1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 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
static 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
= (*(X_CAST(double *, DINFINITY
)));
608 if (guile_Inf
== tmp
)
616 #if defined (HAVE_ISNAN)
619 /* C99 NAN, when available */
623 extern unsigned int DQNAN
[2];
624 guile_NaN
= (*(X_CAST(double *, DQNAN
)));
626 guile_NaN
= guile_Inf
/ guile_Inf
;
632 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
635 #define FUNC_NAME s_scm_inf
637 static int initialized
= 0;
643 return scm_from_double (guile_Inf
);
647 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
650 #define FUNC_NAME s_scm_nan
652 static int initialized
= 0;
658 return scm_from_double (guile_NaN
);
663 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
665 "Return the absolute value of @var{x}.")
670 long int xx
= SCM_I_INUM (x
);
673 else if (SCM_POSFIXABLE (-xx
))
674 return SCM_I_MAKINUM (-xx
);
676 return scm_i_long2big (-xx
);
678 else if (SCM_BIGP (x
))
680 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
682 return scm_i_clonebig (x
, 0);
686 else if (SCM_REALP (x
))
688 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
689 double xx
= SCM_REAL_VALUE (x
);
691 return scm_from_double (-xx
);
695 else if (SCM_FRACTIONP (x
))
697 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
699 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
700 SCM_FRACTION_DENOMINATOR (x
));
703 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
708 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
709 /* "Return the quotient of the numbers @var{x} and @var{y}."
712 scm_quotient (SCM x
, SCM y
)
716 long xx
= SCM_I_INUM (x
);
719 long yy
= SCM_I_INUM (y
);
721 scm_num_overflow (s_quotient
);
726 return SCM_I_MAKINUM (z
);
728 return scm_i_long2big (z
);
731 else if (SCM_BIGP (y
))
733 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
734 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
735 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
737 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
738 scm_remember_upto_here_1 (y
);
739 return SCM_I_MAKINUM (-1);
742 return SCM_I_MAKINUM (0);
745 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
747 else if (SCM_BIGP (x
))
751 long yy
= SCM_I_INUM (y
);
753 scm_num_overflow (s_quotient
);
758 SCM result
= scm_i_mkbig ();
761 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
764 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
767 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
768 scm_remember_upto_here_1 (x
);
769 return scm_i_normbig (result
);
772 else if (SCM_BIGP (y
))
774 SCM result
= scm_i_mkbig ();
775 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
778 scm_remember_upto_here_2 (x
, y
);
779 return scm_i_normbig (result
);
782 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
785 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
788 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
789 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
791 * "(remainder 13 4) @result{} 1\n"
792 * "(remainder -13 4) @result{} -1\n"
796 scm_remainder (SCM x
, SCM y
)
802 long yy
= SCM_I_INUM (y
);
804 scm_num_overflow (s_remainder
);
807 long z
= SCM_I_INUM (x
) % yy
;
808 return SCM_I_MAKINUM (z
);
811 else if (SCM_BIGP (y
))
813 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
814 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
815 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
817 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
818 scm_remember_upto_here_1 (y
);
819 return SCM_I_MAKINUM (0);
825 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
827 else if (SCM_BIGP (x
))
831 long yy
= SCM_I_INUM (y
);
833 scm_num_overflow (s_remainder
);
836 SCM result
= scm_i_mkbig ();
839 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
840 scm_remember_upto_here_1 (x
);
841 return scm_i_normbig (result
);
844 else if (SCM_BIGP (y
))
846 SCM result
= scm_i_mkbig ();
847 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
850 scm_remember_upto_here_2 (x
, y
);
851 return scm_i_normbig (result
);
854 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
857 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
861 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
862 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
864 * "(modulo 13 4) @result{} 1\n"
865 * "(modulo -13 4) @result{} 3\n"
869 scm_modulo (SCM x
, SCM y
)
873 long xx
= SCM_I_INUM (x
);
876 long yy
= SCM_I_INUM (y
);
878 scm_num_overflow (s_modulo
);
881 /* C99 specifies that "%" is the remainder corresponding to a
882 quotient rounded towards zero, and that's also traditional
883 for machine division, so z here should be well defined. */
901 return SCM_I_MAKINUM (result
);
904 else if (SCM_BIGP (y
))
906 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
913 SCM pos_y
= scm_i_clonebig (y
, 0);
914 /* do this after the last scm_op */
915 mpz_init_set_si (z_x
, xx
);
916 result
= pos_y
; /* re-use this bignum */
917 mpz_mod (SCM_I_BIG_MPZ (result
),
919 SCM_I_BIG_MPZ (pos_y
));
920 scm_remember_upto_here_1 (pos_y
);
924 result
= scm_i_mkbig ();
925 /* do this after the last scm_op */
926 mpz_init_set_si (z_x
, xx
);
927 mpz_mod (SCM_I_BIG_MPZ (result
),
930 scm_remember_upto_here_1 (y
);
933 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
934 mpz_add (SCM_I_BIG_MPZ (result
),
936 SCM_I_BIG_MPZ (result
));
937 scm_remember_upto_here_1 (y
);
938 /* and do this before the next one */
940 return scm_i_normbig (result
);
944 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
946 else if (SCM_BIGP (x
))
950 long yy
= SCM_I_INUM (y
);
952 scm_num_overflow (s_modulo
);
955 SCM result
= scm_i_mkbig ();
956 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
958 (yy
< 0) ? - yy
: yy
);
959 scm_remember_upto_here_1 (x
);
960 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
961 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
962 SCM_I_BIG_MPZ (result
),
964 return scm_i_normbig (result
);
967 else if (SCM_BIGP (y
))
970 SCM result
= scm_i_mkbig ();
971 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
972 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
973 mpz_mod (SCM_I_BIG_MPZ (result
),
975 SCM_I_BIG_MPZ (pos_y
));
977 scm_remember_upto_here_1 (x
);
978 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
979 mpz_add (SCM_I_BIG_MPZ (result
),
981 SCM_I_BIG_MPZ (result
));
982 scm_remember_upto_here_2 (y
, pos_y
);
983 return scm_i_normbig (result
);
987 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
990 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
993 SCM_GPROC1 (s_gcd
, "gcd", scm_tc7_asubr
, scm_gcd
, g_gcd
);
994 /* "Return the greatest common divisor of all arguments.\n"
995 * "If called without arguments, 0 is returned."
998 scm_gcd (SCM x
, SCM y
)
1001 return SCM_UNBNDP (x
) ? SCM_INUM0
: x
;
1003 if (SCM_I_INUMP (x
))
1005 if (SCM_I_INUMP (y
))
1007 long xx
= SCM_I_INUM (x
);
1008 long yy
= SCM_I_INUM (y
);
1009 long u
= xx
< 0 ? -xx
: xx
;
1010 long v
= yy
< 0 ? -yy
: yy
;
1020 /* Determine a common factor 2^k */
1021 while (!(1 & (u
| v
)))
1027 /* Now, any factor 2^n can be eliminated */
1047 return (SCM_POSFIXABLE (result
)
1048 ? SCM_I_MAKINUM (result
)
1049 : scm_i_long2big (result
));
1051 else if (SCM_BIGP (y
))
1057 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1059 else if (SCM_BIGP (x
))
1061 if (SCM_I_INUMP (y
))
1063 unsigned long result
;
1066 yy
= SCM_I_INUM (y
);
1071 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1072 scm_remember_upto_here_1 (x
);
1073 return (SCM_POSFIXABLE (result
)
1074 ? SCM_I_MAKINUM (result
)
1075 : scm_from_ulong (result
));
1077 else if (SCM_BIGP (y
))
1079 SCM result
= scm_i_mkbig ();
1080 mpz_gcd (SCM_I_BIG_MPZ (result
),
1083 scm_remember_upto_here_2 (x
, y
);
1084 return scm_i_normbig (result
);
1087 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1090 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1093 SCM_GPROC1 (s_lcm
, "lcm", scm_tc7_asubr
, scm_lcm
, g_lcm
);
1094 /* "Return the least common multiple of the arguments.\n"
1095 * "If called without arguments, 1 is returned."
1098 scm_lcm (SCM n1
, SCM n2
)
1100 if (SCM_UNBNDP (n2
))
1102 if (SCM_UNBNDP (n1
))
1103 return SCM_I_MAKINUM (1L);
1104 n2
= SCM_I_MAKINUM (1L);
1107 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1108 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1109 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1110 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1112 if (SCM_I_INUMP (n1
))
1114 if (SCM_I_INUMP (n2
))
1116 SCM d
= scm_gcd (n1
, n2
);
1117 if (scm_is_eq (d
, SCM_INUM0
))
1120 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1124 /* inum n1, big n2 */
1127 SCM result
= scm_i_mkbig ();
1128 long nn1
= SCM_I_INUM (n1
);
1129 if (nn1
== 0) return SCM_INUM0
;
1130 if (nn1
< 0) nn1
= - nn1
;
1131 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1132 scm_remember_upto_here_1 (n2
);
1140 if (SCM_I_INUMP (n2
))
1147 SCM result
= scm_i_mkbig ();
1148 mpz_lcm(SCM_I_BIG_MPZ (result
),
1150 SCM_I_BIG_MPZ (n2
));
1151 scm_remember_upto_here_2(n1
, n2
);
1152 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1158 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1163 + + + x (map digit:logand X Y)
1164 + - + x (map digit:logand X (lognot (+ -1 Y)))
1165 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1166 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1171 + + + (map digit:logior X Y)
1172 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1173 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1174 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1179 + + + (map digit:logxor X Y)
1180 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1181 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1182 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1187 + + (any digit:logand X Y)
1188 + - (any digit:logand X (lognot (+ -1 Y)))
1189 - + (any digit:logand (lognot (+ -1 X)) Y)
1194 SCM_DEFINE1 (scm_logand
, "logand", scm_tc7_asubr
,
1196 "Return the bitwise AND of the integer arguments.\n\n"
1198 "(logand) @result{} -1\n"
1199 "(logand 7) @result{} 7\n"
1200 "(logand #b111 #b011 #b001) @result{} 1\n"
1202 #define FUNC_NAME s_scm_logand
1206 if (SCM_UNBNDP (n2
))
1208 if (SCM_UNBNDP (n1
))
1209 return SCM_I_MAKINUM (-1);
1210 else if (!SCM_NUMBERP (n1
))
1211 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1212 else if (SCM_NUMBERP (n1
))
1215 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1218 if (SCM_I_INUMP (n1
))
1220 nn1
= SCM_I_INUM (n1
);
1221 if (SCM_I_INUMP (n2
))
1223 long nn2
= SCM_I_INUM (n2
);
1224 return SCM_I_MAKINUM (nn1
& nn2
);
1226 else if SCM_BIGP (n2
)
1232 SCM result_z
= scm_i_mkbig ();
1234 mpz_init_set_si (nn1_z
, nn1
);
1235 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1236 scm_remember_upto_here_1 (n2
);
1238 return scm_i_normbig (result_z
);
1242 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1244 else if (SCM_BIGP (n1
))
1246 if (SCM_I_INUMP (n2
))
1249 nn1
= SCM_I_INUM (n1
);
1252 else if (SCM_BIGP (n2
))
1254 SCM result_z
= scm_i_mkbig ();
1255 mpz_and (SCM_I_BIG_MPZ (result_z
),
1257 SCM_I_BIG_MPZ (n2
));
1258 scm_remember_upto_here_2 (n1
, n2
);
1259 return scm_i_normbig (result_z
);
1262 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1265 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1270 SCM_DEFINE1 (scm_logior
, "logior", scm_tc7_asubr
,
1272 "Return the bitwise OR of the integer arguments.\n\n"
1274 "(logior) @result{} 0\n"
1275 "(logior 7) @result{} 7\n"
1276 "(logior #b000 #b001 #b011) @result{} 3\n"
1278 #define FUNC_NAME s_scm_logior
1282 if (SCM_UNBNDP (n2
))
1284 if (SCM_UNBNDP (n1
))
1286 else if (SCM_NUMBERP (n1
))
1289 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1292 if (SCM_I_INUMP (n1
))
1294 nn1
= SCM_I_INUM (n1
);
1295 if (SCM_I_INUMP (n2
))
1297 long nn2
= SCM_I_INUM (n2
);
1298 return SCM_I_MAKINUM (nn1
| nn2
);
1300 else if (SCM_BIGP (n2
))
1306 SCM result_z
= scm_i_mkbig ();
1308 mpz_init_set_si (nn1_z
, nn1
);
1309 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1310 scm_remember_upto_here_1 (n2
);
1316 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1318 else if (SCM_BIGP (n1
))
1320 if (SCM_I_INUMP (n2
))
1323 nn1
= SCM_I_INUM (n1
);
1326 else if (SCM_BIGP (n2
))
1328 SCM result_z
= scm_i_mkbig ();
1329 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1331 SCM_I_BIG_MPZ (n2
));
1332 scm_remember_upto_here_2 (n1
, n2
);
1336 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1339 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1344 SCM_DEFINE1 (scm_logxor
, "logxor", scm_tc7_asubr
,
1346 "Return the bitwise XOR of the integer arguments. A bit is\n"
1347 "set in the result if it is set in an odd number of arguments.\n"
1349 "(logxor) @result{} 0\n"
1350 "(logxor 7) @result{} 7\n"
1351 "(logxor #b000 #b001 #b011) @result{} 2\n"
1352 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1354 #define FUNC_NAME s_scm_logxor
1358 if (SCM_UNBNDP (n2
))
1360 if (SCM_UNBNDP (n1
))
1362 else if (SCM_NUMBERP (n1
))
1365 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1368 if (SCM_I_INUMP (n1
))
1370 nn1
= SCM_I_INUM (n1
);
1371 if (SCM_I_INUMP (n2
))
1373 long nn2
= SCM_I_INUM (n2
);
1374 return SCM_I_MAKINUM (nn1
^ nn2
);
1376 else if (SCM_BIGP (n2
))
1380 SCM result_z
= scm_i_mkbig ();
1382 mpz_init_set_si (nn1_z
, nn1
);
1383 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1384 scm_remember_upto_here_1 (n2
);
1386 return scm_i_normbig (result_z
);
1390 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1392 else if (SCM_BIGP (n1
))
1394 if (SCM_I_INUMP (n2
))
1397 nn1
= SCM_I_INUM (n1
);
1400 else if (SCM_BIGP (n2
))
1402 SCM result_z
= scm_i_mkbig ();
1403 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1405 SCM_I_BIG_MPZ (n2
));
1406 scm_remember_upto_here_2 (n1
, n2
);
1407 return scm_i_normbig (result_z
);
1410 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1413 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1418 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1420 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1421 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1422 "without actually calculating the @code{logand}, just testing\n"
1426 "(logtest #b0100 #b1011) @result{} #f\n"
1427 "(logtest #b0100 #b0111) @result{} #t\n"
1429 #define FUNC_NAME s_scm_logtest
1433 if (SCM_I_INUMP (j
))
1435 nj
= SCM_I_INUM (j
);
1436 if (SCM_I_INUMP (k
))
1438 long nk
= SCM_I_INUM (k
);
1439 return scm_from_bool (nj
& nk
);
1441 else if (SCM_BIGP (k
))
1449 mpz_init_set_si (nj_z
, nj
);
1450 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1451 scm_remember_upto_here_1 (k
);
1452 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1458 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1460 else if (SCM_BIGP (j
))
1462 if (SCM_I_INUMP (k
))
1465 nj
= SCM_I_INUM (j
);
1468 else if (SCM_BIGP (k
))
1472 mpz_init (result_z
);
1476 scm_remember_upto_here_2 (j
, k
);
1477 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1478 mpz_clear (result_z
);
1482 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1485 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1490 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1492 "Test whether bit number @var{index} in @var{j} is set.\n"
1493 "@var{index} starts from 0 for the least significant bit.\n"
1496 "(logbit? 0 #b1101) @result{} #t\n"
1497 "(logbit? 1 #b1101) @result{} #f\n"
1498 "(logbit? 2 #b1101) @result{} #t\n"
1499 "(logbit? 3 #b1101) @result{} #t\n"
1500 "(logbit? 4 #b1101) @result{} #f\n"
1502 #define FUNC_NAME s_scm_logbit_p
1504 unsigned long int iindex
;
1505 iindex
= scm_to_ulong (index
);
1507 if (SCM_I_INUMP (j
))
1509 /* bits above what's in an inum follow the sign bit */
1510 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1511 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1513 else if (SCM_BIGP (j
))
1515 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1516 scm_remember_upto_here_1 (j
);
1517 return scm_from_bool (val
);
1520 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1525 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1527 "Return the integer which is the ones-complement of the integer\n"
1531 "(number->string (lognot #b10000000) 2)\n"
1532 " @result{} \"-10000001\"\n"
1533 "(number->string (lognot #b0) 2)\n"
1534 " @result{} \"-1\"\n"
1536 #define FUNC_NAME s_scm_lognot
1538 if (SCM_I_INUMP (n
)) {
1539 /* No overflow here, just need to toggle all the bits making up the inum.
1540 Enhancement: No need to strip the tag and add it back, could just xor
1541 a block of 1 bits, if that worked with the various debug versions of
1543 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1545 } else if (SCM_BIGP (n
)) {
1546 SCM result
= scm_i_mkbig ();
1547 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1548 scm_remember_upto_here_1 (n
);
1552 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1557 /* returns 0 if IN is not an integer. OUT must already be
1560 coerce_to_big (SCM in
, mpz_t out
)
1563 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1564 else if (SCM_I_INUMP (in
))
1565 mpz_set_si (out
, SCM_I_INUM (in
));
1572 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1573 (SCM n
, SCM k
, SCM m
),
1574 "Return @var{n} raised to the integer exponent\n"
1575 "@var{k}, modulo @var{m}.\n"
1578 "(modulo-expt 2 3 5)\n"
1581 #define FUNC_NAME s_scm_modulo_expt
1587 /* There are two classes of error we might encounter --
1588 1) Math errors, which we'll report by calling scm_num_overflow,
1590 2) wrong-type errors, which of course we'll report by calling
1592 We don't report those errors immediately, however; instead we do
1593 some cleanup first. These variables tell us which error (if
1594 any) we should report after cleaning up.
1596 int report_overflow
= 0;
1598 int position_of_wrong_type
= 0;
1599 SCM value_of_wrong_type
= SCM_INUM0
;
1601 SCM result
= SCM_UNDEFINED
;
1607 if (scm_is_eq (m
, SCM_INUM0
))
1609 report_overflow
= 1;
1613 if (!coerce_to_big (n
, n_tmp
))
1615 value_of_wrong_type
= n
;
1616 position_of_wrong_type
= 1;
1620 if (!coerce_to_big (k
, k_tmp
))
1622 value_of_wrong_type
= k
;
1623 position_of_wrong_type
= 2;
1627 if (!coerce_to_big (m
, m_tmp
))
1629 value_of_wrong_type
= m
;
1630 position_of_wrong_type
= 3;
1634 /* if the exponent K is negative, and we simply call mpz_powm, we
1635 will get a divide-by-zero exception when an inverse 1/n mod m
1636 doesn't exist (or is not unique). Since exceptions are hard to
1637 handle, we'll attempt the inversion "by hand" -- that way, we get
1638 a simple failure code, which is easy to handle. */
1640 if (-1 == mpz_sgn (k_tmp
))
1642 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1644 report_overflow
= 1;
1647 mpz_neg (k_tmp
, k_tmp
);
1650 result
= scm_i_mkbig ();
1651 mpz_powm (SCM_I_BIG_MPZ (result
),
1656 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1657 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1664 if (report_overflow
)
1665 scm_num_overflow (FUNC_NAME
);
1667 if (position_of_wrong_type
)
1668 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1669 value_of_wrong_type
);
1671 return scm_i_normbig (result
);
1675 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1677 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1678 "exact integer, @var{n} can be any number.\n"
1680 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1681 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1682 "includes @math{0^0} is 1.\n"
1685 "(integer-expt 2 5) @result{} 32\n"
1686 "(integer-expt -3 3) @result{} -27\n"
1687 "(integer-expt 5 -3) @result{} 1/125\n"
1688 "(integer-expt 0 0) @result{} 1\n"
1690 #define FUNC_NAME s_scm_integer_expt
1693 SCM z_i2
= SCM_BOOL_F
;
1695 SCM acc
= SCM_I_MAKINUM (1L);
1697 /* 0^0 == 1 according to R5RS */
1698 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1699 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1700 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1701 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1703 if (SCM_I_INUMP (k
))
1704 i2
= SCM_I_INUM (k
);
1705 else if (SCM_BIGP (k
))
1707 z_i2
= scm_i_clonebig (k
, 1);
1708 scm_remember_upto_here_1 (k
);
1712 SCM_WRONG_TYPE_ARG (2, k
);
1716 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1718 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1719 n
= scm_divide (n
, SCM_UNDEFINED
);
1723 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1727 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1729 return scm_product (acc
, n
);
1731 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1732 acc
= scm_product (acc
, n
);
1733 n
= scm_product (n
, n
);
1734 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1742 n
= scm_divide (n
, SCM_UNDEFINED
);
1749 return scm_product (acc
, n
);
1751 acc
= scm_product (acc
, n
);
1752 n
= scm_product (n
, n
);
1759 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1761 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1762 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1764 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1765 "@var{cnt} is negative it's a division, rounded towards negative\n"
1766 "infinity. (Note that this is not the same rounding as\n"
1767 "@code{quotient} does.)\n"
1769 "With @var{n} viewed as an infinite precision twos complement,\n"
1770 "@code{ash} means a left shift introducing zero bits, or a right\n"
1771 "shift dropping bits.\n"
1774 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1775 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1777 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1778 "(ash -23 -2) @result{} -6\n"
1780 #define FUNC_NAME s_scm_ash
1783 bits_to_shift
= scm_to_long (cnt
);
1785 if (bits_to_shift
< 0)
1787 /* Shift right by abs(cnt) bits. This is realized as a division
1788 by div:=2^abs(cnt). However, to guarantee the floor
1789 rounding, negative values require some special treatment.
1791 SCM div
= scm_integer_expt (SCM_I_MAKINUM (2),
1792 scm_from_long (-bits_to_shift
));
1794 /* scm_quotient assumes its arguments are integers, but it's legal to (ash 1/2 -1) */
1795 if (scm_is_false (scm_negative_p (n
)))
1796 return scm_quotient (n
, div
);
1798 return scm_sum (SCM_I_MAKINUM (-1L),
1799 scm_quotient (scm_sum (SCM_I_MAKINUM (1L), n
), div
));
1802 /* Shift left is done by multiplication with 2^CNT */
1803 return scm_product (n
, scm_integer_expt (SCM_I_MAKINUM (2), cnt
));
1808 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1809 (SCM n
, SCM start
, SCM end
),
1810 "Return the integer composed of the @var{start} (inclusive)\n"
1811 "through @var{end} (exclusive) bits of @var{n}. The\n"
1812 "@var{start}th bit becomes the 0-th bit in the result.\n"
1815 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1816 " @result{} \"1010\"\n"
1817 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1818 " @result{} \"10110\"\n"
1820 #define FUNC_NAME s_scm_bit_extract
1822 unsigned long int istart
, iend
, bits
;
1823 istart
= scm_to_ulong (start
);
1824 iend
= scm_to_ulong (end
);
1825 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1827 /* how many bits to keep */
1828 bits
= iend
- istart
;
1830 if (SCM_I_INUMP (n
))
1832 long int in
= SCM_I_INUM (n
);
1834 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1835 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1836 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1838 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1840 /* Since we emulate two's complement encoded numbers, this
1841 * special case requires us to produce a result that has
1842 * more bits than can be stored in a fixnum.
1844 SCM result
= scm_i_long2big (in
);
1845 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1850 /* mask down to requisite bits */
1851 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1852 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1854 else if (SCM_BIGP (n
))
1859 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1863 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1864 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1865 such bits into a ulong. */
1866 result
= scm_i_mkbig ();
1867 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1868 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1869 result
= scm_i_normbig (result
);
1871 scm_remember_upto_here_1 (n
);
1875 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1880 static const char scm_logtab
[] = {
1881 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1884 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1886 "Return the number of bits in integer @var{n}. If integer is\n"
1887 "positive, the 1-bits in its binary representation are counted.\n"
1888 "If negative, the 0-bits in its two's-complement binary\n"
1889 "representation are counted. If 0, 0 is returned.\n"
1892 "(logcount #b10101010)\n"
1899 #define FUNC_NAME s_scm_logcount
1901 if (SCM_I_INUMP (n
))
1903 unsigned long int c
= 0;
1904 long int nn
= SCM_I_INUM (n
);
1909 c
+= scm_logtab
[15 & nn
];
1912 return SCM_I_MAKINUM (c
);
1914 else if (SCM_BIGP (n
))
1916 unsigned long count
;
1917 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
1918 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
1920 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
1921 scm_remember_upto_here_1 (n
);
1922 return SCM_I_MAKINUM (count
);
1925 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1930 static const char scm_ilentab
[] = {
1931 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
1935 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
1937 "Return the number of bits necessary to represent @var{n}.\n"
1940 "(integer-length #b10101010)\n"
1942 "(integer-length 0)\n"
1944 "(integer-length #b1111)\n"
1947 #define FUNC_NAME s_scm_integer_length
1949 if (SCM_I_INUMP (n
))
1951 unsigned long int c
= 0;
1953 long int nn
= SCM_I_INUM (n
);
1959 l
= scm_ilentab
[15 & nn
];
1962 return SCM_I_MAKINUM (c
- 4 + l
);
1964 else if (SCM_BIGP (n
))
1966 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
1967 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1968 1 too big, so check for that and adjust. */
1969 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
1970 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
1971 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
1972 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
1974 scm_remember_upto_here_1 (n
);
1975 return SCM_I_MAKINUM (size
);
1978 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1982 /*** NUMBERS -> STRINGS ***/
1983 #define SCM_MAX_DBL_PREC 60
1984 #define SCM_MAX_DBL_RADIX 36
1986 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
1987 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
1988 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
1991 void init_dblprec(int *prec
, int radix
) {
1992 /* determine floating point precision by adding successively
1993 smaller increments to 1.0 until it is considered == 1.0 */
1994 double f
= ((double)1.0)/radix
;
1995 double fsum
= 1.0 + f
;
2000 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2012 void init_fx_radix(double *fx_list
, int radix
)
2014 /* initialize a per-radix list of tolerances. When added
2015 to a number < 1.0, we can determine if we should raund
2016 up and quit converting a number to a string. */
2020 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2021 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2024 /* use this array as a way to generate a single digit */
2025 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2028 idbl2str (double f
, char *a
, int radix
)
2030 int efmt
, dpt
, d
, i
, wp
;
2032 #ifdef DBL_MIN_10_EXP
2035 #endif /* DBL_MIN_10_EXP */
2040 radix
> SCM_MAX_DBL_RADIX
)
2042 /* revert to existing behavior */
2046 wp
= scm_dblprec
[radix
-2];
2047 fx
= fx_per_radix
[radix
-2];
2051 #ifdef HAVE_COPYSIGN
2052 double sgn
= copysign (1.0, f
);
2057 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2063 strcpy (a
, "-inf.0");
2065 strcpy (a
, "+inf.0");
2068 else if (xisnan (f
))
2070 strcpy (a
, "+nan.0");
2080 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2081 make-uniform-vector, from causing infinite loops. */
2082 /* just do the checking...if it passes, we do the conversion for our
2083 radix again below */
2090 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2098 while (f_cpy
> 10.0)
2101 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2122 if (f
+ fx
[wp
] >= radix
)
2129 /* adding 9999 makes this equivalent to abs(x) % 3 */
2130 dpt
= (exp
+ 9999) % 3;
2134 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2156 a
[ch
++] = number_chars
[d
];
2159 if (f
+ fx
[wp
] >= 1.0)
2161 a
[ch
- 1] = number_chars
[d
+1];
2173 if ((dpt
> 4) && (exp
> 6))
2175 d
= (a
[0] == '-' ? 2 : 1);
2176 for (i
= ch
++; i
> d
; i
--)
2189 if (a
[ch
- 1] == '.')
2190 a
[ch
++] = '0'; /* trailing zero */
2199 for (i
= radix
; i
<= exp
; i
*= radix
);
2200 for (i
/= radix
; i
; i
/= radix
)
2202 a
[ch
++] = number_chars
[exp
/ i
];
2211 icmplx2str (double real
, double imag
, char *str
, int radix
)
2215 i
= idbl2str (real
, str
, radix
);
2218 /* Don't output a '+' for negative numbers or for Inf and
2219 NaN. They will provide their own sign. */
2220 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2222 i
+= idbl2str (imag
, &str
[i
], radix
);
2229 iflo2str (SCM flt
, char *str
, int radix
)
2232 if (SCM_REALP (flt
))
2233 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2235 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2240 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2241 characters in the result.
2243 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2245 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2250 return scm_iuint2str (-num
, rad
, p
) + 1;
2253 return scm_iuint2str (num
, rad
, p
);
2256 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2257 characters in the result.
2259 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2261 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2265 scm_t_uintmax n
= num
;
2267 for (n
/= rad
; n
> 0; n
/= rad
)
2277 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2282 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2284 "Return a string holding the external representation of the\n"
2285 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2286 "inexact, a radix of 10 will be used.")
2287 #define FUNC_NAME s_scm_number_to_string
2291 if (SCM_UNBNDP (radix
))
2294 base
= scm_to_signed_integer (radix
, 2, 36);
2296 if (SCM_I_INUMP (n
))
2298 char num_buf
[SCM_INTBUFLEN
];
2299 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2300 return scm_from_locale_stringn (num_buf
, length
);
2302 else if (SCM_BIGP (n
))
2304 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2305 scm_remember_upto_here_1 (n
);
2306 return scm_take_locale_string (str
);
2308 else if (SCM_FRACTIONP (n
))
2310 scm_i_fraction_reduce (n
);
2311 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2312 scm_from_locale_string ("/"),
2313 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2315 else if (SCM_INEXACTP (n
))
2317 char num_buf
[FLOBUFLEN
];
2318 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2321 SCM_WRONG_TYPE_ARG (1, n
);
2326 /* These print routines used to be stubbed here so that scm_repl.c
2327 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2330 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2332 char num_buf
[FLOBUFLEN
];
2333 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2338 scm_i_print_double (double val
, SCM port
)
2340 char num_buf
[FLOBUFLEN
];
2341 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2345 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2348 char num_buf
[FLOBUFLEN
];
2349 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2354 scm_i_print_complex (double real
, double imag
, SCM port
)
2356 char num_buf
[FLOBUFLEN
];
2357 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2361 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2364 scm_i_fraction_reduce (sexp
);
2365 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2366 scm_lfwrite (scm_i_string_chars (str
), scm_i_string_length (str
), port
);
2367 scm_remember_upto_here_1 (str
);
2372 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2374 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2375 scm_remember_upto_here_1 (exp
);
2376 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2380 /*** END nums->strs ***/
2383 /*** STRINGS -> NUMBERS ***/
2385 /* The following functions implement the conversion from strings to numbers.
2386 * The implementation somehow follows the grammar for numbers as it is given
2387 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2388 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2389 * points should be noted about the implementation:
2390 * * Each function keeps a local index variable 'idx' that points at the
2391 * current position within the parsed string. The global index is only
2392 * updated if the function could parse the corresponding syntactic unit
2394 * * Similarly, the functions keep track of indicators of inexactness ('#',
2395 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2396 * global exactness information is only updated after each part has been
2397 * successfully parsed.
2398 * * Sequences of digits are parsed into temporary variables holding fixnums.
2399 * Only if these fixnums would overflow, the result variables are updated
2400 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2401 * the temporary variables holding the fixnums are cleared, and the process
2402 * starts over again. If for example fixnums were able to store five decimal
2403 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2404 * and the result was computed as 12345 * 100000 + 67890. In other words,
2405 * only every five digits two bignum operations were performed.
2408 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2410 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2412 /* In non ASCII-style encodings the following macro might not work. */
2413 #define XDIGIT2UINT(d) \
2414 (isdigit ((int) (unsigned char) d) \
2416 : tolower ((int) (unsigned char) d) - 'a' + 10)
2419 mem2uinteger (const char* mem
, size_t len
, unsigned int *p_idx
,
2420 unsigned int radix
, enum t_exactness
*p_exactness
)
2422 unsigned int idx
= *p_idx
;
2423 unsigned int hash_seen
= 0;
2424 scm_t_bits shift
= 1;
2426 unsigned int digit_value
;
2434 if (!isxdigit ((int) (unsigned char) c
))
2436 digit_value
= XDIGIT2UINT (c
);
2437 if (digit_value
>= radix
)
2441 result
= SCM_I_MAKINUM (digit_value
);
2445 if (isxdigit ((int) (unsigned char) c
))
2449 digit_value
= XDIGIT2UINT (c
);
2450 if (digit_value
>= radix
)
2462 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2464 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2466 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2473 shift
= shift
* radix
;
2474 add
= add
* radix
+ digit_value
;
2479 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2481 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2485 *p_exactness
= INEXACT
;
2491 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2492 * covers the parts of the rules that start at a potential point. The value
2493 * of the digits up to the point have been parsed by the caller and are given
2494 * in variable result. The content of *p_exactness indicates, whether a hash
2495 * has already been seen in the digits before the point.
2498 /* In non ASCII-style encodings the following macro might not work. */
2499 #define DIGIT2UINT(d) ((d) - '0')
2502 mem2decimal_from_point (SCM result
, const char* mem
, size_t len
,
2503 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2505 unsigned int idx
= *p_idx
;
2506 enum t_exactness x
= *p_exactness
;
2511 if (mem
[idx
] == '.')
2513 scm_t_bits shift
= 1;
2515 unsigned int digit_value
;
2516 SCM big_shift
= SCM_I_MAKINUM (1);
2522 if (isdigit ((int) (unsigned char) c
))
2527 digit_value
= DIGIT2UINT (c
);
2538 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2540 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2541 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2543 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2551 add
= add
* 10 + digit_value
;
2557 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2558 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2559 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2562 result
= scm_divide (result
, big_shift
);
2564 /* We've seen a decimal point, thus the value is implicitly inexact. */
2576 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2603 if (!isdigit ((int) (unsigned char) c
))
2607 exponent
= DIGIT2UINT (c
);
2611 if (isdigit ((int) (unsigned char) c
))
2614 if (exponent
<= SCM_MAXEXP
)
2615 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2621 if (exponent
> SCM_MAXEXP
)
2623 size_t exp_len
= idx
- start
;
2624 SCM exp_string
= scm_from_locale_stringn (&mem
[start
], exp_len
);
2625 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2626 scm_out_of_range ("string->number", exp_num
);
2629 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2631 result
= scm_product (result
, e
);
2633 result
= scm_divide2real (result
, e
);
2635 /* We've seen an exponent, thus the value is implicitly inexact. */
2653 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2656 mem2ureal (const char* mem
, size_t len
, unsigned int *p_idx
,
2657 unsigned int radix
, enum t_exactness
*p_exactness
)
2659 unsigned int idx
= *p_idx
;
2665 if (idx
+5 <= len
&& !strncmp (mem
+idx
, "inf.0", 5))
2671 if (idx
+4 < len
&& !strncmp (mem
+idx
, "nan.", 4))
2673 enum t_exactness x
= EXACT
;
2675 /* Cobble up the fractional part. We might want to set the
2676 NaN's mantissa from it. */
2678 mem2uinteger (mem
, len
, &idx
, 10, &x
);
2683 if (mem
[idx
] == '.')
2687 else if (idx
+ 1 == len
)
2689 else if (!isdigit ((int) (unsigned char) mem
[idx
+ 1]))
2692 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
, len
,
2693 p_idx
, p_exactness
);
2697 enum t_exactness x
= EXACT
;
2700 uinteger
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2701 if (scm_is_false (uinteger
))
2706 else if (mem
[idx
] == '/')
2712 divisor
= mem2uinteger (mem
, len
, &idx
, radix
, &x
);
2713 if (scm_is_false (divisor
))
2716 /* both are int/big here, I assume */
2717 result
= scm_i_make_ratio (uinteger
, divisor
);
2719 else if (radix
== 10)
2721 result
= mem2decimal_from_point (uinteger
, mem
, len
, &idx
, &x
);
2722 if (scm_is_false (result
))
2733 /* When returning an inexact zero, make sure it is represented as a
2734 floating point value so that we can change its sign.
2736 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2737 result
= scm_from_double (0.0);
2743 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2746 mem2complex (const char* mem
, size_t len
, unsigned int idx
,
2747 unsigned int radix
, enum t_exactness
*p_exactness
)
2771 ureal
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2772 if (scm_is_false (ureal
))
2774 /* input must be either +i or -i */
2779 if (mem
[idx
] == 'i' || mem
[idx
] == 'I')
2785 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2792 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2793 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2802 /* either +<ureal>i or -<ureal>i */
2809 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2812 /* polar input: <real>@<real>. */
2837 angle
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2838 if (scm_is_false (angle
))
2843 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2844 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2846 result
= scm_make_polar (ureal
, angle
);
2851 /* expecting input matching <real>[+-]<ureal>?i */
2858 int sign
= (c
== '+') ? 1 : -1;
2859 SCM imag
= mem2ureal (mem
, len
, &idx
, radix
, p_exactness
);
2861 if (scm_is_false (imag
))
2862 imag
= SCM_I_MAKINUM (sign
);
2863 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2864 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2868 if (mem
[idx
] != 'i' && mem
[idx
] != 'I')
2875 return scm_make_rectangular (ureal
, imag
);
2884 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
2886 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
2889 scm_i_mem2number (const char* mem
, size_t len
, unsigned int default_radix
)
2891 unsigned int idx
= 0;
2892 unsigned int radix
= NO_RADIX
;
2893 enum t_exactness forced_x
= NO_EXACTNESS
;
2894 enum t_exactness implicit_x
= EXACT
;
2897 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
2898 while (idx
+ 2 < len
&& mem
[idx
] == '#')
2900 switch (mem
[idx
+ 1])
2903 if (radix
!= NO_RADIX
)
2908 if (radix
!= NO_RADIX
)
2913 if (forced_x
!= NO_EXACTNESS
)
2918 if (forced_x
!= NO_EXACTNESS
)
2923 if (radix
!= NO_RADIX
)
2928 if (radix
!= NO_RADIX
)
2938 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2939 if (radix
== NO_RADIX
)
2940 result
= mem2complex (mem
, len
, idx
, default_radix
, &implicit_x
);
2942 result
= mem2complex (mem
, len
, idx
, (unsigned int) radix
, &implicit_x
);
2944 if (scm_is_false (result
))
2950 if (SCM_INEXACTP (result
))
2951 return scm_inexact_to_exact (result
);
2955 if (SCM_INEXACTP (result
))
2958 return scm_exact_to_inexact (result
);
2961 if (implicit_x
== INEXACT
)
2963 if (SCM_INEXACTP (result
))
2966 return scm_exact_to_inexact (result
);
2974 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
2975 (SCM string
, SCM radix
),
2976 "Return a number of the maximally precise representation\n"
2977 "expressed by the given @var{string}. @var{radix} must be an\n"
2978 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
2979 "is a default radix that may be overridden by an explicit radix\n"
2980 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
2981 "supplied, then the default radix is 10. If string is not a\n"
2982 "syntactically valid notation for a number, then\n"
2983 "@code{string->number} returns @code{#f}.")
2984 #define FUNC_NAME s_scm_string_to_number
2988 SCM_VALIDATE_STRING (1, string
);
2990 if (SCM_UNBNDP (radix
))
2993 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
2995 answer
= scm_i_mem2number (scm_i_string_chars (string
),
2996 scm_i_string_length (string
),
2998 scm_remember_upto_here_1 (string
);
3004 /*** END strs->nums ***/
3008 scm_bigequal (SCM x
, SCM y
)
3010 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3011 scm_remember_upto_here_2 (x
, y
);
3012 return scm_from_bool (0 == result
);
3016 scm_real_equalp (SCM x
, SCM y
)
3018 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3022 scm_complex_equalp (SCM x
, SCM y
)
3024 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3025 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3029 scm_i_fraction_equalp (SCM x
, SCM y
)
3031 scm_i_fraction_reduce (x
);
3032 scm_i_fraction_reduce (y
);
3033 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3034 SCM_FRACTION_NUMERATOR (y
)))
3035 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3036 SCM_FRACTION_DENOMINATOR (y
))))
3043 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3045 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3047 #define FUNC_NAME s_scm_number_p
3049 return scm_from_bool (SCM_NUMBERP (x
));
3053 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3055 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3056 "otherwise. Note that the sets of real, rational and integer\n"
3057 "values form subsets of the set of complex numbers, i. e. the\n"
3058 "predicate will also be fulfilled if @var{x} is a real,\n"
3059 "rational or integer number.")
3060 #define FUNC_NAME s_scm_complex_p
3062 /* all numbers are complex. */
3063 return scm_number_p (x
);
3067 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3069 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3070 "otherwise. Note that the set of integer values forms a subset of\n"
3071 "the set of real numbers, i. e. the predicate will also be\n"
3072 "fulfilled if @var{x} is an integer number.")
3073 #define FUNC_NAME s_scm_real_p
3075 /* we can't represent irrational numbers. */
3076 return scm_rational_p (x
);
3080 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3082 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3083 "otherwise. Note that the set of integer values forms a subset of\n"
3084 "the set of rational numbers, i. e. the predicate will also be\n"
3085 "fulfilled if @var{x} is an integer number.")
3086 #define FUNC_NAME s_scm_rational_p
3088 if (SCM_I_INUMP (x
))
3090 else if (SCM_IMP (x
))
3092 else if (SCM_BIGP (x
))
3094 else if (SCM_FRACTIONP (x
))
3096 else if (SCM_REALP (x
))
3097 /* due to their limited precision, all floating point numbers are
3098 rational as well. */
3105 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3107 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3109 #define FUNC_NAME s_scm_integer_p
3112 if (SCM_I_INUMP (x
))
3118 if (!SCM_INEXACTP (x
))
3120 if (SCM_COMPLEXP (x
))
3122 r
= SCM_REAL_VALUE (x
);
3123 /* +/-inf passes r==floor(r), making those #t */
3131 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3133 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3135 #define FUNC_NAME s_scm_inexact_p
3137 if (SCM_INEXACTP (x
))
3139 if (SCM_NUMBERP (x
))
3141 SCM_WRONG_TYPE_ARG (1, x
);
3146 SCM_GPROC1 (s_eq_p
, "=", scm_tc7_rpsubr
, scm_num_eq_p
, g_eq_p
);
3147 /* "Return @code{#t} if all parameters are numerically equal." */
3149 scm_num_eq_p (SCM x
, SCM y
)
3152 if (SCM_I_INUMP (x
))
3154 long xx
= SCM_I_INUM (x
);
3155 if (SCM_I_INUMP (y
))
3157 long yy
= SCM_I_INUM (y
);
3158 return scm_from_bool (xx
== yy
);
3160 else if (SCM_BIGP (y
))
3162 else if (SCM_REALP (y
))
3163 return scm_from_bool ((double) xx
== SCM_REAL_VALUE (y
));
3164 else if (SCM_COMPLEXP (y
))
3165 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3166 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3167 else if (SCM_FRACTIONP (y
))
3170 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3172 else if (SCM_BIGP (x
))
3174 if (SCM_I_INUMP (y
))
3176 else if (SCM_BIGP (y
))
3178 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3179 scm_remember_upto_here_2 (x
, y
);
3180 return scm_from_bool (0 == cmp
);
3182 else if (SCM_REALP (y
))
3185 if (xisnan (SCM_REAL_VALUE (y
)))
3187 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3188 scm_remember_upto_here_1 (x
);
3189 return scm_from_bool (0 == cmp
);
3191 else if (SCM_COMPLEXP (y
))
3194 if (0.0 != SCM_COMPLEX_IMAG (y
))
3196 if (xisnan (SCM_COMPLEX_REAL (y
)))
3198 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3199 scm_remember_upto_here_1 (x
);
3200 return scm_from_bool (0 == cmp
);
3202 else if (SCM_FRACTIONP (y
))
3205 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3207 else if (SCM_REALP (x
))
3209 if (SCM_I_INUMP (y
))
3210 return scm_from_bool (SCM_REAL_VALUE (x
) == (double) SCM_I_INUM (y
));
3211 else if (SCM_BIGP (y
))
3214 if (xisnan (SCM_REAL_VALUE (x
)))
3216 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3217 scm_remember_upto_here_1 (y
);
3218 return scm_from_bool (0 == cmp
);
3220 else if (SCM_REALP (y
))
3221 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3222 else if (SCM_COMPLEXP (y
))
3223 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3224 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3225 else if (SCM_FRACTIONP (y
))
3227 double xx
= SCM_REAL_VALUE (x
);
3231 return scm_from_bool (xx
< 0.0);
3232 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3236 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3238 else if (SCM_COMPLEXP (x
))
3240 if (SCM_I_INUMP (y
))
3241 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3242 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3243 else if (SCM_BIGP (y
))
3246 if (0.0 != SCM_COMPLEX_IMAG (x
))
3248 if (xisnan (SCM_COMPLEX_REAL (x
)))
3250 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3251 scm_remember_upto_here_1 (y
);
3252 return scm_from_bool (0 == cmp
);
3254 else if (SCM_REALP (y
))
3255 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3256 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3257 else if (SCM_COMPLEXP (y
))
3258 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3259 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3260 else if (SCM_FRACTIONP (y
))
3263 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3265 xx
= SCM_COMPLEX_REAL (x
);
3269 return scm_from_bool (xx
< 0.0);
3270 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3274 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3276 else if (SCM_FRACTIONP (x
))
3278 if (SCM_I_INUMP (y
))
3280 else if (SCM_BIGP (y
))
3282 else if (SCM_REALP (y
))
3284 double yy
= SCM_REAL_VALUE (y
);
3288 return scm_from_bool (0.0 < yy
);
3289 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3292 else if (SCM_COMPLEXP (y
))
3295 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3297 yy
= SCM_COMPLEX_REAL (y
);
3301 return scm_from_bool (0.0 < yy
);
3302 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3305 else if (SCM_FRACTIONP (y
))
3306 return scm_i_fraction_equalp (x
, y
);
3308 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARGn
, s_eq_p
);
3311 SCM_WTA_DISPATCH_2 (g_eq_p
, x
, y
, SCM_ARG1
, s_eq_p
);
3315 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3316 done are good for inums, but for bignums an answer can almost always be
3317 had by just examining a few high bits of the operands, as done by GMP in
3318 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3319 of the float exponent to take into account. */
3321 SCM_GPROC1 (s_less_p
, "<", scm_tc7_rpsubr
, scm_less_p
, g_less_p
);
3322 /* "Return @code{#t} if the list of parameters is monotonically\n"
3326 scm_less_p (SCM x
, SCM y
)
3329 if (SCM_I_INUMP (x
))
3331 long xx
= SCM_I_INUM (x
);
3332 if (SCM_I_INUMP (y
))
3334 long yy
= SCM_I_INUM (y
);
3335 return scm_from_bool (xx
< yy
);
3337 else if (SCM_BIGP (y
))
3339 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3340 scm_remember_upto_here_1 (y
);
3341 return scm_from_bool (sgn
> 0);
3343 else if (SCM_REALP (y
))
3344 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3345 else if (SCM_FRACTIONP (y
))
3347 /* "x < a/b" becomes "x*b < a" */
3349 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3350 y
= SCM_FRACTION_NUMERATOR (y
);
3354 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3356 else if (SCM_BIGP (x
))
3358 if (SCM_I_INUMP (y
))
3360 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3361 scm_remember_upto_here_1 (x
);
3362 return scm_from_bool (sgn
< 0);
3364 else if (SCM_BIGP (y
))
3366 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3367 scm_remember_upto_here_2 (x
, y
);
3368 return scm_from_bool (cmp
< 0);
3370 else if (SCM_REALP (y
))
3373 if (xisnan (SCM_REAL_VALUE (y
)))
3375 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3376 scm_remember_upto_here_1 (x
);
3377 return scm_from_bool (cmp
< 0);
3379 else if (SCM_FRACTIONP (y
))
3382 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3384 else if (SCM_REALP (x
))
3386 if (SCM_I_INUMP (y
))
3387 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3388 else if (SCM_BIGP (y
))
3391 if (xisnan (SCM_REAL_VALUE (x
)))
3393 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3394 scm_remember_upto_here_1 (y
);
3395 return scm_from_bool (cmp
> 0);
3397 else if (SCM_REALP (y
))
3398 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3399 else if (SCM_FRACTIONP (y
))
3401 double xx
= SCM_REAL_VALUE (x
);
3405 return scm_from_bool (xx
< 0.0);
3406 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3410 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3412 else if (SCM_FRACTIONP (x
))
3414 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3416 /* "a/b < y" becomes "a < y*b" */
3417 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3418 x
= SCM_FRACTION_NUMERATOR (x
);
3421 else if (SCM_REALP (y
))
3423 double yy
= SCM_REAL_VALUE (y
);
3427 return scm_from_bool (0.0 < yy
);
3428 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3431 else if (SCM_FRACTIONP (y
))
3433 /* "a/b < c/d" becomes "a*d < c*b" */
3434 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3435 SCM_FRACTION_DENOMINATOR (y
));
3436 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3437 SCM_FRACTION_DENOMINATOR (x
));
3443 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARGn
, s_less_p
);
3446 SCM_WTA_DISPATCH_2 (g_less_p
, x
, y
, SCM_ARG1
, s_less_p
);
3450 SCM_GPROC1 (s_scm_gr_p
, ">", scm_tc7_rpsubr
, scm_gr_p
, g_gr_p
);
3451 /* "Return @code{#t} if the list of parameters is monotonically\n"
3454 #define FUNC_NAME s_scm_gr_p
3456 scm_gr_p (SCM x
, SCM y
)
3458 if (!SCM_NUMBERP (x
))
3459 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3460 else if (!SCM_NUMBERP (y
))
3461 SCM_WTA_DISPATCH_2 (g_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3463 return scm_less_p (y
, x
);
3468 SCM_GPROC1 (s_scm_leq_p
, "<=", scm_tc7_rpsubr
, scm_leq_p
, g_leq_p
);
3469 /* "Return @code{#t} if the list of parameters is monotonically\n"
3472 #define FUNC_NAME s_scm_leq_p
3474 scm_leq_p (SCM x
, SCM y
)
3476 if (!SCM_NUMBERP (x
))
3477 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3478 else if (!SCM_NUMBERP (y
))
3479 SCM_WTA_DISPATCH_2 (g_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3480 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3483 return scm_not (scm_less_p (y
, x
));
3488 SCM_GPROC1 (s_scm_geq_p
, ">=", scm_tc7_rpsubr
, scm_geq_p
, g_geq_p
);
3489 /* "Return @code{#t} if the list of parameters is monotonically\n"
3492 #define FUNC_NAME s_scm_geq_p
3494 scm_geq_p (SCM x
, SCM y
)
3496 if (!SCM_NUMBERP (x
))
3497 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3498 else if (!SCM_NUMBERP (y
))
3499 SCM_WTA_DISPATCH_2 (g_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3500 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3503 return scm_not (scm_less_p (x
, y
));
3508 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3509 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3515 if (SCM_I_INUMP (z
))
3516 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3517 else if (SCM_BIGP (z
))
3519 else if (SCM_REALP (z
))
3520 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3521 else if (SCM_COMPLEXP (z
))
3522 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3523 && SCM_COMPLEX_IMAG (z
) == 0.0);
3524 else if (SCM_FRACTIONP (z
))
3527 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3531 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3532 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3536 scm_positive_p (SCM x
)
3538 if (SCM_I_INUMP (x
))
3539 return scm_from_bool (SCM_I_INUM (x
) > 0);
3540 else if (SCM_BIGP (x
))
3542 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3543 scm_remember_upto_here_1 (x
);
3544 return scm_from_bool (sgn
> 0);
3546 else if (SCM_REALP (x
))
3547 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3548 else if (SCM_FRACTIONP (x
))
3549 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3551 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3555 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3556 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3560 scm_negative_p (SCM x
)
3562 if (SCM_I_INUMP (x
))
3563 return scm_from_bool (SCM_I_INUM (x
) < 0);
3564 else if (SCM_BIGP (x
))
3566 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3567 scm_remember_upto_here_1 (x
);
3568 return scm_from_bool (sgn
< 0);
3570 else if (SCM_REALP (x
))
3571 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3572 else if (SCM_FRACTIONP (x
))
3573 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3575 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3579 /* scm_min and scm_max return an inexact when either argument is inexact, as
3580 required by r5rs. On that basis, for exact/inexact combinations the
3581 exact is converted to inexact to compare and possibly return. This is
3582 unlike scm_less_p above which takes some trouble to preserve all bits in
3583 its test, such trouble is not required for min and max. */
3585 SCM_GPROC1 (s_max
, "max", scm_tc7_asubr
, scm_max
, g_max
);
3586 /* "Return the maximum of all parameter values."
3589 scm_max (SCM x
, SCM y
)
3594 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3595 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3598 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3601 if (SCM_I_INUMP (x
))
3603 long xx
= SCM_I_INUM (x
);
3604 if (SCM_I_INUMP (y
))
3606 long yy
= SCM_I_INUM (y
);
3607 return (xx
< yy
) ? y
: x
;
3609 else if (SCM_BIGP (y
))
3611 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3612 scm_remember_upto_here_1 (y
);
3613 return (sgn
< 0) ? x
: y
;
3615 else if (SCM_REALP (y
))
3618 /* if y==NaN then ">" is false and we return NaN */
3619 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3621 else if (SCM_FRACTIONP (y
))
3624 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3627 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3629 else if (SCM_BIGP (x
))
3631 if (SCM_I_INUMP (y
))
3633 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3634 scm_remember_upto_here_1 (x
);
3635 return (sgn
< 0) ? y
: x
;
3637 else if (SCM_BIGP (y
))
3639 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3640 scm_remember_upto_here_2 (x
, y
);
3641 return (cmp
> 0) ? x
: y
;
3643 else if (SCM_REALP (y
))
3645 /* if y==NaN then xx>yy is false, so we return the NaN y */
3648 xx
= scm_i_big2dbl (x
);
3649 yy
= SCM_REAL_VALUE (y
);
3650 return (xx
> yy
? scm_from_double (xx
) : y
);
3652 else if (SCM_FRACTIONP (y
))
3657 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3659 else if (SCM_REALP (x
))
3661 if (SCM_I_INUMP (y
))
3663 double z
= SCM_I_INUM (y
);
3664 /* if x==NaN then "<" is false and we return NaN */
3665 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3667 else if (SCM_BIGP (y
))
3672 else if (SCM_REALP (y
))
3674 /* if x==NaN then our explicit check means we return NaN
3675 if y==NaN then ">" is false and we return NaN
3676 calling isnan is unavoidable, since it's the only way to know
3677 which of x or y causes any compares to be false */
3678 double xx
= SCM_REAL_VALUE (x
);
3679 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3681 else if (SCM_FRACTIONP (y
))
3683 double yy
= scm_i_fraction2double (y
);
3684 double xx
= SCM_REAL_VALUE (x
);
3685 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3688 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3690 else if (SCM_FRACTIONP (x
))
3692 if (SCM_I_INUMP (y
))
3696 else if (SCM_BIGP (y
))
3700 else if (SCM_REALP (y
))
3702 double xx
= scm_i_fraction2double (x
);
3703 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3705 else if (SCM_FRACTIONP (y
))
3710 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3713 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3717 SCM_GPROC1 (s_min
, "min", scm_tc7_asubr
, scm_min
, g_min
);
3718 /* "Return the minium of all parameter values."
3721 scm_min (SCM x
, SCM y
)
3726 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
3727 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3730 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
3733 if (SCM_I_INUMP (x
))
3735 long xx
= SCM_I_INUM (x
);
3736 if (SCM_I_INUMP (y
))
3738 long yy
= SCM_I_INUM (y
);
3739 return (xx
< yy
) ? x
: y
;
3741 else if (SCM_BIGP (y
))
3743 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3744 scm_remember_upto_here_1 (y
);
3745 return (sgn
< 0) ? y
: x
;
3747 else if (SCM_REALP (y
))
3750 /* if y==NaN then "<" is false and we return NaN */
3751 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3753 else if (SCM_FRACTIONP (y
))
3756 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
3759 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3761 else if (SCM_BIGP (x
))
3763 if (SCM_I_INUMP (y
))
3765 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3766 scm_remember_upto_here_1 (x
);
3767 return (sgn
< 0) ? x
: y
;
3769 else if (SCM_BIGP (y
))
3771 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3772 scm_remember_upto_here_2 (x
, y
);
3773 return (cmp
> 0) ? y
: x
;
3775 else if (SCM_REALP (y
))
3777 /* if y==NaN then xx<yy is false, so we return the NaN y */
3780 xx
= scm_i_big2dbl (x
);
3781 yy
= SCM_REAL_VALUE (y
);
3782 return (xx
< yy
? scm_from_double (xx
) : y
);
3784 else if (SCM_FRACTIONP (y
))
3789 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3791 else if (SCM_REALP (x
))
3793 if (SCM_I_INUMP (y
))
3795 double z
= SCM_I_INUM (y
);
3796 /* if x==NaN then "<" is false and we return NaN */
3797 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
3799 else if (SCM_BIGP (y
))
3804 else if (SCM_REALP (y
))
3806 /* if x==NaN then our explicit check means we return NaN
3807 if y==NaN then "<" is false and we return NaN
3808 calling isnan is unavoidable, since it's the only way to know
3809 which of x or y causes any compares to be false */
3810 double xx
= SCM_REAL_VALUE (x
);
3811 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
3813 else if (SCM_FRACTIONP (y
))
3815 double yy
= scm_i_fraction2double (y
);
3816 double xx
= SCM_REAL_VALUE (x
);
3817 return (yy
< xx
) ? scm_from_double (yy
) : x
;
3820 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
3822 else if (SCM_FRACTIONP (x
))
3824 if (SCM_I_INUMP (y
))
3828 else if (SCM_BIGP (y
))
3832 else if (SCM_REALP (y
))
3834 double xx
= scm_i_fraction2double (x
);
3835 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
3837 else if (SCM_FRACTIONP (y
))
3842 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3845 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
3849 SCM_GPROC1 (s_sum
, "+", scm_tc7_asubr
, scm_sum
, g_sum
);
3850 /* "Return the sum of all parameter values. Return 0 if called without\n"
3854 scm_sum (SCM x
, SCM y
)
3858 if (SCM_NUMBERP (x
)) return x
;
3859 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
3860 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
3863 if (SCM_I_INUMP (x
))
3865 if (SCM_I_INUMP (y
))
3867 long xx
= SCM_I_INUM (x
);
3868 long yy
= SCM_I_INUM (y
);
3869 long int z
= xx
+ yy
;
3870 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
3872 else if (SCM_BIGP (y
))
3877 else if (SCM_REALP (y
))
3879 long int xx
= SCM_I_INUM (x
);
3880 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
3882 else if (SCM_COMPLEXP (y
))
3884 long int xx
= SCM_I_INUM (x
);
3885 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
3886 SCM_COMPLEX_IMAG (y
));
3888 else if (SCM_FRACTIONP (y
))
3889 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3890 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3891 SCM_FRACTION_DENOMINATOR (y
));
3893 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3894 } else if (SCM_BIGP (x
))
3896 if (SCM_I_INUMP (y
))
3901 inum
= SCM_I_INUM (y
);
3904 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3907 SCM result
= scm_i_mkbig ();
3908 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
3909 scm_remember_upto_here_1 (x
);
3910 /* we know the result will have to be a bignum */
3913 return scm_i_normbig (result
);
3917 SCM result
= scm_i_mkbig ();
3918 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
3919 scm_remember_upto_here_1 (x
);
3920 /* we know the result will have to be a bignum */
3923 return scm_i_normbig (result
);
3926 else if (SCM_BIGP (y
))
3928 SCM result
= scm_i_mkbig ();
3929 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3930 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3931 mpz_add (SCM_I_BIG_MPZ (result
),
3934 scm_remember_upto_here_2 (x
, y
);
3935 /* we know the result will have to be a bignum */
3938 return scm_i_normbig (result
);
3940 else if (SCM_REALP (y
))
3942 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
3943 scm_remember_upto_here_1 (x
);
3944 return scm_from_double (result
);
3946 else if (SCM_COMPLEXP (y
))
3948 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
3949 + SCM_COMPLEX_REAL (y
));
3950 scm_remember_upto_here_1 (x
);
3951 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
3953 else if (SCM_FRACTIONP (y
))
3954 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
3955 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
3956 SCM_FRACTION_DENOMINATOR (y
));
3958 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3960 else if (SCM_REALP (x
))
3962 if (SCM_I_INUMP (y
))
3963 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
3964 else if (SCM_BIGP (y
))
3966 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
3967 scm_remember_upto_here_1 (y
);
3968 return scm_from_double (result
);
3970 else if (SCM_REALP (y
))
3971 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
3972 else if (SCM_COMPLEXP (y
))
3973 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
3974 SCM_COMPLEX_IMAG (y
));
3975 else if (SCM_FRACTIONP (y
))
3976 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
3978 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
3980 else if (SCM_COMPLEXP (x
))
3982 if (SCM_I_INUMP (y
))
3983 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
3984 SCM_COMPLEX_IMAG (x
));
3985 else if (SCM_BIGP (y
))
3987 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
3988 + SCM_COMPLEX_REAL (x
));
3989 scm_remember_upto_here_1 (y
);
3990 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
3992 else if (SCM_REALP (y
))
3993 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
3994 SCM_COMPLEX_IMAG (x
));
3995 else if (SCM_COMPLEXP (y
))
3996 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
3997 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
3998 else if (SCM_FRACTIONP (y
))
3999 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4000 SCM_COMPLEX_IMAG (x
));
4002 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4004 else if (SCM_FRACTIONP (x
))
4006 if (SCM_I_INUMP (y
))
4007 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4008 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4009 SCM_FRACTION_DENOMINATOR (x
));
4010 else if (SCM_BIGP (y
))
4011 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4012 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4013 SCM_FRACTION_DENOMINATOR (x
));
4014 else if (SCM_REALP (y
))
4015 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4016 else if (SCM_COMPLEXP (y
))
4017 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4018 SCM_COMPLEX_IMAG (y
));
4019 else if (SCM_FRACTIONP (y
))
4020 /* a/b + c/d = (ad + bc) / bd */
4021 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4022 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4023 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4025 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4028 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4032 SCM_GPROC1 (s_difference
, "-", scm_tc7_asubr
, scm_difference
, g_difference
);
4033 /* If called with one argument @var{z1}, -@var{z1} returned. Otherwise
4034 * the sum of all but the first argument are subtracted from the first
4036 #define FUNC_NAME s_difference
4038 scm_difference (SCM x
, SCM y
)
4043 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4045 if (SCM_I_INUMP (x
))
4047 long xx
= -SCM_I_INUM (x
);
4048 if (SCM_FIXABLE (xx
))
4049 return SCM_I_MAKINUM (xx
);
4051 return scm_i_long2big (xx
);
4053 else if (SCM_BIGP (x
))
4054 /* FIXME: do we really need to normalize here? */
4055 return scm_i_normbig (scm_i_clonebig (x
, 0));
4056 else if (SCM_REALP (x
))
4057 return scm_from_double (-SCM_REAL_VALUE (x
));
4058 else if (SCM_COMPLEXP (x
))
4059 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4060 -SCM_COMPLEX_IMAG (x
));
4061 else if (SCM_FRACTIONP (x
))
4062 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4063 SCM_FRACTION_DENOMINATOR (x
));
4065 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4068 if (SCM_I_INUMP (x
))
4070 if (SCM_I_INUMP (y
))
4072 long int xx
= SCM_I_INUM (x
);
4073 long int yy
= SCM_I_INUM (y
);
4074 long int z
= xx
- yy
;
4075 if (SCM_FIXABLE (z
))
4076 return SCM_I_MAKINUM (z
);
4078 return scm_i_long2big (z
);
4080 else if (SCM_BIGP (y
))
4082 /* inum-x - big-y */
4083 long xx
= SCM_I_INUM (x
);
4086 return scm_i_clonebig (y
, 0);
4089 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4090 SCM result
= scm_i_mkbig ();
4093 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4096 /* x - y == -(y + -x) */
4097 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4098 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4100 scm_remember_upto_here_1 (y
);
4102 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4103 /* we know the result will have to be a bignum */
4106 return scm_i_normbig (result
);
4109 else if (SCM_REALP (y
))
4111 long int xx
= SCM_I_INUM (x
);
4112 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4114 else if (SCM_COMPLEXP (y
))
4116 long int xx
= SCM_I_INUM (x
);
4117 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4118 - SCM_COMPLEX_IMAG (y
));
4120 else if (SCM_FRACTIONP (y
))
4121 /* a - b/c = (ac - b) / c */
4122 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4123 SCM_FRACTION_NUMERATOR (y
)),
4124 SCM_FRACTION_DENOMINATOR (y
));
4126 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4128 else if (SCM_BIGP (x
))
4130 if (SCM_I_INUMP (y
))
4132 /* big-x - inum-y */
4133 long yy
= SCM_I_INUM (y
);
4134 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4136 scm_remember_upto_here_1 (x
);
4138 return (SCM_FIXABLE (-yy
) ?
4139 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4142 SCM result
= scm_i_mkbig ();
4145 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4147 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4148 scm_remember_upto_here_1 (x
);
4150 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4151 /* we know the result will have to be a bignum */
4154 return scm_i_normbig (result
);
4157 else if (SCM_BIGP (y
))
4159 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4160 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4161 SCM result
= scm_i_mkbig ();
4162 mpz_sub (SCM_I_BIG_MPZ (result
),
4165 scm_remember_upto_here_2 (x
, y
);
4166 /* we know the result will have to be a bignum */
4167 if ((sgn_x
== 1) && (sgn_y
== -1))
4169 if ((sgn_x
== -1) && (sgn_y
== 1))
4171 return scm_i_normbig (result
);
4173 else if (SCM_REALP (y
))
4175 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4176 scm_remember_upto_here_1 (x
);
4177 return scm_from_double (result
);
4179 else if (SCM_COMPLEXP (y
))
4181 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4182 - SCM_COMPLEX_REAL (y
));
4183 scm_remember_upto_here_1 (x
);
4184 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4186 else if (SCM_FRACTIONP (y
))
4187 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4188 SCM_FRACTION_NUMERATOR (y
)),
4189 SCM_FRACTION_DENOMINATOR (y
));
4190 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4192 else if (SCM_REALP (x
))
4194 if (SCM_I_INUMP (y
))
4195 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4196 else if (SCM_BIGP (y
))
4198 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4199 scm_remember_upto_here_1 (x
);
4200 return scm_from_double (result
);
4202 else if (SCM_REALP (y
))
4203 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4204 else if (SCM_COMPLEXP (y
))
4205 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4206 -SCM_COMPLEX_IMAG (y
));
4207 else if (SCM_FRACTIONP (y
))
4208 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4210 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4212 else if (SCM_COMPLEXP (x
))
4214 if (SCM_I_INUMP (y
))
4215 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4216 SCM_COMPLEX_IMAG (x
));
4217 else if (SCM_BIGP (y
))
4219 double real_part
= (SCM_COMPLEX_REAL (x
)
4220 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4221 scm_remember_upto_here_1 (x
);
4222 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4224 else if (SCM_REALP (y
))
4225 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4226 SCM_COMPLEX_IMAG (x
));
4227 else if (SCM_COMPLEXP (y
))
4228 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4229 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4230 else if (SCM_FRACTIONP (y
))
4231 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4232 SCM_COMPLEX_IMAG (x
));
4234 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4236 else if (SCM_FRACTIONP (x
))
4238 if (SCM_I_INUMP (y
))
4239 /* a/b - c = (a - cb) / b */
4240 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4241 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4242 SCM_FRACTION_DENOMINATOR (x
));
4243 else if (SCM_BIGP (y
))
4244 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4245 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4246 SCM_FRACTION_DENOMINATOR (x
));
4247 else if (SCM_REALP (y
))
4248 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4249 else if (SCM_COMPLEXP (y
))
4250 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4251 -SCM_COMPLEX_IMAG (y
));
4252 else if (SCM_FRACTIONP (y
))
4253 /* a/b - c/d = (ad - bc) / bd */
4254 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4255 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4256 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4258 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4261 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4266 SCM_GPROC1 (s_product
, "*", scm_tc7_asubr
, scm_product
, g_product
);
4267 /* "Return the product of all arguments. If called without arguments,\n"
4271 scm_product (SCM x
, SCM y
)
4276 return SCM_I_MAKINUM (1L);
4277 else if (SCM_NUMBERP (x
))
4280 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4283 if (SCM_I_INUMP (x
))
4288 xx
= SCM_I_INUM (x
);
4292 case 0: return x
; break;
4293 case 1: return y
; break;
4296 if (SCM_I_INUMP (y
))
4298 long yy
= SCM_I_INUM (y
);
4300 SCM k
= SCM_I_MAKINUM (kk
);
4301 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4305 SCM result
= scm_i_long2big (xx
);
4306 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4307 return scm_i_normbig (result
);
4310 else if (SCM_BIGP (y
))
4312 SCM result
= scm_i_mkbig ();
4313 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4314 scm_remember_upto_here_1 (y
);
4317 else if (SCM_REALP (y
))
4318 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4319 else if (SCM_COMPLEXP (y
))
4320 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4321 xx
* SCM_COMPLEX_IMAG (y
));
4322 else if (SCM_FRACTIONP (y
))
4323 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4324 SCM_FRACTION_DENOMINATOR (y
));
4326 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4328 else if (SCM_BIGP (x
))
4330 if (SCM_I_INUMP (y
))
4335 else if (SCM_BIGP (y
))
4337 SCM result
= scm_i_mkbig ();
4338 mpz_mul (SCM_I_BIG_MPZ (result
),
4341 scm_remember_upto_here_2 (x
, y
);
4344 else if (SCM_REALP (y
))
4346 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4347 scm_remember_upto_here_1 (x
);
4348 return scm_from_double (result
);
4350 else if (SCM_COMPLEXP (y
))
4352 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4353 scm_remember_upto_here_1 (x
);
4354 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4355 z
* SCM_COMPLEX_IMAG (y
));
4357 else if (SCM_FRACTIONP (y
))
4358 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4359 SCM_FRACTION_DENOMINATOR (y
));
4361 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4363 else if (SCM_REALP (x
))
4365 if (SCM_I_INUMP (y
))
4366 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4367 else if (SCM_BIGP (y
))
4369 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4370 scm_remember_upto_here_1 (y
);
4371 return scm_from_double (result
);
4373 else if (SCM_REALP (y
))
4374 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4375 else if (SCM_COMPLEXP (y
))
4376 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4377 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4378 else if (SCM_FRACTIONP (y
))
4379 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4381 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4383 else if (SCM_COMPLEXP (x
))
4385 if (SCM_I_INUMP (y
))
4386 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4387 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4388 else if (SCM_BIGP (y
))
4390 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4391 scm_remember_upto_here_1 (y
);
4392 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4393 z
* SCM_COMPLEX_IMAG (x
));
4395 else if (SCM_REALP (y
))
4396 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4397 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4398 else if (SCM_COMPLEXP (y
))
4400 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4401 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4402 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4403 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4405 else if (SCM_FRACTIONP (y
))
4407 double yy
= scm_i_fraction2double (y
);
4408 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4409 yy
* SCM_COMPLEX_IMAG (x
));
4412 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4414 else if (SCM_FRACTIONP (x
))
4416 if (SCM_I_INUMP (y
))
4417 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4418 SCM_FRACTION_DENOMINATOR (x
));
4419 else if (SCM_BIGP (y
))
4420 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4421 SCM_FRACTION_DENOMINATOR (x
));
4422 else if (SCM_REALP (y
))
4423 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4424 else if (SCM_COMPLEXP (y
))
4426 double xx
= scm_i_fraction2double (x
);
4427 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4428 xx
* SCM_COMPLEX_IMAG (y
));
4430 else if (SCM_FRACTIONP (y
))
4431 /* a/b * c/d = ac / bd */
4432 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4433 SCM_FRACTION_NUMERATOR (y
)),
4434 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4435 SCM_FRACTION_DENOMINATOR (y
)));
4437 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4440 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4443 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4444 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4445 #define ALLOW_DIVIDE_BY_ZERO
4446 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4449 /* The code below for complex division is adapted from the GNU
4450 libstdc++, which adapted it from f2c's libF77, and is subject to
4453 /****************************************************************
4454 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4456 Permission to use, copy, modify, and distribute this software
4457 and its documentation for any purpose and without fee is hereby
4458 granted, provided that the above copyright notice appear in all
4459 copies and that both that the copyright notice and this
4460 permission notice and warranty disclaimer appear in supporting
4461 documentation, and that the names of AT&T Bell Laboratories or
4462 Bellcore or any of their entities not be used in advertising or
4463 publicity pertaining to distribution of the software without
4464 specific, written prior permission.
4466 AT&T and Bellcore disclaim all warranties with regard to this
4467 software, including all implied warranties of merchantability
4468 and fitness. In no event shall AT&T or Bellcore be liable for
4469 any special, indirect or consequential damages or any damages
4470 whatsoever resulting from loss of use, data or profits, whether
4471 in an action of contract, negligence or other tortious action,
4472 arising out of or in connection with the use or performance of
4474 ****************************************************************/
4476 SCM_GPROC1 (s_divide
, "/", scm_tc7_asubr
, scm_divide
, g_divide
);
4477 /* Divide the first argument by the product of the remaining
4478 arguments. If called with one argument @var{z1}, 1/@var{z1} is
4480 #define FUNC_NAME s_divide
4482 scm_i_divide (SCM x
, SCM y
, int inexact
)
4489 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4490 else if (SCM_I_INUMP (x
))
4492 long xx
= SCM_I_INUM (x
);
4493 if (xx
== 1 || xx
== -1)
4495 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4497 scm_num_overflow (s_divide
);
4502 return scm_from_double (1.0 / (double) xx
);
4503 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4506 else if (SCM_BIGP (x
))
4509 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4510 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4512 else if (SCM_REALP (x
))
4514 double xx
= SCM_REAL_VALUE (x
);
4515 #ifndef ALLOW_DIVIDE_BY_ZERO
4517 scm_num_overflow (s_divide
);
4520 return scm_from_double (1.0 / xx
);
4522 else if (SCM_COMPLEXP (x
))
4524 double r
= SCM_COMPLEX_REAL (x
);
4525 double i
= SCM_COMPLEX_IMAG (x
);
4529 double d
= i
* (1.0 + t
* t
);
4530 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4535 double d
= r
* (1.0 + t
* t
);
4536 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4539 else if (SCM_FRACTIONP (x
))
4540 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4541 SCM_FRACTION_NUMERATOR (x
));
4543 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4546 if (SCM_I_INUMP (x
))
4548 long xx
= SCM_I_INUM (x
);
4549 if (SCM_I_INUMP (y
))
4551 long yy
= SCM_I_INUM (y
);
4554 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4555 scm_num_overflow (s_divide
);
4557 return scm_from_double ((double) xx
/ (double) yy
);
4560 else if (xx
% yy
!= 0)
4563 return scm_from_double ((double) xx
/ (double) yy
);
4564 else return scm_i_make_ratio (x
, y
);
4569 if (SCM_FIXABLE (z
))
4570 return SCM_I_MAKINUM (z
);
4572 return scm_i_long2big (z
);
4575 else if (SCM_BIGP (y
))
4578 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4579 else return scm_i_make_ratio (x
, y
);
4581 else if (SCM_REALP (y
))
4583 double yy
= SCM_REAL_VALUE (y
);
4584 #ifndef ALLOW_DIVIDE_BY_ZERO
4586 scm_num_overflow (s_divide
);
4589 return scm_from_double ((double) xx
/ yy
);
4591 else if (SCM_COMPLEXP (y
))
4594 complex_div
: /* y _must_ be a complex number */
4596 double r
= SCM_COMPLEX_REAL (y
);
4597 double i
= SCM_COMPLEX_IMAG (y
);
4601 double d
= i
* (1.0 + t
* t
);
4602 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4607 double d
= r
* (1.0 + t
* t
);
4608 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4612 else if (SCM_FRACTIONP (y
))
4613 /* a / b/c = ac / b */
4614 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4615 SCM_FRACTION_NUMERATOR (y
));
4617 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4619 else if (SCM_BIGP (x
))
4621 if (SCM_I_INUMP (y
))
4623 long int yy
= SCM_I_INUM (y
);
4626 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4627 scm_num_overflow (s_divide
);
4629 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4630 scm_remember_upto_here_1 (x
);
4631 return (sgn
== 0) ? scm_nan () : scm_inf ();
4638 /* FIXME: HMM, what are the relative performance issues here?
4639 We need to test. Is it faster on average to test
4640 divisible_p, then perform whichever operation, or is it
4641 faster to perform the integer div opportunistically and
4642 switch to real if there's a remainder? For now we take the
4643 middle ground: test, then if divisible, use the faster div
4646 long abs_yy
= yy
< 0 ? -yy
: yy
;
4647 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
4651 SCM result
= scm_i_mkbig ();
4652 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
4653 scm_remember_upto_here_1 (x
);
4655 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4656 return scm_i_normbig (result
);
4661 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
4662 else return scm_i_make_ratio (x
, y
);
4666 else if (SCM_BIGP (y
))
4668 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
4671 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4672 scm_num_overflow (s_divide
);
4674 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4675 scm_remember_upto_here_1 (x
);
4676 return (sgn
== 0) ? scm_nan () : scm_inf ();
4682 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
4686 SCM result
= scm_i_mkbig ();
4687 mpz_divexact (SCM_I_BIG_MPZ (result
),
4690 scm_remember_upto_here_2 (x
, y
);
4691 return scm_i_normbig (result
);
4697 double dbx
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4698 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4699 scm_remember_upto_here_2 (x
, y
);
4700 return scm_from_double (dbx
/ dby
);
4702 else return scm_i_make_ratio (x
, y
);
4706 else if (SCM_REALP (y
))
4708 double yy
= SCM_REAL_VALUE (y
);
4709 #ifndef ALLOW_DIVIDE_BY_ZERO
4711 scm_num_overflow (s_divide
);
4714 return scm_from_double (scm_i_big2dbl (x
) / yy
);
4716 else if (SCM_COMPLEXP (y
))
4718 a
= scm_i_big2dbl (x
);
4721 else if (SCM_FRACTIONP (y
))
4722 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4723 SCM_FRACTION_NUMERATOR (y
));
4725 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4727 else if (SCM_REALP (x
))
4729 double rx
= SCM_REAL_VALUE (x
);
4730 if (SCM_I_INUMP (y
))
4732 long int yy
= SCM_I_INUM (y
);
4733 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4735 scm_num_overflow (s_divide
);
4738 return scm_from_double (rx
/ (double) yy
);
4740 else if (SCM_BIGP (y
))
4742 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4743 scm_remember_upto_here_1 (y
);
4744 return scm_from_double (rx
/ dby
);
4746 else if (SCM_REALP (y
))
4748 double yy
= SCM_REAL_VALUE (y
);
4749 #ifndef ALLOW_DIVIDE_BY_ZERO
4751 scm_num_overflow (s_divide
);
4754 return scm_from_double (rx
/ yy
);
4756 else if (SCM_COMPLEXP (y
))
4761 else if (SCM_FRACTIONP (y
))
4762 return scm_from_double (rx
/ scm_i_fraction2double (y
));
4764 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4766 else if (SCM_COMPLEXP (x
))
4768 double rx
= SCM_COMPLEX_REAL (x
);
4769 double ix
= SCM_COMPLEX_IMAG (x
);
4770 if (SCM_I_INUMP (y
))
4772 long int yy
= SCM_I_INUM (y
);
4773 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4775 scm_num_overflow (s_divide
);
4780 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
4783 else if (SCM_BIGP (y
))
4785 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4786 scm_remember_upto_here_1 (y
);
4787 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
4789 else if (SCM_REALP (y
))
4791 double yy
= SCM_REAL_VALUE (y
);
4792 #ifndef ALLOW_DIVIDE_BY_ZERO
4794 scm_num_overflow (s_divide
);
4797 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4799 else if (SCM_COMPLEXP (y
))
4801 double ry
= SCM_COMPLEX_REAL (y
);
4802 double iy
= SCM_COMPLEX_IMAG (y
);
4806 double d
= iy
* (1.0 + t
* t
);
4807 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
4812 double d
= ry
* (1.0 + t
* t
);
4813 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
4816 else if (SCM_FRACTIONP (y
))
4818 double yy
= scm_i_fraction2double (y
);
4819 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
4822 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4824 else if (SCM_FRACTIONP (x
))
4826 if (SCM_I_INUMP (y
))
4828 long int yy
= SCM_I_INUM (y
);
4829 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4831 scm_num_overflow (s_divide
);
4834 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4835 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4837 else if (SCM_BIGP (y
))
4839 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
4840 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
4842 else if (SCM_REALP (y
))
4844 double yy
= SCM_REAL_VALUE (y
);
4845 #ifndef ALLOW_DIVIDE_BY_ZERO
4847 scm_num_overflow (s_divide
);
4850 return scm_from_double (scm_i_fraction2double (x
) / yy
);
4852 else if (SCM_COMPLEXP (y
))
4854 a
= scm_i_fraction2double (x
);
4857 else if (SCM_FRACTIONP (y
))
4858 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4859 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
4861 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4864 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
4868 scm_divide (SCM x
, SCM y
)
4870 return scm_i_divide (x
, y
, 0);
4873 static SCM
scm_divide2real (SCM x
, SCM y
)
4875 return scm_i_divide (x
, y
, 1);
4881 scm_asinh (double x
)
4886 #define asinh scm_asinh
4887 return log (x
+ sqrt (x
* x
+ 1));
4890 SCM_GPROC1 (s_asinh
, "$asinh", scm_tc7_dsubr
, (SCM (*)()) asinh
, g_asinh
);
4891 /* "Return the inverse hyperbolic sine of @var{x}."
4896 scm_acosh (double x
)
4901 #define acosh scm_acosh
4902 return log (x
+ sqrt (x
* x
- 1));
4905 SCM_GPROC1 (s_acosh
, "$acosh", scm_tc7_dsubr
, (SCM (*)()) acosh
, g_acosh
);
4906 /* "Return the inverse hyperbolic cosine of @var{x}."
4911 scm_atanh (double x
)
4916 #define atanh scm_atanh
4917 return 0.5 * log ((1 + x
) / (1 - x
));
4920 SCM_GPROC1 (s_atanh
, "$atanh", scm_tc7_dsubr
, (SCM (*)()) atanh
, g_atanh
);
4921 /* "Return the inverse hyperbolic tangent of @var{x}."
4926 scm_c_truncate (double x
)
4937 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
4938 half-way case (ie. when x is an integer plus 0.5) going upwards.
4939 Then half-way cases are identified and adjusted down if the
4940 round-upwards didn't give the desired even integer.
4942 "plus_half == result" identifies a half-way case. If plus_half, which is
4943 x + 0.5, is an integer then x must be an integer plus 0.5.
4945 An odd "result" value is identified with result/2 != floor(result/2).
4946 This is done with plus_half, since that value is ready for use sooner in
4947 a pipelined cpu, and we're already requiring plus_half == result.
4949 Note however that we need to be careful when x is big and already an
4950 integer. In that case "x+0.5" may round to an adjacent integer, causing
4951 us to return such a value, incorrectly. For instance if the hardware is
4952 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
4953 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
4954 returned. Or if the hardware is in round-upwards mode, then other bigger
4955 values like say x == 2^128 will see x+0.5 rounding up to the next higher
4956 representable value, 2^128+2^76 (or whatever), again incorrect.
4958 These bad roundings of x+0.5 are avoided by testing at the start whether
4959 x is already an integer. If it is then clearly that's the desired result
4960 already. And if it's not then the exponent must be small enough to allow
4961 an 0.5 to be represented, and hence added without a bad rounding. */
4964 scm_c_round (double x
)
4966 double plus_half
, result
;
4971 plus_half
= x
+ 0.5;
4972 result
= floor (plus_half
);
4973 /* Adjust so that the rounding is towards even. */
4974 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
4979 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
4981 "Round the number @var{x} towards zero.")
4982 #define FUNC_NAME s_scm_truncate_number
4984 if (scm_is_false (scm_negative_p (x
)))
4985 return scm_floor (x
);
4987 return scm_ceiling (x
);
4991 static SCM exactly_one_half
;
4993 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
4995 "Round the number @var{x} towards the nearest integer. "
4996 "When it is exactly halfway between two integers, "
4997 "round towards the even one.")
4998 #define FUNC_NAME s_scm_round_number
5000 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5002 else if (SCM_REALP (x
))
5003 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5006 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5007 single quotient+remainder division then examining to see which way
5008 the rounding should go. */
5009 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5010 SCM result
= scm_floor (plus_half
);
5011 /* Adjust so that the rounding is towards even. */
5012 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5013 && scm_is_true (scm_odd_p (result
)))
5014 return scm_difference (result
, SCM_I_MAKINUM (1));
5021 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5023 "Round the number @var{x} towards minus infinity.")
5024 #define FUNC_NAME s_scm_floor
5026 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5028 else if (SCM_REALP (x
))
5029 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5030 else if (SCM_FRACTIONP (x
))
5032 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5033 SCM_FRACTION_DENOMINATOR (x
));
5034 if (scm_is_false (scm_negative_p (x
)))
5036 /* For positive x, rounding towards zero is correct. */
5041 /* For negative x, we need to return q-1 unless x is an
5042 integer. But fractions are never integer, per our
5044 return scm_difference (q
, SCM_I_MAKINUM (1));
5048 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5052 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5054 "Round the number @var{x} towards infinity.")
5055 #define FUNC_NAME s_scm_ceiling
5057 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5059 else if (SCM_REALP (x
))
5060 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5061 else if (SCM_FRACTIONP (x
))
5063 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5064 SCM_FRACTION_DENOMINATOR (x
));
5065 if (scm_is_false (scm_positive_p (x
)))
5067 /* For negative x, rounding towards zero is correct. */
5072 /* For positive x, we need to return q+1 unless x is an
5073 integer. But fractions are never integer, per our
5075 return scm_sum (q
, SCM_I_MAKINUM (1));
5079 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5083 SCM_GPROC1 (s_i_sqrt
, "$sqrt", scm_tc7_dsubr
, (SCM (*)()) sqrt
, g_i_sqrt
);
5084 /* "Return the square root of the real number @var{x}."
5086 SCM_GPROC1 (s_i_abs
, "$abs", scm_tc7_dsubr
, (SCM (*)()) fabs
, g_i_abs
);
5087 /* "Return the absolute value of the real number @var{x}."
5089 SCM_GPROC1 (s_i_exp
, "$exp", scm_tc7_dsubr
, (SCM (*)()) exp
, g_i_exp
);
5090 /* "Return the @var{x}th power of e."
5092 SCM_GPROC1 (s_i_log
, "$log", scm_tc7_dsubr
, (SCM (*)()) log
, g_i_log
);
5093 /* "Return the natural logarithm of the real number @var{x}."
5095 SCM_GPROC1 (s_i_sin
, "$sin", scm_tc7_dsubr
, (SCM (*)()) sin
, g_i_sin
);
5096 /* "Return the sine of the real number @var{x}."
5098 SCM_GPROC1 (s_i_cos
, "$cos", scm_tc7_dsubr
, (SCM (*)()) cos
, g_i_cos
);
5099 /* "Return the cosine of the real number @var{x}."
5101 SCM_GPROC1 (s_i_tan
, "$tan", scm_tc7_dsubr
, (SCM (*)()) tan
, g_i_tan
);
5102 /* "Return the tangent of the real number @var{x}."
5104 SCM_GPROC1 (s_i_asin
, "$asin", scm_tc7_dsubr
, (SCM (*)()) asin
, g_i_asin
);
5105 /* "Return the arc sine of the real number @var{x}."
5107 SCM_GPROC1 (s_i_acos
, "$acos", scm_tc7_dsubr
, (SCM (*)()) acos
, g_i_acos
);
5108 /* "Return the arc cosine of the real number @var{x}."
5110 SCM_GPROC1 (s_i_atan
, "$atan", scm_tc7_dsubr
, (SCM (*)()) atan
, g_i_atan
);
5111 /* "Return the arc tangent of the real number @var{x}."
5113 SCM_GPROC1 (s_i_sinh
, "$sinh", scm_tc7_dsubr
, (SCM (*)()) sinh
, g_i_sinh
);
5114 /* "Return the hyperbolic sine of the real number @var{x}."
5116 SCM_GPROC1 (s_i_cosh
, "$cosh", scm_tc7_dsubr
, (SCM (*)()) cosh
, g_i_cosh
);
5117 /* "Return the hyperbolic cosine of the real number @var{x}."
5119 SCM_GPROC1 (s_i_tanh
, "$tanh", scm_tc7_dsubr
, (SCM (*)()) tanh
, g_i_tanh
);
5120 /* "Return the hyperbolic tangent of the real number @var{x}."
5128 static void scm_two_doubles (SCM x
,
5130 const char *sstring
,
5134 scm_two_doubles (SCM x
, SCM y
, const char *sstring
, struct dpair
*xy
)
5136 if (SCM_I_INUMP (x
))
5137 xy
->x
= SCM_I_INUM (x
);
5138 else if (SCM_BIGP (x
))
5139 xy
->x
= scm_i_big2dbl (x
);
5140 else if (SCM_REALP (x
))
5141 xy
->x
= SCM_REAL_VALUE (x
);
5142 else if (SCM_FRACTIONP (x
))
5143 xy
->x
= scm_i_fraction2double (x
);
5145 scm_wrong_type_arg (sstring
, SCM_ARG1
, x
);
5147 if (SCM_I_INUMP (y
))
5148 xy
->y
= SCM_I_INUM (y
);
5149 else if (SCM_BIGP (y
))
5150 xy
->y
= scm_i_big2dbl (y
);
5151 else if (SCM_REALP (y
))
5152 xy
->y
= SCM_REAL_VALUE (y
);
5153 else if (SCM_FRACTIONP (y
))
5154 xy
->y
= scm_i_fraction2double (y
);
5156 scm_wrong_type_arg (sstring
, SCM_ARG2
, y
);
5160 SCM_DEFINE (scm_sys_expt
, "$expt", 2, 0, 0,
5162 "Return @var{x} raised to the power of @var{y}. This\n"
5163 "procedure does not accept complex arguments.")
5164 #define FUNC_NAME s_scm_sys_expt
5167 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5168 return scm_from_double (pow (xy
.x
, xy
.y
));
5173 SCM_DEFINE (scm_sys_atan2
, "$atan2", 2, 0, 0,
5175 "Return the arc tangent of the two arguments @var{x} and\n"
5176 "@var{y}. This is similar to calculating the arc tangent of\n"
5177 "@var{x} / @var{y}, except that the signs of both arguments\n"
5178 "are used to determine the quadrant of the result. This\n"
5179 "procedure does not accept complex arguments.")
5180 #define FUNC_NAME s_scm_sys_atan2
5183 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5184 return scm_from_double (atan2 (xy
.x
, xy
.y
));
5189 scm_c_make_rectangular (double re
, double im
)
5192 return scm_from_double (re
);
5196 SCM_NEWSMOB (z
, scm_tc16_complex
, scm_gc_malloc (sizeof (scm_t_complex
),
5198 SCM_COMPLEX_REAL (z
) = re
;
5199 SCM_COMPLEX_IMAG (z
) = im
;
5204 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5205 (SCM real
, SCM imaginary
),
5206 "Return a complex number constructed of the given @var{real} and\n"
5207 "@var{imaginary} parts.")
5208 #define FUNC_NAME s_scm_make_rectangular
5211 scm_two_doubles (real
, imaginary
, FUNC_NAME
, &xy
);
5212 return scm_c_make_rectangular (xy
.x
, xy
.y
);
5217 scm_c_make_polar (double mag
, double ang
)
5221 sincos (ang
, &s
, &c
);
5226 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5229 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5231 "Return the complex number @var{x} * e^(i * @var{y}).")
5232 #define FUNC_NAME s_scm_make_polar
5235 scm_two_doubles (x
, y
, FUNC_NAME
, &xy
);
5236 return scm_c_make_polar (xy
.x
, xy
.y
);
5241 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5242 /* "Return the real part of the number @var{z}."
5245 scm_real_part (SCM z
)
5247 if (SCM_I_INUMP (z
))
5249 else if (SCM_BIGP (z
))
5251 else if (SCM_REALP (z
))
5253 else if (SCM_COMPLEXP (z
))
5254 return scm_from_double (SCM_COMPLEX_REAL (z
));
5255 else if (SCM_FRACTIONP (z
))
5258 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5262 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5263 /* "Return the imaginary part of the number @var{z}."
5266 scm_imag_part (SCM z
)
5268 if (SCM_I_INUMP (z
))
5270 else if (SCM_BIGP (z
))
5272 else if (SCM_REALP (z
))
5274 else if (SCM_COMPLEXP (z
))
5275 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5276 else if (SCM_FRACTIONP (z
))
5279 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5282 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5283 /* "Return the numerator of the number @var{z}."
5286 scm_numerator (SCM z
)
5288 if (SCM_I_INUMP (z
))
5290 else if (SCM_BIGP (z
))
5292 else if (SCM_FRACTIONP (z
))
5294 scm_i_fraction_reduce (z
);
5295 return SCM_FRACTION_NUMERATOR (z
);
5297 else if (SCM_REALP (z
))
5298 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5300 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5304 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5305 /* "Return the denominator of the number @var{z}."
5308 scm_denominator (SCM z
)
5310 if (SCM_I_INUMP (z
))
5311 return SCM_I_MAKINUM (1);
5312 else if (SCM_BIGP (z
))
5313 return SCM_I_MAKINUM (1);
5314 else if (SCM_FRACTIONP (z
))
5316 scm_i_fraction_reduce (z
);
5317 return SCM_FRACTION_DENOMINATOR (z
);
5319 else if (SCM_REALP (z
))
5320 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5322 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5325 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5326 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5327 * "@code{abs} for real arguments, but also allows complex numbers."
5330 scm_magnitude (SCM z
)
5332 if (SCM_I_INUMP (z
))
5334 long int zz
= SCM_I_INUM (z
);
5337 else if (SCM_POSFIXABLE (-zz
))
5338 return SCM_I_MAKINUM (-zz
);
5340 return scm_i_long2big (-zz
);
5342 else if (SCM_BIGP (z
))
5344 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5345 scm_remember_upto_here_1 (z
);
5347 return scm_i_clonebig (z
, 0);
5351 else if (SCM_REALP (z
))
5352 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5353 else if (SCM_COMPLEXP (z
))
5354 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5355 else if (SCM_FRACTIONP (z
))
5357 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5359 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5360 SCM_FRACTION_DENOMINATOR (z
));
5363 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5367 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5368 /* "Return the angle of the complex number @var{z}."
5373 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5374 scm_flo0 to save allocating a new flonum with scm_from_double each time.
5375 But if atan2 follows the floating point rounding mode, then the value
5376 is not a constant. Maybe it'd be close enough though. */
5377 if (SCM_I_INUMP (z
))
5379 if (SCM_I_INUM (z
) >= 0)
5382 return scm_from_double (atan2 (0.0, -1.0));
5384 else if (SCM_BIGP (z
))
5386 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5387 scm_remember_upto_here_1 (z
);
5389 return scm_from_double (atan2 (0.0, -1.0));
5393 else if (SCM_REALP (z
))
5395 if (SCM_REAL_VALUE (z
) >= 0)
5398 return scm_from_double (atan2 (0.0, -1.0));
5400 else if (SCM_COMPLEXP (z
))
5401 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5402 else if (SCM_FRACTIONP (z
))
5404 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5406 else return scm_from_double (atan2 (0.0, -1.0));
5409 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5413 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5414 /* Convert the number @var{x} to its inexact representation.\n"
5417 scm_exact_to_inexact (SCM z
)
5419 if (SCM_I_INUMP (z
))
5420 return scm_from_double ((double) SCM_I_INUM (z
));
5421 else if (SCM_BIGP (z
))
5422 return scm_from_double (scm_i_big2dbl (z
));
5423 else if (SCM_FRACTIONP (z
))
5424 return scm_from_double (scm_i_fraction2double (z
));
5425 else if (SCM_INEXACTP (z
))
5428 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5432 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5434 "Return an exact number that is numerically closest to @var{z}.")
5435 #define FUNC_NAME s_scm_inexact_to_exact
5437 if (SCM_I_INUMP (z
))
5439 else if (SCM_BIGP (z
))
5441 else if (SCM_REALP (z
))
5443 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5444 SCM_OUT_OF_RANGE (1, z
);
5451 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5452 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5453 scm_i_mpz2num (mpq_denref (frac
)));
5455 /* When scm_i_make_ratio throws, we leak the memory allocated
5462 else if (SCM_FRACTIONP (z
))
5465 SCM_WRONG_TYPE_ARG (1, z
);
5469 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5471 "Return an exact number that is within @var{err} of @var{x}.")
5472 #define FUNC_NAME s_scm_rationalize
5474 if (SCM_I_INUMP (x
))
5476 else if (SCM_BIGP (x
))
5478 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5480 /* Use continued fractions to find closest ratio. All
5481 arithmetic is done with exact numbers.
5484 SCM ex
= scm_inexact_to_exact (x
);
5485 SCM int_part
= scm_floor (ex
);
5486 SCM tt
= SCM_I_MAKINUM (1);
5487 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
5488 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
5492 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
5495 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
5496 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
5498 /* We stop after a million iterations just to be absolutely sure
5499 that we don't go into an infinite loop. The process normally
5500 converges after less than a dozen iterations.
5503 err
= scm_abs (err
);
5504 while (++i
< 1000000)
5506 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
5507 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
5508 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
5510 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
5511 err
))) /* abs(x-a/b) <= err */
5513 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
5514 if (scm_is_false (scm_exact_p (x
))
5515 || scm_is_false (scm_exact_p (err
)))
5516 return scm_exact_to_inexact (res
);
5520 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
5522 tt
= scm_floor (rx
); /* tt = floor (rx) */
5528 scm_num_overflow (s_scm_rationalize
);
5531 SCM_WRONG_TYPE_ARG (1, x
);
5535 /* conversion functions */
5538 scm_is_integer (SCM val
)
5540 return scm_is_true (scm_integer_p (val
));
5544 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
5546 if (SCM_I_INUMP (val
))
5548 scm_t_signed_bits n
= SCM_I_INUM (val
);
5549 return n
>= min
&& n
<= max
;
5551 else if (SCM_BIGP (val
))
5553 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
5555 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
5557 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
5559 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
5560 return n
>= min
&& n
<= max
;
5570 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5571 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5574 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5575 SCM_I_BIG_MPZ (val
));
5577 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
5589 return n
>= min
&& n
<= max
;
5597 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
5599 if (SCM_I_INUMP (val
))
5601 scm_t_signed_bits n
= SCM_I_INUM (val
);
5602 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
5604 else if (SCM_BIGP (val
))
5606 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
5608 else if (max
<= ULONG_MAX
)
5610 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
5612 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
5613 return n
>= min
&& n
<= max
;
5623 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
5626 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
5627 > CHAR_BIT
*sizeof (scm_t_uintmax
))
5630 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
5631 SCM_I_BIG_MPZ (val
));
5633 return n
>= min
&& n
<= max
;
5641 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
5643 scm_error (scm_out_of_range_key
,
5645 "Value out of range ~S to ~S: ~S",
5646 scm_list_3 (min
, max
, bad_val
),
5647 scm_list_1 (bad_val
));
5650 #define TYPE scm_t_intmax
5651 #define TYPE_MIN min
5652 #define TYPE_MAX max
5653 #define SIZEOF_TYPE 0
5654 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
5655 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
5656 #include "libguile/conv-integer.i.c"
5658 #define TYPE scm_t_uintmax
5659 #define TYPE_MIN min
5660 #define TYPE_MAX max
5661 #define SIZEOF_TYPE 0
5662 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
5663 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
5664 #include "libguile/conv-uinteger.i.c"
5666 #define TYPE scm_t_int8
5667 #define TYPE_MIN SCM_T_INT8_MIN
5668 #define TYPE_MAX SCM_T_INT8_MAX
5669 #define SIZEOF_TYPE 1
5670 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
5671 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
5672 #include "libguile/conv-integer.i.c"
5674 #define TYPE scm_t_uint8
5676 #define TYPE_MAX SCM_T_UINT8_MAX
5677 #define SIZEOF_TYPE 1
5678 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
5679 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
5680 #include "libguile/conv-uinteger.i.c"
5682 #define TYPE scm_t_int16
5683 #define TYPE_MIN SCM_T_INT16_MIN
5684 #define TYPE_MAX SCM_T_INT16_MAX
5685 #define SIZEOF_TYPE 2
5686 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
5687 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
5688 #include "libguile/conv-integer.i.c"
5690 #define TYPE scm_t_uint16
5692 #define TYPE_MAX SCM_T_UINT16_MAX
5693 #define SIZEOF_TYPE 2
5694 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
5695 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
5696 #include "libguile/conv-uinteger.i.c"
5698 #define TYPE scm_t_int32
5699 #define TYPE_MIN SCM_T_INT32_MIN
5700 #define TYPE_MAX SCM_T_INT32_MAX
5701 #define SIZEOF_TYPE 4
5702 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
5703 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
5704 #include "libguile/conv-integer.i.c"
5706 #define TYPE scm_t_uint32
5708 #define TYPE_MAX SCM_T_UINT32_MAX
5709 #define SIZEOF_TYPE 4
5710 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
5711 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
5712 #include "libguile/conv-uinteger.i.c"
5714 #if SCM_HAVE_T_INT64
5716 #define TYPE scm_t_int64
5717 #define TYPE_MIN SCM_T_INT64_MIN
5718 #define TYPE_MAX SCM_T_INT64_MAX
5719 #define SIZEOF_TYPE 8
5720 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
5721 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
5722 #include "libguile/conv-integer.i.c"
5724 #define TYPE scm_t_uint64
5726 #define TYPE_MAX SCM_T_UINT64_MAX
5727 #define SIZEOF_TYPE 8
5728 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
5729 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
5730 #include "libguile/conv-uinteger.i.c"
5735 scm_to_mpz (SCM val
, mpz_t rop
)
5737 if (SCM_I_INUMP (val
))
5738 mpz_set_si (rop
, SCM_I_INUM (val
));
5739 else if (SCM_BIGP (val
))
5740 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
5742 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
5746 scm_from_mpz (mpz_t val
)
5748 return scm_i_mpz2num (val
);
5752 scm_is_real (SCM val
)
5754 return scm_is_true (scm_real_p (val
));
5758 scm_is_rational (SCM val
)
5760 return scm_is_true (scm_rational_p (val
));
5764 scm_to_double (SCM val
)
5766 if (SCM_I_INUMP (val
))
5767 return SCM_I_INUM (val
);
5768 else if (SCM_BIGP (val
))
5769 return scm_i_big2dbl (val
);
5770 else if (SCM_FRACTIONP (val
))
5771 return scm_i_fraction2double (val
);
5772 else if (SCM_REALP (val
))
5773 return SCM_REAL_VALUE (val
);
5775 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
5779 scm_from_double (double val
)
5781 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
5782 SCM_REAL_VALUE (z
) = val
;
5786 #if SCM_ENABLE_DISCOURAGED == 1
5789 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
5793 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5797 scm_out_of_range (NULL
, num
);
5800 return scm_to_double (num
);
5804 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
5808 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
5812 scm_out_of_range (NULL
, num
);
5815 return scm_to_double (num
);
5821 scm_is_complex (SCM val
)
5823 return scm_is_true (scm_complex_p (val
));
5827 scm_c_real_part (SCM z
)
5829 if (SCM_COMPLEXP (z
))
5830 return SCM_COMPLEX_REAL (z
);
5833 /* Use the scm_real_part to get proper error checking and
5836 return scm_to_double (scm_real_part (z
));
5841 scm_c_imag_part (SCM z
)
5843 if (SCM_COMPLEXP (z
))
5844 return SCM_COMPLEX_IMAG (z
);
5847 /* Use the scm_imag_part to get proper error checking and
5848 dispatching. The result will almost always be 0.0, but not
5851 return scm_to_double (scm_imag_part (z
));
5856 scm_c_magnitude (SCM z
)
5858 return scm_to_double (scm_magnitude (z
));
5864 return scm_to_double (scm_angle (z
));
5868 scm_is_number (SCM z
)
5870 return scm_is_true (scm_number_p (z
));
5878 mpz_init_set_si (z_negative_one
, -1);
5880 /* It may be possible to tune the performance of some algorithms by using
5881 * the following constants to avoid the creation of bignums. Please, before
5882 * using these values, remember the two rules of program optimization:
5883 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
5884 scm_c_define ("most-positive-fixnum",
5885 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
5886 scm_c_define ("most-negative-fixnum",
5887 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
5889 scm_add_feature ("complex");
5890 scm_add_feature ("inexact");
5891 scm_flo0
= scm_from_double (0.0);
5893 /* determine floating point precision */
5894 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
5896 init_dblprec(&scm_dblprec
[i
-2],i
);
5897 init_fx_radix(fx_per_radix
[i
-2],i
);
5900 /* hard code precision for base 10 if the preprocessor tells us to... */
5901 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
5904 exactly_one_half
= scm_permanent_object (scm_divide (SCM_I_MAKINUM (1),
5905 SCM_I_MAKINUM (2)));
5906 #include "libguile/numbers.x"