conservatory
/
stk
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
stk/Src/number.c

1279 lines
32 KiB
C
Raw Permalink Blame History

/*
*
* n u m b e r . c -- Numbers management
*
*
* Numbers recognized by the interpreter are:
* - integer (which fit in a C long)
* - bignum (arbitrary precision integer)
* - flonum (represented as a C double)
*
* Bignum use the GNU gmp API. However to avoid to fall under the GPL terms
* you can use the FGMP package (FGMP is a public domain implementation of
* a subset of the GNU gmp library with the same API, written by Mark
* Henderson <markh@wimsey.bc.ca>). If your concern is speed, and if the
* GPL is not a problem for you, use the gmp package.
*
*
*****
*
* Copyright © 1993-1999 Erick Gallesio - I3S-CNRS/ESSI <eg@unice.fr>
*
*
* Permission to use, copy, modify, distribute,and license this
* software and its documentation for any purpose is hereby granted,
* provided that existing copyright notices are retained in all
* copies and that this notice is included verbatim in any
* distributions. No written agreement, license, or royalty fee is
* required for any of the authorized uses.
* This software is provided ``AS IS'' without express or implied
* warranty.
*
*
* Author: Erick Gallesio [eg@kaolin.unice.fr]
* Creation date: 12-May-1993 10:34
* Last file update: 3-Sep-1999 20:21 (eg)
*/
#include "stk.h"
#define SMALLNUMP(x) (fabs(x) <= SMALLINT_MAX)
#define ISINT(x) (floor(x) == (x))
/*
* This interface has been added for the Frank Schmuelling complex-number
* package (Frank email is * <frank@zeus.ph1.Uni-Koeln.DE>). Don't rely on
* it, since rationnal and complex numbers will be present in a future
* release. As soon as the full tower will be implemented in STk,
* this door will be closed.
*/
#define do_addition _STk_do_addition
#define do_substract _STk_do_substract
#define do_multiply _STk_do_multiply
#define do_divide _STk_do_divide
#define bignum2double _STk_bignum2double
void _STk_do_multiply(SCM *x, SCM y);
void _STk_do_addition(SCM *x, SCM y);
static double do_compare(SCM x, SCM y);
static SCM makesmallint(long x)
{
SCM z;
#ifndef COMPACT_SMALL_CST
NEWCELL(z, tc_integer);
#endif
SET_INTEGER(z, x);
return z;
}
static SCM makebignum(char *s)
{
SCM z;
NEWCELL(z, tc_bignum);
BIGNUM(z) = must_malloc(sizeof(MP_INT));
mpz_init_set_str(BIGNUM(z), s, 10);
return z;
}
/******************************************************************************
*
* Conversion Functions
*
******************************************************************************/
static int digitp(char c, long base)
{
c = ('0' <= c && c <= '9') ? c - '0':
('a' <= c && c <= 'f') ? c - 'a' + 10:
('A' <= c && c <= 'F') ? c - 'A' + 10:
(c == '#') ? 0 :
100;
return (c < base);
}
static SCM int2bignum(long n)
{
SCM z;
NEWCELL(z, tc_bignum);
BIGNUM(z) = must_malloc(sizeof(MP_INT));
mpz_init_set_si(BIGNUM(z), n);
return z;
}
static SCM uint2bignum(unsigned long n)
{
SCM z;
NEWCELL(z, tc_bignum);
BIGNUM(z) = must_malloc(sizeof(MP_INT));
mpz_init_set_ui(BIGNUM(z), n);
return z;
}
static double bignum2double(MP_INT *bn)
{
char *s= mpz_get_str(NULL, 10, bn);
double d = (double) atof(s);
free(s);
return d;
}
static SCM double2integer(double n) /* small or big depending of n's size */
{
int i, j;
size_t size = 20;
char *tmp = NULL;
SCM z;
if (!ISINT(n)) Err("cannot convert this number to an integer", STk_makenumber(n));
/* Try first to convert n to a long */
if (((double) SMALLINT_MIN <= n) && (n <= (double) SMALLINT_MAX))
return STk_makeinteger((long) n);
/* n doesn't fit in a long => build a bignum. THIS IS VERY INEFFICIENT */
tmp = must_malloc(size);
i = 0;
if (n < 0.0) { tmp[i++] = '-'; n = -n; }
do {
if (i >= size) tmp = must_realloc(tmp, size *= 2);
tmp[i++] = (int) fmod(n, (double) 10) + '0';
n = floor(n / 10.0);
}
while (n > 0.0);
tmp[i] = 0;
/* Reverse the content of string tmp */
for (i=i-1, j=(tmp[0]=='-'); i > j; i--, j++) {
char c = tmp[i];
tmp[i] = tmp[j];
tmp[j] = c;
}
/* tmp contains a textual representation of n. Convert it to a bignum */
z = makebignum(tmp);
if (tmp) free(tmp);
return z;
}
/* Convert a number to a C-string. Result must be freed if != from buffer */
char *STk_number2Cstr(SCM n, long base, char buffer[])
{
char *s = buffer;
switch (TYPE(n)) {
case tc_flonum:
if (base != 10) Err("base must be 10 for this number", n);
sprintf(buffer, "%.15g", FLONM(n));
if (strchr(buffer, '.') == NULL && strchr(buffer, 'e') == NULL)
strcat(buffer, ".0");
return buffer;
case tc_integer:
{
long tmp, val = INTEGER(n);
int u;
if (val < 0) {
val = -val;
*s++ = '-';
}
for (s++, tmp=val; tmp >= base; tmp /= base) s++;
*s = '\000'; tmp = val;
do {
u = tmp % base;
*(--s) = u + ((u < 10) ? '0' : 'a'-10);
tmp /= base;
}
while (tmp);
return buffer;
}
case tc_bignum:
s = must_malloc(mpz_sizeinbase(BIGNUM(n), base) + 2);
s = mpz_get_str(s, base, BIGNUM(n));
return s;
}
return NULL; /* never reached */
}
/* Return #f if the C string doesn't contain a valid number. */
/* UGLY! must be rewritten */
SCM STk_Cstr2number(char *str, long base)
{
int i, adigit=0, isint=1, exact=' ', radix = 0;
register char *p = str;
for (i = 0; i < 2; i++) {
if (*p == '#') {
p += 1;
switch (*p++) {
case 'e': if (exact == ' ') { exact = 'e'; break; } else return Ntruth;
case 'i': if (exact == ' ') { exact = 'i'; break; } else return Ntruth;
case 'b': if (!radix) {base = 2; radix = 1; break;} else return Ntruth;
case 'o': if (!radix) {base = 8; radix = 1; break;} else return Ntruth;
case 'd': if (!radix) {base = 10; radix = 1; break;} else return Ntruth;
case 'x': if (!radix) {base = 16; radix = 1; break;} else return Ntruth;
default: return Ntruth;
}
str += 2;
}
if (*p != '#') break;
}
if (*p == '-' || *p == '+') p+=1;
if (*p == '#') goto End;
while(digitp(*p, base)) { p+=1; adigit=1; if (*p == '#') isint = 0; }
if (*p=='.'){
isint = 0; p += 1;
while(digitp(*p, base)) { p+=1; adigit=1; }
}
if (!adigit) goto End;
if (*p && strchr("eEsSfFdDlL", *p)) {
isint = 0;
p += 1;
if (*p == '-' || *p == '+') p+=1;
if (!digitp(*p, base)) goto End;
p+=1;
while (digitp(*p, base)) p+=1;
}
if (*p) goto End;
/* Now we are sure it is a number. Find the more adequate type */
if (isint) {
MP_INT n;
if (*str == '+') str+=1; /* mpz_init_set_str doesn't recognize +xyz !!! */
if (mpz_init_set_str(&n, str, base) < 0) {
mpz_clear(&n);
return Ntruth;
}
if (mpz_cmp_si(&n, SMALLINT_MIN) >=0 && mpz_cmp_si(&n, SMALLINT_MAX) <= 0) {
long num = mpz_get_si(&n);
mpz_clear(&n);
return (exact == 'i') ? STk_makenumber((double) num) : makesmallint(num);
}
/* It's a bignum */
if (exact == 'i') return STk_makenumber(bignum2double(&n));
else {
SCM z;
NEWCELL(z, tc_bignum);
BIGNUM(z) = must_malloc(sizeof(MP_INT));
mpz_init_set(BIGNUM(z), &n);
mpz_clear(&n);
return z;
}
}
/* It's a float */
if (exact == 'e') return Ntruth; /* e cannot be specified on e float */
if (base == 10) {
/* Replace sharp signs by 0 */
for(p=str; *p; p++)
switch (*p) {
case '#': *p = '0'; break;
case 's': case 'S': case 'f': case 'F':
case 'd': case 'D': case 'l': case 'L': *p = 'e';
}
return STk_makenumber((double) atof(str));
}
End:
return Ntruth;
}
long STk_integer_value(SCM x) /* Returns LONG_MIN if not representable as int */
{
if (INTEGERP(x)) return INTEGER(x);
if (BIGNUMP(x)) {
if (mpz_cmp_si(BIGNUM(x), SMALLINT_MIN)>=0 &&
mpz_cmp_si(BIGNUM(x), SMALLINT_MAX)<=0)
return mpz_get_si(BIGNUM(x));
}
return LONG_MIN;
}
long STk_integer_value_no_overflow(SCM x) /* Returns LONG_MIN if not an integer */
{
if (INTEGERP(x)) return INTEGER(x);
if (BIGNUMP(x)) return mpz_get_si(BIGNUM(x));
return LONG_MIN;
}
int STk_equal_numbers(SCM number1, SCM number2) /* number1 = number2 */
{
return do_compare(number1, number2) == 0;
}
long STk_integer2long(SCM x)
{
if (BIGNUMP(x)) return mpz_get_si(BIGNUM(x));
return INTEGER(x);
}
unsigned long STk_integer2ulong(SCM x)
{
if (BIGNUMP(x)) return mpz_get_ui(BIGNUM(x));
return (unsigned long) INTEGER(x);
}
/******************************************************************************/
SCM STk_makenumber(double x)
{
/* Floats are not stored in a struct obj since this leads to memory consumption
* This memory consumption due to alignment problems.
* For instance on a Sun 4, where double are 8 bytes, a struct obj with a double
* in line will occupy 16 bytes whereas it occupies only 12 bytes if the double
* is mallocated.
* This change (94/08/29) will give worst performances when crunching numbers,
* but use Fortran if this is your job :->
*/
SCM z;
NEWCELL(z,tc_flonum);
FLONM(z) = x;
return z;
}
SCM STk_makeinteger(long x)
{
return (SMALLINT_MIN <= x && x <= SMALLINT_MAX) ? makesmallint(x): int2bignum(x);
}
SCM STk_makeunsigned(unsigned long x)
{
return (x <= SMALLINT_MAX) ? makesmallint(x): uint2bignum(x);
}
/******************************************************************************
*
* Simple operations
*
******************************************************************************/
static SCM clone(SCM number)
{
/* clone a number */
switch (TYPE(number)) {
case tc_integer: return makesmallint(INTEGER(number));
case tc_bignum: {
SCM z;
NEWCELL(z, tc_bignum);
BIGNUM(z) = must_malloc(sizeof(MP_INT));
mpz_init(BIGNUM(z));
mpz_set(BIGNUM(z), BIGNUM(number));
return z;
}
case tc_flonum: return STk_makenumber(FLONM(number));
default: /* This is not a number. Return the parameter unmodified. */
return number; /* Error will be signaled later */
}
}
void _STk_do_multiply(SCM *x, SCM y)
{
switch (TYPE(*x)) {
case tc_integer:
switch (TYPE(y)) {
case tc_integer:
{
double prod;
if (SMALLNUMP(prod=(double) INTEGER(*x) * INTEGER(y)))
SET_INTEGER(*x, (long) prod);
else {
MP_INT tmp;
*x = int2bignum(INTEGER(*x));
mpz_init_set_si(&tmp, INTEGER(y));
mpz_mul(BIGNUM(*x), BIGNUM(*x), &tmp);
mpz_clear(&tmp);
}
}
break;
case tc_bignum:
*x = int2bignum(INTEGER(*x));
mpz_mul(BIGNUM(*x), BIGNUM(*x), BIGNUM(y));
break;
case tc_flonum:
*x = STk_makenumber(INTEGER(*x) * FLONM(y));
break;
}
break;
case tc_bignum:
switch (TYPE(y)) {
case tc_integer:
{
SCM tmp = int2bignum(INTEGER(y));
mpz_mul(BIGNUM(*x), BIGNUM(*x), BIGNUM(tmp));
}
break;
case tc_bignum:
mpz_mul(BIGNUM(*x), BIGNUM(*x), BIGNUM(y));
break;
case tc_flonum:
*x = STk_makenumber(bignum2double(BIGNUM(*x)) * FLONM(y));
break;
}
break;
case tc_flonum:
switch (TYPE(y)) {
case tc_integer:
FLONM(*x) *= INTEGER(y);
break;
case tc_bignum:
FLONM(*x) *= bignum2double(BIGNUM(y));
break;
case tc_flonum:
FLONM(*x) *= FLONM(y);
break;
}
break;
}
}
void _STk_do_addition(SCM *x, SCM y)
{
switch (TYPE(*x)) {
case tc_integer:
switch (TYPE(y)) {
case tc_integer:
{
double add =(double) INTEGER(*x) + INTEGER(y);
if (SMALLNUMP(add))
SET_INTEGER(*x, (long) add);
else {
*x = int2bignum(INTEGER(*x));
mpz_add(BIGNUM(*x), BIGNUM(*x), BIGNUM(int2bignum(INTEGER(y))));
}
}
break;
case tc_bignum:
*x = int2bignum(INTEGER(*x));
mpz_add(BIGNUM(*x), BIGNUM(*x), BIGNUM(y));
break;
case tc_flonum:
*x = STk_makenumber(INTEGER(*x) + FLONM(y));
break;
}
break;
case tc_bignum:
switch (TYPE(y)) {
case tc_integer:
{
SCM tmp = int2bignum(INTEGER(y));
mpz_add(BIGNUM(*x), BIGNUM(*x), BIGNUM(tmp));
}
break;
case tc_bignum:
mpz_add(BIGNUM(*x), BIGNUM(*x), BIGNUM(y));
break;
case tc_flonum:
*x = STk_makenumber(bignum2double(BIGNUM(*x)) + FLONM(y));
break;
}
break;
case tc_flonum:
switch (TYPE(y)) {
case tc_integer:
FLONM(*x) += INTEGER(y);
break;
case tc_bignum:
FLONM(*x) += bignum2double(BIGNUM(y));
break;
case tc_flonum:
FLONM(*x) += FLONM(y);
break;
}
break;
}
}
static void _STk_do_substract(SCM *x, SCM y)
{
switch (TYPE(*x)) {
case tc_integer:
switch (TYPE(y)) {
case tc_integer:
{
double add = (double) INTEGER(*x) - INTEGER(y);
if (SMALLNUMP(add))
SET_INTEGER(*x,(long) add);
else {
*x = int2bignum(INTEGER(*x));
mpz_sub(BIGNUM(*x), BIGNUM(*x), BIGNUM(int2bignum(INTEGER(y))));
}
}
break;
case tc_bignum:
*x = int2bignum(INTEGER(*x));
mpz_sub(BIGNUM(*x), BIGNUM(*x), BIGNUM(y));
break;
case tc_flonum:
*x = STk_makenumber(INTEGER(*x) - FLONM(y));
break;
}
break;
case tc_bignum:
switch (TYPE(y)) {
case tc_integer:
{
SCM tmp = int2bignum(INTEGER(y));
mpz_sub(BIGNUM(*x), BIGNUM(*x), BIGNUM(tmp));
}
break;
case tc_bignum:
mpz_sub(BIGNUM(*x), BIGNUM(*x), BIGNUM(y));
break;
case tc_flonum:
*x = STk_makenumber(bignum2double(BIGNUM(*x)) - FLONM(y));
break;
}
break;
case tc_flonum:
switch (TYPE(y)) {
case tc_integer:
FLONM(*x) -= INTEGER(y);
break;
case tc_bignum:
FLONM(*x) -= bignum2double(BIGNUM(y));
break;
case tc_flonum:
FLONM(*x) -= FLONM(y);
break;
}
break;
}
}
static void _STk_do_divide(SCM *x, SCM y)
{
switch (TYPE(*x)) {
case tc_integer:
switch (TYPE(y)) {
case tc_integer:
{
double div = (double) INTEGER(*x) / INTEGER(y);
if (ISINT(div))
SET_INTEGER(*x, (double) div);
else
*x = STk_makenumber(div);
}
break;
case tc_bignum:
{
MP_INT q, r;
mpz_init(&q); mpz_init(&r);
*x = int2bignum(INTEGER(*x));
mpz_mdivmod(&q, &r, BIGNUM(*x), BIGNUM(y));
if (mpz_cmp_ui(&r, 0L) == 0)
mpz_set(BIGNUM(*x), &q);
else
*x=STk_makenumber(bignum2double(BIGNUM(*x)) /
bignum2double(BIGNUM(y)));
mpz_clear(&q); mpz_clear(&r);
}
break;
case tc_flonum:
*x = STk_makenumber(INTEGER(*x) / FLONM(y));
break;
}
break;
case tc_bignum:
switch (TYPE(y)) {
case tc_integer:
case tc_bignum:
{
MP_INT q, r;
mpz_init(&q); mpz_init(&r);
if (INTEGERP(y))
y = int2bignum(INTEGER(y));
mpz_mdivmod(&q, &r, BIGNUM(*x), BIGNUM(y));
if (mpz_cmp_ui(&r, 0L) == 0)
mpz_set(BIGNUM(*x), &q);
else
*x=STk_makenumber(bignum2double(BIGNUM(*x)) /
bignum2double(BIGNUM(y)));
mpz_clear(&q); mpz_clear(&r);
}
break;
case tc_flonum:
*x = STk_makenumber(bignum2double(BIGNUM(*x)) / FLONM(y));
break;
}
break;
case tc_flonum:
switch (TYPE(y)) {
case tc_integer:
FLONM(*x) /= INTEGER(y);
break;
case tc_bignum:
FLONM(*x) /= bignum2double(BIGNUM(y));
break;
case tc_flonum:
FLONM(*x) /= FLONM(y);
break;
}
break;
}
}
static SCM do_integer_division(SCM x, SCM y, int quotientp)
{
MP_INT q, r;
SCM res;
int exact = 1;
/* Divide x by y and return its quotient or remainder. y is not 0 */
if (FLONUMP(x))
if (ISINT(FLONM(x))) { x = double2integer(FLONM(x)); exact = 0; }
else Err("bad number in an integer division", x);
if (FLONUMP(y))
if (ISINT(FLONM(y))) { y = double2integer(FLONM(y)); exact = 0; }
else Err("bad number in an integer division", y);
/* Here, x and y can only be integer or bignum (not float) */
if (INTEGERP(x))
if (INTEGERP(y)) {
res = makesmallint(quotientp? INTEGER(x)/INTEGER(y):
INTEGER(x)%INTEGER(y));
return exact? res : STk_exact2inexact(res);
}
else
x = int2bignum(INTEGER(x));
else {
/* x is a bignum */
if (INTEGERP(y))
y = int2bignum(INTEGER(y));
}
/* x and y are both bignums */
res = makebignum("0");
mpz_init(&q); mpz_init(&r);
mpz_divmod(&q, &r, BIGNUM(x), BIGNUM(y)); /* use divmod instead of mdivmod */
mpz_set(BIGNUM(res), quotientp? &q : &r); /* Thanks to kerch@parc.xerox.com*/
mpz_clear(&q); mpz_clear(&r);
return exact? res : STk_exact2inexact(res);
}
static double do_compare(SCM x, SCM y)
{
switch (TYPE(x)) {
case tc_integer:
switch (TYPE(y)) {
case tc_integer: return (double) (INTEGER(x)-INTEGER(y));
case tc_bignum: return (double) -(mpz_cmp_si(BIGNUM(y), INTEGER(x)));
case tc_flonum: return (double) INTEGER(x) - FLONM(y);
}
case tc_bignum:
switch (TYPE(y)) {
case tc_integer: return (double) mpz_cmp_si(BIGNUM(x), INTEGER(y));
case tc_bignum: return (double) mpz_cmp(BIGNUM(x), BIGNUM(y));
case tc_flonum: return bignum2double(BIGNUM(x))-FLONM(y);
}
case tc_flonum:
switch (TYPE(y)) {
case tc_integer: return (FLONM(x) - (double) INTEGER(y));
case tc_bignum: return (FLONM(x) - bignum2double(BIGNUM(y)));
case tc_flonum: return (FLONM(x) - FLONM(y));
}
default: return 0.0; /* never reached */
}
}
/******************************************************************************
*
* Scheme primitives and utilities
*
******************************************************************************/
/**** Section 6.5 ****/
PRIMITIVE STk_numberp(SCM x)
{
return NUMBERP(x)? Truth : Ntruth;
}
PRIMITIVE STk_integerp(SCM x)
{
switch (TYPE(x)) {
case tc_integer:
case tc_bignum: return Truth;
case tc_flonum: {
double val = FLONM(x);
return (floor(val) == val)? Truth : Ntruth;
}
default: return Ntruth;
}
}
PRIMITIVE STk_exactp(SCM x)
{
if (NNUMBERP(x)) Err("exact?: Bad number", x);
return EXACTP(x) ? Truth: Ntruth;
}
PRIMITIVE STk_inexactp(SCM x)
{
if (NNUMBERP(x)) Err("inexact?: Bad number", x);
return NEXACTP(x) ? Truth: Ntruth;
}
#define Compare(name, operator) \
PRIMITIVE name(SCM l, SCM env, int from_eval) \
{ \
register SCM tmp1, tmp2; \
\
if (NCONSP(l)) Err("too few parameters", l); \
tmp1 = from_eval? EVALCAR(l): CAR(l); \
if (NNUMBERP(tmp1)) goto Error; \
\
for (l=CDR(l); NNULLP(l); l=CDR(l),tmp1=tmp2) { \
tmp2 = from_eval? EVALCAR(l): CAR(l); \
if (NNUMBERP(tmp2)) goto Error; \
if (do_compare(tmp1, tmp2) operator 0) return Ntruth; \
} \
return Truth; \
Error: \
Err("Bad number in a comparison", l); \
return UNDEFINED; /* never reached */ \
}
Compare(STk_numequal, !=)
Compare(STk_lessp, >=)
Compare(STk_greaterp, <=)
Compare(STk_lessep, >)
Compare(STk_greaterep, <)
PRIMITIVE STk_zerop(SCM n)
{
switch (TYPE(n)) {
case tc_integer: return (INTEGER(n) == 0) ? Truth: Ntruth;
case tc_flonum: return (FLONM(n) == 0.0) ? Truth: Ntruth;
case tc_bignum: return (mpz_cmp_ui(BIGNUM(n), 0L) == 0) ? Truth: Ntruth;
default: Err("zero?: bad number", n);
return UNDEFINED; /* never reached */
}
}
PRIMITIVE STk_positivep(SCM n)
{
switch (TYPE(n)) {
case tc_integer: return (INTEGER(n) > 0) ? Truth: Ntruth;
case tc_flonum: return (FLONM(n) > 0.0) ? Truth: Ntruth;
case tc_bignum: return (mpz_cmp_ui(BIGNUM(n), 0L) > 0) ? Truth: Ntruth;
default: Err("positive?: bad number", n);
return UNDEFINED; /* never reached */
}
}
PRIMITIVE STk_negativep(SCM n)
{
switch (TYPE(n)) {
case tc_integer: return (INTEGER(n) < 0) ? Truth: Ntruth;
case tc_flonum: return (FLONM(n) < 0.0) ? Truth: Ntruth;
case tc_bignum: return (mpz_cmp_ui(BIGNUM(n), 0L) < 0) ? Truth: Ntruth;
default: Err("negative?: bad number", n);
return UNDEFINED; /* never reached */
}
}
PRIMITIVE STk_oddp(SCM n)
{
MP_INT q, r;
long res;
switch (TYPE(n)) {
case tc_integer: return (INTEGER(n)%2) ? Truth: Ntruth;
case tc_bignum: mpz_init(&q), mpz_init(&r);
mpz_divmod_ui(&q, &r, BIGNUM(n), 2L);
res = mpz_cmp_ui(&r, 0L);
mpz_clear(&q); mpz_clear(&r);
return (res != 0) ? Truth: Ntruth;
default: Err("odd?: bad number", n);
return UNDEFINED; /* never reached */
}
}
PRIMITIVE STk_evenp(SCM n)
{
MP_INT q, r;
long res;
switch (TYPE(n)) {
case tc_integer: return (INTEGER(n)%2) ? Ntruth: Truth;
case tc_bignum: mpz_init(&q), mpz_init(&r);
mpz_divmod_ui(&q, &r, BIGNUM(n), 2L);
res = mpz_cmp_ui(&r, 0L);
mpz_clear(&q); mpz_clear(&r);
return (res == 0) ? Truth: Ntruth;
default: Err("even?: bad number", n);
return UNDEFINED; /* never reached */
}
}
PRIMITIVE STk_max(SCM l, SCM env, int from_eval)
{
register SCM tmp, max;
int inexact;
if (NULLP(l)) Err("max: bad number of arguments", NIL);
tmp = from_eval? EVALCAR(l): CAR(l);
if (NNUMBERP(tmp)) goto Error;
inexact = FLONUMP(tmp);
max = tmp;
for (l=CDR(l); NNULLP(l); l=CDR(l)) {
tmp = from_eval? EVALCAR(l): CAR(l);
if (NNUMBERP(tmp)) goto Error;
inexact |= FLONUMP(tmp);
if (do_compare(tmp, max) > 0) max = tmp;
}
return (inexact && EXACTP(max)) ? STk_exact2inexact(max) : max;
Error:
Err("max: bad number", tmp);
return UNDEFINED; /* never reached */
}
PRIMITIVE STk_min(SCM l, SCM env, int from_eval)
{
register SCM tmp, min;
int inexact;
if (NULLP(l)) Err("min: bad number of arguments", NIL);
tmp = from_eval? EVALCAR(l): CAR(l);
if (NNUMBERP(tmp)) goto Error;
inexact = FLONUMP(tmp);
min = tmp;
for (l=CDR(l); NNULLP(l); l=CDR(l)) {
tmp = from_eval? EVALCAR(l): CAR(l);
if (NNUMBERP(tmp)) goto Error;
inexact |= FLONUMP(tmp);
if (do_compare(tmp, min) < 0) min = tmp;
}
return (inexact && EXACTP(min)) ? STk_exact2inexact(min) : min;
Error:
Err("min: bad number", tmp);
return UNDEFINED; /* never reached */
}
PRIMITIVE STk_plus(SCM l, SCM env, int from_eval)
{
SCM tmp, res = makesmallint(0);
for ( ; CONSP(l); l=CDR(l)) {
tmp = from_eval? EVALCAR(l): CAR(l);
if (NNUMBERP(tmp)) Err("+: not a number", tmp);
do_addition(&res, tmp);
}
return res;
}
PRIMITIVE STk_difference(SCM l, SCM env, int from_eval)
{
SCM tmp, res;
if (NULLP(l)) Err("-: no argument given", NIL);
if (NULLP(CDR(l))) l = Cons(makesmallint(0), l); /* (- x) --> (- 0 x) */
tmp = res = clone(EVALCAR(l));
if (NNUMBERP(res)) goto Error;
for (l=CDR(l) ; CONSP(l); l=CDR(l)) {
tmp = from_eval? EVALCAR(l): CAR(l);
if (NNUMBERP(tmp)) goto Error;
do_substract(&res, tmp);
}
return res;
Error:
Err("-: not a number", tmp);
return UNDEFINED; /* never reached */
}
PRIMITIVE STk_times(SCM l, SCM env, int from_eval)
{
SCM tmp, res = makesmallint(1);
for ( ; CONSP(l); l=CDR(l)) {
tmp = from_eval? EVALCAR(l): CAR(l);
if (NNUMBERP(tmp)) Err("*: not a number", tmp);
do_multiply(&res, tmp);
}
return res;
}
PRIMITIVE STk_division(SCM l, SCM env, int from_eval)
{
SCM tmp, res;
if (NULLP(l)) Err("/: no argumenent given", NIL);
if (NULLP(CDR(l))) l = Cons(makesmallint(1), l); /* (/ x) --> (/ 1 x) */
tmp = res = clone(EVALCAR(l));
if (NNUMBERP(res)) goto Error;
for (l=CDR(l) ; CONSP(l); l=CDR(l)) {
tmp = from_eval? EVALCAR(l): CAR(l);
if (NNUMBERP(tmp) || STk_zerop(tmp) == Truth) goto Error;
do_divide(&res, tmp);
}
return res;
Error:
Err("/: not a valid number", tmp);
return UNDEFINED; /* never reached */
}
PRIMITIVE STk_absolute(SCM x)
{
switch (TYPE(x)) {
case tc_integer: return (INTEGER(x) < 0) ? makesmallint(-INTEGER(x)) : x;
case tc_flonum: return (FLONM(x) < 0.0) ? STk_makenumber(-FLONM(x)) : x;
case tc_bignum: if (mpz_cmp_ui(BIGNUM(x), 0L) < 0) {
SCM tmp = clone(x);
mpz_neg(BIGNUM(tmp), BIGNUM(x));
return tmp;
}
return x;
default: Err("abs: bad number", x);
return UNDEFINED; /* never reached */
}
}
PRIMITIVE STk_quotient(SCM n1, SCM n2)
{
if (NNUMBERP(n1) || NNUMBERP(n2))
Err("quotient: bad arguments", Cons(n1, n2));
if (STk_zerop(n2) == Truth) Err("quotient: division by 0", NIL);
return do_integer_division(n1, n2, TRUE);
}
PRIMITIVE STk_remainder(SCM n1, SCM n2)
{
if (NNUMBERP(n1) || NNUMBERP(n2))
Err("remainder: bad arguments", Cons(n1, n2));
if (STk_zerop(n2) == Truth) Err("remainder: division by 0", NIL);
return do_integer_division(n1, n2, FALSE);
}
PRIMITIVE STk_modulo(SCM n1, SCM n2)
{
SCM z;
if (NNUMBERP(n1) || NNUMBERP(n2))
Err("modulo: bad arguments", Cons(n1, n2));
if (STk_zerop(n2) == Truth) Err("modulo: division by 0", NIL);
z = do_integer_division(n1, n2, FALSE);
if (STk_negativep(n1) != STk_negativep(n2) && STk_zerop(z) != Truth)
/*kerch@parc.xerox.com*/
do_addition(&z, n2);
return z;
}
static SCM gcd2(SCM n1, SCM n2)
{
return (STk_zerop(n2) == Truth) ? n1 : gcd2(n2, STk_modulo(n1, n2));
}
PRIMITIVE STk_gcd(SCM l, SCM env, int from_eval)
{
register SCM tmp, res;
if (NULLP(l)) return makesmallint(0L);
tmp = EVALCAR(l);
if (NNUMBERP(tmp)) goto Error;
res = STk_absolute(tmp);
for (l=CDR(l); NNULLP(l); l=CDR(l)) {
tmp = from_eval? EVALCAR(l): CAR(l);
if (NNUMBERP(tmp)) goto Error;
res = STk_absolute(gcd2(res, tmp));
}
return res;
Error:
Err("gcd: bad number", tmp);
return UNDEFINED; /* never reached */
}
PRIMITIVE STk_lcm(SCM l, SCM env, int from_eval)
{
SCM tmp, res, gcd;
if (NULLP(l)) return makesmallint(1L);
tmp = EVALCAR(l);
if (NNUMBERP(tmp)) goto Error;
res = STk_absolute(clone(tmp));
for (l=CDR(l); NNULLP(l); l=CDR(l)) {
tmp=EVALCAR(l);
if (NNUMBERP(tmp)) goto Error;
if (STk_zerop(tmp)==Truth) return tmp;
gcd = gcd2(res, tmp);
do_multiply(&res, tmp);
do_divide(&res, gcd);
res = STk_absolute(res);
}
return res;
Error:
Err("lcm: bad number", tmp);
return UNDEFINED; /* never reached */
}
PRIMITIVE STk_floor(SCM x)
{
if (NNUMBERP(x)) Err("floor: bad number", x);
if (FLONUMP(x)) return STk_makenumber(floor(FLONM(x)));
return x;
}
PRIMITIVE STk_ceiling(SCM x)
{
if (NNUMBERP(x)) Err("ceiling: bad number", x);
if (FLONUMP(x)) return STk_makenumber(ceil(FLONM(x)));
return x;
}
PRIMITIVE STk_truncate(SCM x)
{
if (NNUMBERP(x)) Err("truncate: bad number", x);
if (FLONUMP(x)) {
double d = FLONM(x);
return STk_makenumber(d < 0.0 ? ceil(d) : floor(d));
}
return x;
}
PRIMITIVE STk_round(SCM x)
{
if (NNUMBERP(x)) Err("round: bad number", x);
if (FLONUMP(x)) {
/*
* R4RS states that round must returns the closest integer to x, rounding to
* even when x is halfway between two integers (round to even is for consistency
* with IEEE standard).
* This explains the (much too) complicate computation below
*/
double v = FLONM(x);
double v_plus_0_5 = v + 0.5;
double res = floor(v_plus_0_5);
return STk_makenumber(
(v_plus_0_5 == res && v_plus_0_5/2 != floor(v_plus_0_5/2)) ? res-1: res);
}
return x;
}
#if !defined(WIN32) && !defined(__STDC__)
# define concat(a,b) a/**/b
#else
# define concat(a,b) a##b
#endif
#define transcendental(fct) \
PRIMITIVE concat(STk_,fct)(SCM z) \
{ \
if (NUMBERP(z)) \
return STk_makenumber(fct(FLONM(STk_exact2inexact(z)))); \
Err("transcendental function: bad number", z); \
return UNDEFINED /* never reached */; \
}
transcendental(exp)
transcendental(log)
transcendental(sin)
transcendental(cos)
transcendental(tan)
transcendental(asin)
transcendental(acos)
PRIMITIVE STk_atan(SCM y, SCM x)
{
if (NNUMBERP(y)) Err("atan: bad number", y);
if (x == UNBOUND)
return STk_makenumber(atan(FLONM(STk_exact2inexact(y))));
else {
double fl_x = FLONM(STk_exact2inexact(x));
double fl_y = FLONM(STk_exact2inexact(y));
/* Make a special case for (atan 0 0) -- in this case return 0 */
return STk_makenumber((fl_x == 0 && fl_y == 0) ? 0L : atan2(fl_y, fl_x));
}
}
PRIMITIVE STk_sqrt(SCM z)
{
switch (TYPE(z)) {
case tc_integer: {
double d;
if (INTEGER(z) < 0) goto Error;
d = (double) sqrt((double) INTEGER(z));
return ISINT(d) ? makesmallint((long) d): STk_makenumber(d);
}
case tc_bignum: {
MP_INT root, remainder;
SCM res;
if (mpz_cmp_si(BIGNUM(z), 0L) < 0) goto Error;
mpz_init(&root); mpz_init(&remainder);
mpz_sqrtrem(&root, &remainder, BIGNUM(z));
if (mpz_cmp_si(&remainder, 0L) == 0) {
/* Result is an integer */
res = makebignum("0");
mpz_set(BIGNUM(res), &root);
}
else {
/* Result is a flonum */
res = STk_makenumber(bignum2double(&root));
}
mpz_clear(&root); mpz_clear(&remainder);
return res;
}
case tc_flonum: {
if (FLONM(z) < 0.0) goto Error;
return STk_makenumber((double) sqrt(FLONM(z)));
}
default:
Err("sqrt: bad number", z);
}
Error:
Err("sqrt: number is negative", z);
return UNDEFINED; /* never reached */
}
/*
(define (exact_expt z1 z2)
(define (square x) (* x x))
(cond ((= z2 0) 1)
((negative? z2) (/ 1 (exact_expt z1 (abs z2))))
((even? z2) (square (exact_expt z1 (/ z2 2))))
(else (* z1 (exact_expt z1 (- z2 1))))))
*/
static SCM exact_expt(SCM z1, SCM z2)
{
SCM res;
if (STk_zerop(z2) == Truth) return makesmallint(1);
if (STk_negativep(z2) == Truth) {
/* (/ 1 (expt z1 (abs z2) */
res = makesmallint(1);
do_divide(&res, exact_expt(z1, STk_absolute(z2)));
return res;
}
if (STk_evenp(z2) == Truth) {
/* (square (expt z1 (/ z2 2)) */
res = clone(z2);
do_divide(&res, makesmallint(2));
res = exact_expt(z1, res);
do_multiply(&res, res);
return res;
}
/* (* (expt z1 (- z2 1)) z1) */
res = clone(z2);
do_substract(&res, makesmallint(1));
res = exact_expt(z1, res);
do_multiply(&res, z1);
return res;
}
PRIMITIVE STk_expt(SCM z1, SCM z2)
{
if (NNUMBERP(z1)) Err("expt: bad number", z1);
if (NNUMBERP(z2)) Err("expt: bad number", z2);
return (EXACTP(z1) && EXACTP(z2)) ?
exact_expt(z1, z2) :
STk_makenumber((double) pow(FLONM(STk_exact2inexact(z1)),
FLONM(STk_exact2inexact(z2))));
}
PRIMITIVE STk_exact2inexact(SCM z)
{
switch (TYPE(z)) {
case tc_integer: return STk_makenumber((double) INTEGER(z));
case tc_bignum: return STk_makenumber(bignum2double(BIGNUM(z)));
case tc_flonum: return z;
default: Err("exact->inexact: bad number", z);
return UNDEFINED; /* never reached */
}
}
PRIMITIVE STk_inexact2exact(SCM z)
{
switch (TYPE(z)) {
case tc_integer:
case tc_bignum: return z;
case tc_flonum: return double2integer(floor(FLONM(z)+0.5));
default: Err("inexact->exact: bad number", z);
return UNDEFINED; /* never reached */
}
}
PRIMITIVE STk_string2number(SCM str, SCM base)
{
long b = (base==UNBOUND)? 10L : STk_integer_value(base);
if (NSTRINGP(str)) Err("string->number: not a string", str);
if (b != 2 && b != 8 && b != 10 && b != 16)
Err("string->number: base must be 2, 8, 10 or 16", NIL);
return STk_Cstr2number(CHARS(str), b);
}
PRIMITIVE STk_number2string(SCM n, SCM base)
{
long b = (base==UNBOUND)? 10L : STk_integer_value(base);
char *s, buffer[100];
SCM z;
if (NNUMBERP(n)) Err("number->string: bad number", n);
if (b != 2 && b != 8 && b != 10 && b != 16)
Err("number->string: base must be 2, 8, 10 or 16", NIL);
s = STk_number2Cstr(n, b, buffer);
z = STk_makestring(s);
if (s != buffer) free(s);
return z;
}
PRIMITIVE STk_bignump(SCM n)
{
return BIGNUMP(n) ? Truth: Ntruth;
}
Reference in New Issue View Git Blame Copy Permalink