/*
   DCL (Double precision C Library) - part of Mathlib
*/

#include "mathlib.constants"
#include "mathlib.errors"

/*************************************************************************
 *  Dexp(d) Returns double precision exponential of d
 *
 *  Fortran IV plus user's guide, Digital Equipment Corp. pp B-3
 *
 *  Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *  1968, pp. 96-104.
 *
 *  Inputs greater than ln(MAXPOSDBL) result in overflow.
 *  Inputs less than ln(MINPOSDBL) result in underflow.
 *
 *  The maximum relative error for the approximating polynomial
 *  is 10**(-16.4).  However, this assumes exact arithmetic
 *  in the polynomial evaluation.  Additional rounding and
 *  truncation errors may occur as the argument is reduced
 *  to the range over which the polynomial approximation
 *  is valid, and as the polynomial is evaluated using
 *  finite-precision arithmetic.
 *
 *  Computes exponential from:
 *    exp(X)  =  2**Y  *  2**Z  *  2**W
 *
 *  Where:
 *    Y =  INT ( X * LOG2(e) )
 *    V =  16 * FRAC ( X * LOG2(e))
 *    Z =  (1/16) * INT (V)
 *    W =  (1/16) * FRAC (V)
 *
 *  Note that:
 *    0 =< V < 16
 *    Z = {0, 1/16, 2/16, ...15/16}
 *    0 =< W < 1/16
 *
 *  Then:
 *    2**Z is looked up in a table of 2**0, 2**1/16, ...
 *    2**W is computed from an approximation:
 *    2**W   =  (A + B) / (A - B)
 *
 *    A   =  C + (D * W * W)
 *    B   =  W * (E + (F * W * W))
 *    C   =  20.8137711965230361973
 *    D   =  1.0
 *    E   =  7.2135034108448192083
 *    F   =  0.057761135831801928
 *    Coefficients are from HART, table #1121, pg 206.
 */

static double fpof2tbl [] = {
    1.00000000000000000000,   /* 2 ** 0/16  */
    1.04427378242741384020,   /* 2 ** 1/16  */
    1.09050773266525765930,   /* 2 ** 2/16  */
    1.13878863475669165390,   /* 2 ** 3/16  */
    1.18920711500272106640,   /* 2 ** 4/16  */
    1.24185781207348404890,   /* 2 ** 5/16  */
    1.29683955465100966610,   /* 2 ** 6/16  */
    1.35425554693689272850,   /* 2 ** 7/16  */
    1.41421356237309504880,   /* 2 ** 8/16  */
    1.47682614593949931110,   /* 2 ** 9/16  */
    1.54221082540794082350,   /* 2 ** 10/16 */
    1.61049033194925430820,   /* 2 ** 11/16 */
    1.68179283050742908600,   /* 2 ** 12/16 */
    1.75625216037329948340,   /* 2 ** 13/16 */
    1.83400808640934246360,   /* 2 ** 14/16 */
    1.91520656139714729380    /* 2 ** 15/16 */
};

double Dexp(d)
double d;
{
    int y,index;
    double w,v,a,b,t,wpof2,zpof2;

    double Dint(),Dabs(),Dfrac(),Dpow2();

    if (d > LN_MAXPOSDBL) {
        dcl_error = OVERFLOW;
        return(0.0);
    } else if (d < LN_MINPOSDBL) {
        dcl_error = UNDERFLOW;
        return(0.0);
    } else {
        t = d * LOG2E;
        y = (int) Dint(t);
        v = 16.0 * Dfrac(t);
        index = (int) Dint(Dabs(v));
        if (d < 0.0)
           zpof2 = 1.0 / fpof2tbl[index];
        else
           zpof2 = fpof2tbl[index];
        w = Dfrac(v) / 16.0;
        a = 20.8137711965230361973 + w * w;
        b = w * (7.2135034108448192083 + (0.057761135831801928 * w * w));
        wpof2 = (a + b) / (a - b);
        return ( Dpow2(wpof2*zpof2,y) );
    }
}

/*************************************************************************
 *  Dlog(d) returns the logarithm, base 10, of d
 *
 *  Method : log[base b](d) = log[base e](d) / log[base e](b)
 */
double Dlog(d)
double d;
{
    double Dlogd();

    return( Dlogd(d,10.0) );
}

/*************************************************************************
 *  Dlogd(d1,d2) returns the logarithm, base d2, of the
 *  first double precision argument.
 *
 *  log[d2](d1) = log[e](d1) / log[e](d2)
 */
double Dlogd(d1,d2)
double d1,d2;
{
    double Dln();

    if (d2 <= 0.0) {
        dcl_error = RANGE_ERROR;
        return( 0.0 );
    } else
        return( Dln(d1) / Dln(d2) );
}

/*************************************************************************
 *  Dln(d) Returns double precision natural log of d
 *
 *  Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *  1968, pp. 105-111.
 *
 *  The absolute error in the approximating polynomial
 *  is 10**(-19.38).
 *
 *  This error bound assumes exact arithmetic in the
 *  polynomial evaluation.  Additional rounding and
 *  truncation errors may occur as the argument is reduced
 *  to the range over which the polynomial approximation
 *  is valid, and as the polynomial is evaluated using
 *  finite-precision arithmetic.
 *
 *  Computes ln(x) from:
 *
 *    (1) If argument is zero then flag an error and return minus
 *    infinity (or rather its machine representation).
 *
 *    (2) If argument is negative then flag an error and return minus
 *    infinity.
 *
 *    (3) Given that x = m * 2**k then extract mantissa m and exponent k.
 *
 *    (4) Transform m which is in range [0.5,1.0] to range [1/sqrt(2),
 *    sqrt(2)] by:
 *          s = m * sqrt(2)
 *
 *    (4)  Compute z = (s - 1) / (s + 1)
 *
 *    (5) Now use the approximation from HART page 111 to find ln(s):
 *       ln(s) = z * ( P(z**2) / Q(z**2) )
 *
 *       Where:
 *       P(z**2) =  sum [ Pj * (z**(2*j)) ]
 *       over j = {0,1,2,3}
 *
 *       Q(z**2) =  sum [ Qj * (z**(2*j)) ]
 *       over j = {0,1,2,3}
 *
 *       P0 =  -0.240139179559210509868484e2
 *       P1 =  0.30957292821537650062264e2
 *       P2 =  -0.96376909336868659324e1
 *       P3 =  0.4210873712179797145
 *       Q0 =  -0.120069589779605254717525e2
 *       Q1 =  0.19480966070088973051623e2
 *       Q2 =  -0.89111090279378312337e1
 *       Q3 =  1.0000
 *       (coefficients from HART table #2705 pg 229)
 *
 *  (5) Finally, compute ln(x) from:
 *      ln(x) = (k * ln(2)) - ln(sqrt(2)) + ln(s)
 */

static double dln_pcoeffs[] = {
   -0.24013917955921050986e2,
    0.30957292821537650062e2,
   -0.96376909336868659324e1,
    0.4210873712179797145
};
static double dln_qcoeffs[] = {
   -0.12006958977960525471e2,
    0.19480966070088973051e2,
   -0.89111090279378312337e1,
    1.0000
};

double Dln(d)
double d;
{
    double k, s, z, zt2, pqofz;
    double Dexponent(),Dmantissa(),Dpoly();

    if (d <= 0.0) {
       dcl_error = RANGE_ERROR;
       return( 0.0 );
    } else if (d == 1.0)
       return(0.0);
    else {
       k = Dexponent(d);
       s = SQRT2 * Dmantissa(d);
       z = (s - 1.0) / (s + 1.0);
       zt2 = z * z;
       pqofz = z * ( Dpoly(zt2,3,dln_pcoeffs) / Dpoly(zt2,3,dln_qcoeffs) );
       return ( (k * LN2) - LNSQRT2 + pqofz );
    }
}

/*************************************************************************
 *  Dsqrt(d) return the square root of the d
 *
 *  See Computer Approximations, J.F.Hart et al,
 *  John Wiley & sons, 1968.
 *
 *  The relative error is 10**(-30.1) after three iterations
 *  of Herons iteration for the square root.
 *  This assumes however exact arithmetic in the iterations
 *  and initial approximation. Additional errors may occur due
 *  to truncation, rounding or machine precision limits.
 *
 *  Computes the square root by
 *
 *    (1) Range reduction of the argumentto [0.5,1.0] using the
 *    identity
 *        sqrt(d) = 2**(k/2) * sqrt(x*2**(-k))
 *    where k is the exponent when x is written as a mantissa
 *    times a power of two. It is assumed that the mantissa m is
 *    normalised to the range (0.5 < =m < 1.0)
 *
 *    (2) An approximation to sqrt(m) is obtained from
 *        u = sqrt(m) = (PO + P1*m) / (Q0 + Q1*m)
 *    where PO = 0.594604482
 *          P1 = 2.54164041
 *          Q0 = 2.13725758
 *          Q1 = 1.0
 *    (coefficients from Hart table 360)
 *
 *    (3) A number of iterations of Herons method are performed
 *        y[n+1] = 0.5 * (y[n] + m/y[n])
 *    where y[0] = u = approximate sqrt(m)
 *
 *    (4) If the value of k was odd then y is either multiplied
 *    by the square root of two or divided by the square root of
 *    two for k positive or negative respectively. This rescales
 *    y by multiplying by 2**frac(k/2)
 *
 *    (5) Finally y is rescaled by int(k/2) which is equivalent
 *    to multiplication by 2**int(k/2)
 */

double Dsqrt(d)
double d;
{
    int k, kmod2, rescale, count;
    double m, u, y, ynext;
    double Dint(),Dexponent(),Dmantissa(),Dpow2();

    if (d == 0.0)
        return( 0.0 );
    else if (d < 0.0) {
        dcl_error = SQRT_NEGATIVE;
        return( 0.0 );
    } else {
        k = (int) Dint( Dexponent(d) );
        m = Dmantissa(d);
        u = (0.594604482 + 2.54164041 * m) / (2.13725758 + m);
        for (count = 0, y = u; count < SQRT_ITERATIONS; count++) {
            ynext = 0.5 * (y + (m/y));
            y = ynext;
        }
        rescale = k/2;
        if ((kmod2 = k % 2) < 0)
            y /= SQRT2;
        else if (kmod2 > 0)
            y *= SQRT2;
        return( Dpow2(y,rescale) );
    }
}

/*************************************************************************
 *  Dpow(d1,d2) returns d1 raised to the power of d2.
 *
 *  Computes Dpow(d1,d2) by d1**d2 = exp( d2*ln(d1) )
 */
double Dpow(d1,d2)
double d1,d2;
{
    double Dexp(),Dln();

    if (d2 == 0.0)
        return( 1.0 );
    else if (d1 == 0.0)
        return( 0.0 );
    else if (d1 < 0.0) {
        dcl_error = RANGE_ERROR;
        return( 0.0 );
    }
    d1 = Dln(d1);
    d1 *= d2;
    return( Dexp( d1 ) );
}

/*************************************************************************
 *  Dpow2(d,i) returns d scaled by i**2
 *
 *  Dpow2(d,i) = d * 2**i
 */
double Dpow2(d,i)
double d;
int i;
{
    if (i > 0)
        while (0 < i--)
            d *= 2.0;
    else if (i < 0)
        while (0 > i++)
            d /= 2.0;

    return( d );
}

/*************************************************************************
 *  Duniform() returns a uniformly distributed random number
 *  in the half open interval 0->1
 *
 *  Random number generator - adapted from the FORTRAN version
 *  in "Software Manual for the Elementary Functions"
 *  by W.J. Cody, Jr and William Waite.
 */
double Duniform()
{
    static long int iy = 100001;

    iy *= 125;
    iy %= 2796203;
    return( (double) iy/ 2796203.0 );
}

/*************************************************************************
 *  Dnormal(mean,variance) returns a normally (Gaussian)
 *  distributed random number.
 *  The theory that sample means have a Gaussian probability
 *  distribution is used The numbers are generated by taking
 *  means of sample size 12.
 *
 *  p(x) = [ sum(Ei) - 6 ] * variance + mean
 *  where the Ei's are samples from a uniform distribution
 *  on (0,1).
 */
double Dnormal(mean,variance)
double mean, variance;
{
    int i;
    double sum;
    double Duniform();

    sum = 0.0;
    for (i = 0; i < 12; i++)
        sum += Duniform();
    return( (sum - 6.0) * variance + mean );
}

/*************************************************************************
 *  Dabs(d) returns the absolute value of d
 *  ie : if d < 0 return -d else return d
 */
double Dabs(d)
double d;
{
    return( d < 0 ? -d : d );
}

/*************************************************************************
 *  Dint(d) return the integer part (trunction) of d
 *  eg: Dint(3.4) = 3
 *      Dint(-3.4) = -3
 *
 *  Note : restricted to range of long integers
 */
double Dint(d)
double d;
{
    return( (double) ((long) d) );
}

/*************************************************************************
 *  Dround(d) return the nearest integer to d
 *  eg: Dround(3.4)  = 3
 *      Dround(-3.4) = -3
 *      Dround(3.6)  = 4
 *      Dround(-3.6) = -4
 */
double Dround(d)
double d;
{
    double Dabs(),Dfrac(),Dint(),Dsign();

    return( Dabs( Dfrac(d) ) < 0.50 ? Dint(d) : Dint(d) + Dsign(d) );
}

/*************************************************************************
 *  Dfloor(d) returns the largest integer not greater than d
 *            (round towards negative infinity)
 *  eg: Dfloor(3.4)  = 3
 *      Dfloor(-3.4) = -4
 *      Dfloor(3.0)  = 3.0
 */
double Dfloor(d)
double d;
{
    double Dround();

    return( Dround(d - 0.5) );
}

/*************************************************************************
 *  Dceil(d) returns the smallest integer not less than d
 *           (Round towards positive infinity)
 *  eg: Dceil(3.4)  = 4.0
 *      Dceil(-3.4) = -3.0
 *      Dceil(3.0)  = 3.0
 */
double Dceil(d)
double d;
{
    double Dround();

    return( Dround(d + 0.5) );
}

/*************************************************************************
 *  Dfrac(d) return the fractional part of d
 *  eg : Dfrac(3.4)  = 0.4
 *       Dfrac(-3.4) = - 0.4
 *
 *  d = int(d) + frac(d)
 */
double Dfrac(d)
double d;
{
    double Dint();

    return( d - Dint(d) );
}

/*************************************************************************
 *  Dmod(d1,d2) returns the double precision mudulo of the
 *  two double precision numbers.
 */
double Dmod(d1,d2)
double d1,d2;
{
    double Dint(),Dtrans();

    if (d1 == 0.0)
        return( 0.0 );
    else if (d2 == 0.0) {
        dcl_error = ZERO_DIVIDE;
        return(0.0);
    } else {
        d1 /= d2;
        d1 -= Dtrans(Dint(d1),d1);
        return(d1 *= d2);
    }
}

/*************************************************************************
 *  Dmantissa(d) returns the normalised mantissa of d
 *
 *  ie : d = mantissa * 2 ** exponent
 *  where 1/2 <= mantissa < 1
 *  and the exponent an integer
 */
double Dmantissa(d)
double d;
{
    double mantissa;

    if (d < 0.0)
        mantissa = -d;
    else
        mantissa = d;

    if (mantissa >= 1.0)
        while (mantissa >= 1.0)
            mantissa /= 2.0;
    else
        while (mantissa < 0.5)
            mantissa *= 2.0;

    if (d >= 0)
        return(mantissa);
    else
        return(-mantissa);
}

/*************************************************************************
 *  Dexponent(d) returns the exponent of d
 *
 *  ie : d = mantissa * 2 ** exponent
 *  where 1/2 <= mantissa < 1
 *  and the exponent an integer
 */
double Dexponent(d)
double d;
{
    double exponent = 0.0;

    if (d < 0.0)
        d = -d;

    if (d >= 1.0)
        while (d >= 1.0) {
            d /= 2.0;
            exponent += 1.0;
        }
    else
        while (d < 0.5) {
            d *= 2.0;
            exponent -= 1.0;
        }

    return(exponent);
}

/*************************************************************************
 *  Dmax(d1,d2) returns the maximum of d1 and d2
 */
double Dmax(d1,d2)
double d1,d2;
{
    return( (d1>d2) ? d1 : d2 );
}

/*************************************************************************
 *  Dmin(d1,d2) returns the minimum of d1 and d2
 */
double Dmin(d1,d2)
double d1,d2;
{
    return( (d1<d2) ? d1 : d2 );
}

/*************************************************************************
 *  Dsign(d) returns the sign of d (1 or -1)
 */
double Dsign(d)
double d;
{
    return( (d < 0.0) ? -1.0 : 1.0);
}

/*************************************************************************
 *  Dtrans(d1,d2) returns d1 with the same sign as d2
 */
double Dtrans(d1,d2)
double d1,d2;
{
    if (d1 < 0.0)
        d1 = - d1;
    return( (d2 < 0.0) ? -d1 : d1);
}

/*************************************************************************
 *  Dswap(d1,d2) swaps the two arguments in place
 */
int Dswap(d1,d2)
double *d1,*d2;
{
    double temp;
    temp = *d1;
    *d1 = *d2;
    *d2 = temp;
}

/*************************************************************************
 *    Dpoly(d1,n,d) evaluates a polynomial of degree n at d1
 *    whose coefficients are in the complex array d.
 *
 *    f(d) = a0 + a1 d + a2 d^2 + .... + an d^n
 *    Where the a's are the coefficients in increasing order.
 *    This function evaluates f(d1).
 */
double Dpoly(d1,n,d)
double d1;
int n;
double d[];
{
    double dtmp;
    int i;

    dtmp=d[n];
    for (i = n-1; i >= 0; i--) {
        dtmp *= d1;
        dtmp += d[i];
    }
    return(dtmp);
}

/*************************************************************************
 *  Dsin(d) Returns double precision sine of d
 *
 *  Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *  1968, pp. 112-120.
 *
 *  The dsin and dcos routines are interactive in the sense that
 *  in the process of reducing the argument to the range -PI/4
 *  to PI/4, each may call the other.  Ultimately one or the
 *  other uses a polynomial approximation on the reduced
 *  argument.  The sin approximation has a maximum relative error
 *  of 10**(-17.59) and the cos approximation has a maximum
 *  relative error of 10**(-16.18).
 *
 *  These error bounds assume exact arithmetic
 *  in the polynomial evaluation.  Additional rounding and
 *  truncation errors may occur as the argument is reduced
 *  to the range over which the polynomial approximation
 *  is valid, and as the polynomial is evaluated using
 *  finite-precision arithmetic.
 *
 *  Computes sin(X) from:
 *
 *    (1) Reduce argument X to range -PI to PI.
 *
 *    (2) If X > PI/2 then call sin recursively using relation
 *        sin(X) = -sin(X - PI).
 *
 *    (3) If X < -PI/2 then call sin recursively using relation
 *        sin(X) = -sin(X + PI).
 *
 *    (4) If X > PI/4 then call cos using relation
 *        sin(X) = cos(PI/2 - X).
 *
 *    (5) If X < -PI/4 then call cos using relation
 *        sin(X) = -cos(PI/2 + X).
 *
 *    (6) If X is small enough that polynomial evaluation would cause
 *    underflow then return X, since sin(X) approaches X as X approaches
 *    zero.
 *
 *    (7) By now X has been reduced to range -PI/4 to PI/4 and the
 *    approximation from HART pg. 118 can be used:
 *       sin(X) = Y * ( P(Y) / Q(Y) )
 *
 *       Where:
 *       Y    =  X * (4/PI)
 *
 *       P(Y) =  sum [ Pj * (Y**(2*j)) ]
 *       over j = {0,1,2,3}
 *
 *       Q(Y) =  sum [ Qj * (Y**(2*j)) ]
 *       over j = {0,1,2,3}
 *
 *       P0   =  0.206643433369958582409167054e+7
 *       P1   =  -0.18160398797407332550219213e+6
 *       P2   =  0.359993069496361883172836e+4
 *       P3   =  -0.2010748329458861571949e+2
 *       Q0   =  0.263106591026476989637710307e+7
 *       Q1   =  0.3927024277464900030883986e+5
 *       Q2   =  0.27811919481083844087953e+3
 *       Q3   =  1.0000
 *       (coefficients from HART table #3063 pg 234)
 *
 *  **** NOTE ****   The range reduction relations used in
 *  this routine depend on the final approximation being valid
 *  over the negative argument range in addition to the positive
 *  argument range.  The particular approximation chosen from
 *  HART satisfies this requirement, although not explicitly
 *  stated in the text.  This may not be true of other
 *  approximations given in the reference.
 */

static double dsin_pcoeffs[] = {
    0.20664343336995858240e7,
   -0.18160398797407332550e6,
    0.35999306949636188317e4,
   -0.20107483294588615719e2
};
static double dsin_qcoeffs[] = {
    0.26310659102647698963e7,
    0.39270242774649000308e5,
    0.27811919481083844087e3,
    1.0
};

double Dsin(d)
double d;
{
    double y, yt2;
    double Dmod(),Dcos(),Dpoly();

    if (d < -PI || d > PI) {
       d = Dmod( d , TWOPI );
       if (d > PI)
          d = d - TWOPI;
       else if (d < -PI)
          d = d + TWOPI;
    }
    if (d > PID2)
       return( -(Dsin(d - PI)) );
    else if (d < -PID2)
       return( -(Dsin(d + PI)) );
    else if (d > PID4)
       return( Dcos(PID2 - d) );
    else if (d < -PID4)
       return( -(Dcos(PID2 + d)) );
    else if (d < X6_UNDERFLOWS && d > -X6_UNDERFLOWS)
       return( d );
    else {
       y = d / PID4;
       yt2 = y * y;
       return(y*( Dpoly(yt2,3,dsin_pcoeffs) / Dpoly(yt2,3,dsin_qcoeffs)) );
    }
}

/*************************************************************************
 *  Dcsc(d) returns the inverse of the sine of d
 *
 *  csc(d) = 1 / sin(d)
 */
double Dcsc(d)
double d;
{
    double Dsin();

    d = Dsin(d);
    if (d == 0.0) {
        dcl_error = ZERO_DIVIDE;
        return(0.0);
    } else
        return(1.0 / d);
}

/*************************************************************************
 *  Dcos(d) Returns double precision cosine of d
 *
 *  Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *  1968, pp. 112-120.
 *
 *  The Dsin and Dcos routines are interactive in the sense that
 *  in the process of reducing the argument to the range -PI/4
 *  to PI/4, each may call the other.  Ultimately one or the
 *  other uses a polynomial approximation on the reduced
 *  argument.  The SIN approximation has a maximum relative error
 *  of 10**(-17.59) and the COS approximation has a maximum
 *  relative error of 10**(-16.18).
 *
 *  These error bounds assume exact arithmetic
 *  in the polynomial evaluation.  Additional rounding and
 *  truncation errors may occur as the argument is reduced
 *  to the range over which the polynomial approximation
 *  is valid, and as the polynomial is evaluated using
 *  finite-precision arithmetic.
 *
 *  Computes cos(X) from:
 *
 *  (1) Reduce argument X to range -PI to PI.
 *
 *  (2) If X > PI/2 then call cos recursively using relation
 *      cos(X) = -cos(X - PI).
 *
 *  (3) If X < -PI/2 then call cos recursively using relation
 *      cos(X) = -cos(X + PI).
 *
 *  (4) If X > PI/4 then call sin using relation
 *      cos(X) = sin(PI/2 - X).
 *
 *  (5) If X < -PI/4 then call cos using relation
 *      cos(X) = sin(PI/2 + X).
 *
 *  (6) If X would cause underflow in approx evaluation arithmetic
 *  then return sqrt(1.0 - X**2).
 *
 *  (7) By now X has been reduced to range -PI/4 to PI/4 and the
 *  approximation from HART pg. 119 can be used:
 *   cos(X) = ( P(Y) / Q(Y) )
 *
 *   Where:
 *   Y    = X * (4/PI)
 *
 *   P(Y) =  sum [ Pj * (Y**(2*j)) ] over j = {0,1,2,3}
 *
 *   Q(Y) =  sum [ Qj * (Y**(2*j)) ] over j = {0,1,2,3}
 *
 *   P0   = 0.12905394659037374438571854e+7
 *   P1   = -0.3745670391572320471032359e+6
 *   P2   = 0.134323009865390842853673e+5
 *   P3   = -0.112314508233409330923e+3
 *   Q0   = 0.12905394659037373590295914e+7
 *   Q1   = 0.234677731072458350524124e+5
 *   Q2   = 0.2096951819672630628621e+3
 *   Q3   = 1.0000
 *   (coefficients from HART table #3843 pg 244)
 *
 *  **** NOTE ****    The range reduction relations used in
 *  This routine depend on the final approximation being valid
 *  over the negative argument range in addition to the positive
 *  argument range.  The particular approximation chosen from
 *  HART satisfies this requirement, although not explicitly
 *  stated in the text.  This may not be true of other
 *  approximations given in the reference.
 */

static double dcos_pcoeffs[] = {
    0.12905394659037374438e7,
   -0.37456703915723204710e6,
    0.13432300986539084285e5,
   -0.11231450823340933092e3
};
static double dcos_qcoeffs[] = {
    0.12905394659037373590e7,
    0.23467773107245835052e5,
    0.20969518196726306286e3,
    1.0
};

double Dcos(d)
double d;
{
    double y, yt2;
    double Dmod(),Dsqrt(),Dsin(),Dpoly();

    if (d < -PI || d > PI) {
        d = Dmod( d , TWOPI );
        if (d > PI)
           d = d - TWOPI;
        else if (d < -PI)
           d = d + TWOPI;
    }
    if (d > PID2)
       return( -(Dcos(d - PI)) );
    else if (d < -PID2)
       return( -(Dcos(d + PI)) );
    else if (d > PID4)
       return( Dsin(PID2 - d) );
    else if (d < -PID4)
       return( Dsin(PID2 + d) );
    else if (d < X6_UNDERFLOWS && d > -X6_UNDERFLOWS)
       return( Dsqrt(1.0 - d * d) );
    else {
       y = d / PID4;
       yt2 = y * y;
       return( Dpoly(yt2,3,dcos_pcoeffs) / Dpoly(yt2,3,dcos_qcoeffs));
    }
}

/*************************************************************************
 *  Dsec(d) returns the inverse of the cosine of d
 *
 *  sec(d) = 1 / cos(d)
 */
double Dsec(d)
double d;
{
    double Dcos();

    d = Dcos(d);
    if (d == 0.0) {
        dcl_error = ZERO_DIVIDE;
        return(0.0);
    } else
        return(1.0 / d);
}

/*************************************************************************
 *  Dtan(d) returns the tangent of d
 *
 *  tan(d) = sin(d) / cos(d)
 */
double Dtan(d)
double d;
{
    double sind,cosd;
    double Dsin(),Dcos();

    sind = Dsin(d);
    cosd = Dcos(d);
    if (cosd == 0.0) {
        ccl_error = ZERO_DIVIDE;
        return( 0.0 );
    } else
        return( sind / cosd );
}

/*************************************************************************
 *  Dcot(d) returns the cotangent of d
 */
double Dcot(d)
double d;
{
    double sind,cosd;
    double Dsin(),Dcos();

    sind = Dsin(d);
    cosd = Dcos(d);
    if (sind == 0.0) {
        dcl_error = ZERO_DIVIDE;
        return( 0.0 );
    } else
        return( cosd / sind );
}

/*************************************************************************
 *  Dasin() returns the arc sine of d
 *
 *  Method :
 *    x = abs(d)
 *    asin(d) = atan[ x / sqrt( 1 - x*x ) ] for x > 0.3
 *    asin(d) = atan[ x / sqrt( (1.0 - x) * (1.0 + x) ) ) ] for 0 < x <= 0.3
 */
double Dasin(d)
double d;
{
    double Datan(),Dsqrt();
    double x;

    if (d > 1.0 || d < -1.0) {
        dcl_error = RANGE_ERROR;
        return(0.0);
    } else if (d == 0.0)
        return(0.0);
    else if (d == 1.0)
        return(PID2);
    else if (d == -1.0)
        return(-PID2);
    x = Dabs(d);
    if (x > 0.3)
        return( Datan( x / Dsqrt( (1.0 - x) * (1.0 + x) ) ) );
    else
        return( Datan( x / Dsqrt(1.0 - x * x) ) );
}

/*************************************************************************
 *  Dacsc(d) returns the inverse of the arc sine of d
 *
 *  acsc(d) = asin(1/d)
 */
double Dacsc(d)
double d;
{
    double Dasin();

    if (d == 0.0) {
        dcl_error = ZERO_DIVIDE;
        return(0.0);
    } else
        return( Dasin(1.0 / d) );
}

/*************************************************************************
 *  Dacos(d) returns the arc cosine of the double
 *  precision argument.
 *
 *  acos(d) = atan[ sqrt[ (1 + d) / (1 - d) ] ]
 */
double Dacos(d)
double d;
{
    double y;
    double Datan(),Dsqrt();

    if (d > 1.0 || d < -1.0) {
        dcl_error = RANGE_ERROR;
        return(0.0);
    } else if (d == 0.0)
        return(PID2);
    else if (d == 1.0)
        return(0.0);
    else if (d == -1.0)
        return(PI);
    else
        return( 2.0 * Datan( Dsqrt( (1.0 - d)/(1.0 + d) ) ) );
}

/*************************************************************************
 *  Dasec(d) returns the inverse of the arc cosine of d
 *
 *  asec(d) = acos(1/d)
 */
double Dasec(d)
double d;
{
    double Dacos();

    if (d == 0.0) {
        dcl_error = ZERO_DIVIDE;
        return(0.0);
    } else
        return( Dacos(1.0 / d) );
}

/*************************************************************************
 *  Datan(d) Returns double precision arc tangent of d
 *
 *  Fortran 77 user's guide, Digital Equipment Corp. pp B-3
 *
 *  Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *  1968, pp. 120-130.
 *
 *  The maximum relative error for the approximating polynomial
 *  is 10**(-16.82).  However, this assumes exact arithmetic
 *  in the polynomial evaluation.  Additional rounding and
 *  truncation errors may occur as the argument is reduced to the
 *  range over which the polynomial approximation is valid, and as
 *  the polynomial is evaluated using finite-precision arithmetic.
 *
 *  Computes arctangent(X) from:
 *
 *    (1) If X < 0 then negate X, perform steps 2, 3, and 4,
 *    and negate the returned result. This makes use of the
 *    identity atan(-X) = - atan(X).
 *
 *    (2) If argument X > 1 at this point then test to be sure
 *    that X can be inverted without underflowing.  If not, reduce
 *    X to largest possible number that can be inverted, issue
 *    warning, and continue. Perform steps 3 and 4 with arg = 1/X
 *    and subtract returned result from pi/2. This makes use of the
 *    identity atan(X) = pi/2 - atan(1/X) for X>0.
 *
 *    (3) At this point 0 <= X <= 1.  If X > tan(pi/12) then
 *    perform step 4 with arg = (X*sqrt(3)-1)/(sqrt(3)+X)
 *    and add pi/6 to returned result. This last transformation
 *    maps arguments greater than tan(pi/12) to arguments in
 *    range 0...tan(pi/12).
 *
 *    (4) At this point the argument has been transformed so that
 *    it lies in the range -tan(pi/12)...tan(pi/12). Since very small
 *    arguments may cause underflow in the polynomial evaluation,
 *    a final check is performed.  If the argument is less than the
 *    underflow bound, the function returns x, since atan(X) approaches
 *    asin(X) which approaches X as X goes to zero.
 *
 *    (5) atan(X) is now computed by one of the approximations given
 *    in the cited reference (Hart).  That is:
 *
 *       datan(X) = X * sum [ C[i] * X**(2*i) ]
 *       over i = {0,1,...8}.
 *
 *       Where:
 *       C[0] =   .9999999999999999849899
 *       C[1] =   -.333333333333299308717
 *       C[2] =   .1999999999872944792
 *       C[3] =   -.142857141028255452
 *       C[4] =   .11111097898051048
 *       C[5] =   -.0909037114191074
 *       C[6] =   .0767936869066
 *       C[7] =   -.06483193510303
 *       C[8] =   .0443895157187
 *    (coefficients from HART table #4945 pg 267)
 */

static double datan_coeffs[] = {
    0.9999999999999999849899,
    -0.333333333333299308717,
    0.1999999999872944792,
    -0.142857141028255452,
    0.11111097898051048,
    -0.0909037114191074,
    0.0767936869066,
    -0.06483193510303,
    0.0443895157187
};

double Datan(d)
double d;
{
    double xt2,t1,t2;
    double Datan();Dpoly();

    if (d < 0.0)
       return ( - Datan(-d) );
    else if (d > 1.0)
       return( PID2 - Datan( 1.0 / d ) );
    else if (d > 0.2679491924311227074725) {
       t1 = d * SQRT3 - 1.0;
       t2 = SQRT3 + d;
       return( PID6 + Datan( t1 / t2 ) );
    } else if (d < X16_UNDERFLOW)
       return( d );
    else {
       xt2 = d * d;
       return( d * Dpoly(xt2,8,datan_coeffs) );
    }
}

/*************************************************************************
 *  Dacot(d) returns the inverse of the arc tangent of d
 *
 *  acot(d) = atan(1/d)
 */
double Dacot(d)
double d;
{
    double Dacot();

    if (d == 0.0) {
        dcl_error = ZERO_DIVIDE;
        return(0.0);
    } else
        return( Datan(1.0 / d) );
}

/*************************************************************************
 *  Datan2(d1,d2) returns the double precision arc tangent
 *  of d2 / d1.
 *
 *  atan2 returns the angle in radians between the vector (d1,d2)
 *  and the x axis.
 *
 *  NOTE the order of the arguments.
 */
double Datan2(d1,d2)
double d1, d2;
{
    double Datan();

    if (d2 == 0.0)
        if (d1 < 0.0)
            return( - PI );
        else
            return( 0.0 );
    else if (d1 == 0.0)
        if (d2 < 0.0)
            return( - PID2 );
        else
            return( PID2 );
    else if (d1 > 0.0)
        return( Datan( d2 / d1) );
    else
        if (d2 < 0.0)
            return( Datan( d2 / d1 ) - PI );
        else
            return( Datan( d2 / d1 ) + PI);
}

/*************************************************************************
 *  Dsinh(d) returns the hyperbolic sine of d
 *
 *  sinh(d) = [ exp(d) - exp(-d) ] / 2
 */
double Dsinh(d)
double d;
{
    double Dexp();
    double expd;

    expd = Dexp(d);
    return( 0.5 * ( expd - 1.0 / expd ) );
}

/*************************************************************************
 *  Dcsch(d) returns the inverse of the hyperbolic sine of d
 *
 *  csch(d) = 1 / sinh(d)
 */
double Dcsch(d)
double d;
{
    double Dsinh();

    return(1.0 / Dsinh(d));
}

/*************************************************************************
 *  Dcosh(d) returns the hyperbolic cosine of d
 *
 *  cosh(d) = [ exp(d) + exp(-d) ] / 2
 */
double Dcosh(d)
double d;
{
    double Dexp();
    double expd;

    expd = Dexp(d);
    return( 0.5 * ( expd + 1.0 / expd ) );
}

/*************************************************************************
 *  Dsech(d) returns the inverse of the hyperbolic cosine of d
 *
 *  sech(d) = 1 / cosh(d)
 */
double Dsech(d)
double d;
{
    double Dcosh();

    return(1.0 / Dcosh(d));
}

/*************************************************************************
 *  Dtanh(d) returns the hyperbolic tangent of d
 *
 *    tanh(d) = sinh(d) / cosh(d)
 *            = [ 1 - exp(-2d) ] / [ 1 + exp(-2d) ]
 */
double Dtanh(d)
double d;
{
    double Dexp();
    double exp2d;

    exp2d = Dexp(- 2.0 * d);
    return( (1.0 - exp2d) / (1.0 + exp2d) );
}

/*************************************************************************
 *  Dcoth(d) returns the inverse of the hyperbolic tangent of d
 *
 *  coth(d) = 1 / tanh(d)
 */
double Dcoth(d)
double d;
{
    double Dtanh();

    return(1.0 / Dtanh(d));
}

/*************************************************************************
 *    Dasinh(d) returns the arc hyperbolic sine of d
 *
 *    asinh(d) = ln[ d + sqrt[ d*d + 1 ] ]
 */
double Dasinh(d)
double d;
{
    double Dln(),Dsqrt();

    return( Dln( d + Dsqrt( d * d + 1.0 ) ) );
}

/*************************************************************************
 *  Dacsch(d) returns the arc of the inverse of the
 *  hyperbolic sine of d
 *
 *  acsch(d) = ln[ 1/d + sqrt(1/d^2 + 1) ]
 *  for d != 0
 */
double Dacsch(d)
double d;
{
    double Dln(),Dsqrt();

    if (d == 0.0) {
        dcl_error = ZERO_DIVIDE;
        return(0.0);
    } else
        return( Dln(1.0/d + Dsqrt(1.0 / (d * d) + 1.0)) );
}

/*************************************************************************
 *    Dacosh(d) returns the arc hyperbolic cosine of d
 *
 *    acosh(d) = ln[ d + sqrt[ d*d - 1 ] ]
 *    for d >= 1
 */
double Dacosh(d)
double d;
{
    double Dln(),Dsqrt();

    if (d < 1.0) {
        dcl_error = RANGE_ERROR;
        return(0.0);
    } else
        return( Dln( d + Dsqrt( d * d - 1.0 ) ) );
}

/*************************************************************************
 *  Dasech(d) returns the arc of the inverse of the
 *  hyperbolic cosine of d
 *
 *  asech(d) = ln[ 1/d + sqrt(1/d^2 - 1) ]
 *  for 0 < d <= 1
 */
double Dasech(d)
double d;
{
    double Dln(),Dsqrt();

    if (d <= 0.0 || d > 1.0) {
        dcl_error = RANGE_ERROR;
        return(0.0);
    } else
        return( Dln( 1.0/d + Dsqrt( 1.0 / (d * d) - 1.0 ) ) );
}

/*************************************************************************
 *  Datanh(d) returns the arc hyperbolic tangent of d
 *
 *    atanh(d) = 0.5 * ln[ (1+d)/(1-d) ]
 *    for 0 <= d < 1
 */
double Datanh(d)
double d;
{
    double Dln();

    if (d < 0.0 || d >= 1.0) {
        dcl_error = RANGE_ERROR;
        return(0.0);
    } else
        return( 0.5 * Dln( (1.0 + d) / (1.0 - d) ) );
}

/*************************************************************************
 *  Dacoth(d) returns the arc of the inverse of the
 *  hyperbolic tangent of d
 *
 *  acoth(d) = 0.5 * ln[ (d + 1) / (d - 1) ]
 *  for d^2 > 1
 */
double Dacoth(d)
double d;
{
    double absd;
    double Dln();

    if (d < 0.0)
        absd = -d;
    else
        absd = d;

    if (absd < 1.0) {
        dcl_error = RANGE_ERROR;
        return(0.0);
    } else
        return( 0.5 * Dln( (d + 1.0) / (d - 1.0) ) );
}