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