1 /****************************************************************
2
3 The author of this software is David M. Gay.
4
5 Copyright (C) 1998 by Lucent Technologies
6 All Rights Reserved
7
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
16 permission.
17
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
25 THIS 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 #ifndef MULTIPLE_THREADS
35 char *dtoa_result;
36 #endif
37
38 char *
39 #ifdef KR_headers
rv_alloc(i)40 rv_alloc(i) int i;
41 #else
42 rv_alloc(int i)
43 #endif
44 {
45 int j, k, *r;
46
47 j = sizeof(ULong);
48 for(k = 0;
49 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
50 j <<= 1)
51 k++;
52 r = (int*)Balloc(k);
53 if (r == NULL)
54 return (
55 #ifndef MULTIPLE_THREADS
56 dtoa_result =
57 #endif
58 NULL);
59 *r = k;
60 return
61 #ifndef MULTIPLE_THREADS
62 dtoa_result =
63 #endif
64 (char *)(r+1);
65 }
66
67 char *
68 #ifdef KR_headers
nrv_alloc(s,rve,n)69 nrv_alloc(s, rve, n) char *s, **rve; int n;
70 #else
71 nrv_alloc(char *s, char **rve, int n)
72 #endif
73 {
74 char *rv, *t;
75
76 t = rv = rv_alloc(n);
77 if (t == NULL)
78 return (NULL);
79 while((*t = *s++) !=0)
80 t++;
81 if (rve)
82 *rve = t;
83 return rv;
84 }
85
86 /* freedtoa(s) must be used to free values s returned by dtoa
87 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
88 * but for consistency with earlier versions of dtoa, it is optional
89 * when MULTIPLE_THREADS is not defined.
90 */
91
92 void
93 #ifdef KR_headers
freedtoa(s)94 freedtoa(s) char *s;
95 #else
96 freedtoa(char *s)
97 #endif
98 {
99 Bigint *b = (Bigint *)((int *)s - 1);
100 b->maxwds = 1 << (b->k = *(int*)b);
101 Bfree(b);
102 #ifndef MULTIPLE_THREADS
103 if (s == dtoa_result)
104 dtoa_result = 0;
105 #endif
106 }
107
108 int
quorem(b,S)109 quorem
110 #ifdef KR_headers
111 (b, S) Bigint *b, *S;
112 #else
113 (Bigint *b, Bigint *S)
114 #endif
115 {
116 int n;
117 ULong *bx, *bxe, q, *sx, *sxe;
118 #ifdef ULLong
119 ULLong borrow, carry, y, ys;
120 #else
121 ULong borrow, carry, y, ys;
122 #ifdef Pack_32
123 ULong si, z, zs;
124 #endif
125 #endif
126
127 n = S->wds;
128 #ifdef DEBUG
129 /*debug*/ if (b->wds > n)
130 /*debug*/ Bug("oversize b in quorem");
131 #endif
132 if (b->wds < n)
133 return 0;
134 sx = S->x;
135 sxe = sx + --n;
136 bx = b->x;
137 bxe = bx + n;
138 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
139 #ifdef DEBUG
140 /*debug*/ if (q > 9)
141 /*debug*/ Bug("oversized quotient in quorem");
142 #endif
143 if (q) {
144 borrow = 0;
145 carry = 0;
146 do {
147 #ifdef ULLong
148 ys = *sx++ * (ULLong)q + carry;
149 carry = ys >> 32;
150 y = *bx - (ys & 0xffffffffUL) - borrow;
151 borrow = y >> 32 & 1UL;
152 *bx++ = y & 0xffffffffUL;
153 #else
154 #ifdef Pack_32
155 si = *sx++;
156 ys = (si & 0xffff) * q + carry;
157 zs = (si >> 16) * q + (ys >> 16);
158 carry = zs >> 16;
159 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
160 borrow = (y & 0x10000) >> 16;
161 z = (*bx >> 16) - (zs & 0xffff) - borrow;
162 borrow = (z & 0x10000) >> 16;
163 Storeinc(bx, z, y);
164 #else
165 ys = *sx++ * q + carry;
166 carry = ys >> 16;
167 y = *bx - (ys & 0xffff) - borrow;
168 borrow = (y & 0x10000) >> 16;
169 *bx++ = y & 0xffff;
170 #endif
171 #endif
172 }
173 while(sx <= sxe);
174 if (!*bxe) {
175 bx = b->x;
176 while(--bxe > bx && !*bxe)
177 --n;
178 b->wds = n;
179 }
180 }
181 if (cmp(b, S) >= 0) {
182 q++;
183 borrow = 0;
184 carry = 0;
185 bx = b->x;
186 sx = S->x;
187 do {
188 #ifdef ULLong
189 ys = *sx++ + carry;
190 carry = ys >> 32;
191 y = *bx - (ys & 0xffffffffUL) - borrow;
192 borrow = y >> 32 & 1UL;
193 *bx++ = y & 0xffffffffUL;
194 #else
195 #ifdef Pack_32
196 si = *sx++;
197 ys = (si & 0xffff) + carry;
198 zs = (si >> 16) + (ys >> 16);
199 carry = zs >> 16;
200 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
201 borrow = (y & 0x10000) >> 16;
202 z = (*bx >> 16) - (zs & 0xffff) - borrow;
203 borrow = (z & 0x10000) >> 16;
204 Storeinc(bx, z, y);
205 #else
206 ys = *sx++ + carry;
207 carry = ys >> 16;
208 y = *bx - (ys & 0xffff) - borrow;
209 borrow = (y & 0x10000) >> 16;
210 *bx++ = y & 0xffff;
211 #endif
212 #endif
213 }
214 while(sx <= sxe);
215 bx = b->x;
216 bxe = bx + n;
217 if (!*bxe) {
218 while(--bxe > bx && !*bxe)
219 --n;
220 b->wds = n;
221 }
222 }
223 return q;
224 }
225