/* 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 */