function log1p( x ) {
	var hfsq;
	var hu;
	var y;
	var f;
	var c;
	var s;
	var z;
	var R;
	var u;
	var k;

	if ( x < -1.0 || isnan( x ) ) {
		return NaN;
	}
	if ( x === -1.0 ) {
		return NINF;
	}
	if ( x === PINF ) {
		return x;
	}
	if ( x === 0.0 ) {
		return x; // handle +-0 (IEEE 754-2008 spec)
	}
	// Set y = |x|:
	if ( x < 0.0 ) {
		y = -x;
	} else {
		y = x;
	}
	// Argument reduction...
	k = 1;

	// Check if argument reduction is needed and if we can just return a small value approximation requiring less computation but with equivalent accuracy...
	if ( y < SQRT2M1 ) { // if |x| < sqrt(2)-1 => ~0.41422
		if ( y < SMALL ) { // if |x| < 2**-29
			if ( y < TINY ) { // if |x| < 2**-54
				return x;
			}
			// Use a simple two-term Taylor series...
			return x - ( x*x*0.5 );
		}
		// Check if `f=x` can be represented exactly (no need for correction terms), allowing us to bypass argument reduction...
		if ( x > SQRT2HALFM1 ) { // if x > sqrt(2)/2-1 => ~-0.2929
			// -0.2929 < x < 0.41422
			k = 0;
			f = x; // exact
			hu = 1;
		}
	}
	// Address case where `f` cannot be represented exactly...
	if ( k !== 0 ) {
		if ( y < TWO53 ) {
			u = 1.0 + x;
			hu = getHighWord( u );

			// Bit shift to isolate the exponent and then subtract the bias:
			k = (hu>>20) - FLOAT64_EXPONENT_BIAS;

			// Correction term...
			if ( k > 0 ) { // positive unbiased exponent
				c = 1.0 - (u-x);
			} else { // nonpositive unbiased exponent
				c = x - (u-1.0);
			}
			c /= u;
		} else {
			u = x;
			hu = getHighWord( u );

			// Bit shift to isolate the exponent and then subtract the bias:
			k = (hu>>20) - FLOAT64_EXPONENT_BIAS;

			// Correction term is zero:
			c = 0;
		}
		// Apply a bit mask (0 00000000000 11111111111111111111) to remove the exponent:
		hu &= 0x000fffff; // max value => 1048575

		// Check if u significand is less than sqrt(2) significand => 0x6a09e => 01101010000010011110
		if ( hu < 434334 ) {
			// Normalize u by setting the exponent to 1023 (bias) => 0x3ff00000 => 0 01111111111 00000000000000000000
			u = setHighWord( u, hu|0x3ff00000 );
		} else {
			k += 1;

			// Normalize u/2 by setting the exponent to 1022 (bias-1 => 2**-1 = 1/2) => 0x3fe00000 => 0 01111111110 00000000000000000000
			u = setHighWord( u, hu|0x3fe00000 );

			// Subtract hu significand from next largest hu => 0 00000000001 00000000000000000000 => 0x00100000 => 1048576
			hu = (1048576-hu)>>2;
		}
		f = u - 1.0;
	}
	// Approximation of log1p(f)...
	hfsq = 0.5 * f * f;
	if ( hu === 0 ) { // if |f| < 2**-20
		if ( f === 0.0 ) {
			c += k * LN2_LO;
			return ( k * LN2_HI ) + c;
		}
		R = hfsq * (1.0 - ( TWO_THIRDS*f ) ); // avoid division
		return ( k*LN2_HI ) - ( (R - ( (k*LN2_LO) + c)) - f );
	}
	s = f / (2.0 + f);
	z = s * s;

	R = z * polyval( z );

	if ( k === 0 ) {
		return f - ( hfsq - ( s*(hfsq+R) ) );
	}
	return ( k*LN2_HI ) - ( (hfsq - ( (s*(hfsq+R)) + ((k*LN2_LO) + c))) - f );
}
Evaluates the natural logarithm of \\(1+x\\).

## Method

1.  Argument Reduction: find \\(k\\) and \\(f\\) such that

    ```tex
    1+x = 2^k (1+f)
    ```

    where

    ```tex
    \frac{\sqrt{2}}{2} < 1+f < \sqrt{2}
    ```

    <!-- <note> -->

    If \\(k=0\\), then \\(f=x\\) is exact. However, if \\(k \neq 0\\), then \\(f\\) may not be representable exactly. In that case, a correction term is needed. Let

    ```tex
    u = \operatorname{round}(1+x)
    ```

    and

    ```tex
    c = (1+x) - u
    ```

    then

    ```tex
    \ln (1+x) - \ln u \approx \frac{c}{u}
    ```

    We can thus proceed to compute \\(\ln(u)\\), and add back the correction term \\(c/u\\).

    <!-- </note> -->

    <!-- <note> -->

    When \\(x > 2^{53}\\), one can simply return \\(\ln(x)\\).

    <!-- </note> -->

2.  Approximation of \\(\operatorname{log1p}(f)\\). Let

    ```tex
    s = \frac{f}{2+f}
    ```

    based on

    ```tex
    \begin{align*}
    \ln 1+f &= \ln (1+s) - \ln (1-s) \\
            &= 2s + \frac{2}{3} s^3 + \frac{2}{5} s^5 + ... \\
            &= 2s + sR \\
    \end{align*}
    ```

    We use a special Reme algorithm on \\(\[0,0.1716\]\\) to generate a polynomial of degree \\(14\\) to approximate \\(R\\). The maximum error of this polynomial approximation is bounded by \\(2^{-58.45}\\). In other words,

    ```tex
    R(z) \approx \mathrm{Lp}_1 s^2 + \mathrm{Lp}_2 s^4 + \mathrm{Lp}_3 s^6 + \mathrm{Lp}_4 s^8 + \mathrm{Lp}_5 s^{10} + \mathrm{Lp}_6 s^{12} + \mathrm{Lp}_7 s^{14}
    ```

    and

    ```tex
    | \mathrm{Lp}_1 s^2 + \ldots + \mathrm{Lp}_7 s^14 - R(z) | \leq 2^{-58.45}
    ```

    <!-- <note> -->

    The values of \\(Lp1\\) to \\(Lp7\\) may be found in the source.

    <!-- </note> -->

    Note that

    ```tex
    \begin{align*}
    2s &= f - sf \\
       &= f - \frac{f^2}{2} + s \frac{f^2}{2} \\
    \end{align*}
    ```

    In order to guarantee error in \\(\ln\\) below \\(1\ \mathrm{ulp}\\), we compute the log by

    ```tex
    \operatorname{log1p}(f) = f - \biggl(\frac{f^2}{2} - s\biggl(\frac{f^2}{2}+R\biggr)\biggr)
    ```

3.  Finally,

    ```tex
    \begin{align*}
    \operatorname{log1p}(x) &= k \cdot \mathrm{ln2} + \operatorname{log1p}(f) \\
    &= k \cdot \mathrm{ln2}_{hi}+\biggl(f-\biggl(\frac{f^2}{2}-\biggl(s\biggl(\frac{f^2}{2}+R\biggr)+k \cdot \mathrm{ln2}_{lo}\biggr)\biggr)\biggr) \\
    \end{align*}
    ```

    Here \\(\mathrm{ln2}\\) is split into two floating point numbers:

    ```tex
    \mathrm{ln2}_{hi} + \mathrm{ln2}_{lo}
    ```

    where \\(n \cdot \mathrm{ln2}_{hi}\\) is always exact for \\(|n| < 2000\\).

## Special Cases

-   \\(\operatorname{log1p}(x) = \mathrm{NaN}\\) with signal if \\(x < -1\\) (including \\(-\infty\\))
-   \\(\operatorname{log1p}(+\infty) = +\infty\\)
-   \\(\operatorname{log1p}(-1) = -\infty\\) with signal
-   \\(\operatorname{log1p}(\mathrm{NaN})= \mathrm{NaN}\\) with no signal

## Notes

-   According to an error analysis, the error is always less than \\(1\\) ulp (unit in the last place).

-   The hexadecimal values are the intended ones for the used constants. The decimal values may be used, provided that the compiler will convert from decimal to binary accurately enough to produce the hexadecimal values shown.

-   Assuming \\(\ln(x)\\) is accurate, the following algorithm can be used to evaluate \\(\operatorname{log1p}(x)\\) to within a few ULP:

    ```javascript
    var u = 1.0 + x;
    if ( u === 1.0 ) {
        return x;
    } else {
        return ln(u) * (x/(u-1.0));
    }
    ```

    See HP-15C Advanced Functions Handbook, p.193.

@param {number} x - input value
@returns {number} the natural logarithm of `1+x`

@example
var v = log1p( 4.0 );
// returns ~1.609

@example
var v = log1p( -1.0 );
// returns -Infinity

@example
var v = log1p( 0.0 );
// returns 0.0

@example
var v = log1p( -0.0 );
// returns -0.0

@example
var v = log1p( -2.0 );
// returns NaN

@example
var v = log1p( NaN );
// returns NaN