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_BIGP() are too large to fit in a fixnum.
26 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
27 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
28 * XXX What about infinities? They are equal to their own floor! -mhw
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))
86 /* On some platforms, isinf(x) returns 0, 1 or -1, indicating the sign
87 of the infinity, but other platforms return a boolean only. */
88 #define DOUBLE_IS_POSITIVE_INFINITY(x) (isinf(x) && ((x) > 0))
89 #define DOUBLE_IS_NEGATIVE_INFINITY(x) (isinf(x) && ((x) < 0))
94 Wonder if this might be faster for some of our code? A switch on
95 the numtag would jump directly to the right case, and the
96 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
98 #define SCM_I_NUMTAG_NOTNUM 0
99 #define SCM_I_NUMTAG_INUM 1
100 #define SCM_I_NUMTAG_BIG scm_tc16_big
101 #define SCM_I_NUMTAG_REAL scm_tc16_real
102 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
103 #define SCM_I_NUMTAG(x) \
104 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
105 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
106 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
107 : SCM_I_NUMTAG_NOTNUM)))
109 /* the macro above will not work as is with fractions */
113 static SCM exactly_one_half
;
115 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
117 /* FLOBUFLEN is the maximum number of characters neccessary for the
118 * printed or scm_string representation of an inexact number.
120 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
123 #if !defined (HAVE_ASINH)
124 static double asinh (double x
) { return log (x
+ sqrt (x
* x
+ 1)); }
126 #if !defined (HAVE_ACOSH)
127 static double acosh (double x
) { return log (x
+ sqrt (x
* x
- 1)); }
129 #if !defined (HAVE_ATANH)
130 static double atanh (double x
) { return 0.5 * log ((1 + x
) / (1 - x
)); }
133 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
134 an explicit check. In some future gmp (don't know what version number),
135 mpz_cmp_d is supposed to do this itself. */
137 #define xmpz_cmp_d(z, d) \
138 (isinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
140 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
144 #if defined (GUILE_I)
145 #if HAVE_COMPLEX_DOUBLE
147 /* For an SCM object Z which is a complex number (ie. satisfies
148 SCM_COMPLEXP), return its value as a C level "complex double". */
149 #define SCM_COMPLEX_VALUE(z) \
150 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
152 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
154 /* Convert a C "complex double" to an SCM value. */
156 scm_from_complex_double (complex double z
)
158 return scm_c_make_rectangular (creal (z
), cimag (z
));
161 #endif /* HAVE_COMPLEX_DOUBLE */
166 static mpz_t z_negative_one
;
169 /* Clear the `mpz_t' embedded in bignum PTR. */
171 finalize_bignum (GC_PTR ptr
, GC_PTR data
)
175 bignum
= PTR2SCM (ptr
);
176 mpz_clear (SCM_I_BIG_MPZ (bignum
));
179 /* Return a new uninitialized bignum. */
184 GC_finalization_proc prev_finalizer
;
185 GC_PTR prev_finalizer_data
;
187 /* Allocate one word for the type tag and enough room for an `mpz_t'. */
188 p
= scm_gc_malloc_pointerless (sizeof (scm_t_bits
) + sizeof (mpz_t
),
192 GC_REGISTER_FINALIZER_NO_ORDER (p
, finalize_bignum
, NULL
,
194 &prev_finalizer_data
);
203 /* Return a newly created bignum. */
204 SCM z
= make_bignum ();
205 mpz_init (SCM_I_BIG_MPZ (z
));
210 scm_i_inum2big (scm_t_inum x
)
212 /* Return a newly created bignum initialized to X. */
213 SCM z
= make_bignum ();
214 #if SIZEOF_VOID_P == SIZEOF_LONG
215 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
217 /* Note that in this case, you'll also have to check all mpz_*_ui and
218 mpz_*_si invocations in Guile. */
219 #error creation of mpz not implemented for this inum size
225 scm_i_long2big (long x
)
227 /* Return a newly created bignum initialized to X. */
228 SCM z
= make_bignum ();
229 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
234 scm_i_ulong2big (unsigned long x
)
236 /* Return a newly created bignum initialized to X. */
237 SCM z
= make_bignum ();
238 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
243 scm_i_clonebig (SCM src_big
, int same_sign_p
)
245 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
246 SCM z
= make_bignum ();
247 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
249 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
254 scm_i_bigcmp (SCM x
, SCM y
)
256 /* Return neg if x < y, pos if x > y, and 0 if x == y */
257 /* presume we already know x and y are bignums */
258 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
259 scm_remember_upto_here_2 (x
, y
);
264 scm_i_dbl2big (double d
)
266 /* results are only defined if d is an integer */
267 SCM z
= make_bignum ();
268 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
272 /* Convert a integer in double representation to a SCM number. */
275 scm_i_dbl2num (double u
)
277 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
278 powers of 2, so there's no rounding when making "double" values
279 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
280 get rounded on a 64-bit machine, hence the "+1".
282 The use of floor() to force to an integer value ensures we get a
283 "numerically closest" value without depending on how a
284 double->long cast or how mpz_set_d will round. For reference,
285 double->long probably follows the hardware rounding mode,
286 mpz_set_d truncates towards zero. */
288 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
289 representable as a double? */
291 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
292 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
293 return SCM_I_MAKINUM ((scm_t_inum
) u
);
295 return scm_i_dbl2big (u
);
298 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
299 with R5RS exact->inexact.
301 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
302 (ie. truncate towards zero), then adjust to get the closest double by
303 examining the next lower bit and adding 1 (to the absolute value) if
306 Bignums exactly half way between representable doubles are rounded to the
307 next higher absolute value (ie. away from zero). This seems like an
308 adequate interpretation of R5RS "numerically closest", and it's easier
309 and faster than a full "nearest-even" style.
311 The bit test must be done on the absolute value of the mpz_t, which means
312 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
313 negatives as twos complement.
315 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
316 following the hardware rounding mode, but applied to the absolute value
317 of the mpz_t operand. This is not what we want so we put the high
318 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
319 mpz_get_d is supposed to always truncate towards zero.
321 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
322 is a slowdown. It'd be faster to pick out the relevant high bits with
323 mpz_getlimbn if we could be bothered coding that, and if the new
324 truncating gmp doesn't come out. */
327 scm_i_big2dbl (SCM b
)
332 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
336 /* Current GMP, eg. 4.1.3, force truncation towards zero */
338 if (bits
> DBL_MANT_DIG
)
340 size_t shift
= bits
- DBL_MANT_DIG
;
341 mpz_init2 (tmp
, DBL_MANT_DIG
);
342 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
343 result
= ldexp (mpz_get_d (tmp
), shift
);
348 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
353 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
356 if (bits
> DBL_MANT_DIG
)
358 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
359 /* test bit number "pos" in absolute value */
360 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
361 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
363 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
367 scm_remember_upto_here_1 (b
);
372 scm_i_normbig (SCM b
)
374 /* convert a big back to a fixnum if it'll fit */
375 /* presume b is a bignum */
376 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
378 scm_t_inum val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
379 if (SCM_FIXABLE (val
))
380 b
= SCM_I_MAKINUM (val
);
385 static SCM_C_INLINE_KEYWORD SCM
386 scm_i_mpz2num (mpz_t b
)
388 /* convert a mpz number to a SCM number. */
389 if (mpz_fits_slong_p (b
))
391 scm_t_inum val
= mpz_get_si (b
);
392 if (SCM_FIXABLE (val
))
393 return SCM_I_MAKINUM (val
);
397 SCM z
= make_bignum ();
398 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
403 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
404 static SCM
scm_divide2real (SCM x
, SCM y
);
407 scm_i_make_ratio (SCM numerator
, SCM denominator
)
408 #define FUNC_NAME "make-ratio"
410 /* First make sure the arguments are proper.
412 if (SCM_I_INUMP (denominator
))
414 if (scm_is_eq (denominator
, SCM_INUM0
))
415 scm_num_overflow ("make-ratio");
416 if (scm_is_eq (denominator
, SCM_INUM1
))
421 if (!(SCM_BIGP(denominator
)))
422 SCM_WRONG_TYPE_ARG (2, denominator
);
424 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
425 SCM_WRONG_TYPE_ARG (1, numerator
);
427 /* Then flip signs so that the denominator is positive.
429 if (scm_is_true (scm_negative_p (denominator
)))
431 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
432 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
435 /* Now consider for each of the four fixnum/bignum combinations
436 whether the rational number is really an integer.
438 if (SCM_I_INUMP (numerator
))
440 scm_t_inum x
= SCM_I_INUM (numerator
);
441 if (scm_is_eq (numerator
, SCM_INUM0
))
443 if (SCM_I_INUMP (denominator
))
446 y
= SCM_I_INUM (denominator
);
450 return SCM_I_MAKINUM (x
/ y
);
454 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
455 of that value for the denominator, as a bignum. Apart from
456 that case, abs(bignum) > abs(inum) so inum/bignum is not an
458 if (x
== SCM_MOST_NEGATIVE_FIXNUM
459 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
460 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
461 return SCM_I_MAKINUM(-1);
464 else if (SCM_BIGP (numerator
))
466 if (SCM_I_INUMP (denominator
))
468 scm_t_inum yy
= SCM_I_INUM (denominator
);
469 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
470 return scm_divide (numerator
, denominator
);
474 if (scm_is_eq (numerator
, denominator
))
476 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
477 SCM_I_BIG_MPZ (denominator
)))
478 return scm_divide(numerator
, denominator
);
482 /* No, it's a proper fraction.
485 SCM divisor
= scm_gcd (numerator
, denominator
);
486 if (!(scm_is_eq (divisor
, SCM_INUM1
)))
488 numerator
= scm_divide (numerator
, divisor
);
489 denominator
= scm_divide (denominator
, divisor
);
492 return scm_double_cell (scm_tc16_fraction
,
493 SCM_UNPACK (numerator
),
494 SCM_UNPACK (denominator
), 0);
500 scm_i_fraction2double (SCM z
)
502 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
503 SCM_FRACTION_DENOMINATOR (z
)));
507 double_is_non_negative_zero (double x
)
509 static double zero
= 0.0;
511 return !memcmp (&x
, &zero
, sizeof(double));
514 SCM_PRIMITIVE_GENERIC (scm_exact_p
, "exact?", 1, 0, 0,
516 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
518 #define FUNC_NAME s_scm_exact_p
520 if (SCM_INEXACTP (x
))
522 else if (SCM_NUMBERP (x
))
525 SCM_WTA_DISPATCH_1 (g_scm_exact_p
, x
, 1, s_scm_exact_p
);
530 SCM_PRIMITIVE_GENERIC (scm_inexact_p
, "inexact?", 1, 0, 0,
532 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
534 #define FUNC_NAME s_scm_inexact_p
536 if (SCM_INEXACTP (x
))
538 else if (SCM_NUMBERP (x
))
541 SCM_WTA_DISPATCH_1 (g_scm_inexact_p
, x
, 1, s_scm_inexact_p
);
546 SCM_PRIMITIVE_GENERIC (scm_odd_p
, "odd?", 1, 0, 0,
548 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
550 #define FUNC_NAME s_scm_odd_p
554 scm_t_inum val
= SCM_I_INUM (n
);
555 return scm_from_bool ((val
& 1L) != 0);
557 else if (SCM_BIGP (n
))
559 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
560 scm_remember_upto_here_1 (n
);
561 return scm_from_bool (odd_p
);
563 else if (SCM_REALP (n
))
565 double val
= SCM_REAL_VALUE (n
);
566 if (DOUBLE_IS_FINITE (val
))
568 double rem
= fabs (fmod (val
, 2.0));
575 SCM_WTA_DISPATCH_1 (g_scm_odd_p
, n
, 1, s_scm_odd_p
);
580 SCM_PRIMITIVE_GENERIC (scm_even_p
, "even?", 1, 0, 0,
582 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
584 #define FUNC_NAME s_scm_even_p
588 scm_t_inum val
= SCM_I_INUM (n
);
589 return scm_from_bool ((val
& 1L) == 0);
591 else if (SCM_BIGP (n
))
593 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
594 scm_remember_upto_here_1 (n
);
595 return scm_from_bool (even_p
);
597 else if (SCM_REALP (n
))
599 double val
= SCM_REAL_VALUE (n
);
600 if (DOUBLE_IS_FINITE (val
))
602 double rem
= fabs (fmod (val
, 2.0));
609 SCM_WTA_DISPATCH_1 (g_scm_even_p
, n
, 1, s_scm_even_p
);
613 SCM_PRIMITIVE_GENERIC (scm_finite_p
, "finite?", 1, 0, 0,
615 "Return @code{#t} if the real number @var{x} is neither\n"
616 "infinite nor a NaN, @code{#f} otherwise.")
617 #define FUNC_NAME s_scm_finite_p
620 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
621 else if (scm_is_real (x
))
624 SCM_WTA_DISPATCH_1 (g_scm_finite_p
, x
, 1, s_scm_finite_p
);
628 SCM_PRIMITIVE_GENERIC (scm_inf_p
, "inf?", 1, 0, 0,
630 "Return @code{#t} if the real number @var{x} is @samp{+inf.0} or\n"
631 "@samp{-inf.0}. Otherwise return @code{#f}.")
632 #define FUNC_NAME s_scm_inf_p
635 return scm_from_bool (isinf (SCM_REAL_VALUE (x
)));
636 else if (scm_is_real (x
))
639 SCM_WTA_DISPATCH_1 (g_scm_inf_p
, x
, 1, s_scm_inf_p
);
643 SCM_PRIMITIVE_GENERIC (scm_nan_p
, "nan?", 1, 0, 0,
645 "Return @code{#t} if the real number @var{x} is a NaN,\n"
646 "or @code{#f} otherwise.")
647 #define FUNC_NAME s_scm_nan_p
650 return scm_from_bool (isnan (SCM_REAL_VALUE (x
)));
651 else if (scm_is_real (x
))
654 SCM_WTA_DISPATCH_1 (g_scm_nan_p
, x
, 1, s_scm_nan_p
);
658 /* Guile's idea of infinity. */
659 static double guile_Inf
;
661 /* Guile's idea of not a number. */
662 static double guile_NaN
;
665 guile_ieee_init (void)
667 /* Some version of gcc on some old version of Linux used to crash when
668 trying to make Inf and NaN. */
671 /* C99 INFINITY, when available.
672 FIXME: The standard allows for INFINITY to be something that overflows
673 at compile time. We ought to have a configure test to check for that
674 before trying to use it. (But in practice we believe this is not a
675 problem on any system guile is likely to target.) */
676 guile_Inf
= INFINITY
;
677 #elif defined HAVE_DINFINITY
679 extern unsigned int DINFINITY
[2];
680 guile_Inf
= (*((double *) (DINFINITY
)));
687 if (guile_Inf
== tmp
)
694 /* C99 NAN, when available */
696 #elif defined HAVE_DQNAN
699 extern unsigned int DQNAN
[2];
700 guile_NaN
= (*((double *)(DQNAN
)));
703 guile_NaN
= guile_Inf
/ guile_Inf
;
707 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
710 #define FUNC_NAME s_scm_inf
712 static int initialized
= 0;
718 return scm_from_double (guile_Inf
);
722 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
725 #define FUNC_NAME s_scm_nan
727 static int initialized
= 0;
733 return scm_from_double (guile_NaN
);
738 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
740 "Return the absolute value of @var{x}.")
741 #define FUNC_NAME s_scm_abs
745 scm_t_inum xx
= SCM_I_INUM (x
);
748 else if (SCM_POSFIXABLE (-xx
))
749 return SCM_I_MAKINUM (-xx
);
751 return scm_i_inum2big (-xx
);
753 else if (SCM_LIKELY (SCM_REALP (x
)))
755 double xx
= SCM_REAL_VALUE (x
);
756 /* If x is a NaN then xx<0 is false so we return x unchanged */
758 return scm_from_double (-xx
);
759 /* Handle signed zeroes properly */
760 else if (SCM_UNLIKELY (xx
== 0.0))
765 else if (SCM_BIGP (x
))
767 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
769 return scm_i_clonebig (x
, 0);
773 else if (SCM_FRACTIONP (x
))
775 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
777 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
778 SCM_FRACTION_DENOMINATOR (x
));
781 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
786 SCM_PRIMITIVE_GENERIC (scm_quotient
, "quotient", 2, 0, 0,
788 "Return the quotient of the numbers @var{x} and @var{y}.")
789 #define FUNC_NAME s_scm_quotient
791 if (SCM_LIKELY (SCM_I_INUMP (x
)))
793 scm_t_inum xx
= SCM_I_INUM (x
);
794 if (SCM_LIKELY (SCM_I_INUMP (y
)))
796 scm_t_inum yy
= SCM_I_INUM (y
);
797 if (SCM_UNLIKELY (yy
== 0))
798 scm_num_overflow (s_scm_quotient
);
801 scm_t_inum z
= xx
/ yy
;
802 if (SCM_LIKELY (SCM_FIXABLE (z
)))
803 return SCM_I_MAKINUM (z
);
805 return scm_i_inum2big (z
);
808 else if (SCM_BIGP (y
))
810 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
811 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
812 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
814 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
815 scm_remember_upto_here_1 (y
);
816 return SCM_I_MAKINUM (-1);
822 SCM_WTA_DISPATCH_2 (g_scm_quotient
, x
, y
, SCM_ARG2
, s_scm_quotient
);
824 else if (SCM_BIGP (x
))
826 if (SCM_LIKELY (SCM_I_INUMP (y
)))
828 scm_t_inum yy
= SCM_I_INUM (y
);
829 if (SCM_UNLIKELY (yy
== 0))
830 scm_num_overflow (s_scm_quotient
);
831 else if (SCM_UNLIKELY (yy
== 1))
835 SCM result
= scm_i_mkbig ();
838 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
841 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
844 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
845 scm_remember_upto_here_1 (x
);
846 return scm_i_normbig (result
);
849 else if (SCM_BIGP (y
))
851 SCM result
= scm_i_mkbig ();
852 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
855 scm_remember_upto_here_2 (x
, y
);
856 return scm_i_normbig (result
);
859 SCM_WTA_DISPATCH_2 (g_scm_quotient
, x
, y
, SCM_ARG2
, s_scm_quotient
);
862 SCM_WTA_DISPATCH_2 (g_scm_quotient
, x
, y
, SCM_ARG1
, s_scm_quotient
);
866 SCM_PRIMITIVE_GENERIC (scm_remainder
, "remainder", 2, 0, 0,
868 "Return the remainder of the numbers @var{x} and @var{y}.\n"
870 "(remainder 13 4) @result{} 1\n"
871 "(remainder -13 4) @result{} -1\n"
873 #define FUNC_NAME s_scm_remainder
875 if (SCM_LIKELY (SCM_I_INUMP (x
)))
877 if (SCM_LIKELY (SCM_I_INUMP (y
)))
879 scm_t_inum yy
= SCM_I_INUM (y
);
880 if (SCM_UNLIKELY (yy
== 0))
881 scm_num_overflow (s_scm_remainder
);
884 /* C99 specifies that "%" is the remainder corresponding to a
885 quotient rounded towards zero, and that's also traditional
886 for machine division, so z here should be well defined. */
887 scm_t_inum z
= SCM_I_INUM (x
) % yy
;
888 return SCM_I_MAKINUM (z
);
891 else if (SCM_BIGP (y
))
893 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
894 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
895 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
897 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
898 scm_remember_upto_here_1 (y
);
905 SCM_WTA_DISPATCH_2 (g_scm_remainder
, x
, y
, SCM_ARG2
, s_scm_remainder
);
907 else if (SCM_BIGP (x
))
909 if (SCM_LIKELY (SCM_I_INUMP (y
)))
911 scm_t_inum yy
= SCM_I_INUM (y
);
912 if (SCM_UNLIKELY (yy
== 0))
913 scm_num_overflow (s_scm_remainder
);
916 SCM result
= scm_i_mkbig ();
919 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
920 scm_remember_upto_here_1 (x
);
921 return scm_i_normbig (result
);
924 else if (SCM_BIGP (y
))
926 SCM result
= scm_i_mkbig ();
927 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
930 scm_remember_upto_here_2 (x
, y
);
931 return scm_i_normbig (result
);
934 SCM_WTA_DISPATCH_2 (g_scm_remainder
, x
, y
, SCM_ARG2
, s_scm_remainder
);
937 SCM_WTA_DISPATCH_2 (g_scm_remainder
, x
, y
, SCM_ARG1
, s_scm_remainder
);
942 SCM_PRIMITIVE_GENERIC (scm_modulo
, "modulo", 2, 0, 0,
944 "Return the modulo of the numbers @var{x} and @var{y}.\n"
946 "(modulo 13 4) @result{} 1\n"
947 "(modulo -13 4) @result{} 3\n"
949 #define FUNC_NAME s_scm_modulo
951 if (SCM_LIKELY (SCM_I_INUMP (x
)))
953 scm_t_inum xx
= SCM_I_INUM (x
);
954 if (SCM_LIKELY (SCM_I_INUMP (y
)))
956 scm_t_inum yy
= SCM_I_INUM (y
);
957 if (SCM_UNLIKELY (yy
== 0))
958 scm_num_overflow (s_scm_modulo
);
961 /* C99 specifies that "%" is the remainder corresponding to a
962 quotient rounded towards zero, and that's also traditional
963 for machine division, so z here should be well defined. */
964 scm_t_inum z
= xx
% yy
;
981 return SCM_I_MAKINUM (result
);
984 else if (SCM_BIGP (y
))
986 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
993 SCM pos_y
= scm_i_clonebig (y
, 0);
994 /* do this after the last scm_op */
995 mpz_init_set_si (z_x
, xx
);
996 result
= pos_y
; /* re-use this bignum */
997 mpz_mod (SCM_I_BIG_MPZ (result
),
999 SCM_I_BIG_MPZ (pos_y
));
1000 scm_remember_upto_here_1 (pos_y
);
1004 result
= scm_i_mkbig ();
1005 /* do this after the last scm_op */
1006 mpz_init_set_si (z_x
, xx
);
1007 mpz_mod (SCM_I_BIG_MPZ (result
),
1010 scm_remember_upto_here_1 (y
);
1013 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1014 mpz_add (SCM_I_BIG_MPZ (result
),
1016 SCM_I_BIG_MPZ (result
));
1017 scm_remember_upto_here_1 (y
);
1018 /* and do this before the next one */
1020 return scm_i_normbig (result
);
1024 SCM_WTA_DISPATCH_2 (g_scm_modulo
, x
, y
, SCM_ARG2
, s_scm_modulo
);
1026 else if (SCM_BIGP (x
))
1028 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1030 scm_t_inum yy
= SCM_I_INUM (y
);
1031 if (SCM_UNLIKELY (yy
== 0))
1032 scm_num_overflow (s_scm_modulo
);
1035 SCM result
= scm_i_mkbig ();
1036 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
1038 (yy
< 0) ? - yy
: yy
);
1039 scm_remember_upto_here_1 (x
);
1040 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1041 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
1042 SCM_I_BIG_MPZ (result
),
1044 return scm_i_normbig (result
);
1047 else if (SCM_BIGP (y
))
1049 SCM result
= scm_i_mkbig ();
1050 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
1051 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1052 mpz_mod (SCM_I_BIG_MPZ (result
),
1054 SCM_I_BIG_MPZ (pos_y
));
1056 scm_remember_upto_here_1 (x
);
1057 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1058 mpz_add (SCM_I_BIG_MPZ (result
),
1060 SCM_I_BIG_MPZ (result
));
1061 scm_remember_upto_here_2 (y
, pos_y
);
1062 return scm_i_normbig (result
);
1065 SCM_WTA_DISPATCH_2 (g_scm_modulo
, x
, y
, SCM_ARG2
, s_scm_modulo
);
1068 SCM_WTA_DISPATCH_2 (g_scm_modulo
, x
, y
, SCM_ARG1
, s_scm_modulo
);
1072 static SCM
scm_i_inexact_euclidean_quotient (double x
, double y
);
1073 static SCM
scm_i_slow_exact_euclidean_quotient (SCM x
, SCM y
);
1075 SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient
, "euclidean-quotient", 2, 0, 0,
1077 "Return the integer @var{q} such that\n"
1078 "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1079 "where @math{0 <= @var{r} < abs(@var{y})}.\n"
1081 "(euclidean-quotient 123 10) @result{} 12\n"
1082 "(euclidean-quotient 123 -10) @result{} -12\n"
1083 "(euclidean-quotient -123 10) @result{} -13\n"
1084 "(euclidean-quotient -123 -10) @result{} 13\n"
1085 "(euclidean-quotient -123.2 -63.5) @result{} 2.0\n"
1086 "(euclidean-quotient 16/3 -10/7) @result{} -3\n"
1088 #define FUNC_NAME s_scm_euclidean_quotient
1090 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1092 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1094 scm_t_inum yy
= SCM_I_INUM (y
);
1095 if (SCM_UNLIKELY (yy
== 0))
1096 scm_num_overflow (s_scm_euclidean_quotient
);
1099 scm_t_inum xx
= SCM_I_INUM (x
);
1100 scm_t_inum qq
= xx
/ yy
;
1108 return SCM_I_MAKINUM (qq
);
1111 else if (SCM_BIGP (y
))
1113 if (SCM_I_INUM (x
) >= 0)
1116 return SCM_I_MAKINUM (- mpz_sgn (SCM_I_BIG_MPZ (y
)));
1118 else if (SCM_REALP (y
))
1119 return scm_i_inexact_euclidean_quotient
1120 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1121 else if (SCM_FRACTIONP (y
))
1122 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1124 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1125 s_scm_euclidean_quotient
);
1127 else if (SCM_BIGP (x
))
1129 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1131 scm_t_inum yy
= SCM_I_INUM (y
);
1132 if (SCM_UNLIKELY (yy
== 0))
1133 scm_num_overflow (s_scm_euclidean_quotient
);
1136 SCM q
= scm_i_mkbig ();
1138 mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (x
), yy
);
1141 mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (x
), -yy
);
1142 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1144 scm_remember_upto_here_1 (x
);
1145 return scm_i_normbig (q
);
1148 else if (SCM_BIGP (y
))
1150 SCM q
= scm_i_mkbig ();
1151 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1152 mpz_fdiv_q (SCM_I_BIG_MPZ (q
),
1156 mpz_cdiv_q (SCM_I_BIG_MPZ (q
),
1159 scm_remember_upto_here_2 (x
, y
);
1160 return scm_i_normbig (q
);
1162 else if (SCM_REALP (y
))
1163 return scm_i_inexact_euclidean_quotient
1164 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1165 else if (SCM_FRACTIONP (y
))
1166 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1168 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1169 s_scm_euclidean_quotient
);
1171 else if (SCM_REALP (x
))
1173 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1174 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1175 return scm_i_inexact_euclidean_quotient
1176 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1178 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1179 s_scm_euclidean_quotient
);
1181 else if (SCM_FRACTIONP (x
))
1184 return scm_i_inexact_euclidean_quotient
1185 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1187 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1190 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG1
,
1191 s_scm_euclidean_quotient
);
1196 scm_i_inexact_euclidean_quotient (double x
, double y
)
1198 if (SCM_LIKELY (y
> 0))
1199 return scm_from_double (floor (x
/ y
));
1200 else if (SCM_LIKELY (y
< 0))
1201 return scm_from_double (ceil (x
/ y
));
1203 scm_num_overflow (s_scm_euclidean_quotient
); /* or return a NaN? */
1208 /* Compute exact euclidean_quotient the slow way.
1209 We use this only if both arguments are exact,
1210 and at least one of them is a fraction */
1212 scm_i_slow_exact_euclidean_quotient (SCM x
, SCM y
)
1214 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1215 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG1
,
1216 s_scm_euclidean_quotient
);
1217 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1218 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1219 s_scm_euclidean_quotient
);
1220 else if (scm_is_true (scm_positive_p (y
)))
1221 return scm_floor (scm_divide (x
, y
));
1222 else if (scm_is_true (scm_negative_p (y
)))
1223 return scm_ceiling (scm_divide (x
, y
));
1225 scm_num_overflow (s_scm_euclidean_quotient
);
1228 static SCM
scm_i_inexact_euclidean_remainder (double x
, double y
);
1229 static SCM
scm_i_slow_exact_euclidean_remainder (SCM x
, SCM y
);
1231 SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder
, "euclidean-remainder", 2, 0, 0,
1233 "Return the real number @var{r} such that\n"
1234 "@math{0 <= @var{r} < abs(@var{y})} and\n"
1235 "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1236 "for some integer @var{q}.\n"
1238 "(euclidean-remainder 123 10) @result{} 3\n"
1239 "(euclidean-remainder 123 -10) @result{} 3\n"
1240 "(euclidean-remainder -123 10) @result{} 7\n"
1241 "(euclidean-remainder -123 -10) @result{} 7\n"
1242 "(euclidean-remainder -123.2 -63.5) @result{} 3.8\n"
1243 "(euclidean-remainder 16/3 -10/7) @result{} 22/21\n"
1245 #define FUNC_NAME s_scm_euclidean_remainder
1247 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1249 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1251 scm_t_inum yy
= SCM_I_INUM (y
);
1252 if (SCM_UNLIKELY (yy
== 0))
1253 scm_num_overflow (s_scm_euclidean_remainder
);
1256 scm_t_inum rr
= SCM_I_INUM (x
) % yy
;
1258 return SCM_I_MAKINUM (rr
);
1260 return SCM_I_MAKINUM (rr
+ yy
);
1262 return SCM_I_MAKINUM (rr
- yy
);
1265 else if (SCM_BIGP (y
))
1267 scm_t_inum xx
= SCM_I_INUM (x
);
1270 else if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1272 SCM r
= scm_i_mkbig ();
1273 mpz_sub_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1274 scm_remember_upto_here_1 (y
);
1275 return scm_i_normbig (r
);
1279 SCM r
= scm_i_mkbig ();
1280 mpz_add_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1281 scm_remember_upto_here_1 (y
);
1282 mpz_neg (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (r
));
1283 return scm_i_normbig (r
);
1286 else if (SCM_REALP (y
))
1287 return scm_i_inexact_euclidean_remainder
1288 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1289 else if (SCM_FRACTIONP (y
))
1290 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1292 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1293 s_scm_euclidean_remainder
);
1295 else if (SCM_BIGP (x
))
1297 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1299 scm_t_inum yy
= SCM_I_INUM (y
);
1300 if (SCM_UNLIKELY (yy
== 0))
1301 scm_num_overflow (s_scm_euclidean_remainder
);
1307 rr
= mpz_fdiv_ui (SCM_I_BIG_MPZ (x
), yy
);
1308 scm_remember_upto_here_1 (x
);
1309 return SCM_I_MAKINUM (rr
);
1312 else if (SCM_BIGP (y
))
1314 SCM r
= scm_i_mkbig ();
1315 mpz_mod (SCM_I_BIG_MPZ (r
),
1318 scm_remember_upto_here_2 (x
, y
);
1319 return scm_i_normbig (r
);
1321 else if (SCM_REALP (y
))
1322 return scm_i_inexact_euclidean_remainder
1323 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1324 else if (SCM_FRACTIONP (y
))
1325 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1327 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1328 s_scm_euclidean_remainder
);
1330 else if (SCM_REALP (x
))
1332 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1333 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1334 return scm_i_inexact_euclidean_remainder
1335 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1337 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1338 s_scm_euclidean_remainder
);
1340 else if (SCM_FRACTIONP (x
))
1343 return scm_i_inexact_euclidean_remainder
1344 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1346 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1349 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG1
,
1350 s_scm_euclidean_remainder
);
1355 scm_i_inexact_euclidean_remainder (double x
, double y
)
1359 /* Although it would be more efficient to use fmod here, we can't
1360 because it would in some cases produce results inconsistent with
1361 scm_i_inexact_euclidean_quotient, such that x != q * y + r (not
1362 even close). In particular, when x is very close to a multiple of
1363 y, then r might be either 0.0 or abs(y)-epsilon, but those two
1364 cases must correspond to different choices of q. If r = 0.0 then q
1365 must be x/y, and if r = abs(y) then q must be (x-r)/y. If quotient
1366 chooses one and remainder chooses the other, it would be bad. This
1367 problem was observed with x = 130.0 and y = 10/7. */
1368 if (SCM_LIKELY (y
> 0))
1370 else if (SCM_LIKELY (y
< 0))
1373 scm_num_overflow (s_scm_euclidean_remainder
); /* or return a NaN? */
1376 return scm_from_double (x
- q
* y
);
1379 /* Compute exact euclidean_remainder the slow way.
1380 We use this only if both arguments are exact,
1381 and at least one of them is a fraction */
1383 scm_i_slow_exact_euclidean_remainder (SCM x
, SCM y
)
1387 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1388 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG1
,
1389 s_scm_euclidean_remainder
);
1390 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1391 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1392 s_scm_euclidean_remainder
);
1393 else if (scm_is_true (scm_positive_p (y
)))
1394 q
= scm_floor (scm_divide (x
, y
));
1395 else if (scm_is_true (scm_negative_p (y
)))
1396 q
= scm_ceiling (scm_divide (x
, y
));
1398 scm_num_overflow (s_scm_euclidean_remainder
);
1399 return scm_difference (x
, scm_product (y
, q
));
1403 static SCM
scm_i_inexact_euclidean_divide (double x
, double y
);
1404 static SCM
scm_i_slow_exact_euclidean_divide (SCM x
, SCM y
);
1406 SCM_PRIMITIVE_GENERIC (scm_euclidean_divide
, "euclidean/", 2, 0, 0,
1408 "Return the integer @var{q} and the real number @var{r}\n"
1409 "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1410 "and @math{0 <= @var{r} < abs(@var{y})}.\n"
1412 "(euclidean/ 123 10) @result{} 12 and 3\n"
1413 "(euclidean/ 123 -10) @result{} -12 and 3\n"
1414 "(euclidean/ -123 10) @result{} -13 and 7\n"
1415 "(euclidean/ -123 -10) @result{} 13 and 7\n"
1416 "(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
1417 "(euclidean/ 16/3 -10/7) @result{} -3 and 22/21\n"
1419 #define FUNC_NAME s_scm_euclidean_divide
1421 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1423 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1425 scm_t_inum yy
= SCM_I_INUM (y
);
1426 if (SCM_UNLIKELY (yy
== 0))
1427 scm_num_overflow (s_scm_euclidean_divide
);
1430 scm_t_inum xx
= SCM_I_INUM (x
);
1431 scm_t_inum qq
= xx
/ yy
;
1432 scm_t_inum rr
= xx
- qq
* yy
;
1440 return scm_values (scm_list_2 (SCM_I_MAKINUM (qq
),
1441 SCM_I_MAKINUM (rr
)));
1444 else if (SCM_BIGP (y
))
1446 scm_t_inum xx
= SCM_I_INUM (x
);
1448 return scm_values (scm_list_2 (SCM_INUM0
, x
));
1449 else if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1451 SCM r
= scm_i_mkbig ();
1452 mpz_sub_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1453 scm_remember_upto_here_1 (y
);
1455 (scm_list_2 (SCM_I_MAKINUM (-1), scm_i_normbig (r
)));
1459 SCM r
= scm_i_mkbig ();
1460 mpz_add_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1461 scm_remember_upto_here_1 (y
);
1462 mpz_neg (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (r
));
1463 return scm_values (scm_list_2 (SCM_INUM1
, scm_i_normbig (r
)));
1466 else if (SCM_REALP (y
))
1467 return scm_i_inexact_euclidean_divide
1468 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1469 else if (SCM_FRACTIONP (y
))
1470 return scm_i_slow_exact_euclidean_divide (x
, y
);
1472 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1473 s_scm_euclidean_divide
);
1475 else if (SCM_BIGP (x
))
1477 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1479 scm_t_inum yy
= SCM_I_INUM (y
);
1480 if (SCM_UNLIKELY (yy
== 0))
1481 scm_num_overflow (s_scm_euclidean_divide
);
1484 SCM q
= scm_i_mkbig ();
1485 SCM r
= scm_i_mkbig ();
1487 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1488 SCM_I_BIG_MPZ (x
), yy
);
1491 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1492 SCM_I_BIG_MPZ (x
), -yy
);
1493 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1495 scm_remember_upto_here_1 (x
);
1496 return scm_values (scm_list_2 (scm_i_normbig (q
),
1497 scm_i_normbig (r
)));
1500 else if (SCM_BIGP (y
))
1502 SCM q
= scm_i_mkbig ();
1503 SCM r
= scm_i_mkbig ();
1504 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1505 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1506 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1508 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1509 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1510 scm_remember_upto_here_2 (x
, y
);
1511 return scm_values (scm_list_2 (scm_i_normbig (q
),
1512 scm_i_normbig (r
)));
1514 else if (SCM_REALP (y
))
1515 return scm_i_inexact_euclidean_divide
1516 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1517 else if (SCM_FRACTIONP (y
))
1518 return scm_i_slow_exact_euclidean_divide (x
, y
);
1520 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1521 s_scm_euclidean_divide
);
1523 else if (SCM_REALP (x
))
1525 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1526 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1527 return scm_i_inexact_euclidean_divide
1528 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1530 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1531 s_scm_euclidean_divide
);
1533 else if (SCM_FRACTIONP (x
))
1536 return scm_i_inexact_euclidean_divide
1537 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1539 return scm_i_slow_exact_euclidean_divide (x
, y
);
1542 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG1
,
1543 s_scm_euclidean_divide
);
1548 scm_i_inexact_euclidean_divide (double x
, double y
)
1552 if (SCM_LIKELY (y
> 0))
1554 else if (SCM_LIKELY (y
< 0))
1557 scm_num_overflow (s_scm_euclidean_divide
); /* or return a NaN? */
1561 return scm_values (scm_list_2 (scm_from_double (q
),
1562 scm_from_double (r
)));
1565 /* Compute exact euclidean quotient and remainder the slow way.
1566 We use this only if both arguments are exact,
1567 and at least one of them is a fraction */
1569 scm_i_slow_exact_euclidean_divide (SCM x
, SCM y
)
1573 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1574 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG1
,
1575 s_scm_euclidean_divide
);
1576 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1577 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1578 s_scm_euclidean_divide
);
1579 else if (scm_is_true (scm_positive_p (y
)))
1580 q
= scm_floor (scm_divide (x
, y
));
1581 else if (scm_is_true (scm_negative_p (y
)))
1582 q
= scm_ceiling (scm_divide (x
, y
));
1584 scm_num_overflow (s_scm_euclidean_divide
);
1585 r
= scm_difference (x
, scm_product (q
, y
));
1586 return scm_values (scm_list_2 (q
, r
));
1589 static SCM
scm_i_inexact_centered_quotient (double x
, double y
);
1590 static SCM
scm_i_bigint_centered_quotient (SCM x
, SCM y
);
1591 static SCM
scm_i_slow_exact_centered_quotient (SCM x
, SCM y
);
1593 SCM_PRIMITIVE_GENERIC (scm_centered_quotient
, "centered-quotient", 2, 0, 0,
1595 "Return the integer @var{q} such that\n"
1596 "@math{@var{x} = @var{q}*@var{y} + @var{r}} where\n"
1597 "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
1599 "(centered-quotient 123 10) @result{} 12\n"
1600 "(centered-quotient 123 -10) @result{} -12\n"
1601 "(centered-quotient -123 10) @result{} -12\n"
1602 "(centered-quotient -123 -10) @result{} 12\n"
1603 "(centered-quotient -123.2 -63.5) @result{} 2.0\n"
1604 "(centered-quotient 16/3 -10/7) @result{} -4\n"
1606 #define FUNC_NAME s_scm_centered_quotient
1608 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1610 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1612 scm_t_inum yy
= SCM_I_INUM (y
);
1613 if (SCM_UNLIKELY (yy
== 0))
1614 scm_num_overflow (s_scm_centered_quotient
);
1617 scm_t_inum xx
= SCM_I_INUM (x
);
1618 scm_t_inum qq
= xx
/ yy
;
1619 scm_t_inum rr
= xx
- qq
* yy
;
1620 if (SCM_LIKELY (xx
> 0))
1622 if (SCM_LIKELY (yy
> 0))
1624 if (rr
>= (yy
+ 1) / 2)
1629 if (rr
>= (1 - yy
) / 2)
1635 if (SCM_LIKELY (yy
> 0))
1646 return SCM_I_MAKINUM (qq
);
1649 else if (SCM_BIGP (y
))
1651 /* Pass a denormalized bignum version of x (even though it
1652 can fit in a fixnum) to scm_i_bigint_centered_quotient */
1653 return scm_i_bigint_centered_quotient
1654 (scm_i_long2big (SCM_I_INUM (x
)), y
);
1656 else if (SCM_REALP (y
))
1657 return scm_i_inexact_centered_quotient
1658 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1659 else if (SCM_FRACTIONP (y
))
1660 return scm_i_slow_exact_centered_quotient (x
, y
);
1662 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1663 s_scm_centered_quotient
);
1665 else if (SCM_BIGP (x
))
1667 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1669 scm_t_inum yy
= SCM_I_INUM (y
);
1670 if (SCM_UNLIKELY (yy
== 0))
1671 scm_num_overflow (s_scm_centered_quotient
);
1674 SCM q
= scm_i_mkbig ();
1676 /* Arrange for rr to initially be non-positive,
1677 because that simplifies the test to see
1678 if it is within the needed bounds. */
1681 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
1682 SCM_I_BIG_MPZ (x
), yy
);
1683 scm_remember_upto_here_1 (x
);
1685 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
1686 SCM_I_BIG_MPZ (q
), 1);
1690 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
1691 SCM_I_BIG_MPZ (x
), -yy
);
1692 scm_remember_upto_here_1 (x
);
1693 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1695 mpz_add_ui (SCM_I_BIG_MPZ (q
),
1696 SCM_I_BIG_MPZ (q
), 1);
1698 return scm_i_normbig (q
);
1701 else if (SCM_BIGP (y
))
1702 return scm_i_bigint_centered_quotient (x
, y
);
1703 else if (SCM_REALP (y
))
1704 return scm_i_inexact_centered_quotient
1705 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1706 else if (SCM_FRACTIONP (y
))
1707 return scm_i_slow_exact_centered_quotient (x
, y
);
1709 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1710 s_scm_centered_quotient
);
1712 else if (SCM_REALP (x
))
1714 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1715 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1716 return scm_i_inexact_centered_quotient
1717 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1719 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1720 s_scm_centered_quotient
);
1722 else if (SCM_FRACTIONP (x
))
1725 return scm_i_inexact_centered_quotient
1726 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1728 return scm_i_slow_exact_centered_quotient (x
, y
);
1731 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG1
,
1732 s_scm_centered_quotient
);
1737 scm_i_inexact_centered_quotient (double x
, double y
)
1739 if (SCM_LIKELY (y
> 0))
1740 return scm_from_double (floor (x
/y
+ 0.5));
1741 else if (SCM_LIKELY (y
< 0))
1742 return scm_from_double (ceil (x
/y
- 0.5));
1744 scm_num_overflow (s_scm_centered_quotient
); /* or return a NaN? */
1749 /* Assumes that both x and y are bigints, though
1750 x might be able to fit into a fixnum. */
1752 scm_i_bigint_centered_quotient (SCM x
, SCM y
)
1756 /* Note that x might be small enough to fit into a
1757 fixnum, so we must not let it escape into the wild */
1761 /* min_r will eventually become -abs(y)/2 */
1762 min_r
= scm_i_mkbig ();
1763 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
1764 SCM_I_BIG_MPZ (y
), 1);
1766 /* Arrange for rr to initially be non-positive,
1767 because that simplifies the test to see
1768 if it is within the needed bounds. */
1769 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1771 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1772 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1773 scm_remember_upto_here_2 (x
, y
);
1774 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
1775 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1776 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
1777 SCM_I_BIG_MPZ (q
), 1);
1781 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1782 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1783 scm_remember_upto_here_2 (x
, y
);
1784 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1785 mpz_add_ui (SCM_I_BIG_MPZ (q
),
1786 SCM_I_BIG_MPZ (q
), 1);
1788 scm_remember_upto_here_2 (r
, min_r
);
1789 return scm_i_normbig (q
);
1792 /* Compute exact centered quotient the slow way.
1793 We use this only if both arguments are exact,
1794 and at least one of them is a fraction */
1796 scm_i_slow_exact_centered_quotient (SCM x
, SCM y
)
1798 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1799 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG1
,
1800 s_scm_centered_quotient
);
1801 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1802 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1803 s_scm_centered_quotient
);
1804 else if (scm_is_true (scm_positive_p (y
)))
1805 return scm_floor (scm_sum (scm_divide (x
, y
),
1807 else if (scm_is_true (scm_negative_p (y
)))
1808 return scm_ceiling (scm_difference (scm_divide (x
, y
),
1811 scm_num_overflow (s_scm_centered_quotient
);
1814 static SCM
scm_i_inexact_centered_remainder (double x
, double y
);
1815 static SCM
scm_i_bigint_centered_remainder (SCM x
, SCM y
);
1816 static SCM
scm_i_slow_exact_centered_remainder (SCM x
, SCM y
);
1818 SCM_PRIMITIVE_GENERIC (scm_centered_remainder
, "centered-remainder", 2, 0, 0,
1820 "Return the real number @var{r} such that\n"
1821 "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}\n"
1822 "and @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1823 "for some integer @var{q}.\n"
1825 "(centered-remainder 123 10) @result{} 3\n"
1826 "(centered-remainder 123 -10) @result{} 3\n"
1827 "(centered-remainder -123 10) @result{} -3\n"
1828 "(centered-remainder -123 -10) @result{} -3\n"
1829 "(centered-remainder -123.2 -63.5) @result{} 3.8\n"
1830 "(centered-remainder 16/3 -10/7) @result{} -8/21\n"
1832 #define FUNC_NAME s_scm_centered_remainder
1834 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1836 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1838 scm_t_inum yy
= SCM_I_INUM (y
);
1839 if (SCM_UNLIKELY (yy
== 0))
1840 scm_num_overflow (s_scm_centered_remainder
);
1843 scm_t_inum xx
= SCM_I_INUM (x
);
1844 scm_t_inum rr
= xx
% yy
;
1845 if (SCM_LIKELY (xx
> 0))
1847 if (SCM_LIKELY (yy
> 0))
1849 if (rr
>= (yy
+ 1) / 2)
1854 if (rr
>= (1 - yy
) / 2)
1860 if (SCM_LIKELY (yy
> 0))
1871 return SCM_I_MAKINUM (rr
);
1874 else if (SCM_BIGP (y
))
1876 /* Pass a denormalized bignum version of x (even though it
1877 can fit in a fixnum) to scm_i_bigint_centered_remainder */
1878 return scm_i_bigint_centered_remainder
1879 (scm_i_long2big (SCM_I_INUM (x
)), y
);
1881 else if (SCM_REALP (y
))
1882 return scm_i_inexact_centered_remainder
1883 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1884 else if (SCM_FRACTIONP (y
))
1885 return scm_i_slow_exact_centered_remainder (x
, y
);
1887 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1888 s_scm_centered_remainder
);
1890 else if (SCM_BIGP (x
))
1892 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1894 scm_t_inum yy
= SCM_I_INUM (y
);
1895 if (SCM_UNLIKELY (yy
== 0))
1896 scm_num_overflow (s_scm_centered_remainder
);
1900 /* Arrange for rr to initially be non-positive,
1901 because that simplifies the test to see
1902 if it is within the needed bounds. */
1905 rr
= - mpz_cdiv_ui (SCM_I_BIG_MPZ (x
), yy
);
1906 scm_remember_upto_here_1 (x
);
1912 rr
= - mpz_cdiv_ui (SCM_I_BIG_MPZ (x
), -yy
);
1913 scm_remember_upto_here_1 (x
);
1917 return SCM_I_MAKINUM (rr
);
1920 else if (SCM_BIGP (y
))
1921 return scm_i_bigint_centered_remainder (x
, y
);
1922 else if (SCM_REALP (y
))
1923 return scm_i_inexact_centered_remainder
1924 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1925 else if (SCM_FRACTIONP (y
))
1926 return scm_i_slow_exact_centered_remainder (x
, y
);
1928 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1929 s_scm_centered_remainder
);
1931 else if (SCM_REALP (x
))
1933 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1934 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1935 return scm_i_inexact_centered_remainder
1936 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1938 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1939 s_scm_centered_remainder
);
1941 else if (SCM_FRACTIONP (x
))
1944 return scm_i_inexact_centered_remainder
1945 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1947 return scm_i_slow_exact_centered_remainder (x
, y
);
1950 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG1
,
1951 s_scm_centered_remainder
);
1956 scm_i_inexact_centered_remainder (double x
, double y
)
1960 /* Although it would be more efficient to use fmod here, we can't
1961 because it would in some cases produce results inconsistent with
1962 scm_i_inexact_centered_quotient, such that x != r + q * y (not even
1963 close). In particular, when x-y/2 is very close to a multiple of
1964 y, then r might be either -abs(y/2) or abs(y/2)-epsilon, but those
1965 two cases must correspond to different choices of q. If quotient
1966 chooses one and remainder chooses the other, it would be bad. */
1967 if (SCM_LIKELY (y
> 0))
1968 q
= floor (x
/y
+ 0.5);
1969 else if (SCM_LIKELY (y
< 0))
1970 q
= ceil (x
/y
- 0.5);
1972 scm_num_overflow (s_scm_centered_remainder
); /* or return a NaN? */
1975 return scm_from_double (x
- q
* y
);
1978 /* Assumes that both x and y are bigints, though
1979 x might be able to fit into a fixnum. */
1981 scm_i_bigint_centered_remainder (SCM x
, SCM y
)
1985 /* Note that x might be small enough to fit into a
1986 fixnum, so we must not let it escape into the wild */
1989 /* min_r will eventually become -abs(y)/2 */
1990 min_r
= scm_i_mkbig ();
1991 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
1992 SCM_I_BIG_MPZ (y
), 1);
1994 /* Arrange for rr to initially be non-positive,
1995 because that simplifies the test to see
1996 if it is within the needed bounds. */
1997 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1999 mpz_cdiv_r (SCM_I_BIG_MPZ (r
),
2000 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2001 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
2002 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2003 mpz_add (SCM_I_BIG_MPZ (r
),
2009 mpz_fdiv_r (SCM_I_BIG_MPZ (r
),
2010 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2011 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2012 mpz_sub (SCM_I_BIG_MPZ (r
),
2016 scm_remember_upto_here_2 (x
, y
);
2017 return scm_i_normbig (r
);
2020 /* Compute exact centered_remainder the slow way.
2021 We use this only if both arguments are exact,
2022 and at least one of them is a fraction */
2024 scm_i_slow_exact_centered_remainder (SCM x
, SCM y
)
2028 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
2029 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG1
,
2030 s_scm_centered_remainder
);
2031 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
2032 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
2033 s_scm_centered_remainder
);
2034 else if (scm_is_true (scm_positive_p (y
)))
2035 q
= scm_floor (scm_sum (scm_divide (x
, y
), exactly_one_half
));
2036 else if (scm_is_true (scm_negative_p (y
)))
2037 q
= scm_ceiling (scm_difference (scm_divide (x
, y
), exactly_one_half
));
2039 scm_num_overflow (s_scm_centered_remainder
);
2040 return scm_difference (x
, scm_product (y
, q
));
2044 static SCM
scm_i_inexact_centered_divide (double x
, double y
);
2045 static SCM
scm_i_bigint_centered_divide (SCM x
, SCM y
);
2046 static SCM
scm_i_slow_exact_centered_divide (SCM x
, SCM y
);
2048 SCM_PRIMITIVE_GENERIC (scm_centered_divide
, "centered/", 2, 0, 0,
2050 "Return the integer @var{q} and the real number @var{r}\n"
2051 "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
2052 "and @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
2054 "(centered/ 123 10) @result{} 12 and 3\n"
2055 "(centered/ 123 -10) @result{} -12 and 3\n"
2056 "(centered/ -123 10) @result{} -12 and -3\n"
2057 "(centered/ -123 -10) @result{} 12 and -3\n"
2058 "(centered/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
2059 "(centered/ 16/3 -10/7) @result{} -4 and -8/21\n"
2061 #define FUNC_NAME s_scm_centered_divide
2063 if (SCM_LIKELY (SCM_I_INUMP (x
)))
2065 if (SCM_LIKELY (SCM_I_INUMP (y
)))
2067 scm_t_inum yy
= SCM_I_INUM (y
);
2068 if (SCM_UNLIKELY (yy
== 0))
2069 scm_num_overflow (s_scm_centered_divide
);
2072 scm_t_inum xx
= SCM_I_INUM (x
);
2073 scm_t_inum qq
= xx
/ yy
;
2074 scm_t_inum rr
= xx
- qq
* yy
;
2075 if (SCM_LIKELY (xx
> 0))
2077 if (SCM_LIKELY (yy
> 0))
2079 if (rr
>= (yy
+ 1) / 2)
2084 if (rr
>= (1 - yy
) / 2)
2090 if (SCM_LIKELY (yy
> 0))
2101 return scm_values (scm_list_2 (SCM_I_MAKINUM (qq
),
2102 SCM_I_MAKINUM (rr
)));
2105 else if (SCM_BIGP (y
))
2107 /* Pass a denormalized bignum version of x (even though it
2108 can fit in a fixnum) to scm_i_bigint_centered_divide */
2109 return scm_i_bigint_centered_divide
2110 (scm_i_long2big (SCM_I_INUM (x
)), y
);
2112 else if (SCM_REALP (y
))
2113 return scm_i_inexact_centered_divide
2114 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
2115 else if (SCM_FRACTIONP (y
))
2116 return scm_i_slow_exact_centered_divide (x
, y
);
2118 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2119 s_scm_centered_divide
);
2121 else if (SCM_BIGP (x
))
2123 if (SCM_LIKELY (SCM_I_INUMP (y
)))
2125 scm_t_inum yy
= SCM_I_INUM (y
);
2126 if (SCM_UNLIKELY (yy
== 0))
2127 scm_num_overflow (s_scm_centered_divide
);
2130 SCM q
= scm_i_mkbig ();
2132 /* Arrange for rr to initially be non-positive,
2133 because that simplifies the test to see
2134 if it is within the needed bounds. */
2137 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
2138 SCM_I_BIG_MPZ (x
), yy
);
2139 scm_remember_upto_here_1 (x
);
2142 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
2143 SCM_I_BIG_MPZ (q
), 1);
2149 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
2150 SCM_I_BIG_MPZ (x
), -yy
);
2151 scm_remember_upto_here_1 (x
);
2152 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
2155 mpz_add_ui (SCM_I_BIG_MPZ (q
),
2156 SCM_I_BIG_MPZ (q
), 1);
2160 return scm_values (scm_list_2 (scm_i_normbig (q
),
2161 SCM_I_MAKINUM (rr
)));
2164 else if (SCM_BIGP (y
))
2165 return scm_i_bigint_centered_divide (x
, y
);
2166 else if (SCM_REALP (y
))
2167 return scm_i_inexact_centered_divide
2168 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
2169 else if (SCM_FRACTIONP (y
))
2170 return scm_i_slow_exact_centered_divide (x
, y
);
2172 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2173 s_scm_centered_divide
);
2175 else if (SCM_REALP (x
))
2177 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
2178 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
2179 return scm_i_inexact_centered_divide
2180 (SCM_REAL_VALUE (x
), scm_to_double (y
));
2182 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2183 s_scm_centered_divide
);
2185 else if (SCM_FRACTIONP (x
))
2188 return scm_i_inexact_centered_divide
2189 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
2191 return scm_i_slow_exact_centered_divide (x
, y
);
2194 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG1
,
2195 s_scm_centered_divide
);
2200 scm_i_inexact_centered_divide (double x
, double y
)
2204 if (SCM_LIKELY (y
> 0))
2205 q
= floor (x
/y
+ 0.5);
2206 else if (SCM_LIKELY (y
< 0))
2207 q
= ceil (x
/y
- 0.5);
2209 scm_num_overflow (s_scm_centered_divide
); /* or return a NaN? */
2213 return scm_values (scm_list_2 (scm_from_double (q
),
2214 scm_from_double (r
)));
2217 /* Assumes that both x and y are bigints, though
2218 x might be able to fit into a fixnum. */
2220 scm_i_bigint_centered_divide (SCM x
, SCM y
)
2224 /* Note that x might be small enough to fit into a
2225 fixnum, so we must not let it escape into the wild */
2229 /* min_r will eventually become -abs(y/2) */
2230 min_r
= scm_i_mkbig ();
2231 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
2232 SCM_I_BIG_MPZ (y
), 1);
2234 /* Arrange for rr to initially be non-positive,
2235 because that simplifies the test to see
2236 if it is within the needed bounds. */
2237 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
2239 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
2240 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2241 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
2242 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2244 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
2245 SCM_I_BIG_MPZ (q
), 1);
2246 mpz_add (SCM_I_BIG_MPZ (r
),
2253 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
2254 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2255 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2257 mpz_add_ui (SCM_I_BIG_MPZ (q
),
2258 SCM_I_BIG_MPZ (q
), 1);
2259 mpz_sub (SCM_I_BIG_MPZ (r
),
2264 scm_remember_upto_here_2 (x
, y
);
2265 return scm_values (scm_list_2 (scm_i_normbig (q
),
2266 scm_i_normbig (r
)));
2269 /* Compute exact centered quotient and remainder the slow way.
2270 We use this only if both arguments are exact,
2271 and at least one of them is a fraction */
2273 scm_i_slow_exact_centered_divide (SCM x
, SCM y
)
2277 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
2278 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG1
,
2279 s_scm_centered_divide
);
2280 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
2281 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2282 s_scm_centered_divide
);
2283 else if (scm_is_true (scm_positive_p (y
)))
2284 q
= scm_floor (scm_sum (scm_divide (x
, y
),
2286 else if (scm_is_true (scm_negative_p (y
)))
2287 q
= scm_ceiling (scm_difference (scm_divide (x
, y
),
2290 scm_num_overflow (s_scm_centered_divide
);
2291 r
= scm_difference (x
, scm_product (q
, y
));
2292 return scm_values (scm_list_2 (q
, r
));
2296 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
2297 (SCM x
, SCM y
, SCM rest
),
2298 "Return the greatest common divisor of all parameter values.\n"
2299 "If called without arguments, 0 is returned.")
2300 #define FUNC_NAME s_scm_i_gcd
2302 while (!scm_is_null (rest
))
2303 { x
= scm_gcd (x
, y
);
2305 rest
= scm_cdr (rest
);
2307 return scm_gcd (x
, y
);
2311 #define s_gcd s_scm_i_gcd
2312 #define g_gcd g_scm_i_gcd
2315 scm_gcd (SCM x
, SCM y
)
2318 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
2320 if (SCM_I_INUMP (x
))
2322 if (SCM_I_INUMP (y
))
2324 scm_t_inum xx
= SCM_I_INUM (x
);
2325 scm_t_inum yy
= SCM_I_INUM (y
);
2326 scm_t_inum u
= xx
< 0 ? -xx
: xx
;
2327 scm_t_inum v
= yy
< 0 ? -yy
: yy
;
2337 /* Determine a common factor 2^k */
2338 while (!(1 & (u
| v
)))
2344 /* Now, any factor 2^n can be eliminated */
2364 return (SCM_POSFIXABLE (result
)
2365 ? SCM_I_MAKINUM (result
)
2366 : scm_i_inum2big (result
));
2368 else if (SCM_BIGP (y
))
2374 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
2376 else if (SCM_BIGP (x
))
2378 if (SCM_I_INUMP (y
))
2383 yy
= SCM_I_INUM (y
);
2388 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
2389 scm_remember_upto_here_1 (x
);
2390 return (SCM_POSFIXABLE (result
)
2391 ? SCM_I_MAKINUM (result
)
2392 : scm_from_unsigned_integer (result
));
2394 else if (SCM_BIGP (y
))
2396 SCM result
= scm_i_mkbig ();
2397 mpz_gcd (SCM_I_BIG_MPZ (result
),
2400 scm_remember_upto_here_2 (x
, y
);
2401 return scm_i_normbig (result
);
2404 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
2407 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
2410 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
2411 (SCM x
, SCM y
, SCM rest
),
2412 "Return the least common multiple of the arguments.\n"
2413 "If called without arguments, 1 is returned.")
2414 #define FUNC_NAME s_scm_i_lcm
2416 while (!scm_is_null (rest
))
2417 { x
= scm_lcm (x
, y
);
2419 rest
= scm_cdr (rest
);
2421 return scm_lcm (x
, y
);
2425 #define s_lcm s_scm_i_lcm
2426 #define g_lcm g_scm_i_lcm
2429 scm_lcm (SCM n1
, SCM n2
)
2431 if (SCM_UNBNDP (n2
))
2433 if (SCM_UNBNDP (n1
))
2434 return SCM_I_MAKINUM (1L);
2435 n2
= SCM_I_MAKINUM (1L);
2438 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
2439 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
2440 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
2441 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
2443 if (SCM_I_INUMP (n1
))
2445 if (SCM_I_INUMP (n2
))
2447 SCM d
= scm_gcd (n1
, n2
);
2448 if (scm_is_eq (d
, SCM_INUM0
))
2451 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
2455 /* inum n1, big n2 */
2458 SCM result
= scm_i_mkbig ();
2459 scm_t_inum nn1
= SCM_I_INUM (n1
);
2460 if (nn1
== 0) return SCM_INUM0
;
2461 if (nn1
< 0) nn1
= - nn1
;
2462 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
2463 scm_remember_upto_here_1 (n2
);
2471 if (SCM_I_INUMP (n2
))
2478 SCM result
= scm_i_mkbig ();
2479 mpz_lcm(SCM_I_BIG_MPZ (result
),
2481 SCM_I_BIG_MPZ (n2
));
2482 scm_remember_upto_here_2(n1
, n2
);
2483 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
2489 /* Emulating 2's complement bignums with sign magnitude arithmetic:
2494 + + + x (map digit:logand X Y)
2495 + - + x (map digit:logand X (lognot (+ -1 Y)))
2496 - + + y (map digit:logand (lognot (+ -1 X)) Y)
2497 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
2502 + + + (map digit:logior X Y)
2503 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
2504 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
2505 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
2510 + + + (map digit:logxor X Y)
2511 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
2512 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
2513 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
2518 + + (any digit:logand X Y)
2519 + - (any digit:logand X (lognot (+ -1 Y)))
2520 - + (any digit:logand (lognot (+ -1 X)) Y)
2525 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
2526 (SCM x
, SCM y
, SCM rest
),
2527 "Return the bitwise AND of the integer arguments.\n\n"
2529 "(logand) @result{} -1\n"
2530 "(logand 7) @result{} 7\n"
2531 "(logand #b111 #b011 #b001) @result{} 1\n"
2533 #define FUNC_NAME s_scm_i_logand
2535 while (!scm_is_null (rest
))
2536 { x
= scm_logand (x
, y
);
2538 rest
= scm_cdr (rest
);
2540 return scm_logand (x
, y
);
2544 #define s_scm_logand s_scm_i_logand
2546 SCM
scm_logand (SCM n1
, SCM n2
)
2547 #define FUNC_NAME s_scm_logand
2551 if (SCM_UNBNDP (n2
))
2553 if (SCM_UNBNDP (n1
))
2554 return SCM_I_MAKINUM (-1);
2555 else if (!SCM_NUMBERP (n1
))
2556 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2557 else if (SCM_NUMBERP (n1
))
2560 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2563 if (SCM_I_INUMP (n1
))
2565 nn1
= SCM_I_INUM (n1
);
2566 if (SCM_I_INUMP (n2
))
2568 scm_t_inum nn2
= SCM_I_INUM (n2
);
2569 return SCM_I_MAKINUM (nn1
& nn2
);
2571 else if SCM_BIGP (n2
)
2577 SCM result_z
= scm_i_mkbig ();
2579 mpz_init_set_si (nn1_z
, nn1
);
2580 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2581 scm_remember_upto_here_1 (n2
);
2583 return scm_i_normbig (result_z
);
2587 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2589 else if (SCM_BIGP (n1
))
2591 if (SCM_I_INUMP (n2
))
2594 nn1
= SCM_I_INUM (n1
);
2597 else if (SCM_BIGP (n2
))
2599 SCM result_z
= scm_i_mkbig ();
2600 mpz_and (SCM_I_BIG_MPZ (result_z
),
2602 SCM_I_BIG_MPZ (n2
));
2603 scm_remember_upto_here_2 (n1
, n2
);
2604 return scm_i_normbig (result_z
);
2607 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2610 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2615 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
2616 (SCM x
, SCM y
, SCM rest
),
2617 "Return the bitwise OR of the integer arguments.\n\n"
2619 "(logior) @result{} 0\n"
2620 "(logior 7) @result{} 7\n"
2621 "(logior #b000 #b001 #b011) @result{} 3\n"
2623 #define FUNC_NAME s_scm_i_logior
2625 while (!scm_is_null (rest
))
2626 { x
= scm_logior (x
, y
);
2628 rest
= scm_cdr (rest
);
2630 return scm_logior (x
, y
);
2634 #define s_scm_logior s_scm_i_logior
2636 SCM
scm_logior (SCM n1
, SCM n2
)
2637 #define FUNC_NAME s_scm_logior
2641 if (SCM_UNBNDP (n2
))
2643 if (SCM_UNBNDP (n1
))
2645 else if (SCM_NUMBERP (n1
))
2648 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2651 if (SCM_I_INUMP (n1
))
2653 nn1
= SCM_I_INUM (n1
);
2654 if (SCM_I_INUMP (n2
))
2656 long nn2
= SCM_I_INUM (n2
);
2657 return SCM_I_MAKINUM (nn1
| nn2
);
2659 else if (SCM_BIGP (n2
))
2665 SCM result_z
= scm_i_mkbig ();
2667 mpz_init_set_si (nn1_z
, nn1
);
2668 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2669 scm_remember_upto_here_1 (n2
);
2671 return scm_i_normbig (result_z
);
2675 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2677 else if (SCM_BIGP (n1
))
2679 if (SCM_I_INUMP (n2
))
2682 nn1
= SCM_I_INUM (n1
);
2685 else if (SCM_BIGP (n2
))
2687 SCM result_z
= scm_i_mkbig ();
2688 mpz_ior (SCM_I_BIG_MPZ (result_z
),
2690 SCM_I_BIG_MPZ (n2
));
2691 scm_remember_upto_here_2 (n1
, n2
);
2692 return scm_i_normbig (result_z
);
2695 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2698 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2703 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
2704 (SCM x
, SCM y
, SCM rest
),
2705 "Return the bitwise XOR of the integer arguments. A bit is\n"
2706 "set in the result if it is set in an odd number of arguments.\n"
2708 "(logxor) @result{} 0\n"
2709 "(logxor 7) @result{} 7\n"
2710 "(logxor #b000 #b001 #b011) @result{} 2\n"
2711 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
2713 #define FUNC_NAME s_scm_i_logxor
2715 while (!scm_is_null (rest
))
2716 { x
= scm_logxor (x
, y
);
2718 rest
= scm_cdr (rest
);
2720 return scm_logxor (x
, y
);
2724 #define s_scm_logxor s_scm_i_logxor
2726 SCM
scm_logxor (SCM n1
, SCM n2
)
2727 #define FUNC_NAME s_scm_logxor
2731 if (SCM_UNBNDP (n2
))
2733 if (SCM_UNBNDP (n1
))
2735 else if (SCM_NUMBERP (n1
))
2738 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2741 if (SCM_I_INUMP (n1
))
2743 nn1
= SCM_I_INUM (n1
);
2744 if (SCM_I_INUMP (n2
))
2746 scm_t_inum nn2
= SCM_I_INUM (n2
);
2747 return SCM_I_MAKINUM (nn1
^ nn2
);
2749 else if (SCM_BIGP (n2
))
2753 SCM result_z
= scm_i_mkbig ();
2755 mpz_init_set_si (nn1_z
, nn1
);
2756 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2757 scm_remember_upto_here_1 (n2
);
2759 return scm_i_normbig (result_z
);
2763 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2765 else if (SCM_BIGP (n1
))
2767 if (SCM_I_INUMP (n2
))
2770 nn1
= SCM_I_INUM (n1
);
2773 else if (SCM_BIGP (n2
))
2775 SCM result_z
= scm_i_mkbig ();
2776 mpz_xor (SCM_I_BIG_MPZ (result_z
),
2778 SCM_I_BIG_MPZ (n2
));
2779 scm_remember_upto_here_2 (n1
, n2
);
2780 return scm_i_normbig (result_z
);
2783 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2786 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2791 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
2793 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
2794 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
2795 "without actually calculating the @code{logand}, just testing\n"
2799 "(logtest #b0100 #b1011) @result{} #f\n"
2800 "(logtest #b0100 #b0111) @result{} #t\n"
2802 #define FUNC_NAME s_scm_logtest
2806 if (SCM_I_INUMP (j
))
2808 nj
= SCM_I_INUM (j
);
2809 if (SCM_I_INUMP (k
))
2811 scm_t_inum nk
= SCM_I_INUM (k
);
2812 return scm_from_bool (nj
& nk
);
2814 else if (SCM_BIGP (k
))
2822 mpz_init_set_si (nj_z
, nj
);
2823 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
2824 scm_remember_upto_here_1 (k
);
2825 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
2831 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
2833 else if (SCM_BIGP (j
))
2835 if (SCM_I_INUMP (k
))
2838 nj
= SCM_I_INUM (j
);
2841 else if (SCM_BIGP (k
))
2845 mpz_init (result_z
);
2849 scm_remember_upto_here_2 (j
, k
);
2850 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
2851 mpz_clear (result_z
);
2855 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
2858 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
2863 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
2865 "Test whether bit number @var{index} in @var{j} is set.\n"
2866 "@var{index} starts from 0 for the least significant bit.\n"
2869 "(logbit? 0 #b1101) @result{} #t\n"
2870 "(logbit? 1 #b1101) @result{} #f\n"
2871 "(logbit? 2 #b1101) @result{} #t\n"
2872 "(logbit? 3 #b1101) @result{} #t\n"
2873 "(logbit? 4 #b1101) @result{} #f\n"
2875 #define FUNC_NAME s_scm_logbit_p
2877 unsigned long int iindex
;
2878 iindex
= scm_to_ulong (index
);
2880 if (SCM_I_INUMP (j
))
2882 /* bits above what's in an inum follow the sign bit */
2883 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
2884 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
2886 else if (SCM_BIGP (j
))
2888 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
2889 scm_remember_upto_here_1 (j
);
2890 return scm_from_bool (val
);
2893 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
2898 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
2900 "Return the integer which is the ones-complement of the integer\n"
2904 "(number->string (lognot #b10000000) 2)\n"
2905 " @result{} \"-10000001\"\n"
2906 "(number->string (lognot #b0) 2)\n"
2907 " @result{} \"-1\"\n"
2909 #define FUNC_NAME s_scm_lognot
2911 if (SCM_I_INUMP (n
)) {
2912 /* No overflow here, just need to toggle all the bits making up the inum.
2913 Enhancement: No need to strip the tag and add it back, could just xor
2914 a block of 1 bits, if that worked with the various debug versions of
2916 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
2918 } else if (SCM_BIGP (n
)) {
2919 SCM result
= scm_i_mkbig ();
2920 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
2921 scm_remember_upto_here_1 (n
);
2925 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2930 /* returns 0 if IN is not an integer. OUT must already be
2933 coerce_to_big (SCM in
, mpz_t out
)
2936 mpz_set (out
, SCM_I_BIG_MPZ (in
));
2937 else if (SCM_I_INUMP (in
))
2938 mpz_set_si (out
, SCM_I_INUM (in
));
2945 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
2946 (SCM n
, SCM k
, SCM m
),
2947 "Return @var{n} raised to the integer exponent\n"
2948 "@var{k}, modulo @var{m}.\n"
2951 "(modulo-expt 2 3 5)\n"
2954 #define FUNC_NAME s_scm_modulo_expt
2960 /* There are two classes of error we might encounter --
2961 1) Math errors, which we'll report by calling scm_num_overflow,
2963 2) wrong-type errors, which of course we'll report by calling
2965 We don't report those errors immediately, however; instead we do
2966 some cleanup first. These variables tell us which error (if
2967 any) we should report after cleaning up.
2969 int report_overflow
= 0;
2971 int position_of_wrong_type
= 0;
2972 SCM value_of_wrong_type
= SCM_INUM0
;
2974 SCM result
= SCM_UNDEFINED
;
2980 if (scm_is_eq (m
, SCM_INUM0
))
2982 report_overflow
= 1;
2986 if (!coerce_to_big (n
, n_tmp
))
2988 value_of_wrong_type
= n
;
2989 position_of_wrong_type
= 1;
2993 if (!coerce_to_big (k
, k_tmp
))
2995 value_of_wrong_type
= k
;
2996 position_of_wrong_type
= 2;
3000 if (!coerce_to_big (m
, m_tmp
))
3002 value_of_wrong_type
= m
;
3003 position_of_wrong_type
= 3;
3007 /* if the exponent K is negative, and we simply call mpz_powm, we
3008 will get a divide-by-zero exception when an inverse 1/n mod m
3009 doesn't exist (or is not unique). Since exceptions are hard to
3010 handle, we'll attempt the inversion "by hand" -- that way, we get
3011 a simple failure code, which is easy to handle. */
3013 if (-1 == mpz_sgn (k_tmp
))
3015 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
3017 report_overflow
= 1;
3020 mpz_neg (k_tmp
, k_tmp
);
3023 result
= scm_i_mkbig ();
3024 mpz_powm (SCM_I_BIG_MPZ (result
),
3029 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
3030 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
3037 if (report_overflow
)
3038 scm_num_overflow (FUNC_NAME
);
3040 if (position_of_wrong_type
)
3041 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
3042 value_of_wrong_type
);
3044 return scm_i_normbig (result
);
3048 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
3050 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
3051 "exact integer, @var{n} can be any number.\n"
3053 "Negative @var{k} is supported, and results in\n"
3054 "@math{1/@var{n}^abs(@var{k})} in the usual way.\n"
3055 "@math{@var{n}^0} is 1, as usual, and that\n"
3056 "includes @math{0^0} is 1.\n"
3059 "(integer-expt 2 5) @result{} 32\n"
3060 "(integer-expt -3 3) @result{} -27\n"
3061 "(integer-expt 5 -3) @result{} 1/125\n"
3062 "(integer-expt 0 0) @result{} 1\n"
3064 #define FUNC_NAME s_scm_integer_expt
3067 SCM z_i2
= SCM_BOOL_F
;
3069 SCM acc
= SCM_I_MAKINUM (1L);
3071 /* Specifically refrain from checking the type of the first argument.
3072 This allows us to exponentiate any object that can be multiplied.
3073 If we must raise to a negative power, we must also be able to
3074 take its reciprocal. */
3075 if (!SCM_LIKELY (SCM_I_INUMP (k
)) && !SCM_LIKELY (SCM_BIGP (k
)))
3076 SCM_WRONG_TYPE_ARG (2, k
);
3078 if (SCM_UNLIKELY (scm_is_eq (k
, SCM_INUM0
)))
3079 return SCM_INUM1
; /* n^(exact0) is exact 1, regardless of n */
3080 else if (SCM_UNLIKELY (scm_is_eq (n
, SCM_I_MAKINUM (-1L))))
3081 return scm_is_false (scm_even_p (k
)) ? n
: SCM_INUM1
;
3082 /* The next check is necessary only because R6RS specifies different
3083 behavior for 0^(-k) than for (/ 0). If n is not a scheme number,
3084 we simply skip this case and move on. */
3085 else if (SCM_NUMBERP (n
) && scm_is_true (scm_zero_p (n
)))
3087 /* k cannot be 0 at this point, because we
3088 have already checked for that case above */
3089 if (scm_is_true (scm_positive_p (k
)))
3091 else /* return NaN for (0 ^ k) for negative k per R6RS */
3095 if (SCM_I_INUMP (k
))
3096 i2
= SCM_I_INUM (k
);
3097 else if (SCM_BIGP (k
))
3099 z_i2
= scm_i_clonebig (k
, 1);
3100 scm_remember_upto_here_1 (k
);
3104 SCM_WRONG_TYPE_ARG (2, k
);
3108 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
3110 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
3111 n
= scm_divide (n
, SCM_UNDEFINED
);
3115 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
3119 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
3121 return scm_product (acc
, n
);
3123 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
3124 acc
= scm_product (acc
, n
);
3125 n
= scm_product (n
, n
);
3126 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
3134 n
= scm_divide (n
, SCM_UNDEFINED
);
3141 return scm_product (acc
, n
);
3143 acc
= scm_product (acc
, n
);
3144 n
= scm_product (n
, n
);
3151 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
3153 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
3154 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
3156 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
3157 "@var{cnt} is negative it's a division, rounded towards negative\n"
3158 "infinity. (Note that this is not the same rounding as\n"
3159 "@code{quotient} does.)\n"
3161 "With @var{n} viewed as an infinite precision twos complement,\n"
3162 "@code{ash} means a left shift introducing zero bits, or a right\n"
3163 "shift dropping bits.\n"
3166 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
3167 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
3169 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
3170 "(ash -23 -2) @result{} -6\n"
3172 #define FUNC_NAME s_scm_ash
3175 bits_to_shift
= scm_to_long (cnt
);
3177 if (SCM_I_INUMP (n
))
3179 scm_t_inum nn
= SCM_I_INUM (n
);
3181 if (bits_to_shift
> 0)
3183 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
3184 overflow a non-zero fixnum. For smaller shifts we check the
3185 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
3186 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
3187 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
3193 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
3195 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
3198 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
3202 SCM result
= scm_i_inum2big (nn
);
3203 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
3210 bits_to_shift
= -bits_to_shift
;
3211 if (bits_to_shift
>= SCM_LONG_BIT
)
3212 return (nn
>= 0 ? SCM_INUM0
: SCM_I_MAKINUM(-1));
3214 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
3218 else if (SCM_BIGP (n
))
3222 if (bits_to_shift
== 0)
3225 result
= scm_i_mkbig ();
3226 if (bits_to_shift
>= 0)
3228 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
3234 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
3235 we have to allocate a bignum even if the result is going to be a
3237 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
3239 return scm_i_normbig (result
);
3245 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3251 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
3252 (SCM n
, SCM start
, SCM end
),
3253 "Return the integer composed of the @var{start} (inclusive)\n"
3254 "through @var{end} (exclusive) bits of @var{n}. The\n"
3255 "@var{start}th bit becomes the 0-th bit in the result.\n"
3258 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
3259 " @result{} \"1010\"\n"
3260 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
3261 " @result{} \"10110\"\n"
3263 #define FUNC_NAME s_scm_bit_extract
3265 unsigned long int istart
, iend
, bits
;
3266 istart
= scm_to_ulong (start
);
3267 iend
= scm_to_ulong (end
);
3268 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
3270 /* how many bits to keep */
3271 bits
= iend
- istart
;
3273 if (SCM_I_INUMP (n
))
3275 scm_t_inum in
= SCM_I_INUM (n
);
3277 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
3278 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
3279 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
3281 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
3283 /* Since we emulate two's complement encoded numbers, this
3284 * special case requires us to produce a result that has
3285 * more bits than can be stored in a fixnum.
3287 SCM result
= scm_i_inum2big (in
);
3288 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
3293 /* mask down to requisite bits */
3294 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
3295 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
3297 else if (SCM_BIGP (n
))
3302 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
3306 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
3307 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
3308 such bits into a ulong. */
3309 result
= scm_i_mkbig ();
3310 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
3311 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
3312 result
= scm_i_normbig (result
);
3314 scm_remember_upto_here_1 (n
);
3318 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3323 static const char scm_logtab
[] = {
3324 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
3327 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
3329 "Return the number of bits in integer @var{n}. If integer is\n"
3330 "positive, the 1-bits in its binary representation are counted.\n"
3331 "If negative, the 0-bits in its two's-complement binary\n"
3332 "representation are counted. If 0, 0 is returned.\n"
3335 "(logcount #b10101010)\n"
3342 #define FUNC_NAME s_scm_logcount
3344 if (SCM_I_INUMP (n
))
3346 unsigned long c
= 0;
3347 scm_t_inum nn
= SCM_I_INUM (n
);
3352 c
+= scm_logtab
[15 & nn
];
3355 return SCM_I_MAKINUM (c
);
3357 else if (SCM_BIGP (n
))
3359 unsigned long count
;
3360 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
3361 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
3363 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
3364 scm_remember_upto_here_1 (n
);
3365 return SCM_I_MAKINUM (count
);
3368 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3373 static const char scm_ilentab
[] = {
3374 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
3378 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
3380 "Return the number of bits necessary to represent @var{n}.\n"
3383 "(integer-length #b10101010)\n"
3385 "(integer-length 0)\n"
3387 "(integer-length #b1111)\n"
3390 #define FUNC_NAME s_scm_integer_length
3392 if (SCM_I_INUMP (n
))
3394 unsigned long c
= 0;
3396 scm_t_inum nn
= SCM_I_INUM (n
);
3402 l
= scm_ilentab
[15 & nn
];
3405 return SCM_I_MAKINUM (c
- 4 + l
);
3407 else if (SCM_BIGP (n
))
3409 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
3410 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
3411 1 too big, so check for that and adjust. */
3412 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
3413 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
3414 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
3415 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
3417 scm_remember_upto_here_1 (n
);
3418 return SCM_I_MAKINUM (size
);
3421 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3425 /*** NUMBERS -> STRINGS ***/
3426 #define SCM_MAX_DBL_PREC 60
3427 #define SCM_MAX_DBL_RADIX 36
3429 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
3430 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
3431 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
3434 void init_dblprec(int *prec
, int radix
) {
3435 /* determine floating point precision by adding successively
3436 smaller increments to 1.0 until it is considered == 1.0 */
3437 double f
= ((double)1.0)/radix
;
3438 double fsum
= 1.0 + f
;
3443 if (++(*prec
) > SCM_MAX_DBL_PREC
)
3455 void init_fx_radix(double *fx_list
, int radix
)
3457 /* initialize a per-radix list of tolerances. When added
3458 to a number < 1.0, we can determine if we should raund
3459 up and quit converting a number to a string. */
3463 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
3464 fx_list
[i
] = (fx_list
[i
-1] / radix
);
3467 /* use this array as a way to generate a single digit */
3468 static const char number_chars
[] = "0123456789abcdefghijklmnopqrstuvwxyz";
3471 idbl2str (double f
, char *a
, int radix
)
3473 int efmt
, dpt
, d
, i
, wp
;
3475 #ifdef DBL_MIN_10_EXP
3478 #endif /* DBL_MIN_10_EXP */
3483 radix
> SCM_MAX_DBL_RADIX
)
3485 /* revert to existing behavior */
3489 wp
= scm_dblprec
[radix
-2];
3490 fx
= fx_per_radix
[radix
-2];
3494 #ifdef HAVE_COPYSIGN
3495 double sgn
= copysign (1.0, f
);
3500 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
3506 strcpy (a
, "-inf.0");
3508 strcpy (a
, "+inf.0");
3513 strcpy (a
, "+nan.0");
3523 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
3524 make-uniform-vector, from causing infinite loops. */
3525 /* just do the checking...if it passes, we do the conversion for our
3526 radix again below */
3533 if (exp_cpy
-- < DBL_MIN_10_EXP
)
3541 while (f_cpy
> 10.0)
3544 if (exp_cpy
++ > DBL_MAX_10_EXP
)
3565 if (f
+ fx
[wp
] >= radix
)
3572 /* adding 9999 makes this equivalent to abs(x) % 3 */
3573 dpt
= (exp
+ 9999) % 3;
3577 efmt
= (exp
< -3) || (exp
> wp
+ 2);
3599 a
[ch
++] = number_chars
[d
];
3602 if (f
+ fx
[wp
] >= 1.0)
3604 a
[ch
- 1] = number_chars
[d
+1];
3616 if ((dpt
> 4) && (exp
> 6))
3618 d
= (a
[0] == '-' ? 2 : 1);
3619 for (i
= ch
++; i
> d
; i
--)
3632 if (a
[ch
- 1] == '.')
3633 a
[ch
++] = '0'; /* trailing zero */
3642 for (i
= radix
; i
<= exp
; i
*= radix
);
3643 for (i
/= radix
; i
; i
/= radix
)
3645 a
[ch
++] = number_chars
[exp
/ i
];
3654 icmplx2str (double real
, double imag
, char *str
, int radix
)
3659 i
= idbl2str (real
, str
, radix
);
3660 #ifdef HAVE_COPYSIGN
3661 sgn
= copysign (1.0, imag
);
3665 /* Don't output a '+' for negative numbers or for Inf and
3666 NaN. They will provide their own sign. */
3667 if (sgn
>= 0 && DOUBLE_IS_FINITE (imag
))
3669 i
+= idbl2str (imag
, &str
[i
], radix
);
3675 iflo2str (SCM flt
, char *str
, int radix
)
3678 if (SCM_REALP (flt
))
3679 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
3681 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
3686 /* convert a scm_t_intmax to a string (unterminated). returns the number of
3687 characters in the result.
3689 p is destination: worst case (base 2) is SCM_INTBUFLEN */
3691 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
3696 return scm_iuint2str (-num
, rad
, p
) + 1;
3699 return scm_iuint2str (num
, rad
, p
);
3702 /* convert a scm_t_intmax to a string (unterminated). returns the number of
3703 characters in the result.
3705 p is destination: worst case (base 2) is SCM_INTBUFLEN */
3707 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
3711 scm_t_uintmax n
= num
;
3713 if (rad
< 2 || rad
> 36)
3714 scm_out_of_range ("scm_iuint2str", scm_from_int (rad
));
3716 for (n
/= rad
; n
> 0; n
/= rad
)
3726 p
[i
] = number_chars
[d
];
3731 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
3733 "Return a string holding the external representation of the\n"
3734 "number @var{n} in the given @var{radix}. If @var{n} is\n"
3735 "inexact, a radix of 10 will be used.")
3736 #define FUNC_NAME s_scm_number_to_string
3740 if (SCM_UNBNDP (radix
))
3743 base
= scm_to_signed_integer (radix
, 2, 36);
3745 if (SCM_I_INUMP (n
))
3747 char num_buf
[SCM_INTBUFLEN
];
3748 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
3749 return scm_from_locale_stringn (num_buf
, length
);
3751 else if (SCM_BIGP (n
))
3753 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
3754 scm_remember_upto_here_1 (n
);
3755 return scm_take_locale_string (str
);
3757 else if (SCM_FRACTIONP (n
))
3759 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
3760 scm_from_locale_string ("/"),
3761 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
3763 else if (SCM_INEXACTP (n
))
3765 char num_buf
[FLOBUFLEN
];
3766 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
3769 SCM_WRONG_TYPE_ARG (1, n
);
3774 /* These print routines used to be stubbed here so that scm_repl.c
3775 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
3778 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3780 char num_buf
[FLOBUFLEN
];
3781 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
3786 scm_i_print_double (double val
, SCM port
)
3788 char num_buf
[FLOBUFLEN
];
3789 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
3793 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3796 char num_buf
[FLOBUFLEN
];
3797 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
3802 scm_i_print_complex (double real
, double imag
, SCM port
)
3804 char num_buf
[FLOBUFLEN
];
3805 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
3809 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3812 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
3813 scm_display (str
, port
);
3814 scm_remember_upto_here_1 (str
);
3819 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3821 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
3822 scm_remember_upto_here_1 (exp
);
3823 scm_lfwrite (str
, (size_t) strlen (str
), port
);
3827 /*** END nums->strs ***/
3830 /*** STRINGS -> NUMBERS ***/
3832 /* The following functions implement the conversion from strings to numbers.
3833 * The implementation somehow follows the grammar for numbers as it is given
3834 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
3835 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
3836 * points should be noted about the implementation:
3837 * * Each function keeps a local index variable 'idx' that points at the
3838 * current position within the parsed string. The global index is only
3839 * updated if the function could parse the corresponding syntactic unit
3841 * * Similarly, the functions keep track of indicators of inexactness ('#',
3842 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
3843 * global exactness information is only updated after each part has been
3844 * successfully parsed.
3845 * * Sequences of digits are parsed into temporary variables holding fixnums.
3846 * Only if these fixnums would overflow, the result variables are updated
3847 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
3848 * the temporary variables holding the fixnums are cleared, and the process
3849 * starts over again. If for example fixnums were able to store five decimal
3850 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
3851 * and the result was computed as 12345 * 100000 + 67890. In other words,
3852 * only every five digits two bignum operations were performed.
3855 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
3857 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
3859 /* Caller is responsible for checking that the return value is in range
3860 for the given radix, which should be <= 36. */
3862 char_decimal_value (scm_t_uint32 c
)
3864 /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
3865 that's certainly above any valid decimal, so we take advantage of
3866 that to elide some tests. */
3867 unsigned int d
= (unsigned int) uc_decimal_value (c
);
3869 /* If that failed, try extended hexadecimals, then. Only accept ascii
3874 if (c
>= (scm_t_uint32
) 'a')
3875 d
= c
- (scm_t_uint32
)'a' + 10U;
3881 mem2uinteger (SCM mem
, unsigned int *p_idx
,
3882 unsigned int radix
, enum t_exactness
*p_exactness
)
3884 unsigned int idx
= *p_idx
;
3885 unsigned int hash_seen
= 0;
3886 scm_t_bits shift
= 1;
3888 unsigned int digit_value
;
3891 size_t len
= scm_i_string_length (mem
);
3896 c
= scm_i_string_ref (mem
, idx
);
3897 digit_value
= char_decimal_value (c
);
3898 if (digit_value
>= radix
)
3902 result
= SCM_I_MAKINUM (digit_value
);
3905 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
3915 digit_value
= char_decimal_value (c
);
3916 /* This check catches non-decimals in addition to out-of-range
3918 if (digit_value
>= radix
)
3923 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
3925 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3927 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3934 shift
= shift
* radix
;
3935 add
= add
* radix
+ digit_value
;
3940 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3942 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3946 *p_exactness
= INEXACT
;
3952 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
3953 * covers the parts of the rules that start at a potential point. The value
3954 * of the digits up to the point have been parsed by the caller and are given
3955 * in variable result. The content of *p_exactness indicates, whether a hash
3956 * has already been seen in the digits before the point.
3959 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
3962 mem2decimal_from_point (SCM result
, SCM mem
,
3963 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
3965 unsigned int idx
= *p_idx
;
3966 enum t_exactness x
= *p_exactness
;
3967 size_t len
= scm_i_string_length (mem
);
3972 if (scm_i_string_ref (mem
, idx
) == '.')
3974 scm_t_bits shift
= 1;
3976 unsigned int digit_value
;
3977 SCM big_shift
= SCM_INUM1
;
3982 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
3983 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
3988 digit_value
= DIGIT2UINT (c
);
3999 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
4001 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
4002 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
4004 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
4012 add
= add
* 10 + digit_value
;
4018 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
4019 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
4020 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
4023 result
= scm_divide (result
, big_shift
);
4025 /* We've seen a decimal point, thus the value is implicitly inexact. */
4037 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
4039 switch (scm_i_string_ref (mem
, idx
))
4051 c
= scm_i_string_ref (mem
, idx
);
4059 c
= scm_i_string_ref (mem
, idx
);
4068 c
= scm_i_string_ref (mem
, idx
);
4073 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4077 exponent
= DIGIT2UINT (c
);
4080 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
4081 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4084 if (exponent
<= SCM_MAXEXP
)
4085 exponent
= exponent
* 10 + DIGIT2UINT (c
);
4091 if (exponent
> SCM_MAXEXP
)
4093 size_t exp_len
= idx
- start
;
4094 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
4095 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
4096 scm_out_of_range ("string->number", exp_num
);
4099 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
4101 result
= scm_product (result
, e
);
4103 result
= scm_divide2real (result
, e
);
4105 /* We've seen an exponent, thus the value is implicitly inexact. */
4123 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
4126 mem2ureal (SCM mem
, unsigned int *p_idx
,
4127 unsigned int radix
, enum t_exactness forced_x
)
4129 unsigned int idx
= *p_idx
;
4131 size_t len
= scm_i_string_length (mem
);
4133 /* Start off believing that the number will be exact. This changes
4134 to INEXACT if we see a decimal point or a hash. */
4135 enum t_exactness implicit_x
= EXACT
;
4140 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
4146 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
4148 /* Cobble up the fractional part. We might want to set the
4149 NaN's mantissa from it. */
4151 mem2uinteger (mem
, &idx
, 10, &implicit_x
);
4156 if (scm_i_string_ref (mem
, idx
) == '.')
4160 else if (idx
+ 1 == len
)
4162 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
4165 result
= mem2decimal_from_point (SCM_INUM0
, mem
,
4166 p_idx
, &implicit_x
);
4172 uinteger
= mem2uinteger (mem
, &idx
, radix
, &implicit_x
);
4173 if (scm_is_false (uinteger
))
4178 else if (scm_i_string_ref (mem
, idx
) == '/')
4186 divisor
= mem2uinteger (mem
, &idx
, radix
, &implicit_x
);
4187 if (scm_is_false (divisor
))
4190 /* both are int/big here, I assume */
4191 result
= scm_i_make_ratio (uinteger
, divisor
);
4193 else if (radix
== 10)
4195 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &implicit_x
);
4196 if (scm_is_false (result
))
4208 if (SCM_INEXACTP (result
))
4209 return scm_inexact_to_exact (result
);
4213 if (SCM_INEXACTP (result
))
4216 return scm_exact_to_inexact (result
);
4218 if (implicit_x
== INEXACT
)
4220 if (SCM_INEXACTP (result
))
4223 return scm_exact_to_inexact (result
);
4229 /* We should never get here */
4230 scm_syserror ("mem2ureal");
4234 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
4237 mem2complex (SCM mem
, unsigned int idx
,
4238 unsigned int radix
, enum t_exactness forced_x
)
4243 size_t len
= scm_i_string_length (mem
);
4248 c
= scm_i_string_ref (mem
, idx
);
4263 ureal
= mem2ureal (mem
, &idx
, radix
, forced_x
);
4264 if (scm_is_false (ureal
))
4266 /* input must be either +i or -i */
4271 if (scm_i_string_ref (mem
, idx
) == 'i'
4272 || scm_i_string_ref (mem
, idx
) == 'I')
4278 return scm_make_rectangular (SCM_INUM0
, SCM_I_MAKINUM (sign
));
4285 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
4286 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
4291 c
= scm_i_string_ref (mem
, idx
);
4295 /* either +<ureal>i or -<ureal>i */
4302 return scm_make_rectangular (SCM_INUM0
, ureal
);
4305 /* polar input: <real>@<real>. */
4316 c
= scm_i_string_ref (mem
, idx
);
4334 angle
= mem2ureal (mem
, &idx
, radix
, forced_x
);
4335 if (scm_is_false (angle
))
4340 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
4341 angle
= scm_difference (angle
, SCM_UNDEFINED
);
4343 result
= scm_make_polar (ureal
, angle
);
4348 /* expecting input matching <real>[+-]<ureal>?i */
4355 int sign
= (c
== '+') ? 1 : -1;
4356 SCM imag
= mem2ureal (mem
, &idx
, radix
, forced_x
);
4358 if (scm_is_false (imag
))
4359 imag
= SCM_I_MAKINUM (sign
);
4360 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
4361 imag
= scm_difference (imag
, SCM_UNDEFINED
);
4365 if (scm_i_string_ref (mem
, idx
) != 'i'
4366 && scm_i_string_ref (mem
, idx
) != 'I')
4373 return scm_make_rectangular (ureal
, imag
);
4382 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
4384 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
4387 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
4389 unsigned int idx
= 0;
4390 unsigned int radix
= NO_RADIX
;
4391 enum t_exactness forced_x
= NO_EXACTNESS
;
4392 size_t len
= scm_i_string_length (mem
);
4394 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
4395 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
4397 switch (scm_i_string_ref (mem
, idx
+ 1))
4400 if (radix
!= NO_RADIX
)
4405 if (radix
!= NO_RADIX
)
4410 if (forced_x
!= NO_EXACTNESS
)
4415 if (forced_x
!= NO_EXACTNESS
)
4420 if (radix
!= NO_RADIX
)
4425 if (radix
!= NO_RADIX
)
4435 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
4436 if (radix
== NO_RADIX
)
4437 radix
= default_radix
;
4439 return mem2complex (mem
, idx
, radix
, forced_x
);
4443 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
4444 unsigned int default_radix
)
4446 SCM str
= scm_from_locale_stringn (mem
, len
);
4448 return scm_i_string_to_number (str
, default_radix
);
4452 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
4453 (SCM string
, SCM radix
),
4454 "Return a number of the maximally precise representation\n"
4455 "expressed by the given @var{string}. @var{radix} must be an\n"
4456 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
4457 "is a default radix that may be overridden by an explicit radix\n"
4458 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
4459 "supplied, then the default radix is 10. If string is not a\n"
4460 "syntactically valid notation for a number, then\n"
4461 "@code{string->number} returns @code{#f}.")
4462 #define FUNC_NAME s_scm_string_to_number
4466 SCM_VALIDATE_STRING (1, string
);
4468 if (SCM_UNBNDP (radix
))
4471 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
4473 answer
= scm_i_string_to_number (string
, base
);
4474 scm_remember_upto_here_1 (string
);
4480 /*** END strs->nums ***/
4483 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
4485 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
4487 #define FUNC_NAME s_scm_number_p
4489 return scm_from_bool (SCM_NUMBERP (x
));
4493 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
4495 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
4496 "otherwise. Note that the sets of real, rational and integer\n"
4497 "values form subsets of the set of complex numbers, i. e. the\n"
4498 "predicate will also be fulfilled if @var{x} is a real,\n"
4499 "rational or integer number.")
4500 #define FUNC_NAME s_scm_complex_p
4502 /* all numbers are complex. */
4503 return scm_number_p (x
);
4507 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
4509 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
4510 "otherwise. Note that the set of integer values forms a subset of\n"
4511 "the set of real numbers, i. e. the predicate will also be\n"
4512 "fulfilled if @var{x} is an integer number.")
4513 #define FUNC_NAME s_scm_real_p
4515 return scm_from_bool
4516 (SCM_I_INUMP (x
) || SCM_REALP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
));
4520 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
4522 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
4523 "otherwise. Note that the set of integer values forms a subset of\n"
4524 "the set of rational numbers, i. e. the predicate will also be\n"
4525 "fulfilled if @var{x} is an integer number.")
4526 #define FUNC_NAME s_scm_rational_p
4528 if (SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
))
4530 else if (SCM_REALP (x
))
4531 /* due to their limited precision, finite floating point numbers are
4532 rational as well. (finite means neither infinity nor a NaN) */
4533 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
4539 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
4541 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
4543 #define FUNC_NAME s_scm_integer_p
4545 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4547 else if (SCM_REALP (x
))
4549 double val
= SCM_REAL_VALUE (x
);
4550 return scm_from_bool (!isinf (val
) && (val
== floor (val
)));
4558 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
4559 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
4560 (SCM x
, SCM y
, SCM rest
),
4561 "Return @code{#t} if all parameters are numerically equal.")
4562 #define FUNC_NAME s_scm_i_num_eq_p
4564 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4566 while (!scm_is_null (rest
))
4568 if (scm_is_false (scm_num_eq_p (x
, y
)))
4572 rest
= scm_cdr (rest
);
4574 return scm_num_eq_p (x
, y
);
4578 scm_num_eq_p (SCM x
, SCM y
)
4581 if (SCM_I_INUMP (x
))
4583 scm_t_signed_bits xx
= SCM_I_INUM (x
);
4584 if (SCM_I_INUMP (y
))
4586 scm_t_signed_bits yy
= SCM_I_INUM (y
);
4587 return scm_from_bool (xx
== yy
);
4589 else if (SCM_BIGP (y
))
4591 else if (SCM_REALP (y
))
4593 /* On a 32-bit system an inum fits a double, we can cast the inum
4594 to a double and compare.
4596 But on a 64-bit system an inum is bigger than a double and
4597 casting it to a double (call that dxx) will round. dxx is at
4598 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
4599 an integer and fits a long. So we cast yy to a long and
4600 compare with plain xx.
4602 An alternative (for any size system actually) would be to check
4603 yy is an integer (with floor) and is in range of an inum
4604 (compare against appropriate powers of 2) then test
4605 xx==(scm_t_signed_bits)yy. It's just a matter of which
4606 casts/comparisons might be fastest or easiest for the cpu. */
4608 double yy
= SCM_REAL_VALUE (y
);
4609 return scm_from_bool ((double) xx
== yy
4610 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
4611 || xx
== (scm_t_signed_bits
) yy
));
4613 else if (SCM_COMPLEXP (y
))
4614 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
4615 && (0.0 == SCM_COMPLEX_IMAG (y
)));
4616 else if (SCM_FRACTIONP (y
))
4619 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4621 else if (SCM_BIGP (x
))
4623 if (SCM_I_INUMP (y
))
4625 else if (SCM_BIGP (y
))
4627 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4628 scm_remember_upto_here_2 (x
, y
);
4629 return scm_from_bool (0 == cmp
);
4631 else if (SCM_REALP (y
))
4634 if (isnan (SCM_REAL_VALUE (y
)))
4636 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
4637 scm_remember_upto_here_1 (x
);
4638 return scm_from_bool (0 == cmp
);
4640 else if (SCM_COMPLEXP (y
))
4643 if (0.0 != SCM_COMPLEX_IMAG (y
))
4645 if (isnan (SCM_COMPLEX_REAL (y
)))
4647 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
4648 scm_remember_upto_here_1 (x
);
4649 return scm_from_bool (0 == cmp
);
4651 else if (SCM_FRACTIONP (y
))
4654 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4656 else if (SCM_REALP (x
))
4658 double xx
= SCM_REAL_VALUE (x
);
4659 if (SCM_I_INUMP (y
))
4661 /* see comments with inum/real above */
4662 scm_t_signed_bits yy
= SCM_I_INUM (y
);
4663 return scm_from_bool (xx
== (double) yy
4664 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
4665 || (scm_t_signed_bits
) xx
== yy
));
4667 else if (SCM_BIGP (y
))
4670 if (isnan (SCM_REAL_VALUE (x
)))
4672 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
4673 scm_remember_upto_here_1 (y
);
4674 return scm_from_bool (0 == cmp
);
4676 else if (SCM_REALP (y
))
4677 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
4678 else if (SCM_COMPLEXP (y
))
4679 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
4680 && (0.0 == SCM_COMPLEX_IMAG (y
)));
4681 else if (SCM_FRACTIONP (y
))
4683 double xx
= SCM_REAL_VALUE (x
);
4687 return scm_from_bool (xx
< 0.0);
4688 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4692 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4694 else if (SCM_COMPLEXP (x
))
4696 if (SCM_I_INUMP (y
))
4697 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
4698 && (SCM_COMPLEX_IMAG (x
) == 0.0));
4699 else if (SCM_BIGP (y
))
4702 if (0.0 != SCM_COMPLEX_IMAG (x
))
4704 if (isnan (SCM_COMPLEX_REAL (x
)))
4706 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
4707 scm_remember_upto_here_1 (y
);
4708 return scm_from_bool (0 == cmp
);
4710 else if (SCM_REALP (y
))
4711 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
4712 && (SCM_COMPLEX_IMAG (x
) == 0.0));
4713 else if (SCM_COMPLEXP (y
))
4714 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
4715 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
4716 else if (SCM_FRACTIONP (y
))
4719 if (SCM_COMPLEX_IMAG (x
) != 0.0)
4721 xx
= SCM_COMPLEX_REAL (x
);
4725 return scm_from_bool (xx
< 0.0);
4726 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4730 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4732 else if (SCM_FRACTIONP (x
))
4734 if (SCM_I_INUMP (y
))
4736 else if (SCM_BIGP (y
))
4738 else if (SCM_REALP (y
))
4740 double yy
= SCM_REAL_VALUE (y
);
4744 return scm_from_bool (0.0 < yy
);
4745 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4748 else if (SCM_COMPLEXP (y
))
4751 if (SCM_COMPLEX_IMAG (y
) != 0.0)
4753 yy
= SCM_COMPLEX_REAL (y
);
4757 return scm_from_bool (0.0 < yy
);
4758 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4761 else if (SCM_FRACTIONP (y
))
4762 return scm_i_fraction_equalp (x
, y
);
4764 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4767 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
4771 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
4772 done are good for inums, but for bignums an answer can almost always be
4773 had by just examining a few high bits of the operands, as done by GMP in
4774 mpq_cmp. flonum/frac compares likewise, but with the slight complication
4775 of the float exponent to take into account. */
4777 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
4778 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
4779 (SCM x
, SCM y
, SCM rest
),
4780 "Return @code{#t} if the list of parameters is monotonically\n"
4782 #define FUNC_NAME s_scm_i_num_less_p
4784 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4786 while (!scm_is_null (rest
))
4788 if (scm_is_false (scm_less_p (x
, y
)))
4792 rest
= scm_cdr (rest
);
4794 return scm_less_p (x
, y
);
4798 scm_less_p (SCM x
, SCM y
)
4801 if (SCM_I_INUMP (x
))
4803 scm_t_inum xx
= SCM_I_INUM (x
);
4804 if (SCM_I_INUMP (y
))
4806 scm_t_inum yy
= SCM_I_INUM (y
);
4807 return scm_from_bool (xx
< yy
);
4809 else if (SCM_BIGP (y
))
4811 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4812 scm_remember_upto_here_1 (y
);
4813 return scm_from_bool (sgn
> 0);
4815 else if (SCM_REALP (y
))
4816 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
4817 else if (SCM_FRACTIONP (y
))
4819 /* "x < a/b" becomes "x*b < a" */
4821 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
4822 y
= SCM_FRACTION_NUMERATOR (y
);
4826 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4828 else if (SCM_BIGP (x
))
4830 if (SCM_I_INUMP (y
))
4832 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4833 scm_remember_upto_here_1 (x
);
4834 return scm_from_bool (sgn
< 0);
4836 else if (SCM_BIGP (y
))
4838 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4839 scm_remember_upto_here_2 (x
, y
);
4840 return scm_from_bool (cmp
< 0);
4842 else if (SCM_REALP (y
))
4845 if (isnan (SCM_REAL_VALUE (y
)))
4847 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
4848 scm_remember_upto_here_1 (x
);
4849 return scm_from_bool (cmp
< 0);
4851 else if (SCM_FRACTIONP (y
))
4854 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4856 else if (SCM_REALP (x
))
4858 if (SCM_I_INUMP (y
))
4859 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
4860 else if (SCM_BIGP (y
))
4863 if (isnan (SCM_REAL_VALUE (x
)))
4865 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
4866 scm_remember_upto_here_1 (y
);
4867 return scm_from_bool (cmp
> 0);
4869 else if (SCM_REALP (y
))
4870 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
4871 else if (SCM_FRACTIONP (y
))
4873 double xx
= SCM_REAL_VALUE (x
);
4877 return scm_from_bool (xx
< 0.0);
4878 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4882 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4884 else if (SCM_FRACTIONP (x
))
4886 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
4888 /* "a/b < y" becomes "a < y*b" */
4889 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
4890 x
= SCM_FRACTION_NUMERATOR (x
);
4893 else if (SCM_REALP (y
))
4895 double yy
= SCM_REAL_VALUE (y
);
4899 return scm_from_bool (0.0 < yy
);
4900 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4903 else if (SCM_FRACTIONP (y
))
4905 /* "a/b < c/d" becomes "a*d < c*b" */
4906 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
4907 SCM_FRACTION_DENOMINATOR (y
));
4908 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
4909 SCM_FRACTION_DENOMINATOR (x
));
4915 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4918 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
4922 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
4923 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
4924 (SCM x
, SCM y
, SCM rest
),
4925 "Return @code{#t} if the list of parameters is monotonically\n"
4927 #define FUNC_NAME s_scm_i_num_gr_p
4929 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4931 while (!scm_is_null (rest
))
4933 if (scm_is_false (scm_gr_p (x
, y
)))
4937 rest
= scm_cdr (rest
);
4939 return scm_gr_p (x
, y
);
4942 #define FUNC_NAME s_scm_i_num_gr_p
4944 scm_gr_p (SCM x
, SCM y
)
4946 if (!SCM_NUMBERP (x
))
4947 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
4948 else if (!SCM_NUMBERP (y
))
4949 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
4951 return scm_less_p (y
, x
);
4956 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
4957 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
4958 (SCM x
, SCM y
, SCM rest
),
4959 "Return @code{#t} if the list of parameters is monotonically\n"
4961 #define FUNC_NAME s_scm_i_num_leq_p
4963 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4965 while (!scm_is_null (rest
))
4967 if (scm_is_false (scm_leq_p (x
, y
)))
4971 rest
= scm_cdr (rest
);
4973 return scm_leq_p (x
, y
);
4976 #define FUNC_NAME s_scm_i_num_leq_p
4978 scm_leq_p (SCM x
, SCM y
)
4980 if (!SCM_NUMBERP (x
))
4981 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
4982 else if (!SCM_NUMBERP (y
))
4983 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
4984 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
4987 return scm_not (scm_less_p (y
, x
));
4992 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
4993 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
4994 (SCM x
, SCM y
, SCM rest
),
4995 "Return @code{#t} if the list of parameters is monotonically\n"
4997 #define FUNC_NAME s_scm_i_num_geq_p
4999 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
5001 while (!scm_is_null (rest
))
5003 if (scm_is_false (scm_geq_p (x
, y
)))
5007 rest
= scm_cdr (rest
);
5009 return scm_geq_p (x
, y
);
5012 #define FUNC_NAME s_scm_i_num_geq_p
5014 scm_geq_p (SCM x
, SCM y
)
5016 if (!SCM_NUMBERP (x
))
5017 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
5018 else if (!SCM_NUMBERP (y
))
5019 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
5020 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
5023 return scm_not (scm_less_p (x
, y
));
5028 SCM_PRIMITIVE_GENERIC (scm_zero_p
, "zero?", 1, 0, 0,
5030 "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
5032 #define FUNC_NAME s_scm_zero_p
5034 if (SCM_I_INUMP (z
))
5035 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
5036 else if (SCM_BIGP (z
))
5038 else if (SCM_REALP (z
))
5039 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
5040 else if (SCM_COMPLEXP (z
))
5041 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
5042 && SCM_COMPLEX_IMAG (z
) == 0.0);
5043 else if (SCM_FRACTIONP (z
))
5046 SCM_WTA_DISPATCH_1 (g_scm_zero_p
, z
, SCM_ARG1
, s_scm_zero_p
);
5051 SCM_PRIMITIVE_GENERIC (scm_positive_p
, "positive?", 1, 0, 0,
5053 "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
5055 #define FUNC_NAME s_scm_positive_p
5057 if (SCM_I_INUMP (x
))
5058 return scm_from_bool (SCM_I_INUM (x
) > 0);
5059 else if (SCM_BIGP (x
))
5061 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5062 scm_remember_upto_here_1 (x
);
5063 return scm_from_bool (sgn
> 0);
5065 else if (SCM_REALP (x
))
5066 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
5067 else if (SCM_FRACTIONP (x
))
5068 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
5070 SCM_WTA_DISPATCH_1 (g_scm_positive_p
, x
, SCM_ARG1
, s_scm_positive_p
);
5075 SCM_PRIMITIVE_GENERIC (scm_negative_p
, "negative?", 1, 0, 0,
5077 "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
5079 #define FUNC_NAME s_scm_negative_p
5081 if (SCM_I_INUMP (x
))
5082 return scm_from_bool (SCM_I_INUM (x
) < 0);
5083 else if (SCM_BIGP (x
))
5085 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5086 scm_remember_upto_here_1 (x
);
5087 return scm_from_bool (sgn
< 0);
5089 else if (SCM_REALP (x
))
5090 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
5091 else if (SCM_FRACTIONP (x
))
5092 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
5094 SCM_WTA_DISPATCH_1 (g_scm_negative_p
, x
, SCM_ARG1
, s_scm_negative_p
);
5099 /* scm_min and scm_max return an inexact when either argument is inexact, as
5100 required by r5rs. On that basis, for exact/inexact combinations the
5101 exact is converted to inexact to compare and possibly return. This is
5102 unlike scm_less_p above which takes some trouble to preserve all bits in
5103 its test, such trouble is not required for min and max. */
5105 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
5106 (SCM x
, SCM y
, SCM rest
),
5107 "Return the maximum of all parameter values.")
5108 #define FUNC_NAME s_scm_i_max
5110 while (!scm_is_null (rest
))
5111 { x
= scm_max (x
, y
);
5113 rest
= scm_cdr (rest
);
5115 return scm_max (x
, y
);
5119 #define s_max s_scm_i_max
5120 #define g_max g_scm_i_max
5123 scm_max (SCM x
, SCM y
)
5128 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
5129 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
5132 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
5135 if (SCM_I_INUMP (x
))
5137 scm_t_inum xx
= SCM_I_INUM (x
);
5138 if (SCM_I_INUMP (y
))
5140 scm_t_inum yy
= SCM_I_INUM (y
);
5141 return (xx
< yy
) ? y
: x
;
5143 else if (SCM_BIGP (y
))
5145 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5146 scm_remember_upto_here_1 (y
);
5147 return (sgn
< 0) ? x
: y
;
5149 else if (SCM_REALP (y
))
5152 double yyd
= SCM_REAL_VALUE (y
);
5155 return scm_from_double (xxd
);
5156 /* If y is a NaN, then "==" is false and we return the NaN */
5157 else if (SCM_LIKELY (!(xxd
== yyd
)))
5159 /* Handle signed zeroes properly */
5165 else if (SCM_FRACTIONP (y
))
5168 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
5171 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5173 else if (SCM_BIGP (x
))
5175 if (SCM_I_INUMP (y
))
5177 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5178 scm_remember_upto_here_1 (x
);
5179 return (sgn
< 0) ? y
: x
;
5181 else if (SCM_BIGP (y
))
5183 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
5184 scm_remember_upto_here_2 (x
, y
);
5185 return (cmp
> 0) ? x
: y
;
5187 else if (SCM_REALP (y
))
5189 /* if y==NaN then xx>yy is false, so we return the NaN y */
5192 xx
= scm_i_big2dbl (x
);
5193 yy
= SCM_REAL_VALUE (y
);
5194 return (xx
> yy
? scm_from_double (xx
) : y
);
5196 else if (SCM_FRACTIONP (y
))
5201 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5203 else if (SCM_REALP (x
))
5205 if (SCM_I_INUMP (y
))
5207 scm_t_inum yy
= SCM_I_INUM (y
);
5208 double xxd
= SCM_REAL_VALUE (x
);
5212 return scm_from_double (yyd
);
5213 /* If x is a NaN, then "==" is false and we return the NaN */
5214 else if (SCM_LIKELY (!(xxd
== yyd
)))
5216 /* Handle signed zeroes properly */
5222 else if (SCM_BIGP (y
))
5227 else if (SCM_REALP (y
))
5229 double xx
= SCM_REAL_VALUE (x
);
5230 double yy
= SCM_REAL_VALUE (y
);
5232 /* For purposes of max: +inf.0 > nan > everything else, per R6RS */
5235 else if (SCM_LIKELY (xx
< yy
))
5237 /* If neither (xx > yy) nor (xx < yy), then
5238 either they're equal or one is a NaN */
5239 else if (SCM_UNLIKELY (isnan (xx
)))
5240 return DOUBLE_IS_POSITIVE_INFINITY (yy
) ? y
: x
;
5241 else if (SCM_UNLIKELY (isnan (yy
)))
5242 return DOUBLE_IS_POSITIVE_INFINITY (xx
) ? x
: y
;
5243 /* xx == yy, but handle signed zeroes properly */
5244 else if (double_is_non_negative_zero (yy
))
5249 else if (SCM_FRACTIONP (y
))
5251 double yy
= scm_i_fraction2double (y
);
5252 double xx
= SCM_REAL_VALUE (x
);
5253 return (xx
< yy
) ? scm_from_double (yy
) : x
;
5256 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5258 else if (SCM_FRACTIONP (x
))
5260 if (SCM_I_INUMP (y
))
5264 else if (SCM_BIGP (y
))
5268 else if (SCM_REALP (y
))
5270 double xx
= scm_i_fraction2double (x
);
5271 /* if y==NaN then ">" is false, so we return the NaN y */
5272 return (xx
> SCM_REAL_VALUE (y
)) ? scm_from_double (xx
) : y
;
5274 else if (SCM_FRACTIONP (y
))
5279 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5282 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
5286 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
5287 (SCM x
, SCM y
, SCM rest
),
5288 "Return the minimum of all parameter values.")
5289 #define FUNC_NAME s_scm_i_min
5291 while (!scm_is_null (rest
))
5292 { x
= scm_min (x
, y
);
5294 rest
= scm_cdr (rest
);
5296 return scm_min (x
, y
);
5300 #define s_min s_scm_i_min
5301 #define g_min g_scm_i_min
5304 scm_min (SCM x
, SCM y
)
5309 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
5310 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
5313 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
5316 if (SCM_I_INUMP (x
))
5318 scm_t_inum xx
= SCM_I_INUM (x
);
5319 if (SCM_I_INUMP (y
))
5321 scm_t_inum yy
= SCM_I_INUM (y
);
5322 return (xx
< yy
) ? x
: y
;
5324 else if (SCM_BIGP (y
))
5326 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5327 scm_remember_upto_here_1 (y
);
5328 return (sgn
< 0) ? y
: x
;
5330 else if (SCM_REALP (y
))
5333 /* if y==NaN then "<" is false and we return NaN */
5334 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
5336 else if (SCM_FRACTIONP (y
))
5339 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
5342 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5344 else if (SCM_BIGP (x
))
5346 if (SCM_I_INUMP (y
))
5348 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5349 scm_remember_upto_here_1 (x
);
5350 return (sgn
< 0) ? x
: y
;
5352 else if (SCM_BIGP (y
))
5354 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
5355 scm_remember_upto_here_2 (x
, y
);
5356 return (cmp
> 0) ? y
: x
;
5358 else if (SCM_REALP (y
))
5360 /* if y==NaN then xx<yy is false, so we return the NaN y */
5363 xx
= scm_i_big2dbl (x
);
5364 yy
= SCM_REAL_VALUE (y
);
5365 return (xx
< yy
? scm_from_double (xx
) : y
);
5367 else if (SCM_FRACTIONP (y
))
5372 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5374 else if (SCM_REALP (x
))
5376 if (SCM_I_INUMP (y
))
5378 double z
= SCM_I_INUM (y
);
5379 /* if x==NaN then "<" is false and we return NaN */
5380 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
5382 else if (SCM_BIGP (y
))
5387 else if (SCM_REALP (y
))
5389 double xx
= SCM_REAL_VALUE (x
);
5390 double yy
= SCM_REAL_VALUE (y
);
5392 /* For purposes of min: -inf.0 < nan < everything else, per R6RS */
5395 else if (SCM_LIKELY (xx
> yy
))
5397 /* If neither (xx < yy) nor (xx > yy), then
5398 either they're equal or one is a NaN */
5399 else if (SCM_UNLIKELY (isnan (xx
)))
5400 return DOUBLE_IS_NEGATIVE_INFINITY (yy
) ? y
: x
;
5401 else if (SCM_UNLIKELY (isnan (yy
)))
5402 return DOUBLE_IS_NEGATIVE_INFINITY (xx
) ? x
: y
;
5403 /* xx == yy, but handle signed zeroes properly */
5404 else if (double_is_non_negative_zero (xx
))
5409 else if (SCM_FRACTIONP (y
))
5411 double yy
= scm_i_fraction2double (y
);
5412 double xx
= SCM_REAL_VALUE (x
);
5413 return (yy
< xx
) ? scm_from_double (yy
) : x
;
5416 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5418 else if (SCM_FRACTIONP (x
))
5420 if (SCM_I_INUMP (y
))
5424 else if (SCM_BIGP (y
))
5428 else if (SCM_REALP (y
))
5430 double xx
= scm_i_fraction2double (x
);
5431 /* if y==NaN then "<" is false, so we return the NaN y */
5432 return (xx
< SCM_REAL_VALUE (y
)) ? scm_from_double (xx
) : y
;
5434 else if (SCM_FRACTIONP (y
))
5439 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5442 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
5446 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
5447 (SCM x
, SCM y
, SCM rest
),
5448 "Return the sum of all parameter values. Return 0 if called without\n"
5450 #define FUNC_NAME s_scm_i_sum
5452 while (!scm_is_null (rest
))
5453 { x
= scm_sum (x
, y
);
5455 rest
= scm_cdr (rest
);
5457 return scm_sum (x
, y
);
5461 #define s_sum s_scm_i_sum
5462 #define g_sum g_scm_i_sum
5465 scm_sum (SCM x
, SCM y
)
5467 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5469 if (SCM_NUMBERP (x
)) return x
;
5470 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
5471 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
5474 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5476 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5478 scm_t_inum xx
= SCM_I_INUM (x
);
5479 scm_t_inum yy
= SCM_I_INUM (y
);
5480 scm_t_inum z
= xx
+ yy
;
5481 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_inum2big (z
);
5483 else if (SCM_BIGP (y
))
5488 else if (SCM_REALP (y
))
5490 scm_t_inum xx
= SCM_I_INUM (x
);
5491 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
5493 else if (SCM_COMPLEXP (y
))
5495 scm_t_inum xx
= SCM_I_INUM (x
);
5496 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
5497 SCM_COMPLEX_IMAG (y
));
5499 else if (SCM_FRACTIONP (y
))
5500 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
5501 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
5502 SCM_FRACTION_DENOMINATOR (y
));
5504 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5505 } else if (SCM_BIGP (x
))
5507 if (SCM_I_INUMP (y
))
5512 inum
= SCM_I_INUM (y
);
5515 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5518 SCM result
= scm_i_mkbig ();
5519 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
5520 scm_remember_upto_here_1 (x
);
5521 /* we know the result will have to be a bignum */
5524 return scm_i_normbig (result
);
5528 SCM result
= scm_i_mkbig ();
5529 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
5530 scm_remember_upto_here_1 (x
);
5531 /* we know the result will have to be a bignum */
5534 return scm_i_normbig (result
);
5537 else if (SCM_BIGP (y
))
5539 SCM result
= scm_i_mkbig ();
5540 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5541 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5542 mpz_add (SCM_I_BIG_MPZ (result
),
5545 scm_remember_upto_here_2 (x
, y
);
5546 /* we know the result will have to be a bignum */
5549 return scm_i_normbig (result
);
5551 else if (SCM_REALP (y
))
5553 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
5554 scm_remember_upto_here_1 (x
);
5555 return scm_from_double (result
);
5557 else if (SCM_COMPLEXP (y
))
5559 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
5560 + SCM_COMPLEX_REAL (y
));
5561 scm_remember_upto_here_1 (x
);
5562 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
5564 else if (SCM_FRACTIONP (y
))
5565 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
5566 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
5567 SCM_FRACTION_DENOMINATOR (y
));
5569 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5571 else if (SCM_REALP (x
))
5573 if (SCM_I_INUMP (y
))
5574 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
5575 else if (SCM_BIGP (y
))
5577 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
5578 scm_remember_upto_here_1 (y
);
5579 return scm_from_double (result
);
5581 else if (SCM_REALP (y
))
5582 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
5583 else if (SCM_COMPLEXP (y
))
5584 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
5585 SCM_COMPLEX_IMAG (y
));
5586 else if (SCM_FRACTIONP (y
))
5587 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
5589 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5591 else if (SCM_COMPLEXP (x
))
5593 if (SCM_I_INUMP (y
))
5594 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
5595 SCM_COMPLEX_IMAG (x
));
5596 else if (SCM_BIGP (y
))
5598 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
5599 + SCM_COMPLEX_REAL (x
));
5600 scm_remember_upto_here_1 (y
);
5601 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
5603 else if (SCM_REALP (y
))
5604 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
5605 SCM_COMPLEX_IMAG (x
));
5606 else if (SCM_COMPLEXP (y
))
5607 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
5608 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
5609 else if (SCM_FRACTIONP (y
))
5610 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
5611 SCM_COMPLEX_IMAG (x
));
5613 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5615 else if (SCM_FRACTIONP (x
))
5617 if (SCM_I_INUMP (y
))
5618 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
5619 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
5620 SCM_FRACTION_DENOMINATOR (x
));
5621 else if (SCM_BIGP (y
))
5622 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
5623 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
5624 SCM_FRACTION_DENOMINATOR (x
));
5625 else if (SCM_REALP (y
))
5626 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
5627 else if (SCM_COMPLEXP (y
))
5628 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
5629 SCM_COMPLEX_IMAG (y
));
5630 else if (SCM_FRACTIONP (y
))
5631 /* a/b + c/d = (ad + bc) / bd */
5632 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5633 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
5634 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
5636 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5639 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
5643 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
5645 "Return @math{@var{x}+1}.")
5646 #define FUNC_NAME s_scm_oneplus
5648 return scm_sum (x
, SCM_INUM1
);
5653 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
5654 (SCM x
, SCM y
, SCM rest
),
5655 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
5656 "the sum of all but the first argument are subtracted from the first\n"
5658 #define FUNC_NAME s_scm_i_difference
5660 while (!scm_is_null (rest
))
5661 { x
= scm_difference (x
, y
);
5663 rest
= scm_cdr (rest
);
5665 return scm_difference (x
, y
);
5669 #define s_difference s_scm_i_difference
5670 #define g_difference g_scm_i_difference
5673 scm_difference (SCM x
, SCM y
)
5674 #define FUNC_NAME s_difference
5676 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5679 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
5681 if (SCM_I_INUMP (x
))
5683 scm_t_inum xx
= -SCM_I_INUM (x
);
5684 if (SCM_FIXABLE (xx
))
5685 return SCM_I_MAKINUM (xx
);
5687 return scm_i_inum2big (xx
);
5689 else if (SCM_BIGP (x
))
5690 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
5691 bignum, but negating that gives a fixnum. */
5692 return scm_i_normbig (scm_i_clonebig (x
, 0));
5693 else if (SCM_REALP (x
))
5694 return scm_from_double (-SCM_REAL_VALUE (x
));
5695 else if (SCM_COMPLEXP (x
))
5696 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
5697 -SCM_COMPLEX_IMAG (x
));
5698 else if (SCM_FRACTIONP (x
))
5699 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
5700 SCM_FRACTION_DENOMINATOR (x
));
5702 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
5705 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5707 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5709 scm_t_inum xx
= SCM_I_INUM (x
);
5710 scm_t_inum yy
= SCM_I_INUM (y
);
5711 scm_t_inum z
= xx
- yy
;
5712 if (SCM_FIXABLE (z
))
5713 return SCM_I_MAKINUM (z
);
5715 return scm_i_inum2big (z
);
5717 else if (SCM_BIGP (y
))
5719 /* inum-x - big-y */
5720 scm_t_inum xx
= SCM_I_INUM (x
);
5724 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
5725 bignum, but negating that gives a fixnum. */
5726 return scm_i_normbig (scm_i_clonebig (y
, 0));
5730 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5731 SCM result
= scm_i_mkbig ();
5734 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
5737 /* x - y == -(y + -x) */
5738 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
5739 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5741 scm_remember_upto_here_1 (y
);
5743 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
5744 /* we know the result will have to be a bignum */
5747 return scm_i_normbig (result
);
5750 else if (SCM_REALP (y
))
5752 scm_t_inum xx
= SCM_I_INUM (x
);
5755 * We need to handle x == exact 0
5756 * specially because R6RS states that:
5757 * (- 0.0) ==> -0.0 and
5758 * (- 0.0 0.0) ==> 0.0
5759 * and the scheme compiler changes
5760 * (- 0.0) into (- 0 0.0)
5761 * So we need to treat (- 0 0.0) like (- 0.0).
5762 * At the C level, (-x) is different than (0.0 - x).
5763 * (0.0 - 0.0) ==> 0.0, but (- 0.0) ==> -0.0.
5766 return scm_from_double (- SCM_REAL_VALUE (y
));
5768 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
5770 else if (SCM_COMPLEXP (y
))
5772 scm_t_inum xx
= SCM_I_INUM (x
);
5774 /* We need to handle x == exact 0 specially.
5775 See the comment above (for SCM_REALP (y)) */
5777 return scm_c_make_rectangular (- SCM_COMPLEX_REAL (y
),
5778 - SCM_COMPLEX_IMAG (y
));
5780 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
5781 - SCM_COMPLEX_IMAG (y
));
5783 else if (SCM_FRACTIONP (y
))
5784 /* a - b/c = (ac - b) / c */
5785 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5786 SCM_FRACTION_NUMERATOR (y
)),
5787 SCM_FRACTION_DENOMINATOR (y
));
5789 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5791 else if (SCM_BIGP (x
))
5793 if (SCM_I_INUMP (y
))
5795 /* big-x - inum-y */
5796 scm_t_inum yy
= SCM_I_INUM (y
);
5797 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5799 scm_remember_upto_here_1 (x
);
5801 return (SCM_FIXABLE (-yy
) ?
5802 SCM_I_MAKINUM (-yy
) : scm_from_inum (-yy
));
5805 SCM result
= scm_i_mkbig ();
5808 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
5810 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
5811 scm_remember_upto_here_1 (x
);
5813 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
5814 /* we know the result will have to be a bignum */
5817 return scm_i_normbig (result
);
5820 else if (SCM_BIGP (y
))
5822 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5823 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5824 SCM result
= scm_i_mkbig ();
5825 mpz_sub (SCM_I_BIG_MPZ (result
),
5828 scm_remember_upto_here_2 (x
, y
);
5829 /* we know the result will have to be a bignum */
5830 if ((sgn_x
== 1) && (sgn_y
== -1))
5832 if ((sgn_x
== -1) && (sgn_y
== 1))
5834 return scm_i_normbig (result
);
5836 else if (SCM_REALP (y
))
5838 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
5839 scm_remember_upto_here_1 (x
);
5840 return scm_from_double (result
);
5842 else if (SCM_COMPLEXP (y
))
5844 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
5845 - SCM_COMPLEX_REAL (y
));
5846 scm_remember_upto_here_1 (x
);
5847 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
5849 else if (SCM_FRACTIONP (y
))
5850 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5851 SCM_FRACTION_NUMERATOR (y
)),
5852 SCM_FRACTION_DENOMINATOR (y
));
5853 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5855 else if (SCM_REALP (x
))
5857 if (SCM_I_INUMP (y
))
5858 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
5859 else if (SCM_BIGP (y
))
5861 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
5862 scm_remember_upto_here_1 (x
);
5863 return scm_from_double (result
);
5865 else if (SCM_REALP (y
))
5866 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
5867 else if (SCM_COMPLEXP (y
))
5868 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
5869 -SCM_COMPLEX_IMAG (y
));
5870 else if (SCM_FRACTIONP (y
))
5871 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
5873 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5875 else if (SCM_COMPLEXP (x
))
5877 if (SCM_I_INUMP (y
))
5878 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
5879 SCM_COMPLEX_IMAG (x
));
5880 else if (SCM_BIGP (y
))
5882 double real_part
= (SCM_COMPLEX_REAL (x
)
5883 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
5884 scm_remember_upto_here_1 (x
);
5885 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
5887 else if (SCM_REALP (y
))
5888 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
5889 SCM_COMPLEX_IMAG (x
));
5890 else if (SCM_COMPLEXP (y
))
5891 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
5892 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
5893 else if (SCM_FRACTIONP (y
))
5894 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
5895 SCM_COMPLEX_IMAG (x
));
5897 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5899 else if (SCM_FRACTIONP (x
))
5901 if (SCM_I_INUMP (y
))
5902 /* a/b - c = (a - cb) / b */
5903 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
5904 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
5905 SCM_FRACTION_DENOMINATOR (x
));
5906 else if (SCM_BIGP (y
))
5907 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
5908 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
5909 SCM_FRACTION_DENOMINATOR (x
));
5910 else if (SCM_REALP (y
))
5911 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
5912 else if (SCM_COMPLEXP (y
))
5913 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
5914 -SCM_COMPLEX_IMAG (y
));
5915 else if (SCM_FRACTIONP (y
))
5916 /* a/b - c/d = (ad - bc) / bd */
5917 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5918 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
5919 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
5921 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5924 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
5929 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
5931 "Return @math{@var{x}-1}.")
5932 #define FUNC_NAME s_scm_oneminus
5934 return scm_difference (x
, SCM_INUM1
);
5939 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
5940 (SCM x
, SCM y
, SCM rest
),
5941 "Return the product of all arguments. If called without arguments,\n"
5943 #define FUNC_NAME s_scm_i_product
5945 while (!scm_is_null (rest
))
5946 { x
= scm_product (x
, y
);
5948 rest
= scm_cdr (rest
);
5950 return scm_product (x
, y
);
5954 #define s_product s_scm_i_product
5955 #define g_product g_scm_i_product
5958 scm_product (SCM x
, SCM y
)
5960 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5963 return SCM_I_MAKINUM (1L);
5964 else if (SCM_NUMBERP (x
))
5967 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
5970 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5975 xx
= SCM_I_INUM (x
);
5980 /* exact1 is the universal multiplicative identity */
5984 /* exact0 times a fixnum is exact0: optimize this case */
5985 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5987 /* if the other argument is inexact, the result is inexact,
5988 and we must do the multiplication in order to handle
5989 infinities and NaNs properly. */
5990 else if (SCM_REALP (y
))
5991 return scm_from_double (0.0 * SCM_REAL_VALUE (y
));
5992 else if (SCM_COMPLEXP (y
))
5993 return scm_c_make_rectangular (0.0 * SCM_COMPLEX_REAL (y
),
5994 0.0 * SCM_COMPLEX_IMAG (y
));
5995 /* we've already handled inexact numbers,
5996 so y must be exact, and we return exact0 */
5997 else if (SCM_NUMP (y
))
6000 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6004 * This case is important for more than just optimization.
6005 * It handles the case of negating
6006 * (+ 1 most-positive-fixnum) aka (- most-negative-fixnum),
6007 * which is a bignum that must be changed back into a fixnum.
6008 * Failure to do so will cause the following to return #f:
6009 * (= most-negative-fixnum (* -1 (- most-negative-fixnum)))
6011 return scm_difference(y
, SCM_UNDEFINED
);
6015 if (SCM_LIKELY (SCM_I_INUMP (y
)))
6017 scm_t_inum yy
= SCM_I_INUM (y
);
6018 scm_t_inum kk
= xx
* yy
;
6019 SCM k
= SCM_I_MAKINUM (kk
);
6020 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
6024 SCM result
= scm_i_inum2big (xx
);
6025 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
6026 return scm_i_normbig (result
);
6029 else if (SCM_BIGP (y
))
6031 SCM result
= scm_i_mkbig ();
6032 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
6033 scm_remember_upto_here_1 (y
);
6036 else if (SCM_REALP (y
))
6037 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
6038 else if (SCM_COMPLEXP (y
))
6039 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
6040 xx
* SCM_COMPLEX_IMAG (y
));
6041 else if (SCM_FRACTIONP (y
))
6042 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
6043 SCM_FRACTION_DENOMINATOR (y
));
6045 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6047 else if (SCM_BIGP (x
))
6049 if (SCM_I_INUMP (y
))
6054 else if (SCM_BIGP (y
))
6056 SCM result
= scm_i_mkbig ();
6057 mpz_mul (SCM_I_BIG_MPZ (result
),
6060 scm_remember_upto_here_2 (x
, y
);
6063 else if (SCM_REALP (y
))
6065 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
6066 scm_remember_upto_here_1 (x
);
6067 return scm_from_double (result
);
6069 else if (SCM_COMPLEXP (y
))
6071 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
6072 scm_remember_upto_here_1 (x
);
6073 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
6074 z
* SCM_COMPLEX_IMAG (y
));
6076 else if (SCM_FRACTIONP (y
))
6077 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
6078 SCM_FRACTION_DENOMINATOR (y
));
6080 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6082 else if (SCM_REALP (x
))
6084 if (SCM_I_INUMP (y
))
6089 else if (SCM_BIGP (y
))
6091 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
6092 scm_remember_upto_here_1 (y
);
6093 return scm_from_double (result
);
6095 else if (SCM_REALP (y
))
6096 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
6097 else if (SCM_COMPLEXP (y
))
6098 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
6099 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
6100 else if (SCM_FRACTIONP (y
))
6101 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
6103 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6105 else if (SCM_COMPLEXP (x
))
6107 if (SCM_I_INUMP (y
))
6112 else if (SCM_BIGP (y
))
6114 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6115 scm_remember_upto_here_1 (y
);
6116 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
6117 z
* SCM_COMPLEX_IMAG (x
));
6119 else if (SCM_REALP (y
))
6120 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
6121 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
6122 else if (SCM_COMPLEXP (y
))
6124 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
6125 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
6126 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
6127 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
6129 else if (SCM_FRACTIONP (y
))
6131 double yy
= scm_i_fraction2double (y
);
6132 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
6133 yy
* SCM_COMPLEX_IMAG (x
));
6136 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6138 else if (SCM_FRACTIONP (x
))
6140 if (SCM_I_INUMP (y
))
6141 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
6142 SCM_FRACTION_DENOMINATOR (x
));
6143 else if (SCM_BIGP (y
))
6144 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
6145 SCM_FRACTION_DENOMINATOR (x
));
6146 else if (SCM_REALP (y
))
6147 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
6148 else if (SCM_COMPLEXP (y
))
6150 double xx
= scm_i_fraction2double (x
);
6151 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
6152 xx
* SCM_COMPLEX_IMAG (y
));
6154 else if (SCM_FRACTIONP (y
))
6155 /* a/b * c/d = ac / bd */
6156 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
6157 SCM_FRACTION_NUMERATOR (y
)),
6158 scm_product (SCM_FRACTION_DENOMINATOR (x
),
6159 SCM_FRACTION_DENOMINATOR (y
)));
6161 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6164 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
6167 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
6168 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
6169 #define ALLOW_DIVIDE_BY_ZERO
6170 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
6173 /* The code below for complex division is adapted from the GNU
6174 libstdc++, which adapted it from f2c's libF77, and is subject to
6177 /****************************************************************
6178 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
6180 Permission to use, copy, modify, and distribute this software
6181 and its documentation for any purpose and without fee is hereby
6182 granted, provided that the above copyright notice appear in all
6183 copies and that both that the copyright notice and this
6184 permission notice and warranty disclaimer appear in supporting
6185 documentation, and that the names of AT&T Bell Laboratories or
6186 Bellcore or any of their entities not be used in advertising or
6187 publicity pertaining to distribution of the software without
6188 specific, written prior permission.
6190 AT&T and Bellcore disclaim all warranties with regard to this
6191 software, including all implied warranties of merchantability
6192 and fitness. In no event shall AT&T or Bellcore be liable for
6193 any special, indirect or consequential damages or any damages
6194 whatsoever resulting from loss of use, data or profits, whether
6195 in an action of contract, negligence or other tortious action,
6196 arising out of or in connection with the use or performance of
6198 ****************************************************************/
6200 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
6201 (SCM x
, SCM y
, SCM rest
),
6202 "Divide the first argument by the product of the remaining\n"
6203 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
6205 #define FUNC_NAME s_scm_i_divide
6207 while (!scm_is_null (rest
))
6208 { x
= scm_divide (x
, y
);
6210 rest
= scm_cdr (rest
);
6212 return scm_divide (x
, y
);
6216 #define s_divide s_scm_i_divide
6217 #define g_divide g_scm_i_divide
6220 do_divide (SCM x
, SCM y
, int inexact
)
6221 #define FUNC_NAME s_divide
6225 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
6228 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
6229 else if (SCM_I_INUMP (x
))
6231 scm_t_inum xx
= SCM_I_INUM (x
);
6232 if (xx
== 1 || xx
== -1)
6234 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6236 scm_num_overflow (s_divide
);
6241 return scm_from_double (1.0 / (double) xx
);
6242 else return scm_i_make_ratio (SCM_INUM1
, x
);
6245 else if (SCM_BIGP (x
))
6248 return scm_from_double (1.0 / scm_i_big2dbl (x
));
6249 else return scm_i_make_ratio (SCM_INUM1
, x
);
6251 else if (SCM_REALP (x
))
6253 double xx
= SCM_REAL_VALUE (x
);
6254 #ifndef ALLOW_DIVIDE_BY_ZERO
6256 scm_num_overflow (s_divide
);
6259 return scm_from_double (1.0 / xx
);
6261 else if (SCM_COMPLEXP (x
))
6263 double r
= SCM_COMPLEX_REAL (x
);
6264 double i
= SCM_COMPLEX_IMAG (x
);
6265 if (fabs(r
) <= fabs(i
))
6268 double d
= i
* (1.0 + t
* t
);
6269 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
6274 double d
= r
* (1.0 + t
* t
);
6275 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
6278 else if (SCM_FRACTIONP (x
))
6279 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
6280 SCM_FRACTION_NUMERATOR (x
));
6282 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
6285 if (SCM_LIKELY (SCM_I_INUMP (x
)))
6287 scm_t_inum xx
= SCM_I_INUM (x
);
6288 if (SCM_LIKELY (SCM_I_INUMP (y
)))
6290 scm_t_inum yy
= SCM_I_INUM (y
);
6293 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6294 scm_num_overflow (s_divide
);
6296 return scm_from_double ((double) xx
/ (double) yy
);
6299 else if (xx
% yy
!= 0)
6302 return scm_from_double ((double) xx
/ (double) yy
);
6303 else return scm_i_make_ratio (x
, y
);
6307 scm_t_inum z
= xx
/ yy
;
6308 if (SCM_FIXABLE (z
))
6309 return SCM_I_MAKINUM (z
);
6311 return scm_i_inum2big (z
);
6314 else if (SCM_BIGP (y
))
6317 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
6318 else return scm_i_make_ratio (x
, y
);
6320 else if (SCM_REALP (y
))
6322 double yy
= SCM_REAL_VALUE (y
);
6323 #ifndef ALLOW_DIVIDE_BY_ZERO
6325 scm_num_overflow (s_divide
);
6328 return scm_from_double ((double) xx
/ yy
);
6330 else if (SCM_COMPLEXP (y
))
6333 complex_div
: /* y _must_ be a complex number */
6335 double r
= SCM_COMPLEX_REAL (y
);
6336 double i
= SCM_COMPLEX_IMAG (y
);
6337 if (fabs(r
) <= fabs(i
))
6340 double d
= i
* (1.0 + t
* t
);
6341 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
6346 double d
= r
* (1.0 + t
* t
);
6347 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
6351 else if (SCM_FRACTIONP (y
))
6352 /* a / b/c = ac / b */
6353 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
6354 SCM_FRACTION_NUMERATOR (y
));
6356 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6358 else if (SCM_BIGP (x
))
6360 if (SCM_I_INUMP (y
))
6362 scm_t_inum yy
= SCM_I_INUM (y
);
6365 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6366 scm_num_overflow (s_divide
);
6368 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
6369 scm_remember_upto_here_1 (x
);
6370 return (sgn
== 0) ? scm_nan () : scm_inf ();
6377 /* FIXME: HMM, what are the relative performance issues here?
6378 We need to test. Is it faster on average to test
6379 divisible_p, then perform whichever operation, or is it
6380 faster to perform the integer div opportunistically and
6381 switch to real if there's a remainder? For now we take the
6382 middle ground: test, then if divisible, use the faster div
6385 scm_t_inum abs_yy
= yy
< 0 ? -yy
: yy
;
6386 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
6390 SCM result
= scm_i_mkbig ();
6391 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
6392 scm_remember_upto_here_1 (x
);
6394 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
6395 return scm_i_normbig (result
);
6400 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
6401 else return scm_i_make_ratio (x
, y
);
6405 else if (SCM_BIGP (y
))
6410 /* It's easily possible for the ratio x/y to fit a double
6411 but one or both x and y be too big to fit a double,
6412 hence the use of mpq_get_d rather than converting and
6415 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
6416 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
6417 return scm_from_double (mpq_get_d (q
));
6421 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
6425 SCM result
= scm_i_mkbig ();
6426 mpz_divexact (SCM_I_BIG_MPZ (result
),
6429 scm_remember_upto_here_2 (x
, y
);
6430 return scm_i_normbig (result
);
6433 return scm_i_make_ratio (x
, y
);
6436 else if (SCM_REALP (y
))
6438 double yy
= SCM_REAL_VALUE (y
);
6439 #ifndef ALLOW_DIVIDE_BY_ZERO
6441 scm_num_overflow (s_divide
);
6444 return scm_from_double (scm_i_big2dbl (x
) / yy
);
6446 else if (SCM_COMPLEXP (y
))
6448 a
= scm_i_big2dbl (x
);
6451 else if (SCM_FRACTIONP (y
))
6452 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
6453 SCM_FRACTION_NUMERATOR (y
));
6455 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6457 else if (SCM_REALP (x
))
6459 double rx
= SCM_REAL_VALUE (x
);
6460 if (SCM_I_INUMP (y
))
6462 scm_t_inum yy
= SCM_I_INUM (y
);
6463 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6465 scm_num_overflow (s_divide
);
6468 return scm_from_double (rx
/ (double) yy
);
6470 else if (SCM_BIGP (y
))
6472 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6473 scm_remember_upto_here_1 (y
);
6474 return scm_from_double (rx
/ dby
);
6476 else if (SCM_REALP (y
))
6478 double yy
= SCM_REAL_VALUE (y
);
6479 #ifndef ALLOW_DIVIDE_BY_ZERO
6481 scm_num_overflow (s_divide
);
6484 return scm_from_double (rx
/ yy
);
6486 else if (SCM_COMPLEXP (y
))
6491 else if (SCM_FRACTIONP (y
))
6492 return scm_from_double (rx
/ scm_i_fraction2double (y
));
6494 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6496 else if (SCM_COMPLEXP (x
))
6498 double rx
= SCM_COMPLEX_REAL (x
);
6499 double ix
= SCM_COMPLEX_IMAG (x
);
6500 if (SCM_I_INUMP (y
))
6502 scm_t_inum yy
= SCM_I_INUM (y
);
6503 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6505 scm_num_overflow (s_divide
);
6510 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
6513 else if (SCM_BIGP (y
))
6515 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6516 scm_remember_upto_here_1 (y
);
6517 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
6519 else if (SCM_REALP (y
))
6521 double yy
= SCM_REAL_VALUE (y
);
6522 #ifndef ALLOW_DIVIDE_BY_ZERO
6524 scm_num_overflow (s_divide
);
6527 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
6529 else if (SCM_COMPLEXP (y
))
6531 double ry
= SCM_COMPLEX_REAL (y
);
6532 double iy
= SCM_COMPLEX_IMAG (y
);
6533 if (fabs(ry
) <= fabs(iy
))
6536 double d
= iy
* (1.0 + t
* t
);
6537 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
6542 double d
= ry
* (1.0 + t
* t
);
6543 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
6546 else if (SCM_FRACTIONP (y
))
6548 double yy
= scm_i_fraction2double (y
);
6549 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
6552 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6554 else if (SCM_FRACTIONP (x
))
6556 if (SCM_I_INUMP (y
))
6558 scm_t_inum yy
= SCM_I_INUM (y
);
6559 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6561 scm_num_overflow (s_divide
);
6564 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
6565 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
6567 else if (SCM_BIGP (y
))
6569 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
6570 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
6572 else if (SCM_REALP (y
))
6574 double yy
= SCM_REAL_VALUE (y
);
6575 #ifndef ALLOW_DIVIDE_BY_ZERO
6577 scm_num_overflow (s_divide
);
6580 return scm_from_double (scm_i_fraction2double (x
) / yy
);
6582 else if (SCM_COMPLEXP (y
))
6584 a
= scm_i_fraction2double (x
);
6587 else if (SCM_FRACTIONP (y
))
6588 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
6589 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
6591 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6594 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
6598 scm_divide (SCM x
, SCM y
)
6600 return do_divide (x
, y
, 0);
6603 static SCM
scm_divide2real (SCM x
, SCM y
)
6605 return do_divide (x
, y
, 1);
6611 scm_c_truncate (double x
)
6622 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
6623 half-way case (ie. when x is an integer plus 0.5) going upwards.
6624 Then half-way cases are identified and adjusted down if the
6625 round-upwards didn't give the desired even integer.
6627 "plus_half == result" identifies a half-way case. If plus_half, which is
6628 x + 0.5, is an integer then x must be an integer plus 0.5.
6630 An odd "result" value is identified with result/2 != floor(result/2).
6631 This is done with plus_half, since that value is ready for use sooner in
6632 a pipelined cpu, and we're already requiring plus_half == result.
6634 Note however that we need to be careful when x is big and already an
6635 integer. In that case "x+0.5" may round to an adjacent integer, causing
6636 us to return such a value, incorrectly. For instance if the hardware is
6637 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
6638 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
6639 returned. Or if the hardware is in round-upwards mode, then other bigger
6640 values like say x == 2^128 will see x+0.5 rounding up to the next higher
6641 representable value, 2^128+2^76 (or whatever), again incorrect.
6643 These bad roundings of x+0.5 are avoided by testing at the start whether
6644 x is already an integer. If it is then clearly that's the desired result
6645 already. And if it's not then the exponent must be small enough to allow
6646 an 0.5 to be represented, and hence added without a bad rounding. */
6649 scm_c_round (double x
)
6651 double plus_half
, result
;
6656 plus_half
= x
+ 0.5;
6657 result
= floor (plus_half
);
6658 /* Adjust so that the rounding is towards even. */
6659 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
6664 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
6666 "Round the number @var{x} towards zero.")
6667 #define FUNC_NAME s_scm_truncate_number
6669 if (scm_is_false (scm_negative_p (x
)))
6670 return scm_floor (x
);
6672 return scm_ceiling (x
);
6676 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
6678 "Round the number @var{x} towards the nearest integer. "
6679 "When it is exactly halfway between two integers, "
6680 "round towards the even one.")
6681 #define FUNC_NAME s_scm_round_number
6683 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6685 else if (SCM_REALP (x
))
6686 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
6689 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
6690 single quotient+remainder division then examining to see which way
6691 the rounding should go. */
6692 SCM plus_half
= scm_sum (x
, exactly_one_half
);
6693 SCM result
= scm_floor (plus_half
);
6694 /* Adjust so that the rounding is towards even. */
6695 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
6696 && scm_is_true (scm_odd_p (result
)))
6697 return scm_difference (result
, SCM_INUM1
);
6704 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
6706 "Round the number @var{x} towards minus infinity.")
6707 #define FUNC_NAME s_scm_floor
6709 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6711 else if (SCM_REALP (x
))
6712 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
6713 else if (SCM_FRACTIONP (x
))
6715 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
6716 SCM_FRACTION_DENOMINATOR (x
));
6717 if (scm_is_false (scm_negative_p (x
)))
6719 /* For positive x, rounding towards zero is correct. */
6724 /* For negative x, we need to return q-1 unless x is an
6725 integer. But fractions are never integer, per our
6727 return scm_difference (q
, SCM_INUM1
);
6731 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
6735 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
6737 "Round the number @var{x} towards infinity.")
6738 #define FUNC_NAME s_scm_ceiling
6740 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6742 else if (SCM_REALP (x
))
6743 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
6744 else if (SCM_FRACTIONP (x
))
6746 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
6747 SCM_FRACTION_DENOMINATOR (x
));
6748 if (scm_is_false (scm_positive_p (x
)))
6750 /* For negative x, rounding towards zero is correct. */
6755 /* For positive x, we need to return q+1 unless x is an
6756 integer. But fractions are never integer, per our
6758 return scm_sum (q
, SCM_INUM1
);
6762 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
6766 SCM_PRIMITIVE_GENERIC (scm_expt
, "expt", 2, 0, 0,
6768 "Return @var{x} raised to the power of @var{y}.")
6769 #define FUNC_NAME s_scm_expt
6771 if (scm_is_integer (y
))
6773 if (scm_is_true (scm_exact_p (y
)))
6774 return scm_integer_expt (x
, y
);
6777 /* Here we handle the case where the exponent is an inexact
6778 integer. We make the exponent exact in order to use
6779 scm_integer_expt, and thus avoid the spurious imaginary
6780 parts that may result from round-off errors in the general
6781 e^(y log x) method below (for example when squaring a large
6782 negative number). In this case, we must return an inexact
6783 result for correctness. We also make the base inexact so
6784 that scm_integer_expt will use fast inexact arithmetic
6785 internally. Note that making the base inexact is not
6786 sufficient to guarantee an inexact result, because
6787 scm_integer_expt will return an exact 1 when the exponent
6788 is 0, even if the base is inexact. */
6789 return scm_exact_to_inexact
6790 (scm_integer_expt (scm_exact_to_inexact (x
),
6791 scm_inexact_to_exact (y
)));
6794 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
6796 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
6798 else if (scm_is_complex (x
) && scm_is_complex (y
))
6799 return scm_exp (scm_product (scm_log (x
), y
));
6800 else if (scm_is_complex (x
))
6801 SCM_WTA_DISPATCH_2 (g_scm_expt
, x
, y
, SCM_ARG2
, s_scm_expt
);
6803 SCM_WTA_DISPATCH_2 (g_scm_expt
, x
, y
, SCM_ARG1
, s_scm_expt
);
6807 /* sin/cos/tan/asin/acos/atan
6808 sinh/cosh/tanh/asinh/acosh/atanh
6809 Derived from "Transcen.scm", Complex trancendental functions for SCM.
6810 Written by Jerry D. Hedden, (C) FSF.
6811 See the file `COPYING' for terms applying to this program. */
6813 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
6815 "Compute the sine of @var{z}.")
6816 #define FUNC_NAME s_scm_sin
6818 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6819 return z
; /* sin(exact0) = exact0 */
6820 else if (scm_is_real (z
))
6821 return scm_from_double (sin (scm_to_double (z
)));
6822 else if (SCM_COMPLEXP (z
))
6824 x
= SCM_COMPLEX_REAL (z
);
6825 y
= SCM_COMPLEX_IMAG (z
);
6826 return scm_c_make_rectangular (sin (x
) * cosh (y
),
6827 cos (x
) * sinh (y
));
6830 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
6834 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
6836 "Compute the cosine of @var{z}.")
6837 #define FUNC_NAME s_scm_cos
6839 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6840 return SCM_INUM1
; /* cos(exact0) = exact1 */
6841 else if (scm_is_real (z
))
6842 return scm_from_double (cos (scm_to_double (z
)));
6843 else if (SCM_COMPLEXP (z
))
6845 x
= SCM_COMPLEX_REAL (z
);
6846 y
= SCM_COMPLEX_IMAG (z
);
6847 return scm_c_make_rectangular (cos (x
) * cosh (y
),
6848 -sin (x
) * sinh (y
));
6851 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
6855 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
6857 "Compute the tangent of @var{z}.")
6858 #define FUNC_NAME s_scm_tan
6860 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6861 return z
; /* tan(exact0) = exact0 */
6862 else if (scm_is_real (z
))
6863 return scm_from_double (tan (scm_to_double (z
)));
6864 else if (SCM_COMPLEXP (z
))
6866 x
= 2.0 * SCM_COMPLEX_REAL (z
);
6867 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
6868 w
= cos (x
) + cosh (y
);
6869 #ifndef ALLOW_DIVIDE_BY_ZERO
6871 scm_num_overflow (s_scm_tan
);
6873 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
6876 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
6880 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
6882 "Compute the hyperbolic sine of @var{z}.")
6883 #define FUNC_NAME s_scm_sinh
6885 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6886 return z
; /* sinh(exact0) = exact0 */
6887 else if (scm_is_real (z
))
6888 return scm_from_double (sinh (scm_to_double (z
)));
6889 else if (SCM_COMPLEXP (z
))
6891 x
= SCM_COMPLEX_REAL (z
);
6892 y
= SCM_COMPLEX_IMAG (z
);
6893 return scm_c_make_rectangular (sinh (x
) * cos (y
),
6894 cosh (x
) * sin (y
));
6897 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
6901 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
6903 "Compute the hyperbolic cosine of @var{z}.")
6904 #define FUNC_NAME s_scm_cosh
6906 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6907 return SCM_INUM1
; /* cosh(exact0) = exact1 */
6908 else if (scm_is_real (z
))
6909 return scm_from_double (cosh (scm_to_double (z
)));
6910 else if (SCM_COMPLEXP (z
))
6912 x
= SCM_COMPLEX_REAL (z
);
6913 y
= SCM_COMPLEX_IMAG (z
);
6914 return scm_c_make_rectangular (cosh (x
) * cos (y
),
6915 sinh (x
) * sin (y
));
6918 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
6922 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
6924 "Compute the hyperbolic tangent of @var{z}.")
6925 #define FUNC_NAME s_scm_tanh
6927 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6928 return z
; /* tanh(exact0) = exact0 */
6929 else if (scm_is_real (z
))
6930 return scm_from_double (tanh (scm_to_double (z
)));
6931 else if (SCM_COMPLEXP (z
))
6933 x
= 2.0 * SCM_COMPLEX_REAL (z
);
6934 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
6935 w
= cosh (x
) + cos (y
);
6936 #ifndef ALLOW_DIVIDE_BY_ZERO
6938 scm_num_overflow (s_scm_tanh
);
6940 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
6943 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
6947 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
6949 "Compute the arc sine of @var{z}.")
6950 #define FUNC_NAME s_scm_asin
6952 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6953 return z
; /* asin(exact0) = exact0 */
6954 else if (scm_is_real (z
))
6956 double w
= scm_to_double (z
);
6957 if (w
>= -1.0 && w
<= 1.0)
6958 return scm_from_double (asin (w
));
6960 return scm_product (scm_c_make_rectangular (0, -1),
6961 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
6963 else if (SCM_COMPLEXP (z
))
6965 x
= SCM_COMPLEX_REAL (z
);
6966 y
= SCM_COMPLEX_IMAG (z
);
6967 return scm_product (scm_c_make_rectangular (0, -1),
6968 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
6971 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
6975 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
6977 "Compute the arc cosine of @var{z}.")
6978 #define FUNC_NAME s_scm_acos
6980 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM1
)))
6981 return SCM_INUM0
; /* acos(exact1) = exact0 */
6982 else if (scm_is_real (z
))
6984 double w
= scm_to_double (z
);
6985 if (w
>= -1.0 && w
<= 1.0)
6986 return scm_from_double (acos (w
));
6988 return scm_sum (scm_from_double (acos (0.0)),
6989 scm_product (scm_c_make_rectangular (0, 1),
6990 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
6992 else if (SCM_COMPLEXP (z
))
6994 x
= SCM_COMPLEX_REAL (z
);
6995 y
= SCM_COMPLEX_IMAG (z
);
6996 return scm_sum (scm_from_double (acos (0.0)),
6997 scm_product (scm_c_make_rectangular (0, 1),
6998 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
7001 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
7005 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
7007 "With one argument, compute the arc tangent of @var{z}.\n"
7008 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
7009 "using the sign of @var{z} and @var{y} to determine the quadrant.")
7010 #define FUNC_NAME s_scm_atan
7014 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
7015 return z
; /* atan(exact0) = exact0 */
7016 else if (scm_is_real (z
))
7017 return scm_from_double (atan (scm_to_double (z
)));
7018 else if (SCM_COMPLEXP (z
))
7021 v
= SCM_COMPLEX_REAL (z
);
7022 w
= SCM_COMPLEX_IMAG (z
);
7023 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
7024 scm_c_make_rectangular (v
, w
+ 1.0))),
7025 scm_c_make_rectangular (0, 2));
7028 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
7030 else if (scm_is_real (z
))
7032 if (scm_is_real (y
))
7033 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
7035 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
7038 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
7042 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
7044 "Compute the inverse hyperbolic sine of @var{z}.")
7045 #define FUNC_NAME s_scm_sys_asinh
7047 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
7048 return z
; /* asinh(exact0) = exact0 */
7049 else if (scm_is_real (z
))
7050 return scm_from_double (asinh (scm_to_double (z
)));
7051 else if (scm_is_number (z
))
7052 return scm_log (scm_sum (z
,
7053 scm_sqrt (scm_sum (scm_product (z
, z
),
7056 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
7060 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
7062 "Compute the inverse hyperbolic cosine of @var{z}.")
7063 #define FUNC_NAME s_scm_sys_acosh
7065 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM1
)))
7066 return SCM_INUM0
; /* acosh(exact1) = exact0 */
7067 else if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
7068 return scm_from_double (acosh (scm_to_double (z
)));
7069 else if (scm_is_number (z
))
7070 return scm_log (scm_sum (z
,
7071 scm_sqrt (scm_difference (scm_product (z
, z
),
7074 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
7078 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
7080 "Compute the inverse hyperbolic tangent of @var{z}.")
7081 #define FUNC_NAME s_scm_sys_atanh
7083 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
7084 return z
; /* atanh(exact0) = exact0 */
7085 else if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
7086 return scm_from_double (atanh (scm_to_double (z
)));
7087 else if (scm_is_number (z
))
7088 return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1
, z
),
7089 scm_difference (SCM_INUM1
, z
))),
7092 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
7097 scm_c_make_rectangular (double re
, double im
)
7101 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
7103 SCM_SET_CELL_TYPE (z
, scm_tc16_complex
);
7104 SCM_COMPLEX_REAL (z
) = re
;
7105 SCM_COMPLEX_IMAG (z
) = im
;
7109 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
7110 (SCM real_part
, SCM imaginary_part
),
7111 "Return a complex number constructed of the given @var{real-part} "
7112 "and @var{imaginary-part} parts.")
7113 #define FUNC_NAME s_scm_make_rectangular
7115 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
7116 SCM_ARG1
, FUNC_NAME
, "real");
7117 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
7118 SCM_ARG2
, FUNC_NAME
, "real");
7120 /* Return a real if and only if the imaginary_part is an _exact_ 0 */
7121 if (scm_is_eq (imaginary_part
, SCM_INUM0
))
7124 return scm_c_make_rectangular (scm_to_double (real_part
),
7125 scm_to_double (imaginary_part
));
7130 scm_c_make_polar (double mag
, double ang
)
7134 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
7135 use it on Glibc-based systems that have it (it's a GNU extension). See
7136 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
7138 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
7139 sincos (ang
, &s
, &c
);
7145 /* If s and c are NaNs, this indicates that the angle is a NaN,
7146 infinite, or perhaps simply too large to determine its value
7147 mod 2*pi. However, we know something that the floating-point
7148 implementation doesn't know: We know that s and c are finite.
7149 Therefore, if the magnitude is zero, return a complex zero.
7151 The reason we check for the NaNs instead of using this case
7152 whenever mag == 0.0 is because when the angle is known, we'd
7153 like to return the correct kind of non-real complex zero:
7154 +0.0+0.0i, -0.0+0.0i, -0.0-0.0i, or +0.0-0.0i, depending
7155 on which quadrant the angle is in.
7157 if (SCM_UNLIKELY (isnan(s
)) && isnan(c
) && (mag
== 0.0))
7158 return scm_c_make_rectangular (0.0, 0.0);
7160 return scm_c_make_rectangular (mag
* c
, mag
* s
);
7163 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
7165 "Return the complex number @var{mag} * e^(i * @var{ang}).")
7166 #define FUNC_NAME s_scm_make_polar
7168 SCM_ASSERT_TYPE (scm_is_real (mag
), mag
, SCM_ARG1
, FUNC_NAME
, "real");
7169 SCM_ASSERT_TYPE (scm_is_real (ang
), ang
, SCM_ARG2
, FUNC_NAME
, "real");
7171 /* If mag is exact0, return exact0 */
7172 if (scm_is_eq (mag
, SCM_INUM0
))
7174 /* Return a real if ang is exact0 */
7175 else if (scm_is_eq (ang
, SCM_INUM0
))
7178 return scm_c_make_polar (scm_to_double (mag
), scm_to_double (ang
));
7183 SCM_PRIMITIVE_GENERIC (scm_real_part
, "real-part", 1, 0, 0,
7185 "Return the real part of the number @var{z}.")
7186 #define FUNC_NAME s_scm_real_part
7188 if (SCM_COMPLEXP (z
))
7189 return scm_from_double (SCM_COMPLEX_REAL (z
));
7190 else if (SCM_I_INUMP (z
) || SCM_BIGP (z
) || SCM_REALP (z
) || SCM_FRACTIONP (z
))
7193 SCM_WTA_DISPATCH_1 (g_scm_real_part
, z
, SCM_ARG1
, s_scm_real_part
);
7198 SCM_PRIMITIVE_GENERIC (scm_imag_part
, "imag-part", 1, 0, 0,
7200 "Return the imaginary part of the number @var{z}.")
7201 #define FUNC_NAME s_scm_imag_part
7203 if (SCM_COMPLEXP (z
))
7204 return scm_from_double (SCM_COMPLEX_IMAG (z
));
7205 else if (SCM_I_INUMP (z
) || SCM_REALP (z
) || SCM_BIGP (z
) || SCM_FRACTIONP (z
))
7208 SCM_WTA_DISPATCH_1 (g_scm_imag_part
, z
, SCM_ARG1
, s_scm_imag_part
);
7212 SCM_PRIMITIVE_GENERIC (scm_numerator
, "numerator", 1, 0, 0,
7214 "Return the numerator of the number @var{z}.")
7215 #define FUNC_NAME s_scm_numerator
7217 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7219 else if (SCM_FRACTIONP (z
))
7220 return SCM_FRACTION_NUMERATOR (z
);
7221 else if (SCM_REALP (z
))
7222 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
7224 SCM_WTA_DISPATCH_1 (g_scm_numerator
, z
, SCM_ARG1
, s_scm_numerator
);
7229 SCM_PRIMITIVE_GENERIC (scm_denominator
, "denominator", 1, 0, 0,
7231 "Return the denominator of the number @var{z}.")
7232 #define FUNC_NAME s_scm_denominator
7234 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7236 else if (SCM_FRACTIONP (z
))
7237 return SCM_FRACTION_DENOMINATOR (z
);
7238 else if (SCM_REALP (z
))
7239 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
7241 SCM_WTA_DISPATCH_1 (g_scm_denominator
, z
, SCM_ARG1
, s_scm_denominator
);
7246 SCM_PRIMITIVE_GENERIC (scm_magnitude
, "magnitude", 1, 0, 0,
7248 "Return the magnitude of the number @var{z}. This is the same as\n"
7249 "@code{abs} for real arguments, but also allows complex numbers.")
7250 #define FUNC_NAME s_scm_magnitude
7252 if (SCM_I_INUMP (z
))
7254 scm_t_inum zz
= SCM_I_INUM (z
);
7257 else if (SCM_POSFIXABLE (-zz
))
7258 return SCM_I_MAKINUM (-zz
);
7260 return scm_i_inum2big (-zz
);
7262 else if (SCM_BIGP (z
))
7264 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
7265 scm_remember_upto_here_1 (z
);
7267 return scm_i_clonebig (z
, 0);
7271 else if (SCM_REALP (z
))
7272 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
7273 else if (SCM_COMPLEXP (z
))
7274 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
7275 else if (SCM_FRACTIONP (z
))
7277 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
7279 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
7280 SCM_FRACTION_DENOMINATOR (z
));
7283 SCM_WTA_DISPATCH_1 (g_scm_magnitude
, z
, SCM_ARG1
, s_scm_magnitude
);
7288 SCM_PRIMITIVE_GENERIC (scm_angle
, "angle", 1, 0, 0,
7290 "Return the angle of the complex number @var{z}.")
7291 #define FUNC_NAME s_scm_angle
7293 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
7294 flo0 to save allocating a new flonum with scm_from_double each time.
7295 But if atan2 follows the floating point rounding mode, then the value
7296 is not a constant. Maybe it'd be close enough though. */
7297 if (SCM_I_INUMP (z
))
7299 if (SCM_I_INUM (z
) >= 0)
7302 return scm_from_double (atan2 (0.0, -1.0));
7304 else if (SCM_BIGP (z
))
7306 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
7307 scm_remember_upto_here_1 (z
);
7309 return scm_from_double (atan2 (0.0, -1.0));
7313 else if (SCM_REALP (z
))
7315 if (SCM_REAL_VALUE (z
) >= 0)
7318 return scm_from_double (atan2 (0.0, -1.0));
7320 else if (SCM_COMPLEXP (z
))
7321 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
7322 else if (SCM_FRACTIONP (z
))
7324 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
7326 else return scm_from_double (atan2 (0.0, -1.0));
7329 SCM_WTA_DISPATCH_1 (g_scm_angle
, z
, SCM_ARG1
, s_scm_angle
);
7334 SCM_PRIMITIVE_GENERIC (scm_exact_to_inexact
, "exact->inexact", 1, 0, 0,
7336 "Convert the number @var{z} to its inexact representation.\n")
7337 #define FUNC_NAME s_scm_exact_to_inexact
7339 if (SCM_I_INUMP (z
))
7340 return scm_from_double ((double) SCM_I_INUM (z
));
7341 else if (SCM_BIGP (z
))
7342 return scm_from_double (scm_i_big2dbl (z
));
7343 else if (SCM_FRACTIONP (z
))
7344 return scm_from_double (scm_i_fraction2double (z
));
7345 else if (SCM_INEXACTP (z
))
7348 SCM_WTA_DISPATCH_1 (g_scm_exact_to_inexact
, z
, 1, s_scm_exact_to_inexact
);
7353 SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
7355 "Return an exact number that is numerically closest to @var{z}.")
7356 #define FUNC_NAME s_scm_inexact_to_exact
7358 if (SCM_I_INUMP (z
) || SCM_BIGP (z
) || SCM_FRACTIONP (z
))
7365 val
= SCM_REAL_VALUE (z
);
7366 else if (SCM_COMPLEXP (z
) && SCM_COMPLEX_IMAG (z
) == 0.0)
7367 val
= SCM_COMPLEX_REAL (z
);
7369 SCM_WTA_DISPATCH_1 (g_scm_inexact_to_exact
, z
, 1, s_scm_inexact_to_exact
);
7371 if (!SCM_LIKELY (DOUBLE_IS_FINITE (val
)))
7372 SCM_OUT_OF_RANGE (1, z
);
7379 mpq_set_d (frac
, val
);
7380 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
7381 scm_i_mpz2num (mpq_denref (frac
)));
7383 /* When scm_i_make_ratio throws, we leak the memory allocated
7393 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
7395 "Returns the @emph{simplest} rational number differing\n"
7396 "from @var{x} by no more than @var{eps}.\n"
7398 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
7399 "exact result when both its arguments are exact. Thus, you might need\n"
7400 "to use @code{inexact->exact} on the arguments.\n"
7403 "(rationalize (inexact->exact 1.2) 1/100)\n"
7406 #define FUNC_NAME s_scm_rationalize
7408 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
7409 SCM_ASSERT_TYPE (scm_is_real (eps
), eps
, SCM_ARG2
, FUNC_NAME
, "real");
7410 eps
= scm_abs (eps
);
7411 if (scm_is_false (scm_positive_p (eps
)))
7413 /* eps is either zero or a NaN */
7414 if (scm_is_true (scm_nan_p (eps
)))
7416 else if (SCM_INEXACTP (eps
))
7417 return scm_exact_to_inexact (x
);
7421 else if (scm_is_false (scm_finite_p (eps
)))
7423 if (scm_is_true (scm_finite_p (x
)))
7428 else if (scm_is_false (scm_finite_p (x
))) /* checks for both inf and nan */
7430 else if (scm_is_false (scm_less_p (scm_floor (scm_sum (x
, eps
)),
7431 scm_ceiling (scm_difference (x
, eps
)))))
7433 /* There's an integer within range; we want the one closest to zero */
7434 if (scm_is_false (scm_less_p (eps
, scm_abs (x
))))
7436 /* zero is within range */
7437 if (SCM_INEXACTP (x
) || SCM_INEXACTP (eps
))
7442 else if (scm_is_true (scm_positive_p (x
)))
7443 return scm_ceiling (scm_difference (x
, eps
));
7445 return scm_floor (scm_sum (x
, eps
));
7449 /* Use continued fractions to find closest ratio. All
7450 arithmetic is done with exact numbers.
7453 SCM ex
= scm_inexact_to_exact (x
);
7454 SCM int_part
= scm_floor (ex
);
7456 SCM a1
= SCM_INUM0
, a2
= SCM_INUM1
, a
= SCM_INUM0
;
7457 SCM b1
= SCM_INUM1
, b2
= SCM_INUM0
, b
= SCM_INUM0
;
7461 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
7462 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
7464 /* We stop after a million iterations just to be absolutely sure
7465 that we don't go into an infinite loop. The process normally
7466 converges after less than a dozen iterations.
7469 while (++i
< 1000000)
7471 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
7472 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
7473 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
7475 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
7476 eps
))) /* abs(x-a/b) <= eps */
7478 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
7479 if (SCM_INEXACTP (x
) || SCM_INEXACTP (eps
))
7480 return scm_exact_to_inexact (res
);
7484 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
7486 tt
= scm_floor (rx
); /* tt = floor (rx) */
7492 scm_num_overflow (s_scm_rationalize
);
7497 /* conversion functions */
7500 scm_is_integer (SCM val
)
7502 return scm_is_true (scm_integer_p (val
));
7506 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
7508 if (SCM_I_INUMP (val
))
7510 scm_t_signed_bits n
= SCM_I_INUM (val
);
7511 return n
>= min
&& n
<= max
;
7513 else if (SCM_BIGP (val
))
7515 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
7517 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
7519 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
7521 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
7522 return n
>= min
&& n
<= max
;
7532 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
7533 > CHAR_BIT
*sizeof (scm_t_uintmax
))
7536 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
7537 SCM_I_BIG_MPZ (val
));
7539 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
7551 return n
>= min
&& n
<= max
;
7559 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
7561 if (SCM_I_INUMP (val
))
7563 scm_t_signed_bits n
= SCM_I_INUM (val
);
7564 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
7566 else if (SCM_BIGP (val
))
7568 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
7570 else if (max
<= ULONG_MAX
)
7572 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
7574 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
7575 return n
>= min
&& n
<= max
;
7585 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
7588 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
7589 > CHAR_BIT
*sizeof (scm_t_uintmax
))
7592 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
7593 SCM_I_BIG_MPZ (val
));
7595 return n
>= min
&& n
<= max
;
7603 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
7605 scm_error (scm_out_of_range_key
,
7607 "Value out of range ~S to ~S: ~S",
7608 scm_list_3 (min
, max
, bad_val
),
7609 scm_list_1 (bad_val
));
7612 #define TYPE scm_t_intmax
7613 #define TYPE_MIN min
7614 #define TYPE_MAX max
7615 #define SIZEOF_TYPE 0
7616 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
7617 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
7618 #include "libguile/conv-integer.i.c"
7620 #define TYPE scm_t_uintmax
7621 #define TYPE_MIN min
7622 #define TYPE_MAX max
7623 #define SIZEOF_TYPE 0
7624 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
7625 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
7626 #include "libguile/conv-uinteger.i.c"
7628 #define TYPE scm_t_int8
7629 #define TYPE_MIN SCM_T_INT8_MIN
7630 #define TYPE_MAX SCM_T_INT8_MAX
7631 #define SIZEOF_TYPE 1
7632 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
7633 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
7634 #include "libguile/conv-integer.i.c"
7636 #define TYPE scm_t_uint8
7638 #define TYPE_MAX SCM_T_UINT8_MAX
7639 #define SIZEOF_TYPE 1
7640 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
7641 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
7642 #include "libguile/conv-uinteger.i.c"
7644 #define TYPE scm_t_int16
7645 #define TYPE_MIN SCM_T_INT16_MIN
7646 #define TYPE_MAX SCM_T_INT16_MAX
7647 #define SIZEOF_TYPE 2
7648 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
7649 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
7650 #include "libguile/conv-integer.i.c"
7652 #define TYPE scm_t_uint16
7654 #define TYPE_MAX SCM_T_UINT16_MAX
7655 #define SIZEOF_TYPE 2
7656 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
7657 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
7658 #include "libguile/conv-uinteger.i.c"
7660 #define TYPE scm_t_int32
7661 #define TYPE_MIN SCM_T_INT32_MIN
7662 #define TYPE_MAX SCM_T_INT32_MAX
7663 #define SIZEOF_TYPE 4
7664 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
7665 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
7666 #include "libguile/conv-integer.i.c"
7668 #define TYPE scm_t_uint32
7670 #define TYPE_MAX SCM_T_UINT32_MAX
7671 #define SIZEOF_TYPE 4
7672 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
7673 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
7674 #include "libguile/conv-uinteger.i.c"
7676 #define TYPE scm_t_wchar
7677 #define TYPE_MIN (scm_t_int32)-1
7678 #define TYPE_MAX (scm_t_int32)0x10ffff
7679 #define SIZEOF_TYPE 4
7680 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
7681 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
7682 #include "libguile/conv-integer.i.c"
7684 #define TYPE scm_t_int64
7685 #define TYPE_MIN SCM_T_INT64_MIN
7686 #define TYPE_MAX SCM_T_INT64_MAX
7687 #define SIZEOF_TYPE 8
7688 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
7689 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
7690 #include "libguile/conv-integer.i.c"
7692 #define TYPE scm_t_uint64
7694 #define TYPE_MAX SCM_T_UINT64_MAX
7695 #define SIZEOF_TYPE 8
7696 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
7697 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
7698 #include "libguile/conv-uinteger.i.c"
7701 scm_to_mpz (SCM val
, mpz_t rop
)
7703 if (SCM_I_INUMP (val
))
7704 mpz_set_si (rop
, SCM_I_INUM (val
));
7705 else if (SCM_BIGP (val
))
7706 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
7708 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
7712 scm_from_mpz (mpz_t val
)
7714 return scm_i_mpz2num (val
);
7718 scm_is_real (SCM val
)
7720 return scm_is_true (scm_real_p (val
));
7724 scm_is_rational (SCM val
)
7726 return scm_is_true (scm_rational_p (val
));
7730 scm_to_double (SCM val
)
7732 if (SCM_I_INUMP (val
))
7733 return SCM_I_INUM (val
);
7734 else if (SCM_BIGP (val
))
7735 return scm_i_big2dbl (val
);
7736 else if (SCM_FRACTIONP (val
))
7737 return scm_i_fraction2double (val
);
7738 else if (SCM_REALP (val
))
7739 return SCM_REAL_VALUE (val
);
7741 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
7745 scm_from_double (double val
)
7749 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
7751 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
7752 SCM_REAL_VALUE (z
) = val
;
7757 #if SCM_ENABLE_DEPRECATED == 1
7760 scm_num2float (SCM num
, unsigned long pos
, const char *s_caller
)
7762 scm_c_issue_deprecation_warning
7763 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
7767 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
7771 scm_out_of_range (NULL
, num
);
7774 return scm_to_double (num
);
7778 scm_num2double (SCM num
, unsigned long pos
, const char *s_caller
)
7780 scm_c_issue_deprecation_warning
7781 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
7785 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
7789 scm_out_of_range (NULL
, num
);
7792 return scm_to_double (num
);
7798 scm_is_complex (SCM val
)
7800 return scm_is_true (scm_complex_p (val
));
7804 scm_c_real_part (SCM z
)
7806 if (SCM_COMPLEXP (z
))
7807 return SCM_COMPLEX_REAL (z
);
7810 /* Use the scm_real_part to get proper error checking and
7813 return scm_to_double (scm_real_part (z
));
7818 scm_c_imag_part (SCM z
)
7820 if (SCM_COMPLEXP (z
))
7821 return SCM_COMPLEX_IMAG (z
);
7824 /* Use the scm_imag_part to get proper error checking and
7825 dispatching. The result will almost always be 0.0, but not
7828 return scm_to_double (scm_imag_part (z
));
7833 scm_c_magnitude (SCM z
)
7835 return scm_to_double (scm_magnitude (z
));
7841 return scm_to_double (scm_angle (z
));
7845 scm_is_number (SCM z
)
7847 return scm_is_true (scm_number_p (z
));
7851 /* In the following functions we dispatch to the real-arg funcs like log()
7852 when we know the arg is real, instead of just handing everything to
7853 clog() for instance. This is in case clog() doesn't optimize for a
7854 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
7855 well use it to go straight to the applicable C func. */
7857 SCM_PRIMITIVE_GENERIC (scm_log
, "log", 1, 0, 0,
7859 "Return the natural logarithm of @var{z}.")
7860 #define FUNC_NAME s_scm_log
7862 if (SCM_COMPLEXP (z
))
7864 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
7865 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
7867 double re
= SCM_COMPLEX_REAL (z
);
7868 double im
= SCM_COMPLEX_IMAG (z
);
7869 return scm_c_make_rectangular (log (hypot (re
, im
)),
7873 else if (SCM_NUMBERP (z
))
7875 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
7876 although the value itself overflows. */
7877 double re
= scm_to_double (z
);
7878 double l
= log (fabs (re
));
7880 return scm_from_double (l
);
7882 return scm_c_make_rectangular (l
, M_PI
);
7885 SCM_WTA_DISPATCH_1 (g_scm_log
, z
, 1, s_scm_log
);
7890 SCM_PRIMITIVE_GENERIC (scm_log10
, "log10", 1, 0, 0,
7892 "Return the base 10 logarithm of @var{z}.")
7893 #define FUNC_NAME s_scm_log10
7895 if (SCM_COMPLEXP (z
))
7897 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
7898 clog() and a multiply by M_LOG10E, rather than the fallback
7899 log10+hypot+atan2.) */
7900 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
7901 && defined SCM_COMPLEX_VALUE
7902 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
7904 double re
= SCM_COMPLEX_REAL (z
);
7905 double im
= SCM_COMPLEX_IMAG (z
);
7906 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
7907 M_LOG10E
* atan2 (im
, re
));
7910 else if (SCM_NUMBERP (z
))
7912 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
7913 although the value itself overflows. */
7914 double re
= scm_to_double (z
);
7915 double l
= log10 (fabs (re
));
7917 return scm_from_double (l
);
7919 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
7922 SCM_WTA_DISPATCH_1 (g_scm_log10
, z
, 1, s_scm_log10
);
7927 SCM_PRIMITIVE_GENERIC (scm_exp
, "exp", 1, 0, 0,
7929 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
7930 "base of natural logarithms (2.71828@dots{}).")
7931 #define FUNC_NAME s_scm_exp
7933 if (SCM_COMPLEXP (z
))
7935 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
7936 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
7938 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
7939 SCM_COMPLEX_IMAG (z
));
7942 else if (SCM_NUMBERP (z
))
7944 /* When z is a negative bignum the conversion to double overflows,
7945 giving -infinity, but that's ok, the exp is still 0.0. */
7946 return scm_from_double (exp (scm_to_double (z
)));
7949 SCM_WTA_DISPATCH_1 (g_scm_exp
, z
, 1, s_scm_exp
);
7954 SCM_PRIMITIVE_GENERIC (scm_sqrt
, "sqrt", 1, 0, 0,
7956 "Return the square root of @var{z}. Of the two possible roots\n"
7957 "(positive and negative), the one with the a positive real part\n"
7958 "is returned, or if that's zero then a positive imaginary part.\n"
7962 "(sqrt 9.0) @result{} 3.0\n"
7963 "(sqrt -9.0) @result{} 0.0+3.0i\n"
7964 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
7965 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
7967 #define FUNC_NAME s_scm_sqrt
7969 if (SCM_COMPLEXP (z
))
7971 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
7972 && defined SCM_COMPLEX_VALUE
7973 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (z
)));
7975 double re
= SCM_COMPLEX_REAL (z
);
7976 double im
= SCM_COMPLEX_IMAG (z
);
7977 return scm_c_make_polar (sqrt (hypot (re
, im
)),
7978 0.5 * atan2 (im
, re
));
7981 else if (SCM_NUMBERP (z
))
7983 double xx
= scm_to_double (z
);
7985 return scm_c_make_rectangular (0.0, sqrt (-xx
));
7987 return scm_from_double (sqrt (xx
));
7990 SCM_WTA_DISPATCH_1 (g_scm_sqrt
, z
, 1, s_scm_sqrt
);
8001 mpz_init_set_si (z_negative_one
, -1);
8003 /* It may be possible to tune the performance of some algorithms by using
8004 * the following constants to avoid the creation of bignums. Please, before
8005 * using these values, remember the two rules of program optimization:
8006 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
8007 scm_c_define ("most-positive-fixnum",
8008 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
8009 scm_c_define ("most-negative-fixnum",
8010 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
8012 scm_add_feature ("complex");
8013 scm_add_feature ("inexact");
8014 flo0
= scm_from_double (0.0);
8016 /* determine floating point precision */
8017 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
8019 init_dblprec(&scm_dblprec
[i
-2],i
);
8020 init_fx_radix(fx_per_radix
[i
-2],i
);
8023 /* hard code precision for base 10 if the preprocessor tells us to... */
8024 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
8027 exactly_one_half
= scm_divide (SCM_INUM1
, SCM_I_MAKINUM (2));
8028 #include "libguile/numbers.x"