home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / lisp / interpre / xlispplu / sources / xlmath2.c < prev    next >
C/C++ Source or Header  |  1992-02-03  |  46KB  |  1,720 lines
  1. /* commonmath - xlisp math functions modified and augmented to         */
  2. /* correspond more closely to Common Lisp standard                     */
  3. /*      Modifications by Luke Tierney for                              */
  4. /*      XLISP-STAT 2.0 Copyright (c) 1988 by Luke Tierney              */
  5. /* xlmath - xlisp built-in arithmetic functions                        */
  6. /*      Copyright (c) 1985, by David Michael Betz                      */
  7. /*      All Rights Reserved                                            */
  8. /*      Permission is granted for unrestricted non-commercial use      */
  9.  
  10. #include "xlisp.h"
  11. #include <math.h>
  12.  
  13. #ifdef COMPLX
  14. /* These definititions used instead of those in XLMATH.C */
  15. /* Somewhat cleaned by by Tom Almy */
  16.  
  17. /* external variables */
  18. extern LVAL true;
  19.  
  20. #define IN 0
  21. #define FL 1
  22. #define CI 2
  23. #define CF 3
  24. #ifdef RATIOS
  25. #define RT 4
  26. #endif
  27.  
  28. typedef struct {
  29.   int mode;
  30.   FIXTYPE val, crval, cival;
  31. #ifdef RATIOS
  32.   FIXTYPE num, denom;
  33. #endif
  34.   FLOTYPE fval, cfrval, cfival;
  35. } Number;
  36.  
  37.  
  38. typedef struct {
  39.   double real, imag;
  40. } Complex;
  41.  
  42.  
  43. #ifdef ANSI /* local declarations */
  44. void NEAR checkizero(FIXTYPE iarg);
  45. void NEAR checkfzero(FLOTYPE farg);
  46. void NEAR badiop(void);
  47. void NEAR badfop(void);
  48. #ifdef RATIOS
  49. void NEAR badrop(void);
  50. void NEAR add_ratios(FIXTYPE *n1, FIXTYPE *d1, FIXTYPE n2, FIXTYPE d2);
  51. void NEAR mult_ratios(FIXTYPE *n1, FIXTYPE *d1, FIXTYPE n2, FIXTYPE d2);
  52. #endif
  53. void NEAR badcop(void);
  54. Complex NEAR makecomplex(LVAL x);
  55. double NEAR modulus(Complex c);
  56. Complex NEAR cart2complex(double real, double imag);
  57. LVAL NEAR cvcomplex(Complex c);
  58. Complex NEAR polar2complex(double mod, double phi);
  59. Complex NEAR csqrt(Complex c);
  60. Complex NEAR cexp(Complex c);
  61. Complex NEAR clog(Complex c);
  62. Complex NEAR cmul(Complex c1, Complex c2);
  63. Complex NEAR cexpt(Complex c1, Complex c2);
  64. Complex NEAR cadd(Complex c1, Complex c2);
  65. Complex NEAR csub(Complex c1, Complex c2);
  66. Complex NEAR cdiv(Complex c1, Complex c2);
  67. Complex NEAR csin(Complex c);
  68. Complex NEAR ccos(Complex c);
  69. Complex NEAR ctan(Complex c);
  70. Complex NEAR casin(Complex c);
  71. Complex NEAR cacos(Complex c);
  72. Complex NEAR catan(Complex c);
  73. Complex NEAR catan2(Complex n, Complex d);
  74. LVAL NEAR readnumber(Number *num);
  75. void NEAR setmode(Number *x, int mode);
  76. void NEAR matchmodes(Number *x, Number *y);
  77. LVAL NEAR lispnumber(Number *x);
  78. LVAL NEAR binary(int which);
  79. LVAL NEAR logbinary(int which);
  80. void NEAR get_mod_arg(FLOTYPE *fval, int *mode);
  81. FIXTYPE NEAR xget_gcd(FIXTYPE n, FIXTYPE m);
  82. LVAL NEAR unary(int which);
  83. LVAL NEAR unary2(int which);
  84. LVAL NEAR predicate(int fcn);
  85. LVAL NEAR compare(int fcn);
  86. LVAL NEAR ccompare(int which);
  87. LVAL NEAR xsgetreal(void);
  88. double NEAR logarithm(FLOTYPE x, FLOTYPE base, int base_supplied);
  89. #endif
  90.  
  91.  
  92. /* Error checking and messages */
  93.  
  94. /* checkizero - check for integer division by zero */
  95. LOCAL VOID NEAR checkizero(iarg)
  96.   FIXTYPE iarg;
  97. {
  98.   if (iarg == 0)
  99.   xlfail("illegal zero argument");
  100. }
  101.  
  102. /* checkfzero - check for floating point division by zero or log of zero */
  103. LOCAL VOID NEAR checkfzero(farg)
  104.   FLOTYPE farg;
  105. {
  106.   if (farg == 0.0)
  107.   xlfail("illegal zero argument");
  108. }
  109.  
  110. /* badiop - bad integer operation */
  111. LOCAL VOID NEAR badiop()
  112. {
  113.   xlfail("bad integer operation");
  114. }
  115.  
  116. /* badfop - bad floating point operation */
  117. LOCAL VOID NEAR badfop()
  118. {
  119.   xlfail("bad floating point operation");
  120. }
  121.  
  122. #ifdef RATIOS
  123. /* badrop - bad ratio operation */
  124. LOCAL VOID NEAR badrop()
  125. {
  126.   xlfail("bad ratio operation");
  127. }
  128. #endif
  129.  
  130. /* badcop - bad complex number operation */
  131. LOCAL VOID NEAR badcop()
  132. {
  133.   xlfail("bad complex number operation");
  134. }
  135.  
  136. /* TAA MOD badarg() removed, using preexisting xlbadtype which
  137.    gives same message */
  138.  
  139. /* complex - Complex number functions  */
  140.  
  141. /*TAA MOD--do inline with libm call*/
  142. #define phase(c) ((c).imag==0.0 && (c).real == 0.0 ? 0 : atan2((c).imag,(c).real))
  143.  
  144. /* TAA Mod--do inline, old test for complexp inline as appropriate */
  145. /*  badcomplex() eliminated since it didn't make sense (other numereric
  146.     types were ok), so returned error is now from xlbadtype */
  147.  
  148. #define realpart(x) (getelement(x, 0))
  149.  
  150. #define imagpart(x) (getelement(x, 1))
  151.  
  152.  
  153. LOCAL Complex NEAR makecomplex(x)
  154.         LVAL x;
  155. {
  156.   Complex c;
  157.  
  158.   if (numberp(x)) {
  159.     c.real = makefloat(x);
  160.     c.imag = 0.0;
  161.   }
  162.   else if (complexp(x)) {
  163.     c.real = makefloat(realpart(x));
  164.     c.imag = makefloat(imagpart(x));
  165.   }
  166.   else xlerror("not a number", x);
  167.   return(c);
  168. }
  169.  
  170. LOCAL double NEAR modulus(c)
  171.         Complex c;
  172. {
  173.   return(sqrt(c.real * c.real + c.imag * c.imag));
  174. }
  175.  
  176. LOCAL Complex NEAR cart2complex(real, imag)
  177.         double real, imag;
  178. {
  179.   Complex val;
  180.   val.real = real;
  181.   val.imag = imag;
  182.   return(val);
  183. }
  184.  
  185. LOCAL LVAL NEAR cvcomplex(c)
  186.         Complex c;
  187. {
  188.   return(newdcomplex(c.real, c.imag));
  189. }
  190.  
  191. LOCAL Complex NEAR polar2complex(mod, phi)
  192.         double mod, phi;
  193. {
  194.   Complex val;
  195.   double cs, sn;
  196.  
  197.   if (phi == 0) {
  198.     cs = 1.0;
  199.     sn = 0.0;
  200.   }
  201.   else if (phi == PI / 2) {
  202.     cs = 0.0;
  203.     sn = 1.0;
  204.   }
  205.   else if (phi == PI) {
  206.     cs = -1.0;
  207.     sn = 0.0;
  208.   }
  209.   else if (phi == -PI / 2) {
  210.     cs = 0.0;
  211.     sn = -1.0;
  212.   }
  213.   else {
  214.     cs = cos(phi);
  215.     sn = sin(phi);
  216.   }
  217.   val.real = mod * cs;
  218.   val.imag = mod * sn;
  219.   return(val);
  220. }
  221.  
  222. LOCAL Complex NEAR csqrt(c)
  223.         Complex c;
  224. {
  225.   return(polar2complex(sqrt(modulus(c)), phase(c) / 2));
  226. }
  227.  
  228. LOCAL Complex NEAR cexp(c)
  229.         Complex c;
  230. {
  231.   return(polar2complex(exp(c.real), c.imag));
  232. }
  233.  
  234. LOCAL Complex NEAR clog(c)
  235.         Complex c;
  236. {
  237.   double mod;
  238.  
  239.   mod = modulus(c);
  240.   checkfzero(mod);
  241.   return(cart2complex(log(mod), phase(c)));
  242. }
  243.  
  244. LOCAL Complex NEAR cmul(c1, c2)
  245.         Complex c1, c2;
  246. {
  247. #if 0   /* This is rediculous, says TAA */
  248.         /* why pay the cost for the two conversions? */
  249.   double m1, m2, p1, p2;
  250.  
  251.   m1 = modulus(c1);
  252.   p1 = phase(c1);
  253.   m2 = modulus(c2);
  254.   p2 = phase(c2);
  255.   return(polar2complex(m1 * m2, p1 + p2));
  256. #else
  257.     Complex val;
  258.  
  259.     val.real = c1.real * c2.real - c1.imag * c2.imag;
  260.     val.imag = c1.imag * c2.real + c1.real * c2.imag;
  261.     return val;
  262. #endif
  263. }
  264.  
  265. LOCAL Complex NEAR cexpt(cb, cp)
  266.         Complex cb, cp;
  267. {
  268.   if (modulus(cp) == 0.0) return(cart2complex(1.0, 0.0));
  269.   else  return(cexp(cmul(clog(cb), cp)));
  270. }
  271.  
  272. LOCAL Complex NEAR cadd(c1, c2)
  273.         Complex c1, c2;
  274. {
  275.   return(cart2complex(c1.real + c2.real, c1.imag + c2.imag));
  276. }
  277.  
  278. LOCAL Complex NEAR csub(c1, c2)
  279.         Complex c1, c2;
  280. {
  281.   return(cart2complex(c1.real - c2.real, c1.imag - c2.imag));
  282. }
  283.  
  284. LOCAL Complex NEAR cdiv(c1, c2)
  285.         Complex c1, c2;
  286. {
  287. #if 0
  288.   double m1, m2, p1, p2;
  289.  
  290.   m1 = modulus(c1);
  291.   p1 = phase(c1);
  292.   m2 = modulus(c2);
  293.   p2 = phase(c2);
  294.   checkfzero(m2);
  295.   return(polar2complex(m1 / m2, p1 - p2));
  296. #else /* replace with code that is faster */
  297.     double ratio, temp;
  298.  
  299.     if (fabs(c2.real) > fabs(c2.imag)) {
  300.         ratio = c2.imag / c2.real;
  301.         temp = c2.real + ratio*c2.imag;
  302.         return  cart2complex((c1.real + c1.imag*ratio)/temp,
  303.                              (c1.imag - c1.real*ratio)/temp);
  304.     }
  305.     else {
  306.         checkfzero(c2.imag);
  307.         ratio = c2.real / c2.imag;
  308.         temp = c2.imag + ratio*c2.real;
  309.         return  cart2complex ((c1.real*ratio + c1.imag)/temp,
  310.                               (c1.imag*ratio - c1.real)/temp);
  311.     }
  312. #endif
  313. }
  314.  
  315. LOCAL Complex NEAR csin(c)
  316.         Complex c;
  317. {
  318.   Complex x1, x2, val;
  319.  
  320.   x1 = cart2complex(-c.imag, c.real);
  321.   x2 = cart2complex(c.imag, -c.real);
  322.   val = csub(cexp(x1), cexp(x2));
  323.   return(cart2complex(val.imag / 2.0, -val.real / 2.0));
  324. }
  325.  
  326. LOCAL Complex NEAR ccos(c)
  327.         Complex c;
  328. {
  329.   Complex x1, x2, val;
  330.  
  331.   x1 = cart2complex(-c.imag, c.real);
  332.   x2 = cart2complex(c.imag, -c.real);
  333.   val = cadd(cexp(x1), cexp(x2));
  334.   return(cart2complex(val.real / 2.0, val.imag / 2.0));
  335. }
  336.  
  337. LOCAL Complex NEAR ctan(c)
  338.         Complex c;
  339. {
  340.   Complex e1, e2, val;
  341.  
  342.   e1 = cexp(cart2complex(-c.imag, c.real));
  343.   e2 = cexp(cart2complex(c.imag, -c.real));
  344.   val = cdiv(csub(e1, e2), cadd(e1, e2));
  345.   return(cart2complex(val.imag, -val.real));
  346. }
  347.  
  348. LOCAL Complex NEAR casin(c)
  349.         Complex c;
  350. {
  351.   Complex sx, ix, val;
  352.  
  353.   sx = cmul(c, c);
  354.   sx = csqrt(cart2complex(1.0 - sx.real, - sx.imag));
  355.   ix = cart2complex(-c.imag, c.real);
  356.   val = clog(cadd(ix, sx));
  357.   return(cart2complex(val.imag, -val.real));
  358. }
  359.  
  360. LOCAL Complex NEAR cacos(c)
  361.         Complex c;
  362. {
  363.   Complex sx, val;
  364.  
  365.   sx = cmul(c, c);
  366.   sx = csqrt(cart2complex(1.0 - sx.real, - sx.imag));
  367.   sx = cart2complex(-sx.imag, sx.real);
  368.   val = clog(cadd(c, sx));
  369.   return(cart2complex(val.imag, -val.real));
  370. }
  371.  
  372. LOCAL Complex NEAR catan(c)
  373.         Complex c;
  374. {
  375. #if 0   /* This has been redefined in Jan 1989 */
  376.   Complex sx, ix, val, one;
  377.  
  378.   sx = cmul(c, c);
  379.   sx = cart2complex(1.0 + sx.real, sx.imag);
  380.   one = cart2complex(1.0, 0.0);
  381.   sx = csqrt(cdiv(one, sx));
  382.   ix = cadd(one, cart2complex(-c.imag, c.real));
  383.   val = clog(cmul(ix, sx));
  384.   return(cart2complex(val.imag, -val.real));
  385. #else
  386.   Complex sx, ix;
  387.  
  388.   sx = clog(cart2complex(1.0-c.imag, c.real));
  389.   ix = clog(cart2complex(1.0+c.imag, -c.real));
  390.   sx.real -= ix.real;   /* complex addition w.o. subroutine call */
  391.   sx.imag -= ix.imag;
  392.   return cdiv(sx,cart2complex(0.0, 2.0));
  393. #endif
  394. }
  395.  
  396. LOCAL Complex NEAR catan2(n,d)  /* by Tom Almy, and a kludge at that */
  397.  
  398. Complex n, d;
  399. {
  400.     double tmp;
  401.  
  402.     tmp = modulus(d);
  403.     if (tmp == 0 || modulus(n)/tmp > 1e50 ) { /* worst case */
  404.         tmp = phase(n) - phase(d);
  405.         if (tmp < -PI) tmp += 2.0*PI;
  406.         else if (tmp > PI) tmp -= 2.0*PI;
  407.         return cart2complex(fabs(tmp) > PI/2.0 ? -PI/2.0 : PI/2.0 ,0.0);
  408.     }
  409.     n = cdiv(n,d);  /* best case */
  410.     return (catan(n));
  411. }
  412.  
  413. /* Helper functions */
  414.  
  415. LOCAL LVAL NEAR readnumber(number)
  416.         Number *number;
  417. {
  418.   LVAL arg = xlgetarg(), real, imag;
  419.  
  420.   if (fixp(arg)) {
  421.       number->mode = IN;
  422.       number->val = getfixnum(arg);
  423.   }
  424.   else if (floatp(arg)) {
  425.       number->mode = FL;
  426.       number->fval = getflonum(arg);
  427.   }
  428. #ifdef RATIOS
  429.   else if (ratiop(arg)) {
  430.       number->mode = RT;
  431.       number->num = getnumer(arg);
  432.       number->denom = getdenom(arg);
  433.   }
  434. #endif
  435.   else if (complexp(arg)) {
  436.       real = realpart(arg);
  437.       imag = imagpart(arg);
  438.       if (fixp(real)) {
  439.           number->mode = CI;
  440.           number->crval = getfixnum(real);
  441.           number->cival = getfixnum(imag);
  442.       }
  443.       else {
  444.           number->mode = CF;
  445.           number->cfrval = makefloat(real);
  446.           number->cfival = makefloat(imag);
  447.       }
  448.   }
  449.   else xlerror("not a number", arg);
  450.   return(arg);
  451. }
  452.  
  453. LOCAL VOID NEAR setmode(x, mode)
  454.         Number *x;
  455.         int mode;
  456. {
  457.   switch (mode) {
  458. #ifdef RATIOS
  459.   case RT:
  460.     if (x->mode != IN) return;
  461.     x->mode = mode;
  462.     x->num = x->val;
  463.     x->denom = 1;
  464.     break;
  465. #endif
  466.   case FL:
  467.     if (x->mode == IN) {
  468.         x->mode = mode;
  469.         x->fval = x->val;
  470.     }
  471. #ifdef RATIOS
  472.     else if (x->mode == RT) {
  473.         x->mode = mode;
  474.         x->fval = x->num / (FLOTYPE) x->denom;
  475.     }
  476. #endif
  477.     else return;
  478.     break;
  479.   case CI:
  480.     if (x->mode != IN) return;
  481.     x->mode = mode;
  482.     x->crval = x->val;
  483.     x->cival = 0;
  484.     break;
  485.   case CF:
  486.     switch (x->mode) {
  487.     case IN:
  488.       x->mode = mode;
  489.       x->cfrval = x->val;
  490.       x->cfival = 0.0;
  491.       break;
  492. #ifdef RATIOS
  493.     case RT:
  494.       x->mode = mode;
  495.       x->cfrval = x->num / (FLOTYPE) x->denom;
  496.       x->cfival = 0.0;
  497.       break;
  498. #endif
  499.     case FL:
  500.       x->mode = mode;
  501.       x->cfrval = x->fval;
  502.       x->cfival = 0.0;
  503.       break;
  504.     case CI:
  505.       x->mode = mode;
  506.       x->cfrval = x->crval;
  507.       x->cfival = x->cival;
  508.       break;
  509.     }
  510.     break;
  511.   }
  512. }
  513.  
  514. LOCAL VOID NEAR matchmodes(x, y)
  515.         Number *x, *y;
  516. {
  517.   int mode = x->mode;
  518.   switch (mode) {
  519.   case IN: mode = y->mode; break;
  520.   case FL: if (y->mode == CI || y->mode == CF) mode = CF; break;
  521. #ifdef RATIOS
  522.   case CI: if (y->mode == FL || y->mode == CF || y->mode == RT)
  523.       mode = CF;
  524.   break;
  525. #else
  526.   case CI: if (y->mode == FL || y->mode == CF) mode = CF; break;
  527. #endif
  528.   case CF: break;
  529. #ifdef RATIOS
  530.   case RT: if (y->mode == CI) mode = CF;
  531.            else if (y->mode != IN) mode = y->mode;
  532.            break;
  533. #endif
  534.   }
  535.   if (x->mode != mode) setmode(x, mode);
  536.   if (y->mode != mode) setmode(y, mode);
  537. }
  538.  
  539. LOCAL LVAL NEAR lispnumber(x)
  540.         Number *x;
  541. {
  542.   switch (x->mode) {
  543.   case IN: return(cvfixnum(x->val));
  544.   case FL: return(cvflonum(x->fval));
  545.   case CI: return(newicomplex(x->crval, x->cival));
  546.   case CF: return(newdcomplex(x->cfrval, x->cfival));
  547. #ifdef RATIOS
  548.   case RT: return(cvratio (x->num, x->denom));
  549. #endif
  550.   }
  551.   return NIL; /* avoid warning messages */
  552. }
  553.  
  554. #ifdef RATIOS
  555. LOCAL VOID NEAR add_ratios (n1, d1, n2, d2)
  556.   FIXTYPE *n1, *d1, n2, d2;
  557. {
  558.     FIXTYPE lcm;
  559.  
  560.     lcm = *d1 * (d2 / xget_gcd(*d1, d2));
  561.     *n1 = *n1 * (lcm / *d1)  + n2 * (lcm / d2);
  562.     *d1 = lcm;
  563.     return;
  564. }
  565.  
  566. LOCAL VOID NEAR mult_ratios (n1, d1, n2, d2)
  567.   FIXTYPE *n1, *d1, n2, d2;
  568. {
  569.     FIXTYPE gcd;
  570.  
  571.     if ((*n1 *= n2) == 0) return;
  572.     gcd = xget_gcd (*n1, *d1 *= d2);
  573.     *n1 /= gcd;
  574.     *d1 /= gcd;
  575.     return;
  576. }
  577. #endif
  578.  
  579. LOCAL LVAL NEAR binary(which)
  580.         int which;
  581. {
  582.   LVAL larg;
  583.   Number val, arg;
  584.   FIXTYPE rtemp, itemp;
  585.   FLOTYPE frtemp, fitemp;
  586.  
  587.   if (xlargc == 1 && (which == '-' || which == '/')) {
  588.     val.mode = IN;
  589.     switch (which) {
  590.     case '-': val.val = 0; break;
  591.     case '/': val.val = 1; break;
  592.     }
  593.   }
  594.   else larg = readnumber(&val);
  595.   while (moreargs()) {
  596.     larg = readnumber(&arg);
  597.     matchmodes(&val, &arg);
  598.     switch (which) {
  599.     case '+':
  600.       switch (val.mode) {
  601.       case IN: val.val   += arg.val;  break;
  602.       case FL: val.fval  += arg.fval; break;
  603.       case CI: val.crval += arg.crval;   val.cival += arg.cival;   break;
  604.       case CF: val.cfrval += arg.cfrval; val.cfival += arg.cfival; break;
  605. #ifdef RATIOS
  606.       case RT:
  607.          add_ratios (&val.num, &val.denom, arg.num, arg.denom);
  608.          break;
  609. #endif
  610.       }
  611.       break;
  612.     case '-':
  613.       switch (val.mode) {
  614.       case IN: val.val   -= arg.val;  break;
  615.       case FL: val.fval  -= arg.fval; break;
  616.       case CI: val.crval -= arg.crval;   val.cival -= arg.cival;   break;
  617.       case CF: val.cfrval -= arg.cfrval; val.cfival -= arg.cfival; break;
  618. #ifdef RATIOS
  619.       case RT:
  620.          add_ratios (&val.num, &val.denom, -arg.num, arg.denom);
  621.          break;
  622. #endif
  623.       }
  624.       break;
  625.     case '*':
  626.       switch (val.mode) {
  627.       case IN: val.val   *= arg.val;  break;
  628.       case FL: val.fval  *= arg.fval; break;
  629.       case CI:
  630.         rtemp = val.crval * arg.crval - val.cival * arg.cival;
  631.         itemp = val.cival * arg.crval + val.crval * arg.cival;
  632.         val.crval = rtemp; val.cival = itemp;
  633.         break;
  634.       case CF:
  635.         frtemp = val.cfrval * arg.cfrval - val.cfival * arg.cfival;
  636.         fitemp = val.cfival * arg.cfrval + val.cfrval * arg.cfival;
  637.         val.cfrval = frtemp; val.cfival = fitemp;
  638.         break;
  639. #ifdef RATIOS
  640.       case RT:
  641.          mult_ratios (&val.num, &val.denom, arg.num, arg.denom);
  642.          break;
  643. #endif
  644.       }
  645.       break;
  646.     case '/':
  647.       switch (val.mode) {
  648.       case IN:
  649.         checkizero(arg.val);
  650. #ifdef RATIOS
  651.         setmode(&val, RT);
  652.         val.denom = arg.val;
  653.         break;
  654. #else
  655.         if (val.val % arg.val == 0) {
  656.           val.val /= arg.val;
  657.           break;
  658.         }
  659.         else {
  660.           setmode(&val, FL);
  661.           setmode(&arg, FL);
  662.         }
  663.         /* drop through */
  664. #endif
  665.       case FL:
  666.         checkfzero(arg.fval);
  667.         val.fval /= arg.fval;
  668.         break;
  669.       case CI:
  670.         setmode(&val, CF);
  671.         setmode(&arg, CF);
  672.         /* drop through */
  673.       case CF:
  674. #if 0   /* we can do better */
  675.     { double magn;
  676.         magn = arg.cfrval * arg.cfrval + arg.cfival * arg.cfival;
  677.         checkfzero(magn);
  678.         frtemp = (val.cfrval * arg.cfrval + val.cfival * arg.cfival) / magn;
  679.         fitemp = (val.cfival * arg.cfrval - val.cfrval * arg.cfival) / magn;
  680.         val.cfrval = frtemp; val.cfival = fitemp;
  681.         break;
  682.     }
  683. #else
  684.     { double ratio,temp;
  685.         if (fabs(arg.cfrval) > fabs(arg.cfival)) {
  686.             ratio = arg.cfival / arg.cfrval;
  687.             temp = arg.cfrval + ratio*arg.cfival;
  688.             frtemp = (val.cfrval + val.cfival*ratio)/temp;
  689.             fitemp = (val.cfival - val.cfrval*ratio)/temp;
  690.         }
  691.         else {
  692.             checkfzero(arg.cfival);
  693.             ratio = arg.cfrval / arg.cfival;
  694.             temp = arg.cfival + ratio*arg.cfrval;
  695.             frtemp = (val.cfrval*ratio + val.cfival)/temp;
  696.             fitemp = (val.cfival*ratio - val.cfrval)/temp;
  697.         }
  698.         val.cfrval = frtemp; val.cfival = fitemp;
  699.         break;
  700.     }
  701. #endif
  702. #ifdef RATIOS
  703.       case RT:
  704.          checkizero (arg.num);
  705.          mult_ratios (&val.num, &val.denom, arg.denom, arg.num);
  706.          break;
  707. #endif
  708.       }
  709.       break;
  710.     case 'M':
  711.       switch (val.mode) {
  712.       case IN: val.val  = (val.val > arg.val)   ? val.val  : arg.val;  break;
  713.       case FL: val.fval = (val.fval > arg.fval) ? val.fval : arg.fval; break;
  714. #ifdef RATIOS
  715.       case RT:
  716.          if ((val.num * arg.denom) < (arg.num * val.denom)) {
  717.            val.num = arg.num;
  718.            val.denom = arg.denom;
  719.            }
  720.          break;
  721. #endif
  722.       default: xlbadtype(larg);
  723.       }
  724.       break;
  725.     case 'm':
  726.       switch (val.mode) {
  727.       case IN: val.val  = (val.val < arg.val)   ? val.val  : arg.val;  break;
  728.       case FL: val.fval = (val.fval < arg.fval) ? val.fval : arg.fval; break;
  729. #ifdef RATIOS
  730.       case RT:
  731.          if ((val.num * arg.denom) > (arg.num * val.denom)) {
  732.            val.num = arg.num;
  733.            val.denom = arg.denom;
  734.            }
  735.          break;
  736. #endif
  737.       default: xlbadtype(larg);
  738.       }
  739.       break;
  740.     }
  741.   }
  742.   return(lispnumber(&val));
  743. }
  744.  
  745.  
  746. /* This has been completely rewritten by Tom Almy to handle floating point
  747.    arguments */
  748.  
  749. LVAL xrem()
  750. {
  751.     Number numer, div;
  752.     double ftemp;
  753.  
  754.     readnumber(&numer);
  755.     readnumber(&div);
  756.     xllastarg();
  757.  
  758.     matchmodes(&numer, &div);
  759.  
  760.     switch (numer.mode) {
  761.     case IN:    checkizero(div.val);
  762.                 return (cvfixnum((FIXTYPE) numer.val % div.val));
  763.     case FL:    checkfzero(div.fval);
  764.                 modf(numer.fval / div.fval, &ftemp);
  765.                 return (cvflonum((FLOTYPE)(numer.fval - ftemp*div.fval)));
  766. #ifdef RATIOS
  767.     case RT:    checkizero(div.num);
  768.                 mult_ratios(&numer.num, &numer.denom, div.denom, div.num);
  769.                 numer.num %= numer.denom;   /* get remainder */
  770.                 mult_ratios(&numer.num, &numer.denom, div.num, div.denom);
  771.                 return (cvratio(numer.num, numer.denom));
  772. #endif
  773.     }
  774.     badcop();
  775.     return NIL; /* fool compiler into not giving warning */
  776. }
  777.  
  778.  
  779. LVAL xash() /* arithmetic shift left */
  780. {
  781.     FIXTYPE arg, val;
  782.  
  783.     arg = getfixnum(xlgafixnum());
  784.     val = getfixnum(xlgafixnum());
  785.     xllastarg();
  786.  
  787.     return cvfixnum(val < 0 ? arg >> -val : arg << val);
  788. }
  789.  
  790.  
  791. LOCAL LVAL NEAR logbinary(which)
  792.         int which;
  793. {
  794.   FIXTYPE val, arg; /* TAA Mod -- was int */
  795.  
  796.   switch (which) {
  797.   case '&': val = -1; break;
  798.   case '|': val =  0; break;
  799.   case '^': val =  0; break;
  800.   }
  801.   while (moreargs()) {
  802.     arg = getfixnum(xlgafixnum());
  803.     switch (which) {
  804.     case '&': val &= arg; break;
  805.     case '|': val |= arg; break;
  806.     case '^': val ^= arg; break;
  807.     }
  808.   }
  809.   return(cvfixnum((FIXTYPE) val));
  810. }
  811.  
  812. /* binary functions */
  813. /* TAA fix allowing (+) */
  814. LVAL xadd()    { return (moreargs()?binary('+'):cvfixnum((FIXTYPE)0)); }
  815. LVAL xsub()    { return (binary('-')); } /* - */
  816. /* TAA fix allowing (*) */
  817. LVAL xmul()    { return (moreargs()?binary('*'):cvfixnum((FIXTYPE)1)); }
  818. LVAL xdiv()    { return (binary('/')); } /* / */
  819. LVAL xmin()    { return (binary('m')); } /* min */
  820. LVAL xmax()    { return (binary('M')); } /* max */
  821. LVAL xlogand() { return (logbinary('&')); } /* logand */
  822. LVAL xlogior() { return (logbinary('|')); } /* logior */
  823. LVAL xlogxor() { return (logbinary('^')); } /* logxor */
  824.  
  825.  
  826. /* New by Tom Almy */
  827.  
  828. LVAL xmod()
  829. {
  830.     Number numer, div;
  831.  
  832.     readnumber(&numer);
  833.     readnumber(&div);
  834.     xllastarg();
  835.  
  836.     matchmodes(&numer, &div);
  837.  
  838.     switch (numer.mode) {
  839.     case IN:    checkizero(div.val);
  840.                 return(cvfixnum(numer.val -
  841.                     div.val*(FIXTYPE)floor(numer.val/(double)div.val)));
  842.     case FL:    checkfzero(div.fval);
  843.                 return(cvflonum((FLOTYPE)(numer.fval -
  844.                     div.fval*floor((double)(numer.fval/div.fval)))));
  845. #ifdef RATIOS
  846.     case RT:    checkizero(div.num);
  847.                 mult_ratios(&numer.num, &numer.denom, div.denom, div.num);
  848.                 numer.num = numer.num - numer.denom*
  849.                     (FIXTYPE)floor(numer.num/(double)numer.denom);
  850.                 mult_ratios(&numer.num, &numer.denom, div.num, div.denom);
  851.                 return (cvratio(numer.num, numer.denom));
  852. #endif
  853.     }
  854.     badcop();
  855.     return NIL; /* fool compiler into not giving warning */
  856. }
  857.  
  858. LVAL xexpt()    /* modified for ratios by TAA */
  859. {
  860.   LVAL base, power;
  861.   int bsign, psign;
  862.   FIXTYPE b, p, val;
  863.   FLOTYPE fb, fp, fval;
  864.  
  865.   base = xlgetarg();
  866.   power = xlgetarg();
  867.   xllastarg();
  868.  
  869.   if (fixp(base) && fixp(power)) {
  870.     b = getfixnum(base);
  871.     p = getfixnum(power);
  872.     if (p == 0) return(cvfixnum((FIXTYPE) 1));
  873.     if (b == 0 && p > 0) return(cvfixnum((FIXTYPE) 0));
  874.     checkizero(b);
  875.     if ((bsign = (b < 0)) != 0) b = -b;
  876.     if ((psign = (p < 0)) != 0) p = -p;
  877.  
  878.     fval = floor(pow((double) b, (double) p) + 0.1); /* to get integer right */
  879.     if (bsign && (p & 1)) fval = -fval;
  880.     if (!psign) {
  881.       val = fval;
  882.       if (val == fval) return(cvfixnum(val));
  883.       else return(cvflonum((FLOTYPE) fval));    /* to handle precision for large results */
  884.     }
  885.     else {
  886. #ifdef RATIOS
  887.       val = fval;
  888.       if (val == fval) return (cvratio(1, val));
  889.       else return (cvflonum((FLOTYPE) 1.0 / fval));
  890. #else
  891.       checkfzero(fval);
  892.       return(cvflonum((FLOTYPE) 1.0 / fval));
  893. #endif
  894.     }
  895.   }
  896. #ifdef RATIOS
  897.   else if (ratiop(base) && fixp(power)) {
  898.     FIXTYPE vald;   /* integer of denominator result */
  899.     FLOTYPE fvald;  /* denominator result */
  900.     b = getnumer(base); /* start with just the numerator */
  901.     p = getfixnum(power);
  902.     if (p == 0) return(cvfixnum((FIXTYPE) 1));
  903.     if ((bsign = (b < 0)) != 0) b = -b;
  904.     if ((psign = (p < 0)) != 0) p = -p;
  905.  
  906.     fval = floor(pow((double) b, (double) p) + 0.1); /* to get integer right */
  907.     if (bsign && (p & 1)) fval = -fval;
  908.  
  909.     fvald = floor(pow((double) getdenom(base), (double) p) + 0.1);
  910.  
  911.     val = fval;
  912.     vald = fvald;
  913.  
  914.     if (!psign) {
  915.       if (val == fval && vald == fvald) /* will fit in ratio */
  916.           return (cvratio(val, vald));
  917.       else return (cvflonum(fval/fvald));
  918.     }
  919.     else {
  920.       if (val == fval && vald == fvald) /* will fit in ratio */
  921.           return (cvratio(vald, val));
  922.       else return (cvflonum(fvald/fval));
  923.     }
  924.   }
  925. #endif
  926.   else if (floatp(base) && fixp(power)) {
  927.     fb = getflonum(base);
  928.     p = getfixnum(power);
  929.     if (p == 0) return(cvflonum((FLOTYPE) 1.0)); /* TAA MOD - used to return
  930.                                     fixnum 1, but CL says should be flonum */
  931.     if (fb == 0.0 && p > 0) return(cvflonum((FLOTYPE) 0.0));
  932.     checkfzero(fb);
  933.     if ((bsign = (fb < 0.0)) != 0) fb = -fb;
  934.     if ((psign = (p < 0)) != 0) p = -p;
  935.     fval = pow((double) fb, (double) p);
  936.     if (bsign && (p & 1)) fval = -fval;
  937.     if (!psign) return(cvflonum((FLOTYPE) fval));
  938.     else {
  939.       checkfzero(fval);
  940.       return(cvflonum((FLOTYPE) 1.0 / fval));
  941.     }
  942.   }
  943. #ifdef RATIOS
  944.   else if (numberp(base) && (ratiop(power) || floatp(power)))
  945. #else
  946.   else if (numberp(base) && floatp(power))
  947. #endif
  948.   {
  949.     fb = makefloat(base);
  950. #ifdef RATIOS
  951.     fp = ratiop(power) ?
  952.         getnumer(power)/(FLOTYPE)getdenom(power) : getflonum(power);
  953. #else
  954.     fp = getflonum(power);
  955. #endif
  956.     if (fp == 0.0) return(cvflonum((FLOTYPE) 1.0));
  957.     if (fb == 0.0 && fp > 0.0) return(cvflonum((FLOTYPE) 0.0));
  958.     if (fb < 0.0)
  959.       return(cvcomplex(cexpt(makecomplex(base), makecomplex(power))));
  960.     if ((psign = (fp < 0.0)) != 0) fp = -fp;
  961.     fval = pow((double) fb, (double) fp);
  962.     if (!psign) return(cvflonum((FLOTYPE) fval));
  963.     else {
  964.       checkfzero(fval);
  965.       return(cvflonum((FLOTYPE) 1.0 / fval));
  966.     }
  967.   }
  968.   else if (complexp(base) || complexp(power))
  969.     return(cvcomplex(cexpt(makecomplex(base), makecomplex(power))));
  970.   else xlfail("bad argument type(s)");
  971.   return NIL; /* avoid compiler warnings */
  972. }
  973.  
  974. /* arc tangent -- Tom Almy */
  975. LVAL xatan()
  976. {
  977.     Number numer, denom;
  978.     LVAL lnum, ldenom;
  979.     Complex cnum, cdenom;
  980.  
  981.     lnum = readnumber(&numer);
  982.  
  983.     if (moreargs()) {
  984.         ldenom = readnumber(&denom);
  985.         xllastarg();
  986.         matchmodes(&numer, &denom);
  987.  
  988.         switch (numer.mode) {
  989.             case IN:
  990.                 numer.fval = numer.val; denom.fval = denom.val;
  991.             case FL:
  992.                 return (cvflonum((FLOTYPE)atan2(numer.fval, denom.fval)));
  993. #ifdef RATIOS
  994.             case RT:
  995.                 return (cvflonum((FLOTYPE)atan2(numer.num/(FLOTYPE)numer.denom,
  996.                     denom.num/(FLOTYPE)denom.denom)));
  997. #endif
  998.             default: /* complex */
  999.                 cnum = makecomplex(lnum);
  1000.                 cdenom = makecomplex(ldenom);
  1001.                 return (cvcomplex(catan2(cnum,cdenom)));
  1002.         }
  1003.     }
  1004.     else {
  1005.         switch (numer.mode) {
  1006.             case IN:
  1007.                 numer.fval = numer.val;
  1008.             case FL:
  1009.                 return (cvflonum((FLOTYPE)atan(numer.fval)));
  1010. #ifdef RATIOS
  1011.             case RT:
  1012.                 return (cvflonum((FLOTYPE)atan(numer.num/(FLOTYPE)numer.denom)));
  1013. #endif
  1014.             default: /* complex */
  1015.                 cnum = makecomplex(lnum);
  1016.                 return (cvcomplex(catan(cnum)));
  1017.         }
  1018.     }
  1019. }
  1020.  
  1021.  
  1022.  
  1023. /* two argument logarithm */
  1024. LOCAL double NEAR logarithm(x, base, base_supplied)
  1025.      FLOTYPE x, base;
  1026.      int base_supplied;
  1027. {
  1028.   double lbase;
  1029.   if (x <= 0.0) xlfail("logarithm of a nonpositive number");
  1030.   if (base_supplied) {
  1031.     if (base <= 0.0) xlfail("logarithm to a nonpositive base");
  1032.     else {
  1033.       lbase = log(base);
  1034.       if (lbase == 0.0) xlfail("logarith to a unit base");
  1035.       else return((log(x)/lbase));
  1036.     }
  1037.   }
  1038.   return (log(x));
  1039. }
  1040.  
  1041. LVAL xlog()
  1042. {
  1043.   LVAL arg, base;
  1044.   int base_supplied = FALSE;
  1045.   double fx, fb;
  1046.  
  1047.   arg = xlgetarg();
  1048.   if (moreargs()) {
  1049.     base_supplied = TRUE;
  1050.     base = xlgetarg();
  1051.   }
  1052.   if (base_supplied) {
  1053.     if (numberp(arg) && numberp(base)) {
  1054.       fx = makefloat(arg);
  1055.       fb = makefloat(base);
  1056.       if (fx <= 0.0 || fb <= 0.0)
  1057.         return(cvcomplex(cdiv(clog(makecomplex(arg)), clog(makecomplex(base)))));
  1058.       else return(cvflonum((FLOTYPE) logarithm(fx, fb, TRUE)));
  1059.     }
  1060.     else if ((numberp(arg) && complexp(base))
  1061.              || (complexp(arg) && numberp(base))
  1062.              || (complexp(arg) && complexp(base)))
  1063.       return(cvcomplex(cdiv(clog(makecomplex(arg)), clog(makecomplex(base)))));
  1064.     else xlfail("bad argument type(s)");
  1065.   }
  1066.   else {
  1067.     if (numberp(arg)) {
  1068.       fx = makefloat(arg);
  1069.       if (fx <= 0.0) return(cvcomplex(clog(makecomplex(arg))));
  1070.       else return(cvflonum((FLOTYPE) logarithm(fx, 0.0, FALSE)));
  1071.     }
  1072.     else if (complexp(arg))
  1073.       return(cvcomplex(clog(makecomplex(arg))));
  1074.     else xlfail("bad argument type(s)");
  1075.   }
  1076.   return NIL; /* avoid compiler warnings */
  1077. }
  1078.  
  1079. /* TAA Mod to return FIXTYPE */
  1080. LOCAL FIXTYPE NEAR xget_gcd(n,m)                 /* euclid's algorith */
  1081. FIXTYPE n, m;
  1082. {
  1083.    FIXTYPE r;
  1084.  
  1085.    for (;;) {
  1086.       r = m % n;
  1087.       if (r == (FIXTYPE) 0)
  1088.         break;
  1089.       m = n;
  1090.       n = r;
  1091.    }
  1092.  
  1093.    return (n);
  1094. }
  1095.  
  1096. /* xgcd - greatest common divisor */
  1097. LVAL xgcd()
  1098. {
  1099.   FIXTYPE m,n;
  1100.   LVAL arg;
  1101.  
  1102.   if (!moreargs())                  /* check for identity case */
  1103.     return (cvfixnum((FIXTYPE)0));
  1104.   arg = xlgafixnum();
  1105.   n = getfixnum(arg);
  1106.   if (n < (FIXTYPE)0) n = -n;           /* absolute value */
  1107.   while (moreargs()) {
  1108.     arg = xlgafixnum();
  1109.     m = getfixnum(arg);
  1110.     if (m == 0 || n == 0) xlfail("zero argument");
  1111.     if (m < (FIXTYPE)0) m = -m;     /* absolute value */
  1112.  
  1113.     n = xget_gcd(n,m);
  1114.   }
  1115.   return (cvfixnum(n));
  1116. }
  1117.  
  1118. LVAL xlcm()                         /* added by kcw */
  1119. {
  1120.   LVAL arg;
  1121.   FIXTYPE n, m, t, g;
  1122.  
  1123.   arg = xlgafixnum();
  1124.   n = getfixnum(arg);
  1125.   if (!moreargs())  {
  1126.      if (n < (FIXTYPE) 0) n = -n;
  1127.      return (cvfixnum(n));
  1128.   }
  1129.  
  1130.   while (moreargs())  {
  1131.      arg = xlgafixnum();
  1132.      m = getfixnum(arg);
  1133.      if ((n == (FIXTYPE) 0) || (m == (FIXTYPE) 0))
  1134.         return(cvfixnum(0));
  1135.  
  1136.      t = n * m;
  1137.      g = xget_gcd(n,m);
  1138.  
  1139.      n = (FIXTYPE) t / g;
  1140.   }
  1141.  
  1142.   if (n < (FIXTYPE) 0) n = -n;
  1143.  
  1144.   return (cvfixnum(n));
  1145. }
  1146.  
  1147. #ifndef RANDOM
  1148. LOCAL long rseed=1L;
  1149. #endif
  1150.  
  1151. /* unary - handle unary operations */
  1152. LOCAL LVAL NEAR unary(which)
  1153.         int which;
  1154. {
  1155.   FLOTYPE fval;
  1156.   FIXTYPE ival;
  1157. #ifdef RATIOS
  1158.   FIXTYPE numer, denom;
  1159. #endif
  1160.   Complex cval;
  1161.   LVAL arg, real, imag;
  1162.   int mode;
  1163.  
  1164.   /* get the argument */
  1165.   arg = xlgetarg();
  1166.   if (which == 'F' && moreargs()) { /*TAA MOD to eliminate compiler warnings*/
  1167.       if (floatp(*xlargv)) xlargc--;
  1168.       else xlbadtype(*xlargv);
  1169.   }
  1170.   xllastarg();
  1171.  
  1172.   /* check its type */
  1173.   if (fixp(arg)) {
  1174.     ival = getfixnum(arg);
  1175.     mode = IN;
  1176.   }
  1177.   else if (floatp(arg)) {
  1178.     fval = getflonum(arg);
  1179.     mode = FL;
  1180.   }
  1181. #ifdef RATIOS
  1182.   else if (ratiop(arg)) {
  1183.     numer = getnumer(arg);
  1184.     denom = getdenom(arg);
  1185.     mode = RT;
  1186.   }
  1187. #endif
  1188.   else if (complexp(arg)) {
  1189.     cval = makecomplex(arg);
  1190.     real = realpart(arg);
  1191.     imag = imagpart(arg);
  1192.     if (fixp(realpart(arg))) mode = CI;
  1193.     else mode = CF;
  1194.   }
  1195.   else xlerror("not a number", arg);
  1196.  
  1197.   switch (which) {
  1198.   case '~':
  1199.     if (mode == IN) return(cvfixnum((FIXTYPE) ~ival));
  1200.     else badiop();
  1201.     break;
  1202.   case 'A':
  1203.     switch (mode) {
  1204.     case IN: return(cvfixnum((FIXTYPE) (ival < 0   ? -ival : ival)));
  1205.     case FL: return(cvflonum((FLOTYPE) (fval < 0.0 ? -fval : fval)));
  1206.     case CI:
  1207.     case CF: return(cvflonum((FLOTYPE) modulus(cval)));
  1208. #ifdef RATIOS
  1209.     case RT: return(cvratio(numer<0?-numer:numer,denom));
  1210. #endif
  1211.     }
  1212.     break;
  1213.   case '+':
  1214.     switch (mode) {
  1215.     case IN: return(cvfixnum((FIXTYPE) ival + 1));
  1216.     case FL: return(cvflonum((FLOTYPE) fval + 1.0));
  1217.     case CI: return(newicomplex(getfixnum(real) + 1, getfixnum(imag)));
  1218.     case CF: return(newdcomplex(getflonum(real) + 1.0, getflonum(imag)));
  1219. #ifdef RATIOS
  1220.     case RT: return(cvratio(numer+denom, denom));
  1221. #endif
  1222.     }
  1223.     break;
  1224.   case '-':
  1225.     switch (mode) {
  1226.     case IN: return(cvfixnum((FIXTYPE) ival - 1));
  1227.     case FL: return(cvflonum((FLOTYPE) fval - 1.0));
  1228.     case CI: return(newicomplex(getfixnum(real) - 1, getfixnum(imag)));
  1229.     case CF: return(newdcomplex(getflonum(real) - 1.0, getflonum(imag)));
  1230. #ifdef RATIOS
  1231.     case RT: return(cvratio(numer-denom, denom));
  1232. #endif
  1233.     }
  1234.     break;
  1235.   case 'S':
  1236.     switch (mode) {
  1237.     case IN: return(cvflonum((FLOTYPE) sin((double) ival)));
  1238. #ifdef RATIOS
  1239.     case RT: fval = numer / (FLOTYPE) denom;
  1240. #endif
  1241.     case FL: return(cvflonum((FLOTYPE) sin((double) fval)));
  1242.     case CI:
  1243.     case CF: return(cvcomplex(csin(cval)));
  1244.     }
  1245.   case 'C':
  1246.     switch (mode) {
  1247.     case IN: return(cvflonum((FLOTYPE) cos((double) ival)));
  1248. #ifdef RATIOS
  1249.     case RT: fval = numer / (FLOTYPE) denom;
  1250. #endif
  1251.     case FL: return(cvflonum((FLOTYPE) cos((double) fval)));
  1252.     case CI:
  1253.     case CF: return(cvcomplex(ccos(cval)));
  1254.     }
  1255.   case 'T':
  1256.     switch (mode) {
  1257.     case IN: return(cvflonum((FLOTYPE) tan((double) ival)));
  1258. #ifdef RATIOS
  1259.     case RT: fval = numer / (FLOTYPE) denom;
  1260. #endif
  1261.     case FL: return(cvflonum((FLOTYPE) tan((double) fval)));
  1262.     case CI:
  1263.     case CF: return(cvcomplex(ctan(cval)));
  1264.     }
  1265.   case 'E':
  1266.     switch (mode) {
  1267.     case IN: return(cvflonum((FLOTYPE) exp((double) ival)));
  1268. #ifdef RATIOS
  1269.     case RT: fval = numer / (FLOTYPE) denom;
  1270. #endif
  1271.     case FL: return(cvflonum((FLOTYPE) exp((double) fval)));
  1272.     case CI:
  1273.     case CF: return(cvcomplex(cexp(cval)));
  1274.     }
  1275.     break;
  1276.   case 'R':
  1277.     switch (mode) {
  1278.     case IN:
  1279.       if (ival < 0) return(cvcomplex(csqrt(makecomplex(arg))));
  1280.           else return(cvflonum((FLOTYPE) sqrt((double) ival)));
  1281. #ifdef RATIOS
  1282.     case RT: fval = numer / (FLOTYPE) denom;
  1283. #endif
  1284.     case FL:
  1285.       if (fval < 0) return(cvcomplex(csqrt(makecomplex(arg))));
  1286.       else return(cvflonum((FLOTYPE) sqrt(fval)));
  1287.     case CI:
  1288.     case CF: return(cvcomplex(csqrt(cval)));
  1289.     }
  1290.     break;
  1291.   case 'F':
  1292.     switch (mode) {
  1293.     case IN: return (cvflonum((FLOTYPE) ival));
  1294. #ifdef RATIOS
  1295.     case RT: fval = numer / (FLOTYPE) denom;
  1296. #endif
  1297.     case FL: return (cvflonum((FLOTYPE) fval));
  1298.         default: badcop();
  1299.         }
  1300.         break;
  1301. #ifndef RANDOM
  1302.   case '?':
  1303.     switch (mode) {
  1304.     case IN: return (cvfixnum((FIXTYPE)(rseed=osrand(rseed)) % ival));
  1305.     case FL: badfop();
  1306.         default: badcop();
  1307.         }
  1308.         break;
  1309. #endif
  1310.   case 's':
  1311.     switch (mode) {
  1312. #ifdef RATIOS
  1313.     case RT: fval = numer / (FLOTYPE) denom; goto l1;
  1314. #endif
  1315.     case IN:
  1316.       fval = ival;
  1317.       /* drop through */
  1318.     case FL:
  1319. #ifdef RATIOS
  1320. l1:
  1321. #endif
  1322.       if (fval > 1.0 || fval < -1.0)
  1323.         return(cvcomplex(casin(makecomplex(arg))));
  1324.       else return(cvflonum((FLOTYPE) asin(fval)));
  1325.         case CI:
  1326.         case CF: return(cvcomplex(casin(cval)));
  1327.         }
  1328.         break;
  1329.   case 'c':
  1330.     switch (mode) {
  1331. #ifdef RATIOS
  1332.     case RT: fval = numer / (FLOTYPE) denom; goto l2;
  1333. #endif
  1334.     case IN:
  1335.       fval = ival;
  1336.       /* drop through */
  1337.     case FL:
  1338. #ifdef RATIOS
  1339. l2:
  1340. #endif
  1341.       if (fval > 1.0 || fval < -1.0)
  1342.         return(cvcomplex(cacos(makecomplex(arg))));
  1343.       else return(cvflonum((FLOTYPE) acos(fval)));
  1344.         case CI:
  1345.         case CF: return(cvcomplex(cacos(cval)));
  1346.         }
  1347.         break;
  1348.   case 'P':
  1349.     switch (mode) {
  1350.         case IN: return(cvflonum((FLOTYPE) (ival >= 0) ? 0.0 : PI));
  1351. #ifdef RATIOS
  1352.         case RT: fval = numer / (FLOTYPE) denom;
  1353. #endif
  1354.         case FL: return(cvflonum((FLOTYPE) (fval >= 0.0) ? 0.0 : PI));
  1355.         case CI:
  1356.         case CF: return(cvflonum((FLOTYPE) phase(cval)));
  1357.         }
  1358.         break;
  1359.   default: xlfail("unsupported operation");
  1360.   }
  1361.         return NIL; /* avoid compiler warning */
  1362. }
  1363.  
  1364. LOCAL LVAL NEAR unary2(which)   /* handle truncate, floor, ceiling, and round */
  1365.                             /* 1 or two arguments */
  1366.                             /* By Tom Almy */
  1367. int which;
  1368. {
  1369.     Number numer, denom;
  1370.     LVAL lval;
  1371.  
  1372.     lval = readnumber(&numer);
  1373.  
  1374.     if (moreargs()) {   /* two argument version */
  1375.         readnumber(&denom);
  1376.         xllastarg();
  1377.         matchmodes(&numer, &denom);
  1378.  
  1379.         switch (numer.mode) {
  1380.             case IN:
  1381.                 checkizero(denom.val);
  1382.                 numer.fval = numer.val / (double)denom.val;
  1383.                 break;
  1384.             case FL:
  1385.                 checkfzero(denom.fval);
  1386.                 numer.fval /= denom.fval;
  1387.                 break;
  1388. #ifdef RATIOS
  1389.             case RT:
  1390.                 checkizero(denom.num);
  1391.                 numer.fval = ((FLOTYPE)numer.num * denom.denom) /
  1392.                         ((FLOTYPE)numer.denom * denom.num);
  1393.                 break;
  1394. #endif
  1395.             default: badcop();
  1396.         }
  1397.     }
  1398.     else { /* single argument version */
  1399.         switch (numer.mode) {
  1400.             case IN: return lval;   /* no-operation */
  1401.             case FL: break;         /* continue */
  1402. #ifdef RATIOS
  1403.             case RT: numer.fval = numer.num / (FLOTYPE) numer.denom;
  1404.                 break;
  1405. #endif
  1406.             default: badcop();
  1407.         }
  1408.     }
  1409.     switch (which)  { /* now do it! */
  1410.         case 'I': modf(numer.fval,&numer.fval); break;
  1411.         case '_': numer.fval = floor(numer.fval); break;
  1412.         case '^': numer.fval = ceil(numer.fval); break;
  1413.         case 'r': numer.fval = floor(numer.fval + 0.5); break;
  1414.     }
  1415.     if (fabs(numer.fval) > (double)MAXFIX)
  1416.         return cvflonum((FLOTYPE)numer.fval);
  1417.     return cvfixnum((FIXTYPE)numer.fval);
  1418. }
  1419.  
  1420. /* unary functions */
  1421. LVAL xlognot() { return (unary('~')); } /* lognot */
  1422. LVAL xabs()    { return (unary('A')); } /* abs */
  1423. LVAL xadd1()   { return (unary('+')); } /* 1+ */
  1424. LVAL xsub1()   { return (unary('-')); } /* 1- */
  1425. LVAL xsin()    { return (unary('S')); } /* sin */
  1426. LVAL xcos()    { return (unary('C')); } /* cos */
  1427. LVAL xtan()    { return (unary('T')); } /* tan */
  1428. LVAL xexp()    { return (unary('E')); } /* exp */
  1429. LVAL xsqrt()   { return (unary('R')); } /* sqrt */
  1430. LVAL xfloat()  { return (unary('F')); } /* float */
  1431. #ifndef RANDOM
  1432. LVAL xrand()   { return (unary('?')); } /* random */
  1433. #endif
  1434. LVAL xasin()   { return (unary('s')); } /* asin */
  1435. LVAL xacos()   { return (unary('c')); } /* acos */
  1436. LVAL xphase()  { return (unary('P')); } /* phase */
  1437. LVAL xfix()    { return (unary2('I')); } /* truncate */
  1438. LVAL xfloor()  { return (unary2('_')); } /* floor */
  1439. LVAL xceil()   { return (unary2('^')); } /* ceiling */
  1440. LVAL xround()  { return (unary2('r')); } /* round */
  1441.  
  1442.  
  1443. /* predicate - handle a predicate function */
  1444. LOCAL LVAL NEAR predicate(fcn)
  1445.   int fcn;
  1446. {
  1447.   FLOTYPE fval;
  1448.   FIXTYPE ival;
  1449.   LVAL arg;
  1450.  
  1451.   /* get the argument */
  1452.   arg = xlgetarg();
  1453.   xllastarg();
  1454.  
  1455.   /* check the argument type */
  1456.   if (fixp(arg)) {
  1457.     ival = getfixnum(arg);
  1458.     switch (fcn) {
  1459.     case '-': ival = (ival < 0); break;
  1460.     case 'Z': ival = (ival == 0); break;
  1461.     case '+': ival = (ival > 0); break;
  1462.     case 'E': ival = ((ival & 1) == 0); break;
  1463.     case 'O': ival = ((ival & 1) != 0); break;
  1464.     default:  badiop();
  1465.     }
  1466.   }
  1467.   else if (floatp(arg)) {
  1468.     fval = getflonum(arg);
  1469.     switch (fcn) {
  1470.     case '-': ival = (fval < 0); break;
  1471.     case 'Z': ival = (fval == 0); break;
  1472.     case '+': ival = (fval > 0); break;
  1473.     default:  badfop();
  1474.     }
  1475.   }
  1476. #ifdef RATIOS
  1477.   else if (ratiop(arg)) {
  1478.     switch (fcn) {
  1479.     case '-': ival = (getnumer(arg) < 0); break;
  1480.     case 'Z': ival = FALSE; break;  /* can't be zero! */
  1481.     case '+': ival = (getnumer(arg) > 0); break;
  1482.     default:  badrop();
  1483.     }
  1484.   }
  1485. #endif
  1486.   else
  1487.     xlerror("bad argument type",arg);
  1488.  
  1489.   /* return the result value */
  1490.   return (ival ? true : NIL);
  1491. }
  1492.  
  1493. /* unary predicates */
  1494. LVAL xminusp() { return (predicate('-')); } /* minusp */
  1495. LVAL xzerop()  { return (predicate('Z')); } /* zerop */
  1496. LVAL xplusp()  { return (predicate('+')); } /* plusp */
  1497. LVAL xevenp()  { return (predicate('E')); } /* evenp */
  1498. LVAL xoddp()   { return (predicate('O')); } /* oddp */
  1499.  
  1500.  
  1501. /* compare - common compare function */
  1502. LOCAL LVAL NEAR compare(fcn)
  1503.   int fcn;
  1504. {
  1505.     FIXTYPE icmp,ival,iarg;
  1506.     FLOTYPE fcmp,fval,farg;
  1507.     LVAL arg;
  1508.     int mode;
  1509.  
  1510.     /* get the first argument */
  1511.     arg = xlgetarg();
  1512.  
  1513.     /* set the type of the first argument */
  1514.     if (fixp(arg)) {
  1515.         ival = getfixnum(arg);
  1516.         mode = 'I';
  1517.     }
  1518.     else if (floatp(arg)) {
  1519.         fval = getflonum(arg);
  1520.         mode = 'F';
  1521.     }
  1522. #ifdef RATIOS
  1523.     else if (ratiop(arg)) {
  1524.         fval = getnumer(arg) / (FLOTYPE) getdenom(arg);
  1525.         mode = 'F';
  1526.     }
  1527. #endif
  1528.     else
  1529.         xlerror("bad argument type",arg);
  1530.  
  1531.     /* handle each remaining argument */
  1532.     for (icmp = TRUE; icmp && moreargs(); ival = iarg, fval = farg) {
  1533.  
  1534.         /* get the next argument */
  1535.         arg = xlgetarg();
  1536.  
  1537.         /* check its type */
  1538.         if (fixp(arg)) {
  1539.             switch (mode) {
  1540.             case 'I':
  1541.                 iarg = getfixnum(arg);
  1542.                 break;
  1543.             case 'F':
  1544.                 farg = (FLOTYPE)getfixnum(arg);
  1545.                 break;
  1546.             }
  1547.         }
  1548.         else if (floatp(arg)) {
  1549.             switch (mode) {
  1550.             case 'I':
  1551.                 fval = (FLOTYPE)ival;
  1552.                 farg = getflonum(arg);
  1553.                 mode = 'F';
  1554.                 break;
  1555.             case 'F':
  1556.                 farg = getflonum(arg);
  1557.                 break;
  1558.             }
  1559.         }
  1560. #ifdef RATIOS
  1561.         else if (ratiop(arg)) {
  1562.             switch (mode) {
  1563.             case 'I':
  1564.                 fval = (FLOTYPE)ival;
  1565.             case 'F':
  1566.                 farg = getnumer(arg) / (FLOTYPE) getdenom(arg);
  1567.                 mode = 'F';
  1568.                 break;
  1569.             }
  1570.         }
  1571. #endif
  1572.         else
  1573.             xlerror("bad argument type",arg);
  1574.  
  1575.         /* compute result of the compare */
  1576.         switch (mode) {
  1577.         case 'I':
  1578.             icmp = ival - iarg;
  1579.             switch (fcn) {
  1580.             case '<':   icmp = (icmp < 0); break;
  1581.             case 'L':   icmp = (icmp <= 0); break;
  1582.             case '=':   icmp = (icmp == 0); break;
  1583.             case '#':   icmp = (icmp != 0); break;
  1584.             case 'G':   icmp = (icmp >= 0); break;
  1585.             case '>':   icmp = (icmp > 0); break;
  1586.             }
  1587.             break;
  1588.         case 'F':
  1589.             fcmp = fval - farg;
  1590.             switch (fcn) {
  1591.             case '<':   icmp = (fcmp < 0.0); break;
  1592.             case 'L':   icmp = (fcmp <= 0.0); break;
  1593.             case '=':   icmp = (fcmp == 0.0); break;
  1594.             case '#':   icmp = (fcmp != 0.0); break;
  1595.             case 'G':   icmp = (fcmp >= 0.0); break;
  1596.             case '>':   icmp = (fcmp > 0.0); break;
  1597.             }
  1598.             break;
  1599.         }
  1600.     }
  1601.  
  1602.     /* return the result */
  1603.     return (icmp ? true : NIL);
  1604. }
  1605.  
  1606. LOCAL LVAL NEAR ccompare(which)
  1607.         int which;
  1608. {
  1609.   Number val, arg;
  1610.   int icmp;
  1611.  
  1612.   switch (which) {
  1613.   case '=': icmp = TRUE;  break;
  1614.   case '#': icmp = FALSE; break;
  1615.   }
  1616.   /*larg =*/ readnumber(&val);
  1617.   while (moreargs()) {
  1618.     /*larg =*/ readnumber(&arg);
  1619.     matchmodes(&val, &arg);
  1620.     switch (which) {
  1621.     case '=':
  1622.       switch (val.mode) {
  1623.       case IN: icmp = icmp && val.val  == arg.val;  break;
  1624.       case FL: icmp = icmp && val.fval == arg.fval; break;
  1625.       case CI: icmp = icmp && val.crval==arg.crval && val.cival==arg.cival;
  1626.           break;
  1627.       case CF: icmp = icmp && val.cfrval==arg.cfrval && val.cfival==arg.cfival;
  1628.           break;
  1629. #ifdef RATIOS
  1630.       case RT: icmp = icmp && val.num == arg.num && val.denom == arg.denom;
  1631.           break;
  1632. #endif
  1633.       }
  1634.       break;
  1635.     case '#':
  1636.       switch (val.mode) {
  1637.       case IN: icmp = icmp || val.val  != arg.val;  break;
  1638.       case FL: icmp = icmp || val.fval != arg.fval; break;
  1639.       case CI: icmp = icmp || val.crval!=arg.crval || val.cival!=arg.cival;
  1640.           break;
  1641.       case CF: icmp = icmp || val.cfrval!=arg.cfrval || val.cfival!=arg.cfival;
  1642.           break;
  1643. #ifdef RATIOS
  1644.       case RT: icmp = icmp || val.num != arg.num || val.denom != arg.denom;
  1645.           break;
  1646. #endif
  1647.       }
  1648.       break;
  1649.     }
  1650.   }
  1651.   return((icmp) ? true : NIL);
  1652. }
  1653.  
  1654. /* comparison functions */
  1655. LVAL xlss() { return (compare('<'));  } /* < */
  1656. LVAL xleq() { return (compare('L'));  } /* <= */
  1657. LVAL xequ() { return (ccompare('=')); } /* = */
  1658. LVAL xneq() { return (ccompare('#')); } /* /= */
  1659. LVAL xgeq() { return (compare('G'));  } /* >= */
  1660. LVAL xgtr() { return (compare('>'));  } /* > */
  1661.  
  1662.  
  1663. /***********************************************************************/
  1664. /**                                                                   **/
  1665. /**                     Complex Number Functions                      **/
  1666. /**                                                                   **/
  1667. /***********************************************************************/
  1668.  
  1669. LVAL xcomplex() /* TAA rewrite so (complex 12.0) => #c(12.0 0.0) as required
  1670.                     by CL. */
  1671. {
  1672.     LVAL real, imag;
  1673.  
  1674.     real = xlgetarg();
  1675.  
  1676.     if (moreargs()) {
  1677.         imag = xlgetarg();
  1678.         xllastarg();
  1679.         return (newcomplex(real, imag));
  1680.     }
  1681.     if (fixp(real)) return (real);
  1682.     return (newcomplex(real,cvflonum((FLOTYPE)0.0)));
  1683. }
  1684.  
  1685. LVAL xconjugate()
  1686. {
  1687.   LVAL arg = xlgetarg();
  1688.  
  1689.   xllastarg();
  1690.   if (numberp(arg)) return(arg);
  1691.   if (!complexp(arg)) xlbadtype(arg);
  1692.   if (fixp(realpart(arg)) && fixp(imagpart(arg)))
  1693.     return(newicomplex(getfixnum(realpart(arg)), -getfixnum(imagpart(arg))));
  1694.   else
  1695.     return(newdcomplex(makefloat(realpart(arg)), -makefloat(imagpart(arg))));
  1696. }
  1697.  
  1698. LVAL xrealpart()
  1699. {
  1700.   LVAL arg = xlgetarg();
  1701.  
  1702.   xllastarg();
  1703.   if (fixp(arg) || floatp(arg)) return(arg);
  1704.   if (!complexp(arg)) xlbadtype(arg);
  1705.   return(realpart(arg));
  1706. }
  1707.  
  1708. LVAL ximagpart()
  1709. {
  1710.   LVAL arg = xlgetarg();
  1711.  
  1712.   xllastarg();
  1713.   if (fixp(arg)) return(cvfixnum((FIXTYPE) 0));
  1714.   if (floatp(arg)) return(cvflonum((FLOTYPE) 0.0));
  1715.   if (!complexp(arg)) xlbadtype(arg);
  1716.   return(imagpart(arg));
  1717. }
  1718.  
  1719. #endif
  1720.