1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008, 2009 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 License
9 * as published by the Free Software Foundation; either version 3 of
10 * the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful, but
13 * 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
24 /* General assumptions:
25 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
26 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
27 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
28 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
29 * All objects satisfying SCM_FRACTIONP are never an integer.
34 - see if special casing bignums and reals in integer-exponent when
35 possible (to use mpz_pow and mpf_pow_ui) is faster.
37 - look in to better short-circuiting of common cases in
38 integer-expt and elsewhere.
40 - see if direct mpz operations can help in ash and elsewhere.
56 #include "libguile/_scm.h"
57 #include "libguile/feature.h"
58 #include "libguile/ports.h"
59 #include "libguile/root.h"
60 #include "libguile/smob.h"
61 #include "libguile/strings.h"
63 #include "libguile/validate.h"
64 #include "libguile/numbers.h"
65 #include "libguile/deprecation.h"
67 #include "libguile/eq.h"
69 #include "libguile/discouraged.h"
71 /* values per glibc, if not already defined */
73 #define M_LOG10E 0.43429448190325182765
76 #define M_PI 3.14159265358979323846
82 Wonder if this might be faster for some of our code? A switch on
83 the numtag would jump directly to the right case, and the
84 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
86 #define SCM_I_NUMTAG_NOTNUM 0
87 #define SCM_I_NUMTAG_INUM 1
88 #define SCM_I_NUMTAG_BIG scm_tc16_big
89 #define SCM_I_NUMTAG_REAL scm_tc16_real
90 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
91 #define SCM_I_NUMTAG(x) \
92 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
93 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
94 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
95 : SCM_I_NUMTAG_NOTNUM)))
97 /* the macro above will not work as is with fractions */
100 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
102 /* FLOBUFLEN is the maximum number of characters neccessary for the
103 * printed or scm_string representation of an inexact number.
105 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
108 #if ! defined (HAVE_ISNAN)
113 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
116 #if ! defined (HAVE_ISINF)
121 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
128 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
129 an explicit check. In some future gmp (don't know what version number),
130 mpz_cmp_d is supposed to do this itself. */
132 #define xmpz_cmp_d(z, d) \
133 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
135 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
138 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
139 isinf. It does have finite and isnan though, hence the use of those.
140 fpclass would be a possibility on that system too. */
144 #if defined (HAVE_ISINF)
146 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
147 return (! (finite (x
) || isnan (x
)));
156 #if defined (HAVE_ISNAN)
163 #if defined (GUILE_I)
164 #if HAVE_COMPLEX_DOUBLE
166 /* For an SCM object Z which is a complex number (ie. satisfies
167 SCM_COMPLEXP), return its value as a C level "complex double". */
168 #define SCM_COMPLEX_VALUE(z) \
169 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
171 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
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
: scm_abs (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
;
2736 /* Start off believing that the number will be exact. This changes
2737 to INEXACT if we see a decimal point or a hash. */
2738 enum t_exactness x
= EXACT
;
2743 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2749 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2751 /* Cobble up the fractional part. We might want to set the
2752 NaN's mantissa from it. */
2754 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2759 if (mem
[idx
] == '.')
2763 else if (idx
+ 1 == len
)
2765 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2768 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2775 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2776 if (scm_is_false (uinteger
))
2781 else if (mem
[idx
] == '/')
2787 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2788 if (scm_is_false (divisor
))
2791 /* both are int/big here, I assume */
2792 result
= scm_i_make_ratio (uinteger
, divisor
);
2794 else if (radix
== 10)
2796 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2797 if (scm_is_false (result
))
2806 /* Update *p_exactness if the number just read was inexact. This is
2807 important for complex numbers, so that a complex number is
2808 treated as inexact overall if either its real or imaginary part
2814 /* When returning an inexact zero, make sure it is represented as a
2815 floating point value so that we can change its sign.
2817 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2818 result
= scm_from_double (0.0);
2824 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2827 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2828 unsigned int radix
, enum t_exactness
*p_exactness
)
2852 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2853 if (scm_is_false (ureal
))
2855 /* input must be either +i or -i */
2860 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2866 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2873 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2874 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2883 /* either +<ureal>i or -<ureal>i */
2890 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2893 /* polar input: <real>@<real>. */
2918 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2919 if (scm_is_false (angle
))
2924 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2925 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2927 result
= scm_make_polar (ureal
, angle
);
2932 /* expecting input matching <real>[+-]<ureal>?i */
2939 int sign
= (c
== '+') ? 1 : -1;
2940 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2942 if (scm_is_false (imag
))
2943 imag
= SCM_I_MAKINUM (sign
);
2944 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2945 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2949 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2956 return scm_make_rectangular (ureal
, imag
);
2965 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2967 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2970 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
2971 unsigned int default_radix
)
2973 unsigned int idx
= 0;
2974 unsigned int radix
= NO_RADIX
;
2975 enum t_exactness forced_x
= NO_EXACTNESS
;
2976 enum t_exactness implicit_x
= EXACT
;
2979 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2980 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2982 switch (mem
[idx
+ 1])
2985 if (radix
!= NO_RADIX
)
2990 if (radix
!= NO_RADIX
)
2995 if (forced_x
!= NO_EXACTNESS
)
3000 if (forced_x
!= NO_EXACTNESS
)
3005 if (radix
!= NO_RADIX
)
3010 if (radix
!= NO_RADIX
)
3020 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3021 if (radix
== NO_RADIX
)
3022 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
3024 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
3026 if (scm_is_false (result
))
3032 if (SCM_INEXACTP (result
))
3033 return scm_inexact_to_exact (result
);
3037 if (SCM_INEXACTP (result
))
3040 return scm_exact_to_inexact (result
);
3043 if (implicit_x
== INEXACT
)
3045 if (SCM_INEXACTP (result
))
3048 return scm_exact_to_inexact (result
);
3056 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3057 (SCM string
, SCM radix
),
3058 "Return a number of the maximally precise representation\n"
3059 "expressed by the given @var{string}. @var{radix} must be an\n"
3060 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3061 "is a default radix that may be overridden by an explicit radix\n"
3062 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3063 "supplied, then the default radix is 10. If string is not a\n"
3064 "syntactically valid notation for a number, then\n"
3065 "@code{string->number} returns @code{#f}.")
3066 #define FUNC_NAME s_scm_string_to_number
3070 SCM_VALIDATE_STRING (1, string
);
3072 if (SCM_UNBNDP (radix
))
3075 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3077 answer
= scm_c_locale_stringn_to_number (scm_i_string_chars (string
),
3078 scm_i_string_length (string
),
3080 scm_remember_upto_here_1 (string
);
3086 /*** END strs->nums ***/
3090 scm_bigequal (SCM x
, SCM y
)
3092 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3093 scm_remember_upto_here_2 (x
, y
);
3094 return scm_from_bool (0 == result
);
3098 scm_real_equalp (SCM x
, SCM y
)
3100 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3104 scm_complex_equalp (SCM x
, SCM y
)
3106 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3107 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3111 scm_i_fraction_equalp (SCM x
, SCM y
)
3113 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3114 SCM_FRACTION_NUMERATOR (y
)))
3115 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3116 SCM_FRACTION_DENOMINATOR (y
))))
3123 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3125 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3127 #define FUNC_NAME s_scm_number_p
3129 return scm_from_bool (SCM_NUMBERP (x
));
3133 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3135 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3136 "otherwise. Note that the sets of real, rational and integer\n"
3137 "values form subsets of the set of complex numbers, i. e. the\n"
3138 "predicate will also be fulfilled if @var{x} is a real,\n"
3139 "rational or integer number.")
3140 #define FUNC_NAME s_scm_complex_p
3142 /* all numbers are complex. */
3143 return scm_number_p (x
);
3147 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3149 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3150 "otherwise. Note that the set of integer values forms a subset of\n"
3151 "the set of real numbers, i. e. the predicate will also be\n"
3152 "fulfilled if @var{x} is an integer number.")
3153 #define FUNC_NAME s_scm_real_p
3155 /* we can't represent irrational numbers. */
3156 return scm_rational_p (x
);
3160 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3162 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3163 "otherwise. Note that the set of integer values forms a subset of\n"
3164 "the set of rational numbers, i. e. the predicate will also be\n"
3165 "fulfilled if @var{x} is an integer number.")
3166 #define FUNC_NAME s_scm_rational_p
3168 if (SCM_I_INUMP (x
))
3170 else if (SCM_IMP (x
))
3172 else if (SCM_BIGP (x
))
3174 else if (SCM_FRACTIONP (x
))
3176 else if (SCM_REALP (x
))
3177 /* due to their limited precision, all floating point numbers are
3178 rational as well. */
3185 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3187 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3189 #define FUNC_NAME s_scm_integer_p
3192 if (SCM_I_INUMP (x
))
3198 if (!SCM_INEXACTP (x
))
3200 if (SCM_COMPLEXP (x
))
3202 r
= SCM_REAL_VALUE (x
);
3203 /* +/-inf passes r==floor(r), making those #t */
3211 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3213 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3215 #define FUNC_NAME s_scm_inexact_p
3217 if (SCM_INEXACTP (x
))
3219 if (SCM_NUMBERP (x
))
3221 SCM_WRONG_TYPE_ARG (1, x
);
3226 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3227 /* "Return @code{#t} if all parameters are numerically equal." */
3229 scm_num_eq_p (SCM x
, SCM y
)
3232 if (SCM_I_INUMP (x
))
3234 long xx
= SCM_I_INUM (x
);
3235 if (SCM_I_INUMP (y
))
3237 long yy
= SCM_I_INUM (y
);
3238 return scm_from_bool (xx
== yy
);
3240 else if (SCM_BIGP (y
))
3242 else if (SCM_REALP (y
))
3244 /* On a 32-bit system an inum fits a double, we can cast the inum
3245 to a double and compare.
3247 But on a 64-bit system an inum is bigger than a double and
3248 casting it to a double (call that dxx) will round. dxx is at
3249 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3250 an integer and fits a long. So we cast yy to a long and
3251 compare with plain xx.
3253 An alternative (for any size system actually) would be to check
3254 yy is an integer (with floor) and is in range of an inum
3255 (compare against appropriate powers of 2) then test
3256 xx==(long)yy. It's just a matter of which casts/comparisons
3257 might be fastest or easiest for the cpu. */
3259 double yy
= SCM_REAL_VALUE (y
);
3260 return scm_from_bool ((double) xx
== yy
3261 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3262 || xx
== (long) yy
));
3264 else if (SCM_COMPLEXP (y
))
3265 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3266 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3267 else if (SCM_FRACTIONP (y
))
3270 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3272 else if (SCM_BIGP (x
))
3274 if (SCM_I_INUMP (y
))
3276 else if (SCM_BIGP (y
))
3278 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3279 scm_remember_upto_here_2 (x
, y
);
3280 return scm_from_bool (0 == cmp
);
3282 else if (SCM_REALP (y
))
3285 if (xisnan (SCM_REAL_VALUE (y
)))
3287 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3288 scm_remember_upto_here_1 (x
);
3289 return scm_from_bool (0 == cmp
);
3291 else if (SCM_COMPLEXP (y
))
3294 if (0.0 != SCM_COMPLEX_IMAG (y
))
3296 if (xisnan (SCM_COMPLEX_REAL (y
)))
3298 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3299 scm_remember_upto_here_1 (x
);
3300 return scm_from_bool (0 == cmp
);
3302 else if (SCM_FRACTIONP (y
))
3305 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3307 else if (SCM_REALP (x
))
3309 double xx
= SCM_REAL_VALUE (x
);
3310 if (SCM_I_INUMP (y
))
3312 /* see comments with inum/real above */
3313 long yy
= SCM_I_INUM (y
);
3314 return scm_from_bool (xx
== (double) yy
3315 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3316 || (long) xx
== yy
));
3318 else if (SCM_BIGP (y
))
3321 if (xisnan (SCM_REAL_VALUE (x
)))
3323 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3324 scm_remember_upto_here_1 (y
);
3325 return scm_from_bool (0 == cmp
);
3327 else if (SCM_REALP (y
))
3328 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3329 else if (SCM_COMPLEXP (y
))
3330 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3331 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3332 else if (SCM_FRACTIONP (y
))
3334 double xx
= SCM_REAL_VALUE (x
);
3338 return scm_from_bool (xx
< 0.0);
3339 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3343 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3345 else if (SCM_COMPLEXP (x
))
3347 if (SCM_I_INUMP (y
))
3348 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3349 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3350 else if (SCM_BIGP (y
))
3353 if (0.0 != SCM_COMPLEX_IMAG (x
))
3355 if (xisnan (SCM_COMPLEX_REAL (x
)))
3357 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3358 scm_remember_upto_here_1 (y
);
3359 return scm_from_bool (0 == cmp
);
3361 else if (SCM_REALP (y
))
3362 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3363 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3364 else if (SCM_COMPLEXP (y
))
3365 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3366 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3367 else if (SCM_FRACTIONP (y
))
3370 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3372 xx
= SCM_COMPLEX_REAL (x
);
3376 return scm_from_bool (xx
< 0.0);
3377 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3381 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3383 else if (SCM_FRACTIONP (x
))
3385 if (SCM_I_INUMP (y
))
3387 else if (SCM_BIGP (y
))
3389 else if (SCM_REALP (y
))
3391 double yy
= SCM_REAL_VALUE (y
);
3395 return scm_from_bool (0.0 < yy
);
3396 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3399 else if (SCM_COMPLEXP (y
))
3402 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3404 yy
= SCM_COMPLEX_REAL (y
);
3408 return scm_from_bool (0.0 < yy
);
3409 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3412 else if (SCM_FRACTIONP (y
))
3413 return scm_i_fraction_equalp (x
, y
);
3415 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3418 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3422 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3423 done are good for inums, but for bignums an answer can almost always be
3424 had by just examining a few high bits of the operands, as done by GMP in
3425 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3426 of the float exponent to take into account. */
3428 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3429 /* "Return @code{#t} if the list of parameters is monotonically\n"
3433 scm_less_p (SCM x
, SCM y
)
3436 if (SCM_I_INUMP (x
))
3438 long xx
= SCM_I_INUM (x
);
3439 if (SCM_I_INUMP (y
))
3441 long yy
= SCM_I_INUM (y
);
3442 return scm_from_bool (xx
< yy
);
3444 else if (SCM_BIGP (y
))
3446 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3447 scm_remember_upto_here_1 (y
);
3448 return scm_from_bool (sgn
> 0);
3450 else if (SCM_REALP (y
))
3451 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3452 else if (SCM_FRACTIONP (y
))
3454 /* "x < a/b" becomes "x*b < a" */
3456 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3457 y
= SCM_FRACTION_NUMERATOR (y
);
3461 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3463 else if (SCM_BIGP (x
))
3465 if (SCM_I_INUMP (y
))
3467 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3468 scm_remember_upto_here_1 (x
);
3469 return scm_from_bool (sgn
< 0);
3471 else if (SCM_BIGP (y
))
3473 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3474 scm_remember_upto_here_2 (x
, y
);
3475 return scm_from_bool (cmp
< 0);
3477 else if (SCM_REALP (y
))
3480 if (xisnan (SCM_REAL_VALUE (y
)))
3482 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3483 scm_remember_upto_here_1 (x
);
3484 return scm_from_bool (cmp
< 0);
3486 else if (SCM_FRACTIONP (y
))
3489 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3491 else if (SCM_REALP (x
))
3493 if (SCM_I_INUMP (y
))
3494 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3495 else if (SCM_BIGP (y
))
3498 if (xisnan (SCM_REAL_VALUE (x
)))
3500 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3501 scm_remember_upto_here_1 (y
);
3502 return scm_from_bool (cmp
> 0);
3504 else if (SCM_REALP (y
))
3505 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3506 else if (SCM_FRACTIONP (y
))
3508 double xx
= SCM_REAL_VALUE (x
);
3512 return scm_from_bool (xx
< 0.0);
3513 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3517 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3519 else if (SCM_FRACTIONP (x
))
3521 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3523 /* "a/b < y" becomes "a < y*b" */
3524 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3525 x
= SCM_FRACTION_NUMERATOR (x
);
3528 else if (SCM_REALP (y
))
3530 double yy
= SCM_REAL_VALUE (y
);
3534 return scm_from_bool (0.0 < yy
);
3535 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3538 else if (SCM_FRACTIONP (y
))
3540 /* "a/b < c/d" becomes "a*d < c*b" */
3541 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3542 SCM_FRACTION_DENOMINATOR (y
));
3543 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3544 SCM_FRACTION_DENOMINATOR (x
));
3550 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3553 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3557 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3558 /* "Return @code{#t} if the list of parameters is monotonically\n"
3561 #define FUNC_NAME s_scm_gr_p
3563 scm_gr_p (SCM x
, SCM y
)
3565 if (!SCM_NUMBERP (x
))
3566 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3567 else if (!SCM_NUMBERP (y
))
3568 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3570 return scm_less_p (y
, x
);
3575 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3576 /* "Return @code{#t} if the list of parameters is monotonically\n"
3579 #define FUNC_NAME s_scm_leq_p
3581 scm_leq_p (SCM x
, SCM y
)
3583 if (!SCM_NUMBERP (x
))
3584 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3585 else if (!SCM_NUMBERP (y
))
3586 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3587 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3590 return scm_not (scm_less_p (y
, x
));
3595 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3596 /* "Return @code{#t} if the list of parameters is monotonically\n"
3599 #define FUNC_NAME s_scm_geq_p
3601 scm_geq_p (SCM x
, SCM y
)
3603 if (!SCM_NUMBERP (x
))
3604 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3605 else if (!SCM_NUMBERP (y
))
3606 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3607 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3610 return scm_not (scm_less_p (x
, y
));
3615 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3616 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3622 if (SCM_I_INUMP (z
))
3623 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3624 else if (SCM_BIGP (z
))
3626 else if (SCM_REALP (z
))
3627 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3628 else if (SCM_COMPLEXP (z
))
3629 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3630 && SCM_COMPLEX_IMAG (z
) == 0.0);
3631 else if (SCM_FRACTIONP (z
))
3634 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3638 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3639 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3643 scm_positive_p (SCM x
)
3645 if (SCM_I_INUMP (x
))
3646 return scm_from_bool (SCM_I_INUM (x
) > 0);
3647 else if (SCM_BIGP (x
))
3649 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3650 scm_remember_upto_here_1 (x
);
3651 return scm_from_bool (sgn
> 0);
3653 else if (SCM_REALP (x
))
3654 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3655 else if (SCM_FRACTIONP (x
))
3656 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3658 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3662 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3663 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3667 scm_negative_p (SCM x
)
3669 if (SCM_I_INUMP (x
))
3670 return scm_from_bool (SCM_I_INUM (x
) < 0);
3671 else if (SCM_BIGP (x
))
3673 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3674 scm_remember_upto_here_1 (x
);
3675 return scm_from_bool (sgn
< 0);
3677 else if (SCM_REALP (x
))
3678 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3679 else if (SCM_FRACTIONP (x
))
3680 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3682 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3686 /* scm_min and scm_max return an inexact when either argument is inexact, as
3687 required by r5rs. On that basis, for exact/inexact combinations the
3688 exact is converted to inexact to compare and possibly return. This is
3689 unlike scm_less_p above which takes some trouble to preserve all bits in
3690 its test, such trouble is not required for min and max. */
3692 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3693 /* "Return the maximum of all parameter values."
3696 scm_max (SCM x
, SCM y
)
3701 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3702 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3705 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3708 if (SCM_I_INUMP (x
))
3710 long xx
= SCM_I_INUM (x
);
3711 if (SCM_I_INUMP (y
))
3713 long yy
= SCM_I_INUM (y
);
3714 return (xx
< yy
) ? y
: x
;
3716 else if (SCM_BIGP (y
))
3718 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3719 scm_remember_upto_here_1 (y
);
3720 return (sgn
< 0) ? x
: y
;
3722 else if (SCM_REALP (y
))
3725 /* if y==NaN then ">" is false and we return NaN */
3726 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3728 else if (SCM_FRACTIONP (y
))
3731 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3734 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3736 else if (SCM_BIGP (x
))
3738 if (SCM_I_INUMP (y
))
3740 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3741 scm_remember_upto_here_1 (x
);
3742 return (sgn
< 0) ? y
: x
;
3744 else if (SCM_BIGP (y
))
3746 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3747 scm_remember_upto_here_2 (x
, y
);
3748 return (cmp
> 0) ? x
: y
;
3750 else if (SCM_REALP (y
))
3752 /* if y==NaN then xx>yy is false, so we return the NaN y */
3755 xx
= scm_i_big2dbl (x
);
3756 yy
= SCM_REAL_VALUE (y
);
3757 return (xx
> yy
? scm_from_double (xx
) : y
);
3759 else if (SCM_FRACTIONP (y
))
3764 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3766 else if (SCM_REALP (x
))
3768 if (SCM_I_INUMP (y
))
3770 double z
= SCM_I_INUM (y
);
3771 /* if x==NaN then "<" is false and we return NaN */
3772 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3774 else if (SCM_BIGP (y
))
3779 else if (SCM_REALP (y
))
3781 /* if x==NaN then our explicit check means we return NaN
3782 if y==NaN then ">" is false and we return NaN
3783 calling isnan is unavoidable, since it's the only way to know
3784 which of x or y causes any compares to be false */
3785 double xx
= SCM_REAL_VALUE (x
);
3786 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3788 else if (SCM_FRACTIONP (y
))
3790 double yy
= scm_i_fraction2double (y
);
3791 double xx
= SCM_REAL_VALUE (x
);
3792 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3795 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3797 else if (SCM_FRACTIONP (x
))
3799 if (SCM_I_INUMP (y
))
3803 else if (SCM_BIGP (y
))
3807 else if (SCM_REALP (y
))
3809 double xx
= scm_i_fraction2double (x
);
3810 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3812 else if (SCM_FRACTIONP (y
))
3817 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3820 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3824 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3825 /* "Return the minium of all parameter values."
3828 scm_min (SCM x
, SCM y
)
3833 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3834 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3837 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3840 if (SCM_I_INUMP (x
))
3842 long xx
= SCM_I_INUM (x
);
3843 if (SCM_I_INUMP (y
))
3845 long yy
= SCM_I_INUM (y
);
3846 return (xx
< yy
) ? x
: y
;
3848 else if (SCM_BIGP (y
))
3850 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3851 scm_remember_upto_here_1 (y
);
3852 return (sgn
< 0) ? y
: x
;
3854 else if (SCM_REALP (y
))
3857 /* if y==NaN then "<" is false and we return NaN */
3858 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3860 else if (SCM_FRACTIONP (y
))
3863 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3866 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3868 else if (SCM_BIGP (x
))
3870 if (SCM_I_INUMP (y
))
3872 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3873 scm_remember_upto_here_1 (x
);
3874 return (sgn
< 0) ? x
: y
;
3876 else if (SCM_BIGP (y
))
3878 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3879 scm_remember_upto_here_2 (x
, y
);
3880 return (cmp
> 0) ? y
: x
;
3882 else if (SCM_REALP (y
))
3884 /* if y==NaN then xx<yy is false, so we return the NaN y */
3887 xx
= scm_i_big2dbl (x
);
3888 yy
= SCM_REAL_VALUE (y
);
3889 return (xx
< yy
? scm_from_double (xx
) : y
);
3891 else if (SCM_FRACTIONP (y
))
3896 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3898 else if (SCM_REALP (x
))
3900 if (SCM_I_INUMP (y
))
3902 double z
= SCM_I_INUM (y
);
3903 /* if x==NaN then "<" is false and we return NaN */
3904 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3906 else if (SCM_BIGP (y
))
3911 else if (SCM_REALP (y
))
3913 /* if x==NaN then our explicit check means we return NaN
3914 if y==NaN then "<" is false and we return NaN
3915 calling isnan is unavoidable, since it's the only way to know
3916 which of x or y causes any compares to be false */
3917 double xx
= SCM_REAL_VALUE (x
);
3918 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3920 else if (SCM_FRACTIONP (y
))
3922 double yy
= scm_i_fraction2double (y
);
3923 double xx
= SCM_REAL_VALUE (x
);
3924 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3927 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3929 else if (SCM_FRACTIONP (x
))
3931 if (SCM_I_INUMP (y
))
3935 else if (SCM_BIGP (y
))
3939 else if (SCM_REALP (y
))
3941 double xx
= scm_i_fraction2double (x
);
3942 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3944 else if (SCM_FRACTIONP (y
))
3949 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3952 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3956 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3957 /* "Return the sum of all parameter values. Return 0 if called without\n"
3961 scm_sum (SCM x
, SCM y
)
3963 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
3965 if (SCM_NUMBERP (x
)) return x
;
3966 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3967 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3970 if (SCM_LIKELY (SCM_I_INUMP (x
)))
3972 if (SCM_LIKELY (SCM_I_INUMP (y
)))
3974 long xx
= SCM_I_INUM (x
);
3975 long yy
= SCM_I_INUM (y
);
3976 long int z
= xx
+ yy
;
3977 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3979 else if (SCM_BIGP (y
))
3984 else if (SCM_REALP (y
))
3986 long int xx
= SCM_I_INUM (x
);
3987 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3989 else if (SCM_COMPLEXP (y
))
3991 long int xx
= SCM_I_INUM (x
);
3992 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3993 SCM_COMPLEX_IMAG (y
));
3995 else if (SCM_FRACTIONP (y
))
3996 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3997 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3998 SCM_FRACTION_DENOMINATOR (y
));
4000 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4001 } else if (SCM_BIGP (x
))
4003 if (SCM_I_INUMP (y
))
4008 inum
= SCM_I_INUM (y
);
4011 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4014 SCM result
= scm_i_mkbig ();
4015 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4016 scm_remember_upto_here_1 (x
);
4017 /* we know the result will have to be a bignum */
4020 return scm_i_normbig (result
);
4024 SCM result
= scm_i_mkbig ();
4025 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4026 scm_remember_upto_here_1 (x
);
4027 /* we know the result will have to be a bignum */
4030 return scm_i_normbig (result
);
4033 else if (SCM_BIGP (y
))
4035 SCM result
= scm_i_mkbig ();
4036 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4037 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4038 mpz_add (SCM_I_BIG_MPZ (result
),
4041 scm_remember_upto_here_2 (x
, y
);
4042 /* we know the result will have to be a bignum */
4045 return scm_i_normbig (result
);
4047 else if (SCM_REALP (y
))
4049 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4050 scm_remember_upto_here_1 (x
);
4051 return scm_from_double (result
);
4053 else if (SCM_COMPLEXP (y
))
4055 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4056 + SCM_COMPLEX_REAL (y
));
4057 scm_remember_upto_here_1 (x
);
4058 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4060 else if (SCM_FRACTIONP (y
))
4061 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4062 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4063 SCM_FRACTION_DENOMINATOR (y
));
4065 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4067 else if (SCM_REALP (x
))
4069 if (SCM_I_INUMP (y
))
4070 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4071 else if (SCM_BIGP (y
))
4073 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4074 scm_remember_upto_here_1 (y
);
4075 return scm_from_double (result
);
4077 else if (SCM_REALP (y
))
4078 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4079 else if (SCM_COMPLEXP (y
))
4080 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4081 SCM_COMPLEX_IMAG (y
));
4082 else if (SCM_FRACTIONP (y
))
4083 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4085 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4087 else if (SCM_COMPLEXP (x
))
4089 if (SCM_I_INUMP (y
))
4090 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4091 SCM_COMPLEX_IMAG (x
));
4092 else if (SCM_BIGP (y
))
4094 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4095 + SCM_COMPLEX_REAL (x
));
4096 scm_remember_upto_here_1 (y
);
4097 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4099 else if (SCM_REALP (y
))
4100 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4101 SCM_COMPLEX_IMAG (x
));
4102 else if (SCM_COMPLEXP (y
))
4103 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4104 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4105 else if (SCM_FRACTIONP (y
))
4106 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4107 SCM_COMPLEX_IMAG (x
));
4109 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4111 else if (SCM_FRACTIONP (x
))
4113 if (SCM_I_INUMP (y
))
4114 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4115 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4116 SCM_FRACTION_DENOMINATOR (x
));
4117 else if (SCM_BIGP (y
))
4118 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4119 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4120 SCM_FRACTION_DENOMINATOR (x
));
4121 else if (SCM_REALP (y
))
4122 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4123 else if (SCM_COMPLEXP (y
))
4124 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4125 SCM_COMPLEX_IMAG (y
));
4126 else if (SCM_FRACTIONP (y
))
4127 /* a/b + c/d = (ad + bc) / bd */
4128 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4129 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4130 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4132 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4135 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4139 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4141 "Return @math{@var{x}+1}.")
4142 #define FUNC_NAME s_scm_oneplus
4144 return scm_sum (x
, SCM_I_MAKINUM (1));
4149 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4150 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4151 * the sum of all but the first argument are subtracted from the first
4153 #define FUNC_NAME s_difference
4155 scm_difference (SCM x
, SCM y
)
4157 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4160 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4162 if (SCM_I_INUMP (x
))
4164 long xx
= -SCM_I_INUM (x
);
4165 if (SCM_FIXABLE (xx
))
4166 return SCM_I_MAKINUM (xx
);
4168 return scm_i_long2big (xx
);
4170 else if (SCM_BIGP (x
))
4171 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4172 bignum, but negating that gives a fixnum. */
4173 return scm_i_normbig (scm_i_clonebig (x
, 0));
4174 else if (SCM_REALP (x
))
4175 return scm_from_double (-SCM_REAL_VALUE (x
));
4176 else if (SCM_COMPLEXP (x
))
4177 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4178 -SCM_COMPLEX_IMAG (x
));
4179 else if (SCM_FRACTIONP (x
))
4180 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4181 SCM_FRACTION_DENOMINATOR (x
));
4183 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4186 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4188 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4190 long int xx
= SCM_I_INUM (x
);
4191 long int yy
= SCM_I_INUM (y
);
4192 long int z
= xx
- yy
;
4193 if (SCM_FIXABLE (z
))
4194 return SCM_I_MAKINUM (z
);
4196 return scm_i_long2big (z
);
4198 else if (SCM_BIGP (y
))
4200 /* inum-x - big-y */
4201 long xx
= SCM_I_INUM (x
);
4204 return scm_i_clonebig (y
, 0);
4207 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4208 SCM result
= scm_i_mkbig ();
4211 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4214 /* x - y == -(y + -x) */
4215 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4216 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4218 scm_remember_upto_here_1 (y
);
4220 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4221 /* we know the result will have to be a bignum */
4224 return scm_i_normbig (result
);
4227 else if (SCM_REALP (y
))
4229 long int xx
= SCM_I_INUM (x
);
4230 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4232 else if (SCM_COMPLEXP (y
))
4234 long int xx
= SCM_I_INUM (x
);
4235 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4236 - SCM_COMPLEX_IMAG (y
));
4238 else if (SCM_FRACTIONP (y
))
4239 /* a - b/c = (ac - b) / c */
4240 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4241 SCM_FRACTION_NUMERATOR (y
)),
4242 SCM_FRACTION_DENOMINATOR (y
));
4244 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4246 else if (SCM_BIGP (x
))
4248 if (SCM_I_INUMP (y
))
4250 /* big-x - inum-y */
4251 long yy
= SCM_I_INUM (y
);
4252 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4254 scm_remember_upto_here_1 (x
);
4256 return (SCM_FIXABLE (-yy
) ?
4257 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4260 SCM result
= scm_i_mkbig ();
4263 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4265 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4266 scm_remember_upto_here_1 (x
);
4268 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4269 /* we know the result will have to be a bignum */
4272 return scm_i_normbig (result
);
4275 else if (SCM_BIGP (y
))
4277 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4278 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4279 SCM result
= scm_i_mkbig ();
4280 mpz_sub (SCM_I_BIG_MPZ (result
),
4283 scm_remember_upto_here_2 (x
, y
);
4284 /* we know the result will have to be a bignum */
4285 if ((sgn_x
== 1) && (sgn_y
== -1))
4287 if ((sgn_x
== -1) && (sgn_y
== 1))
4289 return scm_i_normbig (result
);
4291 else if (SCM_REALP (y
))
4293 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4294 scm_remember_upto_here_1 (x
);
4295 return scm_from_double (result
);
4297 else if (SCM_COMPLEXP (y
))
4299 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4300 - SCM_COMPLEX_REAL (y
));
4301 scm_remember_upto_here_1 (x
);
4302 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4304 else if (SCM_FRACTIONP (y
))
4305 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4306 SCM_FRACTION_NUMERATOR (y
)),
4307 SCM_FRACTION_DENOMINATOR (y
));
4308 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4310 else if (SCM_REALP (x
))
4312 if (SCM_I_INUMP (y
))
4313 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4314 else if (SCM_BIGP (y
))
4316 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4317 scm_remember_upto_here_1 (x
);
4318 return scm_from_double (result
);
4320 else if (SCM_REALP (y
))
4321 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4322 else if (SCM_COMPLEXP (y
))
4323 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4324 -SCM_COMPLEX_IMAG (y
));
4325 else if (SCM_FRACTIONP (y
))
4326 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4328 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4330 else if (SCM_COMPLEXP (x
))
4332 if (SCM_I_INUMP (y
))
4333 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4334 SCM_COMPLEX_IMAG (x
));
4335 else if (SCM_BIGP (y
))
4337 double real_part
= (SCM_COMPLEX_REAL (x
)
4338 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4339 scm_remember_upto_here_1 (x
);
4340 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4342 else if (SCM_REALP (y
))
4343 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4344 SCM_COMPLEX_IMAG (x
));
4345 else if (SCM_COMPLEXP (y
))
4346 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4347 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4348 else if (SCM_FRACTIONP (y
))
4349 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4350 SCM_COMPLEX_IMAG (x
));
4352 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4354 else if (SCM_FRACTIONP (x
))
4356 if (SCM_I_INUMP (y
))
4357 /* a/b - c = (a - cb) / b */
4358 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4359 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4360 SCM_FRACTION_DENOMINATOR (x
));
4361 else if (SCM_BIGP (y
))
4362 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4363 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4364 SCM_FRACTION_DENOMINATOR (x
));
4365 else if (SCM_REALP (y
))
4366 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4367 else if (SCM_COMPLEXP (y
))
4368 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4369 -SCM_COMPLEX_IMAG (y
));
4370 else if (SCM_FRACTIONP (y
))
4371 /* a/b - c/d = (ad - bc) / bd */
4372 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4373 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4374 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4376 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4379 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4384 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4386 "Return @math{@var{x}-1}.")
4387 #define FUNC_NAME s_scm_oneminus
4389 return scm_difference (x
, SCM_I_MAKINUM (1));
4394 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4395 /* "Return the product of all arguments. If called without arguments,\n"
4399 scm_product (SCM x
, SCM y
)
4401 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4404 return SCM_I_MAKINUM (1L);
4405 else if (SCM_NUMBERP (x
))
4408 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4411 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4416 xx
= SCM_I_INUM (x
);
4420 case 0: return x
; break;
4421 case 1: return y
; break;
4424 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4426 long yy
= SCM_I_INUM (y
);
4428 SCM k
= SCM_I_MAKINUM (kk
);
4429 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4433 SCM result
= scm_i_long2big (xx
);
4434 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4435 return scm_i_normbig (result
);
4438 else if (SCM_BIGP (y
))
4440 SCM result
= scm_i_mkbig ();
4441 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4442 scm_remember_upto_here_1 (y
);
4445 else if (SCM_REALP (y
))
4446 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4447 else if (SCM_COMPLEXP (y
))
4448 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4449 xx
* SCM_COMPLEX_IMAG (y
));
4450 else if (SCM_FRACTIONP (y
))
4451 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4452 SCM_FRACTION_DENOMINATOR (y
));
4454 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4456 else if (SCM_BIGP (x
))
4458 if (SCM_I_INUMP (y
))
4463 else if (SCM_BIGP (y
))
4465 SCM result
= scm_i_mkbig ();
4466 mpz_mul (SCM_I_BIG_MPZ (result
),
4469 scm_remember_upto_here_2 (x
, y
);
4472 else if (SCM_REALP (y
))
4474 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4475 scm_remember_upto_here_1 (x
);
4476 return scm_from_double (result
);
4478 else if (SCM_COMPLEXP (y
))
4480 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4481 scm_remember_upto_here_1 (x
);
4482 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4483 z
* SCM_COMPLEX_IMAG (y
));
4485 else if (SCM_FRACTIONP (y
))
4486 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4487 SCM_FRACTION_DENOMINATOR (y
));
4489 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4491 else if (SCM_REALP (x
))
4493 if (SCM_I_INUMP (y
))
4495 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4496 if (scm_is_eq (y
, SCM_INUM0
))
4498 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4500 else if (SCM_BIGP (y
))
4502 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4503 scm_remember_upto_here_1 (y
);
4504 return scm_from_double (result
);
4506 else if (SCM_REALP (y
))
4507 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4508 else if (SCM_COMPLEXP (y
))
4509 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4510 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4511 else if (SCM_FRACTIONP (y
))
4512 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4514 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4516 else if (SCM_COMPLEXP (x
))
4518 if (SCM_I_INUMP (y
))
4520 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4521 if (scm_is_eq (y
, SCM_INUM0
))
4523 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4524 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4526 else if (SCM_BIGP (y
))
4528 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4529 scm_remember_upto_here_1 (y
);
4530 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4531 z
* SCM_COMPLEX_IMAG (x
));
4533 else if (SCM_REALP (y
))
4534 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4535 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4536 else if (SCM_COMPLEXP (y
))
4538 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4539 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4540 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4541 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4543 else if (SCM_FRACTIONP (y
))
4545 double yy
= scm_i_fraction2double (y
);
4546 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4547 yy
* SCM_COMPLEX_IMAG (x
));
4550 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4552 else if (SCM_FRACTIONP (x
))
4554 if (SCM_I_INUMP (y
))
4555 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4556 SCM_FRACTION_DENOMINATOR (x
));
4557 else if (SCM_BIGP (y
))
4558 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4559 SCM_FRACTION_DENOMINATOR (x
));
4560 else if (SCM_REALP (y
))
4561 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4562 else if (SCM_COMPLEXP (y
))
4564 double xx
= scm_i_fraction2double (x
);
4565 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4566 xx
* SCM_COMPLEX_IMAG (y
));
4568 else if (SCM_FRACTIONP (y
))
4569 /* a/b * c/d = ac / bd */
4570 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4571 SCM_FRACTION_NUMERATOR (y
)),
4572 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4573 SCM_FRACTION_DENOMINATOR (y
)));
4575 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4578 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4581 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4582 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4583 #define ALLOW_DIVIDE_BY_ZERO
4584 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4587 /* The code below for complex division is adapted from the GNU
4588 libstdc++, which adapted it from f2c's libF77, and is subject to
4591 /****************************************************************
4592 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4594 Permission to use, copy, modify, and distribute this software
4595 and its documentation for any purpose and without fee is hereby
4596 granted, provided that the above copyright notice appear in all
4597 copies and that both that the copyright notice and this
4598 permission notice and warranty disclaimer appear in supporting
4599 documentation, and that the names of AT&T Bell Laboratories or
4600 Bellcore or any of their entities not be used in advertising or
4601 publicity pertaining to distribution of the software without
4602 specific, written prior permission.
4604 AT&T and Bellcore disclaim all warranties with regard to this
4605 software, including all implied warranties of merchantability
4606 and fitness. In no event shall AT&T or Bellcore be liable for
4607 any special, indirect or consequential damages or any damages
4608 whatsoever resulting from loss of use, data or profits, whether
4609 in an action of contract, negligence or other tortious action,
4610 arising out of or in connection with the use or performance of
4612 ****************************************************************/
4614 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4615 /* Divide the first argument by the product of the remaining
4616 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4618 #define FUNC_NAME s_divide
4620 scm_i_divide (SCM x
, SCM y
, int inexact
)
4624 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4627 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4628 else if (SCM_I_INUMP (x
))
4630 long xx
= SCM_I_INUM (x
);
4631 if (xx
== 1 || xx
== -1)
4633 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4635 scm_num_overflow (s_divide
);
4640 return scm_from_double (1.0 / (double) xx
);
4641 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4644 else if (SCM_BIGP (x
))
4647 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4648 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4650 else if (SCM_REALP (x
))
4652 double xx
= SCM_REAL_VALUE (x
);
4653 #ifndef ALLOW_DIVIDE_BY_ZERO
4655 scm_num_overflow (s_divide
);
4658 return scm_from_double (1.0 / xx
);
4660 else if (SCM_COMPLEXP (x
))
4662 double r
= SCM_COMPLEX_REAL (x
);
4663 double i
= SCM_COMPLEX_IMAG (x
);
4664 if (fabs(r
) <= fabs(i
))
4667 double d
= i
* (1.0 + t
* t
);
4668 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4673 double d
= r
* (1.0 + t
* t
);
4674 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4677 else if (SCM_FRACTIONP (x
))
4678 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4679 SCM_FRACTION_NUMERATOR (x
));
4681 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4684 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4686 long xx
= SCM_I_INUM (x
);
4687 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4689 long yy
= SCM_I_INUM (y
);
4692 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4693 scm_num_overflow (s_divide
);
4695 return scm_from_double ((double) xx
/ (double) yy
);
4698 else if (xx
% yy
!= 0)
4701 return scm_from_double ((double) xx
/ (double) yy
);
4702 else return scm_i_make_ratio (x
, y
);
4707 if (SCM_FIXABLE (z
))
4708 return SCM_I_MAKINUM (z
);
4710 return scm_i_long2big (z
);
4713 else if (SCM_BIGP (y
))
4716 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4717 else return scm_i_make_ratio (x
, y
);
4719 else if (SCM_REALP (y
))
4721 double yy
= SCM_REAL_VALUE (y
);
4722 #ifndef ALLOW_DIVIDE_BY_ZERO
4724 scm_num_overflow (s_divide
);
4727 return scm_from_double ((double) xx
/ yy
);
4729 else if (SCM_COMPLEXP (y
))
4732 complex_div
: /* y _must_ be a complex number */
4734 double r
= SCM_COMPLEX_REAL (y
);
4735 double i
= SCM_COMPLEX_IMAG (y
);
4736 if (fabs(r
) <= fabs(i
))
4739 double d
= i
* (1.0 + t
* t
);
4740 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4745 double d
= r
* (1.0 + t
* t
);
4746 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4750 else if (SCM_FRACTIONP (y
))
4751 /* a / b/c = ac / b */
4752 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4753 SCM_FRACTION_NUMERATOR (y
));
4755 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4757 else if (SCM_BIGP (x
))
4759 if (SCM_I_INUMP (y
))
4761 long int yy
= SCM_I_INUM (y
);
4764 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4765 scm_num_overflow (s_divide
);
4767 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4768 scm_remember_upto_here_1 (x
);
4769 return (sgn
== 0) ? scm_nan () : scm_inf ();
4776 /* FIXME: HMM, what are the relative performance issues here?
4777 We need to test. Is it faster on average to test
4778 divisible_p, then perform whichever operation, or is it
4779 faster to perform the integer div opportunistically and
4780 switch to real if there's a remainder? For now we take the
4781 middle ground: test, then if divisible, use the faster div
4784 long abs_yy
= yy
< 0 ? -yy
: yy
;
4785 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4789 SCM result
= scm_i_mkbig ();
4790 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4791 scm_remember_upto_here_1 (x
);
4793 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4794 return scm_i_normbig (result
);
4799 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4800 else return scm_i_make_ratio (x
, y
);
4804 else if (SCM_BIGP (y
))
4806 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4809 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4810 scm_num_overflow (s_divide
);
4812 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4813 scm_remember_upto_here_1 (x
);
4814 return (sgn
== 0) ? scm_nan () : scm_inf ();
4822 /* It's easily possible for the ratio x/y to fit a double
4823 but one or both x and y be too big to fit a double,
4824 hence the use of mpq_get_d rather than converting and
4827 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
4828 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
4829 return scm_from_double (mpq_get_d (q
));
4833 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4837 SCM result
= scm_i_mkbig ();
4838 mpz_divexact (SCM_I_BIG_MPZ (result
),
4841 scm_remember_upto_here_2 (x
, y
);
4842 return scm_i_normbig (result
);
4845 return scm_i_make_ratio (x
, y
);
4849 else if (SCM_REALP (y
))
4851 double yy
= SCM_REAL_VALUE (y
);
4852 #ifndef ALLOW_DIVIDE_BY_ZERO
4854 scm_num_overflow (s_divide
);
4857 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4859 else if (SCM_COMPLEXP (y
))
4861 a
= scm_i_big2dbl (x
);
4864 else if (SCM_FRACTIONP (y
))
4865 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4866 SCM_FRACTION_NUMERATOR (y
));
4868 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4870 else if (SCM_REALP (x
))
4872 double rx
= SCM_REAL_VALUE (x
);
4873 if (SCM_I_INUMP (y
))
4875 long int yy
= SCM_I_INUM (y
);
4876 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4878 scm_num_overflow (s_divide
);
4881 return scm_from_double (rx
/ (double) yy
);
4883 else if (SCM_BIGP (y
))
4885 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4886 scm_remember_upto_here_1 (y
);
4887 return scm_from_double (rx
/ dby
);
4889 else if (SCM_REALP (y
))
4891 double yy
= SCM_REAL_VALUE (y
);
4892 #ifndef ALLOW_DIVIDE_BY_ZERO
4894 scm_num_overflow (s_divide
);
4897 return scm_from_double (rx
/ yy
);
4899 else if (SCM_COMPLEXP (y
))
4904 else if (SCM_FRACTIONP (y
))
4905 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4907 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4909 else if (SCM_COMPLEXP (x
))
4911 double rx
= SCM_COMPLEX_REAL (x
);
4912 double ix
= SCM_COMPLEX_IMAG (x
);
4913 if (SCM_I_INUMP (y
))
4915 long int yy
= SCM_I_INUM (y
);
4916 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4918 scm_num_overflow (s_divide
);
4923 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4926 else if (SCM_BIGP (y
))
4928 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4929 scm_remember_upto_here_1 (y
);
4930 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4932 else if (SCM_REALP (y
))
4934 double yy
= SCM_REAL_VALUE (y
);
4935 #ifndef ALLOW_DIVIDE_BY_ZERO
4937 scm_num_overflow (s_divide
);
4940 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4942 else if (SCM_COMPLEXP (y
))
4944 double ry
= SCM_COMPLEX_REAL (y
);
4945 double iy
= SCM_COMPLEX_IMAG (y
);
4946 if (fabs(ry
) <= fabs(iy
))
4949 double d
= iy
* (1.0 + t
* t
);
4950 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4955 double d
= ry
* (1.0 + t
* t
);
4956 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4959 else if (SCM_FRACTIONP (y
))
4961 double yy
= scm_i_fraction2double (y
);
4962 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4965 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4967 else if (SCM_FRACTIONP (x
))
4969 if (SCM_I_INUMP (y
))
4971 long int yy
= SCM_I_INUM (y
);
4972 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4974 scm_num_overflow (s_divide
);
4977 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4978 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4980 else if (SCM_BIGP (y
))
4982 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4983 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4985 else if (SCM_REALP (y
))
4987 double yy
= SCM_REAL_VALUE (y
);
4988 #ifndef ALLOW_DIVIDE_BY_ZERO
4990 scm_num_overflow (s_divide
);
4993 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4995 else if (SCM_COMPLEXP (y
))
4997 a
= scm_i_fraction2double (x
);
5000 else if (SCM_FRACTIONP (y
))
5001 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5002 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5004 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5007 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5011 scm_divide (SCM x
, SCM y
)
5013 return scm_i_divide (x
, y
, 0);
5016 static SCM
scm_divide2real (SCM x
, SCM y
)
5018 return scm_i_divide (x
, y
, 1);
5024 scm_asinh (double x
)
5029 #define asinh scm_asinh
5030 return log (x
+ sqrt (x
* x
+ 1));
5033 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
5034 /* "Return the inverse hyperbolic sine of @var{x}."
5039 scm_acosh (double x
)
5044 #define acosh scm_acosh
5045 return log (x
+ sqrt (x
* x
- 1));
5048 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
5049 /* "Return the inverse hyperbolic cosine of @var{x}."
5054 scm_atanh (double x
)
5059 #define atanh scm_atanh
5060 return 0.5 * log ((1 + x
) / (1 - x
));
5063 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
5064 /* "Return the inverse hyperbolic tangent of @var{x}."
5069 scm_c_truncate (double x
)
5080 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5081 half-way case (ie. when x is an integer plus 0.5) going upwards.
5082 Then half-way cases are identified and adjusted down if the
5083 round-upwards didn't give the desired even integer.
5085 "plus_half == result" identifies a half-way case. If plus_half, which is
5086 x + 0.5, is an integer then x must be an integer plus 0.5.
5088 An odd "result" value is identified with result/2 != floor(result/2).
5089 This is done with plus_half, since that value is ready for use sooner in
5090 a pipelined cpu, and we're already requiring plus_half == result.
5092 Note however that we need to be careful when x is big and already an
5093 integer. In that case "x+0.5" may round to an adjacent integer, causing
5094 us to return such a value, incorrectly. For instance if the hardware is
5095 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5096 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5097 returned. Or if the hardware is in round-upwards mode, then other bigger
5098 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5099 representable value, 2^128+2^76 (or whatever), again incorrect.
5101 These bad roundings of x+0.5 are avoided by testing at the start whether
5102 x is already an integer. If it is then clearly that's the desired result
5103 already. And if it's not then the exponent must be small enough to allow
5104 an 0.5 to be represented, and hence added without a bad rounding. */
5107 scm_c_round (double x
)
5109 double plus_half
, result
;
5114 plus_half
= x
+ 0.5;
5115 result
= floor (plus_half
);
5116 /* Adjust so that the rounding is towards even. */
5117 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5122 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5124 "Round the number @var{x} towards zero.")
5125 #define FUNC_NAME s_scm_truncate_number
5127 if (scm_is_false (scm_negative_p (x
)))
5128 return scm_floor (x
);
5130 return scm_ceiling (x
);
5134 static SCM exactly_one_half
;
5136 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5138 "Round the number @var{x} towards the nearest integer. "
5139 "When it is exactly halfway between two integers, "
5140 "round towards the even one.")
5141 #define FUNC_NAME s_scm_round_number
5143 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5145 else if (SCM_REALP (x
))
5146 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5149 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5150 single quotient+remainder division then examining to see which way
5151 the rounding should go. */
5152 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5153 SCM result
= scm_floor (plus_half
);
5154 /* Adjust so that the rounding is towards even. */
5155 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5156 && scm_is_true (scm_odd_p (result
)))
5157 return scm_difference (result
, SCM_I_MAKINUM (1));
5164 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5166 "Round the number @var{x} towards minus infinity.")
5167 #define FUNC_NAME s_scm_floor
5169 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5171 else if (SCM_REALP (x
))
5172 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5173 else if (SCM_FRACTIONP (x
))
5175 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5176 SCM_FRACTION_DENOMINATOR (x
));
5177 if (scm_is_false (scm_negative_p (x
)))
5179 /* For positive x, rounding towards zero is correct. */
5184 /* For negative x, we need to return q-1 unless x is an
5185 integer. But fractions are never integer, per our
5187 return scm_difference (q
, SCM_I_MAKINUM (1));
5191 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5195 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5197 "Round the number @var{x} towards infinity.")
5198 #define FUNC_NAME s_scm_ceiling
5200 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5202 else if (SCM_REALP (x
))
5203 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5204 else if (SCM_FRACTIONP (x
))
5206 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5207 SCM_FRACTION_DENOMINATOR (x
));
5208 if (scm_is_false (scm_positive_p (x
)))
5210 /* For negative x, rounding towards zero is correct. */
5215 /* For positive x, we need to return q+1 unless x is an
5216 integer. But fractions are never integer, per our
5218 return scm_sum (q
, SCM_I_MAKINUM (1));
5222 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5226 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5227 /* "Return the square root of the real number @var{x}."
5229 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5230 /* "Return the absolute value of the real number @var{x}."
5232 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5233 /* "Return the @var{x}th power of e."
5235 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5236 /* "Return the natural logarithm of the real number @var{x}."
5238 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5239 /* "Return the sine of the real number @var{x}."
5241 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5242 /* "Return the cosine of the real number @var{x}."
5244 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5245 /* "Return the tangent of the real number @var{x}."
5247 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5248 /* "Return the arc sine of the real number @var{x}."
5250 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5251 /* "Return the arc cosine of the real number @var{x}."
5253 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5254 /* "Return the arc tangent of the real number @var{x}."
5256 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5257 /* "Return the hyperbolic sine of the real number @var{x}."
5259 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5260 /* "Return the hyperbolic cosine of the real number @var{x}."
5262 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5263 /* "Return the hyperbolic tangent of the real number @var{x}."
5271 static void scm_two_doubles (SCM x
,
5273 const char *sstring
,
5277 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5279 if (SCM_I_INUMP (x
))
5280 xy
->x
= SCM_I_INUM (x
);
5281 else if (SCM_BIGP (x
))
5282 xy
->x
= scm_i_big2dbl (x
);
5283 else if (SCM_REALP (x
))
5284 xy
->x
= SCM_REAL_VALUE (x
);
5285 else if (SCM_FRACTIONP (x
))
5286 xy
->x
= scm_i_fraction2double (x
);
5288 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5290 if (SCM_I_INUMP (y
))
5291 xy
->y
= SCM_I_INUM (y
);
5292 else if (SCM_BIGP (y
))
5293 xy
->y
= scm_i_big2dbl (y
);
5294 else if (SCM_REALP (y
))
5295 xy
->y
= SCM_REAL_VALUE (y
);
5296 else if (SCM_FRACTIONP (y
))
5297 xy
->y
= scm_i_fraction2double (y
);
5299 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5303 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5305 "Return @var{x} raised to the power of @var{y}. This\n"
5306 "procedure does not accept complex arguments.")
5307 #define FUNC_NAME s_scm_sys_expt
5310 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5311 return scm_from_double (pow (xy
.x
, xy
.y
));
5316 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5318 "Return the arc tangent of the two arguments @var{x} and\n"
5319 "@var{y}. This is similar to calculating the arc tangent of\n"
5320 "@var{x} / @var{y}, except that the signs of both arguments\n"
5321 "are used to determine the quadrant of the result. This\n"
5322 "procedure does not accept complex arguments.")
5323 #define FUNC_NAME s_scm_sys_atan2
5326 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5327 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5332 scm_c_make_rectangular (double re
, double im
)
5335 return scm_from_double (re
);
5339 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5341 SCM_COMPLEX_REAL (z
) = re
;
5342 SCM_COMPLEX_IMAG (z
) = im
;
5347 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5348 (SCM real_part
, SCM imaginary_part
),
5349 "Return a complex number constructed of the given @var{real-part} "
5350 "and @var{imaginary-part} parts.")
5351 #define FUNC_NAME s_scm_make_rectangular
5354 scm_two_doubles (real_part
, imaginary_part
, FUNC_NAME
, &xy
);
5355 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5360 scm_c_make_polar (double mag
, double ang
)
5364 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5365 use it on Glibc-based systems that have it (it's a GNU extension). See
5366 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5368 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5369 sincos (ang
, &s
, &c
);
5374 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5377 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5379 "Return the complex number @var{x} * e^(i * @var{y}).")
5380 #define FUNC_NAME s_scm_make_polar
5383 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5384 return scm_c_make_polar (xy
.x
, xy
.y
);
5389 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5390 /* "Return the real part of the number @var{z}."
5393 scm_real_part (SCM z
)
5395 if (SCM_I_INUMP (z
))
5397 else if (SCM_BIGP (z
))
5399 else if (SCM_REALP (z
))
5401 else if (SCM_COMPLEXP (z
))
5402 return scm_from_double (SCM_COMPLEX_REAL (z
));
5403 else if (SCM_FRACTIONP (z
))
5406 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5410 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5411 /* "Return the imaginary part of the number @var{z}."
5414 scm_imag_part (SCM z
)
5416 if (SCM_I_INUMP (z
))
5418 else if (SCM_BIGP (z
))
5420 else if (SCM_REALP (z
))
5422 else if (SCM_COMPLEXP (z
))
5423 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5424 else if (SCM_FRACTIONP (z
))
5427 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5430 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5431 /* "Return the numerator of the number @var{z}."
5434 scm_numerator (SCM z
)
5436 if (SCM_I_INUMP (z
))
5438 else if (SCM_BIGP (z
))
5440 else if (SCM_FRACTIONP (z
))
5441 return SCM_FRACTION_NUMERATOR (z
);
5442 else if (SCM_REALP (z
))
5443 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5445 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5449 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5450 /* "Return the denominator of the number @var{z}."
5453 scm_denominator (SCM z
)
5455 if (SCM_I_INUMP (z
))
5456 return SCM_I_MAKINUM (1);
5457 else if (SCM_BIGP (z
))
5458 return SCM_I_MAKINUM (1);
5459 else if (SCM_FRACTIONP (z
))
5460 return SCM_FRACTION_DENOMINATOR (z
);
5461 else if (SCM_REALP (z
))
5462 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5464 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5467 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5468 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5469 * "@code{abs} for real arguments, but also allows complex numbers."
5472 scm_magnitude (SCM z
)
5474 if (SCM_I_INUMP (z
))
5476 long int zz
= SCM_I_INUM (z
);
5479 else if (SCM_POSFIXABLE (-zz
))
5480 return SCM_I_MAKINUM (-zz
);
5482 return scm_i_long2big (-zz
);
5484 else if (SCM_BIGP (z
))
5486 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5487 scm_remember_upto_here_1 (z
);
5489 return scm_i_clonebig (z
, 0);
5493 else if (SCM_REALP (z
))
5494 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5495 else if (SCM_COMPLEXP (z
))
5496 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5497 else if (SCM_FRACTIONP (z
))
5499 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5501 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5502 SCM_FRACTION_DENOMINATOR (z
));
5505 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5509 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5510 /* "Return the angle of the complex number @var{z}."
5515 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5516 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5517 But if atan2 follows the floating point rounding mode, then the value
5518 is not a constant. Maybe it'd be close enough though. */
5519 if (SCM_I_INUMP (z
))
5521 if (SCM_I_INUM (z
) >= 0)
5524 return scm_from_double (atan2 (0.0, -1.0));
5526 else if (SCM_BIGP (z
))
5528 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5529 scm_remember_upto_here_1 (z
);
5531 return scm_from_double (atan2 (0.0, -1.0));
5535 else if (SCM_REALP (z
))
5537 if (SCM_REAL_VALUE (z
) >= 0)
5540 return scm_from_double (atan2 (0.0, -1.0));
5542 else if (SCM_COMPLEXP (z
))
5543 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5544 else if (SCM_FRACTIONP (z
))
5546 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5548 else return scm_from_double (atan2 (0.0, -1.0));
5551 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5555 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5556 /* Convert the number @var{x} to its inexact representation.\n"
5559 scm_exact_to_inexact (SCM z
)
5561 if (SCM_I_INUMP (z
))
5562 return scm_from_double ((double) SCM_I_INUM (z
));
5563 else if (SCM_BIGP (z
))
5564 return scm_from_double (scm_i_big2dbl (z
));
5565 else if (SCM_FRACTIONP (z
))
5566 return scm_from_double (scm_i_fraction2double (z
));
5567 else if (SCM_INEXACTP (z
))
5570 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5574 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5576 "Return an exact number that is numerically closest to @var{z}.")
5577 #define FUNC_NAME s_scm_inexact_to_exact
5579 if (SCM_I_INUMP (z
))
5581 else if (SCM_BIGP (z
))
5583 else if (SCM_REALP (z
))
5585 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5586 SCM_OUT_OF_RANGE (1, z
);
5593 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5594 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5595 scm_i_mpz2num (mpq_denref (frac
)));
5597 /* When scm_i_make_ratio throws, we leak the memory allocated
5604 else if (SCM_FRACTIONP (z
))
5607 SCM_WRONG_TYPE_ARG (1, z
);
5611 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5613 "Returns the @emph{simplest} rational number differing\n"
5614 "from @var{x} by no more than @var{eps}.\n"
5616 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
5617 "exact result when both its arguments are exact. Thus, you might need\n"
5618 "to use @code{inexact->exact} on the arguments.\n"
5621 "(rationalize (inexact->exact 1.2) 1/100)\n"
5624 #define FUNC_NAME s_scm_rationalize
5626 if (SCM_I_INUMP (x
))
5628 else if (SCM_BIGP (x
))
5630 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5632 /* Use continued fractions to find closest ratio. All
5633 arithmetic is done with exact numbers.
5636 SCM ex
= scm_inexact_to_exact (x
);
5637 SCM int_part
= scm_floor (ex
);
5638 SCM tt
= SCM_I_MAKINUM (1);
5639 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5640 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5644 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5647 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5648 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5650 /* We stop after a million iterations just to be absolutely sure
5651 that we don't go into an infinite loop. The process normally
5652 converges after less than a dozen iterations.
5655 eps
= scm_abs (eps
);
5656 while (++i
< 1000000)
5658 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5659 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5660 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5662 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5663 eps
))) /* abs(x-a/b) <= eps */
5665 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5666 if (scm_is_false (scm_exact_p (x
))
5667 || scm_is_false (scm_exact_p (eps
)))
5668 return scm_exact_to_inexact (res
);
5672 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5674 tt
= scm_floor (rx
); /* tt = floor (rx) */
5680 scm_num_overflow (s_scm_rationalize
);
5683 SCM_WRONG_TYPE_ARG (1, x
);
5687 /* conversion functions */
5690 scm_is_integer (SCM val
)
5692 return scm_is_true (scm_integer_p (val
));
5696 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5698 if (SCM_I_INUMP (val
))
5700 scm_t_signed_bits n
= SCM_I_INUM (val
);
5701 return n
>= min
&& n
<= max
;
5703 else if (SCM_BIGP (val
))
5705 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5707 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5709 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5711 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5712 return n
>= min
&& n
<= max
;
5722 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5723 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5726 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5727 SCM_I_BIG_MPZ (val
));
5729 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5741 return n
>= min
&& n
<= max
;
5749 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5751 if (SCM_I_INUMP (val
))
5753 scm_t_signed_bits n
= SCM_I_INUM (val
);
5754 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5756 else if (SCM_BIGP (val
))
5758 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5760 else if (max
<= ULONG_MAX
)
5762 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5764 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5765 return n
>= min
&& n
<= max
;
5775 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5778 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5779 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5782 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5783 SCM_I_BIG_MPZ (val
));
5785 return n
>= min
&& n
<= max
;
5793 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5795 scm_error (scm_out_of_range_key
,
5797 "Value out of range ~S to ~S: ~S",
5798 scm_list_3 (min
, max
, bad_val
),
5799 scm_list_1 (bad_val
));
5802 #define TYPE scm_t_intmax
5803 #define TYPE_MIN min
5804 #define TYPE_MAX max
5805 #define SIZEOF_TYPE 0
5806 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5807 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5808 #include "libguile/conv-integer.i.c"
5810 #define TYPE scm_t_uintmax
5811 #define TYPE_MIN min
5812 #define TYPE_MAX max
5813 #define SIZEOF_TYPE 0
5814 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5815 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5816 #include "libguile/conv-uinteger.i.c"
5818 #define TYPE scm_t_int8
5819 #define TYPE_MIN SCM_T_INT8_MIN
5820 #define TYPE_MAX SCM_T_INT8_MAX
5821 #define SIZEOF_TYPE 1
5822 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5823 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5824 #include "libguile/conv-integer.i.c"
5826 #define TYPE scm_t_uint8
5828 #define TYPE_MAX SCM_T_UINT8_MAX
5829 #define SIZEOF_TYPE 1
5830 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5831 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5832 #include "libguile/conv-uinteger.i.c"
5834 #define TYPE scm_t_int16
5835 #define TYPE_MIN SCM_T_INT16_MIN
5836 #define TYPE_MAX SCM_T_INT16_MAX
5837 #define SIZEOF_TYPE 2
5838 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5839 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5840 #include "libguile/conv-integer.i.c"
5842 #define TYPE scm_t_uint16
5844 #define TYPE_MAX SCM_T_UINT16_MAX
5845 #define SIZEOF_TYPE 2
5846 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5847 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5848 #include "libguile/conv-uinteger.i.c"
5850 #define TYPE scm_t_int32
5851 #define TYPE_MIN SCM_T_INT32_MIN
5852 #define TYPE_MAX SCM_T_INT32_MAX
5853 #define SIZEOF_TYPE 4
5854 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5855 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5856 #include "libguile/conv-integer.i.c"
5858 #define TYPE scm_t_uint32
5860 #define TYPE_MAX SCM_T_UINT32_MAX
5861 #define SIZEOF_TYPE 4
5862 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5863 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5864 #include "libguile/conv-uinteger.i.c"
5866 #if SCM_HAVE_T_INT64
5868 #define TYPE scm_t_int64
5869 #define TYPE_MIN SCM_T_INT64_MIN
5870 #define TYPE_MAX SCM_T_INT64_MAX
5871 #define SIZEOF_TYPE 8
5872 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5873 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5874 #include "libguile/conv-integer.i.c"
5876 #define TYPE scm_t_uint64
5878 #define TYPE_MAX SCM_T_UINT64_MAX
5879 #define SIZEOF_TYPE 8
5880 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5881 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5882 #include "libguile/conv-uinteger.i.c"
5887 scm_to_mpz (SCM val
, mpz_t rop
)
5889 if (SCM_I_INUMP (val
))
5890 mpz_set_si (rop
, SCM_I_INUM (val
));
5891 else if (SCM_BIGP (val
))
5892 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5894 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5898 scm_from_mpz (mpz_t val
)
5900 return scm_i_mpz2num (val
);
5904 scm_is_real (SCM val
)
5906 return scm_is_true (scm_real_p (val
));
5910 scm_is_rational (SCM val
)
5912 return scm_is_true (scm_rational_p (val
));
5916 scm_to_double (SCM val
)
5918 if (SCM_I_INUMP (val
))
5919 return SCM_I_INUM (val
);
5920 else if (SCM_BIGP (val
))
5921 return scm_i_big2dbl (val
);
5922 else if (SCM_FRACTIONP (val
))
5923 return scm_i_fraction2double (val
);
5924 else if (SCM_REALP (val
))
5925 return SCM_REAL_VALUE (val
);
5927 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5931 scm_from_double (double val
)
5933 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5934 SCM_REAL_VALUE (z
) = val
;
5938 #if SCM_ENABLE_DISCOURAGED == 1
5941 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5945 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5949 scm_out_of_range (NULL
, num
);
5952 return scm_to_double (num
);
5956 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5960 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5964 scm_out_of_range (NULL
, num
);
5967 return scm_to_double (num
);
5973 scm_is_complex (SCM val
)
5975 return scm_is_true (scm_complex_p (val
));
5979 scm_c_real_part (SCM z
)
5981 if (SCM_COMPLEXP (z
))
5982 return SCM_COMPLEX_REAL (z
);
5985 /* Use the scm_real_part to get proper error checking and
5988 return scm_to_double (scm_real_part (z
));
5993 scm_c_imag_part (SCM z
)
5995 if (SCM_COMPLEXP (z
))
5996 return SCM_COMPLEX_IMAG (z
);
5999 /* Use the scm_imag_part to get proper error checking and
6000 dispatching. The result will almost always be 0.0, but not
6003 return scm_to_double (scm_imag_part (z
));
6008 scm_c_magnitude (SCM z
)
6010 return scm_to_double (scm_magnitude (z
));
6016 return scm_to_double (scm_angle (z
));
6020 scm_is_number (SCM z
)
6022 return scm_is_true (scm_number_p (z
));
6026 /* In the following functions we dispatch to the real-arg funcs like log()
6027 when we know the arg is real, instead of just handing everything to
6028 clog() for instance. This is in case clog() doesn't optimize for a
6029 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6030 well use it to go straight to the applicable C func. */
6032 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6034 "Return the natural logarithm of @var{z}.")
6035 #define FUNC_NAME s_scm_log
6037 if (SCM_COMPLEXP (z
))
6039 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6040 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6042 double re
= SCM_COMPLEX_REAL (z
);
6043 double im
= SCM_COMPLEX_IMAG (z
);
6044 return scm_c_make_rectangular (log (hypot (re
, im
)),
6050 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6051 although the value itself overflows. */
6052 double re
= scm_to_double (z
);
6053 double l
= log (fabs (re
));
6055 return scm_from_double (l
);
6057 return scm_c_make_rectangular (l
, M_PI
);
6063 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6065 "Return the base 10 logarithm of @var{z}.")
6066 #define FUNC_NAME s_scm_log10
6068 if (SCM_COMPLEXP (z
))
6070 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6071 clog() and a multiply by M_LOG10E, rather than the fallback
6072 log10+hypot+atan2.) */
6073 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6074 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6076 double re
= SCM_COMPLEX_REAL (z
);
6077 double im
= SCM_COMPLEX_IMAG (z
);
6078 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6079 M_LOG10E
* atan2 (im
, re
));
6084 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6085 although the value itself overflows. */
6086 double re
= scm_to_double (z
);
6087 double l
= log10 (fabs (re
));
6089 return scm_from_double (l
);
6091 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6097 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6099 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6100 "base of natural logarithms (2.71828@dots{}).")
6101 #define FUNC_NAME s_scm_exp
6103 if (SCM_COMPLEXP (z
))
6105 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6106 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6108 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6109 SCM_COMPLEX_IMAG (z
));
6114 /* When z is a negative bignum the conversion to double overflows,
6115 giving -infinity, but that's ok, the exp is still 0.0. */
6116 return scm_from_double (exp (scm_to_double (z
)));
6122 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6124 "Return the square root of @var{z}. Of the two possible roots\n"
6125 "(positive and negative), the one with the a positive real part\n"
6126 "is returned, or if that's zero then a positive imaginary part.\n"
6130 "(sqrt 9.0) @result{} 3.0\n"
6131 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6132 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6133 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6135 #define FUNC_NAME s_scm_sqrt
6137 if (SCM_COMPLEXP (x
))
6139 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6140 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6142 double re
= SCM_COMPLEX_REAL (x
);
6143 double im
= SCM_COMPLEX_IMAG (x
);
6144 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6145 0.5 * atan2 (im
, re
));
6150 double xx
= scm_to_double (x
);
6152 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6154 return scm_from_double (sqrt (xx
));
6166 mpz_init_set_si (z_negative_one
, -1);
6168 /* It may be possible to tune the performance of some algorithms by using
6169 * the following constants to avoid the creation of bignums. Please, before
6170 * using these values, remember the two rules of program optimization:
6171 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6172 scm_c_define ("most-positive-fixnum",
6173 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6174 scm_c_define ("most-negative-fixnum",
6175 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6177 scm_add_feature ("complex");
6178 scm_add_feature ("inexact");
6179 scm_flo0
= scm_from_double (0.0);
6181 /* determine floating point precision */
6182 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6184 init_dblprec(&scm_dblprec
[i
-2],i
);
6185 init_fx_radix(fx_per_radix
[i
-2],i
);
6188 /* hard code precision for base 10 if the preprocessor tells us to... */
6189 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6192 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6193 SCM_I_MAKINUM (2)));
6194 #include "libguile/numbers.x"