1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004 Free Software Foundation, Inc.
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 /* General assumptions:
24 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
25 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
26 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
27 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
28 * All objects satisfying SCM_FRACTIONP are never an integer.
33 - see if special casing bignums and reals in integer-exponent when
34 possible (to use mpz_pow and mpf_pow_ui) is faster.
36 - look in to better short-circuiting of common cases in
37 integer-expt and elsewhere.
39 - see if direct mpz operations can help in ash and elsewhere.
43 /* tell glibc (2.3) to give prototype for C99 trunc() */
55 #include "libguile/_scm.h"
56 #include "libguile/feature.h"
57 #include "libguile/ports.h"
58 #include "libguile/root.h"
59 #include "libguile/smob.h"
60 #include "libguile/strings.h"
62 #include "libguile/validate.h"
63 #include "libguile/numbers.h"
64 #include "libguile/deprecation.h"
66 #include "libguile/eq.h"
68 #include "libguile/discouraged.h"
73 Wonder if this might be faster for some of our code? A switch on
74 the numtag would jump directly to the right case, and the
75 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
77 #define SCM_I_NUMTAG_NOTNUM 0
78 #define SCM_I_NUMTAG_INUM 1
79 #define SCM_I_NUMTAG_BIG scm_tc16_big
80 #define SCM_I_NUMTAG_REAL scm_tc16_real
81 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
82 #define SCM_I_NUMTAG(x) \
83 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
84 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
85 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
86 : SCM_I_NUMTAG_NOTNUM)))
88 /* the macro above will not work as is with fractions */
91 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
93 /* FLOBUFLEN is the maximum number of characters neccessary for the
94 * printed or scm_string representation of an inexact number.
96 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
99 #if ! defined (HAVE_ISNAN)
104 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
107 #if ! defined (HAVE_ISINF)
112 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
119 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
120 an explicit check. In some future gmp (don't know what version number),
121 mpz_cmp_d is supposed to do this itself. */
123 #define xmpz_cmp_d(z, d) \
124 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
126 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
129 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
130 isinf. It does have finite and isnan though, hence the use of those.
131 fpclass would be a possibility on that system too. */
135 #if defined (HAVE_ISINF)
137 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
138 return (! (finite (x
) || isnan (x
)));
147 #if defined (HAVE_ISNAN)
156 static mpz_t z_negative_one
;
160 SCM_C_INLINE_KEYWORD SCM
163 /* Return a newly created bignum. */
164 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
165 mpz_init (SCM_I_BIG_MPZ (z
));
169 SCM_C_INLINE_KEYWORD SCM
170 scm_i_long2big (long x
)
172 /* Return a newly created bignum initialized to X. */
173 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
174 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
178 SCM_C_INLINE_KEYWORD SCM
179 scm_i_ulong2big (unsigned long x
)
181 /* Return a newly created bignum initialized to X. */
182 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
183 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
187 SCM_C_INLINE_KEYWORD
static SCM
188 scm_i_clonebig (SCM src_big
, int same_sign_p
)
190 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
191 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
192 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
194 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
198 SCM_C_INLINE_KEYWORD
int
199 scm_i_bigcmp (SCM x
, SCM y
)
201 /* Return neg if x < y, pos if x > y, and 0 if x == y */
202 /* presume we already know x and y are bignums */
203 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
204 scm_remember_upto_here_2 (x
, y
);
208 SCM_C_INLINE_KEYWORD SCM
209 scm_i_dbl2big (double d
)
211 /* results are only defined if d is an integer */
212 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
213 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
217 /* Convert a integer in double representation to a SCM number. */
219 SCM_C_INLINE_KEYWORD SCM
220 scm_i_dbl2num (double u
)
222 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
223 powers of 2, so there's no rounding when making "double" values
224 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
225 get rounded on a 64-bit machine, hence the "+1".
227 The use of floor() to force to an integer value ensures we get a
228 "numerically closest" value without depending on how a
229 double->long cast or how mpz_set_d will round. For reference,
230 double->long probably follows the hardware rounding mode,
231 mpz_set_d truncates towards zero. */
233 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
234 representable as a double? */
236 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
237 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
238 return SCM_I_MAKINUM ((long) u
);
240 return scm_i_dbl2big (u
);
243 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
244 with R5RS exact->inexact.
246 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
247 (ie. truncate towards zero), then adjust to get the closest double by
248 examining the next lower bit and adding 1 (to the absolute value) if
251 Bignums exactly half way between representable doubles are rounded to the
252 next higher absolute value (ie. away from zero). This seems like an
253 adequate interpretation of R5RS "numerically closest", and it's easier
254 and faster than a full "nearest-even" style.
256 The bit test must be done on the absolute value of the mpz_t, which means
257 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
258 negatives as twos complement.
260 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
261 following the hardware rounding mode, but applied to the absolute value
262 of the mpz_t operand. This is not what we want so we put the high
263 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
264 mpz_get_d is supposed to always truncate towards zero.
266 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
267 is a slowdown. It'd be faster to pick out the relevant high bits with
268 mpz_getlimbn if we could be bothered coding that, and if the new
269 truncating gmp doesn't come out. */
272 scm_i_big2dbl (SCM b
)
277 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
281 /* Current GMP, eg. 4.1.3, force truncation towards zero */
283 if (bits
> DBL_MANT_DIG
)
285 size_t shift
= bits
- DBL_MANT_DIG
;
286 mpz_init2 (tmp
, DBL_MANT_DIG
);
287 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
288 result
= ldexp (mpz_get_d (tmp
), shift
);
293 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
298 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
301 if (bits
> DBL_MANT_DIG
)
303 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
304 /* test bit number "pos" in absolute value */
305 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
306 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
308 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
312 scm_remember_upto_here_1 (b
);
316 SCM_C_INLINE_KEYWORD SCM
317 scm_i_normbig (SCM b
)
319 /* convert a big back to a fixnum if it'll fit */
320 /* presume b is a bignum */
321 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
323 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
324 if (SCM_FIXABLE (val
))
325 b
= SCM_I_MAKINUM (val
);
330 static SCM_C_INLINE_KEYWORD SCM
331 scm_i_mpz2num (mpz_t b
)
333 /* convert a mpz number to a SCM number. */
334 if (mpz_fits_slong_p (b
))
336 long val
= mpz_get_si (b
);
337 if (SCM_FIXABLE (val
))
338 return SCM_I_MAKINUM (val
);
342 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
343 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
348 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
349 static SCM
scm_divide2real (SCM x
, SCM y
);
352 scm_i_make_ratio (SCM numerator
, SCM denominator
)
353 #define FUNC_NAME "make-ratio"
355 /* First make sure the arguments are proper.
357 if (SCM_I_INUMP (denominator
))
359 if (scm_is_eq (denominator
, SCM_INUM0
))
360 scm_num_overflow ("make-ratio");
361 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
366 if (!(SCM_BIGP(denominator
)))
367 SCM_WRONG_TYPE_ARG (2, denominator
);
369 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
370 SCM_WRONG_TYPE_ARG (1, numerator
);
372 /* Then flip signs so that the denominator is positive.
374 if (scm_is_true (scm_negative_p (denominator
)))
376 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
377 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
380 /* Now consider for each of the four fixnum/bignum combinations
381 whether the rational number is really an integer.
383 if (SCM_I_INUMP (numerator
))
385 long x
= SCM_I_INUM (numerator
);
386 if (scm_is_eq (numerator
, SCM_INUM0
))
388 if (SCM_I_INUMP (denominator
))
391 y
= SCM_I_INUM (denominator
);
393 return SCM_I_MAKINUM(1);
395 return SCM_I_MAKINUM (x
/ y
);
399 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
400 of that value for the denominator, as a bignum. Apart from
401 that case, abs(bignum) > abs(inum) so inum/bignum is not an
403 if (x
== SCM_MOST_NEGATIVE_FIXNUM
404 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
405 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
406 return SCM_I_MAKINUM(-1);
409 else if (SCM_BIGP (numerator
))
411 if (SCM_I_INUMP (denominator
))
413 long yy
= SCM_I_INUM (denominator
);
414 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
415 return scm_divide (numerator
, denominator
);
419 if (scm_is_eq (numerator
, denominator
))
420 return SCM_I_MAKINUM(1);
421 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
422 SCM_I_BIG_MPZ (denominator
)))
423 return scm_divide(numerator
, denominator
);
427 /* No, it's a proper fraction.
429 return scm_double_cell (scm_tc16_fraction
,
430 SCM_UNPACK (numerator
),
431 SCM_UNPACK (denominator
), 0);
435 static void scm_i_fraction_reduce (SCM z
)
437 if (!(SCM_FRACTION_REDUCED (z
)))
440 divisor
= scm_gcd (SCM_FRACTION_NUMERATOR (z
), SCM_FRACTION_DENOMINATOR (z
));
441 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
444 SCM_FRACTION_SET_NUMERATOR (z
, scm_divide (SCM_FRACTION_NUMERATOR (z
), divisor
));
445 SCM_FRACTION_SET_DENOMINATOR (z
, scm_divide (SCM_FRACTION_DENOMINATOR (z
), divisor
));
447 SCM_FRACTION_REDUCED_SET (z
);
452 scm_i_fraction2double (SCM z
)
454 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
455 SCM_FRACTION_DENOMINATOR (z
)));
458 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
460 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
462 #define FUNC_NAME s_scm_exact_p
468 if (SCM_FRACTIONP (x
))
472 SCM_WRONG_TYPE_ARG (1, x
);
477 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
479 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
481 #define FUNC_NAME s_scm_odd_p
485 long val
= SCM_I_INUM (n
);
486 return scm_from_bool ((val
& 1L) != 0);
488 else if (SCM_BIGP (n
))
490 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
491 scm_remember_upto_here_1 (n
);
492 return scm_from_bool (odd_p
);
494 else if (scm_is_true (scm_inf_p (n
)))
496 else if (SCM_REALP (n
))
498 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
504 SCM_WRONG_TYPE_ARG (1, n
);
507 SCM_WRONG_TYPE_ARG (1, n
);
512 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
514 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
516 #define FUNC_NAME s_scm_even_p
520 long val
= SCM_I_INUM (n
);
521 return scm_from_bool ((val
& 1L) == 0);
523 else if (SCM_BIGP (n
))
525 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
526 scm_remember_upto_here_1 (n
);
527 return scm_from_bool (even_p
);
529 else if (scm_is_true (scm_inf_p (n
)))
531 else if (SCM_REALP (n
))
533 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
539 SCM_WRONG_TYPE_ARG (1, n
);
542 SCM_WRONG_TYPE_ARG (1, n
);
546 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
548 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
549 "or @samp{-inf.0}, @code{#f} otherwise.")
550 #define FUNC_NAME s_scm_inf_p
553 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
554 else if (SCM_COMPLEXP (x
))
555 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
556 || xisinf (SCM_COMPLEX_IMAG (x
)));
562 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
564 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
566 #define FUNC_NAME s_scm_nan_p
569 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
570 else if (SCM_COMPLEXP (n
))
571 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
572 || xisnan (SCM_COMPLEX_IMAG (n
)));
578 /* Guile's idea of infinity. */
579 static double guile_Inf
;
581 /* Guile's idea of not a number. */
582 static double guile_NaN
;
585 guile_ieee_init (void)
587 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
589 /* Some version of gcc on some old version of Linux used to crash when
590 trying to make Inf and NaN. */
593 /* C99 INFINITY, when available.
594 FIXME: The standard allows for INFINITY to be something that overflows
595 at compile time. We ought to have a configure test to check for that
596 before trying to use it. (But in practice we believe this is not a
597 problem on any system guile is likely to target.) */
598 guile_Inf
= INFINITY
;
601 extern unsigned int DINFINITY
[2];
602 guile_Inf
= (*(X_CAST(double *, DINFINITY
)));
609 if (guile_Inf
== tmp
)
617 #if defined (HAVE_ISNAN)
620 /* C99 NAN, when available */
624 extern unsigned int DQNAN
[2];
625 guile_NaN
= (*(X_CAST(double *, DQNAN
)));
627 guile_NaN
= guile_Inf
/ guile_Inf
;
633 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
636 #define FUNC_NAME s_scm_inf
638 static int initialized
= 0;
644 return scm_from_double (guile_Inf
);
648 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
651 #define FUNC_NAME s_scm_nan
653 static int initialized
= 0;
659 return scm_from_double (guile_NaN
);
664 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
666 "Return the absolute value of @var{x}.")
671 long int xx
= SCM_I_INUM (x
);
674 else if (SCM_POSFIXABLE (-xx
))
675 return SCM_I_MAKINUM (-xx
);
677 return scm_i_long2big (-xx
);
679 else if (SCM_BIGP (x
))
681 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
683 return scm_i_clonebig (x
, 0);
687 else if (SCM_REALP (x
))
689 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
690 double xx
= SCM_REAL_VALUE (x
);
692 return scm_from_double (-xx
);
696 else if (SCM_FRACTIONP (x
))
698 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
700 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
701 SCM_FRACTION_DENOMINATOR (x
));
704 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
709 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
710 /* "Return the quotient of the numbers @var{x} and @var{y}."
713 scm_quotient (SCM x
, SCM y
)
717 long xx
= SCM_I_INUM (x
);
720 long yy
= SCM_I_INUM (y
);
722 scm_num_overflow (s_quotient
);
727 return SCM_I_MAKINUM (z
);
729 return scm_i_long2big (z
);
732 else if (SCM_BIGP (y
))
734 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
735 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
736 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
738 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
739 scm_remember_upto_here_1 (y
);
740 return SCM_I_MAKINUM (-1);
743 return SCM_I_MAKINUM (0);
746 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
748 else if (SCM_BIGP (x
))
752 long yy
= SCM_I_INUM (y
);
754 scm_num_overflow (s_quotient
);
759 SCM result
= scm_i_mkbig ();
762 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
765 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
768 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
769 scm_remember_upto_here_1 (x
);
770 return scm_i_normbig (result
);
773 else if (SCM_BIGP (y
))
775 SCM result
= scm_i_mkbig ();
776 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
779 scm_remember_upto_here_2 (x
, y
);
780 return scm_i_normbig (result
);
783 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
786 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
789 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
790 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
792 * "(remainder 13 4) @result{} 1\n"
793 * "(remainder -13 4) @result{} -1\n"
797 scm_remainder (SCM x
, SCM y
)
803 long yy
= SCM_I_INUM (y
);
805 scm_num_overflow (s_remainder
);
808 long z
= SCM_I_INUM (x
) % yy
;
809 return SCM_I_MAKINUM (z
);
812 else if (SCM_BIGP (y
))
814 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
815 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
816 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
818 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
819 scm_remember_upto_here_1 (y
);
820 return SCM_I_MAKINUM (0);
826 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
828 else if (SCM_BIGP (x
))
832 long yy
= SCM_I_INUM (y
);
834 scm_num_overflow (s_remainder
);
837 SCM result
= scm_i_mkbig ();
840 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
841 scm_remember_upto_here_1 (x
);
842 return scm_i_normbig (result
);
845 else if (SCM_BIGP (y
))
847 SCM result
= scm_i_mkbig ();
848 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
851 scm_remember_upto_here_2 (x
, y
);
852 return scm_i_normbig (result
);
855 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
858 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
862 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
863 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
865 * "(modulo 13 4) @result{} 1\n"
866 * "(modulo -13 4) @result{} 3\n"
870 scm_modulo (SCM x
, SCM y
)
874 long xx
= SCM_I_INUM (x
);
877 long yy
= SCM_I_INUM (y
);
879 scm_num_overflow (s_modulo
);
882 /* FIXME: I think this may be a bug on some arches -- results
883 of % with negative second arg are undefined... */
901 return SCM_I_MAKINUM (result
);
904 else if (SCM_BIGP (y
))
906 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
913 SCM pos_y
= scm_i_clonebig (y
, 0);
914 /* do this after the last scm_op */
915 mpz_init_set_si (z_x
, xx
);
916 result
= pos_y
; /* re-use this bignum */
917 mpz_mod (SCM_I_BIG_MPZ (result
),
919 SCM_I_BIG_MPZ (pos_y
));
920 scm_remember_upto_here_1 (pos_y
);
924 result
= scm_i_mkbig ();
925 /* do this after the last scm_op */
926 mpz_init_set_si (z_x
, xx
);
927 mpz_mod (SCM_I_BIG_MPZ (result
),
930 scm_remember_upto_here_1 (y
);
933 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
934 mpz_add (SCM_I_BIG_MPZ (result
),
936 SCM_I_BIG_MPZ (result
));
937 scm_remember_upto_here_1 (y
);
938 /* and do this before the next one */
940 return scm_i_normbig (result
);
944 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
946 else if (SCM_BIGP (x
))
950 long yy
= SCM_I_INUM (y
);
952 scm_num_overflow (s_modulo
);
955 SCM result
= scm_i_mkbig ();
956 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
958 (yy
< 0) ? - yy
: yy
);
959 scm_remember_upto_here_1 (x
);
960 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
961 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
962 SCM_I_BIG_MPZ (result
),
964 return scm_i_normbig (result
);
967 else if (SCM_BIGP (y
))
970 SCM result
= scm_i_mkbig ();
971 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
972 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
973 mpz_mod (SCM_I_BIG_MPZ (result
),
975 SCM_I_BIG_MPZ (pos_y
));
977 scm_remember_upto_here_1 (x
);
978 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
979 mpz_add (SCM_I_BIG_MPZ (result
),
981 SCM_I_BIG_MPZ (result
));
982 scm_remember_upto_here_2 (y
, pos_y
);
983 return scm_i_normbig (result
);
987 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
990 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
993 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
994 /* "Return the greatest common divisor of all arguments.\n"
995 * "If called without arguments, 0 is returned."
998 scm_gcd (SCM x
, SCM y
)
1001 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1003 if (SCM_I_INUMP (x
))
1005 if (SCM_I_INUMP (y
))
1007 long xx
= SCM_I_INUM (x
);
1008 long yy
= SCM_I_INUM (y
);
1009 long u
= xx
< 0 ? -xx
: xx
;
1010 long v
= yy
< 0 ? -yy
: yy
;
1020 /* Determine a common factor 2^k */
1021 while (!(1 & (u
| v
)))
1027 /* Now, any factor 2^n can be eliminated */
1047 return (SCM_POSFIXABLE (result
)
1048 ? SCM_I_MAKINUM (result
)
1049 : scm_i_long2big (result
));
1051 else if (SCM_BIGP (y
))
1057 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1059 else if (SCM_BIGP (x
))
1061 if (SCM_I_INUMP (y
))
1063 unsigned long result
;
1066 yy
= SCM_I_INUM (y
);
1071 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1072 scm_remember_upto_here_1 (x
);
1073 return (SCM_POSFIXABLE (result
)
1074 ? SCM_I_MAKINUM (result
)
1075 : scm_from_ulong (result
));
1077 else if (SCM_BIGP (y
))
1079 SCM result
= scm_i_mkbig ();
1080 mpz_gcd (SCM_I_BIG_MPZ (result
),
1083 scm_remember_upto_here_2 (x
, y
);
1084 return scm_i_normbig (result
);
1087 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1090 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1093 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1094 /* "Return the least common multiple of the arguments.\n"
1095 * "If called without arguments, 1 is returned."
1098 scm_lcm (SCM n1
, SCM n2
)
1100 if (SCM_UNBNDP (n2
))
1102 if (SCM_UNBNDP (n1
))
1103 return SCM_I_MAKINUM (1L);
1104 n2
= SCM_I_MAKINUM (1L);
1107 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1108 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1109 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1110 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1112 if (SCM_I_INUMP (n1
))
1114 if (SCM_I_INUMP (n2
))
1116 SCM d
= scm_gcd (n1
, n2
);
1117 if (scm_is_eq (d
, SCM_INUM0
))
1120 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1124 /* inum n1, big n2 */
1127 SCM result
= scm_i_mkbig ();
1128 long nn1
= SCM_I_INUM (n1
);
1129 if (nn1
== 0) return SCM_INUM0
;
1130 if (nn1
< 0) nn1
= - nn1
;
1131 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1132 scm_remember_upto_here_1 (n2
);
1140 if (SCM_I_INUMP (n2
))
1147 SCM result
= scm_i_mkbig ();
1148 mpz_lcm(SCM_I_BIG_MPZ (result
),
1150 SCM_I_BIG_MPZ (n2
));
1151 scm_remember_upto_here_2(n1
, n2
);
1152 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1158 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1163 + + + x (map digit:logand X Y)
1164 + - + x (map digit:logand X (lognot (+ -1 Y)))
1165 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1166 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1171 + + + (map digit:logior X Y)
1172 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1173 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1174 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1179 + + + (map digit:logxor X Y)
1180 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1181 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1182 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1187 + + (any digit:logand X Y)
1188 + - (any digit:logand X (lognot (+ -1 Y)))
1189 - + (any digit:logand (lognot (+ -1 X)) Y)
1194 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1196 "Return the bitwise AND of the integer arguments.\n\n"
1198 "(logand) @result{} -1\n"
1199 "(logand 7) @result{} 7\n"
1200 "(logand #b111 #b011 #b001) @result{} 1\n"
1202 #define FUNC_NAME s_scm_logand
1206 if (SCM_UNBNDP (n2
))
1208 if (SCM_UNBNDP (n1
))
1209 return SCM_I_MAKINUM (-1);
1210 else if (!SCM_NUMBERP (n1
))
1211 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1212 else if (SCM_NUMBERP (n1
))
1215 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1218 if (SCM_I_INUMP (n1
))
1220 nn1
= SCM_I_INUM (n1
);
1221 if (SCM_I_INUMP (n2
))
1223 long nn2
= SCM_I_INUM (n2
);
1224 return SCM_I_MAKINUM (nn1
& nn2
);
1226 else if SCM_BIGP (n2
)
1232 SCM result_z
= scm_i_mkbig ();
1234 mpz_init_set_si (nn1_z
, nn1
);
1235 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1236 scm_remember_upto_here_1 (n2
);
1238 return scm_i_normbig (result_z
);
1242 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1244 else if (SCM_BIGP (n1
))
1246 if (SCM_I_INUMP (n2
))
1249 nn1
= SCM_I_INUM (n1
);
1252 else if (SCM_BIGP (n2
))
1254 SCM result_z
= scm_i_mkbig ();
1255 mpz_and (SCM_I_BIG_MPZ (result_z
),
1257 SCM_I_BIG_MPZ (n2
));
1258 scm_remember_upto_here_2 (n1
, n2
);
1259 return scm_i_normbig (result_z
);
1262 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1265 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1270 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1272 "Return the bitwise OR of the integer arguments.\n\n"
1274 "(logior) @result{} 0\n"
1275 "(logior 7) @result{} 7\n"
1276 "(logior #b000 #b001 #b011) @result{} 3\n"
1278 #define FUNC_NAME s_scm_logior
1282 if (SCM_UNBNDP (n2
))
1284 if (SCM_UNBNDP (n1
))
1286 else if (SCM_NUMBERP (n1
))
1289 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1292 if (SCM_I_INUMP (n1
))
1294 nn1
= SCM_I_INUM (n1
);
1295 if (SCM_I_INUMP (n2
))
1297 long nn2
= SCM_I_INUM (n2
);
1298 return SCM_I_MAKINUM (nn1
| nn2
);
1300 else if (SCM_BIGP (n2
))
1306 SCM result_z
= scm_i_mkbig ();
1308 mpz_init_set_si (nn1_z
, nn1
);
1309 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1310 scm_remember_upto_here_1 (n2
);
1316 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1318 else if (SCM_BIGP (n1
))
1320 if (SCM_I_INUMP (n2
))
1323 nn1
= SCM_I_INUM (n1
);
1326 else if (SCM_BIGP (n2
))
1328 SCM result_z
= scm_i_mkbig ();
1329 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1331 SCM_I_BIG_MPZ (n2
));
1332 scm_remember_upto_here_2 (n1
, n2
);
1336 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1339 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1344 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1346 "Return the bitwise XOR of the integer arguments. A bit is\n"
1347 "set in the result if it is set in an odd number of arguments.\n"
1349 "(logxor) @result{} 0\n"
1350 "(logxor 7) @result{} 7\n"
1351 "(logxor #b000 #b001 #b011) @result{} 2\n"
1352 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1354 #define FUNC_NAME s_scm_logxor
1358 if (SCM_UNBNDP (n2
))
1360 if (SCM_UNBNDP (n1
))
1362 else if (SCM_NUMBERP (n1
))
1365 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1368 if (SCM_I_INUMP (n1
))
1370 nn1
= SCM_I_INUM (n1
);
1371 if (SCM_I_INUMP (n2
))
1373 long nn2
= SCM_I_INUM (n2
);
1374 return SCM_I_MAKINUM (nn1
^ nn2
);
1376 else if (SCM_BIGP (n2
))
1380 SCM result_z
= scm_i_mkbig ();
1382 mpz_init_set_si (nn1_z
, nn1
);
1383 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1384 scm_remember_upto_here_1 (n2
);
1386 return scm_i_normbig (result_z
);
1390 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1392 else if (SCM_BIGP (n1
))
1394 if (SCM_I_INUMP (n2
))
1397 nn1
= SCM_I_INUM (n1
);
1400 else if (SCM_BIGP (n2
))
1402 SCM result_z
= scm_i_mkbig ();
1403 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1405 SCM_I_BIG_MPZ (n2
));
1406 scm_remember_upto_here_2 (n1
, n2
);
1407 return scm_i_normbig (result_z
);
1410 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1413 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1418 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1421 "(logtest j k) @equiv{} (not (zero? (logand j k)))\n\n"
1422 "(logtest #b0100 #b1011) @result{} #f\n"
1423 "(logtest #b0100 #b0111) @result{} #t\n"
1425 #define FUNC_NAME s_scm_logtest
1429 if (SCM_I_INUMP (j
))
1431 nj
= SCM_I_INUM (j
);
1432 if (SCM_I_INUMP (k
))
1434 long nk
= SCM_I_INUM (k
);
1435 return scm_from_bool (nj
& nk
);
1437 else if (SCM_BIGP (k
))
1445 mpz_init_set_si (nj_z
, nj
);
1446 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1447 scm_remember_upto_here_1 (k
);
1448 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1454 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1456 else if (SCM_BIGP (j
))
1458 if (SCM_I_INUMP (k
))
1461 nj
= SCM_I_INUM (j
);
1464 else if (SCM_BIGP (k
))
1468 mpz_init (result_z
);
1472 scm_remember_upto_here_2 (j
, k
);
1473 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1474 mpz_clear (result_z
);
1478 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1481 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1486 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1489 "(logbit? index j) @equiv{} (logtest (integer-expt 2 index) j)\n\n"
1490 "(logbit? 0 #b1101) @result{} #t\n"
1491 "(logbit? 1 #b1101) @result{} #f\n"
1492 "(logbit? 2 #b1101) @result{} #t\n"
1493 "(logbit? 3 #b1101) @result{} #t\n"
1494 "(logbit? 4 #b1101) @result{} #f\n"
1496 #define FUNC_NAME s_scm_logbit_p
1498 unsigned long int iindex
;
1499 iindex
= scm_to_ulong (index
);
1501 if (SCM_I_INUMP (j
))
1503 /* bits above what's in an inum follow the sign bit */
1504 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1505 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1507 else if (SCM_BIGP (j
))
1509 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1510 scm_remember_upto_here_1 (j
);
1511 return scm_from_bool (val
);
1514 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1519 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1521 "Return the integer which is the ones-complement of the integer\n"
1525 "(number->string (lognot #b10000000) 2)\n"
1526 " @result{} \"-10000001\"\n"
1527 "(number->string (lognot #b0) 2)\n"
1528 " @result{} \"-1\"\n"
1530 #define FUNC_NAME s_scm_lognot
1532 if (SCM_I_INUMP (n
)) {
1533 /* No overflow here, just need to toggle all the bits making up the inum.
1534 Enhancement: No need to strip the tag and add it back, could just xor
1535 a block of 1 bits, if that worked with the various debug versions of
1537 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1539 } else if (SCM_BIGP (n
)) {
1540 SCM result
= scm_i_mkbig ();
1541 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1542 scm_remember_upto_here_1 (n
);
1546 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1551 /* returns 0 if IN is not an integer. OUT must already be
1554 coerce_to_big (SCM in
, mpz_t out
)
1557 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1558 else if (SCM_I_INUMP (in
))
1559 mpz_set_si (out
, SCM_I_INUM (in
));
1566 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1567 (SCM n
, SCM k
, SCM m
),
1568 "Return @var{n} raised to the integer exponent\n"
1569 "@var{k}, modulo @var{m}.\n"
1572 "(modulo-expt 2 3 5)\n"
1575 #define FUNC_NAME s_scm_modulo_expt
1581 /* There are two classes of error we might encounter --
1582 1) Math errors, which we'll report by calling scm_num_overflow,
1584 2) wrong-type errors, which of course we'll report by calling
1586 We don't report those errors immediately, however; instead we do
1587 some cleanup first. These variables tell us which error (if
1588 any) we should report after cleaning up.
1590 int report_overflow
= 0;
1592 int position_of_wrong_type
= 0;
1593 SCM value_of_wrong_type
= SCM_INUM0
;
1595 SCM result
= SCM_UNDEFINED
;
1601 if (scm_is_eq (m
, SCM_INUM0
))
1603 report_overflow
= 1;
1607 if (!coerce_to_big (n
, n_tmp
))
1609 value_of_wrong_type
= n
;
1610 position_of_wrong_type
= 1;
1614 if (!coerce_to_big (k
, k_tmp
))
1616 value_of_wrong_type
= k
;
1617 position_of_wrong_type
= 2;
1621 if (!coerce_to_big (m
, m_tmp
))
1623 value_of_wrong_type
= m
;
1624 position_of_wrong_type
= 3;
1628 /* if the exponent K is negative, and we simply call mpz_powm, we
1629 will get a divide-by-zero exception when an inverse 1/n mod m
1630 doesn't exist (or is not unique). Since exceptions are hard to
1631 handle, we'll attempt the inversion "by hand" -- that way, we get
1632 a simple failure code, which is easy to handle. */
1634 if (-1 == mpz_sgn (k_tmp
))
1636 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1638 report_overflow
= 1;
1641 mpz_neg (k_tmp
, k_tmp
);
1644 result
= scm_i_mkbig ();
1645 mpz_powm (SCM_I_BIG_MPZ (result
),
1650 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1651 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1658 if (report_overflow
)
1659 scm_num_overflow (FUNC_NAME
);
1661 if (position_of_wrong_type
)
1662 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1663 value_of_wrong_type
);
1665 return scm_i_normbig (result
);
1669 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1671 "Return @var{n} raised to the exact 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
);
1702 SCM_WRONG_TYPE_ARG (2, k
);
1706 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1708 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1709 n
= scm_divide (n
, SCM_UNDEFINED
);
1713 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1717 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1719 return scm_product (acc
, n
);
1721 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1722 acc
= scm_product (acc
, n
);
1723 n
= scm_product (n
, n
);
1724 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1732 n
= scm_divide (n
, SCM_UNDEFINED
);
1739 return scm_product (acc
, n
);
1741 acc
= scm_product (acc
, n
);
1742 n
= scm_product (n
, n
);
1749 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1751 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1752 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1754 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1755 "@var{cnt} is negative it's a division, rounded towards negative\n"
1756 "infinity. (Note that this is not the same rounding as\n"
1757 "@code{quotient} does.)\n"
1759 "With @var{n} viewed as an infinite precision twos complement,\n"
1760 "@code{ash} means a left shift introducing zero bits, or a right\n"
1761 "shift dropping bits.\n"
1764 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1765 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1767 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1768 "(ash -23 -2) @result{} -6\n"
1770 #define FUNC_NAME s_scm_ash
1773 bits_to_shift
= scm_to_long (cnt
);
1775 if (bits_to_shift
< 0)
1777 /* Shift right by abs(cnt) bits. This is realized as a division
1778 by div:=2^abs(cnt). However, to guarantee the floor
1779 rounding, negative values require some special treatment.
1781 SCM div
= scm_integer_expt (SCM_I_MAKINUM (2),
1782 scm_from_long (-bits_to_shift
));
1784 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1785 if (scm_is_false (scm_negative_p (n
)))
1786 return scm_quotient (n
, div
);
1788 return scm_sum (SCM_I_MAKINUM (-1L),
1789 scm_quotient (scm_sum (SCM_I_MAKINUM (1L), n
), div
));
1792 /* Shift left is done by multiplication with 2^CNT */
1793 return scm_product (n
, scm_integer_expt (SCM_I_MAKINUM (2), cnt
));
1798 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1799 (SCM n
, SCM start
, SCM end
),
1800 "Return the integer composed of the @var{start} (inclusive)\n"
1801 "through @var{end} (exclusive) bits of @var{n}. The\n"
1802 "@var{start}th bit becomes the 0-th bit in the result.\n"
1805 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1806 " @result{} \"1010\"\n"
1807 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1808 " @result{} \"10110\"\n"
1810 #define FUNC_NAME s_scm_bit_extract
1812 unsigned long int istart
, iend
, bits
;
1813 istart
= scm_to_ulong (start
);
1814 iend
= scm_to_ulong (end
);
1815 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1817 /* how many bits to keep */
1818 bits
= iend
- istart
;
1820 if (SCM_I_INUMP (n
))
1822 long int in
= SCM_I_INUM (n
);
1824 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1825 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1826 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1828 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1830 /* Since we emulate two's complement encoded numbers, this
1831 * special case requires us to produce a result that has
1832 * more bits than can be stored in a fixnum.
1834 SCM result
= scm_i_long2big (in
);
1835 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1840 /* mask down to requisite bits */
1841 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1842 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1844 else if (SCM_BIGP (n
))
1849 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1853 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1854 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1855 such bits into a ulong. */
1856 result
= scm_i_mkbig ();
1857 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1858 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1859 result
= scm_i_normbig (result
);
1861 scm_remember_upto_here_1 (n
);
1865 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1870 static const char scm_logtab
[] = {
1871 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1874 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1876 "Return the number of bits in integer @var{n}. If integer is\n"
1877 "positive, the 1-bits in its binary representation are counted.\n"
1878 "If negative, the 0-bits in its two's-complement binary\n"
1879 "representation are counted. If 0, 0 is returned.\n"
1882 "(logcount #b10101010)\n"
1889 #define FUNC_NAME s_scm_logcount
1891 if (SCM_I_INUMP (n
))
1893 unsigned long int c
= 0;
1894 long int nn
= SCM_I_INUM (n
);
1899 c
+= scm_logtab
[15 & nn
];
1902 return SCM_I_MAKINUM (c
);
1904 else if (SCM_BIGP (n
))
1906 unsigned long count
;
1907 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1908 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1910 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1911 scm_remember_upto_here_1 (n
);
1912 return SCM_I_MAKINUM (count
);
1915 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1920 static const char scm_ilentab
[] = {
1921 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1925 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1927 "Return the number of bits necessary to represent @var{n}.\n"
1930 "(integer-length #b10101010)\n"
1932 "(integer-length 0)\n"
1934 "(integer-length #b1111)\n"
1937 #define FUNC_NAME s_scm_integer_length
1939 if (SCM_I_INUMP (n
))
1941 unsigned long int c
= 0;
1943 long int nn
= SCM_I_INUM (n
);
1949 l
= scm_ilentab
[15 & nn
];
1952 return SCM_I_MAKINUM (c
- 4 + l
);
1954 else if (SCM_BIGP (n
))
1956 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1957 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1958 1 too big, so check for that and adjust. */
1959 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1960 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1961 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1962 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1964 scm_remember_upto_here_1 (n
);
1965 return SCM_I_MAKINUM (size
);
1968 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1972 /*** NUMBERS -> STRINGS ***/
1973 #define SCM_MAX_DBL_PREC 60
1974 #define SCM_MAX_DBL_RADIX 36
1976 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
1977 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
1978 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
1981 void init_dblprec(int *prec
, int radix
) {
1982 /* determine floating point precision by adding successively
1983 smaller increments to 1.0 until it is considered == 1.0 */
1984 double f
= ((double)1.0)/radix
;
1985 double fsum
= 1.0 + f
;
1990 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2002 void init_fx_radix(double *fx_list
, int radix
)
2004 /* initialize a per-radix list of tolerances. When added
2005 to a number < 1.0, we can determine if we should raund
2006 up and quit converting a number to a string. */
2010 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2011 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2014 /* use this array as a way to generate a single digit */
2015 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2018 idbl2str (double f
, char *a
, int radix
)
2020 int efmt
, dpt
, d
, i
, wp
;
2022 #ifdef DBL_MIN_10_EXP
2025 #endif /* DBL_MIN_10_EXP */
2030 radix
> SCM_MAX_DBL_RADIX
)
2032 /* revert to existing behavior */
2036 wp
= scm_dblprec
[radix
-2];
2037 fx
= fx_per_radix
[radix
-2];
2041 #ifdef HAVE_COPYSIGN
2042 double sgn
= copysign (1.0, f
);
2047 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2053 strcpy (a
, "-inf.0");
2055 strcpy (a
, "+inf.0");
2058 else if (xisnan (f
))
2060 strcpy (a
, "+nan.0");
2070 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2071 make-uniform-vector, from causing infinite loops. */
2072 /* just do the checking...if it passes, we do the conversion for our
2073 radix again below */
2080 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2088 while (f_cpy
> 10.0)
2091 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2112 if (f
+ fx
[wp
] >= radix
)
2119 /* adding 9999 makes this equivalent to abs(x) % 3 */
2120 dpt
= (exp
+ 9999) % 3;
2124 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2146 a
[ch
++] = number_chars
[d
];
2149 if (f
+ fx
[wp
] >= 1.0)
2151 a
[ch
- 1] = number_chars
[d
+1];
2163 if ((dpt
> 4) && (exp
> 6))
2165 d
= (a
[0] == '-' ? 2 : 1);
2166 for (i
= ch
++; i
> d
; i
--)
2179 if (a
[ch
- 1] == '.')
2180 a
[ch
++] = '0'; /* trailing zero */
2189 for (i
= radix
; i
<= exp
; i
*= radix
);
2190 for (i
/= radix
; i
; i
/= radix
)
2192 a
[ch
++] = number_chars
[exp
/ i
];
2200 iflo2str (SCM flt
, char *str
, int radix
)
2203 if (SCM_REALP (flt
))
2204 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2207 i
= idbl2str (SCM_COMPLEX_REAL (flt
), str
, radix
);
2208 if (SCM_COMPLEX_IMAG (flt
) != 0.0)
2210 double imag
= SCM_COMPLEX_IMAG (flt
);
2211 /* Don't output a '+' for negative numbers or for Inf and
2212 NaN. They will provide their own sign. */
2213 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2215 i
+= idbl2str (imag
, &str
[i
], radix
);
2222 /* convert a long to a string (unterminated). returns the number of
2223 characters in the result.
2225 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2227 scm_iint2str (long num
, int rad
, char *p
)
2231 unsigned long n
= (num
< 0) ? -num
: num
;
2233 for (n
/= rad
; n
> 0; n
/= rad
)
2250 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2255 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2257 "Return a string holding the external representation of the\n"
2258 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2259 "inexact, a radix of 10 will be used.")
2260 #define FUNC_NAME s_scm_number_to_string
2264 if (SCM_UNBNDP (radix
))
2267 base
= scm_to_signed_integer (radix
, 2, 36);
2269 if (SCM_I_INUMP (n
))
2271 char num_buf
[SCM_INTBUFLEN
];
2272 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2273 return scm_from_locale_stringn (num_buf
, length
);
2275 else if (SCM_BIGP (n
))
2277 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2278 scm_remember_upto_here_1 (n
);
2279 return scm_take_locale_string (str
);
2281 else if (SCM_FRACTIONP (n
))
2283 scm_i_fraction_reduce (n
);
2284 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2285 scm_from_locale_string ("/"),
2286 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2288 else if (SCM_INEXACTP (n
))
2290 char num_buf
[FLOBUFLEN
];
2291 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2294 SCM_WRONG_TYPE_ARG (1, n
);
2299 /* These print routines used to be stubbed here so that scm_repl.c
2300 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2303 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2305 char num_buf
[FLOBUFLEN
];
2306 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2311 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2314 char num_buf
[FLOBUFLEN
];
2315 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2320 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2323 scm_i_fraction_reduce (sexp
);
2324 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2325 scm_lfwrite (scm_i_string_chars (str
), scm_i_string_length (str
), port
);
2326 scm_remember_upto_here_1 (str
);
2331 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2333 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2334 scm_remember_upto_here_1 (exp
);
2335 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2339 /*** END nums->strs ***/
2342 /*** STRINGS -> NUMBERS ***/
2344 /* The following functions implement the conversion from strings to numbers.
2345 * The implementation somehow follows the grammar for numbers as it is given
2346 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2347 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2348 * points should be noted about the implementation:
2349 * * Each function keeps a local index variable 'idx' that points at the
2350 * current position within the parsed string. The global index is only
2351 * updated if the function could parse the corresponding syntactic unit
2353 * * Similarly, the functions keep track of indicators of inexactness ('#',
2354 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2355 * global exactness information is only updated after each part has been
2356 * successfully parsed.
2357 * * Sequences of digits are parsed into temporary variables holding fixnums.
2358 * Only if these fixnums would overflow, the result variables are updated
2359 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2360 * the temporary variables holding the fixnums are cleared, and the process
2361 * starts over again. If for example fixnums were able to store five decimal
2362 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2363 * and the result was computed as 12345 * 100000 + 67890. In other words,
2364 * only every five digits two bignum operations were performed.
2367 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2369 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2371 /* In non ASCII-style encodings the following macro might not work. */
2372 #define XDIGIT2UINT(d) \
2373 (isdigit ((int) (unsigned char) d) \
2375 : tolower ((int) (unsigned char) d) - 'a' + 10)
2378 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2379 unsigned int radix
, enum t_exactness
*p_exactness
)
2381 unsigned int idx
= *p_idx
;
2382 unsigned int hash_seen
= 0;
2383 scm_t_bits shift
= 1;
2385 unsigned int digit_value
;
2393 if (!isxdigit ((int) (unsigned char) c
))
2395 digit_value
= XDIGIT2UINT (c
);
2396 if (digit_value
>= radix
)
2400 result
= SCM_I_MAKINUM (digit_value
);
2404 if (isxdigit ((int) (unsigned char) c
))
2408 digit_value
= XDIGIT2UINT (c
);
2409 if (digit_value
>= radix
)
2421 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2423 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2425 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2432 shift
= shift
* radix
;
2433 add
= add
* radix
+ digit_value
;
2438 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2440 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2444 *p_exactness
= INEXACT
;
2450 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2451 * covers the parts of the rules that start at a potential point. The value
2452 * of the digits up to the point have been parsed by the caller and are given
2453 * in variable result. The content of *p_exactness indicates, whether a hash
2454 * has already been seen in the digits before the point.
2457 /* In non ASCII-style encodings the following macro might not work. */
2458 #define DIGIT2UINT(d) ((d) - '0')
2461 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2462 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2464 unsigned int idx
= *p_idx
;
2465 enum t_exactness x
= *p_exactness
;
2470 if (mem
[idx
] == '.')
2472 scm_t_bits shift
= 1;
2474 unsigned int digit_value
;
2475 SCM big_shift
= SCM_I_MAKINUM (1);
2481 if (isdigit ((int) (unsigned char) c
))
2486 digit_value
= DIGIT2UINT (c
);
2497 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2499 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2500 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2502 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2510 add
= add
* 10 + digit_value
;
2516 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2517 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2518 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2521 result
= scm_divide (result
, big_shift
);
2523 /* We've seen a decimal point, thus the value is implicitly inexact. */
2535 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2562 if (!isdigit ((int) (unsigned char) c
))
2566 exponent
= DIGIT2UINT (c
);
2570 if (isdigit ((int) (unsigned char) c
))
2573 if (exponent
<= SCM_MAXEXP
)
2574 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2580 if (exponent
> SCM_MAXEXP
)
2582 size_t exp_len
= idx
- start
;
2583 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2584 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2585 scm_out_of_range ("string->number", exp_num
);
2588 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2590 result
= scm_product (result
, e
);
2592 result
= scm_divide2real (result
, e
);
2594 /* We've seen an exponent, thus the value is implicitly inexact. */
2612 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2615 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2616 unsigned int radix
, enum t_exactness
*p_exactness
)
2618 unsigned int idx
= *p_idx
;
2624 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2630 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2632 enum t_exactness x
= EXACT
;
2634 /* Cobble up the fractional part. We might want to set the
2635 NaN's mantissa from it. */
2637 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2642 if (mem
[idx
] == '.')
2646 else if (idx
+ 1 == len
)
2648 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2651 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2652 p_idx
, p_exactness
);
2656 enum t_exactness x
= EXACT
;
2659 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2660 if (scm_is_false (uinteger
))
2665 else if (mem
[idx
] == '/')
2671 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2672 if (scm_is_false (divisor
))
2675 /* both are int/big here, I assume */
2676 result
= scm_i_make_ratio (uinteger
, divisor
);
2678 else if (radix
== 10)
2680 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2681 if (scm_is_false (result
))
2692 /* When returning an inexact zero, make sure it is represented as a
2693 floating point value so that we can change its sign.
2695 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2696 result
= scm_from_double (0.0);
2702 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2705 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2706 unsigned int radix
, enum t_exactness
*p_exactness
)
2730 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2731 if (scm_is_false (ureal
))
2733 /* input must be either +i or -i */
2738 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2744 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2751 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2752 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2761 /* either +<ureal>i or -<ureal>i */
2768 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2771 /* polar input: <real>@<real>. */
2796 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2797 if (scm_is_false (angle
))
2802 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2803 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2805 result
= scm_make_polar (ureal
, angle
);
2810 /* expecting input matching <real>[+-]<ureal>?i */
2817 int sign
= (c
== '+') ? 1 : -1;
2818 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2820 if (scm_is_false (imag
))
2821 imag
= SCM_I_MAKINUM (sign
);
2822 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2823 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2827 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2834 return scm_make_rectangular (ureal
, imag
);
2843 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2845 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2848 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2850 unsigned int idx
= 0;
2851 unsigned int radix
= NO_RADIX
;
2852 enum t_exactness forced_x
= NO_EXACTNESS
;
2853 enum t_exactness implicit_x
= EXACT
;
2856 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2857 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2859 switch (mem
[idx
+ 1])
2862 if (radix
!= NO_RADIX
)
2867 if (radix
!= NO_RADIX
)
2872 if (forced_x
!= NO_EXACTNESS
)
2877 if (forced_x
!= NO_EXACTNESS
)
2882 if (radix
!= NO_RADIX
)
2887 if (radix
!= NO_RADIX
)
2897 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2898 if (radix
== NO_RADIX
)
2899 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2901 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2903 if (scm_is_false (result
))
2909 if (SCM_INEXACTP (result
))
2910 return scm_inexact_to_exact (result
);
2914 if (SCM_INEXACTP (result
))
2917 return scm_exact_to_inexact (result
);
2920 if (implicit_x
== INEXACT
)
2922 if (SCM_INEXACTP (result
))
2925 return scm_exact_to_inexact (result
);
2933 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2934 (SCM string
, SCM radix
),
2935 "Return a number of the maximally precise representation\n"
2936 "expressed by the given @var{string}. @var{radix} must be an\n"
2937 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2938 "is a default radix that may be overridden by an explicit radix\n"
2939 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2940 "supplied, then the default radix is 10. If string is not a\n"
2941 "syntactically valid notation for a number, then\n"
2942 "@code{string->number} returns @code{#f}.")
2943 #define FUNC_NAME s_scm_string_to_number
2947 SCM_VALIDATE_STRING (1, string
);
2949 if (SCM_UNBNDP (radix
))
2952 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
2954 answer
= scm_i_mem2number (scm_i_string_chars (string
),
2955 scm_i_string_length (string
),
2957 scm_remember_upto_here_1 (string
);
2963 /*** END strs->nums ***/
2967 scm_bigequal (SCM x
, SCM y
)
2969 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
2970 scm_remember_upto_here_2 (x
, y
);
2971 return scm_from_bool (0 == result
);
2975 scm_real_equalp (SCM x
, SCM y
)
2977 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
2981 scm_complex_equalp (SCM x
, SCM y
)
2983 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
2984 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
2988 scm_i_fraction_equalp (SCM x
, SCM y
)
2990 scm_i_fraction_reduce (x
);
2991 scm_i_fraction_reduce (y
);
2992 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
2993 SCM_FRACTION_NUMERATOR (y
)))
2994 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
2995 SCM_FRACTION_DENOMINATOR (y
))))
3002 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3004 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3006 #define FUNC_NAME s_scm_number_p
3008 return scm_from_bool (SCM_NUMBERP (x
));
3012 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3014 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3015 "otherwise. Note that the sets of real, rational and integer\n"
3016 "values form subsets of the set of complex numbers, i. e. the\n"
3017 "predicate will also be fulfilled if @var{x} is a real,\n"
3018 "rational or integer number.")
3019 #define FUNC_NAME s_scm_complex_p
3021 /* all numbers are complex. */
3022 return scm_number_p (x
);
3026 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3028 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3029 "otherwise. Note that the set of integer values forms a subset of\n"
3030 "the set of real numbers, i. e. the predicate will also be\n"
3031 "fulfilled if @var{x} is an integer number.")
3032 #define FUNC_NAME s_scm_real_p
3034 /* we can't represent irrational numbers. */
3035 return scm_rational_p (x
);
3039 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3041 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3042 "otherwise. Note that the set of integer values forms a subset of\n"
3043 "the set of rational numbers, i. e. the predicate will also be\n"
3044 "fulfilled if @var{x} is an integer number.")
3045 #define FUNC_NAME s_scm_rational_p
3047 if (SCM_I_INUMP (x
))
3049 else if (SCM_IMP (x
))
3051 else if (SCM_BIGP (x
))
3053 else if (SCM_FRACTIONP (x
))
3055 else if (SCM_REALP (x
))
3056 /* due to their limited precision, all floating point numbers are
3057 rational as well. */
3064 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3066 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3068 #define FUNC_NAME s_scm_integer_p
3071 if (SCM_I_INUMP (x
))
3077 if (!SCM_INEXACTP (x
))
3079 if (SCM_COMPLEXP (x
))
3081 r
= SCM_REAL_VALUE (x
);
3082 /* +/-inf passes r==floor(r), making those #t */
3090 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3092 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3094 #define FUNC_NAME s_scm_inexact_p
3096 if (SCM_INEXACTP (x
))
3098 if (SCM_NUMBERP (x
))
3100 SCM_WRONG_TYPE_ARG (1, x
);
3105 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3106 /* "Return @code{#t} if all parameters are numerically equal." */
3108 scm_num_eq_p (SCM x
, SCM y
)
3111 if (SCM_I_INUMP (x
))
3113 long xx
= SCM_I_INUM (x
);
3114 if (SCM_I_INUMP (y
))
3116 long yy
= SCM_I_INUM (y
);
3117 return scm_from_bool (xx
== yy
);
3119 else if (SCM_BIGP (y
))
3121 else if (SCM_REALP (y
))
3122 return scm_from_bool ((double) xx
== SCM_REAL_VALUE (y
));
3123 else if (SCM_COMPLEXP (y
))
3124 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3125 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3126 else if (SCM_FRACTIONP (y
))
3129 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3131 else if (SCM_BIGP (x
))
3133 if (SCM_I_INUMP (y
))
3135 else if (SCM_BIGP (y
))
3137 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3138 scm_remember_upto_here_2 (x
, y
);
3139 return scm_from_bool (0 == cmp
);
3141 else if (SCM_REALP (y
))
3144 if (xisnan (SCM_REAL_VALUE (y
)))
3146 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3147 scm_remember_upto_here_1 (x
);
3148 return scm_from_bool (0 == cmp
);
3150 else if (SCM_COMPLEXP (y
))
3153 if (0.0 != SCM_COMPLEX_IMAG (y
))
3155 if (xisnan (SCM_COMPLEX_REAL (y
)))
3157 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3158 scm_remember_upto_here_1 (x
);
3159 return scm_from_bool (0 == cmp
);
3161 else if (SCM_FRACTIONP (y
))
3164 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3166 else if (SCM_REALP (x
))
3168 if (SCM_I_INUMP (y
))
3169 return scm_from_bool (SCM_REAL_VALUE (x
) == (double) SCM_I_INUM (y
));
3170 else if (SCM_BIGP (y
))
3173 if (xisnan (SCM_REAL_VALUE (x
)))
3175 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3176 scm_remember_upto_here_1 (y
);
3177 return scm_from_bool (0 == cmp
);
3179 else if (SCM_REALP (y
))
3180 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3181 else if (SCM_COMPLEXP (y
))
3182 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3183 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3184 else if (SCM_FRACTIONP (y
))
3186 double xx
= SCM_REAL_VALUE (x
);
3190 return scm_from_bool (xx
< 0.0);
3191 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3195 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3197 else if (SCM_COMPLEXP (x
))
3199 if (SCM_I_INUMP (y
))
3200 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3201 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3202 else if (SCM_BIGP (y
))
3205 if (0.0 != SCM_COMPLEX_IMAG (x
))
3207 if (xisnan (SCM_COMPLEX_REAL (x
)))
3209 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3210 scm_remember_upto_here_1 (y
);
3211 return scm_from_bool (0 == cmp
);
3213 else if (SCM_REALP (y
))
3214 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3215 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3216 else if (SCM_COMPLEXP (y
))
3217 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3218 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3219 else if (SCM_FRACTIONP (y
))
3222 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3224 xx
= SCM_COMPLEX_REAL (x
);
3228 return scm_from_bool (xx
< 0.0);
3229 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3233 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3235 else if (SCM_FRACTIONP (x
))
3237 if (SCM_I_INUMP (y
))
3239 else if (SCM_BIGP (y
))
3241 else if (SCM_REALP (y
))
3243 double yy
= SCM_REAL_VALUE (y
);
3247 return scm_from_bool (0.0 < yy
);
3248 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3251 else if (SCM_COMPLEXP (y
))
3254 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3256 yy
= SCM_COMPLEX_REAL (y
);
3260 return scm_from_bool (0.0 < yy
);
3261 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3264 else if (SCM_FRACTIONP (y
))
3265 return scm_i_fraction_equalp (x
, y
);
3267 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3270 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3274 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3275 done are good for inums, but for bignums an answer can almost always be
3276 had by just examining a few high bits of the operands, as done by GMP in
3277 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3278 of the float exponent to take into account. */
3280 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3281 /* "Return @code{#t} if the list of parameters is monotonically\n"
3285 scm_less_p (SCM x
, SCM y
)
3288 if (SCM_I_INUMP (x
))
3290 long xx
= SCM_I_INUM (x
);
3291 if (SCM_I_INUMP (y
))
3293 long yy
= SCM_I_INUM (y
);
3294 return scm_from_bool (xx
< yy
);
3296 else if (SCM_BIGP (y
))
3298 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3299 scm_remember_upto_here_1 (y
);
3300 return scm_from_bool (sgn
> 0);
3302 else if (SCM_REALP (y
))
3303 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3304 else if (SCM_FRACTIONP (y
))
3306 /* "x < a/b" becomes "x*b < a" */
3308 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3309 y
= SCM_FRACTION_NUMERATOR (y
);
3313 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3315 else if (SCM_BIGP (x
))
3317 if (SCM_I_INUMP (y
))
3319 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3320 scm_remember_upto_here_1 (x
);
3321 return scm_from_bool (sgn
< 0);
3323 else if (SCM_BIGP (y
))
3325 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3326 scm_remember_upto_here_2 (x
, y
);
3327 return scm_from_bool (cmp
< 0);
3329 else if (SCM_REALP (y
))
3332 if (xisnan (SCM_REAL_VALUE (y
)))
3334 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3335 scm_remember_upto_here_1 (x
);
3336 return scm_from_bool (cmp
< 0);
3338 else if (SCM_FRACTIONP (y
))
3341 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3343 else if (SCM_REALP (x
))
3345 if (SCM_I_INUMP (y
))
3346 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3347 else if (SCM_BIGP (y
))
3350 if (xisnan (SCM_REAL_VALUE (x
)))
3352 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3353 scm_remember_upto_here_1 (y
);
3354 return scm_from_bool (cmp
> 0);
3356 else if (SCM_REALP (y
))
3357 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3358 else if (SCM_FRACTIONP (y
))
3360 double xx
= SCM_REAL_VALUE (x
);
3364 return scm_from_bool (xx
< 0.0);
3365 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3369 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3371 else if (SCM_FRACTIONP (x
))
3373 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3375 /* "a/b < y" becomes "a < y*b" */
3376 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3377 x
= SCM_FRACTION_NUMERATOR (x
);
3380 else if (SCM_REALP (y
))
3382 double yy
= SCM_REAL_VALUE (y
);
3386 return scm_from_bool (0.0 < yy
);
3387 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3390 else if (SCM_FRACTIONP (y
))
3392 /* "a/b < c/d" becomes "a*d < c*b" */
3393 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3394 SCM_FRACTION_DENOMINATOR (y
));
3395 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3396 SCM_FRACTION_DENOMINATOR (x
));
3402 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3405 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3409 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3410 /* "Return @code{#t} if the list of parameters is monotonically\n"
3413 #define FUNC_NAME s_scm_gr_p
3415 scm_gr_p (SCM x
, SCM y
)
3417 if (!SCM_NUMBERP (x
))
3418 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3419 else if (!SCM_NUMBERP (y
))
3420 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3422 return scm_less_p (y
, x
);
3427 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3428 /* "Return @code{#t} if the list of parameters is monotonically\n"
3431 #define FUNC_NAME s_scm_leq_p
3433 scm_leq_p (SCM x
, SCM y
)
3435 if (!SCM_NUMBERP (x
))
3436 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3437 else if (!SCM_NUMBERP (y
))
3438 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3439 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3442 return scm_not (scm_less_p (y
, x
));
3447 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3448 /* "Return @code{#t} if the list of parameters is monotonically\n"
3451 #define FUNC_NAME s_scm_geq_p
3453 scm_geq_p (SCM x
, SCM y
)
3455 if (!SCM_NUMBERP (x
))
3456 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3457 else if (!SCM_NUMBERP (y
))
3458 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3459 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3462 return scm_not (scm_less_p (x
, y
));
3467 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3468 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3474 if (SCM_I_INUMP (z
))
3475 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3476 else if (SCM_BIGP (z
))
3478 else if (SCM_REALP (z
))
3479 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3480 else if (SCM_COMPLEXP (z
))
3481 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3482 && SCM_COMPLEX_IMAG (z
) == 0.0);
3483 else if (SCM_FRACTIONP (z
))
3486 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3490 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3491 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3495 scm_positive_p (SCM x
)
3497 if (SCM_I_INUMP (x
))
3498 return scm_from_bool (SCM_I_INUM (x
) > 0);
3499 else if (SCM_BIGP (x
))
3501 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3502 scm_remember_upto_here_1 (x
);
3503 return scm_from_bool (sgn
> 0);
3505 else if (SCM_REALP (x
))
3506 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3507 else if (SCM_FRACTIONP (x
))
3508 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3510 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3514 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3515 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3519 scm_negative_p (SCM x
)
3521 if (SCM_I_INUMP (x
))
3522 return scm_from_bool (SCM_I_INUM (x
) < 0);
3523 else if (SCM_BIGP (x
))
3525 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3526 scm_remember_upto_here_1 (x
);
3527 return scm_from_bool (sgn
< 0);
3529 else if (SCM_REALP (x
))
3530 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3531 else if (SCM_FRACTIONP (x
))
3532 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3534 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3538 /* scm_min and scm_max return an inexact when either argument is inexact, as
3539 required by r5rs. On that basis, for exact/inexact combinations the
3540 exact is converted to inexact to compare and possibly return. This is
3541 unlike scm_less_p above which takes some trouble to preserve all bits in
3542 its test, such trouble is not required for min and max. */
3544 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3545 /* "Return the maximum of all parameter values."
3548 scm_max (SCM x
, SCM y
)
3553 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3554 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3557 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3560 if (SCM_I_INUMP (x
))
3562 long xx
= SCM_I_INUM (x
);
3563 if (SCM_I_INUMP (y
))
3565 long yy
= SCM_I_INUM (y
);
3566 return (xx
< yy
) ? y
: x
;
3568 else if (SCM_BIGP (y
))
3570 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3571 scm_remember_upto_here_1 (y
);
3572 return (sgn
< 0) ? x
: y
;
3574 else if (SCM_REALP (y
))
3577 /* if y==NaN then ">" is false and we return NaN */
3578 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3580 else if (SCM_FRACTIONP (y
))
3583 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3586 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3588 else if (SCM_BIGP (x
))
3590 if (SCM_I_INUMP (y
))
3592 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3593 scm_remember_upto_here_1 (x
);
3594 return (sgn
< 0) ? y
: x
;
3596 else if (SCM_BIGP (y
))
3598 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3599 scm_remember_upto_here_2 (x
, y
);
3600 return (cmp
> 0) ? x
: y
;
3602 else if (SCM_REALP (y
))
3604 /* if y==NaN then xx>yy is false, so we return the NaN y */
3607 xx
= scm_i_big2dbl (x
);
3608 yy
= SCM_REAL_VALUE (y
);
3609 return (xx
> yy
? scm_from_double (xx
) : y
);
3611 else if (SCM_FRACTIONP (y
))
3616 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3618 else if (SCM_REALP (x
))
3620 if (SCM_I_INUMP (y
))
3622 double z
= SCM_I_INUM (y
);
3623 /* if x==NaN then "<" is false and we return NaN */
3624 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3626 else if (SCM_BIGP (y
))
3631 else if (SCM_REALP (y
))
3633 /* if x==NaN then our explicit check means we return NaN
3634 if y==NaN then ">" is false and we return NaN
3635 calling isnan is unavoidable, since it's the only way to know
3636 which of x or y causes any compares to be false */
3637 double xx
= SCM_REAL_VALUE (x
);
3638 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3640 else if (SCM_FRACTIONP (y
))
3642 double yy
= scm_i_fraction2double (y
);
3643 double xx
= SCM_REAL_VALUE (x
);
3644 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3647 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3649 else if (SCM_FRACTIONP (x
))
3651 if (SCM_I_INUMP (y
))
3655 else if (SCM_BIGP (y
))
3659 else if (SCM_REALP (y
))
3661 double xx
= scm_i_fraction2double (x
);
3662 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3664 else if (SCM_FRACTIONP (y
))
3669 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3672 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3676 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3677 /* "Return the minium of all parameter values."
3680 scm_min (SCM x
, SCM y
)
3685 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3686 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3689 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3692 if (SCM_I_INUMP (x
))
3694 long xx
= SCM_I_INUM (x
);
3695 if (SCM_I_INUMP (y
))
3697 long yy
= SCM_I_INUM (y
);
3698 return (xx
< yy
) ? x
: y
;
3700 else if (SCM_BIGP (y
))
3702 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3703 scm_remember_upto_here_1 (y
);
3704 return (sgn
< 0) ? y
: x
;
3706 else if (SCM_REALP (y
))
3709 /* if y==NaN then "<" is false and we return NaN */
3710 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3712 else if (SCM_FRACTIONP (y
))
3715 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3718 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3720 else if (SCM_BIGP (x
))
3722 if (SCM_I_INUMP (y
))
3724 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3725 scm_remember_upto_here_1 (x
);
3726 return (sgn
< 0) ? x
: y
;
3728 else if (SCM_BIGP (y
))
3730 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3731 scm_remember_upto_here_2 (x
, y
);
3732 return (cmp
> 0) ? y
: x
;
3734 else if (SCM_REALP (y
))
3736 /* if y==NaN then xx<yy is false, so we return the NaN y */
3739 xx
= scm_i_big2dbl (x
);
3740 yy
= SCM_REAL_VALUE (y
);
3741 return (xx
< yy
? scm_from_double (xx
) : y
);
3743 else if (SCM_FRACTIONP (y
))
3748 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3750 else if (SCM_REALP (x
))
3752 if (SCM_I_INUMP (y
))
3754 double z
= SCM_I_INUM (y
);
3755 /* if x==NaN then "<" is false and we return NaN */
3756 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3758 else if (SCM_BIGP (y
))
3763 else if (SCM_REALP (y
))
3765 /* if x==NaN then our explicit check means we return NaN
3766 if y==NaN then "<" is false and we return NaN
3767 calling isnan is unavoidable, since it's the only way to know
3768 which of x or y causes any compares to be false */
3769 double xx
= SCM_REAL_VALUE (x
);
3770 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3772 else if (SCM_FRACTIONP (y
))
3774 double yy
= scm_i_fraction2double (y
);
3775 double xx
= SCM_REAL_VALUE (x
);
3776 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3779 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3781 else if (SCM_FRACTIONP (x
))
3783 if (SCM_I_INUMP (y
))
3787 else if (SCM_BIGP (y
))
3791 else if (SCM_REALP (y
))
3793 double xx
= scm_i_fraction2double (x
);
3794 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3796 else if (SCM_FRACTIONP (y
))
3801 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3804 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3808 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3809 /* "Return the sum of all parameter values. Return 0 if called without\n"
3813 scm_sum (SCM x
, SCM y
)
3817 if (SCM_NUMBERP (x
)) return x
;
3818 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3819 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3822 if (SCM_I_INUMP (x
))
3824 if (SCM_I_INUMP (y
))
3826 long xx
= SCM_I_INUM (x
);
3827 long yy
= SCM_I_INUM (y
);
3828 long int z
= xx
+ yy
;
3829 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3831 else if (SCM_BIGP (y
))
3836 else if (SCM_REALP (y
))
3838 long int xx
= SCM_I_INUM (x
);
3839 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3841 else if (SCM_COMPLEXP (y
))
3843 long int xx
= SCM_I_INUM (x
);
3844 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3845 SCM_COMPLEX_IMAG (y
));
3847 else if (SCM_FRACTIONP (y
))
3848 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3849 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3850 SCM_FRACTION_DENOMINATOR (y
));
3852 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3853 } else if (SCM_BIGP (x
))
3855 if (SCM_I_INUMP (y
))
3860 inum
= SCM_I_INUM (y
);
3863 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3866 SCM result
= scm_i_mkbig ();
3867 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3868 scm_remember_upto_here_1 (x
);
3869 /* we know the result will have to be a bignum */
3872 return scm_i_normbig (result
);
3876 SCM result
= scm_i_mkbig ();
3877 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3878 scm_remember_upto_here_1 (x
);
3879 /* we know the result will have to be a bignum */
3882 return scm_i_normbig (result
);
3885 else if (SCM_BIGP (y
))
3887 SCM result
= scm_i_mkbig ();
3888 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3889 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3890 mpz_add (SCM_I_BIG_MPZ (result
),
3893 scm_remember_upto_here_2 (x
, y
);
3894 /* we know the result will have to be a bignum */
3897 return scm_i_normbig (result
);
3899 else if (SCM_REALP (y
))
3901 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3902 scm_remember_upto_here_1 (x
);
3903 return scm_from_double (result
);
3905 else if (SCM_COMPLEXP (y
))
3907 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3908 + SCM_COMPLEX_REAL (y
));
3909 scm_remember_upto_here_1 (x
);
3910 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
3912 else if (SCM_FRACTIONP (y
))
3913 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3914 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3915 SCM_FRACTION_DENOMINATOR (y
));
3917 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3919 else if (SCM_REALP (x
))
3921 if (SCM_I_INUMP (y
))
3922 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
3923 else if (SCM_BIGP (y
))
3925 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3926 scm_remember_upto_here_1 (y
);
3927 return scm_from_double (result
);
3929 else if (SCM_REALP (y
))
3930 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3931 else if (SCM_COMPLEXP (y
))
3932 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3933 SCM_COMPLEX_IMAG (y
));
3934 else if (SCM_FRACTIONP (y
))
3935 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3937 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3939 else if (SCM_COMPLEXP (x
))
3941 if (SCM_I_INUMP (y
))
3942 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
3943 SCM_COMPLEX_IMAG (x
));
3944 else if (SCM_BIGP (y
))
3946 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3947 + SCM_COMPLEX_REAL (x
));
3948 scm_remember_upto_here_1 (y
);
3949 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
3951 else if (SCM_REALP (y
))
3952 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3953 SCM_COMPLEX_IMAG (x
));
3954 else if (SCM_COMPLEXP (y
))
3955 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3956 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3957 else if (SCM_FRACTIONP (y
))
3958 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
3959 SCM_COMPLEX_IMAG (x
));
3961 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3963 else if (SCM_FRACTIONP (x
))
3965 if (SCM_I_INUMP (y
))
3966 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3967 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3968 SCM_FRACTION_DENOMINATOR (x
));
3969 else if (SCM_BIGP (y
))
3970 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
3971 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
3972 SCM_FRACTION_DENOMINATOR (x
));
3973 else if (SCM_REALP (y
))
3974 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
3975 else if (SCM_COMPLEXP (y
))
3976 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
3977 SCM_COMPLEX_IMAG (y
));
3978 else if (SCM_FRACTIONP (y
))
3979 /* a/b + c/d = (ad + bc) / bd */
3980 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
3981 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
3982 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
3984 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3987 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
3991 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
3992 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
3993 * the sum of all but the first argument are subtracted from the first
3995 #define FUNC_NAME s_difference
3997 scm_difference (SCM x
, SCM y
)
4002 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4004 if (SCM_I_INUMP (x
))
4006 long xx
= -SCM_I_INUM (x
);
4007 if (SCM_FIXABLE (xx
))
4008 return SCM_I_MAKINUM (xx
);
4010 return scm_i_long2big (xx
);
4012 else if (SCM_BIGP (x
))
4013 /* FIXME: do we really need to normalize here? */
4014 return scm_i_normbig (scm_i_clonebig (x
, 0));
4015 else if (SCM_REALP (x
))
4016 return scm_from_double (-SCM_REAL_VALUE (x
));
4017 else if (SCM_COMPLEXP (x
))
4018 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4019 -SCM_COMPLEX_IMAG (x
));
4020 else if (SCM_FRACTIONP (x
))
4021 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4022 SCM_FRACTION_DENOMINATOR (x
));
4024 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4027 if (SCM_I_INUMP (x
))
4029 if (SCM_I_INUMP (y
))
4031 long int xx
= SCM_I_INUM (x
);
4032 long int yy
= SCM_I_INUM (y
);
4033 long int z
= xx
- yy
;
4034 if (SCM_FIXABLE (z
))
4035 return SCM_I_MAKINUM (z
);
4037 return scm_i_long2big (z
);
4039 else if (SCM_BIGP (y
))
4041 /* inum-x - big-y */
4042 long xx
= SCM_I_INUM (x
);
4045 return scm_i_clonebig (y
, 0);
4048 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4049 SCM result
= scm_i_mkbig ();
4052 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4055 /* x - y == -(y + -x) */
4056 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4057 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4059 scm_remember_upto_here_1 (y
);
4061 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4062 /* we know the result will have to be a bignum */
4065 return scm_i_normbig (result
);
4068 else if (SCM_REALP (y
))
4070 long int xx
= SCM_I_INUM (x
);
4071 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4073 else if (SCM_COMPLEXP (y
))
4075 long int xx
= SCM_I_INUM (x
);
4076 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4077 - SCM_COMPLEX_IMAG (y
));
4079 else if (SCM_FRACTIONP (y
))
4080 /* a - b/c = (ac - b) / c */
4081 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4082 SCM_FRACTION_NUMERATOR (y
)),
4083 SCM_FRACTION_DENOMINATOR (y
));
4085 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4087 else if (SCM_BIGP (x
))
4089 if (SCM_I_INUMP (y
))
4091 /* big-x - inum-y */
4092 long yy
= SCM_I_INUM (y
);
4093 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4095 scm_remember_upto_here_1 (x
);
4097 return (SCM_FIXABLE (-yy
) ?
4098 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4101 SCM result
= scm_i_mkbig ();
4104 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4106 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4107 scm_remember_upto_here_1 (x
);
4109 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4110 /* we know the result will have to be a bignum */
4113 return scm_i_normbig (result
);
4116 else if (SCM_BIGP (y
))
4118 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4119 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4120 SCM result
= scm_i_mkbig ();
4121 mpz_sub (SCM_I_BIG_MPZ (result
),
4124 scm_remember_upto_here_2 (x
, y
);
4125 /* we know the result will have to be a bignum */
4126 if ((sgn_x
== 1) && (sgn_y
== -1))
4128 if ((sgn_x
== -1) && (sgn_y
== 1))
4130 return scm_i_normbig (result
);
4132 else if (SCM_REALP (y
))
4134 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4135 scm_remember_upto_here_1 (x
);
4136 return scm_from_double (result
);
4138 else if (SCM_COMPLEXP (y
))
4140 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4141 - SCM_COMPLEX_REAL (y
));
4142 scm_remember_upto_here_1 (x
);
4143 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4145 else if (SCM_FRACTIONP (y
))
4146 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4147 SCM_FRACTION_NUMERATOR (y
)),
4148 SCM_FRACTION_DENOMINATOR (y
));
4149 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4151 else if (SCM_REALP (x
))
4153 if (SCM_I_INUMP (y
))
4154 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4155 else if (SCM_BIGP (y
))
4157 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4158 scm_remember_upto_here_1 (x
);
4159 return scm_from_double (result
);
4161 else if (SCM_REALP (y
))
4162 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4163 else if (SCM_COMPLEXP (y
))
4164 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4165 -SCM_COMPLEX_IMAG (y
));
4166 else if (SCM_FRACTIONP (y
))
4167 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4169 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4171 else if (SCM_COMPLEXP (x
))
4173 if (SCM_I_INUMP (y
))
4174 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4175 SCM_COMPLEX_IMAG (x
));
4176 else if (SCM_BIGP (y
))
4178 double real_part
= (SCM_COMPLEX_REAL (x
)
4179 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4180 scm_remember_upto_here_1 (x
);
4181 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4183 else if (SCM_REALP (y
))
4184 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4185 SCM_COMPLEX_IMAG (x
));
4186 else if (SCM_COMPLEXP (y
))
4187 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4188 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4189 else if (SCM_FRACTIONP (y
))
4190 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4191 SCM_COMPLEX_IMAG (x
));
4193 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4195 else if (SCM_FRACTIONP (x
))
4197 if (SCM_I_INUMP (y
))
4198 /* a/b - c = (a - cb) / b */
4199 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4200 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4201 SCM_FRACTION_DENOMINATOR (x
));
4202 else if (SCM_BIGP (y
))
4203 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4204 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4205 SCM_FRACTION_DENOMINATOR (x
));
4206 else if (SCM_REALP (y
))
4207 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4208 else if (SCM_COMPLEXP (y
))
4209 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4210 -SCM_COMPLEX_IMAG (y
));
4211 else if (SCM_FRACTIONP (y
))
4212 /* a/b - c/d = (ad - bc) / bd */
4213 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4214 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4215 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4217 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4220 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4225 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4226 /* "Return the product of all arguments. If called without arguments,\n"
4230 scm_product (SCM x
, SCM y
)
4235 return SCM_I_MAKINUM (1L);
4236 else if (SCM_NUMBERP (x
))
4239 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4242 if (SCM_I_INUMP (x
))
4247 xx
= SCM_I_INUM (x
);
4251 case 0: return x
; break;
4252 case 1: return y
; break;
4255 if (SCM_I_INUMP (y
))
4257 long yy
= SCM_I_INUM (y
);
4259 SCM k
= SCM_I_MAKINUM (kk
);
4260 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4264 SCM result
= scm_i_long2big (xx
);
4265 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4266 return scm_i_normbig (result
);
4269 else if (SCM_BIGP (y
))
4271 SCM result
= scm_i_mkbig ();
4272 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4273 scm_remember_upto_here_1 (y
);
4276 else if (SCM_REALP (y
))
4277 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4278 else if (SCM_COMPLEXP (y
))
4279 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4280 xx
* SCM_COMPLEX_IMAG (y
));
4281 else if (SCM_FRACTIONP (y
))
4282 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4283 SCM_FRACTION_DENOMINATOR (y
));
4285 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4287 else if (SCM_BIGP (x
))
4289 if (SCM_I_INUMP (y
))
4294 else if (SCM_BIGP (y
))
4296 SCM result
= scm_i_mkbig ();
4297 mpz_mul (SCM_I_BIG_MPZ (result
),
4300 scm_remember_upto_here_2 (x
, y
);
4303 else if (SCM_REALP (y
))
4305 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4306 scm_remember_upto_here_1 (x
);
4307 return scm_from_double (result
);
4309 else if (SCM_COMPLEXP (y
))
4311 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4312 scm_remember_upto_here_1 (x
);
4313 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4314 z
* SCM_COMPLEX_IMAG (y
));
4316 else if (SCM_FRACTIONP (y
))
4317 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4318 SCM_FRACTION_DENOMINATOR (y
));
4320 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4322 else if (SCM_REALP (x
))
4324 if (SCM_I_INUMP (y
))
4325 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4326 else if (SCM_BIGP (y
))
4328 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4329 scm_remember_upto_here_1 (y
);
4330 return scm_from_double (result
);
4332 else if (SCM_REALP (y
))
4333 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4334 else if (SCM_COMPLEXP (y
))
4335 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4336 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4337 else if (SCM_FRACTIONP (y
))
4338 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4340 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4342 else if (SCM_COMPLEXP (x
))
4344 if (SCM_I_INUMP (y
))
4345 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4346 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4347 else if (SCM_BIGP (y
))
4349 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4350 scm_remember_upto_here_1 (y
);
4351 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4352 z
* SCM_COMPLEX_IMAG (x
));
4354 else if (SCM_REALP (y
))
4355 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4356 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4357 else if (SCM_COMPLEXP (y
))
4359 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4360 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4361 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4362 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4364 else if (SCM_FRACTIONP (y
))
4366 double yy
= scm_i_fraction2double (y
);
4367 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4368 yy
* SCM_COMPLEX_IMAG (x
));
4371 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4373 else if (SCM_FRACTIONP (x
))
4375 if (SCM_I_INUMP (y
))
4376 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4377 SCM_FRACTION_DENOMINATOR (x
));
4378 else if (SCM_BIGP (y
))
4379 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4380 SCM_FRACTION_DENOMINATOR (x
));
4381 else if (SCM_REALP (y
))
4382 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4383 else if (SCM_COMPLEXP (y
))
4385 double xx
= scm_i_fraction2double (x
);
4386 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4387 xx
* SCM_COMPLEX_IMAG (y
));
4389 else if (SCM_FRACTIONP (y
))
4390 /* a/b * c/d = ac / bd */
4391 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4392 SCM_FRACTION_NUMERATOR (y
)),
4393 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4394 SCM_FRACTION_DENOMINATOR (y
)));
4396 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4399 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4402 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4403 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4404 #define ALLOW_DIVIDE_BY_ZERO
4405 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4408 /* The code below for complex division is adapted from the GNU
4409 libstdc++, which adapted it from f2c's libF77, and is subject to
4412 /****************************************************************
4413 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4415 Permission to use, copy, modify, and distribute this software
4416 and its documentation for any purpose and without fee is hereby
4417 granted, provided that the above copyright notice appear in all
4418 copies and that both that the copyright notice and this
4419 permission notice and warranty disclaimer appear in supporting
4420 documentation, and that the names of AT&T Bell Laboratories or
4421 Bellcore or any of their entities not be used in advertising or
4422 publicity pertaining to distribution of the software without
4423 specific, written prior permission.
4425 AT&T and Bellcore disclaim all warranties with regard to this
4426 software, including all implied warranties of merchantability
4427 and fitness. In no event shall AT&T or Bellcore be liable for
4428 any special, indirect or consequential damages or any damages
4429 whatsoever resulting from loss of use, data or profits, whether
4430 in an action of contract, negligence or other tortious action,
4431 arising out of or in connection with the use or performance of
4433 ****************************************************************/
4435 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4436 /* Divide the first argument by the product of the remaining
4437 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4439 #define FUNC_NAME s_divide
4441 scm_i_divide (SCM x
, SCM y
, int inexact
)
4448 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4449 else if (SCM_I_INUMP (x
))
4451 long xx
= SCM_I_INUM (x
);
4452 if (xx
== 1 || xx
== -1)
4454 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4456 scm_num_overflow (s_divide
);
4461 return scm_from_double (1.0 / (double) xx
);
4462 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4465 else if (SCM_BIGP (x
))
4468 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4469 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4471 else if (SCM_REALP (x
))
4473 double xx
= SCM_REAL_VALUE (x
);
4474 #ifndef ALLOW_DIVIDE_BY_ZERO
4476 scm_num_overflow (s_divide
);
4479 return scm_from_double (1.0 / xx
);
4481 else if (SCM_COMPLEXP (x
))
4483 double r
= SCM_COMPLEX_REAL (x
);
4484 double i
= SCM_COMPLEX_IMAG (x
);
4488 double d
= i
* (1.0 + t
* t
);
4489 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4494 double d
= r
* (1.0 + t
* t
);
4495 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4498 else if (SCM_FRACTIONP (x
))
4499 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4500 SCM_FRACTION_NUMERATOR (x
));
4502 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4505 if (SCM_I_INUMP (x
))
4507 long xx
= SCM_I_INUM (x
);
4508 if (SCM_I_INUMP (y
))
4510 long yy
= SCM_I_INUM (y
);
4513 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4514 scm_num_overflow (s_divide
);
4516 return scm_from_double ((double) xx
/ (double) yy
);
4519 else if (xx
% yy
!= 0)
4522 return scm_from_double ((double) xx
/ (double) yy
);
4523 else return scm_i_make_ratio (x
, y
);
4528 if (SCM_FIXABLE (z
))
4529 return SCM_I_MAKINUM (z
);
4531 return scm_i_long2big (z
);
4534 else if (SCM_BIGP (y
))
4537 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4538 else return scm_i_make_ratio (x
, y
);
4540 else if (SCM_REALP (y
))
4542 double yy
= SCM_REAL_VALUE (y
);
4543 #ifndef ALLOW_DIVIDE_BY_ZERO
4545 scm_num_overflow (s_divide
);
4548 return scm_from_double ((double) xx
/ yy
);
4550 else if (SCM_COMPLEXP (y
))
4553 complex_div
: /* y _must_ be a complex number */
4555 double r
= SCM_COMPLEX_REAL (y
);
4556 double i
= SCM_COMPLEX_IMAG (y
);
4560 double d
= i
* (1.0 + t
* t
);
4561 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4566 double d
= r
* (1.0 + t
* t
);
4567 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4571 else if (SCM_FRACTIONP (y
))
4572 /* a / b/c = ac / b */
4573 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4574 SCM_FRACTION_NUMERATOR (y
));
4576 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4578 else if (SCM_BIGP (x
))
4580 if (SCM_I_INUMP (y
))
4582 long int yy
= SCM_I_INUM (y
);
4585 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4586 scm_num_overflow (s_divide
);
4588 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4589 scm_remember_upto_here_1 (x
);
4590 return (sgn
== 0) ? scm_nan () : scm_inf ();
4597 /* FIXME: HMM, what are the relative performance issues here?
4598 We need to test. Is it faster on average to test
4599 divisible_p, then perform whichever operation, or is it
4600 faster to perform the integer div opportunistically and
4601 switch to real if there's a remainder? For now we take the
4602 middle ground: test, then if divisible, use the faster div
4605 long abs_yy
= yy
< 0 ? -yy
: yy
;
4606 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4610 SCM result
= scm_i_mkbig ();
4611 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4612 scm_remember_upto_here_1 (x
);
4614 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4615 return scm_i_normbig (result
);
4620 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4621 else return scm_i_make_ratio (x
, y
);
4625 else if (SCM_BIGP (y
))
4627 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4630 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4631 scm_num_overflow (s_divide
);
4633 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4634 scm_remember_upto_here_1 (x
);
4635 return (sgn
== 0) ? scm_nan () : scm_inf ();
4641 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4645 SCM result
= scm_i_mkbig ();
4646 mpz_divexact (SCM_I_BIG_MPZ (result
),
4649 scm_remember_upto_here_2 (x
, y
);
4650 return scm_i_normbig (result
);
4656 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4657 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4658 scm_remember_upto_here_2 (x
, y
);
4659 return scm_from_double (dbx
/ dby
);
4661 else return scm_i_make_ratio (x
, y
);
4665 else if (SCM_REALP (y
))
4667 double yy
= SCM_REAL_VALUE (y
);
4668 #ifndef ALLOW_DIVIDE_BY_ZERO
4670 scm_num_overflow (s_divide
);
4673 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4675 else if (SCM_COMPLEXP (y
))
4677 a
= scm_i_big2dbl (x
);
4680 else if (SCM_FRACTIONP (y
))
4681 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4682 SCM_FRACTION_NUMERATOR (y
));
4684 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4686 else if (SCM_REALP (x
))
4688 double rx
= SCM_REAL_VALUE (x
);
4689 if (SCM_I_INUMP (y
))
4691 long int yy
= SCM_I_INUM (y
);
4692 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4694 scm_num_overflow (s_divide
);
4697 return scm_from_double (rx
/ (double) yy
);
4699 else if (SCM_BIGP (y
))
4701 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4702 scm_remember_upto_here_1 (y
);
4703 return scm_from_double (rx
/ dby
);
4705 else if (SCM_REALP (y
))
4707 double yy
= SCM_REAL_VALUE (y
);
4708 #ifndef ALLOW_DIVIDE_BY_ZERO
4710 scm_num_overflow (s_divide
);
4713 return scm_from_double (rx
/ yy
);
4715 else if (SCM_COMPLEXP (y
))
4720 else if (SCM_FRACTIONP (y
))
4721 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4723 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4725 else if (SCM_COMPLEXP (x
))
4727 double rx
= SCM_COMPLEX_REAL (x
);
4728 double ix
= SCM_COMPLEX_IMAG (x
);
4729 if (SCM_I_INUMP (y
))
4731 long int yy
= SCM_I_INUM (y
);
4732 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4734 scm_num_overflow (s_divide
);
4739 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4742 else if (SCM_BIGP (y
))
4744 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4745 scm_remember_upto_here_1 (y
);
4746 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4748 else if (SCM_REALP (y
))
4750 double yy
= SCM_REAL_VALUE (y
);
4751 #ifndef ALLOW_DIVIDE_BY_ZERO
4753 scm_num_overflow (s_divide
);
4756 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4758 else if (SCM_COMPLEXP (y
))
4760 double ry
= SCM_COMPLEX_REAL (y
);
4761 double iy
= SCM_COMPLEX_IMAG (y
);
4765 double d
= iy
* (1.0 + t
* t
);
4766 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4771 double d
= ry
* (1.0 + t
* t
);
4772 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4775 else if (SCM_FRACTIONP (y
))
4777 double yy
= scm_i_fraction2double (y
);
4778 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4781 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4783 else if (SCM_FRACTIONP (x
))
4785 if (SCM_I_INUMP (y
))
4787 long int yy
= SCM_I_INUM (y
);
4788 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4790 scm_num_overflow (s_divide
);
4793 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4794 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4796 else if (SCM_BIGP (y
))
4798 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4799 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4801 else if (SCM_REALP (y
))
4803 double yy
= SCM_REAL_VALUE (y
);
4804 #ifndef ALLOW_DIVIDE_BY_ZERO
4806 scm_num_overflow (s_divide
);
4809 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4811 else if (SCM_COMPLEXP (y
))
4813 a
= scm_i_fraction2double (x
);
4816 else if (SCM_FRACTIONP (y
))
4817 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4818 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4820 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4823 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4827 scm_divide (SCM x
, SCM y
)
4829 return scm_i_divide (x
, y
, 0);
4832 static SCM
scm_divide2real (SCM x
, SCM y
)
4834 return scm_i_divide (x
, y
, 1);
4840 scm_asinh (double x
)
4845 #define asinh scm_asinh
4846 return log (x
+ sqrt (x
* x
+ 1));
4849 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4850 /* "Return the inverse hyperbolic sine of @var{x}."
4855 scm_acosh (double x
)
4860 #define acosh scm_acosh
4861 return log (x
+ sqrt (x
* x
- 1));
4864 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4865 /* "Return the inverse hyperbolic cosine of @var{x}."
4870 scm_atanh (double x
)
4875 #define atanh scm_atanh
4876 return 0.5 * log ((1 + x
) / (1 - x
));
4879 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4880 /* "Return the inverse hyperbolic tangent of @var{x}."
4885 scm_c_truncate (double x
)
4896 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
4897 half-way case (ie. when x is an integer plus 0.5) going upwards.
4898 Then half-way cases are identified and adjusted down if the
4899 round-upwards didn't give the desired even integer.
4901 "plus_half == result" identifies a half-way case. If plus_half, which is
4902 x + 0.5, is an integer then x must be an integer plus 0.5.
4904 An odd "result" value is identified with result/2 != floor(result/2).
4905 This is done with plus_half, since that value is ready for use sooner in
4906 a pipelined cpu, and we're already requiring plus_half == result.
4908 Note however that we need to be careful when x is big and already an
4909 integer. In that case "x+0.5" may round to an adjacent integer, causing
4910 us to return such a value, incorrectly. For instance if the hardware is
4911 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
4912 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
4913 returned. Or if the hardware is in round-upwards mode, then other bigger
4914 values like say x == 2^128 will see x+0.5 rounding up to the next higher
4915 representable value, 2^128+2^76 (or whatever), again incorrect.
4917 These bad roundings of x+0.5 are avoided by testing at the start whether
4918 x is already an integer. If it is then clearly that's the desired result
4919 already. And if it's not then the exponent must be small enough to allow
4920 an 0.5 to be represented, and hence added without a bad rounding. */
4923 scm_c_round (double x
)
4925 double plus_half
, result
;
4930 plus_half
= x
+ 0.5;
4931 result
= floor (plus_half
);
4932 /* Adjust so that the rounding is towards even. */
4933 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4938 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4940 "Round the number @var{x} towards zero.")
4941 #define FUNC_NAME s_scm_truncate_number
4943 if (scm_is_false (scm_negative_p (x
)))
4944 return scm_floor (x
);
4946 return scm_ceiling (x
);
4950 static SCM exactly_one_half
;
4952 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4954 "Round the number @var{x} towards the nearest integer. "
4955 "When it is exactly halfway between two integers, "
4956 "round towards the even one.")
4957 #define FUNC_NAME s_scm_round_number
4959 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4961 else if (SCM_REALP (x
))
4962 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
4965 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
4966 single quotient+remainder division then examining to see which way
4967 the rounding should go. */
4968 SCM plus_half
= scm_sum (x
, exactly_one_half
);
4969 SCM result
= scm_floor (plus_half
);
4970 /* Adjust so that the rounding is towards even. */
4971 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
4972 && scm_is_true (scm_odd_p (result
)))
4973 return scm_difference (result
, SCM_I_MAKINUM (1));
4980 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
4982 "Round the number @var{x} towards minus infinity.")
4983 #define FUNC_NAME s_scm_floor
4985 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
4987 else if (SCM_REALP (x
))
4988 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
4989 else if (SCM_FRACTIONP (x
))
4991 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
4992 SCM_FRACTION_DENOMINATOR (x
));
4993 if (scm_is_false (scm_negative_p (x
)))
4995 /* For positive x, rounding towards zero is correct. */
5000 /* For negative x, we need to return q-1 unless x is an
5001 integer. But fractions are never integer, per our
5003 return scm_difference (q
, SCM_I_MAKINUM (1));
5007 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5011 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5013 "Round the number @var{x} towards infinity.")
5014 #define FUNC_NAME s_scm_ceiling
5016 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5018 else if (SCM_REALP (x
))
5019 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5020 else if (SCM_FRACTIONP (x
))
5022 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5023 SCM_FRACTION_DENOMINATOR (x
));
5024 if (scm_is_false (scm_positive_p (x
)))
5026 /* For negative x, rounding towards zero is correct. */
5031 /* For positive x, we need to return q+1 unless x is an
5032 integer. But fractions are never integer, per our
5034 return scm_sum (q
, SCM_I_MAKINUM (1));
5038 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5042 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5043 /* "Return the square root of the real number @var{x}."
5045 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5046 /* "Return the absolute value of the real number @var{x}."
5048 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5049 /* "Return the @var{x}th power of e."
5051 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5052 /* "Return the natural logarithm of the real number @var{x}."
5054 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5055 /* "Return the sine of the real number @var{x}."
5057 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5058 /* "Return the cosine of the real number @var{x}."
5060 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5061 /* "Return the tangent of the real number @var{x}."
5063 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5064 /* "Return the arc sine of the real number @var{x}."
5066 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5067 /* "Return the arc cosine of the real number @var{x}."
5069 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5070 /* "Return the arc tangent of the real number @var{x}."
5072 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5073 /* "Return the hyperbolic sine of the real number @var{x}."
5075 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5076 /* "Return the hyperbolic cosine of the real number @var{x}."
5078 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5079 /* "Return the hyperbolic tangent of the real number @var{x}."
5087 static void scm_two_doubles (SCM x
,
5089 const char *sstring
,
5093 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5095 if (SCM_I_INUMP (x
))
5096 xy
->x
= SCM_I_INUM (x
);
5097 else if (SCM_BIGP (x
))
5098 xy
->x
= scm_i_big2dbl (x
);
5099 else if (SCM_REALP (x
))
5100 xy
->x
= SCM_REAL_VALUE (x
);
5101 else if (SCM_FRACTIONP (x
))
5102 xy
->x
= scm_i_fraction2double (x
);
5104 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5106 if (SCM_I_INUMP (y
))
5107 xy
->y
= SCM_I_INUM (y
);
5108 else if (SCM_BIGP (y
))
5109 xy
->y
= scm_i_big2dbl (y
);
5110 else if (SCM_REALP (y
))
5111 xy
->y
= SCM_REAL_VALUE (y
);
5112 else if (SCM_FRACTIONP (y
))
5113 xy
->y
= scm_i_fraction2double (y
);
5115 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5119 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5121 "Return @var{x} raised to the power of @var{y}. This\n"
5122 "procedure does not accept complex arguments.")
5123 #define FUNC_NAME s_scm_sys_expt
5126 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5127 return scm_from_double (pow (xy
.x
, xy
.y
));
5132 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5134 "Return the arc tangent of the two arguments @var{x} and\n"
5135 "@var{y}. This is similar to calculating the arc tangent of\n"
5136 "@var{x} / @var{y}, except that the signs of both arguments\n"
5137 "are used to determine the quadrant of the result. This\n"
5138 "procedure does not accept complex arguments.")
5139 #define FUNC_NAME s_scm_sys_atan2
5142 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5143 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5148 scm_c_make_rectangular (double re
, double im
)
5151 return scm_from_double (re
);
5155 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5157 SCM_COMPLEX_REAL (z
) = re
;
5158 SCM_COMPLEX_IMAG (z
) = im
;
5163 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5164 (SCM real
, SCM imaginary
),
5165 "Return a complex number constructed of the given @var{real} and\n"
5166 "@var{imaginary} parts.")
5167 #define FUNC_NAME s_scm_make_rectangular
5170 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5171 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5176 scm_c_make_polar (double mag
, double ang
)
5180 sincos (ang
, &s
, &c
);
5185 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5188 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5190 "Return the complex number @var{x} * e^(i * @var{y}).")
5191 #define FUNC_NAME s_scm_make_polar
5194 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5195 return scm_c_make_polar (xy
.x
, xy
.y
);
5200 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5201 /* "Return the real part of the number @var{z}."
5204 scm_real_part (SCM z
)
5206 if (SCM_I_INUMP (z
))
5208 else if (SCM_BIGP (z
))
5210 else if (SCM_REALP (z
))
5212 else if (SCM_COMPLEXP (z
))
5213 return scm_from_double (SCM_COMPLEX_REAL (z
));
5214 else if (SCM_FRACTIONP (z
))
5217 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5221 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5222 /* "Return the imaginary part of the number @var{z}."
5225 scm_imag_part (SCM z
)
5227 if (SCM_I_INUMP (z
))
5229 else if (SCM_BIGP (z
))
5231 else if (SCM_REALP (z
))
5233 else if (SCM_COMPLEXP (z
))
5234 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5235 else if (SCM_FRACTIONP (z
))
5238 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5241 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5242 /* "Return the numerator of the number @var{z}."
5245 scm_numerator (SCM z
)
5247 if (SCM_I_INUMP (z
))
5249 else if (SCM_BIGP (z
))
5251 else if (SCM_FRACTIONP (z
))
5253 scm_i_fraction_reduce (z
);
5254 return SCM_FRACTION_NUMERATOR (z
);
5256 else if (SCM_REALP (z
))
5257 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5259 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5263 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5264 /* "Return the denominator of the number @var{z}."
5267 scm_denominator (SCM z
)
5269 if (SCM_I_INUMP (z
))
5270 return SCM_I_MAKINUM (1);
5271 else if (SCM_BIGP (z
))
5272 return SCM_I_MAKINUM (1);
5273 else if (SCM_FRACTIONP (z
))
5275 scm_i_fraction_reduce (z
);
5276 return SCM_FRACTION_DENOMINATOR (z
);
5278 else if (SCM_REALP (z
))
5279 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5281 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5284 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5285 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5286 * "@code{abs} for real arguments, but also allows complex numbers."
5289 scm_magnitude (SCM z
)
5291 if (SCM_I_INUMP (z
))
5293 long int zz
= SCM_I_INUM (z
);
5296 else if (SCM_POSFIXABLE (-zz
))
5297 return SCM_I_MAKINUM (-zz
);
5299 return scm_i_long2big (-zz
);
5301 else if (SCM_BIGP (z
))
5303 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5304 scm_remember_upto_here_1 (z
);
5306 return scm_i_clonebig (z
, 0);
5310 else if (SCM_REALP (z
))
5311 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5312 else if (SCM_COMPLEXP (z
))
5313 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5314 else if (SCM_FRACTIONP (z
))
5316 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5318 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5319 SCM_FRACTION_DENOMINATOR (z
));
5322 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5326 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5327 /* "Return the angle of the complex number @var{z}."
5332 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5333 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5334 But if atan2 follows the floating point rounding mode, then the value
5335 is not a constant. Maybe it'd be close enough though. */
5336 if (SCM_I_INUMP (z
))
5338 if (SCM_I_INUM (z
) >= 0)
5341 return scm_from_double (atan2 (0.0, -1.0));
5343 else if (SCM_BIGP (z
))
5345 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5346 scm_remember_upto_here_1 (z
);
5348 return scm_from_double (atan2 (0.0, -1.0));
5352 else if (SCM_REALP (z
))
5354 if (SCM_REAL_VALUE (z
) >= 0)
5357 return scm_from_double (atan2 (0.0, -1.0));
5359 else if (SCM_COMPLEXP (z
))
5360 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5361 else if (SCM_FRACTIONP (z
))
5363 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5365 else return scm_from_double (atan2 (0.0, -1.0));
5368 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5372 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5373 /* Convert the number @var{x} to its inexact representation.\n"
5376 scm_exact_to_inexact (SCM z
)
5378 if (SCM_I_INUMP (z
))
5379 return scm_from_double ((double) SCM_I_INUM (z
));
5380 else if (SCM_BIGP (z
))
5381 return scm_from_double (scm_i_big2dbl (z
));
5382 else if (SCM_FRACTIONP (z
))
5383 return scm_from_double (scm_i_fraction2double (z
));
5384 else if (SCM_INEXACTP (z
))
5387 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5391 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5393 "Return an exact number that is numerically closest to @var{z}.")
5394 #define FUNC_NAME s_scm_inexact_to_exact
5396 if (SCM_I_INUMP (z
))
5398 else if (SCM_BIGP (z
))
5400 else if (SCM_REALP (z
))
5402 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5403 SCM_OUT_OF_RANGE (1, z
);
5410 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5411 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5412 scm_i_mpz2num (mpq_denref (frac
)));
5414 /* When scm_i_make_ratio throws, we leak the memory allocated
5421 else if (SCM_FRACTIONP (z
))
5424 SCM_WRONG_TYPE_ARG (1, z
);
5428 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5430 "Return an exact number that is within @var{err} of @var{x}.")
5431 #define FUNC_NAME s_scm_rationalize
5433 if (SCM_I_INUMP (x
))
5435 else if (SCM_BIGP (x
))
5437 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5439 /* Use continued fractions to find closest ratio. All
5440 arithmetic is done with exact numbers.
5443 SCM ex
= scm_inexact_to_exact (x
);
5444 SCM int_part
= scm_floor (ex
);
5445 SCM tt
= SCM_I_MAKINUM (1);
5446 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5447 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5451 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5454 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5455 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5457 /* We stop after a million iterations just to be absolutely sure
5458 that we don't go into an infinite loop. The process normally
5459 converges after less than a dozen iterations.
5462 err
= scm_abs (err
);
5463 while (++i
< 1000000)
5465 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5466 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5467 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5469 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5470 err
))) /* abs(x-a/b) <= err */
5472 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5473 if (scm_is_false (scm_exact_p (x
))
5474 || scm_is_false (scm_exact_p (err
)))
5475 return scm_exact_to_inexact (res
);
5479 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5481 tt
= scm_floor (rx
); /* tt = floor (rx) */
5487 scm_num_overflow (s_scm_rationalize
);
5490 SCM_WRONG_TYPE_ARG (1, x
);
5494 /* conversion functions */
5497 scm_is_integer (SCM val
)
5499 return scm_is_true (scm_integer_p (val
));
5503 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5505 if (SCM_I_INUMP (val
))
5507 scm_t_signed_bits n
= SCM_I_INUM (val
);
5508 return n
>= min
&& n
<= max
;
5510 else if (SCM_BIGP (val
))
5512 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5514 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5516 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5518 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5519 return n
>= min
&& n
<= max
;
5529 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5530 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5533 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5534 SCM_I_BIG_MPZ (val
));
5536 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5548 return n
>= min
&& n
<= max
;
5556 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5558 if (SCM_I_INUMP (val
))
5560 scm_t_signed_bits n
= SCM_I_INUM (val
);
5561 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5563 else if (SCM_BIGP (val
))
5565 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5567 else if (max
<= ULONG_MAX
)
5569 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5571 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5572 return n
>= min
&& n
<= max
;
5582 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5585 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5586 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5589 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5590 SCM_I_BIG_MPZ (val
));
5592 return n
>= min
&& n
<= max
;
5599 #define TYPE scm_t_intmax
5600 #define TYPE_MIN min
5601 #define TYPE_MAX max
5602 #define SIZEOF_TYPE 0
5603 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5604 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5605 #include "libguile/conv-integer.i.c"
5607 #define TYPE scm_t_uintmax
5608 #define TYPE_MIN min
5609 #define TYPE_MAX max
5610 #define SIZEOF_TYPE 0
5611 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5612 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5613 #include "libguile/conv-uinteger.i.c"
5615 #define TYPE scm_t_int8
5616 #define TYPE_MIN SCM_T_INT8_MIN
5617 #define TYPE_MAX SCM_T_INT8_MAX
5618 #define SIZEOF_TYPE 1
5619 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5620 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5621 #include "libguile/conv-integer.i.c"
5623 #define TYPE scm_t_uint8
5625 #define TYPE_MAX SCM_T_UINT8_MAX
5626 #define SIZEOF_TYPE 1
5627 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5628 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5629 #include "libguile/conv-uinteger.i.c"
5631 #define TYPE scm_t_int16
5632 #define TYPE_MIN SCM_T_INT16_MIN
5633 #define TYPE_MAX SCM_T_INT16_MAX
5634 #define SIZEOF_TYPE 2
5635 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5636 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5637 #include "libguile/conv-integer.i.c"
5639 #define TYPE scm_t_uint16
5641 #define TYPE_MAX SCM_T_UINT16_MAX
5642 #define SIZEOF_TYPE 2
5643 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5644 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5645 #include "libguile/conv-uinteger.i.c"
5647 #define TYPE scm_t_int32
5648 #define TYPE_MIN SCM_T_INT32_MIN
5649 #define TYPE_MAX SCM_T_INT32_MAX
5650 #define SIZEOF_TYPE 4
5651 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5652 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5653 #include "libguile/conv-integer.i.c"
5655 #define TYPE scm_t_uint32
5657 #define TYPE_MAX SCM_T_UINT32_MAX
5658 #define SIZEOF_TYPE 4
5659 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5660 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5661 #include "libguile/conv-uinteger.i.c"
5663 #if SCM_HAVE_T_INT64
5665 #define TYPE scm_t_int64
5666 #define TYPE_MIN SCM_T_INT64_MIN
5667 #define TYPE_MAX SCM_T_INT64_MAX
5668 #define SIZEOF_TYPE 8
5669 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5670 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5671 #include "libguile/conv-integer.i.c"
5673 #define TYPE scm_t_uint64
5675 #define TYPE_MAX SCM_T_UINT64_MAX
5676 #define SIZEOF_TYPE 8
5677 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5678 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5679 #include "libguile/conv-uinteger.i.c"
5684 scm_is_real (SCM val
)
5686 return scm_is_true (scm_real_p (val
));
5690 scm_is_rational (SCM val
)
5692 return scm_is_true (scm_rational_p (val
));
5696 scm_to_double (SCM val
)
5698 if (SCM_I_INUMP (val
))
5699 return SCM_I_INUM (val
);
5700 else if (SCM_BIGP (val
))
5701 return scm_i_big2dbl (val
);
5702 else if (SCM_FRACTIONP (val
))
5703 return scm_i_fraction2double (val
);
5704 else if (SCM_REALP (val
))
5705 return SCM_REAL_VALUE (val
);
5707 scm_wrong_type_arg (NULL
, 0, val
);
5711 scm_from_double (double val
)
5713 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5714 SCM_REAL_VALUE (z
) = val
;
5718 #if SCM_ENABLE_DISCOURAGED == 1
5721 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5725 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5729 scm_out_of_range (NULL
, num
);
5732 return scm_to_double (num
);
5736 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5740 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5744 scm_out_of_range (NULL
, num
);
5747 return scm_to_double (num
);
5753 scm_is_complex (SCM val
)
5755 return scm_is_true (scm_complex_p (val
));
5759 scm_c_real_part (SCM z
)
5761 if (SCM_COMPLEXP (z
))
5762 return SCM_COMPLEX_REAL (z
);
5765 /* Use the scm_real_part to get proper error checking and
5768 return scm_to_double (scm_real_part (z
));
5773 scm_c_imag_part (SCM z
)
5775 if (SCM_COMPLEXP (z
))
5776 return SCM_COMPLEX_IMAG (z
);
5779 /* Use the scm_imag_part to get proper error checking and
5780 dispatching. The result will almost always be 0.0, but not
5783 return scm_to_double (scm_imag_part (z
));
5788 scm_c_magnitude (SCM z
)
5790 return scm_to_double (scm_magnitude (z
));
5796 return scm_to_double (scm_angle (z
));
5800 scm_is_number (SCM z
)
5802 return scm_is_true (scm_number_p (z
));
5810 mpz_init_set_si (z_negative_one
, -1);
5812 /* It may be possible to tune the performance of some algorithms by using
5813 * the following constants to avoid the creation of bignums. Please, before
5814 * using these values, remember the two rules of program optimization:
5815 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5816 scm_c_define ("most-positive-fixnum",
5817 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5818 scm_c_define ("most-negative-fixnum",
5819 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5821 scm_add_feature ("complex");
5822 scm_add_feature ("inexact");
5823 scm_flo0
= scm_from_double (0.0);
5825 /* determine floating point precision */
5826 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
5828 init_dblprec(&scm_dblprec
[i
-2],i
);
5829 init_fx_radix(fx_per_radix
[i
-2],i
);
5832 /* hard code precision for base 10 if the preprocessor tells us to... */
5833 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5836 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
5837 SCM_I_MAKINUM (2)));
5838 #include "libguile/numbers.x"