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