1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998, 1999 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and
9its documentation for any purpose and without fee is hereby
10granted, provided that the above copyright notice appear in all
11copies and that both that the copyright notice and this
12permission notice and warranty disclaimer appear in supporting
13documentation, and that the name of Lucent or any of its entities
14not be used in advertising or publicity pertaining to
15distribution of the software without specific, written prior
16permission.
17
18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25THIS SOFTWARE.
26
27****************************************************************/
28
29/* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
31
32#include "gdtoaimp.h"
33
34 static Bigint *freelist[Kmax+1];
35#ifndef Omit_Private_Memory
36#ifndef PRIVATE_MEM
37#define PRIVATE_MEM 2304
38#endif
39#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
40static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
41#endif
42
43 Bigint *
44Balloc
45#ifdef KR_headers
46 (k) int k;
47#else
48 (int k)
49#endif
50{
51 int x;
52 Bigint *rv;
53#ifndef Omit_Private_Memory
54 unsigned int len;
55#endif
56
57 ACQUIRE_DTOA_LOCK(0);
58 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
59 /* but this case seems very unlikely. */
60 if (k <= Kmax && (rv = freelist[k]) !=0) {
61 freelist[k] = rv->next;
62 }
63 else {
64 x = 1 << k;
65#ifdef Omit_Private_Memory
66 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
67#else
68 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
69 /sizeof(double);
70 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
71 rv = (Bigint*)pmem_next;
72 pmem_next += len;
73 }
74 else
75 rv = (Bigint*)MALLOC(len*sizeof(double));
76#endif
77 rv->k = k;
78 rv->maxwds = x;
79 }
80 FREE_DTOA_LOCK(0);
81 rv->sign = rv->wds = 0;
82 return rv;
83 }
84
85 void
86Bfree
87#ifdef KR_headers
88 (v) Bigint *v;
89#else
90 (Bigint *v)
91#endif
92{
93 if (v) {
94 if (v->k > Kmax)
| 1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998, 1999 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and
9its documentation for any purpose and without fee is hereby
10granted, provided that the above copyright notice appear in all
11copies and that both that the copyright notice and this
12permission notice and warranty disclaimer appear in supporting
13documentation, and that the name of Lucent or any of its entities
14not be used in advertising or publicity pertaining to
15distribution of the software without specific, written prior
16permission.
17
18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25THIS SOFTWARE.
26
27****************************************************************/
28
29/* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
31
32#include "gdtoaimp.h"
33
34 static Bigint *freelist[Kmax+1];
35#ifndef Omit_Private_Memory
36#ifndef PRIVATE_MEM
37#define PRIVATE_MEM 2304
38#endif
39#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
40static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
41#endif
42
43 Bigint *
44Balloc
45#ifdef KR_headers
46 (k) int k;
47#else
48 (int k)
49#endif
50{
51 int x;
52 Bigint *rv;
53#ifndef Omit_Private_Memory
54 unsigned int len;
55#endif
56
57 ACQUIRE_DTOA_LOCK(0);
58 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
59 /* but this case seems very unlikely. */
60 if (k <= Kmax && (rv = freelist[k]) !=0) {
61 freelist[k] = rv->next;
62 }
63 else {
64 x = 1 << k;
65#ifdef Omit_Private_Memory
66 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
67#else
68 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
69 /sizeof(double);
70 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
71 rv = (Bigint*)pmem_next;
72 pmem_next += len;
73 }
74 else
75 rv = (Bigint*)MALLOC(len*sizeof(double));
76#endif
77 rv->k = k;
78 rv->maxwds = x;
79 }
80 FREE_DTOA_LOCK(0);
81 rv->sign = rv->wds = 0;
82 return rv;
83 }
84
85 void
86Bfree
87#ifdef KR_headers
88 (v) Bigint *v;
89#else
90 (Bigint *v)
91#endif
92{
93 if (v) {
94 if (v->k > Kmax)
|
| 95#ifdef FREE
96 FREE((void*)v);
97#else
|
95 free((void*)v);
| 98 free((void*)v);
|
| 99#endif
|
96 else {
97 ACQUIRE_DTOA_LOCK(0);
98 v->next = freelist[v->k];
99 freelist[v->k] = v;
100 FREE_DTOA_LOCK(0);
101 }
102 }
103 }
104
105 int
106lo0bits
107#ifdef KR_headers
108 (y) ULong *y;
109#else
110 (ULong *y)
111#endif
112{
| 100 else {
101 ACQUIRE_DTOA_LOCK(0);
102 v->next = freelist[v->k];
103 freelist[v->k] = v;
104 FREE_DTOA_LOCK(0);
105 }
106 }
107 }
108
109 int
110lo0bits
111#ifdef KR_headers
112 (y) ULong *y;
113#else
114 (ULong *y)
115#endif
116{
|
113 register int k;
114 register ULong x = *y;
| 117 int k;
118 ULong x = *y;
|
115
116 if (x & 7) {
117 if (x & 1)
118 return 0;
119 if (x & 2) {
120 *y = x >> 1;
121 return 1;
122 }
123 *y = x >> 2;
124 return 2;
125 }
126 k = 0;
127 if (!(x & 0xffff)) {
128 k = 16;
129 x >>= 16;
130 }
131 if (!(x & 0xff)) {
132 k += 8;
133 x >>= 8;
134 }
135 if (!(x & 0xf)) {
136 k += 4;
137 x >>= 4;
138 }
139 if (!(x & 0x3)) {
140 k += 2;
141 x >>= 2;
142 }
143 if (!(x & 1)) {
144 k++;
145 x >>= 1;
146 if (!x)
147 return 32;
148 }
149 *y = x;
150 return k;
151 }
152
153 Bigint *
154multadd
155#ifdef KR_headers
156 (b, m, a) Bigint *b; int m, a;
157#else
158 (Bigint *b, int m, int a) /* multiply by m and add a */
159#endif
160{
161 int i, wds;
162#ifdef ULLong
163 ULong *x;
164 ULLong carry, y;
165#else
166 ULong carry, *x, y;
167#ifdef Pack_32
168 ULong xi, z;
169#endif
170#endif
171 Bigint *b1;
172
173 wds = b->wds;
174 x = b->x;
175 i = 0;
176 carry = a;
177 do {
178#ifdef ULLong
179 y = *x * (ULLong)m + carry;
180 carry = y >> 32;
181 *x++ = y & 0xffffffffUL;
182#else
183#ifdef Pack_32
184 xi = *x;
185 y = (xi & 0xffff) * m + carry;
186 z = (xi >> 16) * m + (y >> 16);
187 carry = z >> 16;
188 *x++ = (z << 16) + (y & 0xffff);
189#else
190 y = *x * m + carry;
191 carry = y >> 16;
192 *x++ = y & 0xffff;
193#endif
194#endif
195 }
196 while(++i < wds);
197 if (carry) {
198 if (wds >= b->maxwds) {
199 b1 = Balloc(b->k+1);
200 Bcopy(b1, b);
201 Bfree(b);
202 b = b1;
203 }
204 b->x[wds++] = carry;
205 b->wds = wds;
206 }
207 return b;
208 }
209
210 int
211hi0bits_D2A
212#ifdef KR_headers
| 119
120 if (x & 7) {
121 if (x & 1)
122 return 0;
123 if (x & 2) {
124 *y = x >> 1;
125 return 1;
126 }
127 *y = x >> 2;
128 return 2;
129 }
130 k = 0;
131 if (!(x & 0xffff)) {
132 k = 16;
133 x >>= 16;
134 }
135 if (!(x & 0xff)) {
136 k += 8;
137 x >>= 8;
138 }
139 if (!(x & 0xf)) {
140 k += 4;
141 x >>= 4;
142 }
143 if (!(x & 0x3)) {
144 k += 2;
145 x >>= 2;
146 }
147 if (!(x & 1)) {
148 k++;
149 x >>= 1;
150 if (!x)
151 return 32;
152 }
153 *y = x;
154 return k;
155 }
156
157 Bigint *
158multadd
159#ifdef KR_headers
160 (b, m, a) Bigint *b; int m, a;
161#else
162 (Bigint *b, int m, int a) /* multiply by m and add a */
163#endif
164{
165 int i, wds;
166#ifdef ULLong
167 ULong *x;
168 ULLong carry, y;
169#else
170 ULong carry, *x, y;
171#ifdef Pack_32
172 ULong xi, z;
173#endif
174#endif
175 Bigint *b1;
176
177 wds = b->wds;
178 x = b->x;
179 i = 0;
180 carry = a;
181 do {
182#ifdef ULLong
183 y = *x * (ULLong)m + carry;
184 carry = y >> 32;
185 *x++ = y & 0xffffffffUL;
186#else
187#ifdef Pack_32
188 xi = *x;
189 y = (xi & 0xffff) * m + carry;
190 z = (xi >> 16) * m + (y >> 16);
191 carry = z >> 16;
192 *x++ = (z << 16) + (y & 0xffff);
193#else
194 y = *x * m + carry;
195 carry = y >> 16;
196 *x++ = y & 0xffff;
197#endif
198#endif
199 }
200 while(++i < wds);
201 if (carry) {
202 if (wds >= b->maxwds) {
203 b1 = Balloc(b->k+1);
204 Bcopy(b1, b);
205 Bfree(b);
206 b = b1;
207 }
208 b->x[wds++] = carry;
209 b->wds = wds;
210 }
211 return b;
212 }
213
214 int
215hi0bits_D2A
216#ifdef KR_headers
|
213 (x) register ULong x;
| 217 (x) ULong x;
|
214#else
| 218#else
|
215 (register ULong x)
| 219 (ULong x)
|
216#endif
217{
| 220#endif
221{
|
218 register int k = 0;
| 222 int k = 0;
|
219
220 if (!(x & 0xffff0000)) {
221 k = 16;
222 x <<= 16;
223 }
224 if (!(x & 0xff000000)) {
225 k += 8;
226 x <<= 8;
227 }
228 if (!(x & 0xf0000000)) {
229 k += 4;
230 x <<= 4;
231 }
232 if (!(x & 0xc0000000)) {
233 k += 2;
234 x <<= 2;
235 }
236 if (!(x & 0x80000000)) {
237 k++;
238 if (!(x & 0x40000000))
239 return 32;
240 }
241 return k;
242 }
243
244 Bigint *
245i2b
246#ifdef KR_headers
247 (i) int i;
248#else
249 (int i)
250#endif
251{
252 Bigint *b;
253
254 b = Balloc(1);
255 b->x[0] = i;
256 b->wds = 1;
257 return b;
258 }
259
260 Bigint *
261mult
262#ifdef KR_headers
263 (a, b) Bigint *a, *b;
264#else
265 (Bigint *a, Bigint *b)
266#endif
267{
268 Bigint *c;
269 int k, wa, wb, wc;
270 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
271 ULong y;
272#ifdef ULLong
273 ULLong carry, z;
274#else
275 ULong carry, z;
276#ifdef Pack_32
277 ULong z2;
278#endif
279#endif
280
281 if (a->wds < b->wds) {
282 c = a;
283 a = b;
284 b = c;
285 }
286 k = a->k;
287 wa = a->wds;
288 wb = b->wds;
289 wc = wa + wb;
290 if (wc > a->maxwds)
291 k++;
292 c = Balloc(k);
293 for(x = c->x, xa = x + wc; x < xa; x++)
294 *x = 0;
295 xa = a->x;
296 xae = xa + wa;
297 xb = b->x;
298 xbe = xb + wb;
299 xc0 = c->x;
300#ifdef ULLong
301 for(; xb < xbe; xc0++) {
302 if ( (y = *xb++) !=0) {
303 x = xa;
304 xc = xc0;
305 carry = 0;
306 do {
307 z = *x++ * (ULLong)y + *xc + carry;
308 carry = z >> 32;
309 *xc++ = z & 0xffffffffUL;
310 }
311 while(x < xae);
312 *xc = carry;
313 }
314 }
315#else
316#ifdef Pack_32
317 for(; xb < xbe; xb++, xc0++) {
318 if ( (y = *xb & 0xffff) !=0) {
319 x = xa;
320 xc = xc0;
321 carry = 0;
322 do {
323 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
324 carry = z >> 16;
325 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
326 carry = z2 >> 16;
327 Storeinc(xc, z2, z);
328 }
329 while(x < xae);
330 *xc = carry;
331 }
332 if ( (y = *xb >> 16) !=0) {
333 x = xa;
334 xc = xc0;
335 carry = 0;
336 z2 = *xc;
337 do {
338 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
339 carry = z >> 16;
340 Storeinc(xc, z, z2);
341 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
342 carry = z2 >> 16;
343 }
344 while(x < xae);
345 *xc = z2;
346 }
347 }
348#else
349 for(; xb < xbe; xc0++) {
350 if ( (y = *xb++) !=0) {
351 x = xa;
352 xc = xc0;
353 carry = 0;
354 do {
355 z = *x++ * y + *xc + carry;
356 carry = z >> 16;
357 *xc++ = z & 0xffff;
358 }
359 while(x < xae);
360 *xc = carry;
361 }
362 }
363#endif
364#endif
365 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
366 c->wds = wc;
367 return c;
368 }
369
370 static Bigint *p5s;
371
372 Bigint *
373pow5mult
374#ifdef KR_headers
375 (b, k) Bigint *b; int k;
376#else
377 (Bigint *b, int k)
378#endif
379{
380 Bigint *b1, *p5, *p51;
381 int i;
382 static int p05[3] = { 5, 25, 125 };
383
384 if ( (i = k & 3) !=0)
385 b = multadd(b, p05[i-1], 0);
386
387 if (!(k >>= 2))
388 return b;
389 if ((p5 = p5s) == 0) {
390 /* first time */
391#ifdef MULTIPLE_THREADS
392 ACQUIRE_DTOA_LOCK(1);
393 if (!(p5 = p5s)) {
394 p5 = p5s = i2b(625);
395 p5->next = 0;
396 }
397 FREE_DTOA_LOCK(1);
398#else
399 p5 = p5s = i2b(625);
400 p5->next = 0;
401#endif
402 }
403 for(;;) {
404 if (k & 1) {
405 b1 = mult(b, p5);
406 Bfree(b);
407 b = b1;
408 }
409 if (!(k >>= 1))
410 break;
411 if ((p51 = p5->next) == 0) {
412#ifdef MULTIPLE_THREADS
413 ACQUIRE_DTOA_LOCK(1);
414 if (!(p51 = p5->next)) {
415 p51 = p5->next = mult(p5,p5);
416 p51->next = 0;
417 }
418 FREE_DTOA_LOCK(1);
419#else
420 p51 = p5->next = mult(p5,p5);
421 p51->next = 0;
422#endif
423 }
424 p5 = p51;
425 }
426 return b;
427 }
428
429 Bigint *
430lshift
431#ifdef KR_headers
432 (b, k) Bigint *b; int k;
433#else
434 (Bigint *b, int k)
435#endif
436{
437 int i, k1, n, n1;
438 Bigint *b1;
439 ULong *x, *x1, *xe, z;
440
441 n = k >> kshift;
442 k1 = b->k;
443 n1 = n + b->wds + 1;
444 for(i = b->maxwds; n1 > i; i <<= 1)
445 k1++;
446 b1 = Balloc(k1);
447 x1 = b1->x;
448 for(i = 0; i < n; i++)
449 *x1++ = 0;
450 x = b->x;
451 xe = x + b->wds;
452 if (k &= kmask) {
453#ifdef Pack_32
454 k1 = 32 - k;
455 z = 0;
456 do {
457 *x1++ = *x << k | z;
458 z = *x++ >> k1;
459 }
460 while(x < xe);
461 if ((*x1 = z) !=0)
462 ++n1;
463#else
464 k1 = 16 - k;
465 z = 0;
466 do {
467 *x1++ = *x << k & 0xffff | z;
468 z = *x++ >> k1;
469 }
470 while(x < xe);
471 if (*x1 = z)
472 ++n1;
473#endif
474 }
475 else do
476 *x1++ = *x++;
477 while(x < xe);
478 b1->wds = n1 - 1;
479 Bfree(b);
480 return b1;
481 }
482
483 int
484cmp
485#ifdef KR_headers
486 (a, b) Bigint *a, *b;
487#else
488 (Bigint *a, Bigint *b)
489#endif
490{
491 ULong *xa, *xa0, *xb, *xb0;
492 int i, j;
493
494 i = a->wds;
495 j = b->wds;
496#ifdef DEBUG
497 if (i > 1 && !a->x[i-1])
498 Bug("cmp called with a->x[a->wds-1] == 0");
499 if (j > 1 && !b->x[j-1])
500 Bug("cmp called with b->x[b->wds-1] == 0");
501#endif
502 if (i -= j)
503 return i;
504 xa0 = a->x;
505 xa = xa0 + j;
506 xb0 = b->x;
507 xb = xb0 + j;
508 for(;;) {
509 if (*--xa != *--xb)
510 return *xa < *xb ? -1 : 1;
511 if (xa <= xa0)
512 break;
513 }
514 return 0;
515 }
516
517 Bigint *
518diff
519#ifdef KR_headers
520 (a, b) Bigint *a, *b;
521#else
522 (Bigint *a, Bigint *b)
523#endif
524{
525 Bigint *c;
526 int i, wa, wb;
527 ULong *xa, *xae, *xb, *xbe, *xc;
528#ifdef ULLong
529 ULLong borrow, y;
530#else
531 ULong borrow, y;
532#ifdef Pack_32
533 ULong z;
534#endif
535#endif
536
537 i = cmp(a,b);
538 if (!i) {
539 c = Balloc(0);
540 c->wds = 1;
541 c->x[0] = 0;
542 return c;
543 }
544 if (i < 0) {
545 c = a;
546 a = b;
547 b = c;
548 i = 1;
549 }
550 else
551 i = 0;
552 c = Balloc(a->k);
553 c->sign = i;
554 wa = a->wds;
555 xa = a->x;
556 xae = xa + wa;
557 wb = b->wds;
558 xb = b->x;
559 xbe = xb + wb;
560 xc = c->x;
561 borrow = 0;
562#ifdef ULLong
563 do {
564 y = (ULLong)*xa++ - *xb++ - borrow;
565 borrow = y >> 32 & 1UL;
566 *xc++ = y & 0xffffffffUL;
567 }
568 while(xb < xbe);
569 while(xa < xae) {
570 y = *xa++ - borrow;
571 borrow = y >> 32 & 1UL;
572 *xc++ = y & 0xffffffffUL;
573 }
574#else
575#ifdef Pack_32
576 do {
577 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
578 borrow = (y & 0x10000) >> 16;
579 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
580 borrow = (z & 0x10000) >> 16;
581 Storeinc(xc, z, y);
582 }
583 while(xb < xbe);
584 while(xa < xae) {
585 y = (*xa & 0xffff) - borrow;
586 borrow = (y & 0x10000) >> 16;
587 z = (*xa++ >> 16) - borrow;
588 borrow = (z & 0x10000) >> 16;
589 Storeinc(xc, z, y);
590 }
591#else
592 do {
593 y = *xa++ - *xb++ - borrow;
594 borrow = (y & 0x10000) >> 16;
595 *xc++ = y & 0xffff;
596 }
597 while(xb < xbe);
598 while(xa < xae) {
599 y = *xa++ - borrow;
600 borrow = (y & 0x10000) >> 16;
601 *xc++ = y & 0xffff;
602 }
603#endif
604#endif
605 while(!*--xc)
606 wa--;
607 c->wds = wa;
608 return c;
609 }
610
611 double
612b2d
613#ifdef KR_headers
614 (a, e) Bigint *a; int *e;
615#else
616 (Bigint *a, int *e)
617#endif
618{
619 ULong *xa, *xa0, w, y, z;
620 int k;
| 223
224 if (!(x & 0xffff0000)) {
225 k = 16;
226 x <<= 16;
227 }
228 if (!(x & 0xff000000)) {
229 k += 8;
230 x <<= 8;
231 }
232 if (!(x & 0xf0000000)) {
233 k += 4;
234 x <<= 4;
235 }
236 if (!(x & 0xc0000000)) {
237 k += 2;
238 x <<= 2;
239 }
240 if (!(x & 0x80000000)) {
241 k++;
242 if (!(x & 0x40000000))
243 return 32;
244 }
245 return k;
246 }
247
248 Bigint *
249i2b
250#ifdef KR_headers
251 (i) int i;
252#else
253 (int i)
254#endif
255{
256 Bigint *b;
257
258 b = Balloc(1);
259 b->x[0] = i;
260 b->wds = 1;
261 return b;
262 }
263
264 Bigint *
265mult
266#ifdef KR_headers
267 (a, b) Bigint *a, *b;
268#else
269 (Bigint *a, Bigint *b)
270#endif
271{
272 Bigint *c;
273 int k, wa, wb, wc;
274 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
275 ULong y;
276#ifdef ULLong
277 ULLong carry, z;
278#else
279 ULong carry, z;
280#ifdef Pack_32
281 ULong z2;
282#endif
283#endif
284
285 if (a->wds < b->wds) {
286 c = a;
287 a = b;
288 b = c;
289 }
290 k = a->k;
291 wa = a->wds;
292 wb = b->wds;
293 wc = wa + wb;
294 if (wc > a->maxwds)
295 k++;
296 c = Balloc(k);
297 for(x = c->x, xa = x + wc; x < xa; x++)
298 *x = 0;
299 xa = a->x;
300 xae = xa + wa;
301 xb = b->x;
302 xbe = xb + wb;
303 xc0 = c->x;
304#ifdef ULLong
305 for(; xb < xbe; xc0++) {
306 if ( (y = *xb++) !=0) {
307 x = xa;
308 xc = xc0;
309 carry = 0;
310 do {
311 z = *x++ * (ULLong)y + *xc + carry;
312 carry = z >> 32;
313 *xc++ = z & 0xffffffffUL;
314 }
315 while(x < xae);
316 *xc = carry;
317 }
318 }
319#else
320#ifdef Pack_32
321 for(; xb < xbe; xb++, xc0++) {
322 if ( (y = *xb & 0xffff) !=0) {
323 x = xa;
324 xc = xc0;
325 carry = 0;
326 do {
327 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
328 carry = z >> 16;
329 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
330 carry = z2 >> 16;
331 Storeinc(xc, z2, z);
332 }
333 while(x < xae);
334 *xc = carry;
335 }
336 if ( (y = *xb >> 16) !=0) {
337 x = xa;
338 xc = xc0;
339 carry = 0;
340 z2 = *xc;
341 do {
342 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
343 carry = z >> 16;
344 Storeinc(xc, z, z2);
345 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
346 carry = z2 >> 16;
347 }
348 while(x < xae);
349 *xc = z2;
350 }
351 }
352#else
353 for(; xb < xbe; xc0++) {
354 if ( (y = *xb++) !=0) {
355 x = xa;
356 xc = xc0;
357 carry = 0;
358 do {
359 z = *x++ * y + *xc + carry;
360 carry = z >> 16;
361 *xc++ = z & 0xffff;
362 }
363 while(x < xae);
364 *xc = carry;
365 }
366 }
367#endif
368#endif
369 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
370 c->wds = wc;
371 return c;
372 }
373
374 static Bigint *p5s;
375
376 Bigint *
377pow5mult
378#ifdef KR_headers
379 (b, k) Bigint *b; int k;
380#else
381 (Bigint *b, int k)
382#endif
383{
384 Bigint *b1, *p5, *p51;
385 int i;
386 static int p05[3] = { 5, 25, 125 };
387
388 if ( (i = k & 3) !=0)
389 b = multadd(b, p05[i-1], 0);
390
391 if (!(k >>= 2))
392 return b;
393 if ((p5 = p5s) == 0) {
394 /* first time */
395#ifdef MULTIPLE_THREADS
396 ACQUIRE_DTOA_LOCK(1);
397 if (!(p5 = p5s)) {
398 p5 = p5s = i2b(625);
399 p5->next = 0;
400 }
401 FREE_DTOA_LOCK(1);
402#else
403 p5 = p5s = i2b(625);
404 p5->next = 0;
405#endif
406 }
407 for(;;) {
408 if (k & 1) {
409 b1 = mult(b, p5);
410 Bfree(b);
411 b = b1;
412 }
413 if (!(k >>= 1))
414 break;
415 if ((p51 = p5->next) == 0) {
416#ifdef MULTIPLE_THREADS
417 ACQUIRE_DTOA_LOCK(1);
418 if (!(p51 = p5->next)) {
419 p51 = p5->next = mult(p5,p5);
420 p51->next = 0;
421 }
422 FREE_DTOA_LOCK(1);
423#else
424 p51 = p5->next = mult(p5,p5);
425 p51->next = 0;
426#endif
427 }
428 p5 = p51;
429 }
430 return b;
431 }
432
433 Bigint *
434lshift
435#ifdef KR_headers
436 (b, k) Bigint *b; int k;
437#else
438 (Bigint *b, int k)
439#endif
440{
441 int i, k1, n, n1;
442 Bigint *b1;
443 ULong *x, *x1, *xe, z;
444
445 n = k >> kshift;
446 k1 = b->k;
447 n1 = n + b->wds + 1;
448 for(i = b->maxwds; n1 > i; i <<= 1)
449 k1++;
450 b1 = Balloc(k1);
451 x1 = b1->x;
452 for(i = 0; i < n; i++)
453 *x1++ = 0;
454 x = b->x;
455 xe = x + b->wds;
456 if (k &= kmask) {
457#ifdef Pack_32
458 k1 = 32 - k;
459 z = 0;
460 do {
461 *x1++ = *x << k | z;
462 z = *x++ >> k1;
463 }
464 while(x < xe);
465 if ((*x1 = z) !=0)
466 ++n1;
467#else
468 k1 = 16 - k;
469 z = 0;
470 do {
471 *x1++ = *x << k & 0xffff | z;
472 z = *x++ >> k1;
473 }
474 while(x < xe);
475 if (*x1 = z)
476 ++n1;
477#endif
478 }
479 else do
480 *x1++ = *x++;
481 while(x < xe);
482 b1->wds = n1 - 1;
483 Bfree(b);
484 return b1;
485 }
486
487 int
488cmp
489#ifdef KR_headers
490 (a, b) Bigint *a, *b;
491#else
492 (Bigint *a, Bigint *b)
493#endif
494{
495 ULong *xa, *xa0, *xb, *xb0;
496 int i, j;
497
498 i = a->wds;
499 j = b->wds;
500#ifdef DEBUG
501 if (i > 1 && !a->x[i-1])
502 Bug("cmp called with a->x[a->wds-1] == 0");
503 if (j > 1 && !b->x[j-1])
504 Bug("cmp called with b->x[b->wds-1] == 0");
505#endif
506 if (i -= j)
507 return i;
508 xa0 = a->x;
509 xa = xa0 + j;
510 xb0 = b->x;
511 xb = xb0 + j;
512 for(;;) {
513 if (*--xa != *--xb)
514 return *xa < *xb ? -1 : 1;
515 if (xa <= xa0)
516 break;
517 }
518 return 0;
519 }
520
521 Bigint *
522diff
523#ifdef KR_headers
524 (a, b) Bigint *a, *b;
525#else
526 (Bigint *a, Bigint *b)
527#endif
528{
529 Bigint *c;
530 int i, wa, wb;
531 ULong *xa, *xae, *xb, *xbe, *xc;
532#ifdef ULLong
533 ULLong borrow, y;
534#else
535 ULong borrow, y;
536#ifdef Pack_32
537 ULong z;
538#endif
539#endif
540
541 i = cmp(a,b);
542 if (!i) {
543 c = Balloc(0);
544 c->wds = 1;
545 c->x[0] = 0;
546 return c;
547 }
548 if (i < 0) {
549 c = a;
550 a = b;
551 b = c;
552 i = 1;
553 }
554 else
555 i = 0;
556 c = Balloc(a->k);
557 c->sign = i;
558 wa = a->wds;
559 xa = a->x;
560 xae = xa + wa;
561 wb = b->wds;
562 xb = b->x;
563 xbe = xb + wb;
564 xc = c->x;
565 borrow = 0;
566#ifdef ULLong
567 do {
568 y = (ULLong)*xa++ - *xb++ - borrow;
569 borrow = y >> 32 & 1UL;
570 *xc++ = y & 0xffffffffUL;
571 }
572 while(xb < xbe);
573 while(xa < xae) {
574 y = *xa++ - borrow;
575 borrow = y >> 32 & 1UL;
576 *xc++ = y & 0xffffffffUL;
577 }
578#else
579#ifdef Pack_32
580 do {
581 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
582 borrow = (y & 0x10000) >> 16;
583 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
584 borrow = (z & 0x10000) >> 16;
585 Storeinc(xc, z, y);
586 }
587 while(xb < xbe);
588 while(xa < xae) {
589 y = (*xa & 0xffff) - borrow;
590 borrow = (y & 0x10000) >> 16;
591 z = (*xa++ >> 16) - borrow;
592 borrow = (z & 0x10000) >> 16;
593 Storeinc(xc, z, y);
594 }
595#else
596 do {
597 y = *xa++ - *xb++ - borrow;
598 borrow = (y & 0x10000) >> 16;
599 *xc++ = y & 0xffff;
600 }
601 while(xb < xbe);
602 while(xa < xae) {
603 y = *xa++ - borrow;
604 borrow = (y & 0x10000) >> 16;
605 *xc++ = y & 0xffff;
606 }
607#endif
608#endif
609 while(!*--xc)
610 wa--;
611 c->wds = wa;
612 return c;
613 }
614
615 double
616b2d
617#ifdef KR_headers
618 (a, e) Bigint *a; int *e;
619#else
620 (Bigint *a, int *e)
621#endif
622{
623 ULong *xa, *xa0, w, y, z;
624 int k;
|
621 double d;
| 625 U d;
|
622#ifdef VAX
623 ULong d0, d1;
624#else
| 626#ifdef VAX
627 ULong d0, d1;
628#else
|
625#define d0 word0(d)
626#define d1 word1(d)
| 629#define d0 word0(&d)
630#define d1 word1(&d)
|
627#endif
628
629 xa0 = a->x;
630 xa = xa0 + a->wds;
631 y = *--xa;
632#ifdef DEBUG
633 if (!y) Bug("zero y in b2d");
634#endif
635 k = hi0bits(y);
636 *e = 32 - k;
637#ifdef Pack_32
638 if (k < Ebits) {
| 631#endif
632
633 xa0 = a->x;
634 xa = xa0 + a->wds;
635 y = *--xa;
636#ifdef DEBUG
637 if (!y) Bug("zero y in b2d");
638#endif
639 k = hi0bits(y);
640 *e = 32 - k;
641#ifdef Pack_32
642 if (k < Ebits) {
|
639 d0 = Exp_1 | y >> Ebits - k;
| 643 d0 = Exp_1 | y >> (Ebits - k);
|
640 w = xa > xa0 ? *--xa : 0;
| 644 w = xa > xa0 ? *--xa : 0;
|
641 d1 = y << (32-Ebits) + k | w >> Ebits - k;
| 645 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
|
642 goto ret_d;
643 }
644 z = xa > xa0 ? *--xa : 0;
645 if (k -= Ebits) {
| 646 goto ret_d;
647 }
648 z = xa > xa0 ? *--xa : 0;
649 if (k -= Ebits) {
|
646 d0 = Exp_1 | y << k | z >> 32 - k;
| 650 d0 = Exp_1 | y << k | z >> (32 - k);
|
647 y = xa > xa0 ? *--xa : 0;
| 651 y = xa > xa0 ? *--xa : 0;
|
648 d1 = z << k | y >> 32 - k;
| 652 d1 = z << k | y >> (32 - k);
|
649 }
650 else {
651 d0 = Exp_1 | y;
652 d1 = z;
653 }
654#else
655 if (k < Ebits + 16) {
656 z = xa > xa0 ? *--xa : 0;
657 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
658 w = xa > xa0 ? *--xa : 0;
659 y = xa > xa0 ? *--xa : 0;
660 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
661 goto ret_d;
662 }
663 z = xa > xa0 ? *--xa : 0;
664 w = xa > xa0 ? *--xa : 0;
665 k -= Ebits + 16;
666 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
667 y = xa > xa0 ? *--xa : 0;
668 d1 = w << k + 16 | y << k;
669#endif
670 ret_d:
671#ifdef VAX
| 653 }
654 else {
655 d0 = Exp_1 | y;
656 d1 = z;
657 }
658#else
659 if (k < Ebits + 16) {
660 z = xa > xa0 ? *--xa : 0;
661 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
662 w = xa > xa0 ? *--xa : 0;
663 y = xa > xa0 ? *--xa : 0;
664 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
665 goto ret_d;
666 }
667 z = xa > xa0 ? *--xa : 0;
668 w = xa > xa0 ? *--xa : 0;
669 k -= Ebits + 16;
670 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
671 y = xa > xa0 ? *--xa : 0;
672 d1 = w << k + 16 | y << k;
673#endif
674 ret_d:
675#ifdef VAX
|
672 word0(d) = d0 >> 16 | d0 << 16;
673 word1(d) = d1 >> 16 | d1 << 16;
| 676 word0(&d) = d0 >> 16 | d0 << 16;
677 word1(&d) = d1 >> 16 | d1 << 16;
|
674#endif
| 678#endif
|
675 return dval(d);
| 679 return dval(&d);
|
676 }
677#undef d0
678#undef d1
679
680 Bigint *
681d2b
682#ifdef KR_headers
| 680 }
681#undef d0
682#undef d1
683
684 Bigint *
685d2b
686#ifdef KR_headers
|
683 (d, e, bits) double d; int *e, *bits;
| 687 (dd, e, bits) double dd; int *e, *bits;
|
684#else
| 688#else
|
685 (double d, int *e, int *bits)
| 689 (double dd, int *e, int *bits)
|
686#endif
687{
688 Bigint *b;
| 690#endif
691{
692 Bigint *b;
|
| 693 U d;
|
689#ifndef Sudden_Underflow
690 int i;
691#endif
692 int de, k;
693 ULong *x, y, z;
694#ifdef VAX
695 ULong d0, d1;
| 694#ifndef Sudden_Underflow
695 int i;
696#endif
697 int de, k;
698 ULong *x, y, z;
699#ifdef VAX
700 ULong d0, d1;
|
696 d0 = word0(d) >> 16 | word0(d) << 16;
697 d1 = word1(d) >> 16 | word1(d) << 16;
| |
698#else
| 701#else
|
699#define d0 word0(d)
700#define d1 word1(d)
| 702#define d0 word0(&d)
703#define d1 word1(&d)
|
701#endif
| 704#endif
|
| 705 d.d = dd;
706#ifdef VAX
707 d0 = word0(&d) >> 16 | word0(&d) << 16;
708 d1 = word1(&d) >> 16 | word1(&d) << 16;
709#endif
|
702
703#ifdef Pack_32
704 b = Balloc(1);
705#else
706 b = Balloc(2);
707#endif
708 x = b->x;
709
710 z = d0 & Frac_mask;
711 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
712#ifdef Sudden_Underflow
713 de = (int)(d0 >> Exp_shift);
714#ifndef IBM
715 z |= Exp_msk11;
716#endif
717#else
718 if ( (de = (int)(d0 >> Exp_shift)) !=0)
719 z |= Exp_msk1;
720#endif
721#ifdef Pack_32
722 if ( (y = d1) !=0) {
723 if ( (k = lo0bits(&y)) !=0) {
| 710
711#ifdef Pack_32
712 b = Balloc(1);
713#else
714 b = Balloc(2);
715#endif
716 x = b->x;
717
718 z = d0 & Frac_mask;
719 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
720#ifdef Sudden_Underflow
721 de = (int)(d0 >> Exp_shift);
722#ifndef IBM
723 z |= Exp_msk11;
724#endif
725#else
726 if ( (de = (int)(d0 >> Exp_shift)) !=0)
727 z |= Exp_msk1;
728#endif
729#ifdef Pack_32
730 if ( (y = d1) !=0) {
731 if ( (k = lo0bits(&y)) !=0) {
|
724 x[0] = y | z << 32 - k;
| 732 x[0] = y | z << (32 - k);
|
725 z >>= k;
726 }
727 else
728 x[0] = y;
729#ifndef Sudden_Underflow
730 i =
731#endif
732 b->wds = (x[1] = z) !=0 ? 2 : 1;
733 }
734 else {
| 733 z >>= k;
734 }
735 else
736 x[0] = y;
737#ifndef Sudden_Underflow
738 i =
739#endif
740 b->wds = (x[1] = z) !=0 ? 2 : 1;
741 }
742 else {
|
735#ifdef DEBUG
736 if (!z)
737 Bug("Zero passed to d2b");
738#endif
| |
739 k = lo0bits(&z);
740 x[0] = z;
741#ifndef Sudden_Underflow
742 i =
743#endif
744 b->wds = 1;
745 k += 32;
746 }
747#else
748 if ( (y = d1) !=0) {
749 if ( (k = lo0bits(&y)) !=0)
750 if (k >= 16) {
751 x[0] = y | z << 32 - k & 0xffff;
752 x[1] = z >> k - 16 & 0xffff;
753 x[2] = z >> k;
754 i = 2;
755 }
756 else {
757 x[0] = y & 0xffff;
758 x[1] = y >> 16 | z << 16 - k & 0xffff;
759 x[2] = z >> k & 0xffff;
760 x[3] = z >> k+16;
761 i = 3;
762 }
763 else {
764 x[0] = y & 0xffff;
765 x[1] = y >> 16;
766 x[2] = z & 0xffff;
767 x[3] = z >> 16;
768 i = 3;
769 }
770 }
771 else {
772#ifdef DEBUG
773 if (!z)
774 Bug("Zero passed to d2b");
775#endif
776 k = lo0bits(&z);
777 if (k >= 16) {
778 x[0] = z;
779 i = 0;
780 }
781 else {
782 x[0] = z & 0xffff;
783 x[1] = z >> 16;
784 i = 1;
785 }
786 k += 32;
787 }
788 while(!x[i])
789 --i;
790 b->wds = i + 1;
791#endif
792#ifndef Sudden_Underflow
793 if (de) {
794#endif
795#ifdef IBM
796 *e = (de - Bias - (P-1) << 2) + k;
| 743 k = lo0bits(&z);
744 x[0] = z;
745#ifndef Sudden_Underflow
746 i =
747#endif
748 b->wds = 1;
749 k += 32;
750 }
751#else
752 if ( (y = d1) !=0) {
753 if ( (k = lo0bits(&y)) !=0)
754 if (k >= 16) {
755 x[0] = y | z << 32 - k & 0xffff;
756 x[1] = z >> k - 16 & 0xffff;
757 x[2] = z >> k;
758 i = 2;
759 }
760 else {
761 x[0] = y & 0xffff;
762 x[1] = y >> 16 | z << 16 - k & 0xffff;
763 x[2] = z >> k & 0xffff;
764 x[3] = z >> k+16;
765 i = 3;
766 }
767 else {
768 x[0] = y & 0xffff;
769 x[1] = y >> 16;
770 x[2] = z & 0xffff;
771 x[3] = z >> 16;
772 i = 3;
773 }
774 }
775 else {
776#ifdef DEBUG
777 if (!z)
778 Bug("Zero passed to d2b");
779#endif
780 k = lo0bits(&z);
781 if (k >= 16) {
782 x[0] = z;
783 i = 0;
784 }
785 else {
786 x[0] = z & 0xffff;
787 x[1] = z >> 16;
788 i = 1;
789 }
790 k += 32;
791 }
792 while(!x[i])
793 --i;
794 b->wds = i + 1;
795#endif
796#ifndef Sudden_Underflow
797 if (de) {
798#endif
799#ifdef IBM
800 *e = (de - Bias - (P-1) << 2) + k;
|
797 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
| 801 *bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask);
|
798#else
799 *e = de - Bias - (P-1) + k;
800 *bits = P - k;
801#endif
802#ifndef Sudden_Underflow
803 }
804 else {
805 *e = de - Bias - (P-1) + 1 + k;
806#ifdef Pack_32
807 *bits = 32*i - hi0bits(x[i-1]);
808#else
809 *bits = (i+2)*16 - hi0bits(x[i]);
810#endif
811 }
812#endif
813 return b;
814 }
815#undef d0
816#undef d1
817
818 CONST double
819#ifdef IEEE_Arith
820bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
821CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
822 };
823#else
824#ifdef IBM
825bigtens[] = { 1e16, 1e32, 1e64 };
826CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
827#else
828bigtens[] = { 1e16, 1e32 };
829CONST double tinytens[] = { 1e-16, 1e-32 };
830#endif
831#endif
832
833 CONST double
834tens[] = {
835 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
836 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
837 1e20, 1e21, 1e22
838#ifdef VAX
839 , 1e23, 1e24
840#endif
841 };
842
843 char *
844#ifdef KR_headers
845strcp_D2A(a, b) char *a; char *b;
846#else
847strcp_D2A(char *a, CONST char *b)
848#endif
849{
| 802#else
803 *e = de - Bias - (P-1) + k;
804 *bits = P - k;
805#endif
806#ifndef Sudden_Underflow
807 }
808 else {
809 *e = de - Bias - (P-1) + 1 + k;
810#ifdef Pack_32
811 *bits = 32*i - hi0bits(x[i-1]);
812#else
813 *bits = (i+2)*16 - hi0bits(x[i]);
814#endif
815 }
816#endif
817 return b;
818 }
819#undef d0
820#undef d1
821
822 CONST double
823#ifdef IEEE_Arith
824bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
825CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
826 };
827#else
828#ifdef IBM
829bigtens[] = { 1e16, 1e32, 1e64 };
830CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
831#else
832bigtens[] = { 1e16, 1e32 };
833CONST double tinytens[] = { 1e-16, 1e-32 };
834#endif
835#endif
836
837 CONST double
838tens[] = {
839 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
840 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
841 1e20, 1e21, 1e22
842#ifdef VAX
843 , 1e23, 1e24
844#endif
845 };
846
847 char *
848#ifdef KR_headers
849strcp_D2A(a, b) char *a; char *b;
850#else
851strcp_D2A(char *a, CONST char *b)
852#endif
853{
|
850 while(*a = *b++)
| 854 while((*a = *b++))
|
851 a++;
852 return a;
853 }
854
855#ifdef NO_STRING_H
856
857 Char *
858#ifdef KR_headers
859memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
860#else
861memcpy_D2A(void *a1, void *b1, size_t len)
862#endif
863{
| 855 a++;
856 return a;
857 }
858
859#ifdef NO_STRING_H
860
861 Char *
862#ifdef KR_headers
863memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
864#else
865memcpy_D2A(void *a1, void *b1, size_t len)
866#endif
867{
|
864 register char *a = (char*)a1, *ae = a + len;
865 register char *b = (char*)b1, *a0 = a;
| 868 char *a = (char*)a1, *ae = a + len;
869 char *b = (char*)b1, *a0 = a;
|
866 while(a < ae)
867 *a++ = *b++;
868 return a0;
869 }
870
871#endif /* NO_STRING_H */
| 870 while(a < ae)
871 *a++ = *b++;
872 return a0;
873 }
874
875#endif /* NO_STRING_H */
|