1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008 Free Software Foundation, Inc.
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 /* General assumptions:
24 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
25 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
26 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
27 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
28 * All objects satisfying SCM_FRACTIONP are never an integer.
33 - see if special casing bignums and reals in integer-exponent when
34 possible (to use mpz_pow and mpf_pow_ui) is faster.
36 - look in to better short-circuiting of common cases in
37 integer-expt and elsewhere.
39 - see if direct mpz operations can help in ash and elsewhere.
43 /* tell glibc (2.3) to give prototype for C99 trunc(), csqrt(), etc */
58 #include "libguile/_scm.h"
59 #include "libguile/feature.h"
60 #include "libguile/ports.h"
61 #include "libguile/root.h"
62 #include "libguile/smob.h"
63 #include "libguile/strings.h"
65 #include "libguile/validate.h"
66 #include "libguile/numbers.h"
67 #include "libguile/deprecation.h"
69 #include "libguile/eq.h"
71 #include "libguile/discouraged.h"
73 /* values per glibc, if not already defined */
75 #define M_LOG10E 0.43429448190325182765
78 #define M_PI 3.14159265358979323846
84 Wonder if this might be faster for some of our code? A switch on
85 the numtag would jump directly to the right case, and the
86 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
88 #define SCM_I_NUMTAG_NOTNUM 0
89 #define SCM_I_NUMTAG_INUM 1
90 #define SCM_I_NUMTAG_BIG scm_tc16_big
91 #define SCM_I_NUMTAG_REAL scm_tc16_real
92 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
93 #define SCM_I_NUMTAG(x) \
94 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
95 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
96 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
97 : SCM_I_NUMTAG_NOTNUM)))
99 /* the macro above will not work as is with fractions */
102 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
104 /* FLOBUFLEN is the maximum number of characters neccessary for the
105 * printed or scm_string representation of an inexact number.
107 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
110 #if ! defined (HAVE_ISNAN)
115 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
118 #if ! defined (HAVE_ISINF)
123 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
130 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
131 an explicit check. In some future gmp (don't know what version number),
132 mpz_cmp_d is supposed to do this itself. */
134 #define xmpz_cmp_d(z, d) \
135 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
137 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
140 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
141 isinf. It does have finite and isnan though, hence the use of those.
142 fpclass would be a possibility on that system too. */
146 #if defined (HAVE_ISINF)
148 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
149 return (! (finite (x
) || isnan (x
)));
158 #if defined (HAVE_ISNAN)
165 #if defined (GUILE_I)
166 #if HAVE_COMPLEX_DOUBLE
168 /* For an SCM object Z which is a complex number (ie. satisfies
169 SCM_COMPLEXP), return its value as a C level "complex double". */
170 #define SCM_COMPLEX_VALUE(z) \
171 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
173 /* Convert a C "complex double" to an SCM value. */
175 scm_from_complex_double (complex double z
)
177 return scm_c_make_rectangular (creal (z
), cimag (z
));
180 #endif /* HAVE_COMPLEX_DOUBLE */
185 static mpz_t z_negative_one
;
192 /* Return a newly created bignum. */
193 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
194 mpz_init (SCM_I_BIG_MPZ (z
));
199 scm_i_long2big (long x
)
201 /* Return a newly created bignum initialized to X. */
202 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
203 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
208 scm_i_ulong2big (unsigned long x
)
210 /* Return a newly created bignum initialized to X. */
211 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
212 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
217 scm_i_clonebig (SCM src_big
, int same_sign_p
)
219 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
220 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
221 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
223 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
228 scm_i_bigcmp (SCM x
, SCM y
)
230 /* Return neg if x < y, pos if x > y, and 0 if x == y */
231 /* presume we already know x and y are bignums */
232 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
233 scm_remember_upto_here_2 (x
, y
);
238 scm_i_dbl2big (double d
)
240 /* results are only defined if d is an integer */
241 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
242 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
246 /* Convert a integer in double representation to a SCM number. */
249 scm_i_dbl2num (double u
)
251 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
252 powers of 2, so there's no rounding when making "double" values
253 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
254 get rounded on a 64-bit machine, hence the "+1".
256 The use of floor() to force to an integer value ensures we get a
257 "numerically closest" value without depending on how a
258 double->long cast or how mpz_set_d will round. For reference,
259 double->long probably follows the hardware rounding mode,
260 mpz_set_d truncates towards zero. */
262 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
263 representable as a double? */
265 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
266 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
267 return SCM_I_MAKINUM ((long) u
);
269 return scm_i_dbl2big (u
);
272 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
273 with R5RS exact->inexact.
275 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
276 (ie. truncate towards zero), then adjust to get the closest double by
277 examining the next lower bit and adding 1 (to the absolute value) if
280 Bignums exactly half way between representable doubles are rounded to the
281 next higher absolute value (ie. away from zero). This seems like an
282 adequate interpretation of R5RS "numerically closest", and it's easier
283 and faster than a full "nearest-even" style.
285 The bit test must be done on the absolute value of the mpz_t, which means
286 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
287 negatives as twos complement.
289 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
290 following the hardware rounding mode, but applied to the absolute value
291 of the mpz_t operand. This is not what we want so we put the high
292 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
293 mpz_get_d is supposed to always truncate towards zero.
295 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
296 is a slowdown. It'd be faster to pick out the relevant high bits with
297 mpz_getlimbn if we could be bothered coding that, and if the new
298 truncating gmp doesn't come out. */
301 scm_i_big2dbl (SCM b
)
306 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
310 /* Current GMP, eg. 4.1.3, force truncation towards zero */
312 if (bits
> DBL_MANT_DIG
)
314 size_t shift
= bits
- DBL_MANT_DIG
;
315 mpz_init2 (tmp
, DBL_MANT_DIG
);
316 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
317 result
= ldexp (mpz_get_d (tmp
), shift
);
322 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
327 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
330 if (bits
> DBL_MANT_DIG
)
332 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
333 /* test bit number "pos" in absolute value */
334 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
335 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
337 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
341 scm_remember_upto_here_1 (b
);
346 scm_i_normbig (SCM b
)
348 /* convert a big back to a fixnum if it'll fit */
349 /* presume b is a bignum */
350 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
352 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
353 if (SCM_FIXABLE (val
))
354 b
= SCM_I_MAKINUM (val
);
359 static SCM_C_INLINE_KEYWORD SCM
360 scm_i_mpz2num (mpz_t b
)
362 /* convert a mpz number to a SCM number. */
363 if (mpz_fits_slong_p (b
))
365 long val
= mpz_get_si (b
);
366 if (SCM_FIXABLE (val
))
367 return SCM_I_MAKINUM (val
);
371 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
372 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
377 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
378 static SCM
scm_divide2real (SCM x
, SCM y
);
381 scm_i_make_ratio (SCM numerator
, SCM denominator
)
382 #define FUNC_NAME "make-ratio"
384 /* First make sure the arguments are proper.
386 if (SCM_I_INUMP (denominator
))
388 if (scm_is_eq (denominator
, SCM_INUM0
))
389 scm_num_overflow ("make-ratio");
390 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
395 if (!(SCM_BIGP(denominator
)))
396 SCM_WRONG_TYPE_ARG (2, denominator
);
398 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
399 SCM_WRONG_TYPE_ARG (1, numerator
);
401 /* Then flip signs so that the denominator is positive.
403 if (scm_is_true (scm_negative_p (denominator
)))
405 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
406 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
409 /* Now consider for each of the four fixnum/bignum combinations
410 whether the rational number is really an integer.
412 if (SCM_I_INUMP (numerator
))
414 long x
= SCM_I_INUM (numerator
);
415 if (scm_is_eq (numerator
, SCM_INUM0
))
417 if (SCM_I_INUMP (denominator
))
420 y
= SCM_I_INUM (denominator
);
422 return SCM_I_MAKINUM(1);
424 return SCM_I_MAKINUM (x
/ y
);
428 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
429 of that value for the denominator, as a bignum. Apart from
430 that case, abs(bignum) > abs(inum) so inum/bignum is not an
432 if (x
== SCM_MOST_NEGATIVE_FIXNUM
433 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
434 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
435 return SCM_I_MAKINUM(-1);
438 else if (SCM_BIGP (numerator
))
440 if (SCM_I_INUMP (denominator
))
442 long yy
= SCM_I_INUM (denominator
);
443 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
444 return scm_divide (numerator
, denominator
);
448 if (scm_is_eq (numerator
, denominator
))
449 return SCM_I_MAKINUM(1);
450 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
451 SCM_I_BIG_MPZ (denominator
)))
452 return scm_divide(numerator
, denominator
);
456 /* No, it's a proper fraction.
459 SCM divisor
= scm_gcd (numerator
, denominator
);
460 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
462 numerator
= scm_divide (numerator
, divisor
);
463 denominator
= scm_divide (denominator
, divisor
);
466 return scm_double_cell (scm_tc16_fraction
,
467 SCM_UNPACK (numerator
),
468 SCM_UNPACK (denominator
), 0);
474 scm_i_fraction2double (SCM z
)
476 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
477 SCM_FRACTION_DENOMINATOR (z
)));
480 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
482 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
484 #define FUNC_NAME s_scm_exact_p
490 if (SCM_FRACTIONP (x
))
494 SCM_WRONG_TYPE_ARG (1, x
);
499 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
501 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
503 #define FUNC_NAME s_scm_odd_p
507 long val
= SCM_I_INUM (n
);
508 return scm_from_bool ((val
& 1L) != 0);
510 else if (SCM_BIGP (n
))
512 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
513 scm_remember_upto_here_1 (n
);
514 return scm_from_bool (odd_p
);
516 else if (scm_is_true (scm_inf_p (n
)))
518 else if (SCM_REALP (n
))
520 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
526 SCM_WRONG_TYPE_ARG (1, n
);
529 SCM_WRONG_TYPE_ARG (1, n
);
534 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
536 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
538 #define FUNC_NAME s_scm_even_p
542 long val
= SCM_I_INUM (n
);
543 return scm_from_bool ((val
& 1L) == 0);
545 else if (SCM_BIGP (n
))
547 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
548 scm_remember_upto_here_1 (n
);
549 return scm_from_bool (even_p
);
551 else if (scm_is_true (scm_inf_p (n
)))
553 else if (SCM_REALP (n
))
555 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
561 SCM_WRONG_TYPE_ARG (1, n
);
564 SCM_WRONG_TYPE_ARG (1, n
);
568 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
570 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
571 "or @samp{-inf.0}, @code{#f} otherwise.")
572 #define FUNC_NAME s_scm_inf_p
575 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
576 else if (SCM_COMPLEXP (x
))
577 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
578 || xisinf (SCM_COMPLEX_IMAG (x
)));
584 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
586 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
588 #define FUNC_NAME s_scm_nan_p
591 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
592 else if (SCM_COMPLEXP (n
))
593 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
594 || xisnan (SCM_COMPLEX_IMAG (n
)));
600 /* Guile's idea of infinity. */
601 static double guile_Inf
;
603 /* Guile's idea of not a number. */
604 static double guile_NaN
;
607 guile_ieee_init (void)
609 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
611 /* Some version of gcc on some old version of Linux used to crash when
612 trying to make Inf and NaN. */
615 /* C99 INFINITY, when available.
616 FIXME: The standard allows for INFINITY to be something that overflows
617 at compile time. We ought to have a configure test to check for that
618 before trying to use it. (But in practice we believe this is not a
619 problem on any system guile is likely to target.) */
620 guile_Inf
= INFINITY
;
623 extern unsigned int DINFINITY
[2];
624 guile_Inf
= (*((double *) (DINFINITY
)));
631 if (guile_Inf
== tmp
)
639 #if defined (HAVE_ISNAN)
642 /* C99 NAN, when available */
647 extern unsigned int DQNAN
[2];
648 guile_NaN
= (*((double *)(DQNAN
)));
651 guile_NaN
= guile_Inf
/ guile_Inf
;
657 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
660 #define FUNC_NAME s_scm_inf
662 static int initialized
= 0;
668 return scm_from_double (guile_Inf
);
672 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
675 #define FUNC_NAME s_scm_nan
677 static int initialized
= 0;
683 return scm_from_double (guile_NaN
);
688 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
690 "Return the absolute value of @var{x}.")
695 long int xx
= SCM_I_INUM (x
);
698 else if (SCM_POSFIXABLE (-xx
))
699 return SCM_I_MAKINUM (-xx
);
701 return scm_i_long2big (-xx
);
703 else if (SCM_BIGP (x
))
705 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
707 return scm_i_clonebig (x
, 0);
711 else if (SCM_REALP (x
))
713 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
714 double xx
= SCM_REAL_VALUE (x
);
716 return scm_from_double (-xx
);
720 else if (SCM_FRACTIONP (x
))
722 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
724 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
725 SCM_FRACTION_DENOMINATOR (x
));
728 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
733 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
734 /* "Return the quotient of the numbers @var{x} and @var{y}."
737 scm_quotient (SCM x
, SCM y
)
741 long xx
= SCM_I_INUM (x
);
744 long yy
= SCM_I_INUM (y
);
746 scm_num_overflow (s_quotient
);
751 return SCM_I_MAKINUM (z
);
753 return scm_i_long2big (z
);
756 else if (SCM_BIGP (y
))
758 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
759 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
760 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
762 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
763 scm_remember_upto_here_1 (y
);
764 return SCM_I_MAKINUM (-1);
767 return SCM_I_MAKINUM (0);
770 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
772 else if (SCM_BIGP (x
))
776 long yy
= SCM_I_INUM (y
);
778 scm_num_overflow (s_quotient
);
783 SCM result
= scm_i_mkbig ();
786 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
789 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
792 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
793 scm_remember_upto_here_1 (x
);
794 return scm_i_normbig (result
);
797 else if (SCM_BIGP (y
))
799 SCM result
= scm_i_mkbig ();
800 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
803 scm_remember_upto_here_2 (x
, y
);
804 return scm_i_normbig (result
);
807 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
810 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
813 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
814 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
816 * "(remainder 13 4) @result{} 1\n"
817 * "(remainder -13 4) @result{} -1\n"
821 scm_remainder (SCM x
, SCM y
)
827 long yy
= SCM_I_INUM (y
);
829 scm_num_overflow (s_remainder
);
832 long z
= SCM_I_INUM (x
) % yy
;
833 return SCM_I_MAKINUM (z
);
836 else if (SCM_BIGP (y
))
838 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
839 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
840 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
842 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
843 scm_remember_upto_here_1 (y
);
844 return SCM_I_MAKINUM (0);
850 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
852 else if (SCM_BIGP (x
))
856 long yy
= SCM_I_INUM (y
);
858 scm_num_overflow (s_remainder
);
861 SCM result
= scm_i_mkbig ();
864 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
865 scm_remember_upto_here_1 (x
);
866 return scm_i_normbig (result
);
869 else if (SCM_BIGP (y
))
871 SCM result
= scm_i_mkbig ();
872 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
875 scm_remember_upto_here_2 (x
, y
);
876 return scm_i_normbig (result
);
879 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
882 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
886 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
887 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
889 * "(modulo 13 4) @result{} 1\n"
890 * "(modulo -13 4) @result{} 3\n"
894 scm_modulo (SCM x
, SCM y
)
898 long xx
= SCM_I_INUM (x
);
901 long yy
= SCM_I_INUM (y
);
903 scm_num_overflow (s_modulo
);
906 /* C99 specifies that "%" is the remainder corresponding to a
907 quotient rounded towards zero, and that's also traditional
908 for machine division, so z here should be well defined. */
926 return SCM_I_MAKINUM (result
);
929 else if (SCM_BIGP (y
))
931 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
938 SCM pos_y
= scm_i_clonebig (y
, 0);
939 /* do this after the last scm_op */
940 mpz_init_set_si (z_x
, xx
);
941 result
= pos_y
; /* re-use this bignum */
942 mpz_mod (SCM_I_BIG_MPZ (result
),
944 SCM_I_BIG_MPZ (pos_y
));
945 scm_remember_upto_here_1 (pos_y
);
949 result
= scm_i_mkbig ();
950 /* do this after the last scm_op */
951 mpz_init_set_si (z_x
, xx
);
952 mpz_mod (SCM_I_BIG_MPZ (result
),
955 scm_remember_upto_here_1 (y
);
958 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
959 mpz_add (SCM_I_BIG_MPZ (result
),
961 SCM_I_BIG_MPZ (result
));
962 scm_remember_upto_here_1 (y
);
963 /* and do this before the next one */
965 return scm_i_normbig (result
);
969 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
971 else if (SCM_BIGP (x
))
975 long yy
= SCM_I_INUM (y
);
977 scm_num_overflow (s_modulo
);
980 SCM result
= scm_i_mkbig ();
981 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
983 (yy
< 0) ? - yy
: yy
);
984 scm_remember_upto_here_1 (x
);
985 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
986 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
987 SCM_I_BIG_MPZ (result
),
989 return scm_i_normbig (result
);
992 else if (SCM_BIGP (y
))
995 SCM result
= scm_i_mkbig ();
996 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
997 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
998 mpz_mod (SCM_I_BIG_MPZ (result
),
1000 SCM_I_BIG_MPZ (pos_y
));
1002 scm_remember_upto_here_1 (x
);
1003 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1004 mpz_add (SCM_I_BIG_MPZ (result
),
1006 SCM_I_BIG_MPZ (result
));
1007 scm_remember_upto_here_2 (y
, pos_y
);
1008 return scm_i_normbig (result
);
1012 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1015 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1018 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
1019 /* "Return the greatest common divisor of all arguments.\n"
1020 * "If called without arguments, 0 is returned."
1023 scm_gcd (SCM x
, SCM y
)
1026 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1028 if (SCM_I_INUMP (x
))
1030 if (SCM_I_INUMP (y
))
1032 long xx
= SCM_I_INUM (x
);
1033 long yy
= SCM_I_INUM (y
);
1034 long u
= xx
< 0 ? -xx
: xx
;
1035 long v
= yy
< 0 ? -yy
: yy
;
1045 /* Determine a common factor 2^k */
1046 while (!(1 & (u
| v
)))
1052 /* Now, any factor 2^n can be eliminated */
1072 return (SCM_POSFIXABLE (result
)
1073 ? SCM_I_MAKINUM (result
)
1074 : scm_i_long2big (result
));
1076 else if (SCM_BIGP (y
))
1082 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1084 else if (SCM_BIGP (x
))
1086 if (SCM_I_INUMP (y
))
1088 unsigned long result
;
1091 yy
= SCM_I_INUM (y
);
1096 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1097 scm_remember_upto_here_1 (x
);
1098 return (SCM_POSFIXABLE (result
)
1099 ? SCM_I_MAKINUM (result
)
1100 : scm_from_ulong (result
));
1102 else if (SCM_BIGP (y
))
1104 SCM result
= scm_i_mkbig ();
1105 mpz_gcd (SCM_I_BIG_MPZ (result
),
1108 scm_remember_upto_here_2 (x
, y
);
1109 return scm_i_normbig (result
);
1112 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1115 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1118 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1119 /* "Return the least common multiple of the arguments.\n"
1120 * "If called without arguments, 1 is returned."
1123 scm_lcm (SCM n1
, SCM n2
)
1125 if (SCM_UNBNDP (n2
))
1127 if (SCM_UNBNDP (n1
))
1128 return SCM_I_MAKINUM (1L);
1129 n2
= SCM_I_MAKINUM (1L);
1132 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1133 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1134 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1135 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1137 if (SCM_I_INUMP (n1
))
1139 if (SCM_I_INUMP (n2
))
1141 SCM d
= scm_gcd (n1
, n2
);
1142 if (scm_is_eq (d
, SCM_INUM0
))
1145 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1149 /* inum n1, big n2 */
1152 SCM result
= scm_i_mkbig ();
1153 long nn1
= SCM_I_INUM (n1
);
1154 if (nn1
== 0) return SCM_INUM0
;
1155 if (nn1
< 0) nn1
= - nn1
;
1156 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1157 scm_remember_upto_here_1 (n2
);
1165 if (SCM_I_INUMP (n2
))
1172 SCM result
= scm_i_mkbig ();
1173 mpz_lcm(SCM_I_BIG_MPZ (result
),
1175 SCM_I_BIG_MPZ (n2
));
1176 scm_remember_upto_here_2(n1
, n2
);
1177 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1183 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1188 + + + x (map digit:logand X Y)
1189 + - + x (map digit:logand X (lognot (+ -1 Y)))
1190 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1191 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1196 + + + (map digit:logior X Y)
1197 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1198 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1199 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1204 + + + (map digit:logxor X Y)
1205 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1206 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1207 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1212 + + (any digit:logand X Y)
1213 + - (any digit:logand X (lognot (+ -1 Y)))
1214 - + (any digit:logand (lognot (+ -1 X)) Y)
1219 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1221 "Return the bitwise AND of the integer arguments.\n\n"
1223 "(logand) @result{} -1\n"
1224 "(logand 7) @result{} 7\n"
1225 "(logand #b111 #b011 #b001) @result{} 1\n"
1227 #define FUNC_NAME s_scm_logand
1231 if (SCM_UNBNDP (n2
))
1233 if (SCM_UNBNDP (n1
))
1234 return SCM_I_MAKINUM (-1);
1235 else if (!SCM_NUMBERP (n1
))
1236 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1237 else if (SCM_NUMBERP (n1
))
1240 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1243 if (SCM_I_INUMP (n1
))
1245 nn1
= SCM_I_INUM (n1
);
1246 if (SCM_I_INUMP (n2
))
1248 long nn2
= SCM_I_INUM (n2
);
1249 return SCM_I_MAKINUM (nn1
& nn2
);
1251 else if SCM_BIGP (n2
)
1257 SCM result_z
= scm_i_mkbig ();
1259 mpz_init_set_si (nn1_z
, nn1
);
1260 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1261 scm_remember_upto_here_1 (n2
);
1263 return scm_i_normbig (result_z
);
1267 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1269 else if (SCM_BIGP (n1
))
1271 if (SCM_I_INUMP (n2
))
1274 nn1
= SCM_I_INUM (n1
);
1277 else if (SCM_BIGP (n2
))
1279 SCM result_z
= scm_i_mkbig ();
1280 mpz_and (SCM_I_BIG_MPZ (result_z
),
1282 SCM_I_BIG_MPZ (n2
));
1283 scm_remember_upto_here_2 (n1
, n2
);
1284 return scm_i_normbig (result_z
);
1287 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1290 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1295 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1297 "Return the bitwise OR of the integer arguments.\n\n"
1299 "(logior) @result{} 0\n"
1300 "(logior 7) @result{} 7\n"
1301 "(logior #b000 #b001 #b011) @result{} 3\n"
1303 #define FUNC_NAME s_scm_logior
1307 if (SCM_UNBNDP (n2
))
1309 if (SCM_UNBNDP (n1
))
1311 else if (SCM_NUMBERP (n1
))
1314 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1317 if (SCM_I_INUMP (n1
))
1319 nn1
= SCM_I_INUM (n1
);
1320 if (SCM_I_INUMP (n2
))
1322 long nn2
= SCM_I_INUM (n2
);
1323 return SCM_I_MAKINUM (nn1
| nn2
);
1325 else if (SCM_BIGP (n2
))
1331 SCM result_z
= scm_i_mkbig ();
1333 mpz_init_set_si (nn1_z
, nn1
);
1334 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1335 scm_remember_upto_here_1 (n2
);
1337 return scm_i_normbig (result_z
);
1341 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1343 else if (SCM_BIGP (n1
))
1345 if (SCM_I_INUMP (n2
))
1348 nn1
= SCM_I_INUM (n1
);
1351 else if (SCM_BIGP (n2
))
1353 SCM result_z
= scm_i_mkbig ();
1354 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1356 SCM_I_BIG_MPZ (n2
));
1357 scm_remember_upto_here_2 (n1
, n2
);
1358 return scm_i_normbig (result_z
);
1361 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1364 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1369 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1371 "Return the bitwise XOR of the integer arguments. A bit is\n"
1372 "set in the result if it is set in an odd number of arguments.\n"
1374 "(logxor) @result{} 0\n"
1375 "(logxor 7) @result{} 7\n"
1376 "(logxor #b000 #b001 #b011) @result{} 2\n"
1377 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1379 #define FUNC_NAME s_scm_logxor
1383 if (SCM_UNBNDP (n2
))
1385 if (SCM_UNBNDP (n1
))
1387 else if (SCM_NUMBERP (n1
))
1390 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1393 if (SCM_I_INUMP (n1
))
1395 nn1
= SCM_I_INUM (n1
);
1396 if (SCM_I_INUMP (n2
))
1398 long nn2
= SCM_I_INUM (n2
);
1399 return SCM_I_MAKINUM (nn1
^ nn2
);
1401 else if (SCM_BIGP (n2
))
1405 SCM result_z
= scm_i_mkbig ();
1407 mpz_init_set_si (nn1_z
, nn1
);
1408 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1409 scm_remember_upto_here_1 (n2
);
1411 return scm_i_normbig (result_z
);
1415 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1417 else if (SCM_BIGP (n1
))
1419 if (SCM_I_INUMP (n2
))
1422 nn1
= SCM_I_INUM (n1
);
1425 else if (SCM_BIGP (n2
))
1427 SCM result_z
= scm_i_mkbig ();
1428 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1430 SCM_I_BIG_MPZ (n2
));
1431 scm_remember_upto_here_2 (n1
, n2
);
1432 return scm_i_normbig (result_z
);
1435 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1438 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1443 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1445 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1446 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1447 "without actually calculating the @code{logand}, just testing\n"
1451 "(logtest #b0100 #b1011) @result{} #f\n"
1452 "(logtest #b0100 #b0111) @result{} #t\n"
1454 #define FUNC_NAME s_scm_logtest
1458 if (SCM_I_INUMP (j
))
1460 nj
= SCM_I_INUM (j
);
1461 if (SCM_I_INUMP (k
))
1463 long nk
= SCM_I_INUM (k
);
1464 return scm_from_bool (nj
& nk
);
1466 else if (SCM_BIGP (k
))
1474 mpz_init_set_si (nj_z
, nj
);
1475 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1476 scm_remember_upto_here_1 (k
);
1477 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1483 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1485 else if (SCM_BIGP (j
))
1487 if (SCM_I_INUMP (k
))
1490 nj
= SCM_I_INUM (j
);
1493 else if (SCM_BIGP (k
))
1497 mpz_init (result_z
);
1501 scm_remember_upto_here_2 (j
, k
);
1502 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1503 mpz_clear (result_z
);
1507 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1510 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1515 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1517 "Test whether bit number @var{index} in @var{j} is set.\n"
1518 "@var{index} starts from 0 for the least significant bit.\n"
1521 "(logbit? 0 #b1101) @result{} #t\n"
1522 "(logbit? 1 #b1101) @result{} #f\n"
1523 "(logbit? 2 #b1101) @result{} #t\n"
1524 "(logbit? 3 #b1101) @result{} #t\n"
1525 "(logbit? 4 #b1101) @result{} #f\n"
1527 #define FUNC_NAME s_scm_logbit_p
1529 unsigned long int iindex
;
1530 iindex
= scm_to_ulong (index
);
1532 if (SCM_I_INUMP (j
))
1534 /* bits above what's in an inum follow the sign bit */
1535 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1536 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1538 else if (SCM_BIGP (j
))
1540 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1541 scm_remember_upto_here_1 (j
);
1542 return scm_from_bool (val
);
1545 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1550 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1552 "Return the integer which is the ones-complement of the integer\n"
1556 "(number->string (lognot #b10000000) 2)\n"
1557 " @result{} \"-10000001\"\n"
1558 "(number->string (lognot #b0) 2)\n"
1559 " @result{} \"-1\"\n"
1561 #define FUNC_NAME s_scm_lognot
1563 if (SCM_I_INUMP (n
)) {
1564 /* No overflow here, just need to toggle all the bits making up the inum.
1565 Enhancement: No need to strip the tag and add it back, could just xor
1566 a block of 1 bits, if that worked with the various debug versions of
1568 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1570 } else if (SCM_BIGP (n
)) {
1571 SCM result
= scm_i_mkbig ();
1572 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1573 scm_remember_upto_here_1 (n
);
1577 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1582 /* returns 0 if IN is not an integer. OUT must already be
1585 coerce_to_big (SCM in
, mpz_t out
)
1588 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1589 else if (SCM_I_INUMP (in
))
1590 mpz_set_si (out
, SCM_I_INUM (in
));
1597 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1598 (SCM n
, SCM k
, SCM m
),
1599 "Return @var{n} raised to the integer exponent\n"
1600 "@var{k}, modulo @var{m}.\n"
1603 "(modulo-expt 2 3 5)\n"
1606 #define FUNC_NAME s_scm_modulo_expt
1612 /* There are two classes of error we might encounter --
1613 1) Math errors, which we'll report by calling scm_num_overflow,
1615 2) wrong-type errors, which of course we'll report by calling
1617 We don't report those errors immediately, however; instead we do
1618 some cleanup first. These variables tell us which error (if
1619 any) we should report after cleaning up.
1621 int report_overflow
= 0;
1623 int position_of_wrong_type
= 0;
1624 SCM value_of_wrong_type
= SCM_INUM0
;
1626 SCM result
= SCM_UNDEFINED
;
1632 if (scm_is_eq (m
, SCM_INUM0
))
1634 report_overflow
= 1;
1638 if (!coerce_to_big (n
, n_tmp
))
1640 value_of_wrong_type
= n
;
1641 position_of_wrong_type
= 1;
1645 if (!coerce_to_big (k
, k_tmp
))
1647 value_of_wrong_type
= k
;
1648 position_of_wrong_type
= 2;
1652 if (!coerce_to_big (m
, m_tmp
))
1654 value_of_wrong_type
= m
;
1655 position_of_wrong_type
= 3;
1659 /* if the exponent K is negative, and we simply call mpz_powm, we
1660 will get a divide-by-zero exception when an inverse 1/n mod m
1661 doesn't exist (or is not unique). Since exceptions are hard to
1662 handle, we'll attempt the inversion "by hand" -- that way, we get
1663 a simple failure code, which is easy to handle. */
1665 if (-1 == mpz_sgn (k_tmp
))
1667 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1669 report_overflow
= 1;
1672 mpz_neg (k_tmp
, k_tmp
);
1675 result
= scm_i_mkbig ();
1676 mpz_powm (SCM_I_BIG_MPZ (result
),
1681 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1682 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1689 if (report_overflow
)
1690 scm_num_overflow (FUNC_NAME
);
1692 if (position_of_wrong_type
)
1693 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1694 value_of_wrong_type
);
1696 return scm_i_normbig (result
);
1700 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1702 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1703 "exact integer, @var{n} can be any number.\n"
1705 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1706 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1707 "includes @math{0^0} is 1.\n"
1710 "(integer-expt 2 5) @result{} 32\n"
1711 "(integer-expt -3 3) @result{} -27\n"
1712 "(integer-expt 5 -3) @result{} 1/125\n"
1713 "(integer-expt 0 0) @result{} 1\n"
1715 #define FUNC_NAME s_scm_integer_expt
1718 SCM z_i2
= SCM_BOOL_F
;
1720 SCM acc
= SCM_I_MAKINUM (1L);
1722 /* 0^0 == 1 according to R5RS */
1723 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1724 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1725 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1726 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1728 if (SCM_I_INUMP (k
))
1729 i2
= SCM_I_INUM (k
);
1730 else if (SCM_BIGP (k
))
1732 z_i2
= scm_i_clonebig (k
, 1);
1733 scm_remember_upto_here_1 (k
);
1737 SCM_WRONG_TYPE_ARG (2, k
);
1741 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1743 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1744 n
= scm_divide (n
, SCM_UNDEFINED
);
1748 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1752 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1754 return scm_product (acc
, n
);
1756 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1757 acc
= scm_product (acc
, n
);
1758 n
= scm_product (n
, n
);
1759 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1767 n
= scm_divide (n
, SCM_UNDEFINED
);
1774 return scm_product (acc
, n
);
1776 acc
= scm_product (acc
, n
);
1777 n
= scm_product (n
, n
);
1784 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1786 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1787 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1789 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1790 "@var{cnt} is negative it's a division, rounded towards negative\n"
1791 "infinity. (Note that this is not the same rounding as\n"
1792 "@code{quotient} does.)\n"
1794 "With @var{n} viewed as an infinite precision twos complement,\n"
1795 "@code{ash} means a left shift introducing zero bits, or a right\n"
1796 "shift dropping bits.\n"
1799 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1800 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1802 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1803 "(ash -23 -2) @result{} -6\n"
1805 #define FUNC_NAME s_scm_ash
1808 bits_to_shift
= scm_to_long (cnt
);
1810 if (SCM_I_INUMP (n
))
1812 long nn
= SCM_I_INUM (n
);
1814 if (bits_to_shift
> 0)
1816 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1817 overflow a non-zero fixnum. For smaller shifts we check the
1818 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1819 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1820 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1826 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1828 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1831 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1835 SCM result
= scm_i_long2big (nn
);
1836 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1843 bits_to_shift
= -bits_to_shift
;
1844 if (bits_to_shift
>= SCM_LONG_BIT
)
1845 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1847 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1851 else if (SCM_BIGP (n
))
1855 if (bits_to_shift
== 0)
1858 result
= scm_i_mkbig ();
1859 if (bits_to_shift
>= 0)
1861 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1867 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1868 we have to allocate a bignum even if the result is going to be a
1870 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1872 return scm_i_normbig (result
);
1878 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1884 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1885 (SCM n
, SCM start
, SCM end
),
1886 "Return the integer composed of the @var{start} (inclusive)\n"
1887 "through @var{end} (exclusive) bits of @var{n}. The\n"
1888 "@var{start}th bit becomes the 0-th bit in the result.\n"
1891 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1892 " @result{} \"1010\"\n"
1893 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1894 " @result{} \"10110\"\n"
1896 #define FUNC_NAME s_scm_bit_extract
1898 unsigned long int istart
, iend
, bits
;
1899 istart
= scm_to_ulong (start
);
1900 iend
= scm_to_ulong (end
);
1901 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1903 /* how many bits to keep */
1904 bits
= iend
- istart
;
1906 if (SCM_I_INUMP (n
))
1908 long int in
= SCM_I_INUM (n
);
1910 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1911 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1912 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1914 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1916 /* Since we emulate two's complement encoded numbers, this
1917 * special case requires us to produce a result that has
1918 * more bits than can be stored in a fixnum.
1920 SCM result
= scm_i_long2big (in
);
1921 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1926 /* mask down to requisite bits */
1927 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1928 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1930 else if (SCM_BIGP (n
))
1935 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1939 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1940 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1941 such bits into a ulong. */
1942 result
= scm_i_mkbig ();
1943 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1944 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1945 result
= scm_i_normbig (result
);
1947 scm_remember_upto_here_1 (n
);
1951 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1956 static const char scm_logtab
[] = {
1957 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1960 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1962 "Return the number of bits in integer @var{n}. If integer is\n"
1963 "positive, the 1-bits in its binary representation are counted.\n"
1964 "If negative, the 0-bits in its two's-complement binary\n"
1965 "representation are counted. If 0, 0 is returned.\n"
1968 "(logcount #b10101010)\n"
1975 #define FUNC_NAME s_scm_logcount
1977 if (SCM_I_INUMP (n
))
1979 unsigned long int c
= 0;
1980 long int nn
= SCM_I_INUM (n
);
1985 c
+= scm_logtab
[15 & nn
];
1988 return SCM_I_MAKINUM (c
);
1990 else if (SCM_BIGP (n
))
1992 unsigned long count
;
1993 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1994 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1996 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1997 scm_remember_upto_here_1 (n
);
1998 return SCM_I_MAKINUM (count
);
2001 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2006 static const char scm_ilentab
[] = {
2007 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2011 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2013 "Return the number of bits necessary to represent @var{n}.\n"
2016 "(integer-length #b10101010)\n"
2018 "(integer-length 0)\n"
2020 "(integer-length #b1111)\n"
2023 #define FUNC_NAME s_scm_integer_length
2025 if (SCM_I_INUMP (n
))
2027 unsigned long int c
= 0;
2029 long int nn
= SCM_I_INUM (n
);
2035 l
= scm_ilentab
[15 & nn
];
2038 return SCM_I_MAKINUM (c
- 4 + l
);
2040 else if (SCM_BIGP (n
))
2042 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2043 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2044 1 too big, so check for that and adjust. */
2045 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2046 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2047 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2048 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2050 scm_remember_upto_here_1 (n
);
2051 return SCM_I_MAKINUM (size
);
2054 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2058 /*** NUMBERS -> STRINGS ***/
2059 #define SCM_MAX_DBL_PREC 60
2060 #define SCM_MAX_DBL_RADIX 36
2062 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2063 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2064 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2067 void init_dblprec(int *prec
, int radix
) {
2068 /* determine floating point precision by adding successively
2069 smaller increments to 1.0 until it is considered == 1.0 */
2070 double f
= ((double)1.0)/radix
;
2071 double fsum
= 1.0 + f
;
2076 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2088 void init_fx_radix(double *fx_list
, int radix
)
2090 /* initialize a per-radix list of tolerances. When added
2091 to a number < 1.0, we can determine if we should raund
2092 up and quit converting a number to a string. */
2096 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2097 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2100 /* use this array as a way to generate a single digit */
2101 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2104 idbl2str (double f
, char *a
, int radix
)
2106 int efmt
, dpt
, d
, i
, wp
;
2108 #ifdef DBL_MIN_10_EXP
2111 #endif /* DBL_MIN_10_EXP */
2116 radix
> SCM_MAX_DBL_RADIX
)
2118 /* revert to existing behavior */
2122 wp
= scm_dblprec
[radix
-2];
2123 fx
= fx_per_radix
[radix
-2];
2127 #ifdef HAVE_COPYSIGN
2128 double sgn
= copysign (1.0, f
);
2133 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2139 strcpy (a
, "-inf.0");
2141 strcpy (a
, "+inf.0");
2144 else if (xisnan (f
))
2146 strcpy (a
, "+nan.0");
2156 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2157 make-uniform-vector, from causing infinite loops. */
2158 /* just do the checking...if it passes, we do the conversion for our
2159 radix again below */
2166 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2174 while (f_cpy
> 10.0)
2177 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2198 if (f
+ fx
[wp
] >= radix
)
2205 /* adding 9999 makes this equivalent to abs(x) % 3 */
2206 dpt
= (exp
+ 9999) % 3;
2210 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2232 a
[ch
++] = number_chars
[d
];
2235 if (f
+ fx
[wp
] >= 1.0)
2237 a
[ch
- 1] = number_chars
[d
+1];
2249 if ((dpt
> 4) && (exp
> 6))
2251 d
= (a
[0] == '-' ? 2 : 1);
2252 for (i
= ch
++; i
> d
; i
--)
2265 if (a
[ch
- 1] == '.')
2266 a
[ch
++] = '0'; /* trailing zero */
2275 for (i
= radix
; i
<= exp
; i
*= radix
);
2276 for (i
/= radix
; i
; i
/= radix
)
2278 a
[ch
++] = number_chars
[exp
/ i
];
2287 icmplx2str (double real
, double imag
, char *str
, int radix
)
2291 i
= idbl2str (real
, str
, radix
);
2294 /* Don't output a '+' for negative numbers or for Inf and
2295 NaN. They will provide their own sign. */
2296 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2298 i
+= idbl2str (imag
, &str
[i
], radix
);
2305 iflo2str (SCM flt
, char *str
, int radix
)
2308 if (SCM_REALP (flt
))
2309 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2311 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2316 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2317 characters in the result.
2319 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2321 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2326 return scm_iuint2str (-num
, rad
, p
) + 1;
2329 return scm_iuint2str (num
, rad
, p
);
2332 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2333 characters in the result.
2335 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2337 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2341 scm_t_uintmax n
= num
;
2343 for (n
/= rad
; n
> 0; n
/= rad
)
2353 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2358 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2360 "Return a string holding the external representation of the\n"
2361 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2362 "inexact, a radix of 10 will be used.")
2363 #define FUNC_NAME s_scm_number_to_string
2367 if (SCM_UNBNDP (radix
))
2370 base
= scm_to_signed_integer (radix
, 2, 36);
2372 if (SCM_I_INUMP (n
))
2374 char num_buf
[SCM_INTBUFLEN
];
2375 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2376 return scm_from_locale_stringn (num_buf
, length
);
2378 else if (SCM_BIGP (n
))
2380 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2381 scm_remember_upto_here_1 (n
);
2382 return scm_take_locale_string (str
);
2384 else if (SCM_FRACTIONP (n
))
2386 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2387 scm_from_locale_string ("/"),
2388 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2390 else if (SCM_INEXACTP (n
))
2392 char num_buf
[FLOBUFLEN
];
2393 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2396 SCM_WRONG_TYPE_ARG (1, n
);
2401 /* These print routines used to be stubbed here so that scm_repl.c
2402 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2405 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2407 char num_buf
[FLOBUFLEN
];
2408 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2413 scm_i_print_double (double val
, SCM port
)
2415 char num_buf
[FLOBUFLEN
];
2416 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2420 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2423 char num_buf
[FLOBUFLEN
];
2424 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2429 scm_i_print_complex (double real
, double imag
, SCM port
)
2431 char num_buf
[FLOBUFLEN
];
2432 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2436 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2439 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2440 scm_lfwrite (scm_i_string_chars (str
), scm_i_string_length (str
), port
);
2441 scm_remember_upto_here_1 (str
);
2446 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2448 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2449 scm_remember_upto_here_1 (exp
);
2450 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2454 /*** END nums->strs ***/
2457 /*** STRINGS -> NUMBERS ***/
2459 /* The following functions implement the conversion from strings to numbers.
2460 * The implementation somehow follows the grammar for numbers as it is given
2461 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2462 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2463 * points should be noted about the implementation:
2464 * * Each function keeps a local index variable 'idx' that points at the
2465 * current position within the parsed string. The global index is only
2466 * updated if the function could parse the corresponding syntactic unit
2468 * * Similarly, the functions keep track of indicators of inexactness ('#',
2469 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2470 * global exactness information is only updated after each part has been
2471 * successfully parsed.
2472 * * Sequences of digits are parsed into temporary variables holding fixnums.
2473 * Only if these fixnums would overflow, the result variables are updated
2474 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2475 * the temporary variables holding the fixnums are cleared, and the process
2476 * starts over again. If for example fixnums were able to store five decimal
2477 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2478 * and the result was computed as 12345 * 100000 + 67890. In other words,
2479 * only every five digits two bignum operations were performed.
2482 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2484 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2486 /* In non ASCII-style encodings the following macro might not work. */
2487 #define XDIGIT2UINT(d) \
2488 (isdigit ((int) (unsigned char) d) \
2490 : tolower ((int) (unsigned char) d) - 'a' + 10)
2493 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2494 unsigned int radix
, enum t_exactness
*p_exactness
)
2496 unsigned int idx
= *p_idx
;
2497 unsigned int hash_seen
= 0;
2498 scm_t_bits shift
= 1;
2500 unsigned int digit_value
;
2508 if (!isxdigit ((int) (unsigned char) c
))
2510 digit_value
= XDIGIT2UINT (c
);
2511 if (digit_value
>= radix
)
2515 result
= SCM_I_MAKINUM (digit_value
);
2519 if (isxdigit ((int) (unsigned char) c
))
2523 digit_value
= XDIGIT2UINT (c
);
2524 if (digit_value
>= radix
)
2536 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2538 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2540 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2547 shift
= shift
* radix
;
2548 add
= add
* radix
+ digit_value
;
2553 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2555 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2559 *p_exactness
= INEXACT
;
2565 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2566 * covers the parts of the rules that start at a potential point. The value
2567 * of the digits up to the point have been parsed by the caller and are given
2568 * in variable result. The content of *p_exactness indicates, whether a hash
2569 * has already been seen in the digits before the point.
2572 /* In non ASCII-style encodings the following macro might not work. */
2573 #define DIGIT2UINT(d) ((d) - '0')
2576 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2577 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2579 unsigned int idx
= *p_idx
;
2580 enum t_exactness x
= *p_exactness
;
2585 if (mem
[idx
] == '.')
2587 scm_t_bits shift
= 1;
2589 unsigned int digit_value
;
2590 SCM big_shift
= SCM_I_MAKINUM (1);
2596 if (isdigit ((int) (unsigned char) c
))
2601 digit_value
= DIGIT2UINT (c
);
2612 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2614 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2615 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2617 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2625 add
= add
* 10 + digit_value
;
2631 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2632 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2633 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2636 result
= scm_divide (result
, big_shift
);
2638 /* We've seen a decimal point, thus the value is implicitly inexact. */
2650 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2677 if (!isdigit ((int) (unsigned char) c
))
2681 exponent
= DIGIT2UINT (c
);
2685 if (isdigit ((int) (unsigned char) c
))
2688 if (exponent
<= SCM_MAXEXP
)
2689 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2695 if (exponent
> SCM_MAXEXP
)
2697 size_t exp_len
= idx
- start
;
2698 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2699 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2700 scm_out_of_range ("string->number", exp_num
);
2703 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2705 result
= scm_product (result
, e
);
2707 result
= scm_divide2real (result
, e
);
2709 /* We've seen an exponent, thus the value is implicitly inexact. */
2727 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2730 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2731 unsigned int radix
, enum t_exactness
*p_exactness
)
2733 unsigned int idx
= *p_idx
;
2739 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2745 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2747 enum t_exactness x
= EXACT
;
2749 /* Cobble up the fractional part. We might want to set the
2750 NaN's mantissa from it. */
2752 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2757 if (mem
[idx
] == '.')
2761 else if (idx
+ 1 == len
)
2763 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2766 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2767 p_idx
, p_exactness
);
2771 enum t_exactness x
= EXACT
;
2774 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2775 if (scm_is_false (uinteger
))
2780 else if (mem
[idx
] == '/')
2786 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2787 if (scm_is_false (divisor
))
2790 /* both are int/big here, I assume */
2791 result
= scm_i_make_ratio (uinteger
, divisor
);
2793 else if (radix
== 10)
2795 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2796 if (scm_is_false (result
))
2807 /* When returning an inexact zero, make sure it is represented as a
2808 floating point value so that we can change its sign.
2810 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2811 result
= scm_from_double (0.0);
2817 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2820 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2821 unsigned int radix
, enum t_exactness
*p_exactness
)
2845 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2846 if (scm_is_false (ureal
))
2848 /* input must be either +i or -i */
2853 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2859 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2866 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2867 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2876 /* either +<ureal>i or -<ureal>i */
2883 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2886 /* polar input: <real>@<real>. */
2911 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2912 if (scm_is_false (angle
))
2917 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2918 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2920 result
= scm_make_polar (ureal
, angle
);
2925 /* expecting input matching <real>[+-]<ureal>?i */
2932 int sign
= (c
== '+') ? 1 : -1;
2933 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2935 if (scm_is_false (imag
))
2936 imag
= SCM_I_MAKINUM (sign
);
2937 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2938 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2942 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2949 return scm_make_rectangular (ureal
, imag
);
2958 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2960 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2963 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
2964 unsigned int default_radix
)
2966 unsigned int idx
= 0;
2967 unsigned int radix
= NO_RADIX
;
2968 enum t_exactness forced_x
= NO_EXACTNESS
;
2969 enum t_exactness implicit_x
= EXACT
;
2972 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2973 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2975 switch (mem
[idx
+ 1])
2978 if (radix
!= NO_RADIX
)
2983 if (radix
!= NO_RADIX
)
2988 if (forced_x
!= NO_EXACTNESS
)
2993 if (forced_x
!= NO_EXACTNESS
)
2998 if (radix
!= NO_RADIX
)
3003 if (radix
!= NO_RADIX
)
3013 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3014 if (radix
== NO_RADIX
)
3015 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
3017 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
3019 if (scm_is_false (result
))
3025 if (SCM_INEXACTP (result
))
3026 return scm_inexact_to_exact (result
);
3030 if (SCM_INEXACTP (result
))
3033 return scm_exact_to_inexact (result
);
3036 if (implicit_x
== INEXACT
)
3038 if (SCM_INEXACTP (result
))
3041 return scm_exact_to_inexact (result
);
3049 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3050 (SCM string
, SCM radix
),
3051 "Return a number of the maximally precise representation\n"
3052 "expressed by the given @var{string}. @var{radix} must be an\n"
3053 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3054 "is a default radix that may be overridden by an explicit radix\n"
3055 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3056 "supplied, then the default radix is 10. If string is not a\n"
3057 "syntactically valid notation for a number, then\n"
3058 "@code{string->number} returns @code{#f}.")
3059 #define FUNC_NAME s_scm_string_to_number
3063 SCM_VALIDATE_STRING (1, string
);
3065 if (SCM_UNBNDP (radix
))
3068 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3070 answer
= scm_c_locale_stringn_to_number (scm_i_string_chars (string
),
3071 scm_i_string_length (string
),
3073 scm_remember_upto_here_1 (string
);
3079 /*** END strs->nums ***/
3083 scm_bigequal (SCM x
, SCM y
)
3085 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3086 scm_remember_upto_here_2 (x
, y
);
3087 return scm_from_bool (0 == result
);
3091 scm_real_equalp (SCM x
, SCM y
)
3093 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3097 scm_complex_equalp (SCM x
, SCM y
)
3099 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3100 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3104 scm_i_fraction_equalp (SCM x
, SCM y
)
3106 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3107 SCM_FRACTION_NUMERATOR (y
)))
3108 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3109 SCM_FRACTION_DENOMINATOR (y
))))
3116 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3118 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3120 #define FUNC_NAME s_scm_number_p
3122 return scm_from_bool (SCM_NUMBERP (x
));
3126 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3128 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3129 "otherwise. Note that the sets of real, rational and integer\n"
3130 "values form subsets of the set of complex numbers, i. e. the\n"
3131 "predicate will also be fulfilled if @var{x} is a real,\n"
3132 "rational or integer number.")
3133 #define FUNC_NAME s_scm_complex_p
3135 /* all numbers are complex. */
3136 return scm_number_p (x
);
3140 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3142 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3143 "otherwise. Note that the set of integer values forms a subset of\n"
3144 "the set of real numbers, i. e. the predicate will also be\n"
3145 "fulfilled if @var{x} is an integer number.")
3146 #define FUNC_NAME s_scm_real_p
3148 /* we can't represent irrational numbers. */
3149 return scm_rational_p (x
);
3153 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3155 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3156 "otherwise. Note that the set of integer values forms a subset of\n"
3157 "the set of rational numbers, i. e. the predicate will also be\n"
3158 "fulfilled if @var{x} is an integer number.")
3159 #define FUNC_NAME s_scm_rational_p
3161 if (SCM_I_INUMP (x
))
3163 else if (SCM_IMP (x
))
3165 else if (SCM_BIGP (x
))
3167 else if (SCM_FRACTIONP (x
))
3169 else if (SCM_REALP (x
))
3170 /* due to their limited precision, all floating point numbers are
3171 rational as well. */
3178 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3180 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3182 #define FUNC_NAME s_scm_integer_p
3185 if (SCM_I_INUMP (x
))
3191 if (!SCM_INEXACTP (x
))
3193 if (SCM_COMPLEXP (x
))
3195 r
= SCM_REAL_VALUE (x
);
3196 /* +/-inf passes r==floor(r), making those #t */
3204 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3206 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3208 #define FUNC_NAME s_scm_inexact_p
3210 if (SCM_INEXACTP (x
))
3212 if (SCM_NUMBERP (x
))
3214 SCM_WRONG_TYPE_ARG (1, x
);
3219 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3220 /* "Return @code{#t} if all parameters are numerically equal." */
3222 scm_num_eq_p (SCM x
, SCM y
)
3225 if (SCM_I_INUMP (x
))
3227 long xx
= SCM_I_INUM (x
);
3228 if (SCM_I_INUMP (y
))
3230 long yy
= SCM_I_INUM (y
);
3231 return scm_from_bool (xx
== yy
);
3233 else if (SCM_BIGP (y
))
3235 else if (SCM_REALP (y
))
3237 /* On a 32-bit system an inum fits a double, we can cast the inum
3238 to a double and compare.
3240 But on a 64-bit system an inum is bigger than a double and
3241 casting it to a double (call that dxx) will round. dxx is at
3242 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3243 an integer and fits a long. So we cast yy to a long and
3244 compare with plain xx.
3246 An alternative (for any size system actually) would be to check
3247 yy is an integer (with floor) and is in range of an inum
3248 (compare against appropriate powers of 2) then test
3249 xx==(long)yy. It's just a matter of which casts/comparisons
3250 might be fastest or easiest for the cpu. */
3252 double yy
= SCM_REAL_VALUE (y
);
3253 return scm_from_bool ((double) xx
== yy
3254 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3255 || xx
== (long) yy
));
3257 else if (SCM_COMPLEXP (y
))
3258 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3259 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3260 else if (SCM_FRACTIONP (y
))
3263 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3265 else if (SCM_BIGP (x
))
3267 if (SCM_I_INUMP (y
))
3269 else if (SCM_BIGP (y
))
3271 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3272 scm_remember_upto_here_2 (x
, y
);
3273 return scm_from_bool (0 == cmp
);
3275 else if (SCM_REALP (y
))
3278 if (xisnan (SCM_REAL_VALUE (y
)))
3280 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3281 scm_remember_upto_here_1 (x
);
3282 return scm_from_bool (0 == cmp
);
3284 else if (SCM_COMPLEXP (y
))
3287 if (0.0 != SCM_COMPLEX_IMAG (y
))
3289 if (xisnan (SCM_COMPLEX_REAL (y
)))
3291 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3292 scm_remember_upto_here_1 (x
);
3293 return scm_from_bool (0 == cmp
);
3295 else if (SCM_FRACTIONP (y
))
3298 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3300 else if (SCM_REALP (x
))
3302 double xx
= SCM_REAL_VALUE (x
);
3303 if (SCM_I_INUMP (y
))
3305 /* see comments with inum/real above */
3306 long yy
= SCM_I_INUM (y
);
3307 return scm_from_bool (xx
== (double) yy
3308 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3309 || (long) xx
== yy
));
3311 else if (SCM_BIGP (y
))
3314 if (xisnan (SCM_REAL_VALUE (x
)))
3316 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3317 scm_remember_upto_here_1 (y
);
3318 return scm_from_bool (0 == cmp
);
3320 else if (SCM_REALP (y
))
3321 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3322 else if (SCM_COMPLEXP (y
))
3323 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3324 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3325 else if (SCM_FRACTIONP (y
))
3327 double xx
= SCM_REAL_VALUE (x
);
3331 return scm_from_bool (xx
< 0.0);
3332 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3336 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3338 else if (SCM_COMPLEXP (x
))
3340 if (SCM_I_INUMP (y
))
3341 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3342 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3343 else if (SCM_BIGP (y
))
3346 if (0.0 != SCM_COMPLEX_IMAG (x
))
3348 if (xisnan (SCM_COMPLEX_REAL (x
)))
3350 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3351 scm_remember_upto_here_1 (y
);
3352 return scm_from_bool (0 == cmp
);
3354 else if (SCM_REALP (y
))
3355 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3356 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3357 else if (SCM_COMPLEXP (y
))
3358 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3359 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3360 else if (SCM_FRACTIONP (y
))
3363 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3365 xx
= SCM_COMPLEX_REAL (x
);
3369 return scm_from_bool (xx
< 0.0);
3370 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3374 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3376 else if (SCM_FRACTIONP (x
))
3378 if (SCM_I_INUMP (y
))
3380 else if (SCM_BIGP (y
))
3382 else if (SCM_REALP (y
))
3384 double yy
= SCM_REAL_VALUE (y
);
3388 return scm_from_bool (0.0 < yy
);
3389 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3392 else if (SCM_COMPLEXP (y
))
3395 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3397 yy
= SCM_COMPLEX_REAL (y
);
3401 return scm_from_bool (0.0 < yy
);
3402 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3405 else if (SCM_FRACTIONP (y
))
3406 return scm_i_fraction_equalp (x
, y
);
3408 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3411 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3415 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3416 done are good for inums, but for bignums an answer can almost always be
3417 had by just examining a few high bits of the operands, as done by GMP in
3418 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3419 of the float exponent to take into account. */
3421 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3422 /* "Return @code{#t} if the list of parameters is monotonically\n"
3426 scm_less_p (SCM x
, SCM y
)
3429 if (SCM_I_INUMP (x
))
3431 long xx
= SCM_I_INUM (x
);
3432 if (SCM_I_INUMP (y
))
3434 long yy
= SCM_I_INUM (y
);
3435 return scm_from_bool (xx
< yy
);
3437 else if (SCM_BIGP (y
))
3439 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3440 scm_remember_upto_here_1 (y
);
3441 return scm_from_bool (sgn
> 0);
3443 else if (SCM_REALP (y
))
3444 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3445 else if (SCM_FRACTIONP (y
))
3447 /* "x < a/b" becomes "x*b < a" */
3449 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3450 y
= SCM_FRACTION_NUMERATOR (y
);
3454 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3456 else if (SCM_BIGP (x
))
3458 if (SCM_I_INUMP (y
))
3460 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3461 scm_remember_upto_here_1 (x
);
3462 return scm_from_bool (sgn
< 0);
3464 else if (SCM_BIGP (y
))
3466 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3467 scm_remember_upto_here_2 (x
, y
);
3468 return scm_from_bool (cmp
< 0);
3470 else if (SCM_REALP (y
))
3473 if (xisnan (SCM_REAL_VALUE (y
)))
3475 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3476 scm_remember_upto_here_1 (x
);
3477 return scm_from_bool (cmp
< 0);
3479 else if (SCM_FRACTIONP (y
))
3482 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3484 else if (SCM_REALP (x
))
3486 if (SCM_I_INUMP (y
))
3487 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3488 else if (SCM_BIGP (y
))
3491 if (xisnan (SCM_REAL_VALUE (x
)))
3493 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3494 scm_remember_upto_here_1 (y
);
3495 return scm_from_bool (cmp
> 0);
3497 else if (SCM_REALP (y
))
3498 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3499 else if (SCM_FRACTIONP (y
))
3501 double xx
= SCM_REAL_VALUE (x
);
3505 return scm_from_bool (xx
< 0.0);
3506 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3510 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3512 else if (SCM_FRACTIONP (x
))
3514 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3516 /* "a/b < y" becomes "a < y*b" */
3517 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3518 x
= SCM_FRACTION_NUMERATOR (x
);
3521 else if (SCM_REALP (y
))
3523 double yy
= SCM_REAL_VALUE (y
);
3527 return scm_from_bool (0.0 < yy
);
3528 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3531 else if (SCM_FRACTIONP (y
))
3533 /* "a/b < c/d" becomes "a*d < c*b" */
3534 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3535 SCM_FRACTION_DENOMINATOR (y
));
3536 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3537 SCM_FRACTION_DENOMINATOR (x
));
3543 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3546 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3550 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3551 /* "Return @code{#t} if the list of parameters is monotonically\n"
3554 #define FUNC_NAME s_scm_gr_p
3556 scm_gr_p (SCM x
, SCM y
)
3558 if (!SCM_NUMBERP (x
))
3559 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3560 else if (!SCM_NUMBERP (y
))
3561 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3563 return scm_less_p (y
, x
);
3568 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3569 /* "Return @code{#t} if the list of parameters is monotonically\n"
3572 #define FUNC_NAME s_scm_leq_p
3574 scm_leq_p (SCM x
, SCM y
)
3576 if (!SCM_NUMBERP (x
))
3577 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3578 else if (!SCM_NUMBERP (y
))
3579 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3580 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3583 return scm_not (scm_less_p (y
, x
));
3588 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3589 /* "Return @code{#t} if the list of parameters is monotonically\n"
3592 #define FUNC_NAME s_scm_geq_p
3594 scm_geq_p (SCM x
, SCM y
)
3596 if (!SCM_NUMBERP (x
))
3597 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3598 else if (!SCM_NUMBERP (y
))
3599 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3600 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3603 return scm_not (scm_less_p (x
, y
));
3608 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3609 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3615 if (SCM_I_INUMP (z
))
3616 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3617 else if (SCM_BIGP (z
))
3619 else if (SCM_REALP (z
))
3620 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3621 else if (SCM_COMPLEXP (z
))
3622 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3623 && SCM_COMPLEX_IMAG (z
) == 0.0);
3624 else if (SCM_FRACTIONP (z
))
3627 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3631 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3632 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3636 scm_positive_p (SCM x
)
3638 if (SCM_I_INUMP (x
))
3639 return scm_from_bool (SCM_I_INUM (x
) > 0);
3640 else if (SCM_BIGP (x
))
3642 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3643 scm_remember_upto_here_1 (x
);
3644 return scm_from_bool (sgn
> 0);
3646 else if (SCM_REALP (x
))
3647 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3648 else if (SCM_FRACTIONP (x
))
3649 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3651 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3655 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3656 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3660 scm_negative_p (SCM x
)
3662 if (SCM_I_INUMP (x
))
3663 return scm_from_bool (SCM_I_INUM (x
) < 0);
3664 else if (SCM_BIGP (x
))
3666 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3667 scm_remember_upto_here_1 (x
);
3668 return scm_from_bool (sgn
< 0);
3670 else if (SCM_REALP (x
))
3671 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3672 else if (SCM_FRACTIONP (x
))
3673 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3675 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3679 /* scm_min and scm_max return an inexact when either argument is inexact, as
3680 required by r5rs. On that basis, for exact/inexact combinations the
3681 exact is converted to inexact to compare and possibly return. This is
3682 unlike scm_less_p above which takes some trouble to preserve all bits in
3683 its test, such trouble is not required for min and max. */
3685 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3686 /* "Return the maximum of all parameter values."
3689 scm_max (SCM x
, SCM y
)
3694 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3695 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3698 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3701 if (SCM_I_INUMP (x
))
3703 long xx
= SCM_I_INUM (x
);
3704 if (SCM_I_INUMP (y
))
3706 long yy
= SCM_I_INUM (y
);
3707 return (xx
< yy
) ? y
: x
;
3709 else if (SCM_BIGP (y
))
3711 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3712 scm_remember_upto_here_1 (y
);
3713 return (sgn
< 0) ? x
: y
;
3715 else if (SCM_REALP (y
))
3718 /* if y==NaN then ">" is false and we return NaN */
3719 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3721 else if (SCM_FRACTIONP (y
))
3724 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3727 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3729 else if (SCM_BIGP (x
))
3731 if (SCM_I_INUMP (y
))
3733 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3734 scm_remember_upto_here_1 (x
);
3735 return (sgn
< 0) ? y
: x
;
3737 else if (SCM_BIGP (y
))
3739 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3740 scm_remember_upto_here_2 (x
, y
);
3741 return (cmp
> 0) ? x
: y
;
3743 else if (SCM_REALP (y
))
3745 /* if y==NaN then xx>yy is false, so we return the NaN y */
3748 xx
= scm_i_big2dbl (x
);
3749 yy
= SCM_REAL_VALUE (y
);
3750 return (xx
> yy
? scm_from_double (xx
) : y
);
3752 else if (SCM_FRACTIONP (y
))
3757 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3759 else if (SCM_REALP (x
))
3761 if (SCM_I_INUMP (y
))
3763 double z
= SCM_I_INUM (y
);
3764 /* if x==NaN then "<" is false and we return NaN */
3765 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3767 else if (SCM_BIGP (y
))
3772 else if (SCM_REALP (y
))
3774 /* if x==NaN then our explicit check means we return NaN
3775 if y==NaN then ">" is false and we return NaN
3776 calling isnan is unavoidable, since it's the only way to know
3777 which of x or y causes any compares to be false */
3778 double xx
= SCM_REAL_VALUE (x
);
3779 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3781 else if (SCM_FRACTIONP (y
))
3783 double yy
= scm_i_fraction2double (y
);
3784 double xx
= SCM_REAL_VALUE (x
);
3785 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3788 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3790 else if (SCM_FRACTIONP (x
))
3792 if (SCM_I_INUMP (y
))
3796 else if (SCM_BIGP (y
))
3800 else if (SCM_REALP (y
))
3802 double xx
= scm_i_fraction2double (x
);
3803 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3805 else if (SCM_FRACTIONP (y
))
3810 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3813 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3817 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3818 /* "Return the minium of all parameter values."
3821 scm_min (SCM x
, SCM y
)
3826 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3827 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3830 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3833 if (SCM_I_INUMP (x
))
3835 long xx
= SCM_I_INUM (x
);
3836 if (SCM_I_INUMP (y
))
3838 long yy
= SCM_I_INUM (y
);
3839 return (xx
< yy
) ? x
: y
;
3841 else if (SCM_BIGP (y
))
3843 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3844 scm_remember_upto_here_1 (y
);
3845 return (sgn
< 0) ? y
: x
;
3847 else if (SCM_REALP (y
))
3850 /* if y==NaN then "<" is false and we return NaN */
3851 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3853 else if (SCM_FRACTIONP (y
))
3856 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3859 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3861 else if (SCM_BIGP (x
))
3863 if (SCM_I_INUMP (y
))
3865 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3866 scm_remember_upto_here_1 (x
);
3867 return (sgn
< 0) ? x
: y
;
3869 else if (SCM_BIGP (y
))
3871 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3872 scm_remember_upto_here_2 (x
, y
);
3873 return (cmp
> 0) ? y
: x
;
3875 else if (SCM_REALP (y
))
3877 /* if y==NaN then xx<yy is false, so we return the NaN y */
3880 xx
= scm_i_big2dbl (x
);
3881 yy
= SCM_REAL_VALUE (y
);
3882 return (xx
< yy
? scm_from_double (xx
) : y
);
3884 else if (SCM_FRACTIONP (y
))
3889 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3891 else if (SCM_REALP (x
))
3893 if (SCM_I_INUMP (y
))
3895 double z
= SCM_I_INUM (y
);
3896 /* if x==NaN then "<" is false and we return NaN */
3897 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3899 else if (SCM_BIGP (y
))
3904 else if (SCM_REALP (y
))
3906 /* if x==NaN then our explicit check means we return NaN
3907 if y==NaN then "<" is false and we return NaN
3908 calling isnan is unavoidable, since it's the only way to know
3909 which of x or y causes any compares to be false */
3910 double xx
= SCM_REAL_VALUE (x
);
3911 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3913 else if (SCM_FRACTIONP (y
))
3915 double yy
= scm_i_fraction2double (y
);
3916 double xx
= SCM_REAL_VALUE (x
);
3917 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3920 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3922 else if (SCM_FRACTIONP (x
))
3924 if (SCM_I_INUMP (y
))
3928 else if (SCM_BIGP (y
))
3932 else if (SCM_REALP (y
))
3934 double xx
= scm_i_fraction2double (x
);
3935 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3937 else if (SCM_FRACTIONP (y
))
3942 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3945 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3949 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3950 /* "Return the sum of all parameter values. Return 0 if called without\n"
3954 scm_sum (SCM x
, SCM y
)
3956 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
3958 if (SCM_NUMBERP (x
)) return x
;
3959 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3960 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3963 if (SCM_LIKELY (SCM_I_INUMP (x
)))
3965 if (SCM_LIKELY (SCM_I_INUMP (y
)))
3967 long xx
= SCM_I_INUM (x
);
3968 long yy
= SCM_I_INUM (y
);
3969 long int z
= xx
+ yy
;
3970 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3972 else if (SCM_BIGP (y
))
3977 else if (SCM_REALP (y
))
3979 long int xx
= SCM_I_INUM (x
);
3980 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3982 else if (SCM_COMPLEXP (y
))
3984 long int xx
= SCM_I_INUM (x
);
3985 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3986 SCM_COMPLEX_IMAG (y
));
3988 else if (SCM_FRACTIONP (y
))
3989 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3990 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3991 SCM_FRACTION_DENOMINATOR (y
));
3993 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3994 } else if (SCM_BIGP (x
))
3996 if (SCM_I_INUMP (y
))
4001 inum
= SCM_I_INUM (y
);
4004 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4007 SCM result
= scm_i_mkbig ();
4008 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4009 scm_remember_upto_here_1 (x
);
4010 /* we know the result will have to be a bignum */
4013 return scm_i_normbig (result
);
4017 SCM result
= scm_i_mkbig ();
4018 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4019 scm_remember_upto_here_1 (x
);
4020 /* we know the result will have to be a bignum */
4023 return scm_i_normbig (result
);
4026 else if (SCM_BIGP (y
))
4028 SCM result
= scm_i_mkbig ();
4029 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4030 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4031 mpz_add (SCM_I_BIG_MPZ (result
),
4034 scm_remember_upto_here_2 (x
, y
);
4035 /* we know the result will have to be a bignum */
4038 return scm_i_normbig (result
);
4040 else if (SCM_REALP (y
))
4042 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4043 scm_remember_upto_here_1 (x
);
4044 return scm_from_double (result
);
4046 else if (SCM_COMPLEXP (y
))
4048 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4049 + SCM_COMPLEX_REAL (y
));
4050 scm_remember_upto_here_1 (x
);
4051 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4053 else if (SCM_FRACTIONP (y
))
4054 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4055 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4056 SCM_FRACTION_DENOMINATOR (y
));
4058 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4060 else if (SCM_REALP (x
))
4062 if (SCM_I_INUMP (y
))
4063 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4064 else if (SCM_BIGP (y
))
4066 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4067 scm_remember_upto_here_1 (y
);
4068 return scm_from_double (result
);
4070 else if (SCM_REALP (y
))
4071 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4072 else if (SCM_COMPLEXP (y
))
4073 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4074 SCM_COMPLEX_IMAG (y
));
4075 else if (SCM_FRACTIONP (y
))
4076 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4078 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4080 else if (SCM_COMPLEXP (x
))
4082 if (SCM_I_INUMP (y
))
4083 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4084 SCM_COMPLEX_IMAG (x
));
4085 else if (SCM_BIGP (y
))
4087 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4088 + SCM_COMPLEX_REAL (x
));
4089 scm_remember_upto_here_1 (y
);
4090 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4092 else if (SCM_REALP (y
))
4093 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4094 SCM_COMPLEX_IMAG (x
));
4095 else if (SCM_COMPLEXP (y
))
4096 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4097 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4098 else if (SCM_FRACTIONP (y
))
4099 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4100 SCM_COMPLEX_IMAG (x
));
4102 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4104 else if (SCM_FRACTIONP (x
))
4106 if (SCM_I_INUMP (y
))
4107 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4108 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4109 SCM_FRACTION_DENOMINATOR (x
));
4110 else if (SCM_BIGP (y
))
4111 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4112 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4113 SCM_FRACTION_DENOMINATOR (x
));
4114 else if (SCM_REALP (y
))
4115 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4116 else if (SCM_COMPLEXP (y
))
4117 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4118 SCM_COMPLEX_IMAG (y
));
4119 else if (SCM_FRACTIONP (y
))
4120 /* a/b + c/d = (ad + bc) / bd */
4121 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4122 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4123 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4125 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4128 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4132 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4134 "Return @math{@var{x}+1}.")
4135 #define FUNC_NAME s_scm_oneplus
4137 return scm_sum (x
, SCM_I_MAKINUM (1));
4142 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4143 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4144 * the sum of all but the first argument are subtracted from the first
4146 #define FUNC_NAME s_difference
4148 scm_difference (SCM x
, SCM y
)
4150 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4153 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4155 if (SCM_I_INUMP (x
))
4157 long xx
= -SCM_I_INUM (x
);
4158 if (SCM_FIXABLE (xx
))
4159 return SCM_I_MAKINUM (xx
);
4161 return scm_i_long2big (xx
);
4163 else if (SCM_BIGP (x
))
4164 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4165 bignum, but negating that gives a fixnum. */
4166 return scm_i_normbig (scm_i_clonebig (x
, 0));
4167 else if (SCM_REALP (x
))
4168 return scm_from_double (-SCM_REAL_VALUE (x
));
4169 else if (SCM_COMPLEXP (x
))
4170 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4171 -SCM_COMPLEX_IMAG (x
));
4172 else if (SCM_FRACTIONP (x
))
4173 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4174 SCM_FRACTION_DENOMINATOR (x
));
4176 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4179 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4181 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4183 long int xx
= SCM_I_INUM (x
);
4184 long int yy
= SCM_I_INUM (y
);
4185 long int z
= xx
- yy
;
4186 if (SCM_FIXABLE (z
))
4187 return SCM_I_MAKINUM (z
);
4189 return scm_i_long2big (z
);
4191 else if (SCM_BIGP (y
))
4193 /* inum-x - big-y */
4194 long xx
= SCM_I_INUM (x
);
4197 return scm_i_clonebig (y
, 0);
4200 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4201 SCM result
= scm_i_mkbig ();
4204 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4207 /* x - y == -(y + -x) */
4208 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4209 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4211 scm_remember_upto_here_1 (y
);
4213 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4214 /* we know the result will have to be a bignum */
4217 return scm_i_normbig (result
);
4220 else if (SCM_REALP (y
))
4222 long int xx
= SCM_I_INUM (x
);
4223 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4225 else if (SCM_COMPLEXP (y
))
4227 long int xx
= SCM_I_INUM (x
);
4228 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4229 - SCM_COMPLEX_IMAG (y
));
4231 else if (SCM_FRACTIONP (y
))
4232 /* a - b/c = (ac - b) / c */
4233 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4234 SCM_FRACTION_NUMERATOR (y
)),
4235 SCM_FRACTION_DENOMINATOR (y
));
4237 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4239 else if (SCM_BIGP (x
))
4241 if (SCM_I_INUMP (y
))
4243 /* big-x - inum-y */
4244 long yy
= SCM_I_INUM (y
);
4245 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4247 scm_remember_upto_here_1 (x
);
4249 return (SCM_FIXABLE (-yy
) ?
4250 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4253 SCM result
= scm_i_mkbig ();
4256 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4258 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4259 scm_remember_upto_here_1 (x
);
4261 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4262 /* we know the result will have to be a bignum */
4265 return scm_i_normbig (result
);
4268 else if (SCM_BIGP (y
))
4270 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4271 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4272 SCM result
= scm_i_mkbig ();
4273 mpz_sub (SCM_I_BIG_MPZ (result
),
4276 scm_remember_upto_here_2 (x
, y
);
4277 /* we know the result will have to be a bignum */
4278 if ((sgn_x
== 1) && (sgn_y
== -1))
4280 if ((sgn_x
== -1) && (sgn_y
== 1))
4282 return scm_i_normbig (result
);
4284 else if (SCM_REALP (y
))
4286 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4287 scm_remember_upto_here_1 (x
);
4288 return scm_from_double (result
);
4290 else if (SCM_COMPLEXP (y
))
4292 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4293 - SCM_COMPLEX_REAL (y
));
4294 scm_remember_upto_here_1 (x
);
4295 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4297 else if (SCM_FRACTIONP (y
))
4298 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4299 SCM_FRACTION_NUMERATOR (y
)),
4300 SCM_FRACTION_DENOMINATOR (y
));
4301 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4303 else if (SCM_REALP (x
))
4305 if (SCM_I_INUMP (y
))
4306 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4307 else if (SCM_BIGP (y
))
4309 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4310 scm_remember_upto_here_1 (x
);
4311 return scm_from_double (result
);
4313 else if (SCM_REALP (y
))
4314 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4315 else if (SCM_COMPLEXP (y
))
4316 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4317 -SCM_COMPLEX_IMAG (y
));
4318 else if (SCM_FRACTIONP (y
))
4319 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4321 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4323 else if (SCM_COMPLEXP (x
))
4325 if (SCM_I_INUMP (y
))
4326 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4327 SCM_COMPLEX_IMAG (x
));
4328 else if (SCM_BIGP (y
))
4330 double real_part
= (SCM_COMPLEX_REAL (x
)
4331 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4332 scm_remember_upto_here_1 (x
);
4333 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4335 else if (SCM_REALP (y
))
4336 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4337 SCM_COMPLEX_IMAG (x
));
4338 else if (SCM_COMPLEXP (y
))
4339 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4340 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4341 else if (SCM_FRACTIONP (y
))
4342 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4343 SCM_COMPLEX_IMAG (x
));
4345 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4347 else if (SCM_FRACTIONP (x
))
4349 if (SCM_I_INUMP (y
))
4350 /* a/b - c = (a - cb) / b */
4351 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4352 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4353 SCM_FRACTION_DENOMINATOR (x
));
4354 else if (SCM_BIGP (y
))
4355 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4356 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4357 SCM_FRACTION_DENOMINATOR (x
));
4358 else if (SCM_REALP (y
))
4359 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4360 else if (SCM_COMPLEXP (y
))
4361 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4362 -SCM_COMPLEX_IMAG (y
));
4363 else if (SCM_FRACTIONP (y
))
4364 /* a/b - c/d = (ad - bc) / bd */
4365 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4366 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4367 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4369 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4372 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4377 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4379 "Return @math{@var{x}-1}.")
4380 #define FUNC_NAME s_scm_oneminus
4382 return scm_difference (x
, SCM_I_MAKINUM (1));
4387 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4388 /* "Return the product of all arguments. If called without arguments,\n"
4392 scm_product (SCM x
, SCM y
)
4394 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4397 return SCM_I_MAKINUM (1L);
4398 else if (SCM_NUMBERP (x
))
4401 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4404 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4409 xx
= SCM_I_INUM (x
);
4413 case 0: return x
; break;
4414 case 1: return y
; break;
4417 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4419 long yy
= SCM_I_INUM (y
);
4421 SCM k
= SCM_I_MAKINUM (kk
);
4422 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4426 SCM result
= scm_i_long2big (xx
);
4427 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4428 return scm_i_normbig (result
);
4431 else if (SCM_BIGP (y
))
4433 SCM result
= scm_i_mkbig ();
4434 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4435 scm_remember_upto_here_1 (y
);
4438 else if (SCM_REALP (y
))
4439 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4440 else if (SCM_COMPLEXP (y
))
4441 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4442 xx
* SCM_COMPLEX_IMAG (y
));
4443 else if (SCM_FRACTIONP (y
))
4444 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4445 SCM_FRACTION_DENOMINATOR (y
));
4447 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4449 else if (SCM_BIGP (x
))
4451 if (SCM_I_INUMP (y
))
4456 else if (SCM_BIGP (y
))
4458 SCM result
= scm_i_mkbig ();
4459 mpz_mul (SCM_I_BIG_MPZ (result
),
4462 scm_remember_upto_here_2 (x
, y
);
4465 else if (SCM_REALP (y
))
4467 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4468 scm_remember_upto_here_1 (x
);
4469 return scm_from_double (result
);
4471 else if (SCM_COMPLEXP (y
))
4473 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4474 scm_remember_upto_here_1 (x
);
4475 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4476 z
* SCM_COMPLEX_IMAG (y
));
4478 else if (SCM_FRACTIONP (y
))
4479 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4480 SCM_FRACTION_DENOMINATOR (y
));
4482 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4484 else if (SCM_REALP (x
))
4486 if (SCM_I_INUMP (y
))
4488 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4489 if (scm_is_eq (y
, SCM_INUM0
))
4491 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4493 else if (SCM_BIGP (y
))
4495 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4496 scm_remember_upto_here_1 (y
);
4497 return scm_from_double (result
);
4499 else if (SCM_REALP (y
))
4500 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4501 else if (SCM_COMPLEXP (y
))
4502 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4503 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4504 else if (SCM_FRACTIONP (y
))
4505 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4507 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4509 else if (SCM_COMPLEXP (x
))
4511 if (SCM_I_INUMP (y
))
4513 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4514 if (scm_is_eq (y
, SCM_INUM0
))
4516 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4517 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4519 else if (SCM_BIGP (y
))
4521 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4522 scm_remember_upto_here_1 (y
);
4523 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4524 z
* SCM_COMPLEX_IMAG (x
));
4526 else if (SCM_REALP (y
))
4527 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4528 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4529 else if (SCM_COMPLEXP (y
))
4531 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4532 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4533 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4534 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4536 else if (SCM_FRACTIONP (y
))
4538 double yy
= scm_i_fraction2double (y
);
4539 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4540 yy
* SCM_COMPLEX_IMAG (x
));
4543 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4545 else if (SCM_FRACTIONP (x
))
4547 if (SCM_I_INUMP (y
))
4548 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4549 SCM_FRACTION_DENOMINATOR (x
));
4550 else if (SCM_BIGP (y
))
4551 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4552 SCM_FRACTION_DENOMINATOR (x
));
4553 else if (SCM_REALP (y
))
4554 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4555 else if (SCM_COMPLEXP (y
))
4557 double xx
= scm_i_fraction2double (x
);
4558 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4559 xx
* SCM_COMPLEX_IMAG (y
));
4561 else if (SCM_FRACTIONP (y
))
4562 /* a/b * c/d = ac / bd */
4563 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4564 SCM_FRACTION_NUMERATOR (y
)),
4565 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4566 SCM_FRACTION_DENOMINATOR (y
)));
4568 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4571 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4574 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4575 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4576 #define ALLOW_DIVIDE_BY_ZERO
4577 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4580 /* The code below for complex division is adapted from the GNU
4581 libstdc++, which adapted it from f2c's libF77, and is subject to
4584 /****************************************************************
4585 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4587 Permission to use, copy, modify, and distribute this software
4588 and its documentation for any purpose and without fee is hereby
4589 granted, provided that the above copyright notice appear in all
4590 copies and that both that the copyright notice and this
4591 permission notice and warranty disclaimer appear in supporting
4592 documentation, and that the names of AT&T Bell Laboratories or
4593 Bellcore or any of their entities not be used in advertising or
4594 publicity pertaining to distribution of the software without
4595 specific, written prior permission.
4597 AT&T and Bellcore disclaim all warranties with regard to this
4598 software, including all implied warranties of merchantability
4599 and fitness. In no event shall AT&T or Bellcore be liable for
4600 any special, indirect or consequential damages or any damages
4601 whatsoever resulting from loss of use, data or profits, whether
4602 in an action of contract, negligence or other tortious action,
4603 arising out of or in connection with the use or performance of
4605 ****************************************************************/
4607 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4608 /* Divide the first argument by the product of the remaining
4609 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4611 #define FUNC_NAME s_divide
4613 scm_i_divide (SCM x
, SCM y
, int inexact
)
4617 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4620 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4621 else if (SCM_I_INUMP (x
))
4623 long xx
= SCM_I_INUM (x
);
4624 if (xx
== 1 || xx
== -1)
4626 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4628 scm_num_overflow (s_divide
);
4633 return scm_from_double (1.0 / (double) xx
);
4634 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4637 else if (SCM_BIGP (x
))
4640 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4641 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4643 else if (SCM_REALP (x
))
4645 double xx
= SCM_REAL_VALUE (x
);
4646 #ifndef ALLOW_DIVIDE_BY_ZERO
4648 scm_num_overflow (s_divide
);
4651 return scm_from_double (1.0 / xx
);
4653 else if (SCM_COMPLEXP (x
))
4655 double r
= SCM_COMPLEX_REAL (x
);
4656 double i
= SCM_COMPLEX_IMAG (x
);
4657 if (fabs(r
) <= fabs(i
))
4660 double d
= i
* (1.0 + t
* t
);
4661 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4666 double d
= r
* (1.0 + t
* t
);
4667 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4670 else if (SCM_FRACTIONP (x
))
4671 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4672 SCM_FRACTION_NUMERATOR (x
));
4674 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4677 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4679 long xx
= SCM_I_INUM (x
);
4680 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4682 long yy
= SCM_I_INUM (y
);
4685 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4686 scm_num_overflow (s_divide
);
4688 return scm_from_double ((double) xx
/ (double) yy
);
4691 else if (xx
% yy
!= 0)
4694 return scm_from_double ((double) xx
/ (double) yy
);
4695 else return scm_i_make_ratio (x
, y
);
4700 if (SCM_FIXABLE (z
))
4701 return SCM_I_MAKINUM (z
);
4703 return scm_i_long2big (z
);
4706 else if (SCM_BIGP (y
))
4709 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4710 else return scm_i_make_ratio (x
, y
);
4712 else if (SCM_REALP (y
))
4714 double yy
= SCM_REAL_VALUE (y
);
4715 #ifndef ALLOW_DIVIDE_BY_ZERO
4717 scm_num_overflow (s_divide
);
4720 return scm_from_double ((double) xx
/ yy
);
4722 else if (SCM_COMPLEXP (y
))
4725 complex_div
: /* y _must_ be a complex number */
4727 double r
= SCM_COMPLEX_REAL (y
);
4728 double i
= SCM_COMPLEX_IMAG (y
);
4729 if (fabs(r
) <= fabs(i
))
4732 double d
= i
* (1.0 + t
* t
);
4733 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4738 double d
= r
* (1.0 + t
* t
);
4739 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4743 else if (SCM_FRACTIONP (y
))
4744 /* a / b/c = ac / b */
4745 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4746 SCM_FRACTION_NUMERATOR (y
));
4748 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4750 else if (SCM_BIGP (x
))
4752 if (SCM_I_INUMP (y
))
4754 long int yy
= SCM_I_INUM (y
);
4757 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4758 scm_num_overflow (s_divide
);
4760 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4761 scm_remember_upto_here_1 (x
);
4762 return (sgn
== 0) ? scm_nan () : scm_inf ();
4769 /* FIXME: HMM, what are the relative performance issues here?
4770 We need to test. Is it faster on average to test
4771 divisible_p, then perform whichever operation, or is it
4772 faster to perform the integer div opportunistically and
4773 switch to real if there's a remainder? For now we take the
4774 middle ground: test, then if divisible, use the faster div
4777 long abs_yy
= yy
< 0 ? -yy
: yy
;
4778 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4782 SCM result
= scm_i_mkbig ();
4783 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4784 scm_remember_upto_here_1 (x
);
4786 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4787 return scm_i_normbig (result
);
4792 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4793 else return scm_i_make_ratio (x
, y
);
4797 else if (SCM_BIGP (y
))
4799 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4802 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4803 scm_num_overflow (s_divide
);
4805 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4806 scm_remember_upto_here_1 (x
);
4807 return (sgn
== 0) ? scm_nan () : scm_inf ();
4815 /* It's easily possible for the ratio x/y to fit a double
4816 but one or both x and y be too big to fit a double,
4817 hence the use of mpq_get_d rather than converting and
4820 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
4821 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
4822 return scm_from_double (mpq_get_d (q
));
4826 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4830 SCM result
= scm_i_mkbig ();
4831 mpz_divexact (SCM_I_BIG_MPZ (result
),
4834 scm_remember_upto_here_2 (x
, y
);
4835 return scm_i_normbig (result
);
4838 return scm_i_make_ratio (x
, y
);
4842 else if (SCM_REALP (y
))
4844 double yy
= SCM_REAL_VALUE (y
);
4845 #ifndef ALLOW_DIVIDE_BY_ZERO
4847 scm_num_overflow (s_divide
);
4850 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4852 else if (SCM_COMPLEXP (y
))
4854 a
= scm_i_big2dbl (x
);
4857 else if (SCM_FRACTIONP (y
))
4858 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4859 SCM_FRACTION_NUMERATOR (y
));
4861 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4863 else if (SCM_REALP (x
))
4865 double rx
= SCM_REAL_VALUE (x
);
4866 if (SCM_I_INUMP (y
))
4868 long int yy
= SCM_I_INUM (y
);
4869 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4871 scm_num_overflow (s_divide
);
4874 return scm_from_double (rx
/ (double) yy
);
4876 else if (SCM_BIGP (y
))
4878 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4879 scm_remember_upto_here_1 (y
);
4880 return scm_from_double (rx
/ dby
);
4882 else if (SCM_REALP (y
))
4884 double yy
= SCM_REAL_VALUE (y
);
4885 #ifndef ALLOW_DIVIDE_BY_ZERO
4887 scm_num_overflow (s_divide
);
4890 return scm_from_double (rx
/ yy
);
4892 else if (SCM_COMPLEXP (y
))
4897 else if (SCM_FRACTIONP (y
))
4898 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4900 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4902 else if (SCM_COMPLEXP (x
))
4904 double rx
= SCM_COMPLEX_REAL (x
);
4905 double ix
= SCM_COMPLEX_IMAG (x
);
4906 if (SCM_I_INUMP (y
))
4908 long int yy
= SCM_I_INUM (y
);
4909 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4911 scm_num_overflow (s_divide
);
4916 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4919 else if (SCM_BIGP (y
))
4921 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4922 scm_remember_upto_here_1 (y
);
4923 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4925 else if (SCM_REALP (y
))
4927 double yy
= SCM_REAL_VALUE (y
);
4928 #ifndef ALLOW_DIVIDE_BY_ZERO
4930 scm_num_overflow (s_divide
);
4933 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4935 else if (SCM_COMPLEXP (y
))
4937 double ry
= SCM_COMPLEX_REAL (y
);
4938 double iy
= SCM_COMPLEX_IMAG (y
);
4939 if (fabs(ry
) <= fabs(iy
))
4942 double d
= iy
* (1.0 + t
* t
);
4943 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4948 double d
= ry
* (1.0 + t
* t
);
4949 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4952 else if (SCM_FRACTIONP (y
))
4954 double yy
= scm_i_fraction2double (y
);
4955 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4958 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4960 else if (SCM_FRACTIONP (x
))
4962 if (SCM_I_INUMP (y
))
4964 long int yy
= SCM_I_INUM (y
);
4965 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4967 scm_num_overflow (s_divide
);
4970 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4971 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4973 else if (SCM_BIGP (y
))
4975 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4976 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4978 else if (SCM_REALP (y
))
4980 double yy
= SCM_REAL_VALUE (y
);
4981 #ifndef ALLOW_DIVIDE_BY_ZERO
4983 scm_num_overflow (s_divide
);
4986 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4988 else if (SCM_COMPLEXP (y
))
4990 a
= scm_i_fraction2double (x
);
4993 else if (SCM_FRACTIONP (y
))
4994 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4995 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4997 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5000 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5004 scm_divide (SCM x
, SCM y
)
5006 return scm_i_divide (x
, y
, 0);
5009 static SCM
scm_divide2real (SCM x
, SCM y
)
5011 return scm_i_divide (x
, y
, 1);
5017 scm_asinh (double x
)
5022 #define asinh scm_asinh
5023 return log (x
+ sqrt (x
* x
+ 1));
5026 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
5027 /* "Return the inverse hyperbolic sine of @var{x}."
5032 scm_acosh (double x
)
5037 #define acosh scm_acosh
5038 return log (x
+ sqrt (x
* x
- 1));
5041 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
5042 /* "Return the inverse hyperbolic cosine of @var{x}."
5047 scm_atanh (double x
)
5052 #define atanh scm_atanh
5053 return 0.5 * log ((1 + x
) / (1 - x
));
5056 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
5057 /* "Return the inverse hyperbolic tangent of @var{x}."
5062 scm_c_truncate (double x
)
5073 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5074 half-way case (ie. when x is an integer plus 0.5) going upwards.
5075 Then half-way cases are identified and adjusted down if the
5076 round-upwards didn't give the desired even integer.
5078 "plus_half == result" identifies a half-way case. If plus_half, which is
5079 x + 0.5, is an integer then x must be an integer plus 0.5.
5081 An odd "result" value is identified with result/2 != floor(result/2).
5082 This is done with plus_half, since that value is ready for use sooner in
5083 a pipelined cpu, and we're already requiring plus_half == result.
5085 Note however that we need to be careful when x is big and already an
5086 integer. In that case "x+0.5" may round to an adjacent integer, causing
5087 us to return such a value, incorrectly. For instance if the hardware is
5088 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5089 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5090 returned. Or if the hardware is in round-upwards mode, then other bigger
5091 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5092 representable value, 2^128+2^76 (or whatever), again incorrect.
5094 These bad roundings of x+0.5 are avoided by testing at the start whether
5095 x is already an integer. If it is then clearly that's the desired result
5096 already. And if it's not then the exponent must be small enough to allow
5097 an 0.5 to be represented, and hence added without a bad rounding. */
5100 scm_c_round (double x
)
5102 double plus_half
, result
;
5107 plus_half
= x
+ 0.5;
5108 result
= floor (plus_half
);
5109 /* Adjust so that the rounding is towards even. */
5110 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5115 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5117 "Round the number @var{x} towards zero.")
5118 #define FUNC_NAME s_scm_truncate_number
5120 if (scm_is_false (scm_negative_p (x
)))
5121 return scm_floor (x
);
5123 return scm_ceiling (x
);
5127 static SCM exactly_one_half
;
5129 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5131 "Round the number @var{x} towards the nearest integer. "
5132 "When it is exactly halfway between two integers, "
5133 "round towards the even one.")
5134 #define FUNC_NAME s_scm_round_number
5136 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5138 else if (SCM_REALP (x
))
5139 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5142 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5143 single quotient+remainder division then examining to see which way
5144 the rounding should go. */
5145 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5146 SCM result
= scm_floor (plus_half
);
5147 /* Adjust so that the rounding is towards even. */
5148 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5149 && scm_is_true (scm_odd_p (result
)))
5150 return scm_difference (result
, SCM_I_MAKINUM (1));
5157 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5159 "Round the number @var{x} towards minus infinity.")
5160 #define FUNC_NAME s_scm_floor
5162 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5164 else if (SCM_REALP (x
))
5165 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5166 else if (SCM_FRACTIONP (x
))
5168 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5169 SCM_FRACTION_DENOMINATOR (x
));
5170 if (scm_is_false (scm_negative_p (x
)))
5172 /* For positive x, rounding towards zero is correct. */
5177 /* For negative x, we need to return q-1 unless x is an
5178 integer. But fractions are never integer, per our
5180 return scm_difference (q
, SCM_I_MAKINUM (1));
5184 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5188 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5190 "Round the number @var{x} towards infinity.")
5191 #define FUNC_NAME s_scm_ceiling
5193 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5195 else if (SCM_REALP (x
))
5196 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5197 else if (SCM_FRACTIONP (x
))
5199 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5200 SCM_FRACTION_DENOMINATOR (x
));
5201 if (scm_is_false (scm_positive_p (x
)))
5203 /* For negative x, rounding towards zero is correct. */
5208 /* For positive x, we need to return q+1 unless x is an
5209 integer. But fractions are never integer, per our
5211 return scm_sum (q
, SCM_I_MAKINUM (1));
5215 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5219 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5220 /* "Return the square root of the real number @var{x}."
5222 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5223 /* "Return the absolute value of the real number @var{x}."
5225 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5226 /* "Return the @var{x}th power of e."
5228 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5229 /* "Return the natural logarithm of the real number @var{x}."
5231 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5232 /* "Return the sine of the real number @var{x}."
5234 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5235 /* "Return the cosine of the real number @var{x}."
5237 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5238 /* "Return the tangent of the real number @var{x}."
5240 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5241 /* "Return the arc sine of the real number @var{x}."
5243 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5244 /* "Return the arc cosine of the real number @var{x}."
5246 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5247 /* "Return the arc tangent of the real number @var{x}."
5249 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5250 /* "Return the hyperbolic sine of the real number @var{x}."
5252 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5253 /* "Return the hyperbolic cosine of the real number @var{x}."
5255 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5256 /* "Return the hyperbolic tangent of the real number @var{x}."
5264 static void scm_two_doubles (SCM x
,
5266 const char *sstring
,
5270 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5272 if (SCM_I_INUMP (x
))
5273 xy
->x
= SCM_I_INUM (x
);
5274 else if (SCM_BIGP (x
))
5275 xy
->x
= scm_i_big2dbl (x
);
5276 else if (SCM_REALP (x
))
5277 xy
->x
= SCM_REAL_VALUE (x
);
5278 else if (SCM_FRACTIONP (x
))
5279 xy
->x
= scm_i_fraction2double (x
);
5281 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5283 if (SCM_I_INUMP (y
))
5284 xy
->y
= SCM_I_INUM (y
);
5285 else if (SCM_BIGP (y
))
5286 xy
->y
= scm_i_big2dbl (y
);
5287 else if (SCM_REALP (y
))
5288 xy
->y
= SCM_REAL_VALUE (y
);
5289 else if (SCM_FRACTIONP (y
))
5290 xy
->y
= scm_i_fraction2double (y
);
5292 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5296 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5298 "Return @var{x} raised to the power of @var{y}. This\n"
5299 "procedure does not accept complex arguments.")
5300 #define FUNC_NAME s_scm_sys_expt
5303 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5304 return scm_from_double (pow (xy
.x
, xy
.y
));
5309 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5311 "Return the arc tangent of the two arguments @var{x} and\n"
5312 "@var{y}. This is similar to calculating the arc tangent of\n"
5313 "@var{x} / @var{y}, except that the signs of both arguments\n"
5314 "are used to determine the quadrant of the result. This\n"
5315 "procedure does not accept complex arguments.")
5316 #define FUNC_NAME s_scm_sys_atan2
5319 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5320 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5325 scm_c_make_rectangular (double re
, double im
)
5328 return scm_from_double (re
);
5332 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5334 SCM_COMPLEX_REAL (z
) = re
;
5335 SCM_COMPLEX_IMAG (z
) = im
;
5340 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5341 (SCM real_part
, SCM imaginary_part
),
5342 "Return a complex number constructed of the given @var{real-part} "
5343 "and @var{imaginary-part} parts.")
5344 #define FUNC_NAME s_scm_make_rectangular
5347 scm_two_doubles (real_part
, imaginary_part
, FUNC_NAME
, &xy
);
5348 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5353 scm_c_make_polar (double mag
, double ang
)
5357 sincos (ang
, &s
, &c
);
5362 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5365 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5367 "Return the complex number @var{x} * e^(i * @var{y}).")
5368 #define FUNC_NAME s_scm_make_polar
5371 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5372 return scm_c_make_polar (xy
.x
, xy
.y
);
5377 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5378 /* "Return the real part of the number @var{z}."
5381 scm_real_part (SCM z
)
5383 if (SCM_I_INUMP (z
))
5385 else if (SCM_BIGP (z
))
5387 else if (SCM_REALP (z
))
5389 else if (SCM_COMPLEXP (z
))
5390 return scm_from_double (SCM_COMPLEX_REAL (z
));
5391 else if (SCM_FRACTIONP (z
))
5394 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5398 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5399 /* "Return the imaginary part of the number @var{z}."
5402 scm_imag_part (SCM z
)
5404 if (SCM_I_INUMP (z
))
5406 else if (SCM_BIGP (z
))
5408 else if (SCM_REALP (z
))
5410 else if (SCM_COMPLEXP (z
))
5411 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5412 else if (SCM_FRACTIONP (z
))
5415 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5418 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5419 /* "Return the numerator of the number @var{z}."
5422 scm_numerator (SCM z
)
5424 if (SCM_I_INUMP (z
))
5426 else if (SCM_BIGP (z
))
5428 else if (SCM_FRACTIONP (z
))
5429 return SCM_FRACTION_NUMERATOR (z
);
5430 else if (SCM_REALP (z
))
5431 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5433 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5437 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5438 /* "Return the denominator of the number @var{z}."
5441 scm_denominator (SCM z
)
5443 if (SCM_I_INUMP (z
))
5444 return SCM_I_MAKINUM (1);
5445 else if (SCM_BIGP (z
))
5446 return SCM_I_MAKINUM (1);
5447 else if (SCM_FRACTIONP (z
))
5448 return SCM_FRACTION_DENOMINATOR (z
);
5449 else if (SCM_REALP (z
))
5450 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5452 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5455 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5456 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5457 * "@code{abs} for real arguments, but also allows complex numbers."
5460 scm_magnitude (SCM z
)
5462 if (SCM_I_INUMP (z
))
5464 long int zz
= SCM_I_INUM (z
);
5467 else if (SCM_POSFIXABLE (-zz
))
5468 return SCM_I_MAKINUM (-zz
);
5470 return scm_i_long2big (-zz
);
5472 else if (SCM_BIGP (z
))
5474 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5475 scm_remember_upto_here_1 (z
);
5477 return scm_i_clonebig (z
, 0);
5481 else if (SCM_REALP (z
))
5482 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5483 else if (SCM_COMPLEXP (z
))
5484 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5485 else if (SCM_FRACTIONP (z
))
5487 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5489 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5490 SCM_FRACTION_DENOMINATOR (z
));
5493 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5497 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5498 /* "Return the angle of the complex number @var{z}."
5503 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5504 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5505 But if atan2 follows the floating point rounding mode, then the value
5506 is not a constant. Maybe it'd be close enough though. */
5507 if (SCM_I_INUMP (z
))
5509 if (SCM_I_INUM (z
) >= 0)
5512 return scm_from_double (atan2 (0.0, -1.0));
5514 else if (SCM_BIGP (z
))
5516 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5517 scm_remember_upto_here_1 (z
);
5519 return scm_from_double (atan2 (0.0, -1.0));
5523 else if (SCM_REALP (z
))
5525 if (SCM_REAL_VALUE (z
) >= 0)
5528 return scm_from_double (atan2 (0.0, -1.0));
5530 else if (SCM_COMPLEXP (z
))
5531 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5532 else if (SCM_FRACTIONP (z
))
5534 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5536 else return scm_from_double (atan2 (0.0, -1.0));
5539 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5543 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5544 /* Convert the number @var{x} to its inexact representation.\n"
5547 scm_exact_to_inexact (SCM z
)
5549 if (SCM_I_INUMP (z
))
5550 return scm_from_double ((double) SCM_I_INUM (z
));
5551 else if (SCM_BIGP (z
))
5552 return scm_from_double (scm_i_big2dbl (z
));
5553 else if (SCM_FRACTIONP (z
))
5554 return scm_from_double (scm_i_fraction2double (z
));
5555 else if (SCM_INEXACTP (z
))
5558 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5562 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5564 "Return an exact number that is numerically closest to @var{z}.")
5565 #define FUNC_NAME s_scm_inexact_to_exact
5567 if (SCM_I_INUMP (z
))
5569 else if (SCM_BIGP (z
))
5571 else if (SCM_REALP (z
))
5573 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5574 SCM_OUT_OF_RANGE (1, z
);
5581 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5582 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5583 scm_i_mpz2num (mpq_denref (frac
)));
5585 /* When scm_i_make_ratio throws, we leak the memory allocated
5592 else if (SCM_FRACTIONP (z
))
5595 SCM_WRONG_TYPE_ARG (1, z
);
5599 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5601 "Return an exact number that is within @var{err} of @var{x}.")
5602 #define FUNC_NAME s_scm_rationalize
5604 if (SCM_I_INUMP (x
))
5606 else if (SCM_BIGP (x
))
5608 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5610 /* Use continued fractions to find closest ratio. All
5611 arithmetic is done with exact numbers.
5614 SCM ex
= scm_inexact_to_exact (x
);
5615 SCM int_part
= scm_floor (ex
);
5616 SCM tt
= SCM_I_MAKINUM (1);
5617 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5618 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5622 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5625 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5626 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5628 /* We stop after a million iterations just to be absolutely sure
5629 that we don't go into an infinite loop. The process normally
5630 converges after less than a dozen iterations.
5633 err
= scm_abs (err
);
5634 while (++i
< 1000000)
5636 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5637 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5638 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5640 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5641 err
))) /* abs(x-a/b) <= err */
5643 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5644 if (scm_is_false (scm_exact_p (x
))
5645 || scm_is_false (scm_exact_p (err
)))
5646 return scm_exact_to_inexact (res
);
5650 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5652 tt
= scm_floor (rx
); /* tt = floor (rx) */
5658 scm_num_overflow (s_scm_rationalize
);
5661 SCM_WRONG_TYPE_ARG (1, x
);
5665 /* conversion functions */
5668 scm_is_integer (SCM val
)
5670 return scm_is_true (scm_integer_p (val
));
5674 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5676 if (SCM_I_INUMP (val
))
5678 scm_t_signed_bits n
= SCM_I_INUM (val
);
5679 return n
>= min
&& n
<= max
;
5681 else if (SCM_BIGP (val
))
5683 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5685 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5687 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5689 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5690 return n
>= min
&& n
<= max
;
5700 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5701 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5704 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5705 SCM_I_BIG_MPZ (val
));
5707 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5719 return n
>= min
&& n
<= max
;
5727 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5729 if (SCM_I_INUMP (val
))
5731 scm_t_signed_bits n
= SCM_I_INUM (val
);
5732 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5734 else if (SCM_BIGP (val
))
5736 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5738 else if (max
<= ULONG_MAX
)
5740 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5742 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5743 return n
>= min
&& n
<= max
;
5753 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5756 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5757 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5760 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5761 SCM_I_BIG_MPZ (val
));
5763 return n
>= min
&& n
<= max
;
5771 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5773 scm_error (scm_out_of_range_key
,
5775 "Value out of range ~S to ~S: ~S",
5776 scm_list_3 (min
, max
, bad_val
),
5777 scm_list_1 (bad_val
));
5780 #define TYPE scm_t_intmax
5781 #define TYPE_MIN min
5782 #define TYPE_MAX max
5783 #define SIZEOF_TYPE 0
5784 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5785 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5786 #include "libguile/conv-integer.i.c"
5788 #define TYPE scm_t_uintmax
5789 #define TYPE_MIN min
5790 #define TYPE_MAX max
5791 #define SIZEOF_TYPE 0
5792 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5793 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5794 #include "libguile/conv-uinteger.i.c"
5796 #define TYPE scm_t_int8
5797 #define TYPE_MIN SCM_T_INT8_MIN
5798 #define TYPE_MAX SCM_T_INT8_MAX
5799 #define SIZEOF_TYPE 1
5800 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5801 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5802 #include "libguile/conv-integer.i.c"
5804 #define TYPE scm_t_uint8
5806 #define TYPE_MAX SCM_T_UINT8_MAX
5807 #define SIZEOF_TYPE 1
5808 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5809 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5810 #include "libguile/conv-uinteger.i.c"
5812 #define TYPE scm_t_int16
5813 #define TYPE_MIN SCM_T_INT16_MIN
5814 #define TYPE_MAX SCM_T_INT16_MAX
5815 #define SIZEOF_TYPE 2
5816 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5817 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5818 #include "libguile/conv-integer.i.c"
5820 #define TYPE scm_t_uint16
5822 #define TYPE_MAX SCM_T_UINT16_MAX
5823 #define SIZEOF_TYPE 2
5824 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5825 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5826 #include "libguile/conv-uinteger.i.c"
5828 #define TYPE scm_t_int32
5829 #define TYPE_MIN SCM_T_INT32_MIN
5830 #define TYPE_MAX SCM_T_INT32_MAX
5831 #define SIZEOF_TYPE 4
5832 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5833 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5834 #include "libguile/conv-integer.i.c"
5836 #define TYPE scm_t_uint32
5838 #define TYPE_MAX SCM_T_UINT32_MAX
5839 #define SIZEOF_TYPE 4
5840 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5841 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5842 #include "libguile/conv-uinteger.i.c"
5844 #if SCM_HAVE_T_INT64
5846 #define TYPE scm_t_int64
5847 #define TYPE_MIN SCM_T_INT64_MIN
5848 #define TYPE_MAX SCM_T_INT64_MAX
5849 #define SIZEOF_TYPE 8
5850 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5851 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5852 #include "libguile/conv-integer.i.c"
5854 #define TYPE scm_t_uint64
5856 #define TYPE_MAX SCM_T_UINT64_MAX
5857 #define SIZEOF_TYPE 8
5858 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5859 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5860 #include "libguile/conv-uinteger.i.c"
5865 scm_to_mpz (SCM val
, mpz_t rop
)
5867 if (SCM_I_INUMP (val
))
5868 mpz_set_si (rop
, SCM_I_INUM (val
));
5869 else if (SCM_BIGP (val
))
5870 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5872 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5876 scm_from_mpz (mpz_t val
)
5878 return scm_i_mpz2num (val
);
5882 scm_is_real (SCM val
)
5884 return scm_is_true (scm_real_p (val
));
5888 scm_is_rational (SCM val
)
5890 return scm_is_true (scm_rational_p (val
));
5894 scm_to_double (SCM val
)
5896 if (SCM_I_INUMP (val
))
5897 return SCM_I_INUM (val
);
5898 else if (SCM_BIGP (val
))
5899 return scm_i_big2dbl (val
);
5900 else if (SCM_FRACTIONP (val
))
5901 return scm_i_fraction2double (val
);
5902 else if (SCM_REALP (val
))
5903 return SCM_REAL_VALUE (val
);
5905 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5909 scm_from_double (double val
)
5911 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5912 SCM_REAL_VALUE (z
) = val
;
5916 #if SCM_ENABLE_DISCOURAGED == 1
5919 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5923 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5927 scm_out_of_range (NULL
, num
);
5930 return scm_to_double (num
);
5934 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5938 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5942 scm_out_of_range (NULL
, num
);
5945 return scm_to_double (num
);
5951 scm_is_complex (SCM val
)
5953 return scm_is_true (scm_complex_p (val
));
5957 scm_c_real_part (SCM z
)
5959 if (SCM_COMPLEXP (z
))
5960 return SCM_COMPLEX_REAL (z
);
5963 /* Use the scm_real_part to get proper error checking and
5966 return scm_to_double (scm_real_part (z
));
5971 scm_c_imag_part (SCM z
)
5973 if (SCM_COMPLEXP (z
))
5974 return SCM_COMPLEX_IMAG (z
);
5977 /* Use the scm_imag_part to get proper error checking and
5978 dispatching. The result will almost always be 0.0, but not
5981 return scm_to_double (scm_imag_part (z
));
5986 scm_c_magnitude (SCM z
)
5988 return scm_to_double (scm_magnitude (z
));
5994 return scm_to_double (scm_angle (z
));
5998 scm_is_number (SCM z
)
6000 return scm_is_true (scm_number_p (z
));
6004 /* In the following functions we dispatch to the real-arg funcs like log()
6005 when we know the arg is real, instead of just handing everything to
6006 clog() for instance. This is in case clog() doesn't optimize for a
6007 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6008 well use it to go straight to the applicable C func. */
6010 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6012 "Return the natural logarithm of @var{z}.")
6013 #define FUNC_NAME s_scm_log
6015 if (SCM_COMPLEXP (z
))
6017 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6018 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6020 double re
= SCM_COMPLEX_REAL (z
);
6021 double im
= SCM_COMPLEX_IMAG (z
);
6022 return scm_c_make_rectangular (log (hypot (re
, im
)),
6028 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6029 although the value itself overflows. */
6030 double re
= scm_to_double (z
);
6031 double l
= log (fabs (re
));
6033 return scm_from_double (l
);
6035 return scm_c_make_rectangular (l
, M_PI
);
6041 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6043 "Return the base 10 logarithm of @var{z}.")
6044 #define FUNC_NAME s_scm_log10
6046 if (SCM_COMPLEXP (z
))
6048 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6049 clog() and a multiply by M_LOG10E, rather than the fallback
6050 log10+hypot+atan2.) */
6051 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6052 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6054 double re
= SCM_COMPLEX_REAL (z
);
6055 double im
= SCM_COMPLEX_IMAG (z
);
6056 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6057 M_LOG10E
* atan2 (im
, re
));
6062 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6063 although the value itself overflows. */
6064 double re
= scm_to_double (z
);
6065 double l
= log10 (fabs (re
));
6067 return scm_from_double (l
);
6069 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6075 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6077 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6078 "base of natural logarithms (2.71828@dots{}).")
6079 #define FUNC_NAME s_scm_exp
6081 if (SCM_COMPLEXP (z
))
6083 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6084 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6086 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6087 SCM_COMPLEX_IMAG (z
));
6092 /* When z is a negative bignum the conversion to double overflows,
6093 giving -infinity, but that's ok, the exp is still 0.0. */
6094 return scm_from_double (exp (scm_to_double (z
)));
6100 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6102 "Return the square root of @var{z}. Of the two possible roots\n"
6103 "(positive and negative), the one with the a positive real part\n"
6104 "is returned, or if that's zero then a positive imaginary part.\n"
6108 "(sqrt 9.0) @result{} 3.0\n"
6109 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6110 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6111 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6113 #define FUNC_NAME s_scm_sqrt
6115 if (SCM_COMPLEXP (x
))
6117 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6118 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6120 double re
= SCM_COMPLEX_REAL (x
);
6121 double im
= SCM_COMPLEX_IMAG (x
);
6122 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6123 0.5 * atan2 (im
, re
));
6128 double xx
= scm_to_double (x
);
6130 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6132 return scm_from_double (sqrt (xx
));
6144 mpz_init_set_si (z_negative_one
, -1);
6146 /* It may be possible to tune the performance of some algorithms by using
6147 * the following constants to avoid the creation of bignums. Please, before
6148 * using these values, remember the two rules of program optimization:
6149 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6150 scm_c_define ("most-positive-fixnum",
6151 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6152 scm_c_define ("most-negative-fixnum",
6153 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6155 scm_add_feature ("complex");
6156 scm_add_feature ("inexact");
6157 scm_flo0
= scm_from_double (0.0);
6159 /* determine floating point precision */
6160 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6162 init_dblprec(&scm_dblprec
[i
-2],i
);
6163 init_fx_radix(fx_per_radix
[i
-2],i
);
6166 /* hard code precision for base 10 if the preprocessor tells us to... */
6167 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6170 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6171 SCM_I_MAKINUM (2)));
6172 #include "libguile/numbers.x"