1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public License
9 * as published by the Free Software Foundation; either version 3 of
10 * the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful, but
13 * 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.
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., 51 Franklin Street, Fifth Floor, Boston, MA
24 /* General assumptions:
25 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
26 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
27 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
28 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
29 * All objects satisfying SCM_FRACTIONP are never an integer.
34 - see if special casing bignums and reals in integer-exponent when
35 possible (to use mpz_pow and mpf_pow_ui) is faster.
37 - look in to better short-circuiting of common cases in
38 integer-expt and elsewhere.
40 - see if direct mpz operations can help in ash and elsewhere.
57 #include "libguile/_scm.h"
58 #include "libguile/feature.h"
59 #include "libguile/ports.h"
60 #include "libguile/root.h"
61 #include "libguile/smob.h"
62 #include "libguile/strings.h"
63 #include "libguile/bdw-gc.h"
65 #include "libguile/validate.h"
66 #include "libguile/numbers.h"
67 #include "libguile/deprecation.h"
69 #include "libguile/eq.h"
71 /* values per glibc, if not already defined */
73 #define M_LOG10E 0.43429448190325182765
76 #define M_PI 3.14159265358979323846
79 typedef scm_t_signed_bits scm_t_inum
;
80 #define scm_from_inum(x) (scm_from_signed_integer (x))
85 Wonder if this might be faster for some of our code? A switch on
86 the numtag would jump directly to the right case, and the
87 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
89 #define SCM_I_NUMTAG_NOTNUM 0
90 #define SCM_I_NUMTAG_INUM 1
91 #define SCM_I_NUMTAG_BIG scm_tc16_big
92 #define SCM_I_NUMTAG_REAL scm_tc16_real
93 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
94 #define SCM_I_NUMTAG(x) \
95 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
96 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
97 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
98 : SCM_I_NUMTAG_NOTNUM)))
100 /* the macro above will not work as is with fractions */
105 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
107 /* FLOBUFLEN is the maximum number of characters neccessary for the
108 * printed or scm_string representation of an inexact number.
110 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
113 #if !defined (HAVE_ASINH)
114 static double asinh (double x
) { return log (x
+ sqrt (x
* x
+ 1)); }
116 #if !defined (HAVE_ACOSH)
117 static double acosh (double x
) { return log (x
+ sqrt (x
* x
- 1)); }
119 #if !defined (HAVE_ATANH)
120 static double atanh (double x
) { return 0.5 * log ((1 + x
) / (1 - x
)); }
123 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
124 an explicit check. In some future gmp (don't know what version number),
125 mpz_cmp_d is supposed to do this itself. */
127 #define xmpz_cmp_d(z, d) \
128 (isinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
130 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
134 #if defined (GUILE_I)
135 #if HAVE_COMPLEX_DOUBLE
137 /* For an SCM object Z which is a complex number (ie. satisfies
138 SCM_COMPLEXP), return its value as a C level "complex double". */
139 #define SCM_COMPLEX_VALUE(z) \
140 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
142 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
144 /* Convert a C "complex double" to an SCM value. */
146 scm_from_complex_double (complex double z
)
148 return scm_c_make_rectangular (creal (z
), cimag (z
));
151 #endif /* HAVE_COMPLEX_DOUBLE */
156 static mpz_t z_negative_one
;
159 /* Clear the `mpz_t' embedded in bignum PTR. */
161 finalize_bignum (GC_PTR ptr
, GC_PTR data
)
165 bignum
= PTR2SCM (ptr
);
166 mpz_clear (SCM_I_BIG_MPZ (bignum
));
169 /* Return a new uninitialized bignum. */
174 GC_finalization_proc prev_finalizer
;
175 GC_PTR prev_finalizer_data
;
177 /* Allocate one word for the type tag and enough room for an `mpz_t'. */
178 p
= scm_gc_malloc_pointerless (sizeof (scm_t_bits
) + sizeof (mpz_t
),
182 GC_REGISTER_FINALIZER_NO_ORDER (p
, finalize_bignum
, NULL
,
184 &prev_finalizer_data
);
193 /* Return a newly created bignum. */
194 SCM z
= make_bignum ();
195 mpz_init (SCM_I_BIG_MPZ (z
));
200 scm_i_inum2big (scm_t_inum x
)
202 /* Return a newly created bignum initialized to X. */
203 SCM z
= make_bignum ();
204 #if SIZEOF_VOID_P == SIZEOF_LONG
205 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
207 /* Note that in this case, you'll also have to check all mpz_*_ui and
208 mpz_*_si invocations in Guile. */
209 #error creation of mpz not implemented for this inum size
215 scm_i_long2big (long x
)
217 /* Return a newly created bignum initialized to X. */
218 SCM z
= make_bignum ();
219 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
224 scm_i_ulong2big (unsigned long x
)
226 /* Return a newly created bignum initialized to X. */
227 SCM z
= make_bignum ();
228 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
233 scm_i_clonebig (SCM src_big
, int same_sign_p
)
235 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
236 SCM z
= make_bignum ();
237 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
239 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
244 scm_i_bigcmp (SCM x
, SCM y
)
246 /* Return neg if x < y, pos if x > y, and 0 if x == y */
247 /* presume we already know x and y are bignums */
248 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
249 scm_remember_upto_here_2 (x
, y
);
254 scm_i_dbl2big (double d
)
256 /* results are only defined if d is an integer */
257 SCM z
= make_bignum ();
258 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
262 /* Convert a integer in double representation to a SCM number. */
265 scm_i_dbl2num (double u
)
267 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
268 powers of 2, so there's no rounding when making "double" values
269 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
270 get rounded on a 64-bit machine, hence the "+1".
272 The use of floor() to force to an integer value ensures we get a
273 "numerically closest" value without depending on how a
274 double->long cast or how mpz_set_d will round. For reference,
275 double->long probably follows the hardware rounding mode,
276 mpz_set_d truncates towards zero. */
278 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
279 representable as a double? */
281 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
282 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
283 return SCM_I_MAKINUM ((scm_t_inum
) u
);
285 return scm_i_dbl2big (u
);
288 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
289 with R5RS exact->inexact.
291 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
292 (ie. truncate towards zero), then adjust to get the closest double by
293 examining the next lower bit and adding 1 (to the absolute value) if
296 Bignums exactly half way between representable doubles are rounded to the
297 next higher absolute value (ie. away from zero). This seems like an
298 adequate interpretation of R5RS "numerically closest", and it's easier
299 and faster than a full "nearest-even" style.
301 The bit test must be done on the absolute value of the mpz_t, which means
302 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
303 negatives as twos complement.
305 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
306 following the hardware rounding mode, but applied to the absolute value
307 of the mpz_t operand. This is not what we want so we put the high
308 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
309 mpz_get_d is supposed to always truncate towards zero.
311 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
312 is a slowdown. It'd be faster to pick out the relevant high bits with
313 mpz_getlimbn if we could be bothered coding that, and if the new
314 truncating gmp doesn't come out. */
317 scm_i_big2dbl (SCM b
)
322 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
326 /* Current GMP, eg. 4.1.3, force truncation towards zero */
328 if (bits
> DBL_MANT_DIG
)
330 size_t shift
= bits
- DBL_MANT_DIG
;
331 mpz_init2 (tmp
, DBL_MANT_DIG
);
332 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
333 result
= ldexp (mpz_get_d (tmp
), shift
);
338 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
343 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
346 if (bits
> DBL_MANT_DIG
)
348 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
349 /* test bit number "pos" in absolute value */
350 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
351 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
353 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
357 scm_remember_upto_here_1 (b
);
362 scm_i_normbig (SCM b
)
364 /* convert a big back to a fixnum if it'll fit */
365 /* presume b is a bignum */
366 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
368 scm_t_inum val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
369 if (SCM_FIXABLE (val
))
370 b
= SCM_I_MAKINUM (val
);
375 static SCM_C_INLINE_KEYWORD SCM
376 scm_i_mpz2num (mpz_t b
)
378 /* convert a mpz number to a SCM number. */
379 if (mpz_fits_slong_p (b
))
381 scm_t_inum val
= mpz_get_si (b
);
382 if (SCM_FIXABLE (val
))
383 return SCM_I_MAKINUM (val
);
387 SCM z
= make_bignum ();
388 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
393 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
394 static SCM
scm_divide2real (SCM x
, SCM y
);
397 scm_i_make_ratio (SCM numerator
, SCM denominator
)
398 #define FUNC_NAME "make-ratio"
400 /* First make sure the arguments are proper.
402 if (SCM_I_INUMP (denominator
))
404 if (scm_is_eq (denominator
, SCM_INUM0
))
405 scm_num_overflow ("make-ratio");
406 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
411 if (!(SCM_BIGP(denominator
)))
412 SCM_WRONG_TYPE_ARG (2, denominator
);
414 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
415 SCM_WRONG_TYPE_ARG (1, numerator
);
417 /* Then flip signs so that the denominator is positive.
419 if (scm_is_true (scm_negative_p (denominator
)))
421 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
422 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
425 /* Now consider for each of the four fixnum/bignum combinations
426 whether the rational number is really an integer.
428 if (SCM_I_INUMP (numerator
))
430 scm_t_inum x
= SCM_I_INUM (numerator
);
431 if (scm_is_eq (numerator
, SCM_INUM0
))
433 if (SCM_I_INUMP (denominator
))
436 y
= SCM_I_INUM (denominator
);
438 return SCM_I_MAKINUM(1);
440 return SCM_I_MAKINUM (x
/ y
);
444 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
445 of that value for the denominator, as a bignum. Apart from
446 that case, abs(bignum) > abs(inum) so inum/bignum is not an
448 if (x
== SCM_MOST_NEGATIVE_FIXNUM
449 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
450 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
451 return SCM_I_MAKINUM(-1);
454 else if (SCM_BIGP (numerator
))
456 if (SCM_I_INUMP (denominator
))
458 scm_t_inum yy
= SCM_I_INUM (denominator
);
459 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
460 return scm_divide (numerator
, denominator
);
464 if (scm_is_eq (numerator
, denominator
))
465 return SCM_I_MAKINUM(1);
466 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
467 SCM_I_BIG_MPZ (denominator
)))
468 return scm_divide(numerator
, denominator
);
472 /* No, it's a proper fraction.
475 SCM divisor
= scm_gcd (numerator
, denominator
);
476 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
478 numerator
= scm_divide (numerator
, divisor
);
479 denominator
= scm_divide (denominator
, divisor
);
482 return scm_double_cell (scm_tc16_fraction
,
483 SCM_UNPACK (numerator
),
484 SCM_UNPACK (denominator
), 0);
490 scm_i_fraction2double (SCM z
)
492 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
493 SCM_FRACTION_DENOMINATOR (z
)));
496 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
498 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
500 #define FUNC_NAME s_scm_exact_p
506 if (SCM_FRACTIONP (x
))
510 SCM_WRONG_TYPE_ARG (1, x
);
515 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
517 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
519 #define FUNC_NAME s_scm_odd_p
523 scm_t_inum val
= SCM_I_INUM (n
);
524 return scm_from_bool ((val
& 1L) != 0);
526 else if (SCM_BIGP (n
))
528 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
529 scm_remember_upto_here_1 (n
);
530 return scm_from_bool (odd_p
);
532 else if (scm_is_true (scm_inf_p (n
)))
533 SCM_WRONG_TYPE_ARG (1, n
);
534 else if (SCM_REALP (n
))
536 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
542 SCM_WRONG_TYPE_ARG (1, n
);
545 SCM_WRONG_TYPE_ARG (1, n
);
550 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
552 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
554 #define FUNC_NAME s_scm_even_p
558 scm_t_inum val
= SCM_I_INUM (n
);
559 return scm_from_bool ((val
& 1L) == 0);
561 else if (SCM_BIGP (n
))
563 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
564 scm_remember_upto_here_1 (n
);
565 return scm_from_bool (even_p
);
567 else if (scm_is_true (scm_inf_p (n
)))
568 SCM_WRONG_TYPE_ARG (1, n
);
569 else if (SCM_REALP (n
))
571 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
577 SCM_WRONG_TYPE_ARG (1, n
);
580 SCM_WRONG_TYPE_ARG (1, n
);
584 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
586 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
587 "or @samp{-inf.0}, @code{#f} otherwise.")
588 #define FUNC_NAME s_scm_inf_p
591 return scm_from_bool (isinf (SCM_REAL_VALUE (x
)));
592 else if (SCM_COMPLEXP (x
))
593 return scm_from_bool (isinf (SCM_COMPLEX_REAL (x
))
594 || isinf (SCM_COMPLEX_IMAG (x
)));
600 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
602 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
604 #define FUNC_NAME s_scm_nan_p
607 return scm_from_bool (isnan (SCM_REAL_VALUE (n
)));
608 else if (SCM_COMPLEXP (n
))
609 return scm_from_bool (isnan (SCM_COMPLEX_REAL (n
))
610 || isnan (SCM_COMPLEX_IMAG (n
)));
616 /* Guile's idea of infinity. */
617 static double guile_Inf
;
619 /* Guile's idea of not a number. */
620 static double guile_NaN
;
623 guile_ieee_init (void)
625 /* Some version of gcc on some old version of Linux used to crash when
626 trying to make Inf and NaN. */
629 /* C99 INFINITY, when available.
630 FIXME: The standard allows for INFINITY to be something that overflows
631 at compile time. We ought to have a configure test to check for that
632 before trying to use it. (But in practice we believe this is not a
633 problem on any system guile is likely to target.) */
634 guile_Inf
= INFINITY
;
635 #elif defined HAVE_DINFINITY
637 extern unsigned int DINFINITY
[2];
638 guile_Inf
= (*((double *) (DINFINITY
)));
645 if (guile_Inf
== tmp
)
652 /* C99 NAN, when available */
654 #elif defined HAVE_DQNAN
657 extern unsigned int DQNAN
[2];
658 guile_NaN
= (*((double *)(DQNAN
)));
661 guile_NaN
= guile_Inf
/ guile_Inf
;
665 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
668 #define FUNC_NAME s_scm_inf
670 static int initialized
= 0;
676 return scm_from_double (guile_Inf
);
680 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
683 #define FUNC_NAME s_scm_nan
685 static int initialized
= 0;
691 return scm_from_double (guile_NaN
);
696 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
698 "Return the absolute value of @var{x}.")
703 scm_t_inum xx
= SCM_I_INUM (x
);
706 else if (SCM_POSFIXABLE (-xx
))
707 return SCM_I_MAKINUM (-xx
);
709 return scm_i_inum2big (-xx
);
711 else if (SCM_BIGP (x
))
713 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
715 return scm_i_clonebig (x
, 0);
719 else if (SCM_REALP (x
))
721 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
722 double xx
= SCM_REAL_VALUE (x
);
724 return scm_from_double (-xx
);
728 else if (SCM_FRACTIONP (x
))
730 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
732 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
733 SCM_FRACTION_DENOMINATOR (x
));
736 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
741 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
742 /* "Return the quotient of the numbers @var{x} and @var{y}."
745 scm_quotient (SCM x
, SCM y
)
749 scm_t_inum xx
= SCM_I_INUM (x
);
752 scm_t_inum yy
= SCM_I_INUM (y
);
754 scm_num_overflow (s_quotient
);
757 scm_t_inum z
= xx
/ yy
;
759 return SCM_I_MAKINUM (z
);
761 return scm_i_inum2big (z
);
764 else if (SCM_BIGP (y
))
766 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
767 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
768 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
770 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
771 scm_remember_upto_here_1 (y
);
772 return SCM_I_MAKINUM (-1);
775 return SCM_I_MAKINUM (0);
778 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
780 else if (SCM_BIGP (x
))
784 scm_t_inum yy
= SCM_I_INUM (y
);
786 scm_num_overflow (s_quotient
);
791 SCM result
= scm_i_mkbig ();
794 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
797 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
800 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
801 scm_remember_upto_here_1 (x
);
802 return scm_i_normbig (result
);
805 else if (SCM_BIGP (y
))
807 SCM result
= scm_i_mkbig ();
808 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
811 scm_remember_upto_here_2 (x
, y
);
812 return scm_i_normbig (result
);
815 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
818 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
821 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
822 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
824 * "(remainder 13 4) @result{} 1\n"
825 * "(remainder -13 4) @result{} -1\n"
829 scm_remainder (SCM x
, SCM y
)
835 scm_t_inum yy
= SCM_I_INUM (y
);
837 scm_num_overflow (s_remainder
);
840 scm_t_inum z
= SCM_I_INUM (x
) % yy
;
841 return SCM_I_MAKINUM (z
);
844 else if (SCM_BIGP (y
))
846 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
847 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
848 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
850 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
851 scm_remember_upto_here_1 (y
);
852 return SCM_I_MAKINUM (0);
858 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
860 else if (SCM_BIGP (x
))
864 scm_t_inum yy
= SCM_I_INUM (y
);
866 scm_num_overflow (s_remainder
);
869 SCM result
= scm_i_mkbig ();
872 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
873 scm_remember_upto_here_1 (x
);
874 return scm_i_normbig (result
);
877 else if (SCM_BIGP (y
))
879 SCM result
= scm_i_mkbig ();
880 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
883 scm_remember_upto_here_2 (x
, y
);
884 return scm_i_normbig (result
);
887 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
890 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
894 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
895 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
897 * "(modulo 13 4) @result{} 1\n"
898 * "(modulo -13 4) @result{} 3\n"
902 scm_modulo (SCM x
, SCM y
)
906 scm_t_inum xx
= SCM_I_INUM (x
);
909 scm_t_inum yy
= SCM_I_INUM (y
);
911 scm_num_overflow (s_modulo
);
914 /* C99 specifies that "%" is the remainder corresponding to a
915 quotient rounded towards zero, and that's also traditional
916 for machine division, so z here should be well defined. */
917 scm_t_inum z
= xx
% yy
;
934 return SCM_I_MAKINUM (result
);
937 else if (SCM_BIGP (y
))
939 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
946 SCM pos_y
= scm_i_clonebig (y
, 0);
947 /* do this after the last scm_op */
948 mpz_init_set_si (z_x
, xx
);
949 result
= pos_y
; /* re-use this bignum */
950 mpz_mod (SCM_I_BIG_MPZ (result
),
952 SCM_I_BIG_MPZ (pos_y
));
953 scm_remember_upto_here_1 (pos_y
);
957 result
= scm_i_mkbig ();
958 /* do this after the last scm_op */
959 mpz_init_set_si (z_x
, xx
);
960 mpz_mod (SCM_I_BIG_MPZ (result
),
963 scm_remember_upto_here_1 (y
);
966 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
967 mpz_add (SCM_I_BIG_MPZ (result
),
969 SCM_I_BIG_MPZ (result
));
970 scm_remember_upto_here_1 (y
);
971 /* and do this before the next one */
973 return scm_i_normbig (result
);
977 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
979 else if (SCM_BIGP (x
))
983 scm_t_inum yy
= SCM_I_INUM (y
);
985 scm_num_overflow (s_modulo
);
988 SCM result
= scm_i_mkbig ();
989 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
991 (yy
< 0) ? - yy
: yy
);
992 scm_remember_upto_here_1 (x
);
993 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
994 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
995 SCM_I_BIG_MPZ (result
),
997 return scm_i_normbig (result
);
1000 else if (SCM_BIGP (y
))
1003 SCM result
= scm_i_mkbig ();
1004 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
1005 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1006 mpz_mod (SCM_I_BIG_MPZ (result
),
1008 SCM_I_BIG_MPZ (pos_y
));
1010 scm_remember_upto_here_1 (x
);
1011 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1012 mpz_add (SCM_I_BIG_MPZ (result
),
1014 SCM_I_BIG_MPZ (result
));
1015 scm_remember_upto_here_2 (y
, pos_y
);
1016 return scm_i_normbig (result
);
1020 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1023 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1026 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
1027 (SCM x
, SCM y
, SCM rest
),
1028 "Return the greatest common divisor of all parameter values.\n"
1029 "If called without arguments, 0 is returned.")
1030 #define FUNC_NAME s_scm_i_gcd
1032 while (!scm_is_null (rest
))
1033 { x
= scm_gcd (x
, y
);
1035 rest
= scm_cdr (rest
);
1037 return scm_gcd (x
, y
);
1041 #define s_gcd s_scm_i_gcd
1042 #define g_gcd g_scm_i_gcd
1045 scm_gcd (SCM x
, SCM y
)
1048 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
1050 if (SCM_I_INUMP (x
))
1052 if (SCM_I_INUMP (y
))
1054 scm_t_inum xx
= SCM_I_INUM (x
);
1055 scm_t_inum yy
= SCM_I_INUM (y
);
1056 scm_t_inum u
= xx
< 0 ? -xx
: xx
;
1057 scm_t_inum v
= yy
< 0 ? -yy
: yy
;
1067 /* Determine a common factor 2^k */
1068 while (!(1 & (u
| v
)))
1074 /* Now, any factor 2^n can be eliminated */
1094 return (SCM_POSFIXABLE (result
)
1095 ? SCM_I_MAKINUM (result
)
1096 : scm_i_inum2big (result
));
1098 else if (SCM_BIGP (y
))
1104 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1106 else if (SCM_BIGP (x
))
1108 if (SCM_I_INUMP (y
))
1113 yy
= SCM_I_INUM (y
);
1118 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1119 scm_remember_upto_here_1 (x
);
1120 return (SCM_POSFIXABLE (result
)
1121 ? SCM_I_MAKINUM (result
)
1122 : scm_from_unsigned_integer (result
));
1124 else if (SCM_BIGP (y
))
1126 SCM result
= scm_i_mkbig ();
1127 mpz_gcd (SCM_I_BIG_MPZ (result
),
1130 scm_remember_upto_here_2 (x
, y
);
1131 return scm_i_normbig (result
);
1134 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1137 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1140 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
1141 (SCM x
, SCM y
, SCM rest
),
1142 "Return the least common multiple of the arguments.\n"
1143 "If called without arguments, 1 is returned.")
1144 #define FUNC_NAME s_scm_i_lcm
1146 while (!scm_is_null (rest
))
1147 { x
= scm_lcm (x
, y
);
1149 rest
= scm_cdr (rest
);
1151 return scm_lcm (x
, y
);
1155 #define s_lcm s_scm_i_lcm
1156 #define g_lcm g_scm_i_lcm
1159 scm_lcm (SCM n1
, SCM n2
)
1161 if (SCM_UNBNDP (n2
))
1163 if (SCM_UNBNDP (n1
))
1164 return SCM_I_MAKINUM (1L);
1165 n2
= SCM_I_MAKINUM (1L);
1168 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1169 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1170 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1171 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1173 if (SCM_I_INUMP (n1
))
1175 if (SCM_I_INUMP (n2
))
1177 SCM d
= scm_gcd (n1
, n2
);
1178 if (scm_is_eq (d
, SCM_INUM0
))
1181 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1185 /* inum n1, big n2 */
1188 SCM result
= scm_i_mkbig ();
1189 scm_t_inum nn1
= SCM_I_INUM (n1
);
1190 if (nn1
== 0) return SCM_INUM0
;
1191 if (nn1
< 0) nn1
= - nn1
;
1192 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1193 scm_remember_upto_here_1 (n2
);
1201 if (SCM_I_INUMP (n2
))
1208 SCM result
= scm_i_mkbig ();
1209 mpz_lcm(SCM_I_BIG_MPZ (result
),
1211 SCM_I_BIG_MPZ (n2
));
1212 scm_remember_upto_here_2(n1
, n2
);
1213 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1219 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1224 + + + x (map digit:logand X Y)
1225 + - + x (map digit:logand X (lognot (+ -1 Y)))
1226 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1227 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1232 + + + (map digit:logior X Y)
1233 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1234 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1235 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1240 + + + (map digit:logxor X Y)
1241 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1242 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1243 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1248 + + (any digit:logand X Y)
1249 + - (any digit:logand X (lognot (+ -1 Y)))
1250 - + (any digit:logand (lognot (+ -1 X)) Y)
1255 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
1256 (SCM x
, SCM y
, SCM rest
),
1257 "Return the bitwise AND of the integer arguments.\n\n"
1259 "(logand) @result{} -1\n"
1260 "(logand 7) @result{} 7\n"
1261 "(logand #b111 #b011 #b001) @result{} 1\n"
1263 #define FUNC_NAME s_scm_i_logand
1265 while (!scm_is_null (rest
))
1266 { x
= scm_logand (x
, y
);
1268 rest
= scm_cdr (rest
);
1270 return scm_logand (x
, y
);
1274 #define s_scm_logand s_scm_i_logand
1276 SCM
scm_logand (SCM n1
, SCM n2
)
1277 #define FUNC_NAME s_scm_logand
1281 if (SCM_UNBNDP (n2
))
1283 if (SCM_UNBNDP (n1
))
1284 return SCM_I_MAKINUM (-1);
1285 else if (!SCM_NUMBERP (n1
))
1286 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1287 else if (SCM_NUMBERP (n1
))
1290 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1293 if (SCM_I_INUMP (n1
))
1295 nn1
= SCM_I_INUM (n1
);
1296 if (SCM_I_INUMP (n2
))
1298 scm_t_inum nn2
= SCM_I_INUM (n2
);
1299 return SCM_I_MAKINUM (nn1
& nn2
);
1301 else if SCM_BIGP (n2
)
1307 SCM result_z
= scm_i_mkbig ();
1309 mpz_init_set_si (nn1_z
, nn1
);
1310 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1311 scm_remember_upto_here_1 (n2
);
1313 return scm_i_normbig (result_z
);
1317 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1319 else if (SCM_BIGP (n1
))
1321 if (SCM_I_INUMP (n2
))
1324 nn1
= SCM_I_INUM (n1
);
1327 else if (SCM_BIGP (n2
))
1329 SCM result_z
= scm_i_mkbig ();
1330 mpz_and (SCM_I_BIG_MPZ (result_z
),
1332 SCM_I_BIG_MPZ (n2
));
1333 scm_remember_upto_here_2 (n1
, n2
);
1334 return scm_i_normbig (result_z
);
1337 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1340 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1345 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
1346 (SCM x
, SCM y
, SCM rest
),
1347 "Return the bitwise OR of the integer arguments.\n\n"
1349 "(logior) @result{} 0\n"
1350 "(logior 7) @result{} 7\n"
1351 "(logior #b000 #b001 #b011) @result{} 3\n"
1353 #define FUNC_NAME s_scm_i_logior
1355 while (!scm_is_null (rest
))
1356 { x
= scm_logior (x
, y
);
1358 rest
= scm_cdr (rest
);
1360 return scm_logior (x
, y
);
1364 #define s_scm_logior s_scm_i_logior
1366 SCM
scm_logior (SCM n1
, SCM n2
)
1367 #define FUNC_NAME s_scm_logior
1371 if (SCM_UNBNDP (n2
))
1373 if (SCM_UNBNDP (n1
))
1375 else if (SCM_NUMBERP (n1
))
1378 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1381 if (SCM_I_INUMP (n1
))
1383 nn1
= SCM_I_INUM (n1
);
1384 if (SCM_I_INUMP (n2
))
1386 long nn2
= SCM_I_INUM (n2
);
1387 return SCM_I_MAKINUM (nn1
| nn2
);
1389 else if (SCM_BIGP (n2
))
1395 SCM result_z
= scm_i_mkbig ();
1397 mpz_init_set_si (nn1_z
, nn1
);
1398 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1399 scm_remember_upto_here_1 (n2
);
1401 return scm_i_normbig (result_z
);
1405 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1407 else if (SCM_BIGP (n1
))
1409 if (SCM_I_INUMP (n2
))
1412 nn1
= SCM_I_INUM (n1
);
1415 else if (SCM_BIGP (n2
))
1417 SCM result_z
= scm_i_mkbig ();
1418 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1420 SCM_I_BIG_MPZ (n2
));
1421 scm_remember_upto_here_2 (n1
, n2
);
1422 return scm_i_normbig (result_z
);
1425 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1428 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1433 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
1434 (SCM x
, SCM y
, SCM rest
),
1435 "Return the bitwise XOR of the integer arguments. A bit is\n"
1436 "set in the result if it is set in an odd number of arguments.\n"
1438 "(logxor) @result{} 0\n"
1439 "(logxor 7) @result{} 7\n"
1440 "(logxor #b000 #b001 #b011) @result{} 2\n"
1441 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1443 #define FUNC_NAME s_scm_i_logxor
1445 while (!scm_is_null (rest
))
1446 { x
= scm_logxor (x
, y
);
1448 rest
= scm_cdr (rest
);
1450 return scm_logxor (x
, y
);
1454 #define s_scm_logxor s_scm_i_logxor
1456 SCM
scm_logxor (SCM n1
, SCM n2
)
1457 #define FUNC_NAME s_scm_logxor
1461 if (SCM_UNBNDP (n2
))
1463 if (SCM_UNBNDP (n1
))
1465 else if (SCM_NUMBERP (n1
))
1468 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1471 if (SCM_I_INUMP (n1
))
1473 nn1
= SCM_I_INUM (n1
);
1474 if (SCM_I_INUMP (n2
))
1476 scm_t_inum nn2
= SCM_I_INUM (n2
);
1477 return SCM_I_MAKINUM (nn1
^ nn2
);
1479 else if (SCM_BIGP (n2
))
1483 SCM result_z
= scm_i_mkbig ();
1485 mpz_init_set_si (nn1_z
, nn1
);
1486 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1487 scm_remember_upto_here_1 (n2
);
1489 return scm_i_normbig (result_z
);
1493 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1495 else if (SCM_BIGP (n1
))
1497 if (SCM_I_INUMP (n2
))
1500 nn1
= SCM_I_INUM (n1
);
1503 else if (SCM_BIGP (n2
))
1505 SCM result_z
= scm_i_mkbig ();
1506 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1508 SCM_I_BIG_MPZ (n2
));
1509 scm_remember_upto_here_2 (n1
, n2
);
1510 return scm_i_normbig (result_z
);
1513 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1516 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1521 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1523 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1524 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1525 "without actually calculating the @code{logand}, just testing\n"
1529 "(logtest #b0100 #b1011) @result{} #f\n"
1530 "(logtest #b0100 #b0111) @result{} #t\n"
1532 #define FUNC_NAME s_scm_logtest
1536 if (SCM_I_INUMP (j
))
1538 nj
= SCM_I_INUM (j
);
1539 if (SCM_I_INUMP (k
))
1541 scm_t_inum nk
= SCM_I_INUM (k
);
1542 return scm_from_bool (nj
& nk
);
1544 else if (SCM_BIGP (k
))
1552 mpz_init_set_si (nj_z
, nj
);
1553 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1554 scm_remember_upto_here_1 (k
);
1555 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1561 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1563 else if (SCM_BIGP (j
))
1565 if (SCM_I_INUMP (k
))
1568 nj
= SCM_I_INUM (j
);
1571 else if (SCM_BIGP (k
))
1575 mpz_init (result_z
);
1579 scm_remember_upto_here_2 (j
, k
);
1580 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1581 mpz_clear (result_z
);
1585 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1588 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1593 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1595 "Test whether bit number @var{index} in @var{j} is set.\n"
1596 "@var{index} starts from 0 for the least significant bit.\n"
1599 "(logbit? 0 #b1101) @result{} #t\n"
1600 "(logbit? 1 #b1101) @result{} #f\n"
1601 "(logbit? 2 #b1101) @result{} #t\n"
1602 "(logbit? 3 #b1101) @result{} #t\n"
1603 "(logbit? 4 #b1101) @result{} #f\n"
1605 #define FUNC_NAME s_scm_logbit_p
1607 unsigned long int iindex
;
1608 iindex
= scm_to_ulong (index
);
1610 if (SCM_I_INUMP (j
))
1612 /* bits above what's in an inum follow the sign bit */
1613 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1614 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1616 else if (SCM_BIGP (j
))
1618 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1619 scm_remember_upto_here_1 (j
);
1620 return scm_from_bool (val
);
1623 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1628 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1630 "Return the integer which is the ones-complement of the integer\n"
1634 "(number->string (lognot #b10000000) 2)\n"
1635 " @result{} \"-10000001\"\n"
1636 "(number->string (lognot #b0) 2)\n"
1637 " @result{} \"-1\"\n"
1639 #define FUNC_NAME s_scm_lognot
1641 if (SCM_I_INUMP (n
)) {
1642 /* No overflow here, just need to toggle all the bits making up the inum.
1643 Enhancement: No need to strip the tag and add it back, could just xor
1644 a block of 1 bits, if that worked with the various debug versions of
1646 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1648 } else if (SCM_BIGP (n
)) {
1649 SCM result
= scm_i_mkbig ();
1650 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1651 scm_remember_upto_here_1 (n
);
1655 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1660 /* returns 0 if IN is not an integer. OUT must already be
1663 coerce_to_big (SCM in
, mpz_t out
)
1666 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1667 else if (SCM_I_INUMP (in
))
1668 mpz_set_si (out
, SCM_I_INUM (in
));
1675 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1676 (SCM n
, SCM k
, SCM m
),
1677 "Return @var{n} raised to the integer exponent\n"
1678 "@var{k}, modulo @var{m}.\n"
1681 "(modulo-expt 2 3 5)\n"
1684 #define FUNC_NAME s_scm_modulo_expt
1690 /* There are two classes of error we might encounter --
1691 1) Math errors, which we'll report by calling scm_num_overflow,
1693 2) wrong-type errors, which of course we'll report by calling
1695 We don't report those errors immediately, however; instead we do
1696 some cleanup first. These variables tell us which error (if
1697 any) we should report after cleaning up.
1699 int report_overflow
= 0;
1701 int position_of_wrong_type
= 0;
1702 SCM value_of_wrong_type
= SCM_INUM0
;
1704 SCM result
= SCM_UNDEFINED
;
1710 if (scm_is_eq (m
, SCM_INUM0
))
1712 report_overflow
= 1;
1716 if (!coerce_to_big (n
, n_tmp
))
1718 value_of_wrong_type
= n
;
1719 position_of_wrong_type
= 1;
1723 if (!coerce_to_big (k
, k_tmp
))
1725 value_of_wrong_type
= k
;
1726 position_of_wrong_type
= 2;
1730 if (!coerce_to_big (m
, m_tmp
))
1732 value_of_wrong_type
= m
;
1733 position_of_wrong_type
= 3;
1737 /* if the exponent K is negative, and we simply call mpz_powm, we
1738 will get a divide-by-zero exception when an inverse 1/n mod m
1739 doesn't exist (or is not unique). Since exceptions are hard to
1740 handle, we'll attempt the inversion "by hand" -- that way, we get
1741 a simple failure code, which is easy to handle. */
1743 if (-1 == mpz_sgn (k_tmp
))
1745 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1747 report_overflow
= 1;
1750 mpz_neg (k_tmp
, k_tmp
);
1753 result
= scm_i_mkbig ();
1754 mpz_powm (SCM_I_BIG_MPZ (result
),
1759 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1760 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1767 if (report_overflow
)
1768 scm_num_overflow (FUNC_NAME
);
1770 if (position_of_wrong_type
)
1771 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1772 value_of_wrong_type
);
1774 return scm_i_normbig (result
);
1778 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1780 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1781 "exact integer, @var{n} can be any number.\n"
1783 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1784 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1785 "includes @math{0^0} is 1.\n"
1788 "(integer-expt 2 5) @result{} 32\n"
1789 "(integer-expt -3 3) @result{} -27\n"
1790 "(integer-expt 5 -3) @result{} 1/125\n"
1791 "(integer-expt 0 0) @result{} 1\n"
1793 #define FUNC_NAME s_scm_integer_expt
1796 SCM z_i2
= SCM_BOOL_F
;
1798 SCM acc
= SCM_I_MAKINUM (1L);
1800 SCM_VALIDATE_NUMBER (SCM_ARG1
, n
);
1801 if (!SCM_I_INUMP (k
) && !SCM_BIGP (k
))
1802 SCM_WRONG_TYPE_ARG (2, k
);
1804 if (scm_is_true (scm_zero_p (n
)))
1806 if (scm_is_true (scm_zero_p (k
))) /* 0^0 == 1 per R5RS */
1807 return acc
; /* return exact 1, regardless of n */
1808 else if (scm_is_true (scm_positive_p (k
)))
1810 else /* return NaN for (0 ^ k) for negative k per R6RS */
1813 else if (scm_is_eq (n
, acc
))
1815 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1816 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1818 if (SCM_I_INUMP (k
))
1819 i2
= SCM_I_INUM (k
);
1820 else if (SCM_BIGP (k
))
1822 z_i2
= scm_i_clonebig (k
, 1);
1823 scm_remember_upto_here_1 (k
);
1827 SCM_WRONG_TYPE_ARG (2, k
);
1831 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1833 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1834 n
= scm_divide (n
, SCM_UNDEFINED
);
1838 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1842 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1844 return scm_product (acc
, n
);
1846 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1847 acc
= scm_product (acc
, n
);
1848 n
= scm_product (n
, n
);
1849 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1857 n
= scm_divide (n
, SCM_UNDEFINED
);
1864 return scm_product (acc
, n
);
1866 acc
= scm_product (acc
, n
);
1867 n
= scm_product (n
, n
);
1874 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1876 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1877 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1879 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1880 "@var{cnt} is negative it's a division, rounded towards negative\n"
1881 "infinity. (Note that this is not the same rounding as\n"
1882 "@code{quotient} does.)\n"
1884 "With @var{n} viewed as an infinite precision twos complement,\n"
1885 "@code{ash} means a left shift introducing zero bits, or a right\n"
1886 "shift dropping bits.\n"
1889 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1890 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1892 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1893 "(ash -23 -2) @result{} -6\n"
1895 #define FUNC_NAME s_scm_ash
1898 bits_to_shift
= scm_to_long (cnt
);
1900 if (SCM_I_INUMP (n
))
1902 scm_t_inum nn
= SCM_I_INUM (n
);
1904 if (bits_to_shift
> 0)
1906 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1907 overflow a non-zero fixnum. For smaller shifts we check the
1908 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1909 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1910 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1916 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1918 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1921 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1925 SCM result
= scm_i_inum2big (nn
);
1926 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1933 bits_to_shift
= -bits_to_shift
;
1934 if (bits_to_shift
>= SCM_LONG_BIT
)
1935 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1937 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1941 else if (SCM_BIGP (n
))
1945 if (bits_to_shift
== 0)
1948 result
= scm_i_mkbig ();
1949 if (bits_to_shift
>= 0)
1951 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1957 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1958 we have to allocate a bignum even if the result is going to be a
1960 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1962 return scm_i_normbig (result
);
1968 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1974 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1975 (SCM n
, SCM start
, SCM end
),
1976 "Return the integer composed of the @var{start} (inclusive)\n"
1977 "through @var{end} (exclusive) bits of @var{n}. The\n"
1978 "@var{start}th bit becomes the 0-th bit in the result.\n"
1981 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1982 " @result{} \"1010\"\n"
1983 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1984 " @result{} \"10110\"\n"
1986 #define FUNC_NAME s_scm_bit_extract
1988 unsigned long int istart
, iend
, bits
;
1989 istart
= scm_to_ulong (start
);
1990 iend
= scm_to_ulong (end
);
1991 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1993 /* how many bits to keep */
1994 bits
= iend
- istart
;
1996 if (SCM_I_INUMP (n
))
1998 scm_t_inum in
= SCM_I_INUM (n
);
2000 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
2001 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
2002 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
2004 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
2006 /* Since we emulate two's complement encoded numbers, this
2007 * special case requires us to produce a result that has
2008 * more bits than can be stored in a fixnum.
2010 SCM result
= scm_i_inum2big (in
);
2011 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
2016 /* mask down to requisite bits */
2017 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
2018 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
2020 else if (SCM_BIGP (n
))
2025 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
2029 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
2030 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
2031 such bits into a ulong. */
2032 result
= scm_i_mkbig ();
2033 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
2034 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
2035 result
= scm_i_normbig (result
);
2037 scm_remember_upto_here_1 (n
);
2041 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2046 static const char scm_logtab
[] = {
2047 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
2050 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
2052 "Return the number of bits in integer @var{n}. If integer is\n"
2053 "positive, the 1-bits in its binary representation are counted.\n"
2054 "If negative, the 0-bits in its two's-complement binary\n"
2055 "representation are counted. If 0, 0 is returned.\n"
2058 "(logcount #b10101010)\n"
2065 #define FUNC_NAME s_scm_logcount
2067 if (SCM_I_INUMP (n
))
2069 unsigned long c
= 0;
2070 scm_t_inum nn
= SCM_I_INUM (n
);
2075 c
+= scm_logtab
[15 & nn
];
2078 return SCM_I_MAKINUM (c
);
2080 else if (SCM_BIGP (n
))
2082 unsigned long count
;
2083 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
2084 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
2086 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
2087 scm_remember_upto_here_1 (n
);
2088 return SCM_I_MAKINUM (count
);
2091 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2096 static const char scm_ilentab
[] = {
2097 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2101 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2103 "Return the number of bits necessary to represent @var{n}.\n"
2106 "(integer-length #b10101010)\n"
2108 "(integer-length 0)\n"
2110 "(integer-length #b1111)\n"
2113 #define FUNC_NAME s_scm_integer_length
2115 if (SCM_I_INUMP (n
))
2117 unsigned long c
= 0;
2119 scm_t_inum nn
= SCM_I_INUM (n
);
2125 l
= scm_ilentab
[15 & nn
];
2128 return SCM_I_MAKINUM (c
- 4 + l
);
2130 else if (SCM_BIGP (n
))
2132 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2133 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2134 1 too big, so check for that and adjust. */
2135 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2136 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2137 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2138 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2140 scm_remember_upto_here_1 (n
);
2141 return SCM_I_MAKINUM (size
);
2144 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2148 /*** NUMBERS -> STRINGS ***/
2149 #define SCM_MAX_DBL_PREC 60
2150 #define SCM_MAX_DBL_RADIX 36
2152 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2153 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2154 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2157 void init_dblprec(int *prec
, int radix
) {
2158 /* determine floating point precision by adding successively
2159 smaller increments to 1.0 until it is considered == 1.0 */
2160 double f
= ((double)1.0)/radix
;
2161 double fsum
= 1.0 + f
;
2166 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2178 void init_fx_radix(double *fx_list
, int radix
)
2180 /* initialize a per-radix list of tolerances. When added
2181 to a number < 1.0, we can determine if we should raund
2182 up and quit converting a number to a string. */
2186 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2187 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2190 /* use this array as a way to generate a single digit */
2191 static const char number_chars
[] = "0123456789abcdefghijklmnopqrstuvwxyz";
2194 idbl2str (double f
, char *a
, int radix
)
2196 int efmt
, dpt
, d
, i
, wp
;
2198 #ifdef DBL_MIN_10_EXP
2201 #endif /* DBL_MIN_10_EXP */
2206 radix
> SCM_MAX_DBL_RADIX
)
2208 /* revert to existing behavior */
2212 wp
= scm_dblprec
[radix
-2];
2213 fx
= fx_per_radix
[radix
-2];
2217 #ifdef HAVE_COPYSIGN
2218 double sgn
= copysign (1.0, f
);
2223 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2229 strcpy (a
, "-inf.0");
2231 strcpy (a
, "+inf.0");
2236 strcpy (a
, "+nan.0");
2246 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2247 make-uniform-vector, from causing infinite loops. */
2248 /* just do the checking...if it passes, we do the conversion for our
2249 radix again below */
2256 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2264 while (f_cpy
> 10.0)
2267 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2288 if (f
+ fx
[wp
] >= radix
)
2295 /* adding 9999 makes this equivalent to abs(x) % 3 */
2296 dpt
= (exp
+ 9999) % 3;
2300 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2322 a
[ch
++] = number_chars
[d
];
2325 if (f
+ fx
[wp
] >= 1.0)
2327 a
[ch
- 1] = number_chars
[d
+1];
2339 if ((dpt
> 4) && (exp
> 6))
2341 d
= (a
[0] == '-' ? 2 : 1);
2342 for (i
= ch
++; i
> d
; i
--)
2355 if (a
[ch
- 1] == '.')
2356 a
[ch
++] = '0'; /* trailing zero */
2365 for (i
= radix
; i
<= exp
; i
*= radix
);
2366 for (i
/= radix
; i
; i
/= radix
)
2368 a
[ch
++] = number_chars
[exp
/ i
];
2377 icmplx2str (double real
, double imag
, char *str
, int radix
)
2381 i
= idbl2str (real
, str
, radix
);
2384 /* Don't output a '+' for negative numbers or for Inf and
2385 NaN. They will provide their own sign. */
2386 if (0 <= imag
&& !isinf (imag
) && !isnan (imag
))
2388 i
+= idbl2str (imag
, &str
[i
], radix
);
2395 iflo2str (SCM flt
, char *str
, int radix
)
2398 if (SCM_REALP (flt
))
2399 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2401 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2406 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2407 characters in the result.
2409 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2411 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2416 return scm_iuint2str (-num
, rad
, p
) + 1;
2419 return scm_iuint2str (num
, rad
, p
);
2422 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2423 characters in the result.
2425 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2427 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2431 scm_t_uintmax n
= num
;
2433 if (rad
< 2 || rad
> 36)
2434 scm_out_of_range ("scm_iuint2str", scm_from_int (rad
));
2436 for (n
/= rad
; n
> 0; n
/= rad
)
2446 p
[i
] = number_chars
[d
];
2451 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2453 "Return a string holding the external representation of the\n"
2454 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2455 "inexact, a radix of 10 will be used.")
2456 #define FUNC_NAME s_scm_number_to_string
2460 if (SCM_UNBNDP (radix
))
2463 base
= scm_to_signed_integer (radix
, 2, 36);
2465 if (SCM_I_INUMP (n
))
2467 char num_buf
[SCM_INTBUFLEN
];
2468 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2469 return scm_from_locale_stringn (num_buf
, length
);
2471 else if (SCM_BIGP (n
))
2473 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2474 scm_remember_upto_here_1 (n
);
2475 return scm_take_locale_string (str
);
2477 else if (SCM_FRACTIONP (n
))
2479 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2480 scm_from_locale_string ("/"),
2481 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2483 else if (SCM_INEXACTP (n
))
2485 char num_buf
[FLOBUFLEN
];
2486 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2489 SCM_WRONG_TYPE_ARG (1, n
);
2494 /* These print routines used to be stubbed here so that scm_repl.c
2495 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2498 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2500 char num_buf
[FLOBUFLEN
];
2501 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2506 scm_i_print_double (double val
, SCM port
)
2508 char num_buf
[FLOBUFLEN
];
2509 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2513 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2516 char num_buf
[FLOBUFLEN
];
2517 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2522 scm_i_print_complex (double real
, double imag
, SCM port
)
2524 char num_buf
[FLOBUFLEN
];
2525 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2529 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2532 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2533 scm_lfwrite_str (str
, port
);
2534 scm_remember_upto_here_1 (str
);
2539 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2541 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2542 scm_remember_upto_here_1 (exp
);
2543 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2547 /*** END nums->strs ***/
2550 /*** STRINGS -> NUMBERS ***/
2552 /* The following functions implement the conversion from strings to numbers.
2553 * The implementation somehow follows the grammar for numbers as it is given
2554 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2555 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2556 * points should be noted about the implementation:
2557 * * Each function keeps a local index variable 'idx' that points at the
2558 * current position within the parsed string. The global index is only
2559 * updated if the function could parse the corresponding syntactic unit
2561 * * Similarly, the functions keep track of indicators of inexactness ('#',
2562 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2563 * global exactness information is only updated after each part has been
2564 * successfully parsed.
2565 * * Sequences of digits are parsed into temporary variables holding fixnums.
2566 * Only if these fixnums would overflow, the result variables are updated
2567 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2568 * the temporary variables holding the fixnums are cleared, and the process
2569 * starts over again. If for example fixnums were able to store five decimal
2570 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2571 * and the result was computed as 12345 * 100000 + 67890. In other words,
2572 * only every five digits two bignum operations were performed.
2575 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2577 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2579 /* Caller is responsible for checking that the return value is in range
2580 for the given radix, which should be <= 36. */
2582 char_decimal_value (scm_t_uint32 c
)
2584 /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
2585 that's certainly above any valid decimal, so we take advantage of
2586 that to elide some tests. */
2587 unsigned int d
= (unsigned int) uc_decimal_value (c
);
2589 /* If that failed, try extended hexadecimals, then. Only accept ascii
2594 if (c
>= (scm_t_uint32
) 'a')
2595 d
= c
- (scm_t_uint32
)'a' + 10U;
2601 mem2uinteger (SCM mem
, unsigned int *p_idx
,
2602 unsigned int radix
, enum t_exactness
*p_exactness
)
2604 unsigned int idx
= *p_idx
;
2605 unsigned int hash_seen
= 0;
2606 scm_t_bits shift
= 1;
2608 unsigned int digit_value
;
2611 size_t len
= scm_i_string_length (mem
);
2616 c
= scm_i_string_ref (mem
, idx
);
2617 digit_value
= char_decimal_value (c
);
2618 if (digit_value
>= radix
)
2622 result
= SCM_I_MAKINUM (digit_value
);
2625 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2635 digit_value
= char_decimal_value (c
);
2636 /* This check catches non-decimals in addition to out-of-range
2638 if (digit_value
>= radix
)
2643 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2645 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2647 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2654 shift
= shift
* radix
;
2655 add
= add
* radix
+ digit_value
;
2660 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2662 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2666 *p_exactness
= INEXACT
;
2672 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2673 * covers the parts of the rules that start at a potential point. The value
2674 * of the digits up to the point have been parsed by the caller and are given
2675 * in variable result. The content of *p_exactness indicates, whether a hash
2676 * has already been seen in the digits before the point.
2679 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
2682 mem2decimal_from_point (SCM result
, SCM mem
,
2683 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2685 unsigned int idx
= *p_idx
;
2686 enum t_exactness x
= *p_exactness
;
2687 size_t len
= scm_i_string_length (mem
);
2692 if (scm_i_string_ref (mem
, idx
) == '.')
2694 scm_t_bits shift
= 1;
2696 unsigned int digit_value
;
2697 SCM big_shift
= SCM_I_MAKINUM (1);
2702 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2703 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2708 digit_value
= DIGIT2UINT (c
);
2719 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2721 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2722 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2724 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2732 add
= add
* 10 + digit_value
;
2738 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2739 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2740 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2743 result
= scm_divide (result
, big_shift
);
2745 /* We've seen a decimal point, thus the value is implicitly inexact. */
2757 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2759 switch (scm_i_string_ref (mem
, idx
))
2771 c
= scm_i_string_ref (mem
, idx
);
2779 c
= scm_i_string_ref (mem
, idx
);
2788 c
= scm_i_string_ref (mem
, idx
);
2793 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2797 exponent
= DIGIT2UINT (c
);
2800 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2801 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2804 if (exponent
<= SCM_MAXEXP
)
2805 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2811 if (exponent
> SCM_MAXEXP
)
2813 size_t exp_len
= idx
- start
;
2814 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
2815 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2816 scm_out_of_range ("string->number", exp_num
);
2819 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2821 result
= scm_product (result
, e
);
2823 result
= scm_divide2real (result
, e
);
2825 /* We've seen an exponent, thus the value is implicitly inexact. */
2843 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2846 mem2ureal (SCM mem
, unsigned int *p_idx
,
2847 unsigned int radix
, enum t_exactness
*p_exactness
)
2849 unsigned int idx
= *p_idx
;
2851 size_t len
= scm_i_string_length (mem
);
2853 /* Start off believing that the number will be exact. This changes
2854 to INEXACT if we see a decimal point or a hash. */
2855 enum t_exactness x
= EXACT
;
2860 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
2866 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
2868 /* Cobble up the fractional part. We might want to set the
2869 NaN's mantissa from it. */
2871 mem2uinteger (mem
, &idx
, 10, &x
);
2876 if (scm_i_string_ref (mem
, idx
) == '.')
2880 else if (idx
+ 1 == len
)
2882 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
2885 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
,
2892 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
2893 if (scm_is_false (uinteger
))
2898 else if (scm_i_string_ref (mem
, idx
) == '/')
2906 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
2907 if (scm_is_false (divisor
))
2910 /* both are int/big here, I assume */
2911 result
= scm_i_make_ratio (uinteger
, divisor
);
2913 else if (radix
== 10)
2915 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
2916 if (scm_is_false (result
))
2925 /* Update *p_exactness if the number just read was inexact. This is
2926 important for complex numbers, so that a complex number is
2927 treated as inexact overall if either its real or imaginary part
2933 /* When returning an inexact zero, make sure it is represented as a
2934 floating point value so that we can change its sign.
2936 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2937 result
= scm_from_double (0.0);
2943 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2946 mem2complex (SCM mem
, unsigned int idx
,
2947 unsigned int radix
, enum t_exactness
*p_exactness
)
2952 size_t len
= scm_i_string_length (mem
);
2957 c
= scm_i_string_ref (mem
, idx
);
2972 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2973 if (scm_is_false (ureal
))
2975 /* input must be either +i or -i */
2980 if (scm_i_string_ref (mem
, idx
) == 'i'
2981 || scm_i_string_ref (mem
, idx
) == 'I')
2987 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2994 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2995 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
3000 c
= scm_i_string_ref (mem
, idx
);
3004 /* either +<ureal>i or -<ureal>i */
3011 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
3014 /* polar input: <real>@<real>. */
3025 c
= scm_i_string_ref (mem
, idx
);
3043 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3044 if (scm_is_false (angle
))
3049 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3050 angle
= scm_difference (angle
, SCM_UNDEFINED
);
3052 result
= scm_make_polar (ureal
, angle
);
3057 /* expecting input matching <real>[+-]<ureal>?i */
3064 int sign
= (c
== '+') ? 1 : -1;
3065 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3067 if (scm_is_false (imag
))
3068 imag
= SCM_I_MAKINUM (sign
);
3069 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
3070 imag
= scm_difference (imag
, SCM_UNDEFINED
);
3074 if (scm_i_string_ref (mem
, idx
) != 'i'
3075 && scm_i_string_ref (mem
, idx
) != 'I')
3082 return scm_make_rectangular (ureal
, imag
);
3091 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
3093 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
3096 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
3098 unsigned int idx
= 0;
3099 unsigned int radix
= NO_RADIX
;
3100 enum t_exactness forced_x
= NO_EXACTNESS
;
3101 enum t_exactness implicit_x
= EXACT
;
3103 size_t len
= scm_i_string_length (mem
);
3105 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3106 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
3108 switch (scm_i_string_ref (mem
, idx
+ 1))
3111 if (radix
!= NO_RADIX
)
3116 if (radix
!= NO_RADIX
)
3121 if (forced_x
!= NO_EXACTNESS
)
3126 if (forced_x
!= NO_EXACTNESS
)
3131 if (radix
!= NO_RADIX
)
3136 if (radix
!= NO_RADIX
)
3146 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3147 if (radix
== NO_RADIX
)
3148 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
3150 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
3152 if (scm_is_false (result
))
3158 if (SCM_INEXACTP (result
))
3159 return scm_inexact_to_exact (result
);
3163 if (SCM_INEXACTP (result
))
3166 return scm_exact_to_inexact (result
);
3169 if (implicit_x
== INEXACT
)
3171 if (SCM_INEXACTP (result
))
3174 return scm_exact_to_inexact (result
);
3182 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
3183 unsigned int default_radix
)
3185 SCM str
= scm_from_locale_stringn (mem
, len
);
3187 return scm_i_string_to_number (str
, default_radix
);
3191 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3192 (SCM string
, SCM radix
),
3193 "Return a number of the maximally precise representation\n"
3194 "expressed by the given @var{string}. @var{radix} must be an\n"
3195 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3196 "is a default radix that may be overridden by an explicit radix\n"
3197 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3198 "supplied, then the default radix is 10. If string is not a\n"
3199 "syntactically valid notation for a number, then\n"
3200 "@code{string->number} returns @code{#f}.")
3201 #define FUNC_NAME s_scm_string_to_number
3205 SCM_VALIDATE_STRING (1, string
);
3207 if (SCM_UNBNDP (radix
))
3210 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3212 answer
= scm_i_string_to_number (string
, base
);
3213 scm_remember_upto_here_1 (string
);
3219 /*** END strs->nums ***/
3223 scm_bigequal (SCM x
, SCM y
)
3225 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3226 scm_remember_upto_here_2 (x
, y
);
3227 return scm_from_bool (0 == result
);
3231 scm_real_equalp (SCM x
, SCM y
)
3233 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3237 scm_complex_equalp (SCM x
, SCM y
)
3239 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3240 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3244 scm_i_fraction_equalp (SCM x
, SCM y
)
3246 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3247 SCM_FRACTION_NUMERATOR (y
)))
3248 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3249 SCM_FRACTION_DENOMINATOR (y
))))
3256 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3258 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3260 #define FUNC_NAME s_scm_number_p
3262 return scm_from_bool (SCM_NUMBERP (x
));
3266 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3268 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3269 "otherwise. Note that the sets of real, rational and integer\n"
3270 "values form subsets of the set of complex numbers, i. e. the\n"
3271 "predicate will also be fulfilled if @var{x} is a real,\n"
3272 "rational or integer number.")
3273 #define FUNC_NAME s_scm_complex_p
3275 /* all numbers are complex. */
3276 return scm_number_p (x
);
3280 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3282 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3283 "otherwise. Note that the set of integer values forms a subset of\n"
3284 "the set of real numbers, i. e. the predicate will also be\n"
3285 "fulfilled if @var{x} is an integer number.")
3286 #define FUNC_NAME s_scm_real_p
3288 /* we can't represent irrational numbers. */
3289 return scm_rational_p (x
);
3293 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3295 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3296 "otherwise. Note that the set of integer values forms a subset of\n"
3297 "the set of rational numbers, i. e. the predicate will also be\n"
3298 "fulfilled if @var{x} is an integer number.")
3299 #define FUNC_NAME s_scm_rational_p
3301 if (SCM_I_INUMP (x
))
3303 else if (SCM_IMP (x
))
3305 else if (SCM_BIGP (x
))
3307 else if (SCM_FRACTIONP (x
))
3309 else if (SCM_REALP (x
))
3310 /* due to their limited precision, all floating point numbers are
3311 rational as well. */
3318 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3320 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3322 #define FUNC_NAME s_scm_integer_p
3325 if (SCM_I_INUMP (x
))
3331 if (!SCM_INEXACTP (x
))
3333 if (SCM_COMPLEXP (x
))
3335 r
= SCM_REAL_VALUE (x
);
3345 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3347 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3349 #define FUNC_NAME s_scm_inexact_p
3351 if (SCM_INEXACTP (x
))
3353 if (SCM_NUMBERP (x
))
3355 SCM_WRONG_TYPE_ARG (1, x
);
3360 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
3361 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
3362 (SCM x
, SCM y
, SCM rest
),
3363 "Return @code{#t} if all parameters are numerically equal.")
3364 #define FUNC_NAME s_scm_i_num_eq_p
3366 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3368 while (!scm_is_null (rest
))
3370 if (scm_is_false (scm_num_eq_p (x
, y
)))
3374 rest
= scm_cdr (rest
);
3376 return scm_num_eq_p (x
, y
);
3380 scm_num_eq_p (SCM x
, SCM y
)
3383 if (SCM_I_INUMP (x
))
3385 scm_t_signed_bits xx
= SCM_I_INUM (x
);
3386 if (SCM_I_INUMP (y
))
3388 scm_t_signed_bits yy
= SCM_I_INUM (y
);
3389 return scm_from_bool (xx
== yy
);
3391 else if (SCM_BIGP (y
))
3393 else if (SCM_REALP (y
))
3395 /* On a 32-bit system an inum fits a double, we can cast the inum
3396 to a double and compare.
3398 But on a 64-bit system an inum is bigger than a double and
3399 casting it to a double (call that dxx) will round. dxx is at
3400 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3401 an integer and fits a long. So we cast yy to a long and
3402 compare with plain xx.
3404 An alternative (for any size system actually) would be to check
3405 yy is an integer (with floor) and is in range of an inum
3406 (compare against appropriate powers of 2) then test
3407 xx==(scm_t_signed_bits)yy. It's just a matter of which
3408 casts/comparisons might be fastest or easiest for the cpu. */
3410 double yy
= SCM_REAL_VALUE (y
);
3411 return scm_from_bool ((double) xx
== yy
3412 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3413 || xx
== (scm_t_signed_bits
) yy
));
3415 else if (SCM_COMPLEXP (y
))
3416 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3417 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3418 else if (SCM_FRACTIONP (y
))
3421 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3423 else if (SCM_BIGP (x
))
3425 if (SCM_I_INUMP (y
))
3427 else if (SCM_BIGP (y
))
3429 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3430 scm_remember_upto_here_2 (x
, y
);
3431 return scm_from_bool (0 == cmp
);
3433 else if (SCM_REALP (y
))
3436 if (isnan (SCM_REAL_VALUE (y
)))
3438 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3439 scm_remember_upto_here_1 (x
);
3440 return scm_from_bool (0 == cmp
);
3442 else if (SCM_COMPLEXP (y
))
3445 if (0.0 != SCM_COMPLEX_IMAG (y
))
3447 if (isnan (SCM_COMPLEX_REAL (y
)))
3449 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3450 scm_remember_upto_here_1 (x
);
3451 return scm_from_bool (0 == cmp
);
3453 else if (SCM_FRACTIONP (y
))
3456 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3458 else if (SCM_REALP (x
))
3460 double xx
= SCM_REAL_VALUE (x
);
3461 if (SCM_I_INUMP (y
))
3463 /* see comments with inum/real above */
3464 scm_t_signed_bits yy
= SCM_I_INUM (y
);
3465 return scm_from_bool (xx
== (double) yy
3466 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3467 || (scm_t_signed_bits
) xx
== yy
));
3469 else if (SCM_BIGP (y
))
3472 if (isnan (SCM_REAL_VALUE (x
)))
3474 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3475 scm_remember_upto_here_1 (y
);
3476 return scm_from_bool (0 == cmp
);
3478 else if (SCM_REALP (y
))
3479 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3480 else if (SCM_COMPLEXP (y
))
3481 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3482 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3483 else if (SCM_FRACTIONP (y
))
3485 double xx
= SCM_REAL_VALUE (x
);
3489 return scm_from_bool (xx
< 0.0);
3490 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3494 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3496 else if (SCM_COMPLEXP (x
))
3498 if (SCM_I_INUMP (y
))
3499 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3500 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3501 else if (SCM_BIGP (y
))
3504 if (0.0 != SCM_COMPLEX_IMAG (x
))
3506 if (isnan (SCM_COMPLEX_REAL (x
)))
3508 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3509 scm_remember_upto_here_1 (y
);
3510 return scm_from_bool (0 == cmp
);
3512 else if (SCM_REALP (y
))
3513 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3514 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3515 else if (SCM_COMPLEXP (y
))
3516 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3517 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3518 else if (SCM_FRACTIONP (y
))
3521 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3523 xx
= SCM_COMPLEX_REAL (x
);
3527 return scm_from_bool (xx
< 0.0);
3528 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3532 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3534 else if (SCM_FRACTIONP (x
))
3536 if (SCM_I_INUMP (y
))
3538 else if (SCM_BIGP (y
))
3540 else if (SCM_REALP (y
))
3542 double yy
= SCM_REAL_VALUE (y
);
3546 return scm_from_bool (0.0 < yy
);
3547 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3550 else if (SCM_COMPLEXP (y
))
3553 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3555 yy
= SCM_COMPLEX_REAL (y
);
3559 return scm_from_bool (0.0 < yy
);
3560 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3563 else if (SCM_FRACTIONP (y
))
3564 return scm_i_fraction_equalp (x
, y
);
3566 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3569 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
3573 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3574 done are good for inums, but for bignums an answer can almost always be
3575 had by just examining a few high bits of the operands, as done by GMP in
3576 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3577 of the float exponent to take into account. */
3579 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
3580 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
3581 (SCM x
, SCM y
, SCM rest
),
3582 "Return @code{#t} if the list of parameters is monotonically\n"
3584 #define FUNC_NAME s_scm_i_num_less_p
3586 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3588 while (!scm_is_null (rest
))
3590 if (scm_is_false (scm_less_p (x
, y
)))
3594 rest
= scm_cdr (rest
);
3596 return scm_less_p (x
, y
);
3600 scm_less_p (SCM x
, SCM y
)
3603 if (SCM_I_INUMP (x
))
3605 scm_t_inum xx
= SCM_I_INUM (x
);
3606 if (SCM_I_INUMP (y
))
3608 scm_t_inum yy
= SCM_I_INUM (y
);
3609 return scm_from_bool (xx
< yy
);
3611 else if (SCM_BIGP (y
))
3613 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3614 scm_remember_upto_here_1 (y
);
3615 return scm_from_bool (sgn
> 0);
3617 else if (SCM_REALP (y
))
3618 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3619 else if (SCM_FRACTIONP (y
))
3621 /* "x < a/b" becomes "x*b < a" */
3623 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3624 y
= SCM_FRACTION_NUMERATOR (y
);
3628 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3630 else if (SCM_BIGP (x
))
3632 if (SCM_I_INUMP (y
))
3634 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3635 scm_remember_upto_here_1 (x
);
3636 return scm_from_bool (sgn
< 0);
3638 else if (SCM_BIGP (y
))
3640 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3641 scm_remember_upto_here_2 (x
, y
);
3642 return scm_from_bool (cmp
< 0);
3644 else if (SCM_REALP (y
))
3647 if (isnan (SCM_REAL_VALUE (y
)))
3649 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3650 scm_remember_upto_here_1 (x
);
3651 return scm_from_bool (cmp
< 0);
3653 else if (SCM_FRACTIONP (y
))
3656 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3658 else if (SCM_REALP (x
))
3660 if (SCM_I_INUMP (y
))
3661 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3662 else if (SCM_BIGP (y
))
3665 if (isnan (SCM_REAL_VALUE (x
)))
3667 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3668 scm_remember_upto_here_1 (y
);
3669 return scm_from_bool (cmp
> 0);
3671 else if (SCM_REALP (y
))
3672 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3673 else if (SCM_FRACTIONP (y
))
3675 double xx
= SCM_REAL_VALUE (x
);
3679 return scm_from_bool (xx
< 0.0);
3680 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3684 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3686 else if (SCM_FRACTIONP (x
))
3688 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3690 /* "a/b < y" becomes "a < y*b" */
3691 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3692 x
= SCM_FRACTION_NUMERATOR (x
);
3695 else if (SCM_REALP (y
))
3697 double yy
= SCM_REAL_VALUE (y
);
3701 return scm_from_bool (0.0 < yy
);
3702 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3705 else if (SCM_FRACTIONP (y
))
3707 /* "a/b < c/d" becomes "a*d < c*b" */
3708 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3709 SCM_FRACTION_DENOMINATOR (y
));
3710 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3711 SCM_FRACTION_DENOMINATOR (x
));
3717 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3720 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
3724 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
3725 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
3726 (SCM x
, SCM y
, SCM rest
),
3727 "Return @code{#t} if the list of parameters is monotonically\n"
3729 #define FUNC_NAME s_scm_i_num_gr_p
3731 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3733 while (!scm_is_null (rest
))
3735 if (scm_is_false (scm_gr_p (x
, y
)))
3739 rest
= scm_cdr (rest
);
3741 return scm_gr_p (x
, y
);
3744 #define FUNC_NAME s_scm_i_num_gr_p
3746 scm_gr_p (SCM x
, SCM y
)
3748 if (!SCM_NUMBERP (x
))
3749 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3750 else if (!SCM_NUMBERP (y
))
3751 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3753 return scm_less_p (y
, x
);
3758 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
3759 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
3760 (SCM x
, SCM y
, SCM rest
),
3761 "Return @code{#t} if the list of parameters is monotonically\n"
3763 #define FUNC_NAME s_scm_i_num_leq_p
3765 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3767 while (!scm_is_null (rest
))
3769 if (scm_is_false (scm_leq_p (x
, y
)))
3773 rest
= scm_cdr (rest
);
3775 return scm_leq_p (x
, y
);
3778 #define FUNC_NAME s_scm_i_num_leq_p
3780 scm_leq_p (SCM x
, SCM y
)
3782 if (!SCM_NUMBERP (x
))
3783 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3784 else if (!SCM_NUMBERP (y
))
3785 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3786 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3789 return scm_not (scm_less_p (y
, x
));
3794 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
3795 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
3796 (SCM x
, SCM y
, SCM rest
),
3797 "Return @code{#t} if the list of parameters is monotonically\n"
3799 #define FUNC_NAME s_scm_i_num_geq_p
3801 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3803 while (!scm_is_null (rest
))
3805 if (scm_is_false (scm_geq_p (x
, y
)))
3809 rest
= scm_cdr (rest
);
3811 return scm_geq_p (x
, y
);
3814 #define FUNC_NAME s_scm_i_num_geq_p
3816 scm_geq_p (SCM x
, SCM y
)
3818 if (!SCM_NUMBERP (x
))
3819 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3820 else if (!SCM_NUMBERP (y
))
3821 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3822 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3825 return scm_not (scm_less_p (x
, y
));
3830 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3831 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3837 if (SCM_I_INUMP (z
))
3838 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3839 else if (SCM_BIGP (z
))
3841 else if (SCM_REALP (z
))
3842 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3843 else if (SCM_COMPLEXP (z
))
3844 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3845 && SCM_COMPLEX_IMAG (z
) == 0.0);
3846 else if (SCM_FRACTIONP (z
))
3849 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3853 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3854 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3858 scm_positive_p (SCM x
)
3860 if (SCM_I_INUMP (x
))
3861 return scm_from_bool (SCM_I_INUM (x
) > 0);
3862 else if (SCM_BIGP (x
))
3864 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3865 scm_remember_upto_here_1 (x
);
3866 return scm_from_bool (sgn
> 0);
3868 else if (SCM_REALP (x
))
3869 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3870 else if (SCM_FRACTIONP (x
))
3871 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3873 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3877 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3878 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3882 scm_negative_p (SCM x
)
3884 if (SCM_I_INUMP (x
))
3885 return scm_from_bool (SCM_I_INUM (x
) < 0);
3886 else if (SCM_BIGP (x
))
3888 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3889 scm_remember_upto_here_1 (x
);
3890 return scm_from_bool (sgn
< 0);
3892 else if (SCM_REALP (x
))
3893 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3894 else if (SCM_FRACTIONP (x
))
3895 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3897 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3901 /* scm_min and scm_max return an inexact when either argument is inexact, as
3902 required by r5rs. On that basis, for exact/inexact combinations the
3903 exact is converted to inexact to compare and possibly return. This is
3904 unlike scm_less_p above which takes some trouble to preserve all bits in
3905 its test, such trouble is not required for min and max. */
3907 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
3908 (SCM x
, SCM y
, SCM rest
),
3909 "Return the maximum of all parameter values.")
3910 #define FUNC_NAME s_scm_i_max
3912 while (!scm_is_null (rest
))
3913 { x
= scm_max (x
, y
);
3915 rest
= scm_cdr (rest
);
3917 return scm_max (x
, y
);
3921 #define s_max s_scm_i_max
3922 #define g_max g_scm_i_max
3925 scm_max (SCM x
, SCM y
)
3930 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3931 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3934 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3937 if (SCM_I_INUMP (x
))
3939 scm_t_inum xx
= SCM_I_INUM (x
);
3940 if (SCM_I_INUMP (y
))
3942 scm_t_inum yy
= SCM_I_INUM (y
);
3943 return (xx
< yy
) ? y
: x
;
3945 else if (SCM_BIGP (y
))
3947 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3948 scm_remember_upto_here_1 (y
);
3949 return (sgn
< 0) ? x
: y
;
3951 else if (SCM_REALP (y
))
3954 /* if y==NaN then ">" is false and we return NaN */
3955 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3957 else if (SCM_FRACTIONP (y
))
3960 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3963 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3965 else if (SCM_BIGP (x
))
3967 if (SCM_I_INUMP (y
))
3969 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3970 scm_remember_upto_here_1 (x
);
3971 return (sgn
< 0) ? y
: x
;
3973 else if (SCM_BIGP (y
))
3975 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3976 scm_remember_upto_here_2 (x
, y
);
3977 return (cmp
> 0) ? x
: y
;
3979 else if (SCM_REALP (y
))
3981 /* if y==NaN then xx>yy is false, so we return the NaN y */
3984 xx
= scm_i_big2dbl (x
);
3985 yy
= SCM_REAL_VALUE (y
);
3986 return (xx
> yy
? scm_from_double (xx
) : y
);
3988 else if (SCM_FRACTIONP (y
))
3993 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3995 else if (SCM_REALP (x
))
3997 if (SCM_I_INUMP (y
))
3999 double z
= SCM_I_INUM (y
);
4000 /* if x==NaN then "<" is false and we return NaN */
4001 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
4003 else if (SCM_BIGP (y
))
4008 else if (SCM_REALP (y
))
4010 /* if x==NaN then our explicit check means we return NaN
4011 if y==NaN then ">" is false and we return NaN
4012 calling isnan is unavoidable, since it's the only way to know
4013 which of x or y causes any compares to be false */
4014 double xx
= SCM_REAL_VALUE (x
);
4015 return (isnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
4017 else if (SCM_FRACTIONP (y
))
4019 double yy
= scm_i_fraction2double (y
);
4020 double xx
= SCM_REAL_VALUE (x
);
4021 return (xx
< yy
) ? scm_from_double (yy
) : x
;
4024 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4026 else if (SCM_FRACTIONP (x
))
4028 if (SCM_I_INUMP (y
))
4032 else if (SCM_BIGP (y
))
4036 else if (SCM_REALP (y
))
4038 double xx
= scm_i_fraction2double (x
);
4039 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
4041 else if (SCM_FRACTIONP (y
))
4046 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4049 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
4053 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
4054 (SCM x
, SCM y
, SCM rest
),
4055 "Return the minimum of all parameter values.")
4056 #define FUNC_NAME s_scm_i_min
4058 while (!scm_is_null (rest
))
4059 { x
= scm_min (x
, y
);
4061 rest
= scm_cdr (rest
);
4063 return scm_min (x
, y
);
4067 #define s_min s_scm_i_min
4068 #define g_min g_scm_i_min
4071 scm_min (SCM x
, SCM y
)
4076 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
4077 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
4080 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
4083 if (SCM_I_INUMP (x
))
4085 scm_t_inum xx
= SCM_I_INUM (x
);
4086 if (SCM_I_INUMP (y
))
4088 scm_t_inum yy
= SCM_I_INUM (y
);
4089 return (xx
< yy
) ? x
: y
;
4091 else if (SCM_BIGP (y
))
4093 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4094 scm_remember_upto_here_1 (y
);
4095 return (sgn
< 0) ? y
: x
;
4097 else if (SCM_REALP (y
))
4100 /* if y==NaN then "<" is false and we return NaN */
4101 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
4103 else if (SCM_FRACTIONP (y
))
4106 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
4109 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4111 else if (SCM_BIGP (x
))
4113 if (SCM_I_INUMP (y
))
4115 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4116 scm_remember_upto_here_1 (x
);
4117 return (sgn
< 0) ? x
: y
;
4119 else if (SCM_BIGP (y
))
4121 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4122 scm_remember_upto_here_2 (x
, y
);
4123 return (cmp
> 0) ? y
: x
;
4125 else if (SCM_REALP (y
))
4127 /* if y==NaN then xx<yy is false, so we return the NaN y */
4130 xx
= scm_i_big2dbl (x
);
4131 yy
= SCM_REAL_VALUE (y
);
4132 return (xx
< yy
? scm_from_double (xx
) : y
);
4134 else if (SCM_FRACTIONP (y
))
4139 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4141 else if (SCM_REALP (x
))
4143 if (SCM_I_INUMP (y
))
4145 double z
= SCM_I_INUM (y
);
4146 /* if x==NaN then "<" is false and we return NaN */
4147 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
4149 else if (SCM_BIGP (y
))
4154 else if (SCM_REALP (y
))
4156 /* if x==NaN then our explicit check means we return NaN
4157 if y==NaN then "<" is false and we return NaN
4158 calling isnan is unavoidable, since it's the only way to know
4159 which of x or y causes any compares to be false */
4160 double xx
= SCM_REAL_VALUE (x
);
4161 return (isnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
4163 else if (SCM_FRACTIONP (y
))
4165 double yy
= scm_i_fraction2double (y
);
4166 double xx
= SCM_REAL_VALUE (x
);
4167 return (yy
< xx
) ? scm_from_double (yy
) : x
;
4170 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4172 else if (SCM_FRACTIONP (x
))
4174 if (SCM_I_INUMP (y
))
4178 else if (SCM_BIGP (y
))
4182 else if (SCM_REALP (y
))
4184 double xx
= scm_i_fraction2double (x
);
4185 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
4187 else if (SCM_FRACTIONP (y
))
4192 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4195 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
4199 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
4200 (SCM x
, SCM y
, SCM rest
),
4201 "Return the sum of all parameter values. Return 0 if called without\n"
4203 #define FUNC_NAME s_scm_i_sum
4205 while (!scm_is_null (rest
))
4206 { x
= scm_sum (x
, y
);
4208 rest
= scm_cdr (rest
);
4210 return scm_sum (x
, y
);
4214 #define s_sum s_scm_i_sum
4215 #define g_sum g_scm_i_sum
4218 scm_sum (SCM x
, SCM y
)
4220 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4222 if (SCM_NUMBERP (x
)) return x
;
4223 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4224 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4227 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4229 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4231 scm_t_inum xx
= SCM_I_INUM (x
);
4232 scm_t_inum yy
= SCM_I_INUM (y
);
4233 scm_t_inum z
= xx
+ yy
;
4234 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_inum2big (z
);
4236 else if (SCM_BIGP (y
))
4241 else if (SCM_REALP (y
))
4243 scm_t_inum xx
= SCM_I_INUM (x
);
4244 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4246 else if (SCM_COMPLEXP (y
))
4248 scm_t_inum xx
= SCM_I_INUM (x
);
4249 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4250 SCM_COMPLEX_IMAG (y
));
4252 else if (SCM_FRACTIONP (y
))
4253 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4254 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4255 SCM_FRACTION_DENOMINATOR (y
));
4257 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4258 } else if (SCM_BIGP (x
))
4260 if (SCM_I_INUMP (y
))
4265 inum
= SCM_I_INUM (y
);
4268 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4271 SCM result
= scm_i_mkbig ();
4272 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4273 scm_remember_upto_here_1 (x
);
4274 /* we know the result will have to be a bignum */
4277 return scm_i_normbig (result
);
4281 SCM result
= scm_i_mkbig ();
4282 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4283 scm_remember_upto_here_1 (x
);
4284 /* we know the result will have to be a bignum */
4287 return scm_i_normbig (result
);
4290 else if (SCM_BIGP (y
))
4292 SCM result
= scm_i_mkbig ();
4293 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4294 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4295 mpz_add (SCM_I_BIG_MPZ (result
),
4298 scm_remember_upto_here_2 (x
, y
);
4299 /* we know the result will have to be a bignum */
4302 return scm_i_normbig (result
);
4304 else if (SCM_REALP (y
))
4306 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4307 scm_remember_upto_here_1 (x
);
4308 return scm_from_double (result
);
4310 else if (SCM_COMPLEXP (y
))
4312 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4313 + SCM_COMPLEX_REAL (y
));
4314 scm_remember_upto_here_1 (x
);
4315 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4317 else if (SCM_FRACTIONP (y
))
4318 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4319 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4320 SCM_FRACTION_DENOMINATOR (y
));
4322 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4324 else if (SCM_REALP (x
))
4326 if (SCM_I_INUMP (y
))
4327 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4328 else if (SCM_BIGP (y
))
4330 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4331 scm_remember_upto_here_1 (y
);
4332 return scm_from_double (result
);
4334 else if (SCM_REALP (y
))
4335 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4336 else if (SCM_COMPLEXP (y
))
4337 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4338 SCM_COMPLEX_IMAG (y
));
4339 else if (SCM_FRACTIONP (y
))
4340 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4342 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4344 else if (SCM_COMPLEXP (x
))
4346 if (SCM_I_INUMP (y
))
4347 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4348 SCM_COMPLEX_IMAG (x
));
4349 else if (SCM_BIGP (y
))
4351 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4352 + SCM_COMPLEX_REAL (x
));
4353 scm_remember_upto_here_1 (y
);
4354 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4356 else if (SCM_REALP (y
))
4357 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4358 SCM_COMPLEX_IMAG (x
));
4359 else if (SCM_COMPLEXP (y
))
4360 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4361 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4362 else if (SCM_FRACTIONP (y
))
4363 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4364 SCM_COMPLEX_IMAG (x
));
4366 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4368 else if (SCM_FRACTIONP (x
))
4370 if (SCM_I_INUMP (y
))
4371 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4372 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4373 SCM_FRACTION_DENOMINATOR (x
));
4374 else if (SCM_BIGP (y
))
4375 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4376 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4377 SCM_FRACTION_DENOMINATOR (x
));
4378 else if (SCM_REALP (y
))
4379 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4380 else if (SCM_COMPLEXP (y
))
4381 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4382 SCM_COMPLEX_IMAG (y
));
4383 else if (SCM_FRACTIONP (y
))
4384 /* a/b + c/d = (ad + bc) / bd */
4385 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4386 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4387 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4389 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4392 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4396 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4398 "Return @math{@var{x}+1}.")
4399 #define FUNC_NAME s_scm_oneplus
4401 return scm_sum (x
, SCM_I_MAKINUM (1));
4406 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
4407 (SCM x
, SCM y
, SCM rest
),
4408 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
4409 "the sum of all but the first argument are subtracted from the first\n"
4411 #define FUNC_NAME s_scm_i_difference
4413 while (!scm_is_null (rest
))
4414 { x
= scm_difference (x
, y
);
4416 rest
= scm_cdr (rest
);
4418 return scm_difference (x
, y
);
4422 #define s_difference s_scm_i_difference
4423 #define g_difference g_scm_i_difference
4426 scm_difference (SCM x
, SCM y
)
4427 #define FUNC_NAME s_difference
4429 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4432 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4434 if (SCM_I_INUMP (x
))
4436 scm_t_inum xx
= -SCM_I_INUM (x
);
4437 if (SCM_FIXABLE (xx
))
4438 return SCM_I_MAKINUM (xx
);
4440 return scm_i_inum2big (xx
);
4442 else if (SCM_BIGP (x
))
4443 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4444 bignum, but negating that gives a fixnum. */
4445 return scm_i_normbig (scm_i_clonebig (x
, 0));
4446 else if (SCM_REALP (x
))
4447 return scm_from_double (-SCM_REAL_VALUE (x
));
4448 else if (SCM_COMPLEXP (x
))
4449 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4450 -SCM_COMPLEX_IMAG (x
));
4451 else if (SCM_FRACTIONP (x
))
4452 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4453 SCM_FRACTION_DENOMINATOR (x
));
4455 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4458 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4460 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4462 scm_t_inum xx
= SCM_I_INUM (x
);
4463 scm_t_inum yy
= SCM_I_INUM (y
);
4464 scm_t_inum z
= xx
- yy
;
4465 if (SCM_FIXABLE (z
))
4466 return SCM_I_MAKINUM (z
);
4468 return scm_i_inum2big (z
);
4470 else if (SCM_BIGP (y
))
4472 /* inum-x - big-y */
4473 scm_t_inum xx
= SCM_I_INUM (x
);
4476 return scm_i_clonebig (y
, 0);
4479 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4480 SCM result
= scm_i_mkbig ();
4483 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4486 /* x - y == -(y + -x) */
4487 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4488 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4490 scm_remember_upto_here_1 (y
);
4492 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4493 /* we know the result will have to be a bignum */
4496 return scm_i_normbig (result
);
4499 else if (SCM_REALP (y
))
4501 scm_t_inum xx
= SCM_I_INUM (x
);
4502 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4504 else if (SCM_COMPLEXP (y
))
4506 scm_t_inum xx
= SCM_I_INUM (x
);
4507 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4508 - SCM_COMPLEX_IMAG (y
));
4510 else if (SCM_FRACTIONP (y
))
4511 /* a - b/c = (ac - b) / c */
4512 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4513 SCM_FRACTION_NUMERATOR (y
)),
4514 SCM_FRACTION_DENOMINATOR (y
));
4516 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4518 else if (SCM_BIGP (x
))
4520 if (SCM_I_INUMP (y
))
4522 /* big-x - inum-y */
4523 scm_t_inum yy
= SCM_I_INUM (y
);
4524 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4526 scm_remember_upto_here_1 (x
);
4528 return (SCM_FIXABLE (-yy
) ?
4529 SCM_I_MAKINUM (-yy
) : scm_from_inum (-yy
));
4532 SCM result
= scm_i_mkbig ();
4535 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4537 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4538 scm_remember_upto_here_1 (x
);
4540 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4541 /* we know the result will have to be a bignum */
4544 return scm_i_normbig (result
);
4547 else if (SCM_BIGP (y
))
4549 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4550 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4551 SCM result
= scm_i_mkbig ();
4552 mpz_sub (SCM_I_BIG_MPZ (result
),
4555 scm_remember_upto_here_2 (x
, y
);
4556 /* we know the result will have to be a bignum */
4557 if ((sgn_x
== 1) && (sgn_y
== -1))
4559 if ((sgn_x
== -1) && (sgn_y
== 1))
4561 return scm_i_normbig (result
);
4563 else if (SCM_REALP (y
))
4565 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4566 scm_remember_upto_here_1 (x
);
4567 return scm_from_double (result
);
4569 else if (SCM_COMPLEXP (y
))
4571 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4572 - SCM_COMPLEX_REAL (y
));
4573 scm_remember_upto_here_1 (x
);
4574 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4576 else if (SCM_FRACTIONP (y
))
4577 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4578 SCM_FRACTION_NUMERATOR (y
)),
4579 SCM_FRACTION_DENOMINATOR (y
));
4580 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4582 else if (SCM_REALP (x
))
4584 if (SCM_I_INUMP (y
))
4585 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4586 else if (SCM_BIGP (y
))
4588 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4589 scm_remember_upto_here_1 (x
);
4590 return scm_from_double (result
);
4592 else if (SCM_REALP (y
))
4593 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4594 else if (SCM_COMPLEXP (y
))
4595 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4596 -SCM_COMPLEX_IMAG (y
));
4597 else if (SCM_FRACTIONP (y
))
4598 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4600 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4602 else if (SCM_COMPLEXP (x
))
4604 if (SCM_I_INUMP (y
))
4605 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4606 SCM_COMPLEX_IMAG (x
));
4607 else if (SCM_BIGP (y
))
4609 double real_part
= (SCM_COMPLEX_REAL (x
)
4610 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4611 scm_remember_upto_here_1 (x
);
4612 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4614 else if (SCM_REALP (y
))
4615 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4616 SCM_COMPLEX_IMAG (x
));
4617 else if (SCM_COMPLEXP (y
))
4618 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4619 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4620 else if (SCM_FRACTIONP (y
))
4621 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4622 SCM_COMPLEX_IMAG (x
));
4624 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4626 else if (SCM_FRACTIONP (x
))
4628 if (SCM_I_INUMP (y
))
4629 /* a/b - c = (a - cb) / b */
4630 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4631 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4632 SCM_FRACTION_DENOMINATOR (x
));
4633 else if (SCM_BIGP (y
))
4634 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4635 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4636 SCM_FRACTION_DENOMINATOR (x
));
4637 else if (SCM_REALP (y
))
4638 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4639 else if (SCM_COMPLEXP (y
))
4640 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4641 -SCM_COMPLEX_IMAG (y
));
4642 else if (SCM_FRACTIONP (y
))
4643 /* a/b - c/d = (ad - bc) / bd */
4644 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4645 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4646 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4648 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4651 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4656 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4658 "Return @math{@var{x}-1}.")
4659 #define FUNC_NAME s_scm_oneminus
4661 return scm_difference (x
, SCM_I_MAKINUM (1));
4666 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
4667 (SCM x
, SCM y
, SCM rest
),
4668 "Return the product of all arguments. If called without arguments,\n"
4670 #define FUNC_NAME s_scm_i_product
4672 while (!scm_is_null (rest
))
4673 { x
= scm_product (x
, y
);
4675 rest
= scm_cdr (rest
);
4677 return scm_product (x
, y
);
4681 #define s_product s_scm_i_product
4682 #define g_product g_scm_i_product
4685 scm_product (SCM x
, SCM y
)
4687 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4690 return SCM_I_MAKINUM (1L);
4691 else if (SCM_NUMBERP (x
))
4694 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4697 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4702 xx
= SCM_I_INUM (x
);
4706 case 0: return x
; break;
4707 case 1: return y
; break;
4710 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4712 scm_t_inum yy
= SCM_I_INUM (y
);
4713 scm_t_inum kk
= xx
* yy
;
4714 SCM k
= SCM_I_MAKINUM (kk
);
4715 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4719 SCM result
= scm_i_inum2big (xx
);
4720 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4721 return scm_i_normbig (result
);
4724 else if (SCM_BIGP (y
))
4726 SCM result
= scm_i_mkbig ();
4727 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4728 scm_remember_upto_here_1 (y
);
4731 else if (SCM_REALP (y
))
4732 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4733 else if (SCM_COMPLEXP (y
))
4734 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4735 xx
* SCM_COMPLEX_IMAG (y
));
4736 else if (SCM_FRACTIONP (y
))
4737 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4738 SCM_FRACTION_DENOMINATOR (y
));
4740 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4742 else if (SCM_BIGP (x
))
4744 if (SCM_I_INUMP (y
))
4749 else if (SCM_BIGP (y
))
4751 SCM result
= scm_i_mkbig ();
4752 mpz_mul (SCM_I_BIG_MPZ (result
),
4755 scm_remember_upto_here_2 (x
, y
);
4758 else if (SCM_REALP (y
))
4760 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4761 scm_remember_upto_here_1 (x
);
4762 return scm_from_double (result
);
4764 else if (SCM_COMPLEXP (y
))
4766 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4767 scm_remember_upto_here_1 (x
);
4768 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4769 z
* SCM_COMPLEX_IMAG (y
));
4771 else if (SCM_FRACTIONP (y
))
4772 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4773 SCM_FRACTION_DENOMINATOR (y
));
4775 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4777 else if (SCM_REALP (x
))
4779 if (SCM_I_INUMP (y
))
4781 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4782 if (scm_is_eq (y
, SCM_INUM0
))
4784 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4786 else if (SCM_BIGP (y
))
4788 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4789 scm_remember_upto_here_1 (y
);
4790 return scm_from_double (result
);
4792 else if (SCM_REALP (y
))
4793 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4794 else if (SCM_COMPLEXP (y
))
4795 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4796 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4797 else if (SCM_FRACTIONP (y
))
4798 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4800 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4802 else if (SCM_COMPLEXP (x
))
4804 if (SCM_I_INUMP (y
))
4806 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4807 if (scm_is_eq (y
, SCM_INUM0
))
4809 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4810 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4812 else if (SCM_BIGP (y
))
4814 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4815 scm_remember_upto_here_1 (y
);
4816 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4817 z
* SCM_COMPLEX_IMAG (x
));
4819 else if (SCM_REALP (y
))
4820 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4821 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4822 else if (SCM_COMPLEXP (y
))
4824 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4825 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4826 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4827 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4829 else if (SCM_FRACTIONP (y
))
4831 double yy
= scm_i_fraction2double (y
);
4832 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4833 yy
* SCM_COMPLEX_IMAG (x
));
4836 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4838 else if (SCM_FRACTIONP (x
))
4840 if (SCM_I_INUMP (y
))
4841 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4842 SCM_FRACTION_DENOMINATOR (x
));
4843 else if (SCM_BIGP (y
))
4844 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4845 SCM_FRACTION_DENOMINATOR (x
));
4846 else if (SCM_REALP (y
))
4847 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4848 else if (SCM_COMPLEXP (y
))
4850 double xx
= scm_i_fraction2double (x
);
4851 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4852 xx
* SCM_COMPLEX_IMAG (y
));
4854 else if (SCM_FRACTIONP (y
))
4855 /* a/b * c/d = ac / bd */
4856 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4857 SCM_FRACTION_NUMERATOR (y
)),
4858 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4859 SCM_FRACTION_DENOMINATOR (y
)));
4861 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4864 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4867 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4868 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4869 #define ALLOW_DIVIDE_BY_ZERO
4870 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4873 /* The code below for complex division is adapted from the GNU
4874 libstdc++, which adapted it from f2c's libF77, and is subject to
4877 /****************************************************************
4878 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4880 Permission to use, copy, modify, and distribute this software
4881 and its documentation for any purpose and without fee is hereby
4882 granted, provided that the above copyright notice appear in all
4883 copies and that both that the copyright notice and this
4884 permission notice and warranty disclaimer appear in supporting
4885 documentation, and that the names of AT&T Bell Laboratories or
4886 Bellcore or any of their entities not be used in advertising or
4887 publicity pertaining to distribution of the software without
4888 specific, written prior permission.
4890 AT&T and Bellcore disclaim all warranties with regard to this
4891 software, including all implied warranties of merchantability
4892 and fitness. In no event shall AT&T or Bellcore be liable for
4893 any special, indirect or consequential damages or any damages
4894 whatsoever resulting from loss of use, data or profits, whether
4895 in an action of contract, negligence or other tortious action,
4896 arising out of or in connection with the use or performance of
4898 ****************************************************************/
4900 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
4901 (SCM x
, SCM y
, SCM rest
),
4902 "Divide the first argument by the product of the remaining\n"
4903 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
4905 #define FUNC_NAME s_scm_i_divide
4907 while (!scm_is_null (rest
))
4908 { x
= scm_divide (x
, y
);
4910 rest
= scm_cdr (rest
);
4912 return scm_divide (x
, y
);
4916 #define s_divide s_scm_i_divide
4917 #define g_divide g_scm_i_divide
4920 do_divide (SCM x
, SCM y
, int inexact
)
4921 #define FUNC_NAME s_divide
4925 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4928 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4929 else if (SCM_I_INUMP (x
))
4931 scm_t_inum xx
= SCM_I_INUM (x
);
4932 if (xx
== 1 || xx
== -1)
4934 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4936 scm_num_overflow (s_divide
);
4941 return scm_from_double (1.0 / (double) xx
);
4942 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4945 else if (SCM_BIGP (x
))
4948 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4949 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4951 else if (SCM_REALP (x
))
4953 double xx
= SCM_REAL_VALUE (x
);
4954 #ifndef ALLOW_DIVIDE_BY_ZERO
4956 scm_num_overflow (s_divide
);
4959 return scm_from_double (1.0 / xx
);
4961 else if (SCM_COMPLEXP (x
))
4963 double r
= SCM_COMPLEX_REAL (x
);
4964 double i
= SCM_COMPLEX_IMAG (x
);
4965 if (fabs(r
) <= fabs(i
))
4968 double d
= i
* (1.0 + t
* t
);
4969 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4974 double d
= r
* (1.0 + t
* t
);
4975 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4978 else if (SCM_FRACTIONP (x
))
4979 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4980 SCM_FRACTION_NUMERATOR (x
));
4982 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4985 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4987 scm_t_inum xx
= SCM_I_INUM (x
);
4988 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4990 scm_t_inum yy
= SCM_I_INUM (y
);
4993 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4994 scm_num_overflow (s_divide
);
4996 return scm_from_double ((double) xx
/ (double) yy
);
4999 else if (xx
% yy
!= 0)
5002 return scm_from_double ((double) xx
/ (double) yy
);
5003 else return scm_i_make_ratio (x
, y
);
5007 scm_t_inum z
= xx
/ yy
;
5008 if (SCM_FIXABLE (z
))
5009 return SCM_I_MAKINUM (z
);
5011 return scm_i_inum2big (z
);
5014 else if (SCM_BIGP (y
))
5017 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
5018 else return scm_i_make_ratio (x
, y
);
5020 else if (SCM_REALP (y
))
5022 double yy
= SCM_REAL_VALUE (y
);
5023 #ifndef ALLOW_DIVIDE_BY_ZERO
5025 scm_num_overflow (s_divide
);
5028 return scm_from_double ((double) xx
/ yy
);
5030 else if (SCM_COMPLEXP (y
))
5033 complex_div
: /* y _must_ be a complex number */
5035 double r
= SCM_COMPLEX_REAL (y
);
5036 double i
= SCM_COMPLEX_IMAG (y
);
5037 if (fabs(r
) <= fabs(i
))
5040 double d
= i
* (1.0 + t
* t
);
5041 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
5046 double d
= r
* (1.0 + t
* t
);
5047 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
5051 else if (SCM_FRACTIONP (y
))
5052 /* a / b/c = ac / b */
5053 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5054 SCM_FRACTION_NUMERATOR (y
));
5056 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5058 else if (SCM_BIGP (x
))
5060 if (SCM_I_INUMP (y
))
5062 scm_t_inum yy
= SCM_I_INUM (y
);
5065 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5066 scm_num_overflow (s_divide
);
5068 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5069 scm_remember_upto_here_1 (x
);
5070 return (sgn
== 0) ? scm_nan () : scm_inf ();
5077 /* FIXME: HMM, what are the relative performance issues here?
5078 We need to test. Is it faster on average to test
5079 divisible_p, then perform whichever operation, or is it
5080 faster to perform the integer div opportunistically and
5081 switch to real if there's a remainder? For now we take the
5082 middle ground: test, then if divisible, use the faster div
5085 scm_t_inum abs_yy
= yy
< 0 ? -yy
: yy
;
5086 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
5090 SCM result
= scm_i_mkbig ();
5091 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
5092 scm_remember_upto_here_1 (x
);
5094 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5095 return scm_i_normbig (result
);
5100 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
5101 else return scm_i_make_ratio (x
, y
);
5105 else if (SCM_BIGP (y
))
5107 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
5110 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5111 scm_num_overflow (s_divide
);
5113 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5114 scm_remember_upto_here_1 (x
);
5115 return (sgn
== 0) ? scm_nan () : scm_inf ();
5123 /* It's easily possible for the ratio x/y to fit a double
5124 but one or both x and y be too big to fit a double,
5125 hence the use of mpq_get_d rather than converting and
5128 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
5129 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
5130 return scm_from_double (mpq_get_d (q
));
5134 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
5138 SCM result
= scm_i_mkbig ();
5139 mpz_divexact (SCM_I_BIG_MPZ (result
),
5142 scm_remember_upto_here_2 (x
, y
);
5143 return scm_i_normbig (result
);
5146 return scm_i_make_ratio (x
, y
);
5150 else if (SCM_REALP (y
))
5152 double yy
= SCM_REAL_VALUE (y
);
5153 #ifndef ALLOW_DIVIDE_BY_ZERO
5155 scm_num_overflow (s_divide
);
5158 return scm_from_double (scm_i_big2dbl (x
) / yy
);
5160 else if (SCM_COMPLEXP (y
))
5162 a
= scm_i_big2dbl (x
);
5165 else if (SCM_FRACTIONP (y
))
5166 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5167 SCM_FRACTION_NUMERATOR (y
));
5169 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5171 else if (SCM_REALP (x
))
5173 double rx
= SCM_REAL_VALUE (x
);
5174 if (SCM_I_INUMP (y
))
5176 scm_t_inum yy
= SCM_I_INUM (y
);
5177 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5179 scm_num_overflow (s_divide
);
5182 return scm_from_double (rx
/ (double) yy
);
5184 else if (SCM_BIGP (y
))
5186 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5187 scm_remember_upto_here_1 (y
);
5188 return scm_from_double (rx
/ dby
);
5190 else if (SCM_REALP (y
))
5192 double yy
= SCM_REAL_VALUE (y
);
5193 #ifndef ALLOW_DIVIDE_BY_ZERO
5195 scm_num_overflow (s_divide
);
5198 return scm_from_double (rx
/ yy
);
5200 else if (SCM_COMPLEXP (y
))
5205 else if (SCM_FRACTIONP (y
))
5206 return scm_from_double (rx
/ scm_i_fraction2double (y
));
5208 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5210 else if (SCM_COMPLEXP (x
))
5212 double rx
= SCM_COMPLEX_REAL (x
);
5213 double ix
= SCM_COMPLEX_IMAG (x
);
5214 if (SCM_I_INUMP (y
))
5216 scm_t_inum yy
= SCM_I_INUM (y
);
5217 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5219 scm_num_overflow (s_divide
);
5224 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
5227 else if (SCM_BIGP (y
))
5229 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5230 scm_remember_upto_here_1 (y
);
5231 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
5233 else if (SCM_REALP (y
))
5235 double yy
= SCM_REAL_VALUE (y
);
5236 #ifndef ALLOW_DIVIDE_BY_ZERO
5238 scm_num_overflow (s_divide
);
5241 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5243 else if (SCM_COMPLEXP (y
))
5245 double ry
= SCM_COMPLEX_REAL (y
);
5246 double iy
= SCM_COMPLEX_IMAG (y
);
5247 if (fabs(ry
) <= fabs(iy
))
5250 double d
= iy
* (1.0 + t
* t
);
5251 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5256 double d
= ry
* (1.0 + t
* t
);
5257 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5260 else if (SCM_FRACTIONP (y
))
5262 double yy
= scm_i_fraction2double (y
);
5263 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5266 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5268 else if (SCM_FRACTIONP (x
))
5270 if (SCM_I_INUMP (y
))
5272 scm_t_inum yy
= SCM_I_INUM (y
);
5273 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5275 scm_num_overflow (s_divide
);
5278 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5279 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5281 else if (SCM_BIGP (y
))
5283 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5284 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5286 else if (SCM_REALP (y
))
5288 double yy
= SCM_REAL_VALUE (y
);
5289 #ifndef ALLOW_DIVIDE_BY_ZERO
5291 scm_num_overflow (s_divide
);
5294 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5296 else if (SCM_COMPLEXP (y
))
5298 a
= scm_i_fraction2double (x
);
5301 else if (SCM_FRACTIONP (y
))
5302 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5303 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5305 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5308 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5312 scm_divide (SCM x
, SCM y
)
5314 return do_divide (x
, y
, 0);
5317 static SCM
scm_divide2real (SCM x
, SCM y
)
5319 return do_divide (x
, y
, 1);
5325 scm_c_truncate (double x
)
5336 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5337 half-way case (ie. when x is an integer plus 0.5) going upwards.
5338 Then half-way cases are identified and adjusted down if the
5339 round-upwards didn't give the desired even integer.
5341 "plus_half == result" identifies a half-way case. If plus_half, which is
5342 x + 0.5, is an integer then x must be an integer plus 0.5.
5344 An odd "result" value is identified with result/2 != floor(result/2).
5345 This is done with plus_half, since that value is ready for use sooner in
5346 a pipelined cpu, and we're already requiring plus_half == result.
5348 Note however that we need to be careful when x is big and already an
5349 integer. In that case "x+0.5" may round to an adjacent integer, causing
5350 us to return such a value, incorrectly. For instance if the hardware is
5351 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5352 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5353 returned. Or if the hardware is in round-upwards mode, then other bigger
5354 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5355 representable value, 2^128+2^76 (or whatever), again incorrect.
5357 These bad roundings of x+0.5 are avoided by testing at the start whether
5358 x is already an integer. If it is then clearly that's the desired result
5359 already. And if it's not then the exponent must be small enough to allow
5360 an 0.5 to be represented, and hence added without a bad rounding. */
5363 scm_c_round (double x
)
5365 double plus_half
, result
;
5370 plus_half
= x
+ 0.5;
5371 result
= floor (plus_half
);
5372 /* Adjust so that the rounding is towards even. */
5373 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5378 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5380 "Round the number @var{x} towards zero.")
5381 #define FUNC_NAME s_scm_truncate_number
5383 if (scm_is_false (scm_negative_p (x
)))
5384 return scm_floor (x
);
5386 return scm_ceiling (x
);
5390 static SCM exactly_one_half
;
5392 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5394 "Round the number @var{x} towards the nearest integer. "
5395 "When it is exactly halfway between two integers, "
5396 "round towards the even one.")
5397 #define FUNC_NAME s_scm_round_number
5399 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5401 else if (SCM_REALP (x
))
5402 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5405 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5406 single quotient+remainder division then examining to see which way
5407 the rounding should go. */
5408 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5409 SCM result
= scm_floor (plus_half
);
5410 /* Adjust so that the rounding is towards even. */
5411 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5412 && scm_is_true (scm_odd_p (result
)))
5413 return scm_difference (result
, SCM_I_MAKINUM (1));
5420 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5422 "Round the number @var{x} towards minus infinity.")
5423 #define FUNC_NAME s_scm_floor
5425 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5427 else if (SCM_REALP (x
))
5428 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5429 else if (SCM_FRACTIONP (x
))
5431 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5432 SCM_FRACTION_DENOMINATOR (x
));
5433 if (scm_is_false (scm_negative_p (x
)))
5435 /* For positive x, rounding towards zero is correct. */
5440 /* For negative x, we need to return q-1 unless x is an
5441 integer. But fractions are never integer, per our
5443 return scm_difference (q
, SCM_I_MAKINUM (1));
5447 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5451 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5453 "Round the number @var{x} towards infinity.")
5454 #define FUNC_NAME s_scm_ceiling
5456 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5458 else if (SCM_REALP (x
))
5459 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5460 else if (SCM_FRACTIONP (x
))
5462 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5463 SCM_FRACTION_DENOMINATOR (x
));
5464 if (scm_is_false (scm_positive_p (x
)))
5466 /* For negative x, rounding towards zero is correct. */
5471 /* For positive x, we need to return q+1 unless x is an
5472 integer. But fractions are never integer, per our
5474 return scm_sum (q
, SCM_I_MAKINUM (1));
5478 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5482 /* sin/cos/tan/asin/acos/atan
5483 sinh/cosh/tanh/asinh/acosh/atanh
5484 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5485 Written by Jerry D. Hedden, (C) FSF.
5486 See the file `COPYING' for terms applying to this program. */
5488 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5490 "Return @var{x} raised to the power of @var{y}.")
5491 #define FUNC_NAME s_scm_expt
5493 if (scm_is_integer (y
))
5495 if (scm_is_true (scm_exact_p (y
)))
5496 return scm_integer_expt (x
, y
);
5499 /* Here we handle the case where the exponent is an inexact
5500 integer. We make the exponent exact in order to use
5501 scm_integer_expt, and thus avoid the spurious imaginary
5502 parts that may result from round-off errors in the general
5503 e^(y log x) method below (for example when squaring a large
5504 negative number). In this case, we must return an inexact
5505 result for correctness. We also make the base inexact so
5506 that scm_integer_expt will use fast inexact arithmetic
5507 internally. Note that making the base inexact is not
5508 sufficient to guarantee an inexact result, because
5509 scm_integer_expt will return an exact 1 when the exponent
5510 is 0, even if the base is inexact. */
5511 return scm_exact_to_inexact
5512 (scm_integer_expt (scm_exact_to_inexact (x
),
5513 scm_inexact_to_exact (y
)));
5516 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5518 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5521 return scm_exp (scm_product (scm_log (x
), y
));
5525 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5527 "Compute the sine of @var{z}.")
5528 #define FUNC_NAME s_scm_sin
5530 if (scm_is_real (z
))
5531 return scm_from_double (sin (scm_to_double (z
)));
5532 else if (SCM_COMPLEXP (z
))
5534 x
= SCM_COMPLEX_REAL (z
);
5535 y
= SCM_COMPLEX_IMAG (z
);
5536 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5537 cos (x
) * sinh (y
));
5540 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5544 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5546 "Compute the cosine of @var{z}.")
5547 #define FUNC_NAME s_scm_cos
5549 if (scm_is_real (z
))
5550 return scm_from_double (cos (scm_to_double (z
)));
5551 else if (SCM_COMPLEXP (z
))
5553 x
= SCM_COMPLEX_REAL (z
);
5554 y
= SCM_COMPLEX_IMAG (z
);
5555 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5556 -sin (x
) * sinh (y
));
5559 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5563 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5565 "Compute the tangent of @var{z}.")
5566 #define FUNC_NAME s_scm_tan
5568 if (scm_is_real (z
))
5569 return scm_from_double (tan (scm_to_double (z
)));
5570 else if (SCM_COMPLEXP (z
))
5572 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5573 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5574 w
= cos (x
) + cosh (y
);
5575 #ifndef ALLOW_DIVIDE_BY_ZERO
5577 scm_num_overflow (s_scm_tan
);
5579 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5582 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5586 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5588 "Compute the hyperbolic sine of @var{z}.")
5589 #define FUNC_NAME s_scm_sinh
5591 if (scm_is_real (z
))
5592 return scm_from_double (sinh (scm_to_double (z
)));
5593 else if (SCM_COMPLEXP (z
))
5595 x
= SCM_COMPLEX_REAL (z
);
5596 y
= SCM_COMPLEX_IMAG (z
);
5597 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5598 cosh (x
) * sin (y
));
5601 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5605 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5607 "Compute the hyperbolic cosine of @var{z}.")
5608 #define FUNC_NAME s_scm_cosh
5610 if (scm_is_real (z
))
5611 return scm_from_double (cosh (scm_to_double (z
)));
5612 else if (SCM_COMPLEXP (z
))
5614 x
= SCM_COMPLEX_REAL (z
);
5615 y
= SCM_COMPLEX_IMAG (z
);
5616 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5617 sinh (x
) * sin (y
));
5620 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5624 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5626 "Compute the hyperbolic tangent of @var{z}.")
5627 #define FUNC_NAME s_scm_tanh
5629 if (scm_is_real (z
))
5630 return scm_from_double (tanh (scm_to_double (z
)));
5631 else if (SCM_COMPLEXP (z
))
5633 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5634 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5635 w
= cosh (x
) + cos (y
);
5636 #ifndef ALLOW_DIVIDE_BY_ZERO
5638 scm_num_overflow (s_scm_tanh
);
5640 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5643 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5647 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5649 "Compute the arc sine of @var{z}.")
5650 #define FUNC_NAME s_scm_asin
5652 if (scm_is_real (z
))
5654 double w
= scm_to_double (z
);
5655 if (w
>= -1.0 && w
<= 1.0)
5656 return scm_from_double (asin (w
));
5658 return scm_product (scm_c_make_rectangular (0, -1),
5659 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5661 else if (SCM_COMPLEXP (z
))
5663 x
= SCM_COMPLEX_REAL (z
);
5664 y
= SCM_COMPLEX_IMAG (z
);
5665 return scm_product (scm_c_make_rectangular (0, -1),
5666 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5669 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5673 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5675 "Compute the arc cosine of @var{z}.")
5676 #define FUNC_NAME s_scm_acos
5678 if (scm_is_real (z
))
5680 double w
= scm_to_double (z
);
5681 if (w
>= -1.0 && w
<= 1.0)
5682 return scm_from_double (acos (w
));
5684 return scm_sum (scm_from_double (acos (0.0)),
5685 scm_product (scm_c_make_rectangular (0, 1),
5686 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5688 else if (SCM_COMPLEXP (z
))
5690 x
= SCM_COMPLEX_REAL (z
);
5691 y
= SCM_COMPLEX_IMAG (z
);
5692 return scm_sum (scm_from_double (acos (0.0)),
5693 scm_product (scm_c_make_rectangular (0, 1),
5694 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5697 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5701 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5703 "With one argument, compute the arc tangent of @var{z}.\n"
5704 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5705 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5706 #define FUNC_NAME s_scm_atan
5710 if (scm_is_real (z
))
5711 return scm_from_double (atan (scm_to_double (z
)));
5712 else if (SCM_COMPLEXP (z
))
5715 v
= SCM_COMPLEX_REAL (z
);
5716 w
= SCM_COMPLEX_IMAG (z
);
5717 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5718 scm_c_make_rectangular (v
, w
+ 1.0))),
5719 scm_c_make_rectangular (0, 2));
5722 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5724 else if (scm_is_real (z
))
5726 if (scm_is_real (y
))
5727 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5729 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5732 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5736 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5738 "Compute the inverse hyperbolic sine of @var{z}.")
5739 #define FUNC_NAME s_scm_sys_asinh
5741 if (scm_is_real (z
))
5742 return scm_from_double (asinh (scm_to_double (z
)));
5743 else if (scm_is_number (z
))
5744 return scm_log (scm_sum (z
,
5745 scm_sqrt (scm_sum (scm_product (z
, z
),
5746 SCM_I_MAKINUM (1)))));
5748 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5752 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5754 "Compute the inverse hyperbolic cosine of @var{z}.")
5755 #define FUNC_NAME s_scm_sys_acosh
5757 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5758 return scm_from_double (acosh (scm_to_double (z
)));
5759 else if (scm_is_number (z
))
5760 return scm_log (scm_sum (z
,
5761 scm_sqrt (scm_difference (scm_product (z
, z
),
5762 SCM_I_MAKINUM (1)))));
5764 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5768 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5770 "Compute the inverse hyperbolic tangent of @var{z}.")
5771 #define FUNC_NAME s_scm_sys_atanh
5773 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5774 return scm_from_double (atanh (scm_to_double (z
)));
5775 else if (scm_is_number (z
))
5776 return scm_divide (scm_log (scm_divide (scm_sum (SCM_I_MAKINUM (1), z
),
5777 scm_difference (SCM_I_MAKINUM (1), z
))),
5780 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5785 scm_c_make_rectangular (double re
, double im
)
5788 return scm_from_double (re
);
5793 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5795 SCM_SET_CELL_TYPE (z
, scm_tc16_complex
);
5796 SCM_COMPLEX_REAL (z
) = re
;
5797 SCM_COMPLEX_IMAG (z
) = im
;
5802 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5803 (SCM real_part
, SCM imaginary_part
),
5804 "Return a complex number constructed of the given @var{real-part} "
5805 "and @var{imaginary-part} parts.")
5806 #define FUNC_NAME s_scm_make_rectangular
5808 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5809 SCM_ARG1
, FUNC_NAME
, "real");
5810 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5811 SCM_ARG2
, FUNC_NAME
, "real");
5812 return scm_c_make_rectangular (scm_to_double (real_part
),
5813 scm_to_double (imaginary_part
));
5818 scm_c_make_polar (double mag
, double ang
)
5822 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5823 use it on Glibc-based systems that have it (it's a GNU extension). See
5824 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5826 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5827 sincos (ang
, &s
, &c
);
5832 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5835 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5837 "Return the complex number @var{x} * e^(i * @var{y}).")
5838 #define FUNC_NAME s_scm_make_polar
5840 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5841 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5842 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5847 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5848 /* "Return the real part of the number @var{z}."
5851 scm_real_part (SCM z
)
5853 if (SCM_I_INUMP (z
))
5855 else if (SCM_BIGP (z
))
5857 else if (SCM_REALP (z
))
5859 else if (SCM_COMPLEXP (z
))
5860 return scm_from_double (SCM_COMPLEX_REAL (z
));
5861 else if (SCM_FRACTIONP (z
))
5864 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5868 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5869 /* "Return the imaginary part of the number @var{z}."
5872 scm_imag_part (SCM z
)
5874 if (SCM_I_INUMP (z
))
5876 else if (SCM_BIGP (z
))
5878 else if (SCM_REALP (z
))
5880 else if (SCM_COMPLEXP (z
))
5881 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5882 else if (SCM_FRACTIONP (z
))
5885 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5888 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5889 /* "Return the numerator of the number @var{z}."
5892 scm_numerator (SCM z
)
5894 if (SCM_I_INUMP (z
))
5896 else if (SCM_BIGP (z
))
5898 else if (SCM_FRACTIONP (z
))
5899 return SCM_FRACTION_NUMERATOR (z
);
5900 else if (SCM_REALP (z
))
5901 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5903 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5907 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5908 /* "Return the denominator of the number @var{z}."
5911 scm_denominator (SCM z
)
5913 if (SCM_I_INUMP (z
))
5914 return SCM_I_MAKINUM (1);
5915 else if (SCM_BIGP (z
))
5916 return SCM_I_MAKINUM (1);
5917 else if (SCM_FRACTIONP (z
))
5918 return SCM_FRACTION_DENOMINATOR (z
);
5919 else if (SCM_REALP (z
))
5920 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5922 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5925 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5926 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5927 * "@code{abs} for real arguments, but also allows complex numbers."
5930 scm_magnitude (SCM z
)
5932 if (SCM_I_INUMP (z
))
5934 scm_t_inum zz
= SCM_I_INUM (z
);
5937 else if (SCM_POSFIXABLE (-zz
))
5938 return SCM_I_MAKINUM (-zz
);
5940 return scm_i_inum2big (-zz
);
5942 else if (SCM_BIGP (z
))
5944 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5945 scm_remember_upto_here_1 (z
);
5947 return scm_i_clonebig (z
, 0);
5951 else if (SCM_REALP (z
))
5952 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5953 else if (SCM_COMPLEXP (z
))
5954 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5955 else if (SCM_FRACTIONP (z
))
5957 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5959 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5960 SCM_FRACTION_DENOMINATOR (z
));
5963 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5967 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5968 /* "Return the angle of the complex number @var{z}."
5973 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5974 flo0 to save allocating a new flonum with scm_from_double each time.
5975 But if atan2 follows the floating point rounding mode, then the value
5976 is not a constant. Maybe it'd be close enough though. */
5977 if (SCM_I_INUMP (z
))
5979 if (SCM_I_INUM (z
) >= 0)
5982 return scm_from_double (atan2 (0.0, -1.0));
5984 else if (SCM_BIGP (z
))
5986 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5987 scm_remember_upto_here_1 (z
);
5989 return scm_from_double (atan2 (0.0, -1.0));
5993 else if (SCM_REALP (z
))
5995 if (SCM_REAL_VALUE (z
) >= 0)
5998 return scm_from_double (atan2 (0.0, -1.0));
6000 else if (SCM_COMPLEXP (z
))
6001 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
6002 else if (SCM_FRACTIONP (z
))
6004 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
6006 else return scm_from_double (atan2 (0.0, -1.0));
6009 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
6013 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
6014 /* Convert the number @var{x} to its inexact representation.\n"
6017 scm_exact_to_inexact (SCM z
)
6019 if (SCM_I_INUMP (z
))
6020 return scm_from_double ((double) SCM_I_INUM (z
));
6021 else if (SCM_BIGP (z
))
6022 return scm_from_double (scm_i_big2dbl (z
));
6023 else if (SCM_FRACTIONP (z
))
6024 return scm_from_double (scm_i_fraction2double (z
));
6025 else if (SCM_INEXACTP (z
))
6028 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
6032 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
6034 "Return an exact number that is numerically closest to @var{z}.")
6035 #define FUNC_NAME s_scm_inexact_to_exact
6037 if (SCM_I_INUMP (z
))
6039 else if (SCM_BIGP (z
))
6041 else if (SCM_REALP (z
))
6043 if (isinf (SCM_REAL_VALUE (z
)) || isnan (SCM_REAL_VALUE (z
)))
6044 SCM_OUT_OF_RANGE (1, z
);
6051 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
6052 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
6053 scm_i_mpz2num (mpq_denref (frac
)));
6055 /* When scm_i_make_ratio throws, we leak the memory allocated
6062 else if (SCM_FRACTIONP (z
))
6065 SCM_WRONG_TYPE_ARG (1, z
);
6069 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
6071 "Returns the @emph{simplest} rational number differing\n"
6072 "from @var{x} by no more than @var{eps}.\n"
6074 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
6075 "exact result when both its arguments are exact. Thus, you might need\n"
6076 "to use @code{inexact->exact} on the arguments.\n"
6079 "(rationalize (inexact->exact 1.2) 1/100)\n"
6082 #define FUNC_NAME s_scm_rationalize
6084 if (SCM_I_INUMP (x
))
6086 else if (SCM_BIGP (x
))
6088 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
6090 /* Use continued fractions to find closest ratio. All
6091 arithmetic is done with exact numbers.
6094 SCM ex
= scm_inexact_to_exact (x
);
6095 SCM int_part
= scm_floor (ex
);
6096 SCM tt
= SCM_I_MAKINUM (1);
6097 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
6098 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
6102 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
6105 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
6106 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
6108 /* We stop after a million iterations just to be absolutely sure
6109 that we don't go into an infinite loop. The process normally
6110 converges after less than a dozen iterations.
6113 eps
= scm_abs (eps
);
6114 while (++i
< 1000000)
6116 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
6117 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
6118 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
6120 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
6121 eps
))) /* abs(x-a/b) <= eps */
6123 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
6124 if (scm_is_false (scm_exact_p (x
))
6125 || scm_is_false (scm_exact_p (eps
)))
6126 return scm_exact_to_inexact (res
);
6130 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
6132 tt
= scm_floor (rx
); /* tt = floor (rx) */
6138 scm_num_overflow (s_scm_rationalize
);
6141 SCM_WRONG_TYPE_ARG (1, x
);
6145 /* conversion functions */
6148 scm_is_integer (SCM val
)
6150 return scm_is_true (scm_integer_p (val
));
6154 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
6156 if (SCM_I_INUMP (val
))
6158 scm_t_signed_bits n
= SCM_I_INUM (val
);
6159 return n
>= min
&& n
<= max
;
6161 else if (SCM_BIGP (val
))
6163 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
6165 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
6167 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
6169 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
6170 return n
>= min
&& n
<= max
;
6180 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6181 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6184 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6185 SCM_I_BIG_MPZ (val
));
6187 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
6199 return n
>= min
&& n
<= max
;
6207 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
6209 if (SCM_I_INUMP (val
))
6211 scm_t_signed_bits n
= SCM_I_INUM (val
);
6212 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
6214 else if (SCM_BIGP (val
))
6216 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
6218 else if (max
<= ULONG_MAX
)
6220 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
6222 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
6223 return n
>= min
&& n
<= max
;
6233 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
6236 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6237 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6240 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6241 SCM_I_BIG_MPZ (val
));
6243 return n
>= min
&& n
<= max
;
6251 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
6253 scm_error (scm_out_of_range_key
,
6255 "Value out of range ~S to ~S: ~S",
6256 scm_list_3 (min
, max
, bad_val
),
6257 scm_list_1 (bad_val
));
6260 #define TYPE scm_t_intmax
6261 #define TYPE_MIN min
6262 #define TYPE_MAX max
6263 #define SIZEOF_TYPE 0
6264 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
6265 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
6266 #include "libguile/conv-integer.i.c"
6268 #define TYPE scm_t_uintmax
6269 #define TYPE_MIN min
6270 #define TYPE_MAX max
6271 #define SIZEOF_TYPE 0
6272 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6273 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6274 #include "libguile/conv-uinteger.i.c"
6276 #define TYPE scm_t_int8
6277 #define TYPE_MIN SCM_T_INT8_MIN
6278 #define TYPE_MAX SCM_T_INT8_MAX
6279 #define SIZEOF_TYPE 1
6280 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6281 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6282 #include "libguile/conv-integer.i.c"
6284 #define TYPE scm_t_uint8
6286 #define TYPE_MAX SCM_T_UINT8_MAX
6287 #define SIZEOF_TYPE 1
6288 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6289 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6290 #include "libguile/conv-uinteger.i.c"
6292 #define TYPE scm_t_int16
6293 #define TYPE_MIN SCM_T_INT16_MIN
6294 #define TYPE_MAX SCM_T_INT16_MAX
6295 #define SIZEOF_TYPE 2
6296 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6297 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6298 #include "libguile/conv-integer.i.c"
6300 #define TYPE scm_t_uint16
6302 #define TYPE_MAX SCM_T_UINT16_MAX
6303 #define SIZEOF_TYPE 2
6304 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6305 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6306 #include "libguile/conv-uinteger.i.c"
6308 #define TYPE scm_t_int32
6309 #define TYPE_MIN SCM_T_INT32_MIN
6310 #define TYPE_MAX SCM_T_INT32_MAX
6311 #define SIZEOF_TYPE 4
6312 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6313 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6314 #include "libguile/conv-integer.i.c"
6316 #define TYPE scm_t_uint32
6318 #define TYPE_MAX SCM_T_UINT32_MAX
6319 #define SIZEOF_TYPE 4
6320 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6321 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6322 #include "libguile/conv-uinteger.i.c"
6324 #define TYPE scm_t_wchar
6325 #define TYPE_MIN (scm_t_int32)-1
6326 #define TYPE_MAX (scm_t_int32)0x10ffff
6327 #define SIZEOF_TYPE 4
6328 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6329 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6330 #include "libguile/conv-integer.i.c"
6332 #define TYPE scm_t_int64
6333 #define TYPE_MIN SCM_T_INT64_MIN
6334 #define TYPE_MAX SCM_T_INT64_MAX
6335 #define SIZEOF_TYPE 8
6336 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6337 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6338 #include "libguile/conv-integer.i.c"
6340 #define TYPE scm_t_uint64
6342 #define TYPE_MAX SCM_T_UINT64_MAX
6343 #define SIZEOF_TYPE 8
6344 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6345 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6346 #include "libguile/conv-uinteger.i.c"
6349 scm_to_mpz (SCM val
, mpz_t rop
)
6351 if (SCM_I_INUMP (val
))
6352 mpz_set_si (rop
, SCM_I_INUM (val
));
6353 else if (SCM_BIGP (val
))
6354 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6356 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6360 scm_from_mpz (mpz_t val
)
6362 return scm_i_mpz2num (val
);
6366 scm_is_real (SCM val
)
6368 return scm_is_true (scm_real_p (val
));
6372 scm_is_rational (SCM val
)
6374 return scm_is_true (scm_rational_p (val
));
6378 scm_to_double (SCM val
)
6380 if (SCM_I_INUMP (val
))
6381 return SCM_I_INUM (val
);
6382 else if (SCM_BIGP (val
))
6383 return scm_i_big2dbl (val
);
6384 else if (SCM_FRACTIONP (val
))
6385 return scm_i_fraction2double (val
);
6386 else if (SCM_REALP (val
))
6387 return SCM_REAL_VALUE (val
);
6389 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6393 scm_from_double (double val
)
6397 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
6399 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
6400 SCM_REAL_VALUE (z
) = val
;
6405 #if SCM_ENABLE_DEPRECATED == 1
6408 scm_num2float (SCM num
, unsigned long pos
, const char *s_caller
)
6410 scm_c_issue_deprecation_warning
6411 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
6415 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6419 scm_out_of_range (NULL
, num
);
6422 return scm_to_double (num
);
6426 scm_num2double (SCM num
, unsigned long pos
, const char *s_caller
)
6428 scm_c_issue_deprecation_warning
6429 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
6433 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6437 scm_out_of_range (NULL
, num
);
6440 return scm_to_double (num
);
6446 scm_is_complex (SCM val
)
6448 return scm_is_true (scm_complex_p (val
));
6452 scm_c_real_part (SCM z
)
6454 if (SCM_COMPLEXP (z
))
6455 return SCM_COMPLEX_REAL (z
);
6458 /* Use the scm_real_part to get proper error checking and
6461 return scm_to_double (scm_real_part (z
));
6466 scm_c_imag_part (SCM z
)
6468 if (SCM_COMPLEXP (z
))
6469 return SCM_COMPLEX_IMAG (z
);
6472 /* Use the scm_imag_part to get proper error checking and
6473 dispatching. The result will almost always be 0.0, but not
6476 return scm_to_double (scm_imag_part (z
));
6481 scm_c_magnitude (SCM z
)
6483 return scm_to_double (scm_magnitude (z
));
6489 return scm_to_double (scm_angle (z
));
6493 scm_is_number (SCM z
)
6495 return scm_is_true (scm_number_p (z
));
6499 /* In the following functions we dispatch to the real-arg funcs like log()
6500 when we know the arg is real, instead of just handing everything to
6501 clog() for instance. This is in case clog() doesn't optimize for a
6502 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6503 well use it to go straight to the applicable C func. */
6505 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6507 "Return the natural logarithm of @var{z}.")
6508 #define FUNC_NAME s_scm_log
6510 if (SCM_COMPLEXP (z
))
6512 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6513 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6515 double re
= SCM_COMPLEX_REAL (z
);
6516 double im
= SCM_COMPLEX_IMAG (z
);
6517 return scm_c_make_rectangular (log (hypot (re
, im
)),
6523 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6524 although the value itself overflows. */
6525 double re
= scm_to_double (z
);
6526 double l
= log (fabs (re
));
6528 return scm_from_double (l
);
6530 return scm_c_make_rectangular (l
, M_PI
);
6536 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6538 "Return the base 10 logarithm of @var{z}.")
6539 #define FUNC_NAME s_scm_log10
6541 if (SCM_COMPLEXP (z
))
6543 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6544 clog() and a multiply by M_LOG10E, rather than the fallback
6545 log10+hypot+atan2.) */
6546 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
6547 && defined SCM_COMPLEX_VALUE
6548 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6550 double re
= SCM_COMPLEX_REAL (z
);
6551 double im
= SCM_COMPLEX_IMAG (z
);
6552 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6553 M_LOG10E
* atan2 (im
, re
));
6558 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6559 although the value itself overflows. */
6560 double re
= scm_to_double (z
);
6561 double l
= log10 (fabs (re
));
6563 return scm_from_double (l
);
6565 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6571 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6573 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6574 "base of natural logarithms (2.71828@dots{}).")
6575 #define FUNC_NAME s_scm_exp
6577 if (SCM_COMPLEXP (z
))
6579 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6580 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6582 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6583 SCM_COMPLEX_IMAG (z
));
6588 /* When z is a negative bignum the conversion to double overflows,
6589 giving -infinity, but that's ok, the exp is still 0.0. */
6590 return scm_from_double (exp (scm_to_double (z
)));
6596 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6598 "Return the square root of @var{z}. Of the two possible roots\n"
6599 "(positive and negative), the one with the a positive real part\n"
6600 "is returned, or if that's zero then a positive imaginary part.\n"
6604 "(sqrt 9.0) @result{} 3.0\n"
6605 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6606 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6607 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6609 #define FUNC_NAME s_scm_sqrt
6611 if (SCM_COMPLEXP (x
))
6613 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
6614 && defined SCM_COMPLEX_VALUE
6615 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6617 double re
= SCM_COMPLEX_REAL (x
);
6618 double im
= SCM_COMPLEX_IMAG (x
);
6619 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6620 0.5 * atan2 (im
, re
));
6625 double xx
= scm_to_double (x
);
6627 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6629 return scm_from_double (sqrt (xx
));
6641 mpz_init_set_si (z_negative_one
, -1);
6643 /* It may be possible to tune the performance of some algorithms by using
6644 * the following constants to avoid the creation of bignums. Please, before
6645 * using these values, remember the two rules of program optimization:
6646 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6647 scm_c_define ("most-positive-fixnum",
6648 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6649 scm_c_define ("most-negative-fixnum",
6650 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6652 scm_add_feature ("complex");
6653 scm_add_feature ("inexact");
6654 flo0
= scm_from_double (0.0);
6656 /* determine floating point precision */
6657 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6659 init_dblprec(&scm_dblprec
[i
-2],i
);
6660 init_fx_radix(fx_per_radix
[i
-2],i
);
6663 /* hard code precision for base 10 if the preprocessor tells us to... */
6664 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6667 exactly_one_half
= scm_divide (SCM_I_MAKINUM (1), SCM_I_MAKINUM (2));
6668 #include "libguile/numbers.x"