1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004 Free Software Foundation, Inc.
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 /* General assumptions:
24 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
25 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
26 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
27 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
28 * All objects satisfying SCM_FRACTIONP are never an integer.
33 - see if special casing bignums and reals in integer-exponent when
34 possible (to use mpz_pow and mpf_pow_ui) is faster.
36 - look in to better short-circuiting of common cases in
37 integer-expt and elsewhere.
39 - see if direct mpz operations can help in ash and elsewhere.
43 /* tell glibc (2.3) to give prototype for C99 trunc() */
55 #include "libguile/_scm.h"
56 #include "libguile/feature.h"
57 #include "libguile/ports.h"
58 #include "libguile/root.h"
59 #include "libguile/smob.h"
60 #include "libguile/strings.h"
62 #include "libguile/validate.h"
63 #include "libguile/numbers.h"
64 #include "libguile/deprecation.h"
66 #include "libguile/eq.h"
68 #include "libguile/discouraged.h"
73 Wonder if this might be faster for some of our code? A switch on
74 the numtag would jump directly to the right case, and the
75 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
77 #define SCM_I_NUMTAG_NOTNUM 0
78 #define SCM_I_NUMTAG_INUM 1
79 #define SCM_I_NUMTAG_BIG scm_tc16_big
80 #define SCM_I_NUMTAG_REAL scm_tc16_real
81 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
82 #define SCM_I_NUMTAG(x) \
83 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
84 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
85 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
86 : SCM_I_NUMTAG_NOTNUM)))
88 /* the macro above will not work as is with fractions */
91 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
93 /* FLOBUFLEN is the maximum number of characters neccessary for the
94 * printed or scm_string representation of an inexact number.
96 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
99 #if ! defined (HAVE_ISNAN)
104 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
107 #if ! defined (HAVE_ISINF)
112 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
119 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
120 an explicit check. In some future gmp (don't know what version number),
121 mpz_cmp_d is supposed to do this itself. */
123 #define xmpz_cmp_d(z, d) \
124 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
126 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
129 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
130 isinf. It does have finite and isnan though, hence the use of those.
131 fpclass would be a possibility on that system too. */
135 #if defined (HAVE_ISINF)
137 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
138 return (! (finite (x
) || isnan (x
)));
147 #if defined (HAVE_ISNAN)
156 static mpz_t z_negative_one
;
160 SCM_C_INLINE_KEYWORD SCM
163 /* Return a newly created bignum. */
164 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
165 mpz_init (SCM_I_BIG_MPZ (z
));
169 SCM_C_INLINE_KEYWORD SCM
170 scm_i_long2big (long x
)
172 /* Return a newly created bignum initialized to X. */
173 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
174 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
178 SCM_C_INLINE_KEYWORD SCM
179 scm_i_ulong2big (unsigned long x
)
181 /* Return a newly created bignum initialized to X. */
182 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
183 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
187 SCM_C_INLINE_KEYWORD
static SCM
188 scm_i_clonebig (SCM src_big
, int same_sign_p
)
190 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
191 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
192 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
194 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
198 SCM_C_INLINE_KEYWORD
int
199 scm_i_bigcmp (SCM x
, SCM y
)
201 /* Return neg if x < y, pos if x > y, and 0 if x == y */
202 /* presume we already know x and y are bignums */
203 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
204 scm_remember_upto_here_2 (x
, y
);
208 SCM_C_INLINE_KEYWORD SCM
209 scm_i_dbl2big (double d
)
211 /* results are only defined if d is an integer */
212 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
213 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
217 /* Convert a integer in double representation to a SCM number. */
219 SCM_C_INLINE_KEYWORD SCM
220 scm_i_dbl2num (double u
)
222 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
223 powers of 2, so there's no rounding when making "double" values
224 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
225 get rounded on a 64-bit machine, hence the "+1".
227 The use of floor() to force to an integer value ensures we get a
228 "numerically closest" value without depending on how a
229 double->long cast or how mpz_set_d will round. For reference,
230 double->long probably follows the hardware rounding mode,
231 mpz_set_d truncates towards zero. */
233 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
234 representable as a double? */
236 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
237 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
238 return SCM_I_MAKINUM ((long) u
);
240 return scm_i_dbl2big (u
);
243 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
244 with R5RS exact->inexact.
246 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
247 (ie. truncate towards zero), then adjust to get the closest double by
248 examining the next lower bit and adding 1 (to the absolute value) if
251 Bignums exactly half way between representable doubles are rounded to the
252 next higher absolute value (ie. away from zero). This seems like an
253 adequate interpretation of R5RS "numerically closest", and it's easier
254 and faster than a full "nearest-even" style.
256 The bit test must be done on the absolute value of the mpz_t, which means
257 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
258 negatives as twos complement.
260 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
261 following the hardware rounding mode, but applied to the absolute value
262 of the mpz_t operand. This is not what we want so we put the high
263 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
264 mpz_get_d is supposed to always truncate towards zero.
266 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
267 is a slowdown. It'd be faster to pick out the relevant high bits with
268 mpz_getlimbn if we could be bothered coding that, and if the new
269 truncating gmp doesn't come out. */
272 scm_i_big2dbl (SCM b
)
277 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
281 /* Current GMP, eg. 4.1.3, force truncation towards zero */
283 if (bits
> DBL_MANT_DIG
)
285 size_t shift
= bits
- DBL_MANT_DIG
;
286 mpz_init2 (tmp
, DBL_MANT_DIG
);
287 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
288 result
= ldexp (mpz_get_d (tmp
), shift
);
293 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
298 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
301 if (bits
> DBL_MANT_DIG
)
303 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
304 /* test bit number "pos" in absolute value */
305 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
306 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
308 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
312 scm_remember_upto_here_1 (b
);
316 SCM_C_INLINE_KEYWORD SCM
317 scm_i_normbig (SCM b
)
319 /* convert a big back to a fixnum if it'll fit */
320 /* presume b is a bignum */
321 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
323 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
324 if (SCM_FIXABLE (val
))
325 b
= SCM_I_MAKINUM (val
);
330 static SCM_C_INLINE_KEYWORD SCM
331 scm_i_mpz2num (mpz_t b
)
333 /* convert a mpz number to a SCM number. */
334 if (mpz_fits_slong_p (b
))
336 long val
= mpz_get_si (b
);
337 if (SCM_FIXABLE (val
))
338 return SCM_I_MAKINUM (val
);
342 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
343 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
348 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
349 static SCM
scm_divide2real (SCM x
, SCM y
);
352 scm_i_make_ratio (SCM numerator
, SCM denominator
)
353 #define FUNC_NAME "make-ratio"
355 /* First make sure the arguments are proper.
357 if (SCM_I_INUMP (denominator
))
359 if (scm_is_eq (denominator
, SCM_INUM0
))
360 scm_num_overflow ("make-ratio");
361 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
366 if (!(SCM_BIGP(denominator
)))
367 SCM_WRONG_TYPE_ARG (2, denominator
);
369 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
370 SCM_WRONG_TYPE_ARG (1, numerator
);
372 /* Then flip signs so that the denominator is positive.
374 if (scm_is_true (scm_negative_p (denominator
)))
376 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
377 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
380 /* Now consider for each of the four fixnum/bignum combinations
381 whether the rational number is really an integer.
383 if (SCM_I_INUMP (numerator
))
385 long x
= SCM_I_INUM (numerator
);
386 if (scm_is_eq (numerator
, SCM_INUM0
))
388 if (SCM_I_INUMP (denominator
))
391 y
= SCM_I_INUM (denominator
);
393 return SCM_I_MAKINUM(1);
395 return SCM_I_MAKINUM (x
/ y
);
399 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
400 of that value for the denominator, as a bignum. Apart from
401 that case, abs(bignum) > abs(inum) so inum/bignum is not an
403 if (x
== SCM_MOST_NEGATIVE_FIXNUM
404 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
405 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
406 return SCM_I_MAKINUM(-1);
409 else if (SCM_BIGP (numerator
))
411 if (SCM_I_INUMP (denominator
))
413 long yy
= SCM_I_INUM (denominator
);
414 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
415 return scm_divide (numerator
, denominator
);
419 if (scm_is_eq (numerator
, denominator
))
420 return SCM_I_MAKINUM(1);
421 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
422 SCM_I_BIG_MPZ (denominator
)))
423 return scm_divide(numerator
, denominator
);
427 /* No, it's a proper fraction.
429 return scm_double_cell (scm_tc16_fraction
,
430 SCM_UNPACK (numerator
),
431 SCM_UNPACK (denominator
), 0);
435 static void scm_i_fraction_reduce (SCM z
)
437 if (!(SCM_FRACTION_REDUCED (z
)))
440 divisor
= scm_gcd (SCM_FRACTION_NUMERATOR (z
), SCM_FRACTION_DENOMINATOR (z
));
441 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
444 SCM_FRACTION_SET_NUMERATOR (z
, scm_divide (SCM_FRACTION_NUMERATOR (z
), divisor
));
445 SCM_FRACTION_SET_DENOMINATOR (z
, scm_divide (SCM_FRACTION_DENOMINATOR (z
), divisor
));
447 SCM_FRACTION_REDUCED_SET (z
);
452 scm_i_fraction2double (SCM z
)
454 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
455 SCM_FRACTION_DENOMINATOR (z
)));
458 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
460 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
462 #define FUNC_NAME s_scm_exact_p
468 if (SCM_FRACTIONP (x
))
472 SCM_WRONG_TYPE_ARG (1, x
);
477 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
479 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
481 #define FUNC_NAME s_scm_odd_p
485 long val
= SCM_I_INUM (n
);
486 return scm_from_bool ((val
& 1L) != 0);
488 else if (SCM_BIGP (n
))
490 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
491 scm_remember_upto_here_1 (n
);
492 return scm_from_bool (odd_p
);
494 else if (scm_is_true (scm_inf_p (n
)))
496 else if (SCM_REALP (n
))
498 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
504 SCM_WRONG_TYPE_ARG (1, n
);
507 SCM_WRONG_TYPE_ARG (1, n
);
512 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
514 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
516 #define FUNC_NAME s_scm_even_p
520 long val
= SCM_I_INUM (n
);
521 return scm_from_bool ((val
& 1L) == 0);
523 else if (SCM_BIGP (n
))
525 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
526 scm_remember_upto_here_1 (n
);
527 return scm_from_bool (even_p
);
529 else if (scm_is_true (scm_inf_p (n
)))
531 else if (SCM_REALP (n
))
533 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
539 SCM_WRONG_TYPE_ARG (1, n
);
542 SCM_WRONG_TYPE_ARG (1, n
);
546 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
548 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
549 "or @samp{-inf.0}, @code{#f} otherwise.")
550 #define FUNC_NAME s_scm_inf_p
553 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
554 else if (SCM_COMPLEXP (x
))
555 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
556 || xisinf (SCM_COMPLEX_IMAG (x
)));
562 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
564 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
566 #define FUNC_NAME s_scm_nan_p
569 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
570 else if (SCM_COMPLEXP (n
))
571 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
572 || xisnan (SCM_COMPLEX_IMAG (n
)));
578 /* Guile's idea of infinity. */
579 static double guile_Inf
;
581 /* Guile's idea of not a number. */
582 static double guile_NaN
;
585 guile_ieee_init (void)
587 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
589 /* Some version of gcc on some old version of Linux used to crash when
590 trying to make Inf and NaN. */
593 /* C99 INFINITY, when available.
594 FIXME: The standard allows for INFINITY to be something that overflows
595 at compile time. We ought to have a configure test to check for that
596 before trying to use it. (But in practice we believe this is not a
597 problem on any system guile is likely to target.) */
598 guile_Inf
= INFINITY
;
601 extern unsigned int DINFINITY
[2];
602 guile_Inf
= (*(X_CAST(double *, DINFINITY
)));
609 if (guile_Inf
== tmp
)
617 #if defined (HAVE_ISNAN)
620 /* C99 NAN, when available */
624 extern unsigned int DQNAN
[2];
625 guile_NaN
= (*(X_CAST(double *, DQNAN
)));
627 guile_NaN
= guile_Inf
/ guile_Inf
;
633 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
636 #define FUNC_NAME s_scm_inf
638 static int initialized
= 0;
644 return scm_from_double (guile_Inf
);
648 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
651 #define FUNC_NAME s_scm_nan
653 static int initialized
= 0;
659 return scm_from_double (guile_NaN
);
664 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
666 "Return the absolute value of @var{x}.")
671 long int xx
= SCM_I_INUM (x
);
674 else if (SCM_POSFIXABLE (-xx
))
675 return SCM_I_MAKINUM (-xx
);
677 return scm_i_long2big (-xx
);
679 else if (SCM_BIGP (x
))
681 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
683 return scm_i_clonebig (x
, 0);
687 else if (SCM_REALP (x
))
689 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
690 double xx
= SCM_REAL_VALUE (x
);
692 return scm_from_double (-xx
);
696 else if (SCM_FRACTIONP (x
))
698 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
700 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
701 SCM_FRACTION_DENOMINATOR (x
));
704 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
709 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
710 /* "Return the quotient of the numbers @var{x} and @var{y}."
713 scm_quotient (SCM x
, SCM y
)
717 long xx
= SCM_I_INUM (x
);
720 long yy
= SCM_I_INUM (y
);
722 scm_num_overflow (s_quotient
);
727 return SCM_I_MAKINUM (z
);
729 return scm_i_long2big (z
);
732 else if (SCM_BIGP (y
))
734 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
735 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
736 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
738 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
739 scm_remember_upto_here_1 (y
);
740 return SCM_I_MAKINUM (-1);
743 return SCM_I_MAKINUM (0);
746 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
748 else if (SCM_BIGP (x
))
752 long yy
= SCM_I_INUM (y
);
754 scm_num_overflow (s_quotient
);
759 SCM result
= scm_i_mkbig ();
762 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
765 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
768 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
769 scm_remember_upto_here_1 (x
);
770 return scm_i_normbig (result
);
773 else if (SCM_BIGP (y
))
775 SCM result
= scm_i_mkbig ();
776 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
779 scm_remember_upto_here_2 (x
, y
);
780 return scm_i_normbig (result
);
783 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
786 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
789 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
790 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
792 * "(remainder 13 4) @result{} 1\n"
793 * "(remainder -13 4) @result{} -1\n"
797 scm_remainder (SCM x
, SCM y
)
803 long yy
= SCM_I_INUM (y
);
805 scm_num_overflow (s_remainder
);
808 long z
= SCM_I_INUM (x
) % yy
;
809 return SCM_I_MAKINUM (z
);
812 else if (SCM_BIGP (y
))
814 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
815 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
816 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
818 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
819 scm_remember_upto_here_1 (y
);
820 return SCM_I_MAKINUM (0);
826 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
828 else if (SCM_BIGP (x
))
832 long yy
= SCM_I_INUM (y
);
834 scm_num_overflow (s_remainder
);
837 SCM result
= scm_i_mkbig ();
840 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
841 scm_remember_upto_here_1 (x
);
842 return scm_i_normbig (result
);
845 else if (SCM_BIGP (y
))
847 SCM result
= scm_i_mkbig ();
848 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
851 scm_remember_upto_here_2 (x
, y
);
852 return scm_i_normbig (result
);
855 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
858 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
862 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
863 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
865 * "(modulo 13 4) @result{} 1\n"
866 * "(modulo -13 4) @result{} 3\n"
870 scm_modulo (SCM x
, SCM y
)
874 long xx
= SCM_I_INUM (x
);
877 long yy
= SCM_I_INUM (y
);
879 scm_num_overflow (s_modulo
);
882 /* FIXME: I think this may be a bug on some arches -- results
883 of % with negative second arg are undefined... */
901 return SCM_I_MAKINUM (result
);
904 else if (SCM_BIGP (y
))
906 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
913 SCM pos_y
= scm_i_clonebig (y
, 0);
914 /* do this after the last scm_op */
915 mpz_init_set_si (z_x
, xx
);
916 result
= pos_y
; /* re-use this bignum */
917 mpz_mod (SCM_I_BIG_MPZ (result
),
919 SCM_I_BIG_MPZ (pos_y
));
920 scm_remember_upto_here_1 (pos_y
);
924 result
= scm_i_mkbig ();
925 /* do this after the last scm_op */
926 mpz_init_set_si (z_x
, xx
);
927 mpz_mod (SCM_I_BIG_MPZ (result
),
930 scm_remember_upto_here_1 (y
);
933 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
934 mpz_add (SCM_I_BIG_MPZ (result
),
936 SCM_I_BIG_MPZ (result
));
937 scm_remember_upto_here_1 (y
);
938 /* and do this before the next one */
940 return scm_i_normbig (result
);
944 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
946 else if (SCM_BIGP (x
))
950 long yy
= SCM_I_INUM (y
);
952 scm_num_overflow (s_modulo
);
955 SCM result
= scm_i_mkbig ();
956 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
958 (yy
< 0) ? - yy
: yy
);
959 scm_remember_upto_here_1 (x
);
960 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
961 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
962 SCM_I_BIG_MPZ (result
),
964 return scm_i_normbig (result
);
967 else if (SCM_BIGP (y
))
970 SCM result
= scm_i_mkbig ();
971 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
972 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
973 mpz_mod (SCM_I_BIG_MPZ (result
),
975 SCM_I_BIG_MPZ (pos_y
));
977 scm_remember_upto_here_1 (x
);
978 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
979 mpz_add (SCM_I_BIG_MPZ (result
),
981 SCM_I_BIG_MPZ (result
));
982 scm_remember_upto_here_2 (y
, pos_y
);
983 return scm_i_normbig (result
);
987 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
990 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
993 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
994 /* "Return the greatest common divisor of all arguments.\n"
995 * "If called without arguments, 0 is returned."
998 scm_gcd (SCM x
, SCM y
)
1001 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1003 if (SCM_I_INUMP (x
))
1005 if (SCM_I_INUMP (y
))
1007 long xx
= SCM_I_INUM (x
);
1008 long yy
= SCM_I_INUM (y
);
1009 long u
= xx
< 0 ? -xx
: xx
;
1010 long v
= yy
< 0 ? -yy
: yy
;
1020 /* Determine a common factor 2^k */
1021 while (!(1 & (u
| v
)))
1027 /* Now, any factor 2^n can be eliminated */
1047 return (SCM_POSFIXABLE (result
)
1048 ? SCM_I_MAKINUM (result
)
1049 : scm_i_long2big (result
));
1051 else if (SCM_BIGP (y
))
1057 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1059 else if (SCM_BIGP (x
))
1061 if (SCM_I_INUMP (y
))
1063 unsigned long result
;
1066 yy
= SCM_I_INUM (y
);
1071 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1072 scm_remember_upto_here_1 (x
);
1073 return (SCM_POSFIXABLE (result
)
1074 ? SCM_I_MAKINUM (result
)
1075 : scm_from_ulong (result
));
1077 else if (SCM_BIGP (y
))
1079 SCM result
= scm_i_mkbig ();
1080 mpz_gcd (SCM_I_BIG_MPZ (result
),
1083 scm_remember_upto_here_2 (x
, y
);
1084 return scm_i_normbig (result
);
1087 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1090 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1093 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1094 /* "Return the least common multiple of the arguments.\n"
1095 * "If called without arguments, 1 is returned."
1098 scm_lcm (SCM n1
, SCM n2
)
1100 if (SCM_UNBNDP (n2
))
1102 if (SCM_UNBNDP (n1
))
1103 return SCM_I_MAKINUM (1L);
1104 n2
= SCM_I_MAKINUM (1L);
1107 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1108 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1109 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1110 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1112 if (SCM_I_INUMP (n1
))
1114 if (SCM_I_INUMP (n2
))
1116 SCM d
= scm_gcd (n1
, n2
);
1117 if (scm_is_eq (d
, SCM_INUM0
))
1120 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1124 /* inum n1, big n2 */
1127 SCM result
= scm_i_mkbig ();
1128 long nn1
= SCM_I_INUM (n1
);
1129 if (nn1
== 0) return SCM_INUM0
;
1130 if (nn1
< 0) nn1
= - nn1
;
1131 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1132 scm_remember_upto_here_1 (n2
);
1140 if (SCM_I_INUMP (n2
))
1147 SCM result
= scm_i_mkbig ();
1148 mpz_lcm(SCM_I_BIG_MPZ (result
),
1150 SCM_I_BIG_MPZ (n2
));
1151 scm_remember_upto_here_2(n1
, n2
);
1152 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1158 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1163 + + + x (map digit:logand X Y)
1164 + - + x (map digit:logand X (lognot (+ -1 Y)))
1165 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1166 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1171 + + + (map digit:logior X Y)
1172 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1173 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1174 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1179 + + + (map digit:logxor X Y)
1180 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1181 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1182 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1187 + + (any digit:logand X Y)
1188 + - (any digit:logand X (lognot (+ -1 Y)))
1189 - + (any digit:logand (lognot (+ -1 X)) Y)
1194 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1196 "Return the bitwise AND of the integer arguments.\n\n"
1198 "(logand) @result{} -1\n"
1199 "(logand 7) @result{} 7\n"
1200 "(logand #b111 #b011 #b001) @result{} 1\n"
1202 #define FUNC_NAME s_scm_logand
1206 if (SCM_UNBNDP (n2
))
1208 if (SCM_UNBNDP (n1
))
1209 return SCM_I_MAKINUM (-1);
1210 else if (!SCM_NUMBERP (n1
))
1211 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1212 else if (SCM_NUMBERP (n1
))
1215 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1218 if (SCM_I_INUMP (n1
))
1220 nn1
= SCM_I_INUM (n1
);
1221 if (SCM_I_INUMP (n2
))
1223 long nn2
= SCM_I_INUM (n2
);
1224 return SCM_I_MAKINUM (nn1
& nn2
);
1226 else if SCM_BIGP (n2
)
1232 SCM result_z
= scm_i_mkbig ();
1234 mpz_init_set_si (nn1_z
, nn1
);
1235 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1236 scm_remember_upto_here_1 (n2
);
1238 return scm_i_normbig (result_z
);
1242 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1244 else if (SCM_BIGP (n1
))
1246 if (SCM_I_INUMP (n2
))
1249 nn1
= SCM_I_INUM (n1
);
1252 else if (SCM_BIGP (n2
))
1254 SCM result_z
= scm_i_mkbig ();
1255 mpz_and (SCM_I_BIG_MPZ (result_z
),
1257 SCM_I_BIG_MPZ (n2
));
1258 scm_remember_upto_here_2 (n1
, n2
);
1259 return scm_i_normbig (result_z
);
1262 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1265 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1270 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1272 "Return the bitwise OR of the integer arguments.\n\n"
1274 "(logior) @result{} 0\n"
1275 "(logior 7) @result{} 7\n"
1276 "(logior #b000 #b001 #b011) @result{} 3\n"
1278 #define FUNC_NAME s_scm_logior
1282 if (SCM_UNBNDP (n2
))
1284 if (SCM_UNBNDP (n1
))
1286 else if (SCM_NUMBERP (n1
))
1289 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1292 if (SCM_I_INUMP (n1
))
1294 nn1
= SCM_I_INUM (n1
);
1295 if (SCM_I_INUMP (n2
))
1297 long nn2
= SCM_I_INUM (n2
);
1298 return SCM_I_MAKINUM (nn1
| nn2
);
1300 else if (SCM_BIGP (n2
))
1306 SCM result_z
= scm_i_mkbig ();
1308 mpz_init_set_si (nn1_z
, nn1
);
1309 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1310 scm_remember_upto_here_1 (n2
);
1316 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1318 else if (SCM_BIGP (n1
))
1320 if (SCM_I_INUMP (n2
))
1323 nn1
= SCM_I_INUM (n1
);
1326 else if (SCM_BIGP (n2
))
1328 SCM result_z
= scm_i_mkbig ();
1329 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1331 SCM_I_BIG_MPZ (n2
));
1332 scm_remember_upto_here_2 (n1
, n2
);
1336 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1339 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1344 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1346 "Return the bitwise XOR of the integer arguments. A bit is\n"
1347 "set in the result if it is set in an odd number of arguments.\n"
1349 "(logxor) @result{} 0\n"
1350 "(logxor 7) @result{} 7\n"
1351 "(logxor #b000 #b001 #b011) @result{} 2\n"
1352 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1354 #define FUNC_NAME s_scm_logxor
1358 if (SCM_UNBNDP (n2
))
1360 if (SCM_UNBNDP (n1
))
1362 else if (SCM_NUMBERP (n1
))
1365 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1368 if (SCM_I_INUMP (n1
))
1370 nn1
= SCM_I_INUM (n1
);
1371 if (SCM_I_INUMP (n2
))
1373 long nn2
= SCM_I_INUM (n2
);
1374 return SCM_I_MAKINUM (nn1
^ nn2
);
1376 else if (SCM_BIGP (n2
))
1380 SCM result_z
= scm_i_mkbig ();
1382 mpz_init_set_si (nn1_z
, nn1
);
1383 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1384 scm_remember_upto_here_1 (n2
);
1386 return scm_i_normbig (result_z
);
1390 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1392 else if (SCM_BIGP (n1
))
1394 if (SCM_I_INUMP (n2
))
1397 nn1
= SCM_I_INUM (n1
);
1400 else if (SCM_BIGP (n2
))
1402 SCM result_z
= scm_i_mkbig ();
1403 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1405 SCM_I_BIG_MPZ (n2
));
1406 scm_remember_upto_here_2 (n1
, n2
);
1407 return scm_i_normbig (result_z
);
1410 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1413 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1418 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1421 "(logtest j k) @equiv{} (not (zero? (logand j k)))\n\n"
1422 "(logtest #b0100 #b1011) @result{} #f\n"
1423 "(logtest #b0100 #b0111) @result{} #t\n"
1425 #define FUNC_NAME s_scm_logtest
1429 if (SCM_I_INUMP (j
))
1431 nj
= SCM_I_INUM (j
);
1432 if (SCM_I_INUMP (k
))
1434 long nk
= SCM_I_INUM (k
);
1435 return scm_from_bool (nj
& nk
);
1437 else if (SCM_BIGP (k
))
1445 mpz_init_set_si (nj_z
, nj
);
1446 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1447 scm_remember_upto_here_1 (k
);
1448 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1454 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1456 else if (SCM_BIGP (j
))
1458 if (SCM_I_INUMP (k
))
1461 nj
= SCM_I_INUM (j
);
1464 else if (SCM_BIGP (k
))
1468 mpz_init (result_z
);
1472 scm_remember_upto_here_2 (j
, k
);
1473 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1474 mpz_clear (result_z
);
1478 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1481 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1486 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1489 "(logbit? index j) @equiv{} (logtest (integer-expt 2 index) j)\n\n"
1490 "(logbit? 0 #b1101) @result{} #t\n"
1491 "(logbit? 1 #b1101) @result{} #f\n"
1492 "(logbit? 2 #b1101) @result{} #t\n"
1493 "(logbit? 3 #b1101) @result{} #t\n"
1494 "(logbit? 4 #b1101) @result{} #f\n"
1496 #define FUNC_NAME s_scm_logbit_p
1498 unsigned long int iindex
;
1499 iindex
= scm_to_ulong (index
);
1501 if (SCM_I_INUMP (j
))
1503 /* bits above what's in an inum follow the sign bit */
1504 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1505 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1507 else if (SCM_BIGP (j
))
1509 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1510 scm_remember_upto_here_1 (j
);
1511 return scm_from_bool (val
);
1514 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1519 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1521 "Return the integer which is the ones-complement of the integer\n"
1525 "(number->string (lognot #b10000000) 2)\n"
1526 " @result{} \"-10000001\"\n"
1527 "(number->string (lognot #b0) 2)\n"
1528 " @result{} \"-1\"\n"
1530 #define FUNC_NAME s_scm_lognot
1532 if (SCM_I_INUMP (n
)) {
1533 /* No overflow here, just need to toggle all the bits making up the inum.
1534 Enhancement: No need to strip the tag and add it back, could just xor
1535 a block of 1 bits, if that worked with the various debug versions of
1537 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1539 } else if (SCM_BIGP (n
)) {
1540 SCM result
= scm_i_mkbig ();
1541 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1542 scm_remember_upto_here_1 (n
);
1546 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1551 /* returns 0 if IN is not an integer. OUT must already be
1554 coerce_to_big (SCM in
, mpz_t out
)
1557 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1558 else if (SCM_I_INUMP (in
))
1559 mpz_set_si (out
, SCM_I_INUM (in
));
1566 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1567 (SCM n
, SCM k
, SCM m
),
1568 "Return @var{n} raised to the integer exponent\n"
1569 "@var{k}, modulo @var{m}.\n"
1572 "(modulo-expt 2 3 5)\n"
1575 #define FUNC_NAME s_scm_modulo_expt
1581 /* There are two classes of error we might encounter --
1582 1) Math errors, which we'll report by calling scm_num_overflow,
1584 2) wrong-type errors, which of course we'll report by calling
1586 We don't report those errors immediately, however; instead we do
1587 some cleanup first. These variables tell us which error (if
1588 any) we should report after cleaning up.
1590 int report_overflow
= 0;
1592 int position_of_wrong_type
= 0;
1593 SCM value_of_wrong_type
= SCM_INUM0
;
1595 SCM result
= SCM_UNDEFINED
;
1601 if (scm_is_eq (m
, SCM_INUM0
))
1603 report_overflow
= 1;
1607 if (!coerce_to_big (n
, n_tmp
))
1609 value_of_wrong_type
= n
;
1610 position_of_wrong_type
= 1;
1614 if (!coerce_to_big (k
, k_tmp
))
1616 value_of_wrong_type
= k
;
1617 position_of_wrong_type
= 2;
1621 if (!coerce_to_big (m
, m_tmp
))
1623 value_of_wrong_type
= m
;
1624 position_of_wrong_type
= 3;
1628 /* if the exponent K is negative, and we simply call mpz_powm, we
1629 will get a divide-by-zero exception when an inverse 1/n mod m
1630 doesn't exist (or is not unique). Since exceptions are hard to
1631 handle, we'll attempt the inversion "by hand" -- that way, we get
1632 a simple failure code, which is easy to handle. */
1634 if (-1 == mpz_sgn (k_tmp
))
1636 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1638 report_overflow
= 1;
1641 mpz_neg (k_tmp
, k_tmp
);
1644 result
= scm_i_mkbig ();
1645 mpz_powm (SCM_I_BIG_MPZ (result
),
1650 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1651 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1658 if (report_overflow
)
1659 scm_num_overflow (FUNC_NAME
);
1661 if (position_of_wrong_type
)
1662 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1663 value_of_wrong_type
);
1665 return scm_i_normbig (result
);
1669 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1671 "Return @var{n} raised to the non-negative integer exponent\n"
1675 "(integer-expt 2 5)\n"
1677 "(integer-expt -3 3)\n"
1680 #define FUNC_NAME s_scm_integer_expt
1683 SCM z_i2
= SCM_BOOL_F
;
1685 SCM acc
= SCM_I_MAKINUM (1L);
1687 /* 0^0 == 1 according to R5RS */
1688 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1689 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1690 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1691 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1693 if (SCM_I_INUMP (k
))
1694 i2
= SCM_I_INUM (k
);
1695 else if (SCM_BIGP (k
))
1697 z_i2
= scm_i_clonebig (k
, 1);
1698 scm_remember_upto_here_1 (k
);
1701 else if (SCM_REALP (k
))
1703 double r
= SCM_REAL_VALUE (k
);
1704 if (floor (r
) != r
|| xisinf (r
))
1705 SCM_WRONG_TYPE_ARG (2, k
);
1706 if ((r
> SCM_MOST_POSITIVE_FIXNUM
) || (r
< SCM_MOST_NEGATIVE_FIXNUM
))
1708 z_i2
= scm_i_mkbig ();
1709 mpz_set_d (SCM_I_BIG_MPZ (z_i2
), r
);
1718 SCM_WRONG_TYPE_ARG (2, k
);
1722 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1724 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1725 n
= scm_divide (n
, SCM_UNDEFINED
);
1729 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1733 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1735 return scm_product (acc
, n
);
1737 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1738 acc
= scm_product (acc
, n
);
1739 n
= scm_product (n
, n
);
1740 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1748 n
= scm_divide (n
, SCM_UNDEFINED
);
1755 return scm_product (acc
, n
);
1757 acc
= scm_product (acc
, n
);
1758 n
= scm_product (n
, n
);
1765 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1767 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1768 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1770 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1771 "@var{cnt} is negative it's a division, rounded towards negative\n"
1772 "infinity. (Note that this is not the same rounding as\n"
1773 "@code{quotient} does.)\n"
1775 "With @var{n} viewed as an infinite precision twos complement,\n"
1776 "@code{ash} means a left shift introducing zero bits, or a right\n"
1777 "shift dropping bits.\n"
1780 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1781 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1783 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1784 "(ash -23 -2) @result{} -6\n"
1786 #define FUNC_NAME s_scm_ash
1789 bits_to_shift
= scm_to_long (cnt
);
1791 if (bits_to_shift
< 0)
1793 /* Shift right by abs(cnt) bits. This is realized as a division
1794 by div:=2^abs(cnt). However, to guarantee the floor
1795 rounding, negative values require some special treatment.
1797 SCM div
= scm_integer_expt (SCM_I_MAKINUM (2),
1798 scm_from_long (-bits_to_shift
));
1800 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1801 if (scm_is_false (scm_negative_p (n
)))
1802 return scm_quotient (n
, div
);
1804 return scm_sum (SCM_I_MAKINUM (-1L),
1805 scm_quotient (scm_sum (SCM_I_MAKINUM (1L), n
), div
));
1808 /* Shift left is done by multiplication with 2^CNT */
1809 return scm_product (n
, scm_integer_expt (SCM_I_MAKINUM (2), cnt
));
1814 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1815 (SCM n
, SCM start
, SCM end
),
1816 "Return the integer composed of the @var{start} (inclusive)\n"
1817 "through @var{end} (exclusive) bits of @var{n}. The\n"
1818 "@var{start}th bit becomes the 0-th bit in the result.\n"
1821 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1822 " @result{} \"1010\"\n"
1823 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1824 " @result{} \"10110\"\n"
1826 #define FUNC_NAME s_scm_bit_extract
1828 unsigned long int istart
, iend
, bits
;
1829 istart
= scm_to_ulong (start
);
1830 iend
= scm_to_ulong (end
);
1831 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1833 /* how many bits to keep */
1834 bits
= iend
- istart
;
1836 if (SCM_I_INUMP (n
))
1838 long int in
= SCM_I_INUM (n
);
1840 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1841 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1842 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1844 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1846 /* Since we emulate two's complement encoded numbers, this
1847 * special case requires us to produce a result that has
1848 * more bits than can be stored in a fixnum.
1850 SCM result
= scm_i_long2big (in
);
1851 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1856 /* mask down to requisite bits */
1857 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1858 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1860 else if (SCM_BIGP (n
))
1865 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1869 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1870 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1871 such bits into a ulong. */
1872 result
= scm_i_mkbig ();
1873 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1874 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1875 result
= scm_i_normbig (result
);
1877 scm_remember_upto_here_1 (n
);
1881 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1886 static const char scm_logtab
[] = {
1887 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1890 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1892 "Return the number of bits in integer @var{n}. If integer is\n"
1893 "positive, the 1-bits in its binary representation are counted.\n"
1894 "If negative, the 0-bits in its two's-complement binary\n"
1895 "representation are counted. If 0, 0 is returned.\n"
1898 "(logcount #b10101010)\n"
1905 #define FUNC_NAME s_scm_logcount
1907 if (SCM_I_INUMP (n
))
1909 unsigned long int c
= 0;
1910 long int nn
= SCM_I_INUM (n
);
1915 c
+= scm_logtab
[15 & nn
];
1918 return SCM_I_MAKINUM (c
);
1920 else if (SCM_BIGP (n
))
1922 unsigned long count
;
1923 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1924 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1926 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1927 scm_remember_upto_here_1 (n
);
1928 return SCM_I_MAKINUM (count
);
1931 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1936 static const char scm_ilentab
[] = {
1937 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1941 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1943 "Return the number of bits necessary to represent @var{n}.\n"
1946 "(integer-length #b10101010)\n"
1948 "(integer-length 0)\n"
1950 "(integer-length #b1111)\n"
1953 #define FUNC_NAME s_scm_integer_length
1955 if (SCM_I_INUMP (n
))
1957 unsigned long int c
= 0;
1959 long int nn
= SCM_I_INUM (n
);
1965 l
= scm_ilentab
[15 & nn
];
1968 return SCM_I_MAKINUM (c
- 4 + l
);
1970 else if (SCM_BIGP (n
))
1972 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1973 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1974 1 too big, so check for that and adjust. */
1975 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1976 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1977 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1978 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1980 scm_remember_upto_here_1 (n
);
1981 return SCM_I_MAKINUM (size
);
1984 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1988 /*** NUMBERS -> STRINGS ***/
1989 #define SCM_MAX_DBL_PREC 60
1990 #define SCM_MAX_DBL_RADIX 36
1992 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
1993 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
1994 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
1997 void init_dblprec(int *prec
, int radix
) {
1998 /* determine floating point precision by adding successively
1999 smaller increments to 1.0 until it is considered == 1.0 */
2000 double f
= ((double)1.0)/radix
;
2001 double fsum
= 1.0 + f
;
2006 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2018 void init_fx_radix(double *fx_list
, int radix
)
2020 /* initialize a per-radix list of tolerances. When added
2021 to a number < 1.0, we can determine if we should raund
2022 up and quit converting a number to a string. */
2026 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2027 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2030 /* use this array as a way to generate a single digit */
2031 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2034 idbl2str (double f
, char *a
, int radix
)
2036 int efmt
, dpt
, d
, i
, wp
;
2038 #ifdef DBL_MIN_10_EXP
2041 #endif /* DBL_MIN_10_EXP */
2046 radix
> SCM_MAX_DBL_RADIX
)
2048 /* revert to existing behavior */
2052 wp
= scm_dblprec
[radix
-2];
2053 fx
= fx_per_radix
[radix
-2];
2057 #ifdef HAVE_COPYSIGN
2058 double sgn
= copysign (1.0, f
);
2063 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2069 strcpy (a
, "-inf.0");
2071 strcpy (a
, "+inf.0");
2074 else if (xisnan (f
))
2076 strcpy (a
, "+nan.0");
2086 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2087 make-uniform-vector, from causing infinite loops. */
2088 /* just do the checking...if it passes, we do the conversion for our
2089 radix again below */
2096 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2104 while (f_cpy
> 10.0)
2107 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2128 if (f
+ fx
[wp
] >= radix
)
2135 /* adding 9999 makes this equivalent to abs(x) % 3 */
2136 dpt
= (exp
+ 9999) % 3;
2140 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2162 a
[ch
++] = number_chars
[d
];
2165 if (f
+ fx
[wp
] >= 1.0)
2167 a
[ch
- 1] = number_chars
[d
+1];
2179 if ((dpt
> 4) && (exp
> 6))
2181 d
= (a
[0] == '-' ? 2 : 1);
2182 for (i
= ch
++; i
> d
; i
--)
2195 if (a
[ch
- 1] == '.')
2196 a
[ch
++] = '0'; /* trailing zero */
2205 for (i
= radix
; i
<= exp
; i
*= radix
);
2206 for (i
/= radix
; i
; i
/= radix
)
2208 a
[ch
++] = number_chars
[exp
/ i
];
2216 iflo2str (SCM flt
, char *str
, int radix
)
2219 if (SCM_REALP (flt
))
2220 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2223 i
= idbl2str (SCM_COMPLEX_REAL (flt
), str
, radix
);
2224 if (SCM_COMPLEX_IMAG (flt
) != 0.0)
2226 double imag
= SCM_COMPLEX_IMAG (flt
);
2227 /* Don't output a '+' for negative numbers or for Inf and
2228 NaN. They will provide their own sign. */
2229 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2231 i
+= idbl2str (imag
, &str
[i
], radix
);
2238 /* convert a long to a string (unterminated). returns the number of
2239 characters in the result.
2241 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2243 scm_iint2str (long num
, int rad
, char *p
)
2247 unsigned long n
= (num
< 0) ? -num
: num
;
2249 for (n
/= rad
; n
> 0; n
/= rad
)
2266 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2271 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2273 "Return a string holding the external representation of the\n"
2274 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2275 "inexact, a radix of 10 will be used.")
2276 #define FUNC_NAME s_scm_number_to_string
2280 if (SCM_UNBNDP (radix
))
2283 base
= scm_to_signed_integer (radix
, 2, 36);
2285 if (SCM_I_INUMP (n
))
2287 char num_buf
[SCM_INTBUFLEN
];
2288 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2289 return scm_from_locale_stringn (num_buf
, length
);
2291 else if (SCM_BIGP (n
))
2293 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2294 scm_remember_upto_here_1 (n
);
2295 return scm_take_locale_string (str
);
2297 else if (SCM_FRACTIONP (n
))
2299 scm_i_fraction_reduce (n
);
2300 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2301 scm_from_locale_string ("/"),
2302 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2304 else if (SCM_INEXACTP (n
))
2306 char num_buf
[FLOBUFLEN
];
2307 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2310 SCM_WRONG_TYPE_ARG (1, n
);
2315 /* These print routines used to be stubbed here so that scm_repl.c
2316 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2319 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2321 char num_buf
[FLOBUFLEN
];
2322 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2327 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2330 char num_buf
[FLOBUFLEN
];
2331 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2336 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2339 scm_i_fraction_reduce (sexp
);
2340 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2341 scm_lfwrite (scm_i_string_chars (str
), scm_i_string_length (str
), port
);
2342 scm_remember_upto_here_1 (str
);
2347 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2349 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2350 scm_remember_upto_here_1 (exp
);
2351 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2355 /*** END nums->strs ***/
2358 /*** STRINGS -> NUMBERS ***/
2360 /* The following functions implement the conversion from strings to numbers.
2361 * The implementation somehow follows the grammar for numbers as it is given
2362 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2363 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2364 * points should be noted about the implementation:
2365 * * Each function keeps a local index variable 'idx' that points at the
2366 * current position within the parsed string. The global index is only
2367 * updated if the function could parse the corresponding syntactic unit
2369 * * Similarly, the functions keep track of indicators of inexactness ('#',
2370 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2371 * global exactness information is only updated after each part has been
2372 * successfully parsed.
2373 * * Sequences of digits are parsed into temporary variables holding fixnums.
2374 * Only if these fixnums would overflow, the result variables are updated
2375 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2376 * the temporary variables holding the fixnums are cleared, and the process
2377 * starts over again. If for example fixnums were able to store five decimal
2378 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2379 * and the result was computed as 12345 * 100000 + 67890. In other words,
2380 * only every five digits two bignum operations were performed.
2383 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2385 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2387 /* In non ASCII-style encodings the following macro might not work. */
2388 #define XDIGIT2UINT(d) \
2389 (isdigit ((int) (unsigned char) d) \
2391 : tolower ((int) (unsigned char) d) - 'a' + 10)
2394 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2395 unsigned int radix
, enum t_exactness
*p_exactness
)
2397 unsigned int idx
= *p_idx
;
2398 unsigned int hash_seen
= 0;
2399 scm_t_bits shift
= 1;
2401 unsigned int digit_value
;
2409 if (!isxdigit ((int) (unsigned char) c
))
2411 digit_value
= XDIGIT2UINT (c
);
2412 if (digit_value
>= radix
)
2416 result
= SCM_I_MAKINUM (digit_value
);
2420 if (isxdigit ((int) (unsigned char) c
))
2424 digit_value
= XDIGIT2UINT (c
);
2425 if (digit_value
>= radix
)
2437 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2439 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2441 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2448 shift
= shift
* radix
;
2449 add
= add
* radix
+ digit_value
;
2454 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2456 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2460 *p_exactness
= INEXACT
;
2466 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2467 * covers the parts of the rules that start at a potential point. The value
2468 * of the digits up to the point have been parsed by the caller and are given
2469 * in variable result. The content of *p_exactness indicates, whether a hash
2470 * has already been seen in the digits before the point.
2473 /* In non ASCII-style encodings the following macro might not work. */
2474 #define DIGIT2UINT(d) ((d) - '0')
2477 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2478 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2480 unsigned int idx
= *p_idx
;
2481 enum t_exactness x
= *p_exactness
;
2486 if (mem
[idx
] == '.')
2488 scm_t_bits shift
= 1;
2490 unsigned int digit_value
;
2491 SCM big_shift
= SCM_I_MAKINUM (1);
2497 if (isdigit ((int) (unsigned char) c
))
2502 digit_value
= DIGIT2UINT (c
);
2513 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2515 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2516 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2518 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2526 add
= add
* 10 + digit_value
;
2532 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2533 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2534 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2537 result
= scm_divide (result
, big_shift
);
2539 /* We've seen a decimal point, thus the value is implicitly inexact. */
2551 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2578 if (!isdigit ((int) (unsigned char) c
))
2582 exponent
= DIGIT2UINT (c
);
2586 if (isdigit ((int) (unsigned char) c
))
2589 if (exponent
<= SCM_MAXEXP
)
2590 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2596 if (exponent
> SCM_MAXEXP
)
2598 size_t exp_len
= idx
- start
;
2599 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2600 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2601 scm_out_of_range ("string->number", exp_num
);
2604 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2606 result
= scm_product (result
, e
);
2608 result
= scm_divide2real (result
, e
);
2610 /* We've seen an exponent, thus the value is implicitly inexact. */
2628 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2631 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2632 unsigned int radix
, enum t_exactness
*p_exactness
)
2634 unsigned int idx
= *p_idx
;
2640 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2646 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2648 enum t_exactness x
= EXACT
;
2650 /* Cobble up the fractional part. We might want to set the
2651 NaN's mantissa from it. */
2653 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2658 if (mem
[idx
] == '.')
2662 else if (idx
+ 1 == len
)
2664 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2667 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2668 p_idx
, p_exactness
);
2672 enum t_exactness x
= EXACT
;
2675 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2676 if (scm_is_false (uinteger
))
2681 else if (mem
[idx
] == '/')
2687 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2688 if (scm_is_false (divisor
))
2691 /* both are int/big here, I assume */
2692 result
= scm_i_make_ratio (uinteger
, divisor
);
2694 else if (radix
== 10)
2696 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2697 if (scm_is_false (result
))
2708 /* When returning an inexact zero, make sure it is represented as a
2709 floating point value so that we can change its sign.
2711 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2712 result
= scm_from_double (0.0);
2718 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2721 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2722 unsigned int radix
, enum t_exactness
*p_exactness
)
2746 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2747 if (scm_is_false (ureal
))
2749 /* input must be either +i or -i */
2754 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2760 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2767 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2768 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2777 /* either +<ureal>i or -<ureal>i */
2784 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2787 /* polar input: <real>@<real>. */
2812 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2813 if (scm_is_false (angle
))
2818 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2819 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2821 result
= scm_make_polar (ureal
, angle
);
2826 /* expecting input matching <real>[+-]<ureal>?i */
2833 int sign
= (c
== '+') ? 1 : -1;
2834 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2836 if (scm_is_false (imag
))
2837 imag
= SCM_I_MAKINUM (sign
);
2838 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2839 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2843 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2850 return scm_make_rectangular (ureal
, imag
);
2859 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2861 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2864 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2866 unsigned int idx
= 0;
2867 unsigned int radix
= NO_RADIX
;
2868 enum t_exactness forced_x
= NO_EXACTNESS
;
2869 enum t_exactness implicit_x
= EXACT
;
2872 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2873 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2875 switch (mem
[idx
+ 1])
2878 if (radix
!= NO_RADIX
)
2883 if (radix
!= NO_RADIX
)
2888 if (forced_x
!= NO_EXACTNESS
)
2893 if (forced_x
!= NO_EXACTNESS
)
2898 if (radix
!= NO_RADIX
)
2903 if (radix
!= NO_RADIX
)
2913 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2914 if (radix
== NO_RADIX
)
2915 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2917 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2919 if (scm_is_false (result
))
2925 if (SCM_INEXACTP (result
))
2926 return scm_inexact_to_exact (result
);
2930 if (SCM_INEXACTP (result
))
2933 return scm_exact_to_inexact (result
);
2936 if (implicit_x
== INEXACT
)
2938 if (SCM_INEXACTP (result
))
2941 return scm_exact_to_inexact (result
);
2949 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2950 (SCM string
, SCM radix
),
2951 "Return a number of the maximally precise representation\n"
2952 "expressed by the given @var{string}. @var{radix} must be an\n"
2953 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2954 "is a default radix that may be overridden by an explicit radix\n"
2955 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2956 "supplied, then the default radix is 10. If string is not a\n"
2957 "syntactically valid notation for a number, then\n"
2958 "@code{string->number} returns @code{#f}.")
2959 #define FUNC_NAME s_scm_string_to_number
2963 SCM_VALIDATE_STRING (1, string
);
2965 if (SCM_UNBNDP (radix
))
2968 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
2970 answer
= scm_i_mem2number (scm_i_string_chars (string
),
2971 scm_i_string_length (string
),
2973 scm_remember_upto_here_1 (string
);
2979 /*** END strs->nums ***/
2983 scm_bigequal (SCM x
, SCM y
)
2985 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2986 scm_remember_upto_here_2 (x
, y
);
2987 return scm_from_bool (0 == result
);
2991 scm_real_equalp (SCM x
, SCM y
)
2993 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
2997 scm_complex_equalp (SCM x
, SCM y
)
2999 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3000 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3004 scm_i_fraction_equalp (SCM x
, SCM y
)
3006 scm_i_fraction_reduce (x
);
3007 scm_i_fraction_reduce (y
);
3008 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3009 SCM_FRACTION_NUMERATOR (y
)))
3010 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3011 SCM_FRACTION_DENOMINATOR (y
))))
3018 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3020 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3022 #define FUNC_NAME s_scm_number_p
3024 return scm_from_bool (SCM_NUMBERP (x
));
3028 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3030 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3031 "otherwise. Note that the sets of real, rational and integer\n"
3032 "values form subsets of the set of complex numbers, i. e. the\n"
3033 "predicate will also be fulfilled if @var{x} is a real,\n"
3034 "rational or integer number.")
3035 #define FUNC_NAME s_scm_complex_p
3037 /* all numbers are complex. */
3038 return scm_number_p (x
);
3042 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3044 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3045 "otherwise. Note that the set of integer values forms a subset of\n"
3046 "the set of real numbers, i. e. the predicate will also be\n"
3047 "fulfilled if @var{x} is an integer number.")
3048 #define FUNC_NAME s_scm_real_p
3050 /* we can't represent irrational numbers. */
3051 return scm_rational_p (x
);
3055 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3057 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3058 "otherwise. Note that the set of integer values forms a subset of\n"
3059 "the set of rational numbers, i. e. the predicate will also be\n"
3060 "fulfilled if @var{x} is an integer number.")
3061 #define FUNC_NAME s_scm_rational_p
3063 if (SCM_I_INUMP (x
))
3065 else if (SCM_IMP (x
))
3067 else if (SCM_BIGP (x
))
3069 else if (SCM_FRACTIONP (x
))
3071 else if (SCM_REALP (x
))
3072 /* due to their limited precision, all floating point numbers are
3073 rational as well. */
3080 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3082 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3084 #define FUNC_NAME s_scm_integer_p
3087 if (SCM_I_INUMP (x
))
3093 if (!SCM_INEXACTP (x
))
3095 if (SCM_COMPLEXP (x
))
3097 r
= SCM_REAL_VALUE (x
);
3098 /* +/-inf passes r==floor(r), making those #t */
3106 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3108 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3110 #define FUNC_NAME s_scm_inexact_p
3112 if (SCM_INEXACTP (x
))
3114 if (SCM_NUMBERP (x
))
3116 SCM_WRONG_TYPE_ARG (1, x
);
3121 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3122 /* "Return @code{#t} if all parameters are numerically equal." */
3124 scm_num_eq_p (SCM x
, SCM y
)
3127 if (SCM_I_INUMP (x
))
3129 long xx
= SCM_I_INUM (x
);
3130 if (SCM_I_INUMP (y
))
3132 long yy
= SCM_I_INUM (y
);
3133 return scm_from_bool (xx
== yy
);
3135 else if (SCM_BIGP (y
))
3137 else if (SCM_REALP (y
))
3138 return scm_from_bool ((double) xx
== SCM_REAL_VALUE (y
));
3139 else if (SCM_COMPLEXP (y
))
3140 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3141 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3142 else if (SCM_FRACTIONP (y
))
3145 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3147 else if (SCM_BIGP (x
))
3149 if (SCM_I_INUMP (y
))
3151 else if (SCM_BIGP (y
))
3153 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3154 scm_remember_upto_here_2 (x
, y
);
3155 return scm_from_bool (0 == cmp
);
3157 else if (SCM_REALP (y
))
3160 if (xisnan (SCM_REAL_VALUE (y
)))
3162 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3163 scm_remember_upto_here_1 (x
);
3164 return scm_from_bool (0 == cmp
);
3166 else if (SCM_COMPLEXP (y
))
3169 if (0.0 != SCM_COMPLEX_IMAG (y
))
3171 if (xisnan (SCM_COMPLEX_REAL (y
)))
3173 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3174 scm_remember_upto_here_1 (x
);
3175 return scm_from_bool (0 == cmp
);
3177 else if (SCM_FRACTIONP (y
))
3180 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3182 else if (SCM_REALP (x
))
3184 if (SCM_I_INUMP (y
))
3185 return scm_from_bool (SCM_REAL_VALUE (x
) == (double) SCM_I_INUM (y
));
3186 else if (SCM_BIGP (y
))
3189 if (xisnan (SCM_REAL_VALUE (x
)))
3191 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3192 scm_remember_upto_here_1 (y
);
3193 return scm_from_bool (0 == cmp
);
3195 else if (SCM_REALP (y
))
3196 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3197 else if (SCM_COMPLEXP (y
))
3198 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3199 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3200 else if (SCM_FRACTIONP (y
))
3202 double xx
= SCM_REAL_VALUE (x
);
3206 return scm_from_bool (xx
< 0.0);
3207 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3211 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3213 else if (SCM_COMPLEXP (x
))
3215 if (SCM_I_INUMP (y
))
3216 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3217 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3218 else if (SCM_BIGP (y
))
3221 if (0.0 != SCM_COMPLEX_IMAG (x
))
3223 if (xisnan (SCM_COMPLEX_REAL (x
)))
3225 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3226 scm_remember_upto_here_1 (y
);
3227 return scm_from_bool (0 == cmp
);
3229 else if (SCM_REALP (y
))
3230 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3231 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3232 else if (SCM_COMPLEXP (y
))
3233 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3234 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3235 else if (SCM_FRACTIONP (y
))
3238 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3240 xx
= SCM_COMPLEX_REAL (x
);
3244 return scm_from_bool (xx
< 0.0);
3245 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3249 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3251 else if (SCM_FRACTIONP (x
))
3253 if (SCM_I_INUMP (y
))
3255 else if (SCM_BIGP (y
))
3257 else if (SCM_REALP (y
))
3259 double yy
= SCM_REAL_VALUE (y
);
3263 return scm_from_bool (0.0 < yy
);
3264 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3267 else if (SCM_COMPLEXP (y
))
3270 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3272 yy
= SCM_COMPLEX_REAL (y
);
3276 return scm_from_bool (0.0 < yy
);
3277 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3280 else if (SCM_FRACTIONP (y
))
3281 return scm_i_fraction_equalp (x
, y
);
3283 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3286 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3290 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3291 done are good for inums, but for bignums an answer can almost always be
3292 had by just examining a few high bits of the operands, as done by GMP in
3293 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3294 of the float exponent to take into account. */
3296 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3297 /* "Return @code{#t} if the list of parameters is monotonically\n"
3301 scm_less_p (SCM x
, SCM y
)
3304 if (SCM_I_INUMP (x
))
3306 long xx
= SCM_I_INUM (x
);
3307 if (SCM_I_INUMP (y
))
3309 long yy
= SCM_I_INUM (y
);
3310 return scm_from_bool (xx
< yy
);
3312 else if (SCM_BIGP (y
))
3314 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3315 scm_remember_upto_here_1 (y
);
3316 return scm_from_bool (sgn
> 0);
3318 else if (SCM_REALP (y
))
3319 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3320 else if (SCM_FRACTIONP (y
))
3322 /* "x < a/b" becomes "x*b < a" */
3324 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3325 y
= SCM_FRACTION_NUMERATOR (y
);
3329 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3331 else if (SCM_BIGP (x
))
3333 if (SCM_I_INUMP (y
))
3335 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3336 scm_remember_upto_here_1 (x
);
3337 return scm_from_bool (sgn
< 0);
3339 else if (SCM_BIGP (y
))
3341 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3342 scm_remember_upto_here_2 (x
, y
);
3343 return scm_from_bool (cmp
< 0);
3345 else if (SCM_REALP (y
))
3348 if (xisnan (SCM_REAL_VALUE (y
)))
3350 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3351 scm_remember_upto_here_1 (x
);
3352 return scm_from_bool (cmp
< 0);
3354 else if (SCM_FRACTIONP (y
))
3357 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3359 else if (SCM_REALP (x
))
3361 if (SCM_I_INUMP (y
))
3362 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3363 else if (SCM_BIGP (y
))
3366 if (xisnan (SCM_REAL_VALUE (x
)))
3368 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3369 scm_remember_upto_here_1 (y
);
3370 return scm_from_bool (cmp
> 0);
3372 else if (SCM_REALP (y
))
3373 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3374 else if (SCM_FRACTIONP (y
))
3376 double xx
= SCM_REAL_VALUE (x
);
3380 return scm_from_bool (xx
< 0.0);
3381 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3385 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3387 else if (SCM_FRACTIONP (x
))
3389 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3391 /* "a/b < y" becomes "a < y*b" */
3392 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3393 x
= SCM_FRACTION_NUMERATOR (x
);
3396 else if (SCM_REALP (y
))
3398 double yy
= SCM_REAL_VALUE (y
);
3402 return scm_from_bool (0.0 < yy
);
3403 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3406 else if (SCM_FRACTIONP (y
))
3408 /* "a/b < c/d" becomes "a*d < c*b" */
3409 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3410 SCM_FRACTION_DENOMINATOR (y
));
3411 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3412 SCM_FRACTION_DENOMINATOR (x
));
3418 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3421 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3425 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3426 /* "Return @code{#t} if the list of parameters is monotonically\n"
3429 #define FUNC_NAME s_scm_gr_p
3431 scm_gr_p (SCM x
, SCM y
)
3433 if (!SCM_NUMBERP (x
))
3434 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3435 else if (!SCM_NUMBERP (y
))
3436 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3438 return scm_less_p (y
, x
);
3443 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3444 /* "Return @code{#t} if the list of parameters is monotonically\n"
3447 #define FUNC_NAME s_scm_leq_p
3449 scm_leq_p (SCM x
, SCM y
)
3451 if (!SCM_NUMBERP (x
))
3452 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3453 else if (!SCM_NUMBERP (y
))
3454 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3455 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3458 return scm_not (scm_less_p (y
, x
));
3463 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3464 /* "Return @code{#t} if the list of parameters is monotonically\n"
3467 #define FUNC_NAME s_scm_geq_p
3469 scm_geq_p (SCM x
, SCM y
)
3471 if (!SCM_NUMBERP (x
))
3472 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3473 else if (!SCM_NUMBERP (y
))
3474 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3475 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3478 return scm_not (scm_less_p (x
, y
));
3483 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3484 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3490 if (SCM_I_INUMP (z
))
3491 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3492 else if (SCM_BIGP (z
))
3494 else if (SCM_REALP (z
))
3495 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3496 else if (SCM_COMPLEXP (z
))
3497 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3498 && SCM_COMPLEX_IMAG (z
) == 0.0);
3499 else if (SCM_FRACTIONP (z
))
3502 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3506 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3507 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3511 scm_positive_p (SCM x
)
3513 if (SCM_I_INUMP (x
))
3514 return scm_from_bool (SCM_I_INUM (x
) > 0);
3515 else if (SCM_BIGP (x
))
3517 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3518 scm_remember_upto_here_1 (x
);
3519 return scm_from_bool (sgn
> 0);
3521 else if (SCM_REALP (x
))
3522 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3523 else if (SCM_FRACTIONP (x
))
3524 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3526 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3530 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3531 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3535 scm_negative_p (SCM x
)
3537 if (SCM_I_INUMP (x
))
3538 return scm_from_bool (SCM_I_INUM (x
) < 0);
3539 else if (SCM_BIGP (x
))
3541 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3542 scm_remember_upto_here_1 (x
);
3543 return scm_from_bool (sgn
< 0);
3545 else if (SCM_REALP (x
))
3546 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3547 else if (SCM_FRACTIONP (x
))
3548 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3550 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3554 /* scm_min and scm_max return an inexact when either argument is inexact, as
3555 required by r5rs. On that basis, for exact/inexact combinations the
3556 exact is converted to inexact to compare and possibly return. This is
3557 unlike scm_less_p above which takes some trouble to preserve all bits in
3558 its test, such trouble is not required for min and max. */
3560 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3561 /* "Return the maximum of all parameter values."
3564 scm_max (SCM x
, SCM y
)
3569 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3570 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3573 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3576 if (SCM_I_INUMP (x
))
3578 long xx
= SCM_I_INUM (x
);
3579 if (SCM_I_INUMP (y
))
3581 long yy
= SCM_I_INUM (y
);
3582 return (xx
< yy
) ? y
: x
;
3584 else if (SCM_BIGP (y
))
3586 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3587 scm_remember_upto_here_1 (y
);
3588 return (sgn
< 0) ? x
: y
;
3590 else if (SCM_REALP (y
))
3593 /* if y==NaN then ">" is false and we return NaN */
3594 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3596 else if (SCM_FRACTIONP (y
))
3599 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3602 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3604 else if (SCM_BIGP (x
))
3606 if (SCM_I_INUMP (y
))
3608 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3609 scm_remember_upto_here_1 (x
);
3610 return (sgn
< 0) ? y
: x
;
3612 else if (SCM_BIGP (y
))
3614 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3615 scm_remember_upto_here_2 (x
, y
);
3616 return (cmp
> 0) ? x
: y
;
3618 else if (SCM_REALP (y
))
3620 /* if y==NaN then xx>yy is false, so we return the NaN y */
3623 xx
= scm_i_big2dbl (x
);
3624 yy
= SCM_REAL_VALUE (y
);
3625 return (xx
> yy
? scm_from_double (xx
) : y
);
3627 else if (SCM_FRACTIONP (y
))
3632 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3634 else if (SCM_REALP (x
))
3636 if (SCM_I_INUMP (y
))
3638 double z
= SCM_I_INUM (y
);
3639 /* if x==NaN then "<" is false and we return NaN */
3640 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3642 else if (SCM_BIGP (y
))
3647 else if (SCM_REALP (y
))
3649 /* if x==NaN then our explicit check means we return NaN
3650 if y==NaN then ">" is false and we return NaN
3651 calling isnan is unavoidable, since it's the only way to know
3652 which of x or y causes any compares to be false */
3653 double xx
= SCM_REAL_VALUE (x
);
3654 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3656 else if (SCM_FRACTIONP (y
))
3658 double yy
= scm_i_fraction2double (y
);
3659 double xx
= SCM_REAL_VALUE (x
);
3660 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3663 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3665 else if (SCM_FRACTIONP (x
))
3667 if (SCM_I_INUMP (y
))
3671 else if (SCM_BIGP (y
))
3675 else if (SCM_REALP (y
))
3677 double xx
= scm_i_fraction2double (x
);
3678 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3680 else if (SCM_FRACTIONP (y
))
3685 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3688 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3692 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3693 /* "Return the minium of all parameter values."
3696 scm_min (SCM x
, SCM y
)
3701 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3702 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3705 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3708 if (SCM_I_INUMP (x
))
3710 long xx
= SCM_I_INUM (x
);
3711 if (SCM_I_INUMP (y
))
3713 long yy
= SCM_I_INUM (y
);
3714 return (xx
< yy
) ? x
: y
;
3716 else if (SCM_BIGP (y
))
3718 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3719 scm_remember_upto_here_1 (y
);
3720 return (sgn
< 0) ? y
: x
;
3722 else if (SCM_REALP (y
))
3725 /* if y==NaN then "<" is false and we return NaN */
3726 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3728 else if (SCM_FRACTIONP (y
))
3731 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3734 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3736 else if (SCM_BIGP (x
))
3738 if (SCM_I_INUMP (y
))
3740 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3741 scm_remember_upto_here_1 (x
);
3742 return (sgn
< 0) ? x
: y
;
3744 else if (SCM_BIGP (y
))
3746 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3747 scm_remember_upto_here_2 (x
, y
);
3748 return (cmp
> 0) ? y
: x
;
3750 else if (SCM_REALP (y
))
3752 /* if y==NaN then xx<yy is false, so we return the NaN y */
3755 xx
= scm_i_big2dbl (x
);
3756 yy
= SCM_REAL_VALUE (y
);
3757 return (xx
< yy
? scm_from_double (xx
) : y
);
3759 else if (SCM_FRACTIONP (y
))
3764 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3766 else if (SCM_REALP (x
))
3768 if (SCM_I_INUMP (y
))
3770 double z
= SCM_I_INUM (y
);
3771 /* if x==NaN then "<" is false and we return NaN */
3772 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3774 else if (SCM_BIGP (y
))
3779 else if (SCM_REALP (y
))
3781 /* if x==NaN then our explicit check means we return NaN
3782 if y==NaN then "<" is false and we return NaN
3783 calling isnan is unavoidable, since it's the only way to know
3784 which of x or y causes any compares to be false */
3785 double xx
= SCM_REAL_VALUE (x
);
3786 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3788 else if (SCM_FRACTIONP (y
))
3790 double yy
= scm_i_fraction2double (y
);
3791 double xx
= SCM_REAL_VALUE (x
);
3792 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3795 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3797 else if (SCM_FRACTIONP (x
))
3799 if (SCM_I_INUMP (y
))
3803 else if (SCM_BIGP (y
))
3807 else if (SCM_REALP (y
))
3809 double xx
= scm_i_fraction2double (x
);
3810 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3812 else if (SCM_FRACTIONP (y
))
3817 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3820 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3824 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3825 /* "Return the sum of all parameter values. Return 0 if called without\n"
3829 scm_sum (SCM x
, SCM y
)
3833 if (SCM_NUMBERP (x
)) return x
;
3834 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3835 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3838 if (SCM_I_INUMP (x
))
3840 if (SCM_I_INUMP (y
))
3842 long xx
= SCM_I_INUM (x
);
3843 long yy
= SCM_I_INUM (y
);
3844 long int z
= xx
+ yy
;
3845 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3847 else if (SCM_BIGP (y
))
3852 else if (SCM_REALP (y
))
3854 long int xx
= SCM_I_INUM (x
);
3855 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3857 else if (SCM_COMPLEXP (y
))
3859 long int xx
= SCM_I_INUM (x
);
3860 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3861 SCM_COMPLEX_IMAG (y
));
3863 else if (SCM_FRACTIONP (y
))
3864 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3865 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3866 SCM_FRACTION_DENOMINATOR (y
));
3868 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3869 } else if (SCM_BIGP (x
))
3871 if (SCM_I_INUMP (y
))
3876 inum
= SCM_I_INUM (y
);
3879 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3882 SCM result
= scm_i_mkbig ();
3883 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3884 scm_remember_upto_here_1 (x
);
3885 /* we know the result will have to be a bignum */
3888 return scm_i_normbig (result
);
3892 SCM result
= scm_i_mkbig ();
3893 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3894 scm_remember_upto_here_1 (x
);
3895 /* we know the result will have to be a bignum */
3898 return scm_i_normbig (result
);
3901 else if (SCM_BIGP (y
))
3903 SCM result
= scm_i_mkbig ();
3904 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3905 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3906 mpz_add (SCM_I_BIG_MPZ (result
),
3909 scm_remember_upto_here_2 (x
, y
);
3910 /* we know the result will have to be a bignum */
3913 return scm_i_normbig (result
);
3915 else if (SCM_REALP (y
))
3917 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3918 scm_remember_upto_here_1 (x
);
3919 return scm_from_double (result
);
3921 else if (SCM_COMPLEXP (y
))
3923 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3924 + SCM_COMPLEX_REAL (y
));
3925 scm_remember_upto_here_1 (x
);
3926 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
3928 else if (SCM_FRACTIONP (y
))
3929 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3930 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3931 SCM_FRACTION_DENOMINATOR (y
));
3933 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3935 else if (SCM_REALP (x
))
3937 if (SCM_I_INUMP (y
))
3938 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
3939 else if (SCM_BIGP (y
))
3941 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3942 scm_remember_upto_here_1 (y
);
3943 return scm_from_double (result
);
3945 else if (SCM_REALP (y
))
3946 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3947 else if (SCM_COMPLEXP (y
))
3948 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3949 SCM_COMPLEX_IMAG (y
));
3950 else if (SCM_FRACTIONP (y
))
3951 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3953 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3955 else if (SCM_COMPLEXP (x
))
3957 if (SCM_I_INUMP (y
))
3958 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
3959 SCM_COMPLEX_IMAG (x
));
3960 else if (SCM_BIGP (y
))
3962 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3963 + SCM_COMPLEX_REAL (x
));
3964 scm_remember_upto_here_1 (y
);
3965 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
3967 else if (SCM_REALP (y
))
3968 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3969 SCM_COMPLEX_IMAG (x
));
3970 else if (SCM_COMPLEXP (y
))
3971 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3972 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3973 else if (SCM_FRACTIONP (y
))
3974 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3975 SCM_COMPLEX_IMAG (x
));
3977 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3979 else if (SCM_FRACTIONP (x
))
3981 if (SCM_I_INUMP (y
))
3982 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3983 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3984 SCM_FRACTION_DENOMINATOR (x
));
3985 else if (SCM_BIGP (y
))
3986 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3987 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3988 SCM_FRACTION_DENOMINATOR (x
));
3989 else if (SCM_REALP (y
))
3990 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
3991 else if (SCM_COMPLEXP (y
))
3992 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
3993 SCM_COMPLEX_IMAG (y
));
3994 else if (SCM_FRACTIONP (y
))
3995 /* a/b + c/d = (ad + bc) / bd */
3996 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
3997 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
3998 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4000 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4003 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4007 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4008 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4009 * the sum of all but the first argument are subtracted from the first
4011 #define FUNC_NAME s_difference
4013 scm_difference (SCM x
, SCM y
)
4018 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4020 if (SCM_I_INUMP (x
))
4022 long xx
= -SCM_I_INUM (x
);
4023 if (SCM_FIXABLE (xx
))
4024 return SCM_I_MAKINUM (xx
);
4026 return scm_i_long2big (xx
);
4028 else if (SCM_BIGP (x
))
4029 /* FIXME: do we really need to normalize here? */
4030 return scm_i_normbig (scm_i_clonebig (x
, 0));
4031 else if (SCM_REALP (x
))
4032 return scm_from_double (-SCM_REAL_VALUE (x
));
4033 else if (SCM_COMPLEXP (x
))
4034 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4035 -SCM_COMPLEX_IMAG (x
));
4036 else if (SCM_FRACTIONP (x
))
4037 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4038 SCM_FRACTION_DENOMINATOR (x
));
4040 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4043 if (SCM_I_INUMP (x
))
4045 if (SCM_I_INUMP (y
))
4047 long int xx
= SCM_I_INUM (x
);
4048 long int yy
= SCM_I_INUM (y
);
4049 long int z
= xx
- yy
;
4050 if (SCM_FIXABLE (z
))
4051 return SCM_I_MAKINUM (z
);
4053 return scm_i_long2big (z
);
4055 else if (SCM_BIGP (y
))
4057 /* inum-x - big-y */
4058 long xx
= SCM_I_INUM (x
);
4061 return scm_i_clonebig (y
, 0);
4064 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4065 SCM result
= scm_i_mkbig ();
4068 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4071 /* x - y == -(y + -x) */
4072 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4073 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4075 scm_remember_upto_here_1 (y
);
4077 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4078 /* we know the result will have to be a bignum */
4081 return scm_i_normbig (result
);
4084 else if (SCM_REALP (y
))
4086 long int xx
= SCM_I_INUM (x
);
4087 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4089 else if (SCM_COMPLEXP (y
))
4091 long int xx
= SCM_I_INUM (x
);
4092 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4093 - SCM_COMPLEX_IMAG (y
));
4095 else if (SCM_FRACTIONP (y
))
4096 /* a - b/c = (ac - b) / c */
4097 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4098 SCM_FRACTION_NUMERATOR (y
)),
4099 SCM_FRACTION_DENOMINATOR (y
));
4101 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4103 else if (SCM_BIGP (x
))
4105 if (SCM_I_INUMP (y
))
4107 /* big-x - inum-y */
4108 long yy
= SCM_I_INUM (y
);
4109 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4111 scm_remember_upto_here_1 (x
);
4113 return (SCM_FIXABLE (-yy
) ?
4114 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4117 SCM result
= scm_i_mkbig ();
4120 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4122 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4123 scm_remember_upto_here_1 (x
);
4125 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4126 /* we know the result will have to be a bignum */
4129 return scm_i_normbig (result
);
4132 else if (SCM_BIGP (y
))
4134 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4135 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4136 SCM result
= scm_i_mkbig ();
4137 mpz_sub (SCM_I_BIG_MPZ (result
),
4140 scm_remember_upto_here_2 (x
, y
);
4141 /* we know the result will have to be a bignum */
4142 if ((sgn_x
== 1) && (sgn_y
== -1))
4144 if ((sgn_x
== -1) && (sgn_y
== 1))
4146 return scm_i_normbig (result
);
4148 else if (SCM_REALP (y
))
4150 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4151 scm_remember_upto_here_1 (x
);
4152 return scm_from_double (result
);
4154 else if (SCM_COMPLEXP (y
))
4156 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4157 - SCM_COMPLEX_REAL (y
));
4158 scm_remember_upto_here_1 (x
);
4159 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4161 else if (SCM_FRACTIONP (y
))
4162 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4163 SCM_FRACTION_NUMERATOR (y
)),
4164 SCM_FRACTION_DENOMINATOR (y
));
4165 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4167 else if (SCM_REALP (x
))
4169 if (SCM_I_INUMP (y
))
4170 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4171 else if (SCM_BIGP (y
))
4173 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4174 scm_remember_upto_here_1 (x
);
4175 return scm_from_double (result
);
4177 else if (SCM_REALP (y
))
4178 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4179 else if (SCM_COMPLEXP (y
))
4180 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4181 -SCM_COMPLEX_IMAG (y
));
4182 else if (SCM_FRACTIONP (y
))
4183 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4185 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4187 else if (SCM_COMPLEXP (x
))
4189 if (SCM_I_INUMP (y
))
4190 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4191 SCM_COMPLEX_IMAG (x
));
4192 else if (SCM_BIGP (y
))
4194 double real_part
= (SCM_COMPLEX_REAL (x
)
4195 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4196 scm_remember_upto_here_1 (x
);
4197 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4199 else if (SCM_REALP (y
))
4200 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4201 SCM_COMPLEX_IMAG (x
));
4202 else if (SCM_COMPLEXP (y
))
4203 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4204 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4205 else if (SCM_FRACTIONP (y
))
4206 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4207 SCM_COMPLEX_IMAG (x
));
4209 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4211 else if (SCM_FRACTIONP (x
))
4213 if (SCM_I_INUMP (y
))
4214 /* a/b - c = (a - cb) / b */
4215 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4216 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4217 SCM_FRACTION_DENOMINATOR (x
));
4218 else if (SCM_BIGP (y
))
4219 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4220 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4221 SCM_FRACTION_DENOMINATOR (x
));
4222 else if (SCM_REALP (y
))
4223 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4224 else if (SCM_COMPLEXP (y
))
4225 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4226 -SCM_COMPLEX_IMAG (y
));
4227 else if (SCM_FRACTIONP (y
))
4228 /* a/b - c/d = (ad - bc) / bd */
4229 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4230 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4231 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4233 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4236 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4241 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4242 /* "Return the product of all arguments. If called without arguments,\n"
4246 scm_product (SCM x
, SCM y
)
4251 return SCM_I_MAKINUM (1L);
4252 else if (SCM_NUMBERP (x
))
4255 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4258 if (SCM_I_INUMP (x
))
4263 xx
= SCM_I_INUM (x
);
4267 case 0: return x
; break;
4268 case 1: return y
; break;
4271 if (SCM_I_INUMP (y
))
4273 long yy
= SCM_I_INUM (y
);
4275 SCM k
= SCM_I_MAKINUM (kk
);
4276 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4280 SCM result
= scm_i_long2big (xx
);
4281 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4282 return scm_i_normbig (result
);
4285 else if (SCM_BIGP (y
))
4287 SCM result
= scm_i_mkbig ();
4288 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4289 scm_remember_upto_here_1 (y
);
4292 else if (SCM_REALP (y
))
4293 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4294 else if (SCM_COMPLEXP (y
))
4295 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4296 xx
* SCM_COMPLEX_IMAG (y
));
4297 else if (SCM_FRACTIONP (y
))
4298 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4299 SCM_FRACTION_DENOMINATOR (y
));
4301 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4303 else if (SCM_BIGP (x
))
4305 if (SCM_I_INUMP (y
))
4310 else if (SCM_BIGP (y
))
4312 SCM result
= scm_i_mkbig ();
4313 mpz_mul (SCM_I_BIG_MPZ (result
),
4316 scm_remember_upto_here_2 (x
, y
);
4319 else if (SCM_REALP (y
))
4321 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4322 scm_remember_upto_here_1 (x
);
4323 return scm_from_double (result
);
4325 else if (SCM_COMPLEXP (y
))
4327 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4328 scm_remember_upto_here_1 (x
);
4329 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4330 z
* SCM_COMPLEX_IMAG (y
));
4332 else if (SCM_FRACTIONP (y
))
4333 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4334 SCM_FRACTION_DENOMINATOR (y
));
4336 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4338 else if (SCM_REALP (x
))
4340 if (SCM_I_INUMP (y
))
4341 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4342 else if (SCM_BIGP (y
))
4344 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4345 scm_remember_upto_here_1 (y
);
4346 return scm_from_double (result
);
4348 else if (SCM_REALP (y
))
4349 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4350 else if (SCM_COMPLEXP (y
))
4351 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4352 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4353 else if (SCM_FRACTIONP (y
))
4354 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4356 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4358 else if (SCM_COMPLEXP (x
))
4360 if (SCM_I_INUMP (y
))
4361 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4362 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4363 else if (SCM_BIGP (y
))
4365 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4366 scm_remember_upto_here_1 (y
);
4367 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4368 z
* SCM_COMPLEX_IMAG (x
));
4370 else if (SCM_REALP (y
))
4371 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4372 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4373 else if (SCM_COMPLEXP (y
))
4375 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4376 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4377 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4378 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4380 else if (SCM_FRACTIONP (y
))
4382 double yy
= scm_i_fraction2double (y
);
4383 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4384 yy
* SCM_COMPLEX_IMAG (x
));
4387 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4389 else if (SCM_FRACTIONP (x
))
4391 if (SCM_I_INUMP (y
))
4392 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4393 SCM_FRACTION_DENOMINATOR (x
));
4394 else if (SCM_BIGP (y
))
4395 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4396 SCM_FRACTION_DENOMINATOR (x
));
4397 else if (SCM_REALP (y
))
4398 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4399 else if (SCM_COMPLEXP (y
))
4401 double xx
= scm_i_fraction2double (x
);
4402 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4403 xx
* SCM_COMPLEX_IMAG (y
));
4405 else if (SCM_FRACTIONP (y
))
4406 /* a/b * c/d = ac / bd */
4407 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4408 SCM_FRACTION_NUMERATOR (y
)),
4409 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4410 SCM_FRACTION_DENOMINATOR (y
)));
4412 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4415 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4418 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4419 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4420 #define ALLOW_DIVIDE_BY_ZERO
4421 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4424 /* The code below for complex division is adapted from the GNU
4425 libstdc++, which adapted it from f2c's libF77, and is subject to
4428 /****************************************************************
4429 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4431 Permission to use, copy, modify, and distribute this software
4432 and its documentation for any purpose and without fee is hereby
4433 granted, provided that the above copyright notice appear in all
4434 copies and that both that the copyright notice and this
4435 permission notice and warranty disclaimer appear in supporting
4436 documentation, and that the names of AT&T Bell Laboratories or
4437 Bellcore or any of their entities not be used in advertising or
4438 publicity pertaining to distribution of the software without
4439 specific, written prior permission.
4441 AT&T and Bellcore disclaim all warranties with regard to this
4442 software, including all implied warranties of merchantability
4443 and fitness. In no event shall AT&T or Bellcore be liable for
4444 any special, indirect or consequential damages or any damages
4445 whatsoever resulting from loss of use, data or profits, whether
4446 in an action of contract, negligence or other tortious action,
4447 arising out of or in connection with the use or performance of
4449 ****************************************************************/
4451 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4452 /* Divide the first argument by the product of the remaining
4453 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4455 #define FUNC_NAME s_divide
4457 scm_i_divide (SCM x
, SCM y
, int inexact
)
4464 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4465 else if (SCM_I_INUMP (x
))
4467 long xx
= SCM_I_INUM (x
);
4468 if (xx
== 1 || xx
== -1)
4470 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4472 scm_num_overflow (s_divide
);
4477 return scm_from_double (1.0 / (double) xx
);
4478 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4481 else if (SCM_BIGP (x
))
4484 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4485 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4487 else if (SCM_REALP (x
))
4489 double xx
= SCM_REAL_VALUE (x
);
4490 #ifndef ALLOW_DIVIDE_BY_ZERO
4492 scm_num_overflow (s_divide
);
4495 return scm_from_double (1.0 / xx
);
4497 else if (SCM_COMPLEXP (x
))
4499 double r
= SCM_COMPLEX_REAL (x
);
4500 double i
= SCM_COMPLEX_IMAG (x
);
4504 double d
= i
* (1.0 + t
* t
);
4505 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4510 double d
= r
* (1.0 + t
* t
);
4511 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4514 else if (SCM_FRACTIONP (x
))
4515 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4516 SCM_FRACTION_NUMERATOR (x
));
4518 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4521 if (SCM_I_INUMP (x
))
4523 long xx
= SCM_I_INUM (x
);
4524 if (SCM_I_INUMP (y
))
4526 long yy
= SCM_I_INUM (y
);
4529 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4530 scm_num_overflow (s_divide
);
4532 return scm_from_double ((double) xx
/ (double) yy
);
4535 else if (xx
% yy
!= 0)
4538 return scm_from_double ((double) xx
/ (double) yy
);
4539 else return scm_i_make_ratio (x
, y
);
4544 if (SCM_FIXABLE (z
))
4545 return SCM_I_MAKINUM (z
);
4547 return scm_i_long2big (z
);
4550 else if (SCM_BIGP (y
))
4553 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4554 else return scm_i_make_ratio (x
, y
);
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_from_double ((double) xx
/ yy
);
4566 else if (SCM_COMPLEXP (y
))
4569 complex_div
: /* y _must_ be a complex number */
4571 double r
= SCM_COMPLEX_REAL (y
);
4572 double i
= SCM_COMPLEX_IMAG (y
);
4576 double d
= i
* (1.0 + t
* t
);
4577 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4582 double d
= r
* (1.0 + t
* t
);
4583 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4587 else if (SCM_FRACTIONP (y
))
4588 /* a / b/c = ac / b */
4589 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4590 SCM_FRACTION_NUMERATOR (y
));
4592 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4594 else if (SCM_BIGP (x
))
4596 if (SCM_I_INUMP (y
))
4598 long int yy
= SCM_I_INUM (y
);
4601 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4602 scm_num_overflow (s_divide
);
4604 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4605 scm_remember_upto_here_1 (x
);
4606 return (sgn
== 0) ? scm_nan () : scm_inf ();
4613 /* FIXME: HMM, what are the relative performance issues here?
4614 We need to test. Is it faster on average to test
4615 divisible_p, then perform whichever operation, or is it
4616 faster to perform the integer div opportunistically and
4617 switch to real if there's a remainder? For now we take the
4618 middle ground: test, then if divisible, use the faster div
4621 long abs_yy
= yy
< 0 ? -yy
: yy
;
4622 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4626 SCM result
= scm_i_mkbig ();
4627 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4628 scm_remember_upto_here_1 (x
);
4630 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4631 return scm_i_normbig (result
);
4636 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4637 else return scm_i_make_ratio (x
, y
);
4641 else if (SCM_BIGP (y
))
4643 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4646 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4647 scm_num_overflow (s_divide
);
4649 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4650 scm_remember_upto_here_1 (x
);
4651 return (sgn
== 0) ? scm_nan () : scm_inf ();
4657 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4661 SCM result
= scm_i_mkbig ();
4662 mpz_divexact (SCM_I_BIG_MPZ (result
),
4665 scm_remember_upto_here_2 (x
, y
);
4666 return scm_i_normbig (result
);
4672 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4673 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4674 scm_remember_upto_here_2 (x
, y
);
4675 return scm_from_double (dbx
/ dby
);
4677 else return scm_i_make_ratio (x
, y
);
4681 else if (SCM_REALP (y
))
4683 double yy
= SCM_REAL_VALUE (y
);
4684 #ifndef ALLOW_DIVIDE_BY_ZERO
4686 scm_num_overflow (s_divide
);
4689 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4691 else if (SCM_COMPLEXP (y
))
4693 a
= scm_i_big2dbl (x
);
4696 else if (SCM_FRACTIONP (y
))
4697 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4698 SCM_FRACTION_NUMERATOR (y
));
4700 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4702 else if (SCM_REALP (x
))
4704 double rx
= SCM_REAL_VALUE (x
);
4705 if (SCM_I_INUMP (y
))
4707 long int yy
= SCM_I_INUM (y
);
4708 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4710 scm_num_overflow (s_divide
);
4713 return scm_from_double (rx
/ (double) yy
);
4715 else if (SCM_BIGP (y
))
4717 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4718 scm_remember_upto_here_1 (y
);
4719 return scm_from_double (rx
/ dby
);
4721 else if (SCM_REALP (y
))
4723 double yy
= SCM_REAL_VALUE (y
);
4724 #ifndef ALLOW_DIVIDE_BY_ZERO
4726 scm_num_overflow (s_divide
);
4729 return scm_from_double (rx
/ yy
);
4731 else if (SCM_COMPLEXP (y
))
4736 else if (SCM_FRACTIONP (y
))
4737 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4739 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4741 else if (SCM_COMPLEXP (x
))
4743 double rx
= SCM_COMPLEX_REAL (x
);
4744 double ix
= SCM_COMPLEX_IMAG (x
);
4745 if (SCM_I_INUMP (y
))
4747 long int yy
= SCM_I_INUM (y
);
4748 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4750 scm_num_overflow (s_divide
);
4755 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4758 else if (SCM_BIGP (y
))
4760 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4761 scm_remember_upto_here_1 (y
);
4762 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4764 else if (SCM_REALP (y
))
4766 double yy
= SCM_REAL_VALUE (y
);
4767 #ifndef ALLOW_DIVIDE_BY_ZERO
4769 scm_num_overflow (s_divide
);
4772 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4774 else if (SCM_COMPLEXP (y
))
4776 double ry
= SCM_COMPLEX_REAL (y
);
4777 double iy
= SCM_COMPLEX_IMAG (y
);
4781 double d
= iy
* (1.0 + t
* t
);
4782 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4787 double d
= ry
* (1.0 + t
* t
);
4788 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4791 else if (SCM_FRACTIONP (y
))
4793 double yy
= scm_i_fraction2double (y
);
4794 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4797 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4799 else if (SCM_FRACTIONP (x
))
4801 if (SCM_I_INUMP (y
))
4803 long int yy
= SCM_I_INUM (y
);
4804 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4806 scm_num_overflow (s_divide
);
4809 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4810 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4812 else if (SCM_BIGP (y
))
4814 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4815 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4817 else if (SCM_REALP (y
))
4819 double yy
= SCM_REAL_VALUE (y
);
4820 #ifndef ALLOW_DIVIDE_BY_ZERO
4822 scm_num_overflow (s_divide
);
4825 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4827 else if (SCM_COMPLEXP (y
))
4829 a
= scm_i_fraction2double (x
);
4832 else if (SCM_FRACTIONP (y
))
4833 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4834 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4836 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4839 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4843 scm_divide (SCM x
, SCM y
)
4845 return scm_i_divide (x
, y
, 0);
4848 static SCM
scm_divide2real (SCM x
, SCM y
)
4850 return scm_i_divide (x
, y
, 1);
4856 scm_asinh (double x
)
4861 #define asinh scm_asinh
4862 return log (x
+ sqrt (x
* x
+ 1));
4865 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4866 /* "Return the inverse hyperbolic sine of @var{x}."
4871 scm_acosh (double x
)
4876 #define acosh scm_acosh
4877 return log (x
+ sqrt (x
* x
- 1));
4880 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4881 /* "Return the inverse hyperbolic cosine of @var{x}."
4886 scm_atanh (double x
)
4891 #define atanh scm_atanh
4892 return 0.5 * log ((1 + x
) / (1 - x
));
4895 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4896 /* "Return the inverse hyperbolic tangent of @var{x}."
4901 scm_c_truncate (double x
)
4912 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
4913 half-way case (ie. when x is an integer plus 0.5) going upwards.
4914 Then half-way cases are identified and adjusted down if the
4915 round-upwards didn't give the desired even integer.
4917 "plus_half == result" identifies a half-way case. If plus_half, which is
4918 x + 0.5, is an integer then x must be an integer plus 0.5.
4920 An odd "result" value is identified with result/2 != floor(result/2).
4921 This is done with plus_half, since that value is ready for use sooner in
4922 a pipelined cpu, and we're already requiring plus_half == result.
4924 Note however that we need to be careful when x is big and already an
4925 integer. In that case "x+0.5" may round to an adjacent integer, causing
4926 us to return such a value, incorrectly. For instance if the hardware is
4927 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
4928 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
4929 returned. Or if the hardware is in round-upwards mode, then other bigger
4930 values like say x == 2^128 will see x+0.5 rounding up to the next higher
4931 representable value, 2^128+2^76 (or whatever), again incorrect.
4933 These bad roundings of x+0.5 are avoided by testing at the start whether
4934 x is already an integer. If it is then clearly that's the desired result
4935 already. And if it's not then the exponent must be small enough to allow
4936 an 0.5 to be represented, and hence added without a bad rounding. */
4939 scm_c_round (double x
)
4941 double plus_half
, result
;
4946 plus_half
= x
+ 0.5;
4947 result
= floor (plus_half
);
4948 /* Adjust so that the rounding is towards even. */
4949 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4954 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4956 "Round the number @var{x} towards zero.")
4957 #define FUNC_NAME s_scm_truncate_number
4959 if (scm_is_false (scm_negative_p (x
)))
4960 return scm_floor (x
);
4962 return scm_ceiling (x
);
4966 static SCM exactly_one_half
;
4968 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4970 "Round the number @var{x} towards the nearest integer. "
4971 "When it is exactly halfway between two integers, "
4972 "round towards the even one.")
4973 #define FUNC_NAME s_scm_round_number
4975 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4977 else if (SCM_REALP (x
))
4978 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
4981 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
4982 single quotient+remainder division then examining to see which way
4983 the rounding should go. */
4984 SCM plus_half
= scm_sum (x
, exactly_one_half
);
4985 SCM result
= scm_floor (plus_half
);
4986 /* Adjust so that the rounding is towards even. */
4987 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
4988 && scm_is_true (scm_odd_p (result
)))
4989 return scm_difference (result
, SCM_I_MAKINUM (1));
4996 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
4998 "Round the number @var{x} towards minus infinity.")
4999 #define FUNC_NAME s_scm_floor
5001 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5003 else if (SCM_REALP (x
))
5004 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5005 else if (SCM_FRACTIONP (x
))
5007 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5008 SCM_FRACTION_DENOMINATOR (x
));
5009 if (scm_is_false (scm_negative_p (x
)))
5011 /* For positive x, rounding towards zero is correct. */
5016 /* For negative x, we need to return q-1 unless x is an
5017 integer. But fractions are never integer, per our
5019 return scm_difference (q
, SCM_I_MAKINUM (1));
5023 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5027 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5029 "Round the number @var{x} towards infinity.")
5030 #define FUNC_NAME s_scm_ceiling
5032 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5034 else if (SCM_REALP (x
))
5035 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5036 else if (SCM_FRACTIONP (x
))
5038 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5039 SCM_FRACTION_DENOMINATOR (x
));
5040 if (scm_is_false (scm_positive_p (x
)))
5042 /* For negative x, rounding towards zero is correct. */
5047 /* For positive x, we need to return q+1 unless x is an
5048 integer. But fractions are never integer, per our
5050 return scm_sum (q
, SCM_I_MAKINUM (1));
5054 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5058 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5059 /* "Return the square root of the real number @var{x}."
5061 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5062 /* "Return the absolute value of the real number @var{x}."
5064 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5065 /* "Return the @var{x}th power of e."
5067 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5068 /* "Return the natural logarithm of the real number @var{x}."
5070 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5071 /* "Return the sine of the real number @var{x}."
5073 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5074 /* "Return the cosine of the real number @var{x}."
5076 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5077 /* "Return the tangent of the real number @var{x}."
5079 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5080 /* "Return the arc sine of the real number @var{x}."
5082 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5083 /* "Return the arc cosine of the real number @var{x}."
5085 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5086 /* "Return the arc tangent of the real number @var{x}."
5088 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5089 /* "Return the hyperbolic sine of the real number @var{x}."
5091 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5092 /* "Return the hyperbolic cosine of the real number @var{x}."
5094 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5095 /* "Return the hyperbolic tangent of the real number @var{x}."
5103 static void scm_two_doubles (SCM x
,
5105 const char *sstring
,
5109 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5111 if (SCM_I_INUMP (x
))
5112 xy
->x
= SCM_I_INUM (x
);
5113 else if (SCM_BIGP (x
))
5114 xy
->x
= scm_i_big2dbl (x
);
5115 else if (SCM_REALP (x
))
5116 xy
->x
= SCM_REAL_VALUE (x
);
5117 else if (SCM_FRACTIONP (x
))
5118 xy
->x
= scm_i_fraction2double (x
);
5120 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5122 if (SCM_I_INUMP (y
))
5123 xy
->y
= SCM_I_INUM (y
);
5124 else if (SCM_BIGP (y
))
5125 xy
->y
= scm_i_big2dbl (y
);
5126 else if (SCM_REALP (y
))
5127 xy
->y
= SCM_REAL_VALUE (y
);
5128 else if (SCM_FRACTIONP (y
))
5129 xy
->y
= scm_i_fraction2double (y
);
5131 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5135 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5137 "Return @var{x} raised to the power of @var{y}. This\n"
5138 "procedure does not accept complex arguments.")
5139 #define FUNC_NAME s_scm_sys_expt
5142 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5143 return scm_from_double (pow (xy
.x
, xy
.y
));
5148 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5150 "Return the arc tangent of the two arguments @var{x} and\n"
5151 "@var{y}. This is similar to calculating the arc tangent of\n"
5152 "@var{x} / @var{y}, except that the signs of both arguments\n"
5153 "are used to determine the quadrant of the result. This\n"
5154 "procedure does not accept complex arguments.")
5155 #define FUNC_NAME s_scm_sys_atan2
5158 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5159 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5164 scm_c_make_rectangular (double re
, double im
)
5167 return scm_from_double (re
);
5171 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5173 SCM_COMPLEX_REAL (z
) = re
;
5174 SCM_COMPLEX_IMAG (z
) = im
;
5179 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5180 (SCM real
, SCM imaginary
),
5181 "Return a complex number constructed of the given @var{real} and\n"
5182 "@var{imaginary} parts.")
5183 #define FUNC_NAME s_scm_make_rectangular
5186 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5187 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5192 scm_c_make_polar (double mag
, double ang
)
5196 sincos (ang
, &s
, &c
);
5201 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5204 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5206 "Return the complex number @var{x} * e^(i * @var{y}).")
5207 #define FUNC_NAME s_scm_make_polar
5210 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5211 return scm_c_make_polar (xy
.x
, xy
.y
);
5216 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5217 /* "Return the real part of the number @var{z}."
5220 scm_real_part (SCM z
)
5222 if (SCM_I_INUMP (z
))
5224 else if (SCM_BIGP (z
))
5226 else if (SCM_REALP (z
))
5228 else if (SCM_COMPLEXP (z
))
5229 return scm_from_double (SCM_COMPLEX_REAL (z
));
5230 else if (SCM_FRACTIONP (z
))
5233 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5237 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5238 /* "Return the imaginary part of the number @var{z}."
5241 scm_imag_part (SCM z
)
5243 if (SCM_I_INUMP (z
))
5245 else if (SCM_BIGP (z
))
5247 else if (SCM_REALP (z
))
5249 else if (SCM_COMPLEXP (z
))
5250 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5251 else if (SCM_FRACTIONP (z
))
5254 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5257 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5258 /* "Return the numerator of the number @var{z}."
5261 scm_numerator (SCM z
)
5263 if (SCM_I_INUMP (z
))
5265 else if (SCM_BIGP (z
))
5267 else if (SCM_FRACTIONP (z
))
5269 scm_i_fraction_reduce (z
);
5270 return SCM_FRACTION_NUMERATOR (z
);
5272 else if (SCM_REALP (z
))
5273 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5275 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5279 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5280 /* "Return the denominator of the number @var{z}."
5283 scm_denominator (SCM z
)
5285 if (SCM_I_INUMP (z
))
5286 return SCM_I_MAKINUM (1);
5287 else if (SCM_BIGP (z
))
5288 return SCM_I_MAKINUM (1);
5289 else if (SCM_FRACTIONP (z
))
5291 scm_i_fraction_reduce (z
);
5292 return SCM_FRACTION_DENOMINATOR (z
);
5294 else if (SCM_REALP (z
))
5295 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5297 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5300 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5301 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5302 * "@code{abs} for real arguments, but also allows complex numbers."
5305 scm_magnitude (SCM z
)
5307 if (SCM_I_INUMP (z
))
5309 long int zz
= SCM_I_INUM (z
);
5312 else if (SCM_POSFIXABLE (-zz
))
5313 return SCM_I_MAKINUM (-zz
);
5315 return scm_i_long2big (-zz
);
5317 else if (SCM_BIGP (z
))
5319 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5320 scm_remember_upto_here_1 (z
);
5322 return scm_i_clonebig (z
, 0);
5326 else if (SCM_REALP (z
))
5327 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5328 else if (SCM_COMPLEXP (z
))
5329 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5330 else if (SCM_FRACTIONP (z
))
5332 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5334 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5335 SCM_FRACTION_DENOMINATOR (z
));
5338 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5342 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5343 /* "Return the angle of the complex number @var{z}."
5348 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5349 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5350 But if atan2 follows the floating point rounding mode, then the value
5351 is not a constant. Maybe it'd be close enough though. */
5352 if (SCM_I_INUMP (z
))
5354 if (SCM_I_INUM (z
) >= 0)
5357 return scm_from_double (atan2 (0.0, -1.0));
5359 else if (SCM_BIGP (z
))
5361 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5362 scm_remember_upto_here_1 (z
);
5364 return scm_from_double (atan2 (0.0, -1.0));
5368 else if (SCM_REALP (z
))
5370 if (SCM_REAL_VALUE (z
) >= 0)
5373 return scm_from_double (atan2 (0.0, -1.0));
5375 else if (SCM_COMPLEXP (z
))
5376 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5377 else if (SCM_FRACTIONP (z
))
5379 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5381 else return scm_from_double (atan2 (0.0, -1.0));
5384 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5388 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5389 /* Convert the number @var{x} to its inexact representation.\n"
5392 scm_exact_to_inexact (SCM z
)
5394 if (SCM_I_INUMP (z
))
5395 return scm_from_double ((double) SCM_I_INUM (z
));
5396 else if (SCM_BIGP (z
))
5397 return scm_from_double (scm_i_big2dbl (z
));
5398 else if (SCM_FRACTIONP (z
))
5399 return scm_from_double (scm_i_fraction2double (z
));
5400 else if (SCM_INEXACTP (z
))
5403 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5407 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5409 "Return an exact number that is numerically closest to @var{z}.")
5410 #define FUNC_NAME s_scm_inexact_to_exact
5412 if (SCM_I_INUMP (z
))
5414 else if (SCM_BIGP (z
))
5416 else if (SCM_REALP (z
))
5418 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5419 SCM_OUT_OF_RANGE (1, z
);
5426 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5427 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5428 scm_i_mpz2num (mpq_denref (frac
)));
5430 /* When scm_i_make_ratio throws, we leak the memory allocated
5437 else if (SCM_FRACTIONP (z
))
5440 SCM_WRONG_TYPE_ARG (1, z
);
5444 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5446 "Return an exact number that is within @var{err} of @var{x}.")
5447 #define FUNC_NAME s_scm_rationalize
5449 if (SCM_I_INUMP (x
))
5451 else if (SCM_BIGP (x
))
5453 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5455 /* Use continued fractions to find closest ratio. All
5456 arithmetic is done with exact numbers.
5459 SCM ex
= scm_inexact_to_exact (x
);
5460 SCM int_part
= scm_floor (ex
);
5461 SCM tt
= SCM_I_MAKINUM (1);
5462 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5463 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5467 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5470 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5471 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5473 /* We stop after a million iterations just to be absolutely sure
5474 that we don't go into an infinite loop. The process normally
5475 converges after less than a dozen iterations.
5478 err
= scm_abs (err
);
5479 while (++i
< 1000000)
5481 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5482 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5483 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5485 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5486 err
))) /* abs(x-a/b) <= err */
5488 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5489 if (scm_is_false (scm_exact_p (x
))
5490 || scm_is_false (scm_exact_p (err
)))
5491 return scm_exact_to_inexact (res
);
5495 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5497 tt
= scm_floor (rx
); /* tt = floor (rx) */
5503 scm_num_overflow (s_scm_rationalize
);
5506 SCM_WRONG_TYPE_ARG (1, x
);
5510 /* conversion functions */
5513 scm_is_integer (SCM val
)
5515 return scm_is_true (scm_integer_p (val
));
5519 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5521 if (SCM_I_INUMP (val
))
5523 scm_t_signed_bits n
= SCM_I_INUM (val
);
5524 return n
>= min
&& n
<= max
;
5526 else if (SCM_BIGP (val
))
5528 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5530 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5532 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5534 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5535 return n
>= min
&& n
<= max
;
5545 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5546 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5549 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5550 SCM_I_BIG_MPZ (val
));
5552 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5564 return n
>= min
&& n
<= max
;
5572 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5574 if (SCM_I_INUMP (val
))
5576 scm_t_signed_bits n
= SCM_I_INUM (val
);
5577 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5579 else if (SCM_BIGP (val
))
5581 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5583 else if (max
<= ULONG_MAX
)
5585 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5587 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5588 return n
>= min
&& n
<= max
;
5598 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5601 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5602 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5605 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5606 SCM_I_BIG_MPZ (val
));
5608 return n
>= min
&& n
<= max
;
5615 #define TYPE scm_t_intmax
5616 #define TYPE_MIN min
5617 #define TYPE_MAX max
5618 #define SIZEOF_TYPE 0
5619 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5620 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5621 #include "libguile/conv-integer.i.c"
5623 #define TYPE scm_t_uintmax
5624 #define TYPE_MIN min
5625 #define TYPE_MAX max
5626 #define SIZEOF_TYPE 0
5627 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5628 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5629 #include "libguile/conv-uinteger.i.c"
5631 #define TYPE scm_t_int8
5632 #define TYPE_MIN SCM_T_INT8_MIN
5633 #define TYPE_MAX SCM_T_INT8_MAX
5634 #define SIZEOF_TYPE 1
5635 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5636 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5637 #include "libguile/conv-integer.i.c"
5639 #define TYPE scm_t_uint8
5641 #define TYPE_MAX SCM_T_UINT8_MAX
5642 #define SIZEOF_TYPE 1
5643 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5644 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5645 #include "libguile/conv-uinteger.i.c"
5647 #define TYPE scm_t_int16
5648 #define TYPE_MIN SCM_T_INT16_MIN
5649 #define TYPE_MAX SCM_T_INT16_MAX
5650 #define SIZEOF_TYPE 2
5651 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5652 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5653 #include "libguile/conv-integer.i.c"
5655 #define TYPE scm_t_uint16
5657 #define TYPE_MAX SCM_T_UINT16_MAX
5658 #define SIZEOF_TYPE 2
5659 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5660 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5661 #include "libguile/conv-uinteger.i.c"
5663 #define TYPE scm_t_int32
5664 #define TYPE_MIN SCM_T_INT32_MIN
5665 #define TYPE_MAX SCM_T_INT32_MAX
5666 #define SIZEOF_TYPE 4
5667 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5668 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5669 #include "libguile/conv-integer.i.c"
5671 #define TYPE scm_t_uint32
5673 #define TYPE_MAX SCM_T_UINT32_MAX
5674 #define SIZEOF_TYPE 4
5675 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5676 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5677 #include "libguile/conv-uinteger.i.c"
5679 #if SCM_HAVE_T_INT64
5681 #define TYPE scm_t_int64
5682 #define TYPE_MIN SCM_T_INT64_MIN
5683 #define TYPE_MAX SCM_T_INT64_MAX
5684 #define SIZEOF_TYPE 8
5685 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5686 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5687 #include "libguile/conv-integer.i.c"
5689 #define TYPE scm_t_uint64
5691 #define TYPE_MAX SCM_T_UINT64_MAX
5692 #define SIZEOF_TYPE 8
5693 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5694 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5695 #include "libguile/conv-uinteger.i.c"
5700 scm_is_real (SCM val
)
5702 return scm_is_true (scm_real_p (val
));
5706 scm_is_rational (SCM val
)
5708 return scm_is_true (scm_rational_p (val
));
5712 scm_to_double (SCM val
)
5714 if (SCM_I_INUMP (val
))
5715 return SCM_I_INUM (val
);
5716 else if (SCM_BIGP (val
))
5717 return scm_i_big2dbl (val
);
5718 else if (SCM_FRACTIONP (val
))
5719 return scm_i_fraction2double (val
);
5720 else if (SCM_REALP (val
))
5721 return SCM_REAL_VALUE (val
);
5723 scm_wrong_type_arg (NULL
, 0, val
);
5727 scm_from_double (double val
)
5729 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5730 SCM_REAL_VALUE (z
) = val
;
5734 #if SCM_ENABLE_DISCOURAGED == 1
5737 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5741 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5745 scm_out_of_range (NULL
, num
);
5748 return scm_to_double (num
);
5752 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5756 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5760 scm_out_of_range (NULL
, num
);
5763 return scm_to_double (num
);
5769 scm_is_complex (SCM val
)
5771 return scm_is_true (scm_complex_p (val
));
5775 scm_c_real_part (SCM z
)
5777 if (SCM_COMPLEXP (z
))
5778 return SCM_COMPLEX_REAL (z
);
5781 /* Use the scm_real_part to get proper error checking and
5784 return scm_to_double (scm_real_part (z
));
5789 scm_c_imag_part (SCM z
)
5791 if (SCM_COMPLEXP (z
))
5792 return SCM_COMPLEX_IMAG (z
);
5795 /* Use the scm_imag_part to get proper error checking and
5796 dispatching. The result will almost always be 0.0, but not
5799 return scm_to_double (scm_imag_part (z
));
5804 scm_c_magnitude (SCM z
)
5806 return scm_to_double (scm_magnitude (z
));
5812 return scm_to_double (scm_angle (z
));
5816 scm_is_number (SCM z
)
5818 return scm_is_true (scm_number_p (z
));
5826 mpz_init_set_si (z_negative_one
, -1);
5828 /* It may be possible to tune the performance of some algorithms by using
5829 * the following constants to avoid the creation of bignums. Please, before
5830 * using these values, remember the two rules of program optimization:
5831 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5832 scm_c_define ("most-positive-fixnum",
5833 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5834 scm_c_define ("most-negative-fixnum",
5835 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5837 scm_add_feature ("complex");
5838 scm_add_feature ("inexact");
5839 scm_flo0
= scm_from_double (0.0);
5841 /* determine floating point precision */
5842 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
5844 init_dblprec(&scm_dblprec
[i
-2],i
);
5845 init_fx_radix(fx_per_radix
[i
-2],i
);
5848 /* hard code precision for base 10 if the preprocessor tells us to... */
5849 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5852 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
5853 SCM_I_MAKINUM (2)));
5854 #include "libguile/numbers.x"