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:
3838 * * Each function keeps a local index variable 'idx' that points at the
3839 * current position within the parsed string. The global index is only
3840 * updated if the function could parse the corresponding syntactic unit
3843 * * Similarly, the functions keep track of indicators of inexactness ('#',
3844 * '.' or exponents) using local variables ('hash_seen', 'x').
3846 * * Sequences of digits are parsed into temporary variables holding fixnums.
3847 * Only if these fixnums would overflow, the result variables are updated
3848 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
3849 * the temporary variables holding the fixnums are cleared, and the process
3850 * starts over again. If for example fixnums were able to store five decimal
3851 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
3852 * and the result was computed as 12345 * 100000 + 67890. In other words,
3853 * only every five digits two bignum operations were performed.
3855 * Notes on the handling of exactness specifiers:
3857 * When parsing non-real complex numbers, we apply exactness specifiers on
3858 * per-component basis, as is done in PLT Scheme. For complex numbers
3859 * written in rectangular form, exactness specifiers are applied to the
3860 * real and imaginary parts before calling scm_make_rectangular. For
3861 * complex numbers written in polar form, exactness specifiers are applied
3862 * to the magnitude and angle before calling scm_make_polar.
3864 * There are two kinds of exactness specifiers: forced and implicit. A
3865 * forced exactness specifier is a "#e" or "#i" prefix at the beginning of
3866 * the entire number, and applies to both components of a complex number.
3867 * "#e" causes each component to be made exact, and "#i" causes each
3868 * component to be made inexact. If no forced exactness specifier is
3869 * present, then the exactness of each component is determined
3870 * independently by the presence or absence of a decimal point or hash mark
3871 * within that component. If a decimal point or hash mark is present, the
3872 * component is made inexact, otherwise it is made exact.
3874 * After the exactness specifiers have been applied to each component, they
3875 * are passed to either scm_make_rectangular or scm_make_polar to produce
3876 * the final result. Note that this will result in a real number if the
3877 * imaginary part, magnitude, or angle is an exact 0.
3879 * For example, (string->number "#i5.0+0i") does the equivalent of:
3881 * (make-rectangular (exact->inexact 5) (exact->inexact 0))
3884 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
3886 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
3888 /* Caller is responsible for checking that the return value is in range
3889 for the given radix, which should be <= 36. */
3891 char_decimal_value (scm_t_uint32 c
)
3893 /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
3894 that's certainly above any valid decimal, so we take advantage of
3895 that to elide some tests. */
3896 unsigned int d
= (unsigned int) uc_decimal_value (c
);
3898 /* If that failed, try extended hexadecimals, then. Only accept ascii
3903 if (c
>= (scm_t_uint32
) 'a')
3904 d
= c
- (scm_t_uint32
)'a' + 10U;
3910 mem2uinteger (SCM mem
, unsigned int *p_idx
,
3911 unsigned int radix
, enum t_exactness
*p_exactness
)
3913 unsigned int idx
= *p_idx
;
3914 unsigned int hash_seen
= 0;
3915 scm_t_bits shift
= 1;
3917 unsigned int digit_value
;
3920 size_t len
= scm_i_string_length (mem
);
3925 c
= scm_i_string_ref (mem
, idx
);
3926 digit_value
= char_decimal_value (c
);
3927 if (digit_value
>= radix
)
3931 result
= SCM_I_MAKINUM (digit_value
);
3934 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
3944 digit_value
= char_decimal_value (c
);
3945 /* This check catches non-decimals in addition to out-of-range
3947 if (digit_value
>= radix
)
3952 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
3954 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3956 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3963 shift
= shift
* radix
;
3964 add
= add
* radix
+ digit_value
;
3969 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3971 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3975 *p_exactness
= INEXACT
;
3981 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
3982 * covers the parts of the rules that start at a potential point. The value
3983 * of the digits up to the point have been parsed by the caller and are given
3984 * in variable result. The content of *p_exactness indicates, whether a hash
3985 * has already been seen in the digits before the point.
3988 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
3991 mem2decimal_from_point (SCM result
, SCM mem
,
3992 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
3994 unsigned int idx
= *p_idx
;
3995 enum t_exactness x
= *p_exactness
;
3996 size_t len
= scm_i_string_length (mem
);
4001 if (scm_i_string_ref (mem
, idx
) == '.')
4003 scm_t_bits shift
= 1;
4005 unsigned int digit_value
;
4006 SCM big_shift
= SCM_INUM1
;
4011 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
4012 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4017 digit_value
= DIGIT2UINT (c
);
4028 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
4030 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
4031 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
4033 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
4041 add
= add
* 10 + digit_value
;
4047 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
4048 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
4049 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
4052 result
= scm_divide (result
, big_shift
);
4054 /* We've seen a decimal point, thus the value is implicitly inexact. */
4066 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
4068 switch (scm_i_string_ref (mem
, idx
))
4080 c
= scm_i_string_ref (mem
, idx
);
4088 c
= scm_i_string_ref (mem
, idx
);
4097 c
= scm_i_string_ref (mem
, idx
);
4102 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4106 exponent
= DIGIT2UINT (c
);
4109 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
4110 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4113 if (exponent
<= SCM_MAXEXP
)
4114 exponent
= exponent
* 10 + DIGIT2UINT (c
);
4120 if (exponent
> SCM_MAXEXP
)
4122 size_t exp_len
= idx
- start
;
4123 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
4124 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
4125 scm_out_of_range ("string->number", exp_num
);
4128 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
4130 result
= scm_product (result
, e
);
4132 result
= scm_divide2real (result
, e
);
4134 /* We've seen an exponent, thus the value is implicitly inexact. */
4152 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
4155 mem2ureal (SCM mem
, unsigned int *p_idx
,
4156 unsigned int radix
, enum t_exactness forced_x
)
4158 unsigned int idx
= *p_idx
;
4160 size_t len
= scm_i_string_length (mem
);
4162 /* Start off believing that the number will be exact. This changes
4163 to INEXACT if we see a decimal point or a hash. */
4164 enum t_exactness implicit_x
= EXACT
;
4169 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
4175 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
4177 /* Cobble up the fractional part. We might want to set the
4178 NaN's mantissa from it. */
4180 mem2uinteger (mem
, &idx
, 10, &implicit_x
);
4185 if (scm_i_string_ref (mem
, idx
) == '.')
4189 else if (idx
+ 1 == len
)
4191 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
4194 result
= mem2decimal_from_point (SCM_INUM0
, mem
,
4195 p_idx
, &implicit_x
);
4201 uinteger
= mem2uinteger (mem
, &idx
, radix
, &implicit_x
);
4202 if (scm_is_false (uinteger
))
4207 else if (scm_i_string_ref (mem
, idx
) == '/')
4215 divisor
= mem2uinteger (mem
, &idx
, radix
, &implicit_x
);
4216 if (scm_is_false (divisor
))
4219 /* both are int/big here, I assume */
4220 result
= scm_i_make_ratio (uinteger
, divisor
);
4222 else if (radix
== 10)
4224 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &implicit_x
);
4225 if (scm_is_false (result
))
4237 if (SCM_INEXACTP (result
))
4238 return scm_inexact_to_exact (result
);
4242 if (SCM_INEXACTP (result
))
4245 return scm_exact_to_inexact (result
);
4247 if (implicit_x
== INEXACT
)
4249 if (SCM_INEXACTP (result
))
4252 return scm_exact_to_inexact (result
);
4258 /* We should never get here */
4259 scm_syserror ("mem2ureal");
4263 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
4266 mem2complex (SCM mem
, unsigned int idx
,
4267 unsigned int radix
, enum t_exactness forced_x
)
4272 size_t len
= scm_i_string_length (mem
);
4277 c
= scm_i_string_ref (mem
, idx
);
4292 ureal
= mem2ureal (mem
, &idx
, radix
, forced_x
);
4293 if (scm_is_false (ureal
))
4295 /* input must be either +i or -i */
4300 if (scm_i_string_ref (mem
, idx
) == 'i'
4301 || scm_i_string_ref (mem
, idx
) == 'I')
4307 return scm_make_rectangular (SCM_INUM0
, SCM_I_MAKINUM (sign
));
4314 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
4315 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
4320 c
= scm_i_string_ref (mem
, idx
);
4324 /* either +<ureal>i or -<ureal>i */
4331 return scm_make_rectangular (SCM_INUM0
, ureal
);
4334 /* polar input: <real>@<real>. */
4345 c
= scm_i_string_ref (mem
, idx
);
4363 angle
= mem2ureal (mem
, &idx
, radix
, forced_x
);
4364 if (scm_is_false (angle
))
4369 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
4370 angle
= scm_difference (angle
, SCM_UNDEFINED
);
4372 result
= scm_make_polar (ureal
, angle
);
4377 /* expecting input matching <real>[+-]<ureal>?i */
4384 int sign
= (c
== '+') ? 1 : -1;
4385 SCM imag
= mem2ureal (mem
, &idx
, radix
, forced_x
);
4387 if (scm_is_false (imag
))
4388 imag
= SCM_I_MAKINUM (sign
);
4389 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
4390 imag
= scm_difference (imag
, SCM_UNDEFINED
);
4394 if (scm_i_string_ref (mem
, idx
) != 'i'
4395 && scm_i_string_ref (mem
, idx
) != 'I')
4402 return scm_make_rectangular (ureal
, imag
);
4411 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
4413 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
4416 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
4418 unsigned int idx
= 0;
4419 unsigned int radix
= NO_RADIX
;
4420 enum t_exactness forced_x
= NO_EXACTNESS
;
4421 size_t len
= scm_i_string_length (mem
);
4423 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
4424 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
4426 switch (scm_i_string_ref (mem
, idx
+ 1))
4429 if (radix
!= NO_RADIX
)
4434 if (radix
!= NO_RADIX
)
4439 if (forced_x
!= NO_EXACTNESS
)
4444 if (forced_x
!= NO_EXACTNESS
)
4449 if (radix
!= NO_RADIX
)
4454 if (radix
!= NO_RADIX
)
4464 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
4465 if (radix
== NO_RADIX
)
4466 radix
= default_radix
;
4468 return mem2complex (mem
, idx
, radix
, forced_x
);
4472 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
4473 unsigned int default_radix
)
4475 SCM str
= scm_from_locale_stringn (mem
, len
);
4477 return scm_i_string_to_number (str
, default_radix
);
4481 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
4482 (SCM string
, SCM radix
),
4483 "Return a number of the maximally precise representation\n"
4484 "expressed by the given @var{string}. @var{radix} must be an\n"
4485 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
4486 "is a default radix that may be overridden by an explicit radix\n"
4487 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
4488 "supplied, then the default radix is 10. If string is not a\n"
4489 "syntactically valid notation for a number, then\n"
4490 "@code{string->number} returns @code{#f}.")
4491 #define FUNC_NAME s_scm_string_to_number
4495 SCM_VALIDATE_STRING (1, string
);
4497 if (SCM_UNBNDP (radix
))
4500 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
4502 answer
= scm_i_string_to_number (string
, base
);
4503 scm_remember_upto_here_1 (string
);
4509 /*** END strs->nums ***/
4512 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
4514 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
4516 #define FUNC_NAME s_scm_number_p
4518 return scm_from_bool (SCM_NUMBERP (x
));
4522 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
4524 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
4525 "otherwise. Note that the sets of real, rational and integer\n"
4526 "values form subsets of the set of complex numbers, i. e. the\n"
4527 "predicate will also be fulfilled if @var{x} is a real,\n"
4528 "rational or integer number.")
4529 #define FUNC_NAME s_scm_complex_p
4531 /* all numbers are complex. */
4532 return scm_number_p (x
);
4536 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
4538 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
4539 "otherwise. Note that the set of integer values forms a subset of\n"
4540 "the set of real numbers, i. e. the predicate will also be\n"
4541 "fulfilled if @var{x} is an integer number.")
4542 #define FUNC_NAME s_scm_real_p
4544 return scm_from_bool
4545 (SCM_I_INUMP (x
) || SCM_REALP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
));
4549 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
4551 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
4552 "otherwise. Note that the set of integer values forms a subset of\n"
4553 "the set of rational numbers, i. e. the predicate will also be\n"
4554 "fulfilled if @var{x} is an integer number.")
4555 #define FUNC_NAME s_scm_rational_p
4557 if (SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
))
4559 else if (SCM_REALP (x
))
4560 /* due to their limited precision, finite floating point numbers are
4561 rational as well. (finite means neither infinity nor a NaN) */
4562 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
4568 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
4570 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
4572 #define FUNC_NAME s_scm_integer_p
4574 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4576 else if (SCM_REALP (x
))
4578 double val
= SCM_REAL_VALUE (x
);
4579 return scm_from_bool (!isinf (val
) && (val
== floor (val
)));
4587 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
4588 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
4589 (SCM x
, SCM y
, SCM rest
),
4590 "Return @code{#t} if all parameters are numerically equal.")
4591 #define FUNC_NAME s_scm_i_num_eq_p
4593 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4595 while (!scm_is_null (rest
))
4597 if (scm_is_false (scm_num_eq_p (x
, y
)))
4601 rest
= scm_cdr (rest
);
4603 return scm_num_eq_p (x
, y
);
4607 scm_num_eq_p (SCM x
, SCM y
)
4610 if (SCM_I_INUMP (x
))
4612 scm_t_signed_bits xx
= SCM_I_INUM (x
);
4613 if (SCM_I_INUMP (y
))
4615 scm_t_signed_bits yy
= SCM_I_INUM (y
);
4616 return scm_from_bool (xx
== yy
);
4618 else if (SCM_BIGP (y
))
4620 else if (SCM_REALP (y
))
4622 /* On a 32-bit system an inum fits a double, we can cast the inum
4623 to a double and compare.
4625 But on a 64-bit system an inum is bigger than a double and
4626 casting it to a double (call that dxx) will round. dxx is at
4627 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
4628 an integer and fits a long. So we cast yy to a long and
4629 compare with plain xx.
4631 An alternative (for any size system actually) would be to check
4632 yy is an integer (with floor) and is in range of an inum
4633 (compare against appropriate powers of 2) then test
4634 xx==(scm_t_signed_bits)yy. It's just a matter of which
4635 casts/comparisons might be fastest or easiest for the cpu. */
4637 double yy
= SCM_REAL_VALUE (y
);
4638 return scm_from_bool ((double) xx
== yy
4639 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
4640 || xx
== (scm_t_signed_bits
) yy
));
4642 else if (SCM_COMPLEXP (y
))
4643 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
4644 && (0.0 == SCM_COMPLEX_IMAG (y
)));
4645 else if (SCM_FRACTIONP (y
))
4648 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4650 else if (SCM_BIGP (x
))
4652 if (SCM_I_INUMP (y
))
4654 else if (SCM_BIGP (y
))
4656 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4657 scm_remember_upto_here_2 (x
, y
);
4658 return scm_from_bool (0 == cmp
);
4660 else if (SCM_REALP (y
))
4663 if (isnan (SCM_REAL_VALUE (y
)))
4665 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
4666 scm_remember_upto_here_1 (x
);
4667 return scm_from_bool (0 == cmp
);
4669 else if (SCM_COMPLEXP (y
))
4672 if (0.0 != SCM_COMPLEX_IMAG (y
))
4674 if (isnan (SCM_COMPLEX_REAL (y
)))
4676 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
4677 scm_remember_upto_here_1 (x
);
4678 return scm_from_bool (0 == cmp
);
4680 else if (SCM_FRACTIONP (y
))
4683 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4685 else if (SCM_REALP (x
))
4687 double xx
= SCM_REAL_VALUE (x
);
4688 if (SCM_I_INUMP (y
))
4690 /* see comments with inum/real above */
4691 scm_t_signed_bits yy
= SCM_I_INUM (y
);
4692 return scm_from_bool (xx
== (double) yy
4693 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
4694 || (scm_t_signed_bits
) xx
== yy
));
4696 else if (SCM_BIGP (y
))
4699 if (isnan (SCM_REAL_VALUE (x
)))
4701 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
4702 scm_remember_upto_here_1 (y
);
4703 return scm_from_bool (0 == cmp
);
4705 else if (SCM_REALP (y
))
4706 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
4707 else if (SCM_COMPLEXP (y
))
4708 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
4709 && (0.0 == SCM_COMPLEX_IMAG (y
)));
4710 else if (SCM_FRACTIONP (y
))
4712 double xx
= SCM_REAL_VALUE (x
);
4716 return scm_from_bool (xx
< 0.0);
4717 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4721 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4723 else if (SCM_COMPLEXP (x
))
4725 if (SCM_I_INUMP (y
))
4726 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
4727 && (SCM_COMPLEX_IMAG (x
) == 0.0));
4728 else if (SCM_BIGP (y
))
4731 if (0.0 != SCM_COMPLEX_IMAG (x
))
4733 if (isnan (SCM_COMPLEX_REAL (x
)))
4735 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
4736 scm_remember_upto_here_1 (y
);
4737 return scm_from_bool (0 == cmp
);
4739 else if (SCM_REALP (y
))
4740 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
4741 && (SCM_COMPLEX_IMAG (x
) == 0.0));
4742 else if (SCM_COMPLEXP (y
))
4743 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
4744 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
4745 else if (SCM_FRACTIONP (y
))
4748 if (SCM_COMPLEX_IMAG (x
) != 0.0)
4750 xx
= SCM_COMPLEX_REAL (x
);
4754 return scm_from_bool (xx
< 0.0);
4755 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4759 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4761 else if (SCM_FRACTIONP (x
))
4763 if (SCM_I_INUMP (y
))
4765 else if (SCM_BIGP (y
))
4767 else if (SCM_REALP (y
))
4769 double yy
= SCM_REAL_VALUE (y
);
4773 return scm_from_bool (0.0 < yy
);
4774 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4777 else if (SCM_COMPLEXP (y
))
4780 if (SCM_COMPLEX_IMAG (y
) != 0.0)
4782 yy
= SCM_COMPLEX_REAL (y
);
4786 return scm_from_bool (0.0 < yy
);
4787 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4790 else if (SCM_FRACTIONP (y
))
4791 return scm_i_fraction_equalp (x
, y
);
4793 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4796 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
4800 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
4801 done are good for inums, but for bignums an answer can almost always be
4802 had by just examining a few high bits of the operands, as done by GMP in
4803 mpq_cmp. flonum/frac compares likewise, but with the slight complication
4804 of the float exponent to take into account. */
4806 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
4807 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
4808 (SCM x
, SCM y
, SCM rest
),
4809 "Return @code{#t} if the list of parameters is monotonically\n"
4811 #define FUNC_NAME s_scm_i_num_less_p
4813 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4815 while (!scm_is_null (rest
))
4817 if (scm_is_false (scm_less_p (x
, y
)))
4821 rest
= scm_cdr (rest
);
4823 return scm_less_p (x
, y
);
4827 scm_less_p (SCM x
, SCM y
)
4830 if (SCM_I_INUMP (x
))
4832 scm_t_inum xx
= SCM_I_INUM (x
);
4833 if (SCM_I_INUMP (y
))
4835 scm_t_inum yy
= SCM_I_INUM (y
);
4836 return scm_from_bool (xx
< yy
);
4838 else if (SCM_BIGP (y
))
4840 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4841 scm_remember_upto_here_1 (y
);
4842 return scm_from_bool (sgn
> 0);
4844 else if (SCM_REALP (y
))
4845 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
4846 else if (SCM_FRACTIONP (y
))
4848 /* "x < a/b" becomes "x*b < a" */
4850 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
4851 y
= SCM_FRACTION_NUMERATOR (y
);
4855 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4857 else if (SCM_BIGP (x
))
4859 if (SCM_I_INUMP (y
))
4861 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4862 scm_remember_upto_here_1 (x
);
4863 return scm_from_bool (sgn
< 0);
4865 else if (SCM_BIGP (y
))
4867 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4868 scm_remember_upto_here_2 (x
, y
);
4869 return scm_from_bool (cmp
< 0);
4871 else if (SCM_REALP (y
))
4874 if (isnan (SCM_REAL_VALUE (y
)))
4876 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
4877 scm_remember_upto_here_1 (x
);
4878 return scm_from_bool (cmp
< 0);
4880 else if (SCM_FRACTIONP (y
))
4883 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4885 else if (SCM_REALP (x
))
4887 if (SCM_I_INUMP (y
))
4888 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
4889 else if (SCM_BIGP (y
))
4892 if (isnan (SCM_REAL_VALUE (x
)))
4894 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
4895 scm_remember_upto_here_1 (y
);
4896 return scm_from_bool (cmp
> 0);
4898 else if (SCM_REALP (y
))
4899 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
4900 else if (SCM_FRACTIONP (y
))
4902 double xx
= SCM_REAL_VALUE (x
);
4906 return scm_from_bool (xx
< 0.0);
4907 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4911 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4913 else if (SCM_FRACTIONP (x
))
4915 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
4917 /* "a/b < y" becomes "a < y*b" */
4918 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
4919 x
= SCM_FRACTION_NUMERATOR (x
);
4922 else if (SCM_REALP (y
))
4924 double yy
= SCM_REAL_VALUE (y
);
4928 return scm_from_bool (0.0 < yy
);
4929 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4932 else if (SCM_FRACTIONP (y
))
4934 /* "a/b < c/d" becomes "a*d < c*b" */
4935 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
4936 SCM_FRACTION_DENOMINATOR (y
));
4937 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
4938 SCM_FRACTION_DENOMINATOR (x
));
4944 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4947 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
4951 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
4952 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
4953 (SCM x
, SCM y
, SCM rest
),
4954 "Return @code{#t} if the list of parameters is monotonically\n"
4956 #define FUNC_NAME s_scm_i_num_gr_p
4958 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4960 while (!scm_is_null (rest
))
4962 if (scm_is_false (scm_gr_p (x
, y
)))
4966 rest
= scm_cdr (rest
);
4968 return scm_gr_p (x
, y
);
4971 #define FUNC_NAME s_scm_i_num_gr_p
4973 scm_gr_p (SCM x
, SCM y
)
4975 if (!SCM_NUMBERP (x
))
4976 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
4977 else if (!SCM_NUMBERP (y
))
4978 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
4980 return scm_less_p (y
, x
);
4985 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
4986 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
4987 (SCM x
, SCM y
, SCM rest
),
4988 "Return @code{#t} if the list of parameters is monotonically\n"
4990 #define FUNC_NAME s_scm_i_num_leq_p
4992 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4994 while (!scm_is_null (rest
))
4996 if (scm_is_false (scm_leq_p (x
, y
)))
5000 rest
= scm_cdr (rest
);
5002 return scm_leq_p (x
, y
);
5005 #define FUNC_NAME s_scm_i_num_leq_p
5007 scm_leq_p (SCM x
, SCM y
)
5009 if (!SCM_NUMBERP (x
))
5010 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
5011 else if (!SCM_NUMBERP (y
))
5012 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
5013 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
5016 return scm_not (scm_less_p (y
, x
));
5021 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
5022 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
5023 (SCM x
, SCM y
, SCM rest
),
5024 "Return @code{#t} if the list of parameters is monotonically\n"
5026 #define FUNC_NAME s_scm_i_num_geq_p
5028 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
5030 while (!scm_is_null (rest
))
5032 if (scm_is_false (scm_geq_p (x
, y
)))
5036 rest
= scm_cdr (rest
);
5038 return scm_geq_p (x
, y
);
5041 #define FUNC_NAME s_scm_i_num_geq_p
5043 scm_geq_p (SCM x
, SCM y
)
5045 if (!SCM_NUMBERP (x
))
5046 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
5047 else if (!SCM_NUMBERP (y
))
5048 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
5049 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
5052 return scm_not (scm_less_p (x
, y
));
5057 SCM_PRIMITIVE_GENERIC (scm_zero_p
, "zero?", 1, 0, 0,
5059 "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
5061 #define FUNC_NAME s_scm_zero_p
5063 if (SCM_I_INUMP (z
))
5064 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
5065 else if (SCM_BIGP (z
))
5067 else if (SCM_REALP (z
))
5068 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
5069 else if (SCM_COMPLEXP (z
))
5070 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
5071 && SCM_COMPLEX_IMAG (z
) == 0.0);
5072 else if (SCM_FRACTIONP (z
))
5075 SCM_WTA_DISPATCH_1 (g_scm_zero_p
, z
, SCM_ARG1
, s_scm_zero_p
);
5080 SCM_PRIMITIVE_GENERIC (scm_positive_p
, "positive?", 1, 0, 0,
5082 "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
5084 #define FUNC_NAME s_scm_positive_p
5086 if (SCM_I_INUMP (x
))
5087 return scm_from_bool (SCM_I_INUM (x
) > 0);
5088 else if (SCM_BIGP (x
))
5090 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5091 scm_remember_upto_here_1 (x
);
5092 return scm_from_bool (sgn
> 0);
5094 else if (SCM_REALP (x
))
5095 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
5096 else if (SCM_FRACTIONP (x
))
5097 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
5099 SCM_WTA_DISPATCH_1 (g_scm_positive_p
, x
, SCM_ARG1
, s_scm_positive_p
);
5104 SCM_PRIMITIVE_GENERIC (scm_negative_p
, "negative?", 1, 0, 0,
5106 "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
5108 #define FUNC_NAME s_scm_negative_p
5110 if (SCM_I_INUMP (x
))
5111 return scm_from_bool (SCM_I_INUM (x
) < 0);
5112 else if (SCM_BIGP (x
))
5114 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5115 scm_remember_upto_here_1 (x
);
5116 return scm_from_bool (sgn
< 0);
5118 else if (SCM_REALP (x
))
5119 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
5120 else if (SCM_FRACTIONP (x
))
5121 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
5123 SCM_WTA_DISPATCH_1 (g_scm_negative_p
, x
, SCM_ARG1
, s_scm_negative_p
);
5128 /* scm_min and scm_max return an inexact when either argument is inexact, as
5129 required by r5rs. On that basis, for exact/inexact combinations the
5130 exact is converted to inexact to compare and possibly return. This is
5131 unlike scm_less_p above which takes some trouble to preserve all bits in
5132 its test, such trouble is not required for min and max. */
5134 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
5135 (SCM x
, SCM y
, SCM rest
),
5136 "Return the maximum of all parameter values.")
5137 #define FUNC_NAME s_scm_i_max
5139 while (!scm_is_null (rest
))
5140 { x
= scm_max (x
, y
);
5142 rest
= scm_cdr (rest
);
5144 return scm_max (x
, y
);
5148 #define s_max s_scm_i_max
5149 #define g_max g_scm_i_max
5152 scm_max (SCM x
, SCM y
)
5157 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
5158 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
5161 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
5164 if (SCM_I_INUMP (x
))
5166 scm_t_inum xx
= SCM_I_INUM (x
);
5167 if (SCM_I_INUMP (y
))
5169 scm_t_inum yy
= SCM_I_INUM (y
);
5170 return (xx
< yy
) ? y
: x
;
5172 else if (SCM_BIGP (y
))
5174 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5175 scm_remember_upto_here_1 (y
);
5176 return (sgn
< 0) ? x
: y
;
5178 else if (SCM_REALP (y
))
5181 double yyd
= SCM_REAL_VALUE (y
);
5184 return scm_from_double (xxd
);
5185 /* If y is a NaN, then "==" is false and we return the NaN */
5186 else if (SCM_LIKELY (!(xxd
== yyd
)))
5188 /* Handle signed zeroes properly */
5194 else if (SCM_FRACTIONP (y
))
5197 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
5200 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5202 else if (SCM_BIGP (x
))
5204 if (SCM_I_INUMP (y
))
5206 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5207 scm_remember_upto_here_1 (x
);
5208 return (sgn
< 0) ? y
: x
;
5210 else if (SCM_BIGP (y
))
5212 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
5213 scm_remember_upto_here_2 (x
, y
);
5214 return (cmp
> 0) ? x
: y
;
5216 else if (SCM_REALP (y
))
5218 /* if y==NaN then xx>yy is false, so we return the NaN y */
5221 xx
= scm_i_big2dbl (x
);
5222 yy
= SCM_REAL_VALUE (y
);
5223 return (xx
> yy
? scm_from_double (xx
) : y
);
5225 else if (SCM_FRACTIONP (y
))
5230 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5232 else if (SCM_REALP (x
))
5234 if (SCM_I_INUMP (y
))
5236 scm_t_inum yy
= SCM_I_INUM (y
);
5237 double xxd
= SCM_REAL_VALUE (x
);
5241 return scm_from_double (yyd
);
5242 /* If x is a NaN, then "==" is false and we return the NaN */
5243 else if (SCM_LIKELY (!(xxd
== yyd
)))
5245 /* Handle signed zeroes properly */
5251 else if (SCM_BIGP (y
))
5256 else if (SCM_REALP (y
))
5258 double xx
= SCM_REAL_VALUE (x
);
5259 double yy
= SCM_REAL_VALUE (y
);
5261 /* For purposes of max: +inf.0 > nan > everything else, per R6RS */
5264 else if (SCM_LIKELY (xx
< yy
))
5266 /* If neither (xx > yy) nor (xx < yy), then
5267 either they're equal or one is a NaN */
5268 else if (SCM_UNLIKELY (isnan (xx
)))
5269 return DOUBLE_IS_POSITIVE_INFINITY (yy
) ? y
: x
;
5270 else if (SCM_UNLIKELY (isnan (yy
)))
5271 return DOUBLE_IS_POSITIVE_INFINITY (xx
) ? x
: y
;
5272 /* xx == yy, but handle signed zeroes properly */
5273 else if (double_is_non_negative_zero (yy
))
5278 else if (SCM_FRACTIONP (y
))
5280 double yy
= scm_i_fraction2double (y
);
5281 double xx
= SCM_REAL_VALUE (x
);
5282 return (xx
< yy
) ? scm_from_double (yy
) : x
;
5285 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5287 else if (SCM_FRACTIONP (x
))
5289 if (SCM_I_INUMP (y
))
5293 else if (SCM_BIGP (y
))
5297 else if (SCM_REALP (y
))
5299 double xx
= scm_i_fraction2double (x
);
5300 /* if y==NaN then ">" is false, so we return the NaN y */
5301 return (xx
> SCM_REAL_VALUE (y
)) ? scm_from_double (xx
) : y
;
5303 else if (SCM_FRACTIONP (y
))
5308 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5311 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
5315 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
5316 (SCM x
, SCM y
, SCM rest
),
5317 "Return the minimum of all parameter values.")
5318 #define FUNC_NAME s_scm_i_min
5320 while (!scm_is_null (rest
))
5321 { x
= scm_min (x
, y
);
5323 rest
= scm_cdr (rest
);
5325 return scm_min (x
, y
);
5329 #define s_min s_scm_i_min
5330 #define g_min g_scm_i_min
5333 scm_min (SCM x
, SCM y
)
5338 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
5339 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
5342 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
5345 if (SCM_I_INUMP (x
))
5347 scm_t_inum xx
= SCM_I_INUM (x
);
5348 if (SCM_I_INUMP (y
))
5350 scm_t_inum yy
= SCM_I_INUM (y
);
5351 return (xx
< yy
) ? x
: y
;
5353 else if (SCM_BIGP (y
))
5355 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5356 scm_remember_upto_here_1 (y
);
5357 return (sgn
< 0) ? y
: x
;
5359 else if (SCM_REALP (y
))
5362 /* if y==NaN then "<" is false and we return NaN */
5363 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
5365 else if (SCM_FRACTIONP (y
))
5368 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
5371 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5373 else if (SCM_BIGP (x
))
5375 if (SCM_I_INUMP (y
))
5377 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5378 scm_remember_upto_here_1 (x
);
5379 return (sgn
< 0) ? x
: y
;
5381 else if (SCM_BIGP (y
))
5383 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
5384 scm_remember_upto_here_2 (x
, y
);
5385 return (cmp
> 0) ? y
: x
;
5387 else if (SCM_REALP (y
))
5389 /* if y==NaN then xx<yy is false, so we return the NaN y */
5392 xx
= scm_i_big2dbl (x
);
5393 yy
= SCM_REAL_VALUE (y
);
5394 return (xx
< yy
? scm_from_double (xx
) : y
);
5396 else if (SCM_FRACTIONP (y
))
5401 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5403 else if (SCM_REALP (x
))
5405 if (SCM_I_INUMP (y
))
5407 double z
= SCM_I_INUM (y
);
5408 /* if x==NaN then "<" is false and we return NaN */
5409 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
5411 else if (SCM_BIGP (y
))
5416 else if (SCM_REALP (y
))
5418 double xx
= SCM_REAL_VALUE (x
);
5419 double yy
= SCM_REAL_VALUE (y
);
5421 /* For purposes of min: -inf.0 < nan < everything else, per R6RS */
5424 else if (SCM_LIKELY (xx
> yy
))
5426 /* If neither (xx < yy) nor (xx > yy), then
5427 either they're equal or one is a NaN */
5428 else if (SCM_UNLIKELY (isnan (xx
)))
5429 return DOUBLE_IS_NEGATIVE_INFINITY (yy
) ? y
: x
;
5430 else if (SCM_UNLIKELY (isnan (yy
)))
5431 return DOUBLE_IS_NEGATIVE_INFINITY (xx
) ? x
: y
;
5432 /* xx == yy, but handle signed zeroes properly */
5433 else if (double_is_non_negative_zero (xx
))
5438 else if (SCM_FRACTIONP (y
))
5440 double yy
= scm_i_fraction2double (y
);
5441 double xx
= SCM_REAL_VALUE (x
);
5442 return (yy
< xx
) ? scm_from_double (yy
) : x
;
5445 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5447 else if (SCM_FRACTIONP (x
))
5449 if (SCM_I_INUMP (y
))
5453 else if (SCM_BIGP (y
))
5457 else if (SCM_REALP (y
))
5459 double xx
= scm_i_fraction2double (x
);
5460 /* if y==NaN then "<" is false, so we return the NaN y */
5461 return (xx
< SCM_REAL_VALUE (y
)) ? scm_from_double (xx
) : y
;
5463 else if (SCM_FRACTIONP (y
))
5468 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5471 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
5475 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
5476 (SCM x
, SCM y
, SCM rest
),
5477 "Return the sum of all parameter values. Return 0 if called without\n"
5479 #define FUNC_NAME s_scm_i_sum
5481 while (!scm_is_null (rest
))
5482 { x
= scm_sum (x
, y
);
5484 rest
= scm_cdr (rest
);
5486 return scm_sum (x
, y
);
5490 #define s_sum s_scm_i_sum
5491 #define g_sum g_scm_i_sum
5494 scm_sum (SCM x
, SCM y
)
5496 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5498 if (SCM_NUMBERP (x
)) return x
;
5499 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
5500 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
5503 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5505 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5507 scm_t_inum xx
= SCM_I_INUM (x
);
5508 scm_t_inum yy
= SCM_I_INUM (y
);
5509 scm_t_inum z
= xx
+ yy
;
5510 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_inum2big (z
);
5512 else if (SCM_BIGP (y
))
5517 else if (SCM_REALP (y
))
5519 scm_t_inum xx
= SCM_I_INUM (x
);
5520 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
5522 else if (SCM_COMPLEXP (y
))
5524 scm_t_inum xx
= SCM_I_INUM (x
);
5525 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
5526 SCM_COMPLEX_IMAG (y
));
5528 else if (SCM_FRACTIONP (y
))
5529 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
5530 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
5531 SCM_FRACTION_DENOMINATOR (y
));
5533 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5534 } else if (SCM_BIGP (x
))
5536 if (SCM_I_INUMP (y
))
5541 inum
= SCM_I_INUM (y
);
5544 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5547 SCM result
= scm_i_mkbig ();
5548 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
5549 scm_remember_upto_here_1 (x
);
5550 /* we know the result will have to be a bignum */
5553 return scm_i_normbig (result
);
5557 SCM result
= scm_i_mkbig ();
5558 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
5559 scm_remember_upto_here_1 (x
);
5560 /* we know the result will have to be a bignum */
5563 return scm_i_normbig (result
);
5566 else if (SCM_BIGP (y
))
5568 SCM result
= scm_i_mkbig ();
5569 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5570 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5571 mpz_add (SCM_I_BIG_MPZ (result
),
5574 scm_remember_upto_here_2 (x
, y
);
5575 /* we know the result will have to be a bignum */
5578 return scm_i_normbig (result
);
5580 else if (SCM_REALP (y
))
5582 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
5583 scm_remember_upto_here_1 (x
);
5584 return scm_from_double (result
);
5586 else if (SCM_COMPLEXP (y
))
5588 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
5589 + SCM_COMPLEX_REAL (y
));
5590 scm_remember_upto_here_1 (x
);
5591 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
5593 else if (SCM_FRACTIONP (y
))
5594 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
5595 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
5596 SCM_FRACTION_DENOMINATOR (y
));
5598 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5600 else if (SCM_REALP (x
))
5602 if (SCM_I_INUMP (y
))
5603 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
5604 else if (SCM_BIGP (y
))
5606 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
5607 scm_remember_upto_here_1 (y
);
5608 return scm_from_double (result
);
5610 else if (SCM_REALP (y
))
5611 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
5612 else if (SCM_COMPLEXP (y
))
5613 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
5614 SCM_COMPLEX_IMAG (y
));
5615 else if (SCM_FRACTIONP (y
))
5616 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
5618 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5620 else if (SCM_COMPLEXP (x
))
5622 if (SCM_I_INUMP (y
))
5623 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
5624 SCM_COMPLEX_IMAG (x
));
5625 else if (SCM_BIGP (y
))
5627 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
5628 + SCM_COMPLEX_REAL (x
));
5629 scm_remember_upto_here_1 (y
);
5630 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
5632 else if (SCM_REALP (y
))
5633 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
5634 SCM_COMPLEX_IMAG (x
));
5635 else if (SCM_COMPLEXP (y
))
5636 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
5637 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
5638 else if (SCM_FRACTIONP (y
))
5639 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
5640 SCM_COMPLEX_IMAG (x
));
5642 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5644 else if (SCM_FRACTIONP (x
))
5646 if (SCM_I_INUMP (y
))
5647 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
5648 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
5649 SCM_FRACTION_DENOMINATOR (x
));
5650 else if (SCM_BIGP (y
))
5651 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
5652 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
5653 SCM_FRACTION_DENOMINATOR (x
));
5654 else if (SCM_REALP (y
))
5655 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
5656 else if (SCM_COMPLEXP (y
))
5657 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
5658 SCM_COMPLEX_IMAG (y
));
5659 else if (SCM_FRACTIONP (y
))
5660 /* a/b + c/d = (ad + bc) / bd */
5661 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5662 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
5663 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
5665 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5668 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
5672 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
5674 "Return @math{@var{x}+1}.")
5675 #define FUNC_NAME s_scm_oneplus
5677 return scm_sum (x
, SCM_INUM1
);
5682 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
5683 (SCM x
, SCM y
, SCM rest
),
5684 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
5685 "the sum of all but the first argument are subtracted from the first\n"
5687 #define FUNC_NAME s_scm_i_difference
5689 while (!scm_is_null (rest
))
5690 { x
= scm_difference (x
, y
);
5692 rest
= scm_cdr (rest
);
5694 return scm_difference (x
, y
);
5698 #define s_difference s_scm_i_difference
5699 #define g_difference g_scm_i_difference
5702 scm_difference (SCM x
, SCM y
)
5703 #define FUNC_NAME s_difference
5705 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5708 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
5710 if (SCM_I_INUMP (x
))
5712 scm_t_inum xx
= -SCM_I_INUM (x
);
5713 if (SCM_FIXABLE (xx
))
5714 return SCM_I_MAKINUM (xx
);
5716 return scm_i_inum2big (xx
);
5718 else if (SCM_BIGP (x
))
5719 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
5720 bignum, but negating that gives a fixnum. */
5721 return scm_i_normbig (scm_i_clonebig (x
, 0));
5722 else if (SCM_REALP (x
))
5723 return scm_from_double (-SCM_REAL_VALUE (x
));
5724 else if (SCM_COMPLEXP (x
))
5725 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
5726 -SCM_COMPLEX_IMAG (x
));
5727 else if (SCM_FRACTIONP (x
))
5728 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
5729 SCM_FRACTION_DENOMINATOR (x
));
5731 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
5734 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5736 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5738 scm_t_inum xx
= SCM_I_INUM (x
);
5739 scm_t_inum yy
= SCM_I_INUM (y
);
5740 scm_t_inum z
= xx
- yy
;
5741 if (SCM_FIXABLE (z
))
5742 return SCM_I_MAKINUM (z
);
5744 return scm_i_inum2big (z
);
5746 else if (SCM_BIGP (y
))
5748 /* inum-x - big-y */
5749 scm_t_inum xx
= SCM_I_INUM (x
);
5753 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
5754 bignum, but negating that gives a fixnum. */
5755 return scm_i_normbig (scm_i_clonebig (y
, 0));
5759 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5760 SCM result
= scm_i_mkbig ();
5763 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
5766 /* x - y == -(y + -x) */
5767 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
5768 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5770 scm_remember_upto_here_1 (y
);
5772 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
5773 /* we know the result will have to be a bignum */
5776 return scm_i_normbig (result
);
5779 else if (SCM_REALP (y
))
5781 scm_t_inum xx
= SCM_I_INUM (x
);
5784 * We need to handle x == exact 0
5785 * specially because R6RS states that:
5786 * (- 0.0) ==> -0.0 and
5787 * (- 0.0 0.0) ==> 0.0
5788 * and the scheme compiler changes
5789 * (- 0.0) into (- 0 0.0)
5790 * So we need to treat (- 0 0.0) like (- 0.0).
5791 * At the C level, (-x) is different than (0.0 - x).
5792 * (0.0 - 0.0) ==> 0.0, but (- 0.0) ==> -0.0.
5795 return scm_from_double (- SCM_REAL_VALUE (y
));
5797 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
5799 else if (SCM_COMPLEXP (y
))
5801 scm_t_inum xx
= SCM_I_INUM (x
);
5803 /* We need to handle x == exact 0 specially.
5804 See the comment above (for SCM_REALP (y)) */
5806 return scm_c_make_rectangular (- SCM_COMPLEX_REAL (y
),
5807 - SCM_COMPLEX_IMAG (y
));
5809 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
5810 - SCM_COMPLEX_IMAG (y
));
5812 else if (SCM_FRACTIONP (y
))
5813 /* a - b/c = (ac - b) / c */
5814 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5815 SCM_FRACTION_NUMERATOR (y
)),
5816 SCM_FRACTION_DENOMINATOR (y
));
5818 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5820 else if (SCM_BIGP (x
))
5822 if (SCM_I_INUMP (y
))
5824 /* big-x - inum-y */
5825 scm_t_inum yy
= SCM_I_INUM (y
);
5826 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5828 scm_remember_upto_here_1 (x
);
5830 return (SCM_FIXABLE (-yy
) ?
5831 SCM_I_MAKINUM (-yy
) : scm_from_inum (-yy
));
5834 SCM result
= scm_i_mkbig ();
5837 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
5839 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
5840 scm_remember_upto_here_1 (x
);
5842 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
5843 /* we know the result will have to be a bignum */
5846 return scm_i_normbig (result
);
5849 else if (SCM_BIGP (y
))
5851 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5852 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5853 SCM result
= scm_i_mkbig ();
5854 mpz_sub (SCM_I_BIG_MPZ (result
),
5857 scm_remember_upto_here_2 (x
, y
);
5858 /* we know the result will have to be a bignum */
5859 if ((sgn_x
== 1) && (sgn_y
== -1))
5861 if ((sgn_x
== -1) && (sgn_y
== 1))
5863 return scm_i_normbig (result
);
5865 else if (SCM_REALP (y
))
5867 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
5868 scm_remember_upto_here_1 (x
);
5869 return scm_from_double (result
);
5871 else if (SCM_COMPLEXP (y
))
5873 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
5874 - SCM_COMPLEX_REAL (y
));
5875 scm_remember_upto_here_1 (x
);
5876 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
5878 else if (SCM_FRACTIONP (y
))
5879 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5880 SCM_FRACTION_NUMERATOR (y
)),
5881 SCM_FRACTION_DENOMINATOR (y
));
5882 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5884 else if (SCM_REALP (x
))
5886 if (SCM_I_INUMP (y
))
5887 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
5888 else if (SCM_BIGP (y
))
5890 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
5891 scm_remember_upto_here_1 (x
);
5892 return scm_from_double (result
);
5894 else if (SCM_REALP (y
))
5895 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
5896 else if (SCM_COMPLEXP (y
))
5897 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
5898 -SCM_COMPLEX_IMAG (y
));
5899 else if (SCM_FRACTIONP (y
))
5900 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
5902 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5904 else if (SCM_COMPLEXP (x
))
5906 if (SCM_I_INUMP (y
))
5907 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
5908 SCM_COMPLEX_IMAG (x
));
5909 else if (SCM_BIGP (y
))
5911 double real_part
= (SCM_COMPLEX_REAL (x
)
5912 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
5913 scm_remember_upto_here_1 (x
);
5914 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
5916 else if (SCM_REALP (y
))
5917 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
5918 SCM_COMPLEX_IMAG (x
));
5919 else if (SCM_COMPLEXP (y
))
5920 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
5921 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
5922 else if (SCM_FRACTIONP (y
))
5923 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
5924 SCM_COMPLEX_IMAG (x
));
5926 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5928 else if (SCM_FRACTIONP (x
))
5930 if (SCM_I_INUMP (y
))
5931 /* a/b - c = (a - cb) / b */
5932 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
5933 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
5934 SCM_FRACTION_DENOMINATOR (x
));
5935 else if (SCM_BIGP (y
))
5936 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
5937 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
5938 SCM_FRACTION_DENOMINATOR (x
));
5939 else if (SCM_REALP (y
))
5940 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
5941 else if (SCM_COMPLEXP (y
))
5942 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
5943 -SCM_COMPLEX_IMAG (y
));
5944 else if (SCM_FRACTIONP (y
))
5945 /* a/b - c/d = (ad - bc) / bd */
5946 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5947 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
5948 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
5950 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5953 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
5958 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
5960 "Return @math{@var{x}-1}.")
5961 #define FUNC_NAME s_scm_oneminus
5963 return scm_difference (x
, SCM_INUM1
);
5968 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
5969 (SCM x
, SCM y
, SCM rest
),
5970 "Return the product of all arguments. If called without arguments,\n"
5972 #define FUNC_NAME s_scm_i_product
5974 while (!scm_is_null (rest
))
5975 { x
= scm_product (x
, y
);
5977 rest
= scm_cdr (rest
);
5979 return scm_product (x
, y
);
5983 #define s_product s_scm_i_product
5984 #define g_product g_scm_i_product
5987 scm_product (SCM x
, SCM y
)
5989 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5992 return SCM_I_MAKINUM (1L);
5993 else if (SCM_NUMBERP (x
))
5996 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
5999 if (SCM_LIKELY (SCM_I_INUMP (x
)))
6004 xx
= SCM_I_INUM (x
);
6009 /* exact1 is the universal multiplicative identity */
6013 /* exact0 times a fixnum is exact0: optimize this case */
6014 if (SCM_LIKELY (SCM_I_INUMP (y
)))
6016 /* if the other argument is inexact, the result is inexact,
6017 and we must do the multiplication in order to handle
6018 infinities and NaNs properly. */
6019 else if (SCM_REALP (y
))
6020 return scm_from_double (0.0 * SCM_REAL_VALUE (y
));
6021 else if (SCM_COMPLEXP (y
))
6022 return scm_c_make_rectangular (0.0 * SCM_COMPLEX_REAL (y
),
6023 0.0 * SCM_COMPLEX_IMAG (y
));
6024 /* we've already handled inexact numbers,
6025 so y must be exact, and we return exact0 */
6026 else if (SCM_NUMP (y
))
6029 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6033 * This case is important for more than just optimization.
6034 * It handles the case of negating
6035 * (+ 1 most-positive-fixnum) aka (- most-negative-fixnum),
6036 * which is a bignum that must be changed back into a fixnum.
6037 * Failure to do so will cause the following to return #f:
6038 * (= most-negative-fixnum (* -1 (- most-negative-fixnum)))
6040 return scm_difference(y
, SCM_UNDEFINED
);
6044 if (SCM_LIKELY (SCM_I_INUMP (y
)))
6046 scm_t_inum yy
= SCM_I_INUM (y
);
6047 scm_t_inum kk
= xx
* yy
;
6048 SCM k
= SCM_I_MAKINUM (kk
);
6049 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
6053 SCM result
= scm_i_inum2big (xx
);
6054 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
6055 return scm_i_normbig (result
);
6058 else if (SCM_BIGP (y
))
6060 SCM result
= scm_i_mkbig ();
6061 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
6062 scm_remember_upto_here_1 (y
);
6065 else if (SCM_REALP (y
))
6066 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
6067 else if (SCM_COMPLEXP (y
))
6068 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
6069 xx
* SCM_COMPLEX_IMAG (y
));
6070 else if (SCM_FRACTIONP (y
))
6071 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
6072 SCM_FRACTION_DENOMINATOR (y
));
6074 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6076 else if (SCM_BIGP (x
))
6078 if (SCM_I_INUMP (y
))
6083 else if (SCM_BIGP (y
))
6085 SCM result
= scm_i_mkbig ();
6086 mpz_mul (SCM_I_BIG_MPZ (result
),
6089 scm_remember_upto_here_2 (x
, y
);
6092 else if (SCM_REALP (y
))
6094 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
6095 scm_remember_upto_here_1 (x
);
6096 return scm_from_double (result
);
6098 else if (SCM_COMPLEXP (y
))
6100 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
6101 scm_remember_upto_here_1 (x
);
6102 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
6103 z
* SCM_COMPLEX_IMAG (y
));
6105 else if (SCM_FRACTIONP (y
))
6106 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
6107 SCM_FRACTION_DENOMINATOR (y
));
6109 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6111 else if (SCM_REALP (x
))
6113 if (SCM_I_INUMP (y
))
6118 else if (SCM_BIGP (y
))
6120 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
6121 scm_remember_upto_here_1 (y
);
6122 return scm_from_double (result
);
6124 else if (SCM_REALP (y
))
6125 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
6126 else if (SCM_COMPLEXP (y
))
6127 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
6128 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
6129 else if (SCM_FRACTIONP (y
))
6130 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
6132 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6134 else if (SCM_COMPLEXP (x
))
6136 if (SCM_I_INUMP (y
))
6141 else if (SCM_BIGP (y
))
6143 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6144 scm_remember_upto_here_1 (y
);
6145 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
6146 z
* SCM_COMPLEX_IMAG (x
));
6148 else if (SCM_REALP (y
))
6149 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
6150 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
6151 else if (SCM_COMPLEXP (y
))
6153 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
6154 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
6155 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
6156 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
6158 else if (SCM_FRACTIONP (y
))
6160 double yy
= scm_i_fraction2double (y
);
6161 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
6162 yy
* SCM_COMPLEX_IMAG (x
));
6165 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6167 else if (SCM_FRACTIONP (x
))
6169 if (SCM_I_INUMP (y
))
6170 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
6171 SCM_FRACTION_DENOMINATOR (x
));
6172 else if (SCM_BIGP (y
))
6173 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
6174 SCM_FRACTION_DENOMINATOR (x
));
6175 else if (SCM_REALP (y
))
6176 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
6177 else if (SCM_COMPLEXP (y
))
6179 double xx
= scm_i_fraction2double (x
);
6180 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
6181 xx
* SCM_COMPLEX_IMAG (y
));
6183 else if (SCM_FRACTIONP (y
))
6184 /* a/b * c/d = ac / bd */
6185 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
6186 SCM_FRACTION_NUMERATOR (y
)),
6187 scm_product (SCM_FRACTION_DENOMINATOR (x
),
6188 SCM_FRACTION_DENOMINATOR (y
)));
6190 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6193 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
6196 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
6197 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
6198 #define ALLOW_DIVIDE_BY_ZERO
6199 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
6202 /* The code below for complex division is adapted from the GNU
6203 libstdc++, which adapted it from f2c's libF77, and is subject to
6206 /****************************************************************
6207 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
6209 Permission to use, copy, modify, and distribute this software
6210 and its documentation for any purpose and without fee is hereby
6211 granted, provided that the above copyright notice appear in all
6212 copies and that both that the copyright notice and this
6213 permission notice and warranty disclaimer appear in supporting
6214 documentation, and that the names of AT&T Bell Laboratories or
6215 Bellcore or any of their entities not be used in advertising or
6216 publicity pertaining to distribution of the software without
6217 specific, written prior permission.
6219 AT&T and Bellcore disclaim all warranties with regard to this
6220 software, including all implied warranties of merchantability
6221 and fitness. In no event shall AT&T or Bellcore be liable for
6222 any special, indirect or consequential damages or any damages
6223 whatsoever resulting from loss of use, data or profits, whether
6224 in an action of contract, negligence or other tortious action,
6225 arising out of or in connection with the use or performance of
6227 ****************************************************************/
6229 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
6230 (SCM x
, SCM y
, SCM rest
),
6231 "Divide the first argument by the product of the remaining\n"
6232 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
6234 #define FUNC_NAME s_scm_i_divide
6236 while (!scm_is_null (rest
))
6237 { x
= scm_divide (x
, y
);
6239 rest
= scm_cdr (rest
);
6241 return scm_divide (x
, y
);
6245 #define s_divide s_scm_i_divide
6246 #define g_divide g_scm_i_divide
6249 do_divide (SCM x
, SCM y
, int inexact
)
6250 #define FUNC_NAME s_divide
6254 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
6257 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
6258 else if (SCM_I_INUMP (x
))
6260 scm_t_inum xx
= SCM_I_INUM (x
);
6261 if (xx
== 1 || xx
== -1)
6263 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6265 scm_num_overflow (s_divide
);
6270 return scm_from_double (1.0 / (double) xx
);
6271 else return scm_i_make_ratio (SCM_INUM1
, x
);
6274 else if (SCM_BIGP (x
))
6277 return scm_from_double (1.0 / scm_i_big2dbl (x
));
6278 else return scm_i_make_ratio (SCM_INUM1
, x
);
6280 else if (SCM_REALP (x
))
6282 double xx
= SCM_REAL_VALUE (x
);
6283 #ifndef ALLOW_DIVIDE_BY_ZERO
6285 scm_num_overflow (s_divide
);
6288 return scm_from_double (1.0 / xx
);
6290 else if (SCM_COMPLEXP (x
))
6292 double r
= SCM_COMPLEX_REAL (x
);
6293 double i
= SCM_COMPLEX_IMAG (x
);
6294 if (fabs(r
) <= fabs(i
))
6297 double d
= i
* (1.0 + t
* t
);
6298 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
6303 double d
= r
* (1.0 + t
* t
);
6304 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
6307 else if (SCM_FRACTIONP (x
))
6308 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
6309 SCM_FRACTION_NUMERATOR (x
));
6311 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
6314 if (SCM_LIKELY (SCM_I_INUMP (x
)))
6316 scm_t_inum xx
= SCM_I_INUM (x
);
6317 if (SCM_LIKELY (SCM_I_INUMP (y
)))
6319 scm_t_inum yy
= SCM_I_INUM (y
);
6322 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6323 scm_num_overflow (s_divide
);
6325 return scm_from_double ((double) xx
/ (double) yy
);
6328 else if (xx
% yy
!= 0)
6331 return scm_from_double ((double) xx
/ (double) yy
);
6332 else return scm_i_make_ratio (x
, y
);
6336 scm_t_inum z
= xx
/ yy
;
6337 if (SCM_FIXABLE (z
))
6338 return SCM_I_MAKINUM (z
);
6340 return scm_i_inum2big (z
);
6343 else if (SCM_BIGP (y
))
6346 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
6347 else return scm_i_make_ratio (x
, y
);
6349 else if (SCM_REALP (y
))
6351 double yy
= SCM_REAL_VALUE (y
);
6352 #ifndef ALLOW_DIVIDE_BY_ZERO
6354 scm_num_overflow (s_divide
);
6357 return scm_from_double ((double) xx
/ yy
);
6359 else if (SCM_COMPLEXP (y
))
6362 complex_div
: /* y _must_ be a complex number */
6364 double r
= SCM_COMPLEX_REAL (y
);
6365 double i
= SCM_COMPLEX_IMAG (y
);
6366 if (fabs(r
) <= fabs(i
))
6369 double d
= i
* (1.0 + t
* t
);
6370 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
6375 double d
= r
* (1.0 + t
* t
);
6376 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
6380 else if (SCM_FRACTIONP (y
))
6381 /* a / b/c = ac / b */
6382 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
6383 SCM_FRACTION_NUMERATOR (y
));
6385 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6387 else if (SCM_BIGP (x
))
6389 if (SCM_I_INUMP (y
))
6391 scm_t_inum yy
= SCM_I_INUM (y
);
6394 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6395 scm_num_overflow (s_divide
);
6397 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
6398 scm_remember_upto_here_1 (x
);
6399 return (sgn
== 0) ? scm_nan () : scm_inf ();
6406 /* FIXME: HMM, what are the relative performance issues here?
6407 We need to test. Is it faster on average to test
6408 divisible_p, then perform whichever operation, or is it
6409 faster to perform the integer div opportunistically and
6410 switch to real if there's a remainder? For now we take the
6411 middle ground: test, then if divisible, use the faster div
6414 scm_t_inum abs_yy
= yy
< 0 ? -yy
: yy
;
6415 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
6419 SCM result
= scm_i_mkbig ();
6420 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
6421 scm_remember_upto_here_1 (x
);
6423 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
6424 return scm_i_normbig (result
);
6429 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
6430 else return scm_i_make_ratio (x
, y
);
6434 else if (SCM_BIGP (y
))
6439 /* It's easily possible for the ratio x/y to fit a double
6440 but one or both x and y be too big to fit a double,
6441 hence the use of mpq_get_d rather than converting and
6444 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
6445 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
6446 return scm_from_double (mpq_get_d (q
));
6450 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
6454 SCM result
= scm_i_mkbig ();
6455 mpz_divexact (SCM_I_BIG_MPZ (result
),
6458 scm_remember_upto_here_2 (x
, y
);
6459 return scm_i_normbig (result
);
6462 return scm_i_make_ratio (x
, y
);
6465 else if (SCM_REALP (y
))
6467 double yy
= SCM_REAL_VALUE (y
);
6468 #ifndef ALLOW_DIVIDE_BY_ZERO
6470 scm_num_overflow (s_divide
);
6473 return scm_from_double (scm_i_big2dbl (x
) / yy
);
6475 else if (SCM_COMPLEXP (y
))
6477 a
= scm_i_big2dbl (x
);
6480 else if (SCM_FRACTIONP (y
))
6481 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
6482 SCM_FRACTION_NUMERATOR (y
));
6484 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6486 else if (SCM_REALP (x
))
6488 double rx
= SCM_REAL_VALUE (x
);
6489 if (SCM_I_INUMP (y
))
6491 scm_t_inum yy
= SCM_I_INUM (y
);
6492 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6494 scm_num_overflow (s_divide
);
6497 return scm_from_double (rx
/ (double) yy
);
6499 else if (SCM_BIGP (y
))
6501 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6502 scm_remember_upto_here_1 (y
);
6503 return scm_from_double (rx
/ dby
);
6505 else if (SCM_REALP (y
))
6507 double yy
= SCM_REAL_VALUE (y
);
6508 #ifndef ALLOW_DIVIDE_BY_ZERO
6510 scm_num_overflow (s_divide
);
6513 return scm_from_double (rx
/ yy
);
6515 else if (SCM_COMPLEXP (y
))
6520 else if (SCM_FRACTIONP (y
))
6521 return scm_from_double (rx
/ scm_i_fraction2double (y
));
6523 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6525 else if (SCM_COMPLEXP (x
))
6527 double rx
= SCM_COMPLEX_REAL (x
);
6528 double ix
= SCM_COMPLEX_IMAG (x
);
6529 if (SCM_I_INUMP (y
))
6531 scm_t_inum yy
= SCM_I_INUM (y
);
6532 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6534 scm_num_overflow (s_divide
);
6539 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
6542 else if (SCM_BIGP (y
))
6544 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6545 scm_remember_upto_here_1 (y
);
6546 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
6548 else if (SCM_REALP (y
))
6550 double yy
= SCM_REAL_VALUE (y
);
6551 #ifndef ALLOW_DIVIDE_BY_ZERO
6553 scm_num_overflow (s_divide
);
6556 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
6558 else if (SCM_COMPLEXP (y
))
6560 double ry
= SCM_COMPLEX_REAL (y
);
6561 double iy
= SCM_COMPLEX_IMAG (y
);
6562 if (fabs(ry
) <= fabs(iy
))
6565 double d
= iy
* (1.0 + t
* t
);
6566 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
6571 double d
= ry
* (1.0 + t
* t
);
6572 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
6575 else if (SCM_FRACTIONP (y
))
6577 double yy
= scm_i_fraction2double (y
);
6578 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
6581 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6583 else if (SCM_FRACTIONP (x
))
6585 if (SCM_I_INUMP (y
))
6587 scm_t_inum yy
= SCM_I_INUM (y
);
6588 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6590 scm_num_overflow (s_divide
);
6593 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
6594 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
6596 else if (SCM_BIGP (y
))
6598 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
6599 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
6601 else if (SCM_REALP (y
))
6603 double yy
= SCM_REAL_VALUE (y
);
6604 #ifndef ALLOW_DIVIDE_BY_ZERO
6606 scm_num_overflow (s_divide
);
6609 return scm_from_double (scm_i_fraction2double (x
) / yy
);
6611 else if (SCM_COMPLEXP (y
))
6613 a
= scm_i_fraction2double (x
);
6616 else if (SCM_FRACTIONP (y
))
6617 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
6618 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
6620 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6623 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
6627 scm_divide (SCM x
, SCM y
)
6629 return do_divide (x
, y
, 0);
6632 static SCM
scm_divide2real (SCM x
, SCM y
)
6634 return do_divide (x
, y
, 1);
6640 scm_c_truncate (double x
)
6651 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
6652 half-way case (ie. when x is an integer plus 0.5) going upwards.
6653 Then half-way cases are identified and adjusted down if the
6654 round-upwards didn't give the desired even integer.
6656 "plus_half == result" identifies a half-way case. If plus_half, which is
6657 x + 0.5, is an integer then x must be an integer plus 0.5.
6659 An odd "result" value is identified with result/2 != floor(result/2).
6660 This is done with plus_half, since that value is ready for use sooner in
6661 a pipelined cpu, and we're already requiring plus_half == result.
6663 Note however that we need to be careful when x is big and already an
6664 integer. In that case "x+0.5" may round to an adjacent integer, causing
6665 us to return such a value, incorrectly. For instance if the hardware is
6666 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
6667 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
6668 returned. Or if the hardware is in round-upwards mode, then other bigger
6669 values like say x == 2^128 will see x+0.5 rounding up to the next higher
6670 representable value, 2^128+2^76 (or whatever), again incorrect.
6672 These bad roundings of x+0.5 are avoided by testing at the start whether
6673 x is already an integer. If it is then clearly that's the desired result
6674 already. And if it's not then the exponent must be small enough to allow
6675 an 0.5 to be represented, and hence added without a bad rounding. */
6678 scm_c_round (double x
)
6680 double plus_half
, result
;
6685 plus_half
= x
+ 0.5;
6686 result
= floor (plus_half
);
6687 /* Adjust so that the rounding is towards even. */
6688 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
6693 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
6695 "Round the number @var{x} towards zero.")
6696 #define FUNC_NAME s_scm_truncate_number
6698 if (scm_is_false (scm_negative_p (x
)))
6699 return scm_floor (x
);
6701 return scm_ceiling (x
);
6705 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
6707 "Round the number @var{x} towards the nearest integer. "
6708 "When it is exactly halfway between two integers, "
6709 "round towards the even one.")
6710 #define FUNC_NAME s_scm_round_number
6712 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6714 else if (SCM_REALP (x
))
6715 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
6718 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
6719 single quotient+remainder division then examining to see which way
6720 the rounding should go. */
6721 SCM plus_half
= scm_sum (x
, exactly_one_half
);
6722 SCM result
= scm_floor (plus_half
);
6723 /* Adjust so that the rounding is towards even. */
6724 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
6725 && scm_is_true (scm_odd_p (result
)))
6726 return scm_difference (result
, SCM_INUM1
);
6733 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
6735 "Round the number @var{x} towards minus infinity.")
6736 #define FUNC_NAME s_scm_floor
6738 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6740 else if (SCM_REALP (x
))
6741 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
6742 else if (SCM_FRACTIONP (x
))
6744 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
6745 SCM_FRACTION_DENOMINATOR (x
));
6746 if (scm_is_false (scm_negative_p (x
)))
6748 /* For positive x, rounding towards zero is correct. */
6753 /* For negative x, we need to return q-1 unless x is an
6754 integer. But fractions are never integer, per our
6756 return scm_difference (q
, SCM_INUM1
);
6760 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
6764 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
6766 "Round the number @var{x} towards infinity.")
6767 #define FUNC_NAME s_scm_ceiling
6769 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6771 else if (SCM_REALP (x
))
6772 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
6773 else if (SCM_FRACTIONP (x
))
6775 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
6776 SCM_FRACTION_DENOMINATOR (x
));
6777 if (scm_is_false (scm_positive_p (x
)))
6779 /* For negative x, rounding towards zero is correct. */
6784 /* For positive x, we need to return q+1 unless x is an
6785 integer. But fractions are never integer, per our
6787 return scm_sum (q
, SCM_INUM1
);
6791 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
6795 SCM_PRIMITIVE_GENERIC (scm_expt
, "expt", 2, 0, 0,
6797 "Return @var{x} raised to the power of @var{y}.")
6798 #define FUNC_NAME s_scm_expt
6800 if (scm_is_integer (y
))
6802 if (scm_is_true (scm_exact_p (y
)))
6803 return scm_integer_expt (x
, y
);
6806 /* Here we handle the case where the exponent is an inexact
6807 integer. We make the exponent exact in order to use
6808 scm_integer_expt, and thus avoid the spurious imaginary
6809 parts that may result from round-off errors in the general
6810 e^(y log x) method below (for example when squaring a large
6811 negative number). In this case, we must return an inexact
6812 result for correctness. We also make the base inexact so
6813 that scm_integer_expt will use fast inexact arithmetic
6814 internally. Note that making the base inexact is not
6815 sufficient to guarantee an inexact result, because
6816 scm_integer_expt will return an exact 1 when the exponent
6817 is 0, even if the base is inexact. */
6818 return scm_exact_to_inexact
6819 (scm_integer_expt (scm_exact_to_inexact (x
),
6820 scm_inexact_to_exact (y
)));
6823 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
6825 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
6827 else if (scm_is_complex (x
) && scm_is_complex (y
))
6828 return scm_exp (scm_product (scm_log (x
), y
));
6829 else if (scm_is_complex (x
))
6830 SCM_WTA_DISPATCH_2 (g_scm_expt
, x
, y
, SCM_ARG2
, s_scm_expt
);
6832 SCM_WTA_DISPATCH_2 (g_scm_expt
, x
, y
, SCM_ARG1
, s_scm_expt
);
6836 /* sin/cos/tan/asin/acos/atan
6837 sinh/cosh/tanh/asinh/acosh/atanh
6838 Derived from "Transcen.scm", Complex trancendental functions for SCM.
6839 Written by Jerry D. Hedden, (C) FSF.
6840 See the file `COPYING' for terms applying to this program. */
6842 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
6844 "Compute the sine of @var{z}.")
6845 #define FUNC_NAME s_scm_sin
6847 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6848 return z
; /* sin(exact0) = exact0 */
6849 else if (scm_is_real (z
))
6850 return scm_from_double (sin (scm_to_double (z
)));
6851 else if (SCM_COMPLEXP (z
))
6853 x
= SCM_COMPLEX_REAL (z
);
6854 y
= SCM_COMPLEX_IMAG (z
);
6855 return scm_c_make_rectangular (sin (x
) * cosh (y
),
6856 cos (x
) * sinh (y
));
6859 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
6863 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
6865 "Compute the cosine of @var{z}.")
6866 #define FUNC_NAME s_scm_cos
6868 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6869 return SCM_INUM1
; /* cos(exact0) = exact1 */
6870 else if (scm_is_real (z
))
6871 return scm_from_double (cos (scm_to_double (z
)));
6872 else if (SCM_COMPLEXP (z
))
6874 x
= SCM_COMPLEX_REAL (z
);
6875 y
= SCM_COMPLEX_IMAG (z
);
6876 return scm_c_make_rectangular (cos (x
) * cosh (y
),
6877 -sin (x
) * sinh (y
));
6880 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
6884 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
6886 "Compute the tangent of @var{z}.")
6887 #define FUNC_NAME s_scm_tan
6889 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6890 return z
; /* tan(exact0) = exact0 */
6891 else if (scm_is_real (z
))
6892 return scm_from_double (tan (scm_to_double (z
)));
6893 else if (SCM_COMPLEXP (z
))
6895 x
= 2.0 * SCM_COMPLEX_REAL (z
);
6896 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
6897 w
= cos (x
) + cosh (y
);
6898 #ifndef ALLOW_DIVIDE_BY_ZERO
6900 scm_num_overflow (s_scm_tan
);
6902 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
6905 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
6909 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
6911 "Compute the hyperbolic sine of @var{z}.")
6912 #define FUNC_NAME s_scm_sinh
6914 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6915 return z
; /* sinh(exact0) = exact0 */
6916 else if (scm_is_real (z
))
6917 return scm_from_double (sinh (scm_to_double (z
)));
6918 else if (SCM_COMPLEXP (z
))
6920 x
= SCM_COMPLEX_REAL (z
);
6921 y
= SCM_COMPLEX_IMAG (z
);
6922 return scm_c_make_rectangular (sinh (x
) * cos (y
),
6923 cosh (x
) * sin (y
));
6926 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
6930 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
6932 "Compute the hyperbolic cosine of @var{z}.")
6933 #define FUNC_NAME s_scm_cosh
6935 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6936 return SCM_INUM1
; /* cosh(exact0) = exact1 */
6937 else if (scm_is_real (z
))
6938 return scm_from_double (cosh (scm_to_double (z
)));
6939 else if (SCM_COMPLEXP (z
))
6941 x
= SCM_COMPLEX_REAL (z
);
6942 y
= SCM_COMPLEX_IMAG (z
);
6943 return scm_c_make_rectangular (cosh (x
) * cos (y
),
6944 sinh (x
) * sin (y
));
6947 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
6951 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
6953 "Compute the hyperbolic tangent of @var{z}.")
6954 #define FUNC_NAME s_scm_tanh
6956 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6957 return z
; /* tanh(exact0) = exact0 */
6958 else if (scm_is_real (z
))
6959 return scm_from_double (tanh (scm_to_double (z
)));
6960 else if (SCM_COMPLEXP (z
))
6962 x
= 2.0 * SCM_COMPLEX_REAL (z
);
6963 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
6964 w
= cosh (x
) + cos (y
);
6965 #ifndef ALLOW_DIVIDE_BY_ZERO
6967 scm_num_overflow (s_scm_tanh
);
6969 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
6972 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
6976 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
6978 "Compute the arc sine of @var{z}.")
6979 #define FUNC_NAME s_scm_asin
6981 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6982 return z
; /* asin(exact0) = exact0 */
6983 else if (scm_is_real (z
))
6985 double w
= scm_to_double (z
);
6986 if (w
>= -1.0 && w
<= 1.0)
6987 return scm_from_double (asin (w
));
6989 return 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_product (scm_c_make_rectangular (0, -1),
6997 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
7000 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
7004 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
7006 "Compute the arc cosine of @var{z}.")
7007 #define FUNC_NAME s_scm_acos
7009 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM1
)))
7010 return SCM_INUM0
; /* acos(exact1) = exact0 */
7011 else if (scm_is_real (z
))
7013 double w
= scm_to_double (z
);
7014 if (w
>= -1.0 && w
<= 1.0)
7015 return scm_from_double (acos (w
));
7017 return scm_sum (scm_from_double (acos (0.0)),
7018 scm_product (scm_c_make_rectangular (0, 1),
7019 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
7021 else if (SCM_COMPLEXP (z
))
7023 x
= SCM_COMPLEX_REAL (z
);
7024 y
= SCM_COMPLEX_IMAG (z
);
7025 return scm_sum (scm_from_double (acos (0.0)),
7026 scm_product (scm_c_make_rectangular (0, 1),
7027 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
7030 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
7034 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
7036 "With one argument, compute the arc tangent of @var{z}.\n"
7037 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
7038 "using the sign of @var{z} and @var{y} to determine the quadrant.")
7039 #define FUNC_NAME s_scm_atan
7043 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
7044 return z
; /* atan(exact0) = exact0 */
7045 else if (scm_is_real (z
))
7046 return scm_from_double (atan (scm_to_double (z
)));
7047 else if (SCM_COMPLEXP (z
))
7050 v
= SCM_COMPLEX_REAL (z
);
7051 w
= SCM_COMPLEX_IMAG (z
);
7052 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
7053 scm_c_make_rectangular (v
, w
+ 1.0))),
7054 scm_c_make_rectangular (0, 2));
7057 SCM_WTA_DISPATCH_1 (g_scm_atan
, z
, SCM_ARG1
, s_scm_atan
);
7059 else if (scm_is_real (z
))
7061 if (scm_is_real (y
))
7062 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
7064 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
7067 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
7071 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
7073 "Compute the inverse hyperbolic sine of @var{z}.")
7074 #define FUNC_NAME s_scm_sys_asinh
7076 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
7077 return z
; /* asinh(exact0) = exact0 */
7078 else if (scm_is_real (z
))
7079 return scm_from_double (asinh (scm_to_double (z
)));
7080 else if (scm_is_number (z
))
7081 return scm_log (scm_sum (z
,
7082 scm_sqrt (scm_sum (scm_product (z
, z
),
7085 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
7089 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
7091 "Compute the inverse hyperbolic cosine of @var{z}.")
7092 #define FUNC_NAME s_scm_sys_acosh
7094 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM1
)))
7095 return SCM_INUM0
; /* acosh(exact1) = exact0 */
7096 else if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
7097 return scm_from_double (acosh (scm_to_double (z
)));
7098 else if (scm_is_number (z
))
7099 return scm_log (scm_sum (z
,
7100 scm_sqrt (scm_difference (scm_product (z
, z
),
7103 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
7107 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
7109 "Compute the inverse hyperbolic tangent of @var{z}.")
7110 #define FUNC_NAME s_scm_sys_atanh
7112 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
7113 return z
; /* atanh(exact0) = exact0 */
7114 else if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
7115 return scm_from_double (atanh (scm_to_double (z
)));
7116 else if (scm_is_number (z
))
7117 return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1
, z
),
7118 scm_difference (SCM_INUM1
, z
))),
7121 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
7126 scm_c_make_rectangular (double re
, double im
)
7130 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
7132 SCM_SET_CELL_TYPE (z
, scm_tc16_complex
);
7133 SCM_COMPLEX_REAL (z
) = re
;
7134 SCM_COMPLEX_IMAG (z
) = im
;
7138 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
7139 (SCM real_part
, SCM imaginary_part
),
7140 "Return a complex number constructed of the given @var{real-part} "
7141 "and @var{imaginary-part} parts.")
7142 #define FUNC_NAME s_scm_make_rectangular
7144 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
7145 SCM_ARG1
, FUNC_NAME
, "real");
7146 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
7147 SCM_ARG2
, FUNC_NAME
, "real");
7149 /* Return a real if and only if the imaginary_part is an _exact_ 0 */
7150 if (scm_is_eq (imaginary_part
, SCM_INUM0
))
7153 return scm_c_make_rectangular (scm_to_double (real_part
),
7154 scm_to_double (imaginary_part
));
7159 scm_c_make_polar (double mag
, double ang
)
7163 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
7164 use it on Glibc-based systems that have it (it's a GNU extension). See
7165 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
7167 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
7168 sincos (ang
, &s
, &c
);
7174 /* If s and c are NaNs, this indicates that the angle is a NaN,
7175 infinite, or perhaps simply too large to determine its value
7176 mod 2*pi. However, we know something that the floating-point
7177 implementation doesn't know: We know that s and c are finite.
7178 Therefore, if the magnitude is zero, return a complex zero.
7180 The reason we check for the NaNs instead of using this case
7181 whenever mag == 0.0 is because when the angle is known, we'd
7182 like to return the correct kind of non-real complex zero:
7183 +0.0+0.0i, -0.0+0.0i, -0.0-0.0i, or +0.0-0.0i, depending
7184 on which quadrant the angle is in.
7186 if (SCM_UNLIKELY (isnan(s
)) && isnan(c
) && (mag
== 0.0))
7187 return scm_c_make_rectangular (0.0, 0.0);
7189 return scm_c_make_rectangular (mag
* c
, mag
* s
);
7192 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
7194 "Return the complex number @var{mag} * e^(i * @var{ang}).")
7195 #define FUNC_NAME s_scm_make_polar
7197 SCM_ASSERT_TYPE (scm_is_real (mag
), mag
, SCM_ARG1
, FUNC_NAME
, "real");
7198 SCM_ASSERT_TYPE (scm_is_real (ang
), ang
, SCM_ARG2
, FUNC_NAME
, "real");
7200 /* If mag is exact0, return exact0 */
7201 if (scm_is_eq (mag
, SCM_INUM0
))
7203 /* Return a real if ang is exact0 */
7204 else if (scm_is_eq (ang
, SCM_INUM0
))
7207 return scm_c_make_polar (scm_to_double (mag
), scm_to_double (ang
));
7212 SCM_PRIMITIVE_GENERIC (scm_real_part
, "real-part", 1, 0, 0,
7214 "Return the real part of the number @var{z}.")
7215 #define FUNC_NAME s_scm_real_part
7217 if (SCM_COMPLEXP (z
))
7218 return scm_from_double (SCM_COMPLEX_REAL (z
));
7219 else if (SCM_I_INUMP (z
) || SCM_BIGP (z
) || SCM_REALP (z
) || SCM_FRACTIONP (z
))
7222 SCM_WTA_DISPATCH_1 (g_scm_real_part
, z
, SCM_ARG1
, s_scm_real_part
);
7227 SCM_PRIMITIVE_GENERIC (scm_imag_part
, "imag-part", 1, 0, 0,
7229 "Return the imaginary part of the number @var{z}.")
7230 #define FUNC_NAME s_scm_imag_part
7232 if (SCM_COMPLEXP (z
))
7233 return scm_from_double (SCM_COMPLEX_IMAG (z
));
7234 else if (SCM_I_INUMP (z
) || SCM_REALP (z
) || SCM_BIGP (z
) || SCM_FRACTIONP (z
))
7237 SCM_WTA_DISPATCH_1 (g_scm_imag_part
, z
, SCM_ARG1
, s_scm_imag_part
);
7241 SCM_PRIMITIVE_GENERIC (scm_numerator
, "numerator", 1, 0, 0,
7243 "Return the numerator of the number @var{z}.")
7244 #define FUNC_NAME s_scm_numerator
7246 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7248 else if (SCM_FRACTIONP (z
))
7249 return SCM_FRACTION_NUMERATOR (z
);
7250 else if (SCM_REALP (z
))
7251 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
7253 SCM_WTA_DISPATCH_1 (g_scm_numerator
, z
, SCM_ARG1
, s_scm_numerator
);
7258 SCM_PRIMITIVE_GENERIC (scm_denominator
, "denominator", 1, 0, 0,
7260 "Return the denominator of the number @var{z}.")
7261 #define FUNC_NAME s_scm_denominator
7263 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7265 else if (SCM_FRACTIONP (z
))
7266 return SCM_FRACTION_DENOMINATOR (z
);
7267 else if (SCM_REALP (z
))
7268 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
7270 SCM_WTA_DISPATCH_1 (g_scm_denominator
, z
, SCM_ARG1
, s_scm_denominator
);
7275 SCM_PRIMITIVE_GENERIC (scm_magnitude
, "magnitude", 1, 0, 0,
7277 "Return the magnitude of the number @var{z}. This is the same as\n"
7278 "@code{abs} for real arguments, but also allows complex numbers.")
7279 #define FUNC_NAME s_scm_magnitude
7281 if (SCM_I_INUMP (z
))
7283 scm_t_inum zz
= SCM_I_INUM (z
);
7286 else if (SCM_POSFIXABLE (-zz
))
7287 return SCM_I_MAKINUM (-zz
);
7289 return scm_i_inum2big (-zz
);
7291 else if (SCM_BIGP (z
))
7293 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
7294 scm_remember_upto_here_1 (z
);
7296 return scm_i_clonebig (z
, 0);
7300 else if (SCM_REALP (z
))
7301 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
7302 else if (SCM_COMPLEXP (z
))
7303 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
7304 else if (SCM_FRACTIONP (z
))
7306 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
7308 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
7309 SCM_FRACTION_DENOMINATOR (z
));
7312 SCM_WTA_DISPATCH_1 (g_scm_magnitude
, z
, SCM_ARG1
, s_scm_magnitude
);
7317 SCM_PRIMITIVE_GENERIC (scm_angle
, "angle", 1, 0, 0,
7319 "Return the angle of the complex number @var{z}.")
7320 #define FUNC_NAME s_scm_angle
7322 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
7323 flo0 to save allocating a new flonum with scm_from_double each time.
7324 But if atan2 follows the floating point rounding mode, then the value
7325 is not a constant. Maybe it'd be close enough though. */
7326 if (SCM_I_INUMP (z
))
7328 if (SCM_I_INUM (z
) >= 0)
7331 return scm_from_double (atan2 (0.0, -1.0));
7333 else if (SCM_BIGP (z
))
7335 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
7336 scm_remember_upto_here_1 (z
);
7338 return scm_from_double (atan2 (0.0, -1.0));
7342 else if (SCM_REALP (z
))
7344 if (SCM_REAL_VALUE (z
) >= 0)
7347 return scm_from_double (atan2 (0.0, -1.0));
7349 else if (SCM_COMPLEXP (z
))
7350 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
7351 else if (SCM_FRACTIONP (z
))
7353 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
7355 else return scm_from_double (atan2 (0.0, -1.0));
7358 SCM_WTA_DISPATCH_1 (g_scm_angle
, z
, SCM_ARG1
, s_scm_angle
);
7363 SCM_PRIMITIVE_GENERIC (scm_exact_to_inexact
, "exact->inexact", 1, 0, 0,
7365 "Convert the number @var{z} to its inexact representation.\n")
7366 #define FUNC_NAME s_scm_exact_to_inexact
7368 if (SCM_I_INUMP (z
))
7369 return scm_from_double ((double) SCM_I_INUM (z
));
7370 else if (SCM_BIGP (z
))
7371 return scm_from_double (scm_i_big2dbl (z
));
7372 else if (SCM_FRACTIONP (z
))
7373 return scm_from_double (scm_i_fraction2double (z
));
7374 else if (SCM_INEXACTP (z
))
7377 SCM_WTA_DISPATCH_1 (g_scm_exact_to_inexact
, z
, 1, s_scm_exact_to_inexact
);
7382 SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
7384 "Return an exact number that is numerically closest to @var{z}.")
7385 #define FUNC_NAME s_scm_inexact_to_exact
7387 if (SCM_I_INUMP (z
) || SCM_BIGP (z
) || SCM_FRACTIONP (z
))
7394 val
= SCM_REAL_VALUE (z
);
7395 else if (SCM_COMPLEXP (z
) && SCM_COMPLEX_IMAG (z
) == 0.0)
7396 val
= SCM_COMPLEX_REAL (z
);
7398 SCM_WTA_DISPATCH_1 (g_scm_inexact_to_exact
, z
, 1, s_scm_inexact_to_exact
);
7400 if (!SCM_LIKELY (DOUBLE_IS_FINITE (val
)))
7401 SCM_OUT_OF_RANGE (1, z
);
7408 mpq_set_d (frac
, val
);
7409 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
7410 scm_i_mpz2num (mpq_denref (frac
)));
7412 /* When scm_i_make_ratio throws, we leak the memory allocated
7422 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
7424 "Returns the @emph{simplest} rational number differing\n"
7425 "from @var{x} by no more than @var{eps}.\n"
7427 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
7428 "exact result when both its arguments are exact. Thus, you might need\n"
7429 "to use @code{inexact->exact} on the arguments.\n"
7432 "(rationalize (inexact->exact 1.2) 1/100)\n"
7435 #define FUNC_NAME s_scm_rationalize
7437 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
7438 SCM_ASSERT_TYPE (scm_is_real (eps
), eps
, SCM_ARG2
, FUNC_NAME
, "real");
7439 eps
= scm_abs (eps
);
7440 if (scm_is_false (scm_positive_p (eps
)))
7442 /* eps is either zero or a NaN */
7443 if (scm_is_true (scm_nan_p (eps
)))
7445 else if (SCM_INEXACTP (eps
))
7446 return scm_exact_to_inexact (x
);
7450 else if (scm_is_false (scm_finite_p (eps
)))
7452 if (scm_is_true (scm_finite_p (x
)))
7457 else if (scm_is_false (scm_finite_p (x
))) /* checks for both inf and nan */
7459 else if (scm_is_false (scm_less_p (scm_floor (scm_sum (x
, eps
)),
7460 scm_ceiling (scm_difference (x
, eps
)))))
7462 /* There's an integer within range; we want the one closest to zero */
7463 if (scm_is_false (scm_less_p (eps
, scm_abs (x
))))
7465 /* zero is within range */
7466 if (SCM_INEXACTP (x
) || SCM_INEXACTP (eps
))
7471 else if (scm_is_true (scm_positive_p (x
)))
7472 return scm_ceiling (scm_difference (x
, eps
));
7474 return scm_floor (scm_sum (x
, eps
));
7478 /* Use continued fractions to find closest ratio. All
7479 arithmetic is done with exact numbers.
7482 SCM ex
= scm_inexact_to_exact (x
);
7483 SCM int_part
= scm_floor (ex
);
7485 SCM a1
= SCM_INUM0
, a2
= SCM_INUM1
, a
= SCM_INUM0
;
7486 SCM b1
= SCM_INUM1
, b2
= SCM_INUM0
, b
= SCM_INUM0
;
7490 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
7491 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
7493 /* We stop after a million iterations just to be absolutely sure
7494 that we don't go into an infinite loop. The process normally
7495 converges after less than a dozen iterations.
7498 while (++i
< 1000000)
7500 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
7501 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
7502 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
7504 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
7505 eps
))) /* abs(x-a/b) <= eps */
7507 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
7508 if (SCM_INEXACTP (x
) || SCM_INEXACTP (eps
))
7509 return scm_exact_to_inexact (res
);
7513 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
7515 tt
= scm_floor (rx
); /* tt = floor (rx) */
7521 scm_num_overflow (s_scm_rationalize
);
7526 /* conversion functions */
7529 scm_is_integer (SCM val
)
7531 return scm_is_true (scm_integer_p (val
));
7535 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
7537 if (SCM_I_INUMP (val
))
7539 scm_t_signed_bits n
= SCM_I_INUM (val
);
7540 return n
>= min
&& n
<= max
;
7542 else if (SCM_BIGP (val
))
7544 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
7546 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
7548 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
7550 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
7551 return n
>= min
&& n
<= max
;
7561 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
7562 > CHAR_BIT
*sizeof (scm_t_uintmax
))
7565 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
7566 SCM_I_BIG_MPZ (val
));
7568 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
7580 return n
>= min
&& n
<= max
;
7588 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
7590 if (SCM_I_INUMP (val
))
7592 scm_t_signed_bits n
= SCM_I_INUM (val
);
7593 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
7595 else if (SCM_BIGP (val
))
7597 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
7599 else if (max
<= ULONG_MAX
)
7601 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
7603 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
7604 return n
>= min
&& n
<= max
;
7614 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
7617 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
7618 > CHAR_BIT
*sizeof (scm_t_uintmax
))
7621 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
7622 SCM_I_BIG_MPZ (val
));
7624 return n
>= min
&& n
<= max
;
7632 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
7634 scm_error (scm_out_of_range_key
,
7636 "Value out of range ~S to ~S: ~S",
7637 scm_list_3 (min
, max
, bad_val
),
7638 scm_list_1 (bad_val
));
7641 #define TYPE scm_t_intmax
7642 #define TYPE_MIN min
7643 #define TYPE_MAX max
7644 #define SIZEOF_TYPE 0
7645 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
7646 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
7647 #include "libguile/conv-integer.i.c"
7649 #define TYPE scm_t_uintmax
7650 #define TYPE_MIN min
7651 #define TYPE_MAX max
7652 #define SIZEOF_TYPE 0
7653 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
7654 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
7655 #include "libguile/conv-uinteger.i.c"
7657 #define TYPE scm_t_int8
7658 #define TYPE_MIN SCM_T_INT8_MIN
7659 #define TYPE_MAX SCM_T_INT8_MAX
7660 #define SIZEOF_TYPE 1
7661 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
7662 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
7663 #include "libguile/conv-integer.i.c"
7665 #define TYPE scm_t_uint8
7667 #define TYPE_MAX SCM_T_UINT8_MAX
7668 #define SIZEOF_TYPE 1
7669 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
7670 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
7671 #include "libguile/conv-uinteger.i.c"
7673 #define TYPE scm_t_int16
7674 #define TYPE_MIN SCM_T_INT16_MIN
7675 #define TYPE_MAX SCM_T_INT16_MAX
7676 #define SIZEOF_TYPE 2
7677 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
7678 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
7679 #include "libguile/conv-integer.i.c"
7681 #define TYPE scm_t_uint16
7683 #define TYPE_MAX SCM_T_UINT16_MAX
7684 #define SIZEOF_TYPE 2
7685 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
7686 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
7687 #include "libguile/conv-uinteger.i.c"
7689 #define TYPE scm_t_int32
7690 #define TYPE_MIN SCM_T_INT32_MIN
7691 #define TYPE_MAX SCM_T_INT32_MAX
7692 #define SIZEOF_TYPE 4
7693 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
7694 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
7695 #include "libguile/conv-integer.i.c"
7697 #define TYPE scm_t_uint32
7699 #define TYPE_MAX SCM_T_UINT32_MAX
7700 #define SIZEOF_TYPE 4
7701 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
7702 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
7703 #include "libguile/conv-uinteger.i.c"
7705 #define TYPE scm_t_wchar
7706 #define TYPE_MIN (scm_t_int32)-1
7707 #define TYPE_MAX (scm_t_int32)0x10ffff
7708 #define SIZEOF_TYPE 4
7709 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
7710 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
7711 #include "libguile/conv-integer.i.c"
7713 #define TYPE scm_t_int64
7714 #define TYPE_MIN SCM_T_INT64_MIN
7715 #define TYPE_MAX SCM_T_INT64_MAX
7716 #define SIZEOF_TYPE 8
7717 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
7718 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
7719 #include "libguile/conv-integer.i.c"
7721 #define TYPE scm_t_uint64
7723 #define TYPE_MAX SCM_T_UINT64_MAX
7724 #define SIZEOF_TYPE 8
7725 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
7726 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
7727 #include "libguile/conv-uinteger.i.c"
7730 scm_to_mpz (SCM val
, mpz_t rop
)
7732 if (SCM_I_INUMP (val
))
7733 mpz_set_si (rop
, SCM_I_INUM (val
));
7734 else if (SCM_BIGP (val
))
7735 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
7737 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
7741 scm_from_mpz (mpz_t val
)
7743 return scm_i_mpz2num (val
);
7747 scm_is_real (SCM val
)
7749 return scm_is_true (scm_real_p (val
));
7753 scm_is_rational (SCM val
)
7755 return scm_is_true (scm_rational_p (val
));
7759 scm_to_double (SCM val
)
7761 if (SCM_I_INUMP (val
))
7762 return SCM_I_INUM (val
);
7763 else if (SCM_BIGP (val
))
7764 return scm_i_big2dbl (val
);
7765 else if (SCM_FRACTIONP (val
))
7766 return scm_i_fraction2double (val
);
7767 else if (SCM_REALP (val
))
7768 return SCM_REAL_VALUE (val
);
7770 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
7774 scm_from_double (double val
)
7778 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
7780 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
7781 SCM_REAL_VALUE (z
) = val
;
7786 #if SCM_ENABLE_DEPRECATED == 1
7789 scm_num2float (SCM num
, unsigned long pos
, const char *s_caller
)
7791 scm_c_issue_deprecation_warning
7792 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
7796 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
7800 scm_out_of_range (NULL
, num
);
7803 return scm_to_double (num
);
7807 scm_num2double (SCM num
, unsigned long pos
, const char *s_caller
)
7809 scm_c_issue_deprecation_warning
7810 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
7814 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
7818 scm_out_of_range (NULL
, num
);
7821 return scm_to_double (num
);
7827 scm_is_complex (SCM val
)
7829 return scm_is_true (scm_complex_p (val
));
7833 scm_c_real_part (SCM z
)
7835 if (SCM_COMPLEXP (z
))
7836 return SCM_COMPLEX_REAL (z
);
7839 /* Use the scm_real_part to get proper error checking and
7842 return scm_to_double (scm_real_part (z
));
7847 scm_c_imag_part (SCM z
)
7849 if (SCM_COMPLEXP (z
))
7850 return SCM_COMPLEX_IMAG (z
);
7853 /* Use the scm_imag_part to get proper error checking and
7854 dispatching. The result will almost always be 0.0, but not
7857 return scm_to_double (scm_imag_part (z
));
7862 scm_c_magnitude (SCM z
)
7864 return scm_to_double (scm_magnitude (z
));
7870 return scm_to_double (scm_angle (z
));
7874 scm_is_number (SCM z
)
7876 return scm_is_true (scm_number_p (z
));
7880 /* In the following functions we dispatch to the real-arg funcs like log()
7881 when we know the arg is real, instead of just handing everything to
7882 clog() for instance. This is in case clog() doesn't optimize for a
7883 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
7884 well use it to go straight to the applicable C func. */
7886 SCM_PRIMITIVE_GENERIC (scm_log
, "log", 1, 0, 0,
7888 "Return the natural logarithm of @var{z}.")
7889 #define FUNC_NAME s_scm_log
7891 if (SCM_COMPLEXP (z
))
7893 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
7894 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
7896 double re
= SCM_COMPLEX_REAL (z
);
7897 double im
= SCM_COMPLEX_IMAG (z
);
7898 return scm_c_make_rectangular (log (hypot (re
, im
)),
7902 else if (SCM_NUMBERP (z
))
7904 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
7905 although the value itself overflows. */
7906 double re
= scm_to_double (z
);
7907 double l
= log (fabs (re
));
7909 return scm_from_double (l
);
7911 return scm_c_make_rectangular (l
, M_PI
);
7914 SCM_WTA_DISPATCH_1 (g_scm_log
, z
, 1, s_scm_log
);
7919 SCM_PRIMITIVE_GENERIC (scm_log10
, "log10", 1, 0, 0,
7921 "Return the base 10 logarithm of @var{z}.")
7922 #define FUNC_NAME s_scm_log10
7924 if (SCM_COMPLEXP (z
))
7926 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
7927 clog() and a multiply by M_LOG10E, rather than the fallback
7928 log10+hypot+atan2.) */
7929 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
7930 && defined SCM_COMPLEX_VALUE
7931 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
7933 double re
= SCM_COMPLEX_REAL (z
);
7934 double im
= SCM_COMPLEX_IMAG (z
);
7935 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
7936 M_LOG10E
* atan2 (im
, re
));
7939 else if (SCM_NUMBERP (z
))
7941 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
7942 although the value itself overflows. */
7943 double re
= scm_to_double (z
);
7944 double l
= log10 (fabs (re
));
7946 return scm_from_double (l
);
7948 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
7951 SCM_WTA_DISPATCH_1 (g_scm_log10
, z
, 1, s_scm_log10
);
7956 SCM_PRIMITIVE_GENERIC (scm_exp
, "exp", 1, 0, 0,
7958 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
7959 "base of natural logarithms (2.71828@dots{}).")
7960 #define FUNC_NAME s_scm_exp
7962 if (SCM_COMPLEXP (z
))
7964 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
7965 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
7967 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
7968 SCM_COMPLEX_IMAG (z
));
7971 else if (SCM_NUMBERP (z
))
7973 /* When z is a negative bignum the conversion to double overflows,
7974 giving -infinity, but that's ok, the exp is still 0.0. */
7975 return scm_from_double (exp (scm_to_double (z
)));
7978 SCM_WTA_DISPATCH_1 (g_scm_exp
, z
, 1, s_scm_exp
);
7983 SCM_PRIMITIVE_GENERIC (scm_sqrt
, "sqrt", 1, 0, 0,
7985 "Return the square root of @var{z}. Of the two possible roots\n"
7986 "(positive and negative), the one with positive real part\n"
7987 "is returned, or if that's zero then a positive imaginary part.\n"
7991 "(sqrt 9.0) @result{} 3.0\n"
7992 "(sqrt -9.0) @result{} 0.0+3.0i\n"
7993 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
7994 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
7996 #define FUNC_NAME s_scm_sqrt
7998 if (SCM_COMPLEXP (z
))
8000 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
8001 && defined SCM_COMPLEX_VALUE
8002 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (z
)));
8004 double re
= SCM_COMPLEX_REAL (z
);
8005 double im
= SCM_COMPLEX_IMAG (z
);
8006 return scm_c_make_polar (sqrt (hypot (re
, im
)),
8007 0.5 * atan2 (im
, re
));
8010 else if (SCM_NUMBERP (z
))
8012 double xx
= scm_to_double (z
);
8014 return scm_c_make_rectangular (0.0, sqrt (-xx
));
8016 return scm_from_double (sqrt (xx
));
8019 SCM_WTA_DISPATCH_1 (g_scm_sqrt
, z
, 1, s_scm_sqrt
);
8030 mpz_init_set_si (z_negative_one
, -1);
8032 /* It may be possible to tune the performance of some algorithms by using
8033 * the following constants to avoid the creation of bignums. Please, before
8034 * using these values, remember the two rules of program optimization:
8035 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
8036 scm_c_define ("most-positive-fixnum",
8037 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
8038 scm_c_define ("most-negative-fixnum",
8039 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
8041 scm_add_feature ("complex");
8042 scm_add_feature ("inexact");
8043 flo0
= scm_from_double (0.0);
8045 /* determine floating point precision */
8046 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
8048 init_dblprec(&scm_dblprec
[i
-2],i
);
8049 init_fx_radix(fx_per_radix
[i
-2],i
);
8052 /* hard code precision for base 10 if the preprocessor tells us to... */
8053 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
8056 exactly_one_half
= scm_divide (SCM_INUM1
, SCM_I_MAKINUM (2));
8057 #include "libguile/numbers.x"