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
)));
501 SCM_PRIMITIVE_GENERIC (scm_exact_p
, "exact?", 1, 0, 0,
503 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
505 #define FUNC_NAME s_scm_exact_p
507 if (SCM_INEXACTP (x
))
509 else if (SCM_NUMBERP (x
))
512 SCM_WTA_DISPATCH_1 (g_scm_exact_p
, x
, 1, s_scm_exact_p
);
517 SCM_PRIMITIVE_GENERIC (scm_inexact_p
, "inexact?", 1, 0, 0,
519 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
521 #define FUNC_NAME s_scm_inexact_p
523 if (SCM_INEXACTP (x
))
525 else if (SCM_NUMBERP (x
))
528 SCM_WTA_DISPATCH_1 (g_scm_inexact_p
, x
, 1, s_scm_inexact_p
);
533 SCM_PRIMITIVE_GENERIC (scm_odd_p
, "odd?", 1, 0, 0,
535 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
537 #define FUNC_NAME s_scm_odd_p
541 scm_t_inum val
= SCM_I_INUM (n
);
542 return scm_from_bool ((val
& 1L) != 0);
544 else if (SCM_BIGP (n
))
546 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
547 scm_remember_upto_here_1 (n
);
548 return scm_from_bool (odd_p
);
550 else if (SCM_REALP (n
))
552 double val
= SCM_REAL_VALUE (n
);
553 if (DOUBLE_IS_FINITE (val
))
555 double rem
= fabs (fmod (val
, 2.0));
562 SCM_WTA_DISPATCH_1 (g_scm_odd_p
, n
, 1, s_scm_odd_p
);
567 SCM_PRIMITIVE_GENERIC (scm_even_p
, "even?", 1, 0, 0,
569 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
571 #define FUNC_NAME s_scm_even_p
575 scm_t_inum val
= SCM_I_INUM (n
);
576 return scm_from_bool ((val
& 1L) == 0);
578 else if (SCM_BIGP (n
))
580 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
581 scm_remember_upto_here_1 (n
);
582 return scm_from_bool (even_p
);
584 else if (SCM_REALP (n
))
586 double val
= SCM_REAL_VALUE (n
);
587 if (DOUBLE_IS_FINITE (val
))
589 double rem
= fabs (fmod (val
, 2.0));
596 SCM_WTA_DISPATCH_1 (g_scm_even_p
, n
, 1, s_scm_even_p
);
600 SCM_PRIMITIVE_GENERIC (scm_finite_p
, "finite?", 1, 0, 0,
602 "Return @code{#t} if the real number @var{x} is neither\n"
603 "infinite nor a NaN, @code{#f} otherwise.")
604 #define FUNC_NAME s_scm_finite_p
607 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
608 else if (scm_is_real (x
))
611 SCM_WTA_DISPATCH_1 (g_scm_finite_p
, x
, 1, s_scm_finite_p
);
615 SCM_PRIMITIVE_GENERIC (scm_inf_p
, "inf?", 1, 0, 0,
617 "Return @code{#t} if the real number @var{x} is @samp{+inf.0} or\n"
618 "@samp{-inf.0}. Otherwise return @code{#f}.")
619 #define FUNC_NAME s_scm_inf_p
622 return scm_from_bool (isinf (SCM_REAL_VALUE (x
)));
623 else if (scm_is_real (x
))
626 SCM_WTA_DISPATCH_1 (g_scm_inf_p
, x
, 1, s_scm_inf_p
);
630 SCM_PRIMITIVE_GENERIC (scm_nan_p
, "nan?", 1, 0, 0,
632 "Return @code{#t} if the real number @var{x} is a NaN,\n"
633 "or @code{#f} otherwise.")
634 #define FUNC_NAME s_scm_nan_p
637 return scm_from_bool (isnan (SCM_REAL_VALUE (x
)));
638 else if (scm_is_real (x
))
641 SCM_WTA_DISPATCH_1 (g_scm_nan_p
, x
, 1, s_scm_nan_p
);
645 /* Guile's idea of infinity. */
646 static double guile_Inf
;
648 /* Guile's idea of not a number. */
649 static double guile_NaN
;
652 guile_ieee_init (void)
654 /* Some version of gcc on some old version of Linux used to crash when
655 trying to make Inf and NaN. */
658 /* C99 INFINITY, when available.
659 FIXME: The standard allows for INFINITY to be something that overflows
660 at compile time. We ought to have a configure test to check for that
661 before trying to use it. (But in practice we believe this is not a
662 problem on any system guile is likely to target.) */
663 guile_Inf
= INFINITY
;
664 #elif defined HAVE_DINFINITY
666 extern unsigned int DINFINITY
[2];
667 guile_Inf
= (*((double *) (DINFINITY
)));
674 if (guile_Inf
== tmp
)
681 /* C99 NAN, when available */
683 #elif defined HAVE_DQNAN
686 extern unsigned int DQNAN
[2];
687 guile_NaN
= (*((double *)(DQNAN
)));
690 guile_NaN
= guile_Inf
/ guile_Inf
;
694 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
697 #define FUNC_NAME s_scm_inf
699 static int initialized
= 0;
705 return scm_from_double (guile_Inf
);
709 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
712 #define FUNC_NAME s_scm_nan
714 static int initialized
= 0;
720 return scm_from_double (guile_NaN
);
725 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
727 "Return the absolute value of @var{x}.")
728 #define FUNC_NAME s_scm_abs
732 scm_t_inum xx
= SCM_I_INUM (x
);
735 else if (SCM_POSFIXABLE (-xx
))
736 return SCM_I_MAKINUM (-xx
);
738 return scm_i_inum2big (-xx
);
740 else if (SCM_BIGP (x
))
742 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
744 return scm_i_clonebig (x
, 0);
748 else if (SCM_REALP (x
))
750 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
751 double xx
= SCM_REAL_VALUE (x
);
753 return scm_from_double (-xx
);
757 else if (SCM_FRACTIONP (x
))
759 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
761 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
762 SCM_FRACTION_DENOMINATOR (x
));
765 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
770 SCM_PRIMITIVE_GENERIC (scm_quotient
, "quotient", 2, 0, 0,
772 "Return the quotient of the numbers @var{x} and @var{y}.")
773 #define FUNC_NAME s_scm_quotient
775 if (SCM_LIKELY (SCM_I_INUMP (x
)))
777 scm_t_inum xx
= SCM_I_INUM (x
);
778 if (SCM_LIKELY (SCM_I_INUMP (y
)))
780 scm_t_inum yy
= SCM_I_INUM (y
);
781 if (SCM_UNLIKELY (yy
== 0))
782 scm_num_overflow (s_scm_quotient
);
785 scm_t_inum z
= xx
/ yy
;
786 if (SCM_LIKELY (SCM_FIXABLE (z
)))
787 return SCM_I_MAKINUM (z
);
789 return scm_i_inum2big (z
);
792 else if (SCM_BIGP (y
))
794 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
795 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
796 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
798 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
799 scm_remember_upto_here_1 (y
);
800 return SCM_I_MAKINUM (-1);
806 SCM_WTA_DISPATCH_2 (g_scm_quotient
, x
, y
, SCM_ARG2
, s_scm_quotient
);
808 else if (SCM_BIGP (x
))
810 if (SCM_LIKELY (SCM_I_INUMP (y
)))
812 scm_t_inum yy
= SCM_I_INUM (y
);
813 if (SCM_UNLIKELY (yy
== 0))
814 scm_num_overflow (s_scm_quotient
);
815 else if (SCM_UNLIKELY (yy
== 1))
819 SCM result
= scm_i_mkbig ();
822 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
825 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
828 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
829 scm_remember_upto_here_1 (x
);
830 return scm_i_normbig (result
);
833 else if (SCM_BIGP (y
))
835 SCM result
= scm_i_mkbig ();
836 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
839 scm_remember_upto_here_2 (x
, y
);
840 return scm_i_normbig (result
);
843 SCM_WTA_DISPATCH_2 (g_scm_quotient
, x
, y
, SCM_ARG2
, s_scm_quotient
);
846 SCM_WTA_DISPATCH_2 (g_scm_quotient
, x
, y
, SCM_ARG1
, s_scm_quotient
);
850 SCM_PRIMITIVE_GENERIC (scm_remainder
, "remainder", 2, 0, 0,
852 "Return the remainder of the numbers @var{x} and @var{y}.\n"
854 "(remainder 13 4) @result{} 1\n"
855 "(remainder -13 4) @result{} -1\n"
857 #define FUNC_NAME s_scm_remainder
859 if (SCM_LIKELY (SCM_I_INUMP (x
)))
861 if (SCM_LIKELY (SCM_I_INUMP (y
)))
863 scm_t_inum yy
= SCM_I_INUM (y
);
864 if (SCM_UNLIKELY (yy
== 0))
865 scm_num_overflow (s_scm_remainder
);
868 /* C99 specifies that "%" is the remainder corresponding to a
869 quotient rounded towards zero, and that's also traditional
870 for machine division, so z here should be well defined. */
871 scm_t_inum z
= SCM_I_INUM (x
) % yy
;
872 return SCM_I_MAKINUM (z
);
875 else if (SCM_BIGP (y
))
877 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
878 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
879 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
881 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
882 scm_remember_upto_here_1 (y
);
889 SCM_WTA_DISPATCH_2 (g_scm_remainder
, x
, y
, SCM_ARG2
, s_scm_remainder
);
891 else if (SCM_BIGP (x
))
893 if (SCM_LIKELY (SCM_I_INUMP (y
)))
895 scm_t_inum yy
= SCM_I_INUM (y
);
896 if (SCM_UNLIKELY (yy
== 0))
897 scm_num_overflow (s_scm_remainder
);
900 SCM result
= scm_i_mkbig ();
903 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
904 scm_remember_upto_here_1 (x
);
905 return scm_i_normbig (result
);
908 else if (SCM_BIGP (y
))
910 SCM result
= scm_i_mkbig ();
911 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
914 scm_remember_upto_here_2 (x
, y
);
915 return scm_i_normbig (result
);
918 SCM_WTA_DISPATCH_2 (g_scm_remainder
, x
, y
, SCM_ARG2
, s_scm_remainder
);
921 SCM_WTA_DISPATCH_2 (g_scm_remainder
, x
, y
, SCM_ARG1
, s_scm_remainder
);
926 SCM_PRIMITIVE_GENERIC (scm_modulo
, "modulo", 2, 0, 0,
928 "Return the modulo of the numbers @var{x} and @var{y}.\n"
930 "(modulo 13 4) @result{} 1\n"
931 "(modulo -13 4) @result{} 3\n"
933 #define FUNC_NAME s_scm_modulo
935 if (SCM_LIKELY (SCM_I_INUMP (x
)))
937 scm_t_inum xx
= SCM_I_INUM (x
);
938 if (SCM_LIKELY (SCM_I_INUMP (y
)))
940 scm_t_inum yy
= SCM_I_INUM (y
);
941 if (SCM_UNLIKELY (yy
== 0))
942 scm_num_overflow (s_scm_modulo
);
945 /* C99 specifies that "%" is the remainder corresponding to a
946 quotient rounded towards zero, and that's also traditional
947 for machine division, so z here should be well defined. */
948 scm_t_inum z
= xx
% yy
;
965 return SCM_I_MAKINUM (result
);
968 else if (SCM_BIGP (y
))
970 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
977 SCM pos_y
= scm_i_clonebig (y
, 0);
978 /* do this after the last scm_op */
979 mpz_init_set_si (z_x
, xx
);
980 result
= pos_y
; /* re-use this bignum */
981 mpz_mod (SCM_I_BIG_MPZ (result
),
983 SCM_I_BIG_MPZ (pos_y
));
984 scm_remember_upto_here_1 (pos_y
);
988 result
= scm_i_mkbig ();
989 /* do this after the last scm_op */
990 mpz_init_set_si (z_x
, xx
);
991 mpz_mod (SCM_I_BIG_MPZ (result
),
994 scm_remember_upto_here_1 (y
);
997 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
998 mpz_add (SCM_I_BIG_MPZ (result
),
1000 SCM_I_BIG_MPZ (result
));
1001 scm_remember_upto_here_1 (y
);
1002 /* and do this before the next one */
1004 return scm_i_normbig (result
);
1008 SCM_WTA_DISPATCH_2 (g_scm_modulo
, x
, y
, SCM_ARG2
, s_scm_modulo
);
1010 else if (SCM_BIGP (x
))
1012 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1014 scm_t_inum yy
= SCM_I_INUM (y
);
1015 if (SCM_UNLIKELY (yy
== 0))
1016 scm_num_overflow (s_scm_modulo
);
1019 SCM result
= scm_i_mkbig ();
1020 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
1022 (yy
< 0) ? - yy
: yy
);
1023 scm_remember_upto_here_1 (x
);
1024 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1025 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
1026 SCM_I_BIG_MPZ (result
),
1028 return scm_i_normbig (result
);
1031 else if (SCM_BIGP (y
))
1033 SCM result
= scm_i_mkbig ();
1034 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
1035 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1036 mpz_mod (SCM_I_BIG_MPZ (result
),
1038 SCM_I_BIG_MPZ (pos_y
));
1040 scm_remember_upto_here_1 (x
);
1041 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1042 mpz_add (SCM_I_BIG_MPZ (result
),
1044 SCM_I_BIG_MPZ (result
));
1045 scm_remember_upto_here_2 (y
, pos_y
);
1046 return scm_i_normbig (result
);
1049 SCM_WTA_DISPATCH_2 (g_scm_modulo
, x
, y
, SCM_ARG2
, s_scm_modulo
);
1052 SCM_WTA_DISPATCH_2 (g_scm_modulo
, x
, y
, SCM_ARG1
, s_scm_modulo
);
1056 static SCM
scm_i_inexact_euclidean_quotient (double x
, double y
);
1057 static SCM
scm_i_slow_exact_euclidean_quotient (SCM x
, SCM y
);
1059 SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient
, "euclidean-quotient", 2, 0, 0,
1061 "Return the integer @var{q} such that\n"
1062 "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1063 "where @math{0 <= @var{r} < abs(@var{y})}.\n"
1065 "(euclidean-quotient 123 10) @result{} 12\n"
1066 "(euclidean-quotient 123 -10) @result{} -12\n"
1067 "(euclidean-quotient -123 10) @result{} -13\n"
1068 "(euclidean-quotient -123 -10) @result{} 13\n"
1069 "(euclidean-quotient -123.2 -63.5) @result{} 2.0\n"
1070 "(euclidean-quotient 16/3 -10/7) @result{} -3\n"
1072 #define FUNC_NAME s_scm_euclidean_quotient
1074 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1076 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1078 scm_t_inum yy
= SCM_I_INUM (y
);
1079 if (SCM_UNLIKELY (yy
== 0))
1080 scm_num_overflow (s_scm_euclidean_quotient
);
1083 scm_t_inum xx
= SCM_I_INUM (x
);
1084 scm_t_inum qq
= xx
/ yy
;
1092 return SCM_I_MAKINUM (qq
);
1095 else if (SCM_BIGP (y
))
1097 if (SCM_I_INUM (x
) >= 0)
1100 return SCM_I_MAKINUM (- mpz_sgn (SCM_I_BIG_MPZ (y
)));
1102 else if (SCM_REALP (y
))
1103 return scm_i_inexact_euclidean_quotient
1104 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1105 else if (SCM_FRACTIONP (y
))
1106 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1108 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1109 s_scm_euclidean_quotient
);
1111 else if (SCM_BIGP (x
))
1113 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1115 scm_t_inum yy
= SCM_I_INUM (y
);
1116 if (SCM_UNLIKELY (yy
== 0))
1117 scm_num_overflow (s_scm_euclidean_quotient
);
1120 SCM q
= scm_i_mkbig ();
1122 mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (x
), yy
);
1125 mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (x
), -yy
);
1126 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1128 scm_remember_upto_here_1 (x
);
1129 return scm_i_normbig (q
);
1132 else if (SCM_BIGP (y
))
1134 SCM q
= scm_i_mkbig ();
1135 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1136 mpz_fdiv_q (SCM_I_BIG_MPZ (q
),
1140 mpz_cdiv_q (SCM_I_BIG_MPZ (q
),
1143 scm_remember_upto_here_2 (x
, y
);
1144 return scm_i_normbig (q
);
1146 else if (SCM_REALP (y
))
1147 return scm_i_inexact_euclidean_quotient
1148 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1149 else if (SCM_FRACTIONP (y
))
1150 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1152 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1153 s_scm_euclidean_quotient
);
1155 else if (SCM_REALP (x
))
1157 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1158 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1159 return scm_i_inexact_euclidean_quotient
1160 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1162 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1163 s_scm_euclidean_quotient
);
1165 else if (SCM_FRACTIONP (x
))
1168 return scm_i_inexact_euclidean_quotient
1169 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1171 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1174 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG1
,
1175 s_scm_euclidean_quotient
);
1180 scm_i_inexact_euclidean_quotient (double x
, double y
)
1182 if (SCM_LIKELY (y
> 0))
1183 return scm_from_double (floor (x
/ y
));
1184 else if (SCM_LIKELY (y
< 0))
1185 return scm_from_double (ceil (x
/ y
));
1187 scm_num_overflow (s_scm_euclidean_quotient
); /* or return a NaN? */
1192 /* Compute exact euclidean_quotient the slow way.
1193 We use this only if both arguments are exact,
1194 and at least one of them is a fraction */
1196 scm_i_slow_exact_euclidean_quotient (SCM x
, SCM y
)
1198 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1199 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG1
,
1200 s_scm_euclidean_quotient
);
1201 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1202 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1203 s_scm_euclidean_quotient
);
1204 else if (scm_is_true (scm_positive_p (y
)))
1205 return scm_floor (scm_divide (x
, y
));
1206 else if (scm_is_true (scm_negative_p (y
)))
1207 return scm_ceiling (scm_divide (x
, y
));
1209 scm_num_overflow (s_scm_euclidean_quotient
);
1212 static SCM
scm_i_inexact_euclidean_remainder (double x
, double y
);
1213 static SCM
scm_i_slow_exact_euclidean_remainder (SCM x
, SCM y
);
1215 SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder
, "euclidean-remainder", 2, 0, 0,
1217 "Return the real number @var{r} such that\n"
1218 "@math{0 <= @var{r} < abs(@var{y})} and\n"
1219 "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1220 "for some integer @var{q}.\n"
1222 "(euclidean-remainder 123 10) @result{} 3\n"
1223 "(euclidean-remainder 123 -10) @result{} 3\n"
1224 "(euclidean-remainder -123 10) @result{} 7\n"
1225 "(euclidean-remainder -123 -10) @result{} 7\n"
1226 "(euclidean-remainder -123.2 -63.5) @result{} 3.8\n"
1227 "(euclidean-remainder 16/3 -10/7) @result{} 22/21\n"
1229 #define FUNC_NAME s_scm_euclidean_remainder
1231 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1233 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1235 scm_t_inum yy
= SCM_I_INUM (y
);
1236 if (SCM_UNLIKELY (yy
== 0))
1237 scm_num_overflow (s_scm_euclidean_remainder
);
1240 scm_t_inum rr
= SCM_I_INUM (x
) % yy
;
1242 return SCM_I_MAKINUM (rr
);
1244 return SCM_I_MAKINUM (rr
+ yy
);
1246 return SCM_I_MAKINUM (rr
- yy
);
1249 else if (SCM_BIGP (y
))
1251 scm_t_inum xx
= SCM_I_INUM (x
);
1254 else if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1256 SCM r
= scm_i_mkbig ();
1257 mpz_sub_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1258 scm_remember_upto_here_1 (y
);
1259 return scm_i_normbig (r
);
1263 SCM r
= scm_i_mkbig ();
1264 mpz_add_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1265 scm_remember_upto_here_1 (y
);
1266 mpz_neg (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (r
));
1267 return scm_i_normbig (r
);
1270 else if (SCM_REALP (y
))
1271 return scm_i_inexact_euclidean_remainder
1272 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1273 else if (SCM_FRACTIONP (y
))
1274 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1276 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1277 s_scm_euclidean_remainder
);
1279 else if (SCM_BIGP (x
))
1281 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1283 scm_t_inum yy
= SCM_I_INUM (y
);
1284 if (SCM_UNLIKELY (yy
== 0))
1285 scm_num_overflow (s_scm_euclidean_remainder
);
1291 rr
= mpz_fdiv_ui (SCM_I_BIG_MPZ (x
), yy
);
1292 scm_remember_upto_here_1 (x
);
1293 return SCM_I_MAKINUM (rr
);
1296 else if (SCM_BIGP (y
))
1298 SCM r
= scm_i_mkbig ();
1299 mpz_mod (SCM_I_BIG_MPZ (r
),
1302 scm_remember_upto_here_2 (x
, y
);
1303 return scm_i_normbig (r
);
1305 else if (SCM_REALP (y
))
1306 return scm_i_inexact_euclidean_remainder
1307 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1308 else if (SCM_FRACTIONP (y
))
1309 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1311 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1312 s_scm_euclidean_remainder
);
1314 else if (SCM_REALP (x
))
1316 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1317 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1318 return scm_i_inexact_euclidean_remainder
1319 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1321 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1322 s_scm_euclidean_remainder
);
1324 else if (SCM_FRACTIONP (x
))
1327 return scm_i_inexact_euclidean_remainder
1328 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1330 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1333 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG1
,
1334 s_scm_euclidean_remainder
);
1339 scm_i_inexact_euclidean_remainder (double x
, double y
)
1343 /* Although it would be more efficient to use fmod here, we can't
1344 because it would in some cases produce results inconsistent with
1345 scm_i_inexact_euclidean_quotient, such that x != q * y + r (not
1346 even close). In particular, when x is very close to a multiple of
1347 y, then r might be either 0.0 or abs(y)-epsilon, but those two
1348 cases must correspond to different choices of q. If r = 0.0 then q
1349 must be x/y, and if r = abs(y) then q must be (x-r)/y. If quotient
1350 chooses one and remainder chooses the other, it would be bad. This
1351 problem was observed with x = 130.0 and y = 10/7. */
1352 if (SCM_LIKELY (y
> 0))
1354 else if (SCM_LIKELY (y
< 0))
1357 scm_num_overflow (s_scm_euclidean_remainder
); /* or return a NaN? */
1360 return scm_from_double (x
- q
* y
);
1363 /* Compute exact euclidean_remainder the slow way.
1364 We use this only if both arguments are exact,
1365 and at least one of them is a fraction */
1367 scm_i_slow_exact_euclidean_remainder (SCM x
, SCM y
)
1371 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1372 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG1
,
1373 s_scm_euclidean_remainder
);
1374 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1375 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1376 s_scm_euclidean_remainder
);
1377 else if (scm_is_true (scm_positive_p (y
)))
1378 q
= scm_floor (scm_divide (x
, y
));
1379 else if (scm_is_true (scm_negative_p (y
)))
1380 q
= scm_ceiling (scm_divide (x
, y
));
1382 scm_num_overflow (s_scm_euclidean_remainder
);
1383 return scm_difference (x
, scm_product (y
, q
));
1387 static SCM
scm_i_inexact_euclidean_divide (double x
, double y
);
1388 static SCM
scm_i_slow_exact_euclidean_divide (SCM x
, SCM y
);
1390 SCM_PRIMITIVE_GENERIC (scm_euclidean_divide
, "euclidean/", 2, 0, 0,
1392 "Return the integer @var{q} and the real number @var{r}\n"
1393 "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1394 "and @math{0 <= @var{r} < abs(@var{y})}.\n"
1396 "(euclidean/ 123 10) @result{} 12 and 3\n"
1397 "(euclidean/ 123 -10) @result{} -12 and 3\n"
1398 "(euclidean/ -123 10) @result{} -13 and 7\n"
1399 "(euclidean/ -123 -10) @result{} 13 and 7\n"
1400 "(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
1401 "(euclidean/ 16/3 -10/7) @result{} -3 and 22/21\n"
1403 #define FUNC_NAME s_scm_euclidean_divide
1405 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1407 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1409 scm_t_inum yy
= SCM_I_INUM (y
);
1410 if (SCM_UNLIKELY (yy
== 0))
1411 scm_num_overflow (s_scm_euclidean_divide
);
1414 scm_t_inum xx
= SCM_I_INUM (x
);
1415 scm_t_inum qq
= xx
/ yy
;
1416 scm_t_inum rr
= xx
- qq
* yy
;
1424 return scm_values (scm_list_2 (SCM_I_MAKINUM (qq
),
1425 SCM_I_MAKINUM (rr
)));
1428 else if (SCM_BIGP (y
))
1430 scm_t_inum xx
= SCM_I_INUM (x
);
1432 return scm_values (scm_list_2 (SCM_INUM0
, x
));
1433 else if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1435 SCM r
= scm_i_mkbig ();
1436 mpz_sub_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1437 scm_remember_upto_here_1 (y
);
1439 (scm_list_2 (SCM_I_MAKINUM (-1), scm_i_normbig (r
)));
1443 SCM r
= scm_i_mkbig ();
1444 mpz_add_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1445 scm_remember_upto_here_1 (y
);
1446 mpz_neg (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (r
));
1447 return scm_values (scm_list_2 (SCM_INUM1
, scm_i_normbig (r
)));
1450 else if (SCM_REALP (y
))
1451 return scm_i_inexact_euclidean_divide
1452 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1453 else if (SCM_FRACTIONP (y
))
1454 return scm_i_slow_exact_euclidean_divide (x
, y
);
1456 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1457 s_scm_euclidean_divide
);
1459 else if (SCM_BIGP (x
))
1461 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1463 scm_t_inum yy
= SCM_I_INUM (y
);
1464 if (SCM_UNLIKELY (yy
== 0))
1465 scm_num_overflow (s_scm_euclidean_divide
);
1468 SCM q
= scm_i_mkbig ();
1469 SCM r
= scm_i_mkbig ();
1471 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1472 SCM_I_BIG_MPZ (x
), yy
);
1475 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1476 SCM_I_BIG_MPZ (x
), -yy
);
1477 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1479 scm_remember_upto_here_1 (x
);
1480 return scm_values (scm_list_2 (scm_i_normbig (q
),
1481 scm_i_normbig (r
)));
1484 else if (SCM_BIGP (y
))
1486 SCM q
= scm_i_mkbig ();
1487 SCM r
= scm_i_mkbig ();
1488 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1489 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1490 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1492 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1493 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1494 scm_remember_upto_here_2 (x
, y
);
1495 return scm_values (scm_list_2 (scm_i_normbig (q
),
1496 scm_i_normbig (r
)));
1498 else if (SCM_REALP (y
))
1499 return scm_i_inexact_euclidean_divide
1500 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1501 else if (SCM_FRACTIONP (y
))
1502 return scm_i_slow_exact_euclidean_divide (x
, y
);
1504 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1505 s_scm_euclidean_divide
);
1507 else if (SCM_REALP (x
))
1509 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1510 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1511 return scm_i_inexact_euclidean_divide
1512 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1514 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1515 s_scm_euclidean_divide
);
1517 else if (SCM_FRACTIONP (x
))
1520 return scm_i_inexact_euclidean_divide
1521 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1523 return scm_i_slow_exact_euclidean_divide (x
, y
);
1526 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG1
,
1527 s_scm_euclidean_divide
);
1532 scm_i_inexact_euclidean_divide (double x
, double y
)
1536 if (SCM_LIKELY (y
> 0))
1538 else if (SCM_LIKELY (y
< 0))
1541 scm_num_overflow (s_scm_euclidean_divide
); /* or return a NaN? */
1545 return scm_values (scm_list_2 (scm_from_double (q
),
1546 scm_from_double (r
)));
1549 /* Compute exact euclidean quotient and remainder the slow way.
1550 We use this only if both arguments are exact,
1551 and at least one of them is a fraction */
1553 scm_i_slow_exact_euclidean_divide (SCM x
, SCM y
)
1557 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1558 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG1
,
1559 s_scm_euclidean_divide
);
1560 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1561 SCM_WTA_DISPATCH_2 (g_scm_euclidean_divide
, x
, y
, SCM_ARG2
,
1562 s_scm_euclidean_divide
);
1563 else if (scm_is_true (scm_positive_p (y
)))
1564 q
= scm_floor (scm_divide (x
, y
));
1565 else if (scm_is_true (scm_negative_p (y
)))
1566 q
= scm_ceiling (scm_divide (x
, y
));
1568 scm_num_overflow (s_scm_euclidean_divide
);
1569 r
= scm_difference (x
, scm_product (q
, y
));
1570 return scm_values (scm_list_2 (q
, r
));
1573 static SCM
scm_i_inexact_centered_quotient (double x
, double y
);
1574 static SCM
scm_i_bigint_centered_quotient (SCM x
, SCM y
);
1575 static SCM
scm_i_slow_exact_centered_quotient (SCM x
, SCM y
);
1577 SCM_PRIMITIVE_GENERIC (scm_centered_quotient
, "centered-quotient", 2, 0, 0,
1579 "Return the integer @var{q} such that\n"
1580 "@math{@var{x} = @var{q}*@var{y} + @var{r}} where\n"
1581 "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
1583 "(centered-quotient 123 10) @result{} 12\n"
1584 "(centered-quotient 123 -10) @result{} -12\n"
1585 "(centered-quotient -123 10) @result{} -12\n"
1586 "(centered-quotient -123 -10) @result{} 12\n"
1587 "(centered-quotient -123.2 -63.5) @result{} 2.0\n"
1588 "(centered-quotient 16/3 -10/7) @result{} -4\n"
1590 #define FUNC_NAME s_scm_centered_quotient
1592 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1594 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1596 scm_t_inum yy
= SCM_I_INUM (y
);
1597 if (SCM_UNLIKELY (yy
== 0))
1598 scm_num_overflow (s_scm_centered_quotient
);
1601 scm_t_inum xx
= SCM_I_INUM (x
);
1602 scm_t_inum qq
= xx
/ yy
;
1603 scm_t_inum rr
= xx
- qq
* yy
;
1604 if (SCM_LIKELY (xx
> 0))
1606 if (SCM_LIKELY (yy
> 0))
1608 if (rr
>= (yy
+ 1) / 2)
1613 if (rr
>= (1 - yy
) / 2)
1619 if (SCM_LIKELY (yy
> 0))
1630 return SCM_I_MAKINUM (qq
);
1633 else if (SCM_BIGP (y
))
1635 /* Pass a denormalized bignum version of x (even though it
1636 can fit in a fixnum) to scm_i_bigint_centered_quotient */
1637 return scm_i_bigint_centered_quotient
1638 (scm_i_long2big (SCM_I_INUM (x
)), y
);
1640 else if (SCM_REALP (y
))
1641 return scm_i_inexact_centered_quotient
1642 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1643 else if (SCM_FRACTIONP (y
))
1644 return scm_i_slow_exact_centered_quotient (x
, y
);
1646 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1647 s_scm_centered_quotient
);
1649 else if (SCM_BIGP (x
))
1651 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1653 scm_t_inum yy
= SCM_I_INUM (y
);
1654 if (SCM_UNLIKELY (yy
== 0))
1655 scm_num_overflow (s_scm_centered_quotient
);
1658 SCM q
= scm_i_mkbig ();
1660 /* Arrange for rr to initially be non-positive,
1661 because that simplifies the test to see
1662 if it is within the needed bounds. */
1665 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
1666 SCM_I_BIG_MPZ (x
), yy
);
1667 scm_remember_upto_here_1 (x
);
1669 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
1670 SCM_I_BIG_MPZ (q
), 1);
1674 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
1675 SCM_I_BIG_MPZ (x
), -yy
);
1676 scm_remember_upto_here_1 (x
);
1677 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1679 mpz_add_ui (SCM_I_BIG_MPZ (q
),
1680 SCM_I_BIG_MPZ (q
), 1);
1682 return scm_i_normbig (q
);
1685 else if (SCM_BIGP (y
))
1686 return scm_i_bigint_centered_quotient (x
, y
);
1687 else if (SCM_REALP (y
))
1688 return scm_i_inexact_centered_quotient
1689 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1690 else if (SCM_FRACTIONP (y
))
1691 return scm_i_slow_exact_centered_quotient (x
, y
);
1693 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1694 s_scm_centered_quotient
);
1696 else if (SCM_REALP (x
))
1698 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1699 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1700 return scm_i_inexact_centered_quotient
1701 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1703 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1704 s_scm_centered_quotient
);
1706 else if (SCM_FRACTIONP (x
))
1709 return scm_i_inexact_centered_quotient
1710 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1712 return scm_i_slow_exact_centered_quotient (x
, y
);
1715 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG1
,
1716 s_scm_centered_quotient
);
1721 scm_i_inexact_centered_quotient (double x
, double y
)
1723 if (SCM_LIKELY (y
> 0))
1724 return scm_from_double (floor (x
/y
+ 0.5));
1725 else if (SCM_LIKELY (y
< 0))
1726 return scm_from_double (ceil (x
/y
- 0.5));
1728 scm_num_overflow (s_scm_centered_quotient
); /* or return a NaN? */
1733 /* Assumes that both x and y are bigints, though
1734 x might be able to fit into a fixnum. */
1736 scm_i_bigint_centered_quotient (SCM x
, SCM y
)
1740 /* Note that x might be small enough to fit into a
1741 fixnum, so we must not let it escape into the wild */
1745 /* min_r will eventually become -abs(y)/2 */
1746 min_r
= scm_i_mkbig ();
1747 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
1748 SCM_I_BIG_MPZ (y
), 1);
1750 /* Arrange for rr to initially be non-positive,
1751 because that simplifies the test to see
1752 if it is within the needed bounds. */
1753 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1755 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1756 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1757 scm_remember_upto_here_2 (x
, y
);
1758 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
1759 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1760 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
1761 SCM_I_BIG_MPZ (q
), 1);
1765 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1766 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1767 scm_remember_upto_here_2 (x
, y
);
1768 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1769 mpz_add_ui (SCM_I_BIG_MPZ (q
),
1770 SCM_I_BIG_MPZ (q
), 1);
1772 scm_remember_upto_here_2 (r
, min_r
);
1773 return scm_i_normbig (q
);
1776 /* Compute exact centered quotient the slow way.
1777 We use this only if both arguments are exact,
1778 and at least one of them is a fraction */
1780 scm_i_slow_exact_centered_quotient (SCM x
, SCM y
)
1782 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1783 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG1
,
1784 s_scm_centered_quotient
);
1785 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1786 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1787 s_scm_centered_quotient
);
1788 else if (scm_is_true (scm_positive_p (y
)))
1789 return scm_floor (scm_sum (scm_divide (x
, y
),
1791 else if (scm_is_true (scm_negative_p (y
)))
1792 return scm_ceiling (scm_difference (scm_divide (x
, y
),
1795 scm_num_overflow (s_scm_centered_quotient
);
1798 static SCM
scm_i_inexact_centered_remainder (double x
, double y
);
1799 static SCM
scm_i_bigint_centered_remainder (SCM x
, SCM y
);
1800 static SCM
scm_i_slow_exact_centered_remainder (SCM x
, SCM y
);
1802 SCM_PRIMITIVE_GENERIC (scm_centered_remainder
, "centered-remainder", 2, 0, 0,
1804 "Return the real number @var{r} such that\n"
1805 "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}\n"
1806 "and @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1807 "for some integer @var{q}.\n"
1809 "(centered-remainder 123 10) @result{} 3\n"
1810 "(centered-remainder 123 -10) @result{} 3\n"
1811 "(centered-remainder -123 10) @result{} -3\n"
1812 "(centered-remainder -123 -10) @result{} -3\n"
1813 "(centered-remainder -123.2 -63.5) @result{} 3.8\n"
1814 "(centered-remainder 16/3 -10/7) @result{} -8/21\n"
1816 #define FUNC_NAME s_scm_centered_remainder
1818 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1820 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1822 scm_t_inum yy
= SCM_I_INUM (y
);
1823 if (SCM_UNLIKELY (yy
== 0))
1824 scm_num_overflow (s_scm_centered_remainder
);
1827 scm_t_inum xx
= SCM_I_INUM (x
);
1828 scm_t_inum rr
= xx
% yy
;
1829 if (SCM_LIKELY (xx
> 0))
1831 if (SCM_LIKELY (yy
> 0))
1833 if (rr
>= (yy
+ 1) / 2)
1838 if (rr
>= (1 - yy
) / 2)
1844 if (SCM_LIKELY (yy
> 0))
1855 return SCM_I_MAKINUM (rr
);
1858 else if (SCM_BIGP (y
))
1860 /* Pass a denormalized bignum version of x (even though it
1861 can fit in a fixnum) to scm_i_bigint_centered_remainder */
1862 return scm_i_bigint_centered_remainder
1863 (scm_i_long2big (SCM_I_INUM (x
)), y
);
1865 else if (SCM_REALP (y
))
1866 return scm_i_inexact_centered_remainder
1867 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1868 else if (SCM_FRACTIONP (y
))
1869 return scm_i_slow_exact_centered_remainder (x
, y
);
1871 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1872 s_scm_centered_remainder
);
1874 else if (SCM_BIGP (x
))
1876 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1878 scm_t_inum yy
= SCM_I_INUM (y
);
1879 if (SCM_UNLIKELY (yy
== 0))
1880 scm_num_overflow (s_scm_centered_remainder
);
1884 /* Arrange for rr to initially be non-positive,
1885 because that simplifies the test to see
1886 if it is within the needed bounds. */
1889 rr
= - mpz_cdiv_ui (SCM_I_BIG_MPZ (x
), yy
);
1890 scm_remember_upto_here_1 (x
);
1896 rr
= - mpz_cdiv_ui (SCM_I_BIG_MPZ (x
), -yy
);
1897 scm_remember_upto_here_1 (x
);
1901 return SCM_I_MAKINUM (rr
);
1904 else if (SCM_BIGP (y
))
1905 return scm_i_bigint_centered_remainder (x
, y
);
1906 else if (SCM_REALP (y
))
1907 return scm_i_inexact_centered_remainder
1908 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1909 else if (SCM_FRACTIONP (y
))
1910 return scm_i_slow_exact_centered_remainder (x
, y
);
1912 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1913 s_scm_centered_remainder
);
1915 else if (SCM_REALP (x
))
1917 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1918 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1919 return scm_i_inexact_centered_remainder
1920 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1922 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1923 s_scm_centered_remainder
);
1925 else if (SCM_FRACTIONP (x
))
1928 return scm_i_inexact_centered_remainder
1929 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1931 return scm_i_slow_exact_centered_remainder (x
, y
);
1934 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG1
,
1935 s_scm_centered_remainder
);
1940 scm_i_inexact_centered_remainder (double x
, double y
)
1944 /* Although it would be more efficient to use fmod here, we can't
1945 because it would in some cases produce results inconsistent with
1946 scm_i_inexact_centered_quotient, such that x != r + q * y (not even
1947 close). In particular, when x-y/2 is very close to a multiple of
1948 y, then r might be either -abs(y/2) or abs(y/2)-epsilon, but those
1949 two cases must correspond to different choices of q. If quotient
1950 chooses one and remainder chooses the other, it would be bad. */
1951 if (SCM_LIKELY (y
> 0))
1952 q
= floor (x
/y
+ 0.5);
1953 else if (SCM_LIKELY (y
< 0))
1954 q
= ceil (x
/y
- 0.5);
1956 scm_num_overflow (s_scm_centered_remainder
); /* or return a NaN? */
1959 return scm_from_double (x
- q
* y
);
1962 /* Assumes that both x and y are bigints, though
1963 x might be able to fit into a fixnum. */
1965 scm_i_bigint_centered_remainder (SCM x
, SCM y
)
1969 /* Note that x might be small enough to fit into a
1970 fixnum, so we must not let it escape into the wild */
1973 /* min_r will eventually become -abs(y)/2 */
1974 min_r
= scm_i_mkbig ();
1975 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
1976 SCM_I_BIG_MPZ (y
), 1);
1978 /* Arrange for rr to initially be non-positive,
1979 because that simplifies the test to see
1980 if it is within the needed bounds. */
1981 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1983 mpz_cdiv_r (SCM_I_BIG_MPZ (r
),
1984 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1985 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
1986 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1987 mpz_add (SCM_I_BIG_MPZ (r
),
1993 mpz_fdiv_r (SCM_I_BIG_MPZ (r
),
1994 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1995 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1996 mpz_sub (SCM_I_BIG_MPZ (r
),
2000 scm_remember_upto_here_2 (x
, y
);
2001 return scm_i_normbig (r
);
2004 /* Compute exact centered_remainder the slow way.
2005 We use this only if both arguments are exact,
2006 and at least one of them is a fraction */
2008 scm_i_slow_exact_centered_remainder (SCM x
, SCM y
)
2012 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
2013 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG1
,
2014 s_scm_centered_remainder
);
2015 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
2016 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
2017 s_scm_centered_remainder
);
2018 else if (scm_is_true (scm_positive_p (y
)))
2019 q
= scm_floor (scm_sum (scm_divide (x
, y
), exactly_one_half
));
2020 else if (scm_is_true (scm_negative_p (y
)))
2021 q
= scm_ceiling (scm_difference (scm_divide (x
, y
), exactly_one_half
));
2023 scm_num_overflow (s_scm_centered_remainder
);
2024 return scm_difference (x
, scm_product (y
, q
));
2028 static SCM
scm_i_inexact_centered_divide (double x
, double y
);
2029 static SCM
scm_i_bigint_centered_divide (SCM x
, SCM y
);
2030 static SCM
scm_i_slow_exact_centered_divide (SCM x
, SCM y
);
2032 SCM_PRIMITIVE_GENERIC (scm_centered_divide
, "centered/", 2, 0, 0,
2034 "Return the integer @var{q} and the real number @var{r}\n"
2035 "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
2036 "and @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
2038 "(centered/ 123 10) @result{} 12 and 3\n"
2039 "(centered/ 123 -10) @result{} -12 and 3\n"
2040 "(centered/ -123 10) @result{} -12 and -3\n"
2041 "(centered/ -123 -10) @result{} 12 and -3\n"
2042 "(centered/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
2043 "(centered/ 16/3 -10/7) @result{} -4 and -8/21\n"
2045 #define FUNC_NAME s_scm_centered_divide
2047 if (SCM_LIKELY (SCM_I_INUMP (x
)))
2049 if (SCM_LIKELY (SCM_I_INUMP (y
)))
2051 scm_t_inum yy
= SCM_I_INUM (y
);
2052 if (SCM_UNLIKELY (yy
== 0))
2053 scm_num_overflow (s_scm_centered_divide
);
2056 scm_t_inum xx
= SCM_I_INUM (x
);
2057 scm_t_inum qq
= xx
/ yy
;
2058 scm_t_inum rr
= xx
- qq
* yy
;
2059 if (SCM_LIKELY (xx
> 0))
2061 if (SCM_LIKELY (yy
> 0))
2063 if (rr
>= (yy
+ 1) / 2)
2068 if (rr
>= (1 - yy
) / 2)
2074 if (SCM_LIKELY (yy
> 0))
2085 return scm_values (scm_list_2 (SCM_I_MAKINUM (qq
),
2086 SCM_I_MAKINUM (rr
)));
2089 else if (SCM_BIGP (y
))
2091 /* Pass a denormalized bignum version of x (even though it
2092 can fit in a fixnum) to scm_i_bigint_centered_divide */
2093 return scm_i_bigint_centered_divide
2094 (scm_i_long2big (SCM_I_INUM (x
)), y
);
2096 else if (SCM_REALP (y
))
2097 return scm_i_inexact_centered_divide
2098 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
2099 else if (SCM_FRACTIONP (y
))
2100 return scm_i_slow_exact_centered_divide (x
, y
);
2102 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2103 s_scm_centered_divide
);
2105 else if (SCM_BIGP (x
))
2107 if (SCM_LIKELY (SCM_I_INUMP (y
)))
2109 scm_t_inum yy
= SCM_I_INUM (y
);
2110 if (SCM_UNLIKELY (yy
== 0))
2111 scm_num_overflow (s_scm_centered_divide
);
2114 SCM q
= scm_i_mkbig ();
2116 /* Arrange for rr to initially be non-positive,
2117 because that simplifies the test to see
2118 if it is within the needed bounds. */
2121 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
2122 SCM_I_BIG_MPZ (x
), yy
);
2123 scm_remember_upto_here_1 (x
);
2126 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
2127 SCM_I_BIG_MPZ (q
), 1);
2133 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
2134 SCM_I_BIG_MPZ (x
), -yy
);
2135 scm_remember_upto_here_1 (x
);
2136 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
2139 mpz_add_ui (SCM_I_BIG_MPZ (q
),
2140 SCM_I_BIG_MPZ (q
), 1);
2144 return scm_values (scm_list_2 (scm_i_normbig (q
),
2145 SCM_I_MAKINUM (rr
)));
2148 else if (SCM_BIGP (y
))
2149 return scm_i_bigint_centered_divide (x
, y
);
2150 else if (SCM_REALP (y
))
2151 return scm_i_inexact_centered_divide
2152 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
2153 else if (SCM_FRACTIONP (y
))
2154 return scm_i_slow_exact_centered_divide (x
, y
);
2156 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2157 s_scm_centered_divide
);
2159 else if (SCM_REALP (x
))
2161 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
2162 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
2163 return scm_i_inexact_centered_divide
2164 (SCM_REAL_VALUE (x
), scm_to_double (y
));
2166 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2167 s_scm_centered_divide
);
2169 else if (SCM_FRACTIONP (x
))
2172 return scm_i_inexact_centered_divide
2173 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
2175 return scm_i_slow_exact_centered_divide (x
, y
);
2178 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG1
,
2179 s_scm_centered_divide
);
2184 scm_i_inexact_centered_divide (double x
, double y
)
2188 if (SCM_LIKELY (y
> 0))
2189 q
= floor (x
/y
+ 0.5);
2190 else if (SCM_LIKELY (y
< 0))
2191 q
= ceil (x
/y
- 0.5);
2193 scm_num_overflow (s_scm_centered_divide
); /* or return a NaN? */
2197 return scm_values (scm_list_2 (scm_from_double (q
),
2198 scm_from_double (r
)));
2201 /* Assumes that both x and y are bigints, though
2202 x might be able to fit into a fixnum. */
2204 scm_i_bigint_centered_divide (SCM x
, SCM y
)
2208 /* Note that x might be small enough to fit into a
2209 fixnum, so we must not let it escape into the wild */
2213 /* min_r will eventually become -abs(y/2) */
2214 min_r
= scm_i_mkbig ();
2215 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
2216 SCM_I_BIG_MPZ (y
), 1);
2218 /* Arrange for rr to initially be non-positive,
2219 because that simplifies the test to see
2220 if it is within the needed bounds. */
2221 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
2223 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
2224 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2225 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
2226 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2228 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
2229 SCM_I_BIG_MPZ (q
), 1);
2230 mpz_add (SCM_I_BIG_MPZ (r
),
2237 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
2238 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2239 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2241 mpz_add_ui (SCM_I_BIG_MPZ (q
),
2242 SCM_I_BIG_MPZ (q
), 1);
2243 mpz_sub (SCM_I_BIG_MPZ (r
),
2248 scm_remember_upto_here_2 (x
, y
);
2249 return scm_values (scm_list_2 (scm_i_normbig (q
),
2250 scm_i_normbig (r
)));
2253 /* Compute exact centered quotient and remainder the slow way.
2254 We use this only if both arguments are exact,
2255 and at least one of them is a fraction */
2257 scm_i_slow_exact_centered_divide (SCM x
, SCM y
)
2261 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
2262 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG1
,
2263 s_scm_centered_divide
);
2264 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
2265 SCM_WTA_DISPATCH_2 (g_scm_centered_divide
, x
, y
, SCM_ARG2
,
2266 s_scm_centered_divide
);
2267 else if (scm_is_true (scm_positive_p (y
)))
2268 q
= scm_floor (scm_sum (scm_divide (x
, y
),
2270 else if (scm_is_true (scm_negative_p (y
)))
2271 q
= scm_ceiling (scm_difference (scm_divide (x
, y
),
2274 scm_num_overflow (s_scm_centered_divide
);
2275 r
= scm_difference (x
, scm_product (q
, y
));
2276 return scm_values (scm_list_2 (q
, r
));
2280 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
2281 (SCM x
, SCM y
, SCM rest
),
2282 "Return the greatest common divisor of all parameter values.\n"
2283 "If called without arguments, 0 is returned.")
2284 #define FUNC_NAME s_scm_i_gcd
2286 while (!scm_is_null (rest
))
2287 { x
= scm_gcd (x
, y
);
2289 rest
= scm_cdr (rest
);
2291 return scm_gcd (x
, y
);
2295 #define s_gcd s_scm_i_gcd
2296 #define g_gcd g_scm_i_gcd
2299 scm_gcd (SCM x
, SCM y
)
2302 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
2304 if (SCM_I_INUMP (x
))
2306 if (SCM_I_INUMP (y
))
2308 scm_t_inum xx
= SCM_I_INUM (x
);
2309 scm_t_inum yy
= SCM_I_INUM (y
);
2310 scm_t_inum u
= xx
< 0 ? -xx
: xx
;
2311 scm_t_inum v
= yy
< 0 ? -yy
: yy
;
2321 /* Determine a common factor 2^k */
2322 while (!(1 & (u
| v
)))
2328 /* Now, any factor 2^n can be eliminated */
2348 return (SCM_POSFIXABLE (result
)
2349 ? SCM_I_MAKINUM (result
)
2350 : scm_i_inum2big (result
));
2352 else if (SCM_BIGP (y
))
2358 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
2360 else if (SCM_BIGP (x
))
2362 if (SCM_I_INUMP (y
))
2367 yy
= SCM_I_INUM (y
);
2372 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
2373 scm_remember_upto_here_1 (x
);
2374 return (SCM_POSFIXABLE (result
)
2375 ? SCM_I_MAKINUM (result
)
2376 : scm_from_unsigned_integer (result
));
2378 else if (SCM_BIGP (y
))
2380 SCM result
= scm_i_mkbig ();
2381 mpz_gcd (SCM_I_BIG_MPZ (result
),
2384 scm_remember_upto_here_2 (x
, y
);
2385 return scm_i_normbig (result
);
2388 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
2391 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
2394 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
2395 (SCM x
, SCM y
, SCM rest
),
2396 "Return the least common multiple of the arguments.\n"
2397 "If called without arguments, 1 is returned.")
2398 #define FUNC_NAME s_scm_i_lcm
2400 while (!scm_is_null (rest
))
2401 { x
= scm_lcm (x
, y
);
2403 rest
= scm_cdr (rest
);
2405 return scm_lcm (x
, y
);
2409 #define s_lcm s_scm_i_lcm
2410 #define g_lcm g_scm_i_lcm
2413 scm_lcm (SCM n1
, SCM n2
)
2415 if (SCM_UNBNDP (n2
))
2417 if (SCM_UNBNDP (n1
))
2418 return SCM_I_MAKINUM (1L);
2419 n2
= SCM_I_MAKINUM (1L);
2422 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
2423 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
2424 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
2425 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
2427 if (SCM_I_INUMP (n1
))
2429 if (SCM_I_INUMP (n2
))
2431 SCM d
= scm_gcd (n1
, n2
);
2432 if (scm_is_eq (d
, SCM_INUM0
))
2435 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
2439 /* inum n1, big n2 */
2442 SCM result
= scm_i_mkbig ();
2443 scm_t_inum nn1
= SCM_I_INUM (n1
);
2444 if (nn1
== 0) return SCM_INUM0
;
2445 if (nn1
< 0) nn1
= - nn1
;
2446 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
2447 scm_remember_upto_here_1 (n2
);
2455 if (SCM_I_INUMP (n2
))
2462 SCM result
= scm_i_mkbig ();
2463 mpz_lcm(SCM_I_BIG_MPZ (result
),
2465 SCM_I_BIG_MPZ (n2
));
2466 scm_remember_upto_here_2(n1
, n2
);
2467 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
2473 /* Emulating 2's complement bignums with sign magnitude arithmetic:
2478 + + + x (map digit:logand X Y)
2479 + - + x (map digit:logand X (lognot (+ -1 Y)))
2480 - + + y (map digit:logand (lognot (+ -1 X)) Y)
2481 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
2486 + + + (map digit:logior X Y)
2487 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
2488 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
2489 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
2494 + + + (map digit:logxor X Y)
2495 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
2496 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
2497 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
2502 + + (any digit:logand X Y)
2503 + - (any digit:logand X (lognot (+ -1 Y)))
2504 - + (any digit:logand (lognot (+ -1 X)) Y)
2509 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
2510 (SCM x
, SCM y
, SCM rest
),
2511 "Return the bitwise AND of the integer arguments.\n\n"
2513 "(logand) @result{} -1\n"
2514 "(logand 7) @result{} 7\n"
2515 "(logand #b111 #b011 #b001) @result{} 1\n"
2517 #define FUNC_NAME s_scm_i_logand
2519 while (!scm_is_null (rest
))
2520 { x
= scm_logand (x
, y
);
2522 rest
= scm_cdr (rest
);
2524 return scm_logand (x
, y
);
2528 #define s_scm_logand s_scm_i_logand
2530 SCM
scm_logand (SCM n1
, SCM n2
)
2531 #define FUNC_NAME s_scm_logand
2535 if (SCM_UNBNDP (n2
))
2537 if (SCM_UNBNDP (n1
))
2538 return SCM_I_MAKINUM (-1);
2539 else if (!SCM_NUMBERP (n1
))
2540 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2541 else if (SCM_NUMBERP (n1
))
2544 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2547 if (SCM_I_INUMP (n1
))
2549 nn1
= SCM_I_INUM (n1
);
2550 if (SCM_I_INUMP (n2
))
2552 scm_t_inum nn2
= SCM_I_INUM (n2
);
2553 return SCM_I_MAKINUM (nn1
& nn2
);
2555 else if SCM_BIGP (n2
)
2561 SCM result_z
= scm_i_mkbig ();
2563 mpz_init_set_si (nn1_z
, nn1
);
2564 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2565 scm_remember_upto_here_1 (n2
);
2567 return scm_i_normbig (result_z
);
2571 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2573 else if (SCM_BIGP (n1
))
2575 if (SCM_I_INUMP (n2
))
2578 nn1
= SCM_I_INUM (n1
);
2581 else if (SCM_BIGP (n2
))
2583 SCM result_z
= scm_i_mkbig ();
2584 mpz_and (SCM_I_BIG_MPZ (result_z
),
2586 SCM_I_BIG_MPZ (n2
));
2587 scm_remember_upto_here_2 (n1
, n2
);
2588 return scm_i_normbig (result_z
);
2591 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2594 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2599 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
2600 (SCM x
, SCM y
, SCM rest
),
2601 "Return the bitwise OR of the integer arguments.\n\n"
2603 "(logior) @result{} 0\n"
2604 "(logior 7) @result{} 7\n"
2605 "(logior #b000 #b001 #b011) @result{} 3\n"
2607 #define FUNC_NAME s_scm_i_logior
2609 while (!scm_is_null (rest
))
2610 { x
= scm_logior (x
, y
);
2612 rest
= scm_cdr (rest
);
2614 return scm_logior (x
, y
);
2618 #define s_scm_logior s_scm_i_logior
2620 SCM
scm_logior (SCM n1
, SCM n2
)
2621 #define FUNC_NAME s_scm_logior
2625 if (SCM_UNBNDP (n2
))
2627 if (SCM_UNBNDP (n1
))
2629 else if (SCM_NUMBERP (n1
))
2632 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2635 if (SCM_I_INUMP (n1
))
2637 nn1
= SCM_I_INUM (n1
);
2638 if (SCM_I_INUMP (n2
))
2640 long nn2
= SCM_I_INUM (n2
);
2641 return SCM_I_MAKINUM (nn1
| nn2
);
2643 else if (SCM_BIGP (n2
))
2649 SCM result_z
= scm_i_mkbig ();
2651 mpz_init_set_si (nn1_z
, nn1
);
2652 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2653 scm_remember_upto_here_1 (n2
);
2655 return scm_i_normbig (result_z
);
2659 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2661 else if (SCM_BIGP (n1
))
2663 if (SCM_I_INUMP (n2
))
2666 nn1
= SCM_I_INUM (n1
);
2669 else if (SCM_BIGP (n2
))
2671 SCM result_z
= scm_i_mkbig ();
2672 mpz_ior (SCM_I_BIG_MPZ (result_z
),
2674 SCM_I_BIG_MPZ (n2
));
2675 scm_remember_upto_here_2 (n1
, n2
);
2676 return scm_i_normbig (result_z
);
2679 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2682 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2687 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
2688 (SCM x
, SCM y
, SCM rest
),
2689 "Return the bitwise XOR of the integer arguments. A bit is\n"
2690 "set in the result if it is set in an odd number of arguments.\n"
2692 "(logxor) @result{} 0\n"
2693 "(logxor 7) @result{} 7\n"
2694 "(logxor #b000 #b001 #b011) @result{} 2\n"
2695 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
2697 #define FUNC_NAME s_scm_i_logxor
2699 while (!scm_is_null (rest
))
2700 { x
= scm_logxor (x
, y
);
2702 rest
= scm_cdr (rest
);
2704 return scm_logxor (x
, y
);
2708 #define s_scm_logxor s_scm_i_logxor
2710 SCM
scm_logxor (SCM n1
, SCM n2
)
2711 #define FUNC_NAME s_scm_logxor
2715 if (SCM_UNBNDP (n2
))
2717 if (SCM_UNBNDP (n1
))
2719 else if (SCM_NUMBERP (n1
))
2722 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2725 if (SCM_I_INUMP (n1
))
2727 nn1
= SCM_I_INUM (n1
);
2728 if (SCM_I_INUMP (n2
))
2730 scm_t_inum nn2
= SCM_I_INUM (n2
);
2731 return SCM_I_MAKINUM (nn1
^ nn2
);
2733 else if (SCM_BIGP (n2
))
2737 SCM result_z
= scm_i_mkbig ();
2739 mpz_init_set_si (nn1_z
, nn1
);
2740 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2741 scm_remember_upto_here_1 (n2
);
2743 return scm_i_normbig (result_z
);
2747 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2749 else if (SCM_BIGP (n1
))
2751 if (SCM_I_INUMP (n2
))
2754 nn1
= SCM_I_INUM (n1
);
2757 else if (SCM_BIGP (n2
))
2759 SCM result_z
= scm_i_mkbig ();
2760 mpz_xor (SCM_I_BIG_MPZ (result_z
),
2762 SCM_I_BIG_MPZ (n2
));
2763 scm_remember_upto_here_2 (n1
, n2
);
2764 return scm_i_normbig (result_z
);
2767 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2770 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2775 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
2777 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
2778 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
2779 "without actually calculating the @code{logand}, just testing\n"
2783 "(logtest #b0100 #b1011) @result{} #f\n"
2784 "(logtest #b0100 #b0111) @result{} #t\n"
2786 #define FUNC_NAME s_scm_logtest
2790 if (SCM_I_INUMP (j
))
2792 nj
= SCM_I_INUM (j
);
2793 if (SCM_I_INUMP (k
))
2795 scm_t_inum nk
= SCM_I_INUM (k
);
2796 return scm_from_bool (nj
& nk
);
2798 else if (SCM_BIGP (k
))
2806 mpz_init_set_si (nj_z
, nj
);
2807 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
2808 scm_remember_upto_here_1 (k
);
2809 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
2815 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
2817 else if (SCM_BIGP (j
))
2819 if (SCM_I_INUMP (k
))
2822 nj
= SCM_I_INUM (j
);
2825 else if (SCM_BIGP (k
))
2829 mpz_init (result_z
);
2833 scm_remember_upto_here_2 (j
, k
);
2834 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
2835 mpz_clear (result_z
);
2839 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
2842 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
2847 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
2849 "Test whether bit number @var{index} in @var{j} is set.\n"
2850 "@var{index} starts from 0 for the least significant bit.\n"
2853 "(logbit? 0 #b1101) @result{} #t\n"
2854 "(logbit? 1 #b1101) @result{} #f\n"
2855 "(logbit? 2 #b1101) @result{} #t\n"
2856 "(logbit? 3 #b1101) @result{} #t\n"
2857 "(logbit? 4 #b1101) @result{} #f\n"
2859 #define FUNC_NAME s_scm_logbit_p
2861 unsigned long int iindex
;
2862 iindex
= scm_to_ulong (index
);
2864 if (SCM_I_INUMP (j
))
2866 /* bits above what's in an inum follow the sign bit */
2867 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
2868 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
2870 else if (SCM_BIGP (j
))
2872 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
2873 scm_remember_upto_here_1 (j
);
2874 return scm_from_bool (val
);
2877 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
2882 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
2884 "Return the integer which is the ones-complement of the integer\n"
2888 "(number->string (lognot #b10000000) 2)\n"
2889 " @result{} \"-10000001\"\n"
2890 "(number->string (lognot #b0) 2)\n"
2891 " @result{} \"-1\"\n"
2893 #define FUNC_NAME s_scm_lognot
2895 if (SCM_I_INUMP (n
)) {
2896 /* No overflow here, just need to toggle all the bits making up the inum.
2897 Enhancement: No need to strip the tag and add it back, could just xor
2898 a block of 1 bits, if that worked with the various debug versions of
2900 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
2902 } else if (SCM_BIGP (n
)) {
2903 SCM result
= scm_i_mkbig ();
2904 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
2905 scm_remember_upto_here_1 (n
);
2909 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2914 /* returns 0 if IN is not an integer. OUT must already be
2917 coerce_to_big (SCM in
, mpz_t out
)
2920 mpz_set (out
, SCM_I_BIG_MPZ (in
));
2921 else if (SCM_I_INUMP (in
))
2922 mpz_set_si (out
, SCM_I_INUM (in
));
2929 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
2930 (SCM n
, SCM k
, SCM m
),
2931 "Return @var{n} raised to the integer exponent\n"
2932 "@var{k}, modulo @var{m}.\n"
2935 "(modulo-expt 2 3 5)\n"
2938 #define FUNC_NAME s_scm_modulo_expt
2944 /* There are two classes of error we might encounter --
2945 1) Math errors, which we'll report by calling scm_num_overflow,
2947 2) wrong-type errors, which of course we'll report by calling
2949 We don't report those errors immediately, however; instead we do
2950 some cleanup first. These variables tell us which error (if
2951 any) we should report after cleaning up.
2953 int report_overflow
= 0;
2955 int position_of_wrong_type
= 0;
2956 SCM value_of_wrong_type
= SCM_INUM0
;
2958 SCM result
= SCM_UNDEFINED
;
2964 if (scm_is_eq (m
, SCM_INUM0
))
2966 report_overflow
= 1;
2970 if (!coerce_to_big (n
, n_tmp
))
2972 value_of_wrong_type
= n
;
2973 position_of_wrong_type
= 1;
2977 if (!coerce_to_big (k
, k_tmp
))
2979 value_of_wrong_type
= k
;
2980 position_of_wrong_type
= 2;
2984 if (!coerce_to_big (m
, m_tmp
))
2986 value_of_wrong_type
= m
;
2987 position_of_wrong_type
= 3;
2991 /* if the exponent K is negative, and we simply call mpz_powm, we
2992 will get a divide-by-zero exception when an inverse 1/n mod m
2993 doesn't exist (or is not unique). Since exceptions are hard to
2994 handle, we'll attempt the inversion "by hand" -- that way, we get
2995 a simple failure code, which is easy to handle. */
2997 if (-1 == mpz_sgn (k_tmp
))
2999 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
3001 report_overflow
= 1;
3004 mpz_neg (k_tmp
, k_tmp
);
3007 result
= scm_i_mkbig ();
3008 mpz_powm (SCM_I_BIG_MPZ (result
),
3013 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
3014 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
3021 if (report_overflow
)
3022 scm_num_overflow (FUNC_NAME
);
3024 if (position_of_wrong_type
)
3025 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
3026 value_of_wrong_type
);
3028 return scm_i_normbig (result
);
3032 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
3034 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
3035 "exact integer, @var{n} can be any number.\n"
3037 "Negative @var{k} is supported, and results in\n"
3038 "@math{1/@var{n}^abs(@var{k})} in the usual way.\n"
3039 "@math{@var{n}^0} is 1, as usual, and that\n"
3040 "includes @math{0^0} is 1.\n"
3043 "(integer-expt 2 5) @result{} 32\n"
3044 "(integer-expt -3 3) @result{} -27\n"
3045 "(integer-expt 5 -3) @result{} 1/125\n"
3046 "(integer-expt 0 0) @result{} 1\n"
3048 #define FUNC_NAME s_scm_integer_expt
3051 SCM z_i2
= SCM_BOOL_F
;
3053 SCM acc
= SCM_I_MAKINUM (1L);
3055 /* Specifically refrain from checking the type of the first argument.
3056 This allows us to exponentiate any object that can be multiplied.
3057 If we must raise to a negative power, we must also be able to
3058 take its reciprocal. */
3059 if (!SCM_LIKELY (SCM_I_INUMP (k
)) && !SCM_LIKELY (SCM_BIGP (k
)))
3060 SCM_WRONG_TYPE_ARG (2, k
);
3062 if (SCM_UNLIKELY (scm_is_eq (k
, SCM_INUM0
)))
3063 return SCM_INUM1
; /* n^(exact0) is exact 1, regardless of n */
3064 else if (SCM_UNLIKELY (scm_is_eq (n
, SCM_I_MAKINUM (-1L))))
3065 return scm_is_false (scm_even_p (k
)) ? n
: SCM_INUM1
;
3066 /* The next check is necessary only because R6RS specifies different
3067 behavior for 0^(-k) than for (/ 0). If n is not a scheme number,
3068 we simply skip this case and move on. */
3069 else if (SCM_NUMBERP (n
) && scm_is_true (scm_zero_p (n
)))
3071 /* k cannot be 0 at this point, because we
3072 have already checked for that case above */
3073 if (scm_is_true (scm_positive_p (k
)))
3075 else /* return NaN for (0 ^ k) for negative k per R6RS */
3079 if (SCM_I_INUMP (k
))
3080 i2
= SCM_I_INUM (k
);
3081 else if (SCM_BIGP (k
))
3083 z_i2
= scm_i_clonebig (k
, 1);
3084 scm_remember_upto_here_1 (k
);
3088 SCM_WRONG_TYPE_ARG (2, k
);
3092 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
3094 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
3095 n
= scm_divide (n
, SCM_UNDEFINED
);
3099 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
3103 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
3105 return scm_product (acc
, n
);
3107 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
3108 acc
= scm_product (acc
, n
);
3109 n
= scm_product (n
, n
);
3110 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
3118 n
= scm_divide (n
, SCM_UNDEFINED
);
3125 return scm_product (acc
, n
);
3127 acc
= scm_product (acc
, n
);
3128 n
= scm_product (n
, n
);
3135 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
3137 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
3138 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
3140 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
3141 "@var{cnt} is negative it's a division, rounded towards negative\n"
3142 "infinity. (Note that this is not the same rounding as\n"
3143 "@code{quotient} does.)\n"
3145 "With @var{n} viewed as an infinite precision twos complement,\n"
3146 "@code{ash} means a left shift introducing zero bits, or a right\n"
3147 "shift dropping bits.\n"
3150 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
3151 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
3153 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
3154 "(ash -23 -2) @result{} -6\n"
3156 #define FUNC_NAME s_scm_ash
3159 bits_to_shift
= scm_to_long (cnt
);
3161 if (SCM_I_INUMP (n
))
3163 scm_t_inum nn
= SCM_I_INUM (n
);
3165 if (bits_to_shift
> 0)
3167 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
3168 overflow a non-zero fixnum. For smaller shifts we check the
3169 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
3170 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
3171 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
3177 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
3179 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
3182 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
3186 SCM result
= scm_i_inum2big (nn
);
3187 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
3194 bits_to_shift
= -bits_to_shift
;
3195 if (bits_to_shift
>= SCM_LONG_BIT
)
3196 return (nn
>= 0 ? SCM_INUM0
: SCM_I_MAKINUM(-1));
3198 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
3202 else if (SCM_BIGP (n
))
3206 if (bits_to_shift
== 0)
3209 result
= scm_i_mkbig ();
3210 if (bits_to_shift
>= 0)
3212 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
3218 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
3219 we have to allocate a bignum even if the result is going to be a
3221 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
3223 return scm_i_normbig (result
);
3229 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3235 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
3236 (SCM n
, SCM start
, SCM end
),
3237 "Return the integer composed of the @var{start} (inclusive)\n"
3238 "through @var{end} (exclusive) bits of @var{n}. The\n"
3239 "@var{start}th bit becomes the 0-th bit in the result.\n"
3242 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
3243 " @result{} \"1010\"\n"
3244 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
3245 " @result{} \"10110\"\n"
3247 #define FUNC_NAME s_scm_bit_extract
3249 unsigned long int istart
, iend
, bits
;
3250 istart
= scm_to_ulong (start
);
3251 iend
= scm_to_ulong (end
);
3252 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
3254 /* how many bits to keep */
3255 bits
= iend
- istart
;
3257 if (SCM_I_INUMP (n
))
3259 scm_t_inum in
= SCM_I_INUM (n
);
3261 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
3262 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
3263 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
3265 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
3267 /* Since we emulate two's complement encoded numbers, this
3268 * special case requires us to produce a result that has
3269 * more bits than can be stored in a fixnum.
3271 SCM result
= scm_i_inum2big (in
);
3272 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
3277 /* mask down to requisite bits */
3278 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
3279 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
3281 else if (SCM_BIGP (n
))
3286 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
3290 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
3291 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
3292 such bits into a ulong. */
3293 result
= scm_i_mkbig ();
3294 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
3295 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
3296 result
= scm_i_normbig (result
);
3298 scm_remember_upto_here_1 (n
);
3302 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3307 static const char scm_logtab
[] = {
3308 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
3311 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
3313 "Return the number of bits in integer @var{n}. If integer is\n"
3314 "positive, the 1-bits in its binary representation are counted.\n"
3315 "If negative, the 0-bits in its two's-complement binary\n"
3316 "representation are counted. If 0, 0 is returned.\n"
3319 "(logcount #b10101010)\n"
3326 #define FUNC_NAME s_scm_logcount
3328 if (SCM_I_INUMP (n
))
3330 unsigned long c
= 0;
3331 scm_t_inum nn
= SCM_I_INUM (n
);
3336 c
+= scm_logtab
[15 & nn
];
3339 return SCM_I_MAKINUM (c
);
3341 else if (SCM_BIGP (n
))
3343 unsigned long count
;
3344 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
3345 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
3347 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
3348 scm_remember_upto_here_1 (n
);
3349 return SCM_I_MAKINUM (count
);
3352 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3357 static const char scm_ilentab
[] = {
3358 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
3362 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
3364 "Return the number of bits necessary to represent @var{n}.\n"
3367 "(integer-length #b10101010)\n"
3369 "(integer-length 0)\n"
3371 "(integer-length #b1111)\n"
3374 #define FUNC_NAME s_scm_integer_length
3376 if (SCM_I_INUMP (n
))
3378 unsigned long c
= 0;
3380 scm_t_inum nn
= SCM_I_INUM (n
);
3386 l
= scm_ilentab
[15 & nn
];
3389 return SCM_I_MAKINUM (c
- 4 + l
);
3391 else if (SCM_BIGP (n
))
3393 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
3394 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
3395 1 too big, so check for that and adjust. */
3396 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
3397 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
3398 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
3399 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
3401 scm_remember_upto_here_1 (n
);
3402 return SCM_I_MAKINUM (size
);
3405 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3409 /*** NUMBERS -> STRINGS ***/
3410 #define SCM_MAX_DBL_PREC 60
3411 #define SCM_MAX_DBL_RADIX 36
3413 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
3414 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
3415 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
3418 void init_dblprec(int *prec
, int radix
) {
3419 /* determine floating point precision by adding successively
3420 smaller increments to 1.0 until it is considered == 1.0 */
3421 double f
= ((double)1.0)/radix
;
3422 double fsum
= 1.0 + f
;
3427 if (++(*prec
) > SCM_MAX_DBL_PREC
)
3439 void init_fx_radix(double *fx_list
, int radix
)
3441 /* initialize a per-radix list of tolerances. When added
3442 to a number < 1.0, we can determine if we should raund
3443 up and quit converting a number to a string. */
3447 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
3448 fx_list
[i
] = (fx_list
[i
-1] / radix
);
3451 /* use this array as a way to generate a single digit */
3452 static const char number_chars
[] = "0123456789abcdefghijklmnopqrstuvwxyz";
3455 idbl2str (double f
, char *a
, int radix
)
3457 int efmt
, dpt
, d
, i
, wp
;
3459 #ifdef DBL_MIN_10_EXP
3462 #endif /* DBL_MIN_10_EXP */
3467 radix
> SCM_MAX_DBL_RADIX
)
3469 /* revert to existing behavior */
3473 wp
= scm_dblprec
[radix
-2];
3474 fx
= fx_per_radix
[radix
-2];
3478 #ifdef HAVE_COPYSIGN
3479 double sgn
= copysign (1.0, f
);
3484 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
3490 strcpy (a
, "-inf.0");
3492 strcpy (a
, "+inf.0");
3497 strcpy (a
, "+nan.0");
3507 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
3508 make-uniform-vector, from causing infinite loops. */
3509 /* just do the checking...if it passes, we do the conversion for our
3510 radix again below */
3517 if (exp_cpy
-- < DBL_MIN_10_EXP
)
3525 while (f_cpy
> 10.0)
3528 if (exp_cpy
++ > DBL_MAX_10_EXP
)
3549 if (f
+ fx
[wp
] >= radix
)
3556 /* adding 9999 makes this equivalent to abs(x) % 3 */
3557 dpt
= (exp
+ 9999) % 3;
3561 efmt
= (exp
< -3) || (exp
> wp
+ 2);
3583 a
[ch
++] = number_chars
[d
];
3586 if (f
+ fx
[wp
] >= 1.0)
3588 a
[ch
- 1] = number_chars
[d
+1];
3600 if ((dpt
> 4) && (exp
> 6))
3602 d
= (a
[0] == '-' ? 2 : 1);
3603 for (i
= ch
++; i
> d
; i
--)
3616 if (a
[ch
- 1] == '.')
3617 a
[ch
++] = '0'; /* trailing zero */
3626 for (i
= radix
; i
<= exp
; i
*= radix
);
3627 for (i
/= radix
; i
; i
/= radix
)
3629 a
[ch
++] = number_chars
[exp
/ i
];
3638 icmplx2str (double real
, double imag
, char *str
, int radix
)
3642 i
= idbl2str (real
, str
, radix
);
3645 /* Don't output a '+' for negative numbers or for Inf and
3646 NaN. They will provide their own sign. */
3647 if (0 <= imag
&& !isinf (imag
) && !isnan (imag
))
3649 i
+= idbl2str (imag
, &str
[i
], radix
);
3656 iflo2str (SCM flt
, char *str
, int radix
)
3659 if (SCM_REALP (flt
))
3660 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
3662 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
3667 /* convert a scm_t_intmax to a string (unterminated). returns the number of
3668 characters in the result.
3670 p is destination: worst case (base 2) is SCM_INTBUFLEN */
3672 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
3677 return scm_iuint2str (-num
, rad
, p
) + 1;
3680 return scm_iuint2str (num
, rad
, p
);
3683 /* convert a scm_t_intmax to a string (unterminated). returns the number of
3684 characters in the result.
3686 p is destination: worst case (base 2) is SCM_INTBUFLEN */
3688 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
3692 scm_t_uintmax n
= num
;
3694 if (rad
< 2 || rad
> 36)
3695 scm_out_of_range ("scm_iuint2str", scm_from_int (rad
));
3697 for (n
/= rad
; n
> 0; n
/= rad
)
3707 p
[i
] = number_chars
[d
];
3712 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
3714 "Return a string holding the external representation of the\n"
3715 "number @var{n} in the given @var{radix}. If @var{n} is\n"
3716 "inexact, a radix of 10 will be used.")
3717 #define FUNC_NAME s_scm_number_to_string
3721 if (SCM_UNBNDP (radix
))
3724 base
= scm_to_signed_integer (radix
, 2, 36);
3726 if (SCM_I_INUMP (n
))
3728 char num_buf
[SCM_INTBUFLEN
];
3729 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
3730 return scm_from_locale_stringn (num_buf
, length
);
3732 else if (SCM_BIGP (n
))
3734 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
3735 scm_remember_upto_here_1 (n
);
3736 return scm_take_locale_string (str
);
3738 else if (SCM_FRACTIONP (n
))
3740 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
3741 scm_from_locale_string ("/"),
3742 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
3744 else if (SCM_INEXACTP (n
))
3746 char num_buf
[FLOBUFLEN
];
3747 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
3750 SCM_WRONG_TYPE_ARG (1, n
);
3755 /* These print routines used to be stubbed here so that scm_repl.c
3756 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
3759 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3761 char num_buf
[FLOBUFLEN
];
3762 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
3767 scm_i_print_double (double val
, SCM port
)
3769 char num_buf
[FLOBUFLEN
];
3770 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
3774 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3777 char num_buf
[FLOBUFLEN
];
3778 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
3783 scm_i_print_complex (double real
, double imag
, SCM port
)
3785 char num_buf
[FLOBUFLEN
];
3786 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
3790 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3793 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
3794 scm_display (str
, port
);
3795 scm_remember_upto_here_1 (str
);
3800 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3802 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
3803 scm_remember_upto_here_1 (exp
);
3804 scm_lfwrite (str
, (size_t) strlen (str
), port
);
3808 /*** END nums->strs ***/
3811 /*** STRINGS -> NUMBERS ***/
3813 /* The following functions implement the conversion from strings to numbers.
3814 * The implementation somehow follows the grammar for numbers as it is given
3815 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
3816 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
3817 * points should be noted about the implementation:
3818 * * Each function keeps a local index variable 'idx' that points at the
3819 * current position within the parsed string. The global index is only
3820 * updated if the function could parse the corresponding syntactic unit
3822 * * Similarly, the functions keep track of indicators of inexactness ('#',
3823 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
3824 * global exactness information is only updated after each part has been
3825 * successfully parsed.
3826 * * Sequences of digits are parsed into temporary variables holding fixnums.
3827 * Only if these fixnums would overflow, the result variables are updated
3828 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
3829 * the temporary variables holding the fixnums are cleared, and the process
3830 * starts over again. If for example fixnums were able to store five decimal
3831 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
3832 * and the result was computed as 12345 * 100000 + 67890. In other words,
3833 * only every five digits two bignum operations were performed.
3836 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
3838 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
3840 /* Caller is responsible for checking that the return value is in range
3841 for the given radix, which should be <= 36. */
3843 char_decimal_value (scm_t_uint32 c
)
3845 /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
3846 that's certainly above any valid decimal, so we take advantage of
3847 that to elide some tests. */
3848 unsigned int d
= (unsigned int) uc_decimal_value (c
);
3850 /* If that failed, try extended hexadecimals, then. Only accept ascii
3855 if (c
>= (scm_t_uint32
) 'a')
3856 d
= c
- (scm_t_uint32
)'a' + 10U;
3862 mem2uinteger (SCM mem
, unsigned int *p_idx
,
3863 unsigned int radix
, enum t_exactness
*p_exactness
)
3865 unsigned int idx
= *p_idx
;
3866 unsigned int hash_seen
= 0;
3867 scm_t_bits shift
= 1;
3869 unsigned int digit_value
;
3872 size_t len
= scm_i_string_length (mem
);
3877 c
= scm_i_string_ref (mem
, idx
);
3878 digit_value
= char_decimal_value (c
);
3879 if (digit_value
>= radix
)
3883 result
= SCM_I_MAKINUM (digit_value
);
3886 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
3896 digit_value
= char_decimal_value (c
);
3897 /* This check catches non-decimals in addition to out-of-range
3899 if (digit_value
>= radix
)
3904 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
3906 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3908 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3915 shift
= shift
* radix
;
3916 add
= add
* radix
+ digit_value
;
3921 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3923 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3927 *p_exactness
= INEXACT
;
3933 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
3934 * covers the parts of the rules that start at a potential point. The value
3935 * of the digits up to the point have been parsed by the caller and are given
3936 * in variable result. The content of *p_exactness indicates, whether a hash
3937 * has already been seen in the digits before the point.
3940 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
3943 mem2decimal_from_point (SCM result
, SCM mem
,
3944 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
3946 unsigned int idx
= *p_idx
;
3947 enum t_exactness x
= *p_exactness
;
3948 size_t len
= scm_i_string_length (mem
);
3953 if (scm_i_string_ref (mem
, idx
) == '.')
3955 scm_t_bits shift
= 1;
3957 unsigned int digit_value
;
3958 SCM big_shift
= SCM_INUM1
;
3963 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
3964 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
3969 digit_value
= DIGIT2UINT (c
);
3980 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
3982 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
3983 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3985 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3993 add
= add
* 10 + digit_value
;
3999 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
4000 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
4001 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
4004 result
= scm_divide (result
, big_shift
);
4006 /* We've seen a decimal point, thus the value is implicitly inexact. */
4018 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
4020 switch (scm_i_string_ref (mem
, idx
))
4032 c
= scm_i_string_ref (mem
, idx
);
4040 c
= scm_i_string_ref (mem
, idx
);
4049 c
= scm_i_string_ref (mem
, idx
);
4054 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4058 exponent
= DIGIT2UINT (c
);
4061 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
4062 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4065 if (exponent
<= SCM_MAXEXP
)
4066 exponent
= exponent
* 10 + DIGIT2UINT (c
);
4072 if (exponent
> SCM_MAXEXP
)
4074 size_t exp_len
= idx
- start
;
4075 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
4076 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
4077 scm_out_of_range ("string->number", exp_num
);
4080 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
4082 result
= scm_product (result
, e
);
4084 result
= scm_divide2real (result
, e
);
4086 /* We've seen an exponent, thus the value is implicitly inexact. */
4104 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
4107 mem2ureal (SCM mem
, unsigned int *p_idx
,
4108 unsigned int radix
, enum t_exactness
*p_exactness
)
4110 unsigned int idx
= *p_idx
;
4112 size_t len
= scm_i_string_length (mem
);
4114 /* Start off believing that the number will be exact. This changes
4115 to INEXACT if we see a decimal point or a hash. */
4116 enum t_exactness x
= EXACT
;
4121 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
4127 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
4129 /* Cobble up the fractional part. We might want to set the
4130 NaN's mantissa from it. */
4132 mem2uinteger (mem
, &idx
, 10, &x
);
4137 if (scm_i_string_ref (mem
, idx
) == '.')
4141 else if (idx
+ 1 == len
)
4143 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
4146 result
= mem2decimal_from_point (SCM_INUM0
, mem
,
4153 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
4154 if (scm_is_false (uinteger
))
4159 else if (scm_i_string_ref (mem
, idx
) == '/')
4167 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
4168 if (scm_is_false (divisor
))
4171 /* both are int/big here, I assume */
4172 result
= scm_i_make_ratio (uinteger
, divisor
);
4174 else if (radix
== 10)
4176 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
4177 if (scm_is_false (result
))
4186 /* Update *p_exactness if the number just read was inexact. This is
4187 important for complex numbers, so that a complex number is
4188 treated as inexact overall if either its real or imaginary part
4194 /* When returning an inexact zero, make sure it is represented as a
4195 floating point value so that we can change its sign.
4197 if (scm_is_eq (result
, SCM_INUM0
) && *p_exactness
== INEXACT
)
4198 result
= scm_from_double (0.0);
4204 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
4207 mem2complex (SCM mem
, unsigned int idx
,
4208 unsigned int radix
, enum t_exactness
*p_exactness
)
4213 size_t len
= scm_i_string_length (mem
);
4218 c
= scm_i_string_ref (mem
, idx
);
4233 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
4234 if (scm_is_false (ureal
))
4236 /* input must be either +i or -i */
4241 if (scm_i_string_ref (mem
, idx
) == 'i'
4242 || scm_i_string_ref (mem
, idx
) == 'I')
4248 return scm_make_rectangular (SCM_INUM0
, SCM_I_MAKINUM (sign
));
4255 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
4256 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
4261 c
= scm_i_string_ref (mem
, idx
);
4265 /* either +<ureal>i or -<ureal>i */
4272 return scm_make_rectangular (SCM_INUM0
, ureal
);
4275 /* polar input: <real>@<real>. */
4286 c
= scm_i_string_ref (mem
, idx
);
4304 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
4305 if (scm_is_false (angle
))
4310 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
4311 angle
= scm_difference (angle
, SCM_UNDEFINED
);
4313 result
= scm_make_polar (ureal
, angle
);
4318 /* expecting input matching <real>[+-]<ureal>?i */
4325 int sign
= (c
== '+') ? 1 : -1;
4326 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
4328 if (scm_is_false (imag
))
4329 imag
= SCM_I_MAKINUM (sign
);
4330 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
4331 imag
= scm_difference (imag
, SCM_UNDEFINED
);
4335 if (scm_i_string_ref (mem
, idx
) != 'i'
4336 && scm_i_string_ref (mem
, idx
) != 'I')
4343 return scm_make_rectangular (ureal
, imag
);
4352 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
4354 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
4357 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
4359 unsigned int idx
= 0;
4360 unsigned int radix
= NO_RADIX
;
4361 enum t_exactness forced_x
= NO_EXACTNESS
;
4362 enum t_exactness implicit_x
= EXACT
;
4364 size_t len
= scm_i_string_length (mem
);
4366 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
4367 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
4369 switch (scm_i_string_ref (mem
, idx
+ 1))
4372 if (radix
!= NO_RADIX
)
4377 if (radix
!= NO_RADIX
)
4382 if (forced_x
!= NO_EXACTNESS
)
4387 if (forced_x
!= NO_EXACTNESS
)
4392 if (radix
!= NO_RADIX
)
4397 if (radix
!= NO_RADIX
)
4407 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
4408 if (radix
== NO_RADIX
)
4409 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
4411 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
4413 if (scm_is_false (result
))
4419 if (SCM_INEXACTP (result
))
4420 return scm_inexact_to_exact (result
);
4424 if (SCM_INEXACTP (result
))
4427 return scm_exact_to_inexact (result
);
4430 if (implicit_x
== INEXACT
)
4432 if (SCM_INEXACTP (result
))
4435 return scm_exact_to_inexact (result
);
4443 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
4444 unsigned int default_radix
)
4446 SCM str
= scm_from_locale_stringn (mem
, len
);
4448 return scm_i_string_to_number (str
, default_radix
);
4452 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
4453 (SCM string
, SCM radix
),
4454 "Return a number of the maximally precise representation\n"
4455 "expressed by the given @var{string}. @var{radix} must be an\n"
4456 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
4457 "is a default radix that may be overridden by an explicit radix\n"
4458 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
4459 "supplied, then the default radix is 10. If string is not a\n"
4460 "syntactically valid notation for a number, then\n"
4461 "@code{string->number} returns @code{#f}.")
4462 #define FUNC_NAME s_scm_string_to_number
4466 SCM_VALIDATE_STRING (1, string
);
4468 if (SCM_UNBNDP (radix
))
4471 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
4473 answer
= scm_i_string_to_number (string
, base
);
4474 scm_remember_upto_here_1 (string
);
4480 /*** END strs->nums ***/
4483 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
4485 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
4487 #define FUNC_NAME s_scm_number_p
4489 return scm_from_bool (SCM_NUMBERP (x
));
4493 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
4495 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
4496 "otherwise. Note that the sets of real, rational and integer\n"
4497 "values form subsets of the set of complex numbers, i. e. the\n"
4498 "predicate will also be fulfilled if @var{x} is a real,\n"
4499 "rational or integer number.")
4500 #define FUNC_NAME s_scm_complex_p
4502 /* all numbers are complex. */
4503 return scm_number_p (x
);
4507 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
4509 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
4510 "otherwise. Note that the set of integer values forms a subset of\n"
4511 "the set of real numbers, i. e. the predicate will also be\n"
4512 "fulfilled if @var{x} is an integer number.")
4513 #define FUNC_NAME s_scm_real_p
4515 return scm_from_bool
4516 (SCM_I_INUMP (x
) || SCM_REALP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
));
4520 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
4522 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
4523 "otherwise. Note that the set of integer values forms a subset of\n"
4524 "the set of rational numbers, i. e. the predicate will also be\n"
4525 "fulfilled if @var{x} is an integer number.")
4526 #define FUNC_NAME s_scm_rational_p
4528 if (SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
))
4530 else if (SCM_REALP (x
))
4531 /* due to their limited precision, finite floating point numbers are
4532 rational as well. (finite means neither infinity nor a NaN) */
4533 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
4539 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
4541 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
4543 #define FUNC_NAME s_scm_integer_p
4545 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4547 else if (SCM_REALP (x
))
4549 double val
= SCM_REAL_VALUE (x
);
4550 return scm_from_bool (!isinf (val
) && (val
== floor (val
)));
4558 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
4559 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
4560 (SCM x
, SCM y
, SCM rest
),
4561 "Return @code{#t} if all parameters are numerically equal.")
4562 #define FUNC_NAME s_scm_i_num_eq_p
4564 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4566 while (!scm_is_null (rest
))
4568 if (scm_is_false (scm_num_eq_p (x
, y
)))
4572 rest
= scm_cdr (rest
);
4574 return scm_num_eq_p (x
, y
);
4578 scm_num_eq_p (SCM x
, SCM y
)
4581 if (SCM_I_INUMP (x
))
4583 scm_t_signed_bits xx
= SCM_I_INUM (x
);
4584 if (SCM_I_INUMP (y
))
4586 scm_t_signed_bits yy
= SCM_I_INUM (y
);
4587 return scm_from_bool (xx
== yy
);
4589 else if (SCM_BIGP (y
))
4591 else if (SCM_REALP (y
))
4593 /* On a 32-bit system an inum fits a double, we can cast the inum
4594 to a double and compare.
4596 But on a 64-bit system an inum is bigger than a double and
4597 casting it to a double (call that dxx) will round. dxx is at
4598 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
4599 an integer and fits a long. So we cast yy to a long and
4600 compare with plain xx.
4602 An alternative (for any size system actually) would be to check
4603 yy is an integer (with floor) and is in range of an inum
4604 (compare against appropriate powers of 2) then test
4605 xx==(scm_t_signed_bits)yy. It's just a matter of which
4606 casts/comparisons might be fastest or easiest for the cpu. */
4608 double yy
= SCM_REAL_VALUE (y
);
4609 return scm_from_bool ((double) xx
== yy
4610 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
4611 || xx
== (scm_t_signed_bits
) yy
));
4613 else if (SCM_COMPLEXP (y
))
4614 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
4615 && (0.0 == SCM_COMPLEX_IMAG (y
)));
4616 else if (SCM_FRACTIONP (y
))
4619 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4621 else if (SCM_BIGP (x
))
4623 if (SCM_I_INUMP (y
))
4625 else if (SCM_BIGP (y
))
4627 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4628 scm_remember_upto_here_2 (x
, y
);
4629 return scm_from_bool (0 == cmp
);
4631 else if (SCM_REALP (y
))
4634 if (isnan (SCM_REAL_VALUE (y
)))
4636 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
4637 scm_remember_upto_here_1 (x
);
4638 return scm_from_bool (0 == cmp
);
4640 else if (SCM_COMPLEXP (y
))
4643 if (0.0 != SCM_COMPLEX_IMAG (y
))
4645 if (isnan (SCM_COMPLEX_REAL (y
)))
4647 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
4648 scm_remember_upto_here_1 (x
);
4649 return scm_from_bool (0 == cmp
);
4651 else if (SCM_FRACTIONP (y
))
4654 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4656 else if (SCM_REALP (x
))
4658 double xx
= SCM_REAL_VALUE (x
);
4659 if (SCM_I_INUMP (y
))
4661 /* see comments with inum/real above */
4662 scm_t_signed_bits yy
= SCM_I_INUM (y
);
4663 return scm_from_bool (xx
== (double) yy
4664 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
4665 || (scm_t_signed_bits
) xx
== yy
));
4667 else if (SCM_BIGP (y
))
4670 if (isnan (SCM_REAL_VALUE (x
)))
4672 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
4673 scm_remember_upto_here_1 (y
);
4674 return scm_from_bool (0 == cmp
);
4676 else if (SCM_REALP (y
))
4677 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
4678 else if (SCM_COMPLEXP (y
))
4679 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
4680 && (0.0 == SCM_COMPLEX_IMAG (y
)));
4681 else if (SCM_FRACTIONP (y
))
4683 double xx
= SCM_REAL_VALUE (x
);
4687 return scm_from_bool (xx
< 0.0);
4688 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4692 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4694 else if (SCM_COMPLEXP (x
))
4696 if (SCM_I_INUMP (y
))
4697 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
4698 && (SCM_COMPLEX_IMAG (x
) == 0.0));
4699 else if (SCM_BIGP (y
))
4702 if (0.0 != SCM_COMPLEX_IMAG (x
))
4704 if (isnan (SCM_COMPLEX_REAL (x
)))
4706 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
4707 scm_remember_upto_here_1 (y
);
4708 return scm_from_bool (0 == cmp
);
4710 else if (SCM_REALP (y
))
4711 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
4712 && (SCM_COMPLEX_IMAG (x
) == 0.0));
4713 else if (SCM_COMPLEXP (y
))
4714 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
4715 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
4716 else if (SCM_FRACTIONP (y
))
4719 if (SCM_COMPLEX_IMAG (x
) != 0.0)
4721 xx
= SCM_COMPLEX_REAL (x
);
4725 return scm_from_bool (xx
< 0.0);
4726 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4730 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4732 else if (SCM_FRACTIONP (x
))
4734 if (SCM_I_INUMP (y
))
4736 else if (SCM_BIGP (y
))
4738 else if (SCM_REALP (y
))
4740 double yy
= SCM_REAL_VALUE (y
);
4744 return scm_from_bool (0.0 < yy
);
4745 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4748 else if (SCM_COMPLEXP (y
))
4751 if (SCM_COMPLEX_IMAG (y
) != 0.0)
4753 yy
= SCM_COMPLEX_REAL (y
);
4757 return scm_from_bool (0.0 < yy
);
4758 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4761 else if (SCM_FRACTIONP (y
))
4762 return scm_i_fraction_equalp (x
, y
);
4764 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4767 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
4771 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
4772 done are good for inums, but for bignums an answer can almost always be
4773 had by just examining a few high bits of the operands, as done by GMP in
4774 mpq_cmp. flonum/frac compares likewise, but with the slight complication
4775 of the float exponent to take into account. */
4777 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
4778 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
4779 (SCM x
, SCM y
, SCM rest
),
4780 "Return @code{#t} if the list of parameters is monotonically\n"
4782 #define FUNC_NAME s_scm_i_num_less_p
4784 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4786 while (!scm_is_null (rest
))
4788 if (scm_is_false (scm_less_p (x
, y
)))
4792 rest
= scm_cdr (rest
);
4794 return scm_less_p (x
, y
);
4798 scm_less_p (SCM x
, SCM y
)
4801 if (SCM_I_INUMP (x
))
4803 scm_t_inum xx
= SCM_I_INUM (x
);
4804 if (SCM_I_INUMP (y
))
4806 scm_t_inum yy
= SCM_I_INUM (y
);
4807 return scm_from_bool (xx
< yy
);
4809 else if (SCM_BIGP (y
))
4811 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4812 scm_remember_upto_here_1 (y
);
4813 return scm_from_bool (sgn
> 0);
4815 else if (SCM_REALP (y
))
4816 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
4817 else if (SCM_FRACTIONP (y
))
4819 /* "x < a/b" becomes "x*b < a" */
4821 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
4822 y
= SCM_FRACTION_NUMERATOR (y
);
4826 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4828 else if (SCM_BIGP (x
))
4830 if (SCM_I_INUMP (y
))
4832 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4833 scm_remember_upto_here_1 (x
);
4834 return scm_from_bool (sgn
< 0);
4836 else if (SCM_BIGP (y
))
4838 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4839 scm_remember_upto_here_2 (x
, y
);
4840 return scm_from_bool (cmp
< 0);
4842 else if (SCM_REALP (y
))
4845 if (isnan (SCM_REAL_VALUE (y
)))
4847 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
4848 scm_remember_upto_here_1 (x
);
4849 return scm_from_bool (cmp
< 0);
4851 else if (SCM_FRACTIONP (y
))
4854 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4856 else if (SCM_REALP (x
))
4858 if (SCM_I_INUMP (y
))
4859 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
4860 else if (SCM_BIGP (y
))
4863 if (isnan (SCM_REAL_VALUE (x
)))
4865 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
4866 scm_remember_upto_here_1 (y
);
4867 return scm_from_bool (cmp
> 0);
4869 else if (SCM_REALP (y
))
4870 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
4871 else if (SCM_FRACTIONP (y
))
4873 double xx
= SCM_REAL_VALUE (x
);
4877 return scm_from_bool (xx
< 0.0);
4878 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4882 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4884 else if (SCM_FRACTIONP (x
))
4886 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
4888 /* "a/b < y" becomes "a < y*b" */
4889 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
4890 x
= SCM_FRACTION_NUMERATOR (x
);
4893 else if (SCM_REALP (y
))
4895 double yy
= SCM_REAL_VALUE (y
);
4899 return scm_from_bool (0.0 < yy
);
4900 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4903 else if (SCM_FRACTIONP (y
))
4905 /* "a/b < c/d" becomes "a*d < c*b" */
4906 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
4907 SCM_FRACTION_DENOMINATOR (y
));
4908 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
4909 SCM_FRACTION_DENOMINATOR (x
));
4915 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4918 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
4922 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
4923 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
4924 (SCM x
, SCM y
, SCM rest
),
4925 "Return @code{#t} if the list of parameters is monotonically\n"
4927 #define FUNC_NAME s_scm_i_num_gr_p
4929 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4931 while (!scm_is_null (rest
))
4933 if (scm_is_false (scm_gr_p (x
, y
)))
4937 rest
= scm_cdr (rest
);
4939 return scm_gr_p (x
, y
);
4942 #define FUNC_NAME s_scm_i_num_gr_p
4944 scm_gr_p (SCM x
, SCM y
)
4946 if (!SCM_NUMBERP (x
))
4947 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
4948 else if (!SCM_NUMBERP (y
))
4949 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
4951 return scm_less_p (y
, x
);
4956 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
4957 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
4958 (SCM x
, SCM y
, SCM rest
),
4959 "Return @code{#t} if the list of parameters is monotonically\n"
4961 #define FUNC_NAME s_scm_i_num_leq_p
4963 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4965 while (!scm_is_null (rest
))
4967 if (scm_is_false (scm_leq_p (x
, y
)))
4971 rest
= scm_cdr (rest
);
4973 return scm_leq_p (x
, y
);
4976 #define FUNC_NAME s_scm_i_num_leq_p
4978 scm_leq_p (SCM x
, SCM y
)
4980 if (!SCM_NUMBERP (x
))
4981 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
4982 else if (!SCM_NUMBERP (y
))
4983 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
4984 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
4987 return scm_not (scm_less_p (y
, x
));
4992 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
4993 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
4994 (SCM x
, SCM y
, SCM rest
),
4995 "Return @code{#t} if the list of parameters is monotonically\n"
4997 #define FUNC_NAME s_scm_i_num_geq_p
4999 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
5001 while (!scm_is_null (rest
))
5003 if (scm_is_false (scm_geq_p (x
, y
)))
5007 rest
= scm_cdr (rest
);
5009 return scm_geq_p (x
, y
);
5012 #define FUNC_NAME s_scm_i_num_geq_p
5014 scm_geq_p (SCM x
, SCM y
)
5016 if (!SCM_NUMBERP (x
))
5017 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
5018 else if (!SCM_NUMBERP (y
))
5019 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
5020 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
5023 return scm_not (scm_less_p (x
, y
));
5028 SCM_PRIMITIVE_GENERIC (scm_zero_p
, "zero?", 1, 0, 0,
5030 "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
5032 #define FUNC_NAME s_scm_zero_p
5034 if (SCM_I_INUMP (z
))
5035 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
5036 else if (SCM_BIGP (z
))
5038 else if (SCM_REALP (z
))
5039 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
5040 else if (SCM_COMPLEXP (z
))
5041 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
5042 && SCM_COMPLEX_IMAG (z
) == 0.0);
5043 else if (SCM_FRACTIONP (z
))
5046 SCM_WTA_DISPATCH_1 (g_scm_zero_p
, z
, SCM_ARG1
, s_scm_zero_p
);
5051 SCM_PRIMITIVE_GENERIC (scm_positive_p
, "positive?", 1, 0, 0,
5053 "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
5055 #define FUNC_NAME s_scm_positive_p
5057 if (SCM_I_INUMP (x
))
5058 return scm_from_bool (SCM_I_INUM (x
) > 0);
5059 else if (SCM_BIGP (x
))
5061 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5062 scm_remember_upto_here_1 (x
);
5063 return scm_from_bool (sgn
> 0);
5065 else if (SCM_REALP (x
))
5066 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
5067 else if (SCM_FRACTIONP (x
))
5068 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
5070 SCM_WTA_DISPATCH_1 (g_scm_positive_p
, x
, SCM_ARG1
, s_scm_positive_p
);
5075 SCM_PRIMITIVE_GENERIC (scm_negative_p
, "negative?", 1, 0, 0,
5077 "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
5079 #define FUNC_NAME s_scm_negative_p
5081 if (SCM_I_INUMP (x
))
5082 return scm_from_bool (SCM_I_INUM (x
) < 0);
5083 else if (SCM_BIGP (x
))
5085 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5086 scm_remember_upto_here_1 (x
);
5087 return scm_from_bool (sgn
< 0);
5089 else if (SCM_REALP (x
))
5090 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
5091 else if (SCM_FRACTIONP (x
))
5092 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
5094 SCM_WTA_DISPATCH_1 (g_scm_negative_p
, x
, SCM_ARG1
, s_scm_negative_p
);
5099 /* scm_min and scm_max return an inexact when either argument is inexact, as
5100 required by r5rs. On that basis, for exact/inexact combinations the
5101 exact is converted to inexact to compare and possibly return. This is
5102 unlike scm_less_p above which takes some trouble to preserve all bits in
5103 its test, such trouble is not required for min and max. */
5105 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
5106 (SCM x
, SCM y
, SCM rest
),
5107 "Return the maximum of all parameter values.")
5108 #define FUNC_NAME s_scm_i_max
5110 while (!scm_is_null (rest
))
5111 { x
= scm_max (x
, y
);
5113 rest
= scm_cdr (rest
);
5115 return scm_max (x
, y
);
5119 #define s_max s_scm_i_max
5120 #define g_max g_scm_i_max
5123 scm_max (SCM x
, SCM y
)
5128 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
5129 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
5132 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
5135 if (SCM_I_INUMP (x
))
5137 scm_t_inum xx
= SCM_I_INUM (x
);
5138 if (SCM_I_INUMP (y
))
5140 scm_t_inum yy
= SCM_I_INUM (y
);
5141 return (xx
< yy
) ? y
: x
;
5143 else if (SCM_BIGP (y
))
5145 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5146 scm_remember_upto_here_1 (y
);
5147 return (sgn
< 0) ? x
: y
;
5149 else if (SCM_REALP (y
))
5152 /* if y==NaN then ">" is false and we return NaN */
5153 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
5155 else if (SCM_FRACTIONP (y
))
5158 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
5161 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5163 else if (SCM_BIGP (x
))
5165 if (SCM_I_INUMP (y
))
5167 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5168 scm_remember_upto_here_1 (x
);
5169 return (sgn
< 0) ? y
: x
;
5171 else if (SCM_BIGP (y
))
5173 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
5174 scm_remember_upto_here_2 (x
, y
);
5175 return (cmp
> 0) ? x
: y
;
5177 else if (SCM_REALP (y
))
5179 /* if y==NaN then xx>yy is false, so we return the NaN y */
5182 xx
= scm_i_big2dbl (x
);
5183 yy
= SCM_REAL_VALUE (y
);
5184 return (xx
> yy
? scm_from_double (xx
) : y
);
5186 else if (SCM_FRACTIONP (y
))
5191 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5193 else if (SCM_REALP (x
))
5195 if (SCM_I_INUMP (y
))
5197 double z
= SCM_I_INUM (y
);
5198 /* if x==NaN then "<" is false and we return NaN */
5199 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
5201 else if (SCM_BIGP (y
))
5206 else if (SCM_REALP (y
))
5208 /* if x==NaN then our explicit check means we return NaN
5209 if y==NaN then ">" is false and we return NaN
5210 calling isnan is unavoidable, since it's the only way to know
5211 which of x or y causes any compares to be false */
5212 double xx
= SCM_REAL_VALUE (x
);
5213 return (isnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
5215 else if (SCM_FRACTIONP (y
))
5217 double yy
= scm_i_fraction2double (y
);
5218 double xx
= SCM_REAL_VALUE (x
);
5219 return (xx
< yy
) ? scm_from_double (yy
) : x
;
5222 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5224 else if (SCM_FRACTIONP (x
))
5226 if (SCM_I_INUMP (y
))
5230 else if (SCM_BIGP (y
))
5234 else if (SCM_REALP (y
))
5236 double xx
= scm_i_fraction2double (x
);
5237 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
5239 else if (SCM_FRACTIONP (y
))
5244 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5247 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
5251 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
5252 (SCM x
, SCM y
, SCM rest
),
5253 "Return the minimum of all parameter values.")
5254 #define FUNC_NAME s_scm_i_min
5256 while (!scm_is_null (rest
))
5257 { x
= scm_min (x
, y
);
5259 rest
= scm_cdr (rest
);
5261 return scm_min (x
, y
);
5265 #define s_min s_scm_i_min
5266 #define g_min g_scm_i_min
5269 scm_min (SCM x
, SCM y
)
5274 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
5275 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
5278 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
5281 if (SCM_I_INUMP (x
))
5283 scm_t_inum xx
= SCM_I_INUM (x
);
5284 if (SCM_I_INUMP (y
))
5286 scm_t_inum yy
= SCM_I_INUM (y
);
5287 return (xx
< yy
) ? x
: y
;
5289 else if (SCM_BIGP (y
))
5291 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5292 scm_remember_upto_here_1 (y
);
5293 return (sgn
< 0) ? y
: x
;
5295 else if (SCM_REALP (y
))
5298 /* if y==NaN then "<" is false and we return NaN */
5299 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
5301 else if (SCM_FRACTIONP (y
))
5304 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
5307 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5309 else if (SCM_BIGP (x
))
5311 if (SCM_I_INUMP (y
))
5313 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5314 scm_remember_upto_here_1 (x
);
5315 return (sgn
< 0) ? x
: y
;
5317 else if (SCM_BIGP (y
))
5319 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
5320 scm_remember_upto_here_2 (x
, y
);
5321 return (cmp
> 0) ? y
: x
;
5323 else if (SCM_REALP (y
))
5325 /* if y==NaN then xx<yy is false, so we return the NaN y */
5328 xx
= scm_i_big2dbl (x
);
5329 yy
= SCM_REAL_VALUE (y
);
5330 return (xx
< yy
? scm_from_double (xx
) : y
);
5332 else if (SCM_FRACTIONP (y
))
5337 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5339 else if (SCM_REALP (x
))
5341 if (SCM_I_INUMP (y
))
5343 double z
= SCM_I_INUM (y
);
5344 /* if x==NaN then "<" is false and we return NaN */
5345 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
5347 else if (SCM_BIGP (y
))
5352 else if (SCM_REALP (y
))
5354 /* if x==NaN then our explicit check means we return NaN
5355 if y==NaN then "<" is false and we return NaN
5356 calling isnan is unavoidable, since it's the only way to know
5357 which of x or y causes any compares to be false */
5358 double xx
= SCM_REAL_VALUE (x
);
5359 return (isnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
5361 else if (SCM_FRACTIONP (y
))
5363 double yy
= scm_i_fraction2double (y
);
5364 double xx
= SCM_REAL_VALUE (x
);
5365 return (yy
< xx
) ? scm_from_double (yy
) : x
;
5368 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5370 else if (SCM_FRACTIONP (x
))
5372 if (SCM_I_INUMP (y
))
5376 else if (SCM_BIGP (y
))
5380 else if (SCM_REALP (y
))
5382 double xx
= scm_i_fraction2double (x
);
5383 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
5385 else if (SCM_FRACTIONP (y
))
5390 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5393 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
5397 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
5398 (SCM x
, SCM y
, SCM rest
),
5399 "Return the sum of all parameter values. Return 0 if called without\n"
5401 #define FUNC_NAME s_scm_i_sum
5403 while (!scm_is_null (rest
))
5404 { x
= scm_sum (x
, y
);
5406 rest
= scm_cdr (rest
);
5408 return scm_sum (x
, y
);
5412 #define s_sum s_scm_i_sum
5413 #define g_sum g_scm_i_sum
5416 scm_sum (SCM x
, SCM y
)
5418 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5420 if (SCM_NUMBERP (x
)) return x
;
5421 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
5422 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
5425 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5427 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5429 scm_t_inum xx
= SCM_I_INUM (x
);
5430 scm_t_inum yy
= SCM_I_INUM (y
);
5431 scm_t_inum z
= xx
+ yy
;
5432 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_inum2big (z
);
5434 else if (SCM_BIGP (y
))
5439 else if (SCM_REALP (y
))
5441 scm_t_inum xx
= SCM_I_INUM (x
);
5442 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
5444 else if (SCM_COMPLEXP (y
))
5446 scm_t_inum xx
= SCM_I_INUM (x
);
5447 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
5448 SCM_COMPLEX_IMAG (y
));
5450 else if (SCM_FRACTIONP (y
))
5451 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
5452 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
5453 SCM_FRACTION_DENOMINATOR (y
));
5455 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5456 } else if (SCM_BIGP (x
))
5458 if (SCM_I_INUMP (y
))
5463 inum
= SCM_I_INUM (y
);
5466 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5469 SCM result
= scm_i_mkbig ();
5470 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
5471 scm_remember_upto_here_1 (x
);
5472 /* we know the result will have to be a bignum */
5475 return scm_i_normbig (result
);
5479 SCM result
= scm_i_mkbig ();
5480 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
5481 scm_remember_upto_here_1 (x
);
5482 /* we know the result will have to be a bignum */
5485 return scm_i_normbig (result
);
5488 else if (SCM_BIGP (y
))
5490 SCM result
= scm_i_mkbig ();
5491 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5492 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5493 mpz_add (SCM_I_BIG_MPZ (result
),
5496 scm_remember_upto_here_2 (x
, y
);
5497 /* we know the result will have to be a bignum */
5500 return scm_i_normbig (result
);
5502 else if (SCM_REALP (y
))
5504 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
5505 scm_remember_upto_here_1 (x
);
5506 return scm_from_double (result
);
5508 else if (SCM_COMPLEXP (y
))
5510 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
5511 + SCM_COMPLEX_REAL (y
));
5512 scm_remember_upto_here_1 (x
);
5513 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
5515 else if (SCM_FRACTIONP (y
))
5516 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
5517 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
5518 SCM_FRACTION_DENOMINATOR (y
));
5520 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5522 else if (SCM_REALP (x
))
5524 if (SCM_I_INUMP (y
))
5525 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
5526 else if (SCM_BIGP (y
))
5528 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
5529 scm_remember_upto_here_1 (y
);
5530 return scm_from_double (result
);
5532 else if (SCM_REALP (y
))
5533 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
5534 else if (SCM_COMPLEXP (y
))
5535 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
5536 SCM_COMPLEX_IMAG (y
));
5537 else if (SCM_FRACTIONP (y
))
5538 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
5540 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5542 else if (SCM_COMPLEXP (x
))
5544 if (SCM_I_INUMP (y
))
5545 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
5546 SCM_COMPLEX_IMAG (x
));
5547 else if (SCM_BIGP (y
))
5549 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
5550 + SCM_COMPLEX_REAL (x
));
5551 scm_remember_upto_here_1 (y
);
5552 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
5554 else if (SCM_REALP (y
))
5555 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
5556 SCM_COMPLEX_IMAG (x
));
5557 else if (SCM_COMPLEXP (y
))
5558 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
5559 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
5560 else if (SCM_FRACTIONP (y
))
5561 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
5562 SCM_COMPLEX_IMAG (x
));
5564 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5566 else if (SCM_FRACTIONP (x
))
5568 if (SCM_I_INUMP (y
))
5569 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
5570 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
5571 SCM_FRACTION_DENOMINATOR (x
));
5572 else if (SCM_BIGP (y
))
5573 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
5574 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
5575 SCM_FRACTION_DENOMINATOR (x
));
5576 else if (SCM_REALP (y
))
5577 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
5578 else if (SCM_COMPLEXP (y
))
5579 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
5580 SCM_COMPLEX_IMAG (y
));
5581 else if (SCM_FRACTIONP (y
))
5582 /* a/b + c/d = (ad + bc) / bd */
5583 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5584 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
5585 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
5587 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5590 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
5594 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
5596 "Return @math{@var{x}+1}.")
5597 #define FUNC_NAME s_scm_oneplus
5599 return scm_sum (x
, SCM_INUM1
);
5604 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
5605 (SCM x
, SCM y
, SCM rest
),
5606 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
5607 "the sum of all but the first argument are subtracted from the first\n"
5609 #define FUNC_NAME s_scm_i_difference
5611 while (!scm_is_null (rest
))
5612 { x
= scm_difference (x
, y
);
5614 rest
= scm_cdr (rest
);
5616 return scm_difference (x
, y
);
5620 #define s_difference s_scm_i_difference
5621 #define g_difference g_scm_i_difference
5624 scm_difference (SCM x
, SCM y
)
5625 #define FUNC_NAME s_difference
5627 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5630 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
5632 if (SCM_I_INUMP (x
))
5634 scm_t_inum xx
= -SCM_I_INUM (x
);
5635 if (SCM_FIXABLE (xx
))
5636 return SCM_I_MAKINUM (xx
);
5638 return scm_i_inum2big (xx
);
5640 else if (SCM_BIGP (x
))
5641 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
5642 bignum, but negating that gives a fixnum. */
5643 return scm_i_normbig (scm_i_clonebig (x
, 0));
5644 else if (SCM_REALP (x
))
5645 return scm_from_double (-SCM_REAL_VALUE (x
));
5646 else if (SCM_COMPLEXP (x
))
5647 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
5648 -SCM_COMPLEX_IMAG (x
));
5649 else if (SCM_FRACTIONP (x
))
5650 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
5651 SCM_FRACTION_DENOMINATOR (x
));
5653 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
5656 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5658 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5660 scm_t_inum xx
= SCM_I_INUM (x
);
5661 scm_t_inum yy
= SCM_I_INUM (y
);
5662 scm_t_inum z
= xx
- yy
;
5663 if (SCM_FIXABLE (z
))
5664 return SCM_I_MAKINUM (z
);
5666 return scm_i_inum2big (z
);
5668 else if (SCM_BIGP (y
))
5670 /* inum-x - big-y */
5671 scm_t_inum xx
= SCM_I_INUM (x
);
5675 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
5676 bignum, but negating that gives a fixnum. */
5677 return scm_i_normbig (scm_i_clonebig (y
, 0));
5681 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5682 SCM result
= scm_i_mkbig ();
5685 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
5688 /* x - y == -(y + -x) */
5689 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
5690 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5692 scm_remember_upto_here_1 (y
);
5694 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
5695 /* we know the result will have to be a bignum */
5698 return scm_i_normbig (result
);
5701 else if (SCM_REALP (y
))
5703 scm_t_inum xx
= SCM_I_INUM (x
);
5704 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
5706 else if (SCM_COMPLEXP (y
))
5708 scm_t_inum xx
= SCM_I_INUM (x
);
5709 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
5710 - SCM_COMPLEX_IMAG (y
));
5712 else if (SCM_FRACTIONP (y
))
5713 /* a - b/c = (ac - b) / c */
5714 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5715 SCM_FRACTION_NUMERATOR (y
)),
5716 SCM_FRACTION_DENOMINATOR (y
));
5718 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5720 else if (SCM_BIGP (x
))
5722 if (SCM_I_INUMP (y
))
5724 /* big-x - inum-y */
5725 scm_t_inum yy
= SCM_I_INUM (y
);
5726 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5728 scm_remember_upto_here_1 (x
);
5730 return (SCM_FIXABLE (-yy
) ?
5731 SCM_I_MAKINUM (-yy
) : scm_from_inum (-yy
));
5734 SCM result
= scm_i_mkbig ();
5737 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
5739 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
5740 scm_remember_upto_here_1 (x
);
5742 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
5743 /* we know the result will have to be a bignum */
5746 return scm_i_normbig (result
);
5749 else if (SCM_BIGP (y
))
5751 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5752 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5753 SCM result
= scm_i_mkbig ();
5754 mpz_sub (SCM_I_BIG_MPZ (result
),
5757 scm_remember_upto_here_2 (x
, y
);
5758 /* we know the result will have to be a bignum */
5759 if ((sgn_x
== 1) && (sgn_y
== -1))
5761 if ((sgn_x
== -1) && (sgn_y
== 1))
5763 return scm_i_normbig (result
);
5765 else if (SCM_REALP (y
))
5767 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
5768 scm_remember_upto_here_1 (x
);
5769 return scm_from_double (result
);
5771 else if (SCM_COMPLEXP (y
))
5773 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
5774 - SCM_COMPLEX_REAL (y
));
5775 scm_remember_upto_here_1 (x
);
5776 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
5778 else if (SCM_FRACTIONP (y
))
5779 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5780 SCM_FRACTION_NUMERATOR (y
)),
5781 SCM_FRACTION_DENOMINATOR (y
));
5782 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5784 else if (SCM_REALP (x
))
5786 if (SCM_I_INUMP (y
))
5787 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
5788 else if (SCM_BIGP (y
))
5790 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
5791 scm_remember_upto_here_1 (x
);
5792 return scm_from_double (result
);
5794 else if (SCM_REALP (y
))
5795 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
5796 else if (SCM_COMPLEXP (y
))
5797 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
5798 -SCM_COMPLEX_IMAG (y
));
5799 else if (SCM_FRACTIONP (y
))
5800 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
5802 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5804 else if (SCM_COMPLEXP (x
))
5806 if (SCM_I_INUMP (y
))
5807 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
5808 SCM_COMPLEX_IMAG (x
));
5809 else if (SCM_BIGP (y
))
5811 double real_part
= (SCM_COMPLEX_REAL (x
)
5812 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
5813 scm_remember_upto_here_1 (x
);
5814 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
5816 else if (SCM_REALP (y
))
5817 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
5818 SCM_COMPLEX_IMAG (x
));
5819 else if (SCM_COMPLEXP (y
))
5820 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
5821 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
5822 else if (SCM_FRACTIONP (y
))
5823 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
5824 SCM_COMPLEX_IMAG (x
));
5826 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5828 else if (SCM_FRACTIONP (x
))
5830 if (SCM_I_INUMP (y
))
5831 /* a/b - c = (a - cb) / b */
5832 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
5833 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
5834 SCM_FRACTION_DENOMINATOR (x
));
5835 else if (SCM_BIGP (y
))
5836 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
5837 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
5838 SCM_FRACTION_DENOMINATOR (x
));
5839 else if (SCM_REALP (y
))
5840 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
5841 else if (SCM_COMPLEXP (y
))
5842 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
5843 -SCM_COMPLEX_IMAG (y
));
5844 else if (SCM_FRACTIONP (y
))
5845 /* a/b - c/d = (ad - bc) / bd */
5846 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5847 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
5848 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
5850 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5853 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
5858 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
5860 "Return @math{@var{x}-1}.")
5861 #define FUNC_NAME s_scm_oneminus
5863 return scm_difference (x
, SCM_INUM1
);
5868 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
5869 (SCM x
, SCM y
, SCM rest
),
5870 "Return the product of all arguments. If called without arguments,\n"
5872 #define FUNC_NAME s_scm_i_product
5874 while (!scm_is_null (rest
))
5875 { x
= scm_product (x
, y
);
5877 rest
= scm_cdr (rest
);
5879 return scm_product (x
, y
);
5883 #define s_product s_scm_i_product
5884 #define g_product g_scm_i_product
5887 scm_product (SCM x
, SCM y
)
5889 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5892 return SCM_I_MAKINUM (1L);
5893 else if (SCM_NUMBERP (x
))
5896 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
5899 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5904 xx
= SCM_I_INUM (x
);
5908 case 0: return x
; break;
5909 case 1: return y
; break;
5911 * The following case (x = -1) is important for more than
5912 * just optimization. It handles the case of negating
5913 * (+ 1 most-positive-fixnum) aka (- most-negative-fixnum),
5914 * which is a bignum that must be changed back into a fixnum.
5915 * Failure to do so will cause the following to return #f:
5916 * (= most-negative-fixnum (* -1 (- most-negative-fixnum)))
5919 return scm_difference(y
, SCM_UNDEFINED
);
5923 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5925 scm_t_inum yy
= SCM_I_INUM (y
);
5926 scm_t_inum kk
= xx
* yy
;
5927 SCM k
= SCM_I_MAKINUM (kk
);
5928 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
5932 SCM result
= scm_i_inum2big (xx
);
5933 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
5934 return scm_i_normbig (result
);
5937 else if (SCM_BIGP (y
))
5939 SCM result
= scm_i_mkbig ();
5940 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
5941 scm_remember_upto_here_1 (y
);
5944 else if (SCM_REALP (y
))
5945 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
5946 else if (SCM_COMPLEXP (y
))
5947 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
5948 xx
* SCM_COMPLEX_IMAG (y
));
5949 else if (SCM_FRACTIONP (y
))
5950 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
5951 SCM_FRACTION_DENOMINATOR (y
));
5953 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
5955 else if (SCM_BIGP (x
))
5957 if (SCM_I_INUMP (y
))
5962 else if (SCM_BIGP (y
))
5964 SCM result
= scm_i_mkbig ();
5965 mpz_mul (SCM_I_BIG_MPZ (result
),
5968 scm_remember_upto_here_2 (x
, y
);
5971 else if (SCM_REALP (y
))
5973 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
5974 scm_remember_upto_here_1 (x
);
5975 return scm_from_double (result
);
5977 else if (SCM_COMPLEXP (y
))
5979 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
5980 scm_remember_upto_here_1 (x
);
5981 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
5982 z
* SCM_COMPLEX_IMAG (y
));
5984 else if (SCM_FRACTIONP (y
))
5985 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
5986 SCM_FRACTION_DENOMINATOR (y
));
5988 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
5990 else if (SCM_REALP (x
))
5992 if (SCM_I_INUMP (y
))
5994 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
5995 if (scm_is_eq (y
, SCM_INUM0
))
5997 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
5999 else if (SCM_BIGP (y
))
6001 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
6002 scm_remember_upto_here_1 (y
);
6003 return scm_from_double (result
);
6005 else if (SCM_REALP (y
))
6006 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
6007 else if (SCM_COMPLEXP (y
))
6008 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
6009 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
6010 else if (SCM_FRACTIONP (y
))
6011 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
6013 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6015 else if (SCM_COMPLEXP (x
))
6017 if (SCM_I_INUMP (y
))
6019 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
6020 if (scm_is_eq (y
, SCM_INUM0
))
6022 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
6023 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
6025 else if (SCM_BIGP (y
))
6027 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6028 scm_remember_upto_here_1 (y
);
6029 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
6030 z
* SCM_COMPLEX_IMAG (x
));
6032 else if (SCM_REALP (y
))
6033 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
6034 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
6035 else if (SCM_COMPLEXP (y
))
6037 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
6038 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
6039 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
6040 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
6042 else if (SCM_FRACTIONP (y
))
6044 double yy
= scm_i_fraction2double (y
);
6045 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
6046 yy
* SCM_COMPLEX_IMAG (x
));
6049 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6051 else if (SCM_FRACTIONP (x
))
6053 if (SCM_I_INUMP (y
))
6054 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
6055 SCM_FRACTION_DENOMINATOR (x
));
6056 else if (SCM_BIGP (y
))
6057 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
6058 SCM_FRACTION_DENOMINATOR (x
));
6059 else if (SCM_REALP (y
))
6060 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
6061 else if (SCM_COMPLEXP (y
))
6063 double xx
= scm_i_fraction2double (x
);
6064 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
6065 xx
* SCM_COMPLEX_IMAG (y
));
6067 else if (SCM_FRACTIONP (y
))
6068 /* a/b * c/d = ac / bd */
6069 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
6070 SCM_FRACTION_NUMERATOR (y
)),
6071 scm_product (SCM_FRACTION_DENOMINATOR (x
),
6072 SCM_FRACTION_DENOMINATOR (y
)));
6074 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6077 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
6080 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
6081 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
6082 #define ALLOW_DIVIDE_BY_ZERO
6083 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
6086 /* The code below for complex division is adapted from the GNU
6087 libstdc++, which adapted it from f2c's libF77, and is subject to
6090 /****************************************************************
6091 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
6093 Permission to use, copy, modify, and distribute this software
6094 and its documentation for any purpose and without fee is hereby
6095 granted, provided that the above copyright notice appear in all
6096 copies and that both that the copyright notice and this
6097 permission notice and warranty disclaimer appear in supporting
6098 documentation, and that the names of AT&T Bell Laboratories or
6099 Bellcore or any of their entities not be used in advertising or
6100 publicity pertaining to distribution of the software without
6101 specific, written prior permission.
6103 AT&T and Bellcore disclaim all warranties with regard to this
6104 software, including all implied warranties of merchantability
6105 and fitness. In no event shall AT&T or Bellcore be liable for
6106 any special, indirect or consequential damages or any damages
6107 whatsoever resulting from loss of use, data or profits, whether
6108 in an action of contract, negligence or other tortious action,
6109 arising out of or in connection with the use or performance of
6111 ****************************************************************/
6113 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
6114 (SCM x
, SCM y
, SCM rest
),
6115 "Divide the first argument by the product of the remaining\n"
6116 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
6118 #define FUNC_NAME s_scm_i_divide
6120 while (!scm_is_null (rest
))
6121 { x
= scm_divide (x
, y
);
6123 rest
= scm_cdr (rest
);
6125 return scm_divide (x
, y
);
6129 #define s_divide s_scm_i_divide
6130 #define g_divide g_scm_i_divide
6133 do_divide (SCM x
, SCM y
, int inexact
)
6134 #define FUNC_NAME s_divide
6138 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
6141 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
6142 else if (SCM_I_INUMP (x
))
6144 scm_t_inum xx
= SCM_I_INUM (x
);
6145 if (xx
== 1 || xx
== -1)
6147 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6149 scm_num_overflow (s_divide
);
6154 return scm_from_double (1.0 / (double) xx
);
6155 else return scm_i_make_ratio (SCM_INUM1
, x
);
6158 else if (SCM_BIGP (x
))
6161 return scm_from_double (1.0 / scm_i_big2dbl (x
));
6162 else return scm_i_make_ratio (SCM_INUM1
, x
);
6164 else if (SCM_REALP (x
))
6166 double xx
= SCM_REAL_VALUE (x
);
6167 #ifndef ALLOW_DIVIDE_BY_ZERO
6169 scm_num_overflow (s_divide
);
6172 return scm_from_double (1.0 / xx
);
6174 else if (SCM_COMPLEXP (x
))
6176 double r
= SCM_COMPLEX_REAL (x
);
6177 double i
= SCM_COMPLEX_IMAG (x
);
6178 if (fabs(r
) <= fabs(i
))
6181 double d
= i
* (1.0 + t
* t
);
6182 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
6187 double d
= r
* (1.0 + t
* t
);
6188 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
6191 else if (SCM_FRACTIONP (x
))
6192 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
6193 SCM_FRACTION_NUMERATOR (x
));
6195 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
6198 if (SCM_LIKELY (SCM_I_INUMP (x
)))
6200 scm_t_inum xx
= SCM_I_INUM (x
);
6201 if (SCM_LIKELY (SCM_I_INUMP (y
)))
6203 scm_t_inum yy
= SCM_I_INUM (y
);
6206 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6207 scm_num_overflow (s_divide
);
6209 return scm_from_double ((double) xx
/ (double) yy
);
6212 else if (xx
% yy
!= 0)
6215 return scm_from_double ((double) xx
/ (double) yy
);
6216 else return scm_i_make_ratio (x
, y
);
6220 scm_t_inum z
= xx
/ yy
;
6221 if (SCM_FIXABLE (z
))
6222 return SCM_I_MAKINUM (z
);
6224 return scm_i_inum2big (z
);
6227 else if (SCM_BIGP (y
))
6230 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
6231 else return scm_i_make_ratio (x
, y
);
6233 else if (SCM_REALP (y
))
6235 double yy
= SCM_REAL_VALUE (y
);
6236 #ifndef ALLOW_DIVIDE_BY_ZERO
6238 scm_num_overflow (s_divide
);
6241 return scm_from_double ((double) xx
/ yy
);
6243 else if (SCM_COMPLEXP (y
))
6246 complex_div
: /* y _must_ be a complex number */
6248 double r
= SCM_COMPLEX_REAL (y
);
6249 double i
= SCM_COMPLEX_IMAG (y
);
6250 if (fabs(r
) <= fabs(i
))
6253 double d
= i
* (1.0 + t
* t
);
6254 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
6259 double d
= r
* (1.0 + t
* t
);
6260 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
6264 else if (SCM_FRACTIONP (y
))
6265 /* a / b/c = ac / b */
6266 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
6267 SCM_FRACTION_NUMERATOR (y
));
6269 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6271 else if (SCM_BIGP (x
))
6273 if (SCM_I_INUMP (y
))
6275 scm_t_inum yy
= SCM_I_INUM (y
);
6278 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6279 scm_num_overflow (s_divide
);
6281 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
6282 scm_remember_upto_here_1 (x
);
6283 return (sgn
== 0) ? scm_nan () : scm_inf ();
6290 /* FIXME: HMM, what are the relative performance issues here?
6291 We need to test. Is it faster on average to test
6292 divisible_p, then perform whichever operation, or is it
6293 faster to perform the integer div opportunistically and
6294 switch to real if there's a remainder? For now we take the
6295 middle ground: test, then if divisible, use the faster div
6298 scm_t_inum abs_yy
= yy
< 0 ? -yy
: yy
;
6299 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
6303 SCM result
= scm_i_mkbig ();
6304 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
6305 scm_remember_upto_here_1 (x
);
6307 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
6308 return scm_i_normbig (result
);
6313 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
6314 else return scm_i_make_ratio (x
, y
);
6318 else if (SCM_BIGP (y
))
6323 /* It's easily possible for the ratio x/y to fit a double
6324 but one or both x and y be too big to fit a double,
6325 hence the use of mpq_get_d rather than converting and
6328 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
6329 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
6330 return scm_from_double (mpq_get_d (q
));
6334 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
6338 SCM result
= scm_i_mkbig ();
6339 mpz_divexact (SCM_I_BIG_MPZ (result
),
6342 scm_remember_upto_here_2 (x
, y
);
6343 return scm_i_normbig (result
);
6346 return scm_i_make_ratio (x
, y
);
6349 else if (SCM_REALP (y
))
6351 double yy
= SCM_REAL_VALUE (y
);
6352 #ifndef ALLOW_DIVIDE_BY_ZERO
6354 scm_num_overflow (s_divide
);
6357 return scm_from_double (scm_i_big2dbl (x
) / yy
);
6359 else if (SCM_COMPLEXP (y
))
6361 a
= scm_i_big2dbl (x
);
6364 else if (SCM_FRACTIONP (y
))
6365 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
6366 SCM_FRACTION_NUMERATOR (y
));
6368 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6370 else if (SCM_REALP (x
))
6372 double rx
= SCM_REAL_VALUE (x
);
6373 if (SCM_I_INUMP (y
))
6375 scm_t_inum yy
= SCM_I_INUM (y
);
6376 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6378 scm_num_overflow (s_divide
);
6381 return scm_from_double (rx
/ (double) yy
);
6383 else if (SCM_BIGP (y
))
6385 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6386 scm_remember_upto_here_1 (y
);
6387 return scm_from_double (rx
/ dby
);
6389 else if (SCM_REALP (y
))
6391 double yy
= SCM_REAL_VALUE (y
);
6392 #ifndef ALLOW_DIVIDE_BY_ZERO
6394 scm_num_overflow (s_divide
);
6397 return scm_from_double (rx
/ yy
);
6399 else if (SCM_COMPLEXP (y
))
6404 else if (SCM_FRACTIONP (y
))
6405 return scm_from_double (rx
/ scm_i_fraction2double (y
));
6407 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6409 else if (SCM_COMPLEXP (x
))
6411 double rx
= SCM_COMPLEX_REAL (x
);
6412 double ix
= SCM_COMPLEX_IMAG (x
);
6413 if (SCM_I_INUMP (y
))
6415 scm_t_inum yy
= SCM_I_INUM (y
);
6416 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6418 scm_num_overflow (s_divide
);
6423 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
6426 else if (SCM_BIGP (y
))
6428 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6429 scm_remember_upto_here_1 (y
);
6430 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
6432 else if (SCM_REALP (y
))
6434 double yy
= SCM_REAL_VALUE (y
);
6435 #ifndef ALLOW_DIVIDE_BY_ZERO
6437 scm_num_overflow (s_divide
);
6440 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
6442 else if (SCM_COMPLEXP (y
))
6444 double ry
= SCM_COMPLEX_REAL (y
);
6445 double iy
= SCM_COMPLEX_IMAG (y
);
6446 if (fabs(ry
) <= fabs(iy
))
6449 double d
= iy
* (1.0 + t
* t
);
6450 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
6455 double d
= ry
* (1.0 + t
* t
);
6456 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
6459 else if (SCM_FRACTIONP (y
))
6461 double yy
= scm_i_fraction2double (y
);
6462 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
6465 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6467 else if (SCM_FRACTIONP (x
))
6469 if (SCM_I_INUMP (y
))
6471 scm_t_inum yy
= SCM_I_INUM (y
);
6472 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6474 scm_num_overflow (s_divide
);
6477 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
6478 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
6480 else if (SCM_BIGP (y
))
6482 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
6483 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
6485 else if (SCM_REALP (y
))
6487 double yy
= SCM_REAL_VALUE (y
);
6488 #ifndef ALLOW_DIVIDE_BY_ZERO
6490 scm_num_overflow (s_divide
);
6493 return scm_from_double (scm_i_fraction2double (x
) / yy
);
6495 else if (SCM_COMPLEXP (y
))
6497 a
= scm_i_fraction2double (x
);
6500 else if (SCM_FRACTIONP (y
))
6501 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
6502 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
6504 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6507 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
6511 scm_divide (SCM x
, SCM y
)
6513 return do_divide (x
, y
, 0);
6516 static SCM
scm_divide2real (SCM x
, SCM y
)
6518 return do_divide (x
, y
, 1);
6524 scm_c_truncate (double x
)
6535 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
6536 half-way case (ie. when x is an integer plus 0.5) going upwards.
6537 Then half-way cases are identified and adjusted down if the
6538 round-upwards didn't give the desired even integer.
6540 "plus_half == result" identifies a half-way case. If plus_half, which is
6541 x + 0.5, is an integer then x must be an integer plus 0.5.
6543 An odd "result" value is identified with result/2 != floor(result/2).
6544 This is done with plus_half, since that value is ready for use sooner in
6545 a pipelined cpu, and we're already requiring plus_half == result.
6547 Note however that we need to be careful when x is big and already an
6548 integer. In that case "x+0.5" may round to an adjacent integer, causing
6549 us to return such a value, incorrectly. For instance if the hardware is
6550 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
6551 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
6552 returned. Or if the hardware is in round-upwards mode, then other bigger
6553 values like say x == 2^128 will see x+0.5 rounding up to the next higher
6554 representable value, 2^128+2^76 (or whatever), again incorrect.
6556 These bad roundings of x+0.5 are avoided by testing at the start whether
6557 x is already an integer. If it is then clearly that's the desired result
6558 already. And if it's not then the exponent must be small enough to allow
6559 an 0.5 to be represented, and hence added without a bad rounding. */
6562 scm_c_round (double x
)
6564 double plus_half
, result
;
6569 plus_half
= x
+ 0.5;
6570 result
= floor (plus_half
);
6571 /* Adjust so that the rounding is towards even. */
6572 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
6577 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
6579 "Round the number @var{x} towards zero.")
6580 #define FUNC_NAME s_scm_truncate_number
6582 if (scm_is_false (scm_negative_p (x
)))
6583 return scm_floor (x
);
6585 return scm_ceiling (x
);
6589 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
6591 "Round the number @var{x} towards the nearest integer. "
6592 "When it is exactly halfway between two integers, "
6593 "round towards the even one.")
6594 #define FUNC_NAME s_scm_round_number
6596 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6598 else if (SCM_REALP (x
))
6599 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
6602 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
6603 single quotient+remainder division then examining to see which way
6604 the rounding should go. */
6605 SCM plus_half
= scm_sum (x
, exactly_one_half
);
6606 SCM result
= scm_floor (plus_half
);
6607 /* Adjust so that the rounding is towards even. */
6608 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
6609 && scm_is_true (scm_odd_p (result
)))
6610 return scm_difference (result
, SCM_INUM1
);
6617 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
6619 "Round the number @var{x} towards minus infinity.")
6620 #define FUNC_NAME s_scm_floor
6622 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6624 else if (SCM_REALP (x
))
6625 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
6626 else if (SCM_FRACTIONP (x
))
6628 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
6629 SCM_FRACTION_DENOMINATOR (x
));
6630 if (scm_is_false (scm_negative_p (x
)))
6632 /* For positive x, rounding towards zero is correct. */
6637 /* For negative x, we need to return q-1 unless x is an
6638 integer. But fractions are never integer, per our
6640 return scm_difference (q
, SCM_INUM1
);
6644 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
6648 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
6650 "Round the number @var{x} towards infinity.")
6651 #define FUNC_NAME s_scm_ceiling
6653 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6655 else if (SCM_REALP (x
))
6656 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
6657 else if (SCM_FRACTIONP (x
))
6659 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
6660 SCM_FRACTION_DENOMINATOR (x
));
6661 if (scm_is_false (scm_positive_p (x
)))
6663 /* For negative x, rounding towards zero is correct. */
6668 /* For positive x, we need to return q+1 unless x is an
6669 integer. But fractions are never integer, per our
6671 return scm_sum (q
, SCM_INUM1
);
6675 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
6679 /* sin/cos/tan/asin/acos/atan
6680 sinh/cosh/tanh/asinh/acosh/atanh
6681 Derived from "Transcen.scm", Complex trancendental functions for SCM.
6682 Written by Jerry D. Hedden, (C) FSF.
6683 See the file `COPYING' for terms applying to this program. */
6685 SCM_PRIMITIVE_GENERIC (scm_expt
, "expt", 2, 0, 0,
6687 "Return @var{x} raised to the power of @var{y}.")
6688 #define FUNC_NAME s_scm_expt
6690 if (scm_is_integer (y
))
6692 if (scm_is_true (scm_exact_p (y
)))
6693 return scm_integer_expt (x
, y
);
6696 /* Here we handle the case where the exponent is an inexact
6697 integer. We make the exponent exact in order to use
6698 scm_integer_expt, and thus avoid the spurious imaginary
6699 parts that may result from round-off errors in the general
6700 e^(y log x) method below (for example when squaring a large
6701 negative number). In this case, we must return an inexact
6702 result for correctness. We also make the base inexact so
6703 that scm_integer_expt will use fast inexact arithmetic
6704 internally. Note that making the base inexact is not
6705 sufficient to guarantee an inexact result, because
6706 scm_integer_expt will return an exact 1 when the exponent
6707 is 0, even if the base is inexact. */
6708 return scm_exact_to_inexact
6709 (scm_integer_expt (scm_exact_to_inexact (x
),
6710 scm_inexact_to_exact (y
)));
6713 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
6715 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
6717 else if (scm_is_complex (x
) && scm_is_complex (y
))
6718 return scm_exp (scm_product (scm_log (x
), y
));
6719 else if (scm_is_complex (x
))
6720 SCM_WTA_DISPATCH_2 (g_scm_expt
, x
, y
, SCM_ARG2
, s_scm_expt
);
6722 SCM_WTA_DISPATCH_2 (g_scm_expt
, x
, y
, SCM_ARG1
, s_scm_expt
);
6726 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
6728 "Compute the sine of @var{z}.")
6729 #define FUNC_NAME s_scm_sin
6731 if (scm_is_real (z
))
6732 return scm_from_double (sin (scm_to_double (z
)));
6733 else if (SCM_COMPLEXP (z
))
6735 x
= SCM_COMPLEX_REAL (z
);
6736 y
= SCM_COMPLEX_IMAG (z
);
6737 return scm_c_make_rectangular (sin (x
) * cosh (y
),
6738 cos (x
) * sinh (y
));
6741 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
6745 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
6747 "Compute the cosine of @var{z}.")
6748 #define FUNC_NAME s_scm_cos
6750 if (scm_is_real (z
))
6751 return scm_from_double (cos (scm_to_double (z
)));
6752 else if (SCM_COMPLEXP (z
))
6754 x
= SCM_COMPLEX_REAL (z
);
6755 y
= SCM_COMPLEX_IMAG (z
);
6756 return scm_c_make_rectangular (cos (x
) * cosh (y
),
6757 -sin (x
) * sinh (y
));
6760 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
6764 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
6766 "Compute the tangent of @var{z}.")
6767 #define FUNC_NAME s_scm_tan
6769 if (scm_is_real (z
))
6770 return scm_from_double (tan (scm_to_double (z
)));
6771 else if (SCM_COMPLEXP (z
))
6773 x
= 2.0 * SCM_COMPLEX_REAL (z
);
6774 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
6775 w
= cos (x
) + cosh (y
);
6776 #ifndef ALLOW_DIVIDE_BY_ZERO
6778 scm_num_overflow (s_scm_tan
);
6780 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
6783 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
6787 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
6789 "Compute the hyperbolic sine of @var{z}.")
6790 #define FUNC_NAME s_scm_sinh
6792 if (scm_is_real (z
))
6793 return scm_from_double (sinh (scm_to_double (z
)));
6794 else if (SCM_COMPLEXP (z
))
6796 x
= SCM_COMPLEX_REAL (z
);
6797 y
= SCM_COMPLEX_IMAG (z
);
6798 return scm_c_make_rectangular (sinh (x
) * cos (y
),
6799 cosh (x
) * sin (y
));
6802 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
6806 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
6808 "Compute the hyperbolic cosine of @var{z}.")
6809 #define FUNC_NAME s_scm_cosh
6811 if (scm_is_real (z
))
6812 return scm_from_double (cosh (scm_to_double (z
)));
6813 else if (SCM_COMPLEXP (z
))
6815 x
= SCM_COMPLEX_REAL (z
);
6816 y
= SCM_COMPLEX_IMAG (z
);
6817 return scm_c_make_rectangular (cosh (x
) * cos (y
),
6818 sinh (x
) * sin (y
));
6821 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
6825 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
6827 "Compute the hyperbolic tangent of @var{z}.")
6828 #define FUNC_NAME s_scm_tanh
6830 if (scm_is_real (z
))
6831 return scm_from_double (tanh (scm_to_double (z
)));
6832 else if (SCM_COMPLEXP (z
))
6834 x
= 2.0 * SCM_COMPLEX_REAL (z
);
6835 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
6836 w
= cosh (x
) + cos (y
);
6837 #ifndef ALLOW_DIVIDE_BY_ZERO
6839 scm_num_overflow (s_scm_tanh
);
6841 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
6844 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
6848 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
6850 "Compute the arc sine of @var{z}.")
6851 #define FUNC_NAME s_scm_asin
6853 if (scm_is_real (z
))
6855 double w
= scm_to_double (z
);
6856 if (w
>= -1.0 && w
<= 1.0)
6857 return scm_from_double (asin (w
));
6859 return scm_product (scm_c_make_rectangular (0, -1),
6860 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
6862 else if (SCM_COMPLEXP (z
))
6864 x
= SCM_COMPLEX_REAL (z
);
6865 y
= SCM_COMPLEX_IMAG (z
);
6866 return scm_product (scm_c_make_rectangular (0, -1),
6867 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
6870 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
6874 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
6876 "Compute the arc cosine of @var{z}.")
6877 #define FUNC_NAME s_scm_acos
6879 if (scm_is_real (z
))
6881 double w
= scm_to_double (z
);
6882 if (w
>= -1.0 && w
<= 1.0)
6883 return scm_from_double (acos (w
));
6885 return scm_sum (scm_from_double (acos (0.0)),
6886 scm_product (scm_c_make_rectangular (0, 1),
6887 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
6889 else if (SCM_COMPLEXP (z
))
6891 x
= SCM_COMPLEX_REAL (z
);
6892 y
= SCM_COMPLEX_IMAG (z
);
6893 return scm_sum (scm_from_double (acos (0.0)),
6894 scm_product (scm_c_make_rectangular (0, 1),
6895 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
6898 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
6902 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
6904 "With one argument, compute the arc tangent of @var{z}.\n"
6905 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
6906 "using the sign of @var{z} and @var{y} to determine the quadrant.")
6907 #define FUNC_NAME s_scm_atan
6911 if (scm_is_real (z
))
6912 return scm_from_double (atan (scm_to_double (z
)));
6913 else if (SCM_COMPLEXP (z
))
6916 v
= SCM_COMPLEX_REAL (z
);
6917 w
= SCM_COMPLEX_IMAG (z
);
6918 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
6919 scm_c_make_rectangular (v
, w
+ 1.0))),
6920 scm_c_make_rectangular (0, 2));
6923 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
6925 else if (scm_is_real (z
))
6927 if (scm_is_real (y
))
6928 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
6930 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
6933 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
6937 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
6939 "Compute the inverse hyperbolic sine of @var{z}.")
6940 #define FUNC_NAME s_scm_sys_asinh
6942 if (scm_is_real (z
))
6943 return scm_from_double (asinh (scm_to_double (z
)));
6944 else if (scm_is_number (z
))
6945 return scm_log (scm_sum (z
,
6946 scm_sqrt (scm_sum (scm_product (z
, z
),
6949 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
6953 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
6955 "Compute the inverse hyperbolic cosine of @var{z}.")
6956 #define FUNC_NAME s_scm_sys_acosh
6958 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
6959 return scm_from_double (acosh (scm_to_double (z
)));
6960 else if (scm_is_number (z
))
6961 return scm_log (scm_sum (z
,
6962 scm_sqrt (scm_difference (scm_product (z
, z
),
6965 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
6969 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
6971 "Compute the inverse hyperbolic tangent of @var{z}.")
6972 #define FUNC_NAME s_scm_sys_atanh
6974 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
6975 return scm_from_double (atanh (scm_to_double (z
)));
6976 else if (scm_is_number (z
))
6977 return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1
, z
),
6978 scm_difference (SCM_INUM1
, z
))),
6981 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
6986 scm_c_make_rectangular (double re
, double im
)
6989 return scm_from_double (re
);
6994 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
6996 SCM_SET_CELL_TYPE (z
, scm_tc16_complex
);
6997 SCM_COMPLEX_REAL (z
) = re
;
6998 SCM_COMPLEX_IMAG (z
) = im
;
7003 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
7004 (SCM real_part
, SCM imaginary_part
),
7005 "Return a complex number constructed of the given @var{real-part} "
7006 "and @var{imaginary-part} parts.")
7007 #define FUNC_NAME s_scm_make_rectangular
7009 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
7010 SCM_ARG1
, FUNC_NAME
, "real");
7011 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
7012 SCM_ARG2
, FUNC_NAME
, "real");
7013 return scm_c_make_rectangular (scm_to_double (real_part
),
7014 scm_to_double (imaginary_part
));
7019 scm_c_make_polar (double mag
, double ang
)
7023 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
7024 use it on Glibc-based systems that have it (it's a GNU extension). See
7025 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
7027 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
7028 sincos (ang
, &s
, &c
);
7033 return scm_c_make_rectangular (mag
* c
, mag
* s
);
7036 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
7038 "Return the complex number @var{x} * e^(i * @var{y}).")
7039 #define FUNC_NAME s_scm_make_polar
7041 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
7042 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
7043 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
7048 SCM_PRIMITIVE_GENERIC (scm_real_part
, "real-part", 1, 0, 0,
7050 "Return the real part of the number @var{z}.")
7051 #define FUNC_NAME s_scm_real_part
7053 if (SCM_COMPLEXP (z
))
7054 return scm_from_double (SCM_COMPLEX_REAL (z
));
7055 else if (SCM_I_INUMP (z
) || SCM_BIGP (z
) || SCM_REALP (z
) || SCM_FRACTIONP (z
))
7058 SCM_WTA_DISPATCH_1 (g_scm_real_part
, z
, SCM_ARG1
, s_scm_real_part
);
7063 SCM_PRIMITIVE_GENERIC (scm_imag_part
, "imag-part", 1, 0, 0,
7065 "Return the imaginary part of the number @var{z}.")
7066 #define FUNC_NAME s_scm_imag_part
7068 if (SCM_COMPLEXP (z
))
7069 return scm_from_double (SCM_COMPLEX_IMAG (z
));
7070 else if (SCM_REALP (z
))
7072 else if (SCM_I_INUMP (z
) || SCM_BIGP (z
) || SCM_FRACTIONP (z
))
7075 SCM_WTA_DISPATCH_1 (g_scm_imag_part
, z
, SCM_ARG1
, s_scm_imag_part
);
7079 SCM_PRIMITIVE_GENERIC (scm_numerator
, "numerator", 1, 0, 0,
7081 "Return the numerator of the number @var{z}.")
7082 #define FUNC_NAME s_scm_numerator
7084 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7086 else if (SCM_FRACTIONP (z
))
7087 return SCM_FRACTION_NUMERATOR (z
);
7088 else if (SCM_REALP (z
))
7089 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
7091 SCM_WTA_DISPATCH_1 (g_scm_numerator
, z
, SCM_ARG1
, s_scm_numerator
);
7096 SCM_PRIMITIVE_GENERIC (scm_denominator
, "denominator", 1, 0, 0,
7098 "Return the denominator of the number @var{z}.")
7099 #define FUNC_NAME s_scm_denominator
7101 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7103 else if (SCM_FRACTIONP (z
))
7104 return SCM_FRACTION_DENOMINATOR (z
);
7105 else if (SCM_REALP (z
))
7106 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
7108 SCM_WTA_DISPATCH_1 (g_scm_denominator
, z
, SCM_ARG1
, s_scm_denominator
);
7113 SCM_PRIMITIVE_GENERIC (scm_magnitude
, "magnitude", 1, 0, 0,
7115 "Return the magnitude of the number @var{z}. This is the same as\n"
7116 "@code{abs} for real arguments, but also allows complex numbers.")
7117 #define FUNC_NAME s_scm_magnitude
7119 if (SCM_I_INUMP (z
))
7121 scm_t_inum zz
= SCM_I_INUM (z
);
7124 else if (SCM_POSFIXABLE (-zz
))
7125 return SCM_I_MAKINUM (-zz
);
7127 return scm_i_inum2big (-zz
);
7129 else if (SCM_BIGP (z
))
7131 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
7132 scm_remember_upto_here_1 (z
);
7134 return scm_i_clonebig (z
, 0);
7138 else if (SCM_REALP (z
))
7139 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
7140 else if (SCM_COMPLEXP (z
))
7141 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
7142 else if (SCM_FRACTIONP (z
))
7144 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
7146 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
7147 SCM_FRACTION_DENOMINATOR (z
));
7150 SCM_WTA_DISPATCH_1 (g_scm_magnitude
, z
, SCM_ARG1
, s_scm_magnitude
);
7155 SCM_PRIMITIVE_GENERIC (scm_angle
, "angle", 1, 0, 0,
7157 "Return the angle of the complex number @var{z}.")
7158 #define FUNC_NAME s_scm_angle
7160 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
7161 flo0 to save allocating a new flonum with scm_from_double each time.
7162 But if atan2 follows the floating point rounding mode, then the value
7163 is not a constant. Maybe it'd be close enough though. */
7164 if (SCM_I_INUMP (z
))
7166 if (SCM_I_INUM (z
) >= 0)
7169 return scm_from_double (atan2 (0.0, -1.0));
7171 else if (SCM_BIGP (z
))
7173 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
7174 scm_remember_upto_here_1 (z
);
7176 return scm_from_double (atan2 (0.0, -1.0));
7180 else if (SCM_REALP (z
))
7182 if (SCM_REAL_VALUE (z
) >= 0)
7185 return scm_from_double (atan2 (0.0, -1.0));
7187 else if (SCM_COMPLEXP (z
))
7188 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
7189 else if (SCM_FRACTIONP (z
))
7191 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
7193 else return scm_from_double (atan2 (0.0, -1.0));
7196 SCM_WTA_DISPATCH_1 (g_scm_angle
, z
, SCM_ARG1
, s_scm_angle
);
7201 SCM_PRIMITIVE_GENERIC (scm_exact_to_inexact
, "exact->inexact", 1, 0, 0,
7203 "Convert the number @var{z} to its inexact representation.\n")
7204 #define FUNC_NAME s_scm_exact_to_inexact
7206 if (SCM_I_INUMP (z
))
7207 return scm_from_double ((double) SCM_I_INUM (z
));
7208 else if (SCM_BIGP (z
))
7209 return scm_from_double (scm_i_big2dbl (z
));
7210 else if (SCM_FRACTIONP (z
))
7211 return scm_from_double (scm_i_fraction2double (z
));
7212 else if (SCM_INEXACTP (z
))
7215 SCM_WTA_DISPATCH_1 (g_scm_exact_to_inexact
, z
, 1, s_scm_exact_to_inexact
);
7220 SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
7222 "Return an exact number that is numerically closest to @var{z}.")
7223 #define FUNC_NAME s_scm_inexact_to_exact
7225 if (SCM_I_INUMP (z
) || SCM_BIGP (z
))
7227 else if (SCM_REALP (z
))
7229 if (!DOUBLE_IS_FINITE (SCM_REAL_VALUE (z
)))
7230 SCM_OUT_OF_RANGE (1, z
);
7237 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
7238 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
7239 scm_i_mpz2num (mpq_denref (frac
)));
7241 /* When scm_i_make_ratio throws, we leak the memory allocated
7248 else if (SCM_FRACTIONP (z
))
7251 SCM_WTA_DISPATCH_1 (g_scm_inexact_to_exact
, z
, 1, s_scm_inexact_to_exact
);
7255 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
7257 "Returns the @emph{simplest} rational number differing\n"
7258 "from @var{x} by no more than @var{eps}.\n"
7260 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
7261 "exact result when both its arguments are exact. Thus, you might need\n"
7262 "to use @code{inexact->exact} on the arguments.\n"
7265 "(rationalize (inexact->exact 1.2) 1/100)\n"
7268 #define FUNC_NAME s_scm_rationalize
7270 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
7271 SCM_ASSERT_TYPE (scm_is_real (eps
), eps
, SCM_ARG2
, FUNC_NAME
, "real");
7272 eps
= scm_abs (eps
);
7273 if (scm_is_false (scm_positive_p (eps
)))
7275 /* eps is either zero or a NaN */
7276 if (scm_is_true (scm_nan_p (eps
)))
7278 else if (SCM_INEXACTP (eps
))
7279 return scm_exact_to_inexact (x
);
7283 else if (scm_is_false (scm_finite_p (eps
)))
7285 if (scm_is_true (scm_finite_p (x
)))
7290 else if (scm_is_false (scm_finite_p (x
))) /* checks for both inf and nan */
7292 else if (scm_is_false (scm_less_p (scm_floor (scm_sum (x
, eps
)),
7293 scm_ceiling (scm_difference (x
, eps
)))))
7295 /* There's an integer within range; we want the one closest to zero */
7296 if (scm_is_false (scm_less_p (eps
, scm_abs (x
))))
7298 /* zero is within range */
7299 if (SCM_INEXACTP (x
) || SCM_INEXACTP (eps
))
7304 else if (scm_is_true (scm_positive_p (x
)))
7305 return scm_ceiling (scm_difference (x
, eps
));
7307 return scm_floor (scm_sum (x
, eps
));
7311 /* Use continued fractions to find closest ratio. All
7312 arithmetic is done with exact numbers.
7315 SCM ex
= scm_inexact_to_exact (x
);
7316 SCM int_part
= scm_floor (ex
);
7318 SCM a1
= SCM_INUM0
, a2
= SCM_INUM1
, a
= SCM_INUM0
;
7319 SCM b1
= SCM_INUM1
, b2
= SCM_INUM0
, b
= SCM_INUM0
;
7323 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
7324 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
7326 /* We stop after a million iterations just to be absolutely sure
7327 that we don't go into an infinite loop. The process normally
7328 converges after less than a dozen iterations.
7331 while (++i
< 1000000)
7333 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
7334 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
7335 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
7337 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
7338 eps
))) /* abs(x-a/b) <= eps */
7340 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
7341 if (SCM_INEXACTP (x
) || SCM_INEXACTP (eps
))
7342 return scm_exact_to_inexact (res
);
7346 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
7348 tt
= scm_floor (rx
); /* tt = floor (rx) */
7354 scm_num_overflow (s_scm_rationalize
);
7359 /* conversion functions */
7362 scm_is_integer (SCM val
)
7364 return scm_is_true (scm_integer_p (val
));
7368 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
7370 if (SCM_I_INUMP (val
))
7372 scm_t_signed_bits n
= SCM_I_INUM (val
);
7373 return n
>= min
&& n
<= max
;
7375 else if (SCM_BIGP (val
))
7377 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
7379 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
7381 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
7383 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
7384 return n
>= min
&& n
<= max
;
7394 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
7395 > CHAR_BIT
*sizeof (scm_t_uintmax
))
7398 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
7399 SCM_I_BIG_MPZ (val
));
7401 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
7413 return n
>= min
&& n
<= max
;
7421 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
7423 if (SCM_I_INUMP (val
))
7425 scm_t_signed_bits n
= SCM_I_INUM (val
);
7426 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
7428 else if (SCM_BIGP (val
))
7430 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
7432 else if (max
<= ULONG_MAX
)
7434 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
7436 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
7437 return n
>= min
&& n
<= max
;
7447 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
7450 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
7451 > CHAR_BIT
*sizeof (scm_t_uintmax
))
7454 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
7455 SCM_I_BIG_MPZ (val
));
7457 return n
>= min
&& n
<= max
;
7465 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
7467 scm_error (scm_out_of_range_key
,
7469 "Value out of range ~S to ~S: ~S",
7470 scm_list_3 (min
, max
, bad_val
),
7471 scm_list_1 (bad_val
));
7474 #define TYPE scm_t_intmax
7475 #define TYPE_MIN min
7476 #define TYPE_MAX max
7477 #define SIZEOF_TYPE 0
7478 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
7479 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
7480 #include "libguile/conv-integer.i.c"
7482 #define TYPE scm_t_uintmax
7483 #define TYPE_MIN min
7484 #define TYPE_MAX max
7485 #define SIZEOF_TYPE 0
7486 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
7487 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
7488 #include "libguile/conv-uinteger.i.c"
7490 #define TYPE scm_t_int8
7491 #define TYPE_MIN SCM_T_INT8_MIN
7492 #define TYPE_MAX SCM_T_INT8_MAX
7493 #define SIZEOF_TYPE 1
7494 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
7495 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
7496 #include "libguile/conv-integer.i.c"
7498 #define TYPE scm_t_uint8
7500 #define TYPE_MAX SCM_T_UINT8_MAX
7501 #define SIZEOF_TYPE 1
7502 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
7503 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
7504 #include "libguile/conv-uinteger.i.c"
7506 #define TYPE scm_t_int16
7507 #define TYPE_MIN SCM_T_INT16_MIN
7508 #define TYPE_MAX SCM_T_INT16_MAX
7509 #define SIZEOF_TYPE 2
7510 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
7511 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
7512 #include "libguile/conv-integer.i.c"
7514 #define TYPE scm_t_uint16
7516 #define TYPE_MAX SCM_T_UINT16_MAX
7517 #define SIZEOF_TYPE 2
7518 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
7519 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
7520 #include "libguile/conv-uinteger.i.c"
7522 #define TYPE scm_t_int32
7523 #define TYPE_MIN SCM_T_INT32_MIN
7524 #define TYPE_MAX SCM_T_INT32_MAX
7525 #define SIZEOF_TYPE 4
7526 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
7527 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
7528 #include "libguile/conv-integer.i.c"
7530 #define TYPE scm_t_uint32
7532 #define TYPE_MAX SCM_T_UINT32_MAX
7533 #define SIZEOF_TYPE 4
7534 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
7535 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
7536 #include "libguile/conv-uinteger.i.c"
7538 #define TYPE scm_t_wchar
7539 #define TYPE_MIN (scm_t_int32)-1
7540 #define TYPE_MAX (scm_t_int32)0x10ffff
7541 #define SIZEOF_TYPE 4
7542 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
7543 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
7544 #include "libguile/conv-integer.i.c"
7546 #define TYPE scm_t_int64
7547 #define TYPE_MIN SCM_T_INT64_MIN
7548 #define TYPE_MAX SCM_T_INT64_MAX
7549 #define SIZEOF_TYPE 8
7550 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
7551 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
7552 #include "libguile/conv-integer.i.c"
7554 #define TYPE scm_t_uint64
7556 #define TYPE_MAX SCM_T_UINT64_MAX
7557 #define SIZEOF_TYPE 8
7558 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
7559 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
7560 #include "libguile/conv-uinteger.i.c"
7563 scm_to_mpz (SCM val
, mpz_t rop
)
7565 if (SCM_I_INUMP (val
))
7566 mpz_set_si (rop
, SCM_I_INUM (val
));
7567 else if (SCM_BIGP (val
))
7568 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
7570 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
7574 scm_from_mpz (mpz_t val
)
7576 return scm_i_mpz2num (val
);
7580 scm_is_real (SCM val
)
7582 return scm_is_true (scm_real_p (val
));
7586 scm_is_rational (SCM val
)
7588 return scm_is_true (scm_rational_p (val
));
7592 scm_to_double (SCM val
)
7594 if (SCM_I_INUMP (val
))
7595 return SCM_I_INUM (val
);
7596 else if (SCM_BIGP (val
))
7597 return scm_i_big2dbl (val
);
7598 else if (SCM_FRACTIONP (val
))
7599 return scm_i_fraction2double (val
);
7600 else if (SCM_REALP (val
))
7601 return SCM_REAL_VALUE (val
);
7603 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
7607 scm_from_double (double val
)
7611 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
7613 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
7614 SCM_REAL_VALUE (z
) = val
;
7619 #if SCM_ENABLE_DEPRECATED == 1
7622 scm_num2float (SCM num
, unsigned long pos
, const char *s_caller
)
7624 scm_c_issue_deprecation_warning
7625 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
7629 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
7633 scm_out_of_range (NULL
, num
);
7636 return scm_to_double (num
);
7640 scm_num2double (SCM num
, unsigned long pos
, const char *s_caller
)
7642 scm_c_issue_deprecation_warning
7643 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
7647 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
7651 scm_out_of_range (NULL
, num
);
7654 return scm_to_double (num
);
7660 scm_is_complex (SCM val
)
7662 return scm_is_true (scm_complex_p (val
));
7666 scm_c_real_part (SCM z
)
7668 if (SCM_COMPLEXP (z
))
7669 return SCM_COMPLEX_REAL (z
);
7672 /* Use the scm_real_part to get proper error checking and
7675 return scm_to_double (scm_real_part (z
));
7680 scm_c_imag_part (SCM z
)
7682 if (SCM_COMPLEXP (z
))
7683 return SCM_COMPLEX_IMAG (z
);
7686 /* Use the scm_imag_part to get proper error checking and
7687 dispatching. The result will almost always be 0.0, but not
7690 return scm_to_double (scm_imag_part (z
));
7695 scm_c_magnitude (SCM z
)
7697 return scm_to_double (scm_magnitude (z
));
7703 return scm_to_double (scm_angle (z
));
7707 scm_is_number (SCM z
)
7709 return scm_is_true (scm_number_p (z
));
7713 /* In the following functions we dispatch to the real-arg funcs like log()
7714 when we know the arg is real, instead of just handing everything to
7715 clog() for instance. This is in case clog() doesn't optimize for a
7716 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
7717 well use it to go straight to the applicable C func. */
7719 SCM_PRIMITIVE_GENERIC (scm_log
, "log", 1, 0, 0,
7721 "Return the natural logarithm of @var{z}.")
7722 #define FUNC_NAME s_scm_log
7724 if (SCM_COMPLEXP (z
))
7726 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
7727 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
7729 double re
= SCM_COMPLEX_REAL (z
);
7730 double im
= SCM_COMPLEX_IMAG (z
);
7731 return scm_c_make_rectangular (log (hypot (re
, im
)),
7735 else if (SCM_NUMBERP (z
))
7737 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
7738 although the value itself overflows. */
7739 double re
= scm_to_double (z
);
7740 double l
= log (fabs (re
));
7742 return scm_from_double (l
);
7744 return scm_c_make_rectangular (l
, M_PI
);
7747 SCM_WTA_DISPATCH_1 (g_scm_log
, z
, 1, s_scm_log
);
7752 SCM_PRIMITIVE_GENERIC (scm_log10
, "log10", 1, 0, 0,
7754 "Return the base 10 logarithm of @var{z}.")
7755 #define FUNC_NAME s_scm_log10
7757 if (SCM_COMPLEXP (z
))
7759 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
7760 clog() and a multiply by M_LOG10E, rather than the fallback
7761 log10+hypot+atan2.) */
7762 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
7763 && defined SCM_COMPLEX_VALUE
7764 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
7766 double re
= SCM_COMPLEX_REAL (z
);
7767 double im
= SCM_COMPLEX_IMAG (z
);
7768 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
7769 M_LOG10E
* atan2 (im
, re
));
7772 else if (SCM_NUMBERP (z
))
7774 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
7775 although the value itself overflows. */
7776 double re
= scm_to_double (z
);
7777 double l
= log10 (fabs (re
));
7779 return scm_from_double (l
);
7781 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
7784 SCM_WTA_DISPATCH_1 (g_scm_log10
, z
, 1, s_scm_log10
);
7789 SCM_PRIMITIVE_GENERIC (scm_exp
, "exp", 1, 0, 0,
7791 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
7792 "base of natural logarithms (2.71828@dots{}).")
7793 #define FUNC_NAME s_scm_exp
7795 if (SCM_COMPLEXP (z
))
7797 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
7798 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
7800 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
7801 SCM_COMPLEX_IMAG (z
));
7804 else if (SCM_NUMBERP (z
))
7806 /* When z is a negative bignum the conversion to double overflows,
7807 giving -infinity, but that's ok, the exp is still 0.0. */
7808 return scm_from_double (exp (scm_to_double (z
)));
7811 SCM_WTA_DISPATCH_1 (g_scm_exp
, z
, 1, s_scm_exp
);
7816 SCM_PRIMITIVE_GENERIC (scm_sqrt
, "sqrt", 1, 0, 0,
7818 "Return the square root of @var{z}. Of the two possible roots\n"
7819 "(positive and negative), the one with the a positive real part\n"
7820 "is returned, or if that's zero then a positive imaginary part.\n"
7824 "(sqrt 9.0) @result{} 3.0\n"
7825 "(sqrt -9.0) @result{} 0.0+3.0i\n"
7826 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
7827 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
7829 #define FUNC_NAME s_scm_sqrt
7831 if (SCM_COMPLEXP (z
))
7833 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
7834 && defined SCM_COMPLEX_VALUE
7835 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (z
)));
7837 double re
= SCM_COMPLEX_REAL (z
);
7838 double im
= SCM_COMPLEX_IMAG (z
);
7839 return scm_c_make_polar (sqrt (hypot (re
, im
)),
7840 0.5 * atan2 (im
, re
));
7843 else if (SCM_NUMBERP (z
))
7845 double xx
= scm_to_double (z
);
7847 return scm_c_make_rectangular (0.0, sqrt (-xx
));
7849 return scm_from_double (sqrt (xx
));
7852 SCM_WTA_DISPATCH_1 (g_scm_sqrt
, z
, 1, s_scm_sqrt
);
7863 mpz_init_set_si (z_negative_one
, -1);
7865 /* It may be possible to tune the performance of some algorithms by using
7866 * the following constants to avoid the creation of bignums. Please, before
7867 * using these values, remember the two rules of program optimization:
7868 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
7869 scm_c_define ("most-positive-fixnum",
7870 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
7871 scm_c_define ("most-negative-fixnum",
7872 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
7874 scm_add_feature ("complex");
7875 scm_add_feature ("inexact");
7876 flo0
= scm_from_double (0.0);
7878 /* determine floating point precision */
7879 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
7881 init_dblprec(&scm_dblprec
[i
-2],i
);
7882 init_fx_radix(fx_per_radix
[i
-2],i
);
7885 /* hard code precision for base 10 if the preprocessor tells us to... */
7886 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
7889 exactly_one_half
= scm_divide (SCM_INUM1
, SCM_I_MAKINUM (2));
7890 #include "libguile/numbers.x"