Artifact bbdf8ec14feb719af270c6d935769e1ba2e7fa50213c968839842f7f8d87eae5:
- Executable file
r36/cslbase/arith10.c
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 57225) [annotate] [blame] [check-ins using] [more...]
/* arith10.c Copyright (C) 1990-96 Codemist Ltd */ /* * Arithmetic functions. */ /* Signature: 76591977 07-Mar-2000 */ #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 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)>>2)-2]; 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 */