1 /* number.c -*- mode:c; coding:utf-8; -*-
2 *
3 * Copyright (c) 2010-2021 Takashi Kato <ktakashi@ymail.com>
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 *
16 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
17 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
18 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
19 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
20 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
21 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
22 * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
23 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
24 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
25 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 *
28 * $Id: $
29 */
30 #include <string.h>
31 #include <float.h>
32 #include <math.h>
33 #include <ctype.h>
34 #define LIBSAGITTARIUS_BODY
35 #include "sagittarius/private/number.h"
36 #include "sagittarius/private/bignum.h"
37 #include "sagittarius/private/string.h"
38 #include "sagittarius/private/symbol.h"
39 #include "sagittarius/private/pair.h"
40 #include "sagittarius/private/port.h"
41 #include "sagittarius/private/core.h"
42 #include "sagittarius/private/exceptions.h"
43 #include "sagittarius/private/clos.h"
44 #include "sagittarius/private/unicode.h"
45 #include "sagittarius/private/error.h"
46 #include "sagittarius/private/numconst.h"
47 #include "sagittarius/private/arith.h"
48 #include "sagittarius/private/bits.h"
49 #include "sagittarius/private/values.h"
50 #include "sagittarius/private/vm.h"
51 #include "sagittarius/private/writer.h"
52
53 #ifndef M_PI
54 #define M_PI 3.14159265358979323846
55 #endif
56
57 /* Solaris doesn't have isinf */
58 #if defined(__SVR4) && defined(__sun)
59 # ifndef isinf
60 # include <ieeefp.h>
61 # define isinf(x) (!finite((x)) && (x) == (x))
62 # endif
63 #endif
64
65 #include "number.inc"
66
67 /* wrong type of argument */
68 #define wte4(who, t, g, irr) \
69 Sg_WrongTypeOfArgumentViolation(who, SG_MAKE_STRING(t), g, irr)
70 #define wte(who, t, g) wte4(who, t, g, g)
71
72 struct numread_packet {
73 const SgChar *buffer; /* original buffer */
74 int radix; /* radix */
75 int exactness; /* exactness; see enum below */
76 int padread; /* '#' padding has been read */
77 int strict; /* when true, reports an error if the
78 input violates implementation limitation;
79 otherwise, the routine returns #f. */
80 };
81
82 enum { /* used in the exactness flag */
83 NOEXACT, EXACT, INEXACT
84 };
85
86 #define IS_EXACT(ctx) ((ctx)->exactness == EXACT)
87 #define IS_INEXACT(ctx) ((ctx)->exactness == INEXACT)
88
89 #define RADIX_MIN 2
90 #define RADIX_MAX 36
91 #define P_EXP10 22 /* (floor (* 53 (log 2 5))) */
92 #define MAX_EXPONENT 325
93 #define MAX_EXACT_10_EXP (P_EXP10 + 1)
94
95 static const int64_t iexpt_2n52 = 0x10000000000000LL; // 2^(53-1)
96 static const int64_t iexpt_2n53 = 0x20000000000000LL; // 2^53
97
98 #include "roundeven.inc"
99
100 /* classes */
101 static SgClass *numeric_cpl[] = {
102 /* SG_CLASS_NUMBER, */
103 /* SG_CLASS_COMPLEX, */
104 /* SG_CLASS_REAL, */
105 /* SG_CLASS_RATIONAL, */
106 /* SG_CLASS_INTEGER, */
107 SG_CLASS_RATIONAL,
108 SG_CLASS_REAL,
109 SG_CLASS_COMPLEX,
110 SG_CLASS_NUMBER,
111 SG_CLASS_TOP,
112 NULL
113 };
114
115 static void number_print(SgObject obj, SgPort *port, SgWriteContext *ctx);
116
117 SG_DEFINE_BUILTIN_CLASS(Sg_NumberClass, number_print, NULL, NULL, NULL,
118 numeric_cpl+4);
119 SG_DEFINE_BUILTIN_CLASS(Sg_ComplexClass, number_print, NULL, NULL, NULL,
120 numeric_cpl+3);
121 SG_DEFINE_BUILTIN_CLASS(Sg_RealClass, number_print, NULL, NULL, NULL,
122 numeric_cpl+2);
123 SG_DEFINE_BUILTIN_CLASS(Sg_RationalClass, number_print, NULL, NULL, NULL,
124 numeric_cpl+1);
125 SG_DEFINE_BUILTIN_CLASS(Sg_IntegerClass, number_print, NULL, NULL, NULL,
126 numeric_cpl);
127
ipow(int r,long n)128 static inline unsigned long ipow(int r, long n)
129 {
130 unsigned long k;
131 for (k = 1; n > 0; n--) k *= r;
132 return k;
133 }
134
pow10n(double x,long n)135 static double pow10n(double x, long n)
136 {
137 static double dpow10[] = { 1.0, 1.0e1, 1.0e2, 1.0e3, 1.0e4,
138 1.0e5, 1.0e6, 1.0e7, 1.0e8, 1.0e9,
139 1.0e10, 1.0e11, 1.0e12, 1.0e13, 1.0e14,
140 1.0e15, 1.0e16, 1.0e17, 1.0e18, 1.0e19,
141 1.0e20, 1.0e21, 1.0e22, 1.0e23 };
142 if (n >= 0) {
143 while (n > 23) {
144 x *= 1.0e24;
145 n -= 24;
146 }
147 return x*dpow10[n];
148 } else {
149 while (n < -23) {
150 x /= 1.0e24;
151 n += 24;
152 }
153 return x/dpow10[-n];
154 }
155
156 }
157
int64_to_integer(int64_t value)158 static SgObject int64_to_integer(int64_t value)
159 {
160 if ((value <= SG_INT_MAX) && (value >= SG_INT_MIN)) return SG_MAKE_INT(value);
161 return Sg_MakeBignumFromS64(value);
162 }
163
oprtr_norm_integer(SgObject obj)164 static SgObject oprtr_norm_integer(SgObject obj)
165 {
166 ASSERT(SG_INTP(obj) || SG_BIGNUMP(obj));
167 if (SG_BIGNUMP(obj)) return Sg_BignumToInteger(SG_BIGNUM(obj));
168 return obj;
169 }
170
oprtr_norm_complex(SgObject real,SgObject imag)171 static SgObject oprtr_norm_complex(SgObject real, SgObject imag)
172 {
173 ASSERT(!SG_COMPLEXP(real));
174 ASSERT(!SG_COMPLEXP(imag));
175 if (SG_INTP(imag) && SG_EQ(imag, SG_MAKE_INT(0))) return real;
176 if (SG_BIGNUMP(imag) && SG_BIGNUM_GET_SIGN(SG_BIGNUM(imag)) == 0) return real;
177 if (SG_FLONUMP(real) || SG_FLONUMP(imag)) {
178 return Sg_MakeComplex(Sg_Inexact(real), Sg_Inexact(imag));
179 }
180 return Sg_MakeComplex(real, imag);
181 }
182
oprtr_expt(SgObject lhs,long n)183 static SgObject oprtr_expt(SgObject lhs, long n)
184 {
185 SgObject r;
186 if (n == 0) return SG_MAKE_INT(1);
187 if (n == 1) return lhs;
188 if (n < 0) return Sg_Inverse(oprtr_expt(lhs, -n));
189 if (!SG_COMPLEXP(lhs) && Sg_NegativeP(lhs)) {
190 SgObject ans = oprtr_expt(Sg_Negate(lhs), n);
191 if (n & 1) return Sg_Negate(ans);
192 return ans;
193 }
194 if (lhs == SG_MAKE_INT(0)) return lhs;
195 if (lhs == SG_MAKE_INT(1)) return lhs;
196 if (lhs == SG_MAKE_INT(2)) {
197 if (n + 1 <= SG_INT_SIZE) return SG_MAKE_INT(1UL << n);
198 return Sg_Ash(SG_MAKE_INT(1), n);
199 }
200
201 if (SG_RATIONALP(lhs)) {
202 return Sg_MakeRational(oprtr_expt(SG_RATIONAL(lhs)->numerator, n),
203 oprtr_expt(SG_RATIONAL(lhs)->denominator, n));
204 }
205 /* bignum */
206 if (SG_BIGNUMP(lhs)) {
207 return Sg_BignumExpt(SG_BIGNUM(lhs), n);
208 }
209 r = SG_MAKE_INT(1);
210 while (TRUE) {
211 if (n & 1) {
212 if (r == SG_MAKE_INT(1)) r = lhs;
213 else r = Sg_Mul(r, lhs);
214 if (n == 1) return r;
215 }
216 lhs = Sg_Mul(lhs, lhs);
217 n >>= 1;
218 }
219 }
220
decode_double(double n,int * exp,int * sign)221 static int64_t decode_double(double n, int *exp, int *sign)
222 {
223 union { double f64; uint64_t u64; } datum;
224 uint64_t bits;
225 uint64_t mant_bits;
226 uint32_t sign_bits;
227 uint32_t exp_bits;
228 datum.f64 = n;
229 bits = datum.u64;
230 mant_bits = bits & (iexpt_2n52 - 1);
231 sign_bits = bits >> 63;
232 exp_bits = (bits >> 52) & 0x7ff;
233
234 if (n == 0.0) {
235 *exp = 0;
236 *sign = sign_bits ? -1 : 1;
237 return 0;
238 }
239 if (isnan(n)) {
240 *exp = 972;
241 *sign = 1;
242 return 0x18000000000000LL; /* (uint64_t)0x180000 << 32; */
243 }
244 if (isinf(n)) {
245 *exp = 972;
246 *sign = sign_bits ? -1 : 1;
247 return 0x10000000000000LL; /* (uint64_t)0x100000 << 32; */
248 }
249 ASSERT(exp_bits != 0x7ff);
250 *exp = (exp_bits ? (int)exp_bits - 1023 : -1022) - 52;
251 *sign = sign_bits ? -1 : 1;
252 if (exp_bits) mant_bits |= iexpt_2n52;
253 return mant_bits;
254 }
255
integer_init_n_alloc(int64_t m,int shift_left)256 static SgObject integer_init_n_alloc(int64_t m, int shift_left)
257 {
258 ASSERT(m >= 0);
259 if (m == 0) return SG_MAKE_INT(0);
260 else {
261 SgObject b = Sg_MakeBignumFromS64(m);
262 return Sg_BignumShiftLeft(SG_BIGNUM(b), shift_left);
263 }
264 }
265
nextfloat(double z)266 static double nextfloat(double z)
267 {
268 int k, sign;
269 int64_t m = decode_double(z, &k, &sign);
270 ASSERT(sign >= 0);
271 if (m == iexpt_2n53 - 1) return ldexp((double)iexpt_2n52, k + 1);
272 return ldexp((double)(m + 1), k);
273 }
274
prevfloat(double z)275 static double prevfloat(double z)
276 {
277 int k, sign;
278 int64_t m = decode_double(z, &k, &sign);
279 ASSERT(sign >= 0);
280 if (k > -1074 && m == iexpt_2n52) return ldexp((double)(iexpt_2n53 - 1), k - 1);
281 return ldexp((double)(m - 1), k);
282 }
283
284 /*
285 // Reference:
286 // William D. Clinger.
287 // How to read floating point numbers accurately
288 // Proceedings of the ACM SIGPLAN 1990 conference on Programming language design and implementation, p.92-101, June 1990
289 */
algorithmR(SgObject f,const long e,const double z0)290 static double algorithmR(SgObject f, const long e, const double z0)
291 {
292 double z = z0;
293 SgObject x0, pow10e;
294 if (e >= 0) {
295 x0 = Sg_Mul(f, Sg_Expt(SG_MAKE_INT(10), SG_MAKE_INT(e)));
296 pow10e = SG_UNDEF;
297 } else {
298 x0 = SG_UNDEF;
299 pow10e = Sg_Expt(SG_MAKE_INT(10), SG_MAKE_INT(-e));
300 }
301 while (1) {
302 int k, sign, test;
303 int64_t m;
304 SgObject x, y, D, D2;
305 int negP;
306 if (isinf(z)) return z;
307 m = decode_double(z, &k, &sign);
308 ASSERT(sign >= 0);
309 if (e >= 0) {
310 if (k >= 0) {
311 x = x0;
312 y = integer_init_n_alloc(m, k);
313 } else {
314 x = Sg_Ash(x0, -k);
315 y = int64_to_integer(m);
316 }
317 } else {
318 if (k >= 0) {
319 x = f;
320 y = Sg_Mul(integer_init_n_alloc(m, k), pow10e);
321 } else {
322 x = Sg_Ash(f, -k);
323 y = Sg_Mul(int64_to_integer(m), pow10e);
324 }
325 }
326 D = Sg_Sub(x, y);
327 D2 = Sg_Mul(int64_to_integer(m + m), D);
328 negP = Sg_NegativeP(D);
329 if (negP) {
330 if (SG_BIGNUMP(D2)) SG_BIGNUM_SET_SIGN(SG_BIGNUM(D2), -SG_BIGNUM_GET_SIGN(D2));
331 else D2 = SG_MAKE_INT(-SG_INT_VALUE(D2));
332 }
333 test = Sg_NumCmp(D2, y);
334 if (test < 0) {
335 if (negP &&
336 m == iexpt_2n52 &&
337 k > -1074 &&
338 Sg_NumGt(Sg_Ash(D2, 1), y)) {
339 z = prevfloat(z);
340 continue;
341 }
342 return z;
343 }
344 if (test == 0) {
345 if ((m & 1) == 0) {
346 if (negP && m == iexpt_2n52) {
347 z = prevfloat(z);
348 continue;
349 }
350 return z;
351 }
352 return negP ? prevfloat(z) : nextfloat(z);
353 }
354 z = negP ? prevfloat(z) : nextfloat(z);;
355 }
356 }
357
number_read_error(const SgChar * msg,struct numread_packet * context)358 static SgObject number_read_error(const SgChar *msg,
359 struct numread_packet *context)
360 {
361 if (context->strict) {
362 SgObject irr = Sg_HeapString(context->buffer);
363 SgObject m = Sg_Sprintf(UC("%s: %A"), msg, irr);
364 /* FIXME wasting memory. */
365 SgObject err = Sg_SimpleConditions(Sg_MakeReaderCondition(m));
366 err = Sg_Cons(Sg_MakeIrritantsCondition(irr), err);
367 err = Sg_Cons(Sg_MakeWhoCondition(SG_INTERN("read")), err);
368 err = Sg_Cons(Sg_MakeImplementationRestrictionViolation(), err);
369 Sg_Raise(Sg_Condition(err), FALSE);
370 }
371 return SG_FALSE;
372 }
373
374
375 static long longdigs[RADIX_MAX - RADIX_MIN + 1] = {0};
376 static unsigned long longlimit[RADIX_MAX - RADIX_MIN + 1] = {0};
377 static unsigned long bigdig[RADIX_MAX - RADIX_MIN + 1] = {0};
378
read_uint(const SgChar ** strp,long * lenp,struct numread_packet * ctx,SgObject initval)379 static SgObject read_uint(const SgChar **strp, long *lenp,
380 struct numread_packet *ctx,
381 SgObject initval)
382 {
383 const SgChar *str = *strp;
384 int digread = FALSE;
385 long len = *lenp;
386 int radix = ctx->radix;
387 long digits = 0, diglimit = longdigs[radix - RADIX_MIN];
388 unsigned long limit = longlimit[radix - RADIX_MIN],
389 bdig = bigdig[radix - RADIX_MIN];
390 unsigned long value_int = 0;
391 SgBignum *value_big = NULL;
392 SgChar c;
393 static const char tab[] = "0123456789abcdefghijklmnopqrstuvwxyz";
394 const char *ptab;
395
396 if (!SG_FALSEP(initval)) {
397 if (SG_INTP(initval)) {
398 if ((unsigned long)SG_INT_VALUE(initval) > limit) {
399 value_big = Sg_MakeBignumWithSize(4, SG_INT_VALUE(initval));
400 } else {
401 value_int = SG_INT_VALUE(initval);
402 }
403 } else if (SG_BIGNUMP(initval)) {
404 value_big = SG_BIGNUM(Sg_BignumCopy(SG_BIGNUM(initval)));
405 }
406 digread = TRUE;
407 } else if (*str == '0') {
408 /* Ignore leading 0's, to avoid unnecessary bignum operations. */
409 while (len > 0 && *str == '0') { str++; len--; }
410 digread = TRUE;
411 }
412
413 while (len--) {
414 int digval = -1;
415 c = tolower(*str++);
416 if (ctx->padread) {
417 if (c == '#') digval = 0;
418 else break;
419 } else if (digread && c == '#') {
420 digval = 0;
421 ctx->padread = TRUE;
422 if (ctx->exactness == NOEXACT) {
423 ctx->exactness = INEXACT;
424 }
425 } else {
426 for (ptab = tab; ptab < tab + radix; ptab++) {
427 if (c == *ptab) {
428 digval = (int)(ptab - tab);
429 digread = TRUE;
430 break;
431 }
432 }
433 }
434
435 if (digval < 0) break;
436 value_int = value_int * radix + digval;
437 digits++;
438 if (value_big == NULL) {
439 if (value_int >= limit) {
440 value_big = Sg_MakeBignumWithSize(4, value_int);
441 value_int = digits = 0;
442 }
443 } else if (digits > diglimit) {
444 value_big = Sg_BignumAccMultAddUI(value_big, bdig, value_int);
445 value_int = digits = 0;
446 }
447 }
448 *strp = str - 1;
449 *lenp = len + 1;
450 if (value_big == NULL) return Sg_MakeInteger(value_int);
451 if (digits > 0) {
452 value_big = Sg_BignumAccMultAddUI(value_big,
453 ipow(radix, digits),
454 value_int);
455 }
456 return Sg_NormalizeBignum(value_big);
457 }
458
read_real(const SgChar ** strp,long * lenp,struct numread_packet * ctx)459 static SgObject read_real(const SgChar **strp, long *lenp,
460 struct numread_packet *ctx)
461 {
462 int minusp = FALSE, exp_minusp = FALSE, exp_overflow = FALSE;
463 int sign_seen = FALSE, has_fraction = FALSE, has_exponent = FALSE;
464 long fracdigs = 0;
465 long exponent = 0;
466 SgObject intpart, fraction;
467 const SgChar *mark;
468
469 switch (**strp) {
470 case '-': minusp = TRUE;
471 case '+':
472 (*strp)++; (*lenp)--; sign_seen = TRUE;
473 }
474 if ((*lenp) <= 0) return SG_FALSE;
475 mark = *strp;
476 /* recognize specials */
477 if (sign_seen && (*lenp) >= 5) {
478 if (ustrncmp(*strp, "inf.0", 5) == 0) {
479 (*strp) += 5; (*lenp) -=5;
480 return minusp ? SG_NEGATIVE_INFINITY : SG_POSITIVE_INFINITY;
481 }
482 if (ustrncmp(*strp, "nan.0", 5) == 0) {
483 (*strp) += 5; (*lenp) -=5;
484 return SG_NAN;
485 }
486 }
487 /* read integral part */
488 if (**strp != '.') {
489 intpart = read_uint(strp, lenp, ctx, SG_FALSE);
490 if ((*lenp) <= 0) {
491 if (minusp) intpart = Sg_Negate(intpart);
492 if (IS_INEXACT(ctx)) {
493 return Sg_Inexact(intpart);
494 } else {
495 return intpart;
496 }
497 }
498 if (**strp == '/') {
499 /* possibly rational */
500 SgObject denom;
501 long lensave;
502 if ((*lenp) <= 1 || mark == *strp) return SG_FALSE;
503 (*strp)++; (*lenp)--;
504 lensave = *lenp;
505 denom = read_uint(strp, lenp, ctx, SG_FALSE);
506 if (SG_FALSEP(denom)) return SG_FALSE;
507 if (SG_MAKE_INT(0) == denom) {
508 if (lensave > *lenp) {
509 if (!IS_INEXACT(ctx)) {
510 return number_read_error(UC("exact infinity/nan is not supported."),
511 ctx);
512 }
513 if (!SG_VM_IS_SET_FLAG(Sg_VM(), SG_R6RS_MODE)) {
514 if (SG_MAKE_INT(0) == intpart) return SG_NAN;
515 return minusp ? SG_NEGATIVE_INFINITY : SG_POSITIVE_INFINITY;
516 } else {
517 return SG_FALSE;
518 }
519 } else {
520 return SG_FALSE;
521 }
522 }
523 if (minusp) intpart = Sg_Negate(intpart);
524 if (IS_INEXACT(ctx)) {return Sg_Inexact(Sg_Div(intpart, denom));
525 } else {
526 return Sg_MakeRational(intpart, denom);
527 }
528 }
529 } else {
530 intpart = SG_FALSE; /* indicate there was no intpart */
531 }
532
533 /* read fractional part */
534 if (**strp == '.') {
535 long lensave;
536 if (ctx->radix != 10) {
537 return number_read_error(UC("only 10-based fraction is supported"), ctx);
538 }
539 (*strp)++; (*lenp)--;
540 lensave = *lenp;
541 fraction = read_uint(strp, lenp, ctx, intpart);
542 fracdigs = lensave - *lenp;
543 has_fraction = TRUE;
544 } else {
545 fraction = intpart;
546 }
547
548 if (SG_FALSEP(intpart)) {
549 if (fracdigs == 0) return SG_FALSE; /* input was '.' */
550 }
551 if (mark == *strp) return SG_FALSE;
552 /* read exponent */
553 if (*lenp > 0 && strchr("eEsSfFdDlL", (int)**strp)) {
554 (*strp)++;
555 if (**strp == '|') return SG_FALSE; /* `10.2e|43 */
556 if (--(*lenp) <= 0) return SG_FALSE;
557 switch (**strp) {
558 case '-': exp_minusp = TRUE;
559 case '+':
560 (*strp)++;
561 if (--(*lenp) <= 0) return SG_FALSE;
562 }
563 while (*lenp > 0) {
564 SgChar c = **strp;
565 if (!isdigit(c)) break;
566 (*strp)++, (*lenp)--;
567 if (isdigit(c) && !exp_overflow) {
568 exponent = exponent * 10 + (c - '0');
569 if (exponent >= MAX_EXPONENT) exp_overflow = TRUE;
570 }
571 }
572 if (exp_minusp) exponent = -exponent;
573 has_exponent = TRUE;
574 }
575 /* parse precision */
576 if (**strp == '|') {
577 SgObject pre;
578 (*strp)++, (*lenp)--;
579 /* just ignore */
580 pre = read_uint(strp, lenp, ctx, SG_FALSE);
581 if (!SG_FALSEP(pre)) ctx->exactness = INEXACT;
582 }
583
584 if (exp_overflow && IS_INEXACT(ctx)) {
585 if (exp_minusp) {
586 return Sg_MakeFlonum(0.0);
587 } else {
588 return minusp ? SG_NEGATIVE_INFINITY : SG_POSITIVE_INFINITY;
589 }
590 }
591 if (IS_EXACT(ctx)) {
592 /* explicit exact number. */
593 SgObject e = Sg_Mul(fraction, Sg_Expt(SG_MAKE_INT(10),
594 Sg_MakeInteger(exponent - fracdigs)));
595 if (minusp) return Sg_Negate(e);
596 else return e;
597 } else if (ctx->exactness == NOEXACT &&
598 !has_fraction && !has_exponent) {
599 return (minusp) ? Sg_Negate(intpart) : intpart;
600 } else {
601 double realnum = Sg_GetDouble(fraction);
602 realnum = pow10n(realnum, exponent - fracdigs);
603 if (isinf(realnum)) {
604 return minusp ? SG_NEGATIVE_INFINITY : SG_POSITIVE_INFINITY;
605 }
606 if (realnum > 0.0
607 && (Sg_NumCmp(fraction, SG_2_52) > 0
608 || exponent - fracdigs > MAX_EXACT_10_EXP
609 || exponent - fracdigs < -MAX_EXACT_10_EXP)) {
610 realnum = algorithmR(fraction, exponent - fracdigs, realnum);
611 }
612 if (minusp) realnum = -realnum;
613 return Sg_MakeFlonum(realnum);
614 }
615 }
616 /*
617 * Number Parser
618 *
619 * <number> : <prefix> <complex>
620 * <prefix> : <radix> <exactness> | <exactness> <radix>
621 * <radix> : <empty> | '#b' | '#o' | '#d' | '#x'
622 * <exactness> : <empty> | '#e' | '#i'
623 * <complex> : <real>
624 * | <real> '@' <real>
625 * | <real> '+' <ureal> 'i'
626 * | <real> '-' <ureal> 'i'
627 * | <real> '+' 'i'
628 * | <real> '-' 'i'
629 * | '+' <ureal> 'i'
630 * | '-' <ureal> 'i'
631 * | '+' 'i'
632 * | '-' 'i'
633 * <real> : <sign> <ureal>
634 * <sign> : <empty> | '+' | '-'
635 * <ureal> : <uinteger>
636 * | <uinteger> '/' <uinteger>
637 * | <decimal>
638 * <uinteger> : <digit>+ '#'*
639 * <decimal> : <digit10>+ '#'* <suffix>
640 * | '.' <digit10>+ '#'* <suffix>
641 * | <digit10>+ '.' <digit10>+ '#'* <suffix>
642 * | <digit10>+ '#'+ '.' '#'* <suffix>
643 * <suffix> : <empty> | <exponent-marker> <sign> <digit10>+
644 * <exponent-marker> : 'e' | 's' | 'f' | 'd' | 'l'
645 *
646 * The parser reads characters from on-memory buffer.
647 * Multibyte strings are filtered out in the early stage of
648 * parsing, so the subroutines assume the buffer contains
649 * only ASCII chars.
650 */
651 /* entry point */
read_number(const SgChar * str,long len,int radix,int strict)652 static SgObject read_number(const SgChar *str, long len, int radix, int strict)
653 {
654 struct numread_packet ctx;
655 int radix_seen = 0, exactness_seen = 0, sign_seen = 0;
656 SgObject realpart;
657
658 ctx.buffer = str;
659 ctx.exactness = NOEXACT;
660 ctx.padread = FALSE;
661 ctx.strict = strict;
662
663 if (radix <= 1 || radix > 36) return SG_FALSE;
664 ctx.radix = radix;
665 /* read radix and prefix*/
666 for (; len >= 0; len -= 2) {
667 if (*str != '#') break;
668 str++;
669 switch (*str++) {
670 case 'x':; case 'X':;
671 if (radix_seen) return SG_FALSE;
672 ctx.radix = 16; radix_seen++;
673 continue;
674 case 'o':; case 'O':;
675 if (radix_seen) return SG_FALSE;
676 ctx.radix = 8; radix_seen++;
677 continue;
678 case 'b':; case 'B':;
679 if (radix_seen) return SG_FALSE;
680 ctx.radix = 2; radix_seen++;
681 continue;
682 case 'd':; case 'D':;
683 if (radix_seen) return SG_FALSE;
684 ctx.radix = 10; radix_seen++;
685 continue;
686 case 'e':; case 'E':;
687 if (exactness_seen) return SG_FALSE;
688 ctx.exactness = EXACT; exactness_seen++;
689 continue;
690 case 'i':; case 'I':;
691 if (exactness_seen) return SG_FALSE;
692 ctx.exactness = INEXACT; exactness_seen++;
693 continue;
694 }
695 return SG_FALSE;
696 }
697 if (len <= 0) return SG_FALSE;
698
699 if (*str == '+' || *str == '-') {
700 if (len == 1) return SG_FALSE;
701 if (len == 2 && (str[1] == 'i' || str[1] == 'I')) {
702 if (IS_INEXACT(&ctx)) {
703 return Sg_MakeComplex(Sg_MakeFlonum(0.0),
704 Sg_MakeFlonum((*str == '+') ? 1.0 : -1.0));
705 } else {
706 return Sg_MakeComplex(SG_MAKE_INT(0),
707 SG_MAKE_INT((*str == '+') ? 1 : -1));
708 }
709 }
710 sign_seen = TRUE;
711 }
712
713 realpart = read_real(&str, &len, &ctx);
714 if (SG_FALSEP(realpart) || len == 0) return realpart;
715
716 switch (*str) {
717 case '@':
718 /* polar representation of complex */
719 if (len <= 1) {
720 return SG_FALSE;
721 } else {
722 SgObject angle;
723 str++; len--;
724 angle = read_real(&str, &len, &ctx);
725 if (SG_FALSEP(angle) || len != 0) return SG_FALSE;
726 if (IS_EXACT(&ctx)) {
727 /* follow the decision of Chez, it does make sense.
728 c.f. #e1@1 can't be represent as exact number but
729 it should not return inexact number since the
730 prefix specifies exact number. so just raise
731 an exception.
732 */
733 return number_read_error(UC("cannot represent given number"), &ctx);
734 } else {
735 return Sg_MakeComplexPolar(realpart, angle);
736 }
737 }
738 case '+':
739 case '-':
740 /* rectangular representation of complex */
741 if (len <= 1) {
742 return SG_FALSE;
743 } else if (len == 2 && (str[1] == 'i' || str[1] == 'I')) {
744 if (IS_INEXACT(&ctx)) {
745 return Sg_MakeComplex(realpart,
746 Sg_MakeFlonum((*str == '+') ? 1.0 : -1.0));
747 } else {
748 return Sg_MakeComplex(realpart,
749 SG_MAKE_INT((*str == '+') ? 1 : -1));
750 }
751 } else {
752 SgObject imagpart;
753 ctx.exactness = NOEXACT;
754 imagpart = read_real(&str, &len, &ctx);
755 if (SG_FALSEP(imagpart) || len != 1 || (*str != 'i' && *str != 'I')) {
756 return SG_FALSE;
757 }
758 /* if (Sg_Sign(imagpart) == 0) return realpart; */
759 return Sg_MakeComplex(realpart, imagpart);
760 }
761 case 'i':
762 case 'I':
763 /* '+' <ureal> 'i' or '-' <ureal> 'i' */
764 if (!sign_seen || len != 1) return SG_FALSE;
765 if (Sg_Sign(realpart) == 0) return realpart;
766 else {
767 if (IS_INEXACT(&ctx)) return Sg_MakeComplex(Sg_MakeFlonum(0.0), realpart);
768 else return Sg_MakeComplex(SG_MAKE_INT(0), realpart);
769 }
770 default:
771 return SG_FALSE;
772 }
773 }
774
Sg_MakeInteger(long x)775 SgObject Sg_MakeInteger(long x)
776 {
777 if (x >= SG_INT_MIN && x <= SG_INT_MAX) {
778 return SG_MAKE_INT(x);
779 } else {
780 return Sg_MakeBignumFromSI(x);
781 }
782 }
783
Sg_MakeIntegerU(unsigned long x)784 SgObject Sg_MakeIntegerU(unsigned long x)
785 {
786 if (x <= (unsigned long)SG_INT_MAX) {
787 return SG_MAKE_INT(x);
788 } else {
789 return Sg_MakeBignumFromUI(x);
790 }
791 }
792
Sg_MakeIntegerFromS64(int64_t x)793 SgObject Sg_MakeIntegerFromS64(int64_t x)
794 {
795 if ((x <= SG_INT_MAX) && (x >= SG_INT_MIN)) return SG_MAKE_INT(x);
796 return Sg_MakeBignumFromS64(x);
797 }
798
Sg_MakeIntegerFromU64(uint64_t x)799 SgObject Sg_MakeIntegerFromU64(uint64_t x)
800 {
801 if (x <= SG_INT_MAX) return SG_MAKE_INT(x);
802 return Sg_MakeBignumFromU64(x);
803 }
804
range_err(SgObject obj,int clamp,int * oor)805 static void range_err(SgObject obj, int clamp, int *oor)
806 {
807 if (clamp == SG_CLAMP_NONE && oor != NULL) {
808 *oor = TRUE;
809 } else {
810 /* TODO should we make invalid argument error? */
811 Sg_Error(UC("argument out of range: %S"), obj);
812 }
813 }
814
Sg_GetIntegerClamp(SgObject obj,int clamp,int * oor)815 long Sg_GetIntegerClamp(SgObject obj, int clamp, int *oor)
816 {
817 double v = 0.0;
818 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = FALSE;
819 if (SG_INTP(obj)) return SG_INT_VALUE(obj);
820 else if (SG_BIGNUMP(obj)) {
821 return Sg_BignumToSI(SG_BIGNUM(obj), clamp, oor);
822 }
823 else if (SG_FLONUMP(obj)) {
824 v = SG_FLONUM_VALUE(obj);
825 goto flonum;
826 }
827 else if (SG_RATIONALP(obj)) {
828 v = Sg_GetDouble(obj);
829 goto flonum;
830 }
831 else {
832 goto err;
833 }
834 flonum:
835 if (v > (double)LONG_MAX) {
836 if (clamp & SG_CLAMP_HI) return LONG_MAX;
837 else goto err;
838 }
839 if (v < (double)LONG_MIN) {
840 if (clamp & SG_CLAMP_LO) return LONG_MIN;
841 else goto err;
842 }
843 return (long)v;
844 err:
845 range_err(obj, clamp, oor);
846 return 0;
847 }
848
Sg_GetUIntegerClamp(SgObject obj,int clamp,int * oor)849 unsigned long Sg_GetUIntegerClamp(SgObject obj, int clamp, int *oor)
850 {
851 double v = 0.0;
852 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = FALSE;
853 if (SG_INTP(obj)) {
854 if (SG_INT_VALUE(obj) < 0) {
855 if (clamp & SG_CLAMP_LO) return 0;
856 else goto err;
857 }
858 return SG_INT_VALUE(obj);
859 }
860 else if (SG_BIGNUMP(obj)) {
861 return Sg_BignumToUI(SG_BIGNUM(obj), clamp, oor);
862 }
863 else if (SG_FLONUMP(obj)) {
864 v = SG_FLONUM_VALUE(obj);
865 goto flonum;
866 }
867 else if (SG_RATIONALP(obj)) {
868 v = Sg_GetDouble(obj);
869 goto flonum;
870 }
871 else {
872 goto err;
873 }
874 flonum:
875 if (v > (double)ULONG_MAX) {
876 if (clamp & SG_CLAMP_HI) return ULONG_MAX;
877 else goto err;
878 }
879 if (v < 0.0) {
880 if (clamp & SG_CLAMP_LO) return 0;
881 else goto err;
882 }
883 return (unsigned long)v;
884 err:
885 range_err(obj, clamp, oor);
886 return 0;
887 }
888
889 #if SIZEOF_LONG < 8
890 #define SG_SET_INT64_MAX(v64) \
891 ((v64) = ((((int64_t)LONG_MAX)<<32) + (int64_t)ULONG_MAX))
892 #define SG_SET_INT64_MIN(v64) \
893 ((v64) = (((int64_t)LONG_MAX + 1) << 32))
894 #define SG_SET_UINT64_MAX(v64) \
895 ((v64) = ((((uint64_t)ULONG_MAX) << 32) + (uint64_t)ULONG_MAX))
896
Sg_GetIntegerS64Clamp(SgObject obj,int clamp,int * oor)897 int64_t Sg_GetIntegerS64Clamp(SgObject obj, int clamp, int *oor)
898 {
899 int64_t r = 0;
900 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = FALSE;
901 if (SG_INTP(obj)) return (int64_t)SG_INT_VALUE(obj);
902 if (SG_BIGNUMP(obj)) {
903 return Sg_BignumToS64(SG_BIGNUM(obj), clamp, oor);
904 }
905 if (SG_RATIONALP(obj)) {
906 obj = Sg_Inexact(obj);
907 /* fall through */
908 }
909 if (SG_FLONUMP(obj)) {
910 int64_t maxval, minval;
911 double v;
912 SG_SET_INT64_MAX(maxval);
913 SG_SET_INT64_MIN(minval);
914 v = SG_FLONUM_VALUE(obj);
915 if (v > (double)maxval) {
916 if (!(clamp & SG_CLAMP_HI)) goto err;
917 return maxval;
918 } else if (v < (double)minval) {
919 if (!(clamp & SG_CLAMP_LO)) goto err;
920 return minval;
921 } else {
922 return (int64_t)v;
923 }
924 }
925 err:
926 range_err(obj, clamp, oor);
927 return r;
928 }
929
Sg_GetIntegerU64Clamp(SgObject obj,int clamp,int * oor)930 uint64_t Sg_GetIntegerU64Clamp(SgObject obj, int clamp, int *oor)
931 {
932 uint64_t r = 0;
933 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = FALSE;
934 if (SG_INTP(obj)) {
935 long v = SG_INT_VALUE(obj);
936 if (v < 0) {
937 if (!(clamp & SG_CLAMP_LO)) goto err;
938 return 0;
939 } else {
940 return (uint64_t)v;
941 }
942 }
943 if (SG_BIGNUMP(obj)) {
944 return Sg_BignumToU64(SG_BIGNUM(obj), clamp, oor);
945 }
946 if (SG_RATIONALP(obj)) {
947 obj = Sg_Inexact(obj);
948 /* fall through */
949 }
950 if (SG_FLONUMP(obj)) {
951 double v = SG_FLONUM_VALUE(obj);
952 uint64_t maxval;
953 if (v < 0) {
954 if (!(clamp & SG_CLAMP_LO)) goto err;
955 return 0;
956 }
957
958 SG_SET_UINT64_MAX(maxval);
959
960 if (v > (double)maxval) {
961 if (!(clamp & SG_CLAMP_HI)) goto err;
962 return maxval;
963 } else {
964 return (uint64_t)v;
965 }
966 }
967 err:
968 range_err(obj, clamp, oor);
969 return r;
970 }
971 #endif
972
973
make_rational(SgObject nume,SgObject deno)974 static SgRational* make_rational(SgObject nume, SgObject deno)
975 {
976 SgRational *z;
977 z = SG_NEW(SgRational);
978 SG_SET_CLASS(z, SG_CLASS_RATIONAL);
979 z->numerator = nume;
980 z->denominator = deno;
981 return z;
982 }
983
Sg_MakeRational(SgObject numerator,SgObject denominator)984 SgObject Sg_MakeRational(SgObject numerator, SgObject denominator)
985 {
986 SgRational *z;
987 if(!Sg_ExactP(numerator)) {
988 Sg_AssertionViolation(
989 SG_FALSE,
990 Sg_Sprintf(UC("numerator must be an exact integer, but got %S"),
991 numerator),
992 numerator);
993 }
994 if(!Sg_ExactP(denominator)) {
995 Sg_AssertionViolation(
996 SG_FALSE,
997 Sg_Sprintf(UC("denominator must be an exact integer, but got %S"),
998 denominator),
999 denominator);
1000 }
1001 if (denominator == SG_MAKE_INT(0)) {
1002 Sg_AssertionViolation(SG_FALSE, SG_MAKE_STRING("undefined for 0"),
1003 SG_LIST2(numerator, denominator));
1004 }
1005
1006 if (denominator == SG_MAKE_INT(1)) return numerator;
1007 if (numerator == SG_MAKE_INT(0)) return SG_MAKE_INT(0);
1008
1009 z = make_rational(numerator, denominator);
1010 return Sg_ReduceRational(z);
1011 }
1012
Sg_ReduceRational(SgObject rational)1013 SgObject Sg_ReduceRational(SgObject rational)
1014 {
1015 SgObject numer, denom;
1016 SgObject common;
1017 int negated = FALSE;
1018
1019 if (SG_INTP(rational) || SG_BIGNUMP(rational)) return rational;
1020 if (!SG_RATIONALP(rational)) {
1021 wte(SG_FALSE, "exect rational number", rational);
1022 }
1023 numer = SG_RATIONAL(rational)->numerator;
1024 denom = SG_RATIONAL(rational)->denominator;
1025
1026 if (Sg_Sign(denom) < 0) {
1027 numer = Sg_Negate(numer);
1028 denom = Sg_Negate(denom);
1029 negated = TRUE;
1030 }
1031
1032 if (denom == SG_MAKE_INT(1)) return numer;
1033 if (denom == SG_MAKE_INT(0)) {
1034 int s = Sg_Sign(numer);
1035 if (s > 0) return SG_POSITIVE_INFINITY;
1036 if (s < 0) return SG_NEGATIVE_INFINITY;
1037 return SG_NAN;
1038 }
1039
1040 common = Sg_Gcd(numer, denom);
1041 if (SG_MAKE_INT(1) == common) {
1042 if (negated) {
1043 return make_rational(numer, denom);
1044 } else {
1045 return rational;
1046 }
1047 } else {
1048 numer = Sg_Quotient(numer, common, NULL);
1049 denom = Sg_Quotient(denom, common, NULL);
1050 if (SG_EQ(denom, SG_MAKE_INT(1))) {
1051 return numer;
1052 } else {
1053 return make_rational(numer, denom);
1054 }
1055 }
1056 }
1057
1058 #define set_nume_deno(r, n, d) \
1059 do { \
1060 if (SG_RATIONALP(r)) { \
1061 (n) = SG_RATIONAL(r)->numerator; \
1062 (d) = SG_RATIONAL(r)->denominator; \
1063 } else { \
1064 (n) = (r); \
1065 (d) = SG_MAKE_INT(1); \
1066 } \
1067 } while (0)
1068
1069
Sg_RationalAddSub(SgObject x,SgObject y,int subtract)1070 SgObject Sg_RationalAddSub(SgObject x, SgObject y, int subtract)
1071 {
1072 SgObject nx, dx, ny, dy;
1073 SgObject gcd, fx, fy, nr, dr;
1074
1075 set_nume_deno(x, nx, dx);
1076 set_nume_deno(y, ny, dy);
1077
1078 if (Sg_NumEq(dx, dy)) {
1079 dr = dx;
1080 goto finish;
1081 }
1082
1083 if (SG_MAKE_INT(1) == dx || SG_MAKE_INT(1) == dy) gcd = SG_MAKE_INT(1);
1084 else gcd = Sg_Gcd(dx, dy);
1085 if (Sg_NumEq(dx, gcd)) {
1086 /* only factor x */
1087 nx = Sg_Mul(Sg_Quotient(dy, dx, NULL), nx);
1088 dr = dy;
1089 goto finish;
1090 }
1091 if (Sg_NumEq(dy, gcd)) {
1092 /* only factor y */
1093 ny = Sg_Mul(Sg_Quotient(dx, dy, NULL), ny);
1094 dr = dx;
1095 goto finish;
1096 }
1097 /* general case */
1098 fx = Sg_Quotient(dx, gcd, NULL);
1099 fy = Sg_Quotient(dy, gcd, NULL);
1100 nx = Sg_Mul(nx, fy);
1101 ny = Sg_Mul(ny, fx);
1102 dr = Sg_Mul(dx, fy);
1103 finish:
1104 nr = (subtract ? Sg_Sub(nx, ny) : Sg_Add(nx, ny));
1105 return Sg_MakeRational(nr, dr);
1106 }
1107
Sg_RationalMulDiv(SgObject x,SgObject y,int divide)1108 SgObject Sg_RationalMulDiv(SgObject x, SgObject y, int divide)
1109 {
1110 SgObject nx, dx, ny, dy;
1111
1112 set_nume_deno(x, nx, dx);
1113 set_nume_deno(y, ny, dy);
1114
1115 if (divide) {
1116 /* swap */
1117 SgObject t = ny; ny = dy; dy = t;
1118 }
1119 return Sg_MakeRational(Sg_Mul(nx, ny),
1120 Sg_Mul(dx, dy));
1121 }
1122
1123 #define Sg_RationalAdd(x, y) Sg_RationalAddSub(x, y, FALSE)
1124 #define Sg_RationalSub(x, y) Sg_RationalAddSub(x, y, TRUE)
1125 #define Sg_RationalMul(x, y) Sg_RationalMulDiv(x, y, FALSE)
1126 #define Sg_RationalDiv(x, y) Sg_RationalMulDiv(x, y, TRUE)
1127
1128
make_flonum(double d)1129 static SgFlonum* make_flonum(double d)
1130 {
1131 SgFlonum *f;
1132 #ifdef USE_IMMEDIATE_FLONUM
1133 SgIFlonum ifl;
1134 /* check if flonum can be immediate */
1135 # if SIZEOF_VOIDP == 8
1136 ifl.f = d;
1137 if ((ifl.i & SG_IFLONUM_MASK) == 0) {
1138 ifl.i |= SG_IFLONUM_TAG;
1139 return SG_OBJ(ifl.i);
1140 }
1141 # else
1142 if (FLT_MIN <= d && d <= FLT_MAX) {
1143 ifl.f = (float)d;
1144 /* for MSVC we need to make scope... */
1145 {
1146 /* To keep calculation better
1147 TODO:
1148 this actually does not allow to make most of flonums immediate
1149 value except the numbers which have its fraction multiple of 5.
1150 So, this makes only gambit benchmarks sumfp and fibfp a bit faster.
1151
1152 IEEE single precision
1153 s eeeeeeee fffffffffffffffffffffff
1154 s = sign
1155 e = exponent ( 8 bits)
1156 f = fraction (23 bits)
1157
1158 */
1159 int e = ifl.i >> 23;
1160 e &= 0xFF; /* drop sign bit */
1161 if ((ifl.i & SG_IFLONUM_MASK) == 0 &&
1162 /* To make as less error as possible, we allow
1163 only exponent less than 15 or all 1. */
1164 (e == 0xFF || e <= 0xE)) {
1165 ifl.i |= SG_IFLONUM_TAG;
1166 return SG_OBJ(ifl.i);
1167 }
1168 }
1169 }
1170 # endif /* SIZEOF_VOIDP == 8 */
1171 #endif /* USE_IMMEDIATE_FLONUM */
1172 f = SG_NEW_ATOMIC(SgFlonum);
1173 SG_SET_CLASS(f, SG_CLASS_REAL);
1174 f->value = d;
1175 return SG_OBJ(f);
1176 }
1177
Sg_MakeFlonum(double d)1178 SgObject Sg_MakeFlonum(double d)
1179 {
1180 if (d == 0.0) {
1181 union { double f64; int64_t i64; } datum;
1182 datum.f64 = d;
1183 if (datum.i64 < 0) {
1184 return SG_FL_NEGATIVE_ZERO;
1185 } else {
1186 return SG_FL_POSITIVE_ZERO;
1187 }
1188 } else if (isnan(d)) return SG_NAN;
1189 else return make_flonum(d);
1190 }
1191
1192 #ifdef USE_IMMEDIATE_FLONUM
Sg_FlonumValue(SgObject obj)1193 double Sg_FlonumValue(SgObject obj)
1194 {
1195 /* MSVC does not allow me to cast */
1196 #ifndef __GNUC__
1197 if (SG_IFLONUMP(obj)) {
1198 void *p = (void *)((uintptr_t)obj&~SG_IFLONUM_MASK);
1199 return (double)((SgIFlonum *)&p)->f;
1200 } else {
1201 return ((SgFlonum *)(obj))->value;
1202 }
1203 #else
1204 return SG_FLONUM_VALUE(obj);
1205 #endif
1206 }
1207 #endif /* USE_IMMEDIATE_FLONUM */
1208
make_complex(SgObject real,SgObject imag)1209 static inline SgObject make_complex(SgObject real, SgObject imag)
1210 {
1211 SgComplex *c;
1212 ASSERT(!SG_COMPLEXP(real));
1213 ASSERT(!SG_COMPLEXP(imag));
1214 c = SG_NEW(SgComplex);
1215 SG_SET_CLASS(c, SG_CLASS_COMPLEX);
1216 c->real = real;
1217 c->imag = imag;
1218 return SG_OBJ(c);
1219 }
1220
Sg_MakeComplex(SgObject real,SgObject imag)1221 SgObject Sg_MakeComplex(SgObject real, SgObject imag)
1222 {
1223 if (!SG_FLONUMP(imag) && Sg_ZeroP(imag)) return real;
1224 else if (SG_FLONUMP(real) || SG_FLONUMP(imag)) {
1225 return make_complex(Sg_Inexact(real), Sg_Inexact(imag));
1226 }
1227 else
1228 return make_complex(real, imag);
1229 }
1230
Sg_MakeComplexPolar(SgObject magnitude,SgObject angle)1231 SgObject Sg_MakeComplexPolar(SgObject magnitude, SgObject angle)
1232 {
1233 double r, a;
1234 if (angle == SG_MAKE_INT(0)) return magnitude;
1235 r = Sg_GetDouble(magnitude);
1236 a = Sg_GetDouble(angle);
1237 return Sg_MakeComplex(Sg_MakeFlonum(r * cos(a)), Sg_MakeFlonum(r * sin(a)));
1238 }
1239
Sg_GetDouble(SgObject obj)1240 double Sg_GetDouble(SgObject obj)
1241 {
1242 if (SG_FLONUMP(obj)) return SG_FLONUM_VALUE(obj);
1243 else if (SG_INTP(obj)) return (double)SG_INT_VALUE(obj);
1244 else if (SG_BIGNUMP(obj)) return Sg_BignumToDouble(SG_BIGNUM(obj));
1245 else if (SG_RATIONAL(obj)) return Sg_RationalToDouble(obj);
1246 else if (SG_COMPLEXP(obj)) {
1247 SgComplex *c = SG_COMPLEX(obj);
1248 if (Sg_ZeroP(c->imag)) return Sg_GetDouble(c->real);
1249 else return 0.0;
1250 }
1251 else return 0.0; /* should this be error? */
1252 }
1253
Sg_DecodeFlonum(double d,int * exp,int * sign)1254 SgObject Sg_DecodeFlonum(double d, int *exp, int *sign)
1255 {
1256 SgObject f;
1257 int exp0, sign0;
1258 int64_t mant = decode_double(d, &exp0, &sign0);
1259 f = int64_to_integer(mant);
1260 *exp = exp0;
1261 *sign = sign0;
1262 return f;
1263 }
1264
Sg_RationalToDouble(SgRational * obj)1265 double Sg_RationalToDouble(SgRational *obj)
1266 {
1267 const int BITSIZE_TH = 96;
1268 double nume = Sg_GetDouble(obj->numerator);
1269 double deno = Sg_GetDouble(obj->denominator);
1270 if (isinf(nume) || isinf(deno)) {
1271 if (isinf(nume) && isinf(deno)) {
1272 long nume_bitsize = Sg_BignumBitSize(obj->numerator);
1273 long deno_bitsize = Sg_BignumBitSize(obj->denominator);
1274 long shift = (nume_bitsize > deno_bitsize)
1275 ? nume_bitsize - BITSIZE_TH : deno_bitsize - BITSIZE_TH;
1276 if (shift < 1) shift = 1;
1277 nume = Sg_GetDouble(Sg_BignumShiftRight(obj->numerator, shift));
1278 deno = Sg_GetDouble(Sg_BignumShiftRight(obj->denominator, shift));
1279 } else if (isinf(deno)) {
1280 long deno_bitsize = Sg_BignumBitSize(obj->denominator);
1281 long shift = deno_bitsize - BITSIZE_TH;
1282 if (shift < 1) shift = 1;
1283 nume = ldexp(nume, (int)-shift);
1284 deno = Sg_GetDouble(Sg_BignumShiftRight(obj->denominator, shift));
1285 } else {
1286 long nume_bitsize = Sg_BignumBitSize(obj->numerator);
1287 long shift = nume_bitsize - BITSIZE_TH;
1288 if (shift < 1) shift = 1;
1289 nume = Sg_GetDouble(Sg_BignumShiftRight(obj->numerator, shift));
1290 deno = ldexp(deno, (int)-shift);
1291 }
1292 }
1293 return nume / deno;
1294 }
1295
Sg_Numerator(SgObject x)1296 SgObject Sg_Numerator(SgObject x)
1297 {
1298 int inexact = FALSE;
1299 SgObject obj;
1300 if (!SG_NUMBERP(x)) wte(SG_INTERN("numerator"), "number", x);
1301 if (SG_FLONUMP(x)) {
1302 double d = SG_FLONUM_VALUE(x);
1303 if (d == 0.0) return x;
1304 if (isinf(d) || isnan(d)) return x;
1305 inexact = TRUE;
1306 }
1307 obj = Sg_Exact(x);
1308 if (SG_RATIONALP(obj)) {
1309 if (inexact) return Sg_Inexact(SG_RATIONAL(obj)->numerator);
1310 return SG_RATIONAL(obj)->numerator;
1311 }
1312 return inexact ? Sg_Inexact(obj) : obj;
1313 }
1314
Sg_Denominator(SgObject x)1315 SgObject Sg_Denominator(SgObject x)
1316 {
1317 int inexact = FALSE;
1318 SgObject obj;
1319 if (!SG_NUMBERP(x)) wte(SG_INTERN("denominator"), "number", x);
1320 if (SG_FLONUMP(x)) {
1321 double d = SG_FLONUM_VALUE(x);
1322 if (isinf(d)) return Sg_MakeFlonum(1.0);
1323 if (isnan(d)) return SG_NAN;
1324 inexact = TRUE;
1325 }
1326 obj = Sg_Exact(x);
1327 if (SG_RATIONALP(obj)) {
1328 if (inexact) return Sg_Inexact(SG_RATIONAL(obj)->denominator);
1329 return SG_RATIONAL(obj)->denominator;
1330 }
1331 return inexact ? Sg_MakeFlonum(1.0) : SG_MAKE_INT(1);
1332 }
1333
rationalize_rec(SgObject bottom,SgObject top)1334 static inline SgObject rationalize_rec(SgObject bottom, SgObject top)
1335 {
1336 if (Sg_NumCmp(bottom, top) == 0) return bottom;
1337 else {
1338 SgObject x = Sg_Round(bottom, SG_ROUND_CEIL);
1339 if (Sg_NumCmp(x, top) < 0) return x;
1340 else {
1341 SgObject one = SG_MAKE_INT(1);
1342 SgObject a = Sg_Sub(x, one);
1343 return Sg_Add(a, Sg_Div(one,
1344 rationalize_rec(Sg_Div(one,
1345 Sg_Sub(top, a)),
1346 Sg_Div(one,
1347 Sg_Sub(bottom, a)))));
1348 }
1349 }
1350 }
1351
Sg_Rationalize(SgObject x,SgObject e)1352 SgObject Sg_Rationalize(SgObject x, SgObject e)
1353 {
1354 if (Sg_InfiniteP(e)) {
1355 if (Sg_InfiniteP(x)) return SG_NAN;
1356 else return Sg_MakeFlonum(0.0);
1357 } else if (Sg_ZeroP(x)) return x;
1358 else if (Sg_NumCmp(x, e) == 0) return Sg_Sub(x, e);
1359 else if (Sg_NegativeP(x)) return Sg_Negate(Sg_Rationalize(Sg_Negate(x), e));
1360 else {
1361 SgObject bottom, top;
1362 e = Sg_Abs(e);
1363 bottom = Sg_Sub(x, e);
1364 top = Sg_Add(x, e);
1365 return rationalize_rec(bottom, top);
1366 }
1367 }
1368
Sg_StringToNumber(SgString * str,int radix,int strict)1369 SgObject Sg_StringToNumber(SgString *str, int radix, int strict)
1370 {
1371 return read_number(str->value, str->size, radix, strict);
1372 }
1373
Sg_ZeroP(SgObject obj)1374 int Sg_ZeroP(SgObject obj)
1375 {
1376 if (SG_INTP(obj)) return obj == SG_MAKE_INT(0);
1377 if (SG_FLONUMP(obj)) return SG_FLONUM_VALUE(obj) == 0.0;
1378 if (SG_BIGNUMP(obj)) {
1379 /* ASSERT(SG_BIGNUM_GET_SIGN(obj) != 0); */
1380 return SG_BIGNUM_GET_SIGN(obj) == 0;
1381 }
1382 if (SG_RATIONALP(obj)) {
1383 ASSERT(Sg_ZeroP(SG_RATIONAL(obj)->numerator) == FALSE);
1384 return FALSE;
1385 }
1386 if (SG_COMPLEXP(obj)) {
1387 SgComplex *c = SG_COMPLEX(obj);
1388 return Sg_ZeroP(c->real) && Sg_ZeroP(c->imag);
1389 }
1390 wte(SG_INTERN("zero?"), "number", obj);
1391 return -1; /* dummy */
1392 }
1393
Sg_IntegerP(SgObject obj)1394 int Sg_IntegerP(SgObject obj)
1395 {
1396 if (SG_INTP(obj) || SG_BIGNUMP(obj)) return TRUE;
1397 if (SG_RATIONALP(obj)) return FALSE; /* normalized ratnum never be integer */
1398 if (SG_FLONUMP(obj)) {
1399 double d = SG_FLONUM_VALUE(obj);
1400 double f, i;
1401 if (Sg_InfiniteP(obj)) return FALSE;
1402 if ((f = modf(d, &i)) == 0.0) return TRUE;
1403 return FALSE;
1404 }
1405 if (SG_COMPLEXP(obj)) return FALSE;
1406 return FALSE;
1407 }
1408
Sg_OddP(SgObject obj)1409 int Sg_OddP(SgObject obj)
1410 {
1411 if (SG_INTP(obj)) {
1412 return (SG_INT_VALUE(obj) & 1);
1413 }
1414 if (SG_BIGNUMP(obj)) {
1415 return (SG_BIGNUM(obj)->elements[0] & 1);
1416 }
1417 if (SG_FLONUMP(obj) && Sg_IntegerP(obj)) {
1418 return (fmod(SG_FLONUM_VALUE(obj), 2.0) != 0.0);
1419 }
1420 wte(SG_INTERN("odd?"), "integer", obj);
1421 return FALSE; /* dummy */
1422 }
1423
Sg_FiniteP(SgObject obj)1424 int Sg_FiniteP(SgObject obj)
1425 {
1426 return !Sg_InfiniteP(obj) && !Sg_NanP(obj);
1427 }
1428
Sg_InfiniteP(SgObject obj)1429 int Sg_InfiniteP(SgObject obj)
1430 {
1431 if (SG_FLONUMP(obj)) {
1432 double v = SG_FLONUM_VALUE(obj);
1433 return isinf(v);
1434 } else if (SG_COMPLEXP(obj)) {
1435 SgObject r = SG_COMPLEX(obj)->real;
1436 SgObject i = SG_COMPLEX(obj)->imag;
1437 return Sg_InfiniteP(r) || Sg_InfiniteP(i);
1438 } else if (!SG_NUMBERP(obj)) {
1439 wte(SG_INTERN("infinite?"), "number", obj);
1440 }
1441 return FALSE;
1442 }
1443
Sg_NanP(SgObject obj)1444 int Sg_NanP(SgObject obj)
1445 {
1446 if (SG_FLONUMP(obj)) {
1447 double v = SG_FLONUM_VALUE(obj);
1448 return isnan(v);
1449 } else if (SG_COMPLEXP(obj)) {
1450 SgObject r = SG_COMPLEX(obj)->real;
1451 SgObject i = SG_COMPLEX(obj)->imag;
1452 return Sg_NanP(r) || Sg_NanP(i);
1453 } else if (!SG_NUMBERP(obj)) {
1454 wte(SG_INTERN("nan?"), "number", obj);
1455 }
1456 return FALSE;
1457 }
1458
1459 /* TODO: the name is conflicted */
Sg_RationalP(SgObject n)1460 int Sg_RationalP(SgObject n)
1461 {
1462 if (SG_EXACT_INTP(n) || SG_RATIONALP(n)) return TRUE;
1463 if (SG_FLONUMP(n)) {
1464 if (Sg_InfiniteP(n)) return FALSE;
1465 if (Sg_NanP(n)) return FALSE;
1466 return TRUE;
1467 }
1468 return FALSE;
1469 }
1470
Sg_RealValuedP(SgObject n)1471 int Sg_RealValuedP(SgObject n)
1472 {
1473 if (SG_REALP(n)) return TRUE;
1474 if (SG_COMPLEXP(n)) {
1475 return Sg_ZeroP(SG_COMPLEX(n)->imag);
1476 }
1477 return FALSE;
1478 }
1479
Sg_RationalValuedP(SgObject n)1480 int Sg_RationalValuedP(SgObject n)
1481 {
1482 if (SG_RATIONALP(n)) return TRUE;
1483 if (SG_EXACT_INTP(n)) return TRUE;
1484 if (SG_FLONUMP(n)) {
1485 if (Sg_InfiniteP(n)) return FALSE;
1486 if (Sg_NanP(n)) return FALSE;
1487 return TRUE;
1488 }
1489 if (SG_COMPLEXP(n)) {
1490 return Sg_ZeroP(SG_COMPLEX(n)->imag) &&
1491 Sg_RationalP(SG_COMPLEX(n)->real);
1492 }
1493 return FALSE;
1494 }
1495
Sg_IntegerValuedP(SgObject n)1496 int Sg_IntegerValuedP(SgObject n)
1497 {
1498 if (Sg_IntegerP(n)) return TRUE;
1499 if (SG_COMPLEXP(n)) {
1500 return Sg_ZeroP(SG_COMPLEX(n)->imag) &&
1501 Sg_IntegerValuedP(SG_COMPLEX(n)->real);
1502 }
1503 return FALSE;
1504 }
1505
Sg_Negate(SgObject obj)1506 SgObject Sg_Negate(SgObject obj)
1507 {
1508 if (SG_INTP(obj)) {
1509 long n = SG_INT_VALUE(obj);
1510 if (n == SG_INT_MIN) return Sg_MakeInteger(-n);
1511 return SG_MAKE_INT(-n);
1512 }
1513 if (SG_FLONUMP(obj)) {
1514 return Sg_MakeFlonum(-SG_FLONUM_VALUE(obj));
1515 }
1516 if (SG_BIGNUMP(obj)) {
1517 SgBignum *b = Sg_BignumCopy(obj);
1518 SG_BIGNUM_SET_SIGN(b, -SG_BIGNUM_GET_SIGN(obj));
1519 return Sg_NormalizeBignum(b);
1520 }
1521 if (SG_RATIONALP(obj)) {
1522 SgRational *r = SG_RATIONAL(obj);
1523 return Sg_MakeRational(Sg_Negate(r->numerator), r->denominator);
1524 }
1525 if (SG_COMPLEXP(obj)) {
1526 SgComplex *c = SG_COMPLEX(obj);
1527 return Sg_MakeComplex(Sg_Negate(c->real), Sg_Negate(c->imag));
1528 }
1529 wte(SG_INTERN("negate"), "number", obj);
1530 return SG_UNDEF; /* dummy */
1531 }
1532
Sg_NegativeP(SgObject obj)1533 int Sg_NegativeP(SgObject obj)
1534 {
1535 if (SG_INTP(obj)) return SG_INT_VALUE(obj) < 0;
1536 if (SG_BIGNUMP(obj)) return SG_BIGNUM_GET_SIGN(obj) < 0;
1537 if (SG_FLONUMP(obj)) return SG_FLONUM_VALUE(obj) < 0.0;
1538 if (SG_RATIONALP(obj)) return Sg_NegativeP(SG_RATIONAL(obj)->numerator);
1539 if (SG_COMPLEXP(obj)) return Sg_NegativeP(SG_COMPLEX(obj)->real);
1540 wte(SG_INTERN("negative?"), "number", obj);
1541 return FALSE; /* dummy */
1542 }
1543
Sg_PositiveP(SgObject obj)1544 int Sg_PositiveP(SgObject obj)
1545 {
1546 if (SG_INTP(obj)) return SG_INT_VALUE(obj) > 0;
1547 if (SG_BIGNUMP(obj)) return SG_BIGNUM_GET_SIGN(obj) > 0;
1548 if (SG_FLONUMP(obj)) {
1549 #if __WATCOMC__
1550 /* on Watcom, +nan.0 is bigger than 0 */
1551 if (isnan(SG_FLONUM_VALUE(obj))) return FALSE;
1552 else return SG_FLONUM_VALUE(obj) > 0.0;
1553 #else
1554 return SG_FLONUM_VALUE(obj) > 0.0;
1555 #endif
1556 }
1557 if (SG_RATIONALP(obj)) return Sg_PositiveP(SG_RATIONAL(obj)->numerator);
1558 if (SG_COMPLEXP(obj)) return Sg_PositiveP(SG_COMPLEX(obj)->real);
1559 wte(SG_INTERN("positive?"), "number", obj);
1560 return FALSE; /* dummy */
1561 }
1562
Sg_Exact(SgObject obj)1563 SgObject Sg_Exact(SgObject obj)
1564 {
1565 if (SG_FLONUMP(obj)) {
1566 double d = SG_FLONUM_VALUE(obj);
1567 double f, i;
1568 /* inf.0 or nan.0 doesn't have exact value. */
1569
1570 if (isinf(d) || isnan(d)) {
1571 Sg_AssertionViolation(SG_INTERN("exact"),
1572 Sg_Sprintf(UC("no exact representation for %S"), obj),
1573 SG_LIST1(obj));
1574 }
1575 #if defined(__WATCOMC__) || defined(__FreeBSD__) || defined(__OpenBSD__) || defined(__DragonFly__)
1576 /* Yes, on Watcom/FreeBSD if +inf.0 or +nan.0 is passed to modf,
1577 it returns +nan.0. Sucks!!! */
1578 if ((f = modf(d, &i)) == 0.0 ||
1579 isinf(d) || isnan(d)) {
1580 if (isnan(d)) {
1581 return SG_MAKE_INT(0);
1582 } else if (d < SG_INT_MIN || d > SG_INT_MAX) {
1583 obj = Sg_MakeBignumFromDouble(d);
1584 } else {
1585 obj = SG_MAKE_INT((long)d);
1586 }
1587 }
1588 #else
1589 if ((f = modf(d, &i)) == 0.0) {
1590 if (d < SG_INT_MIN || d > SG_INT_MAX) {
1591 obj = Sg_MakeBignumFromDouble(d);
1592 } else {
1593 obj = SG_MAKE_INT((long)d);
1594 }
1595 }
1596 #endif
1597 else {
1598 SgObject m;
1599 int exp, sign;
1600 m = Sg_DecodeFlonum(d, &exp, &sign);
1601 ASSERT(exp < 0); /* exp >= case should be handled above */
1602 obj = Sg_Div(m, Sg_Ash(SG_MAKE_INT(1), -exp));
1603 if (sign < 0) obj = Sg_Negate(obj);
1604 }
1605 return obj;
1606 } else if (SG_COMPLEXP(obj)) {
1607 SgComplex *c = SG_COMPLEX(obj);
1608 if (SG_FLONUMP(c->real) || SG_FLONUMP(c->imag)) {
1609 return oprtr_norm_complex(Sg_Exact(c->real), Sg_Exact(c->imag));
1610 }
1611 return obj;
1612 }
1613 if (SG_INTP(obj) || SG_BIGNUMP(obj) || SG_RATIONALP(obj)) return obj;
1614 wte(SG_INTERN("exact"), "number", obj);
1615 return SG_UNDEF; /* dummy */
1616 }
1617
Sg_Inexact(SgObject obj)1618 SgObject Sg_Inexact(SgObject obj)
1619 {
1620 if (SG_INTP(obj)) {
1621 double z = (double)SG_INT_VALUE(obj);
1622 return Sg_MakeFlonum(z);
1623 } else if (SG_BIGNUMP(obj)) {
1624 double z = Sg_BignumToDouble(SG_BIGNUM(obj));
1625 return Sg_MakeFlonum(z);
1626 } else if (SG_RATIONALP(obj)) {
1627 double z = Sg_GetDouble(obj);
1628 return Sg_MakeFlonum(z);
1629 } else if (SG_FLONUMP(obj)) {
1630 return obj;
1631 } else if (SG_COMPLEXP(obj)) {
1632 SgComplex *c = SG_COMPLEX(obj);
1633 if (SG_FLONUMP(c->real) & SG_FLONUMP(c->imag)) return obj;
1634 return Sg_MakeComplex(Sg_Inexact(c->real), Sg_Inexact(c->imag));
1635 } else {
1636 wte(SG_INTERN("inexact"), "number", obj);
1637 }
1638 return SG_UNDEF; /* dummy */
1639 }
1640
1641
Sg_ExactP(SgObject obj)1642 int Sg_ExactP(SgObject obj)
1643 {
1644 if (SG_INTP(obj) || SG_BIGNUMP(obj) || SG_RATIONALP(obj)) return TRUE;
1645 if (SG_FLONUMP(obj)) return FALSE;
1646 if (SG_COMPLEXP(obj)) {
1647 SgComplex *c = SG_COMPLEX(obj);
1648 return Sg_ExactP(c->real) && Sg_ExactP(c->imag);
1649 }
1650 wte(SG_INTERN("exact?"), "number", obj);
1651 return FALSE;
1652 }
1653
Sg_InexactP(SgObject obj)1654 int Sg_InexactP(SgObject obj)
1655 {
1656 if (SG_INTP(obj) || SG_BIGNUMP(obj) || SG_RATIONALP(obj)) return FALSE;
1657 if (SG_FLONUMP(obj)) return TRUE;
1658 if (SG_COMPLEXP(obj)) {
1659 SgComplex *c = SG_COMPLEX(obj);
1660 /* TODO is it correct? */
1661 return Sg_InexactP(c->real) || Sg_InexactP(c->imag);
1662 }
1663 wte(SG_INTERN("inexact?"), "number", obj);
1664 return FALSE;
1665 }
1666
Sg_Inverse(SgObject obj)1667 SgObject Sg_Inverse(SgObject obj)
1668 {
1669 if (SG_INTP(obj)) {
1670 if (SG_INT_VALUE(obj) == 0) wte(SG_INTERN("inverse"), "non 0 number", obj);
1671 if (SG_INT_VALUE(obj) > 0) {
1672 if (SG_INT_VALUE(obj) == 1) return obj;
1673 return Sg_MakeRational(SG_MAKE_INT(1), obj);
1674 }
1675 if (obj == SG_MAKE_INT(-1)) return obj;
1676 return Sg_MakeRational(SG_MAKE_INT(-1), Sg_Negate(obj));
1677 }
1678 if (SG_FLONUMP(obj)) return Sg_MakeFlonum(1.0 / SG_FLONUM_VALUE(obj));
1679 if (SG_BIGNUMP(obj)) {
1680 if (SG_BIGNUM_GET_SIGN(obj) == 0) wte(SG_INTERN("inverse"), "non 0 number", obj);
1681 if (SG_BIGNUM_GET_SIGN(obj) > 0) {
1682 return Sg_MakeRational(SG_MAKE_INT(1), obj);
1683 }
1684 return Sg_MakeRational(SG_MAKE_INT(-1), Sg_Negate(obj));
1685 }
1686 if (SG_RATIONALP(obj)) {
1687 if (!Sg_NegativeP(SG_RATIONAL(obj)->numerator)) {
1688 if (SG_RATIONAL(obj)->numerator == SG_MAKE_INT(1)) {
1689 return oprtr_norm_integer(SG_RATIONAL(obj)->denominator);
1690 }
1691 return Sg_MakeRational(SG_RATIONAL(obj)->denominator, SG_RATIONAL(obj)->numerator);
1692 }
1693 if (SG_RATIONAL(obj)->numerator == SG_MAKE_INT(-1)) {
1694 return Sg_Negate(SG_RATIONAL(obj)->denominator);
1695 }
1696 return Sg_MakeRational(Sg_Negate(SG_RATIONAL(obj)->denominator),
1697 Sg_Negate(SG_RATIONAL(obj)->numerator));
1698 }
1699 if (SG_COMPLEXP(obj)) return Sg_Div(SG_MAKE_INT(1), obj);
1700 wte(SG_INTERN("inverse"), "number", obj);
1701 return SG_UNDEF; /* dummy */
1702 }
1703
integer_length_rec(SgObject n)1704 static inline long integer_length_rec(SgObject n)
1705 {
1706 long n2;
1707 SgObject n3;
1708 if (SG_INTP(n)) {
1709 switch (SG_INT_VALUE(n)) {
1710 case 0: case -1:
1711 return 0;
1712 case 1: case -2:
1713 return 1;
1714 case 2: case 3: case -3: case -4:
1715 return 2;
1716 case 4: case 5: case 6: case 7:
1717 case -5: case -6: case -7: case -8:
1718 return 3;
1719 }
1720 /* fall through */
1721 }
1722 n2 = integer_length_rec(Sg_Ash(n, -4));
1723 n3 = Sg_Add(SG_MAKE_INT(4),
1724 SG_MAKE_INT(n2));
1725 ASSERT(SG_INTP(n3));
1726 return SG_INT_VALUE(n3);
1727 }
1728
Sg_IntegerLength(SgObject n)1729 long Sg_IntegerLength(SgObject n)
1730 {
1731 if (!Sg_IntegerP(n)) wte(SG_INTERN("integer-length"), "integer", n);
1732 return integer_length_rec(n);
1733 }
1734
Sg_Ash(SgObject x,long count)1735 SgObject Sg_Ash(SgObject x, long count)
1736 {
1737 if (SG_INTP(x)) {
1738 long ix = SG_INT_VALUE(x);
1739 if (x == 0) {
1740 return SG_MAKE_INT(0);
1741 } else if (count <= -(SIZEOF_LONG * 8)) {
1742 ix = (ix < 0) ? -1 : 0;
1743 return Sg_MakeInteger(ix);
1744 } else if (count < 0) {
1745 if (ix < 0) {
1746 ix = ~((~ix) >> (-count));
1747 } else {
1748 ix >>= -count;
1749 }
1750 return Sg_MakeInteger(ix);
1751 } else if (count < SG_INT_SIZE) {
1752 if (ix < 0) {
1753 if (-ix < (SG_INT_MAX >> count)) {
1754 ix <<= count;
1755 return Sg_MakeInteger(ix);
1756 }
1757 } else {
1758 if (ix < (SG_INT_MAX >> count)) {
1759 ix <<= count;
1760 return Sg_MakeInteger(ix);
1761 }
1762 }
1763 }
1764 return Sg_BignumAshSI(ix, count);
1765 } else if (SG_BIGNUMP(x)) {
1766 return Sg_BignumAsh(SG_BIGNUM(x), count);
1767 }
1768 wte(SG_INTERN("bitwise-arithmetic-shift"), "exact integer", x);
1769 return SG_UNDEF; /* dummy */
1770 }
1771
Sg_LogNot(SgObject x)1772 SgObject Sg_LogNot(SgObject x)
1773 {
1774 if (!SG_EXACT_INTP(x)) wte(SG_INTERN("bitwise-not"), "exact integer", x);
1775 if (SG_INTP(x)) {
1776 /* this won't cause an overflow */
1777 return SG_MAKE_INT(~SG_INT_VALUE(x));
1778 } else {
1779 return Sg_Negate(Sg_BignumAddSI(SG_BIGNUM(x), 1));
1780 }
1781 }
1782
Sg_LogAnd(SgObject x,SgObject y)1783 SgObject Sg_LogAnd(SgObject x, SgObject y)
1784 {
1785 if (!SG_EXACT_INTP(x)) wte(SG_INTERN("bitwise-and"), "exact integer", x);
1786 if (!SG_EXACT_INTP(y)) wte(SG_INTERN("bitwise-and"), "exact integer", y);
1787 if (SG_INTP(x)) {
1788 if (SG_INTP(y)) {
1789 return SG_MAKE_INT(SG_INT_VALUE(x) & SG_INT_VALUE(y));
1790 } else if (SG_INT_VALUE(x) >= 0) {
1791 if (SG_BIGNUM_GET_SIGN(y) > 0) {
1792 return Sg_MakeInteger(SG_INT_VALUE(x) & SG_BIGNUM(y)->elements[0]);
1793 } else if (SG_BIGNUM_GET_SIGN(y) == 0) {
1794 return SG_MAKE_INT(0);
1795 }
1796 }
1797 /* x = Sg_MakeBignumFromSI(SG_INT_VALUE(x)); */
1798 return Sg_BignumLogAndSI(SG_BIGNUM(y), SG_INT_VALUE(x));
1799 } else if (SG_INTP(y)) {
1800 if (SG_INT_VALUE(y) >= 0) {
1801 if (SG_BIGNUM_GET_SIGN(x) > 0) {
1802 return Sg_MakeInteger(SG_INT_VALUE(y) & SG_BIGNUM(x)->elements[0]);
1803 } else if (SG_BIGNUM_GET_SIGN(x) == 0) {
1804 return SG_MAKE_INT(0);
1805 }
1806 }
1807 /* y = Sg_MakeBignumFromSI(SG_INT_VALUE(y)); */
1808 return Sg_BignumLogAndSI(SG_BIGNUM(x), SG_INT_VALUE(y));
1809 }
1810 return Sg_BignumLogAnd(SG_BIGNUM(x), SG_BIGNUM(y));
1811 }
1812
Sg_LogIor(SgObject x,SgObject y)1813 SgObject Sg_LogIor(SgObject x, SgObject y)
1814 {
1815 if (!SG_EXACT_INTP(x)) wte(SG_INTERN("bitwise-ior"), "exact integer", x);
1816 if (!SG_EXACT_INTP(y)) wte(SG_INTERN("bitwise-ior"), "exact integer", y);
1817 if (SG_INTP(x)) {
1818 if (SG_INTP(y)) {
1819 return SG_MAKE_INT(SG_INT_VALUE(x) | SG_INT_VALUE(y));
1820 } else {
1821 /* x = Sg_MakeBignumFromSI(SG_INT_VALUE(x)); */
1822 return Sg_BignumLogIorSI(SG_BIGNUM(y), SG_INT_VALUE(x));
1823 }
1824 } else {
1825 if (SG_INTP(y)) {
1826 /* y = Sg_MakeBignumFromSI(SG_INT_VALUE(y)); */
1827 return Sg_BignumLogIorSI(SG_BIGNUM(x), SG_INT_VALUE(y));
1828 }
1829 }
1830 return Sg_BignumLogIor(SG_BIGNUM(x), SG_BIGNUM(y));
1831 }
1832
Sg_LogXor(SgObject x,SgObject y)1833 SgObject Sg_LogXor(SgObject x, SgObject y)
1834 {
1835 if (!SG_EXACT_INTP(x)) wte(SG_INTERN("bitwise-xor"), "exact integer", x);
1836 if (!SG_EXACT_INTP(y)) wte(SG_INTERN("bitwise-xor"), "exact integer", y);
1837 if (SG_INTP(x)) {
1838 if (SG_INTP(y)) {
1839 return SG_MAKE_INT(SG_INT_VALUE(x) ^ SG_INT_VALUE(y));
1840 } else {
1841 /* x = Sg_MakeBignumFromSI(SG_INT_VALUE(x)); */
1842 return Sg_BignumLogXorSI(y, SG_INT_VALUE(x));
1843 }
1844 } else {
1845 if (SG_INTP(y)) {
1846 /* y = Sg_MakeBignumFromSI(SG_INT_VALUE(y)); */
1847 return Sg_BignumLogXorSI(x, SG_INT_VALUE(y));
1848 }
1849 }
1850 return Sg_BignumLogXor(SG_BIGNUM(x), SG_BIGNUM(y));
1851 }
1852
Sg_BitCount(SgObject x)1853 long Sg_BitCount(SgObject x)
1854 {
1855 if (!SG_EXACT_INTP(x)) wte(SG_INTERN("bitwise-bit-count"), "exact integer", x);
1856 if (SG_INTP(x)) {
1857 long n = SG_INT_VALUE(x);
1858 if (n >= 0) {
1859 return nbits(n);
1860 } else {
1861 return ~nbits(~n);
1862 }
1863 } else {
1864 return Sg_BignumBitCount(SG_BIGNUM(x));
1865 }
1866 }
1867
Sg_BitSize(SgObject x)1868 long Sg_BitSize(SgObject x)
1869 {
1870 if (!SG_EXACT_INTP(x)) wte(SG_INTERN("bitwise-length"), "exact integer", x);
1871 if (SG_INTP(x)) {
1872 long n = SG_INT_VALUE(x), n2;
1873 if (n == 0) return 0;
1874 n2 = (n < 0) ? ~n : n;
1875 return WORD_BITS - nlz(n2);
1876 } else {
1877 if (SG_BIGNUM_GET_SIGN(x) > 0) return Sg_BignumBitSize(SG_BIGNUM(x));
1878 else return Sg_BitSize(Sg_LogNot(x));
1879 }
1880 }
1881
Sg_FirstBitSet(SgObject x)1882 long Sg_FirstBitSet(SgObject x)
1883 {
1884 if (!SG_EXACT_INTP(x)) {
1885 wte(SG_INTERN("bitwise-first-bit-set"), "exact integer", x);
1886 }
1887 if (SG_INTP(x)) {
1888 long n = SG_INT_VALUE(x);
1889 int bit;
1890 if (n == 0) return -1;
1891 bit = 0;
1892 bit += ntz(n);
1893 return bit;
1894 } else {
1895 return Sg_BignumFirstBitSet(SG_BIGNUM(x));
1896 }
1897 }
1898
Sg_BitSetP(SgObject x,long n)1899 int Sg_BitSetP(SgObject x, long n)
1900
1901 {
1902 if (!SG_EXACT_INTP(x)) {
1903 wte(SG_INTERN("bitwise-bit-set?"), "exact integer", x);
1904 }
1905 if (n < 0) {
1906 wte(SG_INTERN("bitwise-bit-set?"), "non negative integer", SG_MAKE_INT(n));
1907 }
1908
1909 if (SG_INTP(x)) {
1910 long ix = SG_INT_VALUE(x);
1911 if (n >= SG_INT_SIZE) {
1912 return (ix < 0);
1913 }
1914 return (ix>>n)&1;
1915 }
1916 return Sg_BignumBitSetP(SG_BIGNUM(x), n);
1917 }
1918
Sg_Add(SgObject x,SgObject y)1919 SgObject Sg_Add(SgObject x, SgObject y)
1920 {
1921 if (SG_INTP(x)) {
1922 if (SG_INTP(y)) {
1923 long r = SG_INT_VALUE(x) + SG_INT_VALUE(y);
1924 return Sg_MakeInteger(r);
1925 }
1926 else if (SG_BIGNUMP(y)) {
1927 if (x == SG_MAKE_INT(0)) return y;
1928 else return Sg_BignumAddSI(SG_BIGNUM(y), SG_INT_VALUE(x));
1929 }
1930 else if (SG_RATIONALP(y)) {
1931 if (x == SG_MAKE_INT(0)) return y;
1932 else return Sg_RationalAdd(SG_RATIONAL(y), x);
1933 }
1934 else if (SG_FLONUMP(y)) {
1935 return Sg_MakeFlonum((double)SG_INT_VALUE(x) + SG_FLONUM_VALUE(y));
1936 }
1937 else if (SG_COMPLEXP(y)) {
1938 if (x == SG_MAKE_INT(0)) return y;
1939 else return Sg_MakeComplex(Sg_Add(SG_COMPLEX(y)->real, x),
1940 SG_COMPLEX(y)->imag);
1941 }
1942 }
1943 else if (SG_FLONUMP(x)) {
1944 if (SG_INTP(y)) {
1945 if (y == SG_MAKE_INT(0)) return x;
1946 else {
1947 double z = SG_FLONUM_VALUE(x) + (double)SG_INT_VALUE(y);
1948 return Sg_MakeFlonum(z);
1949 }
1950 }
1951 else if (SG_BIGNUMP(y) || SG_RATIONALP(y)) {
1952 return Sg_MakeFlonum(SG_FLONUM_VALUE(x) + Sg_GetDouble(y));
1953 }
1954 else if (SG_FLONUMP(y)) {
1955 return Sg_MakeFlonum(SG_FLONUM_VALUE(x) + SG_FLONUM_VALUE(y));
1956 }
1957 else if (SG_COMPLEXP(y)) {
1958 return Sg_MakeComplex(Sg_Add(SG_COMPLEX(y)->real, x),
1959 SG_COMPLEX(y)->imag);
1960 }
1961 }
1962 else if (SG_BIGNUMP(x)) {
1963 if (SG_INTP(y)) {
1964 if (y == SG_MAKE_INT(0)) return x;
1965 else return Sg_BignumAddSI(SG_BIGNUM(x), SG_INT_VALUE(y));
1966 }
1967 else if (SG_BIGNUMP(y)) {
1968 return Sg_BignumAdd(SG_BIGNUM(x), SG_BIGNUM(y));
1969 }
1970 else if (SG_RATIONALP(y)) {
1971 return Sg_RationalAdd(x, y);
1972 }
1973 else if (SG_FLONUMP(y)) {
1974 return Sg_MakeFlonum(Sg_BignumToDouble(x) + SG_FLONUM_VALUE(y));
1975 }
1976 else if (SG_COMPLEXP(y)) {
1977 return Sg_MakeComplex(Sg_Add(SG_COMPLEX(y)->real, x), SG_COMPLEX(y)->imag);
1978 }
1979 }
1980 else if (SG_RATIONALP(x)) {
1981 if (SG_INTP(y)) {
1982 if (y == SG_MAKE_INT(0)) return x;
1983 else return Sg_RationalAdd(x, y);
1984 }
1985 else if (SG_BIGNUMP(y) || SG_RATIONALP(y)) {
1986 return Sg_RationalAdd(x, y);
1987 }
1988 else if (SG_FLONUMP(y)) {
1989 return Sg_MakeFlonum(Sg_GetDouble(x) + SG_FLONUM_VALUE(y));
1990 }
1991 else if (SG_COMPLEXP(y)) {
1992 return Sg_MakeComplex(Sg_Add(SG_COMPLEX(y)->real, x),
1993 SG_COMPLEX(y)->imag);
1994 }
1995 }
1996 else if (SG_COMPLEXP(x)) {
1997 if (SG_INTP(y)) {
1998 if (y == SG_MAKE_INT(0)) return x;
1999 else return Sg_MakeComplex(Sg_Add(SG_COMPLEX(x)->real, y),
2000 SG_COMPLEX(x)->imag);
2001 }
2002 else if (SG_BIGNUMP(y) || SG_RATIONALP(y)) {
2003 return Sg_MakeComplex(Sg_Add(SG_COMPLEX(x)->real, y),
2004 SG_COMPLEX(x)->imag);
2005
2006 }
2007 else if (SG_FLONUMP(y)) {
2008 return Sg_MakeComplex(Sg_Add(SG_COMPLEX(x)->real, y),
2009 Sg_Inexact(SG_COMPLEX(x)->imag));
2010 }
2011 else if (SG_COMPLEXP(y)) {
2012 SgObject real = Sg_Add(SG_COMPLEX(x)->real, SG_COMPLEX(y)->real);
2013 SgObject imag = Sg_Add(SG_COMPLEX(x)->imag, SG_COMPLEX(y)->imag);
2014 return oprtr_norm_complex(real, imag);
2015 }
2016 }
2017 wte(SG_INTERN("+"), "number", SG_LIST2(x, y));
2018 return SG_UNDEF; /* dummy */
2019 }
2020
Sg_Sub(SgObject x,SgObject y)2021 SgObject Sg_Sub(SgObject x, SgObject y)
2022 {
2023 if (SG_INTP(x)) {
2024 if (SG_INTP(y)) {
2025 long r = SG_INT_VALUE(x) - SG_INT_VALUE(y);
2026 return Sg_MakeInteger(r);
2027 }
2028 else if (SG_BIGNUMP(y)) {
2029 SgObject big = Sg_MakeBignumFromSI(SG_INT_VALUE(x));
2030 return Sg_BignumSub(SG_BIGNUM(big), SG_BIGNUM(y));
2031 }
2032 else if (SG_RATIONALP(y)) {
2033 if (x == SG_MAKE_INT(0)) return Sg_Negate(y);
2034 else return Sg_RationalSub(x, y);
2035 }
2036 else if (SG_FLONUMP(y)) {
2037 return Sg_MakeFlonum((double)SG_INT_VALUE(x) - SG_FLONUM_VALUE(y));
2038 }
2039 else if (SG_COMPLEXP(y)) {
2040 SgObject img = SG_COMPLEX(y)->imag;
2041 if (x == SG_MAKE_INT(0)) return Sg_Negate(y);
2042 /* imag part must be negate. */
2043 else return Sg_MakeComplex(Sg_Sub(x, SG_COMPLEX(y)->real),
2044 Sg_Negate(img));
2045 }
2046 }
2047 else if (SG_FLONUMP(x)) {
2048 if (SG_INTP(y)) {
2049 if (y == SG_MAKE_INT(0)) return x;
2050 else {
2051 double z = SG_FLONUM_VALUE(x) - (double)SG_INT_VALUE(y);
2052 return Sg_MakeFlonum(z);
2053 }
2054 }
2055 else if (SG_BIGNUMP(y) || SG_RATIONALP(y)) {
2056 return Sg_MakeFlonum(SG_FLONUM_VALUE(x) - Sg_GetDouble(y));
2057 }
2058 else if (SG_FLONUMP(y)) {
2059 return Sg_MakeFlonum(SG_FLONUM_VALUE(x) - SG_FLONUM_VALUE(y));
2060 }
2061 else if (SG_COMPLEXP(y)) {
2062 /* imag part must be negate. */
2063 SgObject img = SG_COMPLEX(y)->imag;
2064 return Sg_MakeComplex(Sg_Sub(x, SG_COMPLEX(y)->real), Sg_Negate(img));
2065 }
2066 }
2067 else if (SG_BIGNUMP(x)) {
2068 if (SG_INTP(y)) {
2069 if (y == SG_MAKE_INT(0)) return x;
2070 else return Sg_BignumSubSI(SG_BIGNUM(x), SG_INT_VALUE(y));
2071 }
2072 else if (SG_BIGNUMP(y)) {
2073 return Sg_BignumSub(SG_BIGNUM(x), SG_BIGNUM(y));
2074 }
2075 else if (SG_RATIONALP(y)) {
2076 return Sg_RationalSub(x, y);
2077 }
2078 else if (SG_FLONUMP(y)) {
2079 return Sg_MakeFlonum(Sg_BignumToDouble(x) - SG_FLONUM_VALUE(y));
2080 }
2081 else if (SG_COMPLEXP(y)) {
2082 /* imag part must be negate. */
2083 SgObject img = SG_COMPLEX(y)->imag;
2084 return Sg_MakeComplex(Sg_Sub(x, SG_COMPLEX(y)->real), Sg_Negate(img));
2085 }
2086 }
2087 else if (SG_RATIONALP(x)) {
2088 if (SG_INTP(y)) {
2089 if (y == SG_MAKE_INT(0)) return x;
2090 else return Sg_RationalSub(x, y);
2091 }
2092 else if (SG_BIGNUMP(y) || SG_RATIONALP(y)) {
2093 return Sg_RationalSub(x, y);
2094 }
2095 else if (SG_FLONUMP(y)) {
2096 return Sg_MakeFlonum(Sg_GetDouble(x) - SG_FLONUM_VALUE(y));
2097 }
2098 else if (SG_COMPLEXP(y)) {
2099 /* imag part must be negate. */
2100 SgObject img = SG_COMPLEX(y)->imag;
2101 return Sg_MakeComplex(Sg_Sub(x, SG_COMPLEX(y)->real), Sg_Negate(img));
2102 }
2103 }
2104 else if (SG_COMPLEXP(x)) {
2105 if (SG_INTP(y)) {
2106 if (y == SG_MAKE_INT(0)) return x;
2107 else return Sg_MakeComplex(Sg_Sub(SG_COMPLEX(x)->real, y),
2108 SG_COMPLEX(x)->imag);
2109 }
2110 else if (SG_BIGNUMP(y) || SG_RATIONALP(y)) {
2111 return Sg_MakeComplex(Sg_Sub(SG_COMPLEX(x)->real, y),
2112 SG_COMPLEX(x)->imag);
2113
2114 }
2115 else if (SG_FLONUMP(y)) {
2116 return Sg_MakeComplex(Sg_Sub(SG_COMPLEX(x)->real, y),
2117 Sg_Inexact(SG_COMPLEX(x)->imag));
2118 }
2119 else if (SG_COMPLEXP(y)) {
2120 SgObject real = Sg_Sub(SG_COMPLEX(x)->real, SG_COMPLEX(y)->real);
2121 SgObject imag = Sg_Sub(SG_COMPLEX(x)->imag, SG_COMPLEX(y)->imag);
2122 return oprtr_norm_complex(real, imag);
2123 }
2124 }
2125 wte(SG_INTERN("-"), "number", SG_LIST2(x, y));
2126 return SG_UNDEF; /* dummy */
2127 }
2128
Sg_Mul(SgObject x,SgObject y)2129 SgObject Sg_Mul(SgObject x, SgObject y)
2130 {
2131 if (SG_INTP(x)) {
2132 if (SG_INTP(y)) {
2133 long v0 = SG_INT_VALUE(x);
2134 long v1 = SG_INT_VALUE(y);
2135 udlong k = v0 * v1;
2136 if ((v1 != 0 && k / v1 != v0) || !(k >= SG_INT_MIN && k <= SG_INT_MAX)) {
2137 SgObject big = Sg_MakeBignumFromSI(v0);
2138 return Sg_BignumMulSI(SG_BIGNUM(big), v1);
2139 } else
2140 return Sg_MakeInteger((long)k);
2141 } else if (SG_BIGNUMP(y)) {
2142 if (x == SG_MAKE_INT(0)) return x;
2143 else if (x == SG_MAKE_INT(1)) return y;
2144 else return Sg_BignumMulSI(SG_BIGNUM(y), SG_INT_VALUE(x));
2145 } else if (SG_RATIONALP(y)) {
2146 if (x == SG_MAKE_INT(0)) return x;
2147 else if (x == SG_MAKE_INT(1)) return y;
2148 else return Sg_RationalMul(x, y);
2149 } else if (SG_FLONUMP(y)) {
2150 if (x == SG_MAKE_INT(0)) return x;
2151 else if (x == SG_MAKE_INT(1)) return y;
2152 else return Sg_MakeFlonum((double)SG_INT_VALUE(x) * SG_FLONUM_VALUE(y));
2153 } else if (SG_COMPLEXP(y)) {
2154 if (x == SG_MAKE_INT(0)) return x;
2155 else if (x == SG_MAKE_INT(1)) return y;
2156 else return Sg_MakeComplex(Sg_Mul(x, SG_COMPLEX(y)->real),
2157 Sg_Mul(x, SG_COMPLEX(y)->imag));
2158 }
2159 }
2160 else if (SG_FLONUMP(x)) {
2161 if (SG_INTP(y)) {
2162 if (x == SG_MAKE_INT(0)) return x;
2163 else if (x == SG_MAKE_INT(1)) return y;
2164 else {
2165 double z = SG_FLONUM_VALUE(x) * (double)SG_INT_VALUE(y);
2166 return Sg_MakeFlonum(z);
2167 }
2168 } else if (SG_BIGNUMP(y) || SG_RATIONALP(y)) {
2169 return Sg_MakeFlonum(SG_FLONUM_VALUE(x) * Sg_GetDouble(y));
2170 } else if (SG_FLONUMP(y)) {
2171 if (SG_FLONUM_VALUE(y) == 1.0) return x;
2172 return Sg_MakeFlonum(SG_FLONUM_VALUE(x) * SG_FLONUM_VALUE(y));
2173 } else if (SG_COMPLEXP(y)) {
2174 return Sg_MakeComplex(Sg_Mul(x, SG_COMPLEX(y)->real),
2175 Sg_Mul(x, SG_COMPLEX(y)->imag));
2176 }
2177 }
2178 else if (SG_BIGNUMP(x)) {
2179 if (SG_INTP(y)) {
2180 if (y == SG_MAKE_INT(0)) return y;
2181 return Sg_BignumMulSI(SG_BIGNUM(x), SG_INT_VALUE(y));
2182 } else if (SG_BIGNUMP(y)) {
2183 /* squaring is faster than multiplication so if given arguments
2184 are the same object, then do it with square. */
2185 if (SG_EQ(x, y)) return Sg_BignumSquare(x);
2186 return Sg_BignumMul(SG_BIGNUM(x), SG_BIGNUM(y));
2187 } else if (SG_RATIONALP(y)) {
2188 return Sg_RationalMul(x, y);
2189 } else if (SG_FLONUMP(y)) {
2190 return Sg_MakeFlonum(Sg_BignumToDouble(x) * SG_FLONUM_VALUE(y));
2191 } else if (SG_COMPLEXP(y)) {
2192 return Sg_MakeComplex(Sg_Mul(x, SG_COMPLEX(y)->real),
2193 Sg_Mul(x, SG_COMPLEX(y)->imag));
2194 }
2195 }
2196 else if (SG_RATIONALP(x)) {
2197 if (SG_INTP(y)) {
2198 if (x == SG_MAKE_INT(0)) return x;
2199 else if (x == SG_MAKE_INT(1)) return y;
2200 else return Sg_RationalMul(x, y);
2201 } else if (SG_BIGNUMP(y) || SG_RATIONALP(y)) {
2202 return Sg_RationalMul(x, y);
2203 } else if (SG_FLONUMP(y)) {
2204 return Sg_MakeFlonum(Sg_GetDouble(x) * SG_FLONUM_VALUE(y));
2205 } else if (SG_COMPLEXP(y)) {
2206 return Sg_MakeComplex(Sg_Mul(x, SG_COMPLEX(y)->real),
2207 Sg_Mul(x, SG_COMPLEX(y)->imag));
2208 }
2209 }
2210 else if (SG_COMPLEXP(x)) {
2211 if (SG_INTP(y)) {
2212 if (x == SG_MAKE_INT(0)) return x;
2213 if (x == SG_MAKE_INT(1)) return y;
2214 return Sg_MakeComplex(Sg_Mul(SG_COMPLEX(x)->real, y),
2215 Sg_Mul(SG_COMPLEX(x)->imag, y));
2216 }
2217 else if (SG_BIGNUMP(y) || SG_RATIONALP(y)) {
2218 return Sg_MakeComplex(Sg_Mul(SG_COMPLEX(x)->real, y),
2219 Sg_Mul(SG_COMPLEX(x)->imag, y));
2220 }
2221 else if (SG_FLONUMP(y)) {
2222 return Sg_MakeComplex(Sg_Mul(SG_COMPLEX(x)->real, y),
2223 Sg_Mul(SG_COMPLEX(x)->imag, y));
2224 }
2225 else if (SG_COMPLEXP(y)) {
2226 SgObject real = Sg_Sub(Sg_Mul(SG_COMPLEX(x)->real, SG_COMPLEX(y)->real),
2227 Sg_Mul(SG_COMPLEX(x)->imag, SG_COMPLEX(y)->imag));
2228 SgObject imag = Sg_Add(Sg_Mul(SG_COMPLEX(x)->imag, SG_COMPLEX(y)->real),
2229 Sg_Mul(SG_COMPLEX(x)->real, SG_COMPLEX(y)->imag));
2230 return oprtr_norm_complex(real, imag);
2231 }
2232 }
2233 wte(SG_INTERN("*"), "number", SG_LIST2(x, y));
2234 return SG_UNDEF; /* dummy */
2235 }
2236
Sg_Div(SgObject x,SgObject y)2237 SgObject Sg_Div(SgObject x, SgObject y)
2238 {
2239 SgObject real, imag;
2240 if (SG_INTP(x)) {
2241 if (SG_INTP(y)) {
2242 if (y == SG_MAKE_INT(0)) goto div_by_zero;
2243 else if (x == SG_MAKE_INT(0)) return x;
2244 else if (y == SG_MAKE_INT(1)) return x;
2245 else return Sg_MakeRational(x, y);
2246 }
2247 else if (SG_BIGNUMP(y)) {
2248 if (x == SG_MAKE_INT(0)) return x;
2249 else return Sg_MakeRational(x, y);
2250 }
2251 else if (SG_RATIONALP(y)) {
2252 return Sg_MakeRational(Sg_Mul(x, SG_RATIONAL(y)->denominator),
2253 SG_RATIONAL(y)->numerator);
2254 }
2255 else if (SG_FLONUMP(y)) {
2256 if (SG_FLONUM_VALUE(y) == 0.0) goto a_normal;
2257 return Sg_MakeFlonum((double)SG_INT_VALUE(x) / SG_FLONUM_VALUE(y));
2258 }
2259 else if (SG_COMPLEXP(y)) {
2260 SgObject r2;
2261 real = SG_COMPLEX(y)->real;
2262 imag = SG_COMPLEX(y)->imag;
2263 r2 = Sg_Add(Sg_Mul(real, real), Sg_Mul(imag, imag));
2264 return Sg_MakeComplex(Sg_Div(Sg_Mul(real, x), r2),
2265 Sg_Negate(Sg_Div(Sg_Mul(x, imag), r2)));
2266 }
2267 }
2268 else if (SG_FLONUMP(x)) {
2269 if (SG_INTP(y)) {
2270 if (y == SG_MAKE_INT(0)) goto a_normal;
2271 else if (y == SG_MAKE_INT(1)) return x;
2272 else return Sg_MakeFlonum(SG_FLONUM_VALUE(x) / SG_INT_VALUE(y));
2273 }
2274 else if (SG_BIGNUMP(y) || SG_RATIONALP(y)) {
2275 return Sg_MakeFlonum(SG_FLONUM_VALUE(x) / Sg_GetDouble(y));
2276 }
2277 else if (SG_FLONUMP(y)) {
2278 if (SG_FLONUM_VALUE(y) == 0.0) goto a_normal;
2279 else return Sg_MakeFlonum(SG_FLONUM_VALUE(x) / SG_FLONUM_VALUE(y));
2280 }
2281 else if (SG_COMPLEXP(y)) {
2282 SgObject r2;
2283 real = SG_COMPLEX(y)->real;
2284 imag = SG_COMPLEX(y)->imag;
2285 r2 = Sg_Add(Sg_Mul(real, real), Sg_Mul(imag, imag));
2286 return Sg_MakeComplex(Sg_Div(Sg_Mul(real, x), r2),
2287 Sg_Negate(Sg_Div(Sg_Mul(x, imag), r2)));
2288 }
2289 }
2290 else if (SG_BIGNUMP(x)) {
2291 if (SG_INTP(y)) {
2292 if (y == SG_MAKE_INT(0)) goto div_by_zero;
2293 else if (y == SG_MAKE_INT(1)) return x;
2294 else return Sg_MakeRational(x, y);
2295 }
2296 else if (SG_BIGNUMP(y)) {
2297 return Sg_MakeRational(x, y);
2298 }
2299 else if (SG_RATIONALP(y)) {
2300 return Sg_MakeRational(Sg_Mul(SG_RATIONAL(y)->denominator, x),
2301 SG_RATIONAL(y)->numerator);
2302 }
2303 else if (SG_FLONUMP(y)) {
2304 if (SG_FLONUM_VALUE(y) == 0.0) goto a_normal;
2305 return Sg_MakeFlonum(Sg_GetDouble(x) / SG_FLONUM_VALUE(y));
2306 }
2307 else if (SG_COMPLEXP(y)) {
2308 SgObject r2;
2309 real = SG_COMPLEX(y)->real;
2310 imag = SG_COMPLEX(y)->imag;
2311 r2 = Sg_Add(Sg_Mul(real, real), Sg_Mul(imag, imag));
2312 return Sg_MakeComplex(Sg_Div(Sg_Mul(real, x), r2),
2313 Sg_Negate(Sg_Div(Sg_Mul(x, imag), r2)));
2314 }
2315 }
2316 else if (SG_RATIONALP(x)) {
2317 if (SG_INTP(y)) {
2318 if (y == SG_MAKE_INT(0)) goto div_by_zero;
2319 else if (y == SG_MAKE_INT(1)) return x;
2320 else return Sg_MakeRational(SG_RATIONAL(x)->numerator,
2321 Sg_Mul(SG_RATIONAL(x)->denominator, y));
2322 }
2323 else if (SG_BIGNUMP(y)) {
2324 return Sg_MakeRational(SG_RATIONAL(x)->numerator,
2325 Sg_Mul(SG_RATIONAL(x)->denominator, y));
2326 }
2327 else if (SG_RATIONALP(y)) {
2328 return Sg_RationalDiv(x, y);
2329 }
2330 else if (SG_FLONUMP(y)) {
2331 return Sg_MakeFlonum(Sg_GetDouble(x) / SG_FLONUM_VALUE(y));
2332 }
2333 else if (SG_COMPLEXP(y)) {
2334 SgObject r2;
2335 real = SG_COMPLEX(y)->real;
2336 imag = SG_COMPLEX(y)->imag;
2337 r2 = Sg_Add(Sg_Mul(real, real), Sg_Mul(imag, imag));
2338 return Sg_MakeComplex(Sg_Div(Sg_Mul(real, x), r2),
2339 Sg_Negate(Sg_Div(Sg_Mul(x, imag), r2)));
2340 }
2341 }
2342 else if (SG_COMPLEXP(x)) {
2343 real = SG_COMPLEX(x)->real;
2344 imag = SG_COMPLEX(x)->imag;
2345 if (SG_INTP(y)) {
2346 return Sg_MakeComplex(Sg_Div(real, y), Sg_Div(imag, y));
2347 }
2348 else if (SG_BIGNUMP(y) || SG_RATIONALP(y)) {
2349 return Sg_MakeComplex(Sg_Div(real, y), Sg_Div(imag, y));
2350 }
2351 else if (SG_FLONUMP(y)) {
2352 return Sg_MakeComplex(Sg_Div(real, y), Sg_Div(imag, y));
2353 }
2354 else if (SG_COMPLEXP(y)) {
2355 SgObject real2 = SG_COMPLEX(y)->real;
2356 SgObject imag2 = SG_COMPLEX(y)->imag;
2357 SgObject r2 = Sg_Add(Sg_Mul(real2, real2), Sg_Mul(imag2, imag2));
2358 SgObject real3 = Sg_Div(Sg_Add(Sg_Mul(real, real2),
2359 Sg_Mul(imag, imag2)), r2);
2360 SgObject imag3 = Sg_Div(Sg_Sub(Sg_Mul(imag, real2),
2361 Sg_Mul(real, imag2)), r2);
2362 return oprtr_norm_complex(real3, imag3);
2363 }
2364 }
2365 wte(SG_INTERN("/"), "number", SG_LIST2(x, y));
2366 return SG_UNDEF; /* dummy */
2367
2368 a_normal:
2369 {
2370 int s = Sg_Sign(x);
2371 if (s == 0) return SG_NAN;
2372 else {
2373 if (SG_INTP(y)) {
2374 only_x:
2375 if (s < 0) return SG_NEGATIVE_INFINITY;
2376 else return SG_POSITIVE_INFINITY;
2377 } else {
2378 union { double f64; int64_t i64; } d;
2379 d.f64 = SG_FLONUM_VALUE(y);
2380 /* if pisitive 0 is given then i64 is 0. */
2381 if (d.i64 < 0)
2382 if (s < 0) return SG_POSITIVE_INFINITY;
2383 else return SG_NEGATIVE_INFINITY;
2384 else goto only_x;
2385 }
2386 }
2387 }
2388 div_by_zero:
2389 Sg_AssertionViolation(SG_INTERN("/"), SG_MAKE_STRING("undefined for 0"),
2390 SG_LIST2(x, y));
2391 return SG_UNDEF;
2392 }
2393
Sg_Quotient(SgObject x,SgObject y,SgObject * rem)2394 SgObject Sg_Quotient(SgObject x, SgObject y, SgObject *rem)
2395 {
2396 SgObject real;
2397 SgObject imag;
2398 double rx, ry;
2399
2400 #define do_complex(t, label) \
2401 real = SG_COMPLEX(t)->real; \
2402 imag = SG_COMPLEX(t)->imag; \
2403 if (Sg_ZeroP(imag)) { \
2404 (t) = real; \
2405 goto label; \
2406 } \
2407 else goto bad_arg; \
2408
2409 start_again:
2410 if (SG_INTP(x)) {
2411 if (SG_INT_VALUE(x) == 0) {
2412 if (rem) *rem = SG_MAKE_INT(0);
2413 return SG_MAKE_INT(0);
2414 }
2415 fixnum_again:
2416 if (SG_INTP(y)) {
2417 long q, r;
2418 if (SG_INT_VALUE(y) == 0) {
2419 goto div_by_zero;
2420 }
2421 q = SG_INT_VALUE(x) / SG_INT_VALUE(y);
2422 if (rem) {
2423 r = SG_INT_VALUE(x) % SG_INT_VALUE(y);
2424 *rem = SG_MAKE_INT(r);
2425 }
2426 return SG_MAKE_INT(q);
2427 }
2428 else if (SG_BIGNUMP(y)) {
2429 SgObject qr = Sg_BignumDivRem(SG_BIGNUM(Sg_MakeBignumFromSI(SG_INT_VALUE(x))),
2430 SG_BIGNUM(y));
2431 if (rem) *rem = SG_CDR(qr);
2432 return SG_CAR(qr);
2433 }
2434 else if (SG_FLONUMP(y)) {
2435 rx = (double)SG_INT_VALUE(x);
2436 ry = SG_FLONUM_VALUE(y);
2437 if (ry != floor(ry)) goto bad_argy;
2438 goto do_flonum;
2439 }
2440 else if (SG_COMPLEXP(y)) {
2441 do_complex(y, fixnum_again);
2442 }
2443 else goto bad_argy;
2444 }
2445 if (SG_BIGNUMP(x)) {
2446 bignum_again:
2447 if (SG_INTP(y)) {
2448 long r;
2449 SgObject q = Sg_BignumDivSI(SG_BIGNUM(x), SG_INT_VALUE(y), &r);
2450 if (rem) *rem = SG_MAKE_INT(r);
2451 return q;
2452 }
2453 else if (SG_BIGNUMP(y)) {
2454 SgObject qr = Sg_BignumDivRem(SG_BIGNUM(x), SG_BIGNUM(y));
2455 if (rem) *rem = SG_CDR(qr);
2456 return SG_CAR(qr);
2457 }
2458 else if (SG_FLONUMP(y)) {
2459 rx = Sg_BignumToDouble(SG_BIGNUM(x));
2460 ry = SG_FLONUM_VALUE(y);
2461 if (ry != floor(ry)) goto bad_argy;
2462 else goto do_flonum;
2463 }
2464 else if (SG_COMPLEXP(y)) {
2465 do_complex(y, bignum_again);
2466 }
2467 goto bad_argy;
2468 }
2469 else if (SG_COMPLEXP(x)) {
2470 do_complex(x, start_again);
2471 }
2472 else if (SG_FLONUMP(x)) {
2473 rx = SG_FLONUM_VALUE(x);
2474 if (rx != floor(rx)) goto bad_arg;
2475 flonum_again:
2476 if (SG_INTP(y)) {
2477 ry = (double)SG_INT_VALUE(y);
2478 } else if (SG_BIGNUMP(y)) {
2479 ry = Sg_BignumToDouble(SG_BIGNUM(y));
2480 } else if (SG_FLONUMP(y)) {
2481 ry = SG_FLONUM_VALUE(y);
2482 if (ry != floor(ry)) goto bad_argy;
2483 } else if (SG_COMPLEXP(y)) {
2484 do_complex(y, flonum_again);
2485 } else {
2486 goto bad_argy;
2487 }
2488 do_flonum:
2489 {
2490 double q;
2491 if (ry == 0.0) goto div_by_zero;
2492 q = (rx*ry > 0)? floor(rx/ry) : ceil(rx/ry);
2493 if (rem) {
2494 double rr = roundeven(rx - q*ry);
2495 *rem = Sg_MakeFlonum(rr);
2496 }
2497 return Sg_MakeFlonum(q);
2498 }
2499 }
2500 goto bad_arg;
2501
2502 div_by_zero:
2503 Sg_AssertionViolation(SG_INTERN("quotient"),
2504 SG_MAKE_STRING("attempt to calculate a quotient by zero"),
2505 SG_LIST2(x, y));
2506 bad_argy:
2507 x = y;
2508 bad_arg:
2509 wte(SG_INTERN("quotient"), "integer", x);
2510 return SG_UNDEF; /* dummy */
2511 #undef do_complex
2512 }
2513
Sg_Modulo(SgObject x,SgObject y,int remp)2514 SgObject Sg_Modulo(SgObject x, SgObject y, int remp)
2515 {
2516 SgObject real, imag;
2517 SgObject bx;
2518 double rx, ry;
2519
2520 #define do_complex(t, label) \
2521 real = SG_COMPLEX(t)->real; \
2522 imag = SG_COMPLEX(t)->imag; \
2523 if (Sg_ZeroP(imag)) { \
2524 (t) = real; \
2525 goto label; \
2526 } \
2527 else goto bad_arg; \
2528
2529 start_again:
2530 if (SG_INTP(x)) {
2531 fixnum_again:
2532 if (SG_INTP(y)) {
2533 long r;
2534 if (SG_EQ((y), SG_MAKE_INT(0))) goto div_by_zero;
2535 r = SG_INT_VALUE(x) % SG_INT_VALUE(y);
2536 if (!remp && r) {
2537 if ((SG_INT_VALUE(x) > 0 && SG_INT_VALUE(y) < 0)
2538 || (SG_INT_VALUE(x) < 0 && SG_INT_VALUE(y) > 0)) {
2539 r += SG_INT_VALUE(y);
2540 }
2541 }
2542 return SG_MAKE_INT(r);
2543 }
2544 else if (SG_BIGNUMP(y)) {
2545 bx = Sg_MakeBignumFromSI(SG_INT_VALUE(x));
2546 goto do_bignumy;
2547 }
2548 else if (SG_FLONUMP(y)) {
2549 rx = (double)SG_INT_VALUE(x);
2550 ry = SG_FLONUM_VALUE(y);
2551 if (ry != floor(ry)) goto bad_argy;
2552 else goto do_flonum;
2553 }
2554 else if (SG_COMPLEXP(y)) {
2555 do_complex(y, fixnum_again);
2556 }
2557 else goto bad_arg;
2558 }
2559 if (SG_BIGNUMP(x)) {
2560 bignum_again:
2561 if (SG_INTP(y)) {
2562 #if 0
2563 long iy = SG_INT_VALUE(y);
2564 long rem;
2565 Sg_BignumDivSI(SG_BIGNUM(x), iy, &rem);
2566 if (!remp
2567 && rem
2568 && ((SG_BIGNUM_GET_SIGN(x) < 0 && iy > 0)
2569 || (SG_BIGNUM_GET_SIGN(x) > 0 && iy < 0))) {
2570 return SG_MAKE_INT(iy + rem);
2571 }
2572 else return SG_MAKE_INT(rem);
2573 #else
2574 if (SG_INT_VALUE(y) == 0) goto div_by_zero;
2575 return Sg_BignumModuloSI(SG_BIGNUM(x), SG_INT_VALUE(y), remp);
2576 #endif
2577 }
2578 else if (SG_BIGNUMP(y)) {
2579 bx = x;
2580 do_bignumy:
2581 return Sg_BignumModulo(bx, SG_BIGNUM(y), remp);
2582 }
2583 else if (SG_FLONUMP(y)) {
2584 rx = Sg_BignumToDouble(SG_BIGNUM(x));
2585 ry = SG_FLONUM_VALUE(y);
2586 if (ry != floor(ry)) goto bad_argy;
2587 else goto do_flonum;
2588 }
2589 else if (SG_COMPLEXP(y)) {
2590 do_complex(y, bignum_again);
2591 }
2592 else goto bad_argy;
2593 }
2594 else if (SG_FLONUMP(x)) {
2595 double rem;
2596 rx = SG_FLONUM_VALUE(x);
2597 flonum_again:
2598 if (rx != floor(rx)) goto bad_arg;
2599 if (SG_INTP(y)) {
2600 ry = (double)SG_INT_VALUE(y);
2601 } else if (SG_BIGNUMP(y)) {
2602 ry = Sg_BignumToDouble(y);
2603 } else if (SG_FLONUMP(y)) {
2604 ry = SG_FLONUM_VALUE(y);
2605 } else if (SG_COMPLEXP(y)) {
2606 do_complex(y, flonum_again);
2607 } else {
2608 goto bad_argy;
2609 }
2610 do_flonum:
2611 if (ry == 0.0) goto div_by_zero;
2612 rem = fmod(rx, ry);
2613 if (!remp) {
2614 if (rem == 0.0) return Sg_MakeFlonum(0.0);
2615 if ((rx > 0 && ry < 0) || (rx < 0 && ry > 0)) {
2616 rem += ry;
2617 }
2618 }
2619 return Sg_MakeFlonum(rem);
2620 }
2621 else if (SG_COMPLEXP(x)) {
2622 do_complex(x, start_again);
2623 }
2624 else goto bad_arg;
2625
2626 div_by_zero:
2627 Sg_AssertionViolation((remp) ? SG_INTERN("remainder") : SG_INTERN("modulo"),
2628 SG_MAKE_STRING("attempt to calculate a remainder/modulo by zero"),
2629 SG_LIST2(x, y));
2630 bad_argy:
2631 x = y;
2632 bad_arg:
2633 wte((remp) ? SG_INTERN("remainder") : SG_INTERN("modulo"),
2634 "integer", x);
2635 return SG_UNDEF; /* dummy */
2636
2637 #undef do_complex
2638 }
2639
cos_pi(double x)2640 static double cos_pi(double x)
2641 {
2642 return cos(M_PI * x);
2643 }
2644
sin_pi(double x)2645 static double sin_pi(double x)
2646 {
2647 return sin(M_PI * x);
2648 }
2649
2650
expt_body(SgObject x,SgObject y)2651 static SgObject expt_body(SgObject x, SgObject y)
2652 {
2653 if (SG_FLONUMP(x) && SG_FLONUM_VALUE(x) == 0.0) {
2654 if (SG_COMPLEXP(y)) {
2655 if (Sg_PositiveP(SG_COMPLEX(y)->real)) return SG_FL_POSITIVE_ZERO;
2656 } else {
2657 if (Sg_PositiveP(y)) {
2658 int exp, sign;
2659 Sg_DecodeFlonum(SG_FLONUM_VALUE(x), &exp, &sign);
2660 if (sign < 0)
2661 return Sg_OddP(y) ? SG_FL_NEGATIVE_ZERO : SG_FL_POSITIVE_ZERO;
2662 return SG_FL_POSITIVE_ZERO;
2663 }
2664 }
2665 }
2666 if (Sg_ExactP(y)) {
2667 if (SG_INTP(y)) {
2668 if (SG_INT_VALUE(y) == 0) {
2669 if (Sg_ExactP(x)) return SG_MAKE_INT(1);
2670 else return Sg_MakeFlonum(1.0);
2671 } else if (SG_FLONUMP(x))
2672 return Sg_MakeFlonum(pow(SG_FLONUM_VALUE(x),
2673 (double)SG_INT_VALUE(y)));
2674 else return oprtr_expt(x, SG_INT_VALUE(y));
2675 }
2676 else if (SG_BIGNUMP(y)) {
2677 if (!Sg_ExactP(x)) {
2678 /* Sg_Error(UC("exact number required, but got %S"), x); */
2679 return SG_UNDEF;
2680 }
2681 else if (SG_REALP(y)) {
2682 double n = Sg_BignumToDouble(SG_BIGNUM(y));
2683 return Sg_MakeFlonum(pow(Sg_GetDouble(x), n));
2684 }
2685 else return Sg_Exp(Sg_Mul(y, Sg_Log(x)));
2686 }
2687 else if (SG_RATIONALP(y)) {
2688 double n = Sg_GetDouble(y);
2689 if (SG_REALP(x) && !Sg_NegativeP(x))
2690 return Sg_MakeFlonum(pow(Sg_GetDouble(x), n));
2691 else
2692 return Sg_Exp(Sg_Mul(y, Sg_Log(x)));
2693 }
2694 else if (SG_COMPLEXP(y)) {
2695 return Sg_Exp(Sg_Mul(y, Sg_Log(x)));
2696 }
2697 else goto bad_arg;
2698 } else {
2699 if (SG_FLONUMP(y)) {
2700 if (SG_REALP(x)) {
2701 double n = SG_FLONUM_VALUE(y);
2702 double m = Sg_GetDouble(x);
2703 if (n < 0 && !Sg_IntegerP(x)) {
2704 /* x^y = exp(y * log(x)) = exp(y*log(|x|))*exp(y*arg(x)*i) */
2705 double mag = exp(n * log(-m));
2706 return Sg_MakeComplex(Sg_MakeFlonum(mag * cos_pi(n)),
2707 Sg_MakeFlonum(mag * sin_pi(n)));
2708 } else return Sg_MakeFlonum(pow(m, n));
2709 }
2710 else return Sg_Exp(Sg_Mul(y, Sg_Log(x)));
2711 } else {
2712 return Sg_Exp(Sg_Mul(y, Sg_Log(x)));
2713 }
2714 }
2715 bad_arg:
2716 /* Sg_Error(UC("real number required, but got %S"), y); */
2717 return SG_UNDEF; /* dummy */
2718 }
2719
Sg_Expt(SgObject x,SgObject y)2720 SgObject Sg_Expt(SgObject x, SgObject y)
2721 {
2722 if (x == SG_MAKE_INT(1)) {
2723 if (Sg_ExactP(y)) return SG_MAKE_INT(1);
2724 else return Sg_MakeFlonum(1.0);
2725 }
2726 else if (x == SG_MAKE_INT(-1) && Sg_ExactP(y)) {
2727 if (Sg_OddP(y)) return SG_MAKE_INT(-1);
2728 else return SG_MAKE_INT(1);
2729 }
2730 else if (x == SG_MAKE_INT(0)) {
2731 if (Sg_ZeroP(y)) {
2732 if (Sg_ExactP(y)) return SG_MAKE_INT(1);
2733 else return Sg_MakeFlonum(1.0);
2734 }
2735 else if (Sg_RealValuedP(y)) {
2736 if (Sg_NegativeP(y)) return Sg_MakeComplex(SG_NAN, SG_NAN);
2737 else if (Sg_ExactP(y)) return SG_MAKE_INT(0);
2738 else return Sg_MakeFlonum(0.0);
2739 } else {
2740 ASSERT(SG_COMPLEXP(y));
2741 if (Sg_PositiveP(SG_COMPLEX(y)->real)) {
2742 if (Sg_ExactP(y)) return SG_MAKE_INT(0);
2743 else return Sg_MakeFlonum(0.0);
2744 }
2745 else return Sg_MakeComplex(SG_NAN, SG_NAN);
2746 }
2747 }
2748 else if (Sg_ExactP(x) && SG_BIGNUMP(y)) {
2749 SgObject msg =
2750 Sg_Sprintf(UC("expt: calculated number is too big to fit into memory %S %S"),
2751 x, y);
2752 Sg_ImplementationRestrictionViolation(SG_INTERN("expt"), msg,
2753 SG_LIST2(x, y));
2754 return SG_UNDEF;
2755 }
2756 else return expt_body(x, y);
2757 }
2758
Sg_Exp(SgObject obj)2759 SgObject Sg_Exp(SgObject obj)
2760 {
2761 if (SG_INTP(obj)) {
2762 if (SG_INT_VALUE(obj) == 0) return SG_MAKE_INT(1);
2763 else return Sg_MakeFlonum(exp((double)SG_INT_VALUE(obj)));
2764 }
2765 else if (SG_COMPLEXP(obj)) {
2766 double real = Sg_GetDouble(SG_COMPLEX(obj)->real);
2767 double imag = Sg_GetDouble(SG_COMPLEX(obj)->imag);
2768 double a = exp(real);
2769 return Sg_MakeComplex(Sg_MakeFlonum(a * cos(imag)),
2770 Sg_MakeFlonum(a * sin(imag)));
2771 }
2772 else if (SG_REALP(obj)) return Sg_MakeFlonum(exp(Sg_GetDouble(obj)));
2773 wte(SG_INTERN("exp"), "real number", obj);
2774 return SG_UNDEF; /* dummy */
2775 }
2776
Sg_Sin(SgObject obj)2777 SgObject Sg_Sin(SgObject obj)
2778 {
2779 if (SG_INTP(obj)) {
2780 if (obj == SG_MAKE_INT(0)) return obj;
2781 else return Sg_MakeFlonum(sin((double)SG_INT_VALUE(obj)));
2782 }
2783 else if (SG_COMPLEXP(obj)) {
2784 double real = Sg_GetDouble(SG_COMPLEX(obj)->real);
2785 double imag = Sg_GetDouble(SG_COMPLEX(obj)->imag);
2786 double e = exp(imag);
2787 double f = 1.0 / e;
2788 return Sg_MakeComplex(Sg_MakeFlonum(0.5 * sin(real) * (e + f)),
2789 Sg_MakeFlonum(0.5 * cos(real) * (e - f)));
2790 }
2791 else if (SG_REALP(obj)) return Sg_MakeFlonum(sin(Sg_GetDouble(obj)));
2792 wte(SG_INTERN("sin"), "number", obj);
2793 return SG_UNDEF; /* dummy */
2794 }
2795
Sg_Cos(SgObject obj)2796 SgObject Sg_Cos(SgObject obj)
2797 {
2798 if (SG_INTP(obj)) {
2799 if (obj == SG_MAKE_INT(0)) return SG_MAKE_INT(1);
2800 else return Sg_MakeFlonum(cos((double)SG_INT_VALUE(obj)));
2801 }
2802 else if (SG_COMPLEXP(obj)) {
2803 double real = Sg_GetDouble(SG_COMPLEX(obj)->real);
2804 double imag = Sg_GetDouble(SG_COMPLEX(obj)->imag);
2805 double e = exp(imag);
2806 double f = 1.0 / e;
2807 return Sg_MakeComplex(Sg_MakeFlonum(0.5 * cos(real) * (f + e)),
2808 Sg_MakeFlonum(0.5 * sin(real) * (f - e)));
2809 }
2810 else if (SG_REALP(obj)) return Sg_MakeFlonum(cos(Sg_GetDouble(obj)));
2811 wte(SG_INTERN("cos"), "number", obj);
2812 return SG_UNDEF; /* dummy */
2813 }
2814
Sg_Tan(SgObject obj)2815 SgObject Sg_Tan(SgObject obj)
2816 {
2817 if (SG_INTP(obj)) {
2818 if (obj == SG_MAKE_INT(0)) return obj;
2819 else return Sg_MakeFlonum(tan((double)SG_INT_VALUE(obj)));
2820 }
2821 else if (SG_COMPLEXP(obj)) {
2822 double real = Sg_GetDouble(SG_COMPLEX(obj)->real);
2823 double imag = Sg_GetDouble(SG_COMPLEX(obj)->imag);
2824 double e = exp(2.0 * imag);
2825 double f = 1.0 / e;
2826 double d = cos(2.0 * real) + 0.5 * (e + f);
2827 return Sg_MakeComplex(Sg_MakeFlonum(sin(2.0 * real) / d),
2828 Sg_MakeFlonum(0.5 * (e - f) / d));
2829 }
2830 else if (SG_REALP(obj)) return Sg_MakeFlonum(tan(Sg_GetDouble(obj)));
2831 wte(SG_INTERN("tan"), "number", obj);
2832 return SG_UNDEF; /* dummy */
2833 }
2834
Sg_Asin(SgObject obj)2835 SgObject Sg_Asin(SgObject obj)
2836 {
2837 SgComplex *cn;
2838 SgObject ans;
2839 if (SG_EQ(obj, SG_MAKE_INT(0))) return obj;
2840
2841 if (SG_REALP(obj) || (SG_COMPLEXP(obj) && Sg_ZeroP(SG_COMPLEX(obj)->imag))) {
2842 double x = Sg_GetDouble(obj);
2843 if (x >= -1.0 && x <= 1.0) return Sg_MakeFlonum(asin(Sg_GetDouble(obj)));
2844 else if (x < 0.0) return Sg_Negate(Sg_Asin(Sg_MakeFlonum(-x)));
2845 cn = Sg_MakeComplex(Sg_MakeFlonum(0.0), Sg_MakeFlonum(x));
2846 } else {
2847 ASSERT(SG_COMPLEXP(obj));
2848 if (Sg_PositiveP(SG_COMPLEX(obj)->imag))
2849 return Sg_Negate(Sg_Asin(Sg_Negate(obj)));
2850 cn = Sg_MakeComplex(Sg_Negate(SG_COMPLEX(obj)->imag),
2851 SG_COMPLEX(obj)->real);
2852 }
2853 ans = Sg_Log(Sg_Add(Sg_Sqrt(Sg_Sub(SG_MAKE_INT(1), Sg_Mul(obj, obj))),
2854 cn));
2855 if (SG_COMPLEXP(ans)) {
2856 return Sg_MakeComplex(Sg_MakeFlonum(Sg_GetDouble(SG_COMPLEX(ans)->imag)),
2857 Sg_MakeFlonum(-Sg_GetDouble(SG_COMPLEX(ans)->real)));
2858 }
2859 else
2860 return Sg_MakeComplex(Sg_MakeFlonum(0.0),
2861 Sg_MakeFlonum(-Sg_GetDouble(ans)));
2862
2863 }
2864
Sg_Acos(SgObject obj)2865 SgObject Sg_Acos(SgObject obj)
2866 {
2867 if (SG_EQ(obj, SG_MAKE_INT(1))) return SG_MAKE_INT(0);
2868
2869 if (SG_REALP(obj) || (SG_COMPLEXP(obj) && Sg_ZeroP(SG_COMPLEX(obj)->imag))) {
2870 double x = Sg_GetDouble(obj);
2871 if (x >= -1.0 && x <= 1.0) return Sg_MakeFlonum(acos(Sg_GetDouble(obj)));
2872 }
2873 return Sg_Sub(Sg_MakeFlonum(M_PI / 2.0), Sg_Asin(obj));
2874 }
2875
Sg_Atan(SgObject obj)2876 SgObject Sg_Atan(SgObject obj)
2877 {
2878 if (SG_INTP(obj)) {
2879 if (obj == SG_MAKE_INT(0)) return obj;
2880 return Sg_MakeFlonum(atan((double)SG_INT_VALUE(obj)));
2881 }
2882 if (SG_COMPLEXP(obj)) {
2883 SgComplex *cn = Sg_MakeComplex(Sg_Negate(SG_COMPLEX(obj)->imag),
2884 SG_COMPLEX(obj)->real);
2885 SgObject ans = Sg_Log(Sg_Div(Sg_Add(SG_MAKE_INT(1), cn),
2886 Sg_Sub(SG_MAKE_INT(1), cn)));
2887 if (SG_COMPLEXP(ans)) {
2888 return Sg_MakeComplex(Sg_MakeFlonum(0.5 * Sg_GetDouble(SG_COMPLEX(ans)->imag)),
2889 Sg_MakeFlonum(-0.5 * Sg_GetDouble(SG_COMPLEX(ans)->real)));
2890 }
2891 return Sg_MakeComplex(Sg_MakeFlonum(0.0),
2892 Sg_MakeFlonum(-0.5 * Sg_GetDouble(ans)));
2893 }
2894 if (SG_REALP(obj)) return Sg_MakeFlonum(atan(Sg_GetDouble(obj)));
2895 wte(SG_INTERN("atan"), "number", obj);
2896 return SG_UNDEF;
2897 }
2898
Sg_Atan2(SgObject x,SgObject y)2899 SgObject Sg_Atan2(SgObject x, SgObject y)
2900 {
2901 if (SG_EQ(x, SG_MAKE_INT(0))) return x;
2902 return Sg_MakeFlonum(atan2(Sg_GetDouble(x), Sg_GetDouble(y)));
2903 }
2904
either_nan_p(SgObject arg0,SgObject arg1)2905 static inline int either_nan_p(SgObject arg0, SgObject arg1)
2906 {
2907 if (SG_FLONUMP(arg0) && isnan(SG_FLONUM_VALUE(arg0))) return TRUE;
2908 if (SG_FLONUMP(arg1) && isnan(SG_FLONUM_VALUE(arg1))) return TRUE;
2909 return FALSE;
2910 }
2911
Sg_NumEq(SgObject x,SgObject y)2912 int Sg_NumEq(SgObject x, SgObject y)
2913 {
2914 /* not to break zero?'s compiler builtin inliner, sucks!! */
2915 if (SG_INTP(y) && SG_EQ(y, SG_MAKE_INT(0)) && Sg_ZeroP(x)) return TRUE;
2916 if (SG_COMPLEXP(x)) {
2917 if (SG_COMPLEXP(y)) {
2918 return ((Sg_NumCmp(SG_COMPLEX(x)->real, SG_COMPLEX(y)->real) == 0)
2919 && (Sg_NumCmp(SG_COMPLEX(x)->imag, SG_COMPLEX(y)->imag) == 0));
2920 }
2921 if (Sg_ZeroP(SG_COMPLEX(x)->imag)) {
2922 return Sg_NumEq(SG_COMPLEX(x)->real, y);
2923 }
2924 return FALSE;
2925 } else {
2926 if (SG_COMPLEXP(y)) {
2927 if (Sg_ZeroP(SG_COMPLEX(y)->imag)) {
2928 return Sg_NumEq(x, SG_COMPLEX(y)->real);
2929 }
2930 return FALSE;
2931 }
2932 if (either_nan_p(x, y)) return FALSE;
2933 return (Sg_NumCmp(x, y) == 0);
2934 }
2935 }
2936
Sg_NumCmp(SgObject x,SgObject y)2937 int Sg_NumCmp(SgObject x, SgObject y)
2938 {
2939 SgObject badnum = SG_FALSE;
2940 /* on WATCOM, somehow NAN is bigger than 0 */
2941 #ifdef __WATCOMC__
2942 #define nan_return(_r) do { if (isnan(_r)) return 0; } while(0)
2943 #else
2944 #define nan_return(_r) /* dummy */
2945 #endif
2946 start_again:
2947 if (SG_INTP(x)) {
2948 int_again:
2949 if (SG_INTP(y)) {
2950 long r = SG_INT_VALUE(x) - SG_INT_VALUE(y);
2951 if (r < 0) return -1;
2952 else if (r > 0) return 1;
2953 else return 0;
2954 }
2955 else if (SG_FLONUMP(y)) {
2956 double r = SG_INT_VALUE(x) - SG_FLONUM_VALUE(y);
2957 nan_return(r);
2958 if (r < 0) return -1;
2959 else if (r > 0) return 1;
2960 else return 0;
2961 }
2962 else if (SG_BIGNUMP(y)) {
2963 return Sg_BignumCmp(SG_BIGNUM(Sg_MakeBignumFromSI(SG_INT_VALUE(x))),
2964 SG_BIGNUM(y));
2965 }
2966 else if (SG_RATIONALP(y)) {
2967 if (SG_MAKE_INT(0) == x) {
2968 return -Sg_Sign(y);
2969 } else {
2970 /* roughly estimate the result by coercing the ratnum to double */
2971 double y2 = Sg_GetDouble(y);
2972 double r = SG_INT_VALUE(x) - y2;
2973 double err = y2 * 2.0e-52;
2974 if (r < -err) return -1;
2975 else if (r > err) return 1;
2976 else return Sg_NumCmp(Sg_Mul(x, SG_RATIONAL(y)->denominator),
2977 SG_RATIONAL(y)->numerator);
2978 }
2979 } else if (SG_COMPLEXP(y)) {
2980 if (Sg_ZeroP(SG_COMPLEX(y)->imag)) {
2981 y = SG_COMPLEX(y)->real;
2982 goto int_again;
2983 }
2984 }
2985 badnum = y;
2986 }
2987 else if (SG_FLONUMP(x)) {
2988 flo_again:
2989 if (SG_INTP(y)) {
2990 double r = SG_FLONUM_VALUE(x) - SG_INT_VALUE(y);
2991 nan_return(r);
2992 if (r < 0) return -1;
2993 else if (r > 0) return 1;
2994 else return 0;
2995 }
2996 else if (SG_FLONUMP(y)) {
2997 double r = SG_FLONUM_VALUE(x) - SG_FLONUM_VALUE(y);
2998 nan_return(r);
2999 if (r < 0) return -1;
3000 else if (r > 0) return 1;
3001 else return 0;
3002 }
3003 else if (SG_BIGNUMP(y) || SG_RATIONALP(y)) {
3004 double r = SG_FLONUM_VALUE(x) - Sg_GetDouble(y);
3005 nan_return(r);
3006 if (r < 0) return -1;
3007 else if (r > 0) return 1;
3008 else {
3009 return Sg_NumCmp(Sg_Exact(x), y);
3010 }
3011 } else if (SG_COMPLEXP(y)) {
3012 if (Sg_ZeroP(SG_COMPLEX(y)->imag)) {
3013 y = SG_COMPLEX(y)->real;
3014 goto flo_again;
3015 }
3016 }
3017 badnum = y;
3018 }
3019 else if (SG_BIGNUMP(x)) {
3020 big_again:
3021 if (SG_INTP(y)) {
3022 return Sg_BignumCmp(SG_BIGNUM(x),
3023 SG_BIGNUM(Sg_MakeBignumFromSI(SG_INT_VALUE(y))));
3024 }
3025 else if (SG_FLONUMP(y)) {
3026 return -Sg_NumCmp(y, x);
3027 }
3028 else if (SG_BIGNUMP(y)) {
3029 return Sg_BignumCmp(SG_BIGNUM(x), SG_BIGNUM(y));
3030 }
3031 else if (SG_RATIONALP(y)) {
3032 SgObject d1 = SG_RATIONAL(y)->denominator;
3033 return Sg_NumCmp(Sg_Mul(x, d1),
3034 SG_RATIONAL(y)->numerator);
3035 }
3036 else if (SG_COMPLEXP(y)) {
3037 if (Sg_ZeroP(SG_COMPLEX(y)->imag)) {
3038 y = SG_COMPLEX(y)->real;
3039 goto big_again;
3040 }
3041 }
3042 badnum = y;
3043 }
3044 else if (SG_RATIONALP(x)) {
3045 rat_again:
3046 if (SG_INTP(y) || SG_BIGNUMP(y) || SG_FLONUMP(y)) {
3047 return -Sg_NumCmp(y, x);
3048 }
3049 if (SG_RATIONALP(y)) {
3050 SgObject nx = SG_RATIONAL(x)->numerator, dx = SG_RATIONAL(x)->denominator;
3051 SgObject ny = SG_RATIONAL(y)->numerator, dy = SG_RATIONAL(y)->denominator;
3052 int sx = Sg_Sign(nx), sy = Sg_Sign(ny), d, n;
3053
3054 if (sx < sy) return -1;
3055 if (sx > sy) return 1;
3056 d = Sg_NumCmp(dx, dy);
3057 if (d == 0) return Sg_NumCmp(nx, ny);
3058 else if ((sx > 0 && sy > 0) || (sx < 0 && sy < 0)) {
3059 n = Sg_NumCmp(nx, ny) * sx;
3060 if (d > 0 && n <= 0) return -sx;
3061 if (d < 0 && n >= 0) return sx;
3062 }
3063 return Sg_NumCmp(Sg_Mul(nx, dy),
3064 Sg_Mul(ny, dx));
3065 } else if (SG_COMPLEXP(y)) {
3066 if (Sg_ZeroP(SG_COMPLEX(y)->imag)) {
3067 y = SG_COMPLEX(y)->real;
3068 goto rat_again;
3069 }
3070 }
3071 badnum = y;
3072 }
3073 else if (SG_COMPLEXP(x)) {
3074 if (Sg_ZeroP(SG_COMPLEX(x)->imag)) {
3075 x = SG_COMPLEX(x)->real;
3076 goto start_again;
3077 }
3078 badnum = x;
3079 }
3080 if (SG_FALSEP(badnum)) badnum = x;
3081
3082 wte4(SG_FALSE, "real number", badnum, SG_LIST2(x, y));
3083 return 0; /* dummy */
3084 }
3085
Sg_NumGt(SgObject x,SgObject y)3086 int Sg_NumGt(SgObject x, SgObject y)
3087 {
3088 if (either_nan_p(x, y)) return FALSE;
3089 return (Sg_NumCmp(x, y) > 0);
3090 }
Sg_NumGe(SgObject x,SgObject y)3091 int Sg_NumGe(SgObject x, SgObject y)
3092 {
3093 if (either_nan_p(x, y)) return FALSE;
3094 return (Sg_NumCmp(x, y) >= 0);
3095 }
3096
Sg_NumLt(SgObject x,SgObject y)3097 int Sg_NumLt(SgObject x, SgObject y)
3098 {
3099 if (either_nan_p(x, y)) return FALSE;
3100 return (Sg_NumCmp(x, y) < 0);
3101 }
Sg_NumLe(SgObject x,SgObject y)3102 int Sg_NumLe(SgObject x, SgObject y)
3103 {
3104 if (either_nan_p(x, y)) return FALSE;
3105 return (Sg_NumCmp(x, y) <= 0);
3106 }
3107
Sg_Abs(SgObject obj)3108 SgObject Sg_Abs(SgObject obj)
3109 {
3110 if (SG_INTP(obj)) {
3111 long v = SG_INT_VALUE(obj);
3112 if (v < 0) {
3113 obj = Sg_MakeInteger(-v);
3114 }
3115 } else if (SG_BIGNUMP(obj)) {
3116 if (SG_BIGNUM_GET_SIGN(obj) < 0) {
3117 obj = Sg_BignumCopy(SG_BIGNUM(obj));
3118 SG_BIGNUM_SET_SIGN(obj, 1);
3119 }
3120 } else if (SG_FLONUMP(obj)) {
3121 double v = SG_FLONUM_VALUE(obj);
3122 if (signbit(v)) return Sg_MakeFlonum(-v);
3123 } else if (SG_RATIONALP(obj)) {
3124 if (Sg_Sign(SG_RATIONAL(obj)->numerator) < 0) {
3125 obj = Sg_MakeRational(Sg_Negate(SG_RATIONAL(obj)->numerator),
3126 SG_RATIONAL(obj)->denominator);
3127 }
3128 } else if (SG_COMPLEXP(obj)) {
3129 SgObject r = SG_COMPLEX(obj)->real;
3130 SgObject i = SG_COMPLEX(obj)->imag;
3131 SgObject a = Sg_Sqrt(Sg_Add(Sg_Mul(r, r), Sg_Mul(i, i)));
3132 return a;
3133 } else {
3134 wte(SG_INTERN("abs"), "number", obj);
3135 }
3136 return obj;
3137 }
3138
Sg_Sqrt(SgObject obj)3139 SgObject Sg_Sqrt(SgObject obj)
3140 {
3141 if (SG_INTP(obj)) {
3142 long value = SG_INT_VALUE(obj);
3143 if (value == 0) return SG_MAKE_INT(0);
3144 if (value > 0) {
3145 double root = sqrt((double)value);
3146 long iroot = (long)floor(root);
3147 if (iroot * iroot == value) return SG_MAKE_INT(iroot);
3148 return Sg_MakeFlonum(root);
3149 } else {
3150 long iroot;
3151 value = -value;
3152 iroot = (long)floor(sqrt((double)value));
3153 if (iroot * iroot == value)
3154 return Sg_MakeComplex(SG_MAKE_INT(0), SG_MAKE_INT(iroot));
3155 return Sg_MakeComplex(Sg_MakeFlonum(0.0),
3156 Sg_MakeFlonum(sqrt((double)value)));
3157 }
3158 }
3159 if (SG_BIGNUMP(obj)) {
3160 return Sg_BignumSqrt(SG_BIGNUM(obj));
3161 }
3162 if (SG_RATIONALP(obj)) {
3163 SgObject nume, deno;
3164 int comp = FALSE;
3165 if (Sg_NegativeP(SG_RATIONAL(obj)->numerator)) {
3166 nume = Sg_Sqrt(Sg_Negate(SG_RATIONAL(obj)->numerator));
3167 comp = TRUE;
3168 } else {
3169 nume = Sg_Sqrt(SG_RATIONAL(obj)->numerator);
3170 }
3171 deno = Sg_Sqrt(SG_RATIONAL(obj)->denominator);
3172 if ((SG_INTP(nume) || SG_BIGNUMP(nume))
3173 && (SG_INTP(deno) || SG_BIGNUMP(deno))) {
3174 if (comp) return Sg_MakeComplex(SG_MAKE_INT(0), Sg_MakeRational(nume, deno));
3175 return Sg_MakeRational(nume, deno);
3176 }
3177 if (comp) return Sg_MakeComplex(Sg_MakeFlonum(0.0), Sg_Div(nume, deno));
3178 return Sg_Div(nume, deno);
3179 }
3180 if (SG_FLONUMP(obj)) {
3181 double s = SG_FLONUM_VALUE(obj);
3182 if (s < 0.0) return Sg_MakeComplex(Sg_MakeFlonum(0.0), Sg_MakeFlonum(sqrt(-s)));
3183 return Sg_MakeFlonum(sqrt(s));
3184 }
3185 if (SG_COMPLEXP(obj)) {
3186 double real, imag, m, x, y, s;
3187 SgComplex *cn = SG_COMPLEX(obj);
3188 if (Sg_ExactP(obj)) {
3189 SgObject m = Sg_Magnitude(cn);
3190 SgObject x = Sg_Div(Sg_Add(cn->real, m), SG_MAKE_INT(2));
3191 SgObject y = Sg_Div(cn->imag, SG_MAKE_INT(2));
3192 SgObject s = Sg_Sqrt(Sg_Div(m, Sg_Add(Sg_Mul(x, x), Sg_Mul(y, y))));
3193 return oprtr_norm_complex(Sg_Mul(x, s), Sg_Mul(y, s));
3194 }
3195 real = Sg_GetDouble(cn->real);
3196 imag = Sg_GetDouble(cn->imag);
3197 m = sqrt(real * real + imag * imag);
3198 x = (real + m) / 2;
3199 y = imag / 2;
3200 s = sqrt(m / (x * x + y * y));
3201 Sg_MakeComplex(Sg_MakeFlonum(x * s), Sg_MakeFlonum(y * s));
3202 }
3203 wte(SG_INTERN("sqrt"), "number", obj);
3204 return SG_UNDEF; /* dummy */
3205 }
3206
exact_integer_sqrt(SgObject k)3207 static inline SgObject exact_integer_sqrt(SgObject k)
3208 {
3209 SgObject ik = Sg_Sqrt(k);
3210 if (Sg_FiniteP(ik)) {
3211 return Sg_Exact(Sg_Round(ik, SG_ROUND_FLOOR));
3212 } else {
3213 long len = Sg_IntegerLength(k);
3214 SgObject quo = Sg_Quotient(SG_MAKE_INT(len), SG_MAKE_INT(2), NULL);
3215 ASSERT(SG_INTP(quo));
3216 return Sg_Ash(SG_MAKE_INT(1), SG_INT_VALUE(quo));
3217 }
3218 }
3219
Sg_ExactIntegerSqrt(SgObject k)3220 SgObject Sg_ExactIntegerSqrt(SgObject k)
3221 {
3222 double d;
3223
3224 if (!SG_EXACT_INTP(k)) {
3225 wte(SG_INTERN("exact-integer-sqrt"), "exact integer", k);
3226 }
3227
3228 d = Sg_GetDouble(k);
3229 if (d < iexpt_2n53) {
3230 double t = floor(sqrt(d));
3231 SgObject s = Sg_Exact(Sg_MakeFlonum(t));
3232 return Sg_Values2(s, Sg_Sub(k, Sg_Mul(s, s)));
3233 } else {
3234 SgObject s = exact_integer_sqrt(k);
3235 while (TRUE) {
3236 SgObject s2 = Sg_Mul(s, s);
3237 if (Sg_NumCmp(k, s2) < 0) {
3238 s = Sg_Quotient(Sg_Add(s2, k), Sg_Mul(SG_MAKE_INT(2), s), NULL);
3239 continue;
3240 } else {
3241 SgObject s2p = Sg_Add(Sg_Add(s2, Sg_Mul(SG_MAKE_INT(2), s)),
3242 SG_MAKE_INT(1));
3243 if (Sg_NumCmp(k, s2p) < 0) {
3244 return Sg_Values2(s, Sg_Sub(k, s2));
3245 } else {
3246 s = Sg_Quotient(Sg_Add(s2, k), Sg_Mul(SG_MAKE_INT(2), s), NULL);
3247 continue;
3248 }
3249 }
3250 }
3251 }
3252 }
3253
Sg_Sign(SgObject obj)3254 int Sg_Sign(SgObject obj)
3255 {
3256 int r = 0;
3257 if (SG_INTP(obj)) {
3258 long v = SG_INT_VALUE(obj);
3259 if (v > 0) r = 1;
3260 else if (v < 0) r = -1;
3261 } else if (SG_BIGNUMP(obj)) {
3262 r = SG_BIGNUM_GET_SIGN(obj);
3263 } else if (SG_FLONUMP(obj)) {
3264 double v = SG_FLONUM_VALUE(obj);
3265 if (v != 0.0) {
3266 r = v > 0.0 ? 1 : -1;
3267 }
3268 } else if (SG_RATIONALP(obj)) {
3269 return Sg_Sign(SG_RATIONAL(obj)->numerator);
3270 } else {
3271 wte(SG_INTERN("sign"), "real number", obj);
3272 }
3273 return r;
3274 }
3275
gcd_fixfix(unsigned long x,unsigned long y)3276 static inline unsigned long gcd_fixfix(unsigned long x, unsigned long y)
3277 {
3278 while (y > 0) {
3279 unsigned long r = x % y;
3280 x = y;
3281 y = r;
3282 }
3283 return x;
3284 }
3285
gcd_floflo(double x,double y)3286 static inline double gcd_floflo(double x, double y)
3287 {
3288 if (x < 0) x = -x;
3289 if (y < 0) y = -y;
3290 if (x < y) { double t = x; x = y; y = t; }
3291
3292 while (y > 0.0) {
3293 double r = fmod(x, y);
3294 x = y;
3295 y = r;
3296 }
3297 return x;
3298 }
3299
gcd_bigfix(SgBignum * x,unsigned long y)3300 static unsigned long gcd_bigfix(SgBignum *x, unsigned long y)
3301 {
3302 long rem;
3303 Sg_BignumDivSI(x, (signed long)y, &rem);
3304 if (rem < 0) rem = -rem;
3305 return gcd_fixfix(y, (unsigned long)rem);
3306 }
3307
Sg_Gcd(SgObject x,SgObject y)3308 SgObject Sg_Gcd(SgObject x, SgObject y)
3309 {
3310 int ox = FALSE, oy = FALSE;
3311 long ix, iy;
3312 unsigned long ux, uy, ur;
3313 if (!Sg_IntegerP(x)) {
3314 wte(SG_INTERN("gcd"), "integer", x);
3315 }
3316 if (!Sg_IntegerP(y)) {
3317 wte(SG_INTERN("gcd"), "integer", y);
3318 }
3319
3320 if (SG_FLONUMP(x) || SG_FLONUMP(y)) {
3321 return Sg_MakeFlonum(gcd_floflo(Sg_GetDouble(x), Sg_GetDouble(y)));
3322 }
3323
3324 if (SG_MAKE_INT(0) == x) return y;
3325 if (SG_MAKE_INT(0) == y) return x;
3326
3327 ix = Sg_GetIntegerClamp(x, SG_CLAMP_NONE, &ox);
3328 iy = Sg_GetIntegerClamp(y, SG_CLAMP_NONE, &oy);
3329
3330 if (!ox && !oy) {
3331 ux = (ix < 0) ? - ix : ix;
3332 uy = (iy < 0) ? - iy : iy;
3333 if (ux >= uy) {
3334 ur = gcd_fixfix(ux, uy);
3335 } else {
3336 ur = gcd_fixfix(uy, ux);
3337 }
3338 return Sg_MakeIntegerU(ur);
3339 }
3340
3341 if (!oy && iy != LONG_MIN) {
3342 ASSERT(SG_BIGNUMP(x));
3343 uy = (iy < 0) ? -iy : iy;
3344 ur = gcd_bigfix(SG_BIGNUM(x), uy);
3345 return Sg_MakeIntegerU(ur);
3346 }
3347
3348 if (!ox && ix != LONG_MIN) {
3349 ASSERT(SG_BIGNUMP(y));
3350 ux = (ix < 0) ? -ix : ix;
3351 ur = gcd_bigfix(SG_BIGNUM(y), ux);
3352 return Sg_MakeIntegerU(ur);
3353 }
3354 ASSERT(SG_BIGNUMP(x) && SG_BIGNUMP(y));
3355 return Sg_BignumGcd(SG_BIGNUM(x), SG_BIGNUM(y));
3356 /*
3357 x = Sg_Abs(x);
3358 y = Sg_Abs(y);
3359 if (Sg_NumCmp(x, y) < 0) {SgObject t = x; x = y; y = t;}
3360
3361 while (!SG_EQ(SG_MAKE_INT(0), y)) {
3362 SgObject r = Sg_Modulo(x, y, TRUE);
3363 x = y;
3364 y = r;
3365 }
3366 return x;
3367 */
3368 }
3369
Sg_Magnitude(SgObject z)3370 SgObject Sg_Magnitude(SgObject z)
3371 {
3372 if (SG_COMPLEXP(z)) {
3373 SgComplex *cn = SG_COMPLEX(z);
3374 if (Sg_ExactP(cn)) {
3375 if (Sg_ZeroP(cn->real)) return Sg_Magnitude(cn->imag);
3376 if (Sg_ZeroP(cn->imag)) return Sg_Magnitude(cn->real);
3377 return Sg_Sqrt(Sg_Add(Sg_Mul(cn->real, cn->real),
3378 Sg_Mul(cn->imag, cn->imag)));
3379 } else {
3380 double real = Sg_GetDouble(cn->real);
3381 double imag = Sg_GetDouble(cn->imag);
3382 double m;
3383 if (isinf(real) || isinf(imag)) return SG_POSITIVE_INFINITY;
3384 m = sqrt(real * real + imag * imag);
3385 if (m < DBL_EPSILON || isinf(m)) return Sg_MakeFlonum(imag / sin(atan2(imag, real)));
3386 return Sg_MakeFlonum(m);
3387 }
3388 }
3389 if (SG_REALP(z)) {
3390 if (Sg_NegativeP(z)) return Sg_Negate(z);
3391 return z;
3392 }
3393 wte(SG_INTERN("magnitude"), "number", z);
3394 return SG_UNDEF; /* dummy */
3395 }
3396
Sg_Angle(SgObject obj)3397 SgObject Sg_Angle(SgObject obj)
3398 {
3399 if (SG_COMPLEXP(obj)) {
3400 double real = Sg_GetDouble(SG_COMPLEX(obj)->real);
3401 double imag = Sg_GetDouble(SG_COMPLEX(obj)->imag);
3402 return Sg_MakeFlonum(atan2(imag, real));
3403 }
3404 if (SG_REALP(obj)) {
3405 if (Sg_NegativeP(obj)) return Sg_MakeFlonum(atan2(0.0, -1.0)); /* pi */
3406 if (SG_FLONUMP(obj)) return Sg_MakeFlonum(0.0);
3407 return SG_MAKE_INT(0);
3408 }
3409 wte(SG_INTERN("angle"), "number", obj);
3410 return SG_UNDEF; /* dummy */
3411 }
3412
log_handle_inf(SgObject n)3413 static inline SgObject log_handle_inf(SgObject n)
3414 {
3415 /* passing argument must be always bignum */
3416 /* Sg_BignumSqrtApprox returns approximate
3417 sqrt and returns always exact number */
3418 SgObject r = Sg_BignumSqrtApprox(n);
3419 return Sg_Add(Sg_Log(r), Sg_Log(r));
3420 }
3421
Sg_Log(SgObject obj)3422 SgObject Sg_Log(SgObject obj)
3423 {
3424 double real;
3425 double imag;
3426
3427 if (SG_INTP(obj)) {
3428 long value = SG_INT_VALUE(obj);
3429 if (value > 0) {
3430 if (value == 1) return SG_MAKE_INT(0);
3431 return Sg_MakeFlonum(log((double)value));
3432 }
3433 real = value;
3434 return Sg_MakeComplex(Sg_MakeFlonum(0.5 * log(real * real)),
3435 Sg_MakeFlonum(atan2(0.0, real)));
3436 }
3437 if (SG_COMPLEXP(obj)) {
3438 SgObject or = SG_COMPLEX(obj)->real, oi = SG_COMPLEX(obj)->imag;
3439
3440 if (!SG_BIGNUMP(or) && !SG_BIGNUMP(oi)) {
3441 real = Sg_GetDouble(or);
3442 imag = Sg_GetDouble(oi);
3443 return Sg_MakeComplex(Sg_MakeFlonum(0.5 * log(real * real + imag * imag)),
3444 Sg_MakeFlonum(atan2(imag, real)));
3445 } else {
3446 /* do the same but with Scheme function */
3447 SgObject r2 = Sg_Mul(or, or);
3448 SgObject i2 = Sg_Mul(oi, oi);
3449 static SgObject HALF = NULL;
3450 if (!HALF) HALF = Sg_MakeFlonum(0.5);
3451 return Sg_MakeComplex(Sg_Mul(HALF, Sg_Log(Sg_Add(r2, i2))),
3452 Sg_Atan2(oi, or));
3453 }
3454 }
3455 if (SG_REALP(obj)) {
3456 real = Sg_GetDouble(obj);
3457 if (isinf(real) && SG_BIGNUMP(obj)) return log_handle_inf(obj);
3458 if (real > 0) return Sg_MakeFlonum(log(real));
3459 imag = atan2(0.0, real);
3460 if (imag == 0.0) return Sg_MakeFlonum(0.5 * log(real * real));
3461 return Sg_MakeComplex(Sg_MakeFlonum(0.5 * log(real * real)),
3462 Sg_MakeFlonum(imag));
3463 }
3464 wte(SG_INTERN("log"), "number", obj);
3465 return SG_UNDEF; /* dummy */
3466 }
3467
Sg_MinMax(SgObject arg0,SgObject args,SgObject * min,SgObject * max)3468 void Sg_MinMax(SgObject arg0, SgObject args, SgObject *min, SgObject *max)
3469 {
3470 #define EXACTP(o) (SG_EXACT_INTP(o) || SG_RATIONALP(o))
3471 int inexact = !EXACTP(arg0);
3472 SgObject mi = arg0;
3473 SgObject ma = arg0;
3474
3475 /* nan check */
3476 if (Sg_NanP(arg0)) {
3477 if (min) *min = arg0;
3478 if (max) *max = arg0;
3479 return;
3480 }
3481
3482 for (;;) {
3483 if (!SG_REALP(arg0))
3484 wte(SG_INTERN("min/max"), "real number", arg0);
3485 if (SG_NULLP(args)) {
3486 if (min) {
3487 if (inexact && EXACTP(mi)) {
3488 *min = Sg_Inexact(mi);
3489 } else {
3490 *min = mi;
3491 }
3492 }
3493 if (max) {
3494 if (inexact && EXACTP(ma)) {
3495 *max = Sg_Inexact(ma);
3496 } else {
3497 *max = ma;
3498 }
3499 }
3500 return;
3501 }
3502 /* nan check */
3503 if (Sg_NanP(SG_CAR(args))) {
3504 if (min) *min = SG_CAR(args);
3505 if (max) *max = SG_CAR(args);
3506 return;
3507 }
3508 if (!EXACTP(SG_CAR(args))) inexact = TRUE;
3509 if (min && Sg_NumCmp(mi, SG_CAR(args)) > 0) {
3510 mi = SG_CAR(args);
3511 }
3512 if (max && Sg_NumCmp(ma, SG_CAR(args)) < 0) {
3513 ma = SG_CAR(args);
3514 }
3515 args = SG_CDR(args);
3516 }
3517 #undef EXACTP
3518 }
3519
Sg_IntegerDiv(SgObject x,SgObject y)3520 SgObject Sg_IntegerDiv(SgObject x, SgObject y)
3521 {
3522 if (!SG_REALP(x)) wte(SG_INTERN("div"), "real number", x);
3523 if (!SG_REALP(y)) wte(SG_INTERN("div"), "real number", y);
3524 if (!Sg_FiniteP(x) || Sg_NanP(x)) {
3525 Sg_AssertionViolation(SG_INTERN("div"),
3526 SG_MAKE_STRING("dividend must be neither infinite nor NaN"), x);
3527 }
3528 if (Sg_ZeroP(y)) {
3529 Sg_AssertionViolation(SG_INTERN("div"),
3530 SG_MAKE_STRING("undefined for 0"), y);
3531 }
3532
3533 if (SG_INTP(x)) {
3534 if (SG_INTP(y)) {
3535 long xx = SG_INT_VALUE(x);
3536 long yy = SG_INT_VALUE(y);
3537 long div;
3538 if (xx == 0) {
3539 div = 0;
3540 } else if (xx > 0) {
3541 div = xx / yy;
3542 } else if (yy > 0) {
3543 div = (xx - yy + 1) / yy;
3544 } else {
3545 div = (xx + yy + 1) / yy;
3546 }
3547 return Sg_MakeInteger(div);
3548 }
3549 }
3550 if (SG_FLONUMP(x) || SG_FLONUMP(y)) {
3551 double xx = Sg_GetDouble(x);
3552 double yy = Sg_GetDouble(y);
3553 return Sg_MakeFlonum((yy > 0.0) ? floor(xx / yy) : - floor(xx / -yy));
3554 }
3555 if (Sg_PositiveP(y)) return Sg_Round(Sg_Div(x, y), SG_ROUND_FLOOR);
3556 return Sg_Negate(Sg_Round(Sg_Div(x, Sg_Negate(y)), SG_ROUND_FLOOR));
3557 }
3558
Sg_IntegerDiv0(SgObject x,SgObject y)3559 SgObject Sg_IntegerDiv0(SgObject x, SgObject y)
3560 {
3561 SgObject div, mod;
3562 if (!SG_REALP(x)) wte(SG_INTERN("div0"), "real number", x);
3563 if (!SG_REALP(y)) wte(SG_INTERN("div0"), "real number", y);
3564 if (!Sg_FiniteP(x) || Sg_NanP(x)) {
3565 Sg_AssertionViolation(SG_INTERN("div0"),
3566 SG_MAKE_STRING("dividend must be neither infinite nor NaN"), x);
3567 }
3568 if (Sg_ZeroP(y)) {
3569 Sg_AssertionViolation(SG_INTERN("div0"),
3570 SG_MAKE_STRING("undefined for 0"), y);
3571 }
3572
3573 div = Sg_IntegerDiv(x, y);
3574 mod = Sg_Sub(x, Sg_Mul(div, y));
3575 if (Sg_NumCmp(mod, Sg_Magnitude(Sg_Div(y, SG_MAKE_INT(2)))) < 0) return div;
3576 if (Sg_PositiveP(y)) return Sg_Add(div, SG_MAKE_INT(1));
3577 return Sg_Sub(div, SG_MAKE_INT(1));
3578
3579 }
3580
Sg_IntegerMod(SgObject x,SgObject y)3581 SgObject Sg_IntegerMod(SgObject x, SgObject y)
3582 {
3583 /* only for bignum */
3584 if (SG_EXACT_INTP(x) && SG_EXACT_INTP(y)) {
3585 /* for performance we need to use Sg_Modulo */
3586 int xsign = Sg_Sign(x), ysign = Sg_Sign(y);
3587 if (ysign == 0) goto err;
3588 if (xsign > 0) {
3589 if (ysign > 0) return Sg_Modulo(x, y, TRUE);
3590 else return Sg_Modulo(x, y, TRUE);
3591 } else if (xsign < 0) {
3592 SgObject r = Sg_Modulo(x, y, TRUE);
3593 if (SG_EQ(SG_MAKE_INT(0), r)) return r;
3594 if (ysign > 0) return Sg_Add(r, y);
3595 else return Sg_Sub(r, y);
3596 } else if (xsign == 0) {
3597 return SG_MAKE_INT(0);
3598 }
3599 } else {
3600 SgObject d = Sg_IntegerDiv(x, y);
3601 SgObject m = Sg_Mul(d, y);
3602 return Sg_Sub(x, m);
3603 }
3604
3605 err:
3606 /* y was 0 */
3607 wte(SG_INTERN("mod"), "non zero number", y);
3608 return SG_UNDEF; /* dummy */
3609 }
3610
Sg_IntegerMod0(SgObject x,SgObject y)3611 SgObject Sg_IntegerMod0(SgObject x, SgObject y)
3612 {
3613 /* only for bignum */
3614 if (SG_EXACT_INTP(x) && SG_EXACT_INTP(y)) {
3615 /* for performance we need to use Sg_Modulo */
3616 int xsign = Sg_Sign(x), ysign = Sg_Sign(y);
3617 SgObject r;
3618 if (ysign == 0) goto err;
3619 r = Sg_Modulo(x, y, TRUE);
3620 if (xsign >= 0) {
3621 if (ysign > 0) {
3622 if (Sg_NumCmp(Sg_Mul(r, SG_MAKE_INT(2)), y) >= 0)
3623 return Sg_Sub(r, y);
3624 else return r;
3625 } else {
3626 if (Sg_NumCmp(Sg_Mul(r, SG_MAKE_INT(2)), Sg_Negate(y)) >= 0)
3627 return Sg_Add(r, y);
3628 else return r;
3629 }
3630 } else {
3631 if (ysign > 0) {
3632 if (Sg_NumCmp(Sg_Mul(r, SG_MAKE_INT(2)), Sg_Negate(y)) < 0)
3633 return Sg_Add(r, y);
3634 else return r;
3635 } else {
3636 if (Sg_NumCmp(Sg_Mul(r, SG_MAKE_INT(2)), y) < 0)
3637 return Sg_Sub(r, y);
3638 else return r;
3639 }
3640 }
3641 } else {
3642 SgObject d = Sg_IntegerDiv0(x, y);
3643 SgObject m = Sg_Mul(d, y);
3644 return Sg_Sub(x, m);
3645 }
3646 err:
3647 /* y was 0 */
3648 wte(SG_INTERN("mod0"), "non zero number", y);
3649 return SG_UNDEF; /* dummy */
3650 }
3651
3652 /* for better performance */
Sg_ModInverse(SgObject x,SgObject m)3653 SgObject Sg_ModInverse(SgObject x, SgObject m)
3654 {
3655 SgObject u1, u3, v1, v3;
3656 int sign = 1;
3657 if (!SG_EXACT_INTP(x)) wte(SG_INTERN("mod-inverse"), "exact integer", x);
3658 if (!SG_EXACT_INTP(m)) wte(SG_INTERN("mod-inverse"), "exact integer", m);
3659 if (Sg_Sign(m) != 1) {
3660 wte(SG_INTERN("mod-inverse"), "positive number", m);
3661 }
3662
3663 if (SG_BIGNUMP(x) && SG_BIGNUMP(m)) {
3664 return Sg_BignumModInverse(x, m);
3665 }
3666 u1 = SG_MAKE_INT(1);
3667 u3 = x;
3668 v1 = SG_MAKE_INT(0);
3669 v3 = m;
3670
3671 while (!SG_EQ(v3, SG_MAKE_INT(0))) {
3672 SgObject t3, w, t1, q;
3673 t3 = Sg_IntegerMod(u3, v3);
3674 q = Sg_Quotient(u3, v3, NULL);
3675 w = Sg_Mul(q, v1);
3676 t1 = Sg_Add(u1, w);
3677 u1 = v1; v1 = t1; u3 = v3; v3 = t3;
3678 sign = -sign;
3679 }
3680 if (sign < 0) {
3681 return Sg_Sub(m, u1);
3682 } else {
3683 return u1;
3684 }
3685 }
3686
Sg_ModExpt(SgObject x,SgObject e,SgObject m)3687 SgObject Sg_ModExpt(SgObject x, SgObject e, SgObject m)
3688 {
3689 SgObject y;
3690 int invertp = FALSE;
3691 if (!SG_EXACT_INTP(x)) wte(SG_INTERN("mod-expt"), "exact integer", x);
3692 if (!SG_EXACT_INTP(e)) wte(SG_INTERN("mod-expt"), "exact integer", e);
3693 if (!SG_EXACT_INTP(m)) wte(SG_INTERN("mod-expt"), "exact integer", m);
3694 if (Sg_Sign(m) <= 0) {
3695 wte(SG_INTERN("mod-inverse"), "positive number", m);
3696 }
3697
3698 /* TODO handle more efficiently */
3699 if (SG_INTP(x)) {
3700 if (SG_INTP(e)) {
3701 if (SG_INTP(m)) {
3702 /* can be done in C */
3703 #if LONG_MAX > 1UL << 32
3704 goto small_entry;
3705 #else
3706 int64_t yy = 1, n = SG_INT_VALUE(e), d = SG_INT_VALUE(m);
3707 int64_t xx = SG_INT_VALUE(x);
3708 if (n < 0) {
3709 n = -n;
3710 invertp = TRUE;
3711 }
3712 while (n > 0) {
3713 if (n % 2) {
3714 yy = (yy * xx) % d;
3715 }
3716 n >>= 1;
3717 if (n > 0) {
3718 xx = (xx * xx) % d;
3719 }
3720 }
3721 y = Sg_MakeIntegerFromS64(yy);
3722 if (invertp) {
3723 return Sg_ModInverse(y, m);
3724 } else {
3725 return y;
3726 }
3727 #endif
3728 }
3729 e = Sg_MakeBignumFromSI(SG_INT_VALUE(e));
3730 }
3731 x = Sg_MakeBignumFromSI(SG_INT_VALUE(x));
3732 if (SG_INTP(m)) {
3733 m = Sg_MakeBignumFromSI(SG_INT_VALUE(m));
3734 }
3735 } else if (SG_INTP(e)) {
3736 if (SG_BIGNUMP(x) && SG_BIGNUMP(m)) {
3737 /* if both are bignum, then the cost of mod and mul would be
3738 more expensive than converting e to bignum */
3739 e = Sg_MakeBignumFromSI(SG_INT_VALUE(e));
3740 } else {
3741 long n;
3742 #if LONG_MAX > 1UL << 32
3743 small_entry:
3744 #endif
3745 n = SG_INT_VALUE(e);
3746 y = SG_MAKE_INT(1);
3747 if (n < 0) {
3748 invertp = TRUE;
3749 n = -n;
3750 }
3751 while (n > 0) {
3752 if (n % 2) {
3753 y = Sg_IntegerMod(Sg_Mul(y, x), m);
3754 }
3755 n >>= 1;
3756 if (n > 0) {
3757 x = Sg_IntegerMod(Sg_Mul(x, x), m);
3758 }
3759 }
3760 if (invertp) {
3761 return Sg_ModInverse(y, m);
3762 } else {
3763 return y;
3764 }
3765 }
3766 } else if (SG_INTP(m)) {
3767 /* both x and e are bignum */
3768 m = Sg_MakeBignumFromSI(SG_INT_VALUE(m));
3769 }
3770 if (!(SG_BIGNUMP(x) && SG_BIGNUMP(e) && SG_BIGNUMP(m))) {
3771 Sg_Error(UC("[internal] bignum required. %S %S %S"), x, e, m);
3772 }
3773 return Sg_BignumModExpt(SG_BIGNUM(x), SG_BIGNUM(e), SG_BIGNUM(m));
3774 }
3775
Sg_Square(SgObject x)3776 SgObject Sg_Square(SgObject x)
3777 {
3778 if (SG_BIGNUMP(x)) {
3779 /* a bit of optimisation */
3780 return Sg_BignumSquare(SG_BIGNUM(x));
3781 }
3782 /* TODO rational and complex but it's rare i think...*/
3783 return Sg_Mul(x, x);
3784 }
3785
Sg_Round(SgObject num,int mode)3786 SgObject Sg_Round(SgObject num, int mode)
3787 {
3788 if (SG_EXACT_INTP(num)) return num;
3789 if (SG_RATIONALP(num)) {
3790 int offset = 0;
3791 SgObject rem;
3792 SgObject quot = Sg_Quotient(SG_RATIONAL(num)->numerator,
3793 SG_RATIONAL(num)->denominator, &rem);
3794 /* this shouldn't happen but just in case */
3795 if (SG_EQ(rem, SG_MAKE_INT(0))) return quot;
3796 switch (mode) {
3797 case SG_ROUND_FLOOR:
3798 offset = (Sg_Sign(num) < 0) ? -1 : 0;
3799 break;
3800 case SG_ROUND_CEIL:
3801 offset = (Sg_Sign(num) < 0) ? 0 : 1;
3802 break;
3803 case SG_ROUND_TRUNC:
3804 offset = 0;
3805 break;
3806 case SG_ROUND_ROUND: {
3807 SgObject rem2 = Sg_Mul(Sg_Abs(rem), SG_MAKE_INT(2));
3808 int cmp = Sg_NumCmp(SG_RATIONAL(num)->denominator, rem2);
3809 if (cmp > 0) {
3810 /* NUM is closer to zero than halfway */
3811 offset = 0;
3812 } else if (cmp < 0) {
3813 /* NUM is further from zero than halfway */
3814 offset = (Sg_Sign(num) < 0) ? -1 : 1;
3815 } else {
3816 /* NUM is exactly the halfway. We round to even */
3817 if (Sg_OddP(quot)) {
3818 offset = (Sg_Sign(num) < 0) ? -1 : 1;
3819 } else {
3820 offset = 0;
3821 }
3822 }
3823 break;
3824 }
3825 default:
3826 Sg_Panic("something screwed up");
3827 }
3828 if (offset == 0) return quot;
3829 else return Sg_Add(quot, SG_MAKE_INT(offset));
3830 }
3831 if (SG_FLONUMP(num)) {
3832 double r = 0.0, v;
3833 v = SG_FLONUM_VALUE(num);
3834 switch (mode) {
3835 case SG_ROUND_FLOOR: r = floor(v); break;
3836 case SG_ROUND_CEIL: r = ceil(v); break;
3837 case SG_ROUND_TRUNC:
3838 r = (v < 0.0) ? ceil(v) : floor(v); break;
3839 case SG_ROUND_ROUND: r = roundeven(v); break;
3840 default:
3841 Sg_Panic("something screwed up");
3842 }
3843 return Sg_MakeFlonum(r);
3844 }
3845 wte(SG_INTERN("round"), "real number", num);
3846 return SG_UNDEF;
3847 }
3848
3849 /* from gauche */
numcmp3(SgObject x,SgObject d,SgObject y)3850 static inline int numcmp3(SgObject x, SgObject d, SgObject y)
3851 {
3852 if (SG_INTP(x) && SG_INTP(d) && SG_INTP(y)) {
3853 long xd = SG_INT_VALUE(x) + SG_INT_VALUE(d);
3854 if (xd < SG_INT_VALUE(y)) return -1;
3855 if (xd > SG_INT_VALUE(y)) return 1;
3856 else return 0;
3857 } else {
3858 SgObject bx = SG_BIGNUMP(x) ? x : Sg_MakeBignumFromSI(SG_INT_VALUE(x));
3859 SgObject bd = SG_BIGNUMP(d) ? d : Sg_MakeBignumFromSI(SG_INT_VALUE(d));
3860 SgObject by = SG_BIGNUMP(y) ? y : Sg_MakeBignumFromSI(SG_INT_VALUE(y));
3861 return Sg_BignumCmp3U(SG_BIGNUM(bx), SG_BIGNUM(bd), SG_BIGNUM(by));
3862 }
3863 }
3864
3865 #ifdef _MSC_VER
3866 /* to shut compiler */
3867 # define strcpy_ strcpy_s
3868 #else
3869 /* this works fine anyway */
3870 # define strcpy_(buf, len, s) strcpy(buf, s)
3871 #endif
3872
double_print(char * buf,int buflen,double val,int plus_sign)3873 static void double_print(char *buf, int buflen, double val, int plus_sign)
3874 {
3875 SgObject f;
3876 int exp, sign;
3877 f = Sg_DecodeFlonum(val, &exp, &sign);
3878 if (val == 0.0) {
3879 if (plus_sign) strcpy_(buf, buflen, "+0.0");
3880 else if (sign < 0) strcpy_(buf, buflen, "-0.0");
3881 else strcpy_(buf, buflen, "0.0");
3882 return;
3883 } else if (isinf(val)) {
3884 if (sign < 0) strcpy_(buf, buflen, "-inf.0");
3885 else strcpy_(buf, buflen, "+inf.0");
3886 return;
3887 } else if (isnan(val)) {
3888 strcpy_(buf, buflen, "+nan.0");
3889 return;
3890 }
3891 if (sign < 0) *buf++ = '-', buflen--;
3892 else if (plus_sign) *buf++ = '+', buflen--;
3893
3894 {
3895 SgObject r, s, mp, mm, q;
3896 int est, tc1, tc2, tc3, digs, point, round;
3897 int mp2 = FALSE, fixup = FALSE;
3898 if (sign < 0) val = -val;
3899 /* initialize r, s, m+ and m- */
3900
3901 round = !Sg_OddP(f);
3902 if (exp >= 0) {
3903 SgObject be = Sg_Ash(SG_MAKE_INT(1), exp);
3904 if (Sg_NumCmp(f, SG_2_52) != 0) {
3905 r = Sg_Ash(f, exp + 1);
3906 s = SG_MAKE_INT(2);
3907 mp2 = FALSE;
3908 mm = be;
3909 } else {
3910 r = Sg_Ash(f, exp + 2);
3911 s = SG_MAKE_INT(4);
3912 mp2 = TRUE;
3913 mm = be;
3914 }
3915 } else {
3916 if (exp == -1023 || Sg_NumCmp(f, SG_2_52) != 0) {
3917 r = Sg_Ash(f, 1);
3918 s = Sg_Ash(SG_MAKE_INT(1), -exp + 1);
3919 mp2 = FALSE;
3920 mm = SG_MAKE_INT(1);
3921 } else {
3922 r = Sg_Ash(f, 2);
3923 s = Sg_Ash(SG_MAKE_INT(1), -exp + 2);
3924 mp2 = TRUE;
3925 mm = SG_MAKE_INT(1);
3926 }
3927 }
3928 /* Sg_Printf(Sg_StandardErrorPort(), UC("exp=%d, round=%d, r=%S, s=%S, mp=%d, mm=%S\n"), */
3929 /* exp, round, r, s, mp, mm); */
3930
3931 /* estimate scale */
3932 est = (int)ceil(log10(val) - 0.1);
3933 if (est >= 0) {
3934 s = Sg_Mul(s, Sg_Expt(SG_MAKE_INT(10), SG_MAKE_INT(est)));
3935 } else {
3936 SgObject scale = Sg_Expt(SG_MAKE_INT(10), SG_MAKE_INT(-est));
3937 r = Sg_Mul(r, scale);
3938 mm = Sg_Mul(mm, scale);
3939 }
3940 /* fixup. avoid calculating m+ for obvious case. */
3941 if (Sg_NumCmp(r, s) >= 0) {
3942 fixup = TRUE;
3943 } else {
3944 mp = (mp2 ? Sg_Ash(mm, 1) : mm);
3945 if (round) {
3946 fixup = (numcmp3(r, mp, s) >= 0);
3947 } else {
3948 fixup = (numcmp3(r, mp, s) > 0);
3949 }
3950 }
3951 if (fixup) {
3952 s = Sg_Mul(s, SG_MAKE_INT(10));
3953 est++;
3954 }
3955 /* Sg_Printf(Sg_StandardErrorPort(), UC("est=%d, r=%S, s=%S, mp=%S, mm=%S\n"), */
3956 /* est, r, s, mp, mm); */
3957
3958 /* determine position of decimal point. */
3959 if (est < 10 && est > -3) {
3960 point = est; est = 1;
3961 } else {
3962 point = 1;
3963 }
3964
3965 /* generate */
3966 if (point <= 0) {
3967 *buf++ = '0'; buflen--;
3968 *buf++ = '.', buflen--;
3969 for (digs = point; digs < 0 && buflen > 5; digs++) {
3970 *buf++ = '0'; buflen--;
3971 }
3972 }
3973 for (digs = 1; buflen > 5; digs++) {
3974 SgObject r10 = Sg_Mul(r, SG_MAKE_INT(10));
3975 q = Sg_Quotient(r10, s, &r);
3976 mm = Sg_Mul(mm, SG_MAKE_INT(10));
3977 mp = (mp2 ? Sg_Ash(mm, 1) : mm);
3978
3979 /* Sg_Printf(Sg_StandardErrorPort(), UC("q=%S, r=%S, mp=%S, mm=%S\n"), */
3980 /* q, r, mp, mm); */
3981
3982 ASSERT(SG_INTP(q));
3983 if (round) {
3984 tc1 = (Sg_NumCmp(r, mm) <= 0);
3985 tc2 = (numcmp3(r, mp, s) >= 0);
3986 } else {
3987 tc1 = (Sg_NumCmp(r, mm) < 0);
3988 tc2 = (numcmp3(r, mp, s) > 0);
3989 }
3990
3991 if (!tc1) {
3992 if (!tc2) {
3993 *buf++ = (char)SG_INT_VALUE(q) + '0';
3994 buflen--;
3995 if (digs == point) *buf++ = '.', buflen--;
3996 continue;
3997 } else {
3998 *buf++ = (char)SG_INT_VALUE(q) + '1';
3999 buflen--;
4000 break;
4001 }
4002 } else {
4003 if (!tc2) {
4004 *buf++ = (char)SG_INT_VALUE(q) + '0';
4005 buflen--;
4006 break;
4007 } else {
4008 tc3 = numcmp3(r, r, s);
4009 if ((round && tc3 <= 0) || (!round && tc3 < 0)) {
4010 *buf++ = (char)SG_INT_VALUE(q) + '0';
4011 buflen--;
4012 break;
4013 } else {
4014 *buf++ = (char)SG_INT_VALUE(q) + '1';
4015 buflen--;
4016 break;
4017 }
4018 }
4019 }
4020 }
4021 if (digs <= point) {
4022 for (; digs < point && buflen > 5; digs++) {
4023 *buf++ = '0', buflen--;
4024 }
4025 *buf++ = '.', buflen--;
4026 *buf++ = '0', buflen--;
4027 }
4028 est--;
4029 if (est != 0) {
4030 *buf++ = 'e', buflen--;
4031 snprintf(buf, buflen, "%d", (int)est);
4032 } else {
4033 *buf++ = 0;
4034 }
4035 }
4036 }
4037
number_print(SgObject obj,SgPort * port,SgWriteContext * ctx)4038 static void number_print(SgObject obj, SgPort *port, SgWriteContext *ctx)
4039 {
4040 SgObject s = Sg_NumberToString(obj, 10, FALSE);
4041 Sg_Puts(port, SG_STRING(s));
4042 }
4043
4044 #define FLT_BUF 50
4045
Sg_NumberToString(SgObject obj,int radix,int use_upper)4046 SgObject Sg_NumberToString(SgObject obj, int radix, int use_upper)
4047 {
4048 SgObject r = SG_NIL;
4049 char buf[FLT_BUF];
4050
4051 if (SG_INTP(obj)) {
4052 char *pbuf = buf;
4053 long value = SG_INT_VALUE(obj);
4054 if (value < 0) {
4055 *pbuf++ = '-';
4056 value = -value;
4057 }
4058 if (radix == 10) {
4059 snprintf(pbuf, FLT_BUF - 1, "%ld" , value);
4060 } else if (radix == 16) {
4061 snprintf(pbuf, FLT_BUF - 1, (use_upper ? "%lX" : "%lx"), value);
4062 } else if (radix == 8) {
4063 snprintf(pbuf, FLT_BUF - 1, "%lo", value);
4064 } else {
4065 /* should use better way */
4066 r = Sg_BignumToString(SG_BIGNUM(Sg_MakeBignumFromSI(SG_INT_VALUE(obj))),
4067 radix, use_upper);
4068 }
4069 if (SG_NULLP(r)) r = Sg_MakeStringC(buf);
4070 } else if (SG_BIGNUMP(obj)) {
4071 r = Sg_BignumToString(SG_BIGNUM(obj), radix, use_upper);
4072 } else if (SG_FLONUMP(obj)) {
4073 double_print(buf, FLT_BUF, SG_FLONUM_VALUE(obj), FALSE);
4074 r = Sg_MakeStringC(buf);
4075 } else if (SG_RATIONALP(obj)) {
4076 SgObject nume, deno;
4077 nume = Sg_NumberToString(SG_RATIONAL(obj)->numerator, radix, use_upper);
4078 deno = Sg_NumberToString(SG_RATIONAL(obj)->denominator, radix, use_upper);
4079 r = Sg_StringAppend2(nume, SG_MAKE_STRING("/"));
4080 r = Sg_StringAppend2(r, deno);
4081 } else if (SG_COMPLEXP(obj)) {
4082 SgObject real, imag;
4083 int needPlus = FALSE;
4084 real = Sg_NumberToString(SG_COMPLEX(obj)->real, radix, use_upper);
4085 imag = Sg_NumberToString(SG_COMPLEX(obj)->imag, radix, use_upper);
4086 needPlus = (SG_STRING_VALUE_AT(imag, 0) != '+' &&
4087 SG_STRING_VALUE_AT(imag, 0) != '-');
4088 if (needPlus) {
4089 r = Sg_StringAppend(SG_LIST4(real, SG_MAKE_STRING("+"),
4090 imag, SG_MAKE_STRING("i")));
4091 } else {
4092 r = Sg_StringAppend(SG_LIST3(real, imag, SG_MAKE_STRING("i")));
4093 }
4094 } else {
4095 wte(SG_INTERN("number->string"), "number", obj);
4096 }
4097 return r;
4098 }
4099
4100 SgObject Sg__ConstObjes[SG_NUM_CONST_OBJS] = {SG_FALSE};
4101
4102 extern void Sg__InitBignum();
Sg__InitNumber()4103 void Sg__InitNumber()
4104 {
4105 /* VC does not have these macros. SUCKS!!
4106 Seems Solaris doesn't have it either...
4107 */
4108 #if !defined(INFINITY) || !defined(NAN)
4109 #undef INFINITY
4110 #undef NAN
4111 static double __d0 = 0.0;
4112 static double __d1 = 1.0;
4113 #define INFINITY (__d1/__d0)
4114 #define NAN (__d0/__d0)
4115 #endif
4116
4117 int radix, i;
4118 unsigned long n;
4119
4120 for (radix = RADIX_MIN; radix <= RADIX_MAX; radix++) {
4121 longlimit[radix - RADIX_MIN] =
4122 (unsigned long)floor((double)LONG_MAX / radix - radix);
4123 /* Find max D where R^(D+1)-1 <= LONG_MAX */
4124 for (i = 0, n = 1; ; i++, n *= radix) {
4125 if (n >= (unsigned long)(LONG_MAX / radix)) {
4126 longdigs[radix - RADIX_MIN] = i - 1;
4127 bigdig[radix - RADIX_MIN] = n;
4128 break;
4129 }
4130 }
4131 }
4132
4133 SG_2_52 = Sg_MakeBignumFromS64(iexpt_2n52);
4134
4135 #define INIT_CONST_FL(o, d) o=make_flonum(d)
4136 INIT_CONST_FL(SG_NAN, NAN);
4137 INIT_CONST_FL(SG_POSITIVE_INFINITY, INFINITY);
4138 INIT_CONST_FL(SG_NEGATIVE_INFINITY, -INFINITY);
4139 INIT_CONST_FL(SG_FL_POSITIVE_ZERO, 0.0);
4140 INIT_CONST_FL(SG_FL_NEGATIVE_ZERO, -0.0);
4141 INIT_CONST_FL(SG_FL_POSITIVE_ONE, 1.0);
4142 INIT_CONST_FL(SG_FL_NEGATIVE_ONE, -1.0);
4143
4144 Sg__InitBignum();
4145 }
4146
4147 /*
4148 end of file
4149 Local Variables:
4150 coding: utf-8-unix
4151 End:
4152 */
4153