1 /****************************************************************
3 The author of this software is David M. Gay.
5 Copyright (C) 1998, 1999 by Lucent Technologies
8 Permission to use, copy, modify, and distribute this software and
9 its documentation for any purpose and without fee is hereby
10 granted, provided that the above copyright notice appear in all
11 copies and that both that the copyright notice and this
12 permission notice and warranty disclaimer appear in supporting
13 documentation, and that the name of Lucent or any of its entities
14 not be used in advertising or publicity pertaining to
15 distribution of the software without specific, written prior
18 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
27 ****************************************************************/
29 /* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
34 static Bigint
*freelist
[Kmax
+1];
35 #ifndef Omit_Private_Memory
37 #define PRIVATE_MEM 2304
39 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
40 static double private_mem
[PRIVATE_mem
], *pmem_next
= private_mem
;
53 #ifndef Omit_Private_Memory
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
;
65 #ifdef Omit_Private_Memory
66 rv
= (Bigint
*)MALLOC(sizeof(Bigint
) + (x
-1)*sizeof(ULong
));
68 len
= (sizeof(Bigint
) + (x
-1)*sizeof(ULong
) + sizeof(double) - 1)
70 if (k
<= Kmax
&& pmem_next
- private_mem
+ len
<= PRIVATE_mem
) {
71 rv
= (Bigint
*)pmem_next
;
75 rv
= (Bigint
*)MALLOC(len
*sizeof(double));
81 rv
->sign
= rv
->wds
= 0;
101 ACQUIRE_DTOA_LOCK(0);
102 v
->next
= freelist
[v
->k
];
160 (b
, m
, a
) Bigint
*b
; int m
, a
;
162 (Bigint
*b
, int m
, int a
) /* multiply by m and add a */
183 y
= *x
* (ULLong
)m
+ carry
;
185 *x
++ = y
& 0xffffffffUL
;
189 y
= (xi
& 0xffff) * m
+ carry
;
190 z
= (xi
>> 16) * m
+ (y
>> 16);
192 *x
++ = (z
<< 16) + (y
& 0xffff);
202 if (wds
>= b
->maxwds
) {
224 if (!(x
& 0xffff0000)) {
228 if (!(x
& 0xff000000)) {
232 if (!(x
& 0xf0000000)) {
236 if (!(x
& 0xc0000000)) {
240 if (!(x
& 0x80000000)) {
242 if (!(x
& 0x40000000))
267 (a
, b
) Bigint
*a
, *b
;
269 (Bigint
*a
, Bigint
*b
)
274 ULong
*x
, *xa
, *xae
, *xb
, *xbe
, *xc
, *xc0
;
285 if (a
->wds
< b
->wds
) {
297 for(x
= c
->x
, xa
= x
+ wc
; x
< xa
; x
++)
305 for(; xb
< xbe
; xc0
++) {
306 if ( (y
= *xb
++) !=0) {
311 z
= *x
++ * (ULLong
)y
+ *xc
+ carry
;
313 *xc
++ = z
& 0xffffffffUL
;
321 for(; xb
< xbe
; xb
++, xc0
++) {
322 if ( (y
= *xb
& 0xffff) !=0) {
327 z
= (*x
& 0xffff) * y
+ (*xc
& 0xffff) + carry
;
329 z2
= (*x
++ >> 16) * y
+ (*xc
>> 16) + carry
;
336 if ( (y
= *xb
>> 16) !=0) {
342 z
= (*x
& 0xffff) * y
+ (*xc
>> 16) + carry
;
345 z2
= (*x
++ >> 16) * y
+ (*xc
& 0xffff) + carry
;
353 for(; xb
< xbe
; xc0
++) {
354 if ( (y
= *xb
++) !=0) {
359 z
= *x
++ * y
+ *xc
+ carry
;
369 for(xc0
= c
->x
, xc
= xc0
+ wc
; wc
> 0 && !*--xc
; --wc
) ;
379 (b
, k
) Bigint
*b
; int k
;
384 Bigint
*b1
, *p5
, *p51
;
386 static int p05
[3] = { 5, 25, 125 };
388 if ( (i
= k
& 3) !=0)
389 b
= multadd(b
, p05
[i
-1], 0);
393 if ((p5
= p5s
) == 0) {
395 #ifdef MULTIPLE_THREADS
396 ACQUIRE_DTOA_LOCK(1);
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
);
424 p51
= p5
->next
= mult(p5
,p5
);
436 (b
, k
) Bigint
*b
; int k
;
443 ULong
*x
, *x1
, *xe
, z
;
448 for(i
= b
->maxwds
; n1
> i
; i
<<= 1)
452 for(i
= 0; i
< n
; i
++)
471 *x1
++ = *x
<< k
& 0xffff | z
;
490 (a
, b
) Bigint
*a
, *b
;
492 (Bigint
*a
, Bigint
*b
)
495 ULong
*xa
, *xa0
, *xb
, *xb0
;
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");
514 return *xa
< *xb
? -1 : 1;
524 (a
, b
) Bigint
*a
, *b
;
526 (Bigint
*a
, Bigint
*b
)
531 ULong
*xa
, *xae
, *xb
, *xbe
, *xc
;
568 y
= (ULLong
)*xa
++ - *xb
++ - borrow
;
569 borrow
= y
>> 32 & 1UL;
570 *xc
++ = y
& 0xffffffffUL
;
575 borrow
= y
>> 32 & 1UL;
576 *xc
++ = y
& 0xffffffffUL
;
581 y
= (*xa
& 0xffff) - (*xb
& 0xffff) - borrow
;
582 borrow
= (y
& 0x10000) >> 16;
583 z
= (*xa
++ >> 16) - (*xb
++ >> 16) - borrow
;
584 borrow
= (z
& 0x10000) >> 16;
589 y
= (*xa
& 0xffff) - borrow
;
590 borrow
= (y
& 0x10000) >> 16;
591 z
= (*xa
++ >> 16) - borrow
;
592 borrow
= (z
& 0x10000) >> 16;
597 y
= *xa
++ - *xb
++ - borrow
;
598 borrow
= (y
& 0x10000) >> 16;
604 borrow
= (y
& 0x10000) >> 16;
618 (a
, e
) Bigint
*a
; int *e
;
623 ULong
*xa
, *xa0
, w
, y
, z
;
637 if (!y
) Bug("zero y in b2d");
643 d0
= Exp_1
| y
>> (Ebits
- k
);
644 w
= xa
> xa0
? *--xa
: 0;
645 d1
= y
<< ((32-Ebits
) + k
) | w
>> (Ebits
- k
);
648 z
= xa
> xa0
? *--xa
: 0;
650 d0
= Exp_1
| y
<< k
| z
>> (32 - k
);
651 y
= xa
> xa0
? *--xa
: 0;
652 d1
= z
<< k
| y
>> (32 - k
);
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
;
667 z
= xa
> xa0
? *--xa
: 0;
668 w
= xa
> xa0
? *--xa
: 0;
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
;
676 word0(&d
) = d0
>> 16 | d0
<< 16;
677 word1(&d
) = d1
>> 16 | d1
<< 16;
687 (dd
, e
, bits
) double dd
; int *e
, *bits
;
689 (double dd
, int *e
, int *bits
)
694 #ifndef Sudden_Underflow
707 d0
= word0(&d
) >> 16 | word0(&d
) << 16;
708 d1
= word1(&d
) >> 16 | word1(&d
) << 16;
719 d0
&= 0x7fffffff; /* clear sign bit, which we ignore */
720 #ifdef Sudden_Underflow
721 de
= (int)(d0
>> Exp_shift
);
726 if ( (de
= (int)(d0
>> Exp_shift
)) !=0)
731 if ( (k
= lo0bits(&y
)) !=0) {
732 x
[0] = y
| z
<< (32 - k
);
737 #ifndef Sudden_Underflow
740 b
->wds
= (x
[1] = z
) !=0 ? 2 : 1;
745 #ifndef Sudden_Underflow
753 if ( (k
= lo0bits(&y
)) !=0)
755 x
[0] = y
| z
<< 32 - k
& 0xffff;
756 x
[1] = z
>> k
- 16 & 0xffff;
762 x
[1] = y
>> 16 | z
<< 16 - k
& 0xffff;
763 x
[2] = z
>> k
& 0xffff;
778 Bug("Zero passed to d2b");
796 #ifndef Sudden_Underflow
800 *e
= (de
- Bias
- (P
-1) << 2) + k
;
801 *bits
= 4*P
+ 8 - k
- hi0bits(word0(&d
) & Frac_mask
);
803 *e
= de
- Bias
- (P
-1) + k
;
806 #ifndef Sudden_Underflow
809 *e
= de
- Bias
- (P
-1) + 1 + k
;
811 *bits
= 32*i
- hi0bits(x
[i
-1]);
813 *bits
= (i
+2)*16 - hi0bits(x
[i
]);
824 bigtens
[] = { 1e16
, 1e32
, 1e64
, 1e128
, 1e256
};
825 CONST
double tinytens
[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
829 bigtens
[] = { 1e16
, 1e32
, 1e64
};
830 CONST
double tinytens
[] = { 1e-16, 1e-32, 1e-64 };
832 bigtens
[] = { 1e16
, 1e32
};
833 CONST
double tinytens
[] = { 1e-16, 1e-32 };
839 1e0
, 1e1
, 1e2
, 1e3
, 1e4
, 1e5
, 1e6
, 1e7
, 1e8
, 1e9
,
840 1e10
, 1e11
, 1e12
, 1e13
, 1e14
, 1e15
, 1e16
, 1e17
, 1e18
, 1e19
,
849 strcp_D2A(a
, b
) char *a
; char *b
;
851 strcp_D2A(char *a
, CONST
char *b
)
863 memcpy_D2A(a
, b
, len
) Char
*a
; Char
*b
; size_t len
;
865 memcpy_D2A(void *a1
, void *b1
, size_t len
)
868 char *a
= (char*)a1
, *ae
= a
+ len
;
869 char *b
= (char*)b1
, *a0
= a
;
875 #endif /* NO_STRING_H */