1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 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.
43 /* tell glibc (2.3) to give prototype for C99 trunc() */
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"
71 Wonder if this might be faster for some of our code? A switch on
72 the numtag would jump directly to the right case, and the
73 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
75 #define SCM_I_NUMTAG_NOTNUM 0
76 #define SCM_I_NUMTAG_INUM 1
77 #define SCM_I_NUMTAG_BIG scm_tc16_big
78 #define SCM_I_NUMTAG_REAL scm_tc16_real
79 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
80 #define SCM_I_NUMTAG(x) \
81 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
82 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
83 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
84 : SCM_I_NUMTAG_NOTNUM)))
86 /* the macro above will not work as is with fractions */
89 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
91 /* FLOBUFLEN is the maximum number of characters neccessary for the
92 * printed or scm_string representation of an inexact number.
94 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
97 #if ! defined (HAVE_ISNAN)
102 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
105 #if ! defined (HAVE_ISINF)
110 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
117 /* mpz_cmp_d only recognises infinities in gmp 4.2 and up.
118 For prior versions use an explicit check here. */
119 #if __GNU_MP_VERSION < 4 \
120 || (__GNU_MP_VERSION == 4 && __GNU_MP_VERSION_MINOR < 2)
121 #define xmpz_cmp_d(z, d) \
122 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
124 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
127 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
128 isinf. It does have finite and isnan though, hence the use of those.
129 fpclass would be a possibility on that system too. */
133 #if defined (HAVE_ISINF)
135 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
136 return (! (finite (x
) || isnan (x
)));
145 #if defined (HAVE_ISNAN)
154 static mpz_t z_negative_one
;
158 SCM_C_INLINE_KEYWORD SCM
161 /* Return a newly created bignum. */
162 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
163 mpz_init (SCM_I_BIG_MPZ (z
));
167 SCM_C_INLINE_KEYWORD
static SCM
168 scm_i_clonebig (SCM src_big
, int same_sign_p
)
170 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
171 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
172 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
174 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
178 SCM_C_INLINE_KEYWORD
int
179 scm_i_bigcmp (SCM x
, SCM y
)
181 /* Return neg if x < y, pos if x > y, and 0 if x == y */
182 /* presume we already know x and y are bignums */
183 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
184 scm_remember_upto_here_2 (x
, y
);
188 SCM_C_INLINE_KEYWORD SCM
189 scm_i_dbl2big (double d
)
191 /* results are only defined if d is an integer */
192 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
193 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
197 /* Convert a integer in double representation to a SCM number. */
199 SCM_C_INLINE_KEYWORD SCM
200 scm_i_dbl2num (double u
)
202 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
203 powers of 2, so there's no rounding when making "double" values
204 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
205 get rounded on a 64-bit machine, hence the "+1".
207 The use of floor() to force to an integer value ensures we get a
208 "numerically closest" value without depending on how a
209 double->long cast or how mpz_set_d will round. For reference,
210 double->long probably follows the hardware rounding mode,
211 mpz_set_d truncates towards zero. */
213 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
214 representable as a double? */
216 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
217 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
218 return SCM_I_MAKINUM ((long) u
);
220 return scm_i_dbl2big (u
);
223 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
224 with R5RS exact->inexact.
226 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
227 (ie. it truncates towards zero), then adjust to get the closest double by
228 examining the next lower bit and adding 1 if necessary.
230 Note that bignums exactly half way between representable doubles are
231 rounded to the next higher absolute value (ie. away from zero). This
232 seems like an adequate interpretation of R5RS "numerically closest", and
233 it's easier and faster than a full "nearest-even" style.
235 The bit test is done on the absolute value of the mpz_t, which means we
236 must use mpz_getlimbn. mpz_tstbit is not right, it treats negatives as
239 Prior to GMP 4.2, the rounding done by mpz_get_d was unspecified. It
240 happened to follow the hardware rounding mode, but on the absolute value
241 of its operand. This is not what we want, so we put the high
242 DBL_MANT_DIG bits into a temporary. This extra init/clear is a slowdown,
243 but doesn't matter too much since it's only for older GMP. */
246 scm_i_big2dbl (SCM b
)
251 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
253 #if __GNU_MP_VERSION < 4 \
254 || (__GNU_MP_VERSION == 4 && __GNU_MP_VERSION_MINOR < 2)
256 /* GMP prior to 4.2, force truncate towards zero */
258 if (bits
> DBL_MANT_DIG
)
260 size_t shift
= bits
- DBL_MANT_DIG
;
261 mpz_init2 (tmp
, DBL_MANT_DIG
);
262 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
263 result
= ldexp (mpz_get_d (tmp
), shift
);
268 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
273 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
276 if (bits
> DBL_MANT_DIG
)
278 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
279 /* test bit number "pos" in absolute value */
280 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
281 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
283 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
287 scm_remember_upto_here_1 (b
);
291 SCM_C_INLINE_KEYWORD SCM
292 scm_i_normbig (SCM b
)
294 /* convert a big back to a fixnum if it'll fit */
295 /* presume b is a bignum */
296 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
298 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
299 if (SCM_FIXABLE (val
))
300 b
= SCM_I_MAKINUM (val
);
305 static SCM_C_INLINE_KEYWORD SCM
306 scm_i_mpz2num (mpz_t b
)
308 /* convert a mpz number to a SCM number. */
309 if (mpz_fits_slong_p (b
))
311 long val
= mpz_get_si (b
);
312 if (SCM_FIXABLE (val
))
313 return SCM_I_MAKINUM (val
);
317 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
318 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
323 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
324 static SCM
scm_divide2real (SCM x
, SCM y
);
327 scm_make_ratio (SCM numerator
, SCM denominator
)
328 #define FUNC_NAME "make-ratio"
330 /* First make sure the arguments are proper.
332 if (SCM_I_INUMP (denominator
))
334 if (scm_is_eq (denominator
, SCM_INUM0
))
335 scm_num_overflow ("make-ratio");
336 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
341 if (!(SCM_BIGP(denominator
)))
342 SCM_WRONG_TYPE_ARG (2, denominator
);
344 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
345 SCM_WRONG_TYPE_ARG (1, numerator
);
347 /* Then flip signs so that the denominator is positive.
349 if (scm_is_true (scm_negative_p (denominator
)))
351 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
352 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
355 /* Now consider for each of the four fixnum/bignum combinations
356 whether the rational number is really an integer.
358 if (SCM_I_INUMP (numerator
))
360 long x
= SCM_I_INUM (numerator
);
361 if (scm_is_eq (numerator
, SCM_INUM0
))
363 if (SCM_I_INUMP (denominator
))
366 y
= SCM_I_INUM (denominator
);
368 return SCM_I_MAKINUM(1);
370 return SCM_I_MAKINUM (x
/ y
);
374 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
375 of that value for the denominator, as a bignum. Apart from
376 that case, abs(bignum) > abs(inum) so inum/bignum is not an
378 if (x
== SCM_MOST_NEGATIVE_FIXNUM
379 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
380 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
381 return SCM_I_MAKINUM(-1);
384 else if (SCM_BIGP (numerator
))
386 if (SCM_I_INUMP (denominator
))
388 long yy
= SCM_I_INUM (denominator
);
389 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
390 return scm_divide (numerator
, denominator
);
394 if (scm_is_eq (numerator
, denominator
))
395 return SCM_I_MAKINUM(1);
396 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
397 SCM_I_BIG_MPZ (denominator
)))
398 return scm_divide(numerator
, denominator
);
402 /* No, it's a proper fraction.
404 return scm_double_cell (scm_tc16_fraction
,
405 SCM_UNPACK (numerator
),
406 SCM_UNPACK (denominator
), 0);
410 static void scm_i_fraction_reduce (SCM z
)
412 if (!(SCM_FRACTION_REDUCED (z
)))
415 divisor
= scm_gcd (SCM_FRACTION_NUMERATOR (z
), SCM_FRACTION_DENOMINATOR (z
));
416 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
419 SCM_FRACTION_SET_NUMERATOR (z
, scm_divide (SCM_FRACTION_NUMERATOR (z
), divisor
));
420 SCM_FRACTION_SET_DENOMINATOR (z
, scm_divide (SCM_FRACTION_DENOMINATOR (z
), divisor
));
422 SCM_FRACTION_REDUCED_SET (z
);
427 scm_i_fraction2double (SCM z
)
429 return scm_num2dbl (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
430 SCM_FRACTION_DENOMINATOR (z
)),
434 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
436 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
438 #define FUNC_NAME s_scm_exact_p
444 if (SCM_FRACTIONP (x
))
448 SCM_WRONG_TYPE_ARG (1, x
);
453 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
455 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
457 #define FUNC_NAME s_scm_odd_p
461 long val
= SCM_I_INUM (n
);
462 return scm_from_bool ((val
& 1L) != 0);
464 else if (SCM_BIGP (n
))
466 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
467 scm_remember_upto_here_1 (n
);
468 return scm_from_bool (odd_p
);
470 else if (scm_is_true (scm_inf_p (n
)))
472 else if (SCM_REALP (n
))
474 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
480 SCM_WRONG_TYPE_ARG (1, n
);
483 SCM_WRONG_TYPE_ARG (1, n
);
488 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
490 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
492 #define FUNC_NAME s_scm_even_p
496 long val
= SCM_I_INUM (n
);
497 return scm_from_bool ((val
& 1L) == 0);
499 else if (SCM_BIGP (n
))
501 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
502 scm_remember_upto_here_1 (n
);
503 return scm_from_bool (even_p
);
505 else if (scm_is_true (scm_inf_p (n
)))
507 else if (SCM_REALP (n
))
509 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
515 SCM_WRONG_TYPE_ARG (1, n
);
518 SCM_WRONG_TYPE_ARG (1, n
);
522 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
524 "Return @code{#t} if @var{n} is infinite, @code{#f}\n"
526 #define FUNC_NAME s_scm_inf_p
529 return scm_from_bool (xisinf (SCM_REAL_VALUE (n
)));
530 else if (SCM_COMPLEXP (n
))
531 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (n
))
532 || xisinf (SCM_COMPLEX_IMAG (n
)));
538 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
540 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
542 #define FUNC_NAME s_scm_nan_p
545 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
546 else if (SCM_COMPLEXP (n
))
547 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
548 || xisnan (SCM_COMPLEX_IMAG (n
)));
554 /* Guile's idea of infinity. */
555 static double guile_Inf
;
557 /* Guile's idea of not a number. */
558 static double guile_NaN
;
561 guile_ieee_init (void)
563 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
565 /* Some version of gcc on some old version of Linux used to crash when
566 trying to make Inf and NaN. */
569 /* C99 INFINITY, when available.
570 FIXME: The standard allows for INFINITY to be something that overflows
571 at compile time. We ought to have a configure test to check for that
572 before trying to use it. (But in practice we believe this is not a
573 problem on any system guile is likely to target.) */
574 guile_Inf
= INFINITY
;
577 extern unsigned int DINFINITY
[2];
578 guile_Inf
= (*(X_CAST(double *, DINFINITY
)));
585 if (guile_Inf
== tmp
)
593 #if defined (HAVE_ISNAN)
596 /* C99 NAN, when available */
600 extern unsigned int DQNAN
[2];
601 guile_NaN
= (*(X_CAST(double *, DQNAN
)));
603 guile_NaN
= guile_Inf
/ guile_Inf
;
609 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
612 #define FUNC_NAME s_scm_inf
614 static int initialized
= 0;
620 return scm_make_real (guile_Inf
);
624 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
627 #define FUNC_NAME s_scm_nan
629 static int initialized
= 0;
635 return scm_make_real (guile_NaN
);
640 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
642 "Return the absolute value of @var{x}.")
647 long int xx
= SCM_I_INUM (x
);
650 else if (SCM_POSFIXABLE (-xx
))
651 return SCM_I_MAKINUM (-xx
);
653 return scm_i_long2big (-xx
);
655 else if (SCM_BIGP (x
))
657 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
659 return scm_i_clonebig (x
, 0);
663 else if (SCM_REALP (x
))
665 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
666 double xx
= SCM_REAL_VALUE (x
);
668 return scm_make_real (-xx
);
672 else if (SCM_FRACTIONP (x
))
674 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
676 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
677 SCM_FRACTION_DENOMINATOR (x
));
680 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
685 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
686 /* "Return the quotient of the numbers @var{x} and @var{y}."
689 scm_quotient (SCM x
, SCM y
)
693 long xx
= SCM_I_INUM (x
);
696 long yy
= SCM_I_INUM (y
);
698 scm_num_overflow (s_quotient
);
703 return SCM_I_MAKINUM (z
);
705 return scm_i_long2big (z
);
708 else if (SCM_BIGP (y
))
710 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
711 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
712 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
714 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
715 scm_remember_upto_here_1 (y
);
716 return SCM_I_MAKINUM (-1);
719 return SCM_I_MAKINUM (0);
722 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
724 else if (SCM_BIGP (x
))
728 long yy
= SCM_I_INUM (y
);
730 scm_num_overflow (s_quotient
);
735 SCM result
= scm_i_mkbig ();
738 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
741 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
744 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
745 scm_remember_upto_here_1 (x
);
746 return scm_i_normbig (result
);
749 else if (SCM_BIGP (y
))
751 SCM result
= scm_i_mkbig ();
752 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
755 scm_remember_upto_here_2 (x
, y
);
756 return scm_i_normbig (result
);
759 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
762 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
765 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
766 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
768 * "(remainder 13 4) @result{} 1\n"
769 * "(remainder -13 4) @result{} -1\n"
773 scm_remainder (SCM x
, SCM y
)
779 long yy
= SCM_I_INUM (y
);
781 scm_num_overflow (s_remainder
);
784 long z
= SCM_I_INUM (x
) % yy
;
785 return SCM_I_MAKINUM (z
);
788 else if (SCM_BIGP (y
))
790 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
791 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
792 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
794 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
795 scm_remember_upto_here_1 (y
);
796 return SCM_I_MAKINUM (0);
802 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
804 else if (SCM_BIGP (x
))
808 long yy
= SCM_I_INUM (y
);
810 scm_num_overflow (s_remainder
);
813 SCM result
= scm_i_mkbig ();
816 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
817 scm_remember_upto_here_1 (x
);
818 return scm_i_normbig (result
);
821 else if (SCM_BIGP (y
))
823 SCM result
= scm_i_mkbig ();
824 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
827 scm_remember_upto_here_2 (x
, y
);
828 return scm_i_normbig (result
);
831 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
834 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
838 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
839 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
841 * "(modulo 13 4) @result{} 1\n"
842 * "(modulo -13 4) @result{} 3\n"
846 scm_modulo (SCM x
, SCM y
)
850 long xx
= SCM_I_INUM (x
);
853 long yy
= SCM_I_INUM (y
);
855 scm_num_overflow (s_modulo
);
858 /* FIXME: I think this may be a bug on some arches -- results
859 of % with negative second arg are undefined... */
877 return SCM_I_MAKINUM (result
);
880 else if (SCM_BIGP (y
))
882 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
889 SCM pos_y
= scm_i_clonebig (y
, 0);
890 /* do this after the last scm_op */
891 mpz_init_set_si (z_x
, xx
);
892 result
= pos_y
; /* re-use this bignum */
893 mpz_mod (SCM_I_BIG_MPZ (result
),
895 SCM_I_BIG_MPZ (pos_y
));
896 scm_remember_upto_here_1 (pos_y
);
900 result
= scm_i_mkbig ();
901 /* do this after the last scm_op */
902 mpz_init_set_si (z_x
, xx
);
903 mpz_mod (SCM_I_BIG_MPZ (result
),
906 scm_remember_upto_here_1 (y
);
909 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
910 mpz_add (SCM_I_BIG_MPZ (result
),
912 SCM_I_BIG_MPZ (result
));
913 scm_remember_upto_here_1 (y
);
914 /* and do this before the next one */
916 return scm_i_normbig (result
);
920 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
922 else if (SCM_BIGP (x
))
926 long yy
= SCM_I_INUM (y
);
928 scm_num_overflow (s_modulo
);
931 SCM result
= scm_i_mkbig ();
932 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
934 (yy
< 0) ? - yy
: yy
);
935 scm_remember_upto_here_1 (x
);
936 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
937 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
938 SCM_I_BIG_MPZ (result
),
940 return scm_i_normbig (result
);
943 else if (SCM_BIGP (y
))
946 SCM result
= scm_i_mkbig ();
947 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
948 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
949 mpz_mod (SCM_I_BIG_MPZ (result
),
951 SCM_I_BIG_MPZ (pos_y
));
953 scm_remember_upto_here_1 (x
);
954 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
955 mpz_add (SCM_I_BIG_MPZ (result
),
957 SCM_I_BIG_MPZ (result
));
958 scm_remember_upto_here_2 (y
, pos_y
);
959 return scm_i_normbig (result
);
963 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
966 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
969 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
970 /* "Return the greatest common divisor of all arguments.\n"
971 * "If called without arguments, 0 is returned."
974 scm_gcd (SCM x
, SCM y
)
977 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
983 long xx
= SCM_I_INUM (x
);
984 long yy
= SCM_I_INUM (y
);
985 long u
= xx
< 0 ? -xx
: xx
;
986 long v
= yy
< 0 ? -yy
: yy
;
996 /* Determine a common factor 2^k */
997 while (!(1 & (u
| v
)))
1003 /* Now, any factor 2^n can be eliminated */
1023 return (SCM_POSFIXABLE (result
)
1024 ? SCM_I_MAKINUM (result
)
1025 : scm_i_long2big (result
));
1027 else if (SCM_BIGP (y
))
1033 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1035 else if (SCM_BIGP (x
))
1037 if (SCM_I_INUMP (y
))
1039 unsigned long result
;
1042 yy
= SCM_I_INUM (y
);
1047 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1048 scm_remember_upto_here_1 (x
);
1049 return (SCM_POSFIXABLE (result
)
1050 ? SCM_I_MAKINUM (result
)
1051 : scm_ulong2num (result
));
1053 else if (SCM_BIGP (y
))
1055 SCM result
= scm_i_mkbig ();
1056 mpz_gcd (SCM_I_BIG_MPZ (result
),
1059 scm_remember_upto_here_2 (x
, y
);
1060 return scm_i_normbig (result
);
1063 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1066 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1069 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1070 /* "Return the least common multiple of the arguments.\n"
1071 * "If called without arguments, 1 is returned."
1074 scm_lcm (SCM n1
, SCM n2
)
1076 if (SCM_UNBNDP (n2
))
1078 if (SCM_UNBNDP (n1
))
1079 return SCM_I_MAKINUM (1L);
1080 n2
= SCM_I_MAKINUM (1L);
1083 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1084 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1085 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1086 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1088 if (SCM_I_INUMP (n1
))
1090 if (SCM_I_INUMP (n2
))
1092 SCM d
= scm_gcd (n1
, n2
);
1093 if (scm_is_eq (d
, SCM_INUM0
))
1096 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1100 /* inum n1, big n2 */
1103 SCM result
= scm_i_mkbig ();
1104 long nn1
= SCM_I_INUM (n1
);
1105 if (nn1
== 0) return SCM_INUM0
;
1106 if (nn1
< 0) nn1
= - nn1
;
1107 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1108 scm_remember_upto_here_1 (n2
);
1116 if (SCM_I_INUMP (n2
))
1123 SCM result
= scm_i_mkbig ();
1124 mpz_lcm(SCM_I_BIG_MPZ (result
),
1126 SCM_I_BIG_MPZ (n2
));
1127 scm_remember_upto_here_2(n1
, n2
);
1128 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1134 #ifndef scm_long2num
1135 #define SCM_LOGOP_RETURN(x) scm_ulong2num(x)
1137 #define SCM_LOGOP_RETURN(x) SCM_I_MAKINUM(x)
1140 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1145 + + + x (map digit:logand X Y)
1146 + - + x (map digit:logand X (lognot (+ -1 Y)))
1147 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1148 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1153 + + + (map digit:logior X Y)
1154 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1155 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1156 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1161 + + + (map digit:logxor X Y)
1162 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1163 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1164 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1169 + + (any digit:logand X Y)
1170 + - (any digit:logand X (lognot (+ -1 Y)))
1171 - + (any digit:logand (lognot (+ -1 X)) Y)
1176 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1178 "Return the bitwise AND of the integer arguments.\n\n"
1180 "(logand) @result{} -1\n"
1181 "(logand 7) @result{} 7\n"
1182 "(logand #b111 #b011 #b001) @result{} 1\n"
1184 #define FUNC_NAME s_scm_logand
1188 if (SCM_UNBNDP (n2
))
1190 if (SCM_UNBNDP (n1
))
1191 return SCM_I_MAKINUM (-1);
1192 else if (!SCM_NUMBERP (n1
))
1193 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1194 else if (SCM_NUMBERP (n1
))
1197 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1200 if (SCM_I_INUMP (n1
))
1202 nn1
= SCM_I_INUM (n1
);
1203 if (SCM_I_INUMP (n2
))
1205 long nn2
= SCM_I_INUM (n2
);
1206 return SCM_I_MAKINUM (nn1
& nn2
);
1208 else if SCM_BIGP (n2
)
1214 SCM result_z
= scm_i_mkbig ();
1216 mpz_init_set_si (nn1_z
, nn1
);
1217 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1218 scm_remember_upto_here_1 (n2
);
1220 return scm_i_normbig (result_z
);
1224 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1226 else if (SCM_BIGP (n1
))
1228 if (SCM_I_INUMP (n2
))
1231 nn1
= SCM_I_INUM (n1
);
1234 else if (SCM_BIGP (n2
))
1236 SCM result_z
= scm_i_mkbig ();
1237 mpz_and (SCM_I_BIG_MPZ (result_z
),
1239 SCM_I_BIG_MPZ (n2
));
1240 scm_remember_upto_here_2 (n1
, n2
);
1241 return scm_i_normbig (result_z
);
1244 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1247 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1252 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1254 "Return the bitwise OR of the integer arguments.\n\n"
1256 "(logior) @result{} 0\n"
1257 "(logior 7) @result{} 7\n"
1258 "(logior #b000 #b001 #b011) @result{} 3\n"
1260 #define FUNC_NAME s_scm_logior
1264 if (SCM_UNBNDP (n2
))
1266 if (SCM_UNBNDP (n1
))
1268 else if (SCM_NUMBERP (n1
))
1271 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1274 if (SCM_I_INUMP (n1
))
1276 nn1
= SCM_I_INUM (n1
);
1277 if (SCM_I_INUMP (n2
))
1279 long nn2
= SCM_I_INUM (n2
);
1280 return SCM_I_MAKINUM (nn1
| nn2
);
1282 else if (SCM_BIGP (n2
))
1288 SCM result_z
= scm_i_mkbig ();
1290 mpz_init_set_si (nn1_z
, nn1
);
1291 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1292 scm_remember_upto_here_1 (n2
);
1298 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1300 else if (SCM_BIGP (n1
))
1302 if (SCM_I_INUMP (n2
))
1305 nn1
= SCM_I_INUM (n1
);
1308 else if (SCM_BIGP (n2
))
1310 SCM result_z
= scm_i_mkbig ();
1311 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1313 SCM_I_BIG_MPZ (n2
));
1314 scm_remember_upto_here_2 (n1
, n2
);
1318 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1321 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1326 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1328 "Return the bitwise XOR of the integer arguments. A bit is\n"
1329 "set in the result if it is set in an odd number of arguments.\n"
1331 "(logxor) @result{} 0\n"
1332 "(logxor 7) @result{} 7\n"
1333 "(logxor #b000 #b001 #b011) @result{} 2\n"
1334 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1336 #define FUNC_NAME s_scm_logxor
1340 if (SCM_UNBNDP (n2
))
1342 if (SCM_UNBNDP (n1
))
1344 else if (SCM_NUMBERP (n1
))
1347 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1350 if (SCM_I_INUMP (n1
))
1352 nn1
= SCM_I_INUM (n1
);
1353 if (SCM_I_INUMP (n2
))
1355 long nn2
= SCM_I_INUM (n2
);
1356 return SCM_I_MAKINUM (nn1
^ nn2
);
1358 else if (SCM_BIGP (n2
))
1362 SCM result_z
= scm_i_mkbig ();
1364 mpz_init_set_si (nn1_z
, nn1
);
1365 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1366 scm_remember_upto_here_1 (n2
);
1368 return scm_i_normbig (result_z
);
1372 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1374 else if (SCM_BIGP (n1
))
1376 if (SCM_I_INUMP (n2
))
1379 nn1
= SCM_I_INUM (n1
);
1382 else if (SCM_BIGP (n2
))
1384 SCM result_z
= scm_i_mkbig ();
1385 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1387 SCM_I_BIG_MPZ (n2
));
1388 scm_remember_upto_here_2 (n1
, n2
);
1389 return scm_i_normbig (result_z
);
1392 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1395 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1400 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1403 "(logtest j k) @equiv{} (not (zero? (logand j k)))\n\n"
1404 "(logtest #b0100 #b1011) @result{} #f\n"
1405 "(logtest #b0100 #b0111) @result{} #t\n"
1407 #define FUNC_NAME s_scm_logtest
1411 if (SCM_I_INUMP (j
))
1413 nj
= SCM_I_INUM (j
);
1414 if (SCM_I_INUMP (k
))
1416 long nk
= SCM_I_INUM (k
);
1417 return scm_from_bool (nj
& nk
);
1419 else if (SCM_BIGP (k
))
1427 mpz_init_set_si (nj_z
, nj
);
1428 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1429 scm_remember_upto_here_1 (k
);
1430 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1436 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1438 else if (SCM_BIGP (j
))
1440 if (SCM_I_INUMP (k
))
1443 nj
= SCM_I_INUM (j
);
1446 else if (SCM_BIGP (k
))
1450 mpz_init (result_z
);
1454 scm_remember_upto_here_2 (j
, k
);
1455 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1456 mpz_clear (result_z
);
1460 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1463 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1468 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1471 "(logbit? index j) @equiv{} (logtest (integer-expt 2 index) j)\n\n"
1472 "(logbit? 0 #b1101) @result{} #t\n"
1473 "(logbit? 1 #b1101) @result{} #f\n"
1474 "(logbit? 2 #b1101) @result{} #t\n"
1475 "(logbit? 3 #b1101) @result{} #t\n"
1476 "(logbit? 4 #b1101) @result{} #f\n"
1478 #define FUNC_NAME s_scm_logbit_p
1480 unsigned long int iindex
;
1481 iindex
= scm_to_ulong (index
);
1483 if (SCM_I_INUMP (j
))
1485 /* bits above what's in an inum follow the sign bit */
1486 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1487 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1489 else if (SCM_BIGP (j
))
1491 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1492 scm_remember_upto_here_1 (j
);
1493 return scm_from_bool (val
);
1496 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1501 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1503 "Return the integer which is the ones-complement of the integer\n"
1507 "(number->string (lognot #b10000000) 2)\n"
1508 " @result{} \"-10000001\"\n"
1509 "(number->string (lognot #b0) 2)\n"
1510 " @result{} \"-1\"\n"
1512 #define FUNC_NAME s_scm_lognot
1514 if (SCM_I_INUMP (n
)) {
1515 /* No overflow here, just need to toggle all the bits making up the inum.
1516 Enhancement: No need to strip the tag and add it back, could just xor
1517 a block of 1 bits, if that worked with the various debug versions of
1519 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1521 } else if (SCM_BIGP (n
)) {
1522 SCM result
= scm_i_mkbig ();
1523 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1524 scm_remember_upto_here_1 (n
);
1528 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1533 /* returns 0 if IN is not an integer. OUT must already be
1536 coerce_to_big (SCM in
, mpz_t out
)
1539 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1540 else if (SCM_I_INUMP (in
))
1541 mpz_set_si (out
, SCM_I_INUM (in
));
1548 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1549 (SCM n
, SCM k
, SCM m
),
1550 "Return @var{n} raised to the integer exponent\n"
1551 "@var{k}, modulo @var{m}.\n"
1554 "(modulo-expt 2 3 5)\n"
1557 #define FUNC_NAME s_scm_modulo_expt
1563 /* There are two classes of error we might encounter --
1564 1) Math errors, which we'll report by calling scm_num_overflow,
1566 2) wrong-type errors, which of course we'll report by calling
1568 We don't report those errors immediately, however; instead we do
1569 some cleanup first. These variables tell us which error (if
1570 any) we should report after cleaning up.
1572 int report_overflow
= 0;
1574 int position_of_wrong_type
= 0;
1575 SCM value_of_wrong_type
= SCM_INUM0
;
1577 SCM result
= SCM_UNDEFINED
;
1583 if (scm_is_eq (m
, SCM_INUM0
))
1585 report_overflow
= 1;
1589 if (!coerce_to_big (n
, n_tmp
))
1591 value_of_wrong_type
= n
;
1592 position_of_wrong_type
= 1;
1596 if (!coerce_to_big (k
, k_tmp
))
1598 value_of_wrong_type
= k
;
1599 position_of_wrong_type
= 2;
1603 if (!coerce_to_big (m
, m_tmp
))
1605 value_of_wrong_type
= m
;
1606 position_of_wrong_type
= 3;
1610 /* if the exponent K is negative, and we simply call mpz_powm, we
1611 will get a divide-by-zero exception when an inverse 1/n mod m
1612 doesn't exist (or is not unique). Since exceptions are hard to
1613 handle, we'll attempt the inversion "by hand" -- that way, we get
1614 a simple failure code, which is easy to handle. */
1616 if (-1 == mpz_sgn (k_tmp
))
1618 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1620 report_overflow
= 1;
1623 mpz_neg (k_tmp
, k_tmp
);
1626 result
= scm_i_mkbig ();
1627 mpz_powm (SCM_I_BIG_MPZ (result
),
1632 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1633 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1640 if (report_overflow
)
1641 scm_num_overflow (FUNC_NAME
);
1643 if (position_of_wrong_type
)
1644 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1645 value_of_wrong_type
);
1647 return scm_i_normbig (result
);
1651 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1653 "Return @var{n} raised to the non-negative integer exponent\n"
1657 "(integer-expt 2 5)\n"
1659 "(integer-expt -3 3)\n"
1662 #define FUNC_NAME s_scm_integer_expt
1665 SCM z_i2
= SCM_BOOL_F
;
1667 SCM acc
= SCM_I_MAKINUM (1L);
1669 /* 0^0 == 1 according to R5RS */
1670 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1671 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1672 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1673 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1675 if (SCM_I_INUMP (k
))
1676 i2
= SCM_I_INUM (k
);
1677 else if (SCM_BIGP (k
))
1679 z_i2
= scm_i_clonebig (k
, 1);
1680 scm_remember_upto_here_1 (k
);
1683 else if (SCM_REALP (k
))
1685 double r
= SCM_REAL_VALUE (k
);
1687 SCM_WRONG_TYPE_ARG (2, k
);
1688 if ((r
> SCM_MOST_POSITIVE_FIXNUM
) || (r
< SCM_MOST_NEGATIVE_FIXNUM
))
1690 z_i2
= scm_i_mkbig ();
1691 mpz_set_d (SCM_I_BIG_MPZ (z_i2
), r
);
1700 SCM_WRONG_TYPE_ARG (2, k
);
1704 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1706 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1707 n
= scm_divide (n
, SCM_UNDEFINED
);
1711 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1715 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1717 return scm_product (acc
, n
);
1719 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1720 acc
= scm_product (acc
, n
);
1721 n
= scm_product (n
, n
);
1722 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1730 n
= scm_divide (n
, SCM_UNDEFINED
);
1737 return scm_product (acc
, n
);
1739 acc
= scm_product (acc
, n
);
1740 n
= scm_product (n
, n
);
1747 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1749 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1750 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1752 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1753 "@var{cnt} is negative it's a division, rounded towards negative\n"
1754 "infinity. (Note that this is not the same rounding as\n"
1755 "@code{quotient} does.)\n"
1757 "With @var{n} viewed as an infinite precision twos complement,\n"
1758 "@code{ash} means a left shift introducing zero bits, or a right\n"
1759 "shift dropping bits.\n"
1762 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1763 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1765 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1766 "(ash -23 -2) @result{} -6\n"
1768 #define FUNC_NAME s_scm_ash
1771 bits_to_shift
= scm_to_long (cnt
);
1773 if (bits_to_shift
< 0)
1775 /* Shift right by abs(cnt) bits. This is realized as a division
1776 by div:=2^abs(cnt). However, to guarantee the floor
1777 rounding, negative values require some special treatment.
1779 SCM div
= scm_integer_expt (SCM_I_MAKINUM (2),
1780 scm_from_long (-bits_to_shift
));
1782 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1783 if (scm_is_false (scm_negative_p (n
)))
1784 return scm_quotient (n
, div
);
1786 return scm_sum (SCM_I_MAKINUM (-1L),
1787 scm_quotient (scm_sum (SCM_I_MAKINUM (1L), n
), div
));
1790 /* Shift left is done by multiplication with 2^CNT */
1791 return scm_product (n
, scm_integer_expt (SCM_I_MAKINUM (2), cnt
));
1796 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1797 (SCM n
, SCM start
, SCM end
),
1798 "Return the integer composed of the @var{start} (inclusive)\n"
1799 "through @var{end} (exclusive) bits of @var{n}. The\n"
1800 "@var{start}th bit becomes the 0-th bit in the result.\n"
1803 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1804 " @result{} \"1010\"\n"
1805 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1806 " @result{} \"10110\"\n"
1808 #define FUNC_NAME s_scm_bit_extract
1810 unsigned long int istart
, iend
, bits
;
1811 istart
= scm_to_ulong (start
);
1812 iend
= scm_to_ulong (end
);
1813 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1815 /* how many bits to keep */
1816 bits
= iend
- istart
;
1818 if (SCM_I_INUMP (n
))
1820 long int in
= SCM_I_INUM (n
);
1822 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1823 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1824 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1826 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1828 /* Since we emulate two's complement encoded numbers, this
1829 * special case requires us to produce a result that has
1830 * more bits than can be stored in a fixnum.
1832 SCM result
= scm_i_long2big (in
);
1833 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1838 /* mask down to requisite bits */
1839 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1840 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1842 else if (SCM_BIGP (n
))
1847 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1851 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1852 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1853 such bits into a ulong. */
1854 result
= scm_i_mkbig ();
1855 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1856 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1857 result
= scm_i_normbig (result
);
1859 scm_remember_upto_here_1 (n
);
1863 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1868 static const char scm_logtab
[] = {
1869 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1872 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1874 "Return the number of bits in integer @var{n}. If integer is\n"
1875 "positive, the 1-bits in its binary representation are counted.\n"
1876 "If negative, the 0-bits in its two's-complement binary\n"
1877 "representation are counted. If 0, 0 is returned.\n"
1880 "(logcount #b10101010)\n"
1887 #define FUNC_NAME s_scm_logcount
1889 if (SCM_I_INUMP (n
))
1891 unsigned long int c
= 0;
1892 long int nn
= SCM_I_INUM (n
);
1897 c
+= scm_logtab
[15 & nn
];
1900 return SCM_I_MAKINUM (c
);
1902 else if (SCM_BIGP (n
))
1904 unsigned long count
;
1905 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1906 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1908 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1909 scm_remember_upto_here_1 (n
);
1910 return SCM_I_MAKINUM (count
);
1913 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1918 static const char scm_ilentab
[] = {
1919 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1923 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1925 "Return the number of bits necessary to represent @var{n}.\n"
1928 "(integer-length #b10101010)\n"
1930 "(integer-length 0)\n"
1932 "(integer-length #b1111)\n"
1935 #define FUNC_NAME s_scm_integer_length
1937 if (SCM_I_INUMP (n
))
1939 unsigned long int c
= 0;
1941 long int nn
= SCM_I_INUM (n
);
1947 l
= scm_ilentab
[15 & nn
];
1950 return SCM_I_MAKINUM (c
- 4 + l
);
1952 else if (SCM_BIGP (n
))
1954 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1955 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1956 1 too big, so check for that and adjust. */
1957 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1958 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1959 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1960 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1962 scm_remember_upto_here_1 (n
);
1963 return SCM_I_MAKINUM (size
);
1966 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1970 /*** NUMBERS -> STRINGS ***/
1971 #define SCM_MAX_DBL_PREC 60
1972 #define SCM_MAX_DBL_RADIX 36
1974 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
1975 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
1976 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
1979 void init_dblprec(int *prec
, int radix
) {
1980 /* determine floating point precision by adding successively
1981 smaller increments to 1.0 until it is considered == 1.0 */
1982 double f
= ((double)1.0)/radix
;
1983 double fsum
= 1.0 + f
;
1988 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2000 void init_fx_radix(double *fx_list
, int radix
)
2002 /* initialize a per-radix list of tolerances. When added
2003 to a number < 1.0, we can determine if we should raund
2004 up and quit converting a number to a string. */
2008 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2009 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2012 /* use this array as a way to generate a single digit */
2013 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2016 idbl2str (double f
, char *a
, int radix
)
2018 int efmt
, dpt
, d
, i
, wp
;
2020 #ifdef DBL_MIN_10_EXP
2023 #endif /* DBL_MIN_10_EXP */
2028 radix
> SCM_MAX_DBL_RADIX
)
2030 /* revert to existing behavior */
2034 wp
= scm_dblprec
[radix
-2];
2035 fx
= fx_per_radix
[radix
-2];
2039 #ifdef HAVE_COPYSIGN
2040 double sgn
= copysign (1.0, f
);
2045 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2051 strcpy (a
, "-inf.0");
2053 strcpy (a
, "+inf.0");
2056 else if (xisnan (f
))
2058 strcpy (a
, "+nan.0");
2068 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2069 make-uniform-vector, from causing infinite loops. */
2070 /* just do the checking...if it passes, we do the conversion for our
2071 radix again below */
2078 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2086 while (f_cpy
> 10.0)
2089 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2110 if (f
+ fx
[wp
] >= radix
)
2117 /* adding 9999 makes this equivalent to abs(x) % 3 */
2118 dpt
= (exp
+ 9999) % 3;
2122 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2144 a
[ch
++] = number_chars
[d
];
2147 if (f
+ fx
[wp
] >= 1.0)
2149 a
[ch
- 1] = number_chars
[d
+1];
2161 if ((dpt
> 4) && (exp
> 6))
2163 d
= (a
[0] == '-' ? 2 : 1);
2164 for (i
= ch
++; i
> d
; i
--)
2177 if (a
[ch
- 1] == '.')
2178 a
[ch
++] = '0'; /* trailing zero */
2187 for (i
= radix
; i
<= exp
; i
*= radix
);
2188 for (i
/= radix
; i
; i
/= radix
)
2190 a
[ch
++] = number_chars
[exp
/ i
];
2198 iflo2str (SCM flt
, char *str
, int radix
)
2201 if (SCM_REALP (flt
))
2202 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2205 i
= idbl2str (SCM_COMPLEX_REAL (flt
), str
, radix
);
2206 if (SCM_COMPLEX_IMAG (flt
) != 0.0)
2208 double imag
= SCM_COMPLEX_IMAG (flt
);
2209 /* Don't output a '+' for negative numbers or for Inf and
2210 NaN. They will provide their own sign. */
2211 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2213 i
+= idbl2str (imag
, &str
[i
], radix
);
2220 /* convert a long to a string (unterminated). returns the number of
2221 characters in the result.
2223 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2225 scm_iint2str (long num
, int rad
, char *p
)
2229 unsigned long n
= (num
< 0) ? -num
: num
;
2231 for (n
/= rad
; n
> 0; n
/= rad
)
2248 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2253 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2255 "Return a string holding the external representation of the\n"
2256 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2257 "inexact, a radix of 10 will be used.")
2258 #define FUNC_NAME s_scm_number_to_string
2262 if (SCM_UNBNDP (radix
))
2265 base
= scm_to_signed_integer (radix
, 2, 36);
2267 if (SCM_I_INUMP (n
))
2269 char num_buf
[SCM_INTBUFLEN
];
2270 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2271 return scm_mem2string (num_buf
, length
);
2273 else if (SCM_BIGP (n
))
2275 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2276 scm_remember_upto_here_1 (n
);
2277 return scm_take0str (str
);
2279 else if (SCM_FRACTIONP (n
))
2281 scm_i_fraction_reduce (n
);
2282 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2283 scm_mem2string ("/", 1),
2284 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2286 else if (SCM_INEXACTP (n
))
2288 char num_buf
[FLOBUFLEN
];
2289 return scm_mem2string (num_buf
, iflo2str (n
, num_buf
, base
));
2292 SCM_WRONG_TYPE_ARG (1, n
);
2297 /* These print routines used to be stubbed here so that scm_repl.c
2298 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2301 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2303 char num_buf
[FLOBUFLEN
];
2304 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2309 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2312 char num_buf
[FLOBUFLEN
];
2313 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2318 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2321 scm_i_fraction_reduce (sexp
);
2322 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2323 scm_lfwrite (SCM_STRING_CHARS (str
), SCM_STRING_LENGTH (str
), port
);
2324 scm_remember_upto_here_1 (str
);
2329 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2331 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2332 scm_remember_upto_here_1 (exp
);
2333 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2337 /*** END nums->strs ***/
2340 /*** STRINGS -> NUMBERS ***/
2342 /* The following functions implement the conversion from strings to numbers.
2343 * The implementation somehow follows the grammar for numbers as it is given
2344 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2345 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2346 * points should be noted about the implementation:
2347 * * Each function keeps a local index variable 'idx' that points at the
2348 * current position within the parsed string. The global index is only
2349 * updated if the function could parse the corresponding syntactic unit
2351 * * Similarly, the functions keep track of indicators of inexactness ('#',
2352 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2353 * global exactness information is only updated after each part has been
2354 * successfully parsed.
2355 * * Sequences of digits are parsed into temporary variables holding fixnums.
2356 * Only if these fixnums would overflow, the result variables are updated
2357 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2358 * the temporary variables holding the fixnums are cleared, and the process
2359 * starts over again. If for example fixnums were able to store five decimal
2360 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2361 * and the result was computed as 12345 * 100000 + 67890. In other words,
2362 * only every five digits two bignum operations were performed.
2365 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2367 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2369 /* In non ASCII-style encodings the following macro might not work. */
2370 #define XDIGIT2UINT(d) \
2371 (isdigit ((int) (unsigned char) d) \
2373 : tolower ((int) (unsigned char) d) - 'a' + 10)
2376 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2377 unsigned int radix
, enum t_exactness
*p_exactness
)
2379 unsigned int idx
= *p_idx
;
2380 unsigned int hash_seen
= 0;
2381 scm_t_bits shift
= 1;
2383 unsigned int digit_value
;
2391 if (!isxdigit ((int) (unsigned char) c
))
2393 digit_value
= XDIGIT2UINT (c
);
2394 if (digit_value
>= radix
)
2398 result
= SCM_I_MAKINUM (digit_value
);
2402 if (isxdigit ((int) (unsigned char) c
))
2406 digit_value
= XDIGIT2UINT (c
);
2407 if (digit_value
>= radix
)
2419 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2421 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2423 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2430 shift
= shift
* radix
;
2431 add
= add
* radix
+ digit_value
;
2436 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2438 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2442 *p_exactness
= INEXACT
;
2448 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2449 * covers the parts of the rules that start at a potential point. The value
2450 * of the digits up to the point have been parsed by the caller and are given
2451 * in variable result. The content of *p_exactness indicates, whether a hash
2452 * has already been seen in the digits before the point.
2455 /* In non ASCII-style encodings the following macro might not work. */
2456 #define DIGIT2UINT(d) ((d) - '0')
2459 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2460 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2462 unsigned int idx
= *p_idx
;
2463 enum t_exactness x
= *p_exactness
;
2468 if (mem
[idx
] == '.')
2470 scm_t_bits shift
= 1;
2472 unsigned int digit_value
;
2473 SCM big_shift
= SCM_I_MAKINUM (1);
2479 if (isdigit ((int) (unsigned char) c
))
2484 digit_value
= DIGIT2UINT (c
);
2495 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2497 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2498 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2500 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2508 add
= add
* 10 + digit_value
;
2514 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2515 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2516 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2519 result
= scm_divide (result
, big_shift
);
2521 /* We've seen a decimal point, thus the value is implicitly inexact. */
2533 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2560 if (!isdigit ((int) (unsigned char) c
))
2564 exponent
= DIGIT2UINT (c
);
2568 if (isdigit ((int) (unsigned char) c
))
2571 if (exponent
<= SCM_MAXEXP
)
2572 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2578 if (exponent
> SCM_MAXEXP
)
2580 size_t exp_len
= idx
- start
;
2581 SCM exp_string
= scm_mem2string (&mem
[start
], exp_len
);
2582 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2583 scm_out_of_range ("string->number", exp_num
);
2586 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2588 result
= scm_product (result
, e
);
2590 result
= scm_divide2real (result
, e
);
2592 /* We've seen an exponent, thus the value is implicitly inexact. */
2610 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2613 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2614 unsigned int radix
, enum t_exactness
*p_exactness
)
2616 unsigned int idx
= *p_idx
;
2622 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2628 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2630 enum t_exactness x
= EXACT
;
2632 /* Cobble up the fractional part. We might want to set the
2633 NaN's mantissa from it. */
2635 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2640 if (mem
[idx
] == '.')
2644 else if (idx
+ 1 == len
)
2646 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2649 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2650 p_idx
, p_exactness
);
2654 enum t_exactness x
= EXACT
;
2657 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2658 if (scm_is_false (uinteger
))
2663 else if (mem
[idx
] == '/')
2669 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2670 if (scm_is_false (divisor
))
2673 /* both are int/big here, I assume */
2674 result
= scm_make_ratio (uinteger
, divisor
);
2676 else if (radix
== 10)
2678 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2679 if (scm_is_false (result
))
2690 /* When returning an inexact zero, make sure it is represented as a
2691 floating point value so that we can change its sign.
2693 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2694 result
= scm_make_real (0.0);
2700 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2703 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2704 unsigned int radix
, enum t_exactness
*p_exactness
)
2728 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2729 if (scm_is_false (ureal
))
2731 /* input must be either +i or -i */
2736 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2742 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2749 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2750 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2759 /* either +<ureal>i or -<ureal>i */
2766 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2769 /* polar input: <real>@<real>. */
2794 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2795 if (scm_is_false (angle
))
2800 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2801 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2803 result
= scm_make_polar (ureal
, angle
);
2808 /* expecting input matching <real>[+-]<ureal>?i */
2815 int sign
= (c
== '+') ? 1 : -1;
2816 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2818 if (scm_is_false (imag
))
2819 imag
= SCM_I_MAKINUM (sign
);
2820 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2821 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2825 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2832 return scm_make_rectangular (ureal
, imag
);
2841 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2843 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2846 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2848 unsigned int idx
= 0;
2849 unsigned int radix
= NO_RADIX
;
2850 enum t_exactness forced_x
= NO_EXACTNESS
;
2851 enum t_exactness implicit_x
= EXACT
;
2854 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2855 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2857 switch (mem
[idx
+ 1])
2860 if (radix
!= NO_RADIX
)
2865 if (radix
!= NO_RADIX
)
2870 if (forced_x
!= NO_EXACTNESS
)
2875 if (forced_x
!= NO_EXACTNESS
)
2880 if (radix
!= NO_RADIX
)
2885 if (radix
!= NO_RADIX
)
2895 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2896 if (radix
== NO_RADIX
)
2897 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2899 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2901 if (scm_is_false (result
))
2907 if (SCM_INEXACTP (result
))
2908 return scm_inexact_to_exact (result
);
2912 if (SCM_INEXACTP (result
))
2915 return scm_exact_to_inexact (result
);
2918 if (implicit_x
== INEXACT
)
2920 if (SCM_INEXACTP (result
))
2923 return scm_exact_to_inexact (result
);
2931 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2932 (SCM string
, SCM radix
),
2933 "Return a number of the maximally precise representation\n"
2934 "expressed by the given @var{string}. @var{radix} must be an\n"
2935 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2936 "is a default radix that may be overridden by an explicit radix\n"
2937 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2938 "supplied, then the default radix is 10. If string is not a\n"
2939 "syntactically valid notation for a number, then\n"
2940 "@code{string->number} returns @code{#f}.")
2941 #define FUNC_NAME s_scm_string_to_number
2945 SCM_VALIDATE_STRING (1, string
);
2947 if (SCM_UNBNDP (radix
))
2950 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
2952 answer
= scm_i_mem2number (SCM_STRING_CHARS (string
),
2953 SCM_STRING_LENGTH (string
),
2955 return scm_return_first (answer
, string
);
2960 /*** END strs->nums ***/
2964 scm_make_real (double x
)
2966 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
2968 SCM_REAL_VALUE (z
) = x
;
2974 scm_make_complex (double x
, double y
)
2977 return scm_make_real (x
);
2981 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
2983 SCM_COMPLEX_REAL (z
) = x
;
2984 SCM_COMPLEX_IMAG (z
) = y
;
2991 scm_bigequal (SCM x
, SCM y
)
2993 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2994 scm_remember_upto_here_2 (x
, y
);
2995 return scm_from_bool (0 == result
);
2999 scm_real_equalp (SCM x
, SCM y
)
3001 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3005 scm_complex_equalp (SCM x
, SCM y
)
3007 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3008 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3012 scm_i_fraction_equalp (SCM x
, SCM y
)
3014 scm_i_fraction_reduce (x
);
3015 scm_i_fraction_reduce (y
);
3016 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3017 SCM_FRACTION_NUMERATOR (y
)))
3018 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3019 SCM_FRACTION_DENOMINATOR (y
))))
3026 SCM_REGISTER_PROC (s_number_p
, "number?", 1, 0, 0, scm_number_p
);
3027 /* "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3028 * "else. Note that the sets of complex, real, rational and\n"
3029 * "integer values form subsets of the set of numbers, i. e. the\n"
3030 * "predicate will be fulfilled for any number."
3032 SCM_DEFINE (scm_number_p
, "complex?", 1, 0, 0,
3034 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3035 "otherwise. Note that the sets of real, rational and integer\n"
3036 "values form subsets of the set of complex numbers, i. e. the\n"
3037 "predicate will also be fulfilled if @var{x} is a real,\n"
3038 "rational or integer number.")
3039 #define FUNC_NAME s_scm_number_p
3041 return scm_from_bool (SCM_NUMBERP (x
));
3046 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3048 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3049 "otherwise. Note that the set of integer values forms a subset of\n"
3050 "the set of real numbers, i. e. the predicate will also be\n"
3051 "fulfilled if @var{x} is an integer number.")
3052 #define FUNC_NAME s_scm_real_p
3054 /* we can't represent irrational numbers. */
3055 return scm_rational_p (x
);
3059 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3061 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3062 "otherwise. Note that the set of integer values forms a subset of\n"
3063 "the set of rational numbers, i. e. the predicate will also be\n"
3064 "fulfilled if @var{x} is an integer number.")
3065 #define FUNC_NAME s_scm_rational_p
3067 if (SCM_I_INUMP (x
))
3069 else if (SCM_IMP (x
))
3071 else if (SCM_BIGP (x
))
3073 else if (SCM_FRACTIONP (x
))
3075 else if (SCM_REALP (x
))
3076 /* due to their limited precision, all floating point numbers are
3077 rational as well. */
3085 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3087 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3089 #define FUNC_NAME s_scm_integer_p
3092 if (SCM_I_INUMP (x
))
3098 if (!SCM_INEXACTP (x
))
3100 if (SCM_COMPLEXP (x
))
3102 r
= SCM_REAL_VALUE (x
);
3110 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3112 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3114 #define FUNC_NAME s_scm_inexact_p
3116 if (SCM_INEXACTP (x
))
3118 if (SCM_NUMBERP (x
))
3120 SCM_WRONG_TYPE_ARG (1, x
);
3125 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3126 /* "Return @code{#t} if all parameters are numerically equal." */
3128 scm_num_eq_p (SCM x
, SCM y
)
3131 if (SCM_I_INUMP (x
))
3133 long xx
= SCM_I_INUM (x
);
3134 if (SCM_I_INUMP (y
))
3136 long yy
= SCM_I_INUM (y
);
3137 return scm_from_bool (xx
== yy
);
3139 else if (SCM_BIGP (y
))
3141 else if (SCM_REALP (y
))
3142 return scm_from_bool ((double) xx
== SCM_REAL_VALUE (y
));
3143 else if (SCM_COMPLEXP (y
))
3144 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3145 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3146 else if (SCM_FRACTIONP (y
))
3149 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3151 else if (SCM_BIGP (x
))
3153 if (SCM_I_INUMP (y
))
3155 else if (SCM_BIGP (y
))
3157 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3158 scm_remember_upto_here_2 (x
, y
);
3159 return scm_from_bool (0 == cmp
);
3161 else if (SCM_REALP (y
))
3164 if (xisnan (SCM_REAL_VALUE (y
)))
3166 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3167 scm_remember_upto_here_1 (x
);
3168 return scm_from_bool (0 == cmp
);
3170 else if (SCM_COMPLEXP (y
))
3173 if (0.0 != SCM_COMPLEX_IMAG (y
))
3175 if (xisnan (SCM_COMPLEX_REAL (y
)))
3177 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3178 scm_remember_upto_here_1 (x
);
3179 return scm_from_bool (0 == cmp
);
3181 else if (SCM_FRACTIONP (y
))
3184 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3186 else if (SCM_REALP (x
))
3188 if (SCM_I_INUMP (y
))
3189 return scm_from_bool (SCM_REAL_VALUE (x
) == (double) SCM_I_INUM (y
));
3190 else if (SCM_BIGP (y
))
3193 if (xisnan (SCM_REAL_VALUE (x
)))
3195 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3196 scm_remember_upto_here_1 (y
);
3197 return scm_from_bool (0 == cmp
);
3199 else if (SCM_REALP (y
))
3200 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3201 else if (SCM_COMPLEXP (y
))
3202 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3203 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3204 else if (SCM_FRACTIONP (y
))
3206 double xx
= SCM_REAL_VALUE (x
);
3210 return scm_from_bool (xx
< 0.0);
3211 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3215 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3217 else if (SCM_COMPLEXP (x
))
3219 if (SCM_I_INUMP (y
))
3220 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3221 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3222 else if (SCM_BIGP (y
))
3225 if (0.0 != SCM_COMPLEX_IMAG (x
))
3227 if (xisnan (SCM_COMPLEX_REAL (x
)))
3229 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3230 scm_remember_upto_here_1 (y
);
3231 return scm_from_bool (0 == cmp
);
3233 else if (SCM_REALP (y
))
3234 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3235 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3236 else if (SCM_COMPLEXP (y
))
3237 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3238 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3239 else if (SCM_FRACTIONP (y
))
3242 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3244 xx
= SCM_COMPLEX_REAL (x
);
3248 return scm_from_bool (xx
< 0.0);
3249 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3253 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3255 else if (SCM_FRACTIONP (x
))
3257 if (SCM_I_INUMP (y
))
3259 else if (SCM_BIGP (y
))
3261 else if (SCM_REALP (y
))
3263 double yy
= SCM_REAL_VALUE (y
);
3267 return scm_from_bool (0.0 < yy
);
3268 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3271 else if (SCM_COMPLEXP (y
))
3274 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3276 yy
= SCM_COMPLEX_REAL (y
);
3280 return scm_from_bool (0.0 < yy
);
3281 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3284 else if (SCM_FRACTIONP (y
))
3285 return scm_i_fraction_equalp (x
, y
);
3287 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3290 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3294 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3295 done are good for inums, but for bignums an answer can almost always be
3296 had by just examining a few high bits of the operands, as done by GMP in
3297 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3298 of the float exponent to take into account. */
3300 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3301 /* "Return @code{#t} if the list of parameters is monotonically\n"
3305 scm_less_p (SCM x
, SCM y
)
3308 if (SCM_I_INUMP (x
))
3310 long xx
= SCM_I_INUM (x
);
3311 if (SCM_I_INUMP (y
))
3313 long yy
= SCM_I_INUM (y
);
3314 return scm_from_bool (xx
< yy
);
3316 else if (SCM_BIGP (y
))
3318 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3319 scm_remember_upto_here_1 (y
);
3320 return scm_from_bool (sgn
> 0);
3322 else if (SCM_REALP (y
))
3323 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3324 else if (SCM_FRACTIONP (y
))
3326 /* "x < a/b" becomes "x*b < a" */
3328 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3329 y
= SCM_FRACTION_NUMERATOR (y
);
3333 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3335 else if (SCM_BIGP (x
))
3337 if (SCM_I_INUMP (y
))
3339 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3340 scm_remember_upto_here_1 (x
);
3341 return scm_from_bool (sgn
< 0);
3343 else if (SCM_BIGP (y
))
3345 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3346 scm_remember_upto_here_2 (x
, y
);
3347 return scm_from_bool (cmp
< 0);
3349 else if (SCM_REALP (y
))
3352 if (xisnan (SCM_REAL_VALUE (y
)))
3354 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3355 scm_remember_upto_here_1 (x
);
3356 return scm_from_bool (cmp
< 0);
3358 else if (SCM_FRACTIONP (y
))
3361 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3363 else if (SCM_REALP (x
))
3365 if (SCM_I_INUMP (y
))
3366 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3367 else if (SCM_BIGP (y
))
3370 if (xisnan (SCM_REAL_VALUE (x
)))
3372 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3373 scm_remember_upto_here_1 (y
);
3374 return scm_from_bool (cmp
> 0);
3376 else if (SCM_REALP (y
))
3377 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3378 else if (SCM_FRACTIONP (y
))
3380 double xx
= SCM_REAL_VALUE (x
);
3384 return scm_from_bool (xx
< 0.0);
3385 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3389 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3391 else if (SCM_FRACTIONP (x
))
3393 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3395 /* "a/b < y" becomes "a < y*b" */
3396 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3397 x
= SCM_FRACTION_NUMERATOR (x
);
3400 else if (SCM_REALP (y
))
3402 double yy
= SCM_REAL_VALUE (y
);
3406 return scm_from_bool (0.0 < yy
);
3407 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3410 else if (SCM_FRACTIONP (y
))
3412 /* "a/b < c/d" becomes "a*d < c*b" */
3413 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3414 SCM_FRACTION_DENOMINATOR (y
));
3415 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3416 SCM_FRACTION_DENOMINATOR (x
));
3422 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3425 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3429 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3430 /* "Return @code{#t} if the list of parameters is monotonically\n"
3433 #define FUNC_NAME s_scm_gr_p
3435 scm_gr_p (SCM x
, SCM y
)
3437 if (!SCM_NUMBERP (x
))
3438 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3439 else if (!SCM_NUMBERP (y
))
3440 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3442 return scm_less_p (y
, x
);
3447 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3448 /* "Return @code{#t} if the list of parameters is monotonically\n"
3451 #define FUNC_NAME s_scm_leq_p
3453 scm_leq_p (SCM x
, SCM y
)
3455 if (!SCM_NUMBERP (x
))
3456 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3457 else if (!SCM_NUMBERP (y
))
3458 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3459 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3462 return scm_not (scm_less_p (y
, x
));
3467 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3468 /* "Return @code{#t} if the list of parameters is monotonically\n"
3471 #define FUNC_NAME s_scm_geq_p
3473 scm_geq_p (SCM x
, SCM y
)
3475 if (!SCM_NUMBERP (x
))
3476 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3477 else if (!SCM_NUMBERP (y
))
3478 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3479 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3482 return scm_not (scm_less_p (x
, y
));
3487 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3488 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3494 if (SCM_I_INUMP (z
))
3495 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3496 else if (SCM_BIGP (z
))
3498 else if (SCM_REALP (z
))
3499 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3500 else if (SCM_COMPLEXP (z
))
3501 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3502 && SCM_COMPLEX_IMAG (z
) == 0.0);
3503 else if (SCM_FRACTIONP (z
))
3506 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3510 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3511 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3515 scm_positive_p (SCM x
)
3517 if (SCM_I_INUMP (x
))
3518 return scm_from_bool (SCM_I_INUM (x
) > 0);
3519 else if (SCM_BIGP (x
))
3521 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3522 scm_remember_upto_here_1 (x
);
3523 return scm_from_bool (sgn
> 0);
3525 else if (SCM_REALP (x
))
3526 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3527 else if (SCM_FRACTIONP (x
))
3528 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3530 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3534 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3535 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3539 scm_negative_p (SCM x
)
3541 if (SCM_I_INUMP (x
))
3542 return scm_from_bool (SCM_I_INUM (x
) < 0);
3543 else if (SCM_BIGP (x
))
3545 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3546 scm_remember_upto_here_1 (x
);
3547 return scm_from_bool (sgn
< 0);
3549 else if (SCM_REALP (x
))
3550 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3551 else if (SCM_FRACTIONP (x
))
3552 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3554 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3558 /* scm_min and scm_max return an inexact when either argument is inexact, as
3559 required by r5rs. On that basis, for exact/inexact combinations the
3560 exact is converted to inexact to compare and possibly return. This is
3561 unlike scm_less_p above which takes some trouble to preserve all bits in
3562 its test, such trouble is not required for min and max. */
3564 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3565 /* "Return the maximum of all parameter values."
3568 scm_max (SCM x
, SCM y
)
3573 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3574 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3577 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3580 if (SCM_I_INUMP (x
))
3582 long xx
= SCM_I_INUM (x
);
3583 if (SCM_I_INUMP (y
))
3585 long yy
= SCM_I_INUM (y
);
3586 return (xx
< yy
) ? y
: x
;
3588 else if (SCM_BIGP (y
))
3590 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3591 scm_remember_upto_here_1 (y
);
3592 return (sgn
< 0) ? x
: y
;
3594 else if (SCM_REALP (y
))
3597 /* if y==NaN then ">" is false and we return NaN */
3598 return (z
> SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3600 else if (SCM_FRACTIONP (y
))
3603 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3606 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3608 else if (SCM_BIGP (x
))
3610 if (SCM_I_INUMP (y
))
3612 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3613 scm_remember_upto_here_1 (x
);
3614 return (sgn
< 0) ? y
: x
;
3616 else if (SCM_BIGP (y
))
3618 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3619 scm_remember_upto_here_2 (x
, y
);
3620 return (cmp
> 0) ? x
: y
;
3622 else if (SCM_REALP (y
))
3624 /* if y==NaN then xx>yy is false, so we return the NaN y */
3627 xx
= scm_i_big2dbl (x
);
3628 yy
= SCM_REAL_VALUE (y
);
3629 return (xx
> yy
? scm_make_real (xx
) : y
);
3631 else if (SCM_FRACTIONP (y
))
3636 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3638 else if (SCM_REALP (x
))
3640 if (SCM_I_INUMP (y
))
3642 double z
= SCM_I_INUM (y
);
3643 /* if x==NaN then "<" is false and we return NaN */
3644 return (SCM_REAL_VALUE (x
) < z
) ? scm_make_real (z
) : x
;
3646 else if (SCM_BIGP (y
))
3651 else if (SCM_REALP (y
))
3653 /* if x==NaN then our explicit check means we return NaN
3654 if y==NaN then ">" is false and we return NaN
3655 calling isnan is unavoidable, since it's the only way to know
3656 which of x or y causes any compares to be false */
3657 double xx
= SCM_REAL_VALUE (x
);
3658 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3660 else if (SCM_FRACTIONP (y
))
3662 double yy
= scm_i_fraction2double (y
);
3663 double xx
= SCM_REAL_VALUE (x
);
3664 return (xx
< yy
) ? scm_make_real (yy
) : x
;
3667 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3669 else if (SCM_FRACTIONP (x
))
3671 if (SCM_I_INUMP (y
))
3675 else if (SCM_BIGP (y
))
3679 else if (SCM_REALP (y
))
3681 double xx
= scm_i_fraction2double (x
);
3682 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_make_real (xx
);
3684 else if (SCM_FRACTIONP (y
))
3689 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3692 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3696 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3697 /* "Return the minium of all parameter values."
3700 scm_min (SCM x
, SCM y
)
3705 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3706 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3709 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3712 if (SCM_I_INUMP (x
))
3714 long xx
= SCM_I_INUM (x
);
3715 if (SCM_I_INUMP (y
))
3717 long yy
= SCM_I_INUM (y
);
3718 return (xx
< yy
) ? x
: y
;
3720 else if (SCM_BIGP (y
))
3722 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3723 scm_remember_upto_here_1 (y
);
3724 return (sgn
< 0) ? y
: x
;
3726 else if (SCM_REALP (y
))
3729 /* if y==NaN then "<" is false and we return NaN */
3730 return (z
< SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3732 else if (SCM_FRACTIONP (y
))
3735 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3738 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3740 else if (SCM_BIGP (x
))
3742 if (SCM_I_INUMP (y
))
3744 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3745 scm_remember_upto_here_1 (x
);
3746 return (sgn
< 0) ? x
: y
;
3748 else if (SCM_BIGP (y
))
3750 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3751 scm_remember_upto_here_2 (x
, y
);
3752 return (cmp
> 0) ? y
: x
;
3754 else if (SCM_REALP (y
))
3756 /* if y==NaN then xx<yy is false, so we return the NaN y */
3759 xx
= scm_i_big2dbl (x
);
3760 yy
= SCM_REAL_VALUE (y
);
3761 return (xx
< yy
? scm_make_real (xx
) : y
);
3763 else if (SCM_FRACTIONP (y
))
3768 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3770 else if (SCM_REALP (x
))
3772 if (SCM_I_INUMP (y
))
3774 double z
= SCM_I_INUM (y
);
3775 /* if x==NaN then "<" is false and we return NaN */
3776 return (z
< SCM_REAL_VALUE (x
)) ? scm_make_real (z
) : x
;
3778 else if (SCM_BIGP (y
))
3783 else if (SCM_REALP (y
))
3785 /* if x==NaN then our explicit check means we return NaN
3786 if y==NaN then "<" is false and we return NaN
3787 calling isnan is unavoidable, since it's the only way to know
3788 which of x or y causes any compares to be false */
3789 double xx
= SCM_REAL_VALUE (x
);
3790 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3792 else if (SCM_FRACTIONP (y
))
3794 double yy
= scm_i_fraction2double (y
);
3795 double xx
= SCM_REAL_VALUE (x
);
3796 return (yy
< xx
) ? scm_make_real (yy
) : x
;
3799 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3801 else if (SCM_FRACTIONP (x
))
3803 if (SCM_I_INUMP (y
))
3807 else if (SCM_BIGP (y
))
3811 else if (SCM_REALP (y
))
3813 double xx
= scm_i_fraction2double (x
);
3814 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_make_real (xx
);
3816 else if (SCM_FRACTIONP (y
))
3821 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3824 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3828 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3829 /* "Return the sum of all parameter values. Return 0 if called without\n"
3833 scm_sum (SCM x
, SCM y
)
3837 if (SCM_NUMBERP (x
)) return x
;
3838 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3839 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3842 if (SCM_I_INUMP (x
))
3844 if (SCM_I_INUMP (y
))
3846 long xx
= SCM_I_INUM (x
);
3847 long yy
= SCM_I_INUM (y
);
3848 long int z
= xx
+ yy
;
3849 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3851 else if (SCM_BIGP (y
))
3856 else if (SCM_REALP (y
))
3858 long int xx
= SCM_I_INUM (x
);
3859 return scm_make_real (xx
+ SCM_REAL_VALUE (y
));
3861 else if (SCM_COMPLEXP (y
))
3863 long int xx
= SCM_I_INUM (x
);
3864 return scm_make_complex (xx
+ SCM_COMPLEX_REAL (y
),
3865 SCM_COMPLEX_IMAG (y
));
3867 else if (SCM_FRACTIONP (y
))
3868 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3869 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3870 SCM_FRACTION_DENOMINATOR (y
));
3872 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3873 } else if (SCM_BIGP (x
))
3875 if (SCM_I_INUMP (y
))
3880 inum
= SCM_I_INUM (y
);
3883 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3886 SCM result
= scm_i_mkbig ();
3887 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3888 scm_remember_upto_here_1 (x
);
3889 /* we know the result will have to be a bignum */
3892 return scm_i_normbig (result
);
3896 SCM result
= scm_i_mkbig ();
3897 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3898 scm_remember_upto_here_1 (x
);
3899 /* we know the result will have to be a bignum */
3902 return scm_i_normbig (result
);
3905 else if (SCM_BIGP (y
))
3907 SCM result
= scm_i_mkbig ();
3908 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3909 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3910 mpz_add (SCM_I_BIG_MPZ (result
),
3913 scm_remember_upto_here_2 (x
, y
);
3914 /* we know the result will have to be a bignum */
3917 return scm_i_normbig (result
);
3919 else if (SCM_REALP (y
))
3921 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3922 scm_remember_upto_here_1 (x
);
3923 return scm_make_real (result
);
3925 else if (SCM_COMPLEXP (y
))
3927 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3928 + SCM_COMPLEX_REAL (y
));
3929 scm_remember_upto_here_1 (x
);
3930 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
3932 else if (SCM_FRACTIONP (y
))
3933 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3934 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3935 SCM_FRACTION_DENOMINATOR (y
));
3937 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3939 else if (SCM_REALP (x
))
3941 if (SCM_I_INUMP (y
))
3942 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
3943 else if (SCM_BIGP (y
))
3945 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3946 scm_remember_upto_here_1 (y
);
3947 return scm_make_real (result
);
3949 else if (SCM_REALP (y
))
3950 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3951 else if (SCM_COMPLEXP (y
))
3952 return scm_make_complex (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3953 SCM_COMPLEX_IMAG (y
));
3954 else if (SCM_FRACTIONP (y
))
3955 return scm_make_real (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3957 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3959 else if (SCM_COMPLEXP (x
))
3961 if (SCM_I_INUMP (y
))
3962 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
3963 SCM_COMPLEX_IMAG (x
));
3964 else if (SCM_BIGP (y
))
3966 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3967 + SCM_COMPLEX_REAL (x
));
3968 scm_remember_upto_here_1 (y
);
3969 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (x
));
3971 else if (SCM_REALP (y
))
3972 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3973 SCM_COMPLEX_IMAG (x
));
3974 else if (SCM_COMPLEXP (y
))
3975 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3976 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3977 else if (SCM_FRACTIONP (y
))
3978 return scm_make_complex (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3979 SCM_COMPLEX_IMAG (x
));
3981 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3983 else if (SCM_FRACTIONP (x
))
3985 if (SCM_I_INUMP (y
))
3986 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3987 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3988 SCM_FRACTION_DENOMINATOR (x
));
3989 else if (SCM_BIGP (y
))
3990 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3991 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3992 SCM_FRACTION_DENOMINATOR (x
));
3993 else if (SCM_REALP (y
))
3994 return scm_make_real (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
3995 else if (SCM_COMPLEXP (y
))
3996 return scm_make_complex (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
3997 SCM_COMPLEX_IMAG (y
));
3998 else if (SCM_FRACTIONP (y
))
3999 /* a/b + c/d = (ad + bc) / bd */
4000 return scm_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4001 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4002 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4004 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4007 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4011 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4012 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4013 * the sum of all but the first argument are subtracted from the first
4015 #define FUNC_NAME s_difference
4017 scm_difference (SCM x
, SCM y
)
4022 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4024 if (SCM_I_INUMP (x
))
4026 long xx
= -SCM_I_INUM (x
);
4027 if (SCM_FIXABLE (xx
))
4028 return SCM_I_MAKINUM (xx
);
4030 return scm_i_long2big (xx
);
4032 else if (SCM_BIGP (x
))
4033 /* FIXME: do we really need to normalize here? */
4034 return scm_i_normbig (scm_i_clonebig (x
, 0));
4035 else if (SCM_REALP (x
))
4036 return scm_make_real (-SCM_REAL_VALUE (x
));
4037 else if (SCM_COMPLEXP (x
))
4038 return scm_make_complex (-SCM_COMPLEX_REAL (x
),
4039 -SCM_COMPLEX_IMAG (x
));
4040 else if (SCM_FRACTIONP (x
))
4041 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4042 SCM_FRACTION_DENOMINATOR (x
));
4044 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4047 if (SCM_I_INUMP (x
))
4049 if (SCM_I_INUMP (y
))
4051 long int xx
= SCM_I_INUM (x
);
4052 long int yy
= SCM_I_INUM (y
);
4053 long int z
= xx
- yy
;
4054 if (SCM_FIXABLE (z
))
4055 return SCM_I_MAKINUM (z
);
4057 return scm_i_long2big (z
);
4059 else if (SCM_BIGP (y
))
4061 /* inum-x - big-y */
4062 long xx
= SCM_I_INUM (x
);
4065 return scm_i_clonebig (y
, 0);
4068 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4069 SCM result
= scm_i_mkbig ();
4072 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4075 /* x - y == -(y + -x) */
4076 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4077 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4079 scm_remember_upto_here_1 (y
);
4081 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4082 /* we know the result will have to be a bignum */
4085 return scm_i_normbig (result
);
4088 else if (SCM_REALP (y
))
4090 long int xx
= SCM_I_INUM (x
);
4091 return scm_make_real (xx
- SCM_REAL_VALUE (y
));
4093 else if (SCM_COMPLEXP (y
))
4095 long int xx
= SCM_I_INUM (x
);
4096 return scm_make_complex (xx
- SCM_COMPLEX_REAL (y
),
4097 - SCM_COMPLEX_IMAG (y
));
4099 else if (SCM_FRACTIONP (y
))
4100 /* a - b/c = (ac - b) / c */
4101 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4102 SCM_FRACTION_NUMERATOR (y
)),
4103 SCM_FRACTION_DENOMINATOR (y
));
4105 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4107 else if (SCM_BIGP (x
))
4109 if (SCM_I_INUMP (y
))
4111 /* big-x - inum-y */
4112 long yy
= SCM_I_INUM (y
);
4113 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4115 scm_remember_upto_here_1 (x
);
4117 return SCM_FIXABLE (-yy
) ? SCM_I_MAKINUM (-yy
) : scm_long2num (-yy
);
4120 SCM result
= scm_i_mkbig ();
4123 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4125 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4126 scm_remember_upto_here_1 (x
);
4128 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4129 /* we know the result will have to be a bignum */
4132 return scm_i_normbig (result
);
4135 else if (SCM_BIGP (y
))
4137 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4138 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4139 SCM result
= scm_i_mkbig ();
4140 mpz_sub (SCM_I_BIG_MPZ (result
),
4143 scm_remember_upto_here_2 (x
, y
);
4144 /* we know the result will have to be a bignum */
4145 if ((sgn_x
== 1) && (sgn_y
== -1))
4147 if ((sgn_x
== -1) && (sgn_y
== 1))
4149 return scm_i_normbig (result
);
4151 else if (SCM_REALP (y
))
4153 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4154 scm_remember_upto_here_1 (x
);
4155 return scm_make_real (result
);
4157 else if (SCM_COMPLEXP (y
))
4159 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4160 - SCM_COMPLEX_REAL (y
));
4161 scm_remember_upto_here_1 (x
);
4162 return scm_make_complex (real_part
, - SCM_COMPLEX_IMAG (y
));
4164 else if (SCM_FRACTIONP (y
))
4165 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4166 SCM_FRACTION_NUMERATOR (y
)),
4167 SCM_FRACTION_DENOMINATOR (y
));
4168 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4170 else if (SCM_REALP (x
))
4172 if (SCM_I_INUMP (y
))
4173 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4174 else if (SCM_BIGP (y
))
4176 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4177 scm_remember_upto_here_1 (x
);
4178 return scm_make_real (result
);
4180 else if (SCM_REALP (y
))
4181 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4182 else if (SCM_COMPLEXP (y
))
4183 return scm_make_complex (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4184 -SCM_COMPLEX_IMAG (y
));
4185 else if (SCM_FRACTIONP (y
))
4186 return scm_make_real (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4188 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4190 else if (SCM_COMPLEXP (x
))
4192 if (SCM_I_INUMP (y
))
4193 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4194 SCM_COMPLEX_IMAG (x
));
4195 else if (SCM_BIGP (y
))
4197 double real_part
= (SCM_COMPLEX_REAL (x
)
4198 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4199 scm_remember_upto_here_1 (x
);
4200 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
4202 else if (SCM_REALP (y
))
4203 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4204 SCM_COMPLEX_IMAG (x
));
4205 else if (SCM_COMPLEXP (y
))
4206 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4207 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4208 else if (SCM_FRACTIONP (y
))
4209 return scm_make_complex (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4210 SCM_COMPLEX_IMAG (x
));
4212 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4214 else if (SCM_FRACTIONP (x
))
4216 if (SCM_I_INUMP (y
))
4217 /* a/b - c = (a - cb) / b */
4218 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4219 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4220 SCM_FRACTION_DENOMINATOR (x
));
4221 else if (SCM_BIGP (y
))
4222 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4223 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4224 SCM_FRACTION_DENOMINATOR (x
));
4225 else if (SCM_REALP (y
))
4226 return scm_make_real (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4227 else if (SCM_COMPLEXP (y
))
4228 return scm_make_complex (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4229 -SCM_COMPLEX_IMAG (y
));
4230 else if (SCM_FRACTIONP (y
))
4231 /* a/b - c/d = (ad - bc) / bd */
4232 return scm_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4233 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4234 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4236 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4239 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4244 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4245 /* "Return the product of all arguments. If called without arguments,\n"
4249 scm_product (SCM x
, SCM y
)
4254 return SCM_I_MAKINUM (1L);
4255 else if (SCM_NUMBERP (x
))
4258 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4261 if (SCM_I_INUMP (x
))
4266 xx
= SCM_I_INUM (x
);
4270 case 0: return x
; break;
4271 case 1: return y
; break;
4274 if (SCM_I_INUMP (y
))
4276 long yy
= SCM_I_INUM (y
);
4278 SCM k
= SCM_I_MAKINUM (kk
);
4279 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4283 SCM result
= scm_i_long2big (xx
);
4284 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4285 return scm_i_normbig (result
);
4288 else if (SCM_BIGP (y
))
4290 SCM result
= scm_i_mkbig ();
4291 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4292 scm_remember_upto_here_1 (y
);
4295 else if (SCM_REALP (y
))
4296 return scm_make_real (xx
* SCM_REAL_VALUE (y
));
4297 else if (SCM_COMPLEXP (y
))
4298 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4299 xx
* SCM_COMPLEX_IMAG (y
));
4300 else if (SCM_FRACTIONP (y
))
4301 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4302 SCM_FRACTION_DENOMINATOR (y
));
4304 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4306 else if (SCM_BIGP (x
))
4308 if (SCM_I_INUMP (y
))
4313 else if (SCM_BIGP (y
))
4315 SCM result
= scm_i_mkbig ();
4316 mpz_mul (SCM_I_BIG_MPZ (result
),
4319 scm_remember_upto_here_2 (x
, y
);
4322 else if (SCM_REALP (y
))
4324 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4325 scm_remember_upto_here_1 (x
);
4326 return scm_make_real (result
);
4328 else if (SCM_COMPLEXP (y
))
4330 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4331 scm_remember_upto_here_1 (x
);
4332 return scm_make_complex (z
* SCM_COMPLEX_REAL (y
),
4333 z
* SCM_COMPLEX_IMAG (y
));
4335 else if (SCM_FRACTIONP (y
))
4336 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4337 SCM_FRACTION_DENOMINATOR (y
));
4339 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4341 else if (SCM_REALP (x
))
4343 if (SCM_I_INUMP (y
))
4344 return scm_make_real (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4345 else if (SCM_BIGP (y
))
4347 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4348 scm_remember_upto_here_1 (y
);
4349 return scm_make_real (result
);
4351 else if (SCM_REALP (y
))
4352 return scm_make_real (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4353 else if (SCM_COMPLEXP (y
))
4354 return scm_make_complex (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4355 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4356 else if (SCM_FRACTIONP (y
))
4357 return scm_make_real (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4359 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4361 else if (SCM_COMPLEXP (x
))
4363 if (SCM_I_INUMP (y
))
4364 return scm_make_complex (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4365 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4366 else if (SCM_BIGP (y
))
4368 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4369 scm_remember_upto_here_1 (y
);
4370 return scm_make_complex (z
* SCM_COMPLEX_REAL (x
),
4371 z
* SCM_COMPLEX_IMAG (x
));
4373 else if (SCM_REALP (y
))
4374 return scm_make_complex (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4375 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4376 else if (SCM_COMPLEXP (y
))
4378 return scm_make_complex (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4379 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4380 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4381 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4383 else if (SCM_FRACTIONP (y
))
4385 double yy
= scm_i_fraction2double (y
);
4386 return scm_make_complex (yy
* SCM_COMPLEX_REAL (x
),
4387 yy
* SCM_COMPLEX_IMAG (x
));
4390 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4392 else if (SCM_FRACTIONP (x
))
4394 if (SCM_I_INUMP (y
))
4395 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4396 SCM_FRACTION_DENOMINATOR (x
));
4397 else if (SCM_BIGP (y
))
4398 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4399 SCM_FRACTION_DENOMINATOR (x
));
4400 else if (SCM_REALP (y
))
4401 return scm_make_real (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4402 else if (SCM_COMPLEXP (y
))
4404 double xx
= scm_i_fraction2double (x
);
4405 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4406 xx
* SCM_COMPLEX_IMAG (y
));
4408 else if (SCM_FRACTIONP (y
))
4409 /* a/b * c/d = ac / bd */
4410 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4411 SCM_FRACTION_NUMERATOR (y
)),
4412 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4413 SCM_FRACTION_DENOMINATOR (y
)));
4415 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4418 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4422 scm_num2dbl (SCM a
, const char *why
)
4423 #define FUNC_NAME why
4425 if (SCM_I_INUMP (a
))
4426 return (double) SCM_I_INUM (a
);
4427 else if (SCM_BIGP (a
))
4429 double result
= mpz_get_d (SCM_I_BIG_MPZ (a
));
4430 scm_remember_upto_here_1 (a
);
4433 else if (SCM_REALP (a
))
4434 return (SCM_REAL_VALUE (a
));
4435 else if (SCM_FRACTIONP (a
))
4436 return scm_i_fraction2double (a
);
4438 SCM_WRONG_TYPE_ARG (SCM_ARGn
, a
);
4442 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4443 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4444 #define ALLOW_DIVIDE_BY_ZERO
4445 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4448 /* The code below for complex division is adapted from the GNU
4449 libstdc++, which adapted it from f2c's libF77, and is subject to
4452 /****************************************************************
4453 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4455 Permission to use, copy, modify, and distribute this software
4456 and its documentation for any purpose and without fee is hereby
4457 granted, provided that the above copyright notice appear in all
4458 copies and that both that the copyright notice and this
4459 permission notice and warranty disclaimer appear in supporting
4460 documentation, and that the names of AT&T Bell Laboratories or
4461 Bellcore or any of their entities not be used in advertising or
4462 publicity pertaining to distribution of the software without
4463 specific, written prior permission.
4465 AT&T and Bellcore disclaim all warranties with regard to this
4466 software, including all implied warranties of merchantability
4467 and fitness. In no event shall AT&T or Bellcore be liable for
4468 any special, indirect or consequential damages or any damages
4469 whatsoever resulting from loss of use, data or profits, whether
4470 in an action of contract, negligence or other tortious action,
4471 arising out of or in connection with the use or performance of
4473 ****************************************************************/
4475 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4476 /* Divide the first argument by the product of the remaining
4477 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4479 #define FUNC_NAME s_divide
4481 scm_i_divide (SCM x
, SCM y
, int inexact
)
4488 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4489 else if (SCM_I_INUMP (x
))
4491 long xx
= SCM_I_INUM (x
);
4492 if (xx
== 1 || xx
== -1)
4494 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4496 scm_num_overflow (s_divide
);
4501 return scm_make_real (1.0 / (double) xx
);
4502 else return scm_make_ratio (SCM_I_MAKINUM(1), x
);
4505 else if (SCM_BIGP (x
))
4508 return scm_make_real (1.0 / scm_i_big2dbl (x
));
4509 else return scm_make_ratio (SCM_I_MAKINUM(1), x
);
4511 else if (SCM_REALP (x
))
4513 double xx
= SCM_REAL_VALUE (x
);
4514 #ifndef ALLOW_DIVIDE_BY_ZERO
4516 scm_num_overflow (s_divide
);
4519 return scm_make_real (1.0 / xx
);
4521 else if (SCM_COMPLEXP (x
))
4523 double r
= SCM_COMPLEX_REAL (x
);
4524 double i
= SCM_COMPLEX_IMAG (x
);
4528 double d
= i
* (1.0 + t
* t
);
4529 return scm_make_complex (t
/ d
, -1.0 / d
);
4534 double d
= r
* (1.0 + t
* t
);
4535 return scm_make_complex (1.0 / d
, -t
/ d
);
4538 else if (SCM_FRACTIONP (x
))
4539 return scm_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4540 SCM_FRACTION_NUMERATOR (x
));
4542 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4545 if (SCM_I_INUMP (x
))
4547 long xx
= SCM_I_INUM (x
);
4548 if (SCM_I_INUMP (y
))
4550 long yy
= SCM_I_INUM (y
);
4553 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4554 scm_num_overflow (s_divide
);
4556 return scm_make_real ((double) xx
/ (double) yy
);
4559 else if (xx
% yy
!= 0)
4562 return scm_make_real ((double) xx
/ (double) yy
);
4563 else return scm_make_ratio (x
, y
);
4568 if (SCM_FIXABLE (z
))
4569 return SCM_I_MAKINUM (z
);
4571 return scm_i_long2big (z
);
4574 else if (SCM_BIGP (y
))
4577 return scm_make_real ((double) xx
/ scm_i_big2dbl (y
));
4578 else return scm_make_ratio (x
, y
);
4580 else if (SCM_REALP (y
))
4582 double yy
= SCM_REAL_VALUE (y
);
4583 #ifndef ALLOW_DIVIDE_BY_ZERO
4585 scm_num_overflow (s_divide
);
4588 return scm_make_real ((double) xx
/ yy
);
4590 else if (SCM_COMPLEXP (y
))
4593 complex_div
: /* y _must_ be a complex number */
4595 double r
= SCM_COMPLEX_REAL (y
);
4596 double i
= SCM_COMPLEX_IMAG (y
);
4600 double d
= i
* (1.0 + t
* t
);
4601 return scm_make_complex ((a
* t
) / d
, -a
/ d
);
4606 double d
= r
* (1.0 + t
* t
);
4607 return scm_make_complex (a
/ d
, -(a
* t
) / d
);
4611 else if (SCM_FRACTIONP (y
))
4612 /* a / b/c = ac / b */
4613 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4614 SCM_FRACTION_NUMERATOR (y
));
4616 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4618 else if (SCM_BIGP (x
))
4620 if (SCM_I_INUMP (y
))
4622 long int yy
= SCM_I_INUM (y
);
4625 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4626 scm_num_overflow (s_divide
);
4628 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4629 scm_remember_upto_here_1 (x
);
4630 return (sgn
== 0) ? scm_nan () : scm_inf ();
4637 /* FIXME: HMM, what are the relative performance issues here?
4638 We need to test. Is it faster on average to test
4639 divisible_p, then perform whichever operation, or is it
4640 faster to perform the integer div opportunistically and
4641 switch to real if there's a remainder? For now we take the
4642 middle ground: test, then if divisible, use the faster div
4645 long abs_yy
= yy
< 0 ? -yy
: yy
;
4646 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4650 SCM result
= scm_i_mkbig ();
4651 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4652 scm_remember_upto_here_1 (x
);
4654 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4655 return scm_i_normbig (result
);
4660 return scm_make_real (scm_i_big2dbl (x
) / (double) yy
);
4661 else return scm_make_ratio (x
, y
);
4665 else if (SCM_BIGP (y
))
4667 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4670 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4671 scm_num_overflow (s_divide
);
4673 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4674 scm_remember_upto_here_1 (x
);
4675 return (sgn
== 0) ? scm_nan () : scm_inf ();
4681 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4685 SCM result
= scm_i_mkbig ();
4686 mpz_divexact (SCM_I_BIG_MPZ (result
),
4689 scm_remember_upto_here_2 (x
, y
);
4690 return scm_i_normbig (result
);
4696 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4697 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4698 scm_remember_upto_here_2 (x
, y
);
4699 return scm_make_real (dbx
/ dby
);
4701 else return scm_make_ratio (x
, y
);
4705 else if (SCM_REALP (y
))
4707 double yy
= SCM_REAL_VALUE (y
);
4708 #ifndef ALLOW_DIVIDE_BY_ZERO
4710 scm_num_overflow (s_divide
);
4713 return scm_make_real (scm_i_big2dbl (x
) / yy
);
4715 else if (SCM_COMPLEXP (y
))
4717 a
= scm_i_big2dbl (x
);
4720 else if (SCM_FRACTIONP (y
))
4721 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4722 SCM_FRACTION_NUMERATOR (y
));
4724 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4726 else if (SCM_REALP (x
))
4728 double rx
= SCM_REAL_VALUE (x
);
4729 if (SCM_I_INUMP (y
))
4731 long int yy
= SCM_I_INUM (y
);
4732 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4734 scm_num_overflow (s_divide
);
4737 return scm_make_real (rx
/ (double) yy
);
4739 else if (SCM_BIGP (y
))
4741 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4742 scm_remember_upto_here_1 (y
);
4743 return scm_make_real (rx
/ dby
);
4745 else if (SCM_REALP (y
))
4747 double yy
= SCM_REAL_VALUE (y
);
4748 #ifndef ALLOW_DIVIDE_BY_ZERO
4750 scm_num_overflow (s_divide
);
4753 return scm_make_real (rx
/ yy
);
4755 else if (SCM_COMPLEXP (y
))
4760 else if (SCM_FRACTIONP (y
))
4761 return scm_make_real (rx
/ scm_i_fraction2double (y
));
4763 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4765 else if (SCM_COMPLEXP (x
))
4767 double rx
= SCM_COMPLEX_REAL (x
);
4768 double ix
= SCM_COMPLEX_IMAG (x
);
4769 if (SCM_I_INUMP (y
))
4771 long int yy
= SCM_I_INUM (y
);
4772 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4774 scm_num_overflow (s_divide
);
4779 return scm_make_complex (rx
/ d
, ix
/ d
);
4782 else if (SCM_BIGP (y
))
4784 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4785 scm_remember_upto_here_1 (y
);
4786 return scm_make_complex (rx
/ dby
, ix
/ dby
);
4788 else if (SCM_REALP (y
))
4790 double yy
= SCM_REAL_VALUE (y
);
4791 #ifndef ALLOW_DIVIDE_BY_ZERO
4793 scm_num_overflow (s_divide
);
4796 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4798 else if (SCM_COMPLEXP (y
))
4800 double ry
= SCM_COMPLEX_REAL (y
);
4801 double iy
= SCM_COMPLEX_IMAG (y
);
4805 double d
= iy
* (1.0 + t
* t
);
4806 return scm_make_complex ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4811 double d
= ry
* (1.0 + t
* t
);
4812 return scm_make_complex ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4815 else if (SCM_FRACTIONP (y
))
4817 double yy
= scm_i_fraction2double (y
);
4818 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4821 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4823 else if (SCM_FRACTIONP (x
))
4825 if (SCM_I_INUMP (y
))
4827 long int yy
= SCM_I_INUM (y
);
4828 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4830 scm_num_overflow (s_divide
);
4833 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4834 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4836 else if (SCM_BIGP (y
))
4838 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4839 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4841 else if (SCM_REALP (y
))
4843 double yy
= SCM_REAL_VALUE (y
);
4844 #ifndef ALLOW_DIVIDE_BY_ZERO
4846 scm_num_overflow (s_divide
);
4849 return scm_make_real (scm_i_fraction2double (x
) / yy
);
4851 else if (SCM_COMPLEXP (y
))
4853 a
= scm_i_fraction2double (x
);
4856 else if (SCM_FRACTIONP (y
))
4857 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4858 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4860 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4863 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4867 scm_divide (SCM x
, SCM y
)
4869 return scm_i_divide (x
, y
, 0);
4872 static SCM
scm_divide2real (SCM x
, SCM y
)
4874 return scm_i_divide (x
, y
, 1);
4880 scm_asinh (double x
)
4885 #define asinh scm_asinh
4886 return log (x
+ sqrt (x
* x
+ 1));
4889 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4890 /* "Return the inverse hyperbolic sine of @var{x}."
4895 scm_acosh (double x
)
4900 #define acosh scm_acosh
4901 return log (x
+ sqrt (x
* x
- 1));
4904 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4905 /* "Return the inverse hyperbolic cosine of @var{x}."
4910 scm_atanh (double x
)
4915 #define atanh scm_atanh
4916 return 0.5 * log ((1 + x
) / (1 - x
));
4919 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4920 /* "Return the inverse hyperbolic tangent of @var{x}."
4924 /* XXX - eventually, we should remove this definition of scm_round and
4925 rename scm_round_number to scm_round. Likewise for scm_truncate
4926 and scm_truncate_number.
4930 scm_truncate (double x
)
4941 /* scm_round is done using floor(x+0.5) to round to nearest and with
4942 half-way case (ie. when x is an integer plus 0.5) going upwards. Then
4943 half-way cases are identified and adjusted down if the round-upwards
4944 didn't give the desired even integer.
4946 "plus_half == result" identifies a half-way case. If plus_half, which is
4947 x + 0.5, is an integer then x must be an integer plus 0.5.
4949 An odd "result" value is identified with result/2 != floor(result/2).
4950 This is done with plus_half, since that value is ready for use sooner in
4951 a pipelined cpu, and we're already requiring plus_half == result.
4953 Note however that we need to be careful when x is big and already an
4954 integer. In that case "x+0.5" may round to an adjacent integer, causing
4955 us to return such a value, incorrectly. For instance if the hardware is
4956 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
4957 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
4958 returned. Or if the hardware is in round-upwards mode, then other bigger
4959 values like say x == 2^128 will see x+0.5 rounding up to the next higher
4960 representable value, 2^128+2^76 (or whatever), again incorrect.
4962 These bad roundings of x+0.5 are avoided by testing at the start whether
4963 x is already an integer. If it is then clearly that's the desired result
4964 already. And if it's not then the exponent must be small enough to allow
4965 an 0.5 to be represented, and hence added without a bad rounding. */
4968 scm_round (double x
)
4970 double plus_half
, result
;
4975 plus_half
= x
+ 0.5;
4976 result
= floor (plus_half
);
4977 /* Adjust so that the scm_round is towards even. */
4978 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4983 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4985 "Round the number @var{x} towards zero.")
4986 #define FUNC_NAME s_scm_truncate_number
4988 if (scm_is_false (scm_negative_p (x
)))
4989 return scm_floor (x
);
4991 return scm_ceiling (x
);
4995 static SCM exactly_one_half
;
4997 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4999 "Round the number @var{x} towards the nearest integer. "
5000 "When it is exactly halfway between two integers, "
5001 "round towards the even one.")
5002 #define FUNC_NAME s_scm_round_number
5004 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5006 else if (SCM_REALP (x
))
5007 return scm_make_real (scm_round (SCM_REAL_VALUE (x
)));
5010 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5011 single quotient+remainder division then examining to see which way
5012 the rounding should go. */
5013 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5014 SCM result
= scm_floor (plus_half
);
5015 /* Adjust so that the scm_round is towards even. */
5016 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5017 && scm_is_true (scm_odd_p (result
)))
5018 return scm_difference (result
, SCM_I_MAKINUM (1));
5025 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5027 "Round the number @var{x} towards minus infinity.")
5028 #define FUNC_NAME s_scm_floor
5030 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5032 else if (SCM_REALP (x
))
5033 return scm_make_real (floor (SCM_REAL_VALUE (x
)));
5034 else if (SCM_FRACTIONP (x
))
5036 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5037 SCM_FRACTION_DENOMINATOR (x
));
5038 if (scm_is_false (scm_negative_p (x
)))
5040 /* For positive x, rounding towards zero is correct. */
5045 /* For negative x, we need to return q-1 unless x is an
5046 integer. But fractions are never integer, per our
5048 return scm_difference (q
, SCM_I_MAKINUM (1));
5052 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5056 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5058 "Round the number @var{x} towards infinity.")
5059 #define FUNC_NAME s_scm_ceiling
5061 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5063 else if (SCM_REALP (x
))
5064 return scm_make_real (ceil (SCM_REAL_VALUE (x
)));
5065 else if (SCM_FRACTIONP (x
))
5067 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5068 SCM_FRACTION_DENOMINATOR (x
));
5069 if (scm_is_false (scm_positive_p (x
)))
5071 /* For negative x, rounding towards zero is correct. */
5076 /* For positive x, we need to return q+1 unless x is an
5077 integer. But fractions are never integer, per our
5079 return scm_sum (q
, SCM_I_MAKINUM (1));
5083 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5087 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5088 /* "Return the square root of the real number @var{x}."
5090 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5091 /* "Return the absolute value of the real number @var{x}."
5093 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5094 /* "Return the @var{x}th power of e."
5096 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5097 /* "Return the natural logarithm of the real number @var{x}."
5099 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5100 /* "Return the sine of the real number @var{x}."
5102 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5103 /* "Return the cosine of the real number @var{x}."
5105 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5106 /* "Return the tangent of the real number @var{x}."
5108 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5109 /* "Return the arc sine of the real number @var{x}."
5111 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5112 /* "Return the arc cosine of the real number @var{x}."
5114 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5115 /* "Return the arc tangent of the real number @var{x}."
5117 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5118 /* "Return the hyperbolic sine of the real number @var{x}."
5120 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5121 /* "Return the hyperbolic cosine of the real number @var{x}."
5123 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5124 /* "Return the hyperbolic tangent of the real number @var{x}."
5132 static void scm_two_doubles (SCM x
,
5134 const char *sstring
,
5138 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5140 if (SCM_I_INUMP (x
))
5141 xy
->x
= SCM_I_INUM (x
);
5142 else if (SCM_BIGP (x
))
5143 xy
->x
= scm_i_big2dbl (x
);
5144 else if (SCM_REALP (x
))
5145 xy
->x
= SCM_REAL_VALUE (x
);
5146 else if (SCM_FRACTIONP (x
))
5147 xy
->x
= scm_i_fraction2double (x
);
5149 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5151 if (SCM_I_INUMP (y
))
5152 xy
->y
= SCM_I_INUM (y
);
5153 else if (SCM_BIGP (y
))
5154 xy
->y
= scm_i_big2dbl (y
);
5155 else if (SCM_REALP (y
))
5156 xy
->y
= SCM_REAL_VALUE (y
);
5157 else if (SCM_FRACTIONP (y
))
5158 xy
->y
= scm_i_fraction2double (y
);
5160 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5164 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5166 "Return @var{x} raised to the power of @var{y}. This\n"
5167 "procedure does not accept complex arguments.")
5168 #define FUNC_NAME s_scm_sys_expt
5171 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5172 return scm_make_real (pow (xy
.x
, xy
.y
));
5177 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5179 "Return the arc tangent of the two arguments @var{x} and\n"
5180 "@var{y}. This is similar to calculating the arc tangent of\n"
5181 "@var{x} / @var{y}, except that the signs of both arguments\n"
5182 "are used to determine the quadrant of the result. This\n"
5183 "procedure does not accept complex arguments.")
5184 #define FUNC_NAME s_scm_sys_atan2
5187 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5188 return scm_make_real (atan2 (xy
.x
, xy
.y
));
5193 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5194 (SCM real
, SCM imaginary
),
5195 "Return a complex number constructed of the given @var{real} and\n"
5196 "@var{imaginary} parts.")
5197 #define FUNC_NAME s_scm_make_rectangular
5200 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5201 return scm_make_complex (xy
.x
, xy
.y
);
5207 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5209 "Return the complex number @var{x} * e^(i * @var{y}).")
5210 #define FUNC_NAME s_scm_make_polar
5214 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5216 sincos (xy
.y
, &s
, &c
);
5221 return scm_make_complex (xy
.x
* c
, xy
.x
* s
);
5226 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5227 /* "Return the real part of the number @var{z}."
5230 scm_real_part (SCM z
)
5232 if (SCM_I_INUMP (z
))
5234 else if (SCM_BIGP (z
))
5236 else if (SCM_REALP (z
))
5238 else if (SCM_COMPLEXP (z
))
5239 return scm_make_real (SCM_COMPLEX_REAL (z
));
5240 else if (SCM_FRACTIONP (z
))
5243 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5247 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5248 /* "Return the imaginary part of the number @var{z}."
5251 scm_imag_part (SCM z
)
5253 if (SCM_I_INUMP (z
))
5255 else if (SCM_BIGP (z
))
5257 else if (SCM_REALP (z
))
5259 else if (SCM_COMPLEXP (z
))
5260 return scm_make_real (SCM_COMPLEX_IMAG (z
));
5261 else if (SCM_FRACTIONP (z
))
5264 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5267 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5268 /* "Return the numerator of the number @var{z}."
5271 scm_numerator (SCM z
)
5273 if (SCM_I_INUMP (z
))
5275 else if (SCM_BIGP (z
))
5277 else if (SCM_FRACTIONP (z
))
5279 scm_i_fraction_reduce (z
);
5280 return SCM_FRACTION_NUMERATOR (z
);
5282 else if (SCM_REALP (z
))
5283 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5285 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5289 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5290 /* "Return the denominator of the number @var{z}."
5293 scm_denominator (SCM z
)
5295 if (SCM_I_INUMP (z
))
5296 return SCM_I_MAKINUM (1);
5297 else if (SCM_BIGP (z
))
5298 return SCM_I_MAKINUM (1);
5299 else if (SCM_FRACTIONP (z
))
5301 scm_i_fraction_reduce (z
);
5302 return SCM_FRACTION_DENOMINATOR (z
);
5304 else if (SCM_REALP (z
))
5305 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5307 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5310 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5311 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5312 * "@code{abs} for real arguments, but also allows complex numbers."
5315 scm_magnitude (SCM z
)
5317 if (SCM_I_INUMP (z
))
5319 long int zz
= SCM_I_INUM (z
);
5322 else if (SCM_POSFIXABLE (-zz
))
5323 return SCM_I_MAKINUM (-zz
);
5325 return scm_i_long2big (-zz
);
5327 else if (SCM_BIGP (z
))
5329 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5330 scm_remember_upto_here_1 (z
);
5332 return scm_i_clonebig (z
, 0);
5336 else if (SCM_REALP (z
))
5337 return scm_make_real (fabs (SCM_REAL_VALUE (z
)));
5338 else if (SCM_COMPLEXP (z
))
5339 return scm_make_real (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5340 else if (SCM_FRACTIONP (z
))
5342 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5344 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5345 SCM_FRACTION_DENOMINATOR (z
));
5348 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5352 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5353 /* "Return the angle of the complex number @var{z}."
5358 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5359 scm_flo0 to save allocating a new flonum with scm_make_real each time.
5360 But if atan2 follows the floating point rounding mode, then the value
5361 is not a constant. Maybe it'd be close enough though. */
5362 if (SCM_I_INUMP (z
))
5364 if (SCM_I_INUM (z
) >= 0)
5367 return scm_make_real (atan2 (0.0, -1.0));
5369 else if (SCM_BIGP (z
))
5371 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5372 scm_remember_upto_here_1 (z
);
5374 return scm_make_real (atan2 (0.0, -1.0));
5378 else if (SCM_REALP (z
))
5380 if (SCM_REAL_VALUE (z
) >= 0)
5383 return scm_make_real (atan2 (0.0, -1.0));
5385 else if (SCM_COMPLEXP (z
))
5386 return scm_make_real (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5387 else if (SCM_FRACTIONP (z
))
5389 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5391 else return scm_make_real (atan2 (0.0, -1.0));
5394 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5398 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5399 /* Convert the number @var{x} to its inexact representation.\n"
5402 scm_exact_to_inexact (SCM z
)
5404 if (SCM_I_INUMP (z
))
5405 return scm_make_real ((double) SCM_I_INUM (z
));
5406 else if (SCM_BIGP (z
))
5407 return scm_make_real (scm_i_big2dbl (z
));
5408 else if (SCM_FRACTIONP (z
))
5409 return scm_make_real (scm_i_fraction2double (z
));
5410 else if (SCM_INEXACTP (z
))
5413 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5417 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5419 "Return an exact number that is numerically closest to @var{z}.")
5420 #define FUNC_NAME s_scm_inexact_to_exact
5422 if (SCM_I_INUMP (z
))
5424 else if (SCM_BIGP (z
))
5426 else if (SCM_REALP (z
))
5428 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5429 SCM_OUT_OF_RANGE (1, z
);
5436 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5437 q
= scm_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5438 scm_i_mpz2num (mpq_denref (frac
)));
5440 /* When scm_make_ratio throws, we leak the memory allocated
5447 else if (SCM_FRACTIONP (z
))
5450 SCM_WRONG_TYPE_ARG (1, z
);
5454 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5456 "Return an exact number that is within @var{err} of @var{x}.")
5457 #define FUNC_NAME s_scm_rationalize
5459 if (SCM_I_INUMP (x
))
5461 else if (SCM_BIGP (x
))
5463 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5465 /* Use continued fractions to find closest ratio. All
5466 arithmetic is done with exact numbers.
5469 SCM ex
= scm_inexact_to_exact (x
);
5470 SCM int_part
= scm_floor (ex
);
5471 SCM tt
= SCM_I_MAKINUM (1);
5472 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5473 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5477 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5480 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5481 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5483 /* We stop after a million iterations just to be absolutely sure
5484 that we don't go into an infinite loop. The process normally
5485 converges after less than a dozen iterations.
5488 err
= scm_abs (err
);
5489 while (++i
< 1000000)
5491 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5492 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5493 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5495 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5496 err
))) /* abs(x-a/b) <= err */
5498 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5499 if (scm_is_false (scm_exact_p (x
))
5500 || scm_is_false (scm_exact_p (err
)))
5501 return scm_exact_to_inexact (res
);
5505 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5507 tt
= scm_floor (rx
); /* tt = floor (rx) */
5513 scm_num_overflow (s_scm_rationalize
);
5516 SCM_WRONG_TYPE_ARG (1, x
);
5520 /* Parameters for creating integer conversion routines.
5522 Define the following preprocessor macros before including
5523 "libguile/num2integral.i.c":
5525 NUM2INTEGRAL - the name of the function for converting from a
5526 Scheme object to the integral type. This function will be
5527 defined when including "num2integral.i.c".
5529 INTEGRAL2NUM - the name of the function for converting from the
5530 integral type to a Scheme object. This function will be defined.
5532 INTEGRAL2BIG - the name of an internal function that createas a
5533 bignum from the integral type. This function will be defined.
5534 The name should start with "scm_i_".
5536 ITYPE - the name of the integral type.
5538 UNSIGNED - Define this to 1 when ITYPE is an unsigned type. Define
5541 UNSIGNED_ITYPE - the name of the the unsigned variant of the
5542 integral type. If you don't define this, it defaults to
5543 "unsigned ITYPE" for signed types and simply "ITYPE" for unsigned
5546 SIZEOF_ITYPE - an expression giving the size of the integral type
5547 in bytes. This expression must be computable by the
5548 preprocessor. (SIZEOF_FOO values are calculated by configure.in
5553 #define NUM2INTEGRAL scm_num2short
5554 #define INTEGRAL2NUM scm_short2num
5555 #define INTEGRAL2BIG scm_i_short2big
5558 #define SIZEOF_ITYPE SIZEOF_SHORT
5559 #include "libguile/num2integral.i.c"
5561 #define NUM2INTEGRAL scm_num2ushort
5562 #define INTEGRAL2NUM scm_ushort2num
5563 #define INTEGRAL2BIG scm_i_ushort2big
5565 #define ITYPE unsigned short
5566 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_SHORT
5567 #include "libguile/num2integral.i.c"
5569 #define NUM2INTEGRAL scm_num2int
5570 #define INTEGRAL2NUM scm_int2num
5571 #define INTEGRAL2BIG scm_i_int2big
5574 #define SIZEOF_ITYPE SIZEOF_INT
5575 #include "libguile/num2integral.i.c"
5577 #define NUM2INTEGRAL scm_num2uint
5578 #define INTEGRAL2NUM scm_uint2num
5579 #define INTEGRAL2BIG scm_i_uint2big
5581 #define ITYPE unsigned int
5582 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_INT
5583 #include "libguile/num2integral.i.c"
5585 #define NUM2INTEGRAL scm_num2long
5586 #define INTEGRAL2NUM scm_long2num
5587 #define INTEGRAL2BIG scm_i_long2big
5590 #define SIZEOF_ITYPE SIZEOF_LONG
5591 #include "libguile/num2integral.i.c"
5593 #define NUM2INTEGRAL scm_num2ulong
5594 #define INTEGRAL2NUM scm_ulong2num
5595 #define INTEGRAL2BIG scm_i_ulong2big
5597 #define ITYPE unsigned long
5598 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG
5599 #include "libguile/num2integral.i.c"
5601 #define NUM2INTEGRAL scm_num2ptrdiff
5602 #define INTEGRAL2NUM scm_ptrdiff2num
5603 #define INTEGRAL2BIG scm_i_ptrdiff2big
5605 #define ITYPE scm_t_ptrdiff
5606 #define UNSIGNED_ITYPE size_t
5607 #define SIZEOF_ITYPE SCM_SIZEOF_SCM_T_PTRDIFF
5608 #include "libguile/num2integral.i.c"
5610 #define NUM2INTEGRAL scm_num2size
5611 #define INTEGRAL2NUM scm_size2num
5612 #define INTEGRAL2BIG scm_i_size2big
5614 #define ITYPE size_t
5615 #define SIZEOF_ITYPE SIZEOF_SIZE_T
5616 #include "libguile/num2integral.i.c"
5618 #if SCM_SIZEOF_LONG_LONG != 0
5620 #define NUM2INTEGRAL scm_num2long_long
5621 #define INTEGRAL2NUM scm_long_long2num
5622 #define INTEGRAL2BIG scm_i_long_long2big
5624 #define ITYPE long long
5625 #define SIZEOF_ITYPE SIZEOF_LONG_LONG
5626 #include "libguile/num2integral.i.c"
5628 #define NUM2INTEGRAL scm_num2ulong_long
5629 #define INTEGRAL2NUM scm_ulong_long2num
5630 #define INTEGRAL2BIG scm_i_ulong_long2big
5632 #define ITYPE unsigned long long
5633 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG_LONG
5634 #include "libguile/num2integral.i.c"
5636 #endif /* SCM_SIZEOF_LONG_LONG != 0 */
5638 #define NUM2FLOAT scm_num2float
5639 #define FLOAT2NUM scm_float2num
5641 #include "libguile/num2float.i.c"
5643 #define NUM2FLOAT scm_num2double
5644 #define FLOAT2NUM scm_double2num
5645 #define FTYPE double
5646 #include "libguile/num2float.i.c"
5648 /* conversion functions */
5651 scm_is_integer (SCM val
)
5653 return scm_is_true (scm_integer_p (val
));
5657 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5659 if (SCM_I_INUMP (val
))
5661 scm_t_signed_bits n
= SCM_I_INUM (val
);
5662 return n
>= min
&& n
<= max
;
5664 else if (SCM_BIGP (val
))
5666 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5668 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5670 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5672 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5673 return n
>= min
&& n
<= max
;
5683 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5684 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5687 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5688 SCM_I_BIG_MPZ (val
));
5690 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5702 return n
>= min
&& n
<= max
;
5710 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5712 if (SCM_I_INUMP (val
))
5714 scm_t_signed_bits n
= SCM_I_INUM (val
);
5715 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5717 else if (SCM_BIGP (val
))
5719 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5721 else if (max
<= ULONG_MAX
)
5723 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5725 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5726 return n
>= min
&& n
<= max
;
5736 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5739 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5740 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5743 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5744 SCM_I_BIG_MPZ (val
));
5746 return n
>= min
&& n
<= max
;
5754 scm_to_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5756 if (SCM_I_INUMP (val
))
5758 scm_t_signed_bits n
= SCM_I_INUM (val
);
5759 if (n
>= min
&& n
<= max
)
5764 scm_out_of_range (NULL
, val
);
5768 else if (SCM_BIGP (val
))
5770 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5772 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5774 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5776 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5777 if (n
>= min
&& n
<= max
)
5790 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5791 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5794 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5795 SCM_I_BIG_MPZ (val
));
5797 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5809 if (n
>= min
&& n
<= max
)
5817 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5823 scm_to_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5825 if (SCM_I_INUMP (val
))
5827 scm_t_signed_bits n
= SCM_I_INUM (val
);
5828 if (n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
)
5833 scm_out_of_range (NULL
, val
);
5837 else if (SCM_BIGP (val
))
5839 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5841 else if (max
<= ULONG_MAX
)
5843 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5845 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5846 if (n
>= min
&& n
<= max
)
5859 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5862 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5863 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5866 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5867 SCM_I_BIG_MPZ (val
));
5869 if (n
>= min
&& n
<= max
)
5877 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5883 scm_from_signed_integer (scm_t_intmax val
)
5885 if (SCM_FIXABLE (val
))
5886 return SCM_I_MAKINUM (val
);
5887 else if (val
>= LONG_MIN
&& val
<= LONG_MAX
)
5889 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
5890 mpz_init_set_si (SCM_I_BIG_MPZ (z
), val
);
5895 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
5896 mpz_init (SCM_I_BIG_MPZ (z
));
5900 mpz_import (SCM_I_BIG_MPZ (z
), 1, 1, sizeof (scm_t_intmax
), 0, 0,
5902 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
5905 mpz_import (SCM_I_BIG_MPZ (z
), 1, 1, sizeof (scm_t_intmax
), 0, 0,
5912 scm_from_unsigned_integer (scm_t_uintmax val
)
5914 if (SCM_POSFIXABLE (val
))
5915 return SCM_I_MAKINUM (val
);
5916 else if (val
<= ULONG_MAX
)
5918 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
5919 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), val
);
5924 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
5925 mpz_init (SCM_I_BIG_MPZ (z
));
5926 mpz_import (SCM_I_BIG_MPZ (z
), 1, 1, sizeof (scm_t_uintmax
), 0, 0,
5933 scm_is_real (SCM val
)
5935 return scm_is_true (scm_real_p (val
));
5939 scm_to_double (SCM val
)
5941 return scm_num2dbl (val
, NULL
);
5945 scm_from_double (double val
)
5947 return scm_make_real (val
);
5953 #define SIZE_MAX ((size_t) (-1))
5956 #define PTRDIFF_MIN \
5957 ((scm_t_ptrdiff) ((scm_t_ptrdiff) 1 \
5958 << ((sizeof (scm_t_ptrdiff) * SCM_CHAR_BIT) - 1)))
5961 #define PTRDIFF_MAX (~ PTRDIFF_MIN)
5964 #define CHECK(type, v) \
5967 if ((v) != scm_num2##type (scm_##type##2num (v), 1, "check_sanity")) \
5987 CHECK (ptrdiff
, -1);
5989 CHECK (short, SHRT_MAX
);
5990 CHECK (short, SHRT_MIN
);
5991 CHECK (ushort
, USHRT_MAX
);
5992 CHECK (int, INT_MAX
);
5993 CHECK (int, INT_MIN
);
5994 CHECK (uint
, UINT_MAX
);
5995 CHECK (long, LONG_MAX
);
5996 CHECK (long, LONG_MIN
);
5997 CHECK (ulong
, ULONG_MAX
);
5998 CHECK (size
, SIZE_MAX
);
5999 CHECK (ptrdiff
, PTRDIFF_MAX
);
6000 CHECK (ptrdiff
, PTRDIFF_MIN
);
6002 #if SCM_SIZEOF_LONG_LONG != 0
6003 CHECK (long_long
, 0LL);
6004 CHECK (ulong_long
, 0ULL);
6005 CHECK (long_long
, -1LL);
6006 CHECK (long_long
, SCM_I_LLONG_MAX
);
6007 CHECK (long_long
, SCM_I_LLONG_MIN
);
6008 CHECK (ulong_long
, SCM_I_ULLONG_MAX
);
6015 scm_internal_catch (SCM_BOOL_T, check_body, &data, check_handler, &data); \
6016 if (scm_is_true (data)) abort();
6019 check_body (void *data
)
6021 SCM num
= *(SCM
*) data
;
6022 scm_num2ulong (num
, 1, NULL
);
6024 return SCM_UNSPECIFIED
;
6028 check_handler (void *data
, SCM tag
, SCM throw_args
)
6030 SCM
*num
= (SCM
*) data
;
6033 return SCM_UNSPECIFIED
;
6036 SCM_DEFINE (scm_sys_check_number_conversions
, "%check-number-conversions", 0, 0, 0,
6038 "Number conversion sanity checking.")
6039 #define FUNC_NAME s_scm_sys_check_number_conversions
6041 SCM data
= SCM_I_MAKINUM (-1);
6043 data
= scm_int2num (INT_MIN
);
6045 data
= scm_ulong2num (ULONG_MAX
);
6046 data
= scm_difference (SCM_INUM0
, data
);
6048 data
= scm_ulong2num (ULONG_MAX
);
6049 data
= scm_sum (SCM_I_MAKINUM (1), data
); data
= scm_difference (SCM_INUM0
, data
);
6051 data
= scm_int2num (-10000); data
= scm_product (data
, data
); data
= scm_product (data
, data
);
6054 return SCM_UNSPECIFIED
;
6065 mpz_init_set_si (z_negative_one
, -1);
6067 /* It may be possible to tune the performance of some algorithms by using
6068 * the following constants to avoid the creation of bignums. Please, before
6069 * using these values, remember the two rules of program optimization:
6070 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6071 scm_c_define ("most-positive-fixnum",
6072 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6073 scm_c_define ("most-negative-fixnum",
6074 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6076 scm_add_feature ("complex");
6077 scm_add_feature ("inexact");
6078 scm_flo0
= scm_make_real (0.0);
6080 /* determine floating point precision */
6081 for(i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6083 init_dblprec(&scm_dblprec
[i
-2],i
);
6084 init_fx_radix(fx_per_radix
[i
-2],i
);
6087 /* hard code precision for base 10 if the preprocessor tells us to... */
6088 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6095 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6096 SCM_I_MAKINUM (2)));
6097 #include "libguile/numbers.x"