unity.c

Go to the documentation of this file.

/*                                                     unity.c
 *
 * Relative error approximations for function arguments near
 * unity.
 *
 *    log1p(x) = log(1+x)
 *    expm1(x) = exp(x) - 1
 *    cosm1(x) = cos(x) - 1
 *    lgam1p(x) = lgam(1+x)
 *
 */
 
/* Scipy changes:
 * - 06-10-2016: added lgam1p
 */
 
#include "mconf.h"
 
extern double MACHEP;
 
 
 
/* log1p(x) = log(1 + x)  */
 
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
 * 1/sqrt(2) <= x < sqrt(2)
 * Theoretical peak relative error = 2.32e-20
 */
static const double LP[] = {
    4.5270000862445199635215E-5,
    4.9854102823193375972212E-1,
    6.5787325942061044846969E0,
    2.9911919328553073277375E1,
    6.0949667980987787057556E1,
    5.7112963590585538103336E1,
    2.0039553499201281259648E1,
};
 
static const double LQ[] = {
    /* 1.0000000000000000000000E0, */
    1.5062909083469192043167E1,
    8.3047565967967209469434E1,
    2.2176239823732856465394E2,
    3.0909872225312059774938E2,
    2.1642788614495947685003E2,
    6.0118660497603843919306E1,
};
 
double log1p(double x)
{
    double z;
 
    z = 1.0 + x;
    if ((z < M_SQRT1_2) || (z > M_SQRT2))
        return (log(z));
    z = x * x;
    z = -0.5 * z + x * (z * polevl(x, LP, 6) / p1evl(x, LQ, 6));
    return (x + z);
}
 
 
/* log(1 + x) - x */
double log1pmx(double x)
{
    if (fabs(x) < 0.5) {
        int n;
        double xfac = x;
        double term;
        double res = 0;
 
        for(n = 2; n < MAXITER; n++) {
            xfac *= -x;
            term = xfac / n;
            res += term;
            if (fabs(term) < MACHEP * fabs(res)) {
                break;
            }
        }
        return res;
    }
    else {
        return log1p(x) - x;
    }
}
 
 
/* expm1(x) = exp(x) - 1  */
 
/*  e^x =  1 + 2x P(x^2)/( Q(x^2) - P(x^2) )
 * -0.5 <= x <= 0.5
 */
 
static double EP[3] = {
    1.2617719307481059087798E-4,
    3.0299440770744196129956E-2,
    9.9999999999999999991025E-1,
};
 
static double EQ[4] = {
    3.0019850513866445504159E-6,
    2.5244834034968410419224E-3,
    2.2726554820815502876593E-1,
    2.0000000000000000000897E0,
};
 
double expm1(double x)
{
    double r, xx;
 
    if (!cephes_isfinite(x)) {
        if (cephes_isnan(x)) {
            return x;
        }
        else if (x > 0) {
            return x;
        }
        else {
            return -1.0;
        }
 
    }
    if ((x < -0.5) || (x > 0.5))
        return (exp(x) - 1.0);
    xx = x * x;
    r = x * polevl(xx, EP, 2);
    r = r / (polevl(xx, EQ, 3) - r);
    return (r + r);
}
 
 
 
/* cosm1(x) = cos(x) - 1  */
 
static double coscof[7] = {
    4.7377507964246204691685E-14,
    -1.1470284843425359765671E-11,
    2.0876754287081521758361E-9,
    -2.7557319214999787979814E-7,
    2.4801587301570552304991E-5,
    -1.3888888888888872993737E-3,
    4.1666666666666666609054E-2,
};
 
double cosm1(double x)
{
    double xx;
 
    if ((x < -M_PI_4) || (x > M_PI_4))
        return (cos(x) - 1.0);
    xx = x * x;
    xx = -0.5 * xx + xx * xx * polevl(xx, coscof, 6);
    return xx;
}
 
 
/* Compute lgam(x + 1) around x = 0 using its Taylor series. */
static double lgam1p_taylor(double x)
{
    int n;
    double xfac, coeff, res;
 
    if (x == 0) {
        return 0;
    }
    res = -SCIPY_EULER * x;
    xfac = -x;
    for (n = 2; n < 42; n++) {
        xfac *= -x;
        coeff = zeta(n, 1) * xfac / n;
        res += coeff;
        if (fabs(coeff) < MACHEP * fabs(res)) {
            break;
        }
    }
    
    return res;
}
 
 
/* Compute lgam(x + 1). */
double lgam1p(double x)
{
    if (fabs(x) <= 0.5) {
        return lgam1p_taylor(x);
    } else if (fabs(x - 1) < 0.5) {
        return log(x) + lgam1p_taylor(x - 1);
    } else {
        return lgam(x + 1);
    }
}