1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public License
9 * as published by the Free Software Foundation; either version 3 of
10 * the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
24 /* General assumptions:
25 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
26 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
27 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
28 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
29 * All objects satisfying SCM_FRACTIONP are never an integer.
34 - see if special casing bignums and reals in integer-exponent when
35 possible (to use mpz_pow and mpf_pow_ui) is faster.
37 - look in to better short-circuiting of common cases in
38 integer-expt and elsewhere.
40 - see if direct mpz operations can help in ash and elsewhere.
57 #include "libguile/_scm.h"
58 #include "libguile/feature.h"
59 #include "libguile/ports.h"
60 #include "libguile/root.h"
61 #include "libguile/smob.h"
62 #include "libguile/strings.h"
63 #include "libguile/bdw-gc.h"
65 #include "libguile/validate.h"
66 #include "libguile/numbers.h"
67 #include "libguile/deprecation.h"
69 #include "libguile/eq.h"
71 /* values per glibc, if not already defined */
73 #define M_LOG10E 0.43429448190325182765
76 #define M_PI 3.14159265358979323846
79 typedef scm_t_signed_bits scm_t_inum
;
80 #define scm_from_inum(x) (scm_from_signed_integer (x))
82 /* Tests to see if a C double is neither infinite nor a NaN.
83 TODO: if it's available, use C99's isfinite(x) instead */
84 #define DOUBLE_IS_FINITE(x) (!isinf(x) && !isnan(x))
89 Wonder if this might be faster for some of our code? A switch on
90 the numtag would jump directly to the right case, and the
91 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
93 #define SCM_I_NUMTAG_NOTNUM 0
94 #define SCM_I_NUMTAG_INUM 1
95 #define SCM_I_NUMTAG_BIG scm_tc16_big
96 #define SCM_I_NUMTAG_REAL scm_tc16_real
97 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
98 #define SCM_I_NUMTAG(x) \
99 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
100 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
101 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
102 : SCM_I_NUMTAG_NOTNUM)))
104 /* the macro above will not work as is with fractions */
108 static SCM exactly_one_half
;
110 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
112 /* FLOBUFLEN is the maximum number of characters neccessary for the
113 * printed or scm_string representation of an inexact number.
115 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
118 #if !defined (HAVE_ASINH)
119 static double asinh (double x
) { return log (x
+ sqrt (x
* x
+ 1)); }
121 #if !defined (HAVE_ACOSH)
122 static double acosh (double x
) { return log (x
+ sqrt (x
* x
- 1)); }
124 #if !defined (HAVE_ATANH)
125 static double atanh (double x
) { return 0.5 * log ((1 + x
) / (1 - x
)); }
128 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
129 an explicit check. In some future gmp (don't know what version number),
130 mpz_cmp_d is supposed to do this itself. */
132 #define xmpz_cmp_d(z, d) \
133 (isinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
135 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
139 #if defined (GUILE_I)
140 #if HAVE_COMPLEX_DOUBLE
142 /* For an SCM object Z which is a complex number (ie. satisfies
143 SCM_COMPLEXP), return its value as a C level "complex double". */
144 #define SCM_COMPLEX_VALUE(z) \
145 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
147 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
149 /* Convert a C "complex double" to an SCM value. */
151 scm_from_complex_double (complex double z
)
153 return scm_c_make_rectangular (creal (z
), cimag (z
));
156 #endif /* HAVE_COMPLEX_DOUBLE */
161 static mpz_t z_negative_one
;
164 /* Clear the `mpz_t' embedded in bignum PTR. */
166 finalize_bignum (GC_PTR ptr
, GC_PTR data
)
170 bignum
= PTR2SCM (ptr
);
171 mpz_clear (SCM_I_BIG_MPZ (bignum
));
174 /* Return a new uninitialized bignum. */
179 GC_finalization_proc prev_finalizer
;
180 GC_PTR prev_finalizer_data
;
182 /* Allocate one word for the type tag and enough room for an `mpz_t'. */
183 p
= scm_gc_malloc_pointerless (sizeof (scm_t_bits
) + sizeof (mpz_t
),
187 GC_REGISTER_FINALIZER_NO_ORDER (p
, finalize_bignum
, NULL
,
189 &prev_finalizer_data
);
198 /* Return a newly created bignum. */
199 SCM z
= make_bignum ();
200 mpz_init (SCM_I_BIG_MPZ (z
));
205 scm_i_inum2big (scm_t_inum x
)
207 /* Return a newly created bignum initialized to X. */
208 SCM z
= make_bignum ();
209 #if SIZEOF_VOID_P == SIZEOF_LONG
210 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
212 /* Note that in this case, you'll also have to check all mpz_*_ui and
213 mpz_*_si invocations in Guile. */
214 #error creation of mpz not implemented for this inum size
220 scm_i_long2big (long x
)
222 /* Return a newly created bignum initialized to X. */
223 SCM z
= make_bignum ();
224 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
229 scm_i_ulong2big (unsigned long x
)
231 /* Return a newly created bignum initialized to X. */
232 SCM z
= make_bignum ();
233 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
238 scm_i_clonebig (SCM src_big
, int same_sign_p
)
240 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
241 SCM z
= make_bignum ();
242 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
244 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
249 scm_i_bigcmp (SCM x
, SCM y
)
251 /* Return neg if x < y, pos if x > y, and 0 if x == y */
252 /* presume we already know x and y are bignums */
253 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
254 scm_remember_upto_here_2 (x
, y
);
259 scm_i_dbl2big (double d
)
261 /* results are only defined if d is an integer */
262 SCM z
= make_bignum ();
263 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
267 /* Convert a integer in double representation to a SCM number. */
270 scm_i_dbl2num (double u
)
272 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
273 powers of 2, so there's no rounding when making "double" values
274 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
275 get rounded on a 64-bit machine, hence the "+1".
277 The use of floor() to force to an integer value ensures we get a
278 "numerically closest" value without depending on how a
279 double->long cast or how mpz_set_d will round. For reference,
280 double->long probably follows the hardware rounding mode,
281 mpz_set_d truncates towards zero. */
283 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
284 representable as a double? */
286 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
287 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
288 return SCM_I_MAKINUM ((scm_t_inum
) u
);
290 return scm_i_dbl2big (u
);
293 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
294 with R5RS exact->inexact.
296 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
297 (ie. truncate towards zero), then adjust to get the closest double by
298 examining the next lower bit and adding 1 (to the absolute value) if
301 Bignums exactly half way between representable doubles are rounded to the
302 next higher absolute value (ie. away from zero). This seems like an
303 adequate interpretation of R5RS "numerically closest", and it's easier
304 and faster than a full "nearest-even" style.
306 The bit test must be done on the absolute value of the mpz_t, which means
307 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
308 negatives as twos complement.
310 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
311 following the hardware rounding mode, but applied to the absolute value
312 of the mpz_t operand. This is not what we want so we put the high
313 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
314 mpz_get_d is supposed to always truncate towards zero.
316 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
317 is a slowdown. It'd be faster to pick out the relevant high bits with
318 mpz_getlimbn if we could be bothered coding that, and if the new
319 truncating gmp doesn't come out. */
322 scm_i_big2dbl (SCM b
)
327 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
331 /* Current GMP, eg. 4.1.3, force truncation towards zero */
333 if (bits
> DBL_MANT_DIG
)
335 size_t shift
= bits
- DBL_MANT_DIG
;
336 mpz_init2 (tmp
, DBL_MANT_DIG
);
337 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
338 result
= ldexp (mpz_get_d (tmp
), shift
);
343 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
348 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
351 if (bits
> DBL_MANT_DIG
)
353 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
354 /* test bit number "pos" in absolute value */
355 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
356 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
358 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
362 scm_remember_upto_here_1 (b
);
367 scm_i_normbig (SCM b
)
369 /* convert a big back to a fixnum if it'll fit */
370 /* presume b is a bignum */
371 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
373 scm_t_inum val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
374 if (SCM_FIXABLE (val
))
375 b
= SCM_I_MAKINUM (val
);
380 static SCM_C_INLINE_KEYWORD SCM
381 scm_i_mpz2num (mpz_t b
)
383 /* convert a mpz number to a SCM number. */
384 if (mpz_fits_slong_p (b
))
386 scm_t_inum val
= mpz_get_si (b
);
387 if (SCM_FIXABLE (val
))
388 return SCM_I_MAKINUM (val
);
392 SCM z
= make_bignum ();
393 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
398 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
399 static SCM
scm_divide2real (SCM x
, SCM y
);
402 scm_i_make_ratio (SCM numerator
, SCM denominator
)
403 #define FUNC_NAME "make-ratio"
405 /* First make sure the arguments are proper.
407 if (SCM_I_INUMP (denominator
))
409 if (scm_is_eq (denominator
, SCM_INUM0
))
410 scm_num_overflow ("make-ratio");
411 if (scm_is_eq (denominator
, SCM_INUM1
))
416 if (!(SCM_BIGP(denominator
)))
417 SCM_WRONG_TYPE_ARG (2, denominator
);
419 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
420 SCM_WRONG_TYPE_ARG (1, numerator
);
422 /* Then flip signs so that the denominator is positive.
424 if (scm_is_true (scm_negative_p (denominator
)))
426 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
427 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
430 /* Now consider for each of the four fixnum/bignum combinations
431 whether the rational number is really an integer.
433 if (SCM_I_INUMP (numerator
))
435 scm_t_inum x
= SCM_I_INUM (numerator
);
436 if (scm_is_eq (numerator
, SCM_INUM0
))
438 if (SCM_I_INUMP (denominator
))
441 y
= SCM_I_INUM (denominator
);
445 return SCM_I_MAKINUM (x
/ y
);
449 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
450 of that value for the denominator, as a bignum. Apart from
451 that case, abs(bignum) > abs(inum) so inum/bignum is not an
453 if (x
== SCM_MOST_NEGATIVE_FIXNUM
454 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
455 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
456 return SCM_I_MAKINUM(-1);
459 else if (SCM_BIGP (numerator
))
461 if (SCM_I_INUMP (denominator
))
463 scm_t_inum yy
= SCM_I_INUM (denominator
);
464 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
465 return scm_divide (numerator
, denominator
);
469 if (scm_is_eq (numerator
, denominator
))
471 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
472 SCM_I_BIG_MPZ (denominator
)))
473 return scm_divide(numerator
, denominator
);
477 /* No, it's a proper fraction.
480 SCM divisor
= scm_gcd (numerator
, denominator
);
481 if (!(scm_is_eq (divisor
, SCM_INUM1
)))
483 numerator
= scm_divide (numerator
, divisor
);
484 denominator
= scm_divide (denominator
, divisor
);
487 return scm_double_cell (scm_tc16_fraction
,
488 SCM_UNPACK (numerator
),
489 SCM_UNPACK (denominator
), 0);
495 scm_i_fraction2double (SCM z
)
497 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
498 SCM_FRACTION_DENOMINATOR (z
)));
502 double_is_non_negative_zero (double x
)
504 static double zero
= 0.0;
506 return !memcmp (&x
, &zero
, sizeof(double));
509 SCM_PRIMITIVE_GENERIC (scm_exact_p
, "exact?", 1, 0, 0,
511 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
513 #define FUNC_NAME s_scm_exact_p
515 if (SCM_INEXACTP (x
))
517 else if (SCM_NUMBERP (x
))
520 SCM_WTA_DISPATCH_1 (g_scm_exact_p
, x
, 1, s_scm_exact_p
);
525 SCM_PRIMITIVE_GENERIC (scm_inexact_p
, "inexact?", 1, 0, 0,
527 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
529 #define FUNC_NAME s_scm_inexact_p
531 if (SCM_INEXACTP (x
))
533 else if (SCM_NUMBERP (x
))
536 SCM_WTA_DISPATCH_1 (g_scm_inexact_p
, x
, 1, s_scm_inexact_p
);
541 SCM_PRIMITIVE_GENERIC (scm_odd_p
, "odd?", 1, 0, 0,
543 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
545 #define FUNC_NAME s_scm_odd_p
549 scm_t_inum val
= SCM_I_INUM (n
);
550 return scm_from_bool ((val
& 1L) != 0);
552 else if (SCM_BIGP (n
))
554 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
555 scm_remember_upto_here_1 (n
);
556 return scm_from_bool (odd_p
);
558 else if (SCM_REALP (n
))
560 double val
= SCM_REAL_VALUE (n
);
561 if (DOUBLE_IS_FINITE (val
))
563 double rem
= fabs (fmod (val
, 2.0));
570 SCM_WTA_DISPATCH_1 (g_scm_odd_p
, n
, 1, s_scm_odd_p
);
575 SCM_PRIMITIVE_GENERIC (scm_even_p
, "even?", 1, 0, 0,
577 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
579 #define FUNC_NAME s_scm_even_p
583 scm_t_inum val
= SCM_I_INUM (n
);
584 return scm_from_bool ((val
& 1L) == 0);
586 else if (SCM_BIGP (n
))
588 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
589 scm_remember_upto_here_1 (n
);
590 return scm_from_bool (even_p
);
592 else if (SCM_REALP (n
))
594 double val
= SCM_REAL_VALUE (n
);
595 if (DOUBLE_IS_FINITE (val
))
597 double rem
= fabs (fmod (val
, 2.0));
604 SCM_WTA_DISPATCH_1 (g_scm_even_p
, n
, 1, s_scm_even_p
);
608 SCM_PRIMITIVE_GENERIC (scm_finite_p
, "finite?", 1, 0, 0,
610 "Return @code{#t} if the real number @var{x} is neither\n"
611 "infinite nor a NaN, @code{#f} otherwise.")
612 #define FUNC_NAME s_scm_finite_p
615 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
616 else if (scm_is_real (x
))
619 SCM_WTA_DISPATCH_1 (g_scm_finite_p
, x
, 1, s_scm_finite_p
);
623 SCM_PRIMITIVE_GENERIC (scm_inf_p
, "inf?", 1, 0, 0,
625 "Return @code{#t} if the real number @var{x} is @samp{+inf.0} or\n"
626 "@samp{-inf.0}. Otherwise return @code{#f}.")
627 #define FUNC_NAME s_scm_inf_p
630 return scm_from_bool (isinf (SCM_REAL_VALUE (x
)));
631 else if (scm_is_real (x
))
634 SCM_WTA_DISPATCH_1 (g_scm_inf_p
, x
, 1, s_scm_inf_p
);
638 SCM_PRIMITIVE_GENERIC (scm_nan_p
, "nan?", 1, 0, 0,
640 "Return @code{#t} if the real number @var{x} is a NaN,\n"
641 "or @code{#f} otherwise.")
642 #define FUNC_NAME s_scm_nan_p
645 return scm_from_bool (isnan (SCM_REAL_VALUE (x
)));
646 else if (scm_is_real (x
))
649 SCM_WTA_DISPATCH_1 (g_scm_nan_p
, x
, 1, s_scm_nan_p
);
653 /* Guile's idea of infinity. */
654 static double guile_Inf
;
656 /* Guile's idea of not a number. */
657 static double guile_NaN
;
660 guile_ieee_init (void)
662 /* Some version of gcc on some old version of Linux used to crash when
663 trying to make Inf and NaN. */
666 /* C99 INFINITY, when available.
667 FIXME: The standard allows for INFINITY to be something that overflows
668 at compile time. We ought to have a configure test to check for that
669 before trying to use it. (But in practice we believe this is not a
670 problem on any system guile is likely to target.) */
671 guile_Inf
= INFINITY
;
672 #elif defined HAVE_DINFINITY
674 extern unsigned int DINFINITY
[2];
675 guile_Inf
= (*((double *) (DINFINITY
)));
682 if (guile_Inf
== tmp
)
689 /* C99 NAN, when available */
691 #elif defined HAVE_DQNAN
694 extern unsigned int DQNAN
[2];
695 guile_NaN
= (*((double *)(DQNAN
)));
698 guile_NaN
= guile_Inf
/ guile_Inf
;
702 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
705 #define FUNC_NAME s_scm_inf
707 static int initialized
= 0;
713 return scm_from_double (guile_Inf
);
717 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
720 #define FUNC_NAME s_scm_nan
722 static int initialized
= 0;
728 return scm_from_double (guile_NaN
);
733 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
735 "Return the absolute value of @var{x}.")
736 #define FUNC_NAME s_scm_abs
740 scm_t_inum xx
= SCM_I_INUM (x
);
743 else if (SCM_POSFIXABLE (-xx
))
744 return SCM_I_MAKINUM (-xx
);
746 return scm_i_inum2big (-xx
);
748 else if (SCM_LIKELY (SCM_REALP (x
)))
750 double xx
= SCM_REAL_VALUE (x
);
751 /* If x is a NaN then xx<0 is false so we return x unchanged */
753 return scm_from_double (-xx
);
754 /* Handle signed zeroes properly */
755 else if (SCM_UNLIKELY (xx
== 0.0))
760 else if (SCM_BIGP (x
))
762 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
764 return scm_i_clonebig (x
, 0);
768 else if (SCM_FRACTIONP (x
))
770 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
772 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
773 SCM_FRACTION_DENOMINATOR (x
));
776 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
781 SCM_PRIMITIVE_GENERIC (scm_quotient
, "quotient", 2, 0, 0,
783 "Return the quotient of the numbers @var{x} and @var{y}.")
784 #define FUNC_NAME s_scm_quotient
786 if (SCM_LIKELY (SCM_I_INUMP (x
)))
788 scm_t_inum xx
= SCM_I_INUM (x
);
789 if (SCM_LIKELY (SCM_I_INUMP (y
)))
791 scm_t_inum yy
= SCM_I_INUM (y
);
792 if (SCM_UNLIKELY (yy
== 0))
793 scm_num_overflow (s_scm_quotient
);
796 scm_t_inum z
= xx
/ yy
;
797 if (SCM_LIKELY (SCM_FIXABLE (z
)))
798 return SCM_I_MAKINUM (z
);
800 return scm_i_inum2big (z
);
803 else if (SCM_BIGP (y
))
805 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
806 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
807 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
809 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
810 scm_remember_upto_here_1 (y
);
811 return SCM_I_MAKINUM (-1);
817 SCM_WTA_DISPATCH_2 (g_scm_quotient
, x
, y
, SCM_ARG2
, s_scm_quotient
);
819 else if (SCM_BIGP (x
))
821 if (SCM_LIKELY (SCM_I_INUMP (y
)))
823 scm_t_inum yy
= SCM_I_INUM (y
);
824 if (SCM_UNLIKELY (yy
== 0))
825 scm_num_overflow (s_scm_quotient
);
826 else if (SCM_UNLIKELY (yy
== 1))
830 SCM result
= scm_i_mkbig ();
833 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
836 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
839 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
840 scm_remember_upto_here_1 (x
);
841 return scm_i_normbig (result
);
844 else if (SCM_BIGP (y
))
846 SCM result
= scm_i_mkbig ();
847 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
850 scm_remember_upto_here_2 (x
, y
);
851 return scm_i_normbig (result
);
854 SCM_WTA_DISPATCH_2 (g_scm_quotient
, x
, y
, SCM_ARG2
, s_scm_quotient
);
857 SCM_WTA_DISPATCH_2 (g_scm_quotient
, x
, y
, SCM_ARG1
, s_scm_quotient
);
861 SCM_PRIMITIVE_GENERIC (scm_remainder
, "remainder", 2, 0, 0,
863 "Return the remainder of the numbers @var{x} and @var{y}.\n"
865 "(remainder 13 4) @result{} 1\n"
866 "(remainder -13 4) @result{} -1\n"
868 #define FUNC_NAME s_scm_remainder
870 if (SCM_LIKELY (SCM_I_INUMP (x
)))
872 if (SCM_LIKELY (SCM_I_INUMP (y
)))
874 scm_t_inum yy
= SCM_I_INUM (y
);
875 if (SCM_UNLIKELY (yy
== 0))
876 scm_num_overflow (s_scm_remainder
);
879 /* C99 specifies that "%" is the remainder corresponding to a
880 quotient rounded towards zero, and that's also traditional
881 for machine division, so z here should be well defined. */
882 scm_t_inum z
= SCM_I_INUM (x
) % yy
;
883 return SCM_I_MAKINUM (z
);
886 else if (SCM_BIGP (y
))
888 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
889 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
890 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
892 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
893 scm_remember_upto_here_1 (y
);
900 SCM_WTA_DISPATCH_2 (g_scm_remainder
, x
, y
, SCM_ARG2
, s_scm_remainder
);
902 else if (SCM_BIGP (x
))
904 if (SCM_LIKELY (SCM_I_INUMP (y
)))
906 scm_t_inum yy
= SCM_I_INUM (y
);
907 if (SCM_UNLIKELY (yy
== 0))
908 scm_num_overflow (s_scm_remainder
);
911 SCM result
= scm_i_mkbig ();
914 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
915 scm_remember_upto_here_1 (x
);
916 return scm_i_normbig (result
);
919 else if (SCM_BIGP (y
))
921 SCM result
= scm_i_mkbig ();
922 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
925 scm_remember_upto_here_2 (x
, y
);
926 return scm_i_normbig (result
);
929 SCM_WTA_DISPATCH_2 (g_scm_remainder
, x
, y
, SCM_ARG2
, s_scm_remainder
);
932 SCM_WTA_DISPATCH_2 (g_scm_remainder
, x
, y
, SCM_ARG1
, s_scm_remainder
);
937 SCM_PRIMITIVE_GENERIC (scm_modulo
, "modulo", 2, 0, 0,
939 "Return the modulo of the numbers @var{x} and @var{y}.\n"
941 "(modulo 13 4) @result{} 1\n"
942 "(modulo -13 4) @result{} 3\n"
944 #define FUNC_NAME s_scm_modulo
946 if (SCM_LIKELY (SCM_I_INUMP (x
)))
948 scm_t_inum xx
= SCM_I_INUM (x
);
949 if (SCM_LIKELY (SCM_I_INUMP (y
)))
951 scm_t_inum yy
= SCM_I_INUM (y
);
952 if (SCM_UNLIKELY (yy
== 0))
953 scm_num_overflow (s_scm_modulo
);
956 /* C99 specifies that "%" is the remainder corresponding to a
957 quotient rounded towards zero, and that's also traditional
958 for machine division, so z here should be well defined. */
959 scm_t_inum z
= xx
% yy
;
976 return SCM_I_MAKINUM (result
);
979 else if (SCM_BIGP (y
))
981 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
988 SCM pos_y
= scm_i_clonebig (y
, 0);
989 /* do this after the last scm_op */
990 mpz_init_set_si (z_x
, xx
);
991 result
= pos_y
; /* re-use this bignum */
992 mpz_mod (SCM_I_BIG_MPZ (result
),
994 SCM_I_BIG_MPZ (pos_y
));
995 scm_remember_upto_here_1 (pos_y
);
999 result
= scm_i_mkbig ();
1000 /* do this after the last scm_op */
1001 mpz_init_set_si (z_x
, xx
);
1002 mpz_mod (SCM_I_BIG_MPZ (result
),
1005 scm_remember_upto_here_1 (y
);
1008 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1009 mpz_add (SCM_I_BIG_MPZ (result
),
1011 SCM_I_BIG_MPZ (result
));
1012 scm_remember_upto_here_1 (y
);
1013 /* and do this before the next one */
1015 return scm_i_normbig (result
);
1019 SCM_WTA_DISPATCH_2 (g_scm_modulo
, x
, y
, SCM_ARG2
, s_scm_modulo
);
1021 else if (SCM_BIGP (x
))
1023 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1025 scm_t_inum yy
= SCM_I_INUM (y
);
1026 if (SCM_UNLIKELY (yy
== 0))
1027 scm_num_overflow (s_scm_modulo
);
1030 SCM result
= scm_i_mkbig ();
1031 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
1033 (yy
< 0) ? - yy
: yy
);
1034 scm_remember_upto_here_1 (x
);
1035 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1036 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
1037 SCM_I_BIG_MPZ (result
),
1039 return scm_i_normbig (result
);
1042 else if (SCM_BIGP (y
))
1044 SCM result
= scm_i_mkbig ();
1045 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
1046 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1047 mpz_mod (SCM_I_BIG_MPZ (result
),
1049 SCM_I_BIG_MPZ (pos_y
));
1051 scm_remember_upto_here_1 (x
);
1052 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1053 mpz_add (SCM_I_BIG_MPZ (result
),
1055 SCM_I_BIG_MPZ (result
));
1056 scm_remember_upto_here_2 (y
, pos_y
);
1057 return scm_i_normbig (result
);
1060 SCM_WTA_DISPATCH_2 (g_scm_modulo
, x
, y
, SCM_ARG2
, s_scm_modulo
);
1063 SCM_WTA_DISPATCH_2 (g_scm_modulo
, x
, y
, SCM_ARG1
, s_scm_modulo
);
1067 static SCM
scm_i_inexact_euclidean_quotient (double x
, double y
);
1068 static SCM
scm_i_slow_exact_euclidean_quotient (SCM x
, SCM y
);
1070 SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient
, "euclidean-quotient", 2, 0, 0,
1072 "Return the integer @var{q} such that\n"
1073 "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1074 "where @math{0 <= @var{r} < abs(@var{y})}.\n"
1076 "(euclidean-quotient 123 10) @result{} 12\n"
1077 "(euclidean-quotient 123 -10) @result{} -12\n"
1078 "(euclidean-quotient -123 10) @result{} -13\n"
1079 "(euclidean-quotient -123 -10) @result{} 13\n"
1080 "(euclidean-quotient -123.2 -63.5) @result{} 2.0\n"
1081 "(euclidean-quotient 16/3 -10/7) @result{} -3\n"
1083 #define FUNC_NAME s_scm_euclidean_quotient
1085 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1087 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1089 scm_t_inum yy
= SCM_I_INUM (y
);
1090 if (SCM_UNLIKELY (yy
== 0))
1091 scm_num_overflow (s_scm_euclidean_quotient
);
1094 scm_t_inum xx
= SCM_I_INUM (x
);
1095 scm_t_inum qq
= xx
/ yy
;
1103 return SCM_I_MAKINUM (qq
);
1106 else if (SCM_BIGP (y
))
1108 if (SCM_I_INUM (x
) >= 0)
1111 return SCM_I_MAKINUM (- mpz_sgn (SCM_I_BIG_MPZ (y
)));
1113 else if (SCM_REALP (y
))
1114 return scm_i_inexact_euclidean_quotient
1115 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1116 else if (SCM_FRACTIONP (y
))
1117 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1119 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1120 s_scm_euclidean_quotient
);
1122 else if (SCM_BIGP (x
))
1124 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1126 scm_t_inum yy
= SCM_I_INUM (y
);
1127 if (SCM_UNLIKELY (yy
== 0))
1128 scm_num_overflow (s_scm_euclidean_quotient
);
1131 SCM q
= scm_i_mkbig ();
1133 mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (x
), yy
);
1136 mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (x
), -yy
);
1137 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1139 scm_remember_upto_here_1 (x
);
1140 return scm_i_normbig (q
);
1143 else if (SCM_BIGP (y
))
1145 SCM q
= scm_i_mkbig ();
1146 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1147 mpz_fdiv_q (SCM_I_BIG_MPZ (q
),
1151 mpz_cdiv_q (SCM_I_BIG_MPZ (q
),
1154 scm_remember_upto_here_2 (x
, y
);
1155 return scm_i_normbig (q
);
1157 else if (SCM_REALP (y
))
1158 return scm_i_inexact_euclidean_quotient
1159 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1160 else if (SCM_FRACTIONP (y
))
1161 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1163 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1164 s_scm_euclidean_quotient
);
1166 else if (SCM_REALP (x
))
1168 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1169 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1170 return scm_i_inexact_euclidean_quotient
1171 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1173 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1174 s_scm_euclidean_quotient
);
1176 else if (SCM_FRACTIONP (x
))
1179 return scm_i_inexact_euclidean_quotient
1180 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1182 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1185 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG1
,
1186 s_scm_euclidean_quotient
);
1191 scm_i_inexact_euclidean_quotient (double x
, double y
)
1193 if (SCM_LIKELY (y
> 0))
1194 return scm_from_double (floor (x
/ y
));
1195 else if (SCM_LIKELY (y
< 0))
1196 return scm_from_double (ceil (x
/ y
));
1198 scm_num_overflow (s_scm_euclidean_quotient
); /* or return a NaN? */
1203 /* Compute exact euclidean_quotient the slow way.
1204 We use this only if both arguments are exact,
1205 and at least one of them is a fraction */
1207 scm_i_slow_exact_euclidean_quotient (SCM x
, SCM y
)
1209 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1210 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG1
,
1211 s_scm_euclidean_quotient
);
1212 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1213 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1214 s_scm_euclidean_quotient
);
1215 else if (scm_is_true (scm_positive_p (y
)))
1216 return scm_floor (scm_divide (x
, y
));
1217 else if (scm_is_true (scm_negative_p (y
)))
1218 return scm_ceiling (scm_divide (x
, y
));
1220 scm_num_overflow (s_scm_euclidean_quotient
);
1223 static SCM
scm_i_inexact_euclidean_remainder (double x
, double y
);
1224 static SCM
scm_i_slow_exact_euclidean_remainder (SCM x
, SCM y
);
1226 SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder
, "euclidean-remainder", 2, 0, 0,
1228 "Return the real number @var{r} such that\n"
1229 "@math{0 <= @var{r} < abs(@var{y})} and\n"
1230 "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1231 "for some integer @var{q}.\n"
1233 "(euclidean-remainder 123 10) @result{} 3\n"
1234 "(euclidean-remainder 123 -10) @result{} 3\n"
1235 "(euclidean-remainder -123 10) @result{} 7\n"
1236 "(euclidean-remainder -123 -10) @result{} 7\n"
1237 "(euclidean-remainder -123.2 -63.5) @result{} 3.8\n"
1238 "(euclidean-remainder 16/3 -10/7) @result{} 22/21\n"
1240 #define FUNC_NAME s_scm_euclidean_remainder
1242 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1244 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1246 scm_t_inum yy
= SCM_I_INUM (y
);
1247 if (SCM_UNLIKELY (yy
== 0))
1248 scm_num_overflow (s_scm_euclidean_remainder
);
1251 scm_t_inum rr
= SCM_I_INUM (x
) % yy
;
1253 return SCM_I_MAKINUM (rr
);
1255 return SCM_I_MAKINUM (rr
+ yy
);
1257 return SCM_I_MAKINUM (rr
- yy
);
1260 else if (SCM_BIGP (y
))
1262 scm_t_inum xx
= SCM_I_INUM (x
);
1265 else if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1267 SCM r
= scm_i_mkbig ();
1268 mpz_sub_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1269 scm_remember_upto_here_1 (y
);
1270 return scm_i_normbig (r
);
1274 SCM r
= scm_i_mkbig ();
1275 mpz_add_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1276 scm_remember_upto_here_1 (y
);
1277 mpz_neg (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (r
));
1278 return scm_i_normbig (r
);
1281 else if (SCM_REALP (y
))
1282 return scm_i_inexact_euclidean_remainder
1283 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1284 else if (SCM_FRACTIONP (y
))
1285 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1287 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1288 s_scm_euclidean_remainder
);
1290 else if (SCM_BIGP (x
))
1292 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1294 scm_t_inum yy
= SCM_I_INUM (y
);
1295 if (SCM_UNLIKELY (yy
== 0))
1296 scm_num_overflow (s_scm_euclidean_remainder
);
1302 rr
= mpz_fdiv_ui (SCM_I_BIG_MPZ (x
), yy
);
1303 scm_remember_upto_here_1 (x
);
1304 return SCM_I_MAKINUM (rr
);
1307 else if (SCM_BIGP (y
))
1309 SCM r
= scm_i_mkbig ();
1310 mpz_mod (SCM_I_BIG_MPZ (r
),
1313 scm_remember_upto_here_2 (x
, y
);
1314 return scm_i_normbig (r
);
1316 else if (SCM_REALP (y
))
1317 return scm_i_inexact_euclidean_remainder
1318 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1319 else if (SCM_FRACTIONP (y
))
1320 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1322 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1323 s_scm_euclidean_remainder
);
1325 else if (SCM_REALP (x
))
1327 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1328 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1329 return scm_i_inexact_euclidean_remainder
1330 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1332 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1333 s_scm_euclidean_remainder
);
1335 else if (SCM_FRACTIONP (x
))
1338 return scm_i_inexact_euclidean_remainder
1339 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1341 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1344 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG1
,
1345 s_scm_euclidean_remainder
);
1350 scm_i_inexact_euclidean_remainder (double x
, double y
)
1354 /* Although it would be more efficient to use fmod here, we can't
1355 because it would in some cases produce results inconsistent with
1356 scm_i_inexact_euclidean_quotient, such that x != q * y + r (not
1357 even close). In particular, when x is very close to a multiple of
1358 y, then r might be either 0.0 or abs(y)-epsilon, but those two
1359 cases must correspond to different choices of q. If r = 0.0 then q
1360 must be x/y, and if r = abs(y) then q must be (x-r)/y. If quotient
1361 chooses one and remainder chooses the other, it would be bad. This
1362 problem was observed with x = 130.0 and y = 10/7. */
1363 if (SCM_LIKELY (y
> 0))
1365 else if (SCM_LIKELY (y
< 0))
1368 scm_num_overflow (s_scm_euclidean_remainder
); /* or return a NaN? */
1371 return scm_from_double (x
- q
* y
);
1374 /* Compute exact euclidean_remainder the slow way.
1375 We use this only if both arguments are exact,
1376 and at least one of them is a fraction */
1378 scm_i_slow_exact_euclidean_remainder (SCM x
, SCM y
)
1382 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1383 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG1
,
1384 s_scm_euclidean_remainder
);
1385 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1386 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1387 s_scm_euclidean_remainder
);
1388 else if (scm_is_true (scm_positive_p (y
)))
1389 q
= scm_floor (scm_divide (x
, y
));
1390 else if (scm_is_true (scm_negative_p (y
)))
1391 q
= scm_ceiling (scm_divide (x
, y
));
1393 scm_num_overflow (s_scm_euclidean_remainder
);
1394 return scm_difference (x
, scm_product (y
, q
));
1398 static SCM
scm_i_inexact_euclidean_divide (double x
, double y
);
1399 static SCM
scm_i_slow_exact_euclidean_divide (SCM x
, SCM y
);
1401 SCM_PRIMITIVE_GENERIC (scm_euclidean_divide
, "euclidean/", 2, 0, 0,
1403 "Return the integer @var{q} and the real number @var{r}\n"
1404 "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1405 "and @math{0 <= @var{r} < abs(@var{y})}.\n"
1407 "(euclidean/ 123 10) @result{} 12 and 3\n"
1408 "(euclidean/ 123 -10) @result{} -12 and 3\n"
1409 "(euclidean/ -123 10) @result{} -13 and 7\n"
1410 "(euclidean/ -123 -10) @result{} 13 and 7\n"
1411 "(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
1412 "(euclidean/ 16/3 -10/7) @result{} -3 and 22/21\n"
1414 #define FUNC_NAME s_scm_euclidean_divide
1416 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1418 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1420 scm_t_inum yy
= SCM_I_INUM (y
);
1421 if (SCM_UNLIKELY (yy
== 0))
1422 scm_num_overflow (s_scm_euclidean_divide
);
1425 scm_t_inum xx
= SCM_I_INUM (x
);
1426 scm_t_inum qq
= xx
/ yy
;
1427 scm_t_inum rr
= xx
- qq
* yy
;
1435 return scm_values (scm_list_2 (SCM_I_MAKINUM (qq
),
1436 SCM_I_MAKINUM (rr
)));
1439 else if (SCM_BIGP (y
))
1441 scm_t_inum xx
= SCM_I_INUM (x
);
1443 return scm_values (scm_list_2 (SCM_INUM0
, x
));
1444 else if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1446 SCM r
= scm_i_mkbig ();
1447 mpz_sub_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1448 scm_remember_upto_here_1 (y
);
1450 (scm_list_2 (SCM_I_MAKINUM (-1), scm_i_normbig (r
)));
1454 SCM r
= scm_i_mkbig ();
1455 mpz_add_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1456 scm_remember_upto_here_1 (y
);
1457 mpz_neg (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (r
));
1458 return scm_values (scm_list_2 (SCM_INUM1
, scm_i_normbig (r
)));
1461 else if (SCM_REALP (y
))
1462 return scm_i_inexact_euclidean_divide
1463 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1464 else if (SCM_FRACTIONP (y
))
1465 return scm_i_slow_exact_euclidean_divide (x
, y
);
1467 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1468 s_scm_euclidean_divide
);
1470 else if (SCM_BIGP (x
))
1472 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1474 scm_t_inum yy
= SCM_I_INUM (y
);
1475 if (SCM_UNLIKELY (yy
== 0))
1476 scm_num_overflow (s_scm_euclidean_divide
);
1479 SCM q
= scm_i_mkbig ();
1480 SCM r
= scm_i_mkbig ();
1482 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1483 SCM_I_BIG_MPZ (x
), yy
);
1486 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1487 SCM_I_BIG_MPZ (x
), -yy
);
1488 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1490 scm_remember_upto_here_1 (x
);
1491 return scm_values (scm_list_2 (scm_i_normbig (q
),
1492 scm_i_normbig (r
)));
1495 else if (SCM_BIGP (y
))
1497 SCM q
= scm_i_mkbig ();
1498 SCM r
= scm_i_mkbig ();
1499 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1500 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1501 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1503 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1504 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1505 scm_remember_upto_here_2 (x
, y
);
1506 return scm_values (scm_list_2 (scm_i_normbig (q
),
1507 scm_i_normbig (r
)));
1509 else if (SCM_REALP (y
))
1510 return scm_i_inexact_euclidean_divide
1511 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1512 else if (SCM_FRACTIONP (y
))
1513 return scm_i_slow_exact_euclidean_divide (x
, y
);
1515 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1516 s_scm_euclidean_divide
);
1518 else if (SCM_REALP (x
))
1520 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1521 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1522 return scm_i_inexact_euclidean_divide
1523 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1525 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1526 s_scm_euclidean_divide
);
1528 else if (SCM_FRACTIONP (x
))
1531 return scm_i_inexact_euclidean_divide
1532 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1534 return scm_i_slow_exact_euclidean_divide (x
, y
);
1537 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG1
,
1538 s_scm_euclidean_divide
);
1543 scm_i_inexact_euclidean_divide (double x
, double y
)
1547 if (SCM_LIKELY (y
> 0))
1549 else if (SCM_LIKELY (y
< 0))
1552 scm_num_overflow (s_scm_euclidean_divide
); /* or return a NaN? */
1556 return scm_values (scm_list_2 (scm_from_double (q
),
1557 scm_from_double (r
)));
1560 /* Compute exact euclidean quotient and remainder the slow way.
1561 We use this only if both arguments are exact,
1562 and at least one of them is a fraction */
1564 scm_i_slow_exact_euclidean_divide (SCM x
, SCM y
)
1568 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1569 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG1
,
1570 s_scm_euclidean_divide
);
1571 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1572 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1573 s_scm_euclidean_divide
);
1574 else if (scm_is_true (scm_positive_p (y
)))
1575 q
= scm_floor (scm_divide (x
, y
));
1576 else if (scm_is_true (scm_negative_p (y
)))
1577 q
= scm_ceiling (scm_divide (x
, y
));
1579 scm_num_overflow (s_scm_euclidean_divide
);
1580 r
= scm_difference (x
, scm_product (q
, y
));
1581 return scm_values (scm_list_2 (q
, r
));
1584 static SCM
scm_i_inexact_centered_quotient (double x
, double y
);
1585 static SCM
scm_i_bigint_centered_quotient (SCM x
, SCM y
);
1586 static SCM
scm_i_slow_exact_centered_quotient (SCM x
, SCM y
);
1588 SCM_PRIMITIVE_GENERIC (scm_centered_quotient
, "centered-quotient", 2, 0, 0,
1590 "Return the integer @var{q} such that\n"
1591 "@math{@var{x} = @var{q}*@var{y} + @var{r}} where\n"
1592 "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
1594 "(centered-quotient 123 10) @result{} 12\n"
1595 "(centered-quotient 123 -10) @result{} -12\n"
1596 "(centered-quotient -123 10) @result{} -12\n"
1597 "(centered-quotient -123 -10) @result{} 12\n"
1598 "(centered-quotient -123.2 -63.5) @result{} 2.0\n"
1599 "(centered-quotient 16/3 -10/7) @result{} -4\n"
1601 #define FUNC_NAME s_scm_centered_quotient
1603 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1605 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1607 scm_t_inum yy
= SCM_I_INUM (y
);
1608 if (SCM_UNLIKELY (yy
== 0))
1609 scm_num_overflow (s_scm_centered_quotient
);
1612 scm_t_inum xx
= SCM_I_INUM (x
);
1613 scm_t_inum qq
= xx
/ yy
;
1614 scm_t_inum rr
= xx
- qq
* yy
;
1615 if (SCM_LIKELY (xx
> 0))
1617 if (SCM_LIKELY (yy
> 0))
1619 if (rr
>= (yy
+ 1) / 2)
1624 if (rr
>= (1 - yy
) / 2)
1630 if (SCM_LIKELY (yy
> 0))
1641 return SCM_I_MAKINUM (qq
);
1644 else if (SCM_BIGP (y
))
1646 /* Pass a denormalized bignum version of x (even though it
1647 can fit in a fixnum) to scm_i_bigint_centered_quotient */
1648 return scm_i_bigint_centered_quotient
1649 (scm_i_long2big (SCM_I_INUM (x
)), y
);
1651 else if (SCM_REALP (y
))
1652 return scm_i_inexact_centered_quotient
1653 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1654 else if (SCM_FRACTIONP (y
))
1655 return scm_i_slow_exact_centered_quotient (x
, y
);
1657 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1658 s_scm_centered_quotient
);
1660 else if (SCM_BIGP (x
))
1662 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1664 scm_t_inum yy
= SCM_I_INUM (y
);
1665 if (SCM_UNLIKELY (yy
== 0))
1666 scm_num_overflow (s_scm_centered_quotient
);
1669 SCM q
= scm_i_mkbig ();
1671 /* Arrange for rr to initially be non-positive,
1672 because that simplifies the test to see
1673 if it is within the needed bounds. */
1676 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
1677 SCM_I_BIG_MPZ (x
), yy
);
1678 scm_remember_upto_here_1 (x
);
1680 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
1681 SCM_I_BIG_MPZ (q
), 1);
1685 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
1686 SCM_I_BIG_MPZ (x
), -yy
);
1687 scm_remember_upto_here_1 (x
);
1688 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1690 mpz_add_ui (SCM_I_BIG_MPZ (q
),
1691 SCM_I_BIG_MPZ (q
), 1);
1693 return scm_i_normbig (q
);
1696 else if (SCM_BIGP (y
))
1697 return scm_i_bigint_centered_quotient (x
, y
);
1698 else if (SCM_REALP (y
))
1699 return scm_i_inexact_centered_quotient
1700 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1701 else if (SCM_FRACTIONP (y
))
1702 return scm_i_slow_exact_centered_quotient (x
, y
);
1704 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1705 s_scm_centered_quotient
);
1707 else if (SCM_REALP (x
))
1709 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1710 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1711 return scm_i_inexact_centered_quotient
1712 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1714 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1715 s_scm_centered_quotient
);
1717 else if (SCM_FRACTIONP (x
))
1720 return scm_i_inexact_centered_quotient
1721 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1723 return scm_i_slow_exact_centered_quotient (x
, y
);
1726 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG1
,
1727 s_scm_centered_quotient
);
1732 scm_i_inexact_centered_quotient (double x
, double y
)
1734 if (SCM_LIKELY (y
> 0))
1735 return scm_from_double (floor (x
/y
+ 0.5));
1736 else if (SCM_LIKELY (y
< 0))
1737 return scm_from_double (ceil (x
/y
- 0.5));
1739 scm_num_overflow (s_scm_centered_quotient
); /* or return a NaN? */
1744 /* Assumes that both x and y are bigints, though
1745 x might be able to fit into a fixnum. */
1747 scm_i_bigint_centered_quotient (SCM x
, SCM y
)
1751 /* Note that x might be small enough to fit into a
1752 fixnum, so we must not let it escape into the wild */
1756 /* min_r will eventually become -abs(y)/2 */
1757 min_r
= scm_i_mkbig ();
1758 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
1759 SCM_I_BIG_MPZ (y
), 1);
1761 /* Arrange for rr to initially be non-positive,
1762 because that simplifies the test to see
1763 if it is within the needed bounds. */
1764 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1766 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1767 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1768 scm_remember_upto_here_2 (x
, y
);
1769 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
1770 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1771 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
1772 SCM_I_BIG_MPZ (q
), 1);
1776 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1777 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1778 scm_remember_upto_here_2 (x
, y
);
1779 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1780 mpz_add_ui (SCM_I_BIG_MPZ (q
),
1781 SCM_I_BIG_MPZ (q
), 1);
1783 scm_remember_upto_here_2 (r
, min_r
);
1784 return scm_i_normbig (q
);
1787 /* Compute exact centered quotient the slow way.
1788 We use this only if both arguments are exact,
1789 and at least one of them is a fraction */
1791 scm_i_slow_exact_centered_quotient (SCM x
, SCM y
)
1793 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1794 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG1
,
1795 s_scm_centered_quotient
);
1796 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1797 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1798 s_scm_centered_quotient
);
1799 else if (scm_is_true (scm_positive_p (y
)))
1800 return scm_floor (scm_sum (scm_divide (x
, y
),
1802 else if (scm_is_true (scm_negative_p (y
)))
1803 return scm_ceiling (scm_difference (scm_divide (x
, y
),
1806 scm_num_overflow (s_scm_centered_quotient
);
1809 static SCM
scm_i_inexact_centered_remainder (double x
, double y
);
1810 static SCM
scm_i_bigint_centered_remainder (SCM x
, SCM y
);
1811 static SCM
scm_i_slow_exact_centered_remainder (SCM x
, SCM y
);
1813 SCM_PRIMITIVE_GENERIC (scm_centered_remainder
, "centered-remainder", 2, 0, 0,
1815 "Return the real number @var{r} such that\n"
1816 "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}\n"
1817 "and @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1818 "for some integer @var{q}.\n"
1820 "(centered-remainder 123 10) @result{} 3\n"
1821 "(centered-remainder 123 -10) @result{} 3\n"
1822 "(centered-remainder -123 10) @result{} -3\n"
1823 "(centered-remainder -123 -10) @result{} -3\n"
1824 "(centered-remainder -123.2 -63.5) @result{} 3.8\n"
1825 "(centered-remainder 16/3 -10/7) @result{} -8/21\n"
1827 #define FUNC_NAME s_scm_centered_remainder
1829 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1831 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1833 scm_t_inum yy
= SCM_I_INUM (y
);
1834 if (SCM_UNLIKELY (yy
== 0))
1835 scm_num_overflow (s_scm_centered_remainder
);
1838 scm_t_inum xx
= SCM_I_INUM (x
);
1839 scm_t_inum rr
= xx
% yy
;
1840 if (SCM_LIKELY (xx
> 0))
1842 if (SCM_LIKELY (yy
> 0))
1844 if (rr
>= (yy
+ 1) / 2)
1849 if (rr
>= (1 - yy
) / 2)
1855 if (SCM_LIKELY (yy
> 0))
1866 return SCM_I_MAKINUM (rr
);
1869 else if (SCM_BIGP (y
))
1871 /* Pass a denormalized bignum version of x (even though it
1872 can fit in a fixnum) to scm_i_bigint_centered_remainder */
1873 return scm_i_bigint_centered_remainder
1874 (scm_i_long2big (SCM_I_INUM (x
)), y
);
1876 else if (SCM_REALP (y
))
1877 return scm_i_inexact_centered_remainder
1878 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1879 else if (SCM_FRACTIONP (y
))
1880 return scm_i_slow_exact_centered_remainder (x
, y
);
1882 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1883 s_scm_centered_remainder
);
1885 else if (SCM_BIGP (x
))
1887 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1889 scm_t_inum yy
= SCM_I_INUM (y
);
1890 if (SCM_UNLIKELY (yy
== 0))
1891 scm_num_overflow (s_scm_centered_remainder
);
1895 /* Arrange for rr to initially be non-positive,
1896 because that simplifies the test to see
1897 if it is within the needed bounds. */
1900 rr
= - mpz_cdiv_ui (SCM_I_BIG_MPZ (x
), yy
);
1901 scm_remember_upto_here_1 (x
);
1907 rr
= - mpz_cdiv_ui (SCM_I_BIG_MPZ (x
), -yy
);
1908 scm_remember_upto_here_1 (x
);
1912 return SCM_I_MAKINUM (rr
);
1915 else if (SCM_BIGP (y
))
1916 return scm_i_bigint_centered_remainder (x
, y
);
1917 else if (SCM_REALP (y
))
1918 return scm_i_inexact_centered_remainder
1919 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1920 else if (SCM_FRACTIONP (y
))
1921 return scm_i_slow_exact_centered_remainder (x
, y
);
1923 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1924 s_scm_centered_remainder
);
1926 else if (SCM_REALP (x
))
1928 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1929 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1930 return scm_i_inexact_centered_remainder
1931 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1933 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1934 s_scm_centered_remainder
);
1936 else if (SCM_FRACTIONP (x
))
1939 return scm_i_inexact_centered_remainder
1940 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1942 return scm_i_slow_exact_centered_remainder (x
, y
);
1945 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG1
,
1946 s_scm_centered_remainder
);
1951 scm_i_inexact_centered_remainder (double x
, double y
)
1955 /* Although it would be more efficient to use fmod here, we can't
1956 because it would in some cases produce results inconsistent with
1957 scm_i_inexact_centered_quotient, such that x != r + q * y (not even
1958 close). In particular, when x-y/2 is very close to a multiple of
1959 y, then r might be either -abs(y/2) or abs(y/2)-epsilon, but those
1960 two cases must correspond to different choices of q. If quotient
1961 chooses one and remainder chooses the other, it would be bad. */
1962 if (SCM_LIKELY (y
> 0))
1963 q
= floor (x
/y
+ 0.5);
1964 else if (SCM_LIKELY (y
< 0))
1965 q
= ceil (x
/y
- 0.5);
1967 scm_num_overflow (s_scm_centered_remainder
); /* or return a NaN? */
1970 return scm_from_double (x
- q
* y
);
1973 /* Assumes that both x and y are bigints, though
1974 x might be able to fit into a fixnum. */
1976 scm_i_bigint_centered_remainder (SCM x
, SCM y
)
1980 /* Note that x might be small enough to fit into a
1981 fixnum, so we must not let it escape into the wild */
1984 /* min_r will eventually become -abs(y)/2 */
1985 min_r
= scm_i_mkbig ();
1986 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
1987 SCM_I_BIG_MPZ (y
), 1);
1989 /* Arrange for rr to initially be non-positive,
1990 because that simplifies the test to see
1991 if it is within the needed bounds. */
1992 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1994 mpz_cdiv_r (SCM_I_BIG_MPZ (r
),
1995 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1996 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
1997 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1998 mpz_add (SCM_I_BIG_MPZ (r
),
2004 mpz_fdiv_r (SCM_I_BIG_MPZ (r
),
2005 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2006 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2007 mpz_sub (SCM_I_BIG_MPZ (r
),
2011 scm_remember_upto_here_2 (x
, y
);
2012 return scm_i_normbig (r
);
2015 /* Compute exact centered_remainder the slow way.
2016 We use this only if both arguments are exact,
2017 and at least one of them is a fraction */
2019 scm_i_slow_exact_centered_remainder (SCM x
, SCM y
)
2023 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
2024 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG1
,
2025 s_scm_centered_remainder
);
2026 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
2027 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
2028 s_scm_centered_remainder
);
2029 else if (scm_is_true (scm_positive_p (y
)))
2030 q
= scm_floor (scm_sum (scm_divide (x
, y
), exactly_one_half
));
2031 else if (scm_is_true (scm_negative_p (y
)))
2032 q
= scm_ceiling (scm_difference (scm_divide (x
, y
), exactly_one_half
));
2034 scm_num_overflow (s_scm_centered_remainder
);
2035 return scm_difference (x
, scm_product (y
, q
));
2039 static SCM
scm_i_inexact_centered_divide (double x
, double y
);
2040 static SCM
scm_i_bigint_centered_divide (SCM x
, SCM y
);
2041 static SCM
scm_i_slow_exact_centered_divide (SCM x
, SCM y
);
2043 SCM_PRIMITIVE_GENERIC (scm_centered_divide
, "centered/", 2, 0, 0,
2045 "Return the integer @var{q} and the real number @var{r}\n"
2046 "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
2047 "and @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
2049 "(centered/ 123 10) @result{} 12 and 3\n"
2050 "(centered/ 123 -10) @result{} -12 and 3\n"
2051 "(centered/ -123 10) @result{} -12 and -3\n"
2052 "(centered/ -123 -10) @result{} 12 and -3\n"
2053 "(centered/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
2054 "(centered/ 16/3 -10/7) @result{} -4 and -8/21\n"
2056 #define FUNC_NAME s_scm_centered_divide
2058 if (SCM_LIKELY (SCM_I_INUMP (x
)))
2060 if (SCM_LIKELY (SCM_I_INUMP (y
)))
2062 scm_t_inum yy
= SCM_I_INUM (y
);
2063 if (SCM_UNLIKELY (yy
== 0))
2064 scm_num_overflow (s_scm_centered_divide
);
2067 scm_t_inum xx
= SCM_I_INUM (x
);
2068 scm_t_inum qq
= xx
/ yy
;
2069 scm_t_inum rr
= xx
- qq
* yy
;
2070 if (SCM_LIKELY (xx
> 0))
2072 if (SCM_LIKELY (yy
> 0))
2074 if (rr
>= (yy
+ 1) / 2)
2079 if (rr
>= (1 - yy
) / 2)
2085 if (SCM_LIKELY (yy
> 0))
2096 return scm_values (scm_list_2 (SCM_I_MAKINUM (qq
),
2097 SCM_I_MAKINUM (rr
)));
2100 else if (SCM_BIGP (y
))
2102 /* Pass a denormalized bignum version of x (even though it
2103 can fit in a fixnum) to scm_i_bigint_centered_divide */
2104 return scm_i_bigint_centered_divide
2105 (scm_i_long2big (SCM_I_INUM (x
)), y
);
2107 else if (SCM_REALP (y
))
2108 return scm_i_inexact_centered_divide
2109 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
2110 else if (SCM_FRACTIONP (y
))
2111 return scm_i_slow_exact_centered_divide (x
, y
);
2113 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2114 s_scm_centered_divide
);
2116 else if (SCM_BIGP (x
))
2118 if (SCM_LIKELY (SCM_I_INUMP (y
)))
2120 scm_t_inum yy
= SCM_I_INUM (y
);
2121 if (SCM_UNLIKELY (yy
== 0))
2122 scm_num_overflow (s_scm_centered_divide
);
2125 SCM q
= scm_i_mkbig ();
2127 /* Arrange for rr to initially be non-positive,
2128 because that simplifies the test to see
2129 if it is within the needed bounds. */
2132 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
2133 SCM_I_BIG_MPZ (x
), yy
);
2134 scm_remember_upto_here_1 (x
);
2137 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
2138 SCM_I_BIG_MPZ (q
), 1);
2144 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
2145 SCM_I_BIG_MPZ (x
), -yy
);
2146 scm_remember_upto_here_1 (x
);
2147 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
2150 mpz_add_ui (SCM_I_BIG_MPZ (q
),
2151 SCM_I_BIG_MPZ (q
), 1);
2155 return scm_values (scm_list_2 (scm_i_normbig (q
),
2156 SCM_I_MAKINUM (rr
)));
2159 else if (SCM_BIGP (y
))
2160 return scm_i_bigint_centered_divide (x
, y
);
2161 else if (SCM_REALP (y
))
2162 return scm_i_inexact_centered_divide
2163 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
2164 else if (SCM_FRACTIONP (y
))
2165 return scm_i_slow_exact_centered_divide (x
, y
);
2167 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2168 s_scm_centered_divide
);
2170 else if (SCM_REALP (x
))
2172 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
2173 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
2174 return scm_i_inexact_centered_divide
2175 (SCM_REAL_VALUE (x
), scm_to_double (y
));
2177 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2178 s_scm_centered_divide
);
2180 else if (SCM_FRACTIONP (x
))
2183 return scm_i_inexact_centered_divide
2184 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
2186 return scm_i_slow_exact_centered_divide (x
, y
);
2189 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG1
,
2190 s_scm_centered_divide
);
2195 scm_i_inexact_centered_divide (double x
, double y
)
2199 if (SCM_LIKELY (y
> 0))
2200 q
= floor (x
/y
+ 0.5);
2201 else if (SCM_LIKELY (y
< 0))
2202 q
= ceil (x
/y
- 0.5);
2204 scm_num_overflow (s_scm_centered_divide
); /* or return a NaN? */
2208 return scm_values (scm_list_2 (scm_from_double (q
),
2209 scm_from_double (r
)));
2212 /* Assumes that both x and y are bigints, though
2213 x might be able to fit into a fixnum. */
2215 scm_i_bigint_centered_divide (SCM x
, SCM y
)
2219 /* Note that x might be small enough to fit into a
2220 fixnum, so we must not let it escape into the wild */
2224 /* min_r will eventually become -abs(y/2) */
2225 min_r
= scm_i_mkbig ();
2226 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
2227 SCM_I_BIG_MPZ (y
), 1);
2229 /* Arrange for rr to initially be non-positive,
2230 because that simplifies the test to see
2231 if it is within the needed bounds. */
2232 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
2234 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
2235 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2236 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
2237 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2239 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
2240 SCM_I_BIG_MPZ (q
), 1);
2241 mpz_add (SCM_I_BIG_MPZ (r
),
2248 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
2249 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2250 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2252 mpz_add_ui (SCM_I_BIG_MPZ (q
),
2253 SCM_I_BIG_MPZ (q
), 1);
2254 mpz_sub (SCM_I_BIG_MPZ (r
),
2259 scm_remember_upto_here_2 (x
, y
);
2260 return scm_values (scm_list_2 (scm_i_normbig (q
),
2261 scm_i_normbig (r
)));
2264 /* Compute exact centered quotient and remainder the slow way.
2265 We use this only if both arguments are exact,
2266 and at least one of them is a fraction */
2268 scm_i_slow_exact_centered_divide (SCM x
, SCM y
)
2272 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
2273 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG1
,
2274 s_scm_centered_divide
);
2275 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
2276 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2277 s_scm_centered_divide
);
2278 else if (scm_is_true (scm_positive_p (y
)))
2279 q
= scm_floor (scm_sum (scm_divide (x
, y
),
2281 else if (scm_is_true (scm_negative_p (y
)))
2282 q
= scm_ceiling (scm_difference (scm_divide (x
, y
),
2285 scm_num_overflow (s_scm_centered_divide
);
2286 r
= scm_difference (x
, scm_product (q
, y
));
2287 return scm_values (scm_list_2 (q
, r
));
2291 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
2292 (SCM x
, SCM y
, SCM rest
),
2293 "Return the greatest common divisor of all parameter values.\n"
2294 "If called without arguments, 0 is returned.")
2295 #define FUNC_NAME s_scm_i_gcd
2297 while (!scm_is_null (rest
))
2298 { x
= scm_gcd (x
, y
);
2300 rest
= scm_cdr (rest
);
2302 return scm_gcd (x
, y
);
2306 #define s_gcd s_scm_i_gcd
2307 #define g_gcd g_scm_i_gcd
2310 scm_gcd (SCM x
, SCM y
)
2313 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
2315 if (SCM_I_INUMP (x
))
2317 if (SCM_I_INUMP (y
))
2319 scm_t_inum xx
= SCM_I_INUM (x
);
2320 scm_t_inum yy
= SCM_I_INUM (y
);
2321 scm_t_inum u
= xx
< 0 ? -xx
: xx
;
2322 scm_t_inum v
= yy
< 0 ? -yy
: yy
;
2332 /* Determine a common factor 2^k */
2333 while (!(1 & (u
| v
)))
2339 /* Now, any factor 2^n can be eliminated */
2359 return (SCM_POSFIXABLE (result
)
2360 ? SCM_I_MAKINUM (result
)
2361 : scm_i_inum2big (result
));
2363 else if (SCM_BIGP (y
))
2369 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
2371 else if (SCM_BIGP (x
))
2373 if (SCM_I_INUMP (y
))
2378 yy
= SCM_I_INUM (y
);
2383 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
2384 scm_remember_upto_here_1 (x
);
2385 return (SCM_POSFIXABLE (result
)
2386 ? SCM_I_MAKINUM (result
)
2387 : scm_from_unsigned_integer (result
));
2389 else if (SCM_BIGP (y
))
2391 SCM result
= scm_i_mkbig ();
2392 mpz_gcd (SCM_I_BIG_MPZ (result
),
2395 scm_remember_upto_here_2 (x
, y
);
2396 return scm_i_normbig (result
);
2399 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
2402 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
2405 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
2406 (SCM x
, SCM y
, SCM rest
),
2407 "Return the least common multiple of the arguments.\n"
2408 "If called without arguments, 1 is returned.")
2409 #define FUNC_NAME s_scm_i_lcm
2411 while (!scm_is_null (rest
))
2412 { x
= scm_lcm (x
, y
);
2414 rest
= scm_cdr (rest
);
2416 return scm_lcm (x
, y
);
2420 #define s_lcm s_scm_i_lcm
2421 #define g_lcm g_scm_i_lcm
2424 scm_lcm (SCM n1
, SCM n2
)
2426 if (SCM_UNBNDP (n2
))
2428 if (SCM_UNBNDP (n1
))
2429 return SCM_I_MAKINUM (1L);
2430 n2
= SCM_I_MAKINUM (1L);
2433 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
2434 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
2435 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
2436 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
2438 if (SCM_I_INUMP (n1
))
2440 if (SCM_I_INUMP (n2
))
2442 SCM d
= scm_gcd (n1
, n2
);
2443 if (scm_is_eq (d
, SCM_INUM0
))
2446 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
2450 /* inum n1, big n2 */
2453 SCM result
= scm_i_mkbig ();
2454 scm_t_inum nn1
= SCM_I_INUM (n1
);
2455 if (nn1
== 0) return SCM_INUM0
;
2456 if (nn1
< 0) nn1
= - nn1
;
2457 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
2458 scm_remember_upto_here_1 (n2
);
2466 if (SCM_I_INUMP (n2
))
2473 SCM result
= scm_i_mkbig ();
2474 mpz_lcm(SCM_I_BIG_MPZ (result
),
2476 SCM_I_BIG_MPZ (n2
));
2477 scm_remember_upto_here_2(n1
, n2
);
2478 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
2484 /* Emulating 2's complement bignums with sign magnitude arithmetic:
2489 + + + x (map digit:logand X Y)
2490 + - + x (map digit:logand X (lognot (+ -1 Y)))
2491 - + + y (map digit:logand (lognot (+ -1 X)) Y)
2492 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
2497 + + + (map digit:logior X Y)
2498 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
2499 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
2500 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
2505 + + + (map digit:logxor X Y)
2506 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
2507 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
2508 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
2513 + + (any digit:logand X Y)
2514 + - (any digit:logand X (lognot (+ -1 Y)))
2515 - + (any digit:logand (lognot (+ -1 X)) Y)
2520 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
2521 (SCM x
, SCM y
, SCM rest
),
2522 "Return the bitwise AND of the integer arguments.\n\n"
2524 "(logand) @result{} -1\n"
2525 "(logand 7) @result{} 7\n"
2526 "(logand #b111 #b011 #b001) @result{} 1\n"
2528 #define FUNC_NAME s_scm_i_logand
2530 while (!scm_is_null (rest
))
2531 { x
= scm_logand (x
, y
);
2533 rest
= scm_cdr (rest
);
2535 return scm_logand (x
, y
);
2539 #define s_scm_logand s_scm_i_logand
2541 SCM
scm_logand (SCM n1
, SCM n2
)
2542 #define FUNC_NAME s_scm_logand
2546 if (SCM_UNBNDP (n2
))
2548 if (SCM_UNBNDP (n1
))
2549 return SCM_I_MAKINUM (-1);
2550 else if (!SCM_NUMBERP (n1
))
2551 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2552 else if (SCM_NUMBERP (n1
))
2555 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2558 if (SCM_I_INUMP (n1
))
2560 nn1
= SCM_I_INUM (n1
);
2561 if (SCM_I_INUMP (n2
))
2563 scm_t_inum nn2
= SCM_I_INUM (n2
);
2564 return SCM_I_MAKINUM (nn1
& nn2
);
2566 else if SCM_BIGP (n2
)
2572 SCM result_z
= scm_i_mkbig ();
2574 mpz_init_set_si (nn1_z
, nn1
);
2575 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2576 scm_remember_upto_here_1 (n2
);
2578 return scm_i_normbig (result_z
);
2582 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2584 else if (SCM_BIGP (n1
))
2586 if (SCM_I_INUMP (n2
))
2589 nn1
= SCM_I_INUM (n1
);
2592 else if (SCM_BIGP (n2
))
2594 SCM result_z
= scm_i_mkbig ();
2595 mpz_and (SCM_I_BIG_MPZ (result_z
),
2597 SCM_I_BIG_MPZ (n2
));
2598 scm_remember_upto_here_2 (n1
, n2
);
2599 return scm_i_normbig (result_z
);
2602 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2605 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2610 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
2611 (SCM x
, SCM y
, SCM rest
),
2612 "Return the bitwise OR of the integer arguments.\n\n"
2614 "(logior) @result{} 0\n"
2615 "(logior 7) @result{} 7\n"
2616 "(logior #b000 #b001 #b011) @result{} 3\n"
2618 #define FUNC_NAME s_scm_i_logior
2620 while (!scm_is_null (rest
))
2621 { x
= scm_logior (x
, y
);
2623 rest
= scm_cdr (rest
);
2625 return scm_logior (x
, y
);
2629 #define s_scm_logior s_scm_i_logior
2631 SCM
scm_logior (SCM n1
, SCM n2
)
2632 #define FUNC_NAME s_scm_logior
2636 if (SCM_UNBNDP (n2
))
2638 if (SCM_UNBNDP (n1
))
2640 else if (SCM_NUMBERP (n1
))
2643 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2646 if (SCM_I_INUMP (n1
))
2648 nn1
= SCM_I_INUM (n1
);
2649 if (SCM_I_INUMP (n2
))
2651 long nn2
= SCM_I_INUM (n2
);
2652 return SCM_I_MAKINUM (nn1
| nn2
);
2654 else if (SCM_BIGP (n2
))
2660 SCM result_z
= scm_i_mkbig ();
2662 mpz_init_set_si (nn1_z
, nn1
);
2663 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2664 scm_remember_upto_here_1 (n2
);
2666 return scm_i_normbig (result_z
);
2670 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2672 else if (SCM_BIGP (n1
))
2674 if (SCM_I_INUMP (n2
))
2677 nn1
= SCM_I_INUM (n1
);
2680 else if (SCM_BIGP (n2
))
2682 SCM result_z
= scm_i_mkbig ();
2683 mpz_ior (SCM_I_BIG_MPZ (result_z
),
2685 SCM_I_BIG_MPZ (n2
));
2686 scm_remember_upto_here_2 (n1
, n2
);
2687 return scm_i_normbig (result_z
);
2690 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2693 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2698 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
2699 (SCM x
, SCM y
, SCM rest
),
2700 "Return the bitwise XOR of the integer arguments. A bit is\n"
2701 "set in the result if it is set in an odd number of arguments.\n"
2703 "(logxor) @result{} 0\n"
2704 "(logxor 7) @result{} 7\n"
2705 "(logxor #b000 #b001 #b011) @result{} 2\n"
2706 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
2708 #define FUNC_NAME s_scm_i_logxor
2710 while (!scm_is_null (rest
))
2711 { x
= scm_logxor (x
, y
);
2713 rest
= scm_cdr (rest
);
2715 return scm_logxor (x
, y
);
2719 #define s_scm_logxor s_scm_i_logxor
2721 SCM
scm_logxor (SCM n1
, SCM n2
)
2722 #define FUNC_NAME s_scm_logxor
2726 if (SCM_UNBNDP (n2
))
2728 if (SCM_UNBNDP (n1
))
2730 else if (SCM_NUMBERP (n1
))
2733 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2736 if (SCM_I_INUMP (n1
))
2738 nn1
= SCM_I_INUM (n1
);
2739 if (SCM_I_INUMP (n2
))
2741 scm_t_inum nn2
= SCM_I_INUM (n2
);
2742 return SCM_I_MAKINUM (nn1
^ nn2
);
2744 else if (SCM_BIGP (n2
))
2748 SCM result_z
= scm_i_mkbig ();
2750 mpz_init_set_si (nn1_z
, nn1
);
2751 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2752 scm_remember_upto_here_1 (n2
);
2754 return scm_i_normbig (result_z
);
2758 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2760 else if (SCM_BIGP (n1
))
2762 if (SCM_I_INUMP (n2
))
2765 nn1
= SCM_I_INUM (n1
);
2768 else if (SCM_BIGP (n2
))
2770 SCM result_z
= scm_i_mkbig ();
2771 mpz_xor (SCM_I_BIG_MPZ (result_z
),
2773 SCM_I_BIG_MPZ (n2
));
2774 scm_remember_upto_here_2 (n1
, n2
);
2775 return scm_i_normbig (result_z
);
2778 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2781 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2786 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
2788 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
2789 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
2790 "without actually calculating the @code{logand}, just testing\n"
2794 "(logtest #b0100 #b1011) @result{} #f\n"
2795 "(logtest #b0100 #b0111) @result{} #t\n"
2797 #define FUNC_NAME s_scm_logtest
2801 if (SCM_I_INUMP (j
))
2803 nj
= SCM_I_INUM (j
);
2804 if (SCM_I_INUMP (k
))
2806 scm_t_inum nk
= SCM_I_INUM (k
);
2807 return scm_from_bool (nj
& nk
);
2809 else if (SCM_BIGP (k
))
2817 mpz_init_set_si (nj_z
, nj
);
2818 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
2819 scm_remember_upto_here_1 (k
);
2820 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
2826 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
2828 else if (SCM_BIGP (j
))
2830 if (SCM_I_INUMP (k
))
2833 nj
= SCM_I_INUM (j
);
2836 else if (SCM_BIGP (k
))
2840 mpz_init (result_z
);
2844 scm_remember_upto_here_2 (j
, k
);
2845 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
2846 mpz_clear (result_z
);
2850 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
2853 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
2858 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
2860 "Test whether bit number @var{index} in @var{j} is set.\n"
2861 "@var{index} starts from 0 for the least significant bit.\n"
2864 "(logbit? 0 #b1101) @result{} #t\n"
2865 "(logbit? 1 #b1101) @result{} #f\n"
2866 "(logbit? 2 #b1101) @result{} #t\n"
2867 "(logbit? 3 #b1101) @result{} #t\n"
2868 "(logbit? 4 #b1101) @result{} #f\n"
2870 #define FUNC_NAME s_scm_logbit_p
2872 unsigned long int iindex
;
2873 iindex
= scm_to_ulong (index
);
2875 if (SCM_I_INUMP (j
))
2877 /* bits above what's in an inum follow the sign bit */
2878 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
2879 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
2881 else if (SCM_BIGP (j
))
2883 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
2884 scm_remember_upto_here_1 (j
);
2885 return scm_from_bool (val
);
2888 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
2893 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
2895 "Return the integer which is the ones-complement of the integer\n"
2899 "(number->string (lognot #b10000000) 2)\n"
2900 " @result{} \"-10000001\"\n"
2901 "(number->string (lognot #b0) 2)\n"
2902 " @result{} \"-1\"\n"
2904 #define FUNC_NAME s_scm_lognot
2906 if (SCM_I_INUMP (n
)) {
2907 /* No overflow here, just need to toggle all the bits making up the inum.
2908 Enhancement: No need to strip the tag and add it back, could just xor
2909 a block of 1 bits, if that worked with the various debug versions of
2911 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
2913 } else if (SCM_BIGP (n
)) {
2914 SCM result
= scm_i_mkbig ();
2915 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
2916 scm_remember_upto_here_1 (n
);
2920 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2925 /* returns 0 if IN is not an integer. OUT must already be
2928 coerce_to_big (SCM in
, mpz_t out
)
2931 mpz_set (out
, SCM_I_BIG_MPZ (in
));
2932 else if (SCM_I_INUMP (in
))
2933 mpz_set_si (out
, SCM_I_INUM (in
));
2940 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
2941 (SCM n
, SCM k
, SCM m
),
2942 "Return @var{n} raised to the integer exponent\n"
2943 "@var{k}, modulo @var{m}.\n"
2946 "(modulo-expt 2 3 5)\n"
2949 #define FUNC_NAME s_scm_modulo_expt
2955 /* There are two classes of error we might encounter --
2956 1) Math errors, which we'll report by calling scm_num_overflow,
2958 2) wrong-type errors, which of course we'll report by calling
2960 We don't report those errors immediately, however; instead we do
2961 some cleanup first. These variables tell us which error (if
2962 any) we should report after cleaning up.
2964 int report_overflow
= 0;
2966 int position_of_wrong_type
= 0;
2967 SCM value_of_wrong_type
= SCM_INUM0
;
2969 SCM result
= SCM_UNDEFINED
;
2975 if (scm_is_eq (m
, SCM_INUM0
))
2977 report_overflow
= 1;
2981 if (!coerce_to_big (n
, n_tmp
))
2983 value_of_wrong_type
= n
;
2984 position_of_wrong_type
= 1;
2988 if (!coerce_to_big (k
, k_tmp
))
2990 value_of_wrong_type
= k
;
2991 position_of_wrong_type
= 2;
2995 if (!coerce_to_big (m
, m_tmp
))
2997 value_of_wrong_type
= m
;
2998 position_of_wrong_type
= 3;
3002 /* if the exponent K is negative, and we simply call mpz_powm, we
3003 will get a divide-by-zero exception when an inverse 1/n mod m
3004 doesn't exist (or is not unique). Since exceptions are hard to
3005 handle, we'll attempt the inversion "by hand" -- that way, we get
3006 a simple failure code, which is easy to handle. */
3008 if (-1 == mpz_sgn (k_tmp
))
3010 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
3012 report_overflow
= 1;
3015 mpz_neg (k_tmp
, k_tmp
);
3018 result
= scm_i_mkbig ();
3019 mpz_powm (SCM_I_BIG_MPZ (result
),
3024 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
3025 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
3032 if (report_overflow
)
3033 scm_num_overflow (FUNC_NAME
);
3035 if (position_of_wrong_type
)
3036 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
3037 value_of_wrong_type
);
3039 return scm_i_normbig (result
);
3043 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
3045 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
3046 "exact integer, @var{n} can be any number.\n"
3048 "Negative @var{k} is supported, and results in\n"
3049 "@math{1/@var{n}^abs(@var{k})} in the usual way.\n"
3050 "@math{@var{n}^0} is 1, as usual, and that\n"
3051 "includes @math{0^0} is 1.\n"
3054 "(integer-expt 2 5) @result{} 32\n"
3055 "(integer-expt -3 3) @result{} -27\n"
3056 "(integer-expt 5 -3) @result{} 1/125\n"
3057 "(integer-expt 0 0) @result{} 1\n"
3059 #define FUNC_NAME s_scm_integer_expt
3062 SCM z_i2
= SCM_BOOL_F
;
3064 SCM acc
= SCM_I_MAKINUM (1L);
3066 /* Specifically refrain from checking the type of the first argument.
3067 This allows us to exponentiate any object that can be multiplied.
3068 If we must raise to a negative power, we must also be able to
3069 take its reciprocal. */
3070 if (!SCM_LIKELY (SCM_I_INUMP (k
)) && !SCM_LIKELY (SCM_BIGP (k
)))
3071 SCM_WRONG_TYPE_ARG (2, k
);
3073 if (SCM_UNLIKELY (scm_is_eq (k
, SCM_INUM0
)))
3074 return SCM_INUM1
; /* n^(exact0) is exact 1, regardless of n */
3075 else if (SCM_UNLIKELY (scm_is_eq (n
, SCM_I_MAKINUM (-1L))))
3076 return scm_is_false (scm_even_p (k
)) ? n
: SCM_INUM1
;
3077 /* The next check is necessary only because R6RS specifies different
3078 behavior for 0^(-k) than for (/ 0). If n is not a scheme number,
3079 we simply skip this case and move on. */
3080 else if (SCM_NUMBERP (n
) && scm_is_true (scm_zero_p (n
)))
3082 /* k cannot be 0 at this point, because we
3083 have already checked for that case above */
3084 if (scm_is_true (scm_positive_p (k
)))
3086 else /* return NaN for (0 ^ k) for negative k per R6RS */
3090 if (SCM_I_INUMP (k
))
3091 i2
= SCM_I_INUM (k
);
3092 else if (SCM_BIGP (k
))
3094 z_i2
= scm_i_clonebig (k
, 1);
3095 scm_remember_upto_here_1 (k
);
3099 SCM_WRONG_TYPE_ARG (2, k
);
3103 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
3105 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
3106 n
= scm_divide (n
, SCM_UNDEFINED
);
3110 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
3114 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
3116 return scm_product (acc
, n
);
3118 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
3119 acc
= scm_product (acc
, n
);
3120 n
= scm_product (n
, n
);
3121 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
3129 n
= scm_divide (n
, SCM_UNDEFINED
);
3136 return scm_product (acc
, n
);
3138 acc
= scm_product (acc
, n
);
3139 n
= scm_product (n
, n
);
3146 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
3148 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
3149 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
3151 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
3152 "@var{cnt} is negative it's a division, rounded towards negative\n"
3153 "infinity. (Note that this is not the same rounding as\n"
3154 "@code{quotient} does.)\n"
3156 "With @var{n} viewed as an infinite precision twos complement,\n"
3157 "@code{ash} means a left shift introducing zero bits, or a right\n"
3158 "shift dropping bits.\n"
3161 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
3162 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
3164 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
3165 "(ash -23 -2) @result{} -6\n"
3167 #define FUNC_NAME s_scm_ash
3170 bits_to_shift
= scm_to_long (cnt
);
3172 if (SCM_I_INUMP (n
))
3174 scm_t_inum nn
= SCM_I_INUM (n
);
3176 if (bits_to_shift
> 0)
3178 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
3179 overflow a non-zero fixnum. For smaller shifts we check the
3180 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
3181 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
3182 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
3188 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
3190 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
3193 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
3197 SCM result
= scm_i_inum2big (nn
);
3198 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
3205 bits_to_shift
= -bits_to_shift
;
3206 if (bits_to_shift
>= SCM_LONG_BIT
)
3207 return (nn
>= 0 ? SCM_INUM0
: SCM_I_MAKINUM(-1));
3209 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
3213 else if (SCM_BIGP (n
))
3217 if (bits_to_shift
== 0)
3220 result
= scm_i_mkbig ();
3221 if (bits_to_shift
>= 0)
3223 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
3229 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
3230 we have to allocate a bignum even if the result is going to be a
3232 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
3234 return scm_i_normbig (result
);
3240 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3246 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
3247 (SCM n
, SCM start
, SCM end
),
3248 "Return the integer composed of the @var{start} (inclusive)\n"
3249 "through @var{end} (exclusive) bits of @var{n}. The\n"
3250 "@var{start}th bit becomes the 0-th bit in the result.\n"
3253 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
3254 " @result{} \"1010\"\n"
3255 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
3256 " @result{} \"10110\"\n"
3258 #define FUNC_NAME s_scm_bit_extract
3260 unsigned long int istart
, iend
, bits
;
3261 istart
= scm_to_ulong (start
);
3262 iend
= scm_to_ulong (end
);
3263 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
3265 /* how many bits to keep */
3266 bits
= iend
- istart
;
3268 if (SCM_I_INUMP (n
))
3270 scm_t_inum in
= SCM_I_INUM (n
);
3272 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
3273 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
3274 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
3276 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
3278 /* Since we emulate two's complement encoded numbers, this
3279 * special case requires us to produce a result that has
3280 * more bits than can be stored in a fixnum.
3282 SCM result
= scm_i_inum2big (in
);
3283 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
3288 /* mask down to requisite bits */
3289 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
3290 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
3292 else if (SCM_BIGP (n
))
3297 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
3301 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
3302 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
3303 such bits into a ulong. */
3304 result
= scm_i_mkbig ();
3305 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
3306 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
3307 result
= scm_i_normbig (result
);
3309 scm_remember_upto_here_1 (n
);
3313 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3318 static const char scm_logtab
[] = {
3319 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
3322 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
3324 "Return the number of bits in integer @var{n}. If integer is\n"
3325 "positive, the 1-bits in its binary representation are counted.\n"
3326 "If negative, the 0-bits in its two's-complement binary\n"
3327 "representation are counted. If 0, 0 is returned.\n"
3330 "(logcount #b10101010)\n"
3337 #define FUNC_NAME s_scm_logcount
3339 if (SCM_I_INUMP (n
))
3341 unsigned long c
= 0;
3342 scm_t_inum nn
= SCM_I_INUM (n
);
3347 c
+= scm_logtab
[15 & nn
];
3350 return SCM_I_MAKINUM (c
);
3352 else if (SCM_BIGP (n
))
3354 unsigned long count
;
3355 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
3356 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
3358 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
3359 scm_remember_upto_here_1 (n
);
3360 return SCM_I_MAKINUM (count
);
3363 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3368 static const char scm_ilentab
[] = {
3369 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
3373 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
3375 "Return the number of bits necessary to represent @var{n}.\n"
3378 "(integer-length #b10101010)\n"
3380 "(integer-length 0)\n"
3382 "(integer-length #b1111)\n"
3385 #define FUNC_NAME s_scm_integer_length
3387 if (SCM_I_INUMP (n
))
3389 unsigned long c
= 0;
3391 scm_t_inum nn
= SCM_I_INUM (n
);
3397 l
= scm_ilentab
[15 & nn
];
3400 return SCM_I_MAKINUM (c
- 4 + l
);
3402 else if (SCM_BIGP (n
))
3404 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
3405 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
3406 1 too big, so check for that and adjust. */
3407 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
3408 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
3409 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
3410 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
3412 scm_remember_upto_here_1 (n
);
3413 return SCM_I_MAKINUM (size
);
3416 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3420 /*** NUMBERS -> STRINGS ***/
3421 #define SCM_MAX_DBL_PREC 60
3422 #define SCM_MAX_DBL_RADIX 36
3424 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
3425 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
3426 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
3429 void init_dblprec(int *prec
, int radix
) {
3430 /* determine floating point precision by adding successively
3431 smaller increments to 1.0 until it is considered == 1.0 */
3432 double f
= ((double)1.0)/radix
;
3433 double fsum
= 1.0 + f
;
3438 if (++(*prec
) > SCM_MAX_DBL_PREC
)
3450 void init_fx_radix(double *fx_list
, int radix
)
3452 /* initialize a per-radix list of tolerances. When added
3453 to a number < 1.0, we can determine if we should raund
3454 up and quit converting a number to a string. */
3458 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
3459 fx_list
[i
] = (fx_list
[i
-1] / radix
);
3462 /* use this array as a way to generate a single digit */
3463 static const char number_chars
[] = "0123456789abcdefghijklmnopqrstuvwxyz";
3466 idbl2str (double f
, char *a
, int radix
)
3468 int efmt
, dpt
, d
, i
, wp
;
3470 #ifdef DBL_MIN_10_EXP
3473 #endif /* DBL_MIN_10_EXP */
3478 radix
> SCM_MAX_DBL_RADIX
)
3480 /* revert to existing behavior */
3484 wp
= scm_dblprec
[radix
-2];
3485 fx
= fx_per_radix
[radix
-2];
3489 #ifdef HAVE_COPYSIGN
3490 double sgn
= copysign (1.0, f
);
3495 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
3501 strcpy (a
, "-inf.0");
3503 strcpy (a
, "+inf.0");
3508 strcpy (a
, "+nan.0");
3518 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
3519 make-uniform-vector, from causing infinite loops. */
3520 /* just do the checking...if it passes, we do the conversion for our
3521 radix again below */
3528 if (exp_cpy
-- < DBL_MIN_10_EXP
)
3536 while (f_cpy
> 10.0)
3539 if (exp_cpy
++ > DBL_MAX_10_EXP
)
3560 if (f
+ fx
[wp
] >= radix
)
3567 /* adding 9999 makes this equivalent to abs(x) % 3 */
3568 dpt
= (exp
+ 9999) % 3;
3572 efmt
= (exp
< -3) || (exp
> wp
+ 2);
3594 a
[ch
++] = number_chars
[d
];
3597 if (f
+ fx
[wp
] >= 1.0)
3599 a
[ch
- 1] = number_chars
[d
+1];
3611 if ((dpt
> 4) && (exp
> 6))
3613 d
= (a
[0] == '-' ? 2 : 1);
3614 for (i
= ch
++; i
> d
; i
--)
3627 if (a
[ch
- 1] == '.')
3628 a
[ch
++] = '0'; /* trailing zero */
3637 for (i
= radix
; i
<= exp
; i
*= radix
);
3638 for (i
/= radix
; i
; i
/= radix
)
3640 a
[ch
++] = number_chars
[exp
/ i
];
3649 icmplx2str (double real
, double imag
, char *str
, int radix
)
3653 i
= idbl2str (real
, str
, radix
);
3656 /* Don't output a '+' for negative numbers or for Inf and
3657 NaN. They will provide their own sign. */
3658 if (0 <= imag
&& !isinf (imag
) && !isnan (imag
))
3660 i
+= idbl2str (imag
, &str
[i
], radix
);
3667 iflo2str (SCM flt
, char *str
, int radix
)
3670 if (SCM_REALP (flt
))
3671 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
3673 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
3678 /* convert a scm_t_intmax to a string (unterminated). returns the number of
3679 characters in the result.
3681 p is destination: worst case (base 2) is SCM_INTBUFLEN */
3683 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
3688 return scm_iuint2str (-num
, rad
, p
) + 1;
3691 return scm_iuint2str (num
, rad
, p
);
3694 /* convert a scm_t_intmax to a string (unterminated). returns the number of
3695 characters in the result.
3697 p is destination: worst case (base 2) is SCM_INTBUFLEN */
3699 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
3703 scm_t_uintmax n
= num
;
3705 if (rad
< 2 || rad
> 36)
3706 scm_out_of_range ("scm_iuint2str", scm_from_int (rad
));
3708 for (n
/= rad
; n
> 0; n
/= rad
)
3718 p
[i
] = number_chars
[d
];
3723 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
3725 "Return a string holding the external representation of the\n"
3726 "number @var{n} in the given @var{radix}. If @var{n} is\n"
3727 "inexact, a radix of 10 will be used.")
3728 #define FUNC_NAME s_scm_number_to_string
3732 if (SCM_UNBNDP (radix
))
3735 base
= scm_to_signed_integer (radix
, 2, 36);
3737 if (SCM_I_INUMP (n
))
3739 char num_buf
[SCM_INTBUFLEN
];
3740 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
3741 return scm_from_locale_stringn (num_buf
, length
);
3743 else if (SCM_BIGP (n
))
3745 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
3746 scm_remember_upto_here_1 (n
);
3747 return scm_take_locale_string (str
);
3749 else if (SCM_FRACTIONP (n
))
3751 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
3752 scm_from_locale_string ("/"),
3753 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
3755 else if (SCM_INEXACTP (n
))
3757 char num_buf
[FLOBUFLEN
];
3758 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
3761 SCM_WRONG_TYPE_ARG (1, n
);
3766 /* These print routines used to be stubbed here so that scm_repl.c
3767 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
3770 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3772 char num_buf
[FLOBUFLEN
];
3773 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
3778 scm_i_print_double (double val
, SCM port
)
3780 char num_buf
[FLOBUFLEN
];
3781 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
3785 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3788 char num_buf
[FLOBUFLEN
];
3789 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
3794 scm_i_print_complex (double real
, double imag
, SCM port
)
3796 char num_buf
[FLOBUFLEN
];
3797 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
3801 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3804 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
3805 scm_display (str
, port
);
3806 scm_remember_upto_here_1 (str
);
3811 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3813 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
3814 scm_remember_upto_here_1 (exp
);
3815 scm_lfwrite (str
, (size_t) strlen (str
), port
);
3819 /*** END nums->strs ***/
3822 /*** STRINGS -> NUMBERS ***/
3824 /* The following functions implement the conversion from strings to numbers.
3825 * The implementation somehow follows the grammar for numbers as it is given
3826 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
3827 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
3828 * points should be noted about the implementation:
3829 * * Each function keeps a local index variable 'idx' that points at the
3830 * current position within the parsed string. The global index is only
3831 * updated if the function could parse the corresponding syntactic unit
3833 * * Similarly, the functions keep track of indicators of inexactness ('#',
3834 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
3835 * global exactness information is only updated after each part has been
3836 * successfully parsed.
3837 * * Sequences of digits are parsed into temporary variables holding fixnums.
3838 * Only if these fixnums would overflow, the result variables are updated
3839 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
3840 * the temporary variables holding the fixnums are cleared, and the process
3841 * starts over again. If for example fixnums were able to store five decimal
3842 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
3843 * and the result was computed as 12345 * 100000 + 67890. In other words,
3844 * only every five digits two bignum operations were performed.
3847 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
3849 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
3851 /* Caller is responsible for checking that the return value is in range
3852 for the given radix, which should be <= 36. */
3854 char_decimal_value (scm_t_uint32 c
)
3856 /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
3857 that's certainly above any valid decimal, so we take advantage of
3858 that to elide some tests. */
3859 unsigned int d
= (unsigned int) uc_decimal_value (c
);
3861 /* If that failed, try extended hexadecimals, then. Only accept ascii
3866 if (c
>= (scm_t_uint32
) 'a')
3867 d
= c
- (scm_t_uint32
)'a' + 10U;
3873 mem2uinteger (SCM mem
, unsigned int *p_idx
,
3874 unsigned int radix
, enum t_exactness
*p_exactness
)
3876 unsigned int idx
= *p_idx
;
3877 unsigned int hash_seen
= 0;
3878 scm_t_bits shift
= 1;
3880 unsigned int digit_value
;
3883 size_t len
= scm_i_string_length (mem
);
3888 c
= scm_i_string_ref (mem
, idx
);
3889 digit_value
= char_decimal_value (c
);
3890 if (digit_value
>= radix
)
3894 result
= SCM_I_MAKINUM (digit_value
);
3897 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
3907 digit_value
= char_decimal_value (c
);
3908 /* This check catches non-decimals in addition to out-of-range
3910 if (digit_value
>= radix
)
3915 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
3917 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3919 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3926 shift
= shift
* radix
;
3927 add
= add
* radix
+ digit_value
;
3932 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3934 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3938 *p_exactness
= INEXACT
;
3944 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
3945 * covers the parts of the rules that start at a potential point. The value
3946 * of the digits up to the point have been parsed by the caller and are given
3947 * in variable result. The content of *p_exactness indicates, whether a hash
3948 * has already been seen in the digits before the point.
3951 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
3954 mem2decimal_from_point (SCM result
, SCM mem
,
3955 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
3957 unsigned int idx
= *p_idx
;
3958 enum t_exactness x
= *p_exactness
;
3959 size_t len
= scm_i_string_length (mem
);
3964 if (scm_i_string_ref (mem
, idx
) == '.')
3966 scm_t_bits shift
= 1;
3968 unsigned int digit_value
;
3969 SCM big_shift
= SCM_INUM1
;
3974 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
3975 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
3980 digit_value
= DIGIT2UINT (c
);
3991 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
3993 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
3994 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3996 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
4004 add
= add
* 10 + digit_value
;
4010 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
4011 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
4012 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
4015 result
= scm_divide (result
, big_shift
);
4017 /* We've seen a decimal point, thus the value is implicitly inexact. */
4029 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
4031 switch (scm_i_string_ref (mem
, idx
))
4043 c
= scm_i_string_ref (mem
, idx
);
4051 c
= scm_i_string_ref (mem
, idx
);
4060 c
= scm_i_string_ref (mem
, idx
);
4065 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4069 exponent
= DIGIT2UINT (c
);
4072 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
4073 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4076 if (exponent
<= SCM_MAXEXP
)
4077 exponent
= exponent
* 10 + DIGIT2UINT (c
);
4083 if (exponent
> SCM_MAXEXP
)
4085 size_t exp_len
= idx
- start
;
4086 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
4087 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
4088 scm_out_of_range ("string->number", exp_num
);
4091 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
4093 result
= scm_product (result
, e
);
4095 result
= scm_divide2real (result
, e
);
4097 /* We've seen an exponent, thus the value is implicitly inexact. */
4115 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
4118 mem2ureal (SCM mem
, unsigned int *p_idx
,
4119 unsigned int radix
, enum t_exactness
*p_exactness
)
4121 unsigned int idx
= *p_idx
;
4123 size_t len
= scm_i_string_length (mem
);
4125 /* Start off believing that the number will be exact. This changes
4126 to INEXACT if we see a decimal point or a hash. */
4127 enum t_exactness x
= EXACT
;
4132 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
4138 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
4140 /* Cobble up the fractional part. We might want to set the
4141 NaN's mantissa from it. */
4143 mem2uinteger (mem
, &idx
, 10, &x
);
4148 if (scm_i_string_ref (mem
, idx
) == '.')
4152 else if (idx
+ 1 == len
)
4154 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
4157 result
= mem2decimal_from_point (SCM_INUM0
, mem
,
4164 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
4165 if (scm_is_false (uinteger
))
4170 else if (scm_i_string_ref (mem
, idx
) == '/')
4178 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
4179 if (scm_is_false (divisor
))
4182 /* both are int/big here, I assume */
4183 result
= scm_i_make_ratio (uinteger
, divisor
);
4185 else if (radix
== 10)
4187 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
4188 if (scm_is_false (result
))
4197 /* Update *p_exactness if the number just read was inexact. This is
4198 important for complex numbers, so that a complex number is
4199 treated as inexact overall if either its real or imaginary part
4205 /* When returning an inexact zero, make sure it is represented as a
4206 floating point value so that we can change its sign.
4208 if (scm_is_eq (result
, SCM_INUM0
) && *p_exactness
== INEXACT
)
4209 result
= scm_from_double (0.0);
4215 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
4218 mem2complex (SCM mem
, unsigned int idx
,
4219 unsigned int radix
, enum t_exactness
*p_exactness
)
4224 size_t len
= scm_i_string_length (mem
);
4229 c
= scm_i_string_ref (mem
, idx
);
4244 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
4245 if (scm_is_false (ureal
))
4247 /* input must be either +i or -i */
4252 if (scm_i_string_ref (mem
, idx
) == 'i'
4253 || scm_i_string_ref (mem
, idx
) == 'I')
4259 return scm_make_rectangular (SCM_INUM0
, SCM_I_MAKINUM (sign
));
4266 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
4267 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
4272 c
= scm_i_string_ref (mem
, idx
);
4276 /* either +<ureal>i or -<ureal>i */
4283 return scm_make_rectangular (SCM_INUM0
, ureal
);
4286 /* polar input: <real>@<real>. */
4297 c
= scm_i_string_ref (mem
, idx
);
4315 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
4316 if (scm_is_false (angle
))
4321 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
4322 angle
= scm_difference (angle
, SCM_UNDEFINED
);
4324 result
= scm_make_polar (ureal
, angle
);
4329 /* expecting input matching <real>[+-]<ureal>?i */
4336 int sign
= (c
== '+') ? 1 : -1;
4337 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
4339 if (scm_is_false (imag
))
4340 imag
= SCM_I_MAKINUM (sign
);
4341 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
4342 imag
= scm_difference (imag
, SCM_UNDEFINED
);
4346 if (scm_i_string_ref (mem
, idx
) != 'i'
4347 && scm_i_string_ref (mem
, idx
) != 'I')
4354 return scm_make_rectangular (ureal
, imag
);
4363 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
4365 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
4368 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
4370 unsigned int idx
= 0;
4371 unsigned int radix
= NO_RADIX
;
4372 enum t_exactness forced_x
= NO_EXACTNESS
;
4373 enum t_exactness implicit_x
= EXACT
;
4375 size_t len
= scm_i_string_length (mem
);
4377 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
4378 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
4380 switch (scm_i_string_ref (mem
, idx
+ 1))
4383 if (radix
!= NO_RADIX
)
4388 if (radix
!= NO_RADIX
)
4393 if (forced_x
!= NO_EXACTNESS
)
4398 if (forced_x
!= NO_EXACTNESS
)
4403 if (radix
!= NO_RADIX
)
4408 if (radix
!= NO_RADIX
)
4418 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
4419 if (radix
== NO_RADIX
)
4420 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
4422 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
4424 if (scm_is_false (result
))
4430 if (SCM_INEXACTP (result
))
4431 return scm_inexact_to_exact (result
);
4435 if (SCM_INEXACTP (result
))
4438 return scm_exact_to_inexact (result
);
4441 if (implicit_x
== INEXACT
)
4443 if (SCM_INEXACTP (result
))
4446 return scm_exact_to_inexact (result
);
4454 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
4455 unsigned int default_radix
)
4457 SCM str
= scm_from_locale_stringn (mem
, len
);
4459 return scm_i_string_to_number (str
, default_radix
);
4463 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
4464 (SCM string
, SCM radix
),
4465 "Return a number of the maximally precise representation\n"
4466 "expressed by the given @var{string}. @var{radix} must be an\n"
4467 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
4468 "is a default radix that may be overridden by an explicit radix\n"
4469 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
4470 "supplied, then the default radix is 10. If string is not a\n"
4471 "syntactically valid notation for a number, then\n"
4472 "@code{string->number} returns @code{#f}.")
4473 #define FUNC_NAME s_scm_string_to_number
4477 SCM_VALIDATE_STRING (1, string
);
4479 if (SCM_UNBNDP (radix
))
4482 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
4484 answer
= scm_i_string_to_number (string
, base
);
4485 scm_remember_upto_here_1 (string
);
4491 /*** END strs->nums ***/
4494 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
4496 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
4498 #define FUNC_NAME s_scm_number_p
4500 return scm_from_bool (SCM_NUMBERP (x
));
4504 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
4506 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
4507 "otherwise. Note that the sets of real, rational and integer\n"
4508 "values form subsets of the set of complex numbers, i. e. the\n"
4509 "predicate will also be fulfilled if @var{x} is a real,\n"
4510 "rational or integer number.")
4511 #define FUNC_NAME s_scm_complex_p
4513 /* all numbers are complex. */
4514 return scm_number_p (x
);
4518 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
4520 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
4521 "otherwise. Note that the set of integer values forms a subset of\n"
4522 "the set of real numbers, i. e. the predicate will also be\n"
4523 "fulfilled if @var{x} is an integer number.")
4524 #define FUNC_NAME s_scm_real_p
4526 return scm_from_bool
4527 (SCM_I_INUMP (x
) || SCM_REALP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
));
4531 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
4533 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
4534 "otherwise. Note that the set of integer values forms a subset of\n"
4535 "the set of rational numbers, i. e. the predicate will also be\n"
4536 "fulfilled if @var{x} is an integer number.")
4537 #define FUNC_NAME s_scm_rational_p
4539 if (SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
))
4541 else if (SCM_REALP (x
))
4542 /* due to their limited precision, finite floating point numbers are
4543 rational as well. (finite means neither infinity nor a NaN) */
4544 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
4550 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
4552 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
4554 #define FUNC_NAME s_scm_integer_p
4556 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4558 else if (SCM_REALP (x
))
4560 double val
= SCM_REAL_VALUE (x
);
4561 return scm_from_bool (!isinf (val
) && (val
== floor (val
)));
4569 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
4570 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
4571 (SCM x
, SCM y
, SCM rest
),
4572 "Return @code{#t} if all parameters are numerically equal.")
4573 #define FUNC_NAME s_scm_i_num_eq_p
4575 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4577 while (!scm_is_null (rest
))
4579 if (scm_is_false (scm_num_eq_p (x
, y
)))
4583 rest
= scm_cdr (rest
);
4585 return scm_num_eq_p (x
, y
);
4589 scm_num_eq_p (SCM x
, SCM y
)
4592 if (SCM_I_INUMP (x
))
4594 scm_t_signed_bits xx
= SCM_I_INUM (x
);
4595 if (SCM_I_INUMP (y
))
4597 scm_t_signed_bits yy
= SCM_I_INUM (y
);
4598 return scm_from_bool (xx
== yy
);
4600 else if (SCM_BIGP (y
))
4602 else if (SCM_REALP (y
))
4604 /* On a 32-bit system an inum fits a double, we can cast the inum
4605 to a double and compare.
4607 But on a 64-bit system an inum is bigger than a double and
4608 casting it to a double (call that dxx) will round. dxx is at
4609 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
4610 an integer and fits a long. So we cast yy to a long and
4611 compare with plain xx.
4613 An alternative (for any size system actually) would be to check
4614 yy is an integer (with floor) and is in range of an inum
4615 (compare against appropriate powers of 2) then test
4616 xx==(scm_t_signed_bits)yy. It's just a matter of which
4617 casts/comparisons might be fastest or easiest for the cpu. */
4619 double yy
= SCM_REAL_VALUE (y
);
4620 return scm_from_bool ((double) xx
== yy
4621 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
4622 || xx
== (scm_t_signed_bits
) yy
));
4624 else if (SCM_COMPLEXP (y
))
4625 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
4626 && (0.0 == SCM_COMPLEX_IMAG (y
)));
4627 else if (SCM_FRACTIONP (y
))
4630 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4632 else if (SCM_BIGP (x
))
4634 if (SCM_I_INUMP (y
))
4636 else if (SCM_BIGP (y
))
4638 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4639 scm_remember_upto_here_2 (x
, y
);
4640 return scm_from_bool (0 == cmp
);
4642 else if (SCM_REALP (y
))
4645 if (isnan (SCM_REAL_VALUE (y
)))
4647 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
4648 scm_remember_upto_here_1 (x
);
4649 return scm_from_bool (0 == cmp
);
4651 else if (SCM_COMPLEXP (y
))
4654 if (0.0 != SCM_COMPLEX_IMAG (y
))
4656 if (isnan (SCM_COMPLEX_REAL (y
)))
4658 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
4659 scm_remember_upto_here_1 (x
);
4660 return scm_from_bool (0 == cmp
);
4662 else if (SCM_FRACTIONP (y
))
4665 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4667 else if (SCM_REALP (x
))
4669 double xx
= SCM_REAL_VALUE (x
);
4670 if (SCM_I_INUMP (y
))
4672 /* see comments with inum/real above */
4673 scm_t_signed_bits yy
= SCM_I_INUM (y
);
4674 return scm_from_bool (xx
== (double) yy
4675 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
4676 || (scm_t_signed_bits
) xx
== yy
));
4678 else if (SCM_BIGP (y
))
4681 if (isnan (SCM_REAL_VALUE (x
)))
4683 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
4684 scm_remember_upto_here_1 (y
);
4685 return scm_from_bool (0 == cmp
);
4687 else if (SCM_REALP (y
))
4688 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
4689 else if (SCM_COMPLEXP (y
))
4690 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
4691 && (0.0 == SCM_COMPLEX_IMAG (y
)));
4692 else if (SCM_FRACTIONP (y
))
4694 double xx
= SCM_REAL_VALUE (x
);
4698 return scm_from_bool (xx
< 0.0);
4699 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4703 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4705 else if (SCM_COMPLEXP (x
))
4707 if (SCM_I_INUMP (y
))
4708 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
4709 && (SCM_COMPLEX_IMAG (x
) == 0.0));
4710 else if (SCM_BIGP (y
))
4713 if (0.0 != SCM_COMPLEX_IMAG (x
))
4715 if (isnan (SCM_COMPLEX_REAL (x
)))
4717 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
4718 scm_remember_upto_here_1 (y
);
4719 return scm_from_bool (0 == cmp
);
4721 else if (SCM_REALP (y
))
4722 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
4723 && (SCM_COMPLEX_IMAG (x
) == 0.0));
4724 else if (SCM_COMPLEXP (y
))
4725 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
4726 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
4727 else if (SCM_FRACTIONP (y
))
4730 if (SCM_COMPLEX_IMAG (x
) != 0.0)
4732 xx
= SCM_COMPLEX_REAL (x
);
4736 return scm_from_bool (xx
< 0.0);
4737 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4741 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4743 else if (SCM_FRACTIONP (x
))
4745 if (SCM_I_INUMP (y
))
4747 else if (SCM_BIGP (y
))
4749 else if (SCM_REALP (y
))
4751 double yy
= SCM_REAL_VALUE (y
);
4755 return scm_from_bool (0.0 < yy
);
4756 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4759 else if (SCM_COMPLEXP (y
))
4762 if (SCM_COMPLEX_IMAG (y
) != 0.0)
4764 yy
= SCM_COMPLEX_REAL (y
);
4768 return scm_from_bool (0.0 < yy
);
4769 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4772 else if (SCM_FRACTIONP (y
))
4773 return scm_i_fraction_equalp (x
, y
);
4775 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4778 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
4782 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
4783 done are good for inums, but for bignums an answer can almost always be
4784 had by just examining a few high bits of the operands, as done by GMP in
4785 mpq_cmp. flonum/frac compares likewise, but with the slight complication
4786 of the float exponent to take into account. */
4788 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
4789 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
4790 (SCM x
, SCM y
, SCM rest
),
4791 "Return @code{#t} if the list of parameters is monotonically\n"
4793 #define FUNC_NAME s_scm_i_num_less_p
4795 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4797 while (!scm_is_null (rest
))
4799 if (scm_is_false (scm_less_p (x
, y
)))
4803 rest
= scm_cdr (rest
);
4805 return scm_less_p (x
, y
);
4809 scm_less_p (SCM x
, SCM y
)
4812 if (SCM_I_INUMP (x
))
4814 scm_t_inum xx
= SCM_I_INUM (x
);
4815 if (SCM_I_INUMP (y
))
4817 scm_t_inum yy
= SCM_I_INUM (y
);
4818 return scm_from_bool (xx
< yy
);
4820 else if (SCM_BIGP (y
))
4822 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4823 scm_remember_upto_here_1 (y
);
4824 return scm_from_bool (sgn
> 0);
4826 else if (SCM_REALP (y
))
4827 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
4828 else if (SCM_FRACTIONP (y
))
4830 /* "x < a/b" becomes "x*b < a" */
4832 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
4833 y
= SCM_FRACTION_NUMERATOR (y
);
4837 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4839 else if (SCM_BIGP (x
))
4841 if (SCM_I_INUMP (y
))
4843 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4844 scm_remember_upto_here_1 (x
);
4845 return scm_from_bool (sgn
< 0);
4847 else if (SCM_BIGP (y
))
4849 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4850 scm_remember_upto_here_2 (x
, y
);
4851 return scm_from_bool (cmp
< 0);
4853 else if (SCM_REALP (y
))
4856 if (isnan (SCM_REAL_VALUE (y
)))
4858 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
4859 scm_remember_upto_here_1 (x
);
4860 return scm_from_bool (cmp
< 0);
4862 else if (SCM_FRACTIONP (y
))
4865 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4867 else if (SCM_REALP (x
))
4869 if (SCM_I_INUMP (y
))
4870 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
4871 else if (SCM_BIGP (y
))
4874 if (isnan (SCM_REAL_VALUE (x
)))
4876 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
4877 scm_remember_upto_here_1 (y
);
4878 return scm_from_bool (cmp
> 0);
4880 else if (SCM_REALP (y
))
4881 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
4882 else if (SCM_FRACTIONP (y
))
4884 double xx
= SCM_REAL_VALUE (x
);
4888 return scm_from_bool (xx
< 0.0);
4889 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4893 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4895 else if (SCM_FRACTIONP (x
))
4897 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
4899 /* "a/b < y" becomes "a < y*b" */
4900 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
4901 x
= SCM_FRACTION_NUMERATOR (x
);
4904 else if (SCM_REALP (y
))
4906 double yy
= SCM_REAL_VALUE (y
);
4910 return scm_from_bool (0.0 < yy
);
4911 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4914 else if (SCM_FRACTIONP (y
))
4916 /* "a/b < c/d" becomes "a*d < c*b" */
4917 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
4918 SCM_FRACTION_DENOMINATOR (y
));
4919 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
4920 SCM_FRACTION_DENOMINATOR (x
));
4926 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4929 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
4933 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
4934 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
4935 (SCM x
, SCM y
, SCM rest
),
4936 "Return @code{#t} if the list of parameters is monotonically\n"
4938 #define FUNC_NAME s_scm_i_num_gr_p
4940 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4942 while (!scm_is_null (rest
))
4944 if (scm_is_false (scm_gr_p (x
, y
)))
4948 rest
= scm_cdr (rest
);
4950 return scm_gr_p (x
, y
);
4953 #define FUNC_NAME s_scm_i_num_gr_p
4955 scm_gr_p (SCM x
, SCM y
)
4957 if (!SCM_NUMBERP (x
))
4958 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
4959 else if (!SCM_NUMBERP (y
))
4960 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
4962 return scm_less_p (y
, x
);
4967 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
4968 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
4969 (SCM x
, SCM y
, SCM rest
),
4970 "Return @code{#t} if the list of parameters is monotonically\n"
4972 #define FUNC_NAME s_scm_i_num_leq_p
4974 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4976 while (!scm_is_null (rest
))
4978 if (scm_is_false (scm_leq_p (x
, y
)))
4982 rest
= scm_cdr (rest
);
4984 return scm_leq_p (x
, y
);
4987 #define FUNC_NAME s_scm_i_num_leq_p
4989 scm_leq_p (SCM x
, SCM y
)
4991 if (!SCM_NUMBERP (x
))
4992 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
4993 else if (!SCM_NUMBERP (y
))
4994 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
4995 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
4998 return scm_not (scm_less_p (y
, x
));
5003 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
5004 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
5005 (SCM x
, SCM y
, SCM rest
),
5006 "Return @code{#t} if the list of parameters is monotonically\n"
5008 #define FUNC_NAME s_scm_i_num_geq_p
5010 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
5012 while (!scm_is_null (rest
))
5014 if (scm_is_false (scm_geq_p (x
, y
)))
5018 rest
= scm_cdr (rest
);
5020 return scm_geq_p (x
, y
);
5023 #define FUNC_NAME s_scm_i_num_geq_p
5025 scm_geq_p (SCM x
, SCM y
)
5027 if (!SCM_NUMBERP (x
))
5028 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
5029 else if (!SCM_NUMBERP (y
))
5030 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
5031 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
5034 return scm_not (scm_less_p (x
, y
));
5039 SCM_PRIMITIVE_GENERIC (scm_zero_p
, "zero?", 1, 0, 0,
5041 "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
5043 #define FUNC_NAME s_scm_zero_p
5045 if (SCM_I_INUMP (z
))
5046 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
5047 else if (SCM_BIGP (z
))
5049 else if (SCM_REALP (z
))
5050 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
5051 else if (SCM_COMPLEXP (z
))
5052 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
5053 && SCM_COMPLEX_IMAG (z
) == 0.0);
5054 else if (SCM_FRACTIONP (z
))
5057 SCM_WTA_DISPATCH_1 (g_scm_zero_p
, z
, SCM_ARG1
, s_scm_zero_p
);
5062 SCM_PRIMITIVE_GENERIC (scm_positive_p
, "positive?", 1, 0, 0,
5064 "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
5066 #define FUNC_NAME s_scm_positive_p
5068 if (SCM_I_INUMP (x
))
5069 return scm_from_bool (SCM_I_INUM (x
) > 0);
5070 else if (SCM_BIGP (x
))
5072 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5073 scm_remember_upto_here_1 (x
);
5074 return scm_from_bool (sgn
> 0);
5076 else if (SCM_REALP (x
))
5077 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
5078 else if (SCM_FRACTIONP (x
))
5079 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
5081 SCM_WTA_DISPATCH_1 (g_scm_positive_p
, x
, SCM_ARG1
, s_scm_positive_p
);
5086 SCM_PRIMITIVE_GENERIC (scm_negative_p
, "negative?", 1, 0, 0,
5088 "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
5090 #define FUNC_NAME s_scm_negative_p
5092 if (SCM_I_INUMP (x
))
5093 return scm_from_bool (SCM_I_INUM (x
) < 0);
5094 else if (SCM_BIGP (x
))
5096 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5097 scm_remember_upto_here_1 (x
);
5098 return scm_from_bool (sgn
< 0);
5100 else if (SCM_REALP (x
))
5101 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
5102 else if (SCM_FRACTIONP (x
))
5103 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
5105 SCM_WTA_DISPATCH_1 (g_scm_negative_p
, x
, SCM_ARG1
, s_scm_negative_p
);
5110 /* scm_min and scm_max return an inexact when either argument is inexact, as
5111 required by r5rs. On that basis, for exact/inexact combinations the
5112 exact is converted to inexact to compare and possibly return. This is
5113 unlike scm_less_p above which takes some trouble to preserve all bits in
5114 its test, such trouble is not required for min and max. */
5116 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
5117 (SCM x
, SCM y
, SCM rest
),
5118 "Return the maximum of all parameter values.")
5119 #define FUNC_NAME s_scm_i_max
5121 while (!scm_is_null (rest
))
5122 { x
= scm_max (x
, y
);
5124 rest
= scm_cdr (rest
);
5126 return scm_max (x
, y
);
5130 #define s_max s_scm_i_max
5131 #define g_max g_scm_i_max
5134 scm_max (SCM x
, SCM y
)
5139 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
5140 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
5143 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
5146 if (SCM_I_INUMP (x
))
5148 scm_t_inum xx
= SCM_I_INUM (x
);
5149 if (SCM_I_INUMP (y
))
5151 scm_t_inum yy
= SCM_I_INUM (y
);
5152 return (xx
< yy
) ? y
: x
;
5154 else if (SCM_BIGP (y
))
5156 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5157 scm_remember_upto_here_1 (y
);
5158 return (sgn
< 0) ? x
: y
;
5160 else if (SCM_REALP (y
))
5163 double yyd
= SCM_REAL_VALUE (y
);
5166 return scm_from_double (xxd
);
5167 /* If y is a NaN, then "==" is false and we return the NaN */
5168 else if (SCM_LIKELY (!(xxd
== yyd
)))
5170 /* Handle signed zeroes properly */
5176 else if (SCM_FRACTIONP (y
))
5179 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
5182 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5184 else if (SCM_BIGP (x
))
5186 if (SCM_I_INUMP (y
))
5188 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5189 scm_remember_upto_here_1 (x
);
5190 return (sgn
< 0) ? y
: x
;
5192 else if (SCM_BIGP (y
))
5194 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
5195 scm_remember_upto_here_2 (x
, y
);
5196 return (cmp
> 0) ? x
: y
;
5198 else if (SCM_REALP (y
))
5200 /* if y==NaN then xx>yy is false, so we return the NaN y */
5203 xx
= scm_i_big2dbl (x
);
5204 yy
= SCM_REAL_VALUE (y
);
5205 return (xx
> yy
? scm_from_double (xx
) : y
);
5207 else if (SCM_FRACTIONP (y
))
5212 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5214 else if (SCM_REALP (x
))
5216 if (SCM_I_INUMP (y
))
5218 scm_t_inum yy
= SCM_I_INUM (y
);
5219 double xxd
= SCM_REAL_VALUE (x
);
5223 return scm_from_double (yyd
);
5224 /* If x is a NaN, then "==" is false and we return the NaN */
5225 else if (SCM_LIKELY (!(xxd
== yyd
)))
5227 /* Handle signed zeroes properly */
5233 else if (SCM_BIGP (y
))
5238 else if (SCM_REALP (y
))
5240 double xx
= SCM_REAL_VALUE (x
);
5241 double yy
= SCM_REAL_VALUE (y
);
5243 /* For purposes of max: +inf.0 > nan > everything else, per R6RS */
5246 else if (SCM_LIKELY (xx
< yy
))
5248 /* If neither (xx > yy) nor (xx < yy), then
5249 either they're equal or one is a NaN */
5250 else if (SCM_UNLIKELY (isnan (xx
)))
5251 return (isinf (yy
) == 1) ? y
: x
;
5252 else if (SCM_UNLIKELY (isnan (yy
)))
5253 return (isinf (xx
) == 1) ? x
: y
;
5254 /* xx == yy, but handle signed zeroes properly */
5255 else if (double_is_non_negative_zero (yy
))
5260 else if (SCM_FRACTIONP (y
))
5262 double yy
= scm_i_fraction2double (y
);
5263 double xx
= SCM_REAL_VALUE (x
);
5264 return (xx
< yy
) ? scm_from_double (yy
) : x
;
5267 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5269 else if (SCM_FRACTIONP (x
))
5271 if (SCM_I_INUMP (y
))
5275 else if (SCM_BIGP (y
))
5279 else if (SCM_REALP (y
))
5281 double xx
= scm_i_fraction2double (x
);
5282 /* if y==NaN then ">" is false, so we return the NaN y */
5283 return (xx
> SCM_REAL_VALUE (y
)) ? scm_from_double (xx
) : y
;
5285 else if (SCM_FRACTIONP (y
))
5290 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5293 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
5297 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
5298 (SCM x
, SCM y
, SCM rest
),
5299 "Return the minimum of all parameter values.")
5300 #define FUNC_NAME s_scm_i_min
5302 while (!scm_is_null (rest
))
5303 { x
= scm_min (x
, y
);
5305 rest
= scm_cdr (rest
);
5307 return scm_min (x
, y
);
5311 #define s_min s_scm_i_min
5312 #define g_min g_scm_i_min
5315 scm_min (SCM x
, SCM y
)
5320 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
5321 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
5324 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
5327 if (SCM_I_INUMP (x
))
5329 scm_t_inum xx
= SCM_I_INUM (x
);
5330 if (SCM_I_INUMP (y
))
5332 scm_t_inum yy
= SCM_I_INUM (y
);
5333 return (xx
< yy
) ? x
: y
;
5335 else if (SCM_BIGP (y
))
5337 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5338 scm_remember_upto_here_1 (y
);
5339 return (sgn
< 0) ? y
: x
;
5341 else if (SCM_REALP (y
))
5344 /* if y==NaN then "<" is false and we return NaN */
5345 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
5347 else if (SCM_FRACTIONP (y
))
5350 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
5353 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5355 else if (SCM_BIGP (x
))
5357 if (SCM_I_INUMP (y
))
5359 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5360 scm_remember_upto_here_1 (x
);
5361 return (sgn
< 0) ? x
: y
;
5363 else if (SCM_BIGP (y
))
5365 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
5366 scm_remember_upto_here_2 (x
, y
);
5367 return (cmp
> 0) ? y
: x
;
5369 else if (SCM_REALP (y
))
5371 /* if y==NaN then xx<yy is false, so we return the NaN y */
5374 xx
= scm_i_big2dbl (x
);
5375 yy
= SCM_REAL_VALUE (y
);
5376 return (xx
< yy
? scm_from_double (xx
) : y
);
5378 else if (SCM_FRACTIONP (y
))
5383 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5385 else if (SCM_REALP (x
))
5387 if (SCM_I_INUMP (y
))
5389 double z
= SCM_I_INUM (y
);
5390 /* if x==NaN then "<" is false and we return NaN */
5391 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
5393 else if (SCM_BIGP (y
))
5398 else if (SCM_REALP (y
))
5400 double xx
= SCM_REAL_VALUE (x
);
5401 double yy
= SCM_REAL_VALUE (y
);
5403 /* For purposes of min: -inf.0 < nan < everything else, per R6RS */
5406 else if (SCM_LIKELY (xx
> yy
))
5408 /* If neither (xx < yy) nor (xx > yy), then
5409 either they're equal or one is a NaN */
5410 else if (SCM_UNLIKELY (isnan (xx
)))
5411 return (isinf (yy
) == -1) ? y
: x
;
5412 else if (SCM_UNLIKELY (isnan (yy
)))
5413 return (isinf (xx
) == -1) ? x
: y
;
5414 /* xx == yy, but handle signed zeroes properly */
5415 else if (double_is_non_negative_zero (xx
))
5420 else if (SCM_FRACTIONP (y
))
5422 double yy
= scm_i_fraction2double (y
);
5423 double xx
= SCM_REAL_VALUE (x
);
5424 return (yy
< xx
) ? scm_from_double (yy
) : x
;
5427 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5429 else if (SCM_FRACTIONP (x
))
5431 if (SCM_I_INUMP (y
))
5435 else if (SCM_BIGP (y
))
5439 else if (SCM_REALP (y
))
5441 double xx
= scm_i_fraction2double (x
);
5442 /* if y==NaN then "<" is false, so we return the NaN y */
5443 return (xx
< SCM_REAL_VALUE (y
)) ? scm_from_double (xx
) : y
;
5445 else if (SCM_FRACTIONP (y
))
5450 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5453 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
5457 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
5458 (SCM x
, SCM y
, SCM rest
),
5459 "Return the sum of all parameter values. Return 0 if called without\n"
5461 #define FUNC_NAME s_scm_i_sum
5463 while (!scm_is_null (rest
))
5464 { x
= scm_sum (x
, y
);
5466 rest
= scm_cdr (rest
);
5468 return scm_sum (x
, y
);
5472 #define s_sum s_scm_i_sum
5473 #define g_sum g_scm_i_sum
5476 scm_sum (SCM x
, SCM y
)
5478 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5480 if (SCM_NUMBERP (x
)) return x
;
5481 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
5482 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
5485 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5487 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5489 scm_t_inum xx
= SCM_I_INUM (x
);
5490 scm_t_inum yy
= SCM_I_INUM (y
);
5491 scm_t_inum z
= xx
+ yy
;
5492 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_inum2big (z
);
5494 else if (SCM_BIGP (y
))
5499 else if (SCM_REALP (y
))
5501 scm_t_inum xx
= SCM_I_INUM (x
);
5502 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
5504 else if (SCM_COMPLEXP (y
))
5506 scm_t_inum xx
= SCM_I_INUM (x
);
5507 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
5508 SCM_COMPLEX_IMAG (y
));
5510 else if (SCM_FRACTIONP (y
))
5511 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
5512 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
5513 SCM_FRACTION_DENOMINATOR (y
));
5515 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5516 } else if (SCM_BIGP (x
))
5518 if (SCM_I_INUMP (y
))
5523 inum
= SCM_I_INUM (y
);
5526 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5529 SCM result
= scm_i_mkbig ();
5530 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
5531 scm_remember_upto_here_1 (x
);
5532 /* we know the result will have to be a bignum */
5535 return scm_i_normbig (result
);
5539 SCM result
= scm_i_mkbig ();
5540 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
5541 scm_remember_upto_here_1 (x
);
5542 /* we know the result will have to be a bignum */
5545 return scm_i_normbig (result
);
5548 else if (SCM_BIGP (y
))
5550 SCM result
= scm_i_mkbig ();
5551 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5552 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5553 mpz_add (SCM_I_BIG_MPZ (result
),
5556 scm_remember_upto_here_2 (x
, y
);
5557 /* we know the result will have to be a bignum */
5560 return scm_i_normbig (result
);
5562 else if (SCM_REALP (y
))
5564 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
5565 scm_remember_upto_here_1 (x
);
5566 return scm_from_double (result
);
5568 else if (SCM_COMPLEXP (y
))
5570 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
5571 + SCM_COMPLEX_REAL (y
));
5572 scm_remember_upto_here_1 (x
);
5573 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
5575 else if (SCM_FRACTIONP (y
))
5576 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
5577 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
5578 SCM_FRACTION_DENOMINATOR (y
));
5580 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5582 else if (SCM_REALP (x
))
5584 if (SCM_I_INUMP (y
))
5585 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
5586 else if (SCM_BIGP (y
))
5588 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
5589 scm_remember_upto_here_1 (y
);
5590 return scm_from_double (result
);
5592 else if (SCM_REALP (y
))
5593 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
5594 else if (SCM_COMPLEXP (y
))
5595 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
5596 SCM_COMPLEX_IMAG (y
));
5597 else if (SCM_FRACTIONP (y
))
5598 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
5600 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5602 else if (SCM_COMPLEXP (x
))
5604 if (SCM_I_INUMP (y
))
5605 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
5606 SCM_COMPLEX_IMAG (x
));
5607 else if (SCM_BIGP (y
))
5609 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
5610 + SCM_COMPLEX_REAL (x
));
5611 scm_remember_upto_here_1 (y
);
5612 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
5614 else if (SCM_REALP (y
))
5615 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
5616 SCM_COMPLEX_IMAG (x
));
5617 else if (SCM_COMPLEXP (y
))
5618 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
5619 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
5620 else if (SCM_FRACTIONP (y
))
5621 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
5622 SCM_COMPLEX_IMAG (x
));
5624 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5626 else if (SCM_FRACTIONP (x
))
5628 if (SCM_I_INUMP (y
))
5629 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
5630 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
5631 SCM_FRACTION_DENOMINATOR (x
));
5632 else if (SCM_BIGP (y
))
5633 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
5634 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
5635 SCM_FRACTION_DENOMINATOR (x
));
5636 else if (SCM_REALP (y
))
5637 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
5638 else if (SCM_COMPLEXP (y
))
5639 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
5640 SCM_COMPLEX_IMAG (y
));
5641 else if (SCM_FRACTIONP (y
))
5642 /* a/b + c/d = (ad + bc) / bd */
5643 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5644 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
5645 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
5647 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5650 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
5654 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
5656 "Return @math{@var{x}+1}.")
5657 #define FUNC_NAME s_scm_oneplus
5659 return scm_sum (x
, SCM_INUM1
);
5664 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
5665 (SCM x
, SCM y
, SCM rest
),
5666 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
5667 "the sum of all but the first argument are subtracted from the first\n"
5669 #define FUNC_NAME s_scm_i_difference
5671 while (!scm_is_null (rest
))
5672 { x
= scm_difference (x
, y
);
5674 rest
= scm_cdr (rest
);
5676 return scm_difference (x
, y
);
5680 #define s_difference s_scm_i_difference
5681 #define g_difference g_scm_i_difference
5684 scm_difference (SCM x
, SCM y
)
5685 #define FUNC_NAME s_difference
5687 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5690 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
5692 if (SCM_I_INUMP (x
))
5694 scm_t_inum xx
= -SCM_I_INUM (x
);
5695 if (SCM_FIXABLE (xx
))
5696 return SCM_I_MAKINUM (xx
);
5698 return scm_i_inum2big (xx
);
5700 else if (SCM_BIGP (x
))
5701 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
5702 bignum, but negating that gives a fixnum. */
5703 return scm_i_normbig (scm_i_clonebig (x
, 0));
5704 else if (SCM_REALP (x
))
5705 return scm_from_double (-SCM_REAL_VALUE (x
));
5706 else if (SCM_COMPLEXP (x
))
5707 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
5708 -SCM_COMPLEX_IMAG (x
));
5709 else if (SCM_FRACTIONP (x
))
5710 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
5711 SCM_FRACTION_DENOMINATOR (x
));
5713 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
5716 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5718 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5720 scm_t_inum xx
= SCM_I_INUM (x
);
5721 scm_t_inum yy
= SCM_I_INUM (y
);
5722 scm_t_inum z
= xx
- yy
;
5723 if (SCM_FIXABLE (z
))
5724 return SCM_I_MAKINUM (z
);
5726 return scm_i_inum2big (z
);
5728 else if (SCM_BIGP (y
))
5730 /* inum-x - big-y */
5731 scm_t_inum xx
= SCM_I_INUM (x
);
5735 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
5736 bignum, but negating that gives a fixnum. */
5737 return scm_i_normbig (scm_i_clonebig (y
, 0));
5741 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5742 SCM result
= scm_i_mkbig ();
5745 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
5748 /* x - y == -(y + -x) */
5749 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
5750 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5752 scm_remember_upto_here_1 (y
);
5754 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
5755 /* we know the result will have to be a bignum */
5758 return scm_i_normbig (result
);
5761 else if (SCM_REALP (y
))
5763 scm_t_inum xx
= SCM_I_INUM (x
);
5766 * We need to handle x == exact 0
5767 * specially because R6RS states that:
5768 * (- 0.0) ==> -0.0 and
5769 * (- 0.0 0.0) ==> 0.0
5770 * and the scheme compiler changes
5771 * (- 0.0) into (- 0 0.0)
5772 * So we need to treat (- 0 0.0) like (- 0.0).
5773 * At the C level, (-x) is different than (0.0 - x).
5774 * (0.0 - 0.0) ==> 0.0, but (- 0.0) ==> -0.0.
5777 return scm_from_double (- SCM_REAL_VALUE (y
));
5779 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
5781 else if (SCM_COMPLEXP (y
))
5783 scm_t_inum xx
= SCM_I_INUM (x
);
5785 /* We need to handle x == exact 0 specially.
5786 See the comment above (for SCM_REALP (y)) */
5788 return scm_c_make_rectangular (- SCM_COMPLEX_REAL (y
),
5789 - SCM_COMPLEX_IMAG (y
));
5791 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
5792 - SCM_COMPLEX_IMAG (y
));
5794 else if (SCM_FRACTIONP (y
))
5795 /* a - b/c = (ac - b) / c */
5796 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5797 SCM_FRACTION_NUMERATOR (y
)),
5798 SCM_FRACTION_DENOMINATOR (y
));
5800 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5802 else if (SCM_BIGP (x
))
5804 if (SCM_I_INUMP (y
))
5806 /* big-x - inum-y */
5807 scm_t_inum yy
= SCM_I_INUM (y
);
5808 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5810 scm_remember_upto_here_1 (x
);
5812 return (SCM_FIXABLE (-yy
) ?
5813 SCM_I_MAKINUM (-yy
) : scm_from_inum (-yy
));
5816 SCM result
= scm_i_mkbig ();
5819 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
5821 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
5822 scm_remember_upto_here_1 (x
);
5824 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
5825 /* we know the result will have to be a bignum */
5828 return scm_i_normbig (result
);
5831 else if (SCM_BIGP (y
))
5833 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5834 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5835 SCM result
= scm_i_mkbig ();
5836 mpz_sub (SCM_I_BIG_MPZ (result
),
5839 scm_remember_upto_here_2 (x
, y
);
5840 /* we know the result will have to be a bignum */
5841 if ((sgn_x
== 1) && (sgn_y
== -1))
5843 if ((sgn_x
== -1) && (sgn_y
== 1))
5845 return scm_i_normbig (result
);
5847 else if (SCM_REALP (y
))
5849 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
5850 scm_remember_upto_here_1 (x
);
5851 return scm_from_double (result
);
5853 else if (SCM_COMPLEXP (y
))
5855 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
5856 - SCM_COMPLEX_REAL (y
));
5857 scm_remember_upto_here_1 (x
);
5858 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
5860 else if (SCM_FRACTIONP (y
))
5861 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5862 SCM_FRACTION_NUMERATOR (y
)),
5863 SCM_FRACTION_DENOMINATOR (y
));
5864 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5866 else if (SCM_REALP (x
))
5868 if (SCM_I_INUMP (y
))
5869 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
5870 else if (SCM_BIGP (y
))
5872 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
5873 scm_remember_upto_here_1 (x
);
5874 return scm_from_double (result
);
5876 else if (SCM_REALP (y
))
5877 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
5878 else if (SCM_COMPLEXP (y
))
5879 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
5880 -SCM_COMPLEX_IMAG (y
));
5881 else if (SCM_FRACTIONP (y
))
5882 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
5884 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5886 else if (SCM_COMPLEXP (x
))
5888 if (SCM_I_INUMP (y
))
5889 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
5890 SCM_COMPLEX_IMAG (x
));
5891 else if (SCM_BIGP (y
))
5893 double real_part
= (SCM_COMPLEX_REAL (x
)
5894 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
5895 scm_remember_upto_here_1 (x
);
5896 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
5898 else if (SCM_REALP (y
))
5899 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
5900 SCM_COMPLEX_IMAG (x
));
5901 else if (SCM_COMPLEXP (y
))
5902 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
5903 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
5904 else if (SCM_FRACTIONP (y
))
5905 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
5906 SCM_COMPLEX_IMAG (x
));
5908 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5910 else if (SCM_FRACTIONP (x
))
5912 if (SCM_I_INUMP (y
))
5913 /* a/b - c = (a - cb) / b */
5914 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
5915 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
5916 SCM_FRACTION_DENOMINATOR (x
));
5917 else if (SCM_BIGP (y
))
5918 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
5919 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
5920 SCM_FRACTION_DENOMINATOR (x
));
5921 else if (SCM_REALP (y
))
5922 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
5923 else if (SCM_COMPLEXP (y
))
5924 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
5925 -SCM_COMPLEX_IMAG (y
));
5926 else if (SCM_FRACTIONP (y
))
5927 /* a/b - c/d = (ad - bc) / bd */
5928 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5929 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
5930 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
5932 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5935 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
5940 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
5942 "Return @math{@var{x}-1}.")
5943 #define FUNC_NAME s_scm_oneminus
5945 return scm_difference (x
, SCM_INUM1
);
5950 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
5951 (SCM x
, SCM y
, SCM rest
),
5952 "Return the product of all arguments. If called without arguments,\n"
5954 #define FUNC_NAME s_scm_i_product
5956 while (!scm_is_null (rest
))
5957 { x
= scm_product (x
, y
);
5959 rest
= scm_cdr (rest
);
5961 return scm_product (x
, y
);
5965 #define s_product s_scm_i_product
5966 #define g_product g_scm_i_product
5969 scm_product (SCM x
, SCM y
)
5971 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5974 return SCM_I_MAKINUM (1L);
5975 else if (SCM_NUMBERP (x
))
5978 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
5981 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5986 xx
= SCM_I_INUM (x
);
5991 /* exact1 is the universal multiplicative identity */
5995 /* exact0 times a fixnum is exact0: optimize this case */
5996 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5998 /* if the other argument is inexact, the result is inexact,
5999 and we must do the multiplication in order to handle
6000 infinities and NaNs properly. */
6001 else if (SCM_REALP (y
))
6002 return scm_from_double (0.0 * SCM_REAL_VALUE (y
));
6003 else if (SCM_COMPLEXP (y
))
6004 return scm_c_make_rectangular (0.0 * SCM_COMPLEX_REAL (y
),
6005 0.0 * SCM_COMPLEX_IMAG (y
));
6006 /* we've already handled inexact numbers,
6007 so y must be exact, and we return exact0 */
6008 else if (SCM_NUMP (y
))
6011 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6015 * This case is important for more than just optimization.
6016 * It handles the case of negating
6017 * (+ 1 most-positive-fixnum) aka (- most-negative-fixnum),
6018 * which is a bignum that must be changed back into a fixnum.
6019 * Failure to do so will cause the following to return #f:
6020 * (= most-negative-fixnum (* -1 (- most-negative-fixnum)))
6022 return scm_difference(y
, SCM_UNDEFINED
);
6026 if (SCM_LIKELY (SCM_I_INUMP (y
)))
6028 scm_t_inum yy
= SCM_I_INUM (y
);
6029 scm_t_inum kk
= xx
* yy
;
6030 SCM k
= SCM_I_MAKINUM (kk
);
6031 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
6035 SCM result
= scm_i_inum2big (xx
);
6036 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
6037 return scm_i_normbig (result
);
6040 else if (SCM_BIGP (y
))
6042 SCM result
= scm_i_mkbig ();
6043 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
6044 scm_remember_upto_here_1 (y
);
6047 else if (SCM_REALP (y
))
6048 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
6049 else if (SCM_COMPLEXP (y
))
6050 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
6051 xx
* SCM_COMPLEX_IMAG (y
));
6052 else if (SCM_FRACTIONP (y
))
6053 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
6054 SCM_FRACTION_DENOMINATOR (y
));
6056 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6058 else if (SCM_BIGP (x
))
6060 if (SCM_I_INUMP (y
))
6065 else if (SCM_BIGP (y
))
6067 SCM result
= scm_i_mkbig ();
6068 mpz_mul (SCM_I_BIG_MPZ (result
),
6071 scm_remember_upto_here_2 (x
, y
);
6074 else if (SCM_REALP (y
))
6076 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
6077 scm_remember_upto_here_1 (x
);
6078 return scm_from_double (result
);
6080 else if (SCM_COMPLEXP (y
))
6082 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
6083 scm_remember_upto_here_1 (x
);
6084 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
6085 z
* SCM_COMPLEX_IMAG (y
));
6087 else if (SCM_FRACTIONP (y
))
6088 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
6089 SCM_FRACTION_DENOMINATOR (y
));
6091 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6093 else if (SCM_REALP (x
))
6095 if (SCM_I_INUMP (y
))
6100 else if (SCM_BIGP (y
))
6102 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
6103 scm_remember_upto_here_1 (y
);
6104 return scm_from_double (result
);
6106 else if (SCM_REALP (y
))
6107 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
6108 else if (SCM_COMPLEXP (y
))
6109 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
6110 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
6111 else if (SCM_FRACTIONP (y
))
6112 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
6114 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6116 else if (SCM_COMPLEXP (x
))
6118 if (SCM_I_INUMP (y
))
6123 else if (SCM_BIGP (y
))
6125 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6126 scm_remember_upto_here_1 (y
);
6127 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
6128 z
* SCM_COMPLEX_IMAG (x
));
6130 else if (SCM_REALP (y
))
6131 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
6132 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
6133 else if (SCM_COMPLEXP (y
))
6135 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
6136 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
6137 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
6138 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
6140 else if (SCM_FRACTIONP (y
))
6142 double yy
= scm_i_fraction2double (y
);
6143 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
6144 yy
* SCM_COMPLEX_IMAG (x
));
6147 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6149 else if (SCM_FRACTIONP (x
))
6151 if (SCM_I_INUMP (y
))
6152 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
6153 SCM_FRACTION_DENOMINATOR (x
));
6154 else if (SCM_BIGP (y
))
6155 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
6156 SCM_FRACTION_DENOMINATOR (x
));
6157 else if (SCM_REALP (y
))
6158 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
6159 else if (SCM_COMPLEXP (y
))
6161 double xx
= scm_i_fraction2double (x
);
6162 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
6163 xx
* SCM_COMPLEX_IMAG (y
));
6165 else if (SCM_FRACTIONP (y
))
6166 /* a/b * c/d = ac / bd */
6167 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
6168 SCM_FRACTION_NUMERATOR (y
)),
6169 scm_product (SCM_FRACTION_DENOMINATOR (x
),
6170 SCM_FRACTION_DENOMINATOR (y
)));
6172 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6175 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
6178 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
6179 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
6180 #define ALLOW_DIVIDE_BY_ZERO
6181 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
6184 /* The code below for complex division is adapted from the GNU
6185 libstdc++, which adapted it from f2c's libF77, and is subject to
6188 /****************************************************************
6189 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
6191 Permission to use, copy, modify, and distribute this software
6192 and its documentation for any purpose and without fee is hereby
6193 granted, provided that the above copyright notice appear in all
6194 copies and that both that the copyright notice and this
6195 permission notice and warranty disclaimer appear in supporting
6196 documentation, and that the names of AT&T Bell Laboratories or
6197 Bellcore or any of their entities not be used in advertising or
6198 publicity pertaining to distribution of the software without
6199 specific, written prior permission.
6201 AT&T and Bellcore disclaim all warranties with regard to this
6202 software, including all implied warranties of merchantability
6203 and fitness. In no event shall AT&T or Bellcore be liable for
6204 any special, indirect or consequential damages or any damages
6205 whatsoever resulting from loss of use, data or profits, whether
6206 in an action of contract, negligence or other tortious action,
6207 arising out of or in connection with the use or performance of
6209 ****************************************************************/
6211 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
6212 (SCM x
, SCM y
, SCM rest
),
6213 "Divide the first argument by the product of the remaining\n"
6214 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
6216 #define FUNC_NAME s_scm_i_divide
6218 while (!scm_is_null (rest
))
6219 { x
= scm_divide (x
, y
);
6221 rest
= scm_cdr (rest
);
6223 return scm_divide (x
, y
);
6227 #define s_divide s_scm_i_divide
6228 #define g_divide g_scm_i_divide
6231 do_divide (SCM x
, SCM y
, int inexact
)
6232 #define FUNC_NAME s_divide
6236 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
6239 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
6240 else if (SCM_I_INUMP (x
))
6242 scm_t_inum xx
= SCM_I_INUM (x
);
6243 if (xx
== 1 || xx
== -1)
6245 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6247 scm_num_overflow (s_divide
);
6252 return scm_from_double (1.0 / (double) xx
);
6253 else return scm_i_make_ratio (SCM_INUM1
, x
);
6256 else if (SCM_BIGP (x
))
6259 return scm_from_double (1.0 / scm_i_big2dbl (x
));
6260 else return scm_i_make_ratio (SCM_INUM1
, x
);
6262 else if (SCM_REALP (x
))
6264 double xx
= SCM_REAL_VALUE (x
);
6265 #ifndef ALLOW_DIVIDE_BY_ZERO
6267 scm_num_overflow (s_divide
);
6270 return scm_from_double (1.0 / xx
);
6272 else if (SCM_COMPLEXP (x
))
6274 double r
= SCM_COMPLEX_REAL (x
);
6275 double i
= SCM_COMPLEX_IMAG (x
);
6276 if (fabs(r
) <= fabs(i
))
6279 double d
= i
* (1.0 + t
* t
);
6280 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
6285 double d
= r
* (1.0 + t
* t
);
6286 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
6289 else if (SCM_FRACTIONP (x
))
6290 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
6291 SCM_FRACTION_NUMERATOR (x
));
6293 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
6296 if (SCM_LIKELY (SCM_I_INUMP (x
)))
6298 scm_t_inum xx
= SCM_I_INUM (x
);
6299 if (SCM_LIKELY (SCM_I_INUMP (y
)))
6301 scm_t_inum yy
= SCM_I_INUM (y
);
6304 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6305 scm_num_overflow (s_divide
);
6307 return scm_from_double ((double) xx
/ (double) yy
);
6310 else if (xx
% yy
!= 0)
6313 return scm_from_double ((double) xx
/ (double) yy
);
6314 else return scm_i_make_ratio (x
, y
);
6318 scm_t_inum z
= xx
/ yy
;
6319 if (SCM_FIXABLE (z
))
6320 return SCM_I_MAKINUM (z
);
6322 return scm_i_inum2big (z
);
6325 else if (SCM_BIGP (y
))
6328 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
6329 else return scm_i_make_ratio (x
, y
);
6331 else if (SCM_REALP (y
))
6333 double yy
= SCM_REAL_VALUE (y
);
6334 #ifndef ALLOW_DIVIDE_BY_ZERO
6336 scm_num_overflow (s_divide
);
6339 return scm_from_double ((double) xx
/ yy
);
6341 else if (SCM_COMPLEXP (y
))
6344 complex_div
: /* y _must_ be a complex number */
6346 double r
= SCM_COMPLEX_REAL (y
);
6347 double i
= SCM_COMPLEX_IMAG (y
);
6348 if (fabs(r
) <= fabs(i
))
6351 double d
= i
* (1.0 + t
* t
);
6352 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
6357 double d
= r
* (1.0 + t
* t
);
6358 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
6362 else if (SCM_FRACTIONP (y
))
6363 /* a / b/c = ac / b */
6364 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
6365 SCM_FRACTION_NUMERATOR (y
));
6367 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6369 else if (SCM_BIGP (x
))
6371 if (SCM_I_INUMP (y
))
6373 scm_t_inum yy
= SCM_I_INUM (y
);
6376 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6377 scm_num_overflow (s_divide
);
6379 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
6380 scm_remember_upto_here_1 (x
);
6381 return (sgn
== 0) ? scm_nan () : scm_inf ();
6388 /* FIXME: HMM, what are the relative performance issues here?
6389 We need to test. Is it faster on average to test
6390 divisible_p, then perform whichever operation, or is it
6391 faster to perform the integer div opportunistically and
6392 switch to real if there's a remainder? For now we take the
6393 middle ground: test, then if divisible, use the faster div
6396 scm_t_inum abs_yy
= yy
< 0 ? -yy
: yy
;
6397 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
6401 SCM result
= scm_i_mkbig ();
6402 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
6403 scm_remember_upto_here_1 (x
);
6405 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
6406 return scm_i_normbig (result
);
6411 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
6412 else return scm_i_make_ratio (x
, y
);
6416 else if (SCM_BIGP (y
))
6421 /* It's easily possible for the ratio x/y to fit a double
6422 but one or both x and y be too big to fit a double,
6423 hence the use of mpq_get_d rather than converting and
6426 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
6427 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
6428 return scm_from_double (mpq_get_d (q
));
6432 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
6436 SCM result
= scm_i_mkbig ();
6437 mpz_divexact (SCM_I_BIG_MPZ (result
),
6440 scm_remember_upto_here_2 (x
, y
);
6441 return scm_i_normbig (result
);
6444 return scm_i_make_ratio (x
, y
);
6447 else if (SCM_REALP (y
))
6449 double yy
= SCM_REAL_VALUE (y
);
6450 #ifndef ALLOW_DIVIDE_BY_ZERO
6452 scm_num_overflow (s_divide
);
6455 return scm_from_double (scm_i_big2dbl (x
) / yy
);
6457 else if (SCM_COMPLEXP (y
))
6459 a
= scm_i_big2dbl (x
);
6462 else if (SCM_FRACTIONP (y
))
6463 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
6464 SCM_FRACTION_NUMERATOR (y
));
6466 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6468 else if (SCM_REALP (x
))
6470 double rx
= SCM_REAL_VALUE (x
);
6471 if (SCM_I_INUMP (y
))
6473 scm_t_inum yy
= SCM_I_INUM (y
);
6474 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6476 scm_num_overflow (s_divide
);
6479 return scm_from_double (rx
/ (double) yy
);
6481 else if (SCM_BIGP (y
))
6483 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6484 scm_remember_upto_here_1 (y
);
6485 return scm_from_double (rx
/ dby
);
6487 else if (SCM_REALP (y
))
6489 double yy
= SCM_REAL_VALUE (y
);
6490 #ifndef ALLOW_DIVIDE_BY_ZERO
6492 scm_num_overflow (s_divide
);
6495 return scm_from_double (rx
/ yy
);
6497 else if (SCM_COMPLEXP (y
))
6502 else if (SCM_FRACTIONP (y
))
6503 return scm_from_double (rx
/ scm_i_fraction2double (y
));
6505 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6507 else if (SCM_COMPLEXP (x
))
6509 double rx
= SCM_COMPLEX_REAL (x
);
6510 double ix
= SCM_COMPLEX_IMAG (x
);
6511 if (SCM_I_INUMP (y
))
6513 scm_t_inum yy
= SCM_I_INUM (y
);
6514 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6516 scm_num_overflow (s_divide
);
6521 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
6524 else if (SCM_BIGP (y
))
6526 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6527 scm_remember_upto_here_1 (y
);
6528 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
6530 else if (SCM_REALP (y
))
6532 double yy
= SCM_REAL_VALUE (y
);
6533 #ifndef ALLOW_DIVIDE_BY_ZERO
6535 scm_num_overflow (s_divide
);
6538 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
6540 else if (SCM_COMPLEXP (y
))
6542 double ry
= SCM_COMPLEX_REAL (y
);
6543 double iy
= SCM_COMPLEX_IMAG (y
);
6544 if (fabs(ry
) <= fabs(iy
))
6547 double d
= iy
* (1.0 + t
* t
);
6548 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
6553 double d
= ry
* (1.0 + t
* t
);
6554 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
6557 else if (SCM_FRACTIONP (y
))
6559 double yy
= scm_i_fraction2double (y
);
6560 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
6563 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6565 else if (SCM_FRACTIONP (x
))
6567 if (SCM_I_INUMP (y
))
6569 scm_t_inum yy
= SCM_I_INUM (y
);
6570 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6572 scm_num_overflow (s_divide
);
6575 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
6576 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
6578 else if (SCM_BIGP (y
))
6580 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
6581 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
6583 else if (SCM_REALP (y
))
6585 double yy
= SCM_REAL_VALUE (y
);
6586 #ifndef ALLOW_DIVIDE_BY_ZERO
6588 scm_num_overflow (s_divide
);
6591 return scm_from_double (scm_i_fraction2double (x
) / yy
);
6593 else if (SCM_COMPLEXP (y
))
6595 a
= scm_i_fraction2double (x
);
6598 else if (SCM_FRACTIONP (y
))
6599 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
6600 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
6602 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6605 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
6609 scm_divide (SCM x
, SCM y
)
6611 return do_divide (x
, y
, 0);
6614 static SCM
scm_divide2real (SCM x
, SCM y
)
6616 return do_divide (x
, y
, 1);
6622 scm_c_truncate (double x
)
6633 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
6634 half-way case (ie. when x is an integer plus 0.5) going upwards.
6635 Then half-way cases are identified and adjusted down if the
6636 round-upwards didn't give the desired even integer.
6638 "plus_half == result" identifies a half-way case. If plus_half, which is
6639 x + 0.5, is an integer then x must be an integer plus 0.5.
6641 An odd "result" value is identified with result/2 != floor(result/2).
6642 This is done with plus_half, since that value is ready for use sooner in
6643 a pipelined cpu, and we're already requiring plus_half == result.
6645 Note however that we need to be careful when x is big and already an
6646 integer. In that case "x+0.5" may round to an adjacent integer, causing
6647 us to return such a value, incorrectly. For instance if the hardware is
6648 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
6649 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
6650 returned. Or if the hardware is in round-upwards mode, then other bigger
6651 values like say x == 2^128 will see x+0.5 rounding up to the next higher
6652 representable value, 2^128+2^76 (or whatever), again incorrect.
6654 These bad roundings of x+0.5 are avoided by testing at the start whether
6655 x is already an integer. If it is then clearly that's the desired result
6656 already. And if it's not then the exponent must be small enough to allow
6657 an 0.5 to be represented, and hence added without a bad rounding. */
6660 scm_c_round (double x
)
6662 double plus_half
, result
;
6667 plus_half
= x
+ 0.5;
6668 result
= floor (plus_half
);
6669 /* Adjust so that the rounding is towards even. */
6670 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
6675 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
6677 "Round the number @var{x} towards zero.")
6678 #define FUNC_NAME s_scm_truncate_number
6680 if (scm_is_false (scm_negative_p (x
)))
6681 return scm_floor (x
);
6683 return scm_ceiling (x
);
6687 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
6689 "Round the number @var{x} towards the nearest integer. "
6690 "When it is exactly halfway between two integers, "
6691 "round towards the even one.")
6692 #define FUNC_NAME s_scm_round_number
6694 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6696 else if (SCM_REALP (x
))
6697 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
6700 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
6701 single quotient+remainder division then examining to see which way
6702 the rounding should go. */
6703 SCM plus_half
= scm_sum (x
, exactly_one_half
);
6704 SCM result
= scm_floor (plus_half
);
6705 /* Adjust so that the rounding is towards even. */
6706 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
6707 && scm_is_true (scm_odd_p (result
)))
6708 return scm_difference (result
, SCM_INUM1
);
6715 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
6717 "Round the number @var{x} towards minus infinity.")
6718 #define FUNC_NAME s_scm_floor
6720 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6722 else if (SCM_REALP (x
))
6723 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
6724 else if (SCM_FRACTIONP (x
))
6726 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
6727 SCM_FRACTION_DENOMINATOR (x
));
6728 if (scm_is_false (scm_negative_p (x
)))
6730 /* For positive x, rounding towards zero is correct. */
6735 /* For negative x, we need to return q-1 unless x is an
6736 integer. But fractions are never integer, per our
6738 return scm_difference (q
, SCM_INUM1
);
6742 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
6746 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
6748 "Round the number @var{x} towards infinity.")
6749 #define FUNC_NAME s_scm_ceiling
6751 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6753 else if (SCM_REALP (x
))
6754 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
6755 else if (SCM_FRACTIONP (x
))
6757 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
6758 SCM_FRACTION_DENOMINATOR (x
));
6759 if (scm_is_false (scm_positive_p (x
)))
6761 /* For negative x, rounding towards zero is correct. */
6766 /* For positive x, we need to return q+1 unless x is an
6767 integer. But fractions are never integer, per our
6769 return scm_sum (q
, SCM_INUM1
);
6773 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
6777 SCM_PRIMITIVE_GENERIC (scm_expt
, "expt", 2, 0, 0,
6779 "Return @var{x} raised to the power of @var{y}.")
6780 #define FUNC_NAME s_scm_expt
6782 if (scm_is_integer (y
))
6784 if (scm_is_true (scm_exact_p (y
)))
6785 return scm_integer_expt (x
, y
);
6788 /* Here we handle the case where the exponent is an inexact
6789 integer. We make the exponent exact in order to use
6790 scm_integer_expt, and thus avoid the spurious imaginary
6791 parts that may result from round-off errors in the general
6792 e^(y log x) method below (for example when squaring a large
6793 negative number). In this case, we must return an inexact
6794 result for correctness. We also make the base inexact so
6795 that scm_integer_expt will use fast inexact arithmetic
6796 internally. Note that making the base inexact is not
6797 sufficient to guarantee an inexact result, because
6798 scm_integer_expt will return an exact 1 when the exponent
6799 is 0, even if the base is inexact. */
6800 return scm_exact_to_inexact
6801 (scm_integer_expt (scm_exact_to_inexact (x
),
6802 scm_inexact_to_exact (y
)));
6805 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
6807 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
6809 else if (scm_is_complex (x
) && scm_is_complex (y
))
6810 return scm_exp (scm_product (scm_log (x
), y
));
6811 else if (scm_is_complex (x
))
6812 SCM_WTA_DISPATCH_2 (g_scm_expt
, x
, y
, SCM_ARG2
, s_scm_expt
);
6814 SCM_WTA_DISPATCH_2 (g_scm_expt
, x
, y
, SCM_ARG1
, s_scm_expt
);
6818 /* sin/cos/tan/asin/acos/atan
6819 sinh/cosh/tanh/asinh/acosh/atanh
6820 Derived from "Transcen.scm", Complex trancendental functions for SCM.
6821 Written by Jerry D. Hedden, (C) FSF.
6822 See the file `COPYING' for terms applying to this program. */
6824 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
6826 "Compute the sine of @var{z}.")
6827 #define FUNC_NAME s_scm_sin
6829 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6830 return z
; /* sin(exact0) = exact0 */
6831 else if (scm_is_real (z
))
6832 return scm_from_double (sin (scm_to_double (z
)));
6833 else if (SCM_COMPLEXP (z
))
6835 x
= SCM_COMPLEX_REAL (z
);
6836 y
= SCM_COMPLEX_IMAG (z
);
6837 return scm_c_make_rectangular (sin (x
) * cosh (y
),
6838 cos (x
) * sinh (y
));
6841 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
6845 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
6847 "Compute the cosine of @var{z}.")
6848 #define FUNC_NAME s_scm_cos
6850 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6851 return SCM_INUM1
; /* cos(exact0) = exact1 */
6852 else if (scm_is_real (z
))
6853 return scm_from_double (cos (scm_to_double (z
)));
6854 else if (SCM_COMPLEXP (z
))
6856 x
= SCM_COMPLEX_REAL (z
);
6857 y
= SCM_COMPLEX_IMAG (z
);
6858 return scm_c_make_rectangular (cos (x
) * cosh (y
),
6859 -sin (x
) * sinh (y
));
6862 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
6866 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
6868 "Compute the tangent of @var{z}.")
6869 #define FUNC_NAME s_scm_tan
6871 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6872 return z
; /* tan(exact0) = exact0 */
6873 else if (scm_is_real (z
))
6874 return scm_from_double (tan (scm_to_double (z
)));
6875 else if (SCM_COMPLEXP (z
))
6877 x
= 2.0 * SCM_COMPLEX_REAL (z
);
6878 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
6879 w
= cos (x
) + cosh (y
);
6880 #ifndef ALLOW_DIVIDE_BY_ZERO
6882 scm_num_overflow (s_scm_tan
);
6884 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
6887 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
6891 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
6893 "Compute the hyperbolic sine of @var{z}.")
6894 #define FUNC_NAME s_scm_sinh
6896 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6897 return z
; /* sinh(exact0) = exact0 */
6898 else if (scm_is_real (z
))
6899 return scm_from_double (sinh (scm_to_double (z
)));
6900 else if (SCM_COMPLEXP (z
))
6902 x
= SCM_COMPLEX_REAL (z
);
6903 y
= SCM_COMPLEX_IMAG (z
);
6904 return scm_c_make_rectangular (sinh (x
) * cos (y
),
6905 cosh (x
) * sin (y
));
6908 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
6912 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
6914 "Compute the hyperbolic cosine of @var{z}.")
6915 #define FUNC_NAME s_scm_cosh
6917 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6918 return SCM_INUM1
; /* cosh(exact0) = exact1 */
6919 else if (scm_is_real (z
))
6920 return scm_from_double (cosh (scm_to_double (z
)));
6921 else if (SCM_COMPLEXP (z
))
6923 x
= SCM_COMPLEX_REAL (z
);
6924 y
= SCM_COMPLEX_IMAG (z
);
6925 return scm_c_make_rectangular (cosh (x
) * cos (y
),
6926 sinh (x
) * sin (y
));
6929 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
6933 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
6935 "Compute the hyperbolic tangent of @var{z}.")
6936 #define FUNC_NAME s_scm_tanh
6938 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6939 return z
; /* tanh(exact0) = exact0 */
6940 else if (scm_is_real (z
))
6941 return scm_from_double (tanh (scm_to_double (z
)));
6942 else if (SCM_COMPLEXP (z
))
6944 x
= 2.0 * SCM_COMPLEX_REAL (z
);
6945 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
6946 w
= cosh (x
) + cos (y
);
6947 #ifndef ALLOW_DIVIDE_BY_ZERO
6949 scm_num_overflow (s_scm_tanh
);
6951 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
6954 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
6958 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
6960 "Compute the arc sine of @var{z}.")
6961 #define FUNC_NAME s_scm_asin
6963 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
6964 return z
; /* asin(exact0) = exact0 */
6965 else if (scm_is_real (z
))
6967 double w
= scm_to_double (z
);
6968 if (w
>= -1.0 && w
<= 1.0)
6969 return scm_from_double (asin (w
));
6971 return scm_product (scm_c_make_rectangular (0, -1),
6972 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
6974 else if (SCM_COMPLEXP (z
))
6976 x
= SCM_COMPLEX_REAL (z
);
6977 y
= SCM_COMPLEX_IMAG (z
);
6978 return scm_product (scm_c_make_rectangular (0, -1),
6979 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
6982 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
6986 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
6988 "Compute the arc cosine of @var{z}.")
6989 #define FUNC_NAME s_scm_acos
6991 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM1
)))
6992 return SCM_INUM0
; /* acos(exact1) = exact0 */
6993 else if (scm_is_real (z
))
6995 double w
= scm_to_double (z
);
6996 if (w
>= -1.0 && w
<= 1.0)
6997 return scm_from_double (acos (w
));
6999 return scm_sum (scm_from_double (acos (0.0)),
7000 scm_product (scm_c_make_rectangular (0, 1),
7001 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
7003 else if (SCM_COMPLEXP (z
))
7005 x
= SCM_COMPLEX_REAL (z
);
7006 y
= SCM_COMPLEX_IMAG (z
);
7007 return scm_sum (scm_from_double (acos (0.0)),
7008 scm_product (scm_c_make_rectangular (0, 1),
7009 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
7012 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
7016 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
7018 "With one argument, compute the arc tangent of @var{z}.\n"
7019 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
7020 "using the sign of @var{z} and @var{y} to determine the quadrant.")
7021 #define FUNC_NAME s_scm_atan
7025 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
7026 return z
; /* atan(exact0) = exact0 */
7027 else if (scm_is_real (z
))
7028 return scm_from_double (atan (scm_to_double (z
)));
7029 else if (SCM_COMPLEXP (z
))
7032 v
= SCM_COMPLEX_REAL (z
);
7033 w
= SCM_COMPLEX_IMAG (z
);
7034 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
7035 scm_c_make_rectangular (v
, w
+ 1.0))),
7036 scm_c_make_rectangular (0, 2));
7039 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
7041 else if (scm_is_real (z
))
7043 if (scm_is_real (y
))
7044 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
7046 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
7049 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
7053 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
7055 "Compute the inverse hyperbolic sine of @var{z}.")
7056 #define FUNC_NAME s_scm_sys_asinh
7058 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
7059 return z
; /* asinh(exact0) = exact0 */
7060 else if (scm_is_real (z
))
7061 return scm_from_double (asinh (scm_to_double (z
)));
7062 else if (scm_is_number (z
))
7063 return scm_log (scm_sum (z
,
7064 scm_sqrt (scm_sum (scm_product (z
, z
),
7067 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
7071 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
7073 "Compute the inverse hyperbolic cosine of @var{z}.")
7074 #define FUNC_NAME s_scm_sys_acosh
7076 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM1
)))
7077 return SCM_INUM0
; /* acosh(exact1) = exact0 */
7078 else if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
7079 return scm_from_double (acosh (scm_to_double (z
)));
7080 else if (scm_is_number (z
))
7081 return scm_log (scm_sum (z
,
7082 scm_sqrt (scm_difference (scm_product (z
, z
),
7085 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
7089 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
7091 "Compute the inverse hyperbolic tangent of @var{z}.")
7092 #define FUNC_NAME s_scm_sys_atanh
7094 if (SCM_UNLIKELY (scm_is_eq (z
, SCM_INUM0
)))
7095 return z
; /* atanh(exact0) = exact0 */
7096 else if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
7097 return scm_from_double (atanh (scm_to_double (z
)));
7098 else if (scm_is_number (z
))
7099 return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1
, z
),
7100 scm_difference (SCM_INUM1
, z
))),
7103 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
7108 scm_c_make_rectangular (double re
, double im
)
7111 return scm_from_double (re
);
7116 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
7118 SCM_SET_CELL_TYPE (z
, scm_tc16_complex
);
7119 SCM_COMPLEX_REAL (z
) = re
;
7120 SCM_COMPLEX_IMAG (z
) = im
;
7125 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
7126 (SCM real_part
, SCM imaginary_part
),
7127 "Return a complex number constructed of the given @var{real-part} "
7128 "and @var{imaginary-part} parts.")
7129 #define FUNC_NAME s_scm_make_rectangular
7131 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
7132 SCM_ARG1
, FUNC_NAME
, "real");
7133 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
7134 SCM_ARG2
, FUNC_NAME
, "real");
7135 return scm_c_make_rectangular (scm_to_double (real_part
),
7136 scm_to_double (imaginary_part
));
7141 scm_c_make_polar (double mag
, double ang
)
7145 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
7146 use it on Glibc-based systems that have it (it's a GNU extension). See
7147 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
7149 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
7150 sincos (ang
, &s
, &c
);
7155 return scm_c_make_rectangular (mag
* c
, mag
* s
);
7158 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
7160 "Return the complex number @var{x} * e^(i * @var{y}).")
7161 #define FUNC_NAME s_scm_make_polar
7163 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
7164 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
7165 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
7170 SCM_PRIMITIVE_GENERIC (scm_real_part
, "real-part", 1, 0, 0,
7172 "Return the real part of the number @var{z}.")
7173 #define FUNC_NAME s_scm_real_part
7175 if (SCM_COMPLEXP (z
))
7176 return scm_from_double (SCM_COMPLEX_REAL (z
));
7177 else if (SCM_I_INUMP (z
) || SCM_BIGP (z
) || SCM_REALP (z
) || SCM_FRACTIONP (z
))
7180 SCM_WTA_DISPATCH_1 (g_scm_real_part
, z
, SCM_ARG1
, s_scm_real_part
);
7185 SCM_PRIMITIVE_GENERIC (scm_imag_part
, "imag-part", 1, 0, 0,
7187 "Return the imaginary part of the number @var{z}.")
7188 #define FUNC_NAME s_scm_imag_part
7190 if (SCM_COMPLEXP (z
))
7191 return scm_from_double (SCM_COMPLEX_IMAG (z
));
7192 else if (SCM_REALP (z
))
7194 else if (SCM_I_INUMP (z
) || SCM_BIGP (z
) || SCM_FRACTIONP (z
))
7197 SCM_WTA_DISPATCH_1 (g_scm_imag_part
, z
, SCM_ARG1
, s_scm_imag_part
);
7201 SCM_PRIMITIVE_GENERIC (scm_numerator
, "numerator", 1, 0, 0,
7203 "Return the numerator of the number @var{z}.")
7204 #define FUNC_NAME s_scm_numerator
7206 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7208 else if (SCM_FRACTIONP (z
))
7209 return SCM_FRACTION_NUMERATOR (z
);
7210 else if (SCM_REALP (z
))
7211 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
7213 SCM_WTA_DISPATCH_1 (g_scm_numerator
, z
, SCM_ARG1
, s_scm_numerator
);
7218 SCM_PRIMITIVE_GENERIC (scm_denominator
, "denominator", 1, 0, 0,
7220 "Return the denominator of the number @var{z}.")
7221 #define FUNC_NAME s_scm_denominator
7223 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7225 else if (SCM_FRACTIONP (z
))
7226 return SCM_FRACTION_DENOMINATOR (z
);
7227 else if (SCM_REALP (z
))
7228 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
7230 SCM_WTA_DISPATCH_1 (g_scm_denominator
, z
, SCM_ARG1
, s_scm_denominator
);
7235 SCM_PRIMITIVE_GENERIC (scm_magnitude
, "magnitude", 1, 0, 0,
7237 "Return the magnitude of the number @var{z}. This is the same as\n"
7238 "@code{abs} for real arguments, but also allows complex numbers.")
7239 #define FUNC_NAME s_scm_magnitude
7241 if (SCM_I_INUMP (z
))
7243 scm_t_inum zz
= SCM_I_INUM (z
);
7246 else if (SCM_POSFIXABLE (-zz
))
7247 return SCM_I_MAKINUM (-zz
);
7249 return scm_i_inum2big (-zz
);
7251 else if (SCM_BIGP (z
))
7253 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
7254 scm_remember_upto_here_1 (z
);
7256 return scm_i_clonebig (z
, 0);
7260 else if (SCM_REALP (z
))
7261 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
7262 else if (SCM_COMPLEXP (z
))
7263 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
7264 else if (SCM_FRACTIONP (z
))
7266 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
7268 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
7269 SCM_FRACTION_DENOMINATOR (z
));
7272 SCM_WTA_DISPATCH_1 (g_scm_magnitude
, z
, SCM_ARG1
, s_scm_magnitude
);
7277 SCM_PRIMITIVE_GENERIC (scm_angle
, "angle", 1, 0, 0,
7279 "Return the angle of the complex number @var{z}.")
7280 #define FUNC_NAME s_scm_angle
7282 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
7283 flo0 to save allocating a new flonum with scm_from_double each time.
7284 But if atan2 follows the floating point rounding mode, then the value
7285 is not a constant. Maybe it'd be close enough though. */
7286 if (SCM_I_INUMP (z
))
7288 if (SCM_I_INUM (z
) >= 0)
7291 return scm_from_double (atan2 (0.0, -1.0));
7293 else if (SCM_BIGP (z
))
7295 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
7296 scm_remember_upto_here_1 (z
);
7298 return scm_from_double (atan2 (0.0, -1.0));
7302 else if (SCM_REALP (z
))
7304 if (SCM_REAL_VALUE (z
) >= 0)
7307 return scm_from_double (atan2 (0.0, -1.0));
7309 else if (SCM_COMPLEXP (z
))
7310 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
7311 else if (SCM_FRACTIONP (z
))
7313 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
7315 else return scm_from_double (atan2 (0.0, -1.0));
7318 SCM_WTA_DISPATCH_1 (g_scm_angle
, z
, SCM_ARG1
, s_scm_angle
);
7323 SCM_PRIMITIVE_GENERIC (scm_exact_to_inexact
, "exact->inexact", 1, 0, 0,
7325 "Convert the number @var{z} to its inexact representation.\n")
7326 #define FUNC_NAME s_scm_exact_to_inexact
7328 if (SCM_I_INUMP (z
))
7329 return scm_from_double ((double) SCM_I_INUM (z
));
7330 else if (SCM_BIGP (z
))
7331 return scm_from_double (scm_i_big2dbl (z
));
7332 else if (SCM_FRACTIONP (z
))
7333 return scm_from_double (scm_i_fraction2double (z
));
7334 else if (SCM_INEXACTP (z
))
7337 SCM_WTA_DISPATCH_1 (g_scm_exact_to_inexact
, z
, 1, s_scm_exact_to_inexact
);
7342 SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
7344 "Return an exact number that is numerically closest to @var{z}.")
7345 #define FUNC_NAME s_scm_inexact_to_exact
7347 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7349 else if (SCM_REALP (z
))
7351 if (!DOUBLE_IS_FINITE (SCM_REAL_VALUE (z
)))
7352 SCM_OUT_OF_RANGE (1, z
);
7359 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
7360 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
7361 scm_i_mpz2num (mpq_denref (frac
)));
7363 /* When scm_i_make_ratio throws, we leak the memory allocated
7370 else if (SCM_FRACTIONP (z
))
7373 SCM_WTA_DISPATCH_1 (g_scm_inexact_to_exact
, z
, 1, s_scm_inexact_to_exact
);
7377 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
7379 "Returns the @emph{simplest} rational number differing\n"
7380 "from @var{x} by no more than @var{eps}.\n"
7382 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
7383 "exact result when both its arguments are exact. Thus, you might need\n"
7384 "to use @code{inexact->exact} on the arguments.\n"
7387 "(rationalize (inexact->exact 1.2) 1/100)\n"
7390 #define FUNC_NAME s_scm_rationalize
7392 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
7393 SCM_ASSERT_TYPE (scm_is_real (eps
), eps
, SCM_ARG2
, FUNC_NAME
, "real");
7394 eps
= scm_abs (eps
);
7395 if (scm_is_false (scm_positive_p (eps
)))
7397 /* eps is either zero or a NaN */
7398 if (scm_is_true (scm_nan_p (eps
)))
7400 else if (SCM_INEXACTP (eps
))
7401 return scm_exact_to_inexact (x
);
7405 else if (scm_is_false (scm_finite_p (eps
)))
7407 if (scm_is_true (scm_finite_p (x
)))
7412 else if (scm_is_false (scm_finite_p (x
))) /* checks for both inf and nan */
7414 else if (scm_is_false (scm_less_p (scm_floor (scm_sum (x
, eps
)),
7415 scm_ceiling (scm_difference (x
, eps
)))))
7417 /* There's an integer within range; we want the one closest to zero */
7418 if (scm_is_false (scm_less_p (eps
, scm_abs (x
))))
7420 /* zero is within range */
7421 if (SCM_INEXACTP (x
) || SCM_INEXACTP (eps
))
7426 else if (scm_is_true (scm_positive_p (x
)))
7427 return scm_ceiling (scm_difference (x
, eps
));
7429 return scm_floor (scm_sum (x
, eps
));
7433 /* Use continued fractions to find closest ratio. All
7434 arithmetic is done with exact numbers.
7437 SCM ex
= scm_inexact_to_exact (x
);
7438 SCM int_part
= scm_floor (ex
);
7440 SCM a1
= SCM_INUM0
, a2
= SCM_INUM1
, a
= SCM_INUM0
;
7441 SCM b1
= SCM_INUM1
, b2
= SCM_INUM0
, b
= SCM_INUM0
;
7445 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
7446 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
7448 /* We stop after a million iterations just to be absolutely sure
7449 that we don't go into an infinite loop. The process normally
7450 converges after less than a dozen iterations.
7453 while (++i
< 1000000)
7455 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
7456 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
7457 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
7459 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
7460 eps
))) /* abs(x-a/b) <= eps */
7462 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
7463 if (SCM_INEXACTP (x
) || SCM_INEXACTP (eps
))
7464 return scm_exact_to_inexact (res
);
7468 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
7470 tt
= scm_floor (rx
); /* tt = floor (rx) */
7476 scm_num_overflow (s_scm_rationalize
);
7481 /* conversion functions */
7484 scm_is_integer (SCM val
)
7486 return scm_is_true (scm_integer_p (val
));
7490 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
7492 if (SCM_I_INUMP (val
))
7494 scm_t_signed_bits n
= SCM_I_INUM (val
);
7495 return n
>= min
&& n
<= max
;
7497 else if (SCM_BIGP (val
))
7499 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
7501 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
7503 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
7505 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
7506 return n
>= min
&& n
<= max
;
7516 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
7517 > CHAR_BIT
*sizeof (scm_t_uintmax
))
7520 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
7521 SCM_I_BIG_MPZ (val
));
7523 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
7535 return n
>= min
&& n
<= max
;
7543 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
7545 if (SCM_I_INUMP (val
))
7547 scm_t_signed_bits n
= SCM_I_INUM (val
);
7548 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
7550 else if (SCM_BIGP (val
))
7552 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
7554 else if (max
<= ULONG_MAX
)
7556 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
7558 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
7559 return n
>= min
&& n
<= max
;
7569 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
7572 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
7573 > CHAR_BIT
*sizeof (scm_t_uintmax
))
7576 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
7577 SCM_I_BIG_MPZ (val
));
7579 return n
>= min
&& n
<= max
;
7587 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
7589 scm_error (scm_out_of_range_key
,
7591 "Value out of range ~S to ~S: ~S",
7592 scm_list_3 (min
, max
, bad_val
),
7593 scm_list_1 (bad_val
));
7596 #define TYPE scm_t_intmax
7597 #define TYPE_MIN min
7598 #define TYPE_MAX max
7599 #define SIZEOF_TYPE 0
7600 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
7601 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
7602 #include "libguile/conv-integer.i.c"
7604 #define TYPE scm_t_uintmax
7605 #define TYPE_MIN min
7606 #define TYPE_MAX max
7607 #define SIZEOF_TYPE 0
7608 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
7609 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
7610 #include "libguile/conv-uinteger.i.c"
7612 #define TYPE scm_t_int8
7613 #define TYPE_MIN SCM_T_INT8_MIN
7614 #define TYPE_MAX SCM_T_INT8_MAX
7615 #define SIZEOF_TYPE 1
7616 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
7617 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
7618 #include "libguile/conv-integer.i.c"
7620 #define TYPE scm_t_uint8
7622 #define TYPE_MAX SCM_T_UINT8_MAX
7623 #define SIZEOF_TYPE 1
7624 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
7625 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
7626 #include "libguile/conv-uinteger.i.c"
7628 #define TYPE scm_t_int16
7629 #define TYPE_MIN SCM_T_INT16_MIN
7630 #define TYPE_MAX SCM_T_INT16_MAX
7631 #define SIZEOF_TYPE 2
7632 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
7633 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
7634 #include "libguile/conv-integer.i.c"
7636 #define TYPE scm_t_uint16
7638 #define TYPE_MAX SCM_T_UINT16_MAX
7639 #define SIZEOF_TYPE 2
7640 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
7641 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
7642 #include "libguile/conv-uinteger.i.c"
7644 #define TYPE scm_t_int32
7645 #define TYPE_MIN SCM_T_INT32_MIN
7646 #define TYPE_MAX SCM_T_INT32_MAX
7647 #define SIZEOF_TYPE 4
7648 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
7649 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
7650 #include "libguile/conv-integer.i.c"
7652 #define TYPE scm_t_uint32
7654 #define TYPE_MAX SCM_T_UINT32_MAX
7655 #define SIZEOF_TYPE 4
7656 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
7657 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
7658 #include "libguile/conv-uinteger.i.c"
7660 #define TYPE scm_t_wchar
7661 #define TYPE_MIN (scm_t_int32)-1
7662 #define TYPE_MAX (scm_t_int32)0x10ffff
7663 #define SIZEOF_TYPE 4
7664 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
7665 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
7666 #include "libguile/conv-integer.i.c"
7668 #define TYPE scm_t_int64
7669 #define TYPE_MIN SCM_T_INT64_MIN
7670 #define TYPE_MAX SCM_T_INT64_MAX
7671 #define SIZEOF_TYPE 8
7672 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
7673 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
7674 #include "libguile/conv-integer.i.c"
7676 #define TYPE scm_t_uint64
7678 #define TYPE_MAX SCM_T_UINT64_MAX
7679 #define SIZEOF_TYPE 8
7680 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
7681 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
7682 #include "libguile/conv-uinteger.i.c"
7685 scm_to_mpz (SCM val
, mpz_t rop
)
7687 if (SCM_I_INUMP (val
))
7688 mpz_set_si (rop
, SCM_I_INUM (val
));
7689 else if (SCM_BIGP (val
))
7690 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
7692 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
7696 scm_from_mpz (mpz_t val
)
7698 return scm_i_mpz2num (val
);
7702 scm_is_real (SCM val
)
7704 return scm_is_true (scm_real_p (val
));
7708 scm_is_rational (SCM val
)
7710 return scm_is_true (scm_rational_p (val
));
7714 scm_to_double (SCM val
)
7716 if (SCM_I_INUMP (val
))
7717 return SCM_I_INUM (val
);
7718 else if (SCM_BIGP (val
))
7719 return scm_i_big2dbl (val
);
7720 else if (SCM_FRACTIONP (val
))
7721 return scm_i_fraction2double (val
);
7722 else if (SCM_REALP (val
))
7723 return SCM_REAL_VALUE (val
);
7725 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
7729 scm_from_double (double val
)
7733 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
7735 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
7736 SCM_REAL_VALUE (z
) = val
;
7741 #if SCM_ENABLE_DEPRECATED == 1
7744 scm_num2float (SCM num
, unsigned long pos
, const char *s_caller
)
7746 scm_c_issue_deprecation_warning
7747 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
7751 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
7755 scm_out_of_range (NULL
, num
);
7758 return scm_to_double (num
);
7762 scm_num2double (SCM num
, unsigned long pos
, const char *s_caller
)
7764 scm_c_issue_deprecation_warning
7765 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
7769 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
7773 scm_out_of_range (NULL
, num
);
7776 return scm_to_double (num
);
7782 scm_is_complex (SCM val
)
7784 return scm_is_true (scm_complex_p (val
));
7788 scm_c_real_part (SCM z
)
7790 if (SCM_COMPLEXP (z
))
7791 return SCM_COMPLEX_REAL (z
);
7794 /* Use the scm_real_part to get proper error checking and
7797 return scm_to_double (scm_real_part (z
));
7802 scm_c_imag_part (SCM z
)
7804 if (SCM_COMPLEXP (z
))
7805 return SCM_COMPLEX_IMAG (z
);
7808 /* Use the scm_imag_part to get proper error checking and
7809 dispatching. The result will almost always be 0.0, but not
7812 return scm_to_double (scm_imag_part (z
));
7817 scm_c_magnitude (SCM z
)
7819 return scm_to_double (scm_magnitude (z
));
7825 return scm_to_double (scm_angle (z
));
7829 scm_is_number (SCM z
)
7831 return scm_is_true (scm_number_p (z
));
7835 /* In the following functions we dispatch to the real-arg funcs like log()
7836 when we know the arg is real, instead of just handing everything to
7837 clog() for instance. This is in case clog() doesn't optimize for a
7838 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
7839 well use it to go straight to the applicable C func. */
7841 SCM_PRIMITIVE_GENERIC (scm_log
, "log", 1, 0, 0,
7843 "Return the natural logarithm of @var{z}.")
7844 #define FUNC_NAME s_scm_log
7846 if (SCM_COMPLEXP (z
))
7848 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
7849 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
7851 double re
= SCM_COMPLEX_REAL (z
);
7852 double im
= SCM_COMPLEX_IMAG (z
);
7853 return scm_c_make_rectangular (log (hypot (re
, im
)),
7857 else if (SCM_NUMBERP (z
))
7859 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
7860 although the value itself overflows. */
7861 double re
= scm_to_double (z
);
7862 double l
= log (fabs (re
));
7864 return scm_from_double (l
);
7866 return scm_c_make_rectangular (l
, M_PI
);
7869 SCM_WTA_DISPATCH_1 (g_scm_log
, z
, 1, s_scm_log
);
7874 SCM_PRIMITIVE_GENERIC (scm_log10
, "log10", 1, 0, 0,
7876 "Return the base 10 logarithm of @var{z}.")
7877 #define FUNC_NAME s_scm_log10
7879 if (SCM_COMPLEXP (z
))
7881 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
7882 clog() and a multiply by M_LOG10E, rather than the fallback
7883 log10+hypot+atan2.) */
7884 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
7885 && defined SCM_COMPLEX_VALUE
7886 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
7888 double re
= SCM_COMPLEX_REAL (z
);
7889 double im
= SCM_COMPLEX_IMAG (z
);
7890 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
7891 M_LOG10E
* atan2 (im
, re
));
7894 else if (SCM_NUMBERP (z
))
7896 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
7897 although the value itself overflows. */
7898 double re
= scm_to_double (z
);
7899 double l
= log10 (fabs (re
));
7901 return scm_from_double (l
);
7903 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
7906 SCM_WTA_DISPATCH_1 (g_scm_log10
, z
, 1, s_scm_log10
);
7911 SCM_PRIMITIVE_GENERIC (scm_exp
, "exp", 1, 0, 0,
7913 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
7914 "base of natural logarithms (2.71828@dots{}).")
7915 #define FUNC_NAME s_scm_exp
7917 if (SCM_COMPLEXP (z
))
7919 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
7920 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
7922 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
7923 SCM_COMPLEX_IMAG (z
));
7926 else if (SCM_NUMBERP (z
))
7928 /* When z is a negative bignum the conversion to double overflows,
7929 giving -infinity, but that's ok, the exp is still 0.0. */
7930 return scm_from_double (exp (scm_to_double (z
)));
7933 SCM_WTA_DISPATCH_1 (g_scm_exp
, z
, 1, s_scm_exp
);
7938 SCM_PRIMITIVE_GENERIC (scm_sqrt
, "sqrt", 1, 0, 0,
7940 "Return the square root of @var{z}. Of the two possible roots\n"
7941 "(positive and negative), the one with the a positive real part\n"
7942 "is returned, or if that's zero then a positive imaginary part.\n"
7946 "(sqrt 9.0) @result{} 3.0\n"
7947 "(sqrt -9.0) @result{} 0.0+3.0i\n"
7948 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
7949 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
7951 #define FUNC_NAME s_scm_sqrt
7953 if (SCM_COMPLEXP (z
))
7955 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
7956 && defined SCM_COMPLEX_VALUE
7957 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (z
)));
7959 double re
= SCM_COMPLEX_REAL (z
);
7960 double im
= SCM_COMPLEX_IMAG (z
);
7961 return scm_c_make_polar (sqrt (hypot (re
, im
)),
7962 0.5 * atan2 (im
, re
));
7965 else if (SCM_NUMBERP (z
))
7967 double xx
= scm_to_double (z
);
7969 return scm_c_make_rectangular (0.0, sqrt (-xx
));
7971 return scm_from_double (sqrt (xx
));
7974 SCM_WTA_DISPATCH_1 (g_scm_sqrt
, z
, 1, s_scm_sqrt
);
7985 mpz_init_set_si (z_negative_one
, -1);
7987 /* It may be possible to tune the performance of some algorithms by using
7988 * the following constants to avoid the creation of bignums. Please, before
7989 * using these values, remember the two rules of program optimization:
7990 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
7991 scm_c_define ("most-positive-fixnum",
7992 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
7993 scm_c_define ("most-negative-fixnum",
7994 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
7996 scm_add_feature ("complex");
7997 scm_add_feature ("inexact");
7998 flo0
= scm_from_double (0.0);
8000 /* determine floating point precision */
8001 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
8003 init_dblprec(&scm_dblprec
[i
-2],i
);
8004 init_fx_radix(fx_per_radix
[i
-2],i
);
8007 /* hard code precision for base 10 if the preprocessor tells us to... */
8008 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
8011 exactly_one_half
= scm_divide (SCM_INUM1
, SCM_I_MAKINUM (2));
8012 #include "libguile/numbers.x"