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
)
330 return scm_divide2real(numerator
, denominator
);
332 #define FUNC_NAME "make-ratio"
333 if (SCM_INUMP (denominator
))
335 if (SCM_EQ_P (denominator
, SCM_INUM0
))
336 scm_num_overflow ("make-ratio");
337 if (SCM_EQ_P (denominator
, SCM_MAKINUM(1)))
342 if (!(SCM_BIGP(denominator
)))
343 SCM_WRONG_TYPE_ARG (2, denominator
);
345 if (SCM_INUMP (numerator
))
347 if (SCM_EQ_P (numerator
, SCM_INUM0
))
349 if (SCM_INUMP (denominator
))
352 x
= SCM_INUM (numerator
);
353 y
= SCM_INUM (denominator
);
355 return SCM_MAKINUM(1);
357 return SCM_MAKINUM (x
/ y
);
359 return scm_double_cell (scm_tc16_fraction
, (scm_t_bits
)SCM_MAKINUM(-x
), (scm_t_bits
)SCM_MAKINUM(-y
), 0);
360 else return scm_double_cell (scm_tc16_fraction
, (scm_t_bits
)numerator
, (scm_t_bits
)denominator
, 0);
364 /* I assume bignums are actually big, so here there's no point in looking for a integer */
365 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (denominator
));
366 if (sgn
< 0) /* if denominator negative, flip signs */
367 return scm_double_cell (scm_tc16_fraction
,
368 (scm_t_bits
)scm_difference (numerator
, SCM_UNDEFINED
),
369 (scm_t_bits
)scm_difference (denominator
, SCM_UNDEFINED
),
371 else return scm_double_cell (scm_tc16_fraction
, (scm_t_bits
)numerator
, (scm_t_bits
)denominator
, 0);
373 /* should this use SCM_UNPACK for the bignums? */
378 if (SCM_BIGP (numerator
))
380 /* can't use scm_divide to find integer here */
381 if (SCM_INUMP (denominator
))
383 long yy
= SCM_INUM (denominator
);
384 long abs_yy
= yy
< 0 ? -yy
: yy
;
385 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), abs_yy
);
387 return scm_divide(numerator
, denominator
);
388 else return scm_double_cell (scm_tc16_fraction
, (scm_t_bits
)numerator
, (scm_t_bits
)denominator
, 0);
392 /* both are bignums */
393 if (SCM_EQ_P (numerator
, denominator
))
394 return SCM_MAKINUM(1);
395 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
396 SCM_I_BIG_MPZ (denominator
));
398 return scm_divide(numerator
, denominator
);
399 else return scm_double_cell (scm_tc16_fraction
, (scm_t_bits
)numerator
, (scm_t_bits
)denominator
, 0);
402 else SCM_WRONG_TYPE_ARG (1, numerator
);
404 return SCM_BOOL_F
; /* won't happen */
409 static void scm_i_fraction_reduce (SCM z
)
411 if (!(SCM_FRACTION_REDUCED (z
)))
414 divisor
= scm_gcd (SCM_FRACTION_NUMERATOR (z
), SCM_FRACTION_DENOMINATOR (z
));
415 if (!(SCM_EQ_P (divisor
, SCM_MAKINUM(1))))
418 SCM_FRACTION_SET_NUMERATOR (z
, scm_divide (SCM_FRACTION_NUMERATOR (z
), divisor
));
419 SCM_FRACTION_SET_DENOMINATOR (z
, scm_divide (SCM_FRACTION_DENOMINATOR (z
), divisor
));
421 SCM_FRACTION_REDUCED_SET (z
);
426 scm_i_fraction2double (SCM z
)
428 return scm_num2dbl (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
429 SCM_FRACTION_DENOMINATOR (z
)),
433 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
435 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
437 #define FUNC_NAME s_scm_exact_p
443 if (SCM_FRACTIONP (x
))
447 SCM_WRONG_TYPE_ARG (1, x
);
452 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
454 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
456 #define FUNC_NAME s_scm_odd_p
460 long val
= SCM_INUM (n
);
461 return SCM_BOOL ((val
& 1L) != 0);
463 else if (SCM_BIGP (n
))
465 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
466 scm_remember_upto_here_1 (n
);
467 return SCM_BOOL (odd_p
);
469 else if (!SCM_FALSEP (scm_inf_p (n
)))
471 else if (SCM_REALP (n
))
473 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
479 SCM_WRONG_TYPE_ARG (1, n
);
482 SCM_WRONG_TYPE_ARG (1, n
);
487 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
489 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
491 #define FUNC_NAME s_scm_even_p
495 long val
= SCM_INUM (n
);
496 return SCM_BOOL ((val
& 1L) == 0);
498 else if (SCM_BIGP (n
))
500 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
501 scm_remember_upto_here_1 (n
);
502 return SCM_BOOL (even_p
);
504 else if (!SCM_FALSEP (scm_inf_p (n
)))
506 else if (SCM_REALP (n
))
508 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
514 SCM_WRONG_TYPE_ARG (1, n
);
517 SCM_WRONG_TYPE_ARG (1, n
);
521 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
523 "Return @code{#t} if @var{n} is infinite, @code{#f}\n"
525 #define FUNC_NAME s_scm_inf_p
528 return SCM_BOOL (xisinf (SCM_REAL_VALUE (n
)));
529 else if (SCM_COMPLEXP (n
))
530 return SCM_BOOL (xisinf (SCM_COMPLEX_REAL (n
))
531 || xisinf (SCM_COMPLEX_IMAG (n
)));
537 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
539 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
541 #define FUNC_NAME s_scm_nan_p
544 return SCM_BOOL (xisnan (SCM_REAL_VALUE (n
)));
545 else if (SCM_COMPLEXP (n
))
546 return SCM_BOOL (xisnan (SCM_COMPLEX_REAL (n
))
547 || xisnan (SCM_COMPLEX_IMAG (n
)));
553 /* Guile's idea of infinity. */
554 static double guile_Inf
;
556 /* Guile's idea of not a number. */
557 static double guile_NaN
;
560 guile_ieee_init (void)
562 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
564 /* Some version of gcc on some old version of Linux used to crash when
565 trying to make Inf and NaN. */
569 guile_Inf
= 1.0 / (tmp
- tmp
);
570 #elif defined (__alpha__) && ! defined (linux)
571 extern unsigned int DINFINITY
[2];
572 guile_Inf
= (*(X_CAST(double *, DINFINITY
)));
579 if (guile_Inf
== tmp
)
587 #if defined (HAVE_ISNAN)
589 #if defined (__alpha__) && ! defined (linux)
590 extern unsigned int DQNAN
[2];
591 guile_NaN
= (*(X_CAST(double *, DQNAN
)));
593 guile_NaN
= guile_Inf
/ guile_Inf
;
599 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
602 #define FUNC_NAME s_scm_inf
604 static int initialized
= 0;
610 return scm_make_real (guile_Inf
);
614 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
617 #define FUNC_NAME s_scm_nan
619 static int initialized
= 0;
625 return scm_make_real (guile_NaN
);
630 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
632 "Return the absolute value of @var{x}.")
637 long int xx
= SCM_INUM (x
);
640 else if (SCM_POSFIXABLE (-xx
))
641 return SCM_MAKINUM (-xx
);
643 return scm_i_long2big (-xx
);
645 else if (SCM_BIGP (x
))
647 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
649 return scm_i_clonebig (x
, 0);
653 else if (SCM_REALP (x
))
654 return scm_make_real (fabs (SCM_REAL_VALUE (x
)));
655 else if (SCM_FRACTIONP (x
))
657 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
659 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
660 SCM_FRACTION_DENOMINATOR (x
));
663 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
668 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
669 /* "Return the quotient of the numbers @var{x} and @var{y}."
672 scm_quotient (SCM x
, SCM y
)
676 long xx
= SCM_INUM (x
);
679 long yy
= SCM_INUM (y
);
681 scm_num_overflow (s_quotient
);
686 return SCM_MAKINUM (z
);
688 return scm_i_long2big (z
);
691 else if (SCM_BIGP (y
))
693 if ((SCM_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
694 && (scm_i_bigcmp (abs_most_negative_fixnum
, y
) == 0))
695 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
696 return SCM_MAKINUM (-1);
698 return SCM_MAKINUM (0);
701 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
703 else if (SCM_BIGP (x
))
707 long yy
= SCM_INUM (y
);
709 scm_num_overflow (s_quotient
);
714 SCM result
= scm_i_mkbig ();
717 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
720 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
723 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
724 scm_remember_upto_here_1 (x
);
725 return scm_i_normbig (result
);
728 else if (SCM_BIGP (y
))
730 SCM result
= scm_i_mkbig ();
731 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
734 scm_remember_upto_here_2 (x
, y
);
735 return scm_i_normbig (result
);
738 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
741 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
744 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
745 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
747 * "(remainder 13 4) @result{} 1\n"
748 * "(remainder -13 4) @result{} -1\n"
752 scm_remainder (SCM x
, SCM y
)
758 long yy
= SCM_INUM (y
);
760 scm_num_overflow (s_remainder
);
763 long z
= SCM_INUM (x
) % yy
;
764 return SCM_MAKINUM (z
);
767 else if (SCM_BIGP (y
))
769 if ((SCM_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
770 && (scm_i_bigcmp (abs_most_negative_fixnum
, y
) == 0))
771 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
772 return SCM_MAKINUM (0);
777 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
779 else if (SCM_BIGP (x
))
783 long yy
= SCM_INUM (y
);
785 scm_num_overflow (s_remainder
);
788 SCM result
= scm_i_mkbig ();
791 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
792 scm_remember_upto_here_1 (x
);
793 return scm_i_normbig (result
);
796 else if (SCM_BIGP (y
))
798 SCM result
= scm_i_mkbig ();
799 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
802 scm_remember_upto_here_2 (x
, y
);
803 return scm_i_normbig (result
);
806 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
809 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
813 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
814 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
816 * "(modulo 13 4) @result{} 1\n"
817 * "(modulo -13 4) @result{} 3\n"
821 scm_modulo (SCM x
, SCM y
)
825 long xx
= SCM_INUM (x
);
828 long yy
= SCM_INUM (y
);
830 scm_num_overflow (s_modulo
);
833 /* FIXME: I think this may be a bug on some arches -- results
834 of % with negative second arg are undefined... */
852 return SCM_MAKINUM (result
);
855 else if (SCM_BIGP (y
))
857 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
860 scm_num_overflow (s_modulo
);
868 SCM pos_y
= scm_i_clonebig (y
, 0);
869 /* do this after the last scm_op */
870 mpz_init_set_si (z_x
, xx
);
871 result
= pos_y
; /* re-use this bignum */
872 mpz_mod (SCM_I_BIG_MPZ (result
),
874 SCM_I_BIG_MPZ (pos_y
));
875 scm_remember_upto_here_1 (pos_y
);
879 result
= scm_i_mkbig ();
880 /* do this after the last scm_op */
881 mpz_init_set_si (z_x
, xx
);
882 mpz_mod (SCM_I_BIG_MPZ (result
),
885 scm_remember_upto_here_1 (y
);
888 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
889 mpz_add (SCM_I_BIG_MPZ (result
),
891 SCM_I_BIG_MPZ (result
));
892 scm_remember_upto_here_1 (y
);
893 /* and do this before the next one */
895 return scm_i_normbig (result
);
899 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
901 else if (SCM_BIGP (x
))
905 long yy
= SCM_INUM (y
);
907 scm_num_overflow (s_modulo
);
910 SCM result
= scm_i_mkbig ();
911 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
913 (yy
< 0) ? - yy
: yy
);
914 scm_remember_upto_here_1 (x
);
915 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
916 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
917 SCM_I_BIG_MPZ (result
),
919 return scm_i_normbig (result
);
922 else if (SCM_BIGP (y
))
924 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
926 scm_num_overflow (s_modulo
);
929 SCM result
= scm_i_mkbig ();
930 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
931 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
932 mpz_mod (SCM_I_BIG_MPZ (result
),
934 SCM_I_BIG_MPZ (pos_y
));
936 scm_remember_upto_here_1 (x
);
937 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
938 mpz_add (SCM_I_BIG_MPZ (result
),
940 SCM_I_BIG_MPZ (result
));
941 scm_remember_upto_here_2 (y
, pos_y
);
942 return scm_i_normbig (result
);
946 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
949 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
952 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
953 /* "Return the greatest common divisor of all arguments.\n"
954 * "If called without arguments, 0 is returned."
957 scm_gcd (SCM x
, SCM y
)
960 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
966 long xx
= SCM_INUM (x
);
967 long yy
= SCM_INUM (y
);
968 long u
= xx
< 0 ? -xx
: xx
;
969 long v
= yy
< 0 ? -yy
: yy
;
979 /* Determine a common factor 2^k */
980 while (!(1 & (u
| v
)))
986 /* Now, any factor 2^n can be eliminated */
1006 return (SCM_POSFIXABLE (result
)
1007 ? SCM_MAKINUM (result
)
1008 : scm_i_long2big (result
));
1010 else if (SCM_BIGP (y
))
1012 SCM result
= scm_i_mkbig ();
1013 SCM mx
= scm_i_mkbig ();
1014 mpz_set_si (SCM_I_BIG_MPZ (mx
), SCM_INUM (x
));
1015 scm_remember_upto_here_1 (x
);
1016 mpz_gcd (SCM_I_BIG_MPZ (result
),
1019 scm_remember_upto_here_2 (mx
, y
);
1020 return scm_i_normbig (result
);
1023 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1025 else if (SCM_BIGP (x
))
1029 unsigned long result
;
1030 long yy
= SCM_INUM (y
);
1035 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1036 scm_remember_upto_here_1 (x
);
1037 return (SCM_POSFIXABLE (result
)
1038 ? SCM_MAKINUM (result
)
1039 : scm_ulong2num (result
));
1041 else if (SCM_BIGP (y
))
1043 SCM result
= scm_i_mkbig ();
1044 mpz_gcd (SCM_I_BIG_MPZ (result
),
1047 scm_remember_upto_here_2 (x
, y
);
1048 return scm_i_normbig (result
);
1051 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1054 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1057 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1058 /* "Return the least common multiple of the arguments.\n"
1059 * "If called without arguments, 1 is returned."
1062 scm_lcm (SCM n1
, SCM n2
)
1064 if (SCM_UNBNDP (n2
))
1066 if (SCM_UNBNDP (n1
))
1067 return SCM_MAKINUM (1L);
1068 n2
= SCM_MAKINUM (1L);
1071 SCM_GASSERT2 (SCM_INUMP (n1
) || SCM_BIGP (n1
),
1072 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1073 SCM_GASSERT2 (SCM_INUMP (n2
) || SCM_BIGP (n2
),
1074 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1080 SCM d
= scm_gcd (n1
, n2
);
1081 if (SCM_EQ_P (d
, SCM_INUM0
))
1084 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1088 /* inum n1, big n2 */
1091 SCM result
= scm_i_mkbig ();
1092 long nn1
= SCM_INUM (n1
);
1093 if (nn1
== 0) return SCM_INUM0
;
1094 if (nn1
< 0) nn1
= - nn1
;
1095 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1096 scm_remember_upto_here_1 (n2
);
1111 SCM result
= scm_i_mkbig ();
1112 mpz_lcm(SCM_I_BIG_MPZ (result
),
1114 SCM_I_BIG_MPZ (n2
));
1115 scm_remember_upto_here_2(n1
, n2
);
1116 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1122 #ifndef scm_long2num
1123 #define SCM_LOGOP_RETURN(x) scm_ulong2num(x)
1125 #define SCM_LOGOP_RETURN(x) SCM_MAKINUM(x)
1128 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1133 + + + x (map digit:logand X Y)
1134 + - + x (map digit:logand X (lognot (+ -1 Y)))
1135 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1136 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1141 + + + (map digit:logior X Y)
1142 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1143 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1144 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1149 + + + (map digit:logxor X Y)
1150 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1151 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1152 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1157 + + (any digit:logand X Y)
1158 + - (any digit:logand X (lognot (+ -1 Y)))
1159 - + (any digit:logand (lognot (+ -1 X)) Y)
1164 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1166 "Return the bitwise AND of the integer arguments.\n\n"
1168 "(logand) @result{} -1\n"
1169 "(logand 7) @result{} 7\n"
1170 "(logand #b111 #b011 #\b001) @result{} 1\n"
1172 #define FUNC_NAME s_scm_logand
1176 if (SCM_UNBNDP (n2
))
1178 if (SCM_UNBNDP (n1
))
1179 return SCM_MAKINUM (-1);
1180 else if (!SCM_NUMBERP (n1
))
1181 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1182 else if (SCM_NUMBERP (n1
))
1185 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1190 nn1
= SCM_INUM (n1
);
1193 long nn2
= SCM_INUM (n2
);
1194 return SCM_MAKINUM (nn1
& nn2
);
1196 else if SCM_BIGP (n2
)
1202 SCM result_z
= scm_i_mkbig ();
1204 mpz_init_set_si (nn1_z
, nn1
);
1205 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1206 scm_remember_upto_here_1 (n2
);
1208 return scm_i_normbig (result_z
);
1212 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1214 else if (SCM_BIGP (n1
))
1219 nn1
= SCM_INUM (n1
);
1222 else if (SCM_BIGP (n2
))
1224 SCM result_z
= scm_i_mkbig ();
1225 mpz_and (SCM_I_BIG_MPZ (result_z
),
1227 SCM_I_BIG_MPZ (n2
));
1228 scm_remember_upto_here_2 (n1
, n2
);
1229 return scm_i_normbig (result_z
);
1232 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1235 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1240 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1242 "Return the bitwise OR of the integer arguments.\n\n"
1244 "(logior) @result{} 0\n"
1245 "(logior 7) @result{} 7\n"
1246 "(logior #b000 #b001 #b011) @result{} 3\n"
1248 #define FUNC_NAME s_scm_logior
1252 if (SCM_UNBNDP (n2
))
1254 if (SCM_UNBNDP (n1
))
1256 else if (SCM_NUMBERP (n1
))
1259 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1264 nn1
= SCM_INUM (n1
);
1267 long nn2
= SCM_INUM (n2
);
1268 return SCM_MAKINUM (nn1
| nn2
);
1270 else if (SCM_BIGP (n2
))
1276 SCM result_z
= scm_i_mkbig ();
1278 mpz_init_set_si (nn1_z
, nn1
);
1279 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1280 scm_remember_upto_here_1 (n2
);
1286 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1288 else if (SCM_BIGP (n1
))
1293 nn1
= SCM_INUM (n1
);
1296 else if (SCM_BIGP (n2
))
1298 SCM result_z
= scm_i_mkbig ();
1299 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1301 SCM_I_BIG_MPZ (n2
));
1302 scm_remember_upto_here_2 (n1
, n2
);
1306 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1309 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1314 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1316 "Return the bitwise XOR of the integer arguments. A bit is\n"
1317 "set in the result if it is set in an odd number of arguments.\n"
1319 "(logxor) @result{} 0\n"
1320 "(logxor 7) @result{} 7\n"
1321 "(logxor #b000 #b001 #b011) @result{} 2\n"
1322 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1324 #define FUNC_NAME s_scm_logxor
1328 if (SCM_UNBNDP (n2
))
1330 if (SCM_UNBNDP (n1
))
1332 else if (SCM_NUMBERP (n1
))
1335 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1340 nn1
= SCM_INUM (n1
);
1343 long nn2
= SCM_INUM (n2
);
1344 return SCM_MAKINUM (nn1
^ nn2
);
1346 else if (SCM_BIGP (n2
))
1350 SCM result_z
= scm_i_mkbig ();
1352 mpz_init_set_si (nn1_z
, nn1
);
1353 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1354 scm_remember_upto_here_1 (n2
);
1356 return scm_i_normbig (result_z
);
1360 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1362 else if (SCM_BIGP (n1
))
1367 nn1
= SCM_INUM (n1
);
1370 else if (SCM_BIGP (n2
))
1372 SCM result_z
= scm_i_mkbig ();
1373 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1375 SCM_I_BIG_MPZ (n2
));
1376 scm_remember_upto_here_2 (n1
, n2
);
1377 return scm_i_normbig (result_z
);
1380 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1383 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1388 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1391 "(logtest j k) @equiv{} (not (zero? (logand j k)))\n\n"
1392 "(logtest #b0100 #b1011) @result{} #f\n"
1393 "(logtest #b0100 #b0111) @result{} #t\n"
1395 #define FUNC_NAME s_scm_logtest
1404 long nk
= SCM_INUM (k
);
1405 return SCM_BOOL (nj
& nk
);
1407 else if (SCM_BIGP (k
))
1415 mpz_init_set_si (nj_z
, nj
);
1416 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1417 scm_remember_upto_here_1 (k
);
1418 result
= SCM_BOOL (mpz_sgn (nj_z
) != 0);
1424 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1426 else if (SCM_BIGP (j
))
1434 else if (SCM_BIGP (k
))
1438 mpz_init (result_z
);
1442 scm_remember_upto_here_2 (j
, k
);
1443 result
= SCM_BOOL (mpz_sgn (result_z
) != 0);
1444 mpz_clear (result_z
);
1448 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1451 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1456 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1459 "(logbit? index j) @equiv{} (logtest (integer-expt 2 index) j)\n\n"
1460 "(logbit? 0 #b1101) @result{} #t\n"
1461 "(logbit? 1 #b1101) @result{} #f\n"
1462 "(logbit? 2 #b1101) @result{} #t\n"
1463 "(logbit? 3 #b1101) @result{} #t\n"
1464 "(logbit? 4 #b1101) @result{} #f\n"
1466 #define FUNC_NAME s_scm_logbit_p
1468 unsigned long int iindex
;
1470 SCM_VALIDATE_INUM_MIN (SCM_ARG1
, index
, 0);
1471 iindex
= (unsigned long int) SCM_INUM (index
);
1474 return SCM_BOOL ((1L << iindex
) & SCM_INUM (j
));
1475 else if (SCM_BIGP (j
))
1477 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1478 scm_remember_upto_here_1 (j
);
1479 return SCM_BOOL (val
);
1482 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1487 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1489 "Return the integer which is the ones-complement of the integer\n"
1493 "(number->string (lognot #b10000000) 2)\n"
1494 " @result{} \"-10000001\"\n"
1495 "(number->string (lognot #b0) 2)\n"
1496 " @result{} \"-1\"\n"
1498 #define FUNC_NAME s_scm_lognot
1500 if (SCM_INUMP (n
)) {
1501 /* No overflow here, just need to toggle all the bits making up the inum.
1502 Enhancement: No need to strip the tag and add it back, could just xor
1503 a block of 1 bits, if that worked with the various debug versions of
1505 return SCM_MAKINUM (~ SCM_INUM (n
));
1507 } else if (SCM_BIGP (n
)) {
1508 SCM result
= scm_i_mkbig ();
1509 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1510 scm_remember_upto_here_1 (n
);
1514 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1519 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1521 "Return @var{n} raised to the non-negative integer exponent\n"
1525 "(integer-expt 2 5)\n"
1527 "(integer-expt -3 3)\n"
1530 #define FUNC_NAME s_scm_integer_expt
1533 SCM z_i2
= SCM_BOOL_F
;
1535 SCM acc
= SCM_MAKINUM (1L);
1537 /* 0^0 == 1 according to R5RS */
1538 if (SCM_EQ_P (n
, SCM_INUM0
) || SCM_EQ_P (n
, acc
))
1539 return SCM_FALSEP (scm_zero_p(k
)) ? n
: acc
;
1540 else if (SCM_EQ_P (n
, SCM_MAKINUM (-1L)))
1541 return SCM_FALSEP (scm_even_p (k
)) ? n
: acc
;
1545 else if (SCM_BIGP (k
))
1547 z_i2
= scm_i_clonebig (k
, 1);
1548 mpz_init_set (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (k
));
1549 scm_remember_upto_here_1 (k
);
1552 else if (SCM_REALP (k
))
1554 double r
= SCM_REAL_VALUE (k
);
1556 SCM_WRONG_TYPE_ARG (2, k
);
1557 if ((r
> SCM_MOST_POSITIVE_FIXNUM
) || (r
< SCM_MOST_NEGATIVE_FIXNUM
))
1559 z_i2
= scm_i_mkbig ();
1560 mpz_init_set_d (SCM_I_BIG_MPZ (z_i2
), r
);
1569 SCM_WRONG_TYPE_ARG (2, k
);
1573 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1575 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1576 n
= scm_divide (n
, SCM_UNDEFINED
);
1580 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1582 mpz_clear (SCM_I_BIG_MPZ (z_i2
));
1585 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1587 mpz_clear (SCM_I_BIG_MPZ (z_i2
));
1588 return scm_product (acc
, n
);
1590 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1591 acc
= scm_product (acc
, n
);
1592 n
= scm_product (n
, n
);
1593 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1601 n
= scm_divide (n
, SCM_UNDEFINED
);
1608 return scm_product (acc
, n
);
1610 acc
= scm_product (acc
, n
);
1611 n
= scm_product (n
, n
);
1618 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1620 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1621 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1623 "This is effectively a multiplication by 2^@var{cnt}}, and when\n"
1624 "@var{cnt} is negative it's a division, rounded towards negative\n"
1625 "infinity. (Note that this is not the same rounding as\n"
1626 "@code{quotient} does.)\n"
1628 "With @var{n} viewed as an infinite precision twos complement,\n"
1629 "@code{ash} means a left shift introducing zero bits, or a right\n"
1630 "shift dropping bits.\n"
1633 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1634 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1636 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1637 "(ash -23 -2) @result{} -6\n"
1639 #define FUNC_NAME s_scm_ash
1643 SCM_VALIDATE_INUM (2, cnt
);
1645 bits_to_shift
= SCM_INUM (cnt
);
1647 if (bits_to_shift
< 0)
1649 /* Shift right by abs(cnt) bits. This is realized as a division
1650 by div:=2^abs(cnt). However, to guarantee the floor
1651 rounding, negative values require some special treatment.
1653 SCM div
= scm_integer_expt (SCM_MAKINUM (2),
1654 SCM_MAKINUM (-bits_to_shift
));
1656 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1657 if (SCM_FALSEP (scm_negative_p (n
)))
1658 return scm_quotient (n
, div
);
1660 return scm_sum (SCM_MAKINUM (-1L),
1661 scm_quotient (scm_sum (SCM_MAKINUM (1L), n
), div
));
1664 /* Shift left is done by multiplication with 2^CNT */
1665 return scm_product (n
, scm_integer_expt (SCM_MAKINUM (2), cnt
));
1670 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1671 (SCM n
, SCM start
, SCM end
),
1672 "Return the integer composed of the @var{start} (inclusive)\n"
1673 "through @var{end} (exclusive) bits of @var{n}. The\n"
1674 "@var{start}th bit becomes the 0-th bit in the result.\n"
1677 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1678 " @result{} \"1010\"\n"
1679 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1680 " @result{} \"10110\"\n"
1682 #define FUNC_NAME s_scm_bit_extract
1684 unsigned long int istart
, iend
;
1685 SCM_VALIDATE_INUM_MIN_COPY (2, start
,0, istart
);
1686 SCM_VALIDATE_INUM_MIN_COPY (3, end
, 0, iend
);
1687 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1691 long int in
= SCM_INUM (n
);
1692 unsigned long int bits
= iend
- istart
;
1694 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1696 /* Since we emulate two's complement encoded numbers, this
1697 * special case requires us to produce a result that has
1698 * more bits than can be stored in a fixnum. Thus, we fall
1699 * back to the more general algorithm that is used for
1705 if (istart
< SCM_I_FIXNUM_BIT
)
1708 if (bits
< SCM_I_FIXNUM_BIT
)
1709 return SCM_MAKINUM (in
& ((1L << bits
) - 1));
1710 else /* we know: in >= 0 */
1711 return SCM_MAKINUM (in
);
1714 return SCM_MAKINUM (-1L & ((1L << bits
) - 1));
1716 return SCM_MAKINUM (0);
1718 else if (SCM_BIGP (n
))
1722 SCM num1
= SCM_MAKINUM (1L);
1723 SCM num2
= SCM_MAKINUM (2L);
1724 SCM bits
= SCM_MAKINUM (iend
- istart
);
1725 SCM mask
= scm_difference (scm_integer_expt (num2
, bits
), num1
);
1726 return scm_logand (mask
, scm_ash (n
, SCM_MAKINUM (-istart
)));
1730 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1734 static const char scm_logtab
[] = {
1735 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1738 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1740 "Return the number of bits in integer @var{n}. If integer is\n"
1741 "positive, the 1-bits in its binary representation are counted.\n"
1742 "If negative, the 0-bits in its two's-complement binary\n"
1743 "representation are counted. If 0, 0 is returned.\n"
1746 "(logcount #b10101010)\n"
1753 #define FUNC_NAME s_scm_logcount
1757 unsigned long int c
= 0;
1758 long int nn
= SCM_INUM (n
);
1763 c
+= scm_logtab
[15 & nn
];
1766 return SCM_MAKINUM (c
);
1768 else if (SCM_BIGP (n
))
1770 unsigned long count
;
1771 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1772 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1774 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1775 scm_remember_upto_here_1 (n
);
1776 return SCM_MAKINUM (count
);
1779 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1784 static const char scm_ilentab
[] = {
1785 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1789 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1791 "Return the number of bits necessary to represent @var{n}.\n"
1794 "(integer-length #b10101010)\n"
1796 "(integer-length 0)\n"
1798 "(integer-length #b1111)\n"
1801 #define FUNC_NAME s_scm_integer_length
1805 unsigned long int c
= 0;
1807 long int nn
= SCM_INUM (n
);
1813 l
= scm_ilentab
[15 & nn
];
1816 return SCM_MAKINUM (c
- 4 + l
);
1818 else if (SCM_BIGP (n
))
1820 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1821 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1822 1 too big, so check for that and adjust. */
1823 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1824 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1825 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1826 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1828 scm_remember_upto_here_1 (n
);
1829 return SCM_MAKINUM (size
);
1832 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1836 /*** NUMBERS -> STRINGS ***/
1838 static const double fx
[] =
1839 { 0.0, 5e-1, 5e-2, 5e-3, 5e-4, 5e-5,
1840 5e-6, 5e-7, 5e-8, 5e-9, 5e-10,
1841 5e-11, 5e-12, 5e-13, 5e-14, 5e-15,
1842 5e-16, 5e-17, 5e-18, 5e-19, 5e-20};
1845 idbl2str (double f
, char *a
)
1847 int efmt
, dpt
, d
, i
, wp
= scm_dblprec
;
1853 #ifdef HAVE_COPYSIGN
1854 double sgn
= copysign (1.0, f
);
1860 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
1866 strcpy (a
, "-inf.0");
1868 strcpy (a
, "+inf.0");
1871 else if (xisnan (f
))
1873 strcpy (a
, "+nan.0");
1883 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
1884 make-uniform-vector, from causing infinite loops. */
1888 if (exp
-- < DBL_MIN_10_EXP
)
1899 if (exp
++ > DBL_MAX_10_EXP
)
1919 if (f
+ fx
[wp
] >= 10.0)
1926 dpt
= (exp
+ 9999) % 3;
1930 efmt
= (exp
< -3) || (exp
> wp
+ 2);
1955 if (f
+ fx
[wp
] >= 1.0)
1969 if ((dpt
> 4) && (exp
> 6))
1971 d
= (a
[0] == '-' ? 2 : 1);
1972 for (i
= ch
++; i
> d
; i
--)
1985 if (a
[ch
- 1] == '.')
1986 a
[ch
++] = '0'; /* trailing zero */
1995 for (i
= 10; i
<= exp
; i
*= 10);
1996 for (i
/= 10; i
; i
/= 10)
1998 a
[ch
++] = exp
/ i
+ '0';
2007 iflo2str (SCM flt
, char *str
)
2010 if (SCM_REALP (flt
))
2011 i
= idbl2str (SCM_REAL_VALUE (flt
), str
);
2014 i
= idbl2str (SCM_COMPLEX_REAL (flt
), str
);
2015 if (SCM_COMPLEX_IMAG (flt
) != 0.0)
2017 double imag
= SCM_COMPLEX_IMAG (flt
);
2018 /* Don't output a '+' for negative numbers or for Inf and
2019 NaN. They will provide their own sign. */
2020 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2022 i
+= idbl2str (imag
, &str
[i
]);
2029 /* convert a long to a string (unterminated). returns the number of
2030 characters in the result.
2032 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2034 scm_iint2str (long num
, int rad
, char *p
)
2038 unsigned long n
= (num
< 0) ? -num
: num
;
2040 for (n
/= rad
; n
> 0; n
/= rad
)
2057 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2062 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2064 "Return a string holding the external representation of the\n"
2065 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2066 "inexact, a radix of 10 will be used.")
2067 #define FUNC_NAME s_scm_number_to_string
2071 if (SCM_UNBNDP (radix
))
2075 SCM_VALIDATE_INUM (2, radix
);
2076 base
= SCM_INUM (radix
);
2077 /* FIXME: ask if range limit was OK, and if so, document */
2078 SCM_ASSERT_RANGE (2, radix
, (base
>= 2) && (base
<= 36));
2083 char num_buf
[SCM_INTBUFLEN
];
2084 size_t length
= scm_iint2str (SCM_INUM (n
), base
, num_buf
);
2085 return scm_mem2string (num_buf
, length
);
2087 else if (SCM_BIGP (n
))
2089 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2090 scm_remember_upto_here_1 (n
);
2091 return scm_take0str (str
);
2093 else if (SCM_FRACTIONP (n
))
2095 scm_i_fraction_reduce (n
);
2096 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2097 scm_mem2string ("/", 1),
2098 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2100 else if (SCM_INEXACTP (n
))
2102 char num_buf
[FLOBUFLEN
];
2103 return scm_mem2string (num_buf
, iflo2str (n
, num_buf
));
2106 SCM_WRONG_TYPE_ARG (1, n
);
2111 /* These print routines used to be stubbed here so that scm_repl.c
2112 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2115 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2117 char num_buf
[FLOBUFLEN
];
2118 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
), port
);
2123 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2126 char num_buf
[FLOBUFLEN
];
2127 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
), port
);
2132 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2135 scm_i_fraction_reduce (sexp
);
2136 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2137 scm_lfwrite (SCM_STRING_CHARS (str
), SCM_STRING_LENGTH (str
), port
);
2138 scm_remember_upto_here_1 (str
);
2143 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2145 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2146 scm_remember_upto_here_1 (exp
);
2147 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2151 /*** END nums->strs ***/
2154 /*** STRINGS -> NUMBERS ***/
2156 /* The following functions implement the conversion from strings to numbers.
2157 * The implementation somehow follows the grammar for numbers as it is given
2158 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2159 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2160 * points should be noted about the implementation:
2161 * * Each function keeps a local index variable 'idx' that points at the
2162 * current position within the parsed string. The global index is only
2163 * updated if the function could parse the corresponding syntactic unit
2165 * * Similarly, the functions keep track of indicators of inexactness ('#',
2166 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2167 * global exactness information is only updated after each part has been
2168 * successfully parsed.
2169 * * Sequences of digits are parsed into temporary variables holding fixnums.
2170 * Only if these fixnums would overflow, the result variables are updated
2171 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2172 * the temporary variables holding the fixnums are cleared, and the process
2173 * starts over again. If for example fixnums were able to store five decimal
2174 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2175 * and the result was computed as 12345 * 100000 + 67890. In other words,
2176 * only every five digits two bignum operations were performed.
2179 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2181 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2183 /* In non ASCII-style encodings the following macro might not work. */
2184 #define XDIGIT2UINT(d) (isdigit (d) ? (d) - '0' : tolower (d) - 'a' + 10)
2187 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2188 unsigned int radix
, enum t_exactness
*p_exactness
)
2190 unsigned int idx
= *p_idx
;
2191 unsigned int hash_seen
= 0;
2192 scm_t_bits shift
= 1;
2194 unsigned int digit_value
;
2204 digit_value
= XDIGIT2UINT (c
);
2205 if (digit_value
>= radix
)
2209 result
= SCM_MAKINUM (digit_value
);
2217 digit_value
= XDIGIT2UINT (c
);
2218 if (digit_value
>= radix
)
2230 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2232 result
= scm_product (result
, SCM_MAKINUM (shift
));
2234 result
= scm_sum (result
, SCM_MAKINUM (add
));
2241 shift
= shift
* radix
;
2242 add
= add
* radix
+ digit_value
;
2247 result
= scm_product (result
, SCM_MAKINUM (shift
));
2249 result
= scm_sum (result
, SCM_MAKINUM (add
));
2253 *p_exactness
= INEXACT
;
2259 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2260 * covers the parts of the rules that start at a potential point. The value
2261 * of the digits up to the point have been parsed by the caller and are given
2262 * in variable result. The content of *p_exactness indicates, whether a hash
2263 * has already been seen in the digits before the point.
2266 /* In non ASCII-style encodings the following macro might not work. */
2267 #define DIGIT2UINT(d) ((d) - '0')
2270 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2271 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2273 unsigned int idx
= *p_idx
;
2274 enum t_exactness x
= *p_exactness
;
2279 if (mem
[idx
] == '.')
2281 scm_t_bits shift
= 1;
2283 unsigned int digit_value
;
2284 SCM big_shift
= SCM_MAKINUM (1);
2295 digit_value
= DIGIT2UINT (c
);
2306 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2308 big_shift
= scm_product (big_shift
, SCM_MAKINUM (shift
));
2309 result
= scm_product (result
, SCM_MAKINUM (shift
));
2311 result
= scm_sum (result
, SCM_MAKINUM (add
));
2319 add
= add
* 10 + digit_value
;
2325 big_shift
= scm_product (big_shift
, SCM_MAKINUM (shift
));
2326 result
= scm_product (result
, SCM_MAKINUM (shift
));
2327 result
= scm_sum (result
, SCM_MAKINUM (add
));
2330 result
= scm_divide (result
, big_shift
);
2332 /* We've seen a decimal point, thus the value is implicitly inexact. */
2344 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2375 exponent
= DIGIT2UINT (c
);
2382 if (exponent
<= SCM_MAXEXP
)
2383 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2389 if (exponent
> SCM_MAXEXP
)
2391 size_t exp_len
= idx
- start
;
2392 SCM exp_string
= scm_mem2string (&mem
[start
], exp_len
);
2393 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2394 scm_out_of_range ("string->number", exp_num
);
2397 e
= scm_integer_expt (SCM_MAKINUM (10), SCM_MAKINUM (exponent
));
2399 result
= scm_product (result
, e
);
2401 result
= scm_divide2real (result
, e
);
2403 /* We've seen an exponent, thus the value is implicitly inexact. */
2421 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2424 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2425 unsigned int radix
, enum t_exactness
*p_exactness
)
2427 unsigned int idx
= *p_idx
;
2433 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2439 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2441 enum t_exactness x
= EXACT
;
2443 /* Cobble up the fractional part. We might want to set the
2444 NaN's mantissa from it. */
2446 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2451 if (mem
[idx
] == '.')
2455 else if (idx
+ 1 == len
)
2457 else if (!isdigit (mem
[idx
+ 1]))
2460 result
= mem2decimal_from_point (SCM_MAKINUM (0), mem
, len
,
2461 p_idx
, p_exactness
);
2465 enum t_exactness x
= EXACT
;
2468 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2469 if (SCM_FALSEP (uinteger
))
2474 else if (mem
[idx
] == '/')
2480 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2481 if (SCM_FALSEP (divisor
))
2484 /* both are int/big here, I assume */
2485 result
= scm_make_ratio (uinteger
, divisor
);
2487 else if (radix
== 10)
2489 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2490 if (SCM_FALSEP (result
))
2501 /* When returning an inexact zero, make sure it is represented as a
2502 floating point value so that we can change its sign.
2504 if (SCM_EQ_P (result
, SCM_MAKINUM(0)) && *p_exactness
== INEXACT
)
2505 result
= scm_make_real (0.0);
2511 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2514 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2515 unsigned int radix
, enum t_exactness
*p_exactness
)
2539 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2540 if (SCM_FALSEP (ureal
))
2542 /* input must be either +i or -i */
2547 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2553 return scm_make_rectangular (SCM_MAKINUM (0), SCM_MAKINUM (sign
));
2560 if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2561 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2570 /* either +<ureal>i or -<ureal>i */
2577 return scm_make_rectangular (SCM_MAKINUM (0), ureal
);
2580 /* polar input: <real>@<real>. */
2605 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2606 if (SCM_FALSEP (angle
))
2611 if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2612 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2614 result
= scm_make_polar (ureal
, angle
);
2619 /* expecting input matching <real>[+-]<ureal>?i */
2626 int sign
= (c
== '+') ? 1 : -1;
2627 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2629 if (SCM_FALSEP (imag
))
2630 imag
= SCM_MAKINUM (sign
);
2631 else if (sign
== -1 && SCM_FALSEP (scm_nan_p (ureal
)))
2632 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2636 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2643 return scm_make_rectangular (ureal
, imag
);
2652 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2654 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2657 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2659 unsigned int idx
= 0;
2660 unsigned int radix
= NO_RADIX
;
2661 enum t_exactness forced_x
= NO_EXACTNESS
;
2662 enum t_exactness implicit_x
= EXACT
;
2665 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2666 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2668 switch (mem
[idx
+ 1])
2671 if (radix
!= NO_RADIX
)
2676 if (radix
!= NO_RADIX
)
2681 if (forced_x
!= NO_EXACTNESS
)
2686 if (forced_x
!= NO_EXACTNESS
)
2691 if (radix
!= NO_RADIX
)
2696 if (radix
!= NO_RADIX
)
2706 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2707 if (radix
== NO_RADIX
)
2708 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2710 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2712 if (SCM_FALSEP (result
))
2718 if (SCM_INEXACTP (result
))
2719 return scm_inexact_to_exact (result
);
2723 if (SCM_INEXACTP (result
))
2726 return scm_exact_to_inexact (result
);
2729 if (implicit_x
== INEXACT
)
2731 if (SCM_INEXACTP (result
))
2734 return scm_exact_to_inexact (result
);
2742 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2743 (SCM string
, SCM radix
),
2744 "Return a number of the maximally precise representation\n"
2745 "expressed by the given @var{string}. @var{radix} must be an\n"
2746 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2747 "is a default radix that may be overridden by an explicit radix\n"
2748 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2749 "supplied, then the default radix is 10. If string is not a\n"
2750 "syntactically valid notation for a number, then\n"
2751 "@code{string->number} returns @code{#f}.")
2752 #define FUNC_NAME s_scm_string_to_number
2756 SCM_VALIDATE_STRING (1, string
);
2757 SCM_VALIDATE_INUM_MIN_DEF_COPY (2, radix
,2,10, base
);
2758 answer
= scm_i_mem2number (SCM_STRING_CHARS (string
),
2759 SCM_STRING_LENGTH (string
),
2761 return scm_return_first (answer
, string
);
2766 /*** END strs->nums ***/
2770 scm_make_real (double x
)
2772 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
2774 SCM_REAL_VALUE (z
) = x
;
2780 scm_make_complex (double x
, double y
)
2783 return scm_make_real (x
);
2787 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
2789 SCM_COMPLEX_REAL (z
) = x
;
2790 SCM_COMPLEX_IMAG (z
) = y
;
2797 scm_bigequal (SCM x
, SCM y
)
2799 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2800 scm_remember_upto_here_2 (x
, y
);
2801 return SCM_BOOL (0 == result
);
2805 scm_real_equalp (SCM x
, SCM y
)
2807 return SCM_BOOL (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
2811 scm_complex_equalp (SCM x
, SCM y
)
2813 return SCM_BOOL (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
2814 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
2818 scm_i_fraction_equalp (SCM x
, SCM y
)
2820 scm_i_fraction_reduce (x
);
2821 scm_i_fraction_reduce (y
);
2822 if (SCM_FALSEP (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
2823 SCM_FRACTION_NUMERATOR (y
)))
2824 || SCM_FALSEP (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
2825 SCM_FRACTION_DENOMINATOR (y
))))
2832 SCM_REGISTER_PROC (s_number_p
, "number?", 1, 0, 0, scm_number_p
);
2833 /* "Return @code{#t} if @var{x} is a number, @code{#f}\n"
2834 * "else. Note that the sets of complex, real, rational and\n"
2835 * "integer values form subsets of the set of numbers, i. e. the\n"
2836 * "predicate will be fulfilled for any number."
2838 SCM_DEFINE (scm_number_p
, "complex?", 1, 0, 0,
2840 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
2841 "otherwise. Note that the sets of real, rational and integer\n"
2842 "values form subsets of the set of complex numbers, i. e. the\n"
2843 "predicate will also be fulfilled if @var{x} is a real,\n"
2844 "rational or integer number.")
2845 #define FUNC_NAME s_scm_number_p
2847 return SCM_BOOL (SCM_NUMBERP (x
));
2852 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
2854 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
2855 "otherwise. Note that the set of integer values forms a subset of\n"
2856 "the set of real numbers, i. e. the predicate will also be\n"
2857 "fulfilled if @var{x} is an integer number.")
2858 #define FUNC_NAME s_scm_real_p
2860 /* we can't represent irrational numbers. */
2861 return scm_rational_p (x
);
2865 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
2867 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
2868 "otherwise. Note that the set of integer values forms a subset of\n"
2869 "the set of rational numbers, i. e. the predicate will also be\n"
2870 "fulfilled if @var{x} is an integer number.")
2871 #define FUNC_NAME s_scm_rational_p
2875 else if (SCM_IMP (x
))
2877 else if (SCM_BIGP (x
))
2879 else if (SCM_FRACTIONP (x
))
2881 else if (SCM_REALP (x
))
2882 /* due to their limited precision, all floating point numbers are
2883 rational as well. */
2891 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
2893 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
2895 #define FUNC_NAME s_scm_integer_p
2904 if (!SCM_INEXACTP (x
))
2906 if (SCM_COMPLEXP (x
))
2908 r
= SCM_REAL_VALUE (x
);
2916 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
2918 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
2920 #define FUNC_NAME s_scm_inexact_p
2922 if (SCM_INEXACTP (x
))
2924 if (SCM_NUMBERP (x
))
2926 SCM_WRONG_TYPE_ARG (1, x
);
2931 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
2932 /* "Return @code{#t} if all parameters are numerically equal." */
2934 scm_num_eq_p (SCM x
, SCM y
)
2938 long xx
= SCM_INUM (x
);
2941 long yy
= SCM_INUM (y
);
2942 return SCM_BOOL (xx
== yy
);
2944 else if (SCM_BIGP (y
))
2946 else if (SCM_REALP (y
))
2947 return SCM_BOOL ((double) xx
== SCM_REAL_VALUE (y
));
2948 else if (SCM_COMPLEXP (y
))
2949 return SCM_BOOL (((double) xx
== SCM_COMPLEX_REAL (y
))
2950 && (0.0 == SCM_COMPLEX_IMAG (y
)));
2951 else if (SCM_FRACTIONP (y
))
2954 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
2956 else if (SCM_BIGP (x
))
2960 else if (SCM_BIGP (y
))
2962 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2963 scm_remember_upto_here_2 (x
, y
);
2964 return SCM_BOOL (0 == cmp
);
2966 else if (SCM_REALP (y
))
2969 if (xisnan (SCM_REAL_VALUE (y
)))
2971 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
2972 scm_remember_upto_here_1 (x
);
2973 return SCM_BOOL (0 == cmp
);
2975 else if (SCM_COMPLEXP (y
))
2978 if (0.0 != SCM_COMPLEX_IMAG (y
))
2980 if (xisnan (SCM_COMPLEX_REAL (y
)))
2982 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
2983 scm_remember_upto_here_1 (x
);
2984 return SCM_BOOL (0 == cmp
);
2986 else if (SCM_FRACTIONP (y
))
2989 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
2991 else if (SCM_REALP (x
))
2994 return SCM_BOOL (SCM_REAL_VALUE (x
) == (double) SCM_INUM (y
));
2995 else if (SCM_BIGP (y
))
2998 if (xisnan (SCM_REAL_VALUE (x
)))
3000 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3001 scm_remember_upto_here_1 (y
);
3002 return SCM_BOOL (0 == cmp
);
3004 else if (SCM_REALP (y
))
3005 return SCM_BOOL (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3006 else if (SCM_COMPLEXP (y
))
3007 return SCM_BOOL ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3008 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3009 else if (SCM_FRACTIONP (y
))
3010 return SCM_BOOL (SCM_REAL_VALUE (x
) == scm_i_fraction2double (y
));
3012 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3014 else if (SCM_COMPLEXP (x
))
3017 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == (double) SCM_INUM (y
))
3018 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3019 else if (SCM_BIGP (y
))
3022 if (0.0 != SCM_COMPLEX_IMAG (x
))
3024 if (xisnan (SCM_COMPLEX_REAL (x
)))
3026 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3027 scm_remember_upto_here_1 (y
);
3028 return SCM_BOOL (0 == cmp
);
3030 else if (SCM_REALP (y
))
3031 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3032 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3033 else if (SCM_COMPLEXP (y
))
3034 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3035 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3036 else if (SCM_FRACTIONP (y
))
3037 return SCM_BOOL ((SCM_COMPLEX_REAL (x
) == scm_i_fraction2double (y
))
3038 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3040 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3042 else if (SCM_FRACTIONP (x
))
3046 else if (SCM_BIGP (y
))
3048 else if (SCM_REALP (y
))
3049 return SCM_BOOL (scm_i_fraction2double (x
) == SCM_REAL_VALUE (y
));
3050 else if (SCM_COMPLEXP (y
))
3051 return SCM_BOOL ((scm_i_fraction2double (x
) == SCM_COMPLEX_REAL (y
))
3052 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3053 else if (SCM_FRACTIONP (y
))
3054 return scm_i_fraction_equalp (x
, y
);
3056 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3059 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3063 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3064 /* "Return @code{#t} if the list of parameters is monotonically\n"
3068 scm_less_p (SCM x
, SCM y
)
3072 long xx
= SCM_INUM (x
);
3075 long yy
= SCM_INUM (y
);
3076 return SCM_BOOL (xx
< yy
);
3078 else if (SCM_BIGP (y
))
3080 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3081 scm_remember_upto_here_1 (y
);
3082 return SCM_BOOL (sgn
> 0);
3084 else if (SCM_REALP (y
))
3085 return SCM_BOOL ((double) xx
< SCM_REAL_VALUE (y
));
3086 else if (SCM_FRACTIONP (y
))
3087 return SCM_BOOL ((double) xx
< scm_i_fraction2double (y
));
3089 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3091 else if (SCM_BIGP (x
))
3095 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3096 scm_remember_upto_here_1 (x
);
3097 return SCM_BOOL (sgn
< 0);
3099 else if (SCM_BIGP (y
))
3101 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3102 scm_remember_upto_here_2 (x
, y
);
3103 return SCM_BOOL (cmp
< 0);
3105 else if (SCM_REALP (y
))
3108 if (xisnan (SCM_REAL_VALUE (y
)))
3110 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3111 scm_remember_upto_here_1 (x
);
3112 return SCM_BOOL (cmp
< 0);
3114 else if (SCM_FRACTIONP (y
))
3117 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), scm_i_fraction2double (y
));
3118 scm_remember_upto_here_1 (x
);
3119 return SCM_BOOL (cmp
< 0);
3122 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3124 else if (SCM_REALP (x
))
3127 return SCM_BOOL (SCM_REAL_VALUE (x
) < (double) SCM_INUM (y
));
3128 else if (SCM_BIGP (y
))
3131 if (xisnan (SCM_REAL_VALUE (x
)))
3133 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3134 scm_remember_upto_here_1 (y
);
3135 return SCM_BOOL (cmp
> 0);
3137 else if (SCM_REALP (y
))
3138 return SCM_BOOL (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3139 else if (SCM_FRACTIONP (y
))
3140 return SCM_BOOL (SCM_REAL_VALUE (x
) < scm_i_fraction2double (y
));
3142 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3144 else if (SCM_FRACTIONP (x
))
3147 return SCM_BOOL (scm_i_fraction2double (x
) < (double) SCM_INUM (y
));
3148 else if (SCM_BIGP (y
))
3151 if (xisnan (SCM_REAL_VALUE (x
)))
3153 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), scm_i_fraction2double (x
));
3154 scm_remember_upto_here_1 (y
);
3155 return SCM_BOOL (cmp
> 0);
3157 else if (SCM_REALP (y
))
3158 return SCM_BOOL (scm_i_fraction2double (x
) < SCM_REAL_VALUE (y
));
3159 else if (SCM_FRACTIONP (y
))
3160 return SCM_BOOL (scm_i_fraction2double (x
) < scm_i_fraction2double (y
));
3162 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3165 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3169 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3170 /* "Return @code{#t} if the list of parameters is monotonically\n"
3173 #define FUNC_NAME s_scm_gr_p
3175 scm_gr_p (SCM x
, SCM y
)
3177 if (!SCM_NUMBERP (x
))
3178 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3179 else if (!SCM_NUMBERP (y
))
3180 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3182 return scm_less_p (y
, x
);
3187 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3188 /* "Return @code{#t} if the list of parameters is monotonically\n"
3191 #define FUNC_NAME s_scm_leq_p
3193 scm_leq_p (SCM x
, SCM y
)
3195 if (!SCM_NUMBERP (x
))
3196 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3197 else if (!SCM_NUMBERP (y
))
3198 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3199 else if (SCM_NFALSEP (scm_nan_p (x
)) || SCM_NFALSEP (scm_nan_p (y
)))
3202 return SCM_BOOL_NOT (scm_less_p (y
, x
));
3207 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3208 /* "Return @code{#t} if the list of parameters is monotonically\n"
3211 #define FUNC_NAME s_scm_geq_p
3213 scm_geq_p (SCM x
, SCM y
)
3215 if (!SCM_NUMBERP (x
))
3216 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3217 else if (!SCM_NUMBERP (y
))
3218 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3219 else if (SCM_NFALSEP (scm_nan_p (x
)) || SCM_NFALSEP (scm_nan_p (y
)))
3222 return SCM_BOOL_NOT (scm_less_p (x
, y
));
3227 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3228 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3235 return SCM_BOOL (SCM_EQ_P (z
, SCM_INUM0
));
3236 else if (SCM_BIGP (z
))
3238 else if (SCM_REALP (z
))
3239 return SCM_BOOL (SCM_REAL_VALUE (z
) == 0.0);
3240 else if (SCM_COMPLEXP (z
))
3241 return SCM_BOOL (SCM_COMPLEX_REAL (z
) == 0.0
3242 && SCM_COMPLEX_IMAG (z
) == 0.0);
3243 else if (SCM_FRACTIONP (z
))
3246 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3250 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3251 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3255 scm_positive_p (SCM x
)
3258 return SCM_BOOL (SCM_INUM (x
) > 0);
3259 else if (SCM_BIGP (x
))
3261 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3262 scm_remember_upto_here_1 (x
);
3263 return SCM_BOOL (sgn
> 0);
3265 else if (SCM_REALP (x
))
3266 return SCM_BOOL(SCM_REAL_VALUE (x
) > 0.0);
3267 else if (SCM_FRACTIONP (x
))
3268 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3270 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3274 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3275 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3279 scm_negative_p (SCM x
)
3282 return SCM_BOOL (SCM_INUM (x
) < 0);
3283 else if (SCM_BIGP (x
))
3285 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3286 scm_remember_upto_here_1 (x
);
3287 return SCM_BOOL (sgn
< 0);
3289 else if (SCM_REALP (x
))
3290 return SCM_BOOL(SCM_REAL_VALUE (x
) < 0.0);
3291 else if (SCM_FRACTIONP (x
))
3292 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3294 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3298 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3299 /* "Return the maximum of all parameter values."
3302 scm_max (SCM x
, SCM y
)
3307 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3308 else if (SCM_NUMBERP (x
))
3311 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3316 long xx
= SCM_INUM (x
);
3319 long yy
= SCM_INUM (y
);
3320 return (xx
< yy
) ? y
: x
;
3322 else if (SCM_BIGP (y
))
3324 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3325 scm_remember_upto_here_1 (y
);
3326 return (sgn
< 0) ? x
: y
;
3328 else if (SCM_REALP (y
))
3331 /* if y==NaN then ">" is false and we return NaN */
3332 return (z
> SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3334 else if (SCM_FRACTIONP (y
))
3337 return (z
> scm_i_fraction2double (y
)) ? x
: y
;
3340 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3342 else if (SCM_BIGP (x
))
3346 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3347 scm_remember_upto_here_1 (x
);
3348 return (sgn
< 0) ? y
: x
;
3350 else if (SCM_BIGP (y
))
3352 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3353 scm_remember_upto_here_2 (x
, y
);
3354 return (cmp
> 0) ? x
: y
;
3356 else if (SCM_REALP (y
))
3358 double yy
= SCM_REAL_VALUE (y
);
3362 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3363 scm_remember_upto_here_1 (x
);
3364 return (cmp
> 0) ? x
: y
;
3366 else if (SCM_FRACTIONP (y
))
3368 double yy
= scm_i_fraction2double (y
);
3370 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3371 scm_remember_upto_here_1 (x
);
3372 return (cmp
> 0) ? x
: y
;
3375 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3377 else if (SCM_REALP (x
))
3381 double z
= SCM_INUM (y
);
3382 /* if x==NaN then "<" is false and we return NaN */
3383 return (SCM_REAL_VALUE (x
) < z
) ? scm_make_real (z
) : x
;
3385 else if (SCM_BIGP (y
))
3387 double xx
= SCM_REAL_VALUE (x
);
3391 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3392 scm_remember_upto_here_1 (y
);
3393 return (cmp
< 0) ? x
: y
;
3395 else if (SCM_REALP (y
))
3397 /* if x==NaN then our explicit check means we return NaN
3398 if y==NaN then ">" is false and we return NaN
3399 calling isnan is unavoidable, since it's the only way to know
3400 which of x or y causes any compares to be false */
3401 double xx
= SCM_REAL_VALUE (x
);
3402 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3404 else if (SCM_FRACTIONP (y
))
3406 double yy
= scm_i_fraction2double (y
);
3407 double xx
= SCM_REAL_VALUE (x
);
3408 return (xx
< yy
) ? scm_make_real (yy
) : x
;
3411 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3413 else if (SCM_FRACTIONP (x
))
3417 double z
= SCM_INUM (y
);
3418 return (scm_i_fraction2double (x
) < z
) ? y
: x
;
3420 else if (SCM_BIGP (y
))
3422 double xx
= scm_i_fraction2double (x
);
3424 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3425 scm_remember_upto_here_1 (y
);
3426 return (cmp
< 0) ? x
: y
;
3428 else if (SCM_REALP (y
))
3430 double xx
= scm_i_fraction2double (x
);
3431 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_make_real (xx
);
3433 else if (SCM_FRACTIONP (y
))
3435 double yy
= scm_i_fraction2double (y
);
3436 double xx
= scm_i_fraction2double (x
);
3437 return (xx
< yy
) ? y
: x
;
3440 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3443 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3447 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3448 /* "Return the minium of all parameter values."
3451 scm_min (SCM x
, SCM y
)
3456 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3457 else if (SCM_NUMBERP (x
))
3460 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3465 long xx
= SCM_INUM (x
);
3468 long yy
= SCM_INUM (y
);
3469 return (xx
< yy
) ? x
: y
;
3471 else if (SCM_BIGP (y
))
3473 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3474 scm_remember_upto_here_1 (y
);
3475 return (sgn
< 0) ? y
: x
;
3477 else if (SCM_REALP (y
))
3480 /* if y==NaN then "<" is false and we return NaN */
3481 return (z
< SCM_REAL_VALUE (y
)) ? scm_make_real (z
) : y
;
3483 else if (SCM_FRACTIONP (y
))
3486 return (z
< scm_i_fraction2double (y
)) ? x
: y
;
3489 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3491 else if (SCM_BIGP (x
))
3495 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3496 scm_remember_upto_here_1 (x
);
3497 return (sgn
< 0) ? x
: y
;
3499 else if (SCM_BIGP (y
))
3501 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3502 scm_remember_upto_here_2 (x
, y
);
3503 return (cmp
> 0) ? y
: x
;
3505 else if (SCM_REALP (y
))
3507 double yy
= SCM_REAL_VALUE (y
);
3511 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3512 scm_remember_upto_here_1 (x
);
3513 return (cmp
> 0) ? y
: x
;
3515 else if (SCM_FRACTIONP (y
))
3517 double yy
= scm_i_fraction2double (y
);
3519 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), yy
);
3520 scm_remember_upto_here_1 (x
);
3521 return (cmp
> 0) ? y
: x
;
3524 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3526 else if (SCM_REALP (x
))
3530 double z
= SCM_INUM (y
);
3531 /* if x==NaN then "<" is false and we return NaN */
3532 return (z
< SCM_REAL_VALUE (x
)) ? scm_make_real (z
) : x
;
3534 else if (SCM_BIGP (y
))
3536 double xx
= SCM_REAL_VALUE (x
);
3540 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3541 scm_remember_upto_here_1 (y
);
3542 return (cmp
< 0) ? y
: x
;
3544 else if (SCM_REALP (y
))
3546 /* if x==NaN then our explicit check means we return NaN
3547 if y==NaN then "<" is false and we return NaN
3548 calling isnan is unavoidable, since it's the only way to know
3549 which of x or y causes any compares to be false */
3550 double xx
= SCM_REAL_VALUE (x
);
3551 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3553 else if (SCM_FRACTIONP (y
))
3555 double yy
= scm_i_fraction2double (y
);
3556 double xx
= SCM_REAL_VALUE (x
);
3557 return (yy
< xx
) ? scm_make_real (yy
) : x
;
3560 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3562 else if (SCM_FRACTIONP (x
))
3566 double z
= SCM_INUM (y
);
3567 return (scm_i_fraction2double (x
) < z
) ? x
: y
;
3569 else if (SCM_BIGP (y
))
3571 double xx
= scm_i_fraction2double (x
);
3573 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), xx
);
3574 scm_remember_upto_here_1 (y
);
3575 return (cmp
< 0) ? y
: x
;
3577 else if (SCM_REALP (y
))
3579 double xx
= scm_i_fraction2double (x
);
3580 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_make_real (xx
);
3582 else if (SCM_FRACTIONP (y
))
3584 double yy
= scm_i_fraction2double (y
);
3585 double xx
= scm_i_fraction2double (x
);
3586 return (xx
< yy
) ? x
: y
;
3589 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3592 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3596 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3597 /* "Return the sum of all parameter values. Return 0 if called without\n"
3601 scm_sum (SCM x
, SCM y
)
3605 if (SCM_NUMBERP (x
)) return x
;
3606 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3607 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3614 long xx
= SCM_INUM (x
);
3615 long yy
= SCM_INUM (y
);
3616 long int z
= xx
+ yy
;
3617 return SCM_FIXABLE (z
) ? SCM_MAKINUM (z
) : scm_i_long2big (z
);
3619 else if (SCM_BIGP (y
))
3624 else if (SCM_REALP (y
))
3626 long int xx
= SCM_INUM (x
);
3627 return scm_make_real (xx
+ SCM_REAL_VALUE (y
));
3629 else if (SCM_COMPLEXP (y
))
3631 long int xx
= SCM_INUM (x
);
3632 return scm_make_complex (xx
+ SCM_COMPLEX_REAL (y
),
3633 SCM_COMPLEX_IMAG (y
));
3635 else if (SCM_FRACTIONP (y
))
3636 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3637 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3638 SCM_FRACTION_DENOMINATOR (y
));
3640 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3641 } else if (SCM_BIGP (x
))
3648 inum
= SCM_INUM (y
);
3651 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3654 SCM result
= scm_i_mkbig ();
3655 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3656 scm_remember_upto_here_1 (x
);
3657 /* we know the result will have to be a bignum */
3660 return scm_i_normbig (result
);
3664 SCM result
= scm_i_mkbig ();
3665 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3666 scm_remember_upto_here_1 (x
);
3667 /* we know the result will have to be a bignum */
3670 return scm_i_normbig (result
);
3673 else if (SCM_BIGP (y
))
3675 SCM result
= scm_i_mkbig ();
3676 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3677 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3678 mpz_add (SCM_I_BIG_MPZ (result
),
3681 scm_remember_upto_here_2 (x
, y
);
3682 /* we know the result will have to be a bignum */
3685 return scm_i_normbig (result
);
3687 else if (SCM_REALP (y
))
3689 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3690 scm_remember_upto_here_1 (x
);
3691 return scm_make_real (result
);
3693 else if (SCM_COMPLEXP (y
))
3695 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3696 + SCM_COMPLEX_REAL (y
));
3697 scm_remember_upto_here_1 (x
);
3698 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
3700 else if (SCM_FRACTIONP (y
))
3701 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3702 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3703 SCM_FRACTION_DENOMINATOR (y
));
3705 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3707 else if (SCM_REALP (x
))
3710 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_INUM (y
));
3711 else if (SCM_BIGP (y
))
3713 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3714 scm_remember_upto_here_1 (y
);
3715 return scm_make_real (result
);
3717 else if (SCM_REALP (y
))
3718 return scm_make_real (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3719 else if (SCM_COMPLEXP (y
))
3720 return scm_make_complex (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3721 SCM_COMPLEX_IMAG (y
));
3722 else if (SCM_FRACTIONP (y
))
3723 return scm_make_real (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3725 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3727 else if (SCM_COMPLEXP (x
))
3730 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_INUM (y
),
3731 SCM_COMPLEX_IMAG (x
));
3732 else if (SCM_BIGP (y
))
3734 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3735 + SCM_COMPLEX_REAL (x
));
3736 scm_remember_upto_here_1 (y
);
3737 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (x
));
3739 else if (SCM_REALP (y
))
3740 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3741 SCM_COMPLEX_IMAG (x
));
3742 else if (SCM_COMPLEXP (y
))
3743 return scm_make_complex (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3744 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3745 else if (SCM_FRACTIONP (y
))
3746 return scm_make_complex (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3747 SCM_COMPLEX_IMAG (x
));
3749 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3751 else if (SCM_FRACTIONP (x
))
3754 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3755 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3756 SCM_FRACTION_DENOMINATOR (x
));
3757 else if (SCM_BIGP (y
))
3758 return scm_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3759 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3760 SCM_FRACTION_DENOMINATOR (x
));
3761 else if (SCM_REALP (y
))
3762 return scm_make_real (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
3763 else if (SCM_COMPLEXP (y
))
3764 return scm_make_complex (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
3765 SCM_COMPLEX_IMAG (y
));
3766 else if (SCM_FRACTIONP (y
))
3767 /* a/b + c/d = (ad + bc) / bd */
3768 return scm_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
3769 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
3770 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
3772 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3775 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
3779 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
3780 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
3781 * the sum of all but the first argument are subtracted from the first
3783 #define FUNC_NAME s_difference
3785 scm_difference (SCM x
, SCM y
)
3790 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
3794 long xx
= -SCM_INUM (x
);
3795 if (SCM_FIXABLE (xx
))
3796 return SCM_MAKINUM (xx
);
3798 return scm_i_long2big (xx
);
3800 else if (SCM_BIGP (x
))
3801 /* FIXME: do we really need to normalize here? */
3802 return scm_i_normbig (scm_i_clonebig (x
, 0));
3803 else if (SCM_REALP (x
))
3804 return scm_make_real (-SCM_REAL_VALUE (x
));
3805 else if (SCM_COMPLEXP (x
))
3806 return scm_make_complex (-SCM_COMPLEX_REAL (x
),
3807 -SCM_COMPLEX_IMAG (x
));
3808 else if (SCM_FRACTIONP (x
))
3809 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
3810 SCM_FRACTION_DENOMINATOR (x
));
3812 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
3819 long int xx
= SCM_INUM (x
);
3820 long int yy
= SCM_INUM (y
);
3821 long int z
= xx
- yy
;
3822 if (SCM_FIXABLE (z
))
3823 return SCM_MAKINUM (z
);
3825 return scm_i_long2big (z
);
3827 else if (SCM_BIGP (y
))
3829 /* inum-x - big-y */
3830 long xx
= SCM_INUM (x
);
3833 return scm_i_clonebig (y
, 0);
3836 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3837 SCM result
= scm_i_mkbig ();
3840 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
3843 /* x - y == -(y + -x) */
3844 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
3845 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
3847 scm_remember_upto_here_1 (y
);
3849 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
3850 /* we know the result will have to be a bignum */
3853 return scm_i_normbig (result
);
3856 else if (SCM_REALP (y
))
3858 long int xx
= SCM_INUM (x
);
3859 return scm_make_real (xx
- SCM_REAL_VALUE (y
));
3861 else if (SCM_COMPLEXP (y
))
3863 long int xx
= SCM_INUM (x
);
3864 return scm_make_complex (xx
- SCM_COMPLEX_REAL (y
),
3865 - SCM_COMPLEX_IMAG (y
));
3867 else if (SCM_FRACTIONP (y
))
3868 /* a - b/c = (ac - b) / c */
3869 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
3870 SCM_FRACTION_NUMERATOR (y
)),
3871 SCM_FRACTION_DENOMINATOR (y
));
3873 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3875 else if (SCM_BIGP (x
))
3879 /* big-x - inum-y */
3880 long yy
= SCM_INUM (y
);
3881 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3883 scm_remember_upto_here_1 (x
);
3885 return SCM_FIXABLE (-yy
) ? SCM_MAKINUM (-yy
) : scm_long2num (-yy
);
3888 SCM result
= scm_i_mkbig ();
3891 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
3893 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
3894 scm_remember_upto_here_1 (x
);
3896 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
3897 /* we know the result will have to be a bignum */
3900 return scm_i_normbig (result
);
3903 else if (SCM_BIGP (y
))
3905 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3906 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3907 SCM result
= scm_i_mkbig ();
3908 mpz_sub (SCM_I_BIG_MPZ (result
),
3911 scm_remember_upto_here_2 (x
, y
);
3912 /* we know the result will have to be a bignum */
3913 if ((sgn_x
== 1) && (sgn_y
== -1))
3915 if ((sgn_x
== -1) && (sgn_y
== 1))
3917 return scm_i_normbig (result
);
3919 else if (SCM_REALP (y
))
3921 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
3922 scm_remember_upto_here_1 (x
);
3923 return scm_make_real (result
);
3925 else if (SCM_COMPLEXP (y
))
3927 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3928 - SCM_COMPLEX_REAL (y
));
3929 scm_remember_upto_here_1 (x
);
3930 return scm_make_complex (real_part
, - SCM_COMPLEX_IMAG (y
));
3932 else if (SCM_FRACTIONP (y
))
3933 return scm_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
3934 SCM_FRACTION_NUMERATOR (y
)),
3935 SCM_FRACTION_DENOMINATOR (y
));
3936 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3938 else if (SCM_REALP (x
))
3941 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_INUM (y
));
3942 else if (SCM_BIGP (y
))
3944 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
3945 scm_remember_upto_here_1 (x
);
3946 return scm_make_real (result
);
3948 else if (SCM_REALP (y
))
3949 return scm_make_real (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
3950 else if (SCM_COMPLEXP (y
))
3951 return scm_make_complex (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
3952 -SCM_COMPLEX_IMAG (y
));
3953 else if (SCM_FRACTIONP (y
))
3954 return scm_make_real (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
3956 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3958 else if (SCM_COMPLEXP (x
))
3961 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_INUM (y
),
3962 SCM_COMPLEX_IMAG (x
));
3963 else if (SCM_BIGP (y
))
3965 double real_part
= (SCM_COMPLEX_REAL (x
)
3966 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
3967 scm_remember_upto_here_1 (x
);
3968 return scm_make_complex (real_part
, SCM_COMPLEX_IMAG (y
));
3970 else if (SCM_REALP (y
))
3971 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
3972 SCM_COMPLEX_IMAG (x
));
3973 else if (SCM_COMPLEXP (y
))
3974 return scm_make_complex (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
3975 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
3976 else if (SCM_FRACTIONP (y
))
3977 return scm_make_complex (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
3978 SCM_COMPLEX_IMAG (x
));
3980 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
3982 else if (SCM_FRACTIONP (x
))
3985 /* a/b - c = (a - cb) / b */
3986 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
3987 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
3988 SCM_FRACTION_DENOMINATOR (x
));
3989 else if (SCM_BIGP (y
))
3990 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
3991 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
3992 SCM_FRACTION_DENOMINATOR (x
));
3993 else if (SCM_REALP (y
))
3994 return scm_make_real (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
3995 else if (SCM_COMPLEXP (y
))
3996 return scm_make_complex (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
3997 -SCM_COMPLEX_IMAG (y
));
3998 else if (SCM_FRACTIONP (y
))
3999 /* a/b - c/d = (ad - bc) / bd */
4000 return scm_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4001 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4002 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4004 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4007 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4012 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4013 /* "Return the product of all arguments. If called without arguments,\n"
4017 scm_product (SCM x
, SCM y
)
4022 return SCM_MAKINUM (1L);
4023 else if (SCM_NUMBERP (x
))
4026 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4038 case 0: return x
; break;
4039 case 1: return y
; break;
4044 long yy
= SCM_INUM (y
);
4046 SCM k
= SCM_MAKINUM (kk
);
4047 if ((kk
== SCM_INUM (k
)) && (kk
/ xx
== yy
))
4051 SCM result
= scm_i_long2big (xx
);
4052 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4053 return scm_i_normbig (result
);
4056 else if (SCM_BIGP (y
))
4058 SCM result
= scm_i_mkbig ();
4059 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4060 scm_remember_upto_here_1 (y
);
4063 else if (SCM_REALP (y
))
4064 return scm_make_real (xx
* SCM_REAL_VALUE (y
));
4065 else if (SCM_COMPLEXP (y
))
4066 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4067 xx
* SCM_COMPLEX_IMAG (y
));
4068 else if (SCM_FRACTIONP (y
))
4069 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4070 SCM_FRACTION_DENOMINATOR (y
));
4072 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4074 else if (SCM_BIGP (x
))
4081 else if (SCM_BIGP (y
))
4083 SCM result
= scm_i_mkbig ();
4084 mpz_mul (SCM_I_BIG_MPZ (result
),
4087 scm_remember_upto_here_2 (x
, y
);
4090 else if (SCM_REALP (y
))
4092 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4093 scm_remember_upto_here_1 (x
);
4094 return scm_make_real (result
);
4096 else if (SCM_COMPLEXP (y
))
4098 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4099 scm_remember_upto_here_1 (x
);
4100 return scm_make_complex (z
* SCM_COMPLEX_REAL (y
),
4101 z
* SCM_COMPLEX_IMAG (y
));
4103 else if (SCM_FRACTIONP (y
))
4104 return scm_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4105 SCM_FRACTION_DENOMINATOR (y
));
4107 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4109 else if (SCM_REALP (x
))
4112 return scm_make_real (SCM_INUM (y
) * SCM_REAL_VALUE (x
));
4113 else if (SCM_BIGP (y
))
4115 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4116 scm_remember_upto_here_1 (y
);
4117 return scm_make_real (result
);
4119 else if (SCM_REALP (y
))
4120 return scm_make_real (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4121 else if (SCM_COMPLEXP (y
))
4122 return scm_make_complex (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4123 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4124 else if (SCM_FRACTIONP (y
))
4125 return scm_make_real (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4127 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4129 else if (SCM_COMPLEXP (x
))
4132 return scm_make_complex (SCM_INUM (y
) * SCM_COMPLEX_REAL (x
),
4133 SCM_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4134 else if (SCM_BIGP (y
))
4136 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4137 scm_remember_upto_here_1 (y
);
4138 return scm_make_complex (z
* SCM_COMPLEX_REAL (x
),
4139 z
* SCM_COMPLEX_IMAG (x
));
4141 else if (SCM_REALP (y
))
4142 return scm_make_complex (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4143 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4144 else if (SCM_COMPLEXP (y
))
4146 return scm_make_complex (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4147 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4148 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4149 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4151 else if (SCM_FRACTIONP (y
))
4153 double yy
= scm_i_fraction2double (y
);
4154 return scm_make_complex (yy
* SCM_COMPLEX_REAL (x
),
4155 yy
* SCM_COMPLEX_IMAG (x
));
4158 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4160 else if (SCM_FRACTIONP (x
))
4163 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4164 SCM_FRACTION_DENOMINATOR (x
));
4165 else if (SCM_BIGP (y
))
4166 return scm_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4167 SCM_FRACTION_DENOMINATOR (x
));
4168 else if (SCM_REALP (y
))
4169 return scm_make_real (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4170 else if (SCM_COMPLEXP (y
))
4172 double xx
= scm_i_fraction2double (x
);
4173 return scm_make_complex (xx
* SCM_COMPLEX_REAL (y
),
4174 xx
* SCM_COMPLEX_IMAG (y
));
4176 else if (SCM_FRACTIONP (y
))
4177 /* a/b * c/d = ac / bd */
4178 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4179 SCM_FRACTION_NUMERATOR (y
)),
4180 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4181 SCM_FRACTION_DENOMINATOR (y
)));
4183 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4186 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4190 scm_num2dbl (SCM a
, const char *why
)
4191 #define FUNC_NAME why
4194 return (double) SCM_INUM (a
);
4195 else if (SCM_BIGP (a
))
4197 double result
= mpz_get_d (SCM_I_BIG_MPZ (a
));
4198 scm_remember_upto_here_1 (a
);
4201 else if (SCM_REALP (a
))
4202 return (SCM_REAL_VALUE (a
));
4203 else if (SCM_FRACTIONP (a
))
4204 return scm_i_fraction2double (a
);
4206 SCM_WRONG_TYPE_ARG (SCM_ARGn
, a
);
4210 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4211 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4212 #define ALLOW_DIVIDE_BY_ZERO
4213 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4216 /* The code below for complex division is adapted from the GNU
4217 libstdc++, which adapted it from f2c's libF77, and is subject to
4220 /****************************************************************
4221 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4223 Permission to use, copy, modify, and distribute this software
4224 and its documentation for any purpose and without fee is hereby
4225 granted, provided that the above copyright notice appear in all
4226 copies and that both that the copyright notice and this
4227 permission notice and warranty disclaimer appear in supporting
4228 documentation, and that the names of AT&T Bell Laboratories or
4229 Bellcore or any of their entities not be used in advertising or
4230 publicity pertaining to distribution of the software without
4231 specific, written prior permission.
4233 AT&T and Bellcore disclaim all warranties with regard to this
4234 software, including all implied warranties of merchantability
4235 and fitness. In no event shall AT&T or Bellcore be liable for
4236 any special, indirect or consequential damages or any damages
4237 whatsoever resulting from loss of use, data or profits, whether
4238 in an action of contract, negligence or other tortious action,
4239 arising out of or in connection with the use or performance of
4241 ****************************************************************/
4243 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4244 /* Divide the first argument by the product of the remaining
4245 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4247 #define FUNC_NAME s_divide
4249 scm_i_divide (SCM x
, SCM y
, int inexact
)
4256 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4257 else if (SCM_INUMP (x
))
4259 long xx
= SCM_INUM (x
);
4260 if (xx
== 1 || xx
== -1)
4262 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4264 scm_num_overflow (s_divide
);
4269 return scm_make_real (1.0 / (double) xx
);
4270 else return scm_make_ratio (SCM_MAKINUM(1), x
);
4273 else if (SCM_BIGP (x
))
4276 return scm_make_real (1.0 / scm_i_big2dbl (x
));
4277 else return scm_make_ratio (SCM_MAKINUM(1), x
);
4279 else if (SCM_REALP (x
))
4281 double xx
= SCM_REAL_VALUE (x
);
4282 #ifndef ALLOW_DIVIDE_BY_ZERO
4284 scm_num_overflow (s_divide
);
4287 return scm_make_real (1.0 / xx
);
4289 else if (SCM_COMPLEXP (x
))
4291 double r
= SCM_COMPLEX_REAL (x
);
4292 double i
= SCM_COMPLEX_IMAG (x
);
4296 double d
= i
* (1.0 + t
* t
);
4297 return scm_make_complex (t
/ d
, -1.0 / d
);
4302 double d
= r
* (1.0 + t
* t
);
4303 return scm_make_complex (1.0 / d
, -t
/ d
);
4306 else if (SCM_FRACTIONP (x
))
4307 return scm_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4308 SCM_FRACTION_NUMERATOR (x
));
4310 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4315 long xx
= SCM_INUM (x
);
4318 long yy
= SCM_INUM (y
);
4321 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4322 scm_num_overflow (s_divide
);
4324 return scm_make_real ((double) xx
/ (double) yy
);
4327 else if (xx
% yy
!= 0)
4330 return scm_make_real ((double) xx
/ (double) yy
);
4331 else return scm_make_ratio (x
, y
);
4336 if (SCM_FIXABLE (z
))
4337 return SCM_MAKINUM (z
);
4339 return scm_i_long2big (z
);
4342 else if (SCM_BIGP (y
))
4345 return scm_make_real ((double) xx
/ scm_i_big2dbl (y
));
4346 else return scm_make_ratio (x
, y
);
4348 else if (SCM_REALP (y
))
4350 double yy
= SCM_REAL_VALUE (y
);
4351 #ifndef ALLOW_DIVIDE_BY_ZERO
4353 scm_num_overflow (s_divide
);
4356 return scm_make_real ((double) xx
/ yy
);
4358 else if (SCM_COMPLEXP (y
))
4361 complex_div
: /* y _must_ be a complex number */
4363 double r
= SCM_COMPLEX_REAL (y
);
4364 double i
= SCM_COMPLEX_IMAG (y
);
4368 double d
= i
* (1.0 + t
* t
);
4369 return scm_make_complex ((a
* t
) / d
, -a
/ d
);
4374 double d
= r
* (1.0 + t
* t
);
4375 return scm_make_complex (a
/ d
, -(a
* t
) / d
);
4379 else if (SCM_FRACTIONP (y
))
4380 /* a / b/c = ac / b */
4381 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4382 SCM_FRACTION_NUMERATOR (y
));
4384 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4386 else if (SCM_BIGP (x
))
4390 long int yy
= SCM_INUM (y
);
4393 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4394 scm_num_overflow (s_divide
);
4396 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4397 scm_remember_upto_here_1 (x
);
4398 return (sgn
== 0) ? scm_nan () : scm_inf ();
4405 /* FIXME: HMM, what are the relative performance issues here?
4406 We need to test. Is it faster on average to test
4407 divisible_p, then perform whichever operation, or is it
4408 faster to perform the integer div opportunistically and
4409 switch to real if there's a remainder? For now we take the
4410 middle ground: test, then if divisible, use the faster div
4413 long abs_yy
= yy
< 0 ? -yy
: yy
;
4414 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4418 SCM result
= scm_i_mkbig ();
4419 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4420 scm_remember_upto_here_1 (x
);
4422 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4423 return scm_i_normbig (result
);
4428 return scm_make_real (scm_i_big2dbl (x
) / (double) yy
);
4429 else return scm_make_ratio (x
, y
);
4433 else if (SCM_BIGP (y
))
4435 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4438 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4439 scm_num_overflow (s_divide
);
4441 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4442 scm_remember_upto_here_1 (x
);
4443 return (sgn
== 0) ? scm_nan () : scm_inf ();
4449 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4453 SCM result
= scm_i_mkbig ();
4454 mpz_divexact (SCM_I_BIG_MPZ (result
),
4457 scm_remember_upto_here_2 (x
, y
);
4458 return scm_i_normbig (result
);
4464 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4465 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4466 scm_remember_upto_here_2 (x
, y
);
4467 return scm_make_real (dbx
/ dby
);
4469 else return scm_make_ratio (x
, y
);
4473 else if (SCM_REALP (y
))
4475 double yy
= SCM_REAL_VALUE (y
);
4476 #ifndef ALLOW_DIVIDE_BY_ZERO
4478 scm_num_overflow (s_divide
);
4481 return scm_make_real (scm_i_big2dbl (x
) / yy
);
4483 else if (SCM_COMPLEXP (y
))
4485 a
= scm_i_big2dbl (x
);
4488 else if (SCM_FRACTIONP (y
))
4489 return scm_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4490 SCM_FRACTION_NUMERATOR (y
));
4492 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4494 else if (SCM_REALP (x
))
4496 double rx
= SCM_REAL_VALUE (x
);
4499 long int yy
= SCM_INUM (y
);
4500 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4502 scm_num_overflow (s_divide
);
4505 return scm_make_real (rx
/ (double) yy
);
4507 else if (SCM_BIGP (y
))
4509 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4510 scm_remember_upto_here_1 (y
);
4511 return scm_make_real (rx
/ dby
);
4513 else if (SCM_REALP (y
))
4515 double yy
= SCM_REAL_VALUE (y
);
4516 #ifndef ALLOW_DIVIDE_BY_ZERO
4518 scm_num_overflow (s_divide
);
4521 return scm_make_real (rx
/ yy
);
4523 else if (SCM_COMPLEXP (y
))
4528 else if (SCM_FRACTIONP (y
))
4529 return scm_make_real (rx
/ scm_i_fraction2double (y
));
4531 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4533 else if (SCM_COMPLEXP (x
))
4535 double rx
= SCM_COMPLEX_REAL (x
);
4536 double ix
= SCM_COMPLEX_IMAG (x
);
4539 long int yy
= SCM_INUM (y
);
4540 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4542 scm_num_overflow (s_divide
);
4547 return scm_make_complex (rx
/ d
, ix
/ d
);
4550 else if (SCM_BIGP (y
))
4552 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4553 scm_remember_upto_here_1 (y
);
4554 return scm_make_complex (rx
/ dby
, ix
/ dby
);
4556 else if (SCM_REALP (y
))
4558 double yy
= SCM_REAL_VALUE (y
);
4559 #ifndef ALLOW_DIVIDE_BY_ZERO
4561 scm_num_overflow (s_divide
);
4564 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4566 else if (SCM_COMPLEXP (y
))
4568 double ry
= SCM_COMPLEX_REAL (y
);
4569 double iy
= SCM_COMPLEX_IMAG (y
);
4573 double d
= iy
* (1.0 + t
* t
);
4574 return scm_make_complex ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4579 double d
= ry
* (1.0 + t
* t
);
4580 return scm_make_complex ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4583 else if (SCM_FRACTIONP (y
))
4585 double yy
= scm_i_fraction2double (y
);
4586 return scm_make_complex (rx
/ yy
, ix
/ yy
);
4589 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4591 else if (SCM_FRACTIONP (x
))
4595 long int yy
= SCM_INUM (y
);
4596 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4598 scm_num_overflow (s_divide
);
4601 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4602 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4604 else if (SCM_BIGP (y
))
4606 return scm_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4607 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4609 else if (SCM_REALP (y
))
4611 double yy
= SCM_REAL_VALUE (y
);
4612 #ifndef ALLOW_DIVIDE_BY_ZERO
4614 scm_num_overflow (s_divide
);
4617 return scm_make_real (scm_i_fraction2double (x
) / yy
);
4619 else if (SCM_COMPLEXP (y
))
4621 a
= scm_i_fraction2double (x
);
4624 else if (SCM_FRACTIONP (y
))
4625 return scm_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4626 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4628 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4631 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4635 scm_divide (SCM x
, SCM y
)
4637 return scm_i_divide (x
, y
, 0);
4640 static SCM
scm_divide2real (SCM x
, SCM y
)
4642 return scm_i_divide (x
, y
, 1);
4648 scm_asinh (double x
)
4653 #define asinh scm_asinh
4654 return log (x
+ sqrt (x
* x
+ 1));
4657 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4658 /* "Return the inverse hyperbolic sine of @var{x}."
4663 scm_acosh (double x
)
4668 #define acosh scm_acosh
4669 return log (x
+ sqrt (x
* x
- 1));
4672 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4673 /* "Return the inverse hyperbolic cosine of @var{x}."
4678 scm_atanh (double x
)
4683 #define atanh scm_atanh
4684 return 0.5 * log ((1 + x
) / (1 - x
));
4687 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4688 /* "Return the inverse hyperbolic tangent of @var{x}."
4692 /* XXX - eventually, we should remove this definition of scm_round and
4693 rename scm_round_number to scm_round. Likewise for scm_truncate
4694 and scm_truncate_number.
4698 scm_truncate (double x
)
4703 #define trunc scm_truncate
4711 scm_round (double x
)
4713 double plus_half
= x
+ 0.5;
4714 double result
= floor (plus_half
);
4715 /* Adjust so that the scm_round is towards even. */
4716 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4721 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4723 "Round the number @var{x} towards zero.")
4724 #define FUNC_NAME s_scm_truncate_number
4726 if (SCM_FALSEP (scm_negative_p (x
)))
4727 return scm_floor (x
);
4729 return scm_ceiling (x
);
4733 static SCM exactly_one_half
;
4735 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4737 "Round the number @var{x} towards the nearest integer. "
4738 "When it is exactly halfway between two integers, "
4739 "round towards the even one.")
4740 #define FUNC_NAME s_scm_round_number
4742 SCM plus_half
= scm_sum (x
, exactly_one_half
);
4743 SCM result
= scm_floor (plus_half
);
4744 /* Adjust so that the scm_round is towards even. */
4745 if (!SCM_FALSEP (scm_num_eq_p (plus_half
, result
))
4746 && !SCM_FALSEP (scm_odd_p (result
)))
4747 return scm_difference (result
, SCM_MAKINUM (1));
4753 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
4755 "Round the number @var{x} towards minus infinity.")
4756 #define FUNC_NAME s_scm_floor
4758 if (SCM_INUMP (x
) || SCM_BIGP (x
))
4760 else if (SCM_REALP (x
))
4761 return scm_make_real (floor (SCM_REAL_VALUE (x
)));
4762 else if (SCM_FRACTIONP (x
))
4764 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
4765 SCM_FRACTION_DENOMINATOR (x
));
4766 if (SCM_FALSEP (scm_negative_p (x
)))
4768 /* For positive x, rounding towards zero is correct. */
4773 /* For negative x, we need to return q-1 unless x is an
4774 integer. But fractions are never integer, per our
4776 return scm_difference (q
, SCM_MAKINUM (1));
4780 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
4784 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
4786 "Round the number @var{x} towards infinity.")
4787 #define FUNC_NAME s_scm_ceiling
4789 if (SCM_INUMP (x
) || SCM_BIGP (x
))
4791 else if (SCM_REALP (x
))
4792 return scm_make_real (ceil (SCM_REAL_VALUE (x
)));
4793 else if (SCM_FRACTIONP (x
))
4795 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
4796 SCM_FRACTION_DENOMINATOR (x
));
4797 if (SCM_FALSEP (scm_positive_p (x
)))
4799 /* For negative x, rounding towards zero is correct. */
4804 /* For positive x, we need to return q+1 unless x is an
4805 integer. But fractions are never integer, per our
4807 return scm_sum (q
, SCM_MAKINUM (1));
4811 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
4815 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
4816 /* "Return the square root of the real number @var{x}."
4818 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
4819 /* "Return the absolute value of the real number @var{x}."
4821 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
4822 /* "Return the @var{x}th power of e."
4824 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
4825 /* "Return the natural logarithm of the real number @var{x}."
4827 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
4828 /* "Return the sine of the real number @var{x}."
4830 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
4831 /* "Return the cosine of the real number @var{x}."
4833 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
4834 /* "Return the tangent of the real number @var{x}."
4836 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
4837 /* "Return the arc sine of the real number @var{x}."
4839 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
4840 /* "Return the arc cosine of the real number @var{x}."
4842 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
4843 /* "Return the arc tangent of the real number @var{x}."
4845 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
4846 /* "Return the hyperbolic sine of the real number @var{x}."
4848 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
4849 /* "Return the hyperbolic cosine of the real number @var{x}."
4851 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
4852 /* "Return the hyperbolic tangent of the real number @var{x}."
4860 static void scm_two_doubles (SCM x
,
4862 const char *sstring
,
4866 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
4869 xy
->x
= SCM_INUM (x
);
4870 else if (SCM_BIGP (x
))
4871 xy
->x
= scm_i_big2dbl (x
);
4872 else if (SCM_REALP (x
))
4873 xy
->x
= SCM_REAL_VALUE (x
);
4874 else if (SCM_FRACTIONP (x
))
4875 xy
->x
= scm_i_fraction2double (x
);
4877 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
4880 xy
->y
= SCM_INUM (y
);
4881 else if (SCM_BIGP (y
))
4882 xy
->y
= scm_i_big2dbl (y
);
4883 else if (SCM_REALP (y
))
4884 xy
->y
= SCM_REAL_VALUE (y
);
4885 else if (SCM_FRACTIONP (y
))
4886 xy
->y
= scm_i_fraction2double (y
);
4888 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
4892 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
4894 "Return @var{x} raised to the power of @var{y}. This\n"
4895 "procedure does not accept complex arguments.")
4896 #define FUNC_NAME s_scm_sys_expt
4899 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4900 return scm_make_real (pow (xy
.x
, xy
.y
));
4905 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
4907 "Return the arc tangent of the two arguments @var{x} and\n"
4908 "@var{y}. This is similar to calculating the arc tangent of\n"
4909 "@var{x} / @var{y}, except that the signs of both arguments\n"
4910 "are used to determine the quadrant of the result. This\n"
4911 "procedure does not accept complex arguments.")
4912 #define FUNC_NAME s_scm_sys_atan2
4915 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4916 return scm_make_real (atan2 (xy
.x
, xy
.y
));
4921 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
4922 (SCM real
, SCM imaginary
),
4923 "Return a complex number constructed of the given @var{real} and\n"
4924 "@var{imaginary} parts.")
4925 #define FUNC_NAME s_scm_make_rectangular
4928 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
4929 return scm_make_complex (xy
.x
, xy
.y
);
4935 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
4937 "Return the complex number @var{x} * e^(i * @var{y}).")
4938 #define FUNC_NAME s_scm_make_polar
4942 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
4944 sincos (xy
.y
, &s
, &c
);
4949 return scm_make_complex (xy
.x
* c
, xy
.x
* s
);
4954 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
4955 /* "Return the real part of the number @var{z}."
4958 scm_real_part (SCM z
)
4962 else if (SCM_BIGP (z
))
4964 else if (SCM_REALP (z
))
4966 else if (SCM_COMPLEXP (z
))
4967 return scm_make_real (SCM_COMPLEX_REAL (z
));
4968 else if (SCM_FRACTIONP (z
))
4969 return scm_make_real (scm_i_fraction2double (z
));
4971 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
4975 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
4976 /* "Return the imaginary part of the number @var{z}."
4979 scm_imag_part (SCM z
)
4983 else if (SCM_BIGP (z
))
4985 else if (SCM_REALP (z
))
4987 else if (SCM_COMPLEXP (z
))
4988 return scm_make_real (SCM_COMPLEX_IMAG (z
));
4989 else if (SCM_FRACTIONP (z
))
4992 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
4995 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
4996 /* "Return the numerator of the number @var{z}."
4999 scm_numerator (SCM z
)
5003 else if (SCM_BIGP (z
))
5005 else if (SCM_FRACTIONP (z
))
5007 scm_i_fraction_reduce (z
);
5008 return SCM_FRACTION_NUMERATOR (z
);
5010 else if (SCM_REALP (z
))
5011 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5013 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5017 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5018 /* "Return the denominator of the number @var{z}."
5021 scm_denominator (SCM z
)
5024 return SCM_MAKINUM (1);
5025 else if (SCM_BIGP (z
))
5026 return SCM_MAKINUM (1);
5027 else if (SCM_FRACTIONP (z
))
5029 scm_i_fraction_reduce (z
);
5030 return SCM_FRACTION_DENOMINATOR (z
);
5032 else if (SCM_REALP (z
))
5033 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5035 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5038 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5039 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5040 * "@code{abs} for real arguments, but also allows complex numbers."
5043 scm_magnitude (SCM z
)
5047 long int zz
= SCM_INUM (z
);
5050 else if (SCM_POSFIXABLE (-zz
))
5051 return SCM_MAKINUM (-zz
);
5053 return scm_i_long2big (-zz
);
5055 else if (SCM_BIGP (z
))
5057 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5058 scm_remember_upto_here_1 (z
);
5060 return scm_i_clonebig (z
, 0);
5064 else if (SCM_REALP (z
))
5065 return scm_make_real (fabs (SCM_REAL_VALUE (z
)));
5066 else if (SCM_COMPLEXP (z
))
5067 return scm_make_real (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5068 else if (SCM_FRACTIONP (z
))
5070 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5072 return scm_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5073 SCM_FRACTION_DENOMINATOR (z
));
5076 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5080 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5081 /* "Return the angle of the complex number @var{z}."
5086 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5087 scm_flo0 to save allocating a new flonum with scm_make_real each time.
5088 But if atan2 follows the floating point rounding mode, then the value
5089 is not a constant. Maybe it'd be close enough though. */
5092 if (SCM_INUM (z
) >= 0)
5095 return scm_make_real (atan2 (0.0, -1.0));
5097 else if (SCM_BIGP (z
))
5099 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5100 scm_remember_upto_here_1 (z
);
5102 return scm_make_real (atan2 (0.0, -1.0));
5106 else if (SCM_REALP (z
))
5108 if (SCM_REAL_VALUE (z
) >= 0)
5111 return scm_make_real (atan2 (0.0, -1.0));
5113 else if (SCM_COMPLEXP (z
))
5114 return scm_make_real (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5115 else if (SCM_FRACTIONP (z
))
5117 if (SCM_FALSEP (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5119 else return scm_make_real (atan2 (0.0, -1.0));
5122 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5126 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5127 /* Convert the number @var{x} to its inexact representation.\n"
5130 scm_exact_to_inexact (SCM z
)
5133 return scm_make_real ((double) SCM_INUM (z
));
5134 else if (SCM_BIGP (z
))
5135 return scm_make_real (scm_i_big2dbl (z
));
5136 else if (SCM_FRACTIONP (z
))
5137 return scm_make_real (scm_i_fraction2double (z
));
5138 else if (SCM_INEXACTP (z
))
5141 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5145 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5147 "Return an exact number that is numerically closest to @var{z}.")
5148 #define FUNC_NAME s_scm_inexact_to_exact
5152 else if (SCM_BIGP (z
))
5154 else if (SCM_REALP (z
))
5156 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5157 SCM_OUT_OF_RANGE (1, z
);
5164 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5165 q
= scm_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5166 scm_i_mpz2num (mpq_denref (frac
)));
5168 /* When scm_make_ratio throws, we leak the memory allocated
5175 else if (SCM_FRACTIONP (z
))
5178 SCM_WRONG_TYPE_ARG (1, z
);
5182 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5184 "Return an exact number that is within @var{err} of @var{x}.")
5185 #define FUNC_NAME s_scm_rationalize
5189 else if (SCM_BIGP (x
))
5191 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5193 /* Use continued fractions to find closest ratio. All
5194 arithmetic is done with exact numbers.
5197 SCM ex
= scm_inexact_to_exact (x
);
5198 SCM int_part
= scm_floor (ex
);
5199 SCM tt
= SCM_MAKINUM (1);
5200 SCM a1
= SCM_MAKINUM (0), a2
= SCM_MAKINUM (1), a
= SCM_MAKINUM (0);
5201 SCM b1
= SCM_MAKINUM (1), b2
= SCM_MAKINUM (0), b
= SCM_MAKINUM (0);
5205 if (!SCM_FALSEP (scm_num_eq_p (ex
, int_part
)))
5208 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5209 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5211 /* We stop after a million iterations just to be absolutely sure
5212 that we don't go into an infinite loop. The process normally
5213 converges after less than a dozen iterations.
5216 err
= scm_abs (err
);
5217 while (++i
< 1000000)
5219 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5220 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5221 if (SCM_FALSEP (scm_zero_p (b
)) && /* b != 0 */
5223 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5224 err
))) /* abs(x-a/b) <= err */
5226 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5227 if (SCM_FALSEP (scm_exact_p (x
))
5228 || SCM_FALSEP (scm_exact_p (err
)))
5229 return scm_exact_to_inexact (res
);
5233 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5235 tt
= scm_floor (rx
); /* tt = floor (rx) */
5241 scm_num_overflow (s_scm_rationalize
);
5244 SCM_WRONG_TYPE_ARG (1, x
);
5248 /* if you need to change this, change test-num2integral.c as well */
5249 #if SCM_SIZEOF_LONG_LONG != 0
5251 # define ULLONG_MAX ((unsigned long long) (-1))
5252 # define LLONG_MAX ((long long) (ULLONG_MAX >> 1))
5253 # define LLONG_MIN (~LLONG_MAX)
5257 /* Parameters for creating integer conversion routines.
5259 Define the following preprocessor macros before including
5260 "libguile/num2integral.i.c":
5262 NUM2INTEGRAL - the name of the function for converting from a
5263 Scheme object to the integral type. This function will be
5264 defined when including "num2integral.i.c".
5266 INTEGRAL2NUM - the name of the function for converting from the
5267 integral type to a Scheme object. This function will be defined.
5269 INTEGRAL2BIG - the name of an internal function that createas a
5270 bignum from the integral type. This function will be defined.
5271 The name should start with "scm_i_".
5273 ITYPE - the name of the integral type.
5275 UNSIGNED - Define this to 1 when ITYPE is an unsigned type. Define
5278 UNSIGNED_ITYPE - the name of the the unsigned variant of the
5279 integral type. If you don't define this, it defaults to
5280 "unsigned ITYPE" for signed types and simply "ITYPE" for unsigned
5283 SIZEOF_ITYPE - an expression giving the size of the integral type
5284 in bytes. This expression must be computable by the
5285 preprocessor. (SIZEOF_FOO values are calculated by configure.in
5290 #define NUM2INTEGRAL scm_num2short
5291 #define INTEGRAL2NUM scm_short2num
5292 #define INTEGRAL2BIG scm_i_short2big
5295 #define SIZEOF_ITYPE SIZEOF_SHORT
5296 #include "libguile/num2integral.i.c"
5298 #define NUM2INTEGRAL scm_num2ushort
5299 #define INTEGRAL2NUM scm_ushort2num
5300 #define INTEGRAL2BIG scm_i_ushort2big
5302 #define ITYPE unsigned short
5303 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_SHORT
5304 #include "libguile/num2integral.i.c"
5306 #define NUM2INTEGRAL scm_num2int
5307 #define INTEGRAL2NUM scm_int2num
5308 #define INTEGRAL2BIG scm_i_int2big
5311 #define SIZEOF_ITYPE SIZEOF_INT
5312 #include "libguile/num2integral.i.c"
5314 #define NUM2INTEGRAL scm_num2uint
5315 #define INTEGRAL2NUM scm_uint2num
5316 #define INTEGRAL2BIG scm_i_uint2big
5318 #define ITYPE unsigned int
5319 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_INT
5320 #include "libguile/num2integral.i.c"
5322 #define NUM2INTEGRAL scm_num2long
5323 #define INTEGRAL2NUM scm_long2num
5324 #define INTEGRAL2BIG scm_i_long2big
5327 #define SIZEOF_ITYPE SIZEOF_LONG
5328 #include "libguile/num2integral.i.c"
5330 #define NUM2INTEGRAL scm_num2ulong
5331 #define INTEGRAL2NUM scm_ulong2num
5332 #define INTEGRAL2BIG scm_i_ulong2big
5334 #define ITYPE unsigned long
5335 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG
5336 #include "libguile/num2integral.i.c"
5338 #define NUM2INTEGRAL scm_num2ptrdiff
5339 #define INTEGRAL2NUM scm_ptrdiff2num
5340 #define INTEGRAL2BIG scm_i_ptrdiff2big
5342 #define ITYPE scm_t_ptrdiff
5343 #define UNSIGNED_ITYPE size_t
5344 #define SIZEOF_ITYPE SCM_SIZEOF_SCM_T_PTRDIFF
5345 #include "libguile/num2integral.i.c"
5347 #define NUM2INTEGRAL scm_num2size
5348 #define INTEGRAL2NUM scm_size2num
5349 #define INTEGRAL2BIG scm_i_size2big
5351 #define ITYPE size_t
5352 #define SIZEOF_ITYPE SIZEOF_SIZE_T
5353 #include "libguile/num2integral.i.c"
5355 #if SCM_SIZEOF_LONG_LONG != 0
5357 #ifndef ULONG_LONG_MAX
5358 #define ULONG_LONG_MAX (~0ULL)
5361 #define NUM2INTEGRAL scm_num2long_long
5362 #define INTEGRAL2NUM scm_long_long2num
5363 #define INTEGRAL2BIG scm_i_long_long2big
5365 #define ITYPE long long
5366 #define SIZEOF_ITYPE SIZEOF_LONG_LONG
5367 #include "libguile/num2integral.i.c"
5369 #define NUM2INTEGRAL scm_num2ulong_long
5370 #define INTEGRAL2NUM scm_ulong_long2num
5371 #define INTEGRAL2BIG scm_i_ulong_long2big
5373 #define ITYPE unsigned long long
5374 #define SIZEOF_ITYPE SIZEOF_UNSIGNED_LONG_LONG
5375 #include "libguile/num2integral.i.c"
5377 #endif /* SCM_SIZEOF_LONG_LONG != 0 */
5379 #define NUM2FLOAT scm_num2float
5380 #define FLOAT2NUM scm_float2num
5382 #include "libguile/num2float.i.c"
5384 #define NUM2FLOAT scm_num2double
5385 #define FLOAT2NUM scm_double2num
5386 #define FTYPE double
5387 #include "libguile/num2float.i.c"
5392 #define SIZE_MAX ((size_t) (-1))
5395 #define PTRDIFF_MIN \
5396 ((scm_t_ptrdiff) ((scm_t_ptrdiff) 1 \
5397 << ((sizeof (scm_t_ptrdiff) * SCM_CHAR_BIT) - 1)))
5400 #define PTRDIFF_MAX (~ PTRDIFF_MIN)
5403 #define CHECK(type, v) \
5406 if ((v) != scm_num2##type (scm_##type##2num (v), 1, "check_sanity")) \
5426 CHECK (ptrdiff
, -1);
5428 CHECK (short, SHRT_MAX
);
5429 CHECK (short, SHRT_MIN
);
5430 CHECK (ushort
, USHRT_MAX
);
5431 CHECK (int, INT_MAX
);
5432 CHECK (int, INT_MIN
);
5433 CHECK (uint
, UINT_MAX
);
5434 CHECK (long, LONG_MAX
);
5435 CHECK (long, LONG_MIN
);
5436 CHECK (ulong
, ULONG_MAX
);
5437 CHECK (size
, SIZE_MAX
);
5438 CHECK (ptrdiff
, PTRDIFF_MAX
);
5439 CHECK (ptrdiff
, PTRDIFF_MIN
);
5441 #if SCM_SIZEOF_LONG_LONG != 0
5442 CHECK (long_long
, 0LL);
5443 CHECK (ulong_long
, 0ULL);
5444 CHECK (long_long
, -1LL);
5445 CHECK (long_long
, LLONG_MAX
);
5446 CHECK (long_long
, LLONG_MIN
);
5447 CHECK (ulong_long
, ULLONG_MAX
);
5454 scm_internal_catch (SCM_BOOL_T, check_body, &data, check_handler, &data); \
5455 if (!SCM_FALSEP (data)) abort();
5458 check_body (void *data
)
5460 SCM num
= *(SCM
*) data
;
5461 scm_num2ulong (num
, 1, NULL
);
5463 return SCM_UNSPECIFIED
;
5467 check_handler (void *data
, SCM tag
, SCM throw_args
)
5469 SCM
*num
= (SCM
*) data
;
5472 return SCM_UNSPECIFIED
;
5475 SCM_DEFINE (scm_sys_check_number_conversions
, "%check-number-conversions", 0, 0, 0,
5477 "Number conversion sanity checking.")
5478 #define FUNC_NAME s_scm_sys_check_number_conversions
5480 SCM data
= SCM_MAKINUM (-1);
5482 data
= scm_int2num (INT_MIN
);
5484 data
= scm_ulong2num (ULONG_MAX
);
5485 data
= scm_difference (SCM_INUM0
, data
);
5487 data
= scm_ulong2num (ULONG_MAX
);
5488 data
= scm_sum (SCM_MAKINUM (1), data
); data
= scm_difference (SCM_INUM0
, data
);
5490 data
= scm_int2num (-10000); data
= scm_product (data
, data
); data
= scm_product (data
, data
);
5493 return SCM_UNSPECIFIED
;
5502 abs_most_negative_fixnum
= scm_i_long2big (- SCM_MOST_NEGATIVE_FIXNUM
);
5503 scm_permanent_object (abs_most_negative_fixnum
);
5505 mpz_init_set_si (z_negative_one
, -1);
5507 /* It may be possible to tune the performance of some algorithms by using
5508 * the following constants to avoid the creation of bignums. Please, before
5509 * using these values, remember the two rules of program optimization:
5510 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5511 scm_c_define ("most-positive-fixnum",
5512 SCM_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5513 scm_c_define ("most-negative-fixnum",
5514 SCM_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5516 scm_add_feature ("complex");
5517 scm_add_feature ("inexact");
5518 scm_flo0
= scm_make_real (0.0);
5520 scm_dblprec
= (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5522 { /* determine floating point precision */
5524 double fsum
= 1.0 + f
;
5527 if (++scm_dblprec
> 20)
5535 scm_dblprec
= scm_dblprec
- 1;
5537 #endif /* DBL_DIG */
5543 exactly_one_half
= scm_permanent_object (scm_divide (SCM_MAKINUM (1),
5545 #include "libguile/numbers.x"