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_DEFINE (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_WRONG_TYPE_ARG (1, x
);
517 SCM_DEFINE (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_WRONG_TYPE_ARG (1, x
);
533 SCM_DEFINE (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_is_true (scm_inf_p (n
)))
551 SCM_WRONG_TYPE_ARG (1, n
);
552 else if (SCM_REALP (n
))
554 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
560 SCM_WRONG_TYPE_ARG (1, n
);
563 SCM_WRONG_TYPE_ARG (1, n
);
568 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
570 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
572 #define FUNC_NAME s_scm_even_p
576 scm_t_inum val
= SCM_I_INUM (n
);
577 return scm_from_bool ((val
& 1L) == 0);
579 else if (SCM_BIGP (n
))
581 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
582 scm_remember_upto_here_1 (n
);
583 return scm_from_bool (even_p
);
585 else if (scm_is_true (scm_inf_p (n
)))
586 SCM_WRONG_TYPE_ARG (1, n
);
587 else if (SCM_REALP (n
))
589 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
595 SCM_WRONG_TYPE_ARG (1, n
);
598 SCM_WRONG_TYPE_ARG (1, n
);
602 SCM_DEFINE (scm_finite_p
, "finite?", 1, 0, 0,
604 "Return @code{#t} if the real number @var{x} is neither\n"
605 "infinite nor a NaN, @code{#f} otherwise.")
606 #define FUNC_NAME s_scm_finite_p
609 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
610 else if (scm_is_real (x
))
613 SCM_WRONG_TYPE_ARG (1, x
);
617 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
619 "Return @code{#t} if the real number @var{x} is @samp{+inf.0} or\n"
620 "@samp{-inf.0}. Otherwise return @code{#f}.")
621 #define FUNC_NAME s_scm_inf_p
624 return scm_from_bool (isinf (SCM_REAL_VALUE (x
)));
625 else if (scm_is_real (x
))
628 SCM_WRONG_TYPE_ARG (1, x
);
632 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
634 "Return @code{#t} if the real number @var{x} is a NaN,\n"
635 "or @code{#f} otherwise.")
636 #define FUNC_NAME s_scm_nan_p
639 return scm_from_bool (isnan (SCM_REAL_VALUE (x
)));
640 else if (scm_is_real (x
))
643 SCM_WRONG_TYPE_ARG (1, x
);
647 /* Guile's idea of infinity. */
648 static double guile_Inf
;
650 /* Guile's idea of not a number. */
651 static double guile_NaN
;
654 guile_ieee_init (void)
656 /* Some version of gcc on some old version of Linux used to crash when
657 trying to make Inf and NaN. */
660 /* C99 INFINITY, when available.
661 FIXME: The standard allows for INFINITY to be something that overflows
662 at compile time. We ought to have a configure test to check for that
663 before trying to use it. (But in practice we believe this is not a
664 problem on any system guile is likely to target.) */
665 guile_Inf
= INFINITY
;
666 #elif defined HAVE_DINFINITY
668 extern unsigned int DINFINITY
[2];
669 guile_Inf
= (*((double *) (DINFINITY
)));
676 if (guile_Inf
== tmp
)
683 /* C99 NAN, when available */
685 #elif defined HAVE_DQNAN
688 extern unsigned int DQNAN
[2];
689 guile_NaN
= (*((double *)(DQNAN
)));
692 guile_NaN
= guile_Inf
/ guile_Inf
;
696 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
699 #define FUNC_NAME s_scm_inf
701 static int initialized
= 0;
707 return scm_from_double (guile_Inf
);
711 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
714 #define FUNC_NAME s_scm_nan
716 static int initialized
= 0;
722 return scm_from_double (guile_NaN
);
727 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
729 "Return the absolute value of @var{x}.")
734 scm_t_inum xx
= SCM_I_INUM (x
);
737 else if (SCM_POSFIXABLE (-xx
))
738 return SCM_I_MAKINUM (-xx
);
740 return scm_i_inum2big (-xx
);
742 else if (SCM_BIGP (x
))
744 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
746 return scm_i_clonebig (x
, 0);
750 else if (SCM_REALP (x
))
752 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
753 double xx
= SCM_REAL_VALUE (x
);
755 return scm_from_double (-xx
);
759 else if (SCM_FRACTIONP (x
))
761 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
763 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
764 SCM_FRACTION_DENOMINATOR (x
));
767 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
772 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
773 /* "Return the quotient of the numbers @var{x} and @var{y}."
776 scm_quotient (SCM x
, SCM y
)
778 if (SCM_LIKELY (SCM_I_INUMP (x
)))
780 scm_t_inum xx
= SCM_I_INUM (x
);
781 if (SCM_LIKELY (SCM_I_INUMP (y
)))
783 scm_t_inum yy
= SCM_I_INUM (y
);
784 if (SCM_UNLIKELY (yy
== 0))
785 scm_num_overflow (s_quotient
);
788 scm_t_inum z
= xx
/ yy
;
789 if (SCM_LIKELY (SCM_FIXABLE (z
)))
790 return SCM_I_MAKINUM (z
);
792 return scm_i_inum2big (z
);
795 else if (SCM_BIGP (y
))
797 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
798 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
799 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
801 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
802 scm_remember_upto_here_1 (y
);
803 return SCM_I_MAKINUM (-1);
809 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
811 else if (SCM_BIGP (x
))
813 if (SCM_LIKELY (SCM_I_INUMP (y
)))
815 scm_t_inum yy
= SCM_I_INUM (y
);
816 if (SCM_UNLIKELY (yy
== 0))
817 scm_num_overflow (s_quotient
);
818 else if (SCM_UNLIKELY (yy
== 1))
822 SCM result
= scm_i_mkbig ();
825 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
828 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
831 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
832 scm_remember_upto_here_1 (x
);
833 return scm_i_normbig (result
);
836 else if (SCM_BIGP (y
))
838 SCM result
= scm_i_mkbig ();
839 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
842 scm_remember_upto_here_2 (x
, y
);
843 return scm_i_normbig (result
);
846 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
849 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
852 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
853 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
855 * "(remainder 13 4) @result{} 1\n"
856 * "(remainder -13 4) @result{} -1\n"
860 scm_remainder (SCM x
, SCM y
)
862 if (SCM_LIKELY (SCM_I_INUMP (x
)))
864 if (SCM_LIKELY (SCM_I_INUMP (y
)))
866 scm_t_inum yy
= SCM_I_INUM (y
);
867 if (SCM_UNLIKELY (yy
== 0))
868 scm_num_overflow (s_remainder
);
871 /* C99 specifies that "%" is the remainder corresponding to a
872 quotient rounded towards zero, and that's also traditional
873 for machine division, so z here should be well defined. */
874 scm_t_inum z
= SCM_I_INUM (x
) % yy
;
875 return SCM_I_MAKINUM (z
);
878 else if (SCM_BIGP (y
))
880 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
881 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
882 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
884 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
885 scm_remember_upto_here_1 (y
);
892 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
894 else if (SCM_BIGP (x
))
896 if (SCM_LIKELY (SCM_I_INUMP (y
)))
898 scm_t_inum yy
= SCM_I_INUM (y
);
899 if (SCM_UNLIKELY (yy
== 0))
900 scm_num_overflow (s_remainder
);
903 SCM result
= scm_i_mkbig ();
906 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
907 scm_remember_upto_here_1 (x
);
908 return scm_i_normbig (result
);
911 else if (SCM_BIGP (y
))
913 SCM result
= scm_i_mkbig ();
914 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
917 scm_remember_upto_here_2 (x
, y
);
918 return scm_i_normbig (result
);
921 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
924 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
928 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
929 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
931 * "(modulo 13 4) @result{} 1\n"
932 * "(modulo -13 4) @result{} 3\n"
936 scm_modulo (SCM x
, SCM y
)
938 if (SCM_LIKELY (SCM_I_INUMP (x
)))
940 scm_t_inum xx
= SCM_I_INUM (x
);
941 if (SCM_LIKELY (SCM_I_INUMP (y
)))
943 scm_t_inum yy
= SCM_I_INUM (y
);
944 if (SCM_UNLIKELY (yy
== 0))
945 scm_num_overflow (s_modulo
);
948 /* C99 specifies that "%" is the remainder corresponding to a
949 quotient rounded towards zero, and that's also traditional
950 for machine division, so z here should be well defined. */
951 scm_t_inum z
= xx
% yy
;
968 return SCM_I_MAKINUM (result
);
971 else if (SCM_BIGP (y
))
973 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
980 SCM pos_y
= scm_i_clonebig (y
, 0);
981 /* do this after the last scm_op */
982 mpz_init_set_si (z_x
, xx
);
983 result
= pos_y
; /* re-use this bignum */
984 mpz_mod (SCM_I_BIG_MPZ (result
),
986 SCM_I_BIG_MPZ (pos_y
));
987 scm_remember_upto_here_1 (pos_y
);
991 result
= scm_i_mkbig ();
992 /* do this after the last scm_op */
993 mpz_init_set_si (z_x
, xx
);
994 mpz_mod (SCM_I_BIG_MPZ (result
),
997 scm_remember_upto_here_1 (y
);
1000 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1001 mpz_add (SCM_I_BIG_MPZ (result
),
1003 SCM_I_BIG_MPZ (result
));
1004 scm_remember_upto_here_1 (y
);
1005 /* and do this before the next one */
1007 return scm_i_normbig (result
);
1011 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1013 else if (SCM_BIGP (x
))
1015 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1017 scm_t_inum yy
= SCM_I_INUM (y
);
1018 if (SCM_UNLIKELY (yy
== 0))
1019 scm_num_overflow (s_modulo
);
1022 SCM result
= scm_i_mkbig ();
1023 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
1025 (yy
< 0) ? - yy
: yy
);
1026 scm_remember_upto_here_1 (x
);
1027 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1028 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
1029 SCM_I_BIG_MPZ (result
),
1031 return scm_i_normbig (result
);
1034 else if (SCM_BIGP (y
))
1036 SCM result
= scm_i_mkbig ();
1037 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
1038 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1039 mpz_mod (SCM_I_BIG_MPZ (result
),
1041 SCM_I_BIG_MPZ (pos_y
));
1043 scm_remember_upto_here_1 (x
);
1044 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1045 mpz_add (SCM_I_BIG_MPZ (result
),
1047 SCM_I_BIG_MPZ (result
));
1048 scm_remember_upto_here_2 (y
, pos_y
);
1049 return scm_i_normbig (result
);
1052 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1055 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1058 static SCM
scm_i_inexact_euclidean_quotient (double x
, double y
);
1059 static SCM
scm_i_slow_exact_euclidean_quotient (SCM x
, SCM y
);
1061 SCM_PRIMITIVE_GENERIC (scm_euclidean_quotient
, "euclidean-quotient", 2, 0, 0,
1063 "Return the integer @var{q} such that\n"
1064 "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1065 "where @math{0 <= @var{r} < abs(@var{y})}.\n"
1067 "(euclidean-quotient 123 10) @result{} 12\n"
1068 "(euclidean-quotient 123 -10) @result{} -12\n"
1069 "(euclidean-quotient -123 10) @result{} -13\n"
1070 "(euclidean-quotient -123 -10) @result{} 13\n"
1071 "(euclidean-quotient -123.2 -63.5) @result{} 2.0\n"
1072 "(euclidean-quotient 16/3 -10/7) @result{} -3\n"
1074 #define FUNC_NAME s_scm_euclidean_quotient
1076 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1078 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1080 scm_t_inum yy
= SCM_I_INUM (y
);
1081 if (SCM_UNLIKELY (yy
== 0))
1082 scm_num_overflow (s_scm_euclidean_quotient
);
1085 scm_t_inum xx
= SCM_I_INUM (x
);
1086 scm_t_inum qq
= xx
/ yy
;
1094 return SCM_I_MAKINUM (qq
);
1097 else if (SCM_BIGP (y
))
1099 if (SCM_I_INUM (x
) >= 0)
1102 return SCM_I_MAKINUM (- mpz_sgn (SCM_I_BIG_MPZ (y
)));
1104 else if (SCM_REALP (y
))
1105 return scm_i_inexact_euclidean_quotient
1106 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1107 else if (SCM_FRACTIONP (y
))
1108 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1110 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1111 s_scm_euclidean_quotient
);
1113 else if (SCM_BIGP (x
))
1115 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1117 scm_t_inum yy
= SCM_I_INUM (y
);
1118 if (SCM_UNLIKELY (yy
== 0))
1119 scm_num_overflow (s_scm_euclidean_quotient
);
1122 SCM q
= scm_i_mkbig ();
1124 mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (x
), yy
);
1127 mpz_fdiv_q_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (x
), -yy
);
1128 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1130 scm_remember_upto_here_1 (x
);
1131 return scm_i_normbig (q
);
1134 else if (SCM_BIGP (y
))
1136 SCM q
= scm_i_mkbig ();
1137 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1138 mpz_fdiv_q (SCM_I_BIG_MPZ (q
),
1142 mpz_cdiv_q (SCM_I_BIG_MPZ (q
),
1145 scm_remember_upto_here_2 (x
, y
);
1146 return scm_i_normbig (q
);
1148 else if (SCM_REALP (y
))
1149 return scm_i_inexact_euclidean_quotient
1150 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1151 else if (SCM_FRACTIONP (y
))
1152 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1154 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1155 s_scm_euclidean_quotient
);
1157 else if (SCM_REALP (x
))
1159 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1160 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1161 return scm_i_inexact_euclidean_quotient
1162 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1164 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1165 s_scm_euclidean_quotient
);
1167 else if (SCM_FRACTIONP (x
))
1170 return scm_i_inexact_euclidean_quotient
1171 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1173 return scm_i_slow_exact_euclidean_quotient (x
, y
);
1176 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG1
,
1177 s_scm_euclidean_quotient
);
1182 scm_i_inexact_euclidean_quotient (double x
, double y
)
1184 if (SCM_LIKELY (y
> 0))
1185 return scm_from_double (floor (x
/ y
));
1186 else if (SCM_LIKELY (y
< 0))
1187 return scm_from_double (ceil (x
/ y
));
1189 scm_num_overflow (s_scm_euclidean_quotient
); /* or return a NaN? */
1194 /* Compute exact euclidean_quotient the slow way.
1195 We use this only if both arguments are exact,
1196 and at least one of them is a fraction */
1198 scm_i_slow_exact_euclidean_quotient (SCM x
, SCM y
)
1200 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1201 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG1
,
1202 s_scm_euclidean_quotient
);
1203 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1204 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quotient
, x
, y
, SCM_ARG2
,
1205 s_scm_euclidean_quotient
);
1206 else if (scm_is_true (scm_positive_p (y
)))
1207 return scm_floor (scm_divide (x
, y
));
1208 else if (scm_is_true (scm_negative_p (y
)))
1209 return scm_ceiling (scm_divide (x
, y
));
1211 scm_num_overflow (s_scm_euclidean_quotient
);
1214 static SCM
scm_i_inexact_euclidean_remainder (double x
, double y
);
1215 static SCM
scm_i_slow_exact_euclidean_remainder (SCM x
, SCM y
);
1217 SCM_PRIMITIVE_GENERIC (scm_euclidean_remainder
, "euclidean-remainder", 2, 0, 0,
1219 "Return the real number @var{r} such that\n"
1220 "@math{0 <= @var{r} < abs(@var{y})} and\n"
1221 "@math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1222 "for some integer @var{q}.\n"
1224 "(euclidean-remainder 123 10) @result{} 3\n"
1225 "(euclidean-remainder 123 -10) @result{} 3\n"
1226 "(euclidean-remainder -123 10) @result{} 7\n"
1227 "(euclidean-remainder -123 -10) @result{} 7\n"
1228 "(euclidean-remainder -123.2 -63.5) @result{} 3.8\n"
1229 "(euclidean-remainder 16/3 -10/7) @result{} 22/21\n"
1231 #define FUNC_NAME s_scm_euclidean_remainder
1233 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1235 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1237 scm_t_inum yy
= SCM_I_INUM (y
);
1238 if (SCM_UNLIKELY (yy
== 0))
1239 scm_num_overflow (s_scm_euclidean_remainder
);
1242 scm_t_inum rr
= SCM_I_INUM (x
) % yy
;
1244 return SCM_I_MAKINUM (rr
);
1246 return SCM_I_MAKINUM (rr
+ yy
);
1248 return SCM_I_MAKINUM (rr
- yy
);
1251 else if (SCM_BIGP (y
))
1253 scm_t_inum xx
= SCM_I_INUM (x
);
1256 else if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1258 SCM r
= scm_i_mkbig ();
1259 mpz_sub_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1260 scm_remember_upto_here_1 (y
);
1261 return scm_i_normbig (r
);
1265 SCM r
= scm_i_mkbig ();
1266 mpz_add_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1267 scm_remember_upto_here_1 (y
);
1268 mpz_neg (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (r
));
1269 return scm_i_normbig (r
);
1272 else if (SCM_REALP (y
))
1273 return scm_i_inexact_euclidean_remainder
1274 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1275 else if (SCM_FRACTIONP (y
))
1276 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1278 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1279 s_scm_euclidean_remainder
);
1281 else if (SCM_BIGP (x
))
1283 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1285 scm_t_inum yy
= SCM_I_INUM (y
);
1286 if (SCM_UNLIKELY (yy
== 0))
1287 scm_num_overflow (s_scm_euclidean_remainder
);
1293 rr
= mpz_fdiv_ui (SCM_I_BIG_MPZ (x
), yy
);
1294 scm_remember_upto_here_1 (x
);
1295 return SCM_I_MAKINUM (rr
);
1298 else if (SCM_BIGP (y
))
1300 SCM r
= scm_i_mkbig ();
1301 mpz_mod (SCM_I_BIG_MPZ (r
),
1304 scm_remember_upto_here_2 (x
, y
);
1305 return scm_i_normbig (r
);
1307 else if (SCM_REALP (y
))
1308 return scm_i_inexact_euclidean_remainder
1309 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1310 else if (SCM_FRACTIONP (y
))
1311 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1313 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1314 s_scm_euclidean_remainder
);
1316 else if (SCM_REALP (x
))
1318 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1319 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1320 return scm_i_inexact_euclidean_remainder
1321 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1323 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1324 s_scm_euclidean_remainder
);
1326 else if (SCM_FRACTIONP (x
))
1329 return scm_i_inexact_euclidean_remainder
1330 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1332 return scm_i_slow_exact_euclidean_remainder (x
, y
);
1335 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG1
,
1336 s_scm_euclidean_remainder
);
1341 scm_i_inexact_euclidean_remainder (double x
, double y
)
1345 /* Although it would be more efficient to use fmod here, we can't
1346 because it would in some cases produce results inconsistent with
1347 scm_i_inexact_euclidean_quotient, such that x != q * y + r (not
1348 even close). In particular, when x is very close to a multiple of
1349 y, then r might be either 0.0 or abs(y)-epsilon, but those two
1350 cases must correspond to different choices of q. If r = 0.0 then q
1351 must be x/y, and if r = abs(y) then q must be (x-r)/y. If quotient
1352 chooses one and remainder chooses the other, it would be bad. This
1353 problem was observed with x = 130.0 and y = 10/7. */
1354 if (SCM_LIKELY (y
> 0))
1356 else if (SCM_LIKELY (y
< 0))
1359 scm_num_overflow (s_scm_euclidean_remainder
); /* or return a NaN? */
1362 return scm_from_double (x
- q
* y
);
1365 /* Compute exact euclidean_remainder the slow way.
1366 We use this only if both arguments are exact,
1367 and at least one of them is a fraction */
1369 scm_i_slow_exact_euclidean_remainder (SCM x
, SCM y
)
1373 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1374 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG1
,
1375 s_scm_euclidean_remainder
);
1376 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1377 SCM_WTA_DISPATCH_2 (g_scm_euclidean_remainder
, x
, y
, SCM_ARG2
,
1378 s_scm_euclidean_remainder
);
1379 else if (scm_is_true (scm_positive_p (y
)))
1380 q
= scm_floor (scm_divide (x
, y
));
1381 else if (scm_is_true (scm_negative_p (y
)))
1382 q
= scm_ceiling (scm_divide (x
, y
));
1384 scm_num_overflow (s_scm_euclidean_remainder
);
1385 return scm_difference (x
, scm_product (y
, q
));
1389 static SCM
scm_i_inexact_euclidean_quo_and_rem (double x
, double y
);
1390 static SCM
scm_i_slow_exact_euclidean_quo_and_rem (SCM x
, SCM y
);
1392 SCM_PRIMITIVE_GENERIC (scm_euclidean_quo_and_rem
, "euclidean/", 2, 0, 0,
1394 "Return the integer @var{q} and the real number @var{r}\n"
1395 "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1396 "and @math{0 <= @var{r} < abs(@var{y})}.\n"
1398 "(euclidean/ 123 10) @result{} 12 and 3\n"
1399 "(euclidean/ 123 -10) @result{} -12 and 3\n"
1400 "(euclidean/ -123 10) @result{} -13 and 7\n"
1401 "(euclidean/ -123 -10) @result{} 13 and 7\n"
1402 "(euclidean/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
1403 "(euclidean/ 16/3 -10/7) @result{} -3 and 22/21\n"
1405 #define FUNC_NAME s_scm_euclidean_quo_and_rem
1407 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1409 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1411 scm_t_inum yy
= SCM_I_INUM (y
);
1412 if (SCM_UNLIKELY (yy
== 0))
1413 scm_num_overflow (s_scm_euclidean_quo_and_rem
);
1416 scm_t_inum xx
= SCM_I_INUM (x
);
1417 scm_t_inum qq
= xx
/ yy
;
1418 scm_t_inum rr
= xx
- qq
* yy
;
1426 return scm_values (scm_list_2 (SCM_I_MAKINUM (qq
),
1427 SCM_I_MAKINUM (rr
)));
1430 else if (SCM_BIGP (y
))
1432 scm_t_inum xx
= SCM_I_INUM (x
);
1434 return scm_values (scm_list_2 (SCM_INUM0
, x
));
1435 else if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1437 SCM r
= scm_i_mkbig ();
1438 mpz_sub_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1439 scm_remember_upto_here_1 (y
);
1441 (scm_list_2 (SCM_I_MAKINUM (-1), scm_i_normbig (r
)));
1445 SCM r
= scm_i_mkbig ();
1446 mpz_add_ui (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (y
), -xx
);
1447 scm_remember_upto_here_1 (y
);
1448 mpz_neg (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (r
));
1449 return scm_values (scm_list_2 (SCM_INUM1
, scm_i_normbig (r
)));
1452 else if (SCM_REALP (y
))
1453 return scm_i_inexact_euclidean_quo_and_rem
1454 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1455 else if (SCM_FRACTIONP (y
))
1456 return scm_i_slow_exact_euclidean_quo_and_rem (x
, y
);
1458 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem
, x
, y
, SCM_ARG2
,
1459 s_scm_euclidean_quo_and_rem
);
1461 else if (SCM_BIGP (x
))
1463 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1465 scm_t_inum yy
= SCM_I_INUM (y
);
1466 if (SCM_UNLIKELY (yy
== 0))
1467 scm_num_overflow (s_scm_euclidean_quo_and_rem
);
1470 SCM q
= scm_i_mkbig ();
1471 SCM r
= scm_i_mkbig ();
1473 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1474 SCM_I_BIG_MPZ (x
), yy
);
1477 mpz_fdiv_qr_ui (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1478 SCM_I_BIG_MPZ (x
), -yy
);
1479 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1481 scm_remember_upto_here_1 (x
);
1482 return scm_values (scm_list_2 (scm_i_normbig (q
),
1483 scm_i_normbig (r
)));
1486 else if (SCM_BIGP (y
))
1488 SCM q
= scm_i_mkbig ();
1489 SCM r
= scm_i_mkbig ();
1490 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1491 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1492 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1494 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1495 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1496 scm_remember_upto_here_2 (x
, y
);
1497 return scm_values (scm_list_2 (scm_i_normbig (q
),
1498 scm_i_normbig (r
)));
1500 else if (SCM_REALP (y
))
1501 return scm_i_inexact_euclidean_quo_and_rem
1502 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1503 else if (SCM_FRACTIONP (y
))
1504 return scm_i_slow_exact_euclidean_quo_and_rem (x
, y
);
1506 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem
, x
, y
, SCM_ARG2
,
1507 s_scm_euclidean_quo_and_rem
);
1509 else if (SCM_REALP (x
))
1511 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1512 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1513 return scm_i_inexact_euclidean_quo_and_rem
1514 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1516 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem
, x
, y
, SCM_ARG2
,
1517 s_scm_euclidean_quo_and_rem
);
1519 else if (SCM_FRACTIONP (x
))
1522 return scm_i_inexact_euclidean_quo_and_rem
1523 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1525 return scm_i_slow_exact_euclidean_quo_and_rem (x
, y
);
1528 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem
, x
, y
, SCM_ARG1
,
1529 s_scm_euclidean_quo_and_rem
);
1534 scm_i_inexact_euclidean_quo_and_rem (double x
, double y
)
1538 if (SCM_LIKELY (y
> 0))
1540 else if (SCM_LIKELY (y
< 0))
1543 scm_num_overflow (s_scm_euclidean_quo_and_rem
); /* or return a NaN? */
1547 return scm_values (scm_list_2 (scm_from_double (q
),
1548 scm_from_double (r
)));
1551 /* Compute exact euclidean quotient and remainder the slow way.
1552 We use this only if both arguments are exact,
1553 and at least one of them is a fraction */
1555 scm_i_slow_exact_euclidean_quo_and_rem (SCM x
, SCM y
)
1559 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1560 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem
, x
, y
, SCM_ARG1
,
1561 s_scm_euclidean_quo_and_rem
);
1562 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1563 SCM_WTA_DISPATCH_2 (g_scm_euclidean_quo_and_rem
, x
, y
, SCM_ARG2
,
1564 s_scm_euclidean_quo_and_rem
);
1565 else if (scm_is_true (scm_positive_p (y
)))
1566 q
= scm_floor (scm_divide (x
, y
));
1567 else if (scm_is_true (scm_negative_p (y
)))
1568 q
= scm_ceiling (scm_divide (x
, y
));
1570 scm_num_overflow (s_scm_euclidean_quo_and_rem
);
1571 r
= scm_difference (x
, scm_product (q
, y
));
1572 return scm_values (scm_list_2 (q
, r
));
1575 static SCM
scm_i_inexact_centered_quotient (double x
, double y
);
1576 static SCM
scm_i_bigint_centered_quotient (SCM x
, SCM y
);
1577 static SCM
scm_i_slow_exact_centered_quotient (SCM x
, SCM y
);
1579 SCM_PRIMITIVE_GENERIC (scm_centered_quotient
, "centered-quotient", 2, 0, 0,
1581 "Return the integer @var{q} such that\n"
1582 "@math{@var{x} = @var{q}*@var{y} + @var{r}} where\n"
1583 "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
1585 "(centered-quotient 123 10) @result{} 12\n"
1586 "(centered-quotient 123 -10) @result{} -12\n"
1587 "(centered-quotient -123 10) @result{} -12\n"
1588 "(centered-quotient -123 -10) @result{} 12\n"
1589 "(centered-quotient -123.2 -63.5) @result{} 2.0\n"
1590 "(centered-quotient 16/3 -10/7) @result{} -4\n"
1592 #define FUNC_NAME s_scm_centered_quotient
1594 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1596 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1598 scm_t_inum yy
= SCM_I_INUM (y
);
1599 if (SCM_UNLIKELY (yy
== 0))
1600 scm_num_overflow (s_scm_centered_quotient
);
1603 scm_t_inum xx
= SCM_I_INUM (x
);
1604 scm_t_inum qq
= xx
/ yy
;
1605 scm_t_inum rr
= xx
- qq
* yy
;
1606 if (SCM_LIKELY (xx
> 0))
1608 if (SCM_LIKELY (yy
> 0))
1610 if (rr
>= (yy
+ 1) / 2)
1615 if (rr
>= (1 - yy
) / 2)
1621 if (SCM_LIKELY (yy
> 0))
1632 return SCM_I_MAKINUM (qq
);
1635 else if (SCM_BIGP (y
))
1637 /* Pass a denormalized bignum version of x (even though it
1638 can fit in a fixnum) to scm_i_bigint_centered_quotient */
1639 return scm_i_bigint_centered_quotient
1640 (scm_i_long2big (SCM_I_INUM (x
)), y
);
1642 else if (SCM_REALP (y
))
1643 return scm_i_inexact_centered_quotient
1644 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1645 else if (SCM_FRACTIONP (y
))
1646 return scm_i_slow_exact_centered_quotient (x
, y
);
1648 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1649 s_scm_centered_quotient
);
1651 else if (SCM_BIGP (x
))
1653 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1655 scm_t_inum yy
= SCM_I_INUM (y
);
1656 if (SCM_UNLIKELY (yy
== 0))
1657 scm_num_overflow (s_scm_centered_quotient
);
1660 SCM q
= scm_i_mkbig ();
1662 /* Arrange for rr to initially be non-positive,
1663 because that simplifies the test to see
1664 if it is within the needed bounds. */
1667 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
1668 SCM_I_BIG_MPZ (x
), yy
);
1669 scm_remember_upto_here_1 (x
);
1671 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
1672 SCM_I_BIG_MPZ (q
), 1);
1676 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
1677 SCM_I_BIG_MPZ (x
), -yy
);
1678 scm_remember_upto_here_1 (x
);
1679 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
1681 mpz_add_ui (SCM_I_BIG_MPZ (q
),
1682 SCM_I_BIG_MPZ (q
), 1);
1684 return scm_i_normbig (q
);
1687 else if (SCM_BIGP (y
))
1688 return scm_i_bigint_centered_quotient (x
, y
);
1689 else if (SCM_REALP (y
))
1690 return scm_i_inexact_centered_quotient
1691 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1692 else if (SCM_FRACTIONP (y
))
1693 return scm_i_slow_exact_centered_quotient (x
, y
);
1695 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1696 s_scm_centered_quotient
);
1698 else if (SCM_REALP (x
))
1700 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1701 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1702 return scm_i_inexact_centered_quotient
1703 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1705 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1706 s_scm_centered_quotient
);
1708 else if (SCM_FRACTIONP (x
))
1711 return scm_i_inexact_centered_quotient
1712 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1714 return scm_i_slow_exact_centered_quotient (x
, y
);
1717 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG1
,
1718 s_scm_centered_quotient
);
1723 scm_i_inexact_centered_quotient (double x
, double y
)
1725 if (SCM_LIKELY (y
> 0))
1726 return scm_from_double (floor (x
/y
+ 0.5));
1727 else if (SCM_LIKELY (y
< 0))
1728 return scm_from_double (ceil (x
/y
- 0.5));
1730 scm_num_overflow (s_scm_centered_quotient
); /* or return a NaN? */
1735 /* Assumes that both x and y are bigints, though
1736 x might be able to fit into a fixnum. */
1738 scm_i_bigint_centered_quotient (SCM x
, SCM y
)
1742 /* Note that x might be small enough to fit into a
1743 fixnum, so we must not let it escape into the wild */
1747 /* min_r will eventually become -abs(y)/2 */
1748 min_r
= scm_i_mkbig ();
1749 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
1750 SCM_I_BIG_MPZ (y
), 1);
1752 /* Arrange for rr to initially be non-positive,
1753 because that simplifies the test to see
1754 if it is within the needed bounds. */
1755 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1757 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1758 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1759 scm_remember_upto_here_2 (x
, y
);
1760 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
1761 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1762 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
1763 SCM_I_BIG_MPZ (q
), 1);
1767 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
1768 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1769 scm_remember_upto_here_2 (x
, y
);
1770 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1771 mpz_add_ui (SCM_I_BIG_MPZ (q
),
1772 SCM_I_BIG_MPZ (q
), 1);
1774 scm_remember_upto_here_2 (r
, min_r
);
1775 return scm_i_normbig (q
);
1778 /* Compute exact centered quotient the slow way.
1779 We use this only if both arguments are exact,
1780 and at least one of them is a fraction */
1782 scm_i_slow_exact_centered_quotient (SCM x
, SCM y
)
1784 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
1785 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG1
,
1786 s_scm_centered_quotient
);
1787 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
1788 SCM_WTA_DISPATCH_2 (g_scm_centered_quotient
, x
, y
, SCM_ARG2
,
1789 s_scm_centered_quotient
);
1790 else if (scm_is_true (scm_positive_p (y
)))
1791 return scm_floor (scm_sum (scm_divide (x
, y
),
1793 else if (scm_is_true (scm_negative_p (y
)))
1794 return scm_ceiling (scm_difference (scm_divide (x
, y
),
1797 scm_num_overflow (s_scm_centered_quotient
);
1800 static SCM
scm_i_inexact_centered_remainder (double x
, double y
);
1801 static SCM
scm_i_bigint_centered_remainder (SCM x
, SCM y
);
1802 static SCM
scm_i_slow_exact_centered_remainder (SCM x
, SCM y
);
1804 SCM_PRIMITIVE_GENERIC (scm_centered_remainder
, "centered-remainder", 2, 0, 0,
1806 "Return the real number @var{r} such that\n"
1807 "@math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}\n"
1808 "and @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
1809 "for some integer @var{q}.\n"
1811 "(centered-remainder 123 10) @result{} 3\n"
1812 "(centered-remainder 123 -10) @result{} 3\n"
1813 "(centered-remainder -123 10) @result{} -3\n"
1814 "(centered-remainder -123 -10) @result{} -3\n"
1815 "(centered-remainder -123.2 -63.5) @result{} 3.8\n"
1816 "(centered-remainder 16/3 -10/7) @result{} -8/21\n"
1818 #define FUNC_NAME s_scm_centered_remainder
1820 if (SCM_LIKELY (SCM_I_INUMP (x
)))
1822 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1824 scm_t_inum yy
= SCM_I_INUM (y
);
1825 if (SCM_UNLIKELY (yy
== 0))
1826 scm_num_overflow (s_scm_centered_remainder
);
1829 scm_t_inum xx
= SCM_I_INUM (x
);
1830 scm_t_inum rr
= xx
% yy
;
1831 if (SCM_LIKELY (xx
> 0))
1833 if (SCM_LIKELY (yy
> 0))
1835 if (rr
>= (yy
+ 1) / 2)
1840 if (rr
>= (1 - yy
) / 2)
1846 if (SCM_LIKELY (yy
> 0))
1857 return SCM_I_MAKINUM (rr
);
1860 else if (SCM_BIGP (y
))
1862 /* Pass a denormalized bignum version of x (even though it
1863 can fit in a fixnum) to scm_i_bigint_centered_remainder */
1864 return scm_i_bigint_centered_remainder
1865 (scm_i_long2big (SCM_I_INUM (x
)), y
);
1867 else if (SCM_REALP (y
))
1868 return scm_i_inexact_centered_remainder
1869 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
1870 else if (SCM_FRACTIONP (y
))
1871 return scm_i_slow_exact_centered_remainder (x
, y
);
1873 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1874 s_scm_centered_remainder
);
1876 else if (SCM_BIGP (x
))
1878 if (SCM_LIKELY (SCM_I_INUMP (y
)))
1880 scm_t_inum yy
= SCM_I_INUM (y
);
1881 if (SCM_UNLIKELY (yy
== 0))
1882 scm_num_overflow (s_scm_centered_remainder
);
1886 /* Arrange for rr to initially be non-positive,
1887 because that simplifies the test to see
1888 if it is within the needed bounds. */
1891 rr
= - mpz_cdiv_ui (SCM_I_BIG_MPZ (x
), yy
);
1892 scm_remember_upto_here_1 (x
);
1898 rr
= - mpz_cdiv_ui (SCM_I_BIG_MPZ (x
), -yy
);
1899 scm_remember_upto_here_1 (x
);
1903 return SCM_I_MAKINUM (rr
);
1906 else if (SCM_BIGP (y
))
1907 return scm_i_bigint_centered_remainder (x
, y
);
1908 else if (SCM_REALP (y
))
1909 return scm_i_inexact_centered_remainder
1910 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
1911 else if (SCM_FRACTIONP (y
))
1912 return scm_i_slow_exact_centered_remainder (x
, y
);
1914 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1915 s_scm_centered_remainder
);
1917 else if (SCM_REALP (x
))
1919 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
1920 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
1921 return scm_i_inexact_centered_remainder
1922 (SCM_REAL_VALUE (x
), scm_to_double (y
));
1924 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
1925 s_scm_centered_remainder
);
1927 else if (SCM_FRACTIONP (x
))
1930 return scm_i_inexact_centered_remainder
1931 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
1933 return scm_i_slow_exact_centered_remainder (x
, y
);
1936 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG1
,
1937 s_scm_centered_remainder
);
1942 scm_i_inexact_centered_remainder (double x
, double y
)
1946 /* Although it would be more efficient to use fmod here, we can't
1947 because it would in some cases produce results inconsistent with
1948 scm_i_inexact_centered_quotient, such that x != r + q * y (not even
1949 close). In particular, when x-y/2 is very close to a multiple of
1950 y, then r might be either -abs(y/2) or abs(y/2)-epsilon, but those
1951 two cases must correspond to different choices of q. If quotient
1952 chooses one and remainder chooses the other, it would be bad. */
1953 if (SCM_LIKELY (y
> 0))
1954 q
= floor (x
/y
+ 0.5);
1955 else if (SCM_LIKELY (y
< 0))
1956 q
= ceil (x
/y
- 0.5);
1958 scm_num_overflow (s_scm_centered_remainder
); /* or return a NaN? */
1961 return scm_from_double (x
- q
* y
);
1964 /* Assumes that both x and y are bigints, though
1965 x might be able to fit into a fixnum. */
1967 scm_i_bigint_centered_remainder (SCM x
, SCM y
)
1971 /* Note that x might be small enough to fit into a
1972 fixnum, so we must not let it escape into the wild */
1975 /* min_r will eventually become -abs(y)/2 */
1976 min_r
= scm_i_mkbig ();
1977 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
1978 SCM_I_BIG_MPZ (y
), 1);
1980 /* Arrange for rr to initially be non-positive,
1981 because that simplifies the test to see
1982 if it is within the needed bounds. */
1983 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
1985 mpz_cdiv_r (SCM_I_BIG_MPZ (r
),
1986 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1987 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
1988 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1989 mpz_add (SCM_I_BIG_MPZ (r
),
1995 mpz_fdiv_r (SCM_I_BIG_MPZ (r
),
1996 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
1997 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
1998 mpz_sub (SCM_I_BIG_MPZ (r
),
2002 scm_remember_upto_here_2 (x
, y
);
2003 return scm_i_normbig (r
);
2006 /* Compute exact centered_remainder the slow way.
2007 We use this only if both arguments are exact,
2008 and at least one of them is a fraction */
2010 scm_i_slow_exact_centered_remainder (SCM x
, SCM y
)
2014 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
2015 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG1
,
2016 s_scm_centered_remainder
);
2017 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
2018 SCM_WTA_DISPATCH_2 (g_scm_centered_remainder
, x
, y
, SCM_ARG2
,
2019 s_scm_centered_remainder
);
2020 else if (scm_is_true (scm_positive_p (y
)))
2021 q
= scm_floor (scm_sum (scm_divide (x
, y
), exactly_one_half
));
2022 else if (scm_is_true (scm_negative_p (y
)))
2023 q
= scm_ceiling (scm_difference (scm_divide (x
, y
), exactly_one_half
));
2025 scm_num_overflow (s_scm_centered_remainder
);
2026 return scm_difference (x
, scm_product (y
, q
));
2030 static SCM
scm_i_inexact_centered_quo_and_rem (double x
, double y
);
2031 static SCM
scm_i_bigint_centered_quo_and_rem (SCM x
, SCM y
);
2032 static SCM
scm_i_slow_exact_centered_quo_and_rem (SCM x
, SCM y
);
2034 SCM_PRIMITIVE_GENERIC (scm_centered_quo_and_rem
, "centered/", 2, 0, 0,
2036 "Return the integer @var{q} and the real number @var{r}\n"
2037 "such that @math{@var{x} = @var{q}*@var{y} + @var{r}}\n"
2038 "and @math{-abs(@var{y}/2) <= @var{r} < abs(@var{y}/2)}.\n"
2040 "(centered/ 123 10) @result{} 12 and 3\n"
2041 "(centered/ 123 -10) @result{} -12 and 3\n"
2042 "(centered/ -123 10) @result{} -12 and -3\n"
2043 "(centered/ -123 -10) @result{} 12 and -3\n"
2044 "(centered/ -123.2 -63.5) @result{} 2.0 and 3.8\n"
2045 "(centered/ 16/3 -10/7) @result{} -4 and -8/21\n"
2047 #define FUNC_NAME s_scm_centered_quo_and_rem
2049 if (SCM_LIKELY (SCM_I_INUMP (x
)))
2051 if (SCM_LIKELY (SCM_I_INUMP (y
)))
2053 scm_t_inum yy
= SCM_I_INUM (y
);
2054 if (SCM_UNLIKELY (yy
== 0))
2055 scm_num_overflow (s_scm_centered_quo_and_rem
);
2058 scm_t_inum xx
= SCM_I_INUM (x
);
2059 scm_t_inum qq
= xx
/ yy
;
2060 scm_t_inum rr
= xx
- qq
* yy
;
2061 if (SCM_LIKELY (xx
> 0))
2063 if (SCM_LIKELY (yy
> 0))
2065 if (rr
>= (yy
+ 1) / 2)
2070 if (rr
>= (1 - yy
) / 2)
2076 if (SCM_LIKELY (yy
> 0))
2087 return scm_values (scm_list_2 (SCM_I_MAKINUM (qq
),
2088 SCM_I_MAKINUM (rr
)));
2091 else if (SCM_BIGP (y
))
2093 /* Pass a denormalized bignum version of x (even though it
2094 can fit in a fixnum) to scm_i_bigint_centered_quo_and_rem */
2095 return scm_i_bigint_centered_quo_and_rem
2096 (scm_i_long2big (SCM_I_INUM (x
)), y
);
2098 else if (SCM_REALP (y
))
2099 return scm_i_inexact_centered_quo_and_rem
2100 (SCM_I_INUM (x
), SCM_REAL_VALUE (y
));
2101 else if (SCM_FRACTIONP (y
))
2102 return scm_i_slow_exact_centered_quo_and_rem (x
, y
);
2104 SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem
, x
, y
, SCM_ARG2
,
2105 s_scm_centered_quo_and_rem
);
2107 else if (SCM_BIGP (x
))
2109 if (SCM_LIKELY (SCM_I_INUMP (y
)))
2111 scm_t_inum yy
= SCM_I_INUM (y
);
2112 if (SCM_UNLIKELY (yy
== 0))
2113 scm_num_overflow (s_scm_centered_quo_and_rem
);
2116 SCM q
= scm_i_mkbig ();
2118 /* Arrange for rr to initially be non-positive,
2119 because that simplifies the test to see
2120 if it is within the needed bounds. */
2123 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
2124 SCM_I_BIG_MPZ (x
), yy
);
2125 scm_remember_upto_here_1 (x
);
2128 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
2129 SCM_I_BIG_MPZ (q
), 1);
2135 rr
= - mpz_cdiv_q_ui (SCM_I_BIG_MPZ (q
),
2136 SCM_I_BIG_MPZ (x
), -yy
);
2137 scm_remember_upto_here_1 (x
);
2138 mpz_neg (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (q
));
2141 mpz_add_ui (SCM_I_BIG_MPZ (q
),
2142 SCM_I_BIG_MPZ (q
), 1);
2146 return scm_values (scm_list_2 (scm_i_normbig (q
),
2147 SCM_I_MAKINUM (rr
)));
2150 else if (SCM_BIGP (y
))
2151 return scm_i_bigint_centered_quo_and_rem (x
, y
);
2152 else if (SCM_REALP (y
))
2153 return scm_i_inexact_centered_quo_and_rem
2154 (scm_i_big2dbl (x
), SCM_REAL_VALUE (y
));
2155 else if (SCM_FRACTIONP (y
))
2156 return scm_i_slow_exact_centered_quo_and_rem (x
, y
);
2158 SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem
, x
, y
, SCM_ARG2
,
2159 s_scm_centered_quo_and_rem
);
2161 else if (SCM_REALP (x
))
2163 if (SCM_REALP (y
) || SCM_I_INUMP (y
) ||
2164 SCM_BIGP (y
) || SCM_FRACTIONP (y
))
2165 return scm_i_inexact_centered_quo_and_rem
2166 (SCM_REAL_VALUE (x
), scm_to_double (y
));
2168 SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem
, x
, y
, SCM_ARG2
,
2169 s_scm_centered_quo_and_rem
);
2171 else if (SCM_FRACTIONP (x
))
2174 return scm_i_inexact_centered_quo_and_rem
2175 (scm_i_fraction2double (x
), SCM_REAL_VALUE (y
));
2177 return scm_i_slow_exact_centered_quo_and_rem (x
, y
);
2180 SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem
, x
, y
, SCM_ARG1
,
2181 s_scm_centered_quo_and_rem
);
2186 scm_i_inexact_centered_quo_and_rem (double x
, double y
)
2190 if (SCM_LIKELY (y
> 0))
2191 q
= floor (x
/y
+ 0.5);
2192 else if (SCM_LIKELY (y
< 0))
2193 q
= ceil (x
/y
- 0.5);
2195 scm_num_overflow (s_scm_centered_quo_and_rem
); /* or return a NaN? */
2199 return scm_values (scm_list_2 (scm_from_double (q
),
2200 scm_from_double (r
)));
2203 /* Assumes that both x and y are bigints, though
2204 x might be able to fit into a fixnum. */
2206 scm_i_bigint_centered_quo_and_rem (SCM x
, SCM y
)
2210 /* Note that x might be small enough to fit into a
2211 fixnum, so we must not let it escape into the wild */
2215 /* min_r will eventually become -abs(y/2) */
2216 min_r
= scm_i_mkbig ();
2217 mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (min_r
),
2218 SCM_I_BIG_MPZ (y
), 1);
2220 /* Arrange for rr to initially be non-positive,
2221 because that simplifies the test to see
2222 if it is within the needed bounds. */
2223 if (mpz_sgn (SCM_I_BIG_MPZ (y
)) > 0)
2225 mpz_cdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
2226 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2227 mpz_neg (SCM_I_BIG_MPZ (min_r
), SCM_I_BIG_MPZ (min_r
));
2228 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2230 mpz_sub_ui (SCM_I_BIG_MPZ (q
),
2231 SCM_I_BIG_MPZ (q
), 1);
2232 mpz_add (SCM_I_BIG_MPZ (r
),
2239 mpz_fdiv_qr (SCM_I_BIG_MPZ (q
), SCM_I_BIG_MPZ (r
),
2240 SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2241 if (mpz_cmp (SCM_I_BIG_MPZ (r
), SCM_I_BIG_MPZ (min_r
)) < 0)
2243 mpz_add_ui (SCM_I_BIG_MPZ (q
),
2244 SCM_I_BIG_MPZ (q
), 1);
2245 mpz_sub (SCM_I_BIG_MPZ (r
),
2250 scm_remember_upto_here_2 (x
, y
);
2251 return scm_values (scm_list_2 (scm_i_normbig (q
),
2252 scm_i_normbig (r
)));
2255 /* Compute exact centered quotient and remainder the slow way.
2256 We use this only if both arguments are exact,
2257 and at least one of them is a fraction */
2259 scm_i_slow_exact_centered_quo_and_rem (SCM x
, SCM y
)
2263 if (!(SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
)))
2264 SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem
, x
, y
, SCM_ARG1
,
2265 s_scm_centered_quo_and_rem
);
2266 else if (!(SCM_I_INUMP (y
) || SCM_BIGP (y
) || SCM_FRACTIONP (y
)))
2267 SCM_WTA_DISPATCH_2 (g_scm_centered_quo_and_rem
, x
, y
, SCM_ARG2
,
2268 s_scm_centered_quo_and_rem
);
2269 else if (scm_is_true (scm_positive_p (y
)))
2270 q
= scm_floor (scm_sum (scm_divide (x
, y
),
2272 else if (scm_is_true (scm_negative_p (y
)))
2273 q
= scm_ceiling (scm_difference (scm_divide (x
, y
),
2276 scm_num_overflow (s_scm_centered_quo_and_rem
);
2277 r
= scm_difference (x
, scm_product (q
, y
));
2278 return scm_values (scm_list_2 (q
, r
));
2282 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
2283 (SCM x
, SCM y
, SCM rest
),
2284 "Return the greatest common divisor of all parameter values.\n"
2285 "If called without arguments, 0 is returned.")
2286 #define FUNC_NAME s_scm_i_gcd
2288 while (!scm_is_null (rest
))
2289 { x
= scm_gcd (x
, y
);
2291 rest
= scm_cdr (rest
);
2293 return scm_gcd (x
, y
);
2297 #define s_gcd s_scm_i_gcd
2298 #define g_gcd g_scm_i_gcd
2301 scm_gcd (SCM x
, SCM y
)
2304 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
2306 if (SCM_I_INUMP (x
))
2308 if (SCM_I_INUMP (y
))
2310 scm_t_inum xx
= SCM_I_INUM (x
);
2311 scm_t_inum yy
= SCM_I_INUM (y
);
2312 scm_t_inum u
= xx
< 0 ? -xx
: xx
;
2313 scm_t_inum v
= yy
< 0 ? -yy
: yy
;
2323 /* Determine a common factor 2^k */
2324 while (!(1 & (u
| v
)))
2330 /* Now, any factor 2^n can be eliminated */
2350 return (SCM_POSFIXABLE (result
)
2351 ? SCM_I_MAKINUM (result
)
2352 : scm_i_inum2big (result
));
2354 else if (SCM_BIGP (y
))
2360 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
2362 else if (SCM_BIGP (x
))
2364 if (SCM_I_INUMP (y
))
2369 yy
= SCM_I_INUM (y
);
2374 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
2375 scm_remember_upto_here_1 (x
);
2376 return (SCM_POSFIXABLE (result
)
2377 ? SCM_I_MAKINUM (result
)
2378 : scm_from_unsigned_integer (result
));
2380 else if (SCM_BIGP (y
))
2382 SCM result
= scm_i_mkbig ();
2383 mpz_gcd (SCM_I_BIG_MPZ (result
),
2386 scm_remember_upto_here_2 (x
, y
);
2387 return scm_i_normbig (result
);
2390 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
2393 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
2396 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
2397 (SCM x
, SCM y
, SCM rest
),
2398 "Return the least common multiple of the arguments.\n"
2399 "If called without arguments, 1 is returned.")
2400 #define FUNC_NAME s_scm_i_lcm
2402 while (!scm_is_null (rest
))
2403 { x
= scm_lcm (x
, y
);
2405 rest
= scm_cdr (rest
);
2407 return scm_lcm (x
, y
);
2411 #define s_lcm s_scm_i_lcm
2412 #define g_lcm g_scm_i_lcm
2415 scm_lcm (SCM n1
, SCM n2
)
2417 if (SCM_UNBNDP (n2
))
2419 if (SCM_UNBNDP (n1
))
2420 return SCM_I_MAKINUM (1L);
2421 n2
= SCM_I_MAKINUM (1L);
2424 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
2425 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
2426 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
2427 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
2429 if (SCM_I_INUMP (n1
))
2431 if (SCM_I_INUMP (n2
))
2433 SCM d
= scm_gcd (n1
, n2
);
2434 if (scm_is_eq (d
, SCM_INUM0
))
2437 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
2441 /* inum n1, big n2 */
2444 SCM result
= scm_i_mkbig ();
2445 scm_t_inum nn1
= SCM_I_INUM (n1
);
2446 if (nn1
== 0) return SCM_INUM0
;
2447 if (nn1
< 0) nn1
= - nn1
;
2448 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
2449 scm_remember_upto_here_1 (n2
);
2457 if (SCM_I_INUMP (n2
))
2464 SCM result
= scm_i_mkbig ();
2465 mpz_lcm(SCM_I_BIG_MPZ (result
),
2467 SCM_I_BIG_MPZ (n2
));
2468 scm_remember_upto_here_2(n1
, n2
);
2469 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
2475 /* Emulating 2's complement bignums with sign magnitude arithmetic:
2480 + + + x (map digit:logand X Y)
2481 + - + x (map digit:logand X (lognot (+ -1 Y)))
2482 - + + y (map digit:logand (lognot (+ -1 X)) Y)
2483 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
2488 + + + (map digit:logior X Y)
2489 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
2490 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
2491 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
2496 + + + (map digit:logxor X Y)
2497 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
2498 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
2499 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
2504 + + (any digit:logand X Y)
2505 + - (any digit:logand X (lognot (+ -1 Y)))
2506 - + (any digit:logand (lognot (+ -1 X)) Y)
2511 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
2512 (SCM x
, SCM y
, SCM rest
),
2513 "Return the bitwise AND of the integer arguments.\n\n"
2515 "(logand) @result{} -1\n"
2516 "(logand 7) @result{} 7\n"
2517 "(logand #b111 #b011 #b001) @result{} 1\n"
2519 #define FUNC_NAME s_scm_i_logand
2521 while (!scm_is_null (rest
))
2522 { x
= scm_logand (x
, y
);
2524 rest
= scm_cdr (rest
);
2526 return scm_logand (x
, y
);
2530 #define s_scm_logand s_scm_i_logand
2532 SCM
scm_logand (SCM n1
, SCM n2
)
2533 #define FUNC_NAME s_scm_logand
2537 if (SCM_UNBNDP (n2
))
2539 if (SCM_UNBNDP (n1
))
2540 return SCM_I_MAKINUM (-1);
2541 else if (!SCM_NUMBERP (n1
))
2542 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2543 else if (SCM_NUMBERP (n1
))
2546 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2549 if (SCM_I_INUMP (n1
))
2551 nn1
= SCM_I_INUM (n1
);
2552 if (SCM_I_INUMP (n2
))
2554 scm_t_inum nn2
= SCM_I_INUM (n2
);
2555 return SCM_I_MAKINUM (nn1
& nn2
);
2557 else if SCM_BIGP (n2
)
2563 SCM result_z
= scm_i_mkbig ();
2565 mpz_init_set_si (nn1_z
, nn1
);
2566 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2567 scm_remember_upto_here_1 (n2
);
2569 return scm_i_normbig (result_z
);
2573 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2575 else if (SCM_BIGP (n1
))
2577 if (SCM_I_INUMP (n2
))
2580 nn1
= SCM_I_INUM (n1
);
2583 else if (SCM_BIGP (n2
))
2585 SCM result_z
= scm_i_mkbig ();
2586 mpz_and (SCM_I_BIG_MPZ (result_z
),
2588 SCM_I_BIG_MPZ (n2
));
2589 scm_remember_upto_here_2 (n1
, n2
);
2590 return scm_i_normbig (result_z
);
2593 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2596 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2601 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
2602 (SCM x
, SCM y
, SCM rest
),
2603 "Return the bitwise OR of the integer arguments.\n\n"
2605 "(logior) @result{} 0\n"
2606 "(logior 7) @result{} 7\n"
2607 "(logior #b000 #b001 #b011) @result{} 3\n"
2609 #define FUNC_NAME s_scm_i_logior
2611 while (!scm_is_null (rest
))
2612 { x
= scm_logior (x
, y
);
2614 rest
= scm_cdr (rest
);
2616 return scm_logior (x
, y
);
2620 #define s_scm_logior s_scm_i_logior
2622 SCM
scm_logior (SCM n1
, SCM n2
)
2623 #define FUNC_NAME s_scm_logior
2627 if (SCM_UNBNDP (n2
))
2629 if (SCM_UNBNDP (n1
))
2631 else if (SCM_NUMBERP (n1
))
2634 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2637 if (SCM_I_INUMP (n1
))
2639 nn1
= SCM_I_INUM (n1
);
2640 if (SCM_I_INUMP (n2
))
2642 long nn2
= SCM_I_INUM (n2
);
2643 return SCM_I_MAKINUM (nn1
| nn2
);
2645 else if (SCM_BIGP (n2
))
2651 SCM result_z
= scm_i_mkbig ();
2653 mpz_init_set_si (nn1_z
, nn1
);
2654 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2655 scm_remember_upto_here_1 (n2
);
2657 return scm_i_normbig (result_z
);
2661 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2663 else if (SCM_BIGP (n1
))
2665 if (SCM_I_INUMP (n2
))
2668 nn1
= SCM_I_INUM (n1
);
2671 else if (SCM_BIGP (n2
))
2673 SCM result_z
= scm_i_mkbig ();
2674 mpz_ior (SCM_I_BIG_MPZ (result_z
),
2676 SCM_I_BIG_MPZ (n2
));
2677 scm_remember_upto_here_2 (n1
, n2
);
2678 return scm_i_normbig (result_z
);
2681 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2684 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2689 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
2690 (SCM x
, SCM y
, SCM rest
),
2691 "Return the bitwise XOR of the integer arguments. A bit is\n"
2692 "set in the result if it is set in an odd number of arguments.\n"
2694 "(logxor) @result{} 0\n"
2695 "(logxor 7) @result{} 7\n"
2696 "(logxor #b000 #b001 #b011) @result{} 2\n"
2697 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
2699 #define FUNC_NAME s_scm_i_logxor
2701 while (!scm_is_null (rest
))
2702 { x
= scm_logxor (x
, y
);
2704 rest
= scm_cdr (rest
);
2706 return scm_logxor (x
, y
);
2710 #define s_scm_logxor s_scm_i_logxor
2712 SCM
scm_logxor (SCM n1
, SCM n2
)
2713 #define FUNC_NAME s_scm_logxor
2717 if (SCM_UNBNDP (n2
))
2719 if (SCM_UNBNDP (n1
))
2721 else if (SCM_NUMBERP (n1
))
2724 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2727 if (SCM_I_INUMP (n1
))
2729 nn1
= SCM_I_INUM (n1
);
2730 if (SCM_I_INUMP (n2
))
2732 scm_t_inum nn2
= SCM_I_INUM (n2
);
2733 return SCM_I_MAKINUM (nn1
^ nn2
);
2735 else if (SCM_BIGP (n2
))
2739 SCM result_z
= scm_i_mkbig ();
2741 mpz_init_set_si (nn1_z
, nn1
);
2742 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
2743 scm_remember_upto_here_1 (n2
);
2745 return scm_i_normbig (result_z
);
2749 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2751 else if (SCM_BIGP (n1
))
2753 if (SCM_I_INUMP (n2
))
2756 nn1
= SCM_I_INUM (n1
);
2759 else if (SCM_BIGP (n2
))
2761 SCM result_z
= scm_i_mkbig ();
2762 mpz_xor (SCM_I_BIG_MPZ (result_z
),
2764 SCM_I_BIG_MPZ (n2
));
2765 scm_remember_upto_here_2 (n1
, n2
);
2766 return scm_i_normbig (result_z
);
2769 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
2772 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
2777 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
2779 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
2780 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
2781 "without actually calculating the @code{logand}, just testing\n"
2785 "(logtest #b0100 #b1011) @result{} #f\n"
2786 "(logtest #b0100 #b0111) @result{} #t\n"
2788 #define FUNC_NAME s_scm_logtest
2792 if (SCM_I_INUMP (j
))
2794 nj
= SCM_I_INUM (j
);
2795 if (SCM_I_INUMP (k
))
2797 scm_t_inum nk
= SCM_I_INUM (k
);
2798 return scm_from_bool (nj
& nk
);
2800 else if (SCM_BIGP (k
))
2808 mpz_init_set_si (nj_z
, nj
);
2809 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
2810 scm_remember_upto_here_1 (k
);
2811 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
2817 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
2819 else if (SCM_BIGP (j
))
2821 if (SCM_I_INUMP (k
))
2824 nj
= SCM_I_INUM (j
);
2827 else if (SCM_BIGP (k
))
2831 mpz_init (result_z
);
2835 scm_remember_upto_here_2 (j
, k
);
2836 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
2837 mpz_clear (result_z
);
2841 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
2844 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
2849 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
2851 "Test whether bit number @var{index} in @var{j} is set.\n"
2852 "@var{index} starts from 0 for the least significant bit.\n"
2855 "(logbit? 0 #b1101) @result{} #t\n"
2856 "(logbit? 1 #b1101) @result{} #f\n"
2857 "(logbit? 2 #b1101) @result{} #t\n"
2858 "(logbit? 3 #b1101) @result{} #t\n"
2859 "(logbit? 4 #b1101) @result{} #f\n"
2861 #define FUNC_NAME s_scm_logbit_p
2863 unsigned long int iindex
;
2864 iindex
= scm_to_ulong (index
);
2866 if (SCM_I_INUMP (j
))
2868 /* bits above what's in an inum follow the sign bit */
2869 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
2870 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
2872 else if (SCM_BIGP (j
))
2874 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
2875 scm_remember_upto_here_1 (j
);
2876 return scm_from_bool (val
);
2879 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
2884 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
2886 "Return the integer which is the ones-complement of the integer\n"
2890 "(number->string (lognot #b10000000) 2)\n"
2891 " @result{} \"-10000001\"\n"
2892 "(number->string (lognot #b0) 2)\n"
2893 " @result{} \"-1\"\n"
2895 #define FUNC_NAME s_scm_lognot
2897 if (SCM_I_INUMP (n
)) {
2898 /* No overflow here, just need to toggle all the bits making up the inum.
2899 Enhancement: No need to strip the tag and add it back, could just xor
2900 a block of 1 bits, if that worked with the various debug versions of
2902 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
2904 } else if (SCM_BIGP (n
)) {
2905 SCM result
= scm_i_mkbig ();
2906 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
2907 scm_remember_upto_here_1 (n
);
2911 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2916 /* returns 0 if IN is not an integer. OUT must already be
2919 coerce_to_big (SCM in
, mpz_t out
)
2922 mpz_set (out
, SCM_I_BIG_MPZ (in
));
2923 else if (SCM_I_INUMP (in
))
2924 mpz_set_si (out
, SCM_I_INUM (in
));
2931 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
2932 (SCM n
, SCM k
, SCM m
),
2933 "Return @var{n} raised to the integer exponent\n"
2934 "@var{k}, modulo @var{m}.\n"
2937 "(modulo-expt 2 3 5)\n"
2940 #define FUNC_NAME s_scm_modulo_expt
2946 /* There are two classes of error we might encounter --
2947 1) Math errors, which we'll report by calling scm_num_overflow,
2949 2) wrong-type errors, which of course we'll report by calling
2951 We don't report those errors immediately, however; instead we do
2952 some cleanup first. These variables tell us which error (if
2953 any) we should report after cleaning up.
2955 int report_overflow
= 0;
2957 int position_of_wrong_type
= 0;
2958 SCM value_of_wrong_type
= SCM_INUM0
;
2960 SCM result
= SCM_UNDEFINED
;
2966 if (scm_is_eq (m
, SCM_INUM0
))
2968 report_overflow
= 1;
2972 if (!coerce_to_big (n
, n_tmp
))
2974 value_of_wrong_type
= n
;
2975 position_of_wrong_type
= 1;
2979 if (!coerce_to_big (k
, k_tmp
))
2981 value_of_wrong_type
= k
;
2982 position_of_wrong_type
= 2;
2986 if (!coerce_to_big (m
, m_tmp
))
2988 value_of_wrong_type
= m
;
2989 position_of_wrong_type
= 3;
2993 /* if the exponent K is negative, and we simply call mpz_powm, we
2994 will get a divide-by-zero exception when an inverse 1/n mod m
2995 doesn't exist (or is not unique). Since exceptions are hard to
2996 handle, we'll attempt the inversion "by hand" -- that way, we get
2997 a simple failure code, which is easy to handle. */
2999 if (-1 == mpz_sgn (k_tmp
))
3001 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
3003 report_overflow
= 1;
3006 mpz_neg (k_tmp
, k_tmp
);
3009 result
= scm_i_mkbig ();
3010 mpz_powm (SCM_I_BIG_MPZ (result
),
3015 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
3016 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
3023 if (report_overflow
)
3024 scm_num_overflow (FUNC_NAME
);
3026 if (position_of_wrong_type
)
3027 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
3028 value_of_wrong_type
);
3030 return scm_i_normbig (result
);
3034 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
3036 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
3037 "exact integer, @var{n} can be any number.\n"
3039 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
3040 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
3041 "includes @math{0^0} is 1.\n"
3044 "(integer-expt 2 5) @result{} 32\n"
3045 "(integer-expt -3 3) @result{} -27\n"
3046 "(integer-expt 5 -3) @result{} 1/125\n"
3047 "(integer-expt 0 0) @result{} 1\n"
3049 #define FUNC_NAME s_scm_integer_expt
3052 SCM z_i2
= SCM_BOOL_F
;
3054 SCM acc
= SCM_I_MAKINUM (1L);
3056 SCM_VALIDATE_NUMBER (SCM_ARG1
, n
);
3057 if (!SCM_I_INUMP (k
) && !SCM_BIGP (k
))
3058 SCM_WRONG_TYPE_ARG (2, k
);
3060 if (scm_is_true (scm_zero_p (n
)))
3062 if (scm_is_true (scm_zero_p (k
))) /* 0^0 == 1 per R5RS */
3063 return acc
; /* return exact 1, regardless of n */
3064 else if (scm_is_true (scm_positive_p (k
)))
3066 else /* return NaN for (0 ^ k) for negative k per R6RS */
3069 else if (scm_is_eq (n
, acc
))
3071 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
3072 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
3074 if (SCM_I_INUMP (k
))
3075 i2
= SCM_I_INUM (k
);
3076 else if (SCM_BIGP (k
))
3078 z_i2
= scm_i_clonebig (k
, 1);
3079 scm_remember_upto_here_1 (k
);
3083 SCM_WRONG_TYPE_ARG (2, k
);
3087 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
3089 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
3090 n
= scm_divide (n
, SCM_UNDEFINED
);
3094 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
3098 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
3100 return scm_product (acc
, n
);
3102 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
3103 acc
= scm_product (acc
, n
);
3104 n
= scm_product (n
, n
);
3105 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
3113 n
= scm_divide (n
, SCM_UNDEFINED
);
3120 return scm_product (acc
, n
);
3122 acc
= scm_product (acc
, n
);
3123 n
= scm_product (n
, n
);
3130 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
3132 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
3133 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
3135 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
3136 "@var{cnt} is negative it's a division, rounded towards negative\n"
3137 "infinity. (Note that this is not the same rounding as\n"
3138 "@code{quotient} does.)\n"
3140 "With @var{n} viewed as an infinite precision twos complement,\n"
3141 "@code{ash} means a left shift introducing zero bits, or a right\n"
3142 "shift dropping bits.\n"
3145 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
3146 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
3148 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
3149 "(ash -23 -2) @result{} -6\n"
3151 #define FUNC_NAME s_scm_ash
3154 bits_to_shift
= scm_to_long (cnt
);
3156 if (SCM_I_INUMP (n
))
3158 scm_t_inum nn
= SCM_I_INUM (n
);
3160 if (bits_to_shift
> 0)
3162 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
3163 overflow a non-zero fixnum. For smaller shifts we check the
3164 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
3165 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
3166 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
3172 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
3174 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
3177 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
3181 SCM result
= scm_i_inum2big (nn
);
3182 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
3189 bits_to_shift
= -bits_to_shift
;
3190 if (bits_to_shift
>= SCM_LONG_BIT
)
3191 return (nn
>= 0 ? SCM_INUM0
: SCM_I_MAKINUM(-1));
3193 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
3197 else if (SCM_BIGP (n
))
3201 if (bits_to_shift
== 0)
3204 result
= scm_i_mkbig ();
3205 if (bits_to_shift
>= 0)
3207 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
3213 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
3214 we have to allocate a bignum even if the result is going to be a
3216 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
3218 return scm_i_normbig (result
);
3224 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3230 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
3231 (SCM n
, SCM start
, SCM end
),
3232 "Return the integer composed of the @var{start} (inclusive)\n"
3233 "through @var{end} (exclusive) bits of @var{n}. The\n"
3234 "@var{start}th bit becomes the 0-th bit in the result.\n"
3237 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
3238 " @result{} \"1010\"\n"
3239 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
3240 " @result{} \"10110\"\n"
3242 #define FUNC_NAME s_scm_bit_extract
3244 unsigned long int istart
, iend
, bits
;
3245 istart
= scm_to_ulong (start
);
3246 iend
= scm_to_ulong (end
);
3247 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
3249 /* how many bits to keep */
3250 bits
= iend
- istart
;
3252 if (SCM_I_INUMP (n
))
3254 scm_t_inum in
= SCM_I_INUM (n
);
3256 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
3257 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
3258 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
3260 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
3262 /* Since we emulate two's complement encoded numbers, this
3263 * special case requires us to produce a result that has
3264 * more bits than can be stored in a fixnum.
3266 SCM result
= scm_i_inum2big (in
);
3267 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
3272 /* mask down to requisite bits */
3273 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
3274 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
3276 else if (SCM_BIGP (n
))
3281 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
3285 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
3286 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
3287 such bits into a ulong. */
3288 result
= scm_i_mkbig ();
3289 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
3290 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
3291 result
= scm_i_normbig (result
);
3293 scm_remember_upto_here_1 (n
);
3297 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3302 static const char scm_logtab
[] = {
3303 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
3306 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
3308 "Return the number of bits in integer @var{n}. If integer is\n"
3309 "positive, the 1-bits in its binary representation are counted.\n"
3310 "If negative, the 0-bits in its two's-complement binary\n"
3311 "representation are counted. If 0, 0 is returned.\n"
3314 "(logcount #b10101010)\n"
3321 #define FUNC_NAME s_scm_logcount
3323 if (SCM_I_INUMP (n
))
3325 unsigned long c
= 0;
3326 scm_t_inum nn
= SCM_I_INUM (n
);
3331 c
+= scm_logtab
[15 & nn
];
3334 return SCM_I_MAKINUM (c
);
3336 else if (SCM_BIGP (n
))
3338 unsigned long count
;
3339 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
3340 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
3342 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
3343 scm_remember_upto_here_1 (n
);
3344 return SCM_I_MAKINUM (count
);
3347 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3352 static const char scm_ilentab
[] = {
3353 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
3357 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
3359 "Return the number of bits necessary to represent @var{n}.\n"
3362 "(integer-length #b10101010)\n"
3364 "(integer-length 0)\n"
3366 "(integer-length #b1111)\n"
3369 #define FUNC_NAME s_scm_integer_length
3371 if (SCM_I_INUMP (n
))
3373 unsigned long c
= 0;
3375 scm_t_inum nn
= SCM_I_INUM (n
);
3381 l
= scm_ilentab
[15 & nn
];
3384 return SCM_I_MAKINUM (c
- 4 + l
);
3386 else if (SCM_BIGP (n
))
3388 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
3389 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
3390 1 too big, so check for that and adjust. */
3391 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
3392 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
3393 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
3394 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
3396 scm_remember_upto_here_1 (n
);
3397 return SCM_I_MAKINUM (size
);
3400 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
3404 /*** NUMBERS -> STRINGS ***/
3405 #define SCM_MAX_DBL_PREC 60
3406 #define SCM_MAX_DBL_RADIX 36
3408 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
3409 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
3410 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
3413 void init_dblprec(int *prec
, int radix
) {
3414 /* determine floating point precision by adding successively
3415 smaller increments to 1.0 until it is considered == 1.0 */
3416 double f
= ((double)1.0)/radix
;
3417 double fsum
= 1.0 + f
;
3422 if (++(*prec
) > SCM_MAX_DBL_PREC
)
3434 void init_fx_radix(double *fx_list
, int radix
)
3436 /* initialize a per-radix list of tolerances. When added
3437 to a number < 1.0, we can determine if we should raund
3438 up and quit converting a number to a string. */
3442 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
3443 fx_list
[i
] = (fx_list
[i
-1] / radix
);
3446 /* use this array as a way to generate a single digit */
3447 static const char number_chars
[] = "0123456789abcdefghijklmnopqrstuvwxyz";
3450 idbl2str (double f
, char *a
, int radix
)
3452 int efmt
, dpt
, d
, i
, wp
;
3454 #ifdef DBL_MIN_10_EXP
3457 #endif /* DBL_MIN_10_EXP */
3462 radix
> SCM_MAX_DBL_RADIX
)
3464 /* revert to existing behavior */
3468 wp
= scm_dblprec
[radix
-2];
3469 fx
= fx_per_radix
[radix
-2];
3473 #ifdef HAVE_COPYSIGN
3474 double sgn
= copysign (1.0, f
);
3479 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
3485 strcpy (a
, "-inf.0");
3487 strcpy (a
, "+inf.0");
3492 strcpy (a
, "+nan.0");
3502 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
3503 make-uniform-vector, from causing infinite loops. */
3504 /* just do the checking...if it passes, we do the conversion for our
3505 radix again below */
3512 if (exp_cpy
-- < DBL_MIN_10_EXP
)
3520 while (f_cpy
> 10.0)
3523 if (exp_cpy
++ > DBL_MAX_10_EXP
)
3544 if (f
+ fx
[wp
] >= radix
)
3551 /* adding 9999 makes this equivalent to abs(x) % 3 */
3552 dpt
= (exp
+ 9999) % 3;
3556 efmt
= (exp
< -3) || (exp
> wp
+ 2);
3578 a
[ch
++] = number_chars
[d
];
3581 if (f
+ fx
[wp
] >= 1.0)
3583 a
[ch
- 1] = number_chars
[d
+1];
3595 if ((dpt
> 4) && (exp
> 6))
3597 d
= (a
[0] == '-' ? 2 : 1);
3598 for (i
= ch
++; i
> d
; i
--)
3611 if (a
[ch
- 1] == '.')
3612 a
[ch
++] = '0'; /* trailing zero */
3621 for (i
= radix
; i
<= exp
; i
*= radix
);
3622 for (i
/= radix
; i
; i
/= radix
)
3624 a
[ch
++] = number_chars
[exp
/ i
];
3633 icmplx2str (double real
, double imag
, char *str
, int radix
)
3637 i
= idbl2str (real
, str
, radix
);
3640 /* Don't output a '+' for negative numbers or for Inf and
3641 NaN. They will provide their own sign. */
3642 if (0 <= imag
&& !isinf (imag
) && !isnan (imag
))
3644 i
+= idbl2str (imag
, &str
[i
], radix
);
3651 iflo2str (SCM flt
, char *str
, int radix
)
3654 if (SCM_REALP (flt
))
3655 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
3657 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
3662 /* convert a scm_t_intmax to a string (unterminated). returns the number of
3663 characters in the result.
3665 p is destination: worst case (base 2) is SCM_INTBUFLEN */
3667 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
3672 return scm_iuint2str (-num
, rad
, p
) + 1;
3675 return scm_iuint2str (num
, rad
, p
);
3678 /* convert a scm_t_intmax to a string (unterminated). returns the number of
3679 characters in the result.
3681 p is destination: worst case (base 2) is SCM_INTBUFLEN */
3683 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
3687 scm_t_uintmax n
= num
;
3689 if (rad
< 2 || rad
> 36)
3690 scm_out_of_range ("scm_iuint2str", scm_from_int (rad
));
3692 for (n
/= rad
; n
> 0; n
/= rad
)
3702 p
[i
] = number_chars
[d
];
3707 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
3709 "Return a string holding the external representation of the\n"
3710 "number @var{n} in the given @var{radix}. If @var{n} is\n"
3711 "inexact, a radix of 10 will be used.")
3712 #define FUNC_NAME s_scm_number_to_string
3716 if (SCM_UNBNDP (radix
))
3719 base
= scm_to_signed_integer (radix
, 2, 36);
3721 if (SCM_I_INUMP (n
))
3723 char num_buf
[SCM_INTBUFLEN
];
3724 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
3725 return scm_from_locale_stringn (num_buf
, length
);
3727 else if (SCM_BIGP (n
))
3729 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
3730 scm_remember_upto_here_1 (n
);
3731 return scm_take_locale_string (str
);
3733 else if (SCM_FRACTIONP (n
))
3735 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
3736 scm_from_locale_string ("/"),
3737 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
3739 else if (SCM_INEXACTP (n
))
3741 char num_buf
[FLOBUFLEN
];
3742 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
3745 SCM_WRONG_TYPE_ARG (1, n
);
3750 /* These print routines used to be stubbed here so that scm_repl.c
3751 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
3754 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3756 char num_buf
[FLOBUFLEN
];
3757 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
3762 scm_i_print_double (double val
, SCM port
)
3764 char num_buf
[FLOBUFLEN
];
3765 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
3769 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3772 char num_buf
[FLOBUFLEN
];
3773 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
3778 scm_i_print_complex (double real
, double imag
, SCM port
)
3780 char num_buf
[FLOBUFLEN
];
3781 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
3785 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3788 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
3789 scm_display (str
, port
);
3790 scm_remember_upto_here_1 (str
);
3795 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
3797 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
3798 scm_remember_upto_here_1 (exp
);
3799 scm_lfwrite (str
, (size_t) strlen (str
), port
);
3803 /*** END nums->strs ***/
3806 /*** STRINGS -> NUMBERS ***/
3808 /* The following functions implement the conversion from strings to numbers.
3809 * The implementation somehow follows the grammar for numbers as it is given
3810 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
3811 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
3812 * points should be noted about the implementation:
3813 * * Each function keeps a local index variable 'idx' that points at the
3814 * current position within the parsed string. The global index is only
3815 * updated if the function could parse the corresponding syntactic unit
3817 * * Similarly, the functions keep track of indicators of inexactness ('#',
3818 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
3819 * global exactness information is only updated after each part has been
3820 * successfully parsed.
3821 * * Sequences of digits are parsed into temporary variables holding fixnums.
3822 * Only if these fixnums would overflow, the result variables are updated
3823 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
3824 * the temporary variables holding the fixnums are cleared, and the process
3825 * starts over again. If for example fixnums were able to store five decimal
3826 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
3827 * and the result was computed as 12345 * 100000 + 67890. In other words,
3828 * only every five digits two bignum operations were performed.
3831 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
3833 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
3835 /* Caller is responsible for checking that the return value is in range
3836 for the given radix, which should be <= 36. */
3838 char_decimal_value (scm_t_uint32 c
)
3840 /* uc_decimal_value returns -1 on error. When cast to an unsigned int,
3841 that's certainly above any valid decimal, so we take advantage of
3842 that to elide some tests. */
3843 unsigned int d
= (unsigned int) uc_decimal_value (c
);
3845 /* If that failed, try extended hexadecimals, then. Only accept ascii
3850 if (c
>= (scm_t_uint32
) 'a')
3851 d
= c
- (scm_t_uint32
)'a' + 10U;
3857 mem2uinteger (SCM mem
, unsigned int *p_idx
,
3858 unsigned int radix
, enum t_exactness
*p_exactness
)
3860 unsigned int idx
= *p_idx
;
3861 unsigned int hash_seen
= 0;
3862 scm_t_bits shift
= 1;
3864 unsigned int digit_value
;
3867 size_t len
= scm_i_string_length (mem
);
3872 c
= scm_i_string_ref (mem
, idx
);
3873 digit_value
= char_decimal_value (c
);
3874 if (digit_value
>= radix
)
3878 result
= SCM_I_MAKINUM (digit_value
);
3881 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
3891 digit_value
= char_decimal_value (c
);
3892 /* This check catches non-decimals in addition to out-of-range
3894 if (digit_value
>= radix
)
3899 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
3901 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3903 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3910 shift
= shift
* radix
;
3911 add
= add
* radix
+ digit_value
;
3916 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3918 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3922 *p_exactness
= INEXACT
;
3928 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
3929 * covers the parts of the rules that start at a potential point. The value
3930 * of the digits up to the point have been parsed by the caller and are given
3931 * in variable result. The content of *p_exactness indicates, whether a hash
3932 * has already been seen in the digits before the point.
3935 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
3938 mem2decimal_from_point (SCM result
, SCM mem
,
3939 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
3941 unsigned int idx
= *p_idx
;
3942 enum t_exactness x
= *p_exactness
;
3943 size_t len
= scm_i_string_length (mem
);
3948 if (scm_i_string_ref (mem
, idx
) == '.')
3950 scm_t_bits shift
= 1;
3952 unsigned int digit_value
;
3953 SCM big_shift
= SCM_INUM1
;
3958 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
3959 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
3964 digit_value
= DIGIT2UINT (c
);
3975 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
3977 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
3978 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3980 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3988 add
= add
* 10 + digit_value
;
3994 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
3995 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
3996 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
3999 result
= scm_divide (result
, big_shift
);
4001 /* We've seen a decimal point, thus the value is implicitly inexact. */
4013 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
4015 switch (scm_i_string_ref (mem
, idx
))
4027 c
= scm_i_string_ref (mem
, idx
);
4035 c
= scm_i_string_ref (mem
, idx
);
4044 c
= scm_i_string_ref (mem
, idx
);
4049 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4053 exponent
= DIGIT2UINT (c
);
4056 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
4057 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
4060 if (exponent
<= SCM_MAXEXP
)
4061 exponent
= exponent
* 10 + DIGIT2UINT (c
);
4067 if (exponent
> SCM_MAXEXP
)
4069 size_t exp_len
= idx
- start
;
4070 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
4071 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
4072 scm_out_of_range ("string->number", exp_num
);
4075 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
4077 result
= scm_product (result
, e
);
4079 result
= scm_divide2real (result
, e
);
4081 /* We've seen an exponent, thus the value is implicitly inexact. */
4099 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
4102 mem2ureal (SCM mem
, unsigned int *p_idx
,
4103 unsigned int radix
, enum t_exactness
*p_exactness
)
4105 unsigned int idx
= *p_idx
;
4107 size_t len
= scm_i_string_length (mem
);
4109 /* Start off believing that the number will be exact. This changes
4110 to INEXACT if we see a decimal point or a hash. */
4111 enum t_exactness x
= EXACT
;
4116 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
4122 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
4124 /* Cobble up the fractional part. We might want to set the
4125 NaN's mantissa from it. */
4127 mem2uinteger (mem
, &idx
, 10, &x
);
4132 if (scm_i_string_ref (mem
, idx
) == '.')
4136 else if (idx
+ 1 == len
)
4138 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
4141 result
= mem2decimal_from_point (SCM_INUM0
, mem
,
4148 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
4149 if (scm_is_false (uinteger
))
4154 else if (scm_i_string_ref (mem
, idx
) == '/')
4162 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
4163 if (scm_is_false (divisor
))
4166 /* both are int/big here, I assume */
4167 result
= scm_i_make_ratio (uinteger
, divisor
);
4169 else if (radix
== 10)
4171 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
4172 if (scm_is_false (result
))
4181 /* Update *p_exactness if the number just read was inexact. This is
4182 important for complex numbers, so that a complex number is
4183 treated as inexact overall if either its real or imaginary part
4189 /* When returning an inexact zero, make sure it is represented as a
4190 floating point value so that we can change its sign.
4192 if (scm_is_eq (result
, SCM_INUM0
) && *p_exactness
== INEXACT
)
4193 result
= scm_from_double (0.0);
4199 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
4202 mem2complex (SCM mem
, unsigned int idx
,
4203 unsigned int radix
, enum t_exactness
*p_exactness
)
4208 size_t len
= scm_i_string_length (mem
);
4213 c
= scm_i_string_ref (mem
, idx
);
4228 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
4229 if (scm_is_false (ureal
))
4231 /* input must be either +i or -i */
4236 if (scm_i_string_ref (mem
, idx
) == 'i'
4237 || scm_i_string_ref (mem
, idx
) == 'I')
4243 return scm_make_rectangular (SCM_INUM0
, SCM_I_MAKINUM (sign
));
4250 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
4251 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
4256 c
= scm_i_string_ref (mem
, idx
);
4260 /* either +<ureal>i or -<ureal>i */
4267 return scm_make_rectangular (SCM_INUM0
, ureal
);
4270 /* polar input: <real>@<real>. */
4281 c
= scm_i_string_ref (mem
, idx
);
4299 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
4300 if (scm_is_false (angle
))
4305 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
4306 angle
= scm_difference (angle
, SCM_UNDEFINED
);
4308 result
= scm_make_polar (ureal
, angle
);
4313 /* expecting input matching <real>[+-]<ureal>?i */
4320 int sign
= (c
== '+') ? 1 : -1;
4321 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
4323 if (scm_is_false (imag
))
4324 imag
= SCM_I_MAKINUM (sign
);
4325 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
4326 imag
= scm_difference (imag
, SCM_UNDEFINED
);
4330 if (scm_i_string_ref (mem
, idx
) != 'i'
4331 && scm_i_string_ref (mem
, idx
) != 'I')
4338 return scm_make_rectangular (ureal
, imag
);
4347 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
4349 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
4352 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
4354 unsigned int idx
= 0;
4355 unsigned int radix
= NO_RADIX
;
4356 enum t_exactness forced_x
= NO_EXACTNESS
;
4357 enum t_exactness implicit_x
= EXACT
;
4359 size_t len
= scm_i_string_length (mem
);
4361 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
4362 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
4364 switch (scm_i_string_ref (mem
, idx
+ 1))
4367 if (radix
!= NO_RADIX
)
4372 if (radix
!= NO_RADIX
)
4377 if (forced_x
!= NO_EXACTNESS
)
4382 if (forced_x
!= NO_EXACTNESS
)
4387 if (radix
!= NO_RADIX
)
4392 if (radix
!= NO_RADIX
)
4402 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
4403 if (radix
== NO_RADIX
)
4404 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
4406 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
4408 if (scm_is_false (result
))
4414 if (SCM_INEXACTP (result
))
4415 return scm_inexact_to_exact (result
);
4419 if (SCM_INEXACTP (result
))
4422 return scm_exact_to_inexact (result
);
4425 if (implicit_x
== INEXACT
)
4427 if (SCM_INEXACTP (result
))
4430 return scm_exact_to_inexact (result
);
4438 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
4439 unsigned int default_radix
)
4441 SCM str
= scm_from_locale_stringn (mem
, len
);
4443 return scm_i_string_to_number (str
, default_radix
);
4447 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
4448 (SCM string
, SCM radix
),
4449 "Return a number of the maximally precise representation\n"
4450 "expressed by the given @var{string}. @var{radix} must be an\n"
4451 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
4452 "is a default radix that may be overridden by an explicit radix\n"
4453 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
4454 "supplied, then the default radix is 10. If string is not a\n"
4455 "syntactically valid notation for a number, then\n"
4456 "@code{string->number} returns @code{#f}.")
4457 #define FUNC_NAME s_scm_string_to_number
4461 SCM_VALIDATE_STRING (1, string
);
4463 if (SCM_UNBNDP (radix
))
4466 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
4468 answer
= scm_i_string_to_number (string
, base
);
4469 scm_remember_upto_here_1 (string
);
4475 /*** END strs->nums ***/
4478 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
4480 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
4482 #define FUNC_NAME s_scm_number_p
4484 return scm_from_bool (SCM_NUMBERP (x
));
4488 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
4490 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
4491 "otherwise. Note that the sets of real, rational and integer\n"
4492 "values form subsets of the set of complex numbers, i. e. the\n"
4493 "predicate will also be fulfilled if @var{x} is a real,\n"
4494 "rational or integer number.")
4495 #define FUNC_NAME s_scm_complex_p
4497 /* all numbers are complex. */
4498 return scm_number_p (x
);
4502 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
4504 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
4505 "otherwise. Note that the set of integer values forms a subset of\n"
4506 "the set of real numbers, i. e. the predicate will also be\n"
4507 "fulfilled if @var{x} is an integer number.")
4508 #define FUNC_NAME s_scm_real_p
4510 return scm_from_bool
4511 (SCM_I_INUMP (x
) || SCM_REALP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
));
4515 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
4517 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
4518 "otherwise. Note that the set of integer values forms a subset of\n"
4519 "the set of rational numbers, i. e. the predicate will also be\n"
4520 "fulfilled if @var{x} is an integer number.")
4521 #define FUNC_NAME s_scm_rational_p
4523 if (SCM_I_INUMP (x
) || SCM_BIGP (x
) || SCM_FRACTIONP (x
))
4525 else if (SCM_REALP (x
))
4526 /* due to their limited precision, finite floating point numbers are
4527 rational as well. (finite means neither infinity nor a NaN) */
4528 return scm_from_bool (DOUBLE_IS_FINITE (SCM_REAL_VALUE (x
)));
4534 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
4536 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
4538 #define FUNC_NAME s_scm_integer_p
4540 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4542 else if (SCM_REALP (x
))
4544 double val
= SCM_REAL_VALUE (x
);
4545 return scm_from_bool (!isinf (val
) && (val
== floor (val
)));
4553 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
4554 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
4555 (SCM x
, SCM y
, SCM rest
),
4556 "Return @code{#t} if all parameters are numerically equal.")
4557 #define FUNC_NAME s_scm_i_num_eq_p
4559 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4561 while (!scm_is_null (rest
))
4563 if (scm_is_false (scm_num_eq_p (x
, y
)))
4567 rest
= scm_cdr (rest
);
4569 return scm_num_eq_p (x
, y
);
4573 scm_num_eq_p (SCM x
, SCM y
)
4576 if (SCM_I_INUMP (x
))
4578 scm_t_signed_bits xx
= SCM_I_INUM (x
);
4579 if (SCM_I_INUMP (y
))
4581 scm_t_signed_bits yy
= SCM_I_INUM (y
);
4582 return scm_from_bool (xx
== yy
);
4584 else if (SCM_BIGP (y
))
4586 else if (SCM_REALP (y
))
4588 /* On a 32-bit system an inum fits a double, we can cast the inum
4589 to a double and compare.
4591 But on a 64-bit system an inum is bigger than a double and
4592 casting it to a double (call that dxx) will round. dxx is at
4593 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
4594 an integer and fits a long. So we cast yy to a long and
4595 compare with plain xx.
4597 An alternative (for any size system actually) would be to check
4598 yy is an integer (with floor) and is in range of an inum
4599 (compare against appropriate powers of 2) then test
4600 xx==(scm_t_signed_bits)yy. It's just a matter of which
4601 casts/comparisons might be fastest or easiest for the cpu. */
4603 double yy
= SCM_REAL_VALUE (y
);
4604 return scm_from_bool ((double) xx
== yy
4605 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
4606 || xx
== (scm_t_signed_bits
) yy
));
4608 else if (SCM_COMPLEXP (y
))
4609 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
4610 && (0.0 == SCM_COMPLEX_IMAG (y
)));
4611 else if (SCM_FRACTIONP (y
))
4614 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4616 else if (SCM_BIGP (x
))
4618 if (SCM_I_INUMP (y
))
4620 else if (SCM_BIGP (y
))
4622 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4623 scm_remember_upto_here_2 (x
, y
);
4624 return scm_from_bool (0 == cmp
);
4626 else if (SCM_REALP (y
))
4629 if (isnan (SCM_REAL_VALUE (y
)))
4631 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
4632 scm_remember_upto_here_1 (x
);
4633 return scm_from_bool (0 == cmp
);
4635 else if (SCM_COMPLEXP (y
))
4638 if (0.0 != SCM_COMPLEX_IMAG (y
))
4640 if (isnan (SCM_COMPLEX_REAL (y
)))
4642 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
4643 scm_remember_upto_here_1 (x
);
4644 return scm_from_bool (0 == cmp
);
4646 else if (SCM_FRACTIONP (y
))
4649 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4651 else if (SCM_REALP (x
))
4653 double xx
= SCM_REAL_VALUE (x
);
4654 if (SCM_I_INUMP (y
))
4656 /* see comments with inum/real above */
4657 scm_t_signed_bits yy
= SCM_I_INUM (y
);
4658 return scm_from_bool (xx
== (double) yy
4659 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
4660 || (scm_t_signed_bits
) xx
== yy
));
4662 else if (SCM_BIGP (y
))
4665 if (isnan (SCM_REAL_VALUE (x
)))
4667 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
4668 scm_remember_upto_here_1 (y
);
4669 return scm_from_bool (0 == cmp
);
4671 else if (SCM_REALP (y
))
4672 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
4673 else if (SCM_COMPLEXP (y
))
4674 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
4675 && (0.0 == SCM_COMPLEX_IMAG (y
)));
4676 else if (SCM_FRACTIONP (y
))
4678 double xx
= SCM_REAL_VALUE (x
);
4682 return scm_from_bool (xx
< 0.0);
4683 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4687 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4689 else if (SCM_COMPLEXP (x
))
4691 if (SCM_I_INUMP (y
))
4692 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
4693 && (SCM_COMPLEX_IMAG (x
) == 0.0));
4694 else if (SCM_BIGP (y
))
4697 if (0.0 != SCM_COMPLEX_IMAG (x
))
4699 if (isnan (SCM_COMPLEX_REAL (x
)))
4701 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
4702 scm_remember_upto_here_1 (y
);
4703 return scm_from_bool (0 == cmp
);
4705 else if (SCM_REALP (y
))
4706 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
4707 && (SCM_COMPLEX_IMAG (x
) == 0.0));
4708 else if (SCM_COMPLEXP (y
))
4709 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
4710 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
4711 else if (SCM_FRACTIONP (y
))
4714 if (SCM_COMPLEX_IMAG (x
) != 0.0)
4716 xx
= SCM_COMPLEX_REAL (x
);
4720 return scm_from_bool (xx
< 0.0);
4721 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4725 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4727 else if (SCM_FRACTIONP (x
))
4729 if (SCM_I_INUMP (y
))
4731 else if (SCM_BIGP (y
))
4733 else if (SCM_REALP (y
))
4735 double yy
= SCM_REAL_VALUE (y
);
4739 return scm_from_bool (0.0 < yy
);
4740 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4743 else if (SCM_COMPLEXP (y
))
4746 if (SCM_COMPLEX_IMAG (y
) != 0.0)
4748 yy
= SCM_COMPLEX_REAL (y
);
4752 return scm_from_bool (0.0 < yy
);
4753 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4756 else if (SCM_FRACTIONP (y
))
4757 return scm_i_fraction_equalp (x
, y
);
4759 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
4762 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
4766 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
4767 done are good for inums, but for bignums an answer can almost always be
4768 had by just examining a few high bits of the operands, as done by GMP in
4769 mpq_cmp. flonum/frac compares likewise, but with the slight complication
4770 of the float exponent to take into account. */
4772 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
4773 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
4774 (SCM x
, SCM y
, SCM rest
),
4775 "Return @code{#t} if the list of parameters is monotonically\n"
4777 #define FUNC_NAME s_scm_i_num_less_p
4779 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4781 while (!scm_is_null (rest
))
4783 if (scm_is_false (scm_less_p (x
, y
)))
4787 rest
= scm_cdr (rest
);
4789 return scm_less_p (x
, y
);
4793 scm_less_p (SCM x
, SCM y
)
4796 if (SCM_I_INUMP (x
))
4798 scm_t_inum xx
= SCM_I_INUM (x
);
4799 if (SCM_I_INUMP (y
))
4801 scm_t_inum yy
= SCM_I_INUM (y
);
4802 return scm_from_bool (xx
< yy
);
4804 else if (SCM_BIGP (y
))
4806 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4807 scm_remember_upto_here_1 (y
);
4808 return scm_from_bool (sgn
> 0);
4810 else if (SCM_REALP (y
))
4811 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
4812 else if (SCM_FRACTIONP (y
))
4814 /* "x < a/b" becomes "x*b < a" */
4816 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
4817 y
= SCM_FRACTION_NUMERATOR (y
);
4821 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4823 else if (SCM_BIGP (x
))
4825 if (SCM_I_INUMP (y
))
4827 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4828 scm_remember_upto_here_1 (x
);
4829 return scm_from_bool (sgn
< 0);
4831 else if (SCM_BIGP (y
))
4833 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4834 scm_remember_upto_here_2 (x
, y
);
4835 return scm_from_bool (cmp
< 0);
4837 else if (SCM_REALP (y
))
4840 if (isnan (SCM_REAL_VALUE (y
)))
4842 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
4843 scm_remember_upto_here_1 (x
);
4844 return scm_from_bool (cmp
< 0);
4846 else if (SCM_FRACTIONP (y
))
4849 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4851 else if (SCM_REALP (x
))
4853 if (SCM_I_INUMP (y
))
4854 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
4855 else if (SCM_BIGP (y
))
4858 if (isnan (SCM_REAL_VALUE (x
)))
4860 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
4861 scm_remember_upto_here_1 (y
);
4862 return scm_from_bool (cmp
> 0);
4864 else if (SCM_REALP (y
))
4865 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
4866 else if (SCM_FRACTIONP (y
))
4868 double xx
= SCM_REAL_VALUE (x
);
4872 return scm_from_bool (xx
< 0.0);
4873 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
4877 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4879 else if (SCM_FRACTIONP (x
))
4881 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
4883 /* "a/b < y" becomes "a < y*b" */
4884 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
4885 x
= SCM_FRACTION_NUMERATOR (x
);
4888 else if (SCM_REALP (y
))
4890 double yy
= SCM_REAL_VALUE (y
);
4894 return scm_from_bool (0.0 < yy
);
4895 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
4898 else if (SCM_FRACTIONP (y
))
4900 /* "a/b < c/d" becomes "a*d < c*b" */
4901 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
4902 SCM_FRACTION_DENOMINATOR (y
));
4903 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
4904 SCM_FRACTION_DENOMINATOR (x
));
4910 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
4913 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
4917 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
4918 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
4919 (SCM x
, SCM y
, SCM rest
),
4920 "Return @code{#t} if the list of parameters is monotonically\n"
4922 #define FUNC_NAME s_scm_i_num_gr_p
4924 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4926 while (!scm_is_null (rest
))
4928 if (scm_is_false (scm_gr_p (x
, y
)))
4932 rest
= scm_cdr (rest
);
4934 return scm_gr_p (x
, y
);
4937 #define FUNC_NAME s_scm_i_num_gr_p
4939 scm_gr_p (SCM x
, SCM y
)
4941 if (!SCM_NUMBERP (x
))
4942 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
4943 else if (!SCM_NUMBERP (y
))
4944 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
4946 return scm_less_p (y
, x
);
4951 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
4952 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
4953 (SCM x
, SCM y
, SCM rest
),
4954 "Return @code{#t} if the list of parameters is monotonically\n"
4956 #define FUNC_NAME s_scm_i_num_leq_p
4958 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4960 while (!scm_is_null (rest
))
4962 if (scm_is_false (scm_leq_p (x
, y
)))
4966 rest
= scm_cdr (rest
);
4968 return scm_leq_p (x
, y
);
4971 #define FUNC_NAME s_scm_i_num_leq_p
4973 scm_leq_p (SCM x
, SCM y
)
4975 if (!SCM_NUMBERP (x
))
4976 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
4977 else if (!SCM_NUMBERP (y
))
4978 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
4979 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
4982 return scm_not (scm_less_p (y
, x
));
4987 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
4988 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
4989 (SCM x
, SCM y
, SCM rest
),
4990 "Return @code{#t} if the list of parameters is monotonically\n"
4992 #define FUNC_NAME s_scm_i_num_geq_p
4994 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
4996 while (!scm_is_null (rest
))
4998 if (scm_is_false (scm_geq_p (x
, y
)))
5002 rest
= scm_cdr (rest
);
5004 return scm_geq_p (x
, y
);
5007 #define FUNC_NAME s_scm_i_num_geq_p
5009 scm_geq_p (SCM x
, SCM y
)
5011 if (!SCM_NUMBERP (x
))
5012 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
5013 else if (!SCM_NUMBERP (y
))
5014 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
5015 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
5018 return scm_not (scm_less_p (x
, y
));
5023 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
5024 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
5030 if (SCM_I_INUMP (z
))
5031 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
5032 else if (SCM_BIGP (z
))
5034 else if (SCM_REALP (z
))
5035 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
5036 else if (SCM_COMPLEXP (z
))
5037 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
5038 && SCM_COMPLEX_IMAG (z
) == 0.0);
5039 else if (SCM_FRACTIONP (z
))
5042 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
5046 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
5047 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
5051 scm_positive_p (SCM x
)
5053 if (SCM_I_INUMP (x
))
5054 return scm_from_bool (SCM_I_INUM (x
) > 0);
5055 else if (SCM_BIGP (x
))
5057 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5058 scm_remember_upto_here_1 (x
);
5059 return scm_from_bool (sgn
> 0);
5061 else if (SCM_REALP (x
))
5062 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
5063 else if (SCM_FRACTIONP (x
))
5064 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
5066 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
5070 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
5071 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
5075 scm_negative_p (SCM x
)
5077 if (SCM_I_INUMP (x
))
5078 return scm_from_bool (SCM_I_INUM (x
) < 0);
5079 else if (SCM_BIGP (x
))
5081 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5082 scm_remember_upto_here_1 (x
);
5083 return scm_from_bool (sgn
< 0);
5085 else if (SCM_REALP (x
))
5086 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
5087 else if (SCM_FRACTIONP (x
))
5088 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
5090 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
5094 /* scm_min and scm_max return an inexact when either argument is inexact, as
5095 required by r5rs. On that basis, for exact/inexact combinations the
5096 exact is converted to inexact to compare and possibly return. This is
5097 unlike scm_less_p above which takes some trouble to preserve all bits in
5098 its test, such trouble is not required for min and max. */
5100 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
5101 (SCM x
, SCM y
, SCM rest
),
5102 "Return the maximum of all parameter values.")
5103 #define FUNC_NAME s_scm_i_max
5105 while (!scm_is_null (rest
))
5106 { x
= scm_max (x
, y
);
5108 rest
= scm_cdr (rest
);
5110 return scm_max (x
, y
);
5114 #define s_max s_scm_i_max
5115 #define g_max g_scm_i_max
5118 scm_max (SCM x
, SCM y
)
5123 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
5124 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
5127 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
5130 if (SCM_I_INUMP (x
))
5132 scm_t_inum xx
= SCM_I_INUM (x
);
5133 if (SCM_I_INUMP (y
))
5135 scm_t_inum yy
= SCM_I_INUM (y
);
5136 return (xx
< yy
) ? y
: x
;
5138 else if (SCM_BIGP (y
))
5140 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5141 scm_remember_upto_here_1 (y
);
5142 return (sgn
< 0) ? x
: y
;
5144 else if (SCM_REALP (y
))
5147 /* if y==NaN then ">" is false and we return NaN */
5148 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
5150 else if (SCM_FRACTIONP (y
))
5153 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
5156 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5158 else if (SCM_BIGP (x
))
5160 if (SCM_I_INUMP (y
))
5162 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5163 scm_remember_upto_here_1 (x
);
5164 return (sgn
< 0) ? y
: x
;
5166 else if (SCM_BIGP (y
))
5168 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
5169 scm_remember_upto_here_2 (x
, y
);
5170 return (cmp
> 0) ? x
: y
;
5172 else if (SCM_REALP (y
))
5174 /* if y==NaN then xx>yy is false, so we return the NaN y */
5177 xx
= scm_i_big2dbl (x
);
5178 yy
= SCM_REAL_VALUE (y
);
5179 return (xx
> yy
? scm_from_double (xx
) : y
);
5181 else if (SCM_FRACTIONP (y
))
5186 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5188 else if (SCM_REALP (x
))
5190 if (SCM_I_INUMP (y
))
5192 double z
= SCM_I_INUM (y
);
5193 /* if x==NaN then "<" is false and we return NaN */
5194 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
5196 else if (SCM_BIGP (y
))
5201 else if (SCM_REALP (y
))
5203 /* if x==NaN then our explicit check means we return NaN
5204 if y==NaN then ">" is false and we return NaN
5205 calling isnan is unavoidable, since it's the only way to know
5206 which of x or y causes any compares to be false */
5207 double xx
= SCM_REAL_VALUE (x
);
5208 return (isnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
5210 else if (SCM_FRACTIONP (y
))
5212 double yy
= scm_i_fraction2double (y
);
5213 double xx
= SCM_REAL_VALUE (x
);
5214 return (xx
< yy
) ? scm_from_double (yy
) : x
;
5217 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5219 else if (SCM_FRACTIONP (x
))
5221 if (SCM_I_INUMP (y
))
5225 else if (SCM_BIGP (y
))
5229 else if (SCM_REALP (y
))
5231 double xx
= scm_i_fraction2double (x
);
5232 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
5234 else if (SCM_FRACTIONP (y
))
5239 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
5242 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
5246 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
5247 (SCM x
, SCM y
, SCM rest
),
5248 "Return the minimum of all parameter values.")
5249 #define FUNC_NAME s_scm_i_min
5251 while (!scm_is_null (rest
))
5252 { x
= scm_min (x
, y
);
5254 rest
= scm_cdr (rest
);
5256 return scm_min (x
, y
);
5260 #define s_min s_scm_i_min
5261 #define g_min g_scm_i_min
5264 scm_min (SCM x
, SCM y
)
5269 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
5270 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
5273 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
5276 if (SCM_I_INUMP (x
))
5278 scm_t_inum xx
= SCM_I_INUM (x
);
5279 if (SCM_I_INUMP (y
))
5281 scm_t_inum yy
= SCM_I_INUM (y
);
5282 return (xx
< yy
) ? x
: y
;
5284 else if (SCM_BIGP (y
))
5286 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5287 scm_remember_upto_here_1 (y
);
5288 return (sgn
< 0) ? y
: x
;
5290 else if (SCM_REALP (y
))
5293 /* if y==NaN then "<" is false and we return NaN */
5294 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
5296 else if (SCM_FRACTIONP (y
))
5299 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
5302 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5304 else if (SCM_BIGP (x
))
5306 if (SCM_I_INUMP (y
))
5308 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5309 scm_remember_upto_here_1 (x
);
5310 return (sgn
< 0) ? x
: y
;
5312 else if (SCM_BIGP (y
))
5314 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
5315 scm_remember_upto_here_2 (x
, y
);
5316 return (cmp
> 0) ? y
: x
;
5318 else if (SCM_REALP (y
))
5320 /* if y==NaN then xx<yy is false, so we return the NaN y */
5323 xx
= scm_i_big2dbl (x
);
5324 yy
= SCM_REAL_VALUE (y
);
5325 return (xx
< yy
? scm_from_double (xx
) : y
);
5327 else if (SCM_FRACTIONP (y
))
5332 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5334 else if (SCM_REALP (x
))
5336 if (SCM_I_INUMP (y
))
5338 double z
= SCM_I_INUM (y
);
5339 /* if x==NaN then "<" is false and we return NaN */
5340 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
5342 else if (SCM_BIGP (y
))
5347 else if (SCM_REALP (y
))
5349 /* if x==NaN then our explicit check means we return NaN
5350 if y==NaN then "<" is false and we return NaN
5351 calling isnan is unavoidable, since it's the only way to know
5352 which of x or y causes any compares to be false */
5353 double xx
= SCM_REAL_VALUE (x
);
5354 return (isnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
5356 else if (SCM_FRACTIONP (y
))
5358 double yy
= scm_i_fraction2double (y
);
5359 double xx
= SCM_REAL_VALUE (x
);
5360 return (yy
< xx
) ? scm_from_double (yy
) : x
;
5363 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5365 else if (SCM_FRACTIONP (x
))
5367 if (SCM_I_INUMP (y
))
5371 else if (SCM_BIGP (y
))
5375 else if (SCM_REALP (y
))
5377 double xx
= scm_i_fraction2double (x
);
5378 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
5380 else if (SCM_FRACTIONP (y
))
5385 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
5388 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
5392 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
5393 (SCM x
, SCM y
, SCM rest
),
5394 "Return the sum of all parameter values. Return 0 if called without\n"
5396 #define FUNC_NAME s_scm_i_sum
5398 while (!scm_is_null (rest
))
5399 { x
= scm_sum (x
, y
);
5401 rest
= scm_cdr (rest
);
5403 return scm_sum (x
, y
);
5407 #define s_sum s_scm_i_sum
5408 #define g_sum g_scm_i_sum
5411 scm_sum (SCM x
, SCM y
)
5413 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5415 if (SCM_NUMBERP (x
)) return x
;
5416 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
5417 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
5420 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5422 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5424 scm_t_inum xx
= SCM_I_INUM (x
);
5425 scm_t_inum yy
= SCM_I_INUM (y
);
5426 scm_t_inum z
= xx
+ yy
;
5427 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_inum2big (z
);
5429 else if (SCM_BIGP (y
))
5434 else if (SCM_REALP (y
))
5436 scm_t_inum xx
= SCM_I_INUM (x
);
5437 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
5439 else if (SCM_COMPLEXP (y
))
5441 scm_t_inum xx
= SCM_I_INUM (x
);
5442 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
5443 SCM_COMPLEX_IMAG (y
));
5445 else if (SCM_FRACTIONP (y
))
5446 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
5447 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
5448 SCM_FRACTION_DENOMINATOR (y
));
5450 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5451 } else if (SCM_BIGP (x
))
5453 if (SCM_I_INUMP (y
))
5458 inum
= SCM_I_INUM (y
);
5461 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5464 SCM result
= scm_i_mkbig ();
5465 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
5466 scm_remember_upto_here_1 (x
);
5467 /* we know the result will have to be a bignum */
5470 return scm_i_normbig (result
);
5474 SCM result
= scm_i_mkbig ();
5475 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
5476 scm_remember_upto_here_1 (x
);
5477 /* we know the result will have to be a bignum */
5480 return scm_i_normbig (result
);
5483 else if (SCM_BIGP (y
))
5485 SCM result
= scm_i_mkbig ();
5486 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5487 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5488 mpz_add (SCM_I_BIG_MPZ (result
),
5491 scm_remember_upto_here_2 (x
, y
);
5492 /* we know the result will have to be a bignum */
5495 return scm_i_normbig (result
);
5497 else if (SCM_REALP (y
))
5499 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
5500 scm_remember_upto_here_1 (x
);
5501 return scm_from_double (result
);
5503 else if (SCM_COMPLEXP (y
))
5505 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
5506 + SCM_COMPLEX_REAL (y
));
5507 scm_remember_upto_here_1 (x
);
5508 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
5510 else if (SCM_FRACTIONP (y
))
5511 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
5512 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
5513 SCM_FRACTION_DENOMINATOR (y
));
5515 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5517 else if (SCM_REALP (x
))
5519 if (SCM_I_INUMP (y
))
5520 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
5521 else if (SCM_BIGP (y
))
5523 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
5524 scm_remember_upto_here_1 (y
);
5525 return scm_from_double (result
);
5527 else if (SCM_REALP (y
))
5528 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
5529 else if (SCM_COMPLEXP (y
))
5530 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
5531 SCM_COMPLEX_IMAG (y
));
5532 else if (SCM_FRACTIONP (y
))
5533 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
5535 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5537 else if (SCM_COMPLEXP (x
))
5539 if (SCM_I_INUMP (y
))
5540 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
5541 SCM_COMPLEX_IMAG (x
));
5542 else if (SCM_BIGP (y
))
5544 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
5545 + SCM_COMPLEX_REAL (x
));
5546 scm_remember_upto_here_1 (y
);
5547 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
5549 else if (SCM_REALP (y
))
5550 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
5551 SCM_COMPLEX_IMAG (x
));
5552 else if (SCM_COMPLEXP (y
))
5553 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
5554 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
5555 else if (SCM_FRACTIONP (y
))
5556 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
5557 SCM_COMPLEX_IMAG (x
));
5559 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5561 else if (SCM_FRACTIONP (x
))
5563 if (SCM_I_INUMP (y
))
5564 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
5565 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
5566 SCM_FRACTION_DENOMINATOR (x
));
5567 else if (SCM_BIGP (y
))
5568 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
5569 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
5570 SCM_FRACTION_DENOMINATOR (x
));
5571 else if (SCM_REALP (y
))
5572 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
5573 else if (SCM_COMPLEXP (y
))
5574 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
5575 SCM_COMPLEX_IMAG (y
));
5576 else if (SCM_FRACTIONP (y
))
5577 /* a/b + c/d = (ad + bc) / bd */
5578 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5579 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
5580 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
5582 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
5585 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
5589 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
5591 "Return @math{@var{x}+1}.")
5592 #define FUNC_NAME s_scm_oneplus
5594 return scm_sum (x
, SCM_INUM1
);
5599 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
5600 (SCM x
, SCM y
, SCM rest
),
5601 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
5602 "the sum of all but the first argument are subtracted from the first\n"
5604 #define FUNC_NAME s_scm_i_difference
5606 while (!scm_is_null (rest
))
5607 { x
= scm_difference (x
, y
);
5609 rest
= scm_cdr (rest
);
5611 return scm_difference (x
, y
);
5615 #define s_difference s_scm_i_difference
5616 #define g_difference g_scm_i_difference
5619 scm_difference (SCM x
, SCM y
)
5620 #define FUNC_NAME s_difference
5622 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5625 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
5627 if (SCM_I_INUMP (x
))
5629 scm_t_inum xx
= -SCM_I_INUM (x
);
5630 if (SCM_FIXABLE (xx
))
5631 return SCM_I_MAKINUM (xx
);
5633 return scm_i_inum2big (xx
);
5635 else if (SCM_BIGP (x
))
5636 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
5637 bignum, but negating that gives a fixnum. */
5638 return scm_i_normbig (scm_i_clonebig (x
, 0));
5639 else if (SCM_REALP (x
))
5640 return scm_from_double (-SCM_REAL_VALUE (x
));
5641 else if (SCM_COMPLEXP (x
))
5642 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
5643 -SCM_COMPLEX_IMAG (x
));
5644 else if (SCM_FRACTIONP (x
))
5645 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
5646 SCM_FRACTION_DENOMINATOR (x
));
5648 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
5651 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5653 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5655 scm_t_inum xx
= SCM_I_INUM (x
);
5656 scm_t_inum yy
= SCM_I_INUM (y
);
5657 scm_t_inum z
= xx
- yy
;
5658 if (SCM_FIXABLE (z
))
5659 return SCM_I_MAKINUM (z
);
5661 return scm_i_inum2big (z
);
5663 else if (SCM_BIGP (y
))
5665 /* inum-x - big-y */
5666 scm_t_inum xx
= SCM_I_INUM (x
);
5670 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
5671 bignum, but negating that gives a fixnum. */
5672 return scm_i_normbig (scm_i_clonebig (y
, 0));
5676 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5677 SCM result
= scm_i_mkbig ();
5680 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
5683 /* x - y == -(y + -x) */
5684 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
5685 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5687 scm_remember_upto_here_1 (y
);
5689 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
5690 /* we know the result will have to be a bignum */
5693 return scm_i_normbig (result
);
5696 else if (SCM_REALP (y
))
5698 scm_t_inum xx
= SCM_I_INUM (x
);
5699 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
5701 else if (SCM_COMPLEXP (y
))
5703 scm_t_inum xx
= SCM_I_INUM (x
);
5704 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
5705 - SCM_COMPLEX_IMAG (y
));
5707 else if (SCM_FRACTIONP (y
))
5708 /* a - b/c = (ac - b) / c */
5709 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5710 SCM_FRACTION_NUMERATOR (y
)),
5711 SCM_FRACTION_DENOMINATOR (y
));
5713 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5715 else if (SCM_BIGP (x
))
5717 if (SCM_I_INUMP (y
))
5719 /* big-x - inum-y */
5720 scm_t_inum yy
= SCM_I_INUM (y
);
5721 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5723 scm_remember_upto_here_1 (x
);
5725 return (SCM_FIXABLE (-yy
) ?
5726 SCM_I_MAKINUM (-yy
) : scm_from_inum (-yy
));
5729 SCM result
= scm_i_mkbig ();
5732 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
5734 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
5735 scm_remember_upto_here_1 (x
);
5737 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
5738 /* we know the result will have to be a bignum */
5741 return scm_i_normbig (result
);
5744 else if (SCM_BIGP (y
))
5746 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5747 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
5748 SCM result
= scm_i_mkbig ();
5749 mpz_sub (SCM_I_BIG_MPZ (result
),
5752 scm_remember_upto_here_2 (x
, y
);
5753 /* we know the result will have to be a bignum */
5754 if ((sgn_x
== 1) && (sgn_y
== -1))
5756 if ((sgn_x
== -1) && (sgn_y
== 1))
5758 return scm_i_normbig (result
);
5760 else if (SCM_REALP (y
))
5762 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
5763 scm_remember_upto_here_1 (x
);
5764 return scm_from_double (result
);
5766 else if (SCM_COMPLEXP (y
))
5768 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
5769 - SCM_COMPLEX_REAL (y
));
5770 scm_remember_upto_here_1 (x
);
5771 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
5773 else if (SCM_FRACTIONP (y
))
5774 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5775 SCM_FRACTION_NUMERATOR (y
)),
5776 SCM_FRACTION_DENOMINATOR (y
));
5777 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5779 else if (SCM_REALP (x
))
5781 if (SCM_I_INUMP (y
))
5782 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
5783 else if (SCM_BIGP (y
))
5785 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
5786 scm_remember_upto_here_1 (x
);
5787 return scm_from_double (result
);
5789 else if (SCM_REALP (y
))
5790 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
5791 else if (SCM_COMPLEXP (y
))
5792 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
5793 -SCM_COMPLEX_IMAG (y
));
5794 else if (SCM_FRACTIONP (y
))
5795 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
5797 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5799 else if (SCM_COMPLEXP (x
))
5801 if (SCM_I_INUMP (y
))
5802 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
5803 SCM_COMPLEX_IMAG (x
));
5804 else if (SCM_BIGP (y
))
5806 double real_part
= (SCM_COMPLEX_REAL (x
)
5807 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
5808 scm_remember_upto_here_1 (x
);
5809 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
5811 else if (SCM_REALP (y
))
5812 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
5813 SCM_COMPLEX_IMAG (x
));
5814 else if (SCM_COMPLEXP (y
))
5815 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
5816 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
5817 else if (SCM_FRACTIONP (y
))
5818 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
5819 SCM_COMPLEX_IMAG (x
));
5821 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5823 else if (SCM_FRACTIONP (x
))
5825 if (SCM_I_INUMP (y
))
5826 /* a/b - c = (a - cb) / b */
5827 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
5828 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
5829 SCM_FRACTION_DENOMINATOR (x
));
5830 else if (SCM_BIGP (y
))
5831 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
5832 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
5833 SCM_FRACTION_DENOMINATOR (x
));
5834 else if (SCM_REALP (y
))
5835 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
5836 else if (SCM_COMPLEXP (y
))
5837 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
5838 -SCM_COMPLEX_IMAG (y
));
5839 else if (SCM_FRACTIONP (y
))
5840 /* a/b - c/d = (ad - bc) / bd */
5841 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5842 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
5843 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
5845 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
5848 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
5853 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
5855 "Return @math{@var{x}-1}.")
5856 #define FUNC_NAME s_scm_oneminus
5858 return scm_difference (x
, SCM_INUM1
);
5863 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
5864 (SCM x
, SCM y
, SCM rest
),
5865 "Return the product of all arguments. If called without arguments,\n"
5867 #define FUNC_NAME s_scm_i_product
5869 while (!scm_is_null (rest
))
5870 { x
= scm_product (x
, y
);
5872 rest
= scm_cdr (rest
);
5874 return scm_product (x
, y
);
5878 #define s_product s_scm_i_product
5879 #define g_product g_scm_i_product
5882 scm_product (SCM x
, SCM y
)
5884 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
5887 return SCM_I_MAKINUM (1L);
5888 else if (SCM_NUMBERP (x
))
5891 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
5894 if (SCM_LIKELY (SCM_I_INUMP (x
)))
5899 xx
= SCM_I_INUM (x
);
5903 case 0: return x
; break;
5904 case 1: return y
; break;
5906 * The following case (x = -1) is important for more than
5907 * just optimization. It handles the case of negating
5908 * (+ 1 most-positive-fixnum) aka (- most-negative-fixnum),
5909 * which is a bignum that must be changed back into a fixnum.
5910 * Failure to do so will cause the following to return #f:
5911 * (= most-negative-fixnum (* -1 (- most-negative-fixnum)))
5914 return scm_difference(y
, SCM_UNDEFINED
);
5918 if (SCM_LIKELY (SCM_I_INUMP (y
)))
5920 scm_t_inum yy
= SCM_I_INUM (y
);
5921 scm_t_inum kk
= xx
* yy
;
5922 SCM k
= SCM_I_MAKINUM (kk
);
5923 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
5927 SCM result
= scm_i_inum2big (xx
);
5928 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
5929 return scm_i_normbig (result
);
5932 else if (SCM_BIGP (y
))
5934 SCM result
= scm_i_mkbig ();
5935 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
5936 scm_remember_upto_here_1 (y
);
5939 else if (SCM_REALP (y
))
5940 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
5941 else if (SCM_COMPLEXP (y
))
5942 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
5943 xx
* SCM_COMPLEX_IMAG (y
));
5944 else if (SCM_FRACTIONP (y
))
5945 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
5946 SCM_FRACTION_DENOMINATOR (y
));
5948 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
5950 else if (SCM_BIGP (x
))
5952 if (SCM_I_INUMP (y
))
5957 else if (SCM_BIGP (y
))
5959 SCM result
= scm_i_mkbig ();
5960 mpz_mul (SCM_I_BIG_MPZ (result
),
5963 scm_remember_upto_here_2 (x
, y
);
5966 else if (SCM_REALP (y
))
5968 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
5969 scm_remember_upto_here_1 (x
);
5970 return scm_from_double (result
);
5972 else if (SCM_COMPLEXP (y
))
5974 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
5975 scm_remember_upto_here_1 (x
);
5976 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
5977 z
* SCM_COMPLEX_IMAG (y
));
5979 else if (SCM_FRACTIONP (y
))
5980 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
5981 SCM_FRACTION_DENOMINATOR (y
));
5983 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
5985 else if (SCM_REALP (x
))
5987 if (SCM_I_INUMP (y
))
5989 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
5990 if (scm_is_eq (y
, SCM_INUM0
))
5992 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
5994 else if (SCM_BIGP (y
))
5996 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
5997 scm_remember_upto_here_1 (y
);
5998 return scm_from_double (result
);
6000 else if (SCM_REALP (y
))
6001 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
6002 else if (SCM_COMPLEXP (y
))
6003 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
6004 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
6005 else if (SCM_FRACTIONP (y
))
6006 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
6008 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6010 else if (SCM_COMPLEXP (x
))
6012 if (SCM_I_INUMP (y
))
6014 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
6015 if (scm_is_eq (y
, SCM_INUM0
))
6017 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
6018 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
6020 else if (SCM_BIGP (y
))
6022 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6023 scm_remember_upto_here_1 (y
);
6024 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
6025 z
* SCM_COMPLEX_IMAG (x
));
6027 else if (SCM_REALP (y
))
6028 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
6029 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
6030 else if (SCM_COMPLEXP (y
))
6032 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
6033 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
6034 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
6035 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
6037 else if (SCM_FRACTIONP (y
))
6039 double yy
= scm_i_fraction2double (y
);
6040 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
6041 yy
* SCM_COMPLEX_IMAG (x
));
6044 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6046 else if (SCM_FRACTIONP (x
))
6048 if (SCM_I_INUMP (y
))
6049 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
6050 SCM_FRACTION_DENOMINATOR (x
));
6051 else if (SCM_BIGP (y
))
6052 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
6053 SCM_FRACTION_DENOMINATOR (x
));
6054 else if (SCM_REALP (y
))
6055 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
6056 else if (SCM_COMPLEXP (y
))
6058 double xx
= scm_i_fraction2double (x
);
6059 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
6060 xx
* SCM_COMPLEX_IMAG (y
));
6062 else if (SCM_FRACTIONP (y
))
6063 /* a/b * c/d = ac / bd */
6064 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
6065 SCM_FRACTION_NUMERATOR (y
)),
6066 scm_product (SCM_FRACTION_DENOMINATOR (x
),
6067 SCM_FRACTION_DENOMINATOR (y
)));
6069 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
6072 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
6075 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
6076 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
6077 #define ALLOW_DIVIDE_BY_ZERO
6078 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
6081 /* The code below for complex division is adapted from the GNU
6082 libstdc++, which adapted it from f2c's libF77, and is subject to
6085 /****************************************************************
6086 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
6088 Permission to use, copy, modify, and distribute this software
6089 and its documentation for any purpose and without fee is hereby
6090 granted, provided that the above copyright notice appear in all
6091 copies and that both that the copyright notice and this
6092 permission notice and warranty disclaimer appear in supporting
6093 documentation, and that the names of AT&T Bell Laboratories or
6094 Bellcore or any of their entities not be used in advertising or
6095 publicity pertaining to distribution of the software without
6096 specific, written prior permission.
6098 AT&T and Bellcore disclaim all warranties with regard to this
6099 software, including all implied warranties of merchantability
6100 and fitness. In no event shall AT&T or Bellcore be liable for
6101 any special, indirect or consequential damages or any damages
6102 whatsoever resulting from loss of use, data or profits, whether
6103 in an action of contract, negligence or other tortious action,
6104 arising out of or in connection with the use or performance of
6106 ****************************************************************/
6108 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
6109 (SCM x
, SCM y
, SCM rest
),
6110 "Divide the first argument by the product of the remaining\n"
6111 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
6113 #define FUNC_NAME s_scm_i_divide
6115 while (!scm_is_null (rest
))
6116 { x
= scm_divide (x
, y
);
6118 rest
= scm_cdr (rest
);
6120 return scm_divide (x
, y
);
6124 #define s_divide s_scm_i_divide
6125 #define g_divide g_scm_i_divide
6128 do_divide (SCM x
, SCM y
, int inexact
)
6129 #define FUNC_NAME s_divide
6133 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
6136 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
6137 else if (SCM_I_INUMP (x
))
6139 scm_t_inum xx
= SCM_I_INUM (x
);
6140 if (xx
== 1 || xx
== -1)
6142 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6144 scm_num_overflow (s_divide
);
6149 return scm_from_double (1.0 / (double) xx
);
6150 else return scm_i_make_ratio (SCM_INUM1
, x
);
6153 else if (SCM_BIGP (x
))
6156 return scm_from_double (1.0 / scm_i_big2dbl (x
));
6157 else return scm_i_make_ratio (SCM_INUM1
, x
);
6159 else if (SCM_REALP (x
))
6161 double xx
= SCM_REAL_VALUE (x
);
6162 #ifndef ALLOW_DIVIDE_BY_ZERO
6164 scm_num_overflow (s_divide
);
6167 return scm_from_double (1.0 / xx
);
6169 else if (SCM_COMPLEXP (x
))
6171 double r
= SCM_COMPLEX_REAL (x
);
6172 double i
= SCM_COMPLEX_IMAG (x
);
6173 if (fabs(r
) <= fabs(i
))
6176 double d
= i
* (1.0 + t
* t
);
6177 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
6182 double d
= r
* (1.0 + t
* t
);
6183 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
6186 else if (SCM_FRACTIONP (x
))
6187 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
6188 SCM_FRACTION_NUMERATOR (x
));
6190 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
6193 if (SCM_LIKELY (SCM_I_INUMP (x
)))
6195 scm_t_inum xx
= SCM_I_INUM (x
);
6196 if (SCM_LIKELY (SCM_I_INUMP (y
)))
6198 scm_t_inum yy
= SCM_I_INUM (y
);
6201 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6202 scm_num_overflow (s_divide
);
6204 return scm_from_double ((double) xx
/ (double) yy
);
6207 else if (xx
% yy
!= 0)
6210 return scm_from_double ((double) xx
/ (double) yy
);
6211 else return scm_i_make_ratio (x
, y
);
6215 scm_t_inum z
= xx
/ yy
;
6216 if (SCM_FIXABLE (z
))
6217 return SCM_I_MAKINUM (z
);
6219 return scm_i_inum2big (z
);
6222 else if (SCM_BIGP (y
))
6225 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
6226 else return scm_i_make_ratio (x
, y
);
6228 else if (SCM_REALP (y
))
6230 double yy
= SCM_REAL_VALUE (y
);
6231 #ifndef ALLOW_DIVIDE_BY_ZERO
6233 scm_num_overflow (s_divide
);
6236 return scm_from_double ((double) xx
/ yy
);
6238 else if (SCM_COMPLEXP (y
))
6241 complex_div
: /* y _must_ be a complex number */
6243 double r
= SCM_COMPLEX_REAL (y
);
6244 double i
= SCM_COMPLEX_IMAG (y
);
6245 if (fabs(r
) <= fabs(i
))
6248 double d
= i
* (1.0 + t
* t
);
6249 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
6254 double d
= r
* (1.0 + t
* t
);
6255 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
6259 else if (SCM_FRACTIONP (y
))
6260 /* a / b/c = ac / b */
6261 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
6262 SCM_FRACTION_NUMERATOR (y
));
6264 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6266 else if (SCM_BIGP (x
))
6268 if (SCM_I_INUMP (y
))
6270 scm_t_inum yy
= SCM_I_INUM (y
);
6273 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6274 scm_num_overflow (s_divide
);
6276 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
6277 scm_remember_upto_here_1 (x
);
6278 return (sgn
== 0) ? scm_nan () : scm_inf ();
6285 /* FIXME: HMM, what are the relative performance issues here?
6286 We need to test. Is it faster on average to test
6287 divisible_p, then perform whichever operation, or is it
6288 faster to perform the integer div opportunistically and
6289 switch to real if there's a remainder? For now we take the
6290 middle ground: test, then if divisible, use the faster div
6293 scm_t_inum abs_yy
= yy
< 0 ? -yy
: yy
;
6294 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
6298 SCM result
= scm_i_mkbig ();
6299 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
6300 scm_remember_upto_here_1 (x
);
6302 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
6303 return scm_i_normbig (result
);
6308 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
6309 else return scm_i_make_ratio (x
, y
);
6313 else if (SCM_BIGP (y
))
6318 /* It's easily possible for the ratio x/y to fit a double
6319 but one or both x and y be too big to fit a double,
6320 hence the use of mpq_get_d rather than converting and
6323 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
6324 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
6325 return scm_from_double (mpq_get_d (q
));
6329 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
6333 SCM result
= scm_i_mkbig ();
6334 mpz_divexact (SCM_I_BIG_MPZ (result
),
6337 scm_remember_upto_here_2 (x
, y
);
6338 return scm_i_normbig (result
);
6341 return scm_i_make_ratio (x
, y
);
6344 else if (SCM_REALP (y
))
6346 double yy
= SCM_REAL_VALUE (y
);
6347 #ifndef ALLOW_DIVIDE_BY_ZERO
6349 scm_num_overflow (s_divide
);
6352 return scm_from_double (scm_i_big2dbl (x
) / yy
);
6354 else if (SCM_COMPLEXP (y
))
6356 a
= scm_i_big2dbl (x
);
6359 else if (SCM_FRACTIONP (y
))
6360 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
6361 SCM_FRACTION_NUMERATOR (y
));
6363 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6365 else if (SCM_REALP (x
))
6367 double rx
= SCM_REAL_VALUE (x
);
6368 if (SCM_I_INUMP (y
))
6370 scm_t_inum yy
= SCM_I_INUM (y
);
6371 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6373 scm_num_overflow (s_divide
);
6376 return scm_from_double (rx
/ (double) yy
);
6378 else if (SCM_BIGP (y
))
6380 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6381 scm_remember_upto_here_1 (y
);
6382 return scm_from_double (rx
/ dby
);
6384 else if (SCM_REALP (y
))
6386 double yy
= SCM_REAL_VALUE (y
);
6387 #ifndef ALLOW_DIVIDE_BY_ZERO
6389 scm_num_overflow (s_divide
);
6392 return scm_from_double (rx
/ yy
);
6394 else if (SCM_COMPLEXP (y
))
6399 else if (SCM_FRACTIONP (y
))
6400 return scm_from_double (rx
/ scm_i_fraction2double (y
));
6402 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6404 else if (SCM_COMPLEXP (x
))
6406 double rx
= SCM_COMPLEX_REAL (x
);
6407 double ix
= SCM_COMPLEX_IMAG (x
);
6408 if (SCM_I_INUMP (y
))
6410 scm_t_inum yy
= SCM_I_INUM (y
);
6411 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6413 scm_num_overflow (s_divide
);
6418 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
6421 else if (SCM_BIGP (y
))
6423 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
6424 scm_remember_upto_here_1 (y
);
6425 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
6427 else if (SCM_REALP (y
))
6429 double yy
= SCM_REAL_VALUE (y
);
6430 #ifndef ALLOW_DIVIDE_BY_ZERO
6432 scm_num_overflow (s_divide
);
6435 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
6437 else if (SCM_COMPLEXP (y
))
6439 double ry
= SCM_COMPLEX_REAL (y
);
6440 double iy
= SCM_COMPLEX_IMAG (y
);
6441 if (fabs(ry
) <= fabs(iy
))
6444 double d
= iy
* (1.0 + t
* t
);
6445 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
6450 double d
= ry
* (1.0 + t
* t
);
6451 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
6454 else if (SCM_FRACTIONP (y
))
6456 double yy
= scm_i_fraction2double (y
);
6457 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
6460 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6462 else if (SCM_FRACTIONP (x
))
6464 if (SCM_I_INUMP (y
))
6466 scm_t_inum yy
= SCM_I_INUM (y
);
6467 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
6469 scm_num_overflow (s_divide
);
6472 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
6473 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
6475 else if (SCM_BIGP (y
))
6477 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
6478 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
6480 else if (SCM_REALP (y
))
6482 double yy
= SCM_REAL_VALUE (y
);
6483 #ifndef ALLOW_DIVIDE_BY_ZERO
6485 scm_num_overflow (s_divide
);
6488 return scm_from_double (scm_i_fraction2double (x
) / yy
);
6490 else if (SCM_COMPLEXP (y
))
6492 a
= scm_i_fraction2double (x
);
6495 else if (SCM_FRACTIONP (y
))
6496 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
6497 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
6499 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
6502 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
6506 scm_divide (SCM x
, SCM y
)
6508 return do_divide (x
, y
, 0);
6511 static SCM
scm_divide2real (SCM x
, SCM y
)
6513 return do_divide (x
, y
, 1);
6519 scm_c_truncate (double x
)
6530 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
6531 half-way case (ie. when x is an integer plus 0.5) going upwards.
6532 Then half-way cases are identified and adjusted down if the
6533 round-upwards didn't give the desired even integer.
6535 "plus_half == result" identifies a half-way case. If plus_half, which is
6536 x + 0.5, is an integer then x must be an integer plus 0.5.
6538 An odd "result" value is identified with result/2 != floor(result/2).
6539 This is done with plus_half, since that value is ready for use sooner in
6540 a pipelined cpu, and we're already requiring plus_half == result.
6542 Note however that we need to be careful when x is big and already an
6543 integer. In that case "x+0.5" may round to an adjacent integer, causing
6544 us to return such a value, incorrectly. For instance if the hardware is
6545 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
6546 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
6547 returned. Or if the hardware is in round-upwards mode, then other bigger
6548 values like say x == 2^128 will see x+0.5 rounding up to the next higher
6549 representable value, 2^128+2^76 (or whatever), again incorrect.
6551 These bad roundings of x+0.5 are avoided by testing at the start whether
6552 x is already an integer. If it is then clearly that's the desired result
6553 already. And if it's not then the exponent must be small enough to allow
6554 an 0.5 to be represented, and hence added without a bad rounding. */
6557 scm_c_round (double x
)
6559 double plus_half
, result
;
6564 plus_half
= x
+ 0.5;
6565 result
= floor (plus_half
);
6566 /* Adjust so that the rounding is towards even. */
6567 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
6572 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
6574 "Round the number @var{x} towards zero.")
6575 #define FUNC_NAME s_scm_truncate_number
6577 if (scm_is_false (scm_negative_p (x
)))
6578 return scm_floor (x
);
6580 return scm_ceiling (x
);
6584 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
6586 "Round the number @var{x} towards the nearest integer. "
6587 "When it is exactly halfway between two integers, "
6588 "round towards the even one.")
6589 #define FUNC_NAME s_scm_round_number
6591 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6593 else if (SCM_REALP (x
))
6594 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
6597 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
6598 single quotient+remainder division then examining to see which way
6599 the rounding should go. */
6600 SCM plus_half
= scm_sum (x
, exactly_one_half
);
6601 SCM result
= scm_floor (plus_half
);
6602 /* Adjust so that the rounding is towards even. */
6603 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
6604 && scm_is_true (scm_odd_p (result
)))
6605 return scm_difference (result
, SCM_INUM1
);
6612 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
6614 "Round the number @var{x} towards minus infinity.")
6615 #define FUNC_NAME s_scm_floor
6617 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6619 else if (SCM_REALP (x
))
6620 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
6621 else if (SCM_FRACTIONP (x
))
6623 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
6624 SCM_FRACTION_DENOMINATOR (x
));
6625 if (scm_is_false (scm_negative_p (x
)))
6627 /* For positive x, rounding towards zero is correct. */
6632 /* For negative x, we need to return q-1 unless x is an
6633 integer. But fractions are never integer, per our
6635 return scm_difference (q
, SCM_INUM1
);
6639 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
6643 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
6645 "Round the number @var{x} towards infinity.")
6646 #define FUNC_NAME s_scm_ceiling
6648 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
6650 else if (SCM_REALP (x
))
6651 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
6652 else if (SCM_FRACTIONP (x
))
6654 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
6655 SCM_FRACTION_DENOMINATOR (x
));
6656 if (scm_is_false (scm_positive_p (x
)))
6658 /* For negative x, rounding towards zero is correct. */
6663 /* For positive x, we need to return q+1 unless x is an
6664 integer. But fractions are never integer, per our
6666 return scm_sum (q
, SCM_INUM1
);
6670 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
6674 /* sin/cos/tan/asin/acos/atan
6675 sinh/cosh/tanh/asinh/acosh/atanh
6676 Derived from "Transcen.scm", Complex trancendental functions for SCM.
6677 Written by Jerry D. Hedden, (C) FSF.
6678 See the file `COPYING' for terms applying to this program. */
6680 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
6682 "Return @var{x} raised to the power of @var{y}.")
6683 #define FUNC_NAME s_scm_expt
6685 if (scm_is_integer (y
))
6687 if (scm_is_true (scm_exact_p (y
)))
6688 return scm_integer_expt (x
, y
);
6691 /* Here we handle the case where the exponent is an inexact
6692 integer. We make the exponent exact in order to use
6693 scm_integer_expt, and thus avoid the spurious imaginary
6694 parts that may result from round-off errors in the general
6695 e^(y log x) method below (for example when squaring a large
6696 negative number). In this case, we must return an inexact
6697 result for correctness. We also make the base inexact so
6698 that scm_integer_expt will use fast inexact arithmetic
6699 internally. Note that making the base inexact is not
6700 sufficient to guarantee an inexact result, because
6701 scm_integer_expt will return an exact 1 when the exponent
6702 is 0, even if the base is inexact. */
6703 return scm_exact_to_inexact
6704 (scm_integer_expt (scm_exact_to_inexact (x
),
6705 scm_inexact_to_exact (y
)));
6708 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
6710 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
6713 return scm_exp (scm_product (scm_log (x
), y
));
6717 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
6719 "Compute the sine of @var{z}.")
6720 #define FUNC_NAME s_scm_sin
6722 if (scm_is_real (z
))
6723 return scm_from_double (sin (scm_to_double (z
)));
6724 else if (SCM_COMPLEXP (z
))
6726 x
= SCM_COMPLEX_REAL (z
);
6727 y
= SCM_COMPLEX_IMAG (z
);
6728 return scm_c_make_rectangular (sin (x
) * cosh (y
),
6729 cos (x
) * sinh (y
));
6732 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
6736 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
6738 "Compute the cosine of @var{z}.")
6739 #define FUNC_NAME s_scm_cos
6741 if (scm_is_real (z
))
6742 return scm_from_double (cos (scm_to_double (z
)));
6743 else if (SCM_COMPLEXP (z
))
6745 x
= SCM_COMPLEX_REAL (z
);
6746 y
= SCM_COMPLEX_IMAG (z
);
6747 return scm_c_make_rectangular (cos (x
) * cosh (y
),
6748 -sin (x
) * sinh (y
));
6751 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
6755 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
6757 "Compute the tangent of @var{z}.")
6758 #define FUNC_NAME s_scm_tan
6760 if (scm_is_real (z
))
6761 return scm_from_double (tan (scm_to_double (z
)));
6762 else if (SCM_COMPLEXP (z
))
6764 x
= 2.0 * SCM_COMPLEX_REAL (z
);
6765 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
6766 w
= cos (x
) + cosh (y
);
6767 #ifndef ALLOW_DIVIDE_BY_ZERO
6769 scm_num_overflow (s_scm_tan
);
6771 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
6774 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
6778 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
6780 "Compute the hyperbolic sine of @var{z}.")
6781 #define FUNC_NAME s_scm_sinh
6783 if (scm_is_real (z
))
6784 return scm_from_double (sinh (scm_to_double (z
)));
6785 else if (SCM_COMPLEXP (z
))
6787 x
= SCM_COMPLEX_REAL (z
);
6788 y
= SCM_COMPLEX_IMAG (z
);
6789 return scm_c_make_rectangular (sinh (x
) * cos (y
),
6790 cosh (x
) * sin (y
));
6793 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
6797 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
6799 "Compute the hyperbolic cosine of @var{z}.")
6800 #define FUNC_NAME s_scm_cosh
6802 if (scm_is_real (z
))
6803 return scm_from_double (cosh (scm_to_double (z
)));
6804 else if (SCM_COMPLEXP (z
))
6806 x
= SCM_COMPLEX_REAL (z
);
6807 y
= SCM_COMPLEX_IMAG (z
);
6808 return scm_c_make_rectangular (cosh (x
) * cos (y
),
6809 sinh (x
) * sin (y
));
6812 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
6816 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
6818 "Compute the hyperbolic tangent of @var{z}.")
6819 #define FUNC_NAME s_scm_tanh
6821 if (scm_is_real (z
))
6822 return scm_from_double (tanh (scm_to_double (z
)));
6823 else if (SCM_COMPLEXP (z
))
6825 x
= 2.0 * SCM_COMPLEX_REAL (z
);
6826 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
6827 w
= cosh (x
) + cos (y
);
6828 #ifndef ALLOW_DIVIDE_BY_ZERO
6830 scm_num_overflow (s_scm_tanh
);
6832 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
6835 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
6839 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
6841 "Compute the arc sine of @var{z}.")
6842 #define FUNC_NAME s_scm_asin
6844 if (scm_is_real (z
))
6846 double w
= scm_to_double (z
);
6847 if (w
>= -1.0 && w
<= 1.0)
6848 return scm_from_double (asin (w
));
6850 return scm_product (scm_c_make_rectangular (0, -1),
6851 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
6853 else if (SCM_COMPLEXP (z
))
6855 x
= SCM_COMPLEX_REAL (z
);
6856 y
= SCM_COMPLEX_IMAG (z
);
6857 return scm_product (scm_c_make_rectangular (0, -1),
6858 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
6861 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
6865 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
6867 "Compute the arc cosine of @var{z}.")
6868 #define FUNC_NAME s_scm_acos
6870 if (scm_is_real (z
))
6872 double w
= scm_to_double (z
);
6873 if (w
>= -1.0 && w
<= 1.0)
6874 return scm_from_double (acos (w
));
6876 return scm_sum (scm_from_double (acos (0.0)),
6877 scm_product (scm_c_make_rectangular (0, 1),
6878 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
6880 else if (SCM_COMPLEXP (z
))
6882 x
= SCM_COMPLEX_REAL (z
);
6883 y
= SCM_COMPLEX_IMAG (z
);
6884 return scm_sum (scm_from_double (acos (0.0)),
6885 scm_product (scm_c_make_rectangular (0, 1),
6886 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
6889 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
6893 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
6895 "With one argument, compute the arc tangent of @var{z}.\n"
6896 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
6897 "using the sign of @var{z} and @var{y} to determine the quadrant.")
6898 #define FUNC_NAME s_scm_atan
6902 if (scm_is_real (z
))
6903 return scm_from_double (atan (scm_to_double (z
)));
6904 else if (SCM_COMPLEXP (z
))
6907 v
= SCM_COMPLEX_REAL (z
);
6908 w
= SCM_COMPLEX_IMAG (z
);
6909 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
6910 scm_c_make_rectangular (v
, w
+ 1.0))),
6911 scm_c_make_rectangular (0, 2));
6914 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
6916 else if (scm_is_real (z
))
6918 if (scm_is_real (y
))
6919 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
6921 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
6924 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
6928 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
6930 "Compute the inverse hyperbolic sine of @var{z}.")
6931 #define FUNC_NAME s_scm_sys_asinh
6933 if (scm_is_real (z
))
6934 return scm_from_double (asinh (scm_to_double (z
)));
6935 else if (scm_is_number (z
))
6936 return scm_log (scm_sum (z
,
6937 scm_sqrt (scm_sum (scm_product (z
, z
),
6940 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
6944 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
6946 "Compute the inverse hyperbolic cosine of @var{z}.")
6947 #define FUNC_NAME s_scm_sys_acosh
6949 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
6950 return scm_from_double (acosh (scm_to_double (z
)));
6951 else if (scm_is_number (z
))
6952 return scm_log (scm_sum (z
,
6953 scm_sqrt (scm_difference (scm_product (z
, z
),
6956 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
6960 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
6962 "Compute the inverse hyperbolic tangent of @var{z}.")
6963 #define FUNC_NAME s_scm_sys_atanh
6965 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
6966 return scm_from_double (atanh (scm_to_double (z
)));
6967 else if (scm_is_number (z
))
6968 return scm_divide (scm_log (scm_divide (scm_sum (SCM_INUM1
, z
),
6969 scm_difference (SCM_INUM1
, z
))),
6972 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
6977 scm_c_make_rectangular (double re
, double im
)
6980 return scm_from_double (re
);
6985 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
6987 SCM_SET_CELL_TYPE (z
, scm_tc16_complex
);
6988 SCM_COMPLEX_REAL (z
) = re
;
6989 SCM_COMPLEX_IMAG (z
) = im
;
6994 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
6995 (SCM real_part
, SCM imaginary_part
),
6996 "Return a complex number constructed of the given @var{real-part} "
6997 "and @var{imaginary-part} parts.")
6998 #define FUNC_NAME s_scm_make_rectangular
7000 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
7001 SCM_ARG1
, FUNC_NAME
, "real");
7002 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
7003 SCM_ARG2
, FUNC_NAME
, "real");
7004 return scm_c_make_rectangular (scm_to_double (real_part
),
7005 scm_to_double (imaginary_part
));
7010 scm_c_make_polar (double mag
, double ang
)
7014 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
7015 use it on Glibc-based systems that have it (it's a GNU extension). See
7016 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
7018 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
7019 sincos (ang
, &s
, &c
);
7024 return scm_c_make_rectangular (mag
* c
, mag
* s
);
7027 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
7029 "Return the complex number @var{x} * e^(i * @var{y}).")
7030 #define FUNC_NAME s_scm_make_polar
7032 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
7033 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
7034 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
7039 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
7040 /* "Return the real part of the number @var{z}."
7043 scm_real_part (SCM z
)
7045 if (SCM_I_INUMP (z
))
7047 else if (SCM_BIGP (z
))
7049 else if (SCM_REALP (z
))
7051 else if (SCM_COMPLEXP (z
))
7052 return scm_from_double (SCM_COMPLEX_REAL (z
));
7053 else if (SCM_FRACTIONP (z
))
7056 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
7060 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
7061 /* "Return the imaginary part of the number @var{z}."
7064 scm_imag_part (SCM z
)
7066 if (SCM_I_INUMP (z
))
7068 else if (SCM_BIGP (z
))
7070 else if (SCM_REALP (z
))
7072 else if (SCM_COMPLEXP (z
))
7073 return scm_from_double (SCM_COMPLEX_IMAG (z
));
7074 else if (SCM_FRACTIONP (z
))
7077 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
7080 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
7081 /* "Return the numerator of the number @var{z}."
7084 scm_numerator (SCM z
)
7086 if (SCM_I_INUMP (z
))
7088 else if (SCM_BIGP (z
))
7090 else if (SCM_FRACTIONP (z
))
7091 return SCM_FRACTION_NUMERATOR (z
);
7092 else if (SCM_REALP (z
))
7093 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
7095 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
7099 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
7100 /* "Return the denominator of the number @var{z}."
7103 scm_denominator (SCM z
)
7105 if (SCM_I_INUMP (z
))
7107 else if (SCM_BIGP (z
))
7109 else if (SCM_FRACTIONP (z
))
7110 return SCM_FRACTION_DENOMINATOR (z
);
7111 else if (SCM_REALP (z
))
7112 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
7114 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
7117 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
7118 /* "Return the magnitude of the number @var{z}. This is the same as\n"
7119 * "@code{abs} for real arguments, but also allows complex numbers."
7122 scm_magnitude (SCM z
)
7124 if (SCM_I_INUMP (z
))
7126 scm_t_inum zz
= SCM_I_INUM (z
);
7129 else if (SCM_POSFIXABLE (-zz
))
7130 return SCM_I_MAKINUM (-zz
);
7132 return scm_i_inum2big (-zz
);
7134 else if (SCM_BIGP (z
))
7136 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
7137 scm_remember_upto_here_1 (z
);
7139 return scm_i_clonebig (z
, 0);
7143 else if (SCM_REALP (z
))
7144 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
7145 else if (SCM_COMPLEXP (z
))
7146 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
7147 else if (SCM_FRACTIONP (z
))
7149 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
7151 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
7152 SCM_FRACTION_DENOMINATOR (z
));
7155 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
7159 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
7160 /* "Return the angle of the complex number @var{z}."
7165 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
7166 flo0 to save allocating a new flonum with scm_from_double each time.
7167 But if atan2 follows the floating point rounding mode, then the value
7168 is not a constant. Maybe it'd be close enough though. */
7169 if (SCM_I_INUMP (z
))
7171 if (SCM_I_INUM (z
) >= 0)
7174 return scm_from_double (atan2 (0.0, -1.0));
7176 else if (SCM_BIGP (z
))
7178 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
7179 scm_remember_upto_here_1 (z
);
7181 return scm_from_double (atan2 (0.0, -1.0));
7185 else if (SCM_REALP (z
))
7187 if (SCM_REAL_VALUE (z
) >= 0)
7190 return scm_from_double (atan2 (0.0, -1.0));
7192 else if (SCM_COMPLEXP (z
))
7193 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
7194 else if (SCM_FRACTIONP (z
))
7196 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
7198 else return scm_from_double (atan2 (0.0, -1.0));
7201 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
7205 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
7206 /* Convert the number @var{x} to its inexact representation.\n"
7209 scm_exact_to_inexact (SCM z
)
7211 if (SCM_I_INUMP (z
))
7212 return scm_from_double ((double) SCM_I_INUM (z
));
7213 else if (SCM_BIGP (z
))
7214 return scm_from_double (scm_i_big2dbl (z
));
7215 else if (SCM_FRACTIONP (z
))
7216 return scm_from_double (scm_i_fraction2double (z
));
7217 else if (SCM_INEXACTP (z
))
7220 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
7224 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
7226 "Return an exact number that is numerically closest to @var{z}.")
7227 #define FUNC_NAME s_scm_inexact_to_exact
7229 if (SCM_I_INUMP (z
))
7231 else if (SCM_BIGP (z
))
7233 else if (SCM_REALP (z
))
7235 if (isinf (SCM_REAL_VALUE (z
)) || isnan (SCM_REAL_VALUE (z
)))
7236 SCM_OUT_OF_RANGE (1, z
);
7243 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
7244 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
7245 scm_i_mpz2num (mpq_denref (frac
)));
7247 /* When scm_i_make_ratio throws, we leak the memory allocated
7254 else if (SCM_FRACTIONP (z
))
7257 SCM_WRONG_TYPE_ARG (1, z
);
7261 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
7263 "Returns the @emph{simplest} rational number differing\n"
7264 "from @var{x} by no more than @var{eps}.\n"
7266 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
7267 "exact result when both its arguments are exact. Thus, you might need\n"
7268 "to use @code{inexact->exact} on the arguments.\n"
7271 "(rationalize (inexact->exact 1.2) 1/100)\n"
7274 #define FUNC_NAME s_scm_rationalize
7276 if (SCM_I_INUMP (x
))
7278 else if (SCM_BIGP (x
))
7280 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
7282 /* Use continued fractions to find closest ratio. All
7283 arithmetic is done with exact numbers.
7286 SCM ex
= scm_inexact_to_exact (x
);
7287 SCM int_part
= scm_floor (ex
);
7289 SCM a1
= SCM_INUM0
, a2
= SCM_INUM1
, a
= SCM_INUM0
;
7290 SCM b1
= SCM_INUM1
, b2
= SCM_INUM0
, b
= SCM_INUM0
;
7294 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
7297 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
7298 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
7300 /* We stop after a million iterations just to be absolutely sure
7301 that we don't go into an infinite loop. The process normally
7302 converges after less than a dozen iterations.
7305 eps
= scm_abs (eps
);
7306 while (++i
< 1000000)
7308 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
7309 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
7310 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
7312 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
7313 eps
))) /* abs(x-a/b) <= eps */
7315 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
7316 if (scm_is_false (scm_exact_p (x
))
7317 || scm_is_false (scm_exact_p (eps
)))
7318 return scm_exact_to_inexact (res
);
7322 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
7324 tt
= scm_floor (rx
); /* tt = floor (rx) */
7330 scm_num_overflow (s_scm_rationalize
);
7333 SCM_WRONG_TYPE_ARG (1, x
);
7337 /* conversion functions */
7340 scm_is_integer (SCM val
)
7342 return scm_is_true (scm_integer_p (val
));
7346 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
7348 if (SCM_I_INUMP (val
))
7350 scm_t_signed_bits n
= SCM_I_INUM (val
);
7351 return n
>= min
&& n
<= max
;
7353 else if (SCM_BIGP (val
))
7355 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
7357 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
7359 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
7361 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
7362 return n
>= min
&& n
<= max
;
7372 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
7373 > CHAR_BIT
*sizeof (scm_t_uintmax
))
7376 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
7377 SCM_I_BIG_MPZ (val
));
7379 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
7391 return n
>= min
&& n
<= max
;
7399 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
7401 if (SCM_I_INUMP (val
))
7403 scm_t_signed_bits n
= SCM_I_INUM (val
);
7404 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
7406 else if (SCM_BIGP (val
))
7408 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
7410 else if (max
<= ULONG_MAX
)
7412 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
7414 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
7415 return n
>= min
&& n
<= max
;
7425 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
7428 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
7429 > CHAR_BIT
*sizeof (scm_t_uintmax
))
7432 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
7433 SCM_I_BIG_MPZ (val
));
7435 return n
>= min
&& n
<= max
;
7443 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
7445 scm_error (scm_out_of_range_key
,
7447 "Value out of range ~S to ~S: ~S",
7448 scm_list_3 (min
, max
, bad_val
),
7449 scm_list_1 (bad_val
));
7452 #define TYPE scm_t_intmax
7453 #define TYPE_MIN min
7454 #define TYPE_MAX max
7455 #define SIZEOF_TYPE 0
7456 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
7457 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
7458 #include "libguile/conv-integer.i.c"
7460 #define TYPE scm_t_uintmax
7461 #define TYPE_MIN min
7462 #define TYPE_MAX max
7463 #define SIZEOF_TYPE 0
7464 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
7465 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
7466 #include "libguile/conv-uinteger.i.c"
7468 #define TYPE scm_t_int8
7469 #define TYPE_MIN SCM_T_INT8_MIN
7470 #define TYPE_MAX SCM_T_INT8_MAX
7471 #define SIZEOF_TYPE 1
7472 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
7473 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
7474 #include "libguile/conv-integer.i.c"
7476 #define TYPE scm_t_uint8
7478 #define TYPE_MAX SCM_T_UINT8_MAX
7479 #define SIZEOF_TYPE 1
7480 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
7481 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
7482 #include "libguile/conv-uinteger.i.c"
7484 #define TYPE scm_t_int16
7485 #define TYPE_MIN SCM_T_INT16_MIN
7486 #define TYPE_MAX SCM_T_INT16_MAX
7487 #define SIZEOF_TYPE 2
7488 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
7489 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
7490 #include "libguile/conv-integer.i.c"
7492 #define TYPE scm_t_uint16
7494 #define TYPE_MAX SCM_T_UINT16_MAX
7495 #define SIZEOF_TYPE 2
7496 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
7497 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
7498 #include "libguile/conv-uinteger.i.c"
7500 #define TYPE scm_t_int32
7501 #define TYPE_MIN SCM_T_INT32_MIN
7502 #define TYPE_MAX SCM_T_INT32_MAX
7503 #define SIZEOF_TYPE 4
7504 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
7505 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
7506 #include "libguile/conv-integer.i.c"
7508 #define TYPE scm_t_uint32
7510 #define TYPE_MAX SCM_T_UINT32_MAX
7511 #define SIZEOF_TYPE 4
7512 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
7513 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
7514 #include "libguile/conv-uinteger.i.c"
7516 #define TYPE scm_t_wchar
7517 #define TYPE_MIN (scm_t_int32)-1
7518 #define TYPE_MAX (scm_t_int32)0x10ffff
7519 #define SIZEOF_TYPE 4
7520 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
7521 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
7522 #include "libguile/conv-integer.i.c"
7524 #define TYPE scm_t_int64
7525 #define TYPE_MIN SCM_T_INT64_MIN
7526 #define TYPE_MAX SCM_T_INT64_MAX
7527 #define SIZEOF_TYPE 8
7528 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
7529 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
7530 #include "libguile/conv-integer.i.c"
7532 #define TYPE scm_t_uint64
7534 #define TYPE_MAX SCM_T_UINT64_MAX
7535 #define SIZEOF_TYPE 8
7536 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
7537 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
7538 #include "libguile/conv-uinteger.i.c"
7541 scm_to_mpz (SCM val
, mpz_t rop
)
7543 if (SCM_I_INUMP (val
))
7544 mpz_set_si (rop
, SCM_I_INUM (val
));
7545 else if (SCM_BIGP (val
))
7546 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
7548 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
7552 scm_from_mpz (mpz_t val
)
7554 return scm_i_mpz2num (val
);
7558 scm_is_real (SCM val
)
7560 return scm_is_true (scm_real_p (val
));
7564 scm_is_rational (SCM val
)
7566 return scm_is_true (scm_rational_p (val
));
7570 scm_to_double (SCM val
)
7572 if (SCM_I_INUMP (val
))
7573 return SCM_I_INUM (val
);
7574 else if (SCM_BIGP (val
))
7575 return scm_i_big2dbl (val
);
7576 else if (SCM_FRACTIONP (val
))
7577 return scm_i_fraction2double (val
);
7578 else if (SCM_REALP (val
))
7579 return SCM_REAL_VALUE (val
);
7581 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
7585 scm_from_double (double val
)
7589 z
= PTR2SCM (scm_gc_malloc_pointerless (sizeof (scm_t_double
), "real"));
7591 SCM_SET_CELL_TYPE (z
, scm_tc16_real
);
7592 SCM_REAL_VALUE (z
) = val
;
7597 #if SCM_ENABLE_DEPRECATED == 1
7600 scm_num2float (SCM num
, unsigned long pos
, const char *s_caller
)
7602 scm_c_issue_deprecation_warning
7603 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
7607 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
7611 scm_out_of_range (NULL
, num
);
7614 return scm_to_double (num
);
7618 scm_num2double (SCM num
, unsigned long pos
, const char *s_caller
)
7620 scm_c_issue_deprecation_warning
7621 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
7625 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
7629 scm_out_of_range (NULL
, num
);
7632 return scm_to_double (num
);
7638 scm_is_complex (SCM val
)
7640 return scm_is_true (scm_complex_p (val
));
7644 scm_c_real_part (SCM z
)
7646 if (SCM_COMPLEXP (z
))
7647 return SCM_COMPLEX_REAL (z
);
7650 /* Use the scm_real_part to get proper error checking and
7653 return scm_to_double (scm_real_part (z
));
7658 scm_c_imag_part (SCM z
)
7660 if (SCM_COMPLEXP (z
))
7661 return SCM_COMPLEX_IMAG (z
);
7664 /* Use the scm_imag_part to get proper error checking and
7665 dispatching. The result will almost always be 0.0, but not
7668 return scm_to_double (scm_imag_part (z
));
7673 scm_c_magnitude (SCM z
)
7675 return scm_to_double (scm_magnitude (z
));
7681 return scm_to_double (scm_angle (z
));
7685 scm_is_number (SCM z
)
7687 return scm_is_true (scm_number_p (z
));
7691 /* In the following functions we dispatch to the real-arg funcs like log()
7692 when we know the arg is real, instead of just handing everything to
7693 clog() for instance. This is in case clog() doesn't optimize for a
7694 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
7695 well use it to go straight to the applicable C func. */
7697 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
7699 "Return the natural logarithm of @var{z}.")
7700 #define FUNC_NAME s_scm_log
7702 if (SCM_COMPLEXP (z
))
7704 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
7705 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
7707 double re
= SCM_COMPLEX_REAL (z
);
7708 double im
= SCM_COMPLEX_IMAG (z
);
7709 return scm_c_make_rectangular (log (hypot (re
, im
)),
7715 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
7716 although the value itself overflows. */
7717 double re
= scm_to_double (z
);
7718 double l
= log (fabs (re
));
7720 return scm_from_double (l
);
7722 return scm_c_make_rectangular (l
, M_PI
);
7728 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
7730 "Return the base 10 logarithm of @var{z}.")
7731 #define FUNC_NAME s_scm_log10
7733 if (SCM_COMPLEXP (z
))
7735 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
7736 clog() and a multiply by M_LOG10E, rather than the fallback
7737 log10+hypot+atan2.) */
7738 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
7739 && defined SCM_COMPLEX_VALUE
7740 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
7742 double re
= SCM_COMPLEX_REAL (z
);
7743 double im
= SCM_COMPLEX_IMAG (z
);
7744 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
7745 M_LOG10E
* atan2 (im
, re
));
7750 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
7751 although the value itself overflows. */
7752 double re
= scm_to_double (z
);
7753 double l
= log10 (fabs (re
));
7755 return scm_from_double (l
);
7757 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
7763 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
7765 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
7766 "base of natural logarithms (2.71828@dots{}).")
7767 #define FUNC_NAME s_scm_exp
7769 if (SCM_COMPLEXP (z
))
7771 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
7772 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
7774 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
7775 SCM_COMPLEX_IMAG (z
));
7780 /* When z is a negative bignum the conversion to double overflows,
7781 giving -infinity, but that's ok, the exp is still 0.0. */
7782 return scm_from_double (exp (scm_to_double (z
)));
7788 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
7790 "Return the square root of @var{z}. Of the two possible roots\n"
7791 "(positive and negative), the one with the a positive real part\n"
7792 "is returned, or if that's zero then a positive imaginary part.\n"
7796 "(sqrt 9.0) @result{} 3.0\n"
7797 "(sqrt -9.0) @result{} 0.0+3.0i\n"
7798 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
7799 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
7801 #define FUNC_NAME s_scm_sqrt
7803 if (SCM_COMPLEXP (x
))
7805 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
7806 && defined SCM_COMPLEX_VALUE
7807 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
7809 double re
= SCM_COMPLEX_REAL (x
);
7810 double im
= SCM_COMPLEX_IMAG (x
);
7811 return scm_c_make_polar (sqrt (hypot (re
, im
)),
7812 0.5 * atan2 (im
, re
));
7817 double xx
= scm_to_double (x
);
7819 return scm_c_make_rectangular (0.0, sqrt (-xx
));
7821 return scm_from_double (sqrt (xx
));
7833 mpz_init_set_si (z_negative_one
, -1);
7835 /* It may be possible to tune the performance of some algorithms by using
7836 * the following constants to avoid the creation of bignums. Please, before
7837 * using these values, remember the two rules of program optimization:
7838 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
7839 scm_c_define ("most-positive-fixnum",
7840 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
7841 scm_c_define ("most-negative-fixnum",
7842 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
7844 scm_add_feature ("complex");
7845 scm_add_feature ("inexact");
7846 flo0
= scm_from_double (0.0);
7848 /* determine floating point precision */
7849 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
7851 init_dblprec(&scm_dblprec
[i
-2],i
);
7852 init_fx_radix(fx_per_radix
[i
-2],i
);
7855 /* hard code precision for base 10 if the preprocessor tells us to... */
7856 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
7859 exactly_one_half
= scm_divide (SCM_INUM1
, SCM_I_MAKINUM (2));
7860 #include "libguile/numbers.x"