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 /* values per glibc, if not already defined */
72 #define M_LOG10E 0.43429448190325182765
75 #define M_PI 3.14159265358979323846
81 Wonder if this might be faster for some of our code? A switch on
82 the numtag would jump directly to the right case, and the
83 SCM_I_NUMTAG code might be faster than repeated SCM_FOOP tests...
85 #define SCM_I_NUMTAG_NOTNUM 0
86 #define SCM_I_NUMTAG_INUM 1
87 #define SCM_I_NUMTAG_BIG scm_tc16_big
88 #define SCM_I_NUMTAG_REAL scm_tc16_real
89 #define SCM_I_NUMTAG_COMPLEX scm_tc16_complex
90 #define SCM_I_NUMTAG(x) \
91 (SCM_I_INUMP(x) ? SCM_I_NUMTAG_INUM \
92 : (SCM_IMP(x) ? SCM_I_NUMTAG_NOTNUM \
93 : (((0xfcff & SCM_CELL_TYPE (x)) == scm_tc7_number) ? SCM_TYP16(x) \
94 : SCM_I_NUMTAG_NOTNUM)))
96 /* the macro above will not work as is with fractions */
101 #define SCM_SWAP(x, y) do { SCM __t = x; x = y; y = __t; } while (0)
103 /* FLOBUFLEN is the maximum number of characters neccessary for the
104 * printed or scm_string representation of an inexact number.
106 #define FLOBUFLEN (40+2*(sizeof(double)/sizeof(char)*SCM_CHAR_BIT*3+9)/10)
109 #if !defined (HAVE_ASINH)
110 static double asinh (double x
) { return log (x
+ sqrt (x
* x
+ 1)); }
112 #if !defined (HAVE_ACOSH)
113 static double acosh (double x
) { return log (x
+ sqrt (x
* x
- 1)); }
115 #if !defined (HAVE_ATANH)
116 static double atanh (double x
) { return 0.5 * log ((1 + x
) / (1 - x
)); }
119 /* mpz_cmp_d in gmp 4.1.3 doesn't recognise infinities, so xmpz_cmp_d uses
120 an explicit check. In some future gmp (don't know what version number),
121 mpz_cmp_d is supposed to do this itself. */
123 #define xmpz_cmp_d(z, d) \
124 (isinf (d) ? (d < 0.0 ? 1 : -1) : mpz_cmp_d (z, d))
126 #define xmpz_cmp_d(z, d) mpz_cmp_d (z, d)
130 #if defined (GUILE_I)
131 #if HAVE_COMPLEX_DOUBLE
133 /* For an SCM object Z which is a complex number (ie. satisfies
134 SCM_COMPLEXP), return its value as a C level "complex double". */
135 #define SCM_COMPLEX_VALUE(z) \
136 (SCM_COMPLEX_REAL (z) + GUILE_I * SCM_COMPLEX_IMAG (z))
138 static inline SCM
scm_from_complex_double (complex double z
) SCM_UNUSED
;
140 /* Convert a C "complex double" to an SCM value. */
142 scm_from_complex_double (complex double z
)
144 return scm_c_make_rectangular (creal (z
), cimag (z
));
147 #endif /* HAVE_COMPLEX_DOUBLE */
152 static mpz_t z_negative_one
;
159 /* Return a newly created bignum. */
160 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
161 mpz_init (SCM_I_BIG_MPZ (z
));
166 scm_i_long2big (long x
)
168 /* Return a newly created bignum initialized to X. */
169 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
170 mpz_init_set_si (SCM_I_BIG_MPZ (z
), x
);
175 scm_i_ulong2big (unsigned long x
)
177 /* Return a newly created bignum initialized to X. */
178 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
179 mpz_init_set_ui (SCM_I_BIG_MPZ (z
), x
);
184 scm_i_clonebig (SCM src_big
, int same_sign_p
)
186 /* Copy src_big's value, negate it if same_sign_p is false, and return. */
187 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
188 mpz_init_set (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (src_big
));
190 mpz_neg (SCM_I_BIG_MPZ (z
), SCM_I_BIG_MPZ (z
));
195 scm_i_bigcmp (SCM x
, SCM y
)
197 /* Return neg if x < y, pos if x > y, and 0 if x == y */
198 /* presume we already know x and y are bignums */
199 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
200 scm_remember_upto_here_2 (x
, y
);
205 scm_i_dbl2big (double d
)
207 /* results are only defined if d is an integer */
208 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
209 mpz_init_set_d (SCM_I_BIG_MPZ (z
), d
);
213 /* Convert a integer in double representation to a SCM number. */
216 scm_i_dbl2num (double u
)
218 /* SCM_MOST_POSITIVE_FIXNUM+1 and SCM_MOST_NEGATIVE_FIXNUM are both
219 powers of 2, so there's no rounding when making "double" values
220 from them. If plain SCM_MOST_POSITIVE_FIXNUM was used it could
221 get rounded on a 64-bit machine, hence the "+1".
223 The use of floor() to force to an integer value ensures we get a
224 "numerically closest" value without depending on how a
225 double->long cast or how mpz_set_d will round. For reference,
226 double->long probably follows the hardware rounding mode,
227 mpz_set_d truncates towards zero. */
229 /* XXX - what happens when SCM_MOST_POSITIVE_FIXNUM etc is not
230 representable as a double? */
232 if (u
< (double) (SCM_MOST_POSITIVE_FIXNUM
+1)
233 && u
>= (double) SCM_MOST_NEGATIVE_FIXNUM
)
234 return SCM_I_MAKINUM ((long) u
);
236 return scm_i_dbl2big (u
);
239 /* scm_i_big2dbl() rounds to the closest representable double, in accordance
240 with R5RS exact->inexact.
242 The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
243 (ie. truncate towards zero), then adjust to get the closest double by
244 examining the next lower bit and adding 1 (to the absolute value) if
247 Bignums exactly half way between representable doubles are rounded to the
248 next higher absolute value (ie. away from zero). This seems like an
249 adequate interpretation of R5RS "numerically closest", and it's easier
250 and faster than a full "nearest-even" style.
252 The bit test must be done on the absolute value of the mpz_t, which means
253 we need to use mpz_getlimbn. mpz_tstbit is not right, it treats
254 negatives as twos complement.
256 In current gmp 4.1.3, mpz_get_d rounding is unspecified. It ends up
257 following the hardware rounding mode, but applied to the absolute value
258 of the mpz_t operand. This is not what we want so we put the high
259 DBL_MANT_DIG bits into a temporary. In some future gmp, don't know when,
260 mpz_get_d is supposed to always truncate towards zero.
262 ENHANCE-ME: The temporary init+clear to force the rounding in gmp 4.1.3
263 is a slowdown. It'd be faster to pick out the relevant high bits with
264 mpz_getlimbn if we could be bothered coding that, and if the new
265 truncating gmp doesn't come out. */
268 scm_i_big2dbl (SCM b
)
273 bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (b
), 2);
277 /* Current GMP, eg. 4.1.3, force truncation towards zero */
279 if (bits
> DBL_MANT_DIG
)
281 size_t shift
= bits
- DBL_MANT_DIG
;
282 mpz_init2 (tmp
, DBL_MANT_DIG
);
283 mpz_tdiv_q_2exp (tmp
, SCM_I_BIG_MPZ (b
), shift
);
284 result
= ldexp (mpz_get_d (tmp
), shift
);
289 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
294 result
= mpz_get_d (SCM_I_BIG_MPZ (b
));
297 if (bits
> DBL_MANT_DIG
)
299 unsigned long pos
= bits
- DBL_MANT_DIG
- 1;
300 /* test bit number "pos" in absolute value */
301 if (mpz_getlimbn (SCM_I_BIG_MPZ (b
), pos
/ GMP_NUMB_BITS
)
302 & ((mp_limb_t
) 1 << (pos
% GMP_NUMB_BITS
)))
304 result
+= ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b
)), pos
+ 1);
308 scm_remember_upto_here_1 (b
);
313 scm_i_normbig (SCM b
)
315 /* convert a big back to a fixnum if it'll fit */
316 /* presume b is a bignum */
317 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (b
)))
319 long val
= mpz_get_si (SCM_I_BIG_MPZ (b
));
320 if (SCM_FIXABLE (val
))
321 b
= SCM_I_MAKINUM (val
);
326 static SCM_C_INLINE_KEYWORD SCM
327 scm_i_mpz2num (mpz_t b
)
329 /* convert a mpz number to a SCM number. */
330 if (mpz_fits_slong_p (b
))
332 long val
= mpz_get_si (b
);
333 if (SCM_FIXABLE (val
))
334 return SCM_I_MAKINUM (val
);
338 SCM z
= scm_double_cell (scm_tc16_big
, 0, 0, 0);
339 mpz_init_set (SCM_I_BIG_MPZ (z
), b
);
344 /* this is needed when we want scm_divide to make a float, not a ratio, even if passed two ints */
345 static SCM
scm_divide2real (SCM x
, SCM y
);
348 scm_i_make_ratio (SCM numerator
, SCM denominator
)
349 #define FUNC_NAME "make-ratio"
351 /* First make sure the arguments are proper.
353 if (SCM_I_INUMP (denominator
))
355 if (scm_is_eq (denominator
, SCM_INUM0
))
356 scm_num_overflow ("make-ratio");
357 if (scm_is_eq (denominator
, SCM_I_MAKINUM(1)))
362 if (!(SCM_BIGP(denominator
)))
363 SCM_WRONG_TYPE_ARG (2, denominator
);
365 if (!SCM_I_INUMP (numerator
) && !SCM_BIGP (numerator
))
366 SCM_WRONG_TYPE_ARG (1, numerator
);
368 /* Then flip signs so that the denominator is positive.
370 if (scm_is_true (scm_negative_p (denominator
)))
372 numerator
= scm_difference (numerator
, SCM_UNDEFINED
);
373 denominator
= scm_difference (denominator
, SCM_UNDEFINED
);
376 /* Now consider for each of the four fixnum/bignum combinations
377 whether the rational number is really an integer.
379 if (SCM_I_INUMP (numerator
))
381 long x
= SCM_I_INUM (numerator
);
382 if (scm_is_eq (numerator
, SCM_INUM0
))
384 if (SCM_I_INUMP (denominator
))
387 y
= SCM_I_INUM (denominator
);
389 return SCM_I_MAKINUM(1);
391 return SCM_I_MAKINUM (x
/ y
);
395 /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
396 of that value for the denominator, as a bignum. Apart from
397 that case, abs(bignum) > abs(inum) so inum/bignum is not an
399 if (x
== SCM_MOST_NEGATIVE_FIXNUM
400 && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator
),
401 - SCM_MOST_NEGATIVE_FIXNUM
) == 0)
402 return SCM_I_MAKINUM(-1);
405 else if (SCM_BIGP (numerator
))
407 if (SCM_I_INUMP (denominator
))
409 long yy
= SCM_I_INUM (denominator
);
410 if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator
), yy
))
411 return scm_divide (numerator
, denominator
);
415 if (scm_is_eq (numerator
, denominator
))
416 return SCM_I_MAKINUM(1);
417 if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator
),
418 SCM_I_BIG_MPZ (denominator
)))
419 return scm_divide(numerator
, denominator
);
423 /* No, it's a proper fraction.
426 SCM divisor
= scm_gcd (numerator
, denominator
);
427 if (!(scm_is_eq (divisor
, SCM_I_MAKINUM(1))))
429 numerator
= scm_divide (numerator
, divisor
);
430 denominator
= scm_divide (denominator
, divisor
);
433 return scm_double_cell (scm_tc16_fraction
,
434 SCM_UNPACK (numerator
),
435 SCM_UNPACK (denominator
), 0);
441 scm_i_fraction2double (SCM z
)
443 return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z
),
444 SCM_FRACTION_DENOMINATOR (z
)));
447 SCM_DEFINE (scm_exact_p
, "exact?", 1, 0, 0,
449 "Return @code{#t} if @var{x} is an exact number, @code{#f}\n"
451 #define FUNC_NAME s_scm_exact_p
457 if (SCM_FRACTIONP (x
))
461 SCM_WRONG_TYPE_ARG (1, x
);
466 SCM_DEFINE (scm_odd_p
, "odd?", 1, 0, 0,
468 "Return @code{#t} if @var{n} is an odd number, @code{#f}\n"
470 #define FUNC_NAME s_scm_odd_p
474 long val
= SCM_I_INUM (n
);
475 return scm_from_bool ((val
& 1L) != 0);
477 else if (SCM_BIGP (n
))
479 int odd_p
= mpz_odd_p (SCM_I_BIG_MPZ (n
));
480 scm_remember_upto_here_1 (n
);
481 return scm_from_bool (odd_p
);
483 else if (scm_is_true (scm_inf_p (n
)))
485 else if (SCM_REALP (n
))
487 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
493 SCM_WRONG_TYPE_ARG (1, n
);
496 SCM_WRONG_TYPE_ARG (1, n
);
501 SCM_DEFINE (scm_even_p
, "even?", 1, 0, 0,
503 "Return @code{#t} if @var{n} is an even number, @code{#f}\n"
505 #define FUNC_NAME s_scm_even_p
509 long val
= SCM_I_INUM (n
);
510 return scm_from_bool ((val
& 1L) == 0);
512 else if (SCM_BIGP (n
))
514 int even_p
= mpz_even_p (SCM_I_BIG_MPZ (n
));
515 scm_remember_upto_here_1 (n
);
516 return scm_from_bool (even_p
);
518 else if (scm_is_true (scm_inf_p (n
)))
520 else if (SCM_REALP (n
))
522 double rem
= fabs (fmod (SCM_REAL_VALUE(n
), 2.0));
528 SCM_WRONG_TYPE_ARG (1, n
);
531 SCM_WRONG_TYPE_ARG (1, n
);
535 SCM_DEFINE (scm_inf_p
, "inf?", 1, 0, 0,
537 "Return @code{#t} if @var{x} is either @samp{+inf.0}\n"
538 "or @samp{-inf.0}, @code{#f} otherwise.")
539 #define FUNC_NAME s_scm_inf_p
542 return scm_from_bool (isinf (SCM_REAL_VALUE (x
)));
543 else if (SCM_COMPLEXP (x
))
544 return scm_from_bool (isinf (SCM_COMPLEX_REAL (x
))
545 || isinf (SCM_COMPLEX_IMAG (x
)));
551 SCM_DEFINE (scm_nan_p
, "nan?", 1, 0, 0,
553 "Return @code{#t} if @var{n} is a NaN, @code{#f}\n"
555 #define FUNC_NAME s_scm_nan_p
558 return scm_from_bool (isnan (SCM_REAL_VALUE (n
)));
559 else if (SCM_COMPLEXP (n
))
560 return scm_from_bool (isnan (SCM_COMPLEX_REAL (n
))
561 || isnan (SCM_COMPLEX_IMAG (n
)));
567 /* Guile's idea of infinity. */
568 static double guile_Inf
;
570 /* Guile's idea of not a number. */
571 static double guile_NaN
;
574 guile_ieee_init (void)
576 /* Some version of gcc on some old version of Linux used to crash when
577 trying to make Inf and NaN. */
580 /* C99 INFINITY, when available.
581 FIXME: The standard allows for INFINITY to be something that overflows
582 at compile time. We ought to have a configure test to check for that
583 before trying to use it. (But in practice we believe this is not a
584 problem on any system guile is likely to target.) */
585 guile_Inf
= INFINITY
;
586 #elif defined HAVE_DINFINITY
588 extern unsigned int DINFINITY
[2];
589 guile_Inf
= (*((double *) (DINFINITY
)));
596 if (guile_Inf
== tmp
)
603 /* C99 NAN, when available */
605 #elif defined HAVE_DQNAN
608 extern unsigned int DQNAN
[2];
609 guile_NaN
= (*((double *)(DQNAN
)));
612 guile_NaN
= guile_Inf
/ guile_Inf
;
616 SCM_DEFINE (scm_inf
, "inf", 0, 0, 0,
619 #define FUNC_NAME s_scm_inf
621 static int initialized
= 0;
627 return scm_from_double (guile_Inf
);
631 SCM_DEFINE (scm_nan
, "nan", 0, 0, 0,
634 #define FUNC_NAME s_scm_nan
636 static int initialized
= 0;
642 return scm_from_double (guile_NaN
);
647 SCM_PRIMITIVE_GENERIC (scm_abs
, "abs", 1, 0, 0,
649 "Return the absolute value of @var{x}.")
654 long int xx
= SCM_I_INUM (x
);
657 else if (SCM_POSFIXABLE (-xx
))
658 return SCM_I_MAKINUM (-xx
);
660 return scm_i_long2big (-xx
);
662 else if (SCM_BIGP (x
))
664 const int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
666 return scm_i_clonebig (x
, 0);
670 else if (SCM_REALP (x
))
672 /* note that if x is a NaN then xx<0 is false so we return x unchanged */
673 double xx
= SCM_REAL_VALUE (x
);
675 return scm_from_double (-xx
);
679 else if (SCM_FRACTIONP (x
))
681 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x
))))
683 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
684 SCM_FRACTION_DENOMINATOR (x
));
687 SCM_WTA_DISPATCH_1 (g_scm_abs
, x
, 1, s_scm_abs
);
692 SCM_GPROC (s_quotient
, "quotient", 2, 0, 0, scm_quotient
, g_quotient
);
693 /* "Return the quotient of the numbers @var{x} and @var{y}."
696 scm_quotient (SCM x
, SCM y
)
700 long xx
= SCM_I_INUM (x
);
703 long yy
= SCM_I_INUM (y
);
705 scm_num_overflow (s_quotient
);
710 return SCM_I_MAKINUM (z
);
712 return scm_i_long2big (z
);
715 else if (SCM_BIGP (y
))
717 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
718 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
719 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
721 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
722 scm_remember_upto_here_1 (y
);
723 return SCM_I_MAKINUM (-1);
726 return SCM_I_MAKINUM (0);
729 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
731 else if (SCM_BIGP (x
))
735 long yy
= SCM_I_INUM (y
);
737 scm_num_overflow (s_quotient
);
742 SCM result
= scm_i_mkbig ();
745 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
),
748 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
751 mpz_tdiv_q_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
752 scm_remember_upto_here_1 (x
);
753 return scm_i_normbig (result
);
756 else if (SCM_BIGP (y
))
758 SCM result
= scm_i_mkbig ();
759 mpz_tdiv_q (SCM_I_BIG_MPZ (result
),
762 scm_remember_upto_here_2 (x
, y
);
763 return scm_i_normbig (result
);
766 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG2
, s_quotient
);
769 SCM_WTA_DISPATCH_2 (g_quotient
, x
, y
, SCM_ARG1
, s_quotient
);
772 SCM_GPROC (s_remainder
, "remainder", 2, 0, 0, scm_remainder
, g_remainder
);
773 /* "Return the remainder of the numbers @var{x} and @var{y}.\n"
775 * "(remainder 13 4) @result{} 1\n"
776 * "(remainder -13 4) @result{} -1\n"
780 scm_remainder (SCM x
, SCM y
)
786 long yy
= SCM_I_INUM (y
);
788 scm_num_overflow (s_remainder
);
791 long z
= SCM_I_INUM (x
) % yy
;
792 return SCM_I_MAKINUM (z
);
795 else if (SCM_BIGP (y
))
797 if ((SCM_I_INUM (x
) == SCM_MOST_NEGATIVE_FIXNUM
)
798 && (mpz_cmp_ui (SCM_I_BIG_MPZ (y
),
799 - SCM_MOST_NEGATIVE_FIXNUM
) == 0))
801 /* Special case: x == fixnum-min && y == abs (fixnum-min) */
802 scm_remember_upto_here_1 (y
);
803 return SCM_I_MAKINUM (0);
809 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
811 else if (SCM_BIGP (x
))
815 long yy
= SCM_I_INUM (y
);
817 scm_num_overflow (s_remainder
);
820 SCM result
= scm_i_mkbig ();
823 mpz_tdiv_r_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ(x
), yy
);
824 scm_remember_upto_here_1 (x
);
825 return scm_i_normbig (result
);
828 else if (SCM_BIGP (y
))
830 SCM result
= scm_i_mkbig ();
831 mpz_tdiv_r (SCM_I_BIG_MPZ (result
),
834 scm_remember_upto_here_2 (x
, y
);
835 return scm_i_normbig (result
);
838 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG2
, s_remainder
);
841 SCM_WTA_DISPATCH_2 (g_remainder
, x
, y
, SCM_ARG1
, s_remainder
);
845 SCM_GPROC (s_modulo
, "modulo", 2, 0, 0, scm_modulo
, g_modulo
);
846 /* "Return the modulo of the numbers @var{x} and @var{y}.\n"
848 * "(modulo 13 4) @result{} 1\n"
849 * "(modulo -13 4) @result{} 3\n"
853 scm_modulo (SCM x
, SCM y
)
857 long xx
= SCM_I_INUM (x
);
860 long yy
= SCM_I_INUM (y
);
862 scm_num_overflow (s_modulo
);
865 /* C99 specifies that "%" is the remainder corresponding to a
866 quotient rounded towards zero, and that's also traditional
867 for machine division, so z here should be well defined. */
885 return SCM_I_MAKINUM (result
);
888 else if (SCM_BIGP (y
))
890 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
897 SCM pos_y
= scm_i_clonebig (y
, 0);
898 /* do this after the last scm_op */
899 mpz_init_set_si (z_x
, xx
);
900 result
= pos_y
; /* re-use this bignum */
901 mpz_mod (SCM_I_BIG_MPZ (result
),
903 SCM_I_BIG_MPZ (pos_y
));
904 scm_remember_upto_here_1 (pos_y
);
908 result
= scm_i_mkbig ();
909 /* do this after the last scm_op */
910 mpz_init_set_si (z_x
, xx
);
911 mpz_mod (SCM_I_BIG_MPZ (result
),
914 scm_remember_upto_here_1 (y
);
917 if ((sgn_y
< 0) && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
918 mpz_add (SCM_I_BIG_MPZ (result
),
920 SCM_I_BIG_MPZ (result
));
921 scm_remember_upto_here_1 (y
);
922 /* and do this before the next one */
924 return scm_i_normbig (result
);
928 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
930 else if (SCM_BIGP (x
))
934 long yy
= SCM_I_INUM (y
);
936 scm_num_overflow (s_modulo
);
939 SCM result
= scm_i_mkbig ();
940 mpz_mod_ui (SCM_I_BIG_MPZ (result
),
942 (yy
< 0) ? - yy
: yy
);
943 scm_remember_upto_here_1 (x
);
944 if ((yy
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
945 mpz_sub_ui (SCM_I_BIG_MPZ (result
),
946 SCM_I_BIG_MPZ (result
),
948 return scm_i_normbig (result
);
951 else if (SCM_BIGP (y
))
954 SCM result
= scm_i_mkbig ();
955 int y_sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
956 SCM pos_y
= scm_i_clonebig (y
, y_sgn
>= 0);
957 mpz_mod (SCM_I_BIG_MPZ (result
),
959 SCM_I_BIG_MPZ (pos_y
));
961 scm_remember_upto_here_1 (x
);
962 if ((y_sgn
< 0) && (mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0))
963 mpz_add (SCM_I_BIG_MPZ (result
),
965 SCM_I_BIG_MPZ (result
));
966 scm_remember_upto_here_2 (y
, pos_y
);
967 return scm_i_normbig (result
);
971 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG2
, s_modulo
);
974 SCM_WTA_DISPATCH_2 (g_modulo
, x
, y
, SCM_ARG1
, s_modulo
);
977 SCM_PRIMITIVE_GENERIC (scm_i_gcd
, "gcd", 0, 2, 1,
978 (SCM x
, SCM y
, SCM rest
),
979 "Return the greatest common divisor of all parameter values.\n"
980 "If called without arguments, 0 is returned.")
981 #define FUNC_NAME s_scm_i_gcd
983 while (!scm_is_null (rest
))
984 { x
= scm_gcd (x
, y
);
986 rest
= scm_cdr (rest
);
988 return scm_gcd (x
, y
);
992 #define s_gcd s_scm_i_gcd
993 #define g_gcd g_scm_i_gcd
996 scm_gcd (SCM x
, SCM y
)
999 return SCM_UNBNDP (x
) ? SCM_INUM0
: scm_abs (x
);
1001 if (SCM_I_INUMP (x
))
1003 if (SCM_I_INUMP (y
))
1005 long xx
= SCM_I_INUM (x
);
1006 long yy
= SCM_I_INUM (y
);
1007 long u
= xx
< 0 ? -xx
: xx
;
1008 long v
= yy
< 0 ? -yy
: yy
;
1018 /* Determine a common factor 2^k */
1019 while (!(1 & (u
| v
)))
1025 /* Now, any factor 2^n can be eliminated */
1045 return (SCM_POSFIXABLE (result
)
1046 ? SCM_I_MAKINUM (result
)
1047 : scm_i_long2big (result
));
1049 else if (SCM_BIGP (y
))
1055 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1057 else if (SCM_BIGP (x
))
1059 if (SCM_I_INUMP (y
))
1061 unsigned long result
;
1064 yy
= SCM_I_INUM (y
);
1069 result
= mpz_gcd_ui (NULL
, SCM_I_BIG_MPZ (x
), yy
);
1070 scm_remember_upto_here_1 (x
);
1071 return (SCM_POSFIXABLE (result
)
1072 ? SCM_I_MAKINUM (result
)
1073 : scm_from_ulong (result
));
1075 else if (SCM_BIGP (y
))
1077 SCM result
= scm_i_mkbig ();
1078 mpz_gcd (SCM_I_BIG_MPZ (result
),
1081 scm_remember_upto_here_2 (x
, y
);
1082 return scm_i_normbig (result
);
1085 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG2
, s_gcd
);
1088 SCM_WTA_DISPATCH_2 (g_gcd
, x
, y
, SCM_ARG1
, s_gcd
);
1091 SCM_PRIMITIVE_GENERIC (scm_i_lcm
, "lcm", 0, 2, 1,
1092 (SCM x
, SCM y
, SCM rest
),
1093 "Return the least common multiple of the arguments.\n"
1094 "If called without arguments, 1 is returned.")
1095 #define FUNC_NAME s_scm_i_lcm
1097 while (!scm_is_null (rest
))
1098 { x
= scm_lcm (x
, y
);
1100 rest
= scm_cdr (rest
);
1102 return scm_lcm (x
, y
);
1106 #define s_lcm s_scm_i_lcm
1107 #define g_lcm g_scm_i_lcm
1110 scm_lcm (SCM n1
, SCM n2
)
1112 if (SCM_UNBNDP (n2
))
1114 if (SCM_UNBNDP (n1
))
1115 return SCM_I_MAKINUM (1L);
1116 n2
= SCM_I_MAKINUM (1L);
1119 SCM_GASSERT2 (SCM_I_INUMP (n1
) || SCM_BIGP (n1
),
1120 g_lcm
, n1
, n2
, SCM_ARG1
, s_lcm
);
1121 SCM_GASSERT2 (SCM_I_INUMP (n2
) || SCM_BIGP (n2
),
1122 g_lcm
, n1
, n2
, SCM_ARGn
, s_lcm
);
1124 if (SCM_I_INUMP (n1
))
1126 if (SCM_I_INUMP (n2
))
1128 SCM d
= scm_gcd (n1
, n2
);
1129 if (scm_is_eq (d
, SCM_INUM0
))
1132 return scm_abs (scm_product (n1
, scm_quotient (n2
, d
)));
1136 /* inum n1, big n2 */
1139 SCM result
= scm_i_mkbig ();
1140 long nn1
= SCM_I_INUM (n1
);
1141 if (nn1
== 0) return SCM_INUM0
;
1142 if (nn1
< 0) nn1
= - nn1
;
1143 mpz_lcm_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n2
), nn1
);
1144 scm_remember_upto_here_1 (n2
);
1152 if (SCM_I_INUMP (n2
))
1159 SCM result
= scm_i_mkbig ();
1160 mpz_lcm(SCM_I_BIG_MPZ (result
),
1162 SCM_I_BIG_MPZ (n2
));
1163 scm_remember_upto_here_2(n1
, n2
);
1164 /* shouldn't need to normalize b/c lcm of 2 bigs should be big */
1170 /* Emulating 2's complement bignums with sign magnitude arithmetic:
1175 + + + x (map digit:logand X Y)
1176 + - + x (map digit:logand X (lognot (+ -1 Y)))
1177 - + + y (map digit:logand (lognot (+ -1 X)) Y)
1178 - - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
1183 + + + (map digit:logior X Y)
1184 + - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
1185 - + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
1186 - - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
1191 + + + (map digit:logxor X Y)
1192 + - - (+ 1 (map digit:logxor X (+ -1 Y)))
1193 - + - (+ 1 (map digit:logxor (+ -1 X) Y))
1194 - - + (map digit:logxor (+ -1 X) (+ -1 Y))
1199 + + (any digit:logand X Y)
1200 + - (any digit:logand X (lognot (+ -1 Y)))
1201 - + (any digit:logand (lognot (+ -1 X)) Y)
1206 SCM_DEFINE (scm_i_logand
, "logand", 0, 2, 1,
1207 (SCM x
, SCM y
, SCM rest
),
1208 "Return the bitwise AND of the integer arguments.\n\n"
1210 "(logand) @result{} -1\n"
1211 "(logand 7) @result{} 7\n"
1212 "(logand #b111 #b011 #b001) @result{} 1\n"
1214 #define FUNC_NAME s_scm_i_logand
1216 while (!scm_is_null (rest
))
1217 { x
= scm_logand (x
, y
);
1219 rest
= scm_cdr (rest
);
1221 return scm_logand (x
, y
);
1225 #define s_scm_logand s_scm_i_logand
1227 SCM
scm_logand (SCM n1
, SCM n2
)
1228 #define FUNC_NAME s_scm_logand
1232 if (SCM_UNBNDP (n2
))
1234 if (SCM_UNBNDP (n1
))
1235 return SCM_I_MAKINUM (-1);
1236 else if (!SCM_NUMBERP (n1
))
1237 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1238 else if (SCM_NUMBERP (n1
))
1241 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1244 if (SCM_I_INUMP (n1
))
1246 nn1
= SCM_I_INUM (n1
);
1247 if (SCM_I_INUMP (n2
))
1249 long nn2
= SCM_I_INUM (n2
);
1250 return SCM_I_MAKINUM (nn1
& nn2
);
1252 else if SCM_BIGP (n2
)
1258 SCM result_z
= scm_i_mkbig ();
1260 mpz_init_set_si (nn1_z
, nn1
);
1261 mpz_and (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1262 scm_remember_upto_here_1 (n2
);
1264 return scm_i_normbig (result_z
);
1268 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1270 else if (SCM_BIGP (n1
))
1272 if (SCM_I_INUMP (n2
))
1275 nn1
= SCM_I_INUM (n1
);
1278 else if (SCM_BIGP (n2
))
1280 SCM result_z
= scm_i_mkbig ();
1281 mpz_and (SCM_I_BIG_MPZ (result_z
),
1283 SCM_I_BIG_MPZ (n2
));
1284 scm_remember_upto_here_2 (n1
, n2
);
1285 return scm_i_normbig (result_z
);
1288 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1291 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1296 SCM_DEFINE (scm_i_logior
, "logior", 0, 2, 1,
1297 (SCM x
, SCM y
, SCM rest
),
1298 "Return the bitwise OR of the integer arguments.\n\n"
1300 "(logior) @result{} 0\n"
1301 "(logior 7) @result{} 7\n"
1302 "(logior #b000 #b001 #b011) @result{} 3\n"
1304 #define FUNC_NAME s_scm_i_logior
1306 while (!scm_is_null (rest
))
1307 { x
= scm_logior (x
, y
);
1309 rest
= scm_cdr (rest
);
1311 return scm_logior (x
, y
);
1315 #define s_scm_logior s_scm_i_logior
1317 SCM
scm_logior (SCM n1
, SCM n2
)
1318 #define FUNC_NAME s_scm_logior
1322 if (SCM_UNBNDP (n2
))
1324 if (SCM_UNBNDP (n1
))
1326 else if (SCM_NUMBERP (n1
))
1329 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1332 if (SCM_I_INUMP (n1
))
1334 nn1
= SCM_I_INUM (n1
);
1335 if (SCM_I_INUMP (n2
))
1337 long nn2
= SCM_I_INUM (n2
);
1338 return SCM_I_MAKINUM (nn1
| nn2
);
1340 else if (SCM_BIGP (n2
))
1346 SCM result_z
= scm_i_mkbig ();
1348 mpz_init_set_si (nn1_z
, nn1
);
1349 mpz_ior (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1350 scm_remember_upto_here_1 (n2
);
1352 return scm_i_normbig (result_z
);
1356 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1358 else if (SCM_BIGP (n1
))
1360 if (SCM_I_INUMP (n2
))
1363 nn1
= SCM_I_INUM (n1
);
1366 else if (SCM_BIGP (n2
))
1368 SCM result_z
= scm_i_mkbig ();
1369 mpz_ior (SCM_I_BIG_MPZ (result_z
),
1371 SCM_I_BIG_MPZ (n2
));
1372 scm_remember_upto_here_2 (n1
, n2
);
1373 return scm_i_normbig (result_z
);
1376 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1379 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1384 SCM_DEFINE (scm_i_logxor
, "logxor", 0, 2, 1,
1385 (SCM x
, SCM y
, SCM rest
),
1386 "Return the bitwise XOR of the integer arguments. A bit is\n"
1387 "set in the result if it is set in an odd number of arguments.\n"
1389 "(logxor) @result{} 0\n"
1390 "(logxor 7) @result{} 7\n"
1391 "(logxor #b000 #b001 #b011) @result{} 2\n"
1392 "(logxor #b000 #b001 #b011 #b011) @result{} 1\n"
1394 #define FUNC_NAME s_scm_i_logxor
1396 while (!scm_is_null (rest
))
1397 { x
= scm_logxor (x
, y
);
1399 rest
= scm_cdr (rest
);
1401 return scm_logxor (x
, y
);
1405 #define s_scm_logxor s_scm_i_logxor
1407 SCM
scm_logxor (SCM n1
, SCM n2
)
1408 #define FUNC_NAME s_scm_logxor
1412 if (SCM_UNBNDP (n2
))
1414 if (SCM_UNBNDP (n1
))
1416 else if (SCM_NUMBERP (n1
))
1419 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1422 if (SCM_I_INUMP (n1
))
1424 nn1
= SCM_I_INUM (n1
);
1425 if (SCM_I_INUMP (n2
))
1427 long nn2
= SCM_I_INUM (n2
);
1428 return SCM_I_MAKINUM (nn1
^ nn2
);
1430 else if (SCM_BIGP (n2
))
1434 SCM result_z
= scm_i_mkbig ();
1436 mpz_init_set_si (nn1_z
, nn1
);
1437 mpz_xor (SCM_I_BIG_MPZ (result_z
), nn1_z
, SCM_I_BIG_MPZ (n2
));
1438 scm_remember_upto_here_1 (n2
);
1440 return scm_i_normbig (result_z
);
1444 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1446 else if (SCM_BIGP (n1
))
1448 if (SCM_I_INUMP (n2
))
1451 nn1
= SCM_I_INUM (n1
);
1454 else if (SCM_BIGP (n2
))
1456 SCM result_z
= scm_i_mkbig ();
1457 mpz_xor (SCM_I_BIG_MPZ (result_z
),
1459 SCM_I_BIG_MPZ (n2
));
1460 scm_remember_upto_here_2 (n1
, n2
);
1461 return scm_i_normbig (result_z
);
1464 SCM_WRONG_TYPE_ARG (SCM_ARG2
, n2
);
1467 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n1
);
1472 SCM_DEFINE (scm_logtest
, "logtest", 2, 0, 0,
1474 "Test whether @var{j} and @var{k} have any 1 bits in common.\n"
1475 "This is equivalent to @code{(not (zero? (logand j k)))}, but\n"
1476 "without actually calculating the @code{logand}, just testing\n"
1480 "(logtest #b0100 #b1011) @result{} #f\n"
1481 "(logtest #b0100 #b0111) @result{} #t\n"
1483 #define FUNC_NAME s_scm_logtest
1487 if (SCM_I_INUMP (j
))
1489 nj
= SCM_I_INUM (j
);
1490 if (SCM_I_INUMP (k
))
1492 long nk
= SCM_I_INUM (k
);
1493 return scm_from_bool (nj
& nk
);
1495 else if (SCM_BIGP (k
))
1503 mpz_init_set_si (nj_z
, nj
);
1504 mpz_and (nj_z
, nj_z
, SCM_I_BIG_MPZ (k
));
1505 scm_remember_upto_here_1 (k
);
1506 result
= scm_from_bool (mpz_sgn (nj_z
) != 0);
1512 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1514 else if (SCM_BIGP (j
))
1516 if (SCM_I_INUMP (k
))
1519 nj
= SCM_I_INUM (j
);
1522 else if (SCM_BIGP (k
))
1526 mpz_init (result_z
);
1530 scm_remember_upto_here_2 (j
, k
);
1531 result
= scm_from_bool (mpz_sgn (result_z
) != 0);
1532 mpz_clear (result_z
);
1536 SCM_WRONG_TYPE_ARG (SCM_ARG2
, k
);
1539 SCM_WRONG_TYPE_ARG (SCM_ARG1
, j
);
1544 SCM_DEFINE (scm_logbit_p
, "logbit?", 2, 0, 0,
1546 "Test whether bit number @var{index} in @var{j} is set.\n"
1547 "@var{index} starts from 0 for the least significant bit.\n"
1550 "(logbit? 0 #b1101) @result{} #t\n"
1551 "(logbit? 1 #b1101) @result{} #f\n"
1552 "(logbit? 2 #b1101) @result{} #t\n"
1553 "(logbit? 3 #b1101) @result{} #t\n"
1554 "(logbit? 4 #b1101) @result{} #f\n"
1556 #define FUNC_NAME s_scm_logbit_p
1558 unsigned long int iindex
;
1559 iindex
= scm_to_ulong (index
);
1561 if (SCM_I_INUMP (j
))
1563 /* bits above what's in an inum follow the sign bit */
1564 iindex
= min (iindex
, SCM_LONG_BIT
- 1);
1565 return scm_from_bool ((1L << iindex
) & SCM_I_INUM (j
));
1567 else if (SCM_BIGP (j
))
1569 int val
= mpz_tstbit (SCM_I_BIG_MPZ (j
), iindex
);
1570 scm_remember_upto_here_1 (j
);
1571 return scm_from_bool (val
);
1574 SCM_WRONG_TYPE_ARG (SCM_ARG2
, j
);
1579 SCM_DEFINE (scm_lognot
, "lognot", 1, 0, 0,
1581 "Return the integer which is the ones-complement of the integer\n"
1585 "(number->string (lognot #b10000000) 2)\n"
1586 " @result{} \"-10000001\"\n"
1587 "(number->string (lognot #b0) 2)\n"
1588 " @result{} \"-1\"\n"
1590 #define FUNC_NAME s_scm_lognot
1592 if (SCM_I_INUMP (n
)) {
1593 /* No overflow here, just need to toggle all the bits making up the inum.
1594 Enhancement: No need to strip the tag and add it back, could just xor
1595 a block of 1 bits, if that worked with the various debug versions of
1597 return SCM_I_MAKINUM (~ SCM_I_INUM (n
));
1599 } else if (SCM_BIGP (n
)) {
1600 SCM result
= scm_i_mkbig ();
1601 mpz_com (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
));
1602 scm_remember_upto_here_1 (n
);
1606 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1611 /* returns 0 if IN is not an integer. OUT must already be
1614 coerce_to_big (SCM in
, mpz_t out
)
1617 mpz_set (out
, SCM_I_BIG_MPZ (in
));
1618 else if (SCM_I_INUMP (in
))
1619 mpz_set_si (out
, SCM_I_INUM (in
));
1626 SCM_DEFINE (scm_modulo_expt
, "modulo-expt", 3, 0, 0,
1627 (SCM n
, SCM k
, SCM m
),
1628 "Return @var{n} raised to the integer exponent\n"
1629 "@var{k}, modulo @var{m}.\n"
1632 "(modulo-expt 2 3 5)\n"
1635 #define FUNC_NAME s_scm_modulo_expt
1641 /* There are two classes of error we might encounter --
1642 1) Math errors, which we'll report by calling scm_num_overflow,
1644 2) wrong-type errors, which of course we'll report by calling
1646 We don't report those errors immediately, however; instead we do
1647 some cleanup first. These variables tell us which error (if
1648 any) we should report after cleaning up.
1650 int report_overflow
= 0;
1652 int position_of_wrong_type
= 0;
1653 SCM value_of_wrong_type
= SCM_INUM0
;
1655 SCM result
= SCM_UNDEFINED
;
1661 if (scm_is_eq (m
, SCM_INUM0
))
1663 report_overflow
= 1;
1667 if (!coerce_to_big (n
, n_tmp
))
1669 value_of_wrong_type
= n
;
1670 position_of_wrong_type
= 1;
1674 if (!coerce_to_big (k
, k_tmp
))
1676 value_of_wrong_type
= k
;
1677 position_of_wrong_type
= 2;
1681 if (!coerce_to_big (m
, m_tmp
))
1683 value_of_wrong_type
= m
;
1684 position_of_wrong_type
= 3;
1688 /* if the exponent K is negative, and we simply call mpz_powm, we
1689 will get a divide-by-zero exception when an inverse 1/n mod m
1690 doesn't exist (or is not unique). Since exceptions are hard to
1691 handle, we'll attempt the inversion "by hand" -- that way, we get
1692 a simple failure code, which is easy to handle. */
1694 if (-1 == mpz_sgn (k_tmp
))
1696 if (!mpz_invert (n_tmp
, n_tmp
, m_tmp
))
1698 report_overflow
= 1;
1701 mpz_neg (k_tmp
, k_tmp
);
1704 result
= scm_i_mkbig ();
1705 mpz_powm (SCM_I_BIG_MPZ (result
),
1710 if (mpz_sgn (m_tmp
) < 0 && mpz_sgn (SCM_I_BIG_MPZ (result
)) != 0)
1711 mpz_add (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), m_tmp
);
1718 if (report_overflow
)
1719 scm_num_overflow (FUNC_NAME
);
1721 if (position_of_wrong_type
)
1722 SCM_WRONG_TYPE_ARG (position_of_wrong_type
,
1723 value_of_wrong_type
);
1725 return scm_i_normbig (result
);
1729 SCM_DEFINE (scm_integer_expt
, "integer-expt", 2, 0, 0,
1731 "Return @var{n} raised to the power @var{k}. @var{k} must be an\n"
1732 "exact integer, @var{n} can be any number.\n"
1734 "Negative @var{k} is supported, and results in @math{1/n^abs(k)}\n"
1735 "in the usual way. @math{@var{n}^0} is 1, as usual, and that\n"
1736 "includes @math{0^0} is 1.\n"
1739 "(integer-expt 2 5) @result{} 32\n"
1740 "(integer-expt -3 3) @result{} -27\n"
1741 "(integer-expt 5 -3) @result{} 1/125\n"
1742 "(integer-expt 0 0) @result{} 1\n"
1744 #define FUNC_NAME s_scm_integer_expt
1747 SCM z_i2
= SCM_BOOL_F
;
1749 SCM acc
= SCM_I_MAKINUM (1L);
1751 SCM_VALIDATE_NUMBER (SCM_ARG1
, n
);
1753 /* 0^0 == 1 according to R5RS */
1754 if (scm_is_eq (n
, SCM_INUM0
) || scm_is_eq (n
, acc
))
1755 return scm_is_false (scm_zero_p(k
)) ? n
: acc
;
1756 else if (scm_is_eq (n
, SCM_I_MAKINUM (-1L)))
1757 return scm_is_false (scm_even_p (k
)) ? n
: acc
;
1759 if (SCM_I_INUMP (k
))
1760 i2
= SCM_I_INUM (k
);
1761 else if (SCM_BIGP (k
))
1763 z_i2
= scm_i_clonebig (k
, 1);
1764 scm_remember_upto_here_1 (k
);
1768 SCM_WRONG_TYPE_ARG (2, k
);
1772 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == -1)
1774 mpz_neg (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
));
1775 n
= scm_divide (n
, SCM_UNDEFINED
);
1779 if (mpz_sgn(SCM_I_BIG_MPZ (z_i2
)) == 0)
1783 if (mpz_cmp_ui(SCM_I_BIG_MPZ (z_i2
), 1) == 0)
1785 return scm_product (acc
, n
);
1787 if (mpz_tstbit(SCM_I_BIG_MPZ (z_i2
), 0))
1788 acc
= scm_product (acc
, n
);
1789 n
= scm_product (n
, n
);
1790 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (z_i2
), SCM_I_BIG_MPZ (z_i2
), 1);
1798 n
= scm_divide (n
, SCM_UNDEFINED
);
1805 return scm_product (acc
, n
);
1807 acc
= scm_product (acc
, n
);
1808 n
= scm_product (n
, n
);
1815 SCM_DEFINE (scm_ash
, "ash", 2, 0, 0,
1817 "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
1818 "if @var{cnt} is negative. This is an ``arithmetic'' shift.\n"
1820 "This is effectively a multiplication by 2^@var{cnt}, and when\n"
1821 "@var{cnt} is negative it's a division, rounded towards negative\n"
1822 "infinity. (Note that this is not the same rounding as\n"
1823 "@code{quotient} does.)\n"
1825 "With @var{n} viewed as an infinite precision twos complement,\n"
1826 "@code{ash} means a left shift introducing zero bits, or a right\n"
1827 "shift dropping bits.\n"
1830 "(number->string (ash #b1 3) 2) @result{} \"1000\"\n"
1831 "(number->string (ash #b1010 -1) 2) @result{} \"101\"\n"
1833 ";; -23 is bits ...11101001, -6 is bits ...111010\n"
1834 "(ash -23 -2) @result{} -6\n"
1836 #define FUNC_NAME s_scm_ash
1839 bits_to_shift
= scm_to_long (cnt
);
1841 if (SCM_I_INUMP (n
))
1843 long nn
= SCM_I_INUM (n
);
1845 if (bits_to_shift
> 0)
1847 /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
1848 overflow a non-zero fixnum. For smaller shifts we check the
1849 bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
1850 all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
1851 Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
1857 if (bits_to_shift
< SCM_I_FIXNUM_BIT
-1
1859 (SCM_SRS (nn
, (SCM_I_FIXNUM_BIT
-1 - bits_to_shift
)) + 1)
1862 return SCM_I_MAKINUM (nn
<< bits_to_shift
);
1866 SCM result
= scm_i_long2big (nn
);
1867 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1874 bits_to_shift
= -bits_to_shift
;
1875 if (bits_to_shift
>= SCM_LONG_BIT
)
1876 return (nn
>= 0 ? SCM_I_MAKINUM (0) : SCM_I_MAKINUM(-1));
1878 return SCM_I_MAKINUM (SCM_SRS (nn
, bits_to_shift
));
1882 else if (SCM_BIGP (n
))
1886 if (bits_to_shift
== 0)
1889 result
= scm_i_mkbig ();
1890 if (bits_to_shift
>= 0)
1892 mpz_mul_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1898 /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
1899 we have to allocate a bignum even if the result is going to be a
1901 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (n
),
1903 return scm_i_normbig (result
);
1909 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1915 SCM_DEFINE (scm_bit_extract
, "bit-extract", 3, 0, 0,
1916 (SCM n
, SCM start
, SCM end
),
1917 "Return the integer composed of the @var{start} (inclusive)\n"
1918 "through @var{end} (exclusive) bits of @var{n}. The\n"
1919 "@var{start}th bit becomes the 0-th bit in the result.\n"
1922 "(number->string (bit-extract #b1101101010 0 4) 2)\n"
1923 " @result{} \"1010\"\n"
1924 "(number->string (bit-extract #b1101101010 4 9) 2)\n"
1925 " @result{} \"10110\"\n"
1927 #define FUNC_NAME s_scm_bit_extract
1929 unsigned long int istart
, iend
, bits
;
1930 istart
= scm_to_ulong (start
);
1931 iend
= scm_to_ulong (end
);
1932 SCM_ASSERT_RANGE (3, end
, (iend
>= istart
));
1934 /* how many bits to keep */
1935 bits
= iend
- istart
;
1937 if (SCM_I_INUMP (n
))
1939 long int in
= SCM_I_INUM (n
);
1941 /* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
1942 SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "in". */
1943 in
= SCM_SRS (in
, min (istart
, SCM_I_FIXNUM_BIT
-1));
1945 if (in
< 0 && bits
>= SCM_I_FIXNUM_BIT
)
1947 /* Since we emulate two's complement encoded numbers, this
1948 * special case requires us to produce a result that has
1949 * more bits than can be stored in a fixnum.
1951 SCM result
= scm_i_long2big (in
);
1952 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
),
1957 /* mask down to requisite bits */
1958 bits
= min (bits
, SCM_I_FIXNUM_BIT
);
1959 return SCM_I_MAKINUM (in
& ((1L << bits
) - 1));
1961 else if (SCM_BIGP (n
))
1966 result
= SCM_I_MAKINUM (mpz_tstbit (SCM_I_BIG_MPZ (n
), istart
));
1970 /* ENHANCE-ME: It'd be nice not to allocate a new bignum when
1971 bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
1972 such bits into a ulong. */
1973 result
= scm_i_mkbig ();
1974 mpz_fdiv_q_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(n
), istart
);
1975 mpz_fdiv_r_2exp (SCM_I_BIG_MPZ(result
), SCM_I_BIG_MPZ(result
), bits
);
1976 result
= scm_i_normbig (result
);
1978 scm_remember_upto_here_1 (n
);
1982 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
1987 static const char scm_logtab
[] = {
1988 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
1991 SCM_DEFINE (scm_logcount
, "logcount", 1, 0, 0,
1993 "Return the number of bits in integer @var{n}. If integer is\n"
1994 "positive, the 1-bits in its binary representation are counted.\n"
1995 "If negative, the 0-bits in its two's-complement binary\n"
1996 "representation are counted. If 0, 0 is returned.\n"
1999 "(logcount #b10101010)\n"
2006 #define FUNC_NAME s_scm_logcount
2008 if (SCM_I_INUMP (n
))
2010 unsigned long int c
= 0;
2011 long int nn
= SCM_I_INUM (n
);
2016 c
+= scm_logtab
[15 & nn
];
2019 return SCM_I_MAKINUM (c
);
2021 else if (SCM_BIGP (n
))
2023 unsigned long count
;
2024 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) >= 0)
2025 count
= mpz_popcount (SCM_I_BIG_MPZ (n
));
2027 count
= mpz_hamdist (SCM_I_BIG_MPZ (n
), z_negative_one
);
2028 scm_remember_upto_here_1 (n
);
2029 return SCM_I_MAKINUM (count
);
2032 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2037 static const char scm_ilentab
[] = {
2038 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
2042 SCM_DEFINE (scm_integer_length
, "integer-length", 1, 0, 0,
2044 "Return the number of bits necessary to represent @var{n}.\n"
2047 "(integer-length #b10101010)\n"
2049 "(integer-length 0)\n"
2051 "(integer-length #b1111)\n"
2054 #define FUNC_NAME s_scm_integer_length
2056 if (SCM_I_INUMP (n
))
2058 unsigned long int c
= 0;
2060 long int nn
= SCM_I_INUM (n
);
2066 l
= scm_ilentab
[15 & nn
];
2069 return SCM_I_MAKINUM (c
- 4 + l
);
2071 else if (SCM_BIGP (n
))
2073 /* mpz_sizeinbase looks at the absolute value of negatives, whereas we
2074 want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
2075 1 too big, so check for that and adjust. */
2076 size_t size
= mpz_sizeinbase (SCM_I_BIG_MPZ (n
), 2);
2077 if (mpz_sgn (SCM_I_BIG_MPZ (n
)) < 0
2078 && mpz_scan0 (SCM_I_BIG_MPZ (n
), /* no 0 bits above the lowest 1 */
2079 mpz_scan1 (SCM_I_BIG_MPZ (n
), 0)) == ULONG_MAX
)
2081 scm_remember_upto_here_1 (n
);
2082 return SCM_I_MAKINUM (size
);
2085 SCM_WRONG_TYPE_ARG (SCM_ARG1
, n
);
2089 /*** NUMBERS -> STRINGS ***/
2090 #define SCM_MAX_DBL_PREC 60
2091 #define SCM_MAX_DBL_RADIX 36
2093 /* this is an array starting with radix 2, and ending with radix SCM_MAX_DBL_RADIX */
2094 static int scm_dblprec
[SCM_MAX_DBL_RADIX
- 1];
2095 static double fx_per_radix
[SCM_MAX_DBL_RADIX
- 1][SCM_MAX_DBL_PREC
];
2098 void init_dblprec(int *prec
, int radix
) {
2099 /* determine floating point precision by adding successively
2100 smaller increments to 1.0 until it is considered == 1.0 */
2101 double f
= ((double)1.0)/radix
;
2102 double fsum
= 1.0 + f
;
2107 if (++(*prec
) > SCM_MAX_DBL_PREC
)
2119 void init_fx_radix(double *fx_list
, int radix
)
2121 /* initialize a per-radix list of tolerances. When added
2122 to a number < 1.0, we can determine if we should raund
2123 up and quit converting a number to a string. */
2127 for( i
= 2 ; i
< SCM_MAX_DBL_PREC
; ++i
)
2128 fx_list
[i
] = (fx_list
[i
-1] / radix
);
2131 /* use this array as a way to generate a single digit */
2132 static const char*number_chars
="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2135 idbl2str (double f
, char *a
, int radix
)
2137 int efmt
, dpt
, d
, i
, wp
;
2139 #ifdef DBL_MIN_10_EXP
2142 #endif /* DBL_MIN_10_EXP */
2147 radix
> SCM_MAX_DBL_RADIX
)
2149 /* revert to existing behavior */
2153 wp
= scm_dblprec
[radix
-2];
2154 fx
= fx_per_radix
[radix
-2];
2158 #ifdef HAVE_COPYSIGN
2159 double sgn
= copysign (1.0, f
);
2164 goto zero
; /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
2170 strcpy (a
, "-inf.0");
2172 strcpy (a
, "+inf.0");
2177 strcpy (a
, "+nan.0");
2187 #ifdef DBL_MIN_10_EXP /* Prevent unnormalized values, as from
2188 make-uniform-vector, from causing infinite loops. */
2189 /* just do the checking...if it passes, we do the conversion for our
2190 radix again below */
2197 if (exp_cpy
-- < DBL_MIN_10_EXP
)
2205 while (f_cpy
> 10.0)
2208 if (exp_cpy
++ > DBL_MAX_10_EXP
)
2229 if (f
+ fx
[wp
] >= radix
)
2236 /* adding 9999 makes this equivalent to abs(x) % 3 */
2237 dpt
= (exp
+ 9999) % 3;
2241 efmt
= (exp
< -3) || (exp
> wp
+ 2);
2263 a
[ch
++] = number_chars
[d
];
2266 if (f
+ fx
[wp
] >= 1.0)
2268 a
[ch
- 1] = number_chars
[d
+1];
2280 if ((dpt
> 4) && (exp
> 6))
2282 d
= (a
[0] == '-' ? 2 : 1);
2283 for (i
= ch
++; i
> d
; i
--)
2296 if (a
[ch
- 1] == '.')
2297 a
[ch
++] = '0'; /* trailing zero */
2306 for (i
= radix
; i
<= exp
; i
*= radix
);
2307 for (i
/= radix
; i
; i
/= radix
)
2309 a
[ch
++] = number_chars
[exp
/ i
];
2318 icmplx2str (double real
, double imag
, char *str
, int radix
)
2322 i
= idbl2str (real
, str
, radix
);
2325 /* Don't output a '+' for negative numbers or for Inf and
2326 NaN. They will provide their own sign. */
2327 if (0 <= imag
&& !isinf (imag
) && !isnan (imag
))
2329 i
+= idbl2str (imag
, &str
[i
], radix
);
2336 iflo2str (SCM flt
, char *str
, int radix
)
2339 if (SCM_REALP (flt
))
2340 i
= idbl2str (SCM_REAL_VALUE (flt
), str
, radix
);
2342 i
= icmplx2str (SCM_COMPLEX_REAL (flt
), SCM_COMPLEX_IMAG (flt
),
2347 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2348 characters in the result.
2350 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2352 scm_iint2str (scm_t_intmax num
, int rad
, char *p
)
2357 return scm_iuint2str (-num
, rad
, p
) + 1;
2360 return scm_iuint2str (num
, rad
, p
);
2363 /* convert a scm_t_intmax to a string (unterminated). returns the number of
2364 characters in the result.
2366 p is destination: worst case (base 2) is SCM_INTBUFLEN */
2368 scm_iuint2str (scm_t_uintmax num
, int rad
, char *p
)
2372 scm_t_uintmax n
= num
;
2374 for (n
/= rad
; n
> 0; n
/= rad
)
2384 p
[i
] = d
+ ((d
< 10) ? '0' : 'a' - 10);
2389 SCM_DEFINE (scm_number_to_string
, "number->string", 1, 1, 0,
2391 "Return a string holding the external representation of the\n"
2392 "number @var{n} in the given @var{radix}. If @var{n} is\n"
2393 "inexact, a radix of 10 will be used.")
2394 #define FUNC_NAME s_scm_number_to_string
2398 if (SCM_UNBNDP (radix
))
2401 base
= scm_to_signed_integer (radix
, 2, 36);
2403 if (SCM_I_INUMP (n
))
2405 char num_buf
[SCM_INTBUFLEN
];
2406 size_t length
= scm_iint2str (SCM_I_INUM (n
), base
, num_buf
);
2407 return scm_from_locale_stringn (num_buf
, length
);
2409 else if (SCM_BIGP (n
))
2411 char *str
= mpz_get_str (NULL
, base
, SCM_I_BIG_MPZ (n
));
2412 scm_remember_upto_here_1 (n
);
2413 return scm_take_locale_string (str
);
2415 else if (SCM_FRACTIONP (n
))
2417 return scm_string_append (scm_list_3 (scm_number_to_string (SCM_FRACTION_NUMERATOR (n
), radix
),
2418 scm_from_locale_string ("/"),
2419 scm_number_to_string (SCM_FRACTION_DENOMINATOR (n
), radix
)));
2421 else if (SCM_INEXACTP (n
))
2423 char num_buf
[FLOBUFLEN
];
2424 return scm_from_locale_stringn (num_buf
, iflo2str (n
, num_buf
, base
));
2427 SCM_WRONG_TYPE_ARG (1, n
);
2432 /* These print routines used to be stubbed here so that scm_repl.c
2433 wouldn't need SCM_BIGDIG conditionals (pre GMP) */
2436 scm_print_real (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2438 char num_buf
[FLOBUFLEN
];
2439 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2444 scm_i_print_double (double val
, SCM port
)
2446 char num_buf
[FLOBUFLEN
];
2447 scm_lfwrite (num_buf
, idbl2str (val
, num_buf
, 10), port
);
2451 scm_print_complex (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2454 char num_buf
[FLOBUFLEN
];
2455 scm_lfwrite (num_buf
, iflo2str (sexp
, num_buf
, 10), port
);
2460 scm_i_print_complex (double real
, double imag
, SCM port
)
2462 char num_buf
[FLOBUFLEN
];
2463 scm_lfwrite (num_buf
, icmplx2str (real
, imag
, num_buf
, 10), port
);
2467 scm_i_print_fraction (SCM sexp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2470 str
= scm_number_to_string (sexp
, SCM_UNDEFINED
);
2471 scm_lfwrite_str (str
, port
);
2472 scm_remember_upto_here_1 (str
);
2477 scm_bigprint (SCM exp
, SCM port
, scm_print_state
*pstate SCM_UNUSED
)
2479 char *str
= mpz_get_str (NULL
, 10, SCM_I_BIG_MPZ (exp
));
2480 scm_remember_upto_here_1 (exp
);
2481 scm_lfwrite (str
, (size_t) strlen (str
), port
);
2485 /*** END nums->strs ***/
2488 /*** STRINGS -> NUMBERS ***/
2490 /* The following functions implement the conversion from strings to numbers.
2491 * The implementation somehow follows the grammar for numbers as it is given
2492 * in R5RS. Thus, the functions resemble syntactic units (<ureal R>,
2493 * <uinteger R>, ...) that are used to build up numbers in the grammar. Some
2494 * points should be noted about the implementation:
2495 * * Each function keeps a local index variable 'idx' that points at the
2496 * current position within the parsed string. The global index is only
2497 * updated if the function could parse the corresponding syntactic unit
2499 * * Similarly, the functions keep track of indicators of inexactness ('#',
2500 * '.' or exponents) using local variables ('hash_seen', 'x'). Again, the
2501 * global exactness information is only updated after each part has been
2502 * successfully parsed.
2503 * * Sequences of digits are parsed into temporary variables holding fixnums.
2504 * Only if these fixnums would overflow, the result variables are updated
2505 * using the standard functions scm_add, scm_product, scm_divide etc. Then,
2506 * the temporary variables holding the fixnums are cleared, and the process
2507 * starts over again. If for example fixnums were able to store five decimal
2508 * digits, a number 1234567890 would be parsed in two parts 12345 and 67890,
2509 * and the result was computed as 12345 * 100000 + 67890. In other words,
2510 * only every five digits two bignum operations were performed.
2513 enum t_exactness
{NO_EXACTNESS
, INEXACT
, EXACT
};
2515 /* R5RS, section 7.1.1, lexical structure of numbers: <uinteger R>. */
2517 /* In non ASCII-style encodings the following macro might not work. */
2518 #define XDIGIT2UINT(d) \
2519 (uc_is_property_decimal_digit ((int) (unsigned char) d) \
2521 : uc_tolower ((int) (unsigned char) d) - 'a' + 10)
2524 mem2uinteger (SCM mem
, unsigned int *p_idx
,
2525 unsigned int radix
, enum t_exactness
*p_exactness
)
2527 unsigned int idx
= *p_idx
;
2528 unsigned int hash_seen
= 0;
2529 scm_t_bits shift
= 1;
2531 unsigned int digit_value
;
2534 size_t len
= scm_i_string_length (mem
);
2539 c
= scm_i_string_ref (mem
, idx
);
2540 if (!uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2542 digit_value
= XDIGIT2UINT (c
);
2543 if (digit_value
>= radix
)
2547 result
= SCM_I_MAKINUM (digit_value
);
2550 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2551 if (uc_is_property_ascii_hex_digit ((scm_t_uint32
) c
))
2555 digit_value
= XDIGIT2UINT (c
);
2556 if (digit_value
>= radix
)
2568 if (SCM_MOST_POSITIVE_FIXNUM
/ radix
< shift
)
2570 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2572 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2579 shift
= shift
* radix
;
2580 add
= add
* radix
+ digit_value
;
2585 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2587 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2591 *p_exactness
= INEXACT
;
2597 /* R5RS, section 7.1.1, lexical structure of numbers: <decimal 10>. Only
2598 * covers the parts of the rules that start at a potential point. The value
2599 * of the digits up to the point have been parsed by the caller and are given
2600 * in variable result. The content of *p_exactness indicates, whether a hash
2601 * has already been seen in the digits before the point.
2604 #define DIGIT2UINT(d) (uc_numeric_value(d).numerator)
2607 mem2decimal_from_point (SCM result
, SCM mem
,
2608 unsigned int *p_idx
, enum t_exactness
*p_exactness
)
2610 unsigned int idx
= *p_idx
;
2611 enum t_exactness x
= *p_exactness
;
2612 size_t len
= scm_i_string_length (mem
);
2617 if (scm_i_string_ref (mem
, idx
) == '.')
2619 scm_t_bits shift
= 1;
2621 unsigned int digit_value
;
2622 SCM big_shift
= SCM_I_MAKINUM (1);
2627 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2628 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2633 digit_value
= DIGIT2UINT (c
);
2644 if (SCM_MOST_POSITIVE_FIXNUM
/ 10 < shift
)
2646 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2647 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2649 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2657 add
= add
* 10 + digit_value
;
2663 big_shift
= scm_product (big_shift
, SCM_I_MAKINUM (shift
));
2664 result
= scm_product (result
, SCM_I_MAKINUM (shift
));
2665 result
= scm_sum (result
, SCM_I_MAKINUM (add
));
2668 result
= scm_divide (result
, big_shift
);
2670 /* We've seen a decimal point, thus the value is implicitly inexact. */
2682 /* R5RS, section 7.1.1, lexical structure of numbers: <suffix> */
2684 switch (scm_i_string_ref (mem
, idx
))
2696 c
= scm_i_string_ref (mem
, idx
);
2704 c
= scm_i_string_ref (mem
, idx
);
2713 c
= scm_i_string_ref (mem
, idx
);
2718 if (!uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2722 exponent
= DIGIT2UINT (c
);
2725 scm_t_wchar c
= scm_i_string_ref (mem
, idx
);
2726 if (uc_is_property_decimal_digit ((scm_t_uint32
) c
))
2729 if (exponent
<= SCM_MAXEXP
)
2730 exponent
= exponent
* 10 + DIGIT2UINT (c
);
2736 if (exponent
> SCM_MAXEXP
)
2738 size_t exp_len
= idx
- start
;
2739 SCM exp_string
= scm_i_substring_copy (mem
, start
, start
+ exp_len
);
2740 SCM exp_num
= scm_string_to_number (exp_string
, SCM_UNDEFINED
);
2741 scm_out_of_range ("string->number", exp_num
);
2744 e
= scm_integer_expt (SCM_I_MAKINUM (10), SCM_I_MAKINUM (exponent
));
2746 result
= scm_product (result
, e
);
2748 result
= scm_divide2real (result
, e
);
2750 /* We've seen an exponent, thus the value is implicitly inexact. */
2768 /* R5RS, section 7.1.1, lexical structure of numbers: <ureal R> */
2771 mem2ureal (SCM mem
, unsigned int *p_idx
,
2772 unsigned int radix
, enum t_exactness
*p_exactness
)
2774 unsigned int idx
= *p_idx
;
2776 size_t len
= scm_i_string_length (mem
);
2778 /* Start off believing that the number will be exact. This changes
2779 to INEXACT if we see a decimal point or a hash. */
2780 enum t_exactness x
= EXACT
;
2785 if (idx
+5 <= len
&& !scm_i_string_strcmp (mem
, idx
, "inf.0"))
2791 if (idx
+4 < len
&& !scm_i_string_strcmp (mem
, idx
, "nan."))
2793 /* Cobble up the fractional part. We might want to set the
2794 NaN's mantissa from it. */
2796 mem2uinteger (mem
, &idx
, 10, &x
);
2801 if (scm_i_string_ref (mem
, idx
) == '.')
2805 else if (idx
+ 1 == len
)
2807 else if (!uc_is_property_decimal_digit ((scm_t_uint32
) scm_i_string_ref (mem
, idx
+1)))
2810 result
= mem2decimal_from_point (SCM_I_MAKINUM (0), mem
,
2817 uinteger
= mem2uinteger (mem
, &idx
, radix
, &x
);
2818 if (scm_is_false (uinteger
))
2823 else if (scm_i_string_ref (mem
, idx
) == '/')
2831 divisor
= mem2uinteger (mem
, &idx
, radix
, &x
);
2832 if (scm_is_false (divisor
))
2835 /* both are int/big here, I assume */
2836 result
= scm_i_make_ratio (uinteger
, divisor
);
2838 else if (radix
== 10)
2840 result
= mem2decimal_from_point (uinteger
, mem
, &idx
, &x
);
2841 if (scm_is_false (result
))
2850 /* Update *p_exactness if the number just read was inexact. This is
2851 important for complex numbers, so that a complex number is
2852 treated as inexact overall if either its real or imaginary part
2858 /* When returning an inexact zero, make sure it is represented as a
2859 floating point value so that we can change its sign.
2861 if (scm_is_eq (result
, SCM_I_MAKINUM(0)) && *p_exactness
== INEXACT
)
2862 result
= scm_from_double (0.0);
2868 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
2871 mem2complex (SCM mem
, unsigned int idx
,
2872 unsigned int radix
, enum t_exactness
*p_exactness
)
2877 size_t len
= scm_i_string_length (mem
);
2882 c
= scm_i_string_ref (mem
, idx
);
2897 ureal
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2898 if (scm_is_false (ureal
))
2900 /* input must be either +i or -i */
2905 if (scm_i_string_ref (mem
, idx
) == 'i'
2906 || scm_i_string_ref (mem
, idx
) == 'I')
2912 return scm_make_rectangular (SCM_I_MAKINUM (0), SCM_I_MAKINUM (sign
));
2919 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2920 ureal
= scm_difference (ureal
, SCM_UNDEFINED
);
2925 c
= scm_i_string_ref (mem
, idx
);
2929 /* either +<ureal>i or -<ureal>i */
2936 return scm_make_rectangular (SCM_I_MAKINUM (0), ureal
);
2939 /* polar input: <real>@<real>. */
2950 c
= scm_i_string_ref (mem
, idx
);
2968 angle
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2969 if (scm_is_false (angle
))
2974 if (sign
== -1 && scm_is_false (scm_nan_p (ureal
)))
2975 angle
= scm_difference (angle
, SCM_UNDEFINED
);
2977 result
= scm_make_polar (ureal
, angle
);
2982 /* expecting input matching <real>[+-]<ureal>?i */
2989 int sign
= (c
== '+') ? 1 : -1;
2990 SCM imag
= mem2ureal (mem
, &idx
, radix
, p_exactness
);
2992 if (scm_is_false (imag
))
2993 imag
= SCM_I_MAKINUM (sign
);
2994 else if (sign
== -1 && scm_is_false (scm_nan_p (imag
)))
2995 imag
= scm_difference (imag
, SCM_UNDEFINED
);
2999 if (scm_i_string_ref (mem
, idx
) != 'i'
3000 && scm_i_string_ref (mem
, idx
) != 'I')
3007 return scm_make_rectangular (ureal
, imag
);
3016 /* R5RS, section 7.1.1, lexical structure of numbers: <number> */
3018 enum t_radix
{NO_RADIX
=0, DUAL
=2, OCT
=8, DEC
=10, HEX
=16};
3021 scm_i_string_to_number (SCM mem
, unsigned int default_radix
)
3023 unsigned int idx
= 0;
3024 unsigned int radix
= NO_RADIX
;
3025 enum t_exactness forced_x
= NO_EXACTNESS
;
3026 enum t_exactness implicit_x
= EXACT
;
3028 size_t len
= scm_i_string_length (mem
);
3030 /* R5RS, section 7.1.1, lexical structure of numbers: <prefix R> */
3031 while (idx
+ 2 < len
&& scm_i_string_ref (mem
, idx
) == '#')
3033 switch (scm_i_string_ref (mem
, idx
+ 1))
3036 if (radix
!= NO_RADIX
)
3041 if (radix
!= NO_RADIX
)
3046 if (forced_x
!= NO_EXACTNESS
)
3051 if (forced_x
!= NO_EXACTNESS
)
3056 if (radix
!= NO_RADIX
)
3061 if (radix
!= NO_RADIX
)
3071 /* R5RS, section 7.1.1, lexical structure of numbers: <complex R> */
3072 if (radix
== NO_RADIX
)
3073 result
= mem2complex (mem
, idx
, default_radix
, &implicit_x
);
3075 result
= mem2complex (mem
, idx
, (unsigned int) radix
, &implicit_x
);
3077 if (scm_is_false (result
))
3083 if (SCM_INEXACTP (result
))
3084 return scm_inexact_to_exact (result
);
3088 if (SCM_INEXACTP (result
))
3091 return scm_exact_to_inexact (result
);
3094 if (implicit_x
== INEXACT
)
3096 if (SCM_INEXACTP (result
))
3099 return scm_exact_to_inexact (result
);
3107 scm_c_locale_stringn_to_number (const char* mem
, size_t len
,
3108 unsigned int default_radix
)
3110 SCM str
= scm_from_locale_stringn (mem
, len
);
3112 return scm_i_string_to_number (str
, default_radix
);
3116 SCM_DEFINE (scm_string_to_number
, "string->number", 1, 1, 0,
3117 (SCM string
, SCM radix
),
3118 "Return a number of the maximally precise representation\n"
3119 "expressed by the given @var{string}. @var{radix} must be an\n"
3120 "exact integer, either 2, 8, 10, or 16. If supplied, @var{radix}\n"
3121 "is a default radix that may be overridden by an explicit radix\n"
3122 "prefix in @var{string} (e.g. \"#o177\"). If @var{radix} is not\n"
3123 "supplied, then the default radix is 10. If string is not a\n"
3124 "syntactically valid notation for a number, then\n"
3125 "@code{string->number} returns @code{#f}.")
3126 #define FUNC_NAME s_scm_string_to_number
3130 SCM_VALIDATE_STRING (1, string
);
3132 if (SCM_UNBNDP (radix
))
3135 base
= scm_to_unsigned_integer (radix
, 2, INT_MAX
);
3137 answer
= scm_i_string_to_number (string
, base
);
3138 scm_remember_upto_here_1 (string
);
3144 /*** END strs->nums ***/
3148 scm_bigequal (SCM x
, SCM y
)
3150 int result
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3151 scm_remember_upto_here_2 (x
, y
);
3152 return scm_from_bool (0 == result
);
3156 scm_real_equalp (SCM x
, SCM y
)
3158 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3162 scm_complex_equalp (SCM x
, SCM y
)
3164 return scm_from_bool (SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
)
3165 && SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
));
3169 scm_i_fraction_equalp (SCM x
, SCM y
)
3171 if (scm_is_false (scm_equal_p (SCM_FRACTION_NUMERATOR (x
),
3172 SCM_FRACTION_NUMERATOR (y
)))
3173 || scm_is_false (scm_equal_p (SCM_FRACTION_DENOMINATOR (x
),
3174 SCM_FRACTION_DENOMINATOR (y
))))
3181 SCM_DEFINE (scm_number_p
, "number?", 1, 0, 0,
3183 "Return @code{#t} if @var{x} is a number, @code{#f}\n"
3185 #define FUNC_NAME s_scm_number_p
3187 return scm_from_bool (SCM_NUMBERP (x
));
3191 SCM_DEFINE (scm_complex_p
, "complex?", 1, 0, 0,
3193 "Return @code{#t} if @var{x} is a complex number, @code{#f}\n"
3194 "otherwise. Note that the sets of real, rational and integer\n"
3195 "values form subsets of the set of complex numbers, i. e. the\n"
3196 "predicate will also be fulfilled if @var{x} is a real,\n"
3197 "rational or integer number.")
3198 #define FUNC_NAME s_scm_complex_p
3200 /* all numbers are complex. */
3201 return scm_number_p (x
);
3205 SCM_DEFINE (scm_real_p
, "real?", 1, 0, 0,
3207 "Return @code{#t} if @var{x} is a real number, @code{#f}\n"
3208 "otherwise. Note that the set of integer values forms a subset of\n"
3209 "the set of real numbers, i. e. the predicate will also be\n"
3210 "fulfilled if @var{x} is an integer number.")
3211 #define FUNC_NAME s_scm_real_p
3213 /* we can't represent irrational numbers. */
3214 return scm_rational_p (x
);
3218 SCM_DEFINE (scm_rational_p
, "rational?", 1, 0, 0,
3220 "Return @code{#t} if @var{x} is a rational number, @code{#f}\n"
3221 "otherwise. Note that the set of integer values forms a subset of\n"
3222 "the set of rational numbers, i. e. the predicate will also be\n"
3223 "fulfilled if @var{x} is an integer number.")
3224 #define FUNC_NAME s_scm_rational_p
3226 if (SCM_I_INUMP (x
))
3228 else if (SCM_IMP (x
))
3230 else if (SCM_BIGP (x
))
3232 else if (SCM_FRACTIONP (x
))
3234 else if (SCM_REALP (x
))
3235 /* due to their limited precision, all floating point numbers are
3236 rational as well. */
3243 SCM_DEFINE (scm_integer_p
, "integer?", 1, 0, 0,
3245 "Return @code{#t} if @var{x} is an integer number, @code{#f}\n"
3247 #define FUNC_NAME s_scm_integer_p
3250 if (SCM_I_INUMP (x
))
3256 if (!SCM_INEXACTP (x
))
3258 if (SCM_COMPLEXP (x
))
3260 r
= SCM_REAL_VALUE (x
);
3261 /* +/-inf passes r==floor(r), making those #t */
3269 SCM_DEFINE (scm_inexact_p
, "inexact?", 1, 0, 0,
3271 "Return @code{#t} if @var{x} is an inexact number, @code{#f}\n"
3273 #define FUNC_NAME s_scm_inexact_p
3275 if (SCM_INEXACTP (x
))
3277 if (SCM_NUMBERP (x
))
3279 SCM_WRONG_TYPE_ARG (1, x
);
3284 SCM
scm_i_num_eq_p (SCM
, SCM
, SCM
);
3285 SCM_PRIMITIVE_GENERIC (scm_i_num_eq_p
, "=", 0, 2, 1,
3286 (SCM x
, SCM y
, SCM rest
),
3287 "Return @code{#t} if all parameters are numerically equal.")
3288 #define FUNC_NAME s_scm_i_num_eq_p
3290 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3292 while (!scm_is_null (rest
))
3294 if (scm_is_false (scm_num_eq_p (x
, y
)))
3298 rest
= scm_cdr (rest
);
3300 return scm_num_eq_p (x
, y
);
3304 scm_num_eq_p (SCM x
, SCM y
)
3307 if (SCM_I_INUMP (x
))
3309 long xx
= SCM_I_INUM (x
);
3310 if (SCM_I_INUMP (y
))
3312 long yy
= SCM_I_INUM (y
);
3313 return scm_from_bool (xx
== yy
);
3315 else if (SCM_BIGP (y
))
3317 else if (SCM_REALP (y
))
3319 /* On a 32-bit system an inum fits a double, we can cast the inum
3320 to a double and compare.
3322 But on a 64-bit system an inum is bigger than a double and
3323 casting it to a double (call that dxx) will round. dxx is at
3324 worst 1 bigger or smaller than xx, so if dxx==yy we know yy is
3325 an integer and fits a long. So we cast yy to a long and
3326 compare with plain xx.
3328 An alternative (for any size system actually) would be to check
3329 yy is an integer (with floor) and is in range of an inum
3330 (compare against appropriate powers of 2) then test
3331 xx==(long)yy. It's just a matter of which casts/comparisons
3332 might be fastest or easiest for the cpu. */
3334 double yy
= SCM_REAL_VALUE (y
);
3335 return scm_from_bool ((double) xx
== yy
3336 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3337 || xx
== (long) yy
));
3339 else if (SCM_COMPLEXP (y
))
3340 return scm_from_bool (((double) xx
== SCM_COMPLEX_REAL (y
))
3341 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3342 else if (SCM_FRACTIONP (y
))
3345 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3347 else if (SCM_BIGP (x
))
3349 if (SCM_I_INUMP (y
))
3351 else if (SCM_BIGP (y
))
3353 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3354 scm_remember_upto_here_2 (x
, y
);
3355 return scm_from_bool (0 == cmp
);
3357 else if (SCM_REALP (y
))
3360 if (isnan (SCM_REAL_VALUE (y
)))
3362 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3363 scm_remember_upto_here_1 (x
);
3364 return scm_from_bool (0 == cmp
);
3366 else if (SCM_COMPLEXP (y
))
3369 if (0.0 != SCM_COMPLEX_IMAG (y
))
3371 if (isnan (SCM_COMPLEX_REAL (y
)))
3373 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_COMPLEX_REAL (y
));
3374 scm_remember_upto_here_1 (x
);
3375 return scm_from_bool (0 == cmp
);
3377 else if (SCM_FRACTIONP (y
))
3380 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3382 else if (SCM_REALP (x
))
3384 double xx
= SCM_REAL_VALUE (x
);
3385 if (SCM_I_INUMP (y
))
3387 /* see comments with inum/real above */
3388 long yy
= SCM_I_INUM (y
);
3389 return scm_from_bool (xx
== (double) yy
3390 && (DBL_MANT_DIG
>= SCM_I_FIXNUM_BIT
-1
3391 || (long) xx
== yy
));
3393 else if (SCM_BIGP (y
))
3396 if (isnan (SCM_REAL_VALUE (x
)))
3398 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3399 scm_remember_upto_here_1 (y
);
3400 return scm_from_bool (0 == cmp
);
3402 else if (SCM_REALP (y
))
3403 return scm_from_bool (SCM_REAL_VALUE (x
) == SCM_REAL_VALUE (y
));
3404 else if (SCM_COMPLEXP (y
))
3405 return scm_from_bool ((SCM_REAL_VALUE (x
) == SCM_COMPLEX_REAL (y
))
3406 && (0.0 == SCM_COMPLEX_IMAG (y
)));
3407 else if (SCM_FRACTIONP (y
))
3409 double xx
= SCM_REAL_VALUE (x
);
3413 return scm_from_bool (xx
< 0.0);
3414 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3418 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3420 else if (SCM_COMPLEXP (x
))
3422 if (SCM_I_INUMP (y
))
3423 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == (double) SCM_I_INUM (y
))
3424 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3425 else if (SCM_BIGP (y
))
3428 if (0.0 != SCM_COMPLEX_IMAG (x
))
3430 if (isnan (SCM_COMPLEX_REAL (x
)))
3432 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_COMPLEX_REAL (x
));
3433 scm_remember_upto_here_1 (y
);
3434 return scm_from_bool (0 == cmp
);
3436 else if (SCM_REALP (y
))
3437 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_REAL_VALUE (y
))
3438 && (SCM_COMPLEX_IMAG (x
) == 0.0));
3439 else if (SCM_COMPLEXP (y
))
3440 return scm_from_bool ((SCM_COMPLEX_REAL (x
) == SCM_COMPLEX_REAL (y
))
3441 && (SCM_COMPLEX_IMAG (x
) == SCM_COMPLEX_IMAG (y
)));
3442 else if (SCM_FRACTIONP (y
))
3445 if (SCM_COMPLEX_IMAG (x
) != 0.0)
3447 xx
= SCM_COMPLEX_REAL (x
);
3451 return scm_from_bool (xx
< 0.0);
3452 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3456 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3458 else if (SCM_FRACTIONP (x
))
3460 if (SCM_I_INUMP (y
))
3462 else if (SCM_BIGP (y
))
3464 else if (SCM_REALP (y
))
3466 double yy
= SCM_REAL_VALUE (y
);
3470 return scm_from_bool (0.0 < yy
);
3471 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3474 else if (SCM_COMPLEXP (y
))
3477 if (SCM_COMPLEX_IMAG (y
) != 0.0)
3479 yy
= SCM_COMPLEX_REAL (y
);
3483 return scm_from_bool (0.0 < yy
);
3484 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3487 else if (SCM_FRACTIONP (y
))
3488 return scm_i_fraction_equalp (x
, y
);
3490 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARGn
, s_scm_i_num_eq_p
);
3493 SCM_WTA_DISPATCH_2 (g_scm_i_num_eq_p
, x
, y
, SCM_ARG1
, s_scm_i_num_eq_p
);
3497 /* OPTIMIZE-ME: For int/frac and frac/frac compares, the multiplications
3498 done are good for inums, but for bignums an answer can almost always be
3499 had by just examining a few high bits of the operands, as done by GMP in
3500 mpq_cmp. flonum/frac compares likewise, but with the slight complication
3501 of the float exponent to take into account. */
3503 SCM_INTERNAL SCM
scm_i_num_less_p (SCM
, SCM
, SCM
);
3504 SCM_PRIMITIVE_GENERIC (scm_i_num_less_p
, "<", 0, 2, 1,
3505 (SCM x
, SCM y
, SCM rest
),
3506 "Return @code{#t} if the list of parameters is monotonically\n"
3508 #define FUNC_NAME s_scm_i_num_less_p
3510 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3512 while (!scm_is_null (rest
))
3514 if (scm_is_false (scm_less_p (x
, y
)))
3518 rest
= scm_cdr (rest
);
3520 return scm_less_p (x
, y
);
3524 scm_less_p (SCM x
, SCM y
)
3527 if (SCM_I_INUMP (x
))
3529 long xx
= SCM_I_INUM (x
);
3530 if (SCM_I_INUMP (y
))
3532 long yy
= SCM_I_INUM (y
);
3533 return scm_from_bool (xx
< yy
);
3535 else if (SCM_BIGP (y
))
3537 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3538 scm_remember_upto_here_1 (y
);
3539 return scm_from_bool (sgn
> 0);
3541 else if (SCM_REALP (y
))
3542 return scm_from_bool ((double) xx
< SCM_REAL_VALUE (y
));
3543 else if (SCM_FRACTIONP (y
))
3545 /* "x < a/b" becomes "x*b < a" */
3547 x
= scm_product (x
, SCM_FRACTION_DENOMINATOR (y
));
3548 y
= SCM_FRACTION_NUMERATOR (y
);
3552 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3554 else if (SCM_BIGP (x
))
3556 if (SCM_I_INUMP (y
))
3558 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3559 scm_remember_upto_here_1 (x
);
3560 return scm_from_bool (sgn
< 0);
3562 else if (SCM_BIGP (y
))
3564 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3565 scm_remember_upto_here_2 (x
, y
);
3566 return scm_from_bool (cmp
< 0);
3568 else if (SCM_REALP (y
))
3571 if (isnan (SCM_REAL_VALUE (y
)))
3573 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (x
), SCM_REAL_VALUE (y
));
3574 scm_remember_upto_here_1 (x
);
3575 return scm_from_bool (cmp
< 0);
3577 else if (SCM_FRACTIONP (y
))
3580 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3582 else if (SCM_REALP (x
))
3584 if (SCM_I_INUMP (y
))
3585 return scm_from_bool (SCM_REAL_VALUE (x
) < (double) SCM_I_INUM (y
));
3586 else if (SCM_BIGP (y
))
3589 if (isnan (SCM_REAL_VALUE (x
)))
3591 cmp
= xmpz_cmp_d (SCM_I_BIG_MPZ (y
), SCM_REAL_VALUE (x
));
3592 scm_remember_upto_here_1 (y
);
3593 return scm_from_bool (cmp
> 0);
3595 else if (SCM_REALP (y
))
3596 return scm_from_bool (SCM_REAL_VALUE (x
) < SCM_REAL_VALUE (y
));
3597 else if (SCM_FRACTIONP (y
))
3599 double xx
= SCM_REAL_VALUE (x
);
3603 return scm_from_bool (xx
< 0.0);
3604 x
= scm_inexact_to_exact (x
); /* with x as frac or int */
3608 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3610 else if (SCM_FRACTIONP (x
))
3612 if (SCM_I_INUMP (y
) || SCM_BIGP (y
))
3614 /* "a/b < y" becomes "a < y*b" */
3615 y
= scm_product (y
, SCM_FRACTION_DENOMINATOR (x
));
3616 x
= SCM_FRACTION_NUMERATOR (x
);
3619 else if (SCM_REALP (y
))
3621 double yy
= SCM_REAL_VALUE (y
);
3625 return scm_from_bool (0.0 < yy
);
3626 y
= scm_inexact_to_exact (y
); /* with y as frac or int */
3629 else if (SCM_FRACTIONP (y
))
3631 /* "a/b < c/d" becomes "a*d < c*b" */
3632 SCM new_x
= scm_product (SCM_FRACTION_NUMERATOR (x
),
3633 SCM_FRACTION_DENOMINATOR (y
));
3634 SCM new_y
= scm_product (SCM_FRACTION_NUMERATOR (y
),
3635 SCM_FRACTION_DENOMINATOR (x
));
3641 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARGn
, s_scm_i_num_less_p
);
3644 SCM_WTA_DISPATCH_2 (g_scm_i_num_less_p
, x
, y
, SCM_ARG1
, s_scm_i_num_less_p
);
3648 SCM
scm_i_num_gr_p (SCM
, SCM
, SCM
);
3649 SCM_PRIMITIVE_GENERIC (scm_i_num_gr_p
, ">", 0, 2, 1,
3650 (SCM x
, SCM y
, SCM rest
),
3651 "Return @code{#t} if the list of parameters is monotonically\n"
3653 #define FUNC_NAME s_scm_i_num_gr_p
3655 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3657 while (!scm_is_null (rest
))
3659 if (scm_is_false (scm_gr_p (x
, y
)))
3663 rest
= scm_cdr (rest
);
3665 return scm_gr_p (x
, y
);
3668 #define FUNC_NAME s_scm_i_num_gr_p
3670 scm_gr_p (SCM x
, SCM y
)
3672 if (!SCM_NUMBERP (x
))
3673 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3674 else if (!SCM_NUMBERP (y
))
3675 SCM_WTA_DISPATCH_2 (g_scm_i_num_gr_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3677 return scm_less_p (y
, x
);
3682 SCM
scm_i_num_leq_p (SCM
, SCM
, SCM
);
3683 SCM_PRIMITIVE_GENERIC (scm_i_num_leq_p
, "<=", 0, 2, 1,
3684 (SCM x
, SCM y
, SCM rest
),
3685 "Return @code{#t} if the list of parameters is monotonically\n"
3687 #define FUNC_NAME s_scm_i_num_leq_p
3689 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3691 while (!scm_is_null (rest
))
3693 if (scm_is_false (scm_leq_p (x
, y
)))
3697 rest
= scm_cdr (rest
);
3699 return scm_leq_p (x
, y
);
3702 #define FUNC_NAME s_scm_i_num_leq_p
3704 scm_leq_p (SCM x
, SCM y
)
3706 if (!SCM_NUMBERP (x
))
3707 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3708 else if (!SCM_NUMBERP (y
))
3709 SCM_WTA_DISPATCH_2 (g_scm_i_num_leq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3710 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3713 return scm_not (scm_less_p (y
, x
));
3718 SCM
scm_i_num_geq_p (SCM
, SCM
, SCM
);
3719 SCM_PRIMITIVE_GENERIC (scm_i_num_geq_p
, ">=", 0, 2, 1,
3720 (SCM x
, SCM y
, SCM rest
),
3721 "Return @code{#t} if the list of parameters is monotonically\n"
3723 #define FUNC_NAME s_scm_i_num_geq_p
3725 if (SCM_UNBNDP (x
) || SCM_UNBNDP (y
))
3727 while (!scm_is_null (rest
))
3729 if (scm_is_false (scm_geq_p (x
, y
)))
3733 rest
= scm_cdr (rest
);
3735 return scm_geq_p (x
, y
);
3738 #define FUNC_NAME s_scm_i_num_geq_p
3740 scm_geq_p (SCM x
, SCM y
)
3742 if (!SCM_NUMBERP (x
))
3743 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG1
, FUNC_NAME
);
3744 else if (!SCM_NUMBERP (y
))
3745 SCM_WTA_DISPATCH_2 (g_scm_i_num_geq_p
, x
, y
, SCM_ARG2
, FUNC_NAME
);
3746 else if (scm_is_true (scm_nan_p (x
)) || scm_is_true (scm_nan_p (y
)))
3749 return scm_not (scm_less_p (x
, y
));
3754 SCM_GPROC (s_zero_p
, "zero?", 1, 0, 0, scm_zero_p
, g_zero_p
);
3755 /* "Return @code{#t} if @var{z} is an exact or inexact number equal to\n"
3761 if (SCM_I_INUMP (z
))
3762 return scm_from_bool (scm_is_eq (z
, SCM_INUM0
));
3763 else if (SCM_BIGP (z
))
3765 else if (SCM_REALP (z
))
3766 return scm_from_bool (SCM_REAL_VALUE (z
) == 0.0);
3767 else if (SCM_COMPLEXP (z
))
3768 return scm_from_bool (SCM_COMPLEX_REAL (z
) == 0.0
3769 && SCM_COMPLEX_IMAG (z
) == 0.0);
3770 else if (SCM_FRACTIONP (z
))
3773 SCM_WTA_DISPATCH_1 (g_zero_p
, z
, SCM_ARG1
, s_zero_p
);
3777 SCM_GPROC (s_positive_p
, "positive?", 1, 0, 0, scm_positive_p
, g_positive_p
);
3778 /* "Return @code{#t} if @var{x} is an exact or inexact number greater than\n"
3782 scm_positive_p (SCM x
)
3784 if (SCM_I_INUMP (x
))
3785 return scm_from_bool (SCM_I_INUM (x
) > 0);
3786 else if (SCM_BIGP (x
))
3788 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3789 scm_remember_upto_here_1 (x
);
3790 return scm_from_bool (sgn
> 0);
3792 else if (SCM_REALP (x
))
3793 return scm_from_bool(SCM_REAL_VALUE (x
) > 0.0);
3794 else if (SCM_FRACTIONP (x
))
3795 return scm_positive_p (SCM_FRACTION_NUMERATOR (x
));
3797 SCM_WTA_DISPATCH_1 (g_positive_p
, x
, SCM_ARG1
, s_positive_p
);
3801 SCM_GPROC (s_negative_p
, "negative?", 1, 0, 0, scm_negative_p
, g_negative_p
);
3802 /* "Return @code{#t} if @var{x} is an exact or inexact number less than\n"
3806 scm_negative_p (SCM x
)
3808 if (SCM_I_INUMP (x
))
3809 return scm_from_bool (SCM_I_INUM (x
) < 0);
3810 else if (SCM_BIGP (x
))
3812 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3813 scm_remember_upto_here_1 (x
);
3814 return scm_from_bool (sgn
< 0);
3816 else if (SCM_REALP (x
))
3817 return scm_from_bool(SCM_REAL_VALUE (x
) < 0.0);
3818 else if (SCM_FRACTIONP (x
))
3819 return scm_negative_p (SCM_FRACTION_NUMERATOR (x
));
3821 SCM_WTA_DISPATCH_1 (g_negative_p
, x
, SCM_ARG1
, s_negative_p
);
3825 /* scm_min and scm_max return an inexact when either argument is inexact, as
3826 required by r5rs. On that basis, for exact/inexact combinations the
3827 exact is converted to inexact to compare and possibly return. This is
3828 unlike scm_less_p above which takes some trouble to preserve all bits in
3829 its test, such trouble is not required for min and max. */
3831 SCM_PRIMITIVE_GENERIC (scm_i_max
, "max", 0, 2, 1,
3832 (SCM x
, SCM y
, SCM rest
),
3833 "Return the maximum of all parameter values.")
3834 #define FUNC_NAME s_scm_i_max
3836 while (!scm_is_null (rest
))
3837 { x
= scm_max (x
, y
);
3839 rest
= scm_cdr (rest
);
3841 return scm_max (x
, y
);
3845 #define s_max s_scm_i_max
3846 #define g_max g_scm_i_max
3849 scm_max (SCM x
, SCM y
)
3854 SCM_WTA_DISPATCH_0 (g_max
, s_max
);
3855 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
3858 SCM_WTA_DISPATCH_1 (g_max
, x
, SCM_ARG1
, s_max
);
3861 if (SCM_I_INUMP (x
))
3863 long xx
= SCM_I_INUM (x
);
3864 if (SCM_I_INUMP (y
))
3866 long yy
= SCM_I_INUM (y
);
3867 return (xx
< yy
) ? y
: x
;
3869 else if (SCM_BIGP (y
))
3871 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
3872 scm_remember_upto_here_1 (y
);
3873 return (sgn
< 0) ? x
: y
;
3875 else if (SCM_REALP (y
))
3878 /* if y==NaN then ">" is false and we return NaN */
3879 return (z
> SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
3881 else if (SCM_FRACTIONP (y
))
3884 return (scm_is_false (scm_less_p (x
, y
)) ? x
: y
);
3887 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3889 else if (SCM_BIGP (x
))
3891 if (SCM_I_INUMP (y
))
3893 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
3894 scm_remember_upto_here_1 (x
);
3895 return (sgn
< 0) ? y
: x
;
3897 else if (SCM_BIGP (y
))
3899 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
3900 scm_remember_upto_here_2 (x
, y
);
3901 return (cmp
> 0) ? x
: y
;
3903 else if (SCM_REALP (y
))
3905 /* if y==NaN then xx>yy is false, so we return the NaN y */
3908 xx
= scm_i_big2dbl (x
);
3909 yy
= SCM_REAL_VALUE (y
);
3910 return (xx
> yy
? scm_from_double (xx
) : y
);
3912 else if (SCM_FRACTIONP (y
))
3917 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3919 else if (SCM_REALP (x
))
3921 if (SCM_I_INUMP (y
))
3923 double z
= SCM_I_INUM (y
);
3924 /* if x==NaN then "<" is false and we return NaN */
3925 return (SCM_REAL_VALUE (x
) < z
) ? scm_from_double (z
) : x
;
3927 else if (SCM_BIGP (y
))
3932 else if (SCM_REALP (y
))
3934 /* if x==NaN then our explicit check means we return NaN
3935 if y==NaN then ">" is false and we return NaN
3936 calling isnan is unavoidable, since it's the only way to know
3937 which of x or y causes any compares to be false */
3938 double xx
= SCM_REAL_VALUE (x
);
3939 return (isnan (xx
) || xx
> SCM_REAL_VALUE (y
)) ? x
: y
;
3941 else if (SCM_FRACTIONP (y
))
3943 double yy
= scm_i_fraction2double (y
);
3944 double xx
= SCM_REAL_VALUE (x
);
3945 return (xx
< yy
) ? scm_from_double (yy
) : x
;
3948 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3950 else if (SCM_FRACTIONP (x
))
3952 if (SCM_I_INUMP (y
))
3956 else if (SCM_BIGP (y
))
3960 else if (SCM_REALP (y
))
3962 double xx
= scm_i_fraction2double (x
);
3963 return (xx
< SCM_REAL_VALUE (y
)) ? y
: scm_from_double (xx
);
3965 else if (SCM_FRACTIONP (y
))
3970 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARGn
, s_max
);
3973 SCM_WTA_DISPATCH_2 (g_max
, x
, y
, SCM_ARG1
, s_max
);
3977 SCM_PRIMITIVE_GENERIC (scm_i_min
, "min", 0, 2, 1,
3978 (SCM x
, SCM y
, SCM rest
),
3979 "Return the minimum of all parameter values.")
3980 #define FUNC_NAME s_scm_i_min
3982 while (!scm_is_null (rest
))
3983 { x
= scm_min (x
, y
);
3985 rest
= scm_cdr (rest
);
3987 return scm_min (x
, y
);
3991 #define s_min s_scm_i_min
3992 #define g_min g_scm_i_min
3995 scm_min (SCM x
, SCM y
)
4000 SCM_WTA_DISPATCH_0 (g_min
, s_min
);
4001 else if (SCM_I_INUMP(x
) || SCM_BIGP(x
) || SCM_REALP(x
) || SCM_FRACTIONP(x
))
4004 SCM_WTA_DISPATCH_1 (g_min
, x
, SCM_ARG1
, s_min
);
4007 if (SCM_I_INUMP (x
))
4009 long xx
= SCM_I_INUM (x
);
4010 if (SCM_I_INUMP (y
))
4012 long yy
= SCM_I_INUM (y
);
4013 return (xx
< yy
) ? x
: y
;
4015 else if (SCM_BIGP (y
))
4017 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4018 scm_remember_upto_here_1 (y
);
4019 return (sgn
< 0) ? y
: x
;
4021 else if (SCM_REALP (y
))
4024 /* if y==NaN then "<" is false and we return NaN */
4025 return (z
< SCM_REAL_VALUE (y
)) ? scm_from_double (z
) : y
;
4027 else if (SCM_FRACTIONP (y
))
4030 return (scm_is_false (scm_less_p (x
, y
)) ? y
: x
);
4033 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4035 else if (SCM_BIGP (x
))
4037 if (SCM_I_INUMP (y
))
4039 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4040 scm_remember_upto_here_1 (x
);
4041 return (sgn
< 0) ? x
: y
;
4043 else if (SCM_BIGP (y
))
4045 int cmp
= mpz_cmp (SCM_I_BIG_MPZ (x
), SCM_I_BIG_MPZ (y
));
4046 scm_remember_upto_here_2 (x
, y
);
4047 return (cmp
> 0) ? y
: x
;
4049 else if (SCM_REALP (y
))
4051 /* if y==NaN then xx<yy is false, so we return the NaN y */
4054 xx
= scm_i_big2dbl (x
);
4055 yy
= SCM_REAL_VALUE (y
);
4056 return (xx
< yy
? scm_from_double (xx
) : y
);
4058 else if (SCM_FRACTIONP (y
))
4063 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4065 else if (SCM_REALP (x
))
4067 if (SCM_I_INUMP (y
))
4069 double z
= SCM_I_INUM (y
);
4070 /* if x==NaN then "<" is false and we return NaN */
4071 return (z
< SCM_REAL_VALUE (x
)) ? scm_from_double (z
) : x
;
4073 else if (SCM_BIGP (y
))
4078 else if (SCM_REALP (y
))
4080 /* if x==NaN then our explicit check means we return NaN
4081 if y==NaN then "<" is false and we return NaN
4082 calling isnan is unavoidable, since it's the only way to know
4083 which of x or y causes any compares to be false */
4084 double xx
= SCM_REAL_VALUE (x
);
4085 return (isnan (xx
) || xx
< SCM_REAL_VALUE (y
)) ? x
: y
;
4087 else if (SCM_FRACTIONP (y
))
4089 double yy
= scm_i_fraction2double (y
);
4090 double xx
= SCM_REAL_VALUE (x
);
4091 return (yy
< xx
) ? scm_from_double (yy
) : x
;
4094 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4096 else if (SCM_FRACTIONP (x
))
4098 if (SCM_I_INUMP (y
))
4102 else if (SCM_BIGP (y
))
4106 else if (SCM_REALP (y
))
4108 double xx
= scm_i_fraction2double (x
);
4109 return (SCM_REAL_VALUE (y
) < xx
) ? y
: scm_from_double (xx
);
4111 else if (SCM_FRACTIONP (y
))
4116 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARGn
, s_min
);
4119 SCM_WTA_DISPATCH_2 (g_min
, x
, y
, SCM_ARG1
, s_min
);
4123 SCM_PRIMITIVE_GENERIC (scm_i_sum
, "+", 0, 2, 1,
4124 (SCM x
, SCM y
, SCM rest
),
4125 "Return the sum of all parameter values. Return 0 if called without\n"
4127 #define FUNC_NAME s_scm_i_sum
4129 while (!scm_is_null (rest
))
4130 { x
= scm_sum (x
, y
);
4132 rest
= scm_cdr (rest
);
4134 return scm_sum (x
, y
);
4138 #define s_sum s_scm_i_sum
4139 #define g_sum g_scm_i_sum
4142 scm_sum (SCM x
, SCM y
)
4144 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4146 if (SCM_NUMBERP (x
)) return x
;
4147 if (SCM_UNBNDP (x
)) return SCM_INUM0
;
4148 SCM_WTA_DISPATCH_1 (g_sum
, x
, SCM_ARG1
, s_sum
);
4151 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4153 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4155 long xx
= SCM_I_INUM (x
);
4156 long yy
= SCM_I_INUM (y
);
4157 long int z
= xx
+ yy
;
4158 return SCM_FIXABLE (z
) ? SCM_I_MAKINUM (z
) : scm_i_long2big (z
);
4160 else if (SCM_BIGP (y
))
4165 else if (SCM_REALP (y
))
4167 long int xx
= SCM_I_INUM (x
);
4168 return scm_from_double (xx
+ SCM_REAL_VALUE (y
));
4170 else if (SCM_COMPLEXP (y
))
4172 long int xx
= SCM_I_INUM (x
);
4173 return scm_c_make_rectangular (xx
+ SCM_COMPLEX_REAL (y
),
4174 SCM_COMPLEX_IMAG (y
));
4176 else if (SCM_FRACTIONP (y
))
4177 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4178 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4179 SCM_FRACTION_DENOMINATOR (y
));
4181 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4182 } else if (SCM_BIGP (x
))
4184 if (SCM_I_INUMP (y
))
4189 inum
= SCM_I_INUM (y
);
4192 bigsgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4195 SCM result
= scm_i_mkbig ();
4196 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), - inum
);
4197 scm_remember_upto_here_1 (x
);
4198 /* we know the result will have to be a bignum */
4201 return scm_i_normbig (result
);
4205 SCM result
= scm_i_mkbig ();
4206 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), inum
);
4207 scm_remember_upto_here_1 (x
);
4208 /* we know the result will have to be a bignum */
4211 return scm_i_normbig (result
);
4214 else if (SCM_BIGP (y
))
4216 SCM result
= scm_i_mkbig ();
4217 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4218 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4219 mpz_add (SCM_I_BIG_MPZ (result
),
4222 scm_remember_upto_here_2 (x
, y
);
4223 /* we know the result will have to be a bignum */
4226 return scm_i_normbig (result
);
4228 else if (SCM_REALP (y
))
4230 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) + SCM_REAL_VALUE (y
);
4231 scm_remember_upto_here_1 (x
);
4232 return scm_from_double (result
);
4234 else if (SCM_COMPLEXP (y
))
4236 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4237 + SCM_COMPLEX_REAL (y
));
4238 scm_remember_upto_here_1 (x
);
4239 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4241 else if (SCM_FRACTIONP (y
))
4242 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (y
),
4243 scm_product (x
, SCM_FRACTION_DENOMINATOR (y
))),
4244 SCM_FRACTION_DENOMINATOR (y
));
4246 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4248 else if (SCM_REALP (x
))
4250 if (SCM_I_INUMP (y
))
4251 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_I_INUM (y
));
4252 else if (SCM_BIGP (y
))
4254 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) + SCM_REAL_VALUE (x
);
4255 scm_remember_upto_here_1 (y
);
4256 return scm_from_double (result
);
4258 else if (SCM_REALP (y
))
4259 return scm_from_double (SCM_REAL_VALUE (x
) + SCM_REAL_VALUE (y
));
4260 else if (SCM_COMPLEXP (y
))
4261 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) + SCM_COMPLEX_REAL (y
),
4262 SCM_COMPLEX_IMAG (y
));
4263 else if (SCM_FRACTIONP (y
))
4264 return scm_from_double (SCM_REAL_VALUE (x
) + scm_i_fraction2double (y
));
4266 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4268 else if (SCM_COMPLEXP (x
))
4270 if (SCM_I_INUMP (y
))
4271 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_I_INUM (y
),
4272 SCM_COMPLEX_IMAG (x
));
4273 else if (SCM_BIGP (y
))
4275 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (y
))
4276 + SCM_COMPLEX_REAL (x
));
4277 scm_remember_upto_here_1 (y
);
4278 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (x
));
4280 else if (SCM_REALP (y
))
4281 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_REAL_VALUE (y
),
4282 SCM_COMPLEX_IMAG (x
));
4283 else if (SCM_COMPLEXP (y
))
4284 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + SCM_COMPLEX_REAL (y
),
4285 SCM_COMPLEX_IMAG (x
) + SCM_COMPLEX_IMAG (y
));
4286 else if (SCM_FRACTIONP (y
))
4287 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) + scm_i_fraction2double (y
),
4288 SCM_COMPLEX_IMAG (x
));
4290 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4292 else if (SCM_FRACTIONP (x
))
4294 if (SCM_I_INUMP (y
))
4295 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4296 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4297 SCM_FRACTION_DENOMINATOR (x
));
4298 else if (SCM_BIGP (y
))
4299 return scm_i_make_ratio (scm_sum (SCM_FRACTION_NUMERATOR (x
),
4300 scm_product (y
, SCM_FRACTION_DENOMINATOR (x
))),
4301 SCM_FRACTION_DENOMINATOR (x
));
4302 else if (SCM_REALP (y
))
4303 return scm_from_double (SCM_REAL_VALUE (y
) + scm_i_fraction2double (x
));
4304 else if (SCM_COMPLEXP (y
))
4305 return scm_c_make_rectangular (SCM_COMPLEX_REAL (y
) + scm_i_fraction2double (x
),
4306 SCM_COMPLEX_IMAG (y
));
4307 else if (SCM_FRACTIONP (y
))
4308 /* a/b + c/d = (ad + bc) / bd */
4309 return scm_i_make_ratio (scm_sum (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4310 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4311 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4313 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARGn
, s_sum
);
4316 SCM_WTA_DISPATCH_2 (g_sum
, x
, y
, SCM_ARG1
, s_sum
);
4320 SCM_DEFINE (scm_oneplus
, "1+", 1, 0, 0,
4322 "Return @math{@var{x}+1}.")
4323 #define FUNC_NAME s_scm_oneplus
4325 return scm_sum (x
, SCM_I_MAKINUM (1));
4330 SCM_PRIMITIVE_GENERIC (scm_i_difference
, "-", 0, 2, 1,
4331 (SCM x
, SCM y
, SCM rest
),
4332 "If called with one argument @var{z1}, -@var{z1} returned. Otherwise\n"
4333 "the sum of all but the first argument are subtracted from the first\n"
4335 #define FUNC_NAME s_scm_i_difference
4337 while (!scm_is_null (rest
))
4338 { x
= scm_difference (x
, y
);
4340 rest
= scm_cdr (rest
);
4342 return scm_difference (x
, y
);
4346 #define s_difference s_scm_i_difference
4347 #define g_difference g_scm_i_difference
4350 scm_difference (SCM x
, SCM y
)
4351 #define FUNC_NAME s_difference
4353 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4356 SCM_WTA_DISPATCH_0 (g_difference
, s_difference
);
4358 if (SCM_I_INUMP (x
))
4360 long xx
= -SCM_I_INUM (x
);
4361 if (SCM_FIXABLE (xx
))
4362 return SCM_I_MAKINUM (xx
);
4364 return scm_i_long2big (xx
);
4366 else if (SCM_BIGP (x
))
4367 /* Must scm_i_normbig here because -SCM_MOST_NEGATIVE_FIXNUM is a
4368 bignum, but negating that gives a fixnum. */
4369 return scm_i_normbig (scm_i_clonebig (x
, 0));
4370 else if (SCM_REALP (x
))
4371 return scm_from_double (-SCM_REAL_VALUE (x
));
4372 else if (SCM_COMPLEXP (x
))
4373 return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x
),
4374 -SCM_COMPLEX_IMAG (x
));
4375 else if (SCM_FRACTIONP (x
))
4376 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
), SCM_UNDEFINED
),
4377 SCM_FRACTION_DENOMINATOR (x
));
4379 SCM_WTA_DISPATCH_1 (g_difference
, x
, SCM_ARG1
, s_difference
);
4382 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4384 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4386 long int xx
= SCM_I_INUM (x
);
4387 long int yy
= SCM_I_INUM (y
);
4388 long int z
= xx
- yy
;
4389 if (SCM_FIXABLE (z
))
4390 return SCM_I_MAKINUM (z
);
4392 return scm_i_long2big (z
);
4394 else if (SCM_BIGP (y
))
4396 /* inum-x - big-y */
4397 long xx
= SCM_I_INUM (x
);
4400 return scm_i_clonebig (y
, 0);
4403 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4404 SCM result
= scm_i_mkbig ();
4407 mpz_ui_sub (SCM_I_BIG_MPZ (result
), xx
, SCM_I_BIG_MPZ (y
));
4410 /* x - y == -(y + -x) */
4411 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), -xx
);
4412 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
4414 scm_remember_upto_here_1 (y
);
4416 if ((xx
< 0 && (sgn_y
> 0)) || ((xx
> 0) && sgn_y
< 0))
4417 /* we know the result will have to be a bignum */
4420 return scm_i_normbig (result
);
4423 else if (SCM_REALP (y
))
4425 long int xx
= SCM_I_INUM (x
);
4426 return scm_from_double (xx
- SCM_REAL_VALUE (y
));
4428 else if (SCM_COMPLEXP (y
))
4430 long int xx
= SCM_I_INUM (x
);
4431 return scm_c_make_rectangular (xx
- SCM_COMPLEX_REAL (y
),
4432 - SCM_COMPLEX_IMAG (y
));
4434 else if (SCM_FRACTIONP (y
))
4435 /* a - b/c = (ac - b) / c */
4436 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4437 SCM_FRACTION_NUMERATOR (y
)),
4438 SCM_FRACTION_DENOMINATOR (y
));
4440 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4442 else if (SCM_BIGP (x
))
4444 if (SCM_I_INUMP (y
))
4446 /* big-x - inum-y */
4447 long yy
= SCM_I_INUM (y
);
4448 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4450 scm_remember_upto_here_1 (x
);
4452 return (SCM_FIXABLE (-yy
) ?
4453 SCM_I_MAKINUM (-yy
) : scm_from_long (-yy
));
4456 SCM result
= scm_i_mkbig ();
4459 mpz_sub_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), yy
);
4461 mpz_add_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), -yy
);
4462 scm_remember_upto_here_1 (x
);
4464 if ((sgn_x
< 0 && (yy
> 0)) || ((sgn_x
> 0) && yy
< 0))
4465 /* we know the result will have to be a bignum */
4468 return scm_i_normbig (result
);
4471 else if (SCM_BIGP (y
))
4473 int sgn_x
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4474 int sgn_y
= mpz_sgn (SCM_I_BIG_MPZ (y
));
4475 SCM result
= scm_i_mkbig ();
4476 mpz_sub (SCM_I_BIG_MPZ (result
),
4479 scm_remember_upto_here_2 (x
, y
);
4480 /* we know the result will have to be a bignum */
4481 if ((sgn_x
== 1) && (sgn_y
== -1))
4483 if ((sgn_x
== -1) && (sgn_y
== 1))
4485 return scm_i_normbig (result
);
4487 else if (SCM_REALP (y
))
4489 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) - SCM_REAL_VALUE (y
);
4490 scm_remember_upto_here_1 (x
);
4491 return scm_from_double (result
);
4493 else if (SCM_COMPLEXP (y
))
4495 double real_part
= (mpz_get_d (SCM_I_BIG_MPZ (x
))
4496 - SCM_COMPLEX_REAL (y
));
4497 scm_remember_upto_here_1 (x
);
4498 return scm_c_make_rectangular (real_part
, - SCM_COMPLEX_IMAG (y
));
4500 else if (SCM_FRACTIONP (y
))
4501 return scm_i_make_ratio (scm_difference (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4502 SCM_FRACTION_NUMERATOR (y
)),
4503 SCM_FRACTION_DENOMINATOR (y
));
4504 else SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4506 else if (SCM_REALP (x
))
4508 if (SCM_I_INUMP (y
))
4509 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_I_INUM (y
));
4510 else if (SCM_BIGP (y
))
4512 double result
= SCM_REAL_VALUE (x
) - mpz_get_d (SCM_I_BIG_MPZ (y
));
4513 scm_remember_upto_here_1 (x
);
4514 return scm_from_double (result
);
4516 else if (SCM_REALP (y
))
4517 return scm_from_double (SCM_REAL_VALUE (x
) - SCM_REAL_VALUE (y
));
4518 else if (SCM_COMPLEXP (y
))
4519 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) - SCM_COMPLEX_REAL (y
),
4520 -SCM_COMPLEX_IMAG (y
));
4521 else if (SCM_FRACTIONP (y
))
4522 return scm_from_double (SCM_REAL_VALUE (x
) - scm_i_fraction2double (y
));
4524 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4526 else if (SCM_COMPLEXP (x
))
4528 if (SCM_I_INUMP (y
))
4529 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_I_INUM (y
),
4530 SCM_COMPLEX_IMAG (x
));
4531 else if (SCM_BIGP (y
))
4533 double real_part
= (SCM_COMPLEX_REAL (x
)
4534 - mpz_get_d (SCM_I_BIG_MPZ (y
)));
4535 scm_remember_upto_here_1 (x
);
4536 return scm_c_make_rectangular (real_part
, SCM_COMPLEX_IMAG (y
));
4538 else if (SCM_REALP (y
))
4539 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_REAL_VALUE (y
),
4540 SCM_COMPLEX_IMAG (x
));
4541 else if (SCM_COMPLEXP (y
))
4542 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - SCM_COMPLEX_REAL (y
),
4543 SCM_COMPLEX_IMAG (x
) - SCM_COMPLEX_IMAG (y
));
4544 else if (SCM_FRACTIONP (y
))
4545 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) - scm_i_fraction2double (y
),
4546 SCM_COMPLEX_IMAG (x
));
4548 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4550 else if (SCM_FRACTIONP (x
))
4552 if (SCM_I_INUMP (y
))
4553 /* a/b - c = (a - cb) / b */
4554 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4555 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4556 SCM_FRACTION_DENOMINATOR (x
));
4557 else if (SCM_BIGP (y
))
4558 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x
),
4559 scm_product(y
, SCM_FRACTION_DENOMINATOR (x
))),
4560 SCM_FRACTION_DENOMINATOR (x
));
4561 else if (SCM_REALP (y
))
4562 return scm_from_double (scm_i_fraction2double (x
) - SCM_REAL_VALUE (y
));
4563 else if (SCM_COMPLEXP (y
))
4564 return scm_c_make_rectangular (scm_i_fraction2double (x
) - SCM_COMPLEX_REAL (y
),
4565 -SCM_COMPLEX_IMAG (y
));
4566 else if (SCM_FRACTIONP (y
))
4567 /* a/b - c/d = (ad - bc) / bd */
4568 return scm_i_make_ratio (scm_difference (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
4569 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
))),
4570 scm_product (SCM_FRACTION_DENOMINATOR (x
), SCM_FRACTION_DENOMINATOR (y
)));
4572 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARGn
, s_difference
);
4575 SCM_WTA_DISPATCH_2 (g_difference
, x
, y
, SCM_ARG1
, s_difference
);
4580 SCM_DEFINE (scm_oneminus
, "1-", 1, 0, 0,
4582 "Return @math{@var{x}-1}.")
4583 #define FUNC_NAME s_scm_oneminus
4585 return scm_difference (x
, SCM_I_MAKINUM (1));
4590 SCM_PRIMITIVE_GENERIC (scm_i_product
, "*", 0, 2, 1,
4591 (SCM x
, SCM y
, SCM rest
),
4592 "Return the product of all arguments. If called without arguments,\n"
4594 #define FUNC_NAME s_scm_i_product
4596 while (!scm_is_null (rest
))
4597 { x
= scm_product (x
, y
);
4599 rest
= scm_cdr (rest
);
4601 return scm_product (x
, y
);
4605 #define s_product s_scm_i_product
4606 #define g_product g_scm_i_product
4609 scm_product (SCM x
, SCM y
)
4611 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4614 return SCM_I_MAKINUM (1L);
4615 else if (SCM_NUMBERP (x
))
4618 SCM_WTA_DISPATCH_1 (g_product
, x
, SCM_ARG1
, s_product
);
4621 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4626 xx
= SCM_I_INUM (x
);
4630 case 0: return x
; break;
4631 case 1: return y
; break;
4634 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4636 long yy
= SCM_I_INUM (y
);
4638 SCM k
= SCM_I_MAKINUM (kk
);
4639 if ((kk
== SCM_I_INUM (k
)) && (kk
/ xx
== yy
))
4643 SCM result
= scm_i_long2big (xx
);
4644 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
), yy
);
4645 return scm_i_normbig (result
);
4648 else if (SCM_BIGP (y
))
4650 SCM result
= scm_i_mkbig ();
4651 mpz_mul_si (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (y
), xx
);
4652 scm_remember_upto_here_1 (y
);
4655 else if (SCM_REALP (y
))
4656 return scm_from_double (xx
* SCM_REAL_VALUE (y
));
4657 else if (SCM_COMPLEXP (y
))
4658 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4659 xx
* SCM_COMPLEX_IMAG (y
));
4660 else if (SCM_FRACTIONP (y
))
4661 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4662 SCM_FRACTION_DENOMINATOR (y
));
4664 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4666 else if (SCM_BIGP (x
))
4668 if (SCM_I_INUMP (y
))
4673 else if (SCM_BIGP (y
))
4675 SCM result
= scm_i_mkbig ();
4676 mpz_mul (SCM_I_BIG_MPZ (result
),
4679 scm_remember_upto_here_2 (x
, y
);
4682 else if (SCM_REALP (y
))
4684 double result
= mpz_get_d (SCM_I_BIG_MPZ (x
)) * SCM_REAL_VALUE (y
);
4685 scm_remember_upto_here_1 (x
);
4686 return scm_from_double (result
);
4688 else if (SCM_COMPLEXP (y
))
4690 double z
= mpz_get_d (SCM_I_BIG_MPZ (x
));
4691 scm_remember_upto_here_1 (x
);
4692 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (y
),
4693 z
* SCM_COMPLEX_IMAG (y
));
4695 else if (SCM_FRACTIONP (y
))
4696 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_NUMERATOR (y
)),
4697 SCM_FRACTION_DENOMINATOR (y
));
4699 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4701 else if (SCM_REALP (x
))
4703 if (SCM_I_INUMP (y
))
4705 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4706 if (scm_is_eq (y
, SCM_INUM0
))
4708 return scm_from_double (SCM_I_INUM (y
) * SCM_REAL_VALUE (x
));
4710 else if (SCM_BIGP (y
))
4712 double result
= mpz_get_d (SCM_I_BIG_MPZ (y
)) * SCM_REAL_VALUE (x
);
4713 scm_remember_upto_here_1 (y
);
4714 return scm_from_double (result
);
4716 else if (SCM_REALP (y
))
4717 return scm_from_double (SCM_REAL_VALUE (x
) * SCM_REAL_VALUE (y
));
4718 else if (SCM_COMPLEXP (y
))
4719 return scm_c_make_rectangular (SCM_REAL_VALUE (x
) * SCM_COMPLEX_REAL (y
),
4720 SCM_REAL_VALUE (x
) * SCM_COMPLEX_IMAG (y
));
4721 else if (SCM_FRACTIONP (y
))
4722 return scm_from_double (SCM_REAL_VALUE (x
) * scm_i_fraction2double (y
));
4724 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4726 else if (SCM_COMPLEXP (x
))
4728 if (SCM_I_INUMP (y
))
4730 /* inexact*exact0 => exact 0, per R5RS "Exactness" section */
4731 if (scm_is_eq (y
, SCM_INUM0
))
4733 return scm_c_make_rectangular (SCM_I_INUM (y
) * SCM_COMPLEX_REAL (x
),
4734 SCM_I_INUM (y
) * SCM_COMPLEX_IMAG (x
));
4736 else if (SCM_BIGP (y
))
4738 double z
= mpz_get_d (SCM_I_BIG_MPZ (y
));
4739 scm_remember_upto_here_1 (y
);
4740 return scm_c_make_rectangular (z
* SCM_COMPLEX_REAL (x
),
4741 z
* SCM_COMPLEX_IMAG (x
));
4743 else if (SCM_REALP (y
))
4744 return scm_c_make_rectangular (SCM_REAL_VALUE (y
) * SCM_COMPLEX_REAL (x
),
4745 SCM_REAL_VALUE (y
) * SCM_COMPLEX_IMAG (x
));
4746 else if (SCM_COMPLEXP (y
))
4748 return scm_c_make_rectangular (SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_REAL (y
)
4749 - SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_IMAG (y
),
4750 SCM_COMPLEX_REAL (x
) * SCM_COMPLEX_IMAG (y
)
4751 + SCM_COMPLEX_IMAG (x
) * SCM_COMPLEX_REAL (y
));
4753 else if (SCM_FRACTIONP (y
))
4755 double yy
= scm_i_fraction2double (y
);
4756 return scm_c_make_rectangular (yy
* SCM_COMPLEX_REAL (x
),
4757 yy
* SCM_COMPLEX_IMAG (x
));
4760 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4762 else if (SCM_FRACTIONP (x
))
4764 if (SCM_I_INUMP (y
))
4765 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4766 SCM_FRACTION_DENOMINATOR (x
));
4767 else if (SCM_BIGP (y
))
4768 return scm_i_make_ratio (scm_product (y
, SCM_FRACTION_NUMERATOR (x
)),
4769 SCM_FRACTION_DENOMINATOR (x
));
4770 else if (SCM_REALP (y
))
4771 return scm_from_double (scm_i_fraction2double (x
) * SCM_REAL_VALUE (y
));
4772 else if (SCM_COMPLEXP (y
))
4774 double xx
= scm_i_fraction2double (x
);
4775 return scm_c_make_rectangular (xx
* SCM_COMPLEX_REAL (y
),
4776 xx
* SCM_COMPLEX_IMAG (y
));
4778 else if (SCM_FRACTIONP (y
))
4779 /* a/b * c/d = ac / bd */
4780 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
),
4781 SCM_FRACTION_NUMERATOR (y
)),
4782 scm_product (SCM_FRACTION_DENOMINATOR (x
),
4783 SCM_FRACTION_DENOMINATOR (y
)));
4785 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARGn
, s_product
);
4788 SCM_WTA_DISPATCH_2 (g_product
, x
, y
, SCM_ARG1
, s_product
);
4791 #if ((defined (HAVE_ISINF) && defined (HAVE_ISNAN)) \
4792 || (defined (HAVE_FINITE) && defined (HAVE_ISNAN)))
4793 #define ALLOW_DIVIDE_BY_ZERO
4794 /* #define ALLOW_DIVIDE_BY_EXACT_ZERO */
4797 /* The code below for complex division is adapted from the GNU
4798 libstdc++, which adapted it from f2c's libF77, and is subject to
4801 /****************************************************************
4802 Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
4804 Permission to use, copy, modify, and distribute this software
4805 and its documentation for any purpose and without fee is hereby
4806 granted, provided that the above copyright notice appear in all
4807 copies and that both that the copyright notice and this
4808 permission notice and warranty disclaimer appear in supporting
4809 documentation, and that the names of AT&T Bell Laboratories or
4810 Bellcore or any of their entities not be used in advertising or
4811 publicity pertaining to distribution of the software without
4812 specific, written prior permission.
4814 AT&T and Bellcore disclaim all warranties with regard to this
4815 software, including all implied warranties of merchantability
4816 and fitness. In no event shall AT&T or Bellcore be liable for
4817 any special, indirect or consequential damages or any damages
4818 whatsoever resulting from loss of use, data or profits, whether
4819 in an action of contract, negligence or other tortious action,
4820 arising out of or in connection with the use or performance of
4822 ****************************************************************/
4824 SCM_PRIMITIVE_GENERIC (scm_i_divide
, "/", 0, 2, 1,
4825 (SCM x
, SCM y
, SCM rest
),
4826 "Divide the first argument by the product of the remaining\n"
4827 "arguments. If called with one argument @var{z1}, 1/@var{z1} is\n"
4829 #define FUNC_NAME s_scm_i_divide
4831 while (!scm_is_null (rest
))
4832 { x
= scm_divide (x
, y
);
4834 rest
= scm_cdr (rest
);
4836 return scm_divide (x
, y
);
4840 #define s_divide s_scm_i_divide
4841 #define g_divide g_scm_i_divide
4844 do_divide (SCM x
, SCM y
, int inexact
)
4845 #define FUNC_NAME s_divide
4849 if (SCM_UNLIKELY (SCM_UNBNDP (y
)))
4852 SCM_WTA_DISPATCH_0 (g_divide
, s_divide
);
4853 else if (SCM_I_INUMP (x
))
4855 long xx
= SCM_I_INUM (x
);
4856 if (xx
== 1 || xx
== -1)
4858 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4860 scm_num_overflow (s_divide
);
4865 return scm_from_double (1.0 / (double) xx
);
4866 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4869 else if (SCM_BIGP (x
))
4872 return scm_from_double (1.0 / scm_i_big2dbl (x
));
4873 else return scm_i_make_ratio (SCM_I_MAKINUM(1), x
);
4875 else if (SCM_REALP (x
))
4877 double xx
= SCM_REAL_VALUE (x
);
4878 #ifndef ALLOW_DIVIDE_BY_ZERO
4880 scm_num_overflow (s_divide
);
4883 return scm_from_double (1.0 / xx
);
4885 else if (SCM_COMPLEXP (x
))
4887 double r
= SCM_COMPLEX_REAL (x
);
4888 double i
= SCM_COMPLEX_IMAG (x
);
4889 if (fabs(r
) <= fabs(i
))
4892 double d
= i
* (1.0 + t
* t
);
4893 return scm_c_make_rectangular (t
/ d
, -1.0 / d
);
4898 double d
= r
* (1.0 + t
* t
);
4899 return scm_c_make_rectangular (1.0 / d
, -t
/ d
);
4902 else if (SCM_FRACTIONP (x
))
4903 return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x
),
4904 SCM_FRACTION_NUMERATOR (x
));
4906 SCM_WTA_DISPATCH_1 (g_divide
, x
, SCM_ARG1
, s_divide
);
4909 if (SCM_LIKELY (SCM_I_INUMP (x
)))
4911 long xx
= SCM_I_INUM (x
);
4912 if (SCM_LIKELY (SCM_I_INUMP (y
)))
4914 long yy
= SCM_I_INUM (y
);
4917 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4918 scm_num_overflow (s_divide
);
4920 return scm_from_double ((double) xx
/ (double) yy
);
4923 else if (xx
% yy
!= 0)
4926 return scm_from_double ((double) xx
/ (double) yy
);
4927 else return scm_i_make_ratio (x
, y
);
4932 if (SCM_FIXABLE (z
))
4933 return SCM_I_MAKINUM (z
);
4935 return scm_i_long2big (z
);
4938 else if (SCM_BIGP (y
))
4941 return scm_from_double ((double) xx
/ scm_i_big2dbl (y
));
4942 else return scm_i_make_ratio (x
, y
);
4944 else if (SCM_REALP (y
))
4946 double yy
= SCM_REAL_VALUE (y
);
4947 #ifndef ALLOW_DIVIDE_BY_ZERO
4949 scm_num_overflow (s_divide
);
4952 return scm_from_double ((double) xx
/ yy
);
4954 else if (SCM_COMPLEXP (y
))
4957 complex_div
: /* y _must_ be a complex number */
4959 double r
= SCM_COMPLEX_REAL (y
);
4960 double i
= SCM_COMPLEX_IMAG (y
);
4961 if (fabs(r
) <= fabs(i
))
4964 double d
= i
* (1.0 + t
* t
);
4965 return scm_c_make_rectangular ((a
* t
) / d
, -a
/ d
);
4970 double d
= r
* (1.0 + t
* t
);
4971 return scm_c_make_rectangular (a
/ d
, -(a
* t
) / d
);
4975 else if (SCM_FRACTIONP (y
))
4976 /* a / b/c = ac / b */
4977 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
4978 SCM_FRACTION_NUMERATOR (y
));
4980 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
4982 else if (SCM_BIGP (x
))
4984 if (SCM_I_INUMP (y
))
4986 long int yy
= SCM_I_INUM (y
);
4989 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
4990 scm_num_overflow (s_divide
);
4992 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
4993 scm_remember_upto_here_1 (x
);
4994 return (sgn
== 0) ? scm_nan () : scm_inf ();
5001 /* FIXME: HMM, what are the relative performance issues here?
5002 We need to test. Is it faster on average to test
5003 divisible_p, then perform whichever operation, or is it
5004 faster to perform the integer div opportunistically and
5005 switch to real if there's a remainder? For now we take the
5006 middle ground: test, then if divisible, use the faster div
5009 long abs_yy
= yy
< 0 ? -yy
: yy
;
5010 int divisible_p
= mpz_divisible_ui_p (SCM_I_BIG_MPZ (x
), abs_yy
);
5014 SCM result
= scm_i_mkbig ();
5015 mpz_divexact_ui (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (x
), abs_yy
);
5016 scm_remember_upto_here_1 (x
);
5018 mpz_neg (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (result
));
5019 return scm_i_normbig (result
);
5024 return scm_from_double (scm_i_big2dbl (x
) / (double) yy
);
5025 else return scm_i_make_ratio (x
, y
);
5029 else if (SCM_BIGP (y
))
5031 int y_is_zero
= (mpz_sgn (SCM_I_BIG_MPZ (y
)) == 0);
5034 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5035 scm_num_overflow (s_divide
);
5037 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (x
));
5038 scm_remember_upto_here_1 (x
);
5039 return (sgn
== 0) ? scm_nan () : scm_inf ();
5047 /* It's easily possible for the ratio x/y to fit a double
5048 but one or both x and y be too big to fit a double,
5049 hence the use of mpq_get_d rather than converting and
5052 *mpq_numref(q
) = *SCM_I_BIG_MPZ (x
);
5053 *mpq_denref(q
) = *SCM_I_BIG_MPZ (y
);
5054 return scm_from_double (mpq_get_d (q
));
5058 int divisible_p
= mpz_divisible_p (SCM_I_BIG_MPZ (x
),
5062 SCM result
= scm_i_mkbig ();
5063 mpz_divexact (SCM_I_BIG_MPZ (result
),
5066 scm_remember_upto_here_2 (x
, y
);
5067 return scm_i_normbig (result
);
5070 return scm_i_make_ratio (x
, y
);
5074 else if (SCM_REALP (y
))
5076 double yy
= SCM_REAL_VALUE (y
);
5077 #ifndef ALLOW_DIVIDE_BY_ZERO
5079 scm_num_overflow (s_divide
);
5082 return scm_from_double (scm_i_big2dbl (x
) / yy
);
5084 else if (SCM_COMPLEXP (y
))
5086 a
= scm_i_big2dbl (x
);
5089 else if (SCM_FRACTIONP (y
))
5090 return scm_i_make_ratio (scm_product (x
, SCM_FRACTION_DENOMINATOR (y
)),
5091 SCM_FRACTION_NUMERATOR (y
));
5093 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5095 else if (SCM_REALP (x
))
5097 double rx
= SCM_REAL_VALUE (x
);
5098 if (SCM_I_INUMP (y
))
5100 long int yy
= SCM_I_INUM (y
);
5101 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5103 scm_num_overflow (s_divide
);
5106 return scm_from_double (rx
/ (double) yy
);
5108 else if (SCM_BIGP (y
))
5110 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5111 scm_remember_upto_here_1 (y
);
5112 return scm_from_double (rx
/ dby
);
5114 else if (SCM_REALP (y
))
5116 double yy
= SCM_REAL_VALUE (y
);
5117 #ifndef ALLOW_DIVIDE_BY_ZERO
5119 scm_num_overflow (s_divide
);
5122 return scm_from_double (rx
/ yy
);
5124 else if (SCM_COMPLEXP (y
))
5129 else if (SCM_FRACTIONP (y
))
5130 return scm_from_double (rx
/ scm_i_fraction2double (y
));
5132 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5134 else if (SCM_COMPLEXP (x
))
5136 double rx
= SCM_COMPLEX_REAL (x
);
5137 double ix
= SCM_COMPLEX_IMAG (x
);
5138 if (SCM_I_INUMP (y
))
5140 long int yy
= SCM_I_INUM (y
);
5141 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5143 scm_num_overflow (s_divide
);
5148 return scm_c_make_rectangular (rx
/ d
, ix
/ d
);
5151 else if (SCM_BIGP (y
))
5153 double dby
= mpz_get_d (SCM_I_BIG_MPZ (y
));
5154 scm_remember_upto_here_1 (y
);
5155 return scm_c_make_rectangular (rx
/ dby
, ix
/ dby
);
5157 else if (SCM_REALP (y
))
5159 double yy
= SCM_REAL_VALUE (y
);
5160 #ifndef ALLOW_DIVIDE_BY_ZERO
5162 scm_num_overflow (s_divide
);
5165 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5167 else if (SCM_COMPLEXP (y
))
5169 double ry
= SCM_COMPLEX_REAL (y
);
5170 double iy
= SCM_COMPLEX_IMAG (y
);
5171 if (fabs(ry
) <= fabs(iy
))
5174 double d
= iy
* (1.0 + t
* t
);
5175 return scm_c_make_rectangular ((rx
* t
+ ix
) / d
, (ix
* t
- rx
) / d
);
5180 double d
= ry
* (1.0 + t
* t
);
5181 return scm_c_make_rectangular ((rx
+ ix
* t
) / d
, (ix
- rx
* t
) / d
);
5184 else if (SCM_FRACTIONP (y
))
5186 double yy
= scm_i_fraction2double (y
);
5187 return scm_c_make_rectangular (rx
/ yy
, ix
/ yy
);
5190 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5192 else if (SCM_FRACTIONP (x
))
5194 if (SCM_I_INUMP (y
))
5196 long int yy
= SCM_I_INUM (y
);
5197 #ifndef ALLOW_DIVIDE_BY_EXACT_ZERO
5199 scm_num_overflow (s_divide
);
5202 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5203 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5205 else if (SCM_BIGP (y
))
5207 return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x
),
5208 scm_product (SCM_FRACTION_DENOMINATOR (x
), y
));
5210 else if (SCM_REALP (y
))
5212 double yy
= SCM_REAL_VALUE (y
);
5213 #ifndef ALLOW_DIVIDE_BY_ZERO
5215 scm_num_overflow (s_divide
);
5218 return scm_from_double (scm_i_fraction2double (x
) / yy
);
5220 else if (SCM_COMPLEXP (y
))
5222 a
= scm_i_fraction2double (x
);
5225 else if (SCM_FRACTIONP (y
))
5226 return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x
), SCM_FRACTION_DENOMINATOR (y
)),
5227 scm_product (SCM_FRACTION_NUMERATOR (y
), SCM_FRACTION_DENOMINATOR (x
)));
5229 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARGn
, s_divide
);
5232 SCM_WTA_DISPATCH_2 (g_divide
, x
, y
, SCM_ARG1
, s_divide
);
5236 scm_divide (SCM x
, SCM y
)
5238 return do_divide (x
, y
, 0);
5241 static SCM
scm_divide2real (SCM x
, SCM y
)
5243 return do_divide (x
, y
, 1);
5249 scm_c_truncate (double x
)
5260 /* scm_c_round is done using floor(x+0.5) to round to nearest and with
5261 half-way case (ie. when x is an integer plus 0.5) going upwards.
5262 Then half-way cases are identified and adjusted down if the
5263 round-upwards didn't give the desired even integer.
5265 "plus_half == result" identifies a half-way case. If plus_half, which is
5266 x + 0.5, is an integer then x must be an integer plus 0.5.
5268 An odd "result" value is identified with result/2 != floor(result/2).
5269 This is done with plus_half, since that value is ready for use sooner in
5270 a pipelined cpu, and we're already requiring plus_half == result.
5272 Note however that we need to be careful when x is big and already an
5273 integer. In that case "x+0.5" may round to an adjacent integer, causing
5274 us to return such a value, incorrectly. For instance if the hardware is
5275 in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF
5276 (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value
5277 returned. Or if the hardware is in round-upwards mode, then other bigger
5278 values like say x == 2^128 will see x+0.5 rounding up to the next higher
5279 representable value, 2^128+2^76 (or whatever), again incorrect.
5281 These bad roundings of x+0.5 are avoided by testing at the start whether
5282 x is already an integer. If it is then clearly that's the desired result
5283 already. And if it's not then the exponent must be small enough to allow
5284 an 0.5 to be represented, and hence added without a bad rounding. */
5287 scm_c_round (double x
)
5289 double plus_half
, result
;
5294 plus_half
= x
+ 0.5;
5295 result
= floor (plus_half
);
5296 /* Adjust so that the rounding is towards even. */
5297 return ((plus_half
== result
&& plus_half
/ 2 != floor (plus_half
/ 2))
5302 SCM_DEFINE (scm_truncate_number
, "truncate", 1, 0, 0,
5304 "Round the number @var{x} towards zero.")
5305 #define FUNC_NAME s_scm_truncate_number
5307 if (scm_is_false (scm_negative_p (x
)))
5308 return scm_floor (x
);
5310 return scm_ceiling (x
);
5314 static SCM exactly_one_half
;
5316 SCM_DEFINE (scm_round_number
, "round", 1, 0, 0,
5318 "Round the number @var{x} towards the nearest integer. "
5319 "When it is exactly halfway between two integers, "
5320 "round towards the even one.")
5321 #define FUNC_NAME s_scm_round_number
5323 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5325 else if (SCM_REALP (x
))
5326 return scm_from_double (scm_c_round (SCM_REAL_VALUE (x
)));
5329 /* OPTIMIZE-ME: Fraction case could be done more efficiently by a
5330 single quotient+remainder division then examining to see which way
5331 the rounding should go. */
5332 SCM plus_half
= scm_sum (x
, exactly_one_half
);
5333 SCM result
= scm_floor (plus_half
);
5334 /* Adjust so that the rounding is towards even. */
5335 if (scm_is_true (scm_num_eq_p (plus_half
, result
))
5336 && scm_is_true (scm_odd_p (result
)))
5337 return scm_difference (result
, SCM_I_MAKINUM (1));
5344 SCM_PRIMITIVE_GENERIC (scm_floor
, "floor", 1, 0, 0,
5346 "Round the number @var{x} towards minus infinity.")
5347 #define FUNC_NAME s_scm_floor
5349 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5351 else if (SCM_REALP (x
))
5352 return scm_from_double (floor (SCM_REAL_VALUE (x
)));
5353 else if (SCM_FRACTIONP (x
))
5355 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5356 SCM_FRACTION_DENOMINATOR (x
));
5357 if (scm_is_false (scm_negative_p (x
)))
5359 /* For positive x, rounding towards zero is correct. */
5364 /* For negative x, we need to return q-1 unless x is an
5365 integer. But fractions are never integer, per our
5367 return scm_difference (q
, SCM_I_MAKINUM (1));
5371 SCM_WTA_DISPATCH_1 (g_scm_floor
, x
, 1, s_scm_floor
);
5375 SCM_PRIMITIVE_GENERIC (scm_ceiling
, "ceiling", 1, 0, 0,
5377 "Round the number @var{x} towards infinity.")
5378 #define FUNC_NAME s_scm_ceiling
5380 if (SCM_I_INUMP (x
) || SCM_BIGP (x
))
5382 else if (SCM_REALP (x
))
5383 return scm_from_double (ceil (SCM_REAL_VALUE (x
)));
5384 else if (SCM_FRACTIONP (x
))
5386 SCM q
= scm_quotient (SCM_FRACTION_NUMERATOR (x
),
5387 SCM_FRACTION_DENOMINATOR (x
));
5388 if (scm_is_false (scm_positive_p (x
)))
5390 /* For negative x, rounding towards zero is correct. */
5395 /* For positive x, we need to return q+1 unless x is an
5396 integer. But fractions are never integer, per our
5398 return scm_sum (q
, SCM_I_MAKINUM (1));
5402 SCM_WTA_DISPATCH_1 (g_scm_ceiling
, x
, 1, s_scm_ceiling
);
5406 /* sin/cos/tan/asin/acos/atan
5407 sinh/cosh/tanh/asinh/acosh/atanh
5408 Derived from "Transcen.scm", Complex trancendental functions for SCM.
5409 Written by Jerry D. Hedden, (C) FSF.
5410 See the file `COPYING' for terms applying to this program. */
5412 SCM_DEFINE (scm_expt
, "expt", 2, 0, 0,
5414 "Return @var{x} raised to the power of @var{y}.")
5415 #define FUNC_NAME s_scm_expt
5417 if (scm_is_true (scm_exact_p (x
)) && scm_is_integer (y
))
5418 return scm_integer_expt (x
, y
);
5419 else if (scm_is_real (x
) && scm_is_real (y
) && scm_to_double (x
) >= 0.0)
5421 return scm_from_double (pow (scm_to_double (x
), scm_to_double (y
)));
5424 return scm_exp (scm_product (scm_log (x
), y
));
5428 SCM_PRIMITIVE_GENERIC (scm_sin
, "sin", 1, 0, 0,
5430 "Compute the sine of @var{z}.")
5431 #define FUNC_NAME s_scm_sin
5433 if (scm_is_real (z
))
5434 return scm_from_double (sin (scm_to_double (z
)));
5435 else if (SCM_COMPLEXP (z
))
5437 x
= SCM_COMPLEX_REAL (z
);
5438 y
= SCM_COMPLEX_IMAG (z
);
5439 return scm_c_make_rectangular (sin (x
) * cosh (y
),
5440 cos (x
) * sinh (y
));
5443 SCM_WTA_DISPATCH_1 (g_scm_sin
, z
, 1, s_scm_sin
);
5447 SCM_PRIMITIVE_GENERIC (scm_cos
, "cos", 1, 0, 0,
5449 "Compute the cosine of @var{z}.")
5450 #define FUNC_NAME s_scm_cos
5452 if (scm_is_real (z
))
5453 return scm_from_double (cos (scm_to_double (z
)));
5454 else if (SCM_COMPLEXP (z
))
5456 x
= SCM_COMPLEX_REAL (z
);
5457 y
= SCM_COMPLEX_IMAG (z
);
5458 return scm_c_make_rectangular (cos (x
) * cosh (y
),
5459 -sin (x
) * sinh (y
));
5462 SCM_WTA_DISPATCH_1 (g_scm_cos
, z
, 1, s_scm_cos
);
5466 SCM_PRIMITIVE_GENERIC (scm_tan
, "tan", 1, 0, 0,
5468 "Compute the tangent of @var{z}.")
5469 #define FUNC_NAME s_scm_tan
5471 if (scm_is_real (z
))
5472 return scm_from_double (tan (scm_to_double (z
)));
5473 else if (SCM_COMPLEXP (z
))
5475 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5476 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5477 w
= cos (x
) + cosh (y
);
5478 #ifndef ALLOW_DIVIDE_BY_ZERO
5480 scm_num_overflow (s_scm_tan
);
5482 return scm_c_make_rectangular (sin (x
) / w
, sinh (y
) / w
);
5485 SCM_WTA_DISPATCH_1 (g_scm_tan
, z
, 1, s_scm_tan
);
5489 SCM_PRIMITIVE_GENERIC (scm_sinh
, "sinh", 1, 0, 0,
5491 "Compute the hyperbolic sine of @var{z}.")
5492 #define FUNC_NAME s_scm_sinh
5494 if (scm_is_real (z
))
5495 return scm_from_double (sinh (scm_to_double (z
)));
5496 else if (SCM_COMPLEXP (z
))
5498 x
= SCM_COMPLEX_REAL (z
);
5499 y
= SCM_COMPLEX_IMAG (z
);
5500 return scm_c_make_rectangular (sinh (x
) * cos (y
),
5501 cosh (x
) * sin (y
));
5504 SCM_WTA_DISPATCH_1 (g_scm_sinh
, z
, 1, s_scm_sinh
);
5508 SCM_PRIMITIVE_GENERIC (scm_cosh
, "cosh", 1, 0, 0,
5510 "Compute the hyperbolic cosine of @var{z}.")
5511 #define FUNC_NAME s_scm_cosh
5513 if (scm_is_real (z
))
5514 return scm_from_double (cosh (scm_to_double (z
)));
5515 else if (SCM_COMPLEXP (z
))
5517 x
= SCM_COMPLEX_REAL (z
);
5518 y
= SCM_COMPLEX_IMAG (z
);
5519 return scm_c_make_rectangular (cosh (x
) * cos (y
),
5520 sinh (x
) * sin (y
));
5523 SCM_WTA_DISPATCH_1 (g_scm_cosh
, z
, 1, s_scm_cosh
);
5527 SCM_PRIMITIVE_GENERIC (scm_tanh
, "tanh", 1, 0, 0,
5529 "Compute the hyperbolic tangent of @var{z}.")
5530 #define FUNC_NAME s_scm_tanh
5532 if (scm_is_real (z
))
5533 return scm_from_double (tanh (scm_to_double (z
)));
5534 else if (SCM_COMPLEXP (z
))
5536 x
= 2.0 * SCM_COMPLEX_REAL (z
);
5537 y
= 2.0 * SCM_COMPLEX_IMAG (z
);
5538 w
= cosh (x
) + cos (y
);
5539 #ifndef ALLOW_DIVIDE_BY_ZERO
5541 scm_num_overflow (s_scm_tanh
);
5543 return scm_c_make_rectangular (sinh (x
) / w
, sin (y
) / w
);
5546 SCM_WTA_DISPATCH_1 (g_scm_tanh
, z
, 1, s_scm_tanh
);
5550 SCM_PRIMITIVE_GENERIC (scm_asin
, "asin", 1, 0, 0,
5552 "Compute the arc sine of @var{z}.")
5553 #define FUNC_NAME s_scm_asin
5555 if (scm_is_real (z
))
5557 double w
= scm_to_double (z
);
5558 if (w
>= -1.0 && w
<= 1.0)
5559 return scm_from_double (asin (w
));
5561 return scm_product (scm_c_make_rectangular (0, -1),
5562 scm_sys_asinh (scm_c_make_rectangular (0, w
)));
5564 else if (SCM_COMPLEXP (z
))
5566 x
= SCM_COMPLEX_REAL (z
);
5567 y
= SCM_COMPLEX_IMAG (z
);
5568 return scm_product (scm_c_make_rectangular (0, -1),
5569 scm_sys_asinh (scm_c_make_rectangular (-y
, x
)));
5572 SCM_WTA_DISPATCH_1 (g_scm_asin
, z
, 1, s_scm_asin
);
5576 SCM_PRIMITIVE_GENERIC (scm_acos
, "acos", 1, 0, 0,
5578 "Compute the arc cosine of @var{z}.")
5579 #define FUNC_NAME s_scm_acos
5581 if (scm_is_real (z
))
5583 double w
= scm_to_double (z
);
5584 if (w
>= -1.0 && w
<= 1.0)
5585 return scm_from_double (acos (w
));
5587 return scm_sum (scm_from_double (acos (0.0)),
5588 scm_product (scm_c_make_rectangular (0, 1),
5589 scm_sys_asinh (scm_c_make_rectangular (0, w
))));
5591 else if (SCM_COMPLEXP (z
))
5593 x
= SCM_COMPLEX_REAL (z
);
5594 y
= SCM_COMPLEX_IMAG (z
);
5595 return scm_sum (scm_from_double (acos (0.0)),
5596 scm_product (scm_c_make_rectangular (0, 1),
5597 scm_sys_asinh (scm_c_make_rectangular (-y
, x
))));
5600 SCM_WTA_DISPATCH_1 (g_scm_acos
, z
, 1, s_scm_acos
);
5604 SCM_PRIMITIVE_GENERIC (scm_atan
, "atan", 1, 1, 0,
5606 "With one argument, compute the arc tangent of @var{z}.\n"
5607 "If @var{y} is present, compute the arc tangent of @var{z}/@var{y},\n"
5608 "using the sign of @var{z} and @var{y} to determine the quadrant.")
5609 #define FUNC_NAME s_scm_atan
5613 if (scm_is_real (z
))
5614 return scm_from_double (atan (scm_to_double (z
)));
5615 else if (SCM_COMPLEXP (z
))
5618 v
= SCM_COMPLEX_REAL (z
);
5619 w
= SCM_COMPLEX_IMAG (z
);
5620 return scm_divide (scm_log (scm_divide (scm_c_make_rectangular (v
, w
- 1.0),
5621 scm_c_make_rectangular (v
, w
+ 1.0))),
5622 scm_c_make_rectangular (0, 2));
5625 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5627 else if (scm_is_real (z
))
5629 if (scm_is_real (y
))
5630 return scm_from_double (atan2 (scm_to_double (z
), scm_to_double (y
)));
5632 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG2
, s_scm_atan
);
5635 SCM_WTA_DISPATCH_2 (g_scm_atan
, z
, y
, SCM_ARG1
, s_scm_atan
);
5639 SCM_PRIMITIVE_GENERIC (scm_sys_asinh
, "asinh", 1, 0, 0,
5641 "Compute the inverse hyperbolic sine of @var{z}.")
5642 #define FUNC_NAME s_scm_sys_asinh
5644 if (scm_is_real (z
))
5645 return scm_from_double (asinh (scm_to_double (z
)));
5646 else if (scm_is_number (z
))
5647 return scm_log (scm_sum (z
,
5648 scm_sqrt (scm_sum (scm_product (z
, z
),
5649 SCM_I_MAKINUM (1)))));
5651 SCM_WTA_DISPATCH_1 (g_scm_sys_asinh
, z
, 1, s_scm_sys_asinh
);
5655 SCM_PRIMITIVE_GENERIC (scm_sys_acosh
, "acosh", 1, 0, 0,
5657 "Compute the inverse hyperbolic cosine of @var{z}.")
5658 #define FUNC_NAME s_scm_sys_acosh
5660 if (scm_is_real (z
) && scm_to_double (z
) >= 1.0)
5661 return scm_from_double (acosh (scm_to_double (z
)));
5662 else if (scm_is_number (z
))
5663 return scm_log (scm_sum (z
,
5664 scm_sqrt (scm_difference (scm_product (z
, z
),
5665 SCM_I_MAKINUM (1)))));
5667 SCM_WTA_DISPATCH_1 (g_scm_sys_acosh
, z
, 1, s_scm_sys_acosh
);
5671 SCM_PRIMITIVE_GENERIC (scm_sys_atanh
, "atanh", 1, 0, 0,
5673 "Compute the inverse hyperbolic tangent of @var{z}.")
5674 #define FUNC_NAME s_scm_sys_atanh
5676 if (scm_is_real (z
) && scm_to_double (z
) >= -1.0 && scm_to_double (z
) <= 1.0)
5677 return scm_from_double (atanh (scm_to_double (z
)));
5678 else if (scm_is_number (z
))
5679 return scm_divide (scm_log (scm_divide (scm_sum (SCM_I_MAKINUM (1), z
),
5680 scm_difference (SCM_I_MAKINUM (1), z
))),
5683 SCM_WTA_DISPATCH_1 (g_scm_sys_atanh
, z
, 1, s_scm_sys_atanh
);
5688 scm_c_make_rectangular (double re
, double im
)
5691 return scm_from_double (re
);
5695 SCM_NEWSMOB (z
, scm_tc16_complex
,
5696 scm_gc_malloc_pointerless (sizeof (scm_t_complex
),
5698 SCM_COMPLEX_REAL (z
) = re
;
5699 SCM_COMPLEX_IMAG (z
) = im
;
5704 SCM_DEFINE (scm_make_rectangular
, "make-rectangular", 2, 0, 0,
5705 (SCM real_part
, SCM imaginary_part
),
5706 "Return a complex number constructed of the given @var{real-part} "
5707 "and @var{imaginary-part} parts.")
5708 #define FUNC_NAME s_scm_make_rectangular
5710 SCM_ASSERT_TYPE (scm_is_real (real_part
), real_part
,
5711 SCM_ARG1
, FUNC_NAME
, "real");
5712 SCM_ASSERT_TYPE (scm_is_real (imaginary_part
), imaginary_part
,
5713 SCM_ARG2
, FUNC_NAME
, "real");
5714 return scm_c_make_rectangular (scm_to_double (real_part
),
5715 scm_to_double (imaginary_part
));
5720 scm_c_make_polar (double mag
, double ang
)
5724 /* The sincos(3) function is undocumented an broken on Tru64. Thus we only
5725 use it on Glibc-based systems that have it (it's a GNU extension). See
5726 http://lists.gnu.org/archive/html/guile-user/2009-04/msg00033.html for
5728 #if (defined HAVE_SINCOS) && (defined __GLIBC__) && (defined _GNU_SOURCE)
5729 sincos (ang
, &s
, &c
);
5734 return scm_c_make_rectangular (mag
* c
, mag
* s
);
5737 SCM_DEFINE (scm_make_polar
, "make-polar", 2, 0, 0,
5739 "Return the complex number @var{x} * e^(i * @var{y}).")
5740 #define FUNC_NAME s_scm_make_polar
5742 SCM_ASSERT_TYPE (scm_is_real (x
), x
, SCM_ARG1
, FUNC_NAME
, "real");
5743 SCM_ASSERT_TYPE (scm_is_real (y
), y
, SCM_ARG2
, FUNC_NAME
, "real");
5744 return scm_c_make_polar (scm_to_double (x
), scm_to_double (y
));
5749 SCM_GPROC (s_real_part
, "real-part", 1, 0, 0, scm_real_part
, g_real_part
);
5750 /* "Return the real part of the number @var{z}."
5753 scm_real_part (SCM z
)
5755 if (SCM_I_INUMP (z
))
5757 else if (SCM_BIGP (z
))
5759 else if (SCM_REALP (z
))
5761 else if (SCM_COMPLEXP (z
))
5762 return scm_from_double (SCM_COMPLEX_REAL (z
));
5763 else if (SCM_FRACTIONP (z
))
5766 SCM_WTA_DISPATCH_1 (g_real_part
, z
, SCM_ARG1
, s_real_part
);
5770 SCM_GPROC (s_imag_part
, "imag-part", 1, 0, 0, scm_imag_part
, g_imag_part
);
5771 /* "Return the imaginary part of the number @var{z}."
5774 scm_imag_part (SCM z
)
5776 if (SCM_I_INUMP (z
))
5778 else if (SCM_BIGP (z
))
5780 else if (SCM_REALP (z
))
5782 else if (SCM_COMPLEXP (z
))
5783 return scm_from_double (SCM_COMPLEX_IMAG (z
));
5784 else if (SCM_FRACTIONP (z
))
5787 SCM_WTA_DISPATCH_1 (g_imag_part
, z
, SCM_ARG1
, s_imag_part
);
5790 SCM_GPROC (s_numerator
, "numerator", 1, 0, 0, scm_numerator
, g_numerator
);
5791 /* "Return the numerator of the number @var{z}."
5794 scm_numerator (SCM z
)
5796 if (SCM_I_INUMP (z
))
5798 else if (SCM_BIGP (z
))
5800 else if (SCM_FRACTIONP (z
))
5801 return SCM_FRACTION_NUMERATOR (z
);
5802 else if (SCM_REALP (z
))
5803 return scm_exact_to_inexact (scm_numerator (scm_inexact_to_exact (z
)));
5805 SCM_WTA_DISPATCH_1 (g_numerator
, z
, SCM_ARG1
, s_numerator
);
5809 SCM_GPROC (s_denominator
, "denominator", 1, 0, 0, scm_denominator
, g_denominator
);
5810 /* "Return the denominator of the number @var{z}."
5813 scm_denominator (SCM z
)
5815 if (SCM_I_INUMP (z
))
5816 return SCM_I_MAKINUM (1);
5817 else if (SCM_BIGP (z
))
5818 return SCM_I_MAKINUM (1);
5819 else if (SCM_FRACTIONP (z
))
5820 return SCM_FRACTION_DENOMINATOR (z
);
5821 else if (SCM_REALP (z
))
5822 return scm_exact_to_inexact (scm_denominator (scm_inexact_to_exact (z
)));
5824 SCM_WTA_DISPATCH_1 (g_denominator
, z
, SCM_ARG1
, s_denominator
);
5827 SCM_GPROC (s_magnitude
, "magnitude", 1, 0, 0, scm_magnitude
, g_magnitude
);
5828 /* "Return the magnitude of the number @var{z}. This is the same as\n"
5829 * "@code{abs} for real arguments, but also allows complex numbers."
5832 scm_magnitude (SCM z
)
5834 if (SCM_I_INUMP (z
))
5836 long int zz
= SCM_I_INUM (z
);
5839 else if (SCM_POSFIXABLE (-zz
))
5840 return SCM_I_MAKINUM (-zz
);
5842 return scm_i_long2big (-zz
);
5844 else if (SCM_BIGP (z
))
5846 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5847 scm_remember_upto_here_1 (z
);
5849 return scm_i_clonebig (z
, 0);
5853 else if (SCM_REALP (z
))
5854 return scm_from_double (fabs (SCM_REAL_VALUE (z
)));
5855 else if (SCM_COMPLEXP (z
))
5856 return scm_from_double (hypot (SCM_COMPLEX_REAL (z
), SCM_COMPLEX_IMAG (z
)));
5857 else if (SCM_FRACTIONP (z
))
5859 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5861 return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z
), SCM_UNDEFINED
),
5862 SCM_FRACTION_DENOMINATOR (z
));
5865 SCM_WTA_DISPATCH_1 (g_magnitude
, z
, SCM_ARG1
, s_magnitude
);
5869 SCM_GPROC (s_angle
, "angle", 1, 0, 0, scm_angle
, g_angle
);
5870 /* "Return the angle of the complex number @var{z}."
5875 /* atan(0,-1) is pi and it'd be possible to have that as a constant like
5876 flo0 to save allocating a new flonum with scm_from_double each time.
5877 But if atan2 follows the floating point rounding mode, then the value
5878 is not a constant. Maybe it'd be close enough though. */
5879 if (SCM_I_INUMP (z
))
5881 if (SCM_I_INUM (z
) >= 0)
5884 return scm_from_double (atan2 (0.0, -1.0));
5886 else if (SCM_BIGP (z
))
5888 int sgn
= mpz_sgn (SCM_I_BIG_MPZ (z
));
5889 scm_remember_upto_here_1 (z
);
5891 return scm_from_double (atan2 (0.0, -1.0));
5895 else if (SCM_REALP (z
))
5897 if (SCM_REAL_VALUE (z
) >= 0)
5900 return scm_from_double (atan2 (0.0, -1.0));
5902 else if (SCM_COMPLEXP (z
))
5903 return scm_from_double (atan2 (SCM_COMPLEX_IMAG (z
), SCM_COMPLEX_REAL (z
)));
5904 else if (SCM_FRACTIONP (z
))
5906 if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z
))))
5908 else return scm_from_double (atan2 (0.0, -1.0));
5911 SCM_WTA_DISPATCH_1 (g_angle
, z
, SCM_ARG1
, s_angle
);
5915 SCM_GPROC (s_exact_to_inexact
, "exact->inexact", 1, 0, 0, scm_exact_to_inexact
, g_exact_to_inexact
);
5916 /* Convert the number @var{x} to its inexact representation.\n"
5919 scm_exact_to_inexact (SCM z
)
5921 if (SCM_I_INUMP (z
))
5922 return scm_from_double ((double) SCM_I_INUM (z
));
5923 else if (SCM_BIGP (z
))
5924 return scm_from_double (scm_i_big2dbl (z
));
5925 else if (SCM_FRACTIONP (z
))
5926 return scm_from_double (scm_i_fraction2double (z
));
5927 else if (SCM_INEXACTP (z
))
5930 SCM_WTA_DISPATCH_1 (g_exact_to_inexact
, z
, 1, s_exact_to_inexact
);
5934 SCM_DEFINE (scm_inexact_to_exact
, "inexact->exact", 1, 0, 0,
5936 "Return an exact number that is numerically closest to @var{z}.")
5937 #define FUNC_NAME s_scm_inexact_to_exact
5939 if (SCM_I_INUMP (z
))
5941 else if (SCM_BIGP (z
))
5943 else if (SCM_REALP (z
))
5945 if (isinf (SCM_REAL_VALUE (z
)) || isnan (SCM_REAL_VALUE (z
)))
5946 SCM_OUT_OF_RANGE (1, z
);
5953 mpq_set_d (frac
, SCM_REAL_VALUE (z
));
5954 q
= scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac
)),
5955 scm_i_mpz2num (mpq_denref (frac
)));
5957 /* When scm_i_make_ratio throws, we leak the memory allocated
5964 else if (SCM_FRACTIONP (z
))
5967 SCM_WRONG_TYPE_ARG (1, z
);
5971 SCM_DEFINE (scm_rationalize
, "rationalize", 2, 0, 0,
5973 "Returns the @emph{simplest} rational number differing\n"
5974 "from @var{x} by no more than @var{eps}.\n"
5976 "As required by @acronym{R5RS}, @code{rationalize} only returns an\n"
5977 "exact result when both its arguments are exact. Thus, you might need\n"
5978 "to use @code{inexact->exact} on the arguments.\n"
5981 "(rationalize (inexact->exact 1.2) 1/100)\n"
5984 #define FUNC_NAME s_scm_rationalize
5986 if (SCM_I_INUMP (x
))
5988 else if (SCM_BIGP (x
))
5990 else if ((SCM_REALP (x
)) || SCM_FRACTIONP (x
))
5992 /* Use continued fractions to find closest ratio. All
5993 arithmetic is done with exact numbers.
5996 SCM ex
= scm_inexact_to_exact (x
);
5997 SCM int_part
= scm_floor (ex
);
5998 SCM tt
= SCM_I_MAKINUM (1);
5999 SCM a1
= SCM_I_MAKINUM (0), a2
= SCM_I_MAKINUM (1), a
= SCM_I_MAKINUM (0);
6000 SCM b1
= SCM_I_MAKINUM (1), b2
= SCM_I_MAKINUM (0), b
= SCM_I_MAKINUM (0);
6004 if (scm_is_true (scm_num_eq_p (ex
, int_part
)))
6007 ex
= scm_difference (ex
, int_part
); /* x = x-int_part */
6008 rx
= scm_divide (ex
, SCM_UNDEFINED
); /* rx = 1/x */
6010 /* We stop after a million iterations just to be absolutely sure
6011 that we don't go into an infinite loop. The process normally
6012 converges after less than a dozen iterations.
6015 eps
= scm_abs (eps
);
6016 while (++i
< 1000000)
6018 a
= scm_sum (scm_product (a1
, tt
), a2
); /* a = a1*tt + a2 */
6019 b
= scm_sum (scm_product (b1
, tt
), b2
); /* b = b1*tt + b2 */
6020 if (scm_is_false (scm_zero_p (b
)) && /* b != 0 */
6022 (scm_gr_p (scm_abs (scm_difference (ex
, scm_divide (a
, b
))),
6023 eps
))) /* abs(x-a/b) <= eps */
6025 SCM res
= scm_sum (int_part
, scm_divide (a
, b
));
6026 if (scm_is_false (scm_exact_p (x
))
6027 || scm_is_false (scm_exact_p (eps
)))
6028 return scm_exact_to_inexact (res
);
6032 rx
= scm_divide (scm_difference (rx
, tt
), /* rx = 1/(rx - tt) */
6034 tt
= scm_floor (rx
); /* tt = floor (rx) */
6040 scm_num_overflow (s_scm_rationalize
);
6043 SCM_WRONG_TYPE_ARG (1, x
);
6047 /* conversion functions */
6050 scm_is_integer (SCM val
)
6052 return scm_is_true (scm_integer_p (val
));
6056 scm_is_signed_integer (SCM val
, scm_t_intmax min
, scm_t_intmax max
)
6058 if (SCM_I_INUMP (val
))
6060 scm_t_signed_bits n
= SCM_I_INUM (val
);
6061 return n
>= min
&& n
<= max
;
6063 else if (SCM_BIGP (val
))
6065 if (min
>= SCM_MOST_NEGATIVE_FIXNUM
&& max
<= SCM_MOST_POSITIVE_FIXNUM
)
6067 else if (min
>= LONG_MIN
&& max
<= LONG_MAX
)
6069 if (mpz_fits_slong_p (SCM_I_BIG_MPZ (val
)))
6071 long n
= mpz_get_si (SCM_I_BIG_MPZ (val
));
6072 return n
>= min
&& n
<= max
;
6082 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6083 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6086 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6087 SCM_I_BIG_MPZ (val
));
6089 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) >= 0)
6101 return n
>= min
&& n
<= max
;
6109 scm_is_unsigned_integer (SCM val
, scm_t_uintmax min
, scm_t_uintmax max
)
6111 if (SCM_I_INUMP (val
))
6113 scm_t_signed_bits n
= SCM_I_INUM (val
);
6114 return n
>= 0 && ((scm_t_uintmax
)n
) >= min
&& ((scm_t_uintmax
)n
) <= max
;
6116 else if (SCM_BIGP (val
))
6118 if (max
<= SCM_MOST_POSITIVE_FIXNUM
)
6120 else if (max
<= ULONG_MAX
)
6122 if (mpz_fits_ulong_p (SCM_I_BIG_MPZ (val
)))
6124 unsigned long n
= mpz_get_ui (SCM_I_BIG_MPZ (val
));
6125 return n
>= min
&& n
<= max
;
6135 if (mpz_sgn (SCM_I_BIG_MPZ (val
)) < 0)
6138 if (mpz_sizeinbase (SCM_I_BIG_MPZ (val
), 2)
6139 > CHAR_BIT
*sizeof (scm_t_uintmax
))
6142 mpz_export (&n
, &count
, 1, sizeof (scm_t_uintmax
), 0, 0,
6143 SCM_I_BIG_MPZ (val
));
6145 return n
>= min
&& n
<= max
;
6153 scm_i_range_error (SCM bad_val
, SCM min
, SCM max
)
6155 scm_error (scm_out_of_range_key
,
6157 "Value out of range ~S to ~S: ~S",
6158 scm_list_3 (min
, max
, bad_val
),
6159 scm_list_1 (bad_val
));
6162 #define TYPE scm_t_intmax
6163 #define TYPE_MIN min
6164 #define TYPE_MAX max
6165 #define SIZEOF_TYPE 0
6166 #define SCM_TO_TYPE_PROTO(arg) scm_to_signed_integer (arg, scm_t_intmax min, scm_t_intmax max)
6167 #define SCM_FROM_TYPE_PROTO(arg) scm_from_signed_integer (arg)
6168 #include "libguile/conv-integer.i.c"
6170 #define TYPE scm_t_uintmax
6171 #define TYPE_MIN min
6172 #define TYPE_MAX max
6173 #define SIZEOF_TYPE 0
6174 #define SCM_TO_TYPE_PROTO(arg) scm_to_unsigned_integer (arg, scm_t_uintmax min, scm_t_uintmax max)
6175 #define SCM_FROM_TYPE_PROTO(arg) scm_from_unsigned_integer (arg)
6176 #include "libguile/conv-uinteger.i.c"
6178 #define TYPE scm_t_int8
6179 #define TYPE_MIN SCM_T_INT8_MIN
6180 #define TYPE_MAX SCM_T_INT8_MAX
6181 #define SIZEOF_TYPE 1
6182 #define SCM_TO_TYPE_PROTO(arg) scm_to_int8 (arg)
6183 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int8 (arg)
6184 #include "libguile/conv-integer.i.c"
6186 #define TYPE scm_t_uint8
6188 #define TYPE_MAX SCM_T_UINT8_MAX
6189 #define SIZEOF_TYPE 1
6190 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint8 (arg)
6191 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint8 (arg)
6192 #include "libguile/conv-uinteger.i.c"
6194 #define TYPE scm_t_int16
6195 #define TYPE_MIN SCM_T_INT16_MIN
6196 #define TYPE_MAX SCM_T_INT16_MAX
6197 #define SIZEOF_TYPE 2
6198 #define SCM_TO_TYPE_PROTO(arg) scm_to_int16 (arg)
6199 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int16 (arg)
6200 #include "libguile/conv-integer.i.c"
6202 #define TYPE scm_t_uint16
6204 #define TYPE_MAX SCM_T_UINT16_MAX
6205 #define SIZEOF_TYPE 2
6206 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint16 (arg)
6207 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint16 (arg)
6208 #include "libguile/conv-uinteger.i.c"
6210 #define TYPE scm_t_int32
6211 #define TYPE_MIN SCM_T_INT32_MIN
6212 #define TYPE_MAX SCM_T_INT32_MAX
6213 #define SIZEOF_TYPE 4
6214 #define SCM_TO_TYPE_PROTO(arg) scm_to_int32 (arg)
6215 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int32 (arg)
6216 #include "libguile/conv-integer.i.c"
6218 #define TYPE scm_t_uint32
6220 #define TYPE_MAX SCM_T_UINT32_MAX
6221 #define SIZEOF_TYPE 4
6222 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint32 (arg)
6223 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint32 (arg)
6224 #include "libguile/conv-uinteger.i.c"
6226 #define TYPE scm_t_wchar
6227 #define TYPE_MIN (scm_t_int32)-1
6228 #define TYPE_MAX (scm_t_int32)0x10ffff
6229 #define SIZEOF_TYPE 4
6230 #define SCM_TO_TYPE_PROTO(arg) scm_to_wchar (arg)
6231 #define SCM_FROM_TYPE_PROTO(arg) scm_from_wchar (arg)
6232 #include "libguile/conv-integer.i.c"
6234 #define TYPE scm_t_int64
6235 #define TYPE_MIN SCM_T_INT64_MIN
6236 #define TYPE_MAX SCM_T_INT64_MAX
6237 #define SIZEOF_TYPE 8
6238 #define SCM_TO_TYPE_PROTO(arg) scm_to_int64 (arg)
6239 #define SCM_FROM_TYPE_PROTO(arg) scm_from_int64 (arg)
6240 #include "libguile/conv-integer.i.c"
6242 #define TYPE scm_t_uint64
6244 #define TYPE_MAX SCM_T_UINT64_MAX
6245 #define SIZEOF_TYPE 8
6246 #define SCM_TO_TYPE_PROTO(arg) scm_to_uint64 (arg)
6247 #define SCM_FROM_TYPE_PROTO(arg) scm_from_uint64 (arg)
6248 #include "libguile/conv-uinteger.i.c"
6251 scm_to_mpz (SCM val
, mpz_t rop
)
6253 if (SCM_I_INUMP (val
))
6254 mpz_set_si (rop
, SCM_I_INUM (val
));
6255 else if (SCM_BIGP (val
))
6256 mpz_set (rop
, SCM_I_BIG_MPZ (val
));
6258 scm_wrong_type_arg_msg (NULL
, 0, val
, "exact integer");
6262 scm_from_mpz (mpz_t val
)
6264 return scm_i_mpz2num (val
);
6268 scm_is_real (SCM val
)
6270 return scm_is_true (scm_real_p (val
));
6274 scm_is_rational (SCM val
)
6276 return scm_is_true (scm_rational_p (val
));
6280 scm_to_double (SCM val
)
6282 if (SCM_I_INUMP (val
))
6283 return SCM_I_INUM (val
);
6284 else if (SCM_BIGP (val
))
6285 return scm_i_big2dbl (val
);
6286 else if (SCM_FRACTIONP (val
))
6287 return scm_i_fraction2double (val
);
6288 else if (SCM_REALP (val
))
6289 return SCM_REAL_VALUE (val
);
6291 scm_wrong_type_arg_msg (NULL
, 0, val
, "real number");
6295 scm_from_double (double val
)
6297 SCM z
= scm_double_cell (scm_tc16_real
, 0, 0, 0);
6298 SCM_REAL_VALUE (z
) = val
;
6302 #if SCM_ENABLE_DEPRECATED == 1
6305 scm_num2float (SCM num
, unsigned long int pos
, const char *s_caller
)
6307 scm_c_issue_deprecation_warning
6308 ("`scm_num2float' is deprecated. Use scm_to_double instead.");
6312 float res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6316 scm_out_of_range (NULL
, num
);
6319 return scm_to_double (num
);
6323 scm_num2double (SCM num
, unsigned long int pos
, const char *s_caller
)
6325 scm_c_issue_deprecation_warning
6326 ("`scm_num2double' is deprecated. Use scm_to_double instead.");
6330 double res
= mpz_get_d (SCM_I_BIG_MPZ (num
));
6334 scm_out_of_range (NULL
, num
);
6337 return scm_to_double (num
);
6343 scm_is_complex (SCM val
)
6345 return scm_is_true (scm_complex_p (val
));
6349 scm_c_real_part (SCM z
)
6351 if (SCM_COMPLEXP (z
))
6352 return SCM_COMPLEX_REAL (z
);
6355 /* Use the scm_real_part to get proper error checking and
6358 return scm_to_double (scm_real_part (z
));
6363 scm_c_imag_part (SCM z
)
6365 if (SCM_COMPLEXP (z
))
6366 return SCM_COMPLEX_IMAG (z
);
6369 /* Use the scm_imag_part to get proper error checking and
6370 dispatching. The result will almost always be 0.0, but not
6373 return scm_to_double (scm_imag_part (z
));
6378 scm_c_magnitude (SCM z
)
6380 return scm_to_double (scm_magnitude (z
));
6386 return scm_to_double (scm_angle (z
));
6390 scm_is_number (SCM z
)
6392 return scm_is_true (scm_number_p (z
));
6396 /* In the following functions we dispatch to the real-arg funcs like log()
6397 when we know the arg is real, instead of just handing everything to
6398 clog() for instance. This is in case clog() doesn't optimize for a
6399 real-only case, and because we have to test SCM_COMPLEXP anyway so may as
6400 well use it to go straight to the applicable C func. */
6402 SCM_DEFINE (scm_log
, "log", 1, 0, 0,
6404 "Return the natural logarithm of @var{z}.")
6405 #define FUNC_NAME s_scm_log
6407 if (SCM_COMPLEXP (z
))
6409 #if HAVE_COMPLEX_DOUBLE && HAVE_CLOG && defined (SCM_COMPLEX_VALUE)
6410 return scm_from_complex_double (clog (SCM_COMPLEX_VALUE (z
)));
6412 double re
= SCM_COMPLEX_REAL (z
);
6413 double im
= SCM_COMPLEX_IMAG (z
);
6414 return scm_c_make_rectangular (log (hypot (re
, im
)),
6420 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6421 although the value itself overflows. */
6422 double re
= scm_to_double (z
);
6423 double l
= log (fabs (re
));
6425 return scm_from_double (l
);
6427 return scm_c_make_rectangular (l
, M_PI
);
6433 SCM_DEFINE (scm_log10
, "log10", 1, 0, 0,
6435 "Return the base 10 logarithm of @var{z}.")
6436 #define FUNC_NAME s_scm_log10
6438 if (SCM_COMPLEXP (z
))
6440 /* Mingw has clog() but not clog10(). (Maybe it'd be worth using
6441 clog() and a multiply by M_LOG10E, rather than the fallback
6442 log10+hypot+atan2.) */
6443 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_CLOG10 \
6444 && defined SCM_COMPLEX_VALUE
6445 return scm_from_complex_double (clog10 (SCM_COMPLEX_VALUE (z
)));
6447 double re
= SCM_COMPLEX_REAL (z
);
6448 double im
= SCM_COMPLEX_IMAG (z
);
6449 return scm_c_make_rectangular (log10 (hypot (re
, im
)),
6450 M_LOG10E
* atan2 (im
, re
));
6455 /* ENHANCE-ME: When z is a bignum the logarithm will fit a double
6456 although the value itself overflows. */
6457 double re
= scm_to_double (z
);
6458 double l
= log10 (fabs (re
));
6460 return scm_from_double (l
);
6462 return scm_c_make_rectangular (l
, M_LOG10E
* M_PI
);
6468 SCM_DEFINE (scm_exp
, "exp", 1, 0, 0,
6470 "Return @math{e} to the power of @var{z}, where @math{e} is the\n"
6471 "base of natural logarithms (2.71828@dots{}).")
6472 #define FUNC_NAME s_scm_exp
6474 if (SCM_COMPLEXP (z
))
6476 #if HAVE_COMPLEX_DOUBLE && HAVE_CEXP && defined (SCM_COMPLEX_VALUE)
6477 return scm_from_complex_double (cexp (SCM_COMPLEX_VALUE (z
)));
6479 return scm_c_make_polar (exp (SCM_COMPLEX_REAL (z
)),
6480 SCM_COMPLEX_IMAG (z
));
6485 /* When z is a negative bignum the conversion to double overflows,
6486 giving -infinity, but that's ok, the exp is still 0.0. */
6487 return scm_from_double (exp (scm_to_double (z
)));
6493 SCM_DEFINE (scm_sqrt
, "sqrt", 1, 0, 0,
6495 "Return the square root of @var{z}. Of the two possible roots\n"
6496 "(positive and negative), the one with the a positive real part\n"
6497 "is returned, or if that's zero then a positive imaginary part.\n"
6501 "(sqrt 9.0) @result{} 3.0\n"
6502 "(sqrt -9.0) @result{} 0.0+3.0i\n"
6503 "(sqrt 1.0+1.0i) @result{} 1.09868411346781+0.455089860562227i\n"
6504 "(sqrt -1.0-1.0i) @result{} 0.455089860562227-1.09868411346781i\n"
6506 #define FUNC_NAME s_scm_sqrt
6508 if (SCM_COMPLEXP (x
))
6510 #if defined HAVE_COMPLEX_DOUBLE && defined HAVE_USABLE_CSQRT \
6511 && defined SCM_COMPLEX_VALUE
6512 return scm_from_complex_double (csqrt (SCM_COMPLEX_VALUE (x
)));
6514 double re
= SCM_COMPLEX_REAL (x
);
6515 double im
= SCM_COMPLEX_IMAG (x
);
6516 return scm_c_make_polar (sqrt (hypot (re
, im
)),
6517 0.5 * atan2 (im
, re
));
6522 double xx
= scm_to_double (x
);
6524 return scm_c_make_rectangular (0.0, sqrt (-xx
));
6526 return scm_from_double (sqrt (xx
));
6538 mpz_init_set_si (z_negative_one
, -1);
6540 /* It may be possible to tune the performance of some algorithms by using
6541 * the following constants to avoid the creation of bignums. Please, before
6542 * using these values, remember the two rules of program optimization:
6543 * 1st Rule: Don't do it. 2nd Rule (experts only): Don't do it yet. */
6544 scm_c_define ("most-positive-fixnum",
6545 SCM_I_MAKINUM (SCM_MOST_POSITIVE_FIXNUM
));
6546 scm_c_define ("most-negative-fixnum",
6547 SCM_I_MAKINUM (SCM_MOST_NEGATIVE_FIXNUM
));
6549 scm_add_feature ("complex");
6550 scm_add_feature ("inexact");
6551 flo0
= scm_from_double (0.0);
6553 /* determine floating point precision */
6554 for (i
=2; i
<= SCM_MAX_DBL_RADIX
; ++i
)
6556 init_dblprec(&scm_dblprec
[i
-2],i
);
6557 init_fx_radix(fx_per_radix
[i
-2],i
);
6560 /* hard code precision for base 10 if the preprocessor tells us to... */
6561 scm_dblprec
[10-2] = (DBL_DIG
> 20) ? 20 : DBL_DIG
;
6564 exactly_one_half
= scm_divide (SCM_I_MAKINUM (1), SCM_I_MAKINUM (2));
6565 #include "libguile/numbers.x"