1 // This file is part of Golly.
2 // See docs/License.html for the copyright notice.
3 
4 #include "bigint.h"
5 #include <cstdlib>
6 #include <string.h>
7 #include <iostream>
8 #include <cmath>
9 #include <limits.h>
10 #include "util.h"
11 #undef SLOWCHECK
12 using namespace std ;
13 /**
14  *   Static data.
15  */
16 static const int MAX_SIMPLE = 0x3fffffff ;
17 static const int MIN_SIMPLE = -0x40000000 ;
18 char *bigint::printbuf ;
19 int *bigint::work ;
20 int bigint::printbuflen ;
21 int bigint::workarrlen ;
22 char bigint::sepchar = ',' ;
23 int bigint::sepcount = 3 ;
24 /**
25  *   Routines.
26  */
bigint(G_INT64 i)27 bigint::bigint(G_INT64 i) {
28    if (i <= INT_MAX && i >= INT_MIN)
29       fromint((int)i) ;
30    else {
31       v.i = 0 ;
32       vectorize(0) ;
33       v.p[0] = 3 ;
34       v.p[1] = (int)(i & 0x7fffffff) ;
35       v.p[2] = (int)((i >> 31) & 0x7fffffff) ;
36       v.p[3] = 0 ;
37       ripple((int)(i >> 62), 3) ;
38    }
39 }
40 // we can parse ####, 2^###, -#####
41 // AKT: we ignore all non-digits (except for leading '-')
42 // so we can parse strings like "1,234" or "+1.234";
43 // it is up to caller to impose smarter restrictions
bigint(const char * s)44 bigint::bigint(const char *s) {
45    if (*s == '2' && s[1] == '^') {
46       long x = atol(s+2) ;
47       if (x < 31)
48          fromint(1 << x) ;
49       else {
50          int sz = 2 + int((x + 1) / 31) ;
51          int asz = sz ;
52          while (asz & (asz - 1))
53             asz &= asz - 1 ;
54          asz *= 2 ;
55          v.p = new int[asz] ;
56          v.p[0] = sz ;
57          for (int i=1; i<=sz; i++)
58             v.p[i] = 0 ;
59          v.p[sz-1] = 1 << (x % 31) ;
60       }
61    } else {
62       int neg = 0 ;
63       if (*s == '-') {
64          neg = 1 ;
65          s++ ;
66       }
67       fromint(0) ;
68       while (*s) {
69          // AKT: was *s != sepchar
70          if (*s >= '0' && *s <= '9') {
71             mul_smallint(10) ;
72             if (neg)
73                add_smallint('0'-*s) ;
74             else
75                add_smallint(*s-'0') ;
76          }
77          s++ ;
78       }
79    }
80 }
copyarr(int * p)81 static int *copyarr(int *p) {
82    int sz = *p ;
83    while (sz & (sz - 1))
84       sz &= sz - 1 ;
85    sz *= 2 ;
86    int *r = new int[sz] ;
87    memcpy(r, p, sizeof(int) * (p[0] + 1)) ;
88 #ifdef SLOWCHECK
89    for (int i=p[0]+1; i<sz; i++)
90       r[i] = 0xdeadbeef ;
91 #endif
92    return r ;
93 }
bigint(const bigint & a)94 bigint::bigint(const bigint &a) {
95    if (a.v.i & 1)
96       v.i = a.v.i ;
97    else
98       v.p = copyarr(a.v.p) ;
99 }
100 /**
101  *   We special-case the p=0 case so we can "initialize" bigint memory
102  *   with zero.
103  */
operator =(const bigint & b)104 bigint &bigint::operator=(const bigint &b) {
105    if (&b != this) {
106       if (0 == (v.i & 1))
107          if (v.p)
108             delete [] v.p ;
109       if (b.v.i & 1)
110          v.i = b.v.i ;
111       else
112          v.p = copyarr(b.v.p) ;
113    }
114    return *this ;
115 }
~bigint()116 bigint::~bigint() {
117    if (0 == (v.i & 1))
118       delete [] v.p ;
119 }
bigint(const bigint & a,const bigint & b,const bigint & c,const bigint & d)120 bigint::bigint(const bigint &a, const bigint &b, const bigint &c, const bigint &d) {
121    const int checkmask = 0xf0000001 ;
122    if ((a.v.i & checkmask) == 1 && (b.v.i & checkmask) == 1 &&
123        (c.v.i & checkmask) == 1 && (d.v.i & checkmask) == 1) {
124       // hot path
125       v.i = a.v.i + b.v.i + c.v.i + d.v.i - 3 ;
126       return ;
127    }
128    v.i = 1 ;
129    *this = a ;
130    *this += b ;
131    *this += c ;
132    *this += d ;
133 }
tostring(char sep) const134 const char *bigint::tostring(char sep) const {
135    int lenreq = 32 ;
136    if ((v.i & 1) == 0)
137       lenreq = size() * 32 ;
138    if (lenreq > printbuflen) {
139      if (printbuf)
140          delete [] printbuf ;
141       printbuf = new char[2 * lenreq] ;
142       printbuflen = 2 * lenreq ;
143    }
144    int sz = 1 ;
145    if (0 == (v.i & 1))
146       sz = size() ;
147    ensurework(sz) ;
148    int neg = sign() < 0 ;
149    if (v.i & 1) {
150       if (neg)
151          work[0] = -(v.i >> 1) ;
152       else
153          work[0] = v.i >> 1 ;
154    } else {
155       if (neg) {
156          int carry = 1 ;
157          for (int i=0; i+1<sz; i++) {
158             int c = carry + (v.p[i+1] ^ 0x7fffffff) ;
159             work[i] = c & 0x7fffffff ;
160             carry = (c >> 31) & 1 ;
161          }
162          work[sz-1] = carry + ~v.p[sz] ;
163       } else {
164          for (int i=0; i<sz; i++)
165             work[i] = v.p[i+1] ;
166       }
167    }
168    char *p = printbuf ;
169    const int bigradix = 1000000000 ; // 9 digits at a time
170    for (;;) {
171       int allbits = 0 ;
172       int carry = 0 ;
173       int i;
174       for (i=sz-1; i>=0; i--) {
175          G_INT64 c = carry * G_MAKEINT64(0x80000000) + work[i] ;
176          carry = (int)(c % bigradix) ;
177          work[i] = (int)(c / bigradix) ;
178          allbits |= work[i] ;
179       }
180       for (i=0; i<9; i++) { // put the nine digits in
181          *p++ = (char)(carry % 10 + '0') ;
182          carry /= 10 ;
183       }
184       if (allbits == 0)
185          break ;
186    }
187    while (p > printbuf + 1 && *(p-1) == '0')
188       p-- ;
189    char *r = p ;
190    if (neg)
191       *r++ = '-' ;
192    for (int i=(int)(p-printbuf-1); i>=0; i--) {
193       *r++ = printbuf[i] ;
194       if (i && sep && (i % sepcount == 0))
195          *r++ = sep ;
196    }
197    *r++ = 0 ;
198    return p ;
199 }
grow(int osz,int nsz)200 void bigint::grow(int osz, int nsz) {
201    int bdiffs = osz ^ nsz ;
202    if (bdiffs > osz) {
203       while (bdiffs & (bdiffs - 1))
204          bdiffs &= bdiffs - 1 ;
205       int *nv = new int[2*bdiffs] ;
206       for (int i=0; i<=osz; i++)
207          nv[i] = v.p[i] ;
208 #ifdef SLOWCHECK
209       for (int i=osz+1; i<2*bdiffs; i++)
210          nv[i] = 0xdeadbeef ;
211 #endif
212       delete [] v.p ;
213       v.p = nv ;
214    }
215    int av = v.p[osz] ;
216    while (osz < nsz) {
217       v.p[osz] = av & 0x7fffffff ;
218       osz++ ;
219    }
220    v.p[nsz] = av ;
221    v.p[0] = nsz ;
222 }
operator +=(const bigint & a)223 bigint& bigint::operator+=(const bigint &a) {
224    if (a.v.i & 1)
225       add_smallint(a.v.i >> 1) ;
226    else {
227       if (v.i & 1)
228          vectorize(v.i >> 1) ;
229       ripple(a, 0) ;
230    }
231    return *this ;
232 }
operator -=(const bigint & a)233 bigint& bigint::operator-=(const bigint &a) {
234    if (a.v.i & 1)
235       add_smallint(-(a.v.i >> 1)) ;
236    else {
237       if (v.i & 1)
238          vectorize(v.i >> 1) ;
239       ripplesub(a, 1) ;
240    }
241    return *this ;
242 }
sign() const243 int bigint::sign() const {
244    int si = v.i ;
245    if (0 == (si & 1))
246       si = v.p[size()] ;
247    if (si > 0)
248       return 1 ;
249    if (si < 0)
250       return -1 ;
251    return 0 ;
252 }
size() const253 int bigint::size() const {
254    return v.p[0] ;
255 }
add_smallint(int a)256 void bigint::add_smallint(int a) {
257    if (v.i & 1)
258       fromint((v.i >> 1) + a) ;
259    else
260       ripple(a, 1) ;
261 }
262 // do we need to shrink it to keep it canonical?
shrink(int pos)263 void bigint::shrink(int pos) {
264    while (pos > 1 && ((v.p[pos] - v.p[pos-1]) & 0x7fffffff) == 0) {
265       pos-- ;
266       v.p[pos] = v.p[pos+1] ;
267       v.p[0] = pos ;
268    }
269    if (pos == 1) {
270       int c = v.p[1] ;
271       delete [] v.p ;
272       v.i = 1 | (c << 1) ;
273    } else if (pos == 2 && ((v.p[2] ^ v.p[1]) & 0x40000000) == 0) {
274       int c = v.p[1] + (v.p[2] << 31) ;
275       delete [] v.p ;
276       v.i = 1 | (c << 1) ;
277    }
278 }
279 void grow(int osz, int nsz) ;
280 // note:  carry may be any legal 31-bit int, *or* positive 2^30
281 // note:  may only be called on arrayed bigints
ripple(int carry,int pos)282 void bigint::ripple(int carry, int pos) {
283    int c ;
284    int sz = size() ;
285    while (pos < sz) {
286       c = v.p[pos] + (carry & 0x7fffffff) ;
287       carry = ((c >> 31) & 1) + (carry >> 31) ;
288       v.p[pos++] = c & 0x7fffffff;
289    }
290    c = v.p[pos] + carry ;
291    if (c == 0 || c == -1) { // see if we can make it smaller
292       v.p[pos] = c ;
293       shrink(pos) ;
294    } else { // need to extend
295       if (0 == (sz & (sz + 1))) { // grow array (oops)
296          grow(sz, sz+1) ;
297       } else
298          v.p[0] = sz + 1 ;
299       v.p[pos] = c & 0x7fffffff ;
300       v.p[pos+1] = -((c >> 31) & 1) ;
301    }
302 }
ripple(const bigint & a,int carry)303 void bigint::ripple(const bigint &a, int carry) {
304    int asz = a.size() ;
305    int tsz = size() ;
306    int pos = 1 ;
307    if (tsz < asz) { // gotta resize
308       grow(tsz, asz) ;
309       tsz = asz ;
310    }
311    while (pos < asz) {
312       int c = v.p[pos] + a.v.p[pos] + carry ;
313       carry = (c >> 31) & 1 ;
314       v.p[pos++] = c & 0x7fffffff;
315    }
316    ripple(carry + a.v.p[pos], pos) ;
317 }
ripplesub(const bigint & a,int carry)318 void bigint::ripplesub(const bigint &a, int carry) {
319    int asz = a.size() ;
320    int tsz = size() ;
321    int pos = 1 ;
322    if (tsz < asz) { // gotta resize
323       grow(tsz, asz) ;
324       tsz = asz ;
325    }
326    while (pos < asz) {
327       int c = v.p[pos] + (0x7fffffff ^ a.v.p[pos]) + carry ;
328       carry = (c >> 31) & 1 ;
329       v.p[pos++] = c & 0x7fffffff;
330    }
331    ripple(carry + ~a.v.p[pos], pos) ;
332 }
333 // make sure it's in vector form; may leave it not canonical!
vectorize(int i)334 void bigint::vectorize(int i) {
335    v.p = new int[4] ;
336    v.p[0] = 2 ;
337    v.p[1] = i & 0x7fffffff ;
338    if (i < 0)
339       v.p[2] = -1 ;
340    else
341       v.p[2] = 0 ;
342 #ifdef SLOWCHECK
343    v.p[3] = 0xdeadbeef ;
344 #endif
345 }
fromint(int i)346 void bigint::fromint(int i) {
347    if (i <= MAX_SIMPLE && i >= MIN_SIMPLE)
348       v.i = (i << 1) | 1 ;
349    else
350       vectorize(i) ;
351 }
ensurework(int sz) const352 void bigint::ensurework(int sz) const {
353    sz += 3 ;
354    if (sz > workarrlen) {
355       if (work)
356          delete [] work ;
357       work = new int[sz * 2] ;
358       workarrlen = sz * 2 ;
359    }
360 }
mul_smallint(int a)361 void bigint::mul_smallint(int a) {
362    int c ;
363    if (a == 0) {
364       *this = 0 ;
365       return ;
366    }
367    if (v.i & 1) {
368       if ((v.i >> 1) <= MAX_SIMPLE / a) {
369          c = (v.i >> 1) * a ;
370          fromint((v.i >> 1) * a) ;
371          return ;
372       }
373       vectorize(v.i >> 1) ;
374    }
375    int sz = size() ;
376    int carry = 0 ;
377    int pos = 1 ;
378    while (pos < sz) {
379       int clo = (v.p[pos] & 0xffff) * a + carry ;
380       int chi = (v.p[pos] >> 16) * a + (clo >> 16) ;
381       carry = chi >> 15 ;
382       v.p[pos++] = (clo & 0xffff) + ((chi & 0x7fff) << 16) ;
383    }
384    c = v.p[pos] * a + carry ;
385    if (c == 0 || c == -1) { // see if we can make it smaller
386       v.p[pos] = c ;
387       shrink(pos) ;
388    } else { // need to extend
389       if (0 == (sz & (sz + 1))) { // grow array (oops)
390          grow(sz, sz+1) ;
391       } else
392          v.p[0] = sz + 1 ;
393       v.p[pos] = c & 0x7fffffff ;
394       v.p[pos+1] = -((c >> 31) & 1) ;
395    }
396 }
div_smallint(int a)397 void bigint::div_smallint(int a) {
398    if (v.i & 1) {
399       int r = (v.i >> 1) / a ;
400       fromint(r) ;
401       return ;
402    }
403    if (v.p[v.p[0]] < 0)
404       lifefatal("we don't support divsmallint when negative yet") ;
405    int carry = 0 ;
406    int pos = v.p[0] ;
407    while (pos > 0) {
408       G_INT64 t = ((G_INT64)carry << G_MAKEINT64(31)) + v.p[pos] ;
409       carry = (int)(t % a) ;
410       v.p[pos] = (int)(t / a) ;
411       pos-- ;
412    }
413    shrink(v.p[0]) ;
414 }
mod_smallint(int a)415 int bigint::mod_smallint(int a) {
416    if (v.i & 1)
417       return (((v.i >> 1) % a) + a) % a ;
418    int pos = v.p[0] ;
419    int mm = (2 * ((1 << 30) % a) % a) ;
420    int r = 0 ;
421    while (pos > 0) {
422       r = (mm * r + v.p[pos]) % a ;
423       pos-- ;
424    }
425    return (r + a) % a ;
426 }
div2()427 void bigint::div2() {
428    if (v.i & 1) {
429       v.i = ((v.i >> 1) | 1) ;
430       return ;
431    }
432    int sz = v.p[0] ;
433    int carry = -v.p[sz] ;
434    for (int i=sz-1; i>0; i--) {
435       int c = (v.p[i] >> 1) + (carry << 30) ;
436       carry = v.p[i] & 1 ;
437       v.p[i] = c ;
438    }
439    shrink(v.p[0]) ;
440 }
operator >>=(int i)441 bigint& bigint::operator>>=(int i) {
442    if (v.i & 1) {
443       if (i > 31)
444          v.i = ((v.i >> 31) | 1) ;
445       else
446          v.i = ((v.i >> i) | 1) ;
447       return *this ;
448    }
449    int bigsh = i / 31 ;
450    if (bigsh) {
451       if (bigsh >= v.p[0]) {
452          v.p[1] = v.p[v.p[0]] ;
453          v.p[0] = 1 ;
454       } else {
455          for (int j=1; j+bigsh<=v.p[0]; j++)
456             v.p[j] = v.p[j+bigsh] ;
457          v.p[0] -= bigsh ;
458       }
459       i -= bigsh * 31 ;
460    }
461    int carry = v.p[v.p[0]] ;
462    if (i) {
463       for (int j=v.p[0]-1; j>0; j--) {
464          int c = ((v.p[j] >> i) | (carry << (31-i))) & 0x7fffffff ;
465          carry = v.p[j] ;
466          v.p[j] = c ;
467       }
468    }
469    shrink(v.p[0]) ;
470    return *this ;
471 }
operator <<=(int i)472 bigint& bigint::operator<<=(int i) {
473    if (v.i & 1) {
474       if (v.i == 1)
475          return *this ;
476       if (i < 30 && (v.i >> 31) == (v.i >> (31 - i))) {
477          v.i = ((v.i & ~1) << i) | 1 ;
478          return *this ;
479       }
480       vectorize(v.i >> 1) ;
481    }
482    int bigsh = i / 31 ;
483    int nsize = v.p[0] + bigsh + 1 ; // how big we need it to be, worst case
484    grow(v.p[0], nsize) ;
485    if (bigsh) {
486       int j ;
487       for (j=v.p[0]-1; j>bigsh; j--)
488          v.p[j] = v.p[j-bigsh] ;
489       for (j=bigsh; j>0; j--)
490          v.p[j] = 0 ;
491       i -= bigsh * 31 ;
492    }
493    int carry = 0 ;
494    if (i) {
495       for (int j=1; j<v.p[0]; j++) {
496          int c = ((v.p[j] << i) | (carry >> (31-i))) & 0x7fffffff ;
497          carry = v.p[j] ;
498          v.p[j] = c ;
499       }
500    }
501    shrink(v.p[0]) ;
502    return *this ;
503 }
mulpow2(int p)504 void bigint::mulpow2(int p) {
505    if (p > 0)
506       *this <<= p ;
507    else if (p < 0)
508       *this >>= -p ;
509 }
even() const510 int bigint::even() const {
511    if (v.i & 1)
512       return 1-((v.i >> 1) & 1) ;
513    else
514       return 1-(v.p[1] & 1) ;
515 }
odd() const516 int bigint::odd() const {
517    if (v.i & 1)
518       return ((v.i >> 1) & 1) ;
519    else
520       return (v.p[1] & 1) ;
521 }
low31() const522 int bigint::low31() const {
523    if (v.i & 1)
524       return ((v.i >> 1) & 0x7fffffff) ;
525    else
526       return v.p[1] ;
527 }
operator ==(const bigint & b) const528 int bigint::operator==(const bigint &b) const {
529    if ((b.v.i - v.i) & 1)
530       return 0 ;
531    if (v.i & 1)
532       return v.i == b.v.i ;
533    if (b.v.p[0] != v.p[0])
534       return 0 ;
535    return memcmp(v.p, b.v.p, sizeof(int) * (v.p[0] + 1)) == 0 ;
536 }
operator !=(const bigint & b) const537 int bigint::operator!=(const bigint &b) const {
538    return !(*this == b) ;
539 }
operator <(const bigint & b) const540 int bigint::operator<(const bigint &b) const {
541    if (b.v.i & 1)
542       if (v.i & 1)
543          return v.i < b.v.i ;
544       else
545          return v.p[v.p[0]] < 0 ;
546    else
547       if (v.i & 1)
548          return b.v.p[b.v.p[0]] >= 0 ;
549    int d = v.p[v.p[0]] - b.v.p[b.v.p[0]] ;
550    if (d < 0)
551       return 1 ;
552    if (d > 0)
553       return 0 ;
554    if (v.p[0] > b.v.p[0])
555       return v.p[v.p[0]] < 0 ;
556    else if (v.p[0] < b.v.p[0])
557       return v.p[v.p[0]] >= 0 ;
558    for (int i=v.p[0]; i>0; i--)
559       if (v.p[i] < b.v.p[i])
560          return 1 ;
561       else if (v.p[i] > b.v.p[i])
562          return 0 ;
563    return 0 ;
564 }
operator <=(const bigint & b) const565 int bigint::operator<=(const bigint &b) const {
566    if (b.v.i & 1)
567       if (v.i & 1)
568          return v.i <= b.v.i ;
569       else
570          return v.p[v.p[0]] < 0 ;
571    else
572       if (v.i & 1)
573          return b.v.p[b.v.p[0]] >= 0 ;
574    int d = v.p[v.p[0]] - b.v.p[b.v.p[0]] ;
575    if (d < 0)
576       return 1 ;
577    if (d > 0)
578       return 0 ;
579    if (v.p[0] > b.v.p[0])
580       return v.p[v.p[0]] < 0 ;
581    else if (v.p[0] < b.v.p[0])
582       return v.p[v.p[0]] >= 0 ;
583    for (int i=v.p[0]; i>0; i--)
584       if (v.p[i] < b.v.p[i])
585          return 1 ;
586       else if (v.p[i] > b.v.p[i])
587          return 0 ;
588    return 1 ;
589 }
operator >(const bigint & b) const590 int bigint::operator>(const bigint &b) const {
591    return !(*this <= b) ;
592 }
operator >=(const bigint & b) const593 int bigint::operator>=(const bigint &b) const {
594    return !(*this < b) ;
595 }
mybpow(int n)596 static double mybpow(int n) {
597    double r = 1 ;
598    double s = 65536.0 * 32768.0 ;
599    while (n) {
600       if (n & 1)
601          r *= s ;
602       s *= s ;
603       n >>= 1 ;
604    }
605    return r ;
606 }
607 /**
608  *   Turn this bigint into a double.
609  */
todouble() const610 double bigint::todouble() const {
611    if (v.i & 1)
612       return (double)(v.i >> 1) ;
613    double r = 0 ;
614    double m = 1 ;
615    int lim = 1 ;
616    if (v.p[0] > 4) {
617       lim = v.p[0] - 3 ;
618       m = mybpow(lim-1) ;
619    }
620    for (int i=lim; i<=v.p[0]; i++) {
621       r = r + m * v.p[i] ;
622       m *= 65536 * 32768.0 ;
623    }
624    return r ;
625 }
626 /**
627  * Turn this bigint into a double in a way that preserves huge exponents.
628  * Here are some examples:
629  *
630  * To represent
631  * the quantity     We return the value
632  *  27               1.27
633  * -6.02e23         -23.602
634  *  6.02e23          23.602
635  *  9.99e299         299.999
636  *  1.0e300          300.1
637  *  1.0e1000         1000.1       This would normally overflow
638  *  6.7e12345        12345.67
639  */
toscinot() const640 double bigint::toscinot() const
641 {
642   double mant, exponent;
643   double k_1_10 = 0.1;
644   double k_1_10000 = 0.0001;
645   double k_base = 65536.0 * 32768.0;
646 
647   exponent = 0;
648   if (v.i & 1) {
649     /* small integer */
650     mant = (double)(v.i >> 1) ;
651   } else {
652     /* big integer: a string of 31-bit chunks */
653     mant = 0 ;
654     double m = 1 ;
655     for (int i=1; i<=v.p[0]; i++) {
656       mant = mant + m * v.p[i] ;
657       m *= k_base ;
658       while (m >= 100000.0) {
659         m *= k_1_10000;
660         mant *= k_1_10000;
661         exponent += 4.0;
662       }
663     }
664   }
665 
666   /* Add the last few powers of 10 back into the mantissa */
667   while(((mant < 0.5) && (mant > -0.5)) && (exponent > 0.0)) {
668     mant *= 10.0;
669     exponent -= 1.0;
670   }
671 
672   /* Mantissa might be greater than 1 at this point */
673   while((mant >= 1.0) || (mant <= -1.0)) {
674     mant *= k_1_10;
675     exponent += 1.0;
676   }
677 
678   if (mant >= 0) {
679     // Normal case: 6.02e23 -> 23.602
680     return(exponent + mant);
681   }
682   // Negative case: -6.02e23 -> -23.602
683   // In this example mant will be -0.602 and exponent will be 23.0
684   return(mant - exponent);
685 }
686 /**
687  *   Return an int.
688  */
toint() const689 int bigint::toint() const {
690    if (v.i & 1)
691       return (v.i >> 1) ;
692    return (v.p[v.p[0]] << 31) | v.p[1] ;
693 }
694 /**
695  *   How many bits required to represent this, approximately?
696  *   Should overestimate but not by too much.
697  */
bitsreq() const698 int bigint::bitsreq() const {
699    if (v.i & 1)
700       return 31 ;
701    return v.p[0] * 31 ;
702 }
703 /**
704  *   Find the lowest bit set.
705  */
lowbitset() const706 int bigint::lowbitset() const {
707    int r = 0 ;
708    if (v.i & 1) {
709       if (v.i == 1)
710          return -1 ;
711       for (int i=1; i<32; i++)
712          if ((v.i >> i) & 1)
713             return i-1 ;
714    }
715    int o = 1 ;
716    while (v.p[o] == 0) {
717       o++ ;
718       r += 31 ;
719    }
720    for (int i=0; i<32; i++)
721       if ((v.p[o] >> i) & 1)
722          return i + r ;
723    return -1 ;
724 }
725 /**
726  *   Fill a character array with the bits.  Top one is always
727  *   the sign bit.
728  */
tochararr(char * fillme,int n) const729 void bigint::tochararr(char *fillme, int n) const {
730    int at = 0 ;
731    while (n > 0) {
732       int w = 0 ;
733       if (v.i & 1) {
734          if (at == 0)
735             w = v.i >> 1 ;
736          else
737             w = v.i >> 31 ;
738       } else {
739          if (at < v.p[0])
740             w = v.p[at+1] ;
741          else
742             w = v.p[v.p[0]] ;
743       }
744       int lim = 31 ;
745       if (n < 31)
746          lim = n ;
747       for (int i=0; i<lim; i++) {
748          *fillme++ = (char)(w & 1) ;
749          w >>= 1 ;
750       }
751       n -= 31 ;
752       at++ ;
753    }
754 }
755 /**
756  *   Manifests.
757  */
758 const bigint bigint::zero(0) ;
759 const bigint bigint::one(1) ;
760 const bigint bigint::two(2) ;
761 const bigint bigint::three(3) ;
762 const bigint bigint::maxint(INT_MAX) ;
763 const bigint bigint::minint(INT_MIN) ;
764 
765 // most editing operations are limited to absolute coordinates <= 10^9,
766 // partly because getcell and setcell only take int parameters, but mostly
767 // to avoid ridiculously long cut/copy/paste/rotate/etc operations
768 const bigint bigint::min_coord(-1000000000) ;
769 const bigint bigint::max_coord(+1000000000) ;
770