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