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))
82 /* Tests to see if a C double is neither infinite nor a NaN.
83 TODO: if it's available, use C99's isfinite(x) instead */
84 #define DOUBLE_IS_FINITE(x) (!isinf(x) && !isnan(x))
89 Wonder if this might be faster for some of our code? A switch on
90 the numtag would jump directly to the right case, and the
91 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
93 #define SCM_I_NUMTAG_NOTNUM 0
94 #define SCM_I_NUMTAG_INUM 1
95 #define SCM_I_NUMTAG_BIG scm_tc16_big
96 #define SCM_I_NUMTAG_REAL scm_tc16_real
97 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
98 #define SCM_I_NUMTAG(x) \
99 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
100 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
101 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
102 : SCM_I_NUMTAG_NOTNUM)))
104 /* the macro above will not work as is with fractions */
109 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
111 /* FLOBUFLEN is the maximum number of characters neccessary for the
112 * printed or scm_string representation of an inexact number.
114 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
117 #if !defined (HAVE_ASINH)
118 static double asinh (double x
) { return log (x
+ sqrt (x
* x
+ 1)); }
120 #if !defined (HAVE_ACOSH)
121 static double acosh (double x
) { return log (x
+ sqrt (x
* x
- 1)); }
123 #if !defined (HAVE_ATANH)
124 static double atanh (double x
) { return 0.5 * log ((1 + x
) / (1 - x
)); }
127 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
128 an explicit check. In some future gmp (don't know what version number),
129 mpz_cmp_d is supposed to do this itself. */
131 #define xmpz_cmp_d(z, d) \
132 (isinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
134 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
138 #if defined (GUILE_I)
139 #if HAVE_COMPLEX_DOUBLE
141 /* For an SCM object Z which is a complex number (ie. satisfies
142 SCM_COMPLEXP), return its value as a C level "complex double". */
143 #define SCM_COMPLEX_VALUE(z) \
144 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
146 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
148 /* Convert a C "complex double" to an SCM value. */
150 scm_from_complex_double (complex double z
)
152 return scm_c_make_rectangular (creal (z
), cimag (z
));
155 #endif /* HAVE_COMPLEX_DOUBLE */
160 static mpz_t z_negative_one
;
163 /* Clear the `mpz_t' embedded in bignum PTR. */
165 finalize_bignum (GC_PTR ptr
, GC_PTR data
)
169 bignum
= PTR2SCM (ptr
);
170 mpz_clear (SCM_I_BIG_MPZ (bignum
));
173 /* Return a new uninitialized bignum. */
178 GC_finalization_proc prev_finalizer
;
179 GC_PTR prev_finalizer_data
;
181 /* Allocate one word for the type tag and enough room for an `mpz_t'. */
182 p
= scm_gc_malloc_pointerless (sizeof (scm_t_bits
) + sizeof (mpz_t
),
186 GC_REGISTER_FINALIZER_NO_ORDER (p
, finalize_bignum
, NULL
,
188 &prev_finalizer_data
);
197 /* Return a newly created bignum. */
198 SCM z
= make_bignum ();
199 mpz_init (SCM_I_BIG_MPZ (z
));
204 scm_i_inum2big (scm_t_inum x
)
206 /* Return a newly created bignum initialized to X. */
207 SCM z
= make_bignum ();
208 #if SIZEOF_VOID_P == SIZEOF_LONG
209 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
211 /* Note that in this case, you'll also have to check all mpz_*_ui and
212 mpz_*_si invocations in Guile. */
213 #error creation of mpz not implemented for this inum size
219 scm_i_long2big (long x
)
221 /* Return a newly created bignum initialized to X. */
222 SCM z
= make_bignum ();
223 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
228 scm_i_ulong2big (unsigned long x
)
230 /* Return a newly created bignum initialized to X. */
231 SCM z
= make_bignum ();
232 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
237 scm_i_clonebig (SCM src_big
, int same_sign_p
)
239 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
240 SCM z
= make_bignum ();
241 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
243 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
248 scm_i_bigcmp (SCM x
, SCM y
)
250 /* Return neg if x < y, pos if x > y, and 0 if x == y */
251 /* presume we already know x and y are bignums */
252 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
253 scm_remember_upto_here_2 (x
, y
);
258 scm_i_dbl2big (double d
)
260 /* results are only defined if d is an integer */
261 SCM z
= make_bignum ();
262 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
266 /* Convert a integer in double representation to a SCM number. */
269 scm_i_dbl2num (double u
)
271 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
272 powers of 2, so there's no rounding when making "double" values
273 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
274 get rounded on a 64-bit machine, hence the "+1".
276 The use of floor() to force to an integer value ensures we get a
277 "numerically closest" value without depending on how a
278 double->long cast or how mpz_set_d will round. For reference,
279 double->long probably follows the hardware rounding mode,
280 mpz_set_d truncates towards zero. */
282 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
283 representable as a double? */
285 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
286 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
287 return SCM_I_MAKINUM ((scm_t_inum
) u
);
289 return scm_i_dbl2big (u
);
292 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
293 with R5RS exact->inexact.
295 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
296 (ie. truncate towards zero), then adjust to get the closest double by
297 examining the next lower bit and adding 1 (to the absolute value) if
300 Bignums exactly half way between representable doubles are rounded to the
301 next higher absolute value (ie. away from zero). This seems like an
302 adequate interpretation of R5RS "numerically closest", and it's easier
303 and faster than a full "nearest-even" style.
305 The bit test must be done on the absolute value of the mpz_t, which means
306 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
307 negatives as twos complement.
309 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
310 following the hardware rounding mode, but applied to the absolute value
311 of the mpz_t operand. This is not what we want so we put the high
312 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
313 mpz_get_d is supposed to always truncate towards zero.
315 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
316 is a slowdown. It'd be faster to pick out the relevant high bits with
317 mpz_getlimbn if we could be bothered coding that, and if the new
318 truncating gmp doesn't come out. */
321 scm_i_big2dbl (SCM b
)
326 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
330 /* Current GMP, eg. 4.1.3, force truncation towards zero */
332 if (bits
> DBL_MANT_DIG
)
334 size_t shift
= bits
- DBL_MANT_DIG
;
335 mpz_init2 (tmp
, DBL_MANT_DIG
);
336 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
337 result
= ldexp (mpz_get_d (tmp
), shift
);
342 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
347 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
350 if (bits
> DBL_MANT_DIG
)
352 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
353 /* test bit number "pos" in absolute value */
354 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
355 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
357 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
361 scm_remember_upto_here_1 (b
);
366 scm_i_normbig (SCM b
)
368 /* convert a big back to a fixnum if it'll fit */
369 /* presume b is a bignum */
370 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
372 scm_t_inum val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
373 if (SCM_FIXABLE (val
))
374 b
= SCM_I_MAKINUM (val
);
379 static SCM_C_INLINE_KEYWORD SCM
380 scm_i_mpz2num (mpz_t b
)
382 /* convert a mpz number to a SCM number. */
383 if (mpz_fits_slong_p (b
))
385 scm_t_inum val
= mpz_get_si (b
);
386 if (SCM_FIXABLE (val
))
387 return SCM_I_MAKINUM (val
);
391 SCM z
= make_bignum ();
392 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
397 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
398 static SCM
scm_divide2real (SCM x
, SCM y
);
401 scm_i_make_ratio (SCM numerator
, SCM denominator
)
402 #define FUNC_NAME "make-ratio"
404 /* First make sure the arguments are proper.
406 if (SCM_I_INUMP (denominator
))
408 if (scm_is_eq (denominator
, SCM_INUM0
))
409 scm_num_overflow ("make-ratio");
410 if (scm_is_eq (denominator
, SCM_INUM1
))
415 if (!(SCM_BIGP(denominator
)))
416 SCM_WRONG_TYPE_ARG (2, denominator
);
418 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
419 SCM_WRONG_TYPE_ARG (1, numerator
);
421 /* Then flip signs so that the denominator is positive.
423 if (scm_is_true (scm_negative_p (denominator
)))
425 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
426 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
429 /* Now consider for each of the four fixnum/bignum combinations
430 whether the rational number is really an integer.
432 if (SCM_I_INUMP (numerator
))
434 scm_t_inum x
= SCM_I_INUM (numerator
);
435 if (scm_is_eq (numerator
, SCM_INUM0
))
437 if (SCM_I_INUMP (denominator
))
440 y
= SCM_I_INUM (denominator
);
444 return SCM_I_MAKINUM (x
/ y
);
448 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
449 of that value for the denominator, as a bignum. Apart from
450 that case, abs(bignum) > abs(inum) so inum/bignum is not an
452 if (x
== SCM_MOST_NEGATIVE_FIXNUM
453 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
454 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
455 return SCM_I_MAKINUM(-1);
458 else if (SCM_BIGP (numerator
))
460 if (SCM_I_INUMP (denominator
))
462 scm_t_inum yy
= SCM_I_INUM (denominator
);
463 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
464 return scm_divide (numerator
, denominator
);
468 if (scm_is_eq (numerator
, denominator
))
470 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
471 SCM_I_BIG_MPZ (denominator
)))
472 return scm_divide(numerator
, denominator
);
476 /* No, it's a proper fraction.
479 SCM divisor
= scm_gcd (numerator
, denominator
);
480 if (!(scm_is_eq (divisor
, SCM_INUM1
)))
482 numerator
= scm_divide (numerator
, divisor
);
483 denominator
= scm_divide (denominator
, divisor
);
486 return scm_double_cell (scm_tc16_fraction
,
487 SCM_UNPACK (numerator
),
488 SCM_UNPACK (denominator
), 0);
494 scm_i_fraction2double (SCM z
)
496 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
497 SCM_FRACTION_DENOMINATOR (z
)));
500 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
502 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
504 #define FUNC_NAME s_scm_exact_p
506 if (SCM_INEXACTP (x
))
508 else if (SCM_NUMBERP (x
))
511 SCM_WRONG_TYPE_ARG (1, x
);
516 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
518 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
520 #define FUNC_NAME s_scm_inexact_p
522 if (SCM_INEXACTP (x
))
524 else if (SCM_NUMBERP (x
))
527 SCM_WRONG_TYPE_ARG (1, x
);
532 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
534 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
536 #define FUNC_NAME s_scm_odd_p
540 scm_t_inum val
= SCM_I_INUM (n
);
541 return scm_from_bool ((val
& 1L) != 0);
543 else if (SCM_BIGP (n
))
545 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
546 scm_remember_upto_here_1 (n
);
547 return scm_from_bool (odd_p
);
549 else if (scm_is_true (scm_inf_p (n
)))
550 SCM_WRONG_TYPE_ARG (1, n
);
551 else if (SCM_REALP (n
))
553 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
559 SCM_WRONG_TYPE_ARG (1, n
);
562 SCM_WRONG_TYPE_ARG (1, n
);
567 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
569 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
571 #define FUNC_NAME s_scm_even_p
575 scm_t_inum val
= SCM_I_INUM (n
);
576 return scm_from_bool ((val
& 1L) == 0);
578 else if (SCM_BIGP (n
))
580 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
581 scm_remember_upto_here_1 (n
);
582 return scm_from_bool (even_p
);
584 else if (scm_is_true (scm_inf_p (n
)))
585 SCM_WRONG_TYPE_ARG (1, n
);
586 else if (SCM_REALP (n
))
588 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
594 SCM_WRONG_TYPE_ARG (1, n
);
597 SCM_WRONG_TYPE_ARG (1, n
);
601 SCM_DEFINE (scm_finite_p
, "finite?", 1, 0, 0,
603 "Return @code{#t} if the real number @var{x} is neither\n"
604 "infinite nor a NaN, @code{#f} otherwise.")
605 #define FUNC_NAME s_scm_finite_p
608 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
609 else if (scm_is_real (x
))
612 SCM_WRONG_TYPE_ARG (1, x
);
616 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
618 "Return @code{#t} if the real number @var{x} is @samp{+inf.0} or\n"
619 "@samp{-inf.0}. Otherwise return @code{#f}.")
620 #define FUNC_NAME s_scm_inf_p
623 return scm_from_bool (isinf (SCM_REAL_VALUE (x
)));
624 else if (scm_is_real (x
))
627 SCM_WRONG_TYPE_ARG (1, x
);
631 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
633 "Return @code{#t} if the real number @var{x} is a NaN,\n"
634 "or @code{#f} otherwise.")
635 #define FUNC_NAME s_scm_nan_p
638 return scm_from_bool (isnan (SCM_REAL_VALUE (x
)));
639 else if (scm_is_real (x
))
642 SCM_WRONG_TYPE_ARG (1, x
);
646 /* Guile's idea of infinity. */
647 static double guile_Inf
;
649 /* Guile's idea of not a number. */
650 static double guile_NaN
;
653 guile_ieee_init (void)
655 /* Some version of gcc on some old version of Linux used to crash when
656 trying to make Inf and NaN. */
659 /* C99 INFINITY, when available.
660 FIXME: The standard allows for INFINITY to be something that overflows
661 at compile time. We ought to have a configure test to check for that
662 before trying to use it. (But in practice we believe this is not a
663 problem on any system guile is likely to target.) */
664 guile_Inf
= INFINITY
;
665 #elif defined HAVE_DINFINITY
667 extern unsigned int DINFINITY
[2];
668 guile_Inf
= (*((double *) (DINFINITY
)));
675 if (guile_Inf
== tmp
)
682 /* C99 NAN, when available */
684 #elif defined HAVE_DQNAN
687 extern unsigned int DQNAN
[2];
688 guile_NaN
= (*((double *)(DQNAN
)));
691 guile_NaN
= guile_Inf
/ guile_Inf
;
695 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
698 #define FUNC_NAME s_scm_inf
700 static int initialized
= 0;
706 return scm_from_double (guile_Inf
);
710 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
713 #define FUNC_NAME s_scm_nan
715 static int initialized
= 0;
721 return scm_from_double (guile_NaN
);
726 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
728 "Return the absolute value of @var{x}.")
733 scm_t_inum xx
= SCM_I_INUM (x
);
736 else if (SCM_POSFIXABLE (-xx
))
737 return SCM_I_MAKINUM (-xx
);
739 return scm_i_inum2big (-xx
);
741 else if (SCM_BIGP (x
))
743 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
745 return scm_i_clonebig (x
, 0);
749 else if (SCM_REALP (x
))
751 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
752 double xx
= SCM_REAL_VALUE (x
);
754 return scm_from_double (-xx
);
758 else if (SCM_FRACTIONP (x
))
760 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
762 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
763 SCM_FRACTION_DENOMINATOR (x
));
766 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
771 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
772 /* "Return the quotient of the numbers @var{x} and @var{y}."
775 scm_quotient (SCM x
, SCM y
)
777 if (SCM_LIKELY (SCM_I_INUMP (x
)))
779 scm_t_inum xx
= SCM_I_INUM (x
);
780 if (SCM_LIKELY (SCM_I_INUMP (y
)))
782 scm_t_inum yy
= SCM_I_INUM (y
);
783 if (SCM_UNLIKELY (yy
== 0))
784 scm_num_overflow (s_quotient
);
787 scm_t_inum z
= xx
/ yy
;
788 if (SCM_LIKELY (SCM_FIXABLE (z
)))
789 return SCM_I_MAKINUM (z
);
791 return scm_i_inum2big (z
);
794 else if (SCM_BIGP (y
))
796 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
797 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
798 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
800 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
801 scm_remember_upto_here_1 (y
);
802 return SCM_I_MAKINUM (-1);
808 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
810 else if (SCM_BIGP (x
))
812 if (SCM_LIKELY (SCM_I_INUMP (y
)))
814 scm_t_inum yy
= SCM_I_INUM (y
);
815 if (SCM_UNLIKELY (yy
== 0))
816 scm_num_overflow (s_quotient
);
817 else if (SCM_UNLIKELY (yy
== 1))
821 SCM result
= scm_i_mkbig ();
824 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
827 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
830 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
831 scm_remember_upto_here_1 (x
);
832 return scm_i_normbig (result
);
835 else if (SCM_BIGP (y
))
837 SCM result
= scm_i_mkbig ();
838 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
841 scm_remember_upto_here_2 (x
, y
);
842 return scm_i_normbig (result
);
845 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
848 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
851 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
852 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
854 * "(remainder 13 4) @result{} 1\n"
855 * "(remainder -13 4) @result{} -1\n"
859 scm_remainder (SCM x
, SCM y
)
861 if (SCM_LIKELY (SCM_I_INUMP (x
)))
863 if (SCM_LIKELY (SCM_I_INUMP (y
)))
865 scm_t_inum yy
= SCM_I_INUM (y
);
866 if (SCM_UNLIKELY (yy
== 0))
867 scm_num_overflow (s_remainder
);
870 /* C99 specifies that "%" is the remainder corresponding to a
871 quotient rounded towards zero, and that's also traditional
872 for machine division, so z here should be well defined. */
873 scm_t_inum z
= SCM_I_INUM (x
) % yy
;
874 return SCM_I_MAKINUM (z
);
877 else if (SCM_BIGP (y
))
879 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
880 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
881 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
883 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
884 scm_remember_upto_here_1 (y
);
891 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
893 else if (SCM_BIGP (x
))
895 if (SCM_LIKELY (SCM_I_INUMP (y
)))
897 scm_t_inum yy
= SCM_I_INUM (y
);
898 if (SCM_UNLIKELY (yy
== 0))
899 scm_num_overflow (s_remainder
);
902 SCM result
= scm_i_mkbig ();
905 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
906 scm_remember_upto_here_1 (x
);
907 return scm_i_normbig (result
);
910 else if (SCM_BIGP (y
))
912 SCM result
= scm_i_mkbig ();
913 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
916 scm_remember_upto_here_2 (x
, y
);
917 return scm_i_normbig (result
);
920 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
923 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
927 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
928 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
930 * "(modulo 13 4) @result{} 1\n"
931 * "(modulo -13 4) @result{} 3\n"
935 scm_modulo (SCM x
, SCM y
)
937 if (SCM_LIKELY (SCM_I_INUMP (x
)))
939 scm_t_inum xx
= SCM_I_INUM (x
);
940 if (SCM_LIKELY (SCM_I_INUMP (y
)))
942 scm_t_inum yy
= SCM_I_INUM (y
);
943 if (SCM_UNLIKELY (yy
== 0))
944 scm_num_overflow (s_modulo
);
947 /* C99 specifies that "%" is the remainder corresponding to a
948 quotient rounded towards zero, and that's also traditional
949 for machine division, so z here should be well defined. */
950 scm_t_inum z
= xx
% yy
;
967 return SCM_I_MAKINUM (result
);
970 else if (SCM_BIGP (y
))
972 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
979 SCM pos_y
= scm_i_clonebig (y
, 0);
980 /* do this after the last scm_op */
981 mpz_init_set_si (z_x
, xx
);
982 result
= pos_y
; /* re-use this bignum */
983 mpz_mod (SCM_I_BIG_MPZ (result
),
985 SCM_I_BIG_MPZ (pos_y
));
986 scm_remember_upto_here_1 (pos_y
);
990 result
= scm_i_mkbig ();
991 /* do this after the last scm_op */
992 mpz_init_set_si (z_x
, xx
);
993 mpz_mod (SCM_I_BIG_MPZ (result
),
996 scm_remember_upto_here_1 (y
);
999 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1000 mpz_add (SCM_I_BIG_MPZ (result
),
1002 SCM_I_BIG_MPZ (result
));
1003 scm_remember_upto_here_1 (y
);
1004 /* and do this before the next one */
1006 return scm_i_normbig (result
);
1010 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1012 else if (SCM_BIGP (x
))
1014 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1016 scm_t_inum yy
= SCM_I_INUM (y
);
1017 if (SCM_UNLIKELY (yy
== 0))
1018 scm_num_overflow (s_modulo
);
1021 SCM result
= scm_i_mkbig ();
1022 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
1024 (yy
< 0) ? - yy
: yy
);
1025 scm_remember_upto_here_1 (x
);
1026 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1027 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
1028 SCM_I_BIG_MPZ (result
),
1030 return scm_i_normbig (result
);
1033 else if (SCM_BIGP (y
))
1035 SCM result
= scm_i_mkbig ();
1036 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
1037 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1038 mpz_mod (SCM_I_BIG_MPZ (result
),
1040 SCM_I_BIG_MPZ (pos_y
));
1042 scm_remember_upto_here_1 (x
);
1043 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1044 mpz_add (SCM_I_BIG_MPZ (result
),
1046 SCM_I_BIG_MPZ (result
));
1047 scm_remember_upto_here_2 (y
, pos_y
);
1048 return scm_i_normbig (result
);
1051 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1054 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1057 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
1058 (SCM x
, SCM y
, SCM rest
),
1059 "Return the greatest common divisor of all parameter values.\n"
1060 "If called without arguments, 0 is returned.")
1061 #define FUNC_NAME s_scm_i_gcd
1063 while (!scm_is_null (rest
))
1064 { x
= scm_gcd (x
, y
);
1066 rest
= scm_cdr (rest
);
1068 return scm_gcd (x
, y
);
1072 #define s_gcd s_scm_i_gcd
1073 #define g_gcd g_scm_i_gcd
1076 scm_gcd (SCM x
, SCM y
)
1079 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
1081 if (SCM_I_INUMP (x
))
1083 if (SCM_I_INUMP (y
))
1085 scm_t_inum xx
= SCM_I_INUM (x
);
1086 scm_t_inum yy
= SCM_I_INUM (y
);
1087 scm_t_inum u
= xx
< 0 ? -xx
: xx
;
1088 scm_t_inum v
= yy
< 0 ? -yy
: yy
;
1098 /* Determine a common factor 2^k */
1099 while (!(1 & (u
| v
)))
1105 /* Now, any factor 2^n can be eliminated */
1125 return (SCM_POSFIXABLE (result
)
1126 ? SCM_I_MAKINUM (result
)
1127 : scm_i_inum2big (result
));
1129 else if (SCM_BIGP (y
))
1135 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1137 else if (SCM_BIGP (x
))
1139 if (SCM_I_INUMP (y
))
1144 yy
= SCM_I_INUM (y
);
1149 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1150 scm_remember_upto_here_1 (x
);
1151 return (SCM_POSFIXABLE (result
)
1152 ? SCM_I_MAKINUM (result
)
1153 : scm_from_unsigned_integer (result
));
1155 else if (SCM_BIGP (y
))
1157 SCM result
= scm_i_mkbig ();
1158 mpz_gcd (SCM_I_BIG_MPZ (result
),
1161 scm_remember_upto_here_2 (x
, y
);
1162 return scm_i_normbig (result
);
1165 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1168 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1171 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
1172 (SCM x
, SCM y
, SCM rest
),
1173 "Return the least common multiple of the arguments.\n"
1174 "If called without arguments, 1 is returned.")
1175 #define FUNC_NAME s_scm_i_lcm
1177 while (!scm_is_null (rest
))
1178 { x
= scm_lcm (x
, y
);
1180 rest
= scm_cdr (rest
);
1182 return scm_lcm (x
, y
);
1186 #define s_lcm s_scm_i_lcm
1187 #define g_lcm g_scm_i_lcm
1190 scm_lcm (SCM n1
, SCM n2
)
1192 if (SCM_UNBNDP (n2
))
1194 if (SCM_UNBNDP (n1
))
1195 return SCM_I_MAKINUM (1L);
1196 n2
= SCM_I_MAKINUM (1L);
1199 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1200 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1201 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1202 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1204 if (SCM_I_INUMP (n1
))
1206 if (SCM_I_INUMP (n2
))
1208 SCM d
= scm_gcd (n1
, n2
);
1209 if (scm_is_eq (d
, SCM_INUM0
))
1212 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1216 /* inum n1, big n2 */
1219 SCM result
= scm_i_mkbig ();
1220 scm_t_inum nn1
= SCM_I_INUM (n1
);
1221 if (nn1
== 0) return SCM_INUM0
;
1222 if (nn1
< 0) nn1
= - nn1
;
1223 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1224 scm_remember_upto_here_1 (n2
);
1232 if (SCM_I_INUMP (n2
))
1239 SCM result
= scm_i_mkbig ();
1240 mpz_lcm(SCM_I_BIG_MPZ (result
),
1242 SCM_I_BIG_MPZ (n2
));
1243 scm_remember_upto_here_2(n1
, n2
);
1244 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1250 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1255 + + + x (map digit:logand X Y)
1256 + - + x (map digit:logand X (lognot (+ -1 Y)))
1257 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1258 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1263 + + + (map digit:logior X Y)
1264 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1265 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1266 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1271 + + + (map digit:logxor X Y)
1272 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1273 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1274 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1279 + + (any digit:logand X Y)
1280 + - (any digit:logand X (lognot (+ -1 Y)))
1281 - + (any digit:logand (lognot (+ -1 X)) Y)
1286 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
1287 (SCM x
, SCM y
, SCM rest
),
1288 "Return the bitwise AND of the integer arguments.\n\n"
1290 "(logand) @result{} -1\n"
1291 "(logand 7) @result{} 7\n"
1292 "(logand #b111 #b011 #b001) @result{} 1\n"
1294 #define FUNC_NAME s_scm_i_logand
1296 while (!scm_is_null (rest
))
1297 { x
= scm_logand (x
, y
);
1299 rest
= scm_cdr (rest
);
1301 return scm_logand (x
, y
);
1305 #define s_scm_logand s_scm_i_logand
1307 SCM
scm_logand (SCM n1
, SCM n2
)
1308 #define FUNC_NAME s_scm_logand
1312 if (SCM_UNBNDP (n2
))
1314 if (SCM_UNBNDP (n1
))
1315 return SCM_I_MAKINUM (-1);
1316 else if (!SCM_NUMBERP (n1
))
1317 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1318 else if (SCM_NUMBERP (n1
))
1321 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1324 if (SCM_I_INUMP (n1
))
1326 nn1
= SCM_I_INUM (n1
);
1327 if (SCM_I_INUMP (n2
))
1329 scm_t_inum nn2
= SCM_I_INUM (n2
);
1330 return SCM_I_MAKINUM (nn1
& nn2
);
1332 else if SCM_BIGP (n2
)
1338 SCM result_z
= scm_i_mkbig ();
1340 mpz_init_set_si (nn1_z
, nn1
);
1341 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1342 scm_remember_upto_here_1 (n2
);
1344 return scm_i_normbig (result_z
);
1348 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1350 else if (SCM_BIGP (n1
))
1352 if (SCM_I_INUMP (n2
))
1355 nn1
= SCM_I_INUM (n1
);
1358 else if (SCM_BIGP (n2
))
1360 SCM result_z
= scm_i_mkbig ();
1361 mpz_and (SCM_I_BIG_MPZ (result_z
),
1363 SCM_I_BIG_MPZ (n2
));
1364 scm_remember_upto_here_2 (n1
, n2
);
1365 return scm_i_normbig (result_z
);
1368 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1371 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1376 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
1377 (SCM x
, SCM y
, SCM rest
),
1378 "Return the bitwise OR of the integer arguments.\n\n"
1380 "(logior) @result{} 0\n"
1381 "(logior 7) @result{} 7\n"
1382 "(logior #b000 #b001 #b011) @result{} 3\n"
1384 #define FUNC_NAME s_scm_i_logior
1386 while (!scm_is_null (rest
))
1387 { x
= scm_logior (x
, y
);
1389 rest
= scm_cdr (rest
);
1391 return scm_logior (x
, y
);
1395 #define s_scm_logior s_scm_i_logior
1397 SCM
scm_logior (SCM n1
, SCM n2
)
1398 #define FUNC_NAME s_scm_logior
1402 if (SCM_UNBNDP (n2
))
1404 if (SCM_UNBNDP (n1
))
1406 else if (SCM_NUMBERP (n1
))
1409 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1412 if (SCM_I_INUMP (n1
))
1414 nn1
= SCM_I_INUM (n1
);
1415 if (SCM_I_INUMP (n2
))
1417 long nn2
= SCM_I_INUM (n2
);
1418 return SCM_I_MAKINUM (nn1
| nn2
);
1420 else if (SCM_BIGP (n2
))
1426 SCM result_z
= scm_i_mkbig ();
1428 mpz_init_set_si (nn1_z
, nn1
);
1429 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1430 scm_remember_upto_here_1 (n2
);
1432 return scm_i_normbig (result_z
);
1436 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1438 else if (SCM_BIGP (n1
))
1440 if (SCM_I_INUMP (n2
))
1443 nn1
= SCM_I_INUM (n1
);
1446 else if (SCM_BIGP (n2
))
1448 SCM result_z
= scm_i_mkbig ();
1449 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1451 SCM_I_BIG_MPZ (n2
));
1452 scm_remember_upto_here_2 (n1
, n2
);
1453 return scm_i_normbig (result_z
);
1456 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1459 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1464 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
1465 (SCM x
, SCM y
, SCM rest
),
1466 "Return the bitwise XOR of the integer arguments. A bit is\n"
1467 "set in the result if it is set in an odd number of arguments.\n"
1469 "(logxor) @result{} 0\n"
1470 "(logxor 7) @result{} 7\n"
1471 "(logxor #b000 #b001 #b011) @result{} 2\n"
1472 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1474 #define FUNC_NAME s_scm_i_logxor
1476 while (!scm_is_null (rest
))
1477 { x
= scm_logxor (x
, y
);
1479 rest
= scm_cdr (rest
);
1481 return scm_logxor (x
, y
);
1485 #define s_scm_logxor s_scm_i_logxor
1487 SCM
scm_logxor (SCM n1
, SCM n2
)
1488 #define FUNC_NAME s_scm_logxor
1492 if (SCM_UNBNDP (n2
))
1494 if (SCM_UNBNDP (n1
))
1496 else if (SCM_NUMBERP (n1
))
1499 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1502 if (SCM_I_INUMP (n1
))
1504 nn1
= SCM_I_INUM (n1
);
1505 if (SCM_I_INUMP (n2
))
1507 scm_t_inum nn2
= SCM_I_INUM (n2
);
1508 return SCM_I_MAKINUM (nn1
^ nn2
);
1510 else if (SCM_BIGP (n2
))
1514 SCM result_z
= scm_i_mkbig ();
1516 mpz_init_set_si (nn1_z
, nn1
);
1517 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1518 scm_remember_upto_here_1 (n2
);
1520 return scm_i_normbig (result_z
);
1524 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1526 else if (SCM_BIGP (n1
))
1528 if (SCM_I_INUMP (n2
))
1531 nn1
= SCM_I_INUM (n1
);
1534 else if (SCM_BIGP (n2
))
1536 SCM result_z
= scm_i_mkbig ();
1537 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1539 SCM_I_BIG_MPZ (n2
));
1540 scm_remember_upto_here_2 (n1
, n2
);
1541 return scm_i_normbig (result_z
);
1544 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1547 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1552 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1554 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1555 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1556 "without actually calculating the @code{logand}, just testing\n"
1560 "(logtest #b0100 #b1011) @result{} #f\n"
1561 "(logtest #b0100 #b0111) @result{} #t\n"
1563 #define FUNC_NAME s_scm_logtest
1567 if (SCM_I_INUMP (j
))
1569 nj
= SCM_I_INUM (j
);
1570 if (SCM_I_INUMP (k
))
1572 scm_t_inum nk
= SCM_I_INUM (k
);
1573 return scm_from_bool (nj
& nk
);
1575 else if (SCM_BIGP (k
))
1583 mpz_init_set_si (nj_z
, nj
);
1584 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1585 scm_remember_upto_here_1 (k
);
1586 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1592 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1594 else if (SCM_BIGP (j
))
1596 if (SCM_I_INUMP (k
))
1599 nj
= SCM_I_INUM (j
);
1602 else if (SCM_BIGP (k
))
1606 mpz_init (result_z
);
1610 scm_remember_upto_here_2 (j
, k
);
1611 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1612 mpz_clear (result_z
);
1616 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1619 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1624 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1626 "Test whether bit number @var{index} in @var{j} is set.\n"
1627 "@var{index} starts from 0 for the least significant bit.\n"
1630 "(logbit? 0 #b1101) @result{} #t\n"
1631 "(logbit? 1 #b1101) @result{} #f\n"
1632 "(logbit? 2 #b1101) @result{} #t\n"
1633 "(logbit? 3 #b1101) @result{} #t\n"
1634 "(logbit? 4 #b1101) @result{} #f\n"
1636 #define FUNC_NAME s_scm_logbit_p
1638 unsigned long int iindex
;
1639 iindex
= scm_to_ulong (index
);
1641 if (SCM_I_INUMP (j
))
1643 /* bits above what's in an inum follow the sign bit */
1644 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1645 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1647 else if (SCM_BIGP (j
))
1649 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1650 scm_remember_upto_here_1 (j
);
1651 return scm_from_bool (val
);
1654 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1659 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1661 "Return the integer which is the ones-complement of the integer\n"
1665 "(number->string (lognot #b10000000) 2)\n"
1666 " @result{} \"-10000001\"\n"
1667 "(number->string (lognot #b0) 2)\n"
1668 " @result{} \"-1\"\n"
1670 #define FUNC_NAME s_scm_lognot
1672 if (SCM_I_INUMP (n
)) {
1673 /* No overflow here, just need to toggle all the bits making up the inum.
1674 Enhancement: No need to strip the tag and add it back, could just xor
1675 a block of 1 bits, if that worked with the various debug versions of
1677 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1679 } else if (SCM_BIGP (n
)) {
1680 SCM result
= scm_i_mkbig ();
1681 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1682 scm_remember_upto_here_1 (n
);
1686 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1691 /* returns 0 if IN is not an integer. OUT must already be
1694 coerce_to_big (SCM in
, mpz_t out
)
1697 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1698 else if (SCM_I_INUMP (in
))
1699 mpz_set_si (out
, SCM_I_INUM (in
));
1706 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1707 (SCM n
, SCM k
, SCM m
),
1708 "Return @var{n} raised to the integer exponent\n"
1709 "@var{k}, modulo @var{m}.\n"
1712 "(modulo-expt 2 3 5)\n"
1715 #define FUNC_NAME s_scm_modulo_expt
1721 /* There are two classes of error we might encounter --
1722 1) Math errors, which we'll report by calling scm_num_overflow,
1724 2) wrong-type errors, which of course we'll report by calling
1726 We don't report those errors immediately, however; instead we do
1727 some cleanup first. These variables tell us which error (if
1728 any) we should report after cleaning up.
1730 int report_overflow
= 0;
1732 int position_of_wrong_type
= 0;
1733 SCM value_of_wrong_type
= SCM_INUM0
;
1735 SCM result
= SCM_UNDEFINED
;
1741 if (scm_is_eq (m
, SCM_INUM0
))
1743 report_overflow
= 1;
1747 if (!coerce_to_big (n
, n_tmp
))
1749 value_of_wrong_type
= n
;
1750 position_of_wrong_type
= 1;
1754 if (!coerce_to_big (k
, k_tmp
))
1756 value_of_wrong_type
= k
;
1757 position_of_wrong_type
= 2;
1761 if (!coerce_to_big (m
, m_tmp
))
1763 value_of_wrong_type
= m
;
1764 position_of_wrong_type
= 3;
1768 /* if the exponent K is negative, and we simply call mpz_powm, we
1769 will get a divide-by-zero exception when an inverse 1/n mod m
1770 doesn't exist (or is not unique). Since exceptions are hard to
1771 handle, we'll attempt the inversion "by hand" -- that way, we get
1772 a simple failure code, which is easy to handle. */
1774 if (-1 == mpz_sgn (k_tmp
))
1776 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1778 report_overflow
= 1;
1781 mpz_neg (k_tmp
, k_tmp
);
1784 result
= scm_i_mkbig ();
1785 mpz_powm (SCM_I_BIG_MPZ (result
),
1790 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1791 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1798 if (report_overflow
)
1799 scm_num_overflow (FUNC_NAME
);
1801 if (position_of_wrong_type
)
1802 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1803 value_of_wrong_type
);
1805 return scm_i_normbig (result
);
1809 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1811 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1812 "exact integer, @var{n} can be any number.\n"
1814 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1815 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1816 "includes @math{0^0} is 1.\n"
1819 "(integer-expt 2 5) @result{} 32\n"
1820 "(integer-expt -3 3) @result{} -27\n"
1821 "(integer-expt 5 -3) @result{} 1/125\n"
1822 "(integer-expt 0 0) @result{} 1\n"
1824 #define FUNC_NAME s_scm_integer_expt
1827 SCM z_i2
= SCM_BOOL_F
;
1829 SCM acc
= SCM_I_MAKINUM (1L);
1831 SCM_VALIDATE_NUMBER (SCM_ARG1
, n
);
1832 if (!SCM_I_INUMP (k
) && !SCM_BIGP (k
))
1833 SCM_WRONG_TYPE_ARG (2, k
);
1835 if (scm_is_true (scm_zero_p (n
)))
1837 if (scm_is_true (scm_zero_p (k
))) /* 0^0 == 1 per R5RS */
1838 return acc
; /* return exact 1, regardless of n */
1839 else if (scm_is_true (scm_positive_p (k
)))
1841 else /* return NaN for (0 ^ k) for negative k per R6RS */
1844 else if (scm_is_eq (n
, acc
))
1846 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1847 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1849 if (SCM_I_INUMP (k
))
1850 i2
= SCM_I_INUM (k
);
1851 else if (SCM_BIGP (k
))
1853 z_i2
= scm_i_clonebig (k
, 1);
1854 scm_remember_upto_here_1 (k
);
1858 SCM_WRONG_TYPE_ARG (2, k
);
1862 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1864 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1865 n
= scm_divide (n
, SCM_UNDEFINED
);
1869 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1873 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1875 return scm_product (acc
, n
);
1877 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1878 acc
= scm_product (acc
, n
);
1879 n
= scm_product (n
, n
);
1880 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1888 n
= scm_divide (n
, SCM_UNDEFINED
);
1895 return scm_product (acc
, n
);
1897 acc
= scm_product (acc
, n
);
1898 n
= scm_product (n
, n
);
1905 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1907 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1908 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1910 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1911 "@var{cnt} is negative it's a division, rounded towards negative\n"
1912 "infinity. (Note that this is not the same rounding as\n"
1913 "@code{quotient} does.)\n"
1915 "With @var{n} viewed as an infinite precision twos complement,\n"
1916 "@code{ash} means a left shift introducing zero bits, or a right\n"
1917 "shift dropping bits.\n"
1920 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1921 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1923 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1924 "(ash -23 -2) @result{} -6\n"
1926 #define FUNC_NAME s_scm_ash
1929 bits_to_shift
= scm_to_long (cnt
);
1931 if (SCM_I_INUMP (n
))
1933 scm_t_inum nn
= SCM_I_INUM (n
);
1935 if (bits_to_shift
> 0)
1937 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1938 overflow a non-zero fixnum. For smaller shifts we check the
1939 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1940 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1941 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1947 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1949 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1952 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1956 SCM result
= scm_i_inum2big (nn
);
1957 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1964 bits_to_shift
= -bits_to_shift
;
1965 if (bits_to_shift
>= SCM_LONG_BIT
)
1966 return (nn
>= 0 ? SCM_INUM0
: SCM_I_MAKINUM(-1));
1968 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1972 else if (SCM_BIGP (n
))
1976 if (bits_to_shift
== 0)
1979 result
= scm_i_mkbig ();
1980 if (bits_to_shift
>= 0)
1982 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1988 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1989 we have to allocate a bignum even if the result is going to be a
1991 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1993 return scm_i_normbig (result
);
1999 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2005 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
2006 (SCM n
, SCM start
, SCM end
),
2007 "Return the integer composed of the @var{start} (inclusive)\n"
2008 "through @var{end} (exclusive) bits of @var{n}. The\n"
2009 "@var{start}th bit becomes the 0-th bit in the result.\n"
2012 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
2013 " @result{} \"1010\"\n"
2014 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
2015 " @result{} \"10110\"\n"
2017 #define FUNC_NAME s_scm_bit_extract
2019 unsigned long int istart
, iend
, bits
;
2020 istart
= scm_to_ulong (start
);
2021 iend
= scm_to_ulong (end
);
2022 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
2024 /* how many bits to keep */
2025 bits
= iend
- istart
;
2027 if (SCM_I_INUMP (n
))
2029 scm_t_inum in
= SCM_I_INUM (n
);
2031 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
2032 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
2033 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
2035 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
2037 /* Since we emulate two's complement encoded numbers, this
2038 * special case requires us to produce a result that has
2039 * more bits than can be stored in a fixnum.
2041 SCM result
= scm_i_inum2big (in
);
2042 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
2047 /* mask down to requisite bits */
2048 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
2049 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
2051 else if (SCM_BIGP (n
))
2056 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
2060 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
2061 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
2062 such bits into a ulong. */
2063 result
= scm_i_mkbig ();
2064 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
2065 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
2066 result
= scm_i_normbig (result
);
2068 scm_remember_upto_here_1 (n
);
2072 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2077 static const char scm_logtab
[] = {
2078 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
2081 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
2083 "Return the number of bits in integer @var{n}. If integer is\n"
2084 "positive, the 1-bits in its binary representation are counted.\n"
2085 "If negative, the 0-bits in its two's-complement binary\n"
2086 "representation are counted. If 0, 0 is returned.\n"
2089 "(logcount #b10101010)\n"
2096 #define FUNC_NAME s_scm_logcount
2098 if (SCM_I_INUMP (n
))
2100 unsigned long c
= 0;
2101 scm_t_inum nn
= SCM_I_INUM (n
);
2106 c
+= scm_logtab
[15 & nn
];
2109 return SCM_I_MAKINUM (c
);
2111 else if (SCM_BIGP (n
))
2113 unsigned long count
;
2114 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
2115 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
2117 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
2118 scm_remember_upto_here_1 (n
);
2119 return SCM_I_MAKINUM (count
);
2122 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2127 static const char scm_ilentab
[] = {
2128 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2132 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2134 "Return the number of bits necessary to represent @var{n}.\n"
2137 "(integer-length #b10101010)\n"
2139 "(integer-length 0)\n"
2141 "(integer-length #b1111)\n"
2144 #define FUNC_NAME s_scm_integer_length
2146 if (SCM_I_INUMP (n
))
2148 unsigned long c
= 0;
2150 scm_t_inum nn
= SCM_I_INUM (n
);
2156 l
= scm_ilentab
[15 & nn
];
2159 return SCM_I_MAKINUM (c
- 4 + l
);
2161 else if (SCM_BIGP (n
))
2163 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2164 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2165 1 too big, so check for that and adjust. */
2166 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2167 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2168 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2169 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2171 scm_remember_upto_here_1 (n
);
2172 return SCM_I_MAKINUM (size
);
2175 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2179 /*** NUMBERS -> STRINGS ***/
2180 #define SCM_MAX_DBL_PREC 60
2181 #define SCM_MAX_DBL_RADIX 36
2183 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2184 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2185 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2188 void init_dblprec(int *prec
, int radix
) {
2189 /* determine floating point precision by adding successively
2190 smaller increments to 1.0 until it is considered == 1.0 */
2191 double f
= ((double)1.0)/radix
;
2192 double fsum
= 1.0 + f
;
2197 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2209 void init_fx_radix(double *fx_list
, int radix
)
2211 /* initialize a per-radix list of tolerances. When added
2212 to a number < 1.0, we can determine if we should raund
2213 up and quit converting a number to a string. */
2217 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2218 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2221 /* use this array as a way to generate a single digit */
2222 static const char number_chars
[] = "0123456789abcdefghijklmnopqrstuvwxyz";
2225 idbl2str (double f
, char *a
, int radix
)
2227 int efmt
, dpt
, d
, i
, wp
;
2229 #ifdef DBL_MIN_10_EXP
2232 #endif /* DBL_MIN_10_EXP */
2237 radix
> SCM_MAX_DBL_RADIX
)
2239 /* revert to existing behavior */
2243 wp
= scm_dblprec
[radix
-2];
2244 fx
= fx_per_radix
[radix
-2];
2248 #ifdef HAVE_COPYSIGN
2249 double sgn
= copysign (1.0, f
);
2254 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2260 strcpy (a
, "-inf.0");
2262 strcpy (a
, "+inf.0");
2267 strcpy (a
, "+nan.0");
2277 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2278 make-uniform-vector, from causing infinite loops. */
2279 /* just do the checking...if it passes, we do the conversion for our
2280 radix again below */
2287 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2295 while (f_cpy
> 10.0)
2298 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2319 if (f
+ fx
[wp
] >= radix
)
2326 /* adding 9999 makes this equivalent to abs(x) % 3 */
2327 dpt
= (exp
+ 9999) % 3;
2331 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2353 a
[ch
++] = number_chars
[d
];
2356 if (f
+ fx
[wp
] >= 1.0)
2358 a
[ch
- 1] = number_chars
[d
+1];
2370 if ((dpt
> 4) && (exp
> 6))
2372 d
= (a
[0] == '-' ? 2 : 1);
2373 for (i
= ch
++; i
> d
; i
--)
2386 if (a
[ch
- 1] == '.')
2387 a
[ch
++] = '0'; /* trailing zero */
2396 for (i
= radix
; i
<= exp
; i
*= radix
);
2397 for (i
/= radix
; i
; i
/= radix
)
2399 a
[ch
++] = number_chars
[exp
/ i
];
2408 icmplx2str (double real
, double imag
, char *str
, int radix
)
2412 i
= idbl2str (real
, str
, radix
);
2415 /* Don't output a '+' for negative numbers or for Inf and
2416 NaN. They will provide their own sign. */
2417 if (0 <= imag
&& !isinf (imag
) && !isnan (imag
))
2419 i
+= idbl2str (imag
, &str
[i
], radix
);
2426 iflo2str (SCM flt
, char *str
, int radix
)
2429 if (SCM_REALP (flt
))
2430 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2432 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2437 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2438 characters in the result.
2440 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2442 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2447 return scm_iuint2str (-num
, rad
, p
) + 1;
2450 return scm_iuint2str (num
, rad
, p
);
2453 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2454 characters in the result.
2456 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2458 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2462 scm_t_uintmax n
= num
;
2464 if (rad
< 2 || rad
> 36)
2465 scm_out_of_range ("scm_iuint2str", scm_from_int (rad
));
2467 for (n
/= rad
; n
> 0; n
/= rad
)
2477 p
[i
] = number_chars
[d
];
2482 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2484 "Return a string holding the external representation of the\n"
2485 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2486 "inexact, a radix of 10 will be used.")
2487 #define FUNC_NAME s_scm_number_to_string
2491 if (SCM_UNBNDP (radix
))
2494 base
= scm_to_signed_integer (radix
, 2, 36);
2496 if (SCM_I_INUMP (n
))
2498 char num_buf
[SCM_INTBUFLEN
];
2499 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2500 return scm_from_locale_stringn (num_buf
, length
);
2502 else if (SCM_BIGP (n
))
2504 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2505 scm_remember_upto_here_1 (n
);
2506 return scm_take_locale_string (str
);
2508 else if (SCM_FRACTIONP (n
))
2510 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2511 scm_from_locale_string ("/"),
2512 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2514 else if (SCM_INEXACTP (n
))
2516 char num_buf
[FLOBUFLEN
];
2517 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2520 SCM_WRONG_TYPE_ARG (1, n
);
2525 /* These print routines used to be stubbed here so that scm_repl.c
2526 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2529 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2531 char num_buf
[FLOBUFLEN
];
2532 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2537 scm_i_print_double (double val
, SCM port
)
2539 char num_buf
[FLOBUFLEN
];
2540 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2544 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2547 char num_buf
[FLOBUFLEN
];
2548 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2553 scm_i_print_complex (double real
, double imag
, SCM port
)
2555 char num_buf
[FLOBUFLEN
];
2556 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2560 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2563 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2564 scm_display (str
, port
);
2565 scm_remember_upto_here_1 (str
);
2570 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2572 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2573 scm_remember_upto_here_1 (exp
);
2574 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2578 /*** END nums->strs ***/
2581 /*** STRINGS -> NUMBERS ***/
2583 /* The following functions implement the conversion from strings to numbers.
2584 * The implementation somehow follows the grammar for numbers as it is given
2585 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2586 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2587 * points should be noted about the implementation:
2588 * * Each function keeps a local index variable 'idx' that points at the
2589 * current position within the parsed string. The global index is only
2590 * updated if the function could parse the corresponding syntactic unit
2592 * * Similarly, the functions keep track of indicators of inexactness ('#',
2593 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2594 * global exactness information is only updated after each part has been
2595 * successfully parsed.
2596 * * Sequences of digits are parsed into temporary variables holding fixnums.
2597 * Only if these fixnums would overflow, the result variables are updated
2598 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2599 * the temporary variables holding the fixnums are cleared, and the process
2600 * starts over again. If for example fixnums were able to store five decimal
2601 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2602 * and the result was computed as 12345 * 100000 + 67890. In other words,
2603 * only every five digits two bignum operations were performed.
2606 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2608 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2610 /* Caller is responsible for checking that the return value is in range
2611 for the given radix, which should be <= 36. */
2613 char_decimal_value (scm_t_uint32 c
)
2615 /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
2616 that's certainly above any valid decimal, so we take advantage of
2617 that to elide some tests. */
2618 unsigned int d
= (unsigned int) uc_decimal_value (c
);
2620 /* If that failed, try extended hexadecimals, then. Only accept ascii
2625 if (c
>= (scm_t_uint32
) 'a')
2626 d
= c
- (scm_t_uint32
)'a' + 10U;
2632 mem2uinteger (SCM mem
, unsigned int *p_idx
,
2633 unsigned int radix
, enum t_exactness
*p_exactness
)
2635 unsigned int idx
= *p_idx
;
2636 unsigned int hash_seen
= 0;
2637 scm_t_bits shift
= 1;
2639 unsigned int digit_value
;
2642 size_t len
= scm_i_string_length (mem
);
2647 c
= scm_i_string_ref (mem
, idx
);
2648 digit_value
= char_decimal_value (c
);
2649 if (digit_value
>= radix
)
2653 result
= SCM_I_MAKINUM (digit_value
);
2656 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2666 digit_value
= char_decimal_value (c
);
2667 /* This check catches non-decimals in addition to out-of-range
2669 if (digit_value
>= radix
)
2674 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2676 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2678 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2685 shift
= shift
* radix
;
2686 add
= add
* radix
+ digit_value
;
2691 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2693 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2697 *p_exactness
= INEXACT
;
2703 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2704 * covers the parts of the rules that start at a potential point. The value
2705 * of the digits up to the point have been parsed by the caller and are given
2706 * in variable result. The content of *p_exactness indicates, whether a hash
2707 * has already been seen in the digits before the point.
2710 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
2713 mem2decimal_from_point (SCM result
, SCM mem
,
2714 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2716 unsigned int idx
= *p_idx
;
2717 enum t_exactness x
= *p_exactness
;
2718 size_t len
= scm_i_string_length (mem
);
2723 if (scm_i_string_ref (mem
, idx
) == '.')
2725 scm_t_bits shift
= 1;
2727 unsigned int digit_value
;
2728 SCM big_shift
= SCM_INUM1
;
2733 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2734 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2739 digit_value
= DIGIT2UINT (c
);
2750 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2752 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2753 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2755 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2763 add
= add
* 10 + digit_value
;
2769 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2770 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2771 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2774 result
= scm_divide (result
, big_shift
);
2776 /* We've seen a decimal point, thus the value is implicitly inexact. */
2788 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2790 switch (scm_i_string_ref (mem
, idx
))
2802 c
= scm_i_string_ref (mem
, idx
);
2810 c
= scm_i_string_ref (mem
, idx
);
2819 c
= scm_i_string_ref (mem
, idx
);
2824 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2828 exponent
= DIGIT2UINT (c
);
2831 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2832 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2835 if (exponent
<= SCM_MAXEXP
)
2836 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2842 if (exponent
> SCM_MAXEXP
)
2844 size_t exp_len
= idx
- start
;
2845 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
2846 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2847 scm_out_of_range ("string->number", exp_num
);
2850 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2852 result
= scm_product (result
, e
);
2854 result
= scm_divide2real (result
, e
);
2856 /* We've seen an exponent, thus the value is implicitly inexact. */
2874 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2877 mem2ureal (SCM mem
, unsigned int *p_idx
,
2878 unsigned int radix
, enum t_exactness
*p_exactness
)
2880 unsigned int idx
= *p_idx
;
2882 size_t len
= scm_i_string_length (mem
);
2884 /* Start off believing that the number will be exact. This changes
2885 to INEXACT if we see a decimal point or a hash. */
2886 enum t_exactness x
= EXACT
;
2891 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
2897 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
2899 /* Cobble up the fractional part. We might want to set the
2900 NaN's mantissa from it. */
2902 mem2uinteger (mem
, &idx
, 10, &x
);
2907 if (scm_i_string_ref (mem
, idx
) == '.')
2911 else if (idx
+ 1 == len
)
2913 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
2916 result
= mem2decimal_from_point (SCM_INUM0
, mem
,
2923 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
2924 if (scm_is_false (uinteger
))
2929 else if (scm_i_string_ref (mem
, idx
) == '/')
2937 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
2938 if (scm_is_false (divisor
))
2941 /* both are int/big here, I assume */
2942 result
= scm_i_make_ratio (uinteger
, divisor
);
2944 else if (radix
== 10)
2946 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
2947 if (scm_is_false (result
))
2956 /* Update *p_exactness if the number just read was inexact. This is
2957 important for complex numbers, so that a complex number is
2958 treated as inexact overall if either its real or imaginary part
2964 /* When returning an inexact zero, make sure it is represented as a
2965 floating point value so that we can change its sign.
2967 if (scm_is_eq (result
, SCM_INUM0
) && *p_exactness
== INEXACT
)
2968 result
= scm_from_double (0.0);
2974 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2977 mem2complex (SCM mem
, unsigned int idx
,
2978 unsigned int radix
, enum t_exactness
*p_exactness
)
2983 size_t len
= scm_i_string_length (mem
);
2988 c
= scm_i_string_ref (mem
, idx
);
3003 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3004 if (scm_is_false (ureal
))
3006 /* input must be either +i or -i */
3011 if (scm_i_string_ref (mem
, idx
) == 'i'
3012 || scm_i_string_ref (mem
, idx
) == 'I')
3018 return scm_make_rectangular (SCM_INUM0
, SCM_I_MAKINUM (sign
));
3025 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3026 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
3031 c
= scm_i_string_ref (mem
, idx
);
3035 /* either +<ureal>i or -<ureal>i */
3042 return scm_make_rectangular (SCM_INUM0
, ureal
);
3045 /* polar input: <real>@<real>. */
3056 c
= scm_i_string_ref (mem
, idx
);
3074 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3075 if (scm_is_false (angle
))
3080 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3081 angle
= scm_difference (angle
, SCM_UNDEFINED
);
3083 result
= scm_make_polar (ureal
, angle
);
3088 /* expecting input matching <real>[+-]<ureal>?i */
3095 int sign
= (c
== '+') ? 1 : -1;
3096 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3098 if (scm_is_false (imag
))
3099 imag
= SCM_I_MAKINUM (sign
);
3100 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
3101 imag
= scm_difference (imag
, SCM_UNDEFINED
);
3105 if (scm_i_string_ref (mem
, idx
) != 'i'
3106 && scm_i_string_ref (mem
, idx
) != 'I')
3113 return scm_make_rectangular (ureal
, imag
);
3122 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
3124 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
3127 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
3129 unsigned int idx
= 0;
3130 unsigned int radix
= NO_RADIX
;
3131 enum t_exactness forced_x
= NO_EXACTNESS
;
3132 enum t_exactness implicit_x
= EXACT
;
3134 size_t len
= scm_i_string_length (mem
);
3136 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3137 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
3139 switch (scm_i_string_ref (mem
, idx
+ 1))
3142 if (radix
!= NO_RADIX
)
3147 if (radix
!= NO_RADIX
)
3152 if (forced_x
!= NO_EXACTNESS
)
3157 if (forced_x
!= NO_EXACTNESS
)
3162 if (radix
!= NO_RADIX
)
3167 if (radix
!= NO_RADIX
)
3177 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3178 if (radix
== NO_RADIX
)
3179 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
3181 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
3183 if (scm_is_false (result
))
3189 if (SCM_INEXACTP (result
))
3190 return scm_inexact_to_exact (result
);
3194 if (SCM_INEXACTP (result
))
3197 return scm_exact_to_inexact (result
);
3200 if (implicit_x
== INEXACT
)
3202 if (SCM_INEXACTP (result
))
3205 return scm_exact_to_inexact (result
);
3213 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
3214 unsigned int default_radix
)
3216 SCM str
= scm_from_locale_stringn (mem
, len
);
3218 return scm_i_string_to_number (str
, default_radix
);
3222 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3223 (SCM string
, SCM radix
),
3224 "Return a number of the maximally precise representation\n"
3225 "expressed by the given @var{string}. @var{radix} must be an\n"
3226 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3227 "is a default radix that may be overridden by an explicit radix\n"
3228 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3229 "supplied, then the default radix is 10. If string is not a\n"
3230 "syntactically valid notation for a number, then\n"
3231 "@code{string->number} returns @code{#f}.")
3232 #define FUNC_NAME s_scm_string_to_number
3236 SCM_VALIDATE_STRING (1, string
);
3238 if (SCM_UNBNDP (radix
))
3241 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3243 answer
= scm_i_string_to_number (string
, base
);
3244 scm_remember_upto_here_1 (string
);
3250 /*** END strs->nums ***/
3253 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3255 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3257 #define FUNC_NAME s_scm_number_p
3259 return scm_from_bool (SCM_NUMBERP (x
));
3263 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3265 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3266 "otherwise. Note that the sets of real, rational and integer\n"
3267 "values form subsets of the set of complex numbers, i. e. the\n"
3268 "predicate will also be fulfilled if @var{x} is a real,\n"
3269 "rational or integer number.")
3270 #define FUNC_NAME s_scm_complex_p
3272 /* all numbers are complex. */
3273 return scm_number_p (x
);
3277 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3279 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3280 "otherwise. Note that the set of integer values forms a subset of\n"
3281 "the set of real numbers, i. e. the predicate will also be\n"
3282 "fulfilled if @var{x} is an integer number.")
3283 #define FUNC_NAME s_scm_real_p
3285 return scm_from_bool
3286 (SCM_I_INUMP (x
) || SCM_REALP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
));
3290 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3292 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3293 "otherwise. Note that the set of integer values forms a subset of\n"
3294 "the set of rational numbers, i. e. the predicate will also be\n"
3295 "fulfilled if @var{x} is an integer number.")
3296 #define FUNC_NAME s_scm_rational_p
3298 if (SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
))
3300 else if (SCM_REALP (x
))
3301 /* due to their limited precision, finite floating point numbers are
3302 rational as well. (finite means neither infinity nor a NaN) */
3303 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
3309 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3311 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3313 #define FUNC_NAME s_scm_integer_p
3315 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
3317 else if (SCM_REALP (x
))
3319 double val
= SCM_REAL_VALUE (x
);
3320 return scm_from_bool (!isinf (val
) && (val
== floor (val
)));
3328 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
3329 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
3330 (SCM x
, SCM y
, SCM rest
),
3331 "Return @code{#t} if all parameters are numerically equal.")
3332 #define FUNC_NAME s_scm_i_num_eq_p
3334 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3336 while (!scm_is_null (rest
))
3338 if (scm_is_false (scm_num_eq_p (x
, y
)))
3342 rest
= scm_cdr (rest
);
3344 return scm_num_eq_p (x
, y
);
3348 scm_num_eq_p (SCM x
, SCM y
)
3351 if (SCM_I_INUMP (x
))
3353 scm_t_signed_bits xx
= SCM_I_INUM (x
);
3354 if (SCM_I_INUMP (y
))
3356 scm_t_signed_bits yy
= SCM_I_INUM (y
);
3357 return scm_from_bool (xx
== yy
);
3359 else if (SCM_BIGP (y
))
3361 else if (SCM_REALP (y
))
3363 /* On a 32-bit system an inum fits a double, we can cast the inum
3364 to a double and compare.
3366 But on a 64-bit system an inum is bigger than a double and
3367 casting it to a double (call that dxx) will round. dxx is at
3368 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3369 an integer and fits a long. So we cast yy to a long and
3370 compare with plain xx.
3372 An alternative (for any size system actually) would be to check
3373 yy is an integer (with floor) and is in range of an inum
3374 (compare against appropriate powers of 2) then test
3375 xx==(scm_t_signed_bits)yy. It's just a matter of which
3376 casts/comparisons might be fastest or easiest for the cpu. */
3378 double yy
= SCM_REAL_VALUE (y
);
3379 return scm_from_bool ((double) xx
== yy
3380 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3381 || xx
== (scm_t_signed_bits
) yy
));
3383 else if (SCM_COMPLEXP (y
))
3384 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3385 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3386 else if (SCM_FRACTIONP (y
))
3389 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3391 else if (SCM_BIGP (x
))
3393 if (SCM_I_INUMP (y
))
3395 else if (SCM_BIGP (y
))
3397 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3398 scm_remember_upto_here_2 (x
, y
);
3399 return scm_from_bool (0 == cmp
);
3401 else if (SCM_REALP (y
))
3404 if (isnan (SCM_REAL_VALUE (y
)))
3406 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3407 scm_remember_upto_here_1 (x
);
3408 return scm_from_bool (0 == cmp
);
3410 else if (SCM_COMPLEXP (y
))
3413 if (0.0 != SCM_COMPLEX_IMAG (y
))
3415 if (isnan (SCM_COMPLEX_REAL (y
)))
3417 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3418 scm_remember_upto_here_1 (x
);
3419 return scm_from_bool (0 == cmp
);
3421 else if (SCM_FRACTIONP (y
))
3424 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3426 else if (SCM_REALP (x
))
3428 double xx
= SCM_REAL_VALUE (x
);
3429 if (SCM_I_INUMP (y
))
3431 /* see comments with inum/real above */
3432 scm_t_signed_bits yy
= SCM_I_INUM (y
);
3433 return scm_from_bool (xx
== (double) yy
3434 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3435 || (scm_t_signed_bits
) xx
== yy
));
3437 else if (SCM_BIGP (y
))
3440 if (isnan (SCM_REAL_VALUE (x
)))
3442 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3443 scm_remember_upto_here_1 (y
);
3444 return scm_from_bool (0 == cmp
);
3446 else if (SCM_REALP (y
))
3447 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3448 else if (SCM_COMPLEXP (y
))
3449 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3450 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3451 else if (SCM_FRACTIONP (y
))
3453 double xx
= SCM_REAL_VALUE (x
);
3457 return scm_from_bool (xx
< 0.0);
3458 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3462 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3464 else if (SCM_COMPLEXP (x
))
3466 if (SCM_I_INUMP (y
))
3467 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3468 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3469 else if (SCM_BIGP (y
))
3472 if (0.0 != SCM_COMPLEX_IMAG (x
))
3474 if (isnan (SCM_COMPLEX_REAL (x
)))
3476 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3477 scm_remember_upto_here_1 (y
);
3478 return scm_from_bool (0 == cmp
);
3480 else if (SCM_REALP (y
))
3481 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3482 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3483 else if (SCM_COMPLEXP (y
))
3484 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3485 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3486 else if (SCM_FRACTIONP (y
))
3489 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3491 xx
= SCM_COMPLEX_REAL (x
);
3495 return scm_from_bool (xx
< 0.0);
3496 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3500 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3502 else if (SCM_FRACTIONP (x
))
3504 if (SCM_I_INUMP (y
))
3506 else if (SCM_BIGP (y
))
3508 else if (SCM_REALP (y
))
3510 double yy
= SCM_REAL_VALUE (y
);
3514 return scm_from_bool (0.0 < yy
);
3515 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3518 else if (SCM_COMPLEXP (y
))
3521 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3523 yy
= SCM_COMPLEX_REAL (y
);
3527 return scm_from_bool (0.0 < yy
);
3528 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3531 else if (SCM_FRACTIONP (y
))
3532 return scm_i_fraction_equalp (x
, y
);
3534 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3537 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
3541 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3542 done are good for inums, but for bignums an answer can almost always be
3543 had by just examining a few high bits of the operands, as done by GMP in
3544 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3545 of the float exponent to take into account. */
3547 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
3548 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
3549 (SCM x
, SCM y
, SCM rest
),
3550 "Return @code{#t} if the list of parameters is monotonically\n"
3552 #define FUNC_NAME s_scm_i_num_less_p
3554 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3556 while (!scm_is_null (rest
))
3558 if (scm_is_false (scm_less_p (x
, y
)))
3562 rest
= scm_cdr (rest
);
3564 return scm_less_p (x
, y
);
3568 scm_less_p (SCM x
, SCM y
)
3571 if (SCM_I_INUMP (x
))
3573 scm_t_inum xx
= SCM_I_INUM (x
);
3574 if (SCM_I_INUMP (y
))
3576 scm_t_inum yy
= SCM_I_INUM (y
);
3577 return scm_from_bool (xx
< yy
);
3579 else if (SCM_BIGP (y
))
3581 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3582 scm_remember_upto_here_1 (y
);
3583 return scm_from_bool (sgn
> 0);
3585 else if (SCM_REALP (y
))
3586 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3587 else if (SCM_FRACTIONP (y
))
3589 /* "x < a/b" becomes "x*b < a" */
3591 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3592 y
= SCM_FRACTION_NUMERATOR (y
);
3596 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3598 else if (SCM_BIGP (x
))
3600 if (SCM_I_INUMP (y
))
3602 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3603 scm_remember_upto_here_1 (x
);
3604 return scm_from_bool (sgn
< 0);
3606 else if (SCM_BIGP (y
))
3608 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3609 scm_remember_upto_here_2 (x
, y
);
3610 return scm_from_bool (cmp
< 0);
3612 else if (SCM_REALP (y
))
3615 if (isnan (SCM_REAL_VALUE (y
)))
3617 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3618 scm_remember_upto_here_1 (x
);
3619 return scm_from_bool (cmp
< 0);
3621 else if (SCM_FRACTIONP (y
))
3624 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3626 else if (SCM_REALP (x
))
3628 if (SCM_I_INUMP (y
))
3629 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3630 else if (SCM_BIGP (y
))
3633 if (isnan (SCM_REAL_VALUE (x
)))
3635 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3636 scm_remember_upto_here_1 (y
);
3637 return scm_from_bool (cmp
> 0);
3639 else if (SCM_REALP (y
))
3640 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3641 else if (SCM_FRACTIONP (y
))
3643 double xx
= SCM_REAL_VALUE (x
);
3647 return scm_from_bool (xx
< 0.0);
3648 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3652 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3654 else if (SCM_FRACTIONP (x
))
3656 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3658 /* "a/b < y" becomes "a < y*b" */
3659 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3660 x
= SCM_FRACTION_NUMERATOR (x
);
3663 else if (SCM_REALP (y
))
3665 double yy
= SCM_REAL_VALUE (y
);
3669 return scm_from_bool (0.0 < yy
);
3670 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3673 else if (SCM_FRACTIONP (y
))
3675 /* "a/b < c/d" becomes "a*d < c*b" */
3676 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3677 SCM_FRACTION_DENOMINATOR (y
));
3678 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3679 SCM_FRACTION_DENOMINATOR (x
));
3685 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3688 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
3692 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
3693 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
3694 (SCM x
, SCM y
, SCM rest
),
3695 "Return @code{#t} if the list of parameters is monotonically\n"
3697 #define FUNC_NAME s_scm_i_num_gr_p
3699 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3701 while (!scm_is_null (rest
))
3703 if (scm_is_false (scm_gr_p (x
, y
)))
3707 rest
= scm_cdr (rest
);
3709 return scm_gr_p (x
, y
);
3712 #define FUNC_NAME s_scm_i_num_gr_p
3714 scm_gr_p (SCM x
, SCM y
)
3716 if (!SCM_NUMBERP (x
))
3717 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3718 else if (!SCM_NUMBERP (y
))
3719 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3721 return scm_less_p (y
, x
);
3726 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
3727 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
3728 (SCM x
, SCM y
, SCM rest
),
3729 "Return @code{#t} if the list of parameters is monotonically\n"
3731 #define FUNC_NAME s_scm_i_num_leq_p
3733 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3735 while (!scm_is_null (rest
))
3737 if (scm_is_false (scm_leq_p (x
, y
)))
3741 rest
= scm_cdr (rest
);
3743 return scm_leq_p (x
, y
);
3746 #define FUNC_NAME s_scm_i_num_leq_p
3748 scm_leq_p (SCM x
, SCM y
)
3750 if (!SCM_NUMBERP (x
))
3751 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3752 else if (!SCM_NUMBERP (y
))
3753 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3754 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3757 return scm_not (scm_less_p (y
, x
));
3762 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
3763 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
3764 (SCM x
, SCM y
, SCM rest
),
3765 "Return @code{#t} if the list of parameters is monotonically\n"
3767 #define FUNC_NAME s_scm_i_num_geq_p
3769 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3771 while (!scm_is_null (rest
))
3773 if (scm_is_false (scm_geq_p (x
, y
)))
3777 rest
= scm_cdr (rest
);
3779 return scm_geq_p (x
, y
);
3782 #define FUNC_NAME s_scm_i_num_geq_p
3784 scm_geq_p (SCM x
, SCM y
)
3786 if (!SCM_NUMBERP (x
))
3787 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3788 else if (!SCM_NUMBERP (y
))
3789 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3790 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3793 return scm_not (scm_less_p (x
, y
));
3798 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3799 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3805 if (SCM_I_INUMP (z
))
3806 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3807 else if (SCM_BIGP (z
))
3809 else if (SCM_REALP (z
))
3810 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3811 else if (SCM_COMPLEXP (z
))
3812 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3813 && SCM_COMPLEX_IMAG (z
) == 0.0);
3814 else if (SCM_FRACTIONP (z
))
3817 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3821 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3822 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3826 scm_positive_p (SCM x
)
3828 if (SCM_I_INUMP (x
))
3829 return scm_from_bool (SCM_I_INUM (x
) > 0);
3830 else if (SCM_BIGP (x
))
3832 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3833 scm_remember_upto_here_1 (x
);
3834 return scm_from_bool (sgn
> 0);
3836 else if (SCM_REALP (x
))
3837 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3838 else if (SCM_FRACTIONP (x
))
3839 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3841 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3845 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3846 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3850 scm_negative_p (SCM x
)
3852 if (SCM_I_INUMP (x
))
3853 return scm_from_bool (SCM_I_INUM (x
) < 0);
3854 else if (SCM_BIGP (x
))
3856 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3857 scm_remember_upto_here_1 (x
);
3858 return scm_from_bool (sgn
< 0);
3860 else if (SCM_REALP (x
))
3861 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3862 else if (SCM_FRACTIONP (x
))
3863 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3865 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3869 /* scm_min and scm_max return an inexact when either argument is inexact, as
3870 required by r5rs. On that basis, for exact/inexact combinations the
3871 exact is converted to inexact to compare and possibly return. This is
3872 unlike scm_less_p above which takes some trouble to preserve all bits in
3873 its test, such trouble is not required for min and max. */
3875 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
3876 (SCM x
, SCM y
, SCM rest
),
3877 "Return the maximum of all parameter values.")
3878 #define FUNC_NAME s_scm_i_max
3880 while (!scm_is_null (rest
))
3881 { x
= scm_max (x
, y
);
3883 rest
= scm_cdr (rest
);
3885 return scm_max (x
, y
);
3889 #define s_max s_scm_i_max
3890 #define g_max g_scm_i_max
3893 scm_max (SCM x
, SCM y
)
3898 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3899 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3902 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3905 if (SCM_I_INUMP (x
))
3907 scm_t_inum xx
= SCM_I_INUM (x
);
3908 if (SCM_I_INUMP (y
))
3910 scm_t_inum yy
= SCM_I_INUM (y
);
3911 return (xx
< yy
) ? y
: x
;
3913 else if (SCM_BIGP (y
))
3915 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3916 scm_remember_upto_here_1 (y
);
3917 return (sgn
< 0) ? x
: y
;
3919 else if (SCM_REALP (y
))
3922 /* if y==NaN then ">" is false and we return NaN */
3923 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3925 else if (SCM_FRACTIONP (y
))
3928 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3931 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3933 else if (SCM_BIGP (x
))
3935 if (SCM_I_INUMP (y
))
3937 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3938 scm_remember_upto_here_1 (x
);
3939 return (sgn
< 0) ? y
: x
;
3941 else if (SCM_BIGP (y
))
3943 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3944 scm_remember_upto_here_2 (x
, y
);
3945 return (cmp
> 0) ? x
: y
;
3947 else if (SCM_REALP (y
))
3949 /* if y==NaN then xx>yy is false, so we return the NaN y */
3952 xx
= scm_i_big2dbl (x
);
3953 yy
= SCM_REAL_VALUE (y
);
3954 return (xx
> yy
? scm_from_double (xx
) : y
);
3956 else if (SCM_FRACTIONP (y
))
3961 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3963 else if (SCM_REALP (x
))
3965 if (SCM_I_INUMP (y
))
3967 double z
= SCM_I_INUM (y
);
3968 /* if x==NaN then "<" is false and we return NaN */
3969 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3971 else if (SCM_BIGP (y
))
3976 else if (SCM_REALP (y
))
3978 /* if x==NaN then our explicit check means we return NaN
3979 if y==NaN then ">" is false and we return NaN
3980 calling isnan is unavoidable, since it's the only way to know
3981 which of x or y causes any compares to be false */
3982 double xx
= SCM_REAL_VALUE (x
);
3983 return (isnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3985 else if (SCM_FRACTIONP (y
))
3987 double yy
= scm_i_fraction2double (y
);
3988 double xx
= SCM_REAL_VALUE (x
);
3989 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3992 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3994 else if (SCM_FRACTIONP (x
))
3996 if (SCM_I_INUMP (y
))
4000 else if (SCM_BIGP (y
))
4004 else if (SCM_REALP (y
))
4006 double xx
= scm_i_fraction2double (x
);
4007 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
4009 else if (SCM_FRACTIONP (y
))
4014 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4017 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
4021 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
4022 (SCM x
, SCM y
, SCM rest
),
4023 "Return the minimum of all parameter values.")
4024 #define FUNC_NAME s_scm_i_min
4026 while (!scm_is_null (rest
))
4027 { x
= scm_min (x
, y
);
4029 rest
= scm_cdr (rest
);
4031 return scm_min (x
, y
);
4035 #define s_min s_scm_i_min
4036 #define g_min g_scm_i_min
4039 scm_min (SCM x
, SCM y
)
4044 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
4045 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
4048 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
4051 if (SCM_I_INUMP (x
))
4053 scm_t_inum xx
= SCM_I_INUM (x
);
4054 if (SCM_I_INUMP (y
))
4056 scm_t_inum yy
= SCM_I_INUM (y
);
4057 return (xx
< yy
) ? x
: y
;
4059 else if (SCM_BIGP (y
))
4061 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4062 scm_remember_upto_here_1 (y
);
4063 return (sgn
< 0) ? y
: x
;
4065 else if (SCM_REALP (y
))
4068 /* if y==NaN then "<" is false and we return NaN */
4069 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
4071 else if (SCM_FRACTIONP (y
))
4074 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
4077 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4079 else if (SCM_BIGP (x
))
4081 if (SCM_I_INUMP (y
))
4083 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4084 scm_remember_upto_here_1 (x
);
4085 return (sgn
< 0) ? x
: y
;
4087 else if (SCM_BIGP (y
))
4089 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4090 scm_remember_upto_here_2 (x
, y
);
4091 return (cmp
> 0) ? y
: x
;
4093 else if (SCM_REALP (y
))
4095 /* if y==NaN then xx<yy is false, so we return the NaN y */
4098 xx
= scm_i_big2dbl (x
);
4099 yy
= SCM_REAL_VALUE (y
);
4100 return (xx
< yy
? scm_from_double (xx
) : y
);
4102 else if (SCM_FRACTIONP (y
))
4107 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4109 else if (SCM_REALP (x
))
4111 if (SCM_I_INUMP (y
))
4113 double z
= SCM_I_INUM (y
);
4114 /* if x==NaN then "<" is false and we return NaN */
4115 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
4117 else if (SCM_BIGP (y
))
4122 else if (SCM_REALP (y
))
4124 /* if x==NaN then our explicit check means we return NaN
4125 if y==NaN then "<" is false and we return NaN
4126 calling isnan is unavoidable, since it's the only way to know
4127 which of x or y causes any compares to be false */
4128 double xx
= SCM_REAL_VALUE (x
);
4129 return (isnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
4131 else if (SCM_FRACTIONP (y
))
4133 double yy
= scm_i_fraction2double (y
);
4134 double xx
= SCM_REAL_VALUE (x
);
4135 return (yy
< xx
) ? scm_from_double (yy
) : x
;
4138 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4140 else if (SCM_FRACTIONP (x
))
4142 if (SCM_I_INUMP (y
))
4146 else if (SCM_BIGP (y
))
4150 else if (SCM_REALP (y
))
4152 double xx
= scm_i_fraction2double (x
);
4153 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
4155 else if (SCM_FRACTIONP (y
))
4160 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4163 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
4167 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
4168 (SCM x
, SCM y
, SCM rest
),
4169 "Return the sum of all parameter values. Return 0 if called without\n"
4171 #define FUNC_NAME s_scm_i_sum
4173 while (!scm_is_null (rest
))
4174 { x
= scm_sum (x
, y
);
4176 rest
= scm_cdr (rest
);
4178 return scm_sum (x
, y
);
4182 #define s_sum s_scm_i_sum
4183 #define g_sum g_scm_i_sum
4186 scm_sum (SCM x
, SCM y
)
4188 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4190 if (SCM_NUMBERP (x
)) return x
;
4191 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4192 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4195 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4197 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4199 scm_t_inum xx
= SCM_I_INUM (x
);
4200 scm_t_inum yy
= SCM_I_INUM (y
);
4201 scm_t_inum z
= xx
+ yy
;
4202 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_inum2big (z
);
4204 else if (SCM_BIGP (y
))
4209 else if (SCM_REALP (y
))
4211 scm_t_inum xx
= SCM_I_INUM (x
);
4212 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4214 else if (SCM_COMPLEXP (y
))
4216 scm_t_inum xx
= SCM_I_INUM (x
);
4217 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4218 SCM_COMPLEX_IMAG (y
));
4220 else if (SCM_FRACTIONP (y
))
4221 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4222 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4223 SCM_FRACTION_DENOMINATOR (y
));
4225 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4226 } else if (SCM_BIGP (x
))
4228 if (SCM_I_INUMP (y
))
4233 inum
= SCM_I_INUM (y
);
4236 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4239 SCM result
= scm_i_mkbig ();
4240 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4241 scm_remember_upto_here_1 (x
);
4242 /* we know the result will have to be a bignum */
4245 return scm_i_normbig (result
);
4249 SCM result
= scm_i_mkbig ();
4250 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4251 scm_remember_upto_here_1 (x
);
4252 /* we know the result will have to be a bignum */
4255 return scm_i_normbig (result
);
4258 else if (SCM_BIGP (y
))
4260 SCM result
= scm_i_mkbig ();
4261 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4262 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4263 mpz_add (SCM_I_BIG_MPZ (result
),
4266 scm_remember_upto_here_2 (x
, y
);
4267 /* we know the result will have to be a bignum */
4270 return scm_i_normbig (result
);
4272 else if (SCM_REALP (y
))
4274 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4275 scm_remember_upto_here_1 (x
);
4276 return scm_from_double (result
);
4278 else if (SCM_COMPLEXP (y
))
4280 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4281 + SCM_COMPLEX_REAL (y
));
4282 scm_remember_upto_here_1 (x
);
4283 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4285 else if (SCM_FRACTIONP (y
))
4286 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4287 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4288 SCM_FRACTION_DENOMINATOR (y
));
4290 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4292 else if (SCM_REALP (x
))
4294 if (SCM_I_INUMP (y
))
4295 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4296 else if (SCM_BIGP (y
))
4298 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4299 scm_remember_upto_here_1 (y
);
4300 return scm_from_double (result
);
4302 else if (SCM_REALP (y
))
4303 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4304 else if (SCM_COMPLEXP (y
))
4305 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4306 SCM_COMPLEX_IMAG (y
));
4307 else if (SCM_FRACTIONP (y
))
4308 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4310 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4312 else if (SCM_COMPLEXP (x
))
4314 if (SCM_I_INUMP (y
))
4315 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4316 SCM_COMPLEX_IMAG (x
));
4317 else if (SCM_BIGP (y
))
4319 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4320 + SCM_COMPLEX_REAL (x
));
4321 scm_remember_upto_here_1 (y
);
4322 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4324 else if (SCM_REALP (y
))
4325 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4326 SCM_COMPLEX_IMAG (x
));
4327 else if (SCM_COMPLEXP (y
))
4328 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4329 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4330 else if (SCM_FRACTIONP (y
))
4331 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4332 SCM_COMPLEX_IMAG (x
));
4334 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4336 else if (SCM_FRACTIONP (x
))
4338 if (SCM_I_INUMP (y
))
4339 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4340 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4341 SCM_FRACTION_DENOMINATOR (x
));
4342 else if (SCM_BIGP (y
))
4343 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4344 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4345 SCM_FRACTION_DENOMINATOR (x
));
4346 else if (SCM_REALP (y
))
4347 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4348 else if (SCM_COMPLEXP (y
))
4349 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4350 SCM_COMPLEX_IMAG (y
));
4351 else if (SCM_FRACTIONP (y
))
4352 /* a/b + c/d = (ad + bc) / bd */
4353 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4354 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4355 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4357 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4360 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4364 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4366 "Return @math{@var{x}+1}.")
4367 #define FUNC_NAME s_scm_oneplus
4369 return scm_sum (x
, SCM_INUM1
);
4374 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
4375 (SCM x
, SCM y
, SCM rest
),
4376 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
4377 "the sum of all but the first argument are subtracted from the first\n"
4379 #define FUNC_NAME s_scm_i_difference
4381 while (!scm_is_null (rest
))
4382 { x
= scm_difference (x
, y
);
4384 rest
= scm_cdr (rest
);
4386 return scm_difference (x
, y
);
4390 #define s_difference s_scm_i_difference
4391 #define g_difference g_scm_i_difference
4394 scm_difference (SCM x
, SCM y
)
4395 #define FUNC_NAME s_difference
4397 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4400 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4402 if (SCM_I_INUMP (x
))
4404 scm_t_inum xx
= -SCM_I_INUM (x
);
4405 if (SCM_FIXABLE (xx
))
4406 return SCM_I_MAKINUM (xx
);
4408 return scm_i_inum2big (xx
);
4410 else if (SCM_BIGP (x
))
4411 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4412 bignum, but negating that gives a fixnum. */
4413 return scm_i_normbig (scm_i_clonebig (x
, 0));
4414 else if (SCM_REALP (x
))
4415 return scm_from_double (-SCM_REAL_VALUE (x
));
4416 else if (SCM_COMPLEXP (x
))
4417 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4418 -SCM_COMPLEX_IMAG (x
));
4419 else if (SCM_FRACTIONP (x
))
4420 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4421 SCM_FRACTION_DENOMINATOR (x
));
4423 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4426 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4428 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4430 scm_t_inum xx
= SCM_I_INUM (x
);
4431 scm_t_inum yy
= SCM_I_INUM (y
);
4432 scm_t_inum z
= xx
- yy
;
4433 if (SCM_FIXABLE (z
))
4434 return SCM_I_MAKINUM (z
);
4436 return scm_i_inum2big (z
);
4438 else if (SCM_BIGP (y
))
4440 /* inum-x - big-y */
4441 scm_t_inum xx
= SCM_I_INUM (x
);
4445 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4446 bignum, but negating that gives a fixnum. */
4447 return scm_i_normbig (scm_i_clonebig (y
, 0));
4451 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4452 SCM result
= scm_i_mkbig ();
4455 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4458 /* x - y == -(y + -x) */
4459 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4460 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4462 scm_remember_upto_here_1 (y
);
4464 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4465 /* we know the result will have to be a bignum */
4468 return scm_i_normbig (result
);
4471 else if (SCM_REALP (y
))
4473 scm_t_inum xx
= SCM_I_INUM (x
);
4474 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4476 else if (SCM_COMPLEXP (y
))
4478 scm_t_inum xx
= SCM_I_INUM (x
);
4479 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4480 - SCM_COMPLEX_IMAG (y
));
4482 else if (SCM_FRACTIONP (y
))
4483 /* a - b/c = (ac - b) / c */
4484 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4485 SCM_FRACTION_NUMERATOR (y
)),
4486 SCM_FRACTION_DENOMINATOR (y
));
4488 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4490 else if (SCM_BIGP (x
))
4492 if (SCM_I_INUMP (y
))
4494 /* big-x - inum-y */
4495 scm_t_inum yy
= SCM_I_INUM (y
);
4496 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4498 scm_remember_upto_here_1 (x
);
4500 return (SCM_FIXABLE (-yy
) ?
4501 SCM_I_MAKINUM (-yy
) : scm_from_inum (-yy
));
4504 SCM result
= scm_i_mkbig ();
4507 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4509 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4510 scm_remember_upto_here_1 (x
);
4512 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4513 /* we know the result will have to be a bignum */
4516 return scm_i_normbig (result
);
4519 else if (SCM_BIGP (y
))
4521 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4522 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4523 SCM result
= scm_i_mkbig ();
4524 mpz_sub (SCM_I_BIG_MPZ (result
),
4527 scm_remember_upto_here_2 (x
, y
);
4528 /* we know the result will have to be a bignum */
4529 if ((sgn_x
== 1) && (sgn_y
== -1))
4531 if ((sgn_x
== -1) && (sgn_y
== 1))
4533 return scm_i_normbig (result
);
4535 else if (SCM_REALP (y
))
4537 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4538 scm_remember_upto_here_1 (x
);
4539 return scm_from_double (result
);
4541 else if (SCM_COMPLEXP (y
))
4543 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4544 - SCM_COMPLEX_REAL (y
));
4545 scm_remember_upto_here_1 (x
);
4546 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4548 else if (SCM_FRACTIONP (y
))
4549 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4550 SCM_FRACTION_NUMERATOR (y
)),
4551 SCM_FRACTION_DENOMINATOR (y
));
4552 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4554 else if (SCM_REALP (x
))
4556 if (SCM_I_INUMP (y
))
4557 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4558 else if (SCM_BIGP (y
))
4560 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4561 scm_remember_upto_here_1 (x
);
4562 return scm_from_double (result
);
4564 else if (SCM_REALP (y
))
4565 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4566 else if (SCM_COMPLEXP (y
))
4567 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4568 -SCM_COMPLEX_IMAG (y
));
4569 else if (SCM_FRACTIONP (y
))
4570 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4572 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4574 else if (SCM_COMPLEXP (x
))
4576 if (SCM_I_INUMP (y
))
4577 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4578 SCM_COMPLEX_IMAG (x
));
4579 else if (SCM_BIGP (y
))
4581 double real_part
= (SCM_COMPLEX_REAL (x
)
4582 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4583 scm_remember_upto_here_1 (x
);
4584 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4586 else if (SCM_REALP (y
))
4587 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4588 SCM_COMPLEX_IMAG (x
));
4589 else if (SCM_COMPLEXP (y
))
4590 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4591 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4592 else if (SCM_FRACTIONP (y
))
4593 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4594 SCM_COMPLEX_IMAG (x
));
4596 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4598 else if (SCM_FRACTIONP (x
))
4600 if (SCM_I_INUMP (y
))
4601 /* a/b - c = (a - cb) / b */
4602 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4603 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4604 SCM_FRACTION_DENOMINATOR (x
));
4605 else if (SCM_BIGP (y
))
4606 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4607 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4608 SCM_FRACTION_DENOMINATOR (x
));
4609 else if (SCM_REALP (y
))
4610 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4611 else if (SCM_COMPLEXP (y
))
4612 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4613 -SCM_COMPLEX_IMAG (y
));
4614 else if (SCM_FRACTIONP (y
))
4615 /* a/b - c/d = (ad - bc) / bd */
4616 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4617 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4618 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4620 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4623 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4628 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4630 "Return @math{@var{x}-1}.")
4631 #define FUNC_NAME s_scm_oneminus
4633 return scm_difference (x
, SCM_INUM1
);
4638 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
4639 (SCM x
, SCM y
, SCM rest
),
4640 "Return the product of all arguments. If called without arguments,\n"
4642 #define FUNC_NAME s_scm_i_product
4644 while (!scm_is_null (rest
))
4645 { x
= scm_product (x
, y
);
4647 rest
= scm_cdr (rest
);
4649 return scm_product (x
, y
);
4653 #define s_product s_scm_i_product
4654 #define g_product g_scm_i_product
4657 scm_product (SCM x
, SCM y
)
4659 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4662 return SCM_I_MAKINUM (1L);
4663 else if (SCM_NUMBERP (x
))
4666 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4669 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4674 xx
= SCM_I_INUM (x
);
4678 case 0: return x
; break;
4679 case 1: return y
; break;
4681 * The following case (x = -1) is important for more than
4682 * just optimization. It handles the case of negating
4683 * (+ 1 most-positive-fixnum) aka (- most-negative-fixnum),
4684 * which is a bignum that must be changed back into a fixnum.
4685 * Failure to do so will cause the following to return #f:
4686 * (= most-negative-fixnum (* -1 (- most-negative-fixnum)))
4689 return scm_difference(y
, SCM_UNDEFINED
);
4693 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4695 scm_t_inum yy
= SCM_I_INUM (y
);
4696 scm_t_inum kk
= xx
* yy
;
4697 SCM k
= SCM_I_MAKINUM (kk
);
4698 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4702 SCM result
= scm_i_inum2big (xx
);
4703 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4704 return scm_i_normbig (result
);
4707 else if (SCM_BIGP (y
))
4709 SCM result
= scm_i_mkbig ();
4710 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4711 scm_remember_upto_here_1 (y
);
4714 else if (SCM_REALP (y
))
4715 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4716 else if (SCM_COMPLEXP (y
))
4717 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4718 xx
* SCM_COMPLEX_IMAG (y
));
4719 else if (SCM_FRACTIONP (y
))
4720 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4721 SCM_FRACTION_DENOMINATOR (y
));
4723 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4725 else if (SCM_BIGP (x
))
4727 if (SCM_I_INUMP (y
))
4732 else if (SCM_BIGP (y
))
4734 SCM result
= scm_i_mkbig ();
4735 mpz_mul (SCM_I_BIG_MPZ (result
),
4738 scm_remember_upto_here_2 (x
, y
);
4741 else if (SCM_REALP (y
))
4743 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4744 scm_remember_upto_here_1 (x
);
4745 return scm_from_double (result
);
4747 else if (SCM_COMPLEXP (y
))
4749 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4750 scm_remember_upto_here_1 (x
);
4751 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4752 z
* SCM_COMPLEX_IMAG (y
));
4754 else if (SCM_FRACTIONP (y
))
4755 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4756 SCM_FRACTION_DENOMINATOR (y
));
4758 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4760 else if (SCM_REALP (x
))
4762 if (SCM_I_INUMP (y
))
4764 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4765 if (scm_is_eq (y
, SCM_INUM0
))
4767 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4769 else if (SCM_BIGP (y
))
4771 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4772 scm_remember_upto_here_1 (y
);
4773 return scm_from_double (result
);
4775 else if (SCM_REALP (y
))
4776 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4777 else if (SCM_COMPLEXP (y
))
4778 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4779 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4780 else if (SCM_FRACTIONP (y
))
4781 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4783 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4785 else if (SCM_COMPLEXP (x
))
4787 if (SCM_I_INUMP (y
))
4789 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4790 if (scm_is_eq (y
, SCM_INUM0
))
4792 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4793 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4795 else if (SCM_BIGP (y
))
4797 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4798 scm_remember_upto_here_1 (y
);
4799 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4800 z
* SCM_COMPLEX_IMAG (x
));
4802 else if (SCM_REALP (y
))
4803 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4804 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4805 else if (SCM_COMPLEXP (y
))
4807 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4808 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4809 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4810 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4812 else if (SCM_FRACTIONP (y
))
4814 double yy
= scm_i_fraction2double (y
);
4815 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4816 yy
* SCM_COMPLEX_IMAG (x
));
4819 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4821 else if (SCM_FRACTIONP (x
))
4823 if (SCM_I_INUMP (y
))
4824 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4825 SCM_FRACTION_DENOMINATOR (x
));
4826 else if (SCM_BIGP (y
))
4827 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4828 SCM_FRACTION_DENOMINATOR (x
));
4829 else if (SCM_REALP (y
))
4830 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4831 else if (SCM_COMPLEXP (y
))
4833 double xx
= scm_i_fraction2double (x
);
4834 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4835 xx
* SCM_COMPLEX_IMAG (y
));
4837 else if (SCM_FRACTIONP (y
))
4838 /* a/b * c/d = ac / bd */
4839 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4840 SCM_FRACTION_NUMERATOR (y
)),
4841 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4842 SCM_FRACTION_DENOMINATOR (y
)));
4844 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4847 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4850 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4851 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4852 #define ALLOW_DIVIDE_BY_ZERO
4853 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4856 /* The code below for complex division is adapted from the GNU
4857 libstdc++, which adapted it from f2c's libF77, and is subject to
4860 /****************************************************************
4861 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4863 Permission to use, copy, modify, and distribute this software
4864 and its documentation for any purpose and without fee is hereby
4865 granted, provided that the above copyright notice appear in all
4866 copies and that both that the copyright notice and this
4867 permission notice and warranty disclaimer appear in supporting
4868 documentation, and that the names of AT&T Bell Laboratories or
4869 Bellcore or any of their entities not be used in advertising or
4870 publicity pertaining to distribution of the software without
4871 specific, written prior permission.
4873 AT&T and Bellcore disclaim all warranties with regard to this
4874 software, including all implied warranties of merchantability
4875 and fitness. In no event shall AT&T or Bellcore be liable for
4876 any special, indirect or consequential damages or any damages
4877 whatsoever resulting from loss of use, data or profits, whether
4878 in an action of contract, negligence or other tortious action,
4879 arising out of or in connection with the use or performance of
4881 ****************************************************************/
4883 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
4884 (SCM x
, SCM y
, SCM rest
),
4885 "Divide the first argument by the product of the remaining\n"
4886 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
4888 #define FUNC_NAME s_scm_i_divide
4890 while (!scm_is_null (rest
))
4891 { x
= scm_divide (x
, y
);
4893 rest
= scm_cdr (rest
);
4895 return scm_divide (x
, y
);
4899 #define s_divide s_scm_i_divide
4900 #define g_divide g_scm_i_divide
4903 do_divide (SCM x
, SCM y
, int inexact
)
4904 #define FUNC_NAME s_divide
4908 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4911 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4912 else if (SCM_I_INUMP (x
))
4914 scm_t_inum xx
= SCM_I_INUM (x
);
4915 if (xx
== 1 || xx
== -1)
4917 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4919 scm_num_overflow (s_divide
);
4924 return scm_from_double (1.0 / (double) xx
);
4925 else return scm_i_make_ratio (SCM_INUM1
, x
);
4928 else if (SCM_BIGP (x
))
4931 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4932 else return scm_i_make_ratio (SCM_INUM1
, x
);
4934 else if (SCM_REALP (x
))
4936 double xx
= SCM_REAL_VALUE (x
);
4937 #ifndef ALLOW_DIVIDE_BY_ZERO
4939 scm_num_overflow (s_divide
);
4942 return scm_from_double (1.0 / xx
);
4944 else if (SCM_COMPLEXP (x
))
4946 double r
= SCM_COMPLEX_REAL (x
);
4947 double i
= SCM_COMPLEX_IMAG (x
);
4948 if (fabs(r
) <= fabs(i
))
4951 double d
= i
* (1.0 + t
* t
);
4952 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4957 double d
= r
* (1.0 + t
* t
);
4958 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4961 else if (SCM_FRACTIONP (x
))
4962 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4963 SCM_FRACTION_NUMERATOR (x
));
4965 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4968 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4970 scm_t_inum xx
= SCM_I_INUM (x
);
4971 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4973 scm_t_inum yy
= SCM_I_INUM (y
);
4976 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4977 scm_num_overflow (s_divide
);
4979 return scm_from_double ((double) xx
/ (double) yy
);
4982 else if (xx
% yy
!= 0)
4985 return scm_from_double ((double) xx
/ (double) yy
);
4986 else return scm_i_make_ratio (x
, y
);
4990 scm_t_inum z
= xx
/ yy
;
4991 if (SCM_FIXABLE (z
))
4992 return SCM_I_MAKINUM (z
);
4994 return scm_i_inum2big (z
);
4997 else if (SCM_BIGP (y
))
5000 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
5001 else return scm_i_make_ratio (x
, y
);
5003 else if (SCM_REALP (y
))
5005 double yy
= SCM_REAL_VALUE (y
);
5006 #ifndef ALLOW_DIVIDE_BY_ZERO
5008 scm_num_overflow (s_divide
);
5011 return scm_from_double ((double) xx
/ yy
);
5013 else if (SCM_COMPLEXP (y
))
5016 complex_div
: /* y _must_ be a complex number */
5018 double r
= SCM_COMPLEX_REAL (y
);
5019 double i
= SCM_COMPLEX_IMAG (y
);
5020 if (fabs(r
) <= fabs(i
))
5023 double d
= i
* (1.0 + t
* t
);
5024 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
5029 double d
= r
* (1.0 + t
* t
);
5030 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
5034 else if (SCM_FRACTIONP (y
))
5035 /* a / b/c = ac / b */
5036 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5037 SCM_FRACTION_NUMERATOR (y
));
5039 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5041 else if (SCM_BIGP (x
))
5043 if (SCM_I_INUMP (y
))
5045 scm_t_inum yy
= SCM_I_INUM (y
);
5048 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5049 scm_num_overflow (s_divide
);
5051 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5052 scm_remember_upto_here_1 (x
);
5053 return (sgn
== 0) ? scm_nan () : scm_inf ();
5060 /* FIXME: HMM, what are the relative performance issues here?
5061 We need to test. Is it faster on average to test
5062 divisible_p, then perform whichever operation, or is it
5063 faster to perform the integer div opportunistically and
5064 switch to real if there's a remainder? For now we take the
5065 middle ground: test, then if divisible, use the faster div
5068 scm_t_inum abs_yy
= yy
< 0 ? -yy
: yy
;
5069 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
5073 SCM result
= scm_i_mkbig ();
5074 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
5075 scm_remember_upto_here_1 (x
);
5077 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5078 return scm_i_normbig (result
);
5083 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
5084 else return scm_i_make_ratio (x
, y
);
5088 else if (SCM_BIGP (y
))
5093 /* It's easily possible for the ratio x/y to fit a double
5094 but one or both x and y be too big to fit a double,
5095 hence the use of mpq_get_d rather than converting and
5098 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
5099 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
5100 return scm_from_double (mpq_get_d (q
));
5104 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
5108 SCM result
= scm_i_mkbig ();
5109 mpz_divexact (SCM_I_BIG_MPZ (result
),
5112 scm_remember_upto_here_2 (x
, y
);
5113 return scm_i_normbig (result
);
5116 return scm_i_make_ratio (x
, y
);
5119 else if (SCM_REALP (y
))
5121 double yy
= SCM_REAL_VALUE (y
);
5122 #ifndef ALLOW_DIVIDE_BY_ZERO
5124 scm_num_overflow (s_divide
);
5127 return scm_from_double (scm_i_big2dbl (x
) / yy
);
5129 else if (SCM_COMPLEXP (y
))
5131 a
= scm_i_big2dbl (x
);
5134 else if (SCM_FRACTIONP (y
))
5135 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5136 SCM_FRACTION_NUMERATOR (y
));
5138 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5140 else if (SCM_REALP (x
))
5142 double rx
= SCM_REAL_VALUE (x
);
5143 if (SCM_I_INUMP (y
))
5145 scm_t_inum yy
= SCM_I_INUM (y
);
5146 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5148 scm_num_overflow (s_divide
);
5151 return scm_from_double (rx
/ (double) yy
);
5153 else if (SCM_BIGP (y
))
5155 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5156 scm_remember_upto_here_1 (y
);
5157 return scm_from_double (rx
/ dby
);
5159 else if (SCM_REALP (y
))
5161 double yy
= SCM_REAL_VALUE (y
);
5162 #ifndef ALLOW_DIVIDE_BY_ZERO
5164 scm_num_overflow (s_divide
);
5167 return scm_from_double (rx
/ yy
);
5169 else if (SCM_COMPLEXP (y
))
5174 else if (SCM_FRACTIONP (y
))
5175 return scm_from_double (rx
/ scm_i_fraction2double (y
));
5177 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5179 else if (SCM_COMPLEXP (x
))
5181 double rx
= SCM_COMPLEX_REAL (x
);
5182 double ix
= SCM_COMPLEX_IMAG (x
);
5183 if (SCM_I_INUMP (y
))
5185 scm_t_inum yy
= SCM_I_INUM (y
);
5186 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5188 scm_num_overflow (s_divide
);
5193 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
5196 else if (SCM_BIGP (y
))
5198 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5199 scm_remember_upto_here_1 (y
);
5200 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
5202 else if (SCM_REALP (y
))
5204 double yy
= SCM_REAL_VALUE (y
);
5205 #ifndef ALLOW_DIVIDE_BY_ZERO
5207 scm_num_overflow (s_divide
);
5210 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5212 else if (SCM_COMPLEXP (y
))
5214 double ry
= SCM_COMPLEX_REAL (y
);
5215 double iy
= SCM_COMPLEX_IMAG (y
);
5216 if (fabs(ry
) <= fabs(iy
))
5219 double d
= iy
* (1.0 + t
* t
);
5220 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5225 double d
= ry
* (1.0 + t
* t
);
5226 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5229 else if (SCM_FRACTIONP (y
))
5231 double yy
= scm_i_fraction2double (y
);
5232 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5235 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5237 else if (SCM_FRACTIONP (x
))
5239 if (SCM_I_INUMP (y
))
5241 scm_t_inum yy
= SCM_I_INUM (y
);
5242 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5244 scm_num_overflow (s_divide
);
5247 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5248 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5250 else if (SCM_BIGP (y
))
5252 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5253 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5255 else if (SCM_REALP (y
))
5257 double yy
= SCM_REAL_VALUE (y
);
5258 #ifndef ALLOW_DIVIDE_BY_ZERO
5260 scm_num_overflow (s_divide
);
5263 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5265 else if (SCM_COMPLEXP (y
))
5267 a
= scm_i_fraction2double (x
);
5270 else if (SCM_FRACTIONP (y
))
5271 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5272 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5274 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5277 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5281 scm_divide (SCM x
, SCM y
)
5283 return do_divide (x
, y
, 0);
5286 static SCM
scm_divide2real (SCM x
, SCM y
)
5288 return do_divide (x
, y
, 1);
5294 scm_c_truncate (double x
)
5305 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5306 half-way case (ie. when x is an integer plus 0.5) going upwards.
5307 Then half-way cases are identified and adjusted down if the
5308 round-upwards didn't give the desired even integer.
5310 "plus_half == result" identifies a half-way case. If plus_half, which is
5311 x + 0.5, is an integer then x must be an integer plus 0.5.
5313 An odd "result" value is identified with result/2 != floor(result/2).
5314 This is done with plus_half, since that value is ready for use sooner in
5315 a pipelined cpu, and we're already requiring plus_half == result.
5317 Note however that we need to be careful when x is big and already an
5318 integer. In that case "x+0.5" may round to an adjacent integer, causing
5319 us to return such a value, incorrectly. For instance if the hardware is
5320 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5321 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5322 returned. Or if the hardware is in round-upwards mode, then other bigger
5323 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5324 representable value, 2^128+2^76 (or whatever), again incorrect.
5326 These bad roundings of x+0.5 are avoided by testing at the start whether
5327 x is already an integer. If it is then clearly that's the desired result
5328 already. And if it's not then the exponent must be small enough to allow
5329 an 0.5 to be represented, and hence added without a bad rounding. */
5332 scm_c_round (double x
)
5334 double plus_half
, result
;
5339 plus_half
= x
+ 0.5;
5340 result
= floor (plus_half
);
5341 /* Adjust so that the rounding is towards even. */
5342 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5347 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5349 "Round the number @var{x} towards zero.")
5350 #define FUNC_NAME s_scm_truncate_number
5352 if (scm_is_false (scm_negative_p (x
)))
5353 return scm_floor (x
);
5355 return scm_ceiling (x
);
5359 static SCM exactly_one_half
;
5361 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5363 "Round the number @var{x} towards the nearest integer. "
5364 "When it is exactly halfway between two integers, "
5365 "round towards the even one.")
5366 #define FUNC_NAME s_scm_round_number
5368 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5370 else if (SCM_REALP (x
))
5371 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5374 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5375 single quotient+remainder division then examining to see which way
5376 the rounding should go. */
5377 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5378 SCM result
= scm_floor (plus_half
);
5379 /* Adjust so that the rounding is towards even. */
5380 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5381 && scm_is_true (scm_odd_p (result
)))
5382 return scm_difference (result
, SCM_INUM1
);
5389 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5391 "Round the number @var{x} towards minus infinity.")
5392 #define FUNC_NAME s_scm_floor
5394 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5396 else if (SCM_REALP (x
))
5397 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5398 else if (SCM_FRACTIONP (x
))
5400 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5401 SCM_FRACTION_DENOMINATOR (x
));
5402 if (scm_is_false (scm_negative_p (x
)))
5404 /* For positive x, rounding towards zero is correct. */
5409 /* For negative x, we need to return q-1 unless x is an
5410 integer. But fractions are never integer, per our
5412 return scm_difference (q
, SCM_INUM1
);
5416 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5420 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5422 "Round the number @var{x} towards infinity.")
5423 #define FUNC_NAME s_scm_ceiling
5425 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5427 else if (SCM_REALP (x
))
5428 return scm_from_double (ceil (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_positive_p (x
)))
5435 /* For negative x, rounding towards zero is correct. */
5440 /* For positive x, we need to return q+1 unless x is an
5441 integer. But fractions are never integer, per our
5443 return scm_sum (q
, SCM_INUM1
);
5447 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5451 /* sin/cos/tan/asin/acos/atan
5452 sinh/cosh/tanh/asinh/acosh/atanh
5453 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5454 Written by Jerry D. Hedden, (C) FSF.
5455 See the file `COPYING' for terms applying to this program. */
5457 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5459 "Return @var{x} raised to the power of @var{y}.")
5460 #define FUNC_NAME s_scm_expt
5462 if (scm_is_integer (y
))
5464 if (scm_is_true (scm_exact_p (y
)))
5465 return scm_integer_expt (x
, y
);
5468 /* Here we handle the case where the exponent is an inexact
5469 integer. We make the exponent exact in order to use
5470 scm_integer_expt, and thus avoid the spurious imaginary
5471 parts that may result from round-off errors in the general
5472 e^(y log x) method below (for example when squaring a large
5473 negative number). In this case, we must return an inexact
5474 result for correctness. We also make the base inexact so
5475 that scm_integer_expt will use fast inexact arithmetic
5476 internally. Note that making the base inexact is not
5477 sufficient to guarantee an inexact result, because
5478 scm_integer_expt will return an exact 1 when the exponent
5479 is 0, even if the base is inexact. */
5480 return scm_exact_to_inexact
5481 (scm_integer_expt (scm_exact_to_inexact (x
),
5482 scm_inexact_to_exact (y
)));
5485 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5487 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5490 return scm_exp (scm_product (scm_log (x
), y
));
5494 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5496 "Compute the sine of @var{z}.")
5497 #define FUNC_NAME s_scm_sin
5499 if (scm_is_real (z
))
5500 return scm_from_double (sin (scm_to_double (z
)));
5501 else if (SCM_COMPLEXP (z
))
5503 x
= SCM_COMPLEX_REAL (z
);
5504 y
= SCM_COMPLEX_IMAG (z
);
5505 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5506 cos (x
) * sinh (y
));
5509 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5513 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5515 "Compute the cosine of @var{z}.")
5516 #define FUNC_NAME s_scm_cos
5518 if (scm_is_real (z
))
5519 return scm_from_double (cos (scm_to_double (z
)));
5520 else if (SCM_COMPLEXP (z
))
5522 x
= SCM_COMPLEX_REAL (z
);
5523 y
= SCM_COMPLEX_IMAG (z
);
5524 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5525 -sin (x
) * sinh (y
));
5528 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5532 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5534 "Compute the tangent of @var{z}.")
5535 #define FUNC_NAME s_scm_tan
5537 if (scm_is_real (z
))
5538 return scm_from_double (tan (scm_to_double (z
)));
5539 else if (SCM_COMPLEXP (z
))
5541 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5542 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5543 w
= cos (x
) + cosh (y
);
5544 #ifndef ALLOW_DIVIDE_BY_ZERO
5546 scm_num_overflow (s_scm_tan
);
5548 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5551 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5555 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5557 "Compute the hyperbolic sine of @var{z}.")
5558 #define FUNC_NAME s_scm_sinh
5560 if (scm_is_real (z
))
5561 return scm_from_double (sinh (scm_to_double (z
)));
5562 else if (SCM_COMPLEXP (z
))
5564 x
= SCM_COMPLEX_REAL (z
);
5565 y
= SCM_COMPLEX_IMAG (z
);
5566 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5567 cosh (x
) * sin (y
));
5570 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5574 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5576 "Compute the hyperbolic cosine of @var{z}.")
5577 #define FUNC_NAME s_scm_cosh
5579 if (scm_is_real (z
))
5580 return scm_from_double (cosh (scm_to_double (z
)));
5581 else if (SCM_COMPLEXP (z
))
5583 x
= SCM_COMPLEX_REAL (z
);
5584 y
= SCM_COMPLEX_IMAG (z
);
5585 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5586 sinh (x
) * sin (y
));
5589 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5593 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5595 "Compute the hyperbolic tangent of @var{z}.")
5596 #define FUNC_NAME s_scm_tanh
5598 if (scm_is_real (z
))
5599 return scm_from_double (tanh (scm_to_double (z
)));
5600 else if (SCM_COMPLEXP (z
))
5602 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5603 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5604 w
= cosh (x
) + cos (y
);
5605 #ifndef ALLOW_DIVIDE_BY_ZERO
5607 scm_num_overflow (s_scm_tanh
);
5609 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5612 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5616 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5618 "Compute the arc sine of @var{z}.")
5619 #define FUNC_NAME s_scm_asin
5621 if (scm_is_real (z
))
5623 double w
= scm_to_double (z
);
5624 if (w
>= -1.0 && w
<= 1.0)
5625 return scm_from_double (asin (w
));
5627 return scm_product (scm_c_make_rectangular (0, -1),
5628 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5630 else if (SCM_COMPLEXP (z
))
5632 x
= SCM_COMPLEX_REAL (z
);
5633 y
= SCM_COMPLEX_IMAG (z
);
5634 return scm_product (scm_c_make_rectangular (0, -1),
5635 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5638 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5642 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5644 "Compute the arc cosine of @var{z}.")
5645 #define FUNC_NAME s_scm_acos
5647 if (scm_is_real (z
))
5649 double w
= scm_to_double (z
);
5650 if (w
>= -1.0 && w
<= 1.0)
5651 return scm_from_double (acos (w
));
5653 return scm_sum (scm_from_double (acos (0.0)),
5654 scm_product (scm_c_make_rectangular (0, 1),
5655 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5657 else if (SCM_COMPLEXP (z
))
5659 x
= SCM_COMPLEX_REAL (z
);
5660 y
= SCM_COMPLEX_IMAG (z
);
5661 return scm_sum (scm_from_double (acos (0.0)),
5662 scm_product (scm_c_make_rectangular (0, 1),
5663 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5666 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5670 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5672 "With one argument, compute the arc tangent of @var{z}.\n"
5673 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5674 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5675 #define FUNC_NAME s_scm_atan
5679 if (scm_is_real (z
))
5680 return scm_from_double (atan (scm_to_double (z
)));
5681 else if (SCM_COMPLEXP (z
))
5684 v
= SCM_COMPLEX_REAL (z
);
5685 w
= SCM_COMPLEX_IMAG (z
);
5686 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5687 scm_c_make_rectangular (v
, w
+ 1.0))),
5688 scm_c_make_rectangular (0, 2));
5691 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5693 else if (scm_is_real (z
))
5695 if (scm_is_real (y
))
5696 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5698 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5701 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5705 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5707 "Compute the inverse hyperbolic sine of @var{z}.")
5708 #define FUNC_NAME s_scm_sys_asinh
5710 if (scm_is_real (z
))
5711 return scm_from_double (asinh (scm_to_double (z
)));
5712 else if (scm_is_number (z
))
5713 return scm_log (scm_sum (z
,
5714 scm_sqrt (scm_sum (scm_product (z
, z
),
5717 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5721 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5723 "Compute the inverse hyperbolic cosine of @var{z}.")
5724 #define FUNC_NAME s_scm_sys_acosh
5726 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5727 return scm_from_double (acosh (scm_to_double (z
)));
5728 else if (scm_is_number (z
))
5729 return scm_log (scm_sum (z
,
5730 scm_sqrt (scm_difference (scm_product (z
, z
),
5733 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5737 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5739 "Compute the inverse hyperbolic tangent of @var{z}.")
5740 #define FUNC_NAME s_scm_sys_atanh
5742 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5743 return scm_from_double (atanh (scm_to_double (z
)));
5744 else if (scm_is_number (z
))
5745 return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1
, z
),
5746 scm_difference (SCM_INUM1
, z
))),
5749 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5754 scm_c_make_rectangular (double re
, double im
)
5757 return scm_from_double (re
);
5762 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5764 SCM_SET_CELL_TYPE (z
, scm_tc16_complex
);
5765 SCM_COMPLEX_REAL (z
) = re
;
5766 SCM_COMPLEX_IMAG (z
) = im
;
5771 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5772 (SCM real_part
, SCM imaginary_part
),
5773 "Return a complex number constructed of the given @var{real-part} "
5774 "and @var{imaginary-part} parts.")
5775 #define FUNC_NAME s_scm_make_rectangular
5777 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5778 SCM_ARG1
, FUNC_NAME
, "real");
5779 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5780 SCM_ARG2
, FUNC_NAME
, "real");
5781 return scm_c_make_rectangular (scm_to_double (real_part
),
5782 scm_to_double (imaginary_part
));
5787 scm_c_make_polar (double mag
, double ang
)
5791 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5792 use it on Glibc-based systems that have it (it's a GNU extension). See
5793 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5795 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5796 sincos (ang
, &s
, &c
);
5801 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5804 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5806 "Return the complex number @var{x} * e^(i * @var{y}).")
5807 #define FUNC_NAME s_scm_make_polar
5809 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5810 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5811 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5816 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5817 /* "Return the real part of the number @var{z}."
5820 scm_real_part (SCM z
)
5822 if (SCM_I_INUMP (z
))
5824 else if (SCM_BIGP (z
))
5826 else if (SCM_REALP (z
))
5828 else if (SCM_COMPLEXP (z
))
5829 return scm_from_double (SCM_COMPLEX_REAL (z
));
5830 else if (SCM_FRACTIONP (z
))
5833 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5837 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5838 /* "Return the imaginary part of the number @var{z}."
5841 scm_imag_part (SCM z
)
5843 if (SCM_I_INUMP (z
))
5845 else if (SCM_BIGP (z
))
5847 else if (SCM_REALP (z
))
5849 else if (SCM_COMPLEXP (z
))
5850 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5851 else if (SCM_FRACTIONP (z
))
5854 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5857 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5858 /* "Return the numerator of the number @var{z}."
5861 scm_numerator (SCM z
)
5863 if (SCM_I_INUMP (z
))
5865 else if (SCM_BIGP (z
))
5867 else if (SCM_FRACTIONP (z
))
5868 return SCM_FRACTION_NUMERATOR (z
);
5869 else if (SCM_REALP (z
))
5870 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5872 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5876 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5877 /* "Return the denominator of the number @var{z}."
5880 scm_denominator (SCM z
)
5882 if (SCM_I_INUMP (z
))
5884 else if (SCM_BIGP (z
))
5886 else if (SCM_FRACTIONP (z
))
5887 return SCM_FRACTION_DENOMINATOR (z
);
5888 else if (SCM_REALP (z
))
5889 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5891 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5894 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5895 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5896 * "@code{abs} for real arguments, but also allows complex numbers."
5899 scm_magnitude (SCM z
)
5901 if (SCM_I_INUMP (z
))
5903 scm_t_inum zz
= SCM_I_INUM (z
);
5906 else if (SCM_POSFIXABLE (-zz
))
5907 return SCM_I_MAKINUM (-zz
);
5909 return scm_i_inum2big (-zz
);
5911 else if (SCM_BIGP (z
))
5913 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5914 scm_remember_upto_here_1 (z
);
5916 return scm_i_clonebig (z
, 0);
5920 else if (SCM_REALP (z
))
5921 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5922 else if (SCM_COMPLEXP (z
))
5923 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5924 else if (SCM_FRACTIONP (z
))
5926 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5928 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5929 SCM_FRACTION_DENOMINATOR (z
));
5932 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5936 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5937 /* "Return the angle of the complex number @var{z}."
5942 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5943 flo0 to save allocating a new flonum with scm_from_double each time.
5944 But if atan2 follows the floating point rounding mode, then the value
5945 is not a constant. Maybe it'd be close enough though. */
5946 if (SCM_I_INUMP (z
))
5948 if (SCM_I_INUM (z
) >= 0)
5951 return scm_from_double (atan2 (0.0, -1.0));
5953 else if (SCM_BIGP (z
))
5955 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5956 scm_remember_upto_here_1 (z
);
5958 return scm_from_double (atan2 (0.0, -1.0));
5962 else if (SCM_REALP (z
))
5964 if (SCM_REAL_VALUE (z
) >= 0)
5967 return scm_from_double (atan2 (0.0, -1.0));
5969 else if (SCM_COMPLEXP (z
))
5970 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5971 else if (SCM_FRACTIONP (z
))
5973 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5975 else return scm_from_double (atan2 (0.0, -1.0));
5978 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5982 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5983 /* Convert the number @var{x} to its inexact representation.\n"
5986 scm_exact_to_inexact (SCM z
)
5988 if (SCM_I_INUMP (z
))
5989 return scm_from_double ((double) SCM_I_INUM (z
));
5990 else if (SCM_BIGP (z
))
5991 return scm_from_double (scm_i_big2dbl (z
));
5992 else if (SCM_FRACTIONP (z
))
5993 return scm_from_double (scm_i_fraction2double (z
));
5994 else if (SCM_INEXACTP (z
))
5997 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
6001 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
6003 "Return an exact number that is numerically closest to @var{z}.")
6004 #define FUNC_NAME s_scm_inexact_to_exact
6006 if (SCM_I_INUMP (z
))
6008 else if (SCM_BIGP (z
))
6010 else if (SCM_REALP (z
))
6012 if (isinf (SCM_REAL_VALUE (z
)) || isnan (SCM_REAL_VALUE (z
)))
6013 SCM_OUT_OF_RANGE (1, z
);
6020 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
6021 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
6022 scm_i_mpz2num (mpq_denref (frac
)));
6024 /* When scm_i_make_ratio throws, we leak the memory allocated
6031 else if (SCM_FRACTIONP (z
))
6034 SCM_WRONG_TYPE_ARG (1, z
);
6038 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
6040 "Returns the @emph{simplest} rational number differing\n"
6041 "from @var{x} by no more than @var{eps}.\n"
6043 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
6044 "exact result when both its arguments are exact. Thus, you might need\n"
6045 "to use @code{inexact->exact} on the arguments.\n"
6048 "(rationalize (inexact->exact 1.2) 1/100)\n"
6051 #define FUNC_NAME s_scm_rationalize
6053 if (SCM_I_INUMP (x
))
6055 else if (SCM_BIGP (x
))
6057 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
6059 /* Use continued fractions to find closest ratio. All
6060 arithmetic is done with exact numbers.
6063 SCM ex
= scm_inexact_to_exact (x
);
6064 SCM int_part
= scm_floor (ex
);
6066 SCM a1
= SCM_INUM0
, a2
= SCM_INUM1
, a
= SCM_INUM0
;
6067 SCM b1
= SCM_INUM1
, b2
= SCM_INUM0
, b
= SCM_INUM0
;
6071 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
6074 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
6075 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
6077 /* We stop after a million iterations just to be absolutely sure
6078 that we don't go into an infinite loop. The process normally
6079 converges after less than a dozen iterations.
6082 eps
= scm_abs (eps
);
6083 while (++i
< 1000000)
6085 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
6086 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
6087 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
6089 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
6090 eps
))) /* abs(x-a/b) <= eps */
6092 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
6093 if (scm_is_false (scm_exact_p (x
))
6094 || scm_is_false (scm_exact_p (eps
)))
6095 return scm_exact_to_inexact (res
);
6099 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
6101 tt
= scm_floor (rx
); /* tt = floor (rx) */
6107 scm_num_overflow (s_scm_rationalize
);
6110 SCM_WRONG_TYPE_ARG (1, x
);
6114 /* conversion functions */
6117 scm_is_integer (SCM val
)
6119 return scm_is_true (scm_integer_p (val
));
6123 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
6125 if (SCM_I_INUMP (val
))
6127 scm_t_signed_bits n
= SCM_I_INUM (val
);
6128 return n
>= min
&& n
<= max
;
6130 else if (SCM_BIGP (val
))
6132 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
6134 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
6136 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
6138 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
6139 return n
>= min
&& n
<= max
;
6149 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6150 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6153 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6154 SCM_I_BIG_MPZ (val
));
6156 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
6168 return n
>= min
&& n
<= max
;
6176 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
6178 if (SCM_I_INUMP (val
))
6180 scm_t_signed_bits n
= SCM_I_INUM (val
);
6181 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
6183 else if (SCM_BIGP (val
))
6185 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
6187 else if (max
<= ULONG_MAX
)
6189 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
6191 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
6192 return n
>= min
&& n
<= max
;
6202 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
6205 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6206 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6209 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6210 SCM_I_BIG_MPZ (val
));
6212 return n
>= min
&& n
<= max
;
6220 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
6222 scm_error (scm_out_of_range_key
,
6224 "Value out of range ~S to ~S: ~S",
6225 scm_list_3 (min
, max
, bad_val
),
6226 scm_list_1 (bad_val
));
6229 #define TYPE scm_t_intmax
6230 #define TYPE_MIN min
6231 #define TYPE_MAX max
6232 #define SIZEOF_TYPE 0
6233 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
6234 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
6235 #include "libguile/conv-integer.i.c"
6237 #define TYPE scm_t_uintmax
6238 #define TYPE_MIN min
6239 #define TYPE_MAX max
6240 #define SIZEOF_TYPE 0
6241 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6242 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6243 #include "libguile/conv-uinteger.i.c"
6245 #define TYPE scm_t_int8
6246 #define TYPE_MIN SCM_T_INT8_MIN
6247 #define TYPE_MAX SCM_T_INT8_MAX
6248 #define SIZEOF_TYPE 1
6249 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6250 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6251 #include "libguile/conv-integer.i.c"
6253 #define TYPE scm_t_uint8
6255 #define TYPE_MAX SCM_T_UINT8_MAX
6256 #define SIZEOF_TYPE 1
6257 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6258 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6259 #include "libguile/conv-uinteger.i.c"
6261 #define TYPE scm_t_int16
6262 #define TYPE_MIN SCM_T_INT16_MIN
6263 #define TYPE_MAX SCM_T_INT16_MAX
6264 #define SIZEOF_TYPE 2
6265 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6266 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6267 #include "libguile/conv-integer.i.c"
6269 #define TYPE scm_t_uint16
6271 #define TYPE_MAX SCM_T_UINT16_MAX
6272 #define SIZEOF_TYPE 2
6273 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6274 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6275 #include "libguile/conv-uinteger.i.c"
6277 #define TYPE scm_t_int32
6278 #define TYPE_MIN SCM_T_INT32_MIN
6279 #define TYPE_MAX SCM_T_INT32_MAX
6280 #define SIZEOF_TYPE 4
6281 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6282 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6283 #include "libguile/conv-integer.i.c"
6285 #define TYPE scm_t_uint32
6287 #define TYPE_MAX SCM_T_UINT32_MAX
6288 #define SIZEOF_TYPE 4
6289 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6290 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6291 #include "libguile/conv-uinteger.i.c"
6293 #define TYPE scm_t_wchar
6294 #define TYPE_MIN (scm_t_int32)-1
6295 #define TYPE_MAX (scm_t_int32)0x10ffff
6296 #define SIZEOF_TYPE 4
6297 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6298 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6299 #include "libguile/conv-integer.i.c"
6301 #define TYPE scm_t_int64
6302 #define TYPE_MIN SCM_T_INT64_MIN
6303 #define TYPE_MAX SCM_T_INT64_MAX
6304 #define SIZEOF_TYPE 8
6305 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6306 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6307 #include "libguile/conv-integer.i.c"
6309 #define TYPE scm_t_uint64
6311 #define TYPE_MAX SCM_T_UINT64_MAX
6312 #define SIZEOF_TYPE 8
6313 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6314 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6315 #include "libguile/conv-uinteger.i.c"
6318 scm_to_mpz (SCM val
, mpz_t rop
)
6320 if (SCM_I_INUMP (val
))
6321 mpz_set_si (rop
, SCM_I_INUM (val
));
6322 else if (SCM_BIGP (val
))
6323 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6325 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6329 scm_from_mpz (mpz_t val
)
6331 return scm_i_mpz2num (val
);
6335 scm_is_real (SCM val
)
6337 return scm_is_true (scm_real_p (val
));
6341 scm_is_rational (SCM val
)
6343 return scm_is_true (scm_rational_p (val
));
6347 scm_to_double (SCM val
)
6349 if (SCM_I_INUMP (val
))
6350 return SCM_I_INUM (val
);
6351 else if (SCM_BIGP (val
))
6352 return scm_i_big2dbl (val
);
6353 else if (SCM_FRACTIONP (val
))
6354 return scm_i_fraction2double (val
);
6355 else if (SCM_REALP (val
))
6356 return SCM_REAL_VALUE (val
);
6358 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6362 scm_from_double (double val
)
6366 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
6368 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
6369 SCM_REAL_VALUE (z
) = val
;
6374 #if SCM_ENABLE_DEPRECATED == 1
6377 scm_num2float (SCM num
, unsigned long pos
, const char *s_caller
)
6379 scm_c_issue_deprecation_warning
6380 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
6384 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6388 scm_out_of_range (NULL
, num
);
6391 return scm_to_double (num
);
6395 scm_num2double (SCM num
, unsigned long pos
, const char *s_caller
)
6397 scm_c_issue_deprecation_warning
6398 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
6402 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6406 scm_out_of_range (NULL
, num
);
6409 return scm_to_double (num
);
6415 scm_is_complex (SCM val
)
6417 return scm_is_true (scm_complex_p (val
));
6421 scm_c_real_part (SCM z
)
6423 if (SCM_COMPLEXP (z
))
6424 return SCM_COMPLEX_REAL (z
);
6427 /* Use the scm_real_part to get proper error checking and
6430 return scm_to_double (scm_real_part (z
));
6435 scm_c_imag_part (SCM z
)
6437 if (SCM_COMPLEXP (z
))
6438 return SCM_COMPLEX_IMAG (z
);
6441 /* Use the scm_imag_part to get proper error checking and
6442 dispatching. The result will almost always be 0.0, but not
6445 return scm_to_double (scm_imag_part (z
));
6450 scm_c_magnitude (SCM z
)
6452 return scm_to_double (scm_magnitude (z
));
6458 return scm_to_double (scm_angle (z
));
6462 scm_is_number (SCM z
)
6464 return scm_is_true (scm_number_p (z
));
6468 /* In the following functions we dispatch to the real-arg funcs like log()
6469 when we know the arg is real, instead of just handing everything to
6470 clog() for instance. This is in case clog() doesn't optimize for a
6471 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6472 well use it to go straight to the applicable C func. */
6474 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6476 "Return the natural logarithm of @var{z}.")
6477 #define FUNC_NAME s_scm_log
6479 if (SCM_COMPLEXP (z
))
6481 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6482 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6484 double re
= SCM_COMPLEX_REAL (z
);
6485 double im
= SCM_COMPLEX_IMAG (z
);
6486 return scm_c_make_rectangular (log (hypot (re
, im
)),
6492 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6493 although the value itself overflows. */
6494 double re
= scm_to_double (z
);
6495 double l
= log (fabs (re
));
6497 return scm_from_double (l
);
6499 return scm_c_make_rectangular (l
, M_PI
);
6505 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6507 "Return the base 10 logarithm of @var{z}.")
6508 #define FUNC_NAME s_scm_log10
6510 if (SCM_COMPLEXP (z
))
6512 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6513 clog() and a multiply by M_LOG10E, rather than the fallback
6514 log10+hypot+atan2.) */
6515 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
6516 && defined SCM_COMPLEX_VALUE
6517 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6519 double re
= SCM_COMPLEX_REAL (z
);
6520 double im
= SCM_COMPLEX_IMAG (z
);
6521 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6522 M_LOG10E
* atan2 (im
, re
));
6527 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6528 although the value itself overflows. */
6529 double re
= scm_to_double (z
);
6530 double l
= log10 (fabs (re
));
6532 return scm_from_double (l
);
6534 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6540 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6542 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6543 "base of natural logarithms (2.71828@dots{}).")
6544 #define FUNC_NAME s_scm_exp
6546 if (SCM_COMPLEXP (z
))
6548 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6549 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6551 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6552 SCM_COMPLEX_IMAG (z
));
6557 /* When z is a negative bignum the conversion to double overflows,
6558 giving -infinity, but that's ok, the exp is still 0.0. */
6559 return scm_from_double (exp (scm_to_double (z
)));
6565 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6567 "Return the square root of @var{z}. Of the two possible roots\n"
6568 "(positive and negative), the one with the a positive real part\n"
6569 "is returned, or if that's zero then a positive imaginary part.\n"
6573 "(sqrt 9.0) @result{} 3.0\n"
6574 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6575 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6576 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6578 #define FUNC_NAME s_scm_sqrt
6580 if (SCM_COMPLEXP (x
))
6582 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
6583 && defined SCM_COMPLEX_VALUE
6584 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6586 double re
= SCM_COMPLEX_REAL (x
);
6587 double im
= SCM_COMPLEX_IMAG (x
);
6588 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6589 0.5 * atan2 (im
, re
));
6594 double xx
= scm_to_double (x
);
6596 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6598 return scm_from_double (sqrt (xx
));
6610 mpz_init_set_si (z_negative_one
, -1);
6612 /* It may be possible to tune the performance of some algorithms by using
6613 * the following constants to avoid the creation of bignums. Please, before
6614 * using these values, remember the two rules of program optimization:
6615 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6616 scm_c_define ("most-positive-fixnum",
6617 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6618 scm_c_define ("most-negative-fixnum",
6619 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6621 scm_add_feature ("complex");
6622 scm_add_feature ("inexact");
6623 flo0
= scm_from_double (0.0);
6625 /* determine floating point precision */
6626 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6628 init_dblprec(&scm_dblprec
[i
-2],i
);
6629 init_fx_radix(fx_per_radix
[i
-2],i
);
6632 /* hard code precision for base 10 if the preprocessor tells us to... */
6633 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6636 exactly_one_half
= scm_divide (SCM_INUM1
, SCM_I_MAKINUM (2));
6637 #include "libguile/numbers.x"