"Fossies" - the Fresh Open Source Software Archive

Member "scm/scl.c" (29 Jun 2018, 79580 Bytes) of package /linux/privat/scm-5f3.zip:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "scl.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 5f2_vs_5f3.

    1 /* "scl.c" non-IEEE utility functions and non-integer arithmetic.
    2  * Copyright (C) 1990, 1991, 1992, 1993, 1994, 1995, 1997, 2005, 2006, 2013, 2014 Free Software Foundation, Inc.
    3  *
    4  * This program is free software: you can redistribute it and/or modify
    5  * it under the terms of the GNU Lesser General Public License as
    6  * published by the Free Software Foundation, either version 3 of the
    7  * License, or (at your option) any later version.
    8  *
    9  * This program is distributed in the hope that it will be useful, but
   10  * WITHOUT ANY WARRANTY; without even the implied warranty of
   11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   12  * Lesser General Public License for more details.
   13  *
   14  * You should have received a copy of the GNU Lesser General Public
   15  * License along with this program.  If not, see
   16  * <http://www.gnu.org/licenses/>.
   17  */
   18 
   19 /* Authors: Jerry D. Hedden and Aubrey Jaffer */
   20 
   21 #include "scm.h"
   22 
   23 #ifdef FLOATS
   24 # ifndef PLAN9
   25 #  include <math.h>
   26 # endif
   27 
   28 static SCM bigdblop P((int op, SCM b, double re, double im));
   29 static SCM inex_divintbig P((SCM a, SCM b));
   30 # ifndef BIGDIG
   31 static int apx_log10 P((double x));
   32 static double lpow10 P((double x, int n));
   33 # endif
   34 static sizet idbl2str P((double f, char *a));
   35 static sizet pdbl2str P((double f, char *a, sizet ch));
   36 static sizet iflo2str P((SCM flt, char *str));
   37 static void safe_add_1 P((double f, double *fsum));
   38 static long scm_twos_power P((SCM n));
   39 static double mantexp2dbl P((SCM manstr, int expo));
   40 static double pmantexp2dbl P((SCM bmant, int expo));
   41 static SCM ex2in P((SCM z));
   42 static SCM ilog P((unsigned long m, SCM b, SCM k, unsigned long *n));
   43 static double dpows5[23];
   44 
   45 static char s_makrect[] = "make-rectangular", s_makpolar[] = "make-polar",
   46         s_magnitude[] = "magnitude", s_angle[] = "angle",
   47         s_real_part[] = "real-part", s_imag_part[] = "imag-part",
   48         s_in2ex[] = "inexact->exact",s_ex2in[] = "exact->inexact";
   49 
   50 static char s_expt[] = "real-expt", s_atan2[] = "$atan2";
   51 #endif
   52 static char s_memv[] = "memv", s_assv[] = "assv";
   53 
   54 SCM sys_protects[NUM_PROTECTS];
   55 sizet num_protects = NUM_PROTECTS;
   56 
   57 char        s_inexactp[] = "inexact?";
   58 static char     s_zerop[] = "zero?", s_abs[] = "abs",
   59         s_positivep[] = "positive?", s_negativep[] = "negative?";
   60 static char     s_lessp[] = "<", s_grp[] = ">";
   61 static char     s_leqp[] = "<=", s_greqp[] = ">=";
   62 #define s_eqp (&s_leqp[1])
   63 static char     s_max[] = "max", s_min[] = "min";
   64 char        s_sum[] = "+", s_difference[] = "-", s_product[] = "*",
   65         s_divide[] = "/";
   66 static char     s_number2string[] = "number->string",
   67         s_str2number[] = "string->number";
   68 
   69 static char s_list_tail[] = "list-tail";
   70 static char s_str2list[] = "string->list";
   71 static char s_st_copy[] = "string-copy", s_st_fill[] = "string-fill!";
   72 static char s_vect2list[] = "vector->list", s_ve_fill[] = "vector-fill!";
   73 static char str_inf0[] = "inf.0", str_nan0[] = "nan.0", str_0[] = "0.";
   74 static char s_intexpt[] = "integer-expt", s_cintlog[] = "ceiling-integer-log";
   75 static char s_dfloat_parts[] = "double-float-parts";
   76 #define s_intlog (&s_cintlog[8])
   77 
   78 /*** NUMBERS -> STRINGS ***/
   79 #ifdef FLOATS
   80 static int dbl_mant_dig = 0;
   81 static double max_dbl_int;       /* Integers less than or equal to max_dbl_int
   82                     are representable exactly as doubles. */
   83 
   84 int inf2str(f, a)
   85      double f;
   86      char *a;
   87 {
   88   sizet ch = 0;
   89   if (f < 0.0) a[ch++] = '-';
   90   else if (f > 0.0) a[ch++] = '+';
   91   else {
   92     a[ch++] = '+';
   93     strcpy(&a[ch], str_nan0);
   94     return ch+sizeof(str_nan0)-1;
   95   }
   96   strcpy(&a[ch], str_inf0);
   97   return ch+sizeof(str_inf0)-1;
   98 }
   99 
  100 /* idbl2str() is the top level double-to-string conversion */
  101 static sizet idbl2str(f, a)
  102      double f;
  103      char *a;
  104 {
  105   sizet ch = 0;
  106   if (f==0.0) {strcpy(a, str_0); return sizeof(str_0)-1;}
  107   if (f==2*f) return inf2str(f, a);
  108   if (f < 0.0) {f = -f;a[ch++]='-';}
  109   else if (f > 0.0) ;
  110   else return inf2str(f, a);
  111   return pdbl2str(f, a, ch);
  112 }
  113 
  114 static double llog2 = .30102999566398114; /* log10(2) */
  115 
  116 /* There are also temporary strings used in number->string conversion. */
  117 void strrecy(str)
  118      SCM str;
  119 {
  120   if (IMP(str) || !STRINGP(str)) return;
  121   DEFER_INTS;
  122   must_free(CHARS(str), (sizet)LENGTH(str));
  123   CAR(str) = MAKINUM(1);
  124   CDR(str) = INUM0;
  125   ALLOW_INTS;
  126 }
  127 
  128 # ifdef BIGDIG
  129 
  130 /* The useful extent of bignums used in floating-point conversions is */
  131 /* limited.  Recycle them explicitly, rather than waiting for GC. */
  132 
  133 void bigrecy(bgnm)
  134      SCM bgnm;
  135 {
  136   if (IMP(bgnm) || !BIGP(bgnm)) return;
  137   DEFER_INTS;
  138   must_free(CHARS(bgnm), (sizet)NUMDIGS(bgnm)*sizeof(BIGDIG));
  139   CAR(bgnm) = INUM0;
  140   CDR(bgnm) = INUM0;
  141   ALLOW_INTS;
  142 }
  143 
  144 /* can convert to string accurately with bignums */
  145 /* f > 0 */
  146 /* DBL_MIN_EXP = -1021 */
  147 /* dbl_mant_dig = 53 */
  148 static sizet pdbl2str(f, a, ch)
  149      double f;
  150      char *a;
  151      sizet ch;
  152 {
  153   SCM mant, quo, num;
  154   sizet dp = ch;
  155   int e2, point, ndig = dbl_mant_dig;
  156   /* code from scm_dfloat_parts() */
  157   double dman = frexp(f, &e2);
  158 #  ifdef DBL_MIN_EXP
  159   if (e2 < DBL_MIN_EXP)
  160     ndig -= DBL_MIN_EXP - e2;
  161 #  endif
  162   e2 -= ndig;
  163   mant = dbl2big(ldexp(dman, ndig));
  164   point = ceil(e2*llog2);
  165   /* if (scm_verbose > 1) printf("mantissa = %g -> #x%s; e2 = %d -> %d; point = %d; ndig = %d -> %d\n", dman, CHARS(number2string(mant, MAKINUM(16))), e2+ndig, e2, point, dbl_mant_dig, ndig); */
  166   if (e2 >= 0) {
  167     /* try first with starved precision */
  168     {
  169       num = scm_ash(mant, MAKINUM(e2 - point));
  170       quo = scm_round_quotient(num, VELTS(pows5)[(long) point]);
  171       if (pmantexp2dbl(quo, point) != f) {
  172     int po2 = MAKINUM(1)==scm_logcount(mant);
  173     if (quo != num) { bigrecy(quo); quo = num; }
  174     num = scm_ash(quo, MAKINUM(1L));
  175     if (quo != num) bigrecy(quo);
  176     quo = scm_round_quotient(num, VELTS(pows5)[(long) --point]);
  177     if (po2 && pmantexp2dbl(quo,point) != f) {
  178       SCM quo1 = sum(quo, MAKINUM(1L));
  179       if (quo != quo1) { bigrecy(quo); quo = quo1; }
  180       if (pmantexp2dbl(quo,point) != f) {
  181         if (quo != num) { bigrecy(quo); quo = num; };
  182         num = scm_ash(quo,MAKINUM(1L));
  183         quo = scm_round_quotient(num, VELTS(pows5)[(long) --point]);
  184       }
  185     }
  186       }
  187       if (num != quo) bigrecy(num);
  188     }
  189   } else {   /* e2 <= 0 */
  190     /* try first with starved precision */
  191     {
  192       SCM den = scm_ash(MAKINUM(1L), MAKINUM(point - e2));
  193       num = product(mant, VELTS(pows5)[- (long) point]);
  194       quo = scm_round_quotient(num, den);
  195       if (pmantexp2dbl(quo, point) != f) {
  196     int po2 = MAKINUM(1L)==scm_logcount(mant);
  197     if (quo != num) { bigrecy(quo); quo = num; }
  198     point--;
  199     num = product(quo, MAKINUM(10));
  200     /* if (quo != num) bigrecy(quo); */
  201     quo = scm_round_quotient(num, den);
  202     if (po2 && pmantexp2dbl(quo,point) != f) {
  203       SCM quo1 = sum(quo, MAKINUM(1L));
  204       if (quo != quo1) { bigrecy(quo); quo = quo1; }
  205       if (pmantexp2dbl(quo,point) != f) {
  206         if (quo != num) { bigrecy(quo); quo = num; }
  207         point--;
  208         num = product(quo, MAKINUM(10));
  209         if (quo != num) bigrecy(quo);
  210         quo = scm_round_quotient(num, den);
  211       }
  212     }
  213       } else if (quo != num) bigrecy(num);
  214       bigrecy(den);
  215     }
  216   }
  217   if (mant != quo) bigrecy(mant);
  218   a[ch++] = '.';
  219   /* if (sizeof(UBIGLONG)>=sizeof(double)) /\* Is ulong larger than mantissa? *\/ */
  220   /*   ch += iulong2str(num2ulong(quo, (char *)ARG1, s_number2string), 10, &a[ch]); */
  221   /* else */ {
  222     SCM str = number2string(quo, MAKINUM(10));
  223     int len = LENGTH(str), i = 0;
  224     bigrecy(quo);
  225     point += len - 1;
  226     while (i < len) a[ch++] = CHARS(str)[i++];
  227     strrecy(str);
  228   }
  229   a[dp] = a[dp+1]; a[++dp] = '.';
  230 #  ifdef ENGNOT
  231   while ((dp+1 < ch) && (point+9999)%3) { a[dp] = a[dp+1]; a[++dp] = '.'; point--; }
  232 #  endif    /* ENGNOT */
  233   while ('0'==a[--ch]); ch++;
  234   if (point != 0) {
  235     a[ch++] = 'e';
  236     return ch + ilong2str(point, 10, &a[ch]);
  237   } else return ch;
  238 }
  239 # else  /* ~BIGDIG */
  240 
  241 /* DBL2STR_FUZZ is a somewhat arbitrary guard against
  242    round off error in scaling f and fprec. */
  243 #  define DBL2STR_FUZZ 0.9
  244 int dblprec;
  245 /* static double dbl_eps; */
  246 double dbl_prec(x)
  247      double x;
  248 {
  249   int expt;
  250   double frac = frexp(x, &expt);
  251 #  ifdef DBL_MIN_EXP
  252   if (0.0==x || expt < DBL_MIN_EXP) /* gradual underflow */
  253     return ldexp(1.0, - dbl_mant_dig) * ldexp(1.0, DBL_MIN_EXP);
  254 #  endif
  255   if (1.0==frac) return ldexp(1.0, expt - dbl_mant_dig + 1);
  256   return ldexp(1.0, expt - dbl_mant_dig);
  257 }
  258 
  259 static int apx_log10(x)
  260      double x;
  261 {
  262   int expt;
  263   frexp(x, &expt);
  264   expt -= 1;
  265   if (expt >= 0)
  266     return (int)(expt * llog2);
  267   return -((int)(-expt * llog2));
  268 }
  269 
  270 static double p10[] = {1.0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7};
  271 static double lpow10(x, n)
  272      double x;
  273      int n;
  274 {
  275   if (n >= 0) {
  276     while (n > 7) {
  277       x *= 1e8;
  278       n -= 8;
  279     }
  280     return x*p10[n];
  281   }
  282   while (n < -7) {
  283     x /= 1e8;
  284     n += 8;
  285   }
  286   return x/p10[-n];
  287 }
  288 
  289 /* f is finite and positive */
  290 static sizet pdbl2str(f, a, ch)
  291      double f;
  292      char *a;
  293      sizet ch;
  294 {
  295   double fprec = dbl_prec(f);
  296   int efmt, dpt, d, i, exp = apx_log10(f);
  297   f = lpow10(f, -exp);
  298   fprec = lpow10(fprec, -exp);
  299 #  ifdef DBL_MIN_10_EXP /* Prevent random unnormalized values, as from
  300                make-uniform-vector, from causing infinite loops,
  301                but try to print gradually underflowing numbers. */
  302   while (f < 1.0) {
  303     f *= 10.0;
  304     fprec *= 10.0;
  305     if (exp-- < DBL_MIN_10_EXP - DBL_DIG - 1) return inf2str(f, a);
  306   }
  307   while (f > 10.0) {
  308     f /= 10.0;
  309     fprec /= 10.0;
  310     if (exp++ > DBL_MAX_10_EXP) return inf2str(f, a);
  311   }
  312 #  else
  313   while (f < 1.0) {f *= 10.0; fprec *= 10.0; exp--;}
  314   while (f > 10.0) {f /= 10.0; fprec /= 10.0; exp++;}
  315 #  endif
  316   fprec *= 0.5;
  317   if (f+fprec >= 10.0) {f = 1.0; exp++;}
  318  /* zero: */
  319 #  ifdef ENGNOT
  320   dpt = (exp+9999)%3;
  321   exp -= dpt++;
  322   efmt = 1;
  323 #  else
  324   efmt = (exp < -3) || (exp > dblprec+2);
  325   if (!efmt)
  326     if (exp < 0) {
  327       a[ch++] = '0';
  328       a[ch++] = '.';
  329       dpt = exp;
  330       while (++dpt)  a[ch++] = '0';
  331     } else
  332       dpt = exp+1;
  333   else
  334     dpt = 1;
  335 #  endif
  336 
  337   for (i = 30; i--;) {
  338     /* printf("  f = %.20g, fprec = %.20g, i = %d\n", f, fprec, i); */
  339     d = f;
  340     f -= d;
  341     a[ch++] = d+'0';
  342     if (f < fprec && f < DBL2STR_FUZZ*fprec) break;
  343     if ((f + fprec) >= 1.0 && (f + DBL2STR_FUZZ*fprec) >= 1.0) {
  344       a[ch-1]++;
  345       break;
  346     }
  347     f *= 10.0;
  348     fprec *= 10.0;
  349     if (!(--dpt))  a[ch++] = '.';
  350   }
  351 
  352   if (dpt > 0)
  353 #  ifndef ENGNOT
  354     if ((dpt > 4) && (exp > 6)) {
  355       d = (a[0]=='-'?2:1);
  356       for (i = ch++; i > d; i--)
  357     a[i] = a[i-1];
  358       a[d] = '.';
  359       efmt = 1;
  360     } else
  361 #  endif
  362       {
  363     while (--dpt)  a[ch++] = '0';
  364     a[ch++] = '.';
  365       }
  366   if (a[ch-1]=='.')  a[ch++]='0'; /* trailing zero */
  367   if (efmt && exp) {
  368     a[ch++] = 'e';
  369     if (exp < 0) {
  370       exp = -exp;
  371       a[ch++] = '-';
  372     }
  373     for (i = 10; i <= exp; i *= 10);
  374     for (i /= 10; i; i /= 10) {
  375       a[ch++] = exp/i + '0';
  376       exp %= i;
  377     }
  378   }
  379   return ch;
  380 }
  381 # endif /* ~BIGDIG */
  382 
  383 static sizet iflo2str(flt, str)
  384      SCM flt;
  385      char *str;
  386 {
  387   sizet i;
  388 # ifdef SINGLES
  389   if (SINGP(flt)) i = idbl2str(FLO(flt), str);
  390   else
  391 # endif
  392     i = idbl2str(REAL(flt), str);
  393   if (scm_narn==flt) return i;
  394   if (CPLXP(flt)) {
  395     if (!(0 > IMAG(flt))) str[i++] = '+';
  396     i += idbl2str(IMAG(flt), &str[i]);
  397     str[i++] = 'i';
  398   }
  399   return i;
  400 }
  401 #endif              /* FLOATS */
  402 
  403 sizet iulong2str(num, rad, p)
  404      unsigned long num;
  405      int rad;
  406      char *p;
  407 {
  408   sizet j;
  409   register int i = 1, d;
  410   register unsigned long n = num;
  411   for (n /= rad;n > 0;n /= rad) i++;
  412   j = i;
  413   n = num;
  414   while (i--) {
  415     d = n % rad;
  416     n /= rad;
  417     p[i] = d + ((d < 10) ? '0' : 'a' - 10);
  418   }
  419   return j;
  420 }
  421 sizet ilong2str(num, rad, p)
  422      long num;
  423      int rad;
  424      char *p;
  425 {
  426   if ((num < 0) && !(rad < 0)) {
  427     *p++ = '-';
  428     return 1 + iulong2str((unsigned long) -num, rad, p);
  429   }
  430   return iulong2str((unsigned long) num, rad < 0 ? -rad : rad, p);
  431 }
  432 #ifdef BIGDIG
  433 static SCM big2str(b, radix)
  434      SCM b;
  435      register unsigned int radix;
  436 {
  437   SCM t = copybig(b, 0);    /* sign of temp doesn't matter */
  438   register BIGDIG *ds = BDIGITS(t);
  439   sizet i = NUMDIGS(t);
  440   sizet j = radix==16 ? (BITSPERDIG*i)/4+2
  441     : radix >= 10 ? (BITSPERDIG*i*241L)/800+2
  442     : (BITSPERDIG*i)+2;
  443   sizet k = 0;
  444   sizet radct = 0;
  445   sizet ch; /* jeh */
  446   BIGDIG radpow = 1, radmod = 0;
  447   SCM ss = makstr((long)j);
  448   char *s = CHARS(ss), c;
  449   scm_protect_temp(&t);
  450   while ((long) radpow * radix < BIGRAD) {
  451     radpow *= radix;
  452     radct++;
  453   }
  454   s[0] = tc16_bigneg==TYP16(b) ? '-' : '+';
  455   while ((i || radmod) && j) {
  456     if (k==0) {
  457       radmod = (BIGDIG)divbigdig(ds, i, radpow);
  458       k = radct;
  459       if (!ds[i-1]) i--;
  460     }
  461     c = radmod % radix; radmod /= radix; k--;
  462     s[--j] = c < 10 ? c + '0' : c + 'a' - 10;
  463   }
  464   ch = s[0]=='-' ? 1 : 0; /* jeh */
  465   if (ch < j) { /* jeh */
  466     for (i = j;j < LENGTH(ss);j++) s[ch+j-i] = s[j]; /* jeh */
  467     resizuve(ss, (SCM)MAKINUM(ch+LENGTH(ss)-i)); /* jeh */
  468   }
  469   bigrecy(t);
  470   return ss;
  471 }
  472 #endif
  473 SCM number2string(x, radix)
  474      SCM x, radix;
  475 {
  476   if (UNBNDP(radix)) radix=MAKINUM(10L);
  477   else ASRTER(INUMP(radix), radix, ARG2, s_number2string);
  478 #ifdef FLOATS
  479   if (NINUMP(x)) {
  480     char num_buf[FLOBUFLEN];
  481 # ifdef BIGDIG
  482     ASRTGO(NIMP(x), badx);
  483     if (BIGP(x)) return big2str(x, (unsigned int)INUM(radix));
  484 #  ifndef RECKLESS
  485     if (!(INEXP(x)))
  486       badx: wta(x, (char *)ARG1, s_number2string);
  487 #  endif
  488 # else
  489     ASRTER(NIMP(x) && INEXP(x), x, ARG1, s_number2string);
  490 # endif
  491     return makfromstr(num_buf, iflo2str(x, num_buf));
  492   }
  493 #else
  494 # ifdef BIGDIG
  495   if (NINUMP(x)) {
  496     ASRTER(NIMP(x) && BIGP(x), x, ARG1, s_number2string);
  497     return big2str(x, (unsigned int)INUM(radix));
  498   }
  499 # else
  500   ASRTER(INUMP(x), x, ARG1, s_number2string);
  501 # endif
  502 #endif
  503   {
  504     char num_buf[INTBUFLEN];
  505     return makfromstr(num_buf, ilong2str(INUM(x), (int)INUM(radix), num_buf));
  506   }
  507 }
  508 /* These print routines are stubbed here so that repl.c doesn't need
  509    FLOATS or BIGDIGs conditionals */
  510 int floprint(sexp, port, writing)
  511      SCM sexp;
  512      SCM port;
  513      int writing;
  514 {
  515 #ifdef FLOATS
  516   if (!errjmp_bad) {
  517     char num_buf[FLOBUFLEN];
  518     lfwrite(num_buf, (sizet)sizeof(char), iflo2str(sexp, num_buf), port);
  519     return !0;
  520   } else
  521 #endif
  522   scm_ipruk("float", sexp, port);
  523   return !0;
  524 }
  525 int bigprint(exp, port, writing)
  526      SCM exp;
  527      SCM port;
  528      int writing;
  529 {
  530 #ifdef BIGDIG
  531   if (!errjmp_bad) {
  532     exp = big2str(exp, (unsigned int)10);
  533     lfwrite(CHARS(exp), (sizet)sizeof(char), (sizet)LENGTH(exp), port);
  534     return !0;
  535   } else
  536 #endif
  537   scm_ipruk("bignum", exp, port);
  538   return !0;
  539 }
  540 /*** END nums->strs ***/
  541 
  542 /*** STRINGS -> NUMBERS ***/
  543 #ifdef BIGDIG
  544 SCM istr2int(str, len, radix)
  545      char *str;
  546      long len;
  547      register int radix;
  548 {
  549   register sizet k, blen = 1;
  550   sizet i = 0, j;
  551   int c;
  552   SCM res;
  553   register BIGDIG *ds;
  554   register UBIGLONG t2;
  555 
  556   if (0 >= len) return BOOL_F;  /* zero length */
  557   /* Estimate number of digits; will trim during finish */
  558   if (10==radix) j = 1+(84*len)/(BITSPERDIG*25);
  559   else j = (8 < radix) ? 1+(4*len)/BITSPERDIG : 1+(3*len)/BITSPERDIG;
  560   switch (str[0]) {     /* leading sign */
  561   case '-':
  562   case '+': if (++i==len) return BOOL_F; /* bad if lone `+' or `-' */
  563   }
  564   res = mkbig(j, '-'==str[0]);
  565   ds = BDIGITS(res);
  566   /* clear allocated digits */
  567   for (k = j;k--;) ds[k] = 0;
  568   do {
  569     switch (c = str[i++]) {
  570     case DIGITS:
  571       c = c - '0';
  572       goto accumulate;
  573     case 'A': case 'B': case 'C': case 'D': case 'E': case 'F':
  574       c = c-'A'+10;
  575       goto accumulate;
  576     case 'a': case 'b': case 'c': case 'd': case 'e': case 'f':
  577       c = c-'a'+10;
  578     accumulate:
  579       if (c >= radix) return BOOL_F; /* bad digit for radix */
  580       k = 0;
  581       t2 = c;
  582     moretodo:
  583       /* Add t2 into temp bignum */
  584       while (k < blen) {
  585     t2 += ((UBIGLONG)ds[k])*radix;
  586     ds[k++] = BIGLO(t2);
  587     t2 = BIGDN(t2);
  588       }
  589       ASRTER(blen <= j, (SCM)MAKINUM(blen), OVFLOW, "bignum");
  590       if (t2) {blen++; goto moretodo;}
  591       break;
  592     default:
  593       return BOOL_F;        /* not a digit */
  594     }
  595   } while (i < len);
  596   if (blen * BITSPERDIG/CHAR_BIT <= sizeof(SCM))
  597     if (INUMP(res = big2inum(res, blen))) return res;
  598   if (j==blen) return res;
  599   return adjbig(res, blen);
  600 }
  601 #else
  602 SCM istr2int(str, len, radix)
  603      register char *str;
  604      long len;
  605      register int radix;
  606 {
  607   register long n = 0, ln;
  608   register int c;
  609   register int i = 0;
  610   int lead_neg = 0;
  611   if (0 >= len) return BOOL_F;  /* zero length */
  612   switch (*str) {       /* leading sign */
  613   case '-': lead_neg = 1;
  614   case '+': if (++i==len) return BOOL_F; /* bad if lone `+' or `-' */
  615   }
  616 
  617   do {
  618     switch (c = str[i++]) {
  619     case DIGITS:
  620       c = c - '0';
  621       goto accumulate;
  622     case 'A': case 'B': case 'C': case 'D': case 'E': case 'F':
  623       c = c-'A'+10;
  624       goto accumulate;
  625     case 'a': case 'b': case 'c': case 'd': case 'e': case 'f':
  626       c = c-'a'+10;
  627     accumulate:
  628       if (c >= radix) return BOOL_F; /* bad digit for radix */
  629       ln = n;
  630       n = n * radix - c;
  631       if (n > ln
  632 # ifdef hpux
  633       || (-n > -MOST_NEGATIVE_FIXNUM) /* workaround for HP700 cc bug */
  634 # endif
  635       ) goto ovfl;
  636       break;
  637     default:
  638       return BOOL_F;        /* not a digit */
  639     }
  640   } while (i < len);
  641   if (lead_neg) {
  642     if (n < MOST_NEGATIVE_FIXNUM) goto ovfl;
  643   }
  644   else {
  645     if (n < -MOST_POSITIVE_FIXNUM) goto ovfl;
  646     n = -n;
  647   }
  648   return MAKINUM(n);
  649  ovfl:              /* overflow scheme integer */
  650   return BOOL_F;
  651 }
  652 #endif
  653 
  654 #ifdef FLOATS
  655 # ifdef BIGDIG
  656 
  657 /* [Jaffer] */
  658 /* In Clinger's "How to Read Floating-Point Numbers Accurately" */
  659 /* the key issue is that successive rounding does not have the */
  660 /* same effect as one rounding operation.  With bignums it is */
  661 /* a simple matter to accumulate digits without rounding. */
  662 /* The approach here is to compute the binary value without rounding, */
  663 /* then round explicitly when scaling the bigint to fit into the */
  664 /* mantissa. */
  665 
  666 /* manstr is the mantissa with the decimal point removed and point
  667    (the exponent) adjusted appropriately */
  668 double mantexp2dbl(manstr, point)
  669      SCM manstr;
  670      int point;
  671 {
  672   SCM bmant = istr2int(CHARS(manstr), LENGTH(manstr), 10L);
  673   double val = pmantexp2dbl(bmant, point);
  674   bigrecy(bmant);
  675   return val;
  676 }
  677 double pmantexp2dbl(bmant, point)
  678      SCM bmant;
  679      int point;
  680 {
  681   double ans;
  682   if (BOOL_F == bmant) return REAL(scm_narn);
  683   if (point >= 0) {
  684     if (point < 23 && INUM(scm_intlength(bmant)) <= dbl_mant_dig)
  685       return ldexp(num2dbl(bmant,ARG1,s_str2number) * dpows5[point], point);
  686     {
  687       SCM quo, num = product(bmant, VELTS(pows5)[(long) point]);
  688       int bex = INUM(scm_intlength(num)) - dbl_mant_dig;
  689       if (bex > 0) {
  690     SCM den = scm_ash(MAKINUM(1L), MAKINUM(bex));
  691     quo = scm_round_quotient(num, den);
  692     bigrecy(den);
  693       }
  694       else
  695     quo = scm_ash(num, MAKINUM(- bex));
  696       /* quo may not be a bignum */
  697       if (INUMP(quo)) ans = ldexp((double)(INUM(quo)), bex + point);
  698       else {
  699     sizet i = NUMDIGS(quo);
  700     sizet j = (dbl_mant_dig + BITSPERDIG - 1)/BITSPERDIG;
  701     BIGDIG *digits = BDIGITS(quo);
  702     if (j < i) j = i - j;
  703     else j = 0;
  704     ans = 0.0;
  705     while (i-- > j) ans = digits[i] + ldexp(ans, BITSPERDIG);
  706     bex += j * BITSPERDIG;
  707     ans = ldexp(ans, bex + point);
  708       }
  709       if (num != quo) bigrecy(quo);
  710       if ((num != bmant) && (bmant != MAKINUM(1L))) bigrecy(num);
  711       return ans;
  712     }
  713   }
  714   if (-point < 23 && INUM(scm_intlength(bmant)) <= dbl_mant_dig)
  715     return ldexp(num2dbl(bmant,ARG1,s_str2number) / dpows5[-point], point);
  716   {
  717     int maxpow = LENGTH(pows5) - 1;
  718     SCM num, quo, scl = (-point <= maxpow) ?
  719       VELTS(pows5)[(long) -point] :
  720       product(VELTS(pows5)[(long)maxpow], VELTS(pows5)[(long)-point-maxpow]);
  721     int mantlen = dbl_mant_dig;
  722     int bex =           /* bex < 0 */
  723       INUM(scm_intlength(bmant)) - INUM(scm_intlength(scl)) - mantlen;
  724     int tmp = bex + point + 1021 + mantlen;
  725     if (tmp < 0) {
  726       bex -= tmp + 1;
  727       mantlen += tmp;
  728     }
  729     num = scm_ash(bmant, MAKINUM(-bex));
  730     quo = scm_round_quotient(num, scl);
  731     if (INUM(scm_intlength(quo)) > mantlen) {
  732       bex++;            /* too many bits of quotient */
  733       quo = scm_round_quotient(num, scm_ash(scl, MAKINUM(1L)));
  734     }
  735     if (-point > maxpow) bigrecy(scl);
  736     if (num != quo) bigrecy(num);
  737     ans = ldexp(int2dbl(quo), bex + point);
  738     bigrecy(quo);
  739     return ans;
  740   }
  741 }
  742 
  743 # else /* def BIGDIG */
  744 
  745 double mantexp2dbl(manstr, point)
  746      SCM manstr;
  747      int point;
  748 {
  749   register int c, i = 0;
  750   double res = 0.0;
  751   char *str = CHARS(manstr);
  752   int len = LENGTH(manstr);
  753   do {              /* check initial digits */
  754     switch (c = str[i]) {
  755     case DIGITS:
  756       c = c - '0';
  757       res = res * 10 + c;
  758       break;
  759     case 'D': case 'E': case 'F':
  760     case 'd': case 'e': case 'f':
  761     default:
  762       goto out1;
  763     }
  764   } while (++i < len);
  765  out1:
  766   if (point >= 0)
  767     while (point--)  res *= 10.0;
  768   else
  769     while (point++) {
  770 #  ifdef _UNICOS
  771       res *= 0.1;
  772 #  else
  773       res /= 10.0;
  774 #  endif
  775     }
  776   return res;
  777 }
  778 
  779 # endif /* def BIGDIG */
  780 
  781 SCM istr2flo(str, len, radix)
  782      register char *str;
  783      register long len;
  784      register long radix;
  785 {
  786   register int c, i = 0, j = 0;
  787   int lead_sgn = 0;
  788   double res = 0.0, tmp = 0.0;
  789   int flg = 0;
  790   int point = 0;
  791   int shrtp = 0;
  792   SCM second, manstr;
  793   char *mant;
  794 
  795   if (i >= len) return BOOL_F;  /* zero length */
  796 
  797   switch (*str) {       /* leading sign */
  798   case '-': lead_sgn = -1; i++; break;
  799   case '+': lead_sgn = 1; i++; break;
  800   }
  801   if (i==len) return BOOL_F;    /* bad if lone `+' or `-' */
  802 
  803   if (6==len && ('+'==str[0] || '-'==str[0])) {
  804     if (0==strcasecmp(str_inf0, &str[1]))
  805       return makdbl(1./0. * ('+'==str[0] ? 1 : -1), 0.0);
  806     else if (0==strcasecmp(str_nan0, &str[1]))
  807       return scm_narn;
  808   }
  809   if (str[i]=='i' || str[i]=='I') { /* handle `+i' and `-i'   */
  810     if (lead_sgn==0) return BOOL_F; /* must have leading sign */
  811     if (++i < len) return BOOL_F; /* `i' not last character */
  812     return makdbl(0.0, (double)lead_sgn);
  813   }
  814 
  815   manstr = makstr(len);
  816   mant = CHARS(manstr);
  817 
  818   do {              /* check initial digits */
  819     switch (c = str[i]) {
  820     case DIGITS:
  821       c = c - '0';
  822       goto accum1;
  823     case 'D': case 'E': case 'F':
  824       if (radix==10) goto out1; /* must be exponent */
  825     case 'A': case 'B': case 'C':
  826       c = c-'A'+10;
  827       goto accum1;
  828     case 'd': case 'e': case 'f':
  829       if (radix==10) goto out1;
  830     case 'a': case 'b': case 'c':
  831       c = c-'a'+10;
  832     accum1:
  833       if (c >= radix) return BOOL_F; /* bad digit for radix */
  834       mant[j++] = str[i];
  835       res = res * radix + c;
  836       flg = 1;          /* res is valid */
  837       break;
  838     default:
  839       goto out1;
  840     }
  841   } while (++i < len);
  842  out1:
  843 
  844   /* if true, then we did see a digit above, and res is valid */
  845   if (i==len) goto done;
  846 
  847   /* By here, must have seen a digit,
  848      or must have next char be a `.' with radix==10 */
  849   if ((!flg) && (!(str[i]=='.' && radix==10))) return BOOL_F;
  850 
  851   while (str[i]=='#') {     /* optional sharps */
  852     res *= radix;
  853     mant[j++] = '0';
  854     if (++i==len) goto done;
  855   }
  856 
  857   if (str[i]=='/' && i+1 < len) {
  858     while (++i < len) {
  859       switch (c = str[i]) {
  860       case DIGITS:
  861     c = c - '0';
  862     goto accum2;
  863       case 'A': case 'B': case 'C': case 'D': case 'E': case 'F':
  864     c = c-'A'+10;
  865     goto accum2;
  866       case 'a': case 'b': case 'c': case 'd': case 'e': case 'f':
  867     c = c-'a'+10;
  868       accum2:
  869     if (c >= radix) return BOOL_F;
  870     tmp = tmp * radix + c;
  871     break;
  872       default:
  873     goto out2;
  874       }
  875     }
  876   out2:
  877     /*     if (tmp==0.0) return BOOL_F; /\* `slash zero' not allowed *\/ */
  878     if (i < len)
  879       while (str[i]=='#') { /* optional sharps */
  880     tmp *= radix;
  881     mant[j++] = '0';
  882     if (++i==len) break;
  883       }
  884     res /= tmp;
  885     goto done;
  886   }
  887 
  888   if (str[i]=='.') {        /* decimal point notation */
  889     if (radix != 10) return BOOL_F; /* must be radix 10 */
  890     while (++i < len) {
  891       switch (c = str[i]) {
  892       case DIGITS:
  893     point--;
  894     res = res*10.0 + c-'0';
  895     flg = 1;
  896     mant[j++] = str[i];
  897     break;
  898       default:
  899     goto out3;
  900       }
  901     }
  902   out3:
  903     if (!flg) return BOOL_F;    /* no digits before or after decimal point */
  904     if (i==len) goto adjust;
  905     while (str[i]=='#') {   /* ignore remaining sharps */
  906       if (++i==len) goto adjust;
  907     }
  908   }
  909 
  910   switch (str[i]) {     /* exponent */
  911   case 'f': case 'F':
  912   case 's': case 'S':
  913     shrtp = !0;
  914   case 'd': case 'D':
  915   case 'e': case 'E':
  916   case 'l': case 'L': {
  917       int expsgn = 1, expon = 0;
  918       if (radix != 10) return BOOL_F; /* only in radix 10 */
  919       if (++i==len) return BOOL_F; /* bad exponent */
  920       switch (str[i]) {
  921       case '-':  expsgn=(-1);
  922       case '+':  if (++i==len) return BOOL_F; /* bad exponent */
  923       }
  924       if (str[i] < '0' || str[i] > '9') return BOOL_F; /* bad exponent */
  925       do {
  926     switch (c = str[i]) {
  927     case DIGITS:
  928       expon = expon*10 + c-'0';
  929       /*    if (expon > MAXEXP) */
  930       /*      if (1==expsgn || expon > (MAXEXP + dblprec + 1)) */
  931       /*        return BOOL_F; /\* exponent too large *\/ */
  932       break;
  933     default:
  934       goto out4;
  935     }
  936       } while (++i < len);
  937     out4:
  938       point += expsgn*expon;
  939     }
  940   }
  941 
  942  adjust:
  943   mant[j] = 0;
  944   manstr = resizuve(manstr, MAKINUM(j));
  945   if (radix == 10) res = mantexp2dbl(manstr, point);
  946 
  947  done:
  948   /* at this point, we have a legitimate floating point result */
  949   if (lead_sgn==-1)  res = -res;
  950   if (i==len) return shrtp ? makflo(res) : makdbl(res, 0.0);
  951 
  952   if (str[i]=='i' || str[i]=='I') { /* pure imaginary number  */
  953     if (lead_sgn==0) return BOOL_F; /* must have leading sign */
  954     if (++i < len) return BOOL_F; /* `i' not last character */
  955     return makdbl(0.0, res);
  956   }
  957 
  958   switch (str[i++]) {
  959   case '-':  lead_sgn = -1; break;
  960   case '+':  lead_sgn = 1;  break;
  961   case '@': {           /* polar input for complex number */
  962     /* get a `real' for angle */
  963     second = istr2flo(&str[i], (long)(len-i), radix);
  964     if (IMP(second)) return BOOL_F;
  965     if (!(INEXP(second))) return BOOL_F; /* not `real' */
  966     if (CPLXP(second))    return BOOL_F; /* not `real' */
  967     tmp = REALPART(second);
  968     return makdbl(res*cos(tmp), res*sin(tmp));
  969   }
  970   default: return BOOL_F;
  971   }
  972 
  973   /* at this point, last char must be `i' */
  974   if (str[len-1] != 'i' && str[len-1] != 'I') return BOOL_F;
  975   /* handles `x+i' and `x-i' */
  976   if (i==(len-1))  return makdbl(res, (double)lead_sgn);
  977   /* get a `ureal' for complex part */
  978   second = istr2flo(&str[i], (long)((len-i)-1), radix);
  979   if (IMP(second)) return BOOL_F;
  980   if (!(INEXP(second))) return BOOL_F; /* not `ureal' */
  981   if (CPLXP(second))    return BOOL_F; /* not `ureal' */
  982   tmp = REALPART(second);
  983   if (tmp < 0.0)    return BOOL_F; /* not `ureal' */
  984   return makdbl(res, (lead_sgn*tmp));
  985 }
  986 #endif              /* FLOATS */
  987 
  988 
  989 SCM istring2number(str, len, radix)
  990      char *str;
  991      long len;
  992      long radix;
  993 {
  994   int i = 0;
  995   char ex = 0;
  996   char ex_p = 0, rx_p = 0;  /* Only allow 1 exactness and 1 radix prefix */
  997   SCM res;
  998   if (len==1)
  999     if (*str=='+' || *str=='-') /* Catches lone `+' and `-' for speed */
 1000       return BOOL_F;
 1001 
 1002   while ((len-i) >= 2  &&  str[i]=='#' && ++i)
 1003     switch (str[i++]) {
 1004     case 'b': case 'B':  if (rx_p++) return BOOL_F; radix = 2;  break;
 1005     case 'o': case 'O':  if (rx_p++) return BOOL_F; radix = 8;  break;
 1006     case 'd': case 'D':  if (rx_p++) return BOOL_F; radix = 10; break;
 1007     case 'x': case 'X':  if (rx_p++) return BOOL_F; radix = 16; break;
 1008     case 'i': case 'I':  if (ex_p++) return BOOL_F; ex = 2;     break;
 1009     case 'e': case 'E':  if (ex_p++) return BOOL_F; ex = 1;     break;
 1010     default:  return BOOL_F;
 1011     }
 1012 
 1013   switch (ex) {
 1014   case 1:
 1015     return istr2int(&str[i], len-i, radix);
 1016   case 0:
 1017     res = istr2int(&str[i], len-i, radix);
 1018     if (NFALSEP(res)) return res;
 1019 #ifdef FLOATS
 1020   case 2: return istr2flo(&str[i], len-i, radix);
 1021 #endif
 1022   }
 1023   return BOOL_F;
 1024 }
 1025 
 1026 
 1027 SCM string2number(str, radix)
 1028      SCM str, radix;
 1029 {
 1030   if (UNBNDP(radix)) radix=MAKINUM(10L);
 1031   else ASRTER(INUMP(radix), radix, ARG2, s_str2number);
 1032   ASRTER(NIMP(str) && STRINGP(str), str, ARG1, s_str2number);
 1033   return istring2number(CHARS(str), LENGTH(str), INUM(radix));
 1034 }
 1035 /*** END strs->nums ***/
 1036 
 1037 #ifdef FLOATS
 1038 SCM makdbl (x, y)
 1039      double x, y;
 1040 {
 1041   SCM z;
 1042   if ((y==0.0) && (x==0.0)) return flo0;
 1043 # ifndef _MSC_VER
 1044 #  ifndef SINGLESONLY
 1045   if ((y != y) || (x != x) || (y==(2 * y) && (y != 0.0))) return scm_narn;
 1046   if ((x==(2 * x)) && (x != 0.0)) y = 0.0;
 1047 #  endif
 1048 # endif
 1049   DEFER_INTS;
 1050   if (y==0.0) {
 1051 # ifdef SINGLES
 1052     float fx = x;        /* David Yeh <theyeh@uclink.berkeley.edu>
 1053                 changed this so that MSVC works */
 1054 #  ifndef SINGLESONLY
 1055     if ((-FLTMAX < x) && (x < FLTMAX) && ( (double)fx == x) )
 1056 #  endif
 1057       {
 1058     NEWCELL(z);
 1059     CAR(z) = tc_flo;
 1060     FLO(z) = x;
 1061     ALLOW_INTS;
 1062     return z;
 1063       }
 1064 # endif             /* def SINGLES */
 1065     z = must_malloc_cell(1L*sizeof(double), (SCM)tc_dblr, "real");
 1066   }
 1067   else {
 1068     z = must_malloc_cell(2L*sizeof(double), (SCM)tc_dblc, "complex");
 1069     IMAG(z) = y;
 1070   }
 1071   REAL(z) = x;
 1072   ALLOW_INTS;
 1073   return z;
 1074 }
 1075 #endif              /* FLOATS */
 1076 
 1077 #ifndef INUMS_ONLY
 1078 SCM eqv(x, y)
 1079      SCM x, y;
 1080 {
 1081   if (x==y) return BOOL_T;
 1082   if (IMP(x)) return BOOL_F;
 1083   if (IMP(y)) return BOOL_F;
 1084   /* this ensures that types and length are the same. */
 1085   if (CAR(x) != CAR(y)) return BOOL_F;
 1086   if (NUMP(x)) {
 1087 # ifdef BIGDIG
 1088     if (BIGP(x)) return (0==bigcomp(x, y)) ? BOOL_T : BOOL_F;
 1089 # endif
 1090 # ifdef FLOATS
 1091     return floequal(x, y);
 1092 # endif
 1093   }
 1094   return BOOL_F;
 1095 }
 1096 SCM memv(x, lst)            /* m.borza  12.2.91 */
 1097 SCM x, lst;
 1098 {
 1099   for (;NIMP(lst);lst = CDR(lst)) {
 1100     ASRTGO(CONSP(lst), badlst);
 1101     if (NFALSEP(eqv(CAR(lst), x))) return lst;
 1102   }
 1103 # ifndef RECKLESS
 1104   if (!(NULLP(lst)))
 1105     badlst: wta(lst, (char *)ARG2, s_memv);
 1106 # endif
 1107   return BOOL_F;
 1108 }
 1109 SCM assv(x, alist)      /* m.borza  12.2.91 */
 1110 SCM x, alist;
 1111 {
 1112   SCM tmp;
 1113   for (;NIMP(alist);alist = CDR(alist)) {
 1114     ASRTGO(CONSP(alist), badlst);
 1115     tmp = CAR(alist);
 1116     ASRTGO(NIMP(tmp) && CONSP(tmp), badlst);
 1117     if (NFALSEP(eqv(CAR(tmp), x))) return tmp;
 1118   }
 1119 # ifndef RECKLESS
 1120   if (!(NULLP(alist)))
 1121     badlst: wta(alist, (char *)ARG2, s_assv);
 1122 # endif
 1123   return BOOL_F;
 1124 }
 1125 #endif
 1126 
 1127 SCM list_tail(lst, k)
 1128      SCM lst, k;
 1129 {
 1130   register long i;
 1131   ASRTER(INUMP(k), k, ARG2, s_list_tail);
 1132   i = INUM(k);
 1133   while (i-- > 0) {
 1134     ASRTER(NIMP(lst) && CONSP(lst), lst, ARG1, s_list_tail);
 1135     lst = CDR(lst);
 1136   }
 1137   return lst;
 1138 }
 1139 
 1140 SCM string2list(str)
 1141      SCM str;
 1142 {
 1143   long i;
 1144   SCM res = EOL;
 1145   unsigned char *src;
 1146   ASRTER(NIMP(str) && STRINGP(str), str, ARG1, s_str2list);
 1147   src = UCHARS(str);
 1148   for (i = LENGTH(str)-1;i >= 0;i--) res = cons((SCM)MAKICHR(src[i]), res);
 1149   return res;
 1150 }
 1151 SCM string_copy(str)
 1152      SCM str;
 1153 {
 1154   ASRTER(NIMP(str) && STRINGP(str), str, ARG1, s_st_copy);
 1155   return makfromstr(CHARS(str), (sizet)LENGTH(str));
 1156 }
 1157 SCM string_fill(str, chr)
 1158      SCM str, chr;
 1159 {
 1160   register char *dst, c;
 1161   register long k;
 1162   ASRTER(NIMP(str) && STRINGP(str), str, ARG1, s_st_fill);
 1163   ASRTER(ICHRP(chr), chr, ARG2, s_st_fill);
 1164   c = ICHR(chr);
 1165   dst = CHARS(str);
 1166   for (k = LENGTH(str)-1;k >= 0;k--) dst[k] = c;
 1167   return UNSPECIFIED;
 1168 }
 1169 SCM vector2list(v)
 1170      SCM v;
 1171 {
 1172   SCM res = EOL;
 1173   long i;
 1174   SCM *data;
 1175   ASRTER(NIMP(v) && VECTORP(v), v, ARG1, s_vect2list);
 1176   data = VELTS(v);
 1177   for (i = LENGTH(v)-1;i >= 0;i--) res = cons(data[i], res);
 1178   return res;
 1179 }
 1180 SCM vector_fill(v, fill)
 1181      SCM v, fill;
 1182 {
 1183   register long i;
 1184   register SCM *data;
 1185   ASRTER(NIMP(v) && VECTORP(v), v, ARG1, s_ve_fill);
 1186   data = VELTS(v);
 1187   for (i = LENGTH(v)-1;i >= 0;i--) data[i] = fill;
 1188   return UNSPECIFIED;
 1189 }
 1190 static SCM vector_equal(x, y)
 1191      SCM x, y;
 1192 {
 1193   long i;
 1194   for (i = LENGTH(x)-1;i >= 0;i--)
 1195     if (FALSEP(equal(VELTS(x)[i], VELTS(y)[i]))) return BOOL_F;
 1196   return BOOL_T;
 1197 }
 1198 #ifdef BIGDIG
 1199 SCM bigequal(x, y)
 1200      SCM x, y;
 1201 {
 1202   if (0==bigcomp(x, y)) return BOOL_T;
 1203   return BOOL_F;
 1204 }
 1205 #endif
 1206 #ifdef FLOATS
 1207 SCM floequal(x, y)
 1208      SCM x, y;
 1209 {
 1210   if ((REALPART(x) != REALPART(y))) return BOOL_F;
 1211   if (CPLXP(x))
 1212     return (CPLXP(y) && (IMAG(x)==IMAG(y))) ? BOOL_T : BOOL_F;
 1213   return CPLXP(y) ? BOOL_F : BOOL_T;
 1214 }
 1215 #endif
 1216 SCM equal(x, y)
 1217      SCM x, y;
 1218 {
 1219   CHECK_STACK;
 1220  tailrecurse: POLL;
 1221   if (x==y) return BOOL_T;
 1222   if (IMP(x)) return BOOL_F;
 1223   if (IMP(y)) return BOOL_F;
 1224   if (CONSP(x) && CONSP(y)) {
 1225     if (FALSEP(equal(CAR(x), CAR(y)))) return BOOL_F;
 1226     x = CDR(x);
 1227     y = CDR(y);
 1228     goto tailrecurse;
 1229   }
 1230   /* this ensures that types and length are the same. */
 1231   if (CAR(x) != CAR(y)) return BOOL_F;
 1232   switch (TYP7(x)) {
 1233   default: return BOOL_F;
 1234   case tc7_string: return st_equal(x, y);
 1235   case tc7_vector: return vector_equal(x, y);
 1236   case tc7_smob: {
 1237     int i = SMOBNUM(x);
 1238     if (!(i < numsmob)) return BOOL_F;
 1239     if (smobs[i].equalp) return (smobs[i].equalp)(x, y);
 1240     else return BOOL_F;
 1241   }
 1242   case tc7_VfixN8: case tc7_VfixZ8: case tc7_VfixN16: case tc7_VfixZ16:
 1243   case tc7_VfixN32: case tc7_VfixZ32: case tc7_VfixN64: case tc7_VfixZ64:
 1244   case tc7_VfloR32: case tc7_VfloC32: case tc7_VfloC64: case tc7_VfloR64:
 1245   case tc7_Vbool: {
 1246     SCM (*pred)() = smobs[0x0ff & (tc16_array>>8)].equalp;
 1247     if (pred) return (*pred)(x, y);
 1248     else return BOOL_F;
 1249   }
 1250   }
 1251 }
 1252 
 1253 SCM numberp(obj)
 1254      SCM obj;
 1255 {
 1256   if (INUMP(obj)) return BOOL_T;
 1257 #ifdef FLOATS
 1258   if (NIMP(obj) && NUMP(obj)) return BOOL_T;
 1259 #else
 1260 # ifdef BIGDIG
 1261   if (NIMP(obj) && NUMP(obj)) return BOOL_T;
 1262 # endif
 1263 #endif
 1264   return BOOL_F;
 1265 }
 1266 #ifdef FLOATS
 1267 SCM scm_complex_p(obj)
 1268      SCM obj;
 1269 {
 1270   if (obj==scm_narn) return BOOL_F;
 1271   return numberp(obj);
 1272 }
 1273 
 1274 # ifdef BIGDIG
 1275 static char twostab[] = {4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0};
 1276 static long scm_twos_power(n)
 1277      SCM n;
 1278 {
 1279   long d, c = 0;
 1280   int d4;
 1281 #  ifdef BIGDIG
 1282   if (NINUMP(n)) {
 1283     BIGDIG *ds;
 1284     int i = 0;
 1285     ds = BDIGITS(n);
 1286     while (0==(d = ds[i++])) c += BITSPERDIG;
 1287     goto out;
 1288   }
 1289 #  endif
 1290   d = INUM(n);
 1291   if (0==d) return 0;
 1292  out:
 1293   do {
 1294     d4 = 15 & d;
 1295     c += twostab[d4];
 1296     d >>= 4;
 1297   } while (0==d4);
 1298   return c;
 1299 }
 1300 
 1301 int scm_bigdblcomp(b, d)
 1302      SCM b;
 1303      double d;
 1304 {
 1305   int dlen = 99;
 1306   sizet blen = 0;
 1307   int dneg = d < 0 ? 1 : 0;
 1308   int bneg = BIGSIGN(b) ? 1 : 0;
 1309   if (bneg < dneg) return -1;
 1310   if (bneg > dneg) return 1;
 1311   if (!(d==2*d && d != 0.0)) {
 1312     frexp(d, &dlen);
 1313     blen = INUM(scm_intlength(b));
 1314   }
 1315   if (blen > dlen) return dneg ? 1 : -1;
 1316   if (blen < dlen) return dneg ? -1 : 1;
 1317   if ((blen <= dbl_mant_dig) || (blen - scm_twos_power(b)) <= dbl_mant_dig) {
 1318     double bd = int2dbl(b);
 1319     if (bd > d) return -1;
 1320     if (bd < d) return 1;
 1321     return 0;
 1322   }
 1323   return bigcomp(b, dbl2big(d));
 1324 }
 1325 # endif
 1326 SCM realp(x)
 1327      SCM x;
 1328 {
 1329   if (INUMP(x)) return BOOL_T;
 1330   if (IMP(x)) return BOOL_F;
 1331   if (REALP(x)) return BOOL_T;
 1332 # ifdef BIGDIG
 1333   if (BIGP(x)) return BOOL_T;
 1334 # endif
 1335   return BOOL_F;
 1336 }
 1337 SCM scm_rationalp(x)
 1338      SCM x;
 1339 {
 1340   if (INUMP(x)) return BOOL_T;
 1341   if (IMP(x)) return BOOL_F;
 1342   if (REALP(x)) {
 1343     float y = REALPART(x);
 1344     if (y==2*y && y != 0.0) return BOOL_F;
 1345     return BOOL_T;
 1346   }
 1347 # ifdef BIGDIG
 1348   if (BIGP(x)) return BOOL_T;
 1349 # endif
 1350   return BOOL_F;
 1351 }
 1352 SCM intp(x)
 1353      SCM x;
 1354 {
 1355   double r;
 1356   if (INUMP(x)) return BOOL_T;
 1357   if (IMP(x)) return BOOL_F;
 1358 # ifdef BIGDIG
 1359   if (BIGP(x)) return BOOL_T;
 1360 # endif
 1361   if (!INEXP(x)) return BOOL_F;
 1362   if (CPLXP(x)) return BOOL_F;
 1363   r = REALPART(x);
 1364   if (r != floor(r)) return BOOL_F;
 1365   if (r==2*r && r != 0.0) return BOOL_F;
 1366   return BOOL_T;
 1367 }
 1368 #endif              /* FLOATS */
 1369 
 1370 SCM inexactp(x)
 1371      SCM x;
 1372 {
 1373 #ifdef FLOATS
 1374   if (NIMP(x) && INEXP(x)) return BOOL_T;
 1375 #endif
 1376   return BOOL_F;
 1377 }
 1378 
 1379 SCM eqp(x, y)
 1380      SCM x, y;
 1381 {
 1382 #ifdef FLOATS
 1383   SCM t;
 1384   if (NINUMP(x)) {
 1385 # ifdef BIGDIG
 1386 #  ifndef RECKLESS
 1387     if (!(NIMP(x)))
 1388       badx: wta(x, (char *)ARG1, s_eqp);
 1389 #  endif
 1390     if (BIGP(x)) {
 1391       if (INUMP(y)) return BOOL_F;
 1392       ASRTGO(NIMP(y), bady);
 1393       if (BIGP(y)) return (0==bigcomp(x, y)) ? BOOL_T : BOOL_F;
 1394       ASRTGO(INEXP(y), bady);
 1395     bigreal:
 1396       return (REALP(y) && (0==scm_bigdblcomp(x, REALPART(y)))) ?
 1397     BOOL_T : BOOL_F;
 1398     }
 1399     ASRTGO(INEXP(x), badx);
 1400 # else
 1401     ASRTER(NIMP(x) && INEXP(x), x, ARG1, s_eqp);
 1402 # endif
 1403     if (INUMP(y)) {t = x; x = y; y = t; goto realint;}
 1404 # ifdef BIGDIG
 1405     ASRTGO(NIMP(y), bady);
 1406     if (BIGP(y)) {t = x; x = y; y = t; goto bigreal;}
 1407     ASRTGO(INEXP(y), bady);
 1408 # else
 1409     ASRTGO(NIMP(y) && INEXP(y), bady);
 1410 # endif
 1411     if (x==y) return BOOL_T;
 1412     return floequal(x, y);
 1413   }
 1414   if (NINUMP(y)) {
 1415 # ifdef BIGDIG
 1416     ASRTGO(NIMP(y), bady);
 1417     if (BIGP(y)) return BOOL_F;
 1418 #  ifndef RECKLESS
 1419     if (!(INEXP(y)))
 1420       bady: wta(y, (char *)ARG2, s_eqp);
 1421 #  endif
 1422 # else
 1423 #  ifndef RECKLESS
 1424     if (!(NIMP(y) && INEXP(y)))
 1425       bady: wta(y, (char *)ARG2, s_eqp);
 1426 #  endif
 1427 # endif
 1428   realint:
 1429     return (REALP(y) && (((double)INUM(x))==REALPART(y))) ? BOOL_T : BOOL_F;
 1430   }
 1431 #else
 1432 # ifdef BIGDIG
 1433   if (NINUMP(x)) {
 1434     ASRTER(NIMP(x) && BIGP(x), x, ARG1, s_eqp);
 1435     if (INUMP(y)) return BOOL_F;
 1436     ASRTGO(NIMP(y) && BIGP(y), bady);
 1437     return (0==bigcomp(x, y)) ? BOOL_T : BOOL_F;
 1438   }
 1439   if (NINUMP(y)) {
 1440 #  ifndef RECKLESS
 1441     if (!(NIMP(y) && BIGP(y)))
 1442       bady: wta(y, (char *)ARG2, s_eqp);
 1443 #  endif
 1444     return BOOL_F;
 1445   }
 1446 # else
 1447   ASRTER(INUMP(x), x, ARG1, s_eqp);
 1448   ASRTER(INUMP(y), y, ARG2, s_eqp);
 1449 # endif
 1450 #endif
 1451   return ((long)x==(long)y) ? BOOL_T : BOOL_F;
 1452 }
 1453 SCM lessp(x, y)
 1454      SCM x, y;
 1455 {
 1456 #ifdef FLOATS
 1457   if (NINUMP(x)) {
 1458 # ifdef BIGDIG
 1459 #  ifndef RECKLESS
 1460     if (!(NIMP(x)))
 1461       badx: wta(x, (char *)ARG1, s_lessp);
 1462 #  endif
 1463     if (BIGP(x)) {
 1464       if (INUMP(y)) return BIGSIGN(x) ? BOOL_T : BOOL_F;
 1465       ASRTGO(NIMP(y), bady);
 1466       if (BIGP(y)) return (1==bigcomp(x, y)) ? BOOL_T : BOOL_F;
 1467       ASRTGO(REALP(y), bady);
 1468       return (1==scm_bigdblcomp(x, REALPART(y))) ? BOOL_T : BOOL_F;
 1469     }
 1470     ASRTGO(REALP(x), badx);
 1471 # else
 1472     ASRTER(NIMP(x) && REALP(x), x, ARG1, s_lessp);
 1473 # endif
 1474     if (INUMP(y)) return (REALPART(x) < ((double)INUM(y))) ? BOOL_T : BOOL_F;
 1475 # ifdef BIGDIG
 1476     ASRTGO(NIMP(y), bady);
 1477     if (BIGP(y)) return (-1==scm_bigdblcomp(y, REALPART(x))) ? BOOL_T : BOOL_F;
 1478     ASRTGO(REALP(y), bady);
 1479 # else
 1480     ASRTGO(NIMP(y) && REALP(y), bady);
 1481 # endif
 1482     return (REALPART(x) < REALPART(y)) ? BOOL_T : BOOL_F;
 1483   }
 1484   if (NINUMP(y)) {
 1485 # ifdef BIGDIG
 1486     ASRTGO(NIMP(y), bady);
 1487     if (BIGP(y)) return BIGSIGN(y) ? BOOL_F : BOOL_T;
 1488 #  ifndef RECKLESS
 1489     if (!(REALP(y)))
 1490       bady: wta(y, (char *)ARG2, s_lessp);
 1491 #  endif
 1492 # else
 1493 #  ifndef RECKLESS
 1494     if (!(NIMP(y) && REALP(y)))
 1495       bady: wta(y, (char *)ARG2, s_lessp);
 1496 #  endif
 1497 # endif
 1498     return (((double)INUM(x)) < REALPART(y)) ? BOOL_T : BOOL_F;
 1499   }
 1500 #else
 1501 # ifdef BIGDIG
 1502   if (NINUMP(x)) {
 1503     ASRTER(NIMP(x) && BIGP(x), x, ARG1, s_lessp);
 1504     if (INUMP(y)) return BIGSIGN(x) ? BOOL_T : BOOL_F;
 1505     ASRTGO(NIMP(y) && BIGP(y), bady);
 1506     return (1==bigcomp(x, y)) ? BOOL_T : BOOL_F;
 1507   }
 1508   if (NINUMP(y)) {
 1509 #  ifndef RECKLESS
 1510     if (!(NIMP(y) && BIGP(y)))
 1511       bady: wta(y, (char *)ARG2, s_lessp);
 1512 #  endif
 1513     return BIGSIGN(y) ? BOOL_F : BOOL_T;
 1514   }
 1515 # else
 1516   ASRTER(INUMP(x), x, ARG1, s_lessp);
 1517   ASRTER(INUMP(y), y, ARG2, s_lessp);
 1518 # endif
 1519 #endif
 1520   return ((long)x < (long)y) ? BOOL_T : BOOL_F;
 1521 }
 1522 SCM greaterp(x, y)
 1523      SCM x, y;
 1524 {
 1525   return lessp(y, x);
 1526 }
 1527 SCM leqp(x, y)
 1528      SCM x, y;
 1529 {
 1530   return BOOL_NOT(lessp(y, x));
 1531 }
 1532 SCM greqp(x, y)
 1533      SCM x, y;
 1534 {
 1535   return BOOL_NOT(lessp(x, y));
 1536 }
 1537 SCM zerop(z)
 1538      SCM z;
 1539 {
 1540 #ifdef FLOATS
 1541   if (NINUMP(z)) {
 1542 # ifdef BIGDIG
 1543     ASRTGO(NIMP(z), badz);
 1544     if (BIGP(z)) return BOOL_F;
 1545 #  ifndef RECKLESS
 1546     if (!(INEXP(z)))
 1547       badz: wta(z, (char *)ARG1, s_zerop);
 1548 #  endif
 1549 # else
 1550     ASRTER(NIMP(z) && INEXP(z), z, ARG1, s_zerop);
 1551 # endif
 1552     return (z==flo0) ? BOOL_T : BOOL_F;
 1553   }
 1554 #else
 1555 # ifdef BIGDIG
 1556   if (NINUMP(z)) {
 1557     ASRTER(NIMP(z) && BIGP(z), z, ARG1, s_zerop);
 1558     return BOOL_F;
 1559   }
 1560 # else
 1561   ASRTER(INUMP(z), z, ARG1, s_zerop);
 1562 # endif
 1563 #endif
 1564   return (z==INUM0) ? BOOL_T: BOOL_F;
 1565 }
 1566 SCM positivep(x)
 1567      SCM x;
 1568 {
 1569 #ifdef FLOATS
 1570   if (NINUMP(x)) {
 1571 # ifdef BIGDIG
 1572     ASRTGO(NIMP(x), badx);
 1573     if (BIGP(x)) return TYP16(x)==tc16_bigpos ? BOOL_T : BOOL_F;
 1574 #  ifndef RECKLESS
 1575     if (!(REALP(x)))
 1576       badx: wta(x, (char *)ARG1, s_positivep);
 1577 #  endif
 1578 # else
 1579     ASRTER(NIMP(x) && REALP(x), x, ARG1, s_positivep);
 1580 # endif
 1581     return (REALPART(x) > 0.0) ? BOOL_T : BOOL_F;
 1582   }
 1583 #else
 1584 # ifdef BIGDIG
 1585   if (NINUMP(x)) {
 1586     ASRTER(NIMP(x) && BIGP(x), x, ARG1, s_positivep);
 1587     return TYP16(x)==tc16_bigpos ? BOOL_T : BOOL_F;
 1588   }
 1589 # else
 1590   ASRTER(INUMP(x), x, ARG1, s_positivep);
 1591 # endif
 1592 #endif
 1593   return (x > INUM0) ? BOOL_T : BOOL_F;
 1594 }
 1595 SCM negativep(x)
 1596      SCM x;
 1597 {
 1598 #ifdef FLOATS
 1599   if (NINUMP(x)) {
 1600 # ifdef BIGDIG
 1601     ASRTGO(NIMP(x), badx);
 1602     if (BIGP(x)) return TYP16(x)==tc16_bigpos ? BOOL_F : BOOL_T;
 1603 #  ifndef RECKLESS
 1604     if (!(REALP(x)))
 1605       badx: wta(x, (char *)ARG1, s_negativep);
 1606 #  endif
 1607 # else
 1608     ASRTER(NIMP(x) && REALP(x), x, ARG1, s_negativep);
 1609 # endif
 1610     return (REALPART(x) < 0.0) ? BOOL_T : BOOL_F;
 1611   }
 1612 #else
 1613 # ifdef BIGDIG
 1614   if (NINUMP(x)) {
 1615     ASRTER(NIMP(x) && BIGP(x), x, ARG1, s_negativep);
 1616     return (TYP16(x)==tc16_bigneg) ? BOOL_T : BOOL_F;
 1617   }
 1618 # else
 1619   ASRTER(INUMP(x), x, ARG1, s_negativep);
 1620 # endif
 1621 #endif
 1622   return (x < INUM0) ? BOOL_T : BOOL_F;
 1623 }
 1624 
 1625 static char s_exactprob[] = "not representable as inexact";
 1626 SCM scm_max(x, y)
 1627      SCM x, y;
 1628 {
 1629 #ifdef FLOATS
 1630   SCM t;
 1631   double z;
 1632 #endif
 1633   if (UNBNDP(y)) {
 1634 #ifndef RECKLESS
 1635     if (!(NUMBERP(x)))
 1636       badx: wta(x, (char *)ARG1, s_max);
 1637 #endif
 1638     return x;
 1639   }
 1640 #ifdef FLOATS
 1641   if (NINUMP(x)) {
 1642 # ifdef BIGDIG
 1643     ASRTGO(NIMP(x), badx);
 1644     if (BIGP(x)) {
 1645       if (INUMP(y)) return BIGSIGN(x) ? y : x;
 1646       ASRTGO(NIMP(y), bady);
 1647       if (BIGP(y)) return (1==bigcomp(x, y)) ? y : x;
 1648       ASRTGO(REALP(y), bady);
 1649     big_dbl:
 1650       if (-1 != scm_bigdblcomp(x, REALPART(y))) return y;
 1651       z = int2dbl(x);
 1652       ASRTER(0==scm_bigdblcomp(x, z), x, s_exactprob, s_max);
 1653       return makdbl(z, 0.0);
 1654     }
 1655     ASRTGO(REALP(x), badx);
 1656 # else
 1657     ASRTER(NIMP(x) && REALP(x), x, ARG1, s_max);
 1658 # endif
 1659     if (INUMP(y)) return (REALPART(x) < (z = INUM(y))) ? makdbl(z, 0.0) : x;
 1660 # ifdef BIGDIG
 1661     ASRTGO(NIMP(y), bady);
 1662     if (BIGP(y)) {
 1663       t = y; y = x; x = t; goto big_dbl;
 1664     }
 1665     ASRTGO(REALP(y), bady);
 1666 # else
 1667     ASRTGO(NIMP(y) && REALP(y), bady);
 1668 # endif
 1669     return (REALPART(x) < REALPART(y)) ? y : x;
 1670   }
 1671   if (NINUMP(y)) {
 1672 # ifdef BIGDIG
 1673     ASRTGO(NIMP(y), bady);
 1674     if (BIGP(y)) return BIGSIGN(y) ? x : y;
 1675 #  ifndef RECKLESS
 1676     if (!(REALP(y)))
 1677       bady: wta(y, (char *)ARG2, s_max);
 1678 #  endif
 1679 # else
 1680 #  ifndef RECKLESS
 1681     if (!(NIMP(y) && REALP(y)))
 1682       bady: wta(y, (char *)ARG2, s_max);
 1683 #  endif
 1684 # endif
 1685     return ((z = INUM(x)) < REALPART(y)) ? y : makdbl(z, 0.0);
 1686   }
 1687 #else
 1688 # ifdef BIGDIG
 1689   if (NINUMP(x)) {
 1690     ASRTER(NIMP(x) && BIGP(x), x, ARG1, s_max);
 1691     if (INUMP(y)) return BIGSIGN(x) ? y : x;
 1692     ASRTGO(NIMP(y) && BIGP(y), bady);
 1693     return (1==bigcomp(x, y)) ? y : x;
 1694   }
 1695   if (NINUMP(y)) {
 1696 #  ifndef RECKLESS
 1697     if (!(NIMP(y) && BIGP(y)))
 1698       bady: wta(y, (char *)ARG2, s_max);
 1699 #  endif
 1700     return BIGSIGN(y) ? x : y;
 1701   }
 1702 # else
 1703   ASRTER(INUMP(x), x, ARG1, s_max);
 1704   ASRTER(INUMP(y), y, ARG2, s_max);
 1705 # endif
 1706 #endif
 1707   return ((long)x < (long)y) ? y : x;
 1708 }
 1709 
 1710 SCM scm_min(x, y)
 1711      SCM x, y;
 1712 {
 1713 #ifdef FLOATS
 1714   SCM t;
 1715   double z;
 1716 #endif
 1717   if (UNBNDP(y)) {
 1718 #ifndef RECKLESS
 1719     if (!(NUMBERP(x)))
 1720       badx: wta(x, (char *)ARG1, s_min);
 1721 #endif
 1722     return x;
 1723   }
 1724 #ifdef FLOATS
 1725   if (NINUMP(x)) {
 1726 # ifdef BIGDIG
 1727     ASRTGO(NIMP(x), badx);
 1728     if (BIGP(x)) {
 1729       if (INUMP(y)) return BIGSIGN(x) ? x : y;
 1730       ASRTGO(NIMP(y), bady);
 1731       if (BIGP(y)) return (-1==bigcomp(x, y)) ? y : x;
 1732       ASRTGO(REALP(y), bady);
 1733     big_dbl:
 1734       if (1 != scm_bigdblcomp(x, REALPART(y))) return y;
 1735       z = int2dbl(x);
 1736       ASRTER(0==scm_bigdblcomp(x, z), x, s_exactprob, s_min);
 1737       return makdbl(z, 0.0);
 1738     }
 1739     ASRTGO(REALP(x), badx);
 1740 # else
 1741     ASRTER(NIMP(x) && REALP(x), x, ARG1, s_min);
 1742 # endif
 1743     if (INUMP(y)) return (REALPART(x) > (z = INUM(y))) ? makdbl(z, 0.0) : x;
 1744 # ifdef BIGDIG
 1745     ASRTGO(NIMP(y), bady);
 1746     if (BIGP(y)) {
 1747       t = y; y = x; x = t; goto big_dbl;
 1748     }
 1749     ASRTGO(REALP(y), bady);
 1750 # else
 1751     ASRTGO(NIMP(y) && REALP(y), bady);
 1752 # endif
 1753     return (REALPART(x) > REALPART(y)) ? y : x;
 1754   }
 1755   if (NINUMP(y)) {
 1756 # ifdef BIGDIG
 1757     ASRTGO(NIMP(y), bady);
 1758     if (BIGP(y)) return BIGSIGN(y) ? y : x;
 1759 #  ifndef RECKLESS
 1760     if (!(REALP(y)))
 1761       bady: wta(y, (char *)ARG2, s_min);
 1762 #  endif
 1763 # else
 1764 #  ifndef RECKLESS
 1765     if (!(NIMP(y) && REALP(y)))
 1766       bady: wta(y, (char *)ARG2, s_min);
 1767 #  endif
 1768 # endif
 1769     return ((z = INUM(x)) > REALPART(y)) ? y : makdbl(z, 0.0);
 1770   }
 1771 #else
 1772 # ifdef BIGDIG
 1773   if (NINUMP(x)) {
 1774     ASRTER(NIMP(x) && BIGP(x), x, ARG1, s_min);
 1775     if (INUMP(y)) return BIGSIGN(x) ? x : y;
 1776     ASRTGO(NIMP(y) && BIGP(y), bady);
 1777     return (-1==bigcomp(x, y)) ? y : x;
 1778   }
 1779   if (NINUMP(y)) {
 1780 #  ifndef RECKLESS
 1781     if (!(NIMP(y) && BIGP(y)))
 1782       bady: wta(y, (char *)ARG2, s_min);
 1783 #  endif
 1784     return BIGSIGN(y) ? y : x;
 1785   }
 1786 # else
 1787   ASRTER(INUMP(x), x, ARG1, s_min);
 1788   ASRTER(INUMP(y), y, ARG2, s_min);
 1789 # endif
 1790 #endif
 1791   return ((long)x > (long)y) ? y : x;
 1792 }
 1793 
 1794 SCM sum(x, y)
 1795      SCM x, y;
 1796 {
 1797   if (UNBNDP(y)) {
 1798     if (UNBNDP(x)) return INUM0;
 1799 #ifndef RECKLESS
 1800     if (!(NUMBERP(x)))
 1801       badx: wta(x, (char *)ARG1, s_sum);
 1802 #endif
 1803     return x;
 1804   }
 1805 #ifdef FLOATS
 1806   if (NINUMP(x)) {
 1807     SCM t;
 1808 # ifdef BIGDIG
 1809     ASRTGO(NIMP(x), badx);
 1810     if (BIGP(x)) {
 1811       if (INUMP(y)) {t = x; x = y; y = t; goto intbig;}
 1812       ASRTGO(NIMP(y), bady);
 1813       if (BIGP(y)) {
 1814     if (NUMDIGS(x) > NUMDIGS(y)) {t = x; x = y; y = t;}
 1815     return addbig(BDIGITS(x), NUMDIGS(x), BIGSIGN(x), y, 0);
 1816       }
 1817       ASRTGO(INEXP(y), bady);
 1818     bigreal: return makdbl(int2dbl(x)+REALPART(y), CPLXP(y)?IMAG(y):0.0);
 1819     }
 1820     ASRTGO(INEXP(x), badx);
 1821 # else
 1822     ASRTGO(NIMP(x) && INEXP(x), badx);
 1823 # endif
 1824     if (INUMP(y)) {t = x; x = y; y = t; goto intreal;}
 1825 # ifdef BIGDIG
 1826     ASRTGO(NIMP(y), bady);
 1827     if (BIGP(y)) {t = x; x = y; y = t; goto bigreal;}
 1828 #  ifndef RECKLESS
 1829     else if (!(INEXP(y)))
 1830       bady: wta(y, (char *)ARG2, s_sum);
 1831 #  endif
 1832 # else
 1833 #  ifndef RECKLESS
 1834     if (!(NIMP(y) && INEXP(y)))
 1835       bady: wta(y, (char *)ARG2, s_sum);
 1836 #  endif
 1837 # endif
 1838     {
 1839       double i = 0.0;
 1840       if (CPLXP(x)) i = IMAG(x);
 1841       if (CPLXP(y)) i += IMAG(y);
 1842       return makdbl(REALPART(x)+REALPART(y), i);
 1843     }
 1844   }
 1845   if (NINUMP(y)) {
 1846 # ifdef BIGDIG
 1847     ASRTGO(NIMP(y), bady);
 1848     if (BIGP(y))
 1849       intbig: {
 1850 #  ifndef DIGSTOOBIG
 1851       long z = pseudolong(INUM(x));
 1852       return addbig((BIGDIG *)&z, DIGSPERLONG, (x < 0) ? 0x0100 : 0, y, 0);
 1853 #  else
 1854       BIGDIG zdigs[DIGSPERLONG];
 1855       longdigs(INUM(x), zdigs);
 1856       return addbig(zdigs, DIGSPERLONG, (x < 0) ? 0x0100 : 0, y, 0);
 1857 #  endif
 1858     }
 1859     ASRTGO(INEXP(y), bady);
 1860 # else
 1861     ASRTGO(NIMP(y) && INEXP(y), bady);
 1862 # endif
 1863   intreal: return makdbl(INUM(x)+REALPART(y), CPLXP(y)?IMAG(y):0.0);
 1864   }
 1865 #else
 1866 # ifdef BIGDIG
 1867   if (NINUMP(x)) {
 1868     SCM t;
 1869     ASRTGO(NIMP(x) && BIGP(x), badx);
 1870     if (INUMP(y)) {t = x; x = y; y = t; goto intbig;}
 1871     ASRTGO(NIMP(y) && BIGP(y), bady);
 1872     if (NUMDIGS(x) > NUMDIGS(y)) {t = x; x = y; y = t;}
 1873     return addbig(BDIGITS(x), NUMDIGS(x), BIGSIGN(x), y, 0);
 1874   }
 1875   if (NINUMP(y)) {
 1876 #  ifndef RECKLESS
 1877     if (!(NIMP(y) && BIGP(y)))
 1878       bady: wta(y, (char *)ARG2, s_sum);
 1879 #  endif
 1880   intbig: {
 1881 #  ifndef DIGSTOOBIG
 1882       long z = pseudolong(INUM(x));
 1883       return addbig((BIGDIG *)&z, DIGSPERLONG, (x < 0) ? 0x0100 : 0, y, 0);
 1884 #  else
 1885       BIGDIG zdigs[DIGSPERLONG];
 1886       longdigs(INUM(x), zdigs);
 1887       return addbig(zdigs, DIGSPERLONG, (x < 0) ? 0x0100 : 0, y, 0);
 1888 #  endif
 1889     }
 1890   }
 1891 # else
 1892   ASRTGO(INUMP(x), badx);
 1893   ASRTER(INUMP(y), y, ARG2, s_sum);
 1894 # endif
 1895 #endif
 1896   x = INUM(x)+INUM(y);
 1897   if (FIXABLE(x)) return MAKINUM(x);
 1898 #ifdef BIGDIG
 1899   return long2big(x);
 1900 #else
 1901 # ifdef FLOATS
 1902   return makdbl((double)x, 0.0);
 1903 # else
 1904   wta(y, (char *)OVFLOW, s_sum);
 1905 # endif
 1906 #endif
 1907 }
 1908 
 1909 SCM difference(x, y)
 1910      SCM x, y;
 1911 {
 1912 #ifdef FLOATS
 1913   if (NINUMP(x)) {
 1914 # ifndef RECKLESS
 1915     if (!(NIMP(x)))
 1916     badx: wta(x, (char *)ARG1, s_difference);
 1917 # endif
 1918     if (UNBNDP(y)) {
 1919 # ifdef BIGDIG
 1920       if (BIGP(x)) {
 1921     x = copybig(x, !BIGSIGN(x));
 1922     return NUMDIGS(x) * BITSPERDIG/CHAR_BIT <= sizeof(SCM) ?
 1923       big2inum(x, NUMDIGS(x)) : x;
 1924       }
 1925 # endif
 1926       ASRTGO(INEXP(x), badx);
 1927       return makdbl(-REALPART(x), CPLXP(x)?-IMAG(x):0.0);
 1928     }
 1929     if (INUMP(y)) return sum(x, MAKINUM(-INUM(y)));
 1930 # ifdef BIGDIG
 1931     ASRTGO(NIMP(y), bady);
 1932     if (BIGP(x)) {
 1933       if (BIGP(y)) return (NUMDIGS(x) < NUMDIGS(y)) ?
 1934              addbig(BDIGITS(x), NUMDIGS(x), BIGSIGN(x), y, 0x0100) :
 1935              addbig(BDIGITS(y), NUMDIGS(y), BIGSIGN(y) ^ 0x0100, x, 0);
 1936       ASRTGO(INEXP(y), bady);
 1937       return makdbl(int2dbl(x)-REALPART(y), CPLXP(y)?-IMAG(y):0.0);
 1938     }
 1939     ASRTGO(INEXP(x), badx);
 1940     if (BIGP(y)) return makdbl(REALPART(x)-int2dbl(y), CPLXP(x)?IMAG(x):0.0);
 1941     ASRTGO(INEXP(y), bady);
 1942 # else
 1943     ASRTGO(INEXP(x), badx);
 1944     ASRTGO(NIMP(y) && INEXP(y), bady);
 1945 # endif
 1946     if (CPLXP(x)) {
 1947       if (CPLXP(y))
 1948     return makdbl(REAL(x)-REAL(y), IMAG(x)-IMAG(y));
 1949       else
 1950     return makdbl(REAL(x)-REALPART(y), IMAG(x));
 1951     }
 1952     return makdbl(REALPART(x)-REALPART(y), CPLXP(y)?-IMAG(y):0.0);
 1953   }
 1954   if (UNBNDP(y)) {x = -INUM(x); goto checkx;}
 1955   if (NINUMP(y)) {
 1956 # ifdef BIGDIG
 1957     ASRTGO(NIMP(y), bady);
 1958     if (BIGP(y)) {
 1959 #  ifndef DIGSTOOBIG
 1960       long z = pseudolong(INUM(x));
 1961       return addbig((BIGDIG *)&z, DIGSPERLONG, (x < 0) ? 0x0100 : 0, y, 0x0100);
 1962 #  else
 1963       BIGDIG zdigs[DIGSPERLONG];
 1964       longdigs(INUM(x), zdigs);
 1965       return addbig(zdigs, DIGSPERLONG, (x < 0) ? 0x0100 : 0, y, 0x0100);
 1966 #  endif
 1967     }
 1968 #  ifndef RECKLESS
 1969     if (!(INEXP(y)))
 1970     bady: wta(y, (char *)ARG2, s_difference);
 1971 #  endif
 1972 # else
 1973 #  ifndef RECKLESS
 1974     if (!(NIMP(y) && INEXP(y)))
 1975     bady: wta(y, (char *)ARG2, s_difference);
 1976 #  endif
 1977 # endif
 1978     return makdbl(INUM(x)-REALPART(y), CPLXP(y)?-IMAG(y):0.0);
 1979   }
 1980 #else
 1981 # ifdef BIGDIG
 1982   if (NINUMP(x)) {
 1983     ASRTER(NIMP(x) && BIGP(x), x, ARG1, s_difference);
 1984     if (UNBNDP(y)) {
 1985       x = copybig(x, !BIGSIGN(x));
 1986       return NUMDIGS(x) * BITSPERDIG/CHAR_BIT <= sizeof(SCM) ?
 1987     big2inum(x, NUMDIGS(x)) : x;
 1988     }
 1989     if (INUMP(y)) {
 1990 #  ifndef DIGSTOOBIG
 1991       long z = pseudolong(INUM(y));
 1992       return addbig((BIGDIG *)&z, DIGSPERLONG, (y < 0) ? 0 : 0x0100, x, 0);
 1993 #  else
 1994       BIGDIG zdigs[DIGSPERLONG];
 1995       longdigs(INUM(x), zdigs);
 1996       return addbig(zdigs, DIGSPERLONG, (y < 0) ? 0 : 0x0100, x, 0);
 1997 #  endif
 1998     }
 1999     ASRTGO(NIMP(y) && BIGP(y), bady);
 2000     return (NUMDIGS(x) < NUMDIGS(y)) ?
 2001       addbig(BDIGITS(x), NUMDIGS(x), BIGSIGN(x), y, 0x0100) :
 2002       addbig(BDIGITS(y), NUMDIGS(y), BIGSIGN(y) ^ 0x0100, x, 0);
 2003   }
 2004   if (UNBNDP(y)) {x = -INUM(x); goto checkx;}
 2005   if (NINUMP(y)) {
 2006 #  ifndef RECKLESS
 2007     if (!(NIMP(y) && BIGP(y)))
 2008     bady: wta(y, (char *)ARG2, s_difference);
 2009 #  endif
 2010     {
 2011 #  ifndef DIGSTOOBIG
 2012       long z = pseudolong(INUM(x));
 2013       return addbig((BIGDIG *)&z, DIGSPERLONG, (x < 0) ? 0x0100 : 0, y, 0x0100);
 2014 #  else
 2015       BIGDIG zdigs[DIGSPERLONG];
 2016       longdigs(INUM(x), zdigs);
 2017       return addbig(zdigs, DIGSPERLONG, (x < 0) ? 0x0100 : 0, y, 0x0100);
 2018 #  endif
 2019     }
 2020   }
 2021 # else
 2022   ASRTER(INUMP(x), x, ARG1, s_difference);
 2023   if (UNBNDP(y)) {x = -INUM(x); goto checkx;}
 2024   ASRTER(INUMP(y), y, ARG2, s_difference);
 2025 # endif
 2026 #endif
 2027   x = INUM(x)-INUM(y);
 2028  checkx:
 2029   if (FIXABLE(x)) return MAKINUM(x);
 2030 #ifdef BIGDIG
 2031   return long2big(x);
 2032 #else
 2033 # ifdef FLOATS
 2034   return makdbl((double)x, 0.0);
 2035 # else
 2036   wta(y, (char *)OVFLOW, s_difference);
 2037 # endif
 2038 #endif
 2039 }
 2040 
 2041 SCM product(x, y)
 2042      SCM x, y;
 2043 {
 2044   if (UNBNDP(y)) {
 2045     if (UNBNDP(x)) return MAKINUM(1L);
 2046 #ifndef RECKLESS
 2047     if (!(NUMBERP(x)))
 2048       badx: wta(x, (char *)ARG1, s_product);
 2049 #endif
 2050     return x;
 2051   }
 2052 #ifdef FLOATS
 2053   if (NINUMP(x)) {
 2054     SCM t;
 2055 # ifdef BIGDIG
 2056     ASRTGO(NIMP(x), badx);
 2057     if (BIGP(x)) {
 2058       if (INUMP(y)) {t = x; x = y; y = t; goto intbig;}
 2059       ASRTGO(NIMP(y), bady);
 2060       if (BIGP(y)) return mulbig(BDIGITS(x), NUMDIGS(x), BDIGITS(y), NUMDIGS(y),
 2061                  BIGSIGN(x) ^ BIGSIGN(y));
 2062       ASRTGO(INEXP(y), bady);
 2063     bigreal:
 2064       return bigdblop('*', x, REALPART(y), CPLXP(y) ? IMAG(y) : 0.0);
 2065     }
 2066     ASRTGO(INEXP(x), badx);
 2067 # else
 2068     ASRTGO(NIMP(x) && INEXP(x), badx);
 2069 # endif
 2070     if (INUMP(y)) {t = x; x = y; y = t; goto intreal;}
 2071 # ifdef BIGDIG
 2072     ASRTGO(NIMP(y), bady);
 2073     if (BIGP(y)) {t = x; x = y; y = t; goto bigreal;}
 2074 #  ifndef RECKLESS
 2075     else if (!(INEXP(y)))
 2076       bady: wta(y, (char *)ARG2, s_product);
 2077 #  endif
 2078 # else
 2079 #  ifndef RECKLESS
 2080     if (!(NIMP(y) && INEXP(y)))
 2081       bady: wta(y, (char *)ARG2, s_product);
 2082 #  endif
 2083 # endif
 2084     if (CPLXP(x)) {
 2085       if (CPLXP(y))
 2086     return makdbl(REAL(x)*REAL(y)-IMAG(x)*IMAG(y),
 2087               REAL(x)*IMAG(y)+IMAG(x)*REAL(y));
 2088       else
 2089     return makdbl(REAL(x)*REALPART(y), IMAG(x)*REALPART(y));
 2090     }
 2091     return makdbl(REALPART(x)*REALPART(y),
 2092           CPLXP(y)?REALPART(x)*IMAG(y):0.0);
 2093   }
 2094   if (NINUMP(y)) {
 2095 # ifdef BIGDIG
 2096     ASRTGO(NIMP(y), bady);
 2097     if (BIGP(y)) {
 2098     intbig: if (INUM0==x) return x; if (MAKINUM(1L)==x) return y;
 2099     {
 2100 #  ifndef DIGSTOOBIG
 2101       long z = pseudolong(INUM(x));
 2102       return mulbig((BIGDIG *)&z, DIGSPERLONG, BDIGITS(y), NUMDIGS(y),
 2103             BIGSIGN(y) ? (x>0) : (x<0));
 2104 #  else
 2105       BIGDIG zdigs[DIGSPERLONG];
 2106       longdigs(INUM(x), zdigs);
 2107       return mulbig(zdigs, DIGSPERLONG, BDIGITS(y), NUMDIGS(y),
 2108             BIGSIGN(y) ? (x>0) : (x<0));
 2109 #  endif
 2110     }
 2111     }
 2112     ASRTGO(INEXP(y), bady);
 2113 # else
 2114     ASRTGO(NIMP(y) && INEXP(y), bady);
 2115 # endif
 2116   intreal: return makdbl(INUM(x)*REALPART(y), CPLXP(y)?INUM(x)*IMAG(y):0.0);
 2117   }
 2118 #else
 2119 # ifdef BIGDIG
 2120   if (NINUMP(x)) {
 2121     ASRTGO(NIMP(x) && BIGP(x), badx);
 2122     if (INUMP(y)) {SCM t = x; x = y; y = t; goto intbig;}
 2123     ASRTGO(NIMP(y) && BIGP(y), bady);
 2124     return mulbig(BDIGITS(x), NUMDIGS(x), BDIGITS(y), NUMDIGS(y),
 2125           BIGSIGN(x) ^ BIGSIGN(y));
 2126   }
 2127   if (NINUMP(y)) {
 2128 #  ifndef RECKLESS
 2129     if (!(NIMP(y) && BIGP(y)))
 2130       bady: wta(y, (char *)ARG2, s_product);
 2131 #  endif
 2132   intbig: if (INUM0==x) return x; if (MAKINUM(1L)==x) return y;
 2133   {
 2134 #  ifndef DIGSTOOBIG
 2135     long z = pseudolong(INUM(x));
 2136     return mulbig((BIGDIG *)&z, DIGSPERLONG, BDIGITS(y), NUMDIGS(y),
 2137           BIGSIGN(y) ? (x>0) : (x<0));
 2138 #  else
 2139     BIGDIG zdigs[DIGSPERLONG];
 2140     longdigs(INUM(x), zdigs);
 2141     return mulbig(zdigs, DIGSPERLONG, BDIGITS(y), NUMDIGS(y),
 2142           BIGSIGN(y) ? (x>0) : (x<0));
 2143 #  endif
 2144   }
 2145   }
 2146 # else
 2147   ASRTGO(INUMP(x), badx);
 2148   ASRTER(INUMP(y), y, ARG2, s_product);
 2149 # endif
 2150 #endif
 2151   {
 2152     long i, j, k;
 2153     i = INUM(x);
 2154     if (0==i) return x;
 2155     j = INUM(y);
 2156     k = i * j;
 2157     y = MAKINUM(k);
 2158     if (k != INUM(y) || k/i != j)
 2159 #ifdef BIGDIG
 2160       {
 2161     int sgn = (i < 0) ^ (j < 0);
 2162 # ifndef DIGSTOOBIG
 2163     i = pseudolong(i);
 2164     j = pseudolong(j);
 2165     return mulbig((BIGDIG *)&i, DIGSPERLONG,
 2166               (BIGDIG *)&j, DIGSPERLONG, sgn);
 2167 # else /* DIGSTOOBIG */
 2168     BIGDIG idigs[DIGSPERLONG];
 2169     BIGDIG jdigs[DIGSPERLONG];
 2170     longdigs(i, idigs);
 2171     longdigs(j, jdigs);
 2172     return mulbig(idigs, DIGSPERLONG, jdigs, DIGSPERLONG, sgn);
 2173 # endif
 2174       }
 2175 #else
 2176 # ifdef FLOATS
 2177     return makdbl(((double)i)*((double)j), 0.0);
 2178 # else
 2179     wta(y, (char *)OVFLOW, s_product);
 2180 # endif
 2181 #endif
 2182     return y;
 2183   }
 2184 }
 2185   /* Use "Smith's formula" to extend dynamic range */
 2186   /* David Goldberg
 2187      What Every Computer Scientist Should Know About Floating-Point Arithmetic
 2188      http://cch.loria.fr/documentation/IEEE754/ACM/goldberg.pdf */
 2189 SCM divide(x, y)
 2190      SCM x, y;
 2191 {
 2192 #ifdef FLOATS
 2193   double den, a = 1.0;
 2194   if (NINUMP(x)) {
 2195 # ifndef RECKLESS
 2196     if (!(NIMP(x)))
 2197       badx: wta(x, (char *)ARG1, s_divide);
 2198 # endif
 2199     if (UNBNDP(y)) {
 2200 # ifdef BIGDIG
 2201       if (BIGP(x)) return inex_divintbig(MAKINUM(1L), x);
 2202 # endif
 2203       /* reciprocal */
 2204       ASRTGO(INEXP(x), badx);
 2205       if (REALP(x)) return makdbl(1.0/REALPART(x), 0.0);
 2206       {
 2207     y = x;
 2208     a = 1.0;
 2209     goto real_over_complex;
 2210       }
 2211     }
 2212 # ifdef BIGDIG
 2213     if (BIGP(x)) {
 2214       SCM z;
 2215       if (INUMP(y)) {
 2216         z = INUM(y);
 2217         ASRTER(z, y, OVFLOW, s_divide);
 2218     if (1==z) return x;
 2219         if (z < 0) z = -z;
 2220         if (z < BIGRAD) {
 2221           SCM w = copybig(x, BIGSIGN(x) ? (y>0) : (y<0));
 2222       int sts = divbigdig(BDIGITS(w), NUMDIGS(w), (BIGDIG)z);
 2223       if (sts) {
 2224         bigrecy(w);
 2225         return bigdblop('/', x, INUM(y), 0.0);
 2226       }
 2227       else return normbig(w);
 2228     }
 2229 #  ifndef DIGSTOOBIG
 2230         z = pseudolong(z);
 2231         z = divbigbig(BDIGITS(x), NUMDIGS(x), (BIGDIG *)&z, DIGSPERLONG,
 2232                       BIGSIGN(x) ? (y>0) : (y<0), 4);
 2233 #  else
 2234     {
 2235       BIGDIG zdigs[DIGSPERLONG];
 2236       longdigs(z, zdigs);
 2237       z = divbigbig(BDIGITS(x), NUMDIGS(x), zdigs, DIGSPERLONG,
 2238             BIGSIGN(x) ? (y>0) : (y<0), 4);
 2239     }
 2240 #  endif
 2241         return z ? z : bigdblop('/', x, INUM(y), 0.0);
 2242       }
 2243       ASRTGO(NIMP(y), bady);
 2244       if (BIGP(y)) {
 2245     z = divbigbig(BDIGITS(x), NUMDIGS(x), BDIGITS(y), NUMDIGS(y),
 2246               BIGSIGN(x) ^ BIGSIGN(y), 4);
 2247     return z ? z : inex_divintbig(x, y);
 2248       }
 2249       ASRTGO(INEXP(y), bady);
 2250       return bigdblop('/', x, REALPART(y), CPLXP(y) ? IMAG(y) : 0.0);
 2251     }
 2252 # endif
 2253     ASRTGO(INEXP(x), badx);
 2254     if (INUMP(y)) {den = INUM(y); goto basic_div;}
 2255 # ifdef BIGDIG
 2256     ASRTGO(NIMP(y), bady);
 2257     if (BIGP(y)) return bigdblop('\\', y, REALPART(x), CPLXP(x) ? IMAG(x) : 0.0);
 2258     ASRTGO(INEXP(y), bady);
 2259 # else
 2260     ASRTGO(NIMP(y) && INEXP(y), bady);
 2261 # endif
 2262     if (REALP(y)) {
 2263       den = REALPART(y);
 2264     basic_div: return makdbl(REALPART(x)/den, CPLXP(x)?IMAG(x)/den:0.0);
 2265     }
 2266     a = REALPART(x);
 2267     if (REALP(x)) goto real_over_complex;
 2268     /* Both x and y are complex */
 2269     /* Use "Smith's formula" to extend dynamic range */
 2270     {
 2271       double b = IMAG(x);
 2272       double c = REALPART(y);
 2273       double d = IMAG(y);
 2274       if ((d > 0 ? d : -d) < (c > 0 ? c : -c)) {
 2275     double r = d/c;
 2276     double i = c + d*r;
 2277     return makdbl((a + b*r)/i, (b - a*r)/i);
 2278       }
 2279       {
 2280     double r = c/d;
 2281     double i = d + c*r;
 2282     return makdbl((b + a*r)/i, (-a + b*r)/i);
 2283       }
 2284     }
 2285   }
 2286   if (UNBNDP(y)) {
 2287     if ((MAKINUM(1L)==x) || (MAKINUM(-1L)==x)) return x;
 2288     return makdbl(1.0/((double)INUM(x)), 0.0);
 2289   }
 2290   if (NINUMP(y)) {
 2291 # ifdef BIGDIG
 2292     ASRTGO(NIMP(y), bady);
 2293     if (BIGP(y)) return inex_divintbig(x, y); /* bigdblop('\\', y, INUM(x), 0.0); */
 2294 #  ifndef RECKLESS
 2295     if (!(INEXP(y)))
 2296       bady: wta(y, (char *)ARG2, s_divide);
 2297 #  endif
 2298 # else
 2299 #  ifndef RECKLESS
 2300     if (!(NIMP(y) && INEXP(y)))
 2301       bady: wta(y, (char *)ARG2, s_divide);
 2302 #  endif
 2303 # endif
 2304     if (REALP(y)) return makdbl(INUM(x)/REALPART(y), 0.0);
 2305     a = INUM(x);
 2306   real_over_complex:
 2307     /* Both x and y are complex */
 2308     /* Use "Smith's formula" to extend dynamic range */
 2309     {
 2310       double c = REALPART(y);
 2311       double d = IMAG(y);
 2312       if ((d > 0 ? d : -d) < (c > 0 ? c : -c)) {
 2313     double r = d/c;
 2314     double i = c + d*r;
 2315     return makdbl((a)/i, (- a*r)/i);
 2316       }
 2317       {
 2318     double r = c/d;
 2319     double i = d + c*r;
 2320     return makdbl((a*r)/i, (-a)/i);
 2321       }
 2322     }
 2323   }
 2324 #else
 2325 # ifdef BIGDIG
 2326   if (NINUMP(x)) {
 2327     SCM z;
 2328     ASRTER(NIMP(x) && BIGP(x), x, ARG1, s_divide);
 2329     if (UNBNDP(y)) goto ov;
 2330     if (INUMP(y)) {
 2331       z = INUM(y);
 2332       if (!z) goto ov;
 2333       if (1==z) return x;
 2334       if (z < 0) z = -z;
 2335       if (z < BIGRAD) {
 2336         SCM w = copybig(x, BIGSIGN(x) ? (y>0) : (y<0));
 2337     int sts = divbigdig(BDIGITS(w), NUMDIGS(w), (BIGDIG)z);
 2338         if (sts) {
 2339       bigrecy(w);
 2340       goto ov;
 2341     }
 2342         return w;
 2343       }
 2344 #  ifndef DIGSTOOBIG
 2345       z = pseudolong(z);
 2346       z = divbigbig(BDIGITS(x), NUMDIGS(x), (BIGDIG *)&z, DIGSPERLONG,
 2347             BIGSIGN(x) ? (y>0) : (y<0), 4);
 2348 #  else
 2349       {
 2350     BIGDIG zdigs[DIGSPERLONG];
 2351     longdigs(z, zdigs);
 2352     z = divbigbig(BDIGITS(x), NUMDIGS(x), zdigs, DIGSPERLONG,
 2353               BIGSIGN(x) ? (y>0) : (y<0), 4);
 2354       }
 2355 #  endif
 2356     } else {
 2357       ASRTGO(NIMP(y) && BIGP(y), bady);
 2358       z = divbigbig(BDIGITS(x), NUMDIGS(x), BDIGITS(y), NUMDIGS(y),
 2359             BIGSIGN(x) ^ BIGSIGN(y), 4);
 2360     }
 2361     if (!z) goto ov;
 2362     return z;
 2363   }
 2364   if (UNBNDP(y)) {
 2365     if ((MAKINUM(1L)==x) || (MAKINUM(-1L)==x)) return x;
 2366     goto ov;
 2367   }
 2368   if (NINUMP(y)) {
 2369 #  ifndef RECKLESS
 2370     if (!(NIMP(y) && BIGP(y)))
 2371       bady: wta(y, (char *)ARG2, s_divide);
 2372 #  endif
 2373     goto ov;
 2374   }
 2375 # else
 2376   ASRTER(INUMP(x), x, ARG1, s_divide);
 2377   if (UNBNDP(y)) {
 2378     if ((MAKINUM(1L)==x) || (MAKINUM(-1L)==x)) return x;
 2379     goto ov;
 2380   }
 2381   ASRTER(INUMP(y), y, ARG2, s_divide);
 2382 # endif
 2383 #endif
 2384   {
 2385     long z = INUM(y);
 2386     if ((0==z) || INUM(x)%z) goto ov;
 2387     z = INUM(x)/z;
 2388     if (FIXABLE(z)) return MAKINUM(z);
 2389 #ifdef BIGDIG
 2390     return long2big(z);
 2391 #endif
 2392 #ifdef FLOATS
 2393   ov: return makdbl(((double)INUM(x))/((double)INUM(y)), 0.0);
 2394 #else
 2395   ov: wta(x, (char *)OVFLOW, s_divide);
 2396 #endif
 2397   }
 2398 }
 2399 
 2400 SCM ilog(m, b, k, n)
 2401      unsigned long m;
 2402      SCM b, k;
 2403      unsigned long *n;
 2404 {
 2405   /* printf("call ilog %ld ", m); scm_iprin1(b, cur_outp, 0); printf(" "); scm_iprin1(k, cur_outp, 0); printf("\n"); */
 2406   if (BOOL_T==lessp(k, b)) return k;
 2407   *n += m;
 2408   {
 2409     SCM q = ilog(2*m, product(b, b), lquotient(k, b), n);
 2410     /* printf("return ilog "); scm_iprin1(q, cur_outp, 0); printf("\n"); */
 2411     if (BOOL_T==lessp(q, b)) return q;
 2412     *n += m;
 2413     return lquotient(q, b);
 2414   }
 2415 }
 2416 
 2417 SCM scm_intlog(base, k)
 2418      SCM base, k;
 2419 {
 2420   unsigned long n = 1;
 2421   ASRTER(INUMP(base) || (NIMP(base) && BIGP(base)), base, ARG1, s_intlog);
 2422   ASRTER(BOOL_T==lessp(MAKINUM(1L), base), base, OUTOFRANGE, s_intlog);
 2423   ASRTER((INUMP(k) && k > 0) || (NIMP(k) && TYP16(k)==tc16_bigpos), k, ARG2, s_intlog);
 2424   if (BOOL_T==lessp(k, base)) return INUM0;
 2425   ilog(1, base, lquotient(k, base), &n);
 2426   return MAKINUM(n);
 2427 }
 2428 
 2429 #ifdef INUMS_ONLY
 2430 # define eqv eqp
 2431 #endif
 2432 SCM scm_cintlog(base, k)
 2433      SCM base, k;
 2434 {
 2435   SCM il = scm_intlog(base, k);
 2436   return (BOOL_T==eqv(k, scm_intexpt(base, il))) ? il : sum(MAKINUM(1L), il);
 2437 }
 2438 
 2439 SCM scm_intexpt(z1, z2)
 2440      SCM z1, z2;
 2441 {
 2442   SCM acc = MAKINUM(1L);
 2443   long iz2;
 2444 #ifdef FLOATS
 2445   double dacc, dz1;
 2446 #endif
 2447   if (INUM0==z2) return sum(acc, product(z1, INUM0));
 2448   ASRTER(INUMP(z2), z2, ARG2, s_intexpt);
 2449   if (acc==z1) return z1;
 2450   if (MAKINUM(-1L)==z1) return BOOL_F==evenp(z2)?z1:acc;
 2451   iz2 = INUM(z2);
 2452   if (iz2 < 0L) {
 2453     iz2 = -iz2;
 2454     z1 = divide(z1, UNDEFINED);
 2455   }
 2456   if (INUMP(z1)) {
 2457     long tmp, iacc = 1L, iz1 = INUM(z1);
 2458     while (1) {
 2459       if (0L==iz2) {
 2460     acc = long2num(iacc);
 2461     break;
 2462       }
 2463       if (0L==iz1) return z1;
 2464       if (1L==iz2) {
 2465     tmp = iacc*iz1;
 2466     if (tmp/iacc != iz1) {
 2467     overflow:
 2468       z1 = long2num(iz1);
 2469       acc = long2num(iacc);
 2470       ASRTGO(NFALSEP(z1) && NFALSEP(acc), errout);
 2471       goto gencase;
 2472     }
 2473     acc = long2num(tmp);
 2474     break;
 2475       }
 2476       if (iz2 & 1) {
 2477     tmp = iacc*iz1;
 2478     if (tmp/iacc != iz1) goto overflow;
 2479     iacc = tmp;
 2480     iz2 = iz2 - 1L;     /* so jumping to gencase works  */
 2481       }
 2482       tmp = iz1*iz1;
 2483       if (tmp/iz1 != iz1) goto overflow;
 2484       iz1 = tmp;
 2485       iz2 >>= 1;
 2486     }
 2487 #ifndef RECKLESS
 2488     if (FALSEP(acc))
 2489     errout: wta(UNDEFINED, (char *)OVFLOW, s_intexpt);
 2490 #endif
 2491     goto ret;
 2492   }
 2493   ASRTER(NIMP(z1), z1, ARG1, s_intexpt);
 2494 #ifdef FLOATS
 2495   if (REALP(z1)) {
 2496     dz1 = REALPART(z1);
 2497     dacc = 1.0;
 2498     while(1) {
 2499       if (0L==iz2) break;
 2500       if (1L==iz2) {dacc = dacc*dz1; break;}
 2501       if (iz2 & 1) dacc = dacc*dz1;
 2502       dz1 = dz1*dz1;
 2503       iz2 >>= 1;
 2504     }
 2505     return makdbl(dacc, 0.0);
 2506   }
 2507 #endif
 2508  gencase:
 2509   {
 2510     SCM tz1 = z1;
 2511     while(1) {
 2512       SCM tmp;
 2513       if (0L==iz2) break;
 2514       if (1L==iz2) {
 2515     tmp = acc;
 2516     acc = product(tmp, tz1);
 2517 #ifdef BIGDIG
 2518     if (acc != tmp) bigrecy(tmp);
 2519 #endif
 2520     break;
 2521       }
 2522       if (iz2 & 1) {
 2523     tmp = acc;
 2524     acc = product(tmp, tz1);
 2525 #ifdef BIGDIG
 2526     if (acc != tmp) bigrecy(tmp);
 2527 #endif
 2528       }
 2529       tmp = tz1;
 2530       tz1 = product(tmp, tmp);
 2531 #ifdef BIGDIG
 2532       if (tmp != z1 && tmp != acc) bigrecy(tmp);
 2533 #endif
 2534       iz2 >>= 1;
 2535     }
 2536 #ifdef BIGDIG
 2537     if (tz1 != acc && tz1 != z1) bigrecy(tz1);
 2538 #endif
 2539   }
 2540  ret: return acc;
 2541 }
 2542 
 2543 #ifdef FLOATS
 2544 # ifndef HAVE_ATANH
 2545 double asinh(x)
 2546      double x;
 2547 {
 2548   return log(x+sqrt(x*x+1));
 2549 }
 2550 
 2551 double acosh(x)
 2552      double x;
 2553 {
 2554   return log(x+sqrt(x*x-1));
 2555 }
 2556 
 2557 double atanh(x)
 2558      double x;
 2559 {
 2560   return 0.5*log((1+x)/(1-x));
 2561 }
 2562 # endif
 2563 
 2564 double scm_truncate(x)
 2565      double x;
 2566 {
 2567   if (x < 0.0) return -floor(-x);
 2568   return floor(x);
 2569 }
 2570 
 2571 double scm_round(x)
 2572      double x;
 2573 {
 2574   double plus_half = x + 0.5;
 2575   double result = floor(plus_half);
 2576   /* Adjust so that the round is towards even.  */
 2577   return (plus_half==result && plus_half / 2 != floor(plus_half / 2))
 2578     ? result - 1 : result;
 2579 }
 2580 
 2581 struct dpair {double x, y;};
 2582 
 2583 void two_doubles(z1, z2, sstring, xy)
 2584      SCM z1, z2;
 2585      char *sstring;
 2586      struct dpair *xy;
 2587 {
 2588   if (INUMP(z1)) xy->x = INUM(z1);
 2589   else {
 2590 # ifdef BIGDIG
 2591     ASRTGO(NIMP(z1), badz1);
 2592     if (BIGP(z1)) xy->x = int2dbl(z1);
 2593     else {
 2594 #  ifndef RECKLESS
 2595       if (!(REALP(z1)))
 2596     badz1: wta(z1, (char *)ARG1, sstring);
 2597 #  endif
 2598       xy->x = REALPART(z1);}
 2599 # else
 2600     {ASRTER(NIMP(z1) && REALP(z1), z1, ARG1, sstring);
 2601     xy->x = REALPART(z1);}
 2602 # endif
 2603   }
 2604   if (INUMP(z2)) xy->y = INUM(z2);
 2605   else {
 2606 # ifdef BIGDIG
 2607     ASRTGO(NIMP(z2), badz2);
 2608     if (BIGP(z2)) xy->y = int2dbl(z2);
 2609     else {
 2610 #  ifndef RECKLESS
 2611       if (!(REALP(z2)))
 2612     badz2: wta(z2, (char *)ARG2, sstring);
 2613 #  endif
 2614       xy->y = REALPART(z2);}
 2615 # else
 2616     {
 2617       ASRTER(NIMP(z2) && REALP(z2), z2, ARG2, sstring);
 2618       xy->y = REALPART(z2);
 2619     }
 2620 # endif
 2621   }
 2622 }
 2623 
 2624 SCM expt(z1, z2)
 2625      SCM z1, z2;
 2626 {
 2627   struct dpair xy;
 2628   two_doubles(z1, z2, s_expt, &xy);
 2629   return makdbl(pow(xy.x, xy.y), 0.0);
 2630 }
 2631 SCM latan2(z1, z2)
 2632      SCM z1, z2;
 2633 {
 2634   struct dpair xy;
 2635   two_doubles(z1, z2, s_atan2, &xy);
 2636   return makdbl(atan2(xy.x, xy.y), 0.0);
 2637 }
 2638 SCM makrect(z1, z2)
 2639      SCM z1, z2;
 2640 {
 2641   struct dpair xy;
 2642   two_doubles(z1, z2, s_makrect, &xy);
 2643   return makdbl(xy.x, xy.y);
 2644 }
 2645 SCM makpolar(z1, z2)
 2646      SCM z1, z2;
 2647 {
 2648   struct dpair xy;
 2649   two_doubles(z1, z2, s_makpolar, &xy);
 2650   return makdbl(xy.x*cos(xy.y), xy.x*sin(xy.y));
 2651 }
 2652 
 2653 SCM real_part(z)
 2654      SCM z;
 2655 {
 2656   if (NINUMP(z)) {
 2657 # ifdef BIGDIG
 2658     ASRTGO(NIMP(z), badz);
 2659     if (BIGP(z)) return z;
 2660 #  ifndef RECKLESS
 2661     if (!(INEXP(z)))
 2662       badz: wta(z, (char *)ARG1, s_real_part);
 2663 #  endif
 2664 # else
 2665     ASRTER(NIMP(z) && INEXP(z), z, ARG1, s_real_part);
 2666 # endif
 2667     if (CPLXP(z)) return makdbl(REAL(z), 0.0);
 2668   }
 2669   return z;
 2670 }
 2671 SCM imag_part(z)
 2672      SCM z;
 2673 {
 2674   if (INUMP(z)) return INUM0;
 2675 # ifdef BIGDIG
 2676   ASRTGO(NIMP(z), badz);
 2677   if (BIGP(z)) return INUM0;
 2678 #  ifndef RECKLESS
 2679   if (!(INEXP(z)))
 2680     badz: wta(z, (char *)ARG1, s_imag_part);
 2681 #  endif
 2682 # else
 2683   ASRTER(NIMP(z) && INEXP(z), z, ARG1, s_imag_part);
 2684 # endif
 2685   if (CPLXP(z)) return makdbl(IMAG(z), 0.0);
 2686   return flo0;
 2687 }
 2688 
 2689 SCM scm_abs(z)
 2690      SCM z;
 2691 {
 2692   if (INUMP(z)) return scm_iabs(z);
 2693   ASRTGO(NIMP(z), badz);
 2694 # ifdef BIGDIG
 2695   if (BIGP(z)) return scm_iabs(z);
 2696 # endif
 2697   if (!REALP(z))
 2698     badz: wta(z, (char *)ARG1, s_abs);
 2699   return makdbl(fabs(REALPART(z)), 0.0);
 2700 }
 2701 
 2702 SCM scm_magnitude(z)
 2703      SCM z;
 2704 {
 2705   if (INUMP(z)) return scm_iabs(z);
 2706   ASRTGO(NIMP(z), badz);
 2707 # ifdef BIGDIG
 2708   if (BIGP(z)) return scm_iabs(z);
 2709 # endif
 2710   if (!INEXP(z))
 2711     badz: wta(z, (char *)ARG1, s_magnitude);
 2712   if (CPLXP(z))
 2713     {
 2714       double i = IMAG(z), r = REAL(z);
 2715       if (i < 0) i = -i;
 2716       if (r < 0) r = -r;
 2717       if (i < r) {
 2718     double q = i / r;
 2719     return makdbl(r * sqrt(1 + q * q), 0.0);
 2720       }
 2721       if (0.0==i) return i;
 2722       {
 2723     double q = r / i;
 2724     return makdbl(i * sqrt(1 + q * q), 0.0);
 2725       }
 2726     }
 2727   return makdbl(fabs(REALPART(z)), 0.0);
 2728 }
 2729 
 2730 SCM angle(z)
 2731      SCM z;
 2732 {
 2733   double x, y = 0.0;
 2734   if (INUMP(z)) {x = (z>=INUM0) ? 1.0 : -1.0; goto do_angle;}
 2735 # ifdef BIGDIG
 2736   ASRTGO(NIMP(z), badz);
 2737   if (BIGP(z)) {x = (TYP16(z)==tc16_bigpos) ? 1.0 : -1.0; goto do_angle;}
 2738 #  ifndef RECKLESS
 2739   if (!(INEXP(z))) {
 2740     badz: wta(z, (char *)ARG1, s_angle);}
 2741 #  endif
 2742 # else
 2743   ASRTER(NIMP(z) && INEXP(z), z, ARG1, s_angle);
 2744 # endif
 2745   if (REALP(z)) {x = REALPART(z); goto do_angle;}
 2746   x = REAL(z); y = IMAG(z);
 2747 do_angle:
 2748   return makdbl(atan2(y, x), 0.0);
 2749 }
 2750 
 2751 
 2752 SCM ex2in(z)
 2753      SCM z;
 2754 {
 2755   if (INUMP(z)) return makdbl((double)INUM(z), 0.0);
 2756   ASRTGO(NIMP(z), badz);
 2757   if (INEXP(z)) return z;
 2758 # ifdef BIGDIG
 2759   if (BIGP(z)) return makdbl(int2dbl(z), 0.0);
 2760 # endif
 2761  badz: wta(z, (char *)ARG1, s_ex2in);
 2762 }
 2763 SCM in2ex(z)
 2764      SCM z;
 2765 {
 2766   if (INUMP(z)) return z;
 2767 # ifdef BIGDIG
 2768   ASRTGO(NIMP(z), badz);
 2769   if (BIGP(z)) return z;
 2770 #  ifndef RECKLESS
 2771   if (!(REALP(z)))
 2772     badz: wta(z, (char *)ARG1, s_in2ex);
 2773 #  endif
 2774 # else
 2775   ASRTER(NIMP(z) && REALP(z), z, ARG1, s_in2ex);
 2776 # endif
 2777 # ifdef BIGDIG
 2778   {
 2779     double u = floor(REALPART(z)+0.5);
 2780     if ((u <= MOST_POSITIVE_FIXNUM)
 2781 #  ifdef hpux
 2782     && (-u <= -MOST_NEGATIVE_FIXNUM) /* workaround for HP700 cc bug */
 2783 #  endif
 2784     ) {
 2785       SCM ans = MAKINUM((long)u);
 2786       if (INUM(ans)==(long)u) return ans;
 2787     }
 2788     ASRTGO(!((u==2*u) || (u)!=(u)), badz); /* problem? */
 2789     return dbl2big(u);
 2790   }
 2791 # else
 2792   return MAKINUM((long)floor(REALPART(z)+0.5));
 2793 # endif
 2794 }
 2795 #else               /* ~FLOATS */
 2796 static char s_trunc[] = "truncate";
 2797 SCM numident(x)
 2798      SCM x;
 2799 {
 2800 # ifdef BIGDIG
 2801   ASRTER(INUMP(x) || (NIMP(x) && BIGP(x)), x, ARG1, s_trunc);
 2802 # else
 2803   ASRTER(INUMP(x), x, ARG1, s_trunc);
 2804 # endif
 2805   return x;
 2806 }
 2807 #endif              /* FLOATS */
 2808 
 2809 SCM scm_iabs(x)
 2810      SCM x;
 2811 {
 2812 #ifdef BIGDIG
 2813   if (NINUMP(x)) {
 2814     ASRTER(NIMP(x) && BIGP(x), x, ARG1, s_abs);
 2815     if (TYP16(x)==tc16_bigpos) return x;
 2816     return copybig(x, 0);
 2817   }
 2818 #else
 2819   ASRTER(INUMP(x), x, ARG1, s_abs);
 2820 #endif
 2821   if (INUM(x) >= 0) return x;
 2822   x = -INUM(x);
 2823   if (!POSFIXABLE(x))
 2824 #ifdef BIGDIG
 2825     return long2big(x);
 2826 #else
 2827     wta(MAKINUM(-x), (char *)OVFLOW, s_abs);
 2828 #endif
 2829   return MAKINUM(x);
 2830 }
 2831 
 2832 #ifdef BIGDIG
 2833 # ifdef FLOATS
 2834 SCM dbl2big(d)
 2835      double d;          /* must be integer */
 2836 {
 2837   sizet i = 0;
 2838   long c;
 2839   BIGDIG *digits;
 2840   SCM ans;
 2841   double u = (d < 0)?-d:d;
 2842   ASRTER(u == u && u != u/2, makdbl(d, 0.0), OVFLOW, "dbl2big");
 2843   while (0 != floor(u)) {u /= BIGRAD;i++;}
 2844   ans = mkbig(i, d < 0);
 2845   digits = BDIGITS(ans);
 2846   while (i--) {
 2847     u = ldexp(u, BITSPERDIG);
 2848     c = floor(u);
 2849     u -= c;
 2850     digits[i] = c;
 2851   }
 2852   ASRTER(0==u, makdbl(d, 0.0), OVFLOW, "dbl2big");
 2853   return normbig(ans);
 2854 }
 2855 /* This turns out to need explicit rounding for bignums */
 2856 double int2dbl(b)
 2857      SCM b;
 2858 {
 2859   if (INUMP(b)) return (double)(INUM(b));
 2860   {
 2861     SCM num = scm_iabs(b), quo = num;
 2862     int bex = INUM(scm_intlength(num)) - dbl_mant_dig;
 2863     double ans = 0.0;
 2864     if (bex > 0) {
 2865       SCM den = scm_ash(MAKINUM(1L), MAKINUM(bex));
 2866       quo = scm_round_quotient(num, den);
 2867       bigrecy(den);
 2868     }
 2869     /* quo may not be a bignum */
 2870     if (INUMP(quo))
 2871       ans = ldexp((double)(INUM(quo)), bex);
 2872     else {
 2873       sizet i = NUMDIGS(quo);
 2874       sizet j = (dbl_mant_dig + BITSPERDIG - 1)/BITSPERDIG;
 2875       BIGDIG *digits = BDIGITS(quo);
 2876       if (j < i) j = i - j;
 2877       else j = 0;
 2878       while (i-- > j) ans = digits[i] + ldexp(ans, BITSPERDIG);
 2879       bex += j * BITSPERDIG;
 2880       if (bex > 0) ans = ldexp(ans, bex);
 2881     }
 2882     if (quo != num) bigrecy(quo);
 2883     if (num != b) bigrecy(num);
 2884     if (tc16_bigneg==TYP16(b)) return -ans;
 2885     return ans;
 2886   }
 2887 }
 2888 static SCM bigdblop(op, b, re, im)
 2889      int op;
 2890      SCM b;
 2891      double re, im;
 2892 {
 2893   double bm = 0.0;
 2894   int i = 0;
 2895   if (NUMDIGS(b)*BITSPERDIG < DBL_MAX_EXP) {
 2896     bm = int2dbl(b);
 2897   }
 2898   else {
 2899     i = INUM(scm_intlength(b));
 2900     if (i < DBL_MAX_EXP) {
 2901       i = 0;
 2902       bm = int2dbl(b);
 2903     }
 2904     else {
 2905       i = i + 1 - DBL_MAX_EXP;
 2906       bm = ldexp(int2dbl(b), -i);
 2907     }
 2908   }
 2909   switch (op) {
 2910   case '*':
 2911     return makdbl(ldexp(bm*re, i), 0.0==im ? 0.0 : ldexp(bm*im, i));
 2912   case '/':
 2913     if (0.0==im) return makdbl(bm/re, 0.0);
 2914     {
 2915       double d = re*re + im*im;
 2916       return makdbl(ldexp(bm*(re/d), i), ldexp(-bm*(im/d), i));
 2917     }
 2918   case '\\':
 2919     return makdbl(ldexp(re/bm, -i), 0.0==im ? 0.0 : ldexp(im/bm, -i));
 2920   default:
 2921     return UNSPECIFIED;
 2922   }
 2923 }
 2924 /* now able to return unnomalized doubles. */
 2925 static SCM inex_divintbig(a, b)
 2926      SCM a, b;
 2927 {
 2928   double r;
 2929   {
 2930     int sgn = (((INUMP(a) ? (INUM(a) < 0):BIGSIGN(a))==0) ^
 2931            (BIGSIGN(b)==0)) ? -1 : 1;
 2932     SCM ma = scm_abs(a);
 2933     SCM mb = scm_abs(b);
 2934     int la = INUM(scm_intlength(ma));
 2935     int lb = INUM(scm_intlength(mb));
 2936     if (la <= DBL_MAX_EXP && lb <= DBL_MAX_EXP) {
 2937       r = int2dbl(a) / int2dbl(b);
 2938     }
 2939     else if (la > DBL_MAX_EXP && lb > DBL_MAX_EXP) {
 2940       int k = (la > lb ? la : lb) - DBL_MAX_EXP;
 2941       r = sgn *
 2942     int2dbl(scm_ash(ma, MAKINUM(-k))) /
 2943     int2dbl(scm_ash(mb, MAKINUM(-k)));
 2944     } else if (la > lb) {
 2945       int k = la - DBL_MAX_EXP;
 2946       r = sgn * ldexp(int2dbl(scm_ash(ma, MAKINUM(-k))) / int2dbl(mb), k);
 2947     } else {
 2948       int k = lb - DBL_MAX_EXP;
 2949       r = sgn * ldexp(int2dbl(ma) / int2dbl(scm_ash(mb, MAKINUM(-k))), -k);
 2950     }
 2951   }
 2952   return makdbl(r, 0.0);
 2953 }
 2954 
 2955 SCM scm_dfloat_parts(f)
 2956      SCM f;
 2957 {
 2958   int expt, ndig = dbl_mant_dig;
 2959   double mant = frexp(num2dbl(f, (char *)ARG1, s_dfloat_parts), &expt);
 2960 #  ifdef DBL_MIN_EXP
 2961   if (expt < DBL_MIN_EXP)
 2962     ndig -= DBL_MIN_EXP - expt;
 2963 #  endif
 2964   mant *= ldexp(1.0, ndig);
 2965   expt -= ndig;
 2966   return scm_values(dbl2big(mant), MAKINUM(expt), EOL, s_dfloat_parts);
 2967 }
 2968 static char s_make_dfloat[] = "make-double-float";
 2969 SCM scm_make_dfloat(mant, expt)
 2970      SCM mant, expt;
 2971 {
 2972   double dmant = num2dbl(mant, (char *)ARG1, s_make_dfloat);
 2973   int e = INUM(expt);
 2974   ASRTER(INUMP(expt), expt, ARG2, s_make_dfloat);
 2975   ASRTER((dmant < 0 ? -dmant : dmant) <= max_dbl_int, mant,
 2976      OUTOFRANGE, s_make_dfloat);
 2977   return makdbl(ldexp(dmant, e), 0.0);
 2978 }
 2979 # endif
 2980 #endif
 2981 
 2982 unsigned long hasher(obj, n, d)
 2983      SCM obj;
 2984      unsigned long n;
 2985      sizet d;
 2986 {
 2987   switch (7 & PTR2INT(obj)) {
 2988   case 2: case 6:       /* INUMP(obj) */
 2989     return INUM(obj) % n;
 2990   case 4:
 2991     if (ICHRP(obj))
 2992       return (unsigned)(downcase[ICHR(obj)]) % n;
 2993     switch ((int) obj) {
 2994 #ifndef SICP
 2995     case (int) EOL: d = 256; break;
 2996 #endif
 2997     case (int) BOOL_T: d = 257; break;
 2998     case (int) BOOL_F: d = 258; break;
 2999     case (int) EOF_VAL: d = 259; break;
 3000     default: d = 263;       /* perhaps should be error */
 3001     }
 3002     return d % n;
 3003   default: return 263 % n;  /* perhaps should be error */
 3004   case 0:
 3005     switch TYP7(obj) {
 3006     default: return 263 % n;
 3007     case tc7_smob:
 3008       switch TYP16(obj) {
 3009       case tcs_bignums:
 3010       bighash: return INUM(modulo(obj, MAKINUM(n)));
 3011       default: return 263 % n;
 3012 #ifdef FLOATS
 3013       case tc16_flo:
 3014     if (REALP(obj)) {
 3015       double r = REALPART(obj);
 3016       if (floor(r)==r) {
 3017         obj = in2ex(obj);
 3018         if (IMP(obj)) return INUM(obj) % n;
 3019         goto bighash;
 3020       }
 3021     }
 3022     obj = number2string(obj, MAKINUM(10));
 3023 #endif
 3024       }
 3025     case tcs_symbols: case tc7_string:
 3026       return strhash(UCHARS(obj), (sizet) LENGTH(obj), n);
 3027     case tc7_vector: {
 3028       sizet len = LENGTH(obj);
 3029       SCM *data = VELTS(obj);
 3030       if (len>5) {
 3031     sizet i = d/2;
 3032     unsigned long h = 1;
 3033     while (i--) h = ((h<<8) + (hasher(data[h % len], n, 2))) % n;
 3034     return h;
 3035       }
 3036       else {
 3037     sizet i = len;
 3038     unsigned long h = (n)-1;
 3039     while (i--) h = ((h<<8) + (hasher(data[i], n, d/len))) % n;
 3040     return h;
 3041       }
 3042     }
 3043     case tcs_cons_imcar: case tcs_cons_nimcar:
 3044       if (d) return (hasher(CAR(obj), n, d/2)+hasher(CDR(obj), n, d/2)) % n;
 3045       else return 1;
 3046     case tc7_port:
 3047       return ((RDNG & CAR(obj)) ? 260 : 261) % n;
 3048     case tcs_closures: case tc7_contin: case tcs_subrs:
 3049       return 262 % n;
 3050     }
 3051   }
 3052 }
 3053 
 3054 static char s_hashv[] = "hashv", s_hashq[] = "hashq";
 3055 extern char s_obunhash[];
 3056 #define s_hash (&s_obunhash[9])
 3057 
 3058 SCM hash(obj, n)
 3059      SCM obj;
 3060      SCM n;
 3061 {
 3062   ASRTER(INUMP(n) && 0 <= n, n, ARG2, s_hash);
 3063   return MAKINUM(hasher(obj, INUM(n), 10));
 3064 }
 3065 
 3066 SCM hashv(obj, n)
 3067      SCM obj;
 3068      SCM n;
 3069 {
 3070   ASRTER(INUMP(n) && 0 <= n, n, ARG2, s_hashv);
 3071   if (ICHRP(obj)) return MAKINUM((unsigned)(downcase[ICHR(obj)]) % INUM(n));
 3072   if (NIMP(obj) && NUMP(obj)) return MAKINUM(hasher(obj, INUM(n), 10));
 3073   else return MAKINUM(obj % INUM(n));
 3074 }
 3075 
 3076 SCM hashq(obj, n)
 3077      SCM obj;
 3078      SCM n;
 3079 {
 3080   ASRTER(INUMP(n) && 0 <= n, n, ARG2, s_hashq);
 3081   return MAKINUM((((unsigned) obj) >> 1) % INUM(n));
 3082 }
 3083 
 3084 static iproc subr1s[] = {
 3085     {"number?", numberp},
 3086     {s_inexactp, inexactp},
 3087 #ifdef FLOATS
 3088     {"complex?", scm_complex_p},
 3089     {"real?", realp},
 3090     {"rational?", scm_rationalp},
 3091     {"integer?", intp},
 3092     {s_real_part, real_part},
 3093     {s_imag_part, imag_part},
 3094     {s_magnitude, scm_magnitude},
 3095     {s_angle, angle},
 3096     {s_in2ex, in2ex},
 3097     {s_ex2in, ex2in},
 3098     {s_abs, scm_abs},
 3099 # ifdef BIGDIG
 3100     {s_dfloat_parts, scm_dfloat_parts},
 3101 # endif
 3102 #else
 3103     {"complex?", numberp},
 3104     {"real?", numberp},
 3105     {"rational?", numberp},
 3106     {"integer?", exactp},
 3107     {"floor", numident},
 3108     {"ceiling", numident},
 3109     {s_trunc, numident},
 3110     {"round", numident},
 3111     {s_abs, scm_iabs},
 3112 #endif
 3113     {s_zerop, zerop},
 3114     {s_positivep, positivep},
 3115     {s_negativep, negativep},
 3116     {s_str2list, string2list},
 3117     {"list->string", string},
 3118     {s_st_copy, string_copy},
 3119     {"list->vector", vector},
 3120     {s_vect2list, vector2list},
 3121     {0, 0}};
 3122 
 3123 static iproc asubrs[] = {
 3124     {s_difference, difference},
 3125     {s_divide, divide},
 3126     {s_max, scm_max},
 3127     {s_min, scm_min},
 3128     {s_sum, sum},
 3129     {s_product, product},
 3130     {0, 0}};
 3131 
 3132 static iproc subr2s[] = {
 3133 #ifdef FLOATS
 3134     {s_makrect, makrect},
 3135     {s_makpolar, makpolar},
 3136     {s_atan2, latan2},
 3137     {s_expt, expt},
 3138 # ifdef BIGDIG
 3139     {s_make_dfloat, scm_make_dfloat},
 3140 # endif
 3141 #endif
 3142 #ifdef INUMS_ONLY
 3143     {s_memv, memq},
 3144     {s_assv, assq},
 3145 #else
 3146     {s_memv, memv},
 3147     {s_assv, assv},
 3148 #endif
 3149     {s_intexpt, scm_intexpt},
 3150     {s_intlog, scm_intlog},
 3151     {s_cintlog, scm_cintlog},
 3152     {s_list_tail, list_tail},
 3153     {s_ve_fill, vector_fill},
 3154     {s_st_fill, string_fill},
 3155     {s_hash, hash},
 3156     {s_hashv, hashv},
 3157     {s_hashq, hashq},
 3158     {0, 0}};
 3159 
 3160 static iproc subr2os[] = {
 3161     {s_str2number, string2number},
 3162     {s_number2string, number2string},
 3163     {0, 0}};
 3164 
 3165 static iproc rpsubrs[] = {
 3166 #ifdef INUMS_ONLY
 3167     {"eqv?", eq},
 3168 #else
 3169     {"eqv?", eqv},
 3170 #endif
 3171     {s_eqp, eqp},
 3172     {s_lessp, lessp},
 3173     {s_grp, greaterp},
 3174     {s_leqp, leqp},
 3175     {s_greqp, greqp},
 3176     {0, 0}};
 3177 
 3178 #ifdef FLOATS
 3179 static dblproc cxrs[] = {
 3180     {"floor", floor},
 3181     {"ceiling", ceil},
 3182     {"truncate", scm_truncate},
 3183     {"round", scm_round},
 3184     {"$abs", fabs},
 3185     {"real-sqrt", sqrt},
 3186     {"real-exp", exp},
 3187     {"real-ln", log},
 3188     {"real-log10", log10},
 3189     {"real-sin", sin},
 3190     {"real-cos", cos},
 3191     {"real-tan", tan},
 3192     {"real-asin", asin},
 3193     {"real-acos", acos},
 3194     {"real-atan", atan},
 3195     {"real-sinh", sinh},
 3196     {"real-cosh", cosh},
 3197     {"real-tanh", tanh},
 3198     {"real-asinh", asinh},
 3199     {"real-acosh", acosh},
 3200     {"real-atanh", atanh},
 3201     {0, 0}};
 3202 #endif
 3203 
 3204 #ifdef FLOATS
 3205 static void safe_add_1(f, fsum)
 3206      double f, *fsum;
 3207 {
 3208   *fsum = f + 1.0;
 3209 }
 3210 #endif
 3211 
 3212 void init_scl()
 3213 {
 3214   init_iprocs(subr1s, tc7_subr_1);
 3215   init_iprocs(subr2os, tc7_subr_2o);
 3216   init_iprocs(subr2s, tc7_subr_2);
 3217   init_iprocs(asubrs, tc7_asubr);
 3218   init_iprocs(rpsubrs, tc7_rpsubr);
 3219 #ifdef SICP
 3220   add_feature("sicp");
 3221 #endif
 3222 #ifdef FLOATS
 3223   init_iprocs((iproc *)cxrs, tc7_cxr);
 3224 # ifdef SINGLES
 3225   NEWCELL(flo0);
 3226   CAR(flo0) = tc_flo;
 3227   FLO(flo0) = 0.0;
 3228 # else
 3229   DEFER_INTS;
 3230   flo0 = must_malloc_cell(1L*sizeof(double), (SCM)tc_dblr, "real");
 3231   REAL(flo0) = 0.0;
 3232   ALLOW_INTS;
 3233 # endif
 3234 # ifndef _MSC_VER
 3235   DEFER_INTS;
 3236   scm_narn = must_malloc_cell(2L*sizeof(double), (SCM)tc_dblc, "complex");
 3237   REAL(scm_narn) = 0.0/0.0;
 3238   IMAG(scm_narn) = 0.0/0.0;
 3239   ALLOW_INTS;
 3240 # endif
 3241 # ifndef BIGDIG
 3242 #  ifdef DBL_DIG
 3243   dblprec = (DBL_DIG > 20) ? 20 : DBL_DIG;
 3244 #  else
 3245   {             /* determine floating point precision */
 3246     double f = 0.1;
 3247     volatile double fsum = 1.0+f;
 3248     while (fsum != 1.0) {
 3249       f /= 10.0;
 3250       if (++dblprec > 20) break;
 3251       safe_add_1(f, &fsum);
 3252     }
 3253     dblprec = dblprec-1;
 3254   }
 3255 #  endif /* DBL_DIG */
 3256 # else   /* !BIGDIG */
 3257   {
 3258     int idx;
 3259     dpows5[0] = 1.0;
 3260     for (idx = 1; idx < 23; idx++) {
 3261       dpows5[idx] = 5*dpows5[idx-1];
 3262     }
 3263   }
 3264 # endif /* !BIGDIG */
 3265 # ifdef DBL_MANT_DIG
 3266   dbl_mant_dig = DBL_MANT_DIG;
 3267 # else
 3268   {             /* means we #defined it. */
 3269     volatile double fsum = 0.0;
 3270     double eps = 1.0;
 3271     int i = 0;
 3272     while (fsum != 1.0) {
 3273       eps = 0.5 * eps;
 3274       safe_add_1(eps, &fsum);
 3275       i++;
 3276     }
 3277     dbl_mant_dig = i;
 3278   }
 3279 # endif /* DBL_MANT_DIG */
 3280   max_dbl_int = pow(2.0, dbl_mant_dig - 1.0);
 3281   max_dbl_int = max_dbl_int + (max_dbl_int - 1.0);
 3282   /* dbl_eps = ldexp(1.0, - dbl_mant_dig); */
 3283   sysintern("double-float-mantissa-length", MAKINUM(dbl_mant_dig));
 3284   sysintern("bigdbl:powers-of-5",
 3285         pows5 = make_vector(MAKINUM(326), MAKINUM(1)));
 3286 #endif
 3287 }