number.c 77.5 KB
Newer Older
eg's avatar
eg committed
1 2 3 4 5 6 7 8 9
/*
 *
 * n u m b e r . c				-- Numbers management
 *
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
Erick's avatar
.  
Erick committed
10
 *
eg's avatar
eg committed
11 12 13 14
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
Erick's avatar
.  
Erick committed
15
 *
eg's avatar
eg committed
16 17
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
Erick's avatar
.  
Erick committed
18
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
eg's avatar
eg committed
19 20 21 22 23
 * USA.
 *
 *
 *           Author: Erick Gallesio [eg@kaolin.unice.fr]
 *    Creation date: 12-May-1993 10:34
Erick's avatar
.  
Erick committed
24
 * Last file update:  5-May-2011 17:52 (eg)
eg's avatar
eg committed
25 26 27
 */


Erick's avatar
.  
Erick committed
28 29
/* workaround for bad optimisations done in glibc2.1. Thanks to  Andreas Jaeger
 * <aj@suse.de> for it
eg's avatar
eg committed
30
 */
Erick's avatar
.  
Erick committed
31
#define __NO_MATH_INLINES
eg's avatar
eg committed
32 33 34 35 36 37 38 39 40 41 42 43 44

#include <math.h>
#include <ctype.h>
#include "stklos.h"

#undef sinc
#if defined(__linux__) && defined(__alpha__)
#  include <signal.h>
#endif


/* Real precision */
static int real_precision = REAL_FORMAT_SIZE;
Erick's avatar
.  
Erick committed
45
static int log10_maxint;
eg's avatar
eg committed
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94

#define FINITE_REALP(n) ((REAL_VAL(n) != minus_inf) && (REAL_VAL(n) != plus_inf))



/* Declaration of bignums. This is done here instead of stklos.h to avoid
 * to expose the file "gmp.h" in "stklos.h" which is the interface users
 * see to access all the sytem (note that we can also use our version which
 * can be different of the one which is sytem installed, and resolve conflits
 * could be hard).
 */
#include <gmp.h>

struct bignum_obj {
  stk_header header;
  mpz_t val;
};

#define BIGNUM_VAL(p) 	(((struct bignum_obj *) (p))->val)


/*==============================================================================*/

#define MY_PI		3.1415926535897932384626433832795029L  /* pi */

#define BIGNUM_FITS_INTEGER(_bn) (mpz_cmp_si((_bn), INT_MIN_VAL) >= 0 && 	\
			          mpz_cmp_si((_bn), INT_MAX_VAL) <= 0)
#define LONG_FITS_INTEGER(_l)    (INT_MIN_VAL <= (_l) && (_l) <= INT_MAX_VAL)
#define TYPEOF(n)	         (INTP(n)? tc_integer: STYPE(n))


#define MINUS_INF "-inf.0"
#define PLUS_INF  "+inf.0"
#define NaN	  "nan.0"


/* Special IEEE values */
static double plus_inf, minus_inf;
double STk_NaN;

/**** forward declarations ****/
static type_cell convert(SCM *px, SCM *py);

static int zerop(SCM n);
static int negativep(SCM n);
static int positivep(SCM n);
static int isexactp(SCM z);
static SCM gcd2(SCM o1, SCM o2);

95 96
EXTERN_PRIMITIVE("make-rectangular", make_rectangular, subr2, (SCM r, SCM i));
EXTERN_PRIMITIVE("real-part", real_part, subr1, (SCM z));
eg's avatar
eg committed
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
EXTERN_PRIMITIVE("magnitude", magnitude, subr1, (SCM z));
EXTERN_PRIMITIVE("angle", angle, subr1, (SCM z));
EXTERN_PRIMITIVE("sqrt", sqrt, subr1, (SCM z));
EXTERN_PRIMITIVE("exact->inexact", ex2inex, subr1, (SCM z));
EXTERN_PRIMITIVE("inexact->exact", inex2ex, subr1, (SCM z));

#define add2 STk_add2
#define mul2 STk_mul2
#define div2 STk_div2
#define sub2 STk_sub2
#define absolute STk_abs
#define exact2inexact STk_ex2inex
#define inexact2exact STk_inex2ex


static SCM int_quotient(SCM o1, SCM o2);
static SCM my_cos(SCM z);
static SCM my_sin(SCM z);



/******************************************************************************
Erick's avatar
.  
Erick committed
119 120
 *
 * Utilities
eg's avatar
eg committed
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158
 *
 ******************************************************************************/
static void error_bad_number(SCM n)
{
  STk_error("~S is a bad number", n);
}

static void error_at_least_1(void)
{
  STk_error("expects at least one argument");
}

static void error_not_on_a_complex(SCM n)
{
  STk_error("cannot be computed on the complex number ~S", n);
}

static void error_cannot_operate(char *operation, SCM o1, SCM o2)
{
  STk_error("cannot perform %s on ~S and ~S", operation, o1, o2);
}

static void error_divide_by_0(SCM n)
{
  STk_error("cannot divide ~S by 0", n);
}

static void error_incorrect_radix(SCM r)
{
  STk_error("base must be 2, 8, 10 or 16. It was ~S", r);
}


/*
<doc EXT real-precision
 * (real-precision)
 * (real-precision value)
 *
Erick's avatar
.  
Erick committed
159 160
 * This parameter object permits to change the default precision used
 * to print real numbers.
eg's avatar
eg committed
161 162 163 164 165 166 167 168 169 170 171 172
 * @lisp
 * (real-precision)        => 15
 * (define f 0.123456789)
 * (display f)             @print{} 0.123456789
 * (real-precision 3)
 * (display f)             @print{} 0.123
 * @end lisp
doc>
*/
static SCM real_precision_conv(SCM value)
{
  long precision = STk_integer_value(value);
Erick's avatar
.  
Erick committed
173

eg's avatar
eg committed
174 175 176 177 178 179 180 181 182
  if (precision <= 0 || precision > 50)
    STk_error("real precision must be an integer in ]0 50]. It was ~S",
	      value);
  real_precision = (int) precision;
  return value;
}


/******************************************************************************
Erick's avatar
.  
Erick committed
183 184
 *
 * Constructor Functions
eg's avatar
eg committed
185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
 *
 ******************************************************************************/
static Inline SCM Cmake_complex(SCM r, SCM i)
{
  SCM z;

  NEWCELL(z, complex);
  COMPLEX_REAL(z) = r;
  COMPLEX_IMAG(z) = i;
  return z;
}

static Inline SCM make_complex(SCM r, SCM i)
{
  return (zerop(i)) ? r : Cmake_complex(r, i);
}

static Inline SCM make_polar(SCM a, SCM m)
{
  return make_complex(mul2(a, my_cos(m)), mul2(a, my_sin(m)));
}




static Inline SCM Cmake_rational(SCM n, SCM d)
{
  SCM z;

  NEWCELL(z, rational);
  RATIONAL_NUM(z) = n;
  RATIONAL_DEN(z) = d;
  return z;
}


static SCM make_rational(SCM n, SCM d)
{
  SCM gcd;

  if (zerop(d))
    STk_error("cannot make rational with null denominator");
227

eg's avatar
eg committed
228 229 230 231 232
  /* Always keep sign in the denominator */
  if (negativep(d)) {
    n = mul2(n, MAKE_INT(-1));
    d = mul2(d, MAKE_INT(-1));
  }
Erick's avatar
.  
Erick committed
233

eg's avatar
eg committed
234 235 236 237 238 239 240 241 242 243 244
  /* Simplify rational */
  gcd = gcd2(n, d);
  if (gcd != MAKE_INT(1)) {
    if (d == gcd) return int_quotient(n, gcd);
    n = int_quotient(n, gcd);
    d = int_quotient(d, gcd);
  }

  /* Make rational if denominator is not 1 (if n is small, it is not a bignum) */
  if (d ==  MAKE_INT(1))
    return n;
Erick's avatar
.  
Erick committed
245
  else
eg's avatar
eg committed
246 247 248 249 250
    return Cmake_rational(n, d);
}


/******************************************************************************
Erick's avatar
.  
Erick committed
251 252
 *
 * Types declaration
eg's avatar
eg committed
253 254 255 256 257 258 259 260 261 262 263
 *
 ******************************************************************************/
static void print_bignum(SCM n, SCM port, int mode)
{
  char *s;

  s = STk_must_malloc_atomic(mpz_sizeinbase(BIGNUM_VAL(n), 10) + 2);
  mpz_get_str(s, 10, BIGNUM_VAL(n));
  STk_puts(s, port);
  STk_free(s);
}
Erick's avatar
.  
Erick committed
264

eg's avatar
eg committed
265 266 267 268 269 270 271

static void print_rational(SCM n, SCM port, int mode)
{
  STk_print(RATIONAL_NUM(n), port, mode);
  STk_putc('/', port);
  STk_print(RATIONAL_DEN(n), port, mode);
}
Erick's avatar
.  
Erick committed
272

eg's avatar
eg committed
273 274 275
static void print_complex(SCM n, SCM port, int mode)
{
  SCM imag = COMPLEX_IMAG(n);
Erick's avatar
.  
Erick committed
276

eg's avatar
eg committed
277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302
  STk_print(COMPLEX_REAL(n), port, mode);
  if (positivep(imag))
    STk_putc('+', port);
  STk_print(imag, port, mode);
  STk_putc('i', port);
}



static struct extended_type_descr xtype_bignum = {
  "bignum",
  print_bignum
};

static struct extended_type_descr xtype_complex = {
  "complex",
  print_complex
};

static struct extended_type_descr xtype_rational = {
  "rational",
  print_rational
};


/******************************************************************************
Erick's avatar
.  
Erick committed
303 304
 *
 * Conversion Functions
eg's avatar
eg committed
305 306 307 308 309 310
 *
 ******************************************************************************/

static Inline SCM long2scheme_bignum(long x)
{
  SCM z;
Erick's avatar
.  
Erick committed
311

eg's avatar
eg committed
312 313 314 315 316 317 318 319 320 321 322 323 324 325 326
  NEWCELL(z, bignum);
  mpz_init_set_si(BIGNUM_VAL(z), x);
  return z;
}


static Inline SCM long2integer(long x)
{
  return (INT_MIN_VAL<=x && x<=INT_MAX_VAL) ?  MAKE_INT(x): long2scheme_bignum(x);
}


static Inline SCM double2real(double x)
{
  SCM z;
Erick's avatar
.  
Erick committed
327

eg's avatar
eg committed
328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343
  NEWCELL(z, real);
  REAL_VAL(z) = x;
  return z;
}

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;

  /* Try first to convert n to a long */
  if (((double) INT_MIN_VAL <= n) && (n <= (double) INT_MAX_VAL))
    return MAKE_INT((long) n);

Erick's avatar
.  
Erick committed
344
  /* n doesn't fit in a long => build a bignum. THIS IS VERY INEFFICIENT */
eg's avatar
eg committed
345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383
  tmp = STk_must_malloc(size);
  i = 0;
  if (n < 0.0) { tmp[i++] = '-'; n = -n; }
  do {
    if (i >= size) tmp = STk_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 = STk_Cstr2number(tmp, 10L);
  if (tmp) STk_free(tmp);
  return z;
}

static SCM double2rational(double d)
{
  double fraction, i;
  SCM int_part, num, den, res;
  int negative = 0;

  if (d < 0.0) { negative = 1; d = -d; }
  fraction = modf(d, &i);
  int_part = double2integer(i);

  if (!fraction) {
    res = int_part;
  } else {
    num = MAKE_INT(0);
    den = MAKE_INT(1);
Erick's avatar
.  
Erick committed
384

eg's avatar
eg committed
385 386
    while (fraction) {
      num      = mul2(num, MAKE_INT(2));
Erick's avatar
.  
Erick committed
387
      den      = mul2(den, MAKE_INT(2));
eg's avatar
eg committed
388 389 390 391 392 393
      fraction = modf(ldexp(fraction, 1), &i);
      if (i)
	num = add2(num, MAKE_INT(1));
    }
    res = add2(int_part, div2(num, den));
  }
Erick's avatar
.  
Erick committed
394

eg's avatar
eg committed
395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418
  return negative? mul2(res, MAKE_INT(-1)): res;
}

static Inline SCM bignum2scheme_bignum(mpz_t n)
{
  SCM z;

  NEWCELL(z, bignum);
  mpz_init_set(BIGNUM_VAL(z), n);
  return z;
}

static Inline SCM bignum2integer(mpz_t n)
{
  return MAKE_INT(mpz_get_si(n));
}

static Inline SCM bignum2number(mpz_t n)  /* => int or bignum */
{
  return (BIGNUM_FITS_INTEGER(n)) ? bignum2integer(n): bignum2scheme_bignum(n);
}

static Inline double bignum2double(mpz_t n)
{
Erick's avatar
.  
Erick committed
419 420
  /* I do not use the function mpz_get_d since it gives an unspecified value
   * when converting a number which is +inf or -inf
eg's avatar
eg committed
421 422 423 424
   */
 char *s = STk_must_malloc_atomic(mpz_sizeinbase(n, 10) + 2);
 return atof(mpz_get_str(s, 10, n));
}
Erick's avatar
.  
Erick committed
425

eg's avatar
eg committed
426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470

static Inline double scheme_bignum2double(SCM b)
{
  return bignum2double(BIGNUM_VAL(b));
}


static Inline SCM scheme_bignum2real(SCM bn)
{
  return double2real(scheme_bignum2double(bn));
}


static double rational2double(SCM r)
{
  SCM num = RATIONAL_NUM(r);
  SCM den = RATIONAL_DEN(r);

  switch (convert(&num, &den)) {
    case tc_integer: return ((double) INT_VAL(num)) / ((double) INT_VAL(den));
    case tc_bignum:  return scheme_bignum2double(num) / scheme_bignum2double(den);
    default:         STk_panic("bad rational ~S", r);
  }
  return 0.0; /* never reached */
}

static Inline SCM rational2real(SCM r)
{
  return double2real(rational2double(r));
}

static Inline SCM real2integer(SCM r)
{
  double v = REAL_VAL(r);

  if (floor(v) != v) {
    /* This is not an inexact integer (weak test) */
    STk_error("bad number (~s) in an integer division", r);
  }
  return double2integer(v);
}

void STk_double2Cstr(char *buffer, double n)
{
  sprintf(buffer, "%.*g", real_precision, n);
Erick's avatar
.  
Erick committed
471
  if (strchr(buffer, '.') == NULL && strchr(buffer, 'e') == NULL)
eg's avatar
eg committed
472 473 474 475 476 477 478 479 480 481 482 483
    strcat(buffer, ".0");
  /* Treat special cases of +nan.0 and +inf.0 */
  if (isalpha(buffer[0])) {
    if (strcmp(buffer, "inf.0") == 0) strcpy(buffer, "+inf.0");
    if (strcmp(buffer, "nan.0") == 0) strcpy(buffer, "+nan.0");
  }
}

/* Convert a number to a C-string. Result must be freed if != from buffer */
static char *number2Cstr(SCM n, long base, char buffer[])
{
  char *s = buffer;
Erick's avatar
.  
Erick committed
484

eg's avatar
eg committed
485 486 487 488 489
  switch (TYPEOF(n)) {
    case tc_integer:
      {
	long tmp, val = INT_VAL(n);
	int u;
Erick's avatar
.  
Erick committed
490

eg's avatar
eg committed
491 492 493 494 495 496
	if (val < 0) {
	  val  = -val;
	  *s++ = '-';
	}
	/* Find how much digit we need */
	for (s++, tmp=val; tmp >= base; tmp /= base) s++;
Erick's avatar
.  
Erick committed
497

eg's avatar
eg committed
498 499 500 501 502 503 504 505 506 507 508 509 510 511
	*s = '\0'; tmp = val;
	do {
	  u = tmp % base;
	  *(--s) = u + ((u < 10) ? '0' : 'a'-10);
	  tmp   /= base;
	}
	while (tmp);
	return buffer;
      }
    case tc_bignum:
      s = STk_must_malloc(mpz_sizeinbase(BIGNUM_VAL(n), base) + 2);
      s = mpz_get_str(s, base, BIGNUM_VAL(n));
      return s;
    case tc_rational:
Erick's avatar
.  
Erick committed
512
      {
eg's avatar
eg committed
513
	char *s1, *s2, *s3, tmp[100];
Erick's avatar
.  
Erick committed
514

eg's avatar
eg committed
515 516 517 518 519 520 521 522 523 524
	s1 = number2Cstr(RATIONAL_NUM(n), base, buffer);
	s2 = number2Cstr(RATIONAL_DEN(n), base, tmp);
	s3 = STk_must_malloc(strlen(s1) + strlen(s2) + 2);
	sprintf(s3, "%s/%s", s1, s2);
	if (s2!=tmp) STk_free(s2); /*buffer will event. be deallocated by caller*/
	return s3;
      }
    case tc_complex:
      {
	char *s1, *s2, *s3, tmp[100];
Erick's avatar
.  
Erick committed
525

eg's avatar
eg committed
526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554
	s1 = number2Cstr(COMPLEX_REAL(n), base, buffer);
	s2 = number2Cstr(COMPLEX_IMAG(n), base, tmp);
	s3 = STk_must_malloc(strlen(s1) + strlen(s2) + 3);
	sprintf(s3, "%s%s%si", s1, ((*s2 == '-') ? "": "+"), s2);
	if (s2!=tmp) STk_free(s2); /*buffer will event. be deallocated by caller*/
	return s3;
      }
    case tc_real:
      if (base != 10) STk_error("base must be 10 for this number", n);
      STk_double2Cstr(buffer, REAL_VAL(n));
      return buffer;

    default: return NULL; /* never reached */
  }
}


/*===== The general conversion routine ==== */

static type_cell convert(SCM *px, SCM *py)
{
  SCM x = *px;
  SCM y = *py;

  if (TYPEOF(x)==TYPEOF(y)) return(TYPEOF(x)); /* avoid testing on current cases */
  switch (TYPEOF(x)) {
    case tc_complex:
      	    switch (TYPEOF(y)) {
	      case tc_complex: /*already done */
Erick's avatar
.  
Erick committed
555
	      case tc_real:
eg's avatar
eg committed
556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591
	      case tc_rational:
	      case tc_bignum:
	      case tc_integer: 	*py = Cmake_complex(y, MAKE_INT(0)); break;
	      default:		error_bad_number(y);		     break;
	    }
	    break;
    case tc_real:
      	    switch (TYPEOF(y)) {
	      case tc_complex:  *px = Cmake_complex(x, MAKE_INT(0));    break;
	      case tc_real:    /*already done */ ; 	 	        break;
	      case tc_rational: *py = rational2real(y);  	        break;
	      case tc_bignum:   *py = scheme_bignum2real(y); 	        break;
	      case tc_integer:  *py = double2real((double) INT_VAL(y)); break;
	      default:		error_bad_number(y);		        break;
	    }
	    break;
    case tc_rational:
      	    switch (TYPEOF(y)) {
	      case tc_complex:  *px = Cmake_complex(x, MAKE_INT(0));   break;
	      case tc_real:     *px = rational2real(x);   	       break;
	      case tc_rational: /*already done */ ; 	  	       break;
	      case tc_bignum:   /* no break */
	      case tc_integer:  *py = Cmake_rational(y , MAKE_INT(1)); break;
	      default:	        error_bad_number(y);		       break;
	    }
	    break;
    case tc_bignum:
      	    switch (TYPEOF(y)) {
	      case tc_complex:  *px = Cmake_complex(x, MAKE_INT(0));    break;
	      case tc_real:     *px = scheme_bignum2real(x);   	        break;
	      case tc_rational: *px = Cmake_rational(x , MAKE_INT(1));  break;
	      case tc_bignum:   /* already done */		        break;
	      case tc_integer:  *py = long2scheme_bignum(INT_VAL(y));	break;
	      default:	        error_bad_number(y);		        break;
	    }
	    break;
Erick's avatar
.  
Erick committed
592
    case tc_integer:
eg's avatar
eg committed
593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643
      	    switch (TYPEOF(y)) {
	      case tc_complex:  *px = Cmake_complex(x, MAKE_INT(0));    break;
	      case tc_real:     *px = double2real((double) INT_VAL(x)); break;
	      case tc_rational: *px = Cmake_rational(x,  MAKE_INT(1));  break;
	      case tc_bignum:   *px = long2scheme_bignum(INT_VAL(x));   break;
	      case tc_integer:  /* already done */		        break;
	      default:	        error_bad_number(y);		        break;
	    }
	    break;
    default: error_bad_number(x);
  }
  return TYPEOF(*px);
}


long STk_integer_value(SCM x) /* Returns LONG_MIN if not representable as long */
{
  if (INTP(x)) return INT_VAL(x);
  if (BIGNUMP(x)) {
    mpz_t *v = &BIGNUM_VAL(x);
    if (mpz_cmp_si(*v, LONG_MIN) > 0 && mpz_cmp_si(*v, LONG_MAX) <= 0)
      return  mpz_get_si(*v);
  }
  return LONG_MIN;
}

unsigned long STk_uinteger_value(SCM x) /* Returns ULONG_MAX if not an ulong */
{
  if (INTP(x) && (x > 0)) return INT_VAL(x); /* sign(INTEGER_VAL(x)) == sign(x) */
  if (BIGNUMP(x)) {
    mpz_t *v = &BIGNUM_VAL(x);
    if (mpz_cmp_ui(*v, 0) >= 0 && mpz_cmp_ui(*v, ULONG_MAX) < 0)
      return mpz_get_ui(*v);
  }
  return ULONG_MAX;
}


SCM STk_long2integer(long n)
{
  return long2integer(n);
}


SCM STk_ulong2integer(unsigned long n)
{
  if (n <= INT_MAX_VAL) {  /* n  >= 0 since it is an ulong */
    return MAKE_INT(n);
  }
  else {
    SCM z;
Erick's avatar
.  
Erick committed
644

eg's avatar
eg committed
645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660
    NEWCELL(z, bignum);
    mpz_init_set_ui(BIGNUM_VAL(z), n);
    return z;
  }
}


long STk_integer2int32(SCM n, int *overflow)
{
  *overflow = 0;

  if (INTP(n)) {
#if (LONG_MAX == INT32_MAX)  	  /* longs are on 32 bits */
    return INT_VAL(n);
#else				  /* longs are more than 32 bits (probably 64) */
    long val = INT_VAL(n);
Erick's avatar
.  
Erick committed
661

eg's avatar
eg committed
662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682
    if ((- INT32_MAX - 1) <= val && val < INT32_MAX)
      return val;
    else {
      *overflow = 1;
      return 0;
    }
#endif
  }
  if (BIGNUMP(n)) {
    mpz_t *v = &BIGNUM_VAL(n);
    if (mpz_cmp_si(*v, (- INT32_MAX - 1)) >= 0 && mpz_cmp_si(*v, INT32_MAX) <= 0)
      return mpz_get_si(*v);
  }
  *overflow = 1;
  return 0;
}


unsigned long STk_integer2uint32(SCM n, int *overflow)
{
  *overflow = 0;
Erick's avatar
.  
Erick committed
683

eg's avatar
eg committed
684 685
  if (INTP(n)) {
    long val = INT_VAL(n);
Erick's avatar
.  
Erick committed
686

eg's avatar
eg committed
687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727
    if (val >= 0) {
#if (ULONG_MAX == UINT32_MAX)  	/* unsigned longs are on 32 bits */
      return (unsigned long) INT_VAL(n);
#else				/* longs are more than 32 bits (probably 64) */
      if (val < UINT32_MAX)
	return val;
      else {
	*overflow = 1;
	return 0;
      }
#endif
    }
  }
  if (BIGNUMP(n)) {
    mpz_t *v = &BIGNUM_VAL(n);
    if (mpz_cmp_ui(*v, 0) >= 0 && mpz_cmp_ui(*v, UINT32_MAX) <= 0)
      return mpz_get_ui(*v);
  }
  *overflow = 1;
  return 0;
}



SCM STk_double2real(double d)
{
  return double2real(d); /* use the inlined function */
}


double STk_number2double(SCM n)	/* returns NaN if not convertible */
{
  switch (TYPEOF(n)) {
    case tc_real:     return REAL_VAL(n);
    case tc_rational: return REAL_VAL(rational2real(n));
    case tc_bignum:   return REAL_VAL(scheme_bignum2real(n));
    case tc_integer:  return (double) INT_VAL(n);
    default:          return STk_NaN;
  }
}

Erick's avatar
.  
Erick committed
728

eg's avatar
eg committed
729 730

/******************************************************************************
Erick's avatar
.  
Erick committed
731 732
 *
 * Utilities
eg's avatar
eg committed
733 734 735 736 737 738
 *
 ******************************************************************************/

static long do_compare(SCM x, SCM y)
{
  double d1=0, d2=0;
Erick's avatar
.  
Erick committed
739

eg's avatar
eg committed
740
  switch (TYPEOF(x)) {
Erick's avatar
.  
Erick committed
741
    case tc_real:
eg's avatar
eg committed
742 743
      	    switch (TYPEOF(y)) {
	      case tc_complex:  goto general_diff;
Erick's avatar
.  
Erick committed
744
	      case tc_real:     d1 = REAL_VAL(x); d2 = REAL_VAL(y);
eg's avatar
eg committed
745
				goto double_diff;
Erick's avatar
.  
Erick committed
746
	      case tc_rational:
eg's avatar
eg committed
747 748 749 750 751 752 753 754
	      case tc_bignum:   goto general_diff;
	      case tc_integer:  d1 = REAL_VAL(x); d2 =  INT_VAL(y);
				goto double_diff;
	      default:		break;
	    }
    case tc_integer:
      	    switch (TYPEOF(y)) {
	      case tc_complex:  goto general_diff;
Erick's avatar
.  
Erick committed
755
	      case tc_real:	d1 = INT_VAL(x); d2 = REAL_VAL(y);
eg's avatar
eg committed
756
				goto double_diff;
Erick's avatar
.  
Erick committed
757
	      case tc_rational:
eg's avatar
eg committed
758 759 760 761 762 763 764 765 766 767
	      case tc_bignum:   goto general_diff;
	      case tc_integer:  return (INT_VAL(x) - INT_VAL(y));
	      default:		break;
	    }
    case tc_complex:
    case tc_rational:
    case tc_bignum:
      	    switch (TYPEOF(y)) {
	      case tc_complex:
	      case tc_real:
Erick's avatar
.  
Erick committed
768 769
	      case tc_rational:
	      case tc_bignum:
eg's avatar
eg committed
770 771 772 773 774 775 776 777 778 779
	      case tc_integer:  goto general_diff;
	      default:		break;
	    }
    default: break;
  }
  /* if we are here, it s that x and y cannot be compared */
  STk_error("comparison between ~S and ~S impossible", x,  y);
double_diff:
  if (isnan(d1) && isnan(d2))
    return 0;
Erick's avatar
.  
Erick committed
780
  return (d1 == d2) ? 0 : ((d1 < d2)?  -1 : 1);
eg's avatar
eg committed
781 782 783
general_diff:
  {
    SCM d = sub2(x, y);
Erick's avatar
.  
Erick committed
784

eg's avatar
eg committed
785 786 787 788 789 790 791
    if (zerop(d)) return 0;
    /* complex numbers cannot be compared => return always 1 */
    return COMPLEXP(d) ? 1 : (negativep(d) ? -1: 1);
  }
}


Erick's avatar
.  
Erick committed
792
static SCM int_quotient(SCM x, SCM y)
eg's avatar
eg committed
793 794 795 796 797 798 799
/* Specialized version for rationals. Accepts only integer or bignums as params */
{
  mpz_t q, r;

  if (INTP(x)) {
    if (INTP(y))
      return MAKE_INT(INT_VAL(x)/INT_VAL(y));
Erick's avatar
.  
Erick committed
800
    else
801
      x = long2scheme_bignum(INT_VAL(x));
eg's avatar
eg committed
802 803 804 805 806 807
  } else {
    if (INTP(y))
      y = long2scheme_bignum(INT_VAL(y));
  }
  /* Here x and y are both bignum */
  mpz_init(q); mpz_init(r);
Erick's avatar
.  
Erick committed
808
  mpz_tdiv_qr(q, r, BIGNUM_VAL(x), BIGNUM_VAL(y));
eg's avatar
eg committed
809 810 811 812 813 814 815 816 817 818 819 820 821 822 823
  return bignum2number(q);
}

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);
}


/******************************************************************************
Erick's avatar
.  
Erick committed
824 825
 *
 * Number parser
eg's avatar
eg committed
826 827 828 829 830 831 832 833 834 835 836 837 838 839
 *
 ******************************************************************************/

static SCM compute_exact_real(char *s, char *p1, char *p2, char *p3, char *p4)
{
  SCM int_part, fract_part, exp_part;
  mpz_t tmp;

  mpz_init(tmp);
  int_part   = MAKE_INT(0);
  fract_part = MAKE_INT(0);
  exp_part   = MAKE_INT(1);

  /* Representation of the given number (number is '\0' terminated)
Erick's avatar
.  
Erick committed
840
   *
eg's avatar
eg committed
841 842 843 844 845 846
   *        +xxxxxxxxxxxxxxxxx.yyyyyyyyyyyyyE+zzzzz
   *        ^                 ^^            ^^
   *        |		      ||            ||
   *        +-str          p1-++-p2      p3-++-p4
   */

Erick's avatar
.  
Erick committed
847 848
  /* patch the given string so that splitting the various parts of the number
   * is easy
eg's avatar
eg committed
849 850 851 852 853 854 855 856
   */
  if (p1) *p1 = '\0';
  if (p3) *p3 = '\0';

  if (p1) {		/* compute integer part */
    if (mpz_init_set_str(tmp, s, 10L) < 0) return STk_false;
    int_part = bignum2number(tmp);
  }
Erick's avatar
.  
Erick committed
857

eg's avatar
eg committed
858 859 860 861 862 863 864 865 866 867 868
  if (p3 > p2) {	/* compute decimal part as a rational 0.12 => 6/5 */
    SCM num, den;

    if (mpz_init_set_str(tmp, p2, 10L) < 0) return STk_false;
    num = bignum2number(tmp);

    mpz_ui_pow_ui(tmp, 10UL, strlen(p2));
    den = bignum2number(tmp);

    fract_part = make_rational(num, den);
  }
Erick's avatar
.  
Erick committed
869

eg's avatar
eg committed
870 871
  if (p4) {		/* compute exposant as a rational 3 => 1000, -3 => 1/1000 */
    long expo;
Erick's avatar
.  
Erick committed
872

eg's avatar
eg committed
873 874 875 876 877 878 879 880 881
    expo = atoi(p4);
    if (expo > 0) {
      mpz_ui_pow_ui(tmp, 10UL, expo);
      exp_part = bignum2number(tmp);
    } else {
      mpz_ui_pow_ui(tmp, 10UL, -expo);
      exp_part = div2(MAKE_INT(1), bignum2number(tmp));
    }
  }
Erick's avatar
.  
Erick committed
882

eg's avatar
eg committed
883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923
  /* now return (int_part + fract_part) * exp_part */
  return mul2(add2(int_part, fract_part), exp_part);
}

static SCM read_integer_or_real(char *str, long base, char exact_flag, char **end)
{
  int adigit=0, isint=1;
  char saved_char = '\0', *p = str, *p1, *p2, *p3, *p4;
  SCM res;

  /* see function compute_exact_real for the meaning of these pointers */
  p1 = p2 = p3 = p4 = NULL;

  if (*p == '-' || *p == '+') p+=1;
  if (*p == '#') return STk_false;
  while(digitp(*p, base)) { p+=1; adigit=1; if (*p == '#') isint = 0; }

  if (adigit) p1 = p;		/* p1 = end of integral part */

  if (*p=='.') {
    isint = 0; p += 1;
    p2 = p;
    while(digitp(*p, base)) { p+=1; adigit=1; }
    p3 = p;
  }

  if (!adigit) return STk_false;

  if (*p && strchr("eEsSfFdDlL", *p)) {
    isint = 0;
    p += 1;
    p4 = p;
    if (*p == '-' || *p == '+') p+=1;
    if (!digitp(*p, base)) return STk_false;
    p+=1;
    while (digitp(*p, base)) p+=1;
  }
  if (*p) {
    /* Patch the end of the number with a '\0' (will be restored on exit) */
    saved_char = *p;
    *p = '\0';
Erick's avatar
.  
Erick committed
924
  }
eg's avatar
eg committed
925 926

  if (isint) {
Erick's avatar
.  
Erick committed
927 928 929 930 931
    /* We are sure to have an integer. Read it as a bignum and see if we can
     * convert it in smallnum after that. Small numbers (those with few
     * digits expressed in base 10) are not read as bignums.
     * This optimisation is easily missed (e.g. 000000000000000001 will be
     * read as a bignum), but it avoids allocation for current numbers
932
     * represented in an usual form.
eg's avatar
eg committed
933 934
     */
    mpz_t n;
Erick's avatar
.  
Erick committed
935

eg's avatar
eg committed
936
    if (*str == '+') str+=1; /* mpz_init_set_str doesn't recognize +xyz !!! */
937
    if (strlen(str) <= log10_maxint && base == 10) {
Erick's avatar
.  
Erick committed
938 939
      long num = atol(str);

eg's avatar
eg committed
940 941
      res = (exact_flag == 'i') ? double2real((double) num): MAKE_INT(num);
    }
942 943 944 945 946 947 948
    else {
      if (mpz_init_set_str(n, str, base) < 0) {
	/* Bad syntax for a bignum */
	res = STk_false;
      } else if (BIGNUM_FITS_INTEGER(n)) {
	/* Can be represented as a short integer */
	long num = mpz_get_si(n);
Erick's avatar
.  
Erick committed
949

950 951 952 953 954 955 956 957
	res = (exact_flag == 'i') ? double2real((double) num): MAKE_INT(num);
      } else {
	/* It's a bignum */
	res = bignum2scheme_bignum(n);
	if (exact_flag == 'i') res = scheme_bignum2real(res);
      }
      mpz_clear(n);
    }
eg's avatar
eg committed
958 959 960 961
  } else {
    /* Number is a float */
    if (base == 10) {
      /* Replace sharp signs by 0 */
Erick's avatar
.  
Erick committed
962
      for(p=str; *p; p++)
eg's avatar
eg committed
963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986
	switch (*p) {
	  case '#': *p = '0'; break;
	  case 's': case 'S': case 'f': case 'F':
	  case 'd': case 'D': case 'l': case 'L': *p = 'e';
	}
      if (exact_flag == 'e') {
	res = compute_exact_real(str, p1, p2, p3, p4);
      } else {
	res = double2real(strtod(str, &p));
      }
    }
    else
      res = STk_false;
  }

  if (saved_char) *p = saved_char;  /* character which ended the number */
  *end = p; 			    /* position of last analysed character */
  return res;
}


static SCM read_rational(SCM num, char *str, long base, char exact_flag, char **end)
{
  SCM den;
Erick's avatar
.  
Erick committed
987

eg's avatar
eg committed
988 989
  den = read_integer_or_real(str, base, exact_flag, end);
  if (den == STk_false) return STk_false;
Erick's avatar
.  
Erick committed
990

eg's avatar
eg committed
991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019
  if ((TYPEOF(num) == tc_integer || TYPEOF(num) == tc_bignum) &&
      (TYPEOF(den) == tc_integer || TYPEOF(den) == tc_bignum))
    return make_rational(num, den);
  else
    STk_error("cannot make rational with ~S and ~S", num, den);

  return STk_false; 		/* never reached */
}

SCM STk_Cstr2number(char *str, long base)
{
  int i, exact, radix, polar, is_signed;
  char *p = str;
  SCM num1, num2;

  is_signed = 0;

  if ((*str == '-' || *str == '+')) {
    is_signed = 1;

    if (isalpha(str[1])) {
      /* Treat special values "+inf.0" -inf.0 and "NaN" as well as +i and -i */
      if (strcmp(str, MINUS_INF)==0) return double2real(minus_inf);
      if (strcmp(str, PLUS_INF)==0)  return double2real(plus_inf);
      if (strcmp(str+1, NaN)==0)     return double2real(STk_NaN);
      if (strcmp(str, "+i")==0)      return make_complex(MAKE_INT(0), MAKE_INT(+1));
      if (strcmp(str, "-i")==0)      return make_complex(MAKE_INT(0), MAKE_INT(-1));
    }
  }
Erick's avatar
.  
Erick committed
1020

eg's avatar
eg committed
1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037
  exact = ' ', radix = 0;
  for (i = 0; i < 2; i++) {
    if (*p == '#') {
      p += 1;
      switch (*p++) {
        case 'e': if (exact == ' ') { exact = 'e'; break; }  else return STk_false;
	case 'i': if (exact == ' ') { exact = 'i'; break; }  else return STk_false;
	case 'b': if (!radix) {base = 2;  radix = 1; break;} else return STk_false;
	case 'o': if (!radix) {base = 8;  radix = 1; break;} else return STk_false;
	case 'd': if (!radix) {base = 10; radix = 1; break;} else return STk_false;
	case 'x': if (!radix) {base = 16; radix = 1; break;} else return STk_false;
	default:  return STk_false;
      }
      str += 2;
    }
    if (*p != '#') break;
  }
Erick's avatar
.  
Erick committed
1038

eg's avatar
eg committed
1039 1040
  num1 = read_integer_or_real(p, base, exact, &p);
  if (num1 == STk_false) return STk_false;
Erick's avatar
.  
Erick committed
1041

eg's avatar
eg committed
1042 1043 1044 1045 1046 1047 1048 1049
  if (*p == '/')
    num1 = read_rational(num1, p+1, base, exact, &p);

  if ((*p == '+') || (*p == '-') || (*p == '@')) {
    /* Start to read a complex number */
    if (*p == '+' && p[1] == 'i') {
      p   += 2;
      num1 = make_complex(num1, MAKE_INT(1));	/* special case ...+i */
Erick's avatar
.  
Erick committed
1050
    }
eg's avatar
eg committed
1051 1052 1053 1054 1055 1056
    else if (*p == '-' && p[1] == 'i') {
      p    += 2;
      num1  = make_complex(num1, MAKE_INT(-1));	/* special case ...-i */
    }
    else {					/* general case ....[+-@]... */
      polar = (*p == '@') ? (p++,1) : 0;
Erick's avatar
.  
Erick committed
1057

eg's avatar
eg committed
1058 1059
      num2 = read_integer_or_real(p, base, exact, &p);
      if (num2 == STk_false) return STk_false;
Erick's avatar
.  
Erick committed
1060

eg's avatar
eg committed
1061 1062 1063 1064 1065
      if (*p == '/') {
	/* Second member of complex number is a rational */
	num2 = read_rational(num2, p+1, base, exact, &p);
	if (num2 == STk_false) return STk_false;
      }
Erick's avatar
.  
Erick committed
1066

eg's avatar
eg committed
1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080
      if (polar) {
	num1 = make_polar(num1, num2);
      } else {
	if (*p == 'i') {
	  num1 = make_complex(num1, num2);
	  p += 1;
	}
      }
    }
  } else if (*p == 'i' && is_signed) {
    /* We had a number of the form '{+|-}...i' */
    p   += 1;
    num1 = make_complex(MAKE_INT(0), num1);
  }
Erick's avatar
.  
Erick committed
1081

eg's avatar
eg committed
1082 1083 1084 1085 1086 1087
  return (*p) ? STk_false : num1;
}



/******************************************************************************
Erick's avatar
.  
Erick committed
1088 1089
 *
 * 			Scheme primitives and utilities
eg's avatar
eg committed
1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111
 *
 ******************************************************************************/


/*
<doc number? complex? real? rational? integer?
 * (number? obj)
 * (complex? obj)
 * (real? obj)
 * (rational? obj)
 * (integer? obj)
 *
 * These numerical type predicates can be applied to any kind of
 * argument, including non-numbers. They return |#t| if the object is of
 * the named type, and otherwise they return |#f|. In general, if a type
 * predicate is true of a number then all higher type predicates are
 * also true of that number. Consequently, if a type predicate is
 * false of a number, then all lower type predicates are also false of
 * that number.
 *
 * If |z| is an inexact complex number, then |(real? z)| is true if and
 * only if |(zero? (imag-part z))| is true. If |x| is an inexact real
Erick's avatar
.  
Erick committed
1112
 * number, then |(integer? x)| is true if and only if
eg's avatar
eg committed
1113
 * |(and (finite? x) (= x (round x)))|
Erick's avatar
.  
Erick committed
1114 1115
 *
 *
eg's avatar
eg committed
1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133
 * @lisp
 *   (complex? 3+4i)         =>  #t
 *   (complex? 3)            =>  #t
 *   (real? 3)               =>  #t
 *   (real? -2.5+0.0i)       =>  #t
 *   (real? #e1e10)          =>  #t
 *   (rational? 6/10)        =>  #t
 *   (rational? 6/3)         =>  #t
 *   (integer? 3+0i)         =>  #t
 *   (integer? 3.0)          =>  #t
 *   (integer? 3.2)          =>  #f
 *   (integer? 8/4)          =>  #t
 *   (integer? "no")         =>  #f
 *   (complex? +inf.0)       =>  #t
 *   (real? -inf.0)          =>  #t
 *   (rational? +inf.0)      =>  #f
 *   (integer? -inf.0)       =>  #f
 * @end lisp
Erick's avatar
.  
Erick committed
1134
 *
eg's avatar
eg committed
1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164
doc>
 */
DEFINE_PRIMITIVE("number?", numberp, subr1, (SCM x))
{
  switch (TYPEOF (x)) {
    case tc_complex:
    case tc_real:
    case tc_rational:
    case tc_bignum:
    case tc_integer: return STk_true;
    default:	     return STk_false;
  }
}


DEFINE_PRIMITIVE("complex?", complexp, subr1, (SCM x))
{
  return STk_numberp(x);
}


DEFINE_PRIMITIVE("real?", realp, subr1, (SCM x))
{
  switch (TYPEOF(x)) {
    case tc_complex: return MAKE_BOOLEAN(zerop(COMPLEX_IMAG(x)));
    case tc_real:
    case tc_rational:
    case tc_bignum:
    case tc_integer: return STk_true;
    default:	     return STk_false;
Erick's avatar
.  
Erick committed
1165
  }
eg's avatar
eg committed
1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176
}


DEFINE_PRIMITIVE("rational?", rationalp, subr1, (SCM x))
{
  switch (TYPEOF(x)) {
    case tc_real:    return MAKE_BOOLEAN(FINITE_REALP(x));
    case tc_rational:
    case tc_bignum:
    case tc_integer: return STk_true;
    default:	     return STk_false;
Erick's avatar
.  
Erick committed
1177
  }
eg's avatar
eg committed
1178 1179 1180 1181 1182 1183
}


/*
<doc EXT bignum?
 * (bignum? x)
Erick's avatar
.  
Erick committed
1184 1185 1186
 *
 * This predicates returns |#t| if |x| is an integer number too large to be
 * represented with a native integer.
eg's avatar
eg committed
1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211
 * @lisp
 * (bignum? (expt 2 300))     => |#t|   (very likely)
 * (bignum? 12)               => |#f|
 * (bignum? "no")             => |#f|
 * @end lisp
doc>
a*/
DEFINE_PRIMITIVE("bignum?", bignump, subr1, (SCM x))
{
  return MAKE_BOOLEAN(BIGNUMP(x));
}


DEFINE_PRIMITIVE("integer?", integerp, subr1, (SCM x))
{
  switch (TYPEOF(x)){
    case tc_real:    {
      		       double val = REAL_VAL(x);
		       return ((val == minus_inf) || (val == plus_inf)) ?
			         STk_false:
			         MAKE_BOOLEAN(floor(val) == val);
    		     }
    case tc_bignum:
    case tc_integer: return STk_true;
    default:	     return STk_false;
Erick's avatar
.  
Erick committed
1212
  }
eg's avatar
eg committed
1213 1214 1215 1216 1217 1218 1219 1220 1221 1222
}


/*
<doc  exact? inexact?
 * (exact? z)
 * (inexact? z)
 *
 * These numerical predicates provide tests for the exactness of a
 * quantity. For any Scheme number, precisely one of these predicates
Erick's avatar
.  
Erick committed
1223
 * is true.
eg's avatar
eg committed
1224 1225 1226 1227 1228 1229 1230 1231
doc>
 */

static int isexactp(SCM z)
{
  switch (TYPEOF(z)) {
    case tc_complex:  return isexactp(COMPLEX_REAL(z)) && isexactp(COMPLEX_IMAG(z));
    case tc_real:     return FALSE;
Erick's avatar
.  
Erick committed
1232
    case tc_rational:
eg's avatar
eg committed
1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299
    case tc_bignum:
    case tc_integer:  return TRUE;
    default:	      error_bad_number(z);
  }
  return FALSE; /* never reached */
}


DEFINE_PRIMITIVE("exact?", exactp, subr1, (SCM z))
{
  return MAKE_BOOLEAN(isexactp(z));
}


DEFINE_PRIMITIVE("inexact?", inexactp, subr1, (SCM z))
{
  return MAKE_BOOLEAN(!isexactp(z));
}



/*
<doc  = < > <= >=
 * (= z1 z2 z3 ...)
 * (< x1 x2 x3 ...)
 * (> x1 x2 x3 ...)
 * (<= x1 x2 x3 ...)
 * (>= x1 x2 x3 ...)
 *
 * These procedures return |#t| if their arguments are (respectively):
 * equal, monotonically increasing, monotonically decreasing,
 * monotonically nondecreasing, or monotonically nonincreasing.
 * @lisp
 * (= +inf.0 +inf.0)           =>  #t
 * (= -inf.0 +inf.0)           =>  #f
 * (= -inf.0 -inf.0)           =>  #t
 * 
 * For any finite real number x:
 * 
 * (< -inf.0 x +inf.0)         =>  #t
 * (> +inf.0 x -inf.0)         =>  #t
 * @end lisp
doc>
 */
#define COMPARE_NUM(_sname_, _cname_, _max_type_, _operator_)		    \
    DEFINE_PRIMITIVE(_sname_, _cname_, vsubr, (int argc, SCM *argv))	    \
    {									    \
      SCM previous;							    \
									    \
      if (argc == 0) error_at_least_1();				    \
      if (_max_type_(*argv) == STk_false) error_bad_number(*argv);	    \
  									    \
      for (previous = *argv--; --argc; previous = *argv--) {		    \
        if (_max_type_(*argv) == STk_false) error_bad_number(*argv);	    \
	if (do_compare(previous, *argv) _operator_ 0) return STk_false;     \
      }									    \
      return STk_true;							    \
    }


#define COMPARE_NUM2(_prim_, _max_type_, _operator_)		    	    \
    long STk_##_prim_##2(SCM o1, SCM o2)				    \
    {									    \
      if (_max_type_(o1) == STk_false) error_bad_number(o1);		    \
      if (_max_type_(o2) == STk_false) error_bad_number(o2);		    \
      return do_compare(o1, o2) _operator_ 0;				    \
    }
Erick's avatar
.  
Erick committed
1300

eg's avatar
eg committed
1301 1302 1303 1304 1305 1306 1307 1308 1309

COMPARE_NUM("=",  numeq, STk_complexp, !=)
COMPARE_NUM("<",  numlt, STk_realp,    >=)
COMPARE_NUM(">",  numgt, STk_realp,    <=)
COMPARE_NUM("<=", numle, STk_realp,    >)
COMPARE_NUM(">=", numge, STk_realp,    <)


/* Version with only two parameters (used by runtime) STk_numeq2, STk_numgt2 ... */
Erick's avatar
.  
Erick committed
1310
COMPARE_NUM2(numeq,   STk_complexp, ==)
eg's avatar
eg committed
1311 1312 1313 1314 1315 1316 1317 1318
COMPARE_NUM2(numlt,   STk_realp,    <)
COMPARE_NUM2(numgt,   STk_realp,    >)
COMPARE_NUM2(numle,   STk_realp,    <=)
COMPARE_NUM2(numge,   STk_realp,    >=)



/*
Erick's avatar
.  
Erick committed
1319
<doc finite? infinite?  zero? positive? negative? odd? even?
eg's avatar
eg committed
1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341
 * (finite? z)
 * (infinite? z)
 * (zero? z)
 * (positive? x)
 * (negative? x)
 * (odd? n)
 * (even? n)
 *
 * These numerical predicates test a number for a particular property,
 * returning |#t| or |#f|.
 * @lisp
 * (positive? +inf.0)          ==>  #t
 * (negative? -inf.0)          ==>  #t
 * (finite? -inf.0)            ==>  #f
 * (infinite? +inf.0)          ==>  #t
 * @end lisp
doc>
 */
static Inline int is_oddp(SCM n)
{
  switch (TYPEOF(n)) {
    case tc_integer: return INT_VAL(n) & 1;
1342
    case tc_bignum:  return mpz_odd_p(BIGNUM_VAL(n));
eg's avatar
eg committed
1343 1344 1345 1346 1347 1348 1349 1350
    default:	     error_bad_number(n);
  }
  return FALSE;	/* never reached */
}

static int zerop(SCM n)
{
  switch (TYPEOF(n)) {
1351
    case tc_integer:  return (INT_VAL(n) == 0);
eg's avatar
eg committed
1352 1353
    case tc_real:     return (REAL_VAL(n) == 0.0);
    case tc_bignum:   return (mpz_cmp_si(BIGNUM_VAL(n), 0L) == 0);
1354 1355
    case tc_complex:  return zerop(COMPLEX_REAL(n)) && zerop(COMPLEX_IMAG(n));
    case tc_rational: return zerop(RATIONAL_NUM(n));
eg's avatar
eg committed
1356 1357 1358 1359 1360 1361 1362 1363
    default:	      error_bad_number(n);
  }
  return FALSE;	/* never reached */
}

static int positivep(SCM n)
{
  switch (TYPEOF(n)) {
1364
    case tc_integer:  return (INT_VAL(n) > 0);
eg's avatar
eg committed
1365 1366
    case tc_real:     return (REAL_VAL(n) > 0.0);
    case tc_bignum:   return (mpz_cmp_si(BIGNUM_VAL(n), 0L) > 0);
1367
    case tc_rational: return positivep(RATIONAL_NUM(n));
eg's avatar
eg committed
1368 1369 1370 1371 1372 1373 1374 1375 1376 1377
    case tc_complex:  error_not_on_a_complex(n); break;
    default:	      error_bad_number(n);
  }
  return FALSE;	/* never reached */
}


static int negativep(SCM n)
{
  switch (TYPEOF(n)) {
1378
    case tc_integer:  return (INT_VAL(n) < 0);
eg's avatar
eg committed
1379 1380
    case tc_real:     return (REAL_VAL(n) < 0.0);
    case tc_bignum:   return (mpz_cmp_si(BIGNUM_VAL(n), 0L) < 0);
1381
    case tc_rational: return negativep(RATIONAL_NUM(n));
eg's avatar
eg committed
1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392
    case tc_complex:  error_not_on_a_complex(n); break;
    default:	      error_bad_number(n);
  }
  return FALSE;	/* never reached */
}


static int finitep(SCM n)
{
  switch (TYPEOF(n)) {
    case tc_real:     return (FINITE_REALP(n));
Erick's avatar
.  
Erick committed
1393 1394
    case tc_rational:
    case tc_bignum:
eg's avatar
eg committed
1395
    case tc_integer:  return TRUE;
Erick's avatar
.  
Erick committed
1396
    case tc_complex:  return (finitep(COMPLEX_REAL(n)) &&
eg's avatar
eg committed
1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460
			      finitep(COMPLEX_IMAG(n)));
    default:	      error_bad_number(n);
  }
  return FALSE;	/* never reached */
}


DEFINE_PRIMITIVE("finite?", finitep, subr1, (SCM n))
{
  return MAKE_BOOLEAN(finitep(n));
}


DEFINE_PRIMITIVE("infinite?", infinitep, subr1, (SCM n))
{
  return MAKE_BOOLEAN(!finitep(n));
}


DEFINE_PRIMITIVE("zero?", zerop, subr1, (SCM n))
{
  return MAKE_BOOLEAN(zerop(n));
}


DEFINE_PRIMITIVE("positive?", positivep, subr1, (SCM n))
{
  return MAKE_BOOLEAN(positivep(n));
}


DEFINE_PRIMITIVE("negative?", negativep, subr1, (SCM n))
{
  return MAKE_BOOLEAN(negativep(n));
}


DEFINE_PRIMITIVE("odd?", oddp, subr1, (SCM n))
{
  return MAKE_BOOLEAN(is_oddp(n));
}


DEFINE_PRIMITIVE("even?", evenp, subr1, (SCM n))
{
  return MAKE_BOOLEAN(!is_oddp(n));
}


/*
<doc  max min
 * (max x1 x2 ...)
 * (min x1 x2 ...)
 *
 * These procedures return the maximum or minimum of their arguments.
 *
 * @lisp
 * (max 3 4)              =>  4    ; exact
 * (max 3.9 4)            =>  4.0  ; inexact
 * @end lisp
 * For any real number x:
 * @lisp
 * (max +inf.0 x)         =>  +inf.0
 * (min -inf.0 x)         =>  -inf.0
Erick's avatar
.  
Erick committed
1461
 * @end lisp
eg's avatar
eg committed
1462 1463 1464 1465 1466 1467 1468
 * 
 * ,(bold "Note:") If any argument is inexact, then the result will also be
 * inexact
doc>
 */
DEFINE_PRIMITIVE("max", max, vsubr, (int argc, SCM *argv))
{
Erick's avatar
.  
Erick committed
1469
  SCM res;
eg's avatar
eg committed
1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495
  int exactp;

  if (argc == 0) error_at_least_1();
  if (argc == 1) {
    if (STk_realp(*argv) == STk_true) return *argv;
    error_bad_number(*argv);
  }

  exactp = isexactp(*argv);

  for (res = *argv--; --argc; argv--) {
    /* See that the argument is a correct number */
    if (STk_realp(*argv) == STk_false) error_bad_number(*argv);

    /* determine if result should be exact or not */
    if (!isexactp(*argv)) exactp = 0;

    /* compute max */
    if (do_compare(res, *argv) < 0) res = *argv;
  }
  return (!exactp && isexactp(res)) ? exact2inexact(res) : res;
}


DEFINE_PRIMITIVE("min", min, vsubr, (int argc, SCM *argv))
{
Erick's avatar
.  
Erick committed
1496
  SCM res;
eg's avatar
eg committed
1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550
  int exactp;

  if (argc == 0) error_at_least_1();
  if (argc == 1) {
    if (STk_realp(*argv) == STk_true) return *argv;
    error_bad_number(*argv);
  }

  exactp = isexactp(*argv);

  for (res = *argv--; --argc; argv--) {
    /* See that the argument is a correct number */
    if (STk_realp(*argv) == STk_false) error_bad_number(*argv);

    /* determine if result should be exact or not */
    if (!isexactp(*argv)) exactp = 0;

    /* compute max */
    if (do_compare(res, *argv) > 0) res = *argv;
  }
  return (!exactp && isexactp(res)) ? exact2inexact(res) : res;
}


/*
<doc    + *
 * (+ z1 ...)
 * (* z1 ...)
 *
 * These procedures return the sum or product of their arguments.
 * @lisp
 * (+ 3 4)                 =>  7
 * (+ 3)                   =>  3
 * (+)                     =>  0
 * (+ +inf.0 +inf.0)       =>  +inf.0
 * (+ +inf.0 -inf.0)       =>  +nan.0
 * (* 4)                   =>  4
 * (*)                     =>  1
 * (* 5 +inf.0)            =>  +inf.0
 * (* -5 +inf.0)           =>  -inf.0
 * (* +inf.0 +inf.0)       =>  +inf.0
 * (* +inf.0 -inf.0)       =>  -inf.0
 * (* 0 +inf.0)            =>  +nan.0
 * @end lisp
 * ,(bold "Note:") For any finite number z:
 * @lisp
 * (+ +inf.0 z)            =>  +inf.0
 * (+ -inf.0 z)            =>  -inf.0
 * @end lisp
doc>
 */
SCM STk_add2(SCM o1, SCM o2)
{
  switch (convert(&o1, &o2)) {
Erick's avatar
.  
Erick committed
1551
    case tc_bignum:
eg's avatar
eg committed
1552 1553
        {
	  mpz_t add;
Erick's avatar
.  
Erick committed
1554

eg's avatar
eg committed
1555 1556
	  mpz_init(add);
	  mpz_add(add, BIGNUM_VAL(o1), BIGNUM_VAL(o2));
Erick's avatar
.  
Erick committed
1557

eg's avatar
eg committed
1558 1559 1560 1561
	  o1 = bignum2number(add);
	  mpz_clear(add);
	  break;
	}
Erick's avatar
.  
Erick committed
1562
      case tc_integer:
eg's avatar
eg committed
1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 1598 1599 1600 1601
	{
	  long add =  (long) INT_VAL(o1) + INT_VAL(o2);

	  if (LONG_FITS_INTEGER(add))
	    o1 = MAKE_INT(add);
	  else
	    o1 = long2scheme_bignum(add);
	  break;
	}
      case tc_real:
	{
	  o1 = double2real(REAL_VAL(o1) + REAL_VAL(o2));
	  break;
	}
      case tc_complex:
	{
	  o1 = make_complex(add2(COMPLEX_REAL(o1), COMPLEX_REAL(o2)),
			    add2(COMPLEX_IMAG(o1), COMPLEX_IMAG(o2)));
	  break;
	}
      case tc_rational:
	{
	  SCM num1, num2, den;

	  den  = mul2(RATIONAL_DEN(o1), RATIONAL_DEN(o2));
	  num1 = mul2(RATIONAL_NUM(o1), RATIONAL_DEN(o2));
	  num2 = mul2(RATIONAL_NUM(o2), RATIONAL_DEN(o1));

	  o1 = make_rational(add2(num1, num2), den);
	  break;
	}
      default: error_cannot_operate("addition", o1, o2);
  }
  return o1;
}

DEFINE_PRIMITIVE("+", plus, vsubr, (int argc, SCM *argv))
{
  SCM res;
Erick's avatar
.  
Erick committed
1602

eg's avatar
eg committed
1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618
  if (argc == 0) return MAKE_INT(0);
  if (argc == 1) return add2(MAKE_INT(0), *argv);

  for (res = *argv--; --argc; argv--)
    res = add2(res, *argv);

  return res;
}


/***
 *** multiplication
 ***/
SCM STk_mul2(SCM o1, SCM o2)
{
  switch (convert(&o1, &o2)) {
Erick's avatar
.  
Erick committed
1619
    case tc_bignum:
1620
      mult_bignum:
eg's avatar
eg committed
1621
      {
Erick's avatar
.  
Erick committed
1622 1623
	mpz_t prod;

eg's avatar
eg committed
1624 1625
	mpz_init(prod);
	mpz_mul(prod, BIGNUM_VAL(o1), BIGNUM_VAL(o2));
Erick's avatar
.  
Erick committed
1626

eg's avatar
eg committed
1627 1628 1629 1630
	o1 = bignum2number(prod);
	mpz_clear(prod);
	break;
      }
Erick's avatar
.  
Erick committed
1631
    case tc_integer:
eg's avatar
eg committed
1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654
      {
	long int i1 = INT_VAL(o1);
	long int i2 = INT_VAL(o2);

	o1 = MAKE_INT(i1*i2);
	if (i1 != 0 && (INT_VAL(o1) / i1) != i2) {
	  o1 = long2scheme_bignum(i1);
	  o2 = long2scheme_bignum(i2);
	  goto mult_bignum;
	}
	break;
      }
      case tc_real:
	{
	  o1 = double2real(REAL_VAL(o1) * REAL_VAL(o2));
	  break;
	}
      case tc_complex:
	{
	  SCM r1 = COMPLEX_REAL(o1);
	  SCM i1 = COMPLEX_IMAG(o1);
	  SCM r2 = COMPLEX_REAL(o2);
	  SCM i2 = COMPLEX_IMAG(o2);
Erick's avatar
.  
Erick committed
1655

eg's avatar
eg committed
1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678 1679
	  o1 = make_complex(sub2(mul2(r1,r2), mul2(i1, i2)),
			    add2(mul2(r1,i2), mul2(r2, i1)));
	  break;
	}
      case tc_rational:
	{
	  o1 = make_rational(mul2(RATIONAL_NUM(o1), RATIONAL_NUM(o2)),
			     mul2(RATIONAL_DEN(o1), RATIONAL_DEN(o2)));
	  break;
	}
      default: error_cannot_operate("multiplication", o1, o2);
  }
  return o1;
}

DEFINE_PRIMITIVE("*", multiplication, vsubr, (int argc, SCM *argv))
{
  SCM res;

  if (argc == 0) return MAKE_INT(1);
  if (argc == 1) return mul2(MAKE_INT(1), *argv);

  for (res = *argv--; --argc; argv--)
    res = mul2(res, *argv);
Erick's avatar
.  
Erick committed
1680

eg's avatar
eg committed
1681 1682 1683 1684 1685 1686 1687 1688 1689 1690 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 1708 1709
  return res;
}

/*
<doc   - /
 * (- z)
 * (- z1 z2)
 * (/ z)
 * (/ z1 z2 ...)
 *
 * With two or more arguments, these procedures return the difference or quotient
 * of their arguments, associating to the left. With one argument, however,
 * they return the additive or multiplicative inverse of their argument.
 *
 * @lisp
 * (- 3 4)                 =>  -1
 * (- 3 4 5)               =>  -6
 * (- 3)                   =>  -3
 * (- +inf.0 +inf.0)       => +nan.0
 * (/ 3 4 5)               =>  3/20
 * (/ 3)                   =>  1/3
 * (/ 0.0)                 => +inf.0
 * (/ 0)                   => error (division by 0)
 * @end lisp
doc>
 */
SCM STk_sub2(SCM o1, SCM o2)
{
  switch (convert(&o1, &o2)) {
Erick's avatar
.  
Erick committed
1710
    case tc_bignum:
eg's avatar
eg committed
1711
      {
Erick's avatar
.  
Erick committed
1712 1713
	mpz_t sub;

eg's avatar
eg committed
1714 1715 1716 1717 1718 1719 1720
	mpz_init(sub);
	mpz_sub(sub, BIGNUM_VAL(o1), BIGNUM_VAL(o2));

	o1 = bignum2number(sub),
	mpz_clear(sub);
	break;
      }
Erick's avatar
.  
Erick committed
1721
    case tc_integer:
eg's avatar
eg committed
1722 1723 1724 1725 1726 1727 1728 1729 1730 1731 1732 1733 1734 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1750 1751 1752 1753 1754 1755 1756 1757 1758 1759 1760 1761 1762 1763 1764 1765 1766 1767
      {
	long sub = (long) INT_VAL(o1) - INT_VAL(o2);
	if (LONG_FITS_INTEGER(sub))
	  o1 = MAKE_INT(sub);
	else
	  o1 = long2scheme_bignum(sub);
	break;
      }
      case tc_real:
	{
	  o1 = double2real(REAL_VAL(o1) - REAL_VAL(o2));
	  break;
	}
      case tc_complex:
	{
	  o1 = make_complex(sub2(COMPLEX_REAL(o1), COMPLEX_REAL(o2)),
			    sub2(COMPLEX_IMAG(o1), COMPLEX_IMAG(o2)));
	  break;
	}
      case tc_rational:
	{
	  SCM num1, num2, den;

	  den  = mul2(RATIONAL_DEN(o1), RATIONAL_DEN(o2));
	  num1 = mul2(RATIONAL_NUM(o1), RATIONAL_DEN(o2));
	  num2 = mul2(RATIONAL_NUM(o2), RATIONAL_DEN(o1));

	  o1 = make_rational(sub2(num1, num2), den);
	  break;
	}
      default: error_cannot_operate("substraction", o1, o2);
  }
  return o1;
}


DEFINE_PRIMITIVE("-", difference, vsubr, (int argc, SCM *argv))
{
  SCM res;

  if (argc == 0) error_at_least_1();
  if (argc == 1) return sub2(MAKE_INT(0), *argv);

  for (res = *argv-- ; --argc; argv--)
    res = sub2(res, *argv);
  return res;
Erick's avatar
.  
Erick committed
1768
}
eg's avatar
eg committed
1769 1770 1771 1772 1773 1774 1775 1776


/***
 ***   Division
 ***/
SCM STk_div2(SCM o1, SCM o2)
{
  switch (convert(&o1, &o2)) {
Erick's avatar
.  
Erick committed
1777
    case tc_bignum:
eg's avatar
eg committed
1778 1779 1780 1781 1782 1783
    case tc_integer:
      o1 = make_rational(o1, o2);
      break;
    case tc_real:
      {
	double r2 = REAL_VAL(o2);
Erick's avatar
.  
Erick committed
1784

eg's avatar
eg committed
1785 1786 1787 1788 1789 1790 1791 1792
	if (r2 != 1.0)
	  o1 = double2real(REAL_VAL(o1) / r2);
	break;
      }
    case tc_rational:
      o1 =  make_rational(mul2(RATIONAL_NUM(o1), RATIONAL_DEN(o2)),
			  mul2(RATIONAL_DEN(o1), RATIONAL_NUM(o2)));
      break;
Erick's avatar
.  
Erick committed
1793
    case tc_complex:
eg's avatar
eg committed
1794 1795
      {
	SCM tmp, new_r, new_i;
Erick's avatar
.  
Erick committed
1796

eg's avatar
eg committed
1797 1798 1799 1800 1801 1802 1803 1804 1805 1806 1807 1808 1809 1810 1811 1812 1813 1814 1815 1816 1817 1818 1819 1820 1821 1822 1823 1824 1825 1826 1827 1828 1829 1830 1831 1832 1833 1834 1835 1836 1837 1838 1839 1840 1841 1842 1843 1844
	if (!zerop(o1)) {
	  tmp   = add2(mul2(COMPLEX_REAL(o2), COMPLEX_REAL(o2)),
		       mul2(COMPLEX_IMAG(o2), COMPLEX_IMAG(o2)));
	  new_r = div2(add2(mul2(COMPLEX_REAL(o1), COMPLEX_REAL(o2)),
			    mul2(COMPLEX_IMAG(o1), COMPLEX_IMAG(o2))),
		       tmp);
	  new_i = div2(sub2(mul2(COMPLEX_IMAG(o1), COMPLEX_REAL(o2)),
			    mul2(COMPLEX_REAL(o1), COMPLEX_IMAG(o2))),
		       tmp);
	  o1 = make_complex(new_r, new_i);
	}
	break;
      }
    default: error_cannot_operate("division", o1, o2);
  }
  return o1;
}

DEFINE_PRIMITIVE("/", division, vsubr, (int argc, SCM *argv))
{
  SCM res;

  if (argc == 0) error_at_least_1();
  if (argc == 1) return div2(MAKE_INT(1), *argv);

  for (res = *argv--; --argc; argv--)
    res = div2(res, *argv);
  return res;
}


/*
<doc  abs
 * (abs x)
 *
 * |Abs| returns the absolute value of its argument.
 * @lisp
 * (abs -7)                =>  7
 * (abs -inf.0)            => +inf.0
 * @end lisp
doc>
 */
DEFINE_PRIMITIVE("abs", abs, subr1, (SCM x))
{
  switch (TYPEOF(x)) {
    case tc_integer:  return (INT_VAL(x) < 0) ? MAKE_INT(-INT_VAL(x)) : x;
    case tc_bignum:   if (mpz_cmp_ui(BIGNUM_VAL(x), 0L) < 0) {
      			mpz_t tmp;
1845 1846

			mpz_init(tmp);
eg's avatar
eg committed
1847
			mpz_neg(tmp, BIGNUM_VAL(x));
1848 1849
			x = bignum2scheme_bignum(tmp);
			mpz_clear(tmp);
eg's avatar
eg committed
1850 1851 1852 1853 1854 1855
		      }
		      return x;
    case tc_real:     return (REAL_VAL(x) < 0.0) ? double2real(-REAL_VAL(x)) : x;
    case tc_rational: return make_rational(absolute(RATIONAL_NUM(x)),
					   RATIONAL_DEN(x));
    default:	      error_bad_number(x);
Erick's avatar
.  
Erick committed
1856
  }
eg's avatar
eg committed
1857 1858 1859 1860 1861
  return STk_void;	/* never reached */
}


/*
Erick's avatar
.  
Erick committed
1862
<doc  quotient remainder modulo
eg's avatar
eg committed
1863 1864 1865 1866 1867
 * (quotient n1 n2)
 * (remainder n1 n2)
 * (modulo n1 n2)
 *
 * These procedures implement number-theoretic (integer) division. n2 should
Erick's avatar
.  
Erick committed
1868
 * be non-zero. All three procedures return integers.
eg's avatar
eg committed
1869 1870 1871 1872 1873 1874 1875 1876 1877 1878 1879 1880 1881 1882 1883
 * 
 * If |n1/n2| is an integer:
 *
 * @lisp
 * (quotient n1 n2)   => n1/n2
 * (remainder n1 n2)  => 0
 * (modulo n1 n2)     => 0
 * @end lisp
 *
 * If n1/n2 is not an integer:
 *
 * @lisp
 * (quotient n1 n2)   => nq
 * (remainder n1 n2)  => nr
 * (modulo n1 n2)     => nm
Erick's avatar
.  
Erick committed
1884
 * @end lisp
eg's avatar
eg committed
1885
 *
Erick's avatar
.  
Erick committed
1886 1887
 * where |nq| is |n1/n2| rounded towards zero, 0 < abs(nr) < abs(n2),
 * 0 < abs(nm) < abs(n2), |nr| and |nm| differ from n1 by a multiple of n2,
eg's avatar
eg committed
1888 1889
 * |nr| has the same sign as n1, and |nm| has the same sign as n2.
 * 
Erick's avatar
.  
Erick committed
1890
 * From this we can conclude that for integers |n1| and |n2| with |n2| not
eg's avatar
eg committed
1891 1892 1893 1894 1895 1896 1897 1898 1899 1900
 * equal to 0,
 * @lisp
 *  (= n1 (+ (* n2 (quotient n1 n2))
 *           (remainder n1 n2)))   =>  #t
 * @end lisp
 * provided all numbers involved in that computation are exact.
 *
 * @lisp
 * (modulo 13 4)           =>  1
 * (remainder 13 4)        =>  1
Erick's avatar
.  
Erick committed
1901
 *
eg's avatar
eg committed
1902 1903
 * (modulo -13 4)          =>  3
 * (remainder -13 4)       =>  -1
Erick's avatar
.  
Erick committed
1904
 *
eg's avatar
eg committed
1905 1906
 * (modulo 13 -4)          =>  -3
 * (remainder 13 -4)       =>  1
Erick's avatar
.  
Erick committed
1907
 *
eg's avatar
eg committed
1908 1909
 * (modulo -13 -4)         =>  -1
 * (remainder -13 -4)      =>  -1
Erick's avatar
.  
Erick committed
1910
 *
eg's avatar
eg committed
1911 1912 1913 1914 1915 1916 1917 1918
 * (remainder -13 -4.0)    =>  -1.0  ; inexact
 * @end lisp
doc>
 */
static void integer_division(SCM x, SCM y, SCM *quotient, SCM* remainder)
{
  mpz_t q, r;
  int exact = 1;
Erick's avatar
.  
Erick committed
1919

eg's avatar
eg committed
1920 1921 1922 1923 1924 1925
  if (!INTP(x) && !BIGNUMP(x) && !REALP(x)) error_bad_number(x);
  if (!INTP(y) && !BIGNUMP(y) && !REALP(y)) error_bad_number(y);
  if (zerop(y)) 			    error_divide_by_0(x);

  if (REALP(x)) { x = real2integer(x); exact = 0; }
  if (REALP(y)) { y = real2integer(y); exact = 0; }
Erick's avatar
.  
Erick committed
1926

eg's avatar
eg committed
1927 1928 1929 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945 1946 1947 1948 1949 1950 1951
  /* Here, x and y can only be integer or bignum (not real) */
  if (INTP(x)) {
    if (INTP(y)) {
      long int i1 = INT_VAL(x);
      long int i2 = INT_VAL(y);

      if (exact) {
	*quotient  = MAKE_INT(i1 / i2);
	*remainder = MAKE_INT(i1 % i2);
      } else {
	*quotient  = double2real((double) (i1 / i2));
	*remainder = double2real((double) (i1 % i2));
      }
      return;
    }
    else
      x = long2scheme_bignum(INT_VAL(x));
  } else {
    /* x is a bignum */
    if (INTP(y))
      y = long2scheme_bignum(INT_VAL(y));
  }

  /* Here x and y are both bignum */
  mpz_init(q); mpz_init(r);
Erick's avatar
.  
Erick committed
1952
  mpz_tdiv_qr(q,r,BIGNUM_VAL(x),BIGNUM_VAL(y));
eg's avatar
eg committed
1953 1954 1955 1956
  if (exact) {
    *quotient  = bignum2number(q);
    *remainder = bignum2number(r);
  } else {
1957 1958 1959 1960
    /*     *quotient  = double2real(mpz_get_d(q)); */
    /*     *remainder = double2real(mpz_get_d(r)); */
    *quotient  = double2real(bignum2double(q));
    *remainder = double2real(bignum2double(r));
eg's avatar
eg committed
1961 1962 1963 1964 1965 1966 1967
  }
  mpz_clear(q); mpz_clear(r);                          /* //FIXME: TESTER */
}

DEFINE_PRIMITIVE("quotient", quotient, subr2, (SCM n1, SCM n2))
{
  SCM q, r;
Erick's avatar
.  
Erick committed
1968

eg's avatar
eg committed
1969 1970 1971
  integer_division(n1, n2, &q, &r);
  return q;
}
Erick's avatar
.  
Erick committed
1972

eg's avatar
eg committed
1973 1974 1975 1976

DEFINE_PRIMITIVE("remainder", remainder, subr2, (SCM n1, SCM n2))
{
  SCM q, r;
Erick's avatar
.  
Erick committed
1977

eg's avatar
eg committed
1978 1979 1980 1981 1982 1983 1984 1985
  integer_division(n1, n2, &q, &r);
  return r;
}


DEFINE_PRIMITIVE("modulo", modulo, subr2, (SCM n1, SCM n2))
{
  SCM q, r;
Erick's avatar
.  
Erick committed
1986

eg's avatar
eg committed
1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999
  integer_division(n1, n2, &q, &r);
  if (negativep(n1) != negativep(n2) && !zerop(r))
     /*kerch@parc.xerox.com*/
    r = add2(r, n2);
  return r;
}


/*
<doc  gcd lcm
 * (gcd n1 ...)
 * (lcm n1 ...)
 *
Erick's avatar
.  
Erick committed
2000
 * These procedures return the greatest common divisor or least common
eg's avatar
eg committed
2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 * multiple of their arguments. The result is always non-negative.
 *
 * @lisp
 * (gcd 32 -36)            =>  4
 * (gcd)                   =>  0
 * (lcm 32 -36)            =>  288
 * (lcm 32.0 -36)          =>  288.0  ; inexact
 * (lcm)                   =>  1
 * @end lisp
doc>
 */
static SCM gcd2(SCM n1, SCM n2)
{
  return zerop(n2) ? n1 : gcd2(n2, STk_modulo(n1, n2));
}

DEFINE_PRIMITIVE("gcd", gcd, vsubr, (int argc, SCM *argv))
{
  SCM res;
Erick's avatar
.  
Erick committed
2020

eg's avatar
eg committed
2021 2022 2023 2024 2025
  if (argc == 0) return MAKE_INT(0);
  if (argc == 1) return gcd2(*argv, MAKE_INT(0));

  for (res = *argv--; --argc; argv--)
    res = absolute(gcd2(res, *argv));
Erick's avatar
.  
Erick committed
2026

eg's avatar
eg committed
2027 2028 2029 2030 2031 2032
  return res;
}

DEFINE_PRIMITIVE("lcm", lcm, vsubr, (int argc, SCM *argv))
{
  SCM res, gcd;
Erick's avatar
.  
Erick committed
2033

eg's avatar
eg committed
2034 2035 2036 2037 2038 2039 2040 2041 2042 2043 2044 2045 2046 2047 2048 2049 2050 2051 2052 2053 2054 2055 2056 2057 2058 2059 2060 2061 2062 2063
  if (argc == 0) return MAKE_INT(1);
  if (STk_numberp(*argv) == STk_false) error_bad_number(*argv);

  for (res = *argv--; --argc; argv--) {
    gcd = gcd2(res, *argv);
    res = mul2(res,div2(*argv, gcd));
  }
  return absolute(res);
}

/*
<doc  numerator denominator
 * (numerator q)
 * (denominator q)
 *
 * These procedures return the numerator or denominator of their argument; the
 * result is computed as if the argument was represented as a fraction in
 * lowest terms. The denominator is always positive. The denominator of
 * 0 is defined to be 1.
 * @lisp
 * (numerator (/ 6 4))  =>  3
 * (denominator (/ 6 4))  =>  2
 * (denominator
 * (exact->inexact (/ 6 4))) => 2.0
 * @end lisp
doc>
 */
DEFINE_PRIMITIVE("numerator", numerator, subr1, (SCM q))
{
  switch (TYPEOF(q)) {
Erick's avatar
.  
Erick committed
2064
    case tc_real:     return
eg's avatar
eg committed
2065 2066
			absolute(exact2inexact(STk_numerator(inexact2exact(q))));
    case tc_rational: return absolute(RATIONAL_NUM(q));
Erick's avatar
.  
Erick committed
2067
    case tc_bignum:
eg's avatar
eg committed
2068 2069 2070 2071 2072 2073 2074 2075 2076 2077 2078
    case tc_integer:  return absolute(q);
    default:	      error_bad_number(q);
  }
  return STk_void; /* never reached */
}

DEFINE_PRIMITIVE("denominator", denominator, subr1, (SCM q))
{
  switch (TYPEOF(q)) {
    case tc_real:     return exact2inexact(STk_denominator(inexact2exact(q)));
    case tc_rational: return RATIONAL_DEN(q);
Erick's avatar
.  
Erick committed
2079
    case tc_bignum:
eg's avatar
eg committed
2080 2081 2082 2083 2084 2085 2086 2087 2088 2089 2090 2091 2092 2093
    case tc_integer:  return MAKE_INT(1);
    default:	      error_bad_number(q);
  }
  return STk_void; /* never reached */
}

/*
<doc  floor ceiling truncate round
 * (floor x)
 * (ceiling x)
 * (truncate x)
 * (round x)
 *
 * These procedures return integers. |Floor| returns the largest integer not
Erick's avatar
.  
Erick committed
2094
 * larger than |x|. |Ceiling| returns the smallest integer not smaller than |x|.
eg's avatar
eg committed
2095
 * |Truncate| returns the integer closest to |x| whose absolute value is not
Erick's avatar
.  
Erick committed
2096
 * larger than the absolute value of |x|. |Round| returns the closest integer
eg's avatar
eg committed
2097 2098
 * to |x|, rounding to even when |x| is halfway between two integers.
 * 
Erick's avatar
.  
Erick committed
2099
 * ,(bold "Rationale:") |Round| rounds to even for consistency with the default
eg's avatar
eg committed
2100 2101 2102
 * rounding mode specified by the IEEE floating point standard.
 * 
 * ,(bold "Note:") If the argument to one of these procedures is inexact, then the
Erick's avatar
.  
Erick committed
2103
 * result will also be inexact. If an exact value is needed, the result should
eg's avatar
eg committed
2104 2105 2106 2107 2108 2109 2110 2111
 * be passed to the |inexact->exact| procedure.
 *
 * @lisp
 *
 * (floor -4.3)          =>  -5.0
 * (ceiling -4.3)        =>  -4.0
 * (truncate -4.3)       =>  -4.0
 * (round -4.3)          =>  -4.0
Erick's avatar
.  
Erick committed
2112
 *
eg's avatar
eg committed
2113 2114 2115 2116
 * (floor 3.5)           =>  3.0
 * (ceiling 3.5)         =>  4.0
 * (truncate 3.5)        =>  3.0
 * (round 3.5)           =>  4.0  ; inexact
Erick's avatar
.  
Erick committed
2117
 *
eg's avatar
eg committed
2118 2119 2120 2121 2122 2123 2124 2125 2126
 * (round 7/2)           =>  4    ; exact
 * (round 7)             =>  7
 * @end lisp
doc>
 */
DEFINE_PRIMITIVE("floor", floor, subr1, (SCM x))
{
  switch (TYPEOF(x)) {
    case tc_real:     return double2real(floor(REAL_VAL(x)));
Erick's avatar
.  
Erick committed
2127
    case tc_rational: {
eg's avatar
eg committed
2128
      			SCM tmp;
Erick's avatar
.  
Erick committed
2129

eg's avatar
eg committed
2130
			tmp = negativep(RATIONAL_NUM(x)) ?
Erick's avatar
.  
Erick committed
2131
			  	 sub2(RATIONAL_NUM(x),
eg's avatar
eg committed
2132 2133 2134 2135
				      sub2(RATIONAL_DEN(x), MAKE_INT(1))):
				 RATIONAL_NUM(x);
			return STk_quotient(tmp, RATIONAL_DEN(x));
   		      }
Erick's avatar
.  
Erick committed
2136
    case tc_bignum:
eg's avatar