2a833c83c8a6ee83f2a9d6fb556b641adb13a3c1
1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006 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
;
186 SCM_C_INLINE_KEYWORD SCM
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
));
195 SCM_C_INLINE_KEYWORD SCM
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
);
204 SCM_C_INLINE_KEYWORD SCM
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
);
213 SCM_C_INLINE_KEYWORD SCM
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
));
224 SCM_C_INLINE_KEYWORD
int
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
);
234 SCM_C_INLINE_KEYWORD SCM
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. */
245 SCM_C_INLINE_KEYWORD SCM
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
);
342 SCM_C_INLINE_KEYWORD SCM
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
)
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_I_INUMP (x
))
3962 if (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
)
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_I_INUMP (x
))
4178 if (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
)
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_I_INUMP (x
))
4406 xx
= SCM_I_INUM (x
);
4410 case 0: return x
; break;
4411 case 1: return y
; break;
4414 if (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
))
4484 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4485 else if (SCM_BIGP (y
))
4487 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4488 scm_remember_upto_here_1 (y
);
4489 return scm_from_double (result
);
4491 else if (SCM_REALP (y
))
4492 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4493 else if (SCM_COMPLEXP (y
))
4494 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4495 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4496 else if (SCM_FRACTIONP (y
))
4497 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4499 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4501 else if (SCM_COMPLEXP (x
))
4503 if (SCM_I_INUMP (y
))
4504 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4505 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4506 else if (SCM_BIGP (y
))
4508 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4509 scm_remember_upto_here_1 (y
);
4510 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4511 z
* SCM_COMPLEX_IMAG (x
));
4513 else if (SCM_REALP (y
))
4514 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4515 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4516 else if (SCM_COMPLEXP (y
))
4518 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4519 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4520 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4521 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4523 else if (SCM_FRACTIONP (y
))
4525 double yy
= scm_i_fraction2double (y
);
4526 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4527 yy
* SCM_COMPLEX_IMAG (x
));
4530 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4532 else if (SCM_FRACTIONP (x
))
4534 if (SCM_I_INUMP (y
))
4535 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4536 SCM_FRACTION_DENOMINATOR (x
));
4537 else if (SCM_BIGP (y
))
4538 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4539 SCM_FRACTION_DENOMINATOR (x
));
4540 else if (SCM_REALP (y
))
4541 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4542 else if (SCM_COMPLEXP (y
))
4544 double xx
= scm_i_fraction2double (x
);
4545 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4546 xx
* SCM_COMPLEX_IMAG (y
));
4548 else if (SCM_FRACTIONP (y
))
4549 /* a/b * c/d = ac / bd */
4550 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4551 SCM_FRACTION_NUMERATOR (y
)),
4552 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4553 SCM_FRACTION_DENOMINATOR (y
)));
4555 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4558 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4561 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4562 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4563 #define ALLOW_DIVIDE_BY_ZERO
4564 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4567 /* The code below for complex division is adapted from the GNU
4568 libstdc++, which adapted it from f2c's libF77, and is subject to
4571 /****************************************************************
4572 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4574 Permission to use, copy, modify, and distribute this software
4575 and its documentation for any purpose and without fee is hereby
4576 granted, provided that the above copyright notice appear in all
4577 copies and that both that the copyright notice and this
4578 permission notice and warranty disclaimer appear in supporting
4579 documentation, and that the names of AT&T Bell Laboratories or
4580 Bellcore or any of their entities not be used in advertising or
4581 publicity pertaining to distribution of the software without
4582 specific, written prior permission.
4584 AT&T and Bellcore disclaim all warranties with regard to this
4585 software, including all implied warranties of merchantability
4586 and fitness. In no event shall AT&T or Bellcore be liable for
4587 any special, indirect or consequential damages or any damages
4588 whatsoever resulting from loss of use, data or profits, whether
4589 in an action of contract, negligence or other tortious action,
4590 arising out of or in connection with the use or performance of
4592 ****************************************************************/
4594 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4595 /* Divide the first argument by the product of the remaining
4596 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4598 #define FUNC_NAME s_divide
4600 scm_i_divide (SCM x
, SCM y
, int inexact
)
4607 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4608 else if (SCM_I_INUMP (x
))
4610 long xx
= SCM_I_INUM (x
);
4611 if (xx
== 1 || xx
== -1)
4613 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4615 scm_num_overflow (s_divide
);
4620 return scm_from_double (1.0 / (double) xx
);
4621 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4624 else if (SCM_BIGP (x
))
4627 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4628 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4630 else if (SCM_REALP (x
))
4632 double xx
= SCM_REAL_VALUE (x
);
4633 #ifndef ALLOW_DIVIDE_BY_ZERO
4635 scm_num_overflow (s_divide
);
4638 return scm_from_double (1.0 / xx
);
4640 else if (SCM_COMPLEXP (x
))
4642 double r
= SCM_COMPLEX_REAL (x
);
4643 double i
= SCM_COMPLEX_IMAG (x
);
4644 if (fabs(r
) <= fabs(i
))
4647 double d
= i
* (1.0 + t
* t
);
4648 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4653 double d
= r
* (1.0 + t
* t
);
4654 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4657 else if (SCM_FRACTIONP (x
))
4658 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4659 SCM_FRACTION_NUMERATOR (x
));
4661 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4664 if (SCM_I_INUMP (x
))
4666 long xx
= SCM_I_INUM (x
);
4667 if (SCM_I_INUMP (y
))
4669 long yy
= SCM_I_INUM (y
);
4672 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4673 scm_num_overflow (s_divide
);
4675 return scm_from_double ((double) xx
/ (double) yy
);
4678 else if (xx
% yy
!= 0)
4681 return scm_from_double ((double) xx
/ (double) yy
);
4682 else return scm_i_make_ratio (x
, y
);
4687 if (SCM_FIXABLE (z
))
4688 return SCM_I_MAKINUM (z
);
4690 return scm_i_long2big (z
);
4693 else if (SCM_BIGP (y
))
4696 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4697 else return scm_i_make_ratio (x
, y
);
4699 else if (SCM_REALP (y
))
4701 double yy
= SCM_REAL_VALUE (y
);
4702 #ifndef ALLOW_DIVIDE_BY_ZERO
4704 scm_num_overflow (s_divide
);
4707 return scm_from_double ((double) xx
/ yy
);
4709 else if (SCM_COMPLEXP (y
))
4712 complex_div
: /* y _must_ be a complex number */
4714 double r
= SCM_COMPLEX_REAL (y
);
4715 double i
= SCM_COMPLEX_IMAG (y
);
4716 if (fabs(r
) <= fabs(i
))
4719 double d
= i
* (1.0 + t
* t
);
4720 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4725 double d
= r
* (1.0 + t
* t
);
4726 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4730 else if (SCM_FRACTIONP (y
))
4731 /* a / b/c = ac / b */
4732 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4733 SCM_FRACTION_NUMERATOR (y
));
4735 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4737 else if (SCM_BIGP (x
))
4739 if (SCM_I_INUMP (y
))
4741 long int yy
= SCM_I_INUM (y
);
4744 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4745 scm_num_overflow (s_divide
);
4747 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4748 scm_remember_upto_here_1 (x
);
4749 return (sgn
== 0) ? scm_nan () : scm_inf ();
4756 /* FIXME: HMM, what are the relative performance issues here?
4757 We need to test. Is it faster on average to test
4758 divisible_p, then perform whichever operation, or is it
4759 faster to perform the integer div opportunistically and
4760 switch to real if there's a remainder? For now we take the
4761 middle ground: test, then if divisible, use the faster div
4764 long abs_yy
= yy
< 0 ? -yy
: yy
;
4765 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4769 SCM result
= scm_i_mkbig ();
4770 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4771 scm_remember_upto_here_1 (x
);
4773 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4774 return scm_i_normbig (result
);
4779 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4780 else return scm_i_make_ratio (x
, y
);
4784 else if (SCM_BIGP (y
))
4786 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4789 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4790 scm_num_overflow (s_divide
);
4792 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4793 scm_remember_upto_here_1 (x
);
4794 return (sgn
== 0) ? scm_nan () : scm_inf ();
4802 /* It's easily possible for the ratio x/y to fit a double
4803 but one or both x and y be too big to fit a double,
4804 hence the use of mpq_get_d rather than converting and
4807 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
4808 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
4809 return scm_from_double (mpq_get_d (q
));
4813 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4817 SCM result
= scm_i_mkbig ();
4818 mpz_divexact (SCM_I_BIG_MPZ (result
),
4821 scm_remember_upto_here_2 (x
, y
);
4822 return scm_i_normbig (result
);
4825 return scm_i_make_ratio (x
, y
);
4829 else if (SCM_REALP (y
))
4831 double yy
= SCM_REAL_VALUE (y
);
4832 #ifndef ALLOW_DIVIDE_BY_ZERO
4834 scm_num_overflow (s_divide
);
4837 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4839 else if (SCM_COMPLEXP (y
))
4841 a
= scm_i_big2dbl (x
);
4844 else if (SCM_FRACTIONP (y
))
4845 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4846 SCM_FRACTION_NUMERATOR (y
));
4848 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4850 else if (SCM_REALP (x
))
4852 double rx
= SCM_REAL_VALUE (x
);
4853 if (SCM_I_INUMP (y
))
4855 long int yy
= SCM_I_INUM (y
);
4856 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4858 scm_num_overflow (s_divide
);
4861 return scm_from_double (rx
/ (double) yy
);
4863 else if (SCM_BIGP (y
))
4865 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4866 scm_remember_upto_here_1 (y
);
4867 return scm_from_double (rx
/ dby
);
4869 else if (SCM_REALP (y
))
4871 double yy
= SCM_REAL_VALUE (y
);
4872 #ifndef ALLOW_DIVIDE_BY_ZERO
4874 scm_num_overflow (s_divide
);
4877 return scm_from_double (rx
/ yy
);
4879 else if (SCM_COMPLEXP (y
))
4884 else if (SCM_FRACTIONP (y
))
4885 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4887 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4889 else if (SCM_COMPLEXP (x
))
4891 double rx
= SCM_COMPLEX_REAL (x
);
4892 double ix
= SCM_COMPLEX_IMAG (x
);
4893 if (SCM_I_INUMP (y
))
4895 long int yy
= SCM_I_INUM (y
);
4896 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4898 scm_num_overflow (s_divide
);
4903 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4906 else if (SCM_BIGP (y
))
4908 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4909 scm_remember_upto_here_1 (y
);
4910 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4912 else if (SCM_REALP (y
))
4914 double yy
= SCM_REAL_VALUE (y
);
4915 #ifndef ALLOW_DIVIDE_BY_ZERO
4917 scm_num_overflow (s_divide
);
4920 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4922 else if (SCM_COMPLEXP (y
))
4924 double ry
= SCM_COMPLEX_REAL (y
);
4925 double iy
= SCM_COMPLEX_IMAG (y
);
4926 if (fabs(ry
) <= fabs(iy
))
4929 double d
= iy
* (1.0 + t
* t
);
4930 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4935 double d
= ry
* (1.0 + t
* t
);
4936 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4939 else if (SCM_FRACTIONP (y
))
4941 double yy
= scm_i_fraction2double (y
);
4942 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4945 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4947 else if (SCM_FRACTIONP (x
))
4949 if (SCM_I_INUMP (y
))
4951 long int yy
= SCM_I_INUM (y
);
4952 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4954 scm_num_overflow (s_divide
);
4957 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4958 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4960 else if (SCM_BIGP (y
))
4962 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4963 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4965 else if (SCM_REALP (y
))
4967 double yy
= SCM_REAL_VALUE (y
);
4968 #ifndef ALLOW_DIVIDE_BY_ZERO
4970 scm_num_overflow (s_divide
);
4973 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4975 else if (SCM_COMPLEXP (y
))
4977 a
= scm_i_fraction2double (x
);
4980 else if (SCM_FRACTIONP (y
))
4981 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4982 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4984 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4987 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4991 scm_divide (SCM x
, SCM y
)
4993 return scm_i_divide (x
, y
, 0);
4996 static SCM
scm_divide2real (SCM x
, SCM y
)
4998 return scm_i_divide (x
, y
, 1);
5004 scm_asinh (double x
)
5009 #define asinh scm_asinh
5010 return log (x
+ sqrt (x
* x
+ 1));
5013 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
5014 /* "Return the inverse hyperbolic sine of @var{x}."
5019 scm_acosh (double x
)
5024 #define acosh scm_acosh
5025 return log (x
+ sqrt (x
* x
- 1));
5028 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
5029 /* "Return the inverse hyperbolic cosine of @var{x}."
5034 scm_atanh (double x
)
5039 #define atanh scm_atanh
5040 return 0.5 * log ((1 + x
) / (1 - x
));
5043 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
5044 /* "Return the inverse hyperbolic tangent of @var{x}."
5049 scm_c_truncate (double x
)
5060 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5061 half-way case (ie. when x is an integer plus 0.5) going upwards.
5062 Then half-way cases are identified and adjusted down if the
5063 round-upwards didn't give the desired even integer.
5065 "plus_half == result" identifies a half-way case. If plus_half, which is
5066 x + 0.5, is an integer then x must be an integer plus 0.5.
5068 An odd "result" value is identified with result/2 != floor(result/2).
5069 This is done with plus_half, since that value is ready for use sooner in
5070 a pipelined cpu, and we're already requiring plus_half == result.
5072 Note however that we need to be careful when x is big and already an
5073 integer. In that case "x+0.5" may round to an adjacent integer, causing
5074 us to return such a value, incorrectly. For instance if the hardware is
5075 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5076 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5077 returned. Or if the hardware is in round-upwards mode, then other bigger
5078 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5079 representable value, 2^128+2^76 (or whatever), again incorrect.
5081 These bad roundings of x+0.5 are avoided by testing at the start whether
5082 x is already an integer. If it is then clearly that's the desired result
5083 already. And if it's not then the exponent must be small enough to allow
5084 an 0.5 to be represented, and hence added without a bad rounding. */
5087 scm_c_round (double x
)
5089 double plus_half
, result
;
5094 plus_half
= x
+ 0.5;
5095 result
= floor (plus_half
);
5096 /* Adjust so that the rounding is towards even. */
5097 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5102 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5104 "Round the number @var{x} towards zero.")
5105 #define FUNC_NAME s_scm_truncate_number
5107 if (scm_is_false (scm_negative_p (x
)))
5108 return scm_floor (x
);
5110 return scm_ceiling (x
);
5114 static SCM exactly_one_half
;
5116 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5118 "Round the number @var{x} towards the nearest integer. "
5119 "When it is exactly halfway between two integers, "
5120 "round towards the even one.")
5121 #define FUNC_NAME s_scm_round_number
5123 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5125 else if (SCM_REALP (x
))
5126 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5129 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5130 single quotient+remainder division then examining to see which way
5131 the rounding should go. */
5132 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5133 SCM result
= scm_floor (plus_half
);
5134 /* Adjust so that the rounding is towards even. */
5135 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5136 && scm_is_true (scm_odd_p (result
)))
5137 return scm_difference (result
, SCM_I_MAKINUM (1));
5144 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5146 "Round the number @var{x} towards minus infinity.")
5147 #define FUNC_NAME s_scm_floor
5149 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5151 else if (SCM_REALP (x
))
5152 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5153 else if (SCM_FRACTIONP (x
))
5155 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5156 SCM_FRACTION_DENOMINATOR (x
));
5157 if (scm_is_false (scm_negative_p (x
)))
5159 /* For positive x, rounding towards zero is correct. */
5164 /* For negative x, we need to return q-1 unless x is an
5165 integer. But fractions are never integer, per our
5167 return scm_difference (q
, SCM_I_MAKINUM (1));
5171 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5175 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5177 "Round the number @var{x} towards infinity.")
5178 #define FUNC_NAME s_scm_ceiling
5180 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5182 else if (SCM_REALP (x
))
5183 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5184 else if (SCM_FRACTIONP (x
))
5186 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5187 SCM_FRACTION_DENOMINATOR (x
));
5188 if (scm_is_false (scm_positive_p (x
)))
5190 /* For negative x, rounding towards zero is correct. */
5195 /* For positive x, we need to return q+1 unless x is an
5196 integer. But fractions are never integer, per our
5198 return scm_sum (q
, SCM_I_MAKINUM (1));
5202 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5206 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5207 /* "Return the square root of the real number @var{x}."
5209 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5210 /* "Return the absolute value of the real number @var{x}."
5212 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5213 /* "Return the @var{x}th power of e."
5215 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5216 /* "Return the natural logarithm of the real number @var{x}."
5218 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5219 /* "Return the sine of the real number @var{x}."
5221 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5222 /* "Return the cosine of the real number @var{x}."
5224 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5225 /* "Return the tangent of the real number @var{x}."
5227 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5228 /* "Return the arc sine of the real number @var{x}."
5230 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5231 /* "Return the arc cosine of the real number @var{x}."
5233 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5234 /* "Return the arc tangent of the real number @var{x}."
5236 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5237 /* "Return the hyperbolic sine of the real number @var{x}."
5239 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5240 /* "Return the hyperbolic cosine of the real number @var{x}."
5242 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5243 /* "Return the hyperbolic tangent of the real number @var{x}."
5251 static void scm_two_doubles (SCM x
,
5253 const char *sstring
,
5257 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5259 if (SCM_I_INUMP (x
))
5260 xy
->x
= SCM_I_INUM (x
);
5261 else if (SCM_BIGP (x
))
5262 xy
->x
= scm_i_big2dbl (x
);
5263 else if (SCM_REALP (x
))
5264 xy
->x
= SCM_REAL_VALUE (x
);
5265 else if (SCM_FRACTIONP (x
))
5266 xy
->x
= scm_i_fraction2double (x
);
5268 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5270 if (SCM_I_INUMP (y
))
5271 xy
->y
= SCM_I_INUM (y
);
5272 else if (SCM_BIGP (y
))
5273 xy
->y
= scm_i_big2dbl (y
);
5274 else if (SCM_REALP (y
))
5275 xy
->y
= SCM_REAL_VALUE (y
);
5276 else if (SCM_FRACTIONP (y
))
5277 xy
->y
= scm_i_fraction2double (y
);
5279 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5283 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5285 "Return @var{x} raised to the power of @var{y}. This\n"
5286 "procedure does not accept complex arguments.")
5287 #define FUNC_NAME s_scm_sys_expt
5290 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5291 return scm_from_double (pow (xy
.x
, xy
.y
));
5296 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5298 "Return the arc tangent of the two arguments @var{x} and\n"
5299 "@var{y}. This is similar to calculating the arc tangent of\n"
5300 "@var{x} / @var{y}, except that the signs of both arguments\n"
5301 "are used to determine the quadrant of the result. This\n"
5302 "procedure does not accept complex arguments.")
5303 #define FUNC_NAME s_scm_sys_atan2
5306 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5307 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5312 scm_c_make_rectangular (double re
, double im
)
5315 return scm_from_double (re
);
5319 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5321 SCM_COMPLEX_REAL (z
) = re
;
5322 SCM_COMPLEX_IMAG (z
) = im
;
5327 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5328 (SCM real
, SCM imaginary
),
5329 "Return a complex number constructed of the given @var{real} and\n"
5330 "@var{imaginary} parts.")
5331 #define FUNC_NAME s_scm_make_rectangular
5334 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5335 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5340 scm_c_make_polar (double mag
, double ang
)
5344 sincos (ang
, &s
, &c
);
5349 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5352 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5354 "Return the complex number @var{x} * e^(i * @var{y}).")
5355 #define FUNC_NAME s_scm_make_polar
5358 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5359 return scm_c_make_polar (xy
.x
, xy
.y
);
5364 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5365 /* "Return the real part of the number @var{z}."
5368 scm_real_part (SCM z
)
5370 if (SCM_I_INUMP (z
))
5372 else if (SCM_BIGP (z
))
5374 else if (SCM_REALP (z
))
5376 else if (SCM_COMPLEXP (z
))
5377 return scm_from_double (SCM_COMPLEX_REAL (z
));
5378 else if (SCM_FRACTIONP (z
))
5381 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5385 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5386 /* "Return the imaginary part of the number @var{z}."
5389 scm_imag_part (SCM z
)
5391 if (SCM_I_INUMP (z
))
5393 else if (SCM_BIGP (z
))
5395 else if (SCM_REALP (z
))
5397 else if (SCM_COMPLEXP (z
))
5398 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5399 else if (SCM_FRACTIONP (z
))
5402 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5405 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5406 /* "Return the numerator of the number @var{z}."
5409 scm_numerator (SCM z
)
5411 if (SCM_I_INUMP (z
))
5413 else if (SCM_BIGP (z
))
5415 else if (SCM_FRACTIONP (z
))
5416 return SCM_FRACTION_NUMERATOR (z
);
5417 else if (SCM_REALP (z
))
5418 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5420 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5424 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5425 /* "Return the denominator of the number @var{z}."
5428 scm_denominator (SCM z
)
5430 if (SCM_I_INUMP (z
))
5431 return SCM_I_MAKINUM (1);
5432 else if (SCM_BIGP (z
))
5433 return SCM_I_MAKINUM (1);
5434 else if (SCM_FRACTIONP (z
))
5435 return SCM_FRACTION_DENOMINATOR (z
);
5436 else if (SCM_REALP (z
))
5437 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5439 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5442 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5443 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5444 * "@code{abs} for real arguments, but also allows complex numbers."
5447 scm_magnitude (SCM z
)
5449 if (SCM_I_INUMP (z
))
5451 long int zz
= SCM_I_INUM (z
);
5454 else if (SCM_POSFIXABLE (-zz
))
5455 return SCM_I_MAKINUM (-zz
);
5457 return scm_i_long2big (-zz
);
5459 else if (SCM_BIGP (z
))
5461 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5462 scm_remember_upto_here_1 (z
);
5464 return scm_i_clonebig (z
, 0);
5468 else if (SCM_REALP (z
))
5469 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5470 else if (SCM_COMPLEXP (z
))
5471 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5472 else if (SCM_FRACTIONP (z
))
5474 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5476 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5477 SCM_FRACTION_DENOMINATOR (z
));
5480 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5484 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5485 /* "Return the angle of the complex number @var{z}."
5490 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5491 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5492 But if atan2 follows the floating point rounding mode, then the value
5493 is not a constant. Maybe it'd be close enough though. */
5494 if (SCM_I_INUMP (z
))
5496 if (SCM_I_INUM (z
) >= 0)
5499 return scm_from_double (atan2 (0.0, -1.0));
5501 else if (SCM_BIGP (z
))
5503 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5504 scm_remember_upto_here_1 (z
);
5506 return scm_from_double (atan2 (0.0, -1.0));
5510 else if (SCM_REALP (z
))
5512 if (SCM_REAL_VALUE (z
) >= 0)
5515 return scm_from_double (atan2 (0.0, -1.0));
5517 else if (SCM_COMPLEXP (z
))
5518 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5519 else if (SCM_FRACTIONP (z
))
5521 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5523 else return scm_from_double (atan2 (0.0, -1.0));
5526 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5530 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5531 /* Convert the number @var{x} to its inexact representation.\n"
5534 scm_exact_to_inexact (SCM z
)
5536 if (SCM_I_INUMP (z
))
5537 return scm_from_double ((double) SCM_I_INUM (z
));
5538 else if (SCM_BIGP (z
))
5539 return scm_from_double (scm_i_big2dbl (z
));
5540 else if (SCM_FRACTIONP (z
))
5541 return scm_from_double (scm_i_fraction2double (z
));
5542 else if (SCM_INEXACTP (z
))
5545 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5549 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5551 "Return an exact number that is numerically closest to @var{z}.")
5552 #define FUNC_NAME s_scm_inexact_to_exact
5554 if (SCM_I_INUMP (z
))
5556 else if (SCM_BIGP (z
))
5558 else if (SCM_REALP (z
))
5560 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5561 SCM_OUT_OF_RANGE (1, z
);
5568 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5569 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5570 scm_i_mpz2num (mpq_denref (frac
)));
5572 /* When scm_i_make_ratio throws, we leak the memory allocated
5579 else if (SCM_FRACTIONP (z
))
5582 SCM_WRONG_TYPE_ARG (1, z
);
5586 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5588 "Return an exact number that is within @var{err} of @var{x}.")
5589 #define FUNC_NAME s_scm_rationalize
5591 if (SCM_I_INUMP (x
))
5593 else if (SCM_BIGP (x
))
5595 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5597 /* Use continued fractions to find closest ratio. All
5598 arithmetic is done with exact numbers.
5601 SCM ex
= scm_inexact_to_exact (x
);
5602 SCM int_part
= scm_floor (ex
);
5603 SCM tt
= SCM_I_MAKINUM (1);
5604 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5605 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5609 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5612 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5613 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5615 /* We stop after a million iterations just to be absolutely sure
5616 that we don't go into an infinite loop. The process normally
5617 converges after less than a dozen iterations.
5620 err
= scm_abs (err
);
5621 while (++i
< 1000000)
5623 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5624 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5625 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5627 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5628 err
))) /* abs(x-a/b) <= err */
5630 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5631 if (scm_is_false (scm_exact_p (x
))
5632 || scm_is_false (scm_exact_p (err
)))
5633 return scm_exact_to_inexact (res
);
5637 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5639 tt
= scm_floor (rx
); /* tt = floor (rx) */
5645 scm_num_overflow (s_scm_rationalize
);
5648 SCM_WRONG_TYPE_ARG (1, x
);
5652 /* conversion functions */
5655 scm_is_integer (SCM val
)
5657 return scm_is_true (scm_integer_p (val
));
5661 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5663 if (SCM_I_INUMP (val
))
5665 scm_t_signed_bits n
= SCM_I_INUM (val
);
5666 return n
>= min
&& n
<= max
;
5668 else if (SCM_BIGP (val
))
5670 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5672 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5674 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5676 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5677 return n
>= min
&& n
<= max
;
5687 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5688 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5691 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5692 SCM_I_BIG_MPZ (val
));
5694 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5706 return n
>= min
&& n
<= max
;
5714 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5716 if (SCM_I_INUMP (val
))
5718 scm_t_signed_bits n
= SCM_I_INUM (val
);
5719 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5721 else if (SCM_BIGP (val
))
5723 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5725 else if (max
<= ULONG_MAX
)
5727 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5729 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5730 return n
>= min
&& n
<= max
;
5740 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5743 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5744 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5747 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5748 SCM_I_BIG_MPZ (val
));
5750 return n
>= min
&& n
<= max
;
5758 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5760 scm_error (scm_out_of_range_key
,
5762 "Value out of range ~S to ~S: ~S",
5763 scm_list_3 (min
, max
, bad_val
),
5764 scm_list_1 (bad_val
));
5767 #define TYPE scm_t_intmax
5768 #define TYPE_MIN min
5769 #define TYPE_MAX max
5770 #define SIZEOF_TYPE 0
5771 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5772 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5773 #include "libguile/conv-integer.i.c"
5775 #define TYPE scm_t_uintmax
5776 #define TYPE_MIN min
5777 #define TYPE_MAX max
5778 #define SIZEOF_TYPE 0
5779 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5780 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5781 #include "libguile/conv-uinteger.i.c"
5783 #define TYPE scm_t_int8
5784 #define TYPE_MIN SCM_T_INT8_MIN
5785 #define TYPE_MAX SCM_T_INT8_MAX
5786 #define SIZEOF_TYPE 1
5787 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5788 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5789 #include "libguile/conv-integer.i.c"
5791 #define TYPE scm_t_uint8
5793 #define TYPE_MAX SCM_T_UINT8_MAX
5794 #define SIZEOF_TYPE 1
5795 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5796 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5797 #include "libguile/conv-uinteger.i.c"
5799 #define TYPE scm_t_int16
5800 #define TYPE_MIN SCM_T_INT16_MIN
5801 #define TYPE_MAX SCM_T_INT16_MAX
5802 #define SIZEOF_TYPE 2
5803 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5804 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5805 #include "libguile/conv-integer.i.c"
5807 #define TYPE scm_t_uint16
5809 #define TYPE_MAX SCM_T_UINT16_MAX
5810 #define SIZEOF_TYPE 2
5811 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5812 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5813 #include "libguile/conv-uinteger.i.c"
5815 #define TYPE scm_t_int32
5816 #define TYPE_MIN SCM_T_INT32_MIN
5817 #define TYPE_MAX SCM_T_INT32_MAX
5818 #define SIZEOF_TYPE 4
5819 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5820 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5821 #include "libguile/conv-integer.i.c"
5823 #define TYPE scm_t_uint32
5825 #define TYPE_MAX SCM_T_UINT32_MAX
5826 #define SIZEOF_TYPE 4
5827 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5828 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5829 #include "libguile/conv-uinteger.i.c"
5831 #if SCM_HAVE_T_INT64
5833 #define TYPE scm_t_int64
5834 #define TYPE_MIN SCM_T_INT64_MIN
5835 #define TYPE_MAX SCM_T_INT64_MAX
5836 #define SIZEOF_TYPE 8
5837 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5838 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5839 #include "libguile/conv-integer.i.c"
5841 #define TYPE scm_t_uint64
5843 #define TYPE_MAX SCM_T_UINT64_MAX
5844 #define SIZEOF_TYPE 8
5845 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5846 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5847 #include "libguile/conv-uinteger.i.c"
5852 scm_to_mpz (SCM val
, mpz_t rop
)
5854 if (SCM_I_INUMP (val
))
5855 mpz_set_si (rop
, SCM_I_INUM (val
));
5856 else if (SCM_BIGP (val
))
5857 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5859 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5863 scm_from_mpz (mpz_t val
)
5865 return scm_i_mpz2num (val
);
5869 scm_is_real (SCM val
)
5871 return scm_is_true (scm_real_p (val
));
5875 scm_is_rational (SCM val
)
5877 return scm_is_true (scm_rational_p (val
));
5881 scm_to_double (SCM val
)
5883 if (SCM_I_INUMP (val
))
5884 return SCM_I_INUM (val
);
5885 else if (SCM_BIGP (val
))
5886 return scm_i_big2dbl (val
);
5887 else if (SCM_FRACTIONP (val
))
5888 return scm_i_fraction2double (val
);
5889 else if (SCM_REALP (val
))
5890 return SCM_REAL_VALUE (val
);
5892 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5896 scm_from_double (double val
)
5898 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5899 SCM_REAL_VALUE (z
) = val
;
5903 #if SCM_ENABLE_DISCOURAGED == 1
5906 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5910 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5914 scm_out_of_range (NULL
, num
);
5917 return scm_to_double (num
);
5921 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5925 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5929 scm_out_of_range (NULL
, num
);
5932 return scm_to_double (num
);
5938 scm_is_complex (SCM val
)
5940 return scm_is_true (scm_complex_p (val
));
5944 scm_c_real_part (SCM z
)
5946 if (SCM_COMPLEXP (z
))
5947 return SCM_COMPLEX_REAL (z
);
5950 /* Use the scm_real_part to get proper error checking and
5953 return scm_to_double (scm_real_part (z
));
5958 scm_c_imag_part (SCM z
)
5960 if (SCM_COMPLEXP (z
))
5961 return SCM_COMPLEX_IMAG (z
);
5964 /* Use the scm_imag_part to get proper error checking and
5965 dispatching. The result will almost always be 0.0, but not
5968 return scm_to_double (scm_imag_part (z
));
5973 scm_c_magnitude (SCM z
)
5975 return scm_to_double (scm_magnitude (z
));
5981 return scm_to_double (scm_angle (z
));
5985 scm_is_number (SCM z
)
5987 return scm_is_true (scm_number_p (z
));
5991 /* In the following functions we dispatch to the real-arg funcs like log()
5992 when we know the arg is real, instead of just handing everything to
5993 clog() for instance. This is in case clog() doesn't optimize for a
5994 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
5995 well use it to go straight to the applicable C func. */
5997 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
5999 "Return the natural logarithm of @var{z}.")
6000 #define FUNC_NAME s_scm_log
6002 if (SCM_COMPLEXP (z
))
6004 #if HAVE_COMPLEX_DOUBLE
6005 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6007 double re
= SCM_COMPLEX_REAL (z
);
6008 double im
= SCM_COMPLEX_IMAG (z
);
6009 return scm_c_make_rectangular (log (hypot (re
, im
)),
6015 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6016 although the value itself overflows. */
6017 double re
= scm_to_double (z
);
6018 double l
= log (fabs (re
));
6020 return scm_from_double (l
);
6022 return scm_c_make_rectangular (l
, M_PI
);
6028 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6030 "Return the base 10 logarithm of @var{z}.")
6031 #define FUNC_NAME s_scm_log10
6033 if (SCM_COMPLEXP (z
))
6035 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6036 clog() and a multiply by M_LOG10E, rather than the fallback
6037 log10+hypot+atan2.) */
6038 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10
6039 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6041 double re
= SCM_COMPLEX_REAL (z
);
6042 double im
= SCM_COMPLEX_IMAG (z
);
6043 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6044 M_LOG10E
* atan2 (im
, re
));
6049 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6050 although the value itself overflows. */
6051 double re
= scm_to_double (z
);
6052 double l
= log10 (fabs (re
));
6054 return scm_from_double (l
);
6056 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6062 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6064 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6065 "base of natural logarithms (2.71828@dots{}).")
6066 #define FUNC_NAME s_scm_exp
6068 if (SCM_COMPLEXP (z
))
6070 #if HAVE_COMPLEX_DOUBLE
6071 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6073 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6074 SCM_COMPLEX_IMAG (z
));
6079 /* When z is a negative bignum the conversion to double overflows,
6080 giving -infinity, but that's ok, the exp is still 0.0. */
6081 return scm_from_double (exp (scm_to_double (z
)));
6087 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6089 "Return the square root of @var{z}. Of the two possible roots\n"
6090 "(positive and negative), the one with the a positive real part\n"
6091 "is returned, or if that's zero then a positive imaginary part.\n"
6095 "(sqrt 9.0) @result{} 3.0\n"
6096 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6097 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6098 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6100 #define FUNC_NAME s_scm_sqrt
6102 if (SCM_COMPLEXP (x
))
6104 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT
6105 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6107 double re
= SCM_COMPLEX_REAL (x
);
6108 double im
= SCM_COMPLEX_IMAG (x
);
6109 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6110 0.5 * atan2 (im
, re
));
6115 double xx
= scm_to_double (x
);
6117 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6119 return scm_from_double (sqrt (xx
));
6131 mpz_init_set_si (z_negative_one
, -1);
6133 /* It may be possible to tune the performance of some algorithms by using
6134 * the following constants to avoid the creation of bignums. Please, before
6135 * using these values, remember the two rules of program optimization:
6136 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6137 scm_c_define ("most-positive-fixnum",
6138 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6139 scm_c_define ("most-negative-fixnum",
6140 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6142 scm_add_feature ("complex");
6143 scm_add_feature ("inexact");
6144 scm_flo0
= scm_from_double (0.0);
6146 /* determine floating point precision */
6147 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6149 init_dblprec(&scm_dblprec
[i
-2],i
);
6150 init_fx_radix(fx_per_radix
[i
-2],i
);
6153 /* hard code precision for base 10 if the preprocessor tells us to... */
6154 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6157 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6158 SCM_I_MAKINUM (2)));
6159 #include "libguile/numbers.x"