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_BIGP (x
))
750 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
752 return scm_i_clonebig (x
, 0);
756 else if (SCM_REALP (x
))
758 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
759 double xx
= SCM_REAL_VALUE (x
);
761 return scm_from_double (-xx
);
765 else if (SCM_FRACTIONP (x
))
767 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
769 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
770 SCM_FRACTION_DENOMINATOR (x
));
773 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
778 SCM_PRIMITIVE_GENERIC (scm_quotient
, "quotient", 2, 0, 0,
780 "Return the quotient of the numbers @var{x} and @var{y}.")
781 #define FUNC_NAME s_scm_quotient
783 if (SCM_LIKELY (SCM_I_INUMP (x
)))
785 scm_t_inum xx
= SCM_I_INUM (x
);
786 if (SCM_LIKELY (SCM_I_INUMP (y
)))
788 scm_t_inum yy
= SCM_I_INUM (y
);
789 if (SCM_UNLIKELY (yy
== 0))
790 scm_num_overflow (s_scm_quotient
);
793 scm_t_inum z
= xx
/ yy
;
794 if (SCM_LIKELY (SCM_FIXABLE (z
)))
795 return SCM_I_MAKINUM (z
);
797 return scm_i_inum2big (z
);
800 else if (SCM_BIGP (y
))
802 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
803 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
804 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
806 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
807 scm_remember_upto_here_1 (y
);
808 return SCM_I_MAKINUM (-1);
814 SCM_WTA_DISPATCH_2 (g_scm_quotient
, x
, y
, SCM_ARG2
, s_scm_quotient
);
816 else if (SCM_BIGP (x
))
818 if (SCM_LIKELY (SCM_I_INUMP (y
)))
820 scm_t_inum yy
= SCM_I_INUM (y
);
821 if (SCM_UNLIKELY (yy
== 0))
822 scm_num_overflow (s_scm_quotient
);
823 else if (SCM_UNLIKELY (yy
== 1))
827 SCM result
= scm_i_mkbig ();
830 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
833 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
836 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
837 scm_remember_upto_here_1 (x
);
838 return scm_i_normbig (result
);
841 else if (SCM_BIGP (y
))
843 SCM result
= scm_i_mkbig ();
844 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
847 scm_remember_upto_here_2 (x
, y
);
848 return scm_i_normbig (result
);
851 SCM_WTA_DISPATCH_2 (g_scm_quotient
, x
, y
, SCM_ARG2
, s_scm_quotient
);
854 SCM_WTA_DISPATCH_2 (g_scm_quotient
, x
, y
, SCM_ARG1
, s_scm_quotient
);
858 SCM_PRIMITIVE_GENERIC (scm_remainder
, "remainder", 2, 0, 0,
860 "Return the remainder of the numbers @var{x} and @var{y}.\n"
862 "(remainder 13 4) @result{} 1\n"
863 "(remainder -13 4) @result{} -1\n"
865 #define FUNC_NAME s_scm_remainder
867 if (SCM_LIKELY (SCM_I_INUMP (x
)))
869 if (SCM_LIKELY (SCM_I_INUMP (y
)))
871 scm_t_inum yy
= SCM_I_INUM (y
);
872 if (SCM_UNLIKELY (yy
== 0))
873 scm_num_overflow (s_scm_remainder
);
876 /* C99 specifies that "%" is the remainder corresponding to a
877 quotient rounded towards zero, and that's also traditional
878 for machine division, so z here should be well defined. */
879 scm_t_inum z
= SCM_I_INUM (x
) % yy
;
880 return SCM_I_MAKINUM (z
);
883 else if (SCM_BIGP (y
))
885 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
886 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
887 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
889 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
890 scm_remember_upto_here_1 (y
);
897 SCM_WTA_DISPATCH_2 (g_scm_remainder
, x
, y
, SCM_ARG2
, s_scm_remainder
);
899 else if (SCM_BIGP (x
))
901 if (SCM_LIKELY (SCM_I_INUMP (y
)))
903 scm_t_inum yy
= SCM_I_INUM (y
);
904 if (SCM_UNLIKELY (yy
== 0))
905 scm_num_overflow (s_scm_remainder
);
908 SCM result
= scm_i_mkbig ();
911 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
912 scm_remember_upto_here_1 (x
);
913 return scm_i_normbig (result
);
916 else if (SCM_BIGP (y
))
918 SCM result
= scm_i_mkbig ();
919 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
922 scm_remember_upto_here_2 (x
, y
);
923 return scm_i_normbig (result
);
926 SCM_WTA_DISPATCH_2 (g_scm_remainder
, x
, y
, SCM_ARG2
, s_scm_remainder
);
929 SCM_WTA_DISPATCH_2 (g_scm_remainder
, x
, y
, SCM_ARG1
, s_scm_remainder
);
934 SCM_PRIMITIVE_GENERIC (scm_modulo
, "modulo", 2, 0, 0,
936 "Return the modulo of the numbers @var{x} and @var{y}.\n"
938 "(modulo 13 4) @result{} 1\n"
939 "(modulo -13 4) @result{} 3\n"
941 #define FUNC_NAME s_scm_modulo
943 if (SCM_LIKELY (SCM_I_INUMP (x
)))
945 scm_t_inum xx
= SCM_I_INUM (x
);
946 if (SCM_LIKELY (SCM_I_INUMP (y
)))
948 scm_t_inum yy
= SCM_I_INUM (y
);
949 if (SCM_UNLIKELY (yy
== 0))
950 scm_num_overflow (s_scm_modulo
);
953 /* C99 specifies that "%" is the remainder corresponding to a
954 quotient rounded towards zero, and that's also traditional
955 for machine division, so z here should be well defined. */
956 scm_t_inum z
= xx
% yy
;
973 return SCM_I_MAKINUM (result
);
976 else if (SCM_BIGP (y
))
978 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
985 SCM pos_y
= scm_i_clonebig (y
, 0);
986 /* do this after the last scm_op */
987 mpz_init_set_si (z_x
, xx
);
988 result
= pos_y
; /* re-use this bignum */
989 mpz_mod (SCM_I_BIG_MPZ (result
),
991 SCM_I_BIG_MPZ (pos_y
));
992 scm_remember_upto_here_1 (pos_y
);
996 result
= scm_i_mkbig ();
997 /* do this after the last scm_op */
998 mpz_init_set_si (z_x
, xx
);
999 mpz_mod (SCM_I_BIG_MPZ (result
),
1002 scm_remember_upto_here_1 (y
);
1005 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1006 mpz_add (SCM_I_BIG_MPZ (result
),
1008 SCM_I_BIG_MPZ (result
));
1009 scm_remember_upto_here_1 (y
);
1010 /* and do this before the next one */
1012 return scm_i_normbig (result
);
1016 SCM_WTA_DISPATCH_2 (g_scm_modulo
, x
, y
, SCM_ARG2
, s_scm_modulo
);
1018 else if (SCM_BIGP (x
))
1020 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1022 scm_t_inum yy
= SCM_I_INUM (y
);
1023 if (SCM_UNLIKELY (yy
== 0))
1024 scm_num_overflow (s_scm_modulo
);
1027 SCM result
= scm_i_mkbig ();
1028 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
1030 (yy
< 0) ? - yy
: yy
);
1031 scm_remember_upto_here_1 (x
);
1032 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1033 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
1034 SCM_I_BIG_MPZ (result
),
1036 return scm_i_normbig (result
);
1039 else if (SCM_BIGP (y
))
1041 SCM result
= scm_i_mkbig ();
1042 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
1043 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1044 mpz_mod (SCM_I_BIG_MPZ (result
),
1046 SCM_I_BIG_MPZ (pos_y
));
1048 scm_remember_upto_here_1 (x
);
1049 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1050 mpz_add (SCM_I_BIG_MPZ (result
),
1052 SCM_I_BIG_MPZ (result
));
1053 scm_remember_upto_here_2 (y
, pos_y
);
1054 return scm_i_normbig (result
);
1057 SCM_WTA_DISPATCH_2 (g_scm_modulo
, x
, y
, SCM_ARG2
, s_scm_modulo
);
1060 SCM_WTA_DISPATCH_2 (g_scm_modulo
, x
, y
, SCM_ARG1
, s_scm_modulo
);
1064 static SCM
scm_i_inexact_euclidean_quotient (double x
, double y
);
1065 static SCM
scm_i_slow_exact_euclidean_quotient (SCM x
, SCM y
);
1067 SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient
, "euclidean-quotient", 2, 0, 0,
1069 "Return the integer @var{q} such that\n"
1070 "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1071 "where @math{0 <= @var{r} < abs(@var{y})}.\n"
1073 "(euclidean-quotient 123 10) @result{} 12\n"
1074 "(euclidean-quotient 123 -10) @result{} -12\n"
1075 "(euclidean-quotient -123 10) @result{} -13\n"
1076 "(euclidean-quotient -123 -10) @result{} 13\n"
1077 "(euclidean-quotient -123.2 -63.5) @result{} 2.0\n"
1078 "(euclidean-quotient 16/3 -10/7) @result{} -3\n"
1080 #define FUNC_NAME s_scm_euclidean_quotient
1082 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1084 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1086 scm_t_inum yy
= SCM_I_INUM (y
);
1087 if (SCM_UNLIKELY (yy
== 0))
1088 scm_num_overflow (s_scm_euclidean_quotient
);
1091 scm_t_inum xx
= SCM_I_INUM (x
);
1092 scm_t_inum qq
= xx
/ yy
;
1100 return SCM_I_MAKINUM (qq
);
1103 else if (SCM_BIGP (y
))
1105 if (SCM_I_INUM (x
) >= 0)
1108 return SCM_I_MAKINUM (- mpz_sgn (SCM_I_BIG_MPZ (y
)));
1110 else if (SCM_REALP (y
))
1111 return scm_i_inexact_euclidean_quotient
1112 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1113 else if (SCM_FRACTIONP (y
))
1114 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1116 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1117 s_scm_euclidean_quotient
);
1119 else if (SCM_BIGP (x
))
1121 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1123 scm_t_inum yy
= SCM_I_INUM (y
);
1124 if (SCM_UNLIKELY (yy
== 0))
1125 scm_num_overflow (s_scm_euclidean_quotient
);
1128 SCM q
= scm_i_mkbig ();
1130 mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (x
), yy
);
1133 mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (x
), -yy
);
1134 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1136 scm_remember_upto_here_1 (x
);
1137 return scm_i_normbig (q
);
1140 else if (SCM_BIGP (y
))
1142 SCM q
= scm_i_mkbig ();
1143 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1144 mpz_fdiv_q (SCM_I_BIG_MPZ (q
),
1148 mpz_cdiv_q (SCM_I_BIG_MPZ (q
),
1151 scm_remember_upto_here_2 (x
, y
);
1152 return scm_i_normbig (q
);
1154 else if (SCM_REALP (y
))
1155 return scm_i_inexact_euclidean_quotient
1156 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1157 else if (SCM_FRACTIONP (y
))
1158 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1160 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1161 s_scm_euclidean_quotient
);
1163 else if (SCM_REALP (x
))
1165 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1166 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1167 return scm_i_inexact_euclidean_quotient
1168 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1170 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1171 s_scm_euclidean_quotient
);
1173 else if (SCM_FRACTIONP (x
))
1176 return scm_i_inexact_euclidean_quotient
1177 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1179 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1182 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG1
,
1183 s_scm_euclidean_quotient
);
1188 scm_i_inexact_euclidean_quotient (double x
, double y
)
1190 if (SCM_LIKELY (y
> 0))
1191 return scm_from_double (floor (x
/ y
));
1192 else if (SCM_LIKELY (y
< 0))
1193 return scm_from_double (ceil (x
/ y
));
1195 scm_num_overflow (s_scm_euclidean_quotient
); /* or return a NaN? */
1200 /* Compute exact euclidean_quotient the slow way.
1201 We use this only if both arguments are exact,
1202 and at least one of them is a fraction */
1204 scm_i_slow_exact_euclidean_quotient (SCM x
, SCM y
)
1206 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1207 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG1
,
1208 s_scm_euclidean_quotient
);
1209 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1210 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1211 s_scm_euclidean_quotient
);
1212 else if (scm_is_true (scm_positive_p (y
)))
1213 return scm_floor (scm_divide (x
, y
));
1214 else if (scm_is_true (scm_negative_p (y
)))
1215 return scm_ceiling (scm_divide (x
, y
));
1217 scm_num_overflow (s_scm_euclidean_quotient
);
1220 static SCM
scm_i_inexact_euclidean_remainder (double x
, double y
);
1221 static SCM
scm_i_slow_exact_euclidean_remainder (SCM x
, SCM y
);
1223 SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder
, "euclidean-remainder", 2, 0, 0,
1225 "Return the real number @var{r} such that\n"
1226 "@math{0 <= @var{r} < abs(@var{y})} and\n"
1227 "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1228 "for some integer @var{q}.\n"
1230 "(euclidean-remainder 123 10) @result{} 3\n"
1231 "(euclidean-remainder 123 -10) @result{} 3\n"
1232 "(euclidean-remainder -123 10) @result{} 7\n"
1233 "(euclidean-remainder -123 -10) @result{} 7\n"
1234 "(euclidean-remainder -123.2 -63.5) @result{} 3.8\n"
1235 "(euclidean-remainder 16/3 -10/7) @result{} 22/21\n"
1237 #define FUNC_NAME s_scm_euclidean_remainder
1239 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1241 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1243 scm_t_inum yy
= SCM_I_INUM (y
);
1244 if (SCM_UNLIKELY (yy
== 0))
1245 scm_num_overflow (s_scm_euclidean_remainder
);
1248 scm_t_inum rr
= SCM_I_INUM (x
) % yy
;
1250 return SCM_I_MAKINUM (rr
);
1252 return SCM_I_MAKINUM (rr
+ yy
);
1254 return SCM_I_MAKINUM (rr
- yy
);
1257 else if (SCM_BIGP (y
))
1259 scm_t_inum xx
= SCM_I_INUM (x
);
1262 else if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1264 SCM r
= scm_i_mkbig ();
1265 mpz_sub_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1266 scm_remember_upto_here_1 (y
);
1267 return scm_i_normbig (r
);
1271 SCM r
= scm_i_mkbig ();
1272 mpz_add_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1273 scm_remember_upto_here_1 (y
);
1274 mpz_neg (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (r
));
1275 return scm_i_normbig (r
);
1278 else if (SCM_REALP (y
))
1279 return scm_i_inexact_euclidean_remainder
1280 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1281 else if (SCM_FRACTIONP (y
))
1282 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1284 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1285 s_scm_euclidean_remainder
);
1287 else if (SCM_BIGP (x
))
1289 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1291 scm_t_inum yy
= SCM_I_INUM (y
);
1292 if (SCM_UNLIKELY (yy
== 0))
1293 scm_num_overflow (s_scm_euclidean_remainder
);
1299 rr
= mpz_fdiv_ui (SCM_I_BIG_MPZ (x
), yy
);
1300 scm_remember_upto_here_1 (x
);
1301 return SCM_I_MAKINUM (rr
);
1304 else if (SCM_BIGP (y
))
1306 SCM r
= scm_i_mkbig ();
1307 mpz_mod (SCM_I_BIG_MPZ (r
),
1310 scm_remember_upto_here_2 (x
, y
);
1311 return scm_i_normbig (r
);
1313 else if (SCM_REALP (y
))
1314 return scm_i_inexact_euclidean_remainder
1315 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1316 else if (SCM_FRACTIONP (y
))
1317 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1319 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1320 s_scm_euclidean_remainder
);
1322 else if (SCM_REALP (x
))
1324 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1325 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1326 return scm_i_inexact_euclidean_remainder
1327 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1329 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1330 s_scm_euclidean_remainder
);
1332 else if (SCM_FRACTIONP (x
))
1335 return scm_i_inexact_euclidean_remainder
1336 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1338 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1341 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG1
,
1342 s_scm_euclidean_remainder
);
1347 scm_i_inexact_euclidean_remainder (double x
, double y
)
1351 /* Although it would be more efficient to use fmod here, we can't
1352 because it would in some cases produce results inconsistent with
1353 scm_i_inexact_euclidean_quotient, such that x != q * y + r (not
1354 even close). In particular, when x is very close to a multiple of
1355 y, then r might be either 0.0 or abs(y)-epsilon, but those two
1356 cases must correspond to different choices of q. If r = 0.0 then q
1357 must be x/y, and if r = abs(y) then q must be (x-r)/y. If quotient
1358 chooses one and remainder chooses the other, it would be bad. This
1359 problem was observed with x = 130.0 and y = 10/7. */
1360 if (SCM_LIKELY (y
> 0))
1362 else if (SCM_LIKELY (y
< 0))
1365 scm_num_overflow (s_scm_euclidean_remainder
); /* or return a NaN? */
1368 return scm_from_double (x
- q
* y
);
1371 /* Compute exact euclidean_remainder the slow way.
1372 We use this only if both arguments are exact,
1373 and at least one of them is a fraction */
1375 scm_i_slow_exact_euclidean_remainder (SCM x
, SCM y
)
1379 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1380 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG1
,
1381 s_scm_euclidean_remainder
);
1382 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1383 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1384 s_scm_euclidean_remainder
);
1385 else if (scm_is_true (scm_positive_p (y
)))
1386 q
= scm_floor (scm_divide (x
, y
));
1387 else if (scm_is_true (scm_negative_p (y
)))
1388 q
= scm_ceiling (scm_divide (x
, y
));
1390 scm_num_overflow (s_scm_euclidean_remainder
);
1391 return scm_difference (x
, scm_product (y
, q
));
1395 static SCM
scm_i_inexact_euclidean_divide (double x
, double y
);
1396 static SCM
scm_i_slow_exact_euclidean_divide (SCM x
, SCM y
);
1398 SCM_PRIMITIVE_GENERIC (scm_euclidean_divide
, "euclidean/", 2, 0, 0,
1400 "Return the integer @var{q} and the real number @var{r}\n"
1401 "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1402 "and @math{0 <= @var{r} < abs(@var{y})}.\n"
1404 "(euclidean/ 123 10) @result{} 12 and 3\n"
1405 "(euclidean/ 123 -10) @result{} -12 and 3\n"
1406 "(euclidean/ -123 10) @result{} -13 and 7\n"
1407 "(euclidean/ -123 -10) @result{} 13 and 7\n"
1408 "(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
1409 "(euclidean/ 16/3 -10/7) @result{} -3 and 22/21\n"
1411 #define FUNC_NAME s_scm_euclidean_divide
1413 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1415 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1417 scm_t_inum yy
= SCM_I_INUM (y
);
1418 if (SCM_UNLIKELY (yy
== 0))
1419 scm_num_overflow (s_scm_euclidean_divide
);
1422 scm_t_inum xx
= SCM_I_INUM (x
);
1423 scm_t_inum qq
= xx
/ yy
;
1424 scm_t_inum rr
= xx
- qq
* yy
;
1432 return scm_values (scm_list_2 (SCM_I_MAKINUM (qq
),
1433 SCM_I_MAKINUM (rr
)));
1436 else if (SCM_BIGP (y
))
1438 scm_t_inum xx
= SCM_I_INUM (x
);
1440 return scm_values (scm_list_2 (SCM_INUM0
, x
));
1441 else if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1443 SCM r
= scm_i_mkbig ();
1444 mpz_sub_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1445 scm_remember_upto_here_1 (y
);
1447 (scm_list_2 (SCM_I_MAKINUM (-1), scm_i_normbig (r
)));
1451 SCM r
= scm_i_mkbig ();
1452 mpz_add_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1453 scm_remember_upto_here_1 (y
);
1454 mpz_neg (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (r
));
1455 return scm_values (scm_list_2 (SCM_INUM1
, scm_i_normbig (r
)));
1458 else if (SCM_REALP (y
))
1459 return scm_i_inexact_euclidean_divide
1460 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1461 else if (SCM_FRACTIONP (y
))
1462 return scm_i_slow_exact_euclidean_divide (x
, y
);
1464 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1465 s_scm_euclidean_divide
);
1467 else if (SCM_BIGP (x
))
1469 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1471 scm_t_inum yy
= SCM_I_INUM (y
);
1472 if (SCM_UNLIKELY (yy
== 0))
1473 scm_num_overflow (s_scm_euclidean_divide
);
1476 SCM q
= scm_i_mkbig ();
1477 SCM r
= scm_i_mkbig ();
1479 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1480 SCM_I_BIG_MPZ (x
), yy
);
1483 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1484 SCM_I_BIG_MPZ (x
), -yy
);
1485 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1487 scm_remember_upto_here_1 (x
);
1488 return scm_values (scm_list_2 (scm_i_normbig (q
),
1489 scm_i_normbig (r
)));
1492 else if (SCM_BIGP (y
))
1494 SCM q
= scm_i_mkbig ();
1495 SCM r
= scm_i_mkbig ();
1496 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1497 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1498 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1500 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1501 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1502 scm_remember_upto_here_2 (x
, y
);
1503 return scm_values (scm_list_2 (scm_i_normbig (q
),
1504 scm_i_normbig (r
)));
1506 else if (SCM_REALP (y
))
1507 return scm_i_inexact_euclidean_divide
1508 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1509 else if (SCM_FRACTIONP (y
))
1510 return scm_i_slow_exact_euclidean_divide (x
, y
);
1512 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1513 s_scm_euclidean_divide
);
1515 else if (SCM_REALP (x
))
1517 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1518 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1519 return scm_i_inexact_euclidean_divide
1520 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1522 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1523 s_scm_euclidean_divide
);
1525 else if (SCM_FRACTIONP (x
))
1528 return scm_i_inexact_euclidean_divide
1529 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1531 return scm_i_slow_exact_euclidean_divide (x
, y
);
1534 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG1
,
1535 s_scm_euclidean_divide
);
1540 scm_i_inexact_euclidean_divide (double x
, double y
)
1544 if (SCM_LIKELY (y
> 0))
1546 else if (SCM_LIKELY (y
< 0))
1549 scm_num_overflow (s_scm_euclidean_divide
); /* or return a NaN? */
1553 return scm_values (scm_list_2 (scm_from_double (q
),
1554 scm_from_double (r
)));
1557 /* Compute exact euclidean quotient and remainder the slow way.
1558 We use this only if both arguments are exact,
1559 and at least one of them is a fraction */
1561 scm_i_slow_exact_euclidean_divide (SCM x
, SCM y
)
1565 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1566 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG1
,
1567 s_scm_euclidean_divide
);
1568 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1569 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1570 s_scm_euclidean_divide
);
1571 else if (scm_is_true (scm_positive_p (y
)))
1572 q
= scm_floor (scm_divide (x
, y
));
1573 else if (scm_is_true (scm_negative_p (y
)))
1574 q
= scm_ceiling (scm_divide (x
, y
));
1576 scm_num_overflow (s_scm_euclidean_divide
);
1577 r
= scm_difference (x
, scm_product (q
, y
));
1578 return scm_values (scm_list_2 (q
, r
));
1581 static SCM
scm_i_inexact_centered_quotient (double x
, double y
);
1582 static SCM
scm_i_bigint_centered_quotient (SCM x
, SCM y
);
1583 static SCM
scm_i_slow_exact_centered_quotient (SCM x
, SCM y
);
1585 SCM_PRIMITIVE_GENERIC (scm_centered_quotient
, "centered-quotient", 2, 0, 0,
1587 "Return the integer @var{q} such that\n"
1588 "@math{@var{x} = @var{q}*@var{y} + @var{r}} where\n"
1589 "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
1591 "(centered-quotient 123 10) @result{} 12\n"
1592 "(centered-quotient 123 -10) @result{} -12\n"
1593 "(centered-quotient -123 10) @result{} -12\n"
1594 "(centered-quotient -123 -10) @result{} 12\n"
1595 "(centered-quotient -123.2 -63.5) @result{} 2.0\n"
1596 "(centered-quotient 16/3 -10/7) @result{} -4\n"
1598 #define FUNC_NAME s_scm_centered_quotient
1600 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1602 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1604 scm_t_inum yy
= SCM_I_INUM (y
);
1605 if (SCM_UNLIKELY (yy
== 0))
1606 scm_num_overflow (s_scm_centered_quotient
);
1609 scm_t_inum xx
= SCM_I_INUM (x
);
1610 scm_t_inum qq
= xx
/ yy
;
1611 scm_t_inum rr
= xx
- qq
* yy
;
1612 if (SCM_LIKELY (xx
> 0))
1614 if (SCM_LIKELY (yy
> 0))
1616 if (rr
>= (yy
+ 1) / 2)
1621 if (rr
>= (1 - yy
) / 2)
1627 if (SCM_LIKELY (yy
> 0))
1638 return SCM_I_MAKINUM (qq
);
1641 else if (SCM_BIGP (y
))
1643 /* Pass a denormalized bignum version of x (even though it
1644 can fit in a fixnum) to scm_i_bigint_centered_quotient */
1645 return scm_i_bigint_centered_quotient
1646 (scm_i_long2big (SCM_I_INUM (x
)), y
);
1648 else if (SCM_REALP (y
))
1649 return scm_i_inexact_centered_quotient
1650 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1651 else if (SCM_FRACTIONP (y
))
1652 return scm_i_slow_exact_centered_quotient (x
, y
);
1654 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1655 s_scm_centered_quotient
);
1657 else if (SCM_BIGP (x
))
1659 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1661 scm_t_inum yy
= SCM_I_INUM (y
);
1662 if (SCM_UNLIKELY (yy
== 0))
1663 scm_num_overflow (s_scm_centered_quotient
);
1666 SCM q
= scm_i_mkbig ();
1668 /* Arrange for rr to initially be non-positive,
1669 because that simplifies the test to see
1670 if it is within the needed bounds. */
1673 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
1674 SCM_I_BIG_MPZ (x
), yy
);
1675 scm_remember_upto_here_1 (x
);
1677 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
1678 SCM_I_BIG_MPZ (q
), 1);
1682 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
1683 SCM_I_BIG_MPZ (x
), -yy
);
1684 scm_remember_upto_here_1 (x
);
1685 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1687 mpz_add_ui (SCM_I_BIG_MPZ (q
),
1688 SCM_I_BIG_MPZ (q
), 1);
1690 return scm_i_normbig (q
);
1693 else if (SCM_BIGP (y
))
1694 return scm_i_bigint_centered_quotient (x
, y
);
1695 else if (SCM_REALP (y
))
1696 return scm_i_inexact_centered_quotient
1697 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1698 else if (SCM_FRACTIONP (y
))
1699 return scm_i_slow_exact_centered_quotient (x
, y
);
1701 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1702 s_scm_centered_quotient
);
1704 else if (SCM_REALP (x
))
1706 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1707 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1708 return scm_i_inexact_centered_quotient
1709 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1711 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1712 s_scm_centered_quotient
);
1714 else if (SCM_FRACTIONP (x
))
1717 return scm_i_inexact_centered_quotient
1718 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1720 return scm_i_slow_exact_centered_quotient (x
, y
);
1723 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG1
,
1724 s_scm_centered_quotient
);
1729 scm_i_inexact_centered_quotient (double x
, double y
)
1731 if (SCM_LIKELY (y
> 0))
1732 return scm_from_double (floor (x
/y
+ 0.5));
1733 else if (SCM_LIKELY (y
< 0))
1734 return scm_from_double (ceil (x
/y
- 0.5));
1736 scm_num_overflow (s_scm_centered_quotient
); /* or return a NaN? */
1741 /* Assumes that both x and y are bigints, though
1742 x might be able to fit into a fixnum. */
1744 scm_i_bigint_centered_quotient (SCM x
, SCM y
)
1748 /* Note that x might be small enough to fit into a
1749 fixnum, so we must not let it escape into the wild */
1753 /* min_r will eventually become -abs(y)/2 */
1754 min_r
= scm_i_mkbig ();
1755 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
1756 SCM_I_BIG_MPZ (y
), 1);
1758 /* Arrange for rr to initially be non-positive,
1759 because that simplifies the test to see
1760 if it is within the needed bounds. */
1761 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1763 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1764 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1765 scm_remember_upto_here_2 (x
, y
);
1766 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
1767 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1768 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
1769 SCM_I_BIG_MPZ (q
), 1);
1773 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1774 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1775 scm_remember_upto_here_2 (x
, y
);
1776 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1777 mpz_add_ui (SCM_I_BIG_MPZ (q
),
1778 SCM_I_BIG_MPZ (q
), 1);
1780 scm_remember_upto_here_2 (r
, min_r
);
1781 return scm_i_normbig (q
);
1784 /* Compute exact centered quotient the slow way.
1785 We use this only if both arguments are exact,
1786 and at least one of them is a fraction */
1788 scm_i_slow_exact_centered_quotient (SCM x
, SCM y
)
1790 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1791 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG1
,
1792 s_scm_centered_quotient
);
1793 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1794 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1795 s_scm_centered_quotient
);
1796 else if (scm_is_true (scm_positive_p (y
)))
1797 return scm_floor (scm_sum (scm_divide (x
, y
),
1799 else if (scm_is_true (scm_negative_p (y
)))
1800 return scm_ceiling (scm_difference (scm_divide (x
, y
),
1803 scm_num_overflow (s_scm_centered_quotient
);
1806 static SCM
scm_i_inexact_centered_remainder (double x
, double y
);
1807 static SCM
scm_i_bigint_centered_remainder (SCM x
, SCM y
);
1808 static SCM
scm_i_slow_exact_centered_remainder (SCM x
, SCM y
);
1810 SCM_PRIMITIVE_GENERIC (scm_centered_remainder
, "centered-remainder", 2, 0, 0,
1812 "Return the real number @var{r} such that\n"
1813 "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}\n"
1814 "and @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1815 "for some integer @var{q}.\n"
1817 "(centered-remainder 123 10) @result{} 3\n"
1818 "(centered-remainder 123 -10) @result{} 3\n"
1819 "(centered-remainder -123 10) @result{} -3\n"
1820 "(centered-remainder -123 -10) @result{} -3\n"
1821 "(centered-remainder -123.2 -63.5) @result{} 3.8\n"
1822 "(centered-remainder 16/3 -10/7) @result{} -8/21\n"
1824 #define FUNC_NAME s_scm_centered_remainder
1826 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1828 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1830 scm_t_inum yy
= SCM_I_INUM (y
);
1831 if (SCM_UNLIKELY (yy
== 0))
1832 scm_num_overflow (s_scm_centered_remainder
);
1835 scm_t_inum xx
= SCM_I_INUM (x
);
1836 scm_t_inum rr
= xx
% yy
;
1837 if (SCM_LIKELY (xx
> 0))
1839 if (SCM_LIKELY (yy
> 0))
1841 if (rr
>= (yy
+ 1) / 2)
1846 if (rr
>= (1 - yy
) / 2)
1852 if (SCM_LIKELY (yy
> 0))
1863 return SCM_I_MAKINUM (rr
);
1866 else if (SCM_BIGP (y
))
1868 /* Pass a denormalized bignum version of x (even though it
1869 can fit in a fixnum) to scm_i_bigint_centered_remainder */
1870 return scm_i_bigint_centered_remainder
1871 (scm_i_long2big (SCM_I_INUM (x
)), y
);
1873 else if (SCM_REALP (y
))
1874 return scm_i_inexact_centered_remainder
1875 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1876 else if (SCM_FRACTIONP (y
))
1877 return scm_i_slow_exact_centered_remainder (x
, y
);
1879 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1880 s_scm_centered_remainder
);
1882 else if (SCM_BIGP (x
))
1884 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1886 scm_t_inum yy
= SCM_I_INUM (y
);
1887 if (SCM_UNLIKELY (yy
== 0))
1888 scm_num_overflow (s_scm_centered_remainder
);
1892 /* Arrange for rr to initially be non-positive,
1893 because that simplifies the test to see
1894 if it is within the needed bounds. */
1897 rr
= - mpz_cdiv_ui (SCM_I_BIG_MPZ (x
), yy
);
1898 scm_remember_upto_here_1 (x
);
1904 rr
= - mpz_cdiv_ui (SCM_I_BIG_MPZ (x
), -yy
);
1905 scm_remember_upto_here_1 (x
);
1909 return SCM_I_MAKINUM (rr
);
1912 else if (SCM_BIGP (y
))
1913 return scm_i_bigint_centered_remainder (x
, y
);
1914 else if (SCM_REALP (y
))
1915 return scm_i_inexact_centered_remainder
1916 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1917 else if (SCM_FRACTIONP (y
))
1918 return scm_i_slow_exact_centered_remainder (x
, y
);
1920 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1921 s_scm_centered_remainder
);
1923 else if (SCM_REALP (x
))
1925 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1926 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1927 return scm_i_inexact_centered_remainder
1928 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1930 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1931 s_scm_centered_remainder
);
1933 else if (SCM_FRACTIONP (x
))
1936 return scm_i_inexact_centered_remainder
1937 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1939 return scm_i_slow_exact_centered_remainder (x
, y
);
1942 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG1
,
1943 s_scm_centered_remainder
);
1948 scm_i_inexact_centered_remainder (double x
, double y
)
1952 /* Although it would be more efficient to use fmod here, we can't
1953 because it would in some cases produce results inconsistent with
1954 scm_i_inexact_centered_quotient, such that x != r + q * y (not even
1955 close). In particular, when x-y/2 is very close to a multiple of
1956 y, then r might be either -abs(y/2) or abs(y/2)-epsilon, but those
1957 two cases must correspond to different choices of q. If quotient
1958 chooses one and remainder chooses the other, it would be bad. */
1959 if (SCM_LIKELY (y
> 0))
1960 q
= floor (x
/y
+ 0.5);
1961 else if (SCM_LIKELY (y
< 0))
1962 q
= ceil (x
/y
- 0.5);
1964 scm_num_overflow (s_scm_centered_remainder
); /* or return a NaN? */
1967 return scm_from_double (x
- q
* y
);
1970 /* Assumes that both x and y are bigints, though
1971 x might be able to fit into a fixnum. */
1973 scm_i_bigint_centered_remainder (SCM x
, SCM y
)
1977 /* Note that x might be small enough to fit into a
1978 fixnum, so we must not let it escape into the wild */
1981 /* min_r will eventually become -abs(y)/2 */
1982 min_r
= scm_i_mkbig ();
1983 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
1984 SCM_I_BIG_MPZ (y
), 1);
1986 /* Arrange for rr to initially be non-positive,
1987 because that simplifies the test to see
1988 if it is within the needed bounds. */
1989 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1991 mpz_cdiv_r (SCM_I_BIG_MPZ (r
),
1992 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1993 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
1994 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1995 mpz_add (SCM_I_BIG_MPZ (r
),
2001 mpz_fdiv_r (SCM_I_BIG_MPZ (r
),
2002 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2003 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2004 mpz_sub (SCM_I_BIG_MPZ (r
),
2008 scm_remember_upto_here_2 (x
, y
);
2009 return scm_i_normbig (r
);
2012 /* Compute exact centered_remainder the slow way.
2013 We use this only if both arguments are exact,
2014 and at least one of them is a fraction */
2016 scm_i_slow_exact_centered_remainder (SCM x
, SCM y
)
2020 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
2021 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG1
,
2022 s_scm_centered_remainder
);
2023 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
2024 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
2025 s_scm_centered_remainder
);
2026 else if (scm_is_true (scm_positive_p (y
)))
2027 q
= scm_floor (scm_sum (scm_divide (x
, y
), exactly_one_half
));
2028 else if (scm_is_true (scm_negative_p (y
)))
2029 q
= scm_ceiling (scm_difference (scm_divide (x
, y
), exactly_one_half
));
2031 scm_num_overflow (s_scm_centered_remainder
);
2032 return scm_difference (x
, scm_product (y
, q
));
2036 static SCM
scm_i_inexact_centered_divide (double x
, double y
);
2037 static SCM
scm_i_bigint_centered_divide (SCM x
, SCM y
);
2038 static SCM
scm_i_slow_exact_centered_divide (SCM x
, SCM y
);
2040 SCM_PRIMITIVE_GENERIC (scm_centered_divide
, "centered/", 2, 0, 0,
2042 "Return the integer @var{q} and the real number @var{r}\n"
2043 "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
2044 "and @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
2046 "(centered/ 123 10) @result{} 12 and 3\n"
2047 "(centered/ 123 -10) @result{} -12 and 3\n"
2048 "(centered/ -123 10) @result{} -12 and -3\n"
2049 "(centered/ -123 -10) @result{} 12 and -3\n"
2050 "(centered/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
2051 "(centered/ 16/3 -10/7) @result{} -4 and -8/21\n"
2053 #define FUNC_NAME s_scm_centered_divide
2055 if (SCM_LIKELY (SCM_I_INUMP (x
)))
2057 if (SCM_LIKELY (SCM_I_INUMP (y
)))
2059 scm_t_inum yy
= SCM_I_INUM (y
);
2060 if (SCM_UNLIKELY (yy
== 0))
2061 scm_num_overflow (s_scm_centered_divide
);
2064 scm_t_inum xx
= SCM_I_INUM (x
);
2065 scm_t_inum qq
= xx
/ yy
;
2066 scm_t_inum rr
= xx
- qq
* yy
;
2067 if (SCM_LIKELY (xx
> 0))
2069 if (SCM_LIKELY (yy
> 0))
2071 if (rr
>= (yy
+ 1) / 2)
2076 if (rr
>= (1 - yy
) / 2)
2082 if (SCM_LIKELY (yy
> 0))
2093 return scm_values (scm_list_2 (SCM_I_MAKINUM (qq
),
2094 SCM_I_MAKINUM (rr
)));
2097 else if (SCM_BIGP (y
))
2099 /* Pass a denormalized bignum version of x (even though it
2100 can fit in a fixnum) to scm_i_bigint_centered_divide */
2101 return scm_i_bigint_centered_divide
2102 (scm_i_long2big (SCM_I_INUM (x
)), y
);
2104 else if (SCM_REALP (y
))
2105 return scm_i_inexact_centered_divide
2106 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
2107 else if (SCM_FRACTIONP (y
))
2108 return scm_i_slow_exact_centered_divide (x
, y
);
2110 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2111 s_scm_centered_divide
);
2113 else if (SCM_BIGP (x
))
2115 if (SCM_LIKELY (SCM_I_INUMP (y
)))
2117 scm_t_inum yy
= SCM_I_INUM (y
);
2118 if (SCM_UNLIKELY (yy
== 0))
2119 scm_num_overflow (s_scm_centered_divide
);
2122 SCM q
= scm_i_mkbig ();
2124 /* Arrange for rr to initially be non-positive,
2125 because that simplifies the test to see
2126 if it is within the needed bounds. */
2129 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
2130 SCM_I_BIG_MPZ (x
), yy
);
2131 scm_remember_upto_here_1 (x
);
2134 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
2135 SCM_I_BIG_MPZ (q
), 1);
2141 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
2142 SCM_I_BIG_MPZ (x
), -yy
);
2143 scm_remember_upto_here_1 (x
);
2144 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
2147 mpz_add_ui (SCM_I_BIG_MPZ (q
),
2148 SCM_I_BIG_MPZ (q
), 1);
2152 return scm_values (scm_list_2 (scm_i_normbig (q
),
2153 SCM_I_MAKINUM (rr
)));
2156 else if (SCM_BIGP (y
))
2157 return scm_i_bigint_centered_divide (x
, y
);
2158 else if (SCM_REALP (y
))
2159 return scm_i_inexact_centered_divide
2160 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
2161 else if (SCM_FRACTIONP (y
))
2162 return scm_i_slow_exact_centered_divide (x
, y
);
2164 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2165 s_scm_centered_divide
);
2167 else if (SCM_REALP (x
))
2169 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
2170 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
2171 return scm_i_inexact_centered_divide
2172 (SCM_REAL_VALUE (x
), scm_to_double (y
));
2174 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2175 s_scm_centered_divide
);
2177 else if (SCM_FRACTIONP (x
))
2180 return scm_i_inexact_centered_divide
2181 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
2183 return scm_i_slow_exact_centered_divide (x
, y
);
2186 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG1
,
2187 s_scm_centered_divide
);
2192 scm_i_inexact_centered_divide (double x
, double y
)
2196 if (SCM_LIKELY (y
> 0))
2197 q
= floor (x
/y
+ 0.5);
2198 else if (SCM_LIKELY (y
< 0))
2199 q
= ceil (x
/y
- 0.5);
2201 scm_num_overflow (s_scm_centered_divide
); /* or return a NaN? */
2205 return scm_values (scm_list_2 (scm_from_double (q
),
2206 scm_from_double (r
)));
2209 /* Assumes that both x and y are bigints, though
2210 x might be able to fit into a fixnum. */
2212 scm_i_bigint_centered_divide (SCM x
, SCM y
)
2216 /* Note that x might be small enough to fit into a
2217 fixnum, so we must not let it escape into the wild */
2221 /* min_r will eventually become -abs(y/2) */
2222 min_r
= scm_i_mkbig ();
2223 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
2224 SCM_I_BIG_MPZ (y
), 1);
2226 /* Arrange for rr to initially be non-positive,
2227 because that simplifies the test to see
2228 if it is within the needed bounds. */
2229 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
2231 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
2232 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2233 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
2234 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2236 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
2237 SCM_I_BIG_MPZ (q
), 1);
2238 mpz_add (SCM_I_BIG_MPZ (r
),
2245 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
2246 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2247 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2249 mpz_add_ui (SCM_I_BIG_MPZ (q
),
2250 SCM_I_BIG_MPZ (q
), 1);
2251 mpz_sub (SCM_I_BIG_MPZ (r
),
2256 scm_remember_upto_here_2 (x
, y
);
2257 return scm_values (scm_list_2 (scm_i_normbig (q
),
2258 scm_i_normbig (r
)));
2261 /* Compute exact centered quotient and remainder the slow way.
2262 We use this only if both arguments are exact,
2263 and at least one of them is a fraction */
2265 scm_i_slow_exact_centered_divide (SCM x
, SCM y
)
2269 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
2270 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG1
,
2271 s_scm_centered_divide
);
2272 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
2273 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2274 s_scm_centered_divide
);
2275 else if (scm_is_true (scm_positive_p (y
)))
2276 q
= scm_floor (scm_sum (scm_divide (x
, y
),
2278 else if (scm_is_true (scm_negative_p (y
)))
2279 q
= scm_ceiling (scm_difference (scm_divide (x
, y
),
2282 scm_num_overflow (s_scm_centered_divide
);
2283 r
= scm_difference (x
, scm_product (q
, y
));
2284 return scm_values (scm_list_2 (q
, r
));
2288 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
2289 (SCM x
, SCM y
, SCM rest
),
2290 "Return the greatest common divisor of all parameter values.\n"
2291 "If called without arguments, 0 is returned.")
2292 #define FUNC_NAME s_scm_i_gcd
2294 while (!scm_is_null (rest
))
2295 { x
= scm_gcd (x
, y
);
2297 rest
= scm_cdr (rest
);
2299 return scm_gcd (x
, y
);
2303 #define s_gcd s_scm_i_gcd
2304 #define g_gcd g_scm_i_gcd
2307 scm_gcd (SCM x
, SCM y
)
2310 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
2312 if (SCM_I_INUMP (x
))
2314 if (SCM_I_INUMP (y
))
2316 scm_t_inum xx
= SCM_I_INUM (x
);
2317 scm_t_inum yy
= SCM_I_INUM (y
);
2318 scm_t_inum u
= xx
< 0 ? -xx
: xx
;
2319 scm_t_inum v
= yy
< 0 ? -yy
: yy
;
2329 /* Determine a common factor 2^k */
2330 while (!(1 & (u
| v
)))
2336 /* Now, any factor 2^n can be eliminated */
2356 return (SCM_POSFIXABLE (result
)
2357 ? SCM_I_MAKINUM (result
)
2358 : scm_i_inum2big (result
));
2360 else if (SCM_BIGP (y
))
2366 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
2368 else if (SCM_BIGP (x
))
2370 if (SCM_I_INUMP (y
))
2375 yy
= SCM_I_INUM (y
);
2380 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
2381 scm_remember_upto_here_1 (x
);
2382 return (SCM_POSFIXABLE (result
)
2383 ? SCM_I_MAKINUM (result
)
2384 : scm_from_unsigned_integer (result
));
2386 else if (SCM_BIGP (y
))
2388 SCM result
= scm_i_mkbig ();
2389 mpz_gcd (SCM_I_BIG_MPZ (result
),
2392 scm_remember_upto_here_2 (x
, y
);
2393 return scm_i_normbig (result
);
2396 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
2399 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
2402 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
2403 (SCM x
, SCM y
, SCM rest
),
2404 "Return the least common multiple of the arguments.\n"
2405 "If called without arguments, 1 is returned.")
2406 #define FUNC_NAME s_scm_i_lcm
2408 while (!scm_is_null (rest
))
2409 { x
= scm_lcm (x
, y
);
2411 rest
= scm_cdr (rest
);
2413 return scm_lcm (x
, y
);
2417 #define s_lcm s_scm_i_lcm
2418 #define g_lcm g_scm_i_lcm
2421 scm_lcm (SCM n1
, SCM n2
)
2423 if (SCM_UNBNDP (n2
))
2425 if (SCM_UNBNDP (n1
))
2426 return SCM_I_MAKINUM (1L);
2427 n2
= SCM_I_MAKINUM (1L);
2430 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
2431 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
2432 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
2433 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
2435 if (SCM_I_INUMP (n1
))
2437 if (SCM_I_INUMP (n2
))
2439 SCM d
= scm_gcd (n1
, n2
);
2440 if (scm_is_eq (d
, SCM_INUM0
))
2443 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
2447 /* inum n1, big n2 */
2450 SCM result
= scm_i_mkbig ();
2451 scm_t_inum nn1
= SCM_I_INUM (n1
);
2452 if (nn1
== 0) return SCM_INUM0
;
2453 if (nn1
< 0) nn1
= - nn1
;
2454 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
2455 scm_remember_upto_here_1 (n2
);
2463 if (SCM_I_INUMP (n2
))
2470 SCM result
= scm_i_mkbig ();
2471 mpz_lcm(SCM_I_BIG_MPZ (result
),
2473 SCM_I_BIG_MPZ (n2
));
2474 scm_remember_upto_here_2(n1
, n2
);
2475 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
2481 /* Emulating 2's complement bignums with sign magnitude arithmetic:
2486 + + + x (map digit:logand X Y)
2487 + - + x (map digit:logand X (lognot (+ -1 Y)))
2488 - + + y (map digit:logand (lognot (+ -1 X)) Y)
2489 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
2494 + + + (map digit:logior X Y)
2495 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
2496 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
2497 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
2502 + + + (map digit:logxor X Y)
2503 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
2504 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
2505 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
2510 + + (any digit:logand X Y)
2511 + - (any digit:logand X (lognot (+ -1 Y)))
2512 - + (any digit:logand (lognot (+ -1 X)) Y)
2517 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
2518 (SCM x
, SCM y
, SCM rest
),
2519 "Return the bitwise AND of the integer arguments.\n\n"
2521 "(logand) @result{} -1\n"
2522 "(logand 7) @result{} 7\n"
2523 "(logand #b111 #b011 #b001) @result{} 1\n"
2525 #define FUNC_NAME s_scm_i_logand
2527 while (!scm_is_null (rest
))
2528 { x
= scm_logand (x
, y
);
2530 rest
= scm_cdr (rest
);
2532 return scm_logand (x
, y
);
2536 #define s_scm_logand s_scm_i_logand
2538 SCM
scm_logand (SCM n1
, SCM n2
)
2539 #define FUNC_NAME s_scm_logand
2543 if (SCM_UNBNDP (n2
))
2545 if (SCM_UNBNDP (n1
))
2546 return SCM_I_MAKINUM (-1);
2547 else if (!SCM_NUMBERP (n1
))
2548 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2549 else if (SCM_NUMBERP (n1
))
2552 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2555 if (SCM_I_INUMP (n1
))
2557 nn1
= SCM_I_INUM (n1
);
2558 if (SCM_I_INUMP (n2
))
2560 scm_t_inum nn2
= SCM_I_INUM (n2
);
2561 return SCM_I_MAKINUM (nn1
& nn2
);
2563 else if SCM_BIGP (n2
)
2569 SCM result_z
= scm_i_mkbig ();
2571 mpz_init_set_si (nn1_z
, nn1
);
2572 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2573 scm_remember_upto_here_1 (n2
);
2575 return scm_i_normbig (result_z
);
2579 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2581 else if (SCM_BIGP (n1
))
2583 if (SCM_I_INUMP (n2
))
2586 nn1
= SCM_I_INUM (n1
);
2589 else if (SCM_BIGP (n2
))
2591 SCM result_z
= scm_i_mkbig ();
2592 mpz_and (SCM_I_BIG_MPZ (result_z
),
2594 SCM_I_BIG_MPZ (n2
));
2595 scm_remember_upto_here_2 (n1
, n2
);
2596 return scm_i_normbig (result_z
);
2599 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2602 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2607 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
2608 (SCM x
, SCM y
, SCM rest
),
2609 "Return the bitwise OR of the integer arguments.\n\n"
2611 "(logior) @result{} 0\n"
2612 "(logior 7) @result{} 7\n"
2613 "(logior #b000 #b001 #b011) @result{} 3\n"
2615 #define FUNC_NAME s_scm_i_logior
2617 while (!scm_is_null (rest
))
2618 { x
= scm_logior (x
, y
);
2620 rest
= scm_cdr (rest
);
2622 return scm_logior (x
, y
);
2626 #define s_scm_logior s_scm_i_logior
2628 SCM
scm_logior (SCM n1
, SCM n2
)
2629 #define FUNC_NAME s_scm_logior
2633 if (SCM_UNBNDP (n2
))
2635 if (SCM_UNBNDP (n1
))
2637 else if (SCM_NUMBERP (n1
))
2640 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2643 if (SCM_I_INUMP (n1
))
2645 nn1
= SCM_I_INUM (n1
);
2646 if (SCM_I_INUMP (n2
))
2648 long nn2
= SCM_I_INUM (n2
);
2649 return SCM_I_MAKINUM (nn1
| nn2
);
2651 else if (SCM_BIGP (n2
))
2657 SCM result_z
= scm_i_mkbig ();
2659 mpz_init_set_si (nn1_z
, nn1
);
2660 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2661 scm_remember_upto_here_1 (n2
);
2663 return scm_i_normbig (result_z
);
2667 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2669 else if (SCM_BIGP (n1
))
2671 if (SCM_I_INUMP (n2
))
2674 nn1
= SCM_I_INUM (n1
);
2677 else if (SCM_BIGP (n2
))
2679 SCM result_z
= scm_i_mkbig ();
2680 mpz_ior (SCM_I_BIG_MPZ (result_z
),
2682 SCM_I_BIG_MPZ (n2
));
2683 scm_remember_upto_here_2 (n1
, n2
);
2684 return scm_i_normbig (result_z
);
2687 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2690 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2695 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
2696 (SCM x
, SCM y
, SCM rest
),
2697 "Return the bitwise XOR of the integer arguments. A bit is\n"
2698 "set in the result if it is set in an odd number of arguments.\n"
2700 "(logxor) @result{} 0\n"
2701 "(logxor 7) @result{} 7\n"
2702 "(logxor #b000 #b001 #b011) @result{} 2\n"
2703 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
2705 #define FUNC_NAME s_scm_i_logxor
2707 while (!scm_is_null (rest
))
2708 { x
= scm_logxor (x
, y
);
2710 rest
= scm_cdr (rest
);
2712 return scm_logxor (x
, y
);
2716 #define s_scm_logxor s_scm_i_logxor
2718 SCM
scm_logxor (SCM n1
, SCM n2
)
2719 #define FUNC_NAME s_scm_logxor
2723 if (SCM_UNBNDP (n2
))
2725 if (SCM_UNBNDP (n1
))
2727 else if (SCM_NUMBERP (n1
))
2730 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2733 if (SCM_I_INUMP (n1
))
2735 nn1
= SCM_I_INUM (n1
);
2736 if (SCM_I_INUMP (n2
))
2738 scm_t_inum nn2
= SCM_I_INUM (n2
);
2739 return SCM_I_MAKINUM (nn1
^ nn2
);
2741 else if (SCM_BIGP (n2
))
2745 SCM result_z
= scm_i_mkbig ();
2747 mpz_init_set_si (nn1_z
, nn1
);
2748 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2749 scm_remember_upto_here_1 (n2
);
2751 return scm_i_normbig (result_z
);
2755 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2757 else if (SCM_BIGP (n1
))
2759 if (SCM_I_INUMP (n2
))
2762 nn1
= SCM_I_INUM (n1
);
2765 else if (SCM_BIGP (n2
))
2767 SCM result_z
= scm_i_mkbig ();
2768 mpz_xor (SCM_I_BIG_MPZ (result_z
),
2770 SCM_I_BIG_MPZ (n2
));
2771 scm_remember_upto_here_2 (n1
, n2
);
2772 return scm_i_normbig (result_z
);
2775 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2778 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2783 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
2785 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
2786 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
2787 "without actually calculating the @code{logand}, just testing\n"
2791 "(logtest #b0100 #b1011) @result{} #f\n"
2792 "(logtest #b0100 #b0111) @result{} #t\n"
2794 #define FUNC_NAME s_scm_logtest
2798 if (SCM_I_INUMP (j
))
2800 nj
= SCM_I_INUM (j
);
2801 if (SCM_I_INUMP (k
))
2803 scm_t_inum nk
= SCM_I_INUM (k
);
2804 return scm_from_bool (nj
& nk
);
2806 else if (SCM_BIGP (k
))
2814 mpz_init_set_si (nj_z
, nj
);
2815 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
2816 scm_remember_upto_here_1 (k
);
2817 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
2823 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
2825 else if (SCM_BIGP (j
))
2827 if (SCM_I_INUMP (k
))
2830 nj
= SCM_I_INUM (j
);
2833 else if (SCM_BIGP (k
))
2837 mpz_init (result_z
);
2841 scm_remember_upto_here_2 (j
, k
);
2842 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
2843 mpz_clear (result_z
);
2847 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
2850 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
2855 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
2857 "Test whether bit number @var{index} in @var{j} is set.\n"
2858 "@var{index} starts from 0 for the least significant bit.\n"
2861 "(logbit? 0 #b1101) @result{} #t\n"
2862 "(logbit? 1 #b1101) @result{} #f\n"
2863 "(logbit? 2 #b1101) @result{} #t\n"
2864 "(logbit? 3 #b1101) @result{} #t\n"
2865 "(logbit? 4 #b1101) @result{} #f\n"
2867 #define FUNC_NAME s_scm_logbit_p
2869 unsigned long int iindex
;
2870 iindex
= scm_to_ulong (index
);
2872 if (SCM_I_INUMP (j
))
2874 /* bits above what's in an inum follow the sign bit */
2875 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
2876 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
2878 else if (SCM_BIGP (j
))
2880 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
2881 scm_remember_upto_here_1 (j
);
2882 return scm_from_bool (val
);
2885 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
2890 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
2892 "Return the integer which is the ones-complement of the integer\n"
2896 "(number->string (lognot #b10000000) 2)\n"
2897 " @result{} \"-10000001\"\n"
2898 "(number->string (lognot #b0) 2)\n"
2899 " @result{} \"-1\"\n"
2901 #define FUNC_NAME s_scm_lognot
2903 if (SCM_I_INUMP (n
)) {
2904 /* No overflow here, just need to toggle all the bits making up the inum.
2905 Enhancement: No need to strip the tag and add it back, could just xor
2906 a block of 1 bits, if that worked with the various debug versions of
2908 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
2910 } else if (SCM_BIGP (n
)) {
2911 SCM result
= scm_i_mkbig ();
2912 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
2913 scm_remember_upto_here_1 (n
);
2917 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2922 /* returns 0 if IN is not an integer. OUT must already be
2925 coerce_to_big (SCM in
, mpz_t out
)
2928 mpz_set (out
, SCM_I_BIG_MPZ (in
));
2929 else if (SCM_I_INUMP (in
))
2930 mpz_set_si (out
, SCM_I_INUM (in
));
2937 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
2938 (SCM n
, SCM k
, SCM m
),
2939 "Return @var{n} raised to the integer exponent\n"
2940 "@var{k}, modulo @var{m}.\n"
2943 "(modulo-expt 2 3 5)\n"
2946 #define FUNC_NAME s_scm_modulo_expt
2952 /* There are two classes of error we might encounter --
2953 1) Math errors, which we'll report by calling scm_num_overflow,
2955 2) wrong-type errors, which of course we'll report by calling
2957 We don't report those errors immediately, however; instead we do
2958 some cleanup first. These variables tell us which error (if
2959 any) we should report after cleaning up.
2961 int report_overflow
= 0;
2963 int position_of_wrong_type
= 0;
2964 SCM value_of_wrong_type
= SCM_INUM0
;
2966 SCM result
= SCM_UNDEFINED
;
2972 if (scm_is_eq (m
, SCM_INUM0
))
2974 report_overflow
= 1;
2978 if (!coerce_to_big (n
, n_tmp
))
2980 value_of_wrong_type
= n
;
2981 position_of_wrong_type
= 1;
2985 if (!coerce_to_big (k
, k_tmp
))
2987 value_of_wrong_type
= k
;
2988 position_of_wrong_type
= 2;
2992 if (!coerce_to_big (m
, m_tmp
))
2994 value_of_wrong_type
= m
;
2995 position_of_wrong_type
= 3;
2999 /* if the exponent K is negative, and we simply call mpz_powm, we
3000 will get a divide-by-zero exception when an inverse 1/n mod m
3001 doesn't exist (or is not unique). Since exceptions are hard to
3002 handle, we'll attempt the inversion "by hand" -- that way, we get
3003 a simple failure code, which is easy to handle. */
3005 if (-1 == mpz_sgn (k_tmp
))
3007 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
3009 report_overflow
= 1;
3012 mpz_neg (k_tmp
, k_tmp
);
3015 result
= scm_i_mkbig ();
3016 mpz_powm (SCM_I_BIG_MPZ (result
),
3021 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
3022 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
3029 if (report_overflow
)
3030 scm_num_overflow (FUNC_NAME
);
3032 if (position_of_wrong_type
)
3033 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
3034 value_of_wrong_type
);
3036 return scm_i_normbig (result
);
3040 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
3042 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
3043 "exact integer, @var{n} can be any number.\n"
3045 "Negative @var{k} is supported, and results in\n"
3046 "@math{1/@var{n}^abs(@var{k})} in the usual way.\n"
3047 "@math{@var{n}^0} is 1, as usual, and that\n"
3048 "includes @math{0^0} is 1.\n"
3051 "(integer-expt 2 5) @result{} 32\n"
3052 "(integer-expt -3 3) @result{} -27\n"
3053 "(integer-expt 5 -3) @result{} 1/125\n"
3054 "(integer-expt 0 0) @result{} 1\n"
3056 #define FUNC_NAME s_scm_integer_expt
3059 SCM z_i2
= SCM_BOOL_F
;
3061 SCM acc
= SCM_I_MAKINUM (1L);
3063 /* Specifically refrain from checking the type of the first argument.
3064 This allows us to exponentiate any object that can be multiplied.
3065 If we must raise to a negative power, we must also be able to
3066 take its reciprocal. */
3067 if (!SCM_LIKELY (SCM_I_INUMP (k
)) && !SCM_LIKELY (SCM_BIGP (k
)))
3068 SCM_WRONG_TYPE_ARG (2, k
);
3070 if (SCM_UNLIKELY (scm_is_eq (k
, SCM_INUM0
)))
3071 return SCM_INUM1
; /* n^(exact0) is exact 1, regardless of n */
3072 else if (SCM_UNLIKELY (scm_is_eq (n
, SCM_I_MAKINUM (-1L))))
3073 return scm_is_false (scm_even_p (k
)) ? n
: SCM_INUM1
;
3074 /* The next check is necessary only because R6RS specifies different
3075 behavior for 0^(-k) than for (/ 0). If n is not a scheme number,
3076 we simply skip this case and move on. */
3077 else if (SCM_NUMBERP (n
) && scm_is_true (scm_zero_p (n
)))
3079 /* k cannot be 0 at this point, because we
3080 have already checked for that case above */
3081 if (scm_is_true (scm_positive_p (k
)))
3083 else /* return NaN for (0 ^ k) for negative k per R6RS */
3087 if (SCM_I_INUMP (k
))
3088 i2
= SCM_I_INUM (k
);
3089 else if (SCM_BIGP (k
))
3091 z_i2
= scm_i_clonebig (k
, 1);
3092 scm_remember_upto_here_1 (k
);
3096 SCM_WRONG_TYPE_ARG (2, k
);
3100 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
3102 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
3103 n
= scm_divide (n
, SCM_UNDEFINED
);
3107 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
3111 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
3113 return scm_product (acc
, n
);
3115 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
3116 acc
= scm_product (acc
, n
);
3117 n
= scm_product (n
, n
);
3118 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
3126 n
= scm_divide (n
, SCM_UNDEFINED
);
3133 return scm_product (acc
, n
);
3135 acc
= scm_product (acc
, n
);
3136 n
= scm_product (n
, n
);
3143 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
3145 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
3146 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
3148 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
3149 "@var{cnt} is negative it's a division, rounded towards negative\n"
3150 "infinity. (Note that this is not the same rounding as\n"
3151 "@code{quotient} does.)\n"
3153 "With @var{n} viewed as an infinite precision twos complement,\n"
3154 "@code{ash} means a left shift introducing zero bits, or a right\n"
3155 "shift dropping bits.\n"
3158 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
3159 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
3161 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
3162 "(ash -23 -2) @result{} -6\n"
3164 #define FUNC_NAME s_scm_ash
3167 bits_to_shift
= scm_to_long (cnt
);
3169 if (SCM_I_INUMP (n
))
3171 scm_t_inum nn
= SCM_I_INUM (n
);
3173 if (bits_to_shift
> 0)
3175 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
3176 overflow a non-zero fixnum. For smaller shifts we check the
3177 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
3178 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
3179 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
3185 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
3187 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
3190 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
3194 SCM result
= scm_i_inum2big (nn
);
3195 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
3202 bits_to_shift
= -bits_to_shift
;
3203 if (bits_to_shift
>= SCM_LONG_BIT
)
3204 return (nn
>= 0 ? SCM_INUM0
: SCM_I_MAKINUM(-1));
3206 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
3210 else if (SCM_BIGP (n
))
3214 if (bits_to_shift
== 0)
3217 result
= scm_i_mkbig ();
3218 if (bits_to_shift
>= 0)
3220 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
3226 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
3227 we have to allocate a bignum even if the result is going to be a
3229 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
3231 return scm_i_normbig (result
);
3237 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3243 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
3244 (SCM n
, SCM start
, SCM end
),
3245 "Return the integer composed of the @var{start} (inclusive)\n"
3246 "through @var{end} (exclusive) bits of @var{n}. The\n"
3247 "@var{start}th bit becomes the 0-th bit in the result.\n"
3250 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
3251 " @result{} \"1010\"\n"
3252 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
3253 " @result{} \"10110\"\n"
3255 #define FUNC_NAME s_scm_bit_extract
3257 unsigned long int istart
, iend
, bits
;
3258 istart
= scm_to_ulong (start
);
3259 iend
= scm_to_ulong (end
);
3260 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
3262 /* how many bits to keep */
3263 bits
= iend
- istart
;
3265 if (SCM_I_INUMP (n
))
3267 scm_t_inum in
= SCM_I_INUM (n
);
3269 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
3270 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
3271 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
3273 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
3275 /* Since we emulate two's complement encoded numbers, this
3276 * special case requires us to produce a result that has
3277 * more bits than can be stored in a fixnum.
3279 SCM result
= scm_i_inum2big (in
);
3280 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
3285 /* mask down to requisite bits */
3286 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
3287 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
3289 else if (SCM_BIGP (n
))
3294 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
3298 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
3299 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
3300 such bits into a ulong. */
3301 result
= scm_i_mkbig ();
3302 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
3303 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
3304 result
= scm_i_normbig (result
);
3306 scm_remember_upto_here_1 (n
);
3310 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3315 static const char scm_logtab
[] = {
3316 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
3319 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
3321 "Return the number of bits in integer @var{n}. If integer is\n"
3322 "positive, the 1-bits in its binary representation are counted.\n"
3323 "If negative, the 0-bits in its two's-complement binary\n"
3324 "representation are counted. If 0, 0 is returned.\n"
3327 "(logcount #b10101010)\n"
3334 #define FUNC_NAME s_scm_logcount
3336 if (SCM_I_INUMP (n
))
3338 unsigned long c
= 0;
3339 scm_t_inum nn
= SCM_I_INUM (n
);
3344 c
+= scm_logtab
[15 & nn
];
3347 return SCM_I_MAKINUM (c
);
3349 else if (SCM_BIGP (n
))
3351 unsigned long count
;
3352 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
3353 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
3355 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
3356 scm_remember_upto_here_1 (n
);
3357 return SCM_I_MAKINUM (count
);
3360 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3365 static const char scm_ilentab
[] = {
3366 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
3370 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
3372 "Return the number of bits necessary to represent @var{n}.\n"
3375 "(integer-length #b10101010)\n"
3377 "(integer-length 0)\n"
3379 "(integer-length #b1111)\n"
3382 #define FUNC_NAME s_scm_integer_length
3384 if (SCM_I_INUMP (n
))
3386 unsigned long c
= 0;
3388 scm_t_inum nn
= SCM_I_INUM (n
);
3394 l
= scm_ilentab
[15 & nn
];
3397 return SCM_I_MAKINUM (c
- 4 + l
);
3399 else if (SCM_BIGP (n
))
3401 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
3402 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
3403 1 too big, so check for that and adjust. */
3404 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
3405 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
3406 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
3407 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
3409 scm_remember_upto_here_1 (n
);
3410 return SCM_I_MAKINUM (size
);
3413 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3417 /*** NUMBERS -> STRINGS ***/
3418 #define SCM_MAX_DBL_PREC 60
3419 #define SCM_MAX_DBL_RADIX 36
3421 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
3422 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
3423 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
3426 void init_dblprec(int *prec
, int radix
) {
3427 /* determine floating point precision by adding successively
3428 smaller increments to 1.0 until it is considered == 1.0 */
3429 double f
= ((double)1.0)/radix
;
3430 double fsum
= 1.0 + f
;
3435 if (++(*prec
) > SCM_MAX_DBL_PREC
)
3447 void init_fx_radix(double *fx_list
, int radix
)
3449 /* initialize a per-radix list of tolerances. When added
3450 to a number < 1.0, we can determine if we should raund
3451 up and quit converting a number to a string. */
3455 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
3456 fx_list
[i
] = (fx_list
[i
-1] / radix
);
3459 /* use this array as a way to generate a single digit */
3460 static const char number_chars
[] = "0123456789abcdefghijklmnopqrstuvwxyz";
3463 idbl2str (double f
, char *a
, int radix
)
3465 int efmt
, dpt
, d
, i
, wp
;
3467 #ifdef DBL_MIN_10_EXP
3470 #endif /* DBL_MIN_10_EXP */
3475 radix
> SCM_MAX_DBL_RADIX
)
3477 /* revert to existing behavior */
3481 wp
= scm_dblprec
[radix
-2];
3482 fx
= fx_per_radix
[radix
-2];
3486 #ifdef HAVE_COPYSIGN
3487 double sgn
= copysign (1.0, f
);
3492 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
3498 strcpy (a
, "-inf.0");
3500 strcpy (a
, "+inf.0");
3505 strcpy (a
, "+nan.0");
3515 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
3516 make-uniform-vector, from causing infinite loops. */
3517 /* just do the checking...if it passes, we do the conversion for our
3518 radix again below */
3525 if (exp_cpy
-- < DBL_MIN_10_EXP
)
3533 while (f_cpy
> 10.0)
3536 if (exp_cpy
++ > DBL_MAX_10_EXP
)
3557 if (f
+ fx
[wp
] >= radix
)
3564 /* adding 9999 makes this equivalent to abs(x) % 3 */
3565 dpt
= (exp
+ 9999) % 3;
3569 efmt
= (exp
< -3) || (exp
> wp
+ 2);
3591 a
[ch
++] = number_chars
[d
];
3594 if (f
+ fx
[wp
] >= 1.0)
3596 a
[ch
- 1] = number_chars
[d
+1];
3608 if ((dpt
> 4) && (exp
> 6))
3610 d
= (a
[0] == '-' ? 2 : 1);
3611 for (i
= ch
++; i
> d
; i
--)
3624 if (a
[ch
- 1] == '.')
3625 a
[ch
++] = '0'; /* trailing zero */
3634 for (i
= radix
; i
<= exp
; i
*= radix
);
3635 for (i
/= radix
; i
; i
/= radix
)
3637 a
[ch
++] = number_chars
[exp
/ i
];
3646 icmplx2str (double real
, double imag
, char *str
, int radix
)
3650 i
= idbl2str (real
, str
, radix
);
3653 /* Don't output a '+' for negative numbers or for Inf and
3654 NaN. They will provide their own sign. */
3655 if (0 <= imag
&& !isinf (imag
) && !isnan (imag
))
3657 i
+= idbl2str (imag
, &str
[i
], radix
);
3664 iflo2str (SCM flt
, char *str
, int radix
)
3667 if (SCM_REALP (flt
))
3668 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
3670 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
3675 /* convert a scm_t_intmax to a string (unterminated). returns the number of
3676 characters in the result.
3678 p is destination: worst case (base 2) is SCM_INTBUFLEN */
3680 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
3685 return scm_iuint2str (-num
, rad
, p
) + 1;
3688 return scm_iuint2str (num
, rad
, p
);
3691 /* convert a scm_t_intmax to a string (unterminated). returns the number of
3692 characters in the result.
3694 p is destination: worst case (base 2) is SCM_INTBUFLEN */
3696 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
3700 scm_t_uintmax n
= num
;
3702 if (rad
< 2 || rad
> 36)
3703 scm_out_of_range ("scm_iuint2str", scm_from_int (rad
));
3705 for (n
/= rad
; n
> 0; n
/= rad
)
3715 p
[i
] = number_chars
[d
];
3720 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
3722 "Return a string holding the external representation of the\n"
3723 "number @var{n} in the given @var{radix}. If @var{n} is\n"
3724 "inexact, a radix of 10 will be used.")
3725 #define FUNC_NAME s_scm_number_to_string
3729 if (SCM_UNBNDP (radix
))
3732 base
= scm_to_signed_integer (radix
, 2, 36);
3734 if (SCM_I_INUMP (n
))
3736 char num_buf
[SCM_INTBUFLEN
];
3737 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
3738 return scm_from_locale_stringn (num_buf
, length
);
3740 else if (SCM_BIGP (n
))
3742 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
3743 scm_remember_upto_here_1 (n
);
3744 return scm_take_locale_string (str
);
3746 else if (SCM_FRACTIONP (n
))
3748 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
3749 scm_from_locale_string ("/"),
3750 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
3752 else if (SCM_INEXACTP (n
))
3754 char num_buf
[FLOBUFLEN
];
3755 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
3758 SCM_WRONG_TYPE_ARG (1, n
);
3763 /* These print routines used to be stubbed here so that scm_repl.c
3764 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
3767 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3769 char num_buf
[FLOBUFLEN
];
3770 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
3775 scm_i_print_double (double val
, SCM port
)
3777 char num_buf
[FLOBUFLEN
];
3778 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
3782 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3785 char num_buf
[FLOBUFLEN
];
3786 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
3791 scm_i_print_complex (double real
, double imag
, SCM port
)
3793 char num_buf
[FLOBUFLEN
];
3794 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
3798 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3801 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
3802 scm_display (str
, port
);
3803 scm_remember_upto_here_1 (str
);
3808 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3810 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
3811 scm_remember_upto_here_1 (exp
);
3812 scm_lfwrite (str
, (size_t) strlen (str
), port
);
3816 /*** END nums->strs ***/
3819 /*** STRINGS -> NUMBERS ***/
3821 /* The following functions implement the conversion from strings to numbers.
3822 * The implementation somehow follows the grammar for numbers as it is given
3823 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
3824 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
3825 * points should be noted about the implementation:
3826 * * Each function keeps a local index variable 'idx' that points at the
3827 * current position within the parsed string. The global index is only
3828 * updated if the function could parse the corresponding syntactic unit
3830 * * Similarly, the functions keep track of indicators of inexactness ('#',
3831 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
3832 * global exactness information is only updated after each part has been
3833 * successfully parsed.
3834 * * Sequences of digits are parsed into temporary variables holding fixnums.
3835 * Only if these fixnums would overflow, the result variables are updated
3836 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
3837 * the temporary variables holding the fixnums are cleared, and the process
3838 * starts over again. If for example fixnums were able to store five decimal
3839 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
3840 * and the result was computed as 12345 * 100000 + 67890. In other words,
3841 * only every five digits two bignum operations were performed.
3844 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
3846 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
3848 /* Caller is responsible for checking that the return value is in range
3849 for the given radix, which should be <= 36. */
3851 char_decimal_value (scm_t_uint32 c
)
3853 /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
3854 that's certainly above any valid decimal, so we take advantage of
3855 that to elide some tests. */
3856 unsigned int d
= (unsigned int) uc_decimal_value (c
);
3858 /* If that failed, try extended hexadecimals, then. Only accept ascii
3863 if (c
>= (scm_t_uint32
) 'a')
3864 d
= c
- (scm_t_uint32
)'a' + 10U;
3870 mem2uinteger (SCM mem
, unsigned int *p_idx
,
3871 unsigned int radix
, enum t_exactness
*p_exactness
)
3873 unsigned int idx
= *p_idx
;
3874 unsigned int hash_seen
= 0;
3875 scm_t_bits shift
= 1;
3877 unsigned int digit_value
;
3880 size_t len
= scm_i_string_length (mem
);
3885 c
= scm_i_string_ref (mem
, idx
);
3886 digit_value
= char_decimal_value (c
);
3887 if (digit_value
>= radix
)
3891 result
= SCM_I_MAKINUM (digit_value
);
3894 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
3904 digit_value
= char_decimal_value (c
);
3905 /* This check catches non-decimals in addition to out-of-range
3907 if (digit_value
>= radix
)
3912 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
3914 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3916 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3923 shift
= shift
* radix
;
3924 add
= add
* radix
+ digit_value
;
3929 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3931 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3935 *p_exactness
= INEXACT
;
3941 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
3942 * covers the parts of the rules that start at a potential point. The value
3943 * of the digits up to the point have been parsed by the caller and are given
3944 * in variable result. The content of *p_exactness indicates, whether a hash
3945 * has already been seen in the digits before the point.
3948 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
3951 mem2decimal_from_point (SCM result
, SCM mem
,
3952 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
3954 unsigned int idx
= *p_idx
;
3955 enum t_exactness x
= *p_exactness
;
3956 size_t len
= scm_i_string_length (mem
);
3961 if (scm_i_string_ref (mem
, idx
) == '.')
3963 scm_t_bits shift
= 1;
3965 unsigned int digit_value
;
3966 SCM big_shift
= SCM_INUM1
;
3971 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
3972 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
3977 digit_value
= DIGIT2UINT (c
);
3988 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
3990 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
3991 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3993 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
4001 add
= add
* 10 + digit_value
;
4007 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
4008 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
4009 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
4012 result
= scm_divide (result
, big_shift
);
4014 /* We've seen a decimal point, thus the value is implicitly inexact. */
4026 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
4028 switch (scm_i_string_ref (mem
, idx
))
4040 c
= scm_i_string_ref (mem
, idx
);
4048 c
= scm_i_string_ref (mem
, idx
);
4057 c
= scm_i_string_ref (mem
, idx
);
4062 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4066 exponent
= DIGIT2UINT (c
);
4069 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
4070 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4073 if (exponent
<= SCM_MAXEXP
)
4074 exponent
= exponent
* 10 + DIGIT2UINT (c
);
4080 if (exponent
> SCM_MAXEXP
)
4082 size_t exp_len
= idx
- start
;
4083 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
4084 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
4085 scm_out_of_range ("string->number", exp_num
);
4088 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
4090 result
= scm_product (result
, e
);
4092 result
= scm_divide2real (result
, e
);
4094 /* We've seen an exponent, thus the value is implicitly inexact. */
4112 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
4115 mem2ureal (SCM mem
, unsigned int *p_idx
,
4116 unsigned int radix
, enum t_exactness
*p_exactness
)
4118 unsigned int idx
= *p_idx
;
4120 size_t len
= scm_i_string_length (mem
);
4122 /* Start off believing that the number will be exact. This changes
4123 to INEXACT if we see a decimal point or a hash. */
4124 enum t_exactness x
= EXACT
;
4129 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
4135 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
4137 /* Cobble up the fractional part. We might want to set the
4138 NaN's mantissa from it. */
4140 mem2uinteger (mem
, &idx
, 10, &x
);
4145 if (scm_i_string_ref (mem
, idx
) == '.')
4149 else if (idx
+ 1 == len
)
4151 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
4154 result
= mem2decimal_from_point (SCM_INUM0
, mem
,
4161 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
4162 if (scm_is_false (uinteger
))
4167 else if (scm_i_string_ref (mem
, idx
) == '/')
4175 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
4176 if (scm_is_false (divisor
))
4179 /* both are int/big here, I assume */
4180 result
= scm_i_make_ratio (uinteger
, divisor
);
4182 else if (radix
== 10)
4184 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
4185 if (scm_is_false (result
))
4194 /* Update *p_exactness if the number just read was inexact. This is
4195 important for complex numbers, so that a complex number is
4196 treated as inexact overall if either its real or imaginary part
4202 /* When returning an inexact zero, make sure it is represented as a
4203 floating point value so that we can change its sign.
4205 if (scm_is_eq (result
, SCM_INUM0
) && *p_exactness
== INEXACT
)
4206 result
= scm_from_double (0.0);
4212 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
4215 mem2complex (SCM mem
, unsigned int idx
,
4216 unsigned int radix
, enum t_exactness
*p_exactness
)
4221 size_t len
= scm_i_string_length (mem
);
4226 c
= scm_i_string_ref (mem
, idx
);
4241 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
4242 if (scm_is_false (ureal
))
4244 /* input must be either +i or -i */
4249 if (scm_i_string_ref (mem
, idx
) == 'i'
4250 || scm_i_string_ref (mem
, idx
) == 'I')
4256 return scm_make_rectangular (SCM_INUM0
, SCM_I_MAKINUM (sign
));
4263 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
4264 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
4269 c
= scm_i_string_ref (mem
, idx
);
4273 /* either +<ureal>i or -<ureal>i */
4280 return scm_make_rectangular (SCM_INUM0
, ureal
);
4283 /* polar input: <real>@<real>. */
4294 c
= scm_i_string_ref (mem
, idx
);
4312 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
4313 if (scm_is_false (angle
))
4318 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
4319 angle
= scm_difference (angle
, SCM_UNDEFINED
);
4321 result
= scm_make_polar (ureal
, angle
);
4326 /* expecting input matching <real>[+-]<ureal>?i */
4333 int sign
= (c
== '+') ? 1 : -1;
4334 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
4336 if (scm_is_false (imag
))
4337 imag
= SCM_I_MAKINUM (sign
);
4338 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
4339 imag
= scm_difference (imag
, SCM_UNDEFINED
);
4343 if (scm_i_string_ref (mem
, idx
) != 'i'
4344 && scm_i_string_ref (mem
, idx
) != 'I')
4351 return scm_make_rectangular (ureal
, imag
);
4360 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
4362 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
4365 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
4367 unsigned int idx
= 0;
4368 unsigned int radix
= NO_RADIX
;
4369 enum t_exactness forced_x
= NO_EXACTNESS
;
4370 enum t_exactness implicit_x
= EXACT
;
4372 size_t len
= scm_i_string_length (mem
);
4374 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
4375 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
4377 switch (scm_i_string_ref (mem
, idx
+ 1))
4380 if (radix
!= NO_RADIX
)
4385 if (radix
!= NO_RADIX
)
4390 if (forced_x
!= NO_EXACTNESS
)
4395 if (forced_x
!= NO_EXACTNESS
)
4400 if (radix
!= NO_RADIX
)
4405 if (radix
!= NO_RADIX
)
4415 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
4416 if (radix
== NO_RADIX
)
4417 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
4419 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
4421 if (scm_is_false (result
))
4427 if (SCM_INEXACTP (result
))
4428 return scm_inexact_to_exact (result
);
4432 if (SCM_INEXACTP (result
))
4435 return scm_exact_to_inexact (result
);
4438 if (implicit_x
== INEXACT
)
4440 if (SCM_INEXACTP (result
))
4443 return scm_exact_to_inexact (result
);
4451 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
4452 unsigned int default_radix
)
4454 SCM str
= scm_from_locale_stringn (mem
, len
);
4456 return scm_i_string_to_number (str
, default_radix
);
4460 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
4461 (SCM string
, SCM radix
),
4462 "Return a number of the maximally precise representation\n"
4463 "expressed by the given @var{string}. @var{radix} must be an\n"
4464 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
4465 "is a default radix that may be overridden by an explicit radix\n"
4466 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
4467 "supplied, then the default radix is 10. If string is not a\n"
4468 "syntactically valid notation for a number, then\n"
4469 "@code{string->number} returns @code{#f}.")
4470 #define FUNC_NAME s_scm_string_to_number
4474 SCM_VALIDATE_STRING (1, string
);
4476 if (SCM_UNBNDP (radix
))
4479 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
4481 answer
= scm_i_string_to_number (string
, base
);
4482 scm_remember_upto_here_1 (string
);
4488 /*** END strs->nums ***/
4491 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
4493 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
4495 #define FUNC_NAME s_scm_number_p
4497 return scm_from_bool (SCM_NUMBERP (x
));
4501 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
4503 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
4504 "otherwise. Note that the sets of real, rational and integer\n"
4505 "values form subsets of the set of complex numbers, i. e. the\n"
4506 "predicate will also be fulfilled if @var{x} is a real,\n"
4507 "rational or integer number.")
4508 #define FUNC_NAME s_scm_complex_p
4510 /* all numbers are complex. */
4511 return scm_number_p (x
);
4515 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
4517 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
4518 "otherwise. Note that the set of integer values forms a subset of\n"
4519 "the set of real numbers, i. e. the predicate will also be\n"
4520 "fulfilled if @var{x} is an integer number.")
4521 #define FUNC_NAME s_scm_real_p
4523 return scm_from_bool
4524 (SCM_I_INUMP (x
) || SCM_REALP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
));
4528 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
4530 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
4531 "otherwise. Note that the set of integer values forms a subset of\n"
4532 "the set of rational numbers, i. e. the predicate will also be\n"
4533 "fulfilled if @var{x} is an integer number.")
4534 #define FUNC_NAME s_scm_rational_p
4536 if (SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
))
4538 else if (SCM_REALP (x
))
4539 /* due to their limited precision, finite floating point numbers are
4540 rational as well. (finite means neither infinity nor a NaN) */
4541 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
4547 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
4549 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
4551 #define FUNC_NAME s_scm_integer_p
4553 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4555 else if (SCM_REALP (x
))
4557 double val
= SCM_REAL_VALUE (x
);
4558 return scm_from_bool (!isinf (val
) && (val
== floor (val
)));
4566 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
4567 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
4568 (SCM x
, SCM y
, SCM rest
),
4569 "Return @code{#t} if all parameters are numerically equal.")
4570 #define FUNC_NAME s_scm_i_num_eq_p
4572 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4574 while (!scm_is_null (rest
))
4576 if (scm_is_false (scm_num_eq_p (x
, y
)))
4580 rest
= scm_cdr (rest
);
4582 return scm_num_eq_p (x
, y
);
4586 scm_num_eq_p (SCM x
, SCM y
)
4589 if (SCM_I_INUMP (x
))
4591 scm_t_signed_bits xx
= SCM_I_INUM (x
);
4592 if (SCM_I_INUMP (y
))
4594 scm_t_signed_bits yy
= SCM_I_INUM (y
);
4595 return scm_from_bool (xx
== yy
);
4597 else if (SCM_BIGP (y
))
4599 else if (SCM_REALP (y
))
4601 /* On a 32-bit system an inum fits a double, we can cast the inum
4602 to a double and compare.
4604 But on a 64-bit system an inum is bigger than a double and
4605 casting it to a double (call that dxx) will round. dxx is at
4606 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
4607 an integer and fits a long. So we cast yy to a long and
4608 compare with plain xx.
4610 An alternative (for any size system actually) would be to check
4611 yy is an integer (with floor) and is in range of an inum
4612 (compare against appropriate powers of 2) then test
4613 xx==(scm_t_signed_bits)yy. It's just a matter of which
4614 casts/comparisons might be fastest or easiest for the cpu. */
4616 double yy
= SCM_REAL_VALUE (y
);
4617 return scm_from_bool ((double) xx
== yy
4618 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
4619 || xx
== (scm_t_signed_bits
) yy
));
4621 else if (SCM_COMPLEXP (y
))
4622 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
4623 && (0.0 == SCM_COMPLEX_IMAG (y
)));
4624 else if (SCM_FRACTIONP (y
))
4627 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4629 else if (SCM_BIGP (x
))
4631 if (SCM_I_INUMP (y
))
4633 else if (SCM_BIGP (y
))
4635 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4636 scm_remember_upto_here_2 (x
, y
);
4637 return scm_from_bool (0 == cmp
);
4639 else if (SCM_REALP (y
))
4642 if (isnan (SCM_REAL_VALUE (y
)))
4644 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
4645 scm_remember_upto_here_1 (x
);
4646 return scm_from_bool (0 == cmp
);
4648 else if (SCM_COMPLEXP (y
))
4651 if (0.0 != SCM_COMPLEX_IMAG (y
))
4653 if (isnan (SCM_COMPLEX_REAL (y
)))
4655 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
4656 scm_remember_upto_here_1 (x
);
4657 return scm_from_bool (0 == cmp
);
4659 else if (SCM_FRACTIONP (y
))
4662 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4664 else if (SCM_REALP (x
))
4666 double xx
= SCM_REAL_VALUE (x
);
4667 if (SCM_I_INUMP (y
))
4669 /* see comments with inum/real above */
4670 scm_t_signed_bits yy
= SCM_I_INUM (y
);
4671 return scm_from_bool (xx
== (double) yy
4672 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
4673 || (scm_t_signed_bits
) xx
== yy
));
4675 else if (SCM_BIGP (y
))
4678 if (isnan (SCM_REAL_VALUE (x
)))
4680 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
4681 scm_remember_upto_here_1 (y
);
4682 return scm_from_bool (0 == cmp
);
4684 else if (SCM_REALP (y
))
4685 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
4686 else if (SCM_COMPLEXP (y
))
4687 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
4688 && (0.0 == SCM_COMPLEX_IMAG (y
)));
4689 else if (SCM_FRACTIONP (y
))
4691 double xx
= SCM_REAL_VALUE (x
);
4695 return scm_from_bool (xx
< 0.0);
4696 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4700 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4702 else if (SCM_COMPLEXP (x
))
4704 if (SCM_I_INUMP (y
))
4705 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
4706 && (SCM_COMPLEX_IMAG (x
) == 0.0));
4707 else if (SCM_BIGP (y
))
4710 if (0.0 != SCM_COMPLEX_IMAG (x
))
4712 if (isnan (SCM_COMPLEX_REAL (x
)))
4714 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
4715 scm_remember_upto_here_1 (y
);
4716 return scm_from_bool (0 == cmp
);
4718 else if (SCM_REALP (y
))
4719 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
4720 && (SCM_COMPLEX_IMAG (x
) == 0.0));
4721 else if (SCM_COMPLEXP (y
))
4722 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
4723 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
4724 else if (SCM_FRACTIONP (y
))
4727 if (SCM_COMPLEX_IMAG (x
) != 0.0)
4729 xx
= SCM_COMPLEX_REAL (x
);
4733 return scm_from_bool (xx
< 0.0);
4734 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4738 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4740 else if (SCM_FRACTIONP (x
))
4742 if (SCM_I_INUMP (y
))
4744 else if (SCM_BIGP (y
))
4746 else if (SCM_REALP (y
))
4748 double yy
= SCM_REAL_VALUE (y
);
4752 return scm_from_bool (0.0 < yy
);
4753 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4756 else if (SCM_COMPLEXP (y
))
4759 if (SCM_COMPLEX_IMAG (y
) != 0.0)
4761 yy
= SCM_COMPLEX_REAL (y
);
4765 return scm_from_bool (0.0 < yy
);
4766 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4769 else if (SCM_FRACTIONP (y
))
4770 return scm_i_fraction_equalp (x
, y
);
4772 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4775 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
4779 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
4780 done are good for inums, but for bignums an answer can almost always be
4781 had by just examining a few high bits of the operands, as done by GMP in
4782 mpq_cmp. flonum/frac compares likewise, but with the slight complication
4783 of the float exponent to take into account. */
4785 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
4786 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
4787 (SCM x
, SCM y
, SCM rest
),
4788 "Return @code{#t} if the list of parameters is monotonically\n"
4790 #define FUNC_NAME s_scm_i_num_less_p
4792 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4794 while (!scm_is_null (rest
))
4796 if (scm_is_false (scm_less_p (x
, y
)))
4800 rest
= scm_cdr (rest
);
4802 return scm_less_p (x
, y
);
4806 scm_less_p (SCM x
, SCM y
)
4809 if (SCM_I_INUMP (x
))
4811 scm_t_inum xx
= SCM_I_INUM (x
);
4812 if (SCM_I_INUMP (y
))
4814 scm_t_inum yy
= SCM_I_INUM (y
);
4815 return scm_from_bool (xx
< yy
);
4817 else if (SCM_BIGP (y
))
4819 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4820 scm_remember_upto_here_1 (y
);
4821 return scm_from_bool (sgn
> 0);
4823 else if (SCM_REALP (y
))
4824 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
4825 else if (SCM_FRACTIONP (y
))
4827 /* "x < a/b" becomes "x*b < a" */
4829 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
4830 y
= SCM_FRACTION_NUMERATOR (y
);
4834 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4836 else if (SCM_BIGP (x
))
4838 if (SCM_I_INUMP (y
))
4840 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4841 scm_remember_upto_here_1 (x
);
4842 return scm_from_bool (sgn
< 0);
4844 else if (SCM_BIGP (y
))
4846 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4847 scm_remember_upto_here_2 (x
, y
);
4848 return scm_from_bool (cmp
< 0);
4850 else if (SCM_REALP (y
))
4853 if (isnan (SCM_REAL_VALUE (y
)))
4855 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
4856 scm_remember_upto_here_1 (x
);
4857 return scm_from_bool (cmp
< 0);
4859 else if (SCM_FRACTIONP (y
))
4862 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4864 else if (SCM_REALP (x
))
4866 if (SCM_I_INUMP (y
))
4867 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
4868 else if (SCM_BIGP (y
))
4871 if (isnan (SCM_REAL_VALUE (x
)))
4873 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
4874 scm_remember_upto_here_1 (y
);
4875 return scm_from_bool (cmp
> 0);
4877 else if (SCM_REALP (y
))
4878 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
4879 else if (SCM_FRACTIONP (y
))
4881 double xx
= SCM_REAL_VALUE (x
);
4885 return scm_from_bool (xx
< 0.0);
4886 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4890 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4892 else if (SCM_FRACTIONP (x
))
4894 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
4896 /* "a/b < y" becomes "a < y*b" */
4897 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
4898 x
= SCM_FRACTION_NUMERATOR (x
);
4901 else if (SCM_REALP (y
))
4903 double yy
= SCM_REAL_VALUE (y
);
4907 return scm_from_bool (0.0 < yy
);
4908 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4911 else if (SCM_FRACTIONP (y
))
4913 /* "a/b < c/d" becomes "a*d < c*b" */
4914 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
4915 SCM_FRACTION_DENOMINATOR (y
));
4916 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
4917 SCM_FRACTION_DENOMINATOR (x
));
4923 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4926 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
4930 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
4931 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
4932 (SCM x
, SCM y
, SCM rest
),
4933 "Return @code{#t} if the list of parameters is monotonically\n"
4935 #define FUNC_NAME s_scm_i_num_gr_p
4937 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4939 while (!scm_is_null (rest
))
4941 if (scm_is_false (scm_gr_p (x
, y
)))
4945 rest
= scm_cdr (rest
);
4947 return scm_gr_p (x
, y
);
4950 #define FUNC_NAME s_scm_i_num_gr_p
4952 scm_gr_p (SCM x
, SCM y
)
4954 if (!SCM_NUMBERP (x
))
4955 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
4956 else if (!SCM_NUMBERP (y
))
4957 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
4959 return scm_less_p (y
, x
);
4964 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
4965 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
4966 (SCM x
, SCM y
, SCM rest
),
4967 "Return @code{#t} if the list of parameters is monotonically\n"
4969 #define FUNC_NAME s_scm_i_num_leq_p
4971 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4973 while (!scm_is_null (rest
))
4975 if (scm_is_false (scm_leq_p (x
, y
)))
4979 rest
= scm_cdr (rest
);
4981 return scm_leq_p (x
, y
);
4984 #define FUNC_NAME s_scm_i_num_leq_p
4986 scm_leq_p (SCM x
, SCM y
)
4988 if (!SCM_NUMBERP (x
))
4989 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
4990 else if (!SCM_NUMBERP (y
))
4991 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
4992 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
4995 return scm_not (scm_less_p (y
, x
));
5000 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
5001 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
5002 (SCM x
, SCM y
, SCM rest
),
5003 "Return @code{#t} if the list of parameters is monotonically\n"
5005 #define FUNC_NAME s_scm_i_num_geq_p
5007 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
5009 while (!scm_is_null (rest
))
5011 if (scm_is_false (scm_geq_p (x
, y
)))
5015 rest
= scm_cdr (rest
);
5017 return scm_geq_p (x
, y
);
5020 #define FUNC_NAME s_scm_i_num_geq_p
5022 scm_geq_p (SCM x
, SCM y
)
5024 if (!SCM_NUMBERP (x
))
5025 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
5026 else if (!SCM_NUMBERP (y
))
5027 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
5028 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
5031 return scm_not (scm_less_p (x
, y
));
5036 SCM_PRIMITIVE_GENERIC (scm_zero_p
, "zero?", 1, 0, 0,
5038 "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
5040 #define FUNC_NAME s_scm_zero_p
5042 if (SCM_I_INUMP (z
))
5043 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
5044 else if (SCM_BIGP (z
))
5046 else if (SCM_REALP (z
))
5047 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
5048 else if (SCM_COMPLEXP (z
))
5049 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
5050 && SCM_COMPLEX_IMAG (z
) == 0.0);
5051 else if (SCM_FRACTIONP (z
))
5054 SCM_WTA_DISPATCH_1 (g_scm_zero_p
, z
, SCM_ARG1
, s_scm_zero_p
);
5059 SCM_PRIMITIVE_GENERIC (scm_positive_p
, "positive?", 1, 0, 0,
5061 "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
5063 #define FUNC_NAME s_scm_positive_p
5065 if (SCM_I_INUMP (x
))
5066 return scm_from_bool (SCM_I_INUM (x
) > 0);
5067 else if (SCM_BIGP (x
))
5069 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5070 scm_remember_upto_here_1 (x
);
5071 return scm_from_bool (sgn
> 0);
5073 else if (SCM_REALP (x
))
5074 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
5075 else if (SCM_FRACTIONP (x
))
5076 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
5078 SCM_WTA_DISPATCH_1 (g_scm_positive_p
, x
, SCM_ARG1
, s_scm_positive_p
);
5083 SCM_PRIMITIVE_GENERIC (scm_negative_p
, "negative?", 1, 0, 0,
5085 "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
5087 #define FUNC_NAME s_scm_negative_p
5089 if (SCM_I_INUMP (x
))
5090 return scm_from_bool (SCM_I_INUM (x
) < 0);
5091 else if (SCM_BIGP (x
))
5093 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5094 scm_remember_upto_here_1 (x
);
5095 return scm_from_bool (sgn
< 0);
5097 else if (SCM_REALP (x
))
5098 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
5099 else if (SCM_FRACTIONP (x
))
5100 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
5102 SCM_WTA_DISPATCH_1 (g_scm_negative_p
, x
, SCM_ARG1
, s_scm_negative_p
);
5107 /* scm_min and scm_max return an inexact when either argument is inexact, as
5108 required by r5rs. On that basis, for exact/inexact combinations the
5109 exact is converted to inexact to compare and possibly return. This is
5110 unlike scm_less_p above which takes some trouble to preserve all bits in
5111 its test, such trouble is not required for min and max. */
5113 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
5114 (SCM x
, SCM y
, SCM rest
),
5115 "Return the maximum of all parameter values.")
5116 #define FUNC_NAME s_scm_i_max
5118 while (!scm_is_null (rest
))
5119 { x
= scm_max (x
, y
);
5121 rest
= scm_cdr (rest
);
5123 return scm_max (x
, y
);
5127 #define s_max s_scm_i_max
5128 #define g_max g_scm_i_max
5131 scm_max (SCM x
, SCM y
)
5136 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
5137 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
5140 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
5143 if (SCM_I_INUMP (x
))
5145 scm_t_inum xx
= SCM_I_INUM (x
);
5146 if (SCM_I_INUMP (y
))
5148 scm_t_inum yy
= SCM_I_INUM (y
);
5149 return (xx
< yy
) ? y
: x
;
5151 else if (SCM_BIGP (y
))
5153 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5154 scm_remember_upto_here_1 (y
);
5155 return (sgn
< 0) ? x
: y
;
5157 else if (SCM_REALP (y
))
5160 double yyd
= SCM_REAL_VALUE (y
);
5163 return scm_from_double (xxd
);
5164 /* If y is a NaN, then "==" is false and we return the NaN */
5165 else if (SCM_LIKELY (!(xxd
== yyd
)))
5167 /* Handle signed zeroes properly */
5173 else if (SCM_FRACTIONP (y
))
5176 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
5179 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5181 else if (SCM_BIGP (x
))
5183 if (SCM_I_INUMP (y
))
5185 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5186 scm_remember_upto_here_1 (x
);
5187 return (sgn
< 0) ? y
: x
;
5189 else if (SCM_BIGP (y
))
5191 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
5192 scm_remember_upto_here_2 (x
, y
);
5193 return (cmp
> 0) ? x
: y
;
5195 else if (SCM_REALP (y
))
5197 /* if y==NaN then xx>yy is false, so we return the NaN y */
5200 xx
= scm_i_big2dbl (x
);
5201 yy
= SCM_REAL_VALUE (y
);
5202 return (xx
> yy
? scm_from_double (xx
) : y
);
5204 else if (SCM_FRACTIONP (y
))
5209 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5211 else if (SCM_REALP (x
))
5213 if (SCM_I_INUMP (y
))
5215 scm_t_inum yy
= SCM_I_INUM (y
);
5216 double xxd
= SCM_REAL_VALUE (x
);
5220 return scm_from_double (yyd
);
5221 /* If x is a NaN, then "==" is false and we return the NaN */
5222 else if (SCM_LIKELY (!(xxd
== yyd
)))
5224 /* Handle signed zeroes properly */
5230 else if (SCM_BIGP (y
))
5235 else if (SCM_REALP (y
))
5237 double xx
= SCM_REAL_VALUE (x
);
5238 double yy
= SCM_REAL_VALUE (y
);
5240 /* For purposes of max: +inf.0 > nan > everything else, per R6RS */
5243 else if (SCM_LIKELY (xx
< yy
))
5245 /* If neither (xx > yy) nor (xx < yy), then
5246 either they're equal or one is a NaN */
5247 else if (SCM_UNLIKELY (isnan (xx
)))
5248 return (isinf (yy
) == 1) ? y
: x
;
5249 else if (SCM_UNLIKELY (isnan (yy
)))
5250 return (isinf (xx
) == 1) ? x
: y
;
5251 /* xx == yy, but handle signed zeroes properly */
5252 else if (double_is_non_negative_zero (yy
))
5257 else if (SCM_FRACTIONP (y
))
5259 double yy
= scm_i_fraction2double (y
);
5260 double xx
= SCM_REAL_VALUE (x
);
5261 return (xx
< yy
) ? scm_from_double (yy
) : x
;
5264 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5266 else if (SCM_FRACTIONP (x
))
5268 if (SCM_I_INUMP (y
))
5272 else if (SCM_BIGP (y
))
5276 else if (SCM_REALP (y
))
5278 double xx
= scm_i_fraction2double (x
);
5279 /* if y==NaN then ">" is false, so we return the NaN y */
5280 return (xx
> SCM_REAL_VALUE (y
)) ? scm_from_double (xx
) : y
;
5282 else if (SCM_FRACTIONP (y
))
5287 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5290 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
5294 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
5295 (SCM x
, SCM y
, SCM rest
),
5296 "Return the minimum of all parameter values.")
5297 #define FUNC_NAME s_scm_i_min
5299 while (!scm_is_null (rest
))
5300 { x
= scm_min (x
, y
);
5302 rest
= scm_cdr (rest
);
5304 return scm_min (x
, y
);
5308 #define s_min s_scm_i_min
5309 #define g_min g_scm_i_min
5312 scm_min (SCM x
, SCM y
)
5317 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
5318 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
5321 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
5324 if (SCM_I_INUMP (x
))
5326 scm_t_inum xx
= SCM_I_INUM (x
);
5327 if (SCM_I_INUMP (y
))
5329 scm_t_inum yy
= SCM_I_INUM (y
);
5330 return (xx
< yy
) ? x
: y
;
5332 else if (SCM_BIGP (y
))
5334 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5335 scm_remember_upto_here_1 (y
);
5336 return (sgn
< 0) ? y
: x
;
5338 else if (SCM_REALP (y
))
5341 /* if y==NaN then "<" is false and we return NaN */
5342 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
5344 else if (SCM_FRACTIONP (y
))
5347 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
5350 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5352 else if (SCM_BIGP (x
))
5354 if (SCM_I_INUMP (y
))
5356 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5357 scm_remember_upto_here_1 (x
);
5358 return (sgn
< 0) ? x
: y
;
5360 else if (SCM_BIGP (y
))
5362 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
5363 scm_remember_upto_here_2 (x
, y
);
5364 return (cmp
> 0) ? y
: x
;
5366 else if (SCM_REALP (y
))
5368 /* if y==NaN then xx<yy is false, so we return the NaN y */
5371 xx
= scm_i_big2dbl (x
);
5372 yy
= SCM_REAL_VALUE (y
);
5373 return (xx
< yy
? scm_from_double (xx
) : y
);
5375 else if (SCM_FRACTIONP (y
))
5380 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5382 else if (SCM_REALP (x
))
5384 if (SCM_I_INUMP (y
))
5386 double z
= SCM_I_INUM (y
);
5387 /* if x==NaN then "<" is false and we return NaN */
5388 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
5390 else if (SCM_BIGP (y
))
5395 else if (SCM_REALP (y
))
5397 double xx
= SCM_REAL_VALUE (x
);
5398 double yy
= SCM_REAL_VALUE (y
);
5400 /* For purposes of min: -inf.0 < nan < everything else, per R6RS */
5403 else if (SCM_LIKELY (xx
> yy
))
5405 /* If neither (xx < yy) nor (xx > yy), then
5406 either they're equal or one is a NaN */
5407 else if (SCM_UNLIKELY (isnan (xx
)))
5408 return (isinf (yy
) == -1) ? y
: x
;
5409 else if (SCM_UNLIKELY (isnan (yy
)))
5410 return (isinf (xx
) == -1) ? x
: y
;
5411 /* xx == yy, but handle signed zeroes properly */
5412 else if (double_is_non_negative_zero (xx
))
5417 else if (SCM_FRACTIONP (y
))
5419 double yy
= scm_i_fraction2double (y
);
5420 double xx
= SCM_REAL_VALUE (x
);
5421 return (yy
< xx
) ? scm_from_double (yy
) : x
;
5424 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5426 else if (SCM_FRACTIONP (x
))
5428 if (SCM_I_INUMP (y
))
5432 else if (SCM_BIGP (y
))
5436 else if (SCM_REALP (y
))
5438 double xx
= scm_i_fraction2double (x
);
5439 /* if y==NaN then "<" is false, so we return the NaN y */
5440 return (xx
< SCM_REAL_VALUE (y
)) ? scm_from_double (xx
) : y
;
5442 else if (SCM_FRACTIONP (y
))
5447 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5450 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
5454 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
5455 (SCM x
, SCM y
, SCM rest
),
5456 "Return the sum of all parameter values. Return 0 if called without\n"
5458 #define FUNC_NAME s_scm_i_sum
5460 while (!scm_is_null (rest
))
5461 { x
= scm_sum (x
, y
);
5463 rest
= scm_cdr (rest
);
5465 return scm_sum (x
, y
);
5469 #define s_sum s_scm_i_sum
5470 #define g_sum g_scm_i_sum
5473 scm_sum (SCM x
, SCM y
)
5475 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5477 if (SCM_NUMBERP (x
)) return x
;
5478 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
5479 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
5482 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5484 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5486 scm_t_inum xx
= SCM_I_INUM (x
);
5487 scm_t_inum yy
= SCM_I_INUM (y
);
5488 scm_t_inum z
= xx
+ yy
;
5489 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_inum2big (z
);
5491 else if (SCM_BIGP (y
))
5496 else if (SCM_REALP (y
))
5498 scm_t_inum xx
= SCM_I_INUM (x
);
5499 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
5501 else if (SCM_COMPLEXP (y
))
5503 scm_t_inum xx
= SCM_I_INUM (x
);
5504 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
5505 SCM_COMPLEX_IMAG (y
));
5507 else if (SCM_FRACTIONP (y
))
5508 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
5509 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
5510 SCM_FRACTION_DENOMINATOR (y
));
5512 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5513 } else if (SCM_BIGP (x
))
5515 if (SCM_I_INUMP (y
))
5520 inum
= SCM_I_INUM (y
);
5523 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5526 SCM result
= scm_i_mkbig ();
5527 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
5528 scm_remember_upto_here_1 (x
);
5529 /* we know the result will have to be a bignum */
5532 return scm_i_normbig (result
);
5536 SCM result
= scm_i_mkbig ();
5537 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
5538 scm_remember_upto_here_1 (x
);
5539 /* we know the result will have to be a bignum */
5542 return scm_i_normbig (result
);
5545 else if (SCM_BIGP (y
))
5547 SCM result
= scm_i_mkbig ();
5548 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5549 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5550 mpz_add (SCM_I_BIG_MPZ (result
),
5553 scm_remember_upto_here_2 (x
, y
);
5554 /* we know the result will have to be a bignum */
5557 return scm_i_normbig (result
);
5559 else if (SCM_REALP (y
))
5561 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
5562 scm_remember_upto_here_1 (x
);
5563 return scm_from_double (result
);
5565 else if (SCM_COMPLEXP (y
))
5567 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
5568 + SCM_COMPLEX_REAL (y
));
5569 scm_remember_upto_here_1 (x
);
5570 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
5572 else if (SCM_FRACTIONP (y
))
5573 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
5574 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
5575 SCM_FRACTION_DENOMINATOR (y
));
5577 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5579 else if (SCM_REALP (x
))
5581 if (SCM_I_INUMP (y
))
5582 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
5583 else if (SCM_BIGP (y
))
5585 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
5586 scm_remember_upto_here_1 (y
);
5587 return scm_from_double (result
);
5589 else if (SCM_REALP (y
))
5590 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
5591 else if (SCM_COMPLEXP (y
))
5592 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
5593 SCM_COMPLEX_IMAG (y
));
5594 else if (SCM_FRACTIONP (y
))
5595 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
5597 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5599 else if (SCM_COMPLEXP (x
))
5601 if (SCM_I_INUMP (y
))
5602 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
5603 SCM_COMPLEX_IMAG (x
));
5604 else if (SCM_BIGP (y
))
5606 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
5607 + SCM_COMPLEX_REAL (x
));
5608 scm_remember_upto_here_1 (y
);
5609 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
5611 else if (SCM_REALP (y
))
5612 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
5613 SCM_COMPLEX_IMAG (x
));
5614 else if (SCM_COMPLEXP (y
))
5615 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
5616 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
5617 else if (SCM_FRACTIONP (y
))
5618 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
5619 SCM_COMPLEX_IMAG (x
));
5621 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5623 else if (SCM_FRACTIONP (x
))
5625 if (SCM_I_INUMP (y
))
5626 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
5627 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
5628 SCM_FRACTION_DENOMINATOR (x
));
5629 else if (SCM_BIGP (y
))
5630 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
5631 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
5632 SCM_FRACTION_DENOMINATOR (x
));
5633 else if (SCM_REALP (y
))
5634 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
5635 else if (SCM_COMPLEXP (y
))
5636 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
5637 SCM_COMPLEX_IMAG (y
));
5638 else if (SCM_FRACTIONP (y
))
5639 /* a/b + c/d = (ad + bc) / bd */
5640 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5641 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
5642 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
5644 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5647 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
5651 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
5653 "Return @math{@var{x}+1}.")
5654 #define FUNC_NAME s_scm_oneplus
5656 return scm_sum (x
, SCM_INUM1
);
5661 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
5662 (SCM x
, SCM y
, SCM rest
),
5663 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
5664 "the sum of all but the first argument are subtracted from the first\n"
5666 #define FUNC_NAME s_scm_i_difference
5668 while (!scm_is_null (rest
))
5669 { x
= scm_difference (x
, y
);
5671 rest
= scm_cdr (rest
);
5673 return scm_difference (x
, y
);
5677 #define s_difference s_scm_i_difference
5678 #define g_difference g_scm_i_difference
5681 scm_difference (SCM x
, SCM y
)
5682 #define FUNC_NAME s_difference
5684 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5687 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
5689 if (SCM_I_INUMP (x
))
5691 scm_t_inum xx
= -SCM_I_INUM (x
);
5692 if (SCM_FIXABLE (xx
))
5693 return SCM_I_MAKINUM (xx
);
5695 return scm_i_inum2big (xx
);
5697 else if (SCM_BIGP (x
))
5698 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
5699 bignum, but negating that gives a fixnum. */
5700 return scm_i_normbig (scm_i_clonebig (x
, 0));
5701 else if (SCM_REALP (x
))
5702 return scm_from_double (-SCM_REAL_VALUE (x
));
5703 else if (SCM_COMPLEXP (x
))
5704 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
5705 -SCM_COMPLEX_IMAG (x
));
5706 else if (SCM_FRACTIONP (x
))
5707 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
5708 SCM_FRACTION_DENOMINATOR (x
));
5710 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
5713 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5715 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5717 scm_t_inum xx
= SCM_I_INUM (x
);
5718 scm_t_inum yy
= SCM_I_INUM (y
);
5719 scm_t_inum z
= xx
- yy
;
5720 if (SCM_FIXABLE (z
))
5721 return SCM_I_MAKINUM (z
);
5723 return scm_i_inum2big (z
);
5725 else if (SCM_BIGP (y
))
5727 /* inum-x - big-y */
5728 scm_t_inum xx
= SCM_I_INUM (x
);
5732 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
5733 bignum, but negating that gives a fixnum. */
5734 return scm_i_normbig (scm_i_clonebig (y
, 0));
5738 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5739 SCM result
= scm_i_mkbig ();
5742 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
5745 /* x - y == -(y + -x) */
5746 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
5747 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5749 scm_remember_upto_here_1 (y
);
5751 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
5752 /* we know the result will have to be a bignum */
5755 return scm_i_normbig (result
);
5758 else if (SCM_REALP (y
))
5760 scm_t_inum xx
= SCM_I_INUM (x
);
5761 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
5763 else if (SCM_COMPLEXP (y
))
5765 scm_t_inum xx
= SCM_I_INUM (x
);
5766 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
5767 - SCM_COMPLEX_IMAG (y
));
5769 else if (SCM_FRACTIONP (y
))
5770 /* a - b/c = (ac - b) / c */
5771 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5772 SCM_FRACTION_NUMERATOR (y
)),
5773 SCM_FRACTION_DENOMINATOR (y
));
5775 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5777 else if (SCM_BIGP (x
))
5779 if (SCM_I_INUMP (y
))
5781 /* big-x - inum-y */
5782 scm_t_inum yy
= SCM_I_INUM (y
);
5783 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5785 scm_remember_upto_here_1 (x
);
5787 return (SCM_FIXABLE (-yy
) ?
5788 SCM_I_MAKINUM (-yy
) : scm_from_inum (-yy
));
5791 SCM result
= scm_i_mkbig ();
5794 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
5796 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
5797 scm_remember_upto_here_1 (x
);
5799 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
5800 /* we know the result will have to be a bignum */
5803 return scm_i_normbig (result
);
5806 else if (SCM_BIGP (y
))
5808 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5809 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5810 SCM result
= scm_i_mkbig ();
5811 mpz_sub (SCM_I_BIG_MPZ (result
),
5814 scm_remember_upto_here_2 (x
, y
);
5815 /* we know the result will have to be a bignum */
5816 if ((sgn_x
== 1) && (sgn_y
== -1))
5818 if ((sgn_x
== -1) && (sgn_y
== 1))
5820 return scm_i_normbig (result
);
5822 else if (SCM_REALP (y
))
5824 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
5825 scm_remember_upto_here_1 (x
);
5826 return scm_from_double (result
);
5828 else if (SCM_COMPLEXP (y
))
5830 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
5831 - SCM_COMPLEX_REAL (y
));
5832 scm_remember_upto_here_1 (x
);
5833 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
5835 else if (SCM_FRACTIONP (y
))
5836 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5837 SCM_FRACTION_NUMERATOR (y
)),
5838 SCM_FRACTION_DENOMINATOR (y
));
5839 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5841 else if (SCM_REALP (x
))
5843 if (SCM_I_INUMP (y
))
5844 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
5845 else if (SCM_BIGP (y
))
5847 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
5848 scm_remember_upto_here_1 (x
);
5849 return scm_from_double (result
);
5851 else if (SCM_REALP (y
))
5852 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
5853 else if (SCM_COMPLEXP (y
))
5854 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
5855 -SCM_COMPLEX_IMAG (y
));
5856 else if (SCM_FRACTIONP (y
))
5857 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
5859 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5861 else if (SCM_COMPLEXP (x
))
5863 if (SCM_I_INUMP (y
))
5864 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
5865 SCM_COMPLEX_IMAG (x
));
5866 else if (SCM_BIGP (y
))
5868 double real_part
= (SCM_COMPLEX_REAL (x
)
5869 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
5870 scm_remember_upto_here_1 (x
);
5871 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
5873 else if (SCM_REALP (y
))
5874 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
5875 SCM_COMPLEX_IMAG (x
));
5876 else if (SCM_COMPLEXP (y
))
5877 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
5878 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
5879 else if (SCM_FRACTIONP (y
))
5880 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
5881 SCM_COMPLEX_IMAG (x
));
5883 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5885 else if (SCM_FRACTIONP (x
))
5887 if (SCM_I_INUMP (y
))
5888 /* a/b - c = (a - cb) / b */
5889 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
5890 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
5891 SCM_FRACTION_DENOMINATOR (x
));
5892 else if (SCM_BIGP (y
))
5893 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
5894 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
5895 SCM_FRACTION_DENOMINATOR (x
));
5896 else if (SCM_REALP (y
))
5897 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
5898 else if (SCM_COMPLEXP (y
))
5899 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
5900 -SCM_COMPLEX_IMAG (y
));
5901 else if (SCM_FRACTIONP (y
))
5902 /* a/b - c/d = (ad - bc) / bd */
5903 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5904 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
5905 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
5907 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5910 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
5915 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
5917 "Return @math{@var{x}-1}.")
5918 #define FUNC_NAME s_scm_oneminus
5920 return scm_difference (x
, SCM_INUM1
);
5925 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
5926 (SCM x
, SCM y
, SCM rest
),
5927 "Return the product of all arguments. If called without arguments,\n"
5929 #define FUNC_NAME s_scm_i_product
5931 while (!scm_is_null (rest
))
5932 { x
= scm_product (x
, y
);
5934 rest
= scm_cdr (rest
);
5936 return scm_product (x
, y
);
5940 #define s_product s_scm_i_product
5941 #define g_product g_scm_i_product
5944 scm_product (SCM x
, SCM y
)
5946 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5949 return SCM_I_MAKINUM (1L);
5950 else if (SCM_NUMBERP (x
))
5953 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
5956 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5961 xx
= SCM_I_INUM (x
);
5966 /* exact1 is the universal multiplicative identity */
5970 /* exact0 times a fixnum is exact0: optimize this case */
5971 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5973 /* if the other argument is inexact, the result is inexact,
5974 and we must do the multiplication in order to handle
5975 infinities and NaNs properly. */
5976 else if (SCM_REALP (y
))
5977 return scm_from_double (0.0 * SCM_REAL_VALUE (y
));
5978 else if (SCM_COMPLEXP (y
))
5979 return scm_c_make_rectangular (0.0 * SCM_COMPLEX_REAL (y
),
5980 0.0 * SCM_COMPLEX_IMAG (y
));
5981 /* we've already handled inexact numbers,
5982 so y must be exact, and we return exact0 */
5983 else if (SCM_NUMP (y
))
5986 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
5990 * This case is important for more than just optimization.
5991 * It handles the case of negating
5992 * (+ 1 most-positive-fixnum) aka (- most-negative-fixnum),
5993 * which is a bignum that must be changed back into a fixnum.
5994 * Failure to do so will cause the following to return #f:
5995 * (= most-negative-fixnum (* -1 (- most-negative-fixnum)))
5997 return scm_difference(y
, SCM_UNDEFINED
);
6001 if (SCM_LIKELY (SCM_I_INUMP (y
)))
6003 scm_t_inum yy
= SCM_I_INUM (y
);
6004 scm_t_inum kk
= xx
* yy
;
6005 SCM k
= SCM_I_MAKINUM (kk
);
6006 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
6010 SCM result
= scm_i_inum2big (xx
);
6011 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
6012 return scm_i_normbig (result
);
6015 else if (SCM_BIGP (y
))
6017 SCM result
= scm_i_mkbig ();
6018 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
6019 scm_remember_upto_here_1 (y
);
6022 else if (SCM_REALP (y
))
6023 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
6024 else if (SCM_COMPLEXP (y
))
6025 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
6026 xx
* SCM_COMPLEX_IMAG (y
));
6027 else if (SCM_FRACTIONP (y
))
6028 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
6029 SCM_FRACTION_DENOMINATOR (y
));
6031 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6033 else if (SCM_BIGP (x
))
6035 if (SCM_I_INUMP (y
))
6040 else if (SCM_BIGP (y
))
6042 SCM result
= scm_i_mkbig ();
6043 mpz_mul (SCM_I_BIG_MPZ (result
),
6046 scm_remember_upto_here_2 (x
, y
);
6049 else if (SCM_REALP (y
))
6051 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
6052 scm_remember_upto_here_1 (x
);
6053 return scm_from_double (result
);
6055 else if (SCM_COMPLEXP (y
))
6057 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
6058 scm_remember_upto_here_1 (x
);
6059 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
6060 z
* SCM_COMPLEX_IMAG (y
));
6062 else if (SCM_FRACTIONP (y
))
6063 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
6064 SCM_FRACTION_DENOMINATOR (y
));
6066 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6068 else if (SCM_REALP (x
))
6070 if (SCM_I_INUMP (y
))
6075 else if (SCM_BIGP (y
))
6077 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
6078 scm_remember_upto_here_1 (y
);
6079 return scm_from_double (result
);
6081 else if (SCM_REALP (y
))
6082 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
6083 else if (SCM_COMPLEXP (y
))
6084 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
6085 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
6086 else if (SCM_FRACTIONP (y
))
6087 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
6089 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6091 else if (SCM_COMPLEXP (x
))
6093 if (SCM_I_INUMP (y
))
6098 else if (SCM_BIGP (y
))
6100 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6101 scm_remember_upto_here_1 (y
);
6102 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
6103 z
* SCM_COMPLEX_IMAG (x
));
6105 else if (SCM_REALP (y
))
6106 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
6107 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
6108 else if (SCM_COMPLEXP (y
))
6110 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
6111 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
6112 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
6113 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
6115 else if (SCM_FRACTIONP (y
))
6117 double yy
= scm_i_fraction2double (y
);
6118 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
6119 yy
* SCM_COMPLEX_IMAG (x
));
6122 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6124 else if (SCM_FRACTIONP (x
))
6126 if (SCM_I_INUMP (y
))
6127 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
6128 SCM_FRACTION_DENOMINATOR (x
));
6129 else if (SCM_BIGP (y
))
6130 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
6131 SCM_FRACTION_DENOMINATOR (x
));
6132 else if (SCM_REALP (y
))
6133 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
6134 else if (SCM_COMPLEXP (y
))
6136 double xx
= scm_i_fraction2double (x
);
6137 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
6138 xx
* SCM_COMPLEX_IMAG (y
));
6140 else if (SCM_FRACTIONP (y
))
6141 /* a/b * c/d = ac / bd */
6142 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
6143 SCM_FRACTION_NUMERATOR (y
)),
6144 scm_product (SCM_FRACTION_DENOMINATOR (x
),
6145 SCM_FRACTION_DENOMINATOR (y
)));
6147 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6150 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
6153 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
6154 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
6155 #define ALLOW_DIVIDE_BY_ZERO
6156 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
6159 /* The code below for complex division is adapted from the GNU
6160 libstdc++, which adapted it from f2c's libF77, and is subject to
6163 /****************************************************************
6164 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
6166 Permission to use, copy, modify, and distribute this software
6167 and its documentation for any purpose and without fee is hereby
6168 granted, provided that the above copyright notice appear in all
6169 copies and that both that the copyright notice and this
6170 permission notice and warranty disclaimer appear in supporting
6171 documentation, and that the names of AT&T Bell Laboratories or
6172 Bellcore or any of their entities not be used in advertising or
6173 publicity pertaining to distribution of the software without
6174 specific, written prior permission.
6176 AT&T and Bellcore disclaim all warranties with regard to this
6177 software, including all implied warranties of merchantability
6178 and fitness. In no event shall AT&T or Bellcore be liable for
6179 any special, indirect or consequential damages or any damages
6180 whatsoever resulting from loss of use, data or profits, whether
6181 in an action of contract, negligence or other tortious action,
6182 arising out of or in connection with the use or performance of
6184 ****************************************************************/
6186 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
6187 (SCM x
, SCM y
, SCM rest
),
6188 "Divide the first argument by the product of the remaining\n"
6189 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
6191 #define FUNC_NAME s_scm_i_divide
6193 while (!scm_is_null (rest
))
6194 { x
= scm_divide (x
, y
);
6196 rest
= scm_cdr (rest
);
6198 return scm_divide (x
, y
);
6202 #define s_divide s_scm_i_divide
6203 #define g_divide g_scm_i_divide
6206 do_divide (SCM x
, SCM y
, int inexact
)
6207 #define FUNC_NAME s_divide
6211 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
6214 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
6215 else if (SCM_I_INUMP (x
))
6217 scm_t_inum xx
= SCM_I_INUM (x
);
6218 if (xx
== 1 || xx
== -1)
6220 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6222 scm_num_overflow (s_divide
);
6227 return scm_from_double (1.0 / (double) xx
);
6228 else return scm_i_make_ratio (SCM_INUM1
, x
);
6231 else if (SCM_BIGP (x
))
6234 return scm_from_double (1.0 / scm_i_big2dbl (x
));
6235 else return scm_i_make_ratio (SCM_INUM1
, x
);
6237 else if (SCM_REALP (x
))
6239 double xx
= SCM_REAL_VALUE (x
);
6240 #ifndef ALLOW_DIVIDE_BY_ZERO
6242 scm_num_overflow (s_divide
);
6245 return scm_from_double (1.0 / xx
);
6247 else if (SCM_COMPLEXP (x
))
6249 double r
= SCM_COMPLEX_REAL (x
);
6250 double i
= SCM_COMPLEX_IMAG (x
);
6251 if (fabs(r
) <= fabs(i
))
6254 double d
= i
* (1.0 + t
* t
);
6255 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
6260 double d
= r
* (1.0 + t
* t
);
6261 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
6264 else if (SCM_FRACTIONP (x
))
6265 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
6266 SCM_FRACTION_NUMERATOR (x
));
6268 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
6271 if (SCM_LIKELY (SCM_I_INUMP (x
)))
6273 scm_t_inum xx
= SCM_I_INUM (x
);
6274 if (SCM_LIKELY (SCM_I_INUMP (y
)))
6276 scm_t_inum yy
= SCM_I_INUM (y
);
6279 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6280 scm_num_overflow (s_divide
);
6282 return scm_from_double ((double) xx
/ (double) yy
);
6285 else if (xx
% yy
!= 0)
6288 return scm_from_double ((double) xx
/ (double) yy
);
6289 else return scm_i_make_ratio (x
, y
);
6293 scm_t_inum z
= xx
/ yy
;
6294 if (SCM_FIXABLE (z
))
6295 return SCM_I_MAKINUM (z
);
6297 return scm_i_inum2big (z
);
6300 else if (SCM_BIGP (y
))
6303 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
6304 else return scm_i_make_ratio (x
, y
);
6306 else if (SCM_REALP (y
))
6308 double yy
= SCM_REAL_VALUE (y
);
6309 #ifndef ALLOW_DIVIDE_BY_ZERO
6311 scm_num_overflow (s_divide
);
6314 return scm_from_double ((double) xx
/ yy
);
6316 else if (SCM_COMPLEXP (y
))
6319 complex_div
: /* y _must_ be a complex number */
6321 double r
= SCM_COMPLEX_REAL (y
);
6322 double i
= SCM_COMPLEX_IMAG (y
);
6323 if (fabs(r
) <= fabs(i
))
6326 double d
= i
* (1.0 + t
* t
);
6327 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
6332 double d
= r
* (1.0 + t
* t
);
6333 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
6337 else if (SCM_FRACTIONP (y
))
6338 /* a / b/c = ac / b */
6339 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
6340 SCM_FRACTION_NUMERATOR (y
));
6342 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6344 else if (SCM_BIGP (x
))
6346 if (SCM_I_INUMP (y
))
6348 scm_t_inum yy
= SCM_I_INUM (y
);
6351 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6352 scm_num_overflow (s_divide
);
6354 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
6355 scm_remember_upto_here_1 (x
);
6356 return (sgn
== 0) ? scm_nan () : scm_inf ();
6363 /* FIXME: HMM, what are the relative performance issues here?
6364 We need to test. Is it faster on average to test
6365 divisible_p, then perform whichever operation, or is it
6366 faster to perform the integer div opportunistically and
6367 switch to real if there's a remainder? For now we take the
6368 middle ground: test, then if divisible, use the faster div
6371 scm_t_inum abs_yy
= yy
< 0 ? -yy
: yy
;
6372 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
6376 SCM result
= scm_i_mkbig ();
6377 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
6378 scm_remember_upto_here_1 (x
);
6380 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
6381 return scm_i_normbig (result
);
6386 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
6387 else return scm_i_make_ratio (x
, y
);
6391 else if (SCM_BIGP (y
))
6396 /* It's easily possible for the ratio x/y to fit a double
6397 but one or both x and y be too big to fit a double,
6398 hence the use of mpq_get_d rather than converting and
6401 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
6402 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
6403 return scm_from_double (mpq_get_d (q
));
6407 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
6411 SCM result
= scm_i_mkbig ();
6412 mpz_divexact (SCM_I_BIG_MPZ (result
),
6415 scm_remember_upto_here_2 (x
, y
);
6416 return scm_i_normbig (result
);
6419 return scm_i_make_ratio (x
, y
);
6422 else if (SCM_REALP (y
))
6424 double yy
= SCM_REAL_VALUE (y
);
6425 #ifndef ALLOW_DIVIDE_BY_ZERO
6427 scm_num_overflow (s_divide
);
6430 return scm_from_double (scm_i_big2dbl (x
) / yy
);
6432 else if (SCM_COMPLEXP (y
))
6434 a
= scm_i_big2dbl (x
);
6437 else if (SCM_FRACTIONP (y
))
6438 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
6439 SCM_FRACTION_NUMERATOR (y
));
6441 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6443 else if (SCM_REALP (x
))
6445 double rx
= SCM_REAL_VALUE (x
);
6446 if (SCM_I_INUMP (y
))
6448 scm_t_inum yy
= SCM_I_INUM (y
);
6449 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6451 scm_num_overflow (s_divide
);
6454 return scm_from_double (rx
/ (double) yy
);
6456 else if (SCM_BIGP (y
))
6458 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6459 scm_remember_upto_here_1 (y
);
6460 return scm_from_double (rx
/ dby
);
6462 else if (SCM_REALP (y
))
6464 double yy
= SCM_REAL_VALUE (y
);
6465 #ifndef ALLOW_DIVIDE_BY_ZERO
6467 scm_num_overflow (s_divide
);
6470 return scm_from_double (rx
/ yy
);
6472 else if (SCM_COMPLEXP (y
))
6477 else if (SCM_FRACTIONP (y
))
6478 return scm_from_double (rx
/ scm_i_fraction2double (y
));
6480 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6482 else if (SCM_COMPLEXP (x
))
6484 double rx
= SCM_COMPLEX_REAL (x
);
6485 double ix
= SCM_COMPLEX_IMAG (x
);
6486 if (SCM_I_INUMP (y
))
6488 scm_t_inum yy
= SCM_I_INUM (y
);
6489 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6491 scm_num_overflow (s_divide
);
6496 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
6499 else if (SCM_BIGP (y
))
6501 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6502 scm_remember_upto_here_1 (y
);
6503 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
6505 else if (SCM_REALP (y
))
6507 double yy
= SCM_REAL_VALUE (y
);
6508 #ifndef ALLOW_DIVIDE_BY_ZERO
6510 scm_num_overflow (s_divide
);
6513 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
6515 else if (SCM_COMPLEXP (y
))
6517 double ry
= SCM_COMPLEX_REAL (y
);
6518 double iy
= SCM_COMPLEX_IMAG (y
);
6519 if (fabs(ry
) <= fabs(iy
))
6522 double d
= iy
* (1.0 + t
* t
);
6523 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
6528 double d
= ry
* (1.0 + t
* t
);
6529 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
6532 else if (SCM_FRACTIONP (y
))
6534 double yy
= scm_i_fraction2double (y
);
6535 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
6538 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6540 else if (SCM_FRACTIONP (x
))
6542 if (SCM_I_INUMP (y
))
6544 scm_t_inum yy
= SCM_I_INUM (y
);
6545 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6547 scm_num_overflow (s_divide
);
6550 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
6551 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
6553 else if (SCM_BIGP (y
))
6555 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
6556 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
6558 else if (SCM_REALP (y
))
6560 double yy
= SCM_REAL_VALUE (y
);
6561 #ifndef ALLOW_DIVIDE_BY_ZERO
6563 scm_num_overflow (s_divide
);
6566 return scm_from_double (scm_i_fraction2double (x
) / yy
);
6568 else if (SCM_COMPLEXP (y
))
6570 a
= scm_i_fraction2double (x
);
6573 else if (SCM_FRACTIONP (y
))
6574 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
6575 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
6577 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6580 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
6584 scm_divide (SCM x
, SCM y
)
6586 return do_divide (x
, y
, 0);
6589 static SCM
scm_divide2real (SCM x
, SCM y
)
6591 return do_divide (x
, y
, 1);
6597 scm_c_truncate (double x
)
6608 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
6609 half-way case (ie. when x is an integer plus 0.5) going upwards.
6610 Then half-way cases are identified and adjusted down if the
6611 round-upwards didn't give the desired even integer.
6613 "plus_half == result" identifies a half-way case. If plus_half, which is
6614 x + 0.5, is an integer then x must be an integer plus 0.5.
6616 An odd "result" value is identified with result/2 != floor(result/2).
6617 This is done with plus_half, since that value is ready for use sooner in
6618 a pipelined cpu, and we're already requiring plus_half == result.
6620 Note however that we need to be careful when x is big and already an
6621 integer. In that case "x+0.5" may round to an adjacent integer, causing
6622 us to return such a value, incorrectly. For instance if the hardware is
6623 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
6624 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
6625 returned. Or if the hardware is in round-upwards mode, then other bigger
6626 values like say x == 2^128 will see x+0.5 rounding up to the next higher
6627 representable value, 2^128+2^76 (or whatever), again incorrect.
6629 These bad roundings of x+0.5 are avoided by testing at the start whether
6630 x is already an integer. If it is then clearly that's the desired result
6631 already. And if it's not then the exponent must be small enough to allow
6632 an 0.5 to be represented, and hence added without a bad rounding. */
6635 scm_c_round (double x
)
6637 double plus_half
, result
;
6642 plus_half
= x
+ 0.5;
6643 result
= floor (plus_half
);
6644 /* Adjust so that the rounding is towards even. */
6645 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
6650 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
6652 "Round the number @var{x} towards zero.")
6653 #define FUNC_NAME s_scm_truncate_number
6655 if (scm_is_false (scm_negative_p (x
)))
6656 return scm_floor (x
);
6658 return scm_ceiling (x
);
6662 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
6664 "Round the number @var{x} towards the nearest integer. "
6665 "When it is exactly halfway between two integers, "
6666 "round towards the even one.")
6667 #define FUNC_NAME s_scm_round_number
6669 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6671 else if (SCM_REALP (x
))
6672 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
6675 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
6676 single quotient+remainder division then examining to see which way
6677 the rounding should go. */
6678 SCM plus_half
= scm_sum (x
, exactly_one_half
);
6679 SCM result
= scm_floor (plus_half
);
6680 /* Adjust so that the rounding is towards even. */
6681 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
6682 && scm_is_true (scm_odd_p (result
)))
6683 return scm_difference (result
, SCM_INUM1
);
6690 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
6692 "Round the number @var{x} towards minus infinity.")
6693 #define FUNC_NAME s_scm_floor
6695 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6697 else if (SCM_REALP (x
))
6698 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
6699 else if (SCM_FRACTIONP (x
))
6701 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
6702 SCM_FRACTION_DENOMINATOR (x
));
6703 if (scm_is_false (scm_negative_p (x
)))
6705 /* For positive x, rounding towards zero is correct. */
6710 /* For negative x, we need to return q-1 unless x is an
6711 integer. But fractions are never integer, per our
6713 return scm_difference (q
, SCM_INUM1
);
6717 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
6721 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
6723 "Round the number @var{x} towards infinity.")
6724 #define FUNC_NAME s_scm_ceiling
6726 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6728 else if (SCM_REALP (x
))
6729 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
6730 else if (SCM_FRACTIONP (x
))
6732 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
6733 SCM_FRACTION_DENOMINATOR (x
));
6734 if (scm_is_false (scm_positive_p (x
)))
6736 /* For negative x, rounding towards zero is correct. */
6741 /* For positive x, we need to return q+1 unless x is an
6742 integer. But fractions are never integer, per our
6744 return scm_sum (q
, SCM_INUM1
);
6748 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
6752 SCM_PRIMITIVE_GENERIC (scm_expt
, "expt", 2, 0, 0,
6754 "Return @var{x} raised to the power of @var{y}.")
6755 #define FUNC_NAME s_scm_expt
6757 if (scm_is_integer (y
))
6759 if (scm_is_true (scm_exact_p (y
)))
6760 return scm_integer_expt (x
, y
);
6763 /* Here we handle the case where the exponent is an inexact
6764 integer. We make the exponent exact in order to use
6765 scm_integer_expt, and thus avoid the spurious imaginary
6766 parts that may result from round-off errors in the general
6767 e^(y log x) method below (for example when squaring a large
6768 negative number). In this case, we must return an inexact
6769 result for correctness. We also make the base inexact so
6770 that scm_integer_expt will use fast inexact arithmetic
6771 internally. Note that making the base inexact is not
6772 sufficient to guarantee an inexact result, because
6773 scm_integer_expt will return an exact 1 when the exponent
6774 is 0, even if the base is inexact. */
6775 return scm_exact_to_inexact
6776 (scm_integer_expt (scm_exact_to_inexact (x
),
6777 scm_inexact_to_exact (y
)));
6780 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
6782 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
6784 else if (scm_is_complex (x
) && scm_is_complex (y
))
6785 return scm_exp (scm_product (scm_log (x
), y
));
6786 else if (scm_is_complex (x
))
6787 SCM_WTA_DISPATCH_2 (g_scm_expt
, x
, y
, SCM_ARG2
, s_scm_expt
);
6789 SCM_WTA_DISPATCH_2 (g_scm_expt
, x
, y
, SCM_ARG1
, s_scm_expt
);
6793 /* sin/cos/tan/asin/acos/atan
6794 sinh/cosh/tanh/asinh/acosh/atanh
6795 Derived from "Transcen.scm", Complex trancendental functions for SCM.
6796 Written by Jerry D. Hedden, (C) FSF.
6797 See the file `COPYING' for terms applying to this program. */
6799 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
6801 "Compute the sine of @var{z}.")
6802 #define FUNC_NAME s_scm_sin
6804 if (scm_is_real (z
))
6805 return scm_from_double (sin (scm_to_double (z
)));
6806 else if (SCM_COMPLEXP (z
))
6808 x
= SCM_COMPLEX_REAL (z
);
6809 y
= SCM_COMPLEX_IMAG (z
);
6810 return scm_c_make_rectangular (sin (x
) * cosh (y
),
6811 cos (x
) * sinh (y
));
6814 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
6818 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
6820 "Compute the cosine of @var{z}.")
6821 #define FUNC_NAME s_scm_cos
6823 if (scm_is_real (z
))
6824 return scm_from_double (cos (scm_to_double (z
)));
6825 else if (SCM_COMPLEXP (z
))
6827 x
= SCM_COMPLEX_REAL (z
);
6828 y
= SCM_COMPLEX_IMAG (z
);
6829 return scm_c_make_rectangular (cos (x
) * cosh (y
),
6830 -sin (x
) * sinh (y
));
6833 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
6837 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
6839 "Compute the tangent of @var{z}.")
6840 #define FUNC_NAME s_scm_tan
6842 if (scm_is_real (z
))
6843 return scm_from_double (tan (scm_to_double (z
)));
6844 else if (SCM_COMPLEXP (z
))
6846 x
= 2.0 * SCM_COMPLEX_REAL (z
);
6847 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
6848 w
= cos (x
) + cosh (y
);
6849 #ifndef ALLOW_DIVIDE_BY_ZERO
6851 scm_num_overflow (s_scm_tan
);
6853 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
6856 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
6860 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
6862 "Compute the hyperbolic sine of @var{z}.")
6863 #define FUNC_NAME s_scm_sinh
6865 if (scm_is_real (z
))
6866 return scm_from_double (sinh (scm_to_double (z
)));
6867 else if (SCM_COMPLEXP (z
))
6869 x
= SCM_COMPLEX_REAL (z
);
6870 y
= SCM_COMPLEX_IMAG (z
);
6871 return scm_c_make_rectangular (sinh (x
) * cos (y
),
6872 cosh (x
) * sin (y
));
6875 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
6879 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
6881 "Compute the hyperbolic cosine of @var{z}.")
6882 #define FUNC_NAME s_scm_cosh
6884 if (scm_is_real (z
))
6885 return scm_from_double (cosh (scm_to_double (z
)));
6886 else if (SCM_COMPLEXP (z
))
6888 x
= SCM_COMPLEX_REAL (z
);
6889 y
= SCM_COMPLEX_IMAG (z
);
6890 return scm_c_make_rectangular (cosh (x
) * cos (y
),
6891 sinh (x
) * sin (y
));
6894 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
6898 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
6900 "Compute the hyperbolic tangent of @var{z}.")
6901 #define FUNC_NAME s_scm_tanh
6903 if (scm_is_real (z
))
6904 return scm_from_double (tanh (scm_to_double (z
)));
6905 else if (SCM_COMPLEXP (z
))
6907 x
= 2.0 * SCM_COMPLEX_REAL (z
);
6908 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
6909 w
= cosh (x
) + cos (y
);
6910 #ifndef ALLOW_DIVIDE_BY_ZERO
6912 scm_num_overflow (s_scm_tanh
);
6914 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
6917 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
6921 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
6923 "Compute the arc sine of @var{z}.")
6924 #define FUNC_NAME s_scm_asin
6926 if (scm_is_real (z
))
6928 double w
= scm_to_double (z
);
6929 if (w
>= -1.0 && w
<= 1.0)
6930 return scm_from_double (asin (w
));
6932 return scm_product (scm_c_make_rectangular (0, -1),
6933 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
6935 else if (SCM_COMPLEXP (z
))
6937 x
= SCM_COMPLEX_REAL (z
);
6938 y
= SCM_COMPLEX_IMAG (z
);
6939 return scm_product (scm_c_make_rectangular (0, -1),
6940 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
6943 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
6947 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
6949 "Compute the arc cosine of @var{z}.")
6950 #define FUNC_NAME s_scm_acos
6952 if (scm_is_real (z
))
6954 double w
= scm_to_double (z
);
6955 if (w
>= -1.0 && w
<= 1.0)
6956 return scm_from_double (acos (w
));
6958 return scm_sum (scm_from_double (acos (0.0)),
6959 scm_product (scm_c_make_rectangular (0, 1),
6960 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
6962 else if (SCM_COMPLEXP (z
))
6964 x
= SCM_COMPLEX_REAL (z
);
6965 y
= SCM_COMPLEX_IMAG (z
);
6966 return scm_sum (scm_from_double (acos (0.0)),
6967 scm_product (scm_c_make_rectangular (0, 1),
6968 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
6971 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
6975 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
6977 "With one argument, compute the arc tangent of @var{z}.\n"
6978 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
6979 "using the sign of @var{z} and @var{y} to determine the quadrant.")
6980 #define FUNC_NAME s_scm_atan
6984 if (scm_is_real (z
))
6985 return scm_from_double (atan (scm_to_double (z
)));
6986 else if (SCM_COMPLEXP (z
))
6989 v
= SCM_COMPLEX_REAL (z
);
6990 w
= SCM_COMPLEX_IMAG (z
);
6991 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
6992 scm_c_make_rectangular (v
, w
+ 1.0))),
6993 scm_c_make_rectangular (0, 2));
6996 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
6998 else if (scm_is_real (z
))
7000 if (scm_is_real (y
))
7001 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
7003 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
7006 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
7010 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
7012 "Compute the inverse hyperbolic sine of @var{z}.")
7013 #define FUNC_NAME s_scm_sys_asinh
7015 if (scm_is_real (z
))
7016 return scm_from_double (asinh (scm_to_double (z
)));
7017 else if (scm_is_number (z
))
7018 return scm_log (scm_sum (z
,
7019 scm_sqrt (scm_sum (scm_product (z
, z
),
7022 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
7026 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
7028 "Compute the inverse hyperbolic cosine of @var{z}.")
7029 #define FUNC_NAME s_scm_sys_acosh
7031 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
7032 return scm_from_double (acosh (scm_to_double (z
)));
7033 else if (scm_is_number (z
))
7034 return scm_log (scm_sum (z
,
7035 scm_sqrt (scm_difference (scm_product (z
, z
),
7038 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
7042 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
7044 "Compute the inverse hyperbolic tangent of @var{z}.")
7045 #define FUNC_NAME s_scm_sys_atanh
7047 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
7048 return scm_from_double (atanh (scm_to_double (z
)));
7049 else if (scm_is_number (z
))
7050 return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1
, z
),
7051 scm_difference (SCM_INUM1
, z
))),
7054 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
7059 scm_c_make_rectangular (double re
, double im
)
7062 return scm_from_double (re
);
7067 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
7069 SCM_SET_CELL_TYPE (z
, scm_tc16_complex
);
7070 SCM_COMPLEX_REAL (z
) = re
;
7071 SCM_COMPLEX_IMAG (z
) = im
;
7076 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
7077 (SCM real_part
, SCM imaginary_part
),
7078 "Return a complex number constructed of the given @var{real-part} "
7079 "and @var{imaginary-part} parts.")
7080 #define FUNC_NAME s_scm_make_rectangular
7082 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
7083 SCM_ARG1
, FUNC_NAME
, "real");
7084 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
7085 SCM_ARG2
, FUNC_NAME
, "real");
7086 return scm_c_make_rectangular (scm_to_double (real_part
),
7087 scm_to_double (imaginary_part
));
7092 scm_c_make_polar (double mag
, double ang
)
7096 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
7097 use it on Glibc-based systems that have it (it's a GNU extension). See
7098 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
7100 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
7101 sincos (ang
, &s
, &c
);
7106 return scm_c_make_rectangular (mag
* c
, mag
* s
);
7109 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
7111 "Return the complex number @var{x} * e^(i * @var{y}).")
7112 #define FUNC_NAME s_scm_make_polar
7114 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
7115 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
7116 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
7121 SCM_PRIMITIVE_GENERIC (scm_real_part
, "real-part", 1, 0, 0,
7123 "Return the real part of the number @var{z}.")
7124 #define FUNC_NAME s_scm_real_part
7126 if (SCM_COMPLEXP (z
))
7127 return scm_from_double (SCM_COMPLEX_REAL (z
));
7128 else if (SCM_I_INUMP (z
) || SCM_BIGP (z
) || SCM_REALP (z
) || SCM_FRACTIONP (z
))
7131 SCM_WTA_DISPATCH_1 (g_scm_real_part
, z
, SCM_ARG1
, s_scm_real_part
);
7136 SCM_PRIMITIVE_GENERIC (scm_imag_part
, "imag-part", 1, 0, 0,
7138 "Return the imaginary part of the number @var{z}.")
7139 #define FUNC_NAME s_scm_imag_part
7141 if (SCM_COMPLEXP (z
))
7142 return scm_from_double (SCM_COMPLEX_IMAG (z
));
7143 else if (SCM_REALP (z
))
7145 else if (SCM_I_INUMP (z
) || SCM_BIGP (z
) || SCM_FRACTIONP (z
))
7148 SCM_WTA_DISPATCH_1 (g_scm_imag_part
, z
, SCM_ARG1
, s_scm_imag_part
);
7152 SCM_PRIMITIVE_GENERIC (scm_numerator
, "numerator", 1, 0, 0,
7154 "Return the numerator of the number @var{z}.")
7155 #define FUNC_NAME s_scm_numerator
7157 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7159 else if (SCM_FRACTIONP (z
))
7160 return SCM_FRACTION_NUMERATOR (z
);
7161 else if (SCM_REALP (z
))
7162 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
7164 SCM_WTA_DISPATCH_1 (g_scm_numerator
, z
, SCM_ARG1
, s_scm_numerator
);
7169 SCM_PRIMITIVE_GENERIC (scm_denominator
, "denominator", 1, 0, 0,
7171 "Return the denominator of the number @var{z}.")
7172 #define FUNC_NAME s_scm_denominator
7174 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7176 else if (SCM_FRACTIONP (z
))
7177 return SCM_FRACTION_DENOMINATOR (z
);
7178 else if (SCM_REALP (z
))
7179 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
7181 SCM_WTA_DISPATCH_1 (g_scm_denominator
, z
, SCM_ARG1
, s_scm_denominator
);
7186 SCM_PRIMITIVE_GENERIC (scm_magnitude
, "magnitude", 1, 0, 0,
7188 "Return the magnitude of the number @var{z}. This is the same as\n"
7189 "@code{abs} for real arguments, but also allows complex numbers.")
7190 #define FUNC_NAME s_scm_magnitude
7192 if (SCM_I_INUMP (z
))
7194 scm_t_inum zz
= SCM_I_INUM (z
);
7197 else if (SCM_POSFIXABLE (-zz
))
7198 return SCM_I_MAKINUM (-zz
);
7200 return scm_i_inum2big (-zz
);
7202 else if (SCM_BIGP (z
))
7204 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
7205 scm_remember_upto_here_1 (z
);
7207 return scm_i_clonebig (z
, 0);
7211 else if (SCM_REALP (z
))
7212 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
7213 else if (SCM_COMPLEXP (z
))
7214 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
7215 else if (SCM_FRACTIONP (z
))
7217 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
7219 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
7220 SCM_FRACTION_DENOMINATOR (z
));
7223 SCM_WTA_DISPATCH_1 (g_scm_magnitude
, z
, SCM_ARG1
, s_scm_magnitude
);
7228 SCM_PRIMITIVE_GENERIC (scm_angle
, "angle", 1, 0, 0,
7230 "Return the angle of the complex number @var{z}.")
7231 #define FUNC_NAME s_scm_angle
7233 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
7234 flo0 to save allocating a new flonum with scm_from_double each time.
7235 But if atan2 follows the floating point rounding mode, then the value
7236 is not a constant. Maybe it'd be close enough though. */
7237 if (SCM_I_INUMP (z
))
7239 if (SCM_I_INUM (z
) >= 0)
7242 return scm_from_double (atan2 (0.0, -1.0));
7244 else if (SCM_BIGP (z
))
7246 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
7247 scm_remember_upto_here_1 (z
);
7249 return scm_from_double (atan2 (0.0, -1.0));
7253 else if (SCM_REALP (z
))
7255 if (SCM_REAL_VALUE (z
) >= 0)
7258 return scm_from_double (atan2 (0.0, -1.0));
7260 else if (SCM_COMPLEXP (z
))
7261 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
7262 else if (SCM_FRACTIONP (z
))
7264 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
7266 else return scm_from_double (atan2 (0.0, -1.0));
7269 SCM_WTA_DISPATCH_1 (g_scm_angle
, z
, SCM_ARG1
, s_scm_angle
);
7274 SCM_PRIMITIVE_GENERIC (scm_exact_to_inexact
, "exact->inexact", 1, 0, 0,
7276 "Convert the number @var{z} to its inexact representation.\n")
7277 #define FUNC_NAME s_scm_exact_to_inexact
7279 if (SCM_I_INUMP (z
))
7280 return scm_from_double ((double) SCM_I_INUM (z
));
7281 else if (SCM_BIGP (z
))
7282 return scm_from_double (scm_i_big2dbl (z
));
7283 else if (SCM_FRACTIONP (z
))
7284 return scm_from_double (scm_i_fraction2double (z
));
7285 else if (SCM_INEXACTP (z
))
7288 SCM_WTA_DISPATCH_1 (g_scm_exact_to_inexact
, z
, 1, s_scm_exact_to_inexact
);
7293 SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
7295 "Return an exact number that is numerically closest to @var{z}.")
7296 #define FUNC_NAME s_scm_inexact_to_exact
7298 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7300 else if (SCM_REALP (z
))
7302 if (!DOUBLE_IS_FINITE (SCM_REAL_VALUE (z
)))
7303 SCM_OUT_OF_RANGE (1, z
);
7310 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
7311 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
7312 scm_i_mpz2num (mpq_denref (frac
)));
7314 /* When scm_i_make_ratio throws, we leak the memory allocated
7321 else if (SCM_FRACTIONP (z
))
7324 SCM_WTA_DISPATCH_1 (g_scm_inexact_to_exact
, z
, 1, s_scm_inexact_to_exact
);
7328 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
7330 "Returns the @emph{simplest} rational number differing\n"
7331 "from @var{x} by no more than @var{eps}.\n"
7333 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
7334 "exact result when both its arguments are exact. Thus, you might need\n"
7335 "to use @code{inexact->exact} on the arguments.\n"
7338 "(rationalize (inexact->exact 1.2) 1/100)\n"
7341 #define FUNC_NAME s_scm_rationalize
7343 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
7344 SCM_ASSERT_TYPE (scm_is_real (eps
), eps
, SCM_ARG2
, FUNC_NAME
, "real");
7345 eps
= scm_abs (eps
);
7346 if (scm_is_false (scm_positive_p (eps
)))
7348 /* eps is either zero or a NaN */
7349 if (scm_is_true (scm_nan_p (eps
)))
7351 else if (SCM_INEXACTP (eps
))
7352 return scm_exact_to_inexact (x
);
7356 else if (scm_is_false (scm_finite_p (eps
)))
7358 if (scm_is_true (scm_finite_p (x
)))
7363 else if (scm_is_false (scm_finite_p (x
))) /* checks for both inf and nan */
7365 else if (scm_is_false (scm_less_p (scm_floor (scm_sum (x
, eps
)),
7366 scm_ceiling (scm_difference (x
, eps
)))))
7368 /* There's an integer within range; we want the one closest to zero */
7369 if (scm_is_false (scm_less_p (eps
, scm_abs (x
))))
7371 /* zero is within range */
7372 if (SCM_INEXACTP (x
) || SCM_INEXACTP (eps
))
7377 else if (scm_is_true (scm_positive_p (x
)))
7378 return scm_ceiling (scm_difference (x
, eps
));
7380 return scm_floor (scm_sum (x
, eps
));
7384 /* Use continued fractions to find closest ratio. All
7385 arithmetic is done with exact numbers.
7388 SCM ex
= scm_inexact_to_exact (x
);
7389 SCM int_part
= scm_floor (ex
);
7391 SCM a1
= SCM_INUM0
, a2
= SCM_INUM1
, a
= SCM_INUM0
;
7392 SCM b1
= SCM_INUM1
, b2
= SCM_INUM0
, b
= SCM_INUM0
;
7396 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
7397 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
7399 /* We stop after a million iterations just to be absolutely sure
7400 that we don't go into an infinite loop. The process normally
7401 converges after less than a dozen iterations.
7404 while (++i
< 1000000)
7406 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
7407 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
7408 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
7410 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
7411 eps
))) /* abs(x-a/b) <= eps */
7413 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
7414 if (SCM_INEXACTP (x
) || SCM_INEXACTP (eps
))
7415 return scm_exact_to_inexact (res
);
7419 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
7421 tt
= scm_floor (rx
); /* tt = floor (rx) */
7427 scm_num_overflow (s_scm_rationalize
);
7432 /* conversion functions */
7435 scm_is_integer (SCM val
)
7437 return scm_is_true (scm_integer_p (val
));
7441 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
7443 if (SCM_I_INUMP (val
))
7445 scm_t_signed_bits n
= SCM_I_INUM (val
);
7446 return n
>= min
&& n
<= max
;
7448 else if (SCM_BIGP (val
))
7450 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
7452 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
7454 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
7456 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
7457 return n
>= min
&& n
<= max
;
7467 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
7468 > CHAR_BIT
*sizeof (scm_t_uintmax
))
7471 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
7472 SCM_I_BIG_MPZ (val
));
7474 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
7486 return n
>= min
&& n
<= max
;
7494 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
7496 if (SCM_I_INUMP (val
))
7498 scm_t_signed_bits n
= SCM_I_INUM (val
);
7499 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
7501 else if (SCM_BIGP (val
))
7503 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
7505 else if (max
<= ULONG_MAX
)
7507 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
7509 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
7510 return n
>= min
&& n
<= max
;
7520 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
7523 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
7524 > CHAR_BIT
*sizeof (scm_t_uintmax
))
7527 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
7528 SCM_I_BIG_MPZ (val
));
7530 return n
>= min
&& n
<= max
;
7538 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
7540 scm_error (scm_out_of_range_key
,
7542 "Value out of range ~S to ~S: ~S",
7543 scm_list_3 (min
, max
, bad_val
),
7544 scm_list_1 (bad_val
));
7547 #define TYPE scm_t_intmax
7548 #define TYPE_MIN min
7549 #define TYPE_MAX max
7550 #define SIZEOF_TYPE 0
7551 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
7552 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
7553 #include "libguile/conv-integer.i.c"
7555 #define TYPE scm_t_uintmax
7556 #define TYPE_MIN min
7557 #define TYPE_MAX max
7558 #define SIZEOF_TYPE 0
7559 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
7560 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
7561 #include "libguile/conv-uinteger.i.c"
7563 #define TYPE scm_t_int8
7564 #define TYPE_MIN SCM_T_INT8_MIN
7565 #define TYPE_MAX SCM_T_INT8_MAX
7566 #define SIZEOF_TYPE 1
7567 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
7568 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
7569 #include "libguile/conv-integer.i.c"
7571 #define TYPE scm_t_uint8
7573 #define TYPE_MAX SCM_T_UINT8_MAX
7574 #define SIZEOF_TYPE 1
7575 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
7576 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
7577 #include "libguile/conv-uinteger.i.c"
7579 #define TYPE scm_t_int16
7580 #define TYPE_MIN SCM_T_INT16_MIN
7581 #define TYPE_MAX SCM_T_INT16_MAX
7582 #define SIZEOF_TYPE 2
7583 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
7584 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
7585 #include "libguile/conv-integer.i.c"
7587 #define TYPE scm_t_uint16
7589 #define TYPE_MAX SCM_T_UINT16_MAX
7590 #define SIZEOF_TYPE 2
7591 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
7592 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
7593 #include "libguile/conv-uinteger.i.c"
7595 #define TYPE scm_t_int32
7596 #define TYPE_MIN SCM_T_INT32_MIN
7597 #define TYPE_MAX SCM_T_INT32_MAX
7598 #define SIZEOF_TYPE 4
7599 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
7600 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
7601 #include "libguile/conv-integer.i.c"
7603 #define TYPE scm_t_uint32
7605 #define TYPE_MAX SCM_T_UINT32_MAX
7606 #define SIZEOF_TYPE 4
7607 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
7608 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
7609 #include "libguile/conv-uinteger.i.c"
7611 #define TYPE scm_t_wchar
7612 #define TYPE_MIN (scm_t_int32)-1
7613 #define TYPE_MAX (scm_t_int32)0x10ffff
7614 #define SIZEOF_TYPE 4
7615 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
7616 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
7617 #include "libguile/conv-integer.i.c"
7619 #define TYPE scm_t_int64
7620 #define TYPE_MIN SCM_T_INT64_MIN
7621 #define TYPE_MAX SCM_T_INT64_MAX
7622 #define SIZEOF_TYPE 8
7623 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
7624 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
7625 #include "libguile/conv-integer.i.c"
7627 #define TYPE scm_t_uint64
7629 #define TYPE_MAX SCM_T_UINT64_MAX
7630 #define SIZEOF_TYPE 8
7631 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
7632 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
7633 #include "libguile/conv-uinteger.i.c"
7636 scm_to_mpz (SCM val
, mpz_t rop
)
7638 if (SCM_I_INUMP (val
))
7639 mpz_set_si (rop
, SCM_I_INUM (val
));
7640 else if (SCM_BIGP (val
))
7641 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
7643 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
7647 scm_from_mpz (mpz_t val
)
7649 return scm_i_mpz2num (val
);
7653 scm_is_real (SCM val
)
7655 return scm_is_true (scm_real_p (val
));
7659 scm_is_rational (SCM val
)
7661 return scm_is_true (scm_rational_p (val
));
7665 scm_to_double (SCM val
)
7667 if (SCM_I_INUMP (val
))
7668 return SCM_I_INUM (val
);
7669 else if (SCM_BIGP (val
))
7670 return scm_i_big2dbl (val
);
7671 else if (SCM_FRACTIONP (val
))
7672 return scm_i_fraction2double (val
);
7673 else if (SCM_REALP (val
))
7674 return SCM_REAL_VALUE (val
);
7676 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
7680 scm_from_double (double val
)
7684 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
7686 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
7687 SCM_REAL_VALUE (z
) = val
;
7692 #if SCM_ENABLE_DEPRECATED == 1
7695 scm_num2float (SCM num
, unsigned long pos
, const char *s_caller
)
7697 scm_c_issue_deprecation_warning
7698 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
7702 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
7706 scm_out_of_range (NULL
, num
);
7709 return scm_to_double (num
);
7713 scm_num2double (SCM num
, unsigned long pos
, const char *s_caller
)
7715 scm_c_issue_deprecation_warning
7716 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
7720 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
7724 scm_out_of_range (NULL
, num
);
7727 return scm_to_double (num
);
7733 scm_is_complex (SCM val
)
7735 return scm_is_true (scm_complex_p (val
));
7739 scm_c_real_part (SCM z
)
7741 if (SCM_COMPLEXP (z
))
7742 return SCM_COMPLEX_REAL (z
);
7745 /* Use the scm_real_part to get proper error checking and
7748 return scm_to_double (scm_real_part (z
));
7753 scm_c_imag_part (SCM z
)
7755 if (SCM_COMPLEXP (z
))
7756 return SCM_COMPLEX_IMAG (z
);
7759 /* Use the scm_imag_part to get proper error checking and
7760 dispatching. The result will almost always be 0.0, but not
7763 return scm_to_double (scm_imag_part (z
));
7768 scm_c_magnitude (SCM z
)
7770 return scm_to_double (scm_magnitude (z
));
7776 return scm_to_double (scm_angle (z
));
7780 scm_is_number (SCM z
)
7782 return scm_is_true (scm_number_p (z
));
7786 /* In the following functions we dispatch to the real-arg funcs like log()
7787 when we know the arg is real, instead of just handing everything to
7788 clog() for instance. This is in case clog() doesn't optimize for a
7789 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
7790 well use it to go straight to the applicable C func. */
7792 SCM_PRIMITIVE_GENERIC (scm_log
, "log", 1, 0, 0,
7794 "Return the natural logarithm of @var{z}.")
7795 #define FUNC_NAME s_scm_log
7797 if (SCM_COMPLEXP (z
))
7799 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
7800 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
7802 double re
= SCM_COMPLEX_REAL (z
);
7803 double im
= SCM_COMPLEX_IMAG (z
);
7804 return scm_c_make_rectangular (log (hypot (re
, im
)),
7808 else if (SCM_NUMBERP (z
))
7810 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
7811 although the value itself overflows. */
7812 double re
= scm_to_double (z
);
7813 double l
= log (fabs (re
));
7815 return scm_from_double (l
);
7817 return scm_c_make_rectangular (l
, M_PI
);
7820 SCM_WTA_DISPATCH_1 (g_scm_log
, z
, 1, s_scm_log
);
7825 SCM_PRIMITIVE_GENERIC (scm_log10
, "log10", 1, 0, 0,
7827 "Return the base 10 logarithm of @var{z}.")
7828 #define FUNC_NAME s_scm_log10
7830 if (SCM_COMPLEXP (z
))
7832 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
7833 clog() and a multiply by M_LOG10E, rather than the fallback
7834 log10+hypot+atan2.) */
7835 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
7836 && defined SCM_COMPLEX_VALUE
7837 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
7839 double re
= SCM_COMPLEX_REAL (z
);
7840 double im
= SCM_COMPLEX_IMAG (z
);
7841 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
7842 M_LOG10E
* atan2 (im
, re
));
7845 else if (SCM_NUMBERP (z
))
7847 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
7848 although the value itself overflows. */
7849 double re
= scm_to_double (z
);
7850 double l
= log10 (fabs (re
));
7852 return scm_from_double (l
);
7854 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
7857 SCM_WTA_DISPATCH_1 (g_scm_log10
, z
, 1, s_scm_log10
);
7862 SCM_PRIMITIVE_GENERIC (scm_exp
, "exp", 1, 0, 0,
7864 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
7865 "base of natural logarithms (2.71828@dots{}).")
7866 #define FUNC_NAME s_scm_exp
7868 if (SCM_COMPLEXP (z
))
7870 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
7871 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
7873 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
7874 SCM_COMPLEX_IMAG (z
));
7877 else if (SCM_NUMBERP (z
))
7879 /* When z is a negative bignum the conversion to double overflows,
7880 giving -infinity, but that's ok, the exp is still 0.0. */
7881 return scm_from_double (exp (scm_to_double (z
)));
7884 SCM_WTA_DISPATCH_1 (g_scm_exp
, z
, 1, s_scm_exp
);
7889 SCM_PRIMITIVE_GENERIC (scm_sqrt
, "sqrt", 1, 0, 0,
7891 "Return the square root of @var{z}. Of the two possible roots\n"
7892 "(positive and negative), the one with the a positive real part\n"
7893 "is returned, or if that's zero then a positive imaginary part.\n"
7897 "(sqrt 9.0) @result{} 3.0\n"
7898 "(sqrt -9.0) @result{} 0.0+3.0i\n"
7899 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
7900 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
7902 #define FUNC_NAME s_scm_sqrt
7904 if (SCM_COMPLEXP (z
))
7906 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
7907 && defined SCM_COMPLEX_VALUE
7908 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (z
)));
7910 double re
= SCM_COMPLEX_REAL (z
);
7911 double im
= SCM_COMPLEX_IMAG (z
);
7912 return scm_c_make_polar (sqrt (hypot (re
, im
)),
7913 0.5 * atan2 (im
, re
));
7916 else if (SCM_NUMBERP (z
))
7918 double xx
= scm_to_double (z
);
7920 return scm_c_make_rectangular (0.0, sqrt (-xx
));
7922 return scm_from_double (sqrt (xx
));
7925 SCM_WTA_DISPATCH_1 (g_scm_sqrt
, z
, 1, s_scm_sqrt
);
7936 mpz_init_set_si (z_negative_one
, -1);
7938 /* It may be possible to tune the performance of some algorithms by using
7939 * the following constants to avoid the creation of bignums. Please, before
7940 * using these values, remember the two rules of program optimization:
7941 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
7942 scm_c_define ("most-positive-fixnum",
7943 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
7944 scm_c_define ("most-negative-fixnum",
7945 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
7947 scm_add_feature ("complex");
7948 scm_add_feature ("inexact");
7949 flo0
= scm_from_double (0.0);
7951 /* determine floating point precision */
7952 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
7954 init_dblprec(&scm_dblprec
[i
-2],i
);
7955 init_fx_radix(fx_per_radix
[i
-2],i
);
7958 /* hard code precision for base 10 if the preprocessor tells us to... */
7959 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
7962 exactly_one_half
= scm_divide (SCM_INUM1
, SCM_I_MAKINUM (2));
7963 #include "libguile/numbers.x"