a00a5c28b416e6a3cfc86de172e7389a0a2633dc
1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 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(), csqrt(), etc */
58 #include "libguile/_scm.h"
59 #include "libguile/feature.h"
60 #include "libguile/ports.h"
61 #include "libguile/root.h"
62 #include "libguile/smob.h"
63 #include "libguile/strings.h"
65 #include "libguile/validate.h"
66 #include "libguile/numbers.h"
67 #include "libguile/deprecation.h"
69 #include "libguile/eq.h"
71 #include "libguile/discouraged.h"
73 /* values per glibc, if not already defined */
75 #define M_LOG10E 0.43429448190325182765
78 #define M_PI 3.14159265358979323846
84 Wonder if this might be faster for some of our code? A switch on
85 the numtag would jump directly to the right case, and the
86 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
88 #define SCM_I_NUMTAG_NOTNUM 0
89 #define SCM_I_NUMTAG_INUM 1
90 #define SCM_I_NUMTAG_BIG scm_tc16_big
91 #define SCM_I_NUMTAG_REAL scm_tc16_real
92 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
93 #define SCM_I_NUMTAG(x) \
94 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
95 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
96 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
97 : SCM_I_NUMTAG_NOTNUM)))
99 /* the macro above will not work as is with fractions */
102 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
104 /* FLOBUFLEN is the maximum number of characters neccessary for the
105 * printed or scm_string representation of an inexact number.
107 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
110 #if ! defined (HAVE_ISNAN)
115 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
118 #if ! defined (HAVE_ISINF)
123 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
130 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
131 an explicit check. In some future gmp (don't know what version number),
132 mpz_cmp_d is supposed to do this itself. */
134 #define xmpz_cmp_d(z, d) \
135 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
137 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
140 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
141 isinf. It does have finite and isnan though, hence the use of those.
142 fpclass would be a possibility on that system too. */
146 #if defined (HAVE_ISINF)
148 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
149 return (! (finite (x
) || isnan (x
)));
158 #if defined (HAVE_ISNAN)
166 /* For an SCM object Z which is a complex number (ie. satisfies
167 SCM_COMPLEXP), return its value as a C level "complex double". */
168 #define SCM_COMPLEX_VALUE(z) \
169 (SCM_COMPLEX_REAL (z) + _Complex_I * SCM_COMPLEX_IMAG (z))
171 /* Convert a C "complex double" to an SCM value. */
172 #if HAVE_COMPLEX_DOUBLE
174 scm_from_complex_double (complex double z
)
176 return scm_c_make_rectangular (creal (z
), cimag (z
));
178 #endif /* HAVE_COMPLEX_DOUBLE */
182 static mpz_t z_negative_one
;
189 /* Return a newly created bignum. */
190 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
191 mpz_init (SCM_I_BIG_MPZ (z
));
196 scm_i_long2big (long x
)
198 /* Return a newly created bignum initialized to X. */
199 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
200 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
205 scm_i_ulong2big (unsigned long x
)
207 /* Return a newly created bignum initialized to X. */
208 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
209 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
214 scm_i_clonebig (SCM src_big
, int same_sign_p
)
216 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
217 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
218 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
220 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
225 scm_i_bigcmp (SCM x
, SCM y
)
227 /* Return neg if x < y, pos if x > y, and 0 if x == y */
228 /* presume we already know x and y are bignums */
229 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
230 scm_remember_upto_here_2 (x
, y
);
235 scm_i_dbl2big (double d
)
237 /* results are only defined if d is an integer */
238 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
239 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
243 /* Convert a integer in double representation to a SCM number. */
246 scm_i_dbl2num (double u
)
248 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
249 powers of 2, so there's no rounding when making "double" values
250 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
251 get rounded on a 64-bit machine, hence the "+1".
253 The use of floor() to force to an integer value ensures we get a
254 "numerically closest" value without depending on how a
255 double->long cast or how mpz_set_d will round. For reference,
256 double->long probably follows the hardware rounding mode,
257 mpz_set_d truncates towards zero. */
259 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
260 representable as a double? */
262 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
263 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
264 return SCM_I_MAKINUM ((long) u
);
266 return scm_i_dbl2big (u
);
269 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
270 with R5RS exact->inexact.
272 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
273 (ie. truncate towards zero), then adjust to get the closest double by
274 examining the next lower bit and adding 1 (to the absolute value) if
277 Bignums exactly half way between representable doubles are rounded to the
278 next higher absolute value (ie. away from zero). This seems like an
279 adequate interpretation of R5RS "numerically closest", and it's easier
280 and faster than a full "nearest-even" style.
282 The bit test must be done on the absolute value of the mpz_t, which means
283 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
284 negatives as twos complement.
286 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
287 following the hardware rounding mode, but applied to the absolute value
288 of the mpz_t operand. This is not what we want so we put the high
289 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
290 mpz_get_d is supposed to always truncate towards zero.
292 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
293 is a slowdown. It'd be faster to pick out the relevant high bits with
294 mpz_getlimbn if we could be bothered coding that, and if the new
295 truncating gmp doesn't come out. */
298 scm_i_big2dbl (SCM b
)
303 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
307 /* Current GMP, eg. 4.1.3, force truncation towards zero */
309 if (bits
> DBL_MANT_DIG
)
311 size_t shift
= bits
- DBL_MANT_DIG
;
312 mpz_init2 (tmp
, DBL_MANT_DIG
);
313 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
314 result
= ldexp (mpz_get_d (tmp
), shift
);
319 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
324 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
327 if (bits
> DBL_MANT_DIG
)
329 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
330 /* test bit number "pos" in absolute value */
331 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
332 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
334 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
338 scm_remember_upto_here_1 (b
);
343 scm_i_normbig (SCM b
)
345 /* convert a big back to a fixnum if it'll fit */
346 /* presume b is a bignum */
347 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
349 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
350 if (SCM_FIXABLE (val
))
351 b
= SCM_I_MAKINUM (val
);
356 static SCM_C_INLINE_KEYWORD SCM
357 scm_i_mpz2num (mpz_t b
)
359 /* convert a mpz number to a SCM number. */
360 if (mpz_fits_slong_p (b
))
362 long val
= mpz_get_si (b
);
363 if (SCM_FIXABLE (val
))
364 return SCM_I_MAKINUM (val
);
368 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
369 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
374 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
375 static SCM
scm_divide2real (SCM x
, SCM y
);
378 scm_i_make_ratio (SCM numerator
, SCM denominator
)
379 #define FUNC_NAME "make-ratio"
381 /* First make sure the arguments are proper.
383 if (SCM_I_INUMP (denominator
))
385 if (scm_is_eq (denominator
, SCM_INUM0
))
386 scm_num_overflow ("make-ratio");
387 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
392 if (!(SCM_BIGP(denominator
)))
393 SCM_WRONG_TYPE_ARG (2, denominator
);
395 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
396 SCM_WRONG_TYPE_ARG (1, numerator
);
398 /* Then flip signs so that the denominator is positive.
400 if (scm_is_true (scm_negative_p (denominator
)))
402 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
403 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
406 /* Now consider for each of the four fixnum/bignum combinations
407 whether the rational number is really an integer.
409 if (SCM_I_INUMP (numerator
))
411 long x
= SCM_I_INUM (numerator
);
412 if (scm_is_eq (numerator
, SCM_INUM0
))
414 if (SCM_I_INUMP (denominator
))
417 y
= SCM_I_INUM (denominator
);
419 return SCM_I_MAKINUM(1);
421 return SCM_I_MAKINUM (x
/ y
);
425 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
426 of that value for the denominator, as a bignum. Apart from
427 that case, abs(bignum) > abs(inum) so inum/bignum is not an
429 if (x
== SCM_MOST_NEGATIVE_FIXNUM
430 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
431 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
432 return SCM_I_MAKINUM(-1);
435 else if (SCM_BIGP (numerator
))
437 if (SCM_I_INUMP (denominator
))
439 long yy
= SCM_I_INUM (denominator
);
440 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
441 return scm_divide (numerator
, denominator
);
445 if (scm_is_eq (numerator
, denominator
))
446 return SCM_I_MAKINUM(1);
447 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
448 SCM_I_BIG_MPZ (denominator
)))
449 return scm_divide(numerator
, denominator
);
453 /* No, it's a proper fraction.
456 SCM divisor
= scm_gcd (numerator
, denominator
);
457 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
459 numerator
= scm_divide (numerator
, divisor
);
460 denominator
= scm_divide (denominator
, divisor
);
463 return scm_double_cell (scm_tc16_fraction
,
464 SCM_UNPACK (numerator
),
465 SCM_UNPACK (denominator
), 0);
471 scm_i_fraction2double (SCM z
)
473 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
474 SCM_FRACTION_DENOMINATOR (z
)));
477 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
479 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
481 #define FUNC_NAME s_scm_exact_p
487 if (SCM_FRACTIONP (x
))
491 SCM_WRONG_TYPE_ARG (1, x
);
496 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
498 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
500 #define FUNC_NAME s_scm_odd_p
504 long val
= SCM_I_INUM (n
);
505 return scm_from_bool ((val
& 1L) != 0);
507 else if (SCM_BIGP (n
))
509 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
510 scm_remember_upto_here_1 (n
);
511 return scm_from_bool (odd_p
);
513 else if (scm_is_true (scm_inf_p (n
)))
515 else if (SCM_REALP (n
))
517 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
523 SCM_WRONG_TYPE_ARG (1, n
);
526 SCM_WRONG_TYPE_ARG (1, n
);
531 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
533 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
535 #define FUNC_NAME s_scm_even_p
539 long val
= SCM_I_INUM (n
);
540 return scm_from_bool ((val
& 1L) == 0);
542 else if (SCM_BIGP (n
))
544 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
545 scm_remember_upto_here_1 (n
);
546 return scm_from_bool (even_p
);
548 else if (scm_is_true (scm_inf_p (n
)))
550 else if (SCM_REALP (n
))
552 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
558 SCM_WRONG_TYPE_ARG (1, n
);
561 SCM_WRONG_TYPE_ARG (1, n
);
565 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
567 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
568 "or @samp{-inf.0}, @code{#f} otherwise.")
569 #define FUNC_NAME s_scm_inf_p
572 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
573 else if (SCM_COMPLEXP (x
))
574 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
575 || xisinf (SCM_COMPLEX_IMAG (x
)));
581 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
583 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
585 #define FUNC_NAME s_scm_nan_p
588 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
589 else if (SCM_COMPLEXP (n
))
590 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
591 || xisnan (SCM_COMPLEX_IMAG (n
)));
597 /* Guile's idea of infinity. */
598 static double guile_Inf
;
600 /* Guile's idea of not a number. */
601 static double guile_NaN
;
604 guile_ieee_init (void)
606 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
608 /* Some version of gcc on some old version of Linux used to crash when
609 trying to make Inf and NaN. */
612 /* C99 INFINITY, when available.
613 FIXME: The standard allows for INFINITY to be something that overflows
614 at compile time. We ought to have a configure test to check for that
615 before trying to use it. (But in practice we believe this is not a
616 problem on any system guile is likely to target.) */
617 guile_Inf
= INFINITY
;
620 extern unsigned int DINFINITY
[2];
621 guile_Inf
= (*((double *) (DINFINITY
)));
628 if (guile_Inf
== tmp
)
636 #if defined (HAVE_ISNAN)
639 /* C99 NAN, when available */
644 extern unsigned int DQNAN
[2];
645 guile_NaN
= (*((double *)(DQNAN
)));
648 guile_NaN
= guile_Inf
/ guile_Inf
;
654 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
657 #define FUNC_NAME s_scm_inf
659 static int initialized
= 0;
665 return scm_from_double (guile_Inf
);
669 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
672 #define FUNC_NAME s_scm_nan
674 static int initialized
= 0;
680 return scm_from_double (guile_NaN
);
685 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
687 "Return the absolute value of @var{x}.")
692 long int xx
= SCM_I_INUM (x
);
695 else if (SCM_POSFIXABLE (-xx
))
696 return SCM_I_MAKINUM (-xx
);
698 return scm_i_long2big (-xx
);
700 else if (SCM_BIGP (x
))
702 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
704 return scm_i_clonebig (x
, 0);
708 else if (SCM_REALP (x
))
710 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
711 double xx
= SCM_REAL_VALUE (x
);
713 return scm_from_double (-xx
);
717 else if (SCM_FRACTIONP (x
))
719 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
721 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
722 SCM_FRACTION_DENOMINATOR (x
));
725 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
730 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
731 /* "Return the quotient of the numbers @var{x} and @var{y}."
734 scm_quotient (SCM x
, SCM y
)
738 long xx
= SCM_I_INUM (x
);
741 long yy
= SCM_I_INUM (y
);
743 scm_num_overflow (s_quotient
);
748 return SCM_I_MAKINUM (z
);
750 return scm_i_long2big (z
);
753 else if (SCM_BIGP (y
))
755 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
756 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
757 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
759 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
760 scm_remember_upto_here_1 (y
);
761 return SCM_I_MAKINUM (-1);
764 return SCM_I_MAKINUM (0);
767 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
769 else if (SCM_BIGP (x
))
773 long yy
= SCM_I_INUM (y
);
775 scm_num_overflow (s_quotient
);
780 SCM result
= scm_i_mkbig ();
783 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
786 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
789 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
790 scm_remember_upto_here_1 (x
);
791 return scm_i_normbig (result
);
794 else if (SCM_BIGP (y
))
796 SCM result
= scm_i_mkbig ();
797 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
800 scm_remember_upto_here_2 (x
, y
);
801 return scm_i_normbig (result
);
804 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
807 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
810 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
811 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
813 * "(remainder 13 4) @result{} 1\n"
814 * "(remainder -13 4) @result{} -1\n"
818 scm_remainder (SCM x
, SCM y
)
824 long yy
= SCM_I_INUM (y
);
826 scm_num_overflow (s_remainder
);
829 long z
= SCM_I_INUM (x
) % yy
;
830 return SCM_I_MAKINUM (z
);
833 else if (SCM_BIGP (y
))
835 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
836 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
837 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
839 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
840 scm_remember_upto_here_1 (y
);
841 return SCM_I_MAKINUM (0);
847 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
849 else if (SCM_BIGP (x
))
853 long yy
= SCM_I_INUM (y
);
855 scm_num_overflow (s_remainder
);
858 SCM result
= scm_i_mkbig ();
861 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
862 scm_remember_upto_here_1 (x
);
863 return scm_i_normbig (result
);
866 else if (SCM_BIGP (y
))
868 SCM result
= scm_i_mkbig ();
869 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
872 scm_remember_upto_here_2 (x
, y
);
873 return scm_i_normbig (result
);
876 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
879 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
883 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
884 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
886 * "(modulo 13 4) @result{} 1\n"
887 * "(modulo -13 4) @result{} 3\n"
891 scm_modulo (SCM x
, SCM y
)
895 long xx
= SCM_I_INUM (x
);
898 long yy
= SCM_I_INUM (y
);
900 scm_num_overflow (s_modulo
);
903 /* C99 specifies that "%" is the remainder corresponding to a
904 quotient rounded towards zero, and that's also traditional
905 for machine division, so z here should be well defined. */
923 return SCM_I_MAKINUM (result
);
926 else if (SCM_BIGP (y
))
928 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
935 SCM pos_y
= scm_i_clonebig (y
, 0);
936 /* do this after the last scm_op */
937 mpz_init_set_si (z_x
, xx
);
938 result
= pos_y
; /* re-use this bignum */
939 mpz_mod (SCM_I_BIG_MPZ (result
),
941 SCM_I_BIG_MPZ (pos_y
));
942 scm_remember_upto_here_1 (pos_y
);
946 result
= scm_i_mkbig ();
947 /* do this after the last scm_op */
948 mpz_init_set_si (z_x
, xx
);
949 mpz_mod (SCM_I_BIG_MPZ (result
),
952 scm_remember_upto_here_1 (y
);
955 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
956 mpz_add (SCM_I_BIG_MPZ (result
),
958 SCM_I_BIG_MPZ (result
));
959 scm_remember_upto_here_1 (y
);
960 /* and do this before the next one */
962 return scm_i_normbig (result
);
966 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
968 else if (SCM_BIGP (x
))
972 long yy
= SCM_I_INUM (y
);
974 scm_num_overflow (s_modulo
);
977 SCM result
= scm_i_mkbig ();
978 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
980 (yy
< 0) ? - yy
: yy
);
981 scm_remember_upto_here_1 (x
);
982 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
983 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
984 SCM_I_BIG_MPZ (result
),
986 return scm_i_normbig (result
);
989 else if (SCM_BIGP (y
))
992 SCM result
= scm_i_mkbig ();
993 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
994 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
995 mpz_mod (SCM_I_BIG_MPZ (result
),
997 SCM_I_BIG_MPZ (pos_y
));
999 scm_remember_upto_here_1 (x
);
1000 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1001 mpz_add (SCM_I_BIG_MPZ (result
),
1003 SCM_I_BIG_MPZ (result
));
1004 scm_remember_upto_here_2 (y
, pos_y
);
1005 return scm_i_normbig (result
);
1009 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1012 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1015 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
1016 /* "Return the greatest common divisor of all arguments.\n"
1017 * "If called without arguments, 0 is returned."
1020 scm_gcd (SCM x
, SCM y
)
1023 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1025 if (SCM_I_INUMP (x
))
1027 if (SCM_I_INUMP (y
))
1029 long xx
= SCM_I_INUM (x
);
1030 long yy
= SCM_I_INUM (y
);
1031 long u
= xx
< 0 ? -xx
: xx
;
1032 long v
= yy
< 0 ? -yy
: yy
;
1042 /* Determine a common factor 2^k */
1043 while (!(1 & (u
| v
)))
1049 /* Now, any factor 2^n can be eliminated */
1069 return (SCM_POSFIXABLE (result
)
1070 ? SCM_I_MAKINUM (result
)
1071 : scm_i_long2big (result
));
1073 else if (SCM_BIGP (y
))
1079 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1081 else if (SCM_BIGP (x
))
1083 if (SCM_I_INUMP (y
))
1085 unsigned long result
;
1088 yy
= SCM_I_INUM (y
);
1093 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1094 scm_remember_upto_here_1 (x
);
1095 return (SCM_POSFIXABLE (result
)
1096 ? SCM_I_MAKINUM (result
)
1097 : scm_from_ulong (result
));
1099 else if (SCM_BIGP (y
))
1101 SCM result
= scm_i_mkbig ();
1102 mpz_gcd (SCM_I_BIG_MPZ (result
),
1105 scm_remember_upto_here_2 (x
, y
);
1106 return scm_i_normbig (result
);
1109 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1112 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1115 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1116 /* "Return the least common multiple of the arguments.\n"
1117 * "If called without arguments, 1 is returned."
1120 scm_lcm (SCM n1
, SCM n2
)
1122 if (SCM_UNBNDP (n2
))
1124 if (SCM_UNBNDP (n1
))
1125 return SCM_I_MAKINUM (1L);
1126 n2
= SCM_I_MAKINUM (1L);
1129 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1130 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1131 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1132 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1134 if (SCM_I_INUMP (n1
))
1136 if (SCM_I_INUMP (n2
))
1138 SCM d
= scm_gcd (n1
, n2
);
1139 if (scm_is_eq (d
, SCM_INUM0
))
1142 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1146 /* inum n1, big n2 */
1149 SCM result
= scm_i_mkbig ();
1150 long nn1
= SCM_I_INUM (n1
);
1151 if (nn1
== 0) return SCM_INUM0
;
1152 if (nn1
< 0) nn1
= - nn1
;
1153 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1154 scm_remember_upto_here_1 (n2
);
1162 if (SCM_I_INUMP (n2
))
1169 SCM result
= scm_i_mkbig ();
1170 mpz_lcm(SCM_I_BIG_MPZ (result
),
1172 SCM_I_BIG_MPZ (n2
));
1173 scm_remember_upto_here_2(n1
, n2
);
1174 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1180 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1185 + + + x (map digit:logand X Y)
1186 + - + x (map digit:logand X (lognot (+ -1 Y)))
1187 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1188 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1193 + + + (map digit:logior X Y)
1194 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1195 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1196 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1201 + + + (map digit:logxor X Y)
1202 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1203 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1204 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1209 + + (any digit:logand X Y)
1210 + - (any digit:logand X (lognot (+ -1 Y)))
1211 - + (any digit:logand (lognot (+ -1 X)) Y)
1216 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1218 "Return the bitwise AND of the integer arguments.\n\n"
1220 "(logand) @result{} -1\n"
1221 "(logand 7) @result{} 7\n"
1222 "(logand #b111 #b011 #b001) @result{} 1\n"
1224 #define FUNC_NAME s_scm_logand
1228 if (SCM_UNBNDP (n2
))
1230 if (SCM_UNBNDP (n1
))
1231 return SCM_I_MAKINUM (-1);
1232 else if (!SCM_NUMBERP (n1
))
1233 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1234 else if (SCM_NUMBERP (n1
))
1237 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1240 if (SCM_I_INUMP (n1
))
1242 nn1
= SCM_I_INUM (n1
);
1243 if (SCM_I_INUMP (n2
))
1245 long nn2
= SCM_I_INUM (n2
);
1246 return SCM_I_MAKINUM (nn1
& nn2
);
1248 else if SCM_BIGP (n2
)
1254 SCM result_z
= scm_i_mkbig ();
1256 mpz_init_set_si (nn1_z
, nn1
);
1257 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1258 scm_remember_upto_here_1 (n2
);
1260 return scm_i_normbig (result_z
);
1264 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1266 else if (SCM_BIGP (n1
))
1268 if (SCM_I_INUMP (n2
))
1271 nn1
= SCM_I_INUM (n1
);
1274 else if (SCM_BIGP (n2
))
1276 SCM result_z
= scm_i_mkbig ();
1277 mpz_and (SCM_I_BIG_MPZ (result_z
),
1279 SCM_I_BIG_MPZ (n2
));
1280 scm_remember_upto_here_2 (n1
, n2
);
1281 return scm_i_normbig (result_z
);
1284 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1287 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1292 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1294 "Return the bitwise OR of the integer arguments.\n\n"
1296 "(logior) @result{} 0\n"
1297 "(logior 7) @result{} 7\n"
1298 "(logior #b000 #b001 #b011) @result{} 3\n"
1300 #define FUNC_NAME s_scm_logior
1304 if (SCM_UNBNDP (n2
))
1306 if (SCM_UNBNDP (n1
))
1308 else if (SCM_NUMBERP (n1
))
1311 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1314 if (SCM_I_INUMP (n1
))
1316 nn1
= SCM_I_INUM (n1
);
1317 if (SCM_I_INUMP (n2
))
1319 long nn2
= SCM_I_INUM (n2
);
1320 return SCM_I_MAKINUM (nn1
| nn2
);
1322 else if (SCM_BIGP (n2
))
1328 SCM result_z
= scm_i_mkbig ();
1330 mpz_init_set_si (nn1_z
, nn1
);
1331 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1332 scm_remember_upto_here_1 (n2
);
1334 return scm_i_normbig (result_z
);
1338 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1340 else if (SCM_BIGP (n1
))
1342 if (SCM_I_INUMP (n2
))
1345 nn1
= SCM_I_INUM (n1
);
1348 else if (SCM_BIGP (n2
))
1350 SCM result_z
= scm_i_mkbig ();
1351 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1353 SCM_I_BIG_MPZ (n2
));
1354 scm_remember_upto_here_2 (n1
, n2
);
1355 return scm_i_normbig (result_z
);
1358 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1361 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1366 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1368 "Return the bitwise XOR of the integer arguments. A bit is\n"
1369 "set in the result if it is set in an odd number of arguments.\n"
1371 "(logxor) @result{} 0\n"
1372 "(logxor 7) @result{} 7\n"
1373 "(logxor #b000 #b001 #b011) @result{} 2\n"
1374 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1376 #define FUNC_NAME s_scm_logxor
1380 if (SCM_UNBNDP (n2
))
1382 if (SCM_UNBNDP (n1
))
1384 else if (SCM_NUMBERP (n1
))
1387 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1390 if (SCM_I_INUMP (n1
))
1392 nn1
= SCM_I_INUM (n1
);
1393 if (SCM_I_INUMP (n2
))
1395 long nn2
= SCM_I_INUM (n2
);
1396 return SCM_I_MAKINUM (nn1
^ nn2
);
1398 else if (SCM_BIGP (n2
))
1402 SCM result_z
= scm_i_mkbig ();
1404 mpz_init_set_si (nn1_z
, nn1
);
1405 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1406 scm_remember_upto_here_1 (n2
);
1408 return scm_i_normbig (result_z
);
1412 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1414 else if (SCM_BIGP (n1
))
1416 if (SCM_I_INUMP (n2
))
1419 nn1
= SCM_I_INUM (n1
);
1422 else if (SCM_BIGP (n2
))
1424 SCM result_z
= scm_i_mkbig ();
1425 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1427 SCM_I_BIG_MPZ (n2
));
1428 scm_remember_upto_here_2 (n1
, n2
);
1429 return scm_i_normbig (result_z
);
1432 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1435 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1440 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1442 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1443 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1444 "without actually calculating the @code{logand}, just testing\n"
1448 "(logtest #b0100 #b1011) @result{} #f\n"
1449 "(logtest #b0100 #b0111) @result{} #t\n"
1451 #define FUNC_NAME s_scm_logtest
1455 if (SCM_I_INUMP (j
))
1457 nj
= SCM_I_INUM (j
);
1458 if (SCM_I_INUMP (k
))
1460 long nk
= SCM_I_INUM (k
);
1461 return scm_from_bool (nj
& nk
);
1463 else if (SCM_BIGP (k
))
1471 mpz_init_set_si (nj_z
, nj
);
1472 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1473 scm_remember_upto_here_1 (k
);
1474 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1480 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1482 else if (SCM_BIGP (j
))
1484 if (SCM_I_INUMP (k
))
1487 nj
= SCM_I_INUM (j
);
1490 else if (SCM_BIGP (k
))
1494 mpz_init (result_z
);
1498 scm_remember_upto_here_2 (j
, k
);
1499 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1500 mpz_clear (result_z
);
1504 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1507 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1512 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1514 "Test whether bit number @var{index} in @var{j} is set.\n"
1515 "@var{index} starts from 0 for the least significant bit.\n"
1518 "(logbit? 0 #b1101) @result{} #t\n"
1519 "(logbit? 1 #b1101) @result{} #f\n"
1520 "(logbit? 2 #b1101) @result{} #t\n"
1521 "(logbit? 3 #b1101) @result{} #t\n"
1522 "(logbit? 4 #b1101) @result{} #f\n"
1524 #define FUNC_NAME s_scm_logbit_p
1526 unsigned long int iindex
;
1527 iindex
= scm_to_ulong (index
);
1529 if (SCM_I_INUMP (j
))
1531 /* bits above what's in an inum follow the sign bit */
1532 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1533 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1535 else if (SCM_BIGP (j
))
1537 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1538 scm_remember_upto_here_1 (j
);
1539 return scm_from_bool (val
);
1542 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1547 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1549 "Return the integer which is the ones-complement of the integer\n"
1553 "(number->string (lognot #b10000000) 2)\n"
1554 " @result{} \"-10000001\"\n"
1555 "(number->string (lognot #b0) 2)\n"
1556 " @result{} \"-1\"\n"
1558 #define FUNC_NAME s_scm_lognot
1560 if (SCM_I_INUMP (n
)) {
1561 /* No overflow here, just need to toggle all the bits making up the inum.
1562 Enhancement: No need to strip the tag and add it back, could just xor
1563 a block of 1 bits, if that worked with the various debug versions of
1565 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1567 } else if (SCM_BIGP (n
)) {
1568 SCM result
= scm_i_mkbig ();
1569 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1570 scm_remember_upto_here_1 (n
);
1574 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1579 /* returns 0 if IN is not an integer. OUT must already be
1582 coerce_to_big (SCM in
, mpz_t out
)
1585 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1586 else if (SCM_I_INUMP (in
))
1587 mpz_set_si (out
, SCM_I_INUM (in
));
1594 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1595 (SCM n
, SCM k
, SCM m
),
1596 "Return @var{n} raised to the integer exponent\n"
1597 "@var{k}, modulo @var{m}.\n"
1600 "(modulo-expt 2 3 5)\n"
1603 #define FUNC_NAME s_scm_modulo_expt
1609 /* There are two classes of error we might encounter --
1610 1) Math errors, which we'll report by calling scm_num_overflow,
1612 2) wrong-type errors, which of course we'll report by calling
1614 We don't report those errors immediately, however; instead we do
1615 some cleanup first. These variables tell us which error (if
1616 any) we should report after cleaning up.
1618 int report_overflow
= 0;
1620 int position_of_wrong_type
= 0;
1621 SCM value_of_wrong_type
= SCM_INUM0
;
1623 SCM result
= SCM_UNDEFINED
;
1629 if (scm_is_eq (m
, SCM_INUM0
))
1631 report_overflow
= 1;
1635 if (!coerce_to_big (n
, n_tmp
))
1637 value_of_wrong_type
= n
;
1638 position_of_wrong_type
= 1;
1642 if (!coerce_to_big (k
, k_tmp
))
1644 value_of_wrong_type
= k
;
1645 position_of_wrong_type
= 2;
1649 if (!coerce_to_big (m
, m_tmp
))
1651 value_of_wrong_type
= m
;
1652 position_of_wrong_type
= 3;
1656 /* if the exponent K is negative, and we simply call mpz_powm, we
1657 will get a divide-by-zero exception when an inverse 1/n mod m
1658 doesn't exist (or is not unique). Since exceptions are hard to
1659 handle, we'll attempt the inversion "by hand" -- that way, we get
1660 a simple failure code, which is easy to handle. */
1662 if (-1 == mpz_sgn (k_tmp
))
1664 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1666 report_overflow
= 1;
1669 mpz_neg (k_tmp
, k_tmp
);
1672 result
= scm_i_mkbig ();
1673 mpz_powm (SCM_I_BIG_MPZ (result
),
1678 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1679 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1686 if (report_overflow
)
1687 scm_num_overflow (FUNC_NAME
);
1689 if (position_of_wrong_type
)
1690 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1691 value_of_wrong_type
);
1693 return scm_i_normbig (result
);
1697 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1699 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1700 "exact integer, @var{n} can be any number.\n"
1702 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1703 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1704 "includes @math{0^0} is 1.\n"
1707 "(integer-expt 2 5) @result{} 32\n"
1708 "(integer-expt -3 3) @result{} -27\n"
1709 "(integer-expt 5 -3) @result{} 1/125\n"
1710 "(integer-expt 0 0) @result{} 1\n"
1712 #define FUNC_NAME s_scm_integer_expt
1715 SCM z_i2
= SCM_BOOL_F
;
1717 SCM acc
= SCM_I_MAKINUM (1L);
1719 /* 0^0 == 1 according to R5RS */
1720 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1721 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1722 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1723 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1725 if (SCM_I_INUMP (k
))
1726 i2
= SCM_I_INUM (k
);
1727 else if (SCM_BIGP (k
))
1729 z_i2
= scm_i_clonebig (k
, 1);
1730 scm_remember_upto_here_1 (k
);
1734 SCM_WRONG_TYPE_ARG (2, k
);
1738 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1740 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1741 n
= scm_divide (n
, SCM_UNDEFINED
);
1745 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1749 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1751 return scm_product (acc
, n
);
1753 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1754 acc
= scm_product (acc
, n
);
1755 n
= scm_product (n
, n
);
1756 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1764 n
= scm_divide (n
, SCM_UNDEFINED
);
1771 return scm_product (acc
, n
);
1773 acc
= scm_product (acc
, n
);
1774 n
= scm_product (n
, n
);
1781 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1783 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1784 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1786 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1787 "@var{cnt} is negative it's a division, rounded towards negative\n"
1788 "infinity. (Note that this is not the same rounding as\n"
1789 "@code{quotient} does.)\n"
1791 "With @var{n} viewed as an infinite precision twos complement,\n"
1792 "@code{ash} means a left shift introducing zero bits, or a right\n"
1793 "shift dropping bits.\n"
1796 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1797 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1799 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1800 "(ash -23 -2) @result{} -6\n"
1802 #define FUNC_NAME s_scm_ash
1805 bits_to_shift
= scm_to_long (cnt
);
1807 if (SCM_I_INUMP (n
))
1809 long nn
= SCM_I_INUM (n
);
1811 if (bits_to_shift
> 0)
1813 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1814 overflow a non-zero fixnum. For smaller shifts we check the
1815 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1816 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1817 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1823 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1825 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1828 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1832 SCM result
= scm_i_long2big (nn
);
1833 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1840 bits_to_shift
= -bits_to_shift
;
1841 if (bits_to_shift
>= SCM_LONG_BIT
)
1842 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1844 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1848 else if (SCM_BIGP (n
))
1852 if (bits_to_shift
== 0)
1855 result
= scm_i_mkbig ();
1856 if (bits_to_shift
>= 0)
1858 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1864 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1865 we have to allocate a bignum even if the result is going to be a
1867 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1869 return scm_i_normbig (result
);
1875 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1881 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1882 (SCM n
, SCM start
, SCM end
),
1883 "Return the integer composed of the @var{start} (inclusive)\n"
1884 "through @var{end} (exclusive) bits of @var{n}. The\n"
1885 "@var{start}th bit becomes the 0-th bit in the result.\n"
1888 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1889 " @result{} \"1010\"\n"
1890 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1891 " @result{} \"10110\"\n"
1893 #define FUNC_NAME s_scm_bit_extract
1895 unsigned long int istart
, iend
, bits
;
1896 istart
= scm_to_ulong (start
);
1897 iend
= scm_to_ulong (end
);
1898 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1900 /* how many bits to keep */
1901 bits
= iend
- istart
;
1903 if (SCM_I_INUMP (n
))
1905 long int in
= SCM_I_INUM (n
);
1907 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1908 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1909 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1911 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1913 /* Since we emulate two's complement encoded numbers, this
1914 * special case requires us to produce a result that has
1915 * more bits than can be stored in a fixnum.
1917 SCM result
= scm_i_long2big (in
);
1918 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1923 /* mask down to requisite bits */
1924 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1925 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1927 else if (SCM_BIGP (n
))
1932 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1936 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1937 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1938 such bits into a ulong. */
1939 result
= scm_i_mkbig ();
1940 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1941 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1942 result
= scm_i_normbig (result
);
1944 scm_remember_upto_here_1 (n
);
1948 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1953 static const char scm_logtab
[] = {
1954 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1957 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1959 "Return the number of bits in integer @var{n}. If integer is\n"
1960 "positive, the 1-bits in its binary representation are counted.\n"
1961 "If negative, the 0-bits in its two's-complement binary\n"
1962 "representation are counted. If 0, 0 is returned.\n"
1965 "(logcount #b10101010)\n"
1972 #define FUNC_NAME s_scm_logcount
1974 if (SCM_I_INUMP (n
))
1976 unsigned long int c
= 0;
1977 long int nn
= SCM_I_INUM (n
);
1982 c
+= scm_logtab
[15 & nn
];
1985 return SCM_I_MAKINUM (c
);
1987 else if (SCM_BIGP (n
))
1989 unsigned long count
;
1990 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1991 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1993 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1994 scm_remember_upto_here_1 (n
);
1995 return SCM_I_MAKINUM (count
);
1998 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2003 static const char scm_ilentab
[] = {
2004 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2008 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2010 "Return the number of bits necessary to represent @var{n}.\n"
2013 "(integer-length #b10101010)\n"
2015 "(integer-length 0)\n"
2017 "(integer-length #b1111)\n"
2020 #define FUNC_NAME s_scm_integer_length
2022 if (SCM_I_INUMP (n
))
2024 unsigned long int c
= 0;
2026 long int nn
= SCM_I_INUM (n
);
2032 l
= scm_ilentab
[15 & nn
];
2035 return SCM_I_MAKINUM (c
- 4 + l
);
2037 else if (SCM_BIGP (n
))
2039 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2040 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2041 1 too big, so check for that and adjust. */
2042 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2043 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2044 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2045 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2047 scm_remember_upto_here_1 (n
);
2048 return SCM_I_MAKINUM (size
);
2051 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2055 /*** NUMBERS -> STRINGS ***/
2056 #define SCM_MAX_DBL_PREC 60
2057 #define SCM_MAX_DBL_RADIX 36
2059 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2060 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2061 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2064 void init_dblprec(int *prec
, int radix
) {
2065 /* determine floating point precision by adding successively
2066 smaller increments to 1.0 until it is considered == 1.0 */
2067 double f
= ((double)1.0)/radix
;
2068 double fsum
= 1.0 + f
;
2073 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2085 void init_fx_radix(double *fx_list
, int radix
)
2087 /* initialize a per-radix list of tolerances. When added
2088 to a number < 1.0, we can determine if we should raund
2089 up and quit converting a number to a string. */
2093 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2094 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2097 /* use this array as a way to generate a single digit */
2098 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2101 idbl2str (double f
, char *a
, int radix
)
2103 int efmt
, dpt
, d
, i
, wp
;
2105 #ifdef DBL_MIN_10_EXP
2108 #endif /* DBL_MIN_10_EXP */
2113 radix
> SCM_MAX_DBL_RADIX
)
2115 /* revert to existing behavior */
2119 wp
= scm_dblprec
[radix
-2];
2120 fx
= fx_per_radix
[radix
-2];
2124 #ifdef HAVE_COPYSIGN
2125 double sgn
= copysign (1.0, f
);
2130 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2136 strcpy (a
, "-inf.0");
2138 strcpy (a
, "+inf.0");
2141 else if (xisnan (f
))
2143 strcpy (a
, "+nan.0");
2153 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2154 make-uniform-vector, from causing infinite loops. */
2155 /* just do the checking...if it passes, we do the conversion for our
2156 radix again below */
2163 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2171 while (f_cpy
> 10.0)
2174 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2195 if (f
+ fx
[wp
] >= radix
)
2202 /* adding 9999 makes this equivalent to abs(x) % 3 */
2203 dpt
= (exp
+ 9999) % 3;
2207 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2229 a
[ch
++] = number_chars
[d
];
2232 if (f
+ fx
[wp
] >= 1.0)
2234 a
[ch
- 1] = number_chars
[d
+1];
2246 if ((dpt
> 4) && (exp
> 6))
2248 d
= (a
[0] == '-' ? 2 : 1);
2249 for (i
= ch
++; i
> d
; i
--)
2262 if (a
[ch
- 1] == '.')
2263 a
[ch
++] = '0'; /* trailing zero */
2272 for (i
= radix
; i
<= exp
; i
*= radix
);
2273 for (i
/= radix
; i
; i
/= radix
)
2275 a
[ch
++] = number_chars
[exp
/ i
];
2284 icmplx2str (double real
, double imag
, char *str
, int radix
)
2288 i
= idbl2str (real
, str
, radix
);
2291 /* Don't output a '+' for negative numbers or for Inf and
2292 NaN. They will provide their own sign. */
2293 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2295 i
+= idbl2str (imag
, &str
[i
], radix
);
2302 iflo2str (SCM flt
, char *str
, int radix
)
2305 if (SCM_REALP (flt
))
2306 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2308 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2313 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2314 characters in the result.
2316 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2318 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2323 return scm_iuint2str (-num
, rad
, p
) + 1;
2326 return scm_iuint2str (num
, rad
, p
);
2329 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2330 characters in the result.
2332 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2334 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2338 scm_t_uintmax n
= num
;
2340 for (n
/= rad
; n
> 0; n
/= rad
)
2350 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2355 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2357 "Return a string holding the external representation of the\n"
2358 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2359 "inexact, a radix of 10 will be used.")
2360 #define FUNC_NAME s_scm_number_to_string
2364 if (SCM_UNBNDP (radix
))
2367 base
= scm_to_signed_integer (radix
, 2, 36);
2369 if (SCM_I_INUMP (n
))
2371 char num_buf
[SCM_INTBUFLEN
];
2372 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2373 return scm_from_locale_stringn (num_buf
, length
);
2375 else if (SCM_BIGP (n
))
2377 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2378 scm_remember_upto_here_1 (n
);
2379 return scm_take_locale_string (str
);
2381 else if (SCM_FRACTIONP (n
))
2383 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2384 scm_from_locale_string ("/"),
2385 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2387 else if (SCM_INEXACTP (n
))
2389 char num_buf
[FLOBUFLEN
];
2390 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2393 SCM_WRONG_TYPE_ARG (1, n
);
2398 /* These print routines used to be stubbed here so that scm_repl.c
2399 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2402 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2404 char num_buf
[FLOBUFLEN
];
2405 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2410 scm_i_print_double (double val
, SCM port
)
2412 char num_buf
[FLOBUFLEN
];
2413 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2417 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2420 char num_buf
[FLOBUFLEN
];
2421 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2426 scm_i_print_complex (double real
, double imag
, SCM port
)
2428 char num_buf
[FLOBUFLEN
];
2429 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2433 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2436 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2437 scm_lfwrite (scm_i_string_chars (str
), scm_i_string_length (str
), port
);
2438 scm_remember_upto_here_1 (str
);
2443 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2445 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2446 scm_remember_upto_here_1 (exp
);
2447 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2451 /*** END nums->strs ***/
2454 /*** STRINGS -> NUMBERS ***/
2456 /* The following functions implement the conversion from strings to numbers.
2457 * The implementation somehow follows the grammar for numbers as it is given
2458 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2459 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2460 * points should be noted about the implementation:
2461 * * Each function keeps a local index variable 'idx' that points at the
2462 * current position within the parsed string. The global index is only
2463 * updated if the function could parse the corresponding syntactic unit
2465 * * Similarly, the functions keep track of indicators of inexactness ('#',
2466 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2467 * global exactness information is only updated after each part has been
2468 * successfully parsed.
2469 * * Sequences of digits are parsed into temporary variables holding fixnums.
2470 * Only if these fixnums would overflow, the result variables are updated
2471 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2472 * the temporary variables holding the fixnums are cleared, and the process
2473 * starts over again. If for example fixnums were able to store five decimal
2474 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2475 * and the result was computed as 12345 * 100000 + 67890. In other words,
2476 * only every five digits two bignum operations were performed.
2479 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2481 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2483 /* In non ASCII-style encodings the following macro might not work. */
2484 #define XDIGIT2UINT(d) \
2485 (isdigit ((int) (unsigned char) d) \
2487 : tolower ((int) (unsigned char) d) - 'a' + 10)
2490 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2491 unsigned int radix
, enum t_exactness
*p_exactness
)
2493 unsigned int idx
= *p_idx
;
2494 unsigned int hash_seen
= 0;
2495 scm_t_bits shift
= 1;
2497 unsigned int digit_value
;
2505 if (!isxdigit ((int) (unsigned char) c
))
2507 digit_value
= XDIGIT2UINT (c
);
2508 if (digit_value
>= radix
)
2512 result
= SCM_I_MAKINUM (digit_value
);
2516 if (isxdigit ((int) (unsigned char) c
))
2520 digit_value
= XDIGIT2UINT (c
);
2521 if (digit_value
>= radix
)
2533 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2535 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2537 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2544 shift
= shift
* radix
;
2545 add
= add
* radix
+ digit_value
;
2550 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2552 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2556 *p_exactness
= INEXACT
;
2562 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2563 * covers the parts of the rules that start at a potential point. The value
2564 * of the digits up to the point have been parsed by the caller and are given
2565 * in variable result. The content of *p_exactness indicates, whether a hash
2566 * has already been seen in the digits before the point.
2569 /* In non ASCII-style encodings the following macro might not work. */
2570 #define DIGIT2UINT(d) ((d) - '0')
2573 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2574 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2576 unsigned int idx
= *p_idx
;
2577 enum t_exactness x
= *p_exactness
;
2582 if (mem
[idx
] == '.')
2584 scm_t_bits shift
= 1;
2586 unsigned int digit_value
;
2587 SCM big_shift
= SCM_I_MAKINUM (1);
2593 if (isdigit ((int) (unsigned char) c
))
2598 digit_value
= DIGIT2UINT (c
);
2609 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2611 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2612 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2614 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2622 add
= add
* 10 + digit_value
;
2628 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2629 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2630 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2633 result
= scm_divide (result
, big_shift
);
2635 /* We've seen a decimal point, thus the value is implicitly inexact. */
2647 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2674 if (!isdigit ((int) (unsigned char) c
))
2678 exponent
= DIGIT2UINT (c
);
2682 if (isdigit ((int) (unsigned char) c
))
2685 if (exponent
<= SCM_MAXEXP
)
2686 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2692 if (exponent
> SCM_MAXEXP
)
2694 size_t exp_len
= idx
- start
;
2695 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2696 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2697 scm_out_of_range ("string->number", exp_num
);
2700 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2702 result
= scm_product (result
, e
);
2704 result
= scm_divide2real (result
, e
);
2706 /* We've seen an exponent, thus the value is implicitly inexact. */
2724 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2727 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2728 unsigned int radix
, enum t_exactness
*p_exactness
)
2730 unsigned int idx
= *p_idx
;
2736 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2742 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2744 enum t_exactness x
= EXACT
;
2746 /* Cobble up the fractional part. We might want to set the
2747 NaN's mantissa from it. */
2749 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2754 if (mem
[idx
] == '.')
2758 else if (idx
+ 1 == len
)
2760 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2763 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2764 p_idx
, p_exactness
);
2768 enum t_exactness x
= EXACT
;
2771 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2772 if (scm_is_false (uinteger
))
2777 else if (mem
[idx
] == '/')
2783 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2784 if (scm_is_false (divisor
))
2787 /* both are int/big here, I assume */
2788 result
= scm_i_make_ratio (uinteger
, divisor
);
2790 else if (radix
== 10)
2792 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2793 if (scm_is_false (result
))
2804 /* When returning an inexact zero, make sure it is represented as a
2805 floating point value so that we can change its sign.
2807 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2808 result
= scm_from_double (0.0);
2814 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2817 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2818 unsigned int radix
, enum t_exactness
*p_exactness
)
2842 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2843 if (scm_is_false (ureal
))
2845 /* input must be either +i or -i */
2850 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2856 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2863 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2864 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2873 /* either +<ureal>i or -<ureal>i */
2880 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2883 /* polar input: <real>@<real>. */
2908 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2909 if (scm_is_false (angle
))
2914 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2915 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2917 result
= scm_make_polar (ureal
, angle
);
2922 /* expecting input matching <real>[+-]<ureal>?i */
2929 int sign
= (c
== '+') ? 1 : -1;
2930 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2932 if (scm_is_false (imag
))
2933 imag
= SCM_I_MAKINUM (sign
);
2934 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2935 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2939 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2946 return scm_make_rectangular (ureal
, imag
);
2955 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2957 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2960 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
2961 unsigned int default_radix
)
2963 unsigned int idx
= 0;
2964 unsigned int radix
= NO_RADIX
;
2965 enum t_exactness forced_x
= NO_EXACTNESS
;
2966 enum t_exactness implicit_x
= EXACT
;
2969 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2970 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2972 switch (mem
[idx
+ 1])
2975 if (radix
!= NO_RADIX
)
2980 if (radix
!= NO_RADIX
)
2985 if (forced_x
!= NO_EXACTNESS
)
2990 if (forced_x
!= NO_EXACTNESS
)
2995 if (radix
!= NO_RADIX
)
3000 if (radix
!= NO_RADIX
)
3010 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3011 if (radix
== NO_RADIX
)
3012 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
3014 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
3016 if (scm_is_false (result
))
3022 if (SCM_INEXACTP (result
))
3023 return scm_inexact_to_exact (result
);
3027 if (SCM_INEXACTP (result
))
3030 return scm_exact_to_inexact (result
);
3033 if (implicit_x
== INEXACT
)
3035 if (SCM_INEXACTP (result
))
3038 return scm_exact_to_inexact (result
);
3046 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3047 (SCM string
, SCM radix
),
3048 "Return a number of the maximally precise representation\n"
3049 "expressed by the given @var{string}. @var{radix} must be an\n"
3050 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3051 "is a default radix that may be overridden by an explicit radix\n"
3052 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3053 "supplied, then the default radix is 10. If string is not a\n"
3054 "syntactically valid notation for a number, then\n"
3055 "@code{string->number} returns @code{#f}.")
3056 #define FUNC_NAME s_scm_string_to_number
3060 SCM_VALIDATE_STRING (1, string
);
3062 if (SCM_UNBNDP (radix
))
3065 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3067 answer
= scm_c_locale_stringn_to_number (scm_i_string_chars (string
),
3068 scm_i_string_length (string
),
3070 scm_remember_upto_here_1 (string
);
3076 /*** END strs->nums ***/
3080 scm_bigequal (SCM x
, SCM y
)
3082 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3083 scm_remember_upto_here_2 (x
, y
);
3084 return scm_from_bool (0 == result
);
3088 scm_real_equalp (SCM x
, SCM y
)
3090 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3094 scm_complex_equalp (SCM x
, SCM y
)
3096 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3097 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3101 scm_i_fraction_equalp (SCM x
, SCM y
)
3103 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3104 SCM_FRACTION_NUMERATOR (y
)))
3105 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3106 SCM_FRACTION_DENOMINATOR (y
))))
3113 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3115 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3117 #define FUNC_NAME s_scm_number_p
3119 return scm_from_bool (SCM_NUMBERP (x
));
3123 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3125 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3126 "otherwise. Note that the sets of real, rational and integer\n"
3127 "values form subsets of the set of complex numbers, i. e. the\n"
3128 "predicate will also be fulfilled if @var{x} is a real,\n"
3129 "rational or integer number.")
3130 #define FUNC_NAME s_scm_complex_p
3132 /* all numbers are complex. */
3133 return scm_number_p (x
);
3137 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3139 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3140 "otherwise. Note that the set of integer values forms a subset of\n"
3141 "the set of real numbers, i. e. the predicate will also be\n"
3142 "fulfilled if @var{x} is an integer number.")
3143 #define FUNC_NAME s_scm_real_p
3145 /* we can't represent irrational numbers. */
3146 return scm_rational_p (x
);
3150 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3152 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3153 "otherwise. Note that the set of integer values forms a subset of\n"
3154 "the set of rational numbers, i. e. the predicate will also be\n"
3155 "fulfilled if @var{x} is an integer number.")
3156 #define FUNC_NAME s_scm_rational_p
3158 if (SCM_I_INUMP (x
))
3160 else if (SCM_IMP (x
))
3162 else if (SCM_BIGP (x
))
3164 else if (SCM_FRACTIONP (x
))
3166 else if (SCM_REALP (x
))
3167 /* due to their limited precision, all floating point numbers are
3168 rational as well. */
3175 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3177 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3179 #define FUNC_NAME s_scm_integer_p
3182 if (SCM_I_INUMP (x
))
3188 if (!SCM_INEXACTP (x
))
3190 if (SCM_COMPLEXP (x
))
3192 r
= SCM_REAL_VALUE (x
);
3193 /* +/-inf passes r==floor(r), making those #t */
3201 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3203 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3205 #define FUNC_NAME s_scm_inexact_p
3207 if (SCM_INEXACTP (x
))
3209 if (SCM_NUMBERP (x
))
3211 SCM_WRONG_TYPE_ARG (1, x
);
3216 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3217 /* "Return @code{#t} if all parameters are numerically equal." */
3219 scm_num_eq_p (SCM x
, SCM y
)
3222 if (SCM_I_INUMP (x
))
3224 long xx
= SCM_I_INUM (x
);
3225 if (SCM_I_INUMP (y
))
3227 long yy
= SCM_I_INUM (y
);
3228 return scm_from_bool (xx
== yy
);
3230 else if (SCM_BIGP (y
))
3232 else if (SCM_REALP (y
))
3234 /* On a 32-bit system an inum fits a double, we can cast the inum
3235 to a double and compare.
3237 But on a 64-bit system an inum is bigger than a double and
3238 casting it to a double (call that dxx) will round. dxx is at
3239 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3240 an integer and fits a long. So we cast yy to a long and
3241 compare with plain xx.
3243 An alternative (for any size system actually) would be to check
3244 yy is an integer (with floor) and is in range of an inum
3245 (compare against appropriate powers of 2) then test
3246 xx==(long)yy. It's just a matter of which casts/comparisons
3247 might be fastest or easiest for the cpu. */
3249 double yy
= SCM_REAL_VALUE (y
);
3250 return scm_from_bool ((double) xx
== yy
3251 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3252 || xx
== (long) yy
));
3254 else if (SCM_COMPLEXP (y
))
3255 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3256 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3257 else if (SCM_FRACTIONP (y
))
3260 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3262 else if (SCM_BIGP (x
))
3264 if (SCM_I_INUMP (y
))
3266 else if (SCM_BIGP (y
))
3268 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3269 scm_remember_upto_here_2 (x
, y
);
3270 return scm_from_bool (0 == cmp
);
3272 else if (SCM_REALP (y
))
3275 if (xisnan (SCM_REAL_VALUE (y
)))
3277 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3278 scm_remember_upto_here_1 (x
);
3279 return scm_from_bool (0 == cmp
);
3281 else if (SCM_COMPLEXP (y
))
3284 if (0.0 != SCM_COMPLEX_IMAG (y
))
3286 if (xisnan (SCM_COMPLEX_REAL (y
)))
3288 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3289 scm_remember_upto_here_1 (x
);
3290 return scm_from_bool (0 == cmp
);
3292 else if (SCM_FRACTIONP (y
))
3295 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3297 else if (SCM_REALP (x
))
3299 double xx
= SCM_REAL_VALUE (x
);
3300 if (SCM_I_INUMP (y
))
3302 /* see comments with inum/real above */
3303 long yy
= SCM_I_INUM (y
);
3304 return scm_from_bool (xx
== (double) yy
3305 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3306 || (long) xx
== yy
));
3308 else if (SCM_BIGP (y
))
3311 if (xisnan (SCM_REAL_VALUE (x
)))
3313 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3314 scm_remember_upto_here_1 (y
);
3315 return scm_from_bool (0 == cmp
);
3317 else if (SCM_REALP (y
))
3318 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3319 else if (SCM_COMPLEXP (y
))
3320 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3321 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3322 else if (SCM_FRACTIONP (y
))
3324 double xx
= SCM_REAL_VALUE (x
);
3328 return scm_from_bool (xx
< 0.0);
3329 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3333 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3335 else if (SCM_COMPLEXP (x
))
3337 if (SCM_I_INUMP (y
))
3338 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3339 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3340 else if (SCM_BIGP (y
))
3343 if (0.0 != SCM_COMPLEX_IMAG (x
))
3345 if (xisnan (SCM_COMPLEX_REAL (x
)))
3347 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3348 scm_remember_upto_here_1 (y
);
3349 return scm_from_bool (0 == cmp
);
3351 else if (SCM_REALP (y
))
3352 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3353 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3354 else if (SCM_COMPLEXP (y
))
3355 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3356 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3357 else if (SCM_FRACTIONP (y
))
3360 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3362 xx
= SCM_COMPLEX_REAL (x
);
3366 return scm_from_bool (xx
< 0.0);
3367 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3371 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3373 else if (SCM_FRACTIONP (x
))
3375 if (SCM_I_INUMP (y
))
3377 else if (SCM_BIGP (y
))
3379 else if (SCM_REALP (y
))
3381 double yy
= SCM_REAL_VALUE (y
);
3385 return scm_from_bool (0.0 < yy
);
3386 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3389 else if (SCM_COMPLEXP (y
))
3392 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3394 yy
= SCM_COMPLEX_REAL (y
);
3398 return scm_from_bool (0.0 < yy
);
3399 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3402 else if (SCM_FRACTIONP (y
))
3403 return scm_i_fraction_equalp (x
, y
);
3405 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3408 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3412 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3413 done are good for inums, but for bignums an answer can almost always be
3414 had by just examining a few high bits of the operands, as done by GMP in
3415 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3416 of the float exponent to take into account. */
3418 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3419 /* "Return @code{#t} if the list of parameters is monotonically\n"
3423 scm_less_p (SCM x
, SCM y
)
3426 if (SCM_I_INUMP (x
))
3428 long xx
= SCM_I_INUM (x
);
3429 if (SCM_I_INUMP (y
))
3431 long yy
= SCM_I_INUM (y
);
3432 return scm_from_bool (xx
< yy
);
3434 else if (SCM_BIGP (y
))
3436 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3437 scm_remember_upto_here_1 (y
);
3438 return scm_from_bool (sgn
> 0);
3440 else if (SCM_REALP (y
))
3441 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3442 else if (SCM_FRACTIONP (y
))
3444 /* "x < a/b" becomes "x*b < a" */
3446 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3447 y
= SCM_FRACTION_NUMERATOR (y
);
3451 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3453 else if (SCM_BIGP (x
))
3455 if (SCM_I_INUMP (y
))
3457 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3458 scm_remember_upto_here_1 (x
);
3459 return scm_from_bool (sgn
< 0);
3461 else if (SCM_BIGP (y
))
3463 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3464 scm_remember_upto_here_2 (x
, y
);
3465 return scm_from_bool (cmp
< 0);
3467 else if (SCM_REALP (y
))
3470 if (xisnan (SCM_REAL_VALUE (y
)))
3472 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3473 scm_remember_upto_here_1 (x
);
3474 return scm_from_bool (cmp
< 0);
3476 else if (SCM_FRACTIONP (y
))
3479 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3481 else if (SCM_REALP (x
))
3483 if (SCM_I_INUMP (y
))
3484 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3485 else if (SCM_BIGP (y
))
3488 if (xisnan (SCM_REAL_VALUE (x
)))
3490 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3491 scm_remember_upto_here_1 (y
);
3492 return scm_from_bool (cmp
> 0);
3494 else if (SCM_REALP (y
))
3495 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3496 else if (SCM_FRACTIONP (y
))
3498 double xx
= SCM_REAL_VALUE (x
);
3502 return scm_from_bool (xx
< 0.0);
3503 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3507 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3509 else if (SCM_FRACTIONP (x
))
3511 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3513 /* "a/b < y" becomes "a < y*b" */
3514 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3515 x
= SCM_FRACTION_NUMERATOR (x
);
3518 else if (SCM_REALP (y
))
3520 double yy
= SCM_REAL_VALUE (y
);
3524 return scm_from_bool (0.0 < yy
);
3525 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3528 else if (SCM_FRACTIONP (y
))
3530 /* "a/b < c/d" becomes "a*d < c*b" */
3531 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3532 SCM_FRACTION_DENOMINATOR (y
));
3533 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3534 SCM_FRACTION_DENOMINATOR (x
));
3540 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3543 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3547 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3548 /* "Return @code{#t} if the list of parameters is monotonically\n"
3551 #define FUNC_NAME s_scm_gr_p
3553 scm_gr_p (SCM x
, SCM y
)
3555 if (!SCM_NUMBERP (x
))
3556 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3557 else if (!SCM_NUMBERP (y
))
3558 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3560 return scm_less_p (y
, x
);
3565 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3566 /* "Return @code{#t} if the list of parameters is monotonically\n"
3569 #define FUNC_NAME s_scm_leq_p
3571 scm_leq_p (SCM x
, SCM y
)
3573 if (!SCM_NUMBERP (x
))
3574 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3575 else if (!SCM_NUMBERP (y
))
3576 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3577 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3580 return scm_not (scm_less_p (y
, x
));
3585 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3586 /* "Return @code{#t} if the list of parameters is monotonically\n"
3589 #define FUNC_NAME s_scm_geq_p
3591 scm_geq_p (SCM x
, SCM y
)
3593 if (!SCM_NUMBERP (x
))
3594 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3595 else if (!SCM_NUMBERP (y
))
3596 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3597 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3600 return scm_not (scm_less_p (x
, y
));
3605 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3606 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3612 if (SCM_I_INUMP (z
))
3613 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3614 else if (SCM_BIGP (z
))
3616 else if (SCM_REALP (z
))
3617 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3618 else if (SCM_COMPLEXP (z
))
3619 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3620 && SCM_COMPLEX_IMAG (z
) == 0.0);
3621 else if (SCM_FRACTIONP (z
))
3624 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3628 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3629 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3633 scm_positive_p (SCM x
)
3635 if (SCM_I_INUMP (x
))
3636 return scm_from_bool (SCM_I_INUM (x
) > 0);
3637 else if (SCM_BIGP (x
))
3639 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3640 scm_remember_upto_here_1 (x
);
3641 return scm_from_bool (sgn
> 0);
3643 else if (SCM_REALP (x
))
3644 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3645 else if (SCM_FRACTIONP (x
))
3646 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3648 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3652 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3653 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3657 scm_negative_p (SCM x
)
3659 if (SCM_I_INUMP (x
))
3660 return scm_from_bool (SCM_I_INUM (x
) < 0);
3661 else if (SCM_BIGP (x
))
3663 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3664 scm_remember_upto_here_1 (x
);
3665 return scm_from_bool (sgn
< 0);
3667 else if (SCM_REALP (x
))
3668 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3669 else if (SCM_FRACTIONP (x
))
3670 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3672 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3676 /* scm_min and scm_max return an inexact when either argument is inexact, as
3677 required by r5rs. On that basis, for exact/inexact combinations the
3678 exact is converted to inexact to compare and possibly return. This is
3679 unlike scm_less_p above which takes some trouble to preserve all bits in
3680 its test, such trouble is not required for min and max. */
3682 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3683 /* "Return the maximum of all parameter values."
3686 scm_max (SCM x
, SCM y
)
3691 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3692 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3695 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3698 if (SCM_I_INUMP (x
))
3700 long xx
= SCM_I_INUM (x
);
3701 if (SCM_I_INUMP (y
))
3703 long yy
= SCM_I_INUM (y
);
3704 return (xx
< yy
) ? y
: x
;
3706 else if (SCM_BIGP (y
))
3708 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3709 scm_remember_upto_here_1 (y
);
3710 return (sgn
< 0) ? x
: y
;
3712 else if (SCM_REALP (y
))
3715 /* if y==NaN then ">" is false and we return NaN */
3716 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3718 else if (SCM_FRACTIONP (y
))
3721 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3724 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3726 else if (SCM_BIGP (x
))
3728 if (SCM_I_INUMP (y
))
3730 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3731 scm_remember_upto_here_1 (x
);
3732 return (sgn
< 0) ? y
: x
;
3734 else if (SCM_BIGP (y
))
3736 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3737 scm_remember_upto_here_2 (x
, y
);
3738 return (cmp
> 0) ? x
: y
;
3740 else if (SCM_REALP (y
))
3742 /* if y==NaN then xx>yy is false, so we return the NaN y */
3745 xx
= scm_i_big2dbl (x
);
3746 yy
= SCM_REAL_VALUE (y
);
3747 return (xx
> yy
? scm_from_double (xx
) : y
);
3749 else if (SCM_FRACTIONP (y
))
3754 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3756 else if (SCM_REALP (x
))
3758 if (SCM_I_INUMP (y
))
3760 double z
= SCM_I_INUM (y
);
3761 /* if x==NaN then "<" is false and we return NaN */
3762 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3764 else if (SCM_BIGP (y
))
3769 else if (SCM_REALP (y
))
3771 /* if x==NaN then our explicit check means we return NaN
3772 if y==NaN then ">" is false and we return NaN
3773 calling isnan is unavoidable, since it's the only way to know
3774 which of x or y causes any compares to be false */
3775 double xx
= SCM_REAL_VALUE (x
);
3776 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3778 else if (SCM_FRACTIONP (y
))
3780 double yy
= scm_i_fraction2double (y
);
3781 double xx
= SCM_REAL_VALUE (x
);
3782 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3785 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3787 else if (SCM_FRACTIONP (x
))
3789 if (SCM_I_INUMP (y
))
3793 else if (SCM_BIGP (y
))
3797 else if (SCM_REALP (y
))
3799 double xx
= scm_i_fraction2double (x
);
3800 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3802 else if (SCM_FRACTIONP (y
))
3807 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3810 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3814 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3815 /* "Return the minium of all parameter values."
3818 scm_min (SCM x
, SCM y
)
3823 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3824 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3827 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3830 if (SCM_I_INUMP (x
))
3832 long xx
= SCM_I_INUM (x
);
3833 if (SCM_I_INUMP (y
))
3835 long yy
= SCM_I_INUM (y
);
3836 return (xx
< yy
) ? x
: y
;
3838 else if (SCM_BIGP (y
))
3840 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3841 scm_remember_upto_here_1 (y
);
3842 return (sgn
< 0) ? y
: x
;
3844 else if (SCM_REALP (y
))
3847 /* if y==NaN then "<" is false and we return NaN */
3848 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3850 else if (SCM_FRACTIONP (y
))
3853 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3856 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3858 else if (SCM_BIGP (x
))
3860 if (SCM_I_INUMP (y
))
3862 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3863 scm_remember_upto_here_1 (x
);
3864 return (sgn
< 0) ? x
: y
;
3866 else if (SCM_BIGP (y
))
3868 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3869 scm_remember_upto_here_2 (x
, y
);
3870 return (cmp
> 0) ? y
: x
;
3872 else if (SCM_REALP (y
))
3874 /* if y==NaN then xx<yy is false, so we return the NaN y */
3877 xx
= scm_i_big2dbl (x
);
3878 yy
= SCM_REAL_VALUE (y
);
3879 return (xx
< yy
? scm_from_double (xx
) : y
);
3881 else if (SCM_FRACTIONP (y
))
3886 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3888 else if (SCM_REALP (x
))
3890 if (SCM_I_INUMP (y
))
3892 double z
= SCM_I_INUM (y
);
3893 /* if x==NaN then "<" is false and we return NaN */
3894 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3896 else if (SCM_BIGP (y
))
3901 else if (SCM_REALP (y
))
3903 /* if x==NaN then our explicit check means we return NaN
3904 if y==NaN then "<" is false and we return NaN
3905 calling isnan is unavoidable, since it's the only way to know
3906 which of x or y causes any compares to be false */
3907 double xx
= SCM_REAL_VALUE (x
);
3908 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3910 else if (SCM_FRACTIONP (y
))
3912 double yy
= scm_i_fraction2double (y
);
3913 double xx
= SCM_REAL_VALUE (x
);
3914 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3917 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3919 else if (SCM_FRACTIONP (x
))
3921 if (SCM_I_INUMP (y
))
3925 else if (SCM_BIGP (y
))
3929 else if (SCM_REALP (y
))
3931 double xx
= scm_i_fraction2double (x
);
3932 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3934 else if (SCM_FRACTIONP (y
))
3939 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3942 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3946 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3947 /* "Return the sum of all parameter values. Return 0 if called without\n"
3951 scm_sum (SCM x
, SCM y
)
3953 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
3955 if (SCM_NUMBERP (x
)) return x
;
3956 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3957 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3960 if (SCM_LIKELY (SCM_I_INUMP (x
)))
3962 if (SCM_LIKELY (SCM_I_INUMP (y
)))
3964 long xx
= SCM_I_INUM (x
);
3965 long yy
= SCM_I_INUM (y
);
3966 long int z
= xx
+ yy
;
3967 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3969 else if (SCM_BIGP (y
))
3974 else if (SCM_REALP (y
))
3976 long int xx
= SCM_I_INUM (x
);
3977 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3979 else if (SCM_COMPLEXP (y
))
3981 long int xx
= SCM_I_INUM (x
);
3982 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3983 SCM_COMPLEX_IMAG (y
));
3985 else if (SCM_FRACTIONP (y
))
3986 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3987 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3988 SCM_FRACTION_DENOMINATOR (y
));
3990 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3991 } else if (SCM_BIGP (x
))
3993 if (SCM_I_INUMP (y
))
3998 inum
= SCM_I_INUM (y
);
4001 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4004 SCM result
= scm_i_mkbig ();
4005 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4006 scm_remember_upto_here_1 (x
);
4007 /* we know the result will have to be a bignum */
4010 return scm_i_normbig (result
);
4014 SCM result
= scm_i_mkbig ();
4015 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4016 scm_remember_upto_here_1 (x
);
4017 /* we know the result will have to be a bignum */
4020 return scm_i_normbig (result
);
4023 else if (SCM_BIGP (y
))
4025 SCM result
= scm_i_mkbig ();
4026 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4027 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4028 mpz_add (SCM_I_BIG_MPZ (result
),
4031 scm_remember_upto_here_2 (x
, y
);
4032 /* we know the result will have to be a bignum */
4035 return scm_i_normbig (result
);
4037 else if (SCM_REALP (y
))
4039 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4040 scm_remember_upto_here_1 (x
);
4041 return scm_from_double (result
);
4043 else if (SCM_COMPLEXP (y
))
4045 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4046 + SCM_COMPLEX_REAL (y
));
4047 scm_remember_upto_here_1 (x
);
4048 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4050 else if (SCM_FRACTIONP (y
))
4051 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4052 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4053 SCM_FRACTION_DENOMINATOR (y
));
4055 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4057 else if (SCM_REALP (x
))
4059 if (SCM_I_INUMP (y
))
4060 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4061 else if (SCM_BIGP (y
))
4063 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4064 scm_remember_upto_here_1 (y
);
4065 return scm_from_double (result
);
4067 else if (SCM_REALP (y
))
4068 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4069 else if (SCM_COMPLEXP (y
))
4070 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4071 SCM_COMPLEX_IMAG (y
));
4072 else if (SCM_FRACTIONP (y
))
4073 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4075 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4077 else if (SCM_COMPLEXP (x
))
4079 if (SCM_I_INUMP (y
))
4080 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4081 SCM_COMPLEX_IMAG (x
));
4082 else if (SCM_BIGP (y
))
4084 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4085 + SCM_COMPLEX_REAL (x
));
4086 scm_remember_upto_here_1 (y
);
4087 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4089 else if (SCM_REALP (y
))
4090 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4091 SCM_COMPLEX_IMAG (x
));
4092 else if (SCM_COMPLEXP (y
))
4093 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4094 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4095 else if (SCM_FRACTIONP (y
))
4096 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4097 SCM_COMPLEX_IMAG (x
));
4099 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4101 else if (SCM_FRACTIONP (x
))
4103 if (SCM_I_INUMP (y
))
4104 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4105 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4106 SCM_FRACTION_DENOMINATOR (x
));
4107 else if (SCM_BIGP (y
))
4108 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4109 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4110 SCM_FRACTION_DENOMINATOR (x
));
4111 else if (SCM_REALP (y
))
4112 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4113 else if (SCM_COMPLEXP (y
))
4114 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4115 SCM_COMPLEX_IMAG (y
));
4116 else if (SCM_FRACTIONP (y
))
4117 /* a/b + c/d = (ad + bc) / bd */
4118 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4119 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4120 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4122 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4125 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4129 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4131 "Return @math{@var{x}+1}.")
4132 #define FUNC_NAME s_scm_oneplus
4134 return scm_sum (x
, SCM_I_MAKINUM (1));
4139 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4140 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4141 * the sum of all but the first argument are subtracted from the first
4143 #define FUNC_NAME s_difference
4145 scm_difference (SCM x
, SCM y
)
4147 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4150 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4152 if (SCM_I_INUMP (x
))
4154 long xx
= -SCM_I_INUM (x
);
4155 if (SCM_FIXABLE (xx
))
4156 return SCM_I_MAKINUM (xx
);
4158 return scm_i_long2big (xx
);
4160 else if (SCM_BIGP (x
))
4161 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4162 bignum, but negating that gives a fixnum. */
4163 return scm_i_normbig (scm_i_clonebig (x
, 0));
4164 else if (SCM_REALP (x
))
4165 return scm_from_double (-SCM_REAL_VALUE (x
));
4166 else if (SCM_COMPLEXP (x
))
4167 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4168 -SCM_COMPLEX_IMAG (x
));
4169 else if (SCM_FRACTIONP (x
))
4170 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4171 SCM_FRACTION_DENOMINATOR (x
));
4173 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4176 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4178 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4180 long int xx
= SCM_I_INUM (x
);
4181 long int yy
= SCM_I_INUM (y
);
4182 long int z
= xx
- yy
;
4183 if (SCM_FIXABLE (z
))
4184 return SCM_I_MAKINUM (z
);
4186 return scm_i_long2big (z
);
4188 else if (SCM_BIGP (y
))
4190 /* inum-x - big-y */
4191 long xx
= SCM_I_INUM (x
);
4194 return scm_i_clonebig (y
, 0);
4197 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4198 SCM result
= scm_i_mkbig ();
4201 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4204 /* x - y == -(y + -x) */
4205 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4206 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4208 scm_remember_upto_here_1 (y
);
4210 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4211 /* we know the result will have to be a bignum */
4214 return scm_i_normbig (result
);
4217 else if (SCM_REALP (y
))
4219 long int xx
= SCM_I_INUM (x
);
4220 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4222 else if (SCM_COMPLEXP (y
))
4224 long int xx
= SCM_I_INUM (x
);
4225 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4226 - SCM_COMPLEX_IMAG (y
));
4228 else if (SCM_FRACTIONP (y
))
4229 /* a - b/c = (ac - b) / c */
4230 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4231 SCM_FRACTION_NUMERATOR (y
)),
4232 SCM_FRACTION_DENOMINATOR (y
));
4234 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4236 else if (SCM_BIGP (x
))
4238 if (SCM_I_INUMP (y
))
4240 /* big-x - inum-y */
4241 long yy
= SCM_I_INUM (y
);
4242 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4244 scm_remember_upto_here_1 (x
);
4246 return (SCM_FIXABLE (-yy
) ?
4247 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4250 SCM result
= scm_i_mkbig ();
4253 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4255 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4256 scm_remember_upto_here_1 (x
);
4258 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4259 /* we know the result will have to be a bignum */
4262 return scm_i_normbig (result
);
4265 else if (SCM_BIGP (y
))
4267 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4268 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4269 SCM result
= scm_i_mkbig ();
4270 mpz_sub (SCM_I_BIG_MPZ (result
),
4273 scm_remember_upto_here_2 (x
, y
);
4274 /* we know the result will have to be a bignum */
4275 if ((sgn_x
== 1) && (sgn_y
== -1))
4277 if ((sgn_x
== -1) && (sgn_y
== 1))
4279 return scm_i_normbig (result
);
4281 else if (SCM_REALP (y
))
4283 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4284 scm_remember_upto_here_1 (x
);
4285 return scm_from_double (result
);
4287 else if (SCM_COMPLEXP (y
))
4289 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4290 - SCM_COMPLEX_REAL (y
));
4291 scm_remember_upto_here_1 (x
);
4292 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4294 else if (SCM_FRACTIONP (y
))
4295 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4296 SCM_FRACTION_NUMERATOR (y
)),
4297 SCM_FRACTION_DENOMINATOR (y
));
4298 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4300 else if (SCM_REALP (x
))
4302 if (SCM_I_INUMP (y
))
4303 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4304 else if (SCM_BIGP (y
))
4306 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4307 scm_remember_upto_here_1 (x
);
4308 return scm_from_double (result
);
4310 else if (SCM_REALP (y
))
4311 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4312 else if (SCM_COMPLEXP (y
))
4313 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4314 -SCM_COMPLEX_IMAG (y
));
4315 else if (SCM_FRACTIONP (y
))
4316 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4318 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4320 else if (SCM_COMPLEXP (x
))
4322 if (SCM_I_INUMP (y
))
4323 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4324 SCM_COMPLEX_IMAG (x
));
4325 else if (SCM_BIGP (y
))
4327 double real_part
= (SCM_COMPLEX_REAL (x
)
4328 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4329 scm_remember_upto_here_1 (x
);
4330 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4332 else if (SCM_REALP (y
))
4333 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4334 SCM_COMPLEX_IMAG (x
));
4335 else if (SCM_COMPLEXP (y
))
4336 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4337 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4338 else if (SCM_FRACTIONP (y
))
4339 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4340 SCM_COMPLEX_IMAG (x
));
4342 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4344 else if (SCM_FRACTIONP (x
))
4346 if (SCM_I_INUMP (y
))
4347 /* a/b - c = (a - cb) / b */
4348 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4349 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4350 SCM_FRACTION_DENOMINATOR (x
));
4351 else if (SCM_BIGP (y
))
4352 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4353 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4354 SCM_FRACTION_DENOMINATOR (x
));
4355 else if (SCM_REALP (y
))
4356 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4357 else if (SCM_COMPLEXP (y
))
4358 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4359 -SCM_COMPLEX_IMAG (y
));
4360 else if (SCM_FRACTIONP (y
))
4361 /* a/b - c/d = (ad - bc) / bd */
4362 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4363 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4364 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4366 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4369 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4374 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4376 "Return @math{@var{x}-1}.")
4377 #define FUNC_NAME s_scm_oneminus
4379 return scm_difference (x
, SCM_I_MAKINUM (1));
4384 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4385 /* "Return the product of all arguments. If called without arguments,\n"
4389 scm_product (SCM x
, SCM y
)
4391 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4394 return SCM_I_MAKINUM (1L);
4395 else if (SCM_NUMBERP (x
))
4398 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4401 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4406 xx
= SCM_I_INUM (x
);
4410 case 0: return x
; break;
4411 case 1: return y
; break;
4414 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4416 long yy
= SCM_I_INUM (y
);
4418 SCM k
= SCM_I_MAKINUM (kk
);
4419 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4423 SCM result
= scm_i_long2big (xx
);
4424 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4425 return scm_i_normbig (result
);
4428 else if (SCM_BIGP (y
))
4430 SCM result
= scm_i_mkbig ();
4431 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4432 scm_remember_upto_here_1 (y
);
4435 else if (SCM_REALP (y
))
4436 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4437 else if (SCM_COMPLEXP (y
))
4438 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4439 xx
* SCM_COMPLEX_IMAG (y
));
4440 else if (SCM_FRACTIONP (y
))
4441 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4442 SCM_FRACTION_DENOMINATOR (y
));
4444 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4446 else if (SCM_BIGP (x
))
4448 if (SCM_I_INUMP (y
))
4453 else if (SCM_BIGP (y
))
4455 SCM result
= scm_i_mkbig ();
4456 mpz_mul (SCM_I_BIG_MPZ (result
),
4459 scm_remember_upto_here_2 (x
, y
);
4462 else if (SCM_REALP (y
))
4464 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4465 scm_remember_upto_here_1 (x
);
4466 return scm_from_double (result
);
4468 else if (SCM_COMPLEXP (y
))
4470 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4471 scm_remember_upto_here_1 (x
);
4472 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4473 z
* SCM_COMPLEX_IMAG (y
));
4475 else if (SCM_FRACTIONP (y
))
4476 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4477 SCM_FRACTION_DENOMINATOR (y
));
4479 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4481 else if (SCM_REALP (x
))
4483 if (SCM_I_INUMP (y
))
4485 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4486 if (scm_is_eq (y
, SCM_INUM0
))
4488 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4490 else if (SCM_BIGP (y
))
4492 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4493 scm_remember_upto_here_1 (y
);
4494 return scm_from_double (result
);
4496 else if (SCM_REALP (y
))
4497 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4498 else if (SCM_COMPLEXP (y
))
4499 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4500 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4501 else if (SCM_FRACTIONP (y
))
4502 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4504 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4506 else if (SCM_COMPLEXP (x
))
4508 if (SCM_I_INUMP (y
))
4510 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4511 if (scm_is_eq (y
, SCM_INUM0
))
4513 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4514 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4516 else if (SCM_BIGP (y
))
4518 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4519 scm_remember_upto_here_1 (y
);
4520 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4521 z
* SCM_COMPLEX_IMAG (x
));
4523 else if (SCM_REALP (y
))
4524 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4525 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4526 else if (SCM_COMPLEXP (y
))
4528 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4529 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4530 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4531 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4533 else if (SCM_FRACTIONP (y
))
4535 double yy
= scm_i_fraction2double (y
);
4536 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4537 yy
* SCM_COMPLEX_IMAG (x
));
4540 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4542 else if (SCM_FRACTIONP (x
))
4544 if (SCM_I_INUMP (y
))
4545 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4546 SCM_FRACTION_DENOMINATOR (x
));
4547 else if (SCM_BIGP (y
))
4548 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4549 SCM_FRACTION_DENOMINATOR (x
));
4550 else if (SCM_REALP (y
))
4551 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4552 else if (SCM_COMPLEXP (y
))
4554 double xx
= scm_i_fraction2double (x
);
4555 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4556 xx
* SCM_COMPLEX_IMAG (y
));
4558 else if (SCM_FRACTIONP (y
))
4559 /* a/b * c/d = ac / bd */
4560 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4561 SCM_FRACTION_NUMERATOR (y
)),
4562 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4563 SCM_FRACTION_DENOMINATOR (y
)));
4565 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4568 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4571 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4572 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4573 #define ALLOW_DIVIDE_BY_ZERO
4574 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4577 /* The code below for complex division is adapted from the GNU
4578 libstdc++, which adapted it from f2c's libF77, and is subject to
4581 /****************************************************************
4582 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4584 Permission to use, copy, modify, and distribute this software
4585 and its documentation for any purpose and without fee is hereby
4586 granted, provided that the above copyright notice appear in all
4587 copies and that both that the copyright notice and this
4588 permission notice and warranty disclaimer appear in supporting
4589 documentation, and that the names of AT&T Bell Laboratories or
4590 Bellcore or any of their entities not be used in advertising or
4591 publicity pertaining to distribution of the software without
4592 specific, written prior permission.
4594 AT&T and Bellcore disclaim all warranties with regard to this
4595 software, including all implied warranties of merchantability
4596 and fitness. In no event shall AT&T or Bellcore be liable for
4597 any special, indirect or consequential damages or any damages
4598 whatsoever resulting from loss of use, data or profits, whether
4599 in an action of contract, negligence or other tortious action,
4600 arising out of or in connection with the use or performance of
4602 ****************************************************************/
4604 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4605 /* Divide the first argument by the product of the remaining
4606 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4608 #define FUNC_NAME s_divide
4610 scm_i_divide (SCM x
, SCM y
, int inexact
)
4614 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4617 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4618 else if (SCM_I_INUMP (x
))
4620 long xx
= SCM_I_INUM (x
);
4621 if (xx
== 1 || xx
== -1)
4623 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4625 scm_num_overflow (s_divide
);
4630 return scm_from_double (1.0 / (double) xx
);
4631 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4634 else if (SCM_BIGP (x
))
4637 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4638 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4640 else if (SCM_REALP (x
))
4642 double xx
= SCM_REAL_VALUE (x
);
4643 #ifndef ALLOW_DIVIDE_BY_ZERO
4645 scm_num_overflow (s_divide
);
4648 return scm_from_double (1.0 / xx
);
4650 else if (SCM_COMPLEXP (x
))
4652 double r
= SCM_COMPLEX_REAL (x
);
4653 double i
= SCM_COMPLEX_IMAG (x
);
4654 if (fabs(r
) <= fabs(i
))
4657 double d
= i
* (1.0 + t
* t
);
4658 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4663 double d
= r
* (1.0 + t
* t
);
4664 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4667 else if (SCM_FRACTIONP (x
))
4668 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4669 SCM_FRACTION_NUMERATOR (x
));
4671 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4674 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4676 long xx
= SCM_I_INUM (x
);
4677 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4679 long yy
= SCM_I_INUM (y
);
4682 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4683 scm_num_overflow (s_divide
);
4685 return scm_from_double ((double) xx
/ (double) yy
);
4688 else if (xx
% yy
!= 0)
4691 return scm_from_double ((double) xx
/ (double) yy
);
4692 else return scm_i_make_ratio (x
, y
);
4697 if (SCM_FIXABLE (z
))
4698 return SCM_I_MAKINUM (z
);
4700 return scm_i_long2big (z
);
4703 else if (SCM_BIGP (y
))
4706 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4707 else return scm_i_make_ratio (x
, y
);
4709 else if (SCM_REALP (y
))
4711 double yy
= SCM_REAL_VALUE (y
);
4712 #ifndef ALLOW_DIVIDE_BY_ZERO
4714 scm_num_overflow (s_divide
);
4717 return scm_from_double ((double) xx
/ yy
);
4719 else if (SCM_COMPLEXP (y
))
4722 complex_div
: /* y _must_ be a complex number */
4724 double r
= SCM_COMPLEX_REAL (y
);
4725 double i
= SCM_COMPLEX_IMAG (y
);
4726 if (fabs(r
) <= fabs(i
))
4729 double d
= i
* (1.0 + t
* t
);
4730 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4735 double d
= r
* (1.0 + t
* t
);
4736 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4740 else if (SCM_FRACTIONP (y
))
4741 /* a / b/c = ac / b */
4742 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4743 SCM_FRACTION_NUMERATOR (y
));
4745 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4747 else if (SCM_BIGP (x
))
4749 if (SCM_I_INUMP (y
))
4751 long int yy
= SCM_I_INUM (y
);
4754 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4755 scm_num_overflow (s_divide
);
4757 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4758 scm_remember_upto_here_1 (x
);
4759 return (sgn
== 0) ? scm_nan () : scm_inf ();
4766 /* FIXME: HMM, what are the relative performance issues here?
4767 We need to test. Is it faster on average to test
4768 divisible_p, then perform whichever operation, or is it
4769 faster to perform the integer div opportunistically and
4770 switch to real if there's a remainder? For now we take the
4771 middle ground: test, then if divisible, use the faster div
4774 long abs_yy
= yy
< 0 ? -yy
: yy
;
4775 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4779 SCM result
= scm_i_mkbig ();
4780 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4781 scm_remember_upto_here_1 (x
);
4783 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4784 return scm_i_normbig (result
);
4789 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4790 else return scm_i_make_ratio (x
, y
);
4794 else if (SCM_BIGP (y
))
4796 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4799 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4800 scm_num_overflow (s_divide
);
4802 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4803 scm_remember_upto_here_1 (x
);
4804 return (sgn
== 0) ? scm_nan () : scm_inf ();
4812 /* It's easily possible for the ratio x/y to fit a double
4813 but one or both x and y be too big to fit a double,
4814 hence the use of mpq_get_d rather than converting and
4817 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
4818 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
4819 return scm_from_double (mpq_get_d (q
));
4823 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4827 SCM result
= scm_i_mkbig ();
4828 mpz_divexact (SCM_I_BIG_MPZ (result
),
4831 scm_remember_upto_here_2 (x
, y
);
4832 return scm_i_normbig (result
);
4835 return scm_i_make_ratio (x
, y
);
4839 else if (SCM_REALP (y
))
4841 double yy
= SCM_REAL_VALUE (y
);
4842 #ifndef ALLOW_DIVIDE_BY_ZERO
4844 scm_num_overflow (s_divide
);
4847 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4849 else if (SCM_COMPLEXP (y
))
4851 a
= scm_i_big2dbl (x
);
4854 else if (SCM_FRACTIONP (y
))
4855 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4856 SCM_FRACTION_NUMERATOR (y
));
4858 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4860 else if (SCM_REALP (x
))
4862 double rx
= SCM_REAL_VALUE (x
);
4863 if (SCM_I_INUMP (y
))
4865 long int yy
= SCM_I_INUM (y
);
4866 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4868 scm_num_overflow (s_divide
);
4871 return scm_from_double (rx
/ (double) yy
);
4873 else if (SCM_BIGP (y
))
4875 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4876 scm_remember_upto_here_1 (y
);
4877 return scm_from_double (rx
/ dby
);
4879 else if (SCM_REALP (y
))
4881 double yy
= SCM_REAL_VALUE (y
);
4882 #ifndef ALLOW_DIVIDE_BY_ZERO
4884 scm_num_overflow (s_divide
);
4887 return scm_from_double (rx
/ yy
);
4889 else if (SCM_COMPLEXP (y
))
4894 else if (SCM_FRACTIONP (y
))
4895 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4897 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4899 else if (SCM_COMPLEXP (x
))
4901 double rx
= SCM_COMPLEX_REAL (x
);
4902 double ix
= SCM_COMPLEX_IMAG (x
);
4903 if (SCM_I_INUMP (y
))
4905 long int yy
= SCM_I_INUM (y
);
4906 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4908 scm_num_overflow (s_divide
);
4913 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4916 else if (SCM_BIGP (y
))
4918 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4919 scm_remember_upto_here_1 (y
);
4920 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4922 else if (SCM_REALP (y
))
4924 double yy
= SCM_REAL_VALUE (y
);
4925 #ifndef ALLOW_DIVIDE_BY_ZERO
4927 scm_num_overflow (s_divide
);
4930 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4932 else if (SCM_COMPLEXP (y
))
4934 double ry
= SCM_COMPLEX_REAL (y
);
4935 double iy
= SCM_COMPLEX_IMAG (y
);
4936 if (fabs(ry
) <= fabs(iy
))
4939 double d
= iy
* (1.0 + t
* t
);
4940 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4945 double d
= ry
* (1.0 + t
* t
);
4946 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4949 else if (SCM_FRACTIONP (y
))
4951 double yy
= scm_i_fraction2double (y
);
4952 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4955 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4957 else if (SCM_FRACTIONP (x
))
4959 if (SCM_I_INUMP (y
))
4961 long int yy
= SCM_I_INUM (y
);
4962 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4964 scm_num_overflow (s_divide
);
4967 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4968 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4970 else if (SCM_BIGP (y
))
4972 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4973 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4975 else if (SCM_REALP (y
))
4977 double yy
= SCM_REAL_VALUE (y
);
4978 #ifndef ALLOW_DIVIDE_BY_ZERO
4980 scm_num_overflow (s_divide
);
4983 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4985 else if (SCM_COMPLEXP (y
))
4987 a
= scm_i_fraction2double (x
);
4990 else if (SCM_FRACTIONP (y
))
4991 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4992 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4994 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4997 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5001 scm_divide (SCM x
, SCM y
)
5003 return scm_i_divide (x
, y
, 0);
5006 static SCM
scm_divide2real (SCM x
, SCM y
)
5008 return scm_i_divide (x
, y
, 1);
5014 scm_asinh (double x
)
5019 #define asinh scm_asinh
5020 return log (x
+ sqrt (x
* x
+ 1));
5023 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
5024 /* "Return the inverse hyperbolic sine of @var{x}."
5029 scm_acosh (double x
)
5034 #define acosh scm_acosh
5035 return log (x
+ sqrt (x
* x
- 1));
5038 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
5039 /* "Return the inverse hyperbolic cosine of @var{x}."
5044 scm_atanh (double x
)
5049 #define atanh scm_atanh
5050 return 0.5 * log ((1 + x
) / (1 - x
));
5053 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
5054 /* "Return the inverse hyperbolic tangent of @var{x}."
5059 scm_c_truncate (double x
)
5070 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5071 half-way case (ie. when x is an integer plus 0.5) going upwards.
5072 Then half-way cases are identified and adjusted down if the
5073 round-upwards didn't give the desired even integer.
5075 "plus_half == result" identifies a half-way case. If plus_half, which is
5076 x + 0.5, is an integer then x must be an integer plus 0.5.
5078 An odd "result" value is identified with result/2 != floor(result/2).
5079 This is done with plus_half, since that value is ready for use sooner in
5080 a pipelined cpu, and we're already requiring plus_half == result.
5082 Note however that we need to be careful when x is big and already an
5083 integer. In that case "x+0.5" may round to an adjacent integer, causing
5084 us to return such a value, incorrectly. For instance if the hardware is
5085 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5086 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5087 returned. Or if the hardware is in round-upwards mode, then other bigger
5088 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5089 representable value, 2^128+2^76 (or whatever), again incorrect.
5091 These bad roundings of x+0.5 are avoided by testing at the start whether
5092 x is already an integer. If it is then clearly that's the desired result
5093 already. And if it's not then the exponent must be small enough to allow
5094 an 0.5 to be represented, and hence added without a bad rounding. */
5097 scm_c_round (double x
)
5099 double plus_half
, result
;
5104 plus_half
= x
+ 0.5;
5105 result
= floor (plus_half
);
5106 /* Adjust so that the rounding is towards even. */
5107 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5112 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5114 "Round the number @var{x} towards zero.")
5115 #define FUNC_NAME s_scm_truncate_number
5117 if (scm_is_false (scm_negative_p (x
)))
5118 return scm_floor (x
);
5120 return scm_ceiling (x
);
5124 static SCM exactly_one_half
;
5126 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5128 "Round the number @var{x} towards the nearest integer. "
5129 "When it is exactly halfway between two integers, "
5130 "round towards the even one.")
5131 #define FUNC_NAME s_scm_round_number
5133 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5135 else if (SCM_REALP (x
))
5136 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5139 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5140 single quotient+remainder division then examining to see which way
5141 the rounding should go. */
5142 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5143 SCM result
= scm_floor (plus_half
);
5144 /* Adjust so that the rounding is towards even. */
5145 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5146 && scm_is_true (scm_odd_p (result
)))
5147 return scm_difference (result
, SCM_I_MAKINUM (1));
5154 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5156 "Round the number @var{x} towards minus infinity.")
5157 #define FUNC_NAME s_scm_floor
5159 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5161 else if (SCM_REALP (x
))
5162 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5163 else if (SCM_FRACTIONP (x
))
5165 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5166 SCM_FRACTION_DENOMINATOR (x
));
5167 if (scm_is_false (scm_negative_p (x
)))
5169 /* For positive x, rounding towards zero is correct. */
5174 /* For negative x, we need to return q-1 unless x is an
5175 integer. But fractions are never integer, per our
5177 return scm_difference (q
, SCM_I_MAKINUM (1));
5181 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5185 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5187 "Round the number @var{x} towards infinity.")
5188 #define FUNC_NAME s_scm_ceiling
5190 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5192 else if (SCM_REALP (x
))
5193 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5194 else if (SCM_FRACTIONP (x
))
5196 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5197 SCM_FRACTION_DENOMINATOR (x
));
5198 if (scm_is_false (scm_positive_p (x
)))
5200 /* For negative x, rounding towards zero is correct. */
5205 /* For positive x, we need to return q+1 unless x is an
5206 integer. But fractions are never integer, per our
5208 return scm_sum (q
, SCM_I_MAKINUM (1));
5212 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5216 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5217 /* "Return the square root of the real number @var{x}."
5219 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5220 /* "Return the absolute value of the real number @var{x}."
5222 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5223 /* "Return the @var{x}th power of e."
5225 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5226 /* "Return the natural logarithm of the real number @var{x}."
5228 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5229 /* "Return the sine of the real number @var{x}."
5231 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5232 /* "Return the cosine of the real number @var{x}."
5234 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5235 /* "Return the tangent of the real number @var{x}."
5237 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5238 /* "Return the arc sine of the real number @var{x}."
5240 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5241 /* "Return the arc cosine of the real number @var{x}."
5243 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5244 /* "Return the arc tangent of the real number @var{x}."
5246 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5247 /* "Return the hyperbolic sine of the real number @var{x}."
5249 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5250 /* "Return the hyperbolic cosine of the real number @var{x}."
5252 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5253 /* "Return the hyperbolic tangent of the real number @var{x}."
5261 static void scm_two_doubles (SCM x
,
5263 const char *sstring
,
5267 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5269 if (SCM_I_INUMP (x
))
5270 xy
->x
= SCM_I_INUM (x
);
5271 else if (SCM_BIGP (x
))
5272 xy
->x
= scm_i_big2dbl (x
);
5273 else if (SCM_REALP (x
))
5274 xy
->x
= SCM_REAL_VALUE (x
);
5275 else if (SCM_FRACTIONP (x
))
5276 xy
->x
= scm_i_fraction2double (x
);
5278 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5280 if (SCM_I_INUMP (y
))
5281 xy
->y
= SCM_I_INUM (y
);
5282 else if (SCM_BIGP (y
))
5283 xy
->y
= scm_i_big2dbl (y
);
5284 else if (SCM_REALP (y
))
5285 xy
->y
= SCM_REAL_VALUE (y
);
5286 else if (SCM_FRACTIONP (y
))
5287 xy
->y
= scm_i_fraction2double (y
);
5289 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5293 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5295 "Return @var{x} raised to the power of @var{y}. This\n"
5296 "procedure does not accept complex arguments.")
5297 #define FUNC_NAME s_scm_sys_expt
5300 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5301 return scm_from_double (pow (xy
.x
, xy
.y
));
5306 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5308 "Return the arc tangent of the two arguments @var{x} and\n"
5309 "@var{y}. This is similar to calculating the arc tangent of\n"
5310 "@var{x} / @var{y}, except that the signs of both arguments\n"
5311 "are used to determine the quadrant of the result. This\n"
5312 "procedure does not accept complex arguments.")
5313 #define FUNC_NAME s_scm_sys_atan2
5316 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5317 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5322 scm_c_make_rectangular (double re
, double im
)
5325 return scm_from_double (re
);
5329 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5331 SCM_COMPLEX_REAL (z
) = re
;
5332 SCM_COMPLEX_IMAG (z
) = im
;
5337 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5338 (SCM real
, SCM imaginary
),
5339 "Return a complex number constructed of the given @var{real} and\n"
5340 "@var{imaginary} parts.")
5341 #define FUNC_NAME s_scm_make_rectangular
5344 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5345 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5350 scm_c_make_polar (double mag
, double ang
)
5354 sincos (ang
, &s
, &c
);
5359 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5362 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5364 "Return the complex number @var{x} * e^(i * @var{y}).")
5365 #define FUNC_NAME s_scm_make_polar
5368 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5369 return scm_c_make_polar (xy
.x
, xy
.y
);
5374 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5375 /* "Return the real part of the number @var{z}."
5378 scm_real_part (SCM z
)
5380 if (SCM_I_INUMP (z
))
5382 else if (SCM_BIGP (z
))
5384 else if (SCM_REALP (z
))
5386 else if (SCM_COMPLEXP (z
))
5387 return scm_from_double (SCM_COMPLEX_REAL (z
));
5388 else if (SCM_FRACTIONP (z
))
5391 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5395 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5396 /* "Return the imaginary part of the number @var{z}."
5399 scm_imag_part (SCM z
)
5401 if (SCM_I_INUMP (z
))
5403 else if (SCM_BIGP (z
))
5405 else if (SCM_REALP (z
))
5407 else if (SCM_COMPLEXP (z
))
5408 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5409 else if (SCM_FRACTIONP (z
))
5412 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5415 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5416 /* "Return the numerator of the number @var{z}."
5419 scm_numerator (SCM z
)
5421 if (SCM_I_INUMP (z
))
5423 else if (SCM_BIGP (z
))
5425 else if (SCM_FRACTIONP (z
))
5426 return SCM_FRACTION_NUMERATOR (z
);
5427 else if (SCM_REALP (z
))
5428 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5430 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5434 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5435 /* "Return the denominator of the number @var{z}."
5438 scm_denominator (SCM z
)
5440 if (SCM_I_INUMP (z
))
5441 return SCM_I_MAKINUM (1);
5442 else if (SCM_BIGP (z
))
5443 return SCM_I_MAKINUM (1);
5444 else if (SCM_FRACTIONP (z
))
5445 return SCM_FRACTION_DENOMINATOR (z
);
5446 else if (SCM_REALP (z
))
5447 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5449 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5452 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5453 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5454 * "@code{abs} for real arguments, but also allows complex numbers."
5457 scm_magnitude (SCM z
)
5459 if (SCM_I_INUMP (z
))
5461 long int zz
= SCM_I_INUM (z
);
5464 else if (SCM_POSFIXABLE (-zz
))
5465 return SCM_I_MAKINUM (-zz
);
5467 return scm_i_long2big (-zz
);
5469 else if (SCM_BIGP (z
))
5471 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5472 scm_remember_upto_here_1 (z
);
5474 return scm_i_clonebig (z
, 0);
5478 else if (SCM_REALP (z
))
5479 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5480 else if (SCM_COMPLEXP (z
))
5481 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5482 else if (SCM_FRACTIONP (z
))
5484 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5486 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5487 SCM_FRACTION_DENOMINATOR (z
));
5490 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5494 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5495 /* "Return the angle of the complex number @var{z}."
5500 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5501 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5502 But if atan2 follows the floating point rounding mode, then the value
5503 is not a constant. Maybe it'd be close enough though. */
5504 if (SCM_I_INUMP (z
))
5506 if (SCM_I_INUM (z
) >= 0)
5509 return scm_from_double (atan2 (0.0, -1.0));
5511 else if (SCM_BIGP (z
))
5513 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5514 scm_remember_upto_here_1 (z
);
5516 return scm_from_double (atan2 (0.0, -1.0));
5520 else if (SCM_REALP (z
))
5522 if (SCM_REAL_VALUE (z
) >= 0)
5525 return scm_from_double (atan2 (0.0, -1.0));
5527 else if (SCM_COMPLEXP (z
))
5528 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5529 else if (SCM_FRACTIONP (z
))
5531 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5533 else return scm_from_double (atan2 (0.0, -1.0));
5536 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5540 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5541 /* Convert the number @var{x} to its inexact representation.\n"
5544 scm_exact_to_inexact (SCM z
)
5546 if (SCM_I_INUMP (z
))
5547 return scm_from_double ((double) SCM_I_INUM (z
));
5548 else if (SCM_BIGP (z
))
5549 return scm_from_double (scm_i_big2dbl (z
));
5550 else if (SCM_FRACTIONP (z
))
5551 return scm_from_double (scm_i_fraction2double (z
));
5552 else if (SCM_INEXACTP (z
))
5555 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5559 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5561 "Return an exact number that is numerically closest to @var{z}.")
5562 #define FUNC_NAME s_scm_inexact_to_exact
5564 if (SCM_I_INUMP (z
))
5566 else if (SCM_BIGP (z
))
5568 else if (SCM_REALP (z
))
5570 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5571 SCM_OUT_OF_RANGE (1, z
);
5578 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5579 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5580 scm_i_mpz2num (mpq_denref (frac
)));
5582 /* When scm_i_make_ratio throws, we leak the memory allocated
5589 else if (SCM_FRACTIONP (z
))
5592 SCM_WRONG_TYPE_ARG (1, z
);
5596 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5598 "Return an exact number that is within @var{err} of @var{x}.")
5599 #define FUNC_NAME s_scm_rationalize
5601 if (SCM_I_INUMP (x
))
5603 else if (SCM_BIGP (x
))
5605 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5607 /* Use continued fractions to find closest ratio. All
5608 arithmetic is done with exact numbers.
5611 SCM ex
= scm_inexact_to_exact (x
);
5612 SCM int_part
= scm_floor (ex
);
5613 SCM tt
= SCM_I_MAKINUM (1);
5614 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5615 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5619 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5622 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5623 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5625 /* We stop after a million iterations just to be absolutely sure
5626 that we don't go into an infinite loop. The process normally
5627 converges after less than a dozen iterations.
5630 err
= scm_abs (err
);
5631 while (++i
< 1000000)
5633 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5634 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5635 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5637 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5638 err
))) /* abs(x-a/b) <= err */
5640 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5641 if (scm_is_false (scm_exact_p (x
))
5642 || scm_is_false (scm_exact_p (err
)))
5643 return scm_exact_to_inexact (res
);
5647 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5649 tt
= scm_floor (rx
); /* tt = floor (rx) */
5655 scm_num_overflow (s_scm_rationalize
);
5658 SCM_WRONG_TYPE_ARG (1, x
);
5662 /* conversion functions */
5665 scm_is_integer (SCM val
)
5667 return scm_is_true (scm_integer_p (val
));
5671 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5673 if (SCM_I_INUMP (val
))
5675 scm_t_signed_bits n
= SCM_I_INUM (val
);
5676 return n
>= min
&& n
<= max
;
5678 else if (SCM_BIGP (val
))
5680 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5682 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5684 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5686 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5687 return n
>= min
&& n
<= max
;
5697 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5698 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5701 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5702 SCM_I_BIG_MPZ (val
));
5704 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5716 return n
>= min
&& n
<= max
;
5724 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5726 if (SCM_I_INUMP (val
))
5728 scm_t_signed_bits n
= SCM_I_INUM (val
);
5729 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5731 else if (SCM_BIGP (val
))
5733 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5735 else if (max
<= ULONG_MAX
)
5737 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5739 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5740 return n
>= min
&& n
<= max
;
5750 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5753 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5754 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5757 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5758 SCM_I_BIG_MPZ (val
));
5760 return n
>= min
&& n
<= max
;
5768 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5770 scm_error (scm_out_of_range_key
,
5772 "Value out of range ~S to ~S: ~S",
5773 scm_list_3 (min
, max
, bad_val
),
5774 scm_list_1 (bad_val
));
5777 #define TYPE scm_t_intmax
5778 #define TYPE_MIN min
5779 #define TYPE_MAX max
5780 #define SIZEOF_TYPE 0
5781 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5782 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5783 #include "libguile/conv-integer.i.c"
5785 #define TYPE scm_t_uintmax
5786 #define TYPE_MIN min
5787 #define TYPE_MAX max
5788 #define SIZEOF_TYPE 0
5789 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5790 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5791 #include "libguile/conv-uinteger.i.c"
5793 #define TYPE scm_t_int8
5794 #define TYPE_MIN SCM_T_INT8_MIN
5795 #define TYPE_MAX SCM_T_INT8_MAX
5796 #define SIZEOF_TYPE 1
5797 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5798 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5799 #include "libguile/conv-integer.i.c"
5801 #define TYPE scm_t_uint8
5803 #define TYPE_MAX SCM_T_UINT8_MAX
5804 #define SIZEOF_TYPE 1
5805 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5806 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5807 #include "libguile/conv-uinteger.i.c"
5809 #define TYPE scm_t_int16
5810 #define TYPE_MIN SCM_T_INT16_MIN
5811 #define TYPE_MAX SCM_T_INT16_MAX
5812 #define SIZEOF_TYPE 2
5813 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5814 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5815 #include "libguile/conv-integer.i.c"
5817 #define TYPE scm_t_uint16
5819 #define TYPE_MAX SCM_T_UINT16_MAX
5820 #define SIZEOF_TYPE 2
5821 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5822 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5823 #include "libguile/conv-uinteger.i.c"
5825 #define TYPE scm_t_int32
5826 #define TYPE_MIN SCM_T_INT32_MIN
5827 #define TYPE_MAX SCM_T_INT32_MAX
5828 #define SIZEOF_TYPE 4
5829 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5830 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5831 #include "libguile/conv-integer.i.c"
5833 #define TYPE scm_t_uint32
5835 #define TYPE_MAX SCM_T_UINT32_MAX
5836 #define SIZEOF_TYPE 4
5837 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5838 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5839 #include "libguile/conv-uinteger.i.c"
5841 #if SCM_HAVE_T_INT64
5843 #define TYPE scm_t_int64
5844 #define TYPE_MIN SCM_T_INT64_MIN
5845 #define TYPE_MAX SCM_T_INT64_MAX
5846 #define SIZEOF_TYPE 8
5847 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5848 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5849 #include "libguile/conv-integer.i.c"
5851 #define TYPE scm_t_uint64
5853 #define TYPE_MAX SCM_T_UINT64_MAX
5854 #define SIZEOF_TYPE 8
5855 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5856 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5857 #include "libguile/conv-uinteger.i.c"
5862 scm_to_mpz (SCM val
, mpz_t rop
)
5864 if (SCM_I_INUMP (val
))
5865 mpz_set_si (rop
, SCM_I_INUM (val
));
5866 else if (SCM_BIGP (val
))
5867 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5869 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5873 scm_from_mpz (mpz_t val
)
5875 return scm_i_mpz2num (val
);
5879 scm_is_real (SCM val
)
5881 return scm_is_true (scm_real_p (val
));
5885 scm_is_rational (SCM val
)
5887 return scm_is_true (scm_rational_p (val
));
5891 scm_to_double (SCM val
)
5893 if (SCM_I_INUMP (val
))
5894 return SCM_I_INUM (val
);
5895 else if (SCM_BIGP (val
))
5896 return scm_i_big2dbl (val
);
5897 else if (SCM_FRACTIONP (val
))
5898 return scm_i_fraction2double (val
);
5899 else if (SCM_REALP (val
))
5900 return SCM_REAL_VALUE (val
);
5902 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5906 scm_from_double (double val
)
5908 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5909 SCM_REAL_VALUE (z
) = val
;
5913 #if SCM_ENABLE_DISCOURAGED == 1
5916 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5920 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5924 scm_out_of_range (NULL
, num
);
5927 return scm_to_double (num
);
5931 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5935 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5939 scm_out_of_range (NULL
, num
);
5942 return scm_to_double (num
);
5948 scm_is_complex (SCM val
)
5950 return scm_is_true (scm_complex_p (val
));
5954 scm_c_real_part (SCM z
)
5956 if (SCM_COMPLEXP (z
))
5957 return SCM_COMPLEX_REAL (z
);
5960 /* Use the scm_real_part to get proper error checking and
5963 return scm_to_double (scm_real_part (z
));
5968 scm_c_imag_part (SCM z
)
5970 if (SCM_COMPLEXP (z
))
5971 return SCM_COMPLEX_IMAG (z
);
5974 /* Use the scm_imag_part to get proper error checking and
5975 dispatching. The result will almost always be 0.0, but not
5978 return scm_to_double (scm_imag_part (z
));
5983 scm_c_magnitude (SCM z
)
5985 return scm_to_double (scm_magnitude (z
));
5991 return scm_to_double (scm_angle (z
));
5995 scm_is_number (SCM z
)
5997 return scm_is_true (scm_number_p (z
));
6001 /* In the following functions we dispatch to the real-arg funcs like log()
6002 when we know the arg is real, instead of just handing everything to
6003 clog() for instance. This is in case clog() doesn't optimize for a
6004 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6005 well use it to go straight to the applicable C func. */
6007 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6009 "Return the natural logarithm of @var{z}.")
6010 #define FUNC_NAME s_scm_log
6012 if (SCM_COMPLEXP (z
))
6014 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG
6015 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6017 double re
= SCM_COMPLEX_REAL (z
);
6018 double im
= SCM_COMPLEX_IMAG (z
);
6019 return scm_c_make_rectangular (log (hypot (re
, im
)),
6025 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6026 although the value itself overflows. */
6027 double re
= scm_to_double (z
);
6028 double l
= log (fabs (re
));
6030 return scm_from_double (l
);
6032 return scm_c_make_rectangular (l
, M_PI
);
6038 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6040 "Return the base 10 logarithm of @var{z}.")
6041 #define FUNC_NAME s_scm_log10
6043 if (SCM_COMPLEXP (z
))
6045 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6046 clog() and a multiply by M_LOG10E, rather than the fallback
6047 log10+hypot+atan2.) */
6048 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10
6049 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6051 double re
= SCM_COMPLEX_REAL (z
);
6052 double im
= SCM_COMPLEX_IMAG (z
);
6053 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6054 M_LOG10E
* atan2 (im
, re
));
6059 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6060 although the value itself overflows. */
6061 double re
= scm_to_double (z
);
6062 double l
= log10 (fabs (re
));
6064 return scm_from_double (l
);
6066 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6072 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6074 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6075 "base of natural logarithms (2.71828@dots{}).")
6076 #define FUNC_NAME s_scm_exp
6078 if (SCM_COMPLEXP (z
))
6080 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP
6081 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6083 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6084 SCM_COMPLEX_IMAG (z
));
6089 /* When z is a negative bignum the conversion to double overflows,
6090 giving -infinity, but that's ok, the exp is still 0.0. */
6091 return scm_from_double (exp (scm_to_double (z
)));
6097 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6099 "Return the square root of @var{z}. Of the two possible roots\n"
6100 "(positive and negative), the one with the a positive real part\n"
6101 "is returned, or if that's zero then a positive imaginary part.\n"
6105 "(sqrt 9.0) @result{} 3.0\n"
6106 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6107 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6108 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6110 #define FUNC_NAME s_scm_sqrt
6112 if (SCM_COMPLEXP (x
))
6114 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT
6115 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6117 double re
= SCM_COMPLEX_REAL (x
);
6118 double im
= SCM_COMPLEX_IMAG (x
);
6119 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6120 0.5 * atan2 (im
, re
));
6125 double xx
= scm_to_double (x
);
6127 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6129 return scm_from_double (sqrt (xx
));
6141 mpz_init_set_si (z_negative_one
, -1);
6143 /* It may be possible to tune the performance of some algorithms by using
6144 * the following constants to avoid the creation of bignums. Please, before
6145 * using these values, remember the two rules of program optimization:
6146 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6147 scm_c_define ("most-positive-fixnum",
6148 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6149 scm_c_define ("most-negative-fixnum",
6150 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6152 scm_add_feature ("complex");
6153 scm_add_feature ("inexact");
6154 scm_flo0
= scm_from_double (0.0);
6156 /* determine floating point precision */
6157 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6159 init_dblprec(&scm_dblprec
[i
-2],i
);
6160 init_fx_radix(fx_per_radix
[i
-2],i
);
6163 /* hard code precision for base 10 if the preprocessor tells us to... */
6164 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6167 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6168 SCM_I_MAKINUM (2)));
6169 #include "libguile/numbers.x"