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 SCM
168 scm_i_long2big (long x
)
170 /* Return a newly created bignum initialized to X. */
171 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
172 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
176 SCM_C_INLINE_KEYWORD SCM
177 scm_i_ulong2big (unsigned long x
)
179 /* Return a newly created bignum initialized to X. */
180 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
181 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
185 SCM_C_INLINE_KEYWORD
static SCM
186 scm_i_clonebig (SCM src_big
, int same_sign_p
)
188 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
189 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
190 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
192 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
196 SCM_C_INLINE_KEYWORD
int
197 scm_i_bigcmp (SCM x
, SCM y
)
199 /* Return neg if x < y, pos if x > y, and 0 if x == y */
200 /* presume we already know x and y are bignums */
201 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
202 scm_remember_upto_here_2 (x
, y
);
206 SCM_C_INLINE_KEYWORD SCM
207 scm_i_dbl2big (double d
)
209 /* results are only defined if d is an integer */
210 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
211 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
215 /* Convert a integer in double representation to a SCM number. */
217 SCM_C_INLINE_KEYWORD SCM
218 scm_i_dbl2num (double u
)
220 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
221 powers of 2, so there's no rounding when making "double" values
222 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
223 get rounded on a 64-bit machine, hence the "+1".
225 The use of floor() to force to an integer value ensures we get a
226 "numerically closest" value without depending on how a
227 double->long cast or how mpz_set_d will round. For reference,
228 double->long probably follows the hardware rounding mode,
229 mpz_set_d truncates towards zero. */
231 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
232 representable as a double? */
234 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
235 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
236 return SCM_I_MAKINUM ((long) u
);
238 return scm_i_dbl2big (u
);
241 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
242 with R5RS exact->inexact.
244 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
245 (ie. it truncates towards zero), then adjust to get the closest double by
246 examining the next lower bit and adding 1 if necessary.
248 Note that bignums exactly half way between representable doubles are
249 rounded to the next higher absolute value (ie. away from zero). This
250 seems like an adequate interpretation of R5RS "numerically closest", and
251 it's easier and faster than a full "nearest-even" style.
253 The bit test is done on the absolute value of the mpz_t, which means we
254 must use mpz_getlimbn. mpz_tstbit is not right, it treats negatives as
257 Prior to GMP 4.2, the rounding done by mpz_get_d was unspecified. It
258 happened to follow the hardware rounding mode, but on the absolute value
259 of its operand. This is not what we want, so we put the high
260 DBL_MANT_DIG bits into a temporary. This extra init/clear is a slowdown,
261 but doesn't matter too much since it's only for older GMP. */
264 scm_i_big2dbl (SCM b
)
269 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
271 #if __GNU_MP_VERSION < 4 \
272 || (__GNU_MP_VERSION == 4 && __GNU_MP_VERSION_MINOR < 2)
274 /* GMP prior to 4.2, force truncate towards zero */
276 if (bits
> DBL_MANT_DIG
)
278 size_t shift
= bits
- DBL_MANT_DIG
;
279 mpz_init2 (tmp
, DBL_MANT_DIG
);
280 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
281 result
= ldexp (mpz_get_d (tmp
), shift
);
286 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
291 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
294 if (bits
> DBL_MANT_DIG
)
296 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
297 /* test bit number "pos" in absolute value */
298 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
299 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
301 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
305 scm_remember_upto_here_1 (b
);
309 SCM_C_INLINE_KEYWORD SCM
310 scm_i_normbig (SCM b
)
312 /* convert a big back to a fixnum if it'll fit */
313 /* presume b is a bignum */
314 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
316 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
317 if (SCM_FIXABLE (val
))
318 b
= SCM_I_MAKINUM (val
);
323 static SCM_C_INLINE_KEYWORD SCM
324 scm_i_mpz2num (mpz_t b
)
326 /* convert a mpz number to a SCM number. */
327 if (mpz_fits_slong_p (b
))
329 long val
= mpz_get_si (b
);
330 if (SCM_FIXABLE (val
))
331 return SCM_I_MAKINUM (val
);
335 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
336 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
341 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
342 static SCM
scm_divide2real (SCM x
, SCM y
);
345 scm_make_ratio (SCM numerator
, SCM denominator
)
346 #define FUNC_NAME "make-ratio"
348 /* First make sure the arguments are proper.
350 if (SCM_I_INUMP (denominator
))
352 if (scm_is_eq (denominator
, SCM_INUM0
))
353 scm_num_overflow ("make-ratio");
354 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
359 if (!(SCM_BIGP(denominator
)))
360 SCM_WRONG_TYPE_ARG (2, denominator
);
362 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
363 SCM_WRONG_TYPE_ARG (1, numerator
);
365 /* Then flip signs so that the denominator is positive.
367 if (scm_is_true (scm_negative_p (denominator
)))
369 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
370 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
373 /* Now consider for each of the four fixnum/bignum combinations
374 whether the rational number is really an integer.
376 if (SCM_I_INUMP (numerator
))
378 long x
= SCM_I_INUM (numerator
);
379 if (scm_is_eq (numerator
, SCM_INUM0
))
381 if (SCM_I_INUMP (denominator
))
384 y
= SCM_I_INUM (denominator
);
386 return SCM_I_MAKINUM(1);
388 return SCM_I_MAKINUM (x
/ y
);
392 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
393 of that value for the denominator, as a bignum. Apart from
394 that case, abs(bignum) > abs(inum) so inum/bignum is not an
396 if (x
== SCM_MOST_NEGATIVE_FIXNUM
397 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
398 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
399 return SCM_I_MAKINUM(-1);
402 else if (SCM_BIGP (numerator
))
404 if (SCM_I_INUMP (denominator
))
406 long yy
= SCM_I_INUM (denominator
);
407 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
408 return scm_divide (numerator
, denominator
);
412 if (scm_is_eq (numerator
, denominator
))
413 return SCM_I_MAKINUM(1);
414 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
415 SCM_I_BIG_MPZ (denominator
)))
416 return scm_divide(numerator
, denominator
);
420 /* No, it's a proper fraction.
422 return scm_double_cell (scm_tc16_fraction
,
423 SCM_UNPACK (numerator
),
424 SCM_UNPACK (denominator
), 0);
428 static void scm_i_fraction_reduce (SCM z
)
430 if (!(SCM_FRACTION_REDUCED (z
)))
433 divisor
= scm_gcd (SCM_FRACTION_NUMERATOR (z
), SCM_FRACTION_DENOMINATOR (z
));
434 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
437 SCM_FRACTION_SET_NUMERATOR (z
, scm_divide (SCM_FRACTION_NUMERATOR (z
), divisor
));
438 SCM_FRACTION_SET_DENOMINATOR (z
, scm_divide (SCM_FRACTION_DENOMINATOR (z
), divisor
));
440 SCM_FRACTION_REDUCED_SET (z
);
445 scm_i_fraction2double (SCM z
)
447 return scm_num2dbl (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
448 SCM_FRACTION_DENOMINATOR (z
)),
452 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
454 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
456 #define FUNC_NAME s_scm_exact_p
462 if (SCM_FRACTIONP (x
))
466 SCM_WRONG_TYPE_ARG (1, x
);
471 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
473 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
475 #define FUNC_NAME s_scm_odd_p
479 long val
= SCM_I_INUM (n
);
480 return scm_from_bool ((val
& 1L) != 0);
482 else if (SCM_BIGP (n
))
484 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
485 scm_remember_upto_here_1 (n
);
486 return scm_from_bool (odd_p
);
488 else if (scm_is_true (scm_inf_p (n
)))
490 else if (SCM_REALP (n
))
492 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
498 SCM_WRONG_TYPE_ARG (1, n
);
501 SCM_WRONG_TYPE_ARG (1, n
);
506 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
508 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
510 #define FUNC_NAME s_scm_even_p
514 long val
= SCM_I_INUM (n
);
515 return scm_from_bool ((val
& 1L) == 0);
517 else if (SCM_BIGP (n
))
519 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
520 scm_remember_upto_here_1 (n
);
521 return scm_from_bool (even_p
);
523 else if (scm_is_true (scm_inf_p (n
)))
525 else if (SCM_REALP (n
))
527 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
533 SCM_WRONG_TYPE_ARG (1, n
);
536 SCM_WRONG_TYPE_ARG (1, n
);
540 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
542 "Return @code{#t} if @var{n} is infinite, @code{#f}\n"
544 #define FUNC_NAME s_scm_inf_p
547 return scm_from_bool (xisinf (SCM_REAL_VALUE (n
)));
548 else if (SCM_COMPLEXP (n
))
549 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (n
))
550 || xisinf (SCM_COMPLEX_IMAG (n
)));
556 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
558 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
560 #define FUNC_NAME s_scm_nan_p
563 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
564 else if (SCM_COMPLEXP (n
))
565 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
566 || xisnan (SCM_COMPLEX_IMAG (n
)));
572 /* Guile's idea of infinity. */
573 static double guile_Inf
;
575 /* Guile's idea of not a number. */
576 static double guile_NaN
;
579 guile_ieee_init (void)
581 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
583 /* Some version of gcc on some old version of Linux used to crash when
584 trying to make Inf and NaN. */
587 /* C99 INFINITY, when available.
588 FIXME: The standard allows for INFINITY to be something that overflows
589 at compile time. We ought to have a configure test to check for that
590 before trying to use it. (But in practice we believe this is not a
591 problem on any system guile is likely to target.) */
592 guile_Inf
= INFINITY
;
595 extern unsigned int DINFINITY
[2];
596 guile_Inf
= (*(X_CAST(double *, DINFINITY
)));
603 if (guile_Inf
== tmp
)
611 #if defined (HAVE_ISNAN)
614 /* C99 NAN, when available */
618 extern unsigned int DQNAN
[2];
619 guile_NaN
= (*(X_CAST(double *, DQNAN
)));
621 guile_NaN
= guile_Inf
/ guile_Inf
;
627 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
630 #define FUNC_NAME s_scm_inf
632 static int initialized
= 0;
638 return scm_make_real (guile_Inf
);
642 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
645 #define FUNC_NAME s_scm_nan
647 static int initialized
= 0;
653 return scm_make_real (guile_NaN
);
658 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
660 "Return the absolute value of @var{x}.")
665 long int xx
= SCM_I_INUM (x
);
668 else if (SCM_POSFIXABLE (-xx
))
669 return SCM_I_MAKINUM (-xx
);
671 return scm_i_long2big (-xx
);
673 else if (SCM_BIGP (x
))
675 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
677 return scm_i_clonebig (x
, 0);
681 else if (SCM_REALP (x
))
683 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
684 double xx
= SCM_REAL_VALUE (x
);
686 return scm_make_real (-xx
);
690 else if (SCM_FRACTIONP (x
))
692 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
694 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
695 SCM_FRACTION_DENOMINATOR (x
));
698 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
703 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
704 /* "Return the quotient of the numbers @var{x} and @var{y}."
707 scm_quotient (SCM x
, SCM y
)
711 long xx
= SCM_I_INUM (x
);
714 long yy
= SCM_I_INUM (y
);
716 scm_num_overflow (s_quotient
);
721 return SCM_I_MAKINUM (z
);
723 return scm_i_long2big (z
);
726 else if (SCM_BIGP (y
))
728 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
729 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
730 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
732 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
733 scm_remember_upto_here_1 (y
);
734 return SCM_I_MAKINUM (-1);
737 return SCM_I_MAKINUM (0);
740 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
742 else if (SCM_BIGP (x
))
746 long yy
= SCM_I_INUM (y
);
748 scm_num_overflow (s_quotient
);
753 SCM result
= scm_i_mkbig ();
756 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
759 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
762 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
763 scm_remember_upto_here_1 (x
);
764 return scm_i_normbig (result
);
767 else if (SCM_BIGP (y
))
769 SCM result
= scm_i_mkbig ();
770 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
773 scm_remember_upto_here_2 (x
, y
);
774 return scm_i_normbig (result
);
777 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
780 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
783 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
784 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
786 * "(remainder 13 4) @result{} 1\n"
787 * "(remainder -13 4) @result{} -1\n"
791 scm_remainder (SCM x
, SCM y
)
797 long yy
= SCM_I_INUM (y
);
799 scm_num_overflow (s_remainder
);
802 long z
= SCM_I_INUM (x
) % yy
;
803 return SCM_I_MAKINUM (z
);
806 else if (SCM_BIGP (y
))
808 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
809 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
810 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
812 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
813 scm_remember_upto_here_1 (y
);
814 return SCM_I_MAKINUM (0);
820 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
822 else if (SCM_BIGP (x
))
826 long yy
= SCM_I_INUM (y
);
828 scm_num_overflow (s_remainder
);
831 SCM result
= scm_i_mkbig ();
834 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
835 scm_remember_upto_here_1 (x
);
836 return scm_i_normbig (result
);
839 else if (SCM_BIGP (y
))
841 SCM result
= scm_i_mkbig ();
842 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
845 scm_remember_upto_here_2 (x
, y
);
846 return scm_i_normbig (result
);
849 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
852 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
856 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
857 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
859 * "(modulo 13 4) @result{} 1\n"
860 * "(modulo -13 4) @result{} 3\n"
864 scm_modulo (SCM x
, SCM y
)
868 long xx
= SCM_I_INUM (x
);
871 long yy
= SCM_I_INUM (y
);
873 scm_num_overflow (s_modulo
);
876 /* FIXME: I think this may be a bug on some arches -- results
877 of % with negative second arg are undefined... */
895 return SCM_I_MAKINUM (result
);
898 else if (SCM_BIGP (y
))
900 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
907 SCM pos_y
= scm_i_clonebig (y
, 0);
908 /* do this after the last scm_op */
909 mpz_init_set_si (z_x
, xx
);
910 result
= pos_y
; /* re-use this bignum */
911 mpz_mod (SCM_I_BIG_MPZ (result
),
913 SCM_I_BIG_MPZ (pos_y
));
914 scm_remember_upto_here_1 (pos_y
);
918 result
= scm_i_mkbig ();
919 /* do this after the last scm_op */
920 mpz_init_set_si (z_x
, xx
);
921 mpz_mod (SCM_I_BIG_MPZ (result
),
924 scm_remember_upto_here_1 (y
);
927 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
928 mpz_add (SCM_I_BIG_MPZ (result
),
930 SCM_I_BIG_MPZ (result
));
931 scm_remember_upto_here_1 (y
);
932 /* and do this before the next one */
934 return scm_i_normbig (result
);
938 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
940 else if (SCM_BIGP (x
))
944 long yy
= SCM_I_INUM (y
);
946 scm_num_overflow (s_modulo
);
949 SCM result
= scm_i_mkbig ();
950 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
952 (yy
< 0) ? - yy
: yy
);
953 scm_remember_upto_here_1 (x
);
954 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
955 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
956 SCM_I_BIG_MPZ (result
),
958 return scm_i_normbig (result
);
961 else if (SCM_BIGP (y
))
964 SCM result
= scm_i_mkbig ();
965 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
966 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
967 mpz_mod (SCM_I_BIG_MPZ (result
),
969 SCM_I_BIG_MPZ (pos_y
));
971 scm_remember_upto_here_1 (x
);
972 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
973 mpz_add (SCM_I_BIG_MPZ (result
),
975 SCM_I_BIG_MPZ (result
));
976 scm_remember_upto_here_2 (y
, pos_y
);
977 return scm_i_normbig (result
);
981 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
984 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
987 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
988 /* "Return the greatest common divisor of all arguments.\n"
989 * "If called without arguments, 0 is returned."
992 scm_gcd (SCM x
, SCM y
)
995 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1001 long xx
= SCM_I_INUM (x
);
1002 long yy
= SCM_I_INUM (y
);
1003 long u
= xx
< 0 ? -xx
: xx
;
1004 long v
= yy
< 0 ? -yy
: yy
;
1014 /* Determine a common factor 2^k */
1015 while (!(1 & (u
| v
)))
1021 /* Now, any factor 2^n can be eliminated */
1041 return (SCM_POSFIXABLE (result
)
1042 ? SCM_I_MAKINUM (result
)
1043 : scm_i_long2big (result
));
1045 else if (SCM_BIGP (y
))
1051 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1053 else if (SCM_BIGP (x
))
1055 if (SCM_I_INUMP (y
))
1057 unsigned long result
;
1060 yy
= SCM_I_INUM (y
);
1065 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1066 scm_remember_upto_here_1 (x
);
1067 return (SCM_POSFIXABLE (result
)
1068 ? SCM_I_MAKINUM (result
)
1069 : scm_from_ulong (result
));
1071 else if (SCM_BIGP (y
))
1073 SCM result
= scm_i_mkbig ();
1074 mpz_gcd (SCM_I_BIG_MPZ (result
),
1077 scm_remember_upto_here_2 (x
, y
);
1078 return scm_i_normbig (result
);
1081 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1084 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1087 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1088 /* "Return the least common multiple of the arguments.\n"
1089 * "If called without arguments, 1 is returned."
1092 scm_lcm (SCM n1
, SCM n2
)
1094 if (SCM_UNBNDP (n2
))
1096 if (SCM_UNBNDP (n1
))
1097 return SCM_I_MAKINUM (1L);
1098 n2
= SCM_I_MAKINUM (1L);
1101 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1102 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1103 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1104 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1106 if (SCM_I_INUMP (n1
))
1108 if (SCM_I_INUMP (n2
))
1110 SCM d
= scm_gcd (n1
, n2
);
1111 if (scm_is_eq (d
, SCM_INUM0
))
1114 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1118 /* inum n1, big n2 */
1121 SCM result
= scm_i_mkbig ();
1122 long nn1
= SCM_I_INUM (n1
);
1123 if (nn1
== 0) return SCM_INUM0
;
1124 if (nn1
< 0) nn1
= - nn1
;
1125 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1126 scm_remember_upto_here_1 (n2
);
1134 if (SCM_I_INUMP (n2
))
1141 SCM result
= scm_i_mkbig ();
1142 mpz_lcm(SCM_I_BIG_MPZ (result
),
1144 SCM_I_BIG_MPZ (n2
));
1145 scm_remember_upto_here_2(n1
, n2
);
1146 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1152 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1157 + + + x (map digit:logand X Y)
1158 + - + x (map digit:logand X (lognot (+ -1 Y)))
1159 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1160 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1165 + + + (map digit:logior X Y)
1166 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1167 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1168 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1173 + + + (map digit:logxor X Y)
1174 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1175 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1176 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1181 + + (any digit:logand X Y)
1182 + - (any digit:logand X (lognot (+ -1 Y)))
1183 - + (any digit:logand (lognot (+ -1 X)) Y)
1188 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1190 "Return the bitwise AND of the integer arguments.\n\n"
1192 "(logand) @result{} -1\n"
1193 "(logand 7) @result{} 7\n"
1194 "(logand #b111 #b011 #b001) @result{} 1\n"
1196 #define FUNC_NAME s_scm_logand
1200 if (SCM_UNBNDP (n2
))
1202 if (SCM_UNBNDP (n1
))
1203 return SCM_I_MAKINUM (-1);
1204 else if (!SCM_NUMBERP (n1
))
1205 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1206 else if (SCM_NUMBERP (n1
))
1209 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1212 if (SCM_I_INUMP (n1
))
1214 nn1
= SCM_I_INUM (n1
);
1215 if (SCM_I_INUMP (n2
))
1217 long nn2
= SCM_I_INUM (n2
);
1218 return SCM_I_MAKINUM (nn1
& nn2
);
1220 else if SCM_BIGP (n2
)
1226 SCM result_z
= scm_i_mkbig ();
1228 mpz_init_set_si (nn1_z
, nn1
);
1229 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1230 scm_remember_upto_here_1 (n2
);
1232 return scm_i_normbig (result_z
);
1236 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1238 else if (SCM_BIGP (n1
))
1240 if (SCM_I_INUMP (n2
))
1243 nn1
= SCM_I_INUM (n1
);
1246 else if (SCM_BIGP (n2
))
1248 SCM result_z
= scm_i_mkbig ();
1249 mpz_and (SCM_I_BIG_MPZ (result_z
),
1251 SCM_I_BIG_MPZ (n2
));
1252 scm_remember_upto_here_2 (n1
, n2
);
1253 return scm_i_normbig (result_z
);
1256 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1259 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1264 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1266 "Return the bitwise OR of the integer arguments.\n\n"
1268 "(logior) @result{} 0\n"
1269 "(logior 7) @result{} 7\n"
1270 "(logior #b000 #b001 #b011) @result{} 3\n"
1272 #define FUNC_NAME s_scm_logior
1276 if (SCM_UNBNDP (n2
))
1278 if (SCM_UNBNDP (n1
))
1280 else if (SCM_NUMBERP (n1
))
1283 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1286 if (SCM_I_INUMP (n1
))
1288 nn1
= SCM_I_INUM (n1
);
1289 if (SCM_I_INUMP (n2
))
1291 long nn2
= SCM_I_INUM (n2
);
1292 return SCM_I_MAKINUM (nn1
| nn2
);
1294 else if (SCM_BIGP (n2
))
1300 SCM result_z
= scm_i_mkbig ();
1302 mpz_init_set_si (nn1_z
, nn1
);
1303 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1304 scm_remember_upto_here_1 (n2
);
1310 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1312 else if (SCM_BIGP (n1
))
1314 if (SCM_I_INUMP (n2
))
1317 nn1
= SCM_I_INUM (n1
);
1320 else if (SCM_BIGP (n2
))
1322 SCM result_z
= scm_i_mkbig ();
1323 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1325 SCM_I_BIG_MPZ (n2
));
1326 scm_remember_upto_here_2 (n1
, n2
);
1330 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1333 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1338 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1340 "Return the bitwise XOR of the integer arguments. A bit is\n"
1341 "set in the result if it is set in an odd number of arguments.\n"
1343 "(logxor) @result{} 0\n"
1344 "(logxor 7) @result{} 7\n"
1345 "(logxor #b000 #b001 #b011) @result{} 2\n"
1346 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1348 #define FUNC_NAME s_scm_logxor
1352 if (SCM_UNBNDP (n2
))
1354 if (SCM_UNBNDP (n1
))
1356 else if (SCM_NUMBERP (n1
))
1359 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1362 if (SCM_I_INUMP (n1
))
1364 nn1
= SCM_I_INUM (n1
);
1365 if (SCM_I_INUMP (n2
))
1367 long nn2
= SCM_I_INUM (n2
);
1368 return SCM_I_MAKINUM (nn1
^ nn2
);
1370 else if (SCM_BIGP (n2
))
1374 SCM result_z
= scm_i_mkbig ();
1376 mpz_init_set_si (nn1_z
, nn1
);
1377 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1378 scm_remember_upto_here_1 (n2
);
1380 return scm_i_normbig (result_z
);
1384 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1386 else if (SCM_BIGP (n1
))
1388 if (SCM_I_INUMP (n2
))
1391 nn1
= SCM_I_INUM (n1
);
1394 else if (SCM_BIGP (n2
))
1396 SCM result_z
= scm_i_mkbig ();
1397 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1399 SCM_I_BIG_MPZ (n2
));
1400 scm_remember_upto_here_2 (n1
, n2
);
1401 return scm_i_normbig (result_z
);
1404 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1407 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1412 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1415 "(logtest j k) @equiv{} (not (zero? (logand j k)))\n\n"
1416 "(logtest #b0100 #b1011) @result{} #f\n"
1417 "(logtest #b0100 #b0111) @result{} #t\n"
1419 #define FUNC_NAME s_scm_logtest
1423 if (SCM_I_INUMP (j
))
1425 nj
= SCM_I_INUM (j
);
1426 if (SCM_I_INUMP (k
))
1428 long nk
= SCM_I_INUM (k
);
1429 return scm_from_bool (nj
& nk
);
1431 else if (SCM_BIGP (k
))
1439 mpz_init_set_si (nj_z
, nj
);
1440 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1441 scm_remember_upto_here_1 (k
);
1442 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1448 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1450 else if (SCM_BIGP (j
))
1452 if (SCM_I_INUMP (k
))
1455 nj
= SCM_I_INUM (j
);
1458 else if (SCM_BIGP (k
))
1462 mpz_init (result_z
);
1466 scm_remember_upto_here_2 (j
, k
);
1467 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1468 mpz_clear (result_z
);
1472 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1475 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1480 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1483 "(logbit? index j) @equiv{} (logtest (integer-expt 2 index) j)\n\n"
1484 "(logbit? 0 #b1101) @result{} #t\n"
1485 "(logbit? 1 #b1101) @result{} #f\n"
1486 "(logbit? 2 #b1101) @result{} #t\n"
1487 "(logbit? 3 #b1101) @result{} #t\n"
1488 "(logbit? 4 #b1101) @result{} #f\n"
1490 #define FUNC_NAME s_scm_logbit_p
1492 unsigned long int iindex
;
1493 iindex
= scm_to_ulong (index
);
1495 if (SCM_I_INUMP (j
))
1497 /* bits above what's in an inum follow the sign bit */
1498 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1499 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1501 else if (SCM_BIGP (j
))
1503 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1504 scm_remember_upto_here_1 (j
);
1505 return scm_from_bool (val
);
1508 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1513 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1515 "Return the integer which is the ones-complement of the integer\n"
1519 "(number->string (lognot #b10000000) 2)\n"
1520 " @result{} \"-10000001\"\n"
1521 "(number->string (lognot #b0) 2)\n"
1522 " @result{} \"-1\"\n"
1524 #define FUNC_NAME s_scm_lognot
1526 if (SCM_I_INUMP (n
)) {
1527 /* No overflow here, just need to toggle all the bits making up the inum.
1528 Enhancement: No need to strip the tag and add it back, could just xor
1529 a block of 1 bits, if that worked with the various debug versions of
1531 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1533 } else if (SCM_BIGP (n
)) {
1534 SCM result
= scm_i_mkbig ();
1535 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1536 scm_remember_upto_here_1 (n
);
1540 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1545 /* returns 0 if IN is not an integer. OUT must already be
1548 coerce_to_big (SCM in
, mpz_t out
)
1551 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1552 else if (SCM_I_INUMP (in
))
1553 mpz_set_si (out
, SCM_I_INUM (in
));
1560 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1561 (SCM n
, SCM k
, SCM m
),
1562 "Return @var{n} raised to the integer exponent\n"
1563 "@var{k}, modulo @var{m}.\n"
1566 "(modulo-expt 2 3 5)\n"
1569 #define FUNC_NAME s_scm_modulo_expt
1575 /* There are two classes of error we might encounter --
1576 1) Math errors, which we'll report by calling scm_num_overflow,
1578 2) wrong-type errors, which of course we'll report by calling
1580 We don't report those errors immediately, however; instead we do
1581 some cleanup first. These variables tell us which error (if
1582 any) we should report after cleaning up.
1584 int report_overflow
= 0;
1586 int position_of_wrong_type
= 0;
1587 SCM value_of_wrong_type
= SCM_INUM0
;
1589 SCM result
= SCM_UNDEFINED
;
1595 if (scm_is_eq (m
, SCM_INUM0
))
1597 report_overflow
= 1;
1601 if (!coerce_to_big (n
, n_tmp
))
1603 value_of_wrong_type
= n
;
1604 position_of_wrong_type
= 1;
1608 if (!coerce_to_big (k
, k_tmp
))
1610 value_of_wrong_type
= k
;
1611 position_of_wrong_type
= 2;
1615 if (!coerce_to_big (m
, m_tmp
))
1617 value_of_wrong_type
= m
;
1618 position_of_wrong_type
= 3;
1622 /* if the exponent K is negative, and we simply call mpz_powm, we
1623 will get a divide-by-zero exception when an inverse 1/n mod m
1624 doesn't exist (or is not unique). Since exceptions are hard to
1625 handle, we'll attempt the inversion "by hand" -- that way, we get
1626 a simple failure code, which is easy to handle. */
1628 if (-1 == mpz_sgn (k_tmp
))
1630 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1632 report_overflow
= 1;
1635 mpz_neg (k_tmp
, k_tmp
);
1638 result
= scm_i_mkbig ();
1639 mpz_powm (SCM_I_BIG_MPZ (result
),
1644 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1645 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1652 if (report_overflow
)
1653 scm_num_overflow (FUNC_NAME
);
1655 if (position_of_wrong_type
)
1656 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1657 value_of_wrong_type
);
1659 return scm_i_normbig (result
);
1663 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1665 "Return @var{n} raised to the non-negative integer exponent\n"
1669 "(integer-expt 2 5)\n"
1671 "(integer-expt -3 3)\n"
1674 #define FUNC_NAME s_scm_integer_expt
1677 SCM z_i2
= SCM_BOOL_F
;
1679 SCM acc
= SCM_I_MAKINUM (1L);
1681 /* 0^0 == 1 according to R5RS */
1682 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1683 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1684 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1685 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1687 if (SCM_I_INUMP (k
))
1688 i2
= SCM_I_INUM (k
);
1689 else if (SCM_BIGP (k
))
1691 z_i2
= scm_i_clonebig (k
, 1);
1692 scm_remember_upto_here_1 (k
);
1695 else if (SCM_REALP (k
))
1697 double r
= SCM_REAL_VALUE (k
);
1699 SCM_WRONG_TYPE_ARG (2, k
);
1700 if ((r
> SCM_MOST_POSITIVE_FIXNUM
) || (r
< SCM_MOST_NEGATIVE_FIXNUM
))
1702 z_i2
= scm_i_mkbig ();
1703 mpz_set_d (SCM_I_BIG_MPZ (z_i2
), r
);
1712 SCM_WRONG_TYPE_ARG (2, k
);
1716 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1718 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1719 n
= scm_divide (n
, SCM_UNDEFINED
);
1723 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1727 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1729 return scm_product (acc
, n
);
1731 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1732 acc
= scm_product (acc
, n
);
1733 n
= scm_product (n
, n
);
1734 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1742 n
= scm_divide (n
, SCM_UNDEFINED
);
1749 return scm_product (acc
, n
);
1751 acc
= scm_product (acc
, n
);
1752 n
= scm_product (n
, n
);
1759 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1761 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1762 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1764 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1765 "@var{cnt} is negative it's a division, rounded towards negative\n"
1766 "infinity. (Note that this is not the same rounding as\n"
1767 "@code{quotient} does.)\n"
1769 "With @var{n} viewed as an infinite precision twos complement,\n"
1770 "@code{ash} means a left shift introducing zero bits, or a right\n"
1771 "shift dropping bits.\n"
1774 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1775 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1777 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1778 "(ash -23 -2) @result{} -6\n"
1780 #define FUNC_NAME s_scm_ash
1783 bits_to_shift
= scm_to_long (cnt
);
1785 if (bits_to_shift
< 0)
1787 /* Shift right by abs(cnt) bits. This is realized as a division
1788 by div:=2^abs(cnt). However, to guarantee the floor
1789 rounding, negative values require some special treatment.
1791 SCM div
= scm_integer_expt (SCM_I_MAKINUM (2),
1792 scm_from_long (-bits_to_shift
));
1794 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1795 if (scm_is_false (scm_negative_p (n
)))
1796 return scm_quotient (n
, div
);
1798 return scm_sum (SCM_I_MAKINUM (-1L),
1799 scm_quotient (scm_sum (SCM_I_MAKINUM (1L), n
), div
));
1802 /* Shift left is done by multiplication with 2^CNT */
1803 return scm_product (n
, scm_integer_expt (SCM_I_MAKINUM (2), cnt
));
1808 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1809 (SCM n
, SCM start
, SCM end
),
1810 "Return the integer composed of the @var{start} (inclusive)\n"
1811 "through @var{end} (exclusive) bits of @var{n}. The\n"
1812 "@var{start}th bit becomes the 0-th bit in the result.\n"
1815 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1816 " @result{} \"1010\"\n"
1817 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1818 " @result{} \"10110\"\n"
1820 #define FUNC_NAME s_scm_bit_extract
1822 unsigned long int istart
, iend
, bits
;
1823 istart
= scm_to_ulong (start
);
1824 iend
= scm_to_ulong (end
);
1825 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1827 /* how many bits to keep */
1828 bits
= iend
- istart
;
1830 if (SCM_I_INUMP (n
))
1832 long int in
= SCM_I_INUM (n
);
1834 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1835 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1836 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1838 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1840 /* Since we emulate two's complement encoded numbers, this
1841 * special case requires us to produce a result that has
1842 * more bits than can be stored in a fixnum.
1844 SCM result
= scm_i_long2big (in
);
1845 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1850 /* mask down to requisite bits */
1851 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1852 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1854 else if (SCM_BIGP (n
))
1859 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1863 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1864 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1865 such bits into a ulong. */
1866 result
= scm_i_mkbig ();
1867 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1868 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1869 result
= scm_i_normbig (result
);
1871 scm_remember_upto_here_1 (n
);
1875 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1880 static const char scm_logtab
[] = {
1881 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1884 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1886 "Return the number of bits in integer @var{n}. If integer is\n"
1887 "positive, the 1-bits in its binary representation are counted.\n"
1888 "If negative, the 0-bits in its two's-complement binary\n"
1889 "representation are counted. If 0, 0 is returned.\n"
1892 "(logcount #b10101010)\n"
1899 #define FUNC_NAME s_scm_logcount
1901 if (SCM_I_INUMP (n
))
1903 unsigned long int c
= 0;
1904 long int nn
= SCM_I_INUM (n
);
1909 c
+= scm_logtab
[15 & nn
];
1912 return SCM_I_MAKINUM (c
);
1914 else if (SCM_BIGP (n
))
1916 unsigned long count
;
1917 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1918 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1920 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1921 scm_remember_upto_here_1 (n
);
1922 return SCM_I_MAKINUM (count
);
1925 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1930 static const char scm_ilentab
[] = {
1931 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1935 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1937 "Return the number of bits necessary to represent @var{n}.\n"
1940 "(integer-length #b10101010)\n"
1942 "(integer-length 0)\n"
1944 "(integer-length #b1111)\n"
1947 #define FUNC_NAME s_scm_integer_length
1949 if (SCM_I_INUMP (n
))
1951 unsigned long int c
= 0;
1953 long int nn
= SCM_I_INUM (n
);
1959 l
= scm_ilentab
[15 & nn
];
1962 return SCM_I_MAKINUM (c
- 4 + l
);
1964 else if (SCM_BIGP (n
))
1966 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1967 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1968 1 too big, so check for that and adjust. */
1969 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1970 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1971 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1972 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1974 scm_remember_upto_here_1 (n
);
1975 return SCM_I_MAKINUM (size
);
1978 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1982 /*** NUMBERS -> STRINGS ***/
1983 #define SCM_MAX_DBL_PREC 60
1984 #define SCM_MAX_DBL_RADIX 36
1986 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
1987 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
1988 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
1991 void init_dblprec(int *prec
, int radix
) {
1992 /* determine floating point precision by adding successively
1993 smaller increments to 1.0 until it is considered == 1.0 */
1994 double f
= ((double)1.0)/radix
;
1995 double fsum
= 1.0 + f
;
2000 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2012 void init_fx_radix(double *fx_list
, int radix
)
2014 /* initialize a per-radix list of tolerances. When added
2015 to a number < 1.0, we can determine if we should raund
2016 up and quit converting a number to a string. */
2020 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2021 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2024 /* use this array as a way to generate a single digit */
2025 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2028 idbl2str (double f
, char *a
, int radix
)
2030 int efmt
, dpt
, d
, i
, wp
;
2032 #ifdef DBL_MIN_10_EXP
2035 #endif /* DBL_MIN_10_EXP */
2040 radix
> SCM_MAX_DBL_RADIX
)
2042 /* revert to existing behavior */
2046 wp
= scm_dblprec
[radix
-2];
2047 fx
= fx_per_radix
[radix
-2];
2051 #ifdef HAVE_COPYSIGN
2052 double sgn
= copysign (1.0, f
);
2057 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2063 strcpy (a
, "-inf.0");
2065 strcpy (a
, "+inf.0");
2068 else if (xisnan (f
))
2070 strcpy (a
, "+nan.0");
2080 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2081 make-uniform-vector, from causing infinite loops. */
2082 /* just do the checking...if it passes, we do the conversion for our
2083 radix again below */
2090 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2098 while (f_cpy
> 10.0)
2101 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2122 if (f
+ fx
[wp
] >= radix
)
2129 /* adding 9999 makes this equivalent to abs(x) % 3 */
2130 dpt
= (exp
+ 9999) % 3;
2134 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2156 a
[ch
++] = number_chars
[d
];
2159 if (f
+ fx
[wp
] >= 1.0)
2161 a
[ch
- 1] = number_chars
[d
+1];
2173 if ((dpt
> 4) && (exp
> 6))
2175 d
= (a
[0] == '-' ? 2 : 1);
2176 for (i
= ch
++; i
> d
; i
--)
2189 if (a
[ch
- 1] == '.')
2190 a
[ch
++] = '0'; /* trailing zero */
2199 for (i
= radix
; i
<= exp
; i
*= radix
);
2200 for (i
/= radix
; i
; i
/= radix
)
2202 a
[ch
++] = number_chars
[exp
/ i
];
2210 iflo2str (SCM flt
, char *str
, int radix
)
2213 if (SCM_REALP (flt
))
2214 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2217 i
= idbl2str (SCM_COMPLEX_REAL (flt
), str
, radix
);
2218 if (SCM_COMPLEX_IMAG (flt
) != 0.0)
2220 double imag
= SCM_COMPLEX_IMAG (flt
);
2221 /* Don't output a '+' for negative numbers or for Inf and
2222 NaN. They will provide their own sign. */
2223 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2225 i
+= idbl2str (imag
, &str
[i
], radix
);
2232 /* convert a long to a string (unterminated). returns the number of
2233 characters in the result.
2235 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2237 scm_iint2str (long num
, int rad
, char *p
)
2241 unsigned long n
= (num
< 0) ? -num
: num
;
2243 for (n
/= rad
; n
> 0; n
/= rad
)
2260 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2265 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2267 "Return a string holding the external representation of the\n"
2268 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2269 "inexact, a radix of 10 will be used.")
2270 #define FUNC_NAME s_scm_number_to_string
2274 if (SCM_UNBNDP (radix
))
2277 base
= scm_to_signed_integer (radix
, 2, 36);
2279 if (SCM_I_INUMP (n
))
2281 char num_buf
[SCM_INTBUFLEN
];
2282 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2283 return scm_mem2string (num_buf
, length
);
2285 else if (SCM_BIGP (n
))
2287 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2288 scm_remember_upto_here_1 (n
);
2289 return scm_take0str (str
);
2291 else if (SCM_FRACTIONP (n
))
2293 scm_i_fraction_reduce (n
);
2294 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2295 scm_mem2string ("/", 1),
2296 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2298 else if (SCM_INEXACTP (n
))
2300 char num_buf
[FLOBUFLEN
];
2301 return scm_mem2string (num_buf
, iflo2str (n
, num_buf
, base
));
2304 SCM_WRONG_TYPE_ARG (1, n
);
2309 /* These print routines used to be stubbed here so that scm_repl.c
2310 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2313 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2315 char num_buf
[FLOBUFLEN
];
2316 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2321 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2324 char num_buf
[FLOBUFLEN
];
2325 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2330 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2333 scm_i_fraction_reduce (sexp
);
2334 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2335 scm_lfwrite (SCM_STRING_CHARS (str
), SCM_STRING_LENGTH (str
), port
);
2336 scm_remember_upto_here_1 (str
);
2341 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2343 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2344 scm_remember_upto_here_1 (exp
);
2345 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2349 /*** END nums->strs ***/
2352 /*** STRINGS -> NUMBERS ***/
2354 /* The following functions implement the conversion from strings to numbers.
2355 * The implementation somehow follows the grammar for numbers as it is given
2356 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2357 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2358 * points should be noted about the implementation:
2359 * * Each function keeps a local index variable 'idx' that points at the
2360 * current position within the parsed string. The global index is only
2361 * updated if the function could parse the corresponding syntactic unit
2363 * * Similarly, the functions keep track of indicators of inexactness ('#',
2364 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2365 * global exactness information is only updated after each part has been
2366 * successfully parsed.
2367 * * Sequences of digits are parsed into temporary variables holding fixnums.
2368 * Only if these fixnums would overflow, the result variables are updated
2369 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2370 * the temporary variables holding the fixnums are cleared, and the process
2371 * starts over again. If for example fixnums were able to store five decimal
2372 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2373 * and the result was computed as 12345 * 100000 + 67890. In other words,
2374 * only every five digits two bignum operations were performed.
2377 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2379 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2381 /* In non ASCII-style encodings the following macro might not work. */
2382 #define XDIGIT2UINT(d) \
2383 (isdigit ((int) (unsigned char) d) \
2385 : tolower ((int) (unsigned char) d) - 'a' + 10)
2388 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2389 unsigned int radix
, enum t_exactness
*p_exactness
)
2391 unsigned int idx
= *p_idx
;
2392 unsigned int hash_seen
= 0;
2393 scm_t_bits shift
= 1;
2395 unsigned int digit_value
;
2403 if (!isxdigit ((int) (unsigned char) c
))
2405 digit_value
= XDIGIT2UINT (c
);
2406 if (digit_value
>= radix
)
2410 result
= SCM_I_MAKINUM (digit_value
);
2414 if (isxdigit ((int) (unsigned char) c
))
2418 digit_value
= XDIGIT2UINT (c
);
2419 if (digit_value
>= radix
)
2431 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2433 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2435 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2442 shift
= shift
* radix
;
2443 add
= add
* radix
+ digit_value
;
2448 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2450 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2454 *p_exactness
= INEXACT
;
2460 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2461 * covers the parts of the rules that start at a potential point. The value
2462 * of the digits up to the point have been parsed by the caller and are given
2463 * in variable result. The content of *p_exactness indicates, whether a hash
2464 * has already been seen in the digits before the point.
2467 /* In non ASCII-style encodings the following macro might not work. */
2468 #define DIGIT2UINT(d) ((d) - '0')
2471 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2472 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2474 unsigned int idx
= *p_idx
;
2475 enum t_exactness x
= *p_exactness
;
2480 if (mem
[idx
] == '.')
2482 scm_t_bits shift
= 1;
2484 unsigned int digit_value
;
2485 SCM big_shift
= SCM_I_MAKINUM (1);
2491 if (isdigit ((int) (unsigned char) c
))
2496 digit_value
= DIGIT2UINT (c
);
2507 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2509 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2510 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2512 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2520 add
= add
* 10 + digit_value
;
2526 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2527 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2528 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2531 result
= scm_divide (result
, big_shift
);
2533 /* We've seen a decimal point, thus the value is implicitly inexact. */
2545 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2572 if (!isdigit ((int) (unsigned char) c
))
2576 exponent
= DIGIT2UINT (c
);
2580 if (isdigit ((int) (unsigned char) c
))
2583 if (exponent
<= SCM_MAXEXP
)
2584 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2590 if (exponent
> SCM_MAXEXP
)
2592 size_t exp_len
= idx
- start
;
2593 SCM exp_string
= scm_mem2string (&mem
[start
], exp_len
);
2594 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2595 scm_out_of_range ("string->number", exp_num
);
2598 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2600 result
= scm_product (result
, e
);
2602 result
= scm_divide2real (result
, e
);
2604 /* We've seen an exponent, thus the value is implicitly inexact. */
2622 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2625 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2626 unsigned int radix
, enum t_exactness
*p_exactness
)
2628 unsigned int idx
= *p_idx
;
2634 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2640 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2642 enum t_exactness x
= EXACT
;
2644 /* Cobble up the fractional part. We might want to set the
2645 NaN's mantissa from it. */
2647 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2652 if (mem
[idx
] == '.')
2656 else if (idx
+ 1 == len
)
2658 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2661 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2662 p_idx
, p_exactness
);
2666 enum t_exactness x
= EXACT
;
2669 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2670 if (scm_is_false (uinteger
))
2675 else if (mem
[idx
] == '/')
2681 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2682 if (scm_is_false (divisor
))
2685 /* both are int/big here, I assume */
2686 result
= scm_make_ratio (uinteger
, divisor
);
2688 else if (radix
== 10)
2690 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2691 if (scm_is_false (result
))
2702 /* When returning an inexact zero, make sure it is represented as a
2703 floating point value so that we can change its sign.
2705 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2706 result
= scm_make_real (0.0);
2712 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2715 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2716 unsigned int radix
, enum t_exactness
*p_exactness
)
2740 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2741 if (scm_is_false (ureal
))
2743 /* input must be either +i or -i */
2748 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2754 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2761 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2762 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2771 /* either +<ureal>i or -<ureal>i */
2778 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2781 /* polar input: <real>@<real>. */
2806 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2807 if (scm_is_false (angle
))
2812 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2813 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2815 result
= scm_make_polar (ureal
, angle
);
2820 /* expecting input matching <real>[+-]<ureal>?i */
2827 int sign
= (c
== '+') ? 1 : -1;
2828 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2830 if (scm_is_false (imag
))
2831 imag
= SCM_I_MAKINUM (sign
);
2832 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2833 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2837 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2844 return scm_make_rectangular (ureal
, imag
);
2853 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2855 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2858 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2860 unsigned int idx
= 0;
2861 unsigned int radix
= NO_RADIX
;
2862 enum t_exactness forced_x
= NO_EXACTNESS
;
2863 enum t_exactness implicit_x
= EXACT
;
2866 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2867 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2869 switch (mem
[idx
+ 1])
2872 if (radix
!= NO_RADIX
)
2877 if (radix
!= NO_RADIX
)
2882 if (forced_x
!= NO_EXACTNESS
)
2887 if (forced_x
!= NO_EXACTNESS
)
2892 if (radix
!= NO_RADIX
)
2897 if (radix
!= NO_RADIX
)
2907 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2908 if (radix
== NO_RADIX
)
2909 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2911 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2913 if (scm_is_false (result
))
2919 if (SCM_INEXACTP (result
))
2920 return scm_inexact_to_exact (result
);
2924 if (SCM_INEXACTP (result
))
2927 return scm_exact_to_inexact (result
);
2930 if (implicit_x
== INEXACT
)
2932 if (SCM_INEXACTP (result
))
2935 return scm_exact_to_inexact (result
);
2943 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2944 (SCM string
, SCM radix
),
2945 "Return a number of the maximally precise representation\n"
2946 "expressed by the given @var{string}. @var{radix} must be an\n"
2947 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2948 "is a default radix that may be overridden by an explicit radix\n"
2949 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2950 "supplied, then the default radix is 10. If string is not a\n"
2951 "syntactically valid notation for a number, then\n"
2952 "@code{string->number} returns @code{#f}.")
2953 #define FUNC_NAME s_scm_string_to_number
2957 SCM_VALIDATE_STRING (1, string
);
2959 if (SCM_UNBNDP (radix
))
2962 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
2964 answer
= scm_i_mem2number (SCM_STRING_CHARS (string
),
2965 SCM_STRING_LENGTH (string
),
2967 return scm_return_first (answer
, string
);
2972 /*** END strs->nums ***/
2976 scm_make_real (double x
)
2978 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
2980 SCM_REAL_VALUE (z
) = x
;
2986 scm_make_complex (double x
, double y
)
2989 return scm_make_real (x
);
2993 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
2995 SCM_COMPLEX_REAL (z
) = x
;
2996 SCM_COMPLEX_IMAG (z
) = y
;
3003 scm_bigequal (SCM x
, SCM y
)
3005 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3006 scm_remember_upto_here_2 (x
, y
);
3007 return scm_from_bool (0 == result
);
3011 scm_real_equalp (SCM x
, SCM y
)
3013 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3017 scm_complex_equalp (SCM x
, SCM y
)
3019 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3020 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3024 scm_i_fraction_equalp (SCM x
, SCM y
)
3026 scm_i_fraction_reduce (x
);
3027 scm_i_fraction_reduce (y
);
3028 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3029 SCM_FRACTION_NUMERATOR (y
)))
3030 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3031 SCM_FRACTION_DENOMINATOR (y
))))
3038 SCM_REGISTER_PROC (s_number_p
, "number?", 1, 0, 0, scm_number_p
);
3039 /* "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3040 * "else. Note that the sets of complex, real, rational and\n"
3041 * "integer values form subsets of the set of numbers, i. e. the\n"
3042 * "predicate will be fulfilled for any number."
3044 SCM_DEFINE (scm_number_p
, "complex?", 1, 0, 0,
3046 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3047 "otherwise. Note that the sets of real, rational and integer\n"
3048 "values form subsets of the set of complex numbers, i. e. the\n"
3049 "predicate will also be fulfilled if @var{x} is a real,\n"
3050 "rational or integer number.")
3051 #define FUNC_NAME s_scm_number_p
3053 return scm_from_bool (SCM_NUMBERP (x
));
3058 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3060 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3061 "otherwise. Note that the set of integer values forms a subset of\n"
3062 "the set of real numbers, i. e. the predicate will also be\n"
3063 "fulfilled if @var{x} is an integer number.")
3064 #define FUNC_NAME s_scm_real_p
3066 /* we can't represent irrational numbers. */
3067 return scm_rational_p (x
);
3071 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3073 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3074 "otherwise. Note that the set of integer values forms a subset of\n"
3075 "the set of rational numbers, i. e. the predicate will also be\n"
3076 "fulfilled if @var{x} is an integer number.")
3077 #define FUNC_NAME s_scm_rational_p
3079 if (SCM_I_INUMP (x
))
3081 else if (SCM_IMP (x
))
3083 else if (SCM_BIGP (x
))
3085 else if (SCM_FRACTIONP (x
))
3087 else if (SCM_REALP (x
))
3088 /* due to their limited precision, all floating point numbers are
3089 rational as well. */
3097 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3099 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3101 #define FUNC_NAME s_scm_integer_p
3104 if (SCM_I_INUMP (x
))
3110 if (!SCM_INEXACTP (x
))
3112 if (SCM_COMPLEXP (x
))
3114 r
= SCM_REAL_VALUE (x
);
3122 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3124 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3126 #define FUNC_NAME s_scm_inexact_p
3128 if (SCM_INEXACTP (x
))
3130 if (SCM_NUMBERP (x
))
3132 SCM_WRONG_TYPE_ARG (1, x
);
3137 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3138 /* "Return @code{#t} if all parameters are numerically equal." */
3140 scm_num_eq_p (SCM x
, SCM y
)
3143 if (SCM_I_INUMP (x
))
3145 long xx
= SCM_I_INUM (x
);
3146 if (SCM_I_INUMP (y
))
3148 long yy
= SCM_I_INUM (y
);
3149 return scm_from_bool (xx
== yy
);
3151 else if (SCM_BIGP (y
))
3153 else if (SCM_REALP (y
))
3154 return scm_from_bool ((double) xx
== SCM_REAL_VALUE (y
));
3155 else if (SCM_COMPLEXP (y
))
3156 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3157 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3158 else if (SCM_FRACTIONP (y
))
3161 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3163 else if (SCM_BIGP (x
))
3165 if (SCM_I_INUMP (y
))
3167 else if (SCM_BIGP (y
))
3169 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3170 scm_remember_upto_here_2 (x
, y
);
3171 return scm_from_bool (0 == cmp
);
3173 else if (SCM_REALP (y
))
3176 if (xisnan (SCM_REAL_VALUE (y
)))
3178 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3179 scm_remember_upto_here_1 (x
);
3180 return scm_from_bool (0 == cmp
);
3182 else if (SCM_COMPLEXP (y
))
3185 if (0.0 != SCM_COMPLEX_IMAG (y
))
3187 if (xisnan (SCM_COMPLEX_REAL (y
)))
3189 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3190 scm_remember_upto_here_1 (x
);
3191 return scm_from_bool (0 == cmp
);
3193 else if (SCM_FRACTIONP (y
))
3196 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3198 else if (SCM_REALP (x
))
3200 if (SCM_I_INUMP (y
))
3201 return scm_from_bool (SCM_REAL_VALUE (x
) == (double) SCM_I_INUM (y
));
3202 else if (SCM_BIGP (y
))
3205 if (xisnan (SCM_REAL_VALUE (x
)))
3207 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3208 scm_remember_upto_here_1 (y
);
3209 return scm_from_bool (0 == cmp
);
3211 else if (SCM_REALP (y
))
3212 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3213 else if (SCM_COMPLEXP (y
))
3214 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3215 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3216 else if (SCM_FRACTIONP (y
))
3218 double xx
= SCM_REAL_VALUE (x
);
3222 return scm_from_bool (xx
< 0.0);
3223 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3227 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3229 else if (SCM_COMPLEXP (x
))
3231 if (SCM_I_INUMP (y
))
3232 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3233 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3234 else if (SCM_BIGP (y
))
3237 if (0.0 != SCM_COMPLEX_IMAG (x
))
3239 if (xisnan (SCM_COMPLEX_REAL (x
)))
3241 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3242 scm_remember_upto_here_1 (y
);
3243 return scm_from_bool (0 == cmp
);
3245 else if (SCM_REALP (y
))
3246 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3247 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3248 else if (SCM_COMPLEXP (y
))
3249 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3250 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3251 else if (SCM_FRACTIONP (y
))
3254 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3256 xx
= SCM_COMPLEX_REAL (x
);
3260 return scm_from_bool (xx
< 0.0);
3261 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3265 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3267 else if (SCM_FRACTIONP (x
))
3269 if (SCM_I_INUMP (y
))
3271 else if (SCM_BIGP (y
))
3273 else if (SCM_REALP (y
))
3275 double yy
= SCM_REAL_VALUE (y
);
3279 return scm_from_bool (0.0 < yy
);
3280 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3283 else if (SCM_COMPLEXP (y
))
3286 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3288 yy
= SCM_COMPLEX_REAL (y
);
3292 return scm_from_bool (0.0 < yy
);
3293 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3296 else if (SCM_FRACTIONP (y
))
3297 return scm_i_fraction_equalp (x
, y
);
3299 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3302 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3306 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3307 done are good for inums, but for bignums an answer can almost always be
3308 had by just examining a few high bits of the operands, as done by GMP in
3309 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3310 of the float exponent to take into account. */
3312 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3313 /* "Return @code{#t} if the list of parameters is monotonically\n"
3317 scm_less_p (SCM x
, SCM y
)
3320 if (SCM_I_INUMP (x
))
3322 long xx
= SCM_I_INUM (x
);
3323 if (SCM_I_INUMP (y
))
3325 long yy
= SCM_I_INUM (y
);
3326 return scm_from_bool (xx
< yy
);
3328 else if (SCM_BIGP (y
))
3330 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3331 scm_remember_upto_here_1 (y
);
3332 return scm_from_bool (sgn
> 0);
3334 else if (SCM_REALP (y
))
3335 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3336 else if (SCM_FRACTIONP (y
))
3338 /* "x < a/b" becomes "x*b < a" */
3340 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3341 y
= SCM_FRACTION_NUMERATOR (y
);
3345 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3347 else if (SCM_BIGP (x
))
3349 if (SCM_I_INUMP (y
))
3351 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3352 scm_remember_upto_here_1 (x
);
3353 return scm_from_bool (sgn
< 0);
3355 else if (SCM_BIGP (y
))
3357 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3358 scm_remember_upto_here_2 (x
, y
);
3359 return scm_from_bool (cmp
< 0);
3361 else if (SCM_REALP (y
))
3364 if (xisnan (SCM_REAL_VALUE (y
)))
3366 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3367 scm_remember_upto_here_1 (x
);
3368 return scm_from_bool (cmp
< 0);
3370 else if (SCM_FRACTIONP (y
))
3373 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3375 else if (SCM_REALP (x
))
3377 if (SCM_I_INUMP (y
))
3378 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3379 else if (SCM_BIGP (y
))
3382 if (xisnan (SCM_REAL_VALUE (x
)))
3384 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3385 scm_remember_upto_here_1 (y
);
3386 return scm_from_bool (cmp
> 0);
3388 else if (SCM_REALP (y
))
3389 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3390 else if (SCM_FRACTIONP (y
))
3392 double xx
= SCM_REAL_VALUE (x
);
3396 return scm_from_bool (xx
< 0.0);
3397 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3401 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3403 else if (SCM_FRACTIONP (x
))
3405 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3407 /* "a/b < y" becomes "a < y*b" */
3408 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3409 x
= SCM_FRACTION_NUMERATOR (x
);
3412 else if (SCM_REALP (y
))
3414 double yy
= SCM_REAL_VALUE (y
);
3418 return scm_from_bool (0.0 < yy
);
3419 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3422 else if (SCM_FRACTIONP (y
))
3424 /* "a/b < c/d" becomes "a*d < c*b" */
3425 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3426 SCM_FRACTION_DENOMINATOR (y
));
3427 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3428 SCM_FRACTION_DENOMINATOR (x
));
3434 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3437 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3441 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3442 /* "Return @code{#t} if the list of parameters is monotonically\n"
3445 #define FUNC_NAME s_scm_gr_p
3447 scm_gr_p (SCM x
, SCM y
)
3449 if (!SCM_NUMBERP (x
))
3450 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3451 else if (!SCM_NUMBERP (y
))
3452 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3454 return scm_less_p (y
, x
);
3459 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3460 /* "Return @code{#t} if the list of parameters is monotonically\n"
3463 #define FUNC_NAME s_scm_leq_p
3465 scm_leq_p (SCM x
, SCM y
)
3467 if (!SCM_NUMBERP (x
))
3468 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3469 else if (!SCM_NUMBERP (y
))
3470 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3471 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3474 return scm_not (scm_less_p (y
, x
));
3479 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3480 /* "Return @code{#t} if the list of parameters is monotonically\n"
3483 #define FUNC_NAME s_scm_geq_p
3485 scm_geq_p (SCM x
, SCM y
)
3487 if (!SCM_NUMBERP (x
))
3488 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3489 else if (!SCM_NUMBERP (y
))
3490 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3491 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3494 return scm_not (scm_less_p (x
, y
));
3499 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3500 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3506 if (SCM_I_INUMP (z
))
3507 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3508 else if (SCM_BIGP (z
))
3510 else if (SCM_REALP (z
))
3511 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3512 else if (SCM_COMPLEXP (z
))
3513 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3514 && SCM_COMPLEX_IMAG (z
) == 0.0);
3515 else if (SCM_FRACTIONP (z
))
3518 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3522 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3523 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3527 scm_positive_p (SCM x
)
3529 if (SCM_I_INUMP (x
))
3530 return scm_from_bool (SCM_I_INUM (x
) > 0);
3531 else if (SCM_BIGP (x
))
3533 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3534 scm_remember_upto_here_1 (x
);
3535 return scm_from_bool (sgn
> 0);
3537 else if (SCM_REALP (x
))
3538 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3539 else if (SCM_FRACTIONP (x
))
3540 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3542 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3546 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3547 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3551 scm_negative_p (SCM x
)
3553 if (SCM_I_INUMP (x
))
3554 return scm_from_bool (SCM_I_INUM (x
) < 0);
3555 else if (SCM_BIGP (x
))
3557 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3558 scm_remember_upto_here_1 (x
);
3559 return scm_from_bool (sgn
< 0);
3561 else if (SCM_REALP (x
))
3562 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3563 else if (SCM_FRACTIONP (x
))
3564 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3566 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3570 /* scm_min and scm_max return an inexact when either argument is inexact, as
3571 required by r5rs. On that basis, for exact/inexact combinations the
3572 exact is converted to inexact to compare and possibly return. This is
3573 unlike scm_less_p above which takes some trouble to preserve all bits in
3574 its test, such trouble is not required for min and max. */
3576 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3577 /* "Return the maximum of all parameter values."
3580 scm_max (SCM x
, SCM y
)
3585 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3586 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3589 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3592 if (SCM_I_INUMP (x
))
3594 long xx
= SCM_I_INUM (x
);
3595 if (SCM_I_INUMP (y
))
3597 long yy
= SCM_I_INUM (y
);
3598 return (xx
< yy
) ? y
: x
;
3600 else if (SCM_BIGP (y
))
3602 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3603 scm_remember_upto_here_1 (y
);
3604 return (sgn
< 0) ? x
: y
;
3606 else if (SCM_REALP (y
))
3609 /* if y==NaN then ">" is false and we return NaN */
3610 return (z
> SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3612 else if (SCM_FRACTIONP (y
))
3615 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3618 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3620 else if (SCM_BIGP (x
))
3622 if (SCM_I_INUMP (y
))
3624 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3625 scm_remember_upto_here_1 (x
);
3626 return (sgn
< 0) ? y
: x
;
3628 else if (SCM_BIGP (y
))
3630 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3631 scm_remember_upto_here_2 (x
, y
);
3632 return (cmp
> 0) ? x
: y
;
3634 else if (SCM_REALP (y
))
3636 /* if y==NaN then xx>yy is false, so we return the NaN y */
3639 xx
= scm_i_big2dbl (x
);
3640 yy
= SCM_REAL_VALUE (y
);
3641 return (xx
> yy
? scm_make_real (xx
) : y
);
3643 else if (SCM_FRACTIONP (y
))
3648 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3650 else if (SCM_REALP (x
))
3652 if (SCM_I_INUMP (y
))
3654 double z
= SCM_I_INUM (y
);
3655 /* if x==NaN then "<" is false and we return NaN */
3656 return (SCM_REAL_VALUE (x
) < z
) ? scm_make_real (z
) : x
;
3658 else if (SCM_BIGP (y
))
3663 else if (SCM_REALP (y
))
3665 /* if x==NaN then our explicit check means we return NaN
3666 if y==NaN then ">" is false and we return NaN
3667 calling isnan is unavoidable, since it's the only way to know
3668 which of x or y causes any compares to be false */
3669 double xx
= SCM_REAL_VALUE (x
);
3670 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3672 else if (SCM_FRACTIONP (y
))
3674 double yy
= scm_i_fraction2double (y
);
3675 double xx
= SCM_REAL_VALUE (x
);
3676 return (xx
< yy
) ? scm_make_real (yy
) : x
;
3679 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3681 else if (SCM_FRACTIONP (x
))
3683 if (SCM_I_INUMP (y
))
3687 else if (SCM_BIGP (y
))
3691 else if (SCM_REALP (y
))
3693 double xx
= scm_i_fraction2double (x
);
3694 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_make_real (xx
);
3696 else if (SCM_FRACTIONP (y
))
3701 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3704 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3708 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3709 /* "Return the minium of all parameter values."
3712 scm_min (SCM x
, SCM y
)
3717 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3718 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3721 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3724 if (SCM_I_INUMP (x
))
3726 long xx
= SCM_I_INUM (x
);
3727 if (SCM_I_INUMP (y
))
3729 long yy
= SCM_I_INUM (y
);
3730 return (xx
< yy
) ? x
: y
;
3732 else if (SCM_BIGP (y
))
3734 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3735 scm_remember_upto_here_1 (y
);
3736 return (sgn
< 0) ? y
: x
;
3738 else if (SCM_REALP (y
))
3741 /* if y==NaN then "<" is false and we return NaN */
3742 return (z
< SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3744 else if (SCM_FRACTIONP (y
))
3747 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3750 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3752 else if (SCM_BIGP (x
))
3754 if (SCM_I_INUMP (y
))
3756 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3757 scm_remember_upto_here_1 (x
);
3758 return (sgn
< 0) ? x
: y
;
3760 else if (SCM_BIGP (y
))
3762 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3763 scm_remember_upto_here_2 (x
, y
);
3764 return (cmp
> 0) ? y
: x
;
3766 else if (SCM_REALP (y
))
3768 /* if y==NaN then xx<yy is false, so we return the NaN y */
3771 xx
= scm_i_big2dbl (x
);
3772 yy
= SCM_REAL_VALUE (y
);
3773 return (xx
< yy
? scm_make_real (xx
) : y
);
3775 else if (SCM_FRACTIONP (y
))
3780 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3782 else if (SCM_REALP (x
))
3784 if (SCM_I_INUMP (y
))
3786 double z
= SCM_I_INUM (y
);
3787 /* if x==NaN then "<" is false and we return NaN */
3788 return (z
< SCM_REAL_VALUE (x
)) ? scm_make_real (z
) : x
;
3790 else if (SCM_BIGP (y
))
3795 else if (SCM_REALP (y
))
3797 /* if x==NaN then our explicit check means we return NaN
3798 if y==NaN then "<" is false and we return NaN
3799 calling isnan is unavoidable, since it's the only way to know
3800 which of x or y causes any compares to be false */
3801 double xx
= SCM_REAL_VALUE (x
);
3802 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3804 else if (SCM_FRACTIONP (y
))
3806 double yy
= scm_i_fraction2double (y
);
3807 double xx
= SCM_REAL_VALUE (x
);
3808 return (yy
< xx
) ? scm_make_real (yy
) : x
;
3811 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3813 else if (SCM_FRACTIONP (x
))
3815 if (SCM_I_INUMP (y
))
3819 else if (SCM_BIGP (y
))
3823 else if (SCM_REALP (y
))
3825 double xx
= scm_i_fraction2double (x
);
3826 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_make_real (xx
);
3828 else if (SCM_FRACTIONP (y
))
3833 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3836 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3840 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3841 /* "Return the sum of all parameter values. Return 0 if called without\n"
3845 scm_sum (SCM x
, SCM y
)
3849 if (SCM_NUMBERP (x
)) return x
;
3850 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3851 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3854 if (SCM_I_INUMP (x
))
3856 if (SCM_I_INUMP (y
))
3858 long xx
= SCM_I_INUM (x
);
3859 long yy
= SCM_I_INUM (y
);
3860 long int z
= xx
+ yy
;
3861 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3863 else if (SCM_BIGP (y
))
3868 else if (SCM_REALP (y
))
3870 long int xx
= SCM_I_INUM (x
);
3871 return scm_make_real (xx
+ SCM_REAL_VALUE (y
));
3873 else if (SCM_COMPLEXP (y
))
3875 long int xx
= SCM_I_INUM (x
);
3876 return scm_make_complex (xx
+ SCM_COMPLEX_REAL (y
),
3877 SCM_COMPLEX_IMAG (y
));
3879 else if (SCM_FRACTIONP (y
))
3880 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3881 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3882 SCM_FRACTION_DENOMINATOR (y
));
3884 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3885 } else if (SCM_BIGP (x
))
3887 if (SCM_I_INUMP (y
))
3892 inum
= SCM_I_INUM (y
);
3895 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3898 SCM result
= scm_i_mkbig ();
3899 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3900 scm_remember_upto_here_1 (x
);
3901 /* we know the result will have to be a bignum */
3904 return scm_i_normbig (result
);
3908 SCM result
= scm_i_mkbig ();
3909 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3910 scm_remember_upto_here_1 (x
);
3911 /* we know the result will have to be a bignum */
3914 return scm_i_normbig (result
);
3917 else if (SCM_BIGP (y
))
3919 SCM result
= scm_i_mkbig ();
3920 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3921 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3922 mpz_add (SCM_I_BIG_MPZ (result
),
3925 scm_remember_upto_here_2 (x
, y
);
3926 /* we know the result will have to be a bignum */
3929 return scm_i_normbig (result
);
3931 else if (SCM_REALP (y
))
3933 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3934 scm_remember_upto_here_1 (x
);
3935 return scm_make_real (result
);
3937 else if (SCM_COMPLEXP (y
))
3939 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3940 + SCM_COMPLEX_REAL (y
));
3941 scm_remember_upto_here_1 (x
);
3942 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
3944 else if (SCM_FRACTIONP (y
))
3945 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3946 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3947 SCM_FRACTION_DENOMINATOR (y
));
3949 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3951 else if (SCM_REALP (x
))
3953 if (SCM_I_INUMP (y
))
3954 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
3955 else if (SCM_BIGP (y
))
3957 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3958 scm_remember_upto_here_1 (y
);
3959 return scm_make_real (result
);
3961 else if (SCM_REALP (y
))
3962 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3963 else if (SCM_COMPLEXP (y
))
3964 return scm_make_complex (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3965 SCM_COMPLEX_IMAG (y
));
3966 else if (SCM_FRACTIONP (y
))
3967 return scm_make_real (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3969 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3971 else if (SCM_COMPLEXP (x
))
3973 if (SCM_I_INUMP (y
))
3974 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
3975 SCM_COMPLEX_IMAG (x
));
3976 else if (SCM_BIGP (y
))
3978 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3979 + SCM_COMPLEX_REAL (x
));
3980 scm_remember_upto_here_1 (y
);
3981 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (x
));
3983 else if (SCM_REALP (y
))
3984 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3985 SCM_COMPLEX_IMAG (x
));
3986 else if (SCM_COMPLEXP (y
))
3987 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3988 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3989 else if (SCM_FRACTIONP (y
))
3990 return scm_make_complex (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3991 SCM_COMPLEX_IMAG (x
));
3993 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3995 else if (SCM_FRACTIONP (x
))
3997 if (SCM_I_INUMP (y
))
3998 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3999 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4000 SCM_FRACTION_DENOMINATOR (x
));
4001 else if (SCM_BIGP (y
))
4002 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4003 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4004 SCM_FRACTION_DENOMINATOR (x
));
4005 else if (SCM_REALP (y
))
4006 return scm_make_real (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4007 else if (SCM_COMPLEXP (y
))
4008 return scm_make_complex (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4009 SCM_COMPLEX_IMAG (y
));
4010 else if (SCM_FRACTIONP (y
))
4011 /* a/b + c/d = (ad + bc) / bd */
4012 return scm_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4013 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4014 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4016 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4019 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4023 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4024 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4025 * the sum of all but the first argument are subtracted from the first
4027 #define FUNC_NAME s_difference
4029 scm_difference (SCM x
, SCM y
)
4034 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4036 if (SCM_I_INUMP (x
))
4038 long xx
= -SCM_I_INUM (x
);
4039 if (SCM_FIXABLE (xx
))
4040 return SCM_I_MAKINUM (xx
);
4042 return scm_i_long2big (xx
);
4044 else if (SCM_BIGP (x
))
4045 /* FIXME: do we really need to normalize here? */
4046 return scm_i_normbig (scm_i_clonebig (x
, 0));
4047 else if (SCM_REALP (x
))
4048 return scm_make_real (-SCM_REAL_VALUE (x
));
4049 else if (SCM_COMPLEXP (x
))
4050 return scm_make_complex (-SCM_COMPLEX_REAL (x
),
4051 -SCM_COMPLEX_IMAG (x
));
4052 else if (SCM_FRACTIONP (x
))
4053 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4054 SCM_FRACTION_DENOMINATOR (x
));
4056 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4059 if (SCM_I_INUMP (x
))
4061 if (SCM_I_INUMP (y
))
4063 long int xx
= SCM_I_INUM (x
);
4064 long int yy
= SCM_I_INUM (y
);
4065 long int z
= xx
- yy
;
4066 if (SCM_FIXABLE (z
))
4067 return SCM_I_MAKINUM (z
);
4069 return scm_i_long2big (z
);
4071 else if (SCM_BIGP (y
))
4073 /* inum-x - big-y */
4074 long xx
= SCM_I_INUM (x
);
4077 return scm_i_clonebig (y
, 0);
4080 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4081 SCM result
= scm_i_mkbig ();
4084 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4087 /* x - y == -(y + -x) */
4088 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4089 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4091 scm_remember_upto_here_1 (y
);
4093 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4094 /* we know the result will have to be a bignum */
4097 return scm_i_normbig (result
);
4100 else if (SCM_REALP (y
))
4102 long int xx
= SCM_I_INUM (x
);
4103 return scm_make_real (xx
- SCM_REAL_VALUE (y
));
4105 else if (SCM_COMPLEXP (y
))
4107 long int xx
= SCM_I_INUM (x
);
4108 return scm_make_complex (xx
- SCM_COMPLEX_REAL (y
),
4109 - SCM_COMPLEX_IMAG (y
));
4111 else if (SCM_FRACTIONP (y
))
4112 /* a - b/c = (ac - b) / c */
4113 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4114 SCM_FRACTION_NUMERATOR (y
)),
4115 SCM_FRACTION_DENOMINATOR (y
));
4117 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4119 else if (SCM_BIGP (x
))
4121 if (SCM_I_INUMP (y
))
4123 /* big-x - inum-y */
4124 long yy
= SCM_I_INUM (y
);
4125 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4127 scm_remember_upto_here_1 (x
);
4129 return (SCM_FIXABLE (-yy
) ?
4130 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4133 SCM result
= scm_i_mkbig ();
4136 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4138 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4139 scm_remember_upto_here_1 (x
);
4141 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4142 /* we know the result will have to be a bignum */
4145 return scm_i_normbig (result
);
4148 else if (SCM_BIGP (y
))
4150 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4151 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4152 SCM result
= scm_i_mkbig ();
4153 mpz_sub (SCM_I_BIG_MPZ (result
),
4156 scm_remember_upto_here_2 (x
, y
);
4157 /* we know the result will have to be a bignum */
4158 if ((sgn_x
== 1) && (sgn_y
== -1))
4160 if ((sgn_x
== -1) && (sgn_y
== 1))
4162 return scm_i_normbig (result
);
4164 else if (SCM_REALP (y
))
4166 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4167 scm_remember_upto_here_1 (x
);
4168 return scm_make_real (result
);
4170 else if (SCM_COMPLEXP (y
))
4172 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4173 - SCM_COMPLEX_REAL (y
));
4174 scm_remember_upto_here_1 (x
);
4175 return scm_make_complex (real_part
, - SCM_COMPLEX_IMAG (y
));
4177 else if (SCM_FRACTIONP (y
))
4178 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4179 SCM_FRACTION_NUMERATOR (y
)),
4180 SCM_FRACTION_DENOMINATOR (y
));
4181 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4183 else if (SCM_REALP (x
))
4185 if (SCM_I_INUMP (y
))
4186 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4187 else if (SCM_BIGP (y
))
4189 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4190 scm_remember_upto_here_1 (x
);
4191 return scm_make_real (result
);
4193 else if (SCM_REALP (y
))
4194 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4195 else if (SCM_COMPLEXP (y
))
4196 return scm_make_complex (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4197 -SCM_COMPLEX_IMAG (y
));
4198 else if (SCM_FRACTIONP (y
))
4199 return scm_make_real (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4201 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4203 else if (SCM_COMPLEXP (x
))
4205 if (SCM_I_INUMP (y
))
4206 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4207 SCM_COMPLEX_IMAG (x
));
4208 else if (SCM_BIGP (y
))
4210 double real_part
= (SCM_COMPLEX_REAL (x
)
4211 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4212 scm_remember_upto_here_1 (x
);
4213 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
4215 else if (SCM_REALP (y
))
4216 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4217 SCM_COMPLEX_IMAG (x
));
4218 else if (SCM_COMPLEXP (y
))
4219 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4220 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4221 else if (SCM_FRACTIONP (y
))
4222 return scm_make_complex (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4223 SCM_COMPLEX_IMAG (x
));
4225 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4227 else if (SCM_FRACTIONP (x
))
4229 if (SCM_I_INUMP (y
))
4230 /* a/b - c = (a - cb) / b */
4231 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4232 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4233 SCM_FRACTION_DENOMINATOR (x
));
4234 else if (SCM_BIGP (y
))
4235 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4236 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4237 SCM_FRACTION_DENOMINATOR (x
));
4238 else if (SCM_REALP (y
))
4239 return scm_make_real (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4240 else if (SCM_COMPLEXP (y
))
4241 return scm_make_complex (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4242 -SCM_COMPLEX_IMAG (y
));
4243 else if (SCM_FRACTIONP (y
))
4244 /* a/b - c/d = (ad - bc) / bd */
4245 return scm_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4246 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4247 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4249 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4252 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4257 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4258 /* "Return the product of all arguments. If called without arguments,\n"
4262 scm_product (SCM x
, SCM y
)
4267 return SCM_I_MAKINUM (1L);
4268 else if (SCM_NUMBERP (x
))
4271 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4274 if (SCM_I_INUMP (x
))
4279 xx
= SCM_I_INUM (x
);
4283 case 0: return x
; break;
4284 case 1: return y
; break;
4287 if (SCM_I_INUMP (y
))
4289 long yy
= SCM_I_INUM (y
);
4291 SCM k
= SCM_I_MAKINUM (kk
);
4292 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4296 SCM result
= scm_i_long2big (xx
);
4297 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4298 return scm_i_normbig (result
);
4301 else if (SCM_BIGP (y
))
4303 SCM result
= scm_i_mkbig ();
4304 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4305 scm_remember_upto_here_1 (y
);
4308 else if (SCM_REALP (y
))
4309 return scm_make_real (xx
* SCM_REAL_VALUE (y
));
4310 else if (SCM_COMPLEXP (y
))
4311 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4312 xx
* SCM_COMPLEX_IMAG (y
));
4313 else if (SCM_FRACTIONP (y
))
4314 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4315 SCM_FRACTION_DENOMINATOR (y
));
4317 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4319 else if (SCM_BIGP (x
))
4321 if (SCM_I_INUMP (y
))
4326 else if (SCM_BIGP (y
))
4328 SCM result
= scm_i_mkbig ();
4329 mpz_mul (SCM_I_BIG_MPZ (result
),
4332 scm_remember_upto_here_2 (x
, y
);
4335 else if (SCM_REALP (y
))
4337 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4338 scm_remember_upto_here_1 (x
);
4339 return scm_make_real (result
);
4341 else if (SCM_COMPLEXP (y
))
4343 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4344 scm_remember_upto_here_1 (x
);
4345 return scm_make_complex (z
* SCM_COMPLEX_REAL (y
),
4346 z
* SCM_COMPLEX_IMAG (y
));
4348 else if (SCM_FRACTIONP (y
))
4349 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4350 SCM_FRACTION_DENOMINATOR (y
));
4352 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4354 else if (SCM_REALP (x
))
4356 if (SCM_I_INUMP (y
))
4357 return scm_make_real (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4358 else if (SCM_BIGP (y
))
4360 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4361 scm_remember_upto_here_1 (y
);
4362 return scm_make_real (result
);
4364 else if (SCM_REALP (y
))
4365 return scm_make_real (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4366 else if (SCM_COMPLEXP (y
))
4367 return scm_make_complex (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4368 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4369 else if (SCM_FRACTIONP (y
))
4370 return scm_make_real (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4372 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4374 else if (SCM_COMPLEXP (x
))
4376 if (SCM_I_INUMP (y
))
4377 return scm_make_complex (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4378 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4379 else if (SCM_BIGP (y
))
4381 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4382 scm_remember_upto_here_1 (y
);
4383 return scm_make_complex (z
* SCM_COMPLEX_REAL (x
),
4384 z
* SCM_COMPLEX_IMAG (x
));
4386 else if (SCM_REALP (y
))
4387 return scm_make_complex (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4388 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4389 else if (SCM_COMPLEXP (y
))
4391 return scm_make_complex (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4392 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4393 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4394 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4396 else if (SCM_FRACTIONP (y
))
4398 double yy
= scm_i_fraction2double (y
);
4399 return scm_make_complex (yy
* SCM_COMPLEX_REAL (x
),
4400 yy
* SCM_COMPLEX_IMAG (x
));
4403 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4405 else if (SCM_FRACTIONP (x
))
4407 if (SCM_I_INUMP (y
))
4408 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4409 SCM_FRACTION_DENOMINATOR (x
));
4410 else if (SCM_BIGP (y
))
4411 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4412 SCM_FRACTION_DENOMINATOR (x
));
4413 else if (SCM_REALP (y
))
4414 return scm_make_real (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4415 else if (SCM_COMPLEXP (y
))
4417 double xx
= scm_i_fraction2double (x
);
4418 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4419 xx
* SCM_COMPLEX_IMAG (y
));
4421 else if (SCM_FRACTIONP (y
))
4422 /* a/b * c/d = ac / bd */
4423 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4424 SCM_FRACTION_NUMERATOR (y
)),
4425 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4426 SCM_FRACTION_DENOMINATOR (y
)));
4428 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4431 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4435 scm_num2dbl (SCM a
, const char *why
)
4436 #define FUNC_NAME why
4438 if (SCM_I_INUMP (a
))
4439 return (double) SCM_I_INUM (a
);
4440 else if (SCM_BIGP (a
))
4442 double result
= mpz_get_d (SCM_I_BIG_MPZ (a
));
4443 scm_remember_upto_here_1 (a
);
4446 else if (SCM_REALP (a
))
4447 return (SCM_REAL_VALUE (a
));
4448 else if (SCM_FRACTIONP (a
))
4449 return scm_i_fraction2double (a
);
4451 SCM_WRONG_TYPE_ARG (SCM_ARGn
, a
);
4455 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4456 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4457 #define ALLOW_DIVIDE_BY_ZERO
4458 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4461 /* The code below for complex division is adapted from the GNU
4462 libstdc++, which adapted it from f2c's libF77, and is subject to
4465 /****************************************************************
4466 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4468 Permission to use, copy, modify, and distribute this software
4469 and its documentation for any purpose and without fee is hereby
4470 granted, provided that the above copyright notice appear in all
4471 copies and that both that the copyright notice and this
4472 permission notice and warranty disclaimer appear in supporting
4473 documentation, and that the names of AT&T Bell Laboratories or
4474 Bellcore or any of their entities not be used in advertising or
4475 publicity pertaining to distribution of the software without
4476 specific, written prior permission.
4478 AT&T and Bellcore disclaim all warranties with regard to this
4479 software, including all implied warranties of merchantability
4480 and fitness. In no event shall AT&T or Bellcore be liable for
4481 any special, indirect or consequential damages or any damages
4482 whatsoever resulting from loss of use, data or profits, whether
4483 in an action of contract, negligence or other tortious action,
4484 arising out of or in connection with the use or performance of
4486 ****************************************************************/
4488 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4489 /* Divide the first argument by the product of the remaining
4490 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4492 #define FUNC_NAME s_divide
4494 scm_i_divide (SCM x
, SCM y
, int inexact
)
4501 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4502 else if (SCM_I_INUMP (x
))
4504 long xx
= SCM_I_INUM (x
);
4505 if (xx
== 1 || xx
== -1)
4507 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4509 scm_num_overflow (s_divide
);
4514 return scm_make_real (1.0 / (double) xx
);
4515 else return scm_make_ratio (SCM_I_MAKINUM(1), x
);
4518 else if (SCM_BIGP (x
))
4521 return scm_make_real (1.0 / scm_i_big2dbl (x
));
4522 else return scm_make_ratio (SCM_I_MAKINUM(1), x
);
4524 else if (SCM_REALP (x
))
4526 double xx
= SCM_REAL_VALUE (x
);
4527 #ifndef ALLOW_DIVIDE_BY_ZERO
4529 scm_num_overflow (s_divide
);
4532 return scm_make_real (1.0 / xx
);
4534 else if (SCM_COMPLEXP (x
))
4536 double r
= SCM_COMPLEX_REAL (x
);
4537 double i
= SCM_COMPLEX_IMAG (x
);
4541 double d
= i
* (1.0 + t
* t
);
4542 return scm_make_complex (t
/ d
, -1.0 / d
);
4547 double d
= r
* (1.0 + t
* t
);
4548 return scm_make_complex (1.0 / d
, -t
/ d
);
4551 else if (SCM_FRACTIONP (x
))
4552 return scm_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4553 SCM_FRACTION_NUMERATOR (x
));
4555 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4558 if (SCM_I_INUMP (x
))
4560 long xx
= SCM_I_INUM (x
);
4561 if (SCM_I_INUMP (y
))
4563 long yy
= SCM_I_INUM (y
);
4566 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4567 scm_num_overflow (s_divide
);
4569 return scm_make_real ((double) xx
/ (double) yy
);
4572 else if (xx
% yy
!= 0)
4575 return scm_make_real ((double) xx
/ (double) yy
);
4576 else return scm_make_ratio (x
, y
);
4581 if (SCM_FIXABLE (z
))
4582 return SCM_I_MAKINUM (z
);
4584 return scm_i_long2big (z
);
4587 else if (SCM_BIGP (y
))
4590 return scm_make_real ((double) xx
/ scm_i_big2dbl (y
));
4591 else return scm_make_ratio (x
, y
);
4593 else if (SCM_REALP (y
))
4595 double yy
= SCM_REAL_VALUE (y
);
4596 #ifndef ALLOW_DIVIDE_BY_ZERO
4598 scm_num_overflow (s_divide
);
4601 return scm_make_real ((double) xx
/ yy
);
4603 else if (SCM_COMPLEXP (y
))
4606 complex_div
: /* y _must_ be a complex number */
4608 double r
= SCM_COMPLEX_REAL (y
);
4609 double i
= SCM_COMPLEX_IMAG (y
);
4613 double d
= i
* (1.0 + t
* t
);
4614 return scm_make_complex ((a
* t
) / d
, -a
/ d
);
4619 double d
= r
* (1.0 + t
* t
);
4620 return scm_make_complex (a
/ d
, -(a
* t
) / d
);
4624 else if (SCM_FRACTIONP (y
))
4625 /* a / b/c = ac / b */
4626 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4627 SCM_FRACTION_NUMERATOR (y
));
4629 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4631 else if (SCM_BIGP (x
))
4633 if (SCM_I_INUMP (y
))
4635 long int yy
= SCM_I_INUM (y
);
4638 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4639 scm_num_overflow (s_divide
);
4641 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4642 scm_remember_upto_here_1 (x
);
4643 return (sgn
== 0) ? scm_nan () : scm_inf ();
4650 /* FIXME: HMM, what are the relative performance issues here?
4651 We need to test. Is it faster on average to test
4652 divisible_p, then perform whichever operation, or is it
4653 faster to perform the integer div opportunistically and
4654 switch to real if there's a remainder? For now we take the
4655 middle ground: test, then if divisible, use the faster div
4658 long abs_yy
= yy
< 0 ? -yy
: yy
;
4659 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4663 SCM result
= scm_i_mkbig ();
4664 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4665 scm_remember_upto_here_1 (x
);
4667 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4668 return scm_i_normbig (result
);
4673 return scm_make_real (scm_i_big2dbl (x
) / (double) yy
);
4674 else return scm_make_ratio (x
, y
);
4678 else if (SCM_BIGP (y
))
4680 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4683 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4684 scm_num_overflow (s_divide
);
4686 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4687 scm_remember_upto_here_1 (x
);
4688 return (sgn
== 0) ? scm_nan () : scm_inf ();
4694 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4698 SCM result
= scm_i_mkbig ();
4699 mpz_divexact (SCM_I_BIG_MPZ (result
),
4702 scm_remember_upto_here_2 (x
, y
);
4703 return scm_i_normbig (result
);
4709 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4710 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4711 scm_remember_upto_here_2 (x
, y
);
4712 return scm_make_real (dbx
/ dby
);
4714 else return scm_make_ratio (x
, y
);
4718 else if (SCM_REALP (y
))
4720 double yy
= SCM_REAL_VALUE (y
);
4721 #ifndef ALLOW_DIVIDE_BY_ZERO
4723 scm_num_overflow (s_divide
);
4726 return scm_make_real (scm_i_big2dbl (x
) / yy
);
4728 else if (SCM_COMPLEXP (y
))
4730 a
= scm_i_big2dbl (x
);
4733 else if (SCM_FRACTIONP (y
))
4734 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4735 SCM_FRACTION_NUMERATOR (y
));
4737 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4739 else if (SCM_REALP (x
))
4741 double rx
= SCM_REAL_VALUE (x
);
4742 if (SCM_I_INUMP (y
))
4744 long int yy
= SCM_I_INUM (y
);
4745 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4747 scm_num_overflow (s_divide
);
4750 return scm_make_real (rx
/ (double) yy
);
4752 else if (SCM_BIGP (y
))
4754 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4755 scm_remember_upto_here_1 (y
);
4756 return scm_make_real (rx
/ dby
);
4758 else if (SCM_REALP (y
))
4760 double yy
= SCM_REAL_VALUE (y
);
4761 #ifndef ALLOW_DIVIDE_BY_ZERO
4763 scm_num_overflow (s_divide
);
4766 return scm_make_real (rx
/ yy
);
4768 else if (SCM_COMPLEXP (y
))
4773 else if (SCM_FRACTIONP (y
))
4774 return scm_make_real (rx
/ scm_i_fraction2double (y
));
4776 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4778 else if (SCM_COMPLEXP (x
))
4780 double rx
= SCM_COMPLEX_REAL (x
);
4781 double ix
= SCM_COMPLEX_IMAG (x
);
4782 if (SCM_I_INUMP (y
))
4784 long int yy
= SCM_I_INUM (y
);
4785 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4787 scm_num_overflow (s_divide
);
4792 return scm_make_complex (rx
/ d
, ix
/ d
);
4795 else if (SCM_BIGP (y
))
4797 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4798 scm_remember_upto_here_1 (y
);
4799 return scm_make_complex (rx
/ dby
, ix
/ dby
);
4801 else if (SCM_REALP (y
))
4803 double yy
= SCM_REAL_VALUE (y
);
4804 #ifndef ALLOW_DIVIDE_BY_ZERO
4806 scm_num_overflow (s_divide
);
4809 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4811 else if (SCM_COMPLEXP (y
))
4813 double ry
= SCM_COMPLEX_REAL (y
);
4814 double iy
= SCM_COMPLEX_IMAG (y
);
4818 double d
= iy
* (1.0 + t
* t
);
4819 return scm_make_complex ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4824 double d
= ry
* (1.0 + t
* t
);
4825 return scm_make_complex ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4828 else if (SCM_FRACTIONP (y
))
4830 double yy
= scm_i_fraction2double (y
);
4831 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4834 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4836 else if (SCM_FRACTIONP (x
))
4838 if (SCM_I_INUMP (y
))
4840 long int yy
= SCM_I_INUM (y
);
4841 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4843 scm_num_overflow (s_divide
);
4846 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4847 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4849 else if (SCM_BIGP (y
))
4851 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4852 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4854 else if (SCM_REALP (y
))
4856 double yy
= SCM_REAL_VALUE (y
);
4857 #ifndef ALLOW_DIVIDE_BY_ZERO
4859 scm_num_overflow (s_divide
);
4862 return scm_make_real (scm_i_fraction2double (x
) / yy
);
4864 else if (SCM_COMPLEXP (y
))
4866 a
= scm_i_fraction2double (x
);
4869 else if (SCM_FRACTIONP (y
))
4870 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4871 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4873 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4876 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4880 scm_divide (SCM x
, SCM y
)
4882 return scm_i_divide (x
, y
, 0);
4885 static SCM
scm_divide2real (SCM x
, SCM y
)
4887 return scm_i_divide (x
, y
, 1);
4893 scm_asinh (double x
)
4898 #define asinh scm_asinh
4899 return log (x
+ sqrt (x
* x
+ 1));
4902 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4903 /* "Return the inverse hyperbolic sine of @var{x}."
4908 scm_acosh (double x
)
4913 #define acosh scm_acosh
4914 return log (x
+ sqrt (x
* x
- 1));
4917 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4918 /* "Return the inverse hyperbolic cosine of @var{x}."
4923 scm_atanh (double x
)
4928 #define atanh scm_atanh
4929 return 0.5 * log ((1 + x
) / (1 - x
));
4932 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4933 /* "Return the inverse hyperbolic tangent of @var{x}."
4937 /* XXX - eventually, we should remove this definition of scm_round and
4938 rename scm_round_number to scm_round. Likewise for scm_truncate
4939 and scm_truncate_number.
4943 scm_truncate (double x
)
4954 /* scm_round is done using floor(x+0.5) to round to nearest and with
4955 half-way case (ie. when x is an integer plus 0.5) going upwards. Then
4956 half-way cases are identified and adjusted down if the round-upwards
4957 didn't give the desired even integer.
4959 "plus_half == result" identifies a half-way case. If plus_half, which is
4960 x + 0.5, is an integer then x must be an integer plus 0.5.
4962 An odd "result" value is identified with result/2 != floor(result/2).
4963 This is done with plus_half, since that value is ready for use sooner in
4964 a pipelined cpu, and we're already requiring plus_half == result.
4966 Note however that we need to be careful when x is big and already an
4967 integer. In that case "x+0.5" may round to an adjacent integer, causing
4968 us to return such a value, incorrectly. For instance if the hardware is
4969 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
4970 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
4971 returned. Or if the hardware is in round-upwards mode, then other bigger
4972 values like say x == 2^128 will see x+0.5 rounding up to the next higher
4973 representable value, 2^128+2^76 (or whatever), again incorrect.
4975 These bad roundings of x+0.5 are avoided by testing at the start whether
4976 x is already an integer. If it is then clearly that's the desired result
4977 already. And if it's not then the exponent must be small enough to allow
4978 an 0.5 to be represented, and hence added without a bad rounding. */
4981 scm_round (double x
)
4983 double plus_half
, result
;
4988 plus_half
= x
+ 0.5;
4989 result
= floor (plus_half
);
4990 /* Adjust so that the scm_round is towards even. */
4991 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4996 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4998 "Round the number @var{x} towards zero.")
4999 #define FUNC_NAME s_scm_truncate_number
5001 if (scm_is_false (scm_negative_p (x
)))
5002 return scm_floor (x
);
5004 return scm_ceiling (x
);
5008 static SCM exactly_one_half
;
5010 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5012 "Round the number @var{x} towards the nearest integer. "
5013 "When it is exactly halfway between two integers, "
5014 "round towards the even one.")
5015 #define FUNC_NAME s_scm_round_number
5017 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5019 else if (SCM_REALP (x
))
5020 return scm_make_real (scm_round (SCM_REAL_VALUE (x
)));
5023 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5024 single quotient+remainder division then examining to see which way
5025 the rounding should go. */
5026 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5027 SCM result
= scm_floor (plus_half
);
5028 /* Adjust so that the scm_round is towards even. */
5029 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5030 && scm_is_true (scm_odd_p (result
)))
5031 return scm_difference (result
, SCM_I_MAKINUM (1));
5038 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5040 "Round the number @var{x} towards minus infinity.")
5041 #define FUNC_NAME s_scm_floor
5043 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5045 else if (SCM_REALP (x
))
5046 return scm_make_real (floor (SCM_REAL_VALUE (x
)));
5047 else if (SCM_FRACTIONP (x
))
5049 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5050 SCM_FRACTION_DENOMINATOR (x
));
5051 if (scm_is_false (scm_negative_p (x
)))
5053 /* For positive x, rounding towards zero is correct. */
5058 /* For negative x, we need to return q-1 unless x is an
5059 integer. But fractions are never integer, per our
5061 return scm_difference (q
, SCM_I_MAKINUM (1));
5065 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5069 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5071 "Round the number @var{x} towards infinity.")
5072 #define FUNC_NAME s_scm_ceiling
5074 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5076 else if (SCM_REALP (x
))
5077 return scm_make_real (ceil (SCM_REAL_VALUE (x
)));
5078 else if (SCM_FRACTIONP (x
))
5080 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5081 SCM_FRACTION_DENOMINATOR (x
));
5082 if (scm_is_false (scm_positive_p (x
)))
5084 /* For negative x, rounding towards zero is correct. */
5089 /* For positive x, we need to return q+1 unless x is an
5090 integer. But fractions are never integer, per our
5092 return scm_sum (q
, SCM_I_MAKINUM (1));
5096 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5100 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5101 /* "Return the square root of the real number @var{x}."
5103 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5104 /* "Return the absolute value of the real number @var{x}."
5106 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5107 /* "Return the @var{x}th power of e."
5109 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5110 /* "Return the natural logarithm of the real number @var{x}."
5112 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5113 /* "Return the sine of the real number @var{x}."
5115 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5116 /* "Return the cosine of the real number @var{x}."
5118 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5119 /* "Return the tangent of the real number @var{x}."
5121 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5122 /* "Return the arc sine of the real number @var{x}."
5124 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5125 /* "Return the arc cosine of the real number @var{x}."
5127 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5128 /* "Return the arc tangent of the real number @var{x}."
5130 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5131 /* "Return the hyperbolic sine of the real number @var{x}."
5133 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5134 /* "Return the hyperbolic cosine of the real number @var{x}."
5136 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5137 /* "Return the hyperbolic tangent of the real number @var{x}."
5145 static void scm_two_doubles (SCM x
,
5147 const char *sstring
,
5151 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5153 if (SCM_I_INUMP (x
))
5154 xy
->x
= SCM_I_INUM (x
);
5155 else if (SCM_BIGP (x
))
5156 xy
->x
= scm_i_big2dbl (x
);
5157 else if (SCM_REALP (x
))
5158 xy
->x
= SCM_REAL_VALUE (x
);
5159 else if (SCM_FRACTIONP (x
))
5160 xy
->x
= scm_i_fraction2double (x
);
5162 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5164 if (SCM_I_INUMP (y
))
5165 xy
->y
= SCM_I_INUM (y
);
5166 else if (SCM_BIGP (y
))
5167 xy
->y
= scm_i_big2dbl (y
);
5168 else if (SCM_REALP (y
))
5169 xy
->y
= SCM_REAL_VALUE (y
);
5170 else if (SCM_FRACTIONP (y
))
5171 xy
->y
= scm_i_fraction2double (y
);
5173 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5177 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5179 "Return @var{x} raised to the power of @var{y}. This\n"
5180 "procedure does not accept complex arguments.")
5181 #define FUNC_NAME s_scm_sys_expt
5184 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5185 return scm_make_real (pow (xy
.x
, xy
.y
));
5190 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5192 "Return the arc tangent of the two arguments @var{x} and\n"
5193 "@var{y}. This is similar to calculating the arc tangent of\n"
5194 "@var{x} / @var{y}, except that the signs of both arguments\n"
5195 "are used to determine the quadrant of the result. This\n"
5196 "procedure does not accept complex arguments.")
5197 #define FUNC_NAME s_scm_sys_atan2
5200 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5201 return scm_make_real (atan2 (xy
.x
, xy
.y
));
5206 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5207 (SCM real
, SCM imaginary
),
5208 "Return a complex number constructed of the given @var{real} and\n"
5209 "@var{imaginary} parts.")
5210 #define FUNC_NAME s_scm_make_rectangular
5213 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5214 return scm_make_complex (xy
.x
, xy
.y
);
5220 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5222 "Return the complex number @var{x} * e^(i * @var{y}).")
5223 #define FUNC_NAME s_scm_make_polar
5227 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5229 sincos (xy
.y
, &s
, &c
);
5234 return scm_make_complex (xy
.x
* c
, xy
.x
* s
);
5239 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5240 /* "Return the real part of the number @var{z}."
5243 scm_real_part (SCM z
)
5245 if (SCM_I_INUMP (z
))
5247 else if (SCM_BIGP (z
))
5249 else if (SCM_REALP (z
))
5251 else if (SCM_COMPLEXP (z
))
5252 return scm_make_real (SCM_COMPLEX_REAL (z
));
5253 else if (SCM_FRACTIONP (z
))
5256 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5260 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5261 /* "Return the imaginary part of the number @var{z}."
5264 scm_imag_part (SCM z
)
5266 if (SCM_I_INUMP (z
))
5268 else if (SCM_BIGP (z
))
5270 else if (SCM_REALP (z
))
5272 else if (SCM_COMPLEXP (z
))
5273 return scm_make_real (SCM_COMPLEX_IMAG (z
));
5274 else if (SCM_FRACTIONP (z
))
5277 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5280 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5281 /* "Return the numerator of the number @var{z}."
5284 scm_numerator (SCM z
)
5286 if (SCM_I_INUMP (z
))
5288 else if (SCM_BIGP (z
))
5290 else if (SCM_FRACTIONP (z
))
5292 scm_i_fraction_reduce (z
);
5293 return SCM_FRACTION_NUMERATOR (z
);
5295 else if (SCM_REALP (z
))
5296 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5298 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5302 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5303 /* "Return the denominator of the number @var{z}."
5306 scm_denominator (SCM z
)
5308 if (SCM_I_INUMP (z
))
5309 return SCM_I_MAKINUM (1);
5310 else if (SCM_BIGP (z
))
5311 return SCM_I_MAKINUM (1);
5312 else if (SCM_FRACTIONP (z
))
5314 scm_i_fraction_reduce (z
);
5315 return SCM_FRACTION_DENOMINATOR (z
);
5317 else if (SCM_REALP (z
))
5318 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5320 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5323 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5324 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5325 * "@code{abs} for real arguments, but also allows complex numbers."
5328 scm_magnitude (SCM z
)
5330 if (SCM_I_INUMP (z
))
5332 long int zz
= SCM_I_INUM (z
);
5335 else if (SCM_POSFIXABLE (-zz
))
5336 return SCM_I_MAKINUM (-zz
);
5338 return scm_i_long2big (-zz
);
5340 else if (SCM_BIGP (z
))
5342 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5343 scm_remember_upto_here_1 (z
);
5345 return scm_i_clonebig (z
, 0);
5349 else if (SCM_REALP (z
))
5350 return scm_make_real (fabs (SCM_REAL_VALUE (z
)));
5351 else if (SCM_COMPLEXP (z
))
5352 return scm_make_real (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5353 else if (SCM_FRACTIONP (z
))
5355 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5357 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5358 SCM_FRACTION_DENOMINATOR (z
));
5361 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5365 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5366 /* "Return the angle of the complex number @var{z}."
5371 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5372 scm_flo0 to save allocating a new flonum with scm_make_real each time.
5373 But if atan2 follows the floating point rounding mode, then the value
5374 is not a constant. Maybe it'd be close enough though. */
5375 if (SCM_I_INUMP (z
))
5377 if (SCM_I_INUM (z
) >= 0)
5380 return scm_make_real (atan2 (0.0, -1.0));
5382 else if (SCM_BIGP (z
))
5384 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5385 scm_remember_upto_here_1 (z
);
5387 return scm_make_real (atan2 (0.0, -1.0));
5391 else if (SCM_REALP (z
))
5393 if (SCM_REAL_VALUE (z
) >= 0)
5396 return scm_make_real (atan2 (0.0, -1.0));
5398 else if (SCM_COMPLEXP (z
))
5399 return scm_make_real (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5400 else if (SCM_FRACTIONP (z
))
5402 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5404 else return scm_make_real (atan2 (0.0, -1.0));
5407 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5411 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5412 /* Convert the number @var{x} to its inexact representation.\n"
5415 scm_exact_to_inexact (SCM z
)
5417 if (SCM_I_INUMP (z
))
5418 return scm_make_real ((double) SCM_I_INUM (z
));
5419 else if (SCM_BIGP (z
))
5420 return scm_make_real (scm_i_big2dbl (z
));
5421 else if (SCM_FRACTIONP (z
))
5422 return scm_make_real (scm_i_fraction2double (z
));
5423 else if (SCM_INEXACTP (z
))
5426 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5430 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5432 "Return an exact number that is numerically closest to @var{z}.")
5433 #define FUNC_NAME s_scm_inexact_to_exact
5435 if (SCM_I_INUMP (z
))
5437 else if (SCM_BIGP (z
))
5439 else if (SCM_REALP (z
))
5441 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5442 SCM_OUT_OF_RANGE (1, z
);
5449 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5450 q
= scm_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5451 scm_i_mpz2num (mpq_denref (frac
)));
5453 /* When scm_make_ratio throws, we leak the memory allocated
5460 else if (SCM_FRACTIONP (z
))
5463 SCM_WRONG_TYPE_ARG (1, z
);
5467 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5469 "Return an exact number that is within @var{err} of @var{x}.")
5470 #define FUNC_NAME s_scm_rationalize
5472 if (SCM_I_INUMP (x
))
5474 else if (SCM_BIGP (x
))
5476 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5478 /* Use continued fractions to find closest ratio. All
5479 arithmetic is done with exact numbers.
5482 SCM ex
= scm_inexact_to_exact (x
);
5483 SCM int_part
= scm_floor (ex
);
5484 SCM tt
= SCM_I_MAKINUM (1);
5485 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5486 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5490 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5493 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5494 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5496 /* We stop after a million iterations just to be absolutely sure
5497 that we don't go into an infinite loop. The process normally
5498 converges after less than a dozen iterations.
5501 err
= scm_abs (err
);
5502 while (++i
< 1000000)
5504 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5505 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5506 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5508 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5509 err
))) /* abs(x-a/b) <= err */
5511 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5512 if (scm_is_false (scm_exact_p (x
))
5513 || scm_is_false (scm_exact_p (err
)))
5514 return scm_exact_to_inexact (res
);
5518 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5520 tt
= scm_floor (rx
); /* tt = floor (rx) */
5526 scm_num_overflow (s_scm_rationalize
);
5529 SCM_WRONG_TYPE_ARG (1, x
);
5533 #define NUM2FLOAT scm_num2float
5534 #define FLOAT2NUM scm_float2num
5536 #include "libguile/num2float.i.c"
5538 #define NUM2FLOAT scm_num2double
5539 #define FLOAT2NUM scm_double2num
5540 #define FTYPE double
5541 #include "libguile/num2float.i.c"
5543 /* conversion functions */
5546 scm_is_integer (SCM val
)
5548 return scm_is_true (scm_integer_p (val
));
5552 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5554 if (SCM_I_INUMP (val
))
5556 scm_t_signed_bits n
= SCM_I_INUM (val
);
5557 return n
>= min
&& n
<= max
;
5559 else if (SCM_BIGP (val
))
5561 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5563 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5565 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5567 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5568 return n
>= min
&& n
<= max
;
5578 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5579 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5582 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5583 SCM_I_BIG_MPZ (val
));
5585 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5597 return n
>= min
&& n
<= max
;
5605 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5607 if (SCM_I_INUMP (val
))
5609 scm_t_signed_bits n
= SCM_I_INUM (val
);
5610 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5612 else if (SCM_BIGP (val
))
5614 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5616 else if (max
<= ULONG_MAX
)
5618 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5620 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5621 return n
>= min
&& n
<= max
;
5631 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5634 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5635 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5638 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5639 SCM_I_BIG_MPZ (val
));
5641 return n
>= min
&& n
<= max
;
5648 #define TYPE scm_t_intmax
5649 #define TYPE_MIN min
5650 #define TYPE_MAX max
5651 #define SIZEOF_TYPE 0
5652 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5653 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5654 #include "libguile/conv-integer.i.c"
5656 #define TYPE scm_t_uintmax
5657 #define TYPE_MIN min
5658 #define TYPE_MAX max
5659 #define SIZEOF_TYPE 0
5660 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5661 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5662 #include "libguile/conv-uinteger.i.c"
5664 #define TYPE scm_t_int8
5665 #define TYPE_MIN SCM_T_INT8_MIN
5666 #define TYPE_MAX SCM_T_INT8_MAX
5667 #define SIZEOF_TYPE 1
5668 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5669 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5670 #include "libguile/conv-integer.i.c"
5672 #define TYPE scm_t_uint8
5674 #define TYPE_MAX SCM_T_UINT8_MAX
5675 #define SIZEOF_TYPE 1
5676 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5677 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5678 #include "libguile/conv-uinteger.i.c"
5680 #define TYPE scm_t_int16
5681 #define TYPE_MIN SCM_T_INT16_MIN
5682 #define TYPE_MAX SCM_T_INT16_MAX
5683 #define SIZEOF_TYPE 2
5684 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5685 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5686 #include "libguile/conv-integer.i.c"
5688 #define TYPE scm_t_uint16
5690 #define TYPE_MAX SCM_T_UINT16_MAX
5691 #define SIZEOF_TYPE 2
5692 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5693 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5694 #include "libguile/conv-uinteger.i.c"
5696 #define TYPE scm_t_int32
5697 #define TYPE_MIN SCM_T_INT32_MIN
5698 #define TYPE_MAX SCM_T_INT32_MAX
5699 #define SIZEOF_TYPE 4
5700 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5701 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5702 #include "libguile/conv-integer.i.c"
5704 #define TYPE scm_t_uint32
5706 #define TYPE_MAX SCM_T_UINT32_MAX
5707 #define SIZEOF_TYPE 4
5708 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5709 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5710 #include "libguile/conv-uinteger.i.c"
5712 #if SCM_HAVE_T_INT64
5714 #define TYPE scm_t_int64
5715 #define TYPE_MIN SCM_T_INT64_MIN
5716 #define TYPE_MAX SCM_T_INT64_MAX
5717 #define SIZEOF_TYPE 8
5718 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5719 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5720 #include "libguile/conv-integer.i.c"
5722 #define TYPE scm_t_uint64
5724 #define TYPE_MAX SCM_T_UINT64_MAX
5725 #define SIZEOF_TYPE 8
5726 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5727 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5728 #include "libguile/conv-uinteger.i.c"
5733 scm_is_real (SCM val
)
5735 return scm_is_true (scm_real_p (val
));
5739 scm_to_double (SCM val
)
5741 return scm_num2dbl (val
, NULL
);
5745 scm_from_double (double val
)
5747 return scm_make_real (val
);
5755 mpz_init_set_si (z_negative_one
, -1);
5757 /* It may be possible to tune the performance of some algorithms by using
5758 * the following constants to avoid the creation of bignums. Please, before
5759 * using these values, remember the two rules of program optimization:
5760 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5761 scm_c_define ("most-positive-fixnum",
5762 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5763 scm_c_define ("most-negative-fixnum",
5764 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5766 scm_add_feature ("complex");
5767 scm_add_feature ("inexact");
5768 scm_flo0
= scm_make_real (0.0);
5770 /* determine floating point precision */
5771 for(i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
5773 init_dblprec(&scm_dblprec
[i
-2],i
);
5774 init_fx_radix(fx_per_radix
[i
-2],i
);
5777 /* hard code precision for base 10 if the preprocessor tells us to... */
5778 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5785 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
5786 SCM_I_MAKINUM (2)));
5787 #include "libguile/numbers.x"