/*****************************************************************************
*                                                                            *
*   LOG1P.SPL  Copyright (C) 2016 DSP Development Corporation                *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:    Randy Race                                                    *
*                                                                            *
*   Synopsis:  Computes log(1 + x) for small x                               *
*                                                                            *
*   Revisions: 11 May 2016  RRR  Creation                                    *
*                                                                            *
*****************************************************************************/


#if @HELP_LOG1P

    LOG1P

    Purpose: Computes log(1 + x) for small x.

    Syntax:  LOG1P(x)

              x - A scalar or series.

    Returns: A scalar or series, the value of log(1 + x).

    Example:
             log1p(1e-17)

             returns 1.000000E-17

             log(1 + 1e-17)

             returns 0.0

    Example:
             expm1(log1p(1e-17))

             returns 1.000000E-17

             exp(log(1 + 1e-17)) - 1

             returns 0.0

    Remarks:
             For real values, LOG1P is accurate for small x where
             1 + x == 1 for double precision accuracy.

             For small x, LOG1P(x) is approximately x, but LOG(1 + x)
             is approximately 0 for double precision accuracy.

             See EXPM1 to compute EXP(x) - 1 for small x, the inverse of
             LOG1P.

    See Also:
             Exp
             Expm1
             Log
#endif


/* log(1 + x) for small x */
log1p(x)
{
        local a, b, k, y;

        a = 1 + x;
        b = a - 1;

        /* core */
        y = log(a) - (b - x) / a;

        /* infinities */
        y[find(x == -1)]  = -inf;

        y[find(isinf(x) && x > 0)] = inf;

        k    = find(isinf(x) && x < 0);
        y[k] = inf +  i * imag(y[k]);

        return(y);
}