/* xlmath3 - xlisp math functions modified and augmented to correspond */
/* more closely to Common Lisp standard */
/* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */
/* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */
/* You may give out copies of this software; for conditions see the */
/* file COPYING included with this distribution. */
#include "xlisp.h"
#ifdef BIGNUMS
/* This code supports rational, bignum and complex arithmetic. It is
* based on the following assumptions:
*
* 1) Fixnums use the full range of a FIXTYPE in 2's complement, so
* that overflow of addition and subtraction can be detected by
* checking signs.
*
* 2) The range of a FIXTYPE is at least [-2^31,2^31-1]. Floating
* point products of integers with results in this range are exact and
* are coerced exactly to the right integer by an assignment. For
* systems with larger FITYPE ranges, BIGNUM arithmetic is used for
* larger products and the results are demoted to FIXNUM's
*
* 3) SFIXMIN <= -1 and SFIXMAX >= 2 so that FIXZERO, FIXONE, FIXTWO,
* and FIXNEGTWO can be computed withot allocation or GC. Eventually
* they should be pre-computed in xlinit.c.
*
* The general approach used in this implementation is to code fixnum,
* flonum and complex flonum operations for efficiency. Other
* operations are coded for simplicity, with recursion to simpler
* types used wherever possible.
*
* Internal operations like add2 assume arguments are protected. */
#define IN 0
#define BG 1
#define RT 2
#define FL 3
#define CR 4
#define CF 5
#ifndef MAXFLOFIX
#define MAXFLOFIX 2147483647.0
#endif
#ifndef MINFLOFIX
#define MINFLOFIX -2147483648.0
#endif
#define goodisum(x, y, z) (((x) > 0) ? ((y) < (z)) : ! ((y) < (z)))
#define goodidiff(x, y, z) (!(((x < 0) ^ (y < 0)) && ((z < 0) ^ (x < 0))))
#define infixnumrange(x) (MINFLOFIX <= (x) && (x) <= MAXFLOFIX)
#define FIXZERO cvfixnum((FIXTYPE) 0)
#define FIXONE cvfixnum((FIXTYPE) 1)
#define FIXTWO cvfixnum((FIXTYPE) 2)
#define FIXNEGTWO cvfixnum((FIXTYPE) -2)
#define fixzerop(x) (fixp(x) && getfixnum(x) == 0)
#define fixonep(x) (fixp(x) && getfixnum(x) == 1)
#define realpart(x) (complexp(x) ? getreal(x) : x)
#define imagpart(x) (complexp(x) ? getimag(x) : FIXZERO)
#define numerator(x) (ratiop(x) ? getnumer(x) : x)
#define denominator(x) (ratiop(x) ? getdenom(x) : FIXONE)
#define as_flotype(x) (floatp(x) ? getflonum(x) : makefloat(x))
#define as_bignum(x) \
(bignump(x) ? x : fixp(x) ? cvtfixbignum(getfixnum(x)) : xlbadtype(x))
#define cvcomplex(c) newdcomplex((c).r,(c).i)
#define lss2(x,y) (compare2(x,y) < 0)
#define leq2(x,y) (compare2(x,y) <= 0)
#define geq2(x,y) (compare2(x,y) >= 0)
#define gtr2(x,y) (compare2(x,y) > 0)
#define neq2(x,y) (! equ2(x, y))
/* Function prototypes */
LOCAL VOID badzero(V);
LOCAL double logarithm P3H(FLOTYPE, FLOTYPE, int);
LOCAL int zerop P1H(LVAL);
LOCAL int plusp P1H(LVAL);
LOCAL int minusp P1H(LVAL);
LOCAL int evenp P1H(LVAL);
LOCAL int oddp P1H(LVAL);
LOCAL int commontype P2H(LVAL, LVAL);
LOCAL int compare2 P2H(LVAL, LVAL);
LOCAL int equ2 P2H(LVAL, LVAL);
LOCAL LVAL uminus P1H(LVAL);
LOCAL LVAL add2 P2H(LVAL, LVAL);
LOCAL LVAL sub2 P2H(LVAL, LVAL);
LOCAL LVAL mul2 P2H(LVAL, LVAL);
LOCAL LVAL div2 P2H(LVAL, LVAL);
#ifdef XLISP_STAT
LOCAL LVAL fdiv2 P2H(LVAL, LVAL);
#endif /* XLISP_STAT */
LOCAL LVAL min2 P2H(LVAL, LVAL);
LOCAL LVAL max2 P2H(LVAL, LVAL);
LOCAL VOID rational_idiv2 P5H(LVAL, LVAL, LVAL *, LVAL *, int);
LOCAL VOID flonum_idiv2 P5H(FLOTYPE, FLOTYPE, FLOTYPE *, FLOTYPE *, int);
LOCAL VOID idiv2 P5H(LVAL, LVAL, LVAL *, LVAL *, int);
LOCAL LVAL add1 P1H(LVAL);
LOCAL LVAL sub1 P1H(LVAL);
LOCAL LVAL gcd2 P2H(LVAL, LVAL);
#ifdef applec
/* a gross hack */
VOID MODF(double x, double *y)
{
extended value;
modf(x, &value);
*y = value;
}
#else
#define MODF modf
#endif /* applec */
/* Error checking for illegal zero argument */
LOCAL VOID badzero(V) { xlfail("illegal zero argument"); }
#define checkzero(x) if (zerop(x)) badzero();
#define checkizero(x) if (x == 0) badzero();
#define checkfzero(x) if (x == 0.0) badzero();
/* Complex number functions */
double d_sign P2C(double *, a, double *, b)
{
/* from the f2c distribution */
double x;
x = *a >= 0 ? *a : - *a;
return *b >= 0 ? x : -x;
}
double z_abs P1C(dcomplex *, z)
{
/* from the f2c distribution */
double real, imag;
double temp;
real = z->r;
imag = z->i;
if(real < 0)
real = -real;
if(imag < 0)
imag = -imag;
if(imag > real){
temp = real;
real = imag;
imag = temp;
}
if((real+imag) == real)
return(real);
temp = imag/real;
temp = real*sqrt(1.0 + temp*temp); /*overflow!!*/
return(temp);
}
VOID z_div P3C(dcomplex *, c, dcomplex *, a, dcomplex *, b)
{
/* from the f2c distribution */
double ratio, den;
double abr, abi;
dcomplex r;
if( (abr = b->r) < 0.)
abr = - abr;
if( (abi = b->i) < 0.)
abi = - abi;
if( abr <= abi ) {
if(abi == 0)
xlfail("complex division by zero");
ratio = b->r / b->i ;
den = b->i * (1 + ratio*ratio);
r.r = (a->r*ratio + a->i) / den;
r.i = (a->i*ratio - a->r) / den;
}
else {
ratio = b->i / b->r ;
den = b->r * (1 + ratio*ratio);
r.r = (a->r + a->i*ratio) / den;
r.i = (a->i - a->r*ratio) / den;
}
c->r = r.r, c->i = r.i;
}
VOID z_sqrt P2C(dcomplex *, r, dcomplex *, z)
{
/* from the f2c distribution */
double mag;
dcomplex v;
if((mag = z_abs(z)) == 0.0)
v.r = v.i = 0.0;
else if(z->r > 0) {
v.r = sqrt(0.5 * (mag + z->r));
v.i = z->i / v.r / 2;
}
else {
v.i = sqrt(0.5 * (mag - z->r));
if(z->i < 0)
v.i = - v.i;
v.r = z->i / v.i / 2;
}
r->r = v.r, r->i = v.i;
}
LOCAL VOID makecomplex P2C(dcomplex *, c, LVAL, x)
{
if (realp(x)) {
c->r = makefloat(x);
c->i = 0.0;
}
else if (complexp(x)) {
c->r = makefloat(getreal(x));
c->i = makefloat(getimag(x));
}
else xlerror("not a number", x);
}
LOCAL double z_phase P1C(dcomplex *, c)
{
double phi;
if (c->r == 0.0) {
if (c->i > 0.0) phi = PI / 2;
else if (c->i == 0.0) phi = 0.0;
else phi = -PI / 2;
}
else phi = atan2(c->i, c->r);
return(phi);
}
LOCAL VOID z_exp P2C(dcomplex *, r, dcomplex *, z)
{
/* from the f2c distribution */
double expx;
expx = exp(z->r);
r->r = expx * cos(z->i);
r->i = expx * sin(z->i);
}
LOCAL VOID z_log P2C(dcomplex *, r, dcomplex *, z)
{
/* adapted from the f2c distribution */
double mod;
dcomplex v;
mod = z_abs(z);
checkfzero(mod);
v.r = log(mod);
v.i = atan2(z->i, z->r);
r->r = v.r, r->i = v.i;
}
LOCAL VOID z_expt P3C(dcomplex *, r, dcomplex *, a, dcomplex *, b)
{
/* adapted from pow_zz.c in the f2c distribution */
if (b->r == 0.0 && b->i == 0.0)
r->r = 1.0, r->i = 0.0;
else if (a->r == 0.0 && a->i == 0.0)
r->r = 0.0, r->i = 0.0;
else {
double logr, logi, x, y;
logr = log(z_abs(a));
logi = atan2(a->i, a->r);
x = exp(logr * b->r - logi * b->i);
y = logr * b->i + logi * b->r;
r->r = x * cos(y);
r->i = x * sin(y);
}
}
LOCAL VOID z_sin P2C(dcomplex *, r, dcomplex *, c)
{
/* Better code by Hume Smith, 11/92 */
dcomplex val;
val.r = sin(c->r) * cosh(c->i);
val.i = cos(c->r) * sinh(c->i);
r->r = val.r, r->i = val.i;
}
LOCAL VOID z_cos P2C(dcomplex *, r, dcomplex *, c)
{
/* Better code by Hume Smith, 11/92 */
dcomplex val;
val.r = cos(c->r) * cosh(c->i);
val.i = -sin(c->r) * sinh(c->i);
r->r = val.r, r->i = val.i;
}
LOCAL VOID z_tan P2C(dcomplex *, ret_val, dcomplex *, z)
{
/* Adapted from ctan.f in FN from NETLIB. */
/* This function does not attempt to signal loss of precision. */
double x2, y2;
double den, sn2x;
x2 = z->r * 2.0;
y2 = z->i * 2.0;
sn2x = sin(x2);
den = cos(x2) + cosh(y2);
checkfzero(den);
ret_val->r = sn2x / den, ret_val->i = sinh(y2) / den;
return;
}
LOCAL VOID z_asin P2C(dcomplex *, r, dcomplex *, c)
{
/* Improved code by Hume Smith 11/92 */
dcomplex sx, val;
/* compute sx = 1 - c*c */
sx.r = 1.0 - c->r*c->r + c->i*c->i;
sx.i = -2.0*c->r*c->i;
/* compute sx = i*c + sqrt(sx) */
z_sqrt(&sx, &sx);
sx.r -= c->i;
sx.i += c->r;
/* compute val = -i ln(sx) */
z_log(&sx, &sx);
val.r = sx.i;
val.i = -sx.r;
r->r = val.r, r->i = val.i;
}
LOCAL VOID z_acos P2C(dcomplex *, r, dcomplex *, c)
{
/* Improved code by Hume Smith 11/92 */
dcomplex sx, val;
/* compute val = 1 - c*c */
sx.r = 1.0 - c->r*c->r + c->i*c->i;
sx.i = -2.0*c->r*c->i;
z_sqrt(&val, &sx);
/* compute sx = c + i sqrt(val) */
sx.r = c->r - val.i;
sx.i = c->i + val.r;
z_log(&sx, &sx);
/* compute val = -i ln(sx)) */
val.r = sx.i;
val.i = -sx.r;
r->r = val.r, r->i = val.i;
}
LOCAL VOID z_atan P2C(dcomplex *, val, dcomplex *, c)
{
/* This has been redefined in Jan 1989 */
dcomplex sx, ix, v;
sx.r = 1.0 - c->i, sx.i = c->r;
ix.r = 1.0 + c->i, ix.i = -c->r;
z_log(&sx, &sx);
z_log(&ix, &ix);
sx.r -= ix.r, sx.i -= ix.i;
ix.r = 0.0, ix.i = 2.0;
z_div(&v, &sx, &ix);
val->r = v.r, val->i = v.i;
}
LOCAL VOID z_atan2 P3C(dcomplex *, val, dcomplex *, n, dcomplex *, d)
{
dcomplex v;
/* adapted from catan2.f in FN from NETLIB */
if (z_abs(d) != 0.0) {
z_div(&v, n, d);
z_atan(&v, &v);
}
else {
v.i = 0.0;
if (z_abs(n) != 0.0) {
double d__1, d__2;
d__1 = PI * 0.5;
d__2 = n->r;
v.r = d_sign(&d__1, &d__2), v.i = 0.0;
}
else
v.r = 0.0;
}
val->r = v.r, val->i = v.i;
}
/* Internal generic arithmetic functions */
LOCAL int zerop P1C(LVAL, x)
{
switch (ntype(x)) {
case FIXNUM: return (getfixnum(x) == 0);
case RATIO: return FALSE;
case BIGNUM: return zeropbignum(x);
case FLONUM: return (getflonum(x) == 0.0);
case COMPLEX: return (zerop(getreal(x)) && zerop(getimag(x)));
default: return FALSE;
}
}
LOCAL int plusp P1C(LVAL, x)
{
switch (ntype(x)) {
case FIXNUM: return (getfixnum(x) > 0);
case RATIO: return plusp(getnumer(x));
case BIGNUM: return ! getbignumsign(x);
case FLONUM: return (getflonum(x) > 0.0);
default: return FALSE;
}
}
LOCAL int minusp P1C(LVAL, x)
{
switch (ntype(x)) {
case FIXNUM: return (getfixnum(x) < 0);
case RATIO: return minusp(getnumer(x));
case BIGNUM: return getbignumsign(x);
case FLONUM: return (getflonum(x) < 0.0);
default: return FALSE;
}
}
LOCAL int evenp P1C(LVAL, x)
{
switch (ntype(x)) {
case FIXNUM:
return (getfixnum(x) % 2 == 0);
case BIGNUM:
{
FIXTYPE temp;
LVAL bigtwo, r;
xlsave1(bigtwo);
bigtwo = cvtfixbignum((FIXTYPE) 2);
divbignum(x, bigtwo, &r);
xlpop();
return (cvtbigfixnum(r, &temp) && temp == 0);
}
default: return FALSE;
}
}
LOCAL int oddp P1C(LVAL, x) { return integerp(x) && ! evenp(x); }
LOCAL int commontype P2C(LVAL, x, LVAL, y)
{
switch (ntype(x)) {
case FIXNUM:
switch (ntype(y)) {
case FIXNUM: return IN;
case BIGNUM: return BG;
case RATIO: return RT;
case FLONUM: return FL;
case COMPLEX: return floatp(getreal(y)) ? CF : CR;
default: xlbadtype(y);
}
case BIGNUM:
switch (ntype(y)) {
case FIXNUM:
case BIGNUM: return BG;
case RATIO: return RT;
case FLONUM: return FL;
case COMPLEX: return floatp(getreal(y)) ? CF : CR;
default: xlbadtype(y);
}
case RATIO:
switch (ntype(y)) {
case FIXNUM:
case BIGNUM:
case RATIO: return RT;
case FLONUM: return FL;
case COMPLEX: return floatp(getreal(y)) ? CF : CR;
default: xlbadtype(y);
}
case FLONUM:
switch (ntype(y)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM: return FL;
case COMPLEX: return CF;
default: xlbadtype(y);
}
case COMPLEX:
switch (ntype(y)) {
case FIXNUM:
case BIGNUM:
case RATIO: return floatp(getreal(x)) ? CF : CR;
case FLONUM: return CF;
case COMPLEX: return (floatp(getreal(x)) || floatp(getreal(y))) ? CF : CR;
default: xlbadtype(y);
}
default: xlbadtype(x);
}
return 0; /* to keep compilers happy */
}
LOCAL int compare2 P2C(LVAL, x, LVAL, y)
{
switch (commontype(x, y)) {
case IN:
{
FIXTYPE xi, yi;
xi = getfixnum(x);
yi = getfixnum(y);
return xi < yi ? -1 : xi > yi ? 1 : 0;
}
case BG:
{
int val;
xlstkcheck(2);
xlprotect(x);
xlprotect(y);
x = as_bignum(x);
y = as_bignum(y);
val = comparebignum(x, y);
xlpopn(2);
return val;
}
case RT:
{
LVAL nx, dx, ny, dy, a, b;
int val;
xlstkcheck(2);
xlsave(a);
xlsave(b);
nx = numerator(x); dx = denominator(x);
ny = numerator(y); dy = denominator(y);
a = mul2(nx, dy); b = mul2(ny, dx);
val = compare2(a, b);
xlpopn(2);
return val;
}
case FL:
{
FLOTYPE xd, yd;
xd = as_flotype(x);
yd = as_flotype(y);
return xd < yd ? -1 : xd > yd ? 1 : 0;
}
default:
xlbadtype(complexp(x) ? x : y);
return 0; /* keep compilers happy */
}
}
LOCAL int equ2 P2C(LVAL, x, LVAL, y)
{
switch (commontype(x, y)) {
case IN:
case BG:
case RT:
case FL:
return compare2(x, y) == 0;
case CR:
return equ2(realpart(x),realpart(y)) && equ2(imagpart(x),imagpart(y));
case CF:
{
LVAL rx, ix, ry, iy;
FLOTYPE rxd, ixd, ryd, iyd;
rx = realpart(x); rxd = as_flotype(rx);
ix = imagpart(x); ixd = as_flotype(ix);
ry = realpart(y); ryd = as_flotype(ry);
iy = imagpart(y); iyd = as_flotype(iy);
return rxd == ryd && ixd == iyd;
}
default: return 0; /* keep compilers happy */
}
}
LOCAL LVAL uminus P1C(LVAL, x)
{
switch (ntype(x)) {
case FIXNUM:
{
FIXTYPE ix;
ix = getfixnum(x);
if (ix > MINFIX)
return cvfixnum(-ix);
/* else fall through */
}
case BIGNUM:
{
FIXTYPE temp;
xlprot1(x);
x = as_bignum(x);
x = copybignum(x, ! getbignumsign(x));
xlpop();
return cvtbigfixnum(x, &temp) ? cvfixnum(temp) : x;
}
case RATIO:
{
LVAL n, y;
xlsave1(n);
n = getnumer(x); n = uminus(n);
y = div2(n, getdenom(x));
xlpop();
return y;
}
case FLONUM:
return cvflonum(-getflonum(x));
case COMPLEX:
if (floatp(getreal(x))) {
return newdcomplex(-getflonum(getreal(x)), -getflonum(getimag(x)));
}
else {
LVAL r, i, y;
xlstkcheck(2);
xlsave(r);
xlsave(i);
r = getreal(x); i = getimag(x);
r = uminus(r); i = uminus(i);
y = newcomplex(r, i);
xlpopn(2);
return y;
}
default: return xlbadtype(x);
}
}
LOCAL LVAL add2 P2C(LVAL, x, LVAL, y)
{
/* try for quick return */
if (fixzerop(x))
return y;
if (fixzerop(y))
return x;
/* general case */
switch (commontype(x, y)) {
case IN:
{
FIXTYPE xi, yi, zi;
xi = getfixnum(x);
yi = getfixnum(y);
zi = xi + yi;
if (goodisum(xi, yi, zi))
return cvfixnum(zi);
/* else drop through */
}
case BG:
{
LVAL z;
FIXTYPE temp;
xlstkcheck(2);
xlprotect(x);
xlprotect(y);
x = as_bignum(x);
y = as_bignum(y);
z = addsubbignum(x, y, FALSE);
xlpopn(2);
return cvtbigfixnum(z, &temp) ? cvfixnum(temp) : z;
}
case RT:
{
LVAL nx, dx, ny, dy, n, d, z;
xlstkcheck(4);
xlsave(nx);
xlsave(ny);
xlsave(n);
xlsave(d);
nx = numerator(x); dx = denominator(x);
ny = numerator(y); dy = denominator(y);
nx = mul2(nx, dy);
ny = mul2(ny, dx);
n = add2(nx, ny);
d = mul2(dx, dy);
z = div2(n, d);
xlpopn(4);
return z;
}
case FL:
{
FLOTYPE xd, yd;
xd = as_flotype(x);
yd = as_flotype(y);
return cvflonum(xd + yd);
}
case CR:
{
LVAL rx, ix, ry, iy, rz, iz, z;
xlstkcheck(2);
xlsave(rz);
xlsave(iz);
rx = realpart(x); ry = realpart(y); rz = add2(rx, ry);
ix = imagpart(x); iy = imagpart(y); iz = add2(ix, iy);
z = newcomplex(rz, iz);
xlpopn(2);
return z;
}
case CF:
{
LVAL rx, ix, ry, iy;
FLOTYPE rxd, ixd, ryd, iyd;
rx = realpart(x); rxd = as_flotype(rx);
ix = imagpart(x); ixd = as_flotype(ix);
ry = realpart(y); ryd = as_flotype(ry);
iy = imagpart(y); iyd = as_flotype(iy);
return newdcomplex(rxd + ryd, ixd + iyd);
}
default: return NIL; /* to keep compiler happy */
}
}
LOCAL LVAL sub2 P2C(LVAL, x, LVAL, y)
{
/* try for quick exit */
if (fixzerop(y))
return x;
if (fixzerop(x))
return uminus(y);
/* general case */
switch (commontype(x, y)) {
case IN:
{
FIXTYPE xi, yi, zi;
xi = getfixnum(x);
yi = getfixnum(y);
zi = xi - yi;
if (goodidiff(xi, yi, zi))
return cvfixnum(zi);
/* else drop through */
}
case BG:
{
LVAL z;
FIXTYPE temp;
xlstkcheck(2);
xlprotect(x);
xlprotect(y);
x = as_bignum(x);
y = as_bignum(y);
z = addsubbignum(x, y, TRUE);
xlpopn(2);
return cvtbigfixnum(z, &temp) ? cvfixnum(temp) : z;
}
case RT:
{
LVAL nx, dx, ny, dy, n, d, z;
xlstkcheck(4);
xlsave(nx);
xlsave(ny);
xlsave(n);
xlsave(d);
nx = numerator(x); dx = denominator(x);
ny = numerator(y); dy = denominator(y);
nx = mul2(nx, dy);
ny = mul2(ny, dx);
n = sub2(nx, ny);
d = mul2(dx, dy);
z = div2(n, d);
xlpopn(4);
return z;
}
case FL:
{
FLOTYPE xd, yd;
xd = as_flotype(x);
yd = as_flotype(y);
return cvflonum(xd - yd);
}
case CR:
{
LVAL rx, ix, ry, iy, rz, iz, z;
xlstkcheck(2);
xlsave(rz);
xlsave(iz);
rx = realpart(x); ry = realpart(y); rz = sub2(rx, ry);
ix = imagpart(x); iy = imagpart(y); iz = sub2(ix, iy);
z = newcomplex(rz, iz);
xlpopn(2);
return z;
}
case CF:
{
LVAL rx, ix, ry, iy;
FLOTYPE rxd, ixd, ryd, iyd;
rx = realpart(x); rxd = as_flotype(rx);
ix = imagpart(x); ixd = as_flotype(ix);
ry = realpart(y); ryd = as_flotype(ry);
iy = imagpart(y); iyd = as_flotype(iy);
return newdcomplex(rxd - ryd, ixd - iyd);
}
default: return NIL; /* to keep compiler happy */
}
}
LOCAL LVAL mul2 P2C(LVAL, x, LVAL, y)
{
/* try for quick exit */
if (fixonep(x))
return y;
if (fixonep(y))
return x;
if (fixzerop(x) && (rationalp(y) || (complexp(y) && rationalp(getreal(y)))))
return FIXZERO;
if (fixzerop(y) && (rationalp(x) || (complexp(x) && rationalp(getreal(x)))))
return FIXZERO;
/* general case */
switch (commontype(x, y)) {
case IN:
{
double xd, yd, zd;
xd = (double) getfixnum(x);
yd = (double) getfixnum(y);
zd = xd * yd;
if (infixnumrange(zd))
return cvfixnum((FIXTYPE) zd);
/* else drop through */
}
case BG:
{
LVAL z;
FIXTYPE temp;
xlstkcheck(2);
xlprotect(x);
xlprotect(y);
x = as_bignum(x);
y = as_bignum(y);
z = multbignum(x, y);
xlpopn(2);
return cvtbigfixnum(z, &temp) ? cvfixnum(temp) : z;
}
case RT:
{
LVAL nx, dx, ny, dy, n, d, z;
xlstkcheck(2);
xlsave(n);
xlsave(d);
nx = numerator(x); dx = denominator(x);
ny = numerator(y); dy = denominator(y);
n = mul2(nx, ny);
d = mul2(dx, dy);
z = div2(n, d);
xlpopn(2);
return z;
}
case FL:
{
FLOTYPE xd, yd;
xd = as_flotype(x);
yd = as_flotype(y);
return cvflonum(xd * yd);
}
case CR:
{
LVAL rx, ix, ry, iy, tmp1, tmp2, rz, iz, z;
xlstkcheck(4);
xlsave(tmp1);
xlsave(tmp2);
xlsave(rz);
xlsave(iz);
rx = realpart(x); ix = imagpart(x);
ry = realpart(y); iy = imagpart(y);
tmp1 = mul2(rx, ry); tmp2 = mul2(ix, iy); rz = sub2(tmp1, tmp2);
tmp1 = mul2(ix, ry); tmp2 = mul2(rx, iy); iz = add2(tmp1, tmp2);
z = newcomplex(rz, iz);
xlpopn(4);
return z;
}
case CF:
{
LVAL rx, ix, ry, iy;
FLOTYPE rxd, ixd, ryd, iyd;
rx = realpart(x); rxd = as_flotype(rx);
ix = imagpart(x); ixd = as_flotype(ix);
ry = realpart(y); ryd = as_flotype(ry);
iy = imagpart(y); iyd = as_flotype(iy);
return newdcomplex(rxd * ryd - ixd * iyd, ixd * ryd + rxd * iyd);
}
default: return NIL; /* to keep compiler happy */
}
}
LOCAL LVAL div2 P2C(LVAL, x, LVAL, y)
{
checkzero(y);
/* try quick exit */
if (fixonep(y))
return x;
if (fixzerop(x) && (rationalp(y) || (complexp(y) && rationalp(getreal(y)))))
return FIXZERO;
/* general case */
switch (commontype(x, y)) {
case IN:
return cvratio(getfixnum(x), getfixnum(y));
case BG:
{
LVAL z;
xlstkcheck(2);
xlprotect(x);
xlprotect(y);
x = as_bignum(x);
y = as_bignum(y);
z = cvbratio(x, y);
xlpopn(2);
return z;
}
case RT:
{
LVAL nx, dx, ny, dy, n, d, z;
xlstkcheck(2);
xlsave(n);
xlsave(d);
nx = numerator(x); dx = denominator(x);
ny = numerator(y); dy = denominator(y);
n = mul2(nx, dy);
d = mul2(dx, ny);
z = div2(n, d);
xlpopn(2);
return z;
}
case FL:
{
FLOTYPE xd, yd;
xd = as_flotype(x);
yd = as_flotype(y);
return cvflonum(xd / yd);
}
case CR:
{
LVAL rx, ix, ry, iy, tmp1, tmp2, rz, iz, d, z;
xlstkcheck(5);
xlsave(tmp1);
xlsave(tmp2);
xlsave(d);
xlsave(rz);
xlsave(iz);
rx = realpart(x); ix = imagpart(x);
ry = realpart(y); iy = imagpart(y);
tmp1 = mul2(ry, ry); tmp2 = mul2(iy, iy); d = add2(tmp1, tmp2);
tmp1 = mul2(rx, ry); tmp2 = mul2(ix, iy); tmp1 = add2(tmp1, tmp2);
rz = div2(tmp1, d);
tmp1 = mul2(ix, ry); tmp2 = mul2(rx, iy); tmp1 = sub2(tmp1, tmp2);
iz = div2(tmp1, d);
z = newcomplex(rz, iz);
xlpopn(5);
return z;
}
case CF:
{
LVAL rx, ix, ry, iy;
dcomplex cx, cy, cz;
rx = realpart(x); cx.r = as_flotype(rx);
ix = imagpart(x); cx.i = as_flotype(ix);
ry = realpart(y); cy.r = as_flotype(ry);
iy = imagpart(y); cy.i = as_flotype(iy);
z_div(&cz, &cx, &cy);
return newdcomplex(cz.r, cz.i);
}
default: return NIL; /* to keep compiler happy */
}
}
#ifdef XLISP_STAT
LOCAL LVAL fdiv2 P2C(LVAL, x, LVAL, y)
{
switch (commontype(x, y)) {
case IN:
{
FIXTYPE xi, yi;
xi = getfixnum(x);
yi = getfixnum(y);
checkizero(yi);
if (yi == -1 && xi == MINFIX)
return cvratio(xi, yi);
else if (xi % yi == 0)
return cvfixnum(xi / yi);
else
return cvflonum(xi / (FLOTYPE) yi);
}
case BG:
case RT:
{
LVAL z;
z = div2(x, y);
if (ratiop(z))
z = cvflonum(makefloat(z));
return z;
}
case FL:
{
FLOTYPE xd, yd;
xd = as_flotype(x);
yd = as_flotype(y);
checkfzero(yd);
return cvflonum(xd / yd);
}
case CR:
{
LVAL z, rz, iz;
z = div2(x,y);
rz = realpart(z); iz = imagpart(z);
if (ratiop(rz) || ratiop(iz))
z = newdcomplex(makefloat(rz), makefloat(iz));
return z;
}
case CF:
{
LVAL rx, ix, ry, iy;
dcomplex cx, cy, cz;
rx = realpart(x); cx.r = as_flotype(rx);
ix = imagpart(x); cx.i = as_flotype(ix);
ry = realpart(y); cy.r = as_flotype(ry);
iy = imagpart(y); cy.i = as_flotype(iy);
z_div(&cz, &cx, &cy);
return newdcomplex(cz.r, cz.i);
}
default: return NIL; /* to keep compiler happy */
}
}
#endif /* XLISP_STAT */
LOCAL LVAL min2 P2C(LVAL, x, LVAL, y)
{
switch (commontype(x, y)) {
case IN: return (getfixnum(x) <= getfixnum(y)) ? x : y;
case FL: return (as_flotype(x) <= as_flotype(y)) ? x : y;
default: return leq2(x, y) ? x : y;
}
}
LOCAL LVAL max2 P2C(LVAL, x, LVAL, y)
{
switch (commontype(x, y)) {
case IN: return (getfixnum(x) >= getfixnum(y)) ? x : y;
case FL: return (as_flotype(x) >= as_flotype(y)) ? x : y;
default: return geq2(x, y) ? x : y;
}
}
LOCAL VOID rational_idiv2 P5C(LVAL, x, LVAL, y,
LVAL *, pn, LVAL *, pr,
int, which)
{
LVAL n, r;
checkzero(y);
if (fixonep(y) && integerp(x)) { /* check for quick exit */
n = x;
r = FIXZERO;
}
else {
xlstkcheck(2);
xlsave(n);
xlsave(r);
if (integerp(x) && integerp(y)) { /* base case */
LVAL xb, yb;
FIXTYPE temp;
xlstkcheck(2);
xlsave(xb);
xlsave(yb);
xb = as_bignum(x);
yb = as_bignum(y);
n = divbignum(xb, yb, &r);
xlpopn(2);
if (cvtbigfixnum(n, &temp))
n = cvfixnum(temp);
if (cvtbigfixnum(r, &temp))
r = cvfixnum(temp);
}
else {
LVAL xn, xd, yn, yd, tmp1, tmp2;
xn = numerator(x); xd = denominator(x);
yn = numerator(y); yd = denominator(y);
xlstkcheck(2);
xlsave(tmp1);
xlsave(tmp2);
tmp1 = mul2(xn, yd);
tmp2 = mul2(yn, xd);
rational_idiv2(tmp1, tmp2, &n, &r, 'I');
tmp1 = mul2(xd, yd);
r = div2(r, tmp1);
xlpopn(2);
}
/* adjust ratio and remainder if necessary */
if (which != 'I' && ! zerop(r)) {
switch (which) {
case '^':
if ((plusp(r) && plusp(y)) || (minusp(r) && minusp(y))) {
n = add1(n);
r = sub2(r, y);
}
break;
case '_':
if ((minusp(r) && plusp(y)) || (plusp(r) && minusp(y))) {
n = sub1(n);
r = add2(r, y);
}
break;
case 'r': /* rounds away from zero on ties, not towards even number */
{
LVAL cut;
xlsave1(cut);
if (plusp(y)) {
if (plusp(r)) {
cut = div2(y, FIXTWO);
if (geq2(r, cut)) { n = add1(n); r = sub2(r, y); }
}
else {
cut = div2(y, FIXNEGTWO);
if (leq2(r, cut)) { n = sub1(n); r = add2(r, y); }
}
}
else {
if (plusp(r)) {
cut = div2(y, FIXNEGTWO);
if (geq2(r, cut)) { n = sub1(n); r = add2(r, y); }
}
else {
cut = div2(y, FIXTWO);
if (leq2(r, cut)) { n = add1(n); r = sub2(r, y); }
}
}
xlpop();
}
break;
}
}
xlpopn(2);
}
*pn = n;
*pr = r;
}
LOCAL VOID flonum_idiv2 P5C(FLOTYPE, xd, FLOTYPE, yd,
FLOTYPE *, pnd, FLOTYPE *, prd,
int, which)
{
FLOTYPE nd = *pnd, rd = *prd;
checkfzero(yd);
MODF(xd / yd, &nd);
rd = xd - nd * yd;
if (rd != 0.0) {
switch (which) {
case '^':
if ((rd > 0.0 && yd > 0.0) || (rd < 0.0 && yd < 0.0)) {
nd++;
rd = rd - yd;
}
break;
case '_':
if ((rd < 0.0 && yd > 0.0) || (rd > 0.0 && yd < 0.0)) {
nd--;
rd = rd + yd;
}
break;
case 'r': /* rounds away from zero on ties, not towards even number */
if (yd > 0.0) {
if (rd > 0.0) {
if (rd >= 0.5 * yd) { nd++; rd -= yd; }
}
else {
if (rd <= -0.5 * yd) { nd--; rd += yd; }
}
}
else {
if (rd > 0.0) {
if (rd >= -0.5 * yd) { nd--; rd += yd; }
}
else {
if (rd <= 0.5 * yd) { nd++; rd -= yd; }
}
}
break;
}
}
*pnd = nd;
*prd = rd;
}
LOCAL VOID idiv2 P5C(LVAL, x, LVAL, y, LVAL *, pn, LVAL *, pr, int, which)
{
LVAL n, r;
switch (commontype(x, y)) {
case IN:
{
FIXTYPE xi, yi;
xi = getfixnum(x);
yi = getfixnum(y);
if (yi == 1) { /* check for quick exit */
n = x;
r = FIXZERO;
break;
}
else if (xi > MINFIX && xi < MAXFIX && yi > MINFIX && yi < MAXFIX) {
FIXTYPE ni, ri;
checkizero(yi);
ni = xi / yi;
ri = xi % yi;
/* adjust ratio and remainder if necessary */
if (ri != 0) {
switch (which) {
case '^':
if ((ri > 0 && yi > 0) || (ri < 0 && yi < 0)) {
ni++;
ri = ri - yi;
}
break;
case '_':
if ((ri < 0 && yi > 0) || (ri > 0 && yi < 0)) {
ni--;
ri = ri + yi;
}
break;
case 'r':/* rounds away from zero on ties, not towards even number */
if (yi > 0) {
if (ri > 0) {
if (ri >= (yi + 1) / 2) { ni++; ri -= yi; }
}
else {
if (ri <= (-yi - 1) / 2) { ni--; ri += yi; }
}
}
else {
if (ri > 0) {
if (ri >= (-yi + 1) / 2) { ni--; ri += yi; }
}
else {
if (ri <= (yi - 1) / 2) { ni++; ri -= yi; }
}
}
break;
}
}
/* only allocate if actually needed */
if (pn) {
xlsave1(n); /* protect n while allocating r */
n = cvfixnum(ni);
r = cvfixnum(ri);
xlpop();
}
else {
n = NIL;
r = cvfixnum(ri);
}
break;
}
/* else drop through */
}
case BG:
case RT:
rational_idiv2(x, y, &n, &r, which);
break;
case FL:
{
FLOTYPE xd, yd, nd, rd;
xd = as_flotype(x);
yd = as_flotype(y);
flonum_idiv2(xd, yd, &nd, &rd, which);
/* only allocate n if actually needed */
if (pn) {
xlsave1(n); /* protect n while allocating r */
if (infixnumrange(nd))
n = cvfixnum((FIXTYPE) nd);
else {
FIXTYPE temp;
n = cvtflobignum(nd);
if (cvtbigfixnum(n, &temp))
n = cvfixnum(temp); /* not needed for 32-big FIXTYPE's */
}
r = cvflonum(rd);
xlpop();
}
else {
n = NIL;
r = cvflonum(rd);
}
break;
}
default:
xlbadtype(complexp(x) ? x : y);
r = n = NIL; /* fool compiler into not giving warning */
}
if (pn) *pn = n;
if (pr) *pr = r;
}
LOCAL LVAL add1 P1C(LVAL, x)
{
switch (ntype(x)) {
case FIXNUM:
{
FIXTYPE xi, yi;
xi = getfixnum(x);
yi = xi + 1;
if (goodisum(xi, 1, yi))
return cvfixnum(yi);
/* else drop through */
}
case BIGNUM:
{
LVAL y;
FIXTYPE temp;
xlprot1(x);
x = as_bignum(x);
y = addsubbignum(x, n_bigmone, TRUE);
y = cvtbigfixnum(y, &temp) ? cvfixnum(temp) : y;
xlpop();
return y;
}
case RATIO:
{
LVAL nx, dx, y;
xlsave1(nx);
nx = numerator(x); dx = denominator(x);
nx = add2(nx, dx);
y = div2(nx, dx);
xlpop();
return y;
}
case FLONUM:
return cvflonum(as_flotype(x) + 1.0);
case COMPLEX:
if (rationalp(getreal(x)) && rationalp(getimag(x))) {
LVAL ry, y;
xlsave1(ry);
ry = add2(realpart(x), FIXONE);
y = newcomplex(ry, imagpart(x));
xlpop();
return y;
}
else {
LVAL rx, ix;
FLOTYPE ixd, rxd;
rx = realpart(x); rxd = as_flotype(rx);
ix = imagpart(x); ixd = as_flotype(ix);
return newdcomplex(rxd + 1.0, ixd);
}
default: return NIL; /* to keep compiler happy */
}
}
LOCAL LVAL sub1 P1C(LVAL, x)
{
switch (ntype(x)) {
case FIXNUM:
{
FIXTYPE xi, yi;
xi = getfixnum(x);
yi = xi - 1;
if (goodidiff(xi, 1, yi))
return cvfixnum(yi);
/* else drop through */
}
case BIGNUM:
{
LVAL y;
FIXTYPE temp;
xlprot1(x);
x = as_bignum(x);
y = addsubbignum(x, n_bigmone, FALSE);
y = cvtbigfixnum(y, &temp) ? cvfixnum(temp) : y;
xlpop();
return y;
}
case RATIO:
{
LVAL nx, dx, y;
xlsave1(nx);
nx = numerator(x); dx = denominator(x);
nx = sub2(nx, dx);
y = div2(nx, dx);
xlpop();
return y;
}
case FLONUM:
return cvflonum(as_flotype(x) - 1.0);
case COMPLEX:
if (rationalp(getreal(x)) && rationalp(getimag(x))) {
LVAL ry, y;
xlsave1(ry);
ry = sub2(realpart(x), FIXONE);
y = newcomplex(ry, imagpart(x));
xlpop();
return y;
}
else {
LVAL rx, ix;
FLOTYPE ixd, rxd;
rx = realpart(x); rxd = as_flotype(rx);
ix = imagpart(x); ixd = as_flotype(ix);
return newdcomplex(rxd - 1.0, ixd);
}
default: return NIL; /* to keep compiler happy */
}
}
LOCAL LVAL gcd2 P2C(LVAL, n, LVAL, m) /* euclid's algorith */
{
LVAL r;
xlstkcheck(3);
xlprotect(n);
xlprotect(m);
xlsave(r);
for (;;) {
idiv2(m, n, NULL, &r, 'I');
if (zerop(r))
break;
m = n;
n = r;
}
xlpopn(3);
return (n);
}
LVAL xadd(V)
{
LVAL x, y;
xlsave1(x);
x = moreargs() ? xlganumber() : FIXZERO;
while (moreargs()) {
y = xlgetarg();
x = add2(x, y);
}
xlpop();
return x;
}
LVAL xsub(V)
{
LVAL x, y;
x = xlgetarg();
if (moreargs()) {
xlprot1(x);
do {
y = xlgetarg();
x = sub2(x, y);
} while (moreargs());
xlpop();
}
else
x = uminus(x);
return x;
}
LVAL xmul(V)
{
LVAL x, y;
xlsave1(x);
x = moreargs() ? xlganumber() : FIXONE;
while (moreargs()) {
y = xlgetarg();
x = mul2(x, y);
}
xlpop();
return x;
}
LVAL xdiv(V)
{
LVAL x, y;
x = xlgetarg();
if (moreargs()) {
xlprot1(x);
do {
y = xlgetarg();
x = div2(x, y);
} while (moreargs());
xlpop();
}
else
x = div2(FIXONE, x);
return x;
}
#ifdef XLISP_STAT
LVAL xfdiv(V)
{
LVAL x, y;
x = xlgetarg();
if (moreargs()) {
xlprot1(x);
do {
y = xlgetarg();
x = fdiv2(x, y);
} while (moreargs());
xlpop();
}
else
x = div2(FIXONE, x);
return x;
}
#endif /* XLISP_STAT */
LVAL xmin(V)
{
LVAL x, y;
x = xlganumber();
while (moreargs()) {
y = xlgetarg();
x = min2(x, y);
}
return x;
}
LVAL xmax(V)
{
LVAL x, y;
x = xlganumber();
while (moreargs()) {
y = xlgetarg();
x = max2(x, y);
}
return x;
}
#define DEFCMPFUN(name,fun) \
LVAL name(V) { \
LVAL x, y; \
x = xlganumber(); \
while (moreargs()) { \
y = xlgetarg(); \
if (fun(x, y)) x = y; \
else return NIL; \
} \
return s_true; \
}
DEFCMPFUN(xlss, lss2)
DEFCMPFUN(xleq, leq2)
DEFCMPFUN(xequ, equ2)
DEFCMPFUN(xgeq, geq2)
DEFCMPFUN(xgtr, gtr2)
LVAL xneq(V)
{
LVAL x;
int i;
x = xlganumber();
while (moreargs()) {
for (i = 0; i < xlargc; i++)
if (equ2(x, xlargv[i]))
return NIL;
x = xlgetarg();
}
return s_true;
}
#define DEFUNARYFUN(name,fun) \
LVAL name(V) { \
LVAL x = xlgetarg(); \
xllastarg(); \
return fun(x); \
}
DEFUNARYFUN(xadd1, add1)
DEFUNARYFUN(xsub1, sub1)
LVAL xrem(V)
{
LVAL x, y, r;
x = xlgetarg();
y = xlgetarg();
xllastarg();
idiv2(x, y, NULL, &r, 'I');
return r;
}
LVAL xmod(V)
{
LVAL x, y, r;
x = xlgetarg();
y = xlgetarg();
xllastarg();
idiv2(x, y, NULL, &r, '_');
return r;
}
#ifdef MULVALS
#define DEFIDIVFUN(name,w) \
LVAL name (V) { \
LVAL x, y, n, r; \
x = xlgetarg(); \
y = moreargs() ? xlgetarg() : FIXONE; \
xllastarg(); \
idiv2(x, y, &n, &r, w); \
xlresults[0] = n; \
xlresults[1] = r; \
xlnumresults = 2; \
return n; \
}
#define DEFFDIVFUN(name,w) \
LVAL name(V) { \
FLOTYPE xd, yd, nd, rd; \
xd = makefloat(xlgetarg()); \
yd = moreargs() ? makefloat(xlgetarg()) : 1.0; \
xllastarg(); \
flonum_idiv2(xd, yd, &nd, &rd, w); \
xlnumresults = 0; \
xlresults[xlnumresults++] = cvflonum(nd); \
xlresults[xlnumresults++] = cvflonum(rd); \
return xlresults[0]; \
}
#else
#define DEFIDIVFUN(name,w) \
LVAL name(V) { \
LVAL x, y, n, r; \
x = xlgetarg(); \
y = moreargs() ? xlgetarg() : FIXONE; \
xllastarg(); \
idiv2(x, y, &n, &r, w); \
return n; \
}
#define DEFFDIVFUN(name,w) \
LVAL name (V) { \
FLOTYPE xd, yd, nd, rd; \
xd = makefloat(xlgetarg()); \
yd = moreargs() ? makefloat(xlgetarg()) : 1.0; \
xllastarg(); \
flonum_idiv2(xd, yd, &nd, &rd, w); \
return cvflonum(nd); \
}
#endif /* MULVALS */
DEFIDIVFUN(xfix, 'I')
DEFIDIVFUN(xceil, '^')
DEFIDIVFUN(xfloor, '_')
DEFIDIVFUN(xround, 'r')
DEFFDIVFUN(xffix, 'I')
DEFFDIVFUN(xfceil, '^')
DEFFDIVFUN(xffloor, '_')
DEFFDIVFUN(xfround, 'r')
LVAL xexpt(V)
{
LVAL base, power;
base = xlganumber();
power = xlganumber();
xllastarg();
if (fixzerop(power)) {
switch (ntype(base)) {
case FIXNUM:
case BIGNUM:
case RATIO:
return FIXONE;
case FLONUM:
return cvflonum((FLOTYPE) 1.0);
case COMPLEX:
return floatp(getreal(base)) ? newdcomplex(1.0, 0.0) : FIXONE;
}
}
else if (zerop(base) && (zerop(power) || minusp(power)))
badzero();
/* complex float base or complex power */
if ((complexp(base) && floatp(getreal(base))) || complexp(power)) {
dcomplex zb, zp, zv;
makecomplex(&zb, base);
makecomplex(&zp, power);
z_expt(&zv, &zb, &zp);
return(cvcomplex(zv));
}
/* real or complex rational base, integer power */
else if (integerp(power)) {
LVAL val, r;
if (fixp(base) && fixp(power)) {
double fb, fp, fv;
fb = (double) getfixnum(base);
fp = (double) getfixnum(power);
fv = floor(pow(fb, fp > 0.0 ? fp : -fp) + 0.5);
if (infixnumrange(fv)) {
if (fp > 0.0)
return cvfixnum((FIXTYPE) fv);
else
return cvratio(1, (FIXTYPE) fv);
}
}
else if (floatp(base) && fixp(power)) {
double fb, fp, fv;
fb = (double) getflonum(base);
fp = (double) getfixnum(power);
/* test range in case FIXNUMS aren't all exactly represented */
if (infixnumrange(fp)) {
fv = pow(fb, fp);
return cvflonum((FLOTYPE) fv);
}
/* else fall through */
}
if (zerop(base))
return FIXZERO;
xlstkcheck(3);
xlprotect(base);
xlprotect(power);
xlsave(val);
val = FIXONE;
if (minusp(power)) {
base = div2(FIXONE, base);
power = uminus(power);
}
while (TRUE) {
idiv2(power, FIXTWO, &power, &r, 'I');
if (! zerop(r)) {
val = mul2(val, base);
if (zerop(power))
break;
}
base = mul2(base, base);
}
xlpopn(3);
return val;
}
/* real power and float base or rational power and real base */
else {
double fb, fp;
fb = makefloat(base);
fp = makefloat(power);
if (fb == 0.0 && fp > 0.0)
return(cvflonum((FLOTYPE) 0.0));
else if (fb < 0.0) {
dcomplex zb, zp, zv;
makecomplex(&zb, base);
makecomplex(&zp, power);
z_expt(&zv, &zb, &zp);
return(cvcomplex(zv));
}
else {
double fval;
if (fp < 0.0) {
fb = 1.0 / fb;
fp = -fp;
}
fval = pow(fb, fp);
return(cvflonum((FLOTYPE) fval));
}
}
}
#ifdef XLISP_STAT
LVAL xfexpt(V)
{
/* this is a hack but it should do */
if (xlargc == 2 && integerp(xlargv[1]) && minusp(xlargv[1])) {
LVAL b = xlargv[0];
if (rationalp(b))
xlargv[0] = cvflonum((FLOTYPE) makefloat(b));
else if (complexp(b) && rationalp(getreal(b)))
xlargv[0] = newdcomplex(makefloat(getreal(b)), makefloat(getimag(b)));
}
return xexpt();
}
#endif /* XLISP_STAT */
/* arc tangent -- Tom Almy */
LVAL xatan(V)
{
LVAL lnum, ldenom;
dcomplex cnum, cdenom, cval;
lnum = xlgetarg();
if (moreargs()) {
ldenom = xlgetarg();
xllastarg();
switch (commontype(lnum, ldenom)) {
case IN:
case FL:
case RT:
return cvflonum((FLOTYPE)atan2(makefloat(lnum), makefloat(ldenom)));
default: /* complex */
makecomplex(&cnum, lnum);
makecomplex(&cdenom, ldenom);
z_atan2(&cval, &cnum, &cdenom);
return (cvcomplex(cval));
}
}
else {
switch (ntype(lnum)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM:
return cvflonum((FLOTYPE)atan(makefloat(lnum)));
default: /* complex */
makecomplex(&cnum, lnum);
z_atan(&cnum, &cnum);
return (cvcomplex(cnum));
}
}
}
/* two argument logarithm */
LOCAL double logarithm P3C(FLOTYPE, x, FLOTYPE, base, int, base_supplied)
{
double lbase;
if (x <= 0.0) xlfail("logarithm of a nonpositive number");
if (base_supplied) {
if (base <= 0.0) xlfail("logarithm to a nonpositive base");
else {
lbase = log(base);
if (lbase == 0.0) xlfail("logarith to a unit base");
else return((log(x)/lbase));
}
}
return (log(x));
}
LVAL xlog(V)
{
LVAL arg, base;
double fx, fb;
dcomplex zx, zb, zv;
arg = xlgetarg();
if (moreargs()) {
base = xlgetarg();
xllastarg();
if (realp(arg) && realp(base)) {
fx = makefloat(arg);
fb = makefloat(base);
if (fx <= 0.0 || fb <= 0.0) {
makecomplex(&zx, arg);
makecomplex(&zb, base);
z_log(&zx, &zx);
z_log(&zb, &zb);
z_div(&zv, &zx, &zb);
return(cvcomplex(zv));
}
else return(cvflonum((FLOTYPE) logarithm(fx, fb, TRUE)));
}
else if ((realp(arg) && complexp(base))
|| (complexp(arg) && realp(base))
|| (complexp(arg) && complexp(base))) {
makecomplex(&zx, arg);
makecomplex(&zb, base);
z_log(&zx, &zx);
z_log(&zb, &zb);
z_div(&zv, &zx, &zb);
return(cvcomplex(zv));
}
else xlfail("bad argument type(s)");
}
else {
xllastarg();
if (realp(arg)) {
fx = makefloat(arg);
if (fx <= 0.0) {
makecomplex(&zx, arg);
z_log(&zx, &zx);
return(cvcomplex(zx));
}
else return(cvflonum((FLOTYPE) logarithm(fx, 0.0, FALSE)));
}
else if (complexp(arg)) {
makecomplex(&zx, arg);
z_log(&zx, &zx);
return(cvcomplex(zx));
}
else xlfail("bad argument type(s)");
}
return NIL; /* avoid compiler warnings */
}
/* xgcd - greatest common divisor */
LVAL xgcd(V)
{
LVAL m, n;
if (!moreargs()) /* check for identity case */
return FIXZERO;
xlstkcheck(2);
xlsave(m);
xlsave(n);
n = xlgainteger();
if (minusp(n))
n = uminus(n); /* absolute value */
checkzero(n);
while (moreargs()) {
m = xlgainteger();
if (minusp(m))
m = uminus(m); /* absolute value */
checkzero(m);
n = gcd2(n,m);
}
xlpopn(2);
return n;
}
LVAL xlcm(V)
{
LVAL n, m, t, g;
xlstkcheck(4);
xlsave(n);
xlsave(m);
xlsave(t);
xlsave(g);
n = xlgainteger();
while (moreargs()) {
m = xlgainteger();
if (zerop(n) || zerop(m))
return FIXZERO;
t = mul2(n, m);
g = gcd2(n, m);
n = div2(t, g);
}
if (minusp(n))
n = uminus(n);
xlpopn(4);
return n;
}
LVAL xabs(V)
{
LVAL x;
x = xlgetarg();
xllastarg();
switch (ntype(x)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM:
return minusp(x) ? uminus(x) : x;
case COMPLEX:
{
dcomplex zx;
makecomplex(&zx, x);
return cvflonum((FLOTYPE) z_abs(&zx));
}
default:
return xlbadtype(x);
}
}
LVAL xsin(V)
{
LVAL x;
x = xlgetarg();
xllastarg();
switch (ntype(x)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM:
return cvflonum((FLOTYPE) sin(makefloat(x)));
case COMPLEX:
{
dcomplex zx;
makecomplex(&zx, x);
z_sin(&zx, &zx);
return cvcomplex(zx);
}
default:
return xlbadtype(x);
}
}
LVAL xcos(V)
{
LVAL x;
x = xlgetarg();
xllastarg();
switch (ntype(x)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM:
return cvflonum((FLOTYPE) cos(makefloat(x)));
case COMPLEX:
{
dcomplex zx;
makecomplex(&zx, x);
z_cos(&zx, &zx);
return cvcomplex(zx);
}
default:
return xlbadtype(x);
}
}
LVAL xtan(V)
{
LVAL x;
x = xlgetarg();
xllastarg();
switch (ntype(x)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM:
return cvflonum((FLOTYPE) tan(makefloat(x)));
case COMPLEX:
{
dcomplex zx;
makecomplex(&zx, x);
z_tan(&zx, &zx);
return cvcomplex(zx);
}
default:
return xlbadtype(x);
}
}
LVAL xexp(V)
{
LVAL x;
x = xlgetarg();
xllastarg();
switch (ntype(x)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM:
return cvflonum((FLOTYPE) exp(makefloat(x)));
case COMPLEX:
{
dcomplex zx;
makecomplex(&zx, x);
z_exp(&zx, &zx);
return cvcomplex(zx);
}
default:
return xlbadtype(x);
}
}
LVAL xsqrt(V)
{
LVAL x;
x = xlgetarg();
xllastarg();
switch (ntype(x)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM:
{
double dx = makefloat(x);
if (dx >= 0.0)
return cvflonum((FLOTYPE) sqrt(dx));
else
return newdcomplex(0.0, sqrt(-dx));
}
case COMPLEX:
{
dcomplex zx;
makecomplex(&zx, x);
z_sqrt(&zx, &zx);
return cvcomplex(zx);
}
default:
return xlbadtype(x);
}
}
LVAL xfloat(V)
{
LVAL x;
x = xlgetarg();
if (moreargs()) xlgaflonum();
xllastarg();
return floatp(x) ? x : cvflonum((FLOTYPE) makefloat(x));
}
LVAL xasin(V)
{
LVAL x;
x = xlgetarg();
xllastarg();
switch (ntype(x)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM:
{
double dx = makefloat(x);
if (dx <= 1.0 && dx >= -1.0)
return cvflonum((FLOTYPE) asin(dx));
/* else drop through */
}
case COMPLEX:
{
dcomplex zx;
makecomplex(&zx, x);
z_asin(&zx, &zx);
return cvcomplex(zx);
}
default:
return xlbadtype(x);
}
}
LVAL xacos(V)
{
LVAL x;
x = xlgetarg();
xllastarg();
switch (ntype(x)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM:
{
double dx = makefloat(x);
if (dx <= 1.0 && dx >= -1.0)
return cvflonum((FLOTYPE) acos(dx));
/* else drop through */
}
case COMPLEX:
{
dcomplex zx;
makecomplex(&zx, x);
z_acos(&zx, &zx);
return cvcomplex(zx);
}
default:
return xlbadtype(x);
}
}
LVAL xphase(V)
{
LVAL x;
x = xlgetarg();
xllastarg();
switch (ntype(x)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM:
return cvflonum((FLOTYPE) minusp(x) ? PI : 0.0);
case COMPLEX:
{
dcomplex zx;
makecomplex(&zx, x);
return cvflonum((FLOTYPE) z_phase(&zx));
}
default:
return xlbadtype(x);
}
}
LVAL xrational(V)
{
LVAL x;
x = xlgetarg();
xllastarg();
switch (ntype(x)) {
case FIXNUM:
case BIGNUM:
case RATIO:
return x;
case FLONUM:
{
double dummy, original, dx;
int e, sign;
dx = original = getflonum(x);
if (dx < 0.0) {
sign = TRUE;
dx = -dx;
}
else
sign = FALSE;
dx = frexp(dx, &e);
while (MODF(dx, &dummy) != 0.0 && dx <= DBL_MAX / 2.0) {
e--;
dx *= 2.0;
}
if (e >= 0)
return cvtflobignum(original);
else {
LVAL num, denom, val;
xlstkcheck(2);
xlsave(num);
xlsave(denom);
num = cvtflobignum(sign ? -dummy : dummy);
denom = cvtflobignum(ldexp(1.0, -e));
val = cvbratio(num, denom);
xlpopn(2);
return val;
}
}
default:
return xlbadtype(x);
}
}
/* unary predicates */
#define DEFPREDFUN(name,fun) \
LVAL name(V) { \
LVAL x; \
x = xlganumber(); \
xllastarg(); \
return fun(x) ? s_true : NIL; \
}
DEFPREDFUN(xminusp, minusp)
DEFPREDFUN(xzerop, zerop)
DEFPREDFUN(xplusp, plusp)
DEFPREDFUN(xevenp, evenp)
DEFPREDFUN(xoddp, oddp)
/***********************************************************************/
/** **/
/** Complex Number Functions **/
/** **/
/***********************************************************************/
LVAL xcomplex(V) /* TAA rewrite so (complex 12.0) => #c(12.0 0.0) as required
by CL. */
{
LVAL real, imag;
real = xlgetarg();
if (moreargs()) {
imag = xlgetarg();
xllastarg();
return newcomplex(real, imag);
}
if (rationalp(real)) return real;
return newdcomplex(makefloat(real), 0.0);
}
LVAL xconjugate(V)
{
LVAL arg, nimag, val;
arg = xlgetarg();
xllastarg();
switch (ntype(arg)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM:
return arg;
case COMPLEX:
xlsave1(nimag);
nimag = uminus(getimag(arg));
val = newcomplex(getreal(arg), nimag);
xlpop();
return val;
default:
return xlbadtype(arg);
}
}
LVAL xrealpart(V)
{
LVAL arg = xlgetarg();
xllastarg();
switch (ntype(arg)) {
case FIXNUM:
case BIGNUM:
case RATIO:
case FLONUM:
return arg;
case COMPLEX:
return getreal(arg);
default:
xlbadtype(arg);
return(NIL);
}
}
LVAL ximagpart(V)
{
LVAL arg = xlgetarg();
xllastarg();
switch (ntype(arg)) {
case FIXNUM:
case BIGNUM:
case RATIO:
return FIXZERO;
case FLONUM:
return cvflonum((FLOTYPE) 0.0);
case COMPLEX:
return(getimag(arg));
default:
xlbadtype(arg);
return(NIL);
}
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1