1 /* Copyright (C) 1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005, 2006, 2007, 2008, 2009, 2010 Free Software Foundation, Inc.
3 * Portions Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories
4 * and Bellcore. See scm_divide.
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public License
9 * as published by the Free Software Foundation; either version 3 of
10 * the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
24 /* General assumptions:
25 * All objects satisfying SCM_COMPLEXP() have a non-zero complex component.
26 * All objects satisfying SCM_BIGP() are too large to fit in a fixnum.
27 * If an object satisfies integer?, it's either an inum, a bignum, or a real.
28 * If floor (r) == r, r is an int, and mpz_set_d will DTRT.
29 * All objects satisfying SCM_FRACTIONP are never an integer.
34 - see if special casing bignums and reals in integer-exponent when
35 possible (to use mpz_pow and mpf_pow_ui) is faster.
37 - look in to better short-circuiting of common cases in
38 integer-expt and elsewhere.
40 - see if direct mpz operations can help in ash and elsewhere.
57 #include "libguile/_scm.h"
58 #include "libguile/feature.h"
59 #include "libguile/ports.h"
60 #include "libguile/root.h"
61 #include "libguile/smob.h"
62 #include "libguile/strings.h"
64 #include "libguile/validate.h"
65 #include "libguile/numbers.h"
66 #include "libguile/deprecation.h"
68 #include "libguile/eq.h"
70 #include "libguile/discouraged.h"
72 /* values per glibc, if not already defined */
74 #define M_LOG10E 0.43429448190325182765
77 #define M_PI 3.14159265358979323846
83 Wonder if this might be faster for some of our code? A switch on
84 the numtag would jump directly to the right case, and the
85 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
87 #define SCM_I_NUMTAG_NOTNUM 0
88 #define SCM_I_NUMTAG_INUM 1
89 #define SCM_I_NUMTAG_BIG scm_tc16_big
90 #define SCM_I_NUMTAG_REAL scm_tc16_real
91 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
92 #define SCM_I_NUMTAG(x) \
93 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
94 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
95 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
96 : SCM_I_NUMTAG_NOTNUM)))
98 /* the macro above will not work as is with fractions */
103 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
105 /* FLOBUFLEN is the maximum number of characters neccessary for the
106 * printed or scm_string representation of an inexact number.
108 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
111 #if ! defined (HAVE_ISNAN)
116 return (IsNANorINF (x
) && NaN (x
) && ! IsINF (x
)) ? 1 : 0;
119 #if ! defined (HAVE_ISINF)
124 return (IsNANorINF (x
) && IsINF (x
)) ? 1 : 0;
131 #if !defined (HAVE_ASINH)
132 static double asinh (double x
) { return log (x
+ sqrt (x
* x
+ 1)); }
134 #if !defined (HAVE_ACOSH)
135 static double acosh (double x
) { return log (x
+ sqrt (x
* x
- 1)); }
137 #if !defined (HAVE_ATANH)
138 static double atanh (double x
) { return 0.5 * log ((1 + x
) / (1 - x
)); }
141 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
142 an explicit check. In some future gmp (don't know what version number),
143 mpz_cmp_d is supposed to do this itself. */
145 #define xmpz_cmp_d(z, d) \
146 (xisinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
148 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
151 /* For reference, sparc solaris 7 has infinities (IEEE) but doesn't have
152 isinf. It does have finite and isnan though, hence the use of those.
153 fpclass would be a possibility on that system too. */
157 #if defined (HAVE_ISINF)
159 #elif defined (HAVE_FINITE) && defined (HAVE_ISNAN)
160 return (! (finite (x
) || isnan (x
)));
169 #if defined (HAVE_ISNAN)
176 #if defined (GUILE_I)
177 #if HAVE_COMPLEX_DOUBLE
179 /* For an SCM object Z which is a complex number (ie. satisfies
180 SCM_COMPLEXP), return its value as a C level "complex double". */
181 #define SCM_COMPLEX_VALUE(z) \
182 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
184 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
186 /* Convert a C "complex double" to an SCM value. */
188 scm_from_complex_double (complex double z
)
190 return scm_c_make_rectangular (creal (z
), cimag (z
));
193 #endif /* HAVE_COMPLEX_DOUBLE */
198 static mpz_t z_negative_one
;
205 /* Return a newly created bignum. */
206 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
207 mpz_init (SCM_I_BIG_MPZ (z
));
212 scm_i_long2big (long x
)
214 /* Return a newly created bignum initialized to X. */
215 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
216 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
221 scm_i_ulong2big (unsigned long x
)
223 /* Return a newly created bignum initialized to X. */
224 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
225 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
230 scm_i_clonebig (SCM src_big
, int same_sign_p
)
232 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
233 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
234 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
236 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
241 scm_i_bigcmp (SCM x
, SCM y
)
243 /* Return neg if x < y, pos if x > y, and 0 if x == y */
244 /* presume we already know x and y are bignums */
245 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
246 scm_remember_upto_here_2 (x
, y
);
251 scm_i_dbl2big (double d
)
253 /* results are only defined if d is an integer */
254 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
255 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
259 /* Convert a integer in double representation to a SCM number. */
262 scm_i_dbl2num (double u
)
264 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
265 powers of 2, so there's no rounding when making "double" values
266 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
267 get rounded on a 64-bit machine, hence the "+1".
269 The use of floor() to force to an integer value ensures we get a
270 "numerically closest" value without depending on how a
271 double->long cast or how mpz_set_d will round. For reference,
272 double->long probably follows the hardware rounding mode,
273 mpz_set_d truncates towards zero. */
275 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
276 representable as a double? */
278 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
279 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
280 return SCM_I_MAKINUM ((long) u
);
282 return scm_i_dbl2big (u
);
285 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
286 with R5RS exact->inexact.
288 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
289 (ie. truncate towards zero), then adjust to get the closest double by
290 examining the next lower bit and adding 1 (to the absolute value) if
293 Bignums exactly half way between representable doubles are rounded to the
294 next higher absolute value (ie. away from zero). This seems like an
295 adequate interpretation of R5RS "numerically closest", and it's easier
296 and faster than a full "nearest-even" style.
298 The bit test must be done on the absolute value of the mpz_t, which means
299 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
300 negatives as twos complement.
302 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
303 following the hardware rounding mode, but applied to the absolute value
304 of the mpz_t operand. This is not what we want so we put the high
305 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
306 mpz_get_d is supposed to always truncate towards zero.
308 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
309 is a slowdown. It'd be faster to pick out the relevant high bits with
310 mpz_getlimbn if we could be bothered coding that, and if the new
311 truncating gmp doesn't come out. */
314 scm_i_big2dbl (SCM b
)
319 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
323 /* Current GMP, eg. 4.1.3, force truncation towards zero */
325 if (bits
> DBL_MANT_DIG
)
327 size_t shift
= bits
- DBL_MANT_DIG
;
328 mpz_init2 (tmp
, DBL_MANT_DIG
);
329 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
330 result
= ldexp (mpz_get_d (tmp
), shift
);
335 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
340 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
343 if (bits
> DBL_MANT_DIG
)
345 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
346 /* test bit number "pos" in absolute value */
347 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
348 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
350 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
354 scm_remember_upto_here_1 (b
);
359 scm_i_normbig (SCM b
)
361 /* convert a big back to a fixnum if it'll fit */
362 /* presume b is a bignum */
363 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
365 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
366 if (SCM_FIXABLE (val
))
367 b
= SCM_I_MAKINUM (val
);
372 static SCM_C_INLINE_KEYWORD SCM
373 scm_i_mpz2num (mpz_t b
)
375 /* convert a mpz number to a SCM number. */
376 if (mpz_fits_slong_p (b
))
378 long val
= mpz_get_si (b
);
379 if (SCM_FIXABLE (val
))
380 return SCM_I_MAKINUM (val
);
384 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
385 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
390 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
391 static SCM
scm_divide2real (SCM x
, SCM y
);
394 scm_i_make_ratio (SCM numerator
, SCM denominator
)
395 #define FUNC_NAME "make-ratio"
397 /* First make sure the arguments are proper.
399 if (SCM_I_INUMP (denominator
))
401 if (scm_is_eq (denominator
, SCM_INUM0
))
402 scm_num_overflow ("make-ratio");
403 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
408 if (!(SCM_BIGP(denominator
)))
409 SCM_WRONG_TYPE_ARG (2, denominator
);
411 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
412 SCM_WRONG_TYPE_ARG (1, numerator
);
414 /* Then flip signs so that the denominator is positive.
416 if (scm_is_true (scm_negative_p (denominator
)))
418 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
419 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
422 /* Now consider for each of the four fixnum/bignum combinations
423 whether the rational number is really an integer.
425 if (SCM_I_INUMP (numerator
))
427 long x
= SCM_I_INUM (numerator
);
428 if (scm_is_eq (numerator
, SCM_INUM0
))
430 if (SCM_I_INUMP (denominator
))
433 y
= SCM_I_INUM (denominator
);
435 return SCM_I_MAKINUM(1);
437 return SCM_I_MAKINUM (x
/ y
);
441 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
442 of that value for the denominator, as a bignum. Apart from
443 that case, abs(bignum) > abs(inum) so inum/bignum is not an
445 if (x
== SCM_MOST_NEGATIVE_FIXNUM
446 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
447 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
448 return SCM_I_MAKINUM(-1);
451 else if (SCM_BIGP (numerator
))
453 if (SCM_I_INUMP (denominator
))
455 long yy
= SCM_I_INUM (denominator
);
456 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
457 return scm_divide (numerator
, denominator
);
461 if (scm_is_eq (numerator
, denominator
))
462 return SCM_I_MAKINUM(1);
463 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
464 SCM_I_BIG_MPZ (denominator
)))
465 return scm_divide(numerator
, denominator
);
469 /* No, it's a proper fraction.
472 SCM divisor
= scm_gcd (numerator
, denominator
);
473 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
475 numerator
= scm_divide (numerator
, divisor
);
476 denominator
= scm_divide (denominator
, divisor
);
479 return scm_double_cell (scm_tc16_fraction
,
480 SCM_UNPACK (numerator
),
481 SCM_UNPACK (denominator
), 0);
487 scm_i_fraction2double (SCM z
)
489 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
490 SCM_FRACTION_DENOMINATOR (z
)));
493 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
495 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
497 #define FUNC_NAME s_scm_exact_p
503 if (SCM_FRACTIONP (x
))
507 SCM_WRONG_TYPE_ARG (1, x
);
512 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
514 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
516 #define FUNC_NAME s_scm_odd_p
520 long val
= SCM_I_INUM (n
);
521 return scm_from_bool ((val
& 1L) != 0);
523 else if (SCM_BIGP (n
))
525 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
526 scm_remember_upto_here_1 (n
);
527 return scm_from_bool (odd_p
);
529 else if (scm_is_true (scm_inf_p (n
)))
531 else if (SCM_REALP (n
))
533 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
539 SCM_WRONG_TYPE_ARG (1, n
);
542 SCM_WRONG_TYPE_ARG (1, n
);
547 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
549 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
551 #define FUNC_NAME s_scm_even_p
555 long val
= SCM_I_INUM (n
);
556 return scm_from_bool ((val
& 1L) == 0);
558 else if (SCM_BIGP (n
))
560 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
561 scm_remember_upto_here_1 (n
);
562 return scm_from_bool (even_p
);
564 else if (scm_is_true (scm_inf_p (n
)))
566 else if (SCM_REALP (n
))
568 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
574 SCM_WRONG_TYPE_ARG (1, n
);
577 SCM_WRONG_TYPE_ARG (1, n
);
581 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
583 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
584 "or @samp{-inf.0}, @code{#f} otherwise.")
585 #define FUNC_NAME s_scm_inf_p
588 return scm_from_bool (xisinf (SCM_REAL_VALUE (x
)));
589 else if (SCM_COMPLEXP (x
))
590 return scm_from_bool (xisinf (SCM_COMPLEX_REAL (x
))
591 || xisinf (SCM_COMPLEX_IMAG (x
)));
597 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
599 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
601 #define FUNC_NAME s_scm_nan_p
604 return scm_from_bool (xisnan (SCM_REAL_VALUE (n
)));
605 else if (SCM_COMPLEXP (n
))
606 return scm_from_bool (xisnan (SCM_COMPLEX_REAL (n
))
607 || xisnan (SCM_COMPLEX_IMAG (n
)));
613 /* Guile's idea of infinity. */
614 static double guile_Inf
;
616 /* Guile's idea of not a number. */
617 static double guile_NaN
;
620 guile_ieee_init (void)
622 #if defined (HAVE_ISINF) || defined (HAVE_FINITE)
624 /* Some version of gcc on some old version of Linux used to crash when
625 trying to make Inf and NaN. */
628 /* C99 INFINITY, when available.
629 FIXME: The standard allows for INFINITY to be something that overflows
630 at compile time. We ought to have a configure test to check for that
631 before trying to use it. (But in practice we believe this is not a
632 problem on any system guile is likely to target.) */
633 guile_Inf
= INFINITY
;
634 #elif defined HAVE_DINFINITY
636 extern unsigned int DINFINITY
[2];
637 guile_Inf
= (*((double *) (DINFINITY
)));
644 if (guile_Inf
== tmp
)
652 #if defined (HAVE_ISNAN)
655 /* C99 NAN, when available */
657 #elif defined HAVE_DQNAN
660 extern unsigned int DQNAN
[2];
661 guile_NaN
= (*((double *)(DQNAN
)));
664 guile_NaN
= guile_Inf
/ guile_Inf
;
670 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
673 #define FUNC_NAME s_scm_inf
675 static int initialized
= 0;
681 return scm_from_double (guile_Inf
);
685 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
688 #define FUNC_NAME s_scm_nan
690 static int initialized
= 0;
696 return scm_from_double (guile_NaN
);
701 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
703 "Return the absolute value of @var{x}.")
708 long int xx
= SCM_I_INUM (x
);
711 else if (SCM_POSFIXABLE (-xx
))
712 return SCM_I_MAKINUM (-xx
);
714 return scm_i_long2big (-xx
);
716 else if (SCM_BIGP (x
))
718 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
720 return scm_i_clonebig (x
, 0);
724 else if (SCM_REALP (x
))
726 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
727 double xx
= SCM_REAL_VALUE (x
);
729 return scm_from_double (-xx
);
733 else if (SCM_FRACTIONP (x
))
735 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
737 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
738 SCM_FRACTION_DENOMINATOR (x
));
741 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
746 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
747 /* "Return the quotient of the numbers @var{x} and @var{y}."
750 scm_quotient (SCM x
, SCM y
)
754 long xx
= SCM_I_INUM (x
);
757 long yy
= SCM_I_INUM (y
);
759 scm_num_overflow (s_quotient
);
764 return SCM_I_MAKINUM (z
);
766 return scm_i_long2big (z
);
769 else if (SCM_BIGP (y
))
771 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
772 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
773 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
775 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
776 scm_remember_upto_here_1 (y
);
777 return SCM_I_MAKINUM (-1);
780 return SCM_I_MAKINUM (0);
783 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
785 else if (SCM_BIGP (x
))
789 long yy
= SCM_I_INUM (y
);
791 scm_num_overflow (s_quotient
);
796 SCM result
= scm_i_mkbig ();
799 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
802 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
805 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
806 scm_remember_upto_here_1 (x
);
807 return scm_i_normbig (result
);
810 else if (SCM_BIGP (y
))
812 SCM result
= scm_i_mkbig ();
813 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
816 scm_remember_upto_here_2 (x
, y
);
817 return scm_i_normbig (result
);
820 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
823 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
826 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
827 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
829 * "(remainder 13 4) @result{} 1\n"
830 * "(remainder -13 4) @result{} -1\n"
834 scm_remainder (SCM x
, SCM y
)
840 long yy
= SCM_I_INUM (y
);
842 scm_num_overflow (s_remainder
);
845 long z
= SCM_I_INUM (x
) % yy
;
846 return SCM_I_MAKINUM (z
);
849 else if (SCM_BIGP (y
))
851 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
852 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
853 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
855 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
856 scm_remember_upto_here_1 (y
);
857 return SCM_I_MAKINUM (0);
863 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
865 else if (SCM_BIGP (x
))
869 long yy
= SCM_I_INUM (y
);
871 scm_num_overflow (s_remainder
);
874 SCM result
= scm_i_mkbig ();
877 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
878 scm_remember_upto_here_1 (x
);
879 return scm_i_normbig (result
);
882 else if (SCM_BIGP (y
))
884 SCM result
= scm_i_mkbig ();
885 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
888 scm_remember_upto_here_2 (x
, y
);
889 return scm_i_normbig (result
);
892 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
895 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
899 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
900 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
902 * "(modulo 13 4) @result{} 1\n"
903 * "(modulo -13 4) @result{} 3\n"
907 scm_modulo (SCM x
, SCM y
)
911 long xx
= SCM_I_INUM (x
);
914 long yy
= SCM_I_INUM (y
);
916 scm_num_overflow (s_modulo
);
919 /* C99 specifies that "%" is the remainder corresponding to a
920 quotient rounded towards zero, and that's also traditional
921 for machine division, so z here should be well defined. */
939 return SCM_I_MAKINUM (result
);
942 else if (SCM_BIGP (y
))
944 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
951 SCM pos_y
= scm_i_clonebig (y
, 0);
952 /* do this after the last scm_op */
953 mpz_init_set_si (z_x
, xx
);
954 result
= pos_y
; /* re-use this bignum */
955 mpz_mod (SCM_I_BIG_MPZ (result
),
957 SCM_I_BIG_MPZ (pos_y
));
958 scm_remember_upto_here_1 (pos_y
);
962 result
= scm_i_mkbig ();
963 /* do this after the last scm_op */
964 mpz_init_set_si (z_x
, xx
);
965 mpz_mod (SCM_I_BIG_MPZ (result
),
968 scm_remember_upto_here_1 (y
);
971 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
972 mpz_add (SCM_I_BIG_MPZ (result
),
974 SCM_I_BIG_MPZ (result
));
975 scm_remember_upto_here_1 (y
);
976 /* and do this before the next one */
978 return scm_i_normbig (result
);
982 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
984 else if (SCM_BIGP (x
))
988 long yy
= SCM_I_INUM (y
);
990 scm_num_overflow (s_modulo
);
993 SCM result
= scm_i_mkbig ();
994 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
996 (yy
< 0) ? - yy
: yy
);
997 scm_remember_upto_here_1 (x
);
998 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
999 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
1000 SCM_I_BIG_MPZ (result
),
1002 return scm_i_normbig (result
);
1005 else if (SCM_BIGP (y
))
1008 SCM result
= scm_i_mkbig ();
1009 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
1010 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
1011 mpz_mod (SCM_I_BIG_MPZ (result
),
1013 SCM_I_BIG_MPZ (pos_y
));
1015 scm_remember_upto_here_1 (x
);
1016 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
1017 mpz_add (SCM_I_BIG_MPZ (result
),
1019 SCM_I_BIG_MPZ (result
));
1020 scm_remember_upto_here_2 (y
, pos_y
);
1021 return scm_i_normbig (result
);
1025 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
1028 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
1031 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
1032 (SCM x
, SCM y
, SCM rest
),
1033 "Return the greatest common divisor of all parameter values.\n"
1034 "If called without arguments, 0 is returned.")
1035 #define FUNC_NAME s_scm_i_gcd
1037 while (!scm_is_null (rest
))
1038 { x
= scm_gcd (x
, y
);
1040 rest
= scm_cdr (rest
);
1042 return scm_gcd (x
, y
);
1046 #define s_gcd s_scm_i_gcd
1047 #define g_gcd g_scm_i_gcd
1050 scm_gcd (SCM x
, SCM y
)
1053 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
1055 if (SCM_I_INUMP (x
))
1057 if (SCM_I_INUMP (y
))
1059 long xx
= SCM_I_INUM (x
);
1060 long yy
= SCM_I_INUM (y
);
1061 long u
= xx
< 0 ? -xx
: xx
;
1062 long v
= yy
< 0 ? -yy
: yy
;
1072 /* Determine a common factor 2^k */
1073 while (!(1 & (u
| v
)))
1079 /* Now, any factor 2^n can be eliminated */
1099 return (SCM_POSFIXABLE (result
)
1100 ? SCM_I_MAKINUM (result
)
1101 : scm_i_long2big (result
));
1103 else if (SCM_BIGP (y
))
1109 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1111 else if (SCM_BIGP (x
))
1113 if (SCM_I_INUMP (y
))
1115 unsigned long result
;
1118 yy
= SCM_I_INUM (y
);
1123 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1124 scm_remember_upto_here_1 (x
);
1125 return (SCM_POSFIXABLE (result
)
1126 ? SCM_I_MAKINUM (result
)
1127 : scm_from_ulong (result
));
1129 else if (SCM_BIGP (y
))
1131 SCM result
= scm_i_mkbig ();
1132 mpz_gcd (SCM_I_BIG_MPZ (result
),
1135 scm_remember_upto_here_2 (x
, y
);
1136 return scm_i_normbig (result
);
1139 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1142 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1145 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
1146 (SCM x
, SCM y
, SCM rest
),
1147 "Return the least common multiple of the arguments.\n"
1148 "If called without arguments, 1 is returned.")
1149 #define FUNC_NAME s_scm_i_lcm
1151 while (!scm_is_null (rest
))
1152 { x
= scm_lcm (x
, y
);
1154 rest
= scm_cdr (rest
);
1156 return scm_lcm (x
, y
);
1160 #define s_lcm s_scm_i_lcm
1161 #define g_lcm g_scm_i_lcm
1164 scm_lcm (SCM n1
, SCM n2
)
1166 if (SCM_UNBNDP (n2
))
1168 if (SCM_UNBNDP (n1
))
1169 return SCM_I_MAKINUM (1L);
1170 n2
= SCM_I_MAKINUM (1L);
1173 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1174 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1175 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1176 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1178 if (SCM_I_INUMP (n1
))
1180 if (SCM_I_INUMP (n2
))
1182 SCM d
= scm_gcd (n1
, n2
);
1183 if (scm_is_eq (d
, SCM_INUM0
))
1186 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1190 /* inum n1, big n2 */
1193 SCM result
= scm_i_mkbig ();
1194 long nn1
= SCM_I_INUM (n1
);
1195 if (nn1
== 0) return SCM_INUM0
;
1196 if (nn1
< 0) nn1
= - nn1
;
1197 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1198 scm_remember_upto_here_1 (n2
);
1206 if (SCM_I_INUMP (n2
))
1213 SCM result
= scm_i_mkbig ();
1214 mpz_lcm(SCM_I_BIG_MPZ (result
),
1216 SCM_I_BIG_MPZ (n2
));
1217 scm_remember_upto_here_2(n1
, n2
);
1218 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1224 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1229 + + + x (map digit:logand X Y)
1230 + - + x (map digit:logand X (lognot (+ -1 Y)))
1231 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1232 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1237 + + + (map digit:logior X Y)
1238 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1239 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1240 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1245 + + + (map digit:logxor X Y)
1246 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1247 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1248 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1253 + + (any digit:logand X Y)
1254 + - (any digit:logand X (lognot (+ -1 Y)))
1255 - + (any digit:logand (lognot (+ -1 X)) Y)
1260 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
1261 (SCM x
, SCM y
, SCM rest
),
1262 "Return the bitwise AND of the integer arguments.\n\n"
1264 "(logand) @result{} -1\n"
1265 "(logand 7) @result{} 7\n"
1266 "(logand #b111 #b011 #b001) @result{} 1\n"
1268 #define FUNC_NAME s_scm_i_logand
1270 while (!scm_is_null (rest
))
1271 { x
= scm_logand (x
, y
);
1273 rest
= scm_cdr (rest
);
1275 return scm_logand (x
, y
);
1279 #define s_scm_logand s_scm_i_logand
1281 SCM
scm_logand (SCM n1
, SCM n2
)
1282 #define FUNC_NAME s_scm_logand
1286 if (SCM_UNBNDP (n2
))
1288 if (SCM_UNBNDP (n1
))
1289 return SCM_I_MAKINUM (-1);
1290 else if (!SCM_NUMBERP (n1
))
1291 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1292 else if (SCM_NUMBERP (n1
))
1295 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1298 if (SCM_I_INUMP (n1
))
1300 nn1
= SCM_I_INUM (n1
);
1301 if (SCM_I_INUMP (n2
))
1303 long nn2
= SCM_I_INUM (n2
);
1304 return SCM_I_MAKINUM (nn1
& nn2
);
1306 else if SCM_BIGP (n2
)
1312 SCM result_z
= scm_i_mkbig ();
1314 mpz_init_set_si (nn1_z
, nn1
);
1315 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1316 scm_remember_upto_here_1 (n2
);
1318 return scm_i_normbig (result_z
);
1322 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1324 else if (SCM_BIGP (n1
))
1326 if (SCM_I_INUMP (n2
))
1329 nn1
= SCM_I_INUM (n1
);
1332 else if (SCM_BIGP (n2
))
1334 SCM result_z
= scm_i_mkbig ();
1335 mpz_and (SCM_I_BIG_MPZ (result_z
),
1337 SCM_I_BIG_MPZ (n2
));
1338 scm_remember_upto_here_2 (n1
, n2
);
1339 return scm_i_normbig (result_z
);
1342 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1345 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1350 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
1351 (SCM x
, SCM y
, SCM rest
),
1352 "Return the bitwise OR of the integer arguments.\n\n"
1354 "(logior) @result{} 0\n"
1355 "(logior 7) @result{} 7\n"
1356 "(logior #b000 #b001 #b011) @result{} 3\n"
1358 #define FUNC_NAME s_scm_i_logior
1360 while (!scm_is_null (rest
))
1361 { x
= scm_logior (x
, y
);
1363 rest
= scm_cdr (rest
);
1365 return scm_logior (x
, y
);
1369 #define s_scm_logior s_scm_i_logior
1371 SCM
scm_logior (SCM n1
, SCM n2
)
1372 #define FUNC_NAME s_scm_logior
1376 if (SCM_UNBNDP (n2
))
1378 if (SCM_UNBNDP (n1
))
1380 else if (SCM_NUMBERP (n1
))
1383 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1386 if (SCM_I_INUMP (n1
))
1388 nn1
= SCM_I_INUM (n1
);
1389 if (SCM_I_INUMP (n2
))
1391 long nn2
= SCM_I_INUM (n2
);
1392 return SCM_I_MAKINUM (nn1
| nn2
);
1394 else if (SCM_BIGP (n2
))
1400 SCM result_z
= scm_i_mkbig ();
1402 mpz_init_set_si (nn1_z
, nn1
);
1403 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1404 scm_remember_upto_here_1 (n2
);
1406 return scm_i_normbig (result_z
);
1410 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1412 else if (SCM_BIGP (n1
))
1414 if (SCM_I_INUMP (n2
))
1417 nn1
= SCM_I_INUM (n1
);
1420 else if (SCM_BIGP (n2
))
1422 SCM result_z
= scm_i_mkbig ();
1423 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1425 SCM_I_BIG_MPZ (n2
));
1426 scm_remember_upto_here_2 (n1
, n2
);
1427 return scm_i_normbig (result_z
);
1430 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1433 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1438 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
1439 (SCM x
, SCM y
, SCM rest
),
1440 "Return the bitwise XOR of the integer arguments. A bit is\n"
1441 "set in the result if it is set in an odd number of arguments.\n"
1443 "(logxor) @result{} 0\n"
1444 "(logxor 7) @result{} 7\n"
1445 "(logxor #b000 #b001 #b011) @result{} 2\n"
1446 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1448 #define FUNC_NAME s_scm_i_logxor
1450 while (!scm_is_null (rest
))
1451 { x
= scm_logxor (x
, y
);
1453 rest
= scm_cdr (rest
);
1455 return scm_logxor (x
, y
);
1459 #define s_scm_logxor s_scm_i_logxor
1461 SCM
scm_logxor (SCM n1
, SCM n2
)
1462 #define FUNC_NAME s_scm_logxor
1466 if (SCM_UNBNDP (n2
))
1468 if (SCM_UNBNDP (n1
))
1470 else if (SCM_NUMBERP (n1
))
1473 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1476 if (SCM_I_INUMP (n1
))
1478 nn1
= SCM_I_INUM (n1
);
1479 if (SCM_I_INUMP (n2
))
1481 long nn2
= SCM_I_INUM (n2
);
1482 return SCM_I_MAKINUM (nn1
^ nn2
);
1484 else if (SCM_BIGP (n2
))
1488 SCM result_z
= scm_i_mkbig ();
1490 mpz_init_set_si (nn1_z
, nn1
);
1491 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1492 scm_remember_upto_here_1 (n2
);
1494 return scm_i_normbig (result_z
);
1498 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1500 else if (SCM_BIGP (n1
))
1502 if (SCM_I_INUMP (n2
))
1505 nn1
= SCM_I_INUM (n1
);
1508 else if (SCM_BIGP (n2
))
1510 SCM result_z
= scm_i_mkbig ();
1511 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1513 SCM_I_BIG_MPZ (n2
));
1514 scm_remember_upto_here_2 (n1
, n2
);
1515 return scm_i_normbig (result_z
);
1518 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1521 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1526 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1528 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1529 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1530 "without actually calculating the @code{logand}, just testing\n"
1534 "(logtest #b0100 #b1011) @result{} #f\n"
1535 "(logtest #b0100 #b0111) @result{} #t\n"
1537 #define FUNC_NAME s_scm_logtest
1541 if (SCM_I_INUMP (j
))
1543 nj
= SCM_I_INUM (j
);
1544 if (SCM_I_INUMP (k
))
1546 long nk
= SCM_I_INUM (k
);
1547 return scm_from_bool (nj
& nk
);
1549 else if (SCM_BIGP (k
))
1557 mpz_init_set_si (nj_z
, nj
);
1558 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1559 scm_remember_upto_here_1 (k
);
1560 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1566 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1568 else if (SCM_BIGP (j
))
1570 if (SCM_I_INUMP (k
))
1573 nj
= SCM_I_INUM (j
);
1576 else if (SCM_BIGP (k
))
1580 mpz_init (result_z
);
1584 scm_remember_upto_here_2 (j
, k
);
1585 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1586 mpz_clear (result_z
);
1590 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1593 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1598 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1600 "Test whether bit number @var{index} in @var{j} is set.\n"
1601 "@var{index} starts from 0 for the least significant bit.\n"
1604 "(logbit? 0 #b1101) @result{} #t\n"
1605 "(logbit? 1 #b1101) @result{} #f\n"
1606 "(logbit? 2 #b1101) @result{} #t\n"
1607 "(logbit? 3 #b1101) @result{} #t\n"
1608 "(logbit? 4 #b1101) @result{} #f\n"
1610 #define FUNC_NAME s_scm_logbit_p
1612 unsigned long int iindex
;
1613 iindex
= scm_to_ulong (index
);
1615 if (SCM_I_INUMP (j
))
1617 /* bits above what's in an inum follow the sign bit */
1618 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1619 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1621 else if (SCM_BIGP (j
))
1623 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1624 scm_remember_upto_here_1 (j
);
1625 return scm_from_bool (val
);
1628 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1633 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1635 "Return the integer which is the ones-complement of the integer\n"
1639 "(number->string (lognot #b10000000) 2)\n"
1640 " @result{} \"-10000001\"\n"
1641 "(number->string (lognot #b0) 2)\n"
1642 " @result{} \"-1\"\n"
1644 #define FUNC_NAME s_scm_lognot
1646 if (SCM_I_INUMP (n
)) {
1647 /* No overflow here, just need to toggle all the bits making up the inum.
1648 Enhancement: No need to strip the tag and add it back, could just xor
1649 a block of 1 bits, if that worked with the various debug versions of
1651 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1653 } else if (SCM_BIGP (n
)) {
1654 SCM result
= scm_i_mkbig ();
1655 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1656 scm_remember_upto_here_1 (n
);
1660 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1665 /* returns 0 if IN is not an integer. OUT must already be
1668 coerce_to_big (SCM in
, mpz_t out
)
1671 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1672 else if (SCM_I_INUMP (in
))
1673 mpz_set_si (out
, SCM_I_INUM (in
));
1680 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1681 (SCM n
, SCM k
, SCM m
),
1682 "Return @var{n} raised to the integer exponent\n"
1683 "@var{k}, modulo @var{m}.\n"
1686 "(modulo-expt 2 3 5)\n"
1689 #define FUNC_NAME s_scm_modulo_expt
1695 /* There are two classes of error we might encounter --
1696 1) Math errors, which we'll report by calling scm_num_overflow,
1698 2) wrong-type errors, which of course we'll report by calling
1700 We don't report those errors immediately, however; instead we do
1701 some cleanup first. These variables tell us which error (if
1702 any) we should report after cleaning up.
1704 int report_overflow
= 0;
1706 int position_of_wrong_type
= 0;
1707 SCM value_of_wrong_type
= SCM_INUM0
;
1709 SCM result
= SCM_UNDEFINED
;
1715 if (scm_is_eq (m
, SCM_INUM0
))
1717 report_overflow
= 1;
1721 if (!coerce_to_big (n
, n_tmp
))
1723 value_of_wrong_type
= n
;
1724 position_of_wrong_type
= 1;
1728 if (!coerce_to_big (k
, k_tmp
))
1730 value_of_wrong_type
= k
;
1731 position_of_wrong_type
= 2;
1735 if (!coerce_to_big (m
, m_tmp
))
1737 value_of_wrong_type
= m
;
1738 position_of_wrong_type
= 3;
1742 /* if the exponent K is negative, and we simply call mpz_powm, we
1743 will get a divide-by-zero exception when an inverse 1/n mod m
1744 doesn't exist (or is not unique). Since exceptions are hard to
1745 handle, we'll attempt the inversion "by hand" -- that way, we get
1746 a simple failure code, which is easy to handle. */
1748 if (-1 == mpz_sgn (k_tmp
))
1750 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1752 report_overflow
= 1;
1755 mpz_neg (k_tmp
, k_tmp
);
1758 result
= scm_i_mkbig ();
1759 mpz_powm (SCM_I_BIG_MPZ (result
),
1764 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1765 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1772 if (report_overflow
)
1773 scm_num_overflow (FUNC_NAME
);
1775 if (position_of_wrong_type
)
1776 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1777 value_of_wrong_type
);
1779 return scm_i_normbig (result
);
1783 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1785 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1786 "exact integer, @var{n} can be any number.\n"
1788 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1789 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1790 "includes @math{0^0} is 1.\n"
1793 "(integer-expt 2 5) @result{} 32\n"
1794 "(integer-expt -3 3) @result{} -27\n"
1795 "(integer-expt 5 -3) @result{} 1/125\n"
1796 "(integer-expt 0 0) @result{} 1\n"
1798 #define FUNC_NAME s_scm_integer_expt
1801 SCM z_i2
= SCM_BOOL_F
;
1803 SCM acc
= SCM_I_MAKINUM (1L);
1805 /* 0^0 == 1 according to R5RS */
1806 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1807 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1808 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1809 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1811 if (SCM_I_INUMP (k
))
1812 i2
= SCM_I_INUM (k
);
1813 else if (SCM_BIGP (k
))
1815 z_i2
= scm_i_clonebig (k
, 1);
1816 scm_remember_upto_here_1 (k
);
1820 SCM_WRONG_TYPE_ARG (2, k
);
1824 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1826 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1827 n
= scm_divide (n
, SCM_UNDEFINED
);
1831 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1835 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1837 return scm_product (acc
, n
);
1839 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1840 acc
= scm_product (acc
, n
);
1841 n
= scm_product (n
, n
);
1842 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1850 n
= scm_divide (n
, SCM_UNDEFINED
);
1857 return scm_product (acc
, n
);
1859 acc
= scm_product (acc
, n
);
1860 n
= scm_product (n
, n
);
1867 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1869 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1870 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1872 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1873 "@var{cnt} is negative it's a division, rounded towards negative\n"
1874 "infinity. (Note that this is not the same rounding as\n"
1875 "@code{quotient} does.)\n"
1877 "With @var{n} viewed as an infinite precision twos complement,\n"
1878 "@code{ash} means a left shift introducing zero bits, or a right\n"
1879 "shift dropping bits.\n"
1882 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1883 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1885 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1886 "(ash -23 -2) @result{} -6\n"
1888 #define FUNC_NAME s_scm_ash
1891 bits_to_shift
= scm_to_long (cnt
);
1893 if (SCM_I_INUMP (n
))
1895 long nn
= SCM_I_INUM (n
);
1897 if (bits_to_shift
> 0)
1899 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1900 overflow a non-zero fixnum. For smaller shifts we check the
1901 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1902 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1903 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1909 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1911 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1914 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1918 SCM result
= scm_i_long2big (nn
);
1919 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1926 bits_to_shift
= -bits_to_shift
;
1927 if (bits_to_shift
>= SCM_LONG_BIT
)
1928 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1930 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1934 else if (SCM_BIGP (n
))
1938 if (bits_to_shift
== 0)
1941 result
= scm_i_mkbig ();
1942 if (bits_to_shift
>= 0)
1944 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1950 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1951 we have to allocate a bignum even if the result is going to be a
1953 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1955 return scm_i_normbig (result
);
1961 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1967 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1968 (SCM n
, SCM start
, SCM end
),
1969 "Return the integer composed of the @var{start} (inclusive)\n"
1970 "through @var{end} (exclusive) bits of @var{n}. The\n"
1971 "@var{start}th bit becomes the 0-th bit in the result.\n"
1974 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1975 " @result{} \"1010\"\n"
1976 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1977 " @result{} \"10110\"\n"
1979 #define FUNC_NAME s_scm_bit_extract
1981 unsigned long int istart
, iend
, bits
;
1982 istart
= scm_to_ulong (start
);
1983 iend
= scm_to_ulong (end
);
1984 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1986 /* how many bits to keep */
1987 bits
= iend
- istart
;
1989 if (SCM_I_INUMP (n
))
1991 long int in
= SCM_I_INUM (n
);
1993 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1994 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1995 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1997 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1999 /* Since we emulate two's complement encoded numbers, this
2000 * special case requires us to produce a result that has
2001 * more bits than can be stored in a fixnum.
2003 SCM result
= scm_i_long2big (in
);
2004 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
2009 /* mask down to requisite bits */
2010 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
2011 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
2013 else if (SCM_BIGP (n
))
2018 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
2022 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
2023 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
2024 such bits into a ulong. */
2025 result
= scm_i_mkbig ();
2026 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
2027 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
2028 result
= scm_i_normbig (result
);
2030 scm_remember_upto_here_1 (n
);
2034 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2039 static const char scm_logtab
[] = {
2040 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
2043 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
2045 "Return the number of bits in integer @var{n}. If integer is\n"
2046 "positive, the 1-bits in its binary representation are counted.\n"
2047 "If negative, the 0-bits in its two's-complement binary\n"
2048 "representation are counted. If 0, 0 is returned.\n"
2051 "(logcount #b10101010)\n"
2058 #define FUNC_NAME s_scm_logcount
2060 if (SCM_I_INUMP (n
))
2062 unsigned long int c
= 0;
2063 long int nn
= SCM_I_INUM (n
);
2068 c
+= scm_logtab
[15 & nn
];
2071 return SCM_I_MAKINUM (c
);
2073 else if (SCM_BIGP (n
))
2075 unsigned long count
;
2076 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
2077 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
2079 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
2080 scm_remember_upto_here_1 (n
);
2081 return SCM_I_MAKINUM (count
);
2084 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2089 static const char scm_ilentab
[] = {
2090 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2094 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2096 "Return the number of bits necessary to represent @var{n}.\n"
2099 "(integer-length #b10101010)\n"
2101 "(integer-length 0)\n"
2103 "(integer-length #b1111)\n"
2106 #define FUNC_NAME s_scm_integer_length
2108 if (SCM_I_INUMP (n
))
2110 unsigned long int c
= 0;
2112 long int nn
= SCM_I_INUM (n
);
2118 l
= scm_ilentab
[15 & nn
];
2121 return SCM_I_MAKINUM (c
- 4 + l
);
2123 else if (SCM_BIGP (n
))
2125 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2126 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2127 1 too big, so check for that and adjust. */
2128 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2129 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2130 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2131 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2133 scm_remember_upto_here_1 (n
);
2134 return SCM_I_MAKINUM (size
);
2137 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2141 /*** NUMBERS -> STRINGS ***/
2142 #define SCM_MAX_DBL_PREC 60
2143 #define SCM_MAX_DBL_RADIX 36
2145 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2146 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2147 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2150 void init_dblprec(int *prec
, int radix
) {
2151 /* determine floating point precision by adding successively
2152 smaller increments to 1.0 until it is considered == 1.0 */
2153 double f
= ((double)1.0)/radix
;
2154 double fsum
= 1.0 + f
;
2159 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2171 void init_fx_radix(double *fx_list
, int radix
)
2173 /* initialize a per-radix list of tolerances. When added
2174 to a number < 1.0, we can determine if we should raund
2175 up and quit converting a number to a string. */
2179 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2180 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2183 /* use this array as a way to generate a single digit */
2184 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2187 idbl2str (double f
, char *a
, int radix
)
2189 int efmt
, dpt
, d
, i
, wp
;
2191 #ifdef DBL_MIN_10_EXP
2194 #endif /* DBL_MIN_10_EXP */
2199 radix
> SCM_MAX_DBL_RADIX
)
2201 /* revert to existing behavior */
2205 wp
= scm_dblprec
[radix
-2];
2206 fx
= fx_per_radix
[radix
-2];
2210 #ifdef HAVE_COPYSIGN
2211 double sgn
= copysign (1.0, f
);
2216 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2222 strcpy (a
, "-inf.0");
2224 strcpy (a
, "+inf.0");
2227 else if (xisnan (f
))
2229 strcpy (a
, "+nan.0");
2239 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2240 make-uniform-vector, from causing infinite loops. */
2241 /* just do the checking...if it passes, we do the conversion for our
2242 radix again below */
2249 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2257 while (f_cpy
> 10.0)
2260 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2281 if (f
+ fx
[wp
] >= radix
)
2288 /* adding 9999 makes this equivalent to abs(x) % 3 */
2289 dpt
= (exp
+ 9999) % 3;
2293 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2315 a
[ch
++] = number_chars
[d
];
2318 if (f
+ fx
[wp
] >= 1.0)
2320 a
[ch
- 1] = number_chars
[d
+1];
2332 if ((dpt
> 4) && (exp
> 6))
2334 d
= (a
[0] == '-' ? 2 : 1);
2335 for (i
= ch
++; i
> d
; i
--)
2348 if (a
[ch
- 1] == '.')
2349 a
[ch
++] = '0'; /* trailing zero */
2358 for (i
= radix
; i
<= exp
; i
*= radix
);
2359 for (i
/= radix
; i
; i
/= radix
)
2361 a
[ch
++] = number_chars
[exp
/ i
];
2370 icmplx2str (double real
, double imag
, char *str
, int radix
)
2374 i
= idbl2str (real
, str
, radix
);
2377 /* Don't output a '+' for negative numbers or for Inf and
2378 NaN. They will provide their own sign. */
2379 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2381 i
+= idbl2str (imag
, &str
[i
], radix
);
2388 iflo2str (SCM flt
, char *str
, int radix
)
2391 if (SCM_REALP (flt
))
2392 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2394 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2399 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2400 characters in the result.
2402 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2404 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2409 return scm_iuint2str (-num
, rad
, p
) + 1;
2412 return scm_iuint2str (num
, rad
, p
);
2415 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2416 characters in the result.
2418 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2420 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2424 scm_t_uintmax n
= num
;
2426 for (n
/= rad
; n
> 0; n
/= rad
)
2436 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2441 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2443 "Return a string holding the external representation of the\n"
2444 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2445 "inexact, a radix of 10 will be used.")
2446 #define FUNC_NAME s_scm_number_to_string
2450 if (SCM_UNBNDP (radix
))
2453 base
= scm_to_signed_integer (radix
, 2, 36);
2455 if (SCM_I_INUMP (n
))
2457 char num_buf
[SCM_INTBUFLEN
];
2458 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2459 return scm_from_locale_stringn (num_buf
, length
);
2461 else if (SCM_BIGP (n
))
2463 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2464 scm_remember_upto_here_1 (n
);
2465 return scm_take_locale_string (str
);
2467 else if (SCM_FRACTIONP (n
))
2469 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2470 scm_from_locale_string ("/"),
2471 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2473 else if (SCM_INEXACTP (n
))
2475 char num_buf
[FLOBUFLEN
];
2476 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2479 SCM_WRONG_TYPE_ARG (1, n
);
2484 /* These print routines used to be stubbed here so that scm_repl.c
2485 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2488 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2490 char num_buf
[FLOBUFLEN
];
2491 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2496 scm_i_print_double (double val
, SCM port
)
2498 char num_buf
[FLOBUFLEN
];
2499 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2503 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2506 char num_buf
[FLOBUFLEN
];
2507 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2512 scm_i_print_complex (double real
, double imag
, SCM port
)
2514 char num_buf
[FLOBUFLEN
];
2515 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2519 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2522 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2523 scm_lfwrite_str (str
, port
);
2524 scm_remember_upto_here_1 (str
);
2529 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2531 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2532 scm_remember_upto_here_1 (exp
);
2533 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2537 /*** END nums->strs ***/
2540 /*** STRINGS -> NUMBERS ***/
2542 /* The following functions implement the conversion from strings to numbers.
2543 * The implementation somehow follows the grammar for numbers as it is given
2544 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2545 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2546 * points should be noted about the implementation:
2547 * * Each function keeps a local index variable 'idx' that points at the
2548 * current position within the parsed string. The global index is only
2549 * updated if the function could parse the corresponding syntactic unit
2551 * * Similarly, the functions keep track of indicators of inexactness ('#',
2552 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2553 * global exactness information is only updated after each part has been
2554 * successfully parsed.
2555 * * Sequences of digits are parsed into temporary variables holding fixnums.
2556 * Only if these fixnums would overflow, the result variables are updated
2557 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2558 * the temporary variables holding the fixnums are cleared, and the process
2559 * starts over again. If for example fixnums were able to store five decimal
2560 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2561 * and the result was computed as 12345 * 100000 + 67890. In other words,
2562 * only every five digits two bignum operations were performed.
2565 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2567 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2569 /* In non ASCII-style encodings the following macro might not work. */
2570 #define XDIGIT2UINT(d) \
2571 (uc_is_property_decimal_digit ((int) (unsigned char) d) \
2573 : uc_tolower ((int) (unsigned char) d) - 'a' + 10)
2576 mem2uinteger (SCM mem
, unsigned int *p_idx
,
2577 unsigned int radix
, enum t_exactness
*p_exactness
)
2579 unsigned int idx
= *p_idx
;
2580 unsigned int hash_seen
= 0;
2581 scm_t_bits shift
= 1;
2583 unsigned int digit_value
;
2586 size_t len
= scm_i_string_length (mem
);
2591 c
= scm_i_string_ref (mem
, idx
);
2592 if (!uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2594 digit_value
= XDIGIT2UINT (c
);
2595 if (digit_value
>= radix
)
2599 result
= SCM_I_MAKINUM (digit_value
);
2602 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2603 if (uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2607 digit_value
= XDIGIT2UINT (c
);
2608 if (digit_value
>= radix
)
2620 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2622 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2624 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2631 shift
= shift
* radix
;
2632 add
= add
* radix
+ digit_value
;
2637 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2639 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2643 *p_exactness
= INEXACT
;
2649 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2650 * covers the parts of the rules that start at a potential point. The value
2651 * of the digits up to the point have been parsed by the caller and are given
2652 * in variable result. The content of *p_exactness indicates, whether a hash
2653 * has already been seen in the digits before the point.
2656 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
2659 mem2decimal_from_point (SCM result
, SCM mem
,
2660 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2662 unsigned int idx
= *p_idx
;
2663 enum t_exactness x
= *p_exactness
;
2664 size_t len
= scm_i_string_length (mem
);
2669 if (scm_i_string_ref (mem
, idx
) == '.')
2671 scm_t_bits shift
= 1;
2673 unsigned int digit_value
;
2674 SCM big_shift
= SCM_I_MAKINUM (1);
2679 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2680 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2685 digit_value
= DIGIT2UINT (c
);
2696 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2698 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2699 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2701 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2709 add
= add
* 10 + digit_value
;
2715 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2716 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2717 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2720 result
= scm_divide (result
, big_shift
);
2722 /* We've seen a decimal point, thus the value is implicitly inexact. */
2734 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2736 switch (scm_i_string_ref (mem
, idx
))
2748 c
= scm_i_string_ref (mem
, idx
);
2756 c
= scm_i_string_ref (mem
, idx
);
2765 c
= scm_i_string_ref (mem
, idx
);
2770 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2774 exponent
= DIGIT2UINT (c
);
2777 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2778 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2781 if (exponent
<= SCM_MAXEXP
)
2782 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2788 if (exponent
> SCM_MAXEXP
)
2790 size_t exp_len
= idx
- start
;
2791 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
2792 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2793 scm_out_of_range ("string->number", exp_num
);
2796 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2798 result
= scm_product (result
, e
);
2800 result
= scm_divide2real (result
, e
);
2802 /* We've seen an exponent, thus the value is implicitly inexact. */
2820 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2823 mem2ureal (SCM mem
, unsigned int *p_idx
,
2824 unsigned int radix
, enum t_exactness
*p_exactness
)
2826 unsigned int idx
= *p_idx
;
2828 size_t len
= scm_i_string_length (mem
);
2830 /* Start off believing that the number will be exact. This changes
2831 to INEXACT if we see a decimal point or a hash. */
2832 enum t_exactness x
= EXACT
;
2837 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
2843 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
2845 /* Cobble up the fractional part. We might want to set the
2846 NaN's mantissa from it. */
2848 mem2uinteger (mem
, &idx
, 10, &x
);
2853 if (scm_i_string_ref (mem
, idx
) == '.')
2857 else if (idx
+ 1 == len
)
2859 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
2862 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
,
2869 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
2870 if (scm_is_false (uinteger
))
2875 else if (scm_i_string_ref (mem
, idx
) == '/')
2883 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
2884 if (scm_is_false (divisor
))
2887 /* both are int/big here, I assume */
2888 result
= scm_i_make_ratio (uinteger
, divisor
);
2890 else if (radix
== 10)
2892 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
2893 if (scm_is_false (result
))
2902 /* Update *p_exactness if the number just read was inexact. This is
2903 important for complex numbers, so that a complex number is
2904 treated as inexact overall if either its real or imaginary part
2910 /* When returning an inexact zero, make sure it is represented as a
2911 floating point value so that we can change its sign.
2913 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2914 result
= scm_from_double (0.0);
2920 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2923 mem2complex (SCM mem
, unsigned int idx
,
2924 unsigned int radix
, enum t_exactness
*p_exactness
)
2929 size_t len
= scm_i_string_length (mem
);
2934 c
= scm_i_string_ref (mem
, idx
);
2949 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2950 if (scm_is_false (ureal
))
2952 /* input must be either +i or -i */
2957 if (scm_i_string_ref (mem
, idx
) == 'i'
2958 || scm_i_string_ref (mem
, idx
) == 'I')
2964 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2971 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2972 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2977 c
= scm_i_string_ref (mem
, idx
);
2981 /* either +<ureal>i or -<ureal>i */
2988 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2991 /* polar input: <real>@<real>. */
3002 c
= scm_i_string_ref (mem
, idx
);
3020 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3021 if (scm_is_false (angle
))
3026 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3027 angle
= scm_difference (angle
, SCM_UNDEFINED
);
3029 result
= scm_make_polar (ureal
, angle
);
3034 /* expecting input matching <real>[+-]<ureal>?i */
3041 int sign
= (c
== '+') ? 1 : -1;
3042 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3044 if (scm_is_false (imag
))
3045 imag
= SCM_I_MAKINUM (sign
);
3046 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3047 imag
= scm_difference (imag
, SCM_UNDEFINED
);
3051 if (scm_i_string_ref (mem
, idx
) != 'i'
3052 && scm_i_string_ref (mem
, idx
) != 'I')
3059 return scm_make_rectangular (ureal
, imag
);
3068 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
3070 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
3073 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
3075 unsigned int idx
= 0;
3076 unsigned int radix
= NO_RADIX
;
3077 enum t_exactness forced_x
= NO_EXACTNESS
;
3078 enum t_exactness implicit_x
= EXACT
;
3080 size_t len
= scm_i_string_length (mem
);
3082 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3083 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
3085 switch (scm_i_string_ref (mem
, idx
+ 1))
3088 if (radix
!= NO_RADIX
)
3093 if (radix
!= NO_RADIX
)
3098 if (forced_x
!= NO_EXACTNESS
)
3103 if (forced_x
!= NO_EXACTNESS
)
3108 if (radix
!= NO_RADIX
)
3113 if (radix
!= NO_RADIX
)
3123 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3124 if (radix
== NO_RADIX
)
3125 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
3127 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
3129 if (scm_is_false (result
))
3135 if (SCM_INEXACTP (result
))
3136 return scm_inexact_to_exact (result
);
3140 if (SCM_INEXACTP (result
))
3143 return scm_exact_to_inexact (result
);
3146 if (implicit_x
== INEXACT
)
3148 if (SCM_INEXACTP (result
))
3151 return scm_exact_to_inexact (result
);
3159 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
3160 unsigned int default_radix
)
3162 SCM str
= scm_from_locale_stringn (mem
, len
);
3164 return scm_i_string_to_number (str
, default_radix
);
3168 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3169 (SCM string
, SCM radix
),
3170 "Return a number of the maximally precise representation\n"
3171 "expressed by the given @var{string}. @var{radix} must be an\n"
3172 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3173 "is a default radix that may be overridden by an explicit radix\n"
3174 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3175 "supplied, then the default radix is 10. If string is not a\n"
3176 "syntactically valid notation for a number, then\n"
3177 "@code{string->number} returns @code{#f}.")
3178 #define FUNC_NAME s_scm_string_to_number
3182 SCM_VALIDATE_STRING (1, string
);
3184 if (SCM_UNBNDP (radix
))
3187 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3189 answer
= scm_i_string_to_number (string
, base
);
3190 scm_remember_upto_here_1 (string
);
3196 /*** END strs->nums ***/
3200 scm_bigequal (SCM x
, SCM y
)
3202 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3203 scm_remember_upto_here_2 (x
, y
);
3204 return scm_from_bool (0 == result
);
3208 scm_real_equalp (SCM x
, SCM y
)
3210 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3214 scm_complex_equalp (SCM x
, SCM y
)
3216 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3217 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3221 scm_i_fraction_equalp (SCM x
, SCM y
)
3223 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3224 SCM_FRACTION_NUMERATOR (y
)))
3225 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3226 SCM_FRACTION_DENOMINATOR (y
))))
3233 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3235 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3237 #define FUNC_NAME s_scm_number_p
3239 return scm_from_bool (SCM_NUMBERP (x
));
3243 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3245 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3246 "otherwise. Note that the sets of real, rational and integer\n"
3247 "values form subsets of the set of complex numbers, i. e. the\n"
3248 "predicate will also be fulfilled if @var{x} is a real,\n"
3249 "rational or integer number.")
3250 #define FUNC_NAME s_scm_complex_p
3252 /* all numbers are complex. */
3253 return scm_number_p (x
);
3257 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3259 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3260 "otherwise. Note that the set of integer values forms a subset of\n"
3261 "the set of real numbers, i. e. the predicate will also be\n"
3262 "fulfilled if @var{x} is an integer number.")
3263 #define FUNC_NAME s_scm_real_p
3265 /* we can't represent irrational numbers. */
3266 return scm_rational_p (x
);
3270 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3272 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3273 "otherwise. Note that the set of integer values forms a subset of\n"
3274 "the set of rational numbers, i. e. the predicate will also be\n"
3275 "fulfilled if @var{x} is an integer number.")
3276 #define FUNC_NAME s_scm_rational_p
3278 if (SCM_I_INUMP (x
))
3280 else if (SCM_IMP (x
))
3282 else if (SCM_BIGP (x
))
3284 else if (SCM_FRACTIONP (x
))
3286 else if (SCM_REALP (x
))
3287 /* due to their limited precision, all floating point numbers are
3288 rational as well. */
3295 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3297 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3299 #define FUNC_NAME s_scm_integer_p
3302 if (SCM_I_INUMP (x
))
3308 if (!SCM_INEXACTP (x
))
3310 if (SCM_COMPLEXP (x
))
3312 r
= SCM_REAL_VALUE (x
);
3313 /* +/-inf passes r==floor(r), making those #t */
3321 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3323 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3325 #define FUNC_NAME s_scm_inexact_p
3327 if (SCM_INEXACTP (x
))
3329 if (SCM_NUMBERP (x
))
3331 SCM_WRONG_TYPE_ARG (1, x
);
3336 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
3337 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
3338 (SCM x
, SCM y
, SCM rest
),
3339 "Return @code{#t} if all parameters are numerically equal.")
3340 #define FUNC_NAME s_scm_i_num_eq_p
3342 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3344 while (!scm_is_null (rest
))
3346 if (scm_is_false (scm_num_eq_p (x
, y
)))
3350 rest
= scm_cdr (rest
);
3352 return scm_num_eq_p (x
, y
);
3356 scm_num_eq_p (SCM x
, SCM y
)
3359 if (SCM_I_INUMP (x
))
3361 long xx
= SCM_I_INUM (x
);
3362 if (SCM_I_INUMP (y
))
3364 long yy
= SCM_I_INUM (y
);
3365 return scm_from_bool (xx
== yy
);
3367 else if (SCM_BIGP (y
))
3369 else if (SCM_REALP (y
))
3371 /* On a 32-bit system an inum fits a double, we can cast the inum
3372 to a double and compare.
3374 But on a 64-bit system an inum is bigger than a double and
3375 casting it to a double (call that dxx) will round. dxx is at
3376 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3377 an integer and fits a long. So we cast yy to a long and
3378 compare with plain xx.
3380 An alternative (for any size system actually) would be to check
3381 yy is an integer (with floor) and is in range of an inum
3382 (compare against appropriate powers of 2) then test
3383 xx==(long)yy. It's just a matter of which casts/comparisons
3384 might be fastest or easiest for the cpu. */
3386 double yy
= SCM_REAL_VALUE (y
);
3387 return scm_from_bool ((double) xx
== yy
3388 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3389 || xx
== (long) yy
));
3391 else if (SCM_COMPLEXP (y
))
3392 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3393 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3394 else if (SCM_FRACTIONP (y
))
3397 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3399 else if (SCM_BIGP (x
))
3401 if (SCM_I_INUMP (y
))
3403 else if (SCM_BIGP (y
))
3405 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3406 scm_remember_upto_here_2 (x
, y
);
3407 return scm_from_bool (0 == cmp
);
3409 else if (SCM_REALP (y
))
3412 if (xisnan (SCM_REAL_VALUE (y
)))
3414 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3415 scm_remember_upto_here_1 (x
);
3416 return scm_from_bool (0 == cmp
);
3418 else if (SCM_COMPLEXP (y
))
3421 if (0.0 != SCM_COMPLEX_IMAG (y
))
3423 if (xisnan (SCM_COMPLEX_REAL (y
)))
3425 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3426 scm_remember_upto_here_1 (x
);
3427 return scm_from_bool (0 == cmp
);
3429 else if (SCM_FRACTIONP (y
))
3432 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3434 else if (SCM_REALP (x
))
3436 double xx
= SCM_REAL_VALUE (x
);
3437 if (SCM_I_INUMP (y
))
3439 /* see comments with inum/real above */
3440 long yy
= SCM_I_INUM (y
);
3441 return scm_from_bool (xx
== (double) yy
3442 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3443 || (long) xx
== yy
));
3445 else if (SCM_BIGP (y
))
3448 if (xisnan (SCM_REAL_VALUE (x
)))
3450 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3451 scm_remember_upto_here_1 (y
);
3452 return scm_from_bool (0 == cmp
);
3454 else if (SCM_REALP (y
))
3455 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3456 else if (SCM_COMPLEXP (y
))
3457 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3458 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3459 else if (SCM_FRACTIONP (y
))
3461 double xx
= SCM_REAL_VALUE (x
);
3465 return scm_from_bool (xx
< 0.0);
3466 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3470 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3472 else if (SCM_COMPLEXP (x
))
3474 if (SCM_I_INUMP (y
))
3475 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3476 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3477 else if (SCM_BIGP (y
))
3480 if (0.0 != SCM_COMPLEX_IMAG (x
))
3482 if (xisnan (SCM_COMPLEX_REAL (x
)))
3484 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3485 scm_remember_upto_here_1 (y
);
3486 return scm_from_bool (0 == cmp
);
3488 else if (SCM_REALP (y
))
3489 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3490 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3491 else if (SCM_COMPLEXP (y
))
3492 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3493 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3494 else if (SCM_FRACTIONP (y
))
3497 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3499 xx
= SCM_COMPLEX_REAL (x
);
3503 return scm_from_bool (xx
< 0.0);
3504 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3508 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3510 else if (SCM_FRACTIONP (x
))
3512 if (SCM_I_INUMP (y
))
3514 else if (SCM_BIGP (y
))
3516 else if (SCM_REALP (y
))
3518 double yy
= SCM_REAL_VALUE (y
);
3522 return scm_from_bool (0.0 < yy
);
3523 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3526 else if (SCM_COMPLEXP (y
))
3529 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3531 yy
= SCM_COMPLEX_REAL (y
);
3535 return scm_from_bool (0.0 < yy
);
3536 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3539 else if (SCM_FRACTIONP (y
))
3540 return scm_i_fraction_equalp (x
, y
);
3542 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3545 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
3549 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3550 done are good for inums, but for bignums an answer can almost always be
3551 had by just examining a few high bits of the operands, as done by GMP in
3552 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3553 of the float exponent to take into account. */
3555 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
3556 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
3557 (SCM x
, SCM y
, SCM rest
),
3558 "Return @code{#t} if the list of parameters is monotonically\n"
3560 #define FUNC_NAME s_scm_i_num_less_p
3562 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3564 while (!scm_is_null (rest
))
3566 if (scm_is_false (scm_less_p (x
, y
)))
3570 rest
= scm_cdr (rest
);
3572 return scm_less_p (x
, y
);
3576 scm_less_p (SCM x
, SCM y
)
3579 if (SCM_I_INUMP (x
))
3581 long xx
= SCM_I_INUM (x
);
3582 if (SCM_I_INUMP (y
))
3584 long yy
= SCM_I_INUM (y
);
3585 return scm_from_bool (xx
< yy
);
3587 else if (SCM_BIGP (y
))
3589 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3590 scm_remember_upto_here_1 (y
);
3591 return scm_from_bool (sgn
> 0);
3593 else if (SCM_REALP (y
))
3594 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3595 else if (SCM_FRACTIONP (y
))
3597 /* "x < a/b" becomes "x*b < a" */
3599 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3600 y
= SCM_FRACTION_NUMERATOR (y
);
3604 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3606 else if (SCM_BIGP (x
))
3608 if (SCM_I_INUMP (y
))
3610 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3611 scm_remember_upto_here_1 (x
);
3612 return scm_from_bool (sgn
< 0);
3614 else if (SCM_BIGP (y
))
3616 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3617 scm_remember_upto_here_2 (x
, y
);
3618 return scm_from_bool (cmp
< 0);
3620 else if (SCM_REALP (y
))
3623 if (xisnan (SCM_REAL_VALUE (y
)))
3625 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3626 scm_remember_upto_here_1 (x
);
3627 return scm_from_bool (cmp
< 0);
3629 else if (SCM_FRACTIONP (y
))
3632 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3634 else if (SCM_REALP (x
))
3636 if (SCM_I_INUMP (y
))
3637 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3638 else if (SCM_BIGP (y
))
3641 if (xisnan (SCM_REAL_VALUE (x
)))
3643 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3644 scm_remember_upto_here_1 (y
);
3645 return scm_from_bool (cmp
> 0);
3647 else if (SCM_REALP (y
))
3648 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3649 else if (SCM_FRACTIONP (y
))
3651 double xx
= SCM_REAL_VALUE (x
);
3655 return scm_from_bool (xx
< 0.0);
3656 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3660 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3662 else if (SCM_FRACTIONP (x
))
3664 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3666 /* "a/b < y" becomes "a < y*b" */
3667 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3668 x
= SCM_FRACTION_NUMERATOR (x
);
3671 else if (SCM_REALP (y
))
3673 double yy
= SCM_REAL_VALUE (y
);
3677 return scm_from_bool (0.0 < yy
);
3678 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3681 else if (SCM_FRACTIONP (y
))
3683 /* "a/b < c/d" becomes "a*d < c*b" */
3684 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3685 SCM_FRACTION_DENOMINATOR (y
));
3686 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3687 SCM_FRACTION_DENOMINATOR (x
));
3693 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3696 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
3700 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
3701 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
3702 (SCM x
, SCM y
, SCM rest
),
3703 "Return @code{#t} if the list of parameters is monotonically\n"
3705 #define FUNC_NAME s_scm_i_num_gr_p
3707 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3709 while (!scm_is_null (rest
))
3711 if (scm_is_false (scm_gr_p (x
, y
)))
3715 rest
= scm_cdr (rest
);
3717 return scm_gr_p (x
, y
);
3720 #define FUNC_NAME s_scm_i_num_gr_p
3722 scm_gr_p (SCM x
, SCM y
)
3724 if (!SCM_NUMBERP (x
))
3725 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3726 else if (!SCM_NUMBERP (y
))
3727 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3729 return scm_less_p (y
, x
);
3734 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
3735 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
3736 (SCM x
, SCM y
, SCM rest
),
3737 "Return @code{#t} if the list of parameters is monotonically\n"
3739 #define FUNC_NAME s_scm_i_num_leq_p
3741 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3743 while (!scm_is_null (rest
))
3745 if (scm_is_false (scm_leq_p (x
, y
)))
3749 rest
= scm_cdr (rest
);
3751 return scm_leq_p (x
, y
);
3754 #define FUNC_NAME s_scm_i_num_leq_p
3756 scm_leq_p (SCM x
, SCM y
)
3758 if (!SCM_NUMBERP (x
))
3759 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3760 else if (!SCM_NUMBERP (y
))
3761 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3762 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3765 return scm_not (scm_less_p (y
, x
));
3770 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
3771 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
3772 (SCM x
, SCM y
, SCM rest
),
3773 "Return @code{#t} if the list of parameters is monotonically\n"
3775 #define FUNC_NAME s_scm_i_num_geq_p
3777 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3779 while (!scm_is_null (rest
))
3781 if (scm_is_false (scm_geq_p (x
, y
)))
3785 rest
= scm_cdr (rest
);
3787 return scm_geq_p (x
, y
);
3790 #define FUNC_NAME s_scm_i_num_geq_p
3792 scm_geq_p (SCM x
, SCM y
)
3794 if (!SCM_NUMBERP (x
))
3795 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3796 else if (!SCM_NUMBERP (y
))
3797 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3798 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3801 return scm_not (scm_less_p (x
, y
));
3806 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3807 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3813 if (SCM_I_INUMP (z
))
3814 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3815 else if (SCM_BIGP (z
))
3817 else if (SCM_REALP (z
))
3818 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3819 else if (SCM_COMPLEXP (z
))
3820 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3821 && SCM_COMPLEX_IMAG (z
) == 0.0);
3822 else if (SCM_FRACTIONP (z
))
3825 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3829 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3830 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3834 scm_positive_p (SCM x
)
3836 if (SCM_I_INUMP (x
))
3837 return scm_from_bool (SCM_I_INUM (x
) > 0);
3838 else if (SCM_BIGP (x
))
3840 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3841 scm_remember_upto_here_1 (x
);
3842 return scm_from_bool (sgn
> 0);
3844 else if (SCM_REALP (x
))
3845 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3846 else if (SCM_FRACTIONP (x
))
3847 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3849 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3853 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3854 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3858 scm_negative_p (SCM x
)
3860 if (SCM_I_INUMP (x
))
3861 return scm_from_bool (SCM_I_INUM (x
) < 0);
3862 else if (SCM_BIGP (x
))
3864 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3865 scm_remember_upto_here_1 (x
);
3866 return scm_from_bool (sgn
< 0);
3868 else if (SCM_REALP (x
))
3869 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3870 else if (SCM_FRACTIONP (x
))
3871 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3873 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3877 /* scm_min and scm_max return an inexact when either argument is inexact, as
3878 required by r5rs. On that basis, for exact/inexact combinations the
3879 exact is converted to inexact to compare and possibly return. This is
3880 unlike scm_less_p above which takes some trouble to preserve all bits in
3881 its test, such trouble is not required for min and max. */
3883 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
3884 (SCM x
, SCM y
, SCM rest
),
3885 "Return the maximum of all parameter values.")
3886 #define FUNC_NAME s_scm_i_max
3888 while (!scm_is_null (rest
))
3889 { x
= scm_max (x
, y
);
3891 rest
= scm_cdr (rest
);
3893 return scm_max (x
, y
);
3897 #define s_max s_scm_i_max
3898 #define g_max g_scm_i_max
3901 scm_max (SCM x
, SCM y
)
3906 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3907 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3910 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3913 if (SCM_I_INUMP (x
))
3915 long xx
= SCM_I_INUM (x
);
3916 if (SCM_I_INUMP (y
))
3918 long yy
= SCM_I_INUM (y
);
3919 return (xx
< yy
) ? y
: x
;
3921 else if (SCM_BIGP (y
))
3923 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3924 scm_remember_upto_here_1 (y
);
3925 return (sgn
< 0) ? x
: y
;
3927 else if (SCM_REALP (y
))
3930 /* if y==NaN then ">" is false and we return NaN */
3931 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3933 else if (SCM_FRACTIONP (y
))
3936 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3939 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3941 else if (SCM_BIGP (x
))
3943 if (SCM_I_INUMP (y
))
3945 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3946 scm_remember_upto_here_1 (x
);
3947 return (sgn
< 0) ? y
: x
;
3949 else if (SCM_BIGP (y
))
3951 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3952 scm_remember_upto_here_2 (x
, y
);
3953 return (cmp
> 0) ? x
: y
;
3955 else if (SCM_REALP (y
))
3957 /* if y==NaN then xx>yy is false, so we return the NaN y */
3960 xx
= scm_i_big2dbl (x
);
3961 yy
= SCM_REAL_VALUE (y
);
3962 return (xx
> yy
? scm_from_double (xx
) : y
);
3964 else if (SCM_FRACTIONP (y
))
3969 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3971 else if (SCM_REALP (x
))
3973 if (SCM_I_INUMP (y
))
3975 double z
= SCM_I_INUM (y
);
3976 /* if x==NaN then "<" is false and we return NaN */
3977 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3979 else if (SCM_BIGP (y
))
3984 else if (SCM_REALP (y
))
3986 /* if x==NaN then our explicit check means we return NaN
3987 if y==NaN then ">" is false and we return NaN
3988 calling isnan is unavoidable, since it's the only way to know
3989 which of x or y causes any compares to be false */
3990 double xx
= SCM_REAL_VALUE (x
);
3991 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3993 else if (SCM_FRACTIONP (y
))
3995 double yy
= scm_i_fraction2double (y
);
3996 double xx
= SCM_REAL_VALUE (x
);
3997 return (xx
< yy
) ? scm_from_double (yy
) : x
;
4000 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4002 else if (SCM_FRACTIONP (x
))
4004 if (SCM_I_INUMP (y
))
4008 else if (SCM_BIGP (y
))
4012 else if (SCM_REALP (y
))
4014 double xx
= scm_i_fraction2double (x
);
4015 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
4017 else if (SCM_FRACTIONP (y
))
4022 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4025 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
4029 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
4030 (SCM x
, SCM y
, SCM rest
),
4031 "Return the minimum of all parameter values.")
4032 #define FUNC_NAME s_scm_i_min
4034 while (!scm_is_null (rest
))
4035 { x
= scm_min (x
, y
);
4037 rest
= scm_cdr (rest
);
4039 return scm_min (x
, y
);
4043 #define s_min s_scm_i_min
4044 #define g_min g_scm_i_min
4047 scm_min (SCM x
, SCM y
)
4052 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
4053 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
4056 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
4059 if (SCM_I_INUMP (x
))
4061 long xx
= SCM_I_INUM (x
);
4062 if (SCM_I_INUMP (y
))
4064 long yy
= SCM_I_INUM (y
);
4065 return (xx
< yy
) ? x
: y
;
4067 else if (SCM_BIGP (y
))
4069 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4070 scm_remember_upto_here_1 (y
);
4071 return (sgn
< 0) ? y
: x
;
4073 else if (SCM_REALP (y
))
4076 /* if y==NaN then "<" is false and we return NaN */
4077 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
4079 else if (SCM_FRACTIONP (y
))
4082 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
4085 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4087 else if (SCM_BIGP (x
))
4089 if (SCM_I_INUMP (y
))
4091 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4092 scm_remember_upto_here_1 (x
);
4093 return (sgn
< 0) ? x
: y
;
4095 else if (SCM_BIGP (y
))
4097 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4098 scm_remember_upto_here_2 (x
, y
);
4099 return (cmp
> 0) ? y
: x
;
4101 else if (SCM_REALP (y
))
4103 /* if y==NaN then xx<yy is false, so we return the NaN y */
4106 xx
= scm_i_big2dbl (x
);
4107 yy
= SCM_REAL_VALUE (y
);
4108 return (xx
< yy
? scm_from_double (xx
) : y
);
4110 else if (SCM_FRACTIONP (y
))
4115 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4117 else if (SCM_REALP (x
))
4119 if (SCM_I_INUMP (y
))
4121 double z
= SCM_I_INUM (y
);
4122 /* if x==NaN then "<" is false and we return NaN */
4123 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
4125 else if (SCM_BIGP (y
))
4130 else if (SCM_REALP (y
))
4132 /* if x==NaN then our explicit check means we return NaN
4133 if y==NaN then "<" is false and we return NaN
4134 calling isnan is unavoidable, since it's the only way to know
4135 which of x or y causes any compares to be false */
4136 double xx
= SCM_REAL_VALUE (x
);
4137 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
4139 else if (SCM_FRACTIONP (y
))
4141 double yy
= scm_i_fraction2double (y
);
4142 double xx
= SCM_REAL_VALUE (x
);
4143 return (yy
< xx
) ? scm_from_double (yy
) : x
;
4146 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4148 else if (SCM_FRACTIONP (x
))
4150 if (SCM_I_INUMP (y
))
4154 else if (SCM_BIGP (y
))
4158 else if (SCM_REALP (y
))
4160 double xx
= scm_i_fraction2double (x
);
4161 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
4163 else if (SCM_FRACTIONP (y
))
4168 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4171 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
4175 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
4176 (SCM x
, SCM y
, SCM rest
),
4177 "Return the sum of all parameter values. Return 0 if called without\n"
4179 #define FUNC_NAME s_scm_i_sum
4181 while (!scm_is_null (rest
))
4182 { x
= scm_sum (x
, y
);
4184 rest
= scm_cdr (rest
);
4186 return scm_sum (x
, y
);
4190 #define s_sum s_scm_i_sum
4191 #define g_sum g_scm_i_sum
4194 scm_sum (SCM x
, SCM y
)
4196 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4198 if (SCM_NUMBERP (x
)) return x
;
4199 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4200 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4203 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4205 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4207 long xx
= SCM_I_INUM (x
);
4208 long yy
= SCM_I_INUM (y
);
4209 long int z
= xx
+ yy
;
4210 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
4212 else if (SCM_BIGP (y
))
4217 else if (SCM_REALP (y
))
4219 long int xx
= SCM_I_INUM (x
);
4220 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4222 else if (SCM_COMPLEXP (y
))
4224 long int xx
= SCM_I_INUM (x
);
4225 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4226 SCM_COMPLEX_IMAG (y
));
4228 else if (SCM_FRACTIONP (y
))
4229 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4230 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4231 SCM_FRACTION_DENOMINATOR (y
));
4233 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4234 } else if (SCM_BIGP (x
))
4236 if (SCM_I_INUMP (y
))
4241 inum
= SCM_I_INUM (y
);
4244 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4247 SCM result
= scm_i_mkbig ();
4248 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4249 scm_remember_upto_here_1 (x
);
4250 /* we know the result will have to be a bignum */
4253 return scm_i_normbig (result
);
4257 SCM result
= scm_i_mkbig ();
4258 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4259 scm_remember_upto_here_1 (x
);
4260 /* we know the result will have to be a bignum */
4263 return scm_i_normbig (result
);
4266 else if (SCM_BIGP (y
))
4268 SCM result
= scm_i_mkbig ();
4269 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4270 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4271 mpz_add (SCM_I_BIG_MPZ (result
),
4274 scm_remember_upto_here_2 (x
, y
);
4275 /* we know the result will have to be a bignum */
4278 return scm_i_normbig (result
);
4280 else if (SCM_REALP (y
))
4282 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4283 scm_remember_upto_here_1 (x
);
4284 return scm_from_double (result
);
4286 else if (SCM_COMPLEXP (y
))
4288 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4289 + SCM_COMPLEX_REAL (y
));
4290 scm_remember_upto_here_1 (x
);
4291 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4293 else if (SCM_FRACTIONP (y
))
4294 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4295 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4296 SCM_FRACTION_DENOMINATOR (y
));
4298 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4300 else if (SCM_REALP (x
))
4302 if (SCM_I_INUMP (y
))
4303 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4304 else if (SCM_BIGP (y
))
4306 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4307 scm_remember_upto_here_1 (y
);
4308 return scm_from_double (result
);
4310 else if (SCM_REALP (y
))
4311 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4312 else if (SCM_COMPLEXP (y
))
4313 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4314 SCM_COMPLEX_IMAG (y
));
4315 else if (SCM_FRACTIONP (y
))
4316 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4318 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4320 else if (SCM_COMPLEXP (x
))
4322 if (SCM_I_INUMP (y
))
4323 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4324 SCM_COMPLEX_IMAG (x
));
4325 else if (SCM_BIGP (y
))
4327 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4328 + SCM_COMPLEX_REAL (x
));
4329 scm_remember_upto_here_1 (y
);
4330 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4332 else if (SCM_REALP (y
))
4333 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4334 SCM_COMPLEX_IMAG (x
));
4335 else if (SCM_COMPLEXP (y
))
4336 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4337 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4338 else if (SCM_FRACTIONP (y
))
4339 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4340 SCM_COMPLEX_IMAG (x
));
4342 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4344 else if (SCM_FRACTIONP (x
))
4346 if (SCM_I_INUMP (y
))
4347 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4348 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4349 SCM_FRACTION_DENOMINATOR (x
));
4350 else if (SCM_BIGP (y
))
4351 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4352 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4353 SCM_FRACTION_DENOMINATOR (x
));
4354 else if (SCM_REALP (y
))
4355 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4356 else if (SCM_COMPLEXP (y
))
4357 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4358 SCM_COMPLEX_IMAG (y
));
4359 else if (SCM_FRACTIONP (y
))
4360 /* a/b + c/d = (ad + bc) / bd */
4361 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4362 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4363 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4365 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4368 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4372 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4374 "Return @math{@var{x}+1}.")
4375 #define FUNC_NAME s_scm_oneplus
4377 return scm_sum (x
, SCM_I_MAKINUM (1));
4382 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
4383 (SCM x
, SCM y
, SCM rest
),
4384 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
4385 "the sum of all but the first argument are subtracted from the first\n"
4387 #define FUNC_NAME s_scm_i_difference
4389 while (!scm_is_null (rest
))
4390 { x
= scm_difference (x
, y
);
4392 rest
= scm_cdr (rest
);
4394 return scm_difference (x
, y
);
4398 #define s_difference s_scm_i_difference
4399 #define g_difference g_scm_i_difference
4402 scm_difference (SCM x
, SCM y
)
4403 #define FUNC_NAME s_difference
4405 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4408 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4410 if (SCM_I_INUMP (x
))
4412 long xx
= -SCM_I_INUM (x
);
4413 if (SCM_FIXABLE (xx
))
4414 return SCM_I_MAKINUM (xx
);
4416 return scm_i_long2big (xx
);
4418 else if (SCM_BIGP (x
))
4419 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4420 bignum, but negating that gives a fixnum. */
4421 return scm_i_normbig (scm_i_clonebig (x
, 0));
4422 else if (SCM_REALP (x
))
4423 return scm_from_double (-SCM_REAL_VALUE (x
));
4424 else if (SCM_COMPLEXP (x
))
4425 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4426 -SCM_COMPLEX_IMAG (x
));
4427 else if (SCM_FRACTIONP (x
))
4428 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4429 SCM_FRACTION_DENOMINATOR (x
));
4431 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4434 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4436 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4438 long int xx
= SCM_I_INUM (x
);
4439 long int yy
= SCM_I_INUM (y
);
4440 long int z
= xx
- yy
;
4441 if (SCM_FIXABLE (z
))
4442 return SCM_I_MAKINUM (z
);
4444 return scm_i_long2big (z
);
4446 else if (SCM_BIGP (y
))
4448 /* inum-x - big-y */
4449 long xx
= SCM_I_INUM (x
);
4452 return scm_i_clonebig (y
, 0);
4455 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4456 SCM result
= scm_i_mkbig ();
4459 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4462 /* x - y == -(y + -x) */
4463 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4464 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4466 scm_remember_upto_here_1 (y
);
4468 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4469 /* we know the result will have to be a bignum */
4472 return scm_i_normbig (result
);
4475 else if (SCM_REALP (y
))
4477 long int xx
= SCM_I_INUM (x
);
4478 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4480 else if (SCM_COMPLEXP (y
))
4482 long int xx
= SCM_I_INUM (x
);
4483 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4484 - SCM_COMPLEX_IMAG (y
));
4486 else if (SCM_FRACTIONP (y
))
4487 /* a - b/c = (ac - b) / c */
4488 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4489 SCM_FRACTION_NUMERATOR (y
)),
4490 SCM_FRACTION_DENOMINATOR (y
));
4492 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4494 else if (SCM_BIGP (x
))
4496 if (SCM_I_INUMP (y
))
4498 /* big-x - inum-y */
4499 long yy
= SCM_I_INUM (y
);
4500 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4502 scm_remember_upto_here_1 (x
);
4504 return (SCM_FIXABLE (-yy
) ?
4505 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4508 SCM result
= scm_i_mkbig ();
4511 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4513 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4514 scm_remember_upto_here_1 (x
);
4516 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4517 /* we know the result will have to be a bignum */
4520 return scm_i_normbig (result
);
4523 else if (SCM_BIGP (y
))
4525 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4526 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4527 SCM result
= scm_i_mkbig ();
4528 mpz_sub (SCM_I_BIG_MPZ (result
),
4531 scm_remember_upto_here_2 (x
, y
);
4532 /* we know the result will have to be a bignum */
4533 if ((sgn_x
== 1) && (sgn_y
== -1))
4535 if ((sgn_x
== -1) && (sgn_y
== 1))
4537 return scm_i_normbig (result
);
4539 else if (SCM_REALP (y
))
4541 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4542 scm_remember_upto_here_1 (x
);
4543 return scm_from_double (result
);
4545 else if (SCM_COMPLEXP (y
))
4547 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4548 - SCM_COMPLEX_REAL (y
));
4549 scm_remember_upto_here_1 (x
);
4550 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4552 else if (SCM_FRACTIONP (y
))
4553 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4554 SCM_FRACTION_NUMERATOR (y
)),
4555 SCM_FRACTION_DENOMINATOR (y
));
4556 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4558 else if (SCM_REALP (x
))
4560 if (SCM_I_INUMP (y
))
4561 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4562 else if (SCM_BIGP (y
))
4564 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4565 scm_remember_upto_here_1 (x
);
4566 return scm_from_double (result
);
4568 else if (SCM_REALP (y
))
4569 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4570 else if (SCM_COMPLEXP (y
))
4571 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4572 -SCM_COMPLEX_IMAG (y
));
4573 else if (SCM_FRACTIONP (y
))
4574 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4576 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4578 else if (SCM_COMPLEXP (x
))
4580 if (SCM_I_INUMP (y
))
4581 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4582 SCM_COMPLEX_IMAG (x
));
4583 else if (SCM_BIGP (y
))
4585 double real_part
= (SCM_COMPLEX_REAL (x
)
4586 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4587 scm_remember_upto_here_1 (x
);
4588 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4590 else if (SCM_REALP (y
))
4591 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4592 SCM_COMPLEX_IMAG (x
));
4593 else if (SCM_COMPLEXP (y
))
4594 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4595 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4596 else if (SCM_FRACTIONP (y
))
4597 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4598 SCM_COMPLEX_IMAG (x
));
4600 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4602 else if (SCM_FRACTIONP (x
))
4604 if (SCM_I_INUMP (y
))
4605 /* a/b - c = (a - cb) / b */
4606 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4607 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4608 SCM_FRACTION_DENOMINATOR (x
));
4609 else if (SCM_BIGP (y
))
4610 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4611 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4612 SCM_FRACTION_DENOMINATOR (x
));
4613 else if (SCM_REALP (y
))
4614 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4615 else if (SCM_COMPLEXP (y
))
4616 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4617 -SCM_COMPLEX_IMAG (y
));
4618 else if (SCM_FRACTIONP (y
))
4619 /* a/b - c/d = (ad - bc) / bd */
4620 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4621 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4622 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4624 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4627 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4632 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4634 "Return @math{@var{x}-1}.")
4635 #define FUNC_NAME s_scm_oneminus
4637 return scm_difference (x
, SCM_I_MAKINUM (1));
4642 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
4643 (SCM x
, SCM y
, SCM rest
),
4644 "Return the product of all arguments. If called without arguments,\n"
4646 #define FUNC_NAME s_scm_i_product
4648 while (!scm_is_null (rest
))
4649 { x
= scm_product (x
, y
);
4651 rest
= scm_cdr (rest
);
4653 return scm_product (x
, y
);
4657 #define s_product s_scm_i_product
4658 #define g_product g_scm_i_product
4661 scm_product (SCM x
, SCM y
)
4663 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4666 return SCM_I_MAKINUM (1L);
4667 else if (SCM_NUMBERP (x
))
4670 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4673 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4678 xx
= SCM_I_INUM (x
);
4682 case 0: return x
; break;
4683 case 1: return y
; break;
4686 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4688 long yy
= SCM_I_INUM (y
);
4690 SCM k
= SCM_I_MAKINUM (kk
);
4691 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4695 SCM result
= scm_i_long2big (xx
);
4696 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4697 return scm_i_normbig (result
);
4700 else if (SCM_BIGP (y
))
4702 SCM result
= scm_i_mkbig ();
4703 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4704 scm_remember_upto_here_1 (y
);
4707 else if (SCM_REALP (y
))
4708 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4709 else if (SCM_COMPLEXP (y
))
4710 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4711 xx
* SCM_COMPLEX_IMAG (y
));
4712 else if (SCM_FRACTIONP (y
))
4713 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4714 SCM_FRACTION_DENOMINATOR (y
));
4716 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4718 else if (SCM_BIGP (x
))
4720 if (SCM_I_INUMP (y
))
4725 else if (SCM_BIGP (y
))
4727 SCM result
= scm_i_mkbig ();
4728 mpz_mul (SCM_I_BIG_MPZ (result
),
4731 scm_remember_upto_here_2 (x
, y
);
4734 else if (SCM_REALP (y
))
4736 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4737 scm_remember_upto_here_1 (x
);
4738 return scm_from_double (result
);
4740 else if (SCM_COMPLEXP (y
))
4742 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4743 scm_remember_upto_here_1 (x
);
4744 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4745 z
* SCM_COMPLEX_IMAG (y
));
4747 else if (SCM_FRACTIONP (y
))
4748 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4749 SCM_FRACTION_DENOMINATOR (y
));
4751 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4753 else if (SCM_REALP (x
))
4755 if (SCM_I_INUMP (y
))
4757 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4758 if (scm_is_eq (y
, SCM_INUM0
))
4760 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4762 else if (SCM_BIGP (y
))
4764 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4765 scm_remember_upto_here_1 (y
);
4766 return scm_from_double (result
);
4768 else if (SCM_REALP (y
))
4769 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4770 else if (SCM_COMPLEXP (y
))
4771 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4772 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4773 else if (SCM_FRACTIONP (y
))
4774 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4776 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4778 else if (SCM_COMPLEXP (x
))
4780 if (SCM_I_INUMP (y
))
4782 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4783 if (scm_is_eq (y
, SCM_INUM0
))
4785 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4786 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4788 else if (SCM_BIGP (y
))
4790 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4791 scm_remember_upto_here_1 (y
);
4792 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4793 z
* SCM_COMPLEX_IMAG (x
));
4795 else if (SCM_REALP (y
))
4796 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4797 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4798 else if (SCM_COMPLEXP (y
))
4800 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4801 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4802 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4803 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4805 else if (SCM_FRACTIONP (y
))
4807 double yy
= scm_i_fraction2double (y
);
4808 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4809 yy
* SCM_COMPLEX_IMAG (x
));
4812 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4814 else if (SCM_FRACTIONP (x
))
4816 if (SCM_I_INUMP (y
))
4817 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4818 SCM_FRACTION_DENOMINATOR (x
));
4819 else if (SCM_BIGP (y
))
4820 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4821 SCM_FRACTION_DENOMINATOR (x
));
4822 else if (SCM_REALP (y
))
4823 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4824 else if (SCM_COMPLEXP (y
))
4826 double xx
= scm_i_fraction2double (x
);
4827 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4828 xx
* SCM_COMPLEX_IMAG (y
));
4830 else if (SCM_FRACTIONP (y
))
4831 /* a/b * c/d = ac / bd */
4832 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4833 SCM_FRACTION_NUMERATOR (y
)),
4834 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4835 SCM_FRACTION_DENOMINATOR (y
)));
4837 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4840 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4843 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4844 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4845 #define ALLOW_DIVIDE_BY_ZERO
4846 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4849 /* The code below for complex division is adapted from the GNU
4850 libstdc++, which adapted it from f2c's libF77, and is subject to
4853 /****************************************************************
4854 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4856 Permission to use, copy, modify, and distribute this software
4857 and its documentation for any purpose and without fee is hereby
4858 granted, provided that the above copyright notice appear in all
4859 copies and that both that the copyright notice and this
4860 permission notice and warranty disclaimer appear in supporting
4861 documentation, and that the names of AT&T Bell Laboratories or
4862 Bellcore or any of their entities not be used in advertising or
4863 publicity pertaining to distribution of the software without
4864 specific, written prior permission.
4866 AT&T and Bellcore disclaim all warranties with regard to this
4867 software, including all implied warranties of merchantability
4868 and fitness. In no event shall AT&T or Bellcore be liable for
4869 any special, indirect or consequential damages or any damages
4870 whatsoever resulting from loss of use, data or profits, whether
4871 in an action of contract, negligence or other tortious action,
4872 arising out of or in connection with the use or performance of
4874 ****************************************************************/
4876 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
4877 (SCM x
, SCM y
, SCM rest
),
4878 "Divide the first argument by the product of the remaining\n"
4879 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
4881 #define FUNC_NAME s_scm_i_divide
4883 while (!scm_is_null (rest
))
4884 { x
= scm_divide (x
, y
);
4886 rest
= scm_cdr (rest
);
4888 return scm_divide (x
, y
);
4892 #define s_divide s_scm_i_divide
4893 #define g_divide g_scm_i_divide
4896 do_divide (SCM x
, SCM y
, int inexact
)
4897 #define FUNC_NAME s_divide
4901 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4904 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4905 else if (SCM_I_INUMP (x
))
4907 long xx
= SCM_I_INUM (x
);
4908 if (xx
== 1 || xx
== -1)
4910 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4912 scm_num_overflow (s_divide
);
4917 return scm_from_double (1.0 / (double) xx
);
4918 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4921 else if (SCM_BIGP (x
))
4924 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4925 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4927 else if (SCM_REALP (x
))
4929 double xx
= SCM_REAL_VALUE (x
);
4930 #ifndef ALLOW_DIVIDE_BY_ZERO
4932 scm_num_overflow (s_divide
);
4935 return scm_from_double (1.0 / xx
);
4937 else if (SCM_COMPLEXP (x
))
4939 double r
= SCM_COMPLEX_REAL (x
);
4940 double i
= SCM_COMPLEX_IMAG (x
);
4941 if (fabs(r
) <= fabs(i
))
4944 double d
= i
* (1.0 + t
* t
);
4945 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4950 double d
= r
* (1.0 + t
* t
);
4951 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4954 else if (SCM_FRACTIONP (x
))
4955 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4956 SCM_FRACTION_NUMERATOR (x
));
4958 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4961 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4963 long xx
= SCM_I_INUM (x
);
4964 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4966 long yy
= SCM_I_INUM (y
);
4969 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4970 scm_num_overflow (s_divide
);
4972 return scm_from_double ((double) xx
/ (double) yy
);
4975 else if (xx
% yy
!= 0)
4978 return scm_from_double ((double) xx
/ (double) yy
);
4979 else return scm_i_make_ratio (x
, y
);
4984 if (SCM_FIXABLE (z
))
4985 return SCM_I_MAKINUM (z
);
4987 return scm_i_long2big (z
);
4990 else if (SCM_BIGP (y
))
4993 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4994 else return scm_i_make_ratio (x
, y
);
4996 else if (SCM_REALP (y
))
4998 double yy
= SCM_REAL_VALUE (y
);
4999 #ifndef ALLOW_DIVIDE_BY_ZERO
5001 scm_num_overflow (s_divide
);
5004 return scm_from_double ((double) xx
/ yy
);
5006 else if (SCM_COMPLEXP (y
))
5009 complex_div
: /* y _must_ be a complex number */
5011 double r
= SCM_COMPLEX_REAL (y
);
5012 double i
= SCM_COMPLEX_IMAG (y
);
5013 if (fabs(r
) <= fabs(i
))
5016 double d
= i
* (1.0 + t
* t
);
5017 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
5022 double d
= r
* (1.0 + t
* t
);
5023 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
5027 else if (SCM_FRACTIONP (y
))
5028 /* a / b/c = ac / b */
5029 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5030 SCM_FRACTION_NUMERATOR (y
));
5032 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5034 else if (SCM_BIGP (x
))
5036 if (SCM_I_INUMP (y
))
5038 long int yy
= SCM_I_INUM (y
);
5041 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5042 scm_num_overflow (s_divide
);
5044 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5045 scm_remember_upto_here_1 (x
);
5046 return (sgn
== 0) ? scm_nan () : scm_inf ();
5053 /* FIXME: HMM, what are the relative performance issues here?
5054 We need to test. Is it faster on average to test
5055 divisible_p, then perform whichever operation, or is it
5056 faster to perform the integer div opportunistically and
5057 switch to real if there's a remainder? For now we take the
5058 middle ground: test, then if divisible, use the faster div
5061 long abs_yy
= yy
< 0 ? -yy
: yy
;
5062 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
5066 SCM result
= scm_i_mkbig ();
5067 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
5068 scm_remember_upto_here_1 (x
);
5070 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5071 return scm_i_normbig (result
);
5076 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
5077 else return scm_i_make_ratio (x
, y
);
5081 else if (SCM_BIGP (y
))
5083 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
5086 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5087 scm_num_overflow (s_divide
);
5089 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5090 scm_remember_upto_here_1 (x
);
5091 return (sgn
== 0) ? scm_nan () : scm_inf ();
5099 /* It's easily possible for the ratio x/y to fit a double
5100 but one or both x and y be too big to fit a double,
5101 hence the use of mpq_get_d rather than converting and
5104 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
5105 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
5106 return scm_from_double (mpq_get_d (q
));
5110 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
5114 SCM result
= scm_i_mkbig ();
5115 mpz_divexact (SCM_I_BIG_MPZ (result
),
5118 scm_remember_upto_here_2 (x
, y
);
5119 return scm_i_normbig (result
);
5122 return scm_i_make_ratio (x
, y
);
5126 else if (SCM_REALP (y
))
5128 double yy
= SCM_REAL_VALUE (y
);
5129 #ifndef ALLOW_DIVIDE_BY_ZERO
5131 scm_num_overflow (s_divide
);
5134 return scm_from_double (scm_i_big2dbl (x
) / yy
);
5136 else if (SCM_COMPLEXP (y
))
5138 a
= scm_i_big2dbl (x
);
5141 else if (SCM_FRACTIONP (y
))
5142 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5143 SCM_FRACTION_NUMERATOR (y
));
5145 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5147 else if (SCM_REALP (x
))
5149 double rx
= SCM_REAL_VALUE (x
);
5150 if (SCM_I_INUMP (y
))
5152 long int yy
= SCM_I_INUM (y
);
5153 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5155 scm_num_overflow (s_divide
);
5158 return scm_from_double (rx
/ (double) yy
);
5160 else if (SCM_BIGP (y
))
5162 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5163 scm_remember_upto_here_1 (y
);
5164 return scm_from_double (rx
/ dby
);
5166 else if (SCM_REALP (y
))
5168 double yy
= SCM_REAL_VALUE (y
);
5169 #ifndef ALLOW_DIVIDE_BY_ZERO
5171 scm_num_overflow (s_divide
);
5174 return scm_from_double (rx
/ yy
);
5176 else if (SCM_COMPLEXP (y
))
5181 else if (SCM_FRACTIONP (y
))
5182 return scm_from_double (rx
/ scm_i_fraction2double (y
));
5184 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5186 else if (SCM_COMPLEXP (x
))
5188 double rx
= SCM_COMPLEX_REAL (x
);
5189 double ix
= SCM_COMPLEX_IMAG (x
);
5190 if (SCM_I_INUMP (y
))
5192 long int yy
= SCM_I_INUM (y
);
5193 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5195 scm_num_overflow (s_divide
);
5200 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
5203 else if (SCM_BIGP (y
))
5205 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5206 scm_remember_upto_here_1 (y
);
5207 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
5209 else if (SCM_REALP (y
))
5211 double yy
= SCM_REAL_VALUE (y
);
5212 #ifndef ALLOW_DIVIDE_BY_ZERO
5214 scm_num_overflow (s_divide
);
5217 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5219 else if (SCM_COMPLEXP (y
))
5221 double ry
= SCM_COMPLEX_REAL (y
);
5222 double iy
= SCM_COMPLEX_IMAG (y
);
5223 if (fabs(ry
) <= fabs(iy
))
5226 double d
= iy
* (1.0 + t
* t
);
5227 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5232 double d
= ry
* (1.0 + t
* t
);
5233 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5236 else if (SCM_FRACTIONP (y
))
5238 double yy
= scm_i_fraction2double (y
);
5239 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5242 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5244 else if (SCM_FRACTIONP (x
))
5246 if (SCM_I_INUMP (y
))
5248 long int yy
= SCM_I_INUM (y
);
5249 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5251 scm_num_overflow (s_divide
);
5254 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5255 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5257 else if (SCM_BIGP (y
))
5259 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5260 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5262 else if (SCM_REALP (y
))
5264 double yy
= SCM_REAL_VALUE (y
);
5265 #ifndef ALLOW_DIVIDE_BY_ZERO
5267 scm_num_overflow (s_divide
);
5270 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5272 else if (SCM_COMPLEXP (y
))
5274 a
= scm_i_fraction2double (x
);
5277 else if (SCM_FRACTIONP (y
))
5278 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5279 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5281 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5284 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5288 scm_divide (SCM x
, SCM y
)
5290 return do_divide (x
, y
, 0);
5293 static SCM
scm_divide2real (SCM x
, SCM y
)
5295 return do_divide (x
, y
, 1);
5301 scm_c_truncate (double x
)
5312 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5313 half-way case (ie. when x is an integer plus 0.5) going upwards.
5314 Then half-way cases are identified and adjusted down if the
5315 round-upwards didn't give the desired even integer.
5317 "plus_half == result" identifies a half-way case. If plus_half, which is
5318 x + 0.5, is an integer then x must be an integer plus 0.5.
5320 An odd "result" value is identified with result/2 != floor(result/2).
5321 This is done with plus_half, since that value is ready for use sooner in
5322 a pipelined cpu, and we're already requiring plus_half == result.
5324 Note however that we need to be careful when x is big and already an
5325 integer. In that case "x+0.5" may round to an adjacent integer, causing
5326 us to return such a value, incorrectly. For instance if the hardware is
5327 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5328 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5329 returned. Or if the hardware is in round-upwards mode, then other bigger
5330 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5331 representable value, 2^128+2^76 (or whatever), again incorrect.
5333 These bad roundings of x+0.5 are avoided by testing at the start whether
5334 x is already an integer. If it is then clearly that's the desired result
5335 already. And if it's not then the exponent must be small enough to allow
5336 an 0.5 to be represented, and hence added without a bad rounding. */
5339 scm_c_round (double x
)
5341 double plus_half
, result
;
5346 plus_half
= x
+ 0.5;
5347 result
= floor (plus_half
);
5348 /* Adjust so that the rounding is towards even. */
5349 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5354 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5356 "Round the number @var{x} towards zero.")
5357 #define FUNC_NAME s_scm_truncate_number
5359 if (scm_is_false (scm_negative_p (x
)))
5360 return scm_floor (x
);
5362 return scm_ceiling (x
);
5366 static SCM exactly_one_half
;
5368 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5370 "Round the number @var{x} towards the nearest integer. "
5371 "When it is exactly halfway between two integers, "
5372 "round towards the even one.")
5373 #define FUNC_NAME s_scm_round_number
5375 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5377 else if (SCM_REALP (x
))
5378 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5381 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5382 single quotient+remainder division then examining to see which way
5383 the rounding should go. */
5384 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5385 SCM result
= scm_floor (plus_half
);
5386 /* Adjust so that the rounding is towards even. */
5387 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5388 && scm_is_true (scm_odd_p (result
)))
5389 return scm_difference (result
, SCM_I_MAKINUM (1));
5396 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5398 "Round the number @var{x} towards minus infinity.")
5399 #define FUNC_NAME s_scm_floor
5401 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5403 else if (SCM_REALP (x
))
5404 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5405 else if (SCM_FRACTIONP (x
))
5407 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5408 SCM_FRACTION_DENOMINATOR (x
));
5409 if (scm_is_false (scm_negative_p (x
)))
5411 /* For positive x, rounding towards zero is correct. */
5416 /* For negative x, we need to return q-1 unless x is an
5417 integer. But fractions are never integer, per our
5419 return scm_difference (q
, SCM_I_MAKINUM (1));
5423 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5427 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5429 "Round the number @var{x} towards infinity.")
5430 #define FUNC_NAME s_scm_ceiling
5432 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5434 else if (SCM_REALP (x
))
5435 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5436 else if (SCM_FRACTIONP (x
))
5438 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5439 SCM_FRACTION_DENOMINATOR (x
));
5440 if (scm_is_false (scm_positive_p (x
)))
5442 /* For negative x, rounding towards zero is correct. */
5447 /* For positive x, we need to return q+1 unless x is an
5448 integer. But fractions are never integer, per our
5450 return scm_sum (q
, SCM_I_MAKINUM (1));
5454 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5458 /* sin/cos/tan/asin/acos/atan
5459 sinh/cosh/tanh/asinh/acosh/atanh
5460 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5461 Written by Jerry D. Hedden, (C) FSF.
5462 See the file `COPYING' for terms applying to this program. */
5464 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5466 "Return @var{x} raised to the power of @var{y}.")
5467 #define FUNC_NAME s_scm_expt
5469 if ((SCM_I_INUMP (x
) || SCM_BIGP (x
)) && scm_is_integer (y
))
5470 return scm_integer_expt (x
, y
);
5471 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5473 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5476 return scm_exp (scm_product (scm_log (x
), y
));
5480 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5482 "Compute the sine of @var{z}.")
5483 #define FUNC_NAME s_scm_sin
5485 if (scm_is_real (z
))
5486 return scm_from_double (sin (scm_to_double (z
)));
5487 else if (SCM_COMPLEXP (z
))
5489 x
= SCM_COMPLEX_REAL (z
);
5490 y
= SCM_COMPLEX_IMAG (z
);
5491 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5492 cos (x
) * sinh (y
));
5495 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5499 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5501 "Compute the cosine of @var{z}.")
5502 #define FUNC_NAME s_scm_cos
5504 if (scm_is_real (z
))
5505 return scm_from_double (cos (scm_to_double (z
)));
5506 else if (SCM_COMPLEXP (z
))
5508 x
= SCM_COMPLEX_REAL (z
);
5509 y
= SCM_COMPLEX_IMAG (z
);
5510 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5511 -sin (x
) * sinh (y
));
5514 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5518 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5520 "Compute the tangent of @var{z}.")
5521 #define FUNC_NAME s_scm_tan
5523 if (scm_is_real (z
))
5524 return scm_from_double (tan (scm_to_double (z
)));
5525 else if (SCM_COMPLEXP (z
))
5527 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5528 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5529 w
= cos (x
) + cosh (y
);
5530 #ifndef ALLOW_DIVIDE_BY_ZERO
5532 scm_num_overflow (s_scm_tan
);
5534 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5537 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5541 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5543 "Compute the hyperbolic sine of @var{z}.")
5544 #define FUNC_NAME s_scm_sinh
5546 if (scm_is_real (z
))
5547 return scm_from_double (sinh (scm_to_double (z
)));
5548 else if (SCM_COMPLEXP (z
))
5550 x
= SCM_COMPLEX_REAL (z
);
5551 y
= SCM_COMPLEX_IMAG (z
);
5552 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5553 cosh (x
) * sin (y
));
5556 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5560 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5562 "Compute the hyperbolic cosine of @var{z}.")
5563 #define FUNC_NAME s_scm_cosh
5565 if (scm_is_real (z
))
5566 return scm_from_double (cosh (scm_to_double (z
)));
5567 else if (SCM_COMPLEXP (z
))
5569 x
= SCM_COMPLEX_REAL (z
);
5570 y
= SCM_COMPLEX_IMAG (z
);
5571 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5572 sinh (x
) * sin (y
));
5575 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5579 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5581 "Compute the hyperbolic tangent of @var{z}.")
5582 #define FUNC_NAME s_scm_tanh
5584 if (scm_is_real (z
))
5585 return scm_from_double (tanh (scm_to_double (z
)));
5586 else if (SCM_COMPLEXP (z
))
5588 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5589 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5590 w
= cosh (x
) + cos (y
);
5591 #ifndef ALLOW_DIVIDE_BY_ZERO
5593 scm_num_overflow (s_scm_tanh
);
5595 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5598 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5602 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5604 "Compute the arc sine of @var{z}.")
5605 #define FUNC_NAME s_scm_asin
5607 if (scm_is_real (z
))
5609 double w
= scm_to_double (z
);
5610 if (w
>= -1.0 && w
<= 1.0)
5611 return scm_from_double (asin (w
));
5613 return scm_product (scm_c_make_rectangular (0, -1),
5614 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5616 else if (SCM_COMPLEXP (z
))
5618 x
= SCM_COMPLEX_REAL (z
);
5619 y
= SCM_COMPLEX_IMAG (z
);
5620 return scm_product (scm_c_make_rectangular (0, -1),
5621 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5624 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5628 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5630 "Compute the arc cosine of @var{z}.")
5631 #define FUNC_NAME s_scm_acos
5633 if (scm_is_real (z
))
5635 double w
= scm_to_double (z
);
5636 if (w
>= -1.0 && w
<= 1.0)
5637 return scm_from_double (acos (w
));
5639 return scm_sum (scm_from_double (acos (0.0)),
5640 scm_product (scm_c_make_rectangular (0, 1),
5641 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5643 else if (SCM_COMPLEXP (z
))
5645 x
= SCM_COMPLEX_REAL (z
);
5646 y
= SCM_COMPLEX_IMAG (z
);
5647 return scm_sum (scm_from_double (acos (0.0)),
5648 scm_product (scm_c_make_rectangular (0, 1),
5649 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5652 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5656 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5658 "With one argument, compute the arc tangent of @var{z}.\n"
5659 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5660 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5661 #define FUNC_NAME s_scm_atan
5665 if (scm_is_real (z
))
5666 return scm_from_double (atan (scm_to_double (z
)));
5667 else if (SCM_COMPLEXP (z
))
5670 v
= SCM_COMPLEX_REAL (z
);
5671 w
= SCM_COMPLEX_IMAG (z
);
5672 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5673 scm_c_make_rectangular (v
, w
+ 1.0))),
5674 scm_c_make_rectangular (0, 2));
5677 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5679 else if (scm_is_real (z
))
5681 if (scm_is_real (y
))
5682 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5684 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5687 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5691 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5693 "Compute the inverse hyperbolic sine of @var{z}.")
5694 #define FUNC_NAME s_scm_sys_asinh
5696 if (scm_is_real (z
))
5697 return scm_from_double (asinh (scm_to_double (z
)));
5698 else if (scm_is_number (z
))
5699 return scm_log (scm_sum (z
,
5700 scm_sqrt (scm_sum (scm_product (z
, z
),
5701 SCM_I_MAKINUM (1)))));
5703 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5707 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5709 "Compute the inverse hyperbolic cosine of @var{z}.")
5710 #define FUNC_NAME s_scm_sys_acosh
5712 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5713 return scm_from_double (acosh (scm_to_double (z
)));
5714 else if (scm_is_number (z
))
5715 return scm_log (scm_sum (z
,
5716 scm_sqrt (scm_difference (scm_product (z
, z
),
5717 SCM_I_MAKINUM (1)))));
5719 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5723 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5725 "Compute the inverse hyperbolic tangent of @var{z}.")
5726 #define FUNC_NAME s_scm_sys_atanh
5728 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5729 return scm_from_double (atanh (scm_to_double (z
)));
5730 else if (scm_is_number (z
))
5731 return scm_divide (scm_log (scm_divide (scm_sum (SCM_I_MAKINUM (1), z
),
5732 scm_difference (SCM_I_MAKINUM (1), z
))),
5735 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5740 scm_c_make_rectangular (double re
, double im
)
5743 return scm_from_double (re
);
5747 SCM_NEWSMOB (z
, scm_tc16_complex
,
5748 scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5750 SCM_COMPLEX_REAL (z
) = re
;
5751 SCM_COMPLEX_IMAG (z
) = im
;
5756 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5757 (SCM real_part
, SCM imaginary_part
),
5758 "Return a complex number constructed of the given @var{real-part} "
5759 "and @var{imaginary-part} parts.")
5760 #define FUNC_NAME s_scm_make_rectangular
5762 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5763 SCM_ARG1
, FUNC_NAME
, "real");
5764 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5765 SCM_ARG2
, FUNC_NAME
, "real");
5766 return scm_c_make_rectangular (scm_to_double (real_part
),
5767 scm_to_double (imaginary_part
));
5772 scm_c_make_polar (double mag
, double ang
)
5776 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5777 use it on Glibc-based systems that have it (it's a GNU extension). See
5778 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5780 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5781 sincos (ang
, &s
, &c
);
5786 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5789 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5791 "Return the complex number @var{x} * e^(i * @var{y}).")
5792 #define FUNC_NAME s_scm_make_polar
5794 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5795 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5796 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5801 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5802 /* "Return the real part of the number @var{z}."
5805 scm_real_part (SCM z
)
5807 if (SCM_I_INUMP (z
))
5809 else if (SCM_BIGP (z
))
5811 else if (SCM_REALP (z
))
5813 else if (SCM_COMPLEXP (z
))
5814 return scm_from_double (SCM_COMPLEX_REAL (z
));
5815 else if (SCM_FRACTIONP (z
))
5818 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5822 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5823 /* "Return the imaginary part of the number @var{z}."
5826 scm_imag_part (SCM z
)
5828 if (SCM_I_INUMP (z
))
5830 else if (SCM_BIGP (z
))
5832 else if (SCM_REALP (z
))
5834 else if (SCM_COMPLEXP (z
))
5835 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5836 else if (SCM_FRACTIONP (z
))
5839 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5842 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5843 /* "Return the numerator of the number @var{z}."
5846 scm_numerator (SCM z
)
5848 if (SCM_I_INUMP (z
))
5850 else if (SCM_BIGP (z
))
5852 else if (SCM_FRACTIONP (z
))
5853 return SCM_FRACTION_NUMERATOR (z
);
5854 else if (SCM_REALP (z
))
5855 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5857 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5861 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5862 /* "Return the denominator of the number @var{z}."
5865 scm_denominator (SCM z
)
5867 if (SCM_I_INUMP (z
))
5868 return SCM_I_MAKINUM (1);
5869 else if (SCM_BIGP (z
))
5870 return SCM_I_MAKINUM (1);
5871 else if (SCM_FRACTIONP (z
))
5872 return SCM_FRACTION_DENOMINATOR (z
);
5873 else if (SCM_REALP (z
))
5874 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5876 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5879 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5880 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5881 * "@code{abs} for real arguments, but also allows complex numbers."
5884 scm_magnitude (SCM z
)
5886 if (SCM_I_INUMP (z
))
5888 long int zz
= SCM_I_INUM (z
);
5891 else if (SCM_POSFIXABLE (-zz
))
5892 return SCM_I_MAKINUM (-zz
);
5894 return scm_i_long2big (-zz
);
5896 else if (SCM_BIGP (z
))
5898 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5899 scm_remember_upto_here_1 (z
);
5901 return scm_i_clonebig (z
, 0);
5905 else if (SCM_REALP (z
))
5906 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5907 else if (SCM_COMPLEXP (z
))
5908 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5909 else if (SCM_FRACTIONP (z
))
5911 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5913 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5914 SCM_FRACTION_DENOMINATOR (z
));
5917 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5921 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5922 /* "Return the angle of the complex number @var{z}."
5927 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5928 flo0 to save allocating a new flonum with scm_from_double each time.
5929 But if atan2 follows the floating point rounding mode, then the value
5930 is not a constant. Maybe it'd be close enough though. */
5931 if (SCM_I_INUMP (z
))
5933 if (SCM_I_INUM (z
) >= 0)
5936 return scm_from_double (atan2 (0.0, -1.0));
5938 else if (SCM_BIGP (z
))
5940 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5941 scm_remember_upto_here_1 (z
);
5943 return scm_from_double (atan2 (0.0, -1.0));
5947 else if (SCM_REALP (z
))
5949 if (SCM_REAL_VALUE (z
) >= 0)
5952 return scm_from_double (atan2 (0.0, -1.0));
5954 else if (SCM_COMPLEXP (z
))
5955 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5956 else if (SCM_FRACTIONP (z
))
5958 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5960 else return scm_from_double (atan2 (0.0, -1.0));
5963 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5967 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5968 /* Convert the number @var{x} to its inexact representation.\n"
5971 scm_exact_to_inexact (SCM z
)
5973 if (SCM_I_INUMP (z
))
5974 return scm_from_double ((double) SCM_I_INUM (z
));
5975 else if (SCM_BIGP (z
))
5976 return scm_from_double (scm_i_big2dbl (z
));
5977 else if (SCM_FRACTIONP (z
))
5978 return scm_from_double (scm_i_fraction2double (z
));
5979 else if (SCM_INEXACTP (z
))
5982 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5986 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5988 "Return an exact number that is numerically closest to @var{z}.")
5989 #define FUNC_NAME s_scm_inexact_to_exact
5991 if (SCM_I_INUMP (z
))
5993 else if (SCM_BIGP (z
))
5995 else if (SCM_REALP (z
))
5997 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
5998 SCM_OUT_OF_RANGE (1, z
);
6005 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
6006 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
6007 scm_i_mpz2num (mpq_denref (frac
)));
6009 /* When scm_i_make_ratio throws, we leak the memory allocated
6016 else if (SCM_FRACTIONP (z
))
6019 SCM_WRONG_TYPE_ARG (1, z
);
6023 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
6025 "Returns the @emph{simplest} rational number differing\n"
6026 "from @var{x} by no more than @var{eps}.\n"
6028 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
6029 "exact result when both its arguments are exact. Thus, you might need\n"
6030 "to use @code{inexact->exact} on the arguments.\n"
6033 "(rationalize (inexact->exact 1.2) 1/100)\n"
6036 #define FUNC_NAME s_scm_rationalize
6038 if (SCM_I_INUMP (x
))
6040 else if (SCM_BIGP (x
))
6042 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
6044 /* Use continued fractions to find closest ratio. All
6045 arithmetic is done with exact numbers.
6048 SCM ex
= scm_inexact_to_exact (x
);
6049 SCM int_part
= scm_floor (ex
);
6050 SCM tt
= SCM_I_MAKINUM (1);
6051 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
6052 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
6056 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
6059 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
6060 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
6062 /* We stop after a million iterations just to be absolutely sure
6063 that we don't go into an infinite loop. The process normally
6064 converges after less than a dozen iterations.
6067 eps
= scm_abs (eps
);
6068 while (++i
< 1000000)
6070 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
6071 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
6072 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
6074 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
6075 eps
))) /* abs(x-a/b) <= eps */
6077 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
6078 if (scm_is_false (scm_exact_p (x
))
6079 || scm_is_false (scm_exact_p (eps
)))
6080 return scm_exact_to_inexact (res
);
6084 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
6086 tt
= scm_floor (rx
); /* tt = floor (rx) */
6092 scm_num_overflow (s_scm_rationalize
);
6095 SCM_WRONG_TYPE_ARG (1, x
);
6099 /* conversion functions */
6102 scm_is_integer (SCM val
)
6104 return scm_is_true (scm_integer_p (val
));
6108 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
6110 if (SCM_I_INUMP (val
))
6112 scm_t_signed_bits n
= SCM_I_INUM (val
);
6113 return n
>= min
&& n
<= max
;
6115 else if (SCM_BIGP (val
))
6117 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
6119 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
6121 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
6123 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
6124 return n
>= min
&& n
<= max
;
6134 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6135 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6138 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6139 SCM_I_BIG_MPZ (val
));
6141 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
6153 return n
>= min
&& n
<= max
;
6161 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
6163 if (SCM_I_INUMP (val
))
6165 scm_t_signed_bits n
= SCM_I_INUM (val
);
6166 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
6168 else if (SCM_BIGP (val
))
6170 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
6172 else if (max
<= ULONG_MAX
)
6174 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
6176 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
6177 return n
>= min
&& n
<= max
;
6187 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
6190 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6191 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6194 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6195 SCM_I_BIG_MPZ (val
));
6197 return n
>= min
&& n
<= max
;
6205 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
6207 scm_error (scm_out_of_range_key
,
6209 "Value out of range ~S to ~S: ~S",
6210 scm_list_3 (min
, max
, bad_val
),
6211 scm_list_1 (bad_val
));
6214 #define TYPE scm_t_intmax
6215 #define TYPE_MIN min
6216 #define TYPE_MAX max
6217 #define SIZEOF_TYPE 0
6218 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
6219 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
6220 #include "libguile/conv-integer.i.c"
6222 #define TYPE scm_t_uintmax
6223 #define TYPE_MIN min
6224 #define TYPE_MAX max
6225 #define SIZEOF_TYPE 0
6226 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6227 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6228 #include "libguile/conv-uinteger.i.c"
6230 #define TYPE scm_t_int8
6231 #define TYPE_MIN SCM_T_INT8_MIN
6232 #define TYPE_MAX SCM_T_INT8_MAX
6233 #define SIZEOF_TYPE 1
6234 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6235 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6236 #include "libguile/conv-integer.i.c"
6238 #define TYPE scm_t_uint8
6240 #define TYPE_MAX SCM_T_UINT8_MAX
6241 #define SIZEOF_TYPE 1
6242 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6243 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6244 #include "libguile/conv-uinteger.i.c"
6246 #define TYPE scm_t_int16
6247 #define TYPE_MIN SCM_T_INT16_MIN
6248 #define TYPE_MAX SCM_T_INT16_MAX
6249 #define SIZEOF_TYPE 2
6250 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6251 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6252 #include "libguile/conv-integer.i.c"
6254 #define TYPE scm_t_uint16
6256 #define TYPE_MAX SCM_T_UINT16_MAX
6257 #define SIZEOF_TYPE 2
6258 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6259 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6260 #include "libguile/conv-uinteger.i.c"
6262 #define TYPE scm_t_int32
6263 #define TYPE_MIN SCM_T_INT32_MIN
6264 #define TYPE_MAX SCM_T_INT32_MAX
6265 #define SIZEOF_TYPE 4
6266 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6267 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6268 #include "libguile/conv-integer.i.c"
6270 #define TYPE scm_t_uint32
6272 #define TYPE_MAX SCM_T_UINT32_MAX
6273 #define SIZEOF_TYPE 4
6274 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6275 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6276 #include "libguile/conv-uinteger.i.c"
6278 #define TYPE scm_t_wchar
6279 #define TYPE_MIN (scm_t_int32)-1
6280 #define TYPE_MAX (scm_t_int32)0x10ffff
6281 #define SIZEOF_TYPE 4
6282 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6283 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6284 #include "libguile/conv-integer.i.c"
6286 #if SCM_HAVE_T_INT64
6288 #define TYPE scm_t_int64
6289 #define TYPE_MIN SCM_T_INT64_MIN
6290 #define TYPE_MAX SCM_T_INT64_MAX
6291 #define SIZEOF_TYPE 8
6292 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6293 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6294 #include "libguile/conv-integer.i.c"
6296 #define TYPE scm_t_uint64
6298 #define TYPE_MAX SCM_T_UINT64_MAX
6299 #define SIZEOF_TYPE 8
6300 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6301 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6302 #include "libguile/conv-uinteger.i.c"
6307 scm_to_mpz (SCM val
, mpz_t rop
)
6309 if (SCM_I_INUMP (val
))
6310 mpz_set_si (rop
, SCM_I_INUM (val
));
6311 else if (SCM_BIGP (val
))
6312 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6314 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6318 scm_from_mpz (mpz_t val
)
6320 return scm_i_mpz2num (val
);
6324 scm_is_real (SCM val
)
6326 return scm_is_true (scm_real_p (val
));
6330 scm_is_rational (SCM val
)
6332 return scm_is_true (scm_rational_p (val
));
6336 scm_to_double (SCM val
)
6338 if (SCM_I_INUMP (val
))
6339 return SCM_I_INUM (val
);
6340 else if (SCM_BIGP (val
))
6341 return scm_i_big2dbl (val
);
6342 else if (SCM_FRACTIONP (val
))
6343 return scm_i_fraction2double (val
);
6344 else if (SCM_REALP (val
))
6345 return SCM_REAL_VALUE (val
);
6347 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6351 scm_from_double (double val
)
6353 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
6354 SCM_REAL_VALUE (z
) = val
;
6358 #if SCM_ENABLE_DISCOURAGED == 1
6361 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
6365 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6369 scm_out_of_range (NULL
, num
);
6372 return scm_to_double (num
);
6376 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
6380 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6384 scm_out_of_range (NULL
, num
);
6387 return scm_to_double (num
);
6393 scm_is_complex (SCM val
)
6395 return scm_is_true (scm_complex_p (val
));
6399 scm_c_real_part (SCM z
)
6401 if (SCM_COMPLEXP (z
))
6402 return SCM_COMPLEX_REAL (z
);
6405 /* Use the scm_real_part to get proper error checking and
6408 return scm_to_double (scm_real_part (z
));
6413 scm_c_imag_part (SCM z
)
6415 if (SCM_COMPLEXP (z
))
6416 return SCM_COMPLEX_IMAG (z
);
6419 /* Use the scm_imag_part to get proper error checking and
6420 dispatching. The result will almost always be 0.0, but not
6423 return scm_to_double (scm_imag_part (z
));
6428 scm_c_magnitude (SCM z
)
6430 return scm_to_double (scm_magnitude (z
));
6436 return scm_to_double (scm_angle (z
));
6440 scm_is_number (SCM z
)
6442 return scm_is_true (scm_number_p (z
));
6446 /* In the following functions we dispatch to the real-arg funcs like log()
6447 when we know the arg is real, instead of just handing everything to
6448 clog() for instance. This is in case clog() doesn't optimize for a
6449 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6450 well use it to go straight to the applicable C func. */
6452 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6454 "Return the natural logarithm of @var{z}.")
6455 #define FUNC_NAME s_scm_log
6457 if (SCM_COMPLEXP (z
))
6459 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6460 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6462 double re
= SCM_COMPLEX_REAL (z
);
6463 double im
= SCM_COMPLEX_IMAG (z
);
6464 return scm_c_make_rectangular (log (hypot (re
, im
)),
6470 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6471 although the value itself overflows. */
6472 double re
= scm_to_double (z
);
6473 double l
= log (fabs (re
));
6475 return scm_from_double (l
);
6477 return scm_c_make_rectangular (l
, M_PI
);
6483 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6485 "Return the base 10 logarithm of @var{z}.")
6486 #define FUNC_NAME s_scm_log10
6488 if (SCM_COMPLEXP (z
))
6490 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6491 clog() and a multiply by M_LOG10E, rather than the fallback
6492 log10+hypot+atan2.) */
6493 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6494 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6496 double re
= SCM_COMPLEX_REAL (z
);
6497 double im
= SCM_COMPLEX_IMAG (z
);
6498 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6499 M_LOG10E
* atan2 (im
, re
));
6504 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6505 although the value itself overflows. */
6506 double re
= scm_to_double (z
);
6507 double l
= log10 (fabs (re
));
6509 return scm_from_double (l
);
6511 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6517 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6519 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6520 "base of natural logarithms (2.71828@dots{}).")
6521 #define FUNC_NAME s_scm_exp
6523 if (SCM_COMPLEXP (z
))
6525 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6526 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6528 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6529 SCM_COMPLEX_IMAG (z
));
6534 /* When z is a negative bignum the conversion to double overflows,
6535 giving -infinity, but that's ok, the exp is still 0.0. */
6536 return scm_from_double (exp (scm_to_double (z
)));
6542 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6544 "Return the square root of @var{z}. Of the two possible roots\n"
6545 "(positive and negative), the one with the a positive real part\n"
6546 "is returned, or if that's zero then a positive imaginary part.\n"
6550 "(sqrt 9.0) @result{} 3.0\n"
6551 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6552 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6553 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6555 #define FUNC_NAME s_scm_sqrt
6557 if (SCM_COMPLEXP (x
))
6559 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6560 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6562 double re
= SCM_COMPLEX_REAL (x
);
6563 double im
= SCM_COMPLEX_IMAG (x
);
6564 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6565 0.5 * atan2 (im
, re
));
6570 double xx
= scm_to_double (x
);
6572 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6574 return scm_from_double (sqrt (xx
));
6586 mpz_init_set_si (z_negative_one
, -1);
6588 /* It may be possible to tune the performance of some algorithms by using
6589 * the following constants to avoid the creation of bignums. Please, before
6590 * using these values, remember the two rules of program optimization:
6591 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6592 scm_c_define ("most-positive-fixnum",
6593 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6594 scm_c_define ("most-negative-fixnum",
6595 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6597 scm_add_feature ("complex");
6598 scm_add_feature ("inexact");
6599 flo0
= scm_from_double (0.0);
6601 /* determine floating point precision */
6602 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6604 init_dblprec(&scm_dblprec
[i
-2],i
);
6605 init_fx_radix(fx_per_radix
[i
-2],i
);
6608 /* hard code precision for base 10 if the preprocessor tells us to... */
6609 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6612 exactly_one_half
= scm_divide (SCM_I_MAKINUM (1), SCM_I_MAKINUM (2));
6613 #include "libguile/numbers.x"