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


#if @HELP_EXPM1

    EXPM1

    Purpose: Computes exp(x) - 1 for small x.

    Syntax:  EXPM1(x)

              x - A scalar or series.

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

    Example:
             expm1(1e-17)

             returns 1.000000E-17

             exp(1e-17) - 1

             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, EXPM1 is accurate for small x where
             1 + x == 1 for double precision accuracy.

             For small x, EXPM1(x) is approximately x, but EXP(x) - 1
             is approximately 0 for double precision accuracy.

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

    See Also:
             Exp
             Log
             Log1p
#endif


/* compute exp(x) - 1 for small x */
expm1(x)
{
        local y, a, k, range, g_e8;

        y = x;
        a = abs(x);

        range = a >= eps && a <= log(2);

        if (any(range))
        {
                g_e8 = a > 1e-8;
                k    = find(range && g_e8);
                y[k] = exp(x[k]) - 1;

                /* Taylor expansion, more accurate in this range */
                k    = find(range && not(g_e8));
                y[k] = (x[k] / 2 + 1) * x[k];

                /* Newton step for solving log(1 + y) = x for y */
                k     = find(range);
                y[k] -= (1 + y[k]) * (log1p(y[k]) - x[k]);
        }

        /* tiny */
        k    = find(a < eps);
        y[k] = x[k];

        /* negligible cancellation */
        k    = find(a > log(2));
        y[k] = exp(x[k]) - 1;

        return(y);
}