Back to index

plt-scheme  4.2.1
bignum.c
Go to the documentation of this file.
00001 /*
00002   MzScheme
00003   Copyright (c) 2004-2009 PLT Scheme Inc.
00004   Copyright (c) 1995-2001 Matthew Flatt, Scott Owens
00005 
00006     This library is free software; you can redistribute it and/or
00007     modify it under the terms of the GNU Library General Public
00008     License as published by the Free Software Foundation; either
00009     version 2 of the License, or (at your option) any later version.
00010 
00011     This library is distributed in the hope that it will be useful,
00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014     Library General Public License for more details.
00015 
00016     You should have received a copy of the GNU Library General Public
00017     License along with this library; if not, write to the Free
00018     Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
00019     Boston, MA 02110-1301 USA.
00020 
00021   libscheme
00022   Copyright (c) 1994 Brent Benson
00023   All rights reserved.
00024 */
00025 
00026 
00027 /* DANGER! DANGER! DANGER! DANGER! DANGER! DANGER! DANGER! DANGER!
00028 
00029    This code is fragile, due to the Small_Bignum optimization, and
00030    memory subtleties of bignums.
00031 
00032    When allocating a bignum for a small integer, a Small_Bignum is
00033    allocated. The Small_Bignum structure has room for one bigdig.
00034 
00035    The digit array pointer of a Small_Bignum points into the (middle
00036    of the) Small_Bignum record itself. This means:
00037 
00038      1) For all collectors, the digit array point must not be copied
00039         to another Scheme_Bignum record (because it points into the
00040         middle of the Small_Bignum record, and interior pointers are
00041         not allowed).
00042 
00043      2) Since SCHEME_BIGDIG() can return an interior pointer, for
00044         precise GC the code must be careful about putting
00045         SCHEME_BIGDIG() results into local variables. The
00046         SCHEME_BIGDIG_SAFE() macro copies the digit into a local
00047         array if necessary.
00048 
00049    In addition, the precise GC needs to distinguish Scheme_Bignum from
00050    Small_Bignum for computing sizes; the allocated_inline flag does
00051    that.
00052 
00053    Finally, when pointers are sent into GMP when GMP might block or
00054    allocate, then the pointer needs to be immobile (but it can and
00055    should be GCable, in case a break exception escapes). The PROTECT
00056    macros copy an array as necessary to immobile memory in precise
00057    GC mode.
00058 */
00059 
00060 #include "schpriv.h"
00061 #include <ctype.h>
00062 #include <math.h>
00063 #include "gmp/gmp.h"
00064 
00065 /* Used by gmp: */
00066 void scheme_bignum_use_fuel(long n);
00067 
00068 
00069 #if defined(SIXTY_FOUR_BIT_INTEGERS)
00070 # define FIRST_BIT_MASK 0x8000000000000000
00071 # define SECOND_BIT_MASK 0x4000000000000000
00072 # define MAX_TWO_BIT_MASK 0xC000000000000000
00073 # define SMALL_NUM_STR_LEN 19 /* conservatively low is OK */
00074 # define SQRT_BIT_MAX 31
00075 #else
00076 # define FIRST_BIT_MASK 0x80000000
00077 # define SECOND_BIT_MASK 0x40000000
00078 # define MAX_TWO_BIT_MASK 0xC0000000
00079 # define SMALL_NUM_STR_LEN 10 /* conservatively low is OK */
00080 # define SQRT_BIT_MAX 15
00081 #endif
00082 
00083 #if defined(USE_LONG_LONG_FOR_BIGDIG)
00084 # define TOP_BITS_MASK 0xFFFFFFFF00000000
00085 # define BOTTOM_BITS_MASK 0x00000000FFFFFFFF
00086 #endif
00087 
00088 #if defined(SIXTY_FOUR_BIT_INTEGERS) || defined(USE_LONG_LONG_FOR_BIGDIG)
00089 # define BIG_RADIX 18446744073709551616.0 /* = 0x10000000000000000 */
00090 # define ALL_ONES 0xFFFFFFFFFFFFFFFF
00091 # define WORD_SIZE 64
00092 #else
00093 # define BIG_RADIX 4294967296.0 /* = 0x100000000 */
00094 # define ALL_ONES 0xFFFFFFFF
00095 # define WORD_SIZE 32
00096 #endif
00097 
00098 static Scheme_Object *bignum_one;
00099 
00100 #ifdef MZ_PRECISE_GC
00101 # define SAFE_SPACE(var) bigdig var[1];
00102 # define SCHEME_BIGDIG_SAFE(b, s) ((SCHEME_BIGDIG(b) == ((Small_Bignum *) mzALIAS b)->v) ? (s[0] = SCHEME_BIGDIG(b)[0], s) : SCHEME_BIGDIG(b))
00103 
00104 # define PROTECT(digarray, len) digarray = copy_to_protected(digarray, len * sizeof(bigdig), 0);
00105 # define RELEASE(digarray) (free_protected(digarray), digarray = NULL);
00106 
00107 # define PROTECT_RESULT(len) copy_to_protected(NULL, len * sizeof(bigdig), 1);
00108 # define FINISH_RESULT(digarray, len) { bigdig *save = digarray; digarray = (bigdig *)scheme_malloc_atomic(len * sizeof(bigdig)); memcpy(digarray, save, len * sizeof(bigdig)); RELEASE(save); }
00109 # define MALLOC_PROTECT(size) copy_to_protected(NULL, size, 0)
00110 # define FREE_PROTECT(ptr) free_protected(ptr)
00111 
00112 extern void GC_check(void *p);
00113 
00114 #define BIGNUM_CACHE_SIZE 16
00115 static THREAD_LOCAL void *bignum_cache[BIGNUM_CACHE_SIZE];
00116 static THREAD_LOCAL int cache_count;
00117 
00118 static void *copy_to_protected(void *p, long len, int zero)
00119 {
00120   void *r;
00121   long minsz;
00122   
00123   minsz = GC_malloc_stays_put_threshold();
00124   if (minsz >= len + sizeof(long)) {
00125     if (cache_count) {
00126       --cache_count;
00127       r = bignum_cache[cache_count];
00128       bignum_cache[cache_count] = NULL;
00129     } else
00130       r = (char *)scheme_malloc_atomic(minsz);
00131     ((long *)r)[0] = 1;
00132   } else {
00133     r = (char *)scheme_malloc_atomic(len + sizeof(long));
00134     ((long *)r)[0] = 0;
00135   }
00136 
00137   r = (char *)r XFORM_OK_PLUS sizeof(long);
00138 
00139   if (p) memcpy(r, p, len);
00140   if (zero) memset(r, 0, len);
00141   return r;
00142 }
00143 
00144 static void free_protected(void *p)
00145 {
00146   if (((long *)p)[-1]) {
00147     if (cache_count < BIGNUM_CACHE_SIZE) {
00148       bignum_cache[cache_count++] = (char *)p - sizeof(long);
00149     }
00150   }
00151 }
00152 
00153 void scheme_clear_bignum_cache(void)
00154 {
00155   int i;
00156   for (i = 0; i < BIGNUM_CACHE_SIZE; i++) {
00157     bignum_cache[i] = NULL;
00158   }
00159   cache_count = 0;
00160 }
00161 
00162 #else
00163 # define SAFE_SPACE(var) /*empty */
00164 # define SCHEME_BIGDIG_SAFE(b, s) SCHEME_BIGDIG(b)
00165 
00166 # define PROTECT(digarray, len) /* no-op */
00167 #define RELEASE(digarray) /* no-op */
00168 
00169 # define PROTECT_RESULT(len) allocate_bigdig_array(len)
00170 # define FINISH_RESULT(digarray, len) /* no-op */
00171 
00172 # define MALLOC_PROTECT(size) scheme_malloc_atomic(size)
00173 # define FREE_PROTECT(ptr) /* no-op */
00174 
00175 void scheme_clear_bignum_cache(void) { }
00176 #endif
00177 
00178 #ifdef MZ_XFORM
00179 START_XFORM_SKIP;
00180 #endif
00181 
00182 
00183 #define xor(a, b) (!(a) ^ !(b))
00184 
00185 Scheme_Object *scheme_make_small_bignum(long v, Small_Bignum *o)
00186 {
00187   bigdig bv;
00188 
00189   o->o.iso.so.type = scheme_bignum_type;
00190   SCHEME_SET_BIGPOS(&o->o, ((v >= 0) ? 1 : 0));
00191   if (v < 0)
00192     bv = -v;
00193   else
00194     bv = v;
00195 
00196 #if defined(USE_LONG_LONG_FOR_BIGDIG)
00197   bv = bv & BOTTOM_BITS_MASK;
00198 #endif
00199 
00200   if (bv == 0)
00201     SCHEME_BIGLEN(&o->o) = 0;
00202   else
00203     SCHEME_BIGLEN(&o->o) = 1;
00204 
00205   SCHEME_BIGDIG(&o->o) = o->v;
00206 
00207   o->v[0] = bv;
00208 
00209   return (Scheme_Object *) mzALIAS o;
00210 }
00211 
00212 #ifdef MZ_XFORM
00213 END_XFORM_SKIP;
00214 #endif
00215 
00216 Scheme_Object *scheme_make_bignum(long v)
00217 {
00218   Small_Bignum *r;
00219   r = MALLOC_ONE_TAGGED(Small_Bignum);
00220 #if MZ_PRECISE_GC
00221   SCHEME_SET_BIGINLINE(&r->o);
00222 #endif
00223   return scheme_make_small_bignum(v, r);
00224 }
00225 
00226 Scheme_Object *scheme_make_bignum_from_unsigned(unsigned long v)
00227 {
00228   Small_Bignum *r;
00229   r = MALLOC_ONE_TAGGED(Small_Bignum);
00230 #if MZ_PRECISE_GC
00231   SCHEME_SET_BIGINLINE(&r->o);
00232 #endif
00233   r->o.iso.so.type = scheme_bignum_type;
00234   SCHEME_SET_BIGPOS(&r->o, 1);
00235   if (v == 0)
00236     SCHEME_BIGLEN(&r->o) = 0;
00237   else
00238     SCHEME_BIGLEN(&r->o) = 1;
00239 
00240   SCHEME_BIGDIG(&r->o) = r->v;
00241 
00242   r->v[0] = v;
00243 
00244   return (Scheme_Object*) mzALIAS r;
00245 }
00246 
00247 #ifndef NO_LONG_LONG_TYPE
00248 
00249 Scheme_Object *scheme_make_bignum_from_long_long(mzlonglong v)
00250 {
00251 #if defined(SIXTY_FOUR_BIT_INTEGERS)
00252   return scheme_make_bignum(v);
00253 #else
00254   if (v < 0) {
00255     mzlonglong v2;
00256     
00257     v2 = -v;
00258     if (v2 == v) {
00259       /* This is 0xFFFFFFFFFFFFFFFFLL */
00260       Scheme_Object *o;
00261       bigdig *o_digs;
00262       int len;
00263 #if defined(USE_LONG_LONG_FOR_BIGDIG)
00264       len = 1;
00265 #else
00266       len = 2;
00267 #endif
00268 
00269       o = (Scheme_Object *)scheme_malloc_tagged(sizeof(Scheme_Bignum));      
00270       o->type = scheme_bignum_type;
00271       SCHEME_BIGLEN(o) = len;
00272       SCHEME_SET_BIGPOS(o, 0);
00273       o_digs = (bigdig *)scheme_malloc_atomic(sizeof(bigdig) * len);
00274       SCHEME_BIGDIG(o) = o_digs;
00275 
00276       o_digs[0] = 0;
00277       o_digs[1] = ((bigdig)1 << (WORD_SIZE - 1));
00278       
00279       return (Scheme_Object *)o;      
00280     } else {
00281       Scheme_Object *o;
00282       o = scheme_make_bignum_from_unsigned_long_long((umzlonglong)v2);
00283       SCHEME_SET_BIGPOS(o, 0);
00284       return o;
00285     }
00286   } else {
00287     return scheme_make_bignum_from_unsigned_long_long((umzlonglong)v);
00288   }
00289 #endif
00290 }
00291 
00292 Scheme_Object *scheme_make_bignum_from_unsigned_long_long(umzlonglong v)
00293 {
00294 #if defined(SIXTY_FOUR_BIT_INTEGERS)
00295   return scheme_make_bignum_from_unsigned(v);
00296 #else
00297   int just_one;
00298 
00299 #if defined(USE_LONG_LONG_FOR_BIGDIG)
00300   just_one = 1;
00301 #else
00302   just_one = !((v >> 32) & 0xFFFFFFFF);
00303 #endif
00304 
00305   if (just_one) {
00306     Small_Bignum *r;
00307     r = MALLOC_ONE_TAGGED(Small_Bignum);
00308 #if MZ_PRECISE_GC
00309     SCHEME_SET_BIGINLINE(&r->o);
00310 #endif
00311     r->o.iso.so.type = scheme_bignum_type;
00312     SCHEME_SET_BIGPOS(&r->o, 1);
00313     SCHEME_BIGLEN(&r->o) = 1;
00314     
00315     SCHEME_BIGDIG(&r->o) = r->v;
00316     
00317     r->v[0] = (bigdig)v;
00318 
00319     return (Scheme_Object*) mzALIAS r;
00320   } else {
00321     Scheme_Object *o;
00322     bigdig *o_digs;
00323     
00324     o = (Scheme_Object *)scheme_malloc_tagged(sizeof(Scheme_Bignum));
00325     
00326     o->type = scheme_bignum_type;
00327     SCHEME_BIGLEN(o) = 2;
00328     SCHEME_SET_BIGPOS(o, 1);
00329     o_digs = (bigdig *)scheme_malloc_atomic(sizeof(bigdig) * 2);
00330     SCHEME_BIGDIG(o) = o_digs;
00331 
00332     o_digs[1] = (bigdig)((v >> 32) & 0xFFFFFFFF);
00333     o_digs[0] = (bigdig)(v & 0xFFFFFFFF);
00334     
00335     return (Scheme_Object *)o;
00336   }
00337 #endif
00338 }
00339 
00340 #endif
00341 
00342 /*
00343   Should only succeed if the bignum can fit into a signed long.
00344   This means that the bignum must have length 0 or 1 and the top bit
00345   of its bigdig must be zero, unless it is -100...000.
00346 
00347 */
00348 int scheme_bignum_get_int_val(const Scheme_Object *o, long *v)
00349 {
00350   if (SCHEME_BIGLEN(o) > 1) {    /* won't fit in a signed long */
00351     return 0;
00352   } else if (SCHEME_BIGLEN(o) == 0) {
00353     *v = 0;
00354     return 1;
00355 #ifdef USE_LONG_LONG_FOR_BIGDIG
00356   } else if (SCHEME_BIGDIG(o)[0] & TOP_BITS_MASK) {
00357     return 0;
00358 #endif
00359   } else if (SCHEME_BIGDIG(o)[0] == FIRST_BIT_MASK && !SCHEME_BIGPOS(o)) {
00360     /* Special case for the most negative number representable in a signed word */
00361     *v = SCHEME_BIGDIG(o)[0];
00362     return 1;
00363   } else if ((SCHEME_BIGDIG(o)[0] & FIRST_BIT_MASK) != 0) { /* Won't fit into a signed long */
00364     return 0;
00365   } else if (SCHEME_BIGPOS(o)) {
00366     *v = (long)SCHEME_BIGDIG(o)[0];
00367     return 1;
00368   } else {
00369     *v = -((long)SCHEME_BIGDIG(o)[0]);
00370     return 1;
00371   }
00372 }
00373 
00374 int scheme_bignum_get_unsigned_int_val(const Scheme_Object *o, unsigned long *v)
00375      /* We want to return anything that will fit into a `unsigned long'. */
00376 {
00377   if ((SCHEME_BIGLEN(o) > 1) || !SCHEME_BIGPOS(o))
00378     /* Won't fit into word, or not positive */
00379     return 0;
00380   else if (SCHEME_BIGLEN(o) == 0) {
00381     *v = 0;
00382     return 1;
00383 #ifdef USE_LONG_LONG_FOR_BIGDIG
00384   } else if (SCHEME_BIGDIG(o)[0] & TOP_BITS_MASK) {
00385     return 0;
00386 #endif
00387   } else {
00388     *v = SCHEME_BIGDIG(o)[0];
00389     return 1;
00390   }
00391 }
00392 
00393 #if defined(USE_LONG_LONG_FOR_BIGDIG) || defined(SIXTY_FOUR_BIT_INTEGERS)
00394 # define MAX_BN_SIZE_FOR_LL 1
00395 #else
00396 # define MAX_BN_SIZE_FOR_LL 2
00397 #endif
00398 
00399 int scheme_bignum_get_long_long_val(const Scheme_Object *o, mzlonglong *v)
00400 {
00401 #ifdef NO_LONG_LONG_TYPE
00402   return scheme_bignum_get_int_val(o, v);
00403 #else
00404   if (SCHEME_BIGLEN(o) > MAX_BN_SIZE_FOR_LL) { /* won't fit in a signed long long */
00405     return 0;
00406   } else if (SCHEME_BIGLEN(o) == 0) {
00407     *v = 0;
00408     return 1;
00409   } else if (SCHEME_BIGDIG(o)[MAX_BN_SIZE_FOR_LL - 1] == FIRST_BIT_MASK 
00410 # ifndef USE_LONG_LONG_FOR_BIGDIG
00411             && !SCHEME_BIGDIG(o)[0]
00412 # endif
00413             && !SCHEME_BIGPOS(o)) {
00414     /* Special case for the most negative number representable in a signed long long */
00415     mzlonglong v2;
00416     v2 = 1;
00417     v2 = (v2 << 63);
00418     *v = v2;
00419     return 1;
00420   } else if ((SCHEME_BIGDIG(o)[MAX_BN_SIZE_FOR_LL - 1] & FIRST_BIT_MASK) != 0) { /* Won't fit into a signed long long */
00421     return 0;
00422   } else {
00423     mzlonglong v2;
00424     v2 = SCHEME_BIGDIG(o)[0];
00425     if (SCHEME_BIGLEN(o) > 1) {
00426       v2 |= ((mzlonglong)(SCHEME_BIGDIG(o)[1])) << 32;
00427     }
00428     if (!SCHEME_BIGPOS(o)) {
00429       v2 = -v2;
00430     }
00431     *v = v2;
00432     return 1;
00433   }
00434 #endif
00435 }
00436 
00437 int scheme_bignum_get_unsigned_long_long_val(const Scheme_Object *o, umzlonglong *v)
00438 {
00439 #ifdef NO_LONG_LONG_TYPE
00440   return scheme_bignum_get_unsigned_int_val(o, v);
00441 #else
00442   if ((SCHEME_BIGLEN(o) > MAX_BN_SIZE_FOR_LL) || !SCHEME_BIGPOS(o))
00443     /* Won't fit into word, or not positive */
00444     return 0;
00445   else if (SCHEME_BIGLEN(o) == 0) {
00446     *v = 0;
00447     return 1;
00448   } else {
00449     umzlonglong v2;
00450     v2 = SCHEME_BIGDIG(o)[0];
00451     if (SCHEME_BIGLEN(o)) {
00452       v2 |= ((umzlonglong)SCHEME_BIGDIG(o)[1]) << 32;
00453     }
00454     *v = v2;
00455     return 1;
00456   }
00457 #endif
00458 }
00459 
00460 /* If the bignum fits into a scheme integer, return that instead */
00461 Scheme_Object *scheme_bignum_normalize(const Scheme_Object *o)
00462 {
00463   long v;
00464 
00465   if (!SCHEME_BIGNUMP(o))
00466     return (Scheme_Object *) mzALIAS o;
00467 
00468   if (scheme_bignum_get_int_val(o, &v)) {
00469     long t;
00470 
00471     t = v & MAX_TWO_BIT_MASK;
00472     if (t == 0 || t == MAX_TWO_BIT_MASK)
00473       return scheme_make_integer(v);
00474     else
00475       return (Scheme_Object*) mzALIAS o;
00476   } else
00477     return (Scheme_Object*) mzALIAS o;
00478 }
00479 
00480 static Scheme_Object *make_single_bigdig_result(int pos, bigdig d)
00481 {
00482   Small_Bignum *sm, quick;
00483   Scheme_Object *o;
00484 
00485   /* May not need to allocate: */
00486   sm = &quick;
00487   sm->o.iso.so.type = scheme_bignum_type;
00488   SCHEME_SET_BIGPOS(sm, pos);
00489   SCHEME_BIGLEN(sm) = 1;
00490   SCHEME_BIGDIG(sm) = sm->v;
00491   sm->v[0] = d;
00492 
00493   o = scheme_bignum_normalize((Scheme_Object *) mzALIAS sm);
00494   if (SAME_OBJ(o, (Scheme_Object *) mzALIAS sm)) {
00495     sm = MALLOC_ONE_TAGGED(Small_Bignum);
00496     sm->o.iso.so.type = scheme_bignum_type;
00497 #if MZ_PRECISE_GC
00498     SCHEME_SET_BIGINLINE(sm);
00499 #endif
00500     SCHEME_SET_BIGPOS(sm, pos);
00501     SCHEME_BIGLEN(sm) = 1;
00502     SCHEME_BIGDIG(sm) = sm->v;
00503     sm->v[0] = d;
00504     return (Scheme_Object *) mzALIAS sm;
00505   } else
00506     return o;
00507 }
00508 
00509 /*
00510    copy the bignum a, and if msd != 0, concat. it as the most significant
00511    digit
00512 */
00513 static Scheme_Object *bignum_copy(const Scheme_Object *a, long msd)
00514 {
00515   Scheme_Object* o;
00516   int c;
00517   bigdig* o_digs;
00518 
00519   c = SCHEME_BIGLEN(a);
00520   o = (Scheme_Object *)scheme_malloc_tagged(sizeof(Scheme_Bignum));
00521 
00522   o->type = scheme_bignum_type;
00523   SCHEME_BIGLEN(o) = c;
00524   SCHEME_SET_BIGPOS(o, SCHEME_BIGPOS(a));
00525   o_digs = (bigdig *)scheme_malloc_atomic(sizeof(bigdig) * (c + (msd ? 1 : 0)));
00526   SCHEME_BIGDIG(o) = o_digs;
00527 
00528   memcpy(o_digs, SCHEME_BIGDIG(a), sizeof(bigdig) * c);
00529 
00530   if (msd) {
00531     o_digs[c] = msd;
00532     SCHEME_BIGLEN(o) = SCHEME_BIGLEN(o) + 1;
00533   }
00534   return o;
00535 }
00536 
00537 int scheme_bignum_eq(const Scheme_Object *a, const Scheme_Object *b)
00538 {
00539   long a_len, b_len;
00540 
00541   a_len = SCHEME_BIGLEN(a);
00542   b_len = SCHEME_BIGLEN(b);
00543 
00544   if (a_len == 0 && b_len == 0)
00545     return 1;
00546 
00547   if (a_len == b_len && SCHEME_BIGPOS(a) == SCHEME_BIGPOS(b))
00548     /* mpn_cmp doesn't allocate or block: */
00549     return mpn_cmp(SCHEME_BIGDIG(a), SCHEME_BIGDIG(b), b_len) == 0;
00550   else
00551     return 0;
00552 }
00553 
00554 /* - if a < b, 0 if a == b, + if  a > b */
00555 XFORM_NONGCING static int bignum_abs_cmp(const Scheme_Object *a, const Scheme_Object *b)
00556 {
00557   long a_len, b_len;
00558 
00559   a_len = SCHEME_BIGLEN(a);
00560   b_len = SCHEME_BIGLEN(b);
00561 
00562   if (a_len > b_len)
00563     return 1;
00564   else if (a_len < b_len)
00565     return -1;
00566   else if (a_len == 0)
00567     return 0;
00568   else
00569     /* mpn_cmp doesn't allocate or block: */
00570     return mpn_cmp(SCHEME_BIGDIG(a), SCHEME_BIGDIG(b), b_len);
00571 }
00572 
00573 int scheme_bignum_lt(const Scheme_Object *a, const Scheme_Object *b)
00574 {
00575   long a_pos, b_pos;
00576   int res;
00577 
00578   a_pos = SCHEME_BIGPOS(a);
00579   b_pos = SCHEME_BIGPOS(b);
00580 
00581   if (!a_pos && b_pos)
00582     return 1;
00583   else if (a_pos && !b_pos)
00584     return 0;
00585   else
00586     res = bignum_abs_cmp(a, b);
00587   if (!a_pos)
00588     return (res > 0);
00589   else
00590     return (res < 0);
00591 }
00592 
00593 int scheme_bignum_gt(const Scheme_Object *a, const Scheme_Object *b)
00594 {
00595   return scheme_bignum_lt(b, a);
00596 }
00597 
00598 int scheme_bignum_le(const Scheme_Object *a, const Scheme_Object *b)
00599 {
00600   return !scheme_bignum_gt(a, b);
00601 }
00602 
00603 int scheme_bignum_ge(const Scheme_Object *a, const Scheme_Object *b)
00604 {
00605   return !scheme_bignum_lt(a, b);
00606 }
00607 
00608 Scheme_Object *scheme_bignum_negate(const Scheme_Object *n)
00609 {
00610   Scheme_Object *o;
00611   int len;
00612 
00613   len = SCHEME_BIGLEN(n);
00614 
00615   if (SCHEME_BIGDIG(n) == ((Small_Bignum *) mzALIAS n)->v) {
00616     /* Can't share bigdig array when n is a Small_Bignum */
00617     o = (Scheme_Object *)scheme_malloc_tagged(sizeof(Small_Bignum));
00618 #if MZ_PRECISE_GC
00619     SCHEME_SET_BIGINLINE(o);
00620 #endif
00621     ((Small_Bignum *)o)->v[0] = SCHEME_BIGDIG(n)[0];
00622     SCHEME_BIGDIG(o) = ((Small_Bignum *) mzALIAS o)->v;
00623   } else {
00624     o = (Scheme_Object *)MALLOC_ONE_TAGGED(Scheme_Bignum);
00625     SCHEME_BIGDIG(o) = SCHEME_BIGDIG(n);
00626   }
00627 
00628   o->type = scheme_bignum_type;
00629   SCHEME_SET_BIGPOS(o, !SCHEME_BIGPOS(n));
00630   SCHEME_BIGLEN(o) = len;
00631 
00632   return o;
00633 }
00634 
00635 static bigdig* allocate_bigdig_array(int length)
00636 {
00637   int i;
00638   bigdig* res;
00639   if (length > 4096) {
00640     res = (bigdig *)scheme_malloc_fail_ok(scheme_malloc_atomic, length * sizeof(bigdig));
00641   } else {
00642     res = (bigdig *)scheme_malloc_atomic(length * sizeof(bigdig));
00643   }
00644   for(i = 0; i < length; ++i) {
00645     res[i] = 0;
00646   }
00647   return res;
00648 }
00649 
00650 /* We don't want to count leading digits of 0 in the bignum's length */
00651 XFORM_NONGCING static int bigdig_length(bigdig* array, int alloced)
00652 {
00653   alloced--;
00654   while (alloced >= 0 && array[alloced] == 0) {
00655     alloced--;
00656   }
00657   return alloced + 1;
00658 }
00659 
00660 /* if (sub) a - b else a + b */
00661 static Scheme_Object *bignum_add_sub(const Scheme_Object *a, const Scheme_Object *b, int sub)
00662 {
00663   Scheme_Object *o;
00664   long a_size, b_size, max_size;
00665   short a_pos, b_pos;
00666 
00667   bigdig *o_digs, *a_digs, *b_digs;
00668   SAFE_SPACE(asd) SAFE_SPACE(bsd)
00669 
00670   a_size = SCHEME_BIGLEN(a);
00671   b_size = SCHEME_BIGLEN(b);
00672   a_pos = SCHEME_BIGPOS(a);
00673   b_pos = xor(SCHEME_BIGPOS(b), sub);
00674   a_digs = SCHEME_BIGDIG_SAFE(a, asd);
00675   b_digs = SCHEME_BIGDIG_SAFE(b, bsd);
00676 
00677   if (b_size == 0)
00678     return scheme_bignum_normalize(bignum_copy(a, 0));
00679   else if (a_size == 0) {
00680     o = bignum_copy(b, 0);
00681     SCHEME_SET_BIGPOS(o, b_pos);
00682     return scheme_bignum_normalize(o);
00683   }
00684 
00685   o = (Scheme_Object *)scheme_malloc_tagged(sizeof(Scheme_Bignum));
00686   o->type = scheme_bignum_type;
00687 
00688   o_digs = NULL; /* Get rid of erroneous gcc warning */
00689 
00690   max_size = (a_size > b_size) ? a_size : b_size;
00691 
00692   if (a_pos == b_pos) /* addition */
00693   {
00694     int carry;
00695 
00696     o_digs = allocate_bigdig_array(max_size);
00697 
00698     /* mpn_add doesn't allocate or block */
00699     if (a_size > b_size)
00700       carry = mpn_add(o_digs, a_digs, a_size, b_digs, b_size);
00701     else
00702       carry = mpn_add(o_digs, b_digs, b_size, a_digs, a_size);
00703 
00704     SCHEME_SET_BIGPOS(o, a_pos);
00705     SCHEME_BIGLEN(o) = max_size;
00706     SCHEME_BIGDIG(o) = o_digs;
00707     if (carry)
00708       o = bignum_copy(o, 1);
00709   }
00710   else /* subtraction */
00711   {
00712     int sw;
00713     if (a_size > b_size)
00714       sw = 0;
00715     else if (b_size > a_size)
00716       sw = 1;
00717     else
00718     {
00719       int cmp;
00720       cmp = mpn_cmp(a_digs, b_digs, a_size); /* doesn't allocate or block */
00721       if (cmp == 0)
00722        return scheme_make_integer(0);
00723       else if (cmp > 0) /* a > b */
00724        sw = 0;
00725       else
00726        sw = 1;
00727     }
00728     o_digs = allocate_bigdig_array(max_size);
00729 
00730     /* mpn_sub doesn't allocate or block */
00731     if (sw)
00732       mpn_sub(o_digs, b_digs, b_size, a_digs, a_size);
00733     else
00734       mpn_sub(o_digs, a_digs, a_size, b_digs, b_size);
00735 
00736     SCHEME_SET_BIGPOS(o, xor(sw, a_pos));
00737     max_size = bigdig_length(o_digs, max_size);
00738     SCHEME_BIGLEN(o) = max_size;
00739     SCHEME_BIGDIG(o) = o_digs;
00740   }
00741   return scheme_bignum_normalize(o);
00742 }
00743 
00744 Scheme_Object *scheme_bignum_add(const Scheme_Object *a, const Scheme_Object *b)
00745 {
00746   return bignum_add_sub(a, b, 0);
00747 }
00748 
00749 Scheme_Object *scheme_bignum_subtract(const Scheme_Object *a, const Scheme_Object *b)
00750 {
00751   return bignum_add_sub(a, b, 1);
00752 }
00753 
00754 Scheme_Object *scheme_bignum_add1(const Scheme_Object *n)
00755 {
00756   if (!bignum_one) {
00757     REGISTER_SO(bignum_one);
00758     bignum_one = scheme_make_bignum(1);
00759   }
00760 
00761   return bignum_add_sub(n, bignum_one, 0);
00762 }
00763 
00764 Scheme_Object *scheme_bignum_sub1(const Scheme_Object *n)
00765 {
00766   if (!bignum_one) {
00767     REGISTER_SO(bignum_one);
00768     bignum_one = scheme_make_bignum(1);
00769   }
00770 
00771   return bignum_add_sub(n, bignum_one, 1);
00772 }
00773 
00774 /* norm determines if we normalize the result */
00775 static Scheme_Object *bignum_multiply(const Scheme_Object *a, const Scheme_Object *b, int norm)
00776 {
00777   Scheme_Object *o;
00778   long a_size, a_pos, b_size, b_pos, res_size, i, j;
00779   bigdig* o_digs, *a_digs, *b_digs;
00780   SAFE_SPACE(asd) SAFE_SPACE(bsd)
00781 
00782   a_size = SCHEME_BIGLEN(a);
00783   b_size = SCHEME_BIGLEN(b);
00784 
00785   SCHEME_USE_FUEL(a_size);
00786   SCHEME_USE_FUEL(b_size);
00787 
00788   if (a_size == 0 || b_size == 0)
00789   {
00790     if (norm)
00791       return scheme_make_integer(0);
00792     else
00793       return scheme_make_bignum(0);
00794   }
00795 
00796   a_pos = SCHEME_BIGPOS(a);
00797   b_pos = SCHEME_BIGPOS(b);
00798   a_digs = SCHEME_BIGDIG_SAFE(a, asd);
00799   b_digs = SCHEME_BIGDIG_SAFE(b, bsd);
00800 
00801   res_size = a_size + b_size;
00802 
00803   o = (Scheme_Object *)scheme_malloc_tagged(sizeof(Scheme_Bignum));
00804   o->type = scheme_bignum_type;
00805 
00806   o_digs = PROTECT_RESULT(res_size);
00807 
00808   PROTECT(a_digs, a_size);
00809   PROTECT(b_digs, b_size);
00810 
00811   for (i = 0; (a_digs[i] == 0) && i < a_size; i++) {
00812     o_digs[i] = 0;
00813   }
00814   for (j = 0; (b_digs[j] == 0) && j < b_size; j++) {
00815     o_digs[i + j] = 0;
00816   }
00817 
00818   if ((a_size - i) > (b_size - j))
00819     mpn_mul(o_digs XFORM_OK_PLUS i + j, a_digs XFORM_OK_PLUS i, a_size - i, b_digs XFORM_OK_PLUS j, b_size - j);
00820   else
00821     mpn_mul(o_digs XFORM_OK_PLUS i + j, b_digs XFORM_OK_PLUS j, b_size - j, a_digs XFORM_OK_PLUS i, a_size - i);
00822 
00823   RELEASE(a_digs);
00824   RELEASE(b_digs);
00825 
00826   FINISH_RESULT(o_digs, res_size);
00827 
00828   res_size = bigdig_length(o_digs, res_size);
00829   SCHEME_BIGLEN(o) = res_size;
00830   SCHEME_BIGDIG(o) = o_digs;
00831   SCHEME_SET_BIGPOS(o, !xor(a_pos, b_pos));
00832 
00833   return (norm ? scheme_bignum_normalize(o) : o);
00834 }
00835 
00836 Scheme_Object *scheme_bignum_multiply(const Scheme_Object *a, const Scheme_Object *b)
00837 {
00838   return bignum_multiply(a, b, 1);
00839 }
00840 
00841 static Scheme_Object *do_power(const Scheme_Object *a, unsigned long b)
00842 {
00843   Scheme_Object *result;
00844   int i;
00845 
00846   result = scheme_make_integer(1);
00847 
00848   i = sizeof(unsigned long) * 8- 1;
00849   while (!((b >> i) & 0x1) && i >= 0)
00850   {
00851     i = i - 1;
00852   }
00853 
00854   while (i >= 0)
00855   {
00856     result = scheme_bin_mult(result, result);
00857     if ((b >> i) & 0x1)
00858       result = scheme_bin_mult(a, result);
00859     i = i - 1;
00860   }
00861   return result;
00862 }
00863 
00864 Scheme_Object *do_big_power(const Scheme_Object *a, const Scheme_Object *b)
00865 {
00866   /* This is really a fancy way of sleeping, because it's only used
00867      when b is a bignum, which means that we have no chance of actually
00868      reaching the result. But just in case... */
00869   Scheme_Object *result, *v[2];
00870 
00871   result = scheme_make_integer(1);
00872   v[1] = scheme_make_integer(-1);
00873 
00874   while (!scheme_is_zero(b)) {
00875     if (SCHEME_TRUEP(scheme_odd_p(1, (Scheme_Object **)&b)))
00876       result = scheme_bin_mult(a, result);
00877     a = scheme_bin_mult(a, a);
00878 
00879     v[0] = (Scheme_Object *)b;
00880     b = scheme_bitwise_shift(2, v);
00881   }
00882 
00883   return result;
00884 }
00885 
00886 
00887 Scheme_Object *scheme_generic_integer_power(const Scheme_Object *a, const Scheme_Object *b)
00888 {
00889   unsigned long exponent;
00890 
00891   if (scheme_current_thread->constant_folding) {
00892     /* if we're trying to fold a constant, limit the work that we're willing to do at compile time */
00893     GC_CAN_IGNORE const char *too_big = "arguments too big to fold `expt'";
00894     if (SCHEME_BIGNUMP(b)
00895         || (SCHEME_INT_VAL(b) > 10000))
00896       scheme_signal_error(too_big);
00897     else if (SCHEME_BIGNUMP(a)) {
00898       int len = SCHEME_BIGLEN(a);
00899       if ((len > 10000)
00900           || (len * SCHEME_INT_VAL(b)) > 10000)
00901         scheme_signal_error(too_big);
00902     }
00903   }
00904 
00905   if (scheme_get_unsigned_int_val((Scheme_Object *)b, &exponent))
00906     return do_power(a, exponent);
00907   else
00908     return do_big_power(a, b);
00909 }
00910 
00911 Scheme_Object *scheme_bignum_max(const Scheme_Object *a, const Scheme_Object *b)
00912 {
00913   int lt;
00914   lt = scheme_bignum_lt(a, b);
00915   return scheme_bignum_normalize(lt ? b : a);
00916 }
00917 
00918 Scheme_Object *scheme_bignum_min(const Scheme_Object *a, const Scheme_Object *b)
00919 {
00920   int lt;
00921   lt = scheme_bignum_lt(a, b);
00922   return scheme_bignum_normalize(lt ? a : b);
00923 }
00924 
00925 /* op = 0 : &
00926    op = 1 : |
00927    op = 2 : ^
00928 assumes len a >= len b */
00929 static Scheme_Object *do_bitop(const Scheme_Object *a, const Scheme_Object *b, int op)
00930 {
00931   long a_size, b_size, a_pos, b_pos, res_alloc, i;
00932   short res_pos;
00933   bigdig* a_digs, *b_digs, *res_digs, quick_digs[1];
00934   int carry_out_a, carry_out_b, carry_out_res, carry_in_a, carry_in_b, carry_in_res;
00935   Scheme_Object* o;
00936   SAFE_SPACE(asd) SAFE_SPACE(bsd)
00937 
00938   a_size = SCHEME_BIGLEN(a);
00939   b_size = SCHEME_BIGLEN(b);
00940 
00941   if (a_size == 0) /* b_size == 0 too */
00942   {
00943     return scheme_make_integer(0); /* for all 3 ops */
00944   }
00945   else if (b_size == 0)
00946   {
00947     if (op == 0)
00948       return scheme_make_integer(0);
00949     else
00950       return scheme_bignum_normalize(bignum_copy(a, 0));
00951   }
00952 
00953   a_pos = SCHEME_BIGPOS(a);
00954   a_digs = SCHEME_BIGDIG_SAFE(a, asd);
00955   b_pos = SCHEME_BIGPOS(b);
00956   b_digs = SCHEME_BIGDIG_SAFE(b, bsd);
00957 
00958   if (op == 0)
00959   {
00960     res_pos = a_pos || b_pos;
00961     res_alloc = (b_pos ? b_size : a_size);
00962   }
00963   else if (op == 1)
00964   {
00965     res_pos = a_pos && b_pos;
00966     res_alloc = (b_pos ? a_size : b_size);
00967   }
00968   else
00969   {
00970     res_pos = !xor(a_pos, b_pos);
00971     res_alloc = a_size;
00972   }
00973 
00974   if (res_alloc < 2)
00975     res_digs = quick_digs;
00976   else
00977     res_digs = allocate_bigdig_array(res_alloc);
00978 
00979   carry_out_a = carry_out_b = carry_out_res = 1;
00980   carry_in_a = carry_in_b = carry_in_res = 0;
00981 
00982   for (i = 0; i < res_alloc; ++i)
00983   {
00984     bigdig a_val, b_val, res_val;
00985 
00986     a_val = a_digs[i];
00987     if (!a_pos)
00988     {
00989       /* We have to do te operation on the 2's complement of a */
00990       carry_in_a = carry_out_a;
00991       carry_out_a = (carry_in_a == 1 && a_val == 0) ? 1 : 0;
00992       a_val = ~a_val + carry_in_a;
00993     }
00994 
00995     if (i < b_size)
00996     {
00997       b_val = b_digs[i];
00998       if (!b_pos)
00999       {
01000        carry_in_b = carry_out_b;
01001        carry_out_b = (carry_in_b == 1 && b_val == 0) ? 1 : 0;
01002        b_val = ~b_val + carry_in_b;
01003       }
01004     }
01005     else
01006     {
01007       if (b_pos)
01008        b_val = 0;
01009       else
01010        b_val = ALL_ONES;
01011     }
01012 
01013     if (op == 0)
01014       res_val = a_val & b_val;
01015     else if (op == 1)
01016       res_val = a_val | b_val;
01017     else
01018       res_val = a_val ^ b_val;
01019 
01020     if (!res_pos)
01021     {
01022       carry_in_res = carry_out_res;
01023       carry_out_res = (carry_in_res == 1 && res_val == 0) ? 1 : 0;
01024       res_val = ~res_val + carry_in_res;
01025     }
01026 
01027     res_digs[i] = res_val;
01028   }
01029 
01030   if (!res_pos && carry_out_res == 1) {
01031     /* Overflow => we need an extra digit */
01032     res_digs = allocate_bigdig_array(res_alloc + 1);
01033     for (i = 0; i < res_alloc; i++) {
01034       res_digs[i] = 0;
01035     }
01036     res_digs[res_alloc] = 1;
01037     res_alloc = res_alloc + 1;
01038   } else {
01039     res_alloc = bigdig_length(res_digs, res_alloc);
01040   }
01041 
01042   if (!res_alloc) {
01043     return scheme_make_integer(0);
01044   } else if (res_alloc == 1) {
01045     return make_single_bigdig_result(res_pos, res_digs[0]);
01046   } else {
01047     o = (Scheme_Object*)scheme_malloc_tagged(sizeof(Scheme_Bignum));
01048     o->type = scheme_bignum_type;
01049     SCHEME_SET_BIGPOS(o, res_pos);
01050     SCHEME_BIGLEN(o) = res_alloc;
01051     SCHEME_BIGDIG(o) = res_digs;
01052 
01053     return o;
01054   }
01055 }
01056 
01057 Scheme_Object *scheme_bignum_and(const Scheme_Object *a, const Scheme_Object *b)
01058 {
01059   if (SCHEME_BIGLEN(a) > SCHEME_BIGLEN(b))
01060     return do_bitop(a, b, 0);
01061   else
01062     return do_bitop(b, a, 0);
01063 }
01064 
01065 Scheme_Object *scheme_bignum_or(const Scheme_Object *a, const Scheme_Object *b)
01066 {
01067   if (SCHEME_BIGLEN(a) > SCHEME_BIGLEN(b))
01068     return do_bitop(a, b, 1);
01069   else
01070     return do_bitop(b, a, 1);
01071 }
01072 
01073 Scheme_Object *scheme_bignum_xor(const Scheme_Object *a, const Scheme_Object *b)
01074 {
01075    if (SCHEME_BIGLEN(a) > SCHEME_BIGLEN(b))
01076     return do_bitop(a, b, 2);
01077   else
01078     return do_bitop(b, a, 2);
01079 }
01080 
01081 Scheme_Object *scheme_bignum_not(const Scheme_Object *a)
01082 {
01083   Scheme_Object *o;
01084 
01085   o = scheme_bignum_add1(a);
01086 
01087   if (SCHEME_BIGNUMP(o)) {
01088     SCHEME_SET_BIGPOS(o, !SCHEME_BIGPOS(o));
01089     return scheme_bignum_normalize(o);
01090   } else {
01091     return scheme_bin_minus(scheme_make_integer(0), o);
01092   }
01093 }
01094 
01095 Scheme_Object *scheme_bignum_shift(const Scheme_Object *n, long shift)
01096 {
01097   Scheme_Object* o;
01098   bigdig* res_digs, *n_digs, quick_digs[1], shift_out;
01099   long res_alloc, shift_words, shift_bits, i, j, n_size;
01100   SAFE_SPACE(nsd)
01101 
01102   n_size = SCHEME_BIGLEN(n);
01103   if (n_size == 0)
01104     return scheme_make_integer(0);
01105   if (shift == 0) /* no shift */
01106     return scheme_bignum_normalize(bignum_copy(n, 0));
01107 
01108   n_digs = SCHEME_BIGDIG_SAFE(n, nsd);
01109 
01110   if (shift < 0) /* right shift */
01111   {
01112     int shifted_off_one = 0;
01113 
01114     shift = -shift;
01115     shift_words = shift / WORD_SIZE;
01116     shift_bits = shift % WORD_SIZE;
01117 
01118     if (shift_words >= n_size) {
01119       if (SCHEME_BIGPOS(n))
01120        return scheme_make_integer(0);
01121       else
01122        return scheme_make_integer(-1);
01123     }
01124 
01125     res_alloc = n_size - shift_words;
01126     if (shift_bits == 0 && !SCHEME_BIGPOS(n))
01127       res_alloc++;   /* Very unlikely event of a carryout on the later add1 increasing the word size */
01128     if (res_alloc < 2)
01129       res_digs = quick_digs;
01130     else
01131       res_digs = allocate_bigdig_array(res_alloc);
01132 
01133     if (!SCHEME_BIGPOS(n)) {
01134       for(i = 0; i < shift_words; ++i) {
01135        if (n_digs[i] != 0) {
01136          shifted_off_one = 1;
01137          break;
01138        }
01139       }
01140     }
01141 
01142     for(i = 0, j = shift_words; j < n_size; ++i, ++j) {
01143       res_digs[i] = n_digs[j];
01144     }
01145 
01146     if (shift_bits)
01147       shift_out = mpn_rshift(res_digs, res_digs, res_alloc, shift_bits); /* no allocation/blocking */
01148     else
01149       shift_out = 0;
01150 
01151     if (!SCHEME_BIGPOS(n) && (shifted_off_one || shift_out)) {
01152       mpn_add_1(res_digs, res_digs, res_alloc, 1); /* no allocation/blocking */
01153     }
01154   }
01155   else /* left shift */
01156   {
01157     shift_words = shift / WORD_SIZE;
01158     shift_bits = shift % WORD_SIZE;
01159     res_alloc = SCHEME_BIGLEN(n) + shift_words;
01160     if (shift_bits != 0)
01161       ++res_alloc;
01162     if (res_alloc < 2)
01163       res_digs = quick_digs;
01164     else
01165       res_digs = allocate_bigdig_array(res_alloc);
01166 
01167     for (i = 0, j = shift_words; i < SCHEME_BIGLEN(n); ++i, ++j) {
01168       res_digs[j] = n_digs[i];
01169     }
01170 
01171     if (shift_bits != 0)
01172       /* no allocation/blocking */
01173       mpn_lshift(res_digs XFORM_OK_PLUS shift_words, res_digs XFORM_OK_PLUS shift_words, res_alloc - shift_words, shift_bits);
01174 
01175   }
01176 
01177   res_alloc = bigdig_length(res_digs, res_alloc);
01178 
01179   if (res_alloc == 0) {
01180     return scheme_make_integer(0);
01181   } else if (res_alloc == 1) {
01182     return make_single_bigdig_result(SCHEME_BIGPOS(n), res_digs[0]);
01183   } else {
01184     o = (Scheme_Object *)scheme_malloc_tagged(sizeof(Scheme_Bignum));
01185     o->type = scheme_bignum_type;
01186     SCHEME_BIGDIG(o) = res_digs;
01187     SCHEME_BIGLEN(o) = res_alloc;
01188     SCHEME_SET_BIGPOS(o, SCHEME_BIGPOS(n));
01189     return scheme_bignum_normalize(o);
01190   }
01191 }
01192 
01193 
01194 char *scheme_bignum_to_allocated_string(const Scheme_Object *b, int radix, int alloc)
01195 {
01196   Scheme_Object *c;
01197   unsigned char* str, *str2;
01198   int i, slen, start, clen;
01199   bigdig *c_digs;
01200   SAFE_SPACE(csd)
01201 
01202   if (radix != 10 && radix != 2 && radix != 8 && radix != 16)
01203     scheme_raise_exn(MZEXN_FAIL_CONTRACT, "bad bignum radix: %d", radix);
01204 
01205   if (SCHEME_BIGLEN(b) == 0) {
01206     if (alloc) {
01207       str2 = (unsigned char *)scheme_malloc_atomic(2);
01208       str2[0] = '0';
01209       str2[1] = 0;
01210       return (char *)str2;
01211     } else
01212       return "0";
01213   }
01214 
01215   c = bignum_copy(b, 1);  /* mpn_get_string may need a word of scratch space */
01216 
01217   if (radix == 2)
01218     slen = WORD_SIZE * SCHEME_BIGLEN(b) + 2;
01219   else if (radix == 8)
01220     slen = (int)(ceil(WORD_SIZE * SCHEME_BIGLEN(b) / 3.0) + 2);
01221   else if (radix == 16)
01222     slen = WORD_SIZE * SCHEME_BIGLEN(b) / 4 + 2;
01223   else /* (radix == 10) */
01224     slen = (int)(ceil(WORD_SIZE * SCHEME_BIGLEN(b) * 0.30102999566398115)) + 1;
01225 
01226   str = (unsigned char *)MALLOC_PROTECT(slen);
01227 
01228   c_digs = SCHEME_BIGDIG_SAFE(c, csd);
01229   clen = SCHEME_BIGLEN(c);
01230   PROTECT(c_digs, clen);
01231 
01232   slen = mpn_get_str(str, radix, c_digs, SCHEME_BIGLEN(c) - 1);
01233 
01234   RELEASE(c_digs);
01235 
01236 #ifdef MZ_PRECISE_GC
01237   {
01238     unsigned char *save = str;
01239     str = (unsigned char*)scheme_malloc_atomic(slen);
01240     memcpy(str, save, slen);
01241     FREE_PROTECT(save);
01242   }
01243 #endif
01244 
01245   i = 0;
01246   while (i < slen && str[i] == 0) {
01247     ++i;
01248   }
01249 
01250   if (i == slen) {
01251     if (alloc) {
01252       str2 = (unsigned char *)scheme_malloc_atomic(2);
01253       str2[0] = '0';
01254       str2[1] = 0;
01255       return (char *)str2;
01256     } else
01257       return "0";
01258   } else
01259     slen = slen - i + 1 + (SCHEME_BIGPOS(b) ? 0 : 1);
01260 
01261   str2 = (unsigned char *)scheme_malloc_atomic(slen);
01262 
01263   start = i;
01264 
01265   if (!(SCHEME_BIGPOS(b))) {
01266     i = 1;
01267     start--;
01268     str2[0] = '-';
01269   } else
01270     i = 0;
01271 
01272   for (; i < slen - 1; ++i) {
01273     if (str[i + start] < 10)
01274       str2[i] = str[i + start] + '0';
01275     else
01276       str2[i] = str[i + start] + 'a' - 10;
01277   }
01278 
01279   str2[slen - 1] = 0;
01280 
01281   return (char *)str2;
01282 }
01283 
01284 char *scheme_bignum_to_string(const Scheme_Object *b, int radix)
01285 {
01286   return scheme_bignum_to_allocated_string(b, radix, 0);
01287 }
01288 
01289 Scheme_Object *scheme_read_bignum(const mzchar *str, int offset, int radix)
01290 {
01291   int len, negate, stri, alloc, i, test;
01292   Scheme_Object* o;
01293   bigdig* digs;
01294   unsigned char* istring;
01295 
01296   if (radix < 0 || radix > 16) {
01297     return scheme_false;
01298   }
01299 
01300   negate = 0;
01301   stri = offset;
01302   while ((str[stri] == '+') || (str[stri] == '-')) {
01303     if (str[stri] == '-')
01304       negate = !negate;
01305     stri++;
01306   }
01307   len = scheme_char_strlen(str XFORM_OK_PLUS stri);
01308 
01309   if (radix == 10 && (len < SMALL_NUM_STR_LEN)) {
01310     /* try simple fixnum read first */
01311     long fx;
01312     if (!str[stri])
01313       return scheme_false;
01314     for (fx = 0; str[stri]; stri++) {
01315       if (str[stri] < '0' || str[stri] > '9')
01316        return scheme_false;
01317       fx = (fx * 10) + (str[stri] - '0');
01318     }
01319     if (negate)
01320        fx = -fx;
01321     return scheme_make_integer(fx);
01322   }
01323 
01324   /* Convert string of chars to string of bytes: */
01325 
01326   istring = (unsigned char *)MALLOC_PROTECT(len);
01327 
01328   i = stri;
01329   while(str[i] != 0) {
01330     if (str[i] >= '0' && str[i] <= '9')
01331       istring[i - stri] = str[i] - '0';
01332     else if (str[i] >= 'a' && str[i] <= 'z')
01333       istring[i - stri] = str[i] - 'a' + 10;
01334     else if (str[i] >= 'A' && str[i] <= 'Z')
01335       istring[i - stri] = str[i] - 'A' + 10;
01336     else
01337       return scheme_false;
01338 
01339     if (istring[i - stri] >= radix)
01340       return scheme_false;
01341     ++i;
01342   }
01343 
01344   o = (Scheme_Object *)scheme_malloc_tagged(sizeof(Scheme_Bignum));
01345   o->type = scheme_bignum_type;
01346 
01347   alloc = (int)(ceil(len * log((double)radix) / (32 * log((double)2))));
01348 
01349   digs = PROTECT_RESULT(alloc);
01350 
01351   SCHEME_SET_BIGPOS(o, !negate);
01352 
01353   test = mpn_set_str(digs, istring, len, radix);
01354 
01355   FREE_PROTECT(istring);
01356   FINISH_RESULT(digs, alloc);
01357 
01358   alloc = bigdig_length(digs, alloc);
01359   SCHEME_BIGLEN(o) = alloc;
01360   SCHEME_BIGDIG(o) = digs;
01361 
01362   return scheme_bignum_normalize(o);
01363 }
01364 
01365 Scheme_Object *scheme_read_bignum_bytes(const char *str, int offset, int radix)
01366 {
01367   mzchar *us;
01368 
01369   us = scheme_utf8_decode_to_buffer((unsigned char *)str, 
01370                                 strlen(str XFORM_OK_PLUS offset), 
01371                                 NULL, 0);
01372   return scheme_read_bignum(us, 0, radix);
01373 }
01374 
01375 static void bignum_double_inplace(Scheme_Object **_stk_o)
01376 {
01377   int carry, len;
01378 
01379   len = SCHEME_BIGLEN(*_stk_o);
01380 
01381   if (len == 0)
01382     return;
01383 
01384   /* We assume that *_stk_o is not small */
01385   carry = mpn_lshift(SCHEME_BIGDIG(*_stk_o), SCHEME_BIGDIG(*_stk_o), len, 1);
01386 
01387   if (carry)
01388     *_stk_o = bignum_copy(*_stk_o, carry);
01389 }
01390 
01391 static void bignum_add1_inplace(Scheme_Object **_stk_o)
01392 {
01393   int carry, len;
01394 
01395   len = SCHEME_BIGLEN(*_stk_o);
01396 
01397   if (len == 0) {
01398     *_stk_o = bignum_copy(*_stk_o, 1);
01399     return;
01400   }
01401   /* We assume that *_stk_o is not small */
01402   carry = mpn_add_1(SCHEME_BIGDIG(*_stk_o), SCHEME_BIGDIG(*_stk_o), len, 1);
01403 
01404   if (carry)
01405     *_stk_o = bignum_copy(*_stk_o, carry);
01406 }
01407 
01408 #define USE_FLOAT_BITS 53
01409 #define FP_TYPE double
01410 #define IS_FLOAT_INF scheme__is_double_inf
01411 #define SCHEME_BIGNUM_TO_FLOAT_INFO scheme_bignum_to_double_inf_info
01412 #define SCHEME_BIGNUM_TO_FLOAT scheme_bignum_to_double
01413 #define SCHEME_CHECK_FLOAT scheme_check_double
01414 #define SCHEME_BIGNUM_FROM_FLOAT scheme_bignum_from_double
01415 #include "bgnfloat.inc"
01416 
01417 #ifdef MZ_USE_SINGLE_FLOATS
01418 # undef USE_FLOAT_BITS
01419 # undef FP_TYPE
01420 # undef IS_FLOAT_INF
01421 # undef SCHEME_BIGNUM_TO_FLOAT_INFO
01422 # undef SCHEME_BIGNUM_TO_FLOAT
01423 # undef SCHEME_CHECK_FLOAT
01424 # undef SCHEME_BIGNUM_FROM_FLOAT
01425 
01426 # define USE_FLOAT_BITS 24
01427 # define FP_TYPE float
01428 # define IS_FLOAT_INF scheme__is_float_inf
01429 # define SCHEME_BIGNUM_TO_FLOAT_INFO scheme_bignum_to_float_inf_info
01430 # define SCHEME_BIGNUM_TO_FLOAT scheme_bignum_to_float
01431 # define SCHEME_CHECK_FLOAT scheme_check_float
01432 # define SCHEME_BIGNUM_FROM_FLOAT scheme_bignum_from_float
01433 # include "bgnfloat.inc"
01434 #endif
01435 
01436 
01437 void scheme_bignum_divide(const Scheme_Object *n, const Scheme_Object *d,
01438                        Scheme_Object **_stk_qp, Scheme_Object **_stk_rp, int norm)
01439 {
01440   int cmp;
01441 
01442   cmp = bignum_abs_cmp(n, d);
01443 
01444   if (cmp == -1) {
01445     if (_stk_qp)
01446       *_stk_qp = (norm ? scheme_make_integer(0) : scheme_make_bignum(0));
01447     if (_stk_rp)
01448       *_stk_rp = (norm ? scheme_bignum_normalize(bignum_copy(n, 0)) : bignum_copy(n, 0));
01449     return;
01450   } else if (cmp == 0) {
01451     int n_pos, d_pos, res;
01452 
01453     n_pos = SCHEME_BIGPOS(n);
01454     d_pos = SCHEME_BIGPOS(d);
01455 
01456     res = (xor(n_pos, d_pos) ? -1 : 1);
01457 
01458     if (_stk_qp)
01459       *_stk_qp = (norm ? scheme_make_integer(res) : scheme_make_bignum(res));
01460     if (_stk_rp)
01461       *_stk_rp = (norm ? scheme_make_integer(0) : scheme_make_bignum(0));
01462     return;
01463   } else {
01464     int i;
01465     long n_size, d_size, q_alloc, r_alloc, d_pos;
01466     short n_pos;
01467     bigdig *q_digs, *r_digs, *n_digs, *d_digs;
01468     Scheme_Object *q, *r;
01469     SAFE_SPACE(ns) SAFE_SPACE(ds)
01470 
01471     n_size = SCHEME_BIGLEN(n);
01472     d_size = SCHEME_BIGLEN(d);
01473 
01474     q = (Scheme_Object *)scheme_malloc_tagged(sizeof(Scheme_Bignum));
01475     q->type = scheme_bignum_type;
01476     r = (Scheme_Object *)scheme_malloc_tagged(sizeof(Scheme_Bignum));
01477     r->type = scheme_bignum_type;
01478 
01479     q_alloc = n_size - d_size + 1;
01480     r_alloc = d_size;
01481 
01482     q_digs = PROTECT_RESULT(q_alloc);
01483     r_digs = PROTECT_RESULT(r_alloc);
01484 
01485     n_digs = SCHEME_BIGDIG_SAFE(n, ns);
01486     d_digs = SCHEME_BIGDIG_SAFE(d, ds);
01487     PROTECT(n_digs, n_size);
01488     PROTECT(d_digs, d_size);
01489 
01490     for (i = 0; (i < d_size) && (d_digs[i] == 0); i++) {
01491       r_digs[i] = n_digs[i];
01492     }
01493 
01494     mpn_tdiv_qr(q_digs, r_digs XFORM_OK_PLUS i, 0,
01495               n_digs XFORM_OK_PLUS i, n_size - i,
01496               d_digs XFORM_OK_PLUS i, d_size - i);
01497 
01498     RELEASE(d_digs);
01499     RELEASE(n_digs);
01500     FINISH_RESULT(q_digs, q_alloc);
01501     FINISH_RESULT(r_digs, r_alloc);
01502 
01503     n_pos = SCHEME_BIGPOS(n);
01504     d_pos = SCHEME_BIGPOS(d);
01505 
01506     if (_stk_rp) {
01507       SCHEME_BIGDIG(r) = r_digs;
01508       r_alloc = bigdig_length(r_digs, r_alloc);
01509       SCHEME_BIGLEN(r) = r_alloc;
01510       SCHEME_SET_BIGPOS(r, n_pos);
01511       *_stk_rp = (norm ? scheme_bignum_normalize(r) : r);
01512     }
01513     if (_stk_qp) {
01514       SCHEME_BIGDIG(q) = q_digs;
01515       q_alloc = bigdig_length(q_digs, q_alloc);
01516       SCHEME_BIGLEN(q) = q_alloc;
01517       SCHEME_SET_BIGPOS(q, !xor(n_pos, d_pos));
01518       *_stk_qp = (norm ? scheme_bignum_normalize(q) : q);
01519     }
01520   }
01521 }
01522 
01523 static unsigned long fixnum_sqrt(unsigned long n, unsigned long *rem)
01524 {
01525   unsigned long root = 0;
01526   unsigned long square = 0;
01527   unsigned long try_root, try_square;
01528   int i;
01529 
01530   for (i = SQRT_BIT_MAX; i >= 0; i--)
01531   {
01532     try_root = root | (0x1 << i);
01533     try_square = try_root * try_root;
01534     if (try_square <= n)
01535     {
01536       root = try_root;
01537       square = try_square;
01538     }
01539   }
01540   if (rem)
01541     *rem = n - square;
01542   return root;
01543 }
01544 
01545 Scheme_Object *scheme_integer_sqrt(const Scheme_Object *n)
01546 {
01547   return scheme_integer_sqrt_rem(n, NULL);
01548 }
01549 
01550 Scheme_Object *scheme_integer_sqrt_rem(const Scheme_Object *n, Scheme_Object **remainder)
01551 {
01552   Scheme_Object *o;
01553   int rem_size;
01554 
01555   SAFE_SPACE(qsd)
01556 
01557   if (SCHEME_INTP(n)) {
01558     unsigned long root, rem;
01559     root = fixnum_sqrt(SCHEME_INT_VAL(n), &rem);
01560     if (remainder) {
01561       o = scheme_make_integer_value(rem);
01562       *remainder = o;
01563     }
01564     rem_size = (rem == 0 ? 0 : 1);
01565     o = scheme_make_integer(root);
01566   } else {
01567     long n_size, res_alloc, rem_alloc;
01568     bigdig *res_digs, *rem_digs, *sqr_digs;
01569 
01570     n_size = SCHEME_BIGLEN(n);
01571     if (n_size == 0)
01572       return scheme_make_integer(0);
01573     sqr_digs = SCHEME_BIGDIG_SAFE(n, qsd);
01574 
01575     if (n_size & 0x1)
01576       res_alloc = (n_size + 1) >> 1;
01577     else
01578       res_alloc = n_size >> 1;
01579 
01580     res_digs = PROTECT_RESULT(res_alloc);
01581 
01582     if (remainder)
01583     {
01584       rem_alloc = n_size;
01585       rem_digs = PROTECT_RESULT(rem_alloc);
01586     }
01587     else
01588     {
01589       rem_alloc = 0;
01590       rem_digs = NULL;
01591     }
01592 
01593     PROTECT(sqr_digs, n_size);
01594 
01595     rem_size = mpn_sqrtrem(res_digs, rem_digs, sqr_digs, n_size);
01596 
01597     RELEASE(sqr_digs);
01598 
01599     if (remainder || rem_size == 0) {
01600       /* An integer result */
01601       FINISH_RESULT(res_digs, res_alloc);
01602 
01603       if (remainder && rem_size == 0) {
01604        *remainder = scheme_make_integer(0);
01605        RELEASE(rem_digs);
01606       } else if (remainder) {
01607        Scheme_Object *p;
01608        FINISH_RESULT(rem_digs, rem_alloc);
01609        p = (Scheme_Object *)scheme_malloc_tagged(sizeof(Scheme_Bignum));
01610        p->type = scheme_bignum_type;
01611        rem_alloc = bigdig_length(rem_digs, rem_alloc);
01612        SCHEME_BIGLEN(p) = rem_alloc;
01613        SCHEME_BIGDIG(p) = rem_digs;
01614        SCHEME_SET_BIGPOS(p, 1);
01615        o = scheme_bignum_normalize(p);
01616        *remainder = o;
01617       }
01618 
01619       o = (Scheme_Object *)scheme_malloc_tagged(sizeof(Scheme_Bignum));
01620       o->type = scheme_bignum_type;
01621       res_alloc = bigdig_length(res_digs, res_alloc);
01622       SCHEME_BIGLEN(o) = res_alloc;
01623       SCHEME_BIGDIG(o) = res_digs;
01624       SCHEME_SET_BIGPOS(o, 1);
01625       return scheme_bignum_normalize(o);
01626     } else
01627       o = NULL;
01628     RELEASE(res_digs);
01629   }
01630 
01631   if (remainder || rem_size == 0)
01632     return o;
01633   else {
01634     double v;
01635 
01636     if (SCHEME_INTP(n))
01637       v = (double)SCHEME_INT_VAL(n);
01638     else {
01639       v = scheme_bignum_to_double(n);
01640 
01641       if (MZ_IS_POS_INFINITY(v)) {
01642 #ifdef USE_SINGLE_FLOATS_AS_DEFAULT
01643        return scheme_make_float(v);
01644 #else
01645        return scheme_make_double(v);
01646 #endif
01647       }
01648     }
01649 
01650     v = sqrt(v);
01651 
01652 #ifdef USE_SINGLE_FLOATS_AS_DEFAULT
01653     return scheme_make_float(v);
01654 #else
01655     return scheme_make_double(v);
01656 #endif
01657   }
01658 }
01659 
01660 int gcd_calls = 0;
01661 
01662 Scheme_Object *scheme_bignum_gcd(const Scheme_Object *n, const Scheme_Object *d)
01663 {
01664   bigdig *r_digs, *n_digs, *d_digs;
01665   long n_size, d_size, r_alloc, r_size;
01666   int res_double;
01667   Scheme_Object *r;
01668   SAFE_SPACE(ns) SAFE_SPACE(ds)
01669 
01670   if (scheme_bignum_lt(d, n)) {
01671     const Scheme_Object *tmp;
01672     tmp = n;
01673     n = d;
01674     d = tmp;
01675   }
01676 
01677   n_size = SCHEME_BIGLEN(n);
01678   d_size = SCHEME_BIGLEN(d);
01679 
01680   if (!n_size)
01681     return (Scheme_Object *)d;
01682 
01683   r = (Scheme_Object *)scheme_malloc_tagged(sizeof(Scheme_Bignum));
01684   r->type = scheme_bignum_type;
01685 
01686 #ifdef MZ_PRECISE_GC
01687   n_digs = SCHEME_BIGDIG_SAFE(n, ns);
01688   d_digs = SCHEME_BIGDIG_SAFE(d, ds);
01689   PROTECT(n_digs, n_size);
01690   PROTECT(d_digs, d_size);
01691 #else
01692   n_digs = (bigdig *)scheme_malloc_atomic(sizeof(bigdig) * n_size);
01693   d_digs = (bigdig *)scheme_malloc_atomic(sizeof(bigdig) * d_size);
01694   memcpy(n_digs, SCHEME_BIGDIG(n), sizeof(bigdig) * n_size);
01695   memcpy(d_digs, SCHEME_BIGDIG(d), sizeof(bigdig) * d_size);
01696 #endif
01697 
01698   /* GMP wants the first argument to be odd. Compute a shift. */
01699   {
01700     bigdig mask;
01701     int b, w, nz = 0, dz = 0;
01702 
01703     b = 1; w = 0; mask = 0x1;
01704     while (!(n_digs[w] & mask)) {
01705       nz++;
01706       if (b == WORD_SIZE) {
01707        b = 1;
01708        mask = 0x1;
01709        w++;
01710       } else {
01711        b++;
01712        mask = mask << 1;
01713       }
01714     }
01715 
01716     b = 1; w = 0; mask = 0x1;
01717     while ((dz < nz) && !(d_digs[w] & mask)) {
01718       dz++;
01719       if (b == WORD_SIZE) {
01720        b = 1;
01721        mask = 0x1;
01722        w++;
01723       } else {
01724        b++;
01725        mask = mask << 1;
01726       }
01727     }
01728 
01729     if (nz) {
01730       w = nz / WORD_SIZE;
01731       memmove(n_digs, n_digs + w, sizeof(bigdig) * (n_size - w));
01732       n_size -= w;
01733       w = nz & (WORD_SIZE - 1);
01734       if (w)
01735        mpn_rshift(n_digs, n_digs, n_size, w);
01736     }
01737     if (dz) {
01738       w = dz / WORD_SIZE;
01739       memmove(d_digs, d_digs + w, sizeof(bigdig) * (d_size - w));
01740       d_size -= w;
01741       w = dz & (WORD_SIZE - 1);
01742       if (w)
01743        mpn_rshift(d_digs, d_digs, d_size, w);
01744     }
01745 
01746     if (nz < dz)
01747       res_double = nz;
01748     else
01749       res_double = dz;
01750 
01751     /* Most-significant word must be non-zero: */
01752     if (!(n_digs[n_size - 1]))
01753       --n_size;
01754     if (!(d_digs[d_size - 1]))
01755       --d_size;
01756   }
01757 
01758   r_alloc = n_size;
01759 
01760   r_digs = PROTECT_RESULT(r_alloc);
01761 
01762   r_size = mpn_gcd(r_digs, d_digs, d_size, n_digs, n_size);
01763 
01764   RELEASE(d_digs);
01765   RELEASE(n_digs);
01766   FINISH_RESULT(r_digs, r_size);
01767 
01768   SCHEME_BIGDIG(r) = r_digs;
01769   r_alloc = bigdig_length(r_digs, r_size);
01770   SCHEME_BIGLEN(r) = r_alloc;
01771   SCHEME_SET_BIGPOS(r, 1);
01772 
01773   if (res_double)
01774     return scheme_bignum_shift(r, res_double);
01775   else
01776     return scheme_bignum_normalize(r);
01777 }
01778 
01779 /* Used by GMP library (which is not xformed for precise GC): */
01780 void scheme_bignum_use_fuel(long n)
01781 {
01782 #ifdef MZ_PRECISE_GC
01783 # ifndef GC_STACK_CALLEE_RESTORE
01784   char *stupid; /* forces __gc_var_stack__ */
01785 # endif
01786 #endif
01787 
01788   SCHEME_USE_FUEL(n);
01789 
01790 #ifdef MZ_PRECISE_GC
01791 # ifndef GC_STACK_CALLEE_RESTORE
01792   /* Restore variable stack. */
01793   if (!stupid)
01794     GC_variable_stack = (void **)__gc_var_stack__[0];
01795 # endif
01796 #endif
01797 }