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