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 SCM_VALIDATE_NUMBER (SCM_ARG1
, n
);
1807 /* 0^0 == 1 according to R5RS */
1808 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1809 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1810 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1811 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1813 if (SCM_I_INUMP (k
))
1814 i2
= SCM_I_INUM (k
);
1815 else if (SCM_BIGP (k
))
1817 z_i2
= scm_i_clonebig (k
, 1);
1818 scm_remember_upto_here_1 (k
);
1822 SCM_WRONG_TYPE_ARG (2, k
);
1826 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1828 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1829 n
= scm_divide (n
, SCM_UNDEFINED
);
1833 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1837 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1839 return scm_product (acc
, n
);
1841 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1842 acc
= scm_product (acc
, n
);
1843 n
= scm_product (n
, n
);
1844 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1852 n
= scm_divide (n
, SCM_UNDEFINED
);
1859 return scm_product (acc
, n
);
1861 acc
= scm_product (acc
, n
);
1862 n
= scm_product (n
, n
);
1869 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1871 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1872 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1874 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1875 "@var{cnt} is negative it's a division, rounded towards negative\n"
1876 "infinity. (Note that this is not the same rounding as\n"
1877 "@code{quotient} does.)\n"
1879 "With @var{n} viewed as an infinite precision twos complement,\n"
1880 "@code{ash} means a left shift introducing zero bits, or a right\n"
1881 "shift dropping bits.\n"
1884 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1885 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1887 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1888 "(ash -23 -2) @result{} -6\n"
1890 #define FUNC_NAME s_scm_ash
1893 bits_to_shift
= scm_to_long (cnt
);
1895 if (SCM_I_INUMP (n
))
1897 long nn
= SCM_I_INUM (n
);
1899 if (bits_to_shift
> 0)
1901 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1902 overflow a non-zero fixnum. For smaller shifts we check the
1903 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1904 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1905 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1911 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1913 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1916 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1920 SCM result
= scm_i_long2big (nn
);
1921 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1928 bits_to_shift
= -bits_to_shift
;
1929 if (bits_to_shift
>= SCM_LONG_BIT
)
1930 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1932 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1936 else if (SCM_BIGP (n
))
1940 if (bits_to_shift
== 0)
1943 result
= scm_i_mkbig ();
1944 if (bits_to_shift
>= 0)
1946 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1952 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1953 we have to allocate a bignum even if the result is going to be a
1955 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1957 return scm_i_normbig (result
);
1963 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1969 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1970 (SCM n
, SCM start
, SCM end
),
1971 "Return the integer composed of the @var{start} (inclusive)\n"
1972 "through @var{end} (exclusive) bits of @var{n}. The\n"
1973 "@var{start}th bit becomes the 0-th bit in the result.\n"
1976 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1977 " @result{} \"1010\"\n"
1978 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1979 " @result{} \"10110\"\n"
1981 #define FUNC_NAME s_scm_bit_extract
1983 unsigned long int istart
, iend
, bits
;
1984 istart
= scm_to_ulong (start
);
1985 iend
= scm_to_ulong (end
);
1986 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1988 /* how many bits to keep */
1989 bits
= iend
- istart
;
1991 if (SCM_I_INUMP (n
))
1993 long int in
= SCM_I_INUM (n
);
1995 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1996 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1997 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1999 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
2001 /* Since we emulate two's complement encoded numbers, this
2002 * special case requires us to produce a result that has
2003 * more bits than can be stored in a fixnum.
2005 SCM result
= scm_i_long2big (in
);
2006 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
2011 /* mask down to requisite bits */
2012 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
2013 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
2015 else if (SCM_BIGP (n
))
2020 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
2024 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
2025 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
2026 such bits into a ulong. */
2027 result
= scm_i_mkbig ();
2028 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
2029 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
2030 result
= scm_i_normbig (result
);
2032 scm_remember_upto_here_1 (n
);
2036 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2041 static const char scm_logtab
[] = {
2042 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
2045 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
2047 "Return the number of bits in integer @var{n}. If integer is\n"
2048 "positive, the 1-bits in its binary representation are counted.\n"
2049 "If negative, the 0-bits in its two's-complement binary\n"
2050 "representation are counted. If 0, 0 is returned.\n"
2053 "(logcount #b10101010)\n"
2060 #define FUNC_NAME s_scm_logcount
2062 if (SCM_I_INUMP (n
))
2064 unsigned long int c
= 0;
2065 long int nn
= SCM_I_INUM (n
);
2070 c
+= scm_logtab
[15 & nn
];
2073 return SCM_I_MAKINUM (c
);
2075 else if (SCM_BIGP (n
))
2077 unsigned long count
;
2078 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
2079 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
2081 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
2082 scm_remember_upto_here_1 (n
);
2083 return SCM_I_MAKINUM (count
);
2086 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2091 static const char scm_ilentab
[] = {
2092 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2096 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2098 "Return the number of bits necessary to represent @var{n}.\n"
2101 "(integer-length #b10101010)\n"
2103 "(integer-length 0)\n"
2105 "(integer-length #b1111)\n"
2108 #define FUNC_NAME s_scm_integer_length
2110 if (SCM_I_INUMP (n
))
2112 unsigned long int c
= 0;
2114 long int nn
= SCM_I_INUM (n
);
2120 l
= scm_ilentab
[15 & nn
];
2123 return SCM_I_MAKINUM (c
- 4 + l
);
2125 else if (SCM_BIGP (n
))
2127 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2128 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2129 1 too big, so check for that and adjust. */
2130 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2131 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2132 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2133 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2135 scm_remember_upto_here_1 (n
);
2136 return SCM_I_MAKINUM (size
);
2139 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2143 /*** NUMBERS -> STRINGS ***/
2144 #define SCM_MAX_DBL_PREC 60
2145 #define SCM_MAX_DBL_RADIX 36
2147 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2148 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2149 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2152 void init_dblprec(int *prec
, int radix
) {
2153 /* determine floating point precision by adding successively
2154 smaller increments to 1.0 until it is considered == 1.0 */
2155 double f
= ((double)1.0)/radix
;
2156 double fsum
= 1.0 + f
;
2161 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2173 void init_fx_radix(double *fx_list
, int radix
)
2175 /* initialize a per-radix list of tolerances. When added
2176 to a number < 1.0, we can determine if we should raund
2177 up and quit converting a number to a string. */
2181 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2182 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2185 /* use this array as a way to generate a single digit */
2186 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2189 idbl2str (double f
, char *a
, int radix
)
2191 int efmt
, dpt
, d
, i
, wp
;
2193 #ifdef DBL_MIN_10_EXP
2196 #endif /* DBL_MIN_10_EXP */
2201 radix
> SCM_MAX_DBL_RADIX
)
2203 /* revert to existing behavior */
2207 wp
= scm_dblprec
[radix
-2];
2208 fx
= fx_per_radix
[radix
-2];
2212 #ifdef HAVE_COPYSIGN
2213 double sgn
= copysign (1.0, f
);
2218 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2224 strcpy (a
, "-inf.0");
2226 strcpy (a
, "+inf.0");
2229 else if (xisnan (f
))
2231 strcpy (a
, "+nan.0");
2241 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2242 make-uniform-vector, from causing infinite loops. */
2243 /* just do the checking...if it passes, we do the conversion for our
2244 radix again below */
2251 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2259 while (f_cpy
> 10.0)
2262 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2283 if (f
+ fx
[wp
] >= radix
)
2290 /* adding 9999 makes this equivalent to abs(x) % 3 */
2291 dpt
= (exp
+ 9999) % 3;
2295 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2317 a
[ch
++] = number_chars
[d
];
2320 if (f
+ fx
[wp
] >= 1.0)
2322 a
[ch
- 1] = number_chars
[d
+1];
2334 if ((dpt
> 4) && (exp
> 6))
2336 d
= (a
[0] == '-' ? 2 : 1);
2337 for (i
= ch
++; i
> d
; i
--)
2350 if (a
[ch
- 1] == '.')
2351 a
[ch
++] = '0'; /* trailing zero */
2360 for (i
= radix
; i
<= exp
; i
*= radix
);
2361 for (i
/= radix
; i
; i
/= radix
)
2363 a
[ch
++] = number_chars
[exp
/ i
];
2372 icmplx2str (double real
, double imag
, char *str
, int radix
)
2376 i
= idbl2str (real
, str
, radix
);
2379 /* Don't output a '+' for negative numbers or for Inf and
2380 NaN. They will provide their own sign. */
2381 if (0 <= imag
&& !xisinf (imag
) && !xisnan (imag
))
2383 i
+= idbl2str (imag
, &str
[i
], radix
);
2390 iflo2str (SCM flt
, char *str
, int radix
)
2393 if (SCM_REALP (flt
))
2394 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2396 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2401 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2402 characters in the result.
2404 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2406 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2411 return scm_iuint2str (-num
, rad
, p
) + 1;
2414 return scm_iuint2str (num
, rad
, p
);
2417 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2418 characters in the result.
2420 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2422 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2426 scm_t_uintmax n
= num
;
2428 for (n
/= rad
; n
> 0; n
/= rad
)
2438 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2443 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2445 "Return a string holding the external representation of the\n"
2446 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2447 "inexact, a radix of 10 will be used.")
2448 #define FUNC_NAME s_scm_number_to_string
2452 if (SCM_UNBNDP (radix
))
2455 base
= scm_to_signed_integer (radix
, 2, 36);
2457 if (SCM_I_INUMP (n
))
2459 char num_buf
[SCM_INTBUFLEN
];
2460 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2461 return scm_from_locale_stringn (num_buf
, length
);
2463 else if (SCM_BIGP (n
))
2465 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2466 scm_remember_upto_here_1 (n
);
2467 return scm_take_locale_string (str
);
2469 else if (SCM_FRACTIONP (n
))
2471 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2472 scm_from_locale_string ("/"),
2473 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2475 else if (SCM_INEXACTP (n
))
2477 char num_buf
[FLOBUFLEN
];
2478 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2481 SCM_WRONG_TYPE_ARG (1, n
);
2486 /* These print routines used to be stubbed here so that scm_repl.c
2487 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2490 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2492 char num_buf
[FLOBUFLEN
];
2493 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2498 scm_i_print_double (double val
, SCM port
)
2500 char num_buf
[FLOBUFLEN
];
2501 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2505 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2508 char num_buf
[FLOBUFLEN
];
2509 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2514 scm_i_print_complex (double real
, double imag
, SCM port
)
2516 char num_buf
[FLOBUFLEN
];
2517 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2521 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2524 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2525 scm_lfwrite_str (str
, port
);
2526 scm_remember_upto_here_1 (str
);
2531 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2533 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2534 scm_remember_upto_here_1 (exp
);
2535 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2539 /*** END nums->strs ***/
2542 /*** STRINGS -> NUMBERS ***/
2544 /* The following functions implement the conversion from strings to numbers.
2545 * The implementation somehow follows the grammar for numbers as it is given
2546 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2547 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2548 * points should be noted about the implementation:
2549 * * Each function keeps a local index variable 'idx' that points at the
2550 * current position within the parsed string. The global index is only
2551 * updated if the function could parse the corresponding syntactic unit
2553 * * Similarly, the functions keep track of indicators of inexactness ('#',
2554 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2555 * global exactness information is only updated after each part has been
2556 * successfully parsed.
2557 * * Sequences of digits are parsed into temporary variables holding fixnums.
2558 * Only if these fixnums would overflow, the result variables are updated
2559 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2560 * the temporary variables holding the fixnums are cleared, and the process
2561 * starts over again. If for example fixnums were able to store five decimal
2562 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2563 * and the result was computed as 12345 * 100000 + 67890. In other words,
2564 * only every five digits two bignum operations were performed.
2567 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2569 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2571 /* In non ASCII-style encodings the following macro might not work. */
2572 #define XDIGIT2UINT(d) \
2573 (uc_is_property_decimal_digit ((int) (unsigned char) d) \
2575 : uc_tolower ((int) (unsigned char) d) - 'a' + 10)
2578 mem2uinteger (SCM mem
, unsigned int *p_idx
,
2579 unsigned int radix
, enum t_exactness
*p_exactness
)
2581 unsigned int idx
= *p_idx
;
2582 unsigned int hash_seen
= 0;
2583 scm_t_bits shift
= 1;
2585 unsigned int digit_value
;
2588 size_t len
= scm_i_string_length (mem
);
2593 c
= scm_i_string_ref (mem
, idx
);
2594 if (!uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2596 digit_value
= XDIGIT2UINT (c
);
2597 if (digit_value
>= radix
)
2601 result
= SCM_I_MAKINUM (digit_value
);
2604 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2605 if (uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2609 digit_value
= XDIGIT2UINT (c
);
2610 if (digit_value
>= radix
)
2622 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2624 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2626 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2633 shift
= shift
* radix
;
2634 add
= add
* radix
+ digit_value
;
2639 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2641 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2645 *p_exactness
= INEXACT
;
2651 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2652 * covers the parts of the rules that start at a potential point. The value
2653 * of the digits up to the point have been parsed by the caller and are given
2654 * in variable result. The content of *p_exactness indicates, whether a hash
2655 * has already been seen in the digits before the point.
2658 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
2661 mem2decimal_from_point (SCM result
, SCM mem
,
2662 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2664 unsigned int idx
= *p_idx
;
2665 enum t_exactness x
= *p_exactness
;
2666 size_t len
= scm_i_string_length (mem
);
2671 if (scm_i_string_ref (mem
, idx
) == '.')
2673 scm_t_bits shift
= 1;
2675 unsigned int digit_value
;
2676 SCM big_shift
= SCM_I_MAKINUM (1);
2681 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2682 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2687 digit_value
= DIGIT2UINT (c
);
2698 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2700 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2701 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2703 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2711 add
= add
* 10 + digit_value
;
2717 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2718 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2719 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2722 result
= scm_divide (result
, big_shift
);
2724 /* We've seen a decimal point, thus the value is implicitly inexact. */
2736 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2738 switch (scm_i_string_ref (mem
, idx
))
2750 c
= scm_i_string_ref (mem
, idx
);
2758 c
= scm_i_string_ref (mem
, idx
);
2767 c
= scm_i_string_ref (mem
, idx
);
2772 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2776 exponent
= DIGIT2UINT (c
);
2779 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2780 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2783 if (exponent
<= SCM_MAXEXP
)
2784 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2790 if (exponent
> SCM_MAXEXP
)
2792 size_t exp_len
= idx
- start
;
2793 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
2794 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2795 scm_out_of_range ("string->number", exp_num
);
2798 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2800 result
= scm_product (result
, e
);
2802 result
= scm_divide2real (result
, e
);
2804 /* We've seen an exponent, thus the value is implicitly inexact. */
2822 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2825 mem2ureal (SCM mem
, unsigned int *p_idx
,
2826 unsigned int radix
, enum t_exactness
*p_exactness
)
2828 unsigned int idx
= *p_idx
;
2830 size_t len
= scm_i_string_length (mem
);
2832 /* Start off believing that the number will be exact. This changes
2833 to INEXACT if we see a decimal point or a hash. */
2834 enum t_exactness x
= EXACT
;
2839 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
2845 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
2847 /* Cobble up the fractional part. We might want to set the
2848 NaN's mantissa from it. */
2850 mem2uinteger (mem
, &idx
, 10, &x
);
2855 if (scm_i_string_ref (mem
, idx
) == '.')
2859 else if (idx
+ 1 == len
)
2861 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
2864 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
,
2871 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
2872 if (scm_is_false (uinteger
))
2877 else if (scm_i_string_ref (mem
, idx
) == '/')
2885 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
2886 if (scm_is_false (divisor
))
2889 /* both are int/big here, I assume */
2890 result
= scm_i_make_ratio (uinteger
, divisor
);
2892 else if (radix
== 10)
2894 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
2895 if (scm_is_false (result
))
2904 /* Update *p_exactness if the number just read was inexact. This is
2905 important for complex numbers, so that a complex number is
2906 treated as inexact overall if either its real or imaginary part
2912 /* When returning an inexact zero, make sure it is represented as a
2913 floating point value so that we can change its sign.
2915 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2916 result
= scm_from_double (0.0);
2922 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2925 mem2complex (SCM mem
, unsigned int idx
,
2926 unsigned int radix
, enum t_exactness
*p_exactness
)
2931 size_t len
= scm_i_string_length (mem
);
2936 c
= scm_i_string_ref (mem
, idx
);
2951 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2952 if (scm_is_false (ureal
))
2954 /* input must be either +i or -i */
2959 if (scm_i_string_ref (mem
, idx
) == 'i'
2960 || scm_i_string_ref (mem
, idx
) == 'I')
2966 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2973 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2974 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2979 c
= scm_i_string_ref (mem
, idx
);
2983 /* either +<ureal>i or -<ureal>i */
2990 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2993 /* polar input: <real>@<real>. */
3004 c
= scm_i_string_ref (mem
, idx
);
3022 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3023 if (scm_is_false (angle
))
3028 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3029 angle
= scm_difference (angle
, SCM_UNDEFINED
);
3031 result
= scm_make_polar (ureal
, angle
);
3036 /* expecting input matching <real>[+-]<ureal>?i */
3043 int sign
= (c
== '+') ? 1 : -1;
3044 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
3046 if (scm_is_false (imag
))
3047 imag
= SCM_I_MAKINUM (sign
);
3048 else if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
3049 imag
= scm_difference (imag
, SCM_UNDEFINED
);
3053 if (scm_i_string_ref (mem
, idx
) != 'i'
3054 && scm_i_string_ref (mem
, idx
) != 'I')
3061 return scm_make_rectangular (ureal
, imag
);
3070 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
3072 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
3075 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
3077 unsigned int idx
= 0;
3078 unsigned int radix
= NO_RADIX
;
3079 enum t_exactness forced_x
= NO_EXACTNESS
;
3080 enum t_exactness implicit_x
= EXACT
;
3082 size_t len
= scm_i_string_length (mem
);
3084 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3085 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
3087 switch (scm_i_string_ref (mem
, idx
+ 1))
3090 if (radix
!= NO_RADIX
)
3095 if (radix
!= NO_RADIX
)
3100 if (forced_x
!= NO_EXACTNESS
)
3105 if (forced_x
!= NO_EXACTNESS
)
3110 if (radix
!= NO_RADIX
)
3115 if (radix
!= NO_RADIX
)
3125 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3126 if (radix
== NO_RADIX
)
3127 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
3129 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
3131 if (scm_is_false (result
))
3137 if (SCM_INEXACTP (result
))
3138 return scm_inexact_to_exact (result
);
3142 if (SCM_INEXACTP (result
))
3145 return scm_exact_to_inexact (result
);
3148 if (implicit_x
== INEXACT
)
3150 if (SCM_INEXACTP (result
))
3153 return scm_exact_to_inexact (result
);
3161 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
3162 unsigned int default_radix
)
3164 SCM str
= scm_from_locale_stringn (mem
, len
);
3166 return scm_i_string_to_number (str
, default_radix
);
3170 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3171 (SCM string
, SCM radix
),
3172 "Return a number of the maximally precise representation\n"
3173 "expressed by the given @var{string}. @var{radix} must be an\n"
3174 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3175 "is a default radix that may be overridden by an explicit radix\n"
3176 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3177 "supplied, then the default radix is 10. If string is not a\n"
3178 "syntactically valid notation for a number, then\n"
3179 "@code{string->number} returns @code{#f}.")
3180 #define FUNC_NAME s_scm_string_to_number
3184 SCM_VALIDATE_STRING (1, string
);
3186 if (SCM_UNBNDP (radix
))
3189 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3191 answer
= scm_i_string_to_number (string
, base
);
3192 scm_remember_upto_here_1 (string
);
3198 /*** END strs->nums ***/
3202 scm_bigequal (SCM x
, SCM y
)
3204 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3205 scm_remember_upto_here_2 (x
, y
);
3206 return scm_from_bool (0 == result
);
3210 scm_real_equalp (SCM x
, SCM y
)
3212 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3216 scm_complex_equalp (SCM x
, SCM y
)
3218 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3219 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3223 scm_i_fraction_equalp (SCM x
, SCM y
)
3225 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3226 SCM_FRACTION_NUMERATOR (y
)))
3227 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3228 SCM_FRACTION_DENOMINATOR (y
))))
3235 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3237 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3239 #define FUNC_NAME s_scm_number_p
3241 return scm_from_bool (SCM_NUMBERP (x
));
3245 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3247 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3248 "otherwise. Note that the sets of real, rational and integer\n"
3249 "values form subsets of the set of complex numbers, i. e. the\n"
3250 "predicate will also be fulfilled if @var{x} is a real,\n"
3251 "rational or integer number.")
3252 #define FUNC_NAME s_scm_complex_p
3254 /* all numbers are complex. */
3255 return scm_number_p (x
);
3259 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3261 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3262 "otherwise. Note that the set of integer values forms a subset of\n"
3263 "the set of real numbers, i. e. the predicate will also be\n"
3264 "fulfilled if @var{x} is an integer number.")
3265 #define FUNC_NAME s_scm_real_p
3267 /* we can't represent irrational numbers. */
3268 return scm_rational_p (x
);
3272 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3274 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3275 "otherwise. Note that the set of integer values forms a subset of\n"
3276 "the set of rational numbers, i. e. the predicate will also be\n"
3277 "fulfilled if @var{x} is an integer number.")
3278 #define FUNC_NAME s_scm_rational_p
3280 if (SCM_I_INUMP (x
))
3282 else if (SCM_IMP (x
))
3284 else if (SCM_BIGP (x
))
3286 else if (SCM_FRACTIONP (x
))
3288 else if (SCM_REALP (x
))
3289 /* due to their limited precision, all floating point numbers are
3290 rational as well. */
3297 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3299 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3301 #define FUNC_NAME s_scm_integer_p
3304 if (SCM_I_INUMP (x
))
3310 if (!SCM_INEXACTP (x
))
3312 if (SCM_COMPLEXP (x
))
3314 r
= SCM_REAL_VALUE (x
);
3315 /* +/-inf passes r==floor(r), making those #t */
3323 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3325 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3327 #define FUNC_NAME s_scm_inexact_p
3329 if (SCM_INEXACTP (x
))
3331 if (SCM_NUMBERP (x
))
3333 SCM_WRONG_TYPE_ARG (1, x
);
3338 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
3339 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
3340 (SCM x
, SCM y
, SCM rest
),
3341 "Return @code{#t} if all parameters are numerically equal.")
3342 #define FUNC_NAME s_scm_i_num_eq_p
3344 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3346 while (!scm_is_null (rest
))
3348 if (scm_is_false (scm_num_eq_p (x
, y
)))
3352 rest
= scm_cdr (rest
);
3354 return scm_num_eq_p (x
, y
);
3358 scm_num_eq_p (SCM x
, SCM y
)
3361 if (SCM_I_INUMP (x
))
3363 long xx
= SCM_I_INUM (x
);
3364 if (SCM_I_INUMP (y
))
3366 long yy
= SCM_I_INUM (y
);
3367 return scm_from_bool (xx
== yy
);
3369 else if (SCM_BIGP (y
))
3371 else if (SCM_REALP (y
))
3373 /* On a 32-bit system an inum fits a double, we can cast the inum
3374 to a double and compare.
3376 But on a 64-bit system an inum is bigger than a double and
3377 casting it to a double (call that dxx) will round. dxx is at
3378 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3379 an integer and fits a long. So we cast yy to a long and
3380 compare with plain xx.
3382 An alternative (for any size system actually) would be to check
3383 yy is an integer (with floor) and is in range of an inum
3384 (compare against appropriate powers of 2) then test
3385 xx==(long)yy. It's just a matter of which casts/comparisons
3386 might be fastest or easiest for the cpu. */
3388 double yy
= SCM_REAL_VALUE (y
);
3389 return scm_from_bool ((double) xx
== yy
3390 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3391 || xx
== (long) yy
));
3393 else if (SCM_COMPLEXP (y
))
3394 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3395 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3396 else if (SCM_FRACTIONP (y
))
3399 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3401 else if (SCM_BIGP (x
))
3403 if (SCM_I_INUMP (y
))
3405 else if (SCM_BIGP (y
))
3407 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3408 scm_remember_upto_here_2 (x
, y
);
3409 return scm_from_bool (0 == cmp
);
3411 else if (SCM_REALP (y
))
3414 if (xisnan (SCM_REAL_VALUE (y
)))
3416 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3417 scm_remember_upto_here_1 (x
);
3418 return scm_from_bool (0 == cmp
);
3420 else if (SCM_COMPLEXP (y
))
3423 if (0.0 != SCM_COMPLEX_IMAG (y
))
3425 if (xisnan (SCM_COMPLEX_REAL (y
)))
3427 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3428 scm_remember_upto_here_1 (x
);
3429 return scm_from_bool (0 == cmp
);
3431 else if (SCM_FRACTIONP (y
))
3434 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3436 else if (SCM_REALP (x
))
3438 double xx
= SCM_REAL_VALUE (x
);
3439 if (SCM_I_INUMP (y
))
3441 /* see comments with inum/real above */
3442 long yy
= SCM_I_INUM (y
);
3443 return scm_from_bool (xx
== (double) yy
3444 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3445 || (long) xx
== yy
));
3447 else if (SCM_BIGP (y
))
3450 if (xisnan (SCM_REAL_VALUE (x
)))
3452 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3453 scm_remember_upto_here_1 (y
);
3454 return scm_from_bool (0 == cmp
);
3456 else if (SCM_REALP (y
))
3457 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3458 else if (SCM_COMPLEXP (y
))
3459 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3460 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3461 else if (SCM_FRACTIONP (y
))
3463 double xx
= SCM_REAL_VALUE (x
);
3467 return scm_from_bool (xx
< 0.0);
3468 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3472 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3474 else if (SCM_COMPLEXP (x
))
3476 if (SCM_I_INUMP (y
))
3477 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3478 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3479 else if (SCM_BIGP (y
))
3482 if (0.0 != SCM_COMPLEX_IMAG (x
))
3484 if (xisnan (SCM_COMPLEX_REAL (x
)))
3486 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3487 scm_remember_upto_here_1 (y
);
3488 return scm_from_bool (0 == cmp
);
3490 else if (SCM_REALP (y
))
3491 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3492 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3493 else if (SCM_COMPLEXP (y
))
3494 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3495 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3496 else if (SCM_FRACTIONP (y
))
3499 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3501 xx
= SCM_COMPLEX_REAL (x
);
3505 return scm_from_bool (xx
< 0.0);
3506 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3510 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3512 else if (SCM_FRACTIONP (x
))
3514 if (SCM_I_INUMP (y
))
3516 else if (SCM_BIGP (y
))
3518 else if (SCM_REALP (y
))
3520 double yy
= SCM_REAL_VALUE (y
);
3524 return scm_from_bool (0.0 < yy
);
3525 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3528 else if (SCM_COMPLEXP (y
))
3531 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3533 yy
= SCM_COMPLEX_REAL (y
);
3537 return scm_from_bool (0.0 < yy
);
3538 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3541 else if (SCM_FRACTIONP (y
))
3542 return scm_i_fraction_equalp (x
, y
);
3544 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3547 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
3551 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3552 done are good for inums, but for bignums an answer can almost always be
3553 had by just examining a few high bits of the operands, as done by GMP in
3554 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3555 of the float exponent to take into account. */
3557 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
3558 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
3559 (SCM x
, SCM y
, SCM rest
),
3560 "Return @code{#t} if the list of parameters is monotonically\n"
3562 #define FUNC_NAME s_scm_i_num_less_p
3564 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3566 while (!scm_is_null (rest
))
3568 if (scm_is_false (scm_less_p (x
, y
)))
3572 rest
= scm_cdr (rest
);
3574 return scm_less_p (x
, y
);
3578 scm_less_p (SCM x
, SCM y
)
3581 if (SCM_I_INUMP (x
))
3583 long xx
= SCM_I_INUM (x
);
3584 if (SCM_I_INUMP (y
))
3586 long yy
= SCM_I_INUM (y
);
3587 return scm_from_bool (xx
< yy
);
3589 else if (SCM_BIGP (y
))
3591 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3592 scm_remember_upto_here_1 (y
);
3593 return scm_from_bool (sgn
> 0);
3595 else if (SCM_REALP (y
))
3596 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3597 else if (SCM_FRACTIONP (y
))
3599 /* "x < a/b" becomes "x*b < a" */
3601 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3602 y
= SCM_FRACTION_NUMERATOR (y
);
3606 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3608 else if (SCM_BIGP (x
))
3610 if (SCM_I_INUMP (y
))
3612 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3613 scm_remember_upto_here_1 (x
);
3614 return scm_from_bool (sgn
< 0);
3616 else if (SCM_BIGP (y
))
3618 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3619 scm_remember_upto_here_2 (x
, y
);
3620 return scm_from_bool (cmp
< 0);
3622 else if (SCM_REALP (y
))
3625 if (xisnan (SCM_REAL_VALUE (y
)))
3627 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3628 scm_remember_upto_here_1 (x
);
3629 return scm_from_bool (cmp
< 0);
3631 else if (SCM_FRACTIONP (y
))
3634 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3636 else if (SCM_REALP (x
))
3638 if (SCM_I_INUMP (y
))
3639 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3640 else if (SCM_BIGP (y
))
3643 if (xisnan (SCM_REAL_VALUE (x
)))
3645 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3646 scm_remember_upto_here_1 (y
);
3647 return scm_from_bool (cmp
> 0);
3649 else if (SCM_REALP (y
))
3650 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3651 else if (SCM_FRACTIONP (y
))
3653 double xx
= SCM_REAL_VALUE (x
);
3657 return scm_from_bool (xx
< 0.0);
3658 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3662 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3664 else if (SCM_FRACTIONP (x
))
3666 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3668 /* "a/b < y" becomes "a < y*b" */
3669 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3670 x
= SCM_FRACTION_NUMERATOR (x
);
3673 else if (SCM_REALP (y
))
3675 double yy
= SCM_REAL_VALUE (y
);
3679 return scm_from_bool (0.0 < yy
);
3680 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3683 else if (SCM_FRACTIONP (y
))
3685 /* "a/b < c/d" becomes "a*d < c*b" */
3686 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3687 SCM_FRACTION_DENOMINATOR (y
));
3688 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3689 SCM_FRACTION_DENOMINATOR (x
));
3695 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3698 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
3702 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
3703 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
3704 (SCM x
, SCM y
, SCM rest
),
3705 "Return @code{#t} if the list of parameters is monotonically\n"
3707 #define FUNC_NAME s_scm_i_num_gr_p
3709 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3711 while (!scm_is_null (rest
))
3713 if (scm_is_false (scm_gr_p (x
, y
)))
3717 rest
= scm_cdr (rest
);
3719 return scm_gr_p (x
, y
);
3722 #define FUNC_NAME s_scm_i_num_gr_p
3724 scm_gr_p (SCM x
, SCM y
)
3726 if (!SCM_NUMBERP (x
))
3727 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3728 else if (!SCM_NUMBERP (y
))
3729 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3731 return scm_less_p (y
, x
);
3736 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
3737 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
3738 (SCM x
, SCM y
, SCM rest
),
3739 "Return @code{#t} if the list of parameters is monotonically\n"
3741 #define FUNC_NAME s_scm_i_num_leq_p
3743 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3745 while (!scm_is_null (rest
))
3747 if (scm_is_false (scm_leq_p (x
, y
)))
3751 rest
= scm_cdr (rest
);
3753 return scm_leq_p (x
, y
);
3756 #define FUNC_NAME s_scm_i_num_leq_p
3758 scm_leq_p (SCM x
, SCM y
)
3760 if (!SCM_NUMBERP (x
))
3761 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3762 else if (!SCM_NUMBERP (y
))
3763 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3764 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3767 return scm_not (scm_less_p (y
, x
));
3772 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
3773 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
3774 (SCM x
, SCM y
, SCM rest
),
3775 "Return @code{#t} if the list of parameters is monotonically\n"
3777 #define FUNC_NAME s_scm_i_num_geq_p
3779 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3781 while (!scm_is_null (rest
))
3783 if (scm_is_false (scm_geq_p (x
, y
)))
3787 rest
= scm_cdr (rest
);
3789 return scm_geq_p (x
, y
);
3792 #define FUNC_NAME s_scm_i_num_geq_p
3794 scm_geq_p (SCM x
, SCM y
)
3796 if (!SCM_NUMBERP (x
))
3797 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3798 else if (!SCM_NUMBERP (y
))
3799 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3800 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3803 return scm_not (scm_less_p (x
, y
));
3808 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3809 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3815 if (SCM_I_INUMP (z
))
3816 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3817 else if (SCM_BIGP (z
))
3819 else if (SCM_REALP (z
))
3820 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3821 else if (SCM_COMPLEXP (z
))
3822 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3823 && SCM_COMPLEX_IMAG (z
) == 0.0);
3824 else if (SCM_FRACTIONP (z
))
3827 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3831 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3832 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3836 scm_positive_p (SCM x
)
3838 if (SCM_I_INUMP (x
))
3839 return scm_from_bool (SCM_I_INUM (x
) > 0);
3840 else if (SCM_BIGP (x
))
3842 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3843 scm_remember_upto_here_1 (x
);
3844 return scm_from_bool (sgn
> 0);
3846 else if (SCM_REALP (x
))
3847 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3848 else if (SCM_FRACTIONP (x
))
3849 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3851 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3855 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3856 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3860 scm_negative_p (SCM x
)
3862 if (SCM_I_INUMP (x
))
3863 return scm_from_bool (SCM_I_INUM (x
) < 0);
3864 else if (SCM_BIGP (x
))
3866 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3867 scm_remember_upto_here_1 (x
);
3868 return scm_from_bool (sgn
< 0);
3870 else if (SCM_REALP (x
))
3871 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3872 else if (SCM_FRACTIONP (x
))
3873 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3875 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3879 /* scm_min and scm_max return an inexact when either argument is inexact, as
3880 required by r5rs. On that basis, for exact/inexact combinations the
3881 exact is converted to inexact to compare and possibly return. This is
3882 unlike scm_less_p above which takes some trouble to preserve all bits in
3883 its test, such trouble is not required for min and max. */
3885 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
3886 (SCM x
, SCM y
, SCM rest
),
3887 "Return the maximum of all parameter values.")
3888 #define FUNC_NAME s_scm_i_max
3890 while (!scm_is_null (rest
))
3891 { x
= scm_max (x
, y
);
3893 rest
= scm_cdr (rest
);
3895 return scm_max (x
, y
);
3899 #define s_max s_scm_i_max
3900 #define g_max g_scm_i_max
3903 scm_max (SCM x
, SCM y
)
3908 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3909 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3912 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3915 if (SCM_I_INUMP (x
))
3917 long xx
= SCM_I_INUM (x
);
3918 if (SCM_I_INUMP (y
))
3920 long yy
= SCM_I_INUM (y
);
3921 return (xx
< yy
) ? y
: x
;
3923 else if (SCM_BIGP (y
))
3925 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3926 scm_remember_upto_here_1 (y
);
3927 return (sgn
< 0) ? x
: y
;
3929 else if (SCM_REALP (y
))
3932 /* if y==NaN then ">" is false and we return NaN */
3933 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3935 else if (SCM_FRACTIONP (y
))
3938 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3941 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3943 else if (SCM_BIGP (x
))
3945 if (SCM_I_INUMP (y
))
3947 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3948 scm_remember_upto_here_1 (x
);
3949 return (sgn
< 0) ? y
: x
;
3951 else if (SCM_BIGP (y
))
3953 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3954 scm_remember_upto_here_2 (x
, y
);
3955 return (cmp
> 0) ? x
: y
;
3957 else if (SCM_REALP (y
))
3959 /* if y==NaN then xx>yy is false, so we return the NaN y */
3962 xx
= scm_i_big2dbl (x
);
3963 yy
= SCM_REAL_VALUE (y
);
3964 return (xx
> yy
? scm_from_double (xx
) : y
);
3966 else if (SCM_FRACTIONP (y
))
3971 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3973 else if (SCM_REALP (x
))
3975 if (SCM_I_INUMP (y
))
3977 double z
= SCM_I_INUM (y
);
3978 /* if x==NaN then "<" is false and we return NaN */
3979 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3981 else if (SCM_BIGP (y
))
3986 else if (SCM_REALP (y
))
3988 /* if x==NaN then our explicit check means we return NaN
3989 if y==NaN then ">" is false and we return NaN
3990 calling isnan is unavoidable, since it's the only way to know
3991 which of x or y causes any compares to be false */
3992 double xx
= SCM_REAL_VALUE (x
);
3993 return (xisnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3995 else if (SCM_FRACTIONP (y
))
3997 double yy
= scm_i_fraction2double (y
);
3998 double xx
= SCM_REAL_VALUE (x
);
3999 return (xx
< yy
) ? scm_from_double (yy
) : x
;
4002 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4004 else if (SCM_FRACTIONP (x
))
4006 if (SCM_I_INUMP (y
))
4010 else if (SCM_BIGP (y
))
4014 else if (SCM_REALP (y
))
4016 double xx
= scm_i_fraction2double (x
);
4017 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
4019 else if (SCM_FRACTIONP (y
))
4024 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
4027 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
4031 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
4032 (SCM x
, SCM y
, SCM rest
),
4033 "Return the minimum of all parameter values.")
4034 #define FUNC_NAME s_scm_i_min
4036 while (!scm_is_null (rest
))
4037 { x
= scm_min (x
, y
);
4039 rest
= scm_cdr (rest
);
4041 return scm_min (x
, y
);
4045 #define s_min s_scm_i_min
4046 #define g_min g_scm_i_min
4049 scm_min (SCM x
, SCM y
)
4054 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
4055 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
4058 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
4061 if (SCM_I_INUMP (x
))
4063 long xx
= SCM_I_INUM (x
);
4064 if (SCM_I_INUMP (y
))
4066 long yy
= SCM_I_INUM (y
);
4067 return (xx
< yy
) ? x
: y
;
4069 else if (SCM_BIGP (y
))
4071 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4072 scm_remember_upto_here_1 (y
);
4073 return (sgn
< 0) ? y
: x
;
4075 else if (SCM_REALP (y
))
4078 /* if y==NaN then "<" is false and we return NaN */
4079 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
4081 else if (SCM_FRACTIONP (y
))
4084 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
4087 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4089 else if (SCM_BIGP (x
))
4091 if (SCM_I_INUMP (y
))
4093 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4094 scm_remember_upto_here_1 (x
);
4095 return (sgn
< 0) ? x
: y
;
4097 else if (SCM_BIGP (y
))
4099 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4100 scm_remember_upto_here_2 (x
, y
);
4101 return (cmp
> 0) ? y
: x
;
4103 else if (SCM_REALP (y
))
4105 /* if y==NaN then xx<yy is false, so we return the NaN y */
4108 xx
= scm_i_big2dbl (x
);
4109 yy
= SCM_REAL_VALUE (y
);
4110 return (xx
< yy
? scm_from_double (xx
) : y
);
4112 else if (SCM_FRACTIONP (y
))
4117 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4119 else if (SCM_REALP (x
))
4121 if (SCM_I_INUMP (y
))
4123 double z
= SCM_I_INUM (y
);
4124 /* if x==NaN then "<" is false and we return NaN */
4125 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
4127 else if (SCM_BIGP (y
))
4132 else if (SCM_REALP (y
))
4134 /* if x==NaN then our explicit check means we return NaN
4135 if y==NaN then "<" is false and we return NaN
4136 calling isnan is unavoidable, since it's the only way to know
4137 which of x or y causes any compares to be false */
4138 double xx
= SCM_REAL_VALUE (x
);
4139 return (xisnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
4141 else if (SCM_FRACTIONP (y
))
4143 double yy
= scm_i_fraction2double (y
);
4144 double xx
= SCM_REAL_VALUE (x
);
4145 return (yy
< xx
) ? scm_from_double (yy
) : x
;
4148 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4150 else if (SCM_FRACTIONP (x
))
4152 if (SCM_I_INUMP (y
))
4156 else if (SCM_BIGP (y
))
4160 else if (SCM_REALP (y
))
4162 double xx
= scm_i_fraction2double (x
);
4163 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
4165 else if (SCM_FRACTIONP (y
))
4170 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4173 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
4177 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
4178 (SCM x
, SCM y
, SCM rest
),
4179 "Return the sum of all parameter values. Return 0 if called without\n"
4181 #define FUNC_NAME s_scm_i_sum
4183 while (!scm_is_null (rest
))
4184 { x
= scm_sum (x
, y
);
4186 rest
= scm_cdr (rest
);
4188 return scm_sum (x
, y
);
4192 #define s_sum s_scm_i_sum
4193 #define g_sum g_scm_i_sum
4196 scm_sum (SCM x
, SCM y
)
4198 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4200 if (SCM_NUMBERP (x
)) return x
;
4201 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4202 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4205 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4207 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4209 long xx
= SCM_I_INUM (x
);
4210 long yy
= SCM_I_INUM (y
);
4211 long int z
= xx
+ yy
;
4212 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
4214 else if (SCM_BIGP (y
))
4219 else if (SCM_REALP (y
))
4221 long int xx
= SCM_I_INUM (x
);
4222 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4224 else if (SCM_COMPLEXP (y
))
4226 long int xx
= SCM_I_INUM (x
);
4227 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4228 SCM_COMPLEX_IMAG (y
));
4230 else if (SCM_FRACTIONP (y
))
4231 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4232 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4233 SCM_FRACTION_DENOMINATOR (y
));
4235 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4236 } else if (SCM_BIGP (x
))
4238 if (SCM_I_INUMP (y
))
4243 inum
= SCM_I_INUM (y
);
4246 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4249 SCM result
= scm_i_mkbig ();
4250 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4251 scm_remember_upto_here_1 (x
);
4252 /* we know the result will have to be a bignum */
4255 return scm_i_normbig (result
);
4259 SCM result
= scm_i_mkbig ();
4260 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4261 scm_remember_upto_here_1 (x
);
4262 /* we know the result will have to be a bignum */
4265 return scm_i_normbig (result
);
4268 else if (SCM_BIGP (y
))
4270 SCM result
= scm_i_mkbig ();
4271 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4272 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4273 mpz_add (SCM_I_BIG_MPZ (result
),
4276 scm_remember_upto_here_2 (x
, y
);
4277 /* we know the result will have to be a bignum */
4280 return scm_i_normbig (result
);
4282 else if (SCM_REALP (y
))
4284 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4285 scm_remember_upto_here_1 (x
);
4286 return scm_from_double (result
);
4288 else if (SCM_COMPLEXP (y
))
4290 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4291 + SCM_COMPLEX_REAL (y
));
4292 scm_remember_upto_here_1 (x
);
4293 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4295 else if (SCM_FRACTIONP (y
))
4296 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4297 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4298 SCM_FRACTION_DENOMINATOR (y
));
4300 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4302 else if (SCM_REALP (x
))
4304 if (SCM_I_INUMP (y
))
4305 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4306 else if (SCM_BIGP (y
))
4308 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4309 scm_remember_upto_here_1 (y
);
4310 return scm_from_double (result
);
4312 else if (SCM_REALP (y
))
4313 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4314 else if (SCM_COMPLEXP (y
))
4315 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4316 SCM_COMPLEX_IMAG (y
));
4317 else if (SCM_FRACTIONP (y
))
4318 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4320 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4322 else if (SCM_COMPLEXP (x
))
4324 if (SCM_I_INUMP (y
))
4325 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4326 SCM_COMPLEX_IMAG (x
));
4327 else if (SCM_BIGP (y
))
4329 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4330 + SCM_COMPLEX_REAL (x
));
4331 scm_remember_upto_here_1 (y
);
4332 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4334 else if (SCM_REALP (y
))
4335 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4336 SCM_COMPLEX_IMAG (x
));
4337 else if (SCM_COMPLEXP (y
))
4338 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4339 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4340 else if (SCM_FRACTIONP (y
))
4341 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4342 SCM_COMPLEX_IMAG (x
));
4344 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4346 else if (SCM_FRACTIONP (x
))
4348 if (SCM_I_INUMP (y
))
4349 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4350 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4351 SCM_FRACTION_DENOMINATOR (x
));
4352 else if (SCM_BIGP (y
))
4353 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4354 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4355 SCM_FRACTION_DENOMINATOR (x
));
4356 else if (SCM_REALP (y
))
4357 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4358 else if (SCM_COMPLEXP (y
))
4359 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4360 SCM_COMPLEX_IMAG (y
));
4361 else if (SCM_FRACTIONP (y
))
4362 /* a/b + c/d = (ad + bc) / bd */
4363 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4364 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4365 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4367 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4370 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4374 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4376 "Return @math{@var{x}+1}.")
4377 #define FUNC_NAME s_scm_oneplus
4379 return scm_sum (x
, SCM_I_MAKINUM (1));
4384 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
4385 (SCM x
, SCM y
, SCM rest
),
4386 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
4387 "the sum of all but the first argument are subtracted from the first\n"
4389 #define FUNC_NAME s_scm_i_difference
4391 while (!scm_is_null (rest
))
4392 { x
= scm_difference (x
, y
);
4394 rest
= scm_cdr (rest
);
4396 return scm_difference (x
, y
);
4400 #define s_difference s_scm_i_difference
4401 #define g_difference g_scm_i_difference
4404 scm_difference (SCM x
, SCM y
)
4405 #define FUNC_NAME s_difference
4407 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4410 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4412 if (SCM_I_INUMP (x
))
4414 long xx
= -SCM_I_INUM (x
);
4415 if (SCM_FIXABLE (xx
))
4416 return SCM_I_MAKINUM (xx
);
4418 return scm_i_long2big (xx
);
4420 else if (SCM_BIGP (x
))
4421 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4422 bignum, but negating that gives a fixnum. */
4423 return scm_i_normbig (scm_i_clonebig (x
, 0));
4424 else if (SCM_REALP (x
))
4425 return scm_from_double (-SCM_REAL_VALUE (x
));
4426 else if (SCM_COMPLEXP (x
))
4427 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4428 -SCM_COMPLEX_IMAG (x
));
4429 else if (SCM_FRACTIONP (x
))
4430 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4431 SCM_FRACTION_DENOMINATOR (x
));
4433 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4436 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4438 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4440 long int xx
= SCM_I_INUM (x
);
4441 long int yy
= SCM_I_INUM (y
);
4442 long int z
= xx
- yy
;
4443 if (SCM_FIXABLE (z
))
4444 return SCM_I_MAKINUM (z
);
4446 return scm_i_long2big (z
);
4448 else if (SCM_BIGP (y
))
4450 /* inum-x - big-y */
4451 long xx
= SCM_I_INUM (x
);
4454 return scm_i_clonebig (y
, 0);
4457 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4458 SCM result
= scm_i_mkbig ();
4461 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4464 /* x - y == -(y + -x) */
4465 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4466 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4468 scm_remember_upto_here_1 (y
);
4470 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4471 /* we know the result will have to be a bignum */
4474 return scm_i_normbig (result
);
4477 else if (SCM_REALP (y
))
4479 long int xx
= SCM_I_INUM (x
);
4480 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4482 else if (SCM_COMPLEXP (y
))
4484 long int xx
= SCM_I_INUM (x
);
4485 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4486 - SCM_COMPLEX_IMAG (y
));
4488 else if (SCM_FRACTIONP (y
))
4489 /* a - b/c = (ac - b) / c */
4490 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4491 SCM_FRACTION_NUMERATOR (y
)),
4492 SCM_FRACTION_DENOMINATOR (y
));
4494 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4496 else if (SCM_BIGP (x
))
4498 if (SCM_I_INUMP (y
))
4500 /* big-x - inum-y */
4501 long yy
= SCM_I_INUM (y
);
4502 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4504 scm_remember_upto_here_1 (x
);
4506 return (SCM_FIXABLE (-yy
) ?
4507 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4510 SCM result
= scm_i_mkbig ();
4513 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4515 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4516 scm_remember_upto_here_1 (x
);
4518 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4519 /* we know the result will have to be a bignum */
4522 return scm_i_normbig (result
);
4525 else if (SCM_BIGP (y
))
4527 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4528 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4529 SCM result
= scm_i_mkbig ();
4530 mpz_sub (SCM_I_BIG_MPZ (result
),
4533 scm_remember_upto_here_2 (x
, y
);
4534 /* we know the result will have to be a bignum */
4535 if ((sgn_x
== 1) && (sgn_y
== -1))
4537 if ((sgn_x
== -1) && (sgn_y
== 1))
4539 return scm_i_normbig (result
);
4541 else if (SCM_REALP (y
))
4543 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4544 scm_remember_upto_here_1 (x
);
4545 return scm_from_double (result
);
4547 else if (SCM_COMPLEXP (y
))
4549 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4550 - SCM_COMPLEX_REAL (y
));
4551 scm_remember_upto_here_1 (x
);
4552 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4554 else if (SCM_FRACTIONP (y
))
4555 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4556 SCM_FRACTION_NUMERATOR (y
)),
4557 SCM_FRACTION_DENOMINATOR (y
));
4558 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4560 else if (SCM_REALP (x
))
4562 if (SCM_I_INUMP (y
))
4563 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4564 else if (SCM_BIGP (y
))
4566 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4567 scm_remember_upto_here_1 (x
);
4568 return scm_from_double (result
);
4570 else if (SCM_REALP (y
))
4571 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4572 else if (SCM_COMPLEXP (y
))
4573 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4574 -SCM_COMPLEX_IMAG (y
));
4575 else if (SCM_FRACTIONP (y
))
4576 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4578 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4580 else if (SCM_COMPLEXP (x
))
4582 if (SCM_I_INUMP (y
))
4583 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4584 SCM_COMPLEX_IMAG (x
));
4585 else if (SCM_BIGP (y
))
4587 double real_part
= (SCM_COMPLEX_REAL (x
)
4588 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4589 scm_remember_upto_here_1 (x
);
4590 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4592 else if (SCM_REALP (y
))
4593 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4594 SCM_COMPLEX_IMAG (x
));
4595 else if (SCM_COMPLEXP (y
))
4596 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4597 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4598 else if (SCM_FRACTIONP (y
))
4599 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4600 SCM_COMPLEX_IMAG (x
));
4602 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4604 else if (SCM_FRACTIONP (x
))
4606 if (SCM_I_INUMP (y
))
4607 /* a/b - c = (a - cb) / b */
4608 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4609 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4610 SCM_FRACTION_DENOMINATOR (x
));
4611 else if (SCM_BIGP (y
))
4612 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4613 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4614 SCM_FRACTION_DENOMINATOR (x
));
4615 else if (SCM_REALP (y
))
4616 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4617 else if (SCM_COMPLEXP (y
))
4618 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4619 -SCM_COMPLEX_IMAG (y
));
4620 else if (SCM_FRACTIONP (y
))
4621 /* a/b - c/d = (ad - bc) / bd */
4622 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4623 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4624 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4626 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4629 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4634 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4636 "Return @math{@var{x}-1}.")
4637 #define FUNC_NAME s_scm_oneminus
4639 return scm_difference (x
, SCM_I_MAKINUM (1));
4644 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
4645 (SCM x
, SCM y
, SCM rest
),
4646 "Return the product of all arguments. If called without arguments,\n"
4648 #define FUNC_NAME s_scm_i_product
4650 while (!scm_is_null (rest
))
4651 { x
= scm_product (x
, y
);
4653 rest
= scm_cdr (rest
);
4655 return scm_product (x
, y
);
4659 #define s_product s_scm_i_product
4660 #define g_product g_scm_i_product
4663 scm_product (SCM x
, SCM y
)
4665 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4668 return SCM_I_MAKINUM (1L);
4669 else if (SCM_NUMBERP (x
))
4672 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4675 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4680 xx
= SCM_I_INUM (x
);
4684 case 0: return x
; break;
4685 case 1: return y
; break;
4688 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4690 long yy
= SCM_I_INUM (y
);
4692 SCM k
= SCM_I_MAKINUM (kk
);
4693 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4697 SCM result
= scm_i_long2big (xx
);
4698 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4699 return scm_i_normbig (result
);
4702 else if (SCM_BIGP (y
))
4704 SCM result
= scm_i_mkbig ();
4705 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4706 scm_remember_upto_here_1 (y
);
4709 else if (SCM_REALP (y
))
4710 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4711 else if (SCM_COMPLEXP (y
))
4712 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4713 xx
* SCM_COMPLEX_IMAG (y
));
4714 else if (SCM_FRACTIONP (y
))
4715 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4716 SCM_FRACTION_DENOMINATOR (y
));
4718 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4720 else if (SCM_BIGP (x
))
4722 if (SCM_I_INUMP (y
))
4727 else if (SCM_BIGP (y
))
4729 SCM result
= scm_i_mkbig ();
4730 mpz_mul (SCM_I_BIG_MPZ (result
),
4733 scm_remember_upto_here_2 (x
, y
);
4736 else if (SCM_REALP (y
))
4738 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4739 scm_remember_upto_here_1 (x
);
4740 return scm_from_double (result
);
4742 else if (SCM_COMPLEXP (y
))
4744 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4745 scm_remember_upto_here_1 (x
);
4746 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4747 z
* SCM_COMPLEX_IMAG (y
));
4749 else if (SCM_FRACTIONP (y
))
4750 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4751 SCM_FRACTION_DENOMINATOR (y
));
4753 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4755 else if (SCM_REALP (x
))
4757 if (SCM_I_INUMP (y
))
4759 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4760 if (scm_is_eq (y
, SCM_INUM0
))
4762 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4764 else if (SCM_BIGP (y
))
4766 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4767 scm_remember_upto_here_1 (y
);
4768 return scm_from_double (result
);
4770 else if (SCM_REALP (y
))
4771 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4772 else if (SCM_COMPLEXP (y
))
4773 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4774 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4775 else if (SCM_FRACTIONP (y
))
4776 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4778 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4780 else if (SCM_COMPLEXP (x
))
4782 if (SCM_I_INUMP (y
))
4784 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4785 if (scm_is_eq (y
, SCM_INUM0
))
4787 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4788 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4790 else if (SCM_BIGP (y
))
4792 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4793 scm_remember_upto_here_1 (y
);
4794 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4795 z
* SCM_COMPLEX_IMAG (x
));
4797 else if (SCM_REALP (y
))
4798 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4799 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4800 else if (SCM_COMPLEXP (y
))
4802 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4803 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4804 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4805 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4807 else if (SCM_FRACTIONP (y
))
4809 double yy
= scm_i_fraction2double (y
);
4810 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4811 yy
* SCM_COMPLEX_IMAG (x
));
4814 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4816 else if (SCM_FRACTIONP (x
))
4818 if (SCM_I_INUMP (y
))
4819 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4820 SCM_FRACTION_DENOMINATOR (x
));
4821 else if (SCM_BIGP (y
))
4822 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4823 SCM_FRACTION_DENOMINATOR (x
));
4824 else if (SCM_REALP (y
))
4825 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4826 else if (SCM_COMPLEXP (y
))
4828 double xx
= scm_i_fraction2double (x
);
4829 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4830 xx
* SCM_COMPLEX_IMAG (y
));
4832 else if (SCM_FRACTIONP (y
))
4833 /* a/b * c/d = ac / bd */
4834 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4835 SCM_FRACTION_NUMERATOR (y
)),
4836 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4837 SCM_FRACTION_DENOMINATOR (y
)));
4839 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4842 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4845 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4846 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4847 #define ALLOW_DIVIDE_BY_ZERO
4848 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4851 /* The code below for complex division is adapted from the GNU
4852 libstdc++, which adapted it from f2c's libF77, and is subject to
4855 /****************************************************************
4856 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4858 Permission to use, copy, modify, and distribute this software
4859 and its documentation for any purpose and without fee is hereby
4860 granted, provided that the above copyright notice appear in all
4861 copies and that both that the copyright notice and this
4862 permission notice and warranty disclaimer appear in supporting
4863 documentation, and that the names of AT&T Bell Laboratories or
4864 Bellcore or any of their entities not be used in advertising or
4865 publicity pertaining to distribution of the software without
4866 specific, written prior permission.
4868 AT&T and Bellcore disclaim all warranties with regard to this
4869 software, including all implied warranties of merchantability
4870 and fitness. In no event shall AT&T or Bellcore be liable for
4871 any special, indirect or consequential damages or any damages
4872 whatsoever resulting from loss of use, data or profits, whether
4873 in an action of contract, negligence or other tortious action,
4874 arising out of or in connection with the use or performance of
4876 ****************************************************************/
4878 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
4879 (SCM x
, SCM y
, SCM rest
),
4880 "Divide the first argument by the product of the remaining\n"
4881 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
4883 #define FUNC_NAME s_scm_i_divide
4885 while (!scm_is_null (rest
))
4886 { x
= scm_divide (x
, y
);
4888 rest
= scm_cdr (rest
);
4890 return scm_divide (x
, y
);
4894 #define s_divide s_scm_i_divide
4895 #define g_divide g_scm_i_divide
4898 do_divide (SCM x
, SCM y
, int inexact
)
4899 #define FUNC_NAME s_divide
4903 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4906 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4907 else if (SCM_I_INUMP (x
))
4909 long xx
= SCM_I_INUM (x
);
4910 if (xx
== 1 || xx
== -1)
4912 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4914 scm_num_overflow (s_divide
);
4919 return scm_from_double (1.0 / (double) xx
);
4920 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4923 else if (SCM_BIGP (x
))
4926 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4927 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4929 else if (SCM_REALP (x
))
4931 double xx
= SCM_REAL_VALUE (x
);
4932 #ifndef ALLOW_DIVIDE_BY_ZERO
4934 scm_num_overflow (s_divide
);
4937 return scm_from_double (1.0 / xx
);
4939 else if (SCM_COMPLEXP (x
))
4941 double r
= SCM_COMPLEX_REAL (x
);
4942 double i
= SCM_COMPLEX_IMAG (x
);
4943 if (fabs(r
) <= fabs(i
))
4946 double d
= i
* (1.0 + t
* t
);
4947 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4952 double d
= r
* (1.0 + t
* t
);
4953 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4956 else if (SCM_FRACTIONP (x
))
4957 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4958 SCM_FRACTION_NUMERATOR (x
));
4960 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4963 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4965 long xx
= SCM_I_INUM (x
);
4966 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4968 long yy
= SCM_I_INUM (y
);
4971 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4972 scm_num_overflow (s_divide
);
4974 return scm_from_double ((double) xx
/ (double) yy
);
4977 else if (xx
% yy
!= 0)
4980 return scm_from_double ((double) xx
/ (double) yy
);
4981 else return scm_i_make_ratio (x
, y
);
4986 if (SCM_FIXABLE (z
))
4987 return SCM_I_MAKINUM (z
);
4989 return scm_i_long2big (z
);
4992 else if (SCM_BIGP (y
))
4995 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4996 else return scm_i_make_ratio (x
, y
);
4998 else if (SCM_REALP (y
))
5000 double yy
= SCM_REAL_VALUE (y
);
5001 #ifndef ALLOW_DIVIDE_BY_ZERO
5003 scm_num_overflow (s_divide
);
5006 return scm_from_double ((double) xx
/ yy
);
5008 else if (SCM_COMPLEXP (y
))
5011 complex_div
: /* y _must_ be a complex number */
5013 double r
= SCM_COMPLEX_REAL (y
);
5014 double i
= SCM_COMPLEX_IMAG (y
);
5015 if (fabs(r
) <= fabs(i
))
5018 double d
= i
* (1.0 + t
* t
);
5019 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
5024 double d
= r
* (1.0 + t
* t
);
5025 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
5029 else if (SCM_FRACTIONP (y
))
5030 /* a / b/c = ac / b */
5031 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5032 SCM_FRACTION_NUMERATOR (y
));
5034 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5036 else if (SCM_BIGP (x
))
5038 if (SCM_I_INUMP (y
))
5040 long int yy
= SCM_I_INUM (y
);
5043 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5044 scm_num_overflow (s_divide
);
5046 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5047 scm_remember_upto_here_1 (x
);
5048 return (sgn
== 0) ? scm_nan () : scm_inf ();
5055 /* FIXME: HMM, what are the relative performance issues here?
5056 We need to test. Is it faster on average to test
5057 divisible_p, then perform whichever operation, or is it
5058 faster to perform the integer div opportunistically and
5059 switch to real if there's a remainder? For now we take the
5060 middle ground: test, then if divisible, use the faster div
5063 long abs_yy
= yy
< 0 ? -yy
: yy
;
5064 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
5068 SCM result
= scm_i_mkbig ();
5069 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
5070 scm_remember_upto_here_1 (x
);
5072 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5073 return scm_i_normbig (result
);
5078 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
5079 else return scm_i_make_ratio (x
, y
);
5083 else if (SCM_BIGP (y
))
5085 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
5088 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5089 scm_num_overflow (s_divide
);
5091 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5092 scm_remember_upto_here_1 (x
);
5093 return (sgn
== 0) ? scm_nan () : scm_inf ();
5101 /* It's easily possible for the ratio x/y to fit a double
5102 but one or both x and y be too big to fit a double,
5103 hence the use of mpq_get_d rather than converting and
5106 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
5107 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
5108 return scm_from_double (mpq_get_d (q
));
5112 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
5116 SCM result
= scm_i_mkbig ();
5117 mpz_divexact (SCM_I_BIG_MPZ (result
),
5120 scm_remember_upto_here_2 (x
, y
);
5121 return scm_i_normbig (result
);
5124 return scm_i_make_ratio (x
, y
);
5128 else if (SCM_REALP (y
))
5130 double yy
= SCM_REAL_VALUE (y
);
5131 #ifndef ALLOW_DIVIDE_BY_ZERO
5133 scm_num_overflow (s_divide
);
5136 return scm_from_double (scm_i_big2dbl (x
) / yy
);
5138 else if (SCM_COMPLEXP (y
))
5140 a
= scm_i_big2dbl (x
);
5143 else if (SCM_FRACTIONP (y
))
5144 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5145 SCM_FRACTION_NUMERATOR (y
));
5147 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5149 else if (SCM_REALP (x
))
5151 double rx
= SCM_REAL_VALUE (x
);
5152 if (SCM_I_INUMP (y
))
5154 long int yy
= SCM_I_INUM (y
);
5155 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5157 scm_num_overflow (s_divide
);
5160 return scm_from_double (rx
/ (double) yy
);
5162 else if (SCM_BIGP (y
))
5164 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5165 scm_remember_upto_here_1 (y
);
5166 return scm_from_double (rx
/ dby
);
5168 else if (SCM_REALP (y
))
5170 double yy
= SCM_REAL_VALUE (y
);
5171 #ifndef ALLOW_DIVIDE_BY_ZERO
5173 scm_num_overflow (s_divide
);
5176 return scm_from_double (rx
/ yy
);
5178 else if (SCM_COMPLEXP (y
))
5183 else if (SCM_FRACTIONP (y
))
5184 return scm_from_double (rx
/ scm_i_fraction2double (y
));
5186 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5188 else if (SCM_COMPLEXP (x
))
5190 double rx
= SCM_COMPLEX_REAL (x
);
5191 double ix
= SCM_COMPLEX_IMAG (x
);
5192 if (SCM_I_INUMP (y
))
5194 long int yy
= SCM_I_INUM (y
);
5195 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5197 scm_num_overflow (s_divide
);
5202 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
5205 else if (SCM_BIGP (y
))
5207 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5208 scm_remember_upto_here_1 (y
);
5209 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
5211 else if (SCM_REALP (y
))
5213 double yy
= SCM_REAL_VALUE (y
);
5214 #ifndef ALLOW_DIVIDE_BY_ZERO
5216 scm_num_overflow (s_divide
);
5219 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5221 else if (SCM_COMPLEXP (y
))
5223 double ry
= SCM_COMPLEX_REAL (y
);
5224 double iy
= SCM_COMPLEX_IMAG (y
);
5225 if (fabs(ry
) <= fabs(iy
))
5228 double d
= iy
* (1.0 + t
* t
);
5229 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5234 double d
= ry
* (1.0 + t
* t
);
5235 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5238 else if (SCM_FRACTIONP (y
))
5240 double yy
= scm_i_fraction2double (y
);
5241 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5244 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5246 else if (SCM_FRACTIONP (x
))
5248 if (SCM_I_INUMP (y
))
5250 long int yy
= SCM_I_INUM (y
);
5251 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5253 scm_num_overflow (s_divide
);
5256 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5257 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5259 else if (SCM_BIGP (y
))
5261 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5262 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5264 else if (SCM_REALP (y
))
5266 double yy
= SCM_REAL_VALUE (y
);
5267 #ifndef ALLOW_DIVIDE_BY_ZERO
5269 scm_num_overflow (s_divide
);
5272 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5274 else if (SCM_COMPLEXP (y
))
5276 a
= scm_i_fraction2double (x
);
5279 else if (SCM_FRACTIONP (y
))
5280 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5281 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5283 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5286 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5290 scm_divide (SCM x
, SCM y
)
5292 return do_divide (x
, y
, 0);
5295 static SCM
scm_divide2real (SCM x
, SCM y
)
5297 return do_divide (x
, y
, 1);
5303 scm_c_truncate (double x
)
5314 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5315 half-way case (ie. when x is an integer plus 0.5) going upwards.
5316 Then half-way cases are identified and adjusted down if the
5317 round-upwards didn't give the desired even integer.
5319 "plus_half == result" identifies a half-way case. If plus_half, which is
5320 x + 0.5, is an integer then x must be an integer plus 0.5.
5322 An odd "result" value is identified with result/2 != floor(result/2).
5323 This is done with plus_half, since that value is ready for use sooner in
5324 a pipelined cpu, and we're already requiring plus_half == result.
5326 Note however that we need to be careful when x is big and already an
5327 integer. In that case "x+0.5" may round to an adjacent integer, causing
5328 us to return such a value, incorrectly. For instance if the hardware is
5329 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5330 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5331 returned. Or if the hardware is in round-upwards mode, then other bigger
5332 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5333 representable value, 2^128+2^76 (or whatever), again incorrect.
5335 These bad roundings of x+0.5 are avoided by testing at the start whether
5336 x is already an integer. If it is then clearly that's the desired result
5337 already. And if it's not then the exponent must be small enough to allow
5338 an 0.5 to be represented, and hence added without a bad rounding. */
5341 scm_c_round (double x
)
5343 double plus_half
, result
;
5348 plus_half
= x
+ 0.5;
5349 result
= floor (plus_half
);
5350 /* Adjust so that the rounding is towards even. */
5351 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5356 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5358 "Round the number @var{x} towards zero.")
5359 #define FUNC_NAME s_scm_truncate_number
5361 if (scm_is_false (scm_negative_p (x
)))
5362 return scm_floor (x
);
5364 return scm_ceiling (x
);
5368 static SCM exactly_one_half
;
5370 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5372 "Round the number @var{x} towards the nearest integer. "
5373 "When it is exactly halfway between two integers, "
5374 "round towards the even one.")
5375 #define FUNC_NAME s_scm_round_number
5377 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5379 else if (SCM_REALP (x
))
5380 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5383 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5384 single quotient+remainder division then examining to see which way
5385 the rounding should go. */
5386 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5387 SCM result
= scm_floor (plus_half
);
5388 /* Adjust so that the rounding is towards even. */
5389 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5390 && scm_is_true (scm_odd_p (result
)))
5391 return scm_difference (result
, SCM_I_MAKINUM (1));
5398 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5400 "Round the number @var{x} towards minus infinity.")
5401 #define FUNC_NAME s_scm_floor
5403 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5405 else if (SCM_REALP (x
))
5406 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5407 else if (SCM_FRACTIONP (x
))
5409 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5410 SCM_FRACTION_DENOMINATOR (x
));
5411 if (scm_is_false (scm_negative_p (x
)))
5413 /* For positive x, rounding towards zero is correct. */
5418 /* For negative x, we need to return q-1 unless x is an
5419 integer. But fractions are never integer, per our
5421 return scm_difference (q
, SCM_I_MAKINUM (1));
5425 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5429 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5431 "Round the number @var{x} towards infinity.")
5432 #define FUNC_NAME s_scm_ceiling
5434 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5436 else if (SCM_REALP (x
))
5437 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5438 else if (SCM_FRACTIONP (x
))
5440 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5441 SCM_FRACTION_DENOMINATOR (x
));
5442 if (scm_is_false (scm_positive_p (x
)))
5444 /* For negative x, rounding towards zero is correct. */
5449 /* For positive x, we need to return q+1 unless x is an
5450 integer. But fractions are never integer, per our
5452 return scm_sum (q
, SCM_I_MAKINUM (1));
5456 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5460 /* sin/cos/tan/asin/acos/atan
5461 sinh/cosh/tanh/asinh/acosh/atanh
5462 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5463 Written by Jerry D. Hedden, (C) FSF.
5464 See the file `COPYING' for terms applying to this program. */
5466 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5468 "Return @var{x} raised to the power of @var{y}.")
5469 #define FUNC_NAME s_scm_expt
5471 if (scm_is_true (scm_exact_p (x
)) && scm_is_integer (y
))
5472 return scm_integer_expt (x
, y
);
5473 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5475 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5478 return scm_exp (scm_product (scm_log (x
), y
));
5482 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5484 "Compute the sine of @var{z}.")
5485 #define FUNC_NAME s_scm_sin
5487 if (scm_is_real (z
))
5488 return scm_from_double (sin (scm_to_double (z
)));
5489 else if (SCM_COMPLEXP (z
))
5491 x
= SCM_COMPLEX_REAL (z
);
5492 y
= SCM_COMPLEX_IMAG (z
);
5493 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5494 cos (x
) * sinh (y
));
5497 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5501 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5503 "Compute the cosine of @var{z}.")
5504 #define FUNC_NAME s_scm_cos
5506 if (scm_is_real (z
))
5507 return scm_from_double (cos (scm_to_double (z
)));
5508 else if (SCM_COMPLEXP (z
))
5510 x
= SCM_COMPLEX_REAL (z
);
5511 y
= SCM_COMPLEX_IMAG (z
);
5512 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5513 -sin (x
) * sinh (y
));
5516 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5520 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5522 "Compute the tangent of @var{z}.")
5523 #define FUNC_NAME s_scm_tan
5525 if (scm_is_real (z
))
5526 return scm_from_double (tan (scm_to_double (z
)));
5527 else if (SCM_COMPLEXP (z
))
5529 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5530 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5531 w
= cos (x
) + cosh (y
);
5532 #ifndef ALLOW_DIVIDE_BY_ZERO
5534 scm_num_overflow (s_scm_tan
);
5536 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5539 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5543 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5545 "Compute the hyperbolic sine of @var{z}.")
5546 #define FUNC_NAME s_scm_sinh
5548 if (scm_is_real (z
))
5549 return scm_from_double (sinh (scm_to_double (z
)));
5550 else if (SCM_COMPLEXP (z
))
5552 x
= SCM_COMPLEX_REAL (z
);
5553 y
= SCM_COMPLEX_IMAG (z
);
5554 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5555 cosh (x
) * sin (y
));
5558 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5562 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5564 "Compute the hyperbolic cosine of @var{z}.")
5565 #define FUNC_NAME s_scm_cosh
5567 if (scm_is_real (z
))
5568 return scm_from_double (cosh (scm_to_double (z
)));
5569 else if (SCM_COMPLEXP (z
))
5571 x
= SCM_COMPLEX_REAL (z
);
5572 y
= SCM_COMPLEX_IMAG (z
);
5573 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5574 sinh (x
) * sin (y
));
5577 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5581 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5583 "Compute the hyperbolic tangent of @var{z}.")
5584 #define FUNC_NAME s_scm_tanh
5586 if (scm_is_real (z
))
5587 return scm_from_double (tanh (scm_to_double (z
)));
5588 else if (SCM_COMPLEXP (z
))
5590 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5591 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5592 w
= cosh (x
) + cos (y
);
5593 #ifndef ALLOW_DIVIDE_BY_ZERO
5595 scm_num_overflow (s_scm_tanh
);
5597 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5600 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5604 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5606 "Compute the arc sine of @var{z}.")
5607 #define FUNC_NAME s_scm_asin
5609 if (scm_is_real (z
))
5611 double w
= scm_to_double (z
);
5612 if (w
>= -1.0 && w
<= 1.0)
5613 return scm_from_double (asin (w
));
5615 return scm_product (scm_c_make_rectangular (0, -1),
5616 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5618 else if (SCM_COMPLEXP (z
))
5620 x
= SCM_COMPLEX_REAL (z
);
5621 y
= SCM_COMPLEX_IMAG (z
);
5622 return scm_product (scm_c_make_rectangular (0, -1),
5623 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5626 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5630 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5632 "Compute the arc cosine of @var{z}.")
5633 #define FUNC_NAME s_scm_acos
5635 if (scm_is_real (z
))
5637 double w
= scm_to_double (z
);
5638 if (w
>= -1.0 && w
<= 1.0)
5639 return scm_from_double (acos (w
));
5641 return scm_sum (scm_from_double (acos (0.0)),
5642 scm_product (scm_c_make_rectangular (0, 1),
5643 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5645 else if (SCM_COMPLEXP (z
))
5647 x
= SCM_COMPLEX_REAL (z
);
5648 y
= SCM_COMPLEX_IMAG (z
);
5649 return scm_sum (scm_from_double (acos (0.0)),
5650 scm_product (scm_c_make_rectangular (0, 1),
5651 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5654 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5658 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5660 "With one argument, compute the arc tangent of @var{z}.\n"
5661 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5662 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5663 #define FUNC_NAME s_scm_atan
5667 if (scm_is_real (z
))
5668 return scm_from_double (atan (scm_to_double (z
)));
5669 else if (SCM_COMPLEXP (z
))
5672 v
= SCM_COMPLEX_REAL (z
);
5673 w
= SCM_COMPLEX_IMAG (z
);
5674 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5675 scm_c_make_rectangular (v
, w
+ 1.0))),
5676 scm_c_make_rectangular (0, 2));
5679 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5681 else if (scm_is_real (z
))
5683 if (scm_is_real (y
))
5684 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5686 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5689 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5693 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5695 "Compute the inverse hyperbolic sine of @var{z}.")
5696 #define FUNC_NAME s_scm_sys_asinh
5698 if (scm_is_real (z
))
5699 return scm_from_double (asinh (scm_to_double (z
)));
5700 else if (scm_is_number (z
))
5701 return scm_log (scm_sum (z
,
5702 scm_sqrt (scm_sum (scm_product (z
, z
),
5703 SCM_I_MAKINUM (1)))));
5705 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5709 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5711 "Compute the inverse hyperbolic cosine of @var{z}.")
5712 #define FUNC_NAME s_scm_sys_acosh
5714 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5715 return scm_from_double (acosh (scm_to_double (z
)));
5716 else if (scm_is_number (z
))
5717 return scm_log (scm_sum (z
,
5718 scm_sqrt (scm_difference (scm_product (z
, z
),
5719 SCM_I_MAKINUM (1)))));
5721 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5725 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5727 "Compute the inverse hyperbolic tangent of @var{z}.")
5728 #define FUNC_NAME s_scm_sys_atanh
5730 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5731 return scm_from_double (atanh (scm_to_double (z
)));
5732 else if (scm_is_number (z
))
5733 return scm_divide (scm_log (scm_divide (scm_sum (SCM_I_MAKINUM (1), z
),
5734 scm_difference (SCM_I_MAKINUM (1), z
))),
5737 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5742 scm_c_make_rectangular (double re
, double im
)
5745 return scm_from_double (re
);
5749 SCM_NEWSMOB (z
, scm_tc16_complex
,
5750 scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5752 SCM_COMPLEX_REAL (z
) = re
;
5753 SCM_COMPLEX_IMAG (z
) = im
;
5758 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5759 (SCM real_part
, SCM imaginary_part
),
5760 "Return a complex number constructed of the given @var{real-part} "
5761 "and @var{imaginary-part} parts.")
5762 #define FUNC_NAME s_scm_make_rectangular
5764 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5765 SCM_ARG1
, FUNC_NAME
, "real");
5766 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5767 SCM_ARG2
, FUNC_NAME
, "real");
5768 return scm_c_make_rectangular (scm_to_double (real_part
),
5769 scm_to_double (imaginary_part
));
5774 scm_c_make_polar (double mag
, double ang
)
5778 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5779 use it on Glibc-based systems that have it (it's a GNU extension). See
5780 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5782 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5783 sincos (ang
, &s
, &c
);
5788 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5791 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5793 "Return the complex number @var{x} * e^(i * @var{y}).")
5794 #define FUNC_NAME s_scm_make_polar
5796 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5797 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5798 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5803 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5804 /* "Return the real part of the number @var{z}."
5807 scm_real_part (SCM z
)
5809 if (SCM_I_INUMP (z
))
5811 else if (SCM_BIGP (z
))
5813 else if (SCM_REALP (z
))
5815 else if (SCM_COMPLEXP (z
))
5816 return scm_from_double (SCM_COMPLEX_REAL (z
));
5817 else if (SCM_FRACTIONP (z
))
5820 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5824 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5825 /* "Return the imaginary part of the number @var{z}."
5828 scm_imag_part (SCM z
)
5830 if (SCM_I_INUMP (z
))
5832 else if (SCM_BIGP (z
))
5834 else if (SCM_REALP (z
))
5836 else if (SCM_COMPLEXP (z
))
5837 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5838 else if (SCM_FRACTIONP (z
))
5841 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5844 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5845 /* "Return the numerator of the number @var{z}."
5848 scm_numerator (SCM z
)
5850 if (SCM_I_INUMP (z
))
5852 else if (SCM_BIGP (z
))
5854 else if (SCM_FRACTIONP (z
))
5855 return SCM_FRACTION_NUMERATOR (z
);
5856 else if (SCM_REALP (z
))
5857 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5859 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5863 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5864 /* "Return the denominator of the number @var{z}."
5867 scm_denominator (SCM z
)
5869 if (SCM_I_INUMP (z
))
5870 return SCM_I_MAKINUM (1);
5871 else if (SCM_BIGP (z
))
5872 return SCM_I_MAKINUM (1);
5873 else if (SCM_FRACTIONP (z
))
5874 return SCM_FRACTION_DENOMINATOR (z
);
5875 else if (SCM_REALP (z
))
5876 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5878 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5881 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5882 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5883 * "@code{abs} for real arguments, but also allows complex numbers."
5886 scm_magnitude (SCM z
)
5888 if (SCM_I_INUMP (z
))
5890 long int zz
= SCM_I_INUM (z
);
5893 else if (SCM_POSFIXABLE (-zz
))
5894 return SCM_I_MAKINUM (-zz
);
5896 return scm_i_long2big (-zz
);
5898 else if (SCM_BIGP (z
))
5900 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5901 scm_remember_upto_here_1 (z
);
5903 return scm_i_clonebig (z
, 0);
5907 else if (SCM_REALP (z
))
5908 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5909 else if (SCM_COMPLEXP (z
))
5910 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5911 else if (SCM_FRACTIONP (z
))
5913 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5915 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5916 SCM_FRACTION_DENOMINATOR (z
));
5919 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5923 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5924 /* "Return the angle of the complex number @var{z}."
5929 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5930 flo0 to save allocating a new flonum with scm_from_double each time.
5931 But if atan2 follows the floating point rounding mode, then the value
5932 is not a constant. Maybe it'd be close enough though. */
5933 if (SCM_I_INUMP (z
))
5935 if (SCM_I_INUM (z
) >= 0)
5938 return scm_from_double (atan2 (0.0, -1.0));
5940 else if (SCM_BIGP (z
))
5942 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5943 scm_remember_upto_here_1 (z
);
5945 return scm_from_double (atan2 (0.0, -1.0));
5949 else if (SCM_REALP (z
))
5951 if (SCM_REAL_VALUE (z
) >= 0)
5954 return scm_from_double (atan2 (0.0, -1.0));
5956 else if (SCM_COMPLEXP (z
))
5957 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5958 else if (SCM_FRACTIONP (z
))
5960 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5962 else return scm_from_double (atan2 (0.0, -1.0));
5965 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5969 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5970 /* Convert the number @var{x} to its inexact representation.\n"
5973 scm_exact_to_inexact (SCM z
)
5975 if (SCM_I_INUMP (z
))
5976 return scm_from_double ((double) SCM_I_INUM (z
));
5977 else if (SCM_BIGP (z
))
5978 return scm_from_double (scm_i_big2dbl (z
));
5979 else if (SCM_FRACTIONP (z
))
5980 return scm_from_double (scm_i_fraction2double (z
));
5981 else if (SCM_INEXACTP (z
))
5984 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5988 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5990 "Return an exact number that is numerically closest to @var{z}.")
5991 #define FUNC_NAME s_scm_inexact_to_exact
5993 if (SCM_I_INUMP (z
))
5995 else if (SCM_BIGP (z
))
5997 else if (SCM_REALP (z
))
5999 if (xisinf (SCM_REAL_VALUE (z
)) || xisnan (SCM_REAL_VALUE (z
)))
6000 SCM_OUT_OF_RANGE (1, z
);
6007 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
6008 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
6009 scm_i_mpz2num (mpq_denref (frac
)));
6011 /* When scm_i_make_ratio throws, we leak the memory allocated
6018 else if (SCM_FRACTIONP (z
))
6021 SCM_WRONG_TYPE_ARG (1, z
);
6025 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
6027 "Returns the @emph{simplest} rational number differing\n"
6028 "from @var{x} by no more than @var{eps}.\n"
6030 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
6031 "exact result when both its arguments are exact. Thus, you might need\n"
6032 "to use @code{inexact->exact} on the arguments.\n"
6035 "(rationalize (inexact->exact 1.2) 1/100)\n"
6038 #define FUNC_NAME s_scm_rationalize
6040 if (SCM_I_INUMP (x
))
6042 else if (SCM_BIGP (x
))
6044 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
6046 /* Use continued fractions to find closest ratio. All
6047 arithmetic is done with exact numbers.
6050 SCM ex
= scm_inexact_to_exact (x
);
6051 SCM int_part
= scm_floor (ex
);
6052 SCM tt
= SCM_I_MAKINUM (1);
6053 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
6054 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
6058 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
6061 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
6062 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
6064 /* We stop after a million iterations just to be absolutely sure
6065 that we don't go into an infinite loop. The process normally
6066 converges after less than a dozen iterations.
6069 eps
= scm_abs (eps
);
6070 while (++i
< 1000000)
6072 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
6073 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
6074 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
6076 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
6077 eps
))) /* abs(x-a/b) <= eps */
6079 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
6080 if (scm_is_false (scm_exact_p (x
))
6081 || scm_is_false (scm_exact_p (eps
)))
6082 return scm_exact_to_inexact (res
);
6086 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
6088 tt
= scm_floor (rx
); /* tt = floor (rx) */
6094 scm_num_overflow (s_scm_rationalize
);
6097 SCM_WRONG_TYPE_ARG (1, x
);
6101 /* conversion functions */
6104 scm_is_integer (SCM val
)
6106 return scm_is_true (scm_integer_p (val
));
6110 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
6112 if (SCM_I_INUMP (val
))
6114 scm_t_signed_bits n
= SCM_I_INUM (val
);
6115 return n
>= min
&& n
<= max
;
6117 else if (SCM_BIGP (val
))
6119 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
6121 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
6123 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
6125 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
6126 return n
>= min
&& n
<= max
;
6136 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6137 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6140 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6141 SCM_I_BIG_MPZ (val
));
6143 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
6155 return n
>= min
&& n
<= max
;
6163 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
6165 if (SCM_I_INUMP (val
))
6167 scm_t_signed_bits n
= SCM_I_INUM (val
);
6168 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
6170 else if (SCM_BIGP (val
))
6172 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
6174 else if (max
<= ULONG_MAX
)
6176 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
6178 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
6179 return n
>= min
&& n
<= max
;
6189 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
6192 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6193 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6196 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6197 SCM_I_BIG_MPZ (val
));
6199 return n
>= min
&& n
<= max
;
6207 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
6209 scm_error (scm_out_of_range_key
,
6211 "Value out of range ~S to ~S: ~S",
6212 scm_list_3 (min
, max
, bad_val
),
6213 scm_list_1 (bad_val
));
6216 #define TYPE scm_t_intmax
6217 #define TYPE_MIN min
6218 #define TYPE_MAX max
6219 #define SIZEOF_TYPE 0
6220 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
6221 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
6222 #include "libguile/conv-integer.i.c"
6224 #define TYPE scm_t_uintmax
6225 #define TYPE_MIN min
6226 #define TYPE_MAX max
6227 #define SIZEOF_TYPE 0
6228 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6229 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6230 #include "libguile/conv-uinteger.i.c"
6232 #define TYPE scm_t_int8
6233 #define TYPE_MIN SCM_T_INT8_MIN
6234 #define TYPE_MAX SCM_T_INT8_MAX
6235 #define SIZEOF_TYPE 1
6236 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6237 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6238 #include "libguile/conv-integer.i.c"
6240 #define TYPE scm_t_uint8
6242 #define TYPE_MAX SCM_T_UINT8_MAX
6243 #define SIZEOF_TYPE 1
6244 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6245 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6246 #include "libguile/conv-uinteger.i.c"
6248 #define TYPE scm_t_int16
6249 #define TYPE_MIN SCM_T_INT16_MIN
6250 #define TYPE_MAX SCM_T_INT16_MAX
6251 #define SIZEOF_TYPE 2
6252 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6253 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6254 #include "libguile/conv-integer.i.c"
6256 #define TYPE scm_t_uint16
6258 #define TYPE_MAX SCM_T_UINT16_MAX
6259 #define SIZEOF_TYPE 2
6260 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6261 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6262 #include "libguile/conv-uinteger.i.c"
6264 #define TYPE scm_t_int32
6265 #define TYPE_MIN SCM_T_INT32_MIN
6266 #define TYPE_MAX SCM_T_INT32_MAX
6267 #define SIZEOF_TYPE 4
6268 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6269 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6270 #include "libguile/conv-integer.i.c"
6272 #define TYPE scm_t_uint32
6274 #define TYPE_MAX SCM_T_UINT32_MAX
6275 #define SIZEOF_TYPE 4
6276 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6277 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6278 #include "libguile/conv-uinteger.i.c"
6280 #define TYPE scm_t_wchar
6281 #define TYPE_MIN (scm_t_int32)-1
6282 #define TYPE_MAX (scm_t_int32)0x10ffff
6283 #define SIZEOF_TYPE 4
6284 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6285 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6286 #include "libguile/conv-integer.i.c"
6288 #if SCM_HAVE_T_INT64
6290 #define TYPE scm_t_int64
6291 #define TYPE_MIN SCM_T_INT64_MIN
6292 #define TYPE_MAX SCM_T_INT64_MAX
6293 #define SIZEOF_TYPE 8
6294 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6295 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6296 #include "libguile/conv-integer.i.c"
6298 #define TYPE scm_t_uint64
6300 #define TYPE_MAX SCM_T_UINT64_MAX
6301 #define SIZEOF_TYPE 8
6302 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6303 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6304 #include "libguile/conv-uinteger.i.c"
6309 scm_to_mpz (SCM val
, mpz_t rop
)
6311 if (SCM_I_INUMP (val
))
6312 mpz_set_si (rop
, SCM_I_INUM (val
));
6313 else if (SCM_BIGP (val
))
6314 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6316 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6320 scm_from_mpz (mpz_t val
)
6322 return scm_i_mpz2num (val
);
6326 scm_is_real (SCM val
)
6328 return scm_is_true (scm_real_p (val
));
6332 scm_is_rational (SCM val
)
6334 return scm_is_true (scm_rational_p (val
));
6338 scm_to_double (SCM val
)
6340 if (SCM_I_INUMP (val
))
6341 return SCM_I_INUM (val
);
6342 else if (SCM_BIGP (val
))
6343 return scm_i_big2dbl (val
);
6344 else if (SCM_FRACTIONP (val
))
6345 return scm_i_fraction2double (val
);
6346 else if (SCM_REALP (val
))
6347 return SCM_REAL_VALUE (val
);
6349 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6353 scm_from_double (double val
)
6355 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
6356 SCM_REAL_VALUE (z
) = val
;
6360 #if SCM_ENABLE_DISCOURAGED == 1
6363 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
6367 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6371 scm_out_of_range (NULL
, num
);
6374 return scm_to_double (num
);
6378 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
6382 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6386 scm_out_of_range (NULL
, num
);
6389 return scm_to_double (num
);
6395 scm_is_complex (SCM val
)
6397 return scm_is_true (scm_complex_p (val
));
6401 scm_c_real_part (SCM z
)
6403 if (SCM_COMPLEXP (z
))
6404 return SCM_COMPLEX_REAL (z
);
6407 /* Use the scm_real_part to get proper error checking and
6410 return scm_to_double (scm_real_part (z
));
6415 scm_c_imag_part (SCM z
)
6417 if (SCM_COMPLEXP (z
))
6418 return SCM_COMPLEX_IMAG (z
);
6421 /* Use the scm_imag_part to get proper error checking and
6422 dispatching. The result will almost always be 0.0, but not
6425 return scm_to_double (scm_imag_part (z
));
6430 scm_c_magnitude (SCM z
)
6432 return scm_to_double (scm_magnitude (z
));
6438 return scm_to_double (scm_angle (z
));
6442 scm_is_number (SCM z
)
6444 return scm_is_true (scm_number_p (z
));
6448 /* In the following functions we dispatch to the real-arg funcs like log()
6449 when we know the arg is real, instead of just handing everything to
6450 clog() for instance. This is in case clog() doesn't optimize for a
6451 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6452 well use it to go straight to the applicable C func. */
6454 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6456 "Return the natural logarithm of @var{z}.")
6457 #define FUNC_NAME s_scm_log
6459 if (SCM_COMPLEXP (z
))
6461 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6462 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6464 double re
= SCM_COMPLEX_REAL (z
);
6465 double im
= SCM_COMPLEX_IMAG (z
);
6466 return scm_c_make_rectangular (log (hypot (re
, im
)),
6472 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6473 although the value itself overflows. */
6474 double re
= scm_to_double (z
);
6475 double l
= log (fabs (re
));
6477 return scm_from_double (l
);
6479 return scm_c_make_rectangular (l
, M_PI
);
6485 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6487 "Return the base 10 logarithm of @var{z}.")
6488 #define FUNC_NAME s_scm_log10
6490 if (SCM_COMPLEXP (z
))
6492 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6493 clog() and a multiply by M_LOG10E, rather than the fallback
6494 log10+hypot+atan2.) */
6495 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG10 && defined (SCM_COMPLEX_VALUE)
6496 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6498 double re
= SCM_COMPLEX_REAL (z
);
6499 double im
= SCM_COMPLEX_IMAG (z
);
6500 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6501 M_LOG10E
* atan2 (im
, re
));
6506 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6507 although the value itself overflows. */
6508 double re
= scm_to_double (z
);
6509 double l
= log10 (fabs (re
));
6511 return scm_from_double (l
);
6513 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6519 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6521 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6522 "base of natural logarithms (2.71828@dots{}).")
6523 #define FUNC_NAME s_scm_exp
6525 if (SCM_COMPLEXP (z
))
6527 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6528 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6530 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6531 SCM_COMPLEX_IMAG (z
));
6536 /* When z is a negative bignum the conversion to double overflows,
6537 giving -infinity, but that's ok, the exp is still 0.0. */
6538 return scm_from_double (exp (scm_to_double (z
)));
6544 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6546 "Return the square root of @var{z}. Of the two possible roots\n"
6547 "(positive and negative), the one with the a positive real part\n"
6548 "is returned, or if that's zero then a positive imaginary part.\n"
6552 "(sqrt 9.0) @result{} 3.0\n"
6553 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6554 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6555 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6557 #define FUNC_NAME s_scm_sqrt
6559 if (SCM_COMPLEXP (x
))
6561 #if HAVE_COMPLEX_DOUBLE && HAVE_USABLE_CSQRT && defined (SCM_COMPLEX_VALUE)
6562 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6564 double re
= SCM_COMPLEX_REAL (x
);
6565 double im
= SCM_COMPLEX_IMAG (x
);
6566 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6567 0.5 * atan2 (im
, re
));
6572 double xx
= scm_to_double (x
);
6574 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6576 return scm_from_double (sqrt (xx
));
6588 mpz_init_set_si (z_negative_one
, -1);
6590 /* It may be possible to tune the performance of some algorithms by using
6591 * the following constants to avoid the creation of bignums. Please, before
6592 * using these values, remember the two rules of program optimization:
6593 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6594 scm_c_define ("most-positive-fixnum",
6595 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6596 scm_c_define ("most-negative-fixnum",
6597 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6599 scm_add_feature ("complex");
6600 scm_add_feature ("inexact");
6601 flo0
= scm_from_double (0.0);
6603 /* determine floating point precision */
6604 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6606 init_dblprec(&scm_dblprec
[i
-2],i
);
6607 init_fx_radix(fx_per_radix
[i
-2],i
);
6610 /* hard code precision for base 10 if the preprocessor tells us to... */
6611 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6614 exactly_one_half
= scm_divide (SCM_I_MAKINUM (1), SCM_I_MAKINUM (2));
6615 #include "libguile/numbers.x"