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> */
2686 if (!isdigit ((int) (unsigned char) c
))
2690 exponent
= DIGIT2UINT (c
);
2694 if (isdigit ((int) (unsigned char) c
))
2697 if (exponent
<= SCM_MAXEXP
)
2698 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2704 if (exponent
> SCM_MAXEXP
)
2706 size_t exp_len
= idx
- start
;
2707 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2708 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2709 scm_out_of_range ("string->number", exp_num
);
2712 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2714 result
= scm_product (result
, e
);
2716 result
= scm_divide2real (result
, e
);
2718 /* We've seen an exponent, thus the value is implicitly inexact. */
2736 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2739 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2740 unsigned int radix
, enum t_exactness
*p_exactness
)
2742 unsigned int idx
= *p_idx
;
2745 /* Start off believing that the number will be exact. This changes
2746 to INEXACT if we see a decimal point or a hash. */
2747 enum t_exactness x
= EXACT
;
2752 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2758 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2760 /* Cobble up the fractional part. We might want to set the
2761 NaN's mantissa from it. */
2763 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2768 if (mem
[idx
] == '.')
2772 else if (idx
+ 1 == len
)
2774 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2777 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2784 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2785 if (scm_is_false (uinteger
))
2790 else if (mem
[idx
] == '/')
2798 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2799 if (scm_is_false (divisor
))
2802 /* both are int/big here, I assume */
2803 result
= scm_i_make_ratio (uinteger
, divisor
);
2805 else if (radix
== 10)
2807 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2808 if (scm_is_false (result
))
2817 /* Update *p_exactness if the number just read was inexact. This is
2818 important for complex numbers, so that a complex number is
2819 treated as inexact overall if either its real or imaginary part
2825 /* When returning an inexact zero, make sure it is represented as a
2826 floating point value so that we can change its sign.
2828 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2829 result
= scm_from_double (0.0);
2835 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2838 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2839 unsigned int radix
, enum t_exactness
*p_exactness
)
2863 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2864 if (scm_is_false (ureal
))
2866 /* input must be either +i or -i */
2871 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2877 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2884 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2885 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2894 /* either +<ureal>i or -<ureal>i */
2901 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2904 /* polar input: <real>@<real>. */
2933 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2934 if (scm_is_false (angle
))
2939 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2940 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2942 result
= scm_make_polar (ureal
, angle
);
2947 /* expecting input matching <real>[+-]<ureal>?i */
2954 int sign
= (c
== '+') ? 1 : -1;
2955 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2957 if (scm_is_false (imag
))
2958 imag
= SCM_I_MAKINUM (sign
);
2959 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2960 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2964 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2971 return scm_make_rectangular (ureal
, imag
);
2980 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2982 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2985 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
2986 unsigned int default_radix
)
2988 unsigned int idx
= 0;
2989 unsigned int radix
= NO_RADIX
;
2990 enum t_exactness forced_x
= NO_EXACTNESS
;
2991 enum t_exactness implicit_x
= EXACT
;
2994 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2995 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2997 switch (mem
[idx
+ 1])
3000 if (radix
!= NO_RADIX
)
3005 if (radix
!= NO_RADIX
)
3010 if (forced_x
!= NO_EXACTNESS
)
3015 if (forced_x
!= NO_EXACTNESS
)
3020 if (radix
!= NO_RADIX
)
3025 if (radix
!= NO_RADIX
)
3035 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3036 if (radix
== NO_RADIX
)
3037 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
3039 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
3041 if (scm_is_false (result
))
3047 if (SCM_INEXACTP (result
))
3048 return scm_inexact_to_exact (result
);
3052 if (SCM_INEXACTP (result
))
3055 return scm_exact_to_inexact (result
);
3058 if (implicit_x
== INEXACT
)
3060 if (SCM_INEXACTP (result
))
3063 return scm_exact_to_inexact (result
);
3071 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3072 (SCM string
, SCM radix
),
3073 "Return a number of the maximally precise representation\n"
3074 "expressed by the given @var{string}. @var{radix} must be an\n"
3075 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3076 "is a default radix that may be overridden by an explicit radix\n"
3077 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3078 "supplied, then the default radix is 10. If string is not a\n"
3079 "syntactically valid notation for a number, then\n"
3080 "@code{string->number} returns @code{#f}.")
3081 #define FUNC_NAME s_scm_string_to_number
3085 SCM_VALIDATE_STRING (1, string
);
3087 if (SCM_UNBNDP (radix
))
3090 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3092 answer
= scm_c_locale_stringn_to_number (scm_i_string_chars (string
),
3093 scm_i_string_length (string
),
3095 scm_remember_upto_here_1 (string
);
3101 /*** END strs->nums ***/
3105 scm_bigequal (SCM x
, SCM y
)
3107 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3108 scm_remember_upto_here_2 (x
, y
);
3109 return scm_from_bool (0 == result
);
3113 scm_real_equalp (SCM x
, SCM y
)
3115 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3119 scm_complex_equalp (SCM x
, SCM y
)
3121 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3122 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3126 scm_i_fraction_equalp (SCM x
, SCM y
)
3128 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3129 SCM_FRACTION_NUMERATOR (y
)))
3130 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3131 SCM_FRACTION_DENOMINATOR (y
))))
3138 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3140 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3142 #define FUNC_NAME s_scm_number_p
3144 return scm_from_bool (SCM_NUMBERP (x
));
3148 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3150 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3151 "otherwise. Note that the sets of real, rational and integer\n"
3152 "values form subsets of the set of complex numbers, i. e. the\n"
3153 "predicate will also be fulfilled if @var{x} is a real,\n"
3154 "rational or integer number.")
3155 #define FUNC_NAME s_scm_complex_p
3157 /* all numbers are complex. */
3158 return scm_number_p (x
);
3162 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3164 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3165 "otherwise. Note that the set of integer values forms a subset of\n"
3166 "the set of real numbers, i. e. the predicate will also be\n"
3167 "fulfilled if @var{x} is an integer number.")
3168 #define FUNC_NAME s_scm_real_p
3170 /* we can't represent irrational numbers. */
3171 return scm_rational_p (x
);
3175 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3177 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3178 "otherwise. Note that the set of integer values forms a subset of\n"
3179 "the set of rational numbers, i. e. the predicate will also be\n"
3180 "fulfilled if @var{x} is an integer number.")
3181 #define FUNC_NAME s_scm_rational_p
3183 if (SCM_I_INUMP (x
))
3185 else if (SCM_IMP (x
))
3187 else if (SCM_BIGP (x
))
3189 else if (SCM_FRACTIONP (x
))
3191 else if (SCM_REALP (x
))
3192 /* due to their limited precision, all floating point numbers are
3193 rational as well. */
3200 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3202 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3204 #define FUNC_NAME s_scm_integer_p
3207 if (SCM_I_INUMP (x
))
3213 if (!SCM_INEXACTP (x
))
3215 if (SCM_COMPLEXP (x
))
3217 r
= SCM_REAL_VALUE (x
);
3218 /* +/-inf passes r==floor(r), making those #t */
3226 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3228 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3230 #define FUNC_NAME s_scm_inexact_p
3232 if (SCM_INEXACTP (x
))
3234 if (SCM_NUMBERP (x
))
3236 SCM_WRONG_TYPE_ARG (1, x
);
3241 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3242 /* "Return @code{#t} if all parameters are numerically equal." */
3244 scm_num_eq_p (SCM x
, SCM y
)
3247 if (SCM_I_INUMP (x
))
3249 long xx
= SCM_I_INUM (x
);
3250 if (SCM_I_INUMP (y
))
3252 long yy
= SCM_I_INUM (y
);
3253 return scm_from_bool (xx
== yy
);
3255 else if (SCM_BIGP (y
))
3257 else if (SCM_REALP (y
))
3259 /* On a 32-bit system an inum fits a double, we can cast the inum
3260 to a double and compare.
3262 But on a 64-bit system an inum is bigger than a double and
3263 casting it to a double (call that dxx) will round. dxx is at
3264 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3265 an integer and fits a long. So we cast yy to a long and
3266 compare with plain xx.
3268 An alternative (for any size system actually) would be to check
3269 yy is an integer (with floor) and is in range of an inum
3270 (compare against appropriate powers of 2) then test
3271 xx==(long)yy. It's just a matter of which casts/comparisons
3272 might be fastest or easiest for the cpu. */
3274 double yy
= SCM_REAL_VALUE (y
);
3275 return scm_from_bool ((double) xx
== yy
3276 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3277 || xx
== (long) yy
));
3279 else if (SCM_COMPLEXP (y
))
3280 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3281 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3282 else if (SCM_FRACTIONP (y
))
3285 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3287 else if (SCM_BIGP (x
))
3289 if (SCM_I_INUMP (y
))
3291 else if (SCM_BIGP (y
))
3293 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3294 scm_remember_upto_here_2 (x
, y
);
3295 return scm_from_bool (0 == cmp
);
3297 else if (SCM_REALP (y
))
3300 if (xisnan (SCM_REAL_VALUE (y
)))
3302 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3303 scm_remember_upto_here_1 (x
);
3304 return scm_from_bool (0 == cmp
);
3306 else if (SCM_COMPLEXP (y
))
3309 if (0.0 != SCM_COMPLEX_IMAG (y
))
3311 if (xisnan (SCM_COMPLEX_REAL (y
)))
3313 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3314 scm_remember_upto_here_1 (x
);
3315 return scm_from_bool (0 == cmp
);
3317 else if (SCM_FRACTIONP (y
))
3320 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3322 else if (SCM_REALP (x
))
3324 double xx
= SCM_REAL_VALUE (x
);
3325 if (SCM_I_INUMP (y
))
3327 /* see comments with inum/real above */
3328 long yy
= SCM_I_INUM (y
);
3329 return scm_from_bool (xx
== (double) yy
3330 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3331 || (long) xx
== yy
));
3333 else if (SCM_BIGP (y
))
3336 if (xisnan (SCM_REAL_VALUE (x
)))
3338 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3339 scm_remember_upto_here_1 (y
);
3340 return scm_from_bool (0 == cmp
);
3342 else if (SCM_REALP (y
))
3343 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3344 else if (SCM_COMPLEXP (y
))
3345 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3346 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3347 else if (SCM_FRACTIONP (y
))
3349 double xx
= SCM_REAL_VALUE (x
);
3353 return scm_from_bool (xx
< 0.0);
3354 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3358 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3360 else if (SCM_COMPLEXP (x
))
3362 if (SCM_I_INUMP (y
))
3363 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3364 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3365 else if (SCM_BIGP (y
))
3368 if (0.0 != SCM_COMPLEX_IMAG (x
))
3370 if (xisnan (SCM_COMPLEX_REAL (x
)))
3372 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3373 scm_remember_upto_here_1 (y
);
3374 return scm_from_bool (0 == cmp
);
3376 else if (SCM_REALP (y
))
3377 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3378 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3379 else if (SCM_COMPLEXP (y
))
3380 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3381 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3382 else if (SCM_FRACTIONP (y
))
3385 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3387 xx
= SCM_COMPLEX_REAL (x
);
3391 return scm_from_bool (xx
< 0.0);
3392 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3396 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3398 else if (SCM_FRACTIONP (x
))
3400 if (SCM_I_INUMP (y
))
3402 else if (SCM_BIGP (y
))
3404 else if (SCM_REALP (y
))
3406 double yy
= SCM_REAL_VALUE (y
);
3410 return scm_from_bool (0.0 < yy
);
3411 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3414 else if (SCM_COMPLEXP (y
))
3417 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3419 yy
= SCM_COMPLEX_REAL (y
);
3423 return scm_from_bool (0.0 < yy
);
3424 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3427 else if (SCM_FRACTIONP (y
))
3428 return scm_i_fraction_equalp (x
, y
);
3430 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3433 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3437 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3438 done are good for inums, but for bignums an answer can almost always be
3439 had by just examining a few high bits of the operands, as done by GMP in
3440 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3441 of the float exponent to take into account. */
3443 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3444 /* "Return @code{#t} if the list of parameters is monotonically\n"
3448 scm_less_p (SCM x
, SCM y
)
3451 if (SCM_I_INUMP (x
))
3453 long xx
= SCM_I_INUM (x
);
3454 if (SCM_I_INUMP (y
))
3456 long yy
= SCM_I_INUM (y
);
3457 return scm_from_bool (xx
< yy
);
3459 else if (SCM_BIGP (y
))
3461 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3462 scm_remember_upto_here_1 (y
);
3463 return scm_from_bool (sgn
> 0);
3465 else if (SCM_REALP (y
))
3466 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3467 else if (SCM_FRACTIONP (y
))
3469 /* "x < a/b" becomes "x*b < a" */
3471 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3472 y
= SCM_FRACTION_NUMERATOR (y
);
3476 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3478 else if (SCM_BIGP (x
))
3480 if (SCM_I_INUMP (y
))
3482 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3483 scm_remember_upto_here_1 (x
);
3484 return scm_from_bool (sgn
< 0);
3486 else if (SCM_BIGP (y
))
3488 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3489 scm_remember_upto_here_2 (x
, y
);
3490 return scm_from_bool (cmp
< 0);
3492 else if (SCM_REALP (y
))
3495 if (xisnan (SCM_REAL_VALUE (y
)))
3497 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3498 scm_remember_upto_here_1 (x
);
3499 return scm_from_bool (cmp
< 0);
3501 else if (SCM_FRACTIONP (y
))
3504 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3506 else if (SCM_REALP (x
))
3508 if (SCM_I_INUMP (y
))
3509 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3510 else if (SCM_BIGP (y
))
3513 if (xisnan (SCM_REAL_VALUE (x
)))
3515 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3516 scm_remember_upto_here_1 (y
);
3517 return scm_from_bool (cmp
> 0);
3519 else if (SCM_REALP (y
))
3520 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3521 else if (SCM_FRACTIONP (y
))
3523 double xx
= SCM_REAL_VALUE (x
);
3527 return scm_from_bool (xx
< 0.0);
3528 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3532 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3534 else if (SCM_FRACTIONP (x
))
3536 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3538 /* "a/b < y" becomes "a < y*b" */
3539 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3540 x
= SCM_FRACTION_NUMERATOR (x
);
3543 else if (SCM_REALP (y
))
3545 double yy
= SCM_REAL_VALUE (y
);
3549 return scm_from_bool (0.0 < yy
);
3550 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3553 else if (SCM_FRACTIONP (y
))
3555 /* "a/b < c/d" becomes "a*d < c*b" */
3556 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3557 SCM_FRACTION_DENOMINATOR (y
));
3558 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3559 SCM_FRACTION_DENOMINATOR (x
));
3565 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3568 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3572 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3573 /* "Return @code{#t} if the list of parameters is monotonically\n"
3576 #define FUNC_NAME s_scm_gr_p
3578 scm_gr_p (SCM x
, SCM y
)
3580 if (!SCM_NUMBERP (x
))
3581 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3582 else if (!SCM_NUMBERP (y
))
3583 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3585 return scm_less_p (y
, x
);
3590 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3591 /* "Return @code{#t} if the list of parameters is monotonically\n"
3594 #define FUNC_NAME s_scm_leq_p
3596 scm_leq_p (SCM x
, SCM y
)
3598 if (!SCM_NUMBERP (x
))
3599 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3600 else if (!SCM_NUMBERP (y
))
3601 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3602 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3605 return scm_not (scm_less_p (y
, x
));
3610 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3611 /* "Return @code{#t} if the list of parameters is monotonically\n"
3614 #define FUNC_NAME s_scm_geq_p
3616 scm_geq_p (SCM x
, SCM y
)
3618 if (!SCM_NUMBERP (x
))
3619 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3620 else if (!SCM_NUMBERP (y
))
3621 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3622 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3625 return scm_not (scm_less_p (x
, y
));
3630 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3631 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3637 if (SCM_I_INUMP (z
))
3638 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3639 else if (SCM_BIGP (z
))
3641 else if (SCM_REALP (z
))
3642 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3643 else if (SCM_COMPLEXP (z
))
3644 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3645 && SCM_COMPLEX_IMAG (z
) == 0.0);
3646 else if (SCM_FRACTIONP (z
))
3649 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3653 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3654 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3658 scm_positive_p (SCM x
)
3660 if (SCM_I_INUMP (x
))
3661 return scm_from_bool (SCM_I_INUM (x
) > 0);
3662 else if (SCM_BIGP (x
))
3664 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3665 scm_remember_upto_here_1 (x
);
3666 return scm_from_bool (sgn
> 0);
3668 else if (SCM_REALP (x
))
3669 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3670 else if (SCM_FRACTIONP (x
))
3671 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3673 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3677 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3678 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3682 scm_negative_p (SCM x
)
3684 if (SCM_I_INUMP (x
))
3685 return scm_from_bool (SCM_I_INUM (x
) < 0);
3686 else if (SCM_BIGP (x
))
3688 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3689 scm_remember_upto_here_1 (x
);
3690 return scm_from_bool (sgn
< 0);
3692 else if (SCM_REALP (x
))
3693 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3694 else if (SCM_FRACTIONP (x
))
3695 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3697 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3701 /* scm_min and scm_max return an inexact when either argument is inexact, as
3702 required by r5rs. On that basis, for exact/inexact combinations the
3703 exact is converted to inexact to compare and possibly return. This is
3704 unlike scm_less_p above which takes some trouble to preserve all bits in
3705 its test, such trouble is not required for min and max. */
3707 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3708 /* "Return the maximum of all parameter values."
3711 scm_max (SCM x
, SCM y
)
3716 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3717 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3720 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3723 if (SCM_I_INUMP (x
))
3725 long xx
= SCM_I_INUM (x
);
3726 if (SCM_I_INUMP (y
))
3728 long yy
= SCM_I_INUM (y
);
3729 return (xx
< yy
) ? y
: x
;
3731 else if (SCM_BIGP (y
))
3733 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3734 scm_remember_upto_here_1 (y
);
3735 return (sgn
< 0) ? x
: y
;
3737 else if (SCM_REALP (y
))
3740 /* if y==NaN then ">" is false and we return NaN */
3741 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3743 else if (SCM_FRACTIONP (y
))
3746 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3749 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3751 else if (SCM_BIGP (x
))
3753 if (SCM_I_INUMP (y
))
3755 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3756 scm_remember_upto_here_1 (x
);
3757 return (sgn
< 0) ? y
: x
;
3759 else if (SCM_BIGP (y
))
3761 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3762 scm_remember_upto_here_2 (x
, y
);
3763 return (cmp
> 0) ? x
: y
;
3765 else if (SCM_REALP (y
))
3767 /* if y==NaN then xx>yy is false, so we return the NaN y */
3770 xx
= scm_i_big2dbl (x
);
3771 yy
= SCM_REAL_VALUE (y
);
3772 return (xx
> yy
? scm_from_double (xx
) : y
);
3774 else if (SCM_FRACTIONP (y
))
3779 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3781 else if (SCM_REALP (x
))
3783 if (SCM_I_INUMP (y
))
3785 double z
= SCM_I_INUM (y
);
3786 /* if x==NaN then "<" is false and we return NaN */
3787 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3789 else if (SCM_BIGP (y
))
3794 else if (SCM_REALP (y
))
3796 /* if x==NaN then our explicit check means we return NaN
3797 if y==NaN then ">" is false and we return NaN
3798 calling isnan is unavoidable, since it's the only way to know
3799 which of x or y causes any compares to be false */
3800 double xx
= SCM_REAL_VALUE (x
);
3801 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3803 else if (SCM_FRACTIONP (y
))
3805 double yy
= scm_i_fraction2double (y
);
3806 double xx
= SCM_REAL_VALUE (x
);
3807 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3810 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3812 else if (SCM_FRACTIONP (x
))
3814 if (SCM_I_INUMP (y
))
3818 else if (SCM_BIGP (y
))
3822 else if (SCM_REALP (y
))
3824 double xx
= scm_i_fraction2double (x
);
3825 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3827 else if (SCM_FRACTIONP (y
))
3832 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3835 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3839 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3840 /* "Return the minium of all parameter values."
3843 scm_min (SCM x
, SCM y
)
3848 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3849 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3852 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3855 if (SCM_I_INUMP (x
))
3857 long xx
= SCM_I_INUM (x
);
3858 if (SCM_I_INUMP (y
))
3860 long yy
= SCM_I_INUM (y
);
3861 return (xx
< yy
) ? x
: y
;
3863 else if (SCM_BIGP (y
))
3865 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3866 scm_remember_upto_here_1 (y
);
3867 return (sgn
< 0) ? y
: x
;
3869 else if (SCM_REALP (y
))
3872 /* if y==NaN then "<" is false and we return NaN */
3873 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3875 else if (SCM_FRACTIONP (y
))
3878 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3881 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3883 else if (SCM_BIGP (x
))
3885 if (SCM_I_INUMP (y
))
3887 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3888 scm_remember_upto_here_1 (x
);
3889 return (sgn
< 0) ? x
: y
;
3891 else if (SCM_BIGP (y
))
3893 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3894 scm_remember_upto_here_2 (x
, y
);
3895 return (cmp
> 0) ? y
: x
;
3897 else if (SCM_REALP (y
))
3899 /* if y==NaN then xx<yy is false, so we return the NaN y */
3902 xx
= scm_i_big2dbl (x
);
3903 yy
= SCM_REAL_VALUE (y
);
3904 return (xx
< yy
? scm_from_double (xx
) : y
);
3906 else if (SCM_FRACTIONP (y
))
3911 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3913 else if (SCM_REALP (x
))
3915 if (SCM_I_INUMP (y
))
3917 double z
= SCM_I_INUM (y
);
3918 /* if x==NaN then "<" is false and we return NaN */
3919 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3921 else if (SCM_BIGP (y
))
3926 else if (SCM_REALP (y
))
3928 /* if x==NaN then our explicit check means we return NaN
3929 if y==NaN then "<" is false and we return NaN
3930 calling isnan is unavoidable, since it's the only way to know
3931 which of x or y causes any compares to be false */
3932 double xx
= SCM_REAL_VALUE (x
);
3933 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3935 else if (SCM_FRACTIONP (y
))
3937 double yy
= scm_i_fraction2double (y
);
3938 double xx
= SCM_REAL_VALUE (x
);
3939 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3942 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3944 else if (SCM_FRACTIONP (x
))
3946 if (SCM_I_INUMP (y
))
3950 else if (SCM_BIGP (y
))
3954 else if (SCM_REALP (y
))
3956 double xx
= scm_i_fraction2double (x
);
3957 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3959 else if (SCM_FRACTIONP (y
))
3964 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3967 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3971 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3972 /* "Return the sum of all parameter values. Return 0 if called without\n"
3976 scm_sum (SCM x
, SCM y
)
3978 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
3980 if (SCM_NUMBERP (x
)) return x
;
3981 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3982 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3985 if (SCM_LIKELY (SCM_I_INUMP (x
)))
3987 if (SCM_LIKELY (SCM_I_INUMP (y
)))
3989 long xx
= SCM_I_INUM (x
);
3990 long yy
= SCM_I_INUM (y
);
3991 long int z
= xx
+ yy
;
3992 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3994 else if (SCM_BIGP (y
))
3999 else if (SCM_REALP (y
))
4001 long int xx
= SCM_I_INUM (x
);
4002 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4004 else if (SCM_COMPLEXP (y
))
4006 long int xx
= SCM_I_INUM (x
);
4007 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4008 SCM_COMPLEX_IMAG (y
));
4010 else if (SCM_FRACTIONP (y
))
4011 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4012 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4013 SCM_FRACTION_DENOMINATOR (y
));
4015 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4016 } else if (SCM_BIGP (x
))
4018 if (SCM_I_INUMP (y
))
4023 inum
= SCM_I_INUM (y
);
4026 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4029 SCM result
= scm_i_mkbig ();
4030 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4031 scm_remember_upto_here_1 (x
);
4032 /* we know the result will have to be a bignum */
4035 return scm_i_normbig (result
);
4039 SCM result
= scm_i_mkbig ();
4040 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4041 scm_remember_upto_here_1 (x
);
4042 /* we know the result will have to be a bignum */
4045 return scm_i_normbig (result
);
4048 else if (SCM_BIGP (y
))
4050 SCM result
= scm_i_mkbig ();
4051 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4052 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4053 mpz_add (SCM_I_BIG_MPZ (result
),
4056 scm_remember_upto_here_2 (x
, y
);
4057 /* we know the result will have to be a bignum */
4060 return scm_i_normbig (result
);
4062 else if (SCM_REALP (y
))
4064 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4065 scm_remember_upto_here_1 (x
);
4066 return scm_from_double (result
);
4068 else if (SCM_COMPLEXP (y
))
4070 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4071 + SCM_COMPLEX_REAL (y
));
4072 scm_remember_upto_here_1 (x
);
4073 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4075 else if (SCM_FRACTIONP (y
))
4076 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4077 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4078 SCM_FRACTION_DENOMINATOR (y
));
4080 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4082 else if (SCM_REALP (x
))
4084 if (SCM_I_INUMP (y
))
4085 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4086 else if (SCM_BIGP (y
))
4088 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4089 scm_remember_upto_here_1 (y
);
4090 return scm_from_double (result
);
4092 else if (SCM_REALP (y
))
4093 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4094 else if (SCM_COMPLEXP (y
))
4095 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4096 SCM_COMPLEX_IMAG (y
));
4097 else if (SCM_FRACTIONP (y
))
4098 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4100 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4102 else if (SCM_COMPLEXP (x
))
4104 if (SCM_I_INUMP (y
))
4105 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4106 SCM_COMPLEX_IMAG (x
));
4107 else if (SCM_BIGP (y
))
4109 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4110 + SCM_COMPLEX_REAL (x
));
4111 scm_remember_upto_here_1 (y
);
4112 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4114 else if (SCM_REALP (y
))
4115 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4116 SCM_COMPLEX_IMAG (x
));
4117 else if (SCM_COMPLEXP (y
))
4118 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4119 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4120 else if (SCM_FRACTIONP (y
))
4121 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4122 SCM_COMPLEX_IMAG (x
));
4124 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4126 else if (SCM_FRACTIONP (x
))
4128 if (SCM_I_INUMP (y
))
4129 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4130 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4131 SCM_FRACTION_DENOMINATOR (x
));
4132 else if (SCM_BIGP (y
))
4133 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4134 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4135 SCM_FRACTION_DENOMINATOR (x
));
4136 else if (SCM_REALP (y
))
4137 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4138 else if (SCM_COMPLEXP (y
))
4139 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4140 SCM_COMPLEX_IMAG (y
));
4141 else if (SCM_FRACTIONP (y
))
4142 /* a/b + c/d = (ad + bc) / bd */
4143 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4144 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4145 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4147 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4150 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4154 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4156 "Return @math{@var{x}+1}.")
4157 #define FUNC_NAME s_scm_oneplus
4159 return scm_sum (x
, SCM_I_MAKINUM (1));
4164 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4165 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4166 * the sum of all but the first argument are subtracted from the first
4168 #define FUNC_NAME s_difference
4170 scm_difference (SCM x
, SCM y
)
4172 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4175 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4177 if (SCM_I_INUMP (x
))
4179 long xx
= -SCM_I_INUM (x
);
4180 if (SCM_FIXABLE (xx
))
4181 return SCM_I_MAKINUM (xx
);
4183 return scm_i_long2big (xx
);
4185 else if (SCM_BIGP (x
))
4186 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4187 bignum, but negating that gives a fixnum. */
4188 return scm_i_normbig (scm_i_clonebig (x
, 0));
4189 else if (SCM_REALP (x
))
4190 return scm_from_double (-SCM_REAL_VALUE (x
));
4191 else if (SCM_COMPLEXP (x
))
4192 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4193 -SCM_COMPLEX_IMAG (x
));
4194 else if (SCM_FRACTIONP (x
))
4195 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4196 SCM_FRACTION_DENOMINATOR (x
));
4198 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4201 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4203 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4205 long int xx
= SCM_I_INUM (x
);
4206 long int yy
= SCM_I_INUM (y
);
4207 long int z
= xx
- yy
;
4208 if (SCM_FIXABLE (z
))
4209 return SCM_I_MAKINUM (z
);
4211 return scm_i_long2big (z
);
4213 else if (SCM_BIGP (y
))
4215 /* inum-x - big-y */
4216 long xx
= SCM_I_INUM (x
);
4219 return scm_i_clonebig (y
, 0);
4222 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4223 SCM result
= scm_i_mkbig ();
4226 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4229 /* x - y == -(y + -x) */
4230 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4231 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4233 scm_remember_upto_here_1 (y
);
4235 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4236 /* we know the result will have to be a bignum */
4239 return scm_i_normbig (result
);
4242 else if (SCM_REALP (y
))
4244 long int xx
= SCM_I_INUM (x
);
4245 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4247 else if (SCM_COMPLEXP (y
))
4249 long int xx
= SCM_I_INUM (x
);
4250 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4251 - SCM_COMPLEX_IMAG (y
));
4253 else if (SCM_FRACTIONP (y
))
4254 /* a - b/c = (ac - b) / c */
4255 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4256 SCM_FRACTION_NUMERATOR (y
)),
4257 SCM_FRACTION_DENOMINATOR (y
));
4259 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4261 else if (SCM_BIGP (x
))
4263 if (SCM_I_INUMP (y
))
4265 /* big-x - inum-y */
4266 long yy
= SCM_I_INUM (y
);
4267 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4269 scm_remember_upto_here_1 (x
);
4271 return (SCM_FIXABLE (-yy
) ?
4272 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4275 SCM result
= scm_i_mkbig ();
4278 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4280 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4281 scm_remember_upto_here_1 (x
);
4283 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4284 /* we know the result will have to be a bignum */
4287 return scm_i_normbig (result
);
4290 else if (SCM_BIGP (y
))
4292 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4293 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4294 SCM result
= scm_i_mkbig ();
4295 mpz_sub (SCM_I_BIG_MPZ (result
),
4298 scm_remember_upto_here_2 (x
, y
);
4299 /* we know the result will have to be a bignum */
4300 if ((sgn_x
== 1) && (sgn_y
== -1))
4302 if ((sgn_x
== -1) && (sgn_y
== 1))
4304 return scm_i_normbig (result
);
4306 else if (SCM_REALP (y
))
4308 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4309 scm_remember_upto_here_1 (x
);
4310 return scm_from_double (result
);
4312 else if (SCM_COMPLEXP (y
))
4314 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4315 - SCM_COMPLEX_REAL (y
));
4316 scm_remember_upto_here_1 (x
);
4317 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4319 else if (SCM_FRACTIONP (y
))
4320 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4321 SCM_FRACTION_NUMERATOR (y
)),
4322 SCM_FRACTION_DENOMINATOR (y
));
4323 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4325 else if (SCM_REALP (x
))
4327 if (SCM_I_INUMP (y
))
4328 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4329 else if (SCM_BIGP (y
))
4331 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4332 scm_remember_upto_here_1 (x
);
4333 return scm_from_double (result
);
4335 else if (SCM_REALP (y
))
4336 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4337 else if (SCM_COMPLEXP (y
))
4338 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4339 -SCM_COMPLEX_IMAG (y
));
4340 else if (SCM_FRACTIONP (y
))
4341 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4343 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4345 else if (SCM_COMPLEXP (x
))
4347 if (SCM_I_INUMP (y
))
4348 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4349 SCM_COMPLEX_IMAG (x
));
4350 else if (SCM_BIGP (y
))
4352 double real_part
= (SCM_COMPLEX_REAL (x
)
4353 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4354 scm_remember_upto_here_1 (x
);
4355 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4357 else if (SCM_REALP (y
))
4358 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4359 SCM_COMPLEX_IMAG (x
));
4360 else if (SCM_COMPLEXP (y
))
4361 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4362 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4363 else if (SCM_FRACTIONP (y
))
4364 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4365 SCM_COMPLEX_IMAG (x
));
4367 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4369 else if (SCM_FRACTIONP (x
))
4371 if (SCM_I_INUMP (y
))
4372 /* a/b - c = (a - cb) / b */
4373 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4374 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4375 SCM_FRACTION_DENOMINATOR (x
));
4376 else if (SCM_BIGP (y
))
4377 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4378 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4379 SCM_FRACTION_DENOMINATOR (x
));
4380 else if (SCM_REALP (y
))
4381 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4382 else if (SCM_COMPLEXP (y
))
4383 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4384 -SCM_COMPLEX_IMAG (y
));
4385 else if (SCM_FRACTIONP (y
))
4386 /* a/b - c/d = (ad - bc) / bd */
4387 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4388 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4389 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4391 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4394 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4399 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4401 "Return @math{@var{x}-1}.")
4402 #define FUNC_NAME s_scm_oneminus
4404 return scm_difference (x
, SCM_I_MAKINUM (1));
4409 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4410 /* "Return the product of all arguments. If called without arguments,\n"
4414 scm_product (SCM x
, SCM y
)
4416 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4419 return SCM_I_MAKINUM (1L);
4420 else if (SCM_NUMBERP (x
))
4423 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4426 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4431 xx
= SCM_I_INUM (x
);
4435 case 0: return x
; break;
4436 case 1: return y
; break;
4439 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4441 long yy
= SCM_I_INUM (y
);
4443 SCM k
= SCM_I_MAKINUM (kk
);
4444 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4448 SCM result
= scm_i_long2big (xx
);
4449 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4450 return scm_i_normbig (result
);
4453 else if (SCM_BIGP (y
))
4455 SCM result
= scm_i_mkbig ();
4456 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4457 scm_remember_upto_here_1 (y
);
4460 else if (SCM_REALP (y
))
4461 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4462 else if (SCM_COMPLEXP (y
))
4463 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4464 xx
* SCM_COMPLEX_IMAG (y
));
4465 else if (SCM_FRACTIONP (y
))
4466 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4467 SCM_FRACTION_DENOMINATOR (y
));
4469 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4471 else if (SCM_BIGP (x
))
4473 if (SCM_I_INUMP (y
))
4478 else if (SCM_BIGP (y
))
4480 SCM result
= scm_i_mkbig ();
4481 mpz_mul (SCM_I_BIG_MPZ (result
),
4484 scm_remember_upto_here_2 (x
, y
);
4487 else if (SCM_REALP (y
))
4489 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4490 scm_remember_upto_here_1 (x
);
4491 return scm_from_double (result
);
4493 else if (SCM_COMPLEXP (y
))
4495 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4496 scm_remember_upto_here_1 (x
);
4497 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4498 z
* SCM_COMPLEX_IMAG (y
));
4500 else if (SCM_FRACTIONP (y
))
4501 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4502 SCM_FRACTION_DENOMINATOR (y
));
4504 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4506 else if (SCM_REALP (x
))
4508 if (SCM_I_INUMP (y
))
4510 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4511 if (scm_is_eq (y
, SCM_INUM0
))
4513 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4515 else if (SCM_BIGP (y
))
4517 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4518 scm_remember_upto_here_1 (y
);
4519 return scm_from_double (result
);
4521 else if (SCM_REALP (y
))
4522 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4523 else if (SCM_COMPLEXP (y
))
4524 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4525 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4526 else if (SCM_FRACTIONP (y
))
4527 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4529 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4531 else if (SCM_COMPLEXP (x
))
4533 if (SCM_I_INUMP (y
))
4535 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4536 if (scm_is_eq (y
, SCM_INUM0
))
4538 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4539 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4541 else if (SCM_BIGP (y
))
4543 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4544 scm_remember_upto_here_1 (y
);
4545 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4546 z
* SCM_COMPLEX_IMAG (x
));
4548 else if (SCM_REALP (y
))
4549 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4550 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4551 else if (SCM_COMPLEXP (y
))
4553 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4554 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4555 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4556 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4558 else if (SCM_FRACTIONP (y
))
4560 double yy
= scm_i_fraction2double (y
);
4561 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4562 yy
* SCM_COMPLEX_IMAG (x
));
4565 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4567 else if (SCM_FRACTIONP (x
))
4569 if (SCM_I_INUMP (y
))
4570 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4571 SCM_FRACTION_DENOMINATOR (x
));
4572 else if (SCM_BIGP (y
))
4573 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4574 SCM_FRACTION_DENOMINATOR (x
));
4575 else if (SCM_REALP (y
))
4576 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4577 else if (SCM_COMPLEXP (y
))
4579 double xx
= scm_i_fraction2double (x
);
4580 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4581 xx
* SCM_COMPLEX_IMAG (y
));
4583 else if (SCM_FRACTIONP (y
))
4584 /* a/b * c/d = ac / bd */
4585 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4586 SCM_FRACTION_NUMERATOR (y
)),
4587 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4588 SCM_FRACTION_DENOMINATOR (y
)));
4590 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4593 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4596 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4597 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4598 #define ALLOW_DIVIDE_BY_ZERO
4599 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4602 /* The code below for complex division is adapted from the GNU
4603 libstdc++, which adapted it from f2c's libF77, and is subject to
4606 /****************************************************************
4607 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4609 Permission to use, copy, modify, and distribute this software
4610 and its documentation for any purpose and without fee is hereby
4611 granted, provided that the above copyright notice appear in all
4612 copies and that both that the copyright notice and this
4613 permission notice and warranty disclaimer appear in supporting
4614 documentation, and that the names of AT&T Bell Laboratories or
4615 Bellcore or any of their entities not be used in advertising or
4616 publicity pertaining to distribution of the software without
4617 specific, written prior permission.
4619 AT&T and Bellcore disclaim all warranties with regard to this
4620 software, including all implied warranties of merchantability
4621 and fitness. In no event shall AT&T or Bellcore be liable for
4622 any special, indirect or consequential damages or any damages
4623 whatsoever resulting from loss of use, data or profits, whether
4624 in an action of contract, negligence or other tortious action,
4625 arising out of or in connection with the use or performance of
4627 ****************************************************************/
4629 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4630 /* Divide the first argument by the product of the remaining
4631 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4633 #define FUNC_NAME s_divide
4635 scm_i_divide (SCM x
, SCM y
, int inexact
)
4639 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4642 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4643 else if (SCM_I_INUMP (x
))
4645 long xx
= SCM_I_INUM (x
);
4646 if (xx
== 1 || xx
== -1)
4648 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4650 scm_num_overflow (s_divide
);
4655 return scm_from_double (1.0 / (double) xx
);
4656 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4659 else if (SCM_BIGP (x
))
4662 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4663 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4665 else if (SCM_REALP (x
))
4667 double xx
= SCM_REAL_VALUE (x
);
4668 #ifndef ALLOW_DIVIDE_BY_ZERO
4670 scm_num_overflow (s_divide
);
4673 return scm_from_double (1.0 / xx
);
4675 else if (SCM_COMPLEXP (x
))
4677 double r
= SCM_COMPLEX_REAL (x
);
4678 double i
= SCM_COMPLEX_IMAG (x
);
4679 if (fabs(r
) <= fabs(i
))
4682 double d
= i
* (1.0 + t
* t
);
4683 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4688 double d
= r
* (1.0 + t
* t
);
4689 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4692 else if (SCM_FRACTIONP (x
))
4693 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4694 SCM_FRACTION_NUMERATOR (x
));
4696 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4699 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4701 long xx
= SCM_I_INUM (x
);
4702 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4704 long yy
= SCM_I_INUM (y
);
4707 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4708 scm_num_overflow (s_divide
);
4710 return scm_from_double ((double) xx
/ (double) yy
);
4713 else if (xx
% yy
!= 0)
4716 return scm_from_double ((double) xx
/ (double) yy
);
4717 else return scm_i_make_ratio (x
, y
);
4722 if (SCM_FIXABLE (z
))
4723 return SCM_I_MAKINUM (z
);
4725 return scm_i_long2big (z
);
4728 else if (SCM_BIGP (y
))
4731 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4732 else return scm_i_make_ratio (x
, y
);
4734 else if (SCM_REALP (y
))
4736 double yy
= SCM_REAL_VALUE (y
);
4737 #ifndef ALLOW_DIVIDE_BY_ZERO
4739 scm_num_overflow (s_divide
);
4742 return scm_from_double ((double) xx
/ yy
);
4744 else if (SCM_COMPLEXP (y
))
4747 complex_div
: /* y _must_ be a complex number */
4749 double r
= SCM_COMPLEX_REAL (y
);
4750 double i
= SCM_COMPLEX_IMAG (y
);
4751 if (fabs(r
) <= fabs(i
))
4754 double d
= i
* (1.0 + t
* t
);
4755 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4760 double d
= r
* (1.0 + t
* t
);
4761 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4765 else if (SCM_FRACTIONP (y
))
4766 /* a / b/c = ac / b */
4767 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4768 SCM_FRACTION_NUMERATOR (y
));
4770 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4772 else if (SCM_BIGP (x
))
4774 if (SCM_I_INUMP (y
))
4776 long int yy
= SCM_I_INUM (y
);
4779 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4780 scm_num_overflow (s_divide
);
4782 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4783 scm_remember_upto_here_1 (x
);
4784 return (sgn
== 0) ? scm_nan () : scm_inf ();
4791 /* FIXME: HMM, what are the relative performance issues here?
4792 We need to test. Is it faster on average to test
4793 divisible_p, then perform whichever operation, or is it
4794 faster to perform the integer div opportunistically and
4795 switch to real if there's a remainder? For now we take the
4796 middle ground: test, then if divisible, use the faster div
4799 long abs_yy
= yy
< 0 ? -yy
: yy
;
4800 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4804 SCM result
= scm_i_mkbig ();
4805 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4806 scm_remember_upto_here_1 (x
);
4808 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4809 return scm_i_normbig (result
);
4814 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4815 else return scm_i_make_ratio (x
, y
);
4819 else if (SCM_BIGP (y
))
4821 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4824 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4825 scm_num_overflow (s_divide
);
4827 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4828 scm_remember_upto_here_1 (x
);
4829 return (sgn
== 0) ? scm_nan () : scm_inf ();
4837 /* It's easily possible for the ratio x/y to fit a double
4838 but one or both x and y be too big to fit a double,
4839 hence the use of mpq_get_d rather than converting and
4842 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
4843 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
4844 return scm_from_double (mpq_get_d (q
));
4848 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4852 SCM result
= scm_i_mkbig ();
4853 mpz_divexact (SCM_I_BIG_MPZ (result
),
4856 scm_remember_upto_here_2 (x
, y
);
4857 return scm_i_normbig (result
);
4860 return scm_i_make_ratio (x
, y
);
4864 else if (SCM_REALP (y
))
4866 double yy
= SCM_REAL_VALUE (y
);
4867 #ifndef ALLOW_DIVIDE_BY_ZERO
4869 scm_num_overflow (s_divide
);
4872 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4874 else if (SCM_COMPLEXP (y
))
4876 a
= scm_i_big2dbl (x
);
4879 else if (SCM_FRACTIONP (y
))
4880 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4881 SCM_FRACTION_NUMERATOR (y
));
4883 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4885 else if (SCM_REALP (x
))
4887 double rx
= SCM_REAL_VALUE (x
);
4888 if (SCM_I_INUMP (y
))
4890 long int yy
= SCM_I_INUM (y
);
4891 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4893 scm_num_overflow (s_divide
);
4896 return scm_from_double (rx
/ (double) yy
);
4898 else if (SCM_BIGP (y
))
4900 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4901 scm_remember_upto_here_1 (y
);
4902 return scm_from_double (rx
/ dby
);
4904 else if (SCM_REALP (y
))
4906 double yy
= SCM_REAL_VALUE (y
);
4907 #ifndef ALLOW_DIVIDE_BY_ZERO
4909 scm_num_overflow (s_divide
);
4912 return scm_from_double (rx
/ yy
);
4914 else if (SCM_COMPLEXP (y
))
4919 else if (SCM_FRACTIONP (y
))
4920 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4922 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4924 else if (SCM_COMPLEXP (x
))
4926 double rx
= SCM_COMPLEX_REAL (x
);
4927 double ix
= SCM_COMPLEX_IMAG (x
);
4928 if (SCM_I_INUMP (y
))
4930 long int yy
= SCM_I_INUM (y
);
4931 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4933 scm_num_overflow (s_divide
);
4938 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4941 else if (SCM_BIGP (y
))
4943 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4944 scm_remember_upto_here_1 (y
);
4945 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4947 else if (SCM_REALP (y
))
4949 double yy
= SCM_REAL_VALUE (y
);
4950 #ifndef ALLOW_DIVIDE_BY_ZERO
4952 scm_num_overflow (s_divide
);
4955 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4957 else if (SCM_COMPLEXP (y
))
4959 double ry
= SCM_COMPLEX_REAL (y
);
4960 double iy
= SCM_COMPLEX_IMAG (y
);
4961 if (fabs(ry
) <= fabs(iy
))
4964 double d
= iy
* (1.0 + t
* t
);
4965 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4970 double d
= ry
* (1.0 + t
* t
);
4971 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4974 else if (SCM_FRACTIONP (y
))
4976 double yy
= scm_i_fraction2double (y
);
4977 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4980 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4982 else if (SCM_FRACTIONP (x
))
4984 if (SCM_I_INUMP (y
))
4986 long int yy
= SCM_I_INUM (y
);
4987 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4989 scm_num_overflow (s_divide
);
4992 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4993 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4995 else if (SCM_BIGP (y
))
4997 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4998 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5000 else if (SCM_REALP (y
))
5002 double yy
= SCM_REAL_VALUE (y
);
5003 #ifndef ALLOW_DIVIDE_BY_ZERO
5005 scm_num_overflow (s_divide
);
5008 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5010 else if (SCM_COMPLEXP (y
))
5012 a
= scm_i_fraction2double (x
);
5015 else if (SCM_FRACTIONP (y
))
5016 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5017 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5019 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5022 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5026 scm_divide (SCM x
, SCM y
)
5028 return scm_i_divide (x
, y
, 0);
5031 static SCM
scm_divide2real (SCM x
, SCM y
)
5033 return scm_i_divide (x
, y
, 1);
5039 scm_asinh (double x
)
5044 #define asinh scm_asinh
5045 return log (x
+ sqrt (x
* x
+ 1));
5048 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
5049 /* "Return the inverse hyperbolic sine of @var{x}."
5054 scm_acosh (double x
)
5059 #define acosh scm_acosh
5060 return log (x
+ sqrt (x
* x
- 1));
5063 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
5064 /* "Return the inverse hyperbolic cosine of @var{x}."
5069 scm_atanh (double x
)
5074 #define atanh scm_atanh
5075 return 0.5 * log ((1 + x
) / (1 - x
));
5078 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
5079 /* "Return the inverse hyperbolic tangent of @var{x}."
5084 scm_c_truncate (double x
)
5095 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5096 half-way case (ie. when x is an integer plus 0.5) going upwards.
5097 Then half-way cases are identified and adjusted down if the
5098 round-upwards didn't give the desired even integer.
5100 "plus_half == result" identifies a half-way case. If plus_half, which is
5101 x + 0.5, is an integer then x must be an integer plus 0.5.
5103 An odd "result" value is identified with result/2 != floor(result/2).
5104 This is done with plus_half, since that value is ready for use sooner in
5105 a pipelined cpu, and we're already requiring plus_half == result.
5107 Note however that we need to be careful when x is big and already an
5108 integer. In that case "x+0.5" may round to an adjacent integer, causing
5109 us to return such a value, incorrectly. For instance if the hardware is
5110 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5111 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5112 returned. Or if the hardware is in round-upwards mode, then other bigger
5113 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5114 representable value, 2^128+2^76 (or whatever), again incorrect.
5116 These bad roundings of x+0.5 are avoided by testing at the start whether
5117 x is already an integer. If it is then clearly that's the desired result
5118 already. And if it's not then the exponent must be small enough to allow
5119 an 0.5 to be represented, and hence added without a bad rounding. */
5122 scm_c_round (double x
)
5124 double plus_half
, result
;
5129 plus_half
= x
+ 0.5;
5130 result
= floor (plus_half
);
5131 /* Adjust so that the rounding is towards even. */
5132 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5137 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5139 "Round the number @var{x} towards zero.")
5140 #define FUNC_NAME s_scm_truncate_number
5142 if (scm_is_false (scm_negative_p (x
)))
5143 return scm_floor (x
);
5145 return scm_ceiling (x
);
5149 static SCM exactly_one_half
;
5151 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5153 "Round the number @var{x} towards the nearest integer. "
5154 "When it is exactly halfway between two integers, "
5155 "round towards the even one.")
5156 #define FUNC_NAME s_scm_round_number
5158 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5160 else if (SCM_REALP (x
))
5161 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5164 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5165 single quotient+remainder division then examining to see which way
5166 the rounding should go. */
5167 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5168 SCM result
= scm_floor (plus_half
);
5169 /* Adjust so that the rounding is towards even. */
5170 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5171 && scm_is_true (scm_odd_p (result
)))
5172 return scm_difference (result
, SCM_I_MAKINUM (1));
5179 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5181 "Round the number @var{x} towards minus infinity.")
5182 #define FUNC_NAME s_scm_floor
5184 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5186 else if (SCM_REALP (x
))
5187 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5188 else if (SCM_FRACTIONP (x
))
5190 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5191 SCM_FRACTION_DENOMINATOR (x
));
5192 if (scm_is_false (scm_negative_p (x
)))
5194 /* For positive x, rounding towards zero is correct. */
5199 /* For negative x, we need to return q-1 unless x is an
5200 integer. But fractions are never integer, per our
5202 return scm_difference (q
, SCM_I_MAKINUM (1));
5206 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5210 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5212 "Round the number @var{x} towards infinity.")
5213 #define FUNC_NAME s_scm_ceiling
5215 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5217 else if (SCM_REALP (x
))
5218 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5219 else if (SCM_FRACTIONP (x
))
5221 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5222 SCM_FRACTION_DENOMINATOR (x
));
5223 if (scm_is_false (scm_positive_p (x
)))
5225 /* For negative x, rounding towards zero is correct. */
5230 /* For positive x, we need to return q+1 unless x is an
5231 integer. But fractions are never integer, per our
5233 return scm_sum (q
, SCM_I_MAKINUM (1));
5237 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5241 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5242 /* "Return the square root of the real number @var{x}."
5244 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5245 /* "Return the absolute value of the real number @var{x}."
5247 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5248 /* "Return the @var{x}th power of e."
5250 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5251 /* "Return the natural logarithm of the real number @var{x}."
5253 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5254 /* "Return the sine of the real number @var{x}."
5256 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5257 /* "Return the cosine of the real number @var{x}."
5259 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5260 /* "Return the tangent of the real number @var{x}."
5262 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5263 /* "Return the arc sine of the real number @var{x}."
5265 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5266 /* "Return the arc cosine of the real number @var{x}."
5268 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5269 /* "Return the arc tangent of the real number @var{x}."
5271 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5272 /* "Return the hyperbolic sine of the real number @var{x}."
5274 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5275 /* "Return the hyperbolic cosine of the real number @var{x}."
5277 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5278 /* "Return the hyperbolic tangent of the real number @var{x}."
5286 static void scm_two_doubles (SCM x
,
5288 const char *sstring
,
5292 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5294 if (SCM_I_INUMP (x
))
5295 xy
->x
= SCM_I_INUM (x
);
5296 else if (SCM_BIGP (x
))
5297 xy
->x
= scm_i_big2dbl (x
);
5298 else if (SCM_REALP (x
))
5299 xy
->x
= SCM_REAL_VALUE (x
);
5300 else if (SCM_FRACTIONP (x
))
5301 xy
->x
= scm_i_fraction2double (x
);
5303 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5305 if (SCM_I_INUMP (y
))
5306 xy
->y
= SCM_I_INUM (y
);
5307 else if (SCM_BIGP (y
))
5308 xy
->y
= scm_i_big2dbl (y
);
5309 else if (SCM_REALP (y
))
5310 xy
->y
= SCM_REAL_VALUE (y
);
5311 else if (SCM_FRACTIONP (y
))
5312 xy
->y
= scm_i_fraction2double (y
);
5314 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5318 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5320 "Return @var{x} raised to the power of @var{y}. This\n"
5321 "procedure does not accept complex arguments.")
5322 #define FUNC_NAME s_scm_sys_expt
5325 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5326 return scm_from_double (pow (xy
.x
, xy
.y
));
5331 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5333 "Return the arc tangent of the two arguments @var{x} and\n"
5334 "@var{y}. This is similar to calculating the arc tangent of\n"
5335 "@var{x} / @var{y}, except that the signs of both arguments\n"
5336 "are used to determine the quadrant of the result. This\n"
5337 "procedure does not accept complex arguments.")
5338 #define FUNC_NAME s_scm_sys_atan2
5341 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5342 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5347 scm_c_make_rectangular (double re
, double im
)
5350 return scm_from_double (re
);
5354 SCM_NEWSMOB (z
, scm_tc16_complex
,
5355 scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5357 SCM_COMPLEX_REAL (z
) = re
;
5358 SCM_COMPLEX_IMAG (z
) = im
;
5363 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5364 (SCM real_part
, SCM imaginary_part
),
5365 "Return a complex number constructed of the given @var{real-part} "
5366 "and @var{imaginary-part} parts.")
5367 #define FUNC_NAME s_scm_make_rectangular
5370 scm_two_doubles (real_part
, imaginary_part
, FUNC_NAME
, &xy
);
5371 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5376 scm_c_make_polar (double mag
, double ang
)
5380 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5381 use it on Glibc-based systems that have it (it's a GNU extension). See
5382 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5384 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5385 sincos (ang
, &s
, &c
);
5390 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5393 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5395 "Return the complex number @var{x} * e^(i * @var{y}).")
5396 #define FUNC_NAME s_scm_make_polar
5399 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5400 return scm_c_make_polar (xy
.x
, xy
.y
);
5405 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5406 /* "Return the real part of the number @var{z}."
5409 scm_real_part (SCM z
)
5411 if (SCM_I_INUMP (z
))
5413 else if (SCM_BIGP (z
))
5415 else if (SCM_REALP (z
))
5417 else if (SCM_COMPLEXP (z
))
5418 return scm_from_double (SCM_COMPLEX_REAL (z
));
5419 else if (SCM_FRACTIONP (z
))
5422 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5426 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5427 /* "Return the imaginary part of the number @var{z}."
5430 scm_imag_part (SCM z
)
5432 if (SCM_I_INUMP (z
))
5434 else if (SCM_BIGP (z
))
5436 else if (SCM_REALP (z
))
5438 else if (SCM_COMPLEXP (z
))
5439 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5440 else if (SCM_FRACTIONP (z
))
5443 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5446 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5447 /* "Return the numerator of the number @var{z}."
5450 scm_numerator (SCM z
)
5452 if (SCM_I_INUMP (z
))
5454 else if (SCM_BIGP (z
))
5456 else if (SCM_FRACTIONP (z
))
5457 return SCM_FRACTION_NUMERATOR (z
);
5458 else if (SCM_REALP (z
))
5459 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5461 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5465 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5466 /* "Return the denominator of the number @var{z}."
5469 scm_denominator (SCM z
)
5471 if (SCM_I_INUMP (z
))
5472 return SCM_I_MAKINUM (1);
5473 else if (SCM_BIGP (z
))
5474 return SCM_I_MAKINUM (1);
5475 else if (SCM_FRACTIONP (z
))
5476 return SCM_FRACTION_DENOMINATOR (z
);
5477 else if (SCM_REALP (z
))
5478 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5480 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5483 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5484 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5485 * "@code{abs} for real arguments, but also allows complex numbers."
5488 scm_magnitude (SCM z
)
5490 if (SCM_I_INUMP (z
))
5492 long int zz
= SCM_I_INUM (z
);
5495 else if (SCM_POSFIXABLE (-zz
))
5496 return SCM_I_MAKINUM (-zz
);
5498 return scm_i_long2big (-zz
);
5500 else if (SCM_BIGP (z
))
5502 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5503 scm_remember_upto_here_1 (z
);
5505 return scm_i_clonebig (z
, 0);
5509 else if (SCM_REALP (z
))
5510 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5511 else if (SCM_COMPLEXP (z
))
5512 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5513 else if (SCM_FRACTIONP (z
))
5515 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5517 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5518 SCM_FRACTION_DENOMINATOR (z
));
5521 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5525 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5526 /* "Return the angle of the complex number @var{z}."
5531 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5532 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5533 But if atan2 follows the floating point rounding mode, then the value
5534 is not a constant. Maybe it'd be close enough though. */
5535 if (SCM_I_INUMP (z
))
5537 if (SCM_I_INUM (z
) >= 0)
5540 return scm_from_double (atan2 (0.0, -1.0));
5542 else if (SCM_BIGP (z
))
5544 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5545 scm_remember_upto_here_1 (z
);
5547 return scm_from_double (atan2 (0.0, -1.0));
5551 else if (SCM_REALP (z
))
5553 if (SCM_REAL_VALUE (z
) >= 0)
5556 return scm_from_double (atan2 (0.0, -1.0));
5558 else if (SCM_COMPLEXP (z
))
5559 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5560 else if (SCM_FRACTIONP (z
))
5562 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5564 else return scm_from_double (atan2 (0.0, -1.0));
5567 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5571 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5572 /* Convert the number @var{x} to its inexact representation.\n"
5575 scm_exact_to_inexact (SCM z
)
5577 if (SCM_I_INUMP (z
))
5578 return scm_from_double ((double) SCM_I_INUM (z
));
5579 else if (SCM_BIGP (z
))
5580 return scm_from_double (scm_i_big2dbl (z
));
5581 else if (SCM_FRACTIONP (z
))
5582 return scm_from_double (scm_i_fraction2double (z
));
5583 else if (SCM_INEXACTP (z
))
5586 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5590 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5592 "Return an exact number that is numerically closest to @var{z}.")
5593 #define FUNC_NAME s_scm_inexact_to_exact
5595 if (SCM_I_INUMP (z
))
5597 else if (SCM_BIGP (z
))
5599 else if (SCM_REALP (z
))
5601 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5602 SCM_OUT_OF_RANGE (1, z
);
5609 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5610 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5611 scm_i_mpz2num (mpq_denref (frac
)));
5613 /* When scm_i_make_ratio throws, we leak the memory allocated
5620 else if (SCM_FRACTIONP (z
))
5623 SCM_WRONG_TYPE_ARG (1, z
);
5627 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5629 "Returns the @emph{simplest} rational number differing\n"
5630 "from @var{x} by no more than @var{eps}.\n"
5632 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
5633 "exact result when both its arguments are exact. Thus, you might need\n"
5634 "to use @code{inexact->exact} on the arguments.\n"
5637 "(rationalize (inexact->exact 1.2) 1/100)\n"
5640 #define FUNC_NAME s_scm_rationalize
5642 if (SCM_I_INUMP (x
))
5644 else if (SCM_BIGP (x
))
5646 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5648 /* Use continued fractions to find closest ratio. All
5649 arithmetic is done with exact numbers.
5652 SCM ex
= scm_inexact_to_exact (x
);
5653 SCM int_part
= scm_floor (ex
);
5654 SCM tt
= SCM_I_MAKINUM (1);
5655 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5656 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5660 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5663 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5664 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5666 /* We stop after a million iterations just to be absolutely sure
5667 that we don't go into an infinite loop. The process normally
5668 converges after less than a dozen iterations.
5671 eps
= scm_abs (eps
);
5672 while (++i
< 1000000)
5674 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5675 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5676 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5678 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5679 eps
))) /* abs(x-a/b) <= eps */
5681 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5682 if (scm_is_false (scm_exact_p (x
))
5683 || scm_is_false (scm_exact_p (eps
)))
5684 return scm_exact_to_inexact (res
);
5688 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5690 tt
= scm_floor (rx
); /* tt = floor (rx) */
5696 scm_num_overflow (s_scm_rationalize
);
5699 SCM_WRONG_TYPE_ARG (1, x
);
5703 /* conversion functions */
5706 scm_is_integer (SCM val
)
5708 return scm_is_true (scm_integer_p (val
));
5712 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5714 if (SCM_I_INUMP (val
))
5716 scm_t_signed_bits n
= SCM_I_INUM (val
);
5717 return n
>= min
&& n
<= max
;
5719 else if (SCM_BIGP (val
))
5721 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5723 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5725 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5727 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5728 return n
>= min
&& n
<= max
;
5738 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5739 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5742 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5743 SCM_I_BIG_MPZ (val
));
5745 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5757 return n
>= min
&& n
<= max
;
5765 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5767 if (SCM_I_INUMP (val
))
5769 scm_t_signed_bits n
= SCM_I_INUM (val
);
5770 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5772 else if (SCM_BIGP (val
))
5774 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5776 else if (max
<= ULONG_MAX
)
5778 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5780 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5781 return n
>= min
&& n
<= max
;
5791 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5794 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5795 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5798 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5799 SCM_I_BIG_MPZ (val
));
5801 return n
>= min
&& n
<= max
;
5809 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5811 scm_error (scm_out_of_range_key
,
5813 "Value out of range ~S to ~S: ~S",
5814 scm_list_3 (min
, max
, bad_val
),
5815 scm_list_1 (bad_val
));
5818 #define TYPE scm_t_intmax
5819 #define TYPE_MIN min
5820 #define TYPE_MAX max
5821 #define SIZEOF_TYPE 0
5822 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5823 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5824 #include "libguile/conv-integer.i.c"
5826 #define TYPE scm_t_uintmax
5827 #define TYPE_MIN min
5828 #define TYPE_MAX max
5829 #define SIZEOF_TYPE 0
5830 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5831 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5832 #include "libguile/conv-uinteger.i.c"
5834 #define TYPE scm_t_int8
5835 #define TYPE_MIN SCM_T_INT8_MIN
5836 #define TYPE_MAX SCM_T_INT8_MAX
5837 #define SIZEOF_TYPE 1
5838 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5839 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5840 #include "libguile/conv-integer.i.c"
5842 #define TYPE scm_t_uint8
5844 #define TYPE_MAX SCM_T_UINT8_MAX
5845 #define SIZEOF_TYPE 1
5846 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5847 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5848 #include "libguile/conv-uinteger.i.c"
5850 #define TYPE scm_t_int16
5851 #define TYPE_MIN SCM_T_INT16_MIN
5852 #define TYPE_MAX SCM_T_INT16_MAX
5853 #define SIZEOF_TYPE 2
5854 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5855 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5856 #include "libguile/conv-integer.i.c"
5858 #define TYPE scm_t_uint16
5860 #define TYPE_MAX SCM_T_UINT16_MAX
5861 #define SIZEOF_TYPE 2
5862 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5863 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5864 #include "libguile/conv-uinteger.i.c"
5866 #define TYPE scm_t_int32
5867 #define TYPE_MIN SCM_T_INT32_MIN
5868 #define TYPE_MAX SCM_T_INT32_MAX
5869 #define SIZEOF_TYPE 4
5870 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5871 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5872 #include "libguile/conv-integer.i.c"
5874 #define TYPE scm_t_uint32
5876 #define TYPE_MAX SCM_T_UINT32_MAX
5877 #define SIZEOF_TYPE 4
5878 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5879 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5880 #include "libguile/conv-uinteger.i.c"
5882 #define TYPE scm_t_wchar
5883 #define TYPE_MIN (scm_t_int32)-1
5884 #define TYPE_MAX (scm_t_int32)0x10ffff
5885 #define SIZEOF_TYPE 4
5886 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
5887 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
5888 #include "libguile/conv-integer.i.c"
5890 #if SCM_HAVE_T_INT64
5892 #define TYPE scm_t_int64
5893 #define TYPE_MIN SCM_T_INT64_MIN
5894 #define TYPE_MAX SCM_T_INT64_MAX
5895 #define SIZEOF_TYPE 8
5896 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5897 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5898 #include "libguile/conv-integer.i.c"
5900 #define TYPE scm_t_uint64
5902 #define TYPE_MAX SCM_T_UINT64_MAX
5903 #define SIZEOF_TYPE 8
5904 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5905 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5906 #include "libguile/conv-uinteger.i.c"
5911 scm_to_mpz (SCM val
, mpz_t rop
)
5913 if (SCM_I_INUMP (val
))
5914 mpz_set_si (rop
, SCM_I_INUM (val
));
5915 else if (SCM_BIGP (val
))
5916 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5918 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5922 scm_from_mpz (mpz_t val
)
5924 return scm_i_mpz2num (val
);
5928 scm_is_real (SCM val
)
5930 return scm_is_true (scm_real_p (val
));
5934 scm_is_rational (SCM val
)
5936 return scm_is_true (scm_rational_p (val
));
5940 scm_to_double (SCM val
)
5942 if (SCM_I_INUMP (val
))
5943 return SCM_I_INUM (val
);
5944 else if (SCM_BIGP (val
))
5945 return scm_i_big2dbl (val
);
5946 else if (SCM_FRACTIONP (val
))
5947 return scm_i_fraction2double (val
);
5948 else if (SCM_REALP (val
))
5949 return SCM_REAL_VALUE (val
);
5951 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5955 scm_from_double (double val
)
5957 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5958 SCM_REAL_VALUE (z
) = val
;
5962 #if SCM_ENABLE_DISCOURAGED == 1
5965 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5969 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5973 scm_out_of_range (NULL
, num
);
5976 return scm_to_double (num
);
5980 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5984 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5988 scm_out_of_range (NULL
, num
);
5991 return scm_to_double (num
);
5997 scm_is_complex (SCM val
)
5999 return scm_is_true (scm_complex_p (val
));
6003 scm_c_real_part (SCM z
)
6005 if (SCM_COMPLEXP (z
))
6006 return SCM_COMPLEX_REAL (z
);
6009 /* Use the scm_real_part to get proper error checking and
6012 return scm_to_double (scm_real_part (z
));
6017 scm_c_imag_part (SCM z
)
6019 if (SCM_COMPLEXP (z
))
6020 return SCM_COMPLEX_IMAG (z
);
6023 /* Use the scm_imag_part to get proper error checking and
6024 dispatching. The result will almost always be 0.0, but not
6027 return scm_to_double (scm_imag_part (z
));
6032 scm_c_magnitude (SCM z
)
6034 return scm_to_double (scm_magnitude (z
));
6040 return scm_to_double (scm_angle (z
));
6044 scm_is_number (SCM z
)
6046 return scm_is_true (scm_number_p (z
));
6050 /* In the following functions we dispatch to the real-arg funcs like log()
6051 when we know the arg is real, instead of just handing everything to
6052 clog() for instance. This is in case clog() doesn't optimize for a
6053 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6054 well use it to go straight to the applicable C func. */
6056 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6058 "Return the natural logarithm of @var{z}.")
6059 #define FUNC_NAME s_scm_log
6061 if (SCM_COMPLEXP (z
))
6063 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6064 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6066 double re
= SCM_COMPLEX_REAL (z
);
6067 double im
= SCM_COMPLEX_IMAG (z
);
6068 return scm_c_make_rectangular (log (hypot (re
, im
)),
6074 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6075 although the value itself overflows. */
6076 double re
= scm_to_double (z
);
6077 double l
= log (fabs (re
));
6079 return scm_from_double (l
);
6081 return scm_c_make_rectangular (l
, M_PI
);
6087 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6089 "Return the base 10 logarithm of @var{z}.")
6090 #define FUNC_NAME s_scm_log10
6092 if (SCM_COMPLEXP (z
))
6094 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6095 clog() and a multiply by M_LOG10E, rather than the fallback
6096 log10+hypot+atan2.) */
6097 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6098 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6100 double re
= SCM_COMPLEX_REAL (z
);
6101 double im
= SCM_COMPLEX_IMAG (z
);
6102 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6103 M_LOG10E
* atan2 (im
, re
));
6108 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6109 although the value itself overflows. */
6110 double re
= scm_to_double (z
);
6111 double l
= log10 (fabs (re
));
6113 return scm_from_double (l
);
6115 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6121 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6123 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6124 "base of natural logarithms (2.71828@dots{}).")
6125 #define FUNC_NAME s_scm_exp
6127 if (SCM_COMPLEXP (z
))
6129 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6130 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6132 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6133 SCM_COMPLEX_IMAG (z
));
6138 /* When z is a negative bignum the conversion to double overflows,
6139 giving -infinity, but that's ok, the exp is still 0.0. */
6140 return scm_from_double (exp (scm_to_double (z
)));
6146 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6148 "Return the square root of @var{z}. Of the two possible roots\n"
6149 "(positive and negative), the one with the a positive real part\n"
6150 "is returned, or if that's zero then a positive imaginary part.\n"
6154 "(sqrt 9.0) @result{} 3.0\n"
6155 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6156 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6157 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6159 #define FUNC_NAME s_scm_sqrt
6161 if (SCM_COMPLEXP (x
))
6163 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6164 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6166 double re
= SCM_COMPLEX_REAL (x
);
6167 double im
= SCM_COMPLEX_IMAG (x
);
6168 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6169 0.5 * atan2 (im
, re
));
6174 double xx
= scm_to_double (x
);
6176 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6178 return scm_from_double (sqrt (xx
));
6190 mpz_init_set_si (z_negative_one
, -1);
6192 /* It may be possible to tune the performance of some algorithms by using
6193 * the following constants to avoid the creation of bignums. Please, before
6194 * using these values, remember the two rules of program optimization:
6195 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6196 scm_c_define ("most-positive-fixnum",
6197 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6198 scm_c_define ("most-negative-fixnum",
6199 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6201 scm_add_feature ("complex");
6202 scm_add_feature ("inexact");
6203 scm_flo0
= scm_from_double (0.0);
6205 /* determine floating point precision */
6206 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6208 init_dblprec(&scm_dblprec
[i
-2],i
);
6209 init_fx_radix(fx_per_radix
[i
-2],i
);
6212 /* hard code precision for base 10 if the preprocessor tells us to... */
6213 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6216 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6217 SCM_I_MAKINUM (2)));
6218 #include "libguile/numbers.x"