math

/* Copyright 1996 Acorn Computers Ltd
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

/* math.c: ANSI draft (X3J11 May 86) library code, section D.5 */
/* Copyright (C) Codemist Ltd, 1988 */

/* version 0.04b */
/* Nov 87: fix bug in ibm frexp(-ve arg).                                   */

/*
 * This version of the code takes the view that whenever there is an
 * error a NaN should be handed back (as well as errno getting set). The
 * value HUGE_VAL is used, which is not actually a NaN but which will
 * often lead to exponent overflow pretty soon if it is used.  ACN is
 * unclear if this is sensible, and has had a program fall over when
 * atan2(0.0, 0.0) handed back HUGE_VAL rather than some less vicious
 * value (e.g. 0.0).  He can imagine people who expect pow(0.0, 0.0) to
 * be 1.0 (or maybe 0.0, but certainly not HUGE_VAL), and who expect
 * sin(x) to be <= 1.0 in absolute value regardless of anything.  Thus
 * the current state is OK if we are being strict, but mey be unfriendly
 * in some cases?  Thoughts and comments, anybody?
 */

#include "hostsys.h"
#include <limits.h>
#include <errno.h>

/* This file contains code for most of the math routines from <math.h>      */

/* On the ARM some of these routines are implemented as floating point      */
/* opcodes and as such appear in startup.s                                  */

#ifndef NO_FLOATING_POINT

#include <math.h>                          /* for forward references */

#ifdef IBMFLOAT
const double __huge_val = 7.2370055773322621e+75;
#else
/* On the ARM, this has moved into the library's static data area */
/* so that it still works with the Shared C Library module.       */
/* const double __huge_val = 1.79769313486231571e+308; */
#endif

#ifdef IBMFLOAT

double frexp(d, lvn)
double d; int *lvn;
{
    fp_number d1;
    int n;
    if (d==0.0)
/* I worry a little about signed zeros here. I hope that -0.0 == 0.0     */
    {   *lvn = 0;
        return 0.0;
    }
    d1.d = d;
    n = 4*(d1.i.x - 0x40);        /* excess 64 exponent */
    d1.i.x = 0x40;
    d = d1.d;
    /* Note that the following code works for unnormalised numbers, but  */
    /* can then take 55 cycles to converge instead of usual 3 max.       */
    while ((d>=0 ? d : -d) < 0.5) d = d+d, n--;
    /* Now d is most definitely normalised.                              */
    *lvn = n;
    return d;
}


double ldexp(d, n)
double d; int n;
{
    fp_number d1;
    int nx;
    if (d==0.0) return 0.0;         /* special case                      */
    d1.d = d;
    nx = d1.i.x + (n & ~3)/4;
    n &= 3;
#ifndef DO_NOT_SUPPORT_UNNORMALIZED_NUMBERS
    /* The following code gets the msd/expt right for unnormalised nos.  */
    d1.i.x = 0x40;
    d1.d += 0.0;                    /* i.e. normalise                    */
    nx += d1.i.x - 0x40;
#endif
    {   int mhi = d1.i.mhi;
        int nx1 = (mhi & 0x00c00000)==0 ?
                   ((mhi & 0x00200000)==0 ? n - 3 : n - 2) :
                   ((mhi & 0x00800000)==0 ? n - 1 : n);
        if (nx1 > 0) nx++, n -= 4;
/* That just dealt with the fact that in IBM format the exponent is for  */
/* base 16 and so scaling by a power of two can involve a real multiply. */
/* I now know what the true exponent (nx) in the result will be.         */
    }
    if (nx > 0x7f)                  /* Overflow (maybe do a raise() ?)   */
    {   d1.i.x = 0x7f;
        d1.i.mhi = 0xffffff;
        d1.i.mlo = 0xffffffff;
        return d1.d;
    }
    if (nx < 0)                     /* Deal with underflow/unnormalised  */
    {   if (nx <= -14) return 0.0;
        d1.i.x = 0;
        while (nx < 0) d1.d /= 16, nx++;     /* de-normalise gracefully  */
    }
    d1.i.x = nx;
    {   double d2;
        switch (n)
        {
case -3:d2 = 0.125; break;
case -2:d2 = 0.25;  break;
case -1:d2 = 0.5;   break;
default:   /* avoid dataflow whinge */
case 0: d2 = 1.0;   break;
case 1: d2 = 2.0;   break;
case 2: d2 = 4.0;   break;
case 3: d2 = 8.0;   break;
        }
        d1.d *= d2;
    }
    return d1.d;
}


#else       /* Here is the IEEE format stuff */

#ifndef DO_NOT_SUPPORT_UNNORMALIZED_NUMBERS

double frexp(d, lvn)
double d; int *lvn;
{
/* This version works even if d starts off as an unnormalized number in  */
/* the IEEE sense. But in that special case it will be mighty slow!      */
/* By that we mean at most 52 iterations for the smallest number.        */
    fp_number d1;
    int n;
    if (d==0.0)
    {   *lvn = 0;
        return 0.0;
    }
    d1.d = d;
    if ((n = d1.i.x - 0x3fe) == -0x3fe)
    {   int flag;
/* Here d1 has zero in its exponent field - this means that the mantissa */
/* is un-normalized. I have to shift it left (at least one bit) until a  */
/* suitable nonzero bit appears to go in the implicit-bit place in the   */
/* fractional result. For each bit shifted I need to adjust the final    */
/* exponent that will be returned.                                       */
/* I have already tested to see if d was zero so the fllowing loop MUST  */
/* terminate.                                                            */
        do
        {   flag = d1.i.mhi & 0x00080000;
            d1.i.mhi = (d1.i.mhi << 1) | (d1.i.mlo >> 31);
            d1.i.mlo = d1.i.mlo << 1;
            n--;
        } while (flag==0);
    }
    *lvn = n;
    d1.i.x = 0x3fe;
    return d1.d;
}

#else   /* DO_NOT_SUPPORT_UNNORMALIZED_NUMBERS */

double frexp(d, lvn)
double d; int *lvn;
{
    fp_number d1;
    if (d==0.0)
    {   *lvn = 0;
        return 0.0;
    }
    d1.d = d;
    *lvn = d1.i.x - 0x3fe;
    d1.i.x = 0x3fe;
    return d1.d;
}

#endif   /* DO_NOT_SUPPORT_UNNORMALIZED_NUMBERS */

double ldexp(d, n)
double d; int n;
{
    fp_number d1;
    int nx;
    if (d==0.0) return 0.0;        /* special case                       */
    d1.d = d;
    nx = (int) d1.i.x + n;
    if (nx >= 0x7ff)
    { errno = ERANGE;
      return HUGE_VAL;              /* overflow yields 'infinity'        */
    }
/* Maybe I should be prepared to generate un-normalized numbers here, or */
/* even deal with input d un-normalized and n positive yielding a proper */
/* result. All that seems like a lot of work and so I will not even try  */
/* in this version of the code!                                          */
    else if (nx <= 0) return 0.0;  /* deal with underflow                */
    d1.i.x = nx;
    return (d1.d);
}

#endif


#ifdef IBMFLOAT
#define _exp_arg_limit 174.673089501106208
#else
#define _exp_arg_limit 709.78271289338397
#endif

/* machine independent code - but beware the 1e-20's */

#define _pi_       3.14159265358979323846
#define _pi_3      1.04719755119659774615
#define _pi_2      1.57079632679489661923
#define _pi_4      0.78539816339744830962
#define _pi_6      0.52359877559829887038
#define _sqrt_half 0.70710678118654752440

#ifndef HOST_HAS_TRIG
/* The ARM has a load of trig functions etc available as opcodes, so     */
/* has no need for these software versions shown here.                   */

#undef sin
#undef cos
#undef atan

static double _sincos(double x, double y, int sign, int coscase)
{
    int n;
    double xn, f, g, r;
    if (y >= 1.0e9)     /* fail if argument is overlarge                 */
    {   errno = EDOM;
        return -HUGE_VAL;
    }
    n = (int) ((y + _pi_2) / _pi_);
    xn = (double) n;
    if ((n & 1) != 0) sign = -sign;
    if (coscase) xn = xn - 0.5;
/* Note that these days C is REQUIRED to observe the brackets used here. */
/* The compiler is NOT allowed to re-associate the additions.            */
#ifdef IBMFLOAT
    {   double x1 = (double)(int)x;
/* observe that the range check on y assures me that (int)x is OK.       */
        double x2 = x - x1;
        f = ((x1 - xn*3.1416015625) + x2) + xn*8.908910206761537356617e-6;
    }
#else
    f = (x - xn*3.1416015625) + xn*8.908910206761537356617e-6;
#endif
/* I expect that the absolute value of f is less than pi/2 here          */
    if (fabs(f) >= 1.e-10)
#define _sincos_r1  -0.16666666666666665052
#define _sincos_r2   0.83333333333331650315e-2
#define _sincos_r3  -0.19841269841201840457e-3
#define _sincos_r4   0.27557319210152756119e-5
#define _sincos_r5  -0.25052106798274584544e-7
#define _sincos_r6   0.16058936490371589114e-9
#define _sincos_r7  -0.76429178068910467734e-12
#define _sincos_r8   0.27204790957888846175e-14
    {   g = f*f;
        r = ((((((((_sincos_r8) * g + _sincos_r7) * g + _sincos_r6) * g +
                    _sincos_r5) * g + _sincos_r4) * g + _sincos_r3) * g +
                    _sincos_r2) * g + _sincos_r1) * g;
        f += f*r;
    };
    if (sign < 0) return -f;
    else return f;
}

double sin(double x)
{
    if (x < 0.0) return _sincos(-x, -x, -1, 0);
    else return _sincos(x, x, 1, 0);
}

double cos(double x)
{
    if (x < 0.0) return _sincos(-x, _pi_2 - x, 1, 1);
    else return _sincos(x, _pi_2 + x, 1, 1);
}

#ifdef IBMFLOAT
#define _exp_negative_arg_limit -180.218266945585778
#else
/* NB: the following value is a proper limit provided denormalised       */
/* values are not being supported. It would need to be changed if they   */
/* were to start existing.                                               */
#define _exp_negative_arg_limit -708.39641853226408
#endif

double exp(double x)
{
    int n;
    double xn, g, z, gp, q, r;
    if (x > _exp_arg_limit)
    {   errno = ERANGE;
        return HUGE_VAL;
    }
    if (x < _exp_negative_arg_limit) return 0.0;
    if (fabs(x) < 1.e-20) return 1.0;
/* In C the cast (int)x truncates towards zero. Here I want to round.    */
    n = (int)((x >= 0 ? 0.5 : -0.5) + 1.44266950408889634074 * x);
    xn = (double)n;
#ifdef IBMFLOAT
    {   double x1 = (double)(int)x;
        double x2 = x - x1;
        g = ((x1 - xn * 0.693359375) + x2) - xn * (-2.1219444005469058277e-4);
    }
#else
    g = (x - xn * 0.693359375) - xn * (-2.1219444005469058277e-4);
#endif
    z = g * g;
#define  _exp_p0  0.249999999999999993
#define  _exp_p1  0.694360001511792852e-2
#define  _exp_p2  0.165203300268279130e-4
#define  _exp_q0  0.500000000000000000
#define  _exp_q1  0.555538666969001188e-1
#define  _exp_q2  0.495862884905441294e-3
    gp = ((_exp_p2 * z + _exp_p1) * z + _exp_p0) * g;
    q = (_exp_q2 * z + _exp_q1) * z + _exp_q0;
    r = 0.5 + gp / (q - gp);
    return ldexp(r, n + 1);
}

double log(double x)
{
    double f, znum, zden, z, w, r, xn;
    int n;
    if (x <= 0.0)
    {   if (x==0.0)
        {   errno = ERANGE;
            return -HUGE_VAL;
        } else
        { errno = EDOM;
          return -HUGE_VAL;
        }
    }
    f = frexp(x, &n);
    if (f > _sqrt_half)
    {   znum = (f - 0.5) - 0.5;
        zden = f * 0.5 + 0.5;
    }
    else
    {   n -= 1;
        znum = f - 0.5;
        zden = znum * 0.5 + 0.5;
    }
    z = znum / zden;
    w = z * z;
#define _log_a0 -0.64124943423745581147e2
#define _log_a1  0.16383943563021534222e2
#define _log_a2 -0.78956112887491257267e0
#define _log_b0 -0.76949932108494879777e3
#define _log_b1  0.31203222091924532844e3
#define _log_b2 -0.35667977739034646171e2
    r = w * ( ((_log_a2 * w + _log_a1) * w + _log_a0) /
              (((w + _log_b2) * w + _log_b1) * w + _log_b0) );
    r = z + z * r;
    xn = (double)n;
    return (xn*(-2.121944400546905827679e-4) + r) + xn*(355.0/512.0);
}

double log10(double x)
{
    return log(x) * 0.43429448190325182765;  /* log10(e) */
}

double sqrt(double x)
{
    fp_number f;
    double y0;
    int n;
    if (x <= 0.0)
    {   if (x < 0.0) errno = EDOM;
        return -HUGE_VAL;
    }
    f.d = x;
#ifdef IBMFLOAT
    n = f.i.x - 0x40;
    f.i.x = 0x40;
#else
    n = f.i.x - 0x3fe;
    f.i.x = 0x3fe;
#endif
    {   double fd = f.d;
#ifdef IBMFLOAT
        y0 = 0.223607 + 0.894427 * fd;
#else
        y0 = 0.41731 + 0.59016 * fd;
#endif
        y0 = 0.5 * (y0 + fd/y0);
        y0 = 0.5 * (y0 + fd/y0);
        y0 = 0.5 * (y0 + fd/y0);
#ifdef IBMFLOAT
#define __EPS 2.2204460492503131e-16
        y0 = y0 + (0.5 * (fd/y0 - y0) + __EPS/32.0);
#endif
    }
    if (n & 1)
    {
#ifdef IBMFLOAT
        y0 = (y0 + __EPS/8.0) * 0.25;
#else
        y0 *= _sqrt_half;
#endif
        n += 1;
    }
    n /= 2;
    f.d = y0;
    f.i.x += n;
    return f.d;
}

double _tancot(double x, int iflag)
{
    int n;
    double f, g, xnum, xden, y, xn;
    y = fabs(x);
    if (y >= 1.0e9)     /* fail if argument is overlarge                 */
    {   errno = EDOM;
        return -HUGE_VAL;
    }
    n = (int) ((2.0 * y + _pi_2) / _pi_);
    if (x < 0) n = - n;
    xn = (double) n;
#ifdef IBMFLOAT
    {   double x1 = (double)(int)x;
        double x2 = x - x1;
        f = ((x1 - xn*1.57080078125) + x2) + xn*4.454455103380768678308e-6;
    }
#else
    f = (x - xn*1.57080078125) + xn*4.454455103380768678308e-6;
#endif
    if (fabs(f) > 1.e-10)
    {   g = f * f;
#define _tan_p1 -0.13338350006421960681
#define _tan_p2  0.34248878235890589960e-2
#define _tan_p3 -0.17861707342254426711e-4
#define _tan_q0  1.00000000000000000000
#define _tan_q1 -0.46671683339755294240
#define _tan_q2  0.25663832289440112864e-1
#define _tan_q3 -0.31181531907010027307e-3
#define _tan_q4  0.49819433993786512270e-6
        xnum = ((_tan_p3*g + _tan_p2)*g + _tan_p1)*g*f + f;
        xden = (((_tan_q4*g + _tan_q3)*g + _tan_q2)*g + _tan_q1)*g + _tan_q0;
    }
    else
    {   xnum = f;
        xden = 1.0;
    }
/* It is probable that overflow can never occur here, since floating     */
/* point values fall about or more 1.e-16 apart near singularities       */
/* of the tangent function.                                              */
    if (iflag==0)
    {   if ((n & 1) == 0) return xnum / xden;
        else return - xden / xnum;
    }
    else
    {   if ((n & 1) == 0) return xden / xnum;
        else return - xnum / xden;
    }
}

double tan(double x)
{
    return _tancot(x, 0);
}

double _cot(double x)      /* Not specified by ANSI hence the funny name */
{
    if (fabs(x) < 1.0/HUGE_VAL)
    {   errno = ERANGE;
        if (x < 0.0) return -HUGE_VAL;
        else return HUGE_VAL;
    }
    return _tancot(x, 1);
}

double atan(double x)
{
    int n;
    double f;
    const static double a[4] = { 0.0, _pi_6, _pi_2, _pi_3 };
    f = fabs(x);
    if (f > 1.0)
    {   f = 1.0 / f;
        n = 2;
    }
    else n = 0;
#define _two_minus_root_three 0.26794919243112270647
#define _sqrt_three           1.73205080756887729353
#define _sqrt_three_minus_one 0.73205080756887729353
    if (f > _two_minus_root_three)
    {   f = (((_sqrt_three_minus_one*f - 0.5) - 0.5) + f) / (_sqrt_three + f);
        n++;
    }
    if (fabs(f) > 1.e-10)
    {   double g = f * f;
        double r;
#define _atan_p0    -0.13688768894191926929e2
#define _atan_p1    -0.20505855195861651981e2
#define _atan_p2    -0.84946240351320683534e1
#define _atan_p3    -0.83758299368150059274
#define _atan_q0     0.41066306682575781263e2
#define _atan_q1     0.86157349597130242515e2
#define _atan_q2     0.59578436142597344465e2
#define _atan_q3     0.15024001160028576121e2
        r = ((((_atan_p3*g + _atan_p2)*g + _atan_p1)*g + _atan_p0)*g) /
             ((((g + _atan_q3)*g + _atan_q2)*g + _atan_q1)*g + _atan_q0);
        f = f + f * r;
    }
    if (n > 1) f = -f;
    f = f + a[n];
    if (x < 0) return -f;
    else return f;
}

double _asinacos(double x, int flag)
{
    int i;
    double y, g, r;
    const static double a[2] = {0.0, _pi_4 };
    const static double b[2] = {_pi_2, _pi_4 };
    y = fabs(x);
    if (y < 1.e-10) i = flag;
    else
    {   if (y > 0.5)
        {   i = 1 - flag;
            if (y > 1.0)
            {   errno = EDOM;
                return -HUGE_VAL;
            }
            g = ((0.5 - y) + 0.5) * 0.5;
            y = -2.0 * sqrt(g);
        }
        else
        {   i = flag;
            g = y * y;
        }
#define _asin_p1    -0.27368494524164255994e2
#define _asin_p2     0.57208227877891731407e2
#define _asin_p3    -0.39688862997504877339e2
#define _asin_p4     0.10152522233806463645e2
#define _asin_p5    -0.69674573447350646411e0
#define _asin_q0    -0.16421096714498560795e3
#define _asin_q1     0.41714430248260412556e3
#define _asin_q2    -0.38186303361750149284e3
#define _asin_q3     0.15095270841030604719e3
#define _asin_q4    -0.23823859153670238830e2
        r = (((((_asin_p5*g + _asin_p4)*g + _asin_p3)*g +
                              _asin_p2)*g + _asin_p1)*g)      /
             (((((g + _asin_q4)*g + _asin_q3)*g +
                      _asin_q2)*g + _asin_q1)*g + _asin_q0);
        y = y + y*r;
    }
    if (flag==0)
    {   y = (a[i] + y) + a[i];
        if (x<0) y = -y;
    }
    else if (x < 0) y = (b[i] + y) + b[i];
    else y = (a[i] - y) + a[i];
    return y;
}

double asin(double x)
{
    return _asinacos(x, 0);
}

double acos(double x)
{
    return _asinacos(x, 1);
}

double pow(double x, double y)
{
    int sign = 0, m, p, i, mdash, pdash;
    double g, r, z, v, u1, u2;
/* The table a1[] contains properly rounded values for 2**(i/16), and    */
/* a2[] contains differences between the true values of 2**(i/16) and    */
/* the a1[] values for odd i.                                            */
#ifdef IBMFLOAT
    const static double a1[17] = {
/* It is painfully important that the following 17 floating point        */
/* numbers are read in to yield the quantities shown on the right.       */
        1.000000000000000000,    /* 41100000:00000000 */
        0.957603280698573644,    /* 40f5257d:152486cc */
        0.917004043204671229,    /* 40eac0c6:e7dd2439 */
        0.878126080186649740,    /* 40e0ccde:ec2a94e1 */
        0.840896415253714543,    /* 40d744fc:cad69d6b */
        0.805245165974627155,    /* 40ce248c:151f8481 */
        0.771105412703970413,    /* 40c5672a:115506db */
        0.738413072969749659,    /* 40bd08a3:9f580c37 */
        0.707106781186547531,    /* 40b504f3:33f9de65 */
        0.677127773468446367,    /* 40ad583e:ea42a14b */
        0.648419777325504834,    /* 40a5fed6:a9b15139 */
        0.620928906036742028,    /* 409ef532:6091a112 */
        0.594603557501360527,    /* 409837f0:518db8a9 */
        0.569394317378345823,    /* 4091c3d3:73ab11c3 */
        0.545253866332628830,    /* 408b95c1:e3ea8bd7 */
        0.522136891213706919,    /* 4085aac3:67cc487b */
        0.500000000000000000     /* 40800000:00000000 */
        };
    const static double a2[8] = {
        2.4114209503420287e-18,
        9.2291566937243078e-19,
        -1.5241915231122319e-18,
        -3.5421849765286817e-18,
        -3.1286215245415074e-18,
        -4.4654376565694489e-18,
        2.9305146686217562e-18,
        1.1260851040933474e-18
        };
#else /* IEEE format */
    const static double a1[17] = {
/* It is painfully important that the following 17 floating point        */
/* numbers are read in to yield the quantities shown on the right.       */
        1.0,                    /* 3ff00000:00000000 */
        0.9576032806985737,     /* 3feea4af:a2a490da */
        0.91700404320467121,    /* 3fed5818:dcfba487 */
        0.87812608018664973,    /* 3fec199b:dd85529c */
        0.8408964152537145,     /* 3feae89f:995ad3ad */
        0.80524516597462714,    /* 3fe9c491:82a3f090 */
        0.77110541270397037,    /* 3fe8ace5:422aa0db */
        0.73841307296974967,    /* 3fe7a114:73eb0187 */
        0.70710678118654757,    /* 3fe6a09e:667f3bcd */
        0.67712777346844632,    /* 3fe5ab07:dd485429 */
        0.64841977732550482,    /* 3fe4bfda:d5362a27 */
        0.620928906036742,      /* 3fe3dea6:4c123422 */
        0.59460355750136051,    /* 3fe306fe:0a31b715 */
        0.56939431737834578,    /* 3fe2387a:6e756238 */
        0.54525386633262884,    /* 3fe172b8:3c7d517b */
        0.52213689121370688,    /* 3fe0b558:6cf9890f */
        0.5                     /* 3fe00000:00000000 */
        };
    const static double a2[8] = {
        -5.3099730280915798e-17,
        1.4800703477186887e-17,
        1.2353596284702225e-17,
        -1.7419972784343138e-17,
        3.8504741898901863e-17,
        2.3290137959059465e-17,
        4.4563878092065126e-17,
        4.2759448527536718e-17
        };
#endif
    if (y == 1.0) return x;
    if (x <= 0.0)
    {   int ny;
        if (x==0.0)
        {   if (y <= 0.0)
            {   errno = EDOM;
                if (y==0.0) return 1.0;
                return -HUGE_VAL;
            }
            return 0.0;
        }
        if (y < (double)INT_MIN || y > (double)INT_MAX ||
            (double)(ny = (int)y) != y)
        {   errno = EDOM;
            return -HUGE_VAL;
        }
/* Here y is an integer and x is negative.                               */
        x = -x;
        sign = (ny & 1);
    }
    if (y == 2.0 && x < 1.e20) return x*x;  /* special case.             */
    g = frexp(x, &m);
    p = 0;
    if (g <= a1[8]) p = 8;
    if (g <= a1[p+4]) p += 4;
    if (g <= a1[p+2]) p += 2;
    z = ((g - a1[p+1]) - a2[p/2]) / (0.5*g + 0.5*a1[p+1]);
/* Expect abs(z) <= 0.044 here */
    v = z * z;
#define _pow_p1 0.83333333333333211405e-1
#define _pow_p2 0.12500000000503799174e-1
#define _pow_p3 0.22321421285924258967e-2
#define _pow_p4 0.43445775672163119635e-3
    r = (((_pow_p4*v + _pow_p3)*v + _pow_p2)*v + _pow_p1)*v*z;
#define _pow_k 0.44269504088896340736
    r = r + _pow_k * r;
    u2 = (r + z * _pow_k) + z;
#define _reduce(v) ((double)(int)(16.0*(v))*0.0625)
    u1 = (double)(16*m-p-1) * 0.0625;
    {   double y1 = _reduce(y);
        double y2 = y - y1;
        double w = u2*y + u1*y2;
        double w1 = _reduce(w);
        double w2 = w - w1;
        int iw1;
        w = w1 + u1*y1;
        w1 = _reduce(w);
        w2 = w2 + (w - w1);
        w = _reduce(w2);
        iw1 = (int)(16.0*(w1+w));
        w2 = w2 - w;
/* The following values have been determined experimentally, buth their  */
/* values are not very critical.                                         */
#ifdef IBMFLOAT
#  define _negative_pow_limit -4160
#else
#  define _negative_pow_limit -16352
#endif
        if (iw1 < _negative_pow_limit)
        {   errno = ERANGE;         /* Underflow certain                 */
            return 0.0;
        }
        if (w2 > 0.0)
        {   iw1 += 1;
            w2 -= 0.0625;
        }
        if (iw1 < 0) i = 0;
        else i = 1;
        mdash = iw1/16 + i;
        pdash = 16*mdash - iw1;
#define _pow_q1 0.69314718055994529629
#define _pow_q2 0.24022650695909537056
#define _pow_q3 0.55504108664085595326e-1
#define _pow_q4 0.96181290595172416964e-2
#define _pow_q5 0.13333541313585784703e-2
#define _pow_q6 0.15400290440989764601e-3
#define _pow_q7 0.14928852680595608186e-4
        z = ((((((_pow_q7*w2 + _pow_q6)*w2 + _pow_q5)*w2 +
                  _pow_q4)*w2 + _pow_q3)*w2 + _pow_q2)*w2 + _pow_q1)*w2;
        z = a1[pdash] + a1[pdash]*z;
        z = frexp(z, &m);
        mdash += m;
#ifdef IBMFLOAT
        if (mdash > 0x3f*4)
        {   errno = ERANGE;
            if (sign) r = -HUGE_VAL;
            else r = HUGE_VAL;
        }
        else if (mdash <= -0x41*4)
        {   errno = ERANGE;
            r = 0.0;
        }
#else
        if (mdash >= 0x7ff-0x3fe)
        {   errno = ERANGE;
            if (sign) r = -HUGE_VAL;
            else r = HUGE_VAL;
        }
        else if (mdash <= -0x3fe)
        {   errno = ERANGE;
            r = 0.0;
        }
#endif
        else
        {   r = ldexp(z, mdash);
            if (sign) r = -r;
        }
    }
    return r;
}

#endif /* HOST_HAS_TRIG */

double atan2(double y, double x)
{
    if (x==0.0 && y==0.0)
    {   errno = EDOM;
        return -HUGE_VAL;
    }
    if (fabs(x) < fabs(y))
    {   if (fabs(x / y)<1.0e-20)
        {   if (y<0.0) return - _pi_2;
            else return _pi_2;
        }
    }
    y = atan(y / x);
    if (x<0.0)
    {   if (y>0.0) y -= _pi_;
        else y += _pi_;
    }
    return y;
}

#undef fabs

double fabs(double x)
{
    if (x<0.0) return -x;
    else return x;
}

double sinh(double x)
{
    int sign;
    double y;
    if (x<0.0) y = -x, sign = 1; else y = x, sign = 0;
    if (y>1.0)
    {
/* _sinh_lnv is REQUIRED to read in as a number with the lower part of   */
/* its floating point representation zero.                               */
#define    _sinh_lnv     0.69316101074218750000          /* ln(v)        */
#define    _sinh_vm2     0.24999308500451499336          /* 1/v^2        */
#define    _sinh_v2m1    0.13830277879601902638e-4       /* (v/2) - 1    */
        double w = y - _sinh_lnv, z, r;
        if (w>_exp_arg_limit)
        {   errno = ERANGE;
            if (sign) return -HUGE_VAL;
            else return HUGE_VAL;
        }
        z = exp(w);   /* should not overflow!                            */
        if (z < 1.0e10) z = z - _sinh_vm2/z;
        r = z + _sinh_v2m1 * z;
        if (sign) return -r;
        else return r;
    }
    else if (y<=1.0e-10) return x;
    else
    {
#define _sinh_p0    -0.35181283430177117881e6
#define _sinh_p1    -0.11563521196851768270e5
#define _sinh_p2    -0.16375798202630751372e3
#define _sinh_p3    -0.78966127417357099479e0
#define _sinh_q0    -0.21108770058106271242e7
#define _sinh_q1     0.36162723109421836460e5
#define _sinh_q2    -0.27773523119650701667e3
#define _sinh_q3     1.0
        double g = x*x;
        double r;
        /* Use a (minimax) rational approximation. See Cody & Waite.     */
        r = ((((_sinh_p3*g + _sinh_p2)*g + _sinh_p1)*g + _sinh_p0)*g) /
             (((g + _sinh_q2)*g + _sinh_q1)*g + _sinh_q0);
        return x + x*r;
    }
}

double cosh(double x)
{
    if (x<0.0) x = -x;
    if (x>1.0)
    {
        x = x - _sinh_lnv;
        if (x>_exp_arg_limit)
        {   errno = ERANGE;
            return HUGE_VAL;
        }
        x = exp(x);   /* the range check above ensures that this does    */
                      /* not overflow.                                   */
        if (x < 1.0e10) x = x + _sinh_vm2/x;
        /* This very final line might JUST overflow even though the call */
        /* to exp is safe and even though _exp_arg_limit is conservative */
        return x + _sinh_v2m1 * x;
    }
/* This second part is satisfactory, even though it is simple!           */
    x = exp(x);
    return 0.5*(x + 1/x);
}

double tanh(double x)
{
/* The first two exits avoid premature overflow as well as needless use  */
/* of the exp() function.                                                */
    int sign;
    if (x>27.0) return 1.0;         /* here exp(2x) dominates 1.0        */
    else if (x<-27.0) return -1.0;
    if (x<0.0) x = -x, sign = 1;
    else sign = 0;
    if (x>0.549306144334054846)     /* ln(3)/2 is crossover point        */
    {   x = exp(2.0*x);
        x = 1.0 - 2.0/(x + 1.0);
        if (sign) return -x;
        else return x;
    }
#define _tanh_p0    -0.16134119023996228053e4
#define _tanh_p1    -0.99225929672236083313e2
#define _tanh_p2    -0.96437492777225469787e0
#define _tanh_q0     0.48402357071988688686e4
#define _tanh_q1     0.22337720718962312926e4
#define _tanh_q2     0.11274474380534949335e3
#define _tanh_q3     1.0
    if (x>1.0e-10)
    {   double y = x*x;
        /* minimax rational approximation                                */
        y = (((_tanh_p2*y + _tanh_p1)*y + _tanh_p0)*y) /
             (((y + _tanh_q2)*y + _tanh_q1)*y + _tanh_q0);
        x = x + x*y;
    }
    if (sign) return -x;
    else return x;
}

double fmod(double x, double y)
{
/* floating point remainder of (x/y) for integral quotient. Remainder    */
/* has same sign as x.                                                   */
    double q, r;
    if (y==0.0)
    {
      errno = EDOM;
      return -HUGE_VAL;
    }
    if (x==0.0) return 0.0;
    if (y < 0.0) y = -y;
    r = modf(x/y, &q);
    r = x - q * y;
/* The next few lines are an ultra-cautious scheme to ensure that the    */
/* result is less than fabs(y) in value and that it has the sign of x.   */
    if (x > 0.0)
    {   while (r >= y) r -= y;
        while (r < 0.0) r += y;
    }
    else
    {   while (r <= -y) r += y;
        while (r > 0.0) r -= y;
    }
    return r;
}

#ifdef IBMFLOAT

double floor(double d)
{
/* round x down to an integer towards minus infinity.                    */
    fp_number x;
    int exponent, mask, exact;
    if (d == 0.0) return 0.0;
    x.d = d;                            /* pun on union type             */
    if ((exponent = x.i.x - 0x40) < 0)
    {   if (x.i.s) return -1.0;
        else return 0.0;
    }
    else if (exponent >= 56/4) return x.d;
    if (exponent >= 24/4)
    {   mask = ((unsigned) 0xffffffff) >> (4*(exponent - 24));
        exact = x.i.mlo & mask;
        x.i.mlo &= ~mask;
    }
    else
    {   mask = 0xfffff >> (4*exponent);
        exact = (x.i.mhi & mask) | x.i.mlo;
        x.i.mhi &= ~mask;
        x.i.mlo = 0;
    }
    if (exact!=0 && x.i.s) return x.d - 1.0;
    else return x.d;
}

double ceil(double d)
{
/* round x up to an integer towards plus infinity.                       */
    fp_number x;
    int exponent, mask, exact;
    if (d == 0.0) return 0.0;
    x.d = d;                            /* pun on union type             */
    if ((exponent = x.i.x - 0x40) < 0)
    {   if (x.i.s) return 0.0;
        else return 1.0;
    }
    else if (exponent >= 56/4) return x.d;
    if (exponent >= 24/4)
    {   mask = ((unsigned) 0xffffffff) >> (4*(exponent - 24));
        exact = x.i.mlo & mask;
        x.i.mlo &= ~mask;
    }
    else
    {   mask = 0xfffff >> (4*exponent);
        exact = (x.i.mhi & mask) | x.i.mlo;
        x.i.mhi &= ~mask;
        x.i.mlo = 0;
    }
    if (exact!=0 && x.i.s==0) return x.d + 1.0;
    else return x.d;
}

double modf(double value, double *iptr)
{
/* splits value into integral part & fraction (both same sign)           */
    fp_number x;
    int exponent, mask;
    if (value == 0.0)
    {   *iptr = 0.0;
        return 0.0;
    }
    x.d = value;
    if ((exponent = x.i.x - 0x40) < 0)
    {   *iptr = 0.0;
        return value;
    }
    else if (exponent >= 56/4)
    {   *iptr = value;
        return 0.0;
    }
    if (exponent >= 24/4)
    {   mask = ((unsigned) 0xffffffff) >> (4*(exponent - 24));
        x.i.mlo &= ~mask;
    }
    else
    {   mask = 0xffffff >> (4*exponent);
        x.i.mhi &= ~mask;
        x.i.mlo = 0;
    }
    *iptr = x.d;
    return value - x.d;
}

#else /* IBMFLOAT */

double floor(double d)
{
/* round x down to an integer towards minus infinity.                    */
    fp_number x;
    int exponent, mask, exact;
    if (d == 0.0) return 0.0;
    x.d = d;                            /* pun on union type             */
    if ((exponent = x.i.x - 0x3ff) < 0)
    {   if (x.i.s) return -1.0;
        else return 0.0;
    }
    else if (exponent >= 52) return x.d;
    if (exponent >= 20)
    {   mask = ((unsigned) 0xffffffff) >> (exponent - 20);
        exact = x.i.mlo & mask;
        x.i.mlo &= ~mask;
    }
    else
    {   mask = 0xfffff >> exponent;
        exact = (x.i.mhi & mask) | x.i.mlo;
        x.i.mhi &= ~mask;
        x.i.mlo = 0;
    }
    if (exact!=0 && x.i.s) return x.d - 1.0;
    else return x.d;
}

double ceil(double d)
{
/* round x up to an integer towards plus infinity.                       */
    fp_number x;
    int exponent, mask, exact;
    if (d == 0.0) return 0.0;
    x.d = d;                            /* pun on union type             */
    if ((exponent = x.i.x - 0x3ff) < 0)
    {   if (x.i.s) return 0.0;
        else return 1.0;
    }
    else if (exponent >= 52) return x.d;
    if (exponent >= 20)
    {   mask = ((unsigned) 0xffffffff) >> (exponent - 20);
        exact = x.i.mlo & mask;
        x.i.mlo &= ~mask;
    }
    else
    {   mask = 0xfffff >> exponent;
        exact = (x.i.mhi & mask) | x.i.mlo;
        x.i.mhi &= ~mask;
        x.i.mlo = 0;
    }
    if (exact!=0 && x.i.s==0) return x.d + 1.0;
    else return x.d;
}

double modf(double value, double *iptr)
{
/* splits value into integral part & fraction (both same sign)           */
    fp_number x;
    int exponent, mask;
    if (value == 0.0)
    {   *iptr = 0.0;
        return 0.0;
    }
    x.d = value;
    if ((exponent = x.i.x - 0x3ff) < 0)
    {   *iptr = 0.0;
        return value;
    }
    else if (exponent >= 52)
    {   *iptr = value;
        return 0.0;
    }
    if (exponent >= 20)
    {   mask = ((unsigned) 0xffffffff) >> (exponent - 20);
        x.i.mlo &= ~mask;
    }
    else
    {   mask = 0xfffff >> exponent;
        x.i.mhi &= ~mask;
        x.i.mlo = 0;
    }
    *iptr = x.d;
    return value - x.d;
}

#endif /* IBMFLOAT */
#endif /* NO_FLOATING_POINT */

/* end of math.c */