root/src/number.c

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
  1. bad_number_method
  2. Scm_MakeFlonum
  3. Scm_MakeFlonumToNumber
  4. Scm_DecodeFlonum
  5. Scm_MakeComplex
  6. Scm_MakeComplexNormalized
  7. Scm_MakeComplexPolar
  8. Scm_Magnitude
  9. Scm_Angle
  10. Scm_MakeInteger
  11. Scm_MakeIntegerU
  12. Scm_GetIntegerClamp
  13. Scm_GetIntegerUClamp
  14. Scm_GetInteger32Clamp
  15. Scm_GetIntegerU32Clamp
  16. Scm_MakeInteger64
  17. Scm_MakeIntegerU64
  18. Scm_GetInteger64Clamp
  19. Scm_GetIntegerU64Clamp
  20. Scm_GetDouble
  21. Scm_IntegerP
  22. Scm_OddP
  23. Scm_Abs
  24. Scm_Sign
  25. Scm_Negate
  26. Scm_Reciprocal
  27. Scm_ExactToInexact
  28. Scm_InexactToExact
  29. Scm_PromoteToBignum
  30. Scm_PromoteToFlonum
  31. Scm_PromoteToComplex
  32. Scm_Add
  33. Scm_Subtract
  34. Scm_Multiply
  35. Scm_Divide
  36. Scm_Quotient
  37. Scm_Modulo
  38. iexpt10_init
  39. exact_expt
  40. Scm_Expt
  41. Scm_NumEq
  42. Scm_NumCmp
  43. Scm_MinMax
  44. roundeven
  45. Scm_Round
  46. Scm_Ash
  47. Scm_LogNot
  48. Scm_LogAnd
  49. Scm_LogIor
  50. Scm_LogXor
  51. iexpt10
  52. ipow
  53. raise_pow10
  54. numcmp3
  55. double_print
  56. number_print
  57. Scm_NumberToString
  58. Scm_PrintDouble
  59. read_uint
  60. algorithmR
  61. read_real
  62. read_number
  63. numread_error
  64. Scm_StringToNumber
  65. Scm__InitNumber

   1 /*
   2  * number.c - numeric functions
   3  *
   4  *   Copyright (c) 2000-2005 Shiro Kawai, All rights reserved.
   5  * 
   6  *   Redistribution and use in source and binary forms, with or without
   7  *   modification, are permitted provided that the following conditions
   8  *   are met:
   9  * 
  10  *   1. Redistributions of source code must retain the above copyright
  11  *      notice, this list of conditions and the following disclaimer.
  12  *
  13  *   2. Redistributions in binary form must reproduce the above copyright
  14  *      notice, this list of conditions and the following disclaimer in the
  15  *      documentation and/or other materials provided with the distribution.
  16  *
  17  *   3. Neither the name of the authors nor the names of its contributors
  18  *      may be used to endorse or promote products derived from this
  19  *      software without specific prior written permission.
  20  *
  21  *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  22  *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  23  *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  24  *   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  25  *   OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  26  *   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
  27  *   TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  28  *   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  29  *   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  30  *   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  31  *   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  32  *
  33  *  $Id: number.c,v 1.123 2005/10/13 08:14:13 shirok Exp $
  34  */
  35 
  36 #include <math.h>
  37 #include <limits.h>
  38 #include <string.h>
  39 #include <ctype.h>
  40 #include <float.h>
  41 #define LIBGAUCHE_BODY
  42 #include "gauche.h"
  43 #include "gauche/scmconst.h"
  44 
  45 #ifndef M_PI
  46 #define M_PI 3.14159265358979323846
  47 #endif
  48 
  49 #ifdef HAVE_ISNAN
  50 #define SCM_IS_NAN(x)  isnan(x)
  51 #else
  52 #define SCM_IS_NAN(x)  FALSE    /* we don't have a clue */
  53 #endif
  54 
  55 #ifdef HAVE_ISINF
  56 #define SCM_IS_INF(x)  isinf(x)
  57 #else
  58 #define SCM_IS_INF(x)  ((x) != 0 && (x) == (x)/2.0)
  59 #endif
  60 
  61 #define RADIX_MIN 2
  62 #define RADIX_MAX 36
  63 
  64 /* Linux gcc have those, but the declarations aren't included unless
  65    __USE_ISOC9X is defined.  Just in case. */
  66 #ifdef HAVE_TRUNC
  67 extern double trunc(double);
  68 #endif
  69 
  70 #ifdef HAVE_RINT
  71 extern double rint(double);
  72 #define roundeven rint
  73 #else
  74 static double roundeven(double);
  75 #endif
  76 
  77 /*
  78  * Classes of Numeric Tower
  79  */
  80 
  81 static ScmClass *numeric_cpl[] = {
  82     SCM_CLASS_STATIC_PTR(Scm_RealClass),
  83     SCM_CLASS_STATIC_PTR(Scm_ComplexClass),
  84     SCM_CLASS_STATIC_PTR(Scm_NumberClass),
  85     SCM_CLASS_STATIC_PTR(Scm_TopClass),
  86     NULL
  87 };
  88 
  89 static void number_print(ScmObj obj, ScmPort *port, ScmWriteContext *ctx);
  90 
  91 SCM_DEFINE_BUILTIN_CLASS(Scm_NumberClass, number_print, NULL, NULL, NULL,
  92                          numeric_cpl+3);
  93 SCM_DEFINE_BUILTIN_CLASS(Scm_ComplexClass, number_print, NULL, NULL, NULL,
  94                          numeric_cpl+2);
  95 SCM_DEFINE_BUILTIN_CLASS(Scm_RealClass, number_print, NULL, NULL, NULL,
  96                          numeric_cpl+1);
  97 SCM_DEFINE_BUILTIN_CLASS(Scm_IntegerClass, number_print, NULL, NULL, NULL,
  98                          numeric_cpl);
  99 
 100 /*=====================================================================
 101  *  Generic Arithmetic
 102  */
 103 
 104 /* Some arithmetic operations calls the corresponding generic function
 105  * if the operand is not a number.
 106  */
 107 
 108 /* Fallback Gf */
 109 static ScmObj bad_number_method(ScmObj *args, int nargs, ScmGeneric *gf)
 110 {
 111     const char *fn = (const char *)SCM_GENERIC_DATA(gf);
 112     if (nargs == 1) {
 113         Scm_Error("operation %s is not defined on object %S", fn, args[0]);
 114     } else if (nargs == 2) {
 115         Scm_Error("operation %s is not defined between %S and %S",
 116                   fn, args[0], args[1]);
 117     } else {
 118         Scm_Error("generic function for %s is called with args %S",
 119                   fn, Scm_ArrayToList(args, nargs));
 120     }
 121     return SCM_UNDEFINED;
 122 }
 123 static SCM_DEFINE_GENERIC(generic_add, bad_number_method, "+");
 124 static SCM_DEFINE_GENERIC(generic_sub, bad_number_method, "-");
 125 static SCM_DEFINE_GENERIC(generic_mul, bad_number_method, "*");
 126 static SCM_DEFINE_GENERIC(generic_div, bad_number_method, "/");
 127 
 128 /*=====================================================================
 129  *  Flonums
 130  */
 131 
 132 ScmObj Scm_MakeFlonum(double d)
 133 {
 134     ScmFlonum *f = SCM_NEW(ScmFlonum);
 135     SCM_SET_CLASS(f, SCM_CLASS_REAL);
 136     f->value = d;
 137     return SCM_OBJ(f);
 138 }
 139 
 140 ScmObj Scm_MakeFlonumToNumber(double d, int exact)
 141 {
 142     if (exact && !SCM_IS_INF(d)) {
 143         /* see if d can be demoted to integer */
 144         double i, f;
 145         f = modf(d, &i);
 146         if (f == 0.0) {
 147             if (i > SCM_SMALL_INT_MAX || i < SCM_SMALL_INT_MIN) {
 148                 return Scm_MakeBignumFromDouble(i);
 149             } else {
 150                 return SCM_MAKE_INT((long)i);
 151             }
 152         }
 153     }
 154     return Scm_MakeFlonum(d);
 155 }
 156 
 157 /* Decompose flonum D into an integer mantissa F and exponent E, where
 158  *   -1022 <= E <= 1023,
 159  *    0 <= abs(F) < 2^53
 160  *    D = F * 2^(E - 53)
 161  * Some special cases:
 162  *    F = 0, E = 0 if D = 0.0 or -0.0
 163  *    F = #t if D is infinity (positive or negative)
 164  *    F = #f if D is NaN.
 165  * If D is normalized number, F >= 2^52.
 166  *
 167  * Cf. IEEE 754 Reference
 168  * http://babbage.cs.qc.edu/courses/cs341/IEEE-754references.html
 169  */
 170 union ieee_double {
 171     double d;
 172     struct {
 173 #ifdef DOUBLE_ARMENDIAN
 174         /* ARM's mixed endian.  TODO: what if we have LP64 ARM? */
 175         unsigned long mant0:20;
 176         unsigned int exp:11;
 177         unsigned int sign:1;
 178         unsigned long mant1:32;
 179 #else  /*!DOUBLE_ARMENDIAN*/
 180 #ifdef WORDS_BIGENDIAN
 181 #if SIZEOF_LONG >= 8
 182         unsigned int sign:1;
 183         unsigned int exp:11;
 184         unsigned long mant:52;
 185 #else  /*SIZEOF_LONG < 8*/
 186         unsigned int sign:1;
 187         unsigned int exp:11;
 188         unsigned long mant0:20;
 189         unsigned long mant1:32;
 190 #endif /*SIZEOF_LONG < 8*/
 191 #else  /*!WORDS_BIGENDIAN*/
 192 #if SIZEOF_LONG >= 8
 193         unsigned long mant:52;
 194         unsigned int  exp:11;
 195         unsigned int  sign:1;
 196 #else  /*SIZEOF_LONG < 8*/
 197         unsigned long mant1:32;
 198         unsigned long mant0:20;
 199         unsigned int  exp:11;
 200         unsigned int  sign:1;
 201 #endif /*SIZEOF_LONG < 8*/
 202 #endif /*!WORDS_BIGENDIAN*/
 203 #endif /*!DOUBLE_ARMENDIAN*/
 204     } components;
 205 };
 206 
 207 ScmObj Scm_DecodeFlonum(double d, int *exp, int *sign)
 208 {
 209     union ieee_double dd;
 210     ScmObj f;
 211     
 212     dd.d = d;
 213 
 214     /* Check exceptional cases */
 215     if (dd.components.exp == 0x7ff) {
 216         *exp = 0;
 217         if (
 218 #if SIZEOF_LONG >= 8
 219             dd.components.mant == 0
 220 #else  /*SIZEOF_LONG < 8*/
 221             dd.components.mant0 == 0 && dd.components.mant1 == 0
 222 #endif /*SIZEOF_LONG < 8*/
 223             ) {
 224             return SCM_TRUE;  /* infinity */
 225         } else {
 226             return SCM_FALSE; /* NaN */
 227         }
 228     }
 229 
 230     *exp  = (dd.components.exp? dd.components.exp - 0x3ff - 52 : -0x3fe - 52);
 231     *sign = (dd.components.sign? -1 : 1);
 232     
 233 #if SIZEOF_LONG >= 8
 234     {
 235         unsigned long lf = dd.components.mant;
 236         if (dd.components.exp > 0) {
 237             lf += (1L<<52);     /* hidden bit */
 238         }
 239         f = Scm_MakeInteger(lf);
 240     }
 241 #else  /*SIZEOF_LONG < 8*/
 242     {
 243         unsigned long values[2];
 244         values[0] = dd.components.mant1;
 245         values[1] = dd.components.mant0;
 246         if (dd.components.exp > 0) {
 247             values[1] += (1L<<20); /* hidden bit */
 248         }
 249         f = Scm_NormalizeBignum(SCM_BIGNUM(Scm_MakeBignumFromUIArray(1, values, 2)));
 250     }
 251 #endif /*SIZEOF_LONG < 8*/
 252     return f;
 253 }
 254 
 255 /*=======================================================================
 256  *  Complex numbers
 257  */
 258 
 259 ScmObj Scm_MakeComplex(double r, double i)
 260 {
 261     ScmComplex *c = SCM_NEW_ATOMIC(ScmComplex);
 262     SCM_SET_CLASS(c, SCM_CLASS_COMPLEX);
 263     c->real = r;
 264     c->imag = i;
 265     return SCM_OBJ(c);
 266 }
 267 
 268 ScmObj Scm_MakeComplexNormalized(double r, double i)
 269 {
 270     if (i == 0.0) return Scm_MakeFlonum(r);
 271     else          return Scm_MakeComplex(r, i);
 272 }
 273 
 274 ScmObj Scm_MakeComplexPolar(double mag, double angle)
 275 {
 276     double real = mag * cos(angle);
 277     double imag = mag * sin(angle);
 278     if (imag == 0.0) return Scm_MakeFlonum(real);
 279     else             return Scm_MakeComplex(real, imag);
 280 }
 281 
 282 ScmObj Scm_Magnitude(ScmObj z)
 283 {
 284     double m;
 285     if (SCM_REALP(z)) {
 286         m = fabs(Scm_GetDouble(z));
 287     } else if (!SCM_COMPLEXP(z)) {
 288         Scm_Error("number required, but got %S", z);
 289         m = 0.0;                /* dummy */
 290     } else {
 291         double r = SCM_COMPLEX_REAL(z);
 292         double i = SCM_COMPLEX_IMAG(z);
 293         m = sqrt(r*r+i*i);
 294     }
 295     return Scm_MakeFlonum(m);
 296 }
 297 
 298 ScmObj Scm_Angle(ScmObj z)
 299 {
 300     double a;
 301     if (SCM_REALP(z)) {
 302         a = (Scm_Sign(z) < 0)? M_PI : 0.0;
 303     } else if (!SCM_COMPLEXP(z)) {
 304         Scm_Error("number required, but got %S", z);
 305         a = 0.0;                /* dummy */
 306     } else {
 307         double r = SCM_COMPLEX_REAL(z);
 308         double i = SCM_COMPLEX_IMAG(z);
 309         a = atan2(i, r);
 310     }
 311     return Scm_MakeFlonum(a);
 312 }
 313 
 314 /*=======================================================================
 315  *  Coertion
 316  */
 317 
 318 ScmObj Scm_MakeInteger(long i)
 319 {
 320     if (i >= SCM_SMALL_INT_MIN && i <= SCM_SMALL_INT_MAX) {
 321         return SCM_MAKE_INT(i);
 322     } else {
 323         return Scm_MakeBignumFromSI(i);
 324     }
 325 }
 326 
 327 ScmObj Scm_MakeIntegerU(u_long i)
 328 {
 329     if (i <= (u_long)SCM_SMALL_INT_MAX) return SCM_MAKE_INT(i);
 330     else return Scm_MakeBignumFromUI(i);
 331 }
 332 
 333 /* Convert scheme integer to C integer */
 334 long Scm_GetIntegerClamp(ScmObj obj, int clamp, int *oor)
 335 {
 336     if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
 337     if (SCM_INTP(obj)) return SCM_INT_VALUE(obj);
 338     else if (SCM_BIGNUMP(obj)) {
 339         return Scm_BignumToSI(SCM_BIGNUM(obj), clamp, oor);
 340     }
 341     else if (SCM_FLONUMP(obj)) {
 342         double v = SCM_FLONUM_VALUE(obj);
 343         if (v > (double)LONG_MAX) {
 344             if (clamp & SCM_CLAMP_HI) return LONG_MAX;
 345             else goto err;
 346         }
 347         if (v < (double)LONG_MIN) {
 348             if (clamp & SCM_CLAMP_LO) return LONG_MIN;
 349             else goto err;
 350         }
 351         return (long)v;
 352     }
 353   err:
 354     if (clamp == SCM_CLAMP_NONE && oor != NULL) {
 355         *oor = TRUE;
 356     } else {
 357         Scm_Error("argument out of range: %S", obj);
 358     }
 359     return 0;
 360 }
 361 
 362 u_long Scm_GetIntegerUClamp(ScmObj obj, int clamp, int *oor)
 363 {
 364     if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
 365     if (SCM_INTP(obj)) {
 366         if (SCM_INT_VALUE(obj) < 0) {
 367             if (clamp & SCM_CLAMP_LO) return 0;
 368             else goto err;
 369         }
 370         return SCM_INT_VALUE(obj);
 371     } else if (SCM_BIGNUMP(obj)) {
 372         return Scm_BignumToUI(SCM_BIGNUM(obj), clamp, oor);
 373     }
 374     else if (SCM_FLONUMP(obj)) {
 375         double v = SCM_FLONUM_VALUE(obj);
 376         if (v > (double)ULONG_MAX) {
 377             if (clamp & SCM_CLAMP_HI) return ULONG_MAX;
 378             else goto err;
 379         }
 380         if (v < 0.0) {
 381             if (clamp & SCM_CLAMP_LO) return 0;
 382             else goto err;
 383         }
 384         return (u_long)v;
 385     }
 386   err:
 387     if (clamp == SCM_CLAMP_NONE && oor != NULL) {
 388         *oor = TRUE;
 389     } else {
 390         Scm_Error("argument out of range: %S", obj);
 391     }
 392     return 0;
 393 }
 394 
 395 /* 32bit integer specific */
 396 ScmInt32 Scm_GetInteger32Clamp(ScmObj obj, int clamp, int *oor)
 397 {
 398 #if SIZEOF_LONG == 4
 399     return (ScmInt32)Scm_GetIntegerClamp(obj, clamp, oor);
 400 #else  /* SIZEOF_LONG >= 8 */
 401     if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
 402     /* NB: we denote the constant directly here.  (1L<<31) fails on
 403        Alpha machines, since the compiler somehow calculates the constant
 404        in 32bit integer even it has 'L'.  We have to write (1LL<<31), but
 405        I'm afraid that it's not portable. */
 406     if (SCM_INTP(obj)) {
 407         long r = SCM_INT_VALUE(obj);
 408         if (r < -0x80000000L) {
 409             if (clamp & SCM_CLAMP_LO) return -0x80000000L;
 410             goto err;
 411         }
 412         if (r > 0x7fffffffL) {
 413             if (clamp & SCM_CLAMP_HI) return 0x7fffffffL;
 414             goto err;
 415         }
 416         return r;
 417     } else if (SCM_BIGNUMP(obj)) {
 418         if (SCM_BIGNUM_SIGN(obj) < 0) {
 419             if (clamp & SCM_CLAMP_LO) return -0x80000000L;
 420             goto err;
 421         } else {
 422             if (clamp & SCM_CLAMP_HI) return 0x7fffffffL;
 423             goto err;
 424         }
 425     }
 426   err:
 427     if (clamp == SCM_CLAMP_NONE && oor != NULL) {
 428         *oor = TRUE;
 429     } else {
 430         Scm_Error("argument out of range: %S", obj);
 431     }
 432     return 0;
 433 #endif /* SIZEOF_LONG >= 8 */
 434 }
 435 
 436 ScmUInt32 Scm_GetIntegerU32Clamp(ScmObj obj, int clamp, int *oor)
 437 {
 438 #if SIZEOF_LONG == 4
 439     return (ScmUInt32)Scm_GetIntegerUClamp(obj, clamp, oor);
 440 #else  /* SIZEOF_LONG >= 8 */
 441     if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
 442     if (SCM_INTP(obj)) {
 443         long r = SCM_INT_VALUE(obj);
 444         if (r < 0) {
 445             if (clamp & SCM_CLAMP_LO) return 0;
 446             goto err;
 447         }
 448         if (r > 0xffffffffUL) {
 449             if (clamp & SCM_CLAMP_HI) return 0xffffffffUL;
 450             goto err;
 451         }
 452         return r;
 453     } else if (SCM_BIGNUMP(obj)) {
 454         if (SCM_BIGNUM_SIGN(obj) < 0) {
 455             if (clamp & SCM_CLAMP_LO) return 0;
 456             goto err;
 457         } else {
 458             if (clamp & SCM_CLAMP_HI) return 0xffffffffUL;
 459             goto err;
 460         }
 461     }
 462   err:
 463     if (clamp == SCM_CLAMP_NONE && oor != NULL) {
 464         *oor = TRUE;
 465     } else {
 466         Scm_Error("argument out of range: %S", obj);
 467     }
 468     return 0;
 469 #endif /* SIZEOF_LONG >= 8 */
 470 }
 471 
 472 
 473 #if SIZEOF_LONG == 4
 474 /* we need special routines */
 475 ScmObj Scm_MakeInteger64(ScmInt64 i)
 476 {
 477 #if SCM_EMULATE_INT64
 478     u_long val[2];
 479     if (i.hi == 0) return Scm_MakeInteger(i.lo);
 480     val[0] = i.lo;
 481     val[1] = i.hi;
 482     return Scm_MakeBignumFromUIArray(0, val, 2); /* bignum checks sign */
 483 #else /*SCM_EMULATE_INT64*/
 484     u_long val[2];
 485     val[0] = (uint64_t)i & ULONG_MAX;
 486     val[1] = (uint64_t)i >> 32;
 487     if (val[1] == 0 && val[0] <= LONG_MAX) return Scm_MakeInteger(val[0]);
 488     return Scm_NormalizeBignum(SCM_BIGNUM(Scm_MakeBignumFromUIArray(0, val, 2)));
 489 #endif
 490 }
 491 
 492 ScmObj Scm_MakeIntegerU64(ScmUInt64 i)
 493 {
 494 #if SCM_EMULATE_INT64
 495     u_long val[2];
 496     if (i.hi == 0) return Scm_MakeIntegerU(i.lo);
 497     val[0] = i.lo;
 498     val[1] = i.hi;
 499     return Scm_MakeBignumFromUIArray(1, val, 2);
 500 #else /*SCM_EMULATE_INT64*/
 501     u_long val[2];
 502     val[0] = (uint64_t)i & ULONG_MAX;
 503     val[1] = (uint64_t)i >> 32;
 504     if (val[1] == 0) return Scm_MakeIntegerU(val[0]);
 505     return Scm_MakeBignumFromUIArray(1, val, 2);
 506 #endif
 507 }
 508 
 509 ScmInt64 Scm_GetInteger64Clamp(ScmObj obj, int clamp, int *oor)
 510 {
 511 #if SCM_EMULATE_INT64
 512     ScmInt64 r = {0, 0};
 513     if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
 514     if (SCM_INTP(obj)) {
 515         long v = SCM_INT_VALUE(obj);
 516         r.lo = v;
 517         if (v < 0) r.hi = ULONG_MAX;
 518         return r;
 519     }
 520     if (SCM_BIGNUMP(obj)) {
 521         return Scm_BignumToSI64(SCM_BIGNUM(obj), clamp, oor);
 522     }
 523     if (SCM_FLONUMP(obj)) {
 524         if (Scm_NumCmp(obj, SCM_2_63) >= 0) {
 525             if (!(clamp&SCM_CLAMP_HI)) goto err;
 526             SCM_SET_INT64_MAX(r);
 527             return r;
 528         } else if (Scm_NumCmp(obj, SCM_MINUS_2_63) < 0) {
 529             if (!(clamp&SCM_CLAMP_LO)) goto err;
 530             SCM_SET_INT64_MIN(r);
 531             return r;
 532         } else {
 533             ScmObj b = Scm_MakeBignumFromDouble(SCM_FLONUM_VALUE(obj));
 534             return Scm_BignumToSI64(SCM_BIGNUM(b), clamp, oor);
 535         }
 536     }
 537 #else /*!SCM_EMULATE_INT64*/
 538     ScmInt64 r = 0;
 539     if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
 540     if (SCM_INTP(obj)) return (ScmInt64)SCM_INT_VALUE(obj);
 541     if (SCM_BIGNUMP(obj)) {
 542         return Scm_BignumToSI64(SCM_BIGNUM(obj), clamp, oor);
 543     }
 544     if (SCM_FLONUMP(obj)) {
 545         int64_t maxval, minval;
 546         double v;
 547         
 548         SCM_SET_INT64_MAX(maxval);
 549         SCM_SET_INT64_MIN(minval);
 550         v = SCM_FLONUM_VALUE(obj);
 551         if (v > (double)maxval) {
 552             if (!(clamp&SCM_CLAMP_HI)) goto err;
 553             return maxval;
 554         } else if (v < (double)minval) {
 555             if (!(clamp&SCM_CLAMP_LO)) goto err;
 556             return minval;
 557         } else {
 558             return (long)v;
 559         }
 560     }
 561 #endif /*!SCM_EMULATE_INT64*/
 562   err:
 563     if (clamp == SCM_CLAMP_NONE && oor != NULL) {
 564         *oor = TRUE;
 565     } else {
 566         Scm_Error("argument out of range: %S", obj);
 567     }
 568     return r;
 569 }
 570                                
 571 ScmUInt64 Scm_GetIntegerU64Clamp(ScmObj obj, int clamp, int *oor)
 572 {
 573 #if SCM_EMULATE_INT64
 574     ScmUInt64 r = {0, 0};
 575     if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
 576     if (SCM_INTP(obj)) {
 577         long v = SCM_INT_VALUE(obj);
 578         if (v < 0) {
 579             if (!(clamp&SCM_CLAMP_LO)) goto err;
 580         } else {
 581             r.lo = v;
 582         }
 583         return r;
 584     }
 585     if (SCM_BIGNUMP(obj)) {
 586         return Scm_BignumToUI64(SCM_BIGNUM(obj), clamp, oor);
 587     }
 588     if (SCM_FLONUMP(obj)) {
 589         if (Scm_NumCmp(obj, SCM_2_64) >= 0) {
 590             if (!(clamp&SCM_CLAMP_HI)) goto err;
 591             SCM_SET_UINT64_MAX(r);
 592             return r;
 593         } else if (SCM_FLONUM_VALUE(obj) < 0) {
 594             if (!(clamp&SCM_CLAMP_LO)) goto err;
 595             return r;
 596         } else {
 597             ScmObj b = Scm_MakeBignumFromDouble(SCM_FLONUM_VALUE(obj));
 598             return Scm_BignumToUI64(SCM_BIGNUM(b), clamp, oor);
 599         }
 600     }
 601 #else /*!SCM_EMULATE_INT64*/
 602     ScmInt64 r = 0;
 603     if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
 604     if (SCM_INTP(obj)) {
 605         long v = SCM_INT_VALUE(obj);
 606         if (v < 0) {
 607             if (!(clamp&SCM_CLAMP_LO)) goto err;
 608             return 0;
 609         } else {
 610             return (ScmUInt64)v;
 611         }
 612     }
 613     if (SCM_BIGNUMP(obj)) {
 614         return Scm_BignumToUI64(SCM_BIGNUM(obj), clamp, oor);
 615     }
 616     if (SCM_FLONUMP(obj)) {
 617         double v = SCM_FLONUM_VALUE(obj);
 618         uint64_t maxval;
 619 
 620         if (v < 0) {
 621             if (!(clamp&SCM_CLAMP_LO)) goto err;
 622             return 0;
 623         }
 624         SCM_SET_UINT64_MAX(maxval);
 625         if (v > (double)maxval) {
 626             if (!(clamp&SCM_CLAMP_HI)) goto err;
 627             return maxval;
 628         } else {
 629             return (uint32_t)v;
 630         }
 631     }
 632 #endif
 633   err:
 634     if (clamp == SCM_CLAMP_NONE && oor != NULL) {
 635         *oor = TRUE;
 636     } else {
 637         Scm_Error("argument out of range: %S", obj);
 638     }
 639     return r;
 640 }
 641                                
 642 #endif /* SIZEOF_LONG == 4 */
 643 
 644 double Scm_GetDouble(ScmObj obj)
 645 {
 646     if (SCM_FLONUMP(obj)) return SCM_FLONUM_VALUE(obj);
 647     else if (SCM_INTP(obj)) return (double)SCM_INT_VALUE(obj);
 648     else if (SCM_BIGNUMP(obj)) return Scm_BignumToDouble(SCM_BIGNUM(obj));
 649     else return 0.0;
 650 }
 651 
 652 /*
 653  *   Generic Methods
 654  */
 655 
 656 /* Predicates */
 657 
 658 int Scm_IntegerP(ScmObj obj)
 659 {
 660     if (SCM_INTP(obj) || SCM_BIGNUMP(obj)) return TRUE;
 661     if (SCM_FLONUMP(obj)) {
 662         double d = SCM_FLONUM_VALUE(obj);
 663         double f, i;
 664         if ((f = modf(d, &i)) == 0.0) return TRUE;
 665         return FALSE;
 666     }
 667     if (SCM_COMPLEXP(obj)) return FALSE;
 668     Scm_Error("number required, but got %S", obj);
 669     return FALSE;           /* dummy */
 670 }
 671 
 672 int Scm_OddP(ScmObj obj)
 673 {
 674     if (SCM_INTP(obj)) {
 675         return (SCM_INT_VALUE(obj)&1);
 676     }
 677     if (SCM_BIGNUMP(obj)) {
 678         return (SCM_BIGNUM(obj)->values[0] & 1);
 679     }
 680     if (SCM_FLONUMP(obj) && Scm_IntegerP(obj)) {
 681         return (fmod(SCM_FLONUM_VALUE(obj), 2.0) != 0.0);
 682     }
 683     Scm_Error("integer required, but got %S", obj);
 684     return FALSE;       /* dummy */
 685     
 686 }
 687 
 688 /* Unary Operator */
 689 
 690 ScmObj Scm_Abs(ScmObj obj)
 691 {
 692     if (SCM_INTP(obj)) {
 693         long v = SCM_INT_VALUE(obj);
 694         if (v < 0) obj = SCM_MAKE_INT(-v);
 695     } else if (SCM_BIGNUMP(obj)) {
 696         if (SCM_BIGNUM_SIGN(obj) < 0) {
 697             obj = Scm_BignumCopy(SCM_BIGNUM(obj));
 698             SCM_BIGNUM_SIGN(obj) = 1;
 699         }
 700     } else if (SCM_FLONUMP(obj)) {
 701         double v = SCM_FLONUM_VALUE(obj);
 702         if (v < 0) obj = Scm_MakeFlonum(-v);
 703     } else if (SCM_COMPLEXP(obj)) {
 704         double r = SCM_COMPLEX_REAL(obj);
 705         double i = SCM_COMPLEX_IMAG(obj);
 706         double a = sqrt(r*r+i*i);
 707         return Scm_MakeFlonum(a);
 708     } else {
 709         Scm_Error("number required: %S", obj);
 710     }
 711     return obj;
 712 }
 713 
 714 /* Return -1, 0 or 1 when arg is minus, zero or plus, respectively.
 715    used to implement zero?, positive? and negative? */
 716 int Scm_Sign(ScmObj obj)
 717 {
 718     long r = 0;
 719     
 720     if (SCM_INTP(obj)) {
 721         r = SCM_INT_VALUE(obj);
 722         if (r > 0) r = 1;
 723         else if (r < 0) r = -1;
 724     } else if (SCM_BIGNUMP(obj)) {
 725         r = SCM_BIGNUM_SIGN(obj);
 726     } else if (SCM_FLONUMP(obj)) {
 727         double v = SCM_FLONUM_VALUE(obj);
 728         if (v != 0.0) {
 729             r = (v > 0.0)? 1 : -1;
 730         }
 731     } else {
 732         /* NB: zero? can accept a complex number, but it is processed in
 733            the stub function.   see stdlib.stub */
 734         Scm_Error("real number required, but got %S", obj);
 735     }
 736     return r;
 737 }
 738 
 739 ScmObj Scm_Negate(ScmObj obj)
 740 {
 741     if (SCM_INTP(obj)) {
 742         long v = SCM_INT_VALUE(obj);
 743         if (v == SCM_SMALL_INT_MIN) {
 744             obj = Scm_MakeBignumFromSI(-v);
 745         } else {
 746             obj = SCM_MAKE_INT(-v);
 747         }
 748     } else if (SCM_BIGNUMP(obj)) {
 749         obj = Scm_BignumNegate(SCM_BIGNUM(obj));
 750     } else if (SCM_FLONUMP(obj)) {
 751         obj = Scm_MakeFlonum(-SCM_FLONUM_VALUE(obj));
 752     } else if (SCM_COMPLEXP(obj)) {
 753         obj = Scm_MakeComplex(-SCM_COMPLEX_REAL(obj),
 754                               -SCM_COMPLEX_IMAG(obj));
 755     } else {
 756         obj = Scm_Apply(SCM_OBJ(&generic_sub), SCM_LIST1(obj));
 757     }
 758     return obj;
 759 }
 760 
 761 ScmObj Scm_Reciprocal(ScmObj obj)
 762 {
 763     if (SCM_INTP(obj)) {
 764         long val = SCM_INT_VALUE(obj);
 765         obj = Scm_MakeFlonum(1.0/(double)val);
 766     } else if (SCM_BIGNUMP(obj)) {
 767         double val = Scm_BignumToDouble(SCM_BIGNUM(obj));
 768         obj = Scm_MakeFlonum(1.0/val);
 769     } else if (SCM_FLONUMP(obj)) {
 770         double val = SCM_FLONUM_VALUE(obj);
 771         obj = Scm_MakeFlonum(1.0/val);
 772     } else if (SCM_COMPLEXP(obj)) {
 773         double r = SCM_COMPLEX_REAL(obj), r1;
 774         double i = SCM_COMPLEX_IMAG(obj), i1;
 775         double d;
 776         d = r*r + i*i;
 777         r1 = r/d;
 778         i1 = -i/d;
 779         obj = Scm_MakeComplexNormalized(r1, i1);
 780     } else {
 781         obj = Scm_Apply(SCM_OBJ(&generic_div), SCM_LIST1(obj));
 782     }
 783     return obj;
 784 }
 785 
 786 /*
 787  * Conversion operators
 788  */
 789 
 790 ScmObj Scm_ExactToInexact(ScmObj obj)
 791 {
 792     if (SCM_INTP(obj)) {
 793         obj = Scm_MakeFlonum((double)SCM_INT_VALUE(obj));
 794     } else if (SCM_BIGNUMP(obj)) {
 795         obj = Scm_MakeFlonum(Scm_BignumToDouble(SCM_BIGNUM(obj)));
 796     } else if (!SCM_FLONUMP(obj) && !SCM_COMPLEXP(obj)) {
 797         Scm_Error("number required: %S", obj);
 798     }
 799     return obj;
 800 }
 801 
 802 ScmObj Scm_InexactToExact(ScmObj obj)
 803 {
 804     if (SCM_FLONUMP(obj)) {
 805         double d = SCM_FLONUM_VALUE(obj);
 806         if (d < SCM_SMALL_INT_MIN || d > SCM_SMALL_INT_MAX) {
 807             obj = Scm_MakeBignumFromDouble(d);
 808         } else {
 809             obj = SCM_MAKE_INT((long)d);
 810         }
 811     } else if (SCM_COMPLEXP(obj)) {
 812         Scm_Error("exact complex is not supported: %S", obj);
 813     } if (!SCM_INTP(obj) && !SCM_BIGNUMP(obj)) {
 814         Scm_Error("number required: %S", obj);
 815     }
 816     return obj;
 817 }
 818 
 819 /* Type conversion:
 820  *   `promote' means a conversion from lower number class to higher,
 821  *      e.g. fixnum -> bignum -> flonum -> complex.
 822  *   `demote' means a conversion from higher number class to lower,
 823  *      e.g. complex -> flonum -> bignum -> fixnum.
 824  */
 825 
 826 ScmObj Scm_PromoteToBignum(ScmObj obj)
 827 {
 828     if (SCM_INTP(obj)) return Scm_MakeBignumFromSI(SCM_INT_VALUE(obj));
 829     if (SCM_BIGNUMP(obj)) return obj;
 830     Scm_Panic("Scm_PromoteToBignum: can't be here");
 831     return SCM_UNDEFINED;       /* dummy */
 832 }
 833 
 834 ScmObj Scm_PromoteToFlonum(ScmObj obj)
 835 {
 836     if (SCM_INTP(obj)) return Scm_MakeFlonum(SCM_INT_VALUE(obj));
 837     if (SCM_BIGNUMP(obj))
 838         return Scm_MakeFlonum(Scm_BignumToDouble(SCM_BIGNUM(obj)));
 839     if (SCM_FLONUMP(obj)) return obj;
 840     Scm_Panic("Scm_PromoteToFlonum: can't be here");
 841     return SCM_UNDEFINED;       /* dummy */
 842 }
 843 
 844 ScmObj Scm_PromoteToComplex(ScmObj obj)
 845 {
 846     if (SCM_INTP(obj))
 847         return Scm_MakeComplex((double)SCM_INT_VALUE(obj), 0.0);
 848     if (SCM_BIGNUMP(obj))
 849         return Scm_MakeComplex(Scm_BignumToDouble(SCM_BIGNUM(obj)), 0.0);
 850     if (SCM_FLONUMP(obj))
 851         return Scm_MakeComplex(SCM_FLONUM_VALUE(obj), 0.0);
 852     if (SCM_COMPLEXP(obj)) return obj;
 853     Scm_Panic("Scm_PromoteToComplex: can't be here");
 854     return SCM_UNDEFINED;       /* dummy */
 855 }
 856 
 857 /*===============================================================
 858  * Arithmetics
 859  */
 860 /* The code of addition, subtraction, multiplication and division
 861    are somewhat ugly (they use the harmful goto!).  My intention is
 862    to keep intermediate result in C-native types whenever possible,
 863    so that I can avoid boxing/unboxing those numbers. */
 864 
 865 #define APPLY_GENERIC_ARITH(v, gf, arg0, arg1, args)    \
 866   do {                                                  \
 867     v = Scm_Apply(SCM_OBJ(&gf), SCM_LIST2(arg0, arg1)); \
 868     if (SCM_NULLP(args)) return v;                      \
 869     arg1 = SCM_CAR(args);                               \
 870     args = SCM_CDR(args);                               \
 871     goto retry;                                         \
 872   } while (0)
 873 
 874 /*
 875  * Addition and subtraction
 876  */
 877 
 878 ScmObj Scm_Add(ScmObj arg0, ScmObj arg1, ScmObj args)
 879 {
 880     long result_int;
 881     double result_real, result_imag;
 882 
 883   retry:
 884     result_int = 0;
 885     if (SCM_INTP(arg0)) {
 886         result_int = SCM_INT_VALUE(arg0);
 887         for (;;) {
 888             if (SCM_INTP(arg1)) {
 889                 result_int += SCM_INT_VALUE(arg1);
 890                 if (result_int > SCM_SMALL_INT_MAX 
 891                     || result_int < SCM_SMALL_INT_MIN) {
 892                     arg0 = Scm_MakeBignumFromSI(result_int);
 893                     break;
 894                 }
 895             } else if (SCM_BIGNUMP(arg1)) {
 896                 arg0 = Scm_BignumAdd(SCM_BIGNUM(Scm_MakeBignumFromSI(result_int)),
 897                                      SCM_BIGNUM(arg1));
 898                 break;
 899             } else if (SCM_FLONUMP(arg1)) {
 900                 result_real = (double)result_int;
 901                 goto DO_FLONUM;
 902             } else if (SCM_COMPLEXP(arg1)) {
 903                 result_real = (double)result_int;
 904                 result_imag = 0.0;
 905                 goto DO_COMPLEX;
 906             } else {
 907                 APPLY_GENERIC_ARITH(arg0, generic_add,
 908                                     Scm_MakeInteger(result_int),
 909                                     arg1, args);
 910             }
 911             if (!SCM_PAIRP(args)) return Scm_MakeInteger(result_int);
 912             arg1 = SCM_CAR(args);
 913             args = SCM_CDR(args);
 914         }
 915         if (!SCM_PAIRP(args)) return arg0;
 916         arg1 = SCM_CAR(args);
 917         args = SCM_CDR(args);
 918     }
 919     if (SCM_BIGNUMP(arg0)) {
 920         /* See if we should call generic version */
 921         if (SCM_NUMBERP(arg1)) {
 922             return Scm_BignumAddN(SCM_BIGNUM(arg0), Scm_Cons(arg1, args));
 923         } else {
 924             APPLY_GENERIC_ARITH(arg0, generic_add, arg0, arg1, args);
 925         }
 926     }
 927     if (SCM_FLONUMP(arg0)) {
 928         result_real = SCM_FLONUM_VALUE(arg0);
 929       DO_FLONUM:
 930         for (;;) {
 931             if (SCM_INTP(arg1)) {
 932                 result_real += (double)SCM_INT_VALUE(arg1);
 933             } else if (SCM_BIGNUMP(arg1)) {
 934                 result_real += Scm_BignumToDouble(SCM_BIGNUM(arg1));
 935             } else if (SCM_FLONUMP(arg1)) {
 936                 result_real += SCM_FLONUM_VALUE(arg1);
 937             } else if (SCM_COMPLEXP(arg1)) {
 938                 result_imag = 0.0;
 939                 goto DO_COMPLEX;
 940             } else {
 941                 APPLY_GENERIC_ARITH(arg0, generic_add,
 942                                     Scm_MakeFlonum(result_real),
 943                                     arg1, args);
 944             }
 945             if (!SCM_PAIRP(args)) return Scm_MakeFlonum(result_real);
 946             arg1 = SCM_CAR(args);
 947             args = SCM_CDR(args);
 948         }
 949     }
 950     if (SCM_COMPLEXP(arg0)) {
 951         result_real = SCM_COMPLEX_REAL(arg0);
 952         result_imag = SCM_COMPLEX_IMAG(arg0);
 953       DO_COMPLEX:
 954         for (;;) {
 955             if (SCM_INTP(arg1)) {
 956                 result_real += (double)SCM_INT_VALUE(arg1);
 957             } else if (SCM_BIGNUMP(arg1)) {
 958                 result_real += Scm_BignumToDouble(SCM_BIGNUM(arg1));
 959             } else if (SCM_FLONUMP(arg1)) {
 960                 result_real += SCM_FLONUM_VALUE(arg1);
 961             } else if (SCM_COMPLEXP(arg1)) {
 962                 result_real += SCM_COMPLEX_REAL(arg1);
 963                 result_imag += SCM_COMPLEX_IMAG(arg1);
 964             } else {
 965                 APPLY_GENERIC_ARITH(arg0, generic_add,
 966                                     Scm_MakeComplexNormalized(result_real,
 967                                                               result_imag),
 968                                     arg1, args);
 969             }
 970             if (!SCM_PAIRP(args)) {
 971                 return Scm_MakeComplexNormalized(result_real, result_imag);
 972             }
 973             arg1 = SCM_CAR(args);
 974             args = SCM_CDR(args);
 975         }
 976     }
 977     APPLY_GENERIC_ARITH(arg0, generic_add,
 978                         arg0, arg1, args);
 979     return SCM_UNDEFINED;       /* NOTREACHED */
 980 }
 981 
 982 ScmObj Scm_Subtract(ScmObj arg0, ScmObj arg1, ScmObj args)
 983 {
 984     long result_int;
 985     double result_real, result_imag;
 986 
 987   retry:
 988     result_int = 0;
 989     result_real = 0.0;
 990     result_imag = 0.0;
 991     if (SCM_INTP(arg0)) {
 992         result_int = SCM_INT_VALUE(arg0);
 993         for (;;) {
 994             if (SCM_INTP(arg1)) {
 995                 result_int -= SCM_INT_VALUE(arg1);
 996                 if (result_int < SCM_SMALL_INT_MIN
 997                     || result_int > SCM_SMALL_INT_MAX) {
 998                     ScmObj big = Scm_MakeBignumFromSI(result_int);
 999                     return Scm_BignumSubN(SCM_BIGNUM(big), args);
1000                 }
1001             } else if (SCM_BIGNUMP(arg1)) {
1002                 ScmObj big = Scm_MakeBignumFromSI(result_int);
1003                 return Scm_BignumSubN(SCM_BIGNUM(big), Scm_Cons(arg1, args));
1004             } else if (SCM_FLONUMP(arg1)) {
1005                 result_real = (double)result_int;
1006                 goto DO_FLONUM;
1007             } else if (SCM_COMPLEXP(arg1)) {
1008                 result_real = (double)result_int;
1009                 goto DO_COMPLEX;
1010             } else {
1011                 APPLY_GENERIC_ARITH(arg0, generic_sub,
1012                                     Scm_MakeInteger(result_int),
1013                                     arg1, args);
1014             }
1015             if (SCM_NULLP(args)) return SCM_MAKE_INT(result_int);
1016             arg1 = SCM_CAR(args);
1017             args = SCM_CDR(args);
1018         }
1019     }
1020     if (SCM_BIGNUMP(arg0)) {
1021         if (SCM_NUMBERP(arg1)) {
1022             return Scm_BignumSubN(SCM_BIGNUM(arg0), Scm_Cons(arg1, args));
1023         } else {
1024             APPLY_GENERIC_ARITH(arg0, generic_sub, arg0, arg1, args);
1025         }
1026     }
1027     if (SCM_FLONUMP(arg0)) {
1028         result_real = SCM_FLONUM_VALUE(arg0);
1029       DO_FLONUM:
1030         for (;;) {
1031             if (SCM_INTP(arg1)) {
1032                 result_real -= (double)SCM_INT_VALUE(arg1);
1033             } else if (SCM_BIGNUMP(arg1)) {
1034                 result_real -= Scm_BignumToDouble(SCM_BIGNUM(arg1));
1035             } else if (SCM_FLONUMP(arg1)) {
1036                 result_real -= SCM_FLONUM_VALUE(arg1);
1037             } else if (SCM_COMPLEXP(arg1)) {
1038                 goto DO_COMPLEX;
1039             } else {
1040                 APPLY_GENERIC_ARITH(arg0, generic_sub,
1041                                     Scm_MakeFlonum(result_real),
1042                                     arg1, args);
1043             }
1044             if (SCM_NULLP(args))
1045                 return Scm_MakeFlonum(result_real);
1046             arg1 = SCM_CAR(args);
1047             args = SCM_CDR(args);
1048         }
1049     }
1050     if (SCM_COMPLEXP(arg0)) {
1051         result_real = SCM_COMPLEX_REAL(arg0);
1052         result_imag = SCM_COMPLEX_IMAG(arg0);
1053       DO_COMPLEX:
1054         for (;;) {
1055             if (SCM_INTP(arg1)) {
1056                 result_real -= (double)SCM_INT_VALUE(arg1);
1057             } else if (SCM_BIGNUMP(arg1)) {
1058                 result_real -= Scm_BignumToDouble(SCM_BIGNUM(arg1));
1059             } else if (SCM_FLONUMP(arg1)) {
1060                 result_real -= SCM_FLONUM_VALUE(arg1);
1061             } else if (SCM_COMPLEXP(arg1)) {
1062                 result_real -= SCM_COMPLEX_REAL(arg1);
1063                 result_imag -= SCM_COMPLEX_IMAG(arg1);
1064             } else {
1065                 APPLY_GENERIC_ARITH(arg0, generic_sub,
1066                                     Scm_MakeComplexNormalized(result_real,
1067                                                               result_imag),
1068                                     arg1, args);
1069             }
1070             if (SCM_NULLP(args))
1071                 return Scm_MakeComplexNormalized(result_real, result_imag);
1072             arg1 = SCM_CAR(args);
1073             args = SCM_CDR(args);
1074         }
1075     }
1076     APPLY_GENERIC_ARITH(arg0, generic_sub, arg0, arg1, args);
1077     return SCM_UNDEFINED;       /* NOTREACHED */
1078 }
1079 
1080 /*
1081  * Multiplication
1082  */
1083 
1084 ScmObj Scm_Multiply(ScmObj arg0, ScmObj arg1, ScmObj args)
1085 {
1086     long result_int;
1087     double result_real, result_imag;
1088 
1089   retry:
1090     if (SCM_INTP(arg0)) {
1091         result_int = SCM_INT_VALUE(arg0);
1092         for (;;) {
1093             if (SCM_INTP(arg1)) {
1094                 long vv = SCM_INT_VALUE(arg1);
1095                 long k = result_int * vv;
1096                 /* TODO: need a better way to check overflow */
1097                 if ((vv != 0 && k/vv != result_int)
1098                     || k < SCM_SMALL_INT_MIN
1099                     || k > SCM_SMALL_INT_MAX) {
1100                     ScmObj big = Scm_MakeBignumFromSI(result_int);
1101                     arg0 = Scm_BignumMulSI(SCM_BIGNUM(big), vv);
1102                     break;
1103                 }
1104                 result_int = k;
1105             } else if (SCM_BIGNUMP(arg1)) {
1106                 arg0 = Scm_BignumMulSI(SCM_BIGNUM(arg1), result_int);
1107                 break;
1108             } else if (SCM_FLONUMP(arg1)) {
1109                 result_real = (double)result_int;
1110                 goto DO_FLONUM;
1111             } else if (SCM_COMPLEXP(arg1)) {
1112                 result_real = (double)result_int;
1113                 result_imag = 0.0;
1114                 goto DO_COMPLEX;
1115             } else {
1116                 APPLY_GENERIC_ARITH(arg0, generic_mul,
1117                                     Scm_MakeInteger(result_int), arg1, args);
1118             }
1119             if (!SCM_PAIRP(args)) return Scm_MakeInteger(result_int);
1120             arg1 = SCM_CAR(args);
1121             args = SCM_CDR(args);
1122         }
1123         if (!SCM_PAIRP(args)) return arg0;
1124         arg1 = SCM_CAR(args);
1125         args = SCM_CDR(args);
1126         goto retry;
1127     }
1128     if (SCM_BIGNUMP(arg0)) {
1129         return Scm_BignumMulN(SCM_BIGNUM(arg0), Scm_Cons(arg1, args));
1130     }
1131     if (SCM_FLONUMP(arg0)) {
1132         result_real = SCM_FLONUM_VALUE(arg0);
1133       DO_FLONUM:
1134         for (;;) {
1135             if (SCM_INTP(arg1)) {
1136                 result_real *= (double)SCM_INT_VALUE(arg1);
1137             } else if (SCM_BIGNUMP(arg1)) {
1138                 result_real *= Scm_BignumToDouble(SCM_BIGNUM(arg1));
1139             } else if (SCM_FLONUMP(arg1)) {
1140                 result_real *= SCM_FLONUM_VALUE(arg1);
1141             } else if (SCM_COMPLEXP(arg1)) {
1142                 result_imag = 0.0;
1143                 goto DO_COMPLEX;
1144             } else {
1145                 APPLY_GENERIC_ARITH(arg0, generic_mul,
1146                                     Scm_MakeFlonum(result_real), arg1, args);
1147             }
1148             if (!SCM_PAIRP(args)) return Scm_MakeFlonum(result_real);
1149             arg1 = SCM_CAR(args);
1150             args = SCM_CDR(args);
1151         }
1152     }
1153     if (SCM_COMPLEXP(arg0)) {
1154         result_real = SCM_COMPLEX_REAL(arg0);
1155         result_imag = SCM_COMPLEX_IMAG(arg0);
1156       DO_COMPLEX:
1157         for (;;) {
1158             if (SCM_INTP(arg1)) {
1159                 result_real *= (double)SCM_INT_VALUE(arg1);
1160                 result_imag *= (double)SCM_INT_VALUE(arg1);
1161             } else if (SCM_BIGNUMP(arg1)) {
1162                 double dd = Scm_BignumToDouble(SCM_BIGNUM(arg1));
1163                 result_real *= dd;
1164                 result_imag *= dd;
1165             } else if (SCM_FLONUMP(arg1)) {
1166                 result_real *= SCM_FLONUM_VALUE(arg1);
1167                 result_imag *= SCM_FLONUM_VALUE(arg1);
1168             } else if (SCM_COMPLEXP(arg1)) {
1169                 double r = SCM_COMPLEX_REAL(arg1);
1170                 double i = SCM_COMPLEX_IMAG(arg1);
1171                 double t = result_real * r - result_imag * i;
1172                 result_imag   = result_real * i + result_imag * r;
1173                 result_real = t;
1174             } else {
1175                 APPLY_GENERIC_ARITH(arg0, generic_mul,
1176                                     Scm_MakeComplexNormalized(result_real,
1177                                                               result_imag),
1178                                     arg1, args);
1179             }
1180             if (!SCM_PAIRP(args)) {
1181                 return Scm_MakeComplexNormalized(result_real, result_imag);
1182             }
1183             arg1 = SCM_CAR(args);
1184             args = SCM_CDR(args);
1185         }
1186     }
1187     APPLY_GENERIC_ARITH(arg0, generic_mul,
1188                         arg0, arg1, args);
1189     return SCM_UNDEFINED;       /* NOTREACHED */
1190 }
1191 
1192 /*
1193  * Division
1194  */
1195 
1196 ScmObj Scm_Divide(ScmObj arg0, ScmObj arg1, ScmObj args)
1197 {
1198     double result_real = 0.0, result_imag = 0.0;
1199     double div_real = 0.0, div_imag = 0.0;
1200     int exact = TRUE;
1201 
1202   retry:
1203     result_real = result_imag = div_real = div_imag = 0.0;
1204     if (SCM_INTP(arg0)) {
1205         result_real = (double)SCM_INT_VALUE(arg0);
1206         goto DO_FLONUM;
1207     }
1208     if (SCM_BIGNUMP(arg0)) {
1209         /* Try integer division first, and if remainder != 0, shift to
1210            inexact number */
1211         if (SCM_INTP(arg1)) {
1212             long rem;
1213             ScmObj div = Scm_BignumDivSI(SCM_BIGNUM(arg0),
1214                                          SCM_INT_VALUE(arg1),
1215                                          &rem);
1216             if (rem != 0) {
1217                 result_real = Scm_BignumToDouble(SCM_BIGNUM(arg0));
1218                 exact = FALSE;
1219                 goto DO_FLONUM;
1220             }
1221             if (SCM_NULLP(args)) return div;
1222             return Scm_Divide(div, SCM_CAR(args), SCM_CDR(args));
1223         }
1224         if (SCM_BIGNUMP(arg1)) {
1225             ScmObj divrem = Scm_BignumDivRem(SCM_BIGNUM(arg0), SCM_BIGNUM(arg1));
1226             if (SCM_CDR(divrem) != SCM_MAKE_INT(0)) {
1227                 result_real = Scm_BignumToDouble(SCM_BIGNUM(arg0));
1228                 exact = FALSE;
1229                 goto DO_FLONUM;
1230             }
1231             if (SCM_NULLP(args)) return SCM_CAR(divrem);
1232             return Scm_Divide(SCM_CAR(divrem), SCM_CAR(args), SCM_CDR(args));
1233         }
1234         if (SCM_FLONUMP(arg1)) {
1235             exact = FALSE;
1236             result_real = Scm_BignumToDouble(SCM_BIGNUM(arg0));
1237             goto DO_FLONUM;
1238         }
1239         if (SCM_COMPLEXP(arg1)) {
1240             exact = FALSE;
1241             result_real = Scm_BignumToDouble(SCM_BIGNUM(arg0));
1242             goto DO_COMPLEX;
1243         }
1244         APPLY_GENERIC_ARITH(arg0, generic_div, arg0, arg1, args);
1245     }
1246     if (SCM_FLONUMP(arg0)) {
1247         result_real = SCM_FLONUM_VALUE(arg0);
1248         exact = FALSE;
1249       DO_FLONUM:
1250         for (;;) {
1251             if (SCM_INTP(arg1)) {
1252                 div_real = (double)SCM_INT_VALUE(arg1);
1253             } else if (SCM_BIGNUMP(arg1)) {
1254                 div_real = Scm_BignumToDouble(SCM_BIGNUM(arg1));
1255             } else if (SCM_FLONUMP(arg1)) {
1256                 div_real = SCM_FLONUM_VALUE(arg1);
1257                 exact = FALSE;
1258             } else if (SCM_COMPLEXP(arg1)) {
1259                 goto DO_COMPLEX;
1260             } else {
1261                 APPLY_GENERIC_ARITH(arg0, generic_div,
1262                                     Scm_MakeFlonumToNumber(result_real, exact),
1263                                     arg1, args);
1264             }
1265             result_real /= div_real;
1266             if (SCM_NULLP(args))
1267                 return Scm_MakeFlonumToNumber(result_real, exact);
1268             arg1 = SCM_CAR(args);
1269             args = SCM_CDR(args);
1270         }
1271     }
1272     if (SCM_COMPLEXP(arg0)) {
1273         double d, r, i;
1274         result_real = SCM_COMPLEX_REAL(arg0);
1275         result_imag = SCM_COMPLEX_IMAG(arg0);
1276       DO_COMPLEX:
1277         for (;;) {
1278             div_imag = 0.0;
1279             if (SCM_INTP(arg1)) {
1280                 div_real = (double)SCM_INT_VALUE(arg1);
1281             } else if (SCM_BIGNUMP(arg1)) {
1282                 div_real = Scm_BignumToDouble(SCM_BIGNUM(arg1));
1283             } else if (SCM_FLONUMP(arg1)) {
1284                 div_real = SCM_FLONUM_VALUE(arg1);
1285             } else if (SCM_COMPLEXP(arg1)) {
1286                 div_real = SCM_COMPLEX_REAL(arg1);
1287                 div_imag = SCM_COMPLEX_IMAG(arg1);
1288             } else {
1289                 APPLY_GENERIC_ARITH(arg0, generic_div,
1290                                     Scm_MakeComplexNormalized(result_real,
1291                                                               result_imag),
1292                                     arg1, args);
1293             }
1294             d = div_real*div_real + div_imag*div_imag;
1295             r = (result_real*div_real + result_imag*div_imag)/d;
1296             i = (result_imag*div_real - result_real*div_imag)/d;
1297             result_real = r;
1298             result_imag = i;
1299             if (SCM_NULLP(args))
1300                 return Scm_MakeComplexNormalized(result_real, result_imag);
1301             arg1 = SCM_CAR(args);
1302             args = SCM_CDR(args);
1303         }
1304     }
1305     APPLY_GENERIC_ARITH(arg0, generic_div, arg0, arg1, args);
1306     return SCM_UNDEFINED;       /* NOTREACHED */
1307 }
1308 
1309 /*
1310  * Integer division
1311  *   Returns (quotient x y)
1312  *   If rem != NULL, sets *rem to be (remainder x y) as well.
1313  */
1314 ScmObj Scm_Quotient(ScmObj x, ScmObj y, ScmObj *rem)
1315 {
1316     double rx, ry;
1317     if (SCM_INTP(x)) {
1318         if (SCM_INTP(y)) {
1319             long q, r;
1320             if (SCM_INT_VALUE(y) == 0) goto DIVBYZERO;
1321             q = SCM_INT_VALUE(x)/SCM_INT_VALUE(y);
1322             if (rem) {
1323                 r = SCM_INT_VALUE(x)%SCM_INT_VALUE(y);
1324                 *rem = SCM_MAKE_INT(r);
1325             }
1326             return SCM_MAKE_INT(q);
1327         }
1328         if (SCM_BIGNUMP(y)) {
1329             if (rem) *rem = x;
1330             return SCM_MAKE_INT(0);
1331         }
1332         if (SCM_FLONUMP(y)) {
1333             rx = (double)SCM_INT_VALUE(x);
1334             ry = SCM_FLONUM_VALUE(y);
1335             if (ry != floor(ry)) goto BADARGY;
1336             goto DO_FLONUM;
1337         }
1338         goto BADARGY;
1339     } else if (SCM_BIGNUMP(x)) {
1340         if (SCM_INTP(y)) {
1341             long r;
1342             ScmObj q = Scm_BignumDivSI(SCM_BIGNUM(x), SCM_INT_VALUE(y), &r);
1343             if (rem) *rem = SCM_MAKE_INT(r);
1344             return q;
1345         } else if (SCM_BIGNUMP(y)) {
1346             ScmObj qr = Scm_BignumDivRem(SCM_BIGNUM(x), SCM_BIGNUM(y));
1347             if (rem) *rem = SCM_CDR(qr);
1348             return SCM_CAR(qr);
1349         } else if (SCM_FLONUMP(y)) {
1350             rx = Scm_BignumToDouble(SCM_BIGNUM(x));
1351             ry = SCM_FLONUM_VALUE(y);
1352             if (ry != floor(ry)) goto BADARGY;
1353             goto DO_FLONUM;
1354         }
1355         goto BADARGY;
1356     } else if (SCM_FLONUMP(x)) {
1357         rx = SCM_FLONUM_VALUE(x);
1358         if (rx != floor(rx)) goto BADARG;
1359         if (SCM_INTP(y)) {
1360             ry = (double)SCM_INT_VALUE(y);
1361         } else if (SCM_BIGNUMP(y)) {
1362             ry = Scm_BignumToDouble(SCM_BIGNUM(y));
1363         } else if (SCM_FLONUMP(y)) {
1364             ry = SCM_FLONUM_VALUE(y);
1365             if (ry != floor(ry)) goto BADARGY;
1366         } else {
1367             goto BADARGY;
1368         }
1369       DO_FLONUM:
1370         {
1371             double q;
1372             if (ry == 0.0) goto DIVBYZERO;
1373             q = roundeven(rx/ry);
1374             if (rem) {
1375                 double rr = roundeven(rx - q*ry);
1376                 *rem = Scm_MakeFlonum(rr);
1377             }
1378             return Scm_MakeFlonum(q);
1379         }
1380     } else {
1381         goto BADARG;
1382     }
1383   DIVBYZERO:
1384     Scm_Error("attempt to calculate a quotient by zero");
1385   BADARGY:
1386     x = y;
1387   BADARG:
1388     Scm_Error("integer required, but got %S", x);
1389     return SCM_UNDEFINED;       /* dummy */
1390 }
1391 
1392 /* Modulo and Reminder.
1393    TODO: on gcc, % works like reminder.  I'm not sure the exact behavior
1394    of % is defined in ANSI C.  Need to check it later. */
1395 ScmObj Scm_Modulo(ScmObj x, ScmObj y, int remp)
1396 {
1397     double rx, ry;
1398     if (SCM_INTP(x)) {
1399         if (SCM_INTP(y)) {
1400             long r;
1401             if (SCM_INT_VALUE(y) == 0) goto DIVBYZERO;
1402             r = SCM_INT_VALUE(x)%SCM_INT_VALUE(y);
1403             if (!remp && r) {
1404                 if ((SCM_INT_VALUE(x) > 0 && SCM_INT_VALUE(y) < 0)
1405                     || (SCM_INT_VALUE(x) < 0 && SCM_INT_VALUE(y) > 0)) {
1406                     r += SCM_INT_VALUE(y);
1407                 }
1408             }
1409             return SCM_MAKE_INT(r);
1410         }
1411         if (SCM_BIGNUMP(y)) {
1412             if (remp) {
1413                 return x;
1414             } else {
1415                 if ((SCM_INT_VALUE(x) < 0 && SCM_BIGNUM_SIGN(y) > 0)
1416                     || (SCM_INT_VALUE(x) > 0 && SCM_BIGNUM_SIGN(y) < 0)) {
1417                     return Scm_BignumAddSI(SCM_BIGNUM(y), SCM_INT_VALUE(x));
1418                 } else {
1419                     return x;
1420                 }
1421             }
1422         }
1423         rx = (double)SCM_INT_VALUE(x);
1424         if (SCM_FLONUMP(y)) {
1425             ry = SCM_FLONUM_VALUE(y);
1426             if (ry != floor(ry)) goto BADARGY;
1427             goto DO_FLONUM;
1428         }
1429         goto BADARGY;
1430     } else if (SCM_BIGNUMP(x)) {
1431         if (SCM_INTP(y)) {
1432             long iy = SCM_INT_VALUE(y);
1433             long rem;
1434             Scm_BignumDivSI(SCM_BIGNUM(x), iy, &rem);
1435             if (!remp
1436                 && rem
1437                 && ((SCM_BIGNUM_SIGN(x) < 0 && iy > 0)
1438                     || (SCM_BIGNUM_SIGN(x) > 0 && iy < 0))) {
1439                 return SCM_MAKE_INT(iy + rem);
1440             }
1441             return SCM_MAKE_INT(rem);
1442         }
1443         if (SCM_BIGNUMP(y)) {
1444             ScmObj rem = SCM_CDR(Scm_BignumDivRem(SCM_BIGNUM(x), SCM_BIGNUM(y)));
1445             if (!remp
1446                 && (rem != SCM_MAKE_INT(0))
1447                 && (SCM_BIGNUM_SIGN(x) * SCM_BIGNUM_SIGN(y) < 0)) {
1448                 if (SCM_BIGNUMP(rem)) {
1449                     return Scm_BignumAdd(SCM_BIGNUM(y), SCM_BIGNUM(rem));
1450                 } else {
1451                     return Scm_BignumAddSI(SCM_BIGNUM(y), SCM_INT_VALUE(rem));
1452                 }       
1453             }
1454             return rem;
1455         }
1456         rx = Scm_BignumToDouble(SCM_BIGNUM(x));
1457         if (SCM_FLONUMP(y)) {
1458             ry = SCM_FLONUM_VALUE(y);
1459             if (ry != floor(ry)) goto BADARGY;
1460             goto DO_FLONUM;
1461         }
1462         goto BADARGY;
1463     } else if (SCM_FLONUMP(x)) {
1464         double rem;
1465         rx = SCM_FLONUM_VALUE(x);
1466         if (rx != floor(rx)) goto BADARG;
1467         if (SCM_INTP(y)) {
1468             ry = (double)SCM_INT_VALUE(y);
1469         } else if (SCM_BIGNUMP(y)) {
1470             ry = Scm_BignumToDouble(SCM_BIGNUM(y));
1471         } else if (SCM_FLONUMP(y)) {
1472             ry = SCM_FLONUM_VALUE(y);
1473             if (ry != floor(ry)) goto BADARGY;
1474         } else {
1475             goto BADARGY;
1476         }
1477       DO_FLONUM:
1478         if (ry == 0.0) goto DIVBYZERO;
1479         rem = fmod(rx, ry);
1480         if (!remp && rem != 0.0) {
1481             if ((rx > 0 && ry < 0) || (rx < 0 && ry > 0)) {
1482                 rem += ry;
1483             }
1484         }
1485         return Scm_MakeFlonum(rem);
1486     } else {
1487         goto BADARG;
1488     }
1489   DIVBYZERO:
1490     Scm_Error("attempt to take a modulo or remainder by zero");
1491   BADARGY:
1492     x = y;
1493   BADARG:
1494     Scm_Error("integer required, but got %S", x);
1495     return SCM_UNDEFINED;       /* dummy */
1496 }
1497 
1498 /*
1499  * Expt
1500  */
1501 
1502 /* Integer power of 10.  It is extensively used during string->number
1503    and number->string operations.
1504    IEXPT10_TABLESIZ is ceil(-log10(ldexp(1.0, -1022-52))) + 2 */
1505 /* NB: actually we need more margin here to handle denormalized numbers. */
1506 #define IEXPT10_TABLESIZ  341
1507 static ScmObj iexpt10_n[IEXPT10_TABLESIZ] = { NULL };
1508 static int    iexpt10_initialized = FALSE;
1509 
1510 static void iexpt10_init(void)
1511 {
1512     int i;
1513     iexpt10_n[0] = SCM_MAKE_INT(1);
1514     iexpt10_n[1] = SCM_MAKE_INT(10);
1515     iexpt10_n[2] = SCM_MAKE_INT(100);
1516     iexpt10_n[3] = SCM_MAKE_INT(1000);
1517     iexpt10_n[4] = SCM_MAKE_INT(10000);
1518     iexpt10_n[5] = SCM_MAKE_INT(100000);
1519     iexpt10_n[6] = SCM_MAKE_INT(1000000);
1520     for (i=7; i<IEXPT10_TABLESIZ; i++) {
1521         iexpt10_n[i] = Scm_Multiply2(iexpt10_n[i-1], SCM_MAKE_INT(10));
1522     }
1523     iexpt10_initialized = TRUE;
1524 }
1525 
1526 #define IEXPT10_INIT() \
1527     do { if (!iexpt10_initialized) iexpt10_init(); } while (0)
1528 
1529 /* short cut for exact numbers */
1530 static ScmObj exact_expt(ScmObj x, ScmObj y)
1531 {
1532     int sign = Scm_Sign(y);
1533     long iy;
1534     ScmObj r = SCM_MAKE_INT(1);
1535 
1536     if (sign == 0) return r;
1537     if (x == SCM_MAKE_INT(1)) return r;
1538     if (x == SCM_MAKE_INT(-1)) return Scm_OddP(y)? SCM_MAKE_INT(-1) : r;
1539 
1540     if (!SCM_INTP(y)) {
1541         /* who wants such a heavy calculation? */
1542         Scm_Error("exponent too big: %S", y);
1543     }
1544     iy = SCM_INT_VALUE(y);
1545     /* Shortcut for special cases */
1546     if (x == SCM_MAKE_INT(10) && iy > 0 && iy < IEXPT10_TABLESIZ) {
1547         IEXPT10_INIT();
1548         r = iexpt10_n[iy];
1549     } else if (x == SCM_MAKE_INT(2) && iy > 0) {
1550         r = Scm_Ash(SCM_MAKE_INT(1), iy);
1551     } else {
1552         if (iy < 0) iy = -iy;
1553         for (;;) {
1554             if (iy == 0) break;
1555             if (iy == 1) { r = Scm_Multiply2(r, x); break; }
1556             if (iy & 0x01) r = Scm_Multiply2(r, x);
1557             x = Scm_Multiply2(x, x);
1558             iy >>= 1;
1559         }
1560     }
1561     return (sign < 0)? Scm_Reciprocal(r) : r;
1562 }
1563 
1564 ScmObj Scm_Expt(ScmObj x, ScmObj y)
1565 {
1566     double dx, dy;
1567     if (SCM_EXACTP(x) && SCM_EXACTP(y)) return exact_expt(x, y);
1568     if (!SCM_REALP(x)) Scm_Error("real number required, but got %S", x);
1569     if (!SCM_REALP(y)) Scm_Error("real number required, but got %S", y);
1570     dx = Scm_GetDouble(x);
1571     dy = Scm_GetDouble(y);
1572     if (dy == 0.0) {
1573         return Scm_MakeFlonum(1.0);
1574     } else if (dx < 0 && !Scm_IntegerP(y)) {
1575         /* x^y == exp(y * log(x)) = exp(y*log(|x|))*exp(y*arg(x)*i)
1576            if x is a negative real number, arg(x) == pi
1577         */
1578         double mag = exp(dy * log(-dx));
1579         double theta = dy * M_PI;
1580         return Scm_MakeComplexNormalized(mag * cos(theta), mag * sin(theta));
1581     } else {
1582         return Scm_MakeFlonum(pow(dx, dy));
1583     }
1584 }
1585 
1586 /*===============================================================
1587  * Comparison
1588  */
1589 
1590 int Scm_NumEq(ScmObj arg0, ScmObj arg1)
1591 {
1592     if (SCM_COMPLEXP(arg0)) {
1593         if (SCM_COMPLEXP(arg1)) {
1594             return ((SCM_COMPLEX_REAL(arg0) == SCM_COMPLEX_REAL(arg1))
1595                     && (SCM_COMPLEX_IMAG(arg0) == SCM_COMPLEX_IMAG(arg1)));
1596         }
1597         return FALSE;
1598     } else {
1599         if (SCM_COMPLEXP(arg1)) return FALSE;
1600         return (Scm_NumCmp(arg0, arg1) == 0);
1601     }
1602 }
1603 
1604 /* 2-arg comparison */
1605 int Scm_NumCmp(ScmObj arg0, ScmObj arg1)
1606 {
1607     ScmObj badnum;
1608     
1609     if (SCM_INTP(arg0)) {
1610         if (SCM_INTP(arg1)) {
1611             long r = SCM_INT_VALUE(arg0) - SCM_INT_VALUE(arg1);
1612             if (r < 0) return -1;
1613             if (r > 0) return 1;
1614             return 0;
1615         }
1616         if (SCM_FLONUMP(arg1)) {
1617             double r = SCM_INT_VALUE(arg0) - SCM_FLONUM_VALUE(arg1);
1618             if (r < 0) return -1;
1619             if (r > 0) return 1;
1620             return 0;
1621         }
1622         if (SCM_BIGNUMP(arg1))
1623             return Scm_BignumCmp(SCM_BIGNUM(Scm_MakeBignumFromSI(SCM_INT_VALUE(arg0))),
1624                                  SCM_BIGNUM(arg1));
1625         badnum = arg1;
1626     }
1627     else if (SCM_FLONUMP(arg0)) {
1628         if (SCM_INTP(arg1)) {
1629             double r = SCM_FLONUM_VALUE(arg0) - SCM_INT_VALUE(arg1);
1630             if (r < 0) return -1;
1631             if (r > 0) return 1;
1632             return 0;
1633         }
1634         if (SCM_FLONUMP(arg1)) {
1635             double r = SCM_FLONUM_VALUE(arg0) - SCM_FLONUM_VALUE(arg1);
1636             if (r < 0) return -1;
1637             if (r > 0) return 1;
1638             return 0;
1639         }
1640         if (SCM_BIGNUMP(arg1))
1641             return Scm_BignumCmp(SCM_BIGNUM(Scm_MakeBignumFromDouble(SCM_FLONUM_VALUE(arg0))),
1642                                  SCM_BIGNUM(arg1));
1643         badnum = arg1;
1644     }
1645     else if (SCM_BIGNUMP(arg0)) {
1646         if (SCM_INTP(arg1))
1647             return Scm_BignumCmp(SCM_BIGNUM(arg0),
1648                                  SCM_BIGNUM(Scm_MakeBignumFromSI(SCM_INT_VALUE(arg1))));
1649         if (SCM_FLONUMP(arg1))
1650             return Scm_BignumCmp(SCM_BIGNUM(arg0),
1651                                  SCM_BIGNUM(Scm_MakeBignumFromDouble(SCM_FLONUM_VALUE(arg1))));
1652         if (SCM_BIGNUMP(arg1))
1653             return Scm_BignumCmp(SCM_BIGNUM(arg0), SCM_BIGNUM(arg1));
1654         badnum = arg1;
1655     }
1656     else badnum = arg0;
1657     Scm_Error("real number required: %S", badnum);
1658     return 0;                    /* dummy */
1659 }
1660 
1661 void Scm_MinMax(ScmObj arg0, ScmObj args, ScmObj *min, ScmObj *max)
1662 {
1663     int inexact = !SCM_EXACTP(arg0);
1664     ScmObj mi = arg0;
1665     ScmObj ma = arg0;
1666     
1667     for (;;) {
1668         if (!SCM_REALP(arg0))
1669             Scm_Error("real number required, but got %S", arg0);
1670         if (SCM_NULLP(args)) {
1671             if (min) {
1672                 if (inexact && SCM_EXACTP(mi)) {
1673                     *min = Scm_ExactToInexact(mi);
1674                 } else {
1675                     *min = mi;
1676                 }
1677             }
1678             if (max) {
1679                 if (inexact && SCM_EXACTP(ma)) {
1680                     *max = Scm_ExactToInexact(ma);
1681                 } else {
1682                     *max = ma;
1683                 }
1684             }
1685             return;
1686         }
1687         if (!SCM_EXACTP(SCM_CAR(args))) inexact = TRUE;
1688         if (min && Scm_NumCmp(mi, SCM_CAR(args)) > 0) {
1689             mi = SCM_CAR(args);
1690         } 
1691         if (max && Scm_NumCmp(ma, SCM_CAR(args)) < 0) {
1692             ma = SCM_CAR(args);
1693         }
1694         args = SCM_CDR(args);
1695     }
1696 }
1697 
1698 /*===============================================================
1699  * ROUNDING
1700  */
1701 
1702 /* NB: rint() is not in POSIX, so the alternative is provided here.
1703    We don't use round(), for it behaves differently when the
1704    argument is exactly the halfway of two whole numbers. */
1705 #ifdef HAVE_RINT
1706 #define roundeven rint
1707 #else  /* !HAVE_RINT */
1708 static double roundeven(double v)
1709 {
1710     double r;
1711     double frac = modf(v, &r);
1712     if (v > 0.0) {
1713         if (frac > 0.5) r += 1.0;
1714         else if (frac == 0.5) {
1715             if (fmod(r,2.0) != 0.0) r += 1.0;
1716         }
1717     } else {
1718         if (frac < -0.5) r -= 1.0;
1719         else if (frac == -0.5) {
1720             if (fmod(r,2.0) != 0.0) r -= 1.0;
1721         }
1722     }
1723     return r;
1724 }
1725 #endif /* !HAVE_RINT */
1726 
1727 ScmObj Scm_Round(ScmObj num, int mode)
1728 {
1729     double r = 0.0, v;
1730     
1731     if (SCM_EXACTP(num)) return num;
1732     if (!SCM_FLONUMP(num))
1733         Scm_Error("real number required, but got %S", num);
1734     v = SCM_FLONUM_VALUE(num);
1735     switch (mode) {
1736     case SCM_ROUND_FLOOR: r = floor(v); break;
1737     case SCM_ROUND_CEIL:  r = ceil(v); break;
1738     /* trunc is neither in ANSI nor in POSIX. */
1739 #ifdef HAVE_TRUNC
1740     case SCM_ROUND_TRUNC: r = trunc(v); break;
1741 #else
1742     case SCM_ROUND_TRUNC: r = (v < 0.0)? ceil(v) : floor(v); break;
1743 #endif
1744     case SCM_ROUND_ROUND: r = roundeven(v); break;
1745     default: Scm_Panic("something screwed up");
1746     }
1747     return Scm_MakeFlonum(r);
1748 }
1749 
1750 /*===============================================================
1751  * Logical (bitwise) operations
1752  */
1753 
1754 ScmObj Scm_Ash(ScmObj x, int cnt)
1755 {
1756     if (SCM_INTP(x)) {
1757         long ix = SCM_INT_VALUE(x);
1758         if (cnt <= -(SIZEOF_LONG * 8)) {
1759             ix = (ix < 0)? -1 : 0;
1760             return Scm_MakeInteger(ix);
1761         } else if (cnt < 0) {
1762             if (ix < 0) {
1763                 ix = ~((~ix) >> (-cnt));
1764             } else {
1765                 ix >>= -cnt;
1766             }
1767             return Scm_MakeInteger(ix);
1768         } else if (cnt < (SIZEOF_LONG*8-3)) {
1769             if (ix < 0) {
1770                 if (-ix < (SCM_SMALL_INT_MAX >> cnt)) {
1771                     ix <<= cnt;
1772                     return Scm_MakeInteger(ix);
1773                 } 
1774             } else {
1775                 if (ix < (SCM_SMALL_INT_MAX >> cnt)) {
1776                     ix <<= cnt;
1777                     return Scm_MakeInteger(ix);
1778                 } 
1779             }
1780         }
1781         /* Here, we know the result must be a bignum. */
1782         {
1783             ScmObj big = Scm_MakeBignumFromSI(ix);
1784             return Scm_BignumAsh(SCM_BIGNUM(big), cnt);
1785         }
1786     } else if (SCM_BIGNUMP(x)) {
1787         return Scm_BignumAsh(SCM_BIGNUM(x), cnt);
1788     }
1789     Scm_Error("exact integer required, but got %S", x);
1790     return SCM_UNDEFINED;
1791 }
1792 
1793 ScmObj Scm_LogNot(ScmObj x)
1794 {
1795     if (!SCM_EXACTP(x)) Scm_Error("exact integer required, but got %S", x);
1796     if (SCM_INTP(x)) {
1797         /* this won't cause an overflow */
1798         return SCM_MAKE_INT(~SCM_INT_VALUE(x));
1799     } else {
1800         return Scm_Negate(Scm_BignumAddSI(SCM_BIGNUM(x), 1));
1801     }
1802 }
1803 
1804 ScmObj Scm_LogAnd(ScmObj x, ScmObj y)
1805 {
1806     if (!SCM_EXACTP(x)) Scm_Error("exact integer required, but got %S", x);
1807     if (!SCM_EXACTP(y)) Scm_Error("exact integer required, but got %S", y);
1808     if (SCM_INTP(x)) {
1809         if (SCM_INTP(y)) {
1810             return SCM_MAKE_INT(SCM_INT_VALUE(x) & SCM_INT_VALUE(y));
1811         } else if (SCM_INT_VALUE(x) >= 0 && SCM_BIGNUM_SIGN(y) >= 0) {
1812             return Scm_MakeInteger(SCM_INT_VALUE(x)&SCM_BIGNUM(y)->values[0]);
1813         }
1814         x = Scm_MakeBignumFromSI(SCM_INT_VALUE(x));
1815     } else if (SCM_INTP(y)) {
1816         if (SCM_INT_VALUE(y) >= 0 && SCM_BIGNUM_SIGN(x) >= 0) {
1817             return Scm_MakeInteger(SCM_INT_VALUE(y)&SCM_BIGNUM(x)->values[0]);
1818         }
1819         y = Scm_MakeBignumFromSI(SCM_INT_VALUE(y));        
1820     }
1821     return Scm_BignumLogAnd(SCM_BIGNUM(x), SCM_BIGNUM(y));
1822 }
1823 
1824 ScmObj Scm_LogIor(ScmObj x, ScmObj y)
1825 {
1826     if (!SCM_EXACTP(x)) Scm_Error("exact integer required, but got %S", x);
1827     if (!SCM_EXACTP(y)) Scm_Error("exact integer required, but got %S", y);
1828     if (SCM_INTP(x)) {
1829         if (SCM_INTP(y))
1830             return SCM_MAKE_INT(SCM_INT_VALUE(x) | SCM_INT_VALUE(y));
1831         else
1832             x = Scm_MakeBignumFromSI(SCM_INT_VALUE(x));
1833     } else {
1834         if (SCM_INTP(y)) y = Scm_MakeBignumFromSI(SCM_INT_VALUE(y));
1835     }
1836     return Scm_BignumLogIor(SCM_BIGNUM(x), SCM_BIGNUM(y));
1837 }
1838 
1839 
1840 ScmObj Scm_LogXor(ScmObj x, ScmObj y)
1841 {
1842     if (!SCM_EXACTP(x)) Scm_Error("exact integer required, but got %S", x);
1843     if (!SCM_EXACTP(y)) Scm_Error("exact integer required, but got %S", y);
1844     if (SCM_INTP(x)) {
1845         if (SCM_INTP(y))
1846             return SCM_MAKE_INT(SCM_INT_VALUE(x) ^ SCM_INT_VALUE(y));
1847         else
1848             x = Scm_MakeBignumFromSI(SCM_INT_VALUE(x));
1849     } else {
1850         if (SCM_INTP(y)) y = Scm_MakeBignumFromSI(SCM_INT_VALUE(y));
1851     }
1852     return Scm_BignumLogXor(SCM_BIGNUM(x), SCM_BIGNUM(y));
1853 }
1854 
1855 /*===============================================================
1856  * Number I/O
1857  */
1858 
1859 /* contants frequently used in number I/O */
1860 static double dexpt2_minus_52  = 0.0;  /* 2.0^-52 */
1861 static double dexpt2_minus_53  = 0.0;  /* 2.0^-53 */
1862 
1863 /* max N where 10.0^N can be representable exactly in double.
1864    it is max N where N * log2(5) < 53. */
1865 #define MAX_EXACT_10_EXP  23
1866 
1867 /* fast 10^n for limited cases */
1868 static inline ScmObj iexpt10(int e)
1869 {
1870     SCM_ASSERT(e < IEXPT10_TABLESIZ);
1871     return iexpt10_n[e];
1872 }
1873 
1874 /* integer power of R by N, N is rather small.
1875    Assuming everything is in range. */
1876 static inline u_long ipow(int r, int n)
1877 {
1878     u_long k;
1879     for (k=1; n>0; n--) k *= r;
1880     return k;
1881 }
1882 
1883 /* X * 10.0^N by double.
1884    10.0^N can be represented _exactly_ in double-precision floating point
1885    number in the range 0 <= N <= 23.
1886    If N is out of this range, a rounding error occurs, which will be
1887    corrected in the algorithmR routine below. */
1888 static double raise_pow10(double x, int n)
1889 {
1890     static double dpow10[] = { 1.0, 1.0e1, 1.0e2, 1.0e3, 1.0e4,
1891                                1.0e5, 1.0e6, 1.0e7, 1.0e8, 1.0e9,
1892                                1.0e10, 1.0e11, 1.0e12, 1.0e13, 1.0e14,
1893                                1.0e15, 1.0e16, 1.0e17, 1.0e18, 1.0e19,
1894                                1.0e20, 1.0e21, 1.0e22, 1.0e23 };
1895     if (n >= 0) {
1896         while (n > 23) {
1897             x *= 1.0e24;
1898             n -= 24;
1899         }
1900         return x*dpow10[n];
1901     } else {
1902         while (n < -23) {
1903             x /= 1.0e24;
1904             n += 24;
1905         }
1906         return x/dpow10[-n];
1907     }
1908 }
1909 
1910 /*
1911  * Number Printer
1912  *
1913  * This version implements Burger&Dybvig algorithm (Robert G. Burger
1914  * and and R. Kent Dybvig, "Priting Floating-Point Numbers Quickly and 
1915  * Accurately", PLDI '96, pp.108--116, 1996).
1916  */
1917 
1918 /* compare x+d and y */
1919 static inline int numcmp3(ScmObj x, ScmObj d, ScmObj y)
1920 {
1921     ScmObj bx = SCM_BIGNUMP(x)? x : Scm_MakeBignumFromSI(SCM_INT_VALUE(x));
1922     ScmObj bd = SCM_BIGNUMP(d)? d : Scm_MakeBignumFromSI(SCM_INT_VALUE(d));
1923     ScmObj by = SCM_BIGNUMP(y)? y : Scm_MakeBignumFromSI(SCM_INT_VALUE(y));
1924     return Scm_BignumCmp3U(SCM_BIGNUM(bx), SCM_BIGNUM(bd), SCM_BIGNUM(by));
1925 }
1926 
1927 static void double_print(char *buf, int buflen, double val, int plus_sign)
1928 {
1929     /* Handle a few special cases first.
1930        The notation of infinity is provisional; see how srfi-70 becomes. */
1931     if (val == 0.0) {
1932         if (plus_sign) strcpy(buf, "+0.0");
1933         else strcpy(buf, "0.0");
1934         return;
1935     } else if (SCM_IS_INF(val)) {
1936         if (val < 0.0) strcpy(buf, "#i-1/0");
1937         else if (plus_sign) strcpy(buf, "#i+1/0");
1938         else strcpy(buf, "#i1/0");
1939         return;
1940     } else if (SCM_IS_NAN(val)) {
1941         strcpy(buf, "#<nan>");
1942         return;
1943     }
1944     
1945     if (val < 0.0) *buf++ = '-', buflen--;
1946     else if (plus_sign) *buf++ = '+', buflen--;
1947     {
1948         /* variable names follows Burger&Dybvig paper. mp, mm for m+, m-.
1949            note that m+ == m- for most cases, and m+ == 2*m- for the rest.
1950            so we calculate m+ from m- for each iteration, using the flag
1951            mp2 as   m+ = mp? m- : 2*m-. */
1952         ScmObj f, r, s, mp, mm, q;
1953         int exp, sign, est, tc1, tc2, tc3, digs, point, round;
1954         int mp2 = FALSE, fixup = FALSE;
1955 
1956         IEXPT10_INIT();
1957         if (val < 0) val = -val;
1958         
1959         /* initialize r, s, m+ and m- */
1960         f = Scm_DecodeFlonum(val, &exp, &sign);
1961         round = !Scm_OddP(f);
1962         if (exp >= 0) {
1963             ScmObj be = Scm_Ash(SCM_MAKE_INT(1), exp);
1964             if (Scm_NumCmp(f, SCM_2_52) != 0) {
1965                 r = Scm_Ash(f, exp+1);
1966                 s = SCM_MAKE_INT(2);
1967                 mp2= FALSE;
1968                 mm = be;
1969             } else {
1970                 r = Scm_Ash(f, exp+2);
1971                 s = SCM_MAKE_INT(4);
1972                 mp2 = TRUE;
1973                 mm = be;
1974             }
1975         } else {
1976             if (exp == -1023 || Scm_NumCmp(f, SCM_2_52) != 0) {
1977                 r = Scm_Ash(f, 1);
1978                 s = Scm_Ash(SCM_MAKE_INT(1), -exp+1);
1979                 mp2 = FALSE;
1980                 mm = SCM_MAKE_INT(1);
1981             } else {
1982                 r = Scm_Ash(f, 2);
1983                 s = Scm_Ash(SCM_MAKE_INT(1), -exp+2);
1984                 mp2 = TRUE;
1985                 mm = SCM_MAKE_INT(1);
1986             }
1987         }
1988 
1989         /* estimate scale */
1990         est = (int)ceil(log10(val) - 0.1);
1991         if (est >= 0) {
1992             s = Scm_Multiply2(s, iexpt10(est));
1993         } else {
1994             ScmObj scale = iexpt10(-est);
1995             r =  Scm_Multiply2(r, scale);
1996             mm = Scm_Multiply2(mm, scale);
1997         }
1998 
1999         /* fixup.  avoid calculating m+ for obvious case. */
2000         if (Scm_NumCmp(r, s) >= 0) {
2001             fixup = TRUE;
2002         } else {
2003             mp = (mp2? Scm_Ash(mm, 1) : mm);
2004             if (round) {
2005                 fixup = (numcmp3(r, mp, s) >= 0);
2006             } else {
2007                 fixup = (numcmp3(r, mp, s) > 0);
2008             }
2009         }
2010         if (fixup) {
2011             s = Scm_Multiply2(s, SCM_MAKE_INT(10));
2012             est++;
2013         }
2014         
2015         /* Scm_Printf(SCM_CURERR, "est=%d, r=%S, s=%S, mp=%S, mm=%S\n",
2016            est, r, s, mp, mm); */
2017 
2018         /* determine position of decimal point.  we avoid exponential
2019            notation if exponent is small, i.e. 0.9 and 30.0 instead of
2020            9.0e-1 and 3.0e1.   The magic number 10 is arbitrary. */
2021         if (est < 10 && est > -3) {
2022             point = est; est = 1;
2023         } else {
2024             point = 1;
2025         }
2026 
2027         /* generate */
2028         if (point <= 0) {
2029             *buf++ = '0'; buflen--;
2030             *buf++ = '.', buflen--;
2031             for (digs=point;digs<0 && buflen>5;digs++) {
2032                 *buf++ = '0'; buflen--;
2033             }
2034         }
2035         for (digs=1;buflen>5;digs++) {
2036             ScmObj r10 = Scm_Multiply2(r, SCM_MAKE_INT(10));
2037             q = Scm_Quotient(r10, s, &r);
2038             mm = Scm_Multiply2(mm, SCM_MAKE_INT(10));
2039             mp = (mp2? Scm_Ash(mm, 1) : mm);
2040             
2041             /* Scm_Printf(SCM_CURERR, "q=%S, r=%S, mp=%S, mm=%S\n",
2042                q, r, mp, mm);*/
2043 
2044             SCM_ASSERT(SCM_INTP(q));
2045             if (round) {
2046                 tc1 = (Scm_NumCmp(r, mm) <= 0);
2047                 tc2 = (numcmp3(r, mp, s) >= 0);
2048             } else {
2049                 tc1 = (Scm_NumCmp(r, mm) < 0);
2050                 tc2 = (numcmp3(r, mp, s) > 0);
2051             }
2052             if (!tc1) {
2053                 if (!tc2) {
2054                     *buf++ = SCM_INT_VALUE(q) + '0';
2055                     if (digs == point) *buf++ = '.', buflen--;
2056                     continue;
2057                 } else {
2058                     *buf++ = SCM_INT_VALUE(q) + '1';
2059                     break;
2060                 }
2061             } else {
2062                 if (!tc2) {
2063                     *buf++ = SCM_INT_VALUE(q) + '0';
2064                     break;
2065                 } else {
2066                     tc3 = numcmp3(r, r, s); /* r*2 <=> s */
2067                     if ((round && tc3 <= 0) || (!round && tc3 < 0)) {
2068                         *buf++ = SCM_INT_VALUE(q) + '0';
2069                         break;
2070                     } else {
2071                         *buf++ = SCM_INT_VALUE(q) + '1';
2072                         break;
2073                     }
2074                 }
2075             }
2076         }
2077 
2078         if (digs <= point) {
2079             for (;digs<point&&buflen>5;digs++) {
2080                 *buf++ = '0', buflen--;
2081             }
2082             *buf++ = '.';
2083             *buf++ = '0';
2084         }
2085 
2086         /* prints exponent.  we shifted decimal point, so -1. */
2087         est--;
2088         if (est != 0) {
2089             *buf++ = 'e';
2090             sprintf(buf, "%d", (int)est);
2091         } else {
2092             *buf++ = 0;
2093         }
2094     }
2095 }
2096 
2097 static void number_print(ScmObj obj, ScmPort *port, ScmWriteContext *ctx)
2098 {
2099     ScmObj s = Scm_NumberToString(obj, 10, FALSE);
2100     SCM_PUTS(SCM_STRING(s), port);
2101 }
2102 
2103 #define FLT_BUF 50
2104 
2105 ScmObj Scm_NumberToString(ScmObj obj, int radix, int use_upper)
2106 {
2107     ScmObj r = SCM_NIL;
2108     char buf[FLT_BUF];
2109     
2110     if (SCM_INTP(obj)) {
2111         char buf[50], *pbuf = buf;
2112         long value = SCM_INT_VALUE(obj);
2113         if (value < 0) {
2114             *pbuf++ = '-';
2115             value = -value;     /* this won't overflow */
2116         }
2117         if (radix == 10) {
2118             snprintf(pbuf, 49, "%ld", value);
2119         } else if (radix == 16) {
2120             snprintf(pbuf, 49, (use_upper? "%lX" : "%lx"), value);
2121         } else if (radix == 8) {
2122             snprintf(pbuf, 49, "%lo", value);
2123         } else {
2124             /* sloppy way ... */
2125             r = Scm_BignumToString(SCM_BIGNUM(Scm_MakeBignumFromSI(SCM_INT_VALUE(obj))),
2126                                    radix, use_upper);
2127         }
2128         if (r == SCM_NIL) r = SCM_MAKE_STR_COPYING(buf);
2129     } else if (SCM_BIGNUMP(obj)) {
2130         r = Scm_BignumToString(SCM_BIGNUM(obj), radix, use_upper);
2131     } else if (SCM_FLONUMP(obj)) {
2132         double_print(buf, FLT_BUF, SCM_FLONUM_VALUE(obj), FALSE);
2133         r = SCM_MAKE_STR_COPYING(buf);
2134     } else if (SCM_COMPLEXP(obj)) {
2135         ScmObj p = Scm_MakeOutputStringPort(TRUE);
2136         double_print(buf, FLT_BUF, SCM_COMPLEX_REAL(obj), FALSE);
2137         SCM_PUTZ(buf, -1, SCM_PORT(p));
2138         double_print(buf, FLT_BUF, SCM_COMPLEX_IMAG(obj), TRUE);
2139         SCM_PUTZ(buf, -1, SCM_PORT(p));
2140         SCM_PUTC('i', SCM_PORT(p));
2141         r = Scm_GetOutputString(SCM_PORT(p));
2142     } else {
2143         Scm_Error("number required: %S", obj);
2144     }
2145     return r;
2146 }
2147 
2148 /* utility to expose Burger&Dybvig algorithm.  FLAGS is not used yet,
2149    but reserved for future extension. */
2150 void Scm_PrintDouble(ScmPort *port, double d, int flags)
2151 {
2152     char buf[FLT_BUF];
2153     double_print(buf, FLT_BUF, d, FALSE);
2154     Scm_Putz(buf, strlen(buf), port);
2155 }
2156 
2157 /*
2158  * Number Parser
2159  *
2160  *  <number> : <prefix> <complex>
2161  *  <prefix> : <radix> <exactness> | <exactness> <radix>
2162  *  <radix>  : <empty> | '#b' | '#o' | '#d' | '#x'
2163  *  <exactness> : <empty> | '#e' | '#i'
2164  *  <complex> : <real>
2165  *            | <real> '@' <real>
2166  *            | <real> '+' <ureal> 'i'
2167  *            | <real> '-' <ureal> 'i'
2168  *            | <real> '+' 'i'
2169  *            | <real> '-' 'i'
2170  *            | '+' <ureal> 'i'
2171  *            | '-' <ureal> 'i'
2172  *            | '+' 'i'
2173  *            | '-' 'i'
2174  *  <real>   : <sign> <ureal>
2175  *  <sign>   : <empty> | '+' | '-'
2176  *  <ureal>  : <uinteger>
2177  *           | <uinteger> '/' <uinteger>
2178  *           | <decimal>
2179  *  <uinteger> : <digit>+ '#'*
2180  *  <decimal> : <digit10>+ '#'* <suffix>
2181  *            | '.' <digit10>+ '#'* <suffix>
2182  *            | <digit10>+ '.' <digit10>+ '#'* <suffix>
2183  *            | <digit10>+ '#'+ '.' '#'* <suffix>
2184  *  <suffix>  : <empty> | <exponent-marker> <sign> <digit10>+
2185  *  <exponent-marker> : 'e' | 's' | 'f' | 'd' | 'l'
2186  *
2187  * The parser reads characters from on-memory buffer.
2188  * Multibyte strings are filtered out in the early stage of
2189  * parsing, so the subroutines assume the buffer contains
2190  * only ASCII chars.
2191  */
2192 
2193 struct numread_packet {
2194     const char *buffer;         /* original buffer */
2195     int buflen;                 /* original length */
2196     int radix;                  /* radix */
2197     int exactness;              /* exactness; see enum below */
2198     int padread;                /* '#' padding has been read */
2199     int strict;                 /* when true, reports an error if the
2200                                    input violates implementation limitation;
2201                                    otherwise, the routine returns #f. */
2202 };
2203 
2204 enum { /* used in the exactness flag */
2205     NOEXACT, EXACT, INEXACT
2206 };
2207 
2208 /* Max digits D such that all D-digit radix R integers fit in signed
2209    long, i.e. R^(D+1)-1 <= LONG_MAX */
2210 static long longdigs[RADIX_MAX-RADIX_MIN+1] = { 0 };
2211 
2212 /* Max integer I such that reading next digit (in radix R) will overflow
2213    long integer.   floor(LONG_MAX/R - R). */
2214 static u_long longlimit[RADIX_MAX-RADIX_MIN+1] = { 0 };
2215 
2216 /* An integer table of R^D, which is a "big digit" to be added
2217    into bignum. */
2218 static u_long bigdig[RADIX_MAX-RADIX_MIN+1] = { 0 };
2219 
2220 static ScmObj numread_error(const char *msg, struct numread_packet *context);
2221 
2222 /* Returns either small integer or bignum.
2223    initval may be a Scheme integer that will be 'concatenated' before
2224    the integer to be read; it is used to read floating-point number.
2225    Note that value_big may keep denormalized bignum. */
2226 static ScmObj read_uint(const char **strp, int *lenp,
2227                         struct numread_packet *ctx,
2228                         ScmObj initval)
2229 {
2230     const char *str = *strp;
2231     int digread = FALSE;
2232     int len = *lenp;
2233     int radix = ctx->radix;
2234     int digits = 0, diglimit = longdigs[radix-RADIX_MIN];
2235     u_long limit = longlimit[radix-RADIX_MIN], bdig = bigdig[radix-RADIX_MIN];
2236     u_long value_int = 0;
2237     ScmBignum *value_big = NULL;
2238     char c;
2239     static const char tab[] = "0123456789abcdefghijklmnopqrstuvwxyz";
2240     const char *ptab;
2241 
2242     if (!SCM_FALSEP(initval)) {
2243         if (SCM_INTP(initval)) {
2244             if (SCM_INT_VALUE(initval) > limit) {
2245                 value_big = Scm_MakeBignumWithSize(4, SCM_INT_VALUE(initval));
2246             } else {
2247                 value_int = SCM_INT_VALUE(initval);
2248             }
2249         } else if (SCM_BIGNUMP(initval)) {
2250             value_big = SCM_BIGNUM(Scm_BignumCopy(SCM_BIGNUM(initval)));
2251         }
2252         digread = TRUE;
2253     } else if (*str == '0') {
2254         /* Ignore leading 0's, to avoid unnecessary bignum operations. */
2255         while (len > 0 && *str == '0') { str++; len--; }
2256         digread = TRUE;
2257     }
2258 
2259     while (len--) {
2260         int digval = -1;
2261         c = tolower(*str++);
2262         if (ctx->padread) {
2263             if (c == '#') digval = 0;
2264             else break;
2265         } else if (digread && c == '#') {
2266             digval = 0;
2267             ctx->padread = TRUE;
2268             if (ctx->exactness == NOEXACT) {
2269                 ctx->exactness = INEXACT;
2270             }
2271         } else {
2272             for (ptab = tab; ptab < tab+radix; ptab++) {
2273                 if (c == *ptab) {
2274                     digval = ptab-tab;
2275                     digread = TRUE;
2276                     break;
2277                 }
2278             }
2279         }
2280         if (digval < 0) break;
2281         value_int = value_int * radix + digval;
2282         digits++;
2283         if (value_big == NULL) {
2284             if (value_int >= limit) {
2285                 value_big = Scm_MakeBignumWithSize(4, value_int);
2286                 value_int = digits = 0;
2287             }
2288         } else if (digits > diglimit) {
2289             value_big = Scm_BignumAccMultAddUI(value_big, bdig, value_int);
2290             value_int = digits = 0;
2291         }
2292     }
2293     *strp = str-1;
2294     *lenp = len+1;
2295 
2296     if (value_big == NULL) return Scm_MakeInteger(value_int);
2297     if (digits > 0) {
2298         value_big = Scm_BignumAccMultAddUI(value_big, 
2299                                            ipow(radix, digits),
2300                                            value_int);
2301     }
2302     return Scm_NormalizeBignum(SCM_BIGNUM(value_big));
2303 }
2304 
2305 /*
2306  * Find a double number closest to f * 10^e, using z as the starting
2307  * approximation.  The algorithm (and its name) is taken from Will Clinger's
2308  * paper "How to Read Floating Point Numbers Accurately", in the ACM
2309  * SIGPLAN '90, pp.92--101.
2310  * The algorithm is modified to take advantage of coherency between loops.
2311  */
2312 static double algorithmR(ScmObj f, int e, double z)
2313 {
2314     ScmObj m, x, y, abs_d, d2;
2315     int k, s, kprev, sign_d;
2316     m = Scm_DecodeFlonum(z, &k, &s);
2317     IEXPT10_INIT();
2318   retry:
2319     if (k >= 0) {
2320         if (e >= 0) {
2321             x = Scm_Multiply2(f, iexpt10(e));
2322             y = Scm_Ash(m, k);
2323         } else {
2324             x = f;
2325             y = Scm_Ash(Scm_Multiply2(m, iexpt10(-e)), k);
2326         }
2327     } else {
2328         if (e >= 0) {
2329             x = Scm_Ash(Scm_Multiply2(f, iexpt10(e)), -k);
2330             y = m;
2331         } else {
2332             x = Scm_Ash(f, -k);
2333             y = Scm_Multiply2(m, iexpt10(-e));
2334         }
2335     }
2336     kprev = k;
2337 
2338     /* compare  */
2339     for (;;) {
2340         /*Scm_Printf(SCM_CURERR, "z=%.20lg,\nx=%S,\ny=%S\nf=%S\nm=%S\ne=%d, k=%d\n", z, x, y, f, m, e, k);*/
2341         /* compare */
2342         sign_d = Scm_NumCmp(x, y);
2343         abs_d = (sign_d > 0)? Scm_Subtract2(x, y) : Scm_Subtract2(y, x);
2344         d2 = Scm_Ash(Scm_Multiply2(m, abs_d), 1);
2345         switch (Scm_NumCmp(d2, y)) {
2346         case -1: /* d2 < y */
2347             if (Scm_NumCmp(m, SCM_2_52) == 0
2348                 && sign_d < 0
2349                 && Scm_NumCmp(Scm_Ash(d2, 1), y) > 0) {
2350                 goto prevfloat;
2351             } else {
2352                 return ldexp(Scm_GetDouble(m), k);
2353             }
2354         case 0: /* d2 == y */
2355             if (!Scm_OddP(m)) {
2356                 if (Scm_NumCmp(m, SCM_2_52) == 0
2357                     && sign_d < 0) {
2358                     goto prevfloat;
2359                 } else {
2360                     return ldexp(Scm_GetDouble(m), k);
2361                 }
2362             } else if (sign_d < 0) {
2363                 goto prevfloat;
2364             } else {
2365                 goto nextfloat;
2366             }
2367         default:
2368             if (sign_d < 0) goto prevfloat;
2369             else            goto nextfloat;
2370         }
2371       prevfloat:
2372         m = Scm_Subtract2(m, SCM_MAKE_INT(1));
2373         if (k > -1074 && Scm_NumCmp(m, SCM_2_52) < 0) {
2374             m = Scm_Ash(m, 1);
2375             k--;
2376         }
2377         goto next;
2378       nextfloat:
2379         m = Scm_Add2(m, SCM_MAKE_INT(1));
2380         if (Scm_NumCmp(m, SCM_2_53) >= 0) {
2381             m = Scm_Ash(m, -1);
2382             k++;
2383         }
2384         /*FALLTHROUGH*/
2385       next:
2386         if (kprev >= 0) {
2387             if (k >= 0) {
2388                 /* k stays positive. x is invariant */
2389                 if (e >= 0) {
2390                     y = Scm_Ash(m, k);
2391                 } else {
2392                     y = Scm_Ash(Scm_Multiply2(m, iexpt10(-e)), k);
2393                 }
2394             } else {
2395                 /* k turned to negative */
2396                 goto retry;
2397             }
2398         } else {
2399             if (k < 0) {
2400                 /* k stays negative. */
2401                 if (e >= 0) {
2402                     if (k != kprev) x = Scm_Ash(Scm_Multiply2(f, iexpt10(e)), -k);
2403                     y = m;
2404                 } else {
2405                     if (k != kprev) x = Scm_Ash(f, -k);
2406                     y = Scm_Multiply2(m, iexpt10(-e));
2407                 }
2408             } else {
2409                 /* k turned to positive */
2410                 goto retry;
2411             }
2412         }
2413     }
2414     /*NOTREACHED*/
2415 }
2416 
2417 static ScmObj read_real(const char **strp, int *lenp,
2418                         struct numread_packet *ctx)
2419 {
2420     int minusp = FALSE, exp_minusp = FALSE;
2421     int fracdigs = 0;
2422     long exponent = 0;
2423     ScmObj intpart, fraction;
2424 
2425     switch (**strp) {
2426     case '-': minusp = TRUE;
2427         /* FALLTHROUGH */
2428     case '+':
2429         (*strp)++; (*lenp)--;
2430     }
2431     if ((*lenp) <= 0) return SCM_FALSE;
2432 
2433     /* Read integral part */
2434     if (**strp != '.') {
2435         intpart = read_uint(strp, lenp, ctx, SCM_FALSE);
2436         if ((*lenp) <= 0) {
2437             if (minusp) intpart = Scm_Negate(intpart);
2438             if (ctx->exactness == INEXACT) {
2439                 return Scm_ExactToInexact(intpart);
2440             } else {
2441                 return intpart;
2442             }
2443         }
2444         if (**strp == '/') {
2445             /* possibly rational */
2446             ScmObj denom, ratval;
2447             int lensave;
2448             
2449             if ((*lenp) <= 1) return SCM_FALSE;
2450             (*strp)++; (*lenp)--;
2451             lensave = *lenp;
2452             denom = read_uint(strp, lenp, ctx, SCM_FALSE);
2453             if (SCM_FALSEP(denom)) return SCM_FALSE;
2454             if (denom == SCM_MAKE_INT(0)) {
2455                 if (lensave > *lenp) {
2456                     if (minusp) {
2457                         return Scm_MakeFlonum(-1.0/0.0);
2458                     } else {
2459                         return Scm_MakeFlonum(1.0/0.0);
2460                     }
2461                 } else {
2462                     return SCM_FALSE;
2463                 }
2464             }
2465             if (minusp) intpart = Scm_Negate(intpart);
2466             ratval = Scm_Divide2(intpart, denom);
2467 
2468             if (ctx->exactness == EXACT && !Scm_IntegerP(ratval)) {
2469                 return numread_error("(exact non-integral rational number is not supported)",
2470                                      ctx);
2471             }
2472             if (ctx->exactness == INEXACT && !SCM_FLONUMP(ratval)) {
2473                 return Scm_ExactToInexact(ratval);
2474             } else {
2475                 return ratval;
2476             }
2477         }
2478         /* fallthrough */
2479     } else {
2480         intpart = SCM_FALSE; /* indicate there was no intpart */
2481     }
2482 
2483     /* Read fractional part.
2484        At this point, simple integer is already eliminated. */
2485     if (**strp == '.') {
2486         int lensave;
2487         if (ctx->radix != 10) {
2488             return numread_error("(only 10-based fraction is supported)",
2489                                  ctx);
2490         }
2491         (*strp)++; (*lenp)--;
2492         lensave = *lenp;
2493         fraction = read_uint(strp, lenp, ctx, intpart);
2494         fracdigs = lensave - *lenp;
2495     } else {
2496         fraction = intpart;
2497     }
2498 
2499     if (SCM_FALSEP(intpart)) {
2500         if (fracdigs == 0) return SCM_FALSE; /* input was "." */
2501     }
2502 
2503     /* Read exponent.  */
2504     if (*lenp > 0 && strchr("eEsSfFdDlL", (int)**strp)) {
2505         (*strp)++;
2506         if (--(*lenp) <= 0) return SCM_FALSE;
2507         switch (**strp) {
2508         case '-': exp_minusp = TRUE;
2509             /*FALLTHROUGH*/
2510         case '+':
2511             (*strp)++;
2512             if (--(*lenp) <= 0) return SCM_FALSE;
2513         }
2514         while (*lenp > 0) {
2515             int c = **strp;
2516             if (!isdigit(c)) break;
2517             (*strp)++, (*lenp)--;
2518             if (isdigit(c)) {
2519                 exponent = exponent * 10 + (c - '0');
2520                 /* just reject obviously wrong exponent.  more precise
2521                    check will be done later. */
2522                 if (exponent >= LONG_MAX/10 - 10) {
2523                     return numread_error("(exponent of floating-point number out of range)", ctx);
2524                 }
2525             }
2526         }
2527         if (exp_minusp) exponent = -exponent;
2528     }
2529 
2530     /*Scm_Printf(SCM_CURERR, "fraction=%S, exponent=%d\n", fraction, exponent);*/
2531     /* Compose flonum.*/
2532     {
2533         double realnum = Scm_GetDouble(fraction);
2534 
2535         realnum = raise_pow10(realnum, exponent-fracdigs);
2536         if (realnum > 0.0
2537             && (Scm_NumCmp(fraction, SCM_2_52) > 0
2538                 || exponent-fracdigs > MAX_EXACT_10_EXP
2539                 || exponent-fracdigs < -MAX_EXACT_10_EXP)) {
2540             realnum = algorithmR(fraction, exponent-fracdigs, realnum);
2541         }
2542         if (minusp) realnum = -realnum;
2543         /* check exactness */
2544         if (ctx->exactness == EXACT) {
2545             double integ;
2546             if (modf(realnum, &integ) != 0.0) {
2547                 return numread_error("(exact non-integral number is not supported)",
2548                                      ctx);
2549             }
2550             return Scm_InexactToExact(Scm_MakeFlonum(realnum));
2551         }
2552         return Scm_MakeFlonum(realnum);
2553     }
2554 }
2555 
2556 /* Entry point */
2557 static ScmObj read_number(const char *str, int len, int radix, int strict)
2558 {
2559     struct numread_packet ctx;
2560     int radix_seen = 0, exactness_seen = 0, sign_seen = 0;
2561     ScmObj realpart;
2562 
2563     ctx.buffer = str;
2564     ctx.buflen = len;
2565     ctx.exactness = NOEXACT;
2566     ctx.padread = FALSE;
2567     ctx.strict = strict;
2568 
2569 #define CHK_EXACT_COMPLEX()                                                 \
2570     do {                                                                    \
2571         if (ctx.exactness == EXACT) {                                       \
2572             return numread_error("(exact complex number is not supported)", \
2573                                  &ctx);                                     \
2574         }                                                                   \
2575     } while (0)
2576 
2577     /* suggested radix.  may be overridden by prefix. */
2578     if (radix <= 1 || radix > 36) return SCM_FALSE;
2579     ctx.radix = radix;
2580     
2581     /* start from prefix part */
2582     for (; len >= 0; len-=2) {
2583         if (*str != '#') break;
2584         str++;
2585         switch (*str++) {
2586         case 'x':; case 'X':;
2587             if (radix_seen) return SCM_FALSE;
2588             ctx.radix = 16; radix_seen++;
2589             continue;
2590         case 'o':; case 'O':;
2591             if (radix_seen) return SCM_FALSE;
2592             ctx.radix = 8; radix_seen++;
2593             continue;
2594         case 'b':; case 'B':;
2595             if (radix_seen) return SCM_FALSE;
2596             ctx.radix = 2; radix_seen++;
2597             continue;
2598         case 'd':; case 'D':;
2599             if (radix_seen) return SCM_FALSE;
2600             ctx.radix = 10; radix_seen++;
2601             continue;
2602         case 'e':; case 'E':;
2603             if (exactness_seen) return SCM_FALSE;
2604             ctx.exactness = EXACT; exactness_seen++;
2605             continue;
2606         case 'i':; case 'I':;
2607             if (exactness_seen) return SCM_FALSE;
2608             ctx.exactness = INEXACT; exactness_seen++;
2609             continue;
2610         }
2611         return SCM_FALSE;
2612     }
2613     if (len <= 0) return SCM_FALSE;
2614 
2615     /* number body.  need to check the special case of pure imaginary */
2616     if (*str == '+' || *str == '-') {
2617         if (len == 1) return SCM_FALSE;
2618         if (len == 2 && (str[1] == 'i' || str[1] == 'I')) {
2619             CHK_EXACT_COMPLEX();
2620             return Scm_MakeComplex(0.0, (*str == '+')? 1.0 : -1.0);
2621         }
2622         sign_seen = TRUE;
2623     }
2624 
2625     realpart = read_real(&str, &len, &ctx);
2626     if (SCM_FALSEP(realpart) || len == 0) return realpart;
2627 
2628     switch (*str) {
2629     case '@':
2630         /* polar representation of complex*/
2631         if (len <= 1) {
2632             return SCM_FALSE;
2633         } else {
2634             ScmObj angle;
2635             double dmag, dangle;
2636             str++; len--;
2637             angle = read_real(&str, &len, &ctx);
2638             if (SCM_FALSEP(angle) || len != 0) return SCM_FALSE;
2639             CHK_EXACT_COMPLEX();
2640             dmag = Scm_GetDouble(realpart);
2641             dangle = Scm_GetDouble(angle);
2642             return Scm_MakeComplexPolar(dmag, dangle);
2643         }
2644     case '+':;
2645     case '-':
2646         /* rectangular representation of complex */
2647         if (len <= 1) {
2648             return SCM_FALSE;
2649         } else if (len == 2 && str[1] == 'i') {
2650             return Scm_MakeComplex(Scm_GetDouble(realpart),
2651                                    (*str == '+' ? 1.0 : -1.0));
2652         } else {
2653             ScmObj imagpart = read_real(&str, &len, &ctx);
2654             if (SCM_FALSEP(imagpart) || len != 1 || *str != 'i') {
2655                 return SCM_FALSE;
2656             }
2657             CHK_EXACT_COMPLEX();
2658             if (Scm_Sign(imagpart) == 0) return realpart;
2659             return Scm_MakeComplexNormalized(Scm_GetDouble(realpart), 
2660                                              Scm_GetDouble(imagpart));
2661         }
2662     case 'i':
2663         /* '+' <ureal> 'i'  or '-' <ureal> 'i' */
2664         if (!sign_seen || len != 1) return SCM_FALSE;
2665         CHK_EXACT_COMPLEX();
2666         if (Scm_Sign(realpart) == 0) return Scm_MakeFlonum(0.0);
2667         else return Scm_MakeComplex(0.0, Scm_GetDouble(realpart));
2668     default:
2669         return SCM_FALSE;
2670     }
2671 }
2672 
2673 static ScmObj numread_error(const char *msg, struct numread_packet *context)
2674 {
2675     if (context->strict) {
2676         Scm_Error("bad number format %s: %A", msg,
2677                   Scm_MakeString(context->buffer, context->buflen,
2678                                  context->buflen, 0));
2679     }
2680     return SCM_FALSE;
2681 }
2682 
2683 
2684 ScmObj Scm_StringToNumber(ScmString *str, int radix, int strict)
2685 {
2686     u_int len, size;
2687     const char *p = Scm_GetStringContent(str, &size, &len, NULL);
2688     if (size != len) {
2689         /* This can't be a proper number. */
2690         return SCM_FALSE;
2691     } else {
2692         return read_number(p, size, radix, strict);
2693     }
2694 }
2695 
2696 /*
2697  * Initialization
2698  */
2699 
2700 ScmObj Scm__ConstObjs[SCM_NUM_CONST_OBJS] = { SCM_FALSE };
2701 
2702 void Scm__InitNumber(void)
2703 {
2704     ScmModule *mod = Scm_GaucheModule();
2705     int radix, i;
2706     u_long n;
2707     
2708     for (radix = RADIX_MIN; radix <= RADIX_MAX; radix++) {
2709         longlimit[radix-RADIX_MIN] =
2710             (u_long)floor((double)LONG_MAX/radix - radix);
2711         /* Find max D where R^(D+1)-1 <= LONG_MAX */
2712         for (i = 0, n = 1; ; i++, n *= radix) {
2713             if (n >= LONG_MAX/radix) {
2714                 longdigs[radix-RADIX_MIN] = i-1;
2715                 bigdig[radix-RADIX_MIN] = n;
2716                 break;
2717             }
2718         }
2719     }
2720     
2721     SCM_2_63 = Scm_Ash(SCM_MAKE_INT(1), 63);
2722     SCM_2_64 = Scm_Ash(SCM_MAKE_INT(1), 64);
2723     SCM_2_64_MINUS_1 = Scm_Subtract2(SCM_2_64, SCM_MAKE_INT(1));
2724     SCM_2_52 = Scm_Ash(SCM_MAKE_INT(1), 52);
2725     SCM_2_53 = Scm_Ash(SCM_MAKE_INT(1), 53);
2726     SCM_MINUS_2_63 = Scm_Negate(SCM_2_63);
2727     SCM_2_32 = Scm_Ash(SCM_MAKE_INT(1), 32);
2728     SCM_2_31 = Scm_Ash(SCM_MAKE_INT(1), 31);
2729     SCM_MINUS_2_31 = Scm_Negate(SCM_2_31);
2730     
2731     dexpt2_minus_52 = ldexp(1.0, -52);
2732     dexpt2_minus_53 = ldexp(1.0, -53);
2733 
2734     Scm_InitBuiltinGeneric(&generic_add, "object-+", mod);
2735     Scm_InitBuiltinGeneric(&generic_sub, "object--", mod);
2736     Scm_InitBuiltinGeneric(&generic_mul, "object-*", mod);
2737     Scm_InitBuiltinGeneric(&generic_div, "object-/", mod);
2738 }

/* [<][>][^][v][top][bottom][index][help] */