reduce-historical/r37/lisp/csl/cslbase/arith10.c

/*  arith10.c                         Copyright (C) 1990-2002 Codemist Ltd */

/*
 * Arithmetic functions.
 */

/*
 * This code may be used and modified, and redistributed in binary
 * or source form, subject to the "CCL Public License", which should
 * accompany it. This license is a variant on the BSD license, and thus
 * permits use of code derived from this in either open and commercial
 * projects: but it does require that updates to this code be made
 * available back to the originators of the package.
 * Before merging other code in with this or linking this code
 * with other packages or libraries please check that the license terms
 * of the other material are compatible with those of this.
 */


/* Signature: 3f726493 08-Apr-2002 */

#include <stdarg.h>
#include <string.h>
#include <ctype.h>
#include <math.h>

#include "machine.h"
#include "tags.h"
#include "cslerror.h"
#include "externs.h"
#include "arith.h"
#include "entries.h"
#ifdef TIMEOUT
#include "timeout.h"
#endif

/*****************************************************************************/
/***              Transcendental functions etcetera.                       ***/
/*****************************************************************************/


/*
 * Much of the code here is extracted from the portable Fortran library
 * used by Codemist with its Fortran compiler.
 */

#define CSL_log_2 0.6931471805599453094

#ifdef COMMON

static Complex MS_CDECL Cdiv_z(Complex, Complex);

static Complex MS_CDECL cpowi(Complex z, int n)
{
/*
 * Raises a complex number to an integer power by repeated squaring.
 */
    if (n == 0)
    {   Complex one;
        one.real = 1.0;
        one.imag = 0.0;
        return one;
    }
    else if (n < 0)
    {   Complex one;
        one.real = 1.0;
        one.imag = 0.0;
        return Cdiv_z(one, cpowi(z, -n));
    }
    else if (n == 1) return z;
    else
    {   Complex r = cpowi(z, n/2);
        double x1 = r.real, y1 = r.imag;
        double x2, y2;
        if (n & 1)
        {   double x = z.real, y = z.imag;
            x2 = x*x1 - y*y1;
            y2 = x*y1 + y*x1;
        }
        else
        {   x2 = x1;
            y2 = y1;
        }
/*
 * This computes r = z^(n/2) and then evaluates either
 *          r*r
 *  or      (r*z)*r
 * where the the calculation or (r*z) will not overflow unless the
 * final result is going to.  Note that if the sum had been done as
 *          (r*r)*z
 * then in some cases r*r might have overflowed even though r*r*z would
 * not (e.g. r*r giving a just-out-of-range real value, and having
 * magnitude between 1 and sqrt(2), and argument pi/4).
 * Although it is SHODDY, for at least the present I will not bother
 * about overflow in during the process of the final complex
 * multiplication.  After all I suspect ALL complex multiplications
 * generated by the compiler will have the same limitation.
 */
        z.real = x1*x2 - y1*y2;
        z.imag = x1*y2 + y1*x2;
        return z;
    }
}

static Complex MS_CDECL csin(Complex z)
{
    double x = z.real, y = z.imag;
/*
 * sin(x + iy) = sin(x)*cosh(y) + i cos(x)*sinh(y)
 * For smallish y this can be used directly.  For |y| > 50 I will
 * compute sinh and cosh as just +/- exp(|y|)/2
 */
    double s = sin(x), c = cos(x);
    double absy = fabs(y);
    if (absy <= 50.0)
    {   double sh = sinh(y), ch = cosh(y);
        z.real = s * ch;
        z.imag = c * sh;
        return z;
    }
    else
    {
        double w;
        int n = _reduced_exp(absy, &w) - 1;
        z.real = ldexp(s*w, n);
        if (y < 0.0) z.imag = ldexp(-c*w, n);
        else z.imag = ldexp(c*w, n);
        return z;
    }
}


#define CSL_sqrt_starter 0.7071
#define CSL_sqrt_iters   6

static Complex MS_CDECL csqrt(Complex z)
{
    double x = z.real, y = z.imag;
/*
 *
 */
    int j, n;
    double scale;
    double vx, vy;
    if (y == 0.0)
    {   if (x < 0.0)
        {   z.real = 0.0;
            z.imag = sqrt(-x);
        }
        else
        {   z.real = sqrt(x);
            z.imag = 0.0;
        }
        return z;           /* Pure real arguments are easy */
    }
    (void)frexp(y, &n);
/* Exact value returned by frexp not critical */
    if (x != 0.0)
    {   int n1;
        (void)frexp(x, &n1);
        if (n1 > n) n = n1;
    }
    n &= ~1;                /* ensure it is even */
    scale = ldexp(1.0, n);
    x = x / scale;
    y = y / scale; /* Now max(|x|, |y|) is in [0.5, 1) */
/*
 * Generate some kind of starting approximation for a NR iteration.
 * The value selected here here will give a relative error of 1.6e-5
 * after 4 iterations, and so would give about 6.0e-20 after 6, which
 * is more accurate (just) than the machines that I use can cope with.
 * Note that for z near the real axis the starting approximation is
 * either real or (pure) imaginary, thus helping ensure that the near-
 * zero component of the answer comes out to decent accuracy.
 */
    if (y < x)
    {   if (x > -y) vx = CSL_sqrt_starter, vy = 0.0;
        else vx = CSL_sqrt_starter, vy = -CSL_sqrt_starter;
    }
    else if (x > -y) vx = CSL_sqrt_starter, vy = CSL_sqrt_starter;
    else if (y > 0.0) vx = 0.0, vy = CSL_sqrt_starter;
    else vx = 0.0, vy = -CSL_sqrt_starter;
/*
 * For starting values as preconditioned here 6 iterations will converge
 * to about adequate accuracy. For some arguments fewer iterations would
 * suffice, but I am taking the view that the cost of providing a more
 * elaborate end-test might well exceed the cost of the extra iterations
 * that this code performs.  Experiment shows that the worst relative
 * error after 4 iterations is 1.3e-5, so after 6 it will be about
 * 3.0e-20, which is better than machine accuracy.  5 iterations would
 * not be enough.
 */
    for (j=0; j<CSL_sqrt_iters; j++)
    {   double t = vx*vx + vy*vy;
        double qx = (x*vx + y*vy)/t,
               qy = (y*vx - x*vy)/t;
        vx = (vx + qx)*0.5;
        vy = (vy + qy)*0.5;
    }
    n = n/2;
    z.real = ldexp(vx, n);
    z.imag = ldexp(vy, n);
    return z;
}

#undef CSL_sqrt_starter
#undef CSL_sqrt_iters


static Complex MS_CDECL ctan(Complex z)
{
    double x = z.real, y = z.imag;
/*
 * tan(x + iy) = (tan(x) + i tanh(y))/(1 - i tan(x)*tanh(y))
 */
    double t = tan(x), th = tanh(y);
    double t2 = t*t, th2 = th*th;
    /* many risks of premature overflow here */
    double d = 1.0 + t2*th2;
    z.real = t*(1.0-th2)/d;  /* /* if th2 is very near 1.0 this is inaccurate */
    z.imag = th*(1.0+t2)/d;
    return z;
}

static Complex MS_CDECL Cdiv_z(Complex p, Complex q)
{
/*
 * (p/q) as a complex number.  Note abominable issues about scaling so
 * the large values of p and q can be handled properly.
 * The easy scheme for (a + ib)/(c + id) would have been
 *         (ac+bd)/w + i(bc - ad)/w            with w = cc + dd
 */
    double a = p.real, b = p.imag;
    double c = q.real, d = q.imag;
    Complex r;
    if (d == 0.0)      /* dividing by a real number */
    {   r.real = a/c;
        r.imag = b/c;
        return r;
    }
    else if (c == 0.0) /* dividing by a purely imaginary number */
    {   r.real = b/d;
        r.imag = -a/d;
        return r;
    }
    else
    {
        double scalep, scaleq, w;
        int n1, n2, n;
/*
 * I avoid going frexp(0.0, &n1) since the exponent handed back in that
 * case is zero, which is not actually very helpful here, where I would
 * rather it was minus infinity.
 */
        if (a == 0.0)
        {   if (b == 0.0) return p;         /* (0.0, 0.0)/z */
/* Exact results from frexp unimportant */
            else (void)frexp(b, &n1);
        }
        else if (b == 0.0) (void)frexp(a, &n1);
        else
        {   (void)frexp(a, &n1);
            (void)frexp(b, &n2);
            if (n2>n1) n1 = n2;
        }
        n = n1;
        scalep = ldexp(1.0, n1);        /* scale numerator */
        a = a / scalep;
        b = b / scalep;
/* At this stage I know that the denominator has nonzero real & imag parts */
        (void)frexp(c, &n1);
        (void)frexp(d, &n2);
        if (n2>n1) n1 = n2;
        n = n - n1;
        scaleq = ldexp(1.0, n1);        /* scale denominator */
        c = c / scaleq;
        d = d / scaleq;
        w = c*c + d*d;  /* no overflow */
        r.real = ldexp((a*c + b*d)/w, n);   /* rescale final result */
        r.imag = ldexp((b*c - a*d)/w, n);
        return r;
    }
}


static Lisp_Object make_complex_float(Complex v, Lisp_Object a)
/*
 * Here a complex result has been been computed (with double precision values
 * for both real and imag parts).  Squash to the required floating point
 * types and package up as a complex value, where (a) was the original input
 * value and so should defined the type needed.  Note that both
 * components of a will have the same type so only one needs testing.
 * I do the 'onevalue' here.
 */
{
    int32 type;
    Lisp_Object a1, a2, nil;
    a = real_part(a);
    if (is_sfloat(a))
    {   Float_union r, i;
        r.f = (float)v.real;
        i.f = (float)v.imag;
        a1 = make_complex((r.i & ~(int32)0xf) + TAG_SFLOAT,
                          (i.i & ~(int32)0xf) + TAG_SFLOAT);
        errexit();
        return onevalue(a1);
    }
    if (is_bfloat(a)) type = type_of_header(flthdr(a));
    else type = TYPE_SINGLE_FLOAT;
    a1 = make_boxfloat(v.real, type);
    errexit();
    a2 = make_boxfloat(v.imag, type);
    errexit();
    a1 = make_complex(a1, a2);
    errexit();
    return onevalue(a1);
}

#endif

#ifdef __alpha
/*
 * NAG report that they need this for the DEC Alpha... some variant on it
 * may be necessary on other systems too?
 */
#define isnegative(x) ((x) < 0.0 && !isnan(x))
#else
#define isnegative(x) ((x) < 0.0)
#endif

static double MS_CDECL rln(double x)
{
    if (isnegative(x)) x = -x;
    return log(x);
}

static double MS_CDECL iln(double x)
{
    if (isnegative(x)) return _pi;
    else return 0.0;
}

static double MS_CDECL rsqrt(double x)
{
    if (isnegative(x)) return 0.0;
    else return sqrt(x);
}

static double MS_CDECL isqrt(double x)
{
    if (isnegative(x)) return sqrt(-x);
    else return 0.0;
}

static double MS_CDECL rasin(double x)
{
    if (1.0 < x) return _half_pi;
    else if (x <= -1.0) return -_half_pi;
    else return asin(x);
}


static double MS_CDECL iasin(double x)
{
    CSLbool sign;
    if (-1.0 <= x && x <= 1.0) return 0.0;
    else if (x < 0.0) x = -x, sign = YES;
    else sign = NO;
    if (x < 2.0)
    {   x += sqrt(x*x - 1.0);
        x = log(x);             /* /* serious inaccuracy here */
    }
    else if (x < 1.0e9)
    {   x += sqrt(x*x - 1.0);
        x = log(x);
    }
    else x = CSL_log_2 + log(x);
    if (sign) return -x;
    else return x;
}

static double MS_CDECL racos(double x)
{
    if (x <= -1.0) return _pi;
    else if (1.0 <= x) return 0.0;
    else return acos(x);
}

static double MS_CDECL iacos(double x)
{
    CSLbool sign;
    if (x < -1.0) x = -x, sign = YES;
    else if (1.0 < x) sign = NO;
    else return 0.0;
    if (x < 2.0)
    {   x += sqrt(x*x - 1.0);
        x = log(x);             /* /* serious inaccuracy here */
    }
    else if (x < 1.0e9)
    {   x += sqrt(x*x - 1.0);
        x = log(x);
    }
    else x = CSL_log_2 + log(x);
    if (sign) return x;
    else return -x;
}

static double MS_CDECL CSLasinh(double x)
{
    CSLbool sign;
    if (x < 0.0) x = -x, sign = YES;
    else sign = NO;
    if (x < 1.0e-3)
    {   double xx = x*x;
        x = x*(1 - xx*((1.0/6.0) - (3.0/40.0)*xx));
    }
    else if (x < 1.0e9)
    {   x += sqrt(1.0 + x*x);
        x = log(x);
    }
    else x = log(x) + CSL_log_2;
    if (sign) x = -x;
    return x;
}

static double acosh_coeffs[] =
{
 -0.15718655513711019382e-5,          /* x^11 */
 +0.81758779765416234142e-5,          /* x^10 */
 -0.24812280287135584149e-4,          /* x^9 */
 +0.62919005027033514743e-4,          /* x^8 */
 -0.15404104307204835991e-3,          /* x^7 */
 +0.38339903706706128921e-3,          /* x^6 */
 -0.98871347029548821795e-3,          /* x^5 */
 +0.26854094489454297811e-2,          /* x^4 */
 -0.78918167367399344521e-2,          /* x^3 */
 +0.26516504294146930609e-1,          /* x^2 */
 -0.11785113019775570984,             /* x */
 +1.41421356237309504786              /* 1 */

};

static double MS_CDECL racosh(double x)
{
    CSLbool sign;
    if (x < -1.0) x = -x, sign = YES;
    else if (1.0 < x) sign = NO;
    else return 0.0;
    if (x < 1.5)
    {   int i;
        double r = acosh_coeffs[0];
        x = (x - 0.5) - 0.5;
/*
 * This is a minimax approximation to acosh(1+x)/sqrt(x) over the
 * range x=0 to 0.5
 */
        for (i=1; i<=11; i++) r = x*r + acosh_coeffs[i];
        x = sqrt(x)*r;
    }
    else if (x < 1.0e9)
    {   x += sqrt((x - 1.0)*(x + 1.0));
        x = log(x);
    }
    else x = log(x) + CSL_log_2;
    if (sign) return -x;
    else return x;
}

static double MS_CDECL iacosh(double x)
{
    if (1.0 <= x) return 0.0;
    else if (x <= -1.0) return _pi;
    else return acos(x);
}

static double MS_CDECL ratanh(double z)
{
    if (z > -0.01 && z < -0.01)
    {   double zz = z*z;
        return z * (1 + zz*((1.0/3.0) + zz*((1.0/5.0) + zz*(1.0/7.0))));
    }
    z = (1.0 + z) / (1.0 - z);
    if (z < 0.0) z = -z;
    return log(z) / 2.0;
}

static double MS_CDECL iatanh(double x)
{
    if (x < -1.0) return _half_pi;
    else if (1.0 < x) return -_half_pi;
    else return 0.0;
}

/*
 * The functions from here on (for a bit) represent re-work to support the
 * full set of elementary functions that Reduce would (maybe) like. They
 * have not been adjusted to allow use with the software floating point
 * option.
 */

#define n180pi 57.2957795130823208768   /* 180/pi */

#define pi180  0.017453292519943295769  /* pi/180 */

#define sqrthalf 0.7071                 /* sqrt(0.5), low accuracy OK */

static double MS_CDECL racosd(double a)
{
    if (a <= -1.0) return 180.0;
    else if (a < -sqrthalf) return 180.0 - n180pi*acos(-a);
    else if (a < 0.0) return 90.0 + n180pi*asin(-a);
    else if (a < sqrthalf) return 90.0 - n180pi*asin(a);
    else if (a < 1.0) return n180pi*acos(a);
    else return 0.0;
}

static double MS_CDECL iacosd(double a)
/*
 * This version is only good enough for real-mode CSL, not for CCL
 */
{
    if (a >= -1.0 && a <= 1.0) return 0.0;
    else return 1.0;
}

static double MS_CDECL racot(double a)
{
    if (a >= 0.0)
        if (a > 1.0) return atan(1.0/a);
        else return _half_pi - atan(a);
    else if (a < -1.0) return _pi - atan(-1.0/a);
    else return _half_pi + atan(-a);
}

static double MS_CDECL racotd(double a)
{
    if (a >= 0.0)
        if (a > 1.0) return n180pi*atan(1.0/a);
        else return 90.0 - n180pi*atan(a);
    else if (a < -1.0) return 180.0 - n180pi*atan(-1.0/a);
    else return 90.0 + n180pi*atan(-a);
}

static double MS_CDECL racoth(double a)
/*
 * No good in complex case
 */
{
    if (a >= -1.0 && a <= 1.0) return 0.0;
    else return ratanh(1.0/a);
}

static double MS_CDECL iacoth(double a)
{
    if (a >= -1.0 && a <= 1.0) return 1.0;
    else return 0.0;
}

static double MS_CDECL racsc(double a)
{
    if (a > -1.0 && a < 1.0) return 0.0;
    else return asin(1.0/a);
}

static double MS_CDECL iacsc(double a)
{
    if (a > -1.0 && a < 1.0) return 1.0;
    else return 0.0;
}

static double MS_CDECL racscd(double a)
{
    if (a > -1.0 && a < 1.0) return 0.0;
/*
 * I could do better than this, I suspect...
 */
    else return n180pi*asin(1.0/a);
}

static double MS_CDECL iacscd(double a)
{
    if (a > -1.0 && a < 1.0) return 1.0;
    else return 0.0;
}

static double MS_CDECL racsch(double a)
{
    if (a == 0.0) return HUGE_VAL;
    else return CSLasinh(1.0/a);
}

static double MS_CDECL rasec(double a)
{
    if (a > -1.0 && a <= 1.0) return 0.0;
    else return acos(1.0/a);
}

static double MS_CDECL iasec(double a)
{
    if (a > -1.0 && a < 1.0) return 1.0;
    else return 0.0;
}

static double MS_CDECL rasecd(double a)
{
    if (a > -1.0 && a <= 1.0) return 0.0;
/*
 * I could do better than this, I suspect...
 */
    else return n180pi*acos(1.0/a);
}

static double MS_CDECL iasecd(double a)
{
    if (a > -1.0 && a < 1.0) return 1.0;
    else return 0.0;
}

static double MS_CDECL rasech(double a)
{
    if (a <= 0.0 || a >= 1.0) return 0.0;
    else return racosh(1.0/a);
}

static double MS_CDECL iasech(double a)
{
    if (a <= 0.0 || a > 1.0) return 1.0;
    else return 0.0;
}

static double MS_CDECL rasind(double a)
{
    if (a <= -1.0) return -90.0;
    else if (a < -sqrthalf) return -90.0 + n180pi*acos(-a);
    else if (a < sqrthalf) return n180pi*asin(a);
    else if (a < 1.0) return 90.0 - n180pi*acos(a);
    else return 90.0;
}

static double MS_CDECL iasind(double a)
{
    if (a >= -1.0 && a <= 1.0) return 0.0;
    else return 1.0;
}

static double MS_CDECL ratand(double a)
{
    if (a < -1.0) return -90.0 + n180pi*atan(-1.0/a);
    else if (a < 1.0) return n180pi*atan(a);
    else return 90.0 - n180pi*atan(1.0/a);
}

static double MS_CDECL rcbrt(double a)
{
    int xx, x, i, neg = 0;
    double b;
    if (a == 0.0) return 0.0;
    else if (a < 0.0) a = -a, neg = 1;
    b = frexp(a, &xx);  /* end-conditions unimportant */
    x = xx;
/*
 * b is now in the range 0.5 to 1.  The next line produces an
 * approximately minimax linear approximation to the cube root
 % function over the range concerned.
 */
    b = 0.5996 + 0.4081*b;
    while (x % 3 != 0) x--; b *= 1.26;
    b = ldexp(b, x/3);
/*
 * Experiment shows that there are values of the input variable
 * that lead to the last of the following iterations making a
 * (small) change to the result, but almost all of the time I
 * could get away with one less.  Still, I do not consider cbrt
 * important enough to optimise any further (e.g. I could go to
 * use of a minimax rational first approximation...)
 */
    for (i=0; i<6; i++)
        b = (2.0*b + a/(b*b))/3.0;
    if (neg) return -b;
    else return b;
}

static double MS_CDECL rcot(double a)
/*
 * Compare this code with rcotd(). Here I just compute a tangent and
 * form its reciprocal.  What about an arg of pi/2 then, where the
 * tangent calculation might overflow?  Well it should not, since no
 * integer multiple of pi/2 has an exact machine representation, so
 * if you try to compute tan(pi/2) in floating point you should get a
 * large but finite value.  For very large a where the tan() function
 * loses resolution there could still be trouble, which is why I use
 * a slower formula for big a.
 */
{
    if (a > 1000.0 || a < -1000.0) return cos(a)/sin(a);
    else return 1.0/tan(a);
}

static double MS_CDECL arg_reduce_degrees(double a, int *quadrant)
/*
 * Reduce argument to the range -45 to 45, and set quadrant to the
 * relevant quadant.  Returns arg converted to radians.
 */
{
    double w = a / 90.0;
    int32 n = (int)w;
    w = a - 90.0*n;
    while (w < -45.0)
    {   n--;
        w = a - 90.0*n;
    }
    while (w >= 45.0)
    {   n++;
        w = a - 90.0*n;
    }
    *quadrant = (int)(n & 3);
    return pi180*w;
}

static double MS_CDECL rsind(double a)
{
    int quadrant;
    a = arg_reduce_degrees(a, &quadrant);
    switch (quadrant)
    {
default:
case 0: return sin(a);
case 1: return cos(a);
case 2: return sin(-a);
case 3: return -cos(a);
    }
}

static double MS_CDECL rcosd(double a)
{
    int quadrant;
    a = arg_reduce_degrees(a, &quadrant);
    switch (quadrant)
    {
default:
case 0: return cos(a);
case 1: return sin(-a);
case 2: return -cos(a);
case 3: return sin(a);
    }
}

static double MS_CDECL rtand(double a)
{
    int quadrant;
    a = arg_reduce_degrees(a, &quadrant);
    switch (quadrant)
    {
default:
case 0:
case 2: return tan(a);
case 1:
case 3: return 1.0/tan(-a);
    }
}

static double MS_CDECL rcotd(double a)
{
    int quadrant;
    a = arg_reduce_degrees(a, &quadrant);
    switch (quadrant)
    {
default:
case 0:
case 2: return 1.0/tan(a);
case 1:
case 3: return tan(-a);
    }
}


static double MS_CDECL rcoth(double a)
{
    if (a == 0.0) return HUGE_VAL;
    else return 1.0/tanh(a);
}

static double MS_CDECL rcsc(double a)
{
    a = sin(a);
    if (a == 0.0) return HUGE_VAL;
    else return 1.0/a;
}

static double MS_CDECL rcscd(double a)
{
    a = rsind(a);
    if (a == 0.0) return HUGE_VAL;
    else return 1.0/a;
}

static double MS_CDECL rcsch(double a)
{
/*
 * This code is imperfect in that (at least!) exp(-a) can underflow to zero
 * before 2*exp(-a) ought to have.  I will not worry about such refinement
 * here. Much.
 */
    if (a == 0.0) return HUGE_VAL;
    else if (a > 20.0) return 2.0*exp(-a);
    else if (a < -20.0) return -2.0*exp(a);
    else return 1.0/sinh(a);
}

#define CSL_log10 2.302585092994045684

static double MS_CDECL rlog10(double a)
{
    if (a > 0.0) return log(a)/CSL_log10;
    else return 0.0;
}

static double MS_CDECL ilog10(double a)
{
    if (a <= 0) return 1.0;
    else return 0.0;
}

static double MS_CDECL rsec(double a)
{
    a = cos(a);
    if (a == 0.0) return HUGE_VAL;
    else return 1.0/a;
}

static double MS_CDECL rsecd(double a)
{
    a = rcosd(a);
    if (a == 0.0) return HUGE_VAL;
    else return 1.0/a;
}

static double MS_CDECL rsech(double a)
{
/*
 * When |a| is big I ought to return 0.0
 */
    return 1.0/cosh(a);
}

#ifdef COMMON

#define i_times(z) \
    { double temp = z.imag; z.imag = z.real; z.real = -temp; }

#define m_i_times(z) \
    { double temp = z.imag; z.imag = -z.real; z.real = temp; }

#endif

/*
 * The calculations in the next few procedures are numerically
 * crummy, but they should get branch cuts correct.  Re-work later.
 */

#ifdef COMMON

static Complex MS_CDECL casinh(Complex z)
/* log(z + sqrt(1 + z^2)) */
{
    int quadrant = 0;
    Complex w;
    if (z.real < 0.0)
    {   z.real = -z.real;
        z.imag = -z.imag;
        quadrant = 1;
    }
    if (z.imag < 0.0)
    {   z.imag = -z.imag;
        quadrant |= 2;
    }
/* /* The next line can overflow or lose precision */
    w.real = z.real*z.real - z.imag*z.imag + 1.0;
    w.imag = 2*z.real*z.imag;
    w = csqrt(w);
    w.real += z.real;
    w.imag += z.imag;
    w = Cln(w);
    if (quadrant & 2) w.imag = -w.imag;
    if (quadrant & 1) w.real = -w.real, w.imag = -w.imag;
    return w;
}

static Complex MS_CDECL cacosh(Complex z)
/* 2*log(sqrt((z+1)/2) + sqrt((z-1)/2)) */
{
    Complex w1, w2;
    w1.real = (z.real + 1.0)/2.0;
    w2.real = (z.real - 1.0)/2.0;
    w1.imag = w2.imag = z.imag/2.0;
    w1 = csqrt(w1);
    w2 = csqrt(w2);
    w1.real += w2.real;
    w1.imag += w2.imag;
    z = Cln(w1);
    z.real *= 2.0;
    z.imag *= 2.0;
    return z;
}

static Complex MS_CDECL catanh(Complex z)
/* (log(1+z) - log(1-z))/2 */
{
    Complex w1, w2;
    w1.real = 1.0 + z.real;
    w2.real = 1.0 - z.real;
    w1.imag = z.imag;
    w2.imag = -z.imag;
    w1 = Cln(w1);
    w2 = Cln(w2);
    w1.real -= w2.real;
    w1.imag -= w2.imag;
    w1.real /= 2.0;
    w1.imag /= 2.0;
    return w1;
}

static Complex MS_CDECL casin(Complex z)
{
    i_times(z);
    z = casinh(z);
    m_i_times(z);
    return z;
}

static Complex MS_CDECL cacos(Complex z)
{
/*  This is the code I had originally had...
    z = cacosh(z);
    m_i_times(z);
    return z;
*/
/*
 * The following is asserted to behave better. I believe that the
 * calculation (pi/2 - z.real) is guaranteed to introduce severe error
 * when the answer is close to zero, and so this is probably not the ultimate
 * proper formula to use.
 */
    z = casin(z);
    z.real = _half_pi - z.real;
    z.imag = - z.imag;
    return z;
   
}

static Complex MS_CDECL catan(Complex z)
{
    i_times(z);
    z = catanh(z);
    m_i_times(z);
    return z;
}

static Complex MS_CDECL csinh(Complex z)
{
    i_times(z);
    z = csin(z);
    m_i_times(z);
    return z;
}

static Complex MS_CDECL ccosh(Complex z)
{
    i_times(z);
    return Ccos(z);
}

static Complex MS_CDECL ctanh(Complex z)
{
    i_times(z);
    z = ctan(z);
    m_i_times(z);
    return z;
}

/*
 * The next collection of things have not been implemented yet,
 * but might be wanted in a full extended system.... maybe.
 */

static Complex MS_CDECL cacosd(Complex a)
{
    return a;
}


static Complex MS_CDECL cacot(Complex a)
{
    return a;
}


static Complex MS_CDECL cacotd(Complex a)
{
    return a;
}


static Complex MS_CDECL cacoth(Complex a)
{
    return a;
}


static Complex MS_CDECL cacsc(Complex a)
{
    return a;
}


static Complex MS_CDECL cacscd(Complex a)
{
    return a;
}


static Complex MS_CDECL cacsch(Complex a)
{
    return a;
}


static Complex MS_CDECL casec(Complex a)
{
    return a;
}


static Complex MS_CDECL casecd(Complex a)
{
    return a;
}


static Complex MS_CDECL casech(Complex a)
{
    return a;
}


static Complex MS_CDECL casind(Complex a)
{
    return a;
}


static Complex MS_CDECL catand(Complex a)
{
    return a;
}


static Complex MS_CDECL ccbrt(Complex a)
{
    return a;
}


static Complex MS_CDECL ccosd(Complex a)
{
    return a;
}


static Complex MS_CDECL ccot(Complex a)
{
    return a;
}


static Complex MS_CDECL ccotd(Complex a)
{
    return a;
}


static Complex MS_CDECL ccoth(Complex a)
{
    return a;
}


static Complex MS_CDECL ccsc(Complex a)
{
    return a;
}


static Complex MS_CDECL ccscd(Complex a)
{
    return a;
}


static Complex MS_CDECL ccsch(Complex a)
{
    return a;
}


static Complex MS_CDECL clog10(Complex a)
{
    return a;
}


static Complex MS_CDECL csec(Complex a)
{
    return a;
}


static Complex MS_CDECL csecd(Complex a)
{
    return a;
}


static Complex MS_CDECL csech(Complex a)
{
    return a;
}


static Complex MS_CDECL csind(Complex a)
{
    return a;
}


static Complex MS_CDECL ctand(Complex a)
{
    return a;
}

/* end of unimplemented bunch */

#endif
/*
 * Now the Lisp callable entrypoints for the above
 */

typedef double MS_CDECL real_arg_fn(double);
#ifdef COMMON
typedef Complex MS_CDECL complex_arg_fn(Complex);
#endif

typedef struct trigfn
{   double (MS_CDECL *real)(double);
    double (MS_CDECL *imag)(double);
#ifdef COMMON
    Complex (MS_CDECL *complex)(Complex);
#endif
} trigfn_record;

#ifdef NeXT
/*
 * The NeXT header files declare some functions from <math.h> with an
 * extra keyword __const__ in the signature, with a result that placing a
 * pointer to the function in a regularly-typed location leads to a
 * warning. The little wrapper functions here should sort that one out.
 */

static double CSL_atan(double x){ return atan(x); }
#define atan(x) CSL_atan(x) 

static double CSL_cos(double x) { return cos(x);  }
#define cos(x)  CSL_cos(x)  

static double CSL_cosh(double x){ return cosh(x); }
#define cosh(x) CSL_cosh(x) 

static double CSL_exp(double x) { return exp(x);  }
#define exp(x)  CSL_exp(x)  

static double CSL_sin(double x) { return sin(x);  }
#define sin(x)  CSL_sin(x)  

static double CSL_sinh(double x){ return sinh(x); }
#define sinh(x) CSL_sinh(x) 

static double CSL_tan(double x) { return tan(x);  }
#define tan(x)  CSL_tan(x)  

static double CSL_tanh(double x){ return tanh(x); }
#define tanh(x) CSL_tanh(x) 

#endif

#ifdef __LCC__
/*
 * Dec 00: LCC seems to be happier with these...
 */

static double mycos(double a)
{
    return cos(a);
}

static double mysin(double a)
{
    return sin(a);
}

#undef cos
#undef sin
#define cos mycos
#define sin mysin
#endif


#ifdef COMMON
static trigfn_record const trig_functions[] =
{
    {racos,  iacos,  cacos},  /* acos   0  inverse cos, rads, [0, pi) */
    {racosd, iacosd, cacosd}, /* acosd  1  inverse cos, degs, [0, 180) */
    {racosh, iacosh, cacosh}, /* acosh  2  inverse hyperbolic cosine */
    {racot,  0,      cacot},  /* acot   3  inverse cot, rads, (0, pi) */
    {racotd, 0,      cacotd}, /* acotd  4  inverse cot, degs, (0, 180) */
    {racoth, iacoth, cacoth}, /* acoth  5  inverse hyperbolic cotangent */
    {racsc,  iacsc,  cacsc},  /* acsc   6  inverse cosec, [-pi/2, pi/2] */
    {racscd, iacscd, cacscd}, /* acscd  7  inverse cosec, degs, [-90, 90] */
    {racsch, 0,      cacsch}, /* acsch  8  inverse hyperbolic cosecant */
    {rasec,  iasec,  casec},  /* asec   9  inverse sec, rads, [0, pi) */
    {rasecd, iasecd, casecd}, /* asecd  10 inverse sec, degs, [0, 180) */
    {rasech, iasech, casech}, /* asech  11 inverse hyperbolic secant */
    {rasin,  iasin,  casin},  /* asin   12 inverse sin, rads, [-pi/2, pi/2] */
    {rasind, iasind, casind}, /* asind  13 inverse sin, degs, [-90, 90] */
    {CSLasinh, 0,    casinh}, /* asinh  14 inverse hyperbolic sin */
    {atan,   0,      catan},  /* atan   15 1-arg inverse tan, (-pi/2, pi/2) */
    {ratand, 0,      catand}, /* atand  16 inverse tan, degs, (-90, 90) */
    {0,      0,      0},      /* atan2  17 2-arg inverse tan, [0, 2pi) */
    {0,      0,      0},      /* atan2d 18 2-arg inverse tan, degs, [0, 360) */
    {ratanh, iatanh, catanh}, /* atanh  19 inverse hyperbolic tan */
    {rcbrt,  0,      ccbrt},  /* cbrt   20 cube root */
    {cos,    0,      Ccos},   /* cos    21 cosine, rads */
    {rcosd,  0,      ccosd},  /* cosd   22 cosine, degs */
    {cosh,   0,      ccosh},  /* cosh   23 hyperbolic cosine */
    {rcot,   0,      ccot},   /* cot    24 cotangent, rads */
    {rcotd,  0,      ccotd},  /* cotd   25 cotangent, degs */
    {rcoth,  0,      ccoth},  /* coth   26 hyperbolic cotangent */
    {rcsc,   0,      ccsc},   /* csc    27 cosecant, rads */
    {rcscd,  0,      ccscd},  /* cscd   28 cosecant, degs */
    {rcsch,  0,      ccsch},  /* csch   29 hyperbolic cosecant */
    {exp,    0,      Cexp},   /* exp    30 exp(x) = e^z, e approx 2.71828 */
    {0,      0,      0},      /* expt   31 expt(a,b) = a^b */
    {0,      0,      0},      /* hypot  32 hypot(a,b) = sqrt(a^2+b^2) */
    {rln,    iln,    Cln},    /* ln     33 log base e, e approx 2.71828 */
    {0,      0,      0},      /* log    34 2-arg log */
    {rlog10, ilog10, clog10}, /* log10  35 log to base 10 */
    {rsec,   0,      csec},   /* sec    36 secant, rads */
    {rsecd,  0,      csecd},  /* secd   37 secant, degs */
    {rsech,  0,      csech},  /* sech   38 hyperbolic secant */
    {sin,    0,      csin},   /* sin    39 sine, rads */
    {rsind,  0,      csind},  /* sind   40 sine, degs */
    {sinh,   0,      csinh},  /* sinh   41 hyperbolic sine */
    {rsqrt,  isqrt,  csqrt},  /* sqrt   42 square root */
    {tan,    0,      ctan},   /* tan    43 tangent, rads */
    {rtand,  0,      ctand},  /* tand   44 tangent, degs */
    {tanh,   0,      ctanh}   /* tanh   45 hyperbolic tangent */
};

#else
static trigfn_record const trig_functions[] =
{
    {racos,  iacos},   /* acos   0  inverse cosine, rads, [0, pi) */
    {racosd, iacosd},  /* acosd  1  inverse cosine, degs, [0, 180) */
    {racosh, iacosh},  /* acosh  2  inverse hyperbolic cosine */
    {racot,  0},       /* acot   3  inverse cotangent, rads, (0, pi) */
    {racotd, 0},       /* acotd  4  inverse cotangent, degs, (0, 180) */
    {racoth, iacoth},  /* acoth  5  inverse hyperbolic cotangent */
    {racsc,  iacsc},   /* acsc   6  inverse cosecant, rads, [-pi/2, pi/2] */
    {racscd, iacscd},  /* acscd  7  inverse cosecant, degs, [-90, 90] */
    {racsch, 0},       /* acsch  8  inverse hyperbolic cosecant */
    {rasec,  iasec},   /* asec   9  inverse secant, rads, [0, pi) */
    {rasecd, iasecd},  /* asecd  10 inverse secant, degs, [0, 180) */
    {rasech, iasech},  /* asech  11 inverse hyperbolic secant */
    {rasin,  iasin},   /* asin   12 inverse sine, rads, [-pi/2, pi/2] */
    {rasind, iasind},  /* asind  13 inverse sine, degs, [-90, 90] */
    {CSLasinh,  0},    /* asinh  14 inverse hyperbolic sine */
    {atan,   0},       /* atan   15 1-arg inverse tangent, (-pi/2, pi/2) */
    {ratand, 0},       /* atand  16 1-arg inverse tangent, degs, (-90, 90) */
    {0,      0},       /* atan2  17 2-arg inverse tangent, [0, 2pi) */
    {0,      0},       /* atan2d 18 2-arg inverse tangent, degs, [0, 360) */
    {ratanh, iatanh},  /* atanh  19 inverse hyperbolic tangent */
    {rcbrt,  0},       /* cbrt   20 cube root */
    {cos,    0},       /* cos    21 cosine, rads */
    {rcosd,  0},       /* cosd   22 cosine, degs */
    {cosh,   0},       /* cosh   23 hyperbolic cosine */
    {rcot,   0},       /* cot    24 cotangent, rads */
    {rcotd,  0},       /* cotd   25 cotangent, degs */
    {rcoth,  0},       /* coth   26 hyperbolic cotangent */
    {rcsc,   0},       /* csc    27 cosecant, rads */
    {rcscd,  0},       /* cscd   28 cosecant, degs */
    {rcsch,  0},       /* csch   29 hyperbolic cosecant */
    {exp,    0},       /* exp    30 exp(x) = e^z, e approx 2.71828 */
    {0,      0},       /* expt   31 expt(a,b) = a^b */
    {0,      0},       /* hypot  32 hypot(a,b) = sqrt(a^2+b^2) */
    {rln,    iln},     /* ln     33 log base e, e approx 2.71828 */
    {0,      0},       /* log    34 2-arg log. log(a,b) is log of a base b */
    {rlog10, ilog10},  /* log10  35 log to base 10 */
    {rsec,   0},       /* sec    36 secant, rads */
    {rsecd,  0},       /* secd   37 secant, degs */
    {rsech,  0},       /* sech   38 hyperbolic secant */
    {sin,    0},       /* sin    39 sine, rads */
    {rsind,  0},       /* sind   40 sine, degs */
    {sinh,   0},       /* sinh   41 hyperbolic sine */
    {rsqrt,  isqrt},   /* sqrt   42 square root */
    {tan,    0},       /* tan    43 tangent, rads */
    {rtand,  0},       /* tand   44 tangent, degs */
    {tanh,   0}        /* tanh   45 hyperbolic tangent */
};

#endif

static Lisp_Object Ltrigfn(unsigned int which_one, Lisp_Object a)
/*
 * This one piece of code does the type-dispatch for the main collection
 * of elementary functions.
 */
{
    double d;
    Lisp_Object nil = C_nil;
#ifndef COMMON
    int32 restype = TYPE_DOUBLE_FLOAT;
#else
    int32 restype = TYPE_SINGLE_FLOAT;
#endif
    if (which_one > 45) return aerror("trigfn internal error");
    switch ((int)a & TAG_BITS)
    {
case TAG_FIXNUM:
        d = (double)int_of_fixnum(a);
        break;
#ifdef COMMON
case TAG_SFLOAT:
        {   Float_union aa;
            aa.i = a - TAG_SFLOAT;
            d = (double)aa.f;
            restype = 0;
            break;
        }
#endif
case TAG_NUMBERS:
        {   int32 ha = type_of_header(numhdr(a));
            switch (ha)
            {
        case TYPE_BIGNUM:
#ifdef COMMON
        case TYPE_RATNUM:
#endif
                d = float_of_number(a);
                break;
#ifdef COMMON
        case TYPE_COMPLEX_NUM:
                {   Complex c1, c2;
                    c1.real = float_of_number(real_part(a));
                    c1.imag = float_of_number(imag_part(a));
                    c2 = (*trig_functions[which_one].complex)(c1);
/* make_complex_float does the onevalue() for me */
                    return make_complex_float(c2, a);
                }
#endif
        default:
                return aerror1("bad arg for trig function",  a);
            }
            break;
        }
case TAG_BOXFLOAT:
        restype = type_of_header(flthdr(a));
        d = float_of_number(a);
        break;
default:
        return aerror1("bad arg for trig function",  a);
    }
    {   double (MS_CDECL *im)(double) = trig_functions[which_one].imag;
        if (im == 0)
/*
 * I really suspect I should do something writh errno here to
 * keep track of when things go wrong.  Doing so feels fairly
 * messy, but it is necessary if I am to raise exceptions for
 * Lisp if an elementary function leads to overflow.
 */
        {   double (MS_CDECL *rl)(double) = trig_functions[which_one].real;
            if (rl == 0) return aerror("unimplemented trig function");
            d = (*rl)(d);
            a = make_boxfloat(d, restype);
            errexit();
            return onevalue(a);
        }
        else
        {   double c1r, c1i;
            Lisp_Object nil;
#ifdef COMMON
            Lisp_Object rp, ip;
#endif
            double (MS_CDECL *rl)(double) = trig_functions[which_one].real;
            if (rl == 0) return aerror("unimplemented trig function");
            c1r = (*rl)(d);
            c1i = (*im)(d);
/*
 * if the imaginary part of the value is zero then I will return a real
 * answer - this is correct since the original argument was real, but
 * it has to be done by hand here because normally complex values with
 * zero imaginary part remain complex.
 */
            if (c1i == 0.0)
            {
                a = make_boxfloat(c1r, restype);
                errexit();
                return onevalue(a);
            }
#ifdef COMMON
            rp = make_boxfloat(c1r, restype);
            errexit();
            ip = make_boxfloat(c1i, restype);
            errexit();
            a = make_complex(rp, ip);
            errexit();
            return onevalue(a);
#else
            return aerror1("bad arg for trig function",  a);
#endif
        }
    }
}

static Lisp_Object makenum(Lisp_Object a, int32 n)
/*
 * Make the value n, but type-consistent with the object a.  Usually
 * used with n=0 or n=1
 */
{
#ifndef COMMON
    int32 restype = TYPE_DOUBLE_FLOAT;
#else
    int32 restype = TYPE_SINGLE_FLOAT;
#endif
    Lisp_Object nil = C_nil;
    switch ((int)a & TAG_BITS)
    {
case TAG_FIXNUM:
        return fixnum_of_int(n);
#ifdef COMMON
case TAG_SFLOAT:
        {   Float_union aa;
            aa.f = (float)n;
            return (aa.i & ~(int32)0xf) + TAG_SFLOAT;
        }
#endif
case TAG_NUMBERS:
        {   int32 ha = type_of_header(numhdr(a));
            switch (ha)
            {
        case TYPE_BIGNUM:
#ifdef COMMON
        case TYPE_RATNUM:
#endif
                return fixnum_of_int(n);
#ifdef COMMON
        case TYPE_COMPLEX_NUM:
                {   Lisp_Object rr, ii;
                    a = real_part(a);
                    rr = makenum(a, 1);
                    errexit();
                    ii = makenum(a, 0);
                    errexit();
                    a = make_complex(rr, ii);
                    errexit();
                    return onevalue(a);
                }
#endif
        default:
                return aerror1("bad arg for makenumber",  a);
            }
            break;
        }
case TAG_BOXFLOAT:
        restype = type_of_header(flthdr(a));
        a = make_boxfloat((double)n, restype);
        errexit();
        return onevalue(a);
default:
        return aerror1("bad arg for makenumber",  a);
    }
}

static Lisp_Object CSLpowi(Lisp_Object a, unsigned32 n)
/*
 * Raise a to the power n by repeated multiplication. The name is CSLpowi
 * rather than just powi because some miserable C compilers come with an
 * external function called powi in <math.h> and then moan about the
 * name clash.
 */
{
    Lisp_Object nil;
    if (n == 0) return makenum(a, 1); /* value 1 of appropriate type */
    else if (n == 1) return a;
    else if ((n & 1) == 0)
    {   a = CSLpowi(a, n/2);
        errexit();
        return times2(a, a);
    }
    else
    {   Lisp_Object b;
        push(a);
        b = CSLpowi(a, n/2);
        nil = C_nil;
        if (!exception_pending()) b = times2(b, b);
        pop(a);
        errexit();
        return times2(a, b);
    }
}

#ifdef COMMON
static Complex MS_CDECL complex_of_number(Lisp_Object a)
{
    Complex z;
    if (is_numbers(a) && is_complex(a))
    {   z.real = float_of_number(real_part(a));
        z.imag = float_of_number(imag_part(a));
    }
    else
    {   z.real = float_of_number(a);
        z.imag = 0.0;
    }
    return z;
}
#endif

static Lisp_Object Lhypot(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    double u, v, r;
    u = float_of_number(a);
    v = float_of_number(b);
    if (u < 0.0) u = -u;
    if (v < 0.0) v = -v;
    if (u > v) { r = u; u = v; v = r; }
    /* Now 0.0 < u < v */
    if (u + v == v) r = v;  /* u is very small compared to v */
    else
    {   r = u/v;
/*
 * A worry I have is that the multiplication on the following line can
 * overflow, blowing me out of the water.
 */
        r = v * sqrt(1.0 + r*r);
    }
    a = make_boxfloat(r, TYPE_DOUBLE_FLOAT);
    errexit();
    return onevalue(a);
}

Lisp_Object Lexpt(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    double d, e;
    int32 restype, n;
#ifdef COMMON
    Lisp_Object w;
    Complex c1, c2, c3;
#endif
/*
 * I take special action on 1, 0 and -1 raised to a power that is an integer
 * or a bignum. In part this is because raising 1 to a power may be a fairly
 * common case worthy of special care, but the main pressure came because
 * these numbers raised to bignum powers can still have fixnum results, and
 * the value can be computed fast.
 */
    if (a == fixnum_of_int(1) ||
        a == fixnum_of_int(0) ||
        a == fixnum_of_int(-1))
    {   if (is_fixnum(b)) 
        {   n = int_of_fixnum(b);
            switch (int_of_fixnum(a))
            {
        case 1:  return onevalue(a);
        case 0:  if (n < 0) return aerror2("expt", a, b);
                 /* In Common Lisp (expt 0 0) is defined to be 0 */
                 else if (n == 0) return onevalue(fixnum_of_int(1));
                 else return onevalue(a);
        case -1: if (n & 1) return onevalue(a);
                 else return onevalue(fixnum_of_int(1));
            }
        }
        else if (is_numbers(b) && is_bignum(b))
        {   switch (int_of_fixnum(a))
            {
        case 1:  return onevalue(a);
        case 0:  n = bignum_digits(b)[(bignum_length(b)-CELL-4)/4];
                 if (n <= 0) return aerror2("expt", a, b);
                 else return onevalue(a);
        case -1: n = bignum_digits(b)[0];
                 if (n & 1) return onevalue(a);
                 else return onevalue(fixnum_of_int(1));
            }
        }
    }
#ifdef COMMON
/*
 * In a similar vein I will take special action on #C(0 1) and #C(0 -1)
 * raise to integer (including bignum) powers.
 */
    if (is_numbers(a) && is_complex(a) && real_part(a)==fixnum_of_int(0) &&
        (imag_part(a) == fixnum_of_int(1) ||
         imag_part(a) == fixnum_of_int(-1)))
    {   n = -1;
        if (is_fixnum(b)) n = int_of_fixnum(b) & 3;
        else if (is_numbers(b) && is_bignum(b))
            n = bignum_digits(b)[0] & 3;
        switch (n)
        {
    case 0:   return onevalue(fixnum_of_int(1));
    case 2:   return onevalue(fixnum_of_int(-1));
    case 1:
    case 3:   if (int_of_fixnum(imag_part(a)) == 1) n ^= 2;
              a = make_complex(fixnum_of_int(0),
                               fixnum_of_int((n & 2) ? 1 : -1));
              errexit();
              return onevalue(a);
    default:  break;
        }
    }
#endif
    if (is_fixnum(b))   /* bignum exponents would yield silly values! */
    {   n = int_of_fixnum(b);
        if (n < 0)
        {   a = CSLpowi(a, (unsigned32)(-n));
            nil = C_nil;
#ifdef COMMON
            if (!exception_pending()) a = CLquot2(fixnum_of_int(1), a);
#else
            if (!exception_pending()) a = quot2(fixnum_of_int(1), a);
#endif
        }
        else a = CSLpowi(a, (unsigned32)n);
        return onevalue(a);
    }
#ifdef COMMON
    if (is_numbers(a) && is_complex(a)) w = real_part(a);
    else w = a;
    if (is_sfloat(w)) restype = 0;
    else
    if (is_bfloat(w)) restype = type_of_header(flthdr(w));
    else restype = TYPE_SINGLE_FLOAT;
    if (is_numbers(b) && is_complex(b)) w = real_part(b);
    else w = b;
    if (is_bfloat(w))
    {   n = type_of_header(flthdr(w));
        if (restype == 0 || restype == TYPE_SINGLE_FLOAT ||
            (restype == TYPE_DOUBLE_FLOAT && n == TYPE_LONG_FLOAT))
            restype = n;
    }
    else if (restype == 0) restype = TYPE_SINGLE_FLOAT;
#else
    restype = TYPE_DOUBLE_FLOAT;
#endif
#ifdef COMMON
    if ((is_numbers(a) && is_complex(a)) ||
        (is_numbers(b) && is_complex(b)))
    {   c1 = complex_of_number(a);
        c2 = complex_of_number(b);
        c3 = Cpow(c1, c2);
        a = make_boxfloat(c3.real, restype);
        errexit();
        b = make_boxfloat(c3.imag, restype);
        errexit();
        a = make_complex(a, b);
        errexit();
        return onevalue(a);
    }
#endif
    d = float_of_number(a);
    e = float_of_number(b);
    if (d < 0.0)
#ifdef COMMON
    {   c1.real = d; c1.imag = 0.0;
        c2.real = e; c2.imag = 0.0;
        c3 = Cpow(c1, c2);
        a = make_boxfloat(c3.real, restype);
        errexit();
        b = make_boxfloat(c3.imag, restype);
        errexit();
        a = make_complex(a, b);
        errexit();
        return onevalue(a);
    }
#else
        return aerror1("bad arg for expt",  b);
#endif
    d = pow(d, e);
    a = make_boxfloat(d, restype);
    errexit();
    return onevalue(a);
}

Lisp_Object Llog_2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
/*
 * Log with specified base.
 */
{
    push(b);
    a = Ltrigfn(33, a);
    pop(b);
    errexit();
    push(a);
    b = Ltrigfn(33, b);
    pop(a);
    errexit();
    return quot2(a, b);
}

#ifdef COMMON

static Lisp_Object Lisqrt(Lisp_Object nil, Lisp_Object a)
{
    double d;
    CSL_IGNORE(nil);
    switch ((int)a & TAG_BITS)
    {
case TAG_FIXNUM:
        d = (double)int_of_fixnum(a);
        break;
case TAG_NUMBERS:
        {   int32 ha = type_of_header(numhdr(a));
            switch (ha)
            {
        case TYPE_BIGNUM:
                d = float_of_number(a);
                break;
        default:
                return aerror1("bad arg for isqrt",  a);
            }
            break;
        }
default:
        return aerror1("bad arg for isqrt",  a);
    }
    d = sqrt(d);
/* /* This is not anything like good enough yet */
    return onevalue(fixnum_of_int((int32)d));
}

#endif

Lisp_Object Labsval(Lisp_Object nil, Lisp_Object a)
/*
 * I call this Labsval not Labs because a non-case-sensitive linker
 * would confuse Labs with labs, and labs is defined in the C libraries...
 * Of course I do not think that case-insensitive linkers should be allowed
 * to remain in service....
 */
{
    switch ((int)a & TAG_BITS)
    {
case TAG_FIXNUM:
#ifdef COMMON
case TAG_SFLOAT:
#endif
        break;
case TAG_NUMBERS:
        {   int32 ha = type_of_header(numhdr(a));
            switch (ha)
            {
        case TYPE_BIGNUM:
#ifdef COMMON
        case TYPE_RATNUM:
#endif
                break;
#ifdef COMMON
        case TYPE_COMPLEX_NUM:
                {   Complex c1;
                    double d;
                    c1.real = float_of_number(real_part(a));
                    c1.imag = float_of_number(imag_part(a));
                    d = Cabs(c1);
/* /* I wonder if I am allowed to promote short or single values to
      double precision here? */
                    a = make_boxfloat(d, TYPE_DOUBLE_FLOAT);
                    errexit();
                    return onevalue(a);
                }
#endif
        default:
                return aerror1("bad arg for abs",  a);
            }
            break;
        }
case TAG_BOXFLOAT:
        break;
default:
        return aerror1("bad arg for abs",  a);
    }
    if (minusp(a))
    {   nil = C_nil;
        if (!exception_pending()) a = negate(a);
    }
    errexit();
    return onevalue(a);
}

#ifdef COMMON

static Lisp_Object Lphase(Lisp_Object nil, Lisp_Object a)
{
    CSLbool s;
    double d;
    if (is_numbers(a) && is_complex(a))
        return Latan2(nil, imag_part(a), real_part(a));
    s = minusp(a);
    errexit();
    if (s) d = -_pi;
    else d = _pi;
    a = make_boxfloat(d, TYPE_DOUBLE_FLOAT);
    errexit();
    return onevalue(a);
/* /* Wrong precision, I guess */
}

static Lisp_Object Lsignum(Lisp_Object nil, Lisp_Object a)
{
/*
 * This seems an expensive way of doing things - huh? Maybe complex values? */
    CSLbool z;
    Lisp_Object w;
    z = zerop(a);
    nil = C_nil;
    if (z || exception_pending()) return onevalue(a);
    push(a);
    w = Labsval(nil, a);
    pop(a);
    errexit();
    a = quot2(a, w);
    errexit();
    return onevalue(a);
}

static Lisp_Object Lcis(Lisp_Object env, Lisp_Object a)
/*
 * Implement as exp(i*a) - this permits complex args which goes
 * beyond the specification of Common Lisp.
 */
{
    Lisp_Object ii, nil;
    CSL_IGNORE(env);
    push(a);
    ii = make_complex(fixnum_of_int(0), fixnum_of_int(1));
    pop(a);
    errexit();
/*
 * it seems a bit gross to multiply by i by calling times2(), but
 * doing so avoids loads of messy type dispatch code here and
 * I am not over-worried about performance at this level (yet).
 */
    a = times2(a, ii);
    errexit();
    return Ltrigfn(30, a);     /* exp() */
}
#endif

#ifdef COMMON
Lisp_Object Latan(Lisp_Object env, Lisp_Object a)
{
    CSL_IGNORE(env);
    return Ltrigfn(15, a);     /* atan() */
}

Lisp_Object Latan_2(Lisp_Object env, Lisp_Object a, Lisp_Object b)
{
    CSL_IGNORE(env);
    return Latan2(env, a, b);
}
#else
Lisp_Object Latan(Lisp_Object env, Lisp_Object a)
{
    CSL_IGNORE(env);
    return Ltrigfn(15, a);     /* atan() */
}
#endif

Lisp_Object Latan2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    double u, v, r;
    int q = 0;
    v = float_of_number(a);
    u = float_of_number(b);
    if (u < 0.0) u = -u, q = 1;
    if (v < 0.0) v = -v, q |= 2;
    if (v > u) { r = u; u = v; v = r; q |= 4; }
    if (u == 0.0 && v == 0.0) r = 0.0;
    else
    {   r = atan(v/u);
        switch (q)
        {
    default:
    case 0: break;
    case 1: r = _pi - r;
            break;
    case 2: r = -r;
            break;
    case 3: r = -_pi + r;
            break;
    case 4: r = _half_pi - r;
            break;
    case 5: r = _half_pi + r;
            break;
    case 6: r = -_half_pi + r;
            break;
    case 7: r = -_half_pi - r;
            break;
        }
    }
    a = make_boxfloat(r, TYPE_DOUBLE_FLOAT);
    errexit();
    return onevalue(a);
}

Lisp_Object Latan2d(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
{
    double u, v, r;
    int q = 0;
    v = float_of_number(a);
    u = float_of_number(b);
    if (u < 0.0) u = -u, q = 1;
    if (v < 0.0) v = -v, q |= 2;
    if (v > u) { r = u; u = v; v = r; q |= 4; }
    if (u == 0.0 && v == 0.0) r = 0.0;
    else
    {   r = n180pi*atan(v/u);
        switch (q)
        {
    default:
    case 0: break;
    case 1: r = 180.0 - r;
            break;
    case 2: r = -r;
            break;
    case 3: r = -180.0 + r;
            break;
    case 4: r = 90.0 - r;
            break;
    case 5: r = 90.0 + r;
            break;
    case 6: r = -90.0 + r;
            break;
    case 7: r = -90.0 - r;
            break;
        }
    }
    a = make_boxfloat(r, TYPE_DOUBLE_FLOAT);
    errexit();
    return onevalue(a);
}

Lisp_Object Lacos(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(0, a);
}

Lisp_Object Lacosd(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(1, a);
}

Lisp_Object Lacosh(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(2, a);
}

Lisp_Object Lacot(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(3, a);
}

Lisp_Object Lacotd(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(4, a);
}

Lisp_Object Lacoth(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(5, a);
}

Lisp_Object Lacsc(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(6, a);
}

Lisp_Object Lacscd(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(7, a);
}

Lisp_Object Lacsch(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(8, a);
}

Lisp_Object Lasec(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(9, a);
}

Lisp_Object Lasecd(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(10, a);
}

Lisp_Object Lasech(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(11, a);
}

Lisp_Object Lasin(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(12, a);
}

Lisp_Object Lasind(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(13, a);
}

Lisp_Object Lasinh(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(14, a);
}

/* code 15 is for atan */

Lisp_Object Latand(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(16, a);
}

/* code 17 is atan2, 18 is atan2d */

Lisp_Object Latanh(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(19, a);
}

Lisp_Object Lcbrt(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(20, a);
}

Lisp_Object Lcos(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(21, a);
}

Lisp_Object Lcosd(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(22, a);
}

Lisp_Object Lcosh(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(23, a);
}

Lisp_Object Lcot(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(24, a);
}

Lisp_Object Lcotd(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(25, a);
}

Lisp_Object Lcoth(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(26, a);
}

Lisp_Object Lcsc(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(27, a);
}

Lisp_Object Lcscd(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(28, a);
}

Lisp_Object Lcsch(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(29, a);
}

Lisp_Object Lexp(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(30, a);
}

/* 31 is expt, 32 is hypot */

Lisp_Object Lln(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(33, a);
}

/* 34 is 2-arg log */

Lisp_Object Llog10(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(35, a);
}

Lisp_Object Lsec(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(36, a);
}

Lisp_Object Lsecd(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(37, a);
}

Lisp_Object Lsech(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(38, a);
}

Lisp_Object Lsin(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(39, a);
}

Lisp_Object Lsind(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(40, a);
}

Lisp_Object Lsinh(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(41, a);
}

Lisp_Object Lsqrt(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(42, a);
}

Lisp_Object Ltan(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(43, a);
}

Lisp_Object Ltand(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(44, a);
}

Lisp_Object Ltanh(Lisp_Object nil, Lisp_Object a)
{
    CSL_IGNORE(nil);
    return Ltrigfn(45, a);
}

setup_type const arith10_setup[] =
{
    {"abs",                     Labsval, too_many_1, wrong_no_1},
    {"acos",                    Lacos, too_many_1, wrong_no_1},
    {"acosd",                   Lacosd, too_many_1, wrong_no_1},
    {"acosh",                   Lacosh, too_many_1, wrong_no_1},
    {"acot",                    Lacot, too_many_1, wrong_no_1},
    {"acotd",                   Lacotd, too_many_1, wrong_no_1},
    {"acoth",                   Lacoth, too_many_1, wrong_no_1},
    {"acsc",                    Lacsc, too_many_1, wrong_no_1},
    {"acscd",                   Lacscd, too_many_1, wrong_no_1},
    {"acsch",                   Lacsch, too_many_1, wrong_no_1},
    {"asec",                    Lasec, too_many_1, wrong_no_1},
    {"asecd",                   Lasecd, too_many_1, wrong_no_1},
    {"asech",                   Lasech, too_many_1, wrong_no_1},
    {"asin",                    Lasin, too_many_1, wrong_no_1},
    {"asind",                   Lasind, too_many_1, wrong_no_1},
    {"asinh",                   Lasinh, too_many_1, wrong_no_1},
    {"atand",                   Latand, too_many_1, wrong_no_1},
    {"atan2",                   too_few_2, Latan2, wrong_no_2},
    {"atan2d",                  too_few_2, Latan2d, wrong_no_2},
    {"atanh",                   Latanh, too_many_1, wrong_no_1},
    {"cbrt",                    Lcbrt, too_many_1, wrong_no_1},
    {"cos",                     Lcos, too_many_1, wrong_no_1},
    {"cosd",                    Lcosd, too_many_1, wrong_no_1},
    {"cosh",                    Lcosh, too_many_1, wrong_no_1},
    {"cot",                     Lcot, too_many_1, wrong_no_1},
    {"cotd",                    Lcotd, too_many_1, wrong_no_1},
    {"coth",                    Lcoth, too_many_1, wrong_no_1},
    {"csc",                     Lcsc, too_many_1, wrong_no_1},
    {"cscd",                    Lcscd, too_many_1, wrong_no_1},
    {"csch",                    Lcsch, too_many_1, wrong_no_1},
    {"exp",                     Lexp, too_many_1, wrong_no_1},
    {"expt",                    too_few_2, Lexpt, wrong_no_2},
    {"hypot",                   too_few_2, Lhypot, wrong_no_2},
    {"ln",                      Lln, too_many_1, wrong_no_1},
    {"log",                     Lln, Llog_2, wrong_no_2},
    {"log10",                   Llog10, too_many_1, wrong_no_1},
    {"sec",                     Lsec, too_many_1, wrong_no_1},
    {"secd",                    Lsecd, too_many_1, wrong_no_1},
    {"sech",                    Lsech, too_many_1, wrong_no_1},
    {"sin",                     Lsin, too_many_1, wrong_no_1},
    {"sind",                    Lsind, too_many_1, wrong_no_1},
    {"sinh",                    Lsinh, too_many_1, wrong_no_1},
    {"sqrt",                    Lsqrt, too_many_1, wrong_no_1},
    {"tan",                     Ltan, too_many_1, wrong_no_1},
    {"tand",                    Ltand, too_many_1, wrong_no_1},
    {"tanh",                    Ltanh, too_many_1, wrong_no_1},
#ifdef COMMON
    {"cis",                     Lcis, too_many_1, wrong_no_1},
    {"isqrt",                   Lisqrt, too_many_1, wrong_no_1},
    {"phase",                   Lphase, too_many_1, wrong_no_1},
    {"signum",                  Lsignum, too_many_1, wrong_no_1},
    {"atan",                    Latan, Latan_2, wrong_no_1},
#else
    {"atan",                    Latan, too_many_1, wrong_no_1},
    {"logb",                    too_few_2, Llog_2, wrong_no_2},
#endif
    {NULL,                      0, 0, 0}
};

/* end of arith10.c */