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_I_STRING_CHARS (str
), SCM_I_STRING_LENGTH (str
), port
);
2342 scm_remember_upto_here_1 (str
);
2347 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2349 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2350 scm_remember_upto_here_1 (exp
);
2351 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2355 /*** END nums->strs ***/
2358 /*** STRINGS -> NUMBERS ***/
2360 /* The following functions implement the conversion from strings to numbers.
2361 * The implementation somehow follows the grammar for numbers as it is given
2362 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2363 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2364 * points should be noted about the implementation:
2365 * * Each function keeps a local index variable 'idx' that points at the
2366 * current position within the parsed string. The global index is only
2367 * updated if the function could parse the corresponding syntactic unit
2369 * * Similarly, the functions keep track of indicators of inexactness ('#',
2370 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2371 * global exactness information is only updated after each part has been
2372 * successfully parsed.
2373 * * Sequences of digits are parsed into temporary variables holding fixnums.
2374 * Only if these fixnums would overflow, the result variables are updated
2375 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2376 * the temporary variables holding the fixnums are cleared, and the process
2377 * starts over again. If for example fixnums were able to store five decimal
2378 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2379 * and the result was computed as 12345 * 100000 + 67890. In other words,
2380 * only every five digits two bignum operations were performed.
2383 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2385 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2387 /* In non ASCII-style encodings the following macro might not work. */
2388 #define XDIGIT2UINT(d) \
2389 (isdigit ((int) (unsigned char) d) \
2391 : tolower ((int) (unsigned char) d) - 'a' + 10)
2394 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2395 unsigned int radix
, enum t_exactness
*p_exactness
)
2397 unsigned int idx
= *p_idx
;
2398 unsigned int hash_seen
= 0;
2399 scm_t_bits shift
= 1;
2401 unsigned int digit_value
;
2409 if (!isxdigit ((int) (unsigned char) c
))
2411 digit_value
= XDIGIT2UINT (c
);
2412 if (digit_value
>= radix
)
2416 result
= SCM_I_MAKINUM (digit_value
);
2420 if (isxdigit ((int) (unsigned char) c
))
2424 digit_value
= XDIGIT2UINT (c
);
2425 if (digit_value
>= radix
)
2437 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2439 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2441 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2448 shift
= shift
* radix
;
2449 add
= add
* radix
+ digit_value
;
2454 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2456 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2460 *p_exactness
= INEXACT
;
2466 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2467 * covers the parts of the rules that start at a potential point. The value
2468 * of the digits up to the point have been parsed by the caller and are given
2469 * in variable result. The content of *p_exactness indicates, whether a hash
2470 * has already been seen in the digits before the point.
2473 /* In non ASCII-style encodings the following macro might not work. */
2474 #define DIGIT2UINT(d) ((d) - '0')
2477 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2478 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2480 unsigned int idx
= *p_idx
;
2481 enum t_exactness x
= *p_exactness
;
2486 if (mem
[idx
] == '.')
2488 scm_t_bits shift
= 1;
2490 unsigned int digit_value
;
2491 SCM big_shift
= SCM_I_MAKINUM (1);
2497 if (isdigit ((int) (unsigned char) c
))
2502 digit_value
= DIGIT2UINT (c
);
2513 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2515 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2516 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2518 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2526 add
= add
* 10 + digit_value
;
2532 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2533 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2534 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2537 result
= scm_divide (result
, big_shift
);
2539 /* We've seen a decimal point, thus the value is implicitly inexact. */
2551 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2578 if (!isdigit ((int) (unsigned char) c
))
2582 exponent
= DIGIT2UINT (c
);
2586 if (isdigit ((int) (unsigned char) c
))
2589 if (exponent
<= SCM_MAXEXP
)
2590 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2596 if (exponent
> SCM_MAXEXP
)
2598 size_t exp_len
= idx
- start
;
2599 SCM exp_string
= scm_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_I_STRING_CHARS (string
),
2971 SCM_I_STRING_LENGTH (string
),
2973 scm_remember_upto_here_1 (string
);
2979 /*** END strs->nums ***/
2983 scm_bigequal (SCM x
, SCM y
)
2985 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2986 scm_remember_upto_here_2 (x
, y
);
2987 return scm_from_bool (0 == result
);
2991 scm_real_equalp (SCM x
, SCM y
)
2993 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
2997 scm_complex_equalp (SCM x
, SCM y
)
2999 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3000 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3004 scm_i_fraction_equalp (SCM x
, SCM y
)
3006 scm_i_fraction_reduce (x
);
3007 scm_i_fraction_reduce (y
);
3008 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3009 SCM_FRACTION_NUMERATOR (y
)))
3010 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3011 SCM_FRACTION_DENOMINATOR (y
))))
3018 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3020 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3022 #define FUNC_NAME s_scm_number_p
3024 return scm_from_bool (SCM_NUMBERP (x
));
3028 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3030 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3031 "otherwise. Note that the sets of real, rational and integer\n"
3032 "values form subsets of the set of complex numbers, i. e. the\n"
3033 "predicate will also be fulfilled if @var{x} is a real,\n"
3034 "rational or integer number.")
3035 #define FUNC_NAME s_scm_complex_p
3037 /* all numbers are complex. */
3038 return scm_number_p (x
);
3042 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3044 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3045 "otherwise. Note that the set of integer values forms a subset of\n"
3046 "the set of real numbers, i. e. the predicate will also be\n"
3047 "fulfilled if @var{x} is an integer number.")
3048 #define FUNC_NAME s_scm_real_p
3050 /* we can't represent irrational numbers. */
3051 return scm_rational_p (x
);
3055 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3057 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3058 "otherwise. Note that the set of integer values forms a subset of\n"
3059 "the set of rational numbers, i. e. the predicate will also be\n"
3060 "fulfilled if @var{x} is an integer number.")
3061 #define FUNC_NAME s_scm_rational_p
3063 if (SCM_I_INUMP (x
))
3065 else if (SCM_IMP (x
))
3067 else if (SCM_BIGP (x
))
3069 else if (SCM_FRACTIONP (x
))
3071 else if (SCM_REALP (x
))
3072 /* due to their limited precision, all floating point numbers are
3073 rational as well. */
3080 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3082 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3084 #define FUNC_NAME s_scm_integer_p
3087 if (SCM_I_INUMP (x
))
3093 if (!SCM_INEXACTP (x
))
3095 if (SCM_COMPLEXP (x
))
3097 r
= SCM_REAL_VALUE (x
);
3105 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3107 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3109 #define FUNC_NAME s_scm_inexact_p
3111 if (SCM_INEXACTP (x
))
3113 if (SCM_NUMBERP (x
))
3115 SCM_WRONG_TYPE_ARG (1, x
);
3120 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3121 /* "Return @code{#t} if all parameters are numerically equal." */
3123 scm_num_eq_p (SCM x
, SCM y
)
3126 if (SCM_I_INUMP (x
))
3128 long xx
= SCM_I_INUM (x
);
3129 if (SCM_I_INUMP (y
))
3131 long yy
= SCM_I_INUM (y
);
3132 return scm_from_bool (xx
== yy
);
3134 else if (SCM_BIGP (y
))
3136 else if (SCM_REALP (y
))
3137 return scm_from_bool ((double) xx
== SCM_REAL_VALUE (y
));
3138 else if (SCM_COMPLEXP (y
))
3139 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3140 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3141 else if (SCM_FRACTIONP (y
))
3144 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3146 else if (SCM_BIGP (x
))
3148 if (SCM_I_INUMP (y
))
3150 else if (SCM_BIGP (y
))
3152 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3153 scm_remember_upto_here_2 (x
, y
);
3154 return scm_from_bool (0 == cmp
);
3156 else if (SCM_REALP (y
))
3159 if (xisnan (SCM_REAL_VALUE (y
)))
3161 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3162 scm_remember_upto_here_1 (x
);
3163 return scm_from_bool (0 == cmp
);
3165 else if (SCM_COMPLEXP (y
))
3168 if (0.0 != SCM_COMPLEX_IMAG (y
))
3170 if (xisnan (SCM_COMPLEX_REAL (y
)))
3172 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3173 scm_remember_upto_here_1 (x
);
3174 return scm_from_bool (0 == cmp
);
3176 else if (SCM_FRACTIONP (y
))
3179 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3181 else if (SCM_REALP (x
))
3183 if (SCM_I_INUMP (y
))
3184 return scm_from_bool (SCM_REAL_VALUE (x
) == (double) SCM_I_INUM (y
));
3185 else if (SCM_BIGP (y
))
3188 if (xisnan (SCM_REAL_VALUE (x
)))
3190 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3191 scm_remember_upto_here_1 (y
);
3192 return scm_from_bool (0 == cmp
);
3194 else if (SCM_REALP (y
))
3195 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3196 else if (SCM_COMPLEXP (y
))
3197 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3198 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3199 else if (SCM_FRACTIONP (y
))
3201 double xx
= SCM_REAL_VALUE (x
);
3205 return scm_from_bool (xx
< 0.0);
3206 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3210 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3212 else if (SCM_COMPLEXP (x
))
3214 if (SCM_I_INUMP (y
))
3215 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3216 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3217 else if (SCM_BIGP (y
))
3220 if (0.0 != SCM_COMPLEX_IMAG (x
))
3222 if (xisnan (SCM_COMPLEX_REAL (x
)))
3224 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3225 scm_remember_upto_here_1 (y
);
3226 return scm_from_bool (0 == cmp
);
3228 else if (SCM_REALP (y
))
3229 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3230 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3231 else if (SCM_COMPLEXP (y
))
3232 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3233 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3234 else if (SCM_FRACTIONP (y
))
3237 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3239 xx
= SCM_COMPLEX_REAL (x
);
3243 return scm_from_bool (xx
< 0.0);
3244 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3248 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3250 else if (SCM_FRACTIONP (x
))
3252 if (SCM_I_INUMP (y
))
3254 else if (SCM_BIGP (y
))
3256 else if (SCM_REALP (y
))
3258 double yy
= SCM_REAL_VALUE (y
);
3262 return scm_from_bool (0.0 < yy
);
3263 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3266 else if (SCM_COMPLEXP (y
))
3269 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3271 yy
= SCM_COMPLEX_REAL (y
);
3275 return scm_from_bool (0.0 < yy
);
3276 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3279 else if (SCM_FRACTIONP (y
))
3280 return scm_i_fraction_equalp (x
, y
);
3282 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3285 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3289 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3290 done are good for inums, but for bignums an answer can almost always be
3291 had by just examining a few high bits of the operands, as done by GMP in
3292 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3293 of the float exponent to take into account. */
3295 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3296 /* "Return @code{#t} if the list of parameters is monotonically\n"
3300 scm_less_p (SCM x
, SCM y
)
3303 if (SCM_I_INUMP (x
))
3305 long xx
= SCM_I_INUM (x
);
3306 if (SCM_I_INUMP (y
))
3308 long yy
= SCM_I_INUM (y
);
3309 return scm_from_bool (xx
< yy
);
3311 else if (SCM_BIGP (y
))
3313 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3314 scm_remember_upto_here_1 (y
);
3315 return scm_from_bool (sgn
> 0);
3317 else if (SCM_REALP (y
))
3318 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3319 else if (SCM_FRACTIONP (y
))
3321 /* "x < a/b" becomes "x*b < a" */
3323 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3324 y
= SCM_FRACTION_NUMERATOR (y
);
3328 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3330 else if (SCM_BIGP (x
))
3332 if (SCM_I_INUMP (y
))
3334 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3335 scm_remember_upto_here_1 (x
);
3336 return scm_from_bool (sgn
< 0);
3338 else if (SCM_BIGP (y
))
3340 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3341 scm_remember_upto_here_2 (x
, y
);
3342 return scm_from_bool (cmp
< 0);
3344 else if (SCM_REALP (y
))
3347 if (xisnan (SCM_REAL_VALUE (y
)))
3349 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3350 scm_remember_upto_here_1 (x
);
3351 return scm_from_bool (cmp
< 0);
3353 else if (SCM_FRACTIONP (y
))
3356 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3358 else if (SCM_REALP (x
))
3360 if (SCM_I_INUMP (y
))
3361 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3362 else if (SCM_BIGP (y
))
3365 if (xisnan (SCM_REAL_VALUE (x
)))
3367 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3368 scm_remember_upto_here_1 (y
);
3369 return scm_from_bool (cmp
> 0);
3371 else if (SCM_REALP (y
))
3372 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3373 else if (SCM_FRACTIONP (y
))
3375 double xx
= SCM_REAL_VALUE (x
);
3379 return scm_from_bool (xx
< 0.0);
3380 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3384 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3386 else if (SCM_FRACTIONP (x
))
3388 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3390 /* "a/b < y" becomes "a < y*b" */
3391 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3392 x
= SCM_FRACTION_NUMERATOR (x
);
3395 else if (SCM_REALP (y
))
3397 double yy
= SCM_REAL_VALUE (y
);
3401 return scm_from_bool (0.0 < yy
);
3402 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3405 else if (SCM_FRACTIONP (y
))
3407 /* "a/b < c/d" becomes "a*d < c*b" */
3408 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3409 SCM_FRACTION_DENOMINATOR (y
));
3410 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3411 SCM_FRACTION_DENOMINATOR (x
));
3417 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3420 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3424 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3425 /* "Return @code{#t} if the list of parameters is monotonically\n"
3428 #define FUNC_NAME s_scm_gr_p
3430 scm_gr_p (SCM x
, SCM y
)
3432 if (!SCM_NUMBERP (x
))
3433 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3434 else if (!SCM_NUMBERP (y
))
3435 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3437 return scm_less_p (y
, x
);
3442 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3443 /* "Return @code{#t} if the list of parameters is monotonically\n"
3446 #define FUNC_NAME s_scm_leq_p
3448 scm_leq_p (SCM x
, SCM y
)
3450 if (!SCM_NUMBERP (x
))
3451 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3452 else if (!SCM_NUMBERP (y
))
3453 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3454 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3457 return scm_not (scm_less_p (y
, x
));
3462 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3463 /* "Return @code{#t} if the list of parameters is monotonically\n"
3466 #define FUNC_NAME s_scm_geq_p
3468 scm_geq_p (SCM x
, SCM y
)
3470 if (!SCM_NUMBERP (x
))
3471 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3472 else if (!SCM_NUMBERP (y
))
3473 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3474 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3477 return scm_not (scm_less_p (x
, y
));
3482 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3483 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3489 if (SCM_I_INUMP (z
))
3490 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3491 else if (SCM_BIGP (z
))
3493 else if (SCM_REALP (z
))
3494 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3495 else if (SCM_COMPLEXP (z
))
3496 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3497 && SCM_COMPLEX_IMAG (z
) == 0.0);
3498 else if (SCM_FRACTIONP (z
))
3501 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3505 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3506 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3510 scm_positive_p (SCM x
)
3512 if (SCM_I_INUMP (x
))
3513 return scm_from_bool (SCM_I_INUM (x
) > 0);
3514 else if (SCM_BIGP (x
))
3516 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3517 scm_remember_upto_here_1 (x
);
3518 return scm_from_bool (sgn
> 0);
3520 else if (SCM_REALP (x
))
3521 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3522 else if (SCM_FRACTIONP (x
))
3523 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3525 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3529 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3530 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3534 scm_negative_p (SCM x
)
3536 if (SCM_I_INUMP (x
))
3537 return scm_from_bool (SCM_I_INUM (x
) < 0);
3538 else if (SCM_BIGP (x
))
3540 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3541 scm_remember_upto_here_1 (x
);
3542 return scm_from_bool (sgn
< 0);
3544 else if (SCM_REALP (x
))
3545 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3546 else if (SCM_FRACTIONP (x
))
3547 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3549 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3553 /* scm_min and scm_max return an inexact when either argument is inexact, as
3554 required by r5rs. On that basis, for exact/inexact combinations the
3555 exact is converted to inexact to compare and possibly return. This is
3556 unlike scm_less_p above which takes some trouble to preserve all bits in
3557 its test, such trouble is not required for min and max. */
3559 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3560 /* "Return the maximum of all parameter values."
3563 scm_max (SCM x
, SCM y
)
3568 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3569 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3572 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3575 if (SCM_I_INUMP (x
))
3577 long xx
= SCM_I_INUM (x
);
3578 if (SCM_I_INUMP (y
))
3580 long yy
= SCM_I_INUM (y
);
3581 return (xx
< yy
) ? y
: x
;
3583 else if (SCM_BIGP (y
))
3585 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3586 scm_remember_upto_here_1 (y
);
3587 return (sgn
< 0) ? x
: y
;
3589 else if (SCM_REALP (y
))
3592 /* if y==NaN then ">" is false and we return NaN */
3593 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3595 else if (SCM_FRACTIONP (y
))
3598 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3601 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3603 else if (SCM_BIGP (x
))
3605 if (SCM_I_INUMP (y
))
3607 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3608 scm_remember_upto_here_1 (x
);
3609 return (sgn
< 0) ? y
: x
;
3611 else if (SCM_BIGP (y
))
3613 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3614 scm_remember_upto_here_2 (x
, y
);
3615 return (cmp
> 0) ? x
: y
;
3617 else if (SCM_REALP (y
))
3619 /* if y==NaN then xx>yy is false, so we return the NaN y */
3622 xx
= scm_i_big2dbl (x
);
3623 yy
= SCM_REAL_VALUE (y
);
3624 return (xx
> yy
? scm_from_double (xx
) : y
);
3626 else if (SCM_FRACTIONP (y
))
3631 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3633 else if (SCM_REALP (x
))
3635 if (SCM_I_INUMP (y
))
3637 double z
= SCM_I_INUM (y
);
3638 /* if x==NaN then "<" is false and we return NaN */
3639 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3641 else if (SCM_BIGP (y
))
3646 else if (SCM_REALP (y
))
3648 /* if x==NaN then our explicit check means we return NaN
3649 if y==NaN then ">" is false and we return NaN
3650 calling isnan is unavoidable, since it's the only way to know
3651 which of x or y causes any compares to be false */
3652 double xx
= SCM_REAL_VALUE (x
);
3653 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3655 else if (SCM_FRACTIONP (y
))
3657 double yy
= scm_i_fraction2double (y
);
3658 double xx
= SCM_REAL_VALUE (x
);
3659 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3662 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3664 else if (SCM_FRACTIONP (x
))
3666 if (SCM_I_INUMP (y
))
3670 else if (SCM_BIGP (y
))
3674 else if (SCM_REALP (y
))
3676 double xx
= scm_i_fraction2double (x
);
3677 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3679 else if (SCM_FRACTIONP (y
))
3684 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3687 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3691 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3692 /* "Return the minium of all parameter values."
3695 scm_min (SCM x
, SCM y
)
3700 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3701 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3704 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3707 if (SCM_I_INUMP (x
))
3709 long xx
= SCM_I_INUM (x
);
3710 if (SCM_I_INUMP (y
))
3712 long yy
= SCM_I_INUM (y
);
3713 return (xx
< yy
) ? x
: y
;
3715 else if (SCM_BIGP (y
))
3717 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3718 scm_remember_upto_here_1 (y
);
3719 return (sgn
< 0) ? y
: x
;
3721 else if (SCM_REALP (y
))
3724 /* if y==NaN then "<" is false and we return NaN */
3725 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3727 else if (SCM_FRACTIONP (y
))
3730 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3733 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3735 else if (SCM_BIGP (x
))
3737 if (SCM_I_INUMP (y
))
3739 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3740 scm_remember_upto_here_1 (x
);
3741 return (sgn
< 0) ? x
: y
;
3743 else if (SCM_BIGP (y
))
3745 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3746 scm_remember_upto_here_2 (x
, y
);
3747 return (cmp
> 0) ? y
: x
;
3749 else if (SCM_REALP (y
))
3751 /* if y==NaN then xx<yy is false, so we return the NaN y */
3754 xx
= scm_i_big2dbl (x
);
3755 yy
= SCM_REAL_VALUE (y
);
3756 return (xx
< yy
? scm_from_double (xx
) : y
);
3758 else if (SCM_FRACTIONP (y
))
3763 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3765 else if (SCM_REALP (x
))
3767 if (SCM_I_INUMP (y
))
3769 double z
= SCM_I_INUM (y
);
3770 /* if x==NaN then "<" is false and we return NaN */
3771 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3773 else if (SCM_BIGP (y
))
3778 else if (SCM_REALP (y
))
3780 /* if x==NaN then our explicit check means we return NaN
3781 if y==NaN then "<" is false and we return NaN
3782 calling isnan is unavoidable, since it's the only way to know
3783 which of x or y causes any compares to be false */
3784 double xx
= SCM_REAL_VALUE (x
);
3785 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3787 else if (SCM_FRACTIONP (y
))
3789 double yy
= scm_i_fraction2double (y
);
3790 double xx
= SCM_REAL_VALUE (x
);
3791 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3794 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3796 else if (SCM_FRACTIONP (x
))
3798 if (SCM_I_INUMP (y
))
3802 else if (SCM_BIGP (y
))
3806 else if (SCM_REALP (y
))
3808 double xx
= scm_i_fraction2double (x
);
3809 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3811 else if (SCM_FRACTIONP (y
))
3816 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3819 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3823 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3824 /* "Return the sum of all parameter values. Return 0 if called without\n"
3828 scm_sum (SCM x
, SCM y
)
3832 if (SCM_NUMBERP (x
)) return x
;
3833 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3834 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3837 if (SCM_I_INUMP (x
))
3839 if (SCM_I_INUMP (y
))
3841 long xx
= SCM_I_INUM (x
);
3842 long yy
= SCM_I_INUM (y
);
3843 long int z
= xx
+ yy
;
3844 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3846 else if (SCM_BIGP (y
))
3851 else if (SCM_REALP (y
))
3853 long int xx
= SCM_I_INUM (x
);
3854 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3856 else if (SCM_COMPLEXP (y
))
3858 long int xx
= SCM_I_INUM (x
);
3859 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3860 SCM_COMPLEX_IMAG (y
));
3862 else if (SCM_FRACTIONP (y
))
3863 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3864 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3865 SCM_FRACTION_DENOMINATOR (y
));
3867 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3868 } else if (SCM_BIGP (x
))
3870 if (SCM_I_INUMP (y
))
3875 inum
= SCM_I_INUM (y
);
3878 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3881 SCM result
= scm_i_mkbig ();
3882 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3883 scm_remember_upto_here_1 (x
);
3884 /* we know the result will have to be a bignum */
3887 return scm_i_normbig (result
);
3891 SCM result
= scm_i_mkbig ();
3892 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3893 scm_remember_upto_here_1 (x
);
3894 /* we know the result will have to be a bignum */
3897 return scm_i_normbig (result
);
3900 else if (SCM_BIGP (y
))
3902 SCM result
= scm_i_mkbig ();
3903 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3904 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3905 mpz_add (SCM_I_BIG_MPZ (result
),
3908 scm_remember_upto_here_2 (x
, y
);
3909 /* we know the result will have to be a bignum */
3912 return scm_i_normbig (result
);
3914 else if (SCM_REALP (y
))
3916 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3917 scm_remember_upto_here_1 (x
);
3918 return scm_from_double (result
);
3920 else if (SCM_COMPLEXP (y
))
3922 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3923 + SCM_COMPLEX_REAL (y
));
3924 scm_remember_upto_here_1 (x
);
3925 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
3927 else if (SCM_FRACTIONP (y
))
3928 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3929 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3930 SCM_FRACTION_DENOMINATOR (y
));
3932 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3934 else if (SCM_REALP (x
))
3936 if (SCM_I_INUMP (y
))
3937 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
3938 else if (SCM_BIGP (y
))
3940 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3941 scm_remember_upto_here_1 (y
);
3942 return scm_from_double (result
);
3944 else if (SCM_REALP (y
))
3945 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3946 else if (SCM_COMPLEXP (y
))
3947 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3948 SCM_COMPLEX_IMAG (y
));
3949 else if (SCM_FRACTIONP (y
))
3950 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3952 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3954 else if (SCM_COMPLEXP (x
))
3956 if (SCM_I_INUMP (y
))
3957 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
3958 SCM_COMPLEX_IMAG (x
));
3959 else if (SCM_BIGP (y
))
3961 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3962 + SCM_COMPLEX_REAL (x
));
3963 scm_remember_upto_here_1 (y
);
3964 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
3966 else if (SCM_REALP (y
))
3967 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3968 SCM_COMPLEX_IMAG (x
));
3969 else if (SCM_COMPLEXP (y
))
3970 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3971 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3972 else if (SCM_FRACTIONP (y
))
3973 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3974 SCM_COMPLEX_IMAG (x
));
3976 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3978 else if (SCM_FRACTIONP (x
))
3980 if (SCM_I_INUMP (y
))
3981 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3982 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3983 SCM_FRACTION_DENOMINATOR (x
));
3984 else if (SCM_BIGP (y
))
3985 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3986 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3987 SCM_FRACTION_DENOMINATOR (x
));
3988 else if (SCM_REALP (y
))
3989 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
3990 else if (SCM_COMPLEXP (y
))
3991 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
3992 SCM_COMPLEX_IMAG (y
));
3993 else if (SCM_FRACTIONP (y
))
3994 /* a/b + c/d = (ad + bc) / bd */
3995 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
3996 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
3997 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
3999 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4002 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4006 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4007 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4008 * the sum of all but the first argument are subtracted from the first
4010 #define FUNC_NAME s_difference
4012 scm_difference (SCM x
, SCM y
)
4017 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4019 if (SCM_I_INUMP (x
))
4021 long xx
= -SCM_I_INUM (x
);
4022 if (SCM_FIXABLE (xx
))
4023 return SCM_I_MAKINUM (xx
);
4025 return scm_i_long2big (xx
);
4027 else if (SCM_BIGP (x
))
4028 /* FIXME: do we really need to normalize here? */
4029 return scm_i_normbig (scm_i_clonebig (x
, 0));
4030 else if (SCM_REALP (x
))
4031 return scm_from_double (-SCM_REAL_VALUE (x
));
4032 else if (SCM_COMPLEXP (x
))
4033 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4034 -SCM_COMPLEX_IMAG (x
));
4035 else if (SCM_FRACTIONP (x
))
4036 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4037 SCM_FRACTION_DENOMINATOR (x
));
4039 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4042 if (SCM_I_INUMP (x
))
4044 if (SCM_I_INUMP (y
))
4046 long int xx
= SCM_I_INUM (x
);
4047 long int yy
= SCM_I_INUM (y
);
4048 long int z
= xx
- yy
;
4049 if (SCM_FIXABLE (z
))
4050 return SCM_I_MAKINUM (z
);
4052 return scm_i_long2big (z
);
4054 else if (SCM_BIGP (y
))
4056 /* inum-x - big-y */
4057 long xx
= SCM_I_INUM (x
);
4060 return scm_i_clonebig (y
, 0);
4063 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4064 SCM result
= scm_i_mkbig ();
4067 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4070 /* x - y == -(y + -x) */
4071 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4072 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4074 scm_remember_upto_here_1 (y
);
4076 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4077 /* we know the result will have to be a bignum */
4080 return scm_i_normbig (result
);
4083 else if (SCM_REALP (y
))
4085 long int xx
= SCM_I_INUM (x
);
4086 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4088 else if (SCM_COMPLEXP (y
))
4090 long int xx
= SCM_I_INUM (x
);
4091 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4092 - SCM_COMPLEX_IMAG (y
));
4094 else if (SCM_FRACTIONP (y
))
4095 /* a - b/c = (ac - b) / c */
4096 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4097 SCM_FRACTION_NUMERATOR (y
)),
4098 SCM_FRACTION_DENOMINATOR (y
));
4100 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4102 else if (SCM_BIGP (x
))
4104 if (SCM_I_INUMP (y
))
4106 /* big-x - inum-y */
4107 long yy
= SCM_I_INUM (y
);
4108 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4110 scm_remember_upto_here_1 (x
);
4112 return (SCM_FIXABLE (-yy
) ?
4113 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4116 SCM result
= scm_i_mkbig ();
4119 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4121 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4122 scm_remember_upto_here_1 (x
);
4124 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4125 /* we know the result will have to be a bignum */
4128 return scm_i_normbig (result
);
4131 else if (SCM_BIGP (y
))
4133 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4134 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4135 SCM result
= scm_i_mkbig ();
4136 mpz_sub (SCM_I_BIG_MPZ (result
),
4139 scm_remember_upto_here_2 (x
, y
);
4140 /* we know the result will have to be a bignum */
4141 if ((sgn_x
== 1) && (sgn_y
== -1))
4143 if ((sgn_x
== -1) && (sgn_y
== 1))
4145 return scm_i_normbig (result
);
4147 else if (SCM_REALP (y
))
4149 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4150 scm_remember_upto_here_1 (x
);
4151 return scm_from_double (result
);
4153 else if (SCM_COMPLEXP (y
))
4155 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4156 - SCM_COMPLEX_REAL (y
));
4157 scm_remember_upto_here_1 (x
);
4158 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4160 else if (SCM_FRACTIONP (y
))
4161 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4162 SCM_FRACTION_NUMERATOR (y
)),
4163 SCM_FRACTION_DENOMINATOR (y
));
4164 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4166 else if (SCM_REALP (x
))
4168 if (SCM_I_INUMP (y
))
4169 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4170 else if (SCM_BIGP (y
))
4172 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4173 scm_remember_upto_here_1 (x
);
4174 return scm_from_double (result
);
4176 else if (SCM_REALP (y
))
4177 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4178 else if (SCM_COMPLEXP (y
))
4179 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4180 -SCM_COMPLEX_IMAG (y
));
4181 else if (SCM_FRACTIONP (y
))
4182 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4184 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4186 else if (SCM_COMPLEXP (x
))
4188 if (SCM_I_INUMP (y
))
4189 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4190 SCM_COMPLEX_IMAG (x
));
4191 else if (SCM_BIGP (y
))
4193 double real_part
= (SCM_COMPLEX_REAL (x
)
4194 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4195 scm_remember_upto_here_1 (x
);
4196 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4198 else if (SCM_REALP (y
))
4199 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4200 SCM_COMPLEX_IMAG (x
));
4201 else if (SCM_COMPLEXP (y
))
4202 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4203 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4204 else if (SCM_FRACTIONP (y
))
4205 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4206 SCM_COMPLEX_IMAG (x
));
4208 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4210 else if (SCM_FRACTIONP (x
))
4212 if (SCM_I_INUMP (y
))
4213 /* a/b - c = (a - cb) / b */
4214 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4215 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4216 SCM_FRACTION_DENOMINATOR (x
));
4217 else if (SCM_BIGP (y
))
4218 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4219 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4220 SCM_FRACTION_DENOMINATOR (x
));
4221 else if (SCM_REALP (y
))
4222 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4223 else if (SCM_COMPLEXP (y
))
4224 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4225 -SCM_COMPLEX_IMAG (y
));
4226 else if (SCM_FRACTIONP (y
))
4227 /* a/b - c/d = (ad - bc) / bd */
4228 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4229 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4230 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4232 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4235 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4240 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4241 /* "Return the product of all arguments. If called without arguments,\n"
4245 scm_product (SCM x
, SCM y
)
4250 return SCM_I_MAKINUM (1L);
4251 else if (SCM_NUMBERP (x
))
4254 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4257 if (SCM_I_INUMP (x
))
4262 xx
= SCM_I_INUM (x
);
4266 case 0: return x
; break;
4267 case 1: return y
; break;
4270 if (SCM_I_INUMP (y
))
4272 long yy
= SCM_I_INUM (y
);
4274 SCM k
= SCM_I_MAKINUM (kk
);
4275 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4279 SCM result
= scm_i_long2big (xx
);
4280 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4281 return scm_i_normbig (result
);
4284 else if (SCM_BIGP (y
))
4286 SCM result
= scm_i_mkbig ();
4287 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4288 scm_remember_upto_here_1 (y
);
4291 else if (SCM_REALP (y
))
4292 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4293 else if (SCM_COMPLEXP (y
))
4294 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4295 xx
* SCM_COMPLEX_IMAG (y
));
4296 else if (SCM_FRACTIONP (y
))
4297 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4298 SCM_FRACTION_DENOMINATOR (y
));
4300 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4302 else if (SCM_BIGP (x
))
4304 if (SCM_I_INUMP (y
))
4309 else if (SCM_BIGP (y
))
4311 SCM result
= scm_i_mkbig ();
4312 mpz_mul (SCM_I_BIG_MPZ (result
),
4315 scm_remember_upto_here_2 (x
, y
);
4318 else if (SCM_REALP (y
))
4320 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4321 scm_remember_upto_here_1 (x
);
4322 return scm_from_double (result
);
4324 else if (SCM_COMPLEXP (y
))
4326 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4327 scm_remember_upto_here_1 (x
);
4328 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4329 z
* SCM_COMPLEX_IMAG (y
));
4331 else if (SCM_FRACTIONP (y
))
4332 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4333 SCM_FRACTION_DENOMINATOR (y
));
4335 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4337 else if (SCM_REALP (x
))
4339 if (SCM_I_INUMP (y
))
4340 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4341 else if (SCM_BIGP (y
))
4343 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4344 scm_remember_upto_here_1 (y
);
4345 return scm_from_double (result
);
4347 else if (SCM_REALP (y
))
4348 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4349 else if (SCM_COMPLEXP (y
))
4350 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4351 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4352 else if (SCM_FRACTIONP (y
))
4353 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4355 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4357 else if (SCM_COMPLEXP (x
))
4359 if (SCM_I_INUMP (y
))
4360 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4361 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4362 else if (SCM_BIGP (y
))
4364 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4365 scm_remember_upto_here_1 (y
);
4366 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4367 z
* SCM_COMPLEX_IMAG (x
));
4369 else if (SCM_REALP (y
))
4370 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4371 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4372 else if (SCM_COMPLEXP (y
))
4374 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4375 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4376 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4377 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4379 else if (SCM_FRACTIONP (y
))
4381 double yy
= scm_i_fraction2double (y
);
4382 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4383 yy
* SCM_COMPLEX_IMAG (x
));
4386 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4388 else if (SCM_FRACTIONP (x
))
4390 if (SCM_I_INUMP (y
))
4391 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4392 SCM_FRACTION_DENOMINATOR (x
));
4393 else if (SCM_BIGP (y
))
4394 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4395 SCM_FRACTION_DENOMINATOR (x
));
4396 else if (SCM_REALP (y
))
4397 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4398 else if (SCM_COMPLEXP (y
))
4400 double xx
= scm_i_fraction2double (x
);
4401 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4402 xx
* SCM_COMPLEX_IMAG (y
));
4404 else if (SCM_FRACTIONP (y
))
4405 /* a/b * c/d = ac / bd */
4406 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4407 SCM_FRACTION_NUMERATOR (y
)),
4408 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4409 SCM_FRACTION_DENOMINATOR (y
)));
4411 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4414 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4417 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4418 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4419 #define ALLOW_DIVIDE_BY_ZERO
4420 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4423 /* The code below for complex division is adapted from the GNU
4424 libstdc++, which adapted it from f2c's libF77, and is subject to
4427 /****************************************************************
4428 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4430 Permission to use, copy, modify, and distribute this software
4431 and its documentation for any purpose and without fee is hereby
4432 granted, provided that the above copyright notice appear in all
4433 copies and that both that the copyright notice and this
4434 permission notice and warranty disclaimer appear in supporting
4435 documentation, and that the names of AT&T Bell Laboratories or
4436 Bellcore or any of their entities not be used in advertising or
4437 publicity pertaining to distribution of the software without
4438 specific, written prior permission.
4440 AT&T and Bellcore disclaim all warranties with regard to this
4441 software, including all implied warranties of merchantability
4442 and fitness. In no event shall AT&T or Bellcore be liable for
4443 any special, indirect or consequential damages or any damages
4444 whatsoever resulting from loss of use, data or profits, whether
4445 in an action of contract, negligence or other tortious action,
4446 arising out of or in connection with the use or performance of
4448 ****************************************************************/
4450 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4451 /* Divide the first argument by the product of the remaining
4452 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4454 #define FUNC_NAME s_divide
4456 scm_i_divide (SCM x
, SCM y
, int inexact
)
4463 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4464 else if (SCM_I_INUMP (x
))
4466 long xx
= SCM_I_INUM (x
);
4467 if (xx
== 1 || xx
== -1)
4469 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4471 scm_num_overflow (s_divide
);
4476 return scm_from_double (1.0 / (double) xx
);
4477 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4480 else if (SCM_BIGP (x
))
4483 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4484 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4486 else if (SCM_REALP (x
))
4488 double xx
= SCM_REAL_VALUE (x
);
4489 #ifndef ALLOW_DIVIDE_BY_ZERO
4491 scm_num_overflow (s_divide
);
4494 return scm_from_double (1.0 / xx
);
4496 else if (SCM_COMPLEXP (x
))
4498 double r
= SCM_COMPLEX_REAL (x
);
4499 double i
= SCM_COMPLEX_IMAG (x
);
4503 double d
= i
* (1.0 + t
* t
);
4504 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4509 double d
= r
* (1.0 + t
* t
);
4510 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4513 else if (SCM_FRACTIONP (x
))
4514 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4515 SCM_FRACTION_NUMERATOR (x
));
4517 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4520 if (SCM_I_INUMP (x
))
4522 long xx
= SCM_I_INUM (x
);
4523 if (SCM_I_INUMP (y
))
4525 long yy
= SCM_I_INUM (y
);
4528 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4529 scm_num_overflow (s_divide
);
4531 return scm_from_double ((double) xx
/ (double) yy
);
4534 else if (xx
% yy
!= 0)
4537 return scm_from_double ((double) xx
/ (double) yy
);
4538 else return scm_i_make_ratio (x
, y
);
4543 if (SCM_FIXABLE (z
))
4544 return SCM_I_MAKINUM (z
);
4546 return scm_i_long2big (z
);
4549 else if (SCM_BIGP (y
))
4552 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4553 else return scm_i_make_ratio (x
, y
);
4555 else if (SCM_REALP (y
))
4557 double yy
= SCM_REAL_VALUE (y
);
4558 #ifndef ALLOW_DIVIDE_BY_ZERO
4560 scm_num_overflow (s_divide
);
4563 return scm_from_double ((double) xx
/ yy
);
4565 else if (SCM_COMPLEXP (y
))
4568 complex_div
: /* y _must_ be a complex number */
4570 double r
= SCM_COMPLEX_REAL (y
);
4571 double i
= SCM_COMPLEX_IMAG (y
);
4575 double d
= i
* (1.0 + t
* t
);
4576 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4581 double d
= r
* (1.0 + t
* t
);
4582 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4586 else if (SCM_FRACTIONP (y
))
4587 /* a / b/c = ac / b */
4588 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4589 SCM_FRACTION_NUMERATOR (y
));
4591 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4593 else if (SCM_BIGP (x
))
4595 if (SCM_I_INUMP (y
))
4597 long int yy
= SCM_I_INUM (y
);
4600 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4601 scm_num_overflow (s_divide
);
4603 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4604 scm_remember_upto_here_1 (x
);
4605 return (sgn
== 0) ? scm_nan () : scm_inf ();
4612 /* FIXME: HMM, what are the relative performance issues here?
4613 We need to test. Is it faster on average to test
4614 divisible_p, then perform whichever operation, or is it
4615 faster to perform the integer div opportunistically and
4616 switch to real if there's a remainder? For now we take the
4617 middle ground: test, then if divisible, use the faster div
4620 long abs_yy
= yy
< 0 ? -yy
: yy
;
4621 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4625 SCM result
= scm_i_mkbig ();
4626 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4627 scm_remember_upto_here_1 (x
);
4629 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4630 return scm_i_normbig (result
);
4635 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4636 else return scm_i_make_ratio (x
, y
);
4640 else if (SCM_BIGP (y
))
4642 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4645 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4646 scm_num_overflow (s_divide
);
4648 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4649 scm_remember_upto_here_1 (x
);
4650 return (sgn
== 0) ? scm_nan () : scm_inf ();
4656 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4660 SCM result
= scm_i_mkbig ();
4661 mpz_divexact (SCM_I_BIG_MPZ (result
),
4664 scm_remember_upto_here_2 (x
, y
);
4665 return scm_i_normbig (result
);
4671 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4672 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4673 scm_remember_upto_here_2 (x
, y
);
4674 return scm_from_double (dbx
/ dby
);
4676 else return scm_i_make_ratio (x
, y
);
4680 else if (SCM_REALP (y
))
4682 double yy
= SCM_REAL_VALUE (y
);
4683 #ifndef ALLOW_DIVIDE_BY_ZERO
4685 scm_num_overflow (s_divide
);
4688 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4690 else if (SCM_COMPLEXP (y
))
4692 a
= scm_i_big2dbl (x
);
4695 else if (SCM_FRACTIONP (y
))
4696 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4697 SCM_FRACTION_NUMERATOR (y
));
4699 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4701 else if (SCM_REALP (x
))
4703 double rx
= SCM_REAL_VALUE (x
);
4704 if (SCM_I_INUMP (y
))
4706 long int yy
= SCM_I_INUM (y
);
4707 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4709 scm_num_overflow (s_divide
);
4712 return scm_from_double (rx
/ (double) yy
);
4714 else if (SCM_BIGP (y
))
4716 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4717 scm_remember_upto_here_1 (y
);
4718 return scm_from_double (rx
/ dby
);
4720 else if (SCM_REALP (y
))
4722 double yy
= SCM_REAL_VALUE (y
);
4723 #ifndef ALLOW_DIVIDE_BY_ZERO
4725 scm_num_overflow (s_divide
);
4728 return scm_from_double (rx
/ yy
);
4730 else if (SCM_COMPLEXP (y
))
4735 else if (SCM_FRACTIONP (y
))
4736 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4738 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4740 else if (SCM_COMPLEXP (x
))
4742 double rx
= SCM_COMPLEX_REAL (x
);
4743 double ix
= SCM_COMPLEX_IMAG (x
);
4744 if (SCM_I_INUMP (y
))
4746 long int yy
= SCM_I_INUM (y
);
4747 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4749 scm_num_overflow (s_divide
);
4754 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4757 else if (SCM_BIGP (y
))
4759 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4760 scm_remember_upto_here_1 (y
);
4761 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4763 else if (SCM_REALP (y
))
4765 double yy
= SCM_REAL_VALUE (y
);
4766 #ifndef ALLOW_DIVIDE_BY_ZERO
4768 scm_num_overflow (s_divide
);
4771 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4773 else if (SCM_COMPLEXP (y
))
4775 double ry
= SCM_COMPLEX_REAL (y
);
4776 double iy
= SCM_COMPLEX_IMAG (y
);
4780 double d
= iy
* (1.0 + t
* t
);
4781 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4786 double d
= ry
* (1.0 + t
* t
);
4787 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4790 else if (SCM_FRACTIONP (y
))
4792 double yy
= scm_i_fraction2double (y
);
4793 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4796 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4798 else if (SCM_FRACTIONP (x
))
4800 if (SCM_I_INUMP (y
))
4802 long int yy
= SCM_I_INUM (y
);
4803 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4805 scm_num_overflow (s_divide
);
4808 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4809 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4811 else if (SCM_BIGP (y
))
4813 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4814 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4816 else if (SCM_REALP (y
))
4818 double yy
= SCM_REAL_VALUE (y
);
4819 #ifndef ALLOW_DIVIDE_BY_ZERO
4821 scm_num_overflow (s_divide
);
4824 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4826 else if (SCM_COMPLEXP (y
))
4828 a
= scm_i_fraction2double (x
);
4831 else if (SCM_FRACTIONP (y
))
4832 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4833 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4835 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4838 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4842 scm_divide (SCM x
, SCM y
)
4844 return scm_i_divide (x
, y
, 0);
4847 static SCM
scm_divide2real (SCM x
, SCM y
)
4849 return scm_i_divide (x
, y
, 1);
4855 scm_asinh (double x
)
4860 #define asinh scm_asinh
4861 return log (x
+ sqrt (x
* x
+ 1));
4864 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4865 /* "Return the inverse hyperbolic sine of @var{x}."
4870 scm_acosh (double x
)
4875 #define acosh scm_acosh
4876 return log (x
+ sqrt (x
* x
- 1));
4879 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4880 /* "Return the inverse hyperbolic cosine of @var{x}."
4885 scm_atanh (double x
)
4890 #define atanh scm_atanh
4891 return 0.5 * log ((1 + x
) / (1 - x
));
4894 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4895 /* "Return the inverse hyperbolic tangent of @var{x}."
4900 scm_c_truncate (double x
)
4911 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
4912 half-way case (ie. when x is an integer plus 0.5) going upwards.
4913 Then half-way cases are identified and adjusted down if the
4914 round-upwards didn't give the desired even integer.
4916 "plus_half == result" identifies a half-way case. If plus_half, which is
4917 x + 0.5, is an integer then x must be an integer plus 0.5.
4919 An odd "result" value is identified with result/2 != floor(result/2).
4920 This is done with plus_half, since that value is ready for use sooner in
4921 a pipelined cpu, and we're already requiring plus_half == result.
4923 Note however that we need to be careful when x is big and already an
4924 integer. In that case "x+0.5" may round to an adjacent integer, causing
4925 us to return such a value, incorrectly. For instance if the hardware is
4926 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
4927 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
4928 returned. Or if the hardware is in round-upwards mode, then other bigger
4929 values like say x == 2^128 will see x+0.5 rounding up to the next higher
4930 representable value, 2^128+2^76 (or whatever), again incorrect.
4932 These bad roundings of x+0.5 are avoided by testing at the start whether
4933 x is already an integer. If it is then clearly that's the desired result
4934 already. And if it's not then the exponent must be small enough to allow
4935 an 0.5 to be represented, and hence added without a bad rounding. */
4938 scm_c_round (double x
)
4940 double plus_half
, result
;
4945 plus_half
= x
+ 0.5;
4946 result
= floor (plus_half
);
4947 /* Adjust so that the rounding is towards even. */
4948 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4953 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4955 "Round the number @var{x} towards zero.")
4956 #define FUNC_NAME s_scm_truncate_number
4958 if (scm_is_false (scm_negative_p (x
)))
4959 return scm_floor (x
);
4961 return scm_ceiling (x
);
4965 static SCM exactly_one_half
;
4967 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4969 "Round the number @var{x} towards the nearest integer. "
4970 "When it is exactly halfway between two integers, "
4971 "round towards the even one.")
4972 #define FUNC_NAME s_scm_round_number
4974 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4976 else if (SCM_REALP (x
))
4977 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
4980 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
4981 single quotient+remainder division then examining to see which way
4982 the rounding should go. */
4983 SCM plus_half
= scm_sum (x
, exactly_one_half
);
4984 SCM result
= scm_floor (plus_half
);
4985 /* Adjust so that the rounding is towards even. */
4986 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
4987 && scm_is_true (scm_odd_p (result
)))
4988 return scm_difference (result
, SCM_I_MAKINUM (1));
4995 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
4997 "Round the number @var{x} towards minus infinity.")
4998 #define FUNC_NAME s_scm_floor
5000 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5002 else if (SCM_REALP (x
))
5003 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5004 else if (SCM_FRACTIONP (x
))
5006 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5007 SCM_FRACTION_DENOMINATOR (x
));
5008 if (scm_is_false (scm_negative_p (x
)))
5010 /* For positive x, rounding towards zero is correct. */
5015 /* For negative x, we need to return q-1 unless x is an
5016 integer. But fractions are never integer, per our
5018 return scm_difference (q
, SCM_I_MAKINUM (1));
5022 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5026 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5028 "Round the number @var{x} towards infinity.")
5029 #define FUNC_NAME s_scm_ceiling
5031 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5033 else if (SCM_REALP (x
))
5034 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5035 else if (SCM_FRACTIONP (x
))
5037 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5038 SCM_FRACTION_DENOMINATOR (x
));
5039 if (scm_is_false (scm_positive_p (x
)))
5041 /* For negative x, rounding towards zero is correct. */
5046 /* For positive x, we need to return q+1 unless x is an
5047 integer. But fractions are never integer, per our
5049 return scm_sum (q
, SCM_I_MAKINUM (1));
5053 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5057 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5058 /* "Return the square root of the real number @var{x}."
5060 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5061 /* "Return the absolute value of the real number @var{x}."
5063 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5064 /* "Return the @var{x}th power of e."
5066 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5067 /* "Return the natural logarithm of the real number @var{x}."
5069 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5070 /* "Return the sine of the real number @var{x}."
5072 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5073 /* "Return the cosine of the real number @var{x}."
5075 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5076 /* "Return the tangent of the real number @var{x}."
5078 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5079 /* "Return the arc sine of the real number @var{x}."
5081 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5082 /* "Return the arc cosine of the real number @var{x}."
5084 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5085 /* "Return the arc tangent of the real number @var{x}."
5087 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5088 /* "Return the hyperbolic sine of the real number @var{x}."
5090 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5091 /* "Return the hyperbolic cosine of the real number @var{x}."
5093 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5094 /* "Return the hyperbolic tangent of the real number @var{x}."
5102 static void scm_two_doubles (SCM x
,
5104 const char *sstring
,
5108 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5110 if (SCM_I_INUMP (x
))
5111 xy
->x
= SCM_I_INUM (x
);
5112 else if (SCM_BIGP (x
))
5113 xy
->x
= scm_i_big2dbl (x
);
5114 else if (SCM_REALP (x
))
5115 xy
->x
= SCM_REAL_VALUE (x
);
5116 else if (SCM_FRACTIONP (x
))
5117 xy
->x
= scm_i_fraction2double (x
);
5119 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5121 if (SCM_I_INUMP (y
))
5122 xy
->y
= SCM_I_INUM (y
);
5123 else if (SCM_BIGP (y
))
5124 xy
->y
= scm_i_big2dbl (y
);
5125 else if (SCM_REALP (y
))
5126 xy
->y
= SCM_REAL_VALUE (y
);
5127 else if (SCM_FRACTIONP (y
))
5128 xy
->y
= scm_i_fraction2double (y
);
5130 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5134 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5136 "Return @var{x} raised to the power of @var{y}. This\n"
5137 "procedure does not accept complex arguments.")
5138 #define FUNC_NAME s_scm_sys_expt
5141 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5142 return scm_from_double (pow (xy
.x
, xy
.y
));
5147 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5149 "Return the arc tangent of the two arguments @var{x} and\n"
5150 "@var{y}. This is similar to calculating the arc tangent of\n"
5151 "@var{x} / @var{y}, except that the signs of both arguments\n"
5152 "are used to determine the quadrant of the result. This\n"
5153 "procedure does not accept complex arguments.")
5154 #define FUNC_NAME s_scm_sys_atan2
5157 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5158 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5163 scm_c_make_rectangular (double re
, double im
)
5166 return scm_from_double (re
);
5170 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5172 SCM_COMPLEX_REAL (z
) = re
;
5173 SCM_COMPLEX_IMAG (z
) = im
;
5178 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5179 (SCM real
, SCM imaginary
),
5180 "Return a complex number constructed of the given @var{real} and\n"
5181 "@var{imaginary} parts.")
5182 #define FUNC_NAME s_scm_make_rectangular
5185 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5186 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5191 scm_c_make_polar (double mag
, double ang
)
5195 sincos (ang
, &s
, &c
);
5200 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5203 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5205 "Return the complex number @var{x} * e^(i * @var{y}).")
5206 #define FUNC_NAME s_scm_make_polar
5209 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5210 return scm_c_make_polar (xy
.x
, xy
.y
);
5215 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5216 /* "Return the real part of the number @var{z}."
5219 scm_real_part (SCM z
)
5221 if (SCM_I_INUMP (z
))
5223 else if (SCM_BIGP (z
))
5225 else if (SCM_REALP (z
))
5227 else if (SCM_COMPLEXP (z
))
5228 return scm_from_double (SCM_COMPLEX_REAL (z
));
5229 else if (SCM_FRACTIONP (z
))
5232 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5236 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5237 /* "Return the imaginary part of the number @var{z}."
5240 scm_imag_part (SCM z
)
5242 if (SCM_I_INUMP (z
))
5244 else if (SCM_BIGP (z
))
5246 else if (SCM_REALP (z
))
5248 else if (SCM_COMPLEXP (z
))
5249 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5250 else if (SCM_FRACTIONP (z
))
5253 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5256 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5257 /* "Return the numerator of the number @var{z}."
5260 scm_numerator (SCM z
)
5262 if (SCM_I_INUMP (z
))
5264 else if (SCM_BIGP (z
))
5266 else if (SCM_FRACTIONP (z
))
5268 scm_i_fraction_reduce (z
);
5269 return SCM_FRACTION_NUMERATOR (z
);
5271 else if (SCM_REALP (z
))
5272 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5274 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5278 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5279 /* "Return the denominator of the number @var{z}."
5282 scm_denominator (SCM z
)
5284 if (SCM_I_INUMP (z
))
5285 return SCM_I_MAKINUM (1);
5286 else if (SCM_BIGP (z
))
5287 return SCM_I_MAKINUM (1);
5288 else if (SCM_FRACTIONP (z
))
5290 scm_i_fraction_reduce (z
);
5291 return SCM_FRACTION_DENOMINATOR (z
);
5293 else if (SCM_REALP (z
))
5294 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5296 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5299 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5300 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5301 * "@code{abs} for real arguments, but also allows complex numbers."
5304 scm_magnitude (SCM z
)
5306 if (SCM_I_INUMP (z
))
5308 long int zz
= SCM_I_INUM (z
);
5311 else if (SCM_POSFIXABLE (-zz
))
5312 return SCM_I_MAKINUM (-zz
);
5314 return scm_i_long2big (-zz
);
5316 else if (SCM_BIGP (z
))
5318 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5319 scm_remember_upto_here_1 (z
);
5321 return scm_i_clonebig (z
, 0);
5325 else if (SCM_REALP (z
))
5326 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5327 else if (SCM_COMPLEXP (z
))
5328 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5329 else if (SCM_FRACTIONP (z
))
5331 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5333 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5334 SCM_FRACTION_DENOMINATOR (z
));
5337 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5341 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5342 /* "Return the angle of the complex number @var{z}."
5347 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5348 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5349 But if atan2 follows the floating point rounding mode, then the value
5350 is not a constant. Maybe it'd be close enough though. */
5351 if (SCM_I_INUMP (z
))
5353 if (SCM_I_INUM (z
) >= 0)
5356 return scm_from_double (atan2 (0.0, -1.0));
5358 else if (SCM_BIGP (z
))
5360 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5361 scm_remember_upto_here_1 (z
);
5363 return scm_from_double (atan2 (0.0, -1.0));
5367 else if (SCM_REALP (z
))
5369 if (SCM_REAL_VALUE (z
) >= 0)
5372 return scm_from_double (atan2 (0.0, -1.0));
5374 else if (SCM_COMPLEXP (z
))
5375 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5376 else if (SCM_FRACTIONP (z
))
5378 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5380 else return scm_from_double (atan2 (0.0, -1.0));
5383 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5387 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5388 /* Convert the number @var{x} to its inexact representation.\n"
5391 scm_exact_to_inexact (SCM z
)
5393 if (SCM_I_INUMP (z
))
5394 return scm_from_double ((double) SCM_I_INUM (z
));
5395 else if (SCM_BIGP (z
))
5396 return scm_from_double (scm_i_big2dbl (z
));
5397 else if (SCM_FRACTIONP (z
))
5398 return scm_from_double (scm_i_fraction2double (z
));
5399 else if (SCM_INEXACTP (z
))
5402 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5406 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5408 "Return an exact number that is numerically closest to @var{z}.")
5409 #define FUNC_NAME s_scm_inexact_to_exact
5411 if (SCM_I_INUMP (z
))
5413 else if (SCM_BIGP (z
))
5415 else if (SCM_REALP (z
))
5417 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5418 SCM_OUT_OF_RANGE (1, z
);
5425 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5426 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5427 scm_i_mpz2num (mpq_denref (frac
)));
5429 /* When scm_i_make_ratio throws, we leak the memory allocated
5436 else if (SCM_FRACTIONP (z
))
5439 SCM_WRONG_TYPE_ARG (1, z
);
5443 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5445 "Return an exact number that is within @var{err} of @var{x}.")
5446 #define FUNC_NAME s_scm_rationalize
5448 if (SCM_I_INUMP (x
))
5450 else if (SCM_BIGP (x
))
5452 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5454 /* Use continued fractions to find closest ratio. All
5455 arithmetic is done with exact numbers.
5458 SCM ex
= scm_inexact_to_exact (x
);
5459 SCM int_part
= scm_floor (ex
);
5460 SCM tt
= SCM_I_MAKINUM (1);
5461 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5462 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5466 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5469 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5470 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5472 /* We stop after a million iterations just to be absolutely sure
5473 that we don't go into an infinite loop. The process normally
5474 converges after less than a dozen iterations.
5477 err
= scm_abs (err
);
5478 while (++i
< 1000000)
5480 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5481 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5482 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5484 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5485 err
))) /* abs(x-a/b) <= err */
5487 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5488 if (scm_is_false (scm_exact_p (x
))
5489 || scm_is_false (scm_exact_p (err
)))
5490 return scm_exact_to_inexact (res
);
5494 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5496 tt
= scm_floor (rx
); /* tt = floor (rx) */
5502 scm_num_overflow (s_scm_rationalize
);
5505 SCM_WRONG_TYPE_ARG (1, x
);
5509 /* conversion functions */
5512 scm_is_integer (SCM val
)
5514 return scm_is_true (scm_integer_p (val
));
5518 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5520 if (SCM_I_INUMP (val
))
5522 scm_t_signed_bits n
= SCM_I_INUM (val
);
5523 return n
>= min
&& n
<= max
;
5525 else if (SCM_BIGP (val
))
5527 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5529 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5531 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5533 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5534 return n
>= min
&& n
<= max
;
5544 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5545 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5548 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5549 SCM_I_BIG_MPZ (val
));
5551 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5563 return n
>= min
&& n
<= max
;
5571 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5573 if (SCM_I_INUMP (val
))
5575 scm_t_signed_bits n
= SCM_I_INUM (val
);
5576 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5578 else if (SCM_BIGP (val
))
5580 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5582 else if (max
<= ULONG_MAX
)
5584 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5586 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5587 return n
>= min
&& n
<= max
;
5597 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5600 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5601 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5604 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5605 SCM_I_BIG_MPZ (val
));
5607 return n
>= min
&& n
<= max
;
5614 #define TYPE scm_t_intmax
5615 #define TYPE_MIN min
5616 #define TYPE_MAX max
5617 #define SIZEOF_TYPE 0
5618 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5619 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5620 #include "libguile/conv-integer.i.c"
5622 #define TYPE scm_t_uintmax
5623 #define TYPE_MIN min
5624 #define TYPE_MAX max
5625 #define SIZEOF_TYPE 0
5626 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5627 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5628 #include "libguile/conv-uinteger.i.c"
5630 #define TYPE scm_t_int8
5631 #define TYPE_MIN SCM_T_INT8_MIN
5632 #define TYPE_MAX SCM_T_INT8_MAX
5633 #define SIZEOF_TYPE 1
5634 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5635 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5636 #include "libguile/conv-integer.i.c"
5638 #define TYPE scm_t_uint8
5640 #define TYPE_MAX SCM_T_UINT8_MAX
5641 #define SIZEOF_TYPE 1
5642 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5643 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5644 #include "libguile/conv-uinteger.i.c"
5646 #define TYPE scm_t_int16
5647 #define TYPE_MIN SCM_T_INT16_MIN
5648 #define TYPE_MAX SCM_T_INT16_MAX
5649 #define SIZEOF_TYPE 2
5650 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5651 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5652 #include "libguile/conv-integer.i.c"
5654 #define TYPE scm_t_uint16
5656 #define TYPE_MAX SCM_T_UINT16_MAX
5657 #define SIZEOF_TYPE 2
5658 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5659 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5660 #include "libguile/conv-uinteger.i.c"
5662 #define TYPE scm_t_int32
5663 #define TYPE_MIN SCM_T_INT32_MIN
5664 #define TYPE_MAX SCM_T_INT32_MAX
5665 #define SIZEOF_TYPE 4
5666 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5667 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5668 #include "libguile/conv-integer.i.c"
5670 #define TYPE scm_t_uint32
5672 #define TYPE_MAX SCM_T_UINT32_MAX
5673 #define SIZEOF_TYPE 4
5674 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5675 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5676 #include "libguile/conv-uinteger.i.c"
5678 #if SCM_HAVE_T_INT64
5680 #define TYPE scm_t_int64
5681 #define TYPE_MIN SCM_T_INT64_MIN
5682 #define TYPE_MAX SCM_T_INT64_MAX
5683 #define SIZEOF_TYPE 8
5684 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5685 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5686 #include "libguile/conv-integer.i.c"
5688 #define TYPE scm_t_uint64
5690 #define TYPE_MAX SCM_T_UINT64_MAX
5691 #define SIZEOF_TYPE 8
5692 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5693 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5694 #include "libguile/conv-uinteger.i.c"
5699 scm_is_real (SCM val
)
5701 return scm_is_true (scm_real_p (val
));
5705 scm_is_rational (SCM val
)
5707 return scm_is_true (scm_rational_p (val
));
5711 scm_to_double (SCM val
)
5713 if (SCM_I_INUMP (val
))
5714 return SCM_I_INUM (val
);
5715 else if (SCM_BIGP (val
))
5716 return scm_i_big2dbl (val
);
5717 else if (SCM_FRACTIONP (val
))
5718 return scm_i_fraction2double (val
);
5719 else if (SCM_REALP (val
))
5720 return SCM_REAL_VALUE (val
);
5722 scm_wrong_type_arg (NULL
, 0, val
);
5726 scm_from_double (double val
)
5728 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5729 SCM_REAL_VALUE (z
) = val
;
5733 #if SCM_ENABLE_DISCOURAGED == 1
5736 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5740 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5744 scm_out_of_range (NULL
, num
);
5747 return scm_to_double (num
);
5751 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5755 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5759 scm_out_of_range (NULL
, num
);
5762 return scm_to_double (num
);
5768 scm_is_complex (SCM val
)
5770 return scm_is_true (scm_complex_p (val
));
5774 scm_c_real_part (SCM z
)
5776 if (SCM_COMPLEXP (z
))
5777 return SCM_COMPLEX_REAL (z
);
5780 /* Use the scm_real_part to get proper error checking and
5783 return scm_to_double (scm_real_part (z
));
5788 scm_c_imag_part (SCM z
)
5790 if (SCM_COMPLEXP (z
))
5791 return SCM_COMPLEX_IMAG (z
);
5794 /* Use the scm_imag_part to get proper error checking and
5795 dispatching. The result will almost always be 0.0, but not
5798 return scm_to_double (scm_imag_part (z
));
5803 scm_c_magnitude (SCM z
)
5805 return scm_to_double (scm_magnitude (z
));
5811 return scm_to_double (scm_angle (z
));
5815 scm_is_number (SCM z
)
5817 return scm_is_true (scm_number_p (z
));
5825 mpz_init_set_si (z_negative_one
, -1);
5827 /* It may be possible to tune the performance of some algorithms by using
5828 * the following constants to avoid the creation of bignums. Please, before
5829 * using these values, remember the two rules of program optimization:
5830 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5831 scm_c_define ("most-positive-fixnum",
5832 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5833 scm_c_define ("most-negative-fixnum",
5834 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5836 scm_add_feature ("complex");
5837 scm_add_feature ("inexact");
5838 scm_flo0
= scm_from_double (0.0);
5840 /* determine floating point precision */
5841 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
5843 init_dblprec(&scm_dblprec
[i
-2],i
);
5844 init_fx_radix(fx_per_radix
[i
-2],i
);
5847 /* hard code precision for base 10 if the preprocessor tells us to... */
5848 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5855 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
5856 SCM_I_MAKINUM (2)));
5857 #include "libguile/numbers.x"