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_INUMP(x) ? SCM_I_NUMTAG_INUM \
82 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
83 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
84 : SCM_I_NUMTAG_NOTNUM)))
86 /* the macro above will not work as is with fractions */
89 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
91 /* FLOBUFLEN is the maximum number of characters neccessary for the
92 * printed or scm_string representation of an inexact number.
94 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
97 #if ! defined (HAVE_ISNAN)
102 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
105 #if ! defined (HAVE_ISINF)
110 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
117 /* mpz_cmp_d only recognises infinities in gmp 4.2 and up.
118 For prior versions use an explicit check here. */
119 #if __GNU_MP_VERSION < 4 \
120 || (__GNU_MP_VERSION == 4 && __GNU_MP_VERSION_MINOR < 2)
121 #define xmpz_cmp_d(z, d) \
122 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
124 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
127 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
128 isinf. It does have finite and isnan though, hence the use of those.
129 fpclass would be a possibility on that system too. */
133 #if defined (HAVE_ISINF)
135 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
136 return (! (finite (x
) || isnan (x
)));
145 #if defined (HAVE_ISNAN)
154 static mpz_t z_negative_one
;
158 SCM_C_INLINE_KEYWORD SCM
161 /* Return a newly created bignum. */
162 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
163 mpz_init (SCM_I_BIG_MPZ (z
));
167 SCM_C_INLINE_KEYWORD
static SCM
168 scm_i_clonebig (SCM src_big
, int same_sign_p
)
170 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
171 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
172 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
174 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
178 SCM_C_INLINE_KEYWORD
int
179 scm_i_bigcmp (SCM x
, SCM y
)
181 /* Return neg if x < y, pos if x > y, and 0 if x == y */
182 /* presume we already know x and y are bignums */
183 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
184 scm_remember_upto_here_2 (x
, y
);
188 SCM_C_INLINE_KEYWORD SCM
189 scm_i_dbl2big (double d
)
191 /* results are only defined if d is an integer */
192 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
193 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
197 /* Convert a integer in double representation to a SCM number. */
199 SCM_C_INLINE_KEYWORD SCM
200 scm_i_dbl2num (double u
)
202 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
203 powers of 2, so there's no rounding when making "double" values
204 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
205 get rounded on a 64-bit machine, hence the "+1".
207 The use of floor() to force to an integer value ensures we get a
208 "numerically closest" value without depending on how a
209 double->long cast or how mpz_set_d will round. For reference,
210 double->long probably follows the hardware rounding mode,
211 mpz_set_d truncates towards zero. */
213 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
214 representable as a double? */
216 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
217 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
218 return SCM_MAKINUM ((long) u
);
220 return scm_i_dbl2big (u
);
223 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
224 with R5RS exact->inexact.
226 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
227 (ie. it truncates towards zero), then adjust to get the closest double by
228 examining the next lower bit and adding 1 if necessary.
230 Note that bignums exactly half way between representable doubles are
231 rounded to the next higher absolute value (ie. away from zero). This
232 seems like an adequate interpretation of R5RS "numerically closest", and
233 it's easier and faster than a full "nearest-even" style.
235 The bit test is done on the absolute value of the mpz_t, which means we
236 must use mpz_getlimbn. mpz_tstbit is not right, it treats negatives as
239 Prior to GMP 4.2, the rounding done by mpz_get_d was unspecified. It
240 happened to follow the hardware rounding mode, but on the absolute value
241 of its operand. This is not what we want, so we put the high
242 DBL_MANT_DIG bits into a temporary. This extra init/clear is a slowdown,
243 but doesn't matter too much since it's only for older GMP. */
246 scm_i_big2dbl (SCM b
)
251 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
253 #if __GNU_MP_VERSION < 4 \
254 || (__GNU_MP_VERSION == 4 && __GNU_MP_VERSION_MINOR < 2)
256 /* GMP prior to 4.2, force truncate towards zero */
258 if (bits
> DBL_MANT_DIG
)
260 size_t shift
= bits
- DBL_MANT_DIG
;
261 mpz_init2 (tmp
, DBL_MANT_DIG
);
262 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
263 result
= ldexp (mpz_get_d (tmp
), shift
);
268 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
273 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
276 if (bits
> DBL_MANT_DIG
)
278 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
279 /* test bit number "pos" in absolute value */
280 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
281 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
283 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
287 scm_remember_upto_here_1 (b
);
291 SCM_C_INLINE_KEYWORD SCM
292 scm_i_normbig (SCM b
)
294 /* convert a big back to a fixnum if it'll fit */
295 /* presume b is a bignum */
296 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
298 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
299 if (SCM_FIXABLE (val
))
300 b
= SCM_MAKINUM (val
);
305 static SCM_C_INLINE_KEYWORD SCM
306 scm_i_mpz2num (mpz_t b
)
308 /* convert a mpz number to a SCM number. */
309 if (mpz_fits_slong_p (b
))
311 long val
= mpz_get_si (b
);
312 if (SCM_FIXABLE (val
))
313 return SCM_MAKINUM (val
);
317 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
318 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
323 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
324 static SCM
scm_divide2real (SCM x
, SCM y
);
327 scm_make_ratio (SCM numerator
, SCM denominator
)
328 #define FUNC_NAME "make-ratio"
330 /* First make sure the arguments are proper.
332 if (SCM_INUMP (denominator
))
334 if (SCM_EQ_P (denominator
, SCM_INUM0
))
335 scm_num_overflow ("make-ratio");
336 if (SCM_EQ_P (denominator
, SCM_MAKINUM(1)))
341 if (!(SCM_BIGP(denominator
)))
342 SCM_WRONG_TYPE_ARG (2, denominator
);
344 if (!SCM_INUMP (numerator
) && !SCM_BIGP (numerator
))
345 SCM_WRONG_TYPE_ARG (1, numerator
);
347 /* Then flip signs so that the denominator is positive.
349 if (SCM_NFALSEP (scm_negative_p (denominator
)))
351 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
352 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
355 /* Now consider for each of the four fixnum/bignum combinations
356 whether the rational number is really an integer.
358 if (SCM_INUMP (numerator
))
360 long x
= SCM_INUM (numerator
);
361 if (SCM_EQ_P (numerator
, SCM_INUM0
))
363 if (SCM_INUMP (denominator
))
366 y
= SCM_INUM (denominator
);
368 return SCM_MAKINUM(1);
370 return SCM_MAKINUM (x
/ y
);
374 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
375 of that value for the denominator, as a bignum. Apart from
376 that case, abs(bignum) > abs(inum) so inum/bignum is not an
378 if (x
== SCM_MOST_NEGATIVE_FIXNUM
379 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
380 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
381 return SCM_MAKINUM(-1);
384 else if (SCM_BIGP (numerator
))
386 if (SCM_INUMP (denominator
))
388 long yy
= SCM_INUM (denominator
);
389 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
390 return scm_divide (numerator
, denominator
);
394 if (SCM_EQ_P (numerator
, denominator
))
395 return SCM_MAKINUM(1);
396 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
397 SCM_I_BIG_MPZ (denominator
)))
398 return scm_divide(numerator
, denominator
);
402 /* No, it's a proper fraction.
404 return scm_double_cell (scm_tc16_fraction
,
405 SCM_UNPACK (numerator
),
406 SCM_UNPACK (denominator
), 0);
410 static void scm_i_fraction_reduce (SCM z
)
412 if (!(SCM_FRACTION_REDUCED (z
)))
415 divisor
= scm_gcd (SCM_FRACTION_NUMERATOR (z
), SCM_FRACTION_DENOMINATOR (z
));
416 if (!(SCM_EQ_P (divisor
, SCM_MAKINUM(1))))
419 SCM_FRACTION_SET_NUMERATOR (z
, scm_divide (SCM_FRACTION_NUMERATOR (z
), divisor
));
420 SCM_FRACTION_SET_DENOMINATOR (z
, scm_divide (SCM_FRACTION_DENOMINATOR (z
), divisor
));
422 SCM_FRACTION_REDUCED_SET (z
);
427 scm_i_fraction2double (SCM z
)
429 return scm_num2dbl (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
430 SCM_FRACTION_DENOMINATOR (z
)),
434 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
436 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
438 #define FUNC_NAME s_scm_exact_p
444 if (SCM_FRACTIONP (x
))
448 SCM_WRONG_TYPE_ARG (1, x
);
453 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
455 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
457 #define FUNC_NAME s_scm_odd_p
461 long val
= SCM_INUM (n
);
462 return SCM_BOOL ((val
& 1L) != 0);
464 else if (SCM_BIGP (n
))
466 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
467 scm_remember_upto_here_1 (n
);
468 return SCM_BOOL (odd_p
);
470 else if (!SCM_FALSEP (scm_inf_p (n
)))
472 else if (SCM_REALP (n
))
474 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
480 SCM_WRONG_TYPE_ARG (1, n
);
483 SCM_WRONG_TYPE_ARG (1, n
);
488 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
490 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
492 #define FUNC_NAME s_scm_even_p
496 long val
= SCM_INUM (n
);
497 return SCM_BOOL ((val
& 1L) == 0);
499 else if (SCM_BIGP (n
))
501 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
502 scm_remember_upto_here_1 (n
);
503 return SCM_BOOL (even_p
);
505 else if (!SCM_FALSEP (scm_inf_p (n
)))
507 else if (SCM_REALP (n
))
509 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
515 SCM_WRONG_TYPE_ARG (1, n
);
518 SCM_WRONG_TYPE_ARG (1, n
);
522 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
524 "Return @code{#t} if @var{n} is infinite, @code{#f}\n"
526 #define FUNC_NAME s_scm_inf_p
529 return SCM_BOOL (xisinf (SCM_REAL_VALUE (n
)));
530 else if (SCM_COMPLEXP (n
))
531 return SCM_BOOL (xisinf (SCM_COMPLEX_REAL (n
))
532 || xisinf (SCM_COMPLEX_IMAG (n
)));
538 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
540 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
542 #define FUNC_NAME s_scm_nan_p
545 return SCM_BOOL (xisnan (SCM_REAL_VALUE (n
)));
546 else if (SCM_COMPLEXP (n
))
547 return SCM_BOOL (xisnan (SCM_COMPLEX_REAL (n
))
548 || xisnan (SCM_COMPLEX_IMAG (n
)));
554 /* Guile's idea of infinity. */
555 static double guile_Inf
;
557 /* Guile's idea of not a number. */
558 static double guile_NaN
;
561 guile_ieee_init (void)
563 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
565 /* Some version of gcc on some old version of Linux used to crash when
566 trying to make Inf and NaN. */
569 /* C99 INFINITY, when available.
570 FIXME: The standard allows for INFINITY to be something that overflows
571 at compile time. We ought to have a configure test to check for that
572 before trying to use it. (But in practice we believe this is not a
573 problem on any system guile is likely to target.) */
574 guile_Inf
= INFINITY
;
577 extern unsigned int DINFINITY
[2];
578 guile_Inf
= (*(X_CAST(double *, DINFINITY
)));
585 if (guile_Inf
== tmp
)
593 #if defined (HAVE_ISNAN)
596 /* C99 NAN, when available */
600 extern unsigned int DQNAN
[2];
601 guile_NaN
= (*(X_CAST(double *, DQNAN
)));
603 guile_NaN
= guile_Inf
/ guile_Inf
;
609 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
612 #define FUNC_NAME s_scm_inf
614 static int initialized
= 0;
620 return scm_make_real (guile_Inf
);
624 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
627 #define FUNC_NAME s_scm_nan
629 static int initialized
= 0;
635 return scm_make_real (guile_NaN
);
640 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
642 "Return the absolute value of @var{x}.")
647 long int xx
= SCM_INUM (x
);
650 else if (SCM_POSFIXABLE (-xx
))
651 return SCM_MAKINUM (-xx
);
653 return scm_i_long2big (-xx
);
655 else if (SCM_BIGP (x
))
657 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
659 return scm_i_clonebig (x
, 0);
663 else if (SCM_REALP (x
))
665 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
666 double xx
= SCM_REAL_VALUE (x
);
668 return scm_make_real (-xx
);
672 else if (SCM_FRACTIONP (x
))
674 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
676 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
677 SCM_FRACTION_DENOMINATOR (x
));
680 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
685 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
686 /* "Return the quotient of the numbers @var{x} and @var{y}."
689 scm_quotient (SCM x
, SCM y
)
693 long xx
= SCM_INUM (x
);
696 long yy
= SCM_INUM (y
);
698 scm_num_overflow (s_quotient
);
703 return SCM_MAKINUM (z
);
705 return scm_i_long2big (z
);
708 else if (SCM_BIGP (y
))
710 if ((SCM_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
711 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
712 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
714 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
715 scm_remember_upto_here_1 (y
);
716 return SCM_MAKINUM (-1);
719 return SCM_MAKINUM (0);
722 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
724 else if (SCM_BIGP (x
))
728 long yy
= SCM_INUM (y
);
730 scm_num_overflow (s_quotient
);
735 SCM result
= scm_i_mkbig ();
738 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
741 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
744 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
745 scm_remember_upto_here_1 (x
);
746 return scm_i_normbig (result
);
749 else if (SCM_BIGP (y
))
751 SCM result
= scm_i_mkbig ();
752 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
755 scm_remember_upto_here_2 (x
, y
);
756 return scm_i_normbig (result
);
759 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
762 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
765 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
766 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
768 * "(remainder 13 4) @result{} 1\n"
769 * "(remainder -13 4) @result{} -1\n"
773 scm_remainder (SCM x
, SCM y
)
779 long yy
= SCM_INUM (y
);
781 scm_num_overflow (s_remainder
);
784 long z
= SCM_INUM (x
) % yy
;
785 return SCM_MAKINUM (z
);
788 else if (SCM_BIGP (y
))
790 if ((SCM_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
791 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
792 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
794 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
795 scm_remember_upto_here_1 (y
);
796 return SCM_MAKINUM (0);
802 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
804 else if (SCM_BIGP (x
))
808 long yy
= SCM_INUM (y
);
810 scm_num_overflow (s_remainder
);
813 SCM result
= scm_i_mkbig ();
816 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
817 scm_remember_upto_here_1 (x
);
818 return scm_i_normbig (result
);
821 else if (SCM_BIGP (y
))
823 SCM result
= scm_i_mkbig ();
824 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
827 scm_remember_upto_here_2 (x
, y
);
828 return scm_i_normbig (result
);
831 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
834 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
838 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
839 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
841 * "(modulo 13 4) @result{} 1\n"
842 * "(modulo -13 4) @result{} 3\n"
846 scm_modulo (SCM x
, SCM y
)
850 long xx
= SCM_INUM (x
);
853 long yy
= SCM_INUM (y
);
855 scm_num_overflow (s_modulo
);
858 /* FIXME: I think this may be a bug on some arches -- results
859 of % with negative second arg are undefined... */
877 return SCM_MAKINUM (result
);
880 else if (SCM_BIGP (y
))
882 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
889 SCM pos_y
= scm_i_clonebig (y
, 0);
890 /* do this after the last scm_op */
891 mpz_init_set_si (z_x
, xx
);
892 result
= pos_y
; /* re-use this bignum */
893 mpz_mod (SCM_I_BIG_MPZ (result
),
895 SCM_I_BIG_MPZ (pos_y
));
896 scm_remember_upto_here_1 (pos_y
);
900 result
= scm_i_mkbig ();
901 /* do this after the last scm_op */
902 mpz_init_set_si (z_x
, xx
);
903 mpz_mod (SCM_I_BIG_MPZ (result
),
906 scm_remember_upto_here_1 (y
);
909 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
910 mpz_add (SCM_I_BIG_MPZ (result
),
912 SCM_I_BIG_MPZ (result
));
913 scm_remember_upto_here_1 (y
);
914 /* and do this before the next one */
916 return scm_i_normbig (result
);
920 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
922 else if (SCM_BIGP (x
))
926 long yy
= SCM_INUM (y
);
928 scm_num_overflow (s_modulo
);
931 SCM result
= scm_i_mkbig ();
932 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
934 (yy
< 0) ? - yy
: yy
);
935 scm_remember_upto_here_1 (x
);
936 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
937 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
938 SCM_I_BIG_MPZ (result
),
940 return scm_i_normbig (result
);
943 else if (SCM_BIGP (y
))
946 SCM result
= scm_i_mkbig ();
947 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
948 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
949 mpz_mod (SCM_I_BIG_MPZ (result
),
951 SCM_I_BIG_MPZ (pos_y
));
953 scm_remember_upto_here_1 (x
);
954 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
955 mpz_add (SCM_I_BIG_MPZ (result
),
957 SCM_I_BIG_MPZ (result
));
958 scm_remember_upto_here_2 (y
, pos_y
);
959 return scm_i_normbig (result
);
963 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
966 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
969 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
970 /* "Return the greatest common divisor of all arguments.\n"
971 * "If called without arguments, 0 is returned."
974 scm_gcd (SCM x
, SCM y
)
977 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
983 long xx
= SCM_INUM (x
);
984 long yy
= SCM_INUM (y
);
985 long u
= xx
< 0 ? -xx
: xx
;
986 long v
= yy
< 0 ? -yy
: yy
;
996 /* Determine a common factor 2^k */
997 while (!(1 & (u
| v
)))
1003 /* Now, any factor 2^n can be eliminated */
1023 return (SCM_POSFIXABLE (result
)
1024 ? SCM_MAKINUM (result
)
1025 : scm_i_long2big (result
));
1027 else if (SCM_BIGP (y
))
1033 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1035 else if (SCM_BIGP (x
))
1039 unsigned long result
;
1047 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1048 scm_remember_upto_here_1 (x
);
1049 return (SCM_POSFIXABLE (result
)
1050 ? SCM_MAKINUM (result
)
1051 : scm_ulong2num (result
));
1053 else if (SCM_BIGP (y
))
1055 SCM result
= scm_i_mkbig ();
1056 mpz_gcd (SCM_I_BIG_MPZ (result
),
1059 scm_remember_upto_here_2 (x
, y
);
1060 return scm_i_normbig (result
);
1063 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1066 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1069 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1070 /* "Return the least common multiple of the arguments.\n"
1071 * "If called without arguments, 1 is returned."
1074 scm_lcm (SCM n1
, SCM n2
)
1076 if (SCM_UNBNDP (n2
))
1078 if (SCM_UNBNDP (n1
))
1079 return SCM_MAKINUM (1L);
1080 n2
= SCM_MAKINUM (1L);
1083 SCM_GASSERT2 (SCM_INUMP (n1
) || SCM_BIGP (n1
),
1084 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1085 SCM_GASSERT2 (SCM_INUMP (n2
) || SCM_BIGP (n2
),
1086 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1092 SCM d
= scm_gcd (n1
, n2
);
1093 if (SCM_EQ_P (d
, SCM_INUM0
))
1096 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1100 /* inum n1, big n2 */
1103 SCM result
= scm_i_mkbig ();
1104 long nn1
= SCM_INUM (n1
);
1105 if (nn1
== 0) return SCM_INUM0
;
1106 if (nn1
< 0) nn1
= - nn1
;
1107 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1108 scm_remember_upto_here_1 (n2
);
1123 SCM result
= scm_i_mkbig ();
1124 mpz_lcm(SCM_I_BIG_MPZ (result
),
1126 SCM_I_BIG_MPZ (n2
));
1127 scm_remember_upto_here_2(n1
, n2
);
1128 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1134 #ifndef scm_long2num
1135 #define SCM_LOGOP_RETURN(x) scm_ulong2num(x)
1137 #define SCM_LOGOP_RETURN(x) SCM_MAKINUM(x)
1140 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1145 + + + x (map digit:logand X Y)
1146 + - + x (map digit:logand X (lognot (+ -1 Y)))
1147 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1148 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1153 + + + (map digit:logior X Y)
1154 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1155 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1156 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1161 + + + (map digit:logxor X Y)
1162 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1163 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1164 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1169 + + (any digit:logand X Y)
1170 + - (any digit:logand X (lognot (+ -1 Y)))
1171 - + (any digit:logand (lognot (+ -1 X)) Y)
1176 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1178 "Return the bitwise AND of the integer arguments.\n\n"
1180 "(logand) @result{} -1\n"
1181 "(logand 7) @result{} 7\n"
1182 "(logand #b111 #b011 #b001) @result{} 1\n"
1184 #define FUNC_NAME s_scm_logand
1188 if (SCM_UNBNDP (n2
))
1190 if (SCM_UNBNDP (n1
))
1191 return SCM_MAKINUM (-1);
1192 else if (!SCM_NUMBERP (n1
))
1193 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1194 else if (SCM_NUMBERP (n1
))
1197 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1202 nn1
= SCM_INUM (n1
);
1205 long nn2
= SCM_INUM (n2
);
1206 return SCM_MAKINUM (nn1
& nn2
);
1208 else if SCM_BIGP (n2
)
1214 SCM result_z
= scm_i_mkbig ();
1216 mpz_init_set_si (nn1_z
, nn1
);
1217 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1218 scm_remember_upto_here_1 (n2
);
1220 return scm_i_normbig (result_z
);
1224 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1226 else if (SCM_BIGP (n1
))
1231 nn1
= SCM_INUM (n1
);
1234 else if (SCM_BIGP (n2
))
1236 SCM result_z
= scm_i_mkbig ();
1237 mpz_and (SCM_I_BIG_MPZ (result_z
),
1239 SCM_I_BIG_MPZ (n2
));
1240 scm_remember_upto_here_2 (n1
, n2
);
1241 return scm_i_normbig (result_z
);
1244 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1247 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1252 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1254 "Return the bitwise OR of the integer arguments.\n\n"
1256 "(logior) @result{} 0\n"
1257 "(logior 7) @result{} 7\n"
1258 "(logior #b000 #b001 #b011) @result{} 3\n"
1260 #define FUNC_NAME s_scm_logior
1264 if (SCM_UNBNDP (n2
))
1266 if (SCM_UNBNDP (n1
))
1268 else if (SCM_NUMBERP (n1
))
1271 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1276 nn1
= SCM_INUM (n1
);
1279 long nn2
= SCM_INUM (n2
);
1280 return SCM_MAKINUM (nn1
| nn2
);
1282 else if (SCM_BIGP (n2
))
1288 SCM result_z
= scm_i_mkbig ();
1290 mpz_init_set_si (nn1_z
, nn1
);
1291 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1292 scm_remember_upto_here_1 (n2
);
1298 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1300 else if (SCM_BIGP (n1
))
1305 nn1
= SCM_INUM (n1
);
1308 else if (SCM_BIGP (n2
))
1310 SCM result_z
= scm_i_mkbig ();
1311 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1313 SCM_I_BIG_MPZ (n2
));
1314 scm_remember_upto_here_2 (n1
, n2
);
1318 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1321 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1326 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1328 "Return the bitwise XOR of the integer arguments. A bit is\n"
1329 "set in the result if it is set in an odd number of arguments.\n"
1331 "(logxor) @result{} 0\n"
1332 "(logxor 7) @result{} 7\n"
1333 "(logxor #b000 #b001 #b011) @result{} 2\n"
1334 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1336 #define FUNC_NAME s_scm_logxor
1340 if (SCM_UNBNDP (n2
))
1342 if (SCM_UNBNDP (n1
))
1344 else if (SCM_NUMBERP (n1
))
1347 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1352 nn1
= SCM_INUM (n1
);
1355 long nn2
= SCM_INUM (n2
);
1356 return SCM_MAKINUM (nn1
^ nn2
);
1358 else if (SCM_BIGP (n2
))
1362 SCM result_z
= scm_i_mkbig ();
1364 mpz_init_set_si (nn1_z
, nn1
);
1365 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1366 scm_remember_upto_here_1 (n2
);
1368 return scm_i_normbig (result_z
);
1372 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1374 else if (SCM_BIGP (n1
))
1379 nn1
= SCM_INUM (n1
);
1382 else if (SCM_BIGP (n2
))
1384 SCM result_z
= scm_i_mkbig ();
1385 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1387 SCM_I_BIG_MPZ (n2
));
1388 scm_remember_upto_here_2 (n1
, n2
);
1389 return scm_i_normbig (result_z
);
1392 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1395 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1400 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1403 "(logtest j k) @equiv{} (not (zero? (logand j k)))\n\n"
1404 "(logtest #b0100 #b1011) @result{} #f\n"
1405 "(logtest #b0100 #b0111) @result{} #t\n"
1407 #define FUNC_NAME s_scm_logtest
1416 long nk
= SCM_INUM (k
);
1417 return SCM_BOOL (nj
& nk
);
1419 else if (SCM_BIGP (k
))
1427 mpz_init_set_si (nj_z
, nj
);
1428 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1429 scm_remember_upto_here_1 (k
);
1430 result
= SCM_BOOL (mpz_sgn (nj_z
) != 0);
1436 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1438 else if (SCM_BIGP (j
))
1446 else if (SCM_BIGP (k
))
1450 mpz_init (result_z
);
1454 scm_remember_upto_here_2 (j
, k
);
1455 result
= SCM_BOOL (mpz_sgn (result_z
) != 0);
1456 mpz_clear (result_z
);
1460 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1463 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1468 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1471 "(logbit? index j) @equiv{} (logtest (integer-expt 2 index) j)\n\n"
1472 "(logbit? 0 #b1101) @result{} #t\n"
1473 "(logbit? 1 #b1101) @result{} #f\n"
1474 "(logbit? 2 #b1101) @result{} #t\n"
1475 "(logbit? 3 #b1101) @result{} #t\n"
1476 "(logbit? 4 #b1101) @result{} #f\n"
1478 #define FUNC_NAME s_scm_logbit_p
1480 unsigned long int iindex
;
1482 SCM_VALIDATE_INUM_MIN (SCM_ARG1
, index
, 0);
1483 iindex
= (unsigned long int) SCM_INUM (index
);
1487 /* bits above what's in an inum follow the sign bit */
1488 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1489 return SCM_BOOL ((1L << iindex
) & SCM_INUM (j
));
1491 else if (SCM_BIGP (j
))
1493 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1494 scm_remember_upto_here_1 (j
);
1495 return SCM_BOOL (val
);
1498 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1503 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1505 "Return the integer which is the ones-complement of the integer\n"
1509 "(number->string (lognot #b10000000) 2)\n"
1510 " @result{} \"-10000001\"\n"
1511 "(number->string (lognot #b0) 2)\n"
1512 " @result{} \"-1\"\n"
1514 #define FUNC_NAME s_scm_lognot
1516 if (SCM_INUMP (n
)) {
1517 /* No overflow here, just need to toggle all the bits making up the inum.
1518 Enhancement: No need to strip the tag and add it back, could just xor
1519 a block of 1 bits, if that worked with the various debug versions of
1521 return SCM_MAKINUM (~ SCM_INUM (n
));
1523 } else if (SCM_BIGP (n
)) {
1524 SCM result
= scm_i_mkbig ();
1525 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1526 scm_remember_upto_here_1 (n
);
1530 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1535 /* returns 0 if IN is not an integer. OUT must already be
1538 coerce_to_big (SCM in
, mpz_t out
)
1541 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1542 else if (SCM_INUMP (in
))
1543 mpz_set_si (out
, SCM_INUM (in
));
1550 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1551 (SCM n
, SCM k
, SCM m
),
1552 "Return @var{n} raised to the integer exponent\n"
1553 "@var{k}, modulo @var{m}.\n"
1556 "(modulo-expt 2 3 5)\n"
1559 #define FUNC_NAME s_scm_modulo_expt
1565 /* There are two classes of error we might encounter --
1566 1) Math errors, which we'll report by calling scm_num_overflow,
1568 2) wrong-type errors, which of course we'll report by calling
1570 We don't report those errors immediately, however; instead we do
1571 some cleanup first. These variables tell us which error (if
1572 any) we should report after cleaning up.
1574 int report_overflow
= 0;
1576 int position_of_wrong_type
= 0;
1577 SCM value_of_wrong_type
= SCM_INUM0
;
1579 SCM result
= SCM_UNDEFINED
;
1585 if (SCM_EQ_P (m
, SCM_INUM0
))
1587 report_overflow
= 1;
1591 if (!coerce_to_big (n
, n_tmp
))
1593 value_of_wrong_type
= n
;
1594 position_of_wrong_type
= 1;
1598 if (!coerce_to_big (k
, k_tmp
))
1600 value_of_wrong_type
= k
;
1601 position_of_wrong_type
= 2;
1605 if (!coerce_to_big (m
, m_tmp
))
1607 value_of_wrong_type
= m
;
1608 position_of_wrong_type
= 3;
1612 /* if the exponent K is negative, and we simply call mpz_powm, we
1613 will get a divide-by-zero exception when an inverse 1/n mod m
1614 doesn't exist (or is not unique). Since exceptions are hard to
1615 handle, we'll attempt the inversion "by hand" -- that way, we get
1616 a simple failure code, which is easy to handle. */
1618 if (-1 == mpz_sgn (k_tmp
))
1620 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1622 report_overflow
= 1;
1625 mpz_neg (k_tmp
, k_tmp
);
1628 result
= scm_i_mkbig ();
1629 mpz_powm (SCM_I_BIG_MPZ (result
),
1634 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1635 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1642 if (report_overflow
)
1643 scm_num_overflow (FUNC_NAME
);
1645 if (position_of_wrong_type
)
1646 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1647 value_of_wrong_type
);
1649 return scm_i_normbig (result
);
1653 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1655 "Return @var{n} raised to the non-negative integer exponent\n"
1659 "(integer-expt 2 5)\n"
1661 "(integer-expt -3 3)\n"
1664 #define FUNC_NAME s_scm_integer_expt
1667 SCM z_i2
= SCM_BOOL_F
;
1669 SCM acc
= SCM_MAKINUM (1L);
1671 /* 0^0 == 1 according to R5RS */
1672 if (SCM_EQ_P (n
, SCM_INUM0
) || SCM_EQ_P (n
, acc
))
1673 return SCM_FALSEP (scm_zero_p(k
)) ? n
: acc
;
1674 else if (SCM_EQ_P (n
, SCM_MAKINUM (-1L)))
1675 return SCM_FALSEP (scm_even_p (k
)) ? n
: acc
;
1679 else if (SCM_BIGP (k
))
1681 z_i2
= scm_i_clonebig (k
, 1);
1682 scm_remember_upto_here_1 (k
);
1685 else if (SCM_REALP (k
))
1687 double r
= SCM_REAL_VALUE (k
);
1689 SCM_WRONG_TYPE_ARG (2, k
);
1690 if ((r
> SCM_MOST_POSITIVE_FIXNUM
) || (r
< SCM_MOST_NEGATIVE_FIXNUM
))
1692 z_i2
= scm_i_mkbig ();
1693 mpz_set_d (SCM_I_BIG_MPZ (z_i2
), r
);
1702 SCM_WRONG_TYPE_ARG (2, k
);
1706 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1708 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1709 n
= scm_divide (n
, SCM_UNDEFINED
);
1713 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1717 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1719 return scm_product (acc
, n
);
1721 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1722 acc
= scm_product (acc
, n
);
1723 n
= scm_product (n
, n
);
1724 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1732 n
= scm_divide (n
, SCM_UNDEFINED
);
1739 return scm_product (acc
, n
);
1741 acc
= scm_product (acc
, n
);
1742 n
= scm_product (n
, n
);
1749 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1751 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1752 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1754 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1755 "@var{cnt} is negative it's a division, rounded towards negative\n"
1756 "infinity. (Note that this is not the same rounding as\n"
1757 "@code{quotient} does.)\n"
1759 "With @var{n} viewed as an infinite precision twos complement,\n"
1760 "@code{ash} means a left shift introducing zero bits, or a right\n"
1761 "shift dropping bits.\n"
1764 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1765 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1767 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1768 "(ash -23 -2) @result{} -6\n"
1770 #define FUNC_NAME s_scm_ash
1774 SCM_VALIDATE_INUM (2, cnt
);
1776 bits_to_shift
= SCM_INUM (cnt
);
1778 if (bits_to_shift
< 0)
1780 /* Shift right by abs(cnt) bits. This is realized as a division
1781 by div:=2^abs(cnt). However, to guarantee the floor
1782 rounding, negative values require some special treatment.
1784 SCM div
= scm_integer_expt (SCM_MAKINUM (2),
1785 SCM_MAKINUM (-bits_to_shift
));
1787 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1788 if (SCM_FALSEP (scm_negative_p (n
)))
1789 return scm_quotient (n
, div
);
1791 return scm_sum (SCM_MAKINUM (-1L),
1792 scm_quotient (scm_sum (SCM_MAKINUM (1L), n
), div
));
1795 /* Shift left is done by multiplication with 2^CNT */
1796 return scm_product (n
, scm_integer_expt (SCM_MAKINUM (2), cnt
));
1801 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1802 (SCM n
, SCM start
, SCM end
),
1803 "Return the integer composed of the @var{start} (inclusive)\n"
1804 "through @var{end} (exclusive) bits of @var{n}. The\n"
1805 "@var{start}th bit becomes the 0-th bit in the result.\n"
1808 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1809 " @result{} \"1010\"\n"
1810 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1811 " @result{} \"10110\"\n"
1813 #define FUNC_NAME s_scm_bit_extract
1815 unsigned long int istart
, iend
, bits
;
1816 SCM_VALIDATE_INUM_MIN_COPY (2, start
,0, istart
);
1817 SCM_VALIDATE_INUM_MIN_COPY (3, end
, 0, iend
);
1818 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1820 /* how many bits to keep */
1821 bits
= iend
- istart
;
1825 long int in
= SCM_INUM (n
);
1827 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1828 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1829 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1831 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1833 /* Since we emulate two's complement encoded numbers, this
1834 * special case requires us to produce a result that has
1835 * more bits than can be stored in a fixnum.
1837 SCM result
= scm_i_long2big (in
);
1838 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1843 /* mask down to requisite bits */
1844 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1845 return SCM_MAKINUM (in
& ((1L << bits
) - 1));
1847 else if (SCM_BIGP (n
))
1852 result
= SCM_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1856 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1857 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1858 such bits into a ulong. */
1859 result
= scm_i_mkbig ();
1860 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1861 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1862 result
= scm_i_normbig (result
);
1864 scm_remember_upto_here_1 (n
);
1868 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1873 static const char scm_logtab
[] = {
1874 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1877 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1879 "Return the number of bits in integer @var{n}. If integer is\n"
1880 "positive, the 1-bits in its binary representation are counted.\n"
1881 "If negative, the 0-bits in its two's-complement binary\n"
1882 "representation are counted. If 0, 0 is returned.\n"
1885 "(logcount #b10101010)\n"
1892 #define FUNC_NAME s_scm_logcount
1896 unsigned long int c
= 0;
1897 long int nn
= SCM_INUM (n
);
1902 c
+= scm_logtab
[15 & nn
];
1905 return SCM_MAKINUM (c
);
1907 else if (SCM_BIGP (n
))
1909 unsigned long count
;
1910 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1911 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1913 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1914 scm_remember_upto_here_1 (n
);
1915 return SCM_MAKINUM (count
);
1918 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1923 static const char scm_ilentab
[] = {
1924 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1928 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1930 "Return the number of bits necessary to represent @var{n}.\n"
1933 "(integer-length #b10101010)\n"
1935 "(integer-length 0)\n"
1937 "(integer-length #b1111)\n"
1940 #define FUNC_NAME s_scm_integer_length
1944 unsigned long int c
= 0;
1946 long int nn
= SCM_INUM (n
);
1952 l
= scm_ilentab
[15 & nn
];
1955 return SCM_MAKINUM (c
- 4 + l
);
1957 else if (SCM_BIGP (n
))
1959 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1960 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1961 1 too big, so check for that and adjust. */
1962 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1963 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1964 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1965 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1967 scm_remember_upto_here_1 (n
);
1968 return SCM_MAKINUM (size
);
1971 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1975 /*** NUMBERS -> STRINGS ***/
1976 #define SCM_MAX_DBL_PREC 60
1977 #define SCM_MAX_DBL_RADIX 36
1979 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
1980 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
1981 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
1984 void init_dblprec(int *prec
, int radix
) {
1985 /* determine floating point precision by adding successively
1986 smaller increments to 1.0 until it is considered == 1.0 */
1987 double f
= ((double)1.0)/radix
;
1988 double fsum
= 1.0 + f
;
1993 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2005 void init_fx_radix(double *fx_list
, int radix
)
2007 /* initialize a per-radix list of tolerances. When added
2008 to a number < 1.0, we can determine if we should raund
2009 up and quit converting a number to a string. */
2013 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2014 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2017 /* use this array as a way to generate a single digit */
2018 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2021 idbl2str (double f
, char *a
, int radix
)
2023 int efmt
, dpt
, d
, i
, wp
;
2025 #ifdef DBL_MIN_10_EXP
2028 #endif /* DBL_MIN_10_EXP */
2033 radix
> SCM_MAX_DBL_RADIX
)
2035 /* revert to existing behavior */
2039 wp
= scm_dblprec
[radix
-2];
2040 fx
= fx_per_radix
[radix
-2];
2044 #ifdef HAVE_COPYSIGN
2045 double sgn
= copysign (1.0, f
);
2050 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2056 strcpy (a
, "-inf.0");
2058 strcpy (a
, "+inf.0");
2061 else if (xisnan (f
))
2063 strcpy (a
, "+nan.0");
2073 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2074 make-uniform-vector, from causing infinite loops. */
2075 /* just do the checking...if it passes, we do the conversion for our
2076 radix again below */
2083 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2091 while (f_cpy
> 10.0)
2094 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2115 if (f
+ fx
[wp
] >= radix
)
2122 /* adding 9999 makes this equivalent to abs(x) % 3 */
2123 dpt
= (exp
+ 9999) % 3;
2127 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2149 a
[ch
++] = number_chars
[d
];
2152 if (f
+ fx
[wp
] >= 1.0)
2154 a
[ch
- 1] = number_chars
[d
+1];
2166 if ((dpt
> 4) && (exp
> 6))
2168 d
= (a
[0] == '-' ? 2 : 1);
2169 for (i
= ch
++; i
> d
; i
--)
2182 if (a
[ch
- 1] == '.')
2183 a
[ch
++] = '0'; /* trailing zero */
2192 for (i
= radix
; i
<= exp
; i
*= radix
);
2193 for (i
/= radix
; i
; i
/= radix
)
2195 a
[ch
++] = number_chars
[exp
/ i
];
2203 iflo2str (SCM flt
, char *str
, int radix
)
2206 if (SCM_REALP (flt
))
2207 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2210 i
= idbl2str (SCM_COMPLEX_REAL (flt
), str
, radix
);
2211 if (SCM_COMPLEX_IMAG (flt
) != 0.0)
2213 double imag
= SCM_COMPLEX_IMAG (flt
);
2214 /* Don't output a '+' for negative numbers or for Inf and
2215 NaN. They will provide their own sign. */
2216 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2218 i
+= idbl2str (imag
, &str
[i
], radix
);
2225 /* convert a long to a string (unterminated). returns the number of
2226 characters in the result.
2228 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2230 scm_iint2str (long num
, int rad
, char *p
)
2234 unsigned long n
= (num
< 0) ? -num
: num
;
2236 for (n
/= rad
; n
> 0; n
/= rad
)
2253 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2258 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2260 "Return a string holding the external representation of the\n"
2261 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2262 "inexact, a radix of 10 will be used.")
2263 #define FUNC_NAME s_scm_number_to_string
2267 if (SCM_UNBNDP (radix
))
2271 SCM_VALIDATE_INUM (2, radix
);
2272 base
= SCM_INUM (radix
);
2273 /* FIXME: ask if range limit was OK, and if so, document */
2274 SCM_ASSERT_RANGE (2, radix
, (base
>= 2) && (base
<= 36));
2279 char num_buf
[SCM_INTBUFLEN
];
2280 size_t length
= scm_iint2str (SCM_INUM (n
), base
, num_buf
);
2281 return scm_mem2string (num_buf
, length
);
2283 else if (SCM_BIGP (n
))
2285 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2286 scm_remember_upto_here_1 (n
);
2287 return scm_take0str (str
);
2289 else if (SCM_FRACTIONP (n
))
2291 scm_i_fraction_reduce (n
);
2292 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2293 scm_mem2string ("/", 1),
2294 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2296 else if (SCM_INEXACTP (n
))
2298 char num_buf
[FLOBUFLEN
];
2299 return scm_mem2string (num_buf
, iflo2str (n
, num_buf
, base
));
2302 SCM_WRONG_TYPE_ARG (1, n
);
2307 /* These print routines used to be stubbed here so that scm_repl.c
2308 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2311 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2313 char num_buf
[FLOBUFLEN
];
2314 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2319 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2322 char num_buf
[FLOBUFLEN
];
2323 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2328 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2331 scm_i_fraction_reduce (sexp
);
2332 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2333 scm_lfwrite (SCM_STRING_CHARS (str
), SCM_STRING_LENGTH (str
), port
);
2334 scm_remember_upto_here_1 (str
);
2339 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2341 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2342 scm_remember_upto_here_1 (exp
);
2343 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2347 /*** END nums->strs ***/
2350 /*** STRINGS -> NUMBERS ***/
2352 /* The following functions implement the conversion from strings to numbers.
2353 * The implementation somehow follows the grammar for numbers as it is given
2354 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2355 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2356 * points should be noted about the implementation:
2357 * * Each function keeps a local index variable 'idx' that points at the
2358 * current position within the parsed string. The global index is only
2359 * updated if the function could parse the corresponding syntactic unit
2361 * * Similarly, the functions keep track of indicators of inexactness ('#',
2362 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2363 * global exactness information is only updated after each part has been
2364 * successfully parsed.
2365 * * Sequences of digits are parsed into temporary variables holding fixnums.
2366 * Only if these fixnums would overflow, the result variables are updated
2367 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2368 * the temporary variables holding the fixnums are cleared, and the process
2369 * starts over again. If for example fixnums were able to store five decimal
2370 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2371 * and the result was computed as 12345 * 100000 + 67890. In other words,
2372 * only every five digits two bignum operations were performed.
2375 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2377 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2379 /* In non ASCII-style encodings the following macro might not work. */
2380 #define XDIGIT2UINT(d) \
2381 (isdigit ((int) (unsigned char) d) \
2383 : tolower ((int) (unsigned char) d) - 'a' + 10)
2386 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2387 unsigned int radix
, enum t_exactness
*p_exactness
)
2389 unsigned int idx
= *p_idx
;
2390 unsigned int hash_seen
= 0;
2391 scm_t_bits shift
= 1;
2393 unsigned int digit_value
;
2401 if (!isxdigit ((int) (unsigned char) c
))
2403 digit_value
= XDIGIT2UINT (c
);
2404 if (digit_value
>= radix
)
2408 result
= SCM_MAKINUM (digit_value
);
2412 if (isxdigit ((int) (unsigned char) c
))
2416 digit_value
= XDIGIT2UINT (c
);
2417 if (digit_value
>= radix
)
2429 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2431 result
= scm_product (result
, SCM_MAKINUM (shift
));
2433 result
= scm_sum (result
, SCM_MAKINUM (add
));
2440 shift
= shift
* radix
;
2441 add
= add
* radix
+ digit_value
;
2446 result
= scm_product (result
, SCM_MAKINUM (shift
));
2448 result
= scm_sum (result
, SCM_MAKINUM (add
));
2452 *p_exactness
= INEXACT
;
2458 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2459 * covers the parts of the rules that start at a potential point. The value
2460 * of the digits up to the point have been parsed by the caller and are given
2461 * in variable result. The content of *p_exactness indicates, whether a hash
2462 * has already been seen in the digits before the point.
2465 /* In non ASCII-style encodings the following macro might not work. */
2466 #define DIGIT2UINT(d) ((d) - '0')
2469 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2470 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2472 unsigned int idx
= *p_idx
;
2473 enum t_exactness x
= *p_exactness
;
2478 if (mem
[idx
] == '.')
2480 scm_t_bits shift
= 1;
2482 unsigned int digit_value
;
2483 SCM big_shift
= SCM_MAKINUM (1);
2489 if (isdigit ((int) (unsigned char) c
))
2494 digit_value
= DIGIT2UINT (c
);
2505 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2507 big_shift
= scm_product (big_shift
, SCM_MAKINUM (shift
));
2508 result
= scm_product (result
, SCM_MAKINUM (shift
));
2510 result
= scm_sum (result
, SCM_MAKINUM (add
));
2518 add
= add
* 10 + digit_value
;
2524 big_shift
= scm_product (big_shift
, SCM_MAKINUM (shift
));
2525 result
= scm_product (result
, SCM_MAKINUM (shift
));
2526 result
= scm_sum (result
, SCM_MAKINUM (add
));
2529 result
= scm_divide (result
, big_shift
);
2531 /* We've seen a decimal point, thus the value is implicitly inexact. */
2543 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2570 if (!isdigit ((int) (unsigned char) c
))
2574 exponent
= DIGIT2UINT (c
);
2578 if (isdigit ((int) (unsigned char) c
))
2581 if (exponent
<= SCM_MAXEXP
)
2582 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2588 if (exponent
> SCM_MAXEXP
)
2590 size_t exp_len
= idx
- start
;
2591 SCM exp_string
= scm_mem2string (&mem
[start
], exp_len
);
2592 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2593 scm_out_of_range ("string->number", exp_num
);
2596 e
= scm_integer_expt (SCM_MAKINUM (10), SCM_MAKINUM (exponent
));
2598 result
= scm_product (result
, e
);
2600 result
= scm_divide2real (result
, e
);
2602 /* We've seen an exponent, thus the value is implicitly inexact. */
2620 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2623 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2624 unsigned int radix
, enum t_exactness
*p_exactness
)
2626 unsigned int idx
= *p_idx
;
2632 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2638 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2640 enum t_exactness x
= EXACT
;
2642 /* Cobble up the fractional part. We might want to set the
2643 NaN's mantissa from it. */
2645 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2650 if (mem
[idx
] == '.')
2654 else if (idx
+ 1 == len
)
2656 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2659 result
= mem2decimal_from_point (SCM_MAKINUM (0), mem
, len
,
2660 p_idx
, p_exactness
);
2664 enum t_exactness x
= EXACT
;
2667 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2668 if (SCM_FALSEP (uinteger
))
2673 else if (mem
[idx
] == '/')
2679 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2680 if (SCM_FALSEP (divisor
))
2683 /* both are int/big here, I assume */
2684 result
= scm_make_ratio (uinteger
, divisor
);
2686 else if (radix
== 10)
2688 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2689 if (SCM_FALSEP (result
))
2700 /* When returning an inexact zero, make sure it is represented as a
2701 floating point value so that we can change its sign.
2703 if (SCM_EQ_P (result
, SCM_MAKINUM(0)) && *p_exactness
== INEXACT
)
2704 result
= scm_make_real (0.0);
2710 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2713 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2714 unsigned int radix
, enum t_exactness
*p_exactness
)
2738 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2739 if (SCM_FALSEP (ureal
))
2741 /* input must be either +i or -i */
2746 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2752 return scm_make_rectangular (SCM_MAKINUM (0), SCM_MAKINUM (sign
));
2759 if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2760 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2769 /* either +<ureal>i or -<ureal>i */
2776 return scm_make_rectangular (SCM_MAKINUM (0), ureal
);
2779 /* polar input: <real>@<real>. */
2804 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2805 if (SCM_FALSEP (angle
))
2810 if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2811 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2813 result
= scm_make_polar (ureal
, angle
);
2818 /* expecting input matching <real>[+-]<ureal>?i */
2825 int sign
= (c
== '+') ? 1 : -1;
2826 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2828 if (SCM_FALSEP (imag
))
2829 imag
= SCM_MAKINUM (sign
);
2830 else if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2831 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2835 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2842 return scm_make_rectangular (ureal
, imag
);
2851 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2853 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2856 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2858 unsigned int idx
= 0;
2859 unsigned int radix
= NO_RADIX
;
2860 enum t_exactness forced_x
= NO_EXACTNESS
;
2861 enum t_exactness implicit_x
= EXACT
;
2864 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2865 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2867 switch (mem
[idx
+ 1])
2870 if (radix
!= NO_RADIX
)
2875 if (radix
!= NO_RADIX
)
2880 if (forced_x
!= NO_EXACTNESS
)
2885 if (forced_x
!= NO_EXACTNESS
)
2890 if (radix
!= NO_RADIX
)
2895 if (radix
!= NO_RADIX
)
2905 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2906 if (radix
== NO_RADIX
)
2907 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2909 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2911 if (SCM_FALSEP (result
))
2917 if (SCM_INEXACTP (result
))
2918 return scm_inexact_to_exact (result
);
2922 if (SCM_INEXACTP (result
))
2925 return scm_exact_to_inexact (result
);
2928 if (implicit_x
== INEXACT
)
2930 if (SCM_INEXACTP (result
))
2933 return scm_exact_to_inexact (result
);
2941 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2942 (SCM string
, SCM radix
),
2943 "Return a number of the maximally precise representation\n"
2944 "expressed by the given @var{string}. @var{radix} must be an\n"
2945 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2946 "is a default radix that may be overridden by an explicit radix\n"
2947 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2948 "supplied, then the default radix is 10. If string is not a\n"
2949 "syntactically valid notation for a number, then\n"
2950 "@code{string->number} returns @code{#f}.")
2951 #define FUNC_NAME s_scm_string_to_number
2955 SCM_VALIDATE_STRING (1, string
);
2956 SCM_VALIDATE_INUM_MIN_DEF_COPY (2, radix
,2,10, base
);
2957 answer
= scm_i_mem2number (SCM_STRING_CHARS (string
),
2958 SCM_STRING_LENGTH (string
),
2960 return scm_return_first (answer
, string
);
2965 /*** END strs->nums ***/
2969 scm_make_real (double x
)
2971 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
2973 SCM_REAL_VALUE (z
) = x
;
2979 scm_make_complex (double x
, double y
)
2982 return scm_make_real (x
);
2986 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
2988 SCM_COMPLEX_REAL (z
) = x
;
2989 SCM_COMPLEX_IMAG (z
) = y
;
2996 scm_bigequal (SCM x
, SCM y
)
2998 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2999 scm_remember_upto_here_2 (x
, y
);
3000 return SCM_BOOL (0 == result
);
3004 scm_real_equalp (SCM x
, SCM y
)
3006 return SCM_BOOL (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3010 scm_complex_equalp (SCM x
, SCM y
)
3012 return SCM_BOOL (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3013 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3017 scm_i_fraction_equalp (SCM x
, SCM y
)
3019 scm_i_fraction_reduce (x
);
3020 scm_i_fraction_reduce (y
);
3021 if (SCM_FALSEP (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3022 SCM_FRACTION_NUMERATOR (y
)))
3023 || SCM_FALSEP (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3024 SCM_FRACTION_DENOMINATOR (y
))))
3031 SCM_REGISTER_PROC (s_number_p
, "number?", 1, 0, 0, scm_number_p
);
3032 /* "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3033 * "else. Note that the sets of complex, real, rational and\n"
3034 * "integer values form subsets of the set of numbers, i. e. the\n"
3035 * "predicate will be fulfilled for any number."
3037 SCM_DEFINE (scm_number_p
, "complex?", 1, 0, 0,
3039 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3040 "otherwise. Note that the sets of real, rational and integer\n"
3041 "values form subsets of the set of complex numbers, i. e. the\n"
3042 "predicate will also be fulfilled if @var{x} is a real,\n"
3043 "rational or integer number.")
3044 #define FUNC_NAME s_scm_number_p
3046 return SCM_BOOL (SCM_NUMBERP (x
));
3051 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3053 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3054 "otherwise. Note that the set of integer values forms a subset of\n"
3055 "the set of real numbers, i. e. the predicate will also be\n"
3056 "fulfilled if @var{x} is an integer number.")
3057 #define FUNC_NAME s_scm_real_p
3059 /* we can't represent irrational numbers. */
3060 return scm_rational_p (x
);
3064 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3066 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3067 "otherwise. Note that the set of integer values forms a subset of\n"
3068 "the set of rational numbers, i. e. the predicate will also be\n"
3069 "fulfilled if @var{x} is an integer number.")
3070 #define FUNC_NAME s_scm_rational_p
3074 else if (SCM_IMP (x
))
3076 else if (SCM_BIGP (x
))
3078 else if (SCM_FRACTIONP (x
))
3080 else if (SCM_REALP (x
))
3081 /* due to their limited precision, all floating point numbers are
3082 rational as well. */
3090 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3092 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3094 #define FUNC_NAME s_scm_integer_p
3103 if (!SCM_INEXACTP (x
))
3105 if (SCM_COMPLEXP (x
))
3107 r
= SCM_REAL_VALUE (x
);
3115 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3117 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3119 #define FUNC_NAME s_scm_inexact_p
3121 if (SCM_INEXACTP (x
))
3123 if (SCM_NUMBERP (x
))
3125 SCM_WRONG_TYPE_ARG (1, x
);
3130 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3131 /* "Return @code{#t} if all parameters are numerically equal." */
3133 scm_num_eq_p (SCM x
, SCM y
)
3138 long xx
= SCM_INUM (x
);
3141 long yy
= SCM_INUM (y
);
3142 return SCM_BOOL (xx
== yy
);
3144 else if (SCM_BIGP (y
))
3146 else if (SCM_REALP (y
))
3147 return SCM_BOOL ((double) xx
== SCM_REAL_VALUE (y
));
3148 else if (SCM_COMPLEXP (y
))
3149 return SCM_BOOL (((double) xx
== SCM_COMPLEX_REAL (y
))
3150 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3151 else if (SCM_FRACTIONP (y
))
3154 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3156 else if (SCM_BIGP (x
))
3160 else if (SCM_BIGP (y
))
3162 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3163 scm_remember_upto_here_2 (x
, y
);
3164 return SCM_BOOL (0 == cmp
);
3166 else if (SCM_REALP (y
))
3169 if (xisnan (SCM_REAL_VALUE (y
)))
3171 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3172 scm_remember_upto_here_1 (x
);
3173 return SCM_BOOL (0 == cmp
);
3175 else if (SCM_COMPLEXP (y
))
3178 if (0.0 != SCM_COMPLEX_IMAG (y
))
3180 if (xisnan (SCM_COMPLEX_REAL (y
)))
3182 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3183 scm_remember_upto_here_1 (x
);
3184 return SCM_BOOL (0 == cmp
);
3186 else if (SCM_FRACTIONP (y
))
3189 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3191 else if (SCM_REALP (x
))
3194 return SCM_BOOL (SCM_REAL_VALUE (x
) == (double) SCM_INUM (y
));
3195 else if (SCM_BIGP (y
))
3198 if (xisnan (SCM_REAL_VALUE (x
)))
3200 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3201 scm_remember_upto_here_1 (y
);
3202 return SCM_BOOL (0 == cmp
);
3204 else if (SCM_REALP (y
))
3205 return SCM_BOOL (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3206 else if (SCM_COMPLEXP (y
))
3207 return SCM_BOOL ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3208 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3209 else if (SCM_FRACTIONP (y
))
3211 double xx
= SCM_REAL_VALUE (x
);
3215 return SCM_BOOL (xx
< 0.0);
3216 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3220 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3222 else if (SCM_COMPLEXP (x
))
3225 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == (double) SCM_INUM (y
))
3226 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3227 else if (SCM_BIGP (y
))
3230 if (0.0 != SCM_COMPLEX_IMAG (x
))
3232 if (xisnan (SCM_COMPLEX_REAL (x
)))
3234 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3235 scm_remember_upto_here_1 (y
);
3236 return SCM_BOOL (0 == cmp
);
3238 else if (SCM_REALP (y
))
3239 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3240 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3241 else if (SCM_COMPLEXP (y
))
3242 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3243 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3244 else if (SCM_FRACTIONP (y
))
3247 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3249 xx
= SCM_COMPLEX_REAL (x
);
3253 return SCM_BOOL (xx
< 0.0);
3254 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3258 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3260 else if (SCM_FRACTIONP (x
))
3264 else if (SCM_BIGP (y
))
3266 else if (SCM_REALP (y
))
3268 double yy
= SCM_REAL_VALUE (y
);
3272 return SCM_BOOL (0.0 < yy
);
3273 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3276 else if (SCM_COMPLEXP (y
))
3279 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3281 yy
= SCM_COMPLEX_REAL (y
);
3285 return SCM_BOOL (0.0 < yy
);
3286 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3289 else if (SCM_FRACTIONP (y
))
3290 return scm_i_fraction_equalp (x
, y
);
3292 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3295 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3299 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3300 done are good for inums, but for bignums an answer can almost always be
3301 had by just examining a few high bits of the operands, as done by GMP in
3302 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3303 of the float exponent to take into account. */
3305 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3306 /* "Return @code{#t} if the list of parameters is monotonically\n"
3310 scm_less_p (SCM x
, SCM y
)
3315 long xx
= SCM_INUM (x
);
3318 long yy
= SCM_INUM (y
);
3319 return SCM_BOOL (xx
< yy
);
3321 else if (SCM_BIGP (y
))
3323 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3324 scm_remember_upto_here_1 (y
);
3325 return SCM_BOOL (sgn
> 0);
3327 else if (SCM_REALP (y
))
3328 return SCM_BOOL ((double) xx
< SCM_REAL_VALUE (y
));
3329 else if (SCM_FRACTIONP (y
))
3331 /* "x < a/b" becomes "x*b < a" */
3333 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3334 y
= SCM_FRACTION_NUMERATOR (y
);
3338 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3340 else if (SCM_BIGP (x
))
3344 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3345 scm_remember_upto_here_1 (x
);
3346 return SCM_BOOL (sgn
< 0);
3348 else if (SCM_BIGP (y
))
3350 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3351 scm_remember_upto_here_2 (x
, y
);
3352 return SCM_BOOL (cmp
< 0);
3354 else if (SCM_REALP (y
))
3357 if (xisnan (SCM_REAL_VALUE (y
)))
3359 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3360 scm_remember_upto_here_1 (x
);
3361 return SCM_BOOL (cmp
< 0);
3363 else if (SCM_FRACTIONP (y
))
3366 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3368 else if (SCM_REALP (x
))
3371 return SCM_BOOL (SCM_REAL_VALUE (x
) < (double) SCM_INUM (y
));
3372 else if (SCM_BIGP (y
))
3375 if (xisnan (SCM_REAL_VALUE (x
)))
3377 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3378 scm_remember_upto_here_1 (y
);
3379 return SCM_BOOL (cmp
> 0);
3381 else if (SCM_REALP (y
))
3382 return SCM_BOOL (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3383 else if (SCM_FRACTIONP (y
))
3385 double xx
= SCM_REAL_VALUE (x
);
3389 return SCM_BOOL (xx
< 0.0);
3390 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3394 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3396 else if (SCM_FRACTIONP (x
))
3398 if (SCM_INUMP (y
) || SCM_BIGP (y
))
3400 /* "a/b < y" becomes "a < y*b" */
3401 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3402 x
= SCM_FRACTION_NUMERATOR (x
);
3405 else if (SCM_REALP (y
))
3407 double yy
= SCM_REAL_VALUE (y
);
3411 return SCM_BOOL (0.0 < yy
);
3412 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3415 else if (SCM_FRACTIONP (y
))
3417 /* "a/b < c/d" becomes "a*d < c*b" */
3418 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3419 SCM_FRACTION_DENOMINATOR (y
));
3420 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3421 SCM_FRACTION_DENOMINATOR (x
));
3427 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3430 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3434 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3435 /* "Return @code{#t} if the list of parameters is monotonically\n"
3438 #define FUNC_NAME s_scm_gr_p
3440 scm_gr_p (SCM x
, SCM y
)
3442 if (!SCM_NUMBERP (x
))
3443 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3444 else if (!SCM_NUMBERP (y
))
3445 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3447 return scm_less_p (y
, x
);
3452 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3453 /* "Return @code{#t} if the list of parameters is monotonically\n"
3456 #define FUNC_NAME s_scm_leq_p
3458 scm_leq_p (SCM x
, SCM y
)
3460 if (!SCM_NUMBERP (x
))
3461 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3462 else if (!SCM_NUMBERP (y
))
3463 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3464 else if (SCM_NFALSEP (scm_nan_p (x
)) || SCM_NFALSEP (scm_nan_p (y
)))
3467 return SCM_BOOL_NOT (scm_less_p (y
, x
));
3472 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3473 /* "Return @code{#t} if the list of parameters is monotonically\n"
3476 #define FUNC_NAME s_scm_geq_p
3478 scm_geq_p (SCM x
, SCM y
)
3480 if (!SCM_NUMBERP (x
))
3481 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3482 else if (!SCM_NUMBERP (y
))
3483 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3484 else if (SCM_NFALSEP (scm_nan_p (x
)) || SCM_NFALSEP (scm_nan_p (y
)))
3487 return SCM_BOOL_NOT (scm_less_p (x
, y
));
3492 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3493 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3500 return SCM_BOOL (SCM_EQ_P (z
, SCM_INUM0
));
3501 else if (SCM_BIGP (z
))
3503 else if (SCM_REALP (z
))
3504 return SCM_BOOL (SCM_REAL_VALUE (z
) == 0.0);
3505 else if (SCM_COMPLEXP (z
))
3506 return SCM_BOOL (SCM_COMPLEX_REAL (z
) == 0.0
3507 && SCM_COMPLEX_IMAG (z
) == 0.0);
3508 else if (SCM_FRACTIONP (z
))
3511 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3515 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3516 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3520 scm_positive_p (SCM x
)
3523 return SCM_BOOL (SCM_INUM (x
) > 0);
3524 else if (SCM_BIGP (x
))
3526 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3527 scm_remember_upto_here_1 (x
);
3528 return SCM_BOOL (sgn
> 0);
3530 else if (SCM_REALP (x
))
3531 return SCM_BOOL(SCM_REAL_VALUE (x
) > 0.0);
3532 else if (SCM_FRACTIONP (x
))
3533 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3535 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3539 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3540 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3544 scm_negative_p (SCM x
)
3547 return SCM_BOOL (SCM_INUM (x
) < 0);
3548 else if (SCM_BIGP (x
))
3550 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3551 scm_remember_upto_here_1 (x
);
3552 return SCM_BOOL (sgn
< 0);
3554 else if (SCM_REALP (x
))
3555 return SCM_BOOL(SCM_REAL_VALUE (x
) < 0.0);
3556 else if (SCM_FRACTIONP (x
))
3557 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3559 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3563 /* scm_min and scm_max return an inexact when either argument is inexact, as
3564 required by r5rs. On that basis, for exact/inexact combinations the
3565 exact is converted to inexact to compare and possibly return. This is
3566 unlike scm_less_p above which takes some trouble to preserve all bits in
3567 its test, such trouble is not required for min and max. */
3569 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3570 /* "Return the maximum of all parameter values."
3573 scm_max (SCM x
, SCM y
)
3578 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3579 else if (SCM_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3582 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3587 long xx
= SCM_INUM (x
);
3590 long yy
= SCM_INUM (y
);
3591 return (xx
< yy
) ? y
: x
;
3593 else if (SCM_BIGP (y
))
3595 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3596 scm_remember_upto_here_1 (y
);
3597 return (sgn
< 0) ? x
: y
;
3599 else if (SCM_REALP (y
))
3602 /* if y==NaN then ">" is false and we return NaN */
3603 return (z
> SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3605 else if (SCM_FRACTIONP (y
))
3608 return (SCM_FALSEP (scm_less_p (x
, y
)) ? x
: y
);
3611 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3613 else if (SCM_BIGP (x
))
3617 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3618 scm_remember_upto_here_1 (x
);
3619 return (sgn
< 0) ? y
: x
;
3621 else if (SCM_BIGP (y
))
3623 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3624 scm_remember_upto_here_2 (x
, y
);
3625 return (cmp
> 0) ? x
: y
;
3627 else if (SCM_REALP (y
))
3629 /* if y==NaN then xx>yy is false, so we return the NaN y */
3632 xx
= scm_i_big2dbl (x
);
3633 yy
= SCM_REAL_VALUE (y
);
3634 return (xx
> yy
? scm_make_real (xx
) : y
);
3636 else if (SCM_FRACTIONP (y
))
3641 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3643 else if (SCM_REALP (x
))
3647 double z
= SCM_INUM (y
);
3648 /* if x==NaN then "<" is false and we return NaN */
3649 return (SCM_REAL_VALUE (x
) < z
) ? scm_make_real (z
) : x
;
3651 else if (SCM_BIGP (y
))
3656 else if (SCM_REALP (y
))
3658 /* if x==NaN then our explicit check means we return NaN
3659 if y==NaN then ">" is false and we return NaN
3660 calling isnan is unavoidable, since it's the only way to know
3661 which of x or y causes any compares to be false */
3662 double xx
= SCM_REAL_VALUE (x
);
3663 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3665 else if (SCM_FRACTIONP (y
))
3667 double yy
= scm_i_fraction2double (y
);
3668 double xx
= SCM_REAL_VALUE (x
);
3669 return (xx
< yy
) ? scm_make_real (yy
) : x
;
3672 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3674 else if (SCM_FRACTIONP (x
))
3680 else if (SCM_BIGP (y
))
3684 else if (SCM_REALP (y
))
3686 double xx
= scm_i_fraction2double (x
);
3687 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_make_real (xx
);
3689 else if (SCM_FRACTIONP (y
))
3694 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3697 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3701 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3702 /* "Return the minium of all parameter values."
3705 scm_min (SCM x
, SCM y
)
3710 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3711 else if (SCM_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3714 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3719 long xx
= SCM_INUM (x
);
3722 long yy
= SCM_INUM (y
);
3723 return (xx
< yy
) ? x
: y
;
3725 else if (SCM_BIGP (y
))
3727 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3728 scm_remember_upto_here_1 (y
);
3729 return (sgn
< 0) ? y
: x
;
3731 else if (SCM_REALP (y
))
3734 /* if y==NaN then "<" is false and we return NaN */
3735 return (z
< SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3737 else if (SCM_FRACTIONP (y
))
3740 return (SCM_FALSEP (scm_less_p (x
, y
)) ? y
: x
);
3743 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3745 else if (SCM_BIGP (x
))
3749 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3750 scm_remember_upto_here_1 (x
);
3751 return (sgn
< 0) ? x
: y
;
3753 else if (SCM_BIGP (y
))
3755 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3756 scm_remember_upto_here_2 (x
, y
);
3757 return (cmp
> 0) ? y
: x
;
3759 else if (SCM_REALP (y
))
3761 /* if y==NaN then xx<yy is false, so we return the NaN y */
3764 xx
= scm_i_big2dbl (x
);
3765 yy
= SCM_REAL_VALUE (y
);
3766 return (xx
< yy
? scm_make_real (xx
) : y
);
3768 else if (SCM_FRACTIONP (y
))
3773 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3775 else if (SCM_REALP (x
))
3779 double z
= SCM_INUM (y
);
3780 /* if x==NaN then "<" is false and we return NaN */
3781 return (z
< SCM_REAL_VALUE (x
)) ? scm_make_real (z
) : x
;
3783 else if (SCM_BIGP (y
))
3788 else if (SCM_REALP (y
))
3790 /* if x==NaN then our explicit check means we return NaN
3791 if y==NaN then "<" is false and we return NaN
3792 calling isnan is unavoidable, since it's the only way to know
3793 which of x or y causes any compares to be false */
3794 double xx
= SCM_REAL_VALUE (x
);
3795 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3797 else if (SCM_FRACTIONP (y
))
3799 double yy
= scm_i_fraction2double (y
);
3800 double xx
= SCM_REAL_VALUE (x
);
3801 return (yy
< xx
) ? scm_make_real (yy
) : x
;
3804 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3806 else if (SCM_FRACTIONP (x
))
3812 else if (SCM_BIGP (y
))
3816 else if (SCM_REALP (y
))
3818 double xx
= scm_i_fraction2double (x
);
3819 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_make_real (xx
);
3821 else if (SCM_FRACTIONP (y
))
3826 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3829 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3833 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3834 /* "Return the sum of all parameter values. Return 0 if called without\n"
3838 scm_sum (SCM x
, SCM y
)
3842 if (SCM_NUMBERP (x
)) return x
;
3843 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3844 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3851 long xx
= SCM_INUM (x
);
3852 long yy
= SCM_INUM (y
);
3853 long int z
= xx
+ yy
;
3854 return SCM_FIXABLE (z
) ? SCM_MAKINUM (z
) : scm_i_long2big (z
);
3856 else if (SCM_BIGP (y
))
3861 else if (SCM_REALP (y
))
3863 long int xx
= SCM_INUM (x
);
3864 return scm_make_real (xx
+ SCM_REAL_VALUE (y
));
3866 else if (SCM_COMPLEXP (y
))
3868 long int xx
= SCM_INUM (x
);
3869 return scm_make_complex (xx
+ SCM_COMPLEX_REAL (y
),
3870 SCM_COMPLEX_IMAG (y
));
3872 else if (SCM_FRACTIONP (y
))
3873 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3874 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3875 SCM_FRACTION_DENOMINATOR (y
));
3877 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3878 } else if (SCM_BIGP (x
))
3885 inum
= SCM_INUM (y
);
3888 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3891 SCM result
= scm_i_mkbig ();
3892 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3893 scm_remember_upto_here_1 (x
);
3894 /* we know the result will have to be a bignum */
3897 return scm_i_normbig (result
);
3901 SCM result
= scm_i_mkbig ();
3902 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3903 scm_remember_upto_here_1 (x
);
3904 /* we know the result will have to be a bignum */
3907 return scm_i_normbig (result
);
3910 else if (SCM_BIGP (y
))
3912 SCM result
= scm_i_mkbig ();
3913 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3914 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3915 mpz_add (SCM_I_BIG_MPZ (result
),
3918 scm_remember_upto_here_2 (x
, y
);
3919 /* we know the result will have to be a bignum */
3922 return scm_i_normbig (result
);
3924 else if (SCM_REALP (y
))
3926 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3927 scm_remember_upto_here_1 (x
);
3928 return scm_make_real (result
);
3930 else if (SCM_COMPLEXP (y
))
3932 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3933 + SCM_COMPLEX_REAL (y
));
3934 scm_remember_upto_here_1 (x
);
3935 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
3937 else if (SCM_FRACTIONP (y
))
3938 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3939 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3940 SCM_FRACTION_DENOMINATOR (y
));
3942 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3944 else if (SCM_REALP (x
))
3947 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_INUM (y
));
3948 else if (SCM_BIGP (y
))
3950 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3951 scm_remember_upto_here_1 (y
);
3952 return scm_make_real (result
);
3954 else if (SCM_REALP (y
))
3955 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3956 else if (SCM_COMPLEXP (y
))
3957 return scm_make_complex (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3958 SCM_COMPLEX_IMAG (y
));
3959 else if (SCM_FRACTIONP (y
))
3960 return scm_make_real (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3962 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3964 else if (SCM_COMPLEXP (x
))
3967 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_INUM (y
),
3968 SCM_COMPLEX_IMAG (x
));
3969 else if (SCM_BIGP (y
))
3971 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3972 + SCM_COMPLEX_REAL (x
));
3973 scm_remember_upto_here_1 (y
);
3974 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (x
));
3976 else if (SCM_REALP (y
))
3977 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3978 SCM_COMPLEX_IMAG (x
));
3979 else if (SCM_COMPLEXP (y
))
3980 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3981 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3982 else if (SCM_FRACTIONP (y
))
3983 return scm_make_complex (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3984 SCM_COMPLEX_IMAG (x
));
3986 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3988 else if (SCM_FRACTIONP (x
))
3991 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3992 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3993 SCM_FRACTION_DENOMINATOR (x
));
3994 else if (SCM_BIGP (y
))
3995 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3996 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3997 SCM_FRACTION_DENOMINATOR (x
));
3998 else if (SCM_REALP (y
))
3999 return scm_make_real (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4000 else if (SCM_COMPLEXP (y
))
4001 return scm_make_complex (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4002 SCM_COMPLEX_IMAG (y
));
4003 else if (SCM_FRACTIONP (y
))
4004 /* a/b + c/d = (ad + bc) / bd */
4005 return scm_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4006 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4007 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4009 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4012 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4016 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4017 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4018 * the sum of all but the first argument are subtracted from the first
4020 #define FUNC_NAME s_difference
4022 scm_difference (SCM x
, SCM y
)
4027 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4031 long xx
= -SCM_INUM (x
);
4032 if (SCM_FIXABLE (xx
))
4033 return SCM_MAKINUM (xx
);
4035 return scm_i_long2big (xx
);
4037 else if (SCM_BIGP (x
))
4038 /* FIXME: do we really need to normalize here? */
4039 return scm_i_normbig (scm_i_clonebig (x
, 0));
4040 else if (SCM_REALP (x
))
4041 return scm_make_real (-SCM_REAL_VALUE (x
));
4042 else if (SCM_COMPLEXP (x
))
4043 return scm_make_complex (-SCM_COMPLEX_REAL (x
),
4044 -SCM_COMPLEX_IMAG (x
));
4045 else if (SCM_FRACTIONP (x
))
4046 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4047 SCM_FRACTION_DENOMINATOR (x
));
4049 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4056 long int xx
= SCM_INUM (x
);
4057 long int yy
= SCM_INUM (y
);
4058 long int z
= xx
- yy
;
4059 if (SCM_FIXABLE (z
))
4060 return SCM_MAKINUM (z
);
4062 return scm_i_long2big (z
);
4064 else if (SCM_BIGP (y
))
4066 /* inum-x - big-y */
4067 long xx
= SCM_INUM (x
);
4070 return scm_i_clonebig (y
, 0);
4073 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4074 SCM result
= scm_i_mkbig ();
4077 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4080 /* x - y == -(y + -x) */
4081 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4082 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4084 scm_remember_upto_here_1 (y
);
4086 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4087 /* we know the result will have to be a bignum */
4090 return scm_i_normbig (result
);
4093 else if (SCM_REALP (y
))
4095 long int xx
= SCM_INUM (x
);
4096 return scm_make_real (xx
- SCM_REAL_VALUE (y
));
4098 else if (SCM_COMPLEXP (y
))
4100 long int xx
= SCM_INUM (x
);
4101 return scm_make_complex (xx
- SCM_COMPLEX_REAL (y
),
4102 - SCM_COMPLEX_IMAG (y
));
4104 else if (SCM_FRACTIONP (y
))
4105 /* a - b/c = (ac - b) / c */
4106 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4107 SCM_FRACTION_NUMERATOR (y
)),
4108 SCM_FRACTION_DENOMINATOR (y
));
4110 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4112 else if (SCM_BIGP (x
))
4116 /* big-x - inum-y */
4117 long yy
= SCM_INUM (y
);
4118 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4120 scm_remember_upto_here_1 (x
);
4122 return SCM_FIXABLE (-yy
) ? SCM_MAKINUM (-yy
) : scm_long2num (-yy
);
4125 SCM result
= scm_i_mkbig ();
4128 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4130 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4131 scm_remember_upto_here_1 (x
);
4133 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4134 /* we know the result will have to be a bignum */
4137 return scm_i_normbig (result
);
4140 else if (SCM_BIGP (y
))
4142 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4143 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4144 SCM result
= scm_i_mkbig ();
4145 mpz_sub (SCM_I_BIG_MPZ (result
),
4148 scm_remember_upto_here_2 (x
, y
);
4149 /* we know the result will have to be a bignum */
4150 if ((sgn_x
== 1) && (sgn_y
== -1))
4152 if ((sgn_x
== -1) && (sgn_y
== 1))
4154 return scm_i_normbig (result
);
4156 else if (SCM_REALP (y
))
4158 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4159 scm_remember_upto_here_1 (x
);
4160 return scm_make_real (result
);
4162 else if (SCM_COMPLEXP (y
))
4164 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4165 - SCM_COMPLEX_REAL (y
));
4166 scm_remember_upto_here_1 (x
);
4167 return scm_make_complex (real_part
, - SCM_COMPLEX_IMAG (y
));
4169 else if (SCM_FRACTIONP (y
))
4170 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4171 SCM_FRACTION_NUMERATOR (y
)),
4172 SCM_FRACTION_DENOMINATOR (y
));
4173 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4175 else if (SCM_REALP (x
))
4178 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_INUM (y
));
4179 else if (SCM_BIGP (y
))
4181 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4182 scm_remember_upto_here_1 (x
);
4183 return scm_make_real (result
);
4185 else if (SCM_REALP (y
))
4186 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4187 else if (SCM_COMPLEXP (y
))
4188 return scm_make_complex (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4189 -SCM_COMPLEX_IMAG (y
));
4190 else if (SCM_FRACTIONP (y
))
4191 return scm_make_real (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4193 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4195 else if (SCM_COMPLEXP (x
))
4198 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_INUM (y
),
4199 SCM_COMPLEX_IMAG (x
));
4200 else if (SCM_BIGP (y
))
4202 double real_part
= (SCM_COMPLEX_REAL (x
)
4203 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4204 scm_remember_upto_here_1 (x
);
4205 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
4207 else if (SCM_REALP (y
))
4208 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4209 SCM_COMPLEX_IMAG (x
));
4210 else if (SCM_COMPLEXP (y
))
4211 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4212 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4213 else if (SCM_FRACTIONP (y
))
4214 return scm_make_complex (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4215 SCM_COMPLEX_IMAG (x
));
4217 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4219 else if (SCM_FRACTIONP (x
))
4222 /* a/b - c = (a - cb) / b */
4223 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4224 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4225 SCM_FRACTION_DENOMINATOR (x
));
4226 else if (SCM_BIGP (y
))
4227 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4228 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4229 SCM_FRACTION_DENOMINATOR (x
));
4230 else if (SCM_REALP (y
))
4231 return scm_make_real (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4232 else if (SCM_COMPLEXP (y
))
4233 return scm_make_complex (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4234 -SCM_COMPLEX_IMAG (y
));
4235 else if (SCM_FRACTIONP (y
))
4236 /* a/b - c/d = (ad - bc) / bd */
4237 return scm_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4238 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4239 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4241 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4244 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4249 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4250 /* "Return the product of all arguments. If called without arguments,\n"
4254 scm_product (SCM x
, SCM y
)
4259 return SCM_MAKINUM (1L);
4260 else if (SCM_NUMBERP (x
))
4263 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4275 case 0: return x
; break;
4276 case 1: return y
; break;
4281 long yy
= SCM_INUM (y
);
4283 SCM k
= SCM_MAKINUM (kk
);
4284 if ((kk
== SCM_INUM (k
)) && (kk
/ xx
== yy
))
4288 SCM result
= scm_i_long2big (xx
);
4289 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4290 return scm_i_normbig (result
);
4293 else if (SCM_BIGP (y
))
4295 SCM result
= scm_i_mkbig ();
4296 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4297 scm_remember_upto_here_1 (y
);
4300 else if (SCM_REALP (y
))
4301 return scm_make_real (xx
* SCM_REAL_VALUE (y
));
4302 else if (SCM_COMPLEXP (y
))
4303 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4304 xx
* SCM_COMPLEX_IMAG (y
));
4305 else if (SCM_FRACTIONP (y
))
4306 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4307 SCM_FRACTION_DENOMINATOR (y
));
4309 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4311 else if (SCM_BIGP (x
))
4318 else if (SCM_BIGP (y
))
4320 SCM result
= scm_i_mkbig ();
4321 mpz_mul (SCM_I_BIG_MPZ (result
),
4324 scm_remember_upto_here_2 (x
, y
);
4327 else if (SCM_REALP (y
))
4329 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4330 scm_remember_upto_here_1 (x
);
4331 return scm_make_real (result
);
4333 else if (SCM_COMPLEXP (y
))
4335 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4336 scm_remember_upto_here_1 (x
);
4337 return scm_make_complex (z
* SCM_COMPLEX_REAL (y
),
4338 z
* SCM_COMPLEX_IMAG (y
));
4340 else if (SCM_FRACTIONP (y
))
4341 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4342 SCM_FRACTION_DENOMINATOR (y
));
4344 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4346 else if (SCM_REALP (x
))
4349 return scm_make_real (SCM_INUM (y
) * SCM_REAL_VALUE (x
));
4350 else if (SCM_BIGP (y
))
4352 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4353 scm_remember_upto_here_1 (y
);
4354 return scm_make_real (result
);
4356 else if (SCM_REALP (y
))
4357 return scm_make_real (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4358 else if (SCM_COMPLEXP (y
))
4359 return scm_make_complex (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4360 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4361 else if (SCM_FRACTIONP (y
))
4362 return scm_make_real (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4364 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4366 else if (SCM_COMPLEXP (x
))
4369 return scm_make_complex (SCM_INUM (y
) * SCM_COMPLEX_REAL (x
),
4370 SCM_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4371 else if (SCM_BIGP (y
))
4373 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4374 scm_remember_upto_here_1 (y
);
4375 return scm_make_complex (z
* SCM_COMPLEX_REAL (x
),
4376 z
* SCM_COMPLEX_IMAG (x
));
4378 else if (SCM_REALP (y
))
4379 return scm_make_complex (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4380 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4381 else if (SCM_COMPLEXP (y
))
4383 return scm_make_complex (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4384 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4385 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4386 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4388 else if (SCM_FRACTIONP (y
))
4390 double yy
= scm_i_fraction2double (y
);
4391 return scm_make_complex (yy
* SCM_COMPLEX_REAL (x
),
4392 yy
* SCM_COMPLEX_IMAG (x
));
4395 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4397 else if (SCM_FRACTIONP (x
))
4400 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4401 SCM_FRACTION_DENOMINATOR (x
));
4402 else if (SCM_BIGP (y
))
4403 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4404 SCM_FRACTION_DENOMINATOR (x
));
4405 else if (SCM_REALP (y
))
4406 return scm_make_real (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4407 else if (SCM_COMPLEXP (y
))
4409 double xx
= scm_i_fraction2double (x
);
4410 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4411 xx
* SCM_COMPLEX_IMAG (y
));
4413 else if (SCM_FRACTIONP (y
))
4414 /* a/b * c/d = ac / bd */
4415 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4416 SCM_FRACTION_NUMERATOR (y
)),
4417 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4418 SCM_FRACTION_DENOMINATOR (y
)));
4420 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4423 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4427 scm_num2dbl (SCM a
, const char *why
)
4428 #define FUNC_NAME why
4431 return (double) SCM_INUM (a
);
4432 else if (SCM_BIGP (a
))
4434 double result
= mpz_get_d (SCM_I_BIG_MPZ (a
));
4435 scm_remember_upto_here_1 (a
);
4438 else if (SCM_REALP (a
))
4439 return (SCM_REAL_VALUE (a
));
4440 else if (SCM_FRACTIONP (a
))
4441 return scm_i_fraction2double (a
);
4443 SCM_WRONG_TYPE_ARG (SCM_ARGn
, a
);
4447 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4448 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4449 #define ALLOW_DIVIDE_BY_ZERO
4450 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4453 /* The code below for complex division is adapted from the GNU
4454 libstdc++, which adapted it from f2c's libF77, and is subject to
4457 /****************************************************************
4458 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4460 Permission to use, copy, modify, and distribute this software
4461 and its documentation for any purpose and without fee is hereby
4462 granted, provided that the above copyright notice appear in all
4463 copies and that both that the copyright notice and this
4464 permission notice and warranty disclaimer appear in supporting
4465 documentation, and that the names of AT&T Bell Laboratories or
4466 Bellcore or any of their entities not be used in advertising or
4467 publicity pertaining to distribution of the software without
4468 specific, written prior permission.
4470 AT&T and Bellcore disclaim all warranties with regard to this
4471 software, including all implied warranties of merchantability
4472 and fitness. In no event shall AT&T or Bellcore be liable for
4473 any special, indirect or consequential damages or any damages
4474 whatsoever resulting from loss of use, data or profits, whether
4475 in an action of contract, negligence or other tortious action,
4476 arising out of or in connection with the use or performance of
4478 ****************************************************************/
4480 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4481 /* Divide the first argument by the product of the remaining
4482 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4484 #define FUNC_NAME s_divide
4486 scm_i_divide (SCM x
, SCM y
, int inexact
)
4493 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4494 else if (SCM_INUMP (x
))
4496 long xx
= SCM_INUM (x
);
4497 if (xx
== 1 || xx
== -1)
4499 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4501 scm_num_overflow (s_divide
);
4506 return scm_make_real (1.0 / (double) xx
);
4507 else return scm_make_ratio (SCM_MAKINUM(1), x
);
4510 else if (SCM_BIGP (x
))
4513 return scm_make_real (1.0 / scm_i_big2dbl (x
));
4514 else return scm_make_ratio (SCM_MAKINUM(1), x
);
4516 else if (SCM_REALP (x
))
4518 double xx
= SCM_REAL_VALUE (x
);
4519 #ifndef ALLOW_DIVIDE_BY_ZERO
4521 scm_num_overflow (s_divide
);
4524 return scm_make_real (1.0 / xx
);
4526 else if (SCM_COMPLEXP (x
))
4528 double r
= SCM_COMPLEX_REAL (x
);
4529 double i
= SCM_COMPLEX_IMAG (x
);
4533 double d
= i
* (1.0 + t
* t
);
4534 return scm_make_complex (t
/ d
, -1.0 / d
);
4539 double d
= r
* (1.0 + t
* t
);
4540 return scm_make_complex (1.0 / d
, -t
/ d
);
4543 else if (SCM_FRACTIONP (x
))
4544 return scm_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4545 SCM_FRACTION_NUMERATOR (x
));
4547 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4552 long xx
= SCM_INUM (x
);
4555 long yy
= SCM_INUM (y
);
4558 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4559 scm_num_overflow (s_divide
);
4561 return scm_make_real ((double) xx
/ (double) yy
);
4564 else if (xx
% yy
!= 0)
4567 return scm_make_real ((double) xx
/ (double) yy
);
4568 else return scm_make_ratio (x
, y
);
4573 if (SCM_FIXABLE (z
))
4574 return SCM_MAKINUM (z
);
4576 return scm_i_long2big (z
);
4579 else if (SCM_BIGP (y
))
4582 return scm_make_real ((double) xx
/ scm_i_big2dbl (y
));
4583 else return scm_make_ratio (x
, y
);
4585 else if (SCM_REALP (y
))
4587 double yy
= SCM_REAL_VALUE (y
);
4588 #ifndef ALLOW_DIVIDE_BY_ZERO
4590 scm_num_overflow (s_divide
);
4593 return scm_make_real ((double) xx
/ yy
);
4595 else if (SCM_COMPLEXP (y
))
4598 complex_div
: /* y _must_ be a complex number */
4600 double r
= SCM_COMPLEX_REAL (y
);
4601 double i
= SCM_COMPLEX_IMAG (y
);
4605 double d
= i
* (1.0 + t
* t
);
4606 return scm_make_complex ((a
* t
) / d
, -a
/ d
);
4611 double d
= r
* (1.0 + t
* t
);
4612 return scm_make_complex (a
/ d
, -(a
* t
) / d
);
4616 else if (SCM_FRACTIONP (y
))
4617 /* a / b/c = ac / b */
4618 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4619 SCM_FRACTION_NUMERATOR (y
));
4621 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4623 else if (SCM_BIGP (x
))
4627 long int yy
= SCM_INUM (y
);
4630 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4631 scm_num_overflow (s_divide
);
4633 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4634 scm_remember_upto_here_1 (x
);
4635 return (sgn
== 0) ? scm_nan () : scm_inf ();
4642 /* FIXME: HMM, what are the relative performance issues here?
4643 We need to test. Is it faster on average to test
4644 divisible_p, then perform whichever operation, or is it
4645 faster to perform the integer div opportunistically and
4646 switch to real if there's a remainder? For now we take the
4647 middle ground: test, then if divisible, use the faster div
4650 long abs_yy
= yy
< 0 ? -yy
: yy
;
4651 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4655 SCM result
= scm_i_mkbig ();
4656 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4657 scm_remember_upto_here_1 (x
);
4659 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4660 return scm_i_normbig (result
);
4665 return scm_make_real (scm_i_big2dbl (x
) / (double) yy
);
4666 else return scm_make_ratio (x
, y
);
4670 else if (SCM_BIGP (y
))
4672 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4675 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4676 scm_num_overflow (s_divide
);
4678 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4679 scm_remember_upto_here_1 (x
);
4680 return (sgn
== 0) ? scm_nan () : scm_inf ();
4686 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4690 SCM result
= scm_i_mkbig ();
4691 mpz_divexact (SCM_I_BIG_MPZ (result
),
4694 scm_remember_upto_here_2 (x
, y
);
4695 return scm_i_normbig (result
);
4701 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4702 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4703 scm_remember_upto_here_2 (x
, y
);
4704 return scm_make_real (dbx
/ dby
);
4706 else return scm_make_ratio (x
, y
);
4710 else if (SCM_REALP (y
))
4712 double yy
= SCM_REAL_VALUE (y
);
4713 #ifndef ALLOW_DIVIDE_BY_ZERO
4715 scm_num_overflow (s_divide
);
4718 return scm_make_real (scm_i_big2dbl (x
) / yy
);
4720 else if (SCM_COMPLEXP (y
))
4722 a
= scm_i_big2dbl (x
);
4725 else if (SCM_FRACTIONP (y
))
4726 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4727 SCM_FRACTION_NUMERATOR (y
));
4729 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4731 else if (SCM_REALP (x
))
4733 double rx
= SCM_REAL_VALUE (x
);
4736 long int yy
= SCM_INUM (y
);
4737 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4739 scm_num_overflow (s_divide
);
4742 return scm_make_real (rx
/ (double) yy
);
4744 else if (SCM_BIGP (y
))
4746 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4747 scm_remember_upto_here_1 (y
);
4748 return scm_make_real (rx
/ dby
);
4750 else if (SCM_REALP (y
))
4752 double yy
= SCM_REAL_VALUE (y
);
4753 #ifndef ALLOW_DIVIDE_BY_ZERO
4755 scm_num_overflow (s_divide
);
4758 return scm_make_real (rx
/ yy
);
4760 else if (SCM_COMPLEXP (y
))
4765 else if (SCM_FRACTIONP (y
))
4766 return scm_make_real (rx
/ scm_i_fraction2double (y
));
4768 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4770 else if (SCM_COMPLEXP (x
))
4772 double rx
= SCM_COMPLEX_REAL (x
);
4773 double ix
= SCM_COMPLEX_IMAG (x
);
4776 long int yy
= SCM_INUM (y
);
4777 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4779 scm_num_overflow (s_divide
);
4784 return scm_make_complex (rx
/ d
, ix
/ d
);
4787 else if (SCM_BIGP (y
))
4789 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4790 scm_remember_upto_here_1 (y
);
4791 return scm_make_complex (rx
/ dby
, ix
/ dby
);
4793 else if (SCM_REALP (y
))
4795 double yy
= SCM_REAL_VALUE (y
);
4796 #ifndef ALLOW_DIVIDE_BY_ZERO
4798 scm_num_overflow (s_divide
);
4801 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4803 else if (SCM_COMPLEXP (y
))
4805 double ry
= SCM_COMPLEX_REAL (y
);
4806 double iy
= SCM_COMPLEX_IMAG (y
);
4810 double d
= iy
* (1.0 + t
* t
);
4811 return scm_make_complex ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4816 double d
= ry
* (1.0 + t
* t
);
4817 return scm_make_complex ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4820 else if (SCM_FRACTIONP (y
))
4822 double yy
= scm_i_fraction2double (y
);
4823 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4826 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4828 else if (SCM_FRACTIONP (x
))
4832 long int yy
= SCM_INUM (y
);
4833 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4835 scm_num_overflow (s_divide
);
4838 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4839 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4841 else if (SCM_BIGP (y
))
4843 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4844 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4846 else if (SCM_REALP (y
))
4848 double yy
= SCM_REAL_VALUE (y
);
4849 #ifndef ALLOW_DIVIDE_BY_ZERO
4851 scm_num_overflow (s_divide
);
4854 return scm_make_real (scm_i_fraction2double (x
) / yy
);
4856 else if (SCM_COMPLEXP (y
))
4858 a
= scm_i_fraction2double (x
);
4861 else if (SCM_FRACTIONP (y
))
4862 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4863 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4865 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4868 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4872 scm_divide (SCM x
, SCM y
)
4874 return scm_i_divide (x
, y
, 0);
4877 static SCM
scm_divide2real (SCM x
, SCM y
)
4879 return scm_i_divide (x
, y
, 1);
4885 scm_asinh (double x
)
4890 #define asinh scm_asinh
4891 return log (x
+ sqrt (x
* x
+ 1));
4894 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4895 /* "Return the inverse hyperbolic sine of @var{x}."
4900 scm_acosh (double x
)
4905 #define acosh scm_acosh
4906 return log (x
+ sqrt (x
* x
- 1));
4909 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4910 /* "Return the inverse hyperbolic cosine of @var{x}."
4915 scm_atanh (double x
)
4920 #define atanh scm_atanh
4921 return 0.5 * log ((1 + x
) / (1 - x
));
4924 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4925 /* "Return the inverse hyperbolic tangent of @var{x}."
4929 /* XXX - eventually, we should remove this definition of scm_round and
4930 rename scm_round_number to scm_round. Likewise for scm_truncate
4931 and scm_truncate_number.
4935 scm_truncate (double x
)
4940 #define trunc scm_truncate
4947 /* scm_round is done using floor(x+0.5) to round to nearest and with
4948 half-way case (ie. when x is an integer plus 0.5) going upwards. Then
4949 half-way cases are identified and adjusted down if the round-upwards
4950 didn't give the desired even integer.
4952 "plus_half == result" identifies a half-way case. If plus_half, which is
4953 x + 0.5, is an integer then x must be an integer plus 0.5.
4955 An odd "result" value is identified with result/2 != floor(result/2).
4956 This is done with plus_half, since that value is ready for use sooner in
4957 a pipelined cpu, and we're already requiring plus_half == result.
4959 Note however that we need to be careful when x is big and already an
4960 integer. In that case "x+0.5" may round to an adjacent integer, causing
4961 us to return such a value, incorrectly. For instance if the hardware is
4962 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
4963 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
4964 returned. Or if the hardware is in round-upwards mode, then other bigger
4965 values like say x == 2^128 will see x+0.5 rounding up to the next higher
4966 representable value, 2^128+2^76 (or whatever), again incorrect.
4968 These bad roundings of x+0.5 are avoided by testing at the start whether
4969 x is already an integer. If it is then clearly that's the desired result
4970 already. And if it's not then the exponent must be small enough to allow
4971 an 0.5 to be represented, and hence added without a bad rounding. */
4974 scm_round (double x
)
4976 double plus_half
, result
;
4981 plus_half
= x
+ 0.5;
4982 result
= floor (plus_half
);
4983 /* Adjust so that the scm_round is towards even. */
4984 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4989 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4991 "Round the number @var{x} towards zero.")
4992 #define FUNC_NAME s_scm_truncate_number
4994 if (SCM_FALSEP (scm_negative_p (x
)))
4995 return scm_floor (x
);
4997 return scm_ceiling (x
);
5001 static SCM exactly_one_half
;
5003 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5005 "Round the number @var{x} towards the nearest integer. "
5006 "When it is exactly halfway between two integers, "
5007 "round towards the even one.")
5008 #define FUNC_NAME s_scm_round_number
5010 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5011 SCM result
= scm_floor (plus_half
);
5012 /* Adjust so that the scm_round is towards even. */
5013 if (!SCM_FALSEP (scm_num_eq_p (plus_half
, result
))
5014 && !SCM_FALSEP (scm_odd_p (result
)))
5015 return scm_difference (result
, SCM_MAKINUM (1));
5021 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5023 "Round the number @var{x} towards minus infinity.")
5024 #define FUNC_NAME s_scm_floor
5026 if (SCM_INUMP (x
) || SCM_BIGP (x
))
5028 else if (SCM_REALP (x
))
5029 return scm_make_real (floor (SCM_REAL_VALUE (x
)));
5030 else if (SCM_FRACTIONP (x
))
5032 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5033 SCM_FRACTION_DENOMINATOR (x
));
5034 if (SCM_FALSEP (scm_negative_p (x
)))
5036 /* For positive x, rounding towards zero is correct. */
5041 /* For negative x, we need to return q-1 unless x is an
5042 integer. But fractions are never integer, per our
5044 return scm_difference (q
, SCM_MAKINUM (1));
5048 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5052 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5054 "Round the number @var{x} towards infinity.")
5055 #define FUNC_NAME s_scm_ceiling
5057 if (SCM_INUMP (x
) || SCM_BIGP (x
))
5059 else if (SCM_REALP (x
))
5060 return scm_make_real (ceil (SCM_REAL_VALUE (x
)));
5061 else if (SCM_FRACTIONP (x
))
5063 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5064 SCM_FRACTION_DENOMINATOR (x
));
5065 if (SCM_FALSEP (scm_positive_p (x
)))
5067 /* For negative x, rounding towards zero is correct. */
5072 /* For positive x, we need to return q+1 unless x is an
5073 integer. But fractions are never integer, per our
5075 return scm_sum (q
, SCM_MAKINUM (1));
5079 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5083 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5084 /* "Return the square root of the real number @var{x}."
5086 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5087 /* "Return the absolute value of the real number @var{x}."
5089 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5090 /* "Return the @var{x}th power of e."
5092 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5093 /* "Return the natural logarithm of the real number @var{x}."
5095 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5096 /* "Return the sine of the real number @var{x}."
5098 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5099 /* "Return the cosine of the real number @var{x}."
5101 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5102 /* "Return the tangent of the real number @var{x}."
5104 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5105 /* "Return the arc sine of the real number @var{x}."
5107 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5108 /* "Return the arc cosine of the real number @var{x}."
5110 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5111 /* "Return the arc tangent of the real number @var{x}."
5113 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5114 /* "Return the hyperbolic sine of the real number @var{x}."
5116 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5117 /* "Return the hyperbolic cosine of the real number @var{x}."
5119 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5120 /* "Return the hyperbolic tangent of the real number @var{x}."
5128 static void scm_two_doubles (SCM x
,
5130 const char *sstring
,
5134 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5137 xy
->x
= SCM_INUM (x
);
5138 else if (SCM_BIGP (x
))
5139 xy
->x
= scm_i_big2dbl (x
);
5140 else if (SCM_REALP (x
))
5141 xy
->x
= SCM_REAL_VALUE (x
);
5142 else if (SCM_FRACTIONP (x
))
5143 xy
->x
= scm_i_fraction2double (x
);
5145 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5148 xy
->y
= SCM_INUM (y
);
5149 else if (SCM_BIGP (y
))
5150 xy
->y
= scm_i_big2dbl (y
);
5151 else if (SCM_REALP (y
))
5152 xy
->y
= SCM_REAL_VALUE (y
);
5153 else if (SCM_FRACTIONP (y
))
5154 xy
->y
= scm_i_fraction2double (y
);
5156 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5160 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5162 "Return @var{x} raised to the power of @var{y}. This\n"
5163 "procedure does not accept complex arguments.")
5164 #define FUNC_NAME s_scm_sys_expt
5167 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5168 return scm_make_real (pow (xy
.x
, xy
.y
));
5173 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5175 "Return the arc tangent of the two arguments @var{x} and\n"
5176 "@var{y}. This is similar to calculating the arc tangent of\n"
5177 "@var{x} / @var{y}, except that the signs of both arguments\n"
5178 "are used to determine the quadrant of the result. This\n"
5179 "procedure does not accept complex arguments.")
5180 #define FUNC_NAME s_scm_sys_atan2
5183 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5184 return scm_make_real (atan2 (xy
.x
, xy
.y
));
5189 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5190 (SCM real
, SCM imaginary
),
5191 "Return a complex number constructed of the given @var{real} and\n"
5192 "@var{imaginary} parts.")
5193 #define FUNC_NAME s_scm_make_rectangular
5196 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5197 return scm_make_complex (xy
.x
, xy
.y
);
5203 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5205 "Return the complex number @var{x} * e^(i * @var{y}).")
5206 #define FUNC_NAME s_scm_make_polar
5210 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5212 sincos (xy
.y
, &s
, &c
);
5217 return scm_make_complex (xy
.x
* c
, xy
.x
* s
);
5222 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5223 /* "Return the real part of the number @var{z}."
5226 scm_real_part (SCM z
)
5230 else if (SCM_BIGP (z
))
5232 else if (SCM_REALP (z
))
5234 else if (SCM_COMPLEXP (z
))
5235 return scm_make_real (SCM_COMPLEX_REAL (z
));
5236 else if (SCM_FRACTIONP (z
))
5239 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5243 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5244 /* "Return the imaginary part of the number @var{z}."
5247 scm_imag_part (SCM z
)
5251 else if (SCM_BIGP (z
))
5253 else if (SCM_REALP (z
))
5255 else if (SCM_COMPLEXP (z
))
5256 return scm_make_real (SCM_COMPLEX_IMAG (z
));
5257 else if (SCM_FRACTIONP (z
))
5260 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5263 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5264 /* "Return the numerator of the number @var{z}."
5267 scm_numerator (SCM z
)
5271 else if (SCM_BIGP (z
))
5273 else if (SCM_FRACTIONP (z
))
5275 scm_i_fraction_reduce (z
);
5276 return SCM_FRACTION_NUMERATOR (z
);
5278 else if (SCM_REALP (z
))
5279 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5281 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5285 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5286 /* "Return the denominator of the number @var{z}."
5289 scm_denominator (SCM z
)
5292 return SCM_MAKINUM (1);
5293 else if (SCM_BIGP (z
))
5294 return SCM_MAKINUM (1);
5295 else if (SCM_FRACTIONP (z
))
5297 scm_i_fraction_reduce (z
);
5298 return SCM_FRACTION_DENOMINATOR (z
);
5300 else if (SCM_REALP (z
))
5301 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5303 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5306 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5307 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5308 * "@code{abs} for real arguments, but also allows complex numbers."
5311 scm_magnitude (SCM z
)
5315 long int zz
= SCM_INUM (z
);
5318 else if (SCM_POSFIXABLE (-zz
))
5319 return SCM_MAKINUM (-zz
);
5321 return scm_i_long2big (-zz
);
5323 else if (SCM_BIGP (z
))
5325 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5326 scm_remember_upto_here_1 (z
);
5328 return scm_i_clonebig (z
, 0);
5332 else if (SCM_REALP (z
))
5333 return scm_make_real (fabs (SCM_REAL_VALUE (z
)));
5334 else if (SCM_COMPLEXP (z
))
5335 return scm_make_real (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5336 else if (SCM_FRACTIONP (z
))
5338 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5340 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5341 SCM_FRACTION_DENOMINATOR (z
));
5344 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5348 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5349 /* "Return the angle of the complex number @var{z}."
5354 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5355 scm_flo0 to save allocating a new flonum with scm_make_real each time.
5356 But if atan2 follows the floating point rounding mode, then the value
5357 is not a constant. Maybe it'd be close enough though. */
5360 if (SCM_INUM (z
) >= 0)
5363 return scm_make_real (atan2 (0.0, -1.0));
5365 else if (SCM_BIGP (z
))
5367 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5368 scm_remember_upto_here_1 (z
);
5370 return scm_make_real (atan2 (0.0, -1.0));
5374 else if (SCM_REALP (z
))
5376 if (SCM_REAL_VALUE (z
) >= 0)
5379 return scm_make_real (atan2 (0.0, -1.0));
5381 else if (SCM_COMPLEXP (z
))
5382 return scm_make_real (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5383 else if (SCM_FRACTIONP (z
))
5385 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5387 else return scm_make_real (atan2 (0.0, -1.0));
5390 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5394 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5395 /* Convert the number @var{x} to its inexact representation.\n"
5398 scm_exact_to_inexact (SCM z
)
5401 return scm_make_real ((double) SCM_INUM (z
));
5402 else if (SCM_BIGP (z
))
5403 return scm_make_real (scm_i_big2dbl (z
));
5404 else if (SCM_FRACTIONP (z
))
5405 return scm_make_real (scm_i_fraction2double (z
));
5406 else if (SCM_INEXACTP (z
))
5409 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5413 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5415 "Return an exact number that is numerically closest to @var{z}.")
5416 #define FUNC_NAME s_scm_inexact_to_exact
5420 else if (SCM_BIGP (z
))
5422 else if (SCM_REALP (z
))
5424 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5425 SCM_OUT_OF_RANGE (1, z
);
5432 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5433 q
= scm_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5434 scm_i_mpz2num (mpq_denref (frac
)));
5436 /* When scm_make_ratio throws, we leak the memory allocated
5443 else if (SCM_FRACTIONP (z
))
5446 SCM_WRONG_TYPE_ARG (1, z
);
5450 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5452 "Return an exact number that is within @var{err} of @var{x}.")
5453 #define FUNC_NAME s_scm_rationalize
5457 else if (SCM_BIGP (x
))
5459 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5461 /* Use continued fractions to find closest ratio. All
5462 arithmetic is done with exact numbers.
5465 SCM ex
= scm_inexact_to_exact (x
);
5466 SCM int_part
= scm_floor (ex
);
5467 SCM tt
= SCM_MAKINUM (1);
5468 SCM a1
= SCM_MAKINUM (0), a2
= SCM_MAKINUM (1), a
= SCM_MAKINUM (0);
5469 SCM b1
= SCM_MAKINUM (1), b2
= SCM_MAKINUM (0), b
= SCM_MAKINUM (0);
5473 if (!SCM_FALSEP (scm_num_eq_p (ex
, int_part
)))
5476 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5477 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5479 /* We stop after a million iterations just to be absolutely sure
5480 that we don't go into an infinite loop. The process normally
5481 converges after less than a dozen iterations.
5484 err
= scm_abs (err
);
5485 while (++i
< 1000000)
5487 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5488 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5489 if (SCM_FALSEP (scm_zero_p (b
)) && /* b != 0 */
5491 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5492 err
))) /* abs(x-a/b) <= err */
5494 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5495 if (SCM_FALSEP (scm_exact_p (x
))
5496 || SCM_FALSEP (scm_exact_p (err
)))
5497 return scm_exact_to_inexact (res
);
5501 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5503 tt
= scm_floor (rx
); /* tt = floor (rx) */
5509 scm_num_overflow (s_scm_rationalize
);
5512 SCM_WRONG_TYPE_ARG (1, x
);
5516 /* if you need to change this, change test-num2integral.c as well */
5517 #if SCM_SIZEOF_LONG_LONG != 0
5519 # define ULLONG_MAX ((unsigned long long) (-1))
5520 # define LLONG_MAX ((long long) (ULLONG_MAX >> 1))
5521 # define LLONG_MIN (~LLONG_MAX)
5525 /* Parameters for creating integer conversion routines.
5527 Define the following preprocessor macros before including
5528 "libguile/num2integral.i.c":
5530 NUM2INTEGRAL - the name of the function for converting from a
5531 Scheme object to the integral type. This function will be
5532 defined when including "num2integral.i.c".
5534 INTEGRAL2NUM - the name of the function for converting from the
5535 integral type to a Scheme object. This function will be defined.
5537 INTEGRAL2BIG - the name of an internal function that createas a
5538 bignum from the integral type. This function will be defined.
5539 The name should start with "scm_i_".
5541 ITYPE - the name of the integral type.
5543 UNSIGNED - Define this to 1 when ITYPE is an unsigned type. Define
5546 UNSIGNED_ITYPE - the name of the the unsigned variant of the
5547 integral type. If you don't define this, it defaults to
5548 "unsigned ITYPE" for signed types and simply "ITYPE" for unsigned
5551 SIZEOF_ITYPE - an expression giving the size of the integral type
5552 in bytes. This expression must be computable by the
5553 preprocessor. (SIZEOF_FOO values are calculated by configure.in
5558 #define NUM2INTEGRAL scm_num2short
5559 #define INTEGRAL2NUM scm_short2num
5560 #define INTEGRAL2BIG scm_i_short2big
5563 #define SIZEOF_ITYPE SIZEOF_SHORT
5564 #include "libguile/num2integral.i.c"
5566 #define NUM2INTEGRAL scm_num2ushort
5567 #define INTEGRAL2NUM scm_ushort2num
5568 #define INTEGRAL2BIG scm_i_ushort2big
5570 #define ITYPE unsigned short
5571 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_SHORT
5572 #include "libguile/num2integral.i.c"
5574 #define NUM2INTEGRAL scm_num2int
5575 #define INTEGRAL2NUM scm_int2num
5576 #define INTEGRAL2BIG scm_i_int2big
5579 #define SIZEOF_ITYPE SIZEOF_INT
5580 #include "libguile/num2integral.i.c"
5582 #define NUM2INTEGRAL scm_num2uint
5583 #define INTEGRAL2NUM scm_uint2num
5584 #define INTEGRAL2BIG scm_i_uint2big
5586 #define ITYPE unsigned int
5587 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_INT
5588 #include "libguile/num2integral.i.c"
5590 #define NUM2INTEGRAL scm_num2long
5591 #define INTEGRAL2NUM scm_long2num
5592 #define INTEGRAL2BIG scm_i_long2big
5595 #define SIZEOF_ITYPE SIZEOF_LONG
5596 #include "libguile/num2integral.i.c"
5598 #define NUM2INTEGRAL scm_num2ulong
5599 #define INTEGRAL2NUM scm_ulong2num
5600 #define INTEGRAL2BIG scm_i_ulong2big
5602 #define ITYPE unsigned long
5603 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG
5604 #include "libguile/num2integral.i.c"
5606 #define NUM2INTEGRAL scm_num2ptrdiff
5607 #define INTEGRAL2NUM scm_ptrdiff2num
5608 #define INTEGRAL2BIG scm_i_ptrdiff2big
5610 #define ITYPE scm_t_ptrdiff
5611 #define UNSIGNED_ITYPE size_t
5612 #define SIZEOF_ITYPE SCM_SIZEOF_SCM_T_PTRDIFF
5613 #include "libguile/num2integral.i.c"
5615 #define NUM2INTEGRAL scm_num2size
5616 #define INTEGRAL2NUM scm_size2num
5617 #define INTEGRAL2BIG scm_i_size2big
5619 #define ITYPE size_t
5620 #define SIZEOF_ITYPE SIZEOF_SIZE_T
5621 #include "libguile/num2integral.i.c"
5623 #if SCM_SIZEOF_LONG_LONG != 0
5625 #ifndef ULONG_LONG_MAX
5626 #define ULONG_LONG_MAX (~0ULL)
5629 #define NUM2INTEGRAL scm_num2long_long
5630 #define INTEGRAL2NUM scm_long_long2num
5631 #define INTEGRAL2BIG scm_i_long_long2big
5633 #define ITYPE long long
5634 #define SIZEOF_ITYPE SIZEOF_LONG_LONG
5635 #include "libguile/num2integral.i.c"
5637 #define NUM2INTEGRAL scm_num2ulong_long
5638 #define INTEGRAL2NUM scm_ulong_long2num
5639 #define INTEGRAL2BIG scm_i_ulong_long2big
5641 #define ITYPE unsigned long long
5642 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG_LONG
5643 #include "libguile/num2integral.i.c"
5645 #endif /* SCM_SIZEOF_LONG_LONG != 0 */
5647 #define NUM2FLOAT scm_num2float
5648 #define FLOAT2NUM scm_float2num
5650 #include "libguile/num2float.i.c"
5652 #define NUM2FLOAT scm_num2double
5653 #define FLOAT2NUM scm_double2num
5654 #define FTYPE double
5655 #include "libguile/num2float.i.c"
5660 #define SIZE_MAX ((size_t) (-1))
5663 #define PTRDIFF_MIN \
5664 ((scm_t_ptrdiff) ((scm_t_ptrdiff) 1 \
5665 << ((sizeof (scm_t_ptrdiff) * SCM_CHAR_BIT) - 1)))
5668 #define PTRDIFF_MAX (~ PTRDIFF_MIN)
5671 #define CHECK(type, v) \
5674 if ((v) != scm_num2##type (scm_##type##2num (v), 1, "check_sanity")) \
5694 CHECK (ptrdiff
, -1);
5696 CHECK (short, SHRT_MAX
);
5697 CHECK (short, SHRT_MIN
);
5698 CHECK (ushort
, USHRT_MAX
);
5699 CHECK (int, INT_MAX
);
5700 CHECK (int, INT_MIN
);
5701 CHECK (uint
, UINT_MAX
);
5702 CHECK (long, LONG_MAX
);
5703 CHECK (long, LONG_MIN
);
5704 CHECK (ulong
, ULONG_MAX
);
5705 CHECK (size
, SIZE_MAX
);
5706 CHECK (ptrdiff
, PTRDIFF_MAX
);
5707 CHECK (ptrdiff
, PTRDIFF_MIN
);
5709 #if SCM_SIZEOF_LONG_LONG != 0
5710 CHECK (long_long
, 0LL);
5711 CHECK (ulong_long
, 0ULL);
5712 CHECK (long_long
, -1LL);
5713 CHECK (long_long
, LLONG_MAX
);
5714 CHECK (long_long
, LLONG_MIN
);
5715 CHECK (ulong_long
, ULLONG_MAX
);
5722 scm_internal_catch (SCM_BOOL_T, check_body, &data, check_handler, &data); \
5723 if (!SCM_FALSEP (data)) abort();
5726 check_body (void *data
)
5728 SCM num
= *(SCM
*) data
;
5729 scm_num2ulong (num
, 1, NULL
);
5731 return SCM_UNSPECIFIED
;
5735 check_handler (void *data
, SCM tag
, SCM throw_args
)
5737 SCM
*num
= (SCM
*) data
;
5740 return SCM_UNSPECIFIED
;
5743 SCM_DEFINE (scm_sys_check_number_conversions
, "%check-number-conversions", 0, 0, 0,
5745 "Number conversion sanity checking.")
5746 #define FUNC_NAME s_scm_sys_check_number_conversions
5748 SCM data
= SCM_MAKINUM (-1);
5750 data
= scm_int2num (INT_MIN
);
5752 data
= scm_ulong2num (ULONG_MAX
);
5753 data
= scm_difference (SCM_INUM0
, data
);
5755 data
= scm_ulong2num (ULONG_MAX
);
5756 data
= scm_sum (SCM_MAKINUM (1), data
); data
= scm_difference (SCM_INUM0
, data
);
5758 data
= scm_int2num (-10000); data
= scm_product (data
, data
); data
= scm_product (data
, data
);
5761 return SCM_UNSPECIFIED
;
5772 mpz_init_set_si (z_negative_one
, -1);
5774 /* It may be possible to tune the performance of some algorithms by using
5775 * the following constants to avoid the creation of bignums. Please, before
5776 * using these values, remember the two rules of program optimization:
5777 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5778 scm_c_define ("most-positive-fixnum",
5779 SCM_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5780 scm_c_define ("most-negative-fixnum",
5781 SCM_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5783 scm_add_feature ("complex");
5784 scm_add_feature ("inexact");
5785 scm_flo0
= scm_make_real (0.0);
5787 /* determine floating point precision */
5788 for(i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
5790 init_dblprec(&scm_dblprec
[i
-2],i
);
5791 init_fx_radix(fx_per_radix
[i
-2],i
);
5794 /* hard code precision for base 10 if the preprocessor tells us to... */
5795 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5802 exactly_one_half
= scm_permanent_object (scm_divide (SCM_MAKINUM (1),
5804 #include "libguile/numbers.x"