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{n} is infinite, @code{#f}\n"
550 #define FUNC_NAME s_scm_inf_p
553 return scm_from_bool (xisinf (SCM_REAL_VALUE (n
)));
554 else if (SCM_COMPLEXP (n
))
555 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (n
))
556 || xisinf (SCM_COMPLEX_IMAG (n
)));
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
);
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_mem2string (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_take0str (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_mem2string ("/", 1),
2302 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2304 else if (SCM_INEXACTP (n
))
2306 char num_buf
[FLOBUFLEN
];
2307 return scm_mem2string (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_STRING_CHARS (str
), SCM_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_mem2string (&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_STRING_CHARS (string
),
2971 SCM_STRING_LENGTH (string
),
2973 return scm_return_first (answer
, string
);
2978 /*** END strs->nums ***/
2982 scm_bigequal (SCM x
, SCM y
)
2984 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2985 scm_remember_upto_here_2 (x
, y
);
2986 return scm_from_bool (0 == result
);
2990 scm_real_equalp (SCM x
, SCM y
)
2992 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
2996 scm_complex_equalp (SCM x
, SCM y
)
2998 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
2999 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3003 scm_i_fraction_equalp (SCM x
, SCM y
)
3005 scm_i_fraction_reduce (x
);
3006 scm_i_fraction_reduce (y
);
3007 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3008 SCM_FRACTION_NUMERATOR (y
)))
3009 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3010 SCM_FRACTION_DENOMINATOR (y
))))
3017 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3019 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3021 #define FUNC_NAME s_scm_number_p
3023 return scm_from_bool (SCM_NUMBERP (x
));
3027 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3029 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3030 "otherwise. Note that the sets of real, rational and integer\n"
3031 "values form subsets of the set of complex numbers, i. e. the\n"
3032 "predicate will also be fulfilled if @var{x} is a real,\n"
3033 "rational or integer number.")
3034 #define FUNC_NAME s_scm_complex_p
3036 /* all numbers are complex. */
3037 return scm_number_p (x
);
3041 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3043 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3044 "otherwise. Note that the set of integer values forms a subset of\n"
3045 "the set of real numbers, i. e. the predicate will also be\n"
3046 "fulfilled if @var{x} is an integer number.")
3047 #define FUNC_NAME s_scm_real_p
3049 /* we can't represent irrational numbers. */
3050 return scm_rational_p (x
);
3054 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3056 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3057 "otherwise. Note that the set of integer values forms a subset of\n"
3058 "the set of rational numbers, i. e. the predicate will also be\n"
3059 "fulfilled if @var{x} is an integer number.")
3060 #define FUNC_NAME s_scm_rational_p
3062 if (SCM_I_INUMP (x
))
3064 else if (SCM_IMP (x
))
3066 else if (SCM_BIGP (x
))
3068 else if (SCM_FRACTIONP (x
))
3070 else if (SCM_REALP (x
))
3071 /* due to their limited precision, all floating point numbers are
3072 rational as well. */
3079 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3081 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3083 #define FUNC_NAME s_scm_integer_p
3086 if (SCM_I_INUMP (x
))
3092 if (!SCM_INEXACTP (x
))
3094 if (SCM_COMPLEXP (x
))
3096 r
= SCM_REAL_VALUE (x
);
3104 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3106 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3108 #define FUNC_NAME s_scm_inexact_p
3110 if (SCM_INEXACTP (x
))
3112 if (SCM_NUMBERP (x
))
3114 SCM_WRONG_TYPE_ARG (1, x
);
3119 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3120 /* "Return @code{#t} if all parameters are numerically equal." */
3122 scm_num_eq_p (SCM x
, SCM y
)
3125 if (SCM_I_INUMP (x
))
3127 long xx
= SCM_I_INUM (x
);
3128 if (SCM_I_INUMP (y
))
3130 long yy
= SCM_I_INUM (y
);
3131 return scm_from_bool (xx
== yy
);
3133 else if (SCM_BIGP (y
))
3135 else if (SCM_REALP (y
))
3136 return scm_from_bool ((double) xx
== SCM_REAL_VALUE (y
));
3137 else if (SCM_COMPLEXP (y
))
3138 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3139 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3140 else if (SCM_FRACTIONP (y
))
3143 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3145 else if (SCM_BIGP (x
))
3147 if (SCM_I_INUMP (y
))
3149 else if (SCM_BIGP (y
))
3151 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3152 scm_remember_upto_here_2 (x
, y
);
3153 return scm_from_bool (0 == cmp
);
3155 else if (SCM_REALP (y
))
3158 if (xisnan (SCM_REAL_VALUE (y
)))
3160 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3161 scm_remember_upto_here_1 (x
);
3162 return scm_from_bool (0 == cmp
);
3164 else if (SCM_COMPLEXP (y
))
3167 if (0.0 != SCM_COMPLEX_IMAG (y
))
3169 if (xisnan (SCM_COMPLEX_REAL (y
)))
3171 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3172 scm_remember_upto_here_1 (x
);
3173 return scm_from_bool (0 == cmp
);
3175 else if (SCM_FRACTIONP (y
))
3178 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3180 else if (SCM_REALP (x
))
3182 if (SCM_I_INUMP (y
))
3183 return scm_from_bool (SCM_REAL_VALUE (x
) == (double) SCM_I_INUM (y
));
3184 else if (SCM_BIGP (y
))
3187 if (xisnan (SCM_REAL_VALUE (x
)))
3189 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3190 scm_remember_upto_here_1 (y
);
3191 return scm_from_bool (0 == cmp
);
3193 else if (SCM_REALP (y
))
3194 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3195 else if (SCM_COMPLEXP (y
))
3196 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3197 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3198 else if (SCM_FRACTIONP (y
))
3200 double xx
= SCM_REAL_VALUE (x
);
3204 return scm_from_bool (xx
< 0.0);
3205 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3209 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3211 else if (SCM_COMPLEXP (x
))
3213 if (SCM_I_INUMP (y
))
3214 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3215 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3216 else if (SCM_BIGP (y
))
3219 if (0.0 != SCM_COMPLEX_IMAG (x
))
3221 if (xisnan (SCM_COMPLEX_REAL (x
)))
3223 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3224 scm_remember_upto_here_1 (y
);
3225 return scm_from_bool (0 == cmp
);
3227 else if (SCM_REALP (y
))
3228 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3229 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3230 else if (SCM_COMPLEXP (y
))
3231 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3232 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3233 else if (SCM_FRACTIONP (y
))
3236 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3238 xx
= SCM_COMPLEX_REAL (x
);
3242 return scm_from_bool (xx
< 0.0);
3243 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3247 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3249 else if (SCM_FRACTIONP (x
))
3251 if (SCM_I_INUMP (y
))
3253 else if (SCM_BIGP (y
))
3255 else if (SCM_REALP (y
))
3257 double yy
= SCM_REAL_VALUE (y
);
3261 return scm_from_bool (0.0 < yy
);
3262 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3265 else if (SCM_COMPLEXP (y
))
3268 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3270 yy
= SCM_COMPLEX_REAL (y
);
3274 return scm_from_bool (0.0 < yy
);
3275 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3278 else if (SCM_FRACTIONP (y
))
3279 return scm_i_fraction_equalp (x
, y
);
3281 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3284 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3288 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3289 done are good for inums, but for bignums an answer can almost always be
3290 had by just examining a few high bits of the operands, as done by GMP in
3291 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3292 of the float exponent to take into account. */
3294 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3295 /* "Return @code{#t} if the list of parameters is monotonically\n"
3299 scm_less_p (SCM x
, SCM y
)
3302 if (SCM_I_INUMP (x
))
3304 long xx
= SCM_I_INUM (x
);
3305 if (SCM_I_INUMP (y
))
3307 long yy
= SCM_I_INUM (y
);
3308 return scm_from_bool (xx
< yy
);
3310 else if (SCM_BIGP (y
))
3312 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3313 scm_remember_upto_here_1 (y
);
3314 return scm_from_bool (sgn
> 0);
3316 else if (SCM_REALP (y
))
3317 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3318 else if (SCM_FRACTIONP (y
))
3320 /* "x < a/b" becomes "x*b < a" */
3322 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3323 y
= SCM_FRACTION_NUMERATOR (y
);
3327 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3329 else if (SCM_BIGP (x
))
3331 if (SCM_I_INUMP (y
))
3333 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3334 scm_remember_upto_here_1 (x
);
3335 return scm_from_bool (sgn
< 0);
3337 else if (SCM_BIGP (y
))
3339 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3340 scm_remember_upto_here_2 (x
, y
);
3341 return scm_from_bool (cmp
< 0);
3343 else if (SCM_REALP (y
))
3346 if (xisnan (SCM_REAL_VALUE (y
)))
3348 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3349 scm_remember_upto_here_1 (x
);
3350 return scm_from_bool (cmp
< 0);
3352 else if (SCM_FRACTIONP (y
))
3355 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3357 else if (SCM_REALP (x
))
3359 if (SCM_I_INUMP (y
))
3360 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3361 else if (SCM_BIGP (y
))
3364 if (xisnan (SCM_REAL_VALUE (x
)))
3366 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3367 scm_remember_upto_here_1 (y
);
3368 return scm_from_bool (cmp
> 0);
3370 else if (SCM_REALP (y
))
3371 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3372 else if (SCM_FRACTIONP (y
))
3374 double xx
= SCM_REAL_VALUE (x
);
3378 return scm_from_bool (xx
< 0.0);
3379 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3383 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3385 else if (SCM_FRACTIONP (x
))
3387 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3389 /* "a/b < y" becomes "a < y*b" */
3390 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3391 x
= SCM_FRACTION_NUMERATOR (x
);
3394 else if (SCM_REALP (y
))
3396 double yy
= SCM_REAL_VALUE (y
);
3400 return scm_from_bool (0.0 < yy
);
3401 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3404 else if (SCM_FRACTIONP (y
))
3406 /* "a/b < c/d" becomes "a*d < c*b" */
3407 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3408 SCM_FRACTION_DENOMINATOR (y
));
3409 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3410 SCM_FRACTION_DENOMINATOR (x
));
3416 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3419 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3423 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3424 /* "Return @code{#t} if the list of parameters is monotonically\n"
3427 #define FUNC_NAME s_scm_gr_p
3429 scm_gr_p (SCM x
, SCM y
)
3431 if (!SCM_NUMBERP (x
))
3432 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3433 else if (!SCM_NUMBERP (y
))
3434 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3436 return scm_less_p (y
, x
);
3441 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3442 /* "Return @code{#t} if the list of parameters is monotonically\n"
3445 #define FUNC_NAME s_scm_leq_p
3447 scm_leq_p (SCM x
, SCM y
)
3449 if (!SCM_NUMBERP (x
))
3450 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3451 else if (!SCM_NUMBERP (y
))
3452 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3453 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3456 return scm_not (scm_less_p (y
, x
));
3461 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3462 /* "Return @code{#t} if the list of parameters is monotonically\n"
3465 #define FUNC_NAME s_scm_geq_p
3467 scm_geq_p (SCM x
, SCM y
)
3469 if (!SCM_NUMBERP (x
))
3470 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3471 else if (!SCM_NUMBERP (y
))
3472 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3473 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3476 return scm_not (scm_less_p (x
, y
));
3481 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3482 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3488 if (SCM_I_INUMP (z
))
3489 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3490 else if (SCM_BIGP (z
))
3492 else if (SCM_REALP (z
))
3493 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3494 else if (SCM_COMPLEXP (z
))
3495 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3496 && SCM_COMPLEX_IMAG (z
) == 0.0);
3497 else if (SCM_FRACTIONP (z
))
3500 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3504 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3505 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3509 scm_positive_p (SCM x
)
3511 if (SCM_I_INUMP (x
))
3512 return scm_from_bool (SCM_I_INUM (x
) > 0);
3513 else if (SCM_BIGP (x
))
3515 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3516 scm_remember_upto_here_1 (x
);
3517 return scm_from_bool (sgn
> 0);
3519 else if (SCM_REALP (x
))
3520 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3521 else if (SCM_FRACTIONP (x
))
3522 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3524 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3528 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3529 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3533 scm_negative_p (SCM x
)
3535 if (SCM_I_INUMP (x
))
3536 return scm_from_bool (SCM_I_INUM (x
) < 0);
3537 else if (SCM_BIGP (x
))
3539 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3540 scm_remember_upto_here_1 (x
);
3541 return scm_from_bool (sgn
< 0);
3543 else if (SCM_REALP (x
))
3544 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3545 else if (SCM_FRACTIONP (x
))
3546 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3548 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3552 /* scm_min and scm_max return an inexact when either argument is inexact, as
3553 required by r5rs. On that basis, for exact/inexact combinations the
3554 exact is converted to inexact to compare and possibly return. This is
3555 unlike scm_less_p above which takes some trouble to preserve all bits in
3556 its test, such trouble is not required for min and max. */
3558 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3559 /* "Return the maximum of all parameter values."
3562 scm_max (SCM x
, SCM y
)
3567 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3568 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3571 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3574 if (SCM_I_INUMP (x
))
3576 long xx
= SCM_I_INUM (x
);
3577 if (SCM_I_INUMP (y
))
3579 long yy
= SCM_I_INUM (y
);
3580 return (xx
< yy
) ? y
: x
;
3582 else if (SCM_BIGP (y
))
3584 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3585 scm_remember_upto_here_1 (y
);
3586 return (sgn
< 0) ? x
: y
;
3588 else if (SCM_REALP (y
))
3591 /* if y==NaN then ">" is false and we return NaN */
3592 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3594 else if (SCM_FRACTIONP (y
))
3597 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3600 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3602 else if (SCM_BIGP (x
))
3604 if (SCM_I_INUMP (y
))
3606 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3607 scm_remember_upto_here_1 (x
);
3608 return (sgn
< 0) ? y
: x
;
3610 else if (SCM_BIGP (y
))
3612 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3613 scm_remember_upto_here_2 (x
, y
);
3614 return (cmp
> 0) ? x
: y
;
3616 else if (SCM_REALP (y
))
3618 /* if y==NaN then xx>yy is false, so we return the NaN y */
3621 xx
= scm_i_big2dbl (x
);
3622 yy
= SCM_REAL_VALUE (y
);
3623 return (xx
> yy
? scm_from_double (xx
) : y
);
3625 else if (SCM_FRACTIONP (y
))
3630 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3632 else if (SCM_REALP (x
))
3634 if (SCM_I_INUMP (y
))
3636 double z
= SCM_I_INUM (y
);
3637 /* if x==NaN then "<" is false and we return NaN */
3638 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3640 else if (SCM_BIGP (y
))
3645 else if (SCM_REALP (y
))
3647 /* if x==NaN then our explicit check means we return NaN
3648 if y==NaN then ">" is false and we return NaN
3649 calling isnan is unavoidable, since it's the only way to know
3650 which of x or y causes any compares to be false */
3651 double xx
= SCM_REAL_VALUE (x
);
3652 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3654 else if (SCM_FRACTIONP (y
))
3656 double yy
= scm_i_fraction2double (y
);
3657 double xx
= SCM_REAL_VALUE (x
);
3658 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3661 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3663 else if (SCM_FRACTIONP (x
))
3665 if (SCM_I_INUMP (y
))
3669 else if (SCM_BIGP (y
))
3673 else if (SCM_REALP (y
))
3675 double xx
= scm_i_fraction2double (x
);
3676 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3678 else if (SCM_FRACTIONP (y
))
3683 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3686 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3690 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3691 /* "Return the minium of all parameter values."
3694 scm_min (SCM x
, SCM y
)
3699 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3700 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3703 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3706 if (SCM_I_INUMP (x
))
3708 long xx
= SCM_I_INUM (x
);
3709 if (SCM_I_INUMP (y
))
3711 long yy
= SCM_I_INUM (y
);
3712 return (xx
< yy
) ? x
: y
;
3714 else if (SCM_BIGP (y
))
3716 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3717 scm_remember_upto_here_1 (y
);
3718 return (sgn
< 0) ? y
: x
;
3720 else if (SCM_REALP (y
))
3723 /* if y==NaN then "<" is false and we return NaN */
3724 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3726 else if (SCM_FRACTIONP (y
))
3729 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3732 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3734 else if (SCM_BIGP (x
))
3736 if (SCM_I_INUMP (y
))
3738 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3739 scm_remember_upto_here_1 (x
);
3740 return (sgn
< 0) ? x
: y
;
3742 else if (SCM_BIGP (y
))
3744 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3745 scm_remember_upto_here_2 (x
, y
);
3746 return (cmp
> 0) ? y
: x
;
3748 else if (SCM_REALP (y
))
3750 /* if y==NaN then xx<yy is false, so we return the NaN y */
3753 xx
= scm_i_big2dbl (x
);
3754 yy
= SCM_REAL_VALUE (y
);
3755 return (xx
< yy
? scm_from_double (xx
) : y
);
3757 else if (SCM_FRACTIONP (y
))
3762 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3764 else if (SCM_REALP (x
))
3766 if (SCM_I_INUMP (y
))
3768 double z
= SCM_I_INUM (y
);
3769 /* if x==NaN then "<" is false and we return NaN */
3770 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3772 else if (SCM_BIGP (y
))
3777 else if (SCM_REALP (y
))
3779 /* if x==NaN then our explicit check means we return NaN
3780 if y==NaN then "<" is false and we return NaN
3781 calling isnan is unavoidable, since it's the only way to know
3782 which of x or y causes any compares to be false */
3783 double xx
= SCM_REAL_VALUE (x
);
3784 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3786 else if (SCM_FRACTIONP (y
))
3788 double yy
= scm_i_fraction2double (y
);
3789 double xx
= SCM_REAL_VALUE (x
);
3790 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3793 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3795 else if (SCM_FRACTIONP (x
))
3797 if (SCM_I_INUMP (y
))
3801 else if (SCM_BIGP (y
))
3805 else if (SCM_REALP (y
))
3807 double xx
= scm_i_fraction2double (x
);
3808 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3810 else if (SCM_FRACTIONP (y
))
3815 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3818 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3822 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3823 /* "Return the sum of all parameter values. Return 0 if called without\n"
3827 scm_sum (SCM x
, SCM y
)
3831 if (SCM_NUMBERP (x
)) return x
;
3832 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3833 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3836 if (SCM_I_INUMP (x
))
3838 if (SCM_I_INUMP (y
))
3840 long xx
= SCM_I_INUM (x
);
3841 long yy
= SCM_I_INUM (y
);
3842 long int z
= xx
+ yy
;
3843 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3845 else if (SCM_BIGP (y
))
3850 else if (SCM_REALP (y
))
3852 long int xx
= SCM_I_INUM (x
);
3853 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3855 else if (SCM_COMPLEXP (y
))
3857 long int xx
= SCM_I_INUM (x
);
3858 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3859 SCM_COMPLEX_IMAG (y
));
3861 else if (SCM_FRACTIONP (y
))
3862 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3863 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3864 SCM_FRACTION_DENOMINATOR (y
));
3866 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3867 } else if (SCM_BIGP (x
))
3869 if (SCM_I_INUMP (y
))
3874 inum
= SCM_I_INUM (y
);
3877 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3880 SCM result
= scm_i_mkbig ();
3881 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3882 scm_remember_upto_here_1 (x
);
3883 /* we know the result will have to be a bignum */
3886 return scm_i_normbig (result
);
3890 SCM result
= scm_i_mkbig ();
3891 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3892 scm_remember_upto_here_1 (x
);
3893 /* we know the result will have to be a bignum */
3896 return scm_i_normbig (result
);
3899 else if (SCM_BIGP (y
))
3901 SCM result
= scm_i_mkbig ();
3902 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3903 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3904 mpz_add (SCM_I_BIG_MPZ (result
),
3907 scm_remember_upto_here_2 (x
, y
);
3908 /* we know the result will have to be a bignum */
3911 return scm_i_normbig (result
);
3913 else if (SCM_REALP (y
))
3915 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3916 scm_remember_upto_here_1 (x
);
3917 return scm_from_double (result
);
3919 else if (SCM_COMPLEXP (y
))
3921 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3922 + SCM_COMPLEX_REAL (y
));
3923 scm_remember_upto_here_1 (x
);
3924 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
3926 else if (SCM_FRACTIONP (y
))
3927 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3928 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3929 SCM_FRACTION_DENOMINATOR (y
));
3931 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3933 else if (SCM_REALP (x
))
3935 if (SCM_I_INUMP (y
))
3936 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
3937 else if (SCM_BIGP (y
))
3939 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3940 scm_remember_upto_here_1 (y
);
3941 return scm_from_double (result
);
3943 else if (SCM_REALP (y
))
3944 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3945 else if (SCM_COMPLEXP (y
))
3946 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3947 SCM_COMPLEX_IMAG (y
));
3948 else if (SCM_FRACTIONP (y
))
3949 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3951 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3953 else if (SCM_COMPLEXP (x
))
3955 if (SCM_I_INUMP (y
))
3956 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
3957 SCM_COMPLEX_IMAG (x
));
3958 else if (SCM_BIGP (y
))
3960 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3961 + SCM_COMPLEX_REAL (x
));
3962 scm_remember_upto_here_1 (y
);
3963 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
3965 else if (SCM_REALP (y
))
3966 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3967 SCM_COMPLEX_IMAG (x
));
3968 else if (SCM_COMPLEXP (y
))
3969 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3970 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3971 else if (SCM_FRACTIONP (y
))
3972 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3973 SCM_COMPLEX_IMAG (x
));
3975 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3977 else if (SCM_FRACTIONP (x
))
3979 if (SCM_I_INUMP (y
))
3980 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3981 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3982 SCM_FRACTION_DENOMINATOR (x
));
3983 else if (SCM_BIGP (y
))
3984 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3985 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3986 SCM_FRACTION_DENOMINATOR (x
));
3987 else if (SCM_REALP (y
))
3988 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
3989 else if (SCM_COMPLEXP (y
))
3990 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
3991 SCM_COMPLEX_IMAG (y
));
3992 else if (SCM_FRACTIONP (y
))
3993 /* a/b + c/d = (ad + bc) / bd */
3994 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
3995 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
3996 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
3998 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4001 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4005 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4006 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4007 * the sum of all but the first argument are subtracted from the first
4009 #define FUNC_NAME s_difference
4011 scm_difference (SCM x
, SCM y
)
4016 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4018 if (SCM_I_INUMP (x
))
4020 long xx
= -SCM_I_INUM (x
);
4021 if (SCM_FIXABLE (xx
))
4022 return SCM_I_MAKINUM (xx
);
4024 return scm_i_long2big (xx
);
4026 else if (SCM_BIGP (x
))
4027 /* FIXME: do we really need to normalize here? */
4028 return scm_i_normbig (scm_i_clonebig (x
, 0));
4029 else if (SCM_REALP (x
))
4030 return scm_from_double (-SCM_REAL_VALUE (x
));
4031 else if (SCM_COMPLEXP (x
))
4032 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4033 -SCM_COMPLEX_IMAG (x
));
4034 else if (SCM_FRACTIONP (x
))
4035 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4036 SCM_FRACTION_DENOMINATOR (x
));
4038 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4041 if (SCM_I_INUMP (x
))
4043 if (SCM_I_INUMP (y
))
4045 long int xx
= SCM_I_INUM (x
);
4046 long int yy
= SCM_I_INUM (y
);
4047 long int z
= xx
- yy
;
4048 if (SCM_FIXABLE (z
))
4049 return SCM_I_MAKINUM (z
);
4051 return scm_i_long2big (z
);
4053 else if (SCM_BIGP (y
))
4055 /* inum-x - big-y */
4056 long xx
= SCM_I_INUM (x
);
4059 return scm_i_clonebig (y
, 0);
4062 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4063 SCM result
= scm_i_mkbig ();
4066 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4069 /* x - y == -(y + -x) */
4070 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4071 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4073 scm_remember_upto_here_1 (y
);
4075 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4076 /* we know the result will have to be a bignum */
4079 return scm_i_normbig (result
);
4082 else if (SCM_REALP (y
))
4084 long int xx
= SCM_I_INUM (x
);
4085 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4087 else if (SCM_COMPLEXP (y
))
4089 long int xx
= SCM_I_INUM (x
);
4090 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4091 - SCM_COMPLEX_IMAG (y
));
4093 else if (SCM_FRACTIONP (y
))
4094 /* a - b/c = (ac - b) / c */
4095 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4096 SCM_FRACTION_NUMERATOR (y
)),
4097 SCM_FRACTION_DENOMINATOR (y
));
4099 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4101 else if (SCM_BIGP (x
))
4103 if (SCM_I_INUMP (y
))
4105 /* big-x - inum-y */
4106 long yy
= SCM_I_INUM (y
);
4107 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4109 scm_remember_upto_here_1 (x
);
4111 return (SCM_FIXABLE (-yy
) ?
4112 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4115 SCM result
= scm_i_mkbig ();
4118 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4120 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4121 scm_remember_upto_here_1 (x
);
4123 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4124 /* we know the result will have to be a bignum */
4127 return scm_i_normbig (result
);
4130 else if (SCM_BIGP (y
))
4132 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4133 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4134 SCM result
= scm_i_mkbig ();
4135 mpz_sub (SCM_I_BIG_MPZ (result
),
4138 scm_remember_upto_here_2 (x
, y
);
4139 /* we know the result will have to be a bignum */
4140 if ((sgn_x
== 1) && (sgn_y
== -1))
4142 if ((sgn_x
== -1) && (sgn_y
== 1))
4144 return scm_i_normbig (result
);
4146 else if (SCM_REALP (y
))
4148 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4149 scm_remember_upto_here_1 (x
);
4150 return scm_from_double (result
);
4152 else if (SCM_COMPLEXP (y
))
4154 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4155 - SCM_COMPLEX_REAL (y
));
4156 scm_remember_upto_here_1 (x
);
4157 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4159 else if (SCM_FRACTIONP (y
))
4160 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4161 SCM_FRACTION_NUMERATOR (y
)),
4162 SCM_FRACTION_DENOMINATOR (y
));
4163 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4165 else if (SCM_REALP (x
))
4167 if (SCM_I_INUMP (y
))
4168 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4169 else if (SCM_BIGP (y
))
4171 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4172 scm_remember_upto_here_1 (x
);
4173 return scm_from_double (result
);
4175 else if (SCM_REALP (y
))
4176 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4177 else if (SCM_COMPLEXP (y
))
4178 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4179 -SCM_COMPLEX_IMAG (y
));
4180 else if (SCM_FRACTIONP (y
))
4181 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4183 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4185 else if (SCM_COMPLEXP (x
))
4187 if (SCM_I_INUMP (y
))
4188 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4189 SCM_COMPLEX_IMAG (x
));
4190 else if (SCM_BIGP (y
))
4192 double real_part
= (SCM_COMPLEX_REAL (x
)
4193 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4194 scm_remember_upto_here_1 (x
);
4195 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4197 else if (SCM_REALP (y
))
4198 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4199 SCM_COMPLEX_IMAG (x
));
4200 else if (SCM_COMPLEXP (y
))
4201 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4202 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4203 else if (SCM_FRACTIONP (y
))
4204 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4205 SCM_COMPLEX_IMAG (x
));
4207 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4209 else if (SCM_FRACTIONP (x
))
4211 if (SCM_I_INUMP (y
))
4212 /* a/b - c = (a - cb) / b */
4213 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4214 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4215 SCM_FRACTION_DENOMINATOR (x
));
4216 else if (SCM_BIGP (y
))
4217 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4218 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4219 SCM_FRACTION_DENOMINATOR (x
));
4220 else if (SCM_REALP (y
))
4221 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4222 else if (SCM_COMPLEXP (y
))
4223 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4224 -SCM_COMPLEX_IMAG (y
));
4225 else if (SCM_FRACTIONP (y
))
4226 /* a/b - c/d = (ad - bc) / bd */
4227 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4228 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4229 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4231 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4234 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4239 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4240 /* "Return the product of all arguments. If called without arguments,\n"
4244 scm_product (SCM x
, SCM y
)
4249 return SCM_I_MAKINUM (1L);
4250 else if (SCM_NUMBERP (x
))
4253 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4256 if (SCM_I_INUMP (x
))
4261 xx
= SCM_I_INUM (x
);
4265 case 0: return x
; break;
4266 case 1: return y
; break;
4269 if (SCM_I_INUMP (y
))
4271 long yy
= SCM_I_INUM (y
);
4273 SCM k
= SCM_I_MAKINUM (kk
);
4274 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4278 SCM result
= scm_i_long2big (xx
);
4279 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4280 return scm_i_normbig (result
);
4283 else if (SCM_BIGP (y
))
4285 SCM result
= scm_i_mkbig ();
4286 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4287 scm_remember_upto_here_1 (y
);
4290 else if (SCM_REALP (y
))
4291 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4292 else if (SCM_COMPLEXP (y
))
4293 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4294 xx
* SCM_COMPLEX_IMAG (y
));
4295 else if (SCM_FRACTIONP (y
))
4296 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4297 SCM_FRACTION_DENOMINATOR (y
));
4299 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4301 else if (SCM_BIGP (x
))
4303 if (SCM_I_INUMP (y
))
4308 else if (SCM_BIGP (y
))
4310 SCM result
= scm_i_mkbig ();
4311 mpz_mul (SCM_I_BIG_MPZ (result
),
4314 scm_remember_upto_here_2 (x
, y
);
4317 else if (SCM_REALP (y
))
4319 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4320 scm_remember_upto_here_1 (x
);
4321 return scm_from_double (result
);
4323 else if (SCM_COMPLEXP (y
))
4325 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4326 scm_remember_upto_here_1 (x
);
4327 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4328 z
* SCM_COMPLEX_IMAG (y
));
4330 else if (SCM_FRACTIONP (y
))
4331 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4332 SCM_FRACTION_DENOMINATOR (y
));
4334 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4336 else if (SCM_REALP (x
))
4338 if (SCM_I_INUMP (y
))
4339 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4340 else if (SCM_BIGP (y
))
4342 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4343 scm_remember_upto_here_1 (y
);
4344 return scm_from_double (result
);
4346 else if (SCM_REALP (y
))
4347 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4348 else if (SCM_COMPLEXP (y
))
4349 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4350 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4351 else if (SCM_FRACTIONP (y
))
4352 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4354 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4356 else if (SCM_COMPLEXP (x
))
4358 if (SCM_I_INUMP (y
))
4359 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4360 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4361 else if (SCM_BIGP (y
))
4363 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4364 scm_remember_upto_here_1 (y
);
4365 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4366 z
* SCM_COMPLEX_IMAG (x
));
4368 else if (SCM_REALP (y
))
4369 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4370 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4371 else if (SCM_COMPLEXP (y
))
4373 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4374 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4375 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4376 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4378 else if (SCM_FRACTIONP (y
))
4380 double yy
= scm_i_fraction2double (y
);
4381 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4382 yy
* SCM_COMPLEX_IMAG (x
));
4385 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4387 else if (SCM_FRACTIONP (x
))
4389 if (SCM_I_INUMP (y
))
4390 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4391 SCM_FRACTION_DENOMINATOR (x
));
4392 else if (SCM_BIGP (y
))
4393 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4394 SCM_FRACTION_DENOMINATOR (x
));
4395 else if (SCM_REALP (y
))
4396 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4397 else if (SCM_COMPLEXP (y
))
4399 double xx
= scm_i_fraction2double (x
);
4400 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4401 xx
* SCM_COMPLEX_IMAG (y
));
4403 else if (SCM_FRACTIONP (y
))
4404 /* a/b * c/d = ac / bd */
4405 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4406 SCM_FRACTION_NUMERATOR (y
)),
4407 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4408 SCM_FRACTION_DENOMINATOR (y
)));
4410 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4413 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4416 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4417 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4418 #define ALLOW_DIVIDE_BY_ZERO
4419 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4422 /* The code below for complex division is adapted from the GNU
4423 libstdc++, which adapted it from f2c's libF77, and is subject to
4426 /****************************************************************
4427 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4429 Permission to use, copy, modify, and distribute this software
4430 and its documentation for any purpose and without fee is hereby
4431 granted, provided that the above copyright notice appear in all
4432 copies and that both that the copyright notice and this
4433 permission notice and warranty disclaimer appear in supporting
4434 documentation, and that the names of AT&T Bell Laboratories or
4435 Bellcore or any of their entities not be used in advertising or
4436 publicity pertaining to distribution of the software without
4437 specific, written prior permission.
4439 AT&T and Bellcore disclaim all warranties with regard to this
4440 software, including all implied warranties of merchantability
4441 and fitness. In no event shall AT&T or Bellcore be liable for
4442 any special, indirect or consequential damages or any damages
4443 whatsoever resulting from loss of use, data or profits, whether
4444 in an action of contract, negligence or other tortious action,
4445 arising out of or in connection with the use or performance of
4447 ****************************************************************/
4449 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4450 /* Divide the first argument by the product of the remaining
4451 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4453 #define FUNC_NAME s_divide
4455 scm_i_divide (SCM x
, SCM y
, int inexact
)
4462 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4463 else if (SCM_I_INUMP (x
))
4465 long xx
= SCM_I_INUM (x
);
4466 if (xx
== 1 || xx
== -1)
4468 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4470 scm_num_overflow (s_divide
);
4475 return scm_from_double (1.0 / (double) xx
);
4476 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4479 else if (SCM_BIGP (x
))
4482 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4483 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4485 else if (SCM_REALP (x
))
4487 double xx
= SCM_REAL_VALUE (x
);
4488 #ifndef ALLOW_DIVIDE_BY_ZERO
4490 scm_num_overflow (s_divide
);
4493 return scm_from_double (1.0 / xx
);
4495 else if (SCM_COMPLEXP (x
))
4497 double r
= SCM_COMPLEX_REAL (x
);
4498 double i
= SCM_COMPLEX_IMAG (x
);
4502 double d
= i
* (1.0 + t
* t
);
4503 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4508 double d
= r
* (1.0 + t
* t
);
4509 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4512 else if (SCM_FRACTIONP (x
))
4513 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4514 SCM_FRACTION_NUMERATOR (x
));
4516 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4519 if (SCM_I_INUMP (x
))
4521 long xx
= SCM_I_INUM (x
);
4522 if (SCM_I_INUMP (y
))
4524 long yy
= SCM_I_INUM (y
);
4527 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4528 scm_num_overflow (s_divide
);
4530 return scm_from_double ((double) xx
/ (double) yy
);
4533 else if (xx
% yy
!= 0)
4536 return scm_from_double ((double) xx
/ (double) yy
);
4537 else return scm_i_make_ratio (x
, y
);
4542 if (SCM_FIXABLE (z
))
4543 return SCM_I_MAKINUM (z
);
4545 return scm_i_long2big (z
);
4548 else if (SCM_BIGP (y
))
4551 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4552 else return scm_i_make_ratio (x
, y
);
4554 else if (SCM_REALP (y
))
4556 double yy
= SCM_REAL_VALUE (y
);
4557 #ifndef ALLOW_DIVIDE_BY_ZERO
4559 scm_num_overflow (s_divide
);
4562 return scm_from_double ((double) xx
/ yy
);
4564 else if (SCM_COMPLEXP (y
))
4567 complex_div
: /* y _must_ be a complex number */
4569 double r
= SCM_COMPLEX_REAL (y
);
4570 double i
= SCM_COMPLEX_IMAG (y
);
4574 double d
= i
* (1.0 + t
* t
);
4575 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4580 double d
= r
* (1.0 + t
* t
);
4581 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4585 else if (SCM_FRACTIONP (y
))
4586 /* a / b/c = ac / b */
4587 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4588 SCM_FRACTION_NUMERATOR (y
));
4590 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4592 else if (SCM_BIGP (x
))
4594 if (SCM_I_INUMP (y
))
4596 long int yy
= SCM_I_INUM (y
);
4599 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4600 scm_num_overflow (s_divide
);
4602 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4603 scm_remember_upto_here_1 (x
);
4604 return (sgn
== 0) ? scm_nan () : scm_inf ();
4611 /* FIXME: HMM, what are the relative performance issues here?
4612 We need to test. Is it faster on average to test
4613 divisible_p, then perform whichever operation, or is it
4614 faster to perform the integer div opportunistically and
4615 switch to real if there's a remainder? For now we take the
4616 middle ground: test, then if divisible, use the faster div
4619 long abs_yy
= yy
< 0 ? -yy
: yy
;
4620 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4624 SCM result
= scm_i_mkbig ();
4625 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4626 scm_remember_upto_here_1 (x
);
4628 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4629 return scm_i_normbig (result
);
4634 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4635 else return scm_i_make_ratio (x
, y
);
4639 else if (SCM_BIGP (y
))
4641 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4644 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4645 scm_num_overflow (s_divide
);
4647 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4648 scm_remember_upto_here_1 (x
);
4649 return (sgn
== 0) ? scm_nan () : scm_inf ();
4655 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4659 SCM result
= scm_i_mkbig ();
4660 mpz_divexact (SCM_I_BIG_MPZ (result
),
4663 scm_remember_upto_here_2 (x
, y
);
4664 return scm_i_normbig (result
);
4670 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4671 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4672 scm_remember_upto_here_2 (x
, y
);
4673 return scm_from_double (dbx
/ dby
);
4675 else return scm_i_make_ratio (x
, y
);
4679 else if (SCM_REALP (y
))
4681 double yy
= SCM_REAL_VALUE (y
);
4682 #ifndef ALLOW_DIVIDE_BY_ZERO
4684 scm_num_overflow (s_divide
);
4687 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4689 else if (SCM_COMPLEXP (y
))
4691 a
= scm_i_big2dbl (x
);
4694 else if (SCM_FRACTIONP (y
))
4695 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4696 SCM_FRACTION_NUMERATOR (y
));
4698 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4700 else if (SCM_REALP (x
))
4702 double rx
= SCM_REAL_VALUE (x
);
4703 if (SCM_I_INUMP (y
))
4705 long int yy
= SCM_I_INUM (y
);
4706 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4708 scm_num_overflow (s_divide
);
4711 return scm_from_double (rx
/ (double) yy
);
4713 else if (SCM_BIGP (y
))
4715 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4716 scm_remember_upto_here_1 (y
);
4717 return scm_from_double (rx
/ dby
);
4719 else if (SCM_REALP (y
))
4721 double yy
= SCM_REAL_VALUE (y
);
4722 #ifndef ALLOW_DIVIDE_BY_ZERO
4724 scm_num_overflow (s_divide
);
4727 return scm_from_double (rx
/ yy
);
4729 else if (SCM_COMPLEXP (y
))
4734 else if (SCM_FRACTIONP (y
))
4735 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4737 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4739 else if (SCM_COMPLEXP (x
))
4741 double rx
= SCM_COMPLEX_REAL (x
);
4742 double ix
= SCM_COMPLEX_IMAG (x
);
4743 if (SCM_I_INUMP (y
))
4745 long int yy
= SCM_I_INUM (y
);
4746 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4748 scm_num_overflow (s_divide
);
4753 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4756 else if (SCM_BIGP (y
))
4758 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4759 scm_remember_upto_here_1 (y
);
4760 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4762 else if (SCM_REALP (y
))
4764 double yy
= SCM_REAL_VALUE (y
);
4765 #ifndef ALLOW_DIVIDE_BY_ZERO
4767 scm_num_overflow (s_divide
);
4770 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4772 else if (SCM_COMPLEXP (y
))
4774 double ry
= SCM_COMPLEX_REAL (y
);
4775 double iy
= SCM_COMPLEX_IMAG (y
);
4779 double d
= iy
* (1.0 + t
* t
);
4780 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4785 double d
= ry
* (1.0 + t
* t
);
4786 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4789 else if (SCM_FRACTIONP (y
))
4791 double yy
= scm_i_fraction2double (y
);
4792 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4795 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4797 else if (SCM_FRACTIONP (x
))
4799 if (SCM_I_INUMP (y
))
4801 long int yy
= SCM_I_INUM (y
);
4802 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4804 scm_num_overflow (s_divide
);
4807 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4808 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4810 else if (SCM_BIGP (y
))
4812 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4813 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4815 else if (SCM_REALP (y
))
4817 double yy
= SCM_REAL_VALUE (y
);
4818 #ifndef ALLOW_DIVIDE_BY_ZERO
4820 scm_num_overflow (s_divide
);
4823 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4825 else if (SCM_COMPLEXP (y
))
4827 a
= scm_i_fraction2double (x
);
4830 else if (SCM_FRACTIONP (y
))
4831 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4832 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4834 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4837 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4841 scm_divide (SCM x
, SCM y
)
4843 return scm_i_divide (x
, y
, 0);
4846 static SCM
scm_divide2real (SCM x
, SCM y
)
4848 return scm_i_divide (x
, y
, 1);
4854 scm_asinh (double x
)
4859 #define asinh scm_asinh
4860 return log (x
+ sqrt (x
* x
+ 1));
4863 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4864 /* "Return the inverse hyperbolic sine of @var{x}."
4869 scm_acosh (double x
)
4874 #define acosh scm_acosh
4875 return log (x
+ sqrt (x
* x
- 1));
4878 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4879 /* "Return the inverse hyperbolic cosine of @var{x}."
4884 scm_atanh (double x
)
4889 #define atanh scm_atanh
4890 return 0.5 * log ((1 + x
) / (1 - x
));
4893 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4894 /* "Return the inverse hyperbolic tangent of @var{x}."
4899 scm_c_truncate (double x
)
4910 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
4911 half-way case (ie. when x is an integer plus 0.5) going upwards.
4912 Then half-way cases are identified and adjusted down if the
4913 round-upwards didn't give the desired even integer.
4915 "plus_half == result" identifies a half-way case. If plus_half, which is
4916 x + 0.5, is an integer then x must be an integer plus 0.5.
4918 An odd "result" value is identified with result/2 != floor(result/2).
4919 This is done with plus_half, since that value is ready for use sooner in
4920 a pipelined cpu, and we're already requiring plus_half == result.
4922 Note however that we need to be careful when x is big and already an
4923 integer. In that case "x+0.5" may round to an adjacent integer, causing
4924 us to return such a value, incorrectly. For instance if the hardware is
4925 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
4926 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
4927 returned. Or if the hardware is in round-upwards mode, then other bigger
4928 values like say x == 2^128 will see x+0.5 rounding up to the next higher
4929 representable value, 2^128+2^76 (or whatever), again incorrect.
4931 These bad roundings of x+0.5 are avoided by testing at the start whether
4932 x is already an integer. If it is then clearly that's the desired result
4933 already. And if it's not then the exponent must be small enough to allow
4934 an 0.5 to be represented, and hence added without a bad rounding. */
4937 scm_c_round (double x
)
4939 double plus_half
, result
;
4944 plus_half
= x
+ 0.5;
4945 result
= floor (plus_half
);
4946 /* Adjust so that the rounding is towards even. */
4947 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4952 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4954 "Round the number @var{x} towards zero.")
4955 #define FUNC_NAME s_scm_truncate_number
4957 if (scm_is_false (scm_negative_p (x
)))
4958 return scm_floor (x
);
4960 return scm_ceiling (x
);
4964 static SCM exactly_one_half
;
4966 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4968 "Round the number @var{x} towards the nearest integer. "
4969 "When it is exactly halfway between two integers, "
4970 "round towards the even one.")
4971 #define FUNC_NAME s_scm_round_number
4973 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4975 else if (SCM_REALP (x
))
4976 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
4979 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
4980 single quotient+remainder division then examining to see which way
4981 the rounding should go. */
4982 SCM plus_half
= scm_sum (x
, exactly_one_half
);
4983 SCM result
= scm_floor (plus_half
);
4984 /* Adjust so that the rounding is towards even. */
4985 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
4986 && scm_is_true (scm_odd_p (result
)))
4987 return scm_difference (result
, SCM_I_MAKINUM (1));
4994 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
4996 "Round the number @var{x} towards minus infinity.")
4997 #define FUNC_NAME s_scm_floor
4999 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5001 else if (SCM_REALP (x
))
5002 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5003 else if (SCM_FRACTIONP (x
))
5005 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5006 SCM_FRACTION_DENOMINATOR (x
));
5007 if (scm_is_false (scm_negative_p (x
)))
5009 /* For positive x, rounding towards zero is correct. */
5014 /* For negative x, we need to return q-1 unless x is an
5015 integer. But fractions are never integer, per our
5017 return scm_difference (q
, SCM_I_MAKINUM (1));
5021 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5025 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5027 "Round the number @var{x} towards infinity.")
5028 #define FUNC_NAME s_scm_ceiling
5030 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5032 else if (SCM_REALP (x
))
5033 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5034 else if (SCM_FRACTIONP (x
))
5036 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5037 SCM_FRACTION_DENOMINATOR (x
));
5038 if (scm_is_false (scm_positive_p (x
)))
5040 /* For negative x, rounding towards zero is correct. */
5045 /* For positive x, we need to return q+1 unless x is an
5046 integer. But fractions are never integer, per our
5048 return scm_sum (q
, SCM_I_MAKINUM (1));
5052 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5056 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5057 /* "Return the square root of the real number @var{x}."
5059 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5060 /* "Return the absolute value of the real number @var{x}."
5062 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5063 /* "Return the @var{x}th power of e."
5065 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5066 /* "Return the natural logarithm of the real number @var{x}."
5068 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5069 /* "Return the sine of the real number @var{x}."
5071 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5072 /* "Return the cosine of the real number @var{x}."
5074 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5075 /* "Return the tangent of the real number @var{x}."
5077 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5078 /* "Return the arc sine of the real number @var{x}."
5080 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5081 /* "Return the arc cosine of the real number @var{x}."
5083 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5084 /* "Return the arc tangent of the real number @var{x}."
5086 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5087 /* "Return the hyperbolic sine of the real number @var{x}."
5089 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5090 /* "Return the hyperbolic cosine of the real number @var{x}."
5092 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5093 /* "Return the hyperbolic tangent of the real number @var{x}."
5101 static void scm_two_doubles (SCM x
,
5103 const char *sstring
,
5107 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5109 if (SCM_I_INUMP (x
))
5110 xy
->x
= SCM_I_INUM (x
);
5111 else if (SCM_BIGP (x
))
5112 xy
->x
= scm_i_big2dbl (x
);
5113 else if (SCM_REALP (x
))
5114 xy
->x
= SCM_REAL_VALUE (x
);
5115 else if (SCM_FRACTIONP (x
))
5116 xy
->x
= scm_i_fraction2double (x
);
5118 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5120 if (SCM_I_INUMP (y
))
5121 xy
->y
= SCM_I_INUM (y
);
5122 else if (SCM_BIGP (y
))
5123 xy
->y
= scm_i_big2dbl (y
);
5124 else if (SCM_REALP (y
))
5125 xy
->y
= SCM_REAL_VALUE (y
);
5126 else if (SCM_FRACTIONP (y
))
5127 xy
->y
= scm_i_fraction2double (y
);
5129 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5133 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5135 "Return @var{x} raised to the power of @var{y}. This\n"
5136 "procedure does not accept complex arguments.")
5137 #define FUNC_NAME s_scm_sys_expt
5140 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5141 return scm_from_double (pow (xy
.x
, xy
.y
));
5146 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5148 "Return the arc tangent of the two arguments @var{x} and\n"
5149 "@var{y}. This is similar to calculating the arc tangent of\n"
5150 "@var{x} / @var{y}, except that the signs of both arguments\n"
5151 "are used to determine the quadrant of the result. This\n"
5152 "procedure does not accept complex arguments.")
5153 #define FUNC_NAME s_scm_sys_atan2
5156 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5157 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5162 scm_c_make_rectangular (double re
, double im
)
5165 return scm_from_double (re
);
5169 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5171 SCM_COMPLEX_REAL (z
) = re
;
5172 SCM_COMPLEX_IMAG (z
) = im
;
5177 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5178 (SCM real
, SCM imaginary
),
5179 "Return a complex number constructed of the given @var{real} and\n"
5180 "@var{imaginary} parts.")
5181 #define FUNC_NAME s_scm_make_rectangular
5184 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5185 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5190 scm_c_make_polar (double mag
, double ang
)
5194 sincos (ang
, &s
, &c
);
5199 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5202 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5204 "Return the complex number @var{x} * e^(i * @var{y}).")
5205 #define FUNC_NAME s_scm_make_polar
5208 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5209 return scm_c_make_polar (xy
.x
, xy
.y
);
5214 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5215 /* "Return the real part of the number @var{z}."
5218 scm_real_part (SCM z
)
5220 if (SCM_I_INUMP (z
))
5222 else if (SCM_BIGP (z
))
5224 else if (SCM_REALP (z
))
5226 else if (SCM_COMPLEXP (z
))
5227 return scm_from_double (SCM_COMPLEX_REAL (z
));
5228 else if (SCM_FRACTIONP (z
))
5231 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5235 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5236 /* "Return the imaginary part of the number @var{z}."
5239 scm_imag_part (SCM z
)
5241 if (SCM_I_INUMP (z
))
5243 else if (SCM_BIGP (z
))
5245 else if (SCM_REALP (z
))
5247 else if (SCM_COMPLEXP (z
))
5248 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5249 else if (SCM_FRACTIONP (z
))
5252 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5255 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5256 /* "Return the numerator of the number @var{z}."
5259 scm_numerator (SCM z
)
5261 if (SCM_I_INUMP (z
))
5263 else if (SCM_BIGP (z
))
5265 else if (SCM_FRACTIONP (z
))
5267 scm_i_fraction_reduce (z
);
5268 return SCM_FRACTION_NUMERATOR (z
);
5270 else if (SCM_REALP (z
))
5271 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5273 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5277 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5278 /* "Return the denominator of the number @var{z}."
5281 scm_denominator (SCM z
)
5283 if (SCM_I_INUMP (z
))
5284 return SCM_I_MAKINUM (1);
5285 else if (SCM_BIGP (z
))
5286 return SCM_I_MAKINUM (1);
5287 else if (SCM_FRACTIONP (z
))
5289 scm_i_fraction_reduce (z
);
5290 return SCM_FRACTION_DENOMINATOR (z
);
5292 else if (SCM_REALP (z
))
5293 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5295 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5298 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5299 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5300 * "@code{abs} for real arguments, but also allows complex numbers."
5303 scm_magnitude (SCM z
)
5305 if (SCM_I_INUMP (z
))
5307 long int zz
= SCM_I_INUM (z
);
5310 else if (SCM_POSFIXABLE (-zz
))
5311 return SCM_I_MAKINUM (-zz
);
5313 return scm_i_long2big (-zz
);
5315 else if (SCM_BIGP (z
))
5317 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5318 scm_remember_upto_here_1 (z
);
5320 return scm_i_clonebig (z
, 0);
5324 else if (SCM_REALP (z
))
5325 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5326 else if (SCM_COMPLEXP (z
))
5327 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5328 else if (SCM_FRACTIONP (z
))
5330 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5332 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5333 SCM_FRACTION_DENOMINATOR (z
));
5336 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5340 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5341 /* "Return the angle of the complex number @var{z}."
5346 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5347 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5348 But if atan2 follows the floating point rounding mode, then the value
5349 is not a constant. Maybe it'd be close enough though. */
5350 if (SCM_I_INUMP (z
))
5352 if (SCM_I_INUM (z
) >= 0)
5355 return scm_from_double (atan2 (0.0, -1.0));
5357 else if (SCM_BIGP (z
))
5359 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5360 scm_remember_upto_here_1 (z
);
5362 return scm_from_double (atan2 (0.0, -1.0));
5366 else if (SCM_REALP (z
))
5368 if (SCM_REAL_VALUE (z
) >= 0)
5371 return scm_from_double (atan2 (0.0, -1.0));
5373 else if (SCM_COMPLEXP (z
))
5374 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5375 else if (SCM_FRACTIONP (z
))
5377 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5379 else return scm_from_double (atan2 (0.0, -1.0));
5382 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5386 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5387 /* Convert the number @var{x} to its inexact representation.\n"
5390 scm_exact_to_inexact (SCM z
)
5392 if (SCM_I_INUMP (z
))
5393 return scm_from_double ((double) SCM_I_INUM (z
));
5394 else if (SCM_BIGP (z
))
5395 return scm_from_double (scm_i_big2dbl (z
));
5396 else if (SCM_FRACTIONP (z
))
5397 return scm_from_double (scm_i_fraction2double (z
));
5398 else if (SCM_INEXACTP (z
))
5401 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5405 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5407 "Return an exact number that is numerically closest to @var{z}.")
5408 #define FUNC_NAME s_scm_inexact_to_exact
5410 if (SCM_I_INUMP (z
))
5412 else if (SCM_BIGP (z
))
5414 else if (SCM_REALP (z
))
5416 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5417 SCM_OUT_OF_RANGE (1, z
);
5424 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5425 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5426 scm_i_mpz2num (mpq_denref (frac
)));
5428 /* When scm_i_make_ratio throws, we leak the memory allocated
5435 else if (SCM_FRACTIONP (z
))
5438 SCM_WRONG_TYPE_ARG (1, z
);
5442 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5444 "Return an exact number that is within @var{err} of @var{x}.")
5445 #define FUNC_NAME s_scm_rationalize
5447 if (SCM_I_INUMP (x
))
5449 else if (SCM_BIGP (x
))
5451 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5453 /* Use continued fractions to find closest ratio. All
5454 arithmetic is done with exact numbers.
5457 SCM ex
= scm_inexact_to_exact (x
);
5458 SCM int_part
= scm_floor (ex
);
5459 SCM tt
= SCM_I_MAKINUM (1);
5460 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5461 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5465 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5468 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5469 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5471 /* We stop after a million iterations just to be absolutely sure
5472 that we don't go into an infinite loop. The process normally
5473 converges after less than a dozen iterations.
5476 err
= scm_abs (err
);
5477 while (++i
< 1000000)
5479 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5480 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5481 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5483 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5484 err
))) /* abs(x-a/b) <= err */
5486 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5487 if (scm_is_false (scm_exact_p (x
))
5488 || scm_is_false (scm_exact_p (err
)))
5489 return scm_exact_to_inexact (res
);
5493 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5495 tt
= scm_floor (rx
); /* tt = floor (rx) */
5501 scm_num_overflow (s_scm_rationalize
);
5504 SCM_WRONG_TYPE_ARG (1, x
);
5508 /* conversion functions */
5511 scm_is_integer (SCM val
)
5513 return scm_is_true (scm_integer_p (val
));
5517 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5519 if (SCM_I_INUMP (val
))
5521 scm_t_signed_bits n
= SCM_I_INUM (val
);
5522 return n
>= min
&& n
<= max
;
5524 else if (SCM_BIGP (val
))
5526 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5528 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5530 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5532 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5533 return n
>= min
&& n
<= max
;
5543 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5544 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5547 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5548 SCM_I_BIG_MPZ (val
));
5550 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5562 return n
>= min
&& n
<= max
;
5570 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5572 if (SCM_I_INUMP (val
))
5574 scm_t_signed_bits n
= SCM_I_INUM (val
);
5575 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5577 else if (SCM_BIGP (val
))
5579 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5581 else if (max
<= ULONG_MAX
)
5583 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5585 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5586 return n
>= min
&& n
<= max
;
5596 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5599 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5600 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5603 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5604 SCM_I_BIG_MPZ (val
));
5606 return n
>= min
&& n
<= max
;
5613 #define TYPE scm_t_intmax
5614 #define TYPE_MIN min
5615 #define TYPE_MAX max
5616 #define SIZEOF_TYPE 0
5617 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5618 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5619 #include "libguile/conv-integer.i.c"
5621 #define TYPE scm_t_uintmax
5622 #define TYPE_MIN min
5623 #define TYPE_MAX max
5624 #define SIZEOF_TYPE 0
5625 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5626 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5627 #include "libguile/conv-uinteger.i.c"
5629 #define TYPE scm_t_int8
5630 #define TYPE_MIN SCM_T_INT8_MIN
5631 #define TYPE_MAX SCM_T_INT8_MAX
5632 #define SIZEOF_TYPE 1
5633 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5634 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5635 #include "libguile/conv-integer.i.c"
5637 #define TYPE scm_t_uint8
5639 #define TYPE_MAX SCM_T_UINT8_MAX
5640 #define SIZEOF_TYPE 1
5641 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5642 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5643 #include "libguile/conv-uinteger.i.c"
5645 #define TYPE scm_t_int16
5646 #define TYPE_MIN SCM_T_INT16_MIN
5647 #define TYPE_MAX SCM_T_INT16_MAX
5648 #define SIZEOF_TYPE 2
5649 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5650 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5651 #include "libguile/conv-integer.i.c"
5653 #define TYPE scm_t_uint16
5655 #define TYPE_MAX SCM_T_UINT16_MAX
5656 #define SIZEOF_TYPE 2
5657 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5658 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5659 #include "libguile/conv-uinteger.i.c"
5661 #define TYPE scm_t_int32
5662 #define TYPE_MIN SCM_T_INT32_MIN
5663 #define TYPE_MAX SCM_T_INT32_MAX
5664 #define SIZEOF_TYPE 4
5665 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5666 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5667 #include "libguile/conv-integer.i.c"
5669 #define TYPE scm_t_uint32
5671 #define TYPE_MAX SCM_T_UINT32_MAX
5672 #define SIZEOF_TYPE 4
5673 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5674 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5675 #include "libguile/conv-uinteger.i.c"
5677 #if SCM_HAVE_T_INT64
5679 #define TYPE scm_t_int64
5680 #define TYPE_MIN SCM_T_INT64_MIN
5681 #define TYPE_MAX SCM_T_INT64_MAX
5682 #define SIZEOF_TYPE 8
5683 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5684 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5685 #include "libguile/conv-integer.i.c"
5687 #define TYPE scm_t_uint64
5689 #define TYPE_MAX SCM_T_UINT64_MAX
5690 #define SIZEOF_TYPE 8
5691 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5692 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5693 #include "libguile/conv-uinteger.i.c"
5698 scm_is_real (SCM val
)
5700 return scm_is_true (scm_real_p (val
));
5704 scm_is_rational (SCM val
)
5706 return scm_is_true (scm_rational_p (val
));
5710 scm_to_double (SCM val
)
5712 if (SCM_I_INUMP (val
))
5713 return SCM_I_INUM (val
);
5714 else if (SCM_BIGP (val
))
5715 return scm_i_big2dbl (val
);
5716 else if (SCM_FRACTIONP (val
))
5717 return scm_i_fraction2double (val
);
5718 else if (SCM_REALP (val
))
5719 return SCM_REAL_VALUE (val
);
5721 scm_wrong_type_arg (NULL
, 0, val
);
5725 scm_from_double (double val
)
5727 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5728 SCM_REAL_VALUE (z
) = val
;
5732 #if SCM_ENABLE_DISCOURAGED == 1
5735 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5739 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5743 scm_out_of_range (NULL
, num
);
5746 return scm_to_double (num
);
5750 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5754 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5758 scm_out_of_range (NULL
, num
);
5761 return scm_to_double (num
);
5767 scm_is_complex (SCM val
)
5769 return scm_is_true (scm_complex_p (val
));
5773 scm_c_real_part (SCM z
)
5775 if (SCM_COMPLEXP (z
))
5776 return SCM_COMPLEX_REAL (z
);
5779 /* Use the scm_real_part to get proper error checking and
5782 return scm_to_double (scm_real_part (z
));
5787 scm_c_imag_part (SCM z
)
5789 if (SCM_COMPLEXP (z
))
5790 return SCM_COMPLEX_IMAG (z
);
5793 /* Use the scm_imag_part to get proper error checking and
5794 dispatching. The result will almost always be 0.0, but not
5797 return scm_to_double (scm_imag_part (z
));
5802 scm_c_magnitude (SCM z
)
5804 return scm_to_double (scm_magnitude (z
));
5810 return scm_to_double (scm_angle (z
));
5814 scm_is_number (SCM z
)
5816 return scm_is_true (scm_number_p (z
));
5824 mpz_init_set_si (z_negative_one
, -1);
5826 /* It may be possible to tune the performance of some algorithms by using
5827 * the following constants to avoid the creation of bignums. Please, before
5828 * using these values, remember the two rules of program optimization:
5829 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5830 scm_c_define ("most-positive-fixnum",
5831 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5832 scm_c_define ("most-negative-fixnum",
5833 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5835 scm_add_feature ("complex");
5836 scm_add_feature ("inexact");
5837 scm_flo0
= scm_from_double (0.0);
5839 /* determine floating point precision */
5840 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
5842 init_dblprec(&scm_dblprec
[i
-2],i
);
5843 init_fx_radix(fx_per_radix
[i
-2],i
);
5846 /* hard code precision for base 10 if the preprocessor tells us to... */
5847 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5854 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
5855 SCM_I_MAKINUM (2)));
5856 #include "libguile/numbers.x"