1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008, 2009, 2010 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
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful,
13 * but 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 02110-1301 USA
23 /* General assumptions:
24 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
25 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
26 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
27 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
28 * All objects satisfying SCM_FRACTIONP are never an integer.
33 - see if special casing bignums and reals in integer-exponent when
34 possible (to use mpz_pow and mpf_pow_ui) is faster.
36 - look in to better short-circuiting of common cases in
37 integer-expt and elsewhere.
39 - see if direct mpz operations can help in ash and elsewhere.
55 #include "libguile/_scm.h"
56 #include "libguile/feature.h"
57 #include "libguile/ports.h"
58 #include "libguile/root.h"
59 #include "libguile/smob.h"
60 #include "libguile/strings.h"
62 #include "libguile/validate.h"
63 #include "libguile/numbers.h"
64 #include "libguile/deprecation.h"
66 #include "libguile/eq.h"
68 #include "libguile/discouraged.h"
70 /* values per glibc, if not already defined */
72 #define M_LOG10E 0.43429448190325182765
75 #define M_PI 3.14159265358979323846
81 Wonder if this might be faster for some of our code? A switch on
82 the numtag would jump directly to the right case, and the
83 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
85 #define SCM_I_NUMTAG_NOTNUM 0
86 #define SCM_I_NUMTAG_INUM 1
87 #define SCM_I_NUMTAG_BIG scm_tc16_big
88 #define SCM_I_NUMTAG_REAL scm_tc16_real
89 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
90 #define SCM_I_NUMTAG(x) \
91 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
92 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
93 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
94 : SCM_I_NUMTAG_NOTNUM)))
96 /* the macro above will not work as is with fractions */
99 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
101 /* FLOBUFLEN is the maximum number of characters neccessary for the
102 * printed or scm_string representation of an inexact number.
104 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
107 #if ! defined (HAVE_ISNAN)
112 return (IsNANorINF (x) && NaN (x) && ! IsINF (x)) ? 1 : 0;
115 #if ! defined (HAVE_ISINF)
120 return (IsNANorINF (x) && IsINF (x)) ? 1 : 0;
127 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
128 an explicit check. In some future gmp (don't know what version number),
129 mpz_cmp_d is supposed to do this itself. */
131 #define xmpz_cmp_d(z, d) \
132 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
134 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
137 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
138 isinf. It does have finite and isnan though, hence the use of those.
139 fpclass would be a possibility on that system too. */
143 #if defined (HAVE_ISINF)
145 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
146 return (! (finite (x) || isnan (x)));
155 #if defined (HAVE_ISNAN)
162 #if defined (GUILE_I)
163 #if HAVE_COMPLEX_DOUBLE
165 /* For an SCM object Z which is a complex number (ie. satisfies
166 SCM_COMPLEXP), return its value as a C level "complex double". */
167 #define SCM_COMPLEX_VALUE(z) \
168 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
170 static inline SCM scm_from_complex_double (complex double z) SCM_UNUSED;
172 /* Convert a C "complex double" to an SCM value. */
174 scm_from_complex_double (complex double z)
176 return scm_c_make_rectangular (creal (z), cimag (z));
179 #endif /* HAVE_COMPLEX_DOUBLE */
184 static mpz_t z_negative_one;
191 /* Return a newly created bignum. */
192 SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
193 mpz_init (SCM_I_BIG_MPZ (z));
198 scm_i_long2big (long x)
200 /* Return a newly created bignum initialized to X. */
201 SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
202 mpz_init_set_si (SCM_I_BIG_MPZ (z), x);
207 scm_i_ulong2big (unsigned long x)
209 /* Return a newly created bignum initialized to X. */
210 SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
211 mpz_init_set_ui (SCM_I_BIG_MPZ (z), x);
216 scm_i_clonebig (SCM src_big, int same_sign_p)
218 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
219 SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
220 mpz_init_set (SCM_I_BIG_MPZ (z), SCM_I_BIG_MPZ (src_big));
222 mpz_neg (SCM_I_BIG_MPZ (z), SCM_I_BIG_MPZ (z));
227 scm_i_bigcmp (SCM x, SCM y)
229 /* Return neg if x < y, pos if x > y, and 0 if x == y */
230 /* presume we already know x and y are bignums */
231 int result = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
232 scm_remember_upto_here_2 (x, y);
237 scm_i_dbl2big (double d)
239 /* results are only defined if d is an integer */
240 SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
241 mpz_init_set_d (SCM_I_BIG_MPZ (z), d);
245 /* Convert a integer in double representation to a SCM number. */
248 scm_i_dbl2num (double u)
250 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
251 powers of 2, so there's no rounding when making "double" values
252 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
253 get rounded on a 64-bit machine, hence the "+1".
255 The use of floor() to force to an integer value ensures we get a
256 "numerically closest" value without depending on how a
257 double->long cast or how mpz_set_d will round. For reference,
258 double->long probably follows the hardware rounding mode,
259 mpz_set_d truncates towards zero. */
261 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
262 representable as a double? */
264 if (u < (double) (SCM_MOST_POSITIVE_FIXNUM+1)
265 && u >= (double) SCM_MOST_NEGATIVE_FIXNUM)
266 return SCM_I_MAKINUM ((long) u);
268 return scm_i_dbl2big (u);
271 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
272 with R5RS exact->inexact.
274 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
275 (ie. truncate towards zero), then adjust to get the closest double by
276 examining the next lower bit and adding 1 (to the absolute value) if
279 Bignums exactly half way between representable doubles are rounded to the
280 next higher absolute value (ie. away from zero). This seems like an
281 adequate interpretation of R5RS "numerically closest", and it's easier
282 and faster than a full "nearest-even" style.
284 The bit test must be done on the absolute value of the mpz_t, which means
285 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
286 negatives as twos complement.
288 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
289 following the hardware rounding mode, but applied to the absolute value
290 of the mpz_t operand. This is not what we want so we put the high
291 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
292 mpz_get_d is supposed to always truncate towards zero.
294 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
295 is a slowdown. It'd be faster to pick out the relevant high bits with
296 mpz_getlimbn if we could be bothered coding that, and if the new
297 truncating gmp doesn't come out. */
300 scm_i_big2dbl (SCM b)
305 bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2);
309 /* Current GMP, eg. 4.1.3, force truncation towards zero */
311 if (bits > DBL_MANT_DIG)
313 size_t shift = bits - DBL_MANT_DIG;
314 mpz_init2 (tmp, DBL_MANT_DIG);
315 mpz_tdiv_q_2exp (tmp, SCM_I_BIG_MPZ (b), shift);
316 result = ldexp (mpz_get_d (tmp), shift);
321 result = mpz_get_d (SCM_I_BIG_MPZ (b));
326 result = mpz_get_d (SCM_I_BIG_MPZ (b));
329 if (bits > DBL_MANT_DIG)
331 unsigned long pos = bits - DBL_MANT_DIG - 1;
332 /* test bit number "pos" in absolute value */
333 if (mpz_getlimbn (SCM_I_BIG_MPZ (b), pos / GMP_NUMB_BITS)
334 & ((mp_limb_t) 1 << (pos % GMP_NUMB_BITS)))
336 result += ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b)), pos + 1);
340 scm_remember_upto_here_1 (b);
345 scm_i_normbig (SCM b)
347 /* convert a big back to a fixnum if it'll fit */
348 /* presume b is a bignum */
349 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b)))
351 long val = mpz_get_si (SCM_I_BIG_MPZ (b));
352 if (SCM_FIXABLE (val))
353 b = SCM_I_MAKINUM (val);
358 static SCM_C_INLINE_KEYWORD SCM
359 scm_i_mpz2num (mpz_t b)
361 /* convert a mpz number to a SCM number. */
362 if (mpz_fits_slong_p (b))
364 long val = mpz_get_si (b);
365 if (SCM_FIXABLE (val))
366 return SCM_I_MAKINUM (val);
370 SCM z = scm_double_cell (scm_tc16_big, 0, 0, 0);
371 mpz_init_set (SCM_I_BIG_MPZ (z), b);
376 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
377 static SCM scm_divide2real (SCM x, SCM y);
380 scm_i_make_ratio (SCM numerator, SCM denominator)
381 #define FUNC_NAME "make-ratio"
383 /* First make sure the arguments are proper.
385 if (SCM_I_INUMP (denominator))
387 if (scm_is_eq (denominator, SCM_INUM0))
388 scm_num_overflow ("make-ratio");
389 if (scm_is_eq (denominator, SCM_I_MAKINUM(1)))
394 if (!(SCM_BIGP(denominator)))
395 SCM_WRONG_TYPE_ARG (2, denominator);
397 if (!SCM_I_INUMP (numerator) && !SCM_BIGP (numerator))
398 SCM_WRONG_TYPE_ARG (1, numerator);
400 /* Then flip signs so that the denominator is positive.
402 if (scm_is_true (scm_negative_p (denominator)))
404 numerator = scm_difference (numerator, SCM_UNDEFINED);
405 denominator = scm_difference (denominator, SCM_UNDEFINED);
408 /* Now consider for each of the four fixnum/bignum combinations
409 whether the rational number is really an integer.
411 if (SCM_I_INUMP (numerator))
413 long x = SCM_I_INUM (numerator);
414 if (scm_is_eq (numerator, SCM_INUM0))
416 if (SCM_I_INUMP (denominator))
419 y = SCM_I_INUM (denominator);
421 return SCM_I_MAKINUM(1);
423 return SCM_I_MAKINUM (x / y);
427 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
428 of that value for the denominator, as a bignum. Apart from
429 that case, abs(bignum) > abs(inum) so inum/bignum is not an
431 if (x == SCM_MOST_NEGATIVE_FIXNUM
432 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator),
433 - SCM_MOST_NEGATIVE_FIXNUM) == 0)
434 return SCM_I_MAKINUM(-1);
437 else if (SCM_BIGP (numerator))
439 if (SCM_I_INUMP (denominator))
441 long yy = SCM_I_INUM (denominator);
442 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator), yy))
443 return scm_divide (numerator, denominator);
447 if (scm_is_eq (numerator, denominator))
448 return SCM_I_MAKINUM(1);
449 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator),
450 SCM_I_BIG_MPZ (denominator)))
451 return scm_divide(numerator, denominator);
455 /* No, it's a proper fraction.
458 SCM divisor = scm_gcd (numerator, denominator);
459 if (!(scm_is_eq (divisor, SCM_I_MAKINUM(1))))
461 numerator = scm_divide (numerator, divisor);
462 denominator = scm_divide (denominator, divisor);
465 return scm_double_cell (scm_tc16_fraction,
466 SCM_UNPACK (numerator),
467 SCM_UNPACK (denominator), 0);
473 scm_i_fraction2double (SCM z)
475 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z),
476 SCM_FRACTION_DENOMINATOR (z)));
479 SCM_DEFINE (scm_exact_p, "exact?", 1, 0, 0,
481 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
483 #define FUNC_NAME s_scm_exact_p
489 if (SCM_FRACTIONP (x))
493 SCM_WRONG_TYPE_ARG (1, x);
498 SCM_DEFINE (scm_odd_p, "odd?", 1, 0, 0,
500 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
502 #define FUNC_NAME s_scm_odd_p
506 long val = SCM_I_INUM (n);
507 return scm_from_bool ((val & 1L) != 0);
509 else if (SCM_BIGP (n))
511 int odd_p = mpz_odd_p (SCM_I_BIG_MPZ (n));
512 scm_remember_upto_here_1 (n);
513 return scm_from_bool (odd_p);
515 else if (scm_is_true (scm_inf_p (n)))
517 else if (SCM_REALP (n))
519 double rem = fabs (fmod (SCM_REAL_VALUE(n), 2.0));
525 SCM_WRONG_TYPE_ARG (1, n);
528 SCM_WRONG_TYPE_ARG (1, n);
533 SCM_DEFINE (scm_even_p, "even?", 1, 0, 0,
535 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
537 #define FUNC_NAME s_scm_even_p
541 long val = SCM_I_INUM (n);
542 return scm_from_bool ((val & 1L) == 0);
544 else if (SCM_BIGP (n))
546 int even_p = mpz_even_p (SCM_I_BIG_MPZ (n));
547 scm_remember_upto_here_1 (n);
548 return scm_from_bool (even_p);
550 else if (scm_is_true (scm_inf_p (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);
567 SCM_DEFINE (scm_inf_p, "inf?", 1, 0, 0,
569 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
570 "or @samp{-inf.0}, @code{#f} otherwise.")
571 #define FUNC_NAME s_scm_inf_p
574 return scm_from_bool (xisinf (SCM_REAL_VALUE (x)));
575 else if (SCM_COMPLEXP (x))
576 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x))
577 || xisinf (SCM_COMPLEX_IMAG (x)));
583 SCM_DEFINE (scm_nan_p, "nan?", 1, 0, 0,
585 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
587 #define FUNC_NAME s_scm_nan_p
590 return scm_from_bool (xisnan (SCM_REAL_VALUE (n)));
591 else if (SCM_COMPLEXP (n))
592 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n))
593 || xisnan (SCM_COMPLEX_IMAG (n)));
599 /* Guile's idea of infinity. */
600 static double guile_Inf;
602 /* Guile's idea of not a number. */
603 static double guile_NaN;
606 guile_ieee_init (void)
608 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
610 /* Some version of gcc on some old version of Linux used to crash when
611 trying to make Inf and NaN. */
614 /* C99 INFINITY, when available.
615 FIXME: The standard allows for INFINITY to be something that overflows
616 at compile time. We ought to have a configure test to check for that
617 before trying to use it. (But in practice we believe this is not a
618 problem on any system guile is likely to target.) */
619 guile_Inf = INFINITY;
622 extern unsigned int DINFINITY[2];
628 alias.i[0] = DINFINITY[0];
629 alias.i[1] = DINFINITY[1];
637 if (guile_Inf == tmp)
645 #if defined (HAVE_ISNAN)
647 #if defined __GNUC__ && defined __alpha__ && !defined _IEEE_FP
648 /* On Alpha GCC must be passed `-mieee' to provide proper NaN handling.
649 See http://lists.gnu.org/archive/html/bug-gnulib/2009-05/msg00010.html
651 # error NaN handling will not work when compiling without -mieee
655 /* C99 NAN, when available */
660 extern unsigned int DQNAN[2];
666 alias.i[0] = DQNAN[0];
667 alias.i[1] = DQNAN[1];
671 guile_NaN = guile_Inf / guile_Inf;
677 SCM_DEFINE (scm_inf, "inf", 0, 0, 0,
680 #define FUNC_NAME s_scm_inf
682 static int initialized = 0;
688 return scm_from_double (guile_Inf);
692 SCM_DEFINE (scm_nan, "nan", 0, 0, 0,
695 #define FUNC_NAME s_scm_nan
697 static int initialized = 0;
703 return scm_from_double (guile_NaN);
708 SCM_PRIMITIVE_GENERIC (scm_abs, "abs", 1, 0, 0,
710 "Return the absolute value of @var{x}.")
715 long int xx = SCM_I_INUM (x);
718 else if (SCM_POSFIXABLE (-xx))
719 return SCM_I_MAKINUM (-xx);
721 return scm_i_long2big (-xx);
723 else if (SCM_BIGP (x))
725 const int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
727 return scm_i_clonebig (x, 0);
731 else if (SCM_REALP (x))
733 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
734 double xx = SCM_REAL_VALUE (x);
736 return scm_from_double (-xx);
740 else if (SCM_FRACTIONP (x))
742 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x))))
744 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED),
745 SCM_FRACTION_DENOMINATOR (x));
748 SCM_WTA_DISPATCH_1 (g_scm_abs, x, 1, s_scm_abs);
753 SCM_GPROC (s_quotient, "quotient", 2, 0, 0, scm_quotient, g_quotient);
754 /* "Return the quotient of the numbers @var{x} and @var{y}."
757 scm_quotient (SCM x, SCM y)
761 long xx = SCM_I_INUM (x);
764 long yy = SCM_I_INUM (y);
766 scm_num_overflow (s_quotient);
771 return SCM_I_MAKINUM (z);
773 return scm_i_long2big (z);
776 else if (SCM_BIGP (y))
778 if ((SCM_I_INUM (x) == SCM_MOST_NEGATIVE_FIXNUM)
779 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y),
780 - SCM_MOST_NEGATIVE_FIXNUM) == 0))
782 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
783 scm_remember_upto_here_1 (y);
784 return SCM_I_MAKINUM (-1);
787 return SCM_I_MAKINUM (0);
790 SCM_WTA_DISPATCH_2 (g_quotient, x, y, SCM_ARG2, s_quotient);
792 else if (SCM_BIGP (x))
796 long yy = SCM_I_INUM (y);
798 scm_num_overflow (s_quotient);
803 SCM result = scm_i_mkbig ();
806 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result),
809 mpz_neg (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result));
812 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), yy);
813 scm_remember_upto_here_1 (x);
814 return scm_i_normbig (result);
817 else if (SCM_BIGP (y))
819 SCM result = scm_i_mkbig ();
820 mpz_tdiv_q (SCM_I_BIG_MPZ (result),
823 scm_remember_upto_here_2 (x, y);
824 return scm_i_normbig (result);
827 SCM_WTA_DISPATCH_2 (g_quotient, x, y, SCM_ARG2, s_quotient);
830 SCM_WTA_DISPATCH_2 (g_quotient, x, y, SCM_ARG1, s_quotient);
833 SCM_GPROC (s_remainder, "remainder", 2, 0, 0, scm_remainder, g_remainder);
834 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
836 * "(remainder 13 4) @result{} 1\n"
837 * "(remainder -13 4) @result{} -1\n"
841 scm_remainder (SCM x, SCM y)
847 long yy = SCM_I_INUM (y);
849 scm_num_overflow (s_remainder);
852 long z = SCM_I_INUM (x) % yy;
853 return SCM_I_MAKINUM (z);
856 else if (SCM_BIGP (y))
858 if ((SCM_I_INUM (x) == SCM_MOST_NEGATIVE_FIXNUM)
859 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y),
860 - SCM_MOST_NEGATIVE_FIXNUM) == 0))
862 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
863 scm_remember_upto_here_1 (y);
864 return SCM_I_MAKINUM (0);
870 SCM_WTA_DISPATCH_2 (g_remainder, x, y, SCM_ARG2, s_remainder);
872 else if (SCM_BIGP (x))
876 long yy = SCM_I_INUM (y);
878 scm_num_overflow (s_remainder);
881 SCM result = scm_i_mkbig ();
884 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ(x), yy);
885 scm_remember_upto_here_1 (x);
886 return scm_i_normbig (result);
889 else if (SCM_BIGP (y))
891 SCM result = scm_i_mkbig ();
892 mpz_tdiv_r (SCM_I_BIG_MPZ (result),
895 scm_remember_upto_here_2 (x, y);
896 return scm_i_normbig (result);
899 SCM_WTA_DISPATCH_2 (g_remainder, x, y, SCM_ARG2, s_remainder);
902 SCM_WTA_DISPATCH_2 (g_remainder, x, y, SCM_ARG1, s_remainder);
906 SCM_GPROC (s_modulo, "modulo", 2, 0, 0, scm_modulo, g_modulo);
907 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
909 * "(modulo 13 4) @result{} 1\n"
910 * "(modulo -13 4) @result{} 3\n"
914 scm_modulo (SCM x, SCM y)
918 long xx = SCM_I_INUM (x);
921 long yy = SCM_I_INUM (y);
923 scm_num_overflow (s_modulo);
926 /* C99 specifies that "%" is the remainder corresponding to a
927 quotient rounded towards zero, and that's also traditional
928 for machine division, so z here should be well defined. */
946 return SCM_I_MAKINUM (result);
949 else if (SCM_BIGP (y))
951 int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y));
958 SCM pos_y = scm_i_clonebig (y, 0);
959 /* do this after the last scm_op */
960 mpz_init_set_si (z_x, xx);
961 result = pos_y; /* re-use this bignum */
962 mpz_mod (SCM_I_BIG_MPZ (result),
964 SCM_I_BIG_MPZ (pos_y));
965 scm_remember_upto_here_1 (pos_y);
969 result = scm_i_mkbig ();
970 /* do this after the last scm_op */
971 mpz_init_set_si (z_x, xx);
972 mpz_mod (SCM_I_BIG_MPZ (result),
975 scm_remember_upto_here_1 (y);
978 if ((sgn_y < 0) && mpz_sgn (SCM_I_BIG_MPZ (result)) != 0)
979 mpz_add (SCM_I_BIG_MPZ (result),
981 SCM_I_BIG_MPZ (result));
982 scm_remember_upto_here_1 (y);
983 /* and do this before the next one */
985 return scm_i_normbig (result);
989 SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG2, s_modulo);
991 else if (SCM_BIGP (x))
995 long yy = SCM_I_INUM (y);
997 scm_num_overflow (s_modulo);
1000 SCM result = scm_i_mkbig ();
1001 mpz_mod_ui (SCM_I_BIG_MPZ (result),
1003 (yy < 0) ? - yy : yy);
1004 scm_remember_upto_here_1 (x);
1005 if ((yy < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0))
1006 mpz_sub_ui (SCM_I_BIG_MPZ (result),
1007 SCM_I_BIG_MPZ (result),
1009 return scm_i_normbig (result);
1012 else if (SCM_BIGP (y))
1015 SCM result = scm_i_mkbig ();
1016 int y_sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
1017 SCM pos_y = scm_i_clonebig (y, y_sgn >= 0);
1018 mpz_mod (SCM_I_BIG_MPZ (result),
1020 SCM_I_BIG_MPZ (pos_y));
1022 scm_remember_upto_here_1 (x);
1023 if ((y_sgn < 0) && (mpz_sgn (SCM_I_BIG_MPZ (result)) != 0))
1024 mpz_add (SCM_I_BIG_MPZ (result),
1026 SCM_I_BIG_MPZ (result));
1027 scm_remember_upto_here_2 (y, pos_y);
1028 return scm_i_normbig (result);
1032 SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG2, s_modulo);
1035 SCM_WTA_DISPATCH_2 (g_modulo, x, y, SCM_ARG1, s_modulo);
1038 SCM_GPROC1 (s_gcd, "gcd", scm_tc7_asubr, scm_gcd, g_gcd);
1039 /* "Return the greatest common divisor of all arguments.\n"
1040 * "If called without arguments, 0 is returned."
1043 scm_gcd (SCM x, SCM y)
1046 return SCM_UNBNDP (x) ? SCM_INUM0 : scm_abs (x);
1048 if (SCM_I_INUMP (x))
1050 if (SCM_I_INUMP (y))
1052 long xx = SCM_I_INUM (x);
1053 long yy = SCM_I_INUM (y);
1054 long u = xx < 0 ? -xx : xx;
1055 long v = yy < 0 ? -yy : yy;
1065 /* Determine a common factor 2^k */
1066 while (!(1 & (u | v)))
1072 /* Now, any factor 2^n can be eliminated */
1092 return (SCM_POSFIXABLE (result)
1093 ? SCM_I_MAKINUM (result)
1094 : scm_i_long2big (result));
1096 else if (SCM_BIGP (y))
1102 SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG2, s_gcd);
1104 else if (SCM_BIGP (x))
1106 if (SCM_I_INUMP (y))
1108 unsigned long result;
1111 yy = SCM_I_INUM (y);
1116 result = mpz_gcd_ui (NULL, SCM_I_BIG_MPZ (x), yy);
1117 scm_remember_upto_here_1 (x);
1118 return (SCM_POSFIXABLE (result)
1119 ? SCM_I_MAKINUM (result)
1120 : scm_from_ulong (result));
1122 else if (SCM_BIGP (y))
1124 SCM result = scm_i_mkbig ();
1125 mpz_gcd (SCM_I_BIG_MPZ (result),
1128 scm_remember_upto_here_2 (x, y);
1129 return scm_i_normbig (result);
1132 SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG2, s_gcd);
1135 SCM_WTA_DISPATCH_2 (g_gcd, x, y, SCM_ARG1, s_gcd);
1138 SCM_GPROC1 (s_lcm, "lcm", scm_tc7_asubr, scm_lcm, g_lcm);
1139 /* "Return the least common multiple of the arguments.\n"
1140 * "If called without arguments, 1 is returned."
1143 scm_lcm (SCM n1, SCM n2)
1145 if (SCM_UNBNDP (n2))
1147 if (SCM_UNBNDP (n1))
1148 return SCM_I_MAKINUM (1L);
1149 n2 = SCM_I_MAKINUM (1L);
1152 SCM_GASSERT2 (SCM_I_INUMP (n1) || SCM_BIGP (n1),
1153 g_lcm, n1, n2, SCM_ARG1, s_lcm);
1154 SCM_GASSERT2 (SCM_I_INUMP (n2) || SCM_BIGP (n2),
1155 g_lcm, n1, n2, SCM_ARGn, s_lcm);
1157 if (SCM_I_INUMP (n1))
1159 if (SCM_I_INUMP (n2))
1161 SCM d = scm_gcd (n1, n2);
1162 if (scm_is_eq (d, SCM_INUM0))
1165 return scm_abs (scm_product (n1, scm_quotient (n2, d)));
1169 /* inum n1, big n2 */
1172 SCM result = scm_i_mkbig ();
1173 long nn1 = SCM_I_INUM (n1);
1174 if (nn1 == 0) return SCM_INUM0;
1175 if (nn1 < 0) nn1 = - nn1;
1176 mpz_lcm_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n2), nn1);
1177 scm_remember_upto_here_1 (n2);
1185 if (SCM_I_INUMP (n2))
1192 SCM result = scm_i_mkbig ();
1193 mpz_lcm(SCM_I_BIG_MPZ (result),
1195 SCM_I_BIG_MPZ (n2));
1196 scm_remember_upto_here_2(n1, n2);
1197 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1203 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1208 + + + x (map digit:logand X Y)
1209 + - + x (map digit:logand X (lognot (+ -1 Y)))
1210 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1211 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1216 + + + (map digit:logior X Y)
1217 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1218 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1219 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1224 + + + (map digit:logxor X Y)
1225 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1226 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1227 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1232 + + (any digit:logand X Y)
1233 + - (any digit:logand X (lognot (+ -1 Y)))
1234 - + (any digit:logand (lognot (+ -1 X)) Y)
1239 SCM_DEFINE1 (scm_logand, "logand", scm_tc7_asubr,
1241 "Return the bitwise AND of the integer arguments.\n\n"
1243 "(logand) @result{} -1\n"
1244 "(logand 7) @result{} 7\n"
1245 "(logand #b111 #b011 #b001) @result{} 1\n"
1247 #define FUNC_NAME s_scm_logand
1251 if (SCM_UNBNDP (n2))
1253 if (SCM_UNBNDP (n1))
1254 return SCM_I_MAKINUM (-1);
1255 else if (!SCM_NUMBERP (n1))
1256 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1257 else if (SCM_NUMBERP (n1))
1260 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1263 if (SCM_I_INUMP (n1))
1265 nn1 = SCM_I_INUM (n1);
1266 if (SCM_I_INUMP (n2))
1268 long nn2 = SCM_I_INUM (n2);
1269 return SCM_I_MAKINUM (nn1 & nn2);
1271 else if SCM_BIGP (n2)
1277 SCM result_z = scm_i_mkbig ();
1279 mpz_init_set_si (nn1_z, nn1);
1280 mpz_and (SCM_I_BIG_MPZ (result_z), nn1_z, SCM_I_BIG_MPZ (n2));
1281 scm_remember_upto_here_1 (n2);
1283 return scm_i_normbig (result_z);
1287 SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
1289 else if (SCM_BIGP (n1))
1291 if (SCM_I_INUMP (n2))
1294 nn1 = SCM_I_INUM (n1);
1297 else if (SCM_BIGP (n2))
1299 SCM result_z = scm_i_mkbig ();
1300 mpz_and (SCM_I_BIG_MPZ (result_z),
1302 SCM_I_BIG_MPZ (n2));
1303 scm_remember_upto_here_2 (n1, n2);
1304 return scm_i_normbig (result_z);
1307 SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
1310 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1315 SCM_DEFINE1 (scm_logior, "logior", scm_tc7_asubr,
1317 "Return the bitwise OR of the integer arguments.\n\n"
1319 "(logior) @result{} 0\n"
1320 "(logior 7) @result{} 7\n"
1321 "(logior #b000 #b001 #b011) @result{} 3\n"
1323 #define FUNC_NAME s_scm_logior
1327 if (SCM_UNBNDP (n2))
1329 if (SCM_UNBNDP (n1))
1331 else if (SCM_NUMBERP (n1))
1334 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1337 if (SCM_I_INUMP (n1))
1339 nn1 = SCM_I_INUM (n1);
1340 if (SCM_I_INUMP (n2))
1342 long nn2 = SCM_I_INUM (n2);
1343 return SCM_I_MAKINUM (nn1 | nn2);
1345 else if (SCM_BIGP (n2))
1351 SCM result_z = scm_i_mkbig ();
1353 mpz_init_set_si (nn1_z, nn1);
1354 mpz_ior (SCM_I_BIG_MPZ (result_z), nn1_z, SCM_I_BIG_MPZ (n2));
1355 scm_remember_upto_here_1 (n2);
1357 return scm_i_normbig (result_z);
1361 SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
1363 else if (SCM_BIGP (n1))
1365 if (SCM_I_INUMP (n2))
1368 nn1 = SCM_I_INUM (n1);
1371 else if (SCM_BIGP (n2))
1373 SCM result_z = scm_i_mkbig ();
1374 mpz_ior (SCM_I_BIG_MPZ (result_z),
1376 SCM_I_BIG_MPZ (n2));
1377 scm_remember_upto_here_2 (n1, n2);
1378 return scm_i_normbig (result_z);
1381 SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
1384 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1389 SCM_DEFINE1 (scm_logxor, "logxor", scm_tc7_asubr,
1391 "Return the bitwise XOR of the integer arguments. A bit is\n"
1392 "set in the result if it is set in an odd number of arguments.\n"
1394 "(logxor) @result{} 0\n"
1395 "(logxor 7) @result{} 7\n"
1396 "(logxor #b000 #b001 #b011) @result{} 2\n"
1397 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1399 #define FUNC_NAME s_scm_logxor
1403 if (SCM_UNBNDP (n2))
1405 if (SCM_UNBNDP (n1))
1407 else if (SCM_NUMBERP (n1))
1410 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1413 if (SCM_I_INUMP (n1))
1415 nn1 = SCM_I_INUM (n1);
1416 if (SCM_I_INUMP (n2))
1418 long nn2 = SCM_I_INUM (n2);
1419 return SCM_I_MAKINUM (nn1 ^ nn2);
1421 else if (SCM_BIGP (n2))
1425 SCM result_z = scm_i_mkbig ();
1427 mpz_init_set_si (nn1_z, nn1);
1428 mpz_xor (SCM_I_BIG_MPZ (result_z), nn1_z, SCM_I_BIG_MPZ (n2));
1429 scm_remember_upto_here_1 (n2);
1431 return scm_i_normbig (result_z);
1435 SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
1437 else if (SCM_BIGP (n1))
1439 if (SCM_I_INUMP (n2))
1442 nn1 = SCM_I_INUM (n1);
1445 else if (SCM_BIGP (n2))
1447 SCM result_z = scm_i_mkbig ();
1448 mpz_xor (SCM_I_BIG_MPZ (result_z),
1450 SCM_I_BIG_MPZ (n2));
1451 scm_remember_upto_here_2 (n1, n2);
1452 return scm_i_normbig (result_z);
1455 SCM_WRONG_TYPE_ARG (SCM_ARG2, n2);
1458 SCM_WRONG_TYPE_ARG (SCM_ARG1, n1);
1463 SCM_DEFINE (scm_logtest, "logtest", 2, 0, 0,
1465 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1466 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1467 "without actually calculating the @code{logand}, just testing\n"
1471 "(logtest #b0100 #b1011) @result{} #f\n"
1472 "(logtest #b0100 #b0111) @result{} #t\n"
1474 #define FUNC_NAME s_scm_logtest
1478 if (SCM_I_INUMP (j))
1480 nj = SCM_I_INUM (j);
1481 if (SCM_I_INUMP (k))
1483 long nk = SCM_I_INUM (k);
1484 return scm_from_bool (nj & nk);
1486 else if (SCM_BIGP (k))
1494 mpz_init_set_si (nj_z, nj);
1495 mpz_and (nj_z, nj_z, SCM_I_BIG_MPZ (k));
1496 scm_remember_upto_here_1 (k);
1497 result = scm_from_bool (mpz_sgn (nj_z) != 0);
1503 SCM_WRONG_TYPE_ARG (SCM_ARG2, k);
1505 else if (SCM_BIGP (j))
1507 if (SCM_I_INUMP (k))
1510 nj = SCM_I_INUM (j);
1513 else if (SCM_BIGP (k))
1517 mpz_init (result_z);
1521 scm_remember_upto_here_2 (j, k);
1522 result = scm_from_bool (mpz_sgn (result_z) != 0);
1523 mpz_clear (result_z);
1527 SCM_WRONG_TYPE_ARG (SCM_ARG2, k);
1530 SCM_WRONG_TYPE_ARG (SCM_ARG1, j);
1535 SCM_DEFINE (scm_logbit_p, "logbit?", 2, 0, 0,
1537 "Test whether bit number @var{index} in @var{j} is set.\n"
1538 "@var{index} starts from 0 for the least significant bit.\n"
1541 "(logbit? 0 #b1101) @result{} #t\n"
1542 "(logbit? 1 #b1101) @result{} #f\n"
1543 "(logbit? 2 #b1101) @result{} #t\n"
1544 "(logbit? 3 #b1101) @result{} #t\n"
1545 "(logbit? 4 #b1101) @result{} #f\n"
1547 #define FUNC_NAME s_scm_logbit_p
1549 unsigned long int iindex;
1550 iindex = scm_to_ulong (index);
1552 if (SCM_I_INUMP (j))
1554 /* bits above what's in an inum follow the sign bit */
1555 iindex = min (iindex, SCM_LONG_BIT - 1);
1556 return scm_from_bool ((1L << iindex) & SCM_I_INUM (j));
1558 else if (SCM_BIGP (j))
1560 int val = mpz_tstbit (SCM_I_BIG_MPZ (j), iindex);
1561 scm_remember_upto_here_1 (j);
1562 return scm_from_bool (val);
1565 SCM_WRONG_TYPE_ARG (SCM_ARG2, j);
1570 SCM_DEFINE (scm_lognot, "lognot", 1, 0, 0,
1572 "Return the integer which is the ones-complement of the integer\n"
1576 "(number->string (lognot #b10000000) 2)\n"
1577 " @result{} \"-10000001\"\n"
1578 "(number->string (lognot #b0) 2)\n"
1579 " @result{} \"-1\"\n"
1581 #define FUNC_NAME s_scm_lognot
1583 if (SCM_I_INUMP (n)) {
1584 /* No overflow here, just need to toggle all the bits making up the inum.
1585 Enhancement: No need to strip the tag and add it back, could just xor
1586 a block of 1 bits, if that worked with the various debug versions of
1588 return SCM_I_MAKINUM (~ SCM_I_INUM (n));
1590 } else if (SCM_BIGP (n)) {
1591 SCM result = scm_i_mkbig ();
1592 mpz_com (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n));
1593 scm_remember_upto_here_1 (n);
1597 SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
1602 /* returns 0 if IN is not an integer. OUT must already be
1605 coerce_to_big (SCM in, mpz_t out)
1608 mpz_set (out, SCM_I_BIG_MPZ (in));
1609 else if (SCM_I_INUMP (in))
1610 mpz_set_si (out, SCM_I_INUM (in));
1617 SCM_DEFINE (scm_modulo_expt, "modulo-expt", 3, 0, 0,
1618 (SCM n, SCM k, SCM m),
1619 "Return @var{n} raised to the integer exponent\n"
1620 "@var{k}, modulo @var{m}.\n"
1623 "(modulo-expt 2 3 5)\n"
1626 #define FUNC_NAME s_scm_modulo_expt
1632 /* There are two classes of error we might encounter --
1633 1) Math errors, which we'll report by calling scm_num_overflow,
1635 2) wrong-type errors, which of course we'll report by calling
1637 We don't report those errors immediately, however; instead we do
1638 some cleanup first. These variables tell us which error (if
1639 any) we should report after cleaning up.
1641 int report_overflow = 0;
1643 int position_of_wrong_type = 0;
1644 SCM value_of_wrong_type = SCM_INUM0;
1646 SCM result = SCM_UNDEFINED;
1652 if (scm_is_eq (m, SCM_INUM0))
1654 report_overflow = 1;
1658 if (!coerce_to_big (n, n_tmp))
1660 value_of_wrong_type = n;
1661 position_of_wrong_type = 1;
1665 if (!coerce_to_big (k, k_tmp))
1667 value_of_wrong_type = k;
1668 position_of_wrong_type = 2;
1672 if (!coerce_to_big (m, m_tmp))
1674 value_of_wrong_type = m;
1675 position_of_wrong_type = 3;
1679 /* if the exponent K is negative, and we simply call mpz_powm, we
1680 will get a divide-by-zero exception when an inverse 1/n mod m
1681 doesn't exist (or is not unique). Since exceptions are hard to
1682 handle, we'll attempt the inversion "by hand" -- that way, we get
1683 a simple failure code, which is easy to handle. */
1685 if (-1 == mpz_sgn (k_tmp))
1687 if (!mpz_invert (n_tmp, n_tmp, m_tmp))
1689 report_overflow = 1;
1692 mpz_neg (k_tmp, k_tmp);
1695 result = scm_i_mkbig ();
1696 mpz_powm (SCM_I_BIG_MPZ (result),
1701 if (mpz_sgn (m_tmp) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result)) != 0)
1702 mpz_add (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result), m_tmp);
1709 if (report_overflow)
1710 scm_num_overflow (FUNC_NAME);
1712 if (position_of_wrong_type)
1713 SCM_WRONG_TYPE_ARG (position_of_wrong_type,
1714 value_of_wrong_type);
1716 return scm_i_normbig (result);
1720 SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
1722 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1723 "exact integer, @var{n} can be any number.\n"
1725 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1726 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1727 "includes @math{0^0} is 1.\n"
1730 "(integer-expt 2 5) @result{} 32\n"
1731 "(integer-expt -3 3) @result{} -27\n"
1732 "(integer-expt 5 -3) @result{} 1/125\n"
1733 "(integer-expt 0 0) @result{} 1\n"
1735 #define FUNC_NAME s_scm_integer_expt
1738 SCM z_i2 = SCM_BOOL_F;
1740 SCM acc = SCM_I_MAKINUM (1L);
1742 /* 0^0 == 1 according to R5RS */
1743 if (scm_is_eq (n, SCM_INUM0) || scm_is_eq (n, acc))
1744 return scm_is_false (scm_zero_p(k)) ? n : acc;
1745 else if (scm_is_eq (n, SCM_I_MAKINUM (-1L)))
1746 return scm_is_false (scm_even_p (k)) ? n : acc;
1748 if (SCM_I_INUMP (k))
1749 i2 = SCM_I_INUM (k);
1750 else if (SCM_BIGP (k))
1752 z_i2 = scm_i_clonebig (k, 1);
1753 scm_remember_upto_here_1 (k);
1757 SCM_WRONG_TYPE_ARG (2, k);
1761 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2)) == -1)
1763 mpz_neg (SCM_I_BIG_MPZ (z_i2), SCM_I_BIG_MPZ (z_i2));
1764 n = scm_divide (n, SCM_UNDEFINED);
1768 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2)) == 0)
1772 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2), 1) == 0)
1774 return scm_product (acc, n);
1776 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2), 0))
1777 acc = scm_product (acc, n);
1778 n = scm_product (n, n);
1779 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2), SCM_I_BIG_MPZ (z_i2), 1);
1787 n = scm_divide (n, SCM_UNDEFINED);
1794 return scm_product (acc, n);
1796 acc = scm_product (acc, n);
1797 n = scm_product (n, n);
1804 SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
1806 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1807 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1809 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1810 "@var{cnt} is negative it's a division, rounded towards negative\n"
1811 "infinity. (Note that this is not the same rounding as\n"
1812 "@code{quotient} does.)\n"
1814 "With @var{n} viewed as an infinite precision twos complement,\n"
1815 "@code{ash} means a left shift introducing zero bits, or a right\n"
1816 "shift dropping bits.\n"
1819 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1820 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1822 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1823 "(ash -23 -2) @result{} -6\n"
1825 #define FUNC_NAME s_scm_ash
1828 bits_to_shift = scm_to_long (cnt);
1830 if (SCM_I_INUMP (n))
1832 long nn = SCM_I_INUM (n);
1834 if (bits_to_shift > 0)
1836 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1837 overflow a non-zero fixnum. For smaller shifts we check the
1838 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1839 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1840 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1846 if (bits_to_shift < SCM_I_FIXNUM_BIT-1
1848 (SCM_SRS (nn, (SCM_I_FIXNUM_BIT-1 - bits_to_shift)) + 1)
1851 return SCM_I_MAKINUM (nn << bits_to_shift);
1855 SCM result = scm_i_long2big (nn);
1856 mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
1863 bits_to_shift = -bits_to_shift;
1864 if (bits_to_shift >= SCM_LONG_BIT)
1865 return (nn >= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1867 return SCM_I_MAKINUM (SCM_SRS (nn, bits_to_shift));
1871 else if (SCM_BIGP (n))
1875 if (bits_to_shift == 0)
1878 result = scm_i_mkbig ();
1879 if (bits_to_shift >= 0)
1881 mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
1887 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1888 we have to allocate a bignum even if the result is going to be a
1890 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
1892 return scm_i_normbig (result);
1898 SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
1904 SCM_DEFINE (scm_bit_extract, "bit-extract", 3, 0, 0,
1905 (SCM n, SCM start, SCM end),
1906 "Return the integer composed of the @var{start} (inclusive)\n"
1907 "through @var{end} (exclusive) bits of @var{n}. The\n"
1908 "@var{start}th bit becomes the 0-th bit in the result.\n"
1911 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1912 " @result{} \"1010\"\n"
1913 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1914 " @result{} \"10110\"\n"
1916 #define FUNC_NAME s_scm_bit_extract
1918 unsigned long int istart, iend, bits;
1919 istart = scm_to_ulong (start);
1920 iend = scm_to_ulong (end);
1921 SCM_ASSERT_RANGE (3, end, (iend >= istart));
1923 /* how many bits to keep */
1924 bits = iend - istart;
1926 if (SCM_I_INUMP (n))
1928 long int in = SCM_I_INUM (n);
1930 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1931 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1932 in = SCM_SRS (in, min (istart, SCM_I_FIXNUM_BIT-1));
1934 if (in < 0 && bits >= SCM_I_FIXNUM_BIT)
1936 /* Since we emulate two's complement encoded numbers, this
1937 * special case requires us to produce a result that has
1938 * more bits than can be stored in a fixnum.
1940 SCM result = scm_i_long2big (in);
1941 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
1946 /* mask down to requisite bits */
1947 bits = min (bits, SCM_I_FIXNUM_BIT);
1948 return SCM_I_MAKINUM (in & ((1L << bits) - 1));
1950 else if (SCM_BIGP (n))
1955 result = SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n), istart));
1959 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1960 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1961 such bits into a ulong. */
1962 result = scm_i_mkbig ();
1963 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result), SCM_I_BIG_MPZ(n), istart);
1964 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result), SCM_I_BIG_MPZ(result), bits);
1965 result = scm_i_normbig (result);
1967 scm_remember_upto_here_1 (n);
1971 SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
1976 static const char scm_logtab[] = {
1977 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1980 SCM_DEFINE (scm_logcount, "logcount", 1, 0, 0,
1982 "Return the number of bits in integer @var{n}. If integer is\n"
1983 "positive, the 1-bits in its binary representation are counted.\n"
1984 "If negative, the 0-bits in its two's-complement binary\n"
1985 "representation are counted. If 0, 0 is returned.\n"
1988 "(logcount #b10101010)\n"
1995 #define FUNC_NAME s_scm_logcount
1997 if (SCM_I_INUMP (n))
1999 unsigned long int c = 0;
2000 long int nn = SCM_I_INUM (n);
2005 c += scm_logtab[15 & nn];
2008 return SCM_I_MAKINUM (c);
2010 else if (SCM_BIGP (n))
2012 unsigned long count;
2013 if (mpz_sgn (SCM_I_BIG_MPZ (n)) >= 0)
2014 count = mpz_popcount (SCM_I_BIG_MPZ (n));
2016 count = mpz_hamdist (SCM_I_BIG_MPZ (n), z_negative_one);
2017 scm_remember_upto_here_1 (n);
2018 return SCM_I_MAKINUM (count);
2021 SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
2026 static const char scm_ilentab[] = {
2027 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2031 SCM_DEFINE (scm_integer_length, "integer-length", 1, 0, 0,
2033 "Return the number of bits necessary to represent @var{n}.\n"
2036 "(integer-length #b10101010)\n"
2038 "(integer-length 0)\n"
2040 "(integer-length #b1111)\n"
2043 #define FUNC_NAME s_scm_integer_length
2045 if (SCM_I_INUMP (n))
2047 unsigned long int c = 0;
2049 long int nn = SCM_I_INUM (n);
2055 l = scm_ilentab [15 & nn];
2058 return SCM_I_MAKINUM (c - 4 + l);
2060 else if (SCM_BIGP (n))
2062 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2063 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2064 1 too big, so check for that and adjust. */
2065 size_t size = mpz_sizeinbase (SCM_I_BIG_MPZ (n), 2);
2066 if (mpz_sgn (SCM_I_BIG_MPZ (n)) < 0
2067 && mpz_scan0 (SCM_I_BIG_MPZ (n), /* no 0 bits above the lowest 1 */
2068 mpz_scan1 (SCM_I_BIG_MPZ (n), 0)) == ULONG_MAX)
2070 scm_remember_upto_here_1 (n);
2071 return SCM_I_MAKINUM (size);
2074 SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
2078 /*** NUMBERS -> STRINGS ***/
2079 #define SCM_MAX_DBL_PREC 60
2080 #define SCM_MAX_DBL_RADIX 36
2082 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2083 static int scm_dblprec[SCM_MAX_DBL_RADIX - 1];
2084 static double fx_per_radix[SCM_MAX_DBL_RADIX - 1][SCM_MAX_DBL_PREC];
2087 void init_dblprec(int *prec, int radix) {
2088 /* determine floating point precision by adding successively
2089 smaller increments to 1.0 until it is considered == 1.0 */
2090 double f = ((double)1.0)/radix;
2091 double fsum = 1.0 + f;
2096 if (++(*prec) > SCM_MAX_DBL_PREC)
2108 void init_fx_radix(double *fx_list, int radix)
2110 /* initialize a per-radix list of tolerances. When added
2111 to a number < 1.0, we can determine if we should raund
2112 up and quit converting a number to a string. */
2116 for( i = 2 ; i < SCM_MAX_DBL_PREC; ++i )
2117 fx_list[i] = (fx_list[i-1] / radix);
2120 /* use this array as a way to generate a single digit */
2121 static const char*number_chars="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2124 idbl2str (double f, char *a, int radix)
2126 int efmt, dpt, d, i, wp;
2128 #ifdef DBL_MIN_10_EXP
2131 #endif /* DBL_MIN_10_EXP */
2136 radix > SCM_MAX_DBL_RADIX)
2138 /* revert to existing behavior */
2142 wp = scm_dblprec[radix-2];
2143 fx = fx_per_radix[radix-2];
2147 #ifdef HAVE_COPYSIGN
2148 double sgn = copysign (1.0, f);
2153 goto zero; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2159 strcpy (a, "-inf.0");
2161 strcpy (a, "+inf.0");
2164 else if (xisnan (f))
2166 strcpy (a, "+nan.0");
2176 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2177 make-uniform-vector, from causing infinite loops. */
2178 /* just do the checking...if it passes, we do the conversion for our
2179 radix again below */
2186 if (exp_cpy-- < DBL_MIN_10_EXP)
2194 while (f_cpy > 10.0)
2197 if (exp_cpy++ > DBL_MAX_10_EXP)
2218 if (f + fx[wp] >= radix)
2225 /* adding 9999 makes this equivalent to abs(x) % 3 */
2226 dpt = (exp + 9999) % 3;
2230 efmt = (exp < -3) || (exp > wp + 2);
2252 a[ch++] = number_chars[d];
2255 if (f + fx[wp] >= 1.0)
2257 a[ch - 1] = number_chars[d+1];
2269 if ((dpt > 4) && (exp > 6))
2271 d = (a[0] == '-' ? 2 : 1);
2272 for (i = ch++; i > d; i--)
2285 if (a[ch - 1] == '.')
2286 a[ch++] = '0'; /* trailing zero */
2295 for (i = radix; i <= exp; i *= radix);
2296 for (i /= radix; i; i /= radix)
2298 a[ch++] = number_chars[exp / i];
2307 icmplx2str (double real, double imag, char *str, int radix)
2311 i = idbl2str (real, str, radix);
2314 /* Don't output a '+' for negative numbers or for Inf and
2315 NaN. They will provide their own sign. */
2316 if (0 <= imag && !xisinf (imag) && !xisnan (imag))
2318 i += idbl2str (imag, &str[i], radix);
2325 iflo2str (SCM flt, char *str, int radix)
2328 if (SCM_REALP (flt))
2329 i = idbl2str (SCM_REAL_VALUE (flt), str, radix);
2331 i = icmplx2str (SCM_COMPLEX_REAL (flt), SCM_COMPLEX_IMAG (flt),
2336 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2337 characters in the result.
2339 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2341 scm_iint2str (scm_t_intmax num, int rad, char *p)
2346 return scm_iuint2str (-num, rad, p) + 1;
2349 return scm_iuint2str (num, rad, p);
2352 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2353 characters in the result.
2355 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2357 scm_iuint2str (scm_t_uintmax num, int rad, char *p)
2361 scm_t_uintmax n = num;
2363 for (n /= rad; n > 0; n /= rad)
2373 p[i] = d + ((d < 10) ? '0' : 'a' - 10);
2378 SCM_DEFINE (scm_number_to_string, "number->string", 1, 1, 0,
2380 "Return a string holding the external representation of the\n"
2381 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2382 "inexact, a radix of 10 will be used.")
2383 #define FUNC_NAME s_scm_number_to_string
2387 if (SCM_UNBNDP (radix))
2390 base = scm_to_signed_integer (radix, 2, 36);
2392 if (SCM_I_INUMP (n))
2394 char num_buf [SCM_INTBUFLEN];
2395 size_t length = scm_iint2str (SCM_I_INUM (n), base, num_buf);
2396 return scm_from_locale_stringn (num_buf, length);
2398 else if (SCM_BIGP (n))
2400 char *str = mpz_get_str (NULL, base, SCM_I_BIG_MPZ (n));
2401 scm_remember_upto_here_1 (n);
2402 return scm_take_locale_string (str);
2404 else if (SCM_FRACTIONP (n))
2406 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n), radix),
2407 scm_from_locale_string ("/"),
2408 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n), radix)));
2410 else if (SCM_INEXACTP (n))
2412 char num_buf [FLOBUFLEN];
2413 return scm_from_locale_stringn (num_buf, iflo2str (n, num_buf, base));
2416 SCM_WRONG_TYPE_ARG (1, n);
2421 /* These print routines used to be stubbed here so that scm_repl.c
2422 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2425 scm_print_real (SCM sexp, SCM port, scm_print_state *pstate SCM_UNUSED)
2427 char num_buf[FLOBUFLEN];
2428 scm_lfwrite (num_buf, iflo2str (sexp, num_buf, 10), port);
2433 scm_i_print_double (double val, SCM port)
2435 char num_buf[FLOBUFLEN];
2436 scm_lfwrite (num_buf, idbl2str (val, num_buf, 10), port);
2440 scm_print_complex (SCM sexp, SCM port, scm_print_state *pstate SCM_UNUSED)
2443 char num_buf[FLOBUFLEN];
2444 scm_lfwrite (num_buf, iflo2str (sexp, num_buf, 10), port);
2449 scm_i_print_complex (double real, double imag, SCM port)
2451 char num_buf[FLOBUFLEN];
2452 scm_lfwrite (num_buf, icmplx2str (real, imag, num_buf, 10), port);
2456 scm_i_print_fraction (SCM sexp, SCM port, scm_print_state *pstate SCM_UNUSED)
2459 str = scm_number_to_string (sexp, SCM_UNDEFINED);
2460 scm_lfwrite (scm_i_string_chars (str), scm_i_string_length (str), port);
2461 scm_remember_upto_here_1 (str);
2466 scm_bigprint (SCM exp, SCM port, scm_print_state *pstate SCM_UNUSED)
2468 char *str = mpz_get_str (NULL, 10, SCM_I_BIG_MPZ (exp));
2469 scm_remember_upto_here_1 (exp);
2470 scm_lfwrite (str, (size_t) strlen (str), port);
2474 /*** END nums->strs ***/
2477 /*** STRINGS -> NUMBERS ***/
2479 /* The following functions implement the conversion from strings to numbers.
2480 * The implementation somehow follows the grammar for numbers as it is given
2481 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2482 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2483 * points should be noted about the implementation:
2484 * * Each function keeps a local index variable 'idx' that points at the
2485 * current position within the parsed string. The global index is only
2486 * updated if the function could parse the corresponding syntactic unit
2488 * * Similarly, the functions keep track of indicators of inexactness ('#',
2489 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2490 * global exactness information is only updated after each part has been
2491 * successfully parsed.
2492 * * Sequences of digits are parsed into temporary variables holding fixnums.
2493 * Only if these fixnums would overflow, the result variables are updated
2494 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2495 * the temporary variables holding the fixnums are cleared, and the process
2496 * starts over again. If for example fixnums were able to store five decimal
2497 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2498 * and the result was computed as 12345 * 100000 + 67890. In other words,
2499 * only every five digits two bignum operations were performed.
2502 enum t_exactness {NO_EXACTNESS, INEXACT, EXACT};
2504 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2506 /* In non ASCII-style encodings the following macro might not work. */
2507 #define XDIGIT2UINT(d) \
2508 (isdigit ((int) (unsigned char) d) \
2510 : tolower ((int) (unsigned char) d) - 'a' + 10)
2513 mem2uinteger (const char* mem, size_t len, unsigned int *p_idx,
2514 unsigned int radix, enum t_exactness *p_exactness)
2516 unsigned int idx = *p_idx;
2517 unsigned int hash_seen = 0;
2518 scm_t_bits shift = 1;
2520 unsigned int digit_value;
2528 if (!isxdigit ((int) (unsigned char) c))
2530 digit_value = XDIGIT2UINT (c);
2531 if (digit_value >= radix)
2535 result = SCM_I_MAKINUM (digit_value);
2539 if (isxdigit ((int) (unsigned char) c))
2543 digit_value = XDIGIT2UINT (c);
2544 if (digit_value >= radix)
2556 if (SCM_MOST_POSITIVE_FIXNUM / radix < shift)
2558 result = scm_product (result, SCM_I_MAKINUM (shift));
2560 result = scm_sum (result, SCM_I_MAKINUM (add));
2567 shift = shift * radix;
2568 add = add * radix + digit_value;
2573 result = scm_product (result, SCM_I_MAKINUM (shift));
2575 result = scm_sum (result, SCM_I_MAKINUM (add));
2579 *p_exactness = INEXACT;
2585 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2586 * covers the parts of the rules that start at a potential point. The value
2587 * of the digits up to the point have been parsed by the caller and are given
2588 * in variable result. The content of *p_exactness indicates, whether a hash
2589 * has already been seen in the digits before the point.
2592 /* In non ASCII-style encodings the following macro might not work. */
2593 #define DIGIT2UINT(d) ((d) - '0')
2596 mem2decimal_from_point (SCM result, const char* mem, size_t len,
2597 unsigned int *p_idx, enum t_exactness *p_exactness)
2599 unsigned int idx = *p_idx;
2600 enum t_exactness x = *p_exactness;
2605 if (mem[idx] == '.')
2607 scm_t_bits shift = 1;
2609 unsigned int digit_value;
2610 SCM big_shift = SCM_I_MAKINUM (1);
2616 if (isdigit ((int) (unsigned char) c))
2621 digit_value = DIGIT2UINT (c);
2632 if (SCM_MOST_POSITIVE_FIXNUM / 10 < shift)
2634 big_shift = scm_product (big_shift, SCM_I_MAKINUM (shift));
2635 result = scm_product (result, SCM_I_MAKINUM (shift));
2637 result = scm_sum (result, SCM_I_MAKINUM (add));
2645 add = add * 10 + digit_value;
2651 big_shift = scm_product (big_shift, SCM_I_MAKINUM (shift));
2652 result = scm_product (result, SCM_I_MAKINUM (shift));
2653 result = scm_sum (result, SCM_I_MAKINUM (add));
2656 result = scm_divide (result, big_shift);
2658 /* We've seen a decimal point, thus the value is implicitly inexact. */
2670 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2706 if (!isdigit ((int) (unsigned char) c))
2710 exponent = DIGIT2UINT (c);
2714 if (isdigit ((int) (unsigned char) c))
2717 if (exponent <= SCM_MAXEXP)
2718 exponent = exponent * 10 + DIGIT2UINT (c);
2724 if (exponent > SCM_MAXEXP)
2726 size_t exp_len = idx - start;
2727 SCM exp_string = scm_from_locale_stringn (&mem[start], exp_len);
2728 SCM exp_num = scm_string_to_number (exp_string, SCM_UNDEFINED);
2729 scm_out_of_range ("string->number", exp_num);
2732 e = scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent));
2734 result = scm_product (result, e);
2736 result = scm_divide2real (result, e);
2738 /* We've seen an exponent, thus the value is implicitly inexact. */
2756 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2759 mem2ureal (const char* mem, size_t len, unsigned int *p_idx,
2760 unsigned int radix, enum t_exactness *p_exactness)
2762 unsigned int idx = *p_idx;
2765 /* Start off believing that the number will be exact. This changes
2766 to INEXACT if we see a decimal point or a hash. */
2767 enum t_exactness x = EXACT;
2772 if (idx+5 <= len && !strncmp (mem+idx, "inf.0", 5))
2778 if (idx+4 < len && !strncmp (mem+idx, "nan.", 4))
2780 /* Cobble up the fractional part. We might want to set the
2781 NaN's mantissa from it. */
2783 mem2uinteger (mem, len, &idx, 10, &x);
2788 if (mem[idx] == '.')
2792 else if (idx + 1 == len)
2794 else if (!isdigit ((int) (unsigned char) mem[idx + 1]))
2797 result = mem2decimal_from_point (SCM_I_MAKINUM (0), mem, len,
2804 uinteger = mem2uinteger (mem, len, &idx, radix, &x);
2805 if (scm_is_false (uinteger))
2810 else if (mem[idx] == '/')
2818 divisor = mem2uinteger (mem, len, &idx, radix, &x);
2819 if (scm_is_false (divisor))
2822 /* both are int/big here, I assume */
2823 result = scm_i_make_ratio (uinteger, divisor);
2825 else if (radix == 10)
2827 result = mem2decimal_from_point (uinteger, mem, len, &idx, &x);
2828 if (scm_is_false (result))
2837 /* Update *p_exactness if the number just read was inexact. This is
2838 important for complex numbers, so that a complex number is
2839 treated as inexact overall if either its real or imaginary part
2845 /* When returning an inexact zero, make sure it is represented as a
2846 floating point value so that we can change its sign.
2848 if (scm_is_eq (result, SCM_I_MAKINUM(0)) && *p_exactness == INEXACT)
2849 result = scm_from_double (0.0);
2855 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2858 mem2complex (const char* mem, size_t len, unsigned int idx,
2859 unsigned int radix, enum t_exactness *p_exactness)
2883 ureal = mem2ureal (mem, len, &idx, radix, p_exactness);
2884 if (scm_is_false (ureal))
2886 /* input must be either +i or -i */
2891 if (mem[idx] == 'i' || mem[idx] == 'I')
2897 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign));
2904 if (sign == -1 && scm_is_false (scm_nan_p (ureal)))
2905 ureal = scm_difference (ureal, SCM_UNDEFINED);
2914 /* either +<ureal>i or -<ureal>i */
2921 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal);
2924 /* polar input: <real>@<real>. */
2953 angle = mem2ureal (mem, len, &idx, radix, p_exactness);
2954 if (scm_is_false (angle))
2959 if (sign == -1 && scm_is_false (scm_nan_p (ureal)))
2960 angle = scm_difference (angle, SCM_UNDEFINED);
2962 result = scm_make_polar (ureal, angle);
2967 /* expecting input matching <real>[+-]<ureal>?i */
2974 int sign = (c == '+') ? 1 : -1;
2975 SCM imag = mem2ureal (mem, len, &idx, radix, p_exactness);
2977 if (scm_is_false (imag))
2978 imag = SCM_I_MAKINUM (sign);
2979 else if (sign == -1 && scm_is_false (scm_nan_p (ureal)))
2980 imag = scm_difference (imag, SCM_UNDEFINED);
2984 if (mem[idx] != 'i' && mem[idx] != 'I')
2991 return scm_make_rectangular (ureal, imag);
3000 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
3002 enum t_radix {NO_RADIX=0, DUAL=2, OCT=8, DEC=10, HEX=16};
3005 scm_c_locale_stringn_to_number (const char* mem, size_t len,
3006 unsigned int default_radix)
3008 unsigned int idx = 0;
3009 unsigned int radix = NO_RADIX;
3010 enum t_exactness forced_x = NO_EXACTNESS;
3011 enum t_exactness implicit_x = EXACT;
3014 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3015 while (idx + 2 < len && mem[idx] == '#')
3017 switch (mem[idx + 1])
3020 if (radix != NO_RADIX)
3025 if (radix != NO_RADIX)
3030 if (forced_x != NO_EXACTNESS)
3035 if (forced_x != NO_EXACTNESS)
3040 if (radix != NO_RADIX)
3045 if (radix != NO_RADIX)
3055 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3056 if (radix == NO_RADIX)
3057 result = mem2complex (mem, len, idx, default_radix, &implicit_x);
3059 result = mem2complex (mem, len, idx, (unsigned int) radix, &implicit_x);
3061 if (scm_is_false (result))
3067 if (SCM_INEXACTP (result))
3068 return scm_inexact_to_exact (result);
3072 if (SCM_INEXACTP (result))
3075 return scm_exact_to_inexact (result);
3078 if (implicit_x == INEXACT)
3080 if (SCM_INEXACTP (result))
3083 return scm_exact_to_inexact (result);
3091 SCM_DEFINE (scm_string_to_number, "string->number", 1, 1, 0,
3092 (SCM string, SCM radix),
3093 "Return a number of the maximally precise representation\n"
3094 "expressed by the given @var{string}. @var{radix} must be an\n"
3095 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3096 "is a default radix that may be overridden by an explicit radix\n"
3097 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3098 "supplied, then the default radix is 10. If string is not a\n"
3099 "syntactically valid notation for a number, then\n"
3100 "@code{string->number} returns @code{#f}.")
3101 #define FUNC_NAME s_scm_string_to_number
3105 SCM_VALIDATE_STRING (1, string);
3107 if (SCM_UNBNDP (radix))
3110 base = scm_to_unsigned_integer (radix, 2, INT_MAX);
3112 answer = scm_c_locale_stringn_to_number (scm_i_string_chars (string),
3113 scm_i_string_length (string),
3115 scm_remember_upto_here_1 (string);
3121 /*** END strs->nums ***/
3125 scm_bigequal (SCM x, SCM y)
3127 int result = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
3128 scm_remember_upto_here_2 (x, y);
3129 return scm_from_bool (0 == result);
3133 scm_real_equalp (SCM x, SCM y)
3135 return scm_from_bool (SCM_REAL_VALUE (x) == SCM_REAL_VALUE (y));
3139 scm_complex_equalp (SCM x, SCM y)
3141 return scm_from_bool (SCM_COMPLEX_REAL (x) == SCM_COMPLEX_REAL (y)
3142 && SCM_COMPLEX_IMAG (x) == SCM_COMPLEX_IMAG (y));
3146 scm_i_fraction_equalp (SCM x, SCM y)
3148 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x),
3149 SCM_FRACTION_NUMERATOR (y)))
3150 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x),
3151 SCM_FRACTION_DENOMINATOR (y))))
3158 SCM_DEFINE (scm_number_p, "number?", 1, 0, 0,
3160 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3162 #define FUNC_NAME s_scm_number_p
3164 return scm_from_bool (SCM_NUMBERP (x));
3168 SCM_DEFINE (scm_complex_p, "complex?", 1, 0, 0,
3170 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3171 "otherwise. Note that the sets of real, rational and integer\n"
3172 "values form subsets of the set of complex numbers, i. e. the\n"
3173 "predicate will also be fulfilled if @var{x} is a real,\n"
3174 "rational or integer number.")
3175 #define FUNC_NAME s_scm_complex_p
3177 /* all numbers are complex. */
3178 return scm_number_p (x);
3182 SCM_DEFINE (scm_real_p, "real?", 1, 0, 0,
3184 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3185 "otherwise. Note that the set of integer values forms a subset of\n"
3186 "the set of real numbers, i. e. the predicate will also be\n"
3187 "fulfilled if @var{x} is an integer number.")
3188 #define FUNC_NAME s_scm_real_p
3190 /* we can't represent irrational numbers. */
3191 return scm_rational_p (x);
3195 SCM_DEFINE (scm_rational_p, "rational?", 1, 0, 0,
3197 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3198 "otherwise. Note that the set of integer values forms a subset of\n"
3199 "the set of rational numbers, i. e. the predicate will also be\n"
3200 "fulfilled if @var{x} is an integer number.")
3201 #define FUNC_NAME s_scm_rational_p
3203 if (SCM_I_INUMP (x))
3205 else if (SCM_IMP (x))
3207 else if (SCM_BIGP (x))
3209 else if (SCM_FRACTIONP (x))
3211 else if (SCM_REALP (x))
3212 /* due to their limited precision, all floating point numbers are
3213 rational as well. */
3220 SCM_DEFINE (scm_integer_p, "integer?", 1, 0, 0,
3222 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3224 #define FUNC_NAME s_scm_integer_p
3227 if (SCM_I_INUMP (x))
3233 if (!SCM_INEXACTP (x))
3235 if (SCM_COMPLEXP (x))
3237 r = SCM_REAL_VALUE (x);
3238 /* +/-inf passes r==floor(r), making those #t */
3246 SCM_DEFINE (scm_inexact_p, "inexact?", 1, 0, 0,
3248 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3250 #define FUNC_NAME s_scm_inexact_p
3252 if (SCM_INEXACTP (x))
3254 if (SCM_NUMBERP (x))
3256 SCM_WRONG_TYPE_ARG (1, x);
3261 SCM_GPROC1 (s_eq_p, "=", scm_tc7_rpsubr, scm_num_eq_p, g_eq_p);
3262 /* "Return @code{#t} if all parameters are numerically equal." */
3264 scm_num_eq_p (SCM x, SCM y)
3267 if (SCM_I_INUMP (x))
3269 long xx = SCM_I_INUM (x);
3270 if (SCM_I_INUMP (y))
3272 long yy = SCM_I_INUM (y);
3273 return scm_from_bool (xx == yy);
3275 else if (SCM_BIGP (y))
3277 else if (SCM_REALP (y))
3279 /* On a 32-bit system an inum fits a double, we can cast the inum
3280 to a double and compare.
3282 But on a 64-bit system an inum is bigger than a double and
3283 casting it to a double (call that dxx) will round. dxx is at
3284 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3285 an integer and fits a long. So we cast yy to a long and
3286 compare with plain xx.
3288 An alternative (for any size system actually) would be to check
3289 yy is an integer (with floor) and is in range of an inum
3290 (compare against appropriate powers of 2) then test
3291 xx==(long)yy. It's just a matter of which casts/comparisons
3292 might be fastest or easiest for the cpu. */
3294 double yy = SCM_REAL_VALUE (y);
3295 return scm_from_bool ((double) xx == yy
3296 && (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1
3297 || xx == (long) yy));
3299 else if (SCM_COMPLEXP (y))
3300 return scm_from_bool (((double) xx == SCM_COMPLEX_REAL (y))
3301 && (0.0 == SCM_COMPLEX_IMAG (y)));
3302 else if (SCM_FRACTIONP (y))
3305 SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
3307 else if (SCM_BIGP (x))
3309 if (SCM_I_INUMP (y))
3311 else if (SCM_BIGP (y))
3313 int cmp = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
3314 scm_remember_upto_here_2 (x, y);
3315 return scm_from_bool (0 == cmp);
3317 else if (SCM_REALP (y))
3320 if (xisnan (SCM_REAL_VALUE (y)))
3322 cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (x), SCM_REAL_VALUE (y));
3323 scm_remember_upto_here_1 (x);
3324 return scm_from_bool (0 == cmp);
3326 else if (SCM_COMPLEXP (y))
3329 if (0.0 != SCM_COMPLEX_IMAG (y))
3331 if (xisnan (SCM_COMPLEX_REAL (y)))
3333 cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (x), SCM_COMPLEX_REAL (y));
3334 scm_remember_upto_here_1 (x);
3335 return scm_from_bool (0 == cmp);
3337 else if (SCM_FRACTIONP (y))
3340 SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
3342 else if (SCM_REALP (x))
3344 double xx = SCM_REAL_VALUE (x);
3345 if (SCM_I_INUMP (y))
3347 /* see comments with inum/real above */
3348 long yy = SCM_I_INUM (y);
3349 return scm_from_bool (xx == (double) yy
3350 && (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1
3351 || (long) xx == yy));
3353 else if (SCM_BIGP (y))
3356 if (xisnan (SCM_REAL_VALUE (x)))
3358 cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), SCM_REAL_VALUE (x));
3359 scm_remember_upto_here_1 (y);
3360 return scm_from_bool (0 == cmp);
3362 else if (SCM_REALP (y))
3363 return scm_from_bool (SCM_REAL_VALUE (x) == SCM_REAL_VALUE (y));
3364 else if (SCM_COMPLEXP (y))
3365 return scm_from_bool ((SCM_REAL_VALUE (x) == SCM_COMPLEX_REAL (y))
3366 && (0.0 == SCM_COMPLEX_IMAG (y)));
3367 else if (SCM_FRACTIONP (y))
3369 double xx = SCM_REAL_VALUE (x);
3373 return scm_from_bool (xx < 0.0);
3374 x = scm_inexact_to_exact (x); /* with x as frac or int */
3378 SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
3380 else if (SCM_COMPLEXP (x))
3382 if (SCM_I_INUMP (y))
3383 return scm_from_bool ((SCM_COMPLEX_REAL (x) == (double) SCM_I_INUM (y))
3384 && (SCM_COMPLEX_IMAG (x) == 0.0));
3385 else if (SCM_BIGP (y))
3388 if (0.0 != SCM_COMPLEX_IMAG (x))
3390 if (xisnan (SCM_COMPLEX_REAL (x)))
3392 cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), SCM_COMPLEX_REAL (x));
3393 scm_remember_upto_here_1 (y);
3394 return scm_from_bool (0 == cmp);
3396 else if (SCM_REALP (y))
3397 return scm_from_bool ((SCM_COMPLEX_REAL (x) == SCM_REAL_VALUE (y))
3398 && (SCM_COMPLEX_IMAG (x) == 0.0));
3399 else if (SCM_COMPLEXP (y))
3400 return scm_from_bool ((SCM_COMPLEX_REAL (x) == SCM_COMPLEX_REAL (y))
3401 && (SCM_COMPLEX_IMAG (x) == SCM_COMPLEX_IMAG (y)));
3402 else if (SCM_FRACTIONP (y))
3405 if (SCM_COMPLEX_IMAG (x) != 0.0)
3407 xx = SCM_COMPLEX_REAL (x);
3411 return scm_from_bool (xx < 0.0);
3412 x = scm_inexact_to_exact (x); /* with x as frac or int */
3416 SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
3418 else if (SCM_FRACTIONP (x))
3420 if (SCM_I_INUMP (y))
3422 else if (SCM_BIGP (y))
3424 else if (SCM_REALP (y))
3426 double yy = SCM_REAL_VALUE (y);
3430 return scm_from_bool (0.0 < yy);
3431 y = scm_inexact_to_exact (y); /* with y as frac or int */
3434 else if (SCM_COMPLEXP (y))
3437 if (SCM_COMPLEX_IMAG (y) != 0.0)
3439 yy = SCM_COMPLEX_REAL (y);
3443 return scm_from_bool (0.0 < yy);
3444 y = scm_inexact_to_exact (y); /* with y as frac or int */
3447 else if (SCM_FRACTIONP (y))
3448 return scm_i_fraction_equalp (x, y);
3450 SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARGn, s_eq_p);
3453 SCM_WTA_DISPATCH_2 (g_eq_p, x, y, SCM_ARG1, s_eq_p);
3457 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3458 done are good for inums, but for bignums an answer can almost always be
3459 had by just examining a few high bits of the operands, as done by GMP in
3460 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3461 of the float exponent to take into account. */
3463 SCM_GPROC1 (s_less_p, "<", scm_tc7_rpsubr, scm_less_p, g_less_p);
3464 /* "Return @code{#t} if the list of parameters is monotonically\n"
3468 scm_less_p (SCM x, SCM y)
3471 if (SCM_I_INUMP (x))
3473 long xx = SCM_I_INUM (x);
3474 if (SCM_I_INUMP (y))
3476 long yy = SCM_I_INUM (y);
3477 return scm_from_bool (xx < yy);
3479 else if (SCM_BIGP (y))
3481 int sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
3482 scm_remember_upto_here_1 (y);
3483 return scm_from_bool (sgn > 0);
3485 else if (SCM_REALP (y))
3486 return scm_from_bool ((double) xx < SCM_REAL_VALUE (y));
3487 else if (SCM_FRACTIONP (y))
3489 /* "x < a/b" becomes "x*b < a" */
3491 x = scm_product (x, SCM_FRACTION_DENOMINATOR (y));
3492 y = SCM_FRACTION_NUMERATOR (y);
3496 SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
3498 else if (SCM_BIGP (x))
3500 if (SCM_I_INUMP (y))
3502 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
3503 scm_remember_upto_here_1 (x);
3504 return scm_from_bool (sgn < 0);
3506 else if (SCM_BIGP (y))
3508 int cmp = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
3509 scm_remember_upto_here_2 (x, y);
3510 return scm_from_bool (cmp < 0);
3512 else if (SCM_REALP (y))
3515 if (xisnan (SCM_REAL_VALUE (y)))
3517 cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (x), SCM_REAL_VALUE (y));
3518 scm_remember_upto_here_1 (x);
3519 return scm_from_bool (cmp < 0);
3521 else if (SCM_FRACTIONP (y))
3524 SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
3526 else if (SCM_REALP (x))
3528 if (SCM_I_INUMP (y))
3529 return scm_from_bool (SCM_REAL_VALUE (x) < (double) SCM_I_INUM (y));
3530 else if (SCM_BIGP (y))
3533 if (xisnan (SCM_REAL_VALUE (x)))
3535 cmp = xmpz_cmp_d (SCM_I_BIG_MPZ (y), SCM_REAL_VALUE (x));
3536 scm_remember_upto_here_1 (y);
3537 return scm_from_bool (cmp > 0);
3539 else if (SCM_REALP (y))
3540 return scm_from_bool (SCM_REAL_VALUE (x) < SCM_REAL_VALUE (y));
3541 else if (SCM_FRACTIONP (y))
3543 double xx = SCM_REAL_VALUE (x);
3547 return scm_from_bool (xx < 0.0);
3548 x = scm_inexact_to_exact (x); /* with x as frac or int */
3552 SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
3554 else if (SCM_FRACTIONP (x))
3556 if (SCM_I_INUMP (y) || SCM_BIGP (y))
3558 /* "a/b < y" becomes "a < y*b" */
3559 y = scm_product (y, SCM_FRACTION_DENOMINATOR (x));
3560 x = SCM_FRACTION_NUMERATOR (x);
3563 else if (SCM_REALP (y))
3565 double yy = SCM_REAL_VALUE (y);
3569 return scm_from_bool (0.0 < yy);
3570 y = scm_inexact_to_exact (y); /* with y as frac or int */
3573 else if (SCM_FRACTIONP (y))
3575 /* "a/b < c/d" becomes "a*d < c*b" */
3576 SCM new_x = scm_product (SCM_FRACTION_NUMERATOR (x),
3577 SCM_FRACTION_DENOMINATOR (y));
3578 SCM new_y = scm_product (SCM_FRACTION_NUMERATOR (y),
3579 SCM_FRACTION_DENOMINATOR (x));
3585 SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARGn, s_less_p);
3588 SCM_WTA_DISPATCH_2 (g_less_p, x, y, SCM_ARG1, s_less_p);
3592 SCM_GPROC1 (s_scm_gr_p, ">", scm_tc7_rpsubr, scm_gr_p, g_gr_p);
3593 /* "Return @code{#t} if the list of parameters is monotonically\n"
3596 #define FUNC_NAME s_scm_gr_p
3598 scm_gr_p (SCM x, SCM y)
3600 if (!SCM_NUMBERP (x))
3601 SCM_WTA_DISPATCH_2 (g_gr_p, x, y, SCM_ARG1, FUNC_NAME);
3602 else if (!SCM_NUMBERP (y))
3603 SCM_WTA_DISPATCH_2 (g_gr_p, x, y, SCM_ARG2, FUNC_NAME);
3605 return scm_less_p (y, x);
3610 SCM_GPROC1 (s_scm_leq_p, "<=", scm_tc7_rpsubr, scm_leq_p, g_leq_p);
3611 /* "Return @code{#t} if the list of parameters is monotonically\n"
3614 #define FUNC_NAME s_scm_leq_p
3616 scm_leq_p (SCM x, SCM y)
3618 if (!SCM_NUMBERP (x))
3619 SCM_WTA_DISPATCH_2 (g_leq_p, x, y, SCM_ARG1, FUNC_NAME);
3620 else if (!SCM_NUMBERP (y))
3621 SCM_WTA_DISPATCH_2 (g_leq_p, x, y, SCM_ARG2, FUNC_NAME);
3622 else if (scm_is_true (scm_nan_p (x)) || scm_is_true (scm_nan_p (y)))
3625 return scm_not (scm_less_p (y, x));
3630 SCM_GPROC1 (s_scm_geq_p, ">=", scm_tc7_rpsubr, scm_geq_p, g_geq_p);
3631 /* "Return @code{#t} if the list of parameters is monotonically\n"
3634 #define FUNC_NAME s_scm_geq_p
3636 scm_geq_p (SCM x, SCM y)
3638 if (!SCM_NUMBERP (x))
3639 SCM_WTA_DISPATCH_2 (g_geq_p, x, y, SCM_ARG1, FUNC_NAME);
3640 else if (!SCM_NUMBERP (y))
3641 SCM_WTA_DISPATCH_2 (g_geq_p, x, y, SCM_ARG2, FUNC_NAME);
3642 else if (scm_is_true (scm_nan_p (x)) || scm_is_true (scm_nan_p (y)))
3645 return scm_not (scm_less_p (x, y));
3650 SCM_GPROC (s_zero_p, "zero?", 1, 0, 0, scm_zero_p, g_zero_p);
3651 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3657 if (SCM_I_INUMP (z))
3658 return scm_from_bool (scm_is_eq (z, SCM_INUM0));
3659 else if (SCM_BIGP (z))
3661 else if (SCM_REALP (z))
3662 return scm_from_bool (SCM_REAL_VALUE (z) == 0.0);
3663 else if (SCM_COMPLEXP (z))
3664 return scm_from_bool (SCM_COMPLEX_REAL (z) == 0.0
3665 && SCM_COMPLEX_IMAG (z) == 0.0);
3666 else if (SCM_FRACTIONP (z))
3669 SCM_WTA_DISPATCH_1 (g_zero_p, z, SCM_ARG1, s_zero_p);
3673 SCM_GPROC (s_positive_p, "positive?", 1, 0, 0, scm_positive_p, g_positive_p);
3674 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3678 scm_positive_p (SCM x)
3680 if (SCM_I_INUMP (x))
3681 return scm_from_bool (SCM_I_INUM (x) > 0);
3682 else if (SCM_BIGP (x))
3684 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
3685 scm_remember_upto_here_1 (x);
3686 return scm_from_bool (sgn > 0);
3688 else if (SCM_REALP (x))
3689 return scm_from_bool(SCM_REAL_VALUE (x) > 0.0);
3690 else if (SCM_FRACTIONP (x))
3691 return scm_positive_p (SCM_FRACTION_NUMERATOR (x));
3693 SCM_WTA_DISPATCH_1 (g_positive_p, x, SCM_ARG1, s_positive_p);
3697 SCM_GPROC (s_negative_p, "negative?", 1, 0, 0, scm_negative_p, g_negative_p);
3698 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3702 scm_negative_p (SCM x)
3704 if (SCM_I_INUMP (x))
3705 return scm_from_bool (SCM_I_INUM (x) < 0);
3706 else if (SCM_BIGP (x))
3708 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
3709 scm_remember_upto_here_1 (x);
3710 return scm_from_bool (sgn < 0);
3712 else if (SCM_REALP (x))
3713 return scm_from_bool(SCM_REAL_VALUE (x) < 0.0);
3714 else if (SCM_FRACTIONP (x))
3715 return scm_negative_p (SCM_FRACTION_NUMERATOR (x));
3717 SCM_WTA_DISPATCH_1 (g_negative_p, x, SCM_ARG1, s_negative_p);
3721 /* scm_min and scm_max return an inexact when either argument is inexact, as
3722 required by r5rs. On that basis, for exact/inexact combinations the
3723 exact is converted to inexact to compare and possibly return. This is
3724 unlike scm_less_p above which takes some trouble to preserve all bits in
3725 its test, such trouble is not required for min and max. */
3727 SCM_GPROC1 (s_max, "max", scm_tc7_asubr, scm_max, g_max);
3728 /* "Return the maximum of all parameter values."
3731 scm_max (SCM x, SCM y)
3736 SCM_WTA_DISPATCH_0 (g_max, s_max);
3737 else if (SCM_I_INUMP(x) || SCM_BIGP(x) || SCM_REALP(x) || SCM_FRACTIONP(x))
3740 SCM_WTA_DISPATCH_1 (g_max, x, SCM_ARG1, s_max);
3743 if (SCM_I_INUMP (x))
3745 long xx = SCM_I_INUM (x);
3746 if (SCM_I_INUMP (y))
3748 long yy = SCM_I_INUM (y);
3749 return (xx < yy) ? y : x;
3751 else if (SCM_BIGP (y))
3753 int sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
3754 scm_remember_upto_here_1 (y);
3755 return (sgn < 0) ? x : y;
3757 else if (SCM_REALP (y))
3760 /* if y==NaN then ">" is false and we return NaN */
3761 return (z > SCM_REAL_VALUE (y)) ? scm_from_double (z) : y;
3763 else if (SCM_FRACTIONP (y))
3766 return (scm_is_false (scm_less_p (x, y)) ? x : y);
3769 SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
3771 else if (SCM_BIGP (x))
3773 if (SCM_I_INUMP (y))
3775 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
3776 scm_remember_upto_here_1 (x);
3777 return (sgn < 0) ? y : x;
3779 else if (SCM_BIGP (y))
3781 int cmp = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
3782 scm_remember_upto_here_2 (x, y);
3783 return (cmp > 0) ? x : y;
3785 else if (SCM_REALP (y))
3787 /* if y==NaN then xx>yy is false, so we return the NaN y */
3790 xx = scm_i_big2dbl (x);
3791 yy = SCM_REAL_VALUE (y);
3792 return (xx > yy ? scm_from_double (xx) : y);
3794 else if (SCM_FRACTIONP (y))
3799 SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
3801 else if (SCM_REALP (x))
3803 if (SCM_I_INUMP (y))
3805 double z = SCM_I_INUM (y);
3806 /* if x==NaN then "<" is false and we return NaN */
3807 return (SCM_REAL_VALUE (x) < z) ? scm_from_double (z) : x;
3809 else if (SCM_BIGP (y))
3814 else if (SCM_REALP (y))
3816 /* if x==NaN then our explicit check means we return NaN
3817 if y==NaN then ">" is false and we return NaN
3818 calling isnan is unavoidable, since it's the only way to know
3819 which of x or y causes any compares to be false */
3820 double xx = SCM_REAL_VALUE (x);
3821 return (xisnan (xx) || xx > SCM_REAL_VALUE (y)) ? x : y;
3823 else if (SCM_FRACTIONP (y))
3825 double yy = scm_i_fraction2double (y);
3826 double xx = SCM_REAL_VALUE (x);
3827 return (xx < yy) ? scm_from_double (yy) : x;
3830 SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
3832 else if (SCM_FRACTIONP (x))
3834 if (SCM_I_INUMP (y))
3838 else if (SCM_BIGP (y))
3842 else if (SCM_REALP (y))
3844 double xx = scm_i_fraction2double (x);
3845 return (xx < SCM_REAL_VALUE (y)) ? y : scm_from_double (xx);
3847 else if (SCM_FRACTIONP (y))
3852 SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
3855 SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARG1, s_max);
3859 SCM_GPROC1 (s_min, "min", scm_tc7_asubr, scm_min, g_min);
3860 /* "Return the minium of all parameter values."
3863 scm_min (SCM x, SCM y)
3868 SCM_WTA_DISPATCH_0 (g_min, s_min);
3869 else if (SCM_I_INUMP(x) || SCM_BIGP(x) || SCM_REALP(x) || SCM_FRACTIONP(x))
3872 SCM_WTA_DISPATCH_1 (g_min, x, SCM_ARG1, s_min);
3875 if (SCM_I_INUMP (x))
3877 long xx = SCM_I_INUM (x);
3878 if (SCM_I_INUMP (y))
3880 long yy = SCM_I_INUM (y);
3881 return (xx < yy) ? x : y;
3883 else if (SCM_BIGP (y))
3885 int sgn = mpz_sgn (SCM_I_BIG_MPZ (y));
3886 scm_remember_upto_here_1 (y);
3887 return (sgn < 0) ? y : x;
3889 else if (SCM_REALP (y))
3892 /* if y==NaN then "<" is false and we return NaN */
3893 return (z < SCM_REAL_VALUE (y)) ? scm_from_double (z) : y;
3895 else if (SCM_FRACTIONP (y))
3898 return (scm_is_false (scm_less_p (x, y)) ? y : x);
3901 SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min);
3903 else if (SCM_BIGP (x))
3905 if (SCM_I_INUMP (y))
3907 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
3908 scm_remember_upto_here_1 (x);
3909 return (sgn < 0) ? x : y;
3911 else if (SCM_BIGP (y))
3913 int cmp = mpz_cmp (SCM_I_BIG_MPZ (x), SCM_I_BIG_MPZ (y));
3914 scm_remember_upto_here_2 (x, y);
3915 return (cmp > 0) ? y : x;
3917 else if (SCM_REALP (y))
3919 /* if y==NaN then xx<yy is false, so we return the NaN y */
3922 xx = scm_i_big2dbl (x);
3923 yy = SCM_REAL_VALUE (y);
3924 return (xx < yy ? scm_from_double (xx) : y);
3926 else if (SCM_FRACTIONP (y))
3931 SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min);
3933 else if (SCM_REALP (x))
3935 if (SCM_I_INUMP (y))
3937 double z = SCM_I_INUM (y);
3938 /* if x==NaN then "<" is false and we return NaN */
3939 return (z < SCM_REAL_VALUE (x)) ? scm_from_double (z) : x;
3941 else if (SCM_BIGP (y))
3946 else if (SCM_REALP (y))
3948 /* if x==NaN then our explicit check means we return NaN
3949 if y==NaN then "<" is false and we return NaN
3950 calling isnan is unavoidable, since it's the only way to know
3951 which of x or y causes any compares to be false */
3952 double xx = SCM_REAL_VALUE (x);
3953 return (xisnan (xx) || xx < SCM_REAL_VALUE (y)) ? x : y;
3955 else if (SCM_FRACTIONP (y))
3957 double yy = scm_i_fraction2double (y);
3958 double xx = SCM_REAL_VALUE (x);
3959 return (yy < xx) ? scm_from_double (yy) : x;
3962 SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARGn, s_min);
3964 else if (SCM_FRACTIONP (x))
3966 if (SCM_I_INUMP (y))
3970 else if (SCM_BIGP (y))
3974 else if (SCM_REALP (y))
3976 double xx = scm_i_fraction2double (x);
3977 return (SCM_REAL_VALUE (y) < xx) ? y : scm_from_double (xx);
3979 else if (SCM_FRACTIONP (y))
3984 SCM_WTA_DISPATCH_2 (g_max, x, y, SCM_ARGn, s_max);
3987 SCM_WTA_DISPATCH_2 (g_min, x, y, SCM_ARG1, s_min);
3991 SCM_GPROC1 (s_sum, "+", scm_tc7_asubr, scm_sum, g_sum);
3992 /* "Return the sum of all parameter values. Return 0 if called without\n"
3996 scm_sum (SCM x, SCM y)
3998 if (SCM_UNLIKELY (SCM_UNBNDP (y)))
4000 if (SCM_NUMBERP (x)) return x;
4001 if (SCM_UNBNDP (x)) return SCM_INUM0;
4002 SCM_WTA_DISPATCH_1 (g_sum, x, SCM_ARG1, s_sum);
4005 if (SCM_LIKELY (SCM_I_INUMP (x)))
4007 if (SCM_LIKELY (SCM_I_INUMP (y)))
4009 long xx = SCM_I_INUM (x);
4010 long yy = SCM_I_INUM (y);
4011 long int z = xx + yy;
4012 return SCM_FIXABLE (z) ? SCM_I_MAKINUM (z) : scm_i_long2big (z);
4014 else if (SCM_BIGP (y))
4019 else if (SCM_REALP (y))
4021 long int xx = SCM_I_INUM (x);
4022 return scm_from_double (xx + SCM_REAL_VALUE (y));
4024 else if (SCM_COMPLEXP (y))
4026 long int xx = SCM_I_INUM (x);
4027 return scm_c_make_rectangular (xx + SCM_COMPLEX_REAL (y),
4028 SCM_COMPLEX_IMAG (y));
4030 else if (SCM_FRACTIONP (y))
4031 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y),
4032 scm_product (x, SCM_FRACTION_DENOMINATOR (y))),
4033 SCM_FRACTION_DENOMINATOR (y));
4035 SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
4036 } else if (SCM_BIGP (x))
4038 if (SCM_I_INUMP (y))
4043 inum = SCM_I_INUM (y);
4046 bigsgn = mpz_sgn (SCM_I_BIG_MPZ (x));
4049 SCM result = scm_i_mkbig ();
4050 mpz_sub_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), - inum);
4051 scm_remember_upto_here_1 (x);
4052 /* we know the result will have to be a bignum */
4055 return scm_i_normbig (result);
4059 SCM result = scm_i_mkbig ();
4060 mpz_add_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), inum);
4061 scm_remember_upto_here_1 (x);
4062 /* we know the result will have to be a bignum */
4065 return scm_i_normbig (result);
4068 else if (SCM_BIGP (y))
4070 SCM result = scm_i_mkbig ();
4071 int sgn_x = mpz_sgn (SCM_I_BIG_MPZ (x));
4072 int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y));
4073 mpz_add (SCM_I_BIG_MPZ (result),
4076 scm_remember_upto_here_2 (x, y);
4077 /* we know the result will have to be a bignum */
4080 return scm_i_normbig (result);
4082 else if (SCM_REALP (y))
4084 double result = mpz_get_d (SCM_I_BIG_MPZ (x)) + SCM_REAL_VALUE (y);
4085 scm_remember_upto_here_1 (x);
4086 return scm_from_double (result);
4088 else if (SCM_COMPLEXP (y))
4090 double real_part = (mpz_get_d (SCM_I_BIG_MPZ (x))
4091 + SCM_COMPLEX_REAL (y));
4092 scm_remember_upto_here_1 (x);
4093 return scm_c_make_rectangular (real_part, SCM_COMPLEX_IMAG (y));
4095 else if (SCM_FRACTIONP (y))
4096 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y),
4097 scm_product (x, SCM_FRACTION_DENOMINATOR (y))),
4098 SCM_FRACTION_DENOMINATOR (y));
4100 SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
4102 else if (SCM_REALP (x))
4104 if (SCM_I_INUMP (y))
4105 return scm_from_double (SCM_REAL_VALUE (x) + SCM_I_INUM (y));
4106 else if (SCM_BIGP (y))
4108 double result = mpz_get_d (SCM_I_BIG_MPZ (y)) + SCM_REAL_VALUE (x);
4109 scm_remember_upto_here_1 (y);
4110 return scm_from_double (result);
4112 else if (SCM_REALP (y))
4113 return scm_from_double (SCM_REAL_VALUE (x) + SCM_REAL_VALUE (y));
4114 else if (SCM_COMPLEXP (y))
4115 return scm_c_make_rectangular (SCM_REAL_VALUE (x) + SCM_COMPLEX_REAL (y),
4116 SCM_COMPLEX_IMAG (y));
4117 else if (SCM_FRACTIONP (y))
4118 return scm_from_double (SCM_REAL_VALUE (x) + scm_i_fraction2double (y));
4120 SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
4122 else if (SCM_COMPLEXP (x))
4124 if (SCM_I_INUMP (y))
4125 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) + SCM_I_INUM (y),
4126 SCM_COMPLEX_IMAG (x));
4127 else if (SCM_BIGP (y))
4129 double real_part = (mpz_get_d (SCM_I_BIG_MPZ (y))
4130 + SCM_COMPLEX_REAL (x));
4131 scm_remember_upto_here_1 (y);
4132 return scm_c_make_rectangular (real_part, SCM_COMPLEX_IMAG (x));
4134 else if (SCM_REALP (y))
4135 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) + SCM_REAL_VALUE (y),
4136 SCM_COMPLEX_IMAG (x));
4137 else if (SCM_COMPLEXP (y))
4138 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) + SCM_COMPLEX_REAL (y),
4139 SCM_COMPLEX_IMAG (x) + SCM_COMPLEX_IMAG (y));
4140 else if (SCM_FRACTIONP (y))
4141 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) + scm_i_fraction2double (y),
4142 SCM_COMPLEX_IMAG (x));
4144 SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
4146 else if (SCM_FRACTIONP (x))
4148 if (SCM_I_INUMP (y))
4149 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x),
4150 scm_product (y, SCM_FRACTION_DENOMINATOR (x))),
4151 SCM_FRACTION_DENOMINATOR (x));
4152 else if (SCM_BIGP (y))
4153 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x),
4154 scm_product (y, SCM_FRACTION_DENOMINATOR (x))),
4155 SCM_FRACTION_DENOMINATOR (x));
4156 else if (SCM_REALP (y))
4157 return scm_from_double (SCM_REAL_VALUE (y) + scm_i_fraction2double (x));
4158 else if (SCM_COMPLEXP (y))
4159 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y) + scm_i_fraction2double (x),
4160 SCM_COMPLEX_IMAG (y));
4161 else if (SCM_FRACTIONP (y))
4162 /* a/b + c/d = (ad + bc) / bd */
4163 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x), SCM_FRACTION_DENOMINATOR (y)),
4164 scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x))),
4165 scm_product (SCM_FRACTION_DENOMINATOR (x), SCM_FRACTION_DENOMINATOR (y)));
4167 SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARGn, s_sum);
4170 SCM_WTA_DISPATCH_2 (g_sum, x, y, SCM_ARG1, s_sum);
4174 SCM_DEFINE (scm_oneplus, "1+", 1, 0, 0,
4176 "Return @math{@var{x}+1}.")
4177 #define FUNC_NAME s_scm_oneplus
4179 return scm_sum (x, SCM_I_MAKINUM (1));
4184 SCM_GPROC1 (s_difference, "-", scm_tc7_asubr, scm_difference, g_difference);
4185 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4186 * the sum of all but the first argument are subtracted from the first
4188 #define FUNC_NAME s_difference
4190 scm_difference (SCM x, SCM y)
4192 if (SCM_UNLIKELY (SCM_UNBNDP (y)))
4195 SCM_WTA_DISPATCH_0 (g_difference, s_difference);
4197 if (SCM_I_INUMP (x))
4199 long xx = -SCM_I_INUM (x);
4200 if (SCM_FIXABLE (xx))
4201 return SCM_I_MAKINUM (xx);
4203 return scm_i_long2big (xx);
4205 else if (SCM_BIGP (x))
4206 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4207 bignum, but negating that gives a fixnum. */
4208 return scm_i_normbig (scm_i_clonebig (x, 0));
4209 else if (SCM_REALP (x))
4210 return scm_from_double (-SCM_REAL_VALUE (x));
4211 else if (SCM_COMPLEXP (x))
4212 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x),
4213 -SCM_COMPLEX_IMAG (x));
4214 else if (SCM_FRACTIONP (x))
4215 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED),
4216 SCM_FRACTION_DENOMINATOR (x));
4218 SCM_WTA_DISPATCH_1 (g_difference, x, SCM_ARG1, s_difference);
4221 if (SCM_LIKELY (SCM_I_INUMP (x)))
4223 if (SCM_LIKELY (SCM_I_INUMP (y)))
4225 long int xx = SCM_I_INUM (x);
4226 long int yy = SCM_I_INUM (y);
4227 long int z = xx - yy;
4228 if (SCM_FIXABLE (z))
4229 return SCM_I_MAKINUM (z);
4231 return scm_i_long2big (z);
4233 else if (SCM_BIGP (y))
4235 /* inum-x - big-y */
4236 long xx = SCM_I_INUM (x);
4239 return scm_i_clonebig (y, 0);
4242 int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y));
4243 SCM result = scm_i_mkbig ();
4246 mpz_ui_sub (SCM_I_BIG_MPZ (result), xx, SCM_I_BIG_MPZ (y));
4249 /* x - y == -(y + -x) */
4250 mpz_add_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (y), -xx);
4251 mpz_neg (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result));
4253 scm_remember_upto_here_1 (y);
4255 if ((xx < 0 && (sgn_y > 0)) || ((xx > 0) && sgn_y < 0))
4256 /* we know the result will have to be a bignum */
4259 return scm_i_normbig (result);
4262 else if (SCM_REALP (y))
4264 long int xx = SCM_I_INUM (x);
4265 return scm_from_double (xx - SCM_REAL_VALUE (y));
4267 else if (SCM_COMPLEXP (y))
4269 long int xx = SCM_I_INUM (x);
4270 return scm_c_make_rectangular (xx - SCM_COMPLEX_REAL (y),
4271 - SCM_COMPLEX_IMAG (y));
4273 else if (SCM_FRACTIONP (y))
4274 /* a - b/c = (ac - b) / c */
4275 return scm_i_make_ratio (scm_difference (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
4276 SCM_FRACTION_NUMERATOR (y)),
4277 SCM_FRACTION_DENOMINATOR (y));
4279 SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
4281 else if (SCM_BIGP (x))
4283 if (SCM_I_INUMP (y))
4285 /* big-x - inum-y */
4286 long yy = SCM_I_INUM (y);
4287 int sgn_x = mpz_sgn (SCM_I_BIG_MPZ (x));
4289 scm_remember_upto_here_1 (x);
4291 return (SCM_FIXABLE (-yy) ?
4292 SCM_I_MAKINUM (-yy) : scm_from_long (-yy));
4295 SCM result = scm_i_mkbig ();
4298 mpz_sub_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), yy);
4300 mpz_add_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), -yy);
4301 scm_remember_upto_here_1 (x);
4303 if ((sgn_x < 0 && (yy > 0)) || ((sgn_x > 0) && yy < 0))
4304 /* we know the result will have to be a bignum */
4307 return scm_i_normbig (result);
4310 else if (SCM_BIGP (y))
4312 int sgn_x = mpz_sgn (SCM_I_BIG_MPZ (x));
4313 int sgn_y = mpz_sgn (SCM_I_BIG_MPZ (y));
4314 SCM result = scm_i_mkbig ();
4315 mpz_sub (SCM_I_BIG_MPZ (result),
4318 scm_remember_upto_here_2 (x, y);
4319 /* we know the result will have to be a bignum */
4320 if ((sgn_x == 1) && (sgn_y == -1))
4322 if ((sgn_x == -1) && (sgn_y == 1))
4324 return scm_i_normbig (result);
4326 else if (SCM_REALP (y))
4328 double result = mpz_get_d (SCM_I_BIG_MPZ (x)) - SCM_REAL_VALUE (y);
4329 scm_remember_upto_here_1 (x);
4330 return scm_from_double (result);
4332 else if (SCM_COMPLEXP (y))
4334 double real_part = (mpz_get_d (SCM_I_BIG_MPZ (x))
4335 - SCM_COMPLEX_REAL (y));
4336 scm_remember_upto_here_1 (x);
4337 return scm_c_make_rectangular (real_part, - SCM_COMPLEX_IMAG (y));
4339 else if (SCM_FRACTIONP (y))
4340 return scm_i_make_ratio (scm_difference (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
4341 SCM_FRACTION_NUMERATOR (y)),
4342 SCM_FRACTION_DENOMINATOR (y));
4343 else SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
4345 else if (SCM_REALP (x))
4347 if (SCM_I_INUMP (y))
4348 return scm_from_double (SCM_REAL_VALUE (x) - SCM_I_INUM (y));
4349 else if (SCM_BIGP (y))
4351 double result = SCM_REAL_VALUE (x) - mpz_get_d (SCM_I_BIG_MPZ (y));
4352 scm_remember_upto_here_1 (x);
4353 return scm_from_double (result);
4355 else if (SCM_REALP (y))
4356 return scm_from_double (SCM_REAL_VALUE (x) - SCM_REAL_VALUE (y));
4357 else if (SCM_COMPLEXP (y))
4358 return scm_c_make_rectangular (SCM_REAL_VALUE (x) - SCM_COMPLEX_REAL (y),
4359 -SCM_COMPLEX_IMAG (y));
4360 else if (SCM_FRACTIONP (y))
4361 return scm_from_double (SCM_REAL_VALUE (x) - scm_i_fraction2double (y));
4363 SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
4365 else if (SCM_COMPLEXP (x))
4367 if (SCM_I_INUMP (y))
4368 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) - SCM_I_INUM (y),
4369 SCM_COMPLEX_IMAG (x));
4370 else if (SCM_BIGP (y))
4372 double real_part = (SCM_COMPLEX_REAL (x)
4373 - mpz_get_d (SCM_I_BIG_MPZ (y)));
4374 scm_remember_upto_here_1 (x);
4375 return scm_c_make_rectangular (real_part, SCM_COMPLEX_IMAG (y));
4377 else if (SCM_REALP (y))
4378 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) - SCM_REAL_VALUE (y),
4379 SCM_COMPLEX_IMAG (x));
4380 else if (SCM_COMPLEXP (y))
4381 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) - SCM_COMPLEX_REAL (y),
4382 SCM_COMPLEX_IMAG (x) - SCM_COMPLEX_IMAG (y));
4383 else if (SCM_FRACTIONP (y))
4384 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) - scm_i_fraction2double (y),
4385 SCM_COMPLEX_IMAG (x));
4387 SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
4389 else if (SCM_FRACTIONP (x))
4391 if (SCM_I_INUMP (y))
4392 /* a/b - c = (a - cb) / b */
4393 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x),
4394 scm_product(y, SCM_FRACTION_DENOMINATOR (x))),
4395 SCM_FRACTION_DENOMINATOR (x));
4396 else if (SCM_BIGP (y))
4397 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x),
4398 scm_product(y, SCM_FRACTION_DENOMINATOR (x))),
4399 SCM_FRACTION_DENOMINATOR (x));
4400 else if (SCM_REALP (y))
4401 return scm_from_double (scm_i_fraction2double (x) - SCM_REAL_VALUE (y));
4402 else if (SCM_COMPLEXP (y))
4403 return scm_c_make_rectangular (scm_i_fraction2double (x) - SCM_COMPLEX_REAL (y),
4404 -SCM_COMPLEX_IMAG (y));
4405 else if (SCM_FRACTIONP (y))
4406 /* a/b - c/d = (ad - bc) / bd */
4407 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x), SCM_FRACTION_DENOMINATOR (y)),
4408 scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x))),
4409 scm_product (SCM_FRACTION_DENOMINATOR (x), SCM_FRACTION_DENOMINATOR (y)));
4411 SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARGn, s_difference);
4414 SCM_WTA_DISPATCH_2 (g_difference, x, y, SCM_ARG1, s_difference);
4419 SCM_DEFINE (scm_oneminus, "1-", 1, 0, 0,
4421 "Return @math{@var{x}-1}.")
4422 #define FUNC_NAME s_scm_oneminus
4424 return scm_difference (x, SCM_I_MAKINUM (1));
4429 SCM_GPROC1 (s_product, "*", scm_tc7_asubr, scm_product, g_product);
4430 /* "Return the product of all arguments. If called without arguments,\n"
4434 scm_product (SCM x, SCM y)
4436 if (SCM_UNLIKELY (SCM_UNBNDP (y)))
4439 return SCM_I_MAKINUM (1L);
4440 else if (SCM_NUMBERP (x))
4443 SCM_WTA_DISPATCH_1 (g_product, x, SCM_ARG1, s_product);
4446 if (SCM_LIKELY (SCM_I_INUMP (x)))
4451 xx = SCM_I_INUM (x);
4455 case 0: return x; break;
4456 case 1: return y; break;
4459 if (SCM_LIKELY (SCM_I_INUMP (y)))
4461 long yy = SCM_I_INUM (y);
4463 SCM k = SCM_I_MAKINUM (kk);
4464 if ((kk == SCM_I_INUM (k)) && (kk / xx == yy))
4468 SCM result = scm_i_long2big (xx);
4469 mpz_mul_si (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result), yy);
4470 return scm_i_normbig (result);
4473 else if (SCM_BIGP (y))
4475 SCM result = scm_i_mkbig ();
4476 mpz_mul_si (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (y), xx);
4477 scm_remember_upto_here_1 (y);
4480 else if (SCM_REALP (y))
4481 return scm_from_double (xx * SCM_REAL_VALUE (y));
4482 else if (SCM_COMPLEXP (y))
4483 return scm_c_make_rectangular (xx * SCM_COMPLEX_REAL (y),
4484 xx * SCM_COMPLEX_IMAG (y));
4485 else if (SCM_FRACTIONP (y))
4486 return scm_i_make_ratio (scm_product (x, SCM_FRACTION_NUMERATOR (y)),
4487 SCM_FRACTION_DENOMINATOR (y));
4489 SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
4491 else if (SCM_BIGP (x))
4493 if (SCM_I_INUMP (y))
4498 else if (SCM_BIGP (y))
4500 SCM result = scm_i_mkbig ();
4501 mpz_mul (SCM_I_BIG_MPZ (result),
4504 scm_remember_upto_here_2 (x, y);
4507 else if (SCM_REALP (y))
4509 double result = mpz_get_d (SCM_I_BIG_MPZ (x)) * SCM_REAL_VALUE (y);
4510 scm_remember_upto_here_1 (x);
4511 return scm_from_double (result);
4513 else if (SCM_COMPLEXP (y))
4515 double z = mpz_get_d (SCM_I_BIG_MPZ (x));
4516 scm_remember_upto_here_1 (x);
4517 return scm_c_make_rectangular (z * SCM_COMPLEX_REAL (y),
4518 z * SCM_COMPLEX_IMAG (y));
4520 else if (SCM_FRACTIONP (y))
4521 return scm_i_make_ratio (scm_product (x, SCM_FRACTION_NUMERATOR (y)),
4522 SCM_FRACTION_DENOMINATOR (y));
4524 SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
4526 else if (SCM_REALP (x))
4528 if (SCM_I_INUMP (y))
4530 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4531 if (scm_is_eq (y, SCM_INUM0))
4533 return scm_from_double (SCM_I_INUM (y) * SCM_REAL_VALUE (x));
4535 else if (SCM_BIGP (y))
4537 double result = mpz_get_d (SCM_I_BIG_MPZ (y)) * SCM_REAL_VALUE (x);
4538 scm_remember_upto_here_1 (y);
4539 return scm_from_double (result);
4541 else if (SCM_REALP (y))
4542 return scm_from_double (SCM_REAL_VALUE (x) * SCM_REAL_VALUE (y));
4543 else if (SCM_COMPLEXP (y))
4544 return scm_c_make_rectangular (SCM_REAL_VALUE (x) * SCM_COMPLEX_REAL (y),
4545 SCM_REAL_VALUE (x) * SCM_COMPLEX_IMAG (y));
4546 else if (SCM_FRACTIONP (y))
4547 return scm_from_double (SCM_REAL_VALUE (x) * scm_i_fraction2double (y));
4549 SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
4551 else if (SCM_COMPLEXP (x))
4553 if (SCM_I_INUMP (y))
4555 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4556 if (scm_is_eq (y, SCM_INUM0))
4558 return scm_c_make_rectangular (SCM_I_INUM (y) * SCM_COMPLEX_REAL (x),
4559 SCM_I_INUM (y) * SCM_COMPLEX_IMAG (x));
4561 else if (SCM_BIGP (y))
4563 double z = mpz_get_d (SCM_I_BIG_MPZ (y));
4564 scm_remember_upto_here_1 (y);
4565 return scm_c_make_rectangular (z * SCM_COMPLEX_REAL (x),
4566 z * SCM_COMPLEX_IMAG (x));
4568 else if (SCM_REALP (y))
4569 return scm_c_make_rectangular (SCM_REAL_VALUE (y) * SCM_COMPLEX_REAL (x),
4570 SCM_REAL_VALUE (y) * SCM_COMPLEX_IMAG (x));
4571 else if (SCM_COMPLEXP (y))
4573 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x) * SCM_COMPLEX_REAL (y)
4574 - SCM_COMPLEX_IMAG (x) * SCM_COMPLEX_IMAG (y),
4575 SCM_COMPLEX_REAL (x) * SCM_COMPLEX_IMAG (y)
4576 + SCM_COMPLEX_IMAG (x) * SCM_COMPLEX_REAL (y));
4578 else if (SCM_FRACTIONP (y))
4580 double yy = scm_i_fraction2double (y);
4581 return scm_c_make_rectangular (yy * SCM_COMPLEX_REAL (x),
4582 yy * SCM_COMPLEX_IMAG (x));
4585 SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
4587 else if (SCM_FRACTIONP (x))
4589 if (SCM_I_INUMP (y))
4590 return scm_i_make_ratio (scm_product (y, SCM_FRACTION_NUMERATOR (x)),
4591 SCM_FRACTION_DENOMINATOR (x));
4592 else if (SCM_BIGP (y))
4593 return scm_i_make_ratio (scm_product (y, SCM_FRACTION_NUMERATOR (x)),
4594 SCM_FRACTION_DENOMINATOR (x));
4595 else if (SCM_REALP (y))
4596 return scm_from_double (scm_i_fraction2double (x) * SCM_REAL_VALUE (y));
4597 else if (SCM_COMPLEXP (y))
4599 double xx = scm_i_fraction2double (x);
4600 return scm_c_make_rectangular (xx * SCM_COMPLEX_REAL (y),
4601 xx * SCM_COMPLEX_IMAG (y));
4603 else if (SCM_FRACTIONP (y))
4604 /* a/b * c/d = ac / bd */
4605 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x),
4606 SCM_FRACTION_NUMERATOR (y)),
4607 scm_product (SCM_FRACTION_DENOMINATOR (x),
4608 SCM_FRACTION_DENOMINATOR (y)));
4610 SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARGn, s_product);
4613 SCM_WTA_DISPATCH_2 (g_product, x, y, SCM_ARG1, s_product);
4616 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4617 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4618 #define ALLOW_DIVIDE_BY_ZERO
4619 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4622 /* The code below for complex division is adapted from the GNU
4623 libstdc++, which adapted it from f2c's libF77, and is subject to
4626 /****************************************************************
4627 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4629 Permission to use, copy, modify, and distribute this software
4630 and its documentation for any purpose and without fee is hereby
4631 granted, provided that the above copyright notice appear in all
4632 copies and that both that the copyright notice and this
4633 permission notice and warranty disclaimer appear in supporting
4634 documentation, and that the names of AT&T Bell Laboratories or
4635 Bellcore or any of their entities not be used in advertising or
4636 publicity pertaining to distribution of the software without
4637 specific, written prior permission.
4639 AT&T and Bellcore disclaim all warranties with regard to this
4640 software, including all implied warranties of merchantability
4641 and fitness. In no event shall AT&T or Bellcore be liable for
4642 any special, indirect or consequential damages or any damages
4643 whatsoever resulting from loss of use, data or profits, whether
4644 in an action of contract, negligence or other tortious action,
4645 arising out of or in connection with the use or performance of
4647 ****************************************************************/
4649 SCM_GPROC1 (s_divide, "/", scm_tc7_asubr, scm_divide, g_divide);
4650 /* Divide the first argument by the product of the remaining
4651 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4653 #define FUNC_NAME s_divide
4655 scm_i_divide (SCM x, SCM y, int inexact)
4659 if (SCM_UNLIKELY (SCM_UNBNDP (y)))
4662 SCM_WTA_DISPATCH_0 (g_divide, s_divide);
4663 else if (SCM_I_INUMP (x))
4665 long xx = SCM_I_INUM (x);
4666 if (xx == 1 || xx == -1)
4668 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4670 scm_num_overflow (s_divide);
4675 return scm_from_double (1.0 / (double) xx);
4676 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x);
4679 else if (SCM_BIGP (x))
4682 return scm_from_double (1.0 / scm_i_big2dbl (x));
4683 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x);
4685 else if (SCM_REALP (x))
4687 double xx = SCM_REAL_VALUE (x);
4688 #ifndef ALLOW_DIVIDE_BY_ZERO
4690 scm_num_overflow (s_divide);
4693 return scm_from_double (1.0 / xx);
4695 else if (SCM_COMPLEXP (x))
4697 double r = SCM_COMPLEX_REAL (x);
4698 double i = SCM_COMPLEX_IMAG (x);
4699 if (fabs(r) <= fabs(i))
4702 double d = i * (1.0 + t * t);
4703 return scm_c_make_rectangular (t / d, -1.0 / d);
4708 double d = r * (1.0 + t * t);
4709 return scm_c_make_rectangular (1.0 / d, -t / d);
4712 else if (SCM_FRACTIONP (x))
4713 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x),
4714 SCM_FRACTION_NUMERATOR (x));
4716 SCM_WTA_DISPATCH_1 (g_divide, x, SCM_ARG1, s_divide);
4719 if (SCM_LIKELY (SCM_I_INUMP (x)))
4721 long xx = SCM_I_INUM (x);
4722 if (SCM_LIKELY (SCM_I_INUMP (y)))
4724 long yy = SCM_I_INUM (y);
4727 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4728 scm_num_overflow (s_divide);
4730 return scm_from_double ((double) xx / (double) yy);
4733 else if (xx % yy != 0)
4736 return scm_from_double ((double) xx / (double) yy);
4737 else return scm_i_make_ratio (x, y);
4742 if (SCM_FIXABLE (z))
4743 return SCM_I_MAKINUM (z);
4745 return scm_i_long2big (z);
4748 else if (SCM_BIGP (y))
4751 return scm_from_double ((double) xx / scm_i_big2dbl (y));
4752 else return scm_i_make_ratio (x, y);
4754 else if (SCM_REALP (y))
4756 double yy = SCM_REAL_VALUE (y);
4757 #ifndef ALLOW_DIVIDE_BY_ZERO
4759 scm_num_overflow (s_divide);
4762 return scm_from_double ((double) xx / yy);
4764 else if (SCM_COMPLEXP (y))
4767 complex_div: /* y _must_ be a complex number */
4769 double r = SCM_COMPLEX_REAL (y);
4770 double i = SCM_COMPLEX_IMAG (y);
4771 if (fabs(r) <= fabs(i))
4774 double d = i * (1.0 + t * t);
4775 return scm_c_make_rectangular ((a * t) / d, -a / d);
4780 double d = r * (1.0 + t * t);
4781 return scm_c_make_rectangular (a / d, -(a * t) / d);
4785 else if (SCM_FRACTIONP (y))
4786 /* a / b/c = ac / b */
4787 return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
4788 SCM_FRACTION_NUMERATOR (y));
4790 SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
4792 else if (SCM_BIGP (x))
4794 if (SCM_I_INUMP (y))
4796 long int yy = SCM_I_INUM (y);
4799 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4800 scm_num_overflow (s_divide);
4802 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
4803 scm_remember_upto_here_1 (x);
4804 return (sgn == 0) ? scm_nan () : scm_inf ();
4811 /* FIXME: HMM, what are the relative performance issues here?
4812 We need to test. Is it faster on average to test
4813 divisible_p, then perform whichever operation, or is it
4814 faster to perform the integer div opportunistically and
4815 switch to real if there's a remainder? For now we take the
4816 middle ground: test, then if divisible, use the faster div
4819 long abs_yy = yy < 0 ? -yy : yy;
4820 int divisible_p = mpz_divisible_ui_p (SCM_I_BIG_MPZ (x), abs_yy);
4824 SCM result = scm_i_mkbig ();
4825 mpz_divexact_ui (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (x), abs_yy);
4826 scm_remember_upto_here_1 (x);
4828 mpz_neg (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result));
4829 return scm_i_normbig (result);
4834 return scm_from_double (scm_i_big2dbl (x) / (double) yy);
4835 else return scm_i_make_ratio (x, y);
4839 else if (SCM_BIGP (y))
4841 int y_is_zero = (mpz_sgn (SCM_I_BIG_MPZ (y)) == 0);
4844 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4845 scm_num_overflow (s_divide);
4847 int sgn = mpz_sgn (SCM_I_BIG_MPZ (x));
4848 scm_remember_upto_here_1 (x);
4849 return (sgn == 0) ? scm_nan () : scm_inf ();
4857 /* It's easily possible for the ratio x/y to fit a double
4858 but one or both x and y be too big to fit a double,
4859 hence the use of mpq_get_d rather than converting and
4862 *mpq_numref(q) = *SCM_I_BIG_MPZ (x);
4863 *mpq_denref(q) = *SCM_I_BIG_MPZ (y);
4864 return scm_from_double (mpq_get_d (q));
4868 int divisible_p = mpz_divisible_p (SCM_I_BIG_MPZ (x),
4872 SCM result = scm_i_mkbig ();
4873 mpz_divexact (SCM_I_BIG_MPZ (result),
4876 scm_remember_upto_here_2 (x, y);
4877 return scm_i_normbig (result);
4880 return scm_i_make_ratio (x, y);
4884 else if (SCM_REALP (y))
4886 double yy = SCM_REAL_VALUE (y);
4887 #ifndef ALLOW_DIVIDE_BY_ZERO
4889 scm_num_overflow (s_divide);
4892 return scm_from_double (scm_i_big2dbl (x) / yy);
4894 else if (SCM_COMPLEXP (y))
4896 a = scm_i_big2dbl (x);
4899 else if (SCM_FRACTIONP (y))
4900 return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
4901 SCM_FRACTION_NUMERATOR (y));
4903 SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
4905 else if (SCM_REALP (x))
4907 double rx = SCM_REAL_VALUE (x);
4908 if (SCM_I_INUMP (y))
4910 long int yy = SCM_I_INUM (y);
4911 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4913 scm_num_overflow (s_divide);
4916 return scm_from_double (rx / (double) yy);
4918 else if (SCM_BIGP (y))
4920 double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
4921 scm_remember_upto_here_1 (y);
4922 return scm_from_double (rx / dby);
4924 else if (SCM_REALP (y))
4926 double yy = SCM_REAL_VALUE (y);
4927 #ifndef ALLOW_DIVIDE_BY_ZERO
4929 scm_num_overflow (s_divide);
4932 return scm_from_double (rx / yy);
4934 else if (SCM_COMPLEXP (y))
4939 else if (SCM_FRACTIONP (y))
4940 return scm_from_double (rx / scm_i_fraction2double (y));
4942 SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
4944 else if (SCM_COMPLEXP (x))
4946 double rx = SCM_COMPLEX_REAL (x);
4947 double ix = SCM_COMPLEX_IMAG (x);
4948 if (SCM_I_INUMP (y))
4950 long int yy = SCM_I_INUM (y);
4951 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4953 scm_num_overflow (s_divide);
4958 return scm_c_make_rectangular (rx / d, ix / d);
4961 else if (SCM_BIGP (y))
4963 double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
4964 scm_remember_upto_here_1 (y);
4965 return scm_c_make_rectangular (rx / dby, ix / dby);
4967 else if (SCM_REALP (y))
4969 double yy = SCM_REAL_VALUE (y);
4970 #ifndef ALLOW_DIVIDE_BY_ZERO
4972 scm_num_overflow (s_divide);
4975 return scm_c_make_rectangular (rx / yy, ix / yy);
4977 else if (SCM_COMPLEXP (y))
4979 double ry = SCM_COMPLEX_REAL (y);
4980 double iy = SCM_COMPLEX_IMAG (y);
4981 if (fabs(ry) <= fabs(iy))
4984 double d = iy * (1.0 + t * t);
4985 return scm_c_make_rectangular ((rx * t + ix) / d, (ix * t - rx) / d);
4990 double d = ry * (1.0 + t * t);
4991 return scm_c_make_rectangular ((rx + ix * t) / d, (ix - rx * t) / d);
4994 else if (SCM_FRACTIONP (y))
4996 double yy = scm_i_fraction2double (y);
4997 return scm_c_make_rectangular (rx / yy, ix / yy);
5000 SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
5002 else if (SCM_FRACTIONP (x))
5004 if (SCM_I_INUMP (y))
5006 long int yy = SCM_I_INUM (y);
5007 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5009 scm_num_overflow (s_divide);
5012 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x),
5013 scm_product (SCM_FRACTION_DENOMINATOR (x), y));
5015 else if (SCM_BIGP (y))
5017 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x),
5018 scm_product (SCM_FRACTION_DENOMINATOR (x), y));
5020 else if (SCM_REALP (y))
5022 double yy = SCM_REAL_VALUE (y);
5023 #ifndef ALLOW_DIVIDE_BY_ZERO
5025 scm_num_overflow (s_divide);
5028 return scm_from_double (scm_i_fraction2double (x) / yy);
5030 else if (SCM_COMPLEXP (y))
5032 a = scm_i_fraction2double (x);
5035 else if (SCM_FRACTIONP (y))
5036 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x), SCM_FRACTION_DENOMINATOR (y)),
5037 scm_product (SCM_FRACTION_NUMERATOR (y), SCM_FRACTION_DENOMINATOR (x)));
5039 SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
5042 SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARG1, s_divide);
5046 scm_divide (SCM x, SCM y)
5048 return scm_i_divide (x, y, 0);
5051 static SCM scm_divide2real (SCM x, SCM y)
5053 return scm_i_divide (x, y, 1);
5059 scm_asinh (double x)
5064 #define asinh scm_asinh
5065 return log (x + sqrt (x * x + 1));
5068 SCM_GPROC1 (s_asinh, "$asinh", scm_tc7_dsubr, (SCM (*)()) asinh, g_asinh);
5069 /* "Return the inverse hyperbolic sine of @var{x}."
5074 scm_acosh (double x)
5079 #define acosh scm_acosh
5080 return log (x + sqrt (x * x - 1));
5083 SCM_GPROC1 (s_acosh, "$acosh", scm_tc7_dsubr, (SCM (*)()) acosh, g_acosh);
5084 /* "Return the inverse hyperbolic cosine of @var{x}."
5089 scm_atanh (double x)
5094 #define atanh scm_atanh
5095 return 0.5 * log ((1 + x) / (1 - x));
5098 SCM_GPROC1 (s_atanh, "$atanh", scm_tc7_dsubr, (SCM (*)()) atanh, g_atanh);
5099 /* "Return the inverse hyperbolic tangent of @var{x}."
5104 scm_c_truncate (double x)
5115 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5116 half-way case (ie. when x is an integer plus 0.5) going upwards.
5117 Then half-way cases are identified and adjusted down if the
5118 round-upwards didn't give the desired even integer.
5120 "plus_half == result" identifies a half-way case. If plus_half, which is
5121 x + 0.5, is an integer then x must be an integer plus 0.5.
5123 An odd "result" value is identified with result/2 != floor(result/2).
5124 This is done with plus_half, since that value is ready for use sooner in
5125 a pipelined cpu, and we're already requiring plus_half == result.
5127 Note however that we need to be careful when x is big and already an
5128 integer. In that case "x+0.5" may round to an adjacent integer, causing
5129 us to return such a value, incorrectly. For instance if the hardware is
5130 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5131 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5132 returned. Or if the hardware is in round-upwards mode, then other bigger
5133 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5134 representable value, 2^128+2^76 (or whatever), again incorrect.
5136 These bad roundings of x+0.5 are avoided by testing at the start whether
5137 x is already an integer. If it is then clearly that's the desired result
5138 already. And if it's not then the exponent must be small enough to allow
5139 an 0.5 to be represented, and hence added without a bad rounding. */
5142 scm_c_round (double x)
5144 double plus_half, result;
5149 plus_half = x + 0.5;
5150 result = floor (plus_half);
5151 /* Adjust so that the rounding is towards even. */
5152 return ((plus_half == result && plus_half / 2 != floor (plus_half / 2))
5157 SCM_DEFINE (scm_truncate_number, "truncate", 1, 0, 0,
5159 "Round the number @var{x} towards zero.")
5160 #define FUNC_NAME s_scm_truncate_number
5162 if (scm_is_false (scm_negative_p (x)))
5163 return scm_floor (x);
5165 return scm_ceiling (x);
5169 static SCM exactly_one_half;
5171 SCM_DEFINE (scm_round_number, "round", 1, 0, 0,
5173 "Round the number @var{x} towards the nearest integer. "
5174 "When it is exactly halfway between two integers, "
5175 "round towards the even one.")
5176 #define FUNC_NAME s_scm_round_number
5178 if (SCM_I_INUMP (x) || SCM_BIGP (x))
5180 else if (SCM_REALP (x))
5181 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x)));
5184 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5185 single quotient+remainder division then examining to see which way
5186 the rounding should go. */
5187 SCM plus_half = scm_sum (x, exactly_one_half);
5188 SCM result = scm_floor (plus_half);
5189 /* Adjust so that the rounding is towards even. */
5190 if (scm_is_true (scm_num_eq_p (plus_half, result))
5191 && scm_is_true (scm_odd_p (result)))
5192 return scm_difference (result, SCM_I_MAKINUM (1));
5199 SCM_PRIMITIVE_GENERIC (scm_floor, "floor", 1, 0, 0,
5201 "Round the number @var{x} towards minus infinity.")
5202 #define FUNC_NAME s_scm_floor
5204 if (SCM_I_INUMP (x) || SCM_BIGP (x))
5206 else if (SCM_REALP (x))
5207 return scm_from_double (floor (SCM_REAL_VALUE (x)));
5208 else if (SCM_FRACTIONP (x))
5210 SCM q = scm_quotient (SCM_FRACTION_NUMERATOR (x),
5211 SCM_FRACTION_DENOMINATOR (x));
5212 if (scm_is_false (scm_negative_p (x)))
5214 /* For positive x, rounding towards zero is correct. */
5219 /* For negative x, we need to return q-1 unless x is an
5220 integer. But fractions are never integer, per our
5222 return scm_difference (q, SCM_I_MAKINUM (1));
5226 SCM_WTA_DISPATCH_1 (g_scm_floor, x, 1, s_scm_floor);
5230 SCM_PRIMITIVE_GENERIC (scm_ceiling, "ceiling", 1, 0, 0,
5232 "Round the number @var{x} towards infinity.")
5233 #define FUNC_NAME s_scm_ceiling
5235 if (SCM_I_INUMP (x) || SCM_BIGP (x))
5237 else if (SCM_REALP (x))
5238 return scm_from_double (ceil (SCM_REAL_VALUE (x)));
5239 else if (SCM_FRACTIONP (x))
5241 SCM q = scm_quotient (SCM_FRACTION_NUMERATOR (x),
5242 SCM_FRACTION_DENOMINATOR (x));
5243 if (scm_is_false (scm_positive_p (x)))
5245 /* For negative x, rounding towards zero is correct. */
5250 /* For positive x, we need to return q+1 unless x is an
5251 integer. But fractions are never integer, per our
5253 return scm_sum (q, SCM_I_MAKINUM (1));
5257 SCM_WTA_DISPATCH_1 (g_scm_ceiling, x, 1, s_scm_ceiling);
5261 SCM_GPROC1 (s_i_sqrt, "$sqrt", scm_tc7_dsubr, (SCM (*)()) sqrt, g_i_sqrt);
5262 /* "Return the square root of the real number @var{x}."
5264 SCM_GPROC1 (s_i_abs, "$abs", scm_tc7_dsubr, (SCM (*)()) fabs, g_i_abs);
5265 /* "Return the absolute value of the real number @var{x}."
5267 SCM_GPROC1 (s_i_exp, "$exp", scm_tc7_dsubr, (SCM (*)()) exp, g_i_exp);
5268 /* "Return the @var{x}th power of e."
5270 SCM_GPROC1 (s_i_log, "$log", scm_tc7_dsubr, (SCM (*)()) log, g_i_log);
5271 /* "Return the natural logarithm of the real number @var{x}."
5273 SCM_GPROC1 (s_i_sin, "$sin", scm_tc7_dsubr, (SCM (*)()) sin, g_i_sin);
5274 /* "Return the sine of the real number @var{x}."
5276 SCM_GPROC1 (s_i_cos, "$cos", scm_tc7_dsubr, (SCM (*)()) cos, g_i_cos);
5277 /* "Return the cosine of the real number @var{x}."
5279 SCM_GPROC1 (s_i_tan, "$tan", scm_tc7_dsubr, (SCM (*)()) tan, g_i_tan);
5280 /* "Return the tangent of the real number @var{x}."
5282 SCM_GPROC1 (s_i_asin, "$asin", scm_tc7_dsubr, (SCM (*)()) asin, g_i_asin);
5283 /* "Return the arc sine of the real number @var{x}."
5285 SCM_GPROC1 (s_i_acos, "$acos", scm_tc7_dsubr, (SCM (*)()) acos, g_i_acos);
5286 /* "Return the arc cosine of the real number @var{x}."
5288 SCM_GPROC1 (s_i_atan, "$atan", scm_tc7_dsubr, (SCM (*)()) atan, g_i_atan);
5289 /* "Return the arc tangent of the real number @var{x}."
5291 SCM_GPROC1 (s_i_sinh, "$sinh", scm_tc7_dsubr, (SCM (*)()) sinh, g_i_sinh);
5292 /* "Return the hyperbolic sine of the real number @var{x}."
5294 SCM_GPROC1 (s_i_cosh, "$cosh", scm_tc7_dsubr, (SCM (*)()) cosh, g_i_cosh);
5295 /* "Return the hyperbolic cosine of the real number @var{x}."
5297 SCM_GPROC1 (s_i_tanh, "$tanh", scm_tc7_dsubr, (SCM (*)()) tanh, g_i_tanh);
5298 /* "Return the hyperbolic tangent of the real number @var{x}."
5306 static void scm_two_doubles (SCM x,
5308 const char *sstring,
5312 scm_two_doubles (SCM x, SCM y, const char *sstring, struct dpair *xy)
5314 if (SCM_I_INUMP (x))
5315 xy->x = SCM_I_INUM (x);
5316 else if (SCM_BIGP (x))
5317 xy->x = scm_i_big2dbl (x);
5318 else if (SCM_REALP (x))
5319 xy->x = SCM_REAL_VALUE (x);
5320 else if (SCM_FRACTIONP (x))
5321 xy->x = scm_i_fraction2double (x);
5323 scm_wrong_type_arg (sstring, SCM_ARG1, x);
5325 if (SCM_I_INUMP (y))
5326 xy->y = SCM_I_INUM (y);
5327 else if (SCM_BIGP (y))
5328 xy->y = scm_i_big2dbl (y);
5329 else if (SCM_REALP (y))
5330 xy->y = SCM_REAL_VALUE (y);
5331 else if (SCM_FRACTIONP (y))
5332 xy->y = scm_i_fraction2double (y);
5334 scm_wrong_type_arg (sstring, SCM_ARG2, y);
5338 SCM_DEFINE (scm_sys_expt, "$expt", 2, 0, 0,
5340 "Return @var{x} raised to the power of @var{y}. This\n"
5341 "procedure does not accept complex arguments.")
5342 #define FUNC_NAME s_scm_sys_expt
5345 scm_two_doubles (x, y, FUNC_NAME, &xy);
5346 return scm_from_double (pow (xy.x, xy.y));
5351 SCM_DEFINE (scm_sys_atan2, "$atan2", 2, 0, 0,
5353 "Return the arc tangent of the two arguments @var{x} and\n"
5354 "@var{y}. This is similar to calculating the arc tangent of\n"
5355 "@var{x} / @var{y}, except that the signs of both arguments\n"
5356 "are used to determine the quadrant of the result. This\n"
5357 "procedure does not accept complex arguments.")
5358 #define FUNC_NAME s_scm_sys_atan2
5361 scm_two_doubles (x, y, FUNC_NAME, &xy);
5362 return scm_from_double (atan2 (xy.x, xy.y));
5367 scm_c_make_rectangular (double re, double im)
5370 return scm_from_double (re);
5374 SCM_NEWSMOB (z, scm_tc16_complex, scm_gc_malloc (sizeof (scm_t_complex),
5376 SCM_COMPLEX_REAL (z) = re;
5377 SCM_COMPLEX_IMAG (z) = im;
5382 SCM_DEFINE (scm_make_rectangular, "make-rectangular", 2, 0, 0,
5383 (SCM real_part, SCM imaginary_part),
5384 "Return a complex number constructed of the given @var{real-part} "
5385 "and @var{imaginary-part} parts.")
5386 #define FUNC_NAME s_scm_make_rectangular
5389 scm_two_doubles (real_part, imaginary_part, FUNC_NAME, &xy);
5390 return scm_c_make_rectangular (xy.x, xy.y);
5395 scm_c_make_polar (double mag, double ang)
5399 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5400 use it on Glibc-based systems that have it (it's a GNU extension). See
5401 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5403 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5404 sincos (ang, &s, &c);
5409 return scm_c_make_rectangular (mag * c, mag * s);
5412 SCM_DEFINE (scm_make_polar, "make-polar", 2, 0, 0,
5414 "Return the complex number @var{x} * e^(i * @var{y}).")
5415 #define FUNC_NAME s_scm_make_polar
5418 scm_two_doubles (x, y, FUNC_NAME, &xy);
5419 return scm_c_make_polar (xy.x, xy.y);
5424 SCM_GPROC (s_real_part, "real-part", 1, 0, 0, scm_real_part, g_real_part);
5425 /* "Return the real part of the number @var{z}."
5428 scm_real_part (SCM z)
5430 if (SCM_I_INUMP (z))
5432 else if (SCM_BIGP (z))
5434 else if (SCM_REALP (z))
5436 else if (SCM_COMPLEXP (z))
5437 return scm_from_double (SCM_COMPLEX_REAL (z));
5438 else if (SCM_FRACTIONP (z))
5441 SCM_WTA_DISPATCH_1 (g_real_part, z, SCM_ARG1, s_real_part);
5445 SCM_GPROC (s_imag_part, "imag-part", 1, 0, 0, scm_imag_part, g_imag_part);
5446 /* "Return the imaginary part of the number @var{z}."
5449 scm_imag_part (SCM z)
5451 if (SCM_I_INUMP (z))
5453 else if (SCM_BIGP (z))
5455 else if (SCM_REALP (z))
5457 else if (SCM_COMPLEXP (z))
5458 return scm_from_double (SCM_COMPLEX_IMAG (z));
5459 else if (SCM_FRACTIONP (z))
5462 SCM_WTA_DISPATCH_1 (g_imag_part, z, SCM_ARG1, s_imag_part);
5465 SCM_GPROC (s_numerator, "numerator", 1, 0, 0, scm_numerator, g_numerator);
5466 /* "Return the numerator of the number @var{z}."
5469 scm_numerator (SCM z)
5471 if (SCM_I_INUMP (z))
5473 else if (SCM_BIGP (z))
5475 else if (SCM_FRACTIONP (z))
5476 return SCM_FRACTION_NUMERATOR (z);
5477 else if (SCM_REALP (z))
5478 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z)));
5480 SCM_WTA_DISPATCH_1 (g_numerator, z, SCM_ARG1, s_numerator);
5484 SCM_GPROC (s_denominator, "denominator", 1, 0, 0, scm_denominator, g_denominator);
5485 /* "Return the denominator of the number @var{z}."
5488 scm_denominator (SCM z)
5490 if (SCM_I_INUMP (z))
5491 return SCM_I_MAKINUM (1);
5492 else if (SCM_BIGP (z))
5493 return SCM_I_MAKINUM (1);
5494 else if (SCM_FRACTIONP (z))
5495 return SCM_FRACTION_DENOMINATOR (z);
5496 else if (SCM_REALP (z))
5497 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z)));
5499 SCM_WTA_DISPATCH_1 (g_denominator, z, SCM_ARG1, s_denominator);
5502 SCM_GPROC (s_magnitude, "magnitude", 1, 0, 0, scm_magnitude, g_magnitude);
5503 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5504 * "@code{abs} for real arguments, but also allows complex numbers."
5507 scm_magnitude (SCM z)
5509 if (SCM_I_INUMP (z))
5511 long int zz = SCM_I_INUM (z);
5514 else if (SCM_POSFIXABLE (-zz))
5515 return SCM_I_MAKINUM (-zz);
5517 return scm_i_long2big (-zz);
5519 else if (SCM_BIGP (z))
5521 int sgn = mpz_sgn (SCM_I_BIG_MPZ (z));
5522 scm_remember_upto_here_1 (z);
5524 return scm_i_clonebig (z, 0);
5528 else if (SCM_REALP (z))
5529 return scm_from_double (fabs (SCM_REAL_VALUE (z)));
5530 else if (SCM_COMPLEXP (z))
5531 return scm_from_double (hypot (SCM_COMPLEX_REAL (z), SCM_COMPLEX_IMAG (z)));
5532 else if (SCM_FRACTIONP (z))
5534 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z))))
5536 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z), SCM_UNDEFINED),
5537 SCM_FRACTION_DENOMINATOR (z));
5540 SCM_WTA_DISPATCH_1 (g_magnitude, z, SCM_ARG1, s_magnitude);
5544 SCM_GPROC (s_angle, "angle", 1, 0, 0, scm_angle, g_angle);
5545 /* "Return the angle of the complex number @var{z}."
5550 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5551 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5552 But if atan2 follows the floating point rounding mode, then the value
5553 is not a constant. Maybe it'd be close enough though. */
5554 if (SCM_I_INUMP (z))
5556 if (SCM_I_INUM (z) >= 0)
5559 return scm_from_double (atan2 (0.0, -1.0));
5561 else if (SCM_BIGP (z))
5563 int sgn = mpz_sgn (SCM_I_BIG_MPZ (z));
5564 scm_remember_upto_here_1 (z);
5566 return scm_from_double (atan2 (0.0, -1.0));
5570 else if (SCM_REALP (z))
5572 if (SCM_REAL_VALUE (z) >= 0)
5575 return scm_from_double (atan2 (0.0, -1.0));
5577 else if (SCM_COMPLEXP (z))
5578 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z), SCM_COMPLEX_REAL (z)));
5579 else if (SCM_FRACTIONP (z))
5581 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z))))
5583 else return scm_from_double (atan2 (0.0, -1.0));
5586 SCM_WTA_DISPATCH_1 (g_angle, z, SCM_ARG1, s_angle);
5590 SCM_GPROC (s_exact_to_inexact, "exact->inexact", 1, 0, 0, scm_exact_to_inexact, g_exact_to_inexact);
5591 /* Convert the number @var{x} to its inexact representation.\n"
5594 scm_exact_to_inexact (SCM z)
5596 if (SCM_I_INUMP (z))
5597 return scm_from_double ((double) SCM_I_INUM (z));
5598 else if (SCM_BIGP (z))
5599 return scm_from_double (scm_i_big2dbl (z));
5600 else if (SCM_FRACTIONP (z))
5601 return scm_from_double (scm_i_fraction2double (z));
5602 else if (SCM_INEXACTP (z))
5605 SCM_WTA_DISPATCH_1 (g_exact_to_inexact, z, 1, s_exact_to_inexact);
5609 SCM_DEFINE (scm_inexact_to_exact, "inexact->exact", 1, 0, 0,
5611 "Return an exact number that is numerically closest to @var{z}.")
5612 #define FUNC_NAME s_scm_inexact_to_exact
5614 if (SCM_I_INUMP (z))
5616 else if (SCM_BIGP (z))
5618 else if (SCM_REALP (z))
5620 if (xisinf (SCM_REAL_VALUE (z)) || xisnan (SCM_REAL_VALUE (z)))
5621 SCM_OUT_OF_RANGE (1, z);
5628 mpq_set_d (frac, SCM_REAL_VALUE (z));
5629 q = scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac)),
5630 scm_i_mpz2num (mpq_denref (frac)));
5632 /* When scm_i_make_ratio throws, we leak the memory allocated
5639 else if (SCM_FRACTIONP (z))
5642 SCM_WRONG_TYPE_ARG (1, z);
5646 SCM_DEFINE (scm_rationalize, "rationalize", 2, 0, 0,
5648 "Returns the @emph{simplest} rational number differing\n"
5649 "from @var{x} by no more than @var{eps}.\n"
5651 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
5652 "exact result when both its arguments are exact. Thus, you might need\n"
5653 "to use @code{inexact->exact} on the arguments.\n"
5656 "(rationalize (inexact->exact 1.2) 1/100)\n"
5659 #define FUNC_NAME s_scm_rationalize
5661 if (SCM_I_INUMP (x))
5663 else if (SCM_BIGP (x))
5665 else if ((SCM_REALP (x)) || SCM_FRACTIONP (x))
5667 /* Use continued fractions to find closest ratio. All
5668 arithmetic is done with exact numbers.
5671 SCM ex = scm_inexact_to_exact (x);
5672 SCM int_part = scm_floor (ex);
5673 SCM tt = SCM_I_MAKINUM (1);
5674 SCM a1 = SCM_I_MAKINUM (0), a2 = SCM_I_MAKINUM (1), a = SCM_I_MAKINUM (0);
5675 SCM b1 = SCM_I_MAKINUM (1), b2 = SCM_I_MAKINUM (0), b = SCM_I_MAKINUM (0);
5679 if (scm_is_true (scm_num_eq_p (ex, int_part)))
5682 ex = scm_difference (ex, int_part); /* x = x-int_part */
5683 rx = scm_divide (ex, SCM_UNDEFINED); /* rx = 1/x */
5685 /* We stop after a million iterations just to be absolutely sure
5686 that we don't go into an infinite loop. The process normally
5687 converges after less than a dozen iterations.
5690 eps = scm_abs (eps);
5691 while (++i < 1000000)
5693 a = scm_sum (scm_product (a1, tt), a2); /* a = a1*tt + a2 */
5694 b = scm_sum (scm_product (b1, tt), b2); /* b = b1*tt + b2 */
5695 if (scm_is_false (scm_zero_p (b)) && /* b != 0 */
5697 (scm_gr_p (scm_abs (scm_difference (ex, scm_divide (a, b))),
5698 eps))) /* abs(x-a/b) <= eps */
5700 SCM res = scm_sum (int_part, scm_divide (a, b));
5701 if (scm_is_false (scm_exact_p (x))
5702 || scm_is_false (scm_exact_p (eps)))
5703 return scm_exact_to_inexact (res);
5707 rx = scm_divide (scm_difference (rx, tt), /* rx = 1/(rx - tt) */
5709 tt = scm_floor (rx); /* tt = floor (rx) */
5715 scm_num_overflow (s_scm_rationalize);
5718 SCM_WRONG_TYPE_ARG (1, x);
5722 /* conversion functions */
5725 scm_is_integer (SCM val)
5727 return scm_is_true (scm_integer_p (val));
5731 scm_is_signed_integer (SCM val, scm_t_intmax min, scm_t_intmax max)
5733 if (SCM_I_INUMP (val))
5735 scm_t_signed_bits n = SCM_I_INUM (val);
5736 return n >= min && n <= max;
5738 else if (SCM_BIGP (val))
5740 if (min >= SCM_MOST_NEGATIVE_FIXNUM && max <= SCM_MOST_POSITIVE_FIXNUM)
5742 else if (min >= LONG_MIN && max <= LONG_MAX)
5744 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val)))
5746 long n = mpz_get_si (SCM_I_BIG_MPZ (val));
5747 return n >= min && n <= max;
5757 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val), 2)
5758 > CHAR_BIT*sizeof (scm_t_uintmax))
5761 mpz_export (&n, &count, 1, sizeof (scm_t_uintmax), 0, 0,
5762 SCM_I_BIG_MPZ (val));
5764 if (mpz_sgn (SCM_I_BIG_MPZ (val)) >= 0)
5776 return n >= min && n <= max;
5784 scm_is_unsigned_integer (SCM val, scm_t_uintmax min, scm_t_uintmax max)
5786 if (SCM_I_INUMP (val))
5788 scm_t_signed_bits n = SCM_I_INUM (val);
5789 return n >= 0 && ((scm_t_uintmax)n) >= min && ((scm_t_uintmax)n) <= max;
5791 else if (SCM_BIGP (val))
5793 if (max <= SCM_MOST_POSITIVE_FIXNUM)
5795 else if (max <= ULONG_MAX)
5797 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val)))
5799 unsigned long n = mpz_get_ui (SCM_I_BIG_MPZ (val));
5800 return n >= min && n <= max;
5810 if (mpz_sgn (SCM_I_BIG_MPZ (val)) < 0)
5813 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val), 2)
5814 > CHAR_BIT*sizeof (scm_t_uintmax))
5817 mpz_export (&n, &count, 1, sizeof (scm_t_uintmax), 0, 0,
5818 SCM_I_BIG_MPZ (val));
5820 return n >= min && n <= max;
5828 scm_i_range_error (SCM bad_val, SCM min, SCM max)
5830 scm_error (scm_out_of_range_key,
5832 "Value out of range ~S to ~S: ~S",
5833 scm_list_3 (min, max, bad_val),
5834 scm_list_1 (bad_val));
5837 #define TYPE scm_t_intmax
5838 #define TYPE_MIN min
5839 #define TYPE_MAX max
5840 #define SIZEOF_TYPE 0
5841 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5842 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5843 #include "libguile/conv-integer.i.c"
5845 #define TYPE scm_t_uintmax
5846 #define TYPE_MIN min
5847 #define TYPE_MAX max
5848 #define SIZEOF_TYPE 0
5849 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5850 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5851 #include "libguile/conv-uinteger.i.c"
5853 #define TYPE scm_t_int8
5854 #define TYPE_MIN SCM_T_INT8_MIN
5855 #define TYPE_MAX SCM_T_INT8_MAX
5856 #define SIZEOF_TYPE 1
5857 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5858 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5859 #include "libguile/conv-integer.i.c"
5861 #define TYPE scm_t_uint8
5863 #define TYPE_MAX SCM_T_UINT8_MAX
5864 #define SIZEOF_TYPE 1
5865 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5866 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5867 #include "libguile/conv-uinteger.i.c"
5869 #define TYPE scm_t_int16
5870 #define TYPE_MIN SCM_T_INT16_MIN
5871 #define TYPE_MAX SCM_T_INT16_MAX
5872 #define SIZEOF_TYPE 2
5873 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5874 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5875 #include "libguile/conv-integer.i.c"
5877 #define TYPE scm_t_uint16
5879 #define TYPE_MAX SCM_T_UINT16_MAX
5880 #define SIZEOF_TYPE 2
5881 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5882 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5883 #include "libguile/conv-uinteger.i.c"
5885 #define TYPE scm_t_int32
5886 #define TYPE_MIN SCM_T_INT32_MIN
5887 #define TYPE_MAX SCM_T_INT32_MAX
5888 #define SIZEOF_TYPE 4
5889 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5890 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5891 #include "libguile/conv-integer.i.c"
5893 #define TYPE scm_t_uint32
5895 #define TYPE_MAX SCM_T_UINT32_MAX
5896 #define SIZEOF_TYPE 4
5897 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5898 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5899 #include "libguile/conv-uinteger.i.c"
5901 #define TYPE scm_t_int64
5902 #define TYPE_MIN SCM_T_INT64_MIN
5903 #define TYPE_MAX SCM_T_INT64_MAX
5904 #define SIZEOF_TYPE 8
5905 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5906 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5907 #include "libguile/conv-integer.i.c"
5909 #define TYPE scm_t_uint64
5911 #define TYPE_MAX SCM_T_UINT64_MAX
5912 #define SIZEOF_TYPE 8
5913 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5914 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5915 #include "libguile/conv-uinteger.i.c"
5918 scm_to_mpz (SCM val, mpz_t rop)
5920 if (SCM_I_INUMP (val))
5921 mpz_set_si (rop, SCM_I_INUM (val));
5922 else if (SCM_BIGP (val))
5923 mpz_set (rop, SCM_I_BIG_MPZ (val));
5925 scm_wrong_type_arg_msg (NULL, 0, val, "exact integer");
5929 scm_from_mpz (mpz_t val)
5931 return scm_i_mpz2num (val);
5935 scm_is_real (SCM val)
5937 return scm_is_true (scm_real_p (val));
5941 scm_is_rational (SCM val)
5943 return scm_is_true (scm_rational_p (val));
5947 scm_to_double (SCM val)
5949 if (SCM_I_INUMP (val))
5950 return SCM_I_INUM (val);
5951 else if (SCM_BIGP (val))
5952 return scm_i_big2dbl (val);
5953 else if (SCM_FRACTIONP (val))
5954 return scm_i_fraction2double (val);
5955 else if (SCM_REALP (val))
5956 return SCM_REAL_VALUE (val);
5958 scm_wrong_type_arg_msg (NULL, 0, val, "real number");
5962 scm_from_double (double val)
5964 SCM z = scm_double_cell (scm_tc16_real, 0, 0, 0);
5965 SCM_REAL_VALUE (z) = val;
5969 #if SCM_ENABLE_DISCOURAGED == 1
5972 scm_num2float (SCM num, unsigned long int pos, const char *s_caller)
5976 float res = mpz_get_d (SCM_I_BIG_MPZ (num));
5980 scm_out_of_range (NULL, num);
5983 return scm_to_double (num);
5987 scm_num2double (SCM num, unsigned long int pos, const char *s_caller)
5991 double res = mpz_get_d (SCM_I_BIG_MPZ (num));
5995 scm_out_of_range (NULL, num);
5998 return scm_to_double (num);
6004 scm_is_complex (SCM val)
6006 return scm_is_true (scm_complex_p (val));
6010 scm_c_real_part (SCM z)
6012 if (SCM_COMPLEXP (z))
6013 return SCM_COMPLEX_REAL (z);
6016 /* Use the scm_real_part to get proper error checking and
6019 return scm_to_double (scm_real_part (z));
6024 scm_c_imag_part (SCM z)
6026 if (SCM_COMPLEXP (z))
6027 return SCM_COMPLEX_IMAG (z);
6030 /* Use the scm_imag_part to get proper error checking and
6031 dispatching. The result will almost always be 0.0, but not
6034 return scm_to_double (scm_imag_part (z));
6039 scm_c_magnitude (SCM z)
6041 return scm_to_double (scm_magnitude (z));
6047 return scm_to_double (scm_angle (z));
6051 scm_is_number (SCM z)
6053 return scm_is_true (scm_number_p (z));
6057 /* In the following functions we dispatch to the real-arg funcs like log()
6058 when we know the arg is real, instead of just handing everything to
6059 clog() for instance. This is in case clog() doesn't optimize for a
6060 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6061 well use it to go straight to the applicable C func. */
6063 SCM_DEFINE (scm_log, "log", 1, 0, 0,
6065 "Return the natural logarithm of @var{z}.")
6066 #define FUNC_NAME s_scm_log
6068 if (SCM_COMPLEXP (z))
6070 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6071 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z)));
6073 double re = SCM_COMPLEX_REAL (z);
6074 double im = SCM_COMPLEX_IMAG (z);
6075 return scm_c_make_rectangular (log (hypot (re, im)),
6081 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6082 although the value itself overflows. */
6083 double re = scm_to_double (z);
6084 double l = log (fabs (re));
6086 return scm_from_double (l);
6088 return scm_c_make_rectangular (l, M_PI);
6094 SCM_DEFINE (scm_log10, "log10", 1, 0, 0,
6096 "Return the base 10 logarithm of @var{z}.")
6097 #define FUNC_NAME s_scm_log10
6099 if (SCM_COMPLEXP (z))
6101 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6102 clog() and a multiply by M_LOG10E, rather than the fallback
6103 log10+hypot+atan2.) */
6104 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6105 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z)));
6107 double re = SCM_COMPLEX_REAL (z);
6108 double im = SCM_COMPLEX_IMAG (z);
6109 return scm_c_make_rectangular (log10 (hypot (re, im)),
6110 M_LOG10E * atan2 (im, re));
6115 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6116 although the value itself overflows. */
6117 double re = scm_to_double (z);
6118 double l = log10 (fabs (re));
6120 return scm_from_double (l);
6122 return scm_c_make_rectangular (l, M_LOG10E * M_PI);
6128 SCM_DEFINE (scm_exp, "exp", 1, 0, 0,
6130 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6131 "base of natural logarithms (2.71828@dots{}).")
6132 #define FUNC_NAME s_scm_exp
6134 if (SCM_COMPLEXP (z))
6136 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6137 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z)));
6139 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z)),
6140 SCM_COMPLEX_IMAG (z));
6145 /* When z is a negative bignum the conversion to double overflows,
6146 giving -infinity, but that's ok, the exp is still 0.0. */
6147 return scm_from_double (exp (scm_to_double (z)));
6153 SCM_DEFINE (scm_sqrt, "sqrt", 1, 0, 0,
6155 "Return the square root of @var{z}. Of the two possible roots\n"
6156 "(positive and negative), the one with the a positive real part\n"
6157 "is returned, or if that's zero then a positive imaginary part.\n"
6161 "(sqrt 9.0) @result{} 3.0\n"
6162 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6163 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6164 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6166 #define FUNC_NAME s_scm_sqrt
6168 if (SCM_COMPLEXP (x))
6170 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6171 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x)));
6173 double re = SCM_COMPLEX_REAL (x);
6174 double im = SCM_COMPLEX_IMAG (x);
6175 return scm_c_make_polar (sqrt (hypot (re, im)),
6176 0.5 * atan2 (im, re));
6181 double xx = scm_to_double (x);
6183 return scm_c_make_rectangular (0.0, sqrt (-xx));
6185 return scm_from_double (sqrt (xx));
6197 mpz_init_set_si (z_negative_one, -1);
6199 /* It may be possible to tune the performance of some algorithms by using
6200 * the following constants to avoid the creation of bignums. Please, before
6201 * using these values, remember the two rules of program optimization:
6202 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6203 scm_c_define ("most-positive-fixnum",
6204 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM));
6205 scm_c_define ("most-negative-fixnum",
6206 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM));
6208 scm_add_feature ("complex");
6209 scm_add_feature ("inexact");
6210 scm_flo0 = scm_from_double (0.0);
6212 /* determine floating point precision */
6213 for (i=2; i <= SCM_MAX_DBL_RADIX; ++i)
6215 init_dblprec(&scm_dblprec[i-2],i);
6216 init_fx_radix(fx_per_radix[i-2],i);
6219 /* hard code precision for base 10 if the preprocessor tells us to... */
6220 scm_dblprec[10-2] = (DBL_DIG > 20) ? 20 : DBL_DIG;
6223 exactly_one_half = scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6224 SCM_I_MAKINUM (2)));
6225 #include "libguile/numbers.x"