Add a copyright year.
[bpt/guile.git] / libguile / numbers.c
1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005 Free Software Foundation, Inc.
2 *
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
5 *
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 */
21
22 \f
23 /* General assumptions:
24 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
25 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
26 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
27 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
28 * All objects satisfying SCM_FRACTIONP are never an integer.
29 */
30
31 /* TODO:
32
33 - see if special casing bignums and reals in integer-exponent when
34 possible (to use mpz_pow and mpf_pow_ui) is faster.
35
36 - look in to better short-circuiting of common cases in
37 integer-expt and elsewhere.
38
39 - see if direct mpz operations can help in ash and elsewhere.
40
41 */
42
43 /* tell glibc (2.3) to give prototype for C99 trunc() */
44 #define _GNU_SOURCE
45
46 #if HAVE_CONFIG_H
47 # include <config.h>
48 #endif
49
50 #include <math.h>
51 #include <ctype.h>
52 #include <string.h>
53
54 #include "libguile/_scm.h"
55 #include "libguile/feature.h"
56 #include "libguile/ports.h"
57 #include "libguile/root.h"
58 #include "libguile/smob.h"
59 #include "libguile/strings.h"
60
61 #include "libguile/validate.h"
62 #include "libguile/numbers.h"
63 #include "libguile/deprecation.h"
64
65 #include "libguile/eq.h"
66
67 #include "libguile/discouraged.h"
68
69 \f
70
71 /*
72 Wonder if this might be faster for some of our code? A switch on
73 the numtag would jump directly to the right case, and the
74 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
75
76 #define SCM_I_NUMTAG_NOTNUM 0
77 #define SCM_I_NUMTAG_INUM 1
78 #define SCM_I_NUMTAG_BIG scm_tc16_big
79 #define SCM_I_NUMTAG_REAL scm_tc16_real
80 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
81 #define SCM_I_NUMTAG(x) \
82 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
83 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
84 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
85 : SCM_I_NUMTAG_NOTNUM)))
86 */
87 /* the macro above will not work as is with fractions */
88
89
90 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
91
92 /* FLOBUFLEN is the maximum number of characters neccessary for the
93 * printed or scm_string representation of an inexact number.
94 */
95 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
96
97 #if defined (SCO)
98 #if ! defined (HAVE_ISNAN)
99 #define HAVE_ISNAN
100 static int
101 isnan (double x)
102 {
103 return (IsNANorINF (x) && NaN (x) && ! IsINF (x)) ? 1 : 0;
104 }
105 #endif
106 #if ! defined (HAVE_ISINF)
107 #define HAVE_ISINF
108 static int
109 isinf (double x)
110 {
111 return (IsNANorINF (x) && IsINF (x)) ? 1 : 0;
112 }
113
114 #endif
115 #endif
116
117
118 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
119 an explicit check. In some future gmp (don't know what version number),
120 mpz_cmp_d is supposed to do this itself. */
121 #if 1
122 #define xmpz_cmp_d(z, d) \
123 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
124 #else
125 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
126 #endif
127
128 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
129 isinf. It does have finite and isnan though, hence the use of those.
130 fpclass would be a possibility on that system too. */
131 static int
132 xisinf (double x)
133 {
134 #if defined (HAVE_ISINF)
135 return isinf (x);
136 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
137 return (! (finite (x) || isnan (x)));
138 #else
139 return 0;
140 #endif
141 }
142
143 static int
144 xisnan (double x)
145 {
146 #if defined (HAVE_ISNAN)
147 return isnan (x);
148 #else
149 return 0;
150 #endif
151 }
152
153 \f
154
155 static mpz_t z_negative_one;
156
157 \f
158
159 SCM_C_INLINE_KEYWORD SCM
160 scm_i_mkbig ()
161 {
162 /* Return a newly created bignum. */
163 SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
164 mpz_init (SCM_I_BIG_MPZ (z));
165 return z;
166 }
167
168 SCM_C_INLINE_KEYWORD SCM
169 scm_i_long2big (long x)
170 {
171 /* Return a newly created bignum initialized to X. */
172 SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
173 mpz_init_set_si (SCM_I_BIG_MPZ (z), x);
174 return z;
175 }
176
177 SCM_C_INLINE_KEYWORD SCM
178 scm_i_ulong2big (unsigned long x)
179 {
180 /* Return a newly created bignum initialized to X. */
181 SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
182 mpz_init_set_ui (SCM_I_BIG_MPZ (z), x);
183 return z;
184 }
185
186 SCM_C_INLINE_KEYWORD static SCM
187 scm_i_clonebig (SCM src_big, int same_sign_p)
188 {
189 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
190 SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
191 mpz_init_set (SCM_I_BIG_MPZ (z), SCM_I_BIG_MPZ (src_big));
192 if (!same_sign_p)
193 mpz_neg (SCM_I_BIG_MPZ (z), SCM_I_BIG_MPZ (z));
194 return z;
195 }
196
197 SCM_C_INLINE_KEYWORD int
198 scm_i_bigcmp (SCM x, SCM y)
199 {
200 /* Return neg if x < y, pos if x > y, and 0 if x == y */
201 /* presume we already know x and y are bignums */
202 int result = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
203 scm_remember_upto_here_2 (x, y);
204 return result;
205 }
206
207 SCM_C_INLINE_KEYWORD SCM
208 scm_i_dbl2big (double d)
209 {
210 /* results are only defined if d is an integer */
211 SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
212 mpz_init_set_d (SCM_I_BIG_MPZ (z), d);
213 return z;
214 }
215
216 /* Convert a integer in double representation to a SCM number. */
217
218 SCM_C_INLINE_KEYWORD SCM
219 scm_i_dbl2num (double u)
220 {
221 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
222 powers of 2, so there's no rounding when making "double" values
223 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
224 get rounded on a 64-bit machine, hence the "+1".
225
226 The use of floor() to force to an integer value ensures we get a
227 "numerically closest" value without depending on how a
228 double->long cast or how mpz_set_d will round. For reference,
229 double->long probably follows the hardware rounding mode,
230 mpz_set_d truncates towards zero. */
231
232 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
233 representable as a double? */
234
235 if (u < (double) (SCM_MOST_POSITIVE_FIXNUM+1)
236 && u >= (double) SCM_MOST_NEGATIVE_FIXNUM)
237 return SCM_I_MAKINUM ((long) u);
238 else
239 return scm_i_dbl2big (u);
240 }
241
242 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
243 with R5RS exact->inexact.
244
245 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
246 (ie. truncate towards zero), then adjust to get the closest double by
247 examining the next lower bit and adding 1 (to the absolute value) if
248 necessary.
249
250 Bignums exactly half way between representable doubles are rounded to the
251 next higher absolute value (ie. away from zero). This seems like an
252 adequate interpretation of R5RS "numerically closest", and it's easier
253 and faster than a full "nearest-even" style.
254
255 The bit test must be done on the absolute value of the mpz_t, which means
256 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
257 negatives as twos complement.
258
259 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
260 following the hardware rounding mode, but applied to the absolute value
261 of the mpz_t operand. This is not what we want so we put the high
262 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
263 mpz_get_d is supposed to always truncate towards zero.
264
265 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
266 is a slowdown. It'd be faster to pick out the relevant high bits with
267 mpz_getlimbn if we could be bothered coding that, and if the new
268 truncating gmp doesn't come out. */
269
270 double
271 scm_i_big2dbl (SCM b)
272 {
273 double result;
274 size_t bits;
275
276 bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2);
277
278 #if 1
279 {
280 /* Current GMP, eg. 4.1.3, force truncation towards zero */
281 mpz_t tmp;
282 if (bits > DBL_MANT_DIG)
283 {
284 size_t shift = bits - DBL_MANT_DIG;
285 mpz_init2 (tmp, DBL_MANT_DIG);
286 mpz_tdiv_q_2exp (tmp, SCM_I_BIG_MPZ (b), shift);
287 result = ldexp (mpz_get_d (tmp), shift);
288 mpz_clear (tmp);
289 }
290 else
291 {
292 result = mpz_get_d (SCM_I_BIG_MPZ (b));
293 }
294 }
295 #else
296 /* Future GMP */
297 result = mpz_get_d (SCM_I_BIG_MPZ (b));
298 #endif
299
300 if (bits > DBL_MANT_DIG)
301 {
302 unsigned long pos = bits - DBL_MANT_DIG - 1;
303 /* test bit number "pos" in absolute value */
304 if (mpz_getlimbn (SCM_I_BIG_MPZ (b), pos / GMP_NUMB_BITS)
305 & ((mp_limb_t) 1 << (pos % GMP_NUMB_BITS)))
306 {
307 result += ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b)), pos + 1);
308 }
309 }
310
311 scm_remember_upto_here_1 (b);
312 return result;
313 }
314
315 SCM_C_INLINE_KEYWORD SCM
316 scm_i_normbig (SCM b)
317 {
318 /* convert a big back to a fixnum if it'll fit */
319 /* presume b is a bignum */
320 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b)))
321 {
322 long val = mpz_get_si (SCM_I_BIG_MPZ (b));
323 if (SCM_FIXABLE (val))
324 b = SCM_I_MAKINUM (val);
325 }
326 return b;
327 }
328
329 static SCM_C_INLINE_KEYWORD SCM
330 scm_i_mpz2num (mpz_t b)
331 {
332 /* convert a mpz number to a SCM number. */
333 if (mpz_fits_slong_p (b))
334 {
335 long val = mpz_get_si (b);
336 if (SCM_FIXABLE (val))
337 return SCM_I_MAKINUM (val);
338 }
339
340 {
341 SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
342 mpz_init_set (SCM_I_BIG_MPZ (z), b);
343 return z;
344 }
345 }
346
347 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
348 static SCM scm_divide2real (SCM x, SCM y);
349
350 static SCM
351 scm_i_make_ratio (SCM numerator, SCM denominator)
352 #define FUNC_NAME "make-ratio"
353 {
354 /* First make sure the arguments are proper.
355 */
356 if (SCM_I_INUMP (denominator))
357 {
358 if (scm_is_eq (denominator, SCM_INUM0))
359 scm_num_overflow ("make-ratio");
360 if (scm_is_eq (denominator, SCM_I_MAKINUM(1)))
361 return numerator;
362 }
363 else
364 {
365 if (!(SCM_BIGP(denominator)))
366 SCM_WRONG_TYPE_ARG (2, denominator);
367 }
368 if (!SCM_I_INUMP (numerator) && !SCM_BIGP (numerator))
369 SCM_WRONG_TYPE_ARG (1, numerator);
370
371 /* Then flip signs so that the denominator is positive.
372 */
373 if (scm_is_true (scm_negative_p (denominator)))
374 {
375 numerator = scm_difference (numerator, SCM_UNDEFINED);
376 denominator = scm_difference (denominator, SCM_UNDEFINED);
377 }
378
379 /* Now consider for each of the four fixnum/bignum combinations
380 whether the rational number is really an integer.
381 */
382 if (SCM_I_INUMP (numerator))
383 {
384 long x = SCM_I_INUM (numerator);
385 if (scm_is_eq (numerator, SCM_INUM0))
386 return SCM_INUM0;
387 if (SCM_I_INUMP (denominator))
388 {
389 long y;
390 y = SCM_I_INUM (denominator);
391 if (x == y)
392 return SCM_I_MAKINUM(1);
393 if ((x % y) == 0)
394 return SCM_I_MAKINUM (x / y);
395 }
396 else
397 {
398 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
399 of that value for the denominator, as a bignum. Apart from
400 that case, abs(bignum) > abs(inum) so inum/bignum is not an
401 integer. */
402 if (x == SCM_MOST_NEGATIVE_FIXNUM
403 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator),
404 - SCM_MOST_NEGATIVE_FIXNUM) == 0)
405 return SCM_I_MAKINUM(-1);
406 }
407 }
408 else if (SCM_BIGP (numerator))
409 {
410 if (SCM_I_INUMP (denominator))
411 {
412 long yy = SCM_I_INUM (denominator);
413 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator), yy))
414 return scm_divide (numerator, denominator);
415 }
416 else
417 {
418 if (scm_is_eq (numerator, denominator))
419 return SCM_I_MAKINUM(1);
420 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator),
421 SCM_I_BIG_MPZ (denominator)))
422 return scm_divide(numerator, denominator);
423 }
424 }
425
426 /* No, it's a proper fraction.
427 */
428 return scm_double_cell (scm_tc16_fraction,
429 SCM_UNPACK (numerator),
430 SCM_UNPACK (denominator), 0);
431 }
432 #undef FUNC_NAME
433
434 static void scm_i_fraction_reduce (SCM z)
435 {
436 if (!(SCM_FRACTION_REDUCED (z)))
437 {
438 SCM divisor;
439 divisor = scm_gcd (SCM_FRACTION_NUMERATOR (z), SCM_FRACTION_DENOMINATOR (z));
440 if (!(scm_is_eq (divisor, SCM_I_MAKINUM(1))))
441 {
442 /* is this safe? */
443 SCM_FRACTION_SET_NUMERATOR (z, scm_divide (SCM_FRACTION_NUMERATOR (z), divisor));
444 SCM_FRACTION_SET_DENOMINATOR (z, scm_divide (SCM_FRACTION_DENOMINATOR (z), divisor));
445 }
446 SCM_FRACTION_REDUCED_SET (z);
447 }
448 }
449
450 double
451 scm_i_fraction2double (SCM z)
452 {
453 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z),
454 SCM_FRACTION_DENOMINATOR (z)));
455 }
456
457 SCM_DEFINE (scm_exact_p, "exact?", 1, 0, 0,
458 (SCM x),
459 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
460 "otherwise.")
461 #define FUNC_NAME s_scm_exact_p
462 {
463 if (SCM_I_INUMP (x))
464 return SCM_BOOL_T;
465 if (SCM_BIGP (x))
466 return SCM_BOOL_T;
467 if (SCM_FRACTIONP (x))
468 return SCM_BOOL_T;
469 if (SCM_NUMBERP (x))
470 return SCM_BOOL_F;
471 SCM_WRONG_TYPE_ARG (1, x);
472 }
473 #undef FUNC_NAME
474
475
476 SCM_DEFINE (scm_odd_p, "odd?", 1, 0, 0,
477 (SCM n),
478 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
479 "otherwise.")
480 #define FUNC_NAME s_scm_odd_p
481 {
482 if (SCM_I_INUMP (n))
483 {
484 long val = SCM_I_INUM (n);
485 return scm_from_bool ((val & 1L) != 0);
486 }
487 else if (SCM_BIGP (n))
488 {
489 int odd_p = mpz_odd_p (SCM_I_BIG_MPZ (n));
490 scm_remember_upto_here_1 (n);
491 return scm_from_bool (odd_p);
492 }
493 else if (scm_is_true (scm_inf_p (n)))
494 return SCM_BOOL_T;
495 else if (SCM_REALP (n))
496 {
497 double rem = fabs (fmod (SCM_REAL_VALUE(n), 2.0));
498 if (rem == 1.0)
499 return SCM_BOOL_T;
500 else if (rem == 0.0)
501 return SCM_BOOL_F;
502 else
503 SCM_WRONG_TYPE_ARG (1, n);
504 }
505 else
506 SCM_WRONG_TYPE_ARG (1, n);
507 }
508 #undef FUNC_NAME
509
510
511 SCM_DEFINE (scm_even_p, "even?", 1, 0, 0,
512 (SCM n),
513 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
514 "otherwise.")
515 #define FUNC_NAME s_scm_even_p
516 {
517 if (SCM_I_INUMP (n))
518 {
519 long val = SCM_I_INUM (n);
520 return scm_from_bool ((val & 1L) == 0);
521 }
522 else if (SCM_BIGP (n))
523 {
524 int even_p = mpz_even_p (SCM_I_BIG_MPZ (n));
525 scm_remember_upto_here_1 (n);
526 return scm_from_bool (even_p);
527 }
528 else if (scm_is_true (scm_inf_p (n)))
529 return SCM_BOOL_T;
530 else if (SCM_REALP (n))
531 {
532 double rem = fabs (fmod (SCM_REAL_VALUE(n), 2.0));
533 if (rem == 1.0)
534 return SCM_BOOL_F;
535 else if (rem == 0.0)
536 return SCM_BOOL_T;
537 else
538 SCM_WRONG_TYPE_ARG (1, n);
539 }
540 else
541 SCM_WRONG_TYPE_ARG (1, n);
542 }
543 #undef FUNC_NAME
544
545 SCM_DEFINE (scm_inf_p, "inf?", 1, 0, 0,
546 (SCM x),
547 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
548 "or @samp{-inf.0}, @code{#f} otherwise.")
549 #define FUNC_NAME s_scm_inf_p
550 {
551 if (SCM_REALP (x))
552 return scm_from_bool (xisinf (SCM_REAL_VALUE (x)));
553 else if (SCM_COMPLEXP (x))
554 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x))
555 || xisinf (SCM_COMPLEX_IMAG (x)));
556 else
557 return SCM_BOOL_F;
558 }
559 #undef FUNC_NAME
560
561 SCM_DEFINE (scm_nan_p, "nan?", 1, 0, 0,
562 (SCM n),
563 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
564 "otherwise.")
565 #define FUNC_NAME s_scm_nan_p
566 {
567 if (SCM_REALP (n))
568 return scm_from_bool (xisnan (SCM_REAL_VALUE (n)));
569 else if (SCM_COMPLEXP (n))
570 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n))
571 || xisnan (SCM_COMPLEX_IMAG (n)));
572 else
573 return SCM_BOOL_F;
574 }
575 #undef FUNC_NAME
576
577 /* Guile's idea of infinity. */
578 static double guile_Inf;
579
580 /* Guile's idea of not a number. */
581 static double guile_NaN;
582
583 static void
584 guile_ieee_init (void)
585 {
586 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
587
588 /* Some version of gcc on some old version of Linux used to crash when
589 trying to make Inf and NaN. */
590
591 #ifdef INFINITY
592 /* C99 INFINITY, when available.
593 FIXME: The standard allows for INFINITY to be something that overflows
594 at compile time. We ought to have a configure test to check for that
595 before trying to use it. (But in practice we believe this is not a
596 problem on any system guile is likely to target.) */
597 guile_Inf = INFINITY;
598 #elif HAVE_DINFINITY
599 /* OSF */
600 extern unsigned int DINFINITY[2];
601 guile_Inf = (*(X_CAST(double *, DINFINITY)));
602 #else
603 double tmp = 1e+10;
604 guile_Inf = tmp;
605 for (;;)
606 {
607 guile_Inf *= 1e+10;
608 if (guile_Inf == tmp)
609 break;
610 tmp = guile_Inf;
611 }
612 #endif
613
614 #endif
615
616 #if defined (HAVE_ISNAN)
617
618 #ifdef NAN
619 /* C99 NAN, when available */
620 guile_NaN = NAN;
621 #elif HAVE_DQNAN
622 /* OSF */
623 extern unsigned int DQNAN[2];
624 guile_NaN = (*(X_CAST(double *, DQNAN)));
625 #else
626 guile_NaN = guile_Inf / guile_Inf;
627 #endif
628
629 #endif
630 }
631
632 SCM_DEFINE (scm_inf, "inf", 0, 0, 0,
633 (void),
634 "Return Inf.")
635 #define FUNC_NAME s_scm_inf
636 {
637 static int initialized = 0;
638 if (! initialized)
639 {
640 guile_ieee_init ();
641 initialized = 1;
642 }
643 return scm_from_double (guile_Inf);
644 }
645 #undef FUNC_NAME
646
647 SCM_DEFINE (scm_nan, "nan", 0, 0, 0,
648 (void),
649 "Return NaN.")
650 #define FUNC_NAME s_scm_nan
651 {
652 static int initialized = 0;
653 if (!initialized)
654 {
655 guile_ieee_init ();
656 initialized = 1;
657 }
658 return scm_from_double (guile_NaN);
659 }
660 #undef FUNC_NAME
661
662
663 SCM_PRIMITIVE_GENERIC (scm_abs, "abs", 1, 0, 0,
664 (SCM x),
665 "Return the absolute value of @var{x}.")
666 #define FUNC_NAME
667 {
668 if (SCM_I_INUMP (x))
669 {
670 long int xx = SCM_I_INUM (x);
671 if (xx >= 0)
672 return x;
673 else if (SCM_POSFIXABLE (-xx))
674 return SCM_I_MAKINUM (-xx);
675 else
676 return scm_i_long2big (-xx);
677 }
678 else if (SCM_BIGP (x))
679 {
680 const int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
681 if (sgn < 0)
682 return scm_i_clonebig (x, 0);
683 else
684 return x;
685 }
686 else if (SCM_REALP (x))
687 {
688 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
689 double xx = SCM_REAL_VALUE (x);
690 if (xx < 0.0)
691 return scm_from_double (-xx);
692 else
693 return x;
694 }
695 else if (SCM_FRACTIONP (x))
696 {
697 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x))))
698 return x;
699 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED),
700 SCM_FRACTION_DENOMINATOR (x));
701 }
702 else
703 SCM_WTA_DISPATCH_1 (g_scm_abs, x, 1, s_scm_abs);
704 }
705 #undef FUNC_NAME
706
707
708 SCM_GPROC (s_quotient, "quotient", 2, 0, 0, scm_quotient, g_quotient);
709 /* "Return the quotient of the numbers @var{x} and @var{y}."
710 */
711 SCM
712 scm_quotient (SCM x, SCM y)
713 {
714 if (SCM_I_INUMP (x))
715 {
716 long xx = SCM_I_INUM (x);
717 if (SCM_I_INUMP (y))
718 {
719 long yy = SCM_I_INUM (y);
720 if (yy == 0)
721 scm_num_overflow (s_quotient);
722 else
723 {
724 long z = xx / yy;
725 if (SCM_FIXABLE (z))
726 return SCM_I_MAKINUM (z);
727 else
728 return scm_i_long2big (z);
729 }
730 }
731 else if (SCM_BIGP (y))
732 {
733 if ((SCM_I_INUM (x) == SCM_MOST_NEGATIVE_FIXNUM)
734 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y),
735 - SCM_MOST_NEGATIVE_FIXNUM) == 0))
736 {
737 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
738 scm_remember_upto_here_1 (y);
739 return SCM_I_MAKINUM (-1);
740 }
741 else
742 return SCM_I_MAKINUM (0);
743 }
744 else
745 SCM_WTA_DISPATCH_2 (g_quotient, x, y, SCM_ARG2, s_quotient);
746 }
747 else if (SCM_BIGP (x))
748 {
749 if (SCM_I_INUMP (y))
750 {
751 long yy = SCM_I_INUM (y);
752 if (yy == 0)
753 scm_num_overflow (s_quotient);
754 else if (yy == 1)
755 return x;
756 else
757 {
758 SCM result = scm_i_mkbig ();
759 if (yy < 0)
760 {
761 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result),
762 SCM_I_BIG_MPZ (x),
763 - yy);
764 mpz_neg (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result));
765 }
766 else
767 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), yy);
768 scm_remember_upto_here_1 (x);
769 return scm_i_normbig (result);
770 }
771 }
772 else if (SCM_BIGP (y))
773 {
774 SCM result = scm_i_mkbig ();
775 mpz_tdiv_q (SCM_I_BIG_MPZ (result),
776 SCM_I_BIG_MPZ (x),
777 SCM_I_BIG_MPZ (y));
778 scm_remember_upto_here_2 (x, y);
779 return scm_i_normbig (result);
780 }
781 else
782 SCM_WTA_DISPATCH_2 (g_quotient, x, y, SCM_ARG2, s_quotient);
783 }
784 else
785 SCM_WTA_DISPATCH_2 (g_quotient, x, y, SCM_ARG1, s_quotient);
786 }
787
788 SCM_GPROC (s_remainder, "remainder", 2, 0, 0, scm_remainder, g_remainder);
789 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
790 * "@lisp\n"
791 * "(remainder 13 4) @result{} 1\n"
792 * "(remainder -13 4) @result{} -1\n"
793 * "@end lisp"
794 */
795 SCM
796 scm_remainder (SCM x, SCM y)
797 {
798 if (SCM_I_INUMP (x))
799 {
800 if (SCM_I_INUMP (y))
801 {
802 long yy = SCM_I_INUM (y);
803 if (yy == 0)
804 scm_num_overflow (s_remainder);
805 else
806 {
807 long z = SCM_I_INUM (x) % yy;
808 return SCM_I_MAKINUM (z);
809 }
810 }
811 else if (SCM_BIGP (y))
812 {
813 if ((SCM_I_INUM (x) == SCM_MOST_NEGATIVE_FIXNUM)
814 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y),
815 - SCM_MOST_NEGATIVE_FIXNUM) == 0))
816 {
817 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
818 scm_remember_upto_here_1 (y);
819 return SCM_I_MAKINUM (0);
820 }
821 else
822 return x;
823 }
824 else
825 SCM_WTA_DISPATCH_2 (g_remainder, x, y, SCM_ARG2, s_remainder);
826 }
827 else if (SCM_BIGP (x))
828 {
829 if (SCM_I_INUMP (y))
830 {
831 long yy = SCM_I_INUM (y);
832 if (yy == 0)
833 scm_num_overflow (s_remainder);
834 else
835 {
836 SCM result = scm_i_mkbig ();
837 if (yy < 0)
838 yy = - yy;
839 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ(x), yy);
840 scm_remember_upto_here_1 (x);
841 return scm_i_normbig (result);
842 }
843 }
844 else if (SCM_BIGP (y))
845 {
846 SCM result = scm_i_mkbig ();
847 mpz_tdiv_r (SCM_I_BIG_MPZ (result),
848 SCM_I_BIG_MPZ (x),
849 SCM_I_BIG_MPZ (y));
850 scm_remember_upto_here_2 (x, y);
851 return scm_i_normbig (result);
852 }
853 else
854 SCM_WTA_DISPATCH_2 (g_remainder, x, y, SCM_ARG2, s_remainder);
855 }
856 else
857 SCM_WTA_DISPATCH_2 (g_remainder, x, y, SCM_ARG1, s_remainder);
858 }
859
860
861 SCM_GPROC (s_modulo, "modulo", 2, 0, 0, scm_modulo, g_modulo);
862 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
863 * "@lisp\n"
864 * "(modulo 13 4) @result{} 1\n"
865 * "(modulo -13 4) @result{} 3\n"
866 * "@end lisp"
867 */
868 SCM
869 scm_modulo (SCM x, SCM y)
870 {
871 if (SCM_I_INUMP (x))
872 {
873 long xx = SCM_I_INUM (x);
874 if (SCM_I_INUMP (y))
875 {
876 long yy = SCM_I_INUM (y);
877 if (yy == 0)
878 scm_num_overflow (s_modulo);
879 else
880 {
881 /* C99 specifies that "%" is the remainder corresponding to a
882 quotient rounded towards zero, and that's also traditional
883 for machine division, so z here should be well defined. */
884 long z = xx % yy;
885 long result;
886
887 if (yy < 0)
888 {
889 if (z > 0)
890 result = z + yy;
891 else
892 result = z;
893 }
894 else
895 {
896 if (z < 0)
897 result = z + yy;
898 else
899 result = z;
900 }
901 return SCM_I_MAKINUM (result);
902 }
903 }
904 else if (SCM_BIGP (y))
905 {
906 int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y));
907 {
908 mpz_t z_x;
909 SCM result;
910
911 if (sgn_y < 0)
912 {
913 SCM pos_y = scm_i_clonebig (y, 0);
914 /* do this after the last scm_op */
915 mpz_init_set_si (z_x, xx);
916 result = pos_y; /* re-use this bignum */
917 mpz_mod (SCM_I_BIG_MPZ (result),
918 z_x,
919 SCM_I_BIG_MPZ (pos_y));
920 scm_remember_upto_here_1 (pos_y);
921 }
922 else
923 {
924 result = scm_i_mkbig ();
925 /* do this after the last scm_op */
926 mpz_init_set_si (z_x, xx);
927 mpz_mod (SCM_I_BIG_MPZ (result),
928 z_x,
929 SCM_I_BIG_MPZ (y));
930 scm_remember_upto_here_1 (y);
931 }
932
933 if ((sgn_y < 0) && mpz_sgn (SCM_I_BIG_MPZ (result)) != 0)
934 mpz_add (SCM_I_BIG_MPZ (result),
935 SCM_I_BIG_MPZ (y),
936 SCM_I_BIG_MPZ (result));
937 scm_remember_upto_here_1 (y);
938 /* and do this before the next one */
939 mpz_clear (z_x);
940 return scm_i_normbig (result);
941 }
942 }
943 else
944 SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG2, s_modulo);
945 }
946 else if (SCM_BIGP (x))
947 {
948 if (SCM_I_INUMP (y))
949 {
950 long yy = SCM_I_INUM (y);
951 if (yy == 0)
952 scm_num_overflow (s_modulo);
953 else
954 {
955 SCM result = scm_i_mkbig ();
956 mpz_mod_ui (SCM_I_BIG_MPZ (result),
957 SCM_I_BIG_MPZ (x),
958 (yy < 0) ? - yy : yy);
959 scm_remember_upto_here_1 (x);
960 if ((yy < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0))
961 mpz_sub_ui (SCM_I_BIG_MPZ (result),
962 SCM_I_BIG_MPZ (result),
963 - yy);
964 return scm_i_normbig (result);
965 }
966 }
967 else if (SCM_BIGP (y))
968 {
969 {
970 SCM result = scm_i_mkbig ();
971 int y_sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
972 SCM pos_y = scm_i_clonebig (y, y_sgn >= 0);
973 mpz_mod (SCM_I_BIG_MPZ (result),
974 SCM_I_BIG_MPZ (x),
975 SCM_I_BIG_MPZ (pos_y));
976
977 scm_remember_upto_here_1 (x);
978 if ((y_sgn < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0))
979 mpz_add (SCM_I_BIG_MPZ (result),
980 SCM_I_BIG_MPZ (y),
981 SCM_I_BIG_MPZ (result));
982 scm_remember_upto_here_2 (y, pos_y);
983 return scm_i_normbig (result);
984 }
985 }
986 else
987 SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG2, s_modulo);
988 }
989 else
990 SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG1, s_modulo);
991 }
992
993 SCM_GPROC1 (s_gcd, "gcd", scm_tc7_asubr, scm_gcd, g_gcd);
994 /* "Return the greatest common divisor of all arguments.\n"
995 * "If called without arguments, 0 is returned."
996 */
997 SCM
998 scm_gcd (SCM x, SCM y)
999 {
1000 if (SCM_UNBNDP (y))
1001 return SCM_UNBNDP (x) ? SCM_INUM0 : x;
1002
1003 if (SCM_I_INUMP (x))
1004 {
1005 if (SCM_I_INUMP (y))
1006 {
1007 long xx = SCM_I_INUM (x);
1008 long yy = SCM_I_INUM (y);
1009 long u = xx < 0 ? -xx : xx;
1010 long v = yy < 0 ? -yy : yy;
1011 long result;
1012 if (xx == 0)
1013 result = v;
1014 else if (yy == 0)
1015 result = u;
1016 else
1017 {
1018 long k = 1;
1019 long t;
1020 /* Determine a common factor 2^k */
1021 while (!(1 & (u | v)))
1022 {
1023 k <<= 1;
1024 u >>= 1;
1025 v >>= 1;
1026 }
1027 /* Now, any factor 2^n can be eliminated */
1028 if (u & 1)
1029 t = -v;
1030 else
1031 {
1032 t = u;
1033 b3:
1034 t = SCM_SRS (t, 1);
1035 }
1036 if (!(1 & t))
1037 goto b3;
1038 if (t > 0)
1039 u = t;
1040 else
1041 v = -t;
1042 t = u - v;
1043 if (t != 0)
1044 goto b3;
1045 result = u * k;
1046 }
1047 return (SCM_POSFIXABLE (result)
1048 ? SCM_I_MAKINUM (result)
1049 : scm_i_long2big (result));
1050 }
1051 else if (SCM_BIGP (y))
1052 {
1053 SCM_SWAP (x, y);
1054 goto big_inum;
1055 }
1056 else
1057 SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG2, s_gcd);
1058 }
1059 else if (SCM_BIGP (x))
1060 {
1061 if (SCM_I_INUMP (y))
1062 {
1063 unsigned long result;
1064 long yy;
1065 big_inum:
1066 yy = SCM_I_INUM (y);
1067 if (yy == 0)
1068 return scm_abs (x);
1069 if (yy < 0)
1070 yy = -yy;
1071 result = mpz_gcd_ui (NULL, SCM_I_BIG_MPZ (x), yy);
1072 scm_remember_upto_here_1 (x);
1073 return (SCM_POSFIXABLE (result)
1074 ? SCM_I_MAKINUM (result)
1075 : scm_from_ulong (result));
1076 }
1077 else if (SCM_BIGP (y))
1078 {
1079 SCM result = scm_i_mkbig ();
1080 mpz_gcd (SCM_I_BIG_MPZ (result),
1081 SCM_I_BIG_MPZ (x),
1082 SCM_I_BIG_MPZ (y));
1083 scm_remember_upto_here_2 (x, y);
1084 return scm_i_normbig (result);
1085 }
1086 else
1087 SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG2, s_gcd);
1088 }
1089 else
1090 SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG1, s_gcd);
1091 }
1092
1093 SCM_GPROC1 (s_lcm, "lcm", scm_tc7_asubr, scm_lcm, g_lcm);
1094 /* "Return the least common multiple of the arguments.\n"
1095 * "If called without arguments, 1 is returned."
1096 */
1097 SCM
1098 scm_lcm (SCM n1, SCM n2)
1099 {
1100 if (SCM_UNBNDP (n2))
1101 {
1102 if (SCM_UNBNDP (n1))
1103 return SCM_I_MAKINUM (1L);
1104 n2 = SCM_I_MAKINUM (1L);
1105 }
1106
1107 SCM_GASSERT2 (SCM_I_INUMP (n1) || SCM_BIGP (n1),
1108 g_lcm, n1, n2, SCM_ARG1, s_lcm);
1109 SCM_GASSERT2 (SCM_I_INUMP (n2) || SCM_BIGP (n2),
1110 g_lcm, n1, n2, SCM_ARGn, s_lcm);
1111
1112 if (SCM_I_INUMP (n1))
1113 {
1114 if (SCM_I_INUMP (n2))
1115 {
1116 SCM d = scm_gcd (n1, n2);
1117 if (scm_is_eq (d, SCM_INUM0))
1118 return d;
1119 else
1120 return scm_abs (scm_product (n1, scm_quotient (n2, d)));
1121 }
1122 else
1123 {
1124 /* inum n1, big n2 */
1125 inumbig:
1126 {
1127 SCM result = scm_i_mkbig ();
1128 long nn1 = SCM_I_INUM (n1);
1129 if (nn1 == 0) return SCM_INUM0;
1130 if (nn1 < 0) nn1 = - nn1;
1131 mpz_lcm_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n2), nn1);
1132 scm_remember_upto_here_1 (n2);
1133 return result;
1134 }
1135 }
1136 }
1137 else
1138 {
1139 /* big n1 */
1140 if (SCM_I_INUMP (n2))
1141 {
1142 SCM_SWAP (n1, n2);
1143 goto inumbig;
1144 }
1145 else
1146 {
1147 SCM result = scm_i_mkbig ();
1148 mpz_lcm(SCM_I_BIG_MPZ (result),
1149 SCM_I_BIG_MPZ (n1),
1150 SCM_I_BIG_MPZ (n2));
1151 scm_remember_upto_here_2(n1, n2);
1152 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1153 return result;
1154 }
1155 }
1156 }
1157
1158 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1159
1160 Logand:
1161 X Y Result Method:
1162 (len)
1163 + + + x (map digit:logand X Y)
1164 + - + x (map digit:logand X (lognot (+ -1 Y)))
1165 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1166 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1167
1168 Logior:
1169 X Y Result Method:
1170
1171 + + + (map digit:logior X Y)
1172 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1173 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1174 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1175
1176 Logxor:
1177 X Y Result Method:
1178
1179 + + + (map digit:logxor X Y)
1180 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1181 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1182 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1183
1184 Logtest:
1185 X Y Result
1186
1187 + + (any digit:logand X Y)
1188 + - (any digit:logand X (lognot (+ -1 Y)))
1189 - + (any digit:logand (lognot (+ -1 X)) Y)
1190 - - #t
1191
1192 */
1193
1194 SCM_DEFINE1 (scm_logand, "logand", scm_tc7_asubr,
1195 (SCM n1, SCM n2),
1196 "Return the bitwise AND of the integer arguments.\n\n"
1197 "@lisp\n"
1198 "(logand) @result{} -1\n"
1199 "(logand 7) @result{} 7\n"
1200 "(logand #b111 #b011 #b001) @result{} 1\n"
1201 "@end lisp")
1202 #define FUNC_NAME s_scm_logand
1203 {
1204 long int nn1;
1205
1206 if (SCM_UNBNDP (n2))
1207 {
1208 if (SCM_UNBNDP (n1))
1209 return SCM_I_MAKINUM (-1);
1210 else if (!SCM_NUMBERP (n1))
1211 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1212 else if (SCM_NUMBERP (n1))
1213 return n1;
1214 else
1215 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1216 }
1217
1218 if (SCM_I_INUMP (n1))
1219 {
1220 nn1 = SCM_I_INUM (n1);
1221 if (SCM_I_INUMP (n2))
1222 {
1223 long nn2 = SCM_I_INUM (n2);
1224 return SCM_I_MAKINUM (nn1 & nn2);
1225 }
1226 else if SCM_BIGP (n2)
1227 {
1228 intbig:
1229 if (n1 == 0)
1230 return SCM_INUM0;
1231 {
1232 SCM result_z = scm_i_mkbig ();
1233 mpz_t nn1_z;
1234 mpz_init_set_si (nn1_z, nn1);
1235 mpz_and (SCM_I_BIG_MPZ (result_z), nn1_z, SCM_I_BIG_MPZ (n2));
1236 scm_remember_upto_here_1 (n2);
1237 mpz_clear (nn1_z);
1238 return scm_i_normbig (result_z);
1239 }
1240 }
1241 else
1242 SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
1243 }
1244 else if (SCM_BIGP (n1))
1245 {
1246 if (SCM_I_INUMP (n2))
1247 {
1248 SCM_SWAP (n1, n2);
1249 nn1 = SCM_I_INUM (n1);
1250 goto intbig;
1251 }
1252 else if (SCM_BIGP (n2))
1253 {
1254 SCM result_z = scm_i_mkbig ();
1255 mpz_and (SCM_I_BIG_MPZ (result_z),
1256 SCM_I_BIG_MPZ (n1),
1257 SCM_I_BIG_MPZ (n2));
1258 scm_remember_upto_here_2 (n1, n2);
1259 return scm_i_normbig (result_z);
1260 }
1261 else
1262 SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
1263 }
1264 else
1265 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1266 }
1267 #undef FUNC_NAME
1268
1269
1270 SCM_DEFINE1 (scm_logior, "logior", scm_tc7_asubr,
1271 (SCM n1, SCM n2),
1272 "Return the bitwise OR of the integer arguments.\n\n"
1273 "@lisp\n"
1274 "(logior) @result{} 0\n"
1275 "(logior 7) @result{} 7\n"
1276 "(logior #b000 #b001 #b011) @result{} 3\n"
1277 "@end lisp")
1278 #define FUNC_NAME s_scm_logior
1279 {
1280 long int nn1;
1281
1282 if (SCM_UNBNDP (n2))
1283 {
1284 if (SCM_UNBNDP (n1))
1285 return SCM_INUM0;
1286 else if (SCM_NUMBERP (n1))
1287 return n1;
1288 else
1289 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1290 }
1291
1292 if (SCM_I_INUMP (n1))
1293 {
1294 nn1 = SCM_I_INUM (n1);
1295 if (SCM_I_INUMP (n2))
1296 {
1297 long nn2 = SCM_I_INUM (n2);
1298 return SCM_I_MAKINUM (nn1 | nn2);
1299 }
1300 else if (SCM_BIGP (n2))
1301 {
1302 intbig:
1303 if (nn1 == 0)
1304 return n2;
1305 {
1306 SCM result_z = scm_i_mkbig ();
1307 mpz_t nn1_z;
1308 mpz_init_set_si (nn1_z, nn1);
1309 mpz_ior (SCM_I_BIG_MPZ (result_z), nn1_z, SCM_I_BIG_MPZ (n2));
1310 scm_remember_upto_here_1 (n2);
1311 mpz_clear (nn1_z);
1312 return result_z;
1313 }
1314 }
1315 else
1316 SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
1317 }
1318 else if (SCM_BIGP (n1))
1319 {
1320 if (SCM_I_INUMP (n2))
1321 {
1322 SCM_SWAP (n1, n2);
1323 nn1 = SCM_I_INUM (n1);
1324 goto intbig;
1325 }
1326 else if (SCM_BIGP (n2))
1327 {
1328 SCM result_z = scm_i_mkbig ();
1329 mpz_ior (SCM_I_BIG_MPZ (result_z),
1330 SCM_I_BIG_MPZ (n1),
1331 SCM_I_BIG_MPZ (n2));
1332 scm_remember_upto_here_2 (n1, n2);
1333 return result_z;
1334 }
1335 else
1336 SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
1337 }
1338 else
1339 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1340 }
1341 #undef FUNC_NAME
1342
1343
1344 SCM_DEFINE1 (scm_logxor, "logxor", scm_tc7_asubr,
1345 (SCM n1, SCM n2),
1346 "Return the bitwise XOR of the integer arguments. A bit is\n"
1347 "set in the result if it is set in an odd number of arguments.\n"
1348 "@lisp\n"
1349 "(logxor) @result{} 0\n"
1350 "(logxor 7) @result{} 7\n"
1351 "(logxor #b000 #b001 #b011) @result{} 2\n"
1352 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1353 "@end lisp")
1354 #define FUNC_NAME s_scm_logxor
1355 {
1356 long int nn1;
1357
1358 if (SCM_UNBNDP (n2))
1359 {
1360 if (SCM_UNBNDP (n1))
1361 return SCM_INUM0;
1362 else if (SCM_NUMBERP (n1))
1363 return n1;
1364 else
1365 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1366 }
1367
1368 if (SCM_I_INUMP (n1))
1369 {
1370 nn1 = SCM_I_INUM (n1);
1371 if (SCM_I_INUMP (n2))
1372 {
1373 long nn2 = SCM_I_INUM (n2);
1374 return SCM_I_MAKINUM (nn1 ^ nn2);
1375 }
1376 else if (SCM_BIGP (n2))
1377 {
1378 intbig:
1379 {
1380 SCM result_z = scm_i_mkbig ();
1381 mpz_t nn1_z;
1382 mpz_init_set_si (nn1_z, nn1);
1383 mpz_xor (SCM_I_BIG_MPZ (result_z), nn1_z, SCM_I_BIG_MPZ (n2));
1384 scm_remember_upto_here_1 (n2);
1385 mpz_clear (nn1_z);
1386 return scm_i_normbig (result_z);
1387 }
1388 }
1389 else
1390 SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
1391 }
1392 else if (SCM_BIGP (n1))
1393 {
1394 if (SCM_I_INUMP (n2))
1395 {
1396 SCM_SWAP (n1, n2);
1397 nn1 = SCM_I_INUM (n1);
1398 goto intbig;
1399 }
1400 else if (SCM_BIGP (n2))
1401 {
1402 SCM result_z = scm_i_mkbig ();
1403 mpz_xor (SCM_I_BIG_MPZ (result_z),
1404 SCM_I_BIG_MPZ (n1),
1405 SCM_I_BIG_MPZ (n2));
1406 scm_remember_upto_here_2 (n1, n2);
1407 return scm_i_normbig (result_z);
1408 }
1409 else
1410 SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
1411 }
1412 else
1413 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1414 }
1415 #undef FUNC_NAME
1416
1417
1418 SCM_DEFINE (scm_logtest, "logtest", 2, 0, 0,
1419 (SCM j, SCM k),
1420 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1421 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1422 "without actually calculating the @code{logand}, just testing\n"
1423 "for non-zero.\n"
1424 "\n"
1425 "@lisp\n"
1426 "(logtest #b0100 #b1011) @result{} #f\n"
1427 "(logtest #b0100 #b0111) @result{} #t\n"
1428 "@end lisp")
1429 #define FUNC_NAME s_scm_logtest
1430 {
1431 long int nj;
1432
1433 if (SCM_I_INUMP (j))
1434 {
1435 nj = SCM_I_INUM (j);
1436 if (SCM_I_INUMP (k))
1437 {
1438 long nk = SCM_I_INUM (k);
1439 return scm_from_bool (nj & nk);
1440 }
1441 else if (SCM_BIGP (k))
1442 {
1443 intbig:
1444 if (nj == 0)
1445 return SCM_BOOL_F;
1446 {
1447 SCM result;
1448 mpz_t nj_z;
1449 mpz_init_set_si (nj_z, nj);
1450 mpz_and (nj_z, nj_z, SCM_I_BIG_MPZ (k));
1451 scm_remember_upto_here_1 (k);
1452 result = scm_from_bool (mpz_sgn (nj_z) != 0);
1453 mpz_clear (nj_z);
1454 return result;
1455 }
1456 }
1457 else
1458 SCM_WRONG_TYPE_ARG (SCM_ARG2, k);
1459 }
1460 else if (SCM_BIGP (j))
1461 {
1462 if (SCM_I_INUMP (k))
1463 {
1464 SCM_SWAP (j, k);
1465 nj = SCM_I_INUM (j);
1466 goto intbig;
1467 }
1468 else if (SCM_BIGP (k))
1469 {
1470 SCM result;
1471 mpz_t result_z;
1472 mpz_init (result_z);
1473 mpz_and (result_z,
1474 SCM_I_BIG_MPZ (j),
1475 SCM_I_BIG_MPZ (k));
1476 scm_remember_upto_here_2 (j, k);
1477 result = scm_from_bool (mpz_sgn (result_z) != 0);
1478 mpz_clear (result_z);
1479 return result;
1480 }
1481 else
1482 SCM_WRONG_TYPE_ARG (SCM_ARG2, k);
1483 }
1484 else
1485 SCM_WRONG_TYPE_ARG (SCM_ARG1, j);
1486 }
1487 #undef FUNC_NAME
1488
1489
1490 SCM_DEFINE (scm_logbit_p, "logbit?", 2, 0, 0,
1491 (SCM index, SCM j),
1492 "Test whether bit number @var{index} in @var{j} is set.\n"
1493 "@var{index} starts from 0 for the least significant bit.\n"
1494 "\n"
1495 "@lisp\n"
1496 "(logbit? 0 #b1101) @result{} #t\n"
1497 "(logbit? 1 #b1101) @result{} #f\n"
1498 "(logbit? 2 #b1101) @result{} #t\n"
1499 "(logbit? 3 #b1101) @result{} #t\n"
1500 "(logbit? 4 #b1101) @result{} #f\n"
1501 "@end lisp")
1502 #define FUNC_NAME s_scm_logbit_p
1503 {
1504 unsigned long int iindex;
1505 iindex = scm_to_ulong (index);
1506
1507 if (SCM_I_INUMP (j))
1508 {
1509 /* bits above what's in an inum follow the sign bit */
1510 iindex = min (iindex, SCM_LONG_BIT - 1);
1511 return scm_from_bool ((1L << iindex) & SCM_I_INUM (j));
1512 }
1513 else if (SCM_BIGP (j))
1514 {
1515 int val = mpz_tstbit (SCM_I_BIG_MPZ (j), iindex);
1516 scm_remember_upto_here_1 (j);
1517 return scm_from_bool (val);
1518 }
1519 else
1520 SCM_WRONG_TYPE_ARG (SCM_ARG2, j);
1521 }
1522 #undef FUNC_NAME
1523
1524
1525 SCM_DEFINE (scm_lognot, "lognot", 1, 0, 0,
1526 (SCM n),
1527 "Return the integer which is the ones-complement of the integer\n"
1528 "argument.\n"
1529 "\n"
1530 "@lisp\n"
1531 "(number->string (lognot #b10000000) 2)\n"
1532 " @result{} \"-10000001\"\n"
1533 "(number->string (lognot #b0) 2)\n"
1534 " @result{} \"-1\"\n"
1535 "@end lisp")
1536 #define FUNC_NAME s_scm_lognot
1537 {
1538 if (SCM_I_INUMP (n)) {
1539 /* No overflow here, just need to toggle all the bits making up the inum.
1540 Enhancement: No need to strip the tag and add it back, could just xor
1541 a block of 1 bits, if that worked with the various debug versions of
1542 the SCM typedef. */
1543 return SCM_I_MAKINUM (~ SCM_I_INUM (n));
1544
1545 } else if (SCM_BIGP (n)) {
1546 SCM result = scm_i_mkbig ();
1547 mpz_com (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n));
1548 scm_remember_upto_here_1 (n);
1549 return result;
1550
1551 } else {
1552 SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
1553 }
1554 }
1555 #undef FUNC_NAME
1556
1557 /* returns 0 if IN is not an integer. OUT must already be
1558 initialized. */
1559 static int
1560 coerce_to_big (SCM in, mpz_t out)
1561 {
1562 if (SCM_BIGP (in))
1563 mpz_set (out, SCM_I_BIG_MPZ (in));
1564 else if (SCM_I_INUMP (in))
1565 mpz_set_si (out, SCM_I_INUM (in));
1566 else
1567 return 0;
1568
1569 return 1;
1570 }
1571
1572 SCM_DEFINE (scm_modulo_expt, "modulo-expt", 3, 0, 0,
1573 (SCM n, SCM k, SCM m),
1574 "Return @var{n} raised to the integer exponent\n"
1575 "@var{k}, modulo @var{m}.\n"
1576 "\n"
1577 "@lisp\n"
1578 "(modulo-expt 2 3 5)\n"
1579 " @result{} 3\n"
1580 "@end lisp")
1581 #define FUNC_NAME s_scm_modulo_expt
1582 {
1583 mpz_t n_tmp;
1584 mpz_t k_tmp;
1585 mpz_t m_tmp;
1586
1587 /* There are two classes of error we might encounter --
1588 1) Math errors, which we'll report by calling scm_num_overflow,
1589 and
1590 2) wrong-type errors, which of course we'll report by calling
1591 SCM_WRONG_TYPE_ARG.
1592 We don't report those errors immediately, however; instead we do
1593 some cleanup first. These variables tell us which error (if
1594 any) we should report after cleaning up.
1595 */
1596 int report_overflow = 0;
1597
1598 int position_of_wrong_type = 0;
1599 SCM value_of_wrong_type = SCM_INUM0;
1600
1601 SCM result = SCM_UNDEFINED;
1602
1603 mpz_init (n_tmp);
1604 mpz_init (k_tmp);
1605 mpz_init (m_tmp);
1606
1607 if (scm_is_eq (m, SCM_INUM0))
1608 {
1609 report_overflow = 1;
1610 goto cleanup;
1611 }
1612
1613 if (!coerce_to_big (n, n_tmp))
1614 {
1615 value_of_wrong_type = n;
1616 position_of_wrong_type = 1;
1617 goto cleanup;
1618 }
1619
1620 if (!coerce_to_big (k, k_tmp))
1621 {
1622 value_of_wrong_type = k;
1623 position_of_wrong_type = 2;
1624 goto cleanup;
1625 }
1626
1627 if (!coerce_to_big (m, m_tmp))
1628 {
1629 value_of_wrong_type = m;
1630 position_of_wrong_type = 3;
1631 goto cleanup;
1632 }
1633
1634 /* if the exponent K is negative, and we simply call mpz_powm, we
1635 will get a divide-by-zero exception when an inverse 1/n mod m
1636 doesn't exist (or is not unique). Since exceptions are hard to
1637 handle, we'll attempt the inversion "by hand" -- that way, we get
1638 a simple failure code, which is easy to handle. */
1639
1640 if (-1 == mpz_sgn (k_tmp))
1641 {
1642 if (!mpz_invert (n_tmp, n_tmp, m_tmp))
1643 {
1644 report_overflow = 1;
1645 goto cleanup;
1646 }
1647 mpz_neg (k_tmp, k_tmp);
1648 }
1649
1650 result = scm_i_mkbig ();
1651 mpz_powm (SCM_I_BIG_MPZ (result),
1652 n_tmp,
1653 k_tmp,
1654 m_tmp);
1655
1656 if (mpz_sgn (m_tmp) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result)) != 0)
1657 mpz_add (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result), m_tmp);
1658
1659 cleanup:
1660 mpz_clear (m_tmp);
1661 mpz_clear (k_tmp);
1662 mpz_clear (n_tmp);
1663
1664 if (report_overflow)
1665 scm_num_overflow (FUNC_NAME);
1666
1667 if (position_of_wrong_type)
1668 SCM_WRONG_TYPE_ARG (position_of_wrong_type,
1669 value_of_wrong_type);
1670
1671 return scm_i_normbig (result);
1672 }
1673 #undef FUNC_NAME
1674
1675 SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
1676 (SCM n, SCM k),
1677 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1678 "exact integer, @var{n} can be any number.\n"
1679 "\n"
1680 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1681 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1682 "includes @math{0^0} is 1.\n"
1683 "\n"
1684 "@lisp\n"
1685 "(integer-expt 2 5) @result{} 32\n"
1686 "(integer-expt -3 3) @result{} -27\n"
1687 "(integer-expt 5 -3) @result{} 1/125\n"
1688 "(integer-expt 0 0) @result{} 1\n"
1689 "@end lisp")
1690 #define FUNC_NAME s_scm_integer_expt
1691 {
1692 long i2 = 0;
1693 SCM z_i2 = SCM_BOOL_F;
1694 int i2_is_big = 0;
1695 SCM acc = SCM_I_MAKINUM (1L);
1696
1697 /* 0^0 == 1 according to R5RS */
1698 if (scm_is_eq (n, SCM_INUM0) || scm_is_eq (n, acc))
1699 return scm_is_false (scm_zero_p(k)) ? n : acc;
1700 else if (scm_is_eq (n, SCM_I_MAKINUM (-1L)))
1701 return scm_is_false (scm_even_p (k)) ? n : acc;
1702
1703 if (SCM_I_INUMP (k))
1704 i2 = SCM_I_INUM (k);
1705 else if (SCM_BIGP (k))
1706 {
1707 z_i2 = scm_i_clonebig (k, 1);
1708 scm_remember_upto_here_1 (k);
1709 i2_is_big = 1;
1710 }
1711 else
1712 SCM_WRONG_TYPE_ARG (2, k);
1713
1714 if (i2_is_big)
1715 {
1716 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2)) == -1)
1717 {
1718 mpz_neg (SCM_I_BIG_MPZ (z_i2), SCM_I_BIG_MPZ (z_i2));
1719 n = scm_divide (n, SCM_UNDEFINED);
1720 }
1721 while (1)
1722 {
1723 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2)) == 0)
1724 {
1725 return acc;
1726 }
1727 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2), 1) == 0)
1728 {
1729 return scm_product (acc, n);
1730 }
1731 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2), 0))
1732 acc = scm_product (acc, n);
1733 n = scm_product (n, n);
1734 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2), SCM_I_BIG_MPZ (z_i2), 1);
1735 }
1736 }
1737 else
1738 {
1739 if (i2 < 0)
1740 {
1741 i2 = -i2;
1742 n = scm_divide (n, SCM_UNDEFINED);
1743 }
1744 while (1)
1745 {
1746 if (0 == i2)
1747 return acc;
1748 if (1 == i2)
1749 return scm_product (acc, n);
1750 if (i2 & 1)
1751 acc = scm_product (acc, n);
1752 n = scm_product (n, n);
1753 i2 >>= 1;
1754 }
1755 }
1756 }
1757 #undef FUNC_NAME
1758
1759 SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
1760 (SCM n, SCM cnt),
1761 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1762 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1763 "\n"
1764 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1765 "@var{cnt} is negative it's a division, rounded towards negative\n"
1766 "infinity. (Note that this is not the same rounding as\n"
1767 "@code{quotient} does.)\n"
1768 "\n"
1769 "With @var{n} viewed as an infinite precision twos complement,\n"
1770 "@code{ash} means a left shift introducing zero bits, or a right\n"
1771 "shift dropping bits.\n"
1772 "\n"
1773 "@lisp\n"
1774 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1775 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1776 "\n"
1777 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1778 "(ash -23 -2) @result{} -6\n"
1779 "@end lisp")
1780 #define FUNC_NAME s_scm_ash
1781 {
1782 long bits_to_shift;
1783 bits_to_shift = scm_to_long (cnt);
1784
1785 if (SCM_I_INUMP (n))
1786 {
1787 long nn = SCM_I_INUM (n);
1788
1789 if (bits_to_shift > 0)
1790 {
1791 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1792 overflow a non-zero fixnum. For smaller shifts we check the
1793 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1794 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1795 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1796 bits_to_shift)". */
1797
1798 if (nn == 0)
1799 return n;
1800
1801 if (bits_to_shift < SCM_I_FIXNUM_BIT-1
1802 && ((unsigned long)
1803 (SCM_SRS (nn, (SCM_I_FIXNUM_BIT-1 - bits_to_shift)) + 1)
1804 <= 1))
1805 {
1806 return SCM_I_MAKINUM (nn << bits_to_shift);
1807 }
1808 else
1809 {
1810 SCM result = scm_i_long2big (nn);
1811 mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
1812 bits_to_shift);
1813 return result;
1814 }
1815 }
1816 else
1817 {
1818 bits_to_shift = -bits_to_shift;
1819 if (bits_to_shift >= SCM_LONG_BIT)
1820 return (nn >= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1821 else
1822 return SCM_I_MAKINUM (SCM_SRS (nn, bits_to_shift));
1823 }
1824
1825 }
1826 else if (SCM_BIGP (n))
1827 {
1828 SCM result;
1829
1830 if (bits_to_shift == 0)
1831 return n;
1832
1833 result = scm_i_mkbig ();
1834 if (bits_to_shift >= 0)
1835 {
1836 mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
1837 bits_to_shift);
1838 return result;
1839 }
1840 else
1841 {
1842 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1843 we have to allocate a bignum even if the result is going to be a
1844 fixnum. */
1845 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
1846 -bits_to_shift);
1847 return scm_i_normbig (result);
1848 }
1849
1850 }
1851 else
1852 {
1853 SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
1854 }
1855 }
1856 #undef FUNC_NAME
1857
1858
1859 SCM_DEFINE (scm_bit_extract, "bit-extract", 3, 0, 0,
1860 (SCM n, SCM start, SCM end),
1861 "Return the integer composed of the @var{start} (inclusive)\n"
1862 "through @var{end} (exclusive) bits of @var{n}. The\n"
1863 "@var{start}th bit becomes the 0-th bit in the result.\n"
1864 "\n"
1865 "@lisp\n"
1866 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1867 " @result{} \"1010\"\n"
1868 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1869 " @result{} \"10110\"\n"
1870 "@end lisp")
1871 #define FUNC_NAME s_scm_bit_extract
1872 {
1873 unsigned long int istart, iend, bits;
1874 istart = scm_to_ulong (start);
1875 iend = scm_to_ulong (end);
1876 SCM_ASSERT_RANGE (3, end, (iend >= istart));
1877
1878 /* how many bits to keep */
1879 bits = iend - istart;
1880
1881 if (SCM_I_INUMP (n))
1882 {
1883 long int in = SCM_I_INUM (n);
1884
1885 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1886 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1887 in = SCM_SRS (in, min (istart, SCM_I_FIXNUM_BIT-1));
1888
1889 if (in < 0 && bits >= SCM_I_FIXNUM_BIT)
1890 {
1891 /* Since we emulate two's complement encoded numbers, this
1892 * special case requires us to produce a result that has
1893 * more bits than can be stored in a fixnum.
1894 */
1895 SCM result = scm_i_long2big (in);
1896 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
1897 bits);
1898 return result;
1899 }
1900
1901 /* mask down to requisite bits */
1902 bits = min (bits, SCM_I_FIXNUM_BIT);
1903 return SCM_I_MAKINUM (in & ((1L << bits) - 1));
1904 }
1905 else if (SCM_BIGP (n))
1906 {
1907 SCM result;
1908 if (bits == 1)
1909 {
1910 result = SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n), istart));
1911 }
1912 else
1913 {
1914 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1915 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1916 such bits into a ulong. */
1917 result = scm_i_mkbig ();
1918 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result), SCM_I_BIG_MPZ(n), istart);
1919 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result), SCM_I_BIG_MPZ(result), bits);
1920 result = scm_i_normbig (result);
1921 }
1922 scm_remember_upto_here_1 (n);
1923 return result;
1924 }
1925 else
1926 SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
1927 }
1928 #undef FUNC_NAME
1929
1930
1931 static const char scm_logtab[] = {
1932 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1933 };
1934
1935 SCM_DEFINE (scm_logcount, "logcount", 1, 0, 0,
1936 (SCM n),
1937 "Return the number of bits in integer @var{n}. If integer is\n"
1938 "positive, the 1-bits in its binary representation are counted.\n"
1939 "If negative, the 0-bits in its two's-complement binary\n"
1940 "representation are counted. If 0, 0 is returned.\n"
1941 "\n"
1942 "@lisp\n"
1943 "(logcount #b10101010)\n"
1944 " @result{} 4\n"
1945 "(logcount 0)\n"
1946 " @result{} 0\n"
1947 "(logcount -2)\n"
1948 " @result{} 1\n"
1949 "@end lisp")
1950 #define FUNC_NAME s_scm_logcount
1951 {
1952 if (SCM_I_INUMP (n))
1953 {
1954 unsigned long int c = 0;
1955 long int nn = SCM_I_INUM (n);
1956 if (nn < 0)
1957 nn = -1 - nn;
1958 while (nn)
1959 {
1960 c += scm_logtab[15 & nn];
1961 nn >>= 4;
1962 }
1963 return SCM_I_MAKINUM (c);
1964 }
1965 else if (SCM_BIGP (n))
1966 {
1967 unsigned long count;
1968 if (mpz_sgn (SCM_I_BIG_MPZ (n)) >= 0)
1969 count = mpz_popcount (SCM_I_BIG_MPZ (n));
1970 else
1971 count = mpz_hamdist (SCM_I_BIG_MPZ (n), z_negative_one);
1972 scm_remember_upto_here_1 (n);
1973 return SCM_I_MAKINUM (count);
1974 }
1975 else
1976 SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
1977 }
1978 #undef FUNC_NAME
1979
1980
1981 static const char scm_ilentab[] = {
1982 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1983 };
1984
1985
1986 SCM_DEFINE (scm_integer_length, "integer-length", 1, 0, 0,
1987 (SCM n),
1988 "Return the number of bits necessary to represent @var{n}.\n"
1989 "\n"
1990 "@lisp\n"
1991 "(integer-length #b10101010)\n"
1992 " @result{} 8\n"
1993 "(integer-length 0)\n"
1994 " @result{} 0\n"
1995 "(integer-length #b1111)\n"
1996 " @result{} 4\n"
1997 "@end lisp")
1998 #define FUNC_NAME s_scm_integer_length
1999 {
2000 if (SCM_I_INUMP (n))
2001 {
2002 unsigned long int c = 0;
2003 unsigned int l = 4;
2004 long int nn = SCM_I_INUM (n);
2005 if (nn < 0)
2006 nn = -1 - nn;
2007 while (nn)
2008 {
2009 c += 4;
2010 l = scm_ilentab [15 & nn];
2011 nn >>= 4;
2012 }
2013 return SCM_I_MAKINUM (c - 4 + l);
2014 }
2015 else if (SCM_BIGP (n))
2016 {
2017 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2018 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2019 1 too big, so check for that and adjust. */
2020 size_t size = mpz_sizeinbase (SCM_I_BIG_MPZ (n), 2);
2021 if (mpz_sgn (SCM_I_BIG_MPZ (n)) < 0
2022 && mpz_scan0 (SCM_I_BIG_MPZ (n), /* no 0 bits above the lowest 1 */
2023 mpz_scan1 (SCM_I_BIG_MPZ (n), 0)) == ULONG_MAX)
2024 size--;
2025 scm_remember_upto_here_1 (n);
2026 return SCM_I_MAKINUM (size);
2027 }
2028 else
2029 SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
2030 }
2031 #undef FUNC_NAME
2032
2033 /*** NUMBERS -> STRINGS ***/
2034 #define SCM_MAX_DBL_PREC 60
2035 #define SCM_MAX_DBL_RADIX 36
2036
2037 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2038 static int scm_dblprec[SCM_MAX_DBL_RADIX - 1];
2039 static double fx_per_radix[SCM_MAX_DBL_RADIX - 1][SCM_MAX_DBL_PREC];
2040
2041 static
2042 void init_dblprec(int *prec, int radix) {
2043 /* determine floating point precision by adding successively
2044 smaller increments to 1.0 until it is considered == 1.0 */
2045 double f = ((double)1.0)/radix;
2046 double fsum = 1.0 + f;
2047
2048 *prec = 0;
2049 while (fsum != 1.0)
2050 {
2051 if (++(*prec) > SCM_MAX_DBL_PREC)
2052 fsum = 1.0;
2053 else
2054 {
2055 f /= radix;
2056 fsum = f + 1.0;
2057 }
2058 }
2059 (*prec) -= 1;
2060 }
2061
2062 static
2063 void init_fx_radix(double *fx_list, int radix)
2064 {
2065 /* initialize a per-radix list of tolerances. When added
2066 to a number < 1.0, we can determine if we should raund
2067 up and quit converting a number to a string. */
2068 int i;
2069 fx_list[0] = 0.0;
2070 fx_list[1] = 0.5;
2071 for( i = 2 ; i < SCM_MAX_DBL_PREC; ++i )
2072 fx_list[i] = (fx_list[i-1] / radix);
2073 }
2074
2075 /* use this array as a way to generate a single digit */
2076 static const char*number_chars="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2077
2078 static size_t
2079 idbl2str (double f, char *a, int radix)
2080 {
2081 int efmt, dpt, d, i, wp;
2082 double *fx;
2083 #ifdef DBL_MIN_10_EXP
2084 double f_cpy;
2085 int exp_cpy;
2086 #endif /* DBL_MIN_10_EXP */
2087 size_t ch = 0;
2088 int exp = 0;
2089
2090 if(radix < 2 ||
2091 radix > SCM_MAX_DBL_RADIX)
2092 {
2093 /* revert to existing behavior */
2094 radix = 10;
2095 }
2096
2097 wp = scm_dblprec[radix-2];
2098 fx = fx_per_radix[radix-2];
2099
2100 if (f == 0.0)
2101 {
2102 #ifdef HAVE_COPYSIGN
2103 double sgn = copysign (1.0, f);
2104
2105 if (sgn < 0.0)
2106 a[ch++] = '-';
2107 #endif
2108 goto zero; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2109 }
2110
2111 if (xisinf (f))
2112 {
2113 if (f < 0)
2114 strcpy (a, "-inf.0");
2115 else
2116 strcpy (a, "+inf.0");
2117 return ch+6;
2118 }
2119 else if (xisnan (f))
2120 {
2121 strcpy (a, "+nan.0");
2122 return ch+6;
2123 }
2124
2125 if (f < 0.0)
2126 {
2127 f = -f;
2128 a[ch++] = '-';
2129 }
2130
2131 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2132 make-uniform-vector, from causing infinite loops. */
2133 /* just do the checking...if it passes, we do the conversion for our
2134 radix again below */
2135 f_cpy = f;
2136 exp_cpy = exp;
2137
2138 while (f_cpy < 1.0)
2139 {
2140 f_cpy *= 10.0;
2141 if (exp_cpy-- < DBL_MIN_10_EXP)
2142 {
2143 a[ch++] = '#';
2144 a[ch++] = '.';
2145 a[ch++] = '#';
2146 return ch;
2147 }
2148 }
2149 while (f_cpy > 10.0)
2150 {
2151 f_cpy *= 0.10;
2152 if (exp_cpy++ > DBL_MAX_10_EXP)
2153 {
2154 a[ch++] = '#';
2155 a[ch++] = '.';
2156 a[ch++] = '#';
2157 return ch;
2158 }
2159 }
2160 #endif
2161
2162 while (f < 1.0)
2163 {
2164 f *= radix;
2165 exp--;
2166 }
2167 while (f > radix)
2168 {
2169 f /= radix;
2170 exp++;
2171 }
2172
2173 if (f + fx[wp] >= radix)
2174 {
2175 f = 1.0;
2176 exp++;
2177 }
2178 zero:
2179 #ifdef ENGNOT
2180 /* adding 9999 makes this equivalent to abs(x) % 3 */
2181 dpt = (exp + 9999) % 3;
2182 exp -= dpt++;
2183 efmt = 1;
2184 #else
2185 efmt = (exp < -3) || (exp > wp + 2);
2186 if (!efmt)
2187 {
2188 if (exp < 0)
2189 {
2190 a[ch++] = '0';
2191 a[ch++] = '.';
2192 dpt = exp;
2193 while (++dpt)
2194 a[ch++] = '0';
2195 }
2196 else
2197 dpt = exp + 1;
2198 }
2199 else
2200 dpt = 1;
2201 #endif
2202
2203 do
2204 {
2205 d = f;
2206 f -= d;
2207 a[ch++] = number_chars[d];
2208 if (f < fx[wp])
2209 break;
2210 if (f + fx[wp] >= 1.0)
2211 {
2212 a[ch - 1] = number_chars[d+1];
2213 break;
2214 }
2215 f *= radix;
2216 if (!(--dpt))
2217 a[ch++] = '.';
2218 }
2219 while (wp--);
2220
2221 if (dpt > 0)
2222 {
2223 #ifndef ENGNOT
2224 if ((dpt > 4) && (exp > 6))
2225 {
2226 d = (a[0] == '-' ? 2 : 1);
2227 for (i = ch++; i > d; i--)
2228 a[i] = a[i - 1];
2229 a[d] = '.';
2230 efmt = 1;
2231 }
2232 else
2233 #endif
2234 {
2235 while (--dpt)
2236 a[ch++] = '0';
2237 a[ch++] = '.';
2238 }
2239 }
2240 if (a[ch - 1] == '.')
2241 a[ch++] = '0'; /* trailing zero */
2242 if (efmt && exp)
2243 {
2244 a[ch++] = 'e';
2245 if (exp < 0)
2246 {
2247 exp = -exp;
2248 a[ch++] = '-';
2249 }
2250 for (i = radix; i <= exp; i *= radix);
2251 for (i /= radix; i; i /= radix)
2252 {
2253 a[ch++] = number_chars[exp / i];
2254 exp %= i;
2255 }
2256 }
2257 return ch;
2258 }
2259
2260
2261 static size_t
2262 icmplx2str (double real, double imag, char *str, int radix)
2263 {
2264 size_t i;
2265
2266 i = idbl2str (real, str, radix);
2267 if (imag != 0.0)
2268 {
2269 /* Don't output a '+' for negative numbers or for Inf and
2270 NaN. They will provide their own sign. */
2271 if (0 <= imag && !xisinf (imag) && !xisnan (imag))
2272 str[i++] = '+';
2273 i += idbl2str (imag, &str[i], radix);
2274 str[i++] = 'i';
2275 }
2276 return i;
2277 }
2278
2279 static size_t
2280 iflo2str (SCM flt, char *str, int radix)
2281 {
2282 size_t i;
2283 if (SCM_REALP (flt))
2284 i = idbl2str (SCM_REAL_VALUE (flt), str, radix);
2285 else
2286 i = icmplx2str (SCM_COMPLEX_REAL (flt), SCM_COMPLEX_IMAG (flt),
2287 str, radix);
2288 return i;
2289 }
2290
2291 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2292 characters in the result.
2293 rad is output base
2294 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2295 size_t
2296 scm_iint2str (scm_t_intmax num, int rad, char *p)
2297 {
2298 if (num < 0)
2299 {
2300 *p++ = '-';
2301 return scm_iuint2str (-num, rad, p) + 1;
2302 }
2303 else
2304 return scm_iuint2str (num, rad, p);
2305 }
2306
2307 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2308 characters in the result.
2309 rad is output base
2310 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2311 size_t
2312 scm_iuint2str (scm_t_uintmax num, int rad, char *p)
2313 {
2314 size_t j = 1;
2315 size_t i;
2316 scm_t_uintmax n = num;
2317
2318 for (n /= rad; n > 0; n /= rad)
2319 j++;
2320
2321 i = j;
2322 n = num;
2323 while (i--)
2324 {
2325 int d = n % rad;
2326
2327 n /= rad;
2328 p[i] = d + ((d < 10) ? '0' : 'a' - 10);
2329 }
2330 return j;
2331 }
2332
2333 SCM_DEFINE (scm_number_to_string, "number->string", 1, 1, 0,
2334 (SCM n, SCM radix),
2335 "Return a string holding the external representation of the\n"
2336 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2337 "inexact, a radix of 10 will be used.")
2338 #define FUNC_NAME s_scm_number_to_string
2339 {
2340 int base;
2341
2342 if (SCM_UNBNDP (radix))
2343 base = 10;
2344 else
2345 base = scm_to_signed_integer (radix, 2, 36);
2346
2347 if (SCM_I_INUMP (n))
2348 {
2349 char num_buf [SCM_INTBUFLEN];
2350 size_t length = scm_iint2str (SCM_I_INUM (n), base, num_buf);
2351 return scm_from_locale_stringn (num_buf, length);
2352 }
2353 else if (SCM_BIGP (n))
2354 {
2355 char *str = mpz_get_str (NULL, base, SCM_I_BIG_MPZ (n));
2356 scm_remember_upto_here_1 (n);
2357 return scm_take_locale_string (str);
2358 }
2359 else if (SCM_FRACTIONP (n))
2360 {
2361 scm_i_fraction_reduce (n);
2362 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n), radix),
2363 scm_from_locale_string ("/"),
2364 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n), radix)));
2365 }
2366 else if (SCM_INEXACTP (n))
2367 {
2368 char num_buf [FLOBUFLEN];
2369 return scm_from_locale_stringn (num_buf, iflo2str (n, num_buf, base));
2370 }
2371 else
2372 SCM_WRONG_TYPE_ARG (1, n);
2373 }
2374 #undef FUNC_NAME
2375
2376
2377 /* These print routines used to be stubbed here so that scm_repl.c
2378 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2379
2380 int
2381 scm_print_real (SCM sexp, SCM port, scm_print_state *pstate SCM_UNUSED)
2382 {
2383 char num_buf[FLOBUFLEN];
2384 scm_lfwrite (num_buf, iflo2str (sexp, num_buf, 10), port);
2385 return !0;
2386 }
2387
2388 void
2389 scm_i_print_double (double val, SCM port)
2390 {
2391 char num_buf[FLOBUFLEN];
2392 scm_lfwrite (num_buf, idbl2str (val, num_buf, 10), port);
2393 }
2394
2395 int
2396 scm_print_complex (SCM sexp, SCM port, scm_print_state *pstate SCM_UNUSED)
2397
2398 {
2399 char num_buf[FLOBUFLEN];
2400 scm_lfwrite (num_buf, iflo2str (sexp, num_buf, 10), port);
2401 return !0;
2402 }
2403
2404 void
2405 scm_i_print_complex (double real, double imag, SCM port)
2406 {
2407 char num_buf[FLOBUFLEN];
2408 scm_lfwrite (num_buf, icmplx2str (real, imag, num_buf, 10), port);
2409 }
2410
2411 int
2412 scm_i_print_fraction (SCM sexp, SCM port, scm_print_state *pstate SCM_UNUSED)
2413 {
2414 SCM str;
2415 scm_i_fraction_reduce (sexp);
2416 str = scm_number_to_string (sexp, SCM_UNDEFINED);
2417 scm_lfwrite (scm_i_string_chars (str), scm_i_string_length (str), port);
2418 scm_remember_upto_here_1 (str);
2419 return !0;
2420 }
2421
2422 int
2423 scm_bigprint (SCM exp, SCM port, scm_print_state *pstate SCM_UNUSED)
2424 {
2425 char *str = mpz_get_str (NULL, 10, SCM_I_BIG_MPZ (exp));
2426 scm_remember_upto_here_1 (exp);
2427 scm_lfwrite (str, (size_t) strlen (str), port);
2428 free (str);
2429 return !0;
2430 }
2431 /*** END nums->strs ***/
2432
2433
2434 /*** STRINGS -> NUMBERS ***/
2435
2436 /* The following functions implement the conversion from strings to numbers.
2437 * The implementation somehow follows the grammar for numbers as it is given
2438 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2439 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2440 * points should be noted about the implementation:
2441 * * Each function keeps a local index variable 'idx' that points at the
2442 * current position within the parsed string. The global index is only
2443 * updated if the function could parse the corresponding syntactic unit
2444 * successfully.
2445 * * Similarly, the functions keep track of indicators of inexactness ('#',
2446 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2447 * global exactness information is only updated after each part has been
2448 * successfully parsed.
2449 * * Sequences of digits are parsed into temporary variables holding fixnums.
2450 * Only if these fixnums would overflow, the result variables are updated
2451 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2452 * the temporary variables holding the fixnums are cleared, and the process
2453 * starts over again. If for example fixnums were able to store five decimal
2454 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2455 * and the result was computed as 12345 * 100000 + 67890. In other words,
2456 * only every five digits two bignum operations were performed.
2457 */
2458
2459 enum t_exactness {NO_EXACTNESS, INEXACT, EXACT};
2460
2461 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2462
2463 /* In non ASCII-style encodings the following macro might not work. */
2464 #define XDIGIT2UINT(d) \
2465 (isdigit ((int) (unsigned char) d) \
2466 ? (d) - '0' \
2467 : tolower ((int) (unsigned char) d) - 'a' + 10)
2468
2469 static SCM
2470 mem2uinteger (const char* mem, size_t len, unsigned int *p_idx,
2471 unsigned int radix, enum t_exactness *p_exactness)
2472 {
2473 unsigned int idx = *p_idx;
2474 unsigned int hash_seen = 0;
2475 scm_t_bits shift = 1;
2476 scm_t_bits add = 0;
2477 unsigned int digit_value;
2478 SCM result;
2479 char c;
2480
2481 if (idx == len)
2482 return SCM_BOOL_F;
2483
2484 c = mem[idx];
2485 if (!isxdigit ((int) (unsigned char) c))
2486 return SCM_BOOL_F;
2487 digit_value = XDIGIT2UINT (c);
2488 if (digit_value >= radix)
2489 return SCM_BOOL_F;
2490
2491 idx++;
2492 result = SCM_I_MAKINUM (digit_value);
2493 while (idx != len)
2494 {
2495 char c = mem[idx];
2496 if (isxdigit ((int) (unsigned char) c))
2497 {
2498 if (hash_seen)
2499 break;
2500 digit_value = XDIGIT2UINT (c);
2501 if (digit_value >= radix)
2502 break;
2503 }
2504 else if (c == '#')
2505 {
2506 hash_seen = 1;
2507 digit_value = 0;
2508 }
2509 else
2510 break;
2511
2512 idx++;
2513 if (SCM_MOST_POSITIVE_FIXNUM / radix < shift)
2514 {
2515 result = scm_product (result, SCM_I_MAKINUM (shift));
2516 if (add > 0)
2517 result = scm_sum (result, SCM_I_MAKINUM (add));
2518
2519 shift = radix;
2520 add = digit_value;
2521 }
2522 else
2523 {
2524 shift = shift * radix;
2525 add = add * radix + digit_value;
2526 }
2527 };
2528
2529 if (shift > 1)
2530 result = scm_product (result, SCM_I_MAKINUM (shift));
2531 if (add > 0)
2532 result = scm_sum (result, SCM_I_MAKINUM (add));
2533
2534 *p_idx = idx;
2535 if (hash_seen)
2536 *p_exactness = INEXACT;
2537
2538 return result;
2539 }
2540
2541
2542 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2543 * covers the parts of the rules that start at a potential point. The value
2544 * of the digits up to the point have been parsed by the caller and are given
2545 * in variable result. The content of *p_exactness indicates, whether a hash
2546 * has already been seen in the digits before the point.
2547 */
2548
2549 /* In non ASCII-style encodings the following macro might not work. */
2550 #define DIGIT2UINT(d) ((d) - '0')
2551
2552 static SCM
2553 mem2decimal_from_point (SCM result, const char* mem, size_t len,
2554 unsigned int *p_idx, enum t_exactness *p_exactness)
2555 {
2556 unsigned int idx = *p_idx;
2557 enum t_exactness x = *p_exactness;
2558
2559 if (idx == len)
2560 return result;
2561
2562 if (mem[idx] == '.')
2563 {
2564 scm_t_bits shift = 1;
2565 scm_t_bits add = 0;
2566 unsigned int digit_value;
2567 SCM big_shift = SCM_I_MAKINUM (1);
2568
2569 idx++;
2570 while (idx != len)
2571 {
2572 char c = mem[idx];
2573 if (isdigit ((int) (unsigned char) c))
2574 {
2575 if (x == INEXACT)
2576 return SCM_BOOL_F;
2577 else
2578 digit_value = DIGIT2UINT (c);
2579 }
2580 else if (c == '#')
2581 {
2582 x = INEXACT;
2583 digit_value = 0;
2584 }
2585 else
2586 break;
2587
2588 idx++;
2589 if (SCM_MOST_POSITIVE_FIXNUM / 10 < shift)
2590 {
2591 big_shift = scm_product (big_shift, SCM_I_MAKINUM (shift));
2592 result = scm_product (result, SCM_I_MAKINUM (shift));
2593 if (add > 0)
2594 result = scm_sum (result, SCM_I_MAKINUM (add));
2595
2596 shift = 10;
2597 add = digit_value;
2598 }
2599 else
2600 {
2601 shift = shift * 10;
2602 add = add * 10 + digit_value;
2603 }
2604 };
2605
2606 if (add > 0)
2607 {
2608 big_shift = scm_product (big_shift, SCM_I_MAKINUM (shift));
2609 result = scm_product (result, SCM_I_MAKINUM (shift));
2610 result = scm_sum (result, SCM_I_MAKINUM (add));
2611 }
2612
2613 result = scm_divide (result, big_shift);
2614
2615 /* We've seen a decimal point, thus the value is implicitly inexact. */
2616 x = INEXACT;
2617 }
2618
2619 if (idx != len)
2620 {
2621 int sign = 1;
2622 unsigned int start;
2623 char c;
2624 int exponent;
2625 SCM e;
2626
2627 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2628
2629 switch (mem[idx])
2630 {
2631 case 'd': case 'D':
2632 case 'e': case 'E':
2633 case 'f': case 'F':
2634 case 'l': case 'L':
2635 case 's': case 'S':
2636 idx++;
2637 start = idx;
2638 c = mem[idx];
2639 if (c == '-')
2640 {
2641 idx++;
2642 sign = -1;
2643 c = mem[idx];
2644 }
2645 else if (c == '+')
2646 {
2647 idx++;
2648 sign = 1;
2649 c = mem[idx];
2650 }
2651 else
2652 sign = 1;
2653
2654 if (!isdigit ((int) (unsigned char) c))
2655 return SCM_BOOL_F;
2656
2657 idx++;
2658 exponent = DIGIT2UINT (c);
2659 while (idx != len)
2660 {
2661 char c = mem[idx];
2662 if (isdigit ((int) (unsigned char) c))
2663 {
2664 idx++;
2665 if (exponent <= SCM_MAXEXP)
2666 exponent = exponent * 10 + DIGIT2UINT (c);
2667 }
2668 else
2669 break;
2670 }
2671
2672 if (exponent > SCM_MAXEXP)
2673 {
2674 size_t exp_len = idx - start;
2675 SCM exp_string = scm_from_locale_stringn (&mem[start], exp_len);
2676 SCM exp_num = scm_string_to_number (exp_string, SCM_UNDEFINED);
2677 scm_out_of_range ("string->number", exp_num);
2678 }
2679
2680 e = scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent));
2681 if (sign == 1)
2682 result = scm_product (result, e);
2683 else
2684 result = scm_divide2real (result, e);
2685
2686 /* We've seen an exponent, thus the value is implicitly inexact. */
2687 x = INEXACT;
2688
2689 break;
2690
2691 default:
2692 break;
2693 }
2694 }
2695
2696 *p_idx = idx;
2697 if (x == INEXACT)
2698 *p_exactness = x;
2699
2700 return result;
2701 }
2702
2703
2704 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2705
2706 static SCM
2707 mem2ureal (const char* mem, size_t len, unsigned int *p_idx,
2708 unsigned int radix, enum t_exactness *p_exactness)
2709 {
2710 unsigned int idx = *p_idx;
2711 SCM result;
2712
2713 if (idx == len)
2714 return SCM_BOOL_F;
2715
2716 if (idx+5 <= len && !strncmp (mem+idx, "inf.0", 5))
2717 {
2718 *p_idx = idx+5;
2719 return scm_inf ();
2720 }
2721
2722 if (idx+4 < len && !strncmp (mem+idx, "nan.", 4))
2723 {
2724 enum t_exactness x = EXACT;
2725
2726 /* Cobble up the fractional part. We might want to set the
2727 NaN's mantissa from it. */
2728 idx += 4;
2729 mem2uinteger (mem, len, &idx, 10, &x);
2730 *p_idx = idx;
2731 return scm_nan ();
2732 }
2733
2734 if (mem[idx] == '.')
2735 {
2736 if (radix != 10)
2737 return SCM_BOOL_F;
2738 else if (idx + 1 == len)
2739 return SCM_BOOL_F;
2740 else if (!isdigit ((int) (unsigned char) mem[idx + 1]))
2741 return SCM_BOOL_F;
2742 else
2743 result = mem2decimal_from_point (SCM_I_MAKINUM (0), mem, len,
2744 p_idx, p_exactness);
2745 }
2746 else
2747 {
2748 enum t_exactness x = EXACT;
2749 SCM uinteger;
2750
2751 uinteger = mem2uinteger (mem, len, &idx, radix, &x);
2752 if (scm_is_false (uinteger))
2753 return SCM_BOOL_F;
2754
2755 if (idx == len)
2756 result = uinteger;
2757 else if (mem[idx] == '/')
2758 {
2759 SCM divisor;
2760
2761 idx++;
2762
2763 divisor = mem2uinteger (mem, len, &idx, radix, &x);
2764 if (scm_is_false (divisor))
2765 return SCM_BOOL_F;
2766
2767 /* both are int/big here, I assume */
2768 result = scm_i_make_ratio (uinteger, divisor);
2769 }
2770 else if (radix == 10)
2771 {
2772 result = mem2decimal_from_point (uinteger, mem, len, &idx, &x);
2773 if (scm_is_false (result))
2774 return SCM_BOOL_F;
2775 }
2776 else
2777 result = uinteger;
2778
2779 *p_idx = idx;
2780 if (x == INEXACT)
2781 *p_exactness = x;
2782 }
2783
2784 /* When returning an inexact zero, make sure it is represented as a
2785 floating point value so that we can change its sign.
2786 */
2787 if (scm_is_eq (result, SCM_I_MAKINUM(0)) && *p_exactness == INEXACT)
2788 result = scm_from_double (0.0);
2789
2790 return result;
2791 }
2792
2793
2794 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2795
2796 static SCM
2797 mem2complex (const char* mem, size_t len, unsigned int idx,
2798 unsigned int radix, enum t_exactness *p_exactness)
2799 {
2800 char c;
2801 int sign = 0;
2802 SCM ureal;
2803
2804 if (idx == len)
2805 return SCM_BOOL_F;
2806
2807 c = mem[idx];
2808 if (c == '+')
2809 {
2810 idx++;
2811 sign = 1;
2812 }
2813 else if (c == '-')
2814 {
2815 idx++;
2816 sign = -1;
2817 }
2818
2819 if (idx == len)
2820 return SCM_BOOL_F;
2821
2822 ureal = mem2ureal (mem, len, &idx, radix, p_exactness);
2823 if (scm_is_false (ureal))
2824 {
2825 /* input must be either +i or -i */
2826
2827 if (sign == 0)
2828 return SCM_BOOL_F;
2829
2830 if (mem[idx] == 'i' || mem[idx] == 'I')
2831 {
2832 idx++;
2833 if (idx != len)
2834 return SCM_BOOL_F;
2835
2836 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign));
2837 }
2838 else
2839 return SCM_BOOL_F;
2840 }
2841 else
2842 {
2843 if (sign == -1 && scm_is_false (scm_nan_p (ureal)))
2844 ureal = scm_difference (ureal, SCM_UNDEFINED);
2845
2846 if (idx == len)
2847 return ureal;
2848
2849 c = mem[idx];
2850 switch (c)
2851 {
2852 case 'i': case 'I':
2853 /* either +<ureal>i or -<ureal>i */
2854
2855 idx++;
2856 if (sign == 0)
2857 return SCM_BOOL_F;
2858 if (idx != len)
2859 return SCM_BOOL_F;
2860 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal);
2861
2862 case '@':
2863 /* polar input: <real>@<real>. */
2864
2865 idx++;
2866 if (idx == len)
2867 return SCM_BOOL_F;
2868 else
2869 {
2870 int sign;
2871 SCM angle;
2872 SCM result;
2873
2874 c = mem[idx];
2875 if (c == '+')
2876 {
2877 idx++;
2878 sign = 1;
2879 }
2880 else if (c == '-')
2881 {
2882 idx++;
2883 sign = -1;
2884 }
2885 else
2886 sign = 1;
2887
2888 angle = mem2ureal (mem, len, &idx, radix, p_exactness);
2889 if (scm_is_false (angle))
2890 return SCM_BOOL_F;
2891 if (idx != len)
2892 return SCM_BOOL_F;
2893
2894 if (sign == -1 && scm_is_false (scm_nan_p (ureal)))
2895 angle = scm_difference (angle, SCM_UNDEFINED);
2896
2897 result = scm_make_polar (ureal, angle);
2898 return result;
2899 }
2900 case '+':
2901 case '-':
2902 /* expecting input matching <real>[+-]<ureal>?i */
2903
2904 idx++;
2905 if (idx == len)
2906 return SCM_BOOL_F;
2907 else
2908 {
2909 int sign = (c == '+') ? 1 : -1;
2910 SCM imag = mem2ureal (mem, len, &idx, radix, p_exactness);
2911
2912 if (scm_is_false (imag))
2913 imag = SCM_I_MAKINUM (sign);
2914 else if (sign == -1 && scm_is_false (scm_nan_p (ureal)))
2915 imag = scm_difference (imag, SCM_UNDEFINED);
2916
2917 if (idx == len)
2918 return SCM_BOOL_F;
2919 if (mem[idx] != 'i' && mem[idx] != 'I')
2920 return SCM_BOOL_F;
2921
2922 idx++;
2923 if (idx != len)
2924 return SCM_BOOL_F;
2925
2926 return scm_make_rectangular (ureal, imag);
2927 }
2928 default:
2929 return SCM_BOOL_F;
2930 }
2931 }
2932 }
2933
2934
2935 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2936
2937 enum t_radix {NO_RADIX=0, DUAL=2, OCT=8, DEC=10, HEX=16};
2938
2939 SCM
2940 scm_i_mem2number (const char* mem, size_t len, unsigned int default_radix)
2941 {
2942 unsigned int idx = 0;
2943 unsigned int radix = NO_RADIX;
2944 enum t_exactness forced_x = NO_EXACTNESS;
2945 enum t_exactness implicit_x = EXACT;
2946 SCM result;
2947
2948 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2949 while (idx + 2 < len && mem[idx] == '#')
2950 {
2951 switch (mem[idx + 1])
2952 {
2953 case 'b': case 'B':
2954 if (radix != NO_RADIX)
2955 return SCM_BOOL_F;
2956 radix = DUAL;
2957 break;
2958 case 'd': case 'D':
2959 if (radix != NO_RADIX)
2960 return SCM_BOOL_F;
2961 radix = DEC;
2962 break;
2963 case 'i': case 'I':
2964 if (forced_x != NO_EXACTNESS)
2965 return SCM_BOOL_F;
2966 forced_x = INEXACT;
2967 break;
2968 case 'e': case 'E':
2969 if (forced_x != NO_EXACTNESS)
2970 return SCM_BOOL_F;
2971 forced_x = EXACT;
2972 break;
2973 case 'o': case 'O':
2974 if (radix != NO_RADIX)
2975 return SCM_BOOL_F;
2976 radix = OCT;
2977 break;
2978 case 'x': case 'X':
2979 if (radix != NO_RADIX)
2980 return SCM_BOOL_F;
2981 radix = HEX;
2982 break;
2983 default:
2984 return SCM_BOOL_F;
2985 }
2986 idx += 2;
2987 }
2988
2989 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2990 if (radix == NO_RADIX)
2991 result = mem2complex (mem, len, idx, default_radix, &implicit_x);
2992 else
2993 result = mem2complex (mem, len, idx, (unsigned int) radix, &implicit_x);
2994
2995 if (scm_is_false (result))
2996 return SCM_BOOL_F;
2997
2998 switch (forced_x)
2999 {
3000 case EXACT:
3001 if (SCM_INEXACTP (result))
3002 return scm_inexact_to_exact (result);
3003 else
3004 return result;
3005 case INEXACT:
3006 if (SCM_INEXACTP (result))
3007 return result;
3008 else
3009 return scm_exact_to_inexact (result);
3010 case NO_EXACTNESS:
3011 default:
3012 if (implicit_x == INEXACT)
3013 {
3014 if (SCM_INEXACTP (result))
3015 return result;
3016 else
3017 return scm_exact_to_inexact (result);
3018 }
3019 else
3020 return result;
3021 }
3022 }
3023
3024
3025 SCM_DEFINE (scm_string_to_number, "string->number", 1, 1, 0,
3026 (SCM string, SCM radix),
3027 "Return a number of the maximally precise representation\n"
3028 "expressed by the given @var{string}. @var{radix} must be an\n"
3029 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3030 "is a default radix that may be overridden by an explicit radix\n"
3031 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3032 "supplied, then the default radix is 10. If string is not a\n"
3033 "syntactically valid notation for a number, then\n"
3034 "@code{string->number} returns @code{#f}.")
3035 #define FUNC_NAME s_scm_string_to_number
3036 {
3037 SCM answer;
3038 unsigned int base;
3039 SCM_VALIDATE_STRING (1, string);
3040
3041 if (SCM_UNBNDP (radix))
3042 base = 10;
3043 else
3044 base = scm_to_unsigned_integer (radix, 2, INT_MAX);
3045
3046 answer = scm_i_mem2number (scm_i_string_chars (string),
3047 scm_i_string_length (string),
3048 base);
3049 scm_remember_upto_here_1 (string);
3050 return answer;
3051 }
3052 #undef FUNC_NAME
3053
3054
3055 /*** END strs->nums ***/
3056
3057
3058 SCM
3059 scm_bigequal (SCM x, SCM y)
3060 {
3061 int result = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
3062 scm_remember_upto_here_2 (x, y);
3063 return scm_from_bool (0 == result);
3064 }
3065
3066 SCM
3067 scm_real_equalp (SCM x, SCM y)
3068 {
3069 return scm_from_bool (SCM_REAL_VALUE (x) == SCM_REAL_VALUE (y));
3070 }
3071
3072 SCM
3073 scm_complex_equalp (SCM x, SCM y)
3074 {
3075 return scm_from_bool (SCM_COMPLEX_REAL (x) == SCM_COMPLEX_REAL (y)
3076 && SCM_COMPLEX_IMAG (x) == SCM_COMPLEX_IMAG (y));
3077 }
3078
3079 SCM
3080 scm_i_fraction_equalp (SCM x, SCM y)
3081 {
3082 scm_i_fraction_reduce (x);
3083 scm_i_fraction_reduce (y);
3084 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x),
3085 SCM_FRACTION_NUMERATOR (y)))
3086 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x),
3087 SCM_FRACTION_DENOMINATOR (y))))
3088 return SCM_BOOL_F;
3089 else
3090 return SCM_BOOL_T;
3091 }
3092
3093
3094 SCM_DEFINE (scm_number_p, "number?", 1, 0, 0,
3095 (SCM x),
3096 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3097 "otherwise.")
3098 #define FUNC_NAME s_scm_number_p
3099 {
3100 return scm_from_bool (SCM_NUMBERP (x));
3101 }
3102 #undef FUNC_NAME
3103
3104 SCM_DEFINE (scm_complex_p, "complex?", 1, 0, 0,
3105 (SCM x),
3106 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3107 "otherwise. Note that the sets of real, rational and integer\n"
3108 "values form subsets of the set of complex numbers, i. e. the\n"
3109 "predicate will also be fulfilled if @var{x} is a real,\n"
3110 "rational or integer number.")
3111 #define FUNC_NAME s_scm_complex_p
3112 {
3113 /* all numbers are complex. */
3114 return scm_number_p (x);
3115 }
3116 #undef FUNC_NAME
3117
3118 SCM_DEFINE (scm_real_p, "real?", 1, 0, 0,
3119 (SCM x),
3120 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3121 "otherwise. Note that the set of integer values forms a subset of\n"
3122 "the set of real numbers, i. e. the predicate will also be\n"
3123 "fulfilled if @var{x} is an integer number.")
3124 #define FUNC_NAME s_scm_real_p
3125 {
3126 /* we can't represent irrational numbers. */
3127 return scm_rational_p (x);
3128 }
3129 #undef FUNC_NAME
3130
3131 SCM_DEFINE (scm_rational_p, "rational?", 1, 0, 0,
3132 (SCM x),
3133 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3134 "otherwise. Note that the set of integer values forms a subset of\n"
3135 "the set of rational numbers, i. e. the predicate will also be\n"
3136 "fulfilled if @var{x} is an integer number.")
3137 #define FUNC_NAME s_scm_rational_p
3138 {
3139 if (SCM_I_INUMP (x))
3140 return SCM_BOOL_T;
3141 else if (SCM_IMP (x))
3142 return SCM_BOOL_F;
3143 else if (SCM_BIGP (x))
3144 return SCM_BOOL_T;
3145 else if (SCM_FRACTIONP (x))
3146 return SCM_BOOL_T;
3147 else if (SCM_REALP (x))
3148 /* due to their limited precision, all floating point numbers are
3149 rational as well. */
3150 return SCM_BOOL_T;
3151 else
3152 return SCM_BOOL_F;
3153 }
3154 #undef FUNC_NAME
3155
3156 SCM_DEFINE (scm_integer_p, "integer?", 1, 0, 0,
3157 (SCM x),
3158 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3159 "else.")
3160 #define FUNC_NAME s_scm_integer_p
3161 {
3162 double r;
3163 if (SCM_I_INUMP (x))
3164 return SCM_BOOL_T;
3165 if (SCM_IMP (x))
3166 return SCM_BOOL_F;
3167 if (SCM_BIGP (x))
3168 return SCM_BOOL_T;
3169 if (!SCM_INEXACTP (x))
3170 return SCM_BOOL_F;
3171 if (SCM_COMPLEXP (x))
3172 return SCM_BOOL_F;
3173 r = SCM_REAL_VALUE (x);
3174 /* +/-inf passes r==floor(r), making those #t */
3175 if (r == floor (r))
3176 return SCM_BOOL_T;
3177 return SCM_BOOL_F;
3178 }
3179 #undef FUNC_NAME
3180
3181
3182 SCM_DEFINE (scm_inexact_p, "inexact?", 1, 0, 0,
3183 (SCM x),
3184 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3185 "else.")
3186 #define FUNC_NAME s_scm_inexact_p
3187 {
3188 if (SCM_INEXACTP (x))
3189 return SCM_BOOL_T;
3190 if (SCM_NUMBERP (x))
3191 return SCM_BOOL_F;
3192 SCM_WRONG_TYPE_ARG (1, x);
3193 }
3194 #undef FUNC_NAME
3195
3196
3197 SCM_GPROC1 (s_eq_p, "=", scm_tc7_rpsubr, scm_num_eq_p, g_eq_p);
3198 /* "Return @code{#t} if all parameters are numerically equal." */
3199 SCM
3200 scm_num_eq_p (SCM x, SCM y)
3201 {
3202 again:
3203 if (SCM_I_INUMP (x))
3204 {
3205 long xx = SCM_I_INUM (x);
3206 if (SCM_I_INUMP (y))
3207 {
3208 long yy = SCM_I_INUM (y);
3209 return scm_from_bool (xx == yy);
3210 }
3211 else if (SCM_BIGP (y))
3212 return SCM_BOOL_F;
3213 else if (SCM_REALP (y))
3214 return scm_from_bool ((double) xx == SCM_REAL_VALUE (y));
3215 else if (SCM_COMPLEXP (y))
3216 return scm_from_bool (((double) xx == SCM_COMPLEX_REAL (y))
3217 && (0.0 == SCM_COMPLEX_IMAG (y)));
3218 else if (SCM_FRACTIONP (y))
3219 return SCM_BOOL_F;
3220 else
3221 SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
3222 }
3223 else if (SCM_BIGP (x))
3224 {
3225 if (SCM_I_INUMP (y))
3226 return SCM_BOOL_F;
3227 else if (SCM_BIGP (y))
3228 {
3229 int cmp = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
3230 scm_remember_upto_here_2 (x, y);
3231 return scm_from_bool (0 == cmp);
3232 }
3233 else if (SCM_REALP (y))
3234 {
3235 int cmp;
3236 if (xisnan (SCM_REAL_VALUE (y)))
3237 return SCM_BOOL_F;
3238 cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (x), SCM_REAL_VALUE (y));
3239 scm_remember_upto_here_1 (x);
3240 return scm_from_bool (0 == cmp);
3241 }
3242 else if (SCM_COMPLEXP (y))
3243 {
3244 int cmp;
3245 if (0.0 != SCM_COMPLEX_IMAG (y))
3246 return SCM_BOOL_F;
3247 if (xisnan (SCM_COMPLEX_REAL (y)))
3248 return SCM_BOOL_F;
3249 cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (x), SCM_COMPLEX_REAL (y));
3250 scm_remember_upto_here_1 (x);
3251 return scm_from_bool (0 == cmp);
3252 }
3253 else if (SCM_FRACTIONP (y))
3254 return SCM_BOOL_F;
3255 else
3256 SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
3257 }
3258 else if (SCM_REALP (x))
3259 {
3260 if (SCM_I_INUMP (y))
3261 return scm_from_bool (SCM_REAL_VALUE (x) == (double) SCM_I_INUM (y));
3262 else if (SCM_BIGP (y))
3263 {
3264 int cmp;
3265 if (xisnan (SCM_REAL_VALUE (x)))
3266 return SCM_BOOL_F;
3267 cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), SCM_REAL_VALUE (x));
3268 scm_remember_upto_here_1 (y);
3269 return scm_from_bool (0 == cmp);
3270 }
3271 else if (SCM_REALP (y))
3272 return scm_from_bool (SCM_REAL_VALUE (x) == SCM_REAL_VALUE (y));
3273 else if (SCM_COMPLEXP (y))
3274 return scm_from_bool ((SCM_REAL_VALUE (x) == SCM_COMPLEX_REAL (y))
3275 && (0.0 == SCM_COMPLEX_IMAG (y)));
3276 else if (SCM_FRACTIONP (y))
3277 {
3278 double xx = SCM_REAL_VALUE (x);
3279 if (xisnan (xx))
3280 return SCM_BOOL_F;
3281 if (xisinf (xx))
3282 return scm_from_bool (xx < 0.0);
3283 x = scm_inexact_to_exact (x); /* with x as frac or int */
3284 goto again;
3285 }
3286 else
3287 SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
3288 }
3289 else if (SCM_COMPLEXP (x))
3290 {
3291 if (SCM_I_INUMP (y))
3292 return scm_from_bool ((SCM_COMPLEX_REAL (x) == (double) SCM_I_INUM (y))
3293 && (SCM_COMPLEX_IMAG (x) == 0.0));
3294 else if (SCM_BIGP (y))
3295 {
3296 int cmp;
3297 if (0.0 != SCM_COMPLEX_IMAG (x))
3298 return SCM_BOOL_F;
3299 if (xisnan (SCM_COMPLEX_REAL (x)))
3300 return SCM_BOOL_F;
3301 cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), SCM_COMPLEX_REAL (x));
3302 scm_remember_upto_here_1 (y);
3303 return scm_from_bool (0 == cmp);
3304 }
3305 else if (SCM_REALP (y))
3306 return scm_from_bool ((SCM_COMPLEX_REAL (x) == SCM_REAL_VALUE (y))
3307 && (SCM_COMPLEX_IMAG (x) == 0.0));
3308 else if (SCM_COMPLEXP (y))
3309 return scm_from_bool ((SCM_COMPLEX_REAL (x) == SCM_COMPLEX_REAL (y))
3310 && (SCM_COMPLEX_IMAG (x) == SCM_COMPLEX_IMAG (y)));
3311 else if (SCM_FRACTIONP (y))
3312 {
3313 double xx;
3314 if (SCM_COMPLEX_IMAG (x) != 0.0)
3315 return SCM_BOOL_F;
3316 xx = SCM_COMPLEX_REAL (x);
3317 if (xisnan (xx))
3318 return SCM_BOOL_F;
3319 if (xisinf (xx))
3320 return scm_from_bool (xx < 0.0);
3321 x = scm_inexact_to_exact (x); /* with x as frac or int */
3322 goto again;
3323 }
3324 else
3325 SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
3326 }
3327 else if (SCM_FRACTIONP (x))
3328 {
3329 if (SCM_I_INUMP (y))
3330 return SCM_BOOL_F;
3331 else if (SCM_BIGP (y))
3332 return SCM_BOOL_F;
3333 else if (SCM_REALP (y))
3334 {
3335 double yy = SCM_REAL_VALUE (y);
3336 if (xisnan (yy))
3337 return SCM_BOOL_F;
3338 if (xisinf (yy))
3339 return scm_from_bool (0.0 < yy);
3340 y = scm_inexact_to_exact (y); /* with y as frac or int */
3341 goto again;
3342 }
3343 else if (SCM_COMPLEXP (y))
3344 {
3345 double yy;
3346 if (SCM_COMPLEX_IMAG (y) != 0.0)
3347 return SCM_BOOL_F;
3348 yy = SCM_COMPLEX_REAL (y);
3349 if (xisnan (yy))
3350 return SCM_BOOL_F;
3351 if (xisinf (yy))
3352 return scm_from_bool (0.0 < yy);
3353 y = scm_inexact_to_exact (y); /* with y as frac or int */
3354 goto again;
3355 }
3356 else if (SCM_FRACTIONP (y))
3357 return scm_i_fraction_equalp (x, y);
3358 else
3359 SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
3360 }
3361 else
3362 SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARG1, s_eq_p);
3363 }
3364
3365
3366 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3367 done are good for inums, but for bignums an answer can almost always be
3368 had by just examining a few high bits of the operands, as done by GMP in
3369 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3370 of the float exponent to take into account. */
3371
3372 SCM_GPROC1 (s_less_p, "<", scm_tc7_rpsubr, scm_less_p, g_less_p);
3373 /* "Return @code{#t} if the list of parameters is monotonically\n"
3374 * "increasing."
3375 */
3376 SCM
3377 scm_less_p (SCM x, SCM y)
3378 {
3379 again:
3380 if (SCM_I_INUMP (x))
3381 {
3382 long xx = SCM_I_INUM (x);
3383 if (SCM_I_INUMP (y))
3384 {
3385 long yy = SCM_I_INUM (y);
3386 return scm_from_bool (xx < yy);
3387 }
3388 else if (SCM_BIGP (y))
3389 {
3390 int sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
3391 scm_remember_upto_here_1 (y);
3392 return scm_from_bool (sgn > 0);
3393 }
3394 else if (SCM_REALP (y))
3395 return scm_from_bool ((double) xx < SCM_REAL_VALUE (y));
3396 else if (SCM_FRACTIONP (y))
3397 {
3398 /* "x < a/b" becomes "x*b < a" */
3399 int_frac:
3400 x = scm_product (x, SCM_FRACTION_DENOMINATOR (y));
3401 y = SCM_FRACTION_NUMERATOR (y);
3402 goto again;
3403 }
3404 else
3405 SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
3406 }
3407 else if (SCM_BIGP (x))
3408 {
3409 if (SCM_I_INUMP (y))
3410 {
3411 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
3412 scm_remember_upto_here_1 (x);
3413 return scm_from_bool (sgn < 0);
3414 }
3415 else if (SCM_BIGP (y))
3416 {
3417 int cmp = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
3418 scm_remember_upto_here_2 (x, y);
3419 return scm_from_bool (cmp < 0);
3420 }
3421 else if (SCM_REALP (y))
3422 {
3423 int cmp;
3424 if (xisnan (SCM_REAL_VALUE (y)))
3425 return SCM_BOOL_F;
3426 cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (x), SCM_REAL_VALUE (y));
3427 scm_remember_upto_here_1 (x);
3428 return scm_from_bool (cmp < 0);
3429 }
3430 else if (SCM_FRACTIONP (y))
3431 goto int_frac;
3432 else
3433 SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
3434 }
3435 else if (SCM_REALP (x))
3436 {
3437 if (SCM_I_INUMP (y))
3438 return scm_from_bool (SCM_REAL_VALUE (x) < (double) SCM_I_INUM (y));
3439 else if (SCM_BIGP (y))
3440 {
3441 int cmp;
3442 if (xisnan (SCM_REAL_VALUE (x)))
3443 return SCM_BOOL_F;
3444 cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), SCM_REAL_VALUE (x));
3445 scm_remember_upto_here_1 (y);
3446 return scm_from_bool (cmp > 0);
3447 }
3448 else if (SCM_REALP (y))
3449 return scm_from_bool (SCM_REAL_VALUE (x) < SCM_REAL_VALUE (y));
3450 else if (SCM_FRACTIONP (y))
3451 {
3452 double xx = SCM_REAL_VALUE (x);
3453 if (xisnan (xx))
3454 return SCM_BOOL_F;
3455 if (xisinf (xx))
3456 return scm_from_bool (xx < 0.0);
3457 x = scm_inexact_to_exact (x); /* with x as frac or int */
3458 goto again;
3459 }
3460 else
3461 SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
3462 }
3463 else if (SCM_FRACTIONP (x))
3464 {
3465 if (SCM_I_INUMP (y) || SCM_BIGP (y))
3466 {
3467 /* "a/b < y" becomes "a < y*b" */
3468 y = scm_product (y, SCM_FRACTION_DENOMINATOR (x));
3469 x = SCM_FRACTION_NUMERATOR (x);
3470 goto again;
3471 }
3472 else if (SCM_REALP (y))
3473 {
3474 double yy = SCM_REAL_VALUE (y);
3475 if (xisnan (yy))
3476 return SCM_BOOL_F;
3477 if (xisinf (yy))
3478 return scm_from_bool (0.0 < yy);
3479 y = scm_inexact_to_exact (y); /* with y as frac or int */
3480 goto again;
3481 }
3482 else if (SCM_FRACTIONP (y))
3483 {
3484 /* "a/b < c/d" becomes "a*d < c*b" */
3485 SCM new_x = scm_product (SCM_FRACTION_NUMERATOR (x),
3486 SCM_FRACTION_DENOMINATOR (y));
3487 SCM new_y = scm_product (SCM_FRACTION_NUMERATOR (y),
3488 SCM_FRACTION_DENOMINATOR (x));
3489 x = new_x;
3490 y = new_y;
3491 goto again;
3492 }
3493 else
3494 SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
3495 }
3496 else
3497 SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARG1, s_less_p);
3498 }
3499
3500
3501 SCM_GPROC1 (s_scm_gr_p, ">", scm_tc7_rpsubr, scm_gr_p, g_gr_p);
3502 /* "Return @code{#t} if the list of parameters is monotonically\n"
3503 * "decreasing."
3504 */
3505 #define FUNC_NAME s_scm_gr_p
3506 SCM
3507 scm_gr_p (SCM x, SCM y)
3508 {
3509 if (!SCM_NUMBERP (x))
3510 SCM_WTA_DISPATCH_2 (g_gr_p, x, y, SCM_ARG1, FUNC_NAME);
3511 else if (!SCM_NUMBERP (y))
3512 SCM_WTA_DISPATCH_2 (g_gr_p, x, y, SCM_ARG2, FUNC_NAME);
3513 else
3514 return scm_less_p (y, x);
3515 }
3516 #undef FUNC_NAME
3517
3518
3519 SCM_GPROC1 (s_scm_leq_p, "<=", scm_tc7_rpsubr, scm_leq_p, g_leq_p);
3520 /* "Return @code{#t} if the list of parameters is monotonically\n"
3521 * "non-decreasing."
3522 */
3523 #define FUNC_NAME s_scm_leq_p
3524 SCM
3525 scm_leq_p (SCM x, SCM y)
3526 {
3527 if (!SCM_NUMBERP (x))
3528 SCM_WTA_DISPATCH_2 (g_leq_p, x, y, SCM_ARG1, FUNC_NAME);
3529 else if (!SCM_NUMBERP (y))
3530 SCM_WTA_DISPATCH_2 (g_leq_p, x, y, SCM_ARG2, FUNC_NAME);
3531 else if (scm_is_true (scm_nan_p (x)) || scm_is_true (scm_nan_p (y)))
3532 return SCM_BOOL_F;
3533 else
3534 return scm_not (scm_less_p (y, x));
3535 }
3536 #undef FUNC_NAME
3537
3538
3539 SCM_GPROC1 (s_scm_geq_p, ">=", scm_tc7_rpsubr, scm_geq_p, g_geq_p);
3540 /* "Return @code{#t} if the list of parameters is monotonically\n"
3541 * "non-increasing."
3542 */
3543 #define FUNC_NAME s_scm_geq_p
3544 SCM
3545 scm_geq_p (SCM x, SCM y)
3546 {
3547 if (!SCM_NUMBERP (x))
3548 SCM_WTA_DISPATCH_2 (g_geq_p, x, y, SCM_ARG1, FUNC_NAME);
3549 else if (!SCM_NUMBERP (y))
3550 SCM_WTA_DISPATCH_2 (g_geq_p, x, y, SCM_ARG2, FUNC_NAME);
3551 else if (scm_is_true (scm_nan_p (x)) || scm_is_true (scm_nan_p (y)))
3552 return SCM_BOOL_F;
3553 else
3554 return scm_not (scm_less_p (x, y));
3555 }
3556 #undef FUNC_NAME
3557
3558
3559 SCM_GPROC (s_zero_p, "zero?", 1, 0, 0, scm_zero_p, g_zero_p);
3560 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3561 * "zero."
3562 */
3563 SCM
3564 scm_zero_p (SCM z)
3565 {
3566 if (SCM_I_INUMP (z))
3567 return scm_from_bool (scm_is_eq (z, SCM_INUM0));
3568 else if (SCM_BIGP (z))
3569 return SCM_BOOL_F;
3570 else if (SCM_REALP (z))
3571 return scm_from_bool (SCM_REAL_VALUE (z) == 0.0);
3572 else if (SCM_COMPLEXP (z))
3573 return scm_from_bool (SCM_COMPLEX_REAL (z) == 0.0
3574 && SCM_COMPLEX_IMAG (z) == 0.0);
3575 else if (SCM_FRACTIONP (z))
3576 return SCM_BOOL_F;
3577 else
3578 SCM_WTA_DISPATCH_1 (g_zero_p, z, SCM_ARG1, s_zero_p);
3579 }
3580
3581
3582 SCM_GPROC (s_positive_p, "positive?", 1, 0, 0, scm_positive_p, g_positive_p);
3583 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3584 * "zero."
3585 */
3586 SCM
3587 scm_positive_p (SCM x)
3588 {
3589 if (SCM_I_INUMP (x))
3590 return scm_from_bool (SCM_I_INUM (x) > 0);
3591 else if (SCM_BIGP (x))
3592 {
3593 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
3594 scm_remember_upto_here_1 (x);
3595 return scm_from_bool (sgn > 0);
3596 }
3597 else if (SCM_REALP (x))
3598 return scm_from_bool(SCM_REAL_VALUE (x) > 0.0);
3599 else if (SCM_FRACTIONP (x))
3600 return scm_positive_p (SCM_FRACTION_NUMERATOR (x));
3601 else
3602 SCM_WTA_DISPATCH_1 (g_positive_p, x, SCM_ARG1, s_positive_p);
3603 }
3604
3605
3606 SCM_GPROC (s_negative_p, "negative?", 1, 0, 0, scm_negative_p, g_negative_p);
3607 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3608 * "zero."
3609 */
3610 SCM
3611 scm_negative_p (SCM x)
3612 {
3613 if (SCM_I_INUMP (x))
3614 return scm_from_bool (SCM_I_INUM (x) < 0);
3615 else if (SCM_BIGP (x))
3616 {
3617 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
3618 scm_remember_upto_here_1 (x);
3619 return scm_from_bool (sgn < 0);
3620 }
3621 else if (SCM_REALP (x))
3622 return scm_from_bool(SCM_REAL_VALUE (x) < 0.0);
3623 else if (SCM_FRACTIONP (x))
3624 return scm_negative_p (SCM_FRACTION_NUMERATOR (x));
3625 else
3626 SCM_WTA_DISPATCH_1 (g_negative_p, x, SCM_ARG1, s_negative_p);
3627 }
3628
3629
3630 /* scm_min and scm_max return an inexact when either argument is inexact, as
3631 required by r5rs. On that basis, for exact/inexact combinations the
3632 exact is converted to inexact to compare and possibly return. This is
3633 unlike scm_less_p above which takes some trouble to preserve all bits in
3634 its test, such trouble is not required for min and max. */
3635
3636 SCM_GPROC1 (s_max, "max", scm_tc7_asubr, scm_max, g_max);
3637 /* "Return the maximum of all parameter values."
3638 */
3639 SCM
3640 scm_max (SCM x, SCM y)
3641 {
3642 if (SCM_UNBNDP (y))
3643 {
3644 if (SCM_UNBNDP (x))
3645 SCM_WTA_DISPATCH_0 (g_max, s_max);
3646 else if (SCM_I_INUMP(x) || SCM_BIGP(x) || SCM_REALP(x) || SCM_FRACTIONP(x))
3647 return x;
3648 else
3649 SCM_WTA_DISPATCH_1 (g_max, x, SCM_ARG1, s_max);
3650 }
3651
3652 if (SCM_I_INUMP (x))
3653 {
3654 long xx = SCM_I_INUM (x);
3655 if (SCM_I_INUMP (y))
3656 {
3657 long yy = SCM_I_INUM (y);
3658 return (xx < yy) ? y : x;
3659 }
3660 else if (SCM_BIGP (y))
3661 {
3662 int sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
3663 scm_remember_upto_here_1 (y);
3664 return (sgn < 0) ? x : y;
3665 }
3666 else if (SCM_REALP (y))
3667 {
3668 double z = xx;
3669 /* if y==NaN then ">" is false and we return NaN */
3670 return (z > SCM_REAL_VALUE (y)) ? scm_from_double (z) : y;
3671 }
3672 else if (SCM_FRACTIONP (y))
3673 {
3674 use_less:
3675 return (scm_is_false (scm_less_p (x, y)) ? x : y);
3676 }
3677 else
3678 SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
3679 }
3680 else if (SCM_BIGP (x))
3681 {
3682 if (SCM_I_INUMP (y))
3683 {
3684 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
3685 scm_remember_upto_here_1 (x);
3686 return (sgn < 0) ? y : x;
3687 }
3688 else if (SCM_BIGP (y))
3689 {
3690 int cmp = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
3691 scm_remember_upto_here_2 (x, y);
3692 return (cmp > 0) ? x : y;
3693 }
3694 else if (SCM_REALP (y))
3695 {
3696 /* if y==NaN then xx>yy is false, so we return the NaN y */
3697 double xx, yy;
3698 big_real:
3699 xx = scm_i_big2dbl (x);
3700 yy = SCM_REAL_VALUE (y);
3701 return (xx > yy ? scm_from_double (xx) : y);
3702 }
3703 else if (SCM_FRACTIONP (y))
3704 {
3705 goto use_less;
3706 }
3707 else
3708 SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
3709 }
3710 else if (SCM_REALP (x))
3711 {
3712 if (SCM_I_INUMP (y))
3713 {
3714 double z = SCM_I_INUM (y);
3715 /* if x==NaN then "<" is false and we return NaN */
3716 return (SCM_REAL_VALUE (x) < z) ? scm_from_double (z) : x;
3717 }
3718 else if (SCM_BIGP (y))
3719 {
3720 SCM_SWAP (x, y);
3721 goto big_real;
3722 }
3723 else if (SCM_REALP (y))
3724 {
3725 /* if x==NaN then our explicit check means we return NaN
3726 if y==NaN then ">" is false and we return NaN
3727 calling isnan is unavoidable, since it's the only way to know
3728 which of x or y causes any compares to be false */
3729 double xx = SCM_REAL_VALUE (x);
3730 return (xisnan (xx) || xx > SCM_REAL_VALUE (y)) ? x : y;
3731 }
3732 else if (SCM_FRACTIONP (y))
3733 {
3734 double yy = scm_i_fraction2double (y);
3735 double xx = SCM_REAL_VALUE (x);
3736 return (xx < yy) ? scm_from_double (yy) : x;
3737 }
3738 else
3739 SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
3740 }
3741 else if (SCM_FRACTIONP (x))
3742 {
3743 if (SCM_I_INUMP (y))
3744 {
3745 goto use_less;
3746 }
3747 else if (SCM_BIGP (y))
3748 {
3749 goto use_less;
3750 }
3751 else if (SCM_REALP (y))
3752 {
3753 double xx = scm_i_fraction2double (x);
3754 return (xx < SCM_REAL_VALUE (y)) ? y : scm_from_double (xx);
3755 }
3756 else if (SCM_FRACTIONP (y))
3757 {
3758 goto use_less;
3759 }
3760 else
3761 SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
3762 }
3763 else
3764 SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARG1, s_max);
3765 }
3766
3767
3768 SCM_GPROC1 (s_min, "min", scm_tc7_asubr, scm_min, g_min);
3769 /* "Return the minium of all parameter values."
3770 */
3771 SCM
3772 scm_min (SCM x, SCM y)
3773 {
3774 if (SCM_UNBNDP (y))
3775 {
3776 if (SCM_UNBNDP (x))
3777 SCM_WTA_DISPATCH_0 (g_min, s_min);
3778 else if (SCM_I_INUMP(x) || SCM_BIGP(x) || SCM_REALP(x) || SCM_FRACTIONP(x))
3779 return x;
3780 else
3781 SCM_WTA_DISPATCH_1 (g_min, x, SCM_ARG1, s_min);
3782 }
3783
3784 if (SCM_I_INUMP (x))
3785 {
3786 long xx = SCM_I_INUM (x);
3787 if (SCM_I_INUMP (y))
3788 {
3789 long yy = SCM_I_INUM (y);
3790 return (xx < yy) ? x : y;
3791 }
3792 else if (SCM_BIGP (y))
3793 {
3794 int sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
3795 scm_remember_upto_here_1 (y);
3796 return (sgn < 0) ? y : x;
3797 }
3798 else if (SCM_REALP (y))
3799 {
3800 double z = xx;
3801 /* if y==NaN then "<" is false and we return NaN */
3802 return (z < SCM_REAL_VALUE (y)) ? scm_from_double (z) : y;
3803 }
3804 else if (SCM_FRACTIONP (y))
3805 {
3806 use_less:
3807 return (scm_is_false (scm_less_p (x, y)) ? y : x);
3808 }
3809 else
3810 SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min);
3811 }
3812 else if (SCM_BIGP (x))
3813 {
3814 if (SCM_I_INUMP (y))
3815 {
3816 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
3817 scm_remember_upto_here_1 (x);
3818 return (sgn < 0) ? x : y;
3819 }
3820 else if (SCM_BIGP (y))
3821 {
3822 int cmp = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
3823 scm_remember_upto_here_2 (x, y);
3824 return (cmp > 0) ? y : x;
3825 }
3826 else if (SCM_REALP (y))
3827 {
3828 /* if y==NaN then xx<yy is false, so we return the NaN y */
3829 double xx, yy;
3830 big_real:
3831 xx = scm_i_big2dbl (x);
3832 yy = SCM_REAL_VALUE (y);
3833 return (xx < yy ? scm_from_double (xx) : y);
3834 }
3835 else if (SCM_FRACTIONP (y))
3836 {
3837 goto use_less;
3838 }
3839 else
3840 SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min);
3841 }
3842 else if (SCM_REALP (x))
3843 {
3844 if (SCM_I_INUMP (y))
3845 {
3846 double z = SCM_I_INUM (y);
3847 /* if x==NaN then "<" is false and we return NaN */
3848 return (z < SCM_REAL_VALUE (x)) ? scm_from_double (z) : x;
3849 }
3850 else if (SCM_BIGP (y))
3851 {
3852 SCM_SWAP (x, y);
3853 goto big_real;
3854 }
3855 else if (SCM_REALP (y))
3856 {
3857 /* if x==NaN then our explicit check means we return NaN
3858 if y==NaN then "<" is false and we return NaN
3859 calling isnan is unavoidable, since it's the only way to know
3860 which of x or y causes any compares to be false */
3861 double xx = SCM_REAL_VALUE (x);
3862 return (xisnan (xx) || xx < SCM_REAL_VALUE (y)) ? x : y;
3863 }
3864 else if (SCM_FRACTIONP (y))
3865 {
3866 double yy = scm_i_fraction2double (y);
3867 double xx = SCM_REAL_VALUE (x);
3868 return (yy < xx) ? scm_from_double (yy) : x;
3869 }
3870 else
3871 SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min);
3872 }
3873 else if (SCM_FRACTIONP (x))
3874 {
3875 if (SCM_I_INUMP (y))
3876 {
3877 goto use_less;
3878 }
3879 else if (SCM_BIGP (y))
3880 {
3881 goto use_less;
3882 }
3883 else if (SCM_REALP (y))
3884 {
3885 double xx = scm_i_fraction2double (x);
3886 return (SCM_REAL_VALUE (y) < xx) ? y : scm_from_double (xx);
3887 }
3888 else if (SCM_FRACTIONP (y))
3889 {
3890 goto use_less;
3891 }
3892 else
3893 SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
3894 }
3895 else
3896 SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARG1, s_min);
3897 }
3898
3899
3900 SCM_GPROC1 (s_sum, "+", scm_tc7_asubr, scm_sum, g_sum);
3901 /* "Return the sum of all parameter values. Return 0 if called without\n"
3902 * "any parameters."
3903 */
3904 SCM
3905 scm_sum (SCM x, SCM y)
3906 {
3907 if (SCM_UNBNDP (y))
3908 {
3909 if (SCM_NUMBERP (x)) return x;
3910 if (SCM_UNBNDP (x)) return SCM_INUM0;
3911 SCM_WTA_DISPATCH_1 (g_sum, x, SCM_ARG1, s_sum);
3912 }
3913
3914 if (SCM_I_INUMP (x))
3915 {
3916 if (SCM_I_INUMP (y))
3917 {
3918 long xx = SCM_I_INUM (x);
3919 long yy = SCM_I_INUM (y);
3920 long int z = xx + yy;
3921 return SCM_FIXABLE (z) ? SCM_I_MAKINUM (z) : scm_i_long2big (z);
3922 }
3923 else if (SCM_BIGP (y))
3924 {
3925 SCM_SWAP (x, y);
3926 goto add_big_inum;
3927 }
3928 else if (SCM_REALP (y))
3929 {
3930 long int xx = SCM_I_INUM (x);
3931 return scm_from_double (xx + SCM_REAL_VALUE (y));
3932 }
3933 else if (SCM_COMPLEXP (y))
3934 {
3935 long int xx = SCM_I_INUM (x);
3936 return scm_c_make_rectangular (xx + SCM_COMPLEX_REAL (y),
3937 SCM_COMPLEX_IMAG (y));
3938 }
3939 else if (SCM_FRACTIONP (y))
3940 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y),
3941 scm_product (x, SCM_FRACTION_DENOMINATOR (y))),
3942 SCM_FRACTION_DENOMINATOR (y));
3943 else
3944 SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
3945 } else if (SCM_BIGP (x))
3946 {
3947 if (SCM_I_INUMP (y))
3948 {
3949 long int inum;
3950 int bigsgn;
3951 add_big_inum:
3952 inum = SCM_I_INUM (y);
3953 if (inum == 0)
3954 return x;
3955 bigsgn = mpz_sgn (SCM_I_BIG_MPZ (x));
3956 if (inum < 0)
3957 {
3958 SCM result = scm_i_mkbig ();
3959 mpz_sub_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), - inum);
3960 scm_remember_upto_here_1 (x);
3961 /* we know the result will have to be a bignum */
3962 if (bigsgn == -1)
3963 return result;
3964 return scm_i_normbig (result);
3965 }
3966 else
3967 {
3968 SCM result = scm_i_mkbig ();
3969 mpz_add_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), inum);
3970 scm_remember_upto_here_1 (x);
3971 /* we know the result will have to be a bignum */
3972 if (bigsgn == 1)
3973 return result;
3974 return scm_i_normbig (result);
3975 }
3976 }
3977 else if (SCM_BIGP (y))
3978 {
3979 SCM result = scm_i_mkbig ();
3980 int sgn_x = mpz_sgn (SCM_I_BIG_MPZ (x));
3981 int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y));
3982 mpz_add (SCM_I_BIG_MPZ (result),
3983 SCM_I_BIG_MPZ (x),
3984 SCM_I_BIG_MPZ (y));
3985 scm_remember_upto_here_2 (x, y);
3986 /* we know the result will have to be a bignum */
3987 if (sgn_x == sgn_y)
3988 return result;
3989 return scm_i_normbig (result);
3990 }
3991 else if (SCM_REALP (y))
3992 {
3993 double result = mpz_get_d (SCM_I_BIG_MPZ (x)) + SCM_REAL_VALUE (y);
3994 scm_remember_upto_here_1 (x);
3995 return scm_from_double (result);
3996 }
3997 else if (SCM_COMPLEXP (y))
3998 {
3999 double real_part = (mpz_get_d (SCM_I_BIG_MPZ (x))
4000 + SCM_COMPLEX_REAL (y));
4001 scm_remember_upto_here_1 (x);
4002 return scm_c_make_rectangular (real_part, SCM_COMPLEX_IMAG (y));
4003 }
4004 else if (SCM_FRACTIONP (y))
4005 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y),
4006 scm_product (x, SCM_FRACTION_DENOMINATOR (y))),
4007 SCM_FRACTION_DENOMINATOR (y));
4008 else
4009 SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
4010 }
4011 else if (SCM_REALP (x))
4012 {
4013 if (SCM_I_INUMP (y))
4014 return scm_from_double (SCM_REAL_VALUE (x) + SCM_I_INUM (y));
4015 else if (SCM_BIGP (y))
4016 {
4017 double result = mpz_get_d (SCM_I_BIG_MPZ (y)) + SCM_REAL_VALUE (x);
4018 scm_remember_upto_here_1 (y);
4019 return scm_from_double (result);
4020 }
4021 else if (SCM_REALP (y))
4022 return scm_from_double (SCM_REAL_VALUE (x) + SCM_REAL_VALUE (y));
4023 else if (SCM_COMPLEXP (y))
4024 return scm_c_make_rectangular (SCM_REAL_VALUE (x) + SCM_COMPLEX_REAL (y),
4025 SCM_COMPLEX_IMAG (y));
4026 else if (SCM_FRACTIONP (y))
4027 return scm_from_double (SCM_REAL_VALUE (x) + scm_i_fraction2double (y));
4028 else
4029 SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
4030 }
4031 else if (SCM_COMPLEXP (x))
4032 {
4033 if (SCM_I_INUMP (y))
4034 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) + SCM_I_INUM (y),
4035 SCM_COMPLEX_IMAG (x));
4036 else if (SCM_BIGP (y))
4037 {
4038 double real_part = (mpz_get_d (SCM_I_BIG_MPZ (y))
4039 + SCM_COMPLEX_REAL (x));
4040 scm_remember_upto_here_1 (y);
4041 return scm_c_make_rectangular (real_part, SCM_COMPLEX_IMAG (x));
4042 }
4043 else if (SCM_REALP (y))
4044 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) + SCM_REAL_VALUE (y),
4045 SCM_COMPLEX_IMAG (x));
4046 else if (SCM_COMPLEXP (y))
4047 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) + SCM_COMPLEX_REAL (y),
4048 SCM_COMPLEX_IMAG (x) + SCM_COMPLEX_IMAG (y));
4049 else if (SCM_FRACTIONP (y))
4050 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) + scm_i_fraction2double (y),
4051 SCM_COMPLEX_IMAG (x));
4052 else
4053 SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
4054 }
4055 else if (SCM_FRACTIONP (x))
4056 {
4057 if (SCM_I_INUMP (y))
4058 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x),
4059 scm_product (y, SCM_FRACTION_DENOMINATOR (x))),
4060 SCM_FRACTION_DENOMINATOR (x));
4061 else if (SCM_BIGP (y))
4062 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x),
4063 scm_product (y, SCM_FRACTION_DENOMINATOR (x))),
4064 SCM_FRACTION_DENOMINATOR (x));
4065 else if (SCM_REALP (y))
4066 return scm_from_double (SCM_REAL_VALUE (y) + scm_i_fraction2double (x));
4067 else if (SCM_COMPLEXP (y))
4068 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y) + scm_i_fraction2double (x),
4069 SCM_COMPLEX_IMAG (y));
4070 else if (SCM_FRACTIONP (y))
4071 /* a/b + c/d = (ad + bc) / bd */
4072 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x), SCM_FRACTION_DENOMINATOR (y)),
4073 scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x))),
4074 scm_product (SCM_FRACTION_DENOMINATOR (x), SCM_FRACTION_DENOMINATOR (y)));
4075 else
4076 SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
4077 }
4078 else
4079 SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARG1, s_sum);
4080 }
4081
4082
4083 SCM_GPROC1 (s_difference, "-", scm_tc7_asubr, scm_difference, g_difference);
4084 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4085 * the sum of all but the first argument are subtracted from the first
4086 * argument. */
4087 #define FUNC_NAME s_difference
4088 SCM
4089 scm_difference (SCM x, SCM y)
4090 {
4091 if (SCM_UNBNDP (y))
4092 {
4093 if (SCM_UNBNDP (x))
4094 SCM_WTA_DISPATCH_0 (g_difference, s_difference);
4095 else
4096 if (SCM_I_INUMP (x))
4097 {
4098 long xx = -SCM_I_INUM (x);
4099 if (SCM_FIXABLE (xx))
4100 return SCM_I_MAKINUM (xx);
4101 else
4102 return scm_i_long2big (xx);
4103 }
4104 else if (SCM_BIGP (x))
4105 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4106 bignum, but negating that gives a fixnum. */
4107 return scm_i_normbig (scm_i_clonebig (x, 0));
4108 else if (SCM_REALP (x))
4109 return scm_from_double (-SCM_REAL_VALUE (x));
4110 else if (SCM_COMPLEXP (x))
4111 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x),
4112 -SCM_COMPLEX_IMAG (x));
4113 else if (SCM_FRACTIONP (x))
4114 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED),
4115 SCM_FRACTION_DENOMINATOR (x));
4116 else
4117 SCM_WTA_DISPATCH_1 (g_difference, x, SCM_ARG1, s_difference);
4118 }
4119
4120 if (SCM_I_INUMP (x))
4121 {
4122 if (SCM_I_INUMP (y))
4123 {
4124 long int xx = SCM_I_INUM (x);
4125 long int yy = SCM_I_INUM (y);
4126 long int z = xx - yy;
4127 if (SCM_FIXABLE (z))
4128 return SCM_I_MAKINUM (z);
4129 else
4130 return scm_i_long2big (z);
4131 }
4132 else if (SCM_BIGP (y))
4133 {
4134 /* inum-x - big-y */
4135 long xx = SCM_I_INUM (x);
4136
4137 if (xx == 0)
4138 return scm_i_clonebig (y, 0);
4139 else
4140 {
4141 int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y));
4142 SCM result = scm_i_mkbig ();
4143
4144 if (xx >= 0)
4145 mpz_ui_sub (SCM_I_BIG_MPZ (result), xx, SCM_I_BIG_MPZ (y));
4146 else
4147 {
4148 /* x - y == -(y + -x) */
4149 mpz_add_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (y), -xx);
4150 mpz_neg (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result));
4151 }
4152 scm_remember_upto_here_1 (y);
4153
4154 if ((xx < 0 && (sgn_y > 0)) || ((xx > 0) && sgn_y < 0))
4155 /* we know the result will have to be a bignum */
4156 return result;
4157 else
4158 return scm_i_normbig (result);
4159 }
4160 }
4161 else if (SCM_REALP (y))
4162 {
4163 long int xx = SCM_I_INUM (x);
4164 return scm_from_double (xx - SCM_REAL_VALUE (y));
4165 }
4166 else if (SCM_COMPLEXP (y))
4167 {
4168 long int xx = SCM_I_INUM (x);
4169 return scm_c_make_rectangular (xx - SCM_COMPLEX_REAL (y),
4170 - SCM_COMPLEX_IMAG (y));
4171 }
4172 else if (SCM_FRACTIONP (y))
4173 /* a - b/c = (ac - b) / c */
4174 return scm_i_make_ratio (scm_difference (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
4175 SCM_FRACTION_NUMERATOR (y)),
4176 SCM_FRACTION_DENOMINATOR (y));
4177 else
4178 SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
4179 }
4180 else if (SCM_BIGP (x))
4181 {
4182 if (SCM_I_INUMP (y))
4183 {
4184 /* big-x - inum-y */
4185 long yy = SCM_I_INUM (y);
4186 int sgn_x = mpz_sgn (SCM_I_BIG_MPZ (x));
4187
4188 scm_remember_upto_here_1 (x);
4189 if (sgn_x == 0)
4190 return (SCM_FIXABLE (-yy) ?
4191 SCM_I_MAKINUM (-yy) : scm_from_long (-yy));
4192 else
4193 {
4194 SCM result = scm_i_mkbig ();
4195
4196 if (yy >= 0)
4197 mpz_sub_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), yy);
4198 else
4199 mpz_add_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), -yy);
4200 scm_remember_upto_here_1 (x);
4201
4202 if ((sgn_x < 0 && (yy > 0)) || ((sgn_x > 0) && yy < 0))
4203 /* we know the result will have to be a bignum */
4204 return result;
4205 else
4206 return scm_i_normbig (result);
4207 }
4208 }
4209 else if (SCM_BIGP (y))
4210 {
4211 int sgn_x = mpz_sgn (SCM_I_BIG_MPZ (x));
4212 int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y));
4213 SCM result = scm_i_mkbig ();
4214 mpz_sub (SCM_I_BIG_MPZ (result),
4215 SCM_I_BIG_MPZ (x),
4216 SCM_I_BIG_MPZ (y));
4217 scm_remember_upto_here_2 (x, y);
4218 /* we know the result will have to be a bignum */
4219 if ((sgn_x == 1) && (sgn_y == -1))
4220 return result;
4221 if ((sgn_x == -1) && (sgn_y == 1))
4222 return result;
4223 return scm_i_normbig (result);
4224 }
4225 else if (SCM_REALP (y))
4226 {
4227 double result = mpz_get_d (SCM_I_BIG_MPZ (x)) - SCM_REAL_VALUE (y);
4228 scm_remember_upto_here_1 (x);
4229 return scm_from_double (result);
4230 }
4231 else if (SCM_COMPLEXP (y))
4232 {
4233 double real_part = (mpz_get_d (SCM_I_BIG_MPZ (x))
4234 - SCM_COMPLEX_REAL (y));
4235 scm_remember_upto_here_1 (x);
4236 return scm_c_make_rectangular (real_part, - SCM_COMPLEX_IMAG (y));
4237 }
4238 else if (SCM_FRACTIONP (y))
4239 return scm_i_make_ratio (scm_difference (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
4240 SCM_FRACTION_NUMERATOR (y)),
4241 SCM_FRACTION_DENOMINATOR (y));
4242 else SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
4243 }
4244 else if (SCM_REALP (x))
4245 {
4246 if (SCM_I_INUMP (y))
4247 return scm_from_double (SCM_REAL_VALUE (x) - SCM_I_INUM (y));
4248 else if (SCM_BIGP (y))
4249 {
4250 double result = SCM_REAL_VALUE (x) - mpz_get_d (SCM_I_BIG_MPZ (y));
4251 scm_remember_upto_here_1 (x);
4252 return scm_from_double (result);
4253 }
4254 else if (SCM_REALP (y))
4255 return scm_from_double (SCM_REAL_VALUE (x) - SCM_REAL_VALUE (y));
4256 else if (SCM_COMPLEXP (y))
4257 return scm_c_make_rectangular (SCM_REAL_VALUE (x) - SCM_COMPLEX_REAL (y),
4258 -SCM_COMPLEX_IMAG (y));
4259 else if (SCM_FRACTIONP (y))
4260 return scm_from_double (SCM_REAL_VALUE (x) - scm_i_fraction2double (y));
4261 else
4262 SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
4263 }
4264 else if (SCM_COMPLEXP (x))
4265 {
4266 if (SCM_I_INUMP (y))
4267 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) - SCM_I_INUM (y),
4268 SCM_COMPLEX_IMAG (x));
4269 else if (SCM_BIGP (y))
4270 {
4271 double real_part = (SCM_COMPLEX_REAL (x)
4272 - mpz_get_d (SCM_I_BIG_MPZ (y)));
4273 scm_remember_upto_here_1 (x);
4274 return scm_c_make_rectangular (real_part, SCM_COMPLEX_IMAG (y));
4275 }
4276 else if (SCM_REALP (y))
4277 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) - SCM_REAL_VALUE (y),
4278 SCM_COMPLEX_IMAG (x));
4279 else if (SCM_COMPLEXP (y))
4280 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) - SCM_COMPLEX_REAL (y),
4281 SCM_COMPLEX_IMAG (x) - SCM_COMPLEX_IMAG (y));
4282 else if (SCM_FRACTIONP (y))
4283 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) - scm_i_fraction2double (y),
4284 SCM_COMPLEX_IMAG (x));
4285 else
4286 SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
4287 }
4288 else if (SCM_FRACTIONP (x))
4289 {
4290 if (SCM_I_INUMP (y))
4291 /* a/b - c = (a - cb) / b */
4292 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x),
4293 scm_product(y, SCM_FRACTION_DENOMINATOR (x))),
4294 SCM_FRACTION_DENOMINATOR (x));
4295 else if (SCM_BIGP (y))
4296 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x),
4297 scm_product(y, SCM_FRACTION_DENOMINATOR (x))),
4298 SCM_FRACTION_DENOMINATOR (x));
4299 else if (SCM_REALP (y))
4300 return scm_from_double (scm_i_fraction2double (x) - SCM_REAL_VALUE (y));
4301 else if (SCM_COMPLEXP (y))
4302 return scm_c_make_rectangular (scm_i_fraction2double (x) - SCM_COMPLEX_REAL (y),
4303 -SCM_COMPLEX_IMAG (y));
4304 else if (SCM_FRACTIONP (y))
4305 /* a/b - c/d = (ad - bc) / bd */
4306 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x), SCM_FRACTION_DENOMINATOR (y)),
4307 scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x))),
4308 scm_product (SCM_FRACTION_DENOMINATOR (x), SCM_FRACTION_DENOMINATOR (y)));
4309 else
4310 SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
4311 }
4312 else
4313 SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARG1, s_difference);
4314 }
4315 #undef FUNC_NAME
4316
4317
4318 SCM_GPROC1 (s_product, "*", scm_tc7_asubr, scm_product, g_product);
4319 /* "Return the product of all arguments. If called without arguments,\n"
4320 * "1 is returned."
4321 */
4322 SCM
4323 scm_product (SCM x, SCM y)
4324 {
4325 if (SCM_UNBNDP (y))
4326 {
4327 if (SCM_UNBNDP (x))
4328 return SCM_I_MAKINUM (1L);
4329 else if (SCM_NUMBERP (x))
4330 return x;
4331 else
4332 SCM_WTA_DISPATCH_1 (g_product, x, SCM_ARG1, s_product);
4333 }
4334
4335 if (SCM_I_INUMP (x))
4336 {
4337 long xx;
4338
4339 intbig:
4340 xx = SCM_I_INUM (x);
4341
4342 switch (xx)
4343 {
4344 case 0: return x; break;
4345 case 1: return y; break;
4346 }
4347
4348 if (SCM_I_INUMP (y))
4349 {
4350 long yy = SCM_I_INUM (y);
4351 long kk = xx * yy;
4352 SCM k = SCM_I_MAKINUM (kk);
4353 if ((kk == SCM_I_INUM (k)) && (kk / xx == yy))
4354 return k;
4355 else
4356 {
4357 SCM result = scm_i_long2big (xx);
4358 mpz_mul_si (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result), yy);
4359 return scm_i_normbig (result);
4360 }
4361 }
4362 else if (SCM_BIGP (y))
4363 {
4364 SCM result = scm_i_mkbig ();
4365 mpz_mul_si (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (y), xx);
4366 scm_remember_upto_here_1 (y);
4367 return result;
4368 }
4369 else if (SCM_REALP (y))
4370 return scm_from_double (xx * SCM_REAL_VALUE (y));
4371 else if (SCM_COMPLEXP (y))
4372 return scm_c_make_rectangular (xx * SCM_COMPLEX_REAL (y),
4373 xx * SCM_COMPLEX_IMAG (y));
4374 else if (SCM_FRACTIONP (y))
4375 return scm_i_make_ratio (scm_product (x, SCM_FRACTION_NUMERATOR (y)),
4376 SCM_FRACTION_DENOMINATOR (y));
4377 else
4378 SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
4379 }
4380 else if (SCM_BIGP (x))
4381 {
4382 if (SCM_I_INUMP (y))
4383 {
4384 SCM_SWAP (x, y);
4385 goto intbig;
4386 }
4387 else if (SCM_BIGP (y))
4388 {
4389 SCM result = scm_i_mkbig ();
4390 mpz_mul (SCM_I_BIG_MPZ (result),
4391 SCM_I_BIG_MPZ (x),
4392 SCM_I_BIG_MPZ (y));
4393 scm_remember_upto_here_2 (x, y);
4394 return result;
4395 }
4396 else if (SCM_REALP (y))
4397 {
4398 double result = mpz_get_d (SCM_I_BIG_MPZ (x)) * SCM_REAL_VALUE (y);
4399 scm_remember_upto_here_1 (x);
4400 return scm_from_double (result);
4401 }
4402 else if (SCM_COMPLEXP (y))
4403 {
4404 double z = mpz_get_d (SCM_I_BIG_MPZ (x));
4405 scm_remember_upto_here_1 (x);
4406 return scm_c_make_rectangular (z * SCM_COMPLEX_REAL (y),
4407 z * SCM_COMPLEX_IMAG (y));
4408 }
4409 else if (SCM_FRACTIONP (y))
4410 return scm_i_make_ratio (scm_product (x, SCM_FRACTION_NUMERATOR (y)),
4411 SCM_FRACTION_DENOMINATOR (y));
4412 else
4413 SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
4414 }
4415 else if (SCM_REALP (x))
4416 {
4417 if (SCM_I_INUMP (y))
4418 return scm_from_double (SCM_I_INUM (y) * SCM_REAL_VALUE (x));
4419 else if (SCM_BIGP (y))
4420 {
4421 double result = mpz_get_d (SCM_I_BIG_MPZ (y)) * SCM_REAL_VALUE (x);
4422 scm_remember_upto_here_1 (y);
4423 return scm_from_double (result);
4424 }
4425 else if (SCM_REALP (y))
4426 return scm_from_double (SCM_REAL_VALUE (x) * SCM_REAL_VALUE (y));
4427 else if (SCM_COMPLEXP (y))
4428 return scm_c_make_rectangular (SCM_REAL_VALUE (x) * SCM_COMPLEX_REAL (y),
4429 SCM_REAL_VALUE (x) * SCM_COMPLEX_IMAG (y));
4430 else if (SCM_FRACTIONP (y))
4431 return scm_from_double (SCM_REAL_VALUE (x) * scm_i_fraction2double (y));
4432 else
4433 SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
4434 }
4435 else if (SCM_COMPLEXP (x))
4436 {
4437 if (SCM_I_INUMP (y))
4438 return scm_c_make_rectangular (SCM_I_INUM (y) * SCM_COMPLEX_REAL (x),
4439 SCM_I_INUM (y) * SCM_COMPLEX_IMAG (x));
4440 else if (SCM_BIGP (y))
4441 {
4442 double z = mpz_get_d (SCM_I_BIG_MPZ (y));
4443 scm_remember_upto_here_1 (y);
4444 return scm_c_make_rectangular (z * SCM_COMPLEX_REAL (x),
4445 z * SCM_COMPLEX_IMAG (x));
4446 }
4447 else if (SCM_REALP (y))
4448 return scm_c_make_rectangular (SCM_REAL_VALUE (y) * SCM_COMPLEX_REAL (x),
4449 SCM_REAL_VALUE (y) * SCM_COMPLEX_IMAG (x));
4450 else if (SCM_COMPLEXP (y))
4451 {
4452 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) * SCM_COMPLEX_REAL (y)
4453 - SCM_COMPLEX_IMAG (x) * SCM_COMPLEX_IMAG (y),
4454 SCM_COMPLEX_REAL (x) * SCM_COMPLEX_IMAG (y)
4455 + SCM_COMPLEX_IMAG (x) * SCM_COMPLEX_REAL (y));
4456 }
4457 else if (SCM_FRACTIONP (y))
4458 {
4459 double yy = scm_i_fraction2double (y);
4460 return scm_c_make_rectangular (yy * SCM_COMPLEX_REAL (x),
4461 yy * SCM_COMPLEX_IMAG (x));
4462 }
4463 else
4464 SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
4465 }
4466 else if (SCM_FRACTIONP (x))
4467 {
4468 if (SCM_I_INUMP (y))
4469 return scm_i_make_ratio (scm_product (y, SCM_FRACTION_NUMERATOR (x)),
4470 SCM_FRACTION_DENOMINATOR (x));
4471 else if (SCM_BIGP (y))
4472 return scm_i_make_ratio (scm_product (y, SCM_FRACTION_NUMERATOR (x)),
4473 SCM_FRACTION_DENOMINATOR (x));
4474 else if (SCM_REALP (y))
4475 return scm_from_double (scm_i_fraction2double (x) * SCM_REAL_VALUE (y));
4476 else if (SCM_COMPLEXP (y))
4477 {
4478 double xx = scm_i_fraction2double (x);
4479 return scm_c_make_rectangular (xx * SCM_COMPLEX_REAL (y),
4480 xx * SCM_COMPLEX_IMAG (y));
4481 }
4482 else if (SCM_FRACTIONP (y))
4483 /* a/b * c/d = ac / bd */
4484 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x),
4485 SCM_FRACTION_NUMERATOR (y)),
4486 scm_product (SCM_FRACTION_DENOMINATOR (x),
4487 SCM_FRACTION_DENOMINATOR (y)));
4488 else
4489 SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
4490 }
4491 else
4492 SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARG1, s_product);
4493 }
4494
4495 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4496 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4497 #define ALLOW_DIVIDE_BY_ZERO
4498 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4499 #endif
4500
4501 /* The code below for complex division is adapted from the GNU
4502 libstdc++, which adapted it from f2c's libF77, and is subject to
4503 this copyright: */
4504
4505 /****************************************************************
4506 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4507
4508 Permission to use, copy, modify, and distribute this software
4509 and its documentation for any purpose and without fee is hereby
4510 granted, provided that the above copyright notice appear in all
4511 copies and that both that the copyright notice and this
4512 permission notice and warranty disclaimer appear in supporting
4513 documentation, and that the names of AT&T Bell Laboratories or
4514 Bellcore or any of their entities not be used in advertising or
4515 publicity pertaining to distribution of the software without
4516 specific, written prior permission.
4517
4518 AT&T and Bellcore disclaim all warranties with regard to this
4519 software, including all implied warranties of merchantability
4520 and fitness. In no event shall AT&T or Bellcore be liable for
4521 any special, indirect or consequential damages or any damages
4522 whatsoever resulting from loss of use, data or profits, whether
4523 in an action of contract, negligence or other tortious action,
4524 arising out of or in connection with the use or performance of
4525 this software.
4526 ****************************************************************/
4527
4528 SCM_GPROC1 (s_divide, "/", scm_tc7_asubr, scm_divide, g_divide);
4529 /* Divide the first argument by the product of the remaining
4530 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4531 returned. */
4532 #define FUNC_NAME s_divide
4533 static SCM
4534 scm_i_divide (SCM x, SCM y, int inexact)
4535 {
4536 double a;
4537
4538 if (SCM_UNBNDP (y))
4539 {
4540 if (SCM_UNBNDP (x))
4541 SCM_WTA_DISPATCH_0 (g_divide, s_divide);
4542 else if (SCM_I_INUMP (x))
4543 {
4544 long xx = SCM_I_INUM (x);
4545 if (xx == 1 || xx == -1)
4546 return x;
4547 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4548 else if (xx == 0)
4549 scm_num_overflow (s_divide);
4550 #endif
4551 else
4552 {
4553 if (inexact)
4554 return scm_from_double (1.0 / (double) xx);
4555 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x);
4556 }
4557 }
4558 else if (SCM_BIGP (x))
4559 {
4560 if (inexact)
4561 return scm_from_double (1.0 / scm_i_big2dbl (x));
4562 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x);
4563 }
4564 else if (SCM_REALP (x))
4565 {
4566 double xx = SCM_REAL_VALUE (x);
4567 #ifndef ALLOW_DIVIDE_BY_ZERO
4568 if (xx == 0.0)
4569 scm_num_overflow (s_divide);
4570 else
4571 #endif
4572 return scm_from_double (1.0 / xx);
4573 }
4574 else if (SCM_COMPLEXP (x))
4575 {
4576 double r = SCM_COMPLEX_REAL (x);
4577 double i = SCM_COMPLEX_IMAG (x);
4578 if (r <= i)
4579 {
4580 double t = r / i;
4581 double d = i * (1.0 + t * t);
4582 return scm_c_make_rectangular (t / d, -1.0 / d);
4583 }
4584 else
4585 {
4586 double t = i / r;
4587 double d = r * (1.0 + t * t);
4588 return scm_c_make_rectangular (1.0 / d, -t / d);
4589 }
4590 }
4591 else if (SCM_FRACTIONP (x))
4592 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x),
4593 SCM_FRACTION_NUMERATOR (x));
4594 else
4595 SCM_WTA_DISPATCH_1 (g_divide, x, SCM_ARG1, s_divide);
4596 }
4597
4598 if (SCM_I_INUMP (x))
4599 {
4600 long xx = SCM_I_INUM (x);
4601 if (SCM_I_INUMP (y))
4602 {
4603 long yy = SCM_I_INUM (y);
4604 if (yy == 0)
4605 {
4606 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4607 scm_num_overflow (s_divide);
4608 #else
4609 return scm_from_double ((double) xx / (double) yy);
4610 #endif
4611 }
4612 else if (xx % yy != 0)
4613 {
4614 if (inexact)
4615 return scm_from_double ((double) xx / (double) yy);
4616 else return scm_i_make_ratio (x, y);
4617 }
4618 else
4619 {
4620 long z = xx / yy;
4621 if (SCM_FIXABLE (z))
4622 return SCM_I_MAKINUM (z);
4623 else
4624 return scm_i_long2big (z);
4625 }
4626 }
4627 else if (SCM_BIGP (y))
4628 {
4629 if (inexact)
4630 return scm_from_double ((double) xx / scm_i_big2dbl (y));
4631 else return scm_i_make_ratio (x, y);
4632 }
4633 else if (SCM_REALP (y))
4634 {
4635 double yy = SCM_REAL_VALUE (y);
4636 #ifndef ALLOW_DIVIDE_BY_ZERO
4637 if (yy == 0.0)
4638 scm_num_overflow (s_divide);
4639 else
4640 #endif
4641 return scm_from_double ((double) xx / yy);
4642 }
4643 else if (SCM_COMPLEXP (y))
4644 {
4645 a = xx;
4646 complex_div: /* y _must_ be a complex number */
4647 {
4648 double r = SCM_COMPLEX_REAL (y);
4649 double i = SCM_COMPLEX_IMAG (y);
4650 if (r <= i)
4651 {
4652 double t = r / i;
4653 double d = i * (1.0 + t * t);
4654 return scm_c_make_rectangular ((a * t) / d, -a / d);
4655 }
4656 else
4657 {
4658 double t = i / r;
4659 double d = r * (1.0 + t * t);
4660 return scm_c_make_rectangular (a / d, -(a * t) / d);
4661 }
4662 }
4663 }
4664 else if (SCM_FRACTIONP (y))
4665 /* a / b/c = ac / b */
4666 return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
4667 SCM_FRACTION_NUMERATOR (y));
4668 else
4669 SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
4670 }
4671 else if (SCM_BIGP (x))
4672 {
4673 if (SCM_I_INUMP (y))
4674 {
4675 long int yy = SCM_I_INUM (y);
4676 if (yy == 0)
4677 {
4678 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4679 scm_num_overflow (s_divide);
4680 #else
4681 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
4682 scm_remember_upto_here_1 (x);
4683 return (sgn == 0) ? scm_nan () : scm_inf ();
4684 #endif
4685 }
4686 else if (yy == 1)
4687 return x;
4688 else
4689 {
4690 /* FIXME: HMM, what are the relative performance issues here?
4691 We need to test. Is it faster on average to test
4692 divisible_p, then perform whichever operation, or is it
4693 faster to perform the integer div opportunistically and
4694 switch to real if there's a remainder? For now we take the
4695 middle ground: test, then if divisible, use the faster div
4696 func. */
4697
4698 long abs_yy = yy < 0 ? -yy : yy;
4699 int divisible_p = mpz_divisible_ui_p (SCM_I_BIG_MPZ (x), abs_yy);
4700
4701 if (divisible_p)
4702 {
4703 SCM result = scm_i_mkbig ();
4704 mpz_divexact_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), abs_yy);
4705 scm_remember_upto_here_1 (x);
4706 if (yy < 0)
4707 mpz_neg (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result));
4708 return scm_i_normbig (result);
4709 }
4710 else
4711 {
4712 if (inexact)
4713 return scm_from_double (scm_i_big2dbl (x) / (double) yy);
4714 else return scm_i_make_ratio (x, y);
4715 }
4716 }
4717 }
4718 else if (SCM_BIGP (y))
4719 {
4720 int y_is_zero = (mpz_sgn (SCM_I_BIG_MPZ (y)) == 0);
4721 if (y_is_zero)
4722 {
4723 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4724 scm_num_overflow (s_divide);
4725 #else
4726 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
4727 scm_remember_upto_here_1 (x);
4728 return (sgn == 0) ? scm_nan () : scm_inf ();
4729 #endif
4730 }
4731 else
4732 {
4733 /* big_x / big_y */
4734 int divisible_p = mpz_divisible_p (SCM_I_BIG_MPZ (x),
4735 SCM_I_BIG_MPZ (y));
4736 if (divisible_p)
4737 {
4738 SCM result = scm_i_mkbig ();
4739 mpz_divexact (SCM_I_BIG_MPZ (result),
4740 SCM_I_BIG_MPZ (x),
4741 SCM_I_BIG_MPZ (y));
4742 scm_remember_upto_here_2 (x, y);
4743 return scm_i_normbig (result);
4744 }
4745 else
4746 {
4747 if (inexact)
4748 {
4749 double dbx = mpz_get_d (SCM_I_BIG_MPZ (x));
4750 double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
4751 scm_remember_upto_here_2 (x, y);
4752 return scm_from_double (dbx / dby);
4753 }
4754 else return scm_i_make_ratio (x, y);
4755 }
4756 }
4757 }
4758 else if (SCM_REALP (y))
4759 {
4760 double yy = SCM_REAL_VALUE (y);
4761 #ifndef ALLOW_DIVIDE_BY_ZERO
4762 if (yy == 0.0)
4763 scm_num_overflow (s_divide);
4764 else
4765 #endif
4766 return scm_from_double (scm_i_big2dbl (x) / yy);
4767 }
4768 else if (SCM_COMPLEXP (y))
4769 {
4770 a = scm_i_big2dbl (x);
4771 goto complex_div;
4772 }
4773 else if (SCM_FRACTIONP (y))
4774 return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
4775 SCM_FRACTION_NUMERATOR (y));
4776 else
4777 SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
4778 }
4779 else if (SCM_REALP (x))
4780 {
4781 double rx = SCM_REAL_VALUE (x);
4782 if (SCM_I_INUMP (y))
4783 {
4784 long int yy = SCM_I_INUM (y);
4785 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4786 if (yy == 0)
4787 scm_num_overflow (s_divide);
4788 else
4789 #endif
4790 return scm_from_double (rx / (double) yy);
4791 }
4792 else if (SCM_BIGP (y))
4793 {
4794 double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
4795 scm_remember_upto_here_1 (y);
4796 return scm_from_double (rx / dby);
4797 }
4798 else if (SCM_REALP (y))
4799 {
4800 double yy = SCM_REAL_VALUE (y);
4801 #ifndef ALLOW_DIVIDE_BY_ZERO
4802 if (yy == 0.0)
4803 scm_num_overflow (s_divide);
4804 else
4805 #endif
4806 return scm_from_double (rx / yy);
4807 }
4808 else if (SCM_COMPLEXP (y))
4809 {
4810 a = rx;
4811 goto complex_div;
4812 }
4813 else if (SCM_FRACTIONP (y))
4814 return scm_from_double (rx / scm_i_fraction2double (y));
4815 else
4816 SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
4817 }
4818 else if (SCM_COMPLEXP (x))
4819 {
4820 double rx = SCM_COMPLEX_REAL (x);
4821 double ix = SCM_COMPLEX_IMAG (x);
4822 if (SCM_I_INUMP (y))
4823 {
4824 long int yy = SCM_I_INUM (y);
4825 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4826 if (yy == 0)
4827 scm_num_overflow (s_divide);
4828 else
4829 #endif
4830 {
4831 double d = yy;
4832 return scm_c_make_rectangular (rx / d, ix / d);
4833 }
4834 }
4835 else if (SCM_BIGP (y))
4836 {
4837 double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
4838 scm_remember_upto_here_1 (y);
4839 return scm_c_make_rectangular (rx / dby, ix / dby);
4840 }
4841 else if (SCM_REALP (y))
4842 {
4843 double yy = SCM_REAL_VALUE (y);
4844 #ifndef ALLOW_DIVIDE_BY_ZERO
4845 if (yy == 0.0)
4846 scm_num_overflow (s_divide);
4847 else
4848 #endif
4849 return scm_c_make_rectangular (rx / yy, ix / yy);
4850 }
4851 else if (SCM_COMPLEXP (y))
4852 {
4853 double ry = SCM_COMPLEX_REAL (y);
4854 double iy = SCM_COMPLEX_IMAG (y);
4855 if (ry <= iy)
4856 {
4857 double t = ry / iy;
4858 double d = iy * (1.0 + t * t);
4859 return scm_c_make_rectangular ((rx * t + ix) / d, (ix * t - rx) / d);
4860 }
4861 else
4862 {
4863 double t = iy / ry;
4864 double d = ry * (1.0 + t * t);
4865 return scm_c_make_rectangular ((rx + ix * t) / d, (ix - rx * t) / d);
4866 }
4867 }
4868 else if (SCM_FRACTIONP (y))
4869 {
4870 double yy = scm_i_fraction2double (y);
4871 return scm_c_make_rectangular (rx / yy, ix / yy);
4872 }
4873 else
4874 SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
4875 }
4876 else if (SCM_FRACTIONP (x))
4877 {
4878 if (SCM_I_INUMP (y))
4879 {
4880 long int yy = SCM_I_INUM (y);
4881 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4882 if (yy == 0)
4883 scm_num_overflow (s_divide);
4884 else
4885 #endif
4886 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x),
4887 scm_product (SCM_FRACTION_DENOMINATOR (x), y));
4888 }
4889 else if (SCM_BIGP (y))
4890 {
4891 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x),
4892 scm_product (SCM_FRACTION_DENOMINATOR (x), y));
4893 }
4894 else if (SCM_REALP (y))
4895 {
4896 double yy = SCM_REAL_VALUE (y);
4897 #ifndef ALLOW_DIVIDE_BY_ZERO
4898 if (yy == 0.0)
4899 scm_num_overflow (s_divide);
4900 else
4901 #endif
4902 return scm_from_double (scm_i_fraction2double (x) / yy);
4903 }
4904 else if (SCM_COMPLEXP (y))
4905 {
4906 a = scm_i_fraction2double (x);
4907 goto complex_div;
4908 }
4909 else if (SCM_FRACTIONP (y))
4910 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x), SCM_FRACTION_DENOMINATOR (y)),
4911 scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x)));
4912 else
4913 SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
4914 }
4915 else
4916 SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARG1, s_divide);
4917 }
4918
4919 SCM
4920 scm_divide (SCM x, SCM y)
4921 {
4922 return scm_i_divide (x, y, 0);
4923 }
4924
4925 static SCM scm_divide2real (SCM x, SCM y)
4926 {
4927 return scm_i_divide (x, y, 1);
4928 }
4929 #undef FUNC_NAME
4930
4931
4932 double
4933 scm_asinh (double x)
4934 {
4935 #if HAVE_ASINH
4936 return asinh (x);
4937 #else
4938 #define asinh scm_asinh
4939 return log (x + sqrt (x * x + 1));
4940 #endif
4941 }
4942 SCM_GPROC1 (s_asinh, "$asinh", scm_tc7_dsubr, (SCM (*)()) asinh, g_asinh);
4943 /* "Return the inverse hyperbolic sine of @var{x}."
4944 */
4945
4946
4947 double
4948 scm_acosh (double x)
4949 {
4950 #if HAVE_ACOSH
4951 return acosh (x);
4952 #else
4953 #define acosh scm_acosh
4954 return log (x + sqrt (x * x - 1));
4955 #endif
4956 }
4957 SCM_GPROC1 (s_acosh, "$acosh", scm_tc7_dsubr, (SCM (*)()) acosh, g_acosh);
4958 /* "Return the inverse hyperbolic cosine of @var{x}."
4959 */
4960
4961
4962 double
4963 scm_atanh (double x)
4964 {
4965 #if HAVE_ATANH
4966 return atanh (x);
4967 #else
4968 #define atanh scm_atanh
4969 return 0.5 * log ((1 + x) / (1 - x));
4970 #endif
4971 }
4972 SCM_GPROC1 (s_atanh, "$atanh", scm_tc7_dsubr, (SCM (*)()) atanh, g_atanh);
4973 /* "Return the inverse hyperbolic tangent of @var{x}."
4974 */
4975
4976
4977 double
4978 scm_c_truncate (double x)
4979 {
4980 #if HAVE_TRUNC
4981 return trunc (x);
4982 #else
4983 if (x < 0.0)
4984 return -floor (-x);
4985 return floor (x);
4986 #endif
4987 }
4988
4989 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
4990 half-way case (ie. when x is an integer plus 0.5) going upwards.
4991 Then half-way cases are identified and adjusted down if the
4992 round-upwards didn't give the desired even integer.
4993
4994 "plus_half == result" identifies a half-way case. If plus_half, which is
4995 x + 0.5, is an integer then x must be an integer plus 0.5.
4996
4997 An odd "result" value is identified with result/2 != floor(result/2).
4998 This is done with plus_half, since that value is ready for use sooner in
4999 a pipelined cpu, and we're already requiring plus_half == result.
5000
5001 Note however that we need to be careful when x is big and already an
5002 integer. In that case "x+0.5" may round to an adjacent integer, causing
5003 us to return such a value, incorrectly. For instance if the hardware is
5004 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5005 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5006 returned. Or if the hardware is in round-upwards mode, then other bigger
5007 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5008 representable value, 2^128+2^76 (or whatever), again incorrect.
5009
5010 These bad roundings of x+0.5 are avoided by testing at the start whether
5011 x is already an integer. If it is then clearly that's the desired result
5012 already. And if it's not then the exponent must be small enough to allow
5013 an 0.5 to be represented, and hence added without a bad rounding. */
5014
5015 double
5016 scm_c_round (double x)
5017 {
5018 double plus_half, result;
5019
5020 if (x == floor (x))
5021 return x;
5022
5023 plus_half = x + 0.5;
5024 result = floor (plus_half);
5025 /* Adjust so that the rounding is towards even. */
5026 return ((plus_half == result && plus_half / 2 != floor (plus_half / 2))
5027 ? result - 1
5028 : result);
5029 }
5030
5031 SCM_DEFINE (scm_truncate_number, "truncate", 1, 0, 0,
5032 (SCM x),
5033 "Round the number @var{x} towards zero.")
5034 #define FUNC_NAME s_scm_truncate_number
5035 {
5036 if (scm_is_false (scm_negative_p (x)))
5037 return scm_floor (x);
5038 else
5039 return scm_ceiling (x);
5040 }
5041 #undef FUNC_NAME
5042
5043 static SCM exactly_one_half;
5044
5045 SCM_DEFINE (scm_round_number, "round", 1, 0, 0,
5046 (SCM x),
5047 "Round the number @var{x} towards the nearest integer. "
5048 "When it is exactly halfway between two integers, "
5049 "round towards the even one.")
5050 #define FUNC_NAME s_scm_round_number
5051 {
5052 if (SCM_I_INUMP (x) || SCM_BIGP (x))
5053 return x;
5054 else if (SCM_REALP (x))
5055 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x)));
5056 else
5057 {
5058 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5059 single quotient+remainder division then examining to see which way
5060 the rounding should go. */
5061 SCM plus_half = scm_sum (x, exactly_one_half);
5062 SCM result = scm_floor (plus_half);
5063 /* Adjust so that the rounding is towards even. */
5064 if (scm_is_true (scm_num_eq_p (plus_half, result))
5065 && scm_is_true (scm_odd_p (result)))
5066 return scm_difference (result, SCM_I_MAKINUM (1));
5067 else
5068 return result;
5069 }
5070 }
5071 #undef FUNC_NAME
5072
5073 SCM_PRIMITIVE_GENERIC (scm_floor, "floor", 1, 0, 0,
5074 (SCM x),
5075 "Round the number @var{x} towards minus infinity.")
5076 #define FUNC_NAME s_scm_floor
5077 {
5078 if (SCM_I_INUMP (x) || SCM_BIGP (x))
5079 return x;
5080 else if (SCM_REALP (x))
5081 return scm_from_double (floor (SCM_REAL_VALUE (x)));
5082 else if (SCM_FRACTIONP (x))
5083 {
5084 SCM q = scm_quotient (SCM_FRACTION_NUMERATOR (x),
5085 SCM_FRACTION_DENOMINATOR (x));
5086 if (scm_is_false (scm_negative_p (x)))
5087 {
5088 /* For positive x, rounding towards zero is correct. */
5089 return q;
5090 }
5091 else
5092 {
5093 /* For negative x, we need to return q-1 unless x is an
5094 integer. But fractions are never integer, per our
5095 assumptions. */
5096 return scm_difference (q, SCM_I_MAKINUM (1));
5097 }
5098 }
5099 else
5100 SCM_WTA_DISPATCH_1 (g_scm_floor, x, 1, s_scm_floor);
5101 }
5102 #undef FUNC_NAME
5103
5104 SCM_PRIMITIVE_GENERIC (scm_ceiling, "ceiling", 1, 0, 0,
5105 (SCM x),
5106 "Round the number @var{x} towards infinity.")
5107 #define FUNC_NAME s_scm_ceiling
5108 {
5109 if (SCM_I_INUMP (x) || SCM_BIGP (x))
5110 return x;
5111 else if (SCM_REALP (x))
5112 return scm_from_double (ceil (SCM_REAL_VALUE (x)));
5113 else if (SCM_FRACTIONP (x))
5114 {
5115 SCM q = scm_quotient (SCM_FRACTION_NUMERATOR (x),
5116 SCM_FRACTION_DENOMINATOR (x));
5117 if (scm_is_false (scm_positive_p (x)))
5118 {
5119 /* For negative x, rounding towards zero is correct. */
5120 return q;
5121 }
5122 else
5123 {
5124 /* For positive x, we need to return q+1 unless x is an
5125 integer. But fractions are never integer, per our
5126 assumptions. */
5127 return scm_sum (q, SCM_I_MAKINUM (1));
5128 }
5129 }
5130 else
5131 SCM_WTA_DISPATCH_1 (g_scm_ceiling, x, 1, s_scm_ceiling);
5132 }
5133 #undef FUNC_NAME
5134
5135 SCM_GPROC1 (s_i_sqrt, "$sqrt", scm_tc7_dsubr, (SCM (*)()) sqrt, g_i_sqrt);
5136 /* "Return the square root of the real number @var{x}."
5137 */
5138 SCM_GPROC1 (s_i_abs, "$abs", scm_tc7_dsubr, (SCM (*)()) fabs, g_i_abs);
5139 /* "Return the absolute value of the real number @var{x}."
5140 */
5141 SCM_GPROC1 (s_i_exp, "$exp", scm_tc7_dsubr, (SCM (*)()) exp, g_i_exp);
5142 /* "Return the @var{x}th power of e."
5143 */
5144 SCM_GPROC1 (s_i_log, "$log", scm_tc7_dsubr, (SCM (*)()) log, g_i_log);
5145 /* "Return the natural logarithm of the real number @var{x}."
5146 */
5147 SCM_GPROC1 (s_i_sin, "$sin", scm_tc7_dsubr, (SCM (*)()) sin, g_i_sin);
5148 /* "Return the sine of the real number @var{x}."
5149 */
5150 SCM_GPROC1 (s_i_cos, "$cos", scm_tc7_dsubr, (SCM (*)()) cos, g_i_cos);
5151 /* "Return the cosine of the real number @var{x}."
5152 */
5153 SCM_GPROC1 (s_i_tan, "$tan", scm_tc7_dsubr, (SCM (*)()) tan, g_i_tan);
5154 /* "Return the tangent of the real number @var{x}."
5155 */
5156 SCM_GPROC1 (s_i_asin, "$asin", scm_tc7_dsubr, (SCM (*)()) asin, g_i_asin);
5157 /* "Return the arc sine of the real number @var{x}."
5158 */
5159 SCM_GPROC1 (s_i_acos, "$acos", scm_tc7_dsubr, (SCM (*)()) acos, g_i_acos);
5160 /* "Return the arc cosine of the real number @var{x}."
5161 */
5162 SCM_GPROC1 (s_i_atan, "$atan", scm_tc7_dsubr, (SCM (*)()) atan, g_i_atan);
5163 /* "Return the arc tangent of the real number @var{x}."
5164 */
5165 SCM_GPROC1 (s_i_sinh, "$sinh", scm_tc7_dsubr, (SCM (*)()) sinh, g_i_sinh);
5166 /* "Return the hyperbolic sine of the real number @var{x}."
5167 */
5168 SCM_GPROC1 (s_i_cosh, "$cosh", scm_tc7_dsubr, (SCM (*)()) cosh, g_i_cosh);
5169 /* "Return the hyperbolic cosine of the real number @var{x}."
5170 */
5171 SCM_GPROC1 (s_i_tanh, "$tanh", scm_tc7_dsubr, (SCM (*)()) tanh, g_i_tanh);
5172 /* "Return the hyperbolic tangent of the real number @var{x}."
5173 */
5174
5175 struct dpair
5176 {
5177 double x, y;
5178 };
5179
5180 static void scm_two_doubles (SCM x,
5181 SCM y,
5182 const char *sstring,
5183 struct dpair * xy);
5184
5185 static void
5186 scm_two_doubles (SCM x, SCM y, const char *sstring, struct dpair *xy)
5187 {
5188 if (SCM_I_INUMP (x))
5189 xy->x = SCM_I_INUM (x);
5190 else if (SCM_BIGP (x))
5191 xy->x = scm_i_big2dbl (x);
5192 else if (SCM_REALP (x))
5193 xy->x = SCM_REAL_VALUE (x);
5194 else if (SCM_FRACTIONP (x))
5195 xy->x = scm_i_fraction2double (x);
5196 else
5197 scm_wrong_type_arg (sstring, SCM_ARG1, x);
5198
5199 if (SCM_I_INUMP (y))
5200 xy->y = SCM_I_INUM (y);
5201 else if (SCM_BIGP (y))
5202 xy->y = scm_i_big2dbl (y);
5203 else if (SCM_REALP (y))
5204 xy->y = SCM_REAL_VALUE (y);
5205 else if (SCM_FRACTIONP (y))
5206 xy->y = scm_i_fraction2double (y);
5207 else
5208 scm_wrong_type_arg (sstring, SCM_ARG2, y);
5209 }
5210
5211
5212 SCM_DEFINE (scm_sys_expt, "$expt", 2, 0, 0,
5213 (SCM x, SCM y),
5214 "Return @var{x} raised to the power of @var{y}. This\n"
5215 "procedure does not accept complex arguments.")
5216 #define FUNC_NAME s_scm_sys_expt
5217 {
5218 struct dpair xy;
5219 scm_two_doubles (x, y, FUNC_NAME, &xy);
5220 return scm_from_double (pow (xy.x, xy.y));
5221 }
5222 #undef FUNC_NAME
5223
5224
5225 SCM_DEFINE (scm_sys_atan2, "$atan2", 2, 0, 0,
5226 (SCM x, SCM y),
5227 "Return the arc tangent of the two arguments @var{x} and\n"
5228 "@var{y}. This is similar to calculating the arc tangent of\n"
5229 "@var{x} / @var{y}, except that the signs of both arguments\n"
5230 "are used to determine the quadrant of the result. This\n"
5231 "procedure does not accept complex arguments.")
5232 #define FUNC_NAME s_scm_sys_atan2
5233 {
5234 struct dpair xy;
5235 scm_two_doubles (x, y, FUNC_NAME, &xy);
5236 return scm_from_double (atan2 (xy.x, xy.y));
5237 }
5238 #undef FUNC_NAME
5239
5240 SCM
5241 scm_c_make_rectangular (double re, double im)
5242 {
5243 if (im == 0.0)
5244 return scm_from_double (re);
5245 else
5246 {
5247 SCM z;
5248 SCM_NEWSMOB (z, scm_tc16_complex, scm_gc_malloc (sizeof (scm_t_complex),
5249 "complex"));
5250 SCM_COMPLEX_REAL (z) = re;
5251 SCM_COMPLEX_IMAG (z) = im;
5252 return z;
5253 }
5254 }
5255
5256 SCM_DEFINE (scm_make_rectangular, "make-rectangular", 2, 0, 0,
5257 (SCM real, SCM imaginary),
5258 "Return a complex number constructed of the given @var{real} and\n"
5259 "@var{imaginary} parts.")
5260 #define FUNC_NAME s_scm_make_rectangular
5261 {
5262 struct dpair xy;
5263 scm_two_doubles (real, imaginary, FUNC_NAME, &xy);
5264 return scm_c_make_rectangular (xy.x, xy.y);
5265 }
5266 #undef FUNC_NAME
5267
5268 SCM
5269 scm_c_make_polar (double mag, double ang)
5270 {
5271 double s, c;
5272 #if HAVE_SINCOS
5273 sincos (ang, &s, &c);
5274 #else
5275 s = sin (ang);
5276 c = cos (ang);
5277 #endif
5278 return scm_c_make_rectangular (mag * c, mag * s);
5279 }
5280
5281 SCM_DEFINE (scm_make_polar, "make-polar", 2, 0, 0,
5282 (SCM x, SCM y),
5283 "Return the complex number @var{x} * e^(i * @var{y}).")
5284 #define FUNC_NAME s_scm_make_polar
5285 {
5286 struct dpair xy;
5287 scm_two_doubles (x, y, FUNC_NAME, &xy);
5288 return scm_c_make_polar (xy.x, xy.y);
5289 }
5290 #undef FUNC_NAME
5291
5292
5293 SCM_GPROC (s_real_part, "real-part", 1, 0, 0, scm_real_part, g_real_part);
5294 /* "Return the real part of the number @var{z}."
5295 */
5296 SCM
5297 scm_real_part (SCM z)
5298 {
5299 if (SCM_I_INUMP (z))
5300 return z;
5301 else if (SCM_BIGP (z))
5302 return z;
5303 else if (SCM_REALP (z))
5304 return z;
5305 else if (SCM_COMPLEXP (z))
5306 return scm_from_double (SCM_COMPLEX_REAL (z));
5307 else if (SCM_FRACTIONP (z))
5308 return z;
5309 else
5310 SCM_WTA_DISPATCH_1 (g_real_part, z, SCM_ARG1, s_real_part);
5311 }
5312
5313
5314 SCM_GPROC (s_imag_part, "imag-part", 1, 0, 0, scm_imag_part, g_imag_part);
5315 /* "Return the imaginary part of the number @var{z}."
5316 */
5317 SCM
5318 scm_imag_part (SCM z)
5319 {
5320 if (SCM_I_INUMP (z))
5321 return SCM_INUM0;
5322 else if (SCM_BIGP (z))
5323 return SCM_INUM0;
5324 else if (SCM_REALP (z))
5325 return scm_flo0;
5326 else if (SCM_COMPLEXP (z))
5327 return scm_from_double (SCM_COMPLEX_IMAG (z));
5328 else if (SCM_FRACTIONP (z))
5329 return SCM_INUM0;
5330 else
5331 SCM_WTA_DISPATCH_1 (g_imag_part, z, SCM_ARG1, s_imag_part);
5332 }
5333
5334 SCM_GPROC (s_numerator, "numerator", 1, 0, 0, scm_numerator, g_numerator);
5335 /* "Return the numerator of the number @var{z}."
5336 */
5337 SCM
5338 scm_numerator (SCM z)
5339 {
5340 if (SCM_I_INUMP (z))
5341 return z;
5342 else if (SCM_BIGP (z))
5343 return z;
5344 else if (SCM_FRACTIONP (z))
5345 {
5346 scm_i_fraction_reduce (z);
5347 return SCM_FRACTION_NUMERATOR (z);
5348 }
5349 else if (SCM_REALP (z))
5350 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z)));
5351 else
5352 SCM_WTA_DISPATCH_1 (g_numerator, z, SCM_ARG1, s_numerator);
5353 }
5354
5355
5356 SCM_GPROC (s_denominator, "denominator", 1, 0, 0, scm_denominator, g_denominator);
5357 /* "Return the denominator of the number @var{z}."
5358 */
5359 SCM
5360 scm_denominator (SCM z)
5361 {
5362 if (SCM_I_INUMP (z))
5363 return SCM_I_MAKINUM (1);
5364 else if (SCM_BIGP (z))
5365 return SCM_I_MAKINUM (1);
5366 else if (SCM_FRACTIONP (z))
5367 {
5368 scm_i_fraction_reduce (z);
5369 return SCM_FRACTION_DENOMINATOR (z);
5370 }
5371 else if (SCM_REALP (z))
5372 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z)));
5373 else
5374 SCM_WTA_DISPATCH_1 (g_denominator, z, SCM_ARG1, s_denominator);
5375 }
5376
5377 SCM_GPROC (s_magnitude, "magnitude", 1, 0, 0, scm_magnitude, g_magnitude);
5378 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5379 * "@code{abs} for real arguments, but also allows complex numbers."
5380 */
5381 SCM
5382 scm_magnitude (SCM z)
5383 {
5384 if (SCM_I_INUMP (z))
5385 {
5386 long int zz = SCM_I_INUM (z);
5387 if (zz >= 0)
5388 return z;
5389 else if (SCM_POSFIXABLE (-zz))
5390 return SCM_I_MAKINUM (-zz);
5391 else
5392 return scm_i_long2big (-zz);
5393 }
5394 else if (SCM_BIGP (z))
5395 {
5396 int sgn = mpz_sgn (SCM_I_BIG_MPZ (z));
5397 scm_remember_upto_here_1 (z);
5398 if (sgn < 0)
5399 return scm_i_clonebig (z, 0);
5400 else
5401 return z;
5402 }
5403 else if (SCM_REALP (z))
5404 return scm_from_double (fabs (SCM_REAL_VALUE (z)));
5405 else if (SCM_COMPLEXP (z))
5406 return scm_from_double (hypot (SCM_COMPLEX_REAL (z), SCM_COMPLEX_IMAG (z)));
5407 else if (SCM_FRACTIONP (z))
5408 {
5409 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z))))
5410 return z;
5411 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z), SCM_UNDEFINED),
5412 SCM_FRACTION_DENOMINATOR (z));
5413 }
5414 else
5415 SCM_WTA_DISPATCH_1 (g_magnitude, z, SCM_ARG1, s_magnitude);
5416 }
5417
5418
5419 SCM_GPROC (s_angle, "angle", 1, 0, 0, scm_angle, g_angle);
5420 /* "Return the angle of the complex number @var{z}."
5421 */
5422 SCM
5423 scm_angle (SCM z)
5424 {
5425 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5426 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5427 But if atan2 follows the floating point rounding mode, then the value
5428 is not a constant. Maybe it'd be close enough though. */
5429 if (SCM_I_INUMP (z))
5430 {
5431 if (SCM_I_INUM (z) >= 0)
5432 return scm_flo0;
5433 else
5434 return scm_from_double (atan2 (0.0, -1.0));
5435 }
5436 else if (SCM_BIGP (z))
5437 {
5438 int sgn = mpz_sgn (SCM_I_BIG_MPZ (z));
5439 scm_remember_upto_here_1 (z);
5440 if (sgn < 0)
5441 return scm_from_double (atan2 (0.0, -1.0));
5442 else
5443 return scm_flo0;
5444 }
5445 else if (SCM_REALP (z))
5446 {
5447 if (SCM_REAL_VALUE (z) >= 0)
5448 return scm_flo0;
5449 else
5450 return scm_from_double (atan2 (0.0, -1.0));
5451 }
5452 else if (SCM_COMPLEXP (z))
5453 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z), SCM_COMPLEX_REAL (z)));
5454 else if (SCM_FRACTIONP (z))
5455 {
5456 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z))))
5457 return scm_flo0;
5458 else return scm_from_double (atan2 (0.0, -1.0));
5459 }
5460 else
5461 SCM_WTA_DISPATCH_1 (g_angle, z, SCM_ARG1, s_angle);
5462 }
5463
5464
5465 SCM_GPROC (s_exact_to_inexact, "exact->inexact", 1, 0, 0, scm_exact_to_inexact, g_exact_to_inexact);
5466 /* Convert the number @var{x} to its inexact representation.\n"
5467 */
5468 SCM
5469 scm_exact_to_inexact (SCM z)
5470 {
5471 if (SCM_I_INUMP (z))
5472 return scm_from_double ((double) SCM_I_INUM (z));
5473 else if (SCM_BIGP (z))
5474 return scm_from_double (scm_i_big2dbl (z));
5475 else if (SCM_FRACTIONP (z))
5476 return scm_from_double (scm_i_fraction2double (z));
5477 else if (SCM_INEXACTP (z))
5478 return z;
5479 else
5480 SCM_WTA_DISPATCH_1 (g_exact_to_inexact, z, 1, s_exact_to_inexact);
5481 }
5482
5483
5484 SCM_DEFINE (scm_inexact_to_exact, "inexact->exact", 1, 0, 0,
5485 (SCM z),
5486 "Return an exact number that is numerically closest to @var{z}.")
5487 #define FUNC_NAME s_scm_inexact_to_exact
5488 {
5489 if (SCM_I_INUMP (z))
5490 return z;
5491 else if (SCM_BIGP (z))
5492 return z;
5493 else if (SCM_REALP (z))
5494 {
5495 if (xisinf (SCM_REAL_VALUE (z)) || xisnan (SCM_REAL_VALUE (z)))
5496 SCM_OUT_OF_RANGE (1, z);
5497 else
5498 {
5499 mpq_t frac;
5500 SCM q;
5501
5502 mpq_init (frac);
5503 mpq_set_d (frac, SCM_REAL_VALUE (z));
5504 q = scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac)),
5505 scm_i_mpz2num (mpq_denref (frac)));
5506
5507 /* When scm_i_make_ratio throws, we leak the memory allocated
5508 for frac...
5509 */
5510 mpq_clear (frac);
5511 return q;
5512 }
5513 }
5514 else if (SCM_FRACTIONP (z))
5515 return z;
5516 else
5517 SCM_WRONG_TYPE_ARG (1, z);
5518 }
5519 #undef FUNC_NAME
5520
5521 SCM_DEFINE (scm_rationalize, "rationalize", 2, 0, 0,
5522 (SCM x, SCM err),
5523 "Return an exact number that is within @var{err} of @var{x}.")
5524 #define FUNC_NAME s_scm_rationalize
5525 {
5526 if (SCM_I_INUMP (x))
5527 return x;
5528 else if (SCM_BIGP (x))
5529 return x;
5530 else if ((SCM_REALP (x)) || SCM_FRACTIONP (x))
5531 {
5532 /* Use continued fractions to find closest ratio. All
5533 arithmetic is done with exact numbers.
5534 */
5535
5536 SCM ex = scm_inexact_to_exact (x);
5537 SCM int_part = scm_floor (ex);
5538 SCM tt = SCM_I_MAKINUM (1);
5539 SCM a1 = SCM_I_MAKINUM (0), a2 = SCM_I_MAKINUM (1), a = SCM_I_MAKINUM (0);
5540 SCM b1 = SCM_I_MAKINUM (1), b2 = SCM_I_MAKINUM (0), b = SCM_I_MAKINUM (0);
5541 SCM rx;
5542 int i = 0;
5543
5544 if (scm_is_true (scm_num_eq_p (ex, int_part)))
5545 return ex;
5546
5547 ex = scm_difference (ex, int_part); /* x = x-int_part */
5548 rx = scm_divide (ex, SCM_UNDEFINED); /* rx = 1/x */
5549
5550 /* We stop after a million iterations just to be absolutely sure
5551 that we don't go into an infinite loop. The process normally
5552 converges after less than a dozen iterations.
5553 */
5554
5555 err = scm_abs (err);
5556 while (++i < 1000000)
5557 {
5558 a = scm_sum (scm_product (a1, tt), a2); /* a = a1*tt + a2 */
5559 b = scm_sum (scm_product (b1, tt), b2); /* b = b1*tt + b2 */
5560 if (scm_is_false (scm_zero_p (b)) && /* b != 0 */
5561 scm_is_false
5562 (scm_gr_p (scm_abs (scm_difference (ex, scm_divide (a, b))),
5563 err))) /* abs(x-a/b) <= err */
5564 {
5565 SCM res = scm_sum (int_part, scm_divide (a, b));
5566 if (scm_is_false (scm_exact_p (x))
5567 || scm_is_false (scm_exact_p (err)))
5568 return scm_exact_to_inexact (res);
5569 else
5570 return res;
5571 }
5572 rx = scm_divide (scm_difference (rx, tt), /* rx = 1/(rx - tt) */
5573 SCM_UNDEFINED);
5574 tt = scm_floor (rx); /* tt = floor (rx) */
5575 a2 = a1;
5576 b2 = b1;
5577 a1 = a;
5578 b1 = b;
5579 }
5580 scm_num_overflow (s_scm_rationalize);
5581 }
5582 else
5583 SCM_WRONG_TYPE_ARG (1, x);
5584 }
5585 #undef FUNC_NAME
5586
5587 /* conversion functions */
5588
5589 int
5590 scm_is_integer (SCM val)
5591 {
5592 return scm_is_true (scm_integer_p (val));
5593 }
5594
5595 int
5596 scm_is_signed_integer (SCM val, scm_t_intmax min, scm_t_intmax max)
5597 {
5598 if (SCM_I_INUMP (val))
5599 {
5600 scm_t_signed_bits n = SCM_I_INUM (val);
5601 return n >= min && n <= max;
5602 }
5603 else if (SCM_BIGP (val))
5604 {
5605 if (min >= SCM_MOST_NEGATIVE_FIXNUM && max <= SCM_MOST_POSITIVE_FIXNUM)
5606 return 0;
5607 else if (min >= LONG_MIN && max <= LONG_MAX)
5608 {
5609 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val)))
5610 {
5611 long n = mpz_get_si (SCM_I_BIG_MPZ (val));
5612 return n >= min && n <= max;
5613 }
5614 else
5615 return 0;
5616 }
5617 else
5618 {
5619 scm_t_intmax n;
5620 size_t count;
5621
5622 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val), 2)
5623 > CHAR_BIT*sizeof (scm_t_uintmax))
5624 return 0;
5625
5626 mpz_export (&n, &count, 1, sizeof (scm_t_uintmax), 0, 0,
5627 SCM_I_BIG_MPZ (val));
5628
5629 if (mpz_sgn (SCM_I_BIG_MPZ (val)) >= 0)
5630 {
5631 if (n < 0)
5632 return 0;
5633 }
5634 else
5635 {
5636 n = -n;
5637 if (n >= 0)
5638 return 0;
5639 }
5640
5641 return n >= min && n <= max;
5642 }
5643 }
5644 else
5645 return 0;
5646 }
5647
5648 int
5649 scm_is_unsigned_integer (SCM val, scm_t_uintmax min, scm_t_uintmax max)
5650 {
5651 if (SCM_I_INUMP (val))
5652 {
5653 scm_t_signed_bits n = SCM_I_INUM (val);
5654 return n >= 0 && ((scm_t_uintmax)n) >= min && ((scm_t_uintmax)n) <= max;
5655 }
5656 else if (SCM_BIGP (val))
5657 {
5658 if (max <= SCM_MOST_POSITIVE_FIXNUM)
5659 return 0;
5660 else if (max <= ULONG_MAX)
5661 {
5662 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val)))
5663 {
5664 unsigned long n = mpz_get_ui (SCM_I_BIG_MPZ (val));
5665 return n >= min && n <= max;
5666 }
5667 else
5668 return 0;
5669 }
5670 else
5671 {
5672 scm_t_uintmax n;
5673 size_t count;
5674
5675 if (mpz_sgn (SCM_I_BIG_MPZ (val)) < 0)
5676 return 0;
5677
5678 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val), 2)
5679 > CHAR_BIT*sizeof (scm_t_uintmax))
5680 return 0;
5681
5682 mpz_export (&n, &count, 1, sizeof (scm_t_uintmax), 0, 0,
5683 SCM_I_BIG_MPZ (val));
5684
5685 return n >= min && n <= max;
5686 }
5687 }
5688 else
5689 return 0;
5690 }
5691
5692 static void
5693 scm_i_range_error (SCM bad_val, SCM min, SCM max)
5694 {
5695 scm_error (scm_out_of_range_key,
5696 NULL,
5697 "Value out of range ~S to ~S: ~S",
5698 scm_list_3 (min, max, bad_val),
5699 scm_list_1 (bad_val));
5700 }
5701
5702 #define TYPE scm_t_intmax
5703 #define TYPE_MIN min
5704 #define TYPE_MAX max
5705 #define SIZEOF_TYPE 0
5706 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5707 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5708 #include "libguile/conv-integer.i.c"
5709
5710 #define TYPE scm_t_uintmax
5711 #define TYPE_MIN min
5712 #define TYPE_MAX max
5713 #define SIZEOF_TYPE 0
5714 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5715 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5716 #include "libguile/conv-uinteger.i.c"
5717
5718 #define TYPE scm_t_int8
5719 #define TYPE_MIN SCM_T_INT8_MIN
5720 #define TYPE_MAX SCM_T_INT8_MAX
5721 #define SIZEOF_TYPE 1
5722 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5723 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5724 #include "libguile/conv-integer.i.c"
5725
5726 #define TYPE scm_t_uint8
5727 #define TYPE_MIN 0
5728 #define TYPE_MAX SCM_T_UINT8_MAX
5729 #define SIZEOF_TYPE 1
5730 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5731 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5732 #include "libguile/conv-uinteger.i.c"
5733
5734 #define TYPE scm_t_int16
5735 #define TYPE_MIN SCM_T_INT16_MIN
5736 #define TYPE_MAX SCM_T_INT16_MAX
5737 #define SIZEOF_TYPE 2
5738 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5739 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5740 #include "libguile/conv-integer.i.c"
5741
5742 #define TYPE scm_t_uint16
5743 #define TYPE_MIN 0
5744 #define TYPE_MAX SCM_T_UINT16_MAX
5745 #define SIZEOF_TYPE 2
5746 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5747 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5748 #include "libguile/conv-uinteger.i.c"
5749
5750 #define TYPE scm_t_int32
5751 #define TYPE_MIN SCM_T_INT32_MIN
5752 #define TYPE_MAX SCM_T_INT32_MAX
5753 #define SIZEOF_TYPE 4
5754 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5755 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5756 #include "libguile/conv-integer.i.c"
5757
5758 #define TYPE scm_t_uint32
5759 #define TYPE_MIN 0
5760 #define TYPE_MAX SCM_T_UINT32_MAX
5761 #define SIZEOF_TYPE 4
5762 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5763 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5764 #include "libguile/conv-uinteger.i.c"
5765
5766 #if SCM_HAVE_T_INT64
5767
5768 #define TYPE scm_t_int64
5769 #define TYPE_MIN SCM_T_INT64_MIN
5770 #define TYPE_MAX SCM_T_INT64_MAX
5771 #define SIZEOF_TYPE 8
5772 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5773 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5774 #include "libguile/conv-integer.i.c"
5775
5776 #define TYPE scm_t_uint64
5777 #define TYPE_MIN 0
5778 #define TYPE_MAX SCM_T_UINT64_MAX
5779 #define SIZEOF_TYPE 8
5780 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5781 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5782 #include "libguile/conv-uinteger.i.c"
5783
5784 #endif
5785
5786 void
5787 scm_to_mpz (SCM val, mpz_t rop)
5788 {
5789 if (SCM_I_INUMP (val))
5790 mpz_set_si (rop, SCM_I_INUM (val));
5791 else if (SCM_BIGP (val))
5792 mpz_set (rop, SCM_I_BIG_MPZ (val));
5793 else
5794 scm_wrong_type_arg_msg (NULL, 0, val, "exact integer");
5795 }
5796
5797 SCM
5798 scm_from_mpz (mpz_t val)
5799 {
5800 return scm_i_mpz2num (val);
5801 }
5802
5803 int
5804 scm_is_real (SCM val)
5805 {
5806 return scm_is_true (scm_real_p (val));
5807 }
5808
5809 int
5810 scm_is_rational (SCM val)
5811 {
5812 return scm_is_true (scm_rational_p (val));
5813 }
5814
5815 double
5816 scm_to_double (SCM val)
5817 {
5818 if (SCM_I_INUMP (val))
5819 return SCM_I_INUM (val);
5820 else if (SCM_BIGP (val))
5821 return scm_i_big2dbl (val);
5822 else if (SCM_FRACTIONP (val))
5823 return scm_i_fraction2double (val);
5824 else if (SCM_REALP (val))
5825 return SCM_REAL_VALUE (val);
5826 else
5827 scm_wrong_type_arg_msg (NULL, 0, val, "real number");
5828 }
5829
5830 SCM
5831 scm_from_double (double val)
5832 {
5833 SCM z = scm_double_cell (scm_tc16_real, 0, 0, 0);
5834 SCM_REAL_VALUE (z) = val;
5835 return z;
5836 }
5837
5838 #if SCM_ENABLE_DISCOURAGED == 1
5839
5840 float
5841 scm_num2float (SCM num, unsigned long int pos, const char *s_caller)
5842 {
5843 if (SCM_BIGP (num))
5844 {
5845 float res = mpz_get_d (SCM_I_BIG_MPZ (num));
5846 if (!xisinf (res))
5847 return res;
5848 else
5849 scm_out_of_range (NULL, num);
5850 }
5851 else
5852 return scm_to_double (num);
5853 }
5854
5855 double
5856 scm_num2double (SCM num, unsigned long int pos, const char *s_caller)
5857 {
5858 if (SCM_BIGP (num))
5859 {
5860 double res = mpz_get_d (SCM_I_BIG_MPZ (num));
5861 if (!xisinf (res))
5862 return res;
5863 else
5864 scm_out_of_range (NULL, num);
5865 }
5866 else
5867 return scm_to_double (num);
5868 }
5869
5870 #endif
5871
5872 int
5873 scm_is_complex (SCM val)
5874 {
5875 return scm_is_true (scm_complex_p (val));
5876 }
5877
5878 double
5879 scm_c_real_part (SCM z)
5880 {
5881 if (SCM_COMPLEXP (z))
5882 return SCM_COMPLEX_REAL (z);
5883 else
5884 {
5885 /* Use the scm_real_part to get proper error checking and
5886 dispatching.
5887 */
5888 return scm_to_double (scm_real_part (z));
5889 }
5890 }
5891
5892 double
5893 scm_c_imag_part (SCM z)
5894 {
5895 if (SCM_COMPLEXP (z))
5896 return SCM_COMPLEX_IMAG (z);
5897 else
5898 {
5899 /* Use the scm_imag_part to get proper error checking and
5900 dispatching. The result will almost always be 0.0, but not
5901 always.
5902 */
5903 return scm_to_double (scm_imag_part (z));
5904 }
5905 }
5906
5907 double
5908 scm_c_magnitude (SCM z)
5909 {
5910 return scm_to_double (scm_magnitude (z));
5911 }
5912
5913 double
5914 scm_c_angle (SCM z)
5915 {
5916 return scm_to_double (scm_angle (z));
5917 }
5918
5919 int
5920 scm_is_number (SCM z)
5921 {
5922 return scm_is_true (scm_number_p (z));
5923 }
5924
5925 void
5926 scm_init_numbers ()
5927 {
5928 int i;
5929
5930 mpz_init_set_si (z_negative_one, -1);
5931
5932 /* It may be possible to tune the performance of some algorithms by using
5933 * the following constants to avoid the creation of bignums. Please, before
5934 * using these values, remember the two rules of program optimization:
5935 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5936 scm_c_define ("most-positive-fixnum",
5937 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM));
5938 scm_c_define ("most-negative-fixnum",
5939 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM));
5940
5941 scm_add_feature ("complex");
5942 scm_add_feature ("inexact");
5943 scm_flo0 = scm_from_double (0.0);
5944
5945 /* determine floating point precision */
5946 for (i=2; i <= SCM_MAX_DBL_RADIX; ++i)
5947 {
5948 init_dblprec(&scm_dblprec[i-2],i);
5949 init_fx_radix(fx_per_radix[i-2],i);
5950 }
5951 #ifdef DBL_DIG
5952 /* hard code precision for base 10 if the preprocessor tells us to... */
5953 scm_dblprec[10-2] = (DBL_DIG > 20) ? 20 : DBL_DIG;
5954 #endif
5955
5956 exactly_one_half = scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
5957 SCM_I_MAKINUM (2)));
5958 #include "libguile/numbers.x"
5959 }
5960
5961 /*
5962 Local Variables:
5963 c-file-style: "gnu"
5964 End:
5965 */