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