1 /*
2 * number.c - numeric functions
3 *
4 * Copyright (c) 2000-2020 Shiro Kawai <shiro@acm.org>
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 #define LIBGAUCHE_BODY
34 #include "gauche.h"
35 #include "gauche/bignum.h"
36 #include "gauche/scmconst.h"
37 #include "gauche/bits.h"
38 #include "gauche/bits_inline.h"
39 #include "gauche/priv/builtin-syms.h"
40 #include "gauche/priv/arith.h"
41 #include "gauche/priv/bytesP.h"
42
43 #include <limits.h>
44 #include <float.h>
45 #include <math.h>
46 #include <ctype.h>
47
48 /*================================================================
49 * Some macros
50 */
51
52 #ifdef HAVE_SUNMATH_H
53 #include "sunmath.h" /* for isinf() on Solaris */
54 #endif /* HAVE_SUNMATH_H */
55
56 #ifndef M_PI
57 #define M_PI 3.14159265358979323846
58 #endif
59
60 #ifndef HAVE_ISINF
61 /* NB: If we inline this, some version of gcc incorrectly assumes
62 the condition would never be satisfied and optimize it away. */
Scm_IsInf(double x)63 int Scm_IsInf(double x)
64 {
65 volatile double xx = x; /* suppress gcc to be too clever */
66 return ((xx) != 0 && (xx) == (xx)/2.0);
67 }
68 #endif
69
70 /* Maximum allowable range of exponent in the number litereal.
71 For flonums, IEEE double can support [-324..308]. For exact
72 numbers we can go further, but it would easily consume huge
73 memory. So I assume it is reasonable to limit its range. */
74 #define MAX_EXPONENT 325
75
76 /* Linux gcc have those, but the declarations aren't included unless
77 __USE_ISOC9X is defined. Just in case. */
78 #ifdef HAVE_TRUNC
79 extern double trunc(double);
80 #endif
81
82 #ifdef HAVE_RINT
83 extern double rint(double);
84 #define roundeven rint
85 #else
86 static double roundeven(double);
87 #endif
88
89 /* Many built-in arithmetic routines come with VM* variants, which
90 are used when the resulting flonum is directly returned to the VM.
91 Such routines share the body with their variants, with an extra
92 flag 'vmp' that indicates whether the call is 'VM' variant or not.
93 RETURN_FLONUM and RETURN_FLOBJ macros are used to hide the different
94 flonum generation in such routines.
95 RETURN_FLONUM takes double arg and returns ScmObj.
96 RETURN_FLOBJ takes ScmObj (ScmFlonum*) and returns ScmObj.
97 */
98 #if GAUCHE_FFX
99 #define RETURN_FLONUM(z) \
100 do { \
101 if (vmp) return Scm_VMReturnFlonum(z); \
102 else return Scm_MakeFlonum(z); \
103 } while (0)
104 #else /*!GAUCHE_FFX*/
105 #define RETURN_FLONUM(z) return Scm_MakeFlonum(z)
106 #endif /*!GAUCHE_FFX*/
107
108 #define DEFINE_DUAL_API1(a, b, kernel) \
109 ScmObj a(ScmObj obj) { return kernel(obj, FALSE); } \
110 ScmObj b(ScmObj obj) { return kernel(obj, TRUE); }
111
112 #define DEFINE_DUAL_API2(a, b, kernel) \
113 ScmObj a(ScmObj obj1, ScmObj obj2) { return kernel(obj1, obj2, FALSE); } \
114 ScmObj b(ScmObj obj1, ScmObj obj2) { return kernel(obj1, obj2, TRUE); }
115
116 /*================================================================
117 * Classes of Numeric Tower
118 */
119
120 static ScmClass *numeric_cpl[] = {
121 SCM_CLASS_STATIC_PTR(Scm_RationalClass),
122 SCM_CLASS_STATIC_PTR(Scm_RealClass),
123 SCM_CLASS_STATIC_PTR(Scm_ComplexClass),
124 SCM_CLASS_STATIC_PTR(Scm_NumberClass),
125 SCM_CLASS_STATIC_PTR(Scm_TopClass),
126 NULL
127 };
128
129 static void number_print(ScmObj obj, ScmPort *port, ScmWriteContext *ctx);
130
131 SCM_DEFINE_BUILTIN_CLASS(Scm_NumberClass, number_print, NULL, NULL, NULL,
132 numeric_cpl+4);
133 SCM_DEFINE_BUILTIN_CLASS(Scm_ComplexClass, number_print, NULL, NULL, NULL,
134 numeric_cpl+3);
135 SCM_DEFINE_BUILTIN_CLASS(Scm_RealClass, number_print, NULL, NULL, NULL,
136 numeric_cpl+2);
137 SCM_DEFINE_BUILTIN_CLASS(Scm_RationalClass, number_print, NULL, NULL, NULL,
138 numeric_cpl+1);
139 SCM_DEFINE_BUILTIN_CLASS(Scm_IntegerClass, number_print, NULL, NULL, NULL,
140 numeric_cpl);
141
142 /*=====================================================================
143 * Generic Arithmetic
144 */
145
146 /* Some arithmetic operations calls the corresponding generic function
147 * if the operand is not a number.
148 */
149
150 /* Fallback Gf */
bad_number_method(ScmObj * args,int nargs,ScmGeneric * gf)151 static ScmObj bad_number_method(ScmObj *args, int nargs, ScmGeneric *gf)
152 {
153 const char *fn = (const char *)SCM_GENERIC_DATA(gf);
154 if (nargs == 1) {
155 Scm_Error("operation %s is not defined on object %S", fn, args[0]);
156 } else if (nargs == 2) {
157 Scm_Error("operation %s is not defined between %S and %S",
158 fn, args[0], args[1]);
159 } else {
160 Scm_Error("generic function for %s is called with args %S",
161 fn, Scm_ArrayToList(args, nargs));
162 }
163 return SCM_UNDEFINED;
164 }
165 static SCM_DEFINE_GENERIC(generic_add, bad_number_method, "+");
166 static SCM_DEFINE_GENERIC(generic_sub, bad_number_method, "-");
167 static SCM_DEFINE_GENERIC(generic_mul, bad_number_method, "*");
168 static SCM_DEFINE_GENERIC(generic_div, bad_number_method, "/");
169
170 /*=====================================================================
171 * IEEE754 double and Endianness
172 */
173
174 /* Structure to extract bits from a double. This info may be provided
175 * by a system header (e.g. ieee754.h) but for the portability we
176 * define by ourselves.
177 */
178 typedef union {
179 double d;
180 struct {
181 #ifdef WORDS_BIGENDIAN
182 #if SIZEOF_LONG >= 8
183 unsigned int sign:1;
184 unsigned int exp:11;
185 unsigned long mant:52;
186 #else /*SIZEOF_LONG < 8*/
187 unsigned int sign:1;
188 unsigned int exp:11;
189 unsigned long mant0:20;
190 unsigned long mant1:32;
191 #endif /*SIZEOF_LONG < 8*/
192 #else /*!WORDS_BIGENDIAN*/
193 #if SIZEOF_LONG >= 8
194 unsigned long mant:52;
195 unsigned int exp:11;
196 unsigned int sign:1;
197 #else /*SIZEOF_LONG < 8*/
198 unsigned long mant1:32;
199 unsigned long mant0:20;
200 unsigned int exp:11;
201 unsigned int sign:1;
202 #endif /*SIZEOF_LONG < 8*/
203 #endif /*!WORDS_BIGENDIAN*/
204 } components;
205 } ScmIEEEDouble;
206
207 #ifdef DOUBLE_ARMENDIAN
208 /* ARM processor may be configured to use a special mixed endian.
209 We check at runtime. */
210 typedef union {
211 double d;
212 struct {
213 unsigned long mant0:20;
214 unsigned int exp:11;
215 unsigned int sign:1;
216 unsigned long mant1:32;
217 } components;
218 } ScmIEEEDoubleARM;
219 #endif /*DOUBLE_ARMENDIAN*/
220
221 /* ARM special handling */
222 #ifdef DOUBLE_ARMENDIAN
223 static int armendian_p = FALSE;
224
225 #define TEST_DBL 1.9999999999999998 /* all '1' bits for mantissa */
226
check_armendian()227 void check_armendian()
228 {
229 ScmIEEEDouble z;
230 z.d = TEST_DBL;
231 if (z.components.exp != 1023) {
232 ScmIEEEDoubleARM z2;
233 z2.d = TEST_DBL;
234 if (z2.components.exp != 1023) {
235 Scm_Panic("Initiaization failed: Cannot determine double's endian "
236 "on this ARM processor.");
237 }
238 armendian_p = TRUE;
239 } else {
240 armendian_p = FALSE;
241 }
242 }
243 #endif /*DOUBLE_ARMENDIAN*/
244
245 /* This is lazyily initialized, for at the time number.c is initialized,
246 not sufficient infrastructure is booted to initialize a parameter.
247 */
248 static ScmPrimitiveParameter *default_endian = NULL;
249 static ScmInternalMutex default_endian_mutex;
250
init_default_endian()251 static void init_default_endian()
252 {
253 if (default_endian == NULL) {
254 SCM_INTERNAL_MUTEX_LOCK(default_endian_mutex);
255 if (default_endian == NULL) {
256 default_endian =
257 Scm_BindPrimitiveParameter(Scm_GaucheModule(), "default-endian",
258 SCM_OBJ(Scm_NativeEndian()), 0);
259
260 }
261 SCM_INTERNAL_MUTEX_UNLOCK(default_endian_mutex);
262 }
263 }
264
Scm_NativeEndian()265 ScmObj Scm_NativeEndian()
266 {
267 #ifdef DOUBLE_ARMENDIAN
268 if (armendian_p) return SCM_SYM_ARM_LITTLE_ENDIAN;
269 #endif /*DOUBLE_ARMENDIAN*/
270 #if WORDS_BIGENDIAN
271 return SCM_SYM_BIG_ENDIAN;
272 #else
273 return SCM_SYM_LITTLE_ENDIAN;
274 #endif
275 }
276
Scm_DefaultEndian(void)277 ScmObj Scm_DefaultEndian(void)
278 {
279 init_default_endian();
280 return Scm_PrimitiveParameterRef(Scm_VM(), default_endian);
281 }
282
Scm_SetDefaultEndian(ScmObj endian)283 void Scm_SetDefaultEndian(ScmObj endian)
284 {
285 /* We trust the caller passes one of symbols big-endian, little-endian
286 or arm-little-endian. */
287 init_default_endian();
288 Scm_PrimitiveParameterSet(Scm_VM(), default_endian, endian);
289 }
290
Scm_IsBE(ScmObj endian)291 int Scm_IsBE(ScmObj endian)
292 {
293 return SCM_IS_BE(endian);
294 }
295
Scm_IsLE(ScmObj endian)296 int Scm_IsLE(ScmObj endian)
297 {
298 return SCM_IS_LE(endian);
299 }
300
Scm_IsArmLE(ScmObj endian)301 int Scm_IsArmLE(ScmObj endian)
302 {
303 return SCM_IS_ARM_LE(endian);
304 }
305
306
307 /*=====================================================================
308 * Flonums
309 */
310
311 #undef COUNT_FLONUM_ALLOC
312
313 #ifdef COUNT_FLONUM_ALLOC /* for benchmarks. usually this should be off. */
314 static u_long flonum_count = 0;
315
report_flonum_count(void * data)316 static void report_flonum_count(void *data)
317 {
318 fprintf(stderr, "allocated flonums = %8d\n", flonum_count);
319 }
320 #endif /*COUNT_FLONUM_ALLOC*/
321
Scm_MakeFlonum(double d)322 ScmObj Scm_MakeFlonum(double d)
323 {
324 ScmFlonum *f = SCM_NEW(ScmFlonum);
325 SCM_FLONUM_VALUE(f) = d;
326 #ifdef COUNT_FLONUM_ALLOC
327 flonum_count++;
328 #endif
329 return SCM_MAKE_FLONUM_MEM(f);
330 }
331
Scm_FlonumIntegerToExact(double d)332 ScmObj Scm_FlonumIntegerToExact(double d) /* d mustn't have fractional part */
333 {
334 #if SIZEOF_LONG >= 8
335 /* On 64bit machine, double can't exactly represent SCM_SMALL_INT_MIN and
336 SCM_SMALL_INT_MAX, and comparing d with them could pass through
337 out-of-range value, so we convert d to long first. */
338 if (LONG_MIN <= d && d <= LONG_MAX) {
339 long n = (long)d;
340 if (SCM_SMALL_INT_MIN <= n && n <= SCM_SMALL_INT_MAX) {
341 return SCM_MAKE_INT(n);
342 }
343 }
344 /* FALLTHROUGH */
345 #else
346 /* On 32bit machine, double has enough precision to cover small int
347 range. */
348 if (SCM_SMALL_INT_MIN <= d && d <= SCM_SMALL_INT_MAX) {
349 return SCM_MAKE_INT((long)d);
350 }
351 #endif
352 return Scm_MakeBignumFromDouble(d);
353 }
354
Scm_MakeFlonumToNumber(double d,int exact)355 ScmObj Scm_MakeFlonumToNumber(double d, int exact)
356 {
357 if (exact && !SCM_IS_INF(d)) {
358 /* see if d can be demoted to integer */
359 double i, f;
360 f = modf(d, &i);
361 if (f == 0.0) return Scm_FlonumIntegerToExact(i);
362 /*FALLTHROUGH*/
363 }
364 return Scm_MakeFlonum(d);
365 }
366
367 /*
368 * Flonum decomposition
369 */
370
decode_double(double d,u_long * mant1 SCM_UNUSED,u_long * mant0,int * exp,int * sign)371 static inline void decode_double(double d,
372 u_long *mant1 SCM_UNUSED,
373 u_long *mant0,
374 int *exp, int *sign)
375 {
376 ScmIEEEDouble dd;
377 #ifdef DOUBLE_ARMENDIAN /* ARM-specific handling */
378 ScmIEEEDoubleARM dd2;
379 if (armendian_p) {
380 dd2.d = d;
381 *mant1 = (u_int)dd2.components.mant1;
382 *mant0 = (u_int)dd2.components.mant0;
383 *exp = dd2.components.exp;
384 *sign = dd2.components.sign;
385 return;
386 }
387 #endif /*DOUBLE_ARMENDIAN*/
388
389 dd.d = d;
390 #if SIZEOF_LONG >= 8
391 *mant0 = dd.components.mant;
392 *exp = dd.components.exp;
393 *sign = dd.components.sign;
394 #else /* SIZEOF_LONG == 4 */
395 *mant1 = (u_int)dd.components.mant1;
396 *mant0 = (u_int)dd.components.mant0;
397 *exp = dd.components.exp;
398 *sign = dd.components.sign;
399 #endif /* SIZEOF_LONG == 4 */
400 }
401
402 /* Decompose flonum D into an integer mantissa F and exponent E, where
403 * -1074 <= E <= 971
404 * 0 <= abs(F) < 2^53
405 * D = F * 2^E
406 * Some special cases:
407 * F = 0, E = 0 if D = 0.0 or -0.0
408 * F = #t if D is infinity (positive or negative)
409 * F = #f if D is NaN.
410 * If D is normalized number (e.g. E > -1074), F >= 2^52.
411 * If D is denormalized number, E == -1074 and F < 2^52
412 * Cf. IEEE 754 Reference
413 * http://babbage.cs.qc.edu/courses/cs341/IEEE-754references.html
414 */
Scm_DecodeFlonum(double d,int * exp,int * sign)415 ScmObj Scm_DecodeFlonum(double d, int *exp, int *sign)
416 {
417 ScmObj f;
418 u_long mant1, mant0;
419 int exp0, sign0;
420
421 decode_double(d, &mant1, &mant0, &exp0, &sign0);
422
423 *sign = (sign0? -1 : 1);
424
425 /* Check exceptional cases */
426 if (exp0 == 0x7ff) {
427 *exp = 0;
428 if (
429 #if SIZEOF_LONG >= 8
430 mant0 == 0
431 #else /*SIZEOF_LONG < 8*/
432 mant0 == 0 && mant1 == 0
433 #endif /*SIZEOF_LONG < 8*/
434 ) {
435 return SCM_TRUE; /* infinity */
436 } else {
437 return SCM_FALSE; /* NaN */
438 }
439 }
440
441 *exp = (exp0? exp0 - 0x3ff - 52 : -0x3fe - 52);
442
443 #if SIZEOF_LONG >= 8
444 {
445 if (exp0 > 0) mant0 += (1L<<52); /* hidden bit */
446 f = Scm_MakeInteger(mant0);
447 }
448 #else /*SIZEOF_LONG < 8*/
449 {
450 u_long values[2];
451 values[0] = mant1;
452 values[1] = mant0;
453 if (exp0 > 0) values[1] += (1L<<20); /* hidden bit */
454 f = Scm_NormalizeBignum(SCM_BIGNUM(Scm_MakeBignumFromUIArray(1, values, 2)));
455 }
456 #endif /*SIZEOF_LONG < 8*/
457 return f;
458 }
459
460 /* returns -1 or 1. Scm_Sign cannot distinguish 0.0 and -0.0. This one can.
461 signbit(3) is in C99. */
Scm_FlonumSign(double d)462 int Scm_FlonumSign(double d)
463 {
464 return signbit(d)? -1 : 1;
465 }
466
467 /* Half float support */
468
Scm_HalfToDouble(ScmHalfFloat v)469 double Scm_HalfToDouble(ScmHalfFloat v)
470 {
471 int e = SCM_HALF_FLOAT_EXPONENT(v);
472 int m = SCM_HALF_FLOAT_MANTISSA(v);
473 int s = SCM_HALF_FLOAT_SIGN_BIT(v);
474 if (e == 31) { /* special */
475 if (m == 0) {
476 if (s) return SCM_DBL_NEGATIVE_INFINITY;
477 else return SCM_DBL_POSITIVE_INFINITY;
478 } else {
479 return SCM_DBL_NAN;
480 }
481 }
482 if (e > 0) { /* normalized */
483 double d = ldexp(1.0 + m/1024.0, e - 15);
484 return s? -d : d;
485 }
486 else { /* denormalized */
487 double d = ldexp(m/1024.0, -14);
488 return s? -d : d;
489 }
490 }
491
Scm_DoubleToHalf(double v)492 ScmHalfFloat Scm_DoubleToHalf(double v)
493 {
494 u_long mant1, mant0;
495 int exp0, sign0;
496
497 decode_double(v, &mant1, &mant0, &exp0, &sign0);
498
499 if (exp0 == 0x7ff) { /* special */
500 if (
501 #if SIZEOF_LONG >= 8
502 mant0 == 0
503 #else /*SIZEOF_LONG < 8*/
504 mant0 == 0 && mant1 == 0
505 #endif /*SIZEOF_LONG < 8*/
506 ) {
507 return sign0? 0xfc00 : 0x7c00;
508 } else {
509 return 0x7fff;
510 }
511 }
512 int e = exp0 - 1023 + 15;
513 if (e >= 31) { /* overflow */
514 return sign0? 0xfc00 : 0x7c00;
515 }
516 /* Calculate required mantissa bits. We need upper 10 bits, unless
517 e < 0, in which case we get denormalized number. */
518 int mbits = 10 + ((e <= 0)? e-1 : 0);
519 if (mbits < -1) { /* underflow (-1 for rounding, see below) */
520 return sign0? 0x8000 : 0x0000;
521 }
522 if (e < 0) e = 0;
523 /* Take the mantissa bits. We take one extra bit to perform
524 roudning. R is used to determine whether lower bits are
525 all 0 or not. */
526 #if SIZEOF_LONG >= 8
527 unsigned long m = mant0 >> (52-mbits-1);
528 unsigned long r = mant0 & ((1UL << (52-mbits-1)) - 1);
529 #else /*SIZEOF_LONG < 8*/
530 unsigned long m = mant0 >> (20-mbits-1);
531 unsigned long r = (mant0 & ((1UL << (20-mbits-1)) - 1))|mant1;
532 #endif /*SIZEOF_LONG < 8*/
533 m += 1<<(mbits+1); /* recover hidden bit */
534
535 if (m%2 == 1) {
536 if (r == 0) {
537 /* half point. we round to even */
538 if (m&2) m += 2;
539 } else {
540 m += 2;
541 }
542 }
543
544 /* drop rounding bits */
545 m >>= 1;
546 if (m >= 0x800) {
547 e += 1;
548 m >>= 1;
549 }
550 if (e == 0 && m >= 0x400) {
551 e += 1;
552 m &= ~0x400;
553 }
554 if (e >= 31) { /* overflow by rounding */
555 return sign0? 0xfc00 : 0x7c00;
556 }
557 /* at this point, normalized numbers should get
558 0x400 <= m <= 0x7ff, e > 0, and denormalized numbers should get
559 0 <= m <= 0x3ff, e == 0. So we don't need to treat denormalized
560 specially. */
561 return (ScmHalfFloat)((sign0? 0x8000 : 0x0000)|(e << 10)|(m & 0x3ff));
562 }
563
564 /* Construct a double directly from the given bit patterns. This is
565 an internal procedure; external procedure should use Scm_EncodeFlonum.
566
567 On 64bit architecture, only mant1 is used for mantissa.
568 On 32bit architecture, mant1 is for lower 32bits of mantissa, and
569 lower 20bits of mant0 is used for higher bits.
570 */
Scm__EncodeDouble(u_long mant1,u_long mant0 SCM_UNUSED,int exp,int signbit)571 double Scm__EncodeDouble(u_long mant1,
572 u_long mant0 SCM_UNUSED,
573 int exp, int signbit)
574 {
575 ScmIEEEDouble dd;
576 #ifdef DOUBLE_ARMENDIAN
577 ScmIEEEDoubleARM dd2;
578 if (armendian_p) {
579 dd2.components.mant1 = mant1;
580 dd2.components.mant0 = mant0;
581 dd2.components.exp = exp;
582 dd2.components.sign = signbit;
583 return dd2.d;
584 }
585 #endif /*DOUBLE_ARMENDIAN*/
586
587 dd.components.exp = exp;
588 dd.components.sign = signbit;
589 #if SIZEOF_LONG >= 8
590 dd.components.mant = mant1;
591 #else /*SIZEOF_LONG==4*/
592 dd.components.mant1 = mant1;
593 dd.components.mant0 = mant0;
594 #endif /*SIZEOF_LONG==4*/
595 return dd.d;
596 }
597
598 /* Inverse of Scm_DecodeFlonum.
599 Returns the double representation of S * F * 2^E,
600 where S = 1 | -1,
601 -1074 <= E <= 971
602 2^52 <= F < 2^53 (if E > -1074)
603 or 0 <= F < 2^52 (if E = -1074)
604 NB: If E < -1074, we scale F accordingly until E becomes -1074.
605 */
Scm_EncodeFlonum(ScmObj mant,int exp,int sign)606 double Scm_EncodeFlonum(ScmObj mant, int exp, int sign)
607 {
608 if (SCM_FALSEP(mant)) return SCM_DBL_NAN;
609 if (SCM_TRUEP(mant)) {
610 if (sign < 0) return SCM_DBL_NEGATIVE_INFINITY;
611 else return SCM_DBL_POSITIVE_INFINITY;
612 }
613
614 int signbit = sign < 0? 1 : 0;
615 if (exp < -1074) {
616 /* scale mantissa */
617 int shift = -1074 - exp;
618 mant = Scm_Ash(mant, -shift);
619 exp = -1074;
620 } else if (exp > 971) {
621 Scm_Error("flonum exponent out of range: %d", exp);
622 }
623
624 uint64_t mant64 = Scm_GetIntegerU64Clamp(mant, SCM_CLAMP_ERROR, NULL);
625 if (!Scm_NumLT(mant, SCM_2_53)) {
626 Scm_Error("flonum mantissa out of range: %S", mant);
627 }
628 if (exp == -1074 && Scm_NumLT(mant, SCM_2_52)) {
629 exp -= 1; /* denormalized range */
630 } else if (Scm_NumLT(mant, SCM_2_52)) {
631 Scm_Error("flonum mantissa out of range: %S", mant);
632 }
633
634 int expfield = exp + 0x3ff + 52;
635 #if SIZEOF_LONG >= 8
636 return Scm__EncodeDouble(mant64, 0, expfield, signbit);
637 #else
638 u_long hi = (mant64 >> 32);
639 u_long lo = (u_long)(mant64 & ULONG_MAX);
640 return Scm__EncodeDouble(lo, hi, expfield, signbit);
641 #endif
642 }
643
644 /*=====================================================================
645 * Ratnums
646 */
647
648 /* possibly returns denomalized number */
Scm_MakeRatnum(ScmObj numer,ScmObj denom)649 ScmObj Scm_MakeRatnum(ScmObj numer, ScmObj denom)
650 {
651 if (!SCM_INTEGERP(numer)) {
652 Scm_Error("numerator must be an exact integer, but got %S", numer);
653 }
654 if (!SCM_INTEGERP(denom)) {
655 Scm_Error("denominator must be an exact integer, but got %S", denom);
656 }
657 if (SCM_EXACT_ZERO_P(denom)) {
658 Scm_Error("attempt to calculate a division by zero");
659 }
660 ScmRatnum *r = SCM_NEW(ScmRatnum);
661 SCM_SET_CLASS(r, SCM_CLASS_RATIONAL);
662 r->numerator = numer;
663 r->denominator = denom;
664 return SCM_OBJ(r);
665 }
666
667 #define ENSURE_RATNUM(integer) \
668 SCM_RATNUM(Scm_MakeRatnum(integer, SCM_MAKE_INT(1)))
669
Scm_MakeRational(ScmObj numer,ScmObj denom)670 ScmObj Scm_MakeRational(ScmObj numer, ScmObj denom)
671 {
672 if (!SCM_INTEGERP(numer)) {
673 Scm_Error("numerator must be an exact integer, but got %S", numer);
674 }
675 if (!SCM_INTEGERP(denom)) {
676 Scm_Error("denominator must be an exact integer, but got %S", denom);
677 }
678 if (SCM_EXACT_ZERO_P(denom)) {
679 Scm_Error("attempt to calculate a division by zero");
680 }
681 if (SCM_EXACT_ONE_P(denom)) return numer;
682 if (SCM_EXACT_ZERO_P(numer)) return SCM_MAKE_INT(0);
683 else return Scm_ReduceRational(Scm_MakeRatnum(numer, denom));
684 }
685
Scm_Numerator(ScmObj n)686 ScmObj Scm_Numerator(ScmObj n)
687 {
688 if (SCM_RATNUMP(n)) return SCM_RATNUM_NUMER(n);
689 if (SCM_EXACTP(n)) return n; /* fixnum or bignum */
690 if (!SCM_REALP(n)) SCM_TYPE_ERROR(n, "real number");
691 return Scm_Inexact(Scm_Numerator(Scm_Exact(n)));
692 }
693
Scm_Denominator(ScmObj n)694 ScmObj Scm_Denominator(ScmObj n)
695 {
696 if (SCM_RATNUMP(n)) return SCM_RATNUM_DENOM(n);
697 if (SCM_INTEGERP(n)) return SCM_MAKE_INT(1);
698 if (!SCM_REALP(n)) SCM_TYPE_ERROR(n, "real number");
699 return Scm_Inexact(Scm_Denominator(Scm_Exact(n)));
700 }
701
Scm_ReduceRational(ScmObj rational)702 ScmObj Scm_ReduceRational(ScmObj rational)
703 {
704 int negated = FALSE;
705
706 if (SCM_INTEGERP(rational)) return rational;
707 if (!SCM_RATNUMP(rational)) {
708 Scm_Error("exact rational number required, but got %S", rational);
709 }
710 ScmObj numer = SCM_RATNUM_NUMER(rational);
711 ScmObj denom = SCM_RATNUM_DENOM(rational);
712
713 if (Scm_Sign(denom) < 0) {
714 numer = Scm_Negate(numer);
715 denom = Scm_Negate(denom);
716 negated = TRUE;
717 }
718
719 /* special cases */
720 if (SCM_EXACT_ONE_P(denom)) return numer;
721 if (SCM_EXACT_ZERO_P(denom)) {
722 int s = Scm_Sign(numer);
723 if (s > 0) return SCM_POSITIVE_INFINITY;
724 if (s < 0) return SCM_NEGATIVE_INFINITY;
725 return SCM_NAN;
726 }
727
728 ScmObj common = Scm_Gcd(numer, denom);
729 if (SCM_EXACT_ONE_P(common)) {
730 if (negated) {
731 return Scm_MakeRatnum(numer, denom);
732 } else {
733 return rational;
734 }
735 } else {
736 numer = Scm_Quotient(numer, common, NULL);
737 denom = Scm_Quotient(denom, common, NULL);
738 if (SCM_EQ(denom, SCM_MAKE_INT(1))) {
739 return numer;
740 } else {
741 return Scm_MakeRatnum(numer, denom);
742 }
743 }
744 }
745
746 /* x, y must be exact numbers */
Scm_RatnumAddSub(ScmObj x,ScmObj y,int subtract)747 ScmObj Scm_RatnumAddSub(ScmObj x, ScmObj y, int subtract)
748 {
749 ScmObj nx = SCM_RATNUMP(x)? SCM_RATNUM_NUMER(x) : x;
750 ScmObj dx = SCM_RATNUMP(x)? SCM_RATNUM_DENOM(x) : SCM_MAKE_INT(1);
751 ScmObj ny = SCM_RATNUMP(y)? SCM_RATNUM_NUMER(y) : y;
752 ScmObj dy = SCM_RATNUMP(y)? SCM_RATNUM_DENOM(y) : SCM_MAKE_INT(1);
753 ScmObj gcd, dr, nr;
754
755 /* shortcut */
756 if (Scm_NumEq(dx, dy)) {
757 dr = dx;
758 goto finish;
759 }
760
761 if (SCM_EXACT_ONE_P(dx)||SCM_EXACT_ONE_P(dx)) gcd = SCM_MAKE_INT(1);
762 else gcd = Scm_Gcd(dx, dy);
763 if (Scm_NumEq(dx, gcd)) {
764 /* only factor x */
765 nx = Scm_Mul(Scm_Quotient(dy, dx, NULL), nx);
766 dr = dy;
767 goto finish;
768 }
769 if (Scm_NumEq(dy, gcd)) {
770 /* only factor y */
771 ny = Scm_Mul(Scm_Quotient(dx, dy, NULL), ny);
772 dr = dx;
773 goto finish;
774 }
775
776 /* general case */
777 ScmObj fx = Scm_Quotient(dx, gcd, NULL);
778 ScmObj fy = Scm_Quotient(dy, gcd, NULL);
779 nx = Scm_Mul(nx, fy);
780 ny = Scm_Mul(ny, fx);
781 dr = Scm_Mul(dx, fy);
782 finish:
783 nr = (subtract? Scm_Sub(nx, ny) : Scm_Add(nx, ny));
784 return Scm_MakeRational(nr, dr);
785 }
786
Scm_RatnumMulDiv(ScmObj x,ScmObj y,int divide)787 ScmObj Scm_RatnumMulDiv(ScmObj x, ScmObj y, int divide)
788 {
789 ScmObj nx = SCM_RATNUMP(x)? SCM_RATNUM_NUMER(x) : x;
790 ScmObj dx = SCM_RATNUMP(x)? SCM_RATNUM_DENOM(x) : SCM_MAKE_INT(1);
791 ScmObj ny = SCM_RATNUMP(y)? SCM_RATNUM_NUMER(y) : y;
792 ScmObj dy = SCM_RATNUMP(y)? SCM_RATNUM_DENOM(y) : SCM_MAKE_INT(1);
793
794 if (divide) {
795 ScmObj t = ny; ny = dy; dy = t;
796 }
797 return Scm_MakeRational(Scm_Mul(nx, ny),
798 Scm_Mul(dx, dy));
799 }
800
801 #define Scm_RatnumAdd(x, y) Scm_RatnumAddSub(x, y, FALSE)
802 #define Scm_RatnumSub(x, y) Scm_RatnumAddSub(x, y, TRUE)
803 #define Scm_RatnumMul(x, y) Scm_RatnumMulDiv(x, y, FALSE)
804 #define Scm_RatnumDiv(x, y) Scm_RatnumMulDiv(x, y, TRUE)
805
806
807 /*=======================================================================
808 * Compnums
809 */
810
Scm_MakeCompnum(double r,double i)811 ScmObj Scm_MakeCompnum(double r, double i)
812 {
813 ScmCompnum *c = SCM_NEW_ATOMIC(ScmCompnum);
814 SCM_SET_CLASS(c, SCM_CLASS_COMPLEX);
815 c->real = r;
816 c->imag = i;
817 return SCM_OBJ(c);
818 }
819
Scm_MakeComplex(double r,double i)820 ScmObj Scm_MakeComplex(double r, double i)
821 {
822 if (i == 0.0) return Scm_MakeFlonum(r);
823 else return Scm_MakeCompnum(r, i);
824 }
825
Scm_MakeComplexPolar(double mag,double angle)826 ScmObj Scm_MakeComplexPolar(double mag, double angle)
827 {
828 double real = mag * cos(angle);
829 double imag = mag * sin(angle);
830 if (imag == 0.0) return Scm_MakeFlonum(real);
831 else return Scm_MakeCompnum(real, imag);
832 }
833
Scm_MakeComplexPolarPI(double mag,double pi_angle)834 ScmObj Scm_MakeComplexPolarPI(double mag, double pi_angle)
835 {
836 double real = mag * Scm_CosPi(pi_angle);
837 double imag = mag * Scm_SinPi(pi_angle);
838 if (imag == 0.0) return Scm_MakeFlonum(real);
839 else return Scm_MakeCompnum(real, imag);
840 }
841
842 /* NB: This isn't called by Scheme's real-part; see libnum.scm */
Scm_RealPart(ScmObj z)843 double Scm_RealPart(ScmObj z)
844 {
845 if (SCM_REALP(z)) {
846 return Scm_GetDouble(z);
847 }
848 if (!SCM_COMPNUMP(z)) {
849 Scm_Error("number required, but got %S", z);
850 return 0.0; /* dummy */
851 }
852 return SCM_COMPNUM_REAL(z);
853 }
854
855 /* NB: This isn't called by Scheme's imag-part; see libnum.scm */
Scm_ImagPart(ScmObj z)856 double Scm_ImagPart(ScmObj z)
857 {
858 if (SCM_COMPNUMP(z)) {
859 return SCM_COMPNUM_IMAG(z);
860 }
861 if (!SCM_REALP(z)) {
862 Scm_Error("number required, but got %S", z);
863 }
864 return 0.0;
865 }
866
Scm_GetHalfComplex(ScmObj z)867 ScmHalfComplex Scm_GetHalfComplex(ScmObj z)
868 {
869 ScmHalfComplex c = {0, 0};
870 if (SCM_COMPNUMP(z)) {
871 c.r = Scm_DoubleToHalf(SCM_COMPNUM_REAL(z));
872 c.i = Scm_DoubleToHalf(SCM_COMPNUM_IMAG(z));
873 } else if (SCM_REALP(z)) {
874 c.r = Scm_DoubleToHalf(Scm_GetDouble(z));
875 c.i = 0;
876 } else {
877 Scm_Error("number required, but got %S", z);
878 }
879 return c;
880 }
881
Scm_GetFloatComplex(ScmObj z)882 ScmFloatComplex Scm_GetFloatComplex(ScmObj z)
883 {
884 ScmFloatComplex c = 0.0f;
885 if (SCM_COMPNUMP(z)) {
886 c = (float)SCM_COMPNUM_REAL(z)
887 + (float)SCM_COMPNUM_IMAG(z) * _Complex_I;
888 } else if (SCM_REALP(z)) {
889 c = (float)Scm_GetDouble(z);
890 } else {
891 Scm_Error("number required, but got %S", z);
892 }
893 return c;
894 }
895
Scm_GetDoubleComplex(ScmObj z)896 ScmDoubleComplex Scm_GetDoubleComplex(ScmObj z)
897 {
898 ScmDoubleComplex c = 0.0f;
899 if (SCM_COMPNUMP(z)) {
900 c = SCM_COMPNUM_REAL(z)
901 + SCM_COMPNUM_IMAG(z) * _Complex_I;
902 } else if (SCM_REALP(z)) {
903 c = Scm_GetDouble(z);
904 } else {
905 Scm_Error("number required, but got %S", z);
906 }
907 return c;
908 }
909
Scm_HalfComplexToComplex(ScmHalfComplex z)910 ScmObj Scm_HalfComplexToComplex(ScmHalfComplex z)
911 {
912 return Scm_MakeComplex(Scm_HalfToDouble(z.r),
913 Scm_HalfToDouble(z.i));
914 }
915
Scm_FloatComplexToComplex(ScmFloatComplex z)916 ScmObj Scm_FloatComplexToComplex(ScmFloatComplex z)
917 {
918 return Scm_MakeComplex((double)crealf(z), (double)cimagf(z));
919 }
920
Scm_DoubleComplexToComplex(ScmDoubleComplex z)921 ScmObj Scm_DoubleComplexToComplex(ScmDoubleComplex z)
922 {
923 return Scm_MakeComplex(creal(z), cimag(z));
924 }
925
926
927 /* NB: This isn't called by Scheme's magnitude; see libnum.scm */
Scm_Magnitude(ScmObj z)928 double Scm_Magnitude(ScmObj z)
929 {
930 if (SCM_REALP(z)) {
931 return fabs(Scm_GetDouble(z));
932 }
933 if (!SCM_COMPNUMP(z)) {
934 Scm_Error("number required, but got %S", z);
935 return 0.0; /* dummy */
936 }
937 double r = SCM_COMPNUM_REAL(z);
938 double i = SCM_COMPNUM_IMAG(z);
939 return sqrt(r*r+i*i);
940 }
941
Scm_Angle(ScmObj z)942 double Scm_Angle(ScmObj z)
943 {
944 if (SCM_REALP(z)) {
945 return (Scm_Sign(z) < 0)? M_PI : 0.0;
946 }
947 if (!SCM_COMPNUMP(z)) {
948 Scm_Error("number required, but got %S", z);
949 return 0.0; /* dummy */
950 }
951 double r = SCM_COMPNUM_REAL(z);
952 double i = SCM_COMPNUM_IMAG(z);
953 return atan2(i, r);
954 }
955
956 /*=======================================================================
957 * Coertion
958 */
959
Scm_MakeInteger(long i)960 ScmObj Scm_MakeInteger(long i)
961 {
962 if (i >= SCM_SMALL_INT_MIN && i <= SCM_SMALL_INT_MAX) {
963 return SCM_MAKE_INT(i);
964 } else {
965 return Scm_MakeBignumFromSI(i);
966 }
967 }
968
Scm_MakeIntegerU(u_long i)969 ScmObj Scm_MakeIntegerU(u_long i)
970 {
971 if (i <= (u_long)SCM_SMALL_INT_MAX) return SCM_MAKE_INT(i);
972 else return Scm_MakeBignumFromUI(i);
973 }
974
range_error(ScmObj obj,int clamp,int * oor)975 static void range_error(ScmObj obj, int clamp, int *oor)
976 {
977 if (clamp != SCM_CLAMP_ERROR && oor != NULL) {
978 *oor = TRUE;
979 } else {
980 Scm_Error("argument out of range: %S", obj);
981 }
982 }
983
984 /* Convert scheme integer to C integer */
Scm_GetIntegerClamp(ScmObj obj,int clamp,int * oor)985 long Scm_GetIntegerClamp(ScmObj obj, int clamp, int *oor)
986 {
987 double v = 0.0;
988 if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
989 if (SCM_INTP(obj)) return SCM_INT_VALUE(obj);
990 else if (SCM_BIGNUMP(obj)) {
991 return Scm_BignumToSI(SCM_BIGNUM(obj), clamp, oor);
992 }
993 else if (SCM_FLONUMP(obj)) {
994 v = SCM_FLONUM_VALUE(obj);
995 goto flonum;
996 }
997 else if (SCM_RATNUMP(obj)) {
998 v = Scm_GetDouble(obj);
999 goto flonum;
1000 }
1001 else {
1002 goto err;
1003 }
1004 flonum:
1005 if (v > (double)LONG_MAX) {
1006 if (clamp & SCM_CLAMP_HI) return LONG_MAX;
1007 else goto err;
1008 }
1009 if (v < (double)LONG_MIN) {
1010 if (clamp & SCM_CLAMP_LO) return LONG_MIN;
1011 else goto err;
1012 }
1013 return (long)v;
1014 err:
1015 range_error(obj, clamp, oor);
1016 return 0;
1017 }
1018
Scm_GetIntegerUClamp(ScmObj obj,int clamp,int * oor)1019 u_long Scm_GetIntegerUClamp(ScmObj obj, int clamp, int *oor)
1020 {
1021 double v = 0.0;
1022
1023 if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
1024 if (SCM_INTP(obj)) {
1025 if (SCM_INT_VALUE(obj) < 0) {
1026 if (clamp & SCM_CLAMP_LO) return 0;
1027 else goto err;
1028 }
1029 return SCM_INT_VALUE(obj);
1030 }
1031 else if (SCM_BIGNUMP(obj)) {
1032 return Scm_BignumToUI(SCM_BIGNUM(obj), clamp, oor);
1033 }
1034 else if (SCM_FLONUMP(obj)) {
1035 v = SCM_FLONUM_VALUE(obj);
1036 goto flonum;
1037 }
1038 else if (SCM_RATNUMP(obj)) {
1039 v = Scm_GetDouble(obj);
1040 goto flonum;
1041 }
1042 else {
1043 goto err;
1044 }
1045 flonum:
1046 if (v > (double)ULONG_MAX) {
1047 if (clamp & SCM_CLAMP_HI) return ULONG_MAX;
1048 else goto err;
1049 }
1050 if (v < 0.0) {
1051 if (clamp & SCM_CLAMP_LO) return 0;
1052 else goto err;
1053 }
1054 return (u_long)v;
1055 err:
1056 range_error(obj, clamp, oor);
1057 return 0;
1058 }
1059
1060 /* 8- and 16-bit integer extraction with range check */
1061 #define SMALL_INT_XTRACT(name, upper, lower) \
1062 name(ScmObj obj, int clamp, int *oor) \
1063 { \
1064 long n = 0; \
1065 if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE; \
1066 if (SCM_INTP(obj)) { \
1067 n = SCM_INT_VALUE(obj); \
1068 } else if (SCM_FLONUMP(obj)) { \
1069 n = (long)SCM_FLONUM_VALUE(obj); \
1070 } else if (SCM_RATNUMP(obj)) { \
1071 n = (long)Scm_GetDouble(obj); \
1072 } else if (SCM_BIGNUMP(obj)) { \
1073 if (Scm_Sign(obj) > 0) { \
1074 if (clamp & SCM_CLAMP_HI) return upper; \
1075 else goto err; \
1076 } else { \
1077 if (clamp & SCM_CLAMP_LO) return lower; \
1078 else goto err; \
1079 } \
1080 } else { \
1081 goto err; \
1082 } \
1083 if (n > upper) { \
1084 if (clamp & SCM_CLAMP_HI) return upper; \
1085 else goto err; \
1086 } \
1087 if (n < lower) { \
1088 if (clamp & SCM_CLAMP_LO) return lower; \
1089 else goto err; \
1090 } \
1091 return n; \
1092 err: \
1093 range_error(obj, clamp, oor); \
1094 return 0; \
1095 }
1096
1097 SMALL_INT_XTRACT(int Scm_GetInteger8Clamp, 127, -128)
1098 SMALL_INT_XTRACT(u_int Scm_GetIntegerU8Clamp, 255, 0)
1099 SMALL_INT_XTRACT(int Scm_GetInteger16Clamp, 32767, -32768)
1100 SMALL_INT_XTRACT(u_int Scm_GetIntegerU16Clamp, 65535, 0)
1101
1102
1103 /* 32bit integer specific */
Scm_GetInteger32Clamp(ScmObj obj,int clamp,int * oor)1104 int32_t Scm_GetInteger32Clamp(ScmObj obj, int clamp, int *oor)
1105 {
1106 #if SIZEOF_LONG == 4
1107 return (int32_t)Scm_GetIntegerClamp(obj, clamp, oor);
1108 #else /* SIZEOF_LONG >= 8 */
1109
1110 if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
1111 /* NB: we denote the constant directly here. (1L<<31) fails on
1112 Alpha machines, since the compiler somehow calculates the constant
1113 in 32bit integer even it has 'L'. We have to write (1LL<<31), but
1114 I'm afraid that it's not portable. */
1115 if (SCM_INTP(obj)) {
1116 long r = SCM_INT_VALUE(obj);
1117 if (r < -0x80000000L) {
1118 if (clamp & SCM_CLAMP_LO) return -0x80000000L;
1119 goto err;
1120 }
1121 if (r > 0x7fffffffL) {
1122 if (clamp & SCM_CLAMP_HI) return 0x7fffffffL;
1123 goto err;
1124 }
1125 return r;
1126 } else if (SCM_BIGNUMP(obj)) {
1127 if (SCM_BIGNUM_SIGN(obj) < 0) {
1128 if (clamp & SCM_CLAMP_LO) return -0x80000000L;
1129 goto err;
1130 } else {
1131 if (clamp & SCM_CLAMP_HI) return 0x7fffffffL;
1132 goto err;
1133 }
1134 }
1135 /*TODO: flonum and ratnum! */
1136 err:
1137 range_error(obj, clamp, oor);
1138 return 0;
1139 #endif /* SIZEOF_LONG >= 8 */
1140 }
1141
Scm_GetIntegerU32Clamp(ScmObj obj,int clamp,int * oor)1142 uint32_t Scm_GetIntegerU32Clamp(ScmObj obj, int clamp, int *oor)
1143 {
1144 #if SIZEOF_LONG == 4
1145 return (uint32_t)Scm_GetIntegerUClamp(obj, clamp, oor);
1146 #else /* SIZEOF_LONG >= 8 */
1147 if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
1148 if (SCM_INTP(obj)) {
1149 long r = SCM_INT_VALUE(obj);
1150 if (r < 0) {
1151 if (clamp & SCM_CLAMP_LO) return 0;
1152 goto err;
1153 }
1154 if (r > 0xffffffffL) {
1155 if (clamp & SCM_CLAMP_HI) return 0xffffffffUL;
1156 goto err;
1157 }
1158 return r;
1159 } else if (SCM_BIGNUMP(obj)) {
1160 if (SCM_BIGNUM_SIGN(obj) < 0) {
1161 if (clamp & SCM_CLAMP_LO) return 0;
1162 goto err;
1163 } else {
1164 if (clamp & SCM_CLAMP_HI) return 0xffffffffUL;
1165 goto err;
1166 }
1167 }
1168 err:
1169 range_error(obj, clamp, oor);
1170 return 0;
1171 #endif /* SIZEOF_LONG >= 8 */
1172 }
1173
1174 /* get an unsigned integer value modulo u_long range.
1175 convenient when you only concern lower bits. */
Scm_GetIntegerUMod(ScmObj obj)1176 u_long Scm_GetIntegerUMod(ScmObj obj)
1177 {
1178 if (SCM_INTP(obj)) return (u_long)SCM_INT_VALUE(obj);
1179 if (SCM_BIGNUMP(obj)) {
1180 if (SCM_BIGNUM_SIZE(obj) == 0) {
1181 /* this shouldn't happen in normalized bignums, but just in case */
1182 return 0;
1183 }
1184 if (SCM_BIGNUM_SIGN(obj) < 0) {
1185 u_long v = SCM_BIGNUM(obj)->values[0];
1186 return ~v + 1;
1187 } else {
1188 return SCM_BIGNUM(obj)->values[0];
1189 }
1190 }
1191 Scm_Error("Exact integer required, but got %S", obj);
1192 return 0; /* dummy */
1193 }
1194
1195 #if SIZEOF_LONG == 4
1196 /* we need special routines */
Scm_MakeInteger64(int64_t i)1197 ScmObj Scm_MakeInteger64(int64_t i)
1198 {
1199 u_long val[2];
1200 val[0] = (u_long)((uint64_t)i & ULONG_MAX);
1201 val[1] = (u_long)((uint64_t)i >> 32);
1202 if (val[1] == 0 && val[0] <= LONG_MAX) return Scm_MakeInteger(val[0]);
1203 return Scm_NormalizeBignum(SCM_BIGNUM(Scm_MakeBignumFromUIArray(0, val, 2)));
1204 }
1205
Scm_MakeIntegerU64(uint64_t i)1206 ScmObj Scm_MakeIntegerU64(uint64_t i)
1207 {
1208 u_long val[2];
1209 val[0] = (u_long)((uint64_t)i & ULONG_MAX);
1210 val[1] = (u_long)((uint64_t)i >> 32);
1211 if (val[1] == 0) return Scm_MakeIntegerU(val[0]);
1212 return Scm_MakeBignumFromUIArray(1, val, 2);
1213 }
1214
Scm_GetInteger64Clamp(ScmObj obj,int clamp,int * oor)1215 int64_t Scm_GetInteger64Clamp(ScmObj obj, int clamp, int *oor)
1216 {
1217 int64_t r = 0;
1218 if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
1219 if (SCM_INTP(obj)) return (int64_t)SCM_INT_VALUE(obj);
1220 if (SCM_BIGNUMP(obj)) {
1221 return Scm_BignumToSI64(SCM_BIGNUM(obj), clamp, oor);
1222 }
1223 if (SCM_RATNUMP(obj)) {
1224 obj = Scm_Inexact(obj);
1225 /* FALLTHROUGH */
1226 }
1227 if (SCM_FLONUMP(obj)) {
1228 int64_t maxval, minval;
1229 double v;
1230
1231 SCM_SET_INT64_MAX(maxval);
1232 SCM_SET_INT64_MIN(minval);
1233 v = SCM_FLONUM_VALUE(obj);
1234 if (v > (double)maxval) {
1235 if (!(clamp&SCM_CLAMP_HI)) goto err;
1236 return maxval;
1237 } else if (v < (double)minval) {
1238 if (!(clamp&SCM_CLAMP_LO)) goto err;
1239 return minval;
1240 } else {
1241 return (int64_t)v;
1242 }
1243 }
1244 err:
1245 range_error(obj, clamp, oor);
1246 return r;
1247 }
1248
Scm_GetIntegerU64Clamp(ScmObj obj,int clamp,int * oor)1249 uint64_t Scm_GetIntegerU64Clamp(ScmObj obj, int clamp, int *oor)
1250 {
1251 uint64_t r = 0;
1252 if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
1253 if (SCM_INTP(obj)) {
1254 long v = SCM_INT_VALUE(obj);
1255 if (v < 0) {
1256 if (!(clamp&SCM_CLAMP_LO)) goto err;
1257 return 0;
1258 } else {
1259 return (uint64_t)v;
1260 }
1261 }
1262 if (SCM_BIGNUMP(obj)) {
1263 return Scm_BignumToUI64(SCM_BIGNUM(obj), clamp, oor);
1264 }
1265 if (SCM_RATNUMP(obj)) {
1266 obj = Scm_Inexact(obj);
1267 /* FALLTHROUGH */
1268 }
1269 if (SCM_FLONUMP(obj)) {
1270 double v = SCM_FLONUM_VALUE(obj);
1271 uint64_t maxval;
1272
1273 if (v < 0) {
1274 if (!(clamp&SCM_CLAMP_LO)) goto err;
1275 return 0;
1276 }
1277 SCM_SET_UINT64_MAX(maxval);
1278 if (v > (double)maxval) {
1279 if (!(clamp&SCM_CLAMP_HI)) goto err;
1280 return maxval;
1281 } else {
1282 return (uint64_t)v;
1283 }
1284 }
1285 err:
1286 range_error(obj, clamp, oor);
1287 return r;
1288 }
1289
1290 #endif /* SIZEOF_LONG == 4 */
1291
1292 #if GAUCHE_API_VERSION < 1000
1293 /* Since we adopt C99, we can simply use cast instead of the funcitons.
1294 These are kept for ABI compatibility until 1.0 release */
Scm_DoubleToInt64(double v)1295 int64_t Scm_DoubleToInt64(double v)
1296 {
1297 return (int64_t)v;
1298 }
1299
Scm_DoubleToUInt64(double v)1300 uint64_t Scm_DoubleToUInt64(double v)
1301 {
1302 return (uint64_t)v;
1303 }
1304
Scm_Int64ToDouble(int64_t v)1305 double Scm_Int64ToDouble(int64_t v)
1306 {
1307 return (double)v;
1308 }
1309
Scm_UInt64ToDouble(uint64_t v)1310 double Scm_UInt64ToDouble(uint64_t v)
1311 {
1312 return (double)v;
1313 }
1314 #endif /*GAUCHE_API_VERSION < 1000*/
1315
1316 /* See if a Scheme integer si can be representable with 53bit mantissa. */
double_precision(ScmObj si,int * hi,int * lo)1317 static int double_precision(ScmObj si, int *hi, int *lo)
1318 {
1319 if (SCM_INTP(si)) {
1320 ScmSmallInt i = SCM_INT_VALUE(si);
1321 if (i < 0) i = -i; /* won't overflow */
1322 ScmBits bi = (ScmBits)i;
1323 *hi = Scm_BitsHighest1(&bi, 0, SCM_WORD_BITS-1);
1324 *lo = Scm_BitsLowest1(&bi, 0, SCM_WORD_BITS-1);
1325 return (*hi - *lo) < 53;
1326 } else {
1327 SCM_ASSERT(SCM_BIGNUMP(si));
1328 const ScmBits *bits = (ScmBits*)SCM_BIGNUM(si)->values;
1329 *hi = Scm_BitsHighest1(bits, 0, SCM_BIGNUM_SIZE(si)*SCM_WORD_BITS-1);
1330 *lo = Scm_BitsLowest1(bits, 0, SCM_BIGNUM_SIZE(si)*SCM_WORD_BITS-1);
1331 return (*hi - *lo) < 53;
1332 }
1333 }
1334
1335 /* Get n-th bit of absolute value of Scheme integer si */
abs_bittest(ScmObj si,int bit)1336 static int abs_bittest(ScmObj si, int bit)
1337 {
1338 if (SCM_INTP(si)) {
1339 ScmSmallInt i = SCM_INT_VALUE(si);
1340 if (i < 0) i = -i;
1341 return i & (1L<<bit);
1342 } else {
1343 SCM_ASSERT(SCM_BIGNUMP(si));
1344 ScmBits *bits = (ScmBits*)SCM_BIGNUM(si)->values;
1345 return SCM_BITS_TEST(bits, bit);
1346 }
1347 }
1348
Scm_GetDouble(ScmObj obj)1349 double Scm_GetDouble(ScmObj obj)
1350 {
1351 if (SCM_FLONUMP(obj)) return SCM_FLONUM_VALUE(obj);
1352 else if (SCM_INTP(obj)) return (double)SCM_INT_VALUE(obj);
1353 else if (SCM_BIGNUMP(obj)) return Scm_BignumToDouble(SCM_BIGNUM(obj));
1354 else if (SCM_RATNUMP(obj)) {
1355 /* This is more subtle than it appears. A naive approach is to
1356 convert numerator and denominator to double, and do flonum division.
1357 However,
1358 - Denominator and/or numerator may exceed FLT_MAX, yielding
1359 infinities in the intermediate results. However, their ratio
1360 can be in a valid flonum range.
1361 - If denominator and/or numerator requires more than 53bit precision,
1362 converting each to double causes rounding, which occurs before
1363 the final division, causing ULP-off error.
1364 E.g. (inexact (/ (+ 1 (* (exact (flonum-epsilon)) 33/100)) 1))
1365 should be 1, but would yield 1.0000000000000002 if the numerator
1366 is rounded up first.
1367 */
1368 ScmObj numer = SCM_RATNUM_NUMER(obj);
1369 ScmObj denom = SCM_RATNUM_DENOM(obj);
1370
1371 int n_hi, n_lo, d_hi, d_lo;
1372 int n_dp = double_precision(numer, &n_hi, &n_lo);
1373 int d_dp = double_precision(denom, &d_hi, &d_lo);
1374
1375 if (!(n_dp && d_dp)) goto fullpath;
1376
1377 double dnumer = Scm_GetDouble(numer);
1378 double ddenom = Scm_GetDouble(denom);
1379
1380 if (!SCM_IS_INF(dnumer) && !SCM_IS_INF(ddenom)) {
1381 /* short path */
1382 volatile double result;
1383 SCM_FP_ENSURE_DOUBLE_PRECISION_BEGIN();
1384 /* It is critical to perform this division in IEEE double (53bit
1385 mantissa) precision, _not_ in x87 extended double precision;
1386 if the latter were used, double-rounding would yield different
1387 results, which makes inexact -> exact -> inexact round-trip
1388 fail. For example, (inexact 3002399751580332/3002399751580331)
1389 should be 1.0000000000000002 (1LSB greater than 1.0), but
1390 extended double precision division yields
1391 1.0000000000000004 (2LSB greater than 1.0).
1392 */
1393 result = dnumer/ddenom;
1394 SCM_FP_ENSURE_DOUBLE_PRECISION_END();
1395 return result;
1396 }
1397 fullpath:;
1398 /* Need more precise but expensive calculation.
1399 We find K such that 2^K * numer >= 2^54 * denom, so that
1400 the integer division yiels more than 53bit integral part.
1401 */
1402 int shift = 0;
1403 if (n_hi - d_hi < 54) {
1404 shift = 54 - (n_hi - d_hi);
1405 numer = Scm_Ash(numer, shift);
1406 }
1407
1408 ScmObj rem;
1409 ScmObj quo = Scm_Quotient(numer, denom, &rem);
1410
1411 /* If shift > 1076, the result would be denomarlized range.
1412 Less than 53bits of quo is used, so we have to mask out the
1413 unncessary digit. (Otherwise, it causes double-rounding.) */
1414 if (shift > 1076) {
1415 ScmObj mask =
1416 Scm_LogNot(Scm_Sub(Scm_Ash(SCM_MAKE_INT(1), shift-1076-1),
1417 SCM_MAKE_INT(1)));
1418 if (Scm_Sign(quo) < 0) {
1419 quo = Scm_Negate(Scm_LogAnd(Scm_Negate(quo), mask));
1420 } else {
1421 quo = Scm_LogAnd(quo, mask);
1422 }
1423 }
1424
1425 int q_hi, q_lo;
1426 if (double_precision(quo, &q_hi, &q_lo)) {
1427 /* Result fits in double precision. */
1428 double dquo = Scm_GetDouble(quo);
1429 return ldexp(dquo, -shift);
1430 }
1431 /* We have to look at the 54th bit and below to decide rounding.
1432 If 54-th bit is 0, we can safely convert it to double, truncating
1433 the lower bits. */
1434 if (!abs_bittest(quo, q_hi - 53)) {
1435 double dquo = Scm_GetDouble(quo);
1436 return ldexp(dquo, -shift);
1437 }
1438 /* If 54-th bit is 1, we see the lower bits. If we have at least
1439 one '1' bit, or remainder isn't zero, we can round up. */
1440 int roundup = FALSE;
1441 ScmObj mask = Scm_Ash(SCM_MAKE_INT(1), q_hi - 53);
1442 if (rem != SCM_MAKE_INT(0)) roundup = TRUE;
1443 else {
1444 ScmObj mask_1 = Scm_Sub(mask, SCM_MAKE_INT(1));
1445 ScmObj q = quo;
1446 if (Scm_Sign(q) < 0) {
1447 q = Scm_Negate(q);
1448 }
1449 if (Scm_LogAnd(q, mask_1) != SCM_MAKE_INT(0)) roundup = TRUE;
1450 }
1451 if (roundup) {
1452 if (Scm_Sign(quo) < 0) {
1453 quo = Scm_Sub(quo, mask);
1454 } else {
1455 quo = Scm_Add(quo, mask);
1456 }
1457 double dquo = Scm_GetDouble(quo);
1458 return ldexp(dquo, -shift);
1459 }
1460 /* We are half-point. See the 53-bit and round to even. */
1461 if (!abs_bittest(quo, q_hi - 52)) {
1462 if (Scm_Sign(quo) < 0) {
1463 quo = Scm_Add(quo, mask);
1464 } else {
1465 quo = Scm_Sub(quo, mask);
1466 }
1467 }
1468 double dquo = Scm_GetDouble(quo);
1469 return ldexp(dquo, -shift);
1470 }
1471 else return 0.0;
1472 }
1473
1474 /*
1475 * Generic Methods
1476 */
1477
1478 /* Predicates */
1479
Scm_IntegerP(ScmObj obj)1480 int Scm_IntegerP(ScmObj obj)
1481 {
1482 if (SCM_INTP(obj) || SCM_BIGNUMP(obj)) return TRUE;
1483 if (SCM_RATNUMP(obj)) return FALSE; /* normalized ratnum never be integer */
1484 if (SCM_FLONUMP(obj)) {
1485 double d = SCM_FLONUM_VALUE(obj);
1486 if (SCM_IS_INF(d) || SCM_IS_NAN(d)) return FALSE;
1487 double i;
1488 double f = modf(d, &i);
1489 if (f == 0.0) return TRUE;
1490 return FALSE;
1491 }
1492 if (SCM_COMPNUMP(obj)) return FALSE;
1493 Scm_Error("number required, but got %S", obj);
1494 return FALSE; /* dummy */
1495 }
1496
Scm_OddP(ScmObj obj)1497 int Scm_OddP(ScmObj obj)
1498 {
1499 if (SCM_INTP(obj)) {
1500 return (SCM_INT_VALUE(obj)&1);
1501 }
1502 if (SCM_BIGNUMP(obj)) {
1503 return (SCM_BIGNUM(obj)->values[0] & 1);
1504 }
1505 if (SCM_FLONUMP(obj) && Scm_IntegerP(obj)) {
1506 return (fmod(SCM_FLONUM_VALUE(obj), 2.0) != 0.0);
1507 }
1508 Scm_Error("integer required, but got %S", obj);
1509 return FALSE; /* dummy */
1510
1511 }
1512
Scm_FiniteP(ScmObj obj)1513 int Scm_FiniteP(ScmObj obj)
1514 {
1515 return !Scm_InfiniteP(obj) && !Scm_NanP(obj);
1516 }
1517
Scm_InfiniteP(ScmObj obj)1518 int Scm_InfiniteP(ScmObj obj)
1519 {
1520 if (SCM_FLONUMP(obj)) {
1521 double v = SCM_FLONUM_VALUE(obj);
1522 return SCM_IS_INF(v);
1523 } else if (SCM_COMPNUMP(obj)) {
1524 double r = SCM_COMPNUM_REAL(obj);
1525 double i = SCM_COMPNUM_IMAG(obj);
1526 return SCM_IS_INF(r) || SCM_IS_INF(i);
1527 } else if (!SCM_NUMBERP(obj)) {
1528 SCM_TYPE_ERROR(obj, "number");
1529 }
1530 return FALSE;
1531 }
1532
Scm_NanP(ScmObj obj)1533 int Scm_NanP(ScmObj obj)
1534 {
1535 if (SCM_FLONUMP(obj)) {
1536 double v = SCM_FLONUM_VALUE(obj);
1537 return SCM_IS_NAN(v);
1538 } else if (SCM_COMPNUMP(obj)) {
1539 double r = SCM_COMPNUM_REAL(obj);
1540 double i = SCM_COMPNUM_IMAG(obj);
1541 return SCM_IS_NAN(r) || SCM_IS_NAN(i);
1542 } else if (!SCM_NUMBERP(obj)) {
1543 SCM_TYPE_ERROR(obj, "number");
1544 }
1545 return FALSE;
1546 }
1547
1548 /* Unary Operator */
1549
scm_abs(ScmObj obj,int vmp)1550 static ScmObj scm_abs(ScmObj obj, int vmp)
1551 {
1552 if (SCM_INTP(obj)) {
1553 long v = SCM_INT_VALUE(obj);
1554 if (v < 0) {
1555 if (v == SCM_SMALL_INT_MIN) {
1556 obj = Scm_MakeBignumFromSI(-v);
1557 } else {
1558 obj = SCM_MAKE_INT(-v);
1559 }
1560 }
1561 } else if (SCM_BIGNUMP(obj)) {
1562 if (SCM_BIGNUM_SIGN(obj) < 0) {
1563 obj = Scm_BignumCopy(SCM_BIGNUM(obj));
1564 SCM_BIGNUM_SIGN(obj) = 1;
1565 }
1566 } else if (SCM_FLONUMP(obj)) {
1567 double v = SCM_FLONUM_VALUE(obj);
1568 if (v < 0) RETURN_FLONUM(-v);
1569 } else if (SCM_RATNUMP(obj)) {
1570 if (Scm_Sign(SCM_RATNUM_NUMER(obj)) < 0) {
1571 obj = Scm_MakeRational(Scm_Negate(SCM_RATNUM_NUMER(obj)),
1572 SCM_RATNUM_DENOM(obj));
1573 }
1574 } else if (SCM_COMPNUMP(obj)) {
1575 double r = SCM_COMPNUM_REAL(obj);
1576 double i = SCM_COMPNUM_IMAG(obj);
1577 double a = sqrt(r*r+i*i);
1578 RETURN_FLONUM(a);
1579 } else {
1580 Scm_Error("number required: %S", obj);
1581 }
1582 return obj;
1583 }
DEFINE_DUAL_API1(Scm_Abs,Scm_VMAbs,scm_abs)1584 DEFINE_DUAL_API1(Scm_Abs, Scm_VMAbs, scm_abs)
1585
1586 /* Return -1, 0 or 1 when arg is minus, zero or plus, respectively.
1587 used to implement zero?, positive? and negative?
1588 NB: This returns 0 for both positive and negative zeros. */
1589 int Scm_Sign(ScmObj obj)
1590 {
1591 if (SCM_INTP(obj)) {
1592 long r = SCM_INT_VALUE(obj);
1593 if (r == 0) return 0;
1594 return (r > 0)? 1 : -1;
1595 }
1596 if (SCM_BIGNUMP(obj)) {
1597 return SCM_BIGNUM_SIGN(obj);
1598 }
1599 if (SCM_FLONUMP(obj)) {
1600 double v = SCM_FLONUM_VALUE(obj);
1601 if (v == 0.0) return 0;
1602 return (v > 0.0)? 1 : -1;
1603 }
1604 if (SCM_RATNUMP(obj)) {
1605 return Scm_Sign(SCM_RATNUM_NUMER(obj));
1606 }
1607 /* NB: zero? can accept a complex number, but it is processed in
1608 the stub function. see libnum.scm */
1609 Scm_Error("real number required, but got %S", obj);
1610 return 0; /* dummy */
1611 }
1612
negate(ScmObj obj,int vmp)1613 static ScmObj negate(ScmObj obj, int vmp)
1614 {
1615 if (SCM_INTP(obj)) {
1616 long v = SCM_INT_VALUE(obj);
1617 if (v == SCM_SMALL_INT_MIN) {
1618 return Scm_MakeBignumFromSI(-v);
1619 } else {
1620 return SCM_MAKE_INT(-v);
1621 }
1622 } else if (SCM_BIGNUMP(obj)) {
1623 return Scm_BignumNegate(SCM_BIGNUM(obj));
1624 } else if (SCM_FLONUMP(obj)) {
1625 double r = -SCM_FLONUM_VALUE(obj);
1626 RETURN_FLONUM(r);
1627 } else if (SCM_RATNUMP(obj)) {
1628 return Scm_MakeRational(Scm_Negate(SCM_RATNUM_NUMER(obj)),
1629 SCM_RATNUM_DENOM(obj));
1630 } else if (SCM_COMPNUMP(obj)) {
1631 return Scm_MakeCompnum(-SCM_COMPNUM_REAL(obj),
1632 -SCM_COMPNUM_IMAG(obj));
1633 } else {
1634 return Scm_ApplyRec(SCM_OBJ(&generic_sub), SCM_LIST1(obj));
1635 }
1636 }
DEFINE_DUAL_API1(Scm_Negate,Scm_VMNegate,negate)1637 DEFINE_DUAL_API1(Scm_Negate, Scm_VMNegate, negate)
1638
1639
1640 static ScmObj reciprocal(ScmObj obj, int vmp)
1641 {
1642 if (SCM_INTP(obj) || SCM_BIGNUMP(obj)) {
1643 return Scm_MakeRational(SCM_MAKE_INT(1), obj);
1644 } else if (SCM_FLONUMP(obj)) {
1645 double val = 1.0/SCM_FLONUM_VALUE(obj);
1646 RETURN_FLONUM(val);
1647 } else if (SCM_RATNUMP(obj)) {
1648 return Scm_MakeRational(SCM_RATNUM_DENOM(obj),
1649 SCM_RATNUM_NUMER(obj));
1650 } else if (SCM_COMPNUMP(obj)) {
1651 double r = SCM_COMPNUM_REAL(obj);
1652 double i = SCM_COMPNUM_IMAG(obj);
1653 double d = r*r + i*i;
1654 double r1 = r/d;
1655 double i1 = -i/d;
1656 return Scm_MakeComplex(r1, i1);
1657 } else {
1658 return Scm_ApplyRec(SCM_OBJ(&generic_div), SCM_LIST1(obj));
1659 }
1660 }
DEFINE_DUAL_API1(Scm_Reciprocal,Scm_VMReciprocal,reciprocal)1661 DEFINE_DUAL_API1(Scm_Reciprocal, Scm_VMReciprocal, reciprocal)
1662
1663
1664 static ScmObj ireciprocal(ScmObj obj, int vmp)
1665 {
1666 if (SCM_EXACT_ZERO_P(obj)) return SCM_POSITIVE_INFINITY;
1667 if (SCM_EXACT_ONE_P(obj)) return obj;
1668 if (SCM_REALP(obj)) {
1669 double z = 1.0/Scm_GetDouble(obj);
1670 RETURN_FLONUM(z);
1671 }
1672 // delegate the rest to exact reciprocal
1673 return reciprocal(obj, vmp);
1674 }
DEFINE_DUAL_API1(Scm_ReciprocalInexact,Scm_VMReciprocalInexact,ireciprocal)1675 DEFINE_DUAL_API1(Scm_ReciprocalInexact, Scm_VMReciprocalInexact, ireciprocal)
1676
1677
1678 /*
1679 * Conversion operators
1680 */
1681
1682 static ScmObj inexact(ScmObj obj, int vmp)
1683 {
1684 if (SCM_INTP(obj)) {
1685 double z = (double)SCM_INT_VALUE(obj);
1686 RETURN_FLONUM(z);
1687 } else if (SCM_BIGNUMP(obj)) {
1688 double z = Scm_BignumToDouble(SCM_BIGNUM(obj));
1689 RETURN_FLONUM(z);
1690 } else if (SCM_RATNUMP(obj)) {
1691 RETURN_FLONUM(Scm_GetDouble(obj));
1692 } else if (!SCM_FLONUMP(obj) && !SCM_COMPNUMP(obj)) {
1693 Scm_Error("number required: %S", obj);
1694 }
1695 return obj;
1696 }
DEFINE_DUAL_API1(Scm_Inexact,Scm_VMInexact,inexact)1697 DEFINE_DUAL_API1(Scm_Inexact, Scm_VMInexact, inexact)
1698
1699
1700 ScmObj Scm_Exact(ScmObj obj)
1701 {
1702 if (SCM_FLONUMP(obj)) {
1703 double d = SCM_FLONUM_VALUE(obj);
1704 double f, i;
1705 if (SCM_IS_NAN(d) || SCM_IS_INF(d)) {
1706 Scm_Error("Exact infinity/nan is not supported: %S", obj);
1707 }
1708 if ((f = modf(d, &i)) == 0.0) {
1709 obj = Scm_FlonumIntegerToExact(i);
1710 } else {
1711 /* We'd find out the simplest rational numebr within the precision
1712 of IEEE double floating point number. The actual code is in
1713 lib/gauche/numerical.scm. */
1714 static ScmObj real_to_rational = SCM_UNDEFINED;
1715 SCM_BIND_PROC(real_to_rational, "real->rational",
1716 Scm_GaucheModule());
1717 obj = Scm_ApplyRec1(real_to_rational, obj);
1718 }
1719 } else if (SCM_COMPNUMP(obj)) {
1720 Scm_Error("exact complex is not supported: %S", obj);
1721 } if (!SCM_EXACTP(obj)) {
1722 Scm_Error("number required: %S", obj);
1723 }
1724 return obj;
1725 }
1726
1727 /*===============================================================
1728 * Arithmetics
1729 */
1730
1731 /* NB: we used to support n-ary operations in C API, expecting
1732 them to be faster since we can carry around the intermediate
1733 results unboxed. The newer versions of compiler, however,
1734 decomposes n-ary arithmetic operations into binary ones
1735 during optimization, so n-ary API hadn't really been used much.
1736 So we dropped them, in favor of simple code. */
1737
1738 /*
1739 * Addition and subtraction
1740 */
1741
scm_add(ScmObj arg0,ScmObj arg1,int vmp)1742 static ScmObj scm_add(ScmObj arg0, ScmObj arg1, int vmp)
1743 {
1744 if (SCM_INTP(arg0)) {
1745 if (SCM_INTP(arg1)) {
1746 long r = SCM_INT_VALUE(arg0) + SCM_INT_VALUE(arg1);
1747 return Scm_MakeInteger(r);
1748 }
1749 if (SCM_BIGNUMP(arg1)) {
1750 if (SCM_EXACT_ZERO_P(arg0)) return arg1;
1751 return Scm_BignumAddSI(SCM_BIGNUM(arg1), SCM_INT_VALUE(arg0));
1752 }
1753 if (SCM_RATNUMP(arg1)) {
1754 if (SCM_EXACT_ZERO_P(arg0)) return arg1;
1755 return Scm_RatnumAdd(arg0, arg1);
1756 }
1757 if (SCM_FLONUMP(arg1)) {
1758 if (SCM_EXACT_ZERO_P(arg0)) return arg1;
1759 double z = (double)SCM_INT_VALUE(arg0) + SCM_FLONUM_VALUE(arg1);
1760 RETURN_FLONUM(z);
1761 }
1762 if (SCM_COMPNUMP(arg1)) {
1763 if (SCM_EXACT_ZERO_P(arg0)) return arg1;
1764 return Scm_MakeComplex((double)SCM_INT_VALUE(arg0)
1765 + SCM_COMPNUM_REAL(arg1),
1766 SCM_COMPNUM_IMAG(arg1));
1767 }
1768 /* fallback to generic */
1769 }
1770 else if (SCM_BIGNUMP(arg0)) {
1771 if (SCM_INTP(arg1)) {
1772 if (SCM_EXACT_ZERO_P(arg1)) return arg0;
1773 return Scm_BignumAddSI(SCM_BIGNUM(arg0), SCM_INT_VALUE(arg1));
1774 }
1775 if (SCM_BIGNUMP(arg1)) {
1776 return Scm_BignumAdd(SCM_BIGNUM(arg0), SCM_BIGNUM(arg1));
1777 }
1778 if (SCM_RATNUMP(arg1)) {
1779 return Scm_RatnumAdd(arg0, arg1);
1780 }
1781 if (SCM_FLONUMP(arg1)) {
1782 double z = Scm_GetDouble(arg0) + SCM_FLONUM_VALUE(arg1);
1783 RETURN_FLONUM(z);
1784 }
1785 if (SCM_COMPNUMP(arg1)) {
1786 return Scm_MakeComplex(Scm_GetDouble(arg0)
1787 + SCM_COMPNUM_REAL(arg1),
1788 SCM_COMPNUM_IMAG(arg1));
1789 }
1790 /* fallback to generic */
1791 }
1792 else if (SCM_RATNUMP(arg0)) {
1793 if (SCM_INTP(arg1)) {
1794 if (SCM_EXACT_ZERO_P(arg1)) return arg0;
1795 return Scm_RatnumAdd(arg0, arg1);
1796 }
1797 if (SCM_BIGNUMP(arg1)||SCM_RATNUMP(arg1)) {
1798 return Scm_RatnumAdd(arg0, arg1);
1799 }
1800 if (SCM_FLONUMP(arg1)) {
1801 double z = Scm_GetDouble(arg0) + SCM_FLONUM_VALUE(arg1);
1802 RETURN_FLONUM(z);
1803 }
1804 if (SCM_COMPNUMP(arg1)) {
1805 return Scm_MakeComplex(Scm_GetDouble(arg0)
1806 + SCM_COMPNUM_REAL(arg1),
1807 SCM_COMPNUM_IMAG(arg1));
1808 }
1809 /* fallback to generic */
1810 }
1811 else if (SCM_FLONUMP(arg0)) {
1812 if (SCM_INTP(arg1)) {
1813 if (SCM_EXACT_ZERO_P(arg1)) return arg0;
1814 double z = SCM_FLONUM_VALUE(arg0) + (double)SCM_INT_VALUE(arg1);
1815 RETURN_FLONUM(z);
1816 }
1817 if (SCM_BIGNUMP(arg1) || SCM_RATNUMP(arg1)) {
1818 double z = SCM_FLONUM_VALUE(arg0) + Scm_GetDouble(arg1);
1819 RETURN_FLONUM(z);
1820 }
1821 if (SCM_FLONUMP(arg1)) {
1822 if (SCM_FLONUM_VALUE(arg0) == 0.0) return arg1;
1823 if (SCM_FLONUM_VALUE(arg1) == 0.0) return arg0;
1824 double z = SCM_FLONUM_VALUE(arg0) + SCM_FLONUM_VALUE(arg1);
1825 RETURN_FLONUM(z);
1826 }
1827 if (SCM_COMPNUMP(arg1)) {
1828 if (SCM_FLONUM_VALUE(arg0) == 0.0) return arg1;
1829 return Scm_MakeComplex(SCM_FLONUM_VALUE(arg0)
1830 + SCM_COMPNUM_REAL(arg1),
1831 SCM_COMPNUM_IMAG(arg1));
1832 }
1833 /* fallback to generic */
1834 }
1835 else if (SCM_COMPNUMP(arg0)) {
1836 if (SCM_INTP(arg1)) {
1837 if (SCM_EXACT_ZERO_P(arg1)) return arg0;
1838 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)
1839 + (double)SCM_INT_VALUE(arg1),
1840 SCM_COMPNUM_IMAG(arg0));
1841 }
1842 if (SCM_BIGNUMP(arg1) || SCM_RATNUMP(arg1)) {
1843 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)
1844 + Scm_GetDouble(arg1),
1845 SCM_COMPNUM_IMAG(arg0));
1846 }
1847 if (SCM_FLONUMP(arg1)) {
1848 if (SCM_FLONUM_VALUE(arg1) == 0.0) return arg0;
1849 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)
1850 + SCM_FLONUM_VALUE(arg1),
1851 SCM_COMPNUM_IMAG(arg0));
1852 }
1853 if (SCM_COMPNUMP(arg1)) {
1854 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)
1855 + SCM_COMPNUM_REAL(arg1),
1856 SCM_COMPNUM_IMAG(arg0)
1857 + SCM_COMPNUM_IMAG(arg1));
1858 }
1859 /* fallback to generic */
1860 }
1861 /* object-+ handling */
1862 SCM_FLONUM_ENSURE_MEM(arg0);
1863 SCM_FLONUM_ENSURE_MEM(arg1);
1864 return Scm_ApplyRec(SCM_OBJ(&generic_add), SCM_LIST2(arg0, arg1));
1865 }
DEFINE_DUAL_API2(Scm_Add,Scm_VMAdd,scm_add)1866 DEFINE_DUAL_API2(Scm_Add, Scm_VMAdd, scm_add)
1867
1868
1869 static ScmObj scm_sub(ScmObj arg0, ScmObj arg1, int vmp)
1870 {
1871 if (SCM_INTP(arg0)) {
1872 if (SCM_INTP(arg1)) {
1873 long r = SCM_INT_VALUE(arg0) - SCM_INT_VALUE(arg1);
1874 return Scm_MakeInteger(r);
1875 }
1876 if (SCM_BIGNUMP(arg1)) {
1877 ScmObj big = Scm_MakeBignumFromSI(SCM_INT_VALUE(arg0));
1878 return Scm_BignumSub(SCM_BIGNUM(big), SCM_BIGNUM(arg1));
1879 }
1880 if (SCM_RATNUMP(arg1)) {
1881 return Scm_RatnumSub(arg0, arg1);
1882 }
1883 if (SCM_FLONUMP(arg1)) {
1884 double z = (double)SCM_INT_VALUE(arg0) - SCM_FLONUM_VALUE(arg1);
1885 RETURN_FLONUM(z);
1886 }
1887 if (SCM_COMPNUMP(arg1)) {
1888 return Scm_MakeComplex((double)SCM_INT_VALUE(arg0)
1889 - SCM_COMPNUM_REAL(arg1),
1890 - SCM_COMPNUM_IMAG(arg1));
1891 }
1892 /* fallback to generic */
1893 }
1894 if (SCM_BIGNUMP(arg0)) {
1895 if (SCM_INTP(arg1)) {
1896 if (SCM_EXACT_ZERO_P(arg1)) return arg0;
1897 return Scm_BignumSubSI(SCM_BIGNUM(arg0), SCM_INT_VALUE(arg1));
1898 }
1899 if (SCM_BIGNUMP(arg1)) {
1900 return Scm_BignumSub(SCM_BIGNUM(arg0), SCM_BIGNUM(arg1));
1901 }
1902 if (SCM_RATNUMP(arg1)) {
1903 return Scm_RatnumSub(arg0, arg1);
1904 }
1905 if (SCM_FLONUMP(arg1)) {
1906 double z = Scm_GetDouble(arg0) - SCM_FLONUM_VALUE(arg1);
1907 RETURN_FLONUM(z);
1908 }
1909 if (SCM_COMPNUMP(arg1)) {
1910 return Scm_MakeComplex(Scm_GetDouble(arg0)
1911 - SCM_COMPNUM_REAL(arg1),
1912 - SCM_COMPNUM_IMAG(arg1));
1913 }
1914 /* fallback to generic */
1915 }
1916 if (SCM_RATNUMP(arg0)) {
1917 if (SCM_INTP(arg1)) {
1918 if (SCM_EXACT_ZERO_P(arg1)) return arg0;
1919 return Scm_RatnumSub(arg0, arg1);
1920 }
1921 if (SCM_BIGNUMP(arg1)||SCM_RATNUMP(arg1)) {
1922 return Scm_RatnumSub(arg0, arg1);
1923 }
1924 if (SCM_FLONUMP(arg1)) {
1925 if (SCM_FLONUM_VALUE(arg1) == 0.0) return arg0;
1926 double z = Scm_GetDouble(arg0) - SCM_FLONUM_VALUE(arg1);
1927 RETURN_FLONUM(z);
1928 }
1929 if (SCM_COMPNUMP(arg1)) {
1930 return Scm_MakeComplex(Scm_GetDouble(arg0)
1931 - SCM_COMPNUM_REAL(arg1),
1932 - SCM_COMPNUM_IMAG(arg1));
1933 }
1934 /* fallback to generic */
1935 }
1936 if (SCM_FLONUMP(arg0)) {
1937 if (SCM_INTP(arg1)) {
1938 if (SCM_EXACT_ZERO_P(arg1)) return arg0;
1939 double z = SCM_FLONUM_VALUE(arg0) - (double)SCM_INT_VALUE(arg1);
1940 RETURN_FLONUM(z);
1941 }
1942 if (SCM_BIGNUMP(arg1) || SCM_RATNUMP(arg1)) {
1943 double z = SCM_FLONUM_VALUE(arg0) - Scm_GetDouble(arg1);
1944 RETURN_FLONUM(z);
1945 }
1946 if (SCM_FLONUMP(arg1)) {
1947 if (SCM_FLONUM_VALUE(arg1) == 0.0) return arg0;
1948 double z = SCM_FLONUM_VALUE(arg0) - SCM_FLONUM_VALUE(arg1);
1949 RETURN_FLONUM(z);
1950 }
1951 if (SCM_COMPNUMP(arg1)) {
1952 return Scm_MakeComplex(SCM_FLONUM_VALUE(arg0)
1953 - SCM_COMPNUM_REAL(arg1),
1954 - SCM_COMPNUM_IMAG(arg1));
1955 }
1956 /* fallback to generic */
1957 }
1958 if (SCM_COMPNUMP(arg0)) {
1959 if (SCM_INTP(arg1)) {
1960 if (SCM_EXACT_ZERO_P(arg1)) return arg0;
1961 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)
1962 - (double)SCM_INT_VALUE(arg1),
1963 SCM_COMPNUM_IMAG(arg0));
1964 }
1965 if (SCM_BIGNUMP(arg1) || SCM_RATNUMP(arg1)) {
1966 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)
1967 - Scm_GetDouble(arg1),
1968 SCM_COMPNUM_IMAG(arg0));
1969 }
1970 if (SCM_FLONUMP(arg1)) {
1971 if (SCM_FLONUM_VALUE(arg1) == 0.0) return arg0;
1972 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)
1973 - Scm_GetDouble(arg1),
1974 SCM_COMPNUM_IMAG(arg0));
1975 }
1976 if (SCM_COMPNUMP(arg1)) {
1977 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)
1978 - SCM_COMPNUM_REAL(arg1),
1979 SCM_COMPNUM_IMAG(arg0)
1980 - SCM_COMPNUM_IMAG(arg1));
1981 }
1982 /* fallback to generic */
1983 }
1984 /* object-- handling */
1985 return Scm_ApplyRec(SCM_OBJ(&generic_sub), SCM_LIST2(arg0, arg1));
1986 }
DEFINE_DUAL_API2(Scm_Sub,Scm_VMSub,scm_sub)1987 DEFINE_DUAL_API2(Scm_Sub, Scm_VMSub, scm_sub)
1988
1989 /*
1990 * Multiplication
1991 */
1992
1993 static ScmObj scm_mul(ScmObj arg0, ScmObj arg1, int vmp)
1994 {
1995 if (SCM_INTP(arg0)) {
1996 if (SCM_INTP(arg1)) {
1997 long v0 = SCM_INT_VALUE(arg0);
1998 long v1 = SCM_INT_VALUE(arg1);
1999 long k;
2000 int ov;
2001 /* Using SMULOV to detect overflow portably. */
2002 SMULOV(k, ov, v0, v1);
2003 if (ov || !SCM_SMALL_INT_FITS(k)) {
2004 ScmObj big = Scm_MakeBignumFromSI(v0);
2005 return Scm_BignumMulSI(SCM_BIGNUM(big), v1);
2006 } else {
2007 return Scm_MakeInteger(k);
2008 }
2009 }
2010 if (SCM_BIGNUMP(arg1)) {
2011 if (SCM_EXACT_ZERO_P(arg0)) return arg0;
2012 if (SCM_EQ(arg0, SCM_MAKE_INT(1))) return arg1;
2013 return Scm_BignumMulSI(SCM_BIGNUM(arg1), SCM_INT_VALUE(arg0));
2014 }
2015 if (SCM_RATNUMP(arg1)) {
2016 if (SCM_EXACT_ZERO_P(arg0)) return arg0;
2017 if (SCM_EQ(arg0, SCM_MAKE_INT(1))) return arg1;
2018 return Scm_RatnumMul(arg0, arg1);
2019 }
2020 if (SCM_FLONUMP(arg1)) {
2021 if (SCM_EXACT_ZERO_P(arg0)) return arg0;
2022 if (SCM_EQ(arg0, SCM_MAKE_INT(1))) return arg1;
2023 double z = (double)SCM_INT_VALUE(arg0) * SCM_FLONUM_VALUE(arg1);
2024 RETURN_FLONUM(z);
2025 }
2026 if (SCM_COMPNUMP(arg1)) {
2027 if (SCM_EXACT_ZERO_P(arg0)) return arg0;
2028 if (SCM_EQ(arg0, SCM_MAKE_INT(1))) return arg1;
2029 return Scm_MakeComplex((double)SCM_INT_VALUE(arg0)
2030 * SCM_COMPNUM_REAL(arg1),
2031 (double)SCM_INT_VALUE(arg0)
2032 * SCM_COMPNUM_IMAG(arg1));
2033 }
2034 /* fallback to generic */
2035 }
2036 if (SCM_BIGNUMP(arg0)) {
2037 if (SCM_INTP(arg1)) {
2038 if (SCM_EXACT_ZERO_P(arg1)) return arg1;
2039 if (SCM_EXACT_ONE_P(arg1)) return arg0;
2040 return Scm_BignumMulSI(SCM_BIGNUM(arg0), SCM_INT_VALUE(arg1));
2041 }
2042 if (SCM_BIGNUMP(arg1)) {
2043 return Scm_BignumMul(SCM_BIGNUM(arg0), SCM_BIGNUM(arg1));
2044 }
2045 if (SCM_RATNUMP(arg1)) {
2046 return Scm_RatnumMul(arg0, arg1);
2047 }
2048 if (SCM_FLONUMP(arg1)) {
2049 double z = Scm_GetDouble(arg0) * SCM_FLONUM_VALUE(arg1);
2050 RETURN_FLONUM(z);
2051 }
2052 if (SCM_COMPNUMP(arg1)) {
2053 double z = Scm_GetDouble(arg0);
2054 return Scm_MakeComplex(z * SCM_COMPNUM_REAL(arg1),
2055 z * SCM_COMPNUM_IMAG(arg1));
2056 }
2057 /* fallback to generic */
2058 }
2059 if (SCM_RATNUMP(arg0)) {
2060 if (SCM_INTP(arg1)) {
2061 if (SCM_EXACT_ZERO_P(arg1)) return arg1;
2062 if (SCM_EXACT_ONE_P(arg1)) return arg0;
2063 return Scm_RatnumMul(arg0, arg1);
2064 }
2065 if (SCM_BIGNUMP(arg1)||SCM_RATNUMP(arg1)) {
2066 return Scm_RatnumMul(arg0, arg1);
2067 }
2068 if (SCM_FLONUMP(arg1)) {
2069 if (SCM_FLONUM_VALUE(arg1) == 0.0) return arg1;
2070 double z = Scm_GetDouble(arg0) * SCM_FLONUM_VALUE(arg1);
2071 RETURN_FLONUM(z);
2072 }
2073 if (SCM_COMPNUMP(arg1)) {
2074 return Scm_MakeComplex(Scm_GetDouble(arg0)
2075 * SCM_COMPNUM_REAL(arg1),
2076 Scm_GetDouble(arg0)
2077 * SCM_COMPNUM_IMAG(arg1));
2078 }
2079 /* fallback to generic */
2080 }
2081 if (SCM_FLONUMP(arg0)) {
2082 if (SCM_INTP(arg1)) {
2083 /* inexact number * exact zero makes exact zero */
2084 if (SCM_EXACT_ZERO_P(arg1)) return arg1;
2085 if (SCM_EXACT_ONE_P(arg1)) return arg0;
2086 double z = SCM_FLONUM_VALUE(arg0) * (double)SCM_INT_VALUE(arg1);
2087 RETURN_FLONUM(z);
2088 }
2089 if (SCM_BIGNUMP(arg1) || SCM_RATNUMP(arg1)) {
2090 RETURN_FLONUM(SCM_FLONUM_VALUE(arg0) * Scm_GetDouble(arg1));
2091 }
2092 if (SCM_FLONUMP(arg1)) {
2093 if (SCM_FLONUM_VALUE(arg1) == 1.0) return arg0;
2094 RETURN_FLONUM(SCM_FLONUM_VALUE(arg0) * SCM_FLONUM_VALUE(arg1));
2095 }
2096 if (SCM_COMPNUMP(arg1)) {
2097 return Scm_MakeComplex(SCM_FLONUM_VALUE(arg0)
2098 * SCM_COMPNUM_REAL(arg1),
2099 SCM_FLONUM_VALUE(arg0)
2100 * SCM_COMPNUM_IMAG(arg1));
2101 }
2102 /* fallback to generic */
2103 }
2104 if (SCM_COMPNUMP(arg0)) {
2105 if (SCM_INTP(arg1)) {
2106 if (SCM_EXACT_ZERO_P(arg1)) return arg1;
2107 if (SCM_EXACT_ONE_P(arg1)) return arg0;
2108 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)
2109 * (double)SCM_INT_VALUE(arg1),
2110 SCM_COMPNUM_IMAG(arg0)
2111 * (double)SCM_INT_VALUE(arg1));
2112 }
2113 if (SCM_BIGNUMP(arg1) || SCM_RATNUMP(arg1)) {
2114 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)
2115 * Scm_GetDouble(arg1),
2116 SCM_COMPNUM_IMAG(arg0)
2117 * Scm_GetDouble(arg1));
2118 }
2119 if (SCM_FLONUMP(arg1)) {
2120 if (SCM_FLONUM_VALUE(arg1) == 1.0) return arg0;
2121 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)
2122 * SCM_FLONUM_VALUE(arg1),
2123 SCM_COMPNUM_IMAG(arg0)
2124 * SCM_FLONUM_VALUE(arg1));
2125 }
2126 if (SCM_COMPNUMP(arg1)) {
2127 double r0 = SCM_COMPNUM_REAL(arg0);
2128 double i0 = SCM_COMPNUM_IMAG(arg0);
2129 double r1 = SCM_COMPNUM_REAL(arg1);
2130 double i1 = SCM_COMPNUM_IMAG(arg1);
2131 return Scm_MakeComplex(r0 * r1 - i0 * i1,
2132 r0 * i1 + r1 * i0);
2133 }
2134 /* fallback to generic */
2135 }
2136 SCM_FLONUM_ENSURE_MEM(arg0);
2137 SCM_FLONUM_ENSURE_MEM(arg1);
2138 return Scm_ApplyRec(SCM_OBJ(&generic_mul), SCM_LIST2(arg0, arg1));
2139 }
DEFINE_DUAL_API2(Scm_Mul,Scm_VMMul,scm_mul)2140 DEFINE_DUAL_API2(Scm_Mul, Scm_VMMul, scm_mul)
2141
2142 /*
2143 * Division
2144 */
2145
2146 /* We have three flavors of division API.
2147 * - Scm_Div : Scheme's `/'. Exact division produces exact result (maybe
2148 * rational)
2149 * - Scm_DivInexact : Scheme's `/.'. The result is always inexact. Fast.
2150 * - Scm_DivCompat : Scheme's `inexact-/'. This is only for the backward
2151 * compatibility, and probably we'll drop this in 1.0.
2152 * It works as Scm_Div, except that when Scm_Div produces ratnum,
2153 * Scm_DivCompat produces flonum.
2154 * There are also two 'VM' API, which can be used if you're returning
2155 * the value directly to the VM.
2156 * - Scm_VMDiv : The 'VM' version of Scm_Div.
2157 * - Scm_VMDivInexact : The 'VM' version of Scm_DivInexact.
2158 *
2159 * All these flavors are handled by a single function scm_div, with
2160 * three flags specifying the behavior.
2161 *
2162 * inexact compat vmp
2163 * Scm_Div FALSE FALSE FALSE
2164 * Scm_DivInexact TRUE FALSE FALSE
2165 * Scm_DivCompat d/c TRUE FALSE
2166 * Scm_VMDiv FALSE FALSE TRUE
2167 * Scm_VMDivInexact TRUE FALSE TRUE
2168 */
2169
2170 static ScmObj
2171 scm_div(ScmObj arg0, ScmObj arg1, int inexact, int compat, int vmp)
2172 {
2173 ScmObj r = SCM_UNBOUND;
2174
2175 #define SIMPLE_RETURN(x) do { r = (x); goto simple_return; } while (0)
2176
2177 if (SCM_INTP(arg0)) {
2178 if (SCM_INTP(arg1)) {
2179 if (SCM_EXACT_ZERO_P(arg1)) {
2180 if (inexact) goto anormal;
2181 else goto div_by_zero;
2182 }
2183 if (SCM_EXACT_ZERO_P(arg0)) SIMPLE_RETURN(arg0);
2184 if (SCM_EXACT_ONE_P(arg1)) SIMPLE_RETURN(arg0);
2185 if (compat) {
2186 if (SCM_INT_VALUE(arg0)%SCM_INT_VALUE(arg1) == 0) {
2187 long q = SCM_INT_VALUE(arg0)/SCM_INT_VALUE(arg1);
2188 return Scm_MakeInteger(q);
2189 } else {
2190 double z = (double)SCM_INT_VALUE(arg0)
2191 / (double)SCM_INT_VALUE(arg1);
2192 RETURN_FLONUM(z);
2193 }
2194 } else if (inexact) {
2195 double z = Scm_GetDouble(arg0)/Scm_GetDouble(arg1);
2196 RETURN_FLONUM(z);
2197 } else {
2198 return Scm_MakeRational(arg0, arg1);
2199 }
2200 }
2201 if (SCM_BIGNUMP(arg1)) {
2202 if (SCM_EXACT_ZERO_P(arg0)) SIMPLE_RETURN(arg0);
2203 goto ratnum_return;
2204 }
2205 if (SCM_RATNUMP(arg1)) {
2206 arg0 = Scm_Mul(arg0, SCM_RATNUM_DENOM(arg1));
2207 arg1 = SCM_RATNUM_NUMER(arg1);
2208 goto ratnum_return;
2209 }
2210 if (SCM_FLONUMP(arg1)) {
2211 if (SCM_FLONUM_VALUE(arg1) == 0.0) goto anormal;
2212 RETURN_FLONUM(SCM_INT_VALUE(arg0)/SCM_FLONUM_VALUE(arg1));
2213 }
2214 if (SCM_COMPNUMP(arg1)) {
2215 goto do_complex;
2216 }
2217 /* fallback to generic */
2218 }
2219 if (SCM_BIGNUMP(arg0)) {
2220 if (SCM_INTP(arg1)) {
2221 if (SCM_EXACT_ZERO_P(arg1)) {
2222 if (inexact) goto anormal;
2223 else goto div_by_zero;
2224 }
2225 if (SCM_EXACT_ONE_P(arg1)) SIMPLE_RETURN(arg0);
2226 goto ratnum_return;
2227 }
2228 if (SCM_BIGNUMP(arg1)) {
2229 goto ratnum_return;
2230 }
2231 if (SCM_RATNUMP(arg1)) {
2232 arg0 = Scm_Mul(arg0, SCM_RATNUM_DENOM(arg1));
2233 arg1 = SCM_RATNUM_NUMER(arg1);
2234 goto ratnum_return;
2235 }
2236 if (SCM_FLONUMP(arg1)) {
2237 if (SCM_FLONUM_VALUE(arg1) == 0.0) goto anormal;
2238 RETURN_FLONUM(Scm_GetDouble(arg0)/SCM_FLONUM_VALUE(arg1));
2239 }
2240 if (SCM_COMPNUMP(arg1)) {
2241 goto do_complex;
2242 }
2243 /* fallback to generic */
2244 }
2245 if (SCM_RATNUMP(arg0)) {
2246 if (SCM_INTP(arg1)) {
2247 if (SCM_EXACT_ZERO_P(arg1)) {
2248 if (inexact) goto anormal;
2249 else goto div_by_zero;
2250 }
2251 if (SCM_EXACT_ONE_P(arg1)) SIMPLE_RETURN(arg0);
2252 arg1 = Scm_Mul(SCM_RATNUM_DENOM(arg0), arg1);
2253 arg0 = SCM_RATNUM_NUMER(arg0);
2254 goto ratnum_return;
2255 }
2256 if (SCM_BIGNUMP(arg1)) {
2257 arg1 = Scm_Mul(SCM_RATNUM_DENOM(arg0), arg1);
2258 arg0 = SCM_RATNUM_NUMER(arg0);
2259 goto ratnum_return;
2260 }
2261 if (SCM_RATNUMP(arg1)) {
2262 if (!compat && !inexact) {
2263 return Scm_RatnumDiv(arg0, arg1);
2264 } else {
2265 ScmObj numer = Scm_Mul(SCM_RATNUM_NUMER(arg0),
2266 SCM_RATNUM_DENOM(arg1));
2267 ScmObj denom = Scm_Mul(SCM_RATNUM_DENOM(arg0),
2268 SCM_RATNUM_NUMER(arg1));
2269 arg0 = numer;
2270 arg1 = denom;
2271 goto ratnum_return;
2272 }
2273 }
2274 if (SCM_FLONUMP(arg1)) {
2275 if (SCM_FLONUM_VALUE(arg1) == 0.0) goto anormal;
2276 RETURN_FLONUM(Scm_GetDouble(arg0)/SCM_FLONUM_VALUE(arg1));
2277 }
2278 if (SCM_COMPNUMP(arg1)) {
2279 goto do_complex;
2280 }
2281 /* fallback to generic */
2282 }
2283 if (SCM_FLONUMP(arg0)) {
2284 if (SCM_INTP(arg1)) {
2285 if (SCM_IS_NAN(SCM_FLONUM_VALUE(arg0))) return SCM_NAN;
2286 if (SCM_EXACT_ZERO_P(arg1)) goto anormal;
2287 if (SCM_EXACT_ONE_P(arg1)) return arg0;
2288 RETURN_FLONUM(SCM_FLONUM_VALUE(arg0)/SCM_INT_VALUE(arg1));
2289 }
2290 if (SCM_BIGNUMP(arg1) || SCM_RATNUMP(arg1)) {
2291 RETURN_FLONUM(SCM_FLONUM_VALUE(arg0)/Scm_GetDouble(arg1));
2292 }
2293 if (SCM_FLONUMP(arg1)) {
2294 if (SCM_IS_NAN(SCM_FLONUM_VALUE(arg0))) return SCM_NAN;
2295 if (SCM_FLONUM_VALUE(arg1) == 0.0) goto anormal;
2296 RETURN_FLONUM(SCM_FLONUM_VALUE(arg0)/SCM_FLONUM_VALUE(arg1));
2297 }
2298 if (SCM_COMPNUMP(arg1)) {
2299 goto do_complex;
2300 }
2301 /* fallback to generic */
2302 }
2303 if (SCM_COMPNUMP(arg0)) {
2304 if (SCM_INTP(arg1)) {
2305 /* NB: Gauche has no exact compnum */
2306 if (SCM_EXACT_ZERO_P(arg1)) goto anormal_comp;
2307 if (SCM_EXACT_ONE_P(arg1)) return arg0;
2308 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)/SCM_INT_VALUE(arg1),
2309 SCM_COMPNUM_IMAG(arg0)/SCM_INT_VALUE(arg1));
2310 }
2311 if (SCM_BIGNUMP(arg1) || SCM_RATNUMP(arg1)) {
2312 double z = Scm_GetDouble(arg1);
2313 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)/z,
2314 SCM_COMPNUM_IMAG(arg0)/z);
2315 }
2316 if (SCM_FLONUMP(arg1)) {
2317 if (SCM_FLONUM_VALUE(arg1) == 0.0) goto anormal_comp;
2318 return Scm_MakeComplex(SCM_COMPNUM_REAL(arg0)/SCM_FLONUM_VALUE(arg1),
2319 SCM_COMPNUM_IMAG(arg0)/SCM_FLONUM_VALUE(arg1));
2320 }
2321 if (SCM_COMPNUMP(arg1)) {
2322 double r0 = SCM_COMPNUM_REAL(arg0);
2323 double i0 = SCM_COMPNUM_IMAG(arg0);
2324 double r1 = SCM_COMPNUM_REAL(arg1);
2325 double i1 = SCM_COMPNUM_IMAG(arg1);
2326 double d = r1*r1+i1*i1;
2327 return Scm_MakeComplex((r0*r1 + i0*i1)/d,
2328 (i0*r1 - r0*i1)/d);
2329 }
2330 /* fallback to generic */
2331 }
2332 SCM_FLONUM_ENSURE_MEM(arg0);
2333 SCM_FLONUM_ENSURE_MEM(arg1);
2334 return Scm_ApplyRec(SCM_OBJ(&generic_div), SCM_LIST2(arg0, arg1));
2335
2336 ratnum_return:
2337 {
2338 /* arg0 and arg1 contains exact numbers.*/
2339 if (compat) goto compat_return;
2340 if (inexact) goto inexact_return;
2341 return Scm_MakeRational(arg0, arg1);
2342 }
2343 compat_return:
2344 {
2345 /* We have exact integer division arg0/arg1 (arg1 != 0).
2346 If it doesn't produce a whole integer, we coerce the
2347 result to flonum. */
2348 ScmObj rem;
2349 ScmObj q = Scm_Quotient(arg0, arg1, &rem);
2350 if (SCM_EXACT_ZERO_P(rem)) {
2351 return q;
2352 }
2353 /*FALLTHROUGH*/
2354 }
2355 inexact_return:
2356 {
2357 double numer = Scm_GetDouble(arg0);
2358 double denom = Scm_GetDouble(arg1);
2359 if (SCM_IS_INF(numer) || SCM_IS_INF(denom)) {
2360 /* special path - we need more sophisticated calculaton. */
2361 ScmObj r = Scm_MakeRational(arg0, arg1);
2362 RETURN_FLONUM(Scm_GetDouble(r));
2363 } else {
2364 RETURN_FLONUM(numer/denom);
2365 }
2366 }
2367 simple_return:
2368 {
2369 if (inexact) return Scm_Inexact(r);
2370 else return r;
2371 }
2372 div_by_zero:
2373 {
2374 Scm_Error("attempt to calculate a division by zero");
2375 }
2376 anormal:
2377 /* real inexact division by zero */
2378 {
2379 int s0 = Scm_Sign(arg0);
2380 int s1 = SCM_FLONUMP(arg1)? Scm_FlonumSign(SCM_FLONUM_VALUE(arg1)) : 1;
2381 if (s0 == 0) return SCM_NAN;
2382 if (s0*s1 < 0) return SCM_NEGATIVE_INFINITY;
2383 else return SCM_POSITIVE_INFINITY;
2384 }
2385 anormal_comp:
2386 /* complex inexact division by zero */
2387 {
2388 double r0 = SCM_COMPNUM_REAL(arg0);
2389 double i0 = SCM_COMPNUM_IMAG(arg0);
2390 int s1 = SCM_FLONUMP(arg1)? Scm_FlonumSign(SCM_FLONUM_VALUE(arg1)) : 1;
2391 double r =
2392 SCM_IS_NAN(r0) ? SCM_DBL_NAN
2393 : ((r0*s1 > 0.0) ? SCM_DBL_POSITIVE_INFINITY
2394 : ((r0*s1 < 0.0) ? SCM_DBL_NEGATIVE_INFINITY
2395 : SCM_DBL_NAN));
2396 double i =
2397 SCM_IS_NAN(i0) ? SCM_DBL_NAN
2398 : ((i0*s1 > 0.0) ? SCM_DBL_POSITIVE_INFINITY
2399 : ((i0*s1 < 0.0) ? SCM_DBL_NEGATIVE_INFINITY
2400 : SCM_DBL_NAN));
2401 return Scm_MakeComplex(r, i);
2402 }
2403 do_complex:
2404 {
2405 double r1 = SCM_COMPNUM_REAL(arg1);
2406 double i1 = SCM_COMPNUM_IMAG(arg1);
2407 double d = r1*r1+i1*i1;
2408 return Scm_MakeComplex(r1 * Scm_GetDouble(arg0) / d,
2409 -i1 * Scm_GetDouble(arg0) / d);
2410 }
2411 }
2412
Scm_Div(ScmObj x,ScmObj y)2413 ScmObj Scm_Div(ScmObj x, ScmObj y)
2414 {
2415 return scm_div(x, y, FALSE, FALSE, FALSE);
2416 }
2417
Scm_DivInexact(ScmObj x,ScmObj y)2418 ScmObj Scm_DivInexact(ScmObj x, ScmObj y)
2419 {
2420 return scm_div(x, y, TRUE, FALSE, FALSE);
2421 }
2422
Scm_DivCompat(ScmObj x,ScmObj y)2423 ScmObj Scm_DivCompat(ScmObj x, ScmObj y)
2424 {
2425 return scm_div(x, y, FALSE, TRUE, FALSE);
2426 }
2427
Scm_VMDiv(ScmObj x,ScmObj y)2428 ScmObj Scm_VMDiv(ScmObj x, ScmObj y)
2429 {
2430 return scm_div(x, y, FALSE, FALSE, TRUE);
2431 }
2432
Scm_VMDivInexact(ScmObj x,ScmObj y)2433 ScmObj Scm_VMDivInexact(ScmObj x, ScmObj y)
2434 {
2435 return scm_div(x, y, TRUE, FALSE, TRUE);
2436 }
2437
2438
2439 /*
2440 * Integer division
2441 * Returns (quotient x y)
2442 * If rem != NULL, sets *rem to be (remainder x y) as well.
2443 * We don't provide Scm_VMQuotient, assuming passing flonums to
2444 * quotient is rare.
2445 */
Scm_Quotient(ScmObj x,ScmObj y,ScmObj * rem)2446 ScmObj Scm_Quotient(ScmObj x, ScmObj y, ScmObj *rem)
2447 {
2448 double rx, ry;
2449
2450 /* Trivial shortcut. This case may seem too specific, but actually
2451 it appears rather often in rational operations. */
2452 if (SCM_EQ(y, SCM_MAKE_INT(1))) {
2453 if (!Scm_IntegerP(x)) goto BADARG;
2454 if (rem) *rem = SCM_MAKE_INT(0);
2455 return x;
2456 }
2457
2458 if (SCM_INTP(x)) {
2459 if (SCM_INTP(y)) {
2460 if (SCM_INT_VALUE(y) == 0) goto DIVBYZERO;
2461 long q = SCM_INT_VALUE(x)/SCM_INT_VALUE(y);
2462 if (rem) {
2463 long r = SCM_INT_VALUE(x)%SCM_INT_VALUE(y);
2464 *rem = SCM_MAKE_INT(r);
2465 }
2466 return SCM_MAKE_INT(q);
2467 }
2468 if (SCM_BIGNUMP(y)) {
2469 if (rem) *rem = x;
2470 return SCM_MAKE_INT(0);
2471 }
2472 if (SCM_FLONUMP(y)) {
2473 rx = (double)SCM_INT_VALUE(x);
2474 ry = SCM_FLONUM_VALUE(y);
2475 if (ry != floor(ry)) goto BADARGY;
2476 goto DO_FLONUM;
2477 }
2478 goto BADARGY;
2479 } else if (SCM_BIGNUMP(x)) {
2480 if (SCM_INTP(y)) {
2481 long r;
2482 ScmObj q = Scm_BignumDivSI(SCM_BIGNUM(x), SCM_INT_VALUE(y), &r);
2483 if (rem) *rem = SCM_MAKE_INT(r);
2484 return q;
2485 } else if (SCM_BIGNUMP(y)) {
2486 ScmObj qr = Scm_BignumDivRem(SCM_BIGNUM(x), SCM_BIGNUM(y));
2487 if (rem) *rem = SCM_CDR(qr);
2488 return SCM_CAR(qr);
2489 } else if (SCM_FLONUMP(y)) {
2490 rx = Scm_BignumToDouble(SCM_BIGNUM(x));
2491 ry = SCM_FLONUM_VALUE(y);
2492 if (ry != floor(ry)) goto BADARGY;
2493 goto DO_FLONUM;
2494 }
2495 goto BADARGY;
2496 } else if (SCM_FLONUMP(x)) {
2497 rx = SCM_FLONUM_VALUE(x);
2498 if (rx != floor(rx)) goto BADARG;
2499 if (SCM_INTP(y)) {
2500 ry = (double)SCM_INT_VALUE(y);
2501 } else if (SCM_BIGNUMP(y)) {
2502 ry = Scm_BignumToDouble(SCM_BIGNUM(y));
2503 } else if (SCM_FLONUMP(y)) {
2504 ry = SCM_FLONUM_VALUE(y);
2505 if (ry != floor(ry)) goto BADARGY;
2506 } else {
2507 goto BADARGY;
2508 }
2509 DO_FLONUM:
2510 {
2511 if (ry == 0.0) goto DIVBYZERO;
2512 double q = (rx*ry > 0)? floor(rx/ry) : ceil(rx/ry);
2513 if (rem) {
2514 double rr = roundeven(rx - q*ry);
2515 *rem = Scm_MakeFlonum(rr);
2516 }
2517 return Scm_MakeFlonum(q);
2518 }
2519 } else {
2520 goto BADARG;
2521 }
2522 DIVBYZERO:
2523 Scm_Error("attempt to calculate a quotient by zero");
2524 BADARGY:
2525 x = y;
2526 BADARG:
2527 SCM_FLONUM_ENSURE_MEM(x);
2528 Scm_Error("integer required, but got %S", x);
2529 return SCM_UNDEFINED; /* dummy */
2530 }
2531
2532 /* Modulo and Reminder.
2533 We don't provide Scm_VMModulo, assuming passing flonums to modulo and
2534 remainder is rare.
2535 TODO: on gcc, % works like reminder. I'm not sure the exact behavior
2536 of % is defined in ANSI C. Need to check it later. */
Scm_Modulo(ScmObj x,ScmObj y,int remp)2537 ScmObj Scm_Modulo(ScmObj x, ScmObj y, int remp)
2538 {
2539 double rx, ry;
2540 if (SCM_INTP(x)) {
2541 if (SCM_INTP(y)) {
2542 if (SCM_INT_VALUE(y) == 0) goto DIVBYZERO;
2543 long r = SCM_INT_VALUE(x)%SCM_INT_VALUE(y);
2544 if (!remp && r) {
2545 if ((SCM_INT_VALUE(x) > 0 && SCM_INT_VALUE(y) < 0)
2546 || (SCM_INT_VALUE(x) < 0 && SCM_INT_VALUE(y) > 0)) {
2547 r += SCM_INT_VALUE(y);
2548 }
2549 }
2550 return SCM_MAKE_INT(r);
2551 }
2552 if (SCM_BIGNUMP(y)) {
2553 if (remp) {
2554 return x;
2555 } else {
2556 if ((SCM_INT_VALUE(x) < 0 && SCM_BIGNUM_SIGN(y) > 0)
2557 || (SCM_INT_VALUE(x) > 0 && SCM_BIGNUM_SIGN(y) < 0)) {
2558 return Scm_BignumAddSI(SCM_BIGNUM(y), SCM_INT_VALUE(x));
2559 } else {
2560 return x;
2561 }
2562 }
2563 }
2564 rx = (double)SCM_INT_VALUE(x);
2565 if (SCM_FLONUMP(y)) {
2566 ry = SCM_FLONUM_VALUE(y);
2567 if (ry != floor(ry)) goto BADARGY;
2568 goto DO_FLONUM;
2569 }
2570 goto BADARGY;
2571 } else if (SCM_BIGNUMP(x)) {
2572 if (SCM_INTP(y)) {
2573 long iy = SCM_INT_VALUE(y);
2574 long rem = Scm_BignumRemSI(SCM_BIGNUM(x), iy);
2575 if (!remp
2576 && rem
2577 && ((SCM_BIGNUM_SIGN(x) < 0 && iy > 0)
2578 || (SCM_BIGNUM_SIGN(x) > 0 && iy < 0))) {
2579 return SCM_MAKE_INT(iy + rem);
2580 }
2581 return SCM_MAKE_INT(rem);
2582 }
2583 if (SCM_BIGNUMP(y)) {
2584 ScmObj rem = SCM_CDR(Scm_BignumDivRem(SCM_BIGNUM(x), SCM_BIGNUM(y)));
2585 if (!remp
2586 && (rem != SCM_MAKE_INT(0))
2587 && (SCM_BIGNUM_SIGN(x) * SCM_BIGNUM_SIGN(y) < 0)) {
2588 if (SCM_BIGNUMP(rem)) {
2589 return Scm_BignumAdd(SCM_BIGNUM(y), SCM_BIGNUM(rem));
2590 } else {
2591 return Scm_BignumAddSI(SCM_BIGNUM(y), SCM_INT_VALUE(rem));
2592 }
2593 }
2594 return rem;
2595 }
2596 rx = Scm_BignumToDouble(SCM_BIGNUM(x));
2597 if (SCM_FLONUMP(y)) {
2598 ry = SCM_FLONUM_VALUE(y);
2599 if (ry != floor(ry)) goto BADARGY;
2600 goto DO_FLONUM;
2601 }
2602 goto BADARGY;
2603 } else if (SCM_FLONUMP(x)) {
2604 double rem;
2605 rx = SCM_FLONUM_VALUE(x);
2606 if (rx != floor(rx)) goto BADARG;
2607 if (SCM_INTP(y)) {
2608 ry = (double)SCM_INT_VALUE(y);
2609 } else if (SCM_BIGNUMP(y)) {
2610 ry = Scm_BignumToDouble(SCM_BIGNUM(y));
2611 } else if (SCM_FLONUMP(y)) {
2612 ry = SCM_FLONUM_VALUE(y);
2613 if (ry != floor(ry)) goto BADARGY;
2614 } else {
2615 goto BADARGY;
2616 }
2617 DO_FLONUM:
2618 if (ry == 0.0) goto DIVBYZERO;
2619 rem = fmod(rx, ry);
2620 if (!remp && rem != 0.0) {
2621 if ((rx > 0 && ry < 0) || (rx < 0 && ry > 0)) {
2622 rem += ry;
2623 }
2624 }
2625 return Scm_MakeFlonum(rem);
2626 } else {
2627 goto BADARG;
2628 }
2629 DIVBYZERO:
2630 Scm_Error("attempt to take a modulo or remainder by zero");
2631 BADARGY:
2632 x = y;
2633 BADARG:
2634 SCM_FLONUM_ENSURE_MEM(x);
2635 Scm_Error("integer required, but got %S", x);
2636 return SCM_UNDEFINED; /* dummy */
2637 }
2638
2639 /*
2640 * Gcd
2641 */
2642
2643 /* assumes x > y >= 0 */
gcd_fixfix(u_long x,u_long y)2644 static u_long gcd_fixfix(u_long x, u_long y)
2645 {
2646 while (y > 0) {
2647 u_long r = x % y;
2648 x = y;
2649 y = r;
2650 }
2651 return x;
2652 }
2653
gcd_floflo(double x,double y)2654 static double gcd_floflo(double x, double y)
2655 {
2656 if (x < 0) x = -x;
2657 if (y < 0) y = -y;
2658 if (x < y) { double t = x; x = y; y = t; }
2659
2660 while (y > 0.0) {
2661 double r = fmod(x, y);
2662 x = y;
2663 y = r;
2664 }
2665 return x;
2666 }
2667
2668 /* assumes y <= LONG_MAX. curiously, the sign of x doesn't matter,
2669 since it only affects the remainder's sign which we adjust afterwards. */
gcd_bigfix(ScmBignum * x,u_long y)2670 static u_long gcd_bigfix(ScmBignum *x, u_long y)
2671 {
2672 long rem;
2673 (void)Scm_BignumDivSI(x, (signed long)y, &rem);
2674 if (rem < 0) rem = -rem;
2675 return gcd_fixfix(y, (u_long)rem);
2676 }
2677
2678 /* We don't provide Scm_VMGcd, assuming passing flonums to gcd is rare. */
Scm_Gcd(ScmObj x,ScmObj y)2679 ScmObj Scm_Gcd(ScmObj x, ScmObj y)
2680 {
2681 if (!Scm_IntegerP(x)) {
2682 Scm_Error("integer required, but got %S", x);
2683 }
2684 if (!Scm_IntegerP(y)) {
2685 Scm_Error("integer required, but got %S", y);
2686 }
2687 if (SCM_FLONUMP(x) || SCM_FLONUMP(y)) {
2688 return Scm_MakeFlonum(gcd_floflo(Scm_GetDouble(x), Scm_GetDouble(y)));
2689 }
2690
2691 if (SCM_EXACT_ZERO_P(x)) return y;
2692 if (SCM_EXACT_ZERO_P(y)) return x;
2693
2694 int ox = FALSE, oy = FALSE;
2695 long ix = Scm_GetIntegerClamp(x, SCM_CLAMP_NONE, &ox);
2696 long iy = Scm_GetIntegerClamp(y, SCM_CLAMP_NONE, &oy);
2697
2698 if (!ox && !oy) {
2699 u_long ux = (ix < 0)? -ix : ix;
2700 u_long uy = (iy < 0)? -iy : iy;
2701 u_long ur = (ux >= uy)? gcd_fixfix(ux, uy) : gcd_fixfix(uy, ux);
2702 return Scm_MakeIntegerU(ur);
2703 }
2704
2705 if (!oy && iy != LONG_MIN) {
2706 /* x overflows long. y doesn't. so we know abs(x) > abs(y)
2707 (abs(x) == abs(y) iff LONG_MAX+1 and y == LONG_MIN, but we've
2708 excluded it above). */
2709 SCM_ASSERT(SCM_BIGNUMP(x));
2710 u_long uy = (iy < 0)? -iy : iy;
2711 u_long ur = gcd_bigfix(SCM_BIGNUM(x), uy);
2712 return Scm_MakeIntegerU(ur);
2713 }
2714
2715 if (!ox && ix != LONG_MIN) {
2716 /* reverse condition of above */
2717 SCM_ASSERT(SCM_BIGNUMP(y));
2718 u_long ux = (ix < 0)? -ix : ix;
2719 u_long ur = gcd_bigfix(SCM_BIGNUM(y), ux);
2720 return Scm_MakeIntegerU(ur);
2721 }
2722
2723 /* Now we need to treat both args as bignums. We could use
2724 Algorithm L in Knuth's TAOCP 4.5.2, but we assume this path
2725 is rarely executed, so we don't bother for now. */
2726 x = Scm_Abs(x);
2727 y = Scm_Abs(y);
2728 if (Scm_NumCmp(x, y) < 0) {ScmObj t = x; x = y; y = t;}
2729
2730 while (!SCM_EXACT_ZERO_P(y)) {
2731 ScmObj r = Scm_Modulo(x, y, TRUE);
2732 x = y;
2733 y = r;
2734 }
2735 return x;
2736 }
2737
2738 /*===============================================================
2739 * Exponential and trigometric functions
2740 */
2741
2742 /* Integer power of 10. It is extensively used during string->number
2743 and number->string operations.
2744 IEXPT10_TABLESIZ is ceil(-log10(ldexp(1.0, -1022-52))) + 2 */
2745 /* NB: actually we need more margin here to handle denormalized numbers. */
2746 #define IEXPT10_TABLESIZ 341
2747 static ScmObj iexpt10_n[IEXPT10_TABLESIZ] = { NULL };
2748 static int iexpt10_initialized = FALSE;
2749
iexpt10_init(void)2750 static void iexpt10_init(void)
2751 {
2752 iexpt10_n[0] = SCM_MAKE_INT(1);
2753 iexpt10_n[1] = SCM_MAKE_INT(10);
2754 iexpt10_n[2] = SCM_MAKE_INT(100);
2755 iexpt10_n[3] = SCM_MAKE_INT(1000);
2756 iexpt10_n[4] = SCM_MAKE_INT(10000);
2757 iexpt10_n[5] = SCM_MAKE_INT(100000);
2758 iexpt10_n[6] = SCM_MAKE_INT(1000000);
2759 for (int i=7; i<IEXPT10_TABLESIZ; i++) {
2760 iexpt10_n[i] = Scm_Mul(iexpt10_n[i-1], SCM_MAKE_INT(10));
2761 }
2762 iexpt10_initialized = TRUE;
2763 }
2764
2765 #define IEXPT10_INIT() \
2766 do { if (!iexpt10_initialized) iexpt10_init(); } while (0)
2767
2768 /* expt(x, y) where x is exact and y is integer */
Scm_ExactIntegerExpt(ScmObj x,ScmObj y)2769 ScmObj Scm_ExactIntegerExpt(ScmObj x, ScmObj y)
2770 {
2771 int sign = Scm_Sign(y);
2772 ScmObj r = SCM_MAKE_INT(1);
2773
2774 if (sign == 0) return r;
2775 if (SCM_EQ(x, SCM_MAKE_INT(1))) return r;
2776 if (SCM_EQ(x, SCM_MAKE_INT(-1))) return Scm_OddP(y)? SCM_MAKE_INT(-1) : r;
2777
2778 if (!SCM_INTP(y)) {
2779 /* who wants such a heavy calculation? */
2780 Scm_Error("exponent too big: %S", y);
2781 }
2782 long iy = SCM_INT_VALUE(y);
2783 /* Shortcut for special cases */
2784 if (SCM_EQ(x, SCM_MAKE_INT(10)) && iy > 0 && iy < IEXPT10_TABLESIZ) {
2785 /* We have a precalculated table for 10^y */
2786 IEXPT10_INIT();
2787 r = iexpt10_n[iy];
2788 } else if (SCM_EQ(x, SCM_MAKE_INT(2)) && iy > 0) {
2789 /* Use shift operation for 2^y, y>0 */
2790 r = Scm_Ash(SCM_MAKE_INT(1), iy);
2791 } else {
2792 /* General case */
2793 if (iy < 0) iy = -iy;
2794 for (;;) {
2795 if (iy == 0) break;
2796 if (iy == 1) { r = Scm_Mul(r, x); break; }
2797 if (iy & 0x01) r = Scm_Mul(r, x);
2798 x = Scm_Mul(x, x);
2799 iy >>= 1;
2800 }
2801 }
2802 return (sign < 0)? Scm_Reciprocal(r) : r;
2803 }
2804
scm_expt(ScmObj x,ScmObj y,int vmp)2805 static ScmObj scm_expt(ScmObj x, ScmObj y, int vmp)
2806 {
2807 /* NB: The exact case is handled by expt in libnum.scm; we check this case
2808 just for the backward compatibility. */
2809 if (SCM_EXACTP(x) && SCM_INTEGERP(y)) return Scm_ExactIntegerExpt(x, y);
2810 if (!SCM_REALP(x)) Scm_Error("real number required, but got %S", x);
2811 if (!SCM_REALP(y)) Scm_Error("real number required, but got %S", y);
2812 double dx = Scm_GetDouble(x);
2813 double dy = Scm_GetDouble(y);
2814 if (dy == 0.0) {
2815 RETURN_FLONUM(1.0);
2816 } else if (dx < 0 && !Scm_IntegerP(y)) {
2817 /* x^y == exp(y * log(x)) = exp(y*log(|x|))*exp(y*arg(x)*i)
2818 if x is a negative real number, arg(x) == pi
2819 */
2820 double mag = exp(dy * log(-dx));
2821 return Scm_MakeComplex(mag * Scm_CosPi(dy), mag * Scm_SinPi(dy));
2822 } else {
2823 RETURN_FLONUM(pow(dx, dy));
2824 }
2825 }
DEFINE_DUAL_API2(Scm_Expt,Scm_VMExpt,scm_expt)2826 DEFINE_DUAL_API2(Scm_Expt, Scm_VMExpt, scm_expt)
2827
2828 /* If num is exact 2^s (s >= 0), returns s. Otherwise returns -1. */
2829 long Scm_TwosPower(ScmObj n)
2830 {
2831 if (SCM_INTP(n)) {
2832 long i = SCM_INT_VALUE(n);
2833 if (i <= 0) return -1;
2834 if ((i<<1) == ((i ^ (i-1)) + 1)) {
2835 return Scm__HighestBitNumber(i);
2836 }
2837 /*FALTHROUGH*/
2838 } else if (SCM_BIGNUMP(n) && SCM_BIGNUM_SIGN(n) > 0) {
2839 ScmBits *b = (ScmBits*)SCM_BIGNUM(n)->values;
2840 int l = SCM_BIGNUM_SIZE(n) * SCM_WORD_BITS;
2841 int c = Scm_BitsLowest1(b, 0, l);
2842 if (c == Scm_BitsHighest1(b, 0, l)) return c;
2843 /*FALTHROUGH*/
2844 }
2845 return -1;
2846 }
2847
2848 /* sinpi(x) = sin(x * pi)
2849 cospi(x) = cos(x * pi)
2850 tanpi(x) = tan(x * pi)
2851
2852 We first reduce input range to -1 <= x <= 1 by trig_pi_reduce_range.
2853 */
trig_pi_reduce_range(double x)2854 static double trig_pi_reduce_range(double x)
2855 {
2856 double xx = fmod(x, 2.0); /* -2.0 < x < 2.0 */
2857 if (xx > 1.0) return xx - 2.0;
2858 if (xx < -1.0) return xx + 2.0;
2859 if (xx == 0.0 && signbit(xx)) return -xx; /* we don't return -0.0 */
2860 else return xx;
2861 }
2862
Scm_SinPi(double x)2863 double Scm_SinPi(double x)
2864 {
2865 double xx = trig_pi_reduce_range(x);
2866 if (xx >= 0) {
2867 if (xx > 0.5) xx = 1 - xx;
2868 if (xx > 0.25) return cos(M_PI*(0.5-xx));
2869 else return sin(M_PI*xx);
2870 } else {
2871 if (xx < -0.5) xx = -1 - xx;
2872 if (xx < -0.25) return -cos(M_PI*(-0.5-xx));
2873 else return sin(M_PI*xx);
2874 }
2875 }
2876
Scm_CosPi(double x)2877 double Scm_CosPi(double x)
2878 {
2879 double xx = fabs(trig_pi_reduce_range(x));
2880 if (xx >= 0.75) return -cos(M_PI*(1-xx));
2881 if (xx > 0.25) return sin(M_PI*(0.5-xx));
2882 else return cos(M_PI*xx);
2883 }
2884
Scm_TanPi(double x)2885 double Scm_TanPi(double x)
2886 {
2887 return Scm_SinPi(x)/Scm_CosPi(x);
2888 }
2889
2890 /*===============================================================
2891 * Comparison
2892 */
2893
nan_p(ScmObj arg)2894 static inline int nan_p(ScmObj arg)
2895 {
2896 return (SCM_FLONUMP(arg) && SCM_IS_NAN(SCM_FLONUM_VALUE(arg)));
2897 }
2898
either_nan_p(ScmObj arg0,ScmObj arg1)2899 static inline int either_nan_p(ScmObj arg0, ScmObj arg1)
2900 {
2901 return (nan_p(arg0) || nan_p(arg1));
2902 }
2903
Scm_NumEq(ScmObj arg0,ScmObj arg1)2904 int Scm_NumEq(ScmObj arg0, ScmObj arg1)
2905 {
2906 if (SCM_COMPNUMP(arg0)) {
2907 if (SCM_COMPNUMP(arg1)) {
2908 return ((SCM_COMPNUM_REAL(arg0) == SCM_COMPNUM_REAL(arg1))
2909 && (SCM_COMPNUM_IMAG(arg0) == SCM_COMPNUM_IMAG(arg1)));
2910 }
2911 return FALSE;
2912 } else {
2913 if (SCM_COMPNUMP(arg1)) return FALSE;
2914 if (either_nan_p(arg0, arg1)) return FALSE;
2915 return (Scm_NumCmp(arg0, arg1) == 0);
2916 }
2917 }
2918
Scm_NumLT(ScmObj arg0,ScmObj arg1)2919 int Scm_NumLT(ScmObj arg0, ScmObj arg1)
2920 {
2921 if (either_nan_p(arg0, arg1)) return FALSE;
2922 return Scm_NumCmp(arg0, arg1) < 0;
2923 }
2924
Scm_NumLE(ScmObj arg0,ScmObj arg1)2925 int Scm_NumLE(ScmObj arg0, ScmObj arg1)
2926 {
2927 if (either_nan_p(arg0, arg1)) return FALSE;
2928 return Scm_NumCmp(arg0, arg1) <= 0;
2929 }
2930
Scm_NumGT(ScmObj arg0,ScmObj arg1)2931 int Scm_NumGT(ScmObj arg0, ScmObj arg1)
2932 {
2933 if (either_nan_p(arg0, arg1)) return FALSE;
2934 return Scm_NumCmp(arg0, arg1) > 0;
2935 }
2936
Scm_NumGE(ScmObj arg0,ScmObj arg1)2937 int Scm_NumGE(ScmObj arg0, ScmObj arg1)
2938 {
2939 if (either_nan_p(arg0, arg1)) return FALSE;
2940 return Scm_NumCmp(arg0, arg1) >= 0;
2941 }
2942
2943 /* 2-arg comparison.
2944 NB: This routine is called from VM or built-in SUBRs when arg0 and/or
2945 arg1 contain register flonums. That's why we insert SCM_FLONUM_ENSURE_MEM
2946 in the path to Scm_Error at the end. For all other paths we know the
2947 register flonum will never leak out. If you make changes here,
2948 keep in mind that args can be register flonums, and make sure to insert
2949 SCM_FLONUM_ENSURE_MEM wherever they can leak out.
2950
2951 Caveat: Scm_NumCmp returns 0 (means equal) when arg0 and/or arg1 is/are NaN.
2952 That's because NaN doesn't make any sense in three-way comparison. The
2953 premise is that NaN has already been filtered out before NumCmp is called.
2954 */
Scm_NumCmp(ScmObj arg0,ScmObj arg1)2955 int Scm_NumCmp(ScmObj arg0, ScmObj arg1)
2956 {
2957 ScmObj badnum;
2958
2959 if (SCM_INTP(arg0)) {
2960 if (SCM_INTP(arg1)) {
2961 long r = SCM_INT_VALUE(arg0) - SCM_INT_VALUE(arg1);
2962 if (r < 0) return -1;
2963 if (r > 0) return 1;
2964 return 0;
2965 }
2966 if (SCM_FLONUMP(arg1)) {
2967 double r = SCM_INT_VALUE(arg0) - SCM_FLONUM_VALUE(arg1);
2968 if (r < 0) return -1;
2969 if (r > 0) return 1;
2970 return 0;
2971 }
2972 if (SCM_BIGNUMP(arg1))
2973 return Scm_BignumCmp(SCM_BIGNUM(Scm_MakeBignumFromSI(SCM_INT_VALUE(arg0))),
2974 SCM_BIGNUM(arg1));
2975 if (SCM_RATNUMP(arg1)) {
2976 if (SCM_EQ(SCM_INT_VALUE(0), arg0)) {
2977 return -Scm_Sign(arg1);
2978 } else {
2979 /* Roughly estimates the result by coercing the RATNUM to
2980 double. We have 53bits of precision (denomalized
2981 numbers only matter when ARG0 is zero, which is excluded
2982 already). */
2983 double y = Scm_GetDouble(arg1);
2984 double r = SCM_INT_VALUE(arg0) - y;
2985 double err = y * 2.0e-52;
2986
2987 if (r < -err) return -1;
2988 if (r > err) return 1;
2989 /* We need more precise comparison. */
2990 return Scm_NumCmp(Scm_Mul(arg0, SCM_RATNUM_DENOM(arg1)),
2991 SCM_RATNUM_NUMER(arg1));
2992 }
2993 }
2994 badnum = arg1;
2995 }
2996 else if (SCM_FLONUMP(arg0)) {
2997 if (SCM_INTP(arg1)) {
2998 double r = SCM_FLONUM_VALUE(arg0) - SCM_INT_VALUE(arg1);
2999 if (r < 0) return -1;
3000 if (r > 0) return 1;
3001 return 0;
3002 }
3003 if (SCM_FLONUMP(arg1)) {
3004 double r = SCM_FLONUM_VALUE(arg0) - SCM_FLONUM_VALUE(arg1);
3005 if (r < 0) return -1;
3006 if (r > 0) return 1;
3007 return 0;
3008 }
3009 if (SCM_BIGNUMP(arg1) || SCM_RATNUMP(arg1)) {
3010 /* NaN is already excluded. We filter out obvious. */
3011 if (SCM_IS_INF(SCM_FLONUM_VALUE(arg0))) {
3012 return Scm_Sign(arg0);
3013 }
3014 /* R7RS requires transitivity in '=', so we can't coerce
3015 arg1 to inexact; we'd rather convert arg0 to exact.
3016 NB: Thus we'll have a case that (= x y) => #f but
3017 (- x y) => 0.0. It's ok for the latter is inexact result. */
3018 return Scm_NumCmp(Scm_Exact(arg0), arg1);
3019 }
3020 badnum = arg1;
3021 }
3022 else if (SCM_BIGNUMP(arg0)) {
3023 if (SCM_INTP(arg1))
3024 return Scm_BignumCmp(SCM_BIGNUM(arg0),
3025 SCM_BIGNUM(Scm_MakeBignumFromSI(SCM_INT_VALUE(arg1))));
3026 if (SCM_FLONUMP(arg1))
3027 return -Scm_NumCmp(arg1, arg0);
3028 if (SCM_BIGNUMP(arg1))
3029 return Scm_BignumCmp(SCM_BIGNUM(arg0), SCM_BIGNUM(arg1));
3030 if (SCM_RATNUMP(arg1)) {
3031 /* we can't coerce to flonum, for it may lose precision. */
3032 ScmObj d1 = SCM_RATNUM_DENOM(arg1);
3033 return Scm_NumCmp(Scm_Mul(arg0, d1),
3034 SCM_RATNUM_NUMER(arg1));
3035 }
3036 badnum = arg1;
3037 }
3038 else if (SCM_RATNUMP(arg0)) {
3039 if (SCM_INTP(arg1) || SCM_BIGNUMP(arg1) || SCM_FLONUMP(arg1)) {
3040 return -Scm_NumCmp(arg1, arg0);
3041 }
3042 if (SCM_RATNUMP(arg1)) {
3043 ScmObj n0 = SCM_RATNUM_NUMER(arg0), d0 = SCM_RATNUM_DENOM(arg0);
3044 ScmObj n1 = SCM_RATNUM_NUMER(arg1), d1 = SCM_RATNUM_DENOM(arg1);
3045 int s0 = Scm_Sign(n0), s1 = Scm_Sign(n1);
3046
3047 /* screen the obvious cases without allocating new numbers */
3048 if (s0 < s1) return -1;
3049 if (s0 > s1) return 1;
3050 int d = Scm_NumCmp(d0, d1);
3051 if (d == 0) return Scm_NumCmp(n0, n1);
3052 if ((s0 > 0 && s1 > 0) || (s0 < 0 && s1 < 0)) {
3053 int n = Scm_NumCmp(n0, n1) * s0;
3054 if (d > 0 && n <= 0) return -s0;
3055 if (d < 0 && n >= 0) return s0;
3056 }
3057
3058 return Scm_NumCmp(Scm_Mul(n0, d1),
3059 Scm_Mul(n1, d0));
3060 }
3061 badnum = arg1;
3062 }
3063 else badnum = arg0;
3064
3065 SCM_FLONUM_ENSURE_MEM(badnum);
3066 Scm_Error("real number required: %S", badnum);
3067 return 0; /* dummy */
3068 }
3069
Scm_MinMax(ScmObj arg0,ScmObj args,ScmObj * min,ScmObj * max)3070 void Scm_MinMax(ScmObj arg0, ScmObj args, ScmObj *min, ScmObj *max)
3071 {
3072 if (!SCM_REALP(arg0)) Scm_Error("real number required, but got %S", arg0);
3073 if (nan_p(arg0)) goto got_nan;
3074 int inexact = !SCM_EXACTP(arg0);
3075 ScmObj mi = arg0;
3076 ScmObj ma = arg0;
3077
3078 for (;SCM_PAIRP(args); args = SCM_CDR(args)) {
3079 ScmObj arg = SCM_CAR(args);
3080 if (!SCM_REALP(arg)) Scm_Error("real number required, but got %S", arg);
3081 if (nan_p(arg)) goto got_nan;
3082 if (!SCM_EXACTP(arg)) inexact = TRUE;
3083 if (min && Scm_NumCmp(mi, arg) > 0) mi = arg;
3084 if (max && Scm_NumCmp(ma, arg) < 0) ma = arg;
3085 }
3086 if (min) {
3087 if (inexact && SCM_EXACTP(mi)) {
3088 *min = Scm_Inexact(mi);
3089 } else {
3090 *min = mi;
3091 }
3092 }
3093 if (max) {
3094 if (inexact && SCM_EXACTP(ma)) {
3095 *max = Scm_Inexact(ma);
3096 } else {
3097 *max = ma;
3098 }
3099 }
3100 return;
3101 got_nan:
3102 if (min) *min = SCM_NAN;
3103 if (max) *max = SCM_NAN;
3104 }
3105
3106 /*===============================================================
3107 * ROUNDING
3108 */
3109
3110 /* NB: rint() is not in POSIX, so an alternative is provided here.
3111 We don't use round(), for it behaves differently when the
3112 argument is exactly the halfway of two whole numbers. */
3113 #ifdef HAVE_RINT
3114 #define roundeven rint
3115 #else /* !HAVE_RINT */
roundeven(double v)3116 static double roundeven(double v)
3117 {
3118 double r;
3119 double frac = modf(v, &r);
3120 if (v > 0.0) {
3121 if (frac > 0.5) r += 1.0;
3122 else if (frac == 0.5) {
3123 if (fmod(r,2.0) != 0.0) r += 1.0;
3124 }
3125 } else {
3126 if (frac < -0.5) r -= 1.0;
3127 else if (frac == -0.5) {
3128 if (fmod(r,2.0) != 0.0) r -= 1.0;
3129 }
3130 }
3131 return r;
3132 }
3133 #endif /* !HAVE_RINT */
3134
scm_round(ScmObj num,int mode,int vmp)3135 static ScmObj scm_round(ScmObj num, int mode, int vmp)
3136 {
3137
3138 if (SCM_INTEGERP(num)) return num;
3139 if (SCM_RATNUMP(num)) {
3140 int offset = 0;
3141 ScmObj rem;
3142 ScmObj quot = Scm_Quotient(SCM_RATNUM_NUMER(num),
3143 SCM_RATNUM_DENOM(num), &rem);
3144 /* this shouldn't happen, but just in case.. */
3145 if (SCM_EXACT_ZERO_P(rem)) return quot;
3146
3147 /* Here we have quotient, which is always closer to zero
3148 than the original value */
3149 switch (mode) {
3150 case SCM_ROUND_FLOOR:
3151 offset = (Scm_Sign(num) < 0)? -1 : 0;
3152 break;
3153 case SCM_ROUND_CEIL:
3154 offset = (Scm_Sign(num) < 0)? 0 : 1;
3155 break;
3156 case SCM_ROUND_TRUNC:
3157 offset = 0;
3158 break;
3159 case SCM_ROUND_ROUND: {
3160 ScmObj rem2 = Scm_Mul(Scm_Abs(rem), SCM_MAKE_INT(2));
3161 int cmp = Scm_NumCmp(SCM_RATNUM_DENOM(num), rem2);
3162
3163 if (cmp > 0) {
3164 /* NUM is closer to zero than halfway */
3165 offset = 0;
3166 } else if (cmp < 0) {
3167 /* NUM is further from zero than halfway */
3168 offset = (Scm_Sign(num) < 0)? -1 : 1;
3169 } else {
3170 /* NUM is exactly the halfway. We round to even */
3171 if (Scm_OddP(quot)) {
3172 offset = (Scm_Sign(num) < 0)? -1 : 1;
3173 } else {
3174 offset = 0;
3175 }
3176 }
3177 break;
3178 }
3179 default: Scm_Panic("something screwed up");
3180 }
3181
3182 if (offset == 0) return quot;
3183 else return scm_add(quot, SCM_MAKE_INT(offset), vmp);
3184 }
3185 if (SCM_FLONUMP(num)) {
3186 double r = 0.0;
3187 double v = SCM_FLONUM_VALUE(num);
3188 switch (mode) {
3189 case SCM_ROUND_FLOOR: r = floor(v); break;
3190 case SCM_ROUND_CEIL: r = ceil(v); break;
3191 /* trunc is neither in ANSI nor in POSIX. */
3192 #ifdef HAVE_TRUNC
3193 case SCM_ROUND_TRUNC: r = trunc(v); break;
3194 #else
3195 case SCM_ROUND_TRUNC: r = (v < 0.0)? ceil(v) : floor(v); break;
3196 #endif
3197 case SCM_ROUND_ROUND: r = roundeven(v); break;
3198 default: Scm_Panic("something screwed up");
3199 }
3200 RETURN_FLONUM(r);
3201 }
3202 Scm_Error("real number required, but got %S", num);
3203 return SCM_UNDEFINED; /* dummy */
3204 }
3205
Scm_Round(ScmObj num,int mode)3206 ScmObj Scm_Round(ScmObj num, int mode)
3207 {
3208 return scm_round(num, mode, FALSE);
3209 }
3210
Scm_VMRound(ScmObj num,int mode)3211 ScmObj Scm_VMRound(ScmObj num, int mode)
3212 {
3213 return scm_round(num, mode, TRUE);
3214 }
3215
3216
Scm_RoundToExact(ScmObj num,int mode)3217 ScmObj Scm_RoundToExact(ScmObj num, int mode)
3218 {
3219 if (SCM_FLONUMP(num)) {
3220 double r = 0.0;
3221 double v = SCM_FLONUM_VALUE(num);
3222 if (SCM_IS_NAN(v) || SCM_IS_INF(v)) {
3223 Scm_Error("Exact infinity/nan is not supported: %S", num);
3224 }
3225 switch (mode) {
3226 case SCM_ROUND_FLOOR: r = floor(v); break;
3227 case SCM_ROUND_CEIL: r = ceil(v); break;
3228 /* trunc is neither in ANSI nor in POSIX. */
3229 #ifdef HAVE_TRUNC
3230 case SCM_ROUND_TRUNC: r = trunc(v); break;
3231 #else
3232 case SCM_ROUND_TRUNC: r = (v < 0.0)? ceil(v) : floor(v); break;
3233 #endif
3234 case SCM_ROUND_ROUND: r = roundeven(v); break;
3235 default: Scm_Panic("something screwed up");
3236 }
3237 return Scm_FlonumIntegerToExact(r);
3238 }
3239 if (SCM_INTEGERP(num)) return num;
3240 if (SCM_RATNUMP(num)) return Scm_Round(num, mode);
3241 Scm_Error("real number required, but got %S", num);
3242 return SCM_UNDEFINED; /* dummy */
3243 }
3244
3245 /*===============================================================
3246 * Logical (bitwise) operations
3247 */
3248
Scm_Ash(ScmObj x,ScmSmallInt cnt)3249 ScmObj Scm_Ash(ScmObj x, ScmSmallInt cnt)
3250 {
3251 /* TODO: This is an arbitrary limit, but we need *some* limit anyway
3252 to prevent a silly mistake from consuming large amount of memory.
3253 Eventually we need a consistent limit on how big a bignum can be. */
3254 if (cnt >= 0x10000000) {
3255 Scm_Error("ash: shift amount too big to handle: %ld", cnt);
3256 }
3257
3258 if (SCM_INTP(x)) {
3259 long ix = SCM_INT_VALUE(x);
3260 if (cnt <= -(SIZEOF_LONG * 8)) {
3261 ix = (ix < 0)? -1 : 0;
3262 return Scm_MakeInteger(ix);
3263 } else if (cnt < 0) {
3264 if (ix < 0) {
3265 ix = ~((~ix) >> (-cnt));
3266 } else {
3267 ix >>= -cnt;
3268 }
3269 return Scm_MakeInteger(ix);
3270 } else if (cnt < SCM_SMALL_INT_SIZE) {
3271 if (ix < 0) {
3272 if (-ix < (SCM_SMALL_INT_MAX >> cnt)) {
3273 ix <<= cnt;
3274 return Scm_MakeInteger(ix);
3275 }
3276 } else {
3277 if (ix < (SCM_SMALL_INT_MAX >> cnt)) {
3278 ix <<= cnt;
3279 return Scm_MakeInteger(ix);
3280 }
3281 }
3282 }
3283 /* Here, we know the result must be a bignum. */
3284 {
3285 ScmObj big = Scm_MakeBignumFromSI(ix);
3286 return Scm_BignumAsh(SCM_BIGNUM(big), cnt);
3287 }
3288 } else if (SCM_BIGNUMP(x)) {
3289 return Scm_BignumAsh(SCM_BIGNUM(x), cnt);
3290 }
3291 Scm_Error("exact integer required, but got %S", x);
3292 return SCM_UNDEFINED;
3293 }
3294
Scm_LogNot(ScmObj x)3295 ScmObj Scm_LogNot(ScmObj x)
3296 {
3297 if (!SCM_INTEGERP(x)) Scm_Error("exact integer required, but got %S", x);
3298 if (SCM_INTP(x)) {
3299 /* this won't cause an overflow */
3300 return SCM_MAKE_INT(~SCM_INT_VALUE(x));
3301 } else {
3302 return Scm_Negate(Scm_BignumAddSI(SCM_BIGNUM(x), 1));
3303 }
3304 }
3305
Scm_LogAnd(ScmObj x,ScmObj y)3306 ScmObj Scm_LogAnd(ScmObj x, ScmObj y)
3307 {
3308 if (!SCM_INTEGERP(x)) Scm_Error("exact integer required, but got %S", x);
3309 if (!SCM_INTEGERP(y)) Scm_Error("exact integer required, but got %S", y);
3310 if (SCM_INTP(x)) {
3311 if (SCM_INTP(y)) {
3312 return SCM_MAKE_INT(SCM_INT_VALUE(x) & SCM_INT_VALUE(y));
3313 } else if (SCM_INT_VALUE(x) >= 0 && SCM_BIGNUM_SIGN(y) >= 0) {
3314 return Scm_MakeInteger(SCM_INT_VALUE(x)&SCM_BIGNUM(y)->values[0]);
3315 }
3316 x = Scm_MakeBignumFromSI(SCM_INT_VALUE(x));
3317 } else if (SCM_INTP(y)) {
3318 if (SCM_INT_VALUE(y) >= 0 && SCM_BIGNUM_SIGN(x) >= 0) {
3319 return Scm_MakeInteger(SCM_INT_VALUE(y)&SCM_BIGNUM(x)->values[0]);
3320 }
3321 y = Scm_MakeBignumFromSI(SCM_INT_VALUE(y));
3322 }
3323 return Scm_BignumLogAnd(SCM_BIGNUM(x), SCM_BIGNUM(y));
3324 }
3325
Scm_LogIor(ScmObj x,ScmObj y)3326 ScmObj Scm_LogIor(ScmObj x, ScmObj y)
3327 {
3328 if (!SCM_INTEGERP(x)) Scm_Error("exact integer required, but got %S", x);
3329 if (!SCM_INTEGERP(y)) Scm_Error("exact integer required, but got %S", y);
3330 if (SCM_INTP(x)) {
3331 if (SCM_INTP(y))
3332 return SCM_MAKE_INT(SCM_INT_VALUE(x) | SCM_INT_VALUE(y));
3333 else
3334 x = Scm_MakeBignumFromSI(SCM_INT_VALUE(x));
3335 } else {
3336 if (SCM_INTP(y)) y = Scm_MakeBignumFromSI(SCM_INT_VALUE(y));
3337 }
3338 return Scm_BignumLogIor(SCM_BIGNUM(x), SCM_BIGNUM(y));
3339 }
3340
3341
Scm_LogXor(ScmObj x,ScmObj y)3342 ScmObj Scm_LogXor(ScmObj x, ScmObj y)
3343 {
3344 if (!SCM_INTEGERP(x)) Scm_Error("exact integer required, but got %S", x);
3345 if (!SCM_INTEGERP(y)) Scm_Error("exact integer required, but got %S", y);
3346 if (SCM_INTP(x)) {
3347 if (SCM_INTP(y))
3348 return SCM_MAKE_INT(SCM_INT_VALUE(x) ^ SCM_INT_VALUE(y));
3349 else
3350 x = Scm_MakeBignumFromSI(SCM_INT_VALUE(x));
3351 } else {
3352 if (SCM_INTP(y)) y = Scm_MakeBignumFromSI(SCM_INT_VALUE(y));
3353 }
3354 return Scm_BignumLogXor(SCM_BIGNUM(x), SCM_BIGNUM(y));
3355 }
3356
3357 /*===============================================================
3358 * Number I/O
3359 */
3360
3361 /* contants frequently used in number I/O */
3362 static double dexpt2_minus_52 = 0.0; /* 2.0^-52 */
3363 static double dexpt2_minus_53 = 0.0; /* 2.0^-53 */
3364
3365 /* max N where 10.0^N can be representable exactly in double.
3366 it is max N where N * log2(5) < 53. */
3367 #define MAX_EXACT_10_EXP 23
3368
3369 /* fast 10^n for limited cases */
iexpt10(int e)3370 static inline ScmObj iexpt10(int e)
3371 {
3372 if (e < IEXPT10_TABLESIZ) {
3373 return iexpt10_n[e];
3374 } else {
3375 /* This recursive case can happen if excessive number of decimal
3376 digits are given below the decimal point, and the exponent
3377 is close to minimum. There may be more efficient way to
3378 prune lower digits so that we can keep scale factor bounded,
3379 since those lower digits are very unlikely to contribute to
3380 the result. But it is difficult to avoid double-rounding
3381 error compeletely, so we take safer path for now. */
3382 return Scm_Mul(iexpt10_n[IEXPT10_TABLESIZ-1],
3383 iexpt10(e - IEXPT10_TABLESIZ + 1));
3384 }
3385 }
3386
3387 /* integer power of R by N, N is rather small.
3388 Assuming everything is in range. */
ipow(int r,int n)3389 static inline u_long ipow(int r, int n)
3390 {
3391 u_long k;
3392 for (k=1; n>0; n--) k *= r;
3393 return k;
3394 }
3395
3396 /* X * 10.0^N by double.
3397 10.0^N can be represented _exactly_ in double-precision floating point
3398 number in the range 0 <= N <= 23.
3399 If N is out of this range, a rounding error occurs, which will be
3400 corrected in the algorithmR routine below. */
raise_pow10(double x,int n)3401 static double raise_pow10(double x, int n)
3402 {
3403 static double dpow10[] = { 1.0, 1.0e1, 1.0e2, 1.0e3, 1.0e4,
3404 1.0e5, 1.0e6, 1.0e7, 1.0e8, 1.0e9,
3405 1.0e10, 1.0e11, 1.0e12, 1.0e13, 1.0e14,
3406 1.0e15, 1.0e16, 1.0e17, 1.0e18, 1.0e19,
3407 1.0e20, 1.0e21, 1.0e22, 1.0e23 };
3408 if (n >= 0) {
3409 while (n > 23) {
3410 x *= 1.0e24;
3411 n -= 24;
3412 }
3413 return x*dpow10[n];
3414 } else {
3415 while (n < -23) {
3416 x /= 1.0e24;
3417 n += 24;
3418 }
3419 return x/dpow10[-n];
3420 }
3421 }
3422
3423 /*
3424 * Number Printer
3425 *
3426 * This version implements Burger&Dybvig algorithm (Robert G. Burger
3427 * and and R. Kent Dybvig, "Priting Floating-Point Numbers Quickly and
3428 * Accurately", PLDI '96, pp.108--116, 1996).
3429 */
3430
3431 /* compare x+d and y. x, d, y are exact positive integers.
3432 this is called in inner loops so we need to be fast. */
numcmp3(ScmObj x,ScmObj d,ScmObj y)3433 static inline int numcmp3(ScmObj x, ScmObj d, ScmObj y)
3434 {
3435 if (SCM_INTP(x) && SCM_INTP(d) && SCM_INTP(y)) {
3436 long xd = SCM_INT_VALUE(x)+SCM_INT_VALUE(d);
3437 if (xd < SCM_INT_VALUE(y)) return -1;
3438 if (xd > SCM_INT_VALUE(y)) return 1;
3439 else return 0;
3440 } else {
3441 ScmObj bx = SCM_BIGNUMP(x)? x : Scm_MakeBignumFromSI(SCM_INT_VALUE(x));
3442 ScmObj bd = SCM_BIGNUMP(d)? d : Scm_MakeBignumFromSI(SCM_INT_VALUE(d));
3443 ScmObj by = SCM_BIGNUMP(y)? y : Scm_MakeBignumFromSI(SCM_INT_VALUE(y));
3444 return Scm_BignumCmp3U(SCM_BIGNUM(bx), SCM_BIGNUM(bd), SCM_BIGNUM(by));
3445 }
3446 }
3447
3448 /* Increment the given decimal number represented as a string.
3449
3450 /------------------ start
3451 | /---- end
3452 v v
3453 buf |- 3 . 1 4 1 5 9 |
3454
3455 In this example, we add 1 to the 10000th ('9') and handle any carryovers.
3456 If the carry over spills on the left, we shift entire string to right.
3457 Returns NUL-terminated string.
3458 */
notational_roundup(const char * src,int start,int end)3459 static char *notational_roundup(const char *src, int start, int end)
3460 {
3461 char *r = SCM_NEW_ATOMIC2(char*, end+2);
3462 memcpy(r, src, end);
3463 r[end] = r[end+1] = '\0';
3464
3465 /* loop as long as we carry over */
3466 char *s = r+start;
3467 char *p = r+end-1;
3468 for (; s <= p; p--) {
3469 if (*p == '.') continue;
3470 if (*p == '9') { *p = '0'; continue; }
3471 else { (*p)++; break; }
3472 }
3473 if (p < s && *s == '0') {
3474 /* spill over */
3475 for (p = r+end-1; s <= p; p--) p[1] = p[0];
3476 *s = '1';
3477 }
3478 return r;
3479 }
3480
3481 /* Notational rounding
3482 Should we use commercial rounding (round half towards infinity)
3483 or banker's rounding? Since notational rounding is inherently
3484 taking convenience than accuracy, we go with commercial rounding.
3485 */
notational_rounding(ScmDString * ds,int numstart,int precision)3486 static void notational_rounding(ScmDString *ds, int numstart, int precision)
3487 {
3488 Scm_DStringWeld(ds);
3489 ScmSmallInt size;
3490 const char *cbuf = Scm_DStringPeek(ds, &size, NULL);
3491 const char *end = cbuf + size;
3492 const char *p = cbuf;
3493 while (p < end && *p != '.') p++;
3494 SCM_ASSERT(p < end);
3495 p++; /* p is on tenths */
3496 if (end-p > precision) {
3497 char c = *(p+precision);
3498 if (c < '5') {
3499 /* round down - we just truncate */
3500 Scm_DStringTruncate(ds, p+precision-cbuf);
3501 } else {
3502 /* round up */
3503 char *nbuf = notational_roundup(cbuf, numstart,
3504 p+precision-cbuf);
3505 Scm_DStringTruncate(ds, 0);
3506 Scm_DStringPutz(ds, nbuf, -1);
3507 }
3508 }
3509 }
3510
3511 /* Increment the last digit and handle all carryovers.
3512 If we always emit highest precision, Burger&Dybvig algorithm strictly
3513 emits digits forward, without backing up to handle carry over.
3514 However, if we round the number before it, there can be a case that
3515 the last digit gets rounded up and its carry-over propagates upward.
3516 In such a case, this is called. The last digit in DS is '9', and we
3517 know we need to round it up.
3518 */
spill_fixup(ScmDString * ds,int numstart)3519 static void spill_fixup(ScmDString *ds, int numstart)
3520 {
3521 ScmSmallInt size;
3522 const char *cbuf = Scm_DStringPeek(ds, &size, NULL);
3523 char *nbuf = notational_roundup(cbuf, numstart, size);
3524 Scm_DStringTruncate(ds, 0);
3525 Scm_DStringPutz(ds, nbuf, -1);
3526 }
3527
3528 /* The main routine to get string representation of double.
3529 Convert VAL to a string and store to BUF, which must have at least FLT_BUF
3530 bytes long.
3531 True for PLUS_SIGN forces adding '+' for nonnegative numbers.
3532 EXP_LO and EXP_HI control when to switch exponential representation.
3533 We use n.nnne+zz representation when zz can be smaller than or equal
3534 to EXP_LO, or greater than or equal to EXP_HI.
3535 PRECISION specifies the number of digits to be printed after
3536 the decimal point; -1 means no limit.
3537 NOTATIONAL flags changes how the rounding with precision is done; if it's
3538 false, we choose closest decimal to the actual number represented by VAL.
3539 If it's true, we first generate optimal decimal notation, then round.
3540 */
print_double(ScmDString * ds,double val,int plus_sign,int precision,int notational,int exp_lo,int exp_hi,int exp_width)3541 static void print_double(ScmDString *ds, double val, int plus_sign,
3542 int precision, int notational,
3543 int exp_lo, int exp_hi, int exp_width)
3544 {
3545 /* Handle a few special cases first. */
3546 if (val == 0.0) {
3547 if (Scm_FlonumSign(val) > 0) {
3548 if (plus_sign) SCM_DSTRING_PUTC(ds, '+');
3549 } else {
3550 SCM_DSTRING_PUTC(ds, '-');
3551 }
3552 if (precision < 0) Scm_DStringPutz(ds, "0.0", 3);
3553 else {
3554 Scm_DStringPutz(ds, "0.", 2);
3555 for (int i=0; i<precision; i++) {
3556 Scm_DStringPutc(ds, '0');
3557 }
3558 }
3559 return;
3560 } else if (SCM_IS_INF(val)) {
3561 if (val < 0.0) Scm_DStringPutz(ds, "-inf.0", 6);
3562 else Scm_DStringPutz(ds, "+inf.0", 6);
3563 return;
3564 } else if (SCM_IS_NAN(val)) {
3565 Scm_DStringPutz(ds, "+nan.0", 6);
3566 return;
3567 }
3568
3569 if (val < 0.0) SCM_DSTRING_PUTC(ds, '-');
3570 else if (plus_sign) SCM_DSTRING_PUTC(ds, '+');
3571
3572 int numstart = Scm_DStringSize(ds); /* remember this for notational rounding */
3573
3574 /* variable names follows Burger&Dybvig paper. mp, mm for m+, m-.
3575 note that m+ == m- for most cases, and m+ == 2*m- for the rest.
3576 so we calculate m+ from m- for each iteration, using the flag
3577 mp2 as m+ = mp? m- : 2*m-. */
3578 ScmObj r, s, mm;
3579 int exp, sign;
3580 int mp2 = FALSE, fixup = FALSE;
3581 int fracdigs = -1; /* Count digits below the decimal point.
3582 Initial value is -1. Once we emit the decimal
3583 point it becomes 0, then we start counting. */
3584
3585 IEXPT10_INIT();
3586 if (val < 0) val = -val;
3587
3588 /* initialize r, s, m+ and m- */
3589 ScmObj f = Scm_DecodeFlonum(val, &exp, &sign);
3590 int round = !Scm_OddP(f);
3591 if (exp >= 0) {
3592 ScmObj be = Scm_Ash(SCM_MAKE_INT(1), exp);
3593 if (Scm_NumCmp(f, SCM_2_52) != 0) {
3594 r = Scm_Ash(f, exp+1);
3595 s = SCM_MAKE_INT(2);
3596 mp2= FALSE;
3597 mm = be;
3598 } else {
3599 r = Scm_Ash(f, exp+2);
3600 s = SCM_MAKE_INT(4);
3601 mp2 = TRUE;
3602 mm = be;
3603 }
3604 } else {
3605 if (exp == -1023 || Scm_NumCmp(f, SCM_2_52) != 0) {
3606 r = Scm_Ash(f, 1);
3607 s = Scm_Ash(SCM_MAKE_INT(1), -exp+1);
3608 mp2 = FALSE;
3609 mm = SCM_MAKE_INT(1);
3610 } else {
3611 r = Scm_Ash(f, 2);
3612 s = Scm_Ash(SCM_MAKE_INT(1), -exp+2);
3613 mp2 = TRUE;
3614 mm = SCM_MAKE_INT(1);
3615 }
3616 }
3617
3618 /* estimate scale */
3619 int est = (int)ceil(log10(val) - 0.1);
3620 if (est >= 0) {
3621 s = Scm_Mul(s, iexpt10(est));
3622 } else {
3623 ScmObj scale = iexpt10(-est);
3624 r = Scm_Mul(r, scale);
3625 mm = Scm_Mul(mm, scale);
3626 }
3627
3628 /* fixup. avoid calculating m+ for obvious case. */
3629 if (Scm_NumCmp(r, s) >= 0) {
3630 fixup = TRUE;
3631 } else {
3632 ScmObj mp = (mp2? Scm_Ash(mm, 1) : mm);
3633 if (round) {
3634 fixup = (numcmp3(r, mp, s) >= 0);
3635 } else {
3636 fixup = (numcmp3(r, mp, s) > 0);
3637 }
3638 }
3639 if (fixup) {
3640 s = Scm_Mul(s, SCM_MAKE_INT(10));
3641 est++;
3642 }
3643
3644 /* Scm_Printf(SCM_CURERR, "est=%d, r=%S, s=%S, mp=%S, mm=%S\n",
3645 est, r, s, mp, mm); */
3646
3647 /* Determine position of decimal point. we avoid exponential
3648 notation if exponent is small, i.e. 0.9 and 30.0 instead of
3649 9.0e-1 and 3.0e1. */
3650 int point;
3651 if (est < exp_hi && est > exp_lo) { point = est; est = 1; }
3652 else { point = 1; }
3653
3654 /* Now, we print XX.YYeZZ, where XX.YY is VAL*10^EST and
3655 ZZ is EST. If EST == 1 we omit exponent part. POINT is
3656 the number of digits in XX part; so POINT=1 for 1.23,
3657 POINT=2 for 12.3 and so on. */
3658
3659 /* If POINT <= 0, we need to emit preceding zeros. */
3660 if (point <= 0) {
3661 for (int digs=point; digs<1 ;digs++) {
3662 /* check the number of digits to be printed */
3663 if (!notational && precision >= 0 && fracdigs >= precision-1) {
3664 if (digs == 0) {
3665 /* whether r/s is greater than 1/2 or not defines the
3666 first number of significant digits */
3667 ScmObj r2 = Scm_Ash(r, 1);
3668 if (Scm_NumCmp(r2, s) > 0) SCM_DSTRING_PUTC(ds, '1');
3669 else SCM_DSTRING_PUTC(ds, '0');
3670 } else {
3671 SCM_DSTRING_PUTC(ds, '0');
3672 }
3673 if (digs == point) SCM_DSTRING_PUTC(ds, '.');
3674 /* no more digits. */
3675 goto show_exponent;
3676 }
3677 SCM_DSTRING_PUTC(ds, '0');
3678 if (digs == point) SCM_DSTRING_PUTC(ds, '.');
3679 fracdigs++;
3680 }
3681 }
3682
3683 /* generate the digits */
3684 int spilled = FALSE;
3685 int digs;
3686 for (digs=1; ; digs++) {
3687 ScmObj r10 = Scm_Mul(r, SCM_MAKE_INT(10));
3688 ScmObj q = Scm_Quotient(r10, s, &r);
3689 ScmObj mp;
3690
3691 /* if we round early, we should extend the boundary
3692 _one digit before the rounded digit_. It's a bit complicated
3693 since when precision == 0 we have to do it before printing
3694 decimal point. */
3695 if (!notational
3696 && ((precision == 0 && digs == point)
3697 || (precision > 0 && fracdigs >= precision-1))) {
3698 mm = mp = Scm_Ash(s, -1);
3699 } else {
3700 mm = Scm_Mul(mm, SCM_MAKE_INT(10));
3701 mp = (mp2? Scm_Ash(mm, 1) : mm);
3702 }
3703
3704 /*Scm_Printf(SCM_CURERR, "q=%S, r=%S, s=%S mp=%S, mm=%S\n",
3705 q, r, s, mp, mm);*/
3706
3707 SCM_ASSERT(SCM_INTP(q));
3708 int tc1, tc2;
3709 if (round) {
3710 tc1 = (Scm_NumCmp(r, mm) <= 0);
3711 tc2 = (numcmp3(r, mp, s) >= 0);
3712 } else {
3713 tc1 = (Scm_NumCmp(r, mm) < 0);
3714 tc2 = (numcmp3(r, mp, s) > 0);
3715 }
3716 if (!tc1) {
3717 if (!tc2) {
3718 ScmChar c = (char)SCM_INT_VALUE(q) + '0';
3719 SCM_DSTRING_PUTC(ds, c);
3720 if (digs == point) SCM_DSTRING_PUTC(ds, '.');
3721 if (digs >= point) fracdigs++;
3722 continue;
3723 } else {
3724 ScmChar c = (char)SCM_INT_VALUE(q) + '1';
3725 if (c > '9') { spilled = TRUE; c = '9'; }
3726 SCM_DSTRING_PUTC(ds, c);
3727 if (digs > point) fracdigs++;
3728 break;
3729 }
3730 } else {
3731 if (!tc2) {
3732 SCM_ASSERT(SCM_INTP(q));
3733 ScmChar c = (char)SCM_INT_VALUE(q) + '0';
3734 SCM_DSTRING_PUTC(ds, c);
3735 if (digs > point) fracdigs++;
3736 break;
3737 } else {
3738 int tc3 = numcmp3(r, r, s); /* r*2 <=> s */
3739 if ((round && tc3 <= 0) || (!round && tc3 < 0)) {
3740 ScmChar c = (char)SCM_INT_VALUE(q) + '0';
3741 SCM_DSTRING_PUTC(ds, c);
3742 if (digs > point) fracdigs++;
3743 break;
3744 } else {
3745 ScmChar c = (char)SCM_INT_VALUE(q) + '1';
3746 if (c > '9') { spilled = TRUE; c = '9'; }
3747 SCM_DSTRING_PUTC(ds, c);
3748 if (digs > point) fracdigs++;
3749 break;
3750 }
3751 }
3752 }
3753 }
3754 if (spilled) spill_fixup(ds, numstart);
3755
3756 /* Notational rounding, if necessary */
3757 if (notational && precision >= 0 && fracdigs > precision) {
3758 notational_rounding(ds, numstart, precision);
3759 }
3760
3761 /* print the trailing zeros if necessary */
3762 if (digs <= point) {
3763 for (;digs < point; digs++) SCM_DSTRING_PUTC(ds, '0');
3764 SCM_DSTRING_PUTC(ds, '.');
3765 if (precision < 0) SCM_DSTRING_PUTC(ds, '0');
3766 }
3767 for (;(digs-point) < precision; digs++) {
3768 SCM_DSTRING_PUTC(ds, '0');
3769 }
3770
3771 show_exponent:
3772 SCM_ASSERT(est < 1000 && est > -1000);
3773 /* prints exponent. we shifted decimal point, so -1. */
3774 est--;
3775 if (est != 0) {
3776 SCM_DSTRING_PUTC(ds, 'e');
3777 if (est < 0) {
3778 Scm_DStringPutc(ds, '-');
3779 est = -est;
3780 }
3781 char zbuf[5]; /* we know est is at most 4 digits */
3782 int echars = sprintf(zbuf, "%d", (int)est);
3783 if (echars < exp_width) {
3784 int fill = exp_width - echars;
3785 while (fill--) {
3786 Scm_DStringPutc(ds, '0');
3787 }
3788 }
3789 Scm_DStringPutz(ds, zbuf, -1);
3790 }
3791 }
3792
number_print(ScmObj obj,ScmPort * port,ScmWriteContext * ctx SCM_UNUSED)3793 static void number_print(ScmObj obj, ScmPort *port,
3794 ScmWriteContext *ctx SCM_UNUSED)
3795 {
3796 Scm_PrintNumber(port, obj, NULL);
3797 }
3798
3799 #define FLT_BUF 65 /* need to hold binary representation of the least fixnum */
3800
3801 static size_t
print_radix_prefix(ScmPort * port,u_long radix)3802 print_radix_prefix(ScmPort *port, u_long radix)
3803 {
3804 char buf[FLT_BUF];
3805 switch (radix) {
3806 case 2: Scm_Putz("#b", 2, port); return 2;
3807 case 8: Scm_Putz("#o", 2, port); return 2;
3808 case 10: Scm_Putz("#d", 2, port); return 2;
3809 case 16: Scm_Putz("#x", 2, port); return 2;
3810 default:
3811 {
3812 int nc = snprintf(buf, sizeof(buf), "#%lur", radix);
3813 Scm_Putz(buf, nc, port);
3814 return nc;
3815 }
3816 }
3817 }
3818
3819 static size_t
print_number(ScmPort * port,ScmObj obj,u_long flags,ScmNumberFormat * fmt)3820 print_number(ScmPort *port, ScmObj obj, u_long flags, ScmNumberFormat *fmt)
3821 {
3822 int use_upper = flags & SCM_NUMBER_FORMAT_USE_UPPER;
3823 int show_plus = flags & SCM_NUMBER_FORMAT_SHOW_PLUS;
3824 int radix = fmt->radix;
3825 int nchars = 0;
3826 char buf[FLT_BUF];
3827
3828 if ((flags & SCM_NUMBER_FORMAT_ALT_RADIX) && SCM_EXACTP(obj)) {
3829 nchars += print_radix_prefix(port, radix);
3830 }
3831
3832 if (SCM_INTP(obj)) {
3833 long value = SCM_INT_VALUE(obj);
3834 if (value < 0) {
3835 SCM_PUTC('-', port);
3836 nchars++;
3837 value = -value; /* this won't overflow */
3838 } else if (show_plus) {
3839 SCM_PUTC('+', port);
3840 nchars++;
3841 }
3842 if (value == 0) { SCM_PUTC('0', port); return nchars+1; }
3843 int i;
3844 for (i = FLT_BUF-1; i >= 0 && value > 0; i--) {
3845 int c = value % radix;
3846 buf[i] = (c<10)?(c+'0'):(use_upper?(c-10+'A'):(c-10+'a'));
3847 value /= radix;
3848 nchars++;
3849 }
3850 Scm_Putz(buf+i+1, FLT_BUF-i-1, port);
3851 return nchars;
3852 } else if (SCM_BIGNUMP(obj)) {
3853 ScmObj s = Scm_BignumToString(SCM_BIGNUM(obj), radix, use_upper);
3854 if(show_plus && Scm_Sign(obj) >= 0) {
3855 Scm_Putc('+', port);
3856 nchars++;
3857 }
3858 Scm_Puts(SCM_STRING(s), port);
3859 return nchars + SCM_STRING_BODY_LENGTH(SCM_STRING_BODY(s));
3860 } else if (SCM_FLONUMP(obj)) {
3861 ScmDString ds;
3862 Scm_DStringInit(&ds);
3863 print_double(&ds, SCM_FLONUM_VALUE(obj), show_plus,
3864 fmt->precision,
3865 fmt->flags&SCM_NUMBER_FORMAT_ROUND_NOTATIONAL,
3866 fmt->exp_lo, fmt->exp_hi, fmt->exp_width);
3867 Scm_Putz(Scm_DStringGetz(&ds), -1, port);
3868 return Scm_DStringSize(&ds);
3869 } else if (SCM_RATNUMP(obj)) {
3870 u_long flags2 = flags & ~SCM_NUMBER_FORMAT_ALT_RADIX;
3871 nchars = print_number(port, SCM_RATNUM_NUMER(obj), flags2, fmt);
3872 Scm_Putc('/', port);
3873 nchars++;
3874 flags2 &= ~SCM_NUMBER_FORMAT_SHOW_PLUS;
3875 nchars += print_number(port, SCM_RATNUM_DENOM(obj), flags2, fmt);
3876 return nchars;
3877 } else if (SCM_COMPNUMP(obj)) {
3878 ScmDString ds;
3879 Scm_DStringInit(&ds);
3880 print_double(&ds, SCM_COMPNUM_REAL(obj), show_plus,
3881 fmt->precision,
3882 fmt->flags&SCM_NUMBER_FORMAT_ROUND_NOTATIONAL,
3883 fmt->exp_lo, fmt->exp_hi, fmt->exp_width);
3884 Scm_Putz(Scm_DStringGetz(&ds), -1, port);
3885 nchars += Scm_DStringSize(&ds);
3886 Scm_DStringTruncate(&ds, 0);
3887 print_double(&ds, SCM_COMPNUM_IMAG(obj), TRUE,
3888 fmt->precision,
3889 fmt->flags&SCM_NUMBER_FORMAT_ROUND_NOTATIONAL,
3890 fmt->exp_lo, fmt->exp_hi, fmt->exp_width);
3891 Scm_Putz(Scm_DStringGetz(&ds), -1, port);
3892 nchars += Scm_DStringSize(&ds);
3893 Scm_Putc('i', port);
3894 return nchars+1;
3895 } else {
3896 Scm_Error("number required: %S", obj);
3897 return 0; /* dummy */
3898 }
3899 }
3900
3901 /* API */
Scm_NumberFormatInit(ScmNumberFormat * fmt)3902 void Scm_NumberFormatInit(ScmNumberFormat* fmt)
3903 {
3904 fmt->flags = 0;
3905 fmt->radix = 10;
3906 fmt->precision = -1;
3907 fmt->exp_lo = -3;
3908 fmt->exp_hi = 10;
3909 fmt->exp_width = 0;
3910 }
3911
3912 /* API */
Scm_NumberToString(ScmObj obj,int radix,u_long flags)3913 ScmObj Scm_NumberToString(ScmObj obj, int radix, u_long flags)
3914 {
3915 if (radix < SCM_RADIX_MIN || radix > SCM_RADIX_MAX)
3916 Scm_Error("radix out of range: %d", radix);
3917 ScmPort *p = SCM_PORT(Scm_MakeOutputStringPort(TRUE));
3918 ScmNumberFormat fmt;
3919 Scm_NumberFormatInit(&fmt);
3920 fmt.flags = flags;
3921 fmt.radix = radix;
3922 Scm_PrintNumber(p, obj, &fmt);
3923 return Scm_GetOutputString(p, 0);
3924 }
3925
3926 /* API. FMT can be NULL. */
Scm_PrintNumber(ScmPort * port,ScmObj n,ScmNumberFormat * fmt)3927 size_t Scm_PrintNumber(ScmPort *port, ScmObj n, ScmNumberFormat *fmt)
3928 {
3929 ScmNumberFormat defaults;
3930 if (fmt == NULL) {
3931 Scm_NumberFormatInit(&defaults);
3932 fmt = &defaults;
3933 }
3934 return print_number(port, n, fmt->flags, fmt);
3935 }
3936
3937 /* API. FMT can be NULL. Utility to expose Burger&Dybvig algorithm. */
Scm_PrintDouble(ScmPort * port,double d,ScmNumberFormat * fmt)3938 size_t Scm_PrintDouble(ScmPort *port, double d, ScmNumberFormat *fmt)
3939 {
3940 ScmNumberFormat defaults;
3941 if (fmt == NULL) {
3942 Scm_NumberFormatInit(&defaults);
3943 fmt = &defaults;
3944 }
3945 ScmDString ds;
3946 Scm_DStringInit(&ds);
3947 print_double(&ds, d,
3948 fmt->flags & SCM_NUMBER_FORMAT_SHOW_PLUS,
3949 fmt->precision,
3950 fmt->flags & SCM_NUMBER_FORMAT_ROUND_NOTATIONAL,
3951 fmt->exp_lo, fmt->exp_hi, fmt->exp_width);
3952 size_t nchars = Scm_DStringSize(&ds);
3953 Scm_Putz(Scm_DStringGetz(&ds), (int)nchars, port);
3954 return nchars;
3955 }
3956
3957 /*
3958 * Number Parser
3959 *
3960 * <number> : <prefix> <complex>
3961 * <prefix> : <radix> <exactness> | <exactness> <radix>
3962 * <radix> : <empty> | '#b' | '#o' | '#d' | '#x'
3963 * <exactness> : <empty> | '#e' | '#i'
3964 * <complex> : <real>
3965 * | <real> '@' <real> ['p' 'i']
3966 * | <real> '+' <ureal> 'i'
3967 * | <real> '-' <ureal> 'i'
3968 * | <real> '+' 'i'
3969 * | <real> '-' 'i'
3970 * | '+' <ureal> 'i'
3971 * | '-' <ureal> 'i'
3972 * | '+' 'i'
3973 * | '-' 'i'
3974 * <real> : <sign> <ureal>
3975 * <sign> : <empty> | '+' | '-'
3976 * <ureal> : <uinteger>
3977 * | <uinteger> '/' <uinteger>
3978 * | <decimal>
3979 * <uinteger> : <digit>+ '#'*
3980 * <decimal> : <digit10>+ '#'* <suffix>
3981 * | '.' <digit10>+ '#'* <suffix>
3982 * | <digit10>+ '.' <digit10>+ '#'* <suffix>
3983 * | <digit10>+ '#'+ '.' '#'* <suffix>
3984 * <suffix> : <empty> | <exponent-marker> <sign> <digit10>+
3985 * <exponent-marker> : 'e' | 's' | 'f' | 'd' | 'l'
3986 *
3987 * The parser reads characters from on-memory buffer.
3988 * Multibyte strings are filtered out in the early stage of
3989 * parsing, so the subroutines assume the buffer contains
3990 * only ASCII chars.
3991 */
3992
3993 struct numread_packet {
3994 const char *buffer; /* original buffer */
3995 int buflen; /* original length */
3996 int radix; /* radix (can be overridden by prefix) */
3997 int noradixprefix; /* do not allow radix prefix */
3998 int exactness; /* exactness; see enum below */
3999 int padread; /* '#' padding has been read */
4000 int explicit; /* explicit prefix is appeared */
4001 int strict; /* reject gauche extension */
4002 int throwerror; /* throws error on parse, instead of
4003 returning #f. */
4004 };
4005
4006 enum { /* used in the exactness flag */
4007 NOEXACT, EXACT, INEXACT
4008 };
4009
4010 /* Max digits D such that all D-digit radix R integers fit in signed
4011 long, i.e. R^(D+1)-1 <= LONG_MAX */
4012 static long longdigs[SCM_RADIX_MAX-SCM_RADIX_MIN+1] = { 0 };
4013
4014 /* Max integer I such that reading next digit (in radix R) will overflow
4015 long integer. floor(LONG_MAX/R - R). */
4016 static u_long longlimit[SCM_RADIX_MAX-SCM_RADIX_MIN+1] = { 0 };
4017
4018 /* An integer table of R^D, which is a "big digit" to be added
4019 into bignum. */
4020 static u_long bigdig[SCM_RADIX_MAX-SCM_RADIX_MIN+1] = { 0 };
4021
4022 static ScmObj numread_error(const char *msg, struct numread_packet *ctx);
4023
4024 /* Returns either small integer or bignum.
4025 initval may be a Scheme integer that will be 'concatenated' before
4026 the integer to be read; it is used to read floating-point number.
4027 Note that value_big may keep denormalized bignum. */
read_uint(const char ** strp,int * lenp,struct numread_packet * ctx,ScmObj initval)4028 static ScmObj read_uint(const char **strp, int *lenp,
4029 struct numread_packet *ctx,
4030 ScmObj initval)
4031 {
4032 const char *str = *strp;
4033 int digread = FALSE;
4034 int len = *lenp;
4035 int radix = ctx->radix;
4036 int digits = 0, diglimit = longdigs[radix-SCM_RADIX_MIN];
4037 u_long limit = longlimit[radix-SCM_RADIX_MIN], bdig = bigdig[radix-SCM_RADIX_MIN];
4038 u_long value_int = 0;
4039 ScmBignum *value_big = NULL;
4040 static const char tab[] = "0123456789abcdefghijklmnopqrstuvwxyz";
4041
4042 if (!SCM_FALSEP(initval)) {
4043 if (SCM_INTP(initval)) {
4044 if ((u_long)SCM_INT_VALUE(initval) > limit) {
4045 value_big = Scm_MakeBignumWithSize(4, SCM_INT_VALUE(initval));
4046 } else {
4047 value_int = SCM_INT_VALUE(initval);
4048 }
4049 } else if (SCM_BIGNUMP(initval)) {
4050 value_big = SCM_BIGNUM(Scm_BignumCopy(SCM_BIGNUM(initval)));
4051 }
4052 digread = TRUE;
4053 } else if (*str == '0') {
4054 /* Ignore leading 0's, to avoid unnecessary bignum operations. */
4055 while (len > 0 && *str == '0') { str++; len--; }
4056 digread = TRUE;
4057 }
4058
4059 while (len--) {
4060 int digval = -1;
4061 char c = tolower(*str++);
4062 if (ctx->explicit && !ctx->strict && c == '_') {
4063 /* Gauche extension - allow '_' in digits for readability
4064 when number is expliticly prefixed. */
4065 continue;
4066 }
4067 if (ctx->padread) {
4068 if (c == '#') digval = 0;
4069 else break;
4070 } else if (digread && c == '#') {
4071 digval = 0;
4072 ctx->padread = TRUE;
4073 if (ctx->exactness == NOEXACT) {
4074 ctx->exactness = INEXACT;
4075 }
4076 } else {
4077 for (const char *ptab = tab; ptab < tab+radix; ptab++) {
4078 if (c == *ptab) {
4079 digval = (int)(ptab-tab);
4080 digread = TRUE;
4081 break;
4082 }
4083 }
4084 }
4085 if (digval < 0) break;
4086 value_int = value_int * radix + digval;
4087 digits++;
4088 if (value_big == NULL) {
4089 if (value_int >= limit) {
4090 value_big = Scm_MakeBignumWithSize(4, value_int);
4091 value_int = digits = 0;
4092 }
4093 } else if (digits > diglimit) {
4094 value_big = Scm_BignumAccMultAddUI(value_big, bdig, value_int);
4095 value_int = digits = 0;
4096 }
4097 }
4098 *strp = str-1;
4099 *lenp = len+1;
4100
4101 if (value_big == NULL) return Scm_MakeInteger(value_int);
4102 if (digits > 0) {
4103 value_big = Scm_BignumAccMultAddUI(value_big,
4104 ipow(radix, digits),
4105 value_int);
4106 }
4107 return Scm_NormalizeBignum(SCM_BIGNUM(value_big));
4108 }
4109
4110 /*
4111 * Find a double number closest to f * 10^e, using z as the starting
4112 * approximation. The algorithm (and its name) is taken from Will Clinger's
4113 * paper "How to Read Floating Point Numbers Accurately", in the ACM
4114 * SIGPLAN '90, pp.92--101.
4115 * The algorithm is modified to take advantage of coherency between loops.
4116 */
algorithmR(ScmObj f,int e,double z)4117 static double algorithmR(ScmObj f, int e, double z)
4118 {
4119 int k, s;
4120 ScmObj m = Scm_DecodeFlonum(z, &k, &s);
4121 ScmObj x, y;
4122 IEXPT10_INIT();
4123 retry:
4124 if (k >= 0) {
4125 if (e >= 0) {
4126 x = Scm_Mul(f, iexpt10(e));
4127 y = Scm_Ash(m, k);
4128 } else {
4129 x = f;
4130 y = Scm_Ash(Scm_Mul(m, iexpt10(-e)), k);
4131 }
4132 } else {
4133 if (e >= 0) {
4134 x = Scm_Ash(Scm_Mul(f, iexpt10(e)), -k);
4135 y = m;
4136 } else {
4137 x = Scm_Ash(f, -k);
4138 y = Scm_Mul(m, iexpt10(-e));
4139 }
4140 }
4141 int kprev = k;
4142
4143 /* compare */
4144 for (;;) {
4145 /*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);*/
4146 /* compare */
4147 int sign_d = Scm_NumCmp(x, y);
4148 ScmObj abs_d = (sign_d > 0)? Scm_Sub(x, y) : Scm_Sub(y, x);
4149 ScmObj d2 = Scm_Ash(Scm_Mul(m, abs_d), 1);
4150 switch (Scm_NumCmp(d2, y)) {
4151 case -1: /* d2 < y */
4152 if (Scm_NumCmp(m, SCM_2_52) == 0
4153 && sign_d < 0
4154 && k > -1074
4155 && Scm_NumCmp(Scm_Ash(d2, 1), y) > 0) {
4156 goto prevfloat;
4157 } else {
4158 return ldexp(Scm_GetDouble(m), k);
4159 }
4160 case 0: /* d2 == y */
4161 if (!Scm_OddP(m)) {
4162 if (Scm_NumCmp(m, SCM_2_52) == 0
4163 && sign_d < 0) {
4164 goto prevfloat;
4165 } else {
4166 return ldexp(Scm_GetDouble(m), k);
4167 }
4168 } else if (sign_d < 0) {
4169 goto prevfloat;
4170 } else {
4171 goto nextfloat;
4172 }
4173 default:
4174 if (sign_d < 0) goto prevfloat;
4175 else goto nextfloat;
4176 }
4177 prevfloat:
4178 m = Scm_Sub(m, SCM_MAKE_INT(1));
4179 if (k > -1074 && Scm_NumCmp(m, SCM_2_52) < 0) {
4180 m = Scm_Ash(m, 1);
4181 k--;
4182 }
4183 goto next;
4184 nextfloat:
4185 m = Scm_Add(m, SCM_MAKE_INT(1));
4186 if (Scm_NumCmp(m, SCM_2_53) >= 0) {
4187 m = Scm_Ash(m, -1);
4188 k++;
4189 }
4190 /*FALLTHROUGH*/
4191 next:
4192 if (kprev >= 0) {
4193 if (k >= 0) {
4194 /* k stays positive. x is invariant */
4195 if (e >= 0) {
4196 y = Scm_Ash(m, k);
4197 } else {
4198 y = Scm_Ash(Scm_Mul(m, iexpt10(-e)), k);
4199 }
4200 } else {
4201 /* k turned to negative */
4202 goto retry;
4203 }
4204 } else {
4205 if (k < 0) {
4206 /* k stays negative. */
4207 if (e >= 0) {
4208 if (k != kprev) x = Scm_Ash(Scm_Mul(f, iexpt10(e)), -k);
4209 y = m;
4210 } else {
4211 if (k != kprev) x = Scm_Ash(f, -k);
4212 y = Scm_Mul(m, iexpt10(-e));
4213 }
4214 } else {
4215 /* k turned to positive */
4216 goto retry;
4217 }
4218 }
4219 }
4220 /*NOTREACHED*/
4221 }
4222
read_real(const char ** strp,int * lenp,struct numread_packet * ctx)4223 static ScmObj read_real(const char **strp, int *lenp,
4224 struct numread_packet *ctx)
4225 {
4226 int minusp = FALSE, exp_minusp = FALSE, exp_overflow = FALSE;
4227 int sign_seen = FALSE;
4228 int fracdigs = 0;
4229 long exponent = 0;
4230 ScmObj intpart, fraction;
4231 const char *mark; /* will point (*strp)[1] if there's a sign,
4232 otherwise (*strp)[0], to check if the
4233 given str has valid content. */
4234
4235 switch (**strp) {
4236 case '-': minusp = TRUE;
4237 /* FALLTHROUGH */
4238 case '+':
4239 (*strp)++; (*lenp)--; sign_seen = TRUE;
4240 }
4241 if ((*lenp) <= 0) return SCM_FALSE;
4242 mark = *strp;
4243
4244 /* Recognize specials */
4245 if (sign_seen && (*lenp) >= 5) {
4246 if (strncasecmp(*strp, "inf.0", 5) == 0) {
4247 (*strp) += 5; (*lenp) -= 5;
4248 return minusp?SCM_NEGATIVE_INFINITY:SCM_POSITIVE_INFINITY;
4249 }
4250 if (strncasecmp(*strp, "nan.0", 5) == 0) {
4251 (*strp) += 5; (*lenp) -= 5;
4252 return SCM_NAN;
4253 }
4254 }
4255
4256 /* Read integral part */
4257 if (**strp != '.') {
4258 intpart = read_uint(strp, lenp, ctx, SCM_FALSE);
4259 if ((*lenp) <= 0) {
4260 if (minusp) intpart = Scm_Negate(intpart);
4261 if (ctx->exactness == INEXACT) {
4262 return Scm_Inexact(intpart);
4263 } else {
4264 return intpart;
4265 }
4266 }
4267 if (**strp == '/') {
4268 /* possibly rational */
4269 ScmObj denom;
4270 int lensave;
4271
4272 if ((*lenp) <= 1 || mark == *strp) return SCM_FALSE;
4273 (*strp)++; (*lenp)--;
4274 lensave = *lenp;
4275 denom = read_uint(strp, lenp, ctx, SCM_FALSE);
4276 if (SCM_FALSEP(denom)) return SCM_FALSE;
4277 if (SCM_EXACT_ZERO_P(denom)) {
4278 if (lensave > *lenp) {
4279 if (ctx->exactness != INEXACT) {
4280 return numread_error("(exact infinity/nan is not supported.)",
4281 ctx);
4282 }
4283 if (SCM_EXACT_ZERO_P(intpart)) return SCM_NAN;
4284 return minusp? SCM_NEGATIVE_INFINITY:SCM_POSITIVE_INFINITY;
4285 } else {
4286 return SCM_FALSE;
4287 }
4288 }
4289 if (minusp) intpart = Scm_Negate(intpart);
4290 if (ctx->exactness == INEXACT) {
4291 return Scm_Inexact(Scm_Div(intpart, denom));
4292 } else {
4293 return Scm_MakeRational(intpart, denom);
4294 }
4295 }
4296 /* fallthrough */
4297 } else {
4298 intpart = SCM_FALSE; /* indicate there was no intpart */
4299 }
4300
4301 /* Read fractional part.
4302 At this point, simple integer is already eliminated. */
4303 if (**strp == '.') {
4304 if (ctx->radix != 10) {
4305 return numread_error("(only 10-based fraction is supported)", ctx);
4306 }
4307 (*strp)++; (*lenp)--;
4308 int lensave = *lenp;
4309 fraction = read_uint(strp, lenp, ctx, intpart);
4310 fracdigs = lensave - *lenp;
4311 } else {
4312 fraction = intpart;
4313 }
4314
4315 if (SCM_FALSEP(intpart)) {
4316 if (fracdigs == 0) return SCM_FALSE; /* input was "." */
4317 }
4318 if (mark == *strp) return SCM_FALSE;
4319
4320 /* Read exponent. */
4321 if (*lenp > 0 && strchr("eEsSfFdDlL", (int)**strp)) {
4322 (*strp)++;
4323 if (--(*lenp) <= 0) return SCM_FALSE;
4324 switch (**strp) {
4325 case '-': exp_minusp = TRUE;
4326 /*FALLTHROUGH*/
4327 case '+':
4328 (*strp)++;
4329 if (--(*lenp) <= 0) return SCM_FALSE;
4330 }
4331 while (*lenp > 0) {
4332 int c = **strp;
4333 if (!isdigit(c)) break;
4334 (*strp)++, (*lenp)--;
4335 if (isdigit(c) && !exp_overflow) {
4336 exponent = exponent * 10 + (c - '0');
4337 /* Check obviously wrong exponent range. More subtle check
4338 will be done later. */
4339 if (exponent >= MAX_EXPONENT) {
4340 exp_overflow = TRUE;
4341 }
4342 }
4343 }
4344 if (exp_minusp) exponent = -exponent;
4345 }
4346 if (exp_overflow) {
4347 if (ctx->exactness == EXACT) {
4348 /* Although we can represent such a number using bignum and
4349 ratnum, such large (or small) exponent is highly unusual
4350 and we assume we can report implementation limitation
4351 violation. */
4352 return numread_error("(such an exact number is out of implementation limitation)",
4353 ctx);
4354 }
4355 if (exp_minusp || SCM_EQ(fraction, SCM_MAKE_INT(0))) {
4356 return Scm_MakeFlonum(minusp? -0.0:0.0);
4357 } else {
4358 return minusp? SCM_NEGATIVE_INFINITY : SCM_POSITIVE_INFINITY;
4359 }
4360 }
4361
4362 /*Scm_Printf(SCM_CURERR, "fraction=%S, exponent=%d\n", fraction, exponent);*/
4363
4364 /* Compose the number. */
4365 if (ctx->exactness == EXACT) {
4366 /* Explicit exact number. We can continue exact arithmetic
4367 (it may end up ratnum) */
4368 ScmObj e = Scm_Mul(fraction,
4369 Scm_ExactIntegerExpt(SCM_MAKE_INT(10),
4370 Scm_MakeInteger(exponent-fracdigs)));
4371 if (minusp) return Scm_Negate(e);
4372 else return e;
4373 }
4374
4375 /* Get double approximaiton of fraction. If fraction >= 2^53 we'll
4376 only get approximation, but the error will be corrected in
4377 AlgorithmR. We have to be careful, however, not to overflow
4378 the following GetDouble call. */
4379 int raise_factor = exponent - fracdigs;
4380 double realnum = Scm_GetDouble(fraction);
4381
4382 if (SCM_IS_INF(realnum)) {
4383 /* We have too many digits to fit in double. We can still get finite
4384 number if raise_factor is small, but we need to calculate realnum
4385 via rational. */
4386 if (raise_factor >= 0) {
4387 /* Exponent is too big. */
4388 return (minusp? SCM_NEGATIVE_INFINITY : SCM_POSITIVE_INFINITY);
4389 }
4390 IEXPT10_INIT();
4391 realnum = Scm_GetDouble(Scm_Div(fraction, iexpt10(-raise_factor)));
4392 } else {
4393 realnum = raise_pow10(realnum, raise_factor);
4394 }
4395
4396 if (SCM_IS_INF(realnum)) {
4397 /* Exponent is too big. */
4398 return (minusp? SCM_NEGATIVE_INFINITY : SCM_POSITIVE_INFINITY);
4399 }
4400
4401 if (realnum > 0.0
4402 && (Scm_NumCmp(fraction, SCM_2_52) > 0
4403 || raise_factor > MAX_EXACT_10_EXP
4404 || raise_factor < -MAX_EXACT_10_EXP)) {
4405 realnum = algorithmR(fraction, raise_factor, realnum);
4406 }
4407 if (minusp) realnum = -realnum;
4408 return Scm_MakeFlonum(realnum);
4409 }
4410
4411 /* Entry point */
read_number(struct numread_packet * ctx)4412 static ScmObj read_number(struct numread_packet *ctx)
4413 {
4414 int radix_seen = 0, exactness_seen = 0, sign_seen = 0;
4415
4416 const char *str = ctx->buffer;
4417 int len = ctx->buflen;
4418
4419 #define CHK_EXACT_COMPLEX() \
4420 do { \
4421 if (ctx->exactness == EXACT) { \
4422 return numread_error("(exact complex number is not supported)", \
4423 ctx); \
4424 } \
4425 } while (0)
4426
4427 /* check suggested radix. */
4428 if (ctx->radix < SCM_RADIX_MIN || ctx->radix > SCM_RADIX_MAX) {
4429 return SCM_FALSE;
4430 }
4431
4432 /* start from prefix part */
4433 while (len >= 0) {
4434 if (*str != '#') break;
4435 str++;
4436 switch (*str++) {
4437 case 'x':; case 'X':;
4438 if (ctx->noradixprefix || radix_seen) return SCM_FALSE;
4439 ctx->radix = 16; radix_seen++;
4440 ctx->explicit = TRUE;
4441 len -= 2;
4442 continue;
4443 case 'o':; case 'O':;
4444 if (ctx->noradixprefix || radix_seen) return SCM_FALSE;
4445 ctx->radix = 8; radix_seen++;
4446 ctx->explicit = TRUE;
4447 len -= 2;
4448 continue;
4449 case 'b':; case 'B':;
4450 if (ctx->noradixprefix || radix_seen) return SCM_FALSE;
4451 ctx->radix = 2; radix_seen++;
4452 ctx->explicit = TRUE;
4453 len -= 2;
4454 continue;
4455 case 'd':; case 'D':;
4456 if (ctx->noradixprefix || radix_seen) return SCM_FALSE;
4457 ctx->radix = 10; radix_seen++;
4458 ctx->explicit = TRUE;
4459 len -= 2;
4460 continue;
4461 case 'e':; case 'E':;
4462 if (exactness_seen) return SCM_FALSE;
4463 ctx->exactness = EXACT; exactness_seen++;
4464 ctx->explicit = TRUE;
4465 len -= 2;
4466 continue;
4467 case 'i':; case 'I':;
4468 if (exactness_seen) return SCM_FALSE;
4469 ctx->exactness = INEXACT; exactness_seen++;
4470 ctx->explicit = TRUE;
4471 len -= 2;
4472 continue;
4473 case '0': case '1': case '2': case '3': case '4':
4474 case '5': case '6': case '7': case '8': case '9':
4475 if (ctx->noradixprefix || ctx->strict || radix_seen) return SCM_FALSE;
4476 else {
4477 ScmSize nread = 0;
4478 long radix = Scm_ParseDigitsAsLong(--str, --len, 10, &nread);
4479 if (radix < SCM_RADIX_MIN || radix > SCM_RADIX_MAX) return SCM_FALSE;
4480 str += nread; len -= nread;
4481 if (len <= 0) return SCM_FALSE;
4482 if (*str != 'r' && *str != 'R') return SCM_FALSE;
4483 str++; len--;
4484 ctx->radix = radix; radix_seen++;
4485 ctx->explicit = TRUE;
4486 continue;
4487 }
4488 }
4489 return SCM_FALSE;
4490 }
4491 if (len <= 0) return SCM_FALSE;
4492
4493 /* number body. need to check the special case of pure imaginary,
4494 and also eliminates some confusing cases. */
4495 if (*str == '+' || *str == '-') {
4496 if (len == 1) return SCM_FALSE;
4497 if (len == 2 && (str[1] == 'i' || str[1] == 'I')) {
4498 CHK_EXACT_COMPLEX();
4499 return Scm_MakeComplex(0.0, (*str == '+')? 1.0 : -1.0);
4500 }
4501 sign_seen = TRUE;
4502 }
4503
4504 ScmObj realpart = read_real(&str, &len, ctx);
4505 if (SCM_FALSEP(realpart) || len == 0) return realpart;
4506
4507 switch (*str) {
4508 case '@':
4509 /* polar representation of complex*/
4510 if (len <= 1) {
4511 return SCM_FALSE;
4512 } else {
4513 str++; len--;
4514 ScmObj angle = read_real(&str, &len, ctx);
4515 if (SCM_FALSEP(angle)) return SCM_FALSE;
4516 /* Gauche extension: X@Ypi */
4517 int pi_angle = (len == 2 && str[0] == 'p' && str[1] == 'i');
4518 if (!pi_angle && len != 0) return SCM_FALSE;
4519 CHK_EXACT_COMPLEX();
4520 double dmag = Scm_GetDouble(realpart);
4521 double dangle = Scm_GetDouble(angle);
4522 return (pi_angle
4523 ? Scm_MakeComplexPolarPI(dmag, dangle)
4524 : Scm_MakeComplexPolar(dmag, dangle));
4525 }
4526 case '+':
4527 case '-':
4528 /* rectangular representation of complex */
4529 if (len <= 1) {
4530 return SCM_FALSE;
4531 } else if (len == 2 && (str[1] == 'i' || str[1] == 'I')) {
4532 return Scm_MakeComplex(Scm_GetDouble(realpart),
4533 (*str == '+' ? 1.0 : -1.0));
4534 } else {
4535 ScmObj imagpart = read_real(&str, &len, ctx);
4536 if (SCM_FALSEP(imagpart) || len != 1
4537 || (*str != 'i' && *str != 'I')) {
4538 return SCM_FALSE;
4539 }
4540 CHK_EXACT_COMPLEX();
4541 if (Scm_Sign(imagpart) == 0) return realpart;
4542 return Scm_MakeComplex(Scm_GetDouble(realpart),
4543 Scm_GetDouble(imagpart));
4544 }
4545 case 'i':
4546 case 'I':
4547 /* '+' <ureal> 'i' or '-' <ureal> 'i' */
4548 if (!sign_seen || len != 1) return SCM_FALSE;
4549 CHK_EXACT_COMPLEX();
4550 if (Scm_Sign(realpart) == 0) return Scm_MakeFlonum(0.0);
4551 else return Scm_MakeComplex(0.0, Scm_GetDouble(realpart));
4552 default:
4553 return SCM_FALSE;
4554 }
4555 }
4556
numread_error(const char * msg,struct numread_packet * ctx)4557 static ScmObj numread_error(const char *msg, struct numread_packet *ctx)
4558 {
4559 if (ctx->throwerror) {
4560 Scm_Error("bad number format %s: %A", msg,
4561 Scm_MakeString(ctx->buffer, ctx->buflen,
4562 ctx->buflen, 0));
4563 }
4564 return SCM_FALSE;
4565 }
4566
4567 /* FLAGS is enum ScmNumberFormatFlags (see number.h). Only some of the
4568 flags are recognized for printing numbers. */
Scm_StringToNumber(ScmString * str,int radix,u_long flags)4569 ScmObj Scm_StringToNumber(ScmString *str, int radix, u_long flags)
4570 {
4571 ScmSmallInt len, size;
4572 const char *p = Scm_GetStringContent(str, &size, &len, NULL);
4573 if (size != len) {
4574 /* This can't be a proper number. */
4575 return SCM_FALSE;
4576 } else {
4577 struct numread_packet ctx;
4578 ctx.buffer = p;
4579 ctx.buflen = size;
4580 ctx.exactness = ((flags&SCM_NUMBER_FORMAT_EXACT)
4581 ? EXACT
4582 : ((flags&SCM_NUMBER_FORMAT_INEXACT)
4583 ? INEXACT
4584 : NOEXACT));
4585 ctx.padread = FALSE;
4586 ctx.explicit = FALSE;
4587 ctx.strict = flags&SCM_NUMBER_FORMAT_STRICT_R7RS;
4588 ctx.throwerror = FALSE;
4589 ctx.radix = radix;
4590 ctx.noradixprefix = flags&SCM_NUMBER_FORMAT_ALT_RADIX;
4591 return read_number(&ctx);
4592 }
4593 }
4594
4595 /*===============================================================
4596 * Initialization
4597 */
4598
4599 ScmObj Scm__ConstObjs[SCM_NUM_CONST_OBJS] = { SCM_FALSE };
4600
Scm__InitNumber(void)4601 void Scm__InitNumber(void)
4602 {
4603 ScmModule *mod = Scm_GaucheModule();
4604
4605 for (int radix = SCM_RADIX_MIN; radix <= SCM_RADIX_MAX; radix++) {
4606 longlimit[radix-SCM_RADIX_MIN] =
4607 (u_long)floor((double)LONG_MAX/radix - radix);
4608 /* Find max D where R^(D+1)-1 <= LONG_MAX */
4609 u_long n = 1;
4610 for (int i = 0; ; i++, n *= radix) {
4611 if (n >= (u_long)(LONG_MAX/radix)) {
4612 longdigs[radix-SCM_RADIX_MIN] = i-1;
4613 bigdig[radix-SCM_RADIX_MIN] = n;
4614 break;
4615 }
4616 }
4617 }
4618
4619 SCM_2_63 = Scm_Ash(SCM_MAKE_INT(1), 63);
4620 SCM_2_64 = Scm_Ash(SCM_MAKE_INT(1), 64);
4621 SCM_2_64_MINUS_1 = Scm_Sub(SCM_2_64, SCM_MAKE_INT(1));
4622 SCM_2_52 = Scm_Ash(SCM_MAKE_INT(1), 52);
4623 SCM_2_53 = Scm_Ash(SCM_MAKE_INT(1), 53);
4624 SCM_MINUS_2_63 = Scm_Negate(SCM_2_63);
4625 SCM_2_32 = Scm_Ash(SCM_MAKE_INT(1), 32);
4626 SCM_2_31 = Scm_Ash(SCM_MAKE_INT(1), 31);
4627 SCM_MINUS_2_31 = Scm_Negate(SCM_2_31);
4628 SCM_MIN_DENORMALIZED_FLONUM_EXACT =
4629 Scm_Reciprocal(Scm_Ash(SCM_MAKE_INT(1), 1075));
4630 SCM_MAX_FINITE_FLONUM_EXACT = Scm_Add(Scm_Sub(Scm_Ash(SCM_MAKE_INT(1), 1024),
4631 Scm_Ash(SCM_MAKE_INT(1), 971)),
4632 Scm_Sub(Scm_Ash(SCM_MAKE_INT(1), 970),
4633 SCM_MAKE_INT(1)));
4634
4635 SCM_POSITIVE_INFINITY = Scm_MakeFlonum(SCM_DBL_POSITIVE_INFINITY);
4636 SCM_NEGATIVE_INFINITY = Scm_MakeFlonum(SCM_DBL_NEGATIVE_INFINITY);
4637 SCM_NAN = Scm_MakeFlonum(SCM_DBL_NAN);
4638
4639 dexpt2_minus_52 = ldexp(1.0, -52);
4640 dexpt2_minus_53 = ldexp(1.0, -53);
4641
4642 Scm_InitBuiltinGeneric(&generic_add, "object-+", mod);
4643 Scm_InitBuiltinGeneric(&generic_sub, "object--", mod);
4644 Scm_InitBuiltinGeneric(&generic_mul, "object-*", mod);
4645 Scm_InitBuiltinGeneric(&generic_div, "object-/", mod);
4646
4647 #ifdef DOUBLE_ARMENDIAN
4648 check_armendian();
4649 #endif /*DOUBLE_ARMENDIAN*/
4650
4651 (void)SCM_INTERNAL_MUTEX_INIT(default_endian_mutex);
4652
4653 #ifdef COUNT_FLONUM_ALLOC
4654 Scm_AddCleanupHandler(report_flonum_count, NULL);
4655 #endif
4656 }
4657