1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006 Free Software Foundation, Inc.
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 /* General assumptions:
24 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
25 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
26 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
27 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
28 * All objects satisfying SCM_FRACTIONP are never an integer.
33 - see if special casing bignums and reals in integer-exponent when
34 possible (to use mpz_pow and mpf_pow_ui) is faster.
36 - look in to better short-circuiting of common cases in
37 integer-expt and elsewhere.
39 - see if direct mpz operations can help in ash and elsewhere.
43 /* tell glibc (2.3) to give prototype for C99 trunc() */
54 #include "libguile/_scm.h"
55 #include "libguile/feature.h"
56 #include "libguile/ports.h"
57 #include "libguile/root.h"
58 #include "libguile/smob.h"
59 #include "libguile/strings.h"
61 #include "libguile/validate.h"
62 #include "libguile/numbers.h"
63 #include "libguile/deprecation.h"
65 #include "libguile/eq.h"
67 #include "libguile/discouraged.h"
72 Wonder if this might be faster for some of our code? A switch on
73 the numtag would jump directly to the right case, and the
74 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
76 #define SCM_I_NUMTAG_NOTNUM 0
77 #define SCM_I_NUMTAG_INUM 1
78 #define SCM_I_NUMTAG_BIG scm_tc16_big
79 #define SCM_I_NUMTAG_REAL scm_tc16_real
80 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
81 #define SCM_I_NUMTAG(x) \
82 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
83 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
84 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
85 : SCM_I_NUMTAG_NOTNUM)))
87 /* the macro above will not work as is with fractions */
90 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
92 /* FLOBUFLEN is the maximum number of characters neccessary for the
93 * printed or scm_string representation of an inexact number.
95 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
98 #if ! defined (HAVE_ISNAN)
103 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
106 #if ! defined (HAVE_ISINF)
111 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
118 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
119 an explicit check. In some future gmp (don't know what version number),
120 mpz_cmp_d is supposed to do this itself. */
122 #define xmpz_cmp_d(z, d) \
123 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
125 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
128 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
129 isinf. It does have finite and isnan though, hence the use of those.
130 fpclass would be a possibility on that system too. */
134 #if defined (HAVE_ISINF)
136 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
137 return (! (finite (x
) || isnan (x
)));
146 #if defined (HAVE_ISNAN)
155 static mpz_t z_negative_one
;
159 SCM_C_INLINE_KEYWORD SCM
162 /* Return a newly created bignum. */
163 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
164 mpz_init (SCM_I_BIG_MPZ (z
));
168 SCM_C_INLINE_KEYWORD SCM
169 scm_i_long2big (long x
)
171 /* Return a newly created bignum initialized to X. */
172 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
173 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
177 SCM_C_INLINE_KEYWORD SCM
178 scm_i_ulong2big (unsigned long x
)
180 /* Return a newly created bignum initialized to X. */
181 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
182 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
186 SCM_C_INLINE_KEYWORD SCM
187 scm_i_clonebig (SCM src_big
, int same_sign_p
)
189 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
190 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
191 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
193 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
197 SCM_C_INLINE_KEYWORD
int
198 scm_i_bigcmp (SCM x
, SCM y
)
200 /* Return neg if x < y, pos if x > y, and 0 if x == y */
201 /* presume we already know x and y are bignums */
202 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
203 scm_remember_upto_here_2 (x
, y
);
207 SCM_C_INLINE_KEYWORD SCM
208 scm_i_dbl2big (double d
)
210 /* results are only defined if d is an integer */
211 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
212 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
216 /* Convert a integer in double representation to a SCM number. */
218 SCM_C_INLINE_KEYWORD SCM
219 scm_i_dbl2num (double u
)
221 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
222 powers of 2, so there's no rounding when making "double" values
223 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
224 get rounded on a 64-bit machine, hence the "+1".
226 The use of floor() to force to an integer value ensures we get a
227 "numerically closest" value without depending on how a
228 double->long cast or how mpz_set_d will round. For reference,
229 double->long probably follows the hardware rounding mode,
230 mpz_set_d truncates towards zero. */
232 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
233 representable as a double? */
235 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
236 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
237 return SCM_I_MAKINUM ((long) u
);
239 return scm_i_dbl2big (u
);
242 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
243 with R5RS exact->inexact.
245 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
246 (ie. truncate towards zero), then adjust to get the closest double by
247 examining the next lower bit and adding 1 (to the absolute value) if
250 Bignums exactly half way between representable doubles are rounded to the
251 next higher absolute value (ie. away from zero). This seems like an
252 adequate interpretation of R5RS "numerically closest", and it's easier
253 and faster than a full "nearest-even" style.
255 The bit test must be done on the absolute value of the mpz_t, which means
256 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
257 negatives as twos complement.
259 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
260 following the hardware rounding mode, but applied to the absolute value
261 of the mpz_t operand. This is not what we want so we put the high
262 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
263 mpz_get_d is supposed to always truncate towards zero.
265 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
266 is a slowdown. It'd be faster to pick out the relevant high bits with
267 mpz_getlimbn if we could be bothered coding that, and if the new
268 truncating gmp doesn't come out. */
271 scm_i_big2dbl (SCM b
)
276 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
280 /* Current GMP, eg. 4.1.3, force truncation towards zero */
282 if (bits
> DBL_MANT_DIG
)
284 size_t shift
= bits
- DBL_MANT_DIG
;
285 mpz_init2 (tmp
, DBL_MANT_DIG
);
286 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
287 result
= ldexp (mpz_get_d (tmp
), shift
);
292 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
297 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
300 if (bits
> DBL_MANT_DIG
)
302 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
303 /* test bit number "pos" in absolute value */
304 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
305 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
307 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
311 scm_remember_upto_here_1 (b
);
315 SCM_C_INLINE_KEYWORD SCM
316 scm_i_normbig (SCM b
)
318 /* convert a big back to a fixnum if it'll fit */
319 /* presume b is a bignum */
320 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
322 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
323 if (SCM_FIXABLE (val
))
324 b
= SCM_I_MAKINUM (val
);
329 static SCM_C_INLINE_KEYWORD SCM
330 scm_i_mpz2num (mpz_t b
)
332 /* convert a mpz number to a SCM number. */
333 if (mpz_fits_slong_p (b
))
335 long val
= mpz_get_si (b
);
336 if (SCM_FIXABLE (val
))
337 return SCM_I_MAKINUM (val
);
341 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
342 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
347 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
348 static SCM
scm_divide2real (SCM x
, SCM y
);
351 scm_i_make_ratio (SCM numerator
, SCM denominator
)
352 #define FUNC_NAME "make-ratio"
354 /* First make sure the arguments are proper.
356 if (SCM_I_INUMP (denominator
))
358 if (scm_is_eq (denominator
, SCM_INUM0
))
359 scm_num_overflow ("make-ratio");
360 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
365 if (!(SCM_BIGP(denominator
)))
366 SCM_WRONG_TYPE_ARG (2, denominator
);
368 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
369 SCM_WRONG_TYPE_ARG (1, numerator
);
371 /* Then flip signs so that the denominator is positive.
373 if (scm_is_true (scm_negative_p (denominator
)))
375 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
376 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
379 /* Now consider for each of the four fixnum/bignum combinations
380 whether the rational number is really an integer.
382 if (SCM_I_INUMP (numerator
))
384 long x
= SCM_I_INUM (numerator
);
385 if (scm_is_eq (numerator
, SCM_INUM0
))
387 if (SCM_I_INUMP (denominator
))
390 y
= SCM_I_INUM (denominator
);
392 return SCM_I_MAKINUM(1);
394 return SCM_I_MAKINUM (x
/ y
);
398 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
399 of that value for the denominator, as a bignum. Apart from
400 that case, abs(bignum) > abs(inum) so inum/bignum is not an
402 if (x
== SCM_MOST_NEGATIVE_FIXNUM
403 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
404 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
405 return SCM_I_MAKINUM(-1);
408 else if (SCM_BIGP (numerator
))
410 if (SCM_I_INUMP (denominator
))
412 long yy
= SCM_I_INUM (denominator
);
413 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
414 return scm_divide (numerator
, denominator
);
418 if (scm_is_eq (numerator
, denominator
))
419 return SCM_I_MAKINUM(1);
420 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
421 SCM_I_BIG_MPZ (denominator
)))
422 return scm_divide(numerator
, denominator
);
426 /* No, it's a proper fraction.
428 return scm_double_cell (scm_tc16_fraction
,
429 SCM_UNPACK (numerator
),
430 SCM_UNPACK (denominator
), 0);
434 static void scm_i_fraction_reduce (SCM z
)
436 if (!(SCM_FRACTION_REDUCED (z
)))
439 divisor
= scm_gcd (SCM_FRACTION_NUMERATOR (z
), SCM_FRACTION_DENOMINATOR (z
));
440 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
443 SCM_FRACTION_SET_NUMERATOR (z
, scm_divide (SCM_FRACTION_NUMERATOR (z
), divisor
));
444 SCM_FRACTION_SET_DENOMINATOR (z
, scm_divide (SCM_FRACTION_DENOMINATOR (z
), divisor
));
446 SCM_FRACTION_REDUCED_SET (z
);
451 scm_i_fraction2double (SCM z
)
453 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
454 SCM_FRACTION_DENOMINATOR (z
)));
457 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
459 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
461 #define FUNC_NAME s_scm_exact_p
467 if (SCM_FRACTIONP (x
))
471 SCM_WRONG_TYPE_ARG (1, x
);
476 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
478 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
480 #define FUNC_NAME s_scm_odd_p
484 long val
= SCM_I_INUM (n
);
485 return scm_from_bool ((val
& 1L) != 0);
487 else if (SCM_BIGP (n
))
489 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
490 scm_remember_upto_here_1 (n
);
491 return scm_from_bool (odd_p
);
493 else if (scm_is_true (scm_inf_p (n
)))
495 else if (SCM_REALP (n
))
497 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
503 SCM_WRONG_TYPE_ARG (1, n
);
506 SCM_WRONG_TYPE_ARG (1, n
);
511 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
513 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
515 #define FUNC_NAME s_scm_even_p
519 long val
= SCM_I_INUM (n
);
520 return scm_from_bool ((val
& 1L) == 0);
522 else if (SCM_BIGP (n
))
524 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
525 scm_remember_upto_here_1 (n
);
526 return scm_from_bool (even_p
);
528 else if (scm_is_true (scm_inf_p (n
)))
530 else if (SCM_REALP (n
))
532 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
538 SCM_WRONG_TYPE_ARG (1, n
);
541 SCM_WRONG_TYPE_ARG (1, n
);
545 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
547 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
548 "or @samp{-inf.0}, @code{#f} otherwise.")
549 #define FUNC_NAME s_scm_inf_p
552 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
553 else if (SCM_COMPLEXP (x
))
554 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
555 || xisinf (SCM_COMPLEX_IMAG (x
)));
561 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
563 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
565 #define FUNC_NAME s_scm_nan_p
568 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
569 else if (SCM_COMPLEXP (n
))
570 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
571 || xisnan (SCM_COMPLEX_IMAG (n
)));
577 /* Guile's idea of infinity. */
578 static double guile_Inf
;
580 /* Guile's idea of not a number. */
581 static double guile_NaN
;
584 guile_ieee_init (void)
586 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
588 /* Some version of gcc on some old version of Linux used to crash when
589 trying to make Inf and NaN. */
592 /* C99 INFINITY, when available.
593 FIXME: The standard allows for INFINITY to be something that overflows
594 at compile time. We ought to have a configure test to check for that
595 before trying to use it. (But in practice we believe this is not a
596 problem on any system guile is likely to target.) */
597 guile_Inf
= INFINITY
;
600 extern unsigned int DINFINITY
[2];
601 guile_Inf
= (*((double *) (DINFINITY
)));
608 if (guile_Inf
== tmp
)
616 #if defined (HAVE_ISNAN)
619 /* C99 NAN, when available */
624 extern unsigned int DQNAN
[2];
625 guile_NaN
= (*((double *)(DQNAN
)));
628 guile_NaN
= guile_Inf
/ guile_Inf
;
634 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
637 #define FUNC_NAME s_scm_inf
639 static int initialized
= 0;
645 return scm_from_double (guile_Inf
);
649 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
652 #define FUNC_NAME s_scm_nan
654 static int initialized
= 0;
660 return scm_from_double (guile_NaN
);
665 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
667 "Return the absolute value of @var{x}.")
672 long int xx
= SCM_I_INUM (x
);
675 else if (SCM_POSFIXABLE (-xx
))
676 return SCM_I_MAKINUM (-xx
);
678 return scm_i_long2big (-xx
);
680 else if (SCM_BIGP (x
))
682 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
684 return scm_i_clonebig (x
, 0);
688 else if (SCM_REALP (x
))
690 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
691 double xx
= SCM_REAL_VALUE (x
);
693 return scm_from_double (-xx
);
697 else if (SCM_FRACTIONP (x
))
699 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
701 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
702 SCM_FRACTION_DENOMINATOR (x
));
705 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
710 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
711 /* "Return the quotient of the numbers @var{x} and @var{y}."
714 scm_quotient (SCM x
, SCM y
)
718 long xx
= SCM_I_INUM (x
);
721 long yy
= SCM_I_INUM (y
);
723 scm_num_overflow (s_quotient
);
728 return SCM_I_MAKINUM (z
);
730 return scm_i_long2big (z
);
733 else if (SCM_BIGP (y
))
735 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
736 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
737 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
739 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
740 scm_remember_upto_here_1 (y
);
741 return SCM_I_MAKINUM (-1);
744 return SCM_I_MAKINUM (0);
747 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
749 else if (SCM_BIGP (x
))
753 long yy
= SCM_I_INUM (y
);
755 scm_num_overflow (s_quotient
);
760 SCM result
= scm_i_mkbig ();
763 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
766 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
769 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
770 scm_remember_upto_here_1 (x
);
771 return scm_i_normbig (result
);
774 else if (SCM_BIGP (y
))
776 SCM result
= scm_i_mkbig ();
777 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
780 scm_remember_upto_here_2 (x
, y
);
781 return scm_i_normbig (result
);
784 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
787 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
790 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
791 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
793 * "(remainder 13 4) @result{} 1\n"
794 * "(remainder -13 4) @result{} -1\n"
798 scm_remainder (SCM x
, SCM y
)
804 long yy
= SCM_I_INUM (y
);
806 scm_num_overflow (s_remainder
);
809 long z
= SCM_I_INUM (x
) % yy
;
810 return SCM_I_MAKINUM (z
);
813 else if (SCM_BIGP (y
))
815 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
816 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
817 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
819 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
820 scm_remember_upto_here_1 (y
);
821 return SCM_I_MAKINUM (0);
827 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
829 else if (SCM_BIGP (x
))
833 long yy
= SCM_I_INUM (y
);
835 scm_num_overflow (s_remainder
);
838 SCM result
= scm_i_mkbig ();
841 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
842 scm_remember_upto_here_1 (x
);
843 return scm_i_normbig (result
);
846 else if (SCM_BIGP (y
))
848 SCM result
= scm_i_mkbig ();
849 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
852 scm_remember_upto_here_2 (x
, y
);
853 return scm_i_normbig (result
);
856 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
859 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
863 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
864 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
866 * "(modulo 13 4) @result{} 1\n"
867 * "(modulo -13 4) @result{} 3\n"
871 scm_modulo (SCM x
, SCM y
)
875 long xx
= SCM_I_INUM (x
);
878 long yy
= SCM_I_INUM (y
);
880 scm_num_overflow (s_modulo
);
883 /* C99 specifies that "%" is the remainder corresponding to a
884 quotient rounded towards zero, and that's also traditional
885 for machine division, so z here should be well defined. */
903 return SCM_I_MAKINUM (result
);
906 else if (SCM_BIGP (y
))
908 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
915 SCM pos_y
= scm_i_clonebig (y
, 0);
916 /* do this after the last scm_op */
917 mpz_init_set_si (z_x
, xx
);
918 result
= pos_y
; /* re-use this bignum */
919 mpz_mod (SCM_I_BIG_MPZ (result
),
921 SCM_I_BIG_MPZ (pos_y
));
922 scm_remember_upto_here_1 (pos_y
);
926 result
= scm_i_mkbig ();
927 /* do this after the last scm_op */
928 mpz_init_set_si (z_x
, xx
);
929 mpz_mod (SCM_I_BIG_MPZ (result
),
932 scm_remember_upto_here_1 (y
);
935 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
936 mpz_add (SCM_I_BIG_MPZ (result
),
938 SCM_I_BIG_MPZ (result
));
939 scm_remember_upto_here_1 (y
);
940 /* and do this before the next one */
942 return scm_i_normbig (result
);
946 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
948 else if (SCM_BIGP (x
))
952 long yy
= SCM_I_INUM (y
);
954 scm_num_overflow (s_modulo
);
957 SCM result
= scm_i_mkbig ();
958 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
960 (yy
< 0) ? - yy
: yy
);
961 scm_remember_upto_here_1 (x
);
962 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
963 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
964 SCM_I_BIG_MPZ (result
),
966 return scm_i_normbig (result
);
969 else if (SCM_BIGP (y
))
972 SCM result
= scm_i_mkbig ();
973 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
974 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
975 mpz_mod (SCM_I_BIG_MPZ (result
),
977 SCM_I_BIG_MPZ (pos_y
));
979 scm_remember_upto_here_1 (x
);
980 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
981 mpz_add (SCM_I_BIG_MPZ (result
),
983 SCM_I_BIG_MPZ (result
));
984 scm_remember_upto_here_2 (y
, pos_y
);
985 return scm_i_normbig (result
);
989 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
992 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
995 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
996 /* "Return the greatest common divisor of all arguments.\n"
997 * "If called without arguments, 0 is returned."
1000 scm_gcd (SCM x
, SCM y
)
1003 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1005 if (SCM_I_INUMP (x
))
1007 if (SCM_I_INUMP (y
))
1009 long xx
= SCM_I_INUM (x
);
1010 long yy
= SCM_I_INUM (y
);
1011 long u
= xx
< 0 ? -xx
: xx
;
1012 long v
= yy
< 0 ? -yy
: yy
;
1022 /* Determine a common factor 2^k */
1023 while (!(1 & (u
| v
)))
1029 /* Now, any factor 2^n can be eliminated */
1049 return (SCM_POSFIXABLE (result
)
1050 ? SCM_I_MAKINUM (result
)
1051 : scm_i_long2big (result
));
1053 else if (SCM_BIGP (y
))
1059 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1061 else if (SCM_BIGP (x
))
1063 if (SCM_I_INUMP (y
))
1065 unsigned long result
;
1068 yy
= SCM_I_INUM (y
);
1073 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1074 scm_remember_upto_here_1 (x
);
1075 return (SCM_POSFIXABLE (result
)
1076 ? SCM_I_MAKINUM (result
)
1077 : scm_from_ulong (result
));
1079 else if (SCM_BIGP (y
))
1081 SCM result
= scm_i_mkbig ();
1082 mpz_gcd (SCM_I_BIG_MPZ (result
),
1085 scm_remember_upto_here_2 (x
, y
);
1086 return scm_i_normbig (result
);
1089 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1092 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1095 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1096 /* "Return the least common multiple of the arguments.\n"
1097 * "If called without arguments, 1 is returned."
1100 scm_lcm (SCM n1
, SCM n2
)
1102 if (SCM_UNBNDP (n2
))
1104 if (SCM_UNBNDP (n1
))
1105 return SCM_I_MAKINUM (1L);
1106 n2
= SCM_I_MAKINUM (1L);
1109 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1110 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1111 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1112 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1114 if (SCM_I_INUMP (n1
))
1116 if (SCM_I_INUMP (n2
))
1118 SCM d
= scm_gcd (n1
, n2
);
1119 if (scm_is_eq (d
, SCM_INUM0
))
1122 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1126 /* inum n1, big n2 */
1129 SCM result
= scm_i_mkbig ();
1130 long nn1
= SCM_I_INUM (n1
);
1131 if (nn1
== 0) return SCM_INUM0
;
1132 if (nn1
< 0) nn1
= - nn1
;
1133 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1134 scm_remember_upto_here_1 (n2
);
1142 if (SCM_I_INUMP (n2
))
1149 SCM result
= scm_i_mkbig ();
1150 mpz_lcm(SCM_I_BIG_MPZ (result
),
1152 SCM_I_BIG_MPZ (n2
));
1153 scm_remember_upto_here_2(n1
, n2
);
1154 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1160 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1165 + + + x (map digit:logand X Y)
1166 + - + x (map digit:logand X (lognot (+ -1 Y)))
1167 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1168 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1173 + + + (map digit:logior X Y)
1174 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1175 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1176 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1181 + + + (map digit:logxor X Y)
1182 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1183 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1184 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1189 + + (any digit:logand X Y)
1190 + - (any digit:logand X (lognot (+ -1 Y)))
1191 - + (any digit:logand (lognot (+ -1 X)) Y)
1196 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1198 "Return the bitwise AND of the integer arguments.\n\n"
1200 "(logand) @result{} -1\n"
1201 "(logand 7) @result{} 7\n"
1202 "(logand #b111 #b011 #b001) @result{} 1\n"
1204 #define FUNC_NAME s_scm_logand
1208 if (SCM_UNBNDP (n2
))
1210 if (SCM_UNBNDP (n1
))
1211 return SCM_I_MAKINUM (-1);
1212 else if (!SCM_NUMBERP (n1
))
1213 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1214 else if (SCM_NUMBERP (n1
))
1217 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1220 if (SCM_I_INUMP (n1
))
1222 nn1
= SCM_I_INUM (n1
);
1223 if (SCM_I_INUMP (n2
))
1225 long nn2
= SCM_I_INUM (n2
);
1226 return SCM_I_MAKINUM (nn1
& nn2
);
1228 else if SCM_BIGP (n2
)
1234 SCM result_z
= scm_i_mkbig ();
1236 mpz_init_set_si (nn1_z
, nn1
);
1237 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1238 scm_remember_upto_here_1 (n2
);
1240 return scm_i_normbig (result_z
);
1244 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1246 else if (SCM_BIGP (n1
))
1248 if (SCM_I_INUMP (n2
))
1251 nn1
= SCM_I_INUM (n1
);
1254 else if (SCM_BIGP (n2
))
1256 SCM result_z
= scm_i_mkbig ();
1257 mpz_and (SCM_I_BIG_MPZ (result_z
),
1259 SCM_I_BIG_MPZ (n2
));
1260 scm_remember_upto_here_2 (n1
, n2
);
1261 return scm_i_normbig (result_z
);
1264 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1267 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1272 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1274 "Return the bitwise OR of the integer arguments.\n\n"
1276 "(logior) @result{} 0\n"
1277 "(logior 7) @result{} 7\n"
1278 "(logior #b000 #b001 #b011) @result{} 3\n"
1280 #define FUNC_NAME s_scm_logior
1284 if (SCM_UNBNDP (n2
))
1286 if (SCM_UNBNDP (n1
))
1288 else if (SCM_NUMBERP (n1
))
1291 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1294 if (SCM_I_INUMP (n1
))
1296 nn1
= SCM_I_INUM (n1
);
1297 if (SCM_I_INUMP (n2
))
1299 long nn2
= SCM_I_INUM (n2
);
1300 return SCM_I_MAKINUM (nn1
| nn2
);
1302 else if (SCM_BIGP (n2
))
1308 SCM result_z
= scm_i_mkbig ();
1310 mpz_init_set_si (nn1_z
, nn1
);
1311 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1312 scm_remember_upto_here_1 (n2
);
1314 return scm_i_normbig (result_z
);
1318 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1320 else if (SCM_BIGP (n1
))
1322 if (SCM_I_INUMP (n2
))
1325 nn1
= SCM_I_INUM (n1
);
1328 else if (SCM_BIGP (n2
))
1330 SCM result_z
= scm_i_mkbig ();
1331 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1333 SCM_I_BIG_MPZ (n2
));
1334 scm_remember_upto_here_2 (n1
, n2
);
1335 return scm_i_normbig (result_z
);
1338 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1341 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1346 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1348 "Return the bitwise XOR of the integer arguments. A bit is\n"
1349 "set in the result if it is set in an odd number of arguments.\n"
1351 "(logxor) @result{} 0\n"
1352 "(logxor 7) @result{} 7\n"
1353 "(logxor #b000 #b001 #b011) @result{} 2\n"
1354 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1356 #define FUNC_NAME s_scm_logxor
1360 if (SCM_UNBNDP (n2
))
1362 if (SCM_UNBNDP (n1
))
1364 else if (SCM_NUMBERP (n1
))
1367 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1370 if (SCM_I_INUMP (n1
))
1372 nn1
= SCM_I_INUM (n1
);
1373 if (SCM_I_INUMP (n2
))
1375 long nn2
= SCM_I_INUM (n2
);
1376 return SCM_I_MAKINUM (nn1
^ nn2
);
1378 else if (SCM_BIGP (n2
))
1382 SCM result_z
= scm_i_mkbig ();
1384 mpz_init_set_si (nn1_z
, nn1
);
1385 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1386 scm_remember_upto_here_1 (n2
);
1388 return scm_i_normbig (result_z
);
1392 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1394 else if (SCM_BIGP (n1
))
1396 if (SCM_I_INUMP (n2
))
1399 nn1
= SCM_I_INUM (n1
);
1402 else if (SCM_BIGP (n2
))
1404 SCM result_z
= scm_i_mkbig ();
1405 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1407 SCM_I_BIG_MPZ (n2
));
1408 scm_remember_upto_here_2 (n1
, n2
);
1409 return scm_i_normbig (result_z
);
1412 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1415 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1420 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1422 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1423 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1424 "without actually calculating the @code{logand}, just testing\n"
1428 "(logtest #b0100 #b1011) @result{} #f\n"
1429 "(logtest #b0100 #b0111) @result{} #t\n"
1431 #define FUNC_NAME s_scm_logtest
1435 if (SCM_I_INUMP (j
))
1437 nj
= SCM_I_INUM (j
);
1438 if (SCM_I_INUMP (k
))
1440 long nk
= SCM_I_INUM (k
);
1441 return scm_from_bool (nj
& nk
);
1443 else if (SCM_BIGP (k
))
1451 mpz_init_set_si (nj_z
, nj
);
1452 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1453 scm_remember_upto_here_1 (k
);
1454 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1460 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1462 else if (SCM_BIGP (j
))
1464 if (SCM_I_INUMP (k
))
1467 nj
= SCM_I_INUM (j
);
1470 else if (SCM_BIGP (k
))
1474 mpz_init (result_z
);
1478 scm_remember_upto_here_2 (j
, k
);
1479 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1480 mpz_clear (result_z
);
1484 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1487 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1492 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1494 "Test whether bit number @var{index} in @var{j} is set.\n"
1495 "@var{index} starts from 0 for the least significant bit.\n"
1498 "(logbit? 0 #b1101) @result{} #t\n"
1499 "(logbit? 1 #b1101) @result{} #f\n"
1500 "(logbit? 2 #b1101) @result{} #t\n"
1501 "(logbit? 3 #b1101) @result{} #t\n"
1502 "(logbit? 4 #b1101) @result{} #f\n"
1504 #define FUNC_NAME s_scm_logbit_p
1506 unsigned long int iindex
;
1507 iindex
= scm_to_ulong (index
);
1509 if (SCM_I_INUMP (j
))
1511 /* bits above what's in an inum follow the sign bit */
1512 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1513 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1515 else if (SCM_BIGP (j
))
1517 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1518 scm_remember_upto_here_1 (j
);
1519 return scm_from_bool (val
);
1522 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1527 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1529 "Return the integer which is the ones-complement of the integer\n"
1533 "(number->string (lognot #b10000000) 2)\n"
1534 " @result{} \"-10000001\"\n"
1535 "(number->string (lognot #b0) 2)\n"
1536 " @result{} \"-1\"\n"
1538 #define FUNC_NAME s_scm_lognot
1540 if (SCM_I_INUMP (n
)) {
1541 /* No overflow here, just need to toggle all the bits making up the inum.
1542 Enhancement: No need to strip the tag and add it back, could just xor
1543 a block of 1 bits, if that worked with the various debug versions of
1545 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1547 } else if (SCM_BIGP (n
)) {
1548 SCM result
= scm_i_mkbig ();
1549 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1550 scm_remember_upto_here_1 (n
);
1554 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1559 /* returns 0 if IN is not an integer. OUT must already be
1562 coerce_to_big (SCM in
, mpz_t out
)
1565 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1566 else if (SCM_I_INUMP (in
))
1567 mpz_set_si (out
, SCM_I_INUM (in
));
1574 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1575 (SCM n
, SCM k
, SCM m
),
1576 "Return @var{n} raised to the integer exponent\n"
1577 "@var{k}, modulo @var{m}.\n"
1580 "(modulo-expt 2 3 5)\n"
1583 #define FUNC_NAME s_scm_modulo_expt
1589 /* There are two classes of error we might encounter --
1590 1) Math errors, which we'll report by calling scm_num_overflow,
1592 2) wrong-type errors, which of course we'll report by calling
1594 We don't report those errors immediately, however; instead we do
1595 some cleanup first. These variables tell us which error (if
1596 any) we should report after cleaning up.
1598 int report_overflow
= 0;
1600 int position_of_wrong_type
= 0;
1601 SCM value_of_wrong_type
= SCM_INUM0
;
1603 SCM result
= SCM_UNDEFINED
;
1609 if (scm_is_eq (m
, SCM_INUM0
))
1611 report_overflow
= 1;
1615 if (!coerce_to_big (n
, n_tmp
))
1617 value_of_wrong_type
= n
;
1618 position_of_wrong_type
= 1;
1622 if (!coerce_to_big (k
, k_tmp
))
1624 value_of_wrong_type
= k
;
1625 position_of_wrong_type
= 2;
1629 if (!coerce_to_big (m
, m_tmp
))
1631 value_of_wrong_type
= m
;
1632 position_of_wrong_type
= 3;
1636 /* if the exponent K is negative, and we simply call mpz_powm, we
1637 will get a divide-by-zero exception when an inverse 1/n mod m
1638 doesn't exist (or is not unique). Since exceptions are hard to
1639 handle, we'll attempt the inversion "by hand" -- that way, we get
1640 a simple failure code, which is easy to handle. */
1642 if (-1 == mpz_sgn (k_tmp
))
1644 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1646 report_overflow
= 1;
1649 mpz_neg (k_tmp
, k_tmp
);
1652 result
= scm_i_mkbig ();
1653 mpz_powm (SCM_I_BIG_MPZ (result
),
1658 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1659 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1666 if (report_overflow
)
1667 scm_num_overflow (FUNC_NAME
);
1669 if (position_of_wrong_type
)
1670 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1671 value_of_wrong_type
);
1673 return scm_i_normbig (result
);
1677 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1679 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1680 "exact integer, @var{n} can be any number.\n"
1682 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1683 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1684 "includes @math{0^0} is 1.\n"
1687 "(integer-expt 2 5) @result{} 32\n"
1688 "(integer-expt -3 3) @result{} -27\n"
1689 "(integer-expt 5 -3) @result{} 1/125\n"
1690 "(integer-expt 0 0) @result{} 1\n"
1692 #define FUNC_NAME s_scm_integer_expt
1695 SCM z_i2
= SCM_BOOL_F
;
1697 SCM acc
= SCM_I_MAKINUM (1L);
1699 /* 0^0 == 1 according to R5RS */
1700 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1701 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1702 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1703 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1705 if (SCM_I_INUMP (k
))
1706 i2
= SCM_I_INUM (k
);
1707 else if (SCM_BIGP (k
))
1709 z_i2
= scm_i_clonebig (k
, 1);
1710 scm_remember_upto_here_1 (k
);
1714 SCM_WRONG_TYPE_ARG (2, k
);
1718 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1720 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1721 n
= scm_divide (n
, SCM_UNDEFINED
);
1725 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1729 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1731 return scm_product (acc
, n
);
1733 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1734 acc
= scm_product (acc
, n
);
1735 n
= scm_product (n
, n
);
1736 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1744 n
= scm_divide (n
, SCM_UNDEFINED
);
1751 return scm_product (acc
, n
);
1753 acc
= scm_product (acc
, n
);
1754 n
= scm_product (n
, n
);
1761 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1763 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1764 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1766 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1767 "@var{cnt} is negative it's a division, rounded towards negative\n"
1768 "infinity. (Note that this is not the same rounding as\n"
1769 "@code{quotient} does.)\n"
1771 "With @var{n} viewed as an infinite precision twos complement,\n"
1772 "@code{ash} means a left shift introducing zero bits, or a right\n"
1773 "shift dropping bits.\n"
1776 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1777 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1779 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1780 "(ash -23 -2) @result{} -6\n"
1782 #define FUNC_NAME s_scm_ash
1785 bits_to_shift
= scm_to_long (cnt
);
1787 if (SCM_I_INUMP (n
))
1789 long nn
= SCM_I_INUM (n
);
1791 if (bits_to_shift
> 0)
1793 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1794 overflow a non-zero fixnum. For smaller shifts we check the
1795 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1796 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1797 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1803 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1805 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1808 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1812 SCM result
= scm_i_long2big (nn
);
1813 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1820 bits_to_shift
= -bits_to_shift
;
1821 if (bits_to_shift
>= SCM_LONG_BIT
)
1822 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1824 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1828 else if (SCM_BIGP (n
))
1832 if (bits_to_shift
== 0)
1835 result
= scm_i_mkbig ();
1836 if (bits_to_shift
>= 0)
1838 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1844 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1845 we have to allocate a bignum even if the result is going to be a
1847 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1849 return scm_i_normbig (result
);
1855 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1861 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1862 (SCM n
, SCM start
, SCM end
),
1863 "Return the integer composed of the @var{start} (inclusive)\n"
1864 "through @var{end} (exclusive) bits of @var{n}. The\n"
1865 "@var{start}th bit becomes the 0-th bit in the result.\n"
1868 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1869 " @result{} \"1010\"\n"
1870 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1871 " @result{} \"10110\"\n"
1873 #define FUNC_NAME s_scm_bit_extract
1875 unsigned long int istart
, iend
, bits
;
1876 istart
= scm_to_ulong (start
);
1877 iend
= scm_to_ulong (end
);
1878 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1880 /* how many bits to keep */
1881 bits
= iend
- istart
;
1883 if (SCM_I_INUMP (n
))
1885 long int in
= SCM_I_INUM (n
);
1887 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1888 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1889 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1891 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1893 /* Since we emulate two's complement encoded numbers, this
1894 * special case requires us to produce a result that has
1895 * more bits than can be stored in a fixnum.
1897 SCM result
= scm_i_long2big (in
);
1898 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1903 /* mask down to requisite bits */
1904 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1905 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1907 else if (SCM_BIGP (n
))
1912 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1916 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1917 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1918 such bits into a ulong. */
1919 result
= scm_i_mkbig ();
1920 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1921 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1922 result
= scm_i_normbig (result
);
1924 scm_remember_upto_here_1 (n
);
1928 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1933 static const char scm_logtab
[] = {
1934 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1937 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1939 "Return the number of bits in integer @var{n}. If integer is\n"
1940 "positive, the 1-bits in its binary representation are counted.\n"
1941 "If negative, the 0-bits in its two's-complement binary\n"
1942 "representation are counted. If 0, 0 is returned.\n"
1945 "(logcount #b10101010)\n"
1952 #define FUNC_NAME s_scm_logcount
1954 if (SCM_I_INUMP (n
))
1956 unsigned long int c
= 0;
1957 long int nn
= SCM_I_INUM (n
);
1962 c
+= scm_logtab
[15 & nn
];
1965 return SCM_I_MAKINUM (c
);
1967 else if (SCM_BIGP (n
))
1969 unsigned long count
;
1970 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1971 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1973 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1974 scm_remember_upto_here_1 (n
);
1975 return SCM_I_MAKINUM (count
);
1978 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1983 static const char scm_ilentab
[] = {
1984 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1988 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1990 "Return the number of bits necessary to represent @var{n}.\n"
1993 "(integer-length #b10101010)\n"
1995 "(integer-length 0)\n"
1997 "(integer-length #b1111)\n"
2000 #define FUNC_NAME s_scm_integer_length
2002 if (SCM_I_INUMP (n
))
2004 unsigned long int c
= 0;
2006 long int nn
= SCM_I_INUM (n
);
2012 l
= scm_ilentab
[15 & nn
];
2015 return SCM_I_MAKINUM (c
- 4 + l
);
2017 else if (SCM_BIGP (n
))
2019 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2020 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2021 1 too big, so check for that and adjust. */
2022 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2023 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2024 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2025 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2027 scm_remember_upto_here_1 (n
);
2028 return SCM_I_MAKINUM (size
);
2031 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2035 /*** NUMBERS -> STRINGS ***/
2036 #define SCM_MAX_DBL_PREC 60
2037 #define SCM_MAX_DBL_RADIX 36
2039 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2040 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2041 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2044 void init_dblprec(int *prec
, int radix
) {
2045 /* determine floating point precision by adding successively
2046 smaller increments to 1.0 until it is considered == 1.0 */
2047 double f
= ((double)1.0)/radix
;
2048 double fsum
= 1.0 + f
;
2053 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2065 void init_fx_radix(double *fx_list
, int radix
)
2067 /* initialize a per-radix list of tolerances. When added
2068 to a number < 1.0, we can determine if we should raund
2069 up and quit converting a number to a string. */
2073 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2074 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2077 /* use this array as a way to generate a single digit */
2078 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2081 idbl2str (double f
, char *a
, int radix
)
2083 int efmt
, dpt
, d
, i
, wp
;
2085 #ifdef DBL_MIN_10_EXP
2088 #endif /* DBL_MIN_10_EXP */
2093 radix
> SCM_MAX_DBL_RADIX
)
2095 /* revert to existing behavior */
2099 wp
= scm_dblprec
[radix
-2];
2100 fx
= fx_per_radix
[radix
-2];
2104 #ifdef HAVE_COPYSIGN
2105 double sgn
= copysign (1.0, f
);
2110 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2116 strcpy (a
, "-inf.0");
2118 strcpy (a
, "+inf.0");
2121 else if (xisnan (f
))
2123 strcpy (a
, "+nan.0");
2133 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2134 make-uniform-vector, from causing infinite loops. */
2135 /* just do the checking...if it passes, we do the conversion for our
2136 radix again below */
2143 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2151 while (f_cpy
> 10.0)
2154 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2175 if (f
+ fx
[wp
] >= radix
)
2182 /* adding 9999 makes this equivalent to abs(x) % 3 */
2183 dpt
= (exp
+ 9999) % 3;
2187 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2209 a
[ch
++] = number_chars
[d
];
2212 if (f
+ fx
[wp
] >= 1.0)
2214 a
[ch
- 1] = number_chars
[d
+1];
2226 if ((dpt
> 4) && (exp
> 6))
2228 d
= (a
[0] == '-' ? 2 : 1);
2229 for (i
= ch
++; i
> d
; i
--)
2242 if (a
[ch
- 1] == '.')
2243 a
[ch
++] = '0'; /* trailing zero */
2252 for (i
= radix
; i
<= exp
; i
*= radix
);
2253 for (i
/= radix
; i
; i
/= radix
)
2255 a
[ch
++] = number_chars
[exp
/ i
];
2264 icmplx2str (double real
, double imag
, char *str
, int radix
)
2268 i
= idbl2str (real
, str
, radix
);
2271 /* Don't output a '+' for negative numbers or for Inf and
2272 NaN. They will provide their own sign. */
2273 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2275 i
+= idbl2str (imag
, &str
[i
], radix
);
2282 iflo2str (SCM flt
, char *str
, int radix
)
2285 if (SCM_REALP (flt
))
2286 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2288 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2293 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2294 characters in the result.
2296 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2298 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2303 return scm_iuint2str (-num
, rad
, p
) + 1;
2306 return scm_iuint2str (num
, rad
, p
);
2309 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2310 characters in the result.
2312 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2314 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2318 scm_t_uintmax n
= num
;
2320 for (n
/= rad
; n
> 0; n
/= rad
)
2330 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2335 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2337 "Return a string holding the external representation of the\n"
2338 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2339 "inexact, a radix of 10 will be used.")
2340 #define FUNC_NAME s_scm_number_to_string
2344 if (SCM_UNBNDP (radix
))
2347 base
= scm_to_signed_integer (radix
, 2, 36);
2349 if (SCM_I_INUMP (n
))
2351 char num_buf
[SCM_INTBUFLEN
];
2352 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2353 return scm_from_locale_stringn (num_buf
, length
);
2355 else if (SCM_BIGP (n
))
2357 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2358 scm_remember_upto_here_1 (n
);
2359 return scm_take_locale_string (str
);
2361 else if (SCM_FRACTIONP (n
))
2363 scm_i_fraction_reduce (n
);
2364 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2365 scm_from_locale_string ("/"),
2366 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2368 else if (SCM_INEXACTP (n
))
2370 char num_buf
[FLOBUFLEN
];
2371 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2374 SCM_WRONG_TYPE_ARG (1, n
);
2379 /* These print routines used to be stubbed here so that scm_repl.c
2380 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2383 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2385 char num_buf
[FLOBUFLEN
];
2386 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2391 scm_i_print_double (double val
, SCM port
)
2393 char num_buf
[FLOBUFLEN
];
2394 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2398 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2401 char num_buf
[FLOBUFLEN
];
2402 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2407 scm_i_print_complex (double real
, double imag
, SCM port
)
2409 char num_buf
[FLOBUFLEN
];
2410 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2414 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2417 scm_i_fraction_reduce (sexp
);
2418 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2419 scm_lfwrite (scm_i_string_chars (str
), scm_i_string_length (str
), port
);
2420 scm_remember_upto_here_1 (str
);
2425 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2427 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2428 scm_remember_upto_here_1 (exp
);
2429 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2433 /*** END nums->strs ***/
2436 /*** STRINGS -> NUMBERS ***/
2438 /* The following functions implement the conversion from strings to numbers.
2439 * The implementation somehow follows the grammar for numbers as it is given
2440 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2441 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2442 * points should be noted about the implementation:
2443 * * Each function keeps a local index variable 'idx' that points at the
2444 * current position within the parsed string. The global index is only
2445 * updated if the function could parse the corresponding syntactic unit
2447 * * Similarly, the functions keep track of indicators of inexactness ('#',
2448 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2449 * global exactness information is only updated after each part has been
2450 * successfully parsed.
2451 * * Sequences of digits are parsed into temporary variables holding fixnums.
2452 * Only if these fixnums would overflow, the result variables are updated
2453 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2454 * the temporary variables holding the fixnums are cleared, and the process
2455 * starts over again. If for example fixnums were able to store five decimal
2456 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2457 * and the result was computed as 12345 * 100000 + 67890. In other words,
2458 * only every five digits two bignum operations were performed.
2461 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2463 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2465 /* In non ASCII-style encodings the following macro might not work. */
2466 #define XDIGIT2UINT(d) \
2467 (isdigit ((int) (unsigned char) d) \
2469 : tolower ((int) (unsigned char) d) - 'a' + 10)
2472 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2473 unsigned int radix
, enum t_exactness
*p_exactness
)
2475 unsigned int idx
= *p_idx
;
2476 unsigned int hash_seen
= 0;
2477 scm_t_bits shift
= 1;
2479 unsigned int digit_value
;
2487 if (!isxdigit ((int) (unsigned char) c
))
2489 digit_value
= XDIGIT2UINT (c
);
2490 if (digit_value
>= radix
)
2494 result
= SCM_I_MAKINUM (digit_value
);
2498 if (isxdigit ((int) (unsigned char) c
))
2502 digit_value
= XDIGIT2UINT (c
);
2503 if (digit_value
>= radix
)
2515 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2517 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2519 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2526 shift
= shift
* radix
;
2527 add
= add
* radix
+ digit_value
;
2532 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2534 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2538 *p_exactness
= INEXACT
;
2544 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2545 * covers the parts of the rules that start at a potential point. The value
2546 * of the digits up to the point have been parsed by the caller and are given
2547 * in variable result. The content of *p_exactness indicates, whether a hash
2548 * has already been seen in the digits before the point.
2551 /* In non ASCII-style encodings the following macro might not work. */
2552 #define DIGIT2UINT(d) ((d) - '0')
2555 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2556 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2558 unsigned int idx
= *p_idx
;
2559 enum t_exactness x
= *p_exactness
;
2564 if (mem
[idx
] == '.')
2566 scm_t_bits shift
= 1;
2568 unsigned int digit_value
;
2569 SCM big_shift
= SCM_I_MAKINUM (1);
2575 if (isdigit ((int) (unsigned char) c
))
2580 digit_value
= DIGIT2UINT (c
);
2591 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2593 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2594 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2596 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2604 add
= add
* 10 + digit_value
;
2610 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2611 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2612 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2615 result
= scm_divide (result
, big_shift
);
2617 /* We've seen a decimal point, thus the value is implicitly inexact. */
2629 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2656 if (!isdigit ((int) (unsigned char) c
))
2660 exponent
= DIGIT2UINT (c
);
2664 if (isdigit ((int) (unsigned char) c
))
2667 if (exponent
<= SCM_MAXEXP
)
2668 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2674 if (exponent
> SCM_MAXEXP
)
2676 size_t exp_len
= idx
- start
;
2677 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2678 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2679 scm_out_of_range ("string->number", exp_num
);
2682 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2684 result
= scm_product (result
, e
);
2686 result
= scm_divide2real (result
, e
);
2688 /* We've seen an exponent, thus the value is implicitly inexact. */
2706 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2709 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2710 unsigned int radix
, enum t_exactness
*p_exactness
)
2712 unsigned int idx
= *p_idx
;
2718 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2724 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2726 enum t_exactness x
= EXACT
;
2728 /* Cobble up the fractional part. We might want to set the
2729 NaN's mantissa from it. */
2731 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2736 if (mem
[idx
] == '.')
2740 else if (idx
+ 1 == len
)
2742 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2745 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2746 p_idx
, p_exactness
);
2750 enum t_exactness x
= EXACT
;
2753 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2754 if (scm_is_false (uinteger
))
2759 else if (mem
[idx
] == '/')
2765 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2766 if (scm_is_false (divisor
))
2769 /* both are int/big here, I assume */
2770 result
= scm_i_make_ratio (uinteger
, divisor
);
2772 else if (radix
== 10)
2774 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2775 if (scm_is_false (result
))
2786 /* When returning an inexact zero, make sure it is represented as a
2787 floating point value so that we can change its sign.
2789 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2790 result
= scm_from_double (0.0);
2796 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2799 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2800 unsigned int radix
, enum t_exactness
*p_exactness
)
2824 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2825 if (scm_is_false (ureal
))
2827 /* input must be either +i or -i */
2832 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2838 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2845 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2846 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2855 /* either +<ureal>i or -<ureal>i */
2862 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2865 /* polar input: <real>@<real>. */
2890 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2891 if (scm_is_false (angle
))
2896 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2897 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2899 result
= scm_make_polar (ureal
, angle
);
2904 /* expecting input matching <real>[+-]<ureal>?i */
2911 int sign
= (c
== '+') ? 1 : -1;
2912 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2914 if (scm_is_false (imag
))
2915 imag
= SCM_I_MAKINUM (sign
);
2916 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2917 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2921 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2928 return scm_make_rectangular (ureal
, imag
);
2937 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2939 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2942 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
2943 unsigned int default_radix
)
2945 unsigned int idx
= 0;
2946 unsigned int radix
= NO_RADIX
;
2947 enum t_exactness forced_x
= NO_EXACTNESS
;
2948 enum t_exactness implicit_x
= EXACT
;
2951 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2952 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2954 switch (mem
[idx
+ 1])
2957 if (radix
!= NO_RADIX
)
2962 if (radix
!= NO_RADIX
)
2967 if (forced_x
!= NO_EXACTNESS
)
2972 if (forced_x
!= NO_EXACTNESS
)
2977 if (radix
!= NO_RADIX
)
2982 if (radix
!= NO_RADIX
)
2992 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2993 if (radix
== NO_RADIX
)
2994 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2996 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2998 if (scm_is_false (result
))
3004 if (SCM_INEXACTP (result
))
3005 return scm_inexact_to_exact (result
);
3009 if (SCM_INEXACTP (result
))
3012 return scm_exact_to_inexact (result
);
3015 if (implicit_x
== INEXACT
)
3017 if (SCM_INEXACTP (result
))
3020 return scm_exact_to_inexact (result
);
3028 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3029 (SCM string
, SCM radix
),
3030 "Return a number of the maximally precise representation\n"
3031 "expressed by the given @var{string}. @var{radix} must be an\n"
3032 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3033 "is a default radix that may be overridden by an explicit radix\n"
3034 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3035 "supplied, then the default radix is 10. If string is not a\n"
3036 "syntactically valid notation for a number, then\n"
3037 "@code{string->number} returns @code{#f}.")
3038 #define FUNC_NAME s_scm_string_to_number
3042 SCM_VALIDATE_STRING (1, string
);
3044 if (SCM_UNBNDP (radix
))
3047 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3049 answer
= scm_c_locale_stringn_to_number (scm_i_string_chars (string
),
3050 scm_i_string_length (string
),
3052 scm_remember_upto_here_1 (string
);
3058 /*** END strs->nums ***/
3062 scm_bigequal (SCM x
, SCM y
)
3064 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3065 scm_remember_upto_here_2 (x
, y
);
3066 return scm_from_bool (0 == result
);
3070 scm_real_equalp (SCM x
, SCM y
)
3072 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3076 scm_complex_equalp (SCM x
, SCM y
)
3078 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3079 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3083 scm_i_fraction_equalp (SCM x
, SCM y
)
3085 scm_i_fraction_reduce (x
);
3086 scm_i_fraction_reduce (y
);
3087 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3088 SCM_FRACTION_NUMERATOR (y
)))
3089 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3090 SCM_FRACTION_DENOMINATOR (y
))))
3097 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3099 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3101 #define FUNC_NAME s_scm_number_p
3103 return scm_from_bool (SCM_NUMBERP (x
));
3107 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3109 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3110 "otherwise. Note that the sets of real, rational and integer\n"
3111 "values form subsets of the set of complex numbers, i. e. the\n"
3112 "predicate will also be fulfilled if @var{x} is a real,\n"
3113 "rational or integer number.")
3114 #define FUNC_NAME s_scm_complex_p
3116 /* all numbers are complex. */
3117 return scm_number_p (x
);
3121 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3123 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3124 "otherwise. Note that the set of integer values forms a subset of\n"
3125 "the set of real numbers, i. e. the predicate will also be\n"
3126 "fulfilled if @var{x} is an integer number.")
3127 #define FUNC_NAME s_scm_real_p
3129 /* we can't represent irrational numbers. */
3130 return scm_rational_p (x
);
3134 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3136 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3137 "otherwise. Note that the set of integer values forms a subset of\n"
3138 "the set of rational numbers, i. e. the predicate will also be\n"
3139 "fulfilled if @var{x} is an integer number.")
3140 #define FUNC_NAME s_scm_rational_p
3142 if (SCM_I_INUMP (x
))
3144 else if (SCM_IMP (x
))
3146 else if (SCM_BIGP (x
))
3148 else if (SCM_FRACTIONP (x
))
3150 else if (SCM_REALP (x
))
3151 /* due to their limited precision, all floating point numbers are
3152 rational as well. */
3159 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3161 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3163 #define FUNC_NAME s_scm_integer_p
3166 if (SCM_I_INUMP (x
))
3172 if (!SCM_INEXACTP (x
))
3174 if (SCM_COMPLEXP (x
))
3176 r
= SCM_REAL_VALUE (x
);
3177 /* +/-inf passes r==floor(r), making those #t */
3185 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3187 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3189 #define FUNC_NAME s_scm_inexact_p
3191 if (SCM_INEXACTP (x
))
3193 if (SCM_NUMBERP (x
))
3195 SCM_WRONG_TYPE_ARG (1, x
);
3200 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3201 /* "Return @code{#t} if all parameters are numerically equal." */
3203 scm_num_eq_p (SCM x
, SCM y
)
3206 if (SCM_I_INUMP (x
))
3208 long xx
= SCM_I_INUM (x
);
3209 if (SCM_I_INUMP (y
))
3211 long yy
= SCM_I_INUM (y
);
3212 return scm_from_bool (xx
== yy
);
3214 else if (SCM_BIGP (y
))
3216 else if (SCM_REALP (y
))
3218 /* On a 32-bit system an inum fits a double, we can cast the inum
3219 to a double and compare.
3221 But on a 64-bit system an inum is bigger than a double and
3222 casting it to a double (call that dxx) will round. dxx is at
3223 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3224 an integer and fits a long. So we cast yy to a long and
3225 compare with plain xx.
3227 An alternative (for any size system actually) would be to check
3228 yy is an integer (with floor) and is in range of an inum
3229 (compare against appropriate powers of 2) then test
3230 xx==(long)yy. It's just a matter of which casts/comparisons
3231 might be fastest or easiest for the cpu. */
3233 double yy
= SCM_REAL_VALUE (y
);
3234 return scm_from_bool ((double) xx
== yy
3235 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3236 || xx
== (long) yy
));
3238 else if (SCM_COMPLEXP (y
))
3239 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3240 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3241 else if (SCM_FRACTIONP (y
))
3244 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3246 else if (SCM_BIGP (x
))
3248 if (SCM_I_INUMP (y
))
3250 else if (SCM_BIGP (y
))
3252 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3253 scm_remember_upto_here_2 (x
, y
);
3254 return scm_from_bool (0 == cmp
);
3256 else if (SCM_REALP (y
))
3259 if (xisnan (SCM_REAL_VALUE (y
)))
3261 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3262 scm_remember_upto_here_1 (x
);
3263 return scm_from_bool (0 == cmp
);
3265 else if (SCM_COMPLEXP (y
))
3268 if (0.0 != SCM_COMPLEX_IMAG (y
))
3270 if (xisnan (SCM_COMPLEX_REAL (y
)))
3272 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3273 scm_remember_upto_here_1 (x
);
3274 return scm_from_bool (0 == cmp
);
3276 else if (SCM_FRACTIONP (y
))
3279 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3281 else if (SCM_REALP (x
))
3283 double xx
= SCM_REAL_VALUE (x
);
3284 if (SCM_I_INUMP (y
))
3286 /* see comments with inum/real above */
3287 long yy
= SCM_I_INUM (y
);
3288 return scm_from_bool (xx
== (double) yy
3289 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3290 || (long) xx
== yy
));
3292 else if (SCM_BIGP (y
))
3295 if (xisnan (SCM_REAL_VALUE (x
)))
3297 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3298 scm_remember_upto_here_1 (y
);
3299 return scm_from_bool (0 == cmp
);
3301 else if (SCM_REALP (y
))
3302 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3303 else if (SCM_COMPLEXP (y
))
3304 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3305 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3306 else if (SCM_FRACTIONP (y
))
3308 double xx
= SCM_REAL_VALUE (x
);
3312 return scm_from_bool (xx
< 0.0);
3313 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3317 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3319 else if (SCM_COMPLEXP (x
))
3321 if (SCM_I_INUMP (y
))
3322 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3323 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3324 else if (SCM_BIGP (y
))
3327 if (0.0 != SCM_COMPLEX_IMAG (x
))
3329 if (xisnan (SCM_COMPLEX_REAL (x
)))
3331 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3332 scm_remember_upto_here_1 (y
);
3333 return scm_from_bool (0 == cmp
);
3335 else if (SCM_REALP (y
))
3336 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3337 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3338 else if (SCM_COMPLEXP (y
))
3339 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3340 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3341 else if (SCM_FRACTIONP (y
))
3344 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3346 xx
= SCM_COMPLEX_REAL (x
);
3350 return scm_from_bool (xx
< 0.0);
3351 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3355 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3357 else if (SCM_FRACTIONP (x
))
3359 if (SCM_I_INUMP (y
))
3361 else if (SCM_BIGP (y
))
3363 else if (SCM_REALP (y
))
3365 double yy
= SCM_REAL_VALUE (y
);
3369 return scm_from_bool (0.0 < yy
);
3370 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3373 else if (SCM_COMPLEXP (y
))
3376 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3378 yy
= SCM_COMPLEX_REAL (y
);
3382 return scm_from_bool (0.0 < yy
);
3383 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3386 else if (SCM_FRACTIONP (y
))
3387 return scm_i_fraction_equalp (x
, y
);
3389 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3392 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3396 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3397 done are good for inums, but for bignums an answer can almost always be
3398 had by just examining a few high bits of the operands, as done by GMP in
3399 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3400 of the float exponent to take into account. */
3402 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3403 /* "Return @code{#t} if the list of parameters is monotonically\n"
3407 scm_less_p (SCM x
, SCM y
)
3410 if (SCM_I_INUMP (x
))
3412 long xx
= SCM_I_INUM (x
);
3413 if (SCM_I_INUMP (y
))
3415 long yy
= SCM_I_INUM (y
);
3416 return scm_from_bool (xx
< yy
);
3418 else if (SCM_BIGP (y
))
3420 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3421 scm_remember_upto_here_1 (y
);
3422 return scm_from_bool (sgn
> 0);
3424 else if (SCM_REALP (y
))
3425 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3426 else if (SCM_FRACTIONP (y
))
3428 /* "x < a/b" becomes "x*b < a" */
3430 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3431 y
= SCM_FRACTION_NUMERATOR (y
);
3435 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3437 else if (SCM_BIGP (x
))
3439 if (SCM_I_INUMP (y
))
3441 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3442 scm_remember_upto_here_1 (x
);
3443 return scm_from_bool (sgn
< 0);
3445 else if (SCM_BIGP (y
))
3447 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3448 scm_remember_upto_here_2 (x
, y
);
3449 return scm_from_bool (cmp
< 0);
3451 else if (SCM_REALP (y
))
3454 if (xisnan (SCM_REAL_VALUE (y
)))
3456 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3457 scm_remember_upto_here_1 (x
);
3458 return scm_from_bool (cmp
< 0);
3460 else if (SCM_FRACTIONP (y
))
3463 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3465 else if (SCM_REALP (x
))
3467 if (SCM_I_INUMP (y
))
3468 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3469 else if (SCM_BIGP (y
))
3472 if (xisnan (SCM_REAL_VALUE (x
)))
3474 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3475 scm_remember_upto_here_1 (y
);
3476 return scm_from_bool (cmp
> 0);
3478 else if (SCM_REALP (y
))
3479 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3480 else if (SCM_FRACTIONP (y
))
3482 double xx
= SCM_REAL_VALUE (x
);
3486 return scm_from_bool (xx
< 0.0);
3487 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3491 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3493 else if (SCM_FRACTIONP (x
))
3495 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3497 /* "a/b < y" becomes "a < y*b" */
3498 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3499 x
= SCM_FRACTION_NUMERATOR (x
);
3502 else if (SCM_REALP (y
))
3504 double yy
= SCM_REAL_VALUE (y
);
3508 return scm_from_bool (0.0 < yy
);
3509 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3512 else if (SCM_FRACTIONP (y
))
3514 /* "a/b < c/d" becomes "a*d < c*b" */
3515 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3516 SCM_FRACTION_DENOMINATOR (y
));
3517 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3518 SCM_FRACTION_DENOMINATOR (x
));
3524 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3527 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3531 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3532 /* "Return @code{#t} if the list of parameters is monotonically\n"
3535 #define FUNC_NAME s_scm_gr_p
3537 scm_gr_p (SCM x
, SCM y
)
3539 if (!SCM_NUMBERP (x
))
3540 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3541 else if (!SCM_NUMBERP (y
))
3542 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3544 return scm_less_p (y
, x
);
3549 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3550 /* "Return @code{#t} if the list of parameters is monotonically\n"
3553 #define FUNC_NAME s_scm_leq_p
3555 scm_leq_p (SCM x
, SCM y
)
3557 if (!SCM_NUMBERP (x
))
3558 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3559 else if (!SCM_NUMBERP (y
))
3560 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3561 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3564 return scm_not (scm_less_p (y
, x
));
3569 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3570 /* "Return @code{#t} if the list of parameters is monotonically\n"
3573 #define FUNC_NAME s_scm_geq_p
3575 scm_geq_p (SCM x
, SCM y
)
3577 if (!SCM_NUMBERP (x
))
3578 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3579 else if (!SCM_NUMBERP (y
))
3580 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3581 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3584 return scm_not (scm_less_p (x
, y
));
3589 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3590 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3596 if (SCM_I_INUMP (z
))
3597 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3598 else if (SCM_BIGP (z
))
3600 else if (SCM_REALP (z
))
3601 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3602 else if (SCM_COMPLEXP (z
))
3603 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3604 && SCM_COMPLEX_IMAG (z
) == 0.0);
3605 else if (SCM_FRACTIONP (z
))
3608 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3612 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3613 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3617 scm_positive_p (SCM x
)
3619 if (SCM_I_INUMP (x
))
3620 return scm_from_bool (SCM_I_INUM (x
) > 0);
3621 else if (SCM_BIGP (x
))
3623 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3624 scm_remember_upto_here_1 (x
);
3625 return scm_from_bool (sgn
> 0);
3627 else if (SCM_REALP (x
))
3628 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3629 else if (SCM_FRACTIONP (x
))
3630 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3632 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3636 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3637 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3641 scm_negative_p (SCM x
)
3643 if (SCM_I_INUMP (x
))
3644 return scm_from_bool (SCM_I_INUM (x
) < 0);
3645 else if (SCM_BIGP (x
))
3647 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3648 scm_remember_upto_here_1 (x
);
3649 return scm_from_bool (sgn
< 0);
3651 else if (SCM_REALP (x
))
3652 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3653 else if (SCM_FRACTIONP (x
))
3654 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3656 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3660 /* scm_min and scm_max return an inexact when either argument is inexact, as
3661 required by r5rs. On that basis, for exact/inexact combinations the
3662 exact is converted to inexact to compare and possibly return. This is
3663 unlike scm_less_p above which takes some trouble to preserve all bits in
3664 its test, such trouble is not required for min and max. */
3666 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3667 /* "Return the maximum of all parameter values."
3670 scm_max (SCM x
, SCM y
)
3675 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3676 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3679 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3682 if (SCM_I_INUMP (x
))
3684 long xx
= SCM_I_INUM (x
);
3685 if (SCM_I_INUMP (y
))
3687 long yy
= SCM_I_INUM (y
);
3688 return (xx
< yy
) ? y
: x
;
3690 else if (SCM_BIGP (y
))
3692 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3693 scm_remember_upto_here_1 (y
);
3694 return (sgn
< 0) ? x
: y
;
3696 else if (SCM_REALP (y
))
3699 /* if y==NaN then ">" is false and we return NaN */
3700 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3702 else if (SCM_FRACTIONP (y
))
3705 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3708 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3710 else if (SCM_BIGP (x
))
3712 if (SCM_I_INUMP (y
))
3714 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3715 scm_remember_upto_here_1 (x
);
3716 return (sgn
< 0) ? y
: x
;
3718 else if (SCM_BIGP (y
))
3720 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3721 scm_remember_upto_here_2 (x
, y
);
3722 return (cmp
> 0) ? x
: y
;
3724 else if (SCM_REALP (y
))
3726 /* if y==NaN then xx>yy is false, so we return the NaN y */
3729 xx
= scm_i_big2dbl (x
);
3730 yy
= SCM_REAL_VALUE (y
);
3731 return (xx
> yy
? scm_from_double (xx
) : y
);
3733 else if (SCM_FRACTIONP (y
))
3738 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3740 else if (SCM_REALP (x
))
3742 if (SCM_I_INUMP (y
))
3744 double z
= SCM_I_INUM (y
);
3745 /* if x==NaN then "<" is false and we return NaN */
3746 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3748 else if (SCM_BIGP (y
))
3753 else if (SCM_REALP (y
))
3755 /* if x==NaN then our explicit check means we return NaN
3756 if y==NaN then ">" is false and we return NaN
3757 calling isnan is unavoidable, since it's the only way to know
3758 which of x or y causes any compares to be false */
3759 double xx
= SCM_REAL_VALUE (x
);
3760 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3762 else if (SCM_FRACTIONP (y
))
3764 double yy
= scm_i_fraction2double (y
);
3765 double xx
= SCM_REAL_VALUE (x
);
3766 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3769 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3771 else if (SCM_FRACTIONP (x
))
3773 if (SCM_I_INUMP (y
))
3777 else if (SCM_BIGP (y
))
3781 else if (SCM_REALP (y
))
3783 double xx
= scm_i_fraction2double (x
);
3784 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3786 else if (SCM_FRACTIONP (y
))
3791 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3794 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3798 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3799 /* "Return the minium of all parameter values."
3802 scm_min (SCM x
, SCM y
)
3807 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3808 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3811 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3814 if (SCM_I_INUMP (x
))
3816 long xx
= SCM_I_INUM (x
);
3817 if (SCM_I_INUMP (y
))
3819 long yy
= SCM_I_INUM (y
);
3820 return (xx
< yy
) ? x
: y
;
3822 else if (SCM_BIGP (y
))
3824 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3825 scm_remember_upto_here_1 (y
);
3826 return (sgn
< 0) ? y
: x
;
3828 else if (SCM_REALP (y
))
3831 /* if y==NaN then "<" is false and we return NaN */
3832 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3834 else if (SCM_FRACTIONP (y
))
3837 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3840 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3842 else if (SCM_BIGP (x
))
3844 if (SCM_I_INUMP (y
))
3846 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3847 scm_remember_upto_here_1 (x
);
3848 return (sgn
< 0) ? x
: y
;
3850 else if (SCM_BIGP (y
))
3852 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3853 scm_remember_upto_here_2 (x
, y
);
3854 return (cmp
> 0) ? y
: x
;
3856 else if (SCM_REALP (y
))
3858 /* if y==NaN then xx<yy is false, so we return the NaN y */
3861 xx
= scm_i_big2dbl (x
);
3862 yy
= SCM_REAL_VALUE (y
);
3863 return (xx
< yy
? scm_from_double (xx
) : y
);
3865 else if (SCM_FRACTIONP (y
))
3870 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3872 else if (SCM_REALP (x
))
3874 if (SCM_I_INUMP (y
))
3876 double z
= SCM_I_INUM (y
);
3877 /* if x==NaN then "<" is false and we return NaN */
3878 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3880 else if (SCM_BIGP (y
))
3885 else if (SCM_REALP (y
))
3887 /* if x==NaN then our explicit check means we return NaN
3888 if y==NaN then "<" is false and we return NaN
3889 calling isnan is unavoidable, since it's the only way to know
3890 which of x or y causes any compares to be false */
3891 double xx
= SCM_REAL_VALUE (x
);
3892 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3894 else if (SCM_FRACTIONP (y
))
3896 double yy
= scm_i_fraction2double (y
);
3897 double xx
= SCM_REAL_VALUE (x
);
3898 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3901 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3903 else if (SCM_FRACTIONP (x
))
3905 if (SCM_I_INUMP (y
))
3909 else if (SCM_BIGP (y
))
3913 else if (SCM_REALP (y
))
3915 double xx
= scm_i_fraction2double (x
);
3916 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3918 else if (SCM_FRACTIONP (y
))
3923 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3926 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3930 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3931 /* "Return the sum of all parameter values. Return 0 if called without\n"
3935 scm_sum (SCM x
, SCM y
)
3939 if (SCM_NUMBERP (x
)) return x
;
3940 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3941 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3944 if (SCM_I_INUMP (x
))
3946 if (SCM_I_INUMP (y
))
3948 long xx
= SCM_I_INUM (x
);
3949 long yy
= SCM_I_INUM (y
);
3950 long int z
= xx
+ yy
;
3951 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3953 else if (SCM_BIGP (y
))
3958 else if (SCM_REALP (y
))
3960 long int xx
= SCM_I_INUM (x
);
3961 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3963 else if (SCM_COMPLEXP (y
))
3965 long int xx
= SCM_I_INUM (x
);
3966 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3967 SCM_COMPLEX_IMAG (y
));
3969 else if (SCM_FRACTIONP (y
))
3970 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3971 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3972 SCM_FRACTION_DENOMINATOR (y
));
3974 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3975 } else if (SCM_BIGP (x
))
3977 if (SCM_I_INUMP (y
))
3982 inum
= SCM_I_INUM (y
);
3985 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3988 SCM result
= scm_i_mkbig ();
3989 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3990 scm_remember_upto_here_1 (x
);
3991 /* we know the result will have to be a bignum */
3994 return scm_i_normbig (result
);
3998 SCM result
= scm_i_mkbig ();
3999 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4000 scm_remember_upto_here_1 (x
);
4001 /* we know the result will have to be a bignum */
4004 return scm_i_normbig (result
);
4007 else if (SCM_BIGP (y
))
4009 SCM result
= scm_i_mkbig ();
4010 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4011 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4012 mpz_add (SCM_I_BIG_MPZ (result
),
4015 scm_remember_upto_here_2 (x
, y
);
4016 /* we know the result will have to be a bignum */
4019 return scm_i_normbig (result
);
4021 else if (SCM_REALP (y
))
4023 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4024 scm_remember_upto_here_1 (x
);
4025 return scm_from_double (result
);
4027 else if (SCM_COMPLEXP (y
))
4029 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4030 + SCM_COMPLEX_REAL (y
));
4031 scm_remember_upto_here_1 (x
);
4032 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4034 else if (SCM_FRACTIONP (y
))
4035 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4036 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4037 SCM_FRACTION_DENOMINATOR (y
));
4039 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4041 else if (SCM_REALP (x
))
4043 if (SCM_I_INUMP (y
))
4044 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4045 else if (SCM_BIGP (y
))
4047 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4048 scm_remember_upto_here_1 (y
);
4049 return scm_from_double (result
);
4051 else if (SCM_REALP (y
))
4052 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4053 else if (SCM_COMPLEXP (y
))
4054 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4055 SCM_COMPLEX_IMAG (y
));
4056 else if (SCM_FRACTIONP (y
))
4057 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4059 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4061 else if (SCM_COMPLEXP (x
))
4063 if (SCM_I_INUMP (y
))
4064 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4065 SCM_COMPLEX_IMAG (x
));
4066 else if (SCM_BIGP (y
))
4068 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4069 + SCM_COMPLEX_REAL (x
));
4070 scm_remember_upto_here_1 (y
);
4071 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4073 else if (SCM_REALP (y
))
4074 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4075 SCM_COMPLEX_IMAG (x
));
4076 else if (SCM_COMPLEXP (y
))
4077 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4078 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4079 else if (SCM_FRACTIONP (y
))
4080 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4081 SCM_COMPLEX_IMAG (x
));
4083 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4085 else if (SCM_FRACTIONP (x
))
4087 if (SCM_I_INUMP (y
))
4088 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4089 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4090 SCM_FRACTION_DENOMINATOR (x
));
4091 else if (SCM_BIGP (y
))
4092 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4093 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4094 SCM_FRACTION_DENOMINATOR (x
));
4095 else if (SCM_REALP (y
))
4096 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4097 else if (SCM_COMPLEXP (y
))
4098 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4099 SCM_COMPLEX_IMAG (y
));
4100 else if (SCM_FRACTIONP (y
))
4101 /* a/b + c/d = (ad + bc) / bd */
4102 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4103 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4104 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4106 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4109 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4113 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4115 "Return @math{@var{x}+1}.")
4116 #define FUNC_NAME s_scm_oneplus
4118 return scm_sum (x
, SCM_I_MAKINUM (1));
4123 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4124 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4125 * the sum of all but the first argument are subtracted from the first
4127 #define FUNC_NAME s_difference
4129 scm_difference (SCM x
, SCM y
)
4134 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4136 if (SCM_I_INUMP (x
))
4138 long xx
= -SCM_I_INUM (x
);
4139 if (SCM_FIXABLE (xx
))
4140 return SCM_I_MAKINUM (xx
);
4142 return scm_i_long2big (xx
);
4144 else if (SCM_BIGP (x
))
4145 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4146 bignum, but negating that gives a fixnum. */
4147 return scm_i_normbig (scm_i_clonebig (x
, 0));
4148 else if (SCM_REALP (x
))
4149 return scm_from_double (-SCM_REAL_VALUE (x
));
4150 else if (SCM_COMPLEXP (x
))
4151 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4152 -SCM_COMPLEX_IMAG (x
));
4153 else if (SCM_FRACTIONP (x
))
4154 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4155 SCM_FRACTION_DENOMINATOR (x
));
4157 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4160 if (SCM_I_INUMP (x
))
4162 if (SCM_I_INUMP (y
))
4164 long int xx
= SCM_I_INUM (x
);
4165 long int yy
= SCM_I_INUM (y
);
4166 long int z
= xx
- yy
;
4167 if (SCM_FIXABLE (z
))
4168 return SCM_I_MAKINUM (z
);
4170 return scm_i_long2big (z
);
4172 else if (SCM_BIGP (y
))
4174 /* inum-x - big-y */
4175 long xx
= SCM_I_INUM (x
);
4178 return scm_i_clonebig (y
, 0);
4181 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4182 SCM result
= scm_i_mkbig ();
4185 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4188 /* x - y == -(y + -x) */
4189 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4190 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4192 scm_remember_upto_here_1 (y
);
4194 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4195 /* we know the result will have to be a bignum */
4198 return scm_i_normbig (result
);
4201 else if (SCM_REALP (y
))
4203 long int xx
= SCM_I_INUM (x
);
4204 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4206 else if (SCM_COMPLEXP (y
))
4208 long int xx
= SCM_I_INUM (x
);
4209 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4210 - SCM_COMPLEX_IMAG (y
));
4212 else if (SCM_FRACTIONP (y
))
4213 /* a - b/c = (ac - b) / c */
4214 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4215 SCM_FRACTION_NUMERATOR (y
)),
4216 SCM_FRACTION_DENOMINATOR (y
));
4218 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4220 else if (SCM_BIGP (x
))
4222 if (SCM_I_INUMP (y
))
4224 /* big-x - inum-y */
4225 long yy
= SCM_I_INUM (y
);
4226 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4228 scm_remember_upto_here_1 (x
);
4230 return (SCM_FIXABLE (-yy
) ?
4231 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4234 SCM result
= scm_i_mkbig ();
4237 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4239 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4240 scm_remember_upto_here_1 (x
);
4242 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4243 /* we know the result will have to be a bignum */
4246 return scm_i_normbig (result
);
4249 else if (SCM_BIGP (y
))
4251 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4252 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4253 SCM result
= scm_i_mkbig ();
4254 mpz_sub (SCM_I_BIG_MPZ (result
),
4257 scm_remember_upto_here_2 (x
, y
);
4258 /* we know the result will have to be a bignum */
4259 if ((sgn_x
== 1) && (sgn_y
== -1))
4261 if ((sgn_x
== -1) && (sgn_y
== 1))
4263 return scm_i_normbig (result
);
4265 else if (SCM_REALP (y
))
4267 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4268 scm_remember_upto_here_1 (x
);
4269 return scm_from_double (result
);
4271 else if (SCM_COMPLEXP (y
))
4273 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4274 - SCM_COMPLEX_REAL (y
));
4275 scm_remember_upto_here_1 (x
);
4276 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4278 else if (SCM_FRACTIONP (y
))
4279 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4280 SCM_FRACTION_NUMERATOR (y
)),
4281 SCM_FRACTION_DENOMINATOR (y
));
4282 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4284 else if (SCM_REALP (x
))
4286 if (SCM_I_INUMP (y
))
4287 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4288 else if (SCM_BIGP (y
))
4290 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4291 scm_remember_upto_here_1 (x
);
4292 return scm_from_double (result
);
4294 else if (SCM_REALP (y
))
4295 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4296 else if (SCM_COMPLEXP (y
))
4297 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4298 -SCM_COMPLEX_IMAG (y
));
4299 else if (SCM_FRACTIONP (y
))
4300 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4302 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4304 else if (SCM_COMPLEXP (x
))
4306 if (SCM_I_INUMP (y
))
4307 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4308 SCM_COMPLEX_IMAG (x
));
4309 else if (SCM_BIGP (y
))
4311 double real_part
= (SCM_COMPLEX_REAL (x
)
4312 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4313 scm_remember_upto_here_1 (x
);
4314 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4316 else if (SCM_REALP (y
))
4317 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4318 SCM_COMPLEX_IMAG (x
));
4319 else if (SCM_COMPLEXP (y
))
4320 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4321 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4322 else if (SCM_FRACTIONP (y
))
4323 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4324 SCM_COMPLEX_IMAG (x
));
4326 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4328 else if (SCM_FRACTIONP (x
))
4330 if (SCM_I_INUMP (y
))
4331 /* a/b - c = (a - cb) / b */
4332 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4333 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4334 SCM_FRACTION_DENOMINATOR (x
));
4335 else if (SCM_BIGP (y
))
4336 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4337 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4338 SCM_FRACTION_DENOMINATOR (x
));
4339 else if (SCM_REALP (y
))
4340 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4341 else if (SCM_COMPLEXP (y
))
4342 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4343 -SCM_COMPLEX_IMAG (y
));
4344 else if (SCM_FRACTIONP (y
))
4345 /* a/b - c/d = (ad - bc) / bd */
4346 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4347 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4348 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4350 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4353 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4358 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4360 "Return @math{@var{x}-1}.")
4361 #define FUNC_NAME s_scm_oneminus
4363 return scm_difference (x
, SCM_I_MAKINUM (1));
4368 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4369 /* "Return the product of all arguments. If called without arguments,\n"
4373 scm_product (SCM x
, SCM y
)
4378 return SCM_I_MAKINUM (1L);
4379 else if (SCM_NUMBERP (x
))
4382 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4385 if (SCM_I_INUMP (x
))
4390 xx
= SCM_I_INUM (x
);
4394 case 0: return x
; break;
4395 case 1: return y
; break;
4398 if (SCM_I_INUMP (y
))
4400 long yy
= SCM_I_INUM (y
);
4402 SCM k
= SCM_I_MAKINUM (kk
);
4403 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4407 SCM result
= scm_i_long2big (xx
);
4408 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4409 return scm_i_normbig (result
);
4412 else if (SCM_BIGP (y
))
4414 SCM result
= scm_i_mkbig ();
4415 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4416 scm_remember_upto_here_1 (y
);
4419 else if (SCM_REALP (y
))
4420 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4421 else if (SCM_COMPLEXP (y
))
4422 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4423 xx
* SCM_COMPLEX_IMAG (y
));
4424 else if (SCM_FRACTIONP (y
))
4425 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4426 SCM_FRACTION_DENOMINATOR (y
));
4428 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4430 else if (SCM_BIGP (x
))
4432 if (SCM_I_INUMP (y
))
4437 else if (SCM_BIGP (y
))
4439 SCM result
= scm_i_mkbig ();
4440 mpz_mul (SCM_I_BIG_MPZ (result
),
4443 scm_remember_upto_here_2 (x
, y
);
4446 else if (SCM_REALP (y
))
4448 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4449 scm_remember_upto_here_1 (x
);
4450 return scm_from_double (result
);
4452 else if (SCM_COMPLEXP (y
))
4454 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4455 scm_remember_upto_here_1 (x
);
4456 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4457 z
* SCM_COMPLEX_IMAG (y
));
4459 else if (SCM_FRACTIONP (y
))
4460 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4461 SCM_FRACTION_DENOMINATOR (y
));
4463 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4465 else if (SCM_REALP (x
))
4467 if (SCM_I_INUMP (y
))
4468 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4469 else if (SCM_BIGP (y
))
4471 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4472 scm_remember_upto_here_1 (y
);
4473 return scm_from_double (result
);
4475 else if (SCM_REALP (y
))
4476 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4477 else if (SCM_COMPLEXP (y
))
4478 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4479 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4480 else if (SCM_FRACTIONP (y
))
4481 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4483 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4485 else if (SCM_COMPLEXP (x
))
4487 if (SCM_I_INUMP (y
))
4488 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4489 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4490 else if (SCM_BIGP (y
))
4492 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4493 scm_remember_upto_here_1 (y
);
4494 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4495 z
* SCM_COMPLEX_IMAG (x
));
4497 else if (SCM_REALP (y
))
4498 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4499 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4500 else if (SCM_COMPLEXP (y
))
4502 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4503 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4504 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4505 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4507 else if (SCM_FRACTIONP (y
))
4509 double yy
= scm_i_fraction2double (y
);
4510 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4511 yy
* SCM_COMPLEX_IMAG (x
));
4514 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4516 else if (SCM_FRACTIONP (x
))
4518 if (SCM_I_INUMP (y
))
4519 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4520 SCM_FRACTION_DENOMINATOR (x
));
4521 else if (SCM_BIGP (y
))
4522 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4523 SCM_FRACTION_DENOMINATOR (x
));
4524 else if (SCM_REALP (y
))
4525 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4526 else if (SCM_COMPLEXP (y
))
4528 double xx
= scm_i_fraction2double (x
);
4529 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4530 xx
* SCM_COMPLEX_IMAG (y
));
4532 else if (SCM_FRACTIONP (y
))
4533 /* a/b * c/d = ac / bd */
4534 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4535 SCM_FRACTION_NUMERATOR (y
)),
4536 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4537 SCM_FRACTION_DENOMINATOR (y
)));
4539 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4542 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4545 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4546 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4547 #define ALLOW_DIVIDE_BY_ZERO
4548 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4551 /* The code below for complex division is adapted from the GNU
4552 libstdc++, which adapted it from f2c's libF77, and is subject to
4555 /****************************************************************
4556 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4558 Permission to use, copy, modify, and distribute this software
4559 and its documentation for any purpose and without fee is hereby
4560 granted, provided that the above copyright notice appear in all
4561 copies and that both that the copyright notice and this
4562 permission notice and warranty disclaimer appear in supporting
4563 documentation, and that the names of AT&T Bell Laboratories or
4564 Bellcore or any of their entities not be used in advertising or
4565 publicity pertaining to distribution of the software without
4566 specific, written prior permission.
4568 AT&T and Bellcore disclaim all warranties with regard to this
4569 software, including all implied warranties of merchantability
4570 and fitness. In no event shall AT&T or Bellcore be liable for
4571 any special, indirect or consequential damages or any damages
4572 whatsoever resulting from loss of use, data or profits, whether
4573 in an action of contract, negligence or other tortious action,
4574 arising out of or in connection with the use or performance of
4576 ****************************************************************/
4578 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4579 /* Divide the first argument by the product of the remaining
4580 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4582 #define FUNC_NAME s_divide
4584 scm_i_divide (SCM x
, SCM y
, int inexact
)
4591 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4592 else if (SCM_I_INUMP (x
))
4594 long xx
= SCM_I_INUM (x
);
4595 if (xx
== 1 || xx
== -1)
4597 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4599 scm_num_overflow (s_divide
);
4604 return scm_from_double (1.0 / (double) xx
);
4605 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4608 else if (SCM_BIGP (x
))
4611 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4612 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4614 else if (SCM_REALP (x
))
4616 double xx
= SCM_REAL_VALUE (x
);
4617 #ifndef ALLOW_DIVIDE_BY_ZERO
4619 scm_num_overflow (s_divide
);
4622 return scm_from_double (1.0 / xx
);
4624 else if (SCM_COMPLEXP (x
))
4626 double r
= SCM_COMPLEX_REAL (x
);
4627 double i
= SCM_COMPLEX_IMAG (x
);
4628 if (fabs(r
) <= fabs(i
))
4631 double d
= i
* (1.0 + t
* t
);
4632 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4637 double d
= r
* (1.0 + t
* t
);
4638 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4641 else if (SCM_FRACTIONP (x
))
4642 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4643 SCM_FRACTION_NUMERATOR (x
));
4645 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4648 if (SCM_I_INUMP (x
))
4650 long xx
= SCM_I_INUM (x
);
4651 if (SCM_I_INUMP (y
))
4653 long yy
= SCM_I_INUM (y
);
4656 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4657 scm_num_overflow (s_divide
);
4659 return scm_from_double ((double) xx
/ (double) yy
);
4662 else if (xx
% yy
!= 0)
4665 return scm_from_double ((double) xx
/ (double) yy
);
4666 else return scm_i_make_ratio (x
, y
);
4671 if (SCM_FIXABLE (z
))
4672 return SCM_I_MAKINUM (z
);
4674 return scm_i_long2big (z
);
4677 else if (SCM_BIGP (y
))
4680 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4681 else return scm_i_make_ratio (x
, y
);
4683 else if (SCM_REALP (y
))
4685 double yy
= SCM_REAL_VALUE (y
);
4686 #ifndef ALLOW_DIVIDE_BY_ZERO
4688 scm_num_overflow (s_divide
);
4691 return scm_from_double ((double) xx
/ yy
);
4693 else if (SCM_COMPLEXP (y
))
4696 complex_div
: /* y _must_ be a complex number */
4698 double r
= SCM_COMPLEX_REAL (y
);
4699 double i
= SCM_COMPLEX_IMAG (y
);
4700 if (fabs(r
) <= fabs(i
))
4703 double d
= i
* (1.0 + t
* t
);
4704 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4709 double d
= r
* (1.0 + t
* t
);
4710 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4714 else if (SCM_FRACTIONP (y
))
4715 /* a / b/c = ac / b */
4716 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4717 SCM_FRACTION_NUMERATOR (y
));
4719 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4721 else if (SCM_BIGP (x
))
4723 if (SCM_I_INUMP (y
))
4725 long int yy
= SCM_I_INUM (y
);
4728 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4729 scm_num_overflow (s_divide
);
4731 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4732 scm_remember_upto_here_1 (x
);
4733 return (sgn
== 0) ? scm_nan () : scm_inf ();
4740 /* FIXME: HMM, what are the relative performance issues here?
4741 We need to test. Is it faster on average to test
4742 divisible_p, then perform whichever operation, or is it
4743 faster to perform the integer div opportunistically and
4744 switch to real if there's a remainder? For now we take the
4745 middle ground: test, then if divisible, use the faster div
4748 long abs_yy
= yy
< 0 ? -yy
: yy
;
4749 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4753 SCM result
= scm_i_mkbig ();
4754 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4755 scm_remember_upto_here_1 (x
);
4757 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4758 return scm_i_normbig (result
);
4763 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4764 else return scm_i_make_ratio (x
, y
);
4768 else if (SCM_BIGP (y
))
4770 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4773 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4774 scm_num_overflow (s_divide
);
4776 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4777 scm_remember_upto_here_1 (x
);
4778 return (sgn
== 0) ? scm_nan () : scm_inf ();
4786 /* It's easily possible for the ratio x/y to fit a double
4787 but one or both x and y be too big to fit a double,
4788 hence the use of mpq_get_d rather than converting and
4791 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
4792 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
4793 return scm_from_double (mpq_get_d (q
));
4797 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4801 SCM result
= scm_i_mkbig ();
4802 mpz_divexact (SCM_I_BIG_MPZ (result
),
4805 scm_remember_upto_here_2 (x
, y
);
4806 return scm_i_normbig (result
);
4809 return scm_i_make_ratio (x
, y
);
4813 else if (SCM_REALP (y
))
4815 double yy
= SCM_REAL_VALUE (y
);
4816 #ifndef ALLOW_DIVIDE_BY_ZERO
4818 scm_num_overflow (s_divide
);
4821 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4823 else if (SCM_COMPLEXP (y
))
4825 a
= scm_i_big2dbl (x
);
4828 else if (SCM_FRACTIONP (y
))
4829 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4830 SCM_FRACTION_NUMERATOR (y
));
4832 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4834 else if (SCM_REALP (x
))
4836 double rx
= SCM_REAL_VALUE (x
);
4837 if (SCM_I_INUMP (y
))
4839 long int yy
= SCM_I_INUM (y
);
4840 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4842 scm_num_overflow (s_divide
);
4845 return scm_from_double (rx
/ (double) yy
);
4847 else if (SCM_BIGP (y
))
4849 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4850 scm_remember_upto_here_1 (y
);
4851 return scm_from_double (rx
/ dby
);
4853 else if (SCM_REALP (y
))
4855 double yy
= SCM_REAL_VALUE (y
);
4856 #ifndef ALLOW_DIVIDE_BY_ZERO
4858 scm_num_overflow (s_divide
);
4861 return scm_from_double (rx
/ yy
);
4863 else if (SCM_COMPLEXP (y
))
4868 else if (SCM_FRACTIONP (y
))
4869 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4871 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4873 else if (SCM_COMPLEXP (x
))
4875 double rx
= SCM_COMPLEX_REAL (x
);
4876 double ix
= SCM_COMPLEX_IMAG (x
);
4877 if (SCM_I_INUMP (y
))
4879 long int yy
= SCM_I_INUM (y
);
4880 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4882 scm_num_overflow (s_divide
);
4887 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4890 else if (SCM_BIGP (y
))
4892 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4893 scm_remember_upto_here_1 (y
);
4894 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4896 else if (SCM_REALP (y
))
4898 double yy
= SCM_REAL_VALUE (y
);
4899 #ifndef ALLOW_DIVIDE_BY_ZERO
4901 scm_num_overflow (s_divide
);
4904 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4906 else if (SCM_COMPLEXP (y
))
4908 double ry
= SCM_COMPLEX_REAL (y
);
4909 double iy
= SCM_COMPLEX_IMAG (y
);
4910 if (fabs(ry
) <= fabs(iy
))
4913 double d
= iy
* (1.0 + t
* t
);
4914 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4919 double d
= ry
* (1.0 + t
* t
);
4920 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4923 else if (SCM_FRACTIONP (y
))
4925 double yy
= scm_i_fraction2double (y
);
4926 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4929 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4931 else if (SCM_FRACTIONP (x
))
4933 if (SCM_I_INUMP (y
))
4935 long int yy
= SCM_I_INUM (y
);
4936 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4938 scm_num_overflow (s_divide
);
4941 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4942 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4944 else if (SCM_BIGP (y
))
4946 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4947 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4949 else if (SCM_REALP (y
))
4951 double yy
= SCM_REAL_VALUE (y
);
4952 #ifndef ALLOW_DIVIDE_BY_ZERO
4954 scm_num_overflow (s_divide
);
4957 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4959 else if (SCM_COMPLEXP (y
))
4961 a
= scm_i_fraction2double (x
);
4964 else if (SCM_FRACTIONP (y
))
4965 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4966 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4968 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4971 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4975 scm_divide (SCM x
, SCM y
)
4977 return scm_i_divide (x
, y
, 0);
4980 static SCM
scm_divide2real (SCM x
, SCM y
)
4982 return scm_i_divide (x
, y
, 1);
4988 scm_asinh (double x
)
4993 #define asinh scm_asinh
4994 return log (x
+ sqrt (x
* x
+ 1));
4997 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4998 /* "Return the inverse hyperbolic sine of @var{x}."
5003 scm_acosh (double x
)
5008 #define acosh scm_acosh
5009 return log (x
+ sqrt (x
* x
- 1));
5012 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
5013 /* "Return the inverse hyperbolic cosine of @var{x}."
5018 scm_atanh (double x
)
5023 #define atanh scm_atanh
5024 return 0.5 * log ((1 + x
) / (1 - x
));
5027 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
5028 /* "Return the inverse hyperbolic tangent of @var{x}."
5033 scm_c_truncate (double x
)
5044 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5045 half-way case (ie. when x is an integer plus 0.5) going upwards.
5046 Then half-way cases are identified and adjusted down if the
5047 round-upwards didn't give the desired even integer.
5049 "plus_half == result" identifies a half-way case. If plus_half, which is
5050 x + 0.5, is an integer then x must be an integer plus 0.5.
5052 An odd "result" value is identified with result/2 != floor(result/2).
5053 This is done with plus_half, since that value is ready for use sooner in
5054 a pipelined cpu, and we're already requiring plus_half == result.
5056 Note however that we need to be careful when x is big and already an
5057 integer. In that case "x+0.5" may round to an adjacent integer, causing
5058 us to return such a value, incorrectly. For instance if the hardware is
5059 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5060 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5061 returned. Or if the hardware is in round-upwards mode, then other bigger
5062 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5063 representable value, 2^128+2^76 (or whatever), again incorrect.
5065 These bad roundings of x+0.5 are avoided by testing at the start whether
5066 x is already an integer. If it is then clearly that's the desired result
5067 already. And if it's not then the exponent must be small enough to allow
5068 an 0.5 to be represented, and hence added without a bad rounding. */
5071 scm_c_round (double x
)
5073 double plus_half
, result
;
5078 plus_half
= x
+ 0.5;
5079 result
= floor (plus_half
);
5080 /* Adjust so that the rounding is towards even. */
5081 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5086 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5088 "Round the number @var{x} towards zero.")
5089 #define FUNC_NAME s_scm_truncate_number
5091 if (scm_is_false (scm_negative_p (x
)))
5092 return scm_floor (x
);
5094 return scm_ceiling (x
);
5098 static SCM exactly_one_half
;
5100 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5102 "Round the number @var{x} towards the nearest integer. "
5103 "When it is exactly halfway between two integers, "
5104 "round towards the even one.")
5105 #define FUNC_NAME s_scm_round_number
5107 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5109 else if (SCM_REALP (x
))
5110 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5113 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5114 single quotient+remainder division then examining to see which way
5115 the rounding should go. */
5116 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5117 SCM result
= scm_floor (plus_half
);
5118 /* Adjust so that the rounding is towards even. */
5119 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5120 && scm_is_true (scm_odd_p (result
)))
5121 return scm_difference (result
, SCM_I_MAKINUM (1));
5128 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5130 "Round the number @var{x} towards minus infinity.")
5131 #define FUNC_NAME s_scm_floor
5133 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5135 else if (SCM_REALP (x
))
5136 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5137 else if (SCM_FRACTIONP (x
))
5139 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5140 SCM_FRACTION_DENOMINATOR (x
));
5141 if (scm_is_false (scm_negative_p (x
)))
5143 /* For positive x, rounding towards zero is correct. */
5148 /* For negative x, we need to return q-1 unless x is an
5149 integer. But fractions are never integer, per our
5151 return scm_difference (q
, SCM_I_MAKINUM (1));
5155 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5159 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5161 "Round the number @var{x} towards infinity.")
5162 #define FUNC_NAME s_scm_ceiling
5164 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5166 else if (SCM_REALP (x
))
5167 return scm_from_double (ceil (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_positive_p (x
)))
5174 /* For negative x, rounding towards zero is correct. */
5179 /* For positive x, we need to return q+1 unless x is an
5180 integer. But fractions are never integer, per our
5182 return scm_sum (q
, SCM_I_MAKINUM (1));
5186 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5190 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5191 /* "Return the square root of the real number @var{x}."
5193 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5194 /* "Return the absolute value of the real number @var{x}."
5196 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5197 /* "Return the @var{x}th power of e."
5199 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5200 /* "Return the natural logarithm of the real number @var{x}."
5202 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5203 /* "Return the sine of the real number @var{x}."
5205 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5206 /* "Return the cosine of the real number @var{x}."
5208 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5209 /* "Return the tangent of the real number @var{x}."
5211 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5212 /* "Return the arc sine of the real number @var{x}."
5214 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5215 /* "Return the arc cosine of the real number @var{x}."
5217 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5218 /* "Return the arc tangent of the real number @var{x}."
5220 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5221 /* "Return the hyperbolic sine of the real number @var{x}."
5223 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5224 /* "Return the hyperbolic cosine of the real number @var{x}."
5226 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5227 /* "Return the hyperbolic tangent of the real number @var{x}."
5235 static void scm_two_doubles (SCM x
,
5237 const char *sstring
,
5241 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5243 if (SCM_I_INUMP (x
))
5244 xy
->x
= SCM_I_INUM (x
);
5245 else if (SCM_BIGP (x
))
5246 xy
->x
= scm_i_big2dbl (x
);
5247 else if (SCM_REALP (x
))
5248 xy
->x
= SCM_REAL_VALUE (x
);
5249 else if (SCM_FRACTIONP (x
))
5250 xy
->x
= scm_i_fraction2double (x
);
5252 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5254 if (SCM_I_INUMP (y
))
5255 xy
->y
= SCM_I_INUM (y
);
5256 else if (SCM_BIGP (y
))
5257 xy
->y
= scm_i_big2dbl (y
);
5258 else if (SCM_REALP (y
))
5259 xy
->y
= SCM_REAL_VALUE (y
);
5260 else if (SCM_FRACTIONP (y
))
5261 xy
->y
= scm_i_fraction2double (y
);
5263 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5267 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5269 "Return @var{x} raised to the power of @var{y}. This\n"
5270 "procedure does not accept complex arguments.")
5271 #define FUNC_NAME s_scm_sys_expt
5274 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5275 return scm_from_double (pow (xy
.x
, xy
.y
));
5280 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5282 "Return the arc tangent of the two arguments @var{x} and\n"
5283 "@var{y}. This is similar to calculating the arc tangent of\n"
5284 "@var{x} / @var{y}, except that the signs of both arguments\n"
5285 "are used to determine the quadrant of the result. This\n"
5286 "procedure does not accept complex arguments.")
5287 #define FUNC_NAME s_scm_sys_atan2
5290 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5291 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5296 scm_c_make_rectangular (double re
, double im
)
5299 return scm_from_double (re
);
5303 SCM_NEWSMOB (z
, scm_tc16_complex
,
5304 scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5306 SCM_COMPLEX_REAL (z
) = re
;
5307 SCM_COMPLEX_IMAG (z
) = im
;
5312 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5313 (SCM real
, SCM imaginary
),
5314 "Return a complex number constructed of the given @var{real} and\n"
5315 "@var{imaginary} parts.")
5316 #define FUNC_NAME s_scm_make_rectangular
5319 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5320 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5325 scm_c_make_polar (double mag
, double ang
)
5329 sincos (ang
, &s
, &c
);
5334 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5337 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5339 "Return the complex number @var{x} * e^(i * @var{y}).")
5340 #define FUNC_NAME s_scm_make_polar
5343 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5344 return scm_c_make_polar (xy
.x
, xy
.y
);
5349 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5350 /* "Return the real part of the number @var{z}."
5353 scm_real_part (SCM z
)
5355 if (SCM_I_INUMP (z
))
5357 else if (SCM_BIGP (z
))
5359 else if (SCM_REALP (z
))
5361 else if (SCM_COMPLEXP (z
))
5362 return scm_from_double (SCM_COMPLEX_REAL (z
));
5363 else if (SCM_FRACTIONP (z
))
5366 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5370 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5371 /* "Return the imaginary part of the number @var{z}."
5374 scm_imag_part (SCM z
)
5376 if (SCM_I_INUMP (z
))
5378 else if (SCM_BIGP (z
))
5380 else if (SCM_REALP (z
))
5382 else if (SCM_COMPLEXP (z
))
5383 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5384 else if (SCM_FRACTIONP (z
))
5387 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5390 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5391 /* "Return the numerator of the number @var{z}."
5394 scm_numerator (SCM z
)
5396 if (SCM_I_INUMP (z
))
5398 else if (SCM_BIGP (z
))
5400 else if (SCM_FRACTIONP (z
))
5402 scm_i_fraction_reduce (z
);
5403 return SCM_FRACTION_NUMERATOR (z
);
5405 else if (SCM_REALP (z
))
5406 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5408 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5412 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5413 /* "Return the denominator of the number @var{z}."
5416 scm_denominator (SCM z
)
5418 if (SCM_I_INUMP (z
))
5419 return SCM_I_MAKINUM (1);
5420 else if (SCM_BIGP (z
))
5421 return SCM_I_MAKINUM (1);
5422 else if (SCM_FRACTIONP (z
))
5424 scm_i_fraction_reduce (z
);
5425 return SCM_FRACTION_DENOMINATOR (z
);
5427 else if (SCM_REALP (z
))
5428 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5430 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5433 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5434 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5435 * "@code{abs} for real arguments, but also allows complex numbers."
5438 scm_magnitude (SCM z
)
5440 if (SCM_I_INUMP (z
))
5442 long int zz
= SCM_I_INUM (z
);
5445 else if (SCM_POSFIXABLE (-zz
))
5446 return SCM_I_MAKINUM (-zz
);
5448 return scm_i_long2big (-zz
);
5450 else if (SCM_BIGP (z
))
5452 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5453 scm_remember_upto_here_1 (z
);
5455 return scm_i_clonebig (z
, 0);
5459 else if (SCM_REALP (z
))
5460 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5461 else if (SCM_COMPLEXP (z
))
5462 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5463 else if (SCM_FRACTIONP (z
))
5465 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5467 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5468 SCM_FRACTION_DENOMINATOR (z
));
5471 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5475 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5476 /* "Return the angle of the complex number @var{z}."
5481 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5482 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5483 But if atan2 follows the floating point rounding mode, then the value
5484 is not a constant. Maybe it'd be close enough though. */
5485 if (SCM_I_INUMP (z
))
5487 if (SCM_I_INUM (z
) >= 0)
5490 return scm_from_double (atan2 (0.0, -1.0));
5492 else if (SCM_BIGP (z
))
5494 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5495 scm_remember_upto_here_1 (z
);
5497 return scm_from_double (atan2 (0.0, -1.0));
5501 else if (SCM_REALP (z
))
5503 if (SCM_REAL_VALUE (z
) >= 0)
5506 return scm_from_double (atan2 (0.0, -1.0));
5508 else if (SCM_COMPLEXP (z
))
5509 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5510 else if (SCM_FRACTIONP (z
))
5512 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5514 else return scm_from_double (atan2 (0.0, -1.0));
5517 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5521 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5522 /* Convert the number @var{x} to its inexact representation.\n"
5525 scm_exact_to_inexact (SCM z
)
5527 if (SCM_I_INUMP (z
))
5528 return scm_from_double ((double) SCM_I_INUM (z
));
5529 else if (SCM_BIGP (z
))
5530 return scm_from_double (scm_i_big2dbl (z
));
5531 else if (SCM_FRACTIONP (z
))
5532 return scm_from_double (scm_i_fraction2double (z
));
5533 else if (SCM_INEXACTP (z
))
5536 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5540 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5542 "Return an exact number that is numerically closest to @var{z}.")
5543 #define FUNC_NAME s_scm_inexact_to_exact
5545 if (SCM_I_INUMP (z
))
5547 else if (SCM_BIGP (z
))
5549 else if (SCM_REALP (z
))
5551 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5552 SCM_OUT_OF_RANGE (1, z
);
5559 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5560 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5561 scm_i_mpz2num (mpq_denref (frac
)));
5563 /* When scm_i_make_ratio throws, we leak the memory allocated
5570 else if (SCM_FRACTIONP (z
))
5573 SCM_WRONG_TYPE_ARG (1, z
);
5577 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5579 "Return an exact number that is within @var{err} of @var{x}.")
5580 #define FUNC_NAME s_scm_rationalize
5582 if (SCM_I_INUMP (x
))
5584 else if (SCM_BIGP (x
))
5586 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5588 /* Use continued fractions to find closest ratio. All
5589 arithmetic is done with exact numbers.
5592 SCM ex
= scm_inexact_to_exact (x
);
5593 SCM int_part
= scm_floor (ex
);
5594 SCM tt
= SCM_I_MAKINUM (1);
5595 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5596 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5600 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5603 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5604 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5606 /* We stop after a million iterations just to be absolutely sure
5607 that we don't go into an infinite loop. The process normally
5608 converges after less than a dozen iterations.
5611 err
= scm_abs (err
);
5612 while (++i
< 1000000)
5614 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5615 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5616 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5618 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5619 err
))) /* abs(x-a/b) <= err */
5621 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5622 if (scm_is_false (scm_exact_p (x
))
5623 || scm_is_false (scm_exact_p (err
)))
5624 return scm_exact_to_inexact (res
);
5628 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5630 tt
= scm_floor (rx
); /* tt = floor (rx) */
5636 scm_num_overflow (s_scm_rationalize
);
5639 SCM_WRONG_TYPE_ARG (1, x
);
5643 /* conversion functions */
5646 scm_is_integer (SCM val
)
5648 return scm_is_true (scm_integer_p (val
));
5652 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5654 if (SCM_I_INUMP (val
))
5656 scm_t_signed_bits n
= SCM_I_INUM (val
);
5657 return n
>= min
&& n
<= max
;
5659 else if (SCM_BIGP (val
))
5661 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5663 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5665 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5667 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5668 return n
>= min
&& n
<= max
;
5678 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5679 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5682 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5683 SCM_I_BIG_MPZ (val
));
5685 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5697 return n
>= min
&& n
<= max
;
5705 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5707 if (SCM_I_INUMP (val
))
5709 scm_t_signed_bits n
= SCM_I_INUM (val
);
5710 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5712 else if (SCM_BIGP (val
))
5714 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5716 else if (max
<= ULONG_MAX
)
5718 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5720 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5721 return n
>= min
&& n
<= max
;
5731 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5734 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5735 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5738 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5739 SCM_I_BIG_MPZ (val
));
5741 return n
>= min
&& n
<= max
;
5749 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5751 scm_error (scm_out_of_range_key
,
5753 "Value out of range ~S to ~S: ~S",
5754 scm_list_3 (min
, max
, bad_val
),
5755 scm_list_1 (bad_val
));
5758 #define TYPE scm_t_intmax
5759 #define TYPE_MIN min
5760 #define TYPE_MAX max
5761 #define SIZEOF_TYPE 0
5762 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5763 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5764 #include "libguile/conv-integer.i.c"
5766 #define TYPE scm_t_uintmax
5767 #define TYPE_MIN min
5768 #define TYPE_MAX max
5769 #define SIZEOF_TYPE 0
5770 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5771 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5772 #include "libguile/conv-uinteger.i.c"
5774 #define TYPE scm_t_int8
5775 #define TYPE_MIN SCM_T_INT8_MIN
5776 #define TYPE_MAX SCM_T_INT8_MAX
5777 #define SIZEOF_TYPE 1
5778 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5779 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5780 #include "libguile/conv-integer.i.c"
5782 #define TYPE scm_t_uint8
5784 #define TYPE_MAX SCM_T_UINT8_MAX
5785 #define SIZEOF_TYPE 1
5786 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5787 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5788 #include "libguile/conv-uinteger.i.c"
5790 #define TYPE scm_t_int16
5791 #define TYPE_MIN SCM_T_INT16_MIN
5792 #define TYPE_MAX SCM_T_INT16_MAX
5793 #define SIZEOF_TYPE 2
5794 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5795 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5796 #include "libguile/conv-integer.i.c"
5798 #define TYPE scm_t_uint16
5800 #define TYPE_MAX SCM_T_UINT16_MAX
5801 #define SIZEOF_TYPE 2
5802 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5803 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5804 #include "libguile/conv-uinteger.i.c"
5806 #define TYPE scm_t_int32
5807 #define TYPE_MIN SCM_T_INT32_MIN
5808 #define TYPE_MAX SCM_T_INT32_MAX
5809 #define SIZEOF_TYPE 4
5810 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5811 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5812 #include "libguile/conv-integer.i.c"
5814 #define TYPE scm_t_uint32
5816 #define TYPE_MAX SCM_T_UINT32_MAX
5817 #define SIZEOF_TYPE 4
5818 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5819 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5820 #include "libguile/conv-uinteger.i.c"
5822 #if SCM_HAVE_T_INT64
5824 #define TYPE scm_t_int64
5825 #define TYPE_MIN SCM_T_INT64_MIN
5826 #define TYPE_MAX SCM_T_INT64_MAX
5827 #define SIZEOF_TYPE 8
5828 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5829 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5830 #include "libguile/conv-integer.i.c"
5832 #define TYPE scm_t_uint64
5834 #define TYPE_MAX SCM_T_UINT64_MAX
5835 #define SIZEOF_TYPE 8
5836 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5837 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5838 #include "libguile/conv-uinteger.i.c"
5843 scm_to_mpz (SCM val
, mpz_t rop
)
5845 if (SCM_I_INUMP (val
))
5846 mpz_set_si (rop
, SCM_I_INUM (val
));
5847 else if (SCM_BIGP (val
))
5848 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5850 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5854 scm_from_mpz (mpz_t val
)
5856 return scm_i_mpz2num (val
);
5860 scm_is_real (SCM val
)
5862 return scm_is_true (scm_real_p (val
));
5866 scm_is_rational (SCM val
)
5868 return scm_is_true (scm_rational_p (val
));
5872 scm_to_double (SCM val
)
5874 if (SCM_I_INUMP (val
))
5875 return SCM_I_INUM (val
);
5876 else if (SCM_BIGP (val
))
5877 return scm_i_big2dbl (val
);
5878 else if (SCM_FRACTIONP (val
))
5879 return scm_i_fraction2double (val
);
5880 else if (SCM_REALP (val
))
5881 return SCM_REAL_VALUE (val
);
5883 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5887 scm_from_double (double val
)
5889 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5890 SCM_REAL_VALUE (z
) = val
;
5894 #if SCM_ENABLE_DISCOURAGED == 1
5897 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5901 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5905 scm_out_of_range (NULL
, num
);
5908 return scm_to_double (num
);
5912 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5916 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5920 scm_out_of_range (NULL
, num
);
5923 return scm_to_double (num
);
5929 scm_is_complex (SCM val
)
5931 return scm_is_true (scm_complex_p (val
));
5935 scm_c_real_part (SCM z
)
5937 if (SCM_COMPLEXP (z
))
5938 return SCM_COMPLEX_REAL (z
);
5941 /* Use the scm_real_part to get proper error checking and
5944 return scm_to_double (scm_real_part (z
));
5949 scm_c_imag_part (SCM z
)
5951 if (SCM_COMPLEXP (z
))
5952 return SCM_COMPLEX_IMAG (z
);
5955 /* Use the scm_imag_part to get proper error checking and
5956 dispatching. The result will almost always be 0.0, but not
5959 return scm_to_double (scm_imag_part (z
));
5964 scm_c_magnitude (SCM z
)
5966 return scm_to_double (scm_magnitude (z
));
5972 return scm_to_double (scm_angle (z
));
5976 scm_is_number (SCM z
)
5978 return scm_is_true (scm_number_p (z
));
5986 mpz_init_set_si (z_negative_one
, -1);
5988 /* It may be possible to tune the performance of some algorithms by using
5989 * the following constants to avoid the creation of bignums. Please, before
5990 * using these values, remember the two rules of program optimization:
5991 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5992 scm_c_define ("most-positive-fixnum",
5993 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5994 scm_c_define ("most-negative-fixnum",
5995 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5997 scm_add_feature ("complex");
5998 scm_add_feature ("inexact");
5999 scm_flo0
= scm_from_double (0.0);
6001 /* determine floating point precision */
6002 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6004 init_dblprec(&scm_dblprec
[i
-2],i
);
6005 init_fx_radix(fx_per_radix
[i
-2],i
);
6008 /* hard code precision for base 10 if the preprocessor tells us to... */
6009 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6012 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
6013 SCM_I_MAKINUM (2)));
6014 #include "libguile/numbers.x"