/* xlmath2 - 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. */
/* 7/92 TAA -- modified so that ratios that "overflowed" became flonums.
also added code so that integers that overflow become flonums when
NOOVFIXNUM is defined. Complex integers fixed as well */
#include "xlisp.h"
/* The stuff enclosed in NOOVFIXNUM defines an attempt to provide
* extended precision by having integer arithmetic overflow into
* floating point when necessary. It is not true extended precision,
* since subsequent results that fit in the fixed point range are not
* demoted. The system assumes two's complement integer arithmetic.
* MAXFIX should contain the largest positive FIXTYPE; MINFIX, the
* most negative FIXTYPE, is computed as ((- MAXFIX) - 1). Multiplication
* is done in floating point and compared to the FIXTYPE range. This
* seems faster than othr methods with floating point hardware.
* Division only requires a check for -MINFIX / -1 (I think).
*/
#ifndef BIGNUMS
#ifdef NOOVFIXNUM
#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) ((FIXTYPE) MINFIX <= (x) && (x) <= (FIXTYPE) MAXFIX)
#endif
#ifndef MINFIX
#define MINFIX ((- MAXFIX) - 1)
#endif
/* These definititions used instead of those in XLMATH.C */
/* Somewhat cleaned by by Tom Almy */
#define IN 0
#define FL 1
#define CI 2
#define CF 3
#ifdef BIGNUMS
#define RT 4
#endif
typedef struct {
int mode;
FIXTYPE val, crval, cival;
#ifdef BIGNUMS
FIXTYPE num, denom;
#endif
FLOTYPE fval, cfrval, cfival;
} Number;
/* Function prototypes */
LOCAL VOID checkizero P1H(FIXTYPE);
LOCAL VOID checkfzero P1H(FLOTYPE);
LOCAL VOID badiop(V);
LOCAL VOID badfop(V);
#ifdef BIGNUMS
LOCAL VOID badrop(V);
LOCAL VOID add_ratios P3H(FIXTYPE, FIXTYPE, Number *);
LOCAL VOID mult_ratios P3H(FIXTYPE, FIXTYPE, Number *);
#endif
LOCAL VOID badcop(V);
LOCAL LVAL readnumber P1H(Number *);
LOCAL VOID setmode P2H(Number *, int);
LOCAL VOID matchmodes P2H(Number *, Number *);
LOCAL LVAL lispnumber P1H(Number *);
LOCAL LVAL binary P1H(int);
LOCAL LVAL logbinary P1H(int);
LOCAL FIXTYPE xget_gcd P2H(FIXTYPE, FIXTYPE);
LOCAL LVAL unary P1H(int);
LOCAL LVAL unary2 P1H(int);
LOCAL LVAL predicate P1H(int);
LOCAL LVAL compare P1H(int);
LOCAL LVAL ccompare P1H(int);
LOCAL double logarithm P3H(FLOTYPE, FLOTYPE, int);
#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 and messages */
/* checkizero - check for integer division by zero */
LOCAL VOID checkizero P1C(FIXTYPE, iarg)
{
if (iarg == 0)
xlfail("illegal zero argument");
}
/* checkfzero - check for floating point division by zero or log of zero */
LOCAL VOID checkfzero P1C(FLOTYPE, farg)
{
if (farg == 0.0)
xlfail("illegal zero argument");
}
/* badiop - bad integer operation */
LOCAL VOID badiop(V)
{
xlfail("bad integer operation");
}
/* badfop - bad floating point operation */
LOCAL VOID badfop(V)
{
xlfail("bad floating point operation");
}
#ifdef BIGNUMS
/* badrop - bad ratio operation */
LOCAL VOID badrop(V)
{
xlfail("bad ratio operation");
}
#endif
/* badcop - bad complex number operation */
LOCAL VOID badcop(V)
{
xlfail("bad complex number operation");
}
/* TAA MOD badarg() removed, using preexisting xlbadtype which
gives same message */
/* complex - Complex number functions */
/* TAA Mod--do inline, old test for complexp inline as appropriate */
/* badcomplex() eliminated since it didn't make sense (other numereric
types were ok), so returned error is now from xlbadtype */
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);
}
#define cvcomplex(c) newdcomplex((c).r,(c).i)
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;
}
/* Helper functions */
LOCAL LVAL readnumber P1C(Number *, number)
{
LVAL arg = xlgetarg(), real, imag;
switch(ntype(arg)) {
case FIXNUM:
number->mode = IN;
number->val = getfixnum(arg);
break;
case FLONUM:
number->mode = FL;
number->fval = getflonum(arg);
break;
#ifdef BIGNUMS
case RATIO:
number->mode = RT;
number->num = getfixnum(getnumer(arg));
number->denom = getfixnum(getdenom(arg));
break;
#endif
case COMPLEX:
real = getreal(arg);
imag = getimag(arg);
if (fixp(real)) {
number->mode = CI;
number->crval = getfixnum(real);
number->cival = getfixnum(imag);
}
else {
number->mode = CF;
number->cfrval = makefloat(real);
number->cfival = makefloat(imag);
}
break;
default:
xlerror("not a number", arg);
}
return(arg);
}
LOCAL VOID setmode P2C(Number *, x, int, mode)
{
switch (mode) {
#ifdef BIGNUMS
case RT:
if (x->mode != IN) return;
x->mode = mode;
x->num = x->val;
x->denom = 1;
break;
#endif
case FL:
if (x->mode == IN) {
x->mode = mode;
x->fval = x->val;
}
#ifdef BIGNUMS
else if (x->mode == RT) {
x->mode = mode;
x->fval = x->num / (FLOTYPE) x->denom;
}
#endif
else return;
break;
case CI:
if (x->mode != IN) return;
x->mode = mode;
x->crval = x->val;
x->cival = 0;
break;
case CF:
switch (x->mode) {
case IN:
x->mode = mode;
x->cfrval = x->val;
x->cfival = 0.0;
break;
#ifdef BIGNUMS
case RT:
x->mode = mode;
x->cfrval = x->num / (FLOTYPE) x->denom;
x->cfival = 0.0;
break;
#endif
case FL:
x->mode = mode;
x->cfrval = x->fval;
x->cfival = 0.0;
break;
case CI:
x->mode = mode;
x->cfrval = x->crval;
x->cfival = x->cival;
break;
}
break;
}
}
LOCAL VOID matchmodes P2C(Number *, x, Number *, y)
{
int mode = x->mode;
switch (mode) {
case IN: mode = y->mode; break;
case FL: if (y->mode == CI || y->mode == CF) mode = CF; break;
#ifdef BIGNUMS
case CI: if (y->mode == FL || y->mode == CF || y->mode == RT)
mode = CF;
break;
#else
case CI: if (y->mode == FL || y->mode == CF) mode = CF; break;
#endif
case CF: break;
#ifdef BIGNUMS
case RT: if (y->mode == CI) mode = CF;
else if (y->mode != IN) mode = y->mode;
break;
#endif
}
if (x->mode != mode) setmode(x, mode);
if (y->mode != mode) setmode(y, mode);
}
LOCAL LVAL lispnumber P1C(Number *, x)
{
switch (x->mode) {
case IN: return(cvfixnum(x->val));
case FL: return(cvflonum(x->fval));
case CI: return(newicomplex(x->crval, x->cival));
case CF: return(newdcomplex(x->cfrval, x->cfival));
#ifdef BIGNUMS
case RT: return(cvratio (x->num, x->denom));
#endif
}
return NIL; /* avoid warning messages */
}
/* TAA Mod to return FIXTYPE */
LOCAL FIXTYPE xget_gcd P2C(FIXTYPE, n, FIXTYPE, m) /* euclid's algorith */
{
FIXTYPE r;
for (;;) {
r = m % n;
if (r == (FIXTYPE) 0)
break;
m = n;
n = r;
}
return (n);
}
#ifdef BIGNUMS
LOCAL VOID add_ratios P3C(FIXTYPE, n1, FIXTYPE, d1, Number *, res)
{
double lcm, nt;
FIXTYPE n2=res->num, d2=res->denom;
lcm = d1 * (d2 /(double)(xget_gcd(d1, d2)));
nt = n1 * (lcm / d1) + n2 * (lcm / d2);
if (fabs(nt) > (double)MAXFIX || fabs(lcm) > (double)MAXFIX) {
res->fval = (n1/(double)d1) + (n2/(double)d2);
res->mode = FL;
}
else {
res->num = (FIXTYPE) nt;
res->denom = (FIXTYPE) lcm;
/* mode is already RT */
}
return;
}
LOCAL VOID mult_ratios P3C(FIXTYPE, n1, FIXTYPE, d1, Number *, res)
{
FIXTYPE gcd;
FIXTYPE n2=res->num, d2=res->denom;
double nt, dt;
/* TAA rewrite from 2.1d -- reduce first, can help fit */
if (n1 == 0 || n2 == 0) { /* handle multiply by zero */
res->val = 0;
res->mode = IN;
return;
}
if ((gcd=xget_gcd(n1, d2)) > 1) {
n1 /= gcd;
d2 /= gcd;
}
if ((gcd=xget_gcd(d1, n2)) > 1) {
d1 /= gcd;
n2 /= gcd;
}
nt = n1 * (double)n2;
dt = d1 * (double)d2;
if (fabs(nt) > (double)MAXFIX || fabs(dt) > (double)MAXFIX) {
res->fval = nt/dt;
res->mode = FL;
}
else {
res->num = (FIXTYPE) nt;
res->denom = (FIXTYPE) dt;
/* mode is already RT */
}
return;
}
#endif
LOCAL LVAL binary P1C(int, which)
{
LVAL larg;
Number val, arg;
FIXTYPE rtemp, itemp;
FLOTYPE frtemp, fitemp;
if (xlargc == 1 && (which == '-' || which == '/'
#ifdef XLISP_STAT
|| which == 'd'
#endif /* XLISP_STAT */
)) {
val.mode = IN;
switch (which) {
case '-': val.val = 0; break;
#ifdef XLISP_STAT
case 'd':
#endif /* XLISP_STAT */
case '/': val.val = 1; break;
}
}
else larg = readnumber(&val);
while (moreargs()) {
larg = readnumber(&arg);
matchmodes(&val, &arg);
switch (which) {
case '+':
switch (val.mode) {
case IN:
#ifdef NOOVFIXNUM
rtemp = val.val + arg.val;
if (goodisum(val.val, arg.val, rtemp))
val.val = rtemp;
else {
val.fval = val.val + (double) arg.val;
val.mode = FL;
}
#else
val.val += arg.val;
#endif /* NOOVFIXNUM */
break;
case FL: val.fval += arg.fval; break;
case CI:
#ifdef NOOVFIXNUM
rtemp = val.crval + arg.crval;
itemp = val.cival + arg.cival;
if (goodisum(val.crval, arg.crval, rtemp) &&
goodisum(val.cival, arg.cival, itemp)) {
val.crval = rtemp;
val.cival = itemp;
}
else {
val.cfrval = val.crval+(double)arg.crval;
val.cfival = val.cival+(double)arg.cival;
val.mode = CF;
}
#else
val.crval += arg.crval;
val.cival += arg.cival;
#endif /* NOOVFIXNUM */
break;
case CF: val.cfrval += arg.cfrval; val.cfival += arg.cfival; break;
#ifdef BIGNUMS
case RT:
add_ratios (arg.num, arg.denom, &val);
break;
#endif
}
break;
case '-':
switch (val.mode) {
case IN:
#ifdef NOOVFIXNUM
rtemp = val.val - arg.val;
if (goodidiff(val.val, arg.val, rtemp))
val.val = rtemp;
else {
val.fval = val.val - (double) arg.val;
val.mode = FL;
}
#else
val.val -= arg.val;
#endif /* NOOVFIXNUM */
break;
case FL: val.fval -= arg.fval; break;
case CI:
#ifdef NOOVFIXNUM
rtemp = val.crval - arg.crval;
itemp = val.cival - arg.cival;
if (goodidiff(val.crval, arg.crval, rtemp) &&
goodidiff(val.cival, arg.cival, itemp)) {
val.crval = rtemp;
val.cival = itemp;
}
else {
val.cfrval = val.crval-(double)arg.crval;
val.cfival = val.cival-(double)arg.cival;
val.mode = CF;
}
#else
val.crval -= arg.crval;
val.cival -= arg.cival;
#endif /* NOOVFIXNUM */
break;
case CF: val.cfrval -= arg.cfrval; val.cfival -= arg.cfival; break;
#ifdef BIGNUMS
case RT:
add_ratios (-arg.num, arg.denom, &val);
break;
#endif
}
break;
case '*':
switch (val.mode) {
case IN:
#ifdef NOOVFIXNUM
frtemp = val.val * (double)arg.val; /* do math in floating point! */
if (infixnumrange(frtemp))
val.val = (FIXTYPE)frtemp;
else {
val.fval = frtemp;
val.mode = FL;
}
#else
val.val *= arg.val;
#endif /* NOOVFIXNUM */
break;
case FL: val.fval *= arg.fval; break;
case CI:
#ifdef NOOVFIXNUM
frtemp = val.crval * (double)arg.crval - val.cival * (double)arg.cival;
fitemp = val.cival * (double)arg.crval + val.crval * (double)arg.cival;
if (infixnumrange(frtemp) && infixnumrange(fitemp)) {
val.crval = (FIXTYPE)frtemp;
val.cival = (FIXTYPE)fitemp;
}
else {
val.cfrval = frtemp;
val.cfival = fitemp;
val.mode = CF;
}
#else
rtemp = val.crval * arg.crval - val.cival * arg.cival;
itemp = val.cival * arg.crval + val.crval * arg.cival;
val.crval = rtemp; val.cival = itemp;
#endif /* NOOVFIXNUM */
break;
case CF:
frtemp = val.cfrval * arg.cfrval - val.cfival * arg.cfival;
fitemp = val.cfival * arg.cfrval + val.cfrval * arg.cfival;
val.cfrval = frtemp; val.cfival = fitemp;
break;
#ifdef BIGNUMS
case RT:
mult_ratios (arg.num, arg.denom, &val);
break;
#endif
}
break;
case '/':
switch (val.mode) {
case IN:
checkizero(arg.val);
#ifdef BIGNUMS
setmode(&val, RT);
val.denom = arg.val;
break;
#else
if (
#ifdef NOOVFIXNUM
(arg.val != -1 || val.val > MINFIX) &&
#endif
val.val % arg.val == 0) {
val.val /= arg.val;
break;
}
else {
setmode(&val, FL);
setmode(&arg, FL);
}
/* drop through */
#endif
case FL:
checkfzero(arg.fval);
val.fval /= arg.fval;
break;
case CI:
setmode(&val, CF);
setmode(&arg, CF);
/* drop through */
case CF:
{
double ratio,temp;
if (fabs(arg.cfrval) > fabs(arg.cfival)) {
ratio = arg.cfival / arg.cfrval;
temp = arg.cfrval + ratio*arg.cfival;
frtemp = (val.cfrval + val.cfival*ratio)/temp;
fitemp = (val.cfival - val.cfrval*ratio)/temp;
}
else {
checkfzero(arg.cfival);
ratio = arg.cfrval / arg.cfival;
temp = arg.cfival + ratio*arg.cfrval;
frtemp = (val.cfrval*ratio + val.cfival)/temp;
fitemp = (val.cfival*ratio - val.cfrval)/temp;
}
val.cfrval = frtemp; val.cfival = fitemp;
break;
}
#ifdef BIGNUMS
case RT:
checkizero (arg.num);
mult_ratios (arg.denom, arg.num, &val);
break;
#endif
}
break;
#ifdef XLISP_STAT
case 'd':
switch (val.mode) {
case IN:
checkizero(arg.val);
if (
#ifdef NOOVFIXNUM
(arg.val != -1 || val.val > MINFIX) &&
#endif
val.val % arg.val == 0) {
val.val /= arg.val;
break;
}
else {
setmode(&val, FL);
setmode(&arg, FL);
}
/* drop through */
case FL:
checkfzero(arg.fval);
val.fval /= arg.fval;
break;
case CI:
setmode(&val, CF);
setmode(&arg, CF);
/* drop through */
case CF:
{
double ratio,temp;
if (fabs(arg.cfrval) > fabs(arg.cfival)) {
ratio = arg.cfival / arg.cfrval;
temp = arg.cfrval + ratio*arg.cfival;
frtemp = (val.cfrval + val.cfival*ratio)/temp;
fitemp = (val.cfival - val.cfrval*ratio)/temp;
}
else {
checkfzero(arg.cfival);
ratio = arg.cfrval / arg.cfival;
temp = arg.cfival + ratio*arg.cfrval;
frtemp = (val.cfrval*ratio + val.cfival)/temp;
fitemp = (val.cfival*ratio - val.cfrval)/temp;
}
val.cfrval = frtemp; val.cfival = fitemp;
break;
}
#ifdef BIGNUMS
case RT:
checkizero (arg.num);
mult_ratios (arg.denom, arg.num, &val);
break;
#endif
}
break;
#endif /* XLISP_STAT */
case 'M':
switch (val.mode) {
case IN: val.val = (val.val > arg.val) ? val.val : arg.val; break;
case FL: val.fval = (val.fval > arg.fval) ? val.fval : arg.fval; break;
#ifdef BIGNUMS
case RT:
if ((val.num * (double)arg.denom) < (arg.num * (double)val.denom)) {
val.num = arg.num;
val.denom = arg.denom;
}
break;
#endif
default: xlbadtype(larg);
}
break;
case 'm':
switch (val.mode) {
case IN: val.val = (val.val < arg.val) ? val.val : arg.val; break;
case FL: val.fval = (val.fval < arg.fval) ? val.fval : arg.fval; break;
#ifdef BIGNUMS
case RT:
if ((val.num * (double)arg.denom) > (arg.num * (double)val.denom)) {
val.num = arg.num;
val.denom = arg.denom;
}
break;
#endif
default: xlbadtype(larg);
}
break;
}
}
return(lispnumber(&val));
}
/* This has been completely rewritten by Tom Almy to handle floating point
arguments */
LVAL xrem(V)
{
Number numer, div;
double ftemp;
readnumber(&numer);
readnumber(&div);
xllastarg();
matchmodes(&numer, &div);
switch (numer.mode) {
case IN:
checkizero(div.val);
return (cvfixnum((FIXTYPE) numer.val % div.val));
case FL:
checkfzero(div.fval);
MODF(numer.fval / div.fval, &ftemp);
return (cvflonum((FLOTYPE)(numer.fval - ftemp*div.fval)));
#ifdef BIGNUMS
case RT:
checkizero(div.num);
mult_ratios(div.denom, div.num, &numer);
if (numer.mode == RT) {
numer.num %= numer.denom; /* get remainder */
mult_ratios(div.num, div.denom, &numer);
}
else if (numer.mode == FL) {
numer.fval = numer.num/numer.denom;
MODF(numer.fval / div.fval, &ftemp);
numer.fval -= ftemp*div.fval;
}
return (lispnumber(&numer));
#endif
}
badcop();
return NIL; /* fool compiler into not giving warning */
}
LVAL xash(V) /* arithmetic shift left */
{
FIXTYPE arg, val;
arg = getfixnum(xlgafixnum());
val = getfixnum(xlgafixnum());
xllastarg();
return cvfixnum(val < 0 ? arg >> -val : arg << val);
}
LOCAL LVAL logbinary P1C(int, which)
{
FIXTYPE val, arg; /* TAA Mod -- was int */
switch (which) {
case '&': val = -1; break;
case '|': val = 0; break;
case '^': val = 0; break;
default: val = 0; /* to keep compiler happy */
}
while (moreargs()) {
arg = getfixnum(xlgafixnum());
switch (which) {
case '&': val &= arg; break;
case '|': val |= arg; break;
case '^': val ^= arg; break;
}
}
return(cvfixnum((FIXTYPE) val));
}
/* binary functions */
/* TAA fix allowing (+) */
LVAL xadd(V) { return (moreargs()?binary('+'):cvfixnum((FIXTYPE)0)); }
LVAL xsub() { return (binary('-')); } /* - */
/* TAA fix allowing (*) */
LVAL xmul(V) { return (moreargs()?binary('*'):cvfixnum((FIXTYPE)1)); }
LVAL xdiv(V) { return (binary('/')); } /* / */
#ifdef XLISP_STAT
LVAL xfdiv(V) { return (binary('d')); } /* / */
#endif /* XLISP_STAT */
LVAL xmin(V) { return (binary('m')); } /* min */
LVAL xmax(V) { return (binary('M')); } /* max */
LVAL xlogand(V) { return (logbinary('&')); } /* logand */
LVAL xlogior(V) { return (logbinary('|')); } /* logior */
LVAL xlogxor(V) { return (logbinary('^')); } /* logxor */
/* New by Tom Almy */
LVAL xmod(V)
{
Number numer, div;
readnumber(&numer);
readnumber(&div);
xllastarg();
matchmodes(&numer, &div);
switch (numer.mode) {
case IN:
checkizero(div.val);
return(cvfixnum(numer.val -
div.val*(FIXTYPE)floor(numer.val/(double)div.val)));
case FL:
checkfzero(div.fval);
return(cvflonum((FLOTYPE)(numer.fval -
div.fval*floor((double)(numer.fval/div.fval)))));
#ifdef BIGNUMS
case RT:
checkizero(div.num);
mult_ratios(div.denom, div.num, &numer);
if (numer.mode == RT) {
numer.num = numer.num - numer.denom*
(FIXTYPE)floor(numer.num/(double)numer.denom);
mult_ratios(div.num, div.denom, &numer);
}
else if (numer.mode == FL) {
numer.fval = numer.num/numer.denom;
numer.fval -= div.fval*floor((double)(numer.fval/div.fval));
}
return (lispnumber(&numer));
#endif
}
badcop();
return NIL; /* fool compiler into not giving warning */
}
LVAL xexpt(V) /* modified for ratios by TAA */
{
LVAL base, power;
int bsign, psign;
FIXTYPE b, p, val;
FLOTYPE fb, fp, fval;
base = xlgetarg();
power = xlgetarg();
xllastarg();
if (fixp(base) && fixp(power)) {
b = getfixnum(base);
p = getfixnum(power);
if (p == 0) return(cvfixnum((FIXTYPE) 1));
if (b == 0 && p > 0) return(cvfixnum((FIXTYPE) 0));
checkizero(b);
if ((bsign = (b < 0)) != 0) b = -b;
if ((psign = (p < 0)) != 0) p = -p;
fval = floor(pow((double) b, (double) p) + 0.1); /* to get integer right */
if (bsign && (p & 1)) fval = -fval;
if (!psign) {
#ifdef NOOVFIXNUM
if (infixnumrange(fval)) {
val = fval;
return(cvfixnum((FIXTYPE) val));
}
#else
val = fval; /* this method causes errors on intel architectures */
if (val == fval) return(cvfixnum(val));
#endif /* NOOVFIXNUM */
else /* to handle precision for large results */
return(cvflonum((FLOTYPE) fval));
}
else {
#ifdef BIGNUMS
val = fval;
if (val == fval) return (cvratio(1, val));
else return (cvflonum((FLOTYPE) 1.0 / fval));
#else
checkfzero(fval);
return(cvflonum((FLOTYPE) 1.0 / fval));
#endif
}
}
#ifdef BIGNUMS
else if (ratiop(base) && fixp(power)) {
FIXTYPE vald; /* integer of denominator result */
FLOTYPE fvald; /* denominator result */
b = getfixnum(getnumer(base)); /* start with just the numerator */
p = getfixnum(power);
if (p == 0) return(cvfixnum((FIXTYPE) 1));
if ((bsign = (b < 0)) != 0) b = -b;
if ((psign = (p < 0)) != 0) p = -p;
fval = floor(pow((double) b, (double) p) + 0.1); /* to get integer right */
if (bsign && (p & 1)) fval = -fval;
fvald = floor(pow((double) getfixnum(getdenom(base)), (double) p) + 0.1);
val = fval;
vald = fvald;
if (!psign) {
if (val == fval && vald == fvald) /* will fit in ratio */
return (cvratio(val, vald));
else return (cvflonum(fval/fvald));
}
else {
if (val == fval && vald == fvald) /* will fit in ratio */
return (cvratio(vald, val));
else return (cvflonum(fvald/fval));
}
}
#endif
else if (floatp(base) && fixp(power)) {
fb = getflonum(base);
p = getfixnum(power);
if (p == 0) return(cvflonum((FLOTYPE) 1.0)); /* TAA MOD - used to return
fixnum 1, but CL says should be flonum */
if (fb == 0.0 && p > 0) return(cvflonum((FLOTYPE) 0.0));
checkfzero(fb);
if ((bsign = (fb < 0.0)) != 0) fb = -fb;
if ((psign = (p < 0)) != 0) p = -p;
fval = pow((double) fb, (double) p);
if (bsign && (p & 1)) fval = -fval;
if (!psign) return(cvflonum((FLOTYPE) fval));
else {
checkfzero(fval);
return(cvflonum((FLOTYPE) 1.0 / fval));
}
}
#ifdef BIGNUMS
else if (realp(base) && (ratiop(power) || floatp(power)))
#else
else if (realp(base) && floatp(power))
#endif
{
fb = makefloat(base);
#ifdef BIGNUMS
fp = ratiop(power) ?
getfixnum(getnumer(power))/(FLOTYPE)getfixnum(getdenom(power)) : getflonum(power);
#else
fp = getflonum(power);
#endif
if (fp == 0.0) return(cvflonum((FLOTYPE) 1.0));
if (fb == 0.0 && fp > 0.0) return(cvflonum((FLOTYPE) 0.0));
if (fb < 0.0) {
dcomplex zb, zp, zv;
makecomplex(&zb, base);
makecomplex(&zp, power);
z_expt(&zv, &zb, &zp);
return(cvcomplex(zv));
}
if ((psign = (fp < 0.0)) != 0) fp = -fp;
fval = pow((double) fb, (double) fp);
if (!psign) return(cvflonum((FLOTYPE) fval));
else {
checkfzero(fval);
return(cvflonum((FLOTYPE) 1.0 / fval));
}
}
else if (complexp(base) || complexp(power)) {
dcomplex zb, zp, zv;
makecomplex(&zb, base);
makecomplex(&zp, power);
z_expt(&zv, &zb, &zp);
return(cvcomplex(zv));
}
else xlfail("bad argument type(s)");
return NIL; /* avoid compiler warnings */
}
/* arc tangent -- Tom Almy */
LVAL xatan(V)
{
Number numer, denom;
LVAL lnum, ldenom;
dcomplex cnum, cdenom, cval;
lnum = readnumber(&numer);
if (moreargs()) {
ldenom = readnumber(&denom);
xllastarg();
matchmodes(&numer, &denom);
switch (numer.mode) {
case IN:
numer.fval = numer.val; denom.fval = denom.val;
case FL:
return (cvflonum((FLOTYPE)atan2(numer.fval, denom.fval)));
#ifdef BIGNUMS
case RT:
return (cvflonum((FLOTYPE)atan2(numer.num/(FLOTYPE)numer.denom,
denom.num/(FLOTYPE)denom.denom)));
#endif
default: /* complex */
makecomplex(&cnum, lnum);
makecomplex(&cdenom, ldenom);
z_atan2(&cval, &cnum, &cdenom);
return (cvcomplex(cval));
}
}
else {
switch (numer.mode) {
case IN:
numer.fval = numer.val;
case FL:
return (cvflonum((FLOTYPE)atan(numer.fval)));
#ifdef BIGNUMS
case RT:
return (cvflonum((FLOTYPE)atan(numer.num/(FLOTYPE)numer.denom)));
#endif
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)
{
FIXTYPE m,n;
LVAL arg;
if (!moreargs()) /* check for identity case */
return (cvfixnum((FIXTYPE)0));
arg = xlgafixnum();
n = getfixnum(arg);
if (n < (FIXTYPE)0) n = -n; /* absolute value */
while (moreargs()) {
arg = xlgafixnum();
m = getfixnum(arg);
if (m == 0 || n == 0) xlfail("zero argument");
if (m < (FIXTYPE)0) m = -m; /* absolute value */
n = xget_gcd(n,m);
}
return (cvfixnum(n));
}
LVAL xlcm(V) /* added by kcw */
{
LVAL arg;
FIXTYPE n, m, g;
double t; /* TAA mod, 7/72, to reduce chance of overflow */
arg = xlgafixnum();
n = getfixnum(arg);
if (!moreargs()) {
if (n < (FIXTYPE) 0) n = -n;
return (cvfixnum(n));
}
while (moreargs()) {
arg = xlgafixnum();
m = getfixnum(arg);
if ((n == (FIXTYPE) 0) || (m == (FIXTYPE) 0))
return(cvfixnum(0));
t = n * (double)m;
g = xget_gcd(n,m);
t = t/g;
if (fabs(t) > (double)MAXFIX) xlfail("lcm larger than maximum fixnum");
n = (FIXTYPE)t;
}
if (n < (FIXTYPE) 0) n = -n;
return (cvfixnum(n));
}
/* unary - handle unary operations */
LOCAL LVAL unary P1C(int, which)
{
FLOTYPE fval;
FIXTYPE ival;
FIXTYPE itemp;
#ifdef BIGNUMS
FIXTYPE numer, denom;
#endif
dcomplex cval;
LVAL arg, real, imag;
int mode;
/* get the argument */
arg = xlgetarg();
if (which == 'F' && moreargs()) xlgaflonum();
xllastarg();
/* check its type */
if (fixp(arg)) {
ival = getfixnum(arg);
mode = IN;
}
else if (floatp(arg)) {
fval = getflonum(arg);
mode = FL;
}
#ifdef BIGNUMS
else if (ratiop(arg)) {
numer = getfixnum(getnumer(arg));
denom = getfixnum(getdenom(arg));
mode = RT;
}
#endif
else if (complexp(arg)) {
makecomplex(&cval, arg);
real = getreal(arg);
imag = getimag(arg);
if (fixp(getreal(arg))) mode = CI;
else mode = CF;
}
else xlerror("not a number", arg);
switch (which) {
case '~':
if (mode == IN) return(cvfixnum((FIXTYPE) ~ival));
else badiop();
break;
case 'A':
switch (mode) {
#ifdef NOOVFIXNUM
case IN:
itemp = ival < 0 ? -ival : ival;
if (0 <= itemp) return(cvfixnum((FIXTYPE) itemp));
else fval = ival;
/* drop through */
#else
case IN: return(cvfixnum((FIXTYPE) (ival < 0 ? -ival : ival)));
#endif /* NOOVFIXNUM */
case FL: return(cvflonum((FLOTYPE) (fval < 0.0 ? -fval : fval)));
case CI:
case CF: return(cvflonum((FLOTYPE) z_abs(&cval)));
#ifdef BIGNUMS
case RT: return(cvratio(numer<0?-numer:numer,denom));
#endif
}
break;
case '+':
switch (mode) {
case IN:
itemp = ival + 1;
#ifdef NOOVFIXNUM
if (itemp < 0 && ival > 0) return(cvflonum(ival + (FLOTYPE) 1.0));
/* drop through */
#endif /* NOOVFIXNUM */
return(cvfixnum((FIXTYPE) ival + 1));
case FL: return(cvflonum((FLOTYPE) fval + 1.0));
case CI:
ival = getfixnum(real);
itemp = ival + 1;
#ifdef NOOVFIXNUM
if (itemp < 0 && ival > 0)
return(newdcomplex(ival + (FLOTYPE) 1.0, (FLOTYPE)getfixnum(imag)));
/* else drop through */
#endif
return(newicomplex(itemp, getfixnum(imag)));
case CF: return(newdcomplex(getflonum(real) + 1.0, getflonum(imag)));
#ifdef BIGNUMS
case RT:
if ((numer+(double)denom) > (double)MAXFIX)
return (cvflonum((FLOTYPE)(1.0+ numer/(double)denom)));
else
return(cvratio(numer+denom, denom));
#endif
}
break;
case '-':
switch (mode) {
case IN:
itemp = ival - 1;
#ifdef NOOVFIXNUM
if (ival < 0 && itemp > 0)
return(cvflonum(ival - (FLOTYPE) 1.0));
/* drop through */
#endif /* NOOVFIXNUM */
return(cvfixnum((FIXTYPE) itemp));
case FL: return(cvflonum((FLOTYPE) fval - 1.0));
case CI:
ival = getfixnum(real);
itemp = ival - 1;
#ifdef NOOVFIXNUM
if (ival < 0 && itemp > 0)
return(newdcomplex(ival - (FLOTYPE) 1.0, (FLOTYPE)getfixnum(imag)));
/* else drop through */
#endif /* NOOVFIXNUM */
return(newicomplex(itemp, getfixnum(imag)));
case CF: return(newdcomplex(getflonum(real) - 1.0, getflonum(imag)));
#ifdef BIGNUMS
case RT:
if ((numer-(double)denom) < (double)-MAXFIX)
return (cvflonum((FLOTYPE)(-1.0 + numer/(double)denom)));
else
return(cvratio(numer-denom, denom));
#endif
}
break;
case 'S':
switch (mode) {
case IN: return(cvflonum((FLOTYPE) sin((double) ival)));
#ifdef BIGNUMS
case RT: fval = numer / (FLOTYPE) denom;
#endif
case FL: return(cvflonum((FLOTYPE) sin((double) fval)));
case CI:
case CF:
z_sin(&cval, &cval);
return(cvcomplex(cval));
}
case 'C':
switch (mode) {
case IN: return(cvflonum((FLOTYPE) cos((double) ival)));
#ifdef BIGNUMS
case RT: fval = numer / (FLOTYPE) denom;
#endif
case FL: return(cvflonum((FLOTYPE) cos((double) fval)));
case CI:
case CF:
z_cos(&cval, &cval);
return(cvcomplex(cval));
}
case 'T':
switch (mode) {
case IN: return(cvflonum((FLOTYPE) tan((double) ival)));
#ifdef BIGNUMS
case RT: fval = numer / (FLOTYPE) denom;
#endif
case FL: return(cvflonum((FLOTYPE) tan((double) fval)));
case CI:
case CF:
z_tan(&cval, &cval);
return(cvcomplex(cval));
}
case 'E':
switch (mode) {
case IN: return(cvflonum((FLOTYPE) exp((double) ival)));
#ifdef BIGNUMS
case RT: fval = numer / (FLOTYPE) denom;
#endif
case FL: return(cvflonum((FLOTYPE) exp((double) fval)));
case CI:
case CF:
z_exp(&cval, &cval);
return(cvcomplex(cval));
}
break;
case 'R':
switch (mode) {
case IN:
if (ival < 0) {
makecomplex(&cval, arg);
z_sqrt(&cval, &cval);
return(cvcomplex(cval));
}
else return(cvflonum((FLOTYPE) sqrt((double) ival)));
#ifdef BIGNUMS
case RT: fval = numer / (FLOTYPE) denom;
#endif
case FL:
if (fval < 0) {
makecomplex(&cval, arg);
z_sqrt(&cval, &cval);
return(cvcomplex(cval));
}
else return(cvflonum((FLOTYPE) sqrt(fval)));
case CI:
case CF:
z_sqrt(&cval, &cval);
return(cvcomplex(cval));
}
break;
case 'F':
switch (mode) {
case IN: return (cvflonum((FLOTYPE) ival));
#ifdef BIGNUMS
case RT: fval = numer / (FLOTYPE) denom;
#endif
case FL: return (cvflonum((FLOTYPE) fval));
default: badcop();
}
break;
case 's':
switch (mode) {
#ifdef BIGNUMS
case RT: fval = numer / (FLOTYPE) denom; goto l1;
#endif
case IN:
fval = ival;
/* drop through */
case FL:
#ifdef BIGNUMS
l1:
#endif
if (fval > 1.0 || fval < -1.0) {
makecomplex(&cval, arg);
z_asin(&cval, &cval);
return(cvcomplex(cval));
}
else return(cvflonum((FLOTYPE) asin(fval)));
case CI:
case CF:
z_asin(&cval, &cval);
return(cvcomplex(cval));
}
break;
case 'c':
switch (mode) {
#ifdef BIGNUMS
case RT: fval = numer / (FLOTYPE) denom; goto l2;
#endif
case IN:
fval = ival;
/* drop through */
case FL:
#ifdef BIGNUMS
l2:
#endif
if (fval > 1.0 || fval < -1.0) {
makecomplex(&cval, arg);
z_acos(&cval, &cval);
return(cvcomplex(cval));
}
else return(cvflonum((FLOTYPE) acos(fval)));
case CI:
case CF:
z_acos(&cval, &cval);
return(cvcomplex(cval));
}
break;
case 'P':
switch (mode) {
case IN: return(cvflonum((FLOTYPE) (ival >= 0) ? 0.0 : PI));
#ifdef BIGNUMS
case RT: fval = numer / (FLOTYPE) denom;
#endif
case FL: return(cvflonum((FLOTYPE) (fval >= 0.0) ? 0.0 : PI));
case CI:
case CF: return(cvflonum((FLOTYPE) z_phase(&cval)));
}
break;
#ifdef BIGNUMS
case 'X':
switch (mode) {
case IN:
case RT: return arg;
case FL:
{
/* NOTE-- This code is specific to 32 bit signed integer math */
/****** FIX THIS? */
double dummy, original=fval;
int e, sign=FALSE;
if (fval < 0.0) {
sign = TRUE;
fval = -fval;
}
fval = frexp(fval, &e);
for (; MODF(fval, &dummy)!=0.0 && e>-30 && fval<=1073741823.; e--)
fval *= 2.0;
if (e >= 0) {
if (fabs(original) > MAXFIX) /* won't fit -- return float */
return (cvflonum((FLOTYPE)original));
else
return(cvfixnum((FIXTYPE)original));
}
return (cvratio((FIXTYPE)(sign? -dummy : dummy),
(FIXTYPE)ldexp(1.0, -e)));
}
case CI:
case CF: badcop();
}
break;
#endif
default: xlfail("unsupported operation");
}
return NIL; /* avoid compiler warning */
}
/**** need to look at thee for overflow handling */
LOCAL LVAL unary2 P1C(int, which) /* handle truncate, floor, ceiling, and round */
/* 1 or two arguments */
/* By Tom Almy */
{
#ifdef MULVALS /* Yep, it all has to be redone for MULVALS! TAA 10/93 */
Number numer, denom;
double quotient;
LVAL lval;
xlnumresults = 2; /* Always the case */
lval = readnumber(&numer);
if (moreargs()) { /* two argument version */
readnumber(&denom);
xllastarg();
matchmodes(&numer, &denom);
switch (numer.mode) {
case IN:
checkizero(denom.val);
quotient = numer.val / (double)denom.val;
break;
case FL:
checkfzero(denom.fval);
quotient = numer.fval/denom.fval;
break;
#ifdef BIGNUMS
case RT:
checkizero(denom.num);
quotient = ((double)numer.num/numer.denom)/
((double)denom.num/denom.denom);
break;
#endif
default: badcop();
}
}
else { /* single argument version, denominator is one */
switch (numer.mode) {
case IN:
xlresults[1] = cvfixnum((FIXTYPE) 0);
return (xlresults[0]=lval); /* done early! */
case FL:
denom.fval = 1.0;
quotient = numer.fval;
break;
#ifdef BIGNUMS
case RT:
denom.mode = RT;
denom.num = denom.denom = 1;
quotient = (double)numer.num/numer.denom;
break;
#endif
default: badcop();
}
}
switch (which) { /* now do it! */
case 'I': MODF(quotient,"ient); break;
case '_': quotient = floor(quotient); break;
case '^': quotient = ceil(quotient); break;
/* round is flakey when value is too big for a long --
then it (probably) rounds the .5 case "randomly" */
case 'r': quotient = (((((long)quotient & 1)==0)^(quotient < 0))?
ceil(quotient - 0.5) :
floor(quotient + 0.5)); break;
}
/* calculate remainder in same type as arguments */
switch (numer.mode) {
#ifdef BIGNUMS
case RT:
/* we can't do it using ratios if the quotient is too big */
if (fabs(quotient) < (double) MAXFIX) {
mult_ratios((FIXTYPE)quotient, (FIXTYPE)1, &denom);
switch (denom.mode) {
case RT: /* we are still a ratio -- good news! */
add_ratios(-denom.num, denom.denom, &numer);
xlresults[1] = cvratio(numer.num, numer.denom);
break;
case IN: /* it has become an integer */
add_ratios(denom.val, (FIXTYPE) 1, &numer);
xlresults[1] = lispnumber(&numer);
break;
case FL: /* became floating after multiply --
can this really happen??*/
xlresults[1]=cvflonum((FLOTYPE)
((double)numer.num/numer.denom
- denom.fval));
break;
}
}
else {
/* quotient is too big -- do it all as floats */
xlresults[1] = cvflonum((FLOTYPE)
((double)numer.num/numer.denom -
quotient*((double)denom.num/denom.denom)));
}
break;
#endif
case IN:
xlresults[1] = cvfixnum(numer.val - (FIXTYPE)quotient*denom.val);
break;
case FL:
xlresults[1] = cvflonum((FLOTYPE)(numer.fval-quotient*denom.fval));
break;
}
/* quotient is (almost) always an integer, but allow for overflow */
return (xlresults[0] = (numer.mode==IN || fabs(quotient)<(double)MAXFIX)?
cvfixnum((FIXTYPE)quotient):
cvflonum((FLOTYPE)quotient));
#else
Number numer, denom;
LVAL lval;
lval = readnumber(&numer);
if (moreargs()) { /* two argument version */
readnumber(&denom);
xllastarg();
matchmodes(&numer, &denom);
switch (numer.mode) {
case IN:
checkizero(denom.val);
numer.fval = numer.val / (double)denom.val;
break;
case FL:
checkfzero(denom.fval);
numer.fval /= denom.fval;
break;
#ifdef BIGNUMS
case RT:
checkizero(denom.num);
numer.fval = ((FLOTYPE)numer.num * denom.denom) /
((FLOTYPE)numer.denom * denom.num);
break;
#endif
default: badcop();
}
}
else { /* single argument version */
switch (numer.mode) {
case IN: return lval; /* no-operation */
case FL: break; /* continue */
#ifdef BIGNUMS
case RT: numer.fval = numer.num / (FLOTYPE) numer.denom;
break;
#endif
default: badcop();
}
}
switch (which) { /* now do it! */
case 'I': MODF(numer.fval,&numer.fval); break;
case '_': numer.fval = floor(numer.fval); break;
case '^': numer.fval = ceil(numer.fval); break;
case 'r': numer.fval = floor(numer.fval + 0.5); break;
}
if (fabs(numer.fval) > (double)MAXFIX)
return cvflonum((FLOTYPE)numer.fval);
return cvfixnum((FIXTYPE)numer.fval);
#endif /* MULVALS */
}
/* unary functions */
LVAL xlognot(V) { return (unary('~')); } /* lognot */
LVAL xabs(V) { return (unary('A')); } /* abs */
LVAL xadd1(V) { return (unary('+')); } /* 1+ */
LVAL xsub1(V) { return (unary('-')); } /* 1- */
LVAL xsin(V) { return (unary('S')); } /* sin */
LVAL xcos(V) { return (unary('C')); } /* cos */
LVAL xtan(V) { return (unary('T')); } /* tan */
LVAL xexp(V) { return (unary('E')); } /* exp */
LVAL xsqrt(V) { return (unary('R')); } /* sqrt */
LVAL xfloat(V) { return (unary('F')); } /* float */
LVAL xasin(V) { return (unary('s')); } /* asin */
LVAL xacos(V) { return (unary('c')); } /* acos */
LVAL xphase(V) { return (unary('P')); } /* phase */
LVAL xfix(V) { return (unary2('I')); } /* truncate */
LVAL xfloor(V) { return (unary2('_')); } /* floor */
LVAL xceil(V) { return (unary2('^')); } /* ceiling */
LVAL xround(V) { return (unary2('r')); } /* round */
#ifdef BIGNUMS /* Added by TAA, 7/92 */
/****** needs vectorizing */
LVAL xrational(V) { return (unary('X')); } /* rational(ize) */
#endif
/* predicate - handle a predicate function */
LOCAL LVAL predicate P1C(int, fcn)
{
FLOTYPE fval;
FIXTYPE ival;
LVAL arg;
/* get the argument */
arg = xlgetarg();
xllastarg();
/* check the argument type */
if (fixp(arg)) {
ival = getfixnum(arg);
switch (fcn) {
case '-': ival = (ival < 0); break;
case 'Z': ival = (ival == 0); break;
case '+': ival = (ival > 0); break;
case 'E': ival = ((ival & 1) == 0); break;
case 'O': ival = ((ival & 1) != 0); break;
default: badiop();
}
}
else if (floatp(arg)) {
fval = getflonum(arg);
switch (fcn) {
case '-': ival = (fval < 0); break;
case 'Z': ival = (fval == 0); break;
case '+': ival = (fval > 0); break;
default: badfop();
}
}
#ifdef BIGNUMS
else if (ratiop(arg)) {
switch (fcn) {
case '-': ival = (getfixnum(getnumer(arg)) < 0); break;
case 'Z': ival = FALSE; break; /* can't be zero! */
case '+': ival = (getfixnum(getnumer(arg)) > 0); break;
default: badrop();
}
}
#endif
else
xlbadtype(arg);
/* return the result value */
return (ival ? s_true : NIL);
}
/* unary predicates */
LVAL xminusp(V) { return (predicate('-')); } /* minusp */
LVAL xzerop(V) { return (predicate('Z')); } /* zerop */
LVAL xplusp(V) { return (predicate('+')); } /* plusp */
LVAL xevenp(V) { return (predicate('E')); } /* evenp */
LVAL xoddp(V) { return (predicate('O')); } /* oddp */
/* compare - common compare function */
LOCAL LVAL compare P1C(int, fcn)
{
FIXTYPE icmp,ival,iarg;
FLOTYPE fcmp,fval,farg;
LVAL arg;
int mode;
/* get the first argument */
arg = xlgetarg();
/* set the type of the first argument */
if (fixp(arg)) {
ival = getfixnum(arg);
mode = 'I';
}
else if (floatp(arg)) {
fval = getflonum(arg);
mode = 'F';
}
#ifdef BIGNUMS
else if (ratiop(arg)) {
fval = getfixnum(getnumer(arg)) / (FLOTYPE) getfixnum(getdenom(arg));
mode = 'F';
}
#endif
else
xlbadtype(arg);
/* handle each remaining argument */
for (icmp = TRUE; icmp && moreargs(); ival = iarg, fval = farg) {
/* get the next argument */
arg = xlgetarg();
/* check its type */
if (fixp(arg)) {
switch (mode) {
case 'I':
iarg = getfixnum(arg);
break;
case 'F':
farg = (FLOTYPE)getfixnum(arg);
break;
}
}
else if (floatp(arg)) {
switch (mode) {
case 'I':
fval = (FLOTYPE)ival;
farg = getflonum(arg);
mode = 'F';
break;
case 'F':
farg = getflonum(arg);
break;
}
}
#ifdef BIGNUMS
else if (ratiop(arg)) {
switch (mode) {
case 'I':
fval = (FLOTYPE)ival;
case 'F':
farg = getfixnum(getnumer(arg)) / (FLOTYPE) getfixnum(getdenom(arg));
mode = 'F';
break;
}
}
#endif
else
xlbadtype(arg);
/* compute result of the compare */
switch (mode) {
case 'I':
icmp = ival - iarg;
switch (fcn) {
case '<': icmp = (icmp < 0); break;
case 'L': icmp = (icmp <= 0); break;
case '=': icmp = (icmp == 0); break;
case '#': icmp = (icmp != 0); break;
case 'G': icmp = (icmp >= 0); break;
case '>': icmp = (icmp > 0); break;
}
break;
case 'F':
fcmp = fval - farg;
switch (fcn) {
case '<': icmp = (fcmp < 0.0); break;
case 'L': icmp = (fcmp <= 0.0); break;
case '=': icmp = (fcmp == 0.0); break;
case '#': icmp = (fcmp != 0.0); break;
case 'G': icmp = (fcmp >= 0.0); break;
case '>': icmp = (fcmp > 0.0); break;
}
break;
}
}
/* return the result */
return (icmp ? s_true : NIL);
}
LOCAL LVAL ccompare P1C(int, which)
{
Number val, arg;
int icmp;
switch (which) {
case '=': icmp = TRUE; break;
case '#': icmp = FALSE; break;
}
readnumber(&val);
while (moreargs()) {
readnumber(&arg);
matchmodes(&val, &arg);
switch (which) {
case '=':
switch (val.mode) {
case IN: icmp = icmp && val.val == arg.val; break;
case FL: icmp = icmp && val.fval == arg.fval; break;
case CI: icmp = icmp && val.crval==arg.crval && val.cival==arg.cival;
break;
case CF: icmp = icmp && val.cfrval==arg.cfrval && val.cfival==arg.cfival;
break;
#ifdef BIGNUMS
case RT: icmp = icmp && val.num == arg.num && val.denom == arg.denom;
break;
#endif
}
break;
case '#':
switch (val.mode) {
case IN: icmp = icmp || val.val != arg.val; break;
case FL: icmp = icmp || val.fval != arg.fval; break;
case CI: icmp = icmp || val.crval!=arg.crval || val.cival!=arg.cival;
break;
case CF: icmp = icmp || val.cfrval!=arg.cfrval || val.cfival!=arg.cfival;
break;
#ifdef BIGNUMS
case RT: icmp = icmp || val.num != arg.num || val.denom != arg.denom;
break;
#endif
}
break;
}
}
return((icmp) ? s_true : NIL);
}
/* comparison functions */
LVAL xlss(V) { return (compare('<')); } /* < */
LVAL xleq(V) { return (compare('L')); } /* <= */
LVAL xequ(V) { return (ccompare('=')); } /* = */
LVAL xneq(V) { return (ccompare('#')); } /* /= */
LVAL xgeq(V) { return (compare('G')); } /* >= */
LVAL xgtr(V) { return (compare('>')); } /* > */
/***********************************************************************/
/** **/
/** 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 (newcomplex(real,cvflonum((FLOTYPE)0.0)));
}
LVAL xconjugate(V)
{
LVAL arg = xlgetarg();
xllastarg();
if (realp(arg)) return(arg);
if (! complexp(arg)) xlbadtype(arg);
if (fixp(getreal(arg)) && fixp(getimag(arg)))
return(newicomplex(getfixnum(getreal(arg)), -getfixnum(getimag(arg))));
else
return(newdcomplex(makefloat(getreal(arg)), -makefloat(getimag(arg))));
}
LVAL xrealpart(V)
{
LVAL arg = xlgetarg();
xllastarg();
if (realp(arg)) return(arg);
if (! complexp(arg)) xlbadtype(arg);
return(getreal(arg));
}
LVAL ximagpart(V)
{
LVAL arg = xlgetarg();
xllastarg();
switch (ntype(arg)) {
case FIXNUM:
#ifdef BIGNUMS
case RATIO:
#endif
return(cvfixnum((FIXTYPE) 0));
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