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
))
645 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
646 double xx
= SCM_REAL_VALUE (x
);
648 return scm_make_real (-xx
);
652 else if (SCM_FRACTIONP (x
))
654 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
656 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
657 SCM_FRACTION_DENOMINATOR (x
));
660 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
665 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
666 /* "Return the quotient of the numbers @var{x} and @var{y}."
669 scm_quotient (SCM x
, SCM y
)
673 long xx
= SCM_INUM (x
);
676 long yy
= SCM_INUM (y
);
678 scm_num_overflow (s_quotient
);
683 return SCM_MAKINUM (z
);
685 return scm_i_long2big (z
);
688 else if (SCM_BIGP (y
))
690 if ((SCM_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
691 && (scm_i_bigcmp (abs_most_negative_fixnum
, y
) == 0))
692 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
693 return SCM_MAKINUM (-1);
695 return SCM_MAKINUM (0);
698 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
700 else if (SCM_BIGP (x
))
704 long yy
= SCM_INUM (y
);
706 scm_num_overflow (s_quotient
);
711 SCM result
= scm_i_mkbig ();
714 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
717 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
720 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
721 scm_remember_upto_here_1 (x
);
722 return scm_i_normbig (result
);
725 else if (SCM_BIGP (y
))
727 SCM result
= scm_i_mkbig ();
728 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
731 scm_remember_upto_here_2 (x
, y
);
732 return scm_i_normbig (result
);
735 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
738 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
741 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
742 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
744 * "(remainder 13 4) @result{} 1\n"
745 * "(remainder -13 4) @result{} -1\n"
749 scm_remainder (SCM x
, SCM y
)
755 long yy
= SCM_INUM (y
);
757 scm_num_overflow (s_remainder
);
760 long z
= SCM_INUM (x
) % yy
;
761 return SCM_MAKINUM (z
);
764 else if (SCM_BIGP (y
))
766 if ((SCM_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
767 && (scm_i_bigcmp (abs_most_negative_fixnum
, y
) == 0))
768 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
769 return SCM_MAKINUM (0);
774 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
776 else if (SCM_BIGP (x
))
780 long yy
= SCM_INUM (y
);
782 scm_num_overflow (s_remainder
);
785 SCM result
= scm_i_mkbig ();
788 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
789 scm_remember_upto_here_1 (x
);
790 return scm_i_normbig (result
);
793 else if (SCM_BIGP (y
))
795 SCM result
= scm_i_mkbig ();
796 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
799 scm_remember_upto_here_2 (x
, y
);
800 return scm_i_normbig (result
);
803 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
806 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
810 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
811 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
813 * "(modulo 13 4) @result{} 1\n"
814 * "(modulo -13 4) @result{} 3\n"
818 scm_modulo (SCM x
, SCM y
)
822 long xx
= SCM_INUM (x
);
825 long yy
= SCM_INUM (y
);
827 scm_num_overflow (s_modulo
);
830 /* FIXME: I think this may be a bug on some arches -- results
831 of % with negative second arg are undefined... */
849 return SCM_MAKINUM (result
);
852 else if (SCM_BIGP (y
))
854 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
857 scm_num_overflow (s_modulo
);
865 SCM pos_y
= scm_i_clonebig (y
, 0);
866 /* do this after the last scm_op */
867 mpz_init_set_si (z_x
, xx
);
868 result
= pos_y
; /* re-use this bignum */
869 mpz_mod (SCM_I_BIG_MPZ (result
),
871 SCM_I_BIG_MPZ (pos_y
));
872 scm_remember_upto_here_1 (pos_y
);
876 result
= scm_i_mkbig ();
877 /* do this after the last scm_op */
878 mpz_init_set_si (z_x
, xx
);
879 mpz_mod (SCM_I_BIG_MPZ (result
),
882 scm_remember_upto_here_1 (y
);
885 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
886 mpz_add (SCM_I_BIG_MPZ (result
),
888 SCM_I_BIG_MPZ (result
));
889 scm_remember_upto_here_1 (y
);
890 /* and do this before the next one */
892 return scm_i_normbig (result
);
896 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
898 else if (SCM_BIGP (x
))
902 long yy
= SCM_INUM (y
);
904 scm_num_overflow (s_modulo
);
907 SCM result
= scm_i_mkbig ();
908 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
910 (yy
< 0) ? - yy
: yy
);
911 scm_remember_upto_here_1 (x
);
912 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
913 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
914 SCM_I_BIG_MPZ (result
),
916 return scm_i_normbig (result
);
919 else if (SCM_BIGP (y
))
921 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
923 scm_num_overflow (s_modulo
);
926 SCM result
= scm_i_mkbig ();
927 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
928 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
929 mpz_mod (SCM_I_BIG_MPZ (result
),
931 SCM_I_BIG_MPZ (pos_y
));
933 scm_remember_upto_here_1 (x
);
934 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
935 mpz_add (SCM_I_BIG_MPZ (result
),
937 SCM_I_BIG_MPZ (result
));
938 scm_remember_upto_here_2 (y
, pos_y
);
939 return scm_i_normbig (result
);
943 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
946 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
949 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
950 /* "Return the greatest common divisor of all arguments.\n"
951 * "If called without arguments, 0 is returned."
954 scm_gcd (SCM x
, SCM y
)
957 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
963 long xx
= SCM_INUM (x
);
964 long yy
= SCM_INUM (y
);
965 long u
= xx
< 0 ? -xx
: xx
;
966 long v
= yy
< 0 ? -yy
: yy
;
976 /* Determine a common factor 2^k */
977 while (!(1 & (u
| v
)))
983 /* Now, any factor 2^n can be eliminated */
1003 return (SCM_POSFIXABLE (result
)
1004 ? SCM_MAKINUM (result
)
1005 : scm_i_long2big (result
));
1007 else if (SCM_BIGP (y
))
1009 SCM result
= scm_i_mkbig ();
1010 SCM mx
= scm_i_mkbig ();
1011 mpz_set_si (SCM_I_BIG_MPZ (mx
), SCM_INUM (x
));
1012 scm_remember_upto_here_1 (x
);
1013 mpz_gcd (SCM_I_BIG_MPZ (result
),
1016 scm_remember_upto_here_2 (mx
, y
);
1017 return scm_i_normbig (result
);
1020 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1022 else if (SCM_BIGP (x
))
1026 unsigned long result
;
1027 long yy
= SCM_INUM (y
);
1032 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1033 scm_remember_upto_here_1 (x
);
1034 return (SCM_POSFIXABLE (result
)
1035 ? SCM_MAKINUM (result
)
1036 : scm_ulong2num (result
));
1038 else if (SCM_BIGP (y
))
1040 SCM result
= scm_i_mkbig ();
1041 mpz_gcd (SCM_I_BIG_MPZ (result
),
1044 scm_remember_upto_here_2 (x
, y
);
1045 return scm_i_normbig (result
);
1048 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1051 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1054 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1055 /* "Return the least common multiple of the arguments.\n"
1056 * "If called without arguments, 1 is returned."
1059 scm_lcm (SCM n1
, SCM n2
)
1061 if (SCM_UNBNDP (n2
))
1063 if (SCM_UNBNDP (n1
))
1064 return SCM_MAKINUM (1L);
1065 n2
= SCM_MAKINUM (1L);
1068 SCM_GASSERT2 (SCM_INUMP (n1
) || SCM_BIGP (n1
),
1069 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1070 SCM_GASSERT2 (SCM_INUMP (n2
) || SCM_BIGP (n2
),
1071 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1077 SCM d
= scm_gcd (n1
, n2
);
1078 if (SCM_EQ_P (d
, SCM_INUM0
))
1081 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1085 /* inum n1, big n2 */
1088 SCM result
= scm_i_mkbig ();
1089 long nn1
= SCM_INUM (n1
);
1090 if (nn1
== 0) return SCM_INUM0
;
1091 if (nn1
< 0) nn1
= - nn1
;
1092 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1093 scm_remember_upto_here_1 (n2
);
1108 SCM result
= scm_i_mkbig ();
1109 mpz_lcm(SCM_I_BIG_MPZ (result
),
1111 SCM_I_BIG_MPZ (n2
));
1112 scm_remember_upto_here_2(n1
, n2
);
1113 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1119 #ifndef scm_long2num
1120 #define SCM_LOGOP_RETURN(x) scm_ulong2num(x)
1122 #define SCM_LOGOP_RETURN(x) SCM_MAKINUM(x)
1125 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1130 + + + x (map digit:logand X Y)
1131 + - + x (map digit:logand X (lognot (+ -1 Y)))
1132 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1133 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1138 + + + (map digit:logior X Y)
1139 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1140 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1141 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1146 + + + (map digit:logxor X Y)
1147 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1148 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1149 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1154 + + (any digit:logand X Y)
1155 + - (any digit:logand X (lognot (+ -1 Y)))
1156 - + (any digit:logand (lognot (+ -1 X)) Y)
1161 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1163 "Return the bitwise AND of the integer arguments.\n\n"
1165 "(logand) @result{} -1\n"
1166 "(logand 7) @result{} 7\n"
1167 "(logand #b111 #b011 #b001) @result{} 1\n"
1169 #define FUNC_NAME s_scm_logand
1173 if (SCM_UNBNDP (n2
))
1175 if (SCM_UNBNDP (n1
))
1176 return SCM_MAKINUM (-1);
1177 else if (!SCM_NUMBERP (n1
))
1178 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1179 else if (SCM_NUMBERP (n1
))
1182 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1187 nn1
= SCM_INUM (n1
);
1190 long nn2
= SCM_INUM (n2
);
1191 return SCM_MAKINUM (nn1
& nn2
);
1193 else if SCM_BIGP (n2
)
1199 SCM result_z
= scm_i_mkbig ();
1201 mpz_init_set_si (nn1_z
, nn1
);
1202 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1203 scm_remember_upto_here_1 (n2
);
1205 return scm_i_normbig (result_z
);
1209 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1211 else if (SCM_BIGP (n1
))
1216 nn1
= SCM_INUM (n1
);
1219 else if (SCM_BIGP (n2
))
1221 SCM result_z
= scm_i_mkbig ();
1222 mpz_and (SCM_I_BIG_MPZ (result_z
),
1224 SCM_I_BIG_MPZ (n2
));
1225 scm_remember_upto_here_2 (n1
, n2
);
1226 return scm_i_normbig (result_z
);
1229 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1232 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1237 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1239 "Return the bitwise OR of the integer arguments.\n\n"
1241 "(logior) @result{} 0\n"
1242 "(logior 7) @result{} 7\n"
1243 "(logior #b000 #b001 #b011) @result{} 3\n"
1245 #define FUNC_NAME s_scm_logior
1249 if (SCM_UNBNDP (n2
))
1251 if (SCM_UNBNDP (n1
))
1253 else if (SCM_NUMBERP (n1
))
1256 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1261 nn1
= SCM_INUM (n1
);
1264 long nn2
= SCM_INUM (n2
);
1265 return SCM_MAKINUM (nn1
| nn2
);
1267 else if (SCM_BIGP (n2
))
1273 SCM result_z
= scm_i_mkbig ();
1275 mpz_init_set_si (nn1_z
, nn1
);
1276 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1277 scm_remember_upto_here_1 (n2
);
1283 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1285 else if (SCM_BIGP (n1
))
1290 nn1
= SCM_INUM (n1
);
1293 else if (SCM_BIGP (n2
))
1295 SCM result_z
= scm_i_mkbig ();
1296 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1298 SCM_I_BIG_MPZ (n2
));
1299 scm_remember_upto_here_2 (n1
, n2
);
1303 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1306 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1311 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1313 "Return the bitwise XOR of the integer arguments. A bit is\n"
1314 "set in the result if it is set in an odd number of arguments.\n"
1316 "(logxor) @result{} 0\n"
1317 "(logxor 7) @result{} 7\n"
1318 "(logxor #b000 #b001 #b011) @result{} 2\n"
1319 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1321 #define FUNC_NAME s_scm_logxor
1325 if (SCM_UNBNDP (n2
))
1327 if (SCM_UNBNDP (n1
))
1329 else if (SCM_NUMBERP (n1
))
1332 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1337 nn1
= SCM_INUM (n1
);
1340 long nn2
= SCM_INUM (n2
);
1341 return SCM_MAKINUM (nn1
^ nn2
);
1343 else if (SCM_BIGP (n2
))
1347 SCM result_z
= scm_i_mkbig ();
1349 mpz_init_set_si (nn1_z
, nn1
);
1350 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1351 scm_remember_upto_here_1 (n2
);
1353 return scm_i_normbig (result_z
);
1357 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1359 else if (SCM_BIGP (n1
))
1364 nn1
= SCM_INUM (n1
);
1367 else if (SCM_BIGP (n2
))
1369 SCM result_z
= scm_i_mkbig ();
1370 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1372 SCM_I_BIG_MPZ (n2
));
1373 scm_remember_upto_here_2 (n1
, n2
);
1374 return scm_i_normbig (result_z
);
1377 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1380 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1385 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1388 "(logtest j k) @equiv{} (not (zero? (logand j k)))\n\n"
1389 "(logtest #b0100 #b1011) @result{} #f\n"
1390 "(logtest #b0100 #b0111) @result{} #t\n"
1392 #define FUNC_NAME s_scm_logtest
1401 long nk
= SCM_INUM (k
);
1402 return SCM_BOOL (nj
& nk
);
1404 else if (SCM_BIGP (k
))
1412 mpz_init_set_si (nj_z
, nj
);
1413 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1414 scm_remember_upto_here_1 (k
);
1415 result
= SCM_BOOL (mpz_sgn (nj_z
) != 0);
1421 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1423 else if (SCM_BIGP (j
))
1431 else if (SCM_BIGP (k
))
1435 mpz_init (result_z
);
1439 scm_remember_upto_here_2 (j
, k
);
1440 result
= SCM_BOOL (mpz_sgn (result_z
) != 0);
1441 mpz_clear (result_z
);
1445 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1448 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1453 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1456 "(logbit? index j) @equiv{} (logtest (integer-expt 2 index) j)\n\n"
1457 "(logbit? 0 #b1101) @result{} #t\n"
1458 "(logbit? 1 #b1101) @result{} #f\n"
1459 "(logbit? 2 #b1101) @result{} #t\n"
1460 "(logbit? 3 #b1101) @result{} #t\n"
1461 "(logbit? 4 #b1101) @result{} #f\n"
1463 #define FUNC_NAME s_scm_logbit_p
1465 unsigned long int iindex
;
1467 SCM_VALIDATE_INUM_MIN (SCM_ARG1
, index
, 0);
1468 iindex
= (unsigned long int) SCM_INUM (index
);
1471 return SCM_BOOL ((1L << iindex
) & SCM_INUM (j
));
1472 else if (SCM_BIGP (j
))
1474 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1475 scm_remember_upto_here_1 (j
);
1476 return SCM_BOOL (val
);
1479 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1484 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1486 "Return the integer which is the ones-complement of the integer\n"
1490 "(number->string (lognot #b10000000) 2)\n"
1491 " @result{} \"-10000001\"\n"
1492 "(number->string (lognot #b0) 2)\n"
1493 " @result{} \"-1\"\n"
1495 #define FUNC_NAME s_scm_lognot
1497 if (SCM_INUMP (n
)) {
1498 /* No overflow here, just need to toggle all the bits making up the inum.
1499 Enhancement: No need to strip the tag and add it back, could just xor
1500 a block of 1 bits, if that worked with the various debug versions of
1502 return SCM_MAKINUM (~ SCM_INUM (n
));
1504 } else if (SCM_BIGP (n
)) {
1505 SCM result
= scm_i_mkbig ();
1506 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1507 scm_remember_upto_here_1 (n
);
1511 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1516 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1518 "Return @var{n} raised to the non-negative integer exponent\n"
1522 "(integer-expt 2 5)\n"
1524 "(integer-expt -3 3)\n"
1527 #define FUNC_NAME s_scm_integer_expt
1530 SCM z_i2
= SCM_BOOL_F
;
1532 SCM acc
= SCM_MAKINUM (1L);
1534 /* 0^0 == 1 according to R5RS */
1535 if (SCM_EQ_P (n
, SCM_INUM0
) || SCM_EQ_P (n
, acc
))
1536 return SCM_FALSEP (scm_zero_p(k
)) ? n
: acc
;
1537 else if (SCM_EQ_P (n
, SCM_MAKINUM (-1L)))
1538 return SCM_FALSEP (scm_even_p (k
)) ? n
: acc
;
1542 else if (SCM_BIGP (k
))
1544 z_i2
= scm_i_clonebig (k
, 1);
1545 scm_remember_upto_here_1 (k
);
1548 else if (SCM_REALP (k
))
1550 double r
= SCM_REAL_VALUE (k
);
1552 SCM_WRONG_TYPE_ARG (2, k
);
1553 if ((r
> SCM_MOST_POSITIVE_FIXNUM
) || (r
< SCM_MOST_NEGATIVE_FIXNUM
))
1555 z_i2
= scm_i_mkbig ();
1556 mpz_set_d (SCM_I_BIG_MPZ (z_i2
), r
);
1565 SCM_WRONG_TYPE_ARG (2, k
);
1569 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1571 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1572 n
= scm_divide (n
, SCM_UNDEFINED
);
1576 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1580 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1582 return scm_product (acc
, n
);
1584 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1585 acc
= scm_product (acc
, n
);
1586 n
= scm_product (n
, n
);
1587 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1595 n
= scm_divide (n
, SCM_UNDEFINED
);
1602 return scm_product (acc
, n
);
1604 acc
= scm_product (acc
, n
);
1605 n
= scm_product (n
, n
);
1612 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1614 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1615 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1617 "This is effectively a multiplication by 2^@var{cnt}}, and when\n"
1618 "@var{cnt} is negative it's a division, rounded towards negative\n"
1619 "infinity. (Note that this is not the same rounding as\n"
1620 "@code{quotient} does.)\n"
1622 "With @var{n} viewed as an infinite precision twos complement,\n"
1623 "@code{ash} means a left shift introducing zero bits, or a right\n"
1624 "shift dropping bits.\n"
1627 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1628 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1630 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1631 "(ash -23 -2) @result{} -6\n"
1633 #define FUNC_NAME s_scm_ash
1637 SCM_VALIDATE_INUM (2, cnt
);
1639 bits_to_shift
= SCM_INUM (cnt
);
1641 if (bits_to_shift
< 0)
1643 /* Shift right by abs(cnt) bits. This is realized as a division
1644 by div:=2^abs(cnt). However, to guarantee the floor
1645 rounding, negative values require some special treatment.
1647 SCM div
= scm_integer_expt (SCM_MAKINUM (2),
1648 SCM_MAKINUM (-bits_to_shift
));
1650 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1651 if (SCM_FALSEP (scm_negative_p (n
)))
1652 return scm_quotient (n
, div
);
1654 return scm_sum (SCM_MAKINUM (-1L),
1655 scm_quotient (scm_sum (SCM_MAKINUM (1L), n
), div
));
1658 /* Shift left is done by multiplication with 2^CNT */
1659 return scm_product (n
, scm_integer_expt (SCM_MAKINUM (2), cnt
));
1664 #define MIN(x,y) ((x) < (y) ? (x) : (y))
1666 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1667 (SCM n
, SCM start
, SCM end
),
1668 "Return the integer composed of the @var{start} (inclusive)\n"
1669 "through @var{end} (exclusive) bits of @var{n}. The\n"
1670 "@var{start}th bit becomes the 0-th bit in the result.\n"
1673 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1674 " @result{} \"1010\"\n"
1675 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1676 " @result{} \"10110\"\n"
1678 #define FUNC_NAME s_scm_bit_extract
1680 unsigned long int istart
, iend
, bits
;
1681 SCM_VALIDATE_INUM_MIN_COPY (2, start
,0, istart
);
1682 SCM_VALIDATE_INUM_MIN_COPY (3, end
, 0, iend
);
1683 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1685 /* how many bits to keep */
1686 bits
= iend
- istart
;
1690 long int in
= SCM_INUM (n
);
1692 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1693 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in".
1694 FIXME: This shift relies on signed right shifts being arithmetic,
1695 which is not guaranteed by C99. */
1696 in
>>= MIN (istart
, SCM_I_FIXNUM_BIT
-1);
1698 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1700 /* Since we emulate two's complement encoded numbers, this
1701 * special case requires us to produce a result that has
1702 * more bits than can be stored in a fixnum.
1704 SCM result
= scm_i_long2big (in
);
1705 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1710 /* mask down to requisite bits */
1711 bits
= MIN (bits
, SCM_I_FIXNUM_BIT
);
1712 return SCM_MAKINUM (in
& ((1L << bits
) - 1));
1714 else if (SCM_BIGP (n
))
1719 result
= SCM_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1723 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1724 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1725 such bits into a ulong. */
1726 result
= scm_i_mkbig ();
1727 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1728 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1729 result
= scm_i_normbig (result
);
1731 scm_remember_upto_here_1 (n
);
1735 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1740 static const char scm_logtab
[] = {
1741 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1744 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1746 "Return the number of bits in integer @var{n}. If integer is\n"
1747 "positive, the 1-bits in its binary representation are counted.\n"
1748 "If negative, the 0-bits in its two's-complement binary\n"
1749 "representation are counted. If 0, 0 is returned.\n"
1752 "(logcount #b10101010)\n"
1759 #define FUNC_NAME s_scm_logcount
1763 unsigned long int c
= 0;
1764 long int nn
= SCM_INUM (n
);
1769 c
+= scm_logtab
[15 & nn
];
1772 return SCM_MAKINUM (c
);
1774 else if (SCM_BIGP (n
))
1776 unsigned long count
;
1777 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1778 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1780 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1781 scm_remember_upto_here_1 (n
);
1782 return SCM_MAKINUM (count
);
1785 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1790 static const char scm_ilentab
[] = {
1791 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1795 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1797 "Return the number of bits necessary to represent @var{n}.\n"
1800 "(integer-length #b10101010)\n"
1802 "(integer-length 0)\n"
1804 "(integer-length #b1111)\n"
1807 #define FUNC_NAME s_scm_integer_length
1811 unsigned long int c
= 0;
1813 long int nn
= SCM_INUM (n
);
1819 l
= scm_ilentab
[15 & nn
];
1822 return SCM_MAKINUM (c
- 4 + l
);
1824 else if (SCM_BIGP (n
))
1826 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1827 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1828 1 too big, so check for that and adjust. */
1829 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1830 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1831 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1832 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1834 scm_remember_upto_here_1 (n
);
1835 return SCM_MAKINUM (size
);
1838 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1842 /*** NUMBERS -> STRINGS ***/
1844 static const double fx
[] =
1845 { 0.0, 5e-1, 5e-2, 5e-3, 5e-4, 5e-5,
1846 5e-6, 5e-7, 5e-8, 5e-9, 5e-10,
1847 5e-11, 5e-12, 5e-13, 5e-14, 5e-15,
1848 5e-16, 5e-17, 5e-18, 5e-19, 5e-20};
1851 idbl2str (double f
, char *a
)
1853 int efmt
, dpt
, d
, i
, wp
= scm_dblprec
;
1859 #ifdef HAVE_COPYSIGN
1860 double sgn
= copysign (1.0, f
);
1866 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
1872 strcpy (a
, "-inf.0");
1874 strcpy (a
, "+inf.0");
1877 else if (xisnan (f
))
1879 strcpy (a
, "+nan.0");
1889 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
1890 make-uniform-vector, from causing infinite loops. */
1894 if (exp
-- < DBL_MIN_10_EXP
)
1905 if (exp
++ > DBL_MAX_10_EXP
)
1925 if (f
+ fx
[wp
] >= 10.0)
1932 dpt
= (exp
+ 9999) % 3;
1936 efmt
= (exp
< -3) || (exp
> wp
+ 2);
1961 if (f
+ fx
[wp
] >= 1.0)
1975 if ((dpt
> 4) && (exp
> 6))
1977 d
= (a
[0] == '-' ? 2 : 1);
1978 for (i
= ch
++; i
> d
; i
--)
1991 if (a
[ch
- 1] == '.')
1992 a
[ch
++] = '0'; /* trailing zero */
2001 for (i
= 10; i
<= exp
; i
*= 10);
2002 for (i
/= 10; i
; i
/= 10)
2004 a
[ch
++] = exp
/ i
+ '0';
2013 iflo2str (SCM flt
, char *str
)
2016 if (SCM_REALP (flt
))
2017 i
= idbl2str (SCM_REAL_VALUE (flt
), str
);
2020 i
= idbl2str (SCM_COMPLEX_REAL (flt
), str
);
2021 if (SCM_COMPLEX_IMAG (flt
) != 0.0)
2023 double imag
= SCM_COMPLEX_IMAG (flt
);
2024 /* Don't output a '+' for negative numbers or for Inf and
2025 NaN. They will provide their own sign. */
2026 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2028 i
+= idbl2str (imag
, &str
[i
]);
2035 /* convert a long to a string (unterminated). returns the number of
2036 characters in the result.
2038 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2040 scm_iint2str (long num
, int rad
, char *p
)
2044 unsigned long n
= (num
< 0) ? -num
: num
;
2046 for (n
/= rad
; n
> 0; n
/= rad
)
2063 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2068 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2070 "Return a string holding the external representation of the\n"
2071 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2072 "inexact, a radix of 10 will be used.")
2073 #define FUNC_NAME s_scm_number_to_string
2077 if (SCM_UNBNDP (radix
))
2081 SCM_VALIDATE_INUM (2, radix
);
2082 base
= SCM_INUM (radix
);
2083 /* FIXME: ask if range limit was OK, and if so, document */
2084 SCM_ASSERT_RANGE (2, radix
, (base
>= 2) && (base
<= 36));
2089 char num_buf
[SCM_INTBUFLEN
];
2090 size_t length
= scm_iint2str (SCM_INUM (n
), base
, num_buf
);
2091 return scm_mem2string (num_buf
, length
);
2093 else if (SCM_BIGP (n
))
2095 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2096 scm_remember_upto_here_1 (n
);
2097 return scm_take0str (str
);
2099 else if (SCM_FRACTIONP (n
))
2101 scm_i_fraction_reduce (n
);
2102 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2103 scm_mem2string ("/", 1),
2104 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2106 else if (SCM_INEXACTP (n
))
2108 char num_buf
[FLOBUFLEN
];
2109 return scm_mem2string (num_buf
, iflo2str (n
, num_buf
));
2112 SCM_WRONG_TYPE_ARG (1, n
);
2117 /* These print routines used to be stubbed here so that scm_repl.c
2118 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2121 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2123 char num_buf
[FLOBUFLEN
];
2124 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
), port
);
2129 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2132 char num_buf
[FLOBUFLEN
];
2133 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
), port
);
2138 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2141 scm_i_fraction_reduce (sexp
);
2142 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2143 scm_lfwrite (SCM_STRING_CHARS (str
), SCM_STRING_LENGTH (str
), port
);
2144 scm_remember_upto_here_1 (str
);
2149 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2151 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2152 scm_remember_upto_here_1 (exp
);
2153 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2157 /*** END nums->strs ***/
2160 /*** STRINGS -> NUMBERS ***/
2162 /* The following functions implement the conversion from strings to numbers.
2163 * The implementation somehow follows the grammar for numbers as it is given
2164 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2165 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2166 * points should be noted about the implementation:
2167 * * Each function keeps a local index variable 'idx' that points at the
2168 * current position within the parsed string. The global index is only
2169 * updated if the function could parse the corresponding syntactic unit
2171 * * Similarly, the functions keep track of indicators of inexactness ('#',
2172 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2173 * global exactness information is only updated after each part has been
2174 * successfully parsed.
2175 * * Sequences of digits are parsed into temporary variables holding fixnums.
2176 * Only if these fixnums would overflow, the result variables are updated
2177 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2178 * the temporary variables holding the fixnums are cleared, and the process
2179 * starts over again. If for example fixnums were able to store five decimal
2180 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2181 * and the result was computed as 12345 * 100000 + 67890. In other words,
2182 * only every five digits two bignum operations were performed.
2185 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2187 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2189 /* In non ASCII-style encodings the following macro might not work. */
2190 #define XDIGIT2UINT(d) (isdigit (d) ? (d) - '0' : tolower (d) - 'a' + 10)
2193 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2194 unsigned int radix
, enum t_exactness
*p_exactness
)
2196 unsigned int idx
= *p_idx
;
2197 unsigned int hash_seen
= 0;
2198 scm_t_bits shift
= 1;
2200 unsigned int digit_value
;
2210 digit_value
= XDIGIT2UINT (c
);
2211 if (digit_value
>= radix
)
2215 result
= SCM_MAKINUM (digit_value
);
2223 digit_value
= XDIGIT2UINT (c
);
2224 if (digit_value
>= radix
)
2236 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2238 result
= scm_product (result
, SCM_MAKINUM (shift
));
2240 result
= scm_sum (result
, SCM_MAKINUM (add
));
2247 shift
= shift
* radix
;
2248 add
= add
* radix
+ digit_value
;
2253 result
= scm_product (result
, SCM_MAKINUM (shift
));
2255 result
= scm_sum (result
, SCM_MAKINUM (add
));
2259 *p_exactness
= INEXACT
;
2265 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2266 * covers the parts of the rules that start at a potential point. The value
2267 * of the digits up to the point have been parsed by the caller and are given
2268 * in variable result. The content of *p_exactness indicates, whether a hash
2269 * has already been seen in the digits before the point.
2272 /* In non ASCII-style encodings the following macro might not work. */
2273 #define DIGIT2UINT(d) ((d) - '0')
2276 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2277 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2279 unsigned int idx
= *p_idx
;
2280 enum t_exactness x
= *p_exactness
;
2285 if (mem
[idx
] == '.')
2287 scm_t_bits shift
= 1;
2289 unsigned int digit_value
;
2290 SCM big_shift
= SCM_MAKINUM (1);
2301 digit_value
= DIGIT2UINT (c
);
2312 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2314 big_shift
= scm_product (big_shift
, SCM_MAKINUM (shift
));
2315 result
= scm_product (result
, SCM_MAKINUM (shift
));
2317 result
= scm_sum (result
, SCM_MAKINUM (add
));
2325 add
= add
* 10 + digit_value
;
2331 big_shift
= scm_product (big_shift
, SCM_MAKINUM (shift
));
2332 result
= scm_product (result
, SCM_MAKINUM (shift
));
2333 result
= scm_sum (result
, SCM_MAKINUM (add
));
2336 result
= scm_divide (result
, big_shift
);
2338 /* We've seen a decimal point, thus the value is implicitly inexact. */
2350 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2381 exponent
= DIGIT2UINT (c
);
2388 if (exponent
<= SCM_MAXEXP
)
2389 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2395 if (exponent
> SCM_MAXEXP
)
2397 size_t exp_len
= idx
- start
;
2398 SCM exp_string
= scm_mem2string (&mem
[start
], exp_len
);
2399 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2400 scm_out_of_range ("string->number", exp_num
);
2403 e
= scm_integer_expt (SCM_MAKINUM (10), SCM_MAKINUM (exponent
));
2405 result
= scm_product (result
, e
);
2407 result
= scm_divide2real (result
, e
);
2409 /* We've seen an exponent, thus the value is implicitly inexact. */
2427 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2430 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2431 unsigned int radix
, enum t_exactness
*p_exactness
)
2433 unsigned int idx
= *p_idx
;
2439 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2445 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2447 enum t_exactness x
= EXACT
;
2449 /* Cobble up the fractional part. We might want to set the
2450 NaN's mantissa from it. */
2452 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2457 if (mem
[idx
] == '.')
2461 else if (idx
+ 1 == len
)
2463 else if (!isdigit (mem
[idx
+ 1]))
2466 result
= mem2decimal_from_point (SCM_MAKINUM (0), mem
, len
,
2467 p_idx
, p_exactness
);
2471 enum t_exactness x
= EXACT
;
2474 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2475 if (SCM_FALSEP (uinteger
))
2480 else if (mem
[idx
] == '/')
2486 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2487 if (SCM_FALSEP (divisor
))
2490 /* both are int/big here, I assume */
2491 result
= scm_make_ratio (uinteger
, divisor
);
2493 else if (radix
== 10)
2495 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2496 if (SCM_FALSEP (result
))
2507 /* When returning an inexact zero, make sure it is represented as a
2508 floating point value so that we can change its sign.
2510 if (SCM_EQ_P (result
, SCM_MAKINUM(0)) && *p_exactness
== INEXACT
)
2511 result
= scm_make_real (0.0);
2517 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2520 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2521 unsigned int radix
, enum t_exactness
*p_exactness
)
2545 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2546 if (SCM_FALSEP (ureal
))
2548 /* input must be either +i or -i */
2553 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2559 return scm_make_rectangular (SCM_MAKINUM (0), SCM_MAKINUM (sign
));
2566 if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2567 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2576 /* either +<ureal>i or -<ureal>i */
2583 return scm_make_rectangular (SCM_MAKINUM (0), ureal
);
2586 /* polar input: <real>@<real>. */
2611 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2612 if (SCM_FALSEP (angle
))
2617 if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2618 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2620 result
= scm_make_polar (ureal
, angle
);
2625 /* expecting input matching <real>[+-]<ureal>?i */
2632 int sign
= (c
== '+') ? 1 : -1;
2633 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2635 if (SCM_FALSEP (imag
))
2636 imag
= SCM_MAKINUM (sign
);
2637 else if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2638 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2642 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2649 return scm_make_rectangular (ureal
, imag
);
2658 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2660 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2663 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2665 unsigned int idx
= 0;
2666 unsigned int radix
= NO_RADIX
;
2667 enum t_exactness forced_x
= NO_EXACTNESS
;
2668 enum t_exactness implicit_x
= EXACT
;
2671 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2672 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2674 switch (mem
[idx
+ 1])
2677 if (radix
!= NO_RADIX
)
2682 if (radix
!= NO_RADIX
)
2687 if (forced_x
!= NO_EXACTNESS
)
2692 if (forced_x
!= NO_EXACTNESS
)
2697 if (radix
!= NO_RADIX
)
2702 if (radix
!= NO_RADIX
)
2712 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2713 if (radix
== NO_RADIX
)
2714 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2716 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2718 if (SCM_FALSEP (result
))
2724 if (SCM_INEXACTP (result
))
2725 return scm_inexact_to_exact (result
);
2729 if (SCM_INEXACTP (result
))
2732 return scm_exact_to_inexact (result
);
2735 if (implicit_x
== INEXACT
)
2737 if (SCM_INEXACTP (result
))
2740 return scm_exact_to_inexact (result
);
2748 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2749 (SCM string
, SCM radix
),
2750 "Return a number of the maximally precise representation\n"
2751 "expressed by the given @var{string}. @var{radix} must be an\n"
2752 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2753 "is a default radix that may be overridden by an explicit radix\n"
2754 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2755 "supplied, then the default radix is 10. If string is not a\n"
2756 "syntactically valid notation for a number, then\n"
2757 "@code{string->number} returns @code{#f}.")
2758 #define FUNC_NAME s_scm_string_to_number
2762 SCM_VALIDATE_STRING (1, string
);
2763 SCM_VALIDATE_INUM_MIN_DEF_COPY (2, radix
,2,10, base
);
2764 answer
= scm_i_mem2number (SCM_STRING_CHARS (string
),
2765 SCM_STRING_LENGTH (string
),
2767 return scm_return_first (answer
, string
);
2772 /*** END strs->nums ***/
2776 scm_make_real (double x
)
2778 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
2780 SCM_REAL_VALUE (z
) = x
;
2786 scm_make_complex (double x
, double y
)
2789 return scm_make_real (x
);
2793 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
2795 SCM_COMPLEX_REAL (z
) = x
;
2796 SCM_COMPLEX_IMAG (z
) = y
;
2803 scm_bigequal (SCM x
, SCM y
)
2805 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2806 scm_remember_upto_here_2 (x
, y
);
2807 return SCM_BOOL (0 == result
);
2811 scm_real_equalp (SCM x
, SCM y
)
2813 return SCM_BOOL (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
2817 scm_complex_equalp (SCM x
, SCM y
)
2819 return SCM_BOOL (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
2820 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
2824 scm_i_fraction_equalp (SCM x
, SCM y
)
2826 scm_i_fraction_reduce (x
);
2827 scm_i_fraction_reduce (y
);
2828 if (SCM_FALSEP (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
2829 SCM_FRACTION_NUMERATOR (y
)))
2830 || SCM_FALSEP (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
2831 SCM_FRACTION_DENOMINATOR (y
))))
2838 SCM_REGISTER_PROC (s_number_p
, "number?", 1, 0, 0, scm_number_p
);
2839 /* "Return @code{#t} if @var{x} is a number, @code{#f}\n"
2840 * "else. Note that the sets of complex, real, rational and\n"
2841 * "integer values form subsets of the set of numbers, i. e. the\n"
2842 * "predicate will be fulfilled for any number."
2844 SCM_DEFINE (scm_number_p
, "complex?", 1, 0, 0,
2846 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
2847 "otherwise. Note that the sets of real, rational and integer\n"
2848 "values form subsets of the set of complex numbers, i. e. the\n"
2849 "predicate will also be fulfilled if @var{x} is a real,\n"
2850 "rational or integer number.")
2851 #define FUNC_NAME s_scm_number_p
2853 return SCM_BOOL (SCM_NUMBERP (x
));
2858 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
2860 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
2861 "otherwise. Note that the set of integer values forms a subset of\n"
2862 "the set of real numbers, i. e. the predicate will also be\n"
2863 "fulfilled if @var{x} is an integer number.")
2864 #define FUNC_NAME s_scm_real_p
2866 /* we can't represent irrational numbers. */
2867 return scm_rational_p (x
);
2871 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
2873 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
2874 "otherwise. Note that the set of integer values forms a subset of\n"
2875 "the set of rational numbers, i. e. the predicate will also be\n"
2876 "fulfilled if @var{x} is an integer number.")
2877 #define FUNC_NAME s_scm_rational_p
2881 else if (SCM_IMP (x
))
2883 else if (SCM_BIGP (x
))
2885 else if (SCM_FRACTIONP (x
))
2887 else if (SCM_REALP (x
))
2888 /* due to their limited precision, all floating point numbers are
2889 rational as well. */
2897 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
2899 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
2901 #define FUNC_NAME s_scm_integer_p
2910 if (!SCM_INEXACTP (x
))
2912 if (SCM_COMPLEXP (x
))
2914 r
= SCM_REAL_VALUE (x
);
2922 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
2924 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
2926 #define FUNC_NAME s_scm_inexact_p
2928 if (SCM_INEXACTP (x
))
2930 if (SCM_NUMBERP (x
))
2932 SCM_WRONG_TYPE_ARG (1, x
);
2937 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
2938 /* "Return @code{#t} if all parameters are numerically equal." */
2940 scm_num_eq_p (SCM x
, SCM y
)
2944 long xx
= SCM_INUM (x
);
2947 long yy
= SCM_INUM (y
);
2948 return SCM_BOOL (xx
== yy
);
2950 else if (SCM_BIGP (y
))
2952 else if (SCM_REALP (y
))
2953 return SCM_BOOL ((double) xx
== SCM_REAL_VALUE (y
));
2954 else if (SCM_COMPLEXP (y
))
2955 return SCM_BOOL (((double) xx
== SCM_COMPLEX_REAL (y
))
2956 && (0.0 == SCM_COMPLEX_IMAG (y
)));
2957 else if (SCM_FRACTIONP (y
))
2960 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
2962 else if (SCM_BIGP (x
))
2966 else if (SCM_BIGP (y
))
2968 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2969 scm_remember_upto_here_2 (x
, y
);
2970 return SCM_BOOL (0 == cmp
);
2972 else if (SCM_REALP (y
))
2975 if (xisnan (SCM_REAL_VALUE (y
)))
2977 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
2978 scm_remember_upto_here_1 (x
);
2979 return SCM_BOOL (0 == cmp
);
2981 else if (SCM_COMPLEXP (y
))
2984 if (0.0 != SCM_COMPLEX_IMAG (y
))
2986 if (xisnan (SCM_COMPLEX_REAL (y
)))
2988 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
2989 scm_remember_upto_here_1 (x
);
2990 return SCM_BOOL (0 == cmp
);
2992 else if (SCM_FRACTIONP (y
))
2995 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
2997 else if (SCM_REALP (x
))
3000 return SCM_BOOL (SCM_REAL_VALUE (x
) == (double) SCM_INUM (y
));
3001 else if (SCM_BIGP (y
))
3004 if (xisnan (SCM_REAL_VALUE (x
)))
3006 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3007 scm_remember_upto_here_1 (y
);
3008 return SCM_BOOL (0 == cmp
);
3010 else if (SCM_REALP (y
))
3011 return SCM_BOOL (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3012 else if (SCM_COMPLEXP (y
))
3013 return SCM_BOOL ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3014 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3015 else if (SCM_FRACTIONP (y
))
3016 return SCM_BOOL (SCM_REAL_VALUE (x
) == scm_i_fraction2double (y
));
3018 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3020 else if (SCM_COMPLEXP (x
))
3023 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == (double) SCM_INUM (y
))
3024 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3025 else if (SCM_BIGP (y
))
3028 if (0.0 != SCM_COMPLEX_IMAG (x
))
3030 if (xisnan (SCM_COMPLEX_REAL (x
)))
3032 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3033 scm_remember_upto_here_1 (y
);
3034 return SCM_BOOL (0 == cmp
);
3036 else if (SCM_REALP (y
))
3037 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3038 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3039 else if (SCM_COMPLEXP (y
))
3040 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3041 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3042 else if (SCM_FRACTIONP (y
))
3043 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == scm_i_fraction2double (y
))
3044 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3046 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3048 else if (SCM_FRACTIONP (x
))
3052 else if (SCM_BIGP (y
))
3054 else if (SCM_REALP (y
))
3055 return SCM_BOOL (scm_i_fraction2double (x
) == SCM_REAL_VALUE (y
));
3056 else if (SCM_COMPLEXP (y
))
3057 return SCM_BOOL ((scm_i_fraction2double (x
) == SCM_COMPLEX_REAL (y
))
3058 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3059 else if (SCM_FRACTIONP (y
))
3060 return scm_i_fraction_equalp (x
, y
);
3062 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3065 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3069 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3070 /* "Return @code{#t} if the list of parameters is monotonically\n"
3074 scm_less_p (SCM x
, SCM y
)
3078 long xx
= SCM_INUM (x
);
3081 long yy
= SCM_INUM (y
);
3082 return SCM_BOOL (xx
< yy
);
3084 else if (SCM_BIGP (y
))
3086 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3087 scm_remember_upto_here_1 (y
);
3088 return SCM_BOOL (sgn
> 0);
3090 else if (SCM_REALP (y
))
3091 return SCM_BOOL ((double) xx
< SCM_REAL_VALUE (y
));
3092 else if (SCM_FRACTIONP (y
))
3093 return SCM_BOOL ((double) xx
< scm_i_fraction2double (y
));
3095 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3097 else if (SCM_BIGP (x
))
3101 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3102 scm_remember_upto_here_1 (x
);
3103 return SCM_BOOL (sgn
< 0);
3105 else if (SCM_BIGP (y
))
3107 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3108 scm_remember_upto_here_2 (x
, y
);
3109 return SCM_BOOL (cmp
< 0);
3111 else if (SCM_REALP (y
))
3114 if (xisnan (SCM_REAL_VALUE (y
)))
3116 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3117 scm_remember_upto_here_1 (x
);
3118 return SCM_BOOL (cmp
< 0);
3120 else if (SCM_FRACTIONP (y
))
3123 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), scm_i_fraction2double (y
));
3124 scm_remember_upto_here_1 (x
);
3125 return SCM_BOOL (cmp
< 0);
3128 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3130 else if (SCM_REALP (x
))
3133 return SCM_BOOL (SCM_REAL_VALUE (x
) < (double) SCM_INUM (y
));
3134 else if (SCM_BIGP (y
))
3137 if (xisnan (SCM_REAL_VALUE (x
)))
3139 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3140 scm_remember_upto_here_1 (y
);
3141 return SCM_BOOL (cmp
> 0);
3143 else if (SCM_REALP (y
))
3144 return SCM_BOOL (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3145 else if (SCM_FRACTIONP (y
))
3146 return SCM_BOOL (SCM_REAL_VALUE (x
) < scm_i_fraction2double (y
));
3148 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3150 else if (SCM_FRACTIONP (x
))
3153 return SCM_BOOL (scm_i_fraction2double (x
) < (double) SCM_INUM (y
));
3154 else if (SCM_BIGP (y
))
3157 if (xisnan (SCM_REAL_VALUE (x
)))
3159 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), scm_i_fraction2double (x
));
3160 scm_remember_upto_here_1 (y
);
3161 return SCM_BOOL (cmp
> 0);
3163 else if (SCM_REALP (y
))
3164 return SCM_BOOL (scm_i_fraction2double (x
) < SCM_REAL_VALUE (y
));
3165 else if (SCM_FRACTIONP (y
))
3166 return SCM_BOOL (scm_i_fraction2double (x
) < scm_i_fraction2double (y
));
3168 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3171 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3175 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3176 /* "Return @code{#t} if the list of parameters is monotonically\n"
3179 #define FUNC_NAME s_scm_gr_p
3181 scm_gr_p (SCM x
, SCM y
)
3183 if (!SCM_NUMBERP (x
))
3184 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3185 else if (!SCM_NUMBERP (y
))
3186 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3188 return scm_less_p (y
, x
);
3193 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3194 /* "Return @code{#t} if the list of parameters is monotonically\n"
3197 #define FUNC_NAME s_scm_leq_p
3199 scm_leq_p (SCM x
, SCM y
)
3201 if (!SCM_NUMBERP (x
))
3202 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3203 else if (!SCM_NUMBERP (y
))
3204 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3205 else if (SCM_NFALSEP (scm_nan_p (x
)) || SCM_NFALSEP (scm_nan_p (y
)))
3208 return SCM_BOOL_NOT (scm_less_p (y
, x
));
3213 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3214 /* "Return @code{#t} if the list of parameters is monotonically\n"
3217 #define FUNC_NAME s_scm_geq_p
3219 scm_geq_p (SCM x
, SCM y
)
3221 if (!SCM_NUMBERP (x
))
3222 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3223 else if (!SCM_NUMBERP (y
))
3224 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3225 else if (SCM_NFALSEP (scm_nan_p (x
)) || SCM_NFALSEP (scm_nan_p (y
)))
3228 return SCM_BOOL_NOT (scm_less_p (x
, y
));
3233 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3234 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3241 return SCM_BOOL (SCM_EQ_P (z
, SCM_INUM0
));
3242 else if (SCM_BIGP (z
))
3244 else if (SCM_REALP (z
))
3245 return SCM_BOOL (SCM_REAL_VALUE (z
) == 0.0);
3246 else if (SCM_COMPLEXP (z
))
3247 return SCM_BOOL (SCM_COMPLEX_REAL (z
) == 0.0
3248 && SCM_COMPLEX_IMAG (z
) == 0.0);
3249 else if (SCM_FRACTIONP (z
))
3252 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3256 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3257 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3261 scm_positive_p (SCM x
)
3264 return SCM_BOOL (SCM_INUM (x
) > 0);
3265 else if (SCM_BIGP (x
))
3267 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3268 scm_remember_upto_here_1 (x
);
3269 return SCM_BOOL (sgn
> 0);
3271 else if (SCM_REALP (x
))
3272 return SCM_BOOL(SCM_REAL_VALUE (x
) > 0.0);
3273 else if (SCM_FRACTIONP (x
))
3274 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3276 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3280 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3281 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3285 scm_negative_p (SCM x
)
3288 return SCM_BOOL (SCM_INUM (x
) < 0);
3289 else if (SCM_BIGP (x
))
3291 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3292 scm_remember_upto_here_1 (x
);
3293 return SCM_BOOL (sgn
< 0);
3295 else if (SCM_REALP (x
))
3296 return SCM_BOOL(SCM_REAL_VALUE (x
) < 0.0);
3297 else if (SCM_FRACTIONP (x
))
3298 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3300 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3304 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3305 /* "Return the maximum of all parameter values."
3308 scm_max (SCM x
, SCM y
)
3313 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3314 else if (SCM_NUMBERP (x
))
3317 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3322 long xx
= SCM_INUM (x
);
3325 long yy
= SCM_INUM (y
);
3326 return (xx
< yy
) ? y
: x
;
3328 else if (SCM_BIGP (y
))
3330 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3331 scm_remember_upto_here_1 (y
);
3332 return (sgn
< 0) ? x
: y
;
3334 else if (SCM_REALP (y
))
3337 /* if y==NaN then ">" is false and we return NaN */
3338 return (z
> SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3340 else if (SCM_FRACTIONP (y
))
3343 return (z
> scm_i_fraction2double (y
)) ? x
: y
;
3346 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3348 else if (SCM_BIGP (x
))
3352 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3353 scm_remember_upto_here_1 (x
);
3354 return (sgn
< 0) ? y
: x
;
3356 else if (SCM_BIGP (y
))
3358 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3359 scm_remember_upto_here_2 (x
, y
);
3360 return (cmp
> 0) ? x
: y
;
3362 else if (SCM_REALP (y
))
3364 double yy
= SCM_REAL_VALUE (y
);
3368 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3369 scm_remember_upto_here_1 (x
);
3370 return (cmp
> 0) ? x
: y
;
3372 else if (SCM_FRACTIONP (y
))
3374 double yy
= scm_i_fraction2double (y
);
3376 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3377 scm_remember_upto_here_1 (x
);
3378 return (cmp
> 0) ? x
: y
;
3381 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3383 else if (SCM_REALP (x
))
3387 double z
= SCM_INUM (y
);
3388 /* if x==NaN then "<" is false and we return NaN */
3389 return (SCM_REAL_VALUE (x
) < z
) ? scm_make_real (z
) : x
;
3391 else if (SCM_BIGP (y
))
3393 double xx
= SCM_REAL_VALUE (x
);
3397 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3398 scm_remember_upto_here_1 (y
);
3399 return (cmp
< 0) ? x
: y
;
3401 else if (SCM_REALP (y
))
3403 /* if x==NaN then our explicit check means we return NaN
3404 if y==NaN then ">" is false and we return NaN
3405 calling isnan is unavoidable, since it's the only way to know
3406 which of x or y causes any compares to be false */
3407 double xx
= SCM_REAL_VALUE (x
);
3408 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3410 else if (SCM_FRACTIONP (y
))
3412 double yy
= scm_i_fraction2double (y
);
3413 double xx
= SCM_REAL_VALUE (x
);
3414 return (xx
< yy
) ? scm_make_real (yy
) : x
;
3417 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3419 else if (SCM_FRACTIONP (x
))
3423 double z
= SCM_INUM (y
);
3424 return (scm_i_fraction2double (x
) < z
) ? y
: x
;
3426 else if (SCM_BIGP (y
))
3428 double xx
= scm_i_fraction2double (x
);
3430 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3431 scm_remember_upto_here_1 (y
);
3432 return (cmp
< 0) ? x
: y
;
3434 else if (SCM_REALP (y
))
3436 double xx
= scm_i_fraction2double (x
);
3437 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_make_real (xx
);
3439 else if (SCM_FRACTIONP (y
))
3441 double yy
= scm_i_fraction2double (y
);
3442 double xx
= scm_i_fraction2double (x
);
3443 return (xx
< yy
) ? y
: x
;
3446 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3449 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3453 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3454 /* "Return the minium of all parameter values."
3457 scm_min (SCM x
, SCM y
)
3462 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3463 else if (SCM_NUMBERP (x
))
3466 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3471 long xx
= SCM_INUM (x
);
3474 long yy
= SCM_INUM (y
);
3475 return (xx
< yy
) ? x
: y
;
3477 else if (SCM_BIGP (y
))
3479 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3480 scm_remember_upto_here_1 (y
);
3481 return (sgn
< 0) ? y
: x
;
3483 else if (SCM_REALP (y
))
3486 /* if y==NaN then "<" is false and we return NaN */
3487 return (z
< SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3489 else if (SCM_FRACTIONP (y
))
3492 return (z
< scm_i_fraction2double (y
)) ? x
: y
;
3495 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3497 else if (SCM_BIGP (x
))
3501 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3502 scm_remember_upto_here_1 (x
);
3503 return (sgn
< 0) ? x
: y
;
3505 else if (SCM_BIGP (y
))
3507 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3508 scm_remember_upto_here_2 (x
, y
);
3509 return (cmp
> 0) ? y
: x
;
3511 else if (SCM_REALP (y
))
3513 double yy
= SCM_REAL_VALUE (y
);
3517 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3518 scm_remember_upto_here_1 (x
);
3519 return (cmp
> 0) ? y
: x
;
3521 else if (SCM_FRACTIONP (y
))
3523 double yy
= scm_i_fraction2double (y
);
3525 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3526 scm_remember_upto_here_1 (x
);
3527 return (cmp
> 0) ? y
: x
;
3530 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3532 else if (SCM_REALP (x
))
3536 double z
= SCM_INUM (y
);
3537 /* if x==NaN then "<" is false and we return NaN */
3538 return (z
< SCM_REAL_VALUE (x
)) ? scm_make_real (z
) : x
;
3540 else if (SCM_BIGP (y
))
3542 double xx
= SCM_REAL_VALUE (x
);
3546 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3547 scm_remember_upto_here_1 (y
);
3548 return (cmp
< 0) ? y
: x
;
3550 else if (SCM_REALP (y
))
3552 /* if x==NaN then our explicit check means we return NaN
3553 if y==NaN then "<" is false and we return NaN
3554 calling isnan is unavoidable, since it's the only way to know
3555 which of x or y causes any compares to be false */
3556 double xx
= SCM_REAL_VALUE (x
);
3557 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3559 else if (SCM_FRACTIONP (y
))
3561 double yy
= scm_i_fraction2double (y
);
3562 double xx
= SCM_REAL_VALUE (x
);
3563 return (yy
< xx
) ? scm_make_real (yy
) : x
;
3566 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3568 else if (SCM_FRACTIONP (x
))
3572 double z
= SCM_INUM (y
);
3573 return (scm_i_fraction2double (x
) < z
) ? x
: y
;
3575 else if (SCM_BIGP (y
))
3577 double xx
= scm_i_fraction2double (x
);
3579 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3580 scm_remember_upto_here_1 (y
);
3581 return (cmp
< 0) ? y
: x
;
3583 else if (SCM_REALP (y
))
3585 double xx
= scm_i_fraction2double (x
);
3586 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_make_real (xx
);
3588 else if (SCM_FRACTIONP (y
))
3590 double yy
= scm_i_fraction2double (y
);
3591 double xx
= scm_i_fraction2double (x
);
3592 return (xx
< yy
) ? x
: y
;
3595 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3598 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3602 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3603 /* "Return the sum of all parameter values. Return 0 if called without\n"
3607 scm_sum (SCM x
, SCM y
)
3611 if (SCM_NUMBERP (x
)) return x
;
3612 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3613 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3620 long xx
= SCM_INUM (x
);
3621 long yy
= SCM_INUM (y
);
3622 long int z
= xx
+ yy
;
3623 return SCM_FIXABLE (z
) ? SCM_MAKINUM (z
) : scm_i_long2big (z
);
3625 else if (SCM_BIGP (y
))
3630 else if (SCM_REALP (y
))
3632 long int xx
= SCM_INUM (x
);
3633 return scm_make_real (xx
+ SCM_REAL_VALUE (y
));
3635 else if (SCM_COMPLEXP (y
))
3637 long int xx
= SCM_INUM (x
);
3638 return scm_make_complex (xx
+ SCM_COMPLEX_REAL (y
),
3639 SCM_COMPLEX_IMAG (y
));
3641 else if (SCM_FRACTIONP (y
))
3642 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3643 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3644 SCM_FRACTION_DENOMINATOR (y
));
3646 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3647 } else if (SCM_BIGP (x
))
3654 inum
= SCM_INUM (y
);
3657 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3660 SCM result
= scm_i_mkbig ();
3661 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3662 scm_remember_upto_here_1 (x
);
3663 /* we know the result will have to be a bignum */
3666 return scm_i_normbig (result
);
3670 SCM result
= scm_i_mkbig ();
3671 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3672 scm_remember_upto_here_1 (x
);
3673 /* we know the result will have to be a bignum */
3676 return scm_i_normbig (result
);
3679 else if (SCM_BIGP (y
))
3681 SCM result
= scm_i_mkbig ();
3682 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3683 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3684 mpz_add (SCM_I_BIG_MPZ (result
),
3687 scm_remember_upto_here_2 (x
, y
);
3688 /* we know the result will have to be a bignum */
3691 return scm_i_normbig (result
);
3693 else if (SCM_REALP (y
))
3695 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3696 scm_remember_upto_here_1 (x
);
3697 return scm_make_real (result
);
3699 else if (SCM_COMPLEXP (y
))
3701 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3702 + SCM_COMPLEX_REAL (y
));
3703 scm_remember_upto_here_1 (x
);
3704 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
3706 else if (SCM_FRACTIONP (y
))
3707 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3708 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3709 SCM_FRACTION_DENOMINATOR (y
));
3711 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3713 else if (SCM_REALP (x
))
3716 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_INUM (y
));
3717 else if (SCM_BIGP (y
))
3719 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3720 scm_remember_upto_here_1 (y
);
3721 return scm_make_real (result
);
3723 else if (SCM_REALP (y
))
3724 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3725 else if (SCM_COMPLEXP (y
))
3726 return scm_make_complex (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3727 SCM_COMPLEX_IMAG (y
));
3728 else if (SCM_FRACTIONP (y
))
3729 return scm_make_real (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3731 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3733 else if (SCM_COMPLEXP (x
))
3736 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_INUM (y
),
3737 SCM_COMPLEX_IMAG (x
));
3738 else if (SCM_BIGP (y
))
3740 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3741 + SCM_COMPLEX_REAL (x
));
3742 scm_remember_upto_here_1 (y
);
3743 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (x
));
3745 else if (SCM_REALP (y
))
3746 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3747 SCM_COMPLEX_IMAG (x
));
3748 else if (SCM_COMPLEXP (y
))
3749 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3750 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3751 else if (SCM_FRACTIONP (y
))
3752 return scm_make_complex (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3753 SCM_COMPLEX_IMAG (x
));
3755 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3757 else if (SCM_FRACTIONP (x
))
3760 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3761 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3762 SCM_FRACTION_DENOMINATOR (x
));
3763 else if (SCM_BIGP (y
))
3764 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3765 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3766 SCM_FRACTION_DENOMINATOR (x
));
3767 else if (SCM_REALP (y
))
3768 return scm_make_real (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
3769 else if (SCM_COMPLEXP (y
))
3770 return scm_make_complex (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
3771 SCM_COMPLEX_IMAG (y
));
3772 else if (SCM_FRACTIONP (y
))
3773 /* a/b + c/d = (ad + bc) / bd */
3774 return scm_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
3775 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
3776 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
3778 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3781 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
3785 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
3786 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
3787 * the sum of all but the first argument are subtracted from the first
3789 #define FUNC_NAME s_difference
3791 scm_difference (SCM x
, SCM y
)
3796 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
3800 long xx
= -SCM_INUM (x
);
3801 if (SCM_FIXABLE (xx
))
3802 return SCM_MAKINUM (xx
);
3804 return scm_i_long2big (xx
);
3806 else if (SCM_BIGP (x
))
3807 /* FIXME: do we really need to normalize here? */
3808 return scm_i_normbig (scm_i_clonebig (x
, 0));
3809 else if (SCM_REALP (x
))
3810 return scm_make_real (-SCM_REAL_VALUE (x
));
3811 else if (SCM_COMPLEXP (x
))
3812 return scm_make_complex (-SCM_COMPLEX_REAL (x
),
3813 -SCM_COMPLEX_IMAG (x
));
3814 else if (SCM_FRACTIONP (x
))
3815 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
3816 SCM_FRACTION_DENOMINATOR (x
));
3818 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
3825 long int xx
= SCM_INUM (x
);
3826 long int yy
= SCM_INUM (y
);
3827 long int z
= xx
- yy
;
3828 if (SCM_FIXABLE (z
))
3829 return SCM_MAKINUM (z
);
3831 return scm_i_long2big (z
);
3833 else if (SCM_BIGP (y
))
3835 /* inum-x - big-y */
3836 long xx
= SCM_INUM (x
);
3839 return scm_i_clonebig (y
, 0);
3842 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3843 SCM result
= scm_i_mkbig ();
3846 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
3849 /* x - y == -(y + -x) */
3850 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
3851 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
3853 scm_remember_upto_here_1 (y
);
3855 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
3856 /* we know the result will have to be a bignum */
3859 return scm_i_normbig (result
);
3862 else if (SCM_REALP (y
))
3864 long int xx
= SCM_INUM (x
);
3865 return scm_make_real (xx
- SCM_REAL_VALUE (y
));
3867 else if (SCM_COMPLEXP (y
))
3869 long int xx
= SCM_INUM (x
);
3870 return scm_make_complex (xx
- SCM_COMPLEX_REAL (y
),
3871 - SCM_COMPLEX_IMAG (y
));
3873 else if (SCM_FRACTIONP (y
))
3874 /* a - b/c = (ac - b) / c */
3875 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
3876 SCM_FRACTION_NUMERATOR (y
)),
3877 SCM_FRACTION_DENOMINATOR (y
));
3879 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3881 else if (SCM_BIGP (x
))
3885 /* big-x - inum-y */
3886 long yy
= SCM_INUM (y
);
3887 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3889 scm_remember_upto_here_1 (x
);
3891 return SCM_FIXABLE (-yy
) ? SCM_MAKINUM (-yy
) : scm_long2num (-yy
);
3894 SCM result
= scm_i_mkbig ();
3897 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
3899 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
3900 scm_remember_upto_here_1 (x
);
3902 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
3903 /* we know the result will have to be a bignum */
3906 return scm_i_normbig (result
);
3909 else if (SCM_BIGP (y
))
3911 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3912 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3913 SCM result
= scm_i_mkbig ();
3914 mpz_sub (SCM_I_BIG_MPZ (result
),
3917 scm_remember_upto_here_2 (x
, y
);
3918 /* we know the result will have to be a bignum */
3919 if ((sgn_x
== 1) && (sgn_y
== -1))
3921 if ((sgn_x
== -1) && (sgn_y
== 1))
3923 return scm_i_normbig (result
);
3925 else if (SCM_REALP (y
))
3927 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
3928 scm_remember_upto_here_1 (x
);
3929 return scm_make_real (result
);
3931 else if (SCM_COMPLEXP (y
))
3933 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3934 - SCM_COMPLEX_REAL (y
));
3935 scm_remember_upto_here_1 (x
);
3936 return scm_make_complex (real_part
, - SCM_COMPLEX_IMAG (y
));
3938 else if (SCM_FRACTIONP (y
))
3939 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
3940 SCM_FRACTION_NUMERATOR (y
)),
3941 SCM_FRACTION_DENOMINATOR (y
));
3942 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3944 else if (SCM_REALP (x
))
3947 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_INUM (y
));
3948 else if (SCM_BIGP (y
))
3950 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
3951 scm_remember_upto_here_1 (x
);
3952 return scm_make_real (result
);
3954 else if (SCM_REALP (y
))
3955 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
3956 else if (SCM_COMPLEXP (y
))
3957 return scm_make_complex (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
3958 -SCM_COMPLEX_IMAG (y
));
3959 else if (SCM_FRACTIONP (y
))
3960 return scm_make_real (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
3962 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3964 else if (SCM_COMPLEXP (x
))
3967 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_INUM (y
),
3968 SCM_COMPLEX_IMAG (x
));
3969 else if (SCM_BIGP (y
))
3971 double real_part
= (SCM_COMPLEX_REAL (x
)
3972 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
3973 scm_remember_upto_here_1 (x
);
3974 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
3976 else if (SCM_REALP (y
))
3977 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
3978 SCM_COMPLEX_IMAG (x
));
3979 else if (SCM_COMPLEXP (y
))
3980 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
3981 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
3982 else if (SCM_FRACTIONP (y
))
3983 return scm_make_complex (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
3984 SCM_COMPLEX_IMAG (x
));
3986 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3988 else if (SCM_FRACTIONP (x
))
3991 /* a/b - c = (a - cb) / b */
3992 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
3993 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
3994 SCM_FRACTION_DENOMINATOR (x
));
3995 else if (SCM_BIGP (y
))
3996 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
3997 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
3998 SCM_FRACTION_DENOMINATOR (x
));
3999 else if (SCM_REALP (y
))
4000 return scm_make_real (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4001 else if (SCM_COMPLEXP (y
))
4002 return scm_make_complex (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4003 -SCM_COMPLEX_IMAG (y
));
4004 else if (SCM_FRACTIONP (y
))
4005 /* a/b - c/d = (ad - bc) / bd */
4006 return scm_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4007 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4008 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4010 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4013 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4018 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4019 /* "Return the product of all arguments. If called without arguments,\n"
4023 scm_product (SCM x
, SCM y
)
4028 return SCM_MAKINUM (1L);
4029 else if (SCM_NUMBERP (x
))
4032 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4044 case 0: return x
; break;
4045 case 1: return y
; break;
4050 long yy
= SCM_INUM (y
);
4052 SCM k
= SCM_MAKINUM (kk
);
4053 if ((kk
== SCM_INUM (k
)) && (kk
/ xx
== yy
))
4057 SCM result
= scm_i_long2big (xx
);
4058 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4059 return scm_i_normbig (result
);
4062 else if (SCM_BIGP (y
))
4064 SCM result
= scm_i_mkbig ();
4065 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4066 scm_remember_upto_here_1 (y
);
4069 else if (SCM_REALP (y
))
4070 return scm_make_real (xx
* SCM_REAL_VALUE (y
));
4071 else if (SCM_COMPLEXP (y
))
4072 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4073 xx
* SCM_COMPLEX_IMAG (y
));
4074 else if (SCM_FRACTIONP (y
))
4075 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4076 SCM_FRACTION_DENOMINATOR (y
));
4078 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4080 else if (SCM_BIGP (x
))
4087 else if (SCM_BIGP (y
))
4089 SCM result
= scm_i_mkbig ();
4090 mpz_mul (SCM_I_BIG_MPZ (result
),
4093 scm_remember_upto_here_2 (x
, y
);
4096 else if (SCM_REALP (y
))
4098 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4099 scm_remember_upto_here_1 (x
);
4100 return scm_make_real (result
);
4102 else if (SCM_COMPLEXP (y
))
4104 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4105 scm_remember_upto_here_1 (x
);
4106 return scm_make_complex (z
* SCM_COMPLEX_REAL (y
),
4107 z
* SCM_COMPLEX_IMAG (y
));
4109 else if (SCM_FRACTIONP (y
))
4110 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4111 SCM_FRACTION_DENOMINATOR (y
));
4113 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4115 else if (SCM_REALP (x
))
4118 return scm_make_real (SCM_INUM (y
) * SCM_REAL_VALUE (x
));
4119 else if (SCM_BIGP (y
))
4121 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4122 scm_remember_upto_here_1 (y
);
4123 return scm_make_real (result
);
4125 else if (SCM_REALP (y
))
4126 return scm_make_real (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4127 else if (SCM_COMPLEXP (y
))
4128 return scm_make_complex (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4129 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4130 else if (SCM_FRACTIONP (y
))
4131 return scm_make_real (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4133 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4135 else if (SCM_COMPLEXP (x
))
4138 return scm_make_complex (SCM_INUM (y
) * SCM_COMPLEX_REAL (x
),
4139 SCM_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4140 else if (SCM_BIGP (y
))
4142 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4143 scm_remember_upto_here_1 (y
);
4144 return scm_make_complex (z
* SCM_COMPLEX_REAL (x
),
4145 z
* SCM_COMPLEX_IMAG (x
));
4147 else if (SCM_REALP (y
))
4148 return scm_make_complex (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4149 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4150 else if (SCM_COMPLEXP (y
))
4152 return scm_make_complex (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4153 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4154 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4155 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4157 else if (SCM_FRACTIONP (y
))
4159 double yy
= scm_i_fraction2double (y
);
4160 return scm_make_complex (yy
* SCM_COMPLEX_REAL (x
),
4161 yy
* SCM_COMPLEX_IMAG (x
));
4164 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4166 else if (SCM_FRACTIONP (x
))
4169 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4170 SCM_FRACTION_DENOMINATOR (x
));
4171 else if (SCM_BIGP (y
))
4172 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4173 SCM_FRACTION_DENOMINATOR (x
));
4174 else if (SCM_REALP (y
))
4175 return scm_make_real (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4176 else if (SCM_COMPLEXP (y
))
4178 double xx
= scm_i_fraction2double (x
);
4179 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4180 xx
* SCM_COMPLEX_IMAG (y
));
4182 else if (SCM_FRACTIONP (y
))
4183 /* a/b * c/d = ac / bd */
4184 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4185 SCM_FRACTION_NUMERATOR (y
)),
4186 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4187 SCM_FRACTION_DENOMINATOR (y
)));
4189 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4192 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4196 scm_num2dbl (SCM a
, const char *why
)
4197 #define FUNC_NAME why
4200 return (double) SCM_INUM (a
);
4201 else if (SCM_BIGP (a
))
4203 double result
= mpz_get_d (SCM_I_BIG_MPZ (a
));
4204 scm_remember_upto_here_1 (a
);
4207 else if (SCM_REALP (a
))
4208 return (SCM_REAL_VALUE (a
));
4209 else if (SCM_FRACTIONP (a
))
4210 return scm_i_fraction2double (a
);
4212 SCM_WRONG_TYPE_ARG (SCM_ARGn
, a
);
4216 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4217 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4218 #define ALLOW_DIVIDE_BY_ZERO
4219 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4222 /* The code below for complex division is adapted from the GNU
4223 libstdc++, which adapted it from f2c's libF77, and is subject to
4226 /****************************************************************
4227 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4229 Permission to use, copy, modify, and distribute this software
4230 and its documentation for any purpose and without fee is hereby
4231 granted, provided that the above copyright notice appear in all
4232 copies and that both that the copyright notice and this
4233 permission notice and warranty disclaimer appear in supporting
4234 documentation, and that the names of AT&T Bell Laboratories or
4235 Bellcore or any of their entities not be used in advertising or
4236 publicity pertaining to distribution of the software without
4237 specific, written prior permission.
4239 AT&T and Bellcore disclaim all warranties with regard to this
4240 software, including all implied warranties of merchantability
4241 and fitness. In no event shall AT&T or Bellcore be liable for
4242 any special, indirect or consequential damages or any damages
4243 whatsoever resulting from loss of use, data or profits, whether
4244 in an action of contract, negligence or other tortious action,
4245 arising out of or in connection with the use or performance of
4247 ****************************************************************/
4249 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4250 /* Divide the first argument by the product of the remaining
4251 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4253 #define FUNC_NAME s_divide
4255 scm_i_divide (SCM x
, SCM y
, int inexact
)
4262 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4263 else if (SCM_INUMP (x
))
4265 long xx
= SCM_INUM (x
);
4266 if (xx
== 1 || xx
== -1)
4268 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4270 scm_num_overflow (s_divide
);
4275 return scm_make_real (1.0 / (double) xx
);
4276 else return scm_make_ratio (SCM_MAKINUM(1), x
);
4279 else if (SCM_BIGP (x
))
4282 return scm_make_real (1.0 / scm_i_big2dbl (x
));
4283 else return scm_make_ratio (SCM_MAKINUM(1), x
);
4285 else if (SCM_REALP (x
))
4287 double xx
= SCM_REAL_VALUE (x
);
4288 #ifndef ALLOW_DIVIDE_BY_ZERO
4290 scm_num_overflow (s_divide
);
4293 return scm_make_real (1.0 / xx
);
4295 else if (SCM_COMPLEXP (x
))
4297 double r
= SCM_COMPLEX_REAL (x
);
4298 double i
= SCM_COMPLEX_IMAG (x
);
4302 double d
= i
* (1.0 + t
* t
);
4303 return scm_make_complex (t
/ d
, -1.0 / d
);
4308 double d
= r
* (1.0 + t
* t
);
4309 return scm_make_complex (1.0 / d
, -t
/ d
);
4312 else if (SCM_FRACTIONP (x
))
4313 return scm_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4314 SCM_FRACTION_NUMERATOR (x
));
4316 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4321 long xx
= SCM_INUM (x
);
4324 long yy
= SCM_INUM (y
);
4327 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4328 scm_num_overflow (s_divide
);
4330 return scm_make_real ((double) xx
/ (double) yy
);
4333 else if (xx
% yy
!= 0)
4336 return scm_make_real ((double) xx
/ (double) yy
);
4337 else return scm_make_ratio (x
, y
);
4342 if (SCM_FIXABLE (z
))
4343 return SCM_MAKINUM (z
);
4345 return scm_i_long2big (z
);
4348 else if (SCM_BIGP (y
))
4351 return scm_make_real ((double) xx
/ scm_i_big2dbl (y
));
4352 else return scm_make_ratio (x
, y
);
4354 else if (SCM_REALP (y
))
4356 double yy
= SCM_REAL_VALUE (y
);
4357 #ifndef ALLOW_DIVIDE_BY_ZERO
4359 scm_num_overflow (s_divide
);
4362 return scm_make_real ((double) xx
/ yy
);
4364 else if (SCM_COMPLEXP (y
))
4367 complex_div
: /* y _must_ be a complex number */
4369 double r
= SCM_COMPLEX_REAL (y
);
4370 double i
= SCM_COMPLEX_IMAG (y
);
4374 double d
= i
* (1.0 + t
* t
);
4375 return scm_make_complex ((a
* t
) / d
, -a
/ d
);
4380 double d
= r
* (1.0 + t
* t
);
4381 return scm_make_complex (a
/ d
, -(a
* t
) / d
);
4385 else if (SCM_FRACTIONP (y
))
4386 /* a / b/c = ac / b */
4387 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4388 SCM_FRACTION_NUMERATOR (y
));
4390 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4392 else if (SCM_BIGP (x
))
4396 long int yy
= SCM_INUM (y
);
4399 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4400 scm_num_overflow (s_divide
);
4402 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4403 scm_remember_upto_here_1 (x
);
4404 return (sgn
== 0) ? scm_nan () : scm_inf ();
4411 /* FIXME: HMM, what are the relative performance issues here?
4412 We need to test. Is it faster on average to test
4413 divisible_p, then perform whichever operation, or is it
4414 faster to perform the integer div opportunistically and
4415 switch to real if there's a remainder? For now we take the
4416 middle ground: test, then if divisible, use the faster div
4419 long abs_yy
= yy
< 0 ? -yy
: yy
;
4420 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4424 SCM result
= scm_i_mkbig ();
4425 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4426 scm_remember_upto_here_1 (x
);
4428 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4429 return scm_i_normbig (result
);
4434 return scm_make_real (scm_i_big2dbl (x
) / (double) yy
);
4435 else return scm_make_ratio (x
, y
);
4439 else if (SCM_BIGP (y
))
4441 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4444 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4445 scm_num_overflow (s_divide
);
4447 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4448 scm_remember_upto_here_1 (x
);
4449 return (sgn
== 0) ? scm_nan () : scm_inf ();
4455 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4459 SCM result
= scm_i_mkbig ();
4460 mpz_divexact (SCM_I_BIG_MPZ (result
),
4463 scm_remember_upto_here_2 (x
, y
);
4464 return scm_i_normbig (result
);
4470 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4471 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4472 scm_remember_upto_here_2 (x
, y
);
4473 return scm_make_real (dbx
/ dby
);
4475 else return scm_make_ratio (x
, y
);
4479 else if (SCM_REALP (y
))
4481 double yy
= SCM_REAL_VALUE (y
);
4482 #ifndef ALLOW_DIVIDE_BY_ZERO
4484 scm_num_overflow (s_divide
);
4487 return scm_make_real (scm_i_big2dbl (x
) / yy
);
4489 else if (SCM_COMPLEXP (y
))
4491 a
= scm_i_big2dbl (x
);
4494 else if (SCM_FRACTIONP (y
))
4495 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4496 SCM_FRACTION_NUMERATOR (y
));
4498 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4500 else if (SCM_REALP (x
))
4502 double rx
= SCM_REAL_VALUE (x
);
4505 long int yy
= SCM_INUM (y
);
4506 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4508 scm_num_overflow (s_divide
);
4511 return scm_make_real (rx
/ (double) yy
);
4513 else if (SCM_BIGP (y
))
4515 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4516 scm_remember_upto_here_1 (y
);
4517 return scm_make_real (rx
/ dby
);
4519 else if (SCM_REALP (y
))
4521 double yy
= SCM_REAL_VALUE (y
);
4522 #ifndef ALLOW_DIVIDE_BY_ZERO
4524 scm_num_overflow (s_divide
);
4527 return scm_make_real (rx
/ yy
);
4529 else if (SCM_COMPLEXP (y
))
4534 else if (SCM_FRACTIONP (y
))
4535 return scm_make_real (rx
/ scm_i_fraction2double (y
));
4537 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4539 else if (SCM_COMPLEXP (x
))
4541 double rx
= SCM_COMPLEX_REAL (x
);
4542 double ix
= SCM_COMPLEX_IMAG (x
);
4545 long int yy
= SCM_INUM (y
);
4546 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4548 scm_num_overflow (s_divide
);
4553 return scm_make_complex (rx
/ d
, ix
/ d
);
4556 else if (SCM_BIGP (y
))
4558 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4559 scm_remember_upto_here_1 (y
);
4560 return scm_make_complex (rx
/ dby
, ix
/ dby
);
4562 else if (SCM_REALP (y
))
4564 double yy
= SCM_REAL_VALUE (y
);
4565 #ifndef ALLOW_DIVIDE_BY_ZERO
4567 scm_num_overflow (s_divide
);
4570 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4572 else if (SCM_COMPLEXP (y
))
4574 double ry
= SCM_COMPLEX_REAL (y
);
4575 double iy
= SCM_COMPLEX_IMAG (y
);
4579 double d
= iy
* (1.0 + t
* t
);
4580 return scm_make_complex ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4585 double d
= ry
* (1.0 + t
* t
);
4586 return scm_make_complex ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4589 else if (SCM_FRACTIONP (y
))
4591 double yy
= scm_i_fraction2double (y
);
4592 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4595 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4597 else if (SCM_FRACTIONP (x
))
4601 long int yy
= SCM_INUM (y
);
4602 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4604 scm_num_overflow (s_divide
);
4607 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4608 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4610 else if (SCM_BIGP (y
))
4612 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4613 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4615 else if (SCM_REALP (y
))
4617 double yy
= SCM_REAL_VALUE (y
);
4618 #ifndef ALLOW_DIVIDE_BY_ZERO
4620 scm_num_overflow (s_divide
);
4623 return scm_make_real (scm_i_fraction2double (x
) / yy
);
4625 else if (SCM_COMPLEXP (y
))
4627 a
= scm_i_fraction2double (x
);
4630 else if (SCM_FRACTIONP (y
))
4631 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4632 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4634 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4637 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4641 scm_divide (SCM x
, SCM y
)
4643 return scm_i_divide (x
, y
, 0);
4646 static SCM
scm_divide2real (SCM x
, SCM y
)
4648 return scm_i_divide (x
, y
, 1);
4654 scm_asinh (double x
)
4659 #define asinh scm_asinh
4660 return log (x
+ sqrt (x
* x
+ 1));
4663 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4664 /* "Return the inverse hyperbolic sine of @var{x}."
4669 scm_acosh (double x
)
4674 #define acosh scm_acosh
4675 return log (x
+ sqrt (x
* x
- 1));
4678 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4679 /* "Return the inverse hyperbolic cosine of @var{x}."
4684 scm_atanh (double x
)
4689 #define atanh scm_atanh
4690 return 0.5 * log ((1 + x
) / (1 - x
));
4693 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4694 /* "Return the inverse hyperbolic tangent of @var{x}."
4698 /* XXX - eventually, we should remove this definition of scm_round and
4699 rename scm_round_number to scm_round. Likewise for scm_truncate
4700 and scm_truncate_number.
4704 scm_truncate (double x
)
4709 #define trunc scm_truncate
4717 scm_round (double x
)
4719 double plus_half
= x
+ 0.5;
4720 double result
= floor (plus_half
);
4721 /* Adjust so that the scm_round is towards even. */
4722 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4727 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4729 "Round the number @var{x} towards zero.")
4730 #define FUNC_NAME s_scm_truncate_number
4732 if (SCM_FALSEP (scm_negative_p (x
)))
4733 return scm_floor (x
);
4735 return scm_ceiling (x
);
4739 static SCM exactly_one_half
;
4741 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4743 "Round the number @var{x} towards the nearest integer. "
4744 "When it is exactly halfway between two integers, "
4745 "round towards the even one.")
4746 #define FUNC_NAME s_scm_round_number
4748 SCM plus_half
= scm_sum (x
, exactly_one_half
);
4749 SCM result
= scm_floor (plus_half
);
4750 /* Adjust so that the scm_round is towards even. */
4751 if (!SCM_FALSEP (scm_num_eq_p (plus_half
, result
))
4752 && !SCM_FALSEP (scm_odd_p (result
)))
4753 return scm_difference (result
, SCM_MAKINUM (1));
4759 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
4761 "Round the number @var{x} towards minus infinity.")
4762 #define FUNC_NAME s_scm_floor
4764 if (SCM_INUMP (x
) || SCM_BIGP (x
))
4766 else if (SCM_REALP (x
))
4767 return scm_make_real (floor (SCM_REAL_VALUE (x
)));
4768 else if (SCM_FRACTIONP (x
))
4770 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
4771 SCM_FRACTION_DENOMINATOR (x
));
4772 if (SCM_FALSEP (scm_negative_p (x
)))
4774 /* For positive x, rounding towards zero is correct. */
4779 /* For negative x, we need to return q-1 unless x is an
4780 integer. But fractions are never integer, per our
4782 return scm_difference (q
, SCM_MAKINUM (1));
4786 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
4790 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
4792 "Round the number @var{x} towards infinity.")
4793 #define FUNC_NAME s_scm_ceiling
4795 if (SCM_INUMP (x
) || SCM_BIGP (x
))
4797 else if (SCM_REALP (x
))
4798 return scm_make_real (ceil (SCM_REAL_VALUE (x
)));
4799 else if (SCM_FRACTIONP (x
))
4801 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
4802 SCM_FRACTION_DENOMINATOR (x
));
4803 if (SCM_FALSEP (scm_positive_p (x
)))
4805 /* For negative x, rounding towards zero is correct. */
4810 /* For positive x, we need to return q+1 unless x is an
4811 integer. But fractions are never integer, per our
4813 return scm_sum (q
, SCM_MAKINUM (1));
4817 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
4821 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
4822 /* "Return the square root of the real number @var{x}."
4824 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
4825 /* "Return the absolute value of the real number @var{x}."
4827 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
4828 /* "Return the @var{x}th power of e."
4830 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
4831 /* "Return the natural logarithm of the real number @var{x}."
4833 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
4834 /* "Return the sine of the real number @var{x}."
4836 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
4837 /* "Return the cosine of the real number @var{x}."
4839 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
4840 /* "Return the tangent of the real number @var{x}."
4842 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
4843 /* "Return the arc sine of the real number @var{x}."
4845 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
4846 /* "Return the arc cosine of the real number @var{x}."
4848 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
4849 /* "Return the arc tangent of the real number @var{x}."
4851 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
4852 /* "Return the hyperbolic sine of the real number @var{x}."
4854 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
4855 /* "Return the hyperbolic cosine of the real number @var{x}."
4857 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
4858 /* "Return the hyperbolic tangent of the real number @var{x}."
4866 static void scm_two_doubles (SCM x
,
4868 const char *sstring
,
4872 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
4875 xy
->x
= SCM_INUM (x
);
4876 else if (SCM_BIGP (x
))
4877 xy
->x
= scm_i_big2dbl (x
);
4878 else if (SCM_REALP (x
))
4879 xy
->x
= SCM_REAL_VALUE (x
);
4880 else if (SCM_FRACTIONP (x
))
4881 xy
->x
= scm_i_fraction2double (x
);
4883 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
4886 xy
->y
= SCM_INUM (y
);
4887 else if (SCM_BIGP (y
))
4888 xy
->y
= scm_i_big2dbl (y
);
4889 else if (SCM_REALP (y
))
4890 xy
->y
= SCM_REAL_VALUE (y
);
4891 else if (SCM_FRACTIONP (y
))
4892 xy
->y
= scm_i_fraction2double (y
);
4894 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
4898 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
4900 "Return @var{x} raised to the power of @var{y}. This\n"
4901 "procedure does not accept complex arguments.")
4902 #define FUNC_NAME s_scm_sys_expt
4905 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4906 return scm_make_real (pow (xy
.x
, xy
.y
));
4911 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
4913 "Return the arc tangent of the two arguments @var{x} and\n"
4914 "@var{y}. This is similar to calculating the arc tangent of\n"
4915 "@var{x} / @var{y}, except that the signs of both arguments\n"
4916 "are used to determine the quadrant of the result. This\n"
4917 "procedure does not accept complex arguments.")
4918 #define FUNC_NAME s_scm_sys_atan2
4921 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4922 return scm_make_real (atan2 (xy
.x
, xy
.y
));
4927 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
4928 (SCM real
, SCM imaginary
),
4929 "Return a complex number constructed of the given @var{real} and\n"
4930 "@var{imaginary} parts.")
4931 #define FUNC_NAME s_scm_make_rectangular
4934 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
4935 return scm_make_complex (xy
.x
, xy
.y
);
4941 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
4943 "Return the complex number @var{x} * e^(i * @var{y}).")
4944 #define FUNC_NAME s_scm_make_polar
4948 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4950 sincos (xy
.y
, &s
, &c
);
4955 return scm_make_complex (xy
.x
* c
, xy
.x
* s
);
4960 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
4961 /* "Return the real part of the number @var{z}."
4964 scm_real_part (SCM z
)
4968 else if (SCM_BIGP (z
))
4970 else if (SCM_REALP (z
))
4972 else if (SCM_COMPLEXP (z
))
4973 return scm_make_real (SCM_COMPLEX_REAL (z
));
4974 else if (SCM_FRACTIONP (z
))
4975 return scm_make_real (scm_i_fraction2double (z
));
4977 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
4981 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
4982 /* "Return the imaginary part of the number @var{z}."
4985 scm_imag_part (SCM z
)
4989 else if (SCM_BIGP (z
))
4991 else if (SCM_REALP (z
))
4993 else if (SCM_COMPLEXP (z
))
4994 return scm_make_real (SCM_COMPLEX_IMAG (z
));
4995 else if (SCM_FRACTIONP (z
))
4998 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5001 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5002 /* "Return the numerator of the number @var{z}."
5005 scm_numerator (SCM z
)
5009 else if (SCM_BIGP (z
))
5011 else if (SCM_FRACTIONP (z
))
5013 scm_i_fraction_reduce (z
);
5014 return SCM_FRACTION_NUMERATOR (z
);
5016 else if (SCM_REALP (z
))
5017 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5019 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5023 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5024 /* "Return the denominator of the number @var{z}."
5027 scm_denominator (SCM z
)
5030 return SCM_MAKINUM (1);
5031 else if (SCM_BIGP (z
))
5032 return SCM_MAKINUM (1);
5033 else if (SCM_FRACTIONP (z
))
5035 scm_i_fraction_reduce (z
);
5036 return SCM_FRACTION_DENOMINATOR (z
);
5038 else if (SCM_REALP (z
))
5039 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5041 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5044 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5045 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5046 * "@code{abs} for real arguments, but also allows complex numbers."
5049 scm_magnitude (SCM z
)
5053 long int zz
= SCM_INUM (z
);
5056 else if (SCM_POSFIXABLE (-zz
))
5057 return SCM_MAKINUM (-zz
);
5059 return scm_i_long2big (-zz
);
5061 else if (SCM_BIGP (z
))
5063 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5064 scm_remember_upto_here_1 (z
);
5066 return scm_i_clonebig (z
, 0);
5070 else if (SCM_REALP (z
))
5071 return scm_make_real (fabs (SCM_REAL_VALUE (z
)));
5072 else if (SCM_COMPLEXP (z
))
5073 return scm_make_real (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5074 else if (SCM_FRACTIONP (z
))
5076 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5078 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5079 SCM_FRACTION_DENOMINATOR (z
));
5082 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5086 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5087 /* "Return the angle of the complex number @var{z}."
5092 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5093 scm_flo0 to save allocating a new flonum with scm_make_real each time.
5094 But if atan2 follows the floating point rounding mode, then the value
5095 is not a constant. Maybe it'd be close enough though. */
5098 if (SCM_INUM (z
) >= 0)
5101 return scm_make_real (atan2 (0.0, -1.0));
5103 else if (SCM_BIGP (z
))
5105 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5106 scm_remember_upto_here_1 (z
);
5108 return scm_make_real (atan2 (0.0, -1.0));
5112 else if (SCM_REALP (z
))
5114 if (SCM_REAL_VALUE (z
) >= 0)
5117 return scm_make_real (atan2 (0.0, -1.0));
5119 else if (SCM_COMPLEXP (z
))
5120 return scm_make_real (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5121 else if (SCM_FRACTIONP (z
))
5123 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5125 else return scm_make_real (atan2 (0.0, -1.0));
5128 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5132 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5133 /* Convert the number @var{x} to its inexact representation.\n"
5136 scm_exact_to_inexact (SCM z
)
5139 return scm_make_real ((double) SCM_INUM (z
));
5140 else if (SCM_BIGP (z
))
5141 return scm_make_real (scm_i_big2dbl (z
));
5142 else if (SCM_FRACTIONP (z
))
5143 return scm_make_real (scm_i_fraction2double (z
));
5144 else if (SCM_INEXACTP (z
))
5147 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5151 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5153 "Return an exact number that is numerically closest to @var{z}.")
5154 #define FUNC_NAME s_scm_inexact_to_exact
5158 else if (SCM_BIGP (z
))
5160 else if (SCM_REALP (z
))
5162 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5163 SCM_OUT_OF_RANGE (1, z
);
5170 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5171 q
= scm_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5172 scm_i_mpz2num (mpq_denref (frac
)));
5174 /* When scm_make_ratio throws, we leak the memory allocated
5181 else if (SCM_FRACTIONP (z
))
5184 SCM_WRONG_TYPE_ARG (1, z
);
5188 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5190 "Return an exact number that is within @var{err} of @var{x}.")
5191 #define FUNC_NAME s_scm_rationalize
5195 else if (SCM_BIGP (x
))
5197 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5199 /* Use continued fractions to find closest ratio. All
5200 arithmetic is done with exact numbers.
5203 SCM ex
= scm_inexact_to_exact (x
);
5204 SCM int_part
= scm_floor (ex
);
5205 SCM tt
= SCM_MAKINUM (1);
5206 SCM a1
= SCM_MAKINUM (0), a2
= SCM_MAKINUM (1), a
= SCM_MAKINUM (0);
5207 SCM b1
= SCM_MAKINUM (1), b2
= SCM_MAKINUM (0), b
= SCM_MAKINUM (0);
5211 if (!SCM_FALSEP (scm_num_eq_p (ex
, int_part
)))
5214 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5215 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5217 /* We stop after a million iterations just to be absolutely sure
5218 that we don't go into an infinite loop. The process normally
5219 converges after less than a dozen iterations.
5222 err
= scm_abs (err
);
5223 while (++i
< 1000000)
5225 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5226 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5227 if (SCM_FALSEP (scm_zero_p (b
)) && /* b != 0 */
5229 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5230 err
))) /* abs(x-a/b) <= err */
5232 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5233 if (SCM_FALSEP (scm_exact_p (x
))
5234 || SCM_FALSEP (scm_exact_p (err
)))
5235 return scm_exact_to_inexact (res
);
5239 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5241 tt
= scm_floor (rx
); /* tt = floor (rx) */
5247 scm_num_overflow (s_scm_rationalize
);
5250 SCM_WRONG_TYPE_ARG (1, x
);
5254 /* if you need to change this, change test-num2integral.c as well */
5255 #if SCM_SIZEOF_LONG_LONG != 0
5257 # define ULLONG_MAX ((unsigned long long) (-1))
5258 # define LLONG_MAX ((long long) (ULLONG_MAX >> 1))
5259 # define LLONG_MIN (~LLONG_MAX)
5263 /* Parameters for creating integer conversion routines.
5265 Define the following preprocessor macros before including
5266 "libguile/num2integral.i.c":
5268 NUM2INTEGRAL - the name of the function for converting from a
5269 Scheme object to the integral type. This function will be
5270 defined when including "num2integral.i.c".
5272 INTEGRAL2NUM - the name of the function for converting from the
5273 integral type to a Scheme object. This function will be defined.
5275 INTEGRAL2BIG - the name of an internal function that createas a
5276 bignum from the integral type. This function will be defined.
5277 The name should start with "scm_i_".
5279 ITYPE - the name of the integral type.
5281 UNSIGNED - Define this to 1 when ITYPE is an unsigned type. Define
5284 UNSIGNED_ITYPE - the name of the the unsigned variant of the
5285 integral type. If you don't define this, it defaults to
5286 "unsigned ITYPE" for signed types and simply "ITYPE" for unsigned
5289 SIZEOF_ITYPE - an expression giving the size of the integral type
5290 in bytes. This expression must be computable by the
5291 preprocessor. (SIZEOF_FOO values are calculated by configure.in
5296 #define NUM2INTEGRAL scm_num2short
5297 #define INTEGRAL2NUM scm_short2num
5298 #define INTEGRAL2BIG scm_i_short2big
5301 #define SIZEOF_ITYPE SIZEOF_SHORT
5302 #include "libguile/num2integral.i.c"
5304 #define NUM2INTEGRAL scm_num2ushort
5305 #define INTEGRAL2NUM scm_ushort2num
5306 #define INTEGRAL2BIG scm_i_ushort2big
5308 #define ITYPE unsigned short
5309 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_SHORT
5310 #include "libguile/num2integral.i.c"
5312 #define NUM2INTEGRAL scm_num2int
5313 #define INTEGRAL2NUM scm_int2num
5314 #define INTEGRAL2BIG scm_i_int2big
5317 #define SIZEOF_ITYPE SIZEOF_INT
5318 #include "libguile/num2integral.i.c"
5320 #define NUM2INTEGRAL scm_num2uint
5321 #define INTEGRAL2NUM scm_uint2num
5322 #define INTEGRAL2BIG scm_i_uint2big
5324 #define ITYPE unsigned int
5325 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_INT
5326 #include "libguile/num2integral.i.c"
5328 #define NUM2INTEGRAL scm_num2long
5329 #define INTEGRAL2NUM scm_long2num
5330 #define INTEGRAL2BIG scm_i_long2big
5333 #define SIZEOF_ITYPE SIZEOF_LONG
5334 #include "libguile/num2integral.i.c"
5336 #define NUM2INTEGRAL scm_num2ulong
5337 #define INTEGRAL2NUM scm_ulong2num
5338 #define INTEGRAL2BIG scm_i_ulong2big
5340 #define ITYPE unsigned long
5341 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG
5342 #include "libguile/num2integral.i.c"
5344 #define NUM2INTEGRAL scm_num2ptrdiff
5345 #define INTEGRAL2NUM scm_ptrdiff2num
5346 #define INTEGRAL2BIG scm_i_ptrdiff2big
5348 #define ITYPE scm_t_ptrdiff
5349 #define UNSIGNED_ITYPE size_t
5350 #define SIZEOF_ITYPE SCM_SIZEOF_SCM_T_PTRDIFF
5351 #include "libguile/num2integral.i.c"
5353 #define NUM2INTEGRAL scm_num2size
5354 #define INTEGRAL2NUM scm_size2num
5355 #define INTEGRAL2BIG scm_i_size2big
5357 #define ITYPE size_t
5358 #define SIZEOF_ITYPE SIZEOF_SIZE_T
5359 #include "libguile/num2integral.i.c"
5361 #if SCM_SIZEOF_LONG_LONG != 0
5363 #ifndef ULONG_LONG_MAX
5364 #define ULONG_LONG_MAX (~0ULL)
5367 #define NUM2INTEGRAL scm_num2long_long
5368 #define INTEGRAL2NUM scm_long_long2num
5369 #define INTEGRAL2BIG scm_i_long_long2big
5371 #define ITYPE long long
5372 #define SIZEOF_ITYPE SIZEOF_LONG_LONG
5373 #include "libguile/num2integral.i.c"
5375 #define NUM2INTEGRAL scm_num2ulong_long
5376 #define INTEGRAL2NUM scm_ulong_long2num
5377 #define INTEGRAL2BIG scm_i_ulong_long2big
5379 #define ITYPE unsigned long long
5380 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG_LONG
5381 #include "libguile/num2integral.i.c"
5383 #endif /* SCM_SIZEOF_LONG_LONG != 0 */
5385 #define NUM2FLOAT scm_num2float
5386 #define FLOAT2NUM scm_float2num
5388 #include "libguile/num2float.i.c"
5390 #define NUM2FLOAT scm_num2double
5391 #define FLOAT2NUM scm_double2num
5392 #define FTYPE double
5393 #include "libguile/num2float.i.c"
5398 #define SIZE_MAX ((size_t) (-1))
5401 #define PTRDIFF_MIN \
5402 ((scm_t_ptrdiff) ((scm_t_ptrdiff) 1 \
5403 << ((sizeof (scm_t_ptrdiff) * SCM_CHAR_BIT) - 1)))
5406 #define PTRDIFF_MAX (~ PTRDIFF_MIN)
5409 #define CHECK(type, v) \
5412 if ((v) != scm_num2##type (scm_##type##2num (v), 1, "check_sanity")) \
5432 CHECK (ptrdiff
, -1);
5434 CHECK (short, SHRT_MAX
);
5435 CHECK (short, SHRT_MIN
);
5436 CHECK (ushort
, USHRT_MAX
);
5437 CHECK (int, INT_MAX
);
5438 CHECK (int, INT_MIN
);
5439 CHECK (uint
, UINT_MAX
);
5440 CHECK (long, LONG_MAX
);
5441 CHECK (long, LONG_MIN
);
5442 CHECK (ulong
, ULONG_MAX
);
5443 CHECK (size
, SIZE_MAX
);
5444 CHECK (ptrdiff
, PTRDIFF_MAX
);
5445 CHECK (ptrdiff
, PTRDIFF_MIN
);
5447 #if SCM_SIZEOF_LONG_LONG != 0
5448 CHECK (long_long
, 0LL);
5449 CHECK (ulong_long
, 0ULL);
5450 CHECK (long_long
, -1LL);
5451 CHECK (long_long
, LLONG_MAX
);
5452 CHECK (long_long
, LLONG_MIN
);
5453 CHECK (ulong_long
, ULLONG_MAX
);
5460 scm_internal_catch (SCM_BOOL_T, check_body, &data, check_handler, &data); \
5461 if (!SCM_FALSEP (data)) abort();
5464 check_body (void *data
)
5466 SCM num
= *(SCM
*) data
;
5467 scm_num2ulong (num
, 1, NULL
);
5469 return SCM_UNSPECIFIED
;
5473 check_handler (void *data
, SCM tag
, SCM throw_args
)
5475 SCM
*num
= (SCM
*) data
;
5478 return SCM_UNSPECIFIED
;
5481 SCM_DEFINE (scm_sys_check_number_conversions
, "%check-number-conversions", 0, 0, 0,
5483 "Number conversion sanity checking.")
5484 #define FUNC_NAME s_scm_sys_check_number_conversions
5486 SCM data
= SCM_MAKINUM (-1);
5488 data
= scm_int2num (INT_MIN
);
5490 data
= scm_ulong2num (ULONG_MAX
);
5491 data
= scm_difference (SCM_INUM0
, data
);
5493 data
= scm_ulong2num (ULONG_MAX
);
5494 data
= scm_sum (SCM_MAKINUM (1), data
); data
= scm_difference (SCM_INUM0
, data
);
5496 data
= scm_int2num (-10000); data
= scm_product (data
, data
); data
= scm_product (data
, data
);
5499 return SCM_UNSPECIFIED
;
5508 abs_most_negative_fixnum
= scm_i_long2big (- SCM_MOST_NEGATIVE_FIXNUM
);
5509 scm_permanent_object (abs_most_negative_fixnum
);
5511 mpz_init_set_si (z_negative_one
, -1);
5513 /* It may be possible to tune the performance of some algorithms by using
5514 * the following constants to avoid the creation of bignums. Please, before
5515 * using these values, remember the two rules of program optimization:
5516 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5517 scm_c_define ("most-positive-fixnum",
5518 SCM_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5519 scm_c_define ("most-negative-fixnum",
5520 SCM_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5522 scm_add_feature ("complex");
5523 scm_add_feature ("inexact");
5524 scm_flo0
= scm_make_real (0.0);
5526 scm_dblprec
= (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5528 { /* determine floating point precision */
5530 double fsum
= 1.0 + f
;
5533 if (++scm_dblprec
> 20)
5541 scm_dblprec
= scm_dblprec
- 1;
5543 #endif /* DBL_DIG */
5549 exactly_one_half
= scm_permanent_object (scm_divide (SCM_MAKINUM (1),
5551 #include "libguile/numbers.x"