1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003 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 (10+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)
130 #if defined (HAVE_ISINF)
132 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
133 return (! (finite (x
) || isnan (x
)));
142 #if defined (HAVE_ISNAN)
151 static SCM abs_most_negative_fixnum
;
152 static mpz_t z_negative_one
;
156 static const char s_bignum
[] = "bignum";
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 if (SCM_EQ_P (numerator
, SCM_INUM0
))
362 if (SCM_INUMP (denominator
))
365 x
= SCM_INUM (numerator
);
366 y
= SCM_INUM (denominator
);
368 return SCM_MAKINUM(1);
370 return SCM_MAKINUM (x
/ y
);
373 else if (SCM_BIGP (numerator
))
375 if (SCM_INUMP (denominator
))
377 long yy
= SCM_INUM (denominator
);
378 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
379 return scm_divide (numerator
, denominator
);
383 if (SCM_EQ_P (numerator
, denominator
))
384 return SCM_MAKINUM(1);
385 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
386 SCM_I_BIG_MPZ (denominator
)))
387 return scm_divide(numerator
, denominator
);
391 /* No, it's a proper fraction.
393 return scm_double_cell (scm_tc16_fraction
,
394 SCM_UNPACK (numerator
),
395 SCM_UNPACK (denominator
), 0);
399 static void scm_i_fraction_reduce (SCM z
)
401 if (!(SCM_FRACTION_REDUCED (z
)))
404 divisor
= scm_gcd (SCM_FRACTION_NUMERATOR (z
), SCM_FRACTION_DENOMINATOR (z
));
405 if (!(SCM_EQ_P (divisor
, SCM_MAKINUM(1))))
408 SCM_FRACTION_SET_NUMERATOR (z
, scm_divide (SCM_FRACTION_NUMERATOR (z
), divisor
));
409 SCM_FRACTION_SET_DENOMINATOR (z
, scm_divide (SCM_FRACTION_DENOMINATOR (z
), divisor
));
411 SCM_FRACTION_REDUCED_SET (z
);
416 scm_i_fraction2double (SCM z
)
418 return scm_num2dbl (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
419 SCM_FRACTION_DENOMINATOR (z
)),
423 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
425 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
427 #define FUNC_NAME s_scm_exact_p
433 if (SCM_FRACTIONP (x
))
437 SCM_WRONG_TYPE_ARG (1, x
);
442 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
444 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
446 #define FUNC_NAME s_scm_odd_p
450 long val
= SCM_INUM (n
);
451 return SCM_BOOL ((val
& 1L) != 0);
453 else if (SCM_BIGP (n
))
455 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
456 scm_remember_upto_here_1 (n
);
457 return SCM_BOOL (odd_p
);
459 else if (!SCM_FALSEP (scm_inf_p (n
)))
461 else if (SCM_REALP (n
))
463 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
469 SCM_WRONG_TYPE_ARG (1, n
);
472 SCM_WRONG_TYPE_ARG (1, n
);
477 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
479 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
481 #define FUNC_NAME s_scm_even_p
485 long val
= SCM_INUM (n
);
486 return SCM_BOOL ((val
& 1L) == 0);
488 else if (SCM_BIGP (n
))
490 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
491 scm_remember_upto_here_1 (n
);
492 return SCM_BOOL (even_p
);
494 else if (!SCM_FALSEP (scm_inf_p (n
)))
496 else if (SCM_REALP (n
))
498 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
504 SCM_WRONG_TYPE_ARG (1, n
);
507 SCM_WRONG_TYPE_ARG (1, n
);
511 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
513 "Return @code{#t} if @var{n} is infinite, @code{#f}\n"
515 #define FUNC_NAME s_scm_inf_p
518 return SCM_BOOL (xisinf (SCM_REAL_VALUE (n
)));
519 else if (SCM_COMPLEXP (n
))
520 return SCM_BOOL (xisinf (SCM_COMPLEX_REAL (n
))
521 || xisinf (SCM_COMPLEX_IMAG (n
)));
527 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
529 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
531 #define FUNC_NAME s_scm_nan_p
534 return SCM_BOOL (xisnan (SCM_REAL_VALUE (n
)));
535 else if (SCM_COMPLEXP (n
))
536 return SCM_BOOL (xisnan (SCM_COMPLEX_REAL (n
))
537 || xisnan (SCM_COMPLEX_IMAG (n
)));
543 /* Guile's idea of infinity. */
544 static double guile_Inf
;
546 /* Guile's idea of not a number. */
547 static double guile_NaN
;
550 guile_ieee_init (void)
552 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
554 /* Some version of gcc on some old version of Linux used to crash when
555 trying to make Inf and NaN. */
559 guile_Inf
= 1.0 / (tmp
- tmp
);
560 #elif defined (__alpha__) && ! defined (linux)
561 extern unsigned int DINFINITY
[2];
562 guile_Inf
= (*(X_CAST(double *, DINFINITY
)));
569 if (guile_Inf
== tmp
)
577 #if defined (HAVE_ISNAN)
579 #if defined (__alpha__) && ! defined (linux)
580 extern unsigned int DQNAN
[2];
581 guile_NaN
= (*(X_CAST(double *, DQNAN
)));
583 guile_NaN
= guile_Inf
/ guile_Inf
;
589 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
592 #define FUNC_NAME s_scm_inf
594 static int initialized
= 0;
600 return scm_make_real (guile_Inf
);
604 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
607 #define FUNC_NAME s_scm_nan
609 static int initialized
= 0;
615 return scm_make_real (guile_NaN
);
620 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
622 "Return the absolute value of @var{x}.")
627 long int xx
= SCM_INUM (x
);
630 else if (SCM_POSFIXABLE (-xx
))
631 return SCM_MAKINUM (-xx
);
633 return scm_i_long2big (-xx
);
635 else if (SCM_BIGP (x
))
637 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
639 return scm_i_clonebig (x
, 0);
643 else if (SCM_REALP (x
))
644 return scm_make_real (fabs (SCM_REAL_VALUE (x
)));
645 else if (SCM_FRACTIONP (x
))
647 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
649 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
650 SCM_FRACTION_DENOMINATOR (x
));
653 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
658 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
659 /* "Return the quotient of the numbers @var{x} and @var{y}."
662 scm_quotient (SCM x
, SCM y
)
666 long xx
= SCM_INUM (x
);
669 long yy
= SCM_INUM (y
);
671 scm_num_overflow (s_quotient
);
676 return SCM_MAKINUM (z
);
678 return scm_i_long2big (z
);
681 else if (SCM_BIGP (y
))
683 if ((SCM_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
684 && (scm_i_bigcmp (abs_most_negative_fixnum
, y
) == 0))
685 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
686 return SCM_MAKINUM (-1);
688 return SCM_MAKINUM (0);
691 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
693 else if (SCM_BIGP (x
))
697 long yy
= SCM_INUM (y
);
699 scm_num_overflow (s_quotient
);
704 SCM result
= scm_i_mkbig ();
707 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
710 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
713 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
714 scm_remember_upto_here_1 (x
);
715 return scm_i_normbig (result
);
718 else if (SCM_BIGP (y
))
720 SCM result
= scm_i_mkbig ();
721 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
724 scm_remember_upto_here_2 (x
, y
);
725 return scm_i_normbig (result
);
728 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
731 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
734 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
735 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
737 * "(remainder 13 4) @result{} 1\n"
738 * "(remainder -13 4) @result{} -1\n"
742 scm_remainder (SCM x
, SCM y
)
748 long yy
= SCM_INUM (y
);
750 scm_num_overflow (s_remainder
);
753 long z
= SCM_INUM (x
) % yy
;
754 return SCM_MAKINUM (z
);
757 else if (SCM_BIGP (y
))
759 if ((SCM_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
760 && (scm_i_bigcmp (abs_most_negative_fixnum
, y
) == 0))
761 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
762 return SCM_MAKINUM (0);
767 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
769 else if (SCM_BIGP (x
))
773 long yy
= SCM_INUM (y
);
775 scm_num_overflow (s_remainder
);
778 SCM result
= scm_i_mkbig ();
781 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
782 scm_remember_upto_here_1 (x
);
783 return scm_i_normbig (result
);
786 else if (SCM_BIGP (y
))
788 SCM result
= scm_i_mkbig ();
789 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
792 scm_remember_upto_here_2 (x
, y
);
793 return scm_i_normbig (result
);
796 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
799 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
803 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
804 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
806 * "(modulo 13 4) @result{} 1\n"
807 * "(modulo -13 4) @result{} 3\n"
811 scm_modulo (SCM x
, SCM y
)
815 long xx
= SCM_INUM (x
);
818 long yy
= SCM_INUM (y
);
820 scm_num_overflow (s_modulo
);
823 /* FIXME: I think this may be a bug on some arches -- results
824 of % with negative second arg are undefined... */
842 return SCM_MAKINUM (result
);
845 else if (SCM_BIGP (y
))
847 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
850 scm_num_overflow (s_modulo
);
858 SCM pos_y
= scm_i_clonebig (y
, 0);
859 /* do this after the last scm_op */
860 mpz_init_set_si (z_x
, xx
);
861 result
= pos_y
; /* re-use this bignum */
862 mpz_mod (SCM_I_BIG_MPZ (result
),
864 SCM_I_BIG_MPZ (pos_y
));
865 scm_remember_upto_here_1 (pos_y
);
869 result
= scm_i_mkbig ();
870 /* do this after the last scm_op */
871 mpz_init_set_si (z_x
, xx
);
872 mpz_mod (SCM_I_BIG_MPZ (result
),
875 scm_remember_upto_here_1 (y
);
878 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
879 mpz_add (SCM_I_BIG_MPZ (result
),
881 SCM_I_BIG_MPZ (result
));
882 scm_remember_upto_here_1 (y
);
883 /* and do this before the next one */
885 return scm_i_normbig (result
);
889 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
891 else if (SCM_BIGP (x
))
895 long yy
= SCM_INUM (y
);
897 scm_num_overflow (s_modulo
);
900 SCM result
= scm_i_mkbig ();
901 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
903 (yy
< 0) ? - yy
: yy
);
904 scm_remember_upto_here_1 (x
);
905 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
906 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
907 SCM_I_BIG_MPZ (result
),
909 return scm_i_normbig (result
);
912 else if (SCM_BIGP (y
))
914 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
916 scm_num_overflow (s_modulo
);
919 SCM result
= scm_i_mkbig ();
920 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
921 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
922 mpz_mod (SCM_I_BIG_MPZ (result
),
924 SCM_I_BIG_MPZ (pos_y
));
926 scm_remember_upto_here_1 (x
);
927 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
928 mpz_add (SCM_I_BIG_MPZ (result
),
930 SCM_I_BIG_MPZ (result
));
931 scm_remember_upto_here_2 (y
, pos_y
);
932 return scm_i_normbig (result
);
936 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
939 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
942 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
943 /* "Return the greatest common divisor of all arguments.\n"
944 * "If called without arguments, 0 is returned."
947 scm_gcd (SCM x
, SCM y
)
950 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
956 long xx
= SCM_INUM (x
);
957 long yy
= SCM_INUM (y
);
958 long u
= xx
< 0 ? -xx
: xx
;
959 long v
= yy
< 0 ? -yy
: yy
;
969 /* Determine a common factor 2^k */
970 while (!(1 & (u
| v
)))
976 /* Now, any factor 2^n can be eliminated */
996 return (SCM_POSFIXABLE (result
)
997 ? SCM_MAKINUM (result
)
998 : scm_i_long2big (result
));
1000 else if (SCM_BIGP (y
))
1002 SCM result
= scm_i_mkbig ();
1003 SCM mx
= scm_i_mkbig ();
1004 mpz_set_si (SCM_I_BIG_MPZ (mx
), SCM_INUM (x
));
1005 scm_remember_upto_here_1 (x
);
1006 mpz_gcd (SCM_I_BIG_MPZ (result
),
1009 scm_remember_upto_here_2 (mx
, y
);
1010 return scm_i_normbig (result
);
1013 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1015 else if (SCM_BIGP (x
))
1019 unsigned long result
;
1020 long yy
= SCM_INUM (y
);
1025 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1026 scm_remember_upto_here_1 (x
);
1027 return (SCM_POSFIXABLE (result
)
1028 ? SCM_MAKINUM (result
)
1029 : scm_ulong2num (result
));
1031 else if (SCM_BIGP (y
))
1033 SCM result
= scm_i_mkbig ();
1034 mpz_gcd (SCM_I_BIG_MPZ (result
),
1037 scm_remember_upto_here_2 (x
, y
);
1038 return scm_i_normbig (result
);
1041 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1044 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1047 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1048 /* "Return the least common multiple of the arguments.\n"
1049 * "If called without arguments, 1 is returned."
1052 scm_lcm (SCM n1
, SCM n2
)
1054 if (SCM_UNBNDP (n2
))
1056 if (SCM_UNBNDP (n1
))
1057 return SCM_MAKINUM (1L);
1058 n2
= SCM_MAKINUM (1L);
1061 SCM_GASSERT2 (SCM_INUMP (n1
) || SCM_BIGP (n1
),
1062 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1063 SCM_GASSERT2 (SCM_INUMP (n2
) || SCM_BIGP (n2
),
1064 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1070 SCM d
= scm_gcd (n1
, n2
);
1071 if (SCM_EQ_P (d
, SCM_INUM0
))
1074 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1078 /* inum n1, big n2 */
1081 SCM result
= scm_i_mkbig ();
1082 long nn1
= SCM_INUM (n1
);
1083 if (nn1
== 0) return SCM_INUM0
;
1084 if (nn1
< 0) nn1
= - nn1
;
1085 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1086 scm_remember_upto_here_1 (n2
);
1101 SCM result
= scm_i_mkbig ();
1102 mpz_lcm(SCM_I_BIG_MPZ (result
),
1104 SCM_I_BIG_MPZ (n2
));
1105 scm_remember_upto_here_2(n1
, n2
);
1106 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1112 #ifndef scm_long2num
1113 #define SCM_LOGOP_RETURN(x) scm_ulong2num(x)
1115 #define SCM_LOGOP_RETURN(x) SCM_MAKINUM(x)
1118 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1123 + + + x (map digit:logand X Y)
1124 + - + x (map digit:logand X (lognot (+ -1 Y)))
1125 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1126 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1131 + + + (map digit:logior X Y)
1132 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1133 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1134 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1139 + + + (map digit:logxor X Y)
1140 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1141 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1142 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1147 + + (any digit:logand X Y)
1148 + - (any digit:logand X (lognot (+ -1 Y)))
1149 - + (any digit:logand (lognot (+ -1 X)) Y)
1154 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1156 "Return the bitwise AND of the integer arguments.\n\n"
1158 "(logand) @result{} -1\n"
1159 "(logand 7) @result{} 7\n"
1160 "(logand #b111 #b011 #\b001) @result{} 1\n"
1162 #define FUNC_NAME s_scm_logand
1166 if (SCM_UNBNDP (n2
))
1168 if (SCM_UNBNDP (n1
))
1169 return SCM_MAKINUM (-1);
1170 else if (!SCM_NUMBERP (n1
))
1171 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1172 else if (SCM_NUMBERP (n1
))
1175 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1180 nn1
= SCM_INUM (n1
);
1183 long nn2
= SCM_INUM (n2
);
1184 return SCM_MAKINUM (nn1
& nn2
);
1186 else if SCM_BIGP (n2
)
1192 SCM result_z
= scm_i_mkbig ();
1194 mpz_init_set_si (nn1_z
, nn1
);
1195 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1196 scm_remember_upto_here_1 (n2
);
1198 return scm_i_normbig (result_z
);
1202 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1204 else if (SCM_BIGP (n1
))
1209 nn1
= SCM_INUM (n1
);
1212 else if (SCM_BIGP (n2
))
1214 SCM result_z
= scm_i_mkbig ();
1215 mpz_and (SCM_I_BIG_MPZ (result_z
),
1217 SCM_I_BIG_MPZ (n2
));
1218 scm_remember_upto_here_2 (n1
, n2
);
1219 return scm_i_normbig (result_z
);
1222 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1225 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1230 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1232 "Return the bitwise OR of the integer arguments.\n\n"
1234 "(logior) @result{} 0\n"
1235 "(logior 7) @result{} 7\n"
1236 "(logior #b000 #b001 #b011) @result{} 3\n"
1238 #define FUNC_NAME s_scm_logior
1242 if (SCM_UNBNDP (n2
))
1244 if (SCM_UNBNDP (n1
))
1246 else if (SCM_NUMBERP (n1
))
1249 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1254 nn1
= SCM_INUM (n1
);
1257 long nn2
= SCM_INUM (n2
);
1258 return SCM_MAKINUM (nn1
| nn2
);
1260 else if (SCM_BIGP (n2
))
1266 SCM result_z
= scm_i_mkbig ();
1268 mpz_init_set_si (nn1_z
, nn1
);
1269 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1270 scm_remember_upto_here_1 (n2
);
1276 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1278 else if (SCM_BIGP (n1
))
1283 nn1
= SCM_INUM (n1
);
1286 else if (SCM_BIGP (n2
))
1288 SCM result_z
= scm_i_mkbig ();
1289 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1291 SCM_I_BIG_MPZ (n2
));
1292 scm_remember_upto_here_2 (n1
, n2
);
1296 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1299 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1304 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1306 "Return the bitwise XOR of the integer arguments. A bit is\n"
1307 "set in the result if it is set in an odd number of arguments.\n"
1309 "(logxor) @result{} 0\n"
1310 "(logxor 7) @result{} 7\n"
1311 "(logxor #b000 #b001 #b011) @result{} 2\n"
1312 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1314 #define FUNC_NAME s_scm_logxor
1318 if (SCM_UNBNDP (n2
))
1320 if (SCM_UNBNDP (n1
))
1322 else if (SCM_NUMBERP (n1
))
1325 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1330 nn1
= SCM_INUM (n1
);
1333 long nn2
= SCM_INUM (n2
);
1334 return SCM_MAKINUM (nn1
^ nn2
);
1336 else if (SCM_BIGP (n2
))
1340 SCM result_z
= scm_i_mkbig ();
1342 mpz_init_set_si (nn1_z
, nn1
);
1343 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1344 scm_remember_upto_here_1 (n2
);
1346 return scm_i_normbig (result_z
);
1350 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1352 else if (SCM_BIGP (n1
))
1357 nn1
= SCM_INUM (n1
);
1360 else if (SCM_BIGP (n2
))
1362 SCM result_z
= scm_i_mkbig ();
1363 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1365 SCM_I_BIG_MPZ (n2
));
1366 scm_remember_upto_here_2 (n1
, n2
);
1367 return scm_i_normbig (result_z
);
1370 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1373 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1378 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1381 "(logtest j k) @equiv{} (not (zero? (logand j k)))\n\n"
1382 "(logtest #b0100 #b1011) @result{} #f\n"
1383 "(logtest #b0100 #b0111) @result{} #t\n"
1385 #define FUNC_NAME s_scm_logtest
1394 long nk
= SCM_INUM (k
);
1395 return SCM_BOOL (nj
& nk
);
1397 else if (SCM_BIGP (k
))
1405 mpz_init_set_si (nj_z
, nj
);
1406 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1407 scm_remember_upto_here_1 (k
);
1408 result
= SCM_BOOL (mpz_sgn (nj_z
) != 0);
1414 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1416 else if (SCM_BIGP (j
))
1424 else if (SCM_BIGP (k
))
1428 mpz_init (result_z
);
1432 scm_remember_upto_here_2 (j
, k
);
1433 result
= SCM_BOOL (mpz_sgn (result_z
) != 0);
1434 mpz_clear (result_z
);
1438 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1441 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1446 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1449 "(logbit? index j) @equiv{} (logtest (integer-expt 2 index) j)\n\n"
1450 "(logbit? 0 #b1101) @result{} #t\n"
1451 "(logbit? 1 #b1101) @result{} #f\n"
1452 "(logbit? 2 #b1101) @result{} #t\n"
1453 "(logbit? 3 #b1101) @result{} #t\n"
1454 "(logbit? 4 #b1101) @result{} #f\n"
1456 #define FUNC_NAME s_scm_logbit_p
1458 unsigned long int iindex
;
1460 SCM_VALIDATE_INUM_MIN (SCM_ARG1
, index
, 0);
1461 iindex
= (unsigned long int) SCM_INUM (index
);
1464 return SCM_BOOL ((1L << iindex
) & SCM_INUM (j
));
1465 else if (SCM_BIGP (j
))
1467 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1468 scm_remember_upto_here_1 (j
);
1469 return SCM_BOOL (val
);
1472 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1477 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1479 "Return the integer which is the ones-complement of the integer\n"
1483 "(number->string (lognot #b10000000) 2)\n"
1484 " @result{} \"-10000001\"\n"
1485 "(number->string (lognot #b0) 2)\n"
1486 " @result{} \"-1\"\n"
1488 #define FUNC_NAME s_scm_lognot
1490 if (SCM_INUMP (n
)) {
1491 /* No overflow here, just need to toggle all the bits making up the inum.
1492 Enhancement: No need to strip the tag and add it back, could just xor
1493 a block of 1 bits, if that worked with the various debug versions of
1495 return SCM_MAKINUM (~ SCM_INUM (n
));
1497 } else if (SCM_BIGP (n
)) {
1498 SCM result
= scm_i_mkbig ();
1499 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1500 scm_remember_upto_here_1 (n
);
1504 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1509 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1511 "Return @var{n} raised to the non-negative integer exponent\n"
1515 "(integer-expt 2 5)\n"
1517 "(integer-expt -3 3)\n"
1520 #define FUNC_NAME s_scm_integer_expt
1523 SCM z_i2
= SCM_BOOL_F
;
1525 SCM acc
= SCM_MAKINUM (1L);
1527 /* 0^0 == 1 according to R5RS */
1528 if (SCM_EQ_P (n
, SCM_INUM0
) || SCM_EQ_P (n
, acc
))
1529 return SCM_FALSEP (scm_zero_p(k
)) ? n
: acc
;
1530 else if (SCM_EQ_P (n
, SCM_MAKINUM (-1L)))
1531 return SCM_FALSEP (scm_even_p (k
)) ? n
: acc
;
1535 else if (SCM_BIGP (k
))
1537 z_i2
= scm_i_clonebig (k
, 1);
1538 scm_remember_upto_here_1 (k
);
1541 else if (SCM_REALP (k
))
1543 double r
= SCM_REAL_VALUE (k
);
1545 SCM_WRONG_TYPE_ARG (2, k
);
1546 if ((r
> SCM_MOST_POSITIVE_FIXNUM
) || (r
< SCM_MOST_NEGATIVE_FIXNUM
))
1548 z_i2
= scm_i_mkbig ();
1549 mpz_set_d (SCM_I_BIG_MPZ (z_i2
), r
);
1558 SCM_WRONG_TYPE_ARG (2, k
);
1562 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1564 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1565 n
= scm_divide (n
, SCM_UNDEFINED
);
1569 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1573 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1575 return scm_product (acc
, n
);
1577 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1578 acc
= scm_product (acc
, n
);
1579 n
= scm_product (n
, n
);
1580 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1588 n
= scm_divide (n
, SCM_UNDEFINED
);
1595 return scm_product (acc
, n
);
1597 acc
= scm_product (acc
, n
);
1598 n
= scm_product (n
, n
);
1605 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1607 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1608 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1610 "This is effectively a multiplication by 2^@var{cnt}}, and when\n"
1611 "@var{cnt} is negative it's a division, rounded towards negative\n"
1612 "infinity. (Note that this is not the same rounding as\n"
1613 "@code{quotient} does.)\n"
1615 "With @var{n} viewed as an infinite precision twos complement,\n"
1616 "@code{ash} means a left shift introducing zero bits, or a right\n"
1617 "shift dropping bits.\n"
1620 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1621 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1623 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1624 "(ash -23 -2) @result{} -6\n"
1626 #define FUNC_NAME s_scm_ash
1630 SCM_VALIDATE_INUM (2, cnt
);
1632 bits_to_shift
= SCM_INUM (cnt
);
1634 if (bits_to_shift
< 0)
1636 /* Shift right by abs(cnt) bits. This is realized as a division
1637 by div:=2^abs(cnt). However, to guarantee the floor
1638 rounding, negative values require some special treatment.
1640 SCM div
= scm_integer_expt (SCM_MAKINUM (2),
1641 SCM_MAKINUM (-bits_to_shift
));
1643 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1644 if (SCM_FALSEP (scm_negative_p (n
)))
1645 return scm_quotient (n
, div
);
1647 return scm_sum (SCM_MAKINUM (-1L),
1648 scm_quotient (scm_sum (SCM_MAKINUM (1L), n
), div
));
1651 /* Shift left is done by multiplication with 2^CNT */
1652 return scm_product (n
, scm_integer_expt (SCM_MAKINUM (2), cnt
));
1657 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1658 (SCM n
, SCM start
, SCM end
),
1659 "Return the integer composed of the @var{start} (inclusive)\n"
1660 "through @var{end} (exclusive) bits of @var{n}. The\n"
1661 "@var{start}th bit becomes the 0-th bit in the result.\n"
1664 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1665 " @result{} \"1010\"\n"
1666 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1667 " @result{} \"10110\"\n"
1669 #define FUNC_NAME s_scm_bit_extract
1671 unsigned long int istart
, iend
;
1672 SCM_VALIDATE_INUM_MIN_COPY (2, start
,0, istart
);
1673 SCM_VALIDATE_INUM_MIN_COPY (3, end
, 0, iend
);
1674 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1678 long int in
= SCM_INUM (n
);
1679 unsigned long int bits
= iend
- istart
;
1681 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1683 /* Since we emulate two's complement encoded numbers, this
1684 * special case requires us to produce a result that has
1685 * more bits than can be stored in a fixnum. Thus, we fall
1686 * back to the more general algorithm that is used for
1692 if (istart
< SCM_I_FIXNUM_BIT
)
1695 if (bits
< SCM_I_FIXNUM_BIT
)
1696 return SCM_MAKINUM (in
& ((1L << bits
) - 1));
1697 else /* we know: in >= 0 */
1698 return SCM_MAKINUM (in
);
1701 return SCM_MAKINUM (-1L & ((1L << bits
) - 1));
1703 return SCM_MAKINUM (0);
1705 else if (SCM_BIGP (n
))
1709 SCM num1
= SCM_MAKINUM (1L);
1710 SCM num2
= SCM_MAKINUM (2L);
1711 SCM bits
= SCM_MAKINUM (iend
- istart
);
1712 SCM mask
= scm_difference (scm_integer_expt (num2
, bits
), num1
);
1713 return scm_logand (mask
, scm_ash (n
, SCM_MAKINUM (-istart
)));
1717 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1721 static const char scm_logtab
[] = {
1722 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1725 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1727 "Return the number of bits in integer @var{n}. If integer is\n"
1728 "positive, the 1-bits in its binary representation are counted.\n"
1729 "If negative, the 0-bits in its two's-complement binary\n"
1730 "representation are counted. If 0, 0 is returned.\n"
1733 "(logcount #b10101010)\n"
1740 #define FUNC_NAME s_scm_logcount
1744 unsigned long int c
= 0;
1745 long int nn
= SCM_INUM (n
);
1750 c
+= scm_logtab
[15 & nn
];
1753 return SCM_MAKINUM (c
);
1755 else if (SCM_BIGP (n
))
1757 unsigned long count
;
1758 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1759 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1761 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1762 scm_remember_upto_here_1 (n
);
1763 return SCM_MAKINUM (count
);
1766 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1771 static const char scm_ilentab
[] = {
1772 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1776 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1778 "Return the number of bits necessary to represent @var{n}.\n"
1781 "(integer-length #b10101010)\n"
1783 "(integer-length 0)\n"
1785 "(integer-length #b1111)\n"
1788 #define FUNC_NAME s_scm_integer_length
1792 unsigned long int c
= 0;
1794 long int nn
= SCM_INUM (n
);
1800 l
= scm_ilentab
[15 & nn
];
1803 return SCM_MAKINUM (c
- 4 + l
);
1805 else if (SCM_BIGP (n
))
1807 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1808 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1809 1 too big, so check for that and adjust. */
1810 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1811 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1812 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1813 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1815 scm_remember_upto_here_1 (n
);
1816 return SCM_MAKINUM (size
);
1819 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1823 /*** NUMBERS -> STRINGS ***/
1825 static const double fx
[] =
1826 { 0.0, 5e-1, 5e-2, 5e-3, 5e-4, 5e-5,
1827 5e-6, 5e-7, 5e-8, 5e-9, 5e-10,
1828 5e-11, 5e-12, 5e-13, 5e-14, 5e-15,
1829 5e-16, 5e-17, 5e-18, 5e-19, 5e-20};
1832 idbl2str (double f
, char *a
)
1834 int efmt
, dpt
, d
, i
, wp
= scm_dblprec
;
1840 #ifdef HAVE_COPYSIGN
1841 double sgn
= copysign (1.0, f
);
1847 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
1853 strcpy (a
, "-inf.0");
1855 strcpy (a
, "+inf.0");
1858 else if (xisnan (f
))
1860 strcpy (a
, "+nan.0");
1870 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
1871 make-uniform-vector, from causing infinite loops. */
1875 if (exp
-- < DBL_MIN_10_EXP
)
1886 if (exp
++ > DBL_MAX_10_EXP
)
1906 if (f
+ fx
[wp
] >= 10.0)
1913 dpt
= (exp
+ 9999) % 3;
1917 efmt
= (exp
< -3) || (exp
> wp
+ 2);
1942 if (f
+ fx
[wp
] >= 1.0)
1956 if ((dpt
> 4) && (exp
> 6))
1958 d
= (a
[0] == '-' ? 2 : 1);
1959 for (i
= ch
++; i
> d
; i
--)
1972 if (a
[ch
- 1] == '.')
1973 a
[ch
++] = '0'; /* trailing zero */
1982 for (i
= 10; i
<= exp
; i
*= 10);
1983 for (i
/= 10; i
; i
/= 10)
1985 a
[ch
++] = exp
/ i
+ '0';
1994 iflo2str (SCM flt
, char *str
)
1997 if (SCM_REALP (flt
))
1998 i
= idbl2str (SCM_REAL_VALUE (flt
), str
);
2001 i
= idbl2str (SCM_COMPLEX_REAL (flt
), str
);
2002 if (SCM_COMPLEX_IMAG (flt
) != 0.0)
2004 double imag
= SCM_COMPLEX_IMAG (flt
);
2005 /* Don't output a '+' for negative numbers or for Inf and
2006 NaN. They will provide their own sign. */
2007 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2009 i
+= idbl2str (imag
, &str
[i
]);
2016 /* convert a long to a string (unterminated). returns the number of
2017 characters in the result.
2019 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2021 scm_iint2str (long num
, int rad
, char *p
)
2025 unsigned long n
= (num
< 0) ? -num
: num
;
2027 for (n
/= rad
; n
> 0; n
/= rad
)
2044 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2049 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2051 "Return a string holding the external representation of the\n"
2052 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2053 "inexact, a radix of 10 will be used.")
2054 #define FUNC_NAME s_scm_number_to_string
2058 if (SCM_UNBNDP (radix
))
2062 SCM_VALIDATE_INUM (2, radix
);
2063 base
= SCM_INUM (radix
);
2064 /* FIXME: ask if range limit was OK, and if so, document */
2065 SCM_ASSERT_RANGE (2, radix
, (base
>= 2) && (base
<= 36));
2070 char num_buf
[SCM_INTBUFLEN
];
2071 size_t length
= scm_iint2str (SCM_INUM (n
), base
, num_buf
);
2072 return scm_mem2string (num_buf
, length
);
2074 else if (SCM_BIGP (n
))
2076 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2077 scm_remember_upto_here_1 (n
);
2078 return scm_take0str (str
);
2080 else if (SCM_FRACTIONP (n
))
2082 scm_i_fraction_reduce (n
);
2083 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2084 scm_mem2string ("/", 1),
2085 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2087 else if (SCM_INEXACTP (n
))
2089 char num_buf
[FLOBUFLEN
];
2090 return scm_mem2string (num_buf
, iflo2str (n
, num_buf
));
2093 SCM_WRONG_TYPE_ARG (1, n
);
2098 /* These print routines used to be stubbed here so that scm_repl.c
2099 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2102 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2104 char num_buf
[FLOBUFLEN
];
2105 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
), port
);
2110 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2113 char num_buf
[FLOBUFLEN
];
2114 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
), port
);
2119 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2122 scm_i_fraction_reduce (sexp
);
2123 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2124 scm_lfwrite (SCM_STRING_CHARS (str
), SCM_STRING_LENGTH (str
), port
);
2125 scm_remember_upto_here_1 (str
);
2130 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2132 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2133 scm_remember_upto_here_1 (exp
);
2134 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2138 /*** END nums->strs ***/
2141 /*** STRINGS -> NUMBERS ***/
2143 /* The following functions implement the conversion from strings to numbers.
2144 * The implementation somehow follows the grammar for numbers as it is given
2145 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2146 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2147 * points should be noted about the implementation:
2148 * * Each function keeps a local index variable 'idx' that points at the
2149 * current position within the parsed string. The global index is only
2150 * updated if the function could parse the corresponding syntactic unit
2152 * * Similarly, the functions keep track of indicators of inexactness ('#',
2153 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2154 * global exactness information is only updated after each part has been
2155 * successfully parsed.
2156 * * Sequences of digits are parsed into temporary variables holding fixnums.
2157 * Only if these fixnums would overflow, the result variables are updated
2158 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2159 * the temporary variables holding the fixnums are cleared, and the process
2160 * starts over again. If for example fixnums were able to store five decimal
2161 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2162 * and the result was computed as 12345 * 100000 + 67890. In other words,
2163 * only every five digits two bignum operations were performed.
2166 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2168 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2170 /* In non ASCII-style encodings the following macro might not work. */
2171 #define XDIGIT2UINT(d) (isdigit (d) ? (d) - '0' : tolower (d) - 'a' + 10)
2174 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2175 unsigned int radix
, enum t_exactness
*p_exactness
)
2177 unsigned int idx
= *p_idx
;
2178 unsigned int hash_seen
= 0;
2179 scm_t_bits shift
= 1;
2181 unsigned int digit_value
;
2191 digit_value
= XDIGIT2UINT (c
);
2192 if (digit_value
>= radix
)
2196 result
= SCM_MAKINUM (digit_value
);
2204 digit_value
= XDIGIT2UINT (c
);
2205 if (digit_value
>= radix
)
2217 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2219 result
= scm_product (result
, SCM_MAKINUM (shift
));
2221 result
= scm_sum (result
, SCM_MAKINUM (add
));
2228 shift
= shift
* radix
;
2229 add
= add
* radix
+ digit_value
;
2234 result
= scm_product (result
, SCM_MAKINUM (shift
));
2236 result
= scm_sum (result
, SCM_MAKINUM (add
));
2240 *p_exactness
= INEXACT
;
2246 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2247 * covers the parts of the rules that start at a potential point. The value
2248 * of the digits up to the point have been parsed by the caller and are given
2249 * in variable result. The content of *p_exactness indicates, whether a hash
2250 * has already been seen in the digits before the point.
2253 /* In non ASCII-style encodings the following macro might not work. */
2254 #define DIGIT2UINT(d) ((d) - '0')
2257 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2258 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2260 unsigned int idx
= *p_idx
;
2261 enum t_exactness x
= *p_exactness
;
2266 if (mem
[idx
] == '.')
2268 scm_t_bits shift
= 1;
2270 unsigned int digit_value
;
2271 SCM big_shift
= SCM_MAKINUM (1);
2282 digit_value
= DIGIT2UINT (c
);
2293 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2295 big_shift
= scm_product (big_shift
, SCM_MAKINUM (shift
));
2296 result
= scm_product (result
, SCM_MAKINUM (shift
));
2298 result
= scm_sum (result
, SCM_MAKINUM (add
));
2306 add
= add
* 10 + digit_value
;
2312 big_shift
= scm_product (big_shift
, SCM_MAKINUM (shift
));
2313 result
= scm_product (result
, SCM_MAKINUM (shift
));
2314 result
= scm_sum (result
, SCM_MAKINUM (add
));
2317 result
= scm_divide (result
, big_shift
);
2319 /* We've seen a decimal point, thus the value is implicitly inexact. */
2331 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2362 exponent
= DIGIT2UINT (c
);
2369 if (exponent
<= SCM_MAXEXP
)
2370 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2376 if (exponent
> SCM_MAXEXP
)
2378 size_t exp_len
= idx
- start
;
2379 SCM exp_string
= scm_mem2string (&mem
[start
], exp_len
);
2380 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2381 scm_out_of_range ("string->number", exp_num
);
2384 e
= scm_integer_expt (SCM_MAKINUM (10), SCM_MAKINUM (exponent
));
2386 result
= scm_product (result
, e
);
2388 result
= scm_divide2real (result
, e
);
2390 /* We've seen an exponent, thus the value is implicitly inexact. */
2408 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2411 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2412 unsigned int radix
, enum t_exactness
*p_exactness
)
2414 unsigned int idx
= *p_idx
;
2420 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2426 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2428 enum t_exactness x
= EXACT
;
2430 /* Cobble up the fractional part. We might want to set the
2431 NaN's mantissa from it. */
2433 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2438 if (mem
[idx
] == '.')
2442 else if (idx
+ 1 == len
)
2444 else if (!isdigit (mem
[idx
+ 1]))
2447 result
= mem2decimal_from_point (SCM_MAKINUM (0), mem
, len
,
2448 p_idx
, p_exactness
);
2452 enum t_exactness x
= EXACT
;
2455 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2456 if (SCM_FALSEP (uinteger
))
2461 else if (mem
[idx
] == '/')
2467 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2468 if (SCM_FALSEP (divisor
))
2471 /* both are int/big here, I assume */
2472 result
= scm_make_ratio (uinteger
, divisor
);
2474 else if (radix
== 10)
2476 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2477 if (SCM_FALSEP (result
))
2488 /* When returning an inexact zero, make sure it is represented as a
2489 floating point value so that we can change its sign.
2491 if (SCM_EQ_P (result
, SCM_MAKINUM(0)) && *p_exactness
== INEXACT
)
2492 result
= scm_make_real (0.0);
2498 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2501 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2502 unsigned int radix
, enum t_exactness
*p_exactness
)
2526 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2527 if (SCM_FALSEP (ureal
))
2529 /* input must be either +i or -i */
2534 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2540 return scm_make_rectangular (SCM_MAKINUM (0), SCM_MAKINUM (sign
));
2547 if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2548 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2557 /* either +<ureal>i or -<ureal>i */
2564 return scm_make_rectangular (SCM_MAKINUM (0), ureal
);
2567 /* polar input: <real>@<real>. */
2592 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2593 if (SCM_FALSEP (angle
))
2598 if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2599 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2601 result
= scm_make_polar (ureal
, angle
);
2606 /* expecting input matching <real>[+-]<ureal>?i */
2613 int sign
= (c
== '+') ? 1 : -1;
2614 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2616 if (SCM_FALSEP (imag
))
2617 imag
= SCM_MAKINUM (sign
);
2618 else if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2619 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2623 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2630 return scm_make_rectangular (ureal
, imag
);
2639 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2641 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2644 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2646 unsigned int idx
= 0;
2647 unsigned int radix
= NO_RADIX
;
2648 enum t_exactness forced_x
= NO_EXACTNESS
;
2649 enum t_exactness implicit_x
= EXACT
;
2652 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2653 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2655 switch (mem
[idx
+ 1])
2658 if (radix
!= NO_RADIX
)
2663 if (radix
!= NO_RADIX
)
2668 if (forced_x
!= NO_EXACTNESS
)
2673 if (forced_x
!= NO_EXACTNESS
)
2678 if (radix
!= NO_RADIX
)
2683 if (radix
!= NO_RADIX
)
2693 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2694 if (radix
== NO_RADIX
)
2695 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2697 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2699 if (SCM_FALSEP (result
))
2705 if (SCM_INEXACTP (result
))
2706 return scm_inexact_to_exact (result
);
2710 if (SCM_INEXACTP (result
))
2713 return scm_exact_to_inexact (result
);
2716 if (implicit_x
== INEXACT
)
2718 if (SCM_INEXACTP (result
))
2721 return scm_exact_to_inexact (result
);
2729 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2730 (SCM string
, SCM radix
),
2731 "Return a number of the maximally precise representation\n"
2732 "expressed by the given @var{string}. @var{radix} must be an\n"
2733 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2734 "is a default radix that may be overridden by an explicit radix\n"
2735 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2736 "supplied, then the default radix is 10. If string is not a\n"
2737 "syntactically valid notation for a number, then\n"
2738 "@code{string->number} returns @code{#f}.")
2739 #define FUNC_NAME s_scm_string_to_number
2743 SCM_VALIDATE_STRING (1, string
);
2744 SCM_VALIDATE_INUM_MIN_DEF_COPY (2, radix
,2,10, base
);
2745 answer
= scm_i_mem2number (SCM_STRING_CHARS (string
),
2746 SCM_STRING_LENGTH (string
),
2748 return scm_return_first (answer
, string
);
2753 /*** END strs->nums ***/
2757 scm_make_real (double x
)
2759 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
2761 SCM_REAL_VALUE (z
) = x
;
2767 scm_make_complex (double x
, double y
)
2770 return scm_make_real (x
);
2774 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
2776 SCM_COMPLEX_REAL (z
) = x
;
2777 SCM_COMPLEX_IMAG (z
) = y
;
2784 scm_bigequal (SCM x
, SCM y
)
2786 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2787 scm_remember_upto_here_2 (x
, y
);
2788 return SCM_BOOL (0 == result
);
2792 scm_real_equalp (SCM x
, SCM y
)
2794 return SCM_BOOL (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
2798 scm_complex_equalp (SCM x
, SCM y
)
2800 return SCM_BOOL (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
2801 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
2805 scm_i_fraction_equalp (SCM x
, SCM y
)
2807 scm_i_fraction_reduce (x
);
2808 scm_i_fraction_reduce (y
);
2809 if (SCM_FALSEP (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
2810 SCM_FRACTION_NUMERATOR (y
)))
2811 || SCM_FALSEP (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
2812 SCM_FRACTION_DENOMINATOR (y
))))
2819 SCM_REGISTER_PROC (s_number_p
, "number?", 1, 0, 0, scm_number_p
);
2820 /* "Return @code{#t} if @var{x} is a number, @code{#f}\n"
2821 * "else. Note that the sets of complex, real, rational and\n"
2822 * "integer values form subsets of the set of numbers, i. e. the\n"
2823 * "predicate will be fulfilled for any number."
2825 SCM_DEFINE (scm_number_p
, "complex?", 1, 0, 0,
2827 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
2828 "otherwise. Note that the sets of real, rational and integer\n"
2829 "values form subsets of the set of complex numbers, i. e. the\n"
2830 "predicate will also be fulfilled if @var{x} is a real,\n"
2831 "rational or integer number.")
2832 #define FUNC_NAME s_scm_number_p
2834 return SCM_BOOL (SCM_NUMBERP (x
));
2839 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
2841 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
2842 "otherwise. Note that the set of integer values forms a subset of\n"
2843 "the set of real numbers, i. e. the predicate will also be\n"
2844 "fulfilled if @var{x} is an integer number.")
2845 #define FUNC_NAME s_scm_real_p
2847 /* we can't represent irrational numbers. */
2848 return scm_rational_p (x
);
2852 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
2854 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
2855 "otherwise. Note that the set of integer values forms a subset of\n"
2856 "the set of rational numbers, i. e. the predicate will also be\n"
2857 "fulfilled if @var{x} is an integer number.")
2858 #define FUNC_NAME s_scm_rational_p
2862 else if (SCM_IMP (x
))
2864 else if (SCM_BIGP (x
))
2866 else if (SCM_FRACTIONP (x
))
2868 else if (SCM_REALP (x
))
2869 /* due to their limited precision, all floating point numbers are
2870 rational as well. */
2878 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
2880 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
2882 #define FUNC_NAME s_scm_integer_p
2891 if (!SCM_INEXACTP (x
))
2893 if (SCM_COMPLEXP (x
))
2895 r
= SCM_REAL_VALUE (x
);
2903 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
2905 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
2907 #define FUNC_NAME s_scm_inexact_p
2909 if (SCM_INEXACTP (x
))
2911 if (SCM_NUMBERP (x
))
2913 SCM_WRONG_TYPE_ARG (1, x
);
2918 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
2919 /* "Return @code{#t} if all parameters are numerically equal." */
2921 scm_num_eq_p (SCM x
, SCM y
)
2925 long xx
= SCM_INUM (x
);
2928 long yy
= SCM_INUM (y
);
2929 return SCM_BOOL (xx
== yy
);
2931 else if (SCM_BIGP (y
))
2933 else if (SCM_REALP (y
))
2934 return SCM_BOOL ((double) xx
== SCM_REAL_VALUE (y
));
2935 else if (SCM_COMPLEXP (y
))
2936 return SCM_BOOL (((double) xx
== SCM_COMPLEX_REAL (y
))
2937 && (0.0 == SCM_COMPLEX_IMAG (y
)));
2938 else if (SCM_FRACTIONP (y
))
2941 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
2943 else if (SCM_BIGP (x
))
2947 else if (SCM_BIGP (y
))
2949 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2950 scm_remember_upto_here_2 (x
, y
);
2951 return SCM_BOOL (0 == cmp
);
2953 else if (SCM_REALP (y
))
2956 if (xisnan (SCM_REAL_VALUE (y
)))
2958 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
2959 scm_remember_upto_here_1 (x
);
2960 return SCM_BOOL (0 == cmp
);
2962 else if (SCM_COMPLEXP (y
))
2965 if (0.0 != SCM_COMPLEX_IMAG (y
))
2967 if (xisnan (SCM_COMPLEX_REAL (y
)))
2969 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
2970 scm_remember_upto_here_1 (x
);
2971 return SCM_BOOL (0 == cmp
);
2973 else if (SCM_FRACTIONP (y
))
2976 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
2978 else if (SCM_REALP (x
))
2981 return SCM_BOOL (SCM_REAL_VALUE (x
) == (double) SCM_INUM (y
));
2982 else if (SCM_BIGP (y
))
2985 if (xisnan (SCM_REAL_VALUE (x
)))
2987 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
2988 scm_remember_upto_here_1 (y
);
2989 return SCM_BOOL (0 == cmp
);
2991 else if (SCM_REALP (y
))
2992 return SCM_BOOL (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
2993 else if (SCM_COMPLEXP (y
))
2994 return SCM_BOOL ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
2995 && (0.0 == SCM_COMPLEX_IMAG (y
)));
2996 else if (SCM_FRACTIONP (y
))
2997 return SCM_BOOL (SCM_REAL_VALUE (x
) == scm_i_fraction2double (y
));
2999 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3001 else if (SCM_COMPLEXP (x
))
3004 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == (double) SCM_INUM (y
))
3005 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3006 else if (SCM_BIGP (y
))
3009 if (0.0 != SCM_COMPLEX_IMAG (x
))
3011 if (xisnan (SCM_COMPLEX_REAL (x
)))
3013 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3014 scm_remember_upto_here_1 (y
);
3015 return SCM_BOOL (0 == cmp
);
3017 else if (SCM_REALP (y
))
3018 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3019 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3020 else if (SCM_COMPLEXP (y
))
3021 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3022 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3023 else if (SCM_FRACTIONP (y
))
3024 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == scm_i_fraction2double (y
))
3025 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3027 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3029 else if (SCM_FRACTIONP (x
))
3033 else if (SCM_BIGP (y
))
3035 else if (SCM_REALP (y
))
3036 return SCM_BOOL (scm_i_fraction2double (x
) == SCM_REAL_VALUE (y
));
3037 else if (SCM_COMPLEXP (y
))
3038 return SCM_BOOL ((scm_i_fraction2double (x
) == SCM_COMPLEX_REAL (y
))
3039 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3040 else if (SCM_FRACTIONP (y
))
3041 return scm_i_fraction_equalp (x
, y
);
3043 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3046 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3050 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3051 /* "Return @code{#t} if the list of parameters is monotonically\n"
3055 scm_less_p (SCM x
, SCM y
)
3059 long xx
= SCM_INUM (x
);
3062 long yy
= SCM_INUM (y
);
3063 return SCM_BOOL (xx
< yy
);
3065 else if (SCM_BIGP (y
))
3067 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3068 scm_remember_upto_here_1 (y
);
3069 return SCM_BOOL (sgn
> 0);
3071 else if (SCM_REALP (y
))
3072 return SCM_BOOL ((double) xx
< SCM_REAL_VALUE (y
));
3073 else if (SCM_FRACTIONP (y
))
3074 return SCM_BOOL ((double) xx
< scm_i_fraction2double (y
));
3076 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3078 else if (SCM_BIGP (x
))
3082 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3083 scm_remember_upto_here_1 (x
);
3084 return SCM_BOOL (sgn
< 0);
3086 else if (SCM_BIGP (y
))
3088 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3089 scm_remember_upto_here_2 (x
, y
);
3090 return SCM_BOOL (cmp
< 0);
3092 else if (SCM_REALP (y
))
3095 if (xisnan (SCM_REAL_VALUE (y
)))
3097 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3098 scm_remember_upto_here_1 (x
);
3099 return SCM_BOOL (cmp
< 0);
3101 else if (SCM_FRACTIONP (y
))
3104 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), scm_i_fraction2double (y
));
3105 scm_remember_upto_here_1 (x
);
3106 return SCM_BOOL (cmp
< 0);
3109 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3111 else if (SCM_REALP (x
))
3114 return SCM_BOOL (SCM_REAL_VALUE (x
) < (double) SCM_INUM (y
));
3115 else if (SCM_BIGP (y
))
3118 if (xisnan (SCM_REAL_VALUE (x
)))
3120 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3121 scm_remember_upto_here_1 (y
);
3122 return SCM_BOOL (cmp
> 0);
3124 else if (SCM_REALP (y
))
3125 return SCM_BOOL (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3126 else if (SCM_FRACTIONP (y
))
3127 return SCM_BOOL (SCM_REAL_VALUE (x
) < scm_i_fraction2double (y
));
3129 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3131 else if (SCM_FRACTIONP (x
))
3134 return SCM_BOOL (scm_i_fraction2double (x
) < (double) SCM_INUM (y
));
3135 else if (SCM_BIGP (y
))
3138 if (xisnan (SCM_REAL_VALUE (x
)))
3140 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), scm_i_fraction2double (x
));
3141 scm_remember_upto_here_1 (y
);
3142 return SCM_BOOL (cmp
> 0);
3144 else if (SCM_REALP (y
))
3145 return SCM_BOOL (scm_i_fraction2double (x
) < SCM_REAL_VALUE (y
));
3146 else if (SCM_FRACTIONP (y
))
3147 return SCM_BOOL (scm_i_fraction2double (x
) < scm_i_fraction2double (y
));
3149 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3152 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3156 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3157 /* "Return @code{#t} if the list of parameters is monotonically\n"
3160 #define FUNC_NAME s_scm_gr_p
3162 scm_gr_p (SCM x
, SCM y
)
3164 if (!SCM_NUMBERP (x
))
3165 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3166 else if (!SCM_NUMBERP (y
))
3167 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3169 return scm_less_p (y
, x
);
3174 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3175 /* "Return @code{#t} if the list of parameters is monotonically\n"
3178 #define FUNC_NAME s_scm_leq_p
3180 scm_leq_p (SCM x
, SCM y
)
3182 if (!SCM_NUMBERP (x
))
3183 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3184 else if (!SCM_NUMBERP (y
))
3185 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3186 else if (SCM_NFALSEP (scm_nan_p (x
)) || SCM_NFALSEP (scm_nan_p (y
)))
3189 return SCM_BOOL_NOT (scm_less_p (y
, x
));
3194 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3195 /* "Return @code{#t} if the list of parameters is monotonically\n"
3198 #define FUNC_NAME s_scm_geq_p
3200 scm_geq_p (SCM x
, SCM y
)
3202 if (!SCM_NUMBERP (x
))
3203 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3204 else if (!SCM_NUMBERP (y
))
3205 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3206 else if (SCM_NFALSEP (scm_nan_p (x
)) || SCM_NFALSEP (scm_nan_p (y
)))
3209 return SCM_BOOL_NOT (scm_less_p (x
, y
));
3214 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3215 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3222 return SCM_BOOL (SCM_EQ_P (z
, SCM_INUM0
));
3223 else if (SCM_BIGP (z
))
3225 else if (SCM_REALP (z
))
3226 return SCM_BOOL (SCM_REAL_VALUE (z
) == 0.0);
3227 else if (SCM_COMPLEXP (z
))
3228 return SCM_BOOL (SCM_COMPLEX_REAL (z
) == 0.0
3229 && SCM_COMPLEX_IMAG (z
) == 0.0);
3230 else if (SCM_FRACTIONP (z
))
3233 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3237 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3238 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3242 scm_positive_p (SCM x
)
3245 return SCM_BOOL (SCM_INUM (x
) > 0);
3246 else if (SCM_BIGP (x
))
3248 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3249 scm_remember_upto_here_1 (x
);
3250 return SCM_BOOL (sgn
> 0);
3252 else if (SCM_REALP (x
))
3253 return SCM_BOOL(SCM_REAL_VALUE (x
) > 0.0);
3254 else if (SCM_FRACTIONP (x
))
3255 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3257 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3261 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3262 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3266 scm_negative_p (SCM x
)
3269 return SCM_BOOL (SCM_INUM (x
) < 0);
3270 else if (SCM_BIGP (x
))
3272 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3273 scm_remember_upto_here_1 (x
);
3274 return SCM_BOOL (sgn
< 0);
3276 else if (SCM_REALP (x
))
3277 return SCM_BOOL(SCM_REAL_VALUE (x
) < 0.0);
3278 else if (SCM_FRACTIONP (x
))
3279 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3281 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3285 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3286 /* "Return the maximum of all parameter values."
3289 scm_max (SCM x
, SCM y
)
3294 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3295 else if (SCM_NUMBERP (x
))
3298 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3303 long xx
= SCM_INUM (x
);
3306 long yy
= SCM_INUM (y
);
3307 return (xx
< yy
) ? y
: x
;
3309 else if (SCM_BIGP (y
))
3311 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3312 scm_remember_upto_here_1 (y
);
3313 return (sgn
< 0) ? x
: y
;
3315 else if (SCM_REALP (y
))
3318 /* if y==NaN then ">" is false and we return NaN */
3319 return (z
> SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3321 else if (SCM_FRACTIONP (y
))
3324 return (z
> scm_i_fraction2double (y
)) ? x
: y
;
3327 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3329 else if (SCM_BIGP (x
))
3333 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3334 scm_remember_upto_here_1 (x
);
3335 return (sgn
< 0) ? y
: x
;
3337 else if (SCM_BIGP (y
))
3339 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3340 scm_remember_upto_here_2 (x
, y
);
3341 return (cmp
> 0) ? x
: y
;
3343 else if (SCM_REALP (y
))
3345 double yy
= SCM_REAL_VALUE (y
);
3349 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3350 scm_remember_upto_here_1 (x
);
3351 return (cmp
> 0) ? x
: y
;
3353 else if (SCM_FRACTIONP (y
))
3355 double yy
= scm_i_fraction2double (y
);
3357 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3358 scm_remember_upto_here_1 (x
);
3359 return (cmp
> 0) ? x
: y
;
3362 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3364 else if (SCM_REALP (x
))
3368 double z
= SCM_INUM (y
);
3369 /* if x==NaN then "<" is false and we return NaN */
3370 return (SCM_REAL_VALUE (x
) < z
) ? scm_make_real (z
) : x
;
3372 else if (SCM_BIGP (y
))
3374 double xx
= SCM_REAL_VALUE (x
);
3378 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3379 scm_remember_upto_here_1 (y
);
3380 return (cmp
< 0) ? x
: y
;
3382 else if (SCM_REALP (y
))
3384 /* if x==NaN then our explicit check means we return NaN
3385 if y==NaN then ">" is false and we return NaN
3386 calling isnan is unavoidable, since it's the only way to know
3387 which of x or y causes any compares to be false */
3388 double xx
= SCM_REAL_VALUE (x
);
3389 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3391 else if (SCM_FRACTIONP (y
))
3393 double yy
= scm_i_fraction2double (y
);
3394 double xx
= SCM_REAL_VALUE (x
);
3395 return (xx
< yy
) ? scm_make_real (yy
) : x
;
3398 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3400 else if (SCM_FRACTIONP (x
))
3404 double z
= SCM_INUM (y
);
3405 return (scm_i_fraction2double (x
) < z
) ? y
: x
;
3407 else if (SCM_BIGP (y
))
3409 double xx
= scm_i_fraction2double (x
);
3411 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3412 scm_remember_upto_here_1 (y
);
3413 return (cmp
< 0) ? x
: y
;
3415 else if (SCM_REALP (y
))
3417 double xx
= scm_i_fraction2double (x
);
3418 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_make_real (xx
);
3420 else if (SCM_FRACTIONP (y
))
3422 double yy
= scm_i_fraction2double (y
);
3423 double xx
= scm_i_fraction2double (x
);
3424 return (xx
< yy
) ? y
: x
;
3427 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3430 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3434 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3435 /* "Return the minium of all parameter values."
3438 scm_min (SCM x
, SCM y
)
3443 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3444 else if (SCM_NUMBERP (x
))
3447 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3452 long xx
= SCM_INUM (x
);
3455 long yy
= SCM_INUM (y
);
3456 return (xx
< yy
) ? x
: y
;
3458 else if (SCM_BIGP (y
))
3460 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3461 scm_remember_upto_here_1 (y
);
3462 return (sgn
< 0) ? y
: x
;
3464 else if (SCM_REALP (y
))
3467 /* if y==NaN then "<" is false and we return NaN */
3468 return (z
< SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3470 else if (SCM_FRACTIONP (y
))
3473 return (z
< scm_i_fraction2double (y
)) ? x
: y
;
3476 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3478 else if (SCM_BIGP (x
))
3482 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3483 scm_remember_upto_here_1 (x
);
3484 return (sgn
< 0) ? x
: y
;
3486 else if (SCM_BIGP (y
))
3488 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3489 scm_remember_upto_here_2 (x
, y
);
3490 return (cmp
> 0) ? y
: x
;
3492 else if (SCM_REALP (y
))
3494 double yy
= SCM_REAL_VALUE (y
);
3498 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3499 scm_remember_upto_here_1 (x
);
3500 return (cmp
> 0) ? y
: x
;
3502 else if (SCM_FRACTIONP (y
))
3504 double yy
= scm_i_fraction2double (y
);
3506 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3507 scm_remember_upto_here_1 (x
);
3508 return (cmp
> 0) ? y
: x
;
3511 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3513 else if (SCM_REALP (x
))
3517 double z
= SCM_INUM (y
);
3518 /* if x==NaN then "<" is false and we return NaN */
3519 return (z
< SCM_REAL_VALUE (x
)) ? scm_make_real (z
) : x
;
3521 else if (SCM_BIGP (y
))
3523 double xx
= SCM_REAL_VALUE (x
);
3527 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3528 scm_remember_upto_here_1 (y
);
3529 return (cmp
< 0) ? y
: x
;
3531 else if (SCM_REALP (y
))
3533 /* if x==NaN then our explicit check means we return NaN
3534 if y==NaN then "<" is false and we return NaN
3535 calling isnan is unavoidable, since it's the only way to know
3536 which of x or y causes any compares to be false */
3537 double xx
= SCM_REAL_VALUE (x
);
3538 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3540 else if (SCM_FRACTIONP (y
))
3542 double yy
= scm_i_fraction2double (y
);
3543 double xx
= SCM_REAL_VALUE (x
);
3544 return (yy
< xx
) ? scm_make_real (yy
) : x
;
3547 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3549 else if (SCM_FRACTIONP (x
))
3553 double z
= SCM_INUM (y
);
3554 return (scm_i_fraction2double (x
) < z
) ? x
: y
;
3556 else if (SCM_BIGP (y
))
3558 double xx
= scm_i_fraction2double (x
);
3560 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3561 scm_remember_upto_here_1 (y
);
3562 return (cmp
< 0) ? y
: x
;
3564 else if (SCM_REALP (y
))
3566 double xx
= scm_i_fraction2double (x
);
3567 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_make_real (xx
);
3569 else if (SCM_FRACTIONP (y
))
3571 double yy
= scm_i_fraction2double (y
);
3572 double xx
= scm_i_fraction2double (x
);
3573 return (xx
< yy
) ? x
: y
;
3576 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3579 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3583 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3584 /* "Return the sum of all parameter values. Return 0 if called without\n"
3588 scm_sum (SCM x
, SCM y
)
3592 if (SCM_NUMBERP (x
)) return x
;
3593 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3594 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3601 long xx
= SCM_INUM (x
);
3602 long yy
= SCM_INUM (y
);
3603 long int z
= xx
+ yy
;
3604 return SCM_FIXABLE (z
) ? SCM_MAKINUM (z
) : scm_i_long2big (z
);
3606 else if (SCM_BIGP (y
))
3611 else if (SCM_REALP (y
))
3613 long int xx
= SCM_INUM (x
);
3614 return scm_make_real (xx
+ SCM_REAL_VALUE (y
));
3616 else if (SCM_COMPLEXP (y
))
3618 long int xx
= SCM_INUM (x
);
3619 return scm_make_complex (xx
+ SCM_COMPLEX_REAL (y
),
3620 SCM_COMPLEX_IMAG (y
));
3622 else if (SCM_FRACTIONP (y
))
3623 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3624 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3625 SCM_FRACTION_DENOMINATOR (y
));
3627 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3628 } else if (SCM_BIGP (x
))
3635 inum
= SCM_INUM (y
);
3638 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3641 SCM result
= scm_i_mkbig ();
3642 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3643 scm_remember_upto_here_1 (x
);
3644 /* we know the result will have to be a bignum */
3647 return scm_i_normbig (result
);
3651 SCM result
= scm_i_mkbig ();
3652 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3653 scm_remember_upto_here_1 (x
);
3654 /* we know the result will have to be a bignum */
3657 return scm_i_normbig (result
);
3660 else if (SCM_BIGP (y
))
3662 SCM result
= scm_i_mkbig ();
3663 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3664 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3665 mpz_add (SCM_I_BIG_MPZ (result
),
3668 scm_remember_upto_here_2 (x
, y
);
3669 /* we know the result will have to be a bignum */
3672 return scm_i_normbig (result
);
3674 else if (SCM_REALP (y
))
3676 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3677 scm_remember_upto_here_1 (x
);
3678 return scm_make_real (result
);
3680 else if (SCM_COMPLEXP (y
))
3682 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3683 + SCM_COMPLEX_REAL (y
));
3684 scm_remember_upto_here_1 (x
);
3685 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
3687 else if (SCM_FRACTIONP (y
))
3688 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3689 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3690 SCM_FRACTION_DENOMINATOR (y
));
3692 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3694 else if (SCM_REALP (x
))
3697 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_INUM (y
));
3698 else if (SCM_BIGP (y
))
3700 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3701 scm_remember_upto_here_1 (y
);
3702 return scm_make_real (result
);
3704 else if (SCM_REALP (y
))
3705 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3706 else if (SCM_COMPLEXP (y
))
3707 return scm_make_complex (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3708 SCM_COMPLEX_IMAG (y
));
3709 else if (SCM_FRACTIONP (y
))
3710 return scm_make_real (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3712 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3714 else if (SCM_COMPLEXP (x
))
3717 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_INUM (y
),
3718 SCM_COMPLEX_IMAG (x
));
3719 else if (SCM_BIGP (y
))
3721 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3722 + SCM_COMPLEX_REAL (x
));
3723 scm_remember_upto_here_1 (y
);
3724 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (x
));
3726 else if (SCM_REALP (y
))
3727 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3728 SCM_COMPLEX_IMAG (x
));
3729 else if (SCM_COMPLEXP (y
))
3730 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3731 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3732 else if (SCM_FRACTIONP (y
))
3733 return scm_make_complex (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3734 SCM_COMPLEX_IMAG (x
));
3736 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3738 else if (SCM_FRACTIONP (x
))
3741 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3742 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3743 SCM_FRACTION_DENOMINATOR (x
));
3744 else if (SCM_BIGP (y
))
3745 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3746 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3747 SCM_FRACTION_DENOMINATOR (x
));
3748 else if (SCM_REALP (y
))
3749 return scm_make_real (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
3750 else if (SCM_COMPLEXP (y
))
3751 return scm_make_complex (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
3752 SCM_COMPLEX_IMAG (y
));
3753 else if (SCM_FRACTIONP (y
))
3754 /* a/b + c/d = (ad + bc) / bd */
3755 return scm_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
3756 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
3757 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
3759 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3762 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
3766 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
3767 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
3768 * the sum of all but the first argument are subtracted from the first
3770 #define FUNC_NAME s_difference
3772 scm_difference (SCM x
, SCM y
)
3777 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
3781 long xx
= -SCM_INUM (x
);
3782 if (SCM_FIXABLE (xx
))
3783 return SCM_MAKINUM (xx
);
3785 return scm_i_long2big (xx
);
3787 else if (SCM_BIGP (x
))
3788 /* FIXME: do we really need to normalize here? */
3789 return scm_i_normbig (scm_i_clonebig (x
, 0));
3790 else if (SCM_REALP (x
))
3791 return scm_make_real (-SCM_REAL_VALUE (x
));
3792 else if (SCM_COMPLEXP (x
))
3793 return scm_make_complex (-SCM_COMPLEX_REAL (x
),
3794 -SCM_COMPLEX_IMAG (x
));
3795 else if (SCM_FRACTIONP (x
))
3796 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
3797 SCM_FRACTION_DENOMINATOR (x
));
3799 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
3806 long int xx
= SCM_INUM (x
);
3807 long int yy
= SCM_INUM (y
);
3808 long int z
= xx
- yy
;
3809 if (SCM_FIXABLE (z
))
3810 return SCM_MAKINUM (z
);
3812 return scm_i_long2big (z
);
3814 else if (SCM_BIGP (y
))
3816 /* inum-x - big-y */
3817 long xx
= SCM_INUM (x
);
3820 return scm_i_clonebig (y
, 0);
3823 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3824 SCM result
= scm_i_mkbig ();
3827 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
3830 /* x - y == -(y + -x) */
3831 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
3832 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
3834 scm_remember_upto_here_1 (y
);
3836 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
3837 /* we know the result will have to be a bignum */
3840 return scm_i_normbig (result
);
3843 else if (SCM_REALP (y
))
3845 long int xx
= SCM_INUM (x
);
3846 return scm_make_real (xx
- SCM_REAL_VALUE (y
));
3848 else if (SCM_COMPLEXP (y
))
3850 long int xx
= SCM_INUM (x
);
3851 return scm_make_complex (xx
- SCM_COMPLEX_REAL (y
),
3852 - SCM_COMPLEX_IMAG (y
));
3854 else if (SCM_FRACTIONP (y
))
3855 /* a - b/c = (ac - b) / c */
3856 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
3857 SCM_FRACTION_NUMERATOR (y
)),
3858 SCM_FRACTION_DENOMINATOR (y
));
3860 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3862 else if (SCM_BIGP (x
))
3866 /* big-x - inum-y */
3867 long yy
= SCM_INUM (y
);
3868 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3870 scm_remember_upto_here_1 (x
);
3872 return SCM_FIXABLE (-yy
) ? SCM_MAKINUM (-yy
) : scm_long2num (-yy
);
3875 SCM result
= scm_i_mkbig ();
3878 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
3880 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
3881 scm_remember_upto_here_1 (x
);
3883 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
3884 /* we know the result will have to be a bignum */
3887 return scm_i_normbig (result
);
3890 else if (SCM_BIGP (y
))
3892 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3893 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3894 SCM result
= scm_i_mkbig ();
3895 mpz_sub (SCM_I_BIG_MPZ (result
),
3898 scm_remember_upto_here_2 (x
, y
);
3899 /* we know the result will have to be a bignum */
3900 if ((sgn_x
== 1) && (sgn_y
== -1))
3902 if ((sgn_x
== -1) && (sgn_y
== 1))
3904 return scm_i_normbig (result
);
3906 else if (SCM_REALP (y
))
3908 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
3909 scm_remember_upto_here_1 (x
);
3910 return scm_make_real (result
);
3912 else if (SCM_COMPLEXP (y
))
3914 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3915 - SCM_COMPLEX_REAL (y
));
3916 scm_remember_upto_here_1 (x
);
3917 return scm_make_complex (real_part
, - SCM_COMPLEX_IMAG (y
));
3919 else if (SCM_FRACTIONP (y
))
3920 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
3921 SCM_FRACTION_NUMERATOR (y
)),
3922 SCM_FRACTION_DENOMINATOR (y
));
3923 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3925 else if (SCM_REALP (x
))
3928 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_INUM (y
));
3929 else if (SCM_BIGP (y
))
3931 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
3932 scm_remember_upto_here_1 (x
);
3933 return scm_make_real (result
);
3935 else if (SCM_REALP (y
))
3936 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
3937 else if (SCM_COMPLEXP (y
))
3938 return scm_make_complex (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
3939 -SCM_COMPLEX_IMAG (y
));
3940 else if (SCM_FRACTIONP (y
))
3941 return scm_make_real (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
3943 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3945 else if (SCM_COMPLEXP (x
))
3948 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_INUM (y
),
3949 SCM_COMPLEX_IMAG (x
));
3950 else if (SCM_BIGP (y
))
3952 double real_part
= (SCM_COMPLEX_REAL (x
)
3953 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
3954 scm_remember_upto_here_1 (x
);
3955 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
3957 else if (SCM_REALP (y
))
3958 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
3959 SCM_COMPLEX_IMAG (x
));
3960 else if (SCM_COMPLEXP (y
))
3961 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
3962 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
3963 else if (SCM_FRACTIONP (y
))
3964 return scm_make_complex (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
3965 SCM_COMPLEX_IMAG (x
));
3967 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3969 else if (SCM_FRACTIONP (x
))
3972 /* a/b - c = (a - cb) / b */
3973 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
3974 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
3975 SCM_FRACTION_DENOMINATOR (x
));
3976 else if (SCM_BIGP (y
))
3977 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
3978 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
3979 SCM_FRACTION_DENOMINATOR (x
));
3980 else if (SCM_REALP (y
))
3981 return scm_make_real (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
3982 else if (SCM_COMPLEXP (y
))
3983 return scm_make_complex (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
3984 -SCM_COMPLEX_IMAG (y
));
3985 else if (SCM_FRACTIONP (y
))
3986 /* a/b - c/d = (ad - bc) / bd */
3987 return scm_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
3988 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
3989 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
3991 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3994 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
3999 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4000 /* "Return the product of all arguments. If called without arguments,\n"
4004 scm_product (SCM x
, SCM y
)
4009 return SCM_MAKINUM (1L);
4010 else if (SCM_NUMBERP (x
))
4013 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4025 case 0: return x
; break;
4026 case 1: return y
; break;
4031 long yy
= SCM_INUM (y
);
4033 SCM k
= SCM_MAKINUM (kk
);
4034 if ((kk
== SCM_INUM (k
)) && (kk
/ xx
== yy
))
4038 SCM result
= scm_i_long2big (xx
);
4039 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4040 return scm_i_normbig (result
);
4043 else if (SCM_BIGP (y
))
4045 SCM result
= scm_i_mkbig ();
4046 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4047 scm_remember_upto_here_1 (y
);
4050 else if (SCM_REALP (y
))
4051 return scm_make_real (xx
* SCM_REAL_VALUE (y
));
4052 else if (SCM_COMPLEXP (y
))
4053 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4054 xx
* SCM_COMPLEX_IMAG (y
));
4055 else if (SCM_FRACTIONP (y
))
4056 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4057 SCM_FRACTION_DENOMINATOR (y
));
4059 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4061 else if (SCM_BIGP (x
))
4068 else if (SCM_BIGP (y
))
4070 SCM result
= scm_i_mkbig ();
4071 mpz_mul (SCM_I_BIG_MPZ (result
),
4074 scm_remember_upto_here_2 (x
, y
);
4077 else if (SCM_REALP (y
))
4079 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4080 scm_remember_upto_here_1 (x
);
4081 return scm_make_real (result
);
4083 else if (SCM_COMPLEXP (y
))
4085 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4086 scm_remember_upto_here_1 (x
);
4087 return scm_make_complex (z
* SCM_COMPLEX_REAL (y
),
4088 z
* SCM_COMPLEX_IMAG (y
));
4090 else if (SCM_FRACTIONP (y
))
4091 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4092 SCM_FRACTION_DENOMINATOR (y
));
4094 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4096 else if (SCM_REALP (x
))
4099 return scm_make_real (SCM_INUM (y
) * SCM_REAL_VALUE (x
));
4100 else if (SCM_BIGP (y
))
4102 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4103 scm_remember_upto_here_1 (y
);
4104 return scm_make_real (result
);
4106 else if (SCM_REALP (y
))
4107 return scm_make_real (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4108 else if (SCM_COMPLEXP (y
))
4109 return scm_make_complex (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4110 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4111 else if (SCM_FRACTIONP (y
))
4112 return scm_make_real (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4114 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4116 else if (SCM_COMPLEXP (x
))
4119 return scm_make_complex (SCM_INUM (y
) * SCM_COMPLEX_REAL (x
),
4120 SCM_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4121 else if (SCM_BIGP (y
))
4123 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4124 scm_remember_upto_here_1 (y
);
4125 return scm_make_complex (z
* SCM_COMPLEX_REAL (x
),
4126 z
* SCM_COMPLEX_IMAG (x
));
4128 else if (SCM_REALP (y
))
4129 return scm_make_complex (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4130 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4131 else if (SCM_COMPLEXP (y
))
4133 return scm_make_complex (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4134 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4135 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4136 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4138 else if (SCM_FRACTIONP (y
))
4140 double yy
= scm_i_fraction2double (y
);
4141 return scm_make_complex (yy
* SCM_COMPLEX_REAL (x
),
4142 yy
* SCM_COMPLEX_IMAG (x
));
4145 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4147 else if (SCM_FRACTIONP (x
))
4150 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4151 SCM_FRACTION_DENOMINATOR (x
));
4152 else if (SCM_BIGP (y
))
4153 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4154 SCM_FRACTION_DENOMINATOR (x
));
4155 else if (SCM_REALP (y
))
4156 return scm_make_real (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4157 else if (SCM_COMPLEXP (y
))
4159 double xx
= scm_i_fraction2double (x
);
4160 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4161 xx
* SCM_COMPLEX_IMAG (y
));
4163 else if (SCM_FRACTIONP (y
))
4164 /* a/b * c/d = ac / bd */
4165 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4166 SCM_FRACTION_NUMERATOR (y
)),
4167 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4168 SCM_FRACTION_DENOMINATOR (y
)));
4170 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4173 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4177 scm_num2dbl (SCM a
, const char *why
)
4178 #define FUNC_NAME why
4181 return (double) SCM_INUM (a
);
4182 else if (SCM_BIGP (a
))
4184 double result
= mpz_get_d (SCM_I_BIG_MPZ (a
));
4185 scm_remember_upto_here_1 (a
);
4188 else if (SCM_REALP (a
))
4189 return (SCM_REAL_VALUE (a
));
4190 else if (SCM_FRACTIONP (a
))
4191 return scm_i_fraction2double (a
);
4193 SCM_WRONG_TYPE_ARG (SCM_ARGn
, a
);
4197 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4198 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4199 #define ALLOW_DIVIDE_BY_ZERO
4200 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4203 /* The code below for complex division is adapted from the GNU
4204 libstdc++, which adapted it from f2c's libF77, and is subject to
4207 /****************************************************************
4208 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4210 Permission to use, copy, modify, and distribute this software
4211 and its documentation for any purpose and without fee is hereby
4212 granted, provided that the above copyright notice appear in all
4213 copies and that both that the copyright notice and this
4214 permission notice and warranty disclaimer appear in supporting
4215 documentation, and that the names of AT&T Bell Laboratories or
4216 Bellcore or any of their entities not be used in advertising or
4217 publicity pertaining to distribution of the software without
4218 specific, written prior permission.
4220 AT&T and Bellcore disclaim all warranties with regard to this
4221 software, including all implied warranties of merchantability
4222 and fitness. In no event shall AT&T or Bellcore be liable for
4223 any special, indirect or consequential damages or any damages
4224 whatsoever resulting from loss of use, data or profits, whether
4225 in an action of contract, negligence or other tortious action,
4226 arising out of or in connection with the use or performance of
4228 ****************************************************************/
4230 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4231 /* Divide the first argument by the product of the remaining
4232 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4234 #define FUNC_NAME s_divide
4236 scm_i_divide (SCM x
, SCM y
, int inexact
)
4243 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4244 else if (SCM_INUMP (x
))
4246 long xx
= SCM_INUM (x
);
4247 if (xx
== 1 || xx
== -1)
4249 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4251 scm_num_overflow (s_divide
);
4256 return scm_make_real (1.0 / (double) xx
);
4257 else return scm_make_ratio (SCM_MAKINUM(1), x
);
4260 else if (SCM_BIGP (x
))
4263 return scm_make_real (1.0 / scm_i_big2dbl (x
));
4264 else return scm_make_ratio (SCM_MAKINUM(1), x
);
4266 else if (SCM_REALP (x
))
4268 double xx
= SCM_REAL_VALUE (x
);
4269 #ifndef ALLOW_DIVIDE_BY_ZERO
4271 scm_num_overflow (s_divide
);
4274 return scm_make_real (1.0 / xx
);
4276 else if (SCM_COMPLEXP (x
))
4278 double r
= SCM_COMPLEX_REAL (x
);
4279 double i
= SCM_COMPLEX_IMAG (x
);
4283 double d
= i
* (1.0 + t
* t
);
4284 return scm_make_complex (t
/ d
, -1.0 / d
);
4289 double d
= r
* (1.0 + t
* t
);
4290 return scm_make_complex (1.0 / d
, -t
/ d
);
4293 else if (SCM_FRACTIONP (x
))
4294 return scm_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4295 SCM_FRACTION_NUMERATOR (x
));
4297 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4302 long xx
= SCM_INUM (x
);
4305 long yy
= SCM_INUM (y
);
4308 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4309 scm_num_overflow (s_divide
);
4311 return scm_make_real ((double) xx
/ (double) yy
);
4314 else if (xx
% yy
!= 0)
4317 return scm_make_real ((double) xx
/ (double) yy
);
4318 else return scm_make_ratio (x
, y
);
4323 if (SCM_FIXABLE (z
))
4324 return SCM_MAKINUM (z
);
4326 return scm_i_long2big (z
);
4329 else if (SCM_BIGP (y
))
4332 return scm_make_real ((double) xx
/ scm_i_big2dbl (y
));
4333 else return scm_make_ratio (x
, y
);
4335 else if (SCM_REALP (y
))
4337 double yy
= SCM_REAL_VALUE (y
);
4338 #ifndef ALLOW_DIVIDE_BY_ZERO
4340 scm_num_overflow (s_divide
);
4343 return scm_make_real ((double) xx
/ yy
);
4345 else if (SCM_COMPLEXP (y
))
4348 complex_div
: /* y _must_ be a complex number */
4350 double r
= SCM_COMPLEX_REAL (y
);
4351 double i
= SCM_COMPLEX_IMAG (y
);
4355 double d
= i
* (1.0 + t
* t
);
4356 return scm_make_complex ((a
* t
) / d
, -a
/ d
);
4361 double d
= r
* (1.0 + t
* t
);
4362 return scm_make_complex (a
/ d
, -(a
* t
) / d
);
4366 else if (SCM_FRACTIONP (y
))
4367 /* a / b/c = ac / b */
4368 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4369 SCM_FRACTION_NUMERATOR (y
));
4371 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4373 else if (SCM_BIGP (x
))
4377 long int yy
= SCM_INUM (y
);
4380 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4381 scm_num_overflow (s_divide
);
4383 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4384 scm_remember_upto_here_1 (x
);
4385 return (sgn
== 0) ? scm_nan () : scm_inf ();
4392 /* FIXME: HMM, what are the relative performance issues here?
4393 We need to test. Is it faster on average to test
4394 divisible_p, then perform whichever operation, or is it
4395 faster to perform the integer div opportunistically and
4396 switch to real if there's a remainder? For now we take the
4397 middle ground: test, then if divisible, use the faster div
4400 long abs_yy
= yy
< 0 ? -yy
: yy
;
4401 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4405 SCM result
= scm_i_mkbig ();
4406 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4407 scm_remember_upto_here_1 (x
);
4409 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4410 return scm_i_normbig (result
);
4415 return scm_make_real (scm_i_big2dbl (x
) / (double) yy
);
4416 else return scm_make_ratio (x
, y
);
4420 else if (SCM_BIGP (y
))
4422 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4425 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4426 scm_num_overflow (s_divide
);
4428 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4429 scm_remember_upto_here_1 (x
);
4430 return (sgn
== 0) ? scm_nan () : scm_inf ();
4436 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4440 SCM result
= scm_i_mkbig ();
4441 mpz_divexact (SCM_I_BIG_MPZ (result
),
4444 scm_remember_upto_here_2 (x
, y
);
4445 return scm_i_normbig (result
);
4451 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4452 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4453 scm_remember_upto_here_2 (x
, y
);
4454 return scm_make_real (dbx
/ dby
);
4456 else return scm_make_ratio (x
, y
);
4460 else if (SCM_REALP (y
))
4462 double yy
= SCM_REAL_VALUE (y
);
4463 #ifndef ALLOW_DIVIDE_BY_ZERO
4465 scm_num_overflow (s_divide
);
4468 return scm_make_real (scm_i_big2dbl (x
) / yy
);
4470 else if (SCM_COMPLEXP (y
))
4472 a
= scm_i_big2dbl (x
);
4475 else if (SCM_FRACTIONP (y
))
4476 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4477 SCM_FRACTION_NUMERATOR (y
));
4479 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4481 else if (SCM_REALP (x
))
4483 double rx
= SCM_REAL_VALUE (x
);
4486 long int yy
= SCM_INUM (y
);
4487 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4489 scm_num_overflow (s_divide
);
4492 return scm_make_real (rx
/ (double) yy
);
4494 else if (SCM_BIGP (y
))
4496 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4497 scm_remember_upto_here_1 (y
);
4498 return scm_make_real (rx
/ dby
);
4500 else if (SCM_REALP (y
))
4502 double yy
= SCM_REAL_VALUE (y
);
4503 #ifndef ALLOW_DIVIDE_BY_ZERO
4505 scm_num_overflow (s_divide
);
4508 return scm_make_real (rx
/ yy
);
4510 else if (SCM_COMPLEXP (y
))
4515 else if (SCM_FRACTIONP (y
))
4516 return scm_make_real (rx
/ scm_i_fraction2double (y
));
4518 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4520 else if (SCM_COMPLEXP (x
))
4522 double rx
= SCM_COMPLEX_REAL (x
);
4523 double ix
= SCM_COMPLEX_IMAG (x
);
4526 long int yy
= SCM_INUM (y
);
4527 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4529 scm_num_overflow (s_divide
);
4534 return scm_make_complex (rx
/ d
, ix
/ d
);
4537 else if (SCM_BIGP (y
))
4539 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4540 scm_remember_upto_here_1 (y
);
4541 return scm_make_complex (rx
/ dby
, ix
/ dby
);
4543 else if (SCM_REALP (y
))
4545 double yy
= SCM_REAL_VALUE (y
);
4546 #ifndef ALLOW_DIVIDE_BY_ZERO
4548 scm_num_overflow (s_divide
);
4551 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4553 else if (SCM_COMPLEXP (y
))
4555 double ry
= SCM_COMPLEX_REAL (y
);
4556 double iy
= SCM_COMPLEX_IMAG (y
);
4560 double d
= iy
* (1.0 + t
* t
);
4561 return scm_make_complex ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4566 double d
= ry
* (1.0 + t
* t
);
4567 return scm_make_complex ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4570 else if (SCM_FRACTIONP (y
))
4572 double yy
= scm_i_fraction2double (y
);
4573 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4576 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4578 else if (SCM_FRACTIONP (x
))
4582 long int yy
= SCM_INUM (y
);
4583 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4585 scm_num_overflow (s_divide
);
4588 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4589 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4591 else if (SCM_BIGP (y
))
4593 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4594 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4596 else if (SCM_REALP (y
))
4598 double yy
= SCM_REAL_VALUE (y
);
4599 #ifndef ALLOW_DIVIDE_BY_ZERO
4601 scm_num_overflow (s_divide
);
4604 return scm_make_real (scm_i_fraction2double (x
) / yy
);
4606 else if (SCM_COMPLEXP (y
))
4608 a
= scm_i_fraction2double (x
);
4611 else if (SCM_FRACTIONP (y
))
4612 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4613 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4615 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4618 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4622 scm_divide (SCM x
, SCM y
)
4624 return scm_i_divide (x
, y
, 0);
4627 static SCM
scm_divide2real (SCM x
, SCM y
)
4629 return scm_i_divide (x
, y
, 1);
4635 scm_asinh (double x
)
4640 #define asinh scm_asinh
4641 return log (x
+ sqrt (x
* x
+ 1));
4644 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4645 /* "Return the inverse hyperbolic sine of @var{x}."
4650 scm_acosh (double x
)
4655 #define acosh scm_acosh
4656 return log (x
+ sqrt (x
* x
- 1));
4659 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4660 /* "Return the inverse hyperbolic cosine of @var{x}."
4665 scm_atanh (double x
)
4670 #define atanh scm_atanh
4671 return 0.5 * log ((1 + x
) / (1 - x
));
4674 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4675 /* "Return the inverse hyperbolic tangent of @var{x}."
4679 /* XXX - eventually, we should remove this definition of scm_round and
4680 rename scm_round_number to scm_round. Likewise for scm_truncate
4681 and scm_truncate_number.
4685 scm_truncate (double x
)
4690 #define trunc scm_truncate
4698 scm_round (double x
)
4700 double plus_half
= x
+ 0.5;
4701 double result
= floor (plus_half
);
4702 /* Adjust so that the scm_round is towards even. */
4703 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4708 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4710 "Round the number @var{x} towards zero.")
4711 #define FUNC_NAME s_scm_truncate_number
4713 if (SCM_FALSEP (scm_negative_p (x
)))
4714 return scm_floor (x
);
4716 return scm_ceiling (x
);
4720 static SCM exactly_one_half
;
4722 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4724 "Round the number @var{x} towards the nearest integer. "
4725 "When it is exactly halfway between two integers, "
4726 "round towards the even one.")
4727 #define FUNC_NAME s_scm_round_number
4729 SCM plus_half
= scm_sum (x
, exactly_one_half
);
4730 SCM result
= scm_floor (plus_half
);
4731 /* Adjust so that the scm_round is towards even. */
4732 if (!SCM_FALSEP (scm_num_eq_p (plus_half
, result
))
4733 && !SCM_FALSEP (scm_odd_p (result
)))
4734 return scm_difference (result
, SCM_MAKINUM (1));
4740 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
4742 "Round the number @var{x} towards minus infinity.")
4743 #define FUNC_NAME s_scm_floor
4745 if (SCM_INUMP (x
) || SCM_BIGP (x
))
4747 else if (SCM_REALP (x
))
4748 return scm_make_real (floor (SCM_REAL_VALUE (x
)));
4749 else if (SCM_FRACTIONP (x
))
4751 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
4752 SCM_FRACTION_DENOMINATOR (x
));
4753 if (SCM_FALSEP (scm_negative_p (x
)))
4755 /* For positive x, rounding towards zero is correct. */
4760 /* For negative x, we need to return q-1 unless x is an
4761 integer. But fractions are never integer, per our
4763 return scm_difference (q
, SCM_MAKINUM (1));
4767 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
4771 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
4773 "Round the number @var{x} towards infinity.")
4774 #define FUNC_NAME s_scm_ceiling
4776 if (SCM_INUMP (x
) || SCM_BIGP (x
))
4778 else if (SCM_REALP (x
))
4779 return scm_make_real (ceil (SCM_REAL_VALUE (x
)));
4780 else if (SCM_FRACTIONP (x
))
4782 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
4783 SCM_FRACTION_DENOMINATOR (x
));
4784 if (SCM_FALSEP (scm_positive_p (x
)))
4786 /* For negative x, rounding towards zero is correct. */
4791 /* For positive x, we need to return q+1 unless x is an
4792 integer. But fractions are never integer, per our
4794 return scm_sum (q
, SCM_MAKINUM (1));
4798 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
4802 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
4803 /* "Return the square root of the real number @var{x}."
4805 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
4806 /* "Return the absolute value of the real number @var{x}."
4808 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
4809 /* "Return the @var{x}th power of e."
4811 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
4812 /* "Return the natural logarithm of the real number @var{x}."
4814 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
4815 /* "Return the sine of the real number @var{x}."
4817 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
4818 /* "Return the cosine of the real number @var{x}."
4820 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
4821 /* "Return the tangent of the real number @var{x}."
4823 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
4824 /* "Return the arc sine of the real number @var{x}."
4826 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
4827 /* "Return the arc cosine of the real number @var{x}."
4829 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
4830 /* "Return the arc tangent of the real number @var{x}."
4832 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
4833 /* "Return the hyperbolic sine of the real number @var{x}."
4835 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
4836 /* "Return the hyperbolic cosine of the real number @var{x}."
4838 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
4839 /* "Return the hyperbolic tangent of the real number @var{x}."
4847 static void scm_two_doubles (SCM x
,
4849 const char *sstring
,
4853 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
4856 xy
->x
= SCM_INUM (x
);
4857 else if (SCM_BIGP (x
))
4858 xy
->x
= scm_i_big2dbl (x
);
4859 else if (SCM_REALP (x
))
4860 xy
->x
= SCM_REAL_VALUE (x
);
4861 else if (SCM_FRACTIONP (x
))
4862 xy
->x
= scm_i_fraction2double (x
);
4864 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
4867 xy
->y
= SCM_INUM (y
);
4868 else if (SCM_BIGP (y
))
4869 xy
->y
= scm_i_big2dbl (y
);
4870 else if (SCM_REALP (y
))
4871 xy
->y
= SCM_REAL_VALUE (y
);
4872 else if (SCM_FRACTIONP (y
))
4873 xy
->y
= scm_i_fraction2double (y
);
4875 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
4879 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
4881 "Return @var{x} raised to the power of @var{y}. This\n"
4882 "procedure does not accept complex arguments.")
4883 #define FUNC_NAME s_scm_sys_expt
4886 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4887 return scm_make_real (pow (xy
.x
, xy
.y
));
4892 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
4894 "Return the arc tangent of the two arguments @var{x} and\n"
4895 "@var{y}. This is similar to calculating the arc tangent of\n"
4896 "@var{x} / @var{y}, except that the signs of both arguments\n"
4897 "are used to determine the quadrant of the result. This\n"
4898 "procedure does not accept complex arguments.")
4899 #define FUNC_NAME s_scm_sys_atan2
4902 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4903 return scm_make_real (atan2 (xy
.x
, xy
.y
));
4908 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
4909 (SCM real
, SCM imaginary
),
4910 "Return a complex number constructed of the given @var{real} and\n"
4911 "@var{imaginary} parts.")
4912 #define FUNC_NAME s_scm_make_rectangular
4915 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
4916 return scm_make_complex (xy
.x
, xy
.y
);
4922 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
4924 "Return the complex number @var{x} * e^(i * @var{y}).")
4925 #define FUNC_NAME s_scm_make_polar
4929 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4931 sincos (xy
.y
, &s
, &c
);
4936 return scm_make_complex (xy
.x
* c
, xy
.x
* s
);
4941 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
4942 /* "Return the real part of the number @var{z}."
4945 scm_real_part (SCM z
)
4949 else if (SCM_BIGP (z
))
4951 else if (SCM_REALP (z
))
4953 else if (SCM_COMPLEXP (z
))
4954 return scm_make_real (SCM_COMPLEX_REAL (z
));
4955 else if (SCM_FRACTIONP (z
))
4956 return scm_make_real (scm_i_fraction2double (z
));
4958 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
4962 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
4963 /* "Return the imaginary part of the number @var{z}."
4966 scm_imag_part (SCM z
)
4970 else if (SCM_BIGP (z
))
4972 else if (SCM_REALP (z
))
4974 else if (SCM_COMPLEXP (z
))
4975 return scm_make_real (SCM_COMPLEX_IMAG (z
));
4976 else if (SCM_FRACTIONP (z
))
4979 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
4982 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
4983 /* "Return the numerator of the number @var{z}."
4986 scm_numerator (SCM z
)
4990 else if (SCM_BIGP (z
))
4992 else if (SCM_FRACTIONP (z
))
4994 scm_i_fraction_reduce (z
);
4995 return SCM_FRACTION_NUMERATOR (z
);
4997 else if (SCM_REALP (z
))
4998 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5000 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5004 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5005 /* "Return the denominator of the number @var{z}."
5008 scm_denominator (SCM z
)
5011 return SCM_MAKINUM (1);
5012 else if (SCM_BIGP (z
))
5013 return SCM_MAKINUM (1);
5014 else if (SCM_FRACTIONP (z
))
5016 scm_i_fraction_reduce (z
);
5017 return SCM_FRACTION_DENOMINATOR (z
);
5019 else if (SCM_REALP (z
))
5020 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5022 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5025 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5026 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5027 * "@code{abs} for real arguments, but also allows complex numbers."
5030 scm_magnitude (SCM z
)
5034 long int zz
= SCM_INUM (z
);
5037 else if (SCM_POSFIXABLE (-zz
))
5038 return SCM_MAKINUM (-zz
);
5040 return scm_i_long2big (-zz
);
5042 else if (SCM_BIGP (z
))
5044 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5045 scm_remember_upto_here_1 (z
);
5047 return scm_i_clonebig (z
, 0);
5051 else if (SCM_REALP (z
))
5052 return scm_make_real (fabs (SCM_REAL_VALUE (z
)));
5053 else if (SCM_COMPLEXP (z
))
5054 return scm_make_real (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5055 else if (SCM_FRACTIONP (z
))
5057 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5059 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5060 SCM_FRACTION_DENOMINATOR (z
));
5063 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5067 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5068 /* "Return the angle of the complex number @var{z}."
5073 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5074 scm_flo0 to save allocating a new flonum with scm_make_real each time.
5075 But if atan2 follows the floating point rounding mode, then the value
5076 is not a constant. Maybe it'd be close enough though. */
5079 if (SCM_INUM (z
) >= 0)
5082 return scm_make_real (atan2 (0.0, -1.0));
5084 else if (SCM_BIGP (z
))
5086 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5087 scm_remember_upto_here_1 (z
);
5089 return scm_make_real (atan2 (0.0, -1.0));
5093 else if (SCM_REALP (z
))
5095 if (SCM_REAL_VALUE (z
) >= 0)
5098 return scm_make_real (atan2 (0.0, -1.0));
5100 else if (SCM_COMPLEXP (z
))
5101 return scm_make_real (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5102 else if (SCM_FRACTIONP (z
))
5104 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5106 else return scm_make_real (atan2 (0.0, -1.0));
5109 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5113 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5114 /* Convert the number @var{x} to its inexact representation.\n"
5117 scm_exact_to_inexact (SCM z
)
5120 return scm_make_real ((double) SCM_INUM (z
));
5121 else if (SCM_BIGP (z
))
5122 return scm_make_real (scm_i_big2dbl (z
));
5123 else if (SCM_FRACTIONP (z
))
5124 return scm_make_real (scm_i_fraction2double (z
));
5125 else if (SCM_INEXACTP (z
))
5128 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5132 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5134 "Return an exact number that is numerically closest to @var{z}.")
5135 #define FUNC_NAME s_scm_inexact_to_exact
5139 else if (SCM_BIGP (z
))
5141 else if (SCM_REALP (z
))
5143 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5144 SCM_OUT_OF_RANGE (1, z
);
5151 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5152 q
= scm_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5153 scm_i_mpz2num (mpq_denref (frac
)));
5155 /* When scm_make_ratio throws, we leak the memory allocated
5162 else if (SCM_FRACTIONP (z
))
5165 SCM_WRONG_TYPE_ARG (1, z
);
5169 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5171 "Return an exact number that is within @var{err} of @var{x}.")
5172 #define FUNC_NAME s_scm_rationalize
5176 else if (SCM_BIGP (x
))
5178 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5180 /* Use continued fractions to find closest ratio. All
5181 arithmetic is done with exact numbers.
5184 SCM ex
= scm_inexact_to_exact (x
);
5185 SCM int_part
= scm_floor (ex
);
5186 SCM tt
= SCM_MAKINUM (1);
5187 SCM a1
= SCM_MAKINUM (0), a2
= SCM_MAKINUM (1), a
= SCM_MAKINUM (0);
5188 SCM b1
= SCM_MAKINUM (1), b2
= SCM_MAKINUM (0), b
= SCM_MAKINUM (0);
5192 if (!SCM_FALSEP (scm_num_eq_p (ex
, int_part
)))
5195 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5196 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5198 /* We stop after a million iterations just to be absolutely sure
5199 that we don't go into an infinite loop. The process normally
5200 converges after less than a dozen iterations.
5203 err
= scm_abs (err
);
5204 while (++i
< 1000000)
5206 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5207 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5208 if (SCM_FALSEP (scm_zero_p (b
)) && /* b != 0 */
5210 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5211 err
))) /* abs(x-a/b) <= err */
5213 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5214 if (SCM_FALSEP (scm_exact_p (x
))
5215 || SCM_FALSEP (scm_exact_p (err
)))
5216 return scm_exact_to_inexact (res
);
5220 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5222 tt
= scm_floor (rx
); /* tt = floor (rx) */
5228 scm_num_overflow (s_scm_rationalize
);
5231 SCM_WRONG_TYPE_ARG (1, x
);
5235 /* if you need to change this, change test-num2integral.c as well */
5236 #if SCM_SIZEOF_LONG_LONG != 0
5238 # define ULLONG_MAX ((unsigned long long) (-1))
5239 # define LLONG_MAX ((long long) (ULLONG_MAX >> 1))
5240 # define LLONG_MIN (~LLONG_MAX)
5244 /* Parameters for creating integer conversion routines.
5246 Define the following preprocessor macros before including
5247 "libguile/num2integral.i.c":
5249 NUM2INTEGRAL - the name of the function for converting from a
5250 Scheme object to the integral type. This function will be
5251 defined when including "num2integral.i.c".
5253 INTEGRAL2NUM - the name of the function for converting from the
5254 integral type to a Scheme object. This function will be defined.
5256 INTEGRAL2BIG - the name of an internal function that createas a
5257 bignum from the integral type. This function will be defined.
5258 The name should start with "scm_i_".
5260 ITYPE - the name of the integral type.
5262 UNSIGNED - Define this to 1 when ITYPE is an unsigned type. Define
5265 UNSIGNED_ITYPE - the name of the the unsigned variant of the
5266 integral type. If you don't define this, it defaults to
5267 "unsigned ITYPE" for signed types and simply "ITYPE" for unsigned
5270 SIZEOF_ITYPE - an expression giving the size of the integral type
5271 in bytes. This expression must be computable by the
5272 preprocessor. (SIZEOF_FOO values are calculated by configure.in
5277 #define NUM2INTEGRAL scm_num2short
5278 #define INTEGRAL2NUM scm_short2num
5279 #define INTEGRAL2BIG scm_i_short2big
5282 #define SIZEOF_ITYPE SIZEOF_SHORT
5283 #include "libguile/num2integral.i.c"
5285 #define NUM2INTEGRAL scm_num2ushort
5286 #define INTEGRAL2NUM scm_ushort2num
5287 #define INTEGRAL2BIG scm_i_ushort2big
5289 #define ITYPE unsigned short
5290 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_SHORT
5291 #include "libguile/num2integral.i.c"
5293 #define NUM2INTEGRAL scm_num2int
5294 #define INTEGRAL2NUM scm_int2num
5295 #define INTEGRAL2BIG scm_i_int2big
5298 #define SIZEOF_ITYPE SIZEOF_INT
5299 #include "libguile/num2integral.i.c"
5301 #define NUM2INTEGRAL scm_num2uint
5302 #define INTEGRAL2NUM scm_uint2num
5303 #define INTEGRAL2BIG scm_i_uint2big
5305 #define ITYPE unsigned int
5306 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_INT
5307 #include "libguile/num2integral.i.c"
5309 #define NUM2INTEGRAL scm_num2long
5310 #define INTEGRAL2NUM scm_long2num
5311 #define INTEGRAL2BIG scm_i_long2big
5314 #define SIZEOF_ITYPE SIZEOF_LONG
5315 #include "libguile/num2integral.i.c"
5317 #define NUM2INTEGRAL scm_num2ulong
5318 #define INTEGRAL2NUM scm_ulong2num
5319 #define INTEGRAL2BIG scm_i_ulong2big
5321 #define ITYPE unsigned long
5322 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG
5323 #include "libguile/num2integral.i.c"
5325 #define NUM2INTEGRAL scm_num2ptrdiff
5326 #define INTEGRAL2NUM scm_ptrdiff2num
5327 #define INTEGRAL2BIG scm_i_ptrdiff2big
5329 #define ITYPE scm_t_ptrdiff
5330 #define UNSIGNED_ITYPE size_t
5331 #define SIZEOF_ITYPE SCM_SIZEOF_SCM_T_PTRDIFF
5332 #include "libguile/num2integral.i.c"
5334 #define NUM2INTEGRAL scm_num2size
5335 #define INTEGRAL2NUM scm_size2num
5336 #define INTEGRAL2BIG scm_i_size2big
5338 #define ITYPE size_t
5339 #define SIZEOF_ITYPE SIZEOF_SIZE_T
5340 #include "libguile/num2integral.i.c"
5342 #if SCM_SIZEOF_LONG_LONG != 0
5344 #ifndef ULONG_LONG_MAX
5345 #define ULONG_LONG_MAX (~0ULL)
5348 #define NUM2INTEGRAL scm_num2long_long
5349 #define INTEGRAL2NUM scm_long_long2num
5350 #define INTEGRAL2BIG scm_i_long_long2big
5352 #define ITYPE long long
5353 #define SIZEOF_ITYPE SIZEOF_LONG_LONG
5354 #include "libguile/num2integral.i.c"
5356 #define NUM2INTEGRAL scm_num2ulong_long
5357 #define INTEGRAL2NUM scm_ulong_long2num
5358 #define INTEGRAL2BIG scm_i_ulong_long2big
5360 #define ITYPE unsigned long long
5361 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG_LONG
5362 #include "libguile/num2integral.i.c"
5364 #endif /* SCM_SIZEOF_LONG_LONG != 0 */
5366 #define NUM2FLOAT scm_num2float
5367 #define FLOAT2NUM scm_float2num
5369 #include "libguile/num2float.i.c"
5371 #define NUM2FLOAT scm_num2double
5372 #define FLOAT2NUM scm_double2num
5373 #define FTYPE double
5374 #include "libguile/num2float.i.c"
5379 #define SIZE_MAX ((size_t) (-1))
5382 #define PTRDIFF_MIN \
5383 ((scm_t_ptrdiff) ((scm_t_ptrdiff) 1 \
5384 << ((sizeof (scm_t_ptrdiff) * SCM_CHAR_BIT) - 1)))
5387 #define PTRDIFF_MAX (~ PTRDIFF_MIN)
5390 #define CHECK(type, v) \
5393 if ((v) != scm_num2##type (scm_##type##2num (v), 1, "check_sanity")) \
5413 CHECK (ptrdiff
, -1);
5415 CHECK (short, SHRT_MAX
);
5416 CHECK (short, SHRT_MIN
);
5417 CHECK (ushort
, USHRT_MAX
);
5418 CHECK (int, INT_MAX
);
5419 CHECK (int, INT_MIN
);
5420 CHECK (uint
, UINT_MAX
);
5421 CHECK (long, LONG_MAX
);
5422 CHECK (long, LONG_MIN
);
5423 CHECK (ulong
, ULONG_MAX
);
5424 CHECK (size
, SIZE_MAX
);
5425 CHECK (ptrdiff
, PTRDIFF_MAX
);
5426 CHECK (ptrdiff
, PTRDIFF_MIN
);
5428 #if SCM_SIZEOF_LONG_LONG != 0
5429 CHECK (long_long
, 0LL);
5430 CHECK (ulong_long
, 0ULL);
5431 CHECK (long_long
, -1LL);
5432 CHECK (long_long
, LLONG_MAX
);
5433 CHECK (long_long
, LLONG_MIN
);
5434 CHECK (ulong_long
, ULLONG_MAX
);
5441 scm_internal_catch (SCM_BOOL_T, check_body, &data, check_handler, &data); \
5442 if (!SCM_FALSEP (data)) abort();
5445 check_body (void *data
)
5447 SCM num
= *(SCM
*) data
;
5448 scm_num2ulong (num
, 1, NULL
);
5450 return SCM_UNSPECIFIED
;
5454 check_handler (void *data
, SCM tag
, SCM throw_args
)
5456 SCM
*num
= (SCM
*) data
;
5459 return SCM_UNSPECIFIED
;
5462 SCM_DEFINE (scm_sys_check_number_conversions
, "%check-number-conversions", 0, 0, 0,
5464 "Number conversion sanity checking.")
5465 #define FUNC_NAME s_scm_sys_check_number_conversions
5467 SCM data
= SCM_MAKINUM (-1);
5469 data
= scm_int2num (INT_MIN
);
5471 data
= scm_ulong2num (ULONG_MAX
);
5472 data
= scm_difference (SCM_INUM0
, data
);
5474 data
= scm_ulong2num (ULONG_MAX
);
5475 data
= scm_sum (SCM_MAKINUM (1), data
); data
= scm_difference (SCM_INUM0
, data
);
5477 data
= scm_int2num (-10000); data
= scm_product (data
, data
); data
= scm_product (data
, data
);
5480 return SCM_UNSPECIFIED
;
5489 abs_most_negative_fixnum
= scm_i_long2big (- SCM_MOST_NEGATIVE_FIXNUM
);
5490 scm_permanent_object (abs_most_negative_fixnum
);
5492 mpz_init_set_si (z_negative_one
, -1);
5494 /* It may be possible to tune the performance of some algorithms by using
5495 * the following constants to avoid the creation of bignums. Please, before
5496 * using these values, remember the two rules of program optimization:
5497 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5498 scm_c_define ("most-positive-fixnum",
5499 SCM_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5500 scm_c_define ("most-negative-fixnum",
5501 SCM_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5503 scm_add_feature ("complex");
5504 scm_add_feature ("inexact");
5505 scm_flo0
= scm_make_real (0.0);
5507 scm_dblprec
= (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5509 { /* determine floating point precision */
5511 double fsum
= 1.0 + f
;
5514 if (++scm_dblprec
> 20)
5522 scm_dblprec
= scm_dblprec
- 1;
5524 #endif /* DBL_DIG */
5530 exactly_one_half
= scm_permanent_object (scm_divide (SCM_MAKINUM (1),
5532 #include "libguile/numbers.x"