1/* BEGIN LICENSE BLOCK 2 * Version: CMPL 1.1 3 * 4 * The contents of this file are subject to the Cisco-style Mozilla Public 5 * License Version 1.1 (the "License"); you may not use this file except 6 * in compliance with the License. You may obtain a copy of the License 7 * at www.eclipse-clp.org/license. 8 * 9 * Software distributed under the License is distributed on an "AS IS" 10 * basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See 11 * the License for the specific language governing rights and limitations 12 * under the License. 13 * 14 * The Original Code is The ECLiPSe Constraint Logic Programming System. 15 * The Initial Developer of the Original Code is Cisco Systems, Inc. 16 * Portions created by the Initial Developer are 17 * Copyright (C) 1993-2006 Cisco Systems, Inc. All Rights Reserved. 18 * 19 * Contributor(s): Joachim Schimpf, ECRC and IC-Parc 20 * 21 * END LICENSE BLOCK */ 22 23/* 24 * IDENTIFICATION bigrat.c 25 * 26 * VERSION $Id: bigrat.c,v 1.10 2013/02/10 03:48:37 jschimpf Exp $ 27 * 28 * AUTHOR Joachim Schimpf 29 * 30 * DESCRIPTION bignum and rational arithmetic 31 * Interface to GNUmp library (version 2 or 3) 32 * 33 * All interfacing via struct tag_desc{} 34 * Keep this file free of built-ins !!! 35 * 36 * Has to be linked with libgmp.{a,so,lib} 37 */ 38 39#include "config.h" 40#include "sepia.h" 41#include "types.h" 42#include "error.h" 43#include "embed.h" 44#include "mem.h" 45#include "dict.h" 46#include "emu_export.h" 47#include "rounding_control.h" 48 49 50#ifndef HAVE_LIBGMP 51 52void 53bigrat_init(void) 54{ 55} 56 57extern int 58ec_double_to_int_or_bignum(double f, pword *pres) 59{ 60 if (MIN_S_WORD_DBL <= (f) && (f) < MAX_S_WORD_1_DBL) 61 { 62 pres->val.nint = (word) (f); 63 pres->tag.kernel = TINT; 64 } 65 else 66 { 67 Bip_Error(ARITH_EXCEPTION); 68 } 69 70 Succeed_; 71} 72 73extern int 74ec_array_to_big(const void *p, int count, int order, int size, int endian, unsigned nails, pword *result) 75{ 76 Bip_Error(ARITH_EXCEPTION); 77} 78 79extern int 80ec_big_to_chunks(pword *pw1, uword chunksize, pword *result) 81{ 82 Bip_Error(ARITH_EXCEPTION); 83} 84 85 86#else 87 88 89#include <math.h> 90#include "gmp.h" 91 92#ifndef GMP_LIMB_BITS 93/* for older gmp versions */ 94#define GMP_LIMB_BITS __GMP_BITS_PER_MP_LIMB 95#endif 96 97#ifdef USING_MPIR 98/* MP_INT and MP_RAT are from gmp 1, and gmp provides the 99 following compatibility macros, but MPIR does not. 100 This is a qick work-around, a proper fix would be to update 101 our code to use the more current GMP defs 102*/ 103typedef __mpz_struct MP_INT; 104typedef __mpq_struct MP_RAT; 105#endif 106 107 108#if defined(HAVE_LONG_LONG) || defined(__GNUC__) 109typedef long long long_long; 110#else 111#ifdef HAVE___INT64 112typedef __int64 long_long; 113#endif 114#endif 115 116 117/* 118#define DEBUG_RAT_ALLOC 119*/ 120 121 122/* 123 * Conversions between ECLiPSe representation and gmp representation 124 * 125 * 126 * |---------------| MP_INT: 127 * | | ----------------- 128 * | array of | | size & sign | 129 * | limbs | | alloc | 130 * | | <------ d | 131 * |--------------| |---------------| ----------------- 132 * | TBIG | |BIGSIGN TBUFFER| 133 * | -----------> | bufsize | 134 * |--------------| |---------------| 135 * pword pwords on global 136 * 137 * When we make an MP_INT/MP_RAT from a TBIG/TRAT we share the limb 138 * array on the global stack. Note that such a MP_INT/MP_RAT may not 139 * be assigned to, otherwise the gmp functions try to free the space 140 * to the heap. The sign must be copied from the TBUFFER header. 141 * The macros/functions that convert MP_INT/MP_RAT to Prolog numbers 142 * always clear (free) the MP_INT/MP_RATs, so they must not be 143 * accessed afterwards. 144 */ 145 146#define RatNegative(pw) ((pw)->val.ptr->tag.kernel & BIGSIGN) 147#define BigZero(pw) (BufferSize(pw) == sizeof(mp_limb_t) && \ 148 ((mp_limb_t *) BufferStart(pw))[0] == 0) 149#define RatZero(pw) BigZero(Numer(pw)->val.ptr) 150#define BigPosMin(pw) (BufferSize(pw) == sizeof(mp_limb_t) && !BigNegative(pw) \ 151 && ((mp_limb_t *) BufferStart(pw))[0] == MIN_S_WORD) 152#define BigOne(pw) (BufferSize(pw) == sizeof(mp_limb_t) && \ 153 ((mp_limb_t *) BufferStart(pw))[0] == 1) 154#define RatIntegral(pw) BigOne(Denom(pw)->val.ptr) 155 156#define Duplicate_Buffer(old, new) { /*allow old==new*/ \ 157 pword *_to, *_from = old; \ 158 int _i = BufferPwords(_from); \ 159 (new) = _to = TG; \ 160 TG += _i; \ 161 Check_Gc; \ 162 do { *_to++ = *_from++; } while (--_i); \ 163 } 164 165#define Negate_Big(pbig) \ 166 (pbig)->tag.kernel ^= BIGSIGN; 167 168#define Make_Big(pw, pbig) \ 169 (pw)->val.ptr = (pbig); \ 170 (pw)->tag.kernel = TBIG; 171 172#define Make_Rat(pw, prat) \ 173 (pw)->val.ptr = (prat); \ 174 (pw)->tag.kernel = TRAT; 175 176#define Push_Rat_Frame() Push_List_Frame() 177#define Numer(prat) (prat) 178#define Denom(prat) ((prat)+1) 179 180#define Big_To_Mpi(pbig, mpi) { \ 181 int _limbs = BufferSize(pbig)/sizeof(mp_limb_t); \ 182 (mpi)->_mp_d = (mp_limb_t *) BufferStart(pbig); \ 183 (mpi)->_mp_alloc = _limbs; \ 184 (mpi)->_mp_size = \ 185 (_limbs == 1 && ((mp_limb_t *) BufferStart(pbig))[0] == 0) ? 0 : \ 186 BigNegative(pbig) ? -_limbs : _limbs; \ 187 } 188 189/* Create a properly normalized TINT or TBIG pword from the MP_INT mpi */ 190/* Clears the MP_INT! */ 191#define Pw_From_Mpi(pw, mpi) _pw_from_mpi(pw, mpi) 192 193/* produces a buffer holding the bignum (know not to be zero) */ 194/* Clears the MP_INT! */ 195#define Push_Big_Mpi_Nonzero(mpi) { \ 196 int _size = (MpiNegative(mpi) ? -(mpi)->_mp_size : (mpi)->_mp_size); \ 197 pword *_pbig = TG; \ 198 mp_limb_t *_from, *_to; \ 199 Push_Buffer(_size * sizeof(mp_limb_t)); \ 200 if (MpiNegative(mpi)) { \ 201 (_pbig)->tag.kernel |= BIGSIGN; \ 202 } \ 203 _from = (mpi)->_mp_d; \ 204 _to = (mp_limb_t *) BufferStart(_pbig); \ 205 do { *_to++ = *_from++; } while (--_size); \ 206 mpz_clear(mpi); \ 207 } 208 209/* produces a buffer holding the bignum (may be zero) */ 210/* Clears the MP_INT! */ 211#define Push_Big_Mpi(mpi) \ 212 if (MpiZero(mpi)) { \ 213 pword *_pbig = TG; \ 214 Push_Buffer(sizeof(mp_limb_t)); \ 215 ((mp_limb_t *) BufferStart(_pbig))[0] = 0; \ 216 mpz_clear(mpi); \ 217 } else { \ 218 Push_Big_Mpi_Nonzero(mpi); \ 219 } 220 221/* produces a bignum buffer holding small integer n */ 222#define Push_Big_Int(neg, n) { \ 223 pword *_pbig = TG; \ 224 Push_Buffer(sizeof(mp_limb_t)); \ 225 if (neg) { \ 226 (_pbig)->tag.kernel |= BIGSIGN; \ 227 ((mp_limb_t *) BufferStart(_pbig))[0] = -(n); \ 228 } else { \ 229 ((mp_limb_t *) BufferStart(_pbig))[0] = (n); \ 230 }} 231 232#define Push_Big_PosInt(n) { \ 233 pword *_pbig = TG; \ 234 Push_Buffer(sizeof(mp_limb_t)); \ 235 ((mp_limb_t *) BufferStart(_pbig))[0] = (n); \ 236 } 237 238/* Produces a bignum buffer holding 64 bit integer n. 239 * For 64 bit architectures an mp_limb_t is 64 bit so 240 * Push_Big_Int/Push_Big_PosInt suffices. 241 */ 242#if GMP_LIMB_BITS == 32 243 244#define Push_Big_Int64(neg, n) { \ 245 pword *_pbig = TG; \ 246 Push_Buffer(2 * sizeof(mp_limb_t)); \ 247 if (neg) { \ 248 (_pbig)->tag.kernel |= BIGSIGN; \ 249 ((mp_limb_t *) BufferStart(_pbig))[0] = ((-n) & 0xFFFFFFFF); \ 250 ((mp_limb_t *) BufferStart(_pbig))[1] = (((-n) >> 32) & 0xFFFFFFFF); \ 251 } else { \ 252 ((mp_limb_t *) BufferStart(_pbig))[0] = (n & 0xFFFFFFFF); \ 253 ((mp_limb_t *) BufferStart(_pbig))[1] = ((n >> 32) & 0xFFFFFFFF); \ 254 }} 255 256#define Push_Big_PosInt64(n) { \ 257 pword *_pbig = TG; \ 258 Push_Buffer( 2 * sizeof(mp_limb_t)); \ 259 ((mp_limb_t *) BufferStart(_pbig))[0] = (n & 0xFFFFFFFF); \ 260 ((mp_limb_t *) BufferStart(_pbig))[1] = ((n >> 32) & 0xFFFFFFFF); \ 261 } 262#elif GMP_LIMB_BITS == 64 263 264#define Push_Big_Int64(neg, n) Push_Big_Int(neg, n) 265#define Push_Big_PosInt64(neg, n) Push_Big_PosInt(neg, n) 266 267#else 268 269PROBLEM: No code to cope with MP_LIMB size for this platform! 270 271#endif 272 273#define Push_Big_Copy(old) { \ 274 pword *_from = old; \ 275 pword *_to = TG; \ 276 int _i = BufferPwords(_from); \ 277 TG += _i; \ 278 Check_Gc; \ 279 do { *_to++ = *_from++; } while (--_i); \ 280 } 281 282#define Normalize_Big(pw, pbig) _pw_from_big(pw, pbig) 283 284#define Rat_To_Mpr(prat, mpr) { \ 285 pword *_pw = (prat)[0].val.ptr; \ 286 Big_To_Mpi(_pw, mpq_numref(mpr)); \ 287 _pw = (prat)[1].val.ptr; \ 288 Big_To_Mpi(_pw, mpq_denref(mpr)); \ 289 } 290 291/* Clears the MP_RAT! */ 292#define Push_Rat_Mpr(mpr) { \ 293 pword *_pw = TG; \ 294 Push_List_Frame(); \ 295 Make_Big(_pw, TG); \ 296 Push_Big_Mpi(mpq_numref(mpr)); \ 297 Make_Big(_pw+1, TG); \ 298 Push_Big_Mpi(mpq_denref(mpr)); \ 299 } 300 301/* Create a TRAT pword from the MP_RAT mpr */ 302/* Clears the MP_RAT! */ 303#define Pw_From_Mpr(pw, mpr) { \ 304 Make_Rat(pw, TG); \ 305 Push_Rat_Mpr(mpr); \ 306 } 307 308 309/*-------------------------------------------------------------------------- 310 * Stuff that should be in the gmp library but isn't 311 * It relies on the internal gmp representation 312 *--------------------------------------------------------------------------*/ 313 314#define MpiNegative(x) ((x)->_mp_size < 0) 315#define MpiZero(x) ((x)->_mp_size == 0) 316 317#define ABS(x) (x >= 0 ? x : -x) 318#ifndef HUGE_VAL 319#define HUGE_VAL (limb_value[1]*limb_value[31]) 320#endif 321 322/* A 1024 bit bignum is the largest representable as a double. 323 * That corresponds to 32 32-bit limbs or 16 64-bit limbs. 324 * mpn_to_double looks at the upper 3 limbs (2 if 64-bit limbs), 325 * which is enough for 52-bit double precision, and multiplies 326 * them with their corresponding limb value. 327 */ 328#define LIMB_VALUE_TABLE_SIZE 32 329static double limb_value[LIMB_VALUE_TABLE_SIZE] = { 3301.0, /* 2^(0*32) = 2^(0*64) */ 3314294967296.0, /* 2^(1*32) */ 3321.8446744073709552e+19, /* 2^(2*32) = 2^(1*64) */ 3337.9228162514264338e+28, /* 2^(3*32) */ 3343.4028236692093846e+38, 3351.4615016373309029e+48, 3366.2771017353866808e+57, 3372.695994666715064e+67, 3381.157920892373162e+77, 3394.9732323640978664e+86, 3402.13598703592091e+96, 3419.173994463960286e+105, 3423.9402006196394479e+115, 3431.6923032801030364e+125, 3447.2683872429560689e+134, 3453.1217485503159922e+144, 3461.3407807929942597e+154, 3475.7586096570152914e+163, 3482.4733040147310453e+173, 3491.0622759856335342e+183, 3504.5624406176221952e+192, 3511.959553324262937e+202, 3528.4162174424773976e+211, 3533.6147378671465184e+221, 3541.5525180923007089e+231, 3556.6680144328798543e+240, 3562.8638903918474961e+250, 3571.2300315572313621e+260, 3585.2829453113566525e+269, 3592.269007733883336e+279, 3609.7453140114e+288, /* 2^(30*32) = 2(15*64) */ 3614.1855804968213567e+298 /* 2^(31*32) */ 362}; 363 364static double 365mpn_to_double(mp_ptr d, mp_size_t size) 366{ 367 double res = 0.0; 368 int i = ABS(size); 369#if GMP_LIMB_BITS == 32 370 switch (i) 371 { 372 case 3: 373 res = limb_value[2] * d[2]; 374 /* fall through */ 375 case 2: 376 res += limb_value[1] * d[1]; 377 /* fall through */ 378 case 1: 379 res += (double) d[0]; 380 /* fall through */ 381 case 0: 382 break; 383 default: 384 if (i-3 < LIMB_VALUE_TABLE_SIZE) 385 res = ( (double) d[i-3] 386 + limb_value[1] * d[i-2] 387 + limb_value[2] * d[i-1] 388 ) * limb_value[i-3]; 389 else 390 res = HUGE_VAL; 391 break; 392 } 393#else 394 switch (i) 395 { 396 case 2: 397 res += limb_value[2] * d[1]; 398 /* fall through */ 399 case 1: 400 res += (double) d[0]; 401 /* fall through */ 402 case 0: 403 break; 404 default: 405 if (2*(i-2) < LIMB_VALUE_TABLE_SIZE) 406 res = ( (double) d[i-2] 407 + limb_value[2] * d[i-1] 408 ) * limb_value[2*(i-2)]; 409 else 410 res = HUGE_VAL; 411 break; 412 } 413#endif 414 return size < 0 ? -res : res; 415} 416 417#if __GNU_MP_VERSION < 3 418static double 419mpz_to_double(MP_INT *z) 420{ 421 return mpn_to_double(z->_mp_d, (mp_size_t) z->_mp_size); 422} 423 424static int 425mpz_getbit(MP_INT *z, uword bit) 426{ 427 428 429 mp_size_t size = z->_mp_size; 430 mp_size_t offs = bit / GMP_LIMB_BITS; 431 432 if (size > 0) 433 return offs >= size ? 0 : (z->_mp_d[offs] >> (bit%GMP_LIMB_BITS)) & 1; 434/* 435 else if (size < 0) 436 { 437 Simulate two's complement arithmetic. Not implemented. 438 } 439*/ 440 else 441 return 0; 442} 443 444static void 445mpz_swap(MP_INT *x, MP_INT *y) 446{ 447 MP_INT z; 448 z = *x; 449 *x= *y; 450 *y= z; 451} 452#else 453#define mpz_to_double(z) mpz_get_d(z) 454#define mpz_getbit(z,i) mpz_tstbit(z,i) 455#endif 456 457/* 458 * Divide two bignums giving a double float result. The naive solution 459 * return mpz_to_double(num) / mpz_to_double(den); 460 * suffers from floating point overflows when the numbers are huge 461 * and is inefficient because it looks at unnecessarily many digits. 462 * 463 * IEEE double precision is 53 bits mantissa and 12 bits signed exponent. 464 * So the largest integer representable with doubles is 1024 bits wide, 465 * of which the first 53 are ones, i.e. it lies between 2^1023 and 2^1024. 466 * If the dividend's MSB is more than 1024 bits higher than the divisor's, 467 * the result will always be floating point infinity (no need to divide). 468 * If we do divide, we first drop excess integer precision by keeping only 469 * DBL_PRECISION_LIMBS and ignoring the lower limbs for both operands 470 * (i.e. we effectively scale the integers down, or right-shift them). 471 */ 472 473#define MIN_LIMB_DIFF (1+1024/GMP_NUMB_BITS) 474#define DBL_PRECISION_LIMBS (1+53/GMP_NUMB_BITS) 475 476static double 477mpz_fdiv(MP_INT *num, MP_INT *den) 478{ 479 mp_ptr longer_d, shorter_d; 480 mp_size_t shorter_size, longer_size, ignored_limbs = 0; 481 int negative, swapped; 482 /* By declaring res volatile we make sure that the result is rounded 483 * to double precision instead of being returned with extended precision 484 * in a floating point register, which can have confusing consequences */ 485 volatile double res; 486 487 shorter_size = num->_mp_size; 488 longer_size = den->_mp_size; 489 negative = 0; 490 if (shorter_size < 0) 491 { 492 shorter_size = -shorter_size; 493 negative = !negative; 494 } 495 if (longer_size < 0) 496 { 497 longer_size = -longer_size; 498 negative = !negative; 499 } 500 if (shorter_size > longer_size) 501 { 502 longer_size = shorter_size; 503 longer_d = num->_mp_d; 504 shorter_size = ABS(den->_mp_size); 505 shorter_d = den->_mp_d; 506 swapped = 1; /* abs(res) > 1 */ 507 } 508 else 509 { 510 longer_d = den->_mp_d; 511 shorter_d = num->_mp_d; 512 swapped = 0; /* abs(res) < 1 */ 513 } 514 515 if (longer_size - shorter_size > MIN_LIMB_DIFF) 516 { 517 res = swapped ? HUGE_VAL : 0.0; 518 } 519 else 520 { 521 /* we ignore limbs that are not significant for the result */ 522 if (longer_size > MIN_LIMB_DIFF) /* more can't be represented */ 523 { 524 ignored_limbs = longer_size - MIN_LIMB_DIFF; 525 longer_size -= ignored_limbs; 526 shorter_size -= ignored_limbs; 527 } 528 if (shorter_size > DBL_PRECISION_LIMBS) /* more exceeds the precision */ 529 { 530 ignored_limbs += shorter_size - DBL_PRECISION_LIMBS; 531 longer_size -= shorter_size - DBL_PRECISION_LIMBS; 532 shorter_size = DBL_PRECISION_LIMBS; 533 } 534 longer_d += ignored_limbs; 535 shorter_d += ignored_limbs; 536 537 res = swapped ? 538 mpn_to_double(longer_d, longer_size) 539 / mpn_to_double(shorter_d, shorter_size): 540 mpn_to_double(shorter_d, shorter_size) 541 / mpn_to_double(longer_d, longer_size); 542 543 } 544 return negative ? -res : res; 545} 546 547 548#define HAVE_MPQ_GET_D (__GNU_MP_VERSION >= 3) 549#define MPQ_GET_D_ROUNDS_TOWARDS_ZERO (__GNU_MP_VERSION > 4 || (__GNU_MP_VERSION == 4 && __GNU_MP_VERSION_MINOR >= 2)) 550 551#if HAVE_MPQ_GET_D && !MPQ_GET_D_ROUNDS_TOWARDS_ZERO 552#define mpq_to_double(q) mpq_get_d(q) 553#else 554static double 555mpq_to_double(MP_RAT *q) 556{ 557 return mpz_fdiv(mpq_numref(q), mpq_denref(q)); 558} 559#endif 560 561static void 562mpq_set_double(MP_RAT *q, double f) 563{ 564 extern double floor(); 565 double fabs = (f < 0.0) ? -f : f; /* get rid of the sign */ 566 double x = fabs; 567 MP_INT na, nb, da, db, big_xi, tmpn, tmpd; 568 569 mpz_init_set_ui(&na, 1L); 570 mpz_init_set_ui(&nb, 0L); 571 mpz_init_set_ui(&da, 0L); 572 mpz_init_set_ui(&db, 1L); 573 mpz_init(&tmpn); 574 mpz_init(&tmpd); 575 mpz_init(&big_xi); 576 577 while (mpz_fdiv(&na, &da) != fabs) 578 { 579 double xi = floor(x); 580 double xf = x - xi; 581 582 mpz_swap(&tmpn, &na); 583 mpz_swap(&tmpd, &da); 584 585 if (x < MAX_S_WORD_1_DBL) 586 { 587 uword int_xi = (uword) xi; 588 mpz_mul_ui(&na, &tmpn, int_xi); 589 mpz_mul_ui(&da, &tmpd, int_xi); 590 } 591 else 592 { 593 mpz_set_d(&big_xi, xi); 594 mpz_mul(&na, &tmpn, &big_xi); 595 mpz_mul(&da, &tmpd, &big_xi); 596 } 597 mpz_add(&na, &na, &nb); 598 mpz_add(&da, &da, &db); 599 mpz_swap(&nb, &tmpn); 600 mpz_swap(&db, &tmpd); 601 602 if (xf == 0.0) 603 break; 604 x = 1.0/xf; 605 } 606 607 if (f < 0.0) /* compute q := [+-]na/da */ 608 mpz_neg(&na, &na); 609 mpq_set_num(q, &na); 610 mpq_set_den(q, &da); 611 mpq_canonicalize(q); 612 613 mpz_clear(&big_xi); /* clean up */ 614 mpz_clear(&tmpd); 615 mpz_clear(&tmpn); 616 mpz_clear(&db); 617 mpz_clear(&da); 618 mpz_clear(&nb); 619 mpz_clear(&na); 620} 621 622 623/*-------------------------------------------------------------------------- 624 * memory allocation 625 *--------------------------------------------------------------------------*/ 626 627#ifdef DEBUG_RAT_ALLOC 628static void * 629_rat_alloc(int size) 630{ 631 void *ptr = hp_alloc_size(size); 632 p_fprintf(current_err_, "alloc %d at 0x%x\n", size, ptr); 633 ec_flush(current_err_); 634 return ptr; 635} 636 637static void * 638_rat_realloc(void *ptr, int oldsize, int newsize) 639{ 640 void *newptr = hp_realloc_size(ptr, oldsize, newsize); 641 p_fprintf(current_err_, "reall %d at 0x%x to %d at 0x%x\n", oldsize, ptr, 642 newsize, newptr); 643 ec_flush(current_err_); 644 return newptr; 645} 646 647static void 648_rat_free(void *ptr, int size) 649{ 650 p_fprintf(current_err_, "free %d at 0x%x\n", size, ptr); 651 ec_flush(current_err_); 652 hp_free_size(ptr, size); 653} 654#endif /* DEBUG_RAT_ALLOC */ 655 656static void 657_pw_from_mpi(pword *pw, MP_INT *mpi) 658{ 659 if (mpi->_mp_size == 1) { 660 if (mpi->_mp_d[0] < MIN_S_WORD) { 661 pw->tag.kernel = TINT; 662 pw->val.nint = mpi->_mp_d[0]; 663 mpz_clear(mpi); 664 return; 665 } 666 } else if (mpi->_mp_size == -1) { 667 if (mpi->_mp_d[0] <= MIN_S_WORD) { 668 pw->tag.kernel = TINT; 669 pw->val.nint = - (word) mpi->_mp_d[0]; 670 mpz_clear(mpi); 671 return; 672 } 673 } else if (mpi->_mp_size == 0) { 674 pw->tag.kernel = TINT; 675 pw->val.nint = 0; 676 mpz_clear(mpi); 677 return; 678 } 679 Make_Big(pw, TG); 680 Push_Big_Mpi_Nonzero(mpi); 681} 682 683static void 684_pw_from_big(pword *pw, pword *pbig) 685{ 686 /* BigZero not handled specially here! */ 687 if (BufferSize(pbig) == sizeof(mp_limb_t)) 688 { 689 mp_limb_t i = ((mp_limb_t *) BufferStart(pbig))[0]; 690 if (BigNegative(pbig)) { 691 if (i <= (mp_limb_t) MIN_S_WORD) { 692 pw->tag.kernel = TINT; 693 pw->val.nint = - (word)i; 694 return; 695 } 696 } else { 697 if (i < (mp_limb_t) MIN_S_WORD) { 698 pw->tag.kernel = TINT; 699 pw->val.nint = i; 700 return; 701 } 702 } 703 } 704 Make_Big(pw,pbig); 705} 706 707 708/*-------------------------------------------------------------------------- 709 * methods for the TBIG type 710 *--------------------------------------------------------------------------*/ 711 712/*ARGSUSED*/ 713static int 714_write_big(int quoted, stream_id stream, value vbig, type tbig) 715{ 716 pword *old_tg = TG; 717 int len; 718 value v; 719 char *s; 720 MP_INT a; 721 722 Big_To_Mpi(vbig.ptr, &a); 723 len = mpz_sizeinbase(&a, 10) + (mpz_sgn(&a) < 0 ? 1 : 0); 724 Make_Stack_String(len, v, s) 725 (void) mpz_get_str(s, 10, &a); 726 (void) ec_outfs(stream, s); /* len is not precise! */ 727 TG = old_tg; 728 Succeed_; 729} 730 731/*ARGSUSED*/ 732static int 733_big_string_size(value vb, type tb, int quoted_or_base) /* the result is not exact, may be greater */ 734{ 735 MP_INT a; 736 Big_To_Mpi(vb.ptr, &a); 737 return mpz_sizeinbase(&a, quoted_or_base < 2 ? 10 : quoted_or_base) 738 + (mpz_sgn(&a) < 0 ? 1 : 0); 739} 740 741/*ARGSUSED*/ 742static int 743_big_to_string(value vb, type tb, char *buf, int quoted_or_base) 744{ 745 char *s = buf; 746 MP_INT a; 747 748 Big_To_Mpi(vb.ptr, &a); 749 (void) mpz_get_str(s, quoted_or_base < 2 ? 10 : quoted_or_base, &a); 750 while (*s) s++; 751 return s - buf; 752} 753 754/*ARGSUSED*/ 755static int 756_compare_big(value v1, value v2) 757{ 758 MP_INT a, b; 759 760 Big_To_Mpi(v1.ptr, &a); 761 Big_To_Mpi(v2.ptr, &b); 762 return mpz_cmp(&a, &b); 763} 764 765/*ARGSUSED*/ 766static int 767_arith_compare_big(value v1, value v2, int *res) 768{ 769 MP_INT a, b; 770 771 Big_To_Mpi(v1.ptr, &a); 772 Big_To_Mpi(v2.ptr, &b); 773 *res = mpz_cmp(&a, &b); 774 Succeed_; 775} 776 777/*ARGSUSED*/ 778static int 779_equal_big(pword *pw1, pword *pw2) 780{ 781 MP_INT a, b; 782 783 Big_To_Mpi(pw1, &a); 784 Big_To_Mpi(pw2, &b); 785 return (mpz_cmp(&a, &b) == 0); 786} 787 788/*ARGSUSED*/ 789static int 790_copy_size_big(value v1, type t1) 791{ 792 return BufferPwords(v1.ptr) * sizeof(pword); 793} 794 795 796/*ARGSUSED*/ 797static pword * 798_copy_heap_big(value v1, type t1, pword *top, pword *dest) 799{ 800 int arity; 801 pword *pw; 802 803 dest->val.ptr = top; 804 dest->tag.kernel = Tag(t1.kernel); 805 arity = BufferPwords(v1.ptr); 806 pw = v1.ptr; 807 do /* copy arity pwords */ 808 *top++ = *pw++; 809 while(--arity > 0); 810 return(top); 811} 812 813 814static int 815_big_from_string(char *s, pword *result, int base) 816{ 817 MP_INT a; 818 if (mpz_init_set_str(&a, s, base)) 819 { Bip_Error(BAD_FORMAT_STRING) } 820 Pw_From_Mpi(result, &a); 821 Succeed_; 822} 823 824 825extern int 826ec_array_to_big(const void *p, int count, int order, int size, int endian, unsigned nails, pword *result) 827{ 828#ifndef _WIN32 829 MP_INT z; 830 mpz_init(&z); 831 mpz_import(&z, count, order, size, endian, nails, (const void *) p); 832 Pw_From_Mpi(result, &z); 833#else 834 /* 835 * while we a use using GMP 3, a replacement for: 836 * mpz_import(&z, count, 1, sizeof(mp_limb_t), 0, 0, (const void *) p); 837 */ 838 int i; 839 mp_limb_t *plimbs; 840 Make_Big(result, TG); 841 plimbs = (mp_limb_t *) BufferStart(TG); 842 Push_Buffer(count * sizeof(mp_limb_t)); 843 for(i=0; i<count; ++i) 844 plimbs[i] = ((mp_limb_t *)p)[count-i-1]; 845#endif 846 Succeed_; 847} 848 849 850/* 851 * Break a bignum into integers of chunksize bits 852 */ 853 854#ifndef ULONG_MAX 855#define ULONG_MAX ((unsigned long) -1) 856#endif 857 858extern int 859ec_big_to_chunks(pword *pw1, uword chunksize, pword *result) 860{ 861 unsigned long pos = 0; 862 unsigned long offset = 0; 863 MP_INT z; 864 865 if (chunksize > SIZEOF_WORD*8) 866 { 867 Bip_Error(RANGE_ERROR) 868 } 869 if (BigNegative(pw1)) 870 { 871 Bip_Error(UNIMPLEMENTED); 872 } 873 Big_To_Mpi(pw1, &z); 874 while (pos < ULONG_MAX) 875 { 876 pword *pw = TG; 877 uword chunk = 0; 878 for(;;) 879 { 880 unsigned long lpos; 881 pos = mpz_scan1(&z, pos); 882 lpos = pos-offset; 883 if (lpos >= chunksize) 884 break; 885 chunk |= 1<<lpos; 886 ++pos; 887 } 888 Make_List(result, pw); 889 Push_List_Frame(); 890 Make_Integer(pw, chunk); 891 result = pw+1; 892 offset += chunksize; 893 } 894 Make_Nil(result); 895 Succeed_; 896} 897 898 899/*-------------------------------------------------------------------------- 900 * methods for the TRAT type 901 *--------------------------------------------------------------------------*/ 902 903/*ARGSUSED*/ 904static int 905_write_rat(int quoted, stream_id stream, value vrat, type trat) 906{ 907 int res = _write_big(quoted, stream, vrat.ptr[0].val, vrat.ptr[0].tag); 908 if (res != PSUCCEED) return res; 909 (void) ec_outfc(stream, '_'); 910#ifdef ALT_RAT_SYNTAX 911 (void) ec_outfc(stream, '/'); 912#endif 913 return _write_big(quoted, stream, vrat.ptr[1].val, vrat.ptr[1].tag); 914} 915 916/*ARGSUSED*/ 917static int 918_rat_string_size(value vr, type tr, int quoted) 919{ 920 MP_INT bign; 921 int len; 922 923 Big_To_Mpi(vr.ptr[0].val.ptr, &bign); 924 len = mpz_sizeinbase(&bign, 10) + (mpz_sgn(&bign) < 0 ? 1 : 0); 925 Big_To_Mpi(vr.ptr[1].val.ptr, &bign); 926 len += mpz_sizeinbase(&bign, 10); 927#ifdef ALT_RAT_SYNTAX 928 return len + 2; /* for the "_/" */ 929#else 930 return len + 1; /* for the "_" */ 931#endif 932} 933 934/*ARGSUSED*/ 935static int 936_rat_to_string(value vr, type tr, char *buf, int quoted) 937{ 938 MP_INT bign; 939 char *s = buf; 940 941 Big_To_Mpi(vr.ptr[0].val.ptr, &bign); 942 (void) mpz_get_str(s, 10, &bign); 943 while (*s) s++; 944 *s++ = '_'; 945#ifdef ALT_RAT_SYNTAX 946 *s++ = '/'; 947#endif 948 Big_To_Mpi(vr.ptr[1].val.ptr, &bign); 949 (void) mpz_get_str(s, 10, &bign); 950 while (*s) s++; 951 return s - buf; 952} 953 954/*ARGSUSED*/ 955static int 956_compare_rat(value v1, value v2) 957{ 958 MP_RAT a, b; 959 960 Rat_To_Mpr(v1.ptr, &a); 961 Rat_To_Mpr(v2.ptr, &b); 962 return mpq_cmp(&a, &b); 963} 964 965/*ARGSUSED*/ 966static int 967_arith_compare_rat(value v1, value v2, int *res) 968{ 969 MP_RAT a, b; 970 971 Rat_To_Mpr(v1.ptr, &a); 972 Rat_To_Mpr(v2.ptr, &b); 973 *res = mpq_cmp(&a, &b); 974 Succeed_; 975} 976 977/*ARGSUSED*/ 978static int 979_equal_rat(pword *pw1, pword *pw2) 980{ 981 MP_RAT a, b; 982 983 Rat_To_Mpr(pw1, &a); 984 Rat_To_Mpr(pw2, &b); 985 return mpq_equal(&a, &b); 986} 987 988/*ARGSUSED*/ 989static int 990_copy_size_rat(value v1, type t1) 991{ 992 return _copy_size_big(v1.ptr[0].val, v1.ptr[0].tag) 993 + _copy_size_big(v1.ptr[1].val, v1.ptr[1].tag) 994 + 2 * sizeof(pword); 995} 996 997/*ARGSUSED*/ 998static pword * 999_copy_heap_rat(value v1, type t1, pword *top, pword *dest) 1000{ 1001 dest->tag.kernel = Tag(t1.kernel); 1002 dest->val.ptr = top; 1003 dest = top; 1004 top += 2; 1005 top = _copy_heap_big(v1.ptr[0].val, v1.ptr[0].tag, top, dest); 1006 dest++; 1007 top = _copy_heap_big(v1.ptr[1].val, v1.ptr[1].tag, top, dest); 1008 return top; 1009} 1010 1011static int 1012_rat_from_string(char *s, /* points to valid rational representation */ 1013 pword *result, 1014 int base) 1015{ 1016 MP_RAT a; 1017 char *s1; 1018 1019 mpq_init(&a); 1020 for (s1=s; *s1 != '_'; s1++) 1021 ; 1022 *s1 = '\0'; /* not quite clean ... */ 1023 if (mpz_set_str(mpq_numref(&a), s, base)) { 1024 *s1++ = '_'; 1025 mpq_clear(&a); 1026 Bip_Error(BAD_FORMAT_STRING) 1027 } 1028 *s1++ = '_'; 1029#ifdef ALT_RAT_SYNTAX 1030 if (*s1 == '/') ++s1; /* allow 3_4 or 3_/4 */ 1031#endif 1032 if (mpz_set_str(mpq_denref(&a), s1, base) || MpiZero(mpq_denref(&a))) { 1033 mpq_clear(&a); 1034 Bip_Error(BAD_FORMAT_STRING) 1035 } 1036 mpq_canonicalize(&a); 1037 Pw_From_Mpr(result, &a); 1038 Succeed_; 1039} 1040 1041 1042/*-------------------------------------------------------------------------- 1043 * type coercion functions 1044 *--------------------------------------------------------------------------*/ 1045 1046/*ARGSUSED*/ 1047static int 1048_unimp_err(void) 1049{ 1050 return UNIMPLEMENTED; 1051} 1052 1053static int 1054_int_big(value in, value *out) /* CAUTION: we allow out == &in */ 1055{ 1056 /* this is the only place where we create unnormalized bignums */ 1057 out->ptr = TG; 1058 Push_Big_Int(in.nint < 0, in.nint); 1059 Succeed_; 1060} 1061 1062static int 1063_big_boxlonglong(long_long in, pword *out) 1064{ 1065 if (in >= MIN_S_WORD && in <= MAX_S_WORD) { 1066 Make_Integer(out, in); 1067 } 1068 else { 1069 Make_Big(out, TG); 1070 Push_Big_Int64(in < 0, in); 1071 } 1072 Succeed_; 1073} 1074 1075static int 1076_big_toclonglong(pword *in, long_long *out) 1077{ 1078 int _limbs = BufferSize(in)/sizeof(mp_limb_t); 1079 switch(_limbs) { 1080 case 1: 1081 *out = ((mp_limb_t *) BufferStart(in))[0]; 1082 if (BigNegative(in)) { 1083 *out = -*out; 1084 } 1085 break; 1086 case 2: 1087 *out = ((mp_limb_t *) BufferStart(in))[1]; 1088 if (BigNegative(in)) { 1089 *out = -((*out << 32) | (((mp_limb_t *) BufferStart(in))[0])); 1090 } 1091 else { 1092 *out = ((*out << 32) | ((mp_limb_t *) BufferStart(in))[0]); 1093 } 1094 break; 1095 default: 1096 Bip_Error(INTEGER_OVERFLOW); 1097 break; 1098 } 1099 1100 Succeed_; 1101} 1102 1103 1104static int 1105_int_rat(value in, value *out) /* CAUTION: we allow out == &in */ 1106{ 1107 pword *pw = TG; 1108 Push_Rat_Frame(); 1109 Make_Big(Numer(pw), TG); 1110 Push_Big_Int(in.nint < 0, in.nint); 1111 Make_Big(Denom(pw), TG); 1112 Push_Big_PosInt(1); 1113 out->ptr = pw; 1114 Succeed_; 1115} 1116 1117static int 1118_int_nicerat(value in, pword *pres) 1119{ 1120 pres->tag.kernel = TRAT; 1121 return _int_rat(in, &pres->val); 1122} 1123 1124static int 1125_big_rat(value in, value *out) /* CAUTION: we allow out == &in */ 1126{ 1127 pword *pw = TG; 1128 Push_Rat_Frame(); 1129 Make_Big(Numer(pw), in.ptr); 1130 Make_Big(Denom(pw), TG); 1131 Push_Big_PosInt(1); 1132 out->ptr = pw; 1133 Succeed_; 1134} 1135 1136static int 1137_big_nicerat(value in, pword *pres) 1138{ 1139 pres->tag.kernel = TRAT; 1140 return _big_rat(in, &pres->val); 1141} 1142 1143static int 1144_big_dbl(value in, value *out) /* CAUTION: we allow out == &in */ 1145{ 1146 MP_INT a; 1147 Big_To_Mpi(in.ptr, &a); 1148 Make_Double_Val(*out, mpz_to_double(&a)); 1149 Succeed_; 1150} 1151 1152static int 1153_rat_dbl(value in, value *out) /* CAUTION: we allow out == &in */ 1154{ 1155 MP_RAT a; 1156 Rat_To_Mpr(in.ptr, &a); 1157 Make_Double_Val(*out, mpq_to_double(&a)); 1158 Succeed_; 1159} 1160 1161static int 1162_rat_ivl(value in, value *out) /* CAUTION: we allow out == &in */ 1163{ 1164 pword *pw; 1165 value dval; 1166 MP_RAT a; 1167 Rat_To_Mpr(in.ptr, &a); 1168 1169/* Disabled, because it leads to machine-dependent test results */ 1170/* #if HAVE_MPQ_GET_D && MPQ_GET_D_ROUNDS_TOWARDS_ZERO */ 1171#if 0 1172 Make_Double_Val(dval, mpq_get_d(&a)); 1173 if (Dbl(dval) < 0.0) { 1174 Push_Interval(pw, down(Dbl(dval)), Dbl(dval)); 1175 } else { 1176 Push_Interval(pw, Dbl(dval), up(Dbl(dval))); 1177 } 1178#else 1179 /* mpq_get_d() rounds to nearest, or we are using our own mpq_to_double() */ 1180 Make_Double_Val(dval, mpq_to_double(&a)); 1181 Push_Interval(pw, down(Dbl(dval)), up(Dbl(dval))); 1182#endif 1183 1184 out->ptr = pw; 1185 Succeed_; 1186} 1187 1188/*ARGSUSED*/ 1189static int 1190_dbl_rat(value in, value *out) /* CAUTION: we allow out == &in */ 1191{ 1192 MP_RAT c; 1193 if (!finite(Dbl(in))) 1194 { Bip_Error(ARITH_EXCEPTION); } 1195 mpq_init(&c); 1196 mpq_set_d(&c, Dbl(in)); 1197 out->ptr = TG; 1198 Push_Rat_Mpr(&c); 1199 Succeed_; 1200} 1201 1202static int 1203_dbl_nicerat(value in, pword *pres) 1204{ 1205 MP_RAT c; 1206 if (!finite(Dbl(in))) 1207 { Bip_Error(ARITH_EXCEPTION); } 1208 mpq_init(&c); 1209 mpq_set_double(&c, Dbl(in)); 1210 Make_Rat(pres, TG); 1211 Push_Rat_Mpr(&c); 1212 Succeed_; 1213} 1214 1215 1216/*-------------------------------------------------------------------------- 1217 * Bignum operations 1218 * 1219 * Integers that fit in a word are normally represented as TINT, so TBIGs 1220 * are really big (> 2147483647 or < -2147483648). 1221 * However, for operations with two input arguments, one of them may be an 1222 * unnormalized TBIG due to type coercion. 1223 * Results always have to be normalized. 1224 *--------------------------------------------------------------------------*/ 1225 1226static int 1227_big_nop(value v1, pword *pres) 1228{ 1229 Make_Big(pres, v1.ptr); 1230 Succeed_; 1231} 1232 1233static int 1234_big_add(value v1, value v2, pword *pres) 1235{ 1236 MP_INT a,b,c; 1237 mpz_init(&c); 1238 Big_To_Mpi(v1.ptr, &a); 1239 Big_To_Mpi(v2.ptr, &b); 1240 mpz_add(&c, &a, &b); 1241 Pw_From_Mpi(pres, &c); 1242 Succeed_; 1243} 1244 1245static int 1246_big_next(value v1, pword *pres) 1247{ 1248 MP_INT a,c; 1249 if (BigNegative(v1.ptr)) { 1250 Fail_; 1251 } 1252 mpz_init(&c); 1253 Big_To_Mpi(v1.ptr, &a); 1254 mpz_add_ui(&c, &a, 1L); 1255 Pw_From_Mpi(pres, &c); 1256 Succeed_; 1257} 1258 1259static int 1260_big_prev(value v1, pword *pres) 1261{ 1262 MP_INT a,c; 1263 if (BigNegative(v1.ptr) /*|| BigZero(v1.ptr)*/) { 1264 Fail_; 1265 } 1266 mpz_init(&c); 1267 Big_To_Mpi(v1.ptr, &a); 1268 mpz_sub_ui(&c, &a, 1L); 1269 Pw_From_Mpi(pres, &c); 1270 Succeed_; 1271} 1272 1273static int 1274_big_sub(value v1, value v2, pword *pres) 1275{ 1276 MP_INT a,b,c; 1277 mpz_init(&c); 1278 Big_To_Mpi(v1.ptr, &a); 1279 Big_To_Mpi(v2.ptr, &b); 1280 mpz_sub(&c, &a, &b); 1281 Pw_From_Mpi(pres, &c); 1282 Succeed_; 1283} 1284 1285static int 1286_big_mul(value v1, value v2, pword *pres) 1287{ 1288 MP_INT a,b,c; 1289 mpz_init(&c); 1290 Big_To_Mpi(v1.ptr, &a); 1291 Big_To_Mpi(v2.ptr, &b); 1292 mpz_mul(&c, &a, &b); 1293 Pw_From_Mpi(pres, &c); 1294 Succeed_; 1295} 1296 1297static int 1298_big_idiv(value v1, value v2, pword *pres) 1299{ 1300 MP_INT a,b,c; 1301 if (BigZero(v2.ptr)) 1302 { Bip_Error(ARITH_EXCEPTION); } 1303 mpz_init(&c); 1304 Big_To_Mpi(v1.ptr, &a); 1305 Big_To_Mpi(v2.ptr, &b); 1306 mpz_tdiv_q(&c, &a, &b); 1307 Pw_From_Mpi(pres, &c); 1308 Succeed_; 1309} 1310 1311static int 1312_big_rem(value v1, value v2, pword *pres) 1313{ 1314 MP_INT a,b,c; 1315 if (BigZero(v2.ptr)) 1316 { Bip_Error(ARITH_EXCEPTION); } 1317 mpz_init(&c); 1318 Big_To_Mpi(v1.ptr, &a); 1319 Big_To_Mpi(v2.ptr, &b); 1320 mpz_tdiv_r(&c, &a, &b); 1321 Pw_From_Mpi(pres, &c); 1322 Succeed_; 1323} 1324 1325static int 1326_big_floordiv(value v1, value v2, pword *pres) 1327{ 1328 MP_INT a,b,c; 1329 if (BigZero(v2.ptr)) 1330 { Bip_Error(ARITH_EXCEPTION); } 1331 mpz_init(&c); 1332 Big_To_Mpi(v1.ptr, &a); 1333 Big_To_Mpi(v2.ptr, &b); 1334 mpz_fdiv_q(&c, &a, &b); 1335 Pw_From_Mpi(pres, &c); 1336 Succeed_; 1337} 1338 1339static int 1340_big_floorrem(value v1, value v2, pword *pres) 1341{ 1342 MP_INT a,b,c; 1343 if (BigZero(v2.ptr)) { 1344 Normalize_Big(pres, v1.ptr); 1345 Succeed_; 1346 } 1347 mpz_init(&c); 1348 Big_To_Mpi(v1.ptr, &a); 1349 Big_To_Mpi(v2.ptr, &b); 1350 mpz_fdiv_r(&c, &a, &b); 1351 Pw_From_Mpi(pres, &c); 1352 Succeed_; 1353} 1354 1355static int 1356_big_pow(value v1, value v2, pword *pres) 1357{ 1358 MP_INT a,c; 1359 mpz_init(&c); 1360 Big_To_Mpi(v1.ptr, &a); 1361 if (v2.nint < 0) /* big x neg_int -> rat */ 1362 { 1363 mpz_pow_ui(&c, &a, (uword) (-v2.nint)); 1364 Make_Rat(pres, TG); 1365 Push_Rat_Frame(); 1366 Make_Big(Numer(pres->val.ptr), TG); 1367 Push_Big_PosInt(1); 1368 Make_Big(Denom(pres->val.ptr), TG); 1369 Push_Big_Mpi(&c); 1370 } 1371 else /* big x int -> big */ 1372 { 1373 mpz_pow_ui(&c, &a, (uword) v2.nint); 1374 Pw_From_Mpi(pres, &c); 1375 } 1376 Succeed_; 1377} 1378 1379static int 1380_big_min(value v1, /* CAUTION: One of the inputs may be unnormalized */ 1381 value v2, 1382 pword *pres) 1383{ 1384 MP_INT a,b; 1385 Big_To_Mpi(v1.ptr, &a); 1386 Big_To_Mpi(v2.ptr, &b); 1387 pres->val.ptr = mpz_cmp(&a, &b) < 0 ? v1.ptr : v2.ptr; 1388 Normalize_Big(pres, pres->val.ptr); 1389 Succeed_; 1390} 1391 1392static int 1393_big_max(value v1, /* CAUTION: One of the inputs may be unnormalized */ 1394 value v2, 1395 pword *pres) 1396{ 1397 MP_INT a,b; 1398 Big_To_Mpi(v1.ptr, &a); 1399 Big_To_Mpi(v2.ptr, &b); 1400 pres->val.ptr = mpz_cmp(&a, &b) > 0 ? v1.ptr : v2.ptr; 1401 Normalize_Big(pres, pres->val.ptr); 1402 Succeed_; 1403} 1404 1405static int 1406_big_neg(value v1, /* can't be zero */ 1407 pword *pres) 1408{ 1409 pword *pw; 1410 if (BigPosMin(v1.ptr)) { 1411 Make_Integer(pres, MIN_S_WORD); 1412 } else { 1413 Duplicate_Buffer(v1.ptr, pw); 1414 Negate_Big(pw); 1415 Make_Big(pres, pw); 1416 } 1417 Succeed_; 1418} 1419 1420static int 1421_big_abs(value v1, pword *pres) 1422{ 1423 if (BigNegative(v1.ptr)) 1424 { 1425 return _big_neg(v1, pres); 1426 } 1427 Make_Big(pres, v1.ptr); 1428 Succeed_; 1429} 1430 1431static int 1432_big_sgn(value v1, /* can't be zero */ 1433 pword *pres) 1434{ 1435 Make_Integer(pres, BigNegative(v1.ptr) ? -1: 1); 1436 Succeed_; 1437} 1438 1439static int 1440_big_and(value v1, value v2, pword *pres) 1441{ 1442 MP_INT a,b,c; 1443 mpz_init(&c); 1444 Big_To_Mpi(v1.ptr, &a); 1445 Big_To_Mpi(v2.ptr, &b); 1446 mpz_and(&c, &a, &b); 1447 Pw_From_Mpi(pres, &c); 1448 Succeed_; 1449} 1450 1451static int 1452_big_or(value v1, value v2, pword *pres) 1453{ 1454 MP_INT a,b,c; 1455 mpz_init(&c); 1456 Big_To_Mpi(v1.ptr, &a); 1457 Big_To_Mpi(v2.ptr, &b); 1458 mpz_ior(&c, &a, &b); 1459 Pw_From_Mpi(pres, &c); /* might be small negative number */ 1460 Succeed_; 1461} 1462 1463static int 1464_big_xor(value v1, value v2, pword *pres) 1465{ 1466#if __GNU_MP_VERSION < 3 1467 Bip_Error(UNIMPLEMENTED); 1468#else 1469 MP_INT a,b,c; 1470 mpz_init(&c); 1471 Big_To_Mpi(v1.ptr, &a); 1472 Big_To_Mpi(v2.ptr, &b); 1473 mpz_xor(&c, &a, &b); 1474 Pw_From_Mpi(pres, &c); 1475 Succeed_; 1476#endif 1477} 1478 1479static int 1480_big_com(value v1, pword *pres) 1481{ 1482 MP_INT a,c; 1483 mpz_init(&c); 1484 Big_To_Mpi(v1.ptr, &a); 1485 mpz_com(&c, &a); 1486 Pw_From_Mpi(pres, &c); 1487 Succeed_; 1488} 1489 1490static int 1491_big_shl(value v1, value v2, pword *pres) /* big x int -> big */ 1492{ 1493 MP_INT a,c; 1494 mpz_init(&c); 1495 Big_To_Mpi(v1.ptr, &a); 1496 if (v2.nint >= 0) 1497 mpz_mul_2exp(&c, &a, (uword) v2.nint); 1498 else 1499#if __GNU_MP_VERSION < 3 1500 mpz_div_2exp(&c, &a, (uword) -v2.nint); 1501#else 1502 mpz_fdiv_q_2exp(&c, &a, (uword) -v2.nint); 1503#endif 1504 Pw_From_Mpi(pres, &c); 1505 Succeed_; 1506} 1507 1508static int 1509_big_shr(value v1, value v2, pword *pres) /* big x int -> big */ 1510{ 1511 MP_INT a,c; 1512 mpz_init(&c); 1513 Big_To_Mpi(v1.ptr, &a); 1514 if (v2.nint >= 0) 1515#if __GNU_MP_VERSION < 3 1516 mpz_div_2exp(&c, &a, (uword) v2.nint); 1517#else 1518 mpz_fdiv_q_2exp(&c, &a, (uword) v2.nint); 1519#endif 1520 else 1521 mpz_mul_2exp(&c, &a, (uword) -v2.nint); 1522 Pw_From_Mpi(pres, &c); 1523 Succeed_; 1524} 1525 1526static int 1527_big_setbit(value v1, value v2, pword *pres) /* big x int -> big */ 1528{ 1529 MP_INT a,c; 1530 Big_To_Mpi(v1.ptr, &a); 1531 mpz_init_set(&c, &a); 1532 mpz_setbit(&c, (uword) v2.nint); 1533 Pw_From_Mpi(pres, &c); 1534 Succeed_; 1535} 1536 1537static int 1538_big_clrbit(value v1, value v2, pword *pres) /* big x int -> big */ 1539{ 1540 MP_INT a,c; 1541 Big_To_Mpi(v1.ptr, &a); 1542 mpz_init_set(&c, &a); 1543 mpz_clrbit(&c, (uword) v2.nint); 1544 Pw_From_Mpi(pres, &c); 1545 Succeed_; 1546} 1547 1548static int 1549_big_getbit(value v1, value v2, pword *pres) /* big x int -> int */ 1550{ 1551 MP_INT a; 1552 if (BigNegative(v1.ptr)) 1553 { Bip_Error(UNIMPLEMENTED); } 1554 Big_To_Mpi(v1.ptr, &a); 1555 Make_Integer(pres, mpz_getbit(&a, v2.nint)); 1556 Succeed_; 1557} 1558 1559static int 1560_big_gcd(value v1, value v2, pword *pres) 1561{ 1562 MP_INT a,b,c; 1563 mpz_init(&c); 1564 Big_To_Mpi(v1.ptr, &a); 1565 Big_To_Mpi(v2.ptr, &b); 1566 mpz_gcd(&c, &a, &b); 1567 Pw_From_Mpi(pres, &c); 1568 Succeed_; 1569} 1570 1571static int 1572_big_gcd_ext(value v1, value v2, pword *ps, pword *pt, pword *pg) 1573{ 1574 MP_INT a,b,g,s,t; 1575 mpz_init(&g); 1576 mpz_init(&s); 1577 mpz_init(&t); 1578 Big_To_Mpi(v1.ptr, &a); 1579 Big_To_Mpi(v2.ptr, &b); 1580 mpz_gcdext(&g, &s, &t, &a, &b); 1581 1582 Pw_From_Mpi(ps, &s); 1583 Pw_From_Mpi(pt, &t); 1584 Pw_From_Mpi(pg, &g); 1585 1586 Succeed_; 1587} 1588 1589 1590static int 1591_big_lcm(value v1, value v2, pword *pres) 1592{ 1593 MP_INT a,b,c; 1594 mpz_init(&c); 1595 Big_To_Mpi(v1.ptr, &a); 1596 Big_To_Mpi(v2.ptr, &b); 1597 mpz_lcm(&c, &a, &b); 1598 Pw_From_Mpi(pres, &c); 1599 Succeed_; 1600} 1601 1602static int 1603_big_powm(value vbase, value vexp, value vmod, pword *pres) 1604{ 1605 MP_INT a,b,c,d; 1606 if (BigNegative(vexp.ptr)) { Bip_Error(UNIMPLEMENTED); } 1607 Big_To_Mpi(vbase.ptr, &a); 1608 Big_To_Mpi(vexp.ptr, &b); 1609 Big_To_Mpi(vmod.ptr, &c); 1610 mpz_init(&d); 1611 mpz_powm(&d, &a, &b, &c); 1612 Pw_From_Mpi(pres, &d); 1613 Succeed_; 1614} 1615 1616static int 1617_big_atan2(value v1, value v2, pword *pres) 1618{ 1619 MP_INT a,b; 1620 Big_To_Mpi(v1.ptr, &a); 1621 Big_To_Mpi(v2.ptr, &b); 1622 Make_Double(pres, Atan2(mpz_to_double(&a), mpz_to_double(&b))); 1623 Succeed_; 1624} 1625 1626 1627/*-------------------------------------------------------------------------- 1628 * Miscellaneous operations 1629 *--------------------------------------------------------------------------*/ 1630 1631static int 1632_int_neg(value v1, pword *pres) 1633{ 1634 if (v1.nint == MIN_S_WORD) 1635 { 1636 Make_Big(pres, TG); 1637 Push_Big_PosInt(MIN_S_WORD); 1638 } 1639 else 1640 { 1641 Make_Integer(pres, -v1.nint); 1642 } 1643 Succeed_; 1644} 1645 1646/* 1647 * We want to provide a double -> TINT or TBIG function with two slightly 1648 * different interfaces: one suitable for the arith_op table's ARITH_FIX 1649 * operator, and one suitable for external code to call which doesn't 1650 * require the double to be in TDBL packaging. 1651 */ 1652 1653#define Dbl_Fix(f, pres) { \ 1654 if (MIN_S_WORD_DBL <= (f) && (f) < MAX_S_WORD_1_DBL) \ 1655 { \ 1656 pres->val.nint = (word) (f); \ 1657 pres->tag.kernel = TINT; \ 1658 } \ 1659 else if (finite(f)) /* convert to a bignum */ \ 1660 { \ 1661 MP_INT c; \ 1662 mpz_init(&c); \ 1663 mpz_set_d(&c, f); \ 1664 Make_Big(pres, TG); \ 1665 Push_Big_Mpi_Nonzero(&c); \ 1666 } \ 1667 else \ 1668 { \ 1669 Bip_Error(ARITH_EXCEPTION); \ 1670 } \ 1671 } 1672 1673static int 1674_dbl_fix(value v1, pword *pres) 1675{ 1676 double f = Dbl(v1); 1677 Dbl_Fix(f, pres); 1678 Succeed_; 1679} 1680 1681static int 1682_dbl_int2(value v1, pword *pres) 1683{ 1684 double ipart; 1685 double fpart = modf(Dbl(v1), &ipart); 1686 1687 if (fpart != 0.0) 1688 { 1689 Bip_Error(ARITH_EXCEPTION); 1690 } 1691 else if (MIN_S_WORD_DBL <= (ipart) && (ipart) < MAX_S_WORD_1_DBL) 1692 { 1693 pres->val.nint = (word) ipart; 1694 pres->tag.kernel = TINT; 1695 } 1696 else if (finite(ipart)) 1697 { 1698 MP_INT c; 1699 mpz_init(&c); 1700 mpz_set_d(&c, ipart); 1701 Make_Big(pres, TG); 1702 Push_Big_Mpi_Nonzero(&c); 1703 } 1704 else 1705 { 1706 Bip_Error(ARITH_EXCEPTION); 1707 } 1708 Succeed_; 1709} 1710 1711extern int 1712ec_double_to_int_or_bignum(double f, pword *pres) 1713{ 1714 Dbl_Fix(f, pres); 1715 Succeed_; 1716} 1717 1718/*-------------------------------------------------------------------------- 1719 * Rational operations 1720 *--------------------------------------------------------------------------*/ 1721 1722static int 1723_rat_nop(value v1, pword *pres) 1724{ 1725 Make_Rat(pres, v1.ptr); 1726 Succeed_; 1727} 1728 1729static int 1730_rat_neg(value v1, pword *pres) 1731{ 1732 if (RatZero(v1.ptr)) 1733 { 1734 Make_Rat(pres, v1.ptr); 1735 } 1736 else 1737 { 1738 pword *pw = TG; 1739 Push_Rat_Frame(); 1740 Make_Big(Numer(pw), TG); 1741 Push_Big_Copy(Numer(v1.ptr)->val.ptr); 1742 Negate_Big(Numer(pw)->val.ptr); 1743 Make_Big(Denom(pw), Denom(v1.ptr)->val.ptr); 1744 Make_Rat(pres, pw); 1745 } 1746 Succeed_; 1747} 1748 1749static int 1750_rat_abs(value v1, pword *pres) 1751{ 1752 if (RatNegative(v1.ptr)) 1753 { 1754 return _rat_neg(v1, pres); 1755 } 1756 Make_Rat(pres, v1.ptr); 1757 Succeed_; 1758} 1759 1760static int 1761_rat_sgn(value v1, pword *pres) 1762{ 1763 Make_Integer(pres, RatNegative(v1.ptr) ? -1: RatZero(v1.ptr)? 0: 1); 1764 Succeed_; 1765} 1766 1767static int 1768_rat_add(value v1, value v2, pword *pres) 1769{ 1770 MP_RAT a,b,c; 1771 mpq_init(&c); 1772 Rat_To_Mpr(v1.ptr, &a); 1773 Rat_To_Mpr(v2.ptr, &b); 1774 mpq_add(&c, &a, &b); 1775 Pw_From_Mpr(pres, &c); 1776 Succeed_; 1777} 1778 1779static int 1780_rat_sub(value v1, value v2, pword *pres) 1781{ 1782 MP_RAT a,b,c; 1783 mpq_init(&c); 1784 Rat_To_Mpr(v1.ptr, &a); 1785 Rat_To_Mpr(v2.ptr, &b); 1786 mpq_sub(&c, &a, &b); 1787 Pw_From_Mpr(pres, &c); 1788 Succeed_; 1789} 1790 1791static int 1792_rat_mul(value v1, value v2, pword *pres) 1793{ 1794 MP_RAT a,b,c; 1795 mpq_init(&c); 1796 Rat_To_Mpr(v1.ptr, &a); 1797 Rat_To_Mpr(v2.ptr, &b); 1798 mpq_mul(&c, &a, &b); 1799 Pw_From_Mpr(pres, &c); 1800 Succeed_; 1801} 1802 1803/* only used if PREFER_RATIONALS */ 1804static int 1805_int_div(value v1, value v2, pword *pres) /* int x int -> rat */ 1806{ 1807 MP_RAT c; 1808 mpz_init_set_si(mpq_numref(&c), v1.nint); 1809 mpz_init_set_si(mpq_denref(&c), v2.nint); 1810 mpq_canonicalize(&c); 1811 Make_Rat(pres, TG); 1812 Push_Rat_Mpr(&c); 1813 Succeed_; 1814} 1815 1816 1817static int 1818_big_div(value v1, value v2, pword *pres) /* big x big -> rat */ 1819{ 1820 MP_INT a,b; 1821 Big_To_Mpi(v1.ptr, &a); 1822 Big_To_Mpi(v2.ptr, &b); 1823 if (GlobalFlags & PREFER_RATIONALS) 1824 { 1825 MP_RAT c; 1826 mpz_init_set(mpq_numref(&c), &a); 1827 mpz_init_set(mpq_denref(&c), &b); 1828 mpq_canonicalize(&c); 1829 Pw_From_Mpr(pres, &c); 1830 } 1831 else 1832 { 1833 double d = mpz_fdiv(&a, &b); 1834 Make_Checked_Float(pres, d); 1835 } 1836 Succeed_; 1837} 1838 1839static int 1840_rat_div(value v1, value v2, pword *pres) 1841{ 1842 MP_RAT a,b,c; 1843 if (RatZero(v2.ptr)) 1844 { Bip_Error(ARITH_EXCEPTION); } 1845 mpq_init(&c); 1846 Rat_To_Mpr(v1.ptr, &a); 1847 Rat_To_Mpr(v2.ptr, &b); 1848 mpq_div(&c, &a, &b); 1849 Pw_From_Mpr(pres, &c); 1850 Succeed_; 1851} 1852 1853static int 1854_rat_pow(value v1, value v2, pword *pres) /* rat x int -> rat */ 1855{ 1856 MP_RAT c; 1857 mpq_init(&c); 1858 if (v2.nint == 0) 1859 { 1860 mpq_set_ui(&c, 1L, 1L); 1861 } 1862 else 1863 { 1864 word exp = v2.nint < 0 ? -v2.nint : v2.nint; 1865 MP_RAT a,b; 1866 Rat_To_Mpr(v1.ptr, &b); 1867 mpq_init(&a); 1868 mpq_set(&a, &b); /* need to copy because we assign to it */ 1869 if (exp % 2) 1870 mpq_set(&c, &a); 1871 else 1872 mpq_set_ui(&c, 1L, 1L); 1873 while ((exp /= 2) != 0) 1874 { 1875 mpq_mul(&a, &a, &a); 1876 if (exp % 2) 1877 mpq_mul(&c, &c, &a); 1878 } 1879 mpq_clear(&a); 1880 if (v2.nint < 0) 1881 mpq_inv(&c, &c); 1882 } 1883 Pw_From_Mpr(pres, &c); 1884 Succeed_; 1885} 1886 1887static int 1888_rat_min(value v1, value v2, pword *pres) 1889{ 1890 MP_RAT a,b; 1891 Rat_To_Mpr(v1.ptr, &a); 1892 Rat_To_Mpr(v2.ptr, &b); 1893 Make_Rat(pres, mpq_cmp(&a, &b) < 0 ? v1.ptr : v2.ptr); 1894 Succeed_; 1895} 1896 1897static int 1898_rat_max(value v1, value v2, pword *pres) 1899{ 1900 MP_RAT a,b; 1901 Rat_To_Mpr(v1.ptr, &a); 1902 Rat_To_Mpr(v2.ptr, &b); 1903 Make_Rat(pres, mpq_cmp(&a, &b) > 0 ? v1.ptr : v2.ptr); 1904 Succeed_; 1905} 1906 1907static int 1908_rat_fix(value v1, pword *pres) 1909{ 1910 MP_INT a,b,c; 1911 Big_To_Mpi(Numer(v1.ptr)->val.ptr, &a); 1912 Big_To_Mpi(Denom(v1.ptr)->val.ptr, &b); 1913 mpz_init(&c); 1914 mpz_tdiv_q(&c, &a, &b); 1915 Pw_From_Mpi(pres, &c); 1916 Succeed_; 1917} 1918 1919static int 1920_rat_int2(value v1, pword *pres) 1921{ 1922 if (!RatIntegral(v1.ptr)) 1923 { 1924 Bip_Error(ARITH_EXCEPTION); 1925 } 1926 Normalize_Big(pres, Numer(v1.ptr)->val.ptr); 1927 Succeed_; 1928} 1929 1930static int 1931_rat_f_c(value v1, pword *pres, void (*div_function) (MP_INT*, const MP_INT*, const MP_INT*)) 1932{ 1933 MP_INT a,b,c; 1934 Big_To_Mpi(Denom(v1.ptr)->val.ptr, &b); 1935 if (mpz_cmp_si(&b, 1L) == 0) 1936 { 1937 Make_Rat(pres, v1.ptr); /* it's already integer */ 1938 } 1939 else 1940 { 1941 Big_To_Mpi(Numer(v1.ptr)->val.ptr, &a); 1942 mpz_init(&c); 1943 (*div_function)(&c, &a, &b); 1944 Make_Rat(pres, TG); 1945 Push_Rat_Frame(); 1946 Make_Big(Numer(pres->val.ptr), TG); 1947 Push_Big_Mpi(&c); 1948 Make_Big(Denom(pres->val.ptr), TG); 1949 Push_Big_PosInt(1); 1950 } 1951 Succeed_; 1952} 1953 1954static int 1955_rat_floor(value v1, pword *pres) 1956{ 1957 return _rat_f_c(v1, pres, mpz_fdiv_q); 1958} 1959 1960static int 1961_rat_ceil(value v1, pword *pres) 1962{ 1963 return _rat_f_c(v1, pres, mpz_cdiv_q); 1964} 1965 1966static int 1967_rat_truncate(value v1, pword *pres) 1968{ 1969 return _rat_f_c(v1, pres, mpz_tdiv_q); 1970} 1971 1972static int 1973_rat_num(value v1, pword *pres) 1974{ 1975 Normalize_Big(pres, Numer(v1.ptr)->val.ptr); 1976 Succeed_; 1977} 1978 1979static int 1980_rat_den(value v1, pword *pres) 1981{ 1982 Normalize_Big(pres, Denom(v1.ptr)->val.ptr); 1983 Succeed_; 1984} 1985 1986static int 1987_rat_atan2(value v1, value v2, pword *pres) 1988{ 1989 MP_RAT a,b; 1990 Rat_To_Mpr(v1.ptr, &a); 1991 Rat_To_Mpr(v2.ptr, &b); 1992 Make_Double(pres, Atan2(mpq_to_double(&a), mpq_to_double(&b))); 1993 Succeed_; 1994} 1995 1996 1997/*-------------------------------------------------------------------------- 1998 * Initialize bignums and rationals 1999 *--------------------------------------------------------------------------*/ 2000 2001void 2002bigrat_init(void) 2003{ 2004#ifdef DEBUG_RAT_ALLOC 2005 mp_set_memory_functions(_rat_alloc, _rat_realloc, _rat_free); 2006#else 2007 mp_set_memory_functions((void*(*)(size_t)) hp_alloc_size, 2008 (void*(*)(void*, size_t, size_t)) hp_realloc_size, 2009 (void(*)(void*, size_t)) hp_free_size); 2010#endif 2011 2012 tag_desc[TINT].coerce_to[TBIG] = _int_big; /* 32 bit integers */ 2013 tag_desc[TINT].coerce_to[TRAT] = _int_rat; 2014 tag_desc[TINT].arith_op[ARITH_DIV] = _int_div; /* may yield rat */ 2015 tag_desc[TINT].arith_op[ARITH_CHGSIGN] = 2016 tag_desc[TINT].arith_op[ARITH_NEG] = _int_neg; /* may yield big */ 2017 tag_desc[TINT].arith_op[ARITH_NICERAT] = _int_nicerat; 2018 2019 tag_desc[TDBL].coerce_to[TRAT] = _dbl_rat; /* double */ 2020 tag_desc[TDBL].arith_op[ARITH_FIX] = _dbl_fix; 2021 tag_desc[TDBL].arith_op[ARITH_INT] = _dbl_int2; 2022 tag_desc[TDBL].arith_op[ARITH_NICERAT] = _dbl_nicerat; 2023 2024 tag_desc[TBIG].from_string = _big_from_string; /* bignums */ 2025 tag_desc[TBIG].write = _write_big; 2026 tag_desc[TBIG].compare = _compare_big; 2027 tag_desc[TBIG].arith_compare = _arith_compare_big; 2028 tag_desc[TBIG].equal = _equal_big; 2029 tag_desc[TBIG].copy_size = _copy_size_big; 2030 tag_desc[TBIG].copy_to_heap = _copy_heap_big; 2031 tag_desc[TBIG].string_size = _big_string_size; 2032 tag_desc[TBIG].to_string = _big_to_string; 2033 tag_desc[TBIG].coerce_to[TRAT] = _big_rat; 2034 tag_desc[TBIG].coerce_to[TDBL] = _big_dbl; 2035 tag_desc[TBIG].arith_op[ARITH_PLUS] = _big_nop; 2036 tag_desc[TBIG].arith_op[ARITH_FLOAT] = _big_nop; /* handled by coercion */ 2037 tag_desc[TBIG].arith_op[ARITH_ROUND] = _big_nop; 2038 tag_desc[TBIG].arith_op[ARITH_FLOOR] = _big_nop; 2039 tag_desc[TBIG].arith_op[ARITH_CEIL] = _big_nop; 2040 tag_desc[TBIG].arith_op[ARITH_TRUNCATE] = _big_nop; 2041 tag_desc[TBIG].arith_op[ARITH_FIX] = _big_nop; 2042 tag_desc[TBIG].arith_op[ARITH_INT] = _big_nop; 2043 tag_desc[TBIG].arith_op[ARITH_MIN] = _big_min; 2044 tag_desc[TBIG].arith_op[ARITH_MAX] = _big_max; 2045 tag_desc[TBIG].arith_op[ARITH_ABS] = _big_abs; 2046 tag_desc[TBIG].arith_op[ARITH_CHGSIGN] = 2047 tag_desc[TBIG].arith_op[ARITH_NEG] = _big_neg; 2048 tag_desc[TBIG].arith_op[ARITH_ADD] = _big_add; 2049 tag_desc[TBIG].arith_op[ARITH_SUB] = _big_sub; 2050 tag_desc[TBIG].arith_op[ARITH_MUL] = _big_mul; 2051 tag_desc[TBIG].arith_op[ARITH_DIV] = _big_div; 2052 tag_desc[TBIG].arith_op[ARITH_IDIV] = _big_idiv; 2053 tag_desc[TBIG].arith_op[ARITH_MOD] = _big_rem; 2054 tag_desc[TBIG].arith_op[ARITH_POW] = _big_pow; 2055 tag_desc[TBIG].arith_op[ARITH_FLOORDIV] = _big_floordiv; 2056 tag_desc[TBIG].arith_op[ARITH_FLOORREM] = _big_floorrem; 2057 tag_desc[TBIG].arith_op[ARITH_AND] = _big_and; 2058 tag_desc[TBIG].arith_op[ARITH_OR] = _big_or; 2059 tag_desc[TBIG].arith_op[ARITH_COM] = _big_com; 2060 tag_desc[TBIG].arith_op[ARITH_XOR] = _big_xor; 2061 tag_desc[TBIG].arith_op[ARITH_SHL] = _big_shl; 2062 tag_desc[TBIG].arith_op[ARITH_SHR] = _big_shr; 2063 tag_desc[TBIG].arith_op[ARITH_SGN] = _big_sgn; 2064 tag_desc[TBIG].arith_op[ARITH_SETBIT] = _big_setbit; 2065 tag_desc[TBIG].arith_op[ARITH_GETBIT] = _big_getbit; 2066 tag_desc[TBIG].arith_op[ARITH_CLRBIT] = _big_clrbit; 2067 tag_desc[TBIG].arith_op[ARITH_BOXLONGLONG] = _big_boxlonglong; 2068 tag_desc[TBIG].arith_op[ARITH_TOCLONGLONG] = _big_toclonglong; 2069 tag_desc[TBIG].arith_op[ARITH_NICERAT] = _big_nicerat; 2070 tag_desc[TBIG].arith_op[ARITH_GCD] = _big_gcd; 2071 tag_desc[TBIG].arith_op[ARITH_GCD_EXT] = _big_gcd_ext; 2072 tag_desc[TBIG].arith_op[ARITH_LCM] = _big_lcm; 2073 tag_desc[TBIG].arith_op[ARITH_POWM] = _big_powm; 2074 tag_desc[TBIG].arith_op[ARITH_NEXT] = _big_next; 2075 tag_desc[TBIG].arith_op[ARITH_PREV] = _big_prev; 2076 tag_desc[TBIG].arith_op[ARITH_ATAN2] = _big_atan2; 2077 2078 tag_desc[TRAT].from_string = _rat_from_string; /* rationals */ 2079 tag_desc[TRAT].write = _write_rat; 2080 tag_desc[TRAT].compare = _compare_rat; 2081 tag_desc[TRAT].arith_compare = _arith_compare_rat; 2082 tag_desc[TRAT].equal = _equal_rat; 2083 tag_desc[TRAT].copy_size = _copy_size_rat; 2084 tag_desc[TRAT].copy_to_heap = _copy_heap_rat; 2085 tag_desc[TRAT].string_size = _rat_string_size; 2086 tag_desc[TRAT].to_string = _rat_to_string; 2087 tag_desc[TRAT].coerce_to[TDBL] = _rat_dbl; 2088 tag_desc[TRAT].arith_op[ARITH_PLUS] = _rat_nop; 2089 tag_desc[TRAT].arith_op[ARITH_FLOAT] = _rat_nop; /* handled by coercion */ 2090 tag_desc[TRAT].arith_op[ARITH_NICERAT] = _rat_nop; 2091 tag_desc[TRAT].coerce_to[TIVL] = _rat_ivl; 2092 tag_desc[TRAT].arith_op[ARITH_CHGSIGN] = 2093 tag_desc[TRAT].arith_op[ARITH_NEG] = _rat_neg; 2094 tag_desc[TRAT].arith_op[ARITH_ABS] = _rat_abs; 2095 tag_desc[TRAT].arith_op[ARITH_MIN] = _rat_min; 2096 tag_desc[TRAT].arith_op[ARITH_MAX] = _rat_max; 2097 tag_desc[TRAT].arith_op[ARITH_ADD] = _rat_add; 2098 tag_desc[TRAT].arith_op[ARITH_SUB] = _rat_sub; 2099 tag_desc[TRAT].arith_op[ARITH_MUL] = _rat_mul; 2100 tag_desc[TRAT].arith_op[ARITH_DIV] = _rat_div; 2101 tag_desc[TRAT].arith_op[ARITH_POW] = _rat_pow; 2102 tag_desc[TRAT].arith_op[ARITH_FLOOR] = _rat_floor; 2103 tag_desc[TRAT].arith_op[ARITH_CEIL] = _rat_ceil; 2104 tag_desc[TRAT].arith_op[ARITH_TRUNCATE] = _rat_truncate; 2105 tag_desc[TRAT].arith_op[ARITH_ROUND] = _unimp_err; 2106 tag_desc[TRAT].arith_op[ARITH_FIX] = _rat_fix; 2107 tag_desc[TRAT].arith_op[ARITH_INT] = _rat_int2; 2108 tag_desc[TRAT].arith_op[ARITH_NUM] = _rat_num; 2109 tag_desc[TRAT].arith_op[ARITH_DEN] = _rat_den; 2110 tag_desc[TRAT].arith_op[ARITH_SGN] = _rat_sgn; 2111 tag_desc[TRAT].arith_op[ARITH_ATAN2] = _rat_atan2; 2112} 2113 2114#endif 2115