1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008 Free Software Foundation, Inc.
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 /* General assumptions:
24 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
25 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
26 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
27 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
28 * All objects satisfying SCM_FRACTIONP are never an integer.
33 - see if special casing bignums and reals in integer-exponent when
34 possible (to use mpz_pow and mpf_pow_ui) is faster.
36 - look in to better short-circuiting of common cases in
37 integer-expt and elsewhere.
39 - see if direct mpz operations can help in ash and elsewhere.
43 /* tell glibc (2.3) to give prototype for C99 trunc(), csqrt(), etc */
58 #include "libguile/_scm.h"
59 #include "libguile/feature.h"
60 #include "libguile/ports.h"
61 #include "libguile/root.h"
62 #include "libguile/smob.h"
63 #include "libguile/strings.h"
65 #include "libguile/validate.h"
66 #include "libguile/numbers.h"
67 #include "libguile/deprecation.h"
69 #include "libguile/eq.h"
71 #include "libguile/discouraged.h"
73 /* values per glibc, if not already defined */
75 #define M_LOG10E 0.43429448190325182765
78 #define M_PI 3.14159265358979323846
84 Wonder if this might be faster for some of our code? A switch on
85 the numtag would jump directly to the right case, and the
86 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
88 #define SCM_I_NUMTAG_NOTNUM 0
89 #define SCM_I_NUMTAG_INUM 1
90 #define SCM_I_NUMTAG_BIG scm_tc16_big
91 #define SCM_I_NUMTAG_REAL scm_tc16_real
92 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
93 #define SCM_I_NUMTAG(x) \
94 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
95 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
96 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
97 : SCM_I_NUMTAG_NOTNUM)))
99 /* the macro above will not work as is with fractions */
102 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
104 /* FLOBUFLEN is the maximum number of characters neccessary for the
105 * printed or scm_string representation of an inexact number.
107 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
110 #if ! defined (HAVE_ISNAN)
115 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
118 #if ! defined (HAVE_ISINF)
123 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
130 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
131 an explicit check. In some future gmp (don't know what version number),
132 mpz_cmp_d is supposed to do this itself. */
134 #define xmpz_cmp_d(z, d) \
135 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
137 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
140 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
141 isinf. It does have finite and isnan though, hence the use of those.
142 fpclass would be a possibility on that system too. */
146 #if defined (HAVE_ISINF)
148 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
149 return (! (finite (x
) || isnan (x
)));
158 #if defined (HAVE_ISNAN)
165 #if defined (GUILE_I)
166 #if HAVE_COMPLEX_DOUBLE
168 /* For an SCM object Z which is a complex number (ie. satisfies
169 SCM_COMPLEXP), return its value as a C level "complex double". */
170 #define SCM_COMPLEX_VALUE(z) \
171 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
173 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
175 /* Convert a C "complex double" to an SCM value. */
177 scm_from_complex_double (complex double z
)
179 return scm_c_make_rectangular (creal (z
), cimag (z
));
182 #endif /* HAVE_COMPLEX_DOUBLE */
187 static mpz_t z_negative_one
;
194 /* Return a newly created bignum. */
195 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
196 mpz_init (SCM_I_BIG_MPZ (z
));
201 scm_i_long2big (long x
)
203 /* Return a newly created bignum initialized to X. */
204 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
205 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
210 scm_i_ulong2big (unsigned long x
)
212 /* Return a newly created bignum initialized to X. */
213 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
214 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
219 scm_i_clonebig (SCM src_big
, int same_sign_p
)
221 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
222 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
223 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
225 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
230 scm_i_bigcmp (SCM x
, SCM y
)
232 /* Return neg if x < y, pos if x > y, and 0 if x == y */
233 /* presume we already know x and y are bignums */
234 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
235 scm_remember_upto_here_2 (x
, y
);
240 scm_i_dbl2big (double d
)
242 /* results are only defined if d is an integer */
243 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
244 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
248 /* Convert a integer in double representation to a SCM number. */
251 scm_i_dbl2num (double u
)
253 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
254 powers of 2, so there's no rounding when making "double" values
255 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
256 get rounded on a 64-bit machine, hence the "+1".
258 The use of floor() to force to an integer value ensures we get a
259 "numerically closest" value without depending on how a
260 double->long cast or how mpz_set_d will round. For reference,
261 double->long probably follows the hardware rounding mode,
262 mpz_set_d truncates towards zero. */
264 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
265 representable as a double? */
267 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
268 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
269 return SCM_I_MAKINUM ((long) u
);
271 return scm_i_dbl2big (u
);
274 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
275 with R5RS exact->inexact.
277 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
278 (ie. truncate towards zero), then adjust to get the closest double by
279 examining the next lower bit and adding 1 (to the absolute value) if
282 Bignums exactly half way between representable doubles are rounded to the
283 next higher absolute value (ie. away from zero). This seems like an
284 adequate interpretation of R5RS "numerically closest", and it's easier
285 and faster than a full "nearest-even" style.
287 The bit test must be done on the absolute value of the mpz_t, which means
288 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
289 negatives as twos complement.
291 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
292 following the hardware rounding mode, but applied to the absolute value
293 of the mpz_t operand. This is not what we want so we put the high
294 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
295 mpz_get_d is supposed to always truncate towards zero.
297 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
298 is a slowdown. It'd be faster to pick out the relevant high bits with
299 mpz_getlimbn if we could be bothered coding that, and if the new
300 truncating gmp doesn't come out. */
303 scm_i_big2dbl (SCM b
)
308 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
312 /* Current GMP, eg. 4.1.3, force truncation towards zero */
314 if (bits
> DBL_MANT_DIG
)
316 size_t shift
= bits
- DBL_MANT_DIG
;
317 mpz_init2 (tmp
, DBL_MANT_DIG
);
318 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
319 result
= ldexp (mpz_get_d (tmp
), shift
);
324 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
329 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
332 if (bits
> DBL_MANT_DIG
)
334 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
335 /* test bit number "pos" in absolute value */
336 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
337 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
339 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
343 scm_remember_upto_here_1 (b
);
348 scm_i_normbig (SCM b
)
350 /* convert a big back to a fixnum if it'll fit */
351 /* presume b is a bignum */
352 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
354 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
355 if (SCM_FIXABLE (val
))
356 b
= SCM_I_MAKINUM (val
);
361 static SCM_C_INLINE_KEYWORD SCM
362 scm_i_mpz2num (mpz_t b
)
364 /* convert a mpz number to a SCM number. */
365 if (mpz_fits_slong_p (b
))
367 long val
= mpz_get_si (b
);
368 if (SCM_FIXABLE (val
))
369 return SCM_I_MAKINUM (val
);
373 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
374 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
379 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
380 static SCM
scm_divide2real (SCM x
, SCM y
);
383 scm_i_make_ratio (SCM numerator
, SCM denominator
)
384 #define FUNC_NAME "make-ratio"
386 /* First make sure the arguments are proper.
388 if (SCM_I_INUMP (denominator
))
390 if (scm_is_eq (denominator
, SCM_INUM0
))
391 scm_num_overflow ("make-ratio");
392 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
397 if (!(SCM_BIGP(denominator
)))
398 SCM_WRONG_TYPE_ARG (2, denominator
);
400 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
401 SCM_WRONG_TYPE_ARG (1, numerator
);
403 /* Then flip signs so that the denominator is positive.
405 if (scm_is_true (scm_negative_p (denominator
)))
407 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
408 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
411 /* Now consider for each of the four fixnum/bignum combinations
412 whether the rational number is really an integer.
414 if (SCM_I_INUMP (numerator
))
416 long x
= SCM_I_INUM (numerator
);
417 if (scm_is_eq (numerator
, SCM_INUM0
))
419 if (SCM_I_INUMP (denominator
))
422 y
= SCM_I_INUM (denominator
);
424 return SCM_I_MAKINUM(1);
426 return SCM_I_MAKINUM (x
/ y
);
430 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
431 of that value for the denominator, as a bignum. Apart from
432 that case, abs(bignum) > abs(inum) so inum/bignum is not an
434 if (x
== SCM_MOST_NEGATIVE_FIXNUM
435 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
436 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
437 return SCM_I_MAKINUM(-1);
440 else if (SCM_BIGP (numerator
))
442 if (SCM_I_INUMP (denominator
))
444 long yy
= SCM_I_INUM (denominator
);
445 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
446 return scm_divide (numerator
, denominator
);
450 if (scm_is_eq (numerator
, denominator
))
451 return SCM_I_MAKINUM(1);
452 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
453 SCM_I_BIG_MPZ (denominator
)))
454 return scm_divide(numerator
, denominator
);
458 /* No, it's a proper fraction.
461 SCM divisor
= scm_gcd (numerator
, denominator
);
462 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
464 numerator
= scm_divide (numerator
, divisor
);
465 denominator
= scm_divide (denominator
, divisor
);
468 return scm_double_cell (scm_tc16_fraction
,
469 SCM_UNPACK (numerator
),
470 SCM_UNPACK (denominator
), 0);
476 scm_i_fraction2double (SCM z
)
478 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
479 SCM_FRACTION_DENOMINATOR (z
)));
482 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
484 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
486 #define FUNC_NAME s_scm_exact_p
492 if (SCM_FRACTIONP (x
))
496 SCM_WRONG_TYPE_ARG (1, x
);
501 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
503 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
505 #define FUNC_NAME s_scm_odd_p
509 long val
= SCM_I_INUM (n
);
510 return scm_from_bool ((val
& 1L) != 0);
512 else if (SCM_BIGP (n
))
514 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
515 scm_remember_upto_here_1 (n
);
516 return scm_from_bool (odd_p
);
518 else if (scm_is_true (scm_inf_p (n
)))
520 else if (SCM_REALP (n
))
522 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
528 SCM_WRONG_TYPE_ARG (1, n
);
531 SCM_WRONG_TYPE_ARG (1, n
);
536 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
538 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
540 #define FUNC_NAME s_scm_even_p
544 long val
= SCM_I_INUM (n
);
545 return scm_from_bool ((val
& 1L) == 0);
547 else if (SCM_BIGP (n
))
549 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
550 scm_remember_upto_here_1 (n
);
551 return scm_from_bool (even_p
);
553 else if (scm_is_true (scm_inf_p (n
)))
555 else if (SCM_REALP (n
))
557 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
563 SCM_WRONG_TYPE_ARG (1, n
);
566 SCM_WRONG_TYPE_ARG (1, n
);
570 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
572 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
573 "or @samp{-inf.0}, @code{#f} otherwise.")
574 #define FUNC_NAME s_scm_inf_p
577 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
578 else if (SCM_COMPLEXP (x
))
579 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
580 || xisinf (SCM_COMPLEX_IMAG (x
)));
586 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
588 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
590 #define FUNC_NAME s_scm_nan_p
593 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
594 else if (SCM_COMPLEXP (n
))
595 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
596 || xisnan (SCM_COMPLEX_IMAG (n
)));
602 /* Guile's idea of infinity. */
603 static double guile_Inf
;
605 /* Guile's idea of not a number. */
606 static double guile_NaN
;
609 guile_ieee_init (void)
611 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
613 /* Some version of gcc on some old version of Linux used to crash when
614 trying to make Inf and NaN. */
617 /* C99 INFINITY, when available.
618 FIXME: The standard allows for INFINITY to be something that overflows
619 at compile time. We ought to have a configure test to check for that
620 before trying to use it. (But in practice we believe this is not a
621 problem on any system guile is likely to target.) */
622 guile_Inf
= INFINITY
;
625 extern unsigned int DINFINITY
[2];
626 guile_Inf
= (*((double *) (DINFINITY
)));
633 if (guile_Inf
== tmp
)
641 #if defined (HAVE_ISNAN)
644 /* C99 NAN, when available */
649 extern unsigned int DQNAN
[2];
650 guile_NaN
= (*((double *)(DQNAN
)));
653 guile_NaN
= guile_Inf
/ guile_Inf
;
659 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
662 #define FUNC_NAME s_scm_inf
664 static int initialized
= 0;
670 return scm_from_double (guile_Inf
);
674 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
677 #define FUNC_NAME s_scm_nan
679 static int initialized
= 0;
685 return scm_from_double (guile_NaN
);
690 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
692 "Return the absolute value of @var{x}.")
697 long int xx
= SCM_I_INUM (x
);
700 else if (SCM_POSFIXABLE (-xx
))
701 return SCM_I_MAKINUM (-xx
);
703 return scm_i_long2big (-xx
);
705 else if (SCM_BIGP (x
))
707 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
709 return scm_i_clonebig (x
, 0);
713 else if (SCM_REALP (x
))
715 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
716 double xx
= SCM_REAL_VALUE (x
);
718 return scm_from_double (-xx
);
722 else if (SCM_FRACTIONP (x
))
724 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
726 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
727 SCM_FRACTION_DENOMINATOR (x
));
730 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
735 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
736 /* "Return the quotient of the numbers @var{x} and @var{y}."
739 scm_quotient (SCM x
, SCM y
)
743 long xx
= SCM_I_INUM (x
);
746 long yy
= SCM_I_INUM (y
);
748 scm_num_overflow (s_quotient
);
753 return SCM_I_MAKINUM (z
);
755 return scm_i_long2big (z
);
758 else if (SCM_BIGP (y
))
760 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
761 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
762 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
764 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
765 scm_remember_upto_here_1 (y
);
766 return SCM_I_MAKINUM (-1);
769 return SCM_I_MAKINUM (0);
772 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
774 else if (SCM_BIGP (x
))
778 long yy
= SCM_I_INUM (y
);
780 scm_num_overflow (s_quotient
);
785 SCM result
= scm_i_mkbig ();
788 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
791 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
794 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
795 scm_remember_upto_here_1 (x
);
796 return scm_i_normbig (result
);
799 else if (SCM_BIGP (y
))
801 SCM result
= scm_i_mkbig ();
802 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
805 scm_remember_upto_here_2 (x
, y
);
806 return scm_i_normbig (result
);
809 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
812 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
815 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
816 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
818 * "(remainder 13 4) @result{} 1\n"
819 * "(remainder -13 4) @result{} -1\n"
823 scm_remainder (SCM x
, SCM y
)
829 long yy
= SCM_I_INUM (y
);
831 scm_num_overflow (s_remainder
);
834 long z
= SCM_I_INUM (x
) % yy
;
835 return SCM_I_MAKINUM (z
);
838 else if (SCM_BIGP (y
))
840 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
841 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
842 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
844 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
845 scm_remember_upto_here_1 (y
);
846 return SCM_I_MAKINUM (0);
852 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
854 else if (SCM_BIGP (x
))
858 long yy
= SCM_I_INUM (y
);
860 scm_num_overflow (s_remainder
);
863 SCM result
= scm_i_mkbig ();
866 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
867 scm_remember_upto_here_1 (x
);
868 return scm_i_normbig (result
);
871 else if (SCM_BIGP (y
))
873 SCM result
= scm_i_mkbig ();
874 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
877 scm_remember_upto_here_2 (x
, y
);
878 return scm_i_normbig (result
);
881 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
884 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
888 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
889 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
891 * "(modulo 13 4) @result{} 1\n"
892 * "(modulo -13 4) @result{} 3\n"
896 scm_modulo (SCM x
, SCM y
)
900 long xx
= SCM_I_INUM (x
);
903 long yy
= SCM_I_INUM (y
);
905 scm_num_overflow (s_modulo
);
908 /* C99 specifies that "%" is the remainder corresponding to a
909 quotient rounded towards zero, and that's also traditional
910 for machine division, so z here should be well defined. */
928 return SCM_I_MAKINUM (result
);
931 else if (SCM_BIGP (y
))
933 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
940 SCM pos_y
= scm_i_clonebig (y
, 0);
941 /* do this after the last scm_op */
942 mpz_init_set_si (z_x
, xx
);
943 result
= pos_y
; /* re-use this bignum */
944 mpz_mod (SCM_I_BIG_MPZ (result
),
946 SCM_I_BIG_MPZ (pos_y
));
947 scm_remember_upto_here_1 (pos_y
);
951 result
= scm_i_mkbig ();
952 /* do this after the last scm_op */
953 mpz_init_set_si (z_x
, xx
);
954 mpz_mod (SCM_I_BIG_MPZ (result
),
957 scm_remember_upto_here_1 (y
);
960 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
961 mpz_add (SCM_I_BIG_MPZ (result
),
963 SCM_I_BIG_MPZ (result
));
964 scm_remember_upto_here_1 (y
);
965 /* and do this before the next one */
967 return scm_i_normbig (result
);
971 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
973 else if (SCM_BIGP (x
))
977 long yy
= SCM_I_INUM (y
);
979 scm_num_overflow (s_modulo
);
982 SCM result
= scm_i_mkbig ();
983 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
985 (yy
< 0) ? - yy
: yy
);
986 scm_remember_upto_here_1 (x
);
987 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
988 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
989 SCM_I_BIG_MPZ (result
),
991 return scm_i_normbig (result
);
994 else if (SCM_BIGP (y
))
997 SCM result
= scm_i_mkbig ();
998 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
999 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1000 mpz_mod (SCM_I_BIG_MPZ (result
),
1002 SCM_I_BIG_MPZ (pos_y
));
1004 scm_remember_upto_here_1 (x
);
1005 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1006 mpz_add (SCM_I_BIG_MPZ (result
),
1008 SCM_I_BIG_MPZ (result
));
1009 scm_remember_upto_here_2 (y
, pos_y
);
1010 return scm_i_normbig (result
);
1014 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1017 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1020 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
1021 /* "Return the greatest common divisor of all arguments.\n"
1022 * "If called without arguments, 0 is returned."
1025 scm_gcd (SCM x
, SCM y
)
1028 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1030 if (SCM_I_INUMP (x
))
1032 if (SCM_I_INUMP (y
))
1034 long xx
= SCM_I_INUM (x
);
1035 long yy
= SCM_I_INUM (y
);
1036 long u
= xx
< 0 ? -xx
: xx
;
1037 long v
= yy
< 0 ? -yy
: yy
;
1047 /* Determine a common factor 2^k */
1048 while (!(1 & (u
| v
)))
1054 /* Now, any factor 2^n can be eliminated */
1074 return (SCM_POSFIXABLE (result
)
1075 ? SCM_I_MAKINUM (result
)
1076 : scm_i_long2big (result
));
1078 else if (SCM_BIGP (y
))
1084 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1086 else if (SCM_BIGP (x
))
1088 if (SCM_I_INUMP (y
))
1090 unsigned long result
;
1093 yy
= SCM_I_INUM (y
);
1098 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1099 scm_remember_upto_here_1 (x
);
1100 return (SCM_POSFIXABLE (result
)
1101 ? SCM_I_MAKINUM (result
)
1102 : scm_from_ulong (result
));
1104 else if (SCM_BIGP (y
))
1106 SCM result
= scm_i_mkbig ();
1107 mpz_gcd (SCM_I_BIG_MPZ (result
),
1110 scm_remember_upto_here_2 (x
, y
);
1111 return scm_i_normbig (result
);
1114 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1117 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1120 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1121 /* "Return the least common multiple of the arguments.\n"
1122 * "If called without arguments, 1 is returned."
1125 scm_lcm (SCM n1
, SCM n2
)
1127 if (SCM_UNBNDP (n2
))
1129 if (SCM_UNBNDP (n1
))
1130 return SCM_I_MAKINUM (1L);
1131 n2
= SCM_I_MAKINUM (1L);
1134 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1135 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1136 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1137 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1139 if (SCM_I_INUMP (n1
))
1141 if (SCM_I_INUMP (n2
))
1143 SCM d
= scm_gcd (n1
, n2
);
1144 if (scm_is_eq (d
, SCM_INUM0
))
1147 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1151 /* inum n1, big n2 */
1154 SCM result
= scm_i_mkbig ();
1155 long nn1
= SCM_I_INUM (n1
);
1156 if (nn1
== 0) return SCM_INUM0
;
1157 if (nn1
< 0) nn1
= - nn1
;
1158 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1159 scm_remember_upto_here_1 (n2
);
1167 if (SCM_I_INUMP (n2
))
1174 SCM result
= scm_i_mkbig ();
1175 mpz_lcm(SCM_I_BIG_MPZ (result
),
1177 SCM_I_BIG_MPZ (n2
));
1178 scm_remember_upto_here_2(n1
, n2
);
1179 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1185 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1190 + + + x (map digit:logand X Y)
1191 + - + x (map digit:logand X (lognot (+ -1 Y)))
1192 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1193 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1198 + + + (map digit:logior X Y)
1199 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1200 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1201 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1206 + + + (map digit:logxor X Y)
1207 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1208 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1209 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1214 + + (any digit:logand X Y)
1215 + - (any digit:logand X (lognot (+ -1 Y)))
1216 - + (any digit:logand (lognot (+ -1 X)) Y)
1221 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1223 "Return the bitwise AND of the integer arguments.\n\n"
1225 "(logand) @result{} -1\n"
1226 "(logand 7) @result{} 7\n"
1227 "(logand #b111 #b011 #b001) @result{} 1\n"
1229 #define FUNC_NAME s_scm_logand
1233 if (SCM_UNBNDP (n2
))
1235 if (SCM_UNBNDP (n1
))
1236 return SCM_I_MAKINUM (-1);
1237 else if (!SCM_NUMBERP (n1
))
1238 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1239 else if (SCM_NUMBERP (n1
))
1242 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1245 if (SCM_I_INUMP (n1
))
1247 nn1
= SCM_I_INUM (n1
);
1248 if (SCM_I_INUMP (n2
))
1250 long nn2
= SCM_I_INUM (n2
);
1251 return SCM_I_MAKINUM (nn1
& nn2
);
1253 else if SCM_BIGP (n2
)
1259 SCM result_z
= scm_i_mkbig ();
1261 mpz_init_set_si (nn1_z
, nn1
);
1262 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1263 scm_remember_upto_here_1 (n2
);
1265 return scm_i_normbig (result_z
);
1269 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1271 else if (SCM_BIGP (n1
))
1273 if (SCM_I_INUMP (n2
))
1276 nn1
= SCM_I_INUM (n1
);
1279 else if (SCM_BIGP (n2
))
1281 SCM result_z
= scm_i_mkbig ();
1282 mpz_and (SCM_I_BIG_MPZ (result_z
),
1284 SCM_I_BIG_MPZ (n2
));
1285 scm_remember_upto_here_2 (n1
, n2
);
1286 return scm_i_normbig (result_z
);
1289 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1292 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1297 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1299 "Return the bitwise OR of the integer arguments.\n\n"
1301 "(logior) @result{} 0\n"
1302 "(logior 7) @result{} 7\n"
1303 "(logior #b000 #b001 #b011) @result{} 3\n"
1305 #define FUNC_NAME s_scm_logior
1309 if (SCM_UNBNDP (n2
))
1311 if (SCM_UNBNDP (n1
))
1313 else if (SCM_NUMBERP (n1
))
1316 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1319 if (SCM_I_INUMP (n1
))
1321 nn1
= SCM_I_INUM (n1
);
1322 if (SCM_I_INUMP (n2
))
1324 long nn2
= SCM_I_INUM (n2
);
1325 return SCM_I_MAKINUM (nn1
| nn2
);
1327 else if (SCM_BIGP (n2
))
1333 SCM result_z
= scm_i_mkbig ();
1335 mpz_init_set_si (nn1_z
, nn1
);
1336 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1337 scm_remember_upto_here_1 (n2
);
1339 return scm_i_normbig (result_z
);
1343 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1345 else if (SCM_BIGP (n1
))
1347 if (SCM_I_INUMP (n2
))
1350 nn1
= SCM_I_INUM (n1
);
1353 else if (SCM_BIGP (n2
))
1355 SCM result_z
= scm_i_mkbig ();
1356 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1358 SCM_I_BIG_MPZ (n2
));
1359 scm_remember_upto_here_2 (n1
, n2
);
1360 return scm_i_normbig (result_z
);
1363 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1366 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1371 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1373 "Return the bitwise XOR of the integer arguments. A bit is\n"
1374 "set in the result if it is set in an odd number of arguments.\n"
1376 "(logxor) @result{} 0\n"
1377 "(logxor 7) @result{} 7\n"
1378 "(logxor #b000 #b001 #b011) @result{} 2\n"
1379 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1381 #define FUNC_NAME s_scm_logxor
1385 if (SCM_UNBNDP (n2
))
1387 if (SCM_UNBNDP (n1
))
1389 else if (SCM_NUMBERP (n1
))
1392 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1395 if (SCM_I_INUMP (n1
))
1397 nn1
= SCM_I_INUM (n1
);
1398 if (SCM_I_INUMP (n2
))
1400 long nn2
= SCM_I_INUM (n2
);
1401 return SCM_I_MAKINUM (nn1
^ nn2
);
1403 else if (SCM_BIGP (n2
))
1407 SCM result_z
= scm_i_mkbig ();
1409 mpz_init_set_si (nn1_z
, nn1
);
1410 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1411 scm_remember_upto_here_1 (n2
);
1413 return scm_i_normbig (result_z
);
1417 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1419 else if (SCM_BIGP (n1
))
1421 if (SCM_I_INUMP (n2
))
1424 nn1
= SCM_I_INUM (n1
);
1427 else if (SCM_BIGP (n2
))
1429 SCM result_z
= scm_i_mkbig ();
1430 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1432 SCM_I_BIG_MPZ (n2
));
1433 scm_remember_upto_here_2 (n1
, n2
);
1434 return scm_i_normbig (result_z
);
1437 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1440 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1445 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1447 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1448 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1449 "without actually calculating the @code{logand}, just testing\n"
1453 "(logtest #b0100 #b1011) @result{} #f\n"
1454 "(logtest #b0100 #b0111) @result{} #t\n"
1456 #define FUNC_NAME s_scm_logtest
1460 if (SCM_I_INUMP (j
))
1462 nj
= SCM_I_INUM (j
);
1463 if (SCM_I_INUMP (k
))
1465 long nk
= SCM_I_INUM (k
);
1466 return scm_from_bool (nj
& nk
);
1468 else if (SCM_BIGP (k
))
1476 mpz_init_set_si (nj_z
, nj
);
1477 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1478 scm_remember_upto_here_1 (k
);
1479 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1485 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1487 else if (SCM_BIGP (j
))
1489 if (SCM_I_INUMP (k
))
1492 nj
= SCM_I_INUM (j
);
1495 else if (SCM_BIGP (k
))
1499 mpz_init (result_z
);
1503 scm_remember_upto_here_2 (j
, k
);
1504 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1505 mpz_clear (result_z
);
1509 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1512 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1517 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1519 "Test whether bit number @var{index} in @var{j} is set.\n"
1520 "@var{index} starts from 0 for the least significant bit.\n"
1523 "(logbit? 0 #b1101) @result{} #t\n"
1524 "(logbit? 1 #b1101) @result{} #f\n"
1525 "(logbit? 2 #b1101) @result{} #t\n"
1526 "(logbit? 3 #b1101) @result{} #t\n"
1527 "(logbit? 4 #b1101) @result{} #f\n"
1529 #define FUNC_NAME s_scm_logbit_p
1531 unsigned long int iindex
;
1532 iindex
= scm_to_ulong (index
);
1534 if (SCM_I_INUMP (j
))
1536 /* bits above what's in an inum follow the sign bit */
1537 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1538 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1540 else if (SCM_BIGP (j
))
1542 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1543 scm_remember_upto_here_1 (j
);
1544 return scm_from_bool (val
);
1547 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1552 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1554 "Return the integer which is the ones-complement of the integer\n"
1558 "(number->string (lognot #b10000000) 2)\n"
1559 " @result{} \"-10000001\"\n"
1560 "(number->string (lognot #b0) 2)\n"
1561 " @result{} \"-1\"\n"
1563 #define FUNC_NAME s_scm_lognot
1565 if (SCM_I_INUMP (n
)) {
1566 /* No overflow here, just need to toggle all the bits making up the inum.
1567 Enhancement: No need to strip the tag and add it back, could just xor
1568 a block of 1 bits, if that worked with the various debug versions of
1570 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1572 } else if (SCM_BIGP (n
)) {
1573 SCM result
= scm_i_mkbig ();
1574 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1575 scm_remember_upto_here_1 (n
);
1579 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1584 /* returns 0 if IN is not an integer. OUT must already be
1587 coerce_to_big (SCM in
, mpz_t out
)
1590 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1591 else if (SCM_I_INUMP (in
))
1592 mpz_set_si (out
, SCM_I_INUM (in
));
1599 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1600 (SCM n
, SCM k
, SCM m
),
1601 "Return @var{n} raised to the integer exponent\n"
1602 "@var{k}, modulo @var{m}.\n"
1605 "(modulo-expt 2 3 5)\n"
1608 #define FUNC_NAME s_scm_modulo_expt
1614 /* There are two classes of error we might encounter --
1615 1) Math errors, which we'll report by calling scm_num_overflow,
1617 2) wrong-type errors, which of course we'll report by calling
1619 We don't report those errors immediately, however; instead we do
1620 some cleanup first. These variables tell us which error (if
1621 any) we should report after cleaning up.
1623 int report_overflow
= 0;
1625 int position_of_wrong_type
= 0;
1626 SCM value_of_wrong_type
= SCM_INUM0
;
1628 SCM result
= SCM_UNDEFINED
;
1634 if (scm_is_eq (m
, SCM_INUM0
))
1636 report_overflow
= 1;
1640 if (!coerce_to_big (n
, n_tmp
))
1642 value_of_wrong_type
= n
;
1643 position_of_wrong_type
= 1;
1647 if (!coerce_to_big (k
, k_tmp
))
1649 value_of_wrong_type
= k
;
1650 position_of_wrong_type
= 2;
1654 if (!coerce_to_big (m
, m_tmp
))
1656 value_of_wrong_type
= m
;
1657 position_of_wrong_type
= 3;
1661 /* if the exponent K is negative, and we simply call mpz_powm, we
1662 will get a divide-by-zero exception when an inverse 1/n mod m
1663 doesn't exist (or is not unique). Since exceptions are hard to
1664 handle, we'll attempt the inversion "by hand" -- that way, we get
1665 a simple failure code, which is easy to handle. */
1667 if (-1 == mpz_sgn (k_tmp
))
1669 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1671 report_overflow
= 1;
1674 mpz_neg (k_tmp
, k_tmp
);
1677 result
= scm_i_mkbig ();
1678 mpz_powm (SCM_I_BIG_MPZ (result
),
1683 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1684 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1691 if (report_overflow
)
1692 scm_num_overflow (FUNC_NAME
);
1694 if (position_of_wrong_type
)
1695 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1696 value_of_wrong_type
);
1698 return scm_i_normbig (result
);
1702 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1704 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1705 "exact integer, @var{n} can be any number.\n"
1707 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1708 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1709 "includes @math{0^0} is 1.\n"
1712 "(integer-expt 2 5) @result{} 32\n"
1713 "(integer-expt -3 3) @result{} -27\n"
1714 "(integer-expt 5 -3) @result{} 1/125\n"
1715 "(integer-expt 0 0) @result{} 1\n"
1717 #define FUNC_NAME s_scm_integer_expt
1720 SCM z_i2
= SCM_BOOL_F
;
1722 SCM acc
= SCM_I_MAKINUM (1L);
1724 /* 0^0 == 1 according to R5RS */
1725 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1726 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1727 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1728 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1730 if (SCM_I_INUMP (k
))
1731 i2
= SCM_I_INUM (k
);
1732 else if (SCM_BIGP (k
))
1734 z_i2
= scm_i_clonebig (k
, 1);
1735 scm_remember_upto_here_1 (k
);
1739 SCM_WRONG_TYPE_ARG (2, k
);
1743 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1745 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1746 n
= scm_divide (n
, SCM_UNDEFINED
);
1750 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1754 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1756 return scm_product (acc
, n
);
1758 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1759 acc
= scm_product (acc
, n
);
1760 n
= scm_product (n
, n
);
1761 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1769 n
= scm_divide (n
, SCM_UNDEFINED
);
1776 return scm_product (acc
, n
);
1778 acc
= scm_product (acc
, n
);
1779 n
= scm_product (n
, n
);
1786 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1788 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1789 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1791 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1792 "@var{cnt} is negative it's a division, rounded towards negative\n"
1793 "infinity. (Note that this is not the same rounding as\n"
1794 "@code{quotient} does.)\n"
1796 "With @var{n} viewed as an infinite precision twos complement,\n"
1797 "@code{ash} means a left shift introducing zero bits, or a right\n"
1798 "shift dropping bits.\n"
1801 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1802 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1804 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1805 "(ash -23 -2) @result{} -6\n"
1807 #define FUNC_NAME s_scm_ash
1810 bits_to_shift
= scm_to_long (cnt
);
1812 if (SCM_I_INUMP (n
))
1814 long nn
= SCM_I_INUM (n
);
1816 if (bits_to_shift
> 0)
1818 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1819 overflow a non-zero fixnum. For smaller shifts we check the
1820 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1821 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1822 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1828 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1830 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1833 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1837 SCM result
= scm_i_long2big (nn
);
1838 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1845 bits_to_shift
= -bits_to_shift
;
1846 if (bits_to_shift
>= SCM_LONG_BIT
)
1847 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1849 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1853 else if (SCM_BIGP (n
))
1857 if (bits_to_shift
== 0)
1860 result
= scm_i_mkbig ();
1861 if (bits_to_shift
>= 0)
1863 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1869 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1870 we have to allocate a bignum even if the result is going to be a
1872 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1874 return scm_i_normbig (result
);
1880 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1886 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1887 (SCM n
, SCM start
, SCM end
),
1888 "Return the integer composed of the @var{start} (inclusive)\n"
1889 "through @var{end} (exclusive) bits of @var{n}. The\n"
1890 "@var{start}th bit becomes the 0-th bit in the result.\n"
1893 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1894 " @result{} \"1010\"\n"
1895 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1896 " @result{} \"10110\"\n"
1898 #define FUNC_NAME s_scm_bit_extract
1900 unsigned long int istart
, iend
, bits
;
1901 istart
= scm_to_ulong (start
);
1902 iend
= scm_to_ulong (end
);
1903 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1905 /* how many bits to keep */
1906 bits
= iend
- istart
;
1908 if (SCM_I_INUMP (n
))
1910 long int in
= SCM_I_INUM (n
);
1912 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1913 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1914 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1916 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1918 /* Since we emulate two's complement encoded numbers, this
1919 * special case requires us to produce a result that has
1920 * more bits than can be stored in a fixnum.
1922 SCM result
= scm_i_long2big (in
);
1923 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1928 /* mask down to requisite bits */
1929 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1930 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1932 else if (SCM_BIGP (n
))
1937 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1941 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1942 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1943 such bits into a ulong. */
1944 result
= scm_i_mkbig ();
1945 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1946 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1947 result
= scm_i_normbig (result
);
1949 scm_remember_upto_here_1 (n
);
1953 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1958 static const char scm_logtab
[] = {
1959 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1962 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1964 "Return the number of bits in integer @var{n}. If integer is\n"
1965 "positive, the 1-bits in its binary representation are counted.\n"
1966 "If negative, the 0-bits in its two's-complement binary\n"
1967 "representation are counted. If 0, 0 is returned.\n"
1970 "(logcount #b10101010)\n"
1977 #define FUNC_NAME s_scm_logcount
1979 if (SCM_I_INUMP (n
))
1981 unsigned long int c
= 0;
1982 long int nn
= SCM_I_INUM (n
);
1987 c
+= scm_logtab
[15 & nn
];
1990 return SCM_I_MAKINUM (c
);
1992 else if (SCM_BIGP (n
))
1994 unsigned long count
;
1995 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1996 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1998 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1999 scm_remember_upto_here_1 (n
);
2000 return SCM_I_MAKINUM (count
);
2003 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2008 static const char scm_ilentab
[] = {
2009 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2013 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2015 "Return the number of bits necessary to represent @var{n}.\n"
2018 "(integer-length #b10101010)\n"
2020 "(integer-length 0)\n"
2022 "(integer-length #b1111)\n"
2025 #define FUNC_NAME s_scm_integer_length
2027 if (SCM_I_INUMP (n
))
2029 unsigned long int c
= 0;
2031 long int nn
= SCM_I_INUM (n
);
2037 l
= scm_ilentab
[15 & nn
];
2040 return SCM_I_MAKINUM (c
- 4 + l
);
2042 else if (SCM_BIGP (n
))
2044 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2045 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2046 1 too big, so check for that and adjust. */
2047 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2048 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2049 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2050 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2052 scm_remember_upto_here_1 (n
);
2053 return SCM_I_MAKINUM (size
);
2056 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2060 /*** NUMBERS -> STRINGS ***/
2061 #define SCM_MAX_DBL_PREC 60
2062 #define SCM_MAX_DBL_RADIX 36
2064 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2065 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2066 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2069 void init_dblprec(int *prec
, int radix
) {
2070 /* determine floating point precision by adding successively
2071 smaller increments to 1.0 until it is considered == 1.0 */
2072 double f
= ((double)1.0)/radix
;
2073 double fsum
= 1.0 + f
;
2078 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2090 void init_fx_radix(double *fx_list
, int radix
)
2092 /* initialize a per-radix list of tolerances. When added
2093 to a number < 1.0, we can determine if we should raund
2094 up and quit converting a number to a string. */
2098 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2099 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2102 /* use this array as a way to generate a single digit */
2103 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2106 idbl2str (double f
, char *a
, int radix
)
2108 int efmt
, dpt
, d
, i
, wp
;
2110 #ifdef DBL_MIN_10_EXP
2113 #endif /* DBL_MIN_10_EXP */
2118 radix
> SCM_MAX_DBL_RADIX
)
2120 /* revert to existing behavior */
2124 wp
= scm_dblprec
[radix
-2];
2125 fx
= fx_per_radix
[radix
-2];
2129 #ifdef HAVE_COPYSIGN
2130 double sgn
= copysign (1.0, f
);
2135 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2141 strcpy (a
, "-inf.0");
2143 strcpy (a
, "+inf.0");
2146 else if (xisnan (f
))
2148 strcpy (a
, "+nan.0");
2158 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2159 make-uniform-vector, from causing infinite loops. */
2160 /* just do the checking...if it passes, we do the conversion for our
2161 radix again below */
2168 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2176 while (f_cpy
> 10.0)
2179 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2200 if (f
+ fx
[wp
] >= radix
)
2207 /* adding 9999 makes this equivalent to abs(x) % 3 */
2208 dpt
= (exp
+ 9999) % 3;
2212 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2234 a
[ch
++] = number_chars
[d
];
2237 if (f
+ fx
[wp
] >= 1.0)
2239 a
[ch
- 1] = number_chars
[d
+1];
2251 if ((dpt
> 4) && (exp
> 6))
2253 d
= (a
[0] == '-' ? 2 : 1);
2254 for (i
= ch
++; i
> d
; i
--)
2267 if (a
[ch
- 1] == '.')
2268 a
[ch
++] = '0'; /* trailing zero */
2277 for (i
= radix
; i
<= exp
; i
*= radix
);
2278 for (i
/= radix
; i
; i
/= radix
)
2280 a
[ch
++] = number_chars
[exp
/ i
];
2289 icmplx2str (double real
, double imag
, char *str
, int radix
)
2293 i
= idbl2str (real
, str
, radix
);
2296 /* Don't output a '+' for negative numbers or for Inf and
2297 NaN. They will provide their own sign. */
2298 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2300 i
+= idbl2str (imag
, &str
[i
], radix
);
2307 iflo2str (SCM flt
, char *str
, int radix
)
2310 if (SCM_REALP (flt
))
2311 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2313 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2318 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2319 characters in the result.
2321 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2323 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2328 return scm_iuint2str (-num
, rad
, p
) + 1;
2331 return scm_iuint2str (num
, rad
, p
);
2334 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2335 characters in the result.
2337 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2339 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2343 scm_t_uintmax n
= num
;
2345 for (n
/= rad
; n
> 0; n
/= rad
)
2355 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2360 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2362 "Return a string holding the external representation of the\n"
2363 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2364 "inexact, a radix of 10 will be used.")
2365 #define FUNC_NAME s_scm_number_to_string
2369 if (SCM_UNBNDP (radix
))
2372 base
= scm_to_signed_integer (radix
, 2, 36);
2374 if (SCM_I_INUMP (n
))
2376 char num_buf
[SCM_INTBUFLEN
];
2377 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2378 return scm_from_locale_stringn (num_buf
, length
);
2380 else if (SCM_BIGP (n
))
2382 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2383 scm_remember_upto_here_1 (n
);
2384 return scm_take_locale_string (str
);
2386 else if (SCM_FRACTIONP (n
))
2388 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2389 scm_from_locale_string ("/"),
2390 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2392 else if (SCM_INEXACTP (n
))
2394 char num_buf
[FLOBUFLEN
];
2395 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2398 SCM_WRONG_TYPE_ARG (1, n
);
2403 /* These print routines used to be stubbed here so that scm_repl.c
2404 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2407 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2409 char num_buf
[FLOBUFLEN
];
2410 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2415 scm_i_print_double (double val
, SCM port
)
2417 char num_buf
[FLOBUFLEN
];
2418 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2422 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2425 char num_buf
[FLOBUFLEN
];
2426 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2431 scm_i_print_complex (double real
, double imag
, SCM port
)
2433 char num_buf
[FLOBUFLEN
];
2434 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2438 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2441 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2442 scm_lfwrite (scm_i_string_chars (str
), scm_i_string_length (str
), port
);
2443 scm_remember_upto_here_1 (str
);
2448 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2450 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2451 scm_remember_upto_here_1 (exp
);
2452 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2456 /*** END nums->strs ***/
2459 /*** STRINGS -> NUMBERS ***/
2461 /* The following functions implement the conversion from strings to numbers.
2462 * The implementation somehow follows the grammar for numbers as it is given
2463 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2464 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2465 * points should be noted about the implementation:
2466 * * Each function keeps a local index variable 'idx' that points at the
2467 * current position within the parsed string. The global index is only
2468 * updated if the function could parse the corresponding syntactic unit
2470 * * Similarly, the functions keep track of indicators of inexactness ('#',
2471 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2472 * global exactness information is only updated after each part has been
2473 * successfully parsed.
2474 * * Sequences of digits are parsed into temporary variables holding fixnums.
2475 * Only if these fixnums would overflow, the result variables are updated
2476 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2477 * the temporary variables holding the fixnums are cleared, and the process
2478 * starts over again. If for example fixnums were able to store five decimal
2479 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2480 * and the result was computed as 12345 * 100000 + 67890. In other words,
2481 * only every five digits two bignum operations were performed.
2484 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2486 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2488 /* In non ASCII-style encodings the following macro might not work. */
2489 #define XDIGIT2UINT(d) \
2490 (isdigit ((int) (unsigned char) d) \
2492 : tolower ((int) (unsigned char) d) - 'a' + 10)
2495 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2496 unsigned int radix
, enum t_exactness
*p_exactness
)
2498 unsigned int idx
= *p_idx
;
2499 unsigned int hash_seen
= 0;
2500 scm_t_bits shift
= 1;
2502 unsigned int digit_value
;
2510 if (!isxdigit ((int) (unsigned char) c
))
2512 digit_value
= XDIGIT2UINT (c
);
2513 if (digit_value
>= radix
)
2517 result
= SCM_I_MAKINUM (digit_value
);
2521 if (isxdigit ((int) (unsigned char) c
))
2525 digit_value
= XDIGIT2UINT (c
);
2526 if (digit_value
>= radix
)
2538 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2540 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2542 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2549 shift
= shift
* radix
;
2550 add
= add
* radix
+ digit_value
;
2555 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2557 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2561 *p_exactness
= INEXACT
;
2567 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2568 * covers the parts of the rules that start at a potential point. The value
2569 * of the digits up to the point have been parsed by the caller and are given
2570 * in variable result. The content of *p_exactness indicates, whether a hash
2571 * has already been seen in the digits before the point.
2574 /* In non ASCII-style encodings the following macro might not work. */
2575 #define DIGIT2UINT(d) ((d) - '0')
2578 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2579 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2581 unsigned int idx
= *p_idx
;
2582 enum t_exactness x
= *p_exactness
;
2587 if (mem
[idx
] == '.')
2589 scm_t_bits shift
= 1;
2591 unsigned int digit_value
;
2592 SCM big_shift
= SCM_I_MAKINUM (1);
2598 if (isdigit ((int) (unsigned char) c
))
2603 digit_value
= DIGIT2UINT (c
);
2614 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2616 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2617 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2619 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2627 add
= add
* 10 + digit_value
;
2633 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2634 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2635 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2638 result
= scm_divide (result
, big_shift
);
2640 /* We've seen a decimal point, thus the value is implicitly inexact. */
2652 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2679 if (!isdigit ((int) (unsigned char) c
))
2683 exponent
= DIGIT2UINT (c
);
2687 if (isdigit ((int) (unsigned char) c
))
2690 if (exponent
<= SCM_MAXEXP
)
2691 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2697 if (exponent
> SCM_MAXEXP
)
2699 size_t exp_len
= idx
- start
;
2700 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2701 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2702 scm_out_of_range ("string->number", exp_num
);
2705 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2707 result
= scm_product (result
, e
);
2709 result
= scm_divide2real (result
, e
);
2711 /* We've seen an exponent, thus the value is implicitly inexact. */
2729 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2732 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2733 unsigned int radix
, enum t_exactness
*p_exactness
)
2735 unsigned int idx
= *p_idx
;
2741 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2747 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2749 enum t_exactness x
= EXACT
;
2751 /* Cobble up the fractional part. We might want to set the
2752 NaN's mantissa from it. */
2754 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2759 if (mem
[idx
] == '.')
2763 else if (idx
+ 1 == len
)
2765 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2768 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2769 p_idx
, p_exactness
);
2773 enum t_exactness x
= EXACT
;
2776 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2777 if (scm_is_false (uinteger
))
2782 else if (mem
[idx
] == '/')
2788 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2789 if (scm_is_false (divisor
))
2792 /* both are int/big here, I assume */
2793 result
= scm_i_make_ratio (uinteger
, divisor
);
2795 else if (radix
== 10)
2797 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2798 if (scm_is_false (result
))
2809 /* When returning an inexact zero, make sure it is represented as a
2810 floating point value so that we can change its sign.
2812 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2813 result
= scm_from_double (0.0);
2819 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2822 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2823 unsigned int radix
, enum t_exactness
*p_exactness
)
2847 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2848 if (scm_is_false (ureal
))
2850 /* input must be either +i or -i */
2855 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2861 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2868 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2869 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2878 /* either +<ureal>i or -<ureal>i */
2885 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2888 /* polar input: <real>@<real>. */
2913 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2914 if (scm_is_false (angle
))
2919 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2920 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2922 result
= scm_make_polar (ureal
, angle
);
2927 /* expecting input matching <real>[+-]<ureal>?i */
2934 int sign
= (c
== '+') ? 1 : -1;
2935 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2937 if (scm_is_false (imag
))
2938 imag
= SCM_I_MAKINUM (sign
);
2939 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2940 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2944 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2951 return scm_make_rectangular (ureal
, imag
);
2960 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2962 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2965 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
2966 unsigned int default_radix
)
2968 unsigned int idx
= 0;
2969 unsigned int radix
= NO_RADIX
;
2970 enum t_exactness forced_x
= NO_EXACTNESS
;
2971 enum t_exactness implicit_x
= EXACT
;
2974 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2975 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2977 switch (mem
[idx
+ 1])
2980 if (radix
!= NO_RADIX
)
2985 if (radix
!= NO_RADIX
)
2990 if (forced_x
!= NO_EXACTNESS
)
2995 if (forced_x
!= NO_EXACTNESS
)
3000 if (radix
!= NO_RADIX
)
3005 if (radix
!= NO_RADIX
)
3015 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3016 if (radix
== NO_RADIX
)
3017 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
3019 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
3021 if (scm_is_false (result
))
3027 if (SCM_INEXACTP (result
))
3028 return scm_inexact_to_exact (result
);
3032 if (SCM_INEXACTP (result
))
3035 return scm_exact_to_inexact (result
);
3038 if (implicit_x
== INEXACT
)
3040 if (SCM_INEXACTP (result
))
3043 return scm_exact_to_inexact (result
);
3051 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3052 (SCM string
, SCM radix
),
3053 "Return a number of the maximally precise representation\n"
3054 "expressed by the given @var{string}. @var{radix} must be an\n"
3055 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3056 "is a default radix that may be overridden by an explicit radix\n"
3057 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3058 "supplied, then the default radix is 10. If string is not a\n"
3059 "syntactically valid notation for a number, then\n"
3060 "@code{string->number} returns @code{#f}.")
3061 #define FUNC_NAME s_scm_string_to_number
3065 SCM_VALIDATE_STRING (1, string
);
3067 if (SCM_UNBNDP (radix
))
3070 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3072 answer
= scm_c_locale_stringn_to_number (scm_i_string_chars (string
),
3073 scm_i_string_length (string
),
3075 scm_remember_upto_here_1 (string
);
3081 /*** END strs->nums ***/
3085 scm_bigequal (SCM x
, SCM y
)
3087 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3088 scm_remember_upto_here_2 (x
, y
);
3089 return scm_from_bool (0 == result
);
3093 scm_real_equalp (SCM x
, SCM y
)
3095 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3099 scm_complex_equalp (SCM x
, SCM y
)
3101 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3102 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3106 scm_i_fraction_equalp (SCM x
, SCM y
)
3108 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3109 SCM_FRACTION_NUMERATOR (y
)))
3110 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3111 SCM_FRACTION_DENOMINATOR (y
))))
3118 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3120 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3122 #define FUNC_NAME s_scm_number_p
3124 return scm_from_bool (SCM_NUMBERP (x
));
3128 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3130 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3131 "otherwise. Note that the sets of real, rational and integer\n"
3132 "values form subsets of the set of complex numbers, i. e. the\n"
3133 "predicate will also be fulfilled if @var{x} is a real,\n"
3134 "rational or integer number.")
3135 #define FUNC_NAME s_scm_complex_p
3137 /* all numbers are complex. */
3138 return scm_number_p (x
);
3142 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3144 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3145 "otherwise. Note that the set of integer values forms a subset of\n"
3146 "the set of real numbers, i. e. the predicate will also be\n"
3147 "fulfilled if @var{x} is an integer number.")
3148 #define FUNC_NAME s_scm_real_p
3150 /* we can't represent irrational numbers. */
3151 return scm_rational_p (x
);
3155 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3157 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3158 "otherwise. Note that the set of integer values forms a subset of\n"
3159 "the set of rational numbers, i. e. the predicate will also be\n"
3160 "fulfilled if @var{x} is an integer number.")
3161 #define FUNC_NAME s_scm_rational_p
3163 if (SCM_I_INUMP (x
))
3165 else if (SCM_IMP (x
))
3167 else if (SCM_BIGP (x
))
3169 else if (SCM_FRACTIONP (x
))
3171 else if (SCM_REALP (x
))
3172 /* due to their limited precision, all floating point numbers are
3173 rational as well. */
3180 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3182 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3184 #define FUNC_NAME s_scm_integer_p
3187 if (SCM_I_INUMP (x
))
3193 if (!SCM_INEXACTP (x
))
3195 if (SCM_COMPLEXP (x
))
3197 r
= SCM_REAL_VALUE (x
);
3198 /* +/-inf passes r==floor(r), making those #t */
3206 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3208 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3210 #define FUNC_NAME s_scm_inexact_p
3212 if (SCM_INEXACTP (x
))
3214 if (SCM_NUMBERP (x
))
3216 SCM_WRONG_TYPE_ARG (1, x
);
3221 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3222 /* "Return @code{#t} if all parameters are numerically equal." */
3224 scm_num_eq_p (SCM x
, SCM y
)
3227 if (SCM_I_INUMP (x
))
3229 long xx
= SCM_I_INUM (x
);
3230 if (SCM_I_INUMP (y
))
3232 long yy
= SCM_I_INUM (y
);
3233 return scm_from_bool (xx
== yy
);
3235 else if (SCM_BIGP (y
))
3237 else if (SCM_REALP (y
))
3239 /* On a 32-bit system an inum fits a double, we can cast the inum
3240 to a double and compare.
3242 But on a 64-bit system an inum is bigger than a double and
3243 casting it to a double (call that dxx) will round. dxx is at
3244 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3245 an integer and fits a long. So we cast yy to a long and
3246 compare with plain xx.
3248 An alternative (for any size system actually) would be to check
3249 yy is an integer (with floor) and is in range of an inum
3250 (compare against appropriate powers of 2) then test
3251 xx==(long)yy. It's just a matter of which casts/comparisons
3252 might be fastest or easiest for the cpu. */
3254 double yy
= SCM_REAL_VALUE (y
);
3255 return scm_from_bool ((double) xx
== yy
3256 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3257 || xx
== (long) yy
));
3259 else if (SCM_COMPLEXP (y
))
3260 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3261 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3262 else if (SCM_FRACTIONP (y
))
3265 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3267 else if (SCM_BIGP (x
))
3269 if (SCM_I_INUMP (y
))
3271 else if (SCM_BIGP (y
))
3273 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3274 scm_remember_upto_here_2 (x
, y
);
3275 return scm_from_bool (0 == cmp
);
3277 else if (SCM_REALP (y
))
3280 if (xisnan (SCM_REAL_VALUE (y
)))
3282 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3283 scm_remember_upto_here_1 (x
);
3284 return scm_from_bool (0 == cmp
);
3286 else if (SCM_COMPLEXP (y
))
3289 if (0.0 != SCM_COMPLEX_IMAG (y
))
3291 if (xisnan (SCM_COMPLEX_REAL (y
)))
3293 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3294 scm_remember_upto_here_1 (x
);
3295 return scm_from_bool (0 == cmp
);
3297 else if (SCM_FRACTIONP (y
))
3300 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3302 else if (SCM_REALP (x
))
3304 double xx
= SCM_REAL_VALUE (x
);
3305 if (SCM_I_INUMP (y
))
3307 /* see comments with inum/real above */
3308 long yy
= SCM_I_INUM (y
);
3309 return scm_from_bool (xx
== (double) yy
3310 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3311 || (long) xx
== yy
));
3313 else if (SCM_BIGP (y
))
3316 if (xisnan (SCM_REAL_VALUE (x
)))
3318 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3319 scm_remember_upto_here_1 (y
);
3320 return scm_from_bool (0 == cmp
);
3322 else if (SCM_REALP (y
))
3323 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3324 else if (SCM_COMPLEXP (y
))
3325 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3326 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3327 else if (SCM_FRACTIONP (y
))
3329 double xx
= SCM_REAL_VALUE (x
);
3333 return scm_from_bool (xx
< 0.0);
3334 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3338 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3340 else if (SCM_COMPLEXP (x
))
3342 if (SCM_I_INUMP (y
))
3343 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3344 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3345 else if (SCM_BIGP (y
))
3348 if (0.0 != SCM_COMPLEX_IMAG (x
))
3350 if (xisnan (SCM_COMPLEX_REAL (x
)))
3352 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3353 scm_remember_upto_here_1 (y
);
3354 return scm_from_bool (0 == cmp
);
3356 else if (SCM_REALP (y
))
3357 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3358 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3359 else if (SCM_COMPLEXP (y
))
3360 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3361 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3362 else if (SCM_FRACTIONP (y
))
3365 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3367 xx
= SCM_COMPLEX_REAL (x
);
3371 return scm_from_bool (xx
< 0.0);
3372 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3376 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3378 else if (SCM_FRACTIONP (x
))
3380 if (SCM_I_INUMP (y
))
3382 else if (SCM_BIGP (y
))
3384 else if (SCM_REALP (y
))
3386 double yy
= SCM_REAL_VALUE (y
);
3390 return scm_from_bool (0.0 < yy
);
3391 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3394 else if (SCM_COMPLEXP (y
))
3397 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3399 yy
= SCM_COMPLEX_REAL (y
);
3403 return scm_from_bool (0.0 < yy
);
3404 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3407 else if (SCM_FRACTIONP (y
))
3408 return scm_i_fraction_equalp (x
, y
);
3410 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3413 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3417 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3418 done are good for inums, but for bignums an answer can almost always be
3419 had by just examining a few high bits of the operands, as done by GMP in
3420 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3421 of the float exponent to take into account. */
3423 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3424 /* "Return @code{#t} if the list of parameters is monotonically\n"
3428 scm_less_p (SCM x
, SCM y
)
3431 if (SCM_I_INUMP (x
))
3433 long xx
= SCM_I_INUM (x
);
3434 if (SCM_I_INUMP (y
))
3436 long yy
= SCM_I_INUM (y
);
3437 return scm_from_bool (xx
< yy
);
3439 else if (SCM_BIGP (y
))
3441 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3442 scm_remember_upto_here_1 (y
);
3443 return scm_from_bool (sgn
> 0);
3445 else if (SCM_REALP (y
))
3446 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3447 else if (SCM_FRACTIONP (y
))
3449 /* "x < a/b" becomes "x*b < a" */
3451 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3452 y
= SCM_FRACTION_NUMERATOR (y
);
3456 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3458 else if (SCM_BIGP (x
))
3460 if (SCM_I_INUMP (y
))
3462 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3463 scm_remember_upto_here_1 (x
);
3464 return scm_from_bool (sgn
< 0);
3466 else if (SCM_BIGP (y
))
3468 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3469 scm_remember_upto_here_2 (x
, y
);
3470 return scm_from_bool (cmp
< 0);
3472 else if (SCM_REALP (y
))
3475 if (xisnan (SCM_REAL_VALUE (y
)))
3477 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3478 scm_remember_upto_here_1 (x
);
3479 return scm_from_bool (cmp
< 0);
3481 else if (SCM_FRACTIONP (y
))
3484 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3486 else if (SCM_REALP (x
))
3488 if (SCM_I_INUMP (y
))
3489 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3490 else if (SCM_BIGP (y
))
3493 if (xisnan (SCM_REAL_VALUE (x
)))
3495 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3496 scm_remember_upto_here_1 (y
);
3497 return scm_from_bool (cmp
> 0);
3499 else if (SCM_REALP (y
))
3500 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3501 else if (SCM_FRACTIONP (y
))
3503 double xx
= SCM_REAL_VALUE (x
);
3507 return scm_from_bool (xx
< 0.0);
3508 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3512 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3514 else if (SCM_FRACTIONP (x
))
3516 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3518 /* "a/b < y" becomes "a < y*b" */
3519 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3520 x
= SCM_FRACTION_NUMERATOR (x
);
3523 else if (SCM_REALP (y
))
3525 double yy
= SCM_REAL_VALUE (y
);
3529 return scm_from_bool (0.0 < yy
);
3530 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3533 else if (SCM_FRACTIONP (y
))
3535 /* "a/b < c/d" becomes "a*d < c*b" */
3536 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3537 SCM_FRACTION_DENOMINATOR (y
));
3538 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3539 SCM_FRACTION_DENOMINATOR (x
));
3545 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3548 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3552 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3553 /* "Return @code{#t} if the list of parameters is monotonically\n"
3556 #define FUNC_NAME s_scm_gr_p
3558 scm_gr_p (SCM x
, SCM y
)
3560 if (!SCM_NUMBERP (x
))
3561 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3562 else if (!SCM_NUMBERP (y
))
3563 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3565 return scm_less_p (y
, x
);
3570 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3571 /* "Return @code{#t} if the list of parameters is monotonically\n"
3574 #define FUNC_NAME s_scm_leq_p
3576 scm_leq_p (SCM x
, SCM y
)
3578 if (!SCM_NUMBERP (x
))
3579 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3580 else if (!SCM_NUMBERP (y
))
3581 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3582 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3585 return scm_not (scm_less_p (y
, x
));
3590 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3591 /* "Return @code{#t} if the list of parameters is monotonically\n"
3594 #define FUNC_NAME s_scm_geq_p
3596 scm_geq_p (SCM x
, SCM y
)
3598 if (!SCM_NUMBERP (x
))
3599 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3600 else if (!SCM_NUMBERP (y
))
3601 SCM_WTA_DISPATCH_2 (g_geq_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 (x
, y
));
3610 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3611 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3617 if (SCM_I_INUMP (z
))
3618 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3619 else if (SCM_BIGP (z
))
3621 else if (SCM_REALP (z
))
3622 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3623 else if (SCM_COMPLEXP (z
))
3624 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3625 && SCM_COMPLEX_IMAG (z
) == 0.0);
3626 else if (SCM_FRACTIONP (z
))
3629 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3633 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3634 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3638 scm_positive_p (SCM x
)
3640 if (SCM_I_INUMP (x
))
3641 return scm_from_bool (SCM_I_INUM (x
) > 0);
3642 else if (SCM_BIGP (x
))
3644 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3645 scm_remember_upto_here_1 (x
);
3646 return scm_from_bool (sgn
> 0);
3648 else if (SCM_REALP (x
))
3649 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3650 else if (SCM_FRACTIONP (x
))
3651 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3653 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3657 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3658 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3662 scm_negative_p (SCM x
)
3664 if (SCM_I_INUMP (x
))
3665 return scm_from_bool (SCM_I_INUM (x
) < 0);
3666 else if (SCM_BIGP (x
))
3668 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3669 scm_remember_upto_here_1 (x
);
3670 return scm_from_bool (sgn
< 0);
3672 else if (SCM_REALP (x
))
3673 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3674 else if (SCM_FRACTIONP (x
))
3675 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3677 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3681 /* scm_min and scm_max return an inexact when either argument is inexact, as
3682 required by r5rs. On that basis, for exact/inexact combinations the
3683 exact is converted to inexact to compare and possibly return. This is
3684 unlike scm_less_p above which takes some trouble to preserve all bits in
3685 its test, such trouble is not required for min and max. */
3687 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3688 /* "Return the maximum of all parameter values."
3691 scm_max (SCM x
, SCM y
)
3696 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3697 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3700 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3703 if (SCM_I_INUMP (x
))
3705 long xx
= SCM_I_INUM (x
);
3706 if (SCM_I_INUMP (y
))
3708 long yy
= SCM_I_INUM (y
);
3709 return (xx
< yy
) ? y
: x
;
3711 else if (SCM_BIGP (y
))
3713 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3714 scm_remember_upto_here_1 (y
);
3715 return (sgn
< 0) ? x
: y
;
3717 else if (SCM_REALP (y
))
3720 /* if y==NaN then ">" is false and we return NaN */
3721 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3723 else if (SCM_FRACTIONP (y
))
3726 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3729 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3731 else if (SCM_BIGP (x
))
3733 if (SCM_I_INUMP (y
))
3735 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3736 scm_remember_upto_here_1 (x
);
3737 return (sgn
< 0) ? y
: x
;
3739 else if (SCM_BIGP (y
))
3741 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3742 scm_remember_upto_here_2 (x
, y
);
3743 return (cmp
> 0) ? x
: y
;
3745 else if (SCM_REALP (y
))
3747 /* if y==NaN then xx>yy is false, so we return the NaN y */
3750 xx
= scm_i_big2dbl (x
);
3751 yy
= SCM_REAL_VALUE (y
);
3752 return (xx
> yy
? scm_from_double (xx
) : y
);
3754 else if (SCM_FRACTIONP (y
))
3759 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3761 else if (SCM_REALP (x
))
3763 if (SCM_I_INUMP (y
))
3765 double z
= SCM_I_INUM (y
);
3766 /* if x==NaN then "<" is false and we return NaN */
3767 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3769 else if (SCM_BIGP (y
))
3774 else if (SCM_REALP (y
))
3776 /* if x==NaN then our explicit check means we return NaN
3777 if y==NaN then ">" is false and we return NaN
3778 calling isnan is unavoidable, since it's the only way to know
3779 which of x or y causes any compares to be false */
3780 double xx
= SCM_REAL_VALUE (x
);
3781 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3783 else if (SCM_FRACTIONP (y
))
3785 double yy
= scm_i_fraction2double (y
);
3786 double xx
= SCM_REAL_VALUE (x
);
3787 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3790 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3792 else if (SCM_FRACTIONP (x
))
3794 if (SCM_I_INUMP (y
))
3798 else if (SCM_BIGP (y
))
3802 else if (SCM_REALP (y
))
3804 double xx
= scm_i_fraction2double (x
);
3805 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3807 else if (SCM_FRACTIONP (y
))
3812 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3815 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3819 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3820 /* "Return the minium of all parameter values."
3823 scm_min (SCM x
, SCM y
)
3828 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3829 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3832 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3835 if (SCM_I_INUMP (x
))
3837 long xx
= SCM_I_INUM (x
);
3838 if (SCM_I_INUMP (y
))
3840 long yy
= SCM_I_INUM (y
);
3841 return (xx
< yy
) ? x
: y
;
3843 else if (SCM_BIGP (y
))
3845 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3846 scm_remember_upto_here_1 (y
);
3847 return (sgn
< 0) ? y
: x
;
3849 else if (SCM_REALP (y
))
3852 /* if y==NaN then "<" is false and we return NaN */
3853 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3855 else if (SCM_FRACTIONP (y
))
3858 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3861 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3863 else if (SCM_BIGP (x
))
3865 if (SCM_I_INUMP (y
))
3867 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3868 scm_remember_upto_here_1 (x
);
3869 return (sgn
< 0) ? x
: y
;
3871 else if (SCM_BIGP (y
))
3873 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3874 scm_remember_upto_here_2 (x
, y
);
3875 return (cmp
> 0) ? y
: x
;
3877 else if (SCM_REALP (y
))
3879 /* if y==NaN then xx<yy is false, so we return the NaN y */
3882 xx
= scm_i_big2dbl (x
);
3883 yy
= SCM_REAL_VALUE (y
);
3884 return (xx
< yy
? scm_from_double (xx
) : y
);
3886 else if (SCM_FRACTIONP (y
))
3891 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3893 else if (SCM_REALP (x
))
3895 if (SCM_I_INUMP (y
))
3897 double z
= SCM_I_INUM (y
);
3898 /* if x==NaN then "<" is false and we return NaN */
3899 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3901 else if (SCM_BIGP (y
))
3906 else if (SCM_REALP (y
))
3908 /* if x==NaN then our explicit check means we return NaN
3909 if y==NaN then "<" is false and we return NaN
3910 calling isnan is unavoidable, since it's the only way to know
3911 which of x or y causes any compares to be false */
3912 double xx
= SCM_REAL_VALUE (x
);
3913 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3915 else if (SCM_FRACTIONP (y
))
3917 double yy
= scm_i_fraction2double (y
);
3918 double xx
= SCM_REAL_VALUE (x
);
3919 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3922 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3924 else if (SCM_FRACTIONP (x
))
3926 if (SCM_I_INUMP (y
))
3930 else if (SCM_BIGP (y
))
3934 else if (SCM_REALP (y
))
3936 double xx
= scm_i_fraction2double (x
);
3937 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3939 else if (SCM_FRACTIONP (y
))
3944 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3947 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3951 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3952 /* "Return the sum of all parameter values. Return 0 if called without\n"
3956 scm_sum (SCM x
, SCM y
)
3958 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
3960 if (SCM_NUMBERP (x
)) return x
;
3961 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3962 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3965 if (SCM_LIKELY (SCM_I_INUMP (x
)))
3967 if (SCM_LIKELY (SCM_I_INUMP (y
)))
3969 long xx
= SCM_I_INUM (x
);
3970 long yy
= SCM_I_INUM (y
);
3971 long int z
= xx
+ yy
;
3972 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3974 else if (SCM_BIGP (y
))
3979 else if (SCM_REALP (y
))
3981 long int xx
= SCM_I_INUM (x
);
3982 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3984 else if (SCM_COMPLEXP (y
))
3986 long int xx
= SCM_I_INUM (x
);
3987 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3988 SCM_COMPLEX_IMAG (y
));
3990 else if (SCM_FRACTIONP (y
))
3991 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3992 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3993 SCM_FRACTION_DENOMINATOR (y
));
3995 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3996 } else if (SCM_BIGP (x
))
3998 if (SCM_I_INUMP (y
))
4003 inum
= SCM_I_INUM (y
);
4006 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4009 SCM result
= scm_i_mkbig ();
4010 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4011 scm_remember_upto_here_1 (x
);
4012 /* we know the result will have to be a bignum */
4015 return scm_i_normbig (result
);
4019 SCM result
= scm_i_mkbig ();
4020 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4021 scm_remember_upto_here_1 (x
);
4022 /* we know the result will have to be a bignum */
4025 return scm_i_normbig (result
);
4028 else if (SCM_BIGP (y
))
4030 SCM result
= scm_i_mkbig ();
4031 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4032 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4033 mpz_add (SCM_I_BIG_MPZ (result
),
4036 scm_remember_upto_here_2 (x
, y
);
4037 /* we know the result will have to be a bignum */
4040 return scm_i_normbig (result
);
4042 else if (SCM_REALP (y
))
4044 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4045 scm_remember_upto_here_1 (x
);
4046 return scm_from_double (result
);
4048 else if (SCM_COMPLEXP (y
))
4050 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4051 + SCM_COMPLEX_REAL (y
));
4052 scm_remember_upto_here_1 (x
);
4053 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4055 else if (SCM_FRACTIONP (y
))
4056 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4057 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4058 SCM_FRACTION_DENOMINATOR (y
));
4060 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4062 else if (SCM_REALP (x
))
4064 if (SCM_I_INUMP (y
))
4065 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4066 else if (SCM_BIGP (y
))
4068 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4069 scm_remember_upto_here_1 (y
);
4070 return scm_from_double (result
);
4072 else if (SCM_REALP (y
))
4073 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4074 else if (SCM_COMPLEXP (y
))
4075 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4076 SCM_COMPLEX_IMAG (y
));
4077 else if (SCM_FRACTIONP (y
))
4078 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4080 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4082 else if (SCM_COMPLEXP (x
))
4084 if (SCM_I_INUMP (y
))
4085 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4086 SCM_COMPLEX_IMAG (x
));
4087 else if (SCM_BIGP (y
))
4089 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4090 + SCM_COMPLEX_REAL (x
));
4091 scm_remember_upto_here_1 (y
);
4092 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4094 else if (SCM_REALP (y
))
4095 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4096 SCM_COMPLEX_IMAG (x
));
4097 else if (SCM_COMPLEXP (y
))
4098 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4099 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4100 else if (SCM_FRACTIONP (y
))
4101 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4102 SCM_COMPLEX_IMAG (x
));
4104 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4106 else if (SCM_FRACTIONP (x
))
4108 if (SCM_I_INUMP (y
))
4109 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4110 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4111 SCM_FRACTION_DENOMINATOR (x
));
4112 else if (SCM_BIGP (y
))
4113 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4114 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4115 SCM_FRACTION_DENOMINATOR (x
));
4116 else if (SCM_REALP (y
))
4117 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4118 else if (SCM_COMPLEXP (y
))
4119 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4120 SCM_COMPLEX_IMAG (y
));
4121 else if (SCM_FRACTIONP (y
))
4122 /* a/b + c/d = (ad + bc) / bd */
4123 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4124 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4125 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4127 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4130 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4134 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4136 "Return @math{@var{x}+1}.")
4137 #define FUNC_NAME s_scm_oneplus
4139 return scm_sum (x
, SCM_I_MAKINUM (1));
4144 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4145 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4146 * the sum of all but the first argument are subtracted from the first
4148 #define FUNC_NAME s_difference
4150 scm_difference (SCM x
, SCM y
)
4152 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4155 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4157 if (SCM_I_INUMP (x
))
4159 long xx
= -SCM_I_INUM (x
);
4160 if (SCM_FIXABLE (xx
))
4161 return SCM_I_MAKINUM (xx
);
4163 return scm_i_long2big (xx
);
4165 else if (SCM_BIGP (x
))
4166 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4167 bignum, but negating that gives a fixnum. */
4168 return scm_i_normbig (scm_i_clonebig (x
, 0));
4169 else if (SCM_REALP (x
))
4170 return scm_from_double (-SCM_REAL_VALUE (x
));
4171 else if (SCM_COMPLEXP (x
))
4172 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4173 -SCM_COMPLEX_IMAG (x
));
4174 else if (SCM_FRACTIONP (x
))
4175 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4176 SCM_FRACTION_DENOMINATOR (x
));
4178 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4181 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4183 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4185 long int xx
= SCM_I_INUM (x
);
4186 long int yy
= SCM_I_INUM (y
);
4187 long int z
= xx
- yy
;
4188 if (SCM_FIXABLE (z
))
4189 return SCM_I_MAKINUM (z
);
4191 return scm_i_long2big (z
);
4193 else if (SCM_BIGP (y
))
4195 /* inum-x - big-y */
4196 long xx
= SCM_I_INUM (x
);
4199 return scm_i_clonebig (y
, 0);
4202 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4203 SCM result
= scm_i_mkbig ();
4206 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4209 /* x - y == -(y + -x) */
4210 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4211 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4213 scm_remember_upto_here_1 (y
);
4215 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4216 /* we know the result will have to be a bignum */
4219 return scm_i_normbig (result
);
4222 else if (SCM_REALP (y
))
4224 long int xx
= SCM_I_INUM (x
);
4225 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4227 else if (SCM_COMPLEXP (y
))
4229 long int xx
= SCM_I_INUM (x
);
4230 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4231 - SCM_COMPLEX_IMAG (y
));
4233 else if (SCM_FRACTIONP (y
))
4234 /* a - b/c = (ac - b) / c */
4235 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4236 SCM_FRACTION_NUMERATOR (y
)),
4237 SCM_FRACTION_DENOMINATOR (y
));
4239 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4241 else if (SCM_BIGP (x
))
4243 if (SCM_I_INUMP (y
))
4245 /* big-x - inum-y */
4246 long yy
= SCM_I_INUM (y
);
4247 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4249 scm_remember_upto_here_1 (x
);
4251 return (SCM_FIXABLE (-yy
) ?
4252 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4255 SCM result
= scm_i_mkbig ();
4258 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4260 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4261 scm_remember_upto_here_1 (x
);
4263 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4264 /* we know the result will have to be a bignum */
4267 return scm_i_normbig (result
);
4270 else if (SCM_BIGP (y
))
4272 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4273 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4274 SCM result
= scm_i_mkbig ();
4275 mpz_sub (SCM_I_BIG_MPZ (result
),
4278 scm_remember_upto_here_2 (x
, y
);
4279 /* we know the result will have to be a bignum */
4280 if ((sgn_x
== 1) && (sgn_y
== -1))
4282 if ((sgn_x
== -1) && (sgn_y
== 1))
4284 return scm_i_normbig (result
);
4286 else if (SCM_REALP (y
))
4288 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4289 scm_remember_upto_here_1 (x
);
4290 return scm_from_double (result
);
4292 else if (SCM_COMPLEXP (y
))
4294 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4295 - SCM_COMPLEX_REAL (y
));
4296 scm_remember_upto_here_1 (x
);
4297 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4299 else if (SCM_FRACTIONP (y
))
4300 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4301 SCM_FRACTION_NUMERATOR (y
)),
4302 SCM_FRACTION_DENOMINATOR (y
));
4303 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4305 else if (SCM_REALP (x
))
4307 if (SCM_I_INUMP (y
))
4308 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4309 else if (SCM_BIGP (y
))
4311 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4312 scm_remember_upto_here_1 (x
);
4313 return scm_from_double (result
);
4315 else if (SCM_REALP (y
))
4316 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4317 else if (SCM_COMPLEXP (y
))
4318 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4319 -SCM_COMPLEX_IMAG (y
));
4320 else if (SCM_FRACTIONP (y
))
4321 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4323 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4325 else if (SCM_COMPLEXP (x
))
4327 if (SCM_I_INUMP (y
))
4328 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4329 SCM_COMPLEX_IMAG (x
));
4330 else if (SCM_BIGP (y
))
4332 double real_part
= (SCM_COMPLEX_REAL (x
)
4333 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4334 scm_remember_upto_here_1 (x
);
4335 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4337 else if (SCM_REALP (y
))
4338 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4339 SCM_COMPLEX_IMAG (x
));
4340 else if (SCM_COMPLEXP (y
))
4341 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4342 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4343 else if (SCM_FRACTIONP (y
))
4344 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4345 SCM_COMPLEX_IMAG (x
));
4347 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4349 else if (SCM_FRACTIONP (x
))
4351 if (SCM_I_INUMP (y
))
4352 /* a/b - c = (a - cb) / b */
4353 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4354 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4355 SCM_FRACTION_DENOMINATOR (x
));
4356 else if (SCM_BIGP (y
))
4357 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4358 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4359 SCM_FRACTION_DENOMINATOR (x
));
4360 else if (SCM_REALP (y
))
4361 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4362 else if (SCM_COMPLEXP (y
))
4363 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4364 -SCM_COMPLEX_IMAG (y
));
4365 else if (SCM_FRACTIONP (y
))
4366 /* a/b - c/d = (ad - bc) / bd */
4367 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4368 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4369 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4371 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4374 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4379 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4381 "Return @math{@var{x}-1}.")
4382 #define FUNC_NAME s_scm_oneminus
4384 return scm_difference (x
, SCM_I_MAKINUM (1));
4389 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4390 /* "Return the product of all arguments. If called without arguments,\n"
4394 scm_product (SCM x
, SCM y
)
4396 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4399 return SCM_I_MAKINUM (1L);
4400 else if (SCM_NUMBERP (x
))
4403 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4406 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4411 xx
= SCM_I_INUM (x
);
4415 case 0: return x
; break;
4416 case 1: return y
; break;
4419 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4421 long yy
= SCM_I_INUM (y
);
4423 SCM k
= SCM_I_MAKINUM (kk
);
4424 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4428 SCM result
= scm_i_long2big (xx
);
4429 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4430 return scm_i_normbig (result
);
4433 else if (SCM_BIGP (y
))
4435 SCM result
= scm_i_mkbig ();
4436 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4437 scm_remember_upto_here_1 (y
);
4440 else if (SCM_REALP (y
))
4441 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4442 else if (SCM_COMPLEXP (y
))
4443 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4444 xx
* SCM_COMPLEX_IMAG (y
));
4445 else if (SCM_FRACTIONP (y
))
4446 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4447 SCM_FRACTION_DENOMINATOR (y
));
4449 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4451 else if (SCM_BIGP (x
))
4453 if (SCM_I_INUMP (y
))
4458 else if (SCM_BIGP (y
))
4460 SCM result
= scm_i_mkbig ();
4461 mpz_mul (SCM_I_BIG_MPZ (result
),
4464 scm_remember_upto_here_2 (x
, y
);
4467 else if (SCM_REALP (y
))
4469 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4470 scm_remember_upto_here_1 (x
);
4471 return scm_from_double (result
);
4473 else if (SCM_COMPLEXP (y
))
4475 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4476 scm_remember_upto_here_1 (x
);
4477 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4478 z
* SCM_COMPLEX_IMAG (y
));
4480 else if (SCM_FRACTIONP (y
))
4481 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4482 SCM_FRACTION_DENOMINATOR (y
));
4484 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4486 else if (SCM_REALP (x
))
4488 if (SCM_I_INUMP (y
))
4490 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4491 if (scm_is_eq (y
, SCM_INUM0
))
4493 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4495 else if (SCM_BIGP (y
))
4497 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4498 scm_remember_upto_here_1 (y
);
4499 return scm_from_double (result
);
4501 else if (SCM_REALP (y
))
4502 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4503 else if (SCM_COMPLEXP (y
))
4504 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4505 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4506 else if (SCM_FRACTIONP (y
))
4507 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4509 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4511 else if (SCM_COMPLEXP (x
))
4513 if (SCM_I_INUMP (y
))
4515 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4516 if (scm_is_eq (y
, SCM_INUM0
))
4518 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4519 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4521 else if (SCM_BIGP (y
))
4523 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4524 scm_remember_upto_here_1 (y
);
4525 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4526 z
* SCM_COMPLEX_IMAG (x
));
4528 else if (SCM_REALP (y
))
4529 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4530 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4531 else if (SCM_COMPLEXP (y
))
4533 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4534 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4535 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4536 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4538 else if (SCM_FRACTIONP (y
))
4540 double yy
= scm_i_fraction2double (y
);
4541 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4542 yy
* SCM_COMPLEX_IMAG (x
));
4545 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4547 else if (SCM_FRACTIONP (x
))
4549 if (SCM_I_INUMP (y
))
4550 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4551 SCM_FRACTION_DENOMINATOR (x
));
4552 else if (SCM_BIGP (y
))
4553 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4554 SCM_FRACTION_DENOMINATOR (x
));
4555 else if (SCM_REALP (y
))
4556 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4557 else if (SCM_COMPLEXP (y
))
4559 double xx
= scm_i_fraction2double (x
);
4560 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4561 xx
* SCM_COMPLEX_IMAG (y
));
4563 else if (SCM_FRACTIONP (y
))
4564 /* a/b * c/d = ac / bd */
4565 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4566 SCM_FRACTION_NUMERATOR (y
)),
4567 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4568 SCM_FRACTION_DENOMINATOR (y
)));
4570 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4573 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4576 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4577 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4578 #define ALLOW_DIVIDE_BY_ZERO
4579 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4582 /* The code below for complex division is adapted from the GNU
4583 libstdc++, which adapted it from f2c's libF77, and is subject to
4586 /****************************************************************
4587 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4589 Permission to use, copy, modify, and distribute this software
4590 and its documentation for any purpose and without fee is hereby
4591 granted, provided that the above copyright notice appear in all
4592 copies and that both that the copyright notice and this
4593 permission notice and warranty disclaimer appear in supporting
4594 documentation, and that the names of AT&T Bell Laboratories or
4595 Bellcore or any of their entities not be used in advertising or
4596 publicity pertaining to distribution of the software without
4597 specific, written prior permission.
4599 AT&T and Bellcore disclaim all warranties with regard to this
4600 software, including all implied warranties of merchantability
4601 and fitness. In no event shall AT&T or Bellcore be liable for
4602 any special, indirect or consequential damages or any damages
4603 whatsoever resulting from loss of use, data or profits, whether
4604 in an action of contract, negligence or other tortious action,
4605 arising out of or in connection with the use or performance of
4607 ****************************************************************/
4609 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4610 /* Divide the first argument by the product of the remaining
4611 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4613 #define FUNC_NAME s_divide
4615 scm_i_divide (SCM x
, SCM y
, int inexact
)
4619 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4622 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4623 else if (SCM_I_INUMP (x
))
4625 long xx
= SCM_I_INUM (x
);
4626 if (xx
== 1 || xx
== -1)
4628 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4630 scm_num_overflow (s_divide
);
4635 return scm_from_double (1.0 / (double) xx
);
4636 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4639 else if (SCM_BIGP (x
))
4642 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4643 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4645 else if (SCM_REALP (x
))
4647 double xx
= SCM_REAL_VALUE (x
);
4648 #ifndef ALLOW_DIVIDE_BY_ZERO
4650 scm_num_overflow (s_divide
);
4653 return scm_from_double (1.0 / xx
);
4655 else if (SCM_COMPLEXP (x
))
4657 double r
= SCM_COMPLEX_REAL (x
);
4658 double i
= SCM_COMPLEX_IMAG (x
);
4659 if (fabs(r
) <= fabs(i
))
4662 double d
= i
* (1.0 + t
* t
);
4663 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4668 double d
= r
* (1.0 + t
* t
);
4669 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4672 else if (SCM_FRACTIONP (x
))
4673 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4674 SCM_FRACTION_NUMERATOR (x
));
4676 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4679 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4681 long xx
= SCM_I_INUM (x
);
4682 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4684 long yy
= SCM_I_INUM (y
);
4687 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4688 scm_num_overflow (s_divide
);
4690 return scm_from_double ((double) xx
/ (double) yy
);
4693 else if (xx
% yy
!= 0)
4696 return scm_from_double ((double) xx
/ (double) yy
);
4697 else return scm_i_make_ratio (x
, y
);
4702 if (SCM_FIXABLE (z
))
4703 return SCM_I_MAKINUM (z
);
4705 return scm_i_long2big (z
);
4708 else if (SCM_BIGP (y
))
4711 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4712 else return scm_i_make_ratio (x
, y
);
4714 else if (SCM_REALP (y
))
4716 double yy
= SCM_REAL_VALUE (y
);
4717 #ifndef ALLOW_DIVIDE_BY_ZERO
4719 scm_num_overflow (s_divide
);
4722 return scm_from_double ((double) xx
/ yy
);
4724 else if (SCM_COMPLEXP (y
))
4727 complex_div
: /* y _must_ be a complex number */
4729 double r
= SCM_COMPLEX_REAL (y
);
4730 double i
= SCM_COMPLEX_IMAG (y
);
4731 if (fabs(r
) <= fabs(i
))
4734 double d
= i
* (1.0 + t
* t
);
4735 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4740 double d
= r
* (1.0 + t
* t
);
4741 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4745 else if (SCM_FRACTIONP (y
))
4746 /* a / b/c = ac / b */
4747 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4748 SCM_FRACTION_NUMERATOR (y
));
4750 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4752 else if (SCM_BIGP (x
))
4754 if (SCM_I_INUMP (y
))
4756 long int yy
= SCM_I_INUM (y
);
4759 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4760 scm_num_overflow (s_divide
);
4762 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4763 scm_remember_upto_here_1 (x
);
4764 return (sgn
== 0) ? scm_nan () : scm_inf ();
4771 /* FIXME: HMM, what are the relative performance issues here?
4772 We need to test. Is it faster on average to test
4773 divisible_p, then perform whichever operation, or is it
4774 faster to perform the integer div opportunistically and
4775 switch to real if there's a remainder? For now we take the
4776 middle ground: test, then if divisible, use the faster div
4779 long abs_yy
= yy
< 0 ? -yy
: yy
;
4780 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4784 SCM result
= scm_i_mkbig ();
4785 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4786 scm_remember_upto_here_1 (x
);
4788 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4789 return scm_i_normbig (result
);
4794 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4795 else return scm_i_make_ratio (x
, y
);
4799 else if (SCM_BIGP (y
))
4801 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4804 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4805 scm_num_overflow (s_divide
);
4807 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4808 scm_remember_upto_here_1 (x
);
4809 return (sgn
== 0) ? scm_nan () : scm_inf ();
4817 /* It's easily possible for the ratio x/y to fit a double
4818 but one or both x and y be too big to fit a double,
4819 hence the use of mpq_get_d rather than converting and
4822 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
4823 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
4824 return scm_from_double (mpq_get_d (q
));
4828 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4832 SCM result
= scm_i_mkbig ();
4833 mpz_divexact (SCM_I_BIG_MPZ (result
),
4836 scm_remember_upto_here_2 (x
, y
);
4837 return scm_i_normbig (result
);
4840 return scm_i_make_ratio (x
, y
);
4844 else if (SCM_REALP (y
))
4846 double yy
= SCM_REAL_VALUE (y
);
4847 #ifndef ALLOW_DIVIDE_BY_ZERO
4849 scm_num_overflow (s_divide
);
4852 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4854 else if (SCM_COMPLEXP (y
))
4856 a
= scm_i_big2dbl (x
);
4859 else if (SCM_FRACTIONP (y
))
4860 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4861 SCM_FRACTION_NUMERATOR (y
));
4863 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4865 else if (SCM_REALP (x
))
4867 double rx
= SCM_REAL_VALUE (x
);
4868 if (SCM_I_INUMP (y
))
4870 long int yy
= SCM_I_INUM (y
);
4871 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4873 scm_num_overflow (s_divide
);
4876 return scm_from_double (rx
/ (double) yy
);
4878 else if (SCM_BIGP (y
))
4880 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4881 scm_remember_upto_here_1 (y
);
4882 return scm_from_double (rx
/ dby
);
4884 else if (SCM_REALP (y
))
4886 double yy
= SCM_REAL_VALUE (y
);
4887 #ifndef ALLOW_DIVIDE_BY_ZERO
4889 scm_num_overflow (s_divide
);
4892 return scm_from_double (rx
/ yy
);
4894 else if (SCM_COMPLEXP (y
))
4899 else if (SCM_FRACTIONP (y
))
4900 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4902 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4904 else if (SCM_COMPLEXP (x
))
4906 double rx
= SCM_COMPLEX_REAL (x
);
4907 double ix
= SCM_COMPLEX_IMAG (x
);
4908 if (SCM_I_INUMP (y
))
4910 long int yy
= SCM_I_INUM (y
);
4911 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4913 scm_num_overflow (s_divide
);
4918 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4921 else if (SCM_BIGP (y
))
4923 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4924 scm_remember_upto_here_1 (y
);
4925 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4927 else if (SCM_REALP (y
))
4929 double yy
= SCM_REAL_VALUE (y
);
4930 #ifndef ALLOW_DIVIDE_BY_ZERO
4932 scm_num_overflow (s_divide
);
4935 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4937 else if (SCM_COMPLEXP (y
))
4939 double ry
= SCM_COMPLEX_REAL (y
);
4940 double iy
= SCM_COMPLEX_IMAG (y
);
4941 if (fabs(ry
) <= fabs(iy
))
4944 double d
= iy
* (1.0 + t
* t
);
4945 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4950 double d
= ry
* (1.0 + t
* t
);
4951 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4954 else if (SCM_FRACTIONP (y
))
4956 double yy
= scm_i_fraction2double (y
);
4957 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4960 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4962 else if (SCM_FRACTIONP (x
))
4964 if (SCM_I_INUMP (y
))
4966 long int yy
= SCM_I_INUM (y
);
4967 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4969 scm_num_overflow (s_divide
);
4972 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4973 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4975 else if (SCM_BIGP (y
))
4977 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4978 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4980 else if (SCM_REALP (y
))
4982 double yy
= SCM_REAL_VALUE (y
);
4983 #ifndef ALLOW_DIVIDE_BY_ZERO
4985 scm_num_overflow (s_divide
);
4988 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4990 else if (SCM_COMPLEXP (y
))
4992 a
= scm_i_fraction2double (x
);
4995 else if (SCM_FRACTIONP (y
))
4996 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4997 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4999 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5002 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5006 scm_divide (SCM x
, SCM y
)
5008 return scm_i_divide (x
, y
, 0);
5011 static SCM
scm_divide2real (SCM x
, SCM y
)
5013 return scm_i_divide (x
, y
, 1);
5019 scm_asinh (double x
)
5024 #define asinh scm_asinh
5025 return log (x
+ sqrt (x
* x
+ 1));
5028 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
5029 /* "Return the inverse hyperbolic sine of @var{x}."
5034 scm_acosh (double x
)
5039 #define acosh scm_acosh
5040 return log (x
+ sqrt (x
* x
- 1));
5043 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
5044 /* "Return the inverse hyperbolic cosine of @var{x}."
5049 scm_atanh (double x
)
5054 #define atanh scm_atanh
5055 return 0.5 * log ((1 + x
) / (1 - x
));
5058 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
5059 /* "Return the inverse hyperbolic tangent of @var{x}."
5064 scm_c_truncate (double x
)
5075 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5076 half-way case (ie. when x is an integer plus 0.5) going upwards.
5077 Then half-way cases are identified and adjusted down if the
5078 round-upwards didn't give the desired even integer.
5080 "plus_half == result" identifies a half-way case. If plus_half, which is
5081 x + 0.5, is an integer then x must be an integer plus 0.5.
5083 An odd "result" value is identified with result/2 != floor(result/2).
5084 This is done with plus_half, since that value is ready for use sooner in
5085 a pipelined cpu, and we're already requiring plus_half == result.
5087 Note however that we need to be careful when x is big and already an
5088 integer. In that case "x+0.5" may round to an adjacent integer, causing
5089 us to return such a value, incorrectly. For instance if the hardware is
5090 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5091 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5092 returned. Or if the hardware is in round-upwards mode, then other bigger
5093 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5094 representable value, 2^128+2^76 (or whatever), again incorrect.
5096 These bad roundings of x+0.5 are avoided by testing at the start whether
5097 x is already an integer. If it is then clearly that's the desired result
5098 already. And if it's not then the exponent must be small enough to allow
5099 an 0.5 to be represented, and hence added without a bad rounding. */
5102 scm_c_round (double x
)
5104 double plus_half
, result
;
5109 plus_half
= x
+ 0.5;
5110 result
= floor (plus_half
);
5111 /* Adjust so that the rounding is towards even. */
5112 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5117 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5119 "Round the number @var{x} towards zero.")
5120 #define FUNC_NAME s_scm_truncate_number
5122 if (scm_is_false (scm_negative_p (x
)))
5123 return scm_floor (x
);
5125 return scm_ceiling (x
);
5129 static SCM exactly_one_half
;
5131 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5133 "Round the number @var{x} towards the nearest integer. "
5134 "When it is exactly halfway between two integers, "
5135 "round towards the even one.")
5136 #define FUNC_NAME s_scm_round_number
5138 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5140 else if (SCM_REALP (x
))
5141 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5144 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5145 single quotient+remainder division then examining to see which way
5146 the rounding should go. */
5147 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5148 SCM result
= scm_floor (plus_half
);
5149 /* Adjust so that the rounding is towards even. */
5150 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5151 && scm_is_true (scm_odd_p (result
)))
5152 return scm_difference (result
, SCM_I_MAKINUM (1));
5159 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5161 "Round the number @var{x} towards minus infinity.")
5162 #define FUNC_NAME s_scm_floor
5164 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5166 else if (SCM_REALP (x
))
5167 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5168 else if (SCM_FRACTIONP (x
))
5170 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5171 SCM_FRACTION_DENOMINATOR (x
));
5172 if (scm_is_false (scm_negative_p (x
)))
5174 /* For positive x, rounding towards zero is correct. */
5179 /* For negative x, we need to return q-1 unless x is an
5180 integer. But fractions are never integer, per our
5182 return scm_difference (q
, SCM_I_MAKINUM (1));
5186 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5190 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5192 "Round the number @var{x} towards infinity.")
5193 #define FUNC_NAME s_scm_ceiling
5195 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5197 else if (SCM_REALP (x
))
5198 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5199 else if (SCM_FRACTIONP (x
))
5201 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5202 SCM_FRACTION_DENOMINATOR (x
));
5203 if (scm_is_false (scm_positive_p (x
)))
5205 /* For negative x, rounding towards zero is correct. */
5210 /* For positive x, we need to return q+1 unless x is an
5211 integer. But fractions are never integer, per our
5213 return scm_sum (q
, SCM_I_MAKINUM (1));
5217 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5221 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5222 /* "Return the square root of the real number @var{x}."
5224 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5225 /* "Return the absolute value of the real number @var{x}."
5227 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5228 /* "Return the @var{x}th power of e."
5230 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5231 /* "Return the natural logarithm of the real number @var{x}."
5233 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5234 /* "Return the sine of the real number @var{x}."
5236 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5237 /* "Return the cosine of the real number @var{x}."
5239 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5240 /* "Return the tangent of the real number @var{x}."
5242 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5243 /* "Return the arc sine of the real number @var{x}."
5245 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5246 /* "Return the arc cosine of the real number @var{x}."
5248 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5249 /* "Return the arc tangent of the real number @var{x}."
5251 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5252 /* "Return the hyperbolic sine of the real number @var{x}."
5254 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5255 /* "Return the hyperbolic cosine of the real number @var{x}."
5257 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5258 /* "Return the hyperbolic tangent of the real number @var{x}."
5266 static void scm_two_doubles (SCM x
,
5268 const char *sstring
,
5272 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5274 if (SCM_I_INUMP (x
))
5275 xy
->x
= SCM_I_INUM (x
);
5276 else if (SCM_BIGP (x
))
5277 xy
->x
= scm_i_big2dbl (x
);
5278 else if (SCM_REALP (x
))
5279 xy
->x
= SCM_REAL_VALUE (x
);
5280 else if (SCM_FRACTIONP (x
))
5281 xy
->x
= scm_i_fraction2double (x
);
5283 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5285 if (SCM_I_INUMP (y
))
5286 xy
->y
= SCM_I_INUM (y
);
5287 else if (SCM_BIGP (y
))
5288 xy
->y
= scm_i_big2dbl (y
);
5289 else if (SCM_REALP (y
))
5290 xy
->y
= SCM_REAL_VALUE (y
);
5291 else if (SCM_FRACTIONP (y
))
5292 xy
->y
= scm_i_fraction2double (y
);
5294 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5298 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5300 "Return @var{x} raised to the power of @var{y}. This\n"
5301 "procedure does not accept complex arguments.")
5302 #define FUNC_NAME s_scm_sys_expt
5305 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5306 return scm_from_double (pow (xy
.x
, xy
.y
));
5311 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5313 "Return the arc tangent of the two arguments @var{x} and\n"
5314 "@var{y}. This is similar to calculating the arc tangent of\n"
5315 "@var{x} / @var{y}, except that the signs of both arguments\n"
5316 "are used to determine the quadrant of the result. This\n"
5317 "procedure does not accept complex arguments.")
5318 #define FUNC_NAME s_scm_sys_atan2
5321 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5322 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5327 scm_c_make_rectangular (double re
, double im
)
5330 return scm_from_double (re
);
5334 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5336 SCM_COMPLEX_REAL (z
) = re
;
5337 SCM_COMPLEX_IMAG (z
) = im
;
5342 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5343 (SCM real_part
, SCM imaginary_part
),
5344 "Return a complex number constructed of the given @var{real-part} "
5345 "and @var{imaginary-part} parts.")
5346 #define FUNC_NAME s_scm_make_rectangular
5349 scm_two_doubles (real_part
, imaginary_part
, FUNC_NAME
, &xy
);
5350 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5355 scm_c_make_polar (double mag
, double ang
)
5359 sincos (ang
, &s
, &c
);
5364 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5367 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5369 "Return the complex number @var{x} * e^(i * @var{y}).")
5370 #define FUNC_NAME s_scm_make_polar
5373 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5374 return scm_c_make_polar (xy
.x
, xy
.y
);
5379 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5380 /* "Return the real part of the number @var{z}."
5383 scm_real_part (SCM z
)
5385 if (SCM_I_INUMP (z
))
5387 else if (SCM_BIGP (z
))
5389 else if (SCM_REALP (z
))
5391 else if (SCM_COMPLEXP (z
))
5392 return scm_from_double (SCM_COMPLEX_REAL (z
));
5393 else if (SCM_FRACTIONP (z
))
5396 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5400 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5401 /* "Return the imaginary part of the number @var{z}."
5404 scm_imag_part (SCM z
)
5406 if (SCM_I_INUMP (z
))
5408 else if (SCM_BIGP (z
))
5410 else if (SCM_REALP (z
))
5412 else if (SCM_COMPLEXP (z
))
5413 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5414 else if (SCM_FRACTIONP (z
))
5417 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5420 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5421 /* "Return the numerator of the number @var{z}."
5424 scm_numerator (SCM z
)
5426 if (SCM_I_INUMP (z
))
5428 else if (SCM_BIGP (z
))
5430 else if (SCM_FRACTIONP (z
))
5431 return SCM_FRACTION_NUMERATOR (z
);
5432 else if (SCM_REALP (z
))
5433 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5435 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5439 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5440 /* "Return the denominator of the number @var{z}."
5443 scm_denominator (SCM z
)
5445 if (SCM_I_INUMP (z
))
5446 return SCM_I_MAKINUM (1);
5447 else if (SCM_BIGP (z
))
5448 return SCM_I_MAKINUM (1);
5449 else if (SCM_FRACTIONP (z
))
5450 return SCM_FRACTION_DENOMINATOR (z
);
5451 else if (SCM_REALP (z
))
5452 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5454 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5457 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5458 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5459 * "@code{abs} for real arguments, but also allows complex numbers."
5462 scm_magnitude (SCM z
)
5464 if (SCM_I_INUMP (z
))
5466 long int zz
= SCM_I_INUM (z
);
5469 else if (SCM_POSFIXABLE (-zz
))
5470 return SCM_I_MAKINUM (-zz
);
5472 return scm_i_long2big (-zz
);
5474 else if (SCM_BIGP (z
))
5476 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5477 scm_remember_upto_here_1 (z
);
5479 return scm_i_clonebig (z
, 0);
5483 else if (SCM_REALP (z
))
5484 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5485 else if (SCM_COMPLEXP (z
))
5486 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5487 else if (SCM_FRACTIONP (z
))
5489 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5491 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5492 SCM_FRACTION_DENOMINATOR (z
));
5495 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5499 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5500 /* "Return the angle of the complex number @var{z}."
5505 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5506 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5507 But if atan2 follows the floating point rounding mode, then the value
5508 is not a constant. Maybe it'd be close enough though. */
5509 if (SCM_I_INUMP (z
))
5511 if (SCM_I_INUM (z
) >= 0)
5514 return scm_from_double (atan2 (0.0, -1.0));
5516 else if (SCM_BIGP (z
))
5518 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5519 scm_remember_upto_here_1 (z
);
5521 return scm_from_double (atan2 (0.0, -1.0));
5525 else if (SCM_REALP (z
))
5527 if (SCM_REAL_VALUE (z
) >= 0)
5530 return scm_from_double (atan2 (0.0, -1.0));
5532 else if (SCM_COMPLEXP (z
))
5533 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5534 else if (SCM_FRACTIONP (z
))
5536 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5538 else return scm_from_double (atan2 (0.0, -1.0));
5541 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5545 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5546 /* Convert the number @var{x} to its inexact representation.\n"
5549 scm_exact_to_inexact (SCM z
)
5551 if (SCM_I_INUMP (z
))
5552 return scm_from_double ((double) SCM_I_INUM (z
));
5553 else if (SCM_BIGP (z
))
5554 return scm_from_double (scm_i_big2dbl (z
));
5555 else if (SCM_FRACTIONP (z
))
5556 return scm_from_double (scm_i_fraction2double (z
));
5557 else if (SCM_INEXACTP (z
))
5560 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5564 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5566 "Return an exact number that is numerically closest to @var{z}.")
5567 #define FUNC_NAME s_scm_inexact_to_exact
5569 if (SCM_I_INUMP (z
))
5571 else if (SCM_BIGP (z
))
5573 else if (SCM_REALP (z
))
5575 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5576 SCM_OUT_OF_RANGE (1, z
);
5583 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5584 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5585 scm_i_mpz2num (mpq_denref (frac
)));
5587 /* When scm_i_make_ratio throws, we leak the memory allocated
5594 else if (SCM_FRACTIONP (z
))
5597 SCM_WRONG_TYPE_ARG (1, z
);
5601 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5603 "Return an exact number that is within @var{err} of @var{x}.")
5604 #define FUNC_NAME s_scm_rationalize
5606 if (SCM_I_INUMP (x
))
5608 else if (SCM_BIGP (x
))
5610 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5612 /* Use continued fractions to find closest ratio. All
5613 arithmetic is done with exact numbers.
5616 SCM ex
= scm_inexact_to_exact (x
);
5617 SCM int_part
= scm_floor (ex
);
5618 SCM tt
= SCM_I_MAKINUM (1);
5619 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5620 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5624 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5627 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5628 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5630 /* We stop after a million iterations just to be absolutely sure
5631 that we don't go into an infinite loop. The process normally
5632 converges after less than a dozen iterations.
5635 err
= scm_abs (err
);
5636 while (++i
< 1000000)
5638 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5639 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5640 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5642 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5643 err
))) /* abs(x-a/b) <= err */
5645 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5646 if (scm_is_false (scm_exact_p (x
))
5647 || scm_is_false (scm_exact_p (err
)))
5648 return scm_exact_to_inexact (res
);
5652 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5654 tt
= scm_floor (rx
); /* tt = floor (rx) */
5660 scm_num_overflow (s_scm_rationalize
);
5663 SCM_WRONG_TYPE_ARG (1, x
);
5667 /* conversion functions */
5670 scm_is_integer (SCM val
)
5672 return scm_is_true (scm_integer_p (val
));
5676 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5678 if (SCM_I_INUMP (val
))
5680 scm_t_signed_bits n
= SCM_I_INUM (val
);
5681 return n
>= min
&& n
<= max
;
5683 else if (SCM_BIGP (val
))
5685 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5687 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5689 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5691 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5692 return n
>= min
&& n
<= max
;
5702 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5703 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5706 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5707 SCM_I_BIG_MPZ (val
));
5709 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5721 return n
>= min
&& n
<= max
;
5729 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5731 if (SCM_I_INUMP (val
))
5733 scm_t_signed_bits n
= SCM_I_INUM (val
);
5734 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5736 else if (SCM_BIGP (val
))
5738 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5740 else if (max
<= ULONG_MAX
)
5742 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5744 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5745 return n
>= min
&& n
<= max
;
5755 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5758 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5759 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5762 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5763 SCM_I_BIG_MPZ (val
));
5765 return n
>= min
&& n
<= max
;
5773 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5775 scm_error (scm_out_of_range_key
,
5777 "Value out of range ~S to ~S: ~S",
5778 scm_list_3 (min
, max
, bad_val
),
5779 scm_list_1 (bad_val
));
5782 #define TYPE scm_t_intmax
5783 #define TYPE_MIN min
5784 #define TYPE_MAX max
5785 #define SIZEOF_TYPE 0
5786 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5787 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5788 #include "libguile/conv-integer.i.c"
5790 #define TYPE scm_t_uintmax
5791 #define TYPE_MIN min
5792 #define TYPE_MAX max
5793 #define SIZEOF_TYPE 0
5794 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5795 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5796 #include "libguile/conv-uinteger.i.c"
5798 #define TYPE scm_t_int8
5799 #define TYPE_MIN SCM_T_INT8_MIN
5800 #define TYPE_MAX SCM_T_INT8_MAX
5801 #define SIZEOF_TYPE 1
5802 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5803 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5804 #include "libguile/conv-integer.i.c"
5806 #define TYPE scm_t_uint8
5808 #define TYPE_MAX SCM_T_UINT8_MAX
5809 #define SIZEOF_TYPE 1
5810 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5811 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5812 #include "libguile/conv-uinteger.i.c"
5814 #define TYPE scm_t_int16
5815 #define TYPE_MIN SCM_T_INT16_MIN
5816 #define TYPE_MAX SCM_T_INT16_MAX
5817 #define SIZEOF_TYPE 2
5818 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5819 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5820 #include "libguile/conv-integer.i.c"
5822 #define TYPE scm_t_uint16
5824 #define TYPE_MAX SCM_T_UINT16_MAX
5825 #define SIZEOF_TYPE 2
5826 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5827 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5828 #include "libguile/conv-uinteger.i.c"
5830 #define TYPE scm_t_int32
5831 #define TYPE_MIN SCM_T_INT32_MIN
5832 #define TYPE_MAX SCM_T_INT32_MAX
5833 #define SIZEOF_TYPE 4
5834 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5835 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5836 #include "libguile/conv-integer.i.c"
5838 #define TYPE scm_t_uint32
5840 #define TYPE_MAX SCM_T_UINT32_MAX
5841 #define SIZEOF_TYPE 4
5842 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5843 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5844 #include "libguile/conv-uinteger.i.c"
5846 #if SCM_HAVE_T_INT64
5848 #define TYPE scm_t_int64
5849 #define TYPE_MIN SCM_T_INT64_MIN
5850 #define TYPE_MAX SCM_T_INT64_MAX
5851 #define SIZEOF_TYPE 8
5852 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5853 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5854 #include "libguile/conv-integer.i.c"
5856 #define TYPE scm_t_uint64
5858 #define TYPE_MAX SCM_T_UINT64_MAX
5859 #define SIZEOF_TYPE 8
5860 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5861 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5862 #include "libguile/conv-uinteger.i.c"
5867 scm_to_mpz (SCM val
, mpz_t rop
)
5869 if (SCM_I_INUMP (val
))
5870 mpz_set_si (rop
, SCM_I_INUM (val
));
5871 else if (SCM_BIGP (val
))
5872 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5874 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5878 scm_from_mpz (mpz_t val
)
5880 return scm_i_mpz2num (val
);
5884 scm_is_real (SCM val
)
5886 return scm_is_true (scm_real_p (val
));
5890 scm_is_rational (SCM val
)
5892 return scm_is_true (scm_rational_p (val
));
5896 scm_to_double (SCM val
)
5898 if (SCM_I_INUMP (val
))
5899 return SCM_I_INUM (val
);
5900 else if (SCM_BIGP (val
))
5901 return scm_i_big2dbl (val
);
5902 else if (SCM_FRACTIONP (val
))
5903 return scm_i_fraction2double (val
);
5904 else if (SCM_REALP (val
))
5905 return SCM_REAL_VALUE (val
);
5907 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5911 scm_from_double (double val
)
5913 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5914 SCM_REAL_VALUE (z
) = val
;
5918 #if SCM_ENABLE_DISCOURAGED == 1
5921 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5925 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5929 scm_out_of_range (NULL
, num
);
5932 return scm_to_double (num
);
5936 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5940 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5944 scm_out_of_range (NULL
, num
);
5947 return scm_to_double (num
);
5953 scm_is_complex (SCM val
)
5955 return scm_is_true (scm_complex_p (val
));
5959 scm_c_real_part (SCM z
)
5961 if (SCM_COMPLEXP (z
))
5962 return SCM_COMPLEX_REAL (z
);
5965 /* Use the scm_real_part to get proper error checking and
5968 return scm_to_double (scm_real_part (z
));
5973 scm_c_imag_part (SCM z
)
5975 if (SCM_COMPLEXP (z
))
5976 return SCM_COMPLEX_IMAG (z
);
5979 /* Use the scm_imag_part to get proper error checking and
5980 dispatching. The result will almost always be 0.0, but not
5983 return scm_to_double (scm_imag_part (z
));
5988 scm_c_magnitude (SCM z
)
5990 return scm_to_double (scm_magnitude (z
));
5996 return scm_to_double (scm_angle (z
));
6000 scm_is_number (SCM z
)
6002 return scm_is_true (scm_number_p (z
));
6006 /* In the following functions we dispatch to the real-arg funcs like log()
6007 when we know the arg is real, instead of just handing everything to
6008 clog() for instance. This is in case clog() doesn't optimize for a
6009 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6010 well use it to go straight to the applicable C func. */
6012 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6014 "Return the natural logarithm of @var{z}.")
6015 #define FUNC_NAME s_scm_log
6017 if (SCM_COMPLEXP (z
))
6019 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6020 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6022 double re
= SCM_COMPLEX_REAL (z
);
6023 double im
= SCM_COMPLEX_IMAG (z
);
6024 return scm_c_make_rectangular (log (hypot (re
, im
)),
6030 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6031 although the value itself overflows. */
6032 double re
= scm_to_double (z
);
6033 double l
= log (fabs (re
));
6035 return scm_from_double (l
);
6037 return scm_c_make_rectangular (l
, M_PI
);
6043 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6045 "Return the base 10 logarithm of @var{z}.")
6046 #define FUNC_NAME s_scm_log10
6048 if (SCM_COMPLEXP (z
))
6050 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6051 clog() and a multiply by M_LOG10E, rather than the fallback
6052 log10+hypot+atan2.) */
6053 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6054 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6056 double re
= SCM_COMPLEX_REAL (z
);
6057 double im
= SCM_COMPLEX_IMAG (z
);
6058 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6059 M_LOG10E
* atan2 (im
, re
));
6064 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6065 although the value itself overflows. */
6066 double re
= scm_to_double (z
);
6067 double l
= log10 (fabs (re
));
6069 return scm_from_double (l
);
6071 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6077 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6079 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6080 "base of natural logarithms (2.71828@dots{}).")
6081 #define FUNC_NAME s_scm_exp
6083 if (SCM_COMPLEXP (z
))
6085 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6086 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6088 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6089 SCM_COMPLEX_IMAG (z
));
6094 /* When z is a negative bignum the conversion to double overflows,
6095 giving -infinity, but that's ok, the exp is still 0.0. */
6096 return scm_from_double (exp (scm_to_double (z
)));
6102 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6104 "Return the square root of @var{z}. Of the two possible roots\n"
6105 "(positive and negative), the one with the a positive real part\n"
6106 "is returned, or if that's zero then a positive imaginary part.\n"
6110 "(sqrt 9.0) @result{} 3.0\n"
6111 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6112 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6113 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6115 #define FUNC_NAME s_scm_sqrt
6117 if (SCM_COMPLEXP (x
))
6119 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6120 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6122 double re
= SCM_COMPLEX_REAL (x
);
6123 double im
= SCM_COMPLEX_IMAG (x
);
6124 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6125 0.5 * atan2 (im
, re
));
6130 double xx
= scm_to_double (x
);
6132 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6134 return scm_from_double (sqrt (xx
));
6146 mpz_init_set_si (z_negative_one
, -1);
6148 /* It may be possible to tune the performance of some algorithms by using
6149 * the following constants to avoid the creation of bignums. Please, before
6150 * using these values, remember the two rules of program optimization:
6151 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6152 scm_c_define ("most-positive-fixnum",
6153 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6154 scm_c_define ("most-negative-fixnum",
6155 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6157 scm_add_feature ("complex");
6158 scm_add_feature ("inexact");
6159 scm_flo0
= scm_from_double (0.0);
6161 /* determine floating point precision */
6162 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6164 init_dblprec(&scm_dblprec
[i
-2],i
);
6165 init_fx_radix(fx_per_radix
[i
-2],i
);
6168 /* hard code precision for base 10 if the preprocessor tells us to... */
6169 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6172 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6173 SCM_I_MAKINUM (2)));
6174 #include "libguile/numbers.x"