source: vendor/FreeBSD-libc/current/gdtoa/misc.c

Last change on this file was 1862, checked in by bird, 20 years ago

Initial revision
  • Property cvs2svn:cvs-rev set to 1.1
  • Property svn:eol-style set to native
  • Property svn:executable set to *
File size: 13.8 KB
Line 
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
30 David M. Gay
31 Bell Laboratories, Room 2C-463
32 600 Mountain Avenue
33 Murray Hill, NJ 07974-0636
34 U.S.A.
35 dmg@bell-labs.com
36 */
37
38#include "gdtoaimp.h"
39
40 static Bigint *freelist[Kmax+1];
41#ifndef Omit_Private_Memory
42#ifndef PRIVATE_MEM
43#define PRIVATE_MEM 2304
44#endif
45#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
46static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
47#endif
48
49 Bigint *
50Balloc
51#ifdef KR_headers
52 (k) int k;
53#else
54 (int k)
55#endif
56{
57 int x;
58 Bigint *rv;
59#ifndef Omit_Private_Memory
60 unsigned int len;
61#endif
62
63 ACQUIRE_DTOA_LOCK(0);
64 if ( (rv = freelist[k]) !=0) {
65 freelist[k] = rv->next;
66 }
67 else {
68 x = 1 << k;
69#ifdef Omit_Private_Memory
70 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
71#else
72 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
73 /sizeof(double);
74 if (pmem_next - private_mem + len <= PRIVATE_mem) {
75 rv = (Bigint*)pmem_next;
76 pmem_next += len;
77 }
78 else
79 rv = (Bigint*)MALLOC(len*sizeof(double));
80#endif
81 rv->k = k;
82 rv->maxwds = x;
83 }
84 FREE_DTOA_LOCK(0);
85 rv->sign = rv->wds = 0;
86 return rv;
87 }
88
89 void
90Bfree
91#ifdef KR_headers
92 (v) Bigint *v;
93#else
94 (Bigint *v)
95#endif
96{
97 if (v) {
98 ACQUIRE_DTOA_LOCK(0);
99 v->next = freelist[v->k];
100 freelist[v->k] = v;
101 FREE_DTOA_LOCK(0);
102 }
103 }
104
105 int
106lo0bits
107#ifdef KR_headers
108 (y) ULong *y;
109#else
110 (ULong *y)
111#endif
112{
113 register int k;
114 register 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 & 1)
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
212#ifdef KR_headers
213 (x) register ULong x;
214#else
215 (register ULong x)
216#endif
217{
218 register 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;
621 double d;
622#ifdef VAX
623 ULong d0, d1;
624#else
625#define d0 word0(d)
626#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) {
639 d0 = Exp_1 | y >> Ebits - k;
640 w = xa > xa0 ? *--xa : 0;
641 d1 = y << (32-Ebits) + k | w >> Ebits - k;
642 goto ret_d;
643 }
644 z = xa > xa0 ? *--xa : 0;
645 if (k -= Ebits) {
646 d0 = Exp_1 | y << k | z >> 32 - k;
647 y = xa > xa0 ? *--xa : 0;
648 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
672 word0(d) = d0 >> 16 | d0 << 16;
673 word1(d) = d1 >> 16 | d1 << 16;
674#endif
675 return dval(d);
676 }
677#undef d0
678#undef d1
679
680 Bigint *
681d2b
682#ifdef KR_headers
683 (d, e, bits) double d; int *e, *bits;
684#else
685 (double d, int *e, int *bits)
686#endif
687{
688 Bigint *b;
689 int de, i, k;
690 ULong *x, y, z;
691#ifdef VAX
692 ULong d0, d1;
693 d0 = word0(d) >> 16 | word0(d) << 16;
694 d1 = word1(d) >> 16 | word1(d) << 16;
695#else
696#define d0 word0(d)
697#define d1 word1(d)
698#endif
699
700#ifdef Pack_32
701 b = Balloc(1);
702#else
703 b = Balloc(2);
704#endif
705 x = b->x;
706
707 z = d0 & Frac_mask;
708 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
709#ifdef Sudden_Underflow
710 de = (int)(d0 >> Exp_shift);
711#ifndef IBM
712 z |= Exp_msk11;
713#endif
714#else
715 if ( (de = (int)(d0 >> Exp_shift)) !=0)
716 z |= Exp_msk1;
717#endif
718#ifdef Pack_32
719 if ( (y = d1) !=0) {
720 if ( (k = lo0bits(&y)) !=0) {
721 x[0] = y | z << 32 - k;
722 z >>= k;
723 }
724 else
725 x[0] = y;
726 i = 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 i = b->wds = 1;
736 k += 32;
737 }
738#else
739 if ( (y = d1) !=0) {
740 if ( (k = lo0bits(&y)) !=0)
741 if (k >= 16) {
742 x[0] = y | z << 32 - k & 0xffff;
743 x[1] = z >> k - 16 & 0xffff;
744 x[2] = z >> k;
745 i = 2;
746 }
747 else {
748 x[0] = y & 0xffff;
749 x[1] = y >> 16 | z << 16 - k & 0xffff;
750 x[2] = z >> k & 0xffff;
751 x[3] = z >> k+16;
752 i = 3;
753 }
754 else {
755 x[0] = y & 0xffff;
756 x[1] = y >> 16;
757 x[2] = z & 0xffff;
758 x[3] = z >> 16;
759 i = 3;
760 }
761 }
762 else {
763#ifdef DEBUG
764 if (!z)
765 Bug("Zero passed to d2b");
766#endif
767 k = lo0bits(&z);
768 if (k >= 16) {
769 x[0] = z;
770 i = 0;
771 }
772 else {
773 x[0] = z & 0xffff;
774 x[1] = z >> 16;
775 i = 1;
776 }
777 k += 32;
778 }
779 while(!x[i])
780 --i;
781 b->wds = i + 1;
782#endif
783#ifndef Sudden_Underflow
784 if (de) {
785#endif
786#ifdef IBM
787 *e = (de - Bias - (P-1) << 2) + k;
788 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
789#else
790 *e = de - Bias - (P-1) + k;
791 *bits = P - k;
792#endif
793#ifndef Sudden_Underflow
794 }
795 else {
796 *e = de - Bias - (P-1) + 1 + k;
797#ifdef Pack_32
798 *bits = 32*i - hi0bits(x[i-1]);
799#else
800 *bits = (i+2)*16 - hi0bits(x[i]);
801#endif
802 }
803#endif
804 return b;
805 }
806#undef d0
807#undef d1
808
809 CONST double
810#ifdef IEEE_Arith
811bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
812CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
813 };
814#else
815#ifdef IBM
816bigtens[] = { 1e16, 1e32, 1e64 };
817CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
818#else
819bigtens[] = { 1e16, 1e32 };
820CONST double tinytens[] = { 1e-16, 1e-32 };
821#endif
822#endif
823
824 CONST double
825tens[] = {
826 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
827 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
828 1e20, 1e21, 1e22
829#ifdef VAX
830 , 1e23, 1e24
831#endif
832 };
833
834 char *
835#ifdef KR_headers
836strcp_D2A(a, b) char *a; char *b;
837#else
838strcp_D2A(char *a, CONST char *b)
839#endif
840{
841 while(*a = *b++)
842 a++;
843 return a;
844 }
845
846#ifdef NO_STRING_H
847
848 Char *
849#ifdef KR_headers
850memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
851#else
852memcpy_D2A(void *a1, void *b1, size_t len)
853#endif
854{
855 register char *a = (char*)a1, *ae = a + len;
856 register char *b = (char*)b1, *a0 = a;
857 while(a < ae)
858 *a++ = *b++;
859 return a0;
860 }
861
862#endif /* NO_STRING_H */
Note: See TracBrowser for help on using the repository browser.