1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License, Version 1.0 only
6 * (the "License"). You may not use this file except in compliance
7 * with the License.
8 *
9 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
10 * or http://www.opensolaris.org/os/licensing.
11 * See the License for the specific language governing permissions
12 * and limitations under the License.
13 *
14 * When distributing Covered Code, include this CDDL HEADER in each
15 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
16 * If applicable, add the following below this CDDL HEADER, with the
17 * fields enclosed by brackets "[]" replaced with your own identifying
18 * information: Portions Copyright [yyyy] [name of copyright owner]
19 *
20 * CDDL HEADER END
21 */
22 #pragma ident "%Z%%M% %I% %E% SMI"
23
24 /*
25 * Copyright (c) 1988 by Sun Microsystems, Inc.
26 */
27
28 #include "_Qquad.h"
29 #include "_Qglobals.h"
30
31 void
_fp_div(px,py,pz)32 _fp_div(px, py, pz)
33 unpacked *px, *py, *pz;
34
35 {
36 unsigned r[4],*y,q,c;
37 int n;
38
39 *pz = *px;
40 pz->sign = px->sign ^ py->sign;
41
42 if ((py->fpclass == fp_quiet) || (py->fpclass == fp_signaling)) {
43 *pz = *py;
44 return;
45 }
46 switch (px->fpclass) {
47 case fp_quiet:
48 case fp_signaling:
49 return;
50 case fp_zero:
51 case fp_infinity:
52 if (px->fpclass == py->fpclass) { /* 0/0 or inf/inf */
53 fpu_error_nan(pz);
54 pz->fpclass = fp_quiet;
55 }
56 return;
57 case fp_normal:
58 switch (py->fpclass) {
59 case fp_zero: /* number/0 */
60 fpu_set_exception(fp_division);
61 pz->fpclass = fp_infinity;
62 return;
63 case fp_infinity: /* number/inf */
64 pz->fpclass = fp_zero;
65 return;
66 }
67 }
68
69 /* Now x and y are both normal or subnormal. */
70
71 r[0] = px->significand[0];
72 r[1] = px->significand[1];
73 r[2] = px->significand[2];
74 r[3] = px->significand[3];
75 y = py->significand;
76
77 if(fpu_cmpli(r,y,4)>=0)
78 pz->exponent = px->exponent - py->exponent;
79 else
80 pz->exponent = px->exponent - py->exponent - 1;
81
82 q=0;
83 while(q<0x10000) { /* generate quo[0] */
84 q<<=1;
85 if(fpu_cmpli(r,y,4)>=0) {
86 q += 1; /* if r>y do r-=y and q+=1 */
87 c = 0;
88 c = fpu_sub3wc(&r[3],r[3],y[3],c);
89 c = fpu_sub3wc(&r[2],r[2],y[2],c);
90 c = fpu_sub3wc(&r[1],r[1],y[1],c);
91 c = fpu_sub3wc(&r[0],r[0],y[0],c);
92 }
93 r[0] = (r[0]<<1)|((r[1]&0x80000000)>>31); /* r << 1 */
94 r[1] = (r[1]<<1)|((r[2]&0x80000000)>>31);
95 r[2] = (r[2]<<1)|((r[3]&0x80000000)>>31);
96 r[3] = (r[3]<<1);
97 }
98 pz->significand[0]=q;
99 q=0; /* generate quo[1] */
100 n = 32;
101 while(n--) {
102 q<<=1;
103 if(fpu_cmpli(r,y,4)>=0) {
104 q += 1; /* if r>y do r-=y and q+=1 */
105 c = 0;
106 c = fpu_sub3wc(&r[3],r[3],y[3],c);
107 c = fpu_sub3wc(&r[2],r[2],y[2],c);
108 c = fpu_sub3wc(&r[1],r[1],y[1],c);
109 c = fpu_sub3wc(&r[0],r[0],y[0],c);
110 }
111 r[0] = (r[0]<<1)|((r[1]&0x80000000)>>31); /* r << 1 */
112 r[1] = (r[1]<<1)|((r[2]&0x80000000)>>31);
113 r[2] = (r[2]<<1)|((r[3]&0x80000000)>>31);
114 r[3] = (r[3]<<1);
115 }
116 pz->significand[1] = q;
117 q=0; /* generate quo[2] */
118 n = 32;
119 while(n--) {
120 q<<=1;
121 if(fpu_cmpli(r,y,4)>=0) {
122 q += 1; /* if r>y do r-=y and q+=1 */
123 c = 0;
124 c = fpu_sub3wc(&r[3],r[3],y[3],c);
125 c = fpu_sub3wc(&r[2],r[2],y[2],c);
126 c = fpu_sub3wc(&r[1],r[1],y[1],c);
127 c = fpu_sub3wc(&r[0],r[0],y[0],c);
128 }
129 r[0] = (r[0]<<1)|((r[1]&0x80000000)>>31); /* r << 1 */
130 r[1] = (r[1]<<1)|((r[2]&0x80000000)>>31);
131 r[2] = (r[2]<<1)|((r[3]&0x80000000)>>31);
132 r[3] = (r[3]<<1);
133 }
134 pz->significand[2] = q;
135 q=0; /* generate quo[3] */
136 n = 32;
137 while(n--) {
138 q<<=1;
139 if(fpu_cmpli(r,y,4)>=0) {
140 q += 1; /* if r>y do r-=y and q+=1 */
141 c = 0;
142 c = fpu_sub3wc(&r[3],r[3],y[3],c);
143 c = fpu_sub3wc(&r[2],r[2],y[2],c);
144 c = fpu_sub3wc(&r[1],r[1],y[1],c);
145 c = fpu_sub3wc(&r[0],r[0],y[0],c);
146 }
147 r[0] = (r[0]<<1)|((r[1]&0x80000000)>>31); /* r << 1 */
148 r[1] = (r[1]<<1)|((r[2]&0x80000000)>>31);
149 r[2] = (r[2]<<1)|((r[3]&0x80000000)>>31);
150 r[3] = (r[3]<<1);
151 }
152 pz->significand[3] = q;
153 if((r[0]|r[1]|r[2]|r[3])==0) pz->sticky = pz->rounded = 0;
154 else {
155 pz->sticky = 1; /* half way case won't occur */
156 if(fpu_cmpli(r,y,4)>=0) pz->rounded = 1;
157 }
158 }
159
160 void
_fp_sqrt(px,pz)161 _fp_sqrt(px, pz)
162 unpacked *px, *pz;
163
164 { /* *pz gets sqrt(*px) */
165
166 unsigned *x,r,c,q,t[4],s[4];
167 *pz = *px;
168 switch (px->fpclass) {
169 case fp_quiet:
170 case fp_signaling:
171 case fp_zero:
172 return;
173 case fp_infinity:
174 if (px->sign == 1) { /* sqrt(-inf) */
175 fpu_error_nan(pz);
176 pz->fpclass = fp_quiet;
177 }
178 return;
179 case fp_normal:
180 if (px->sign == 1) { /* sqrt(-norm) */
181 fpu_error_nan(pz);
182 pz->fpclass = fp_quiet;
183 return;
184 }
185 }
186
187 /* Now x is normal. */
188 x = px->significand;
189 if (px->exponent & 1) { /* sqrt(1.f * 2**odd) = sqrt (2.+2f) *
190 * 2**(odd-1)/2 */
191 pz->exponent = (px->exponent - 1) / 2;
192 x[0] = (x[0]<<1)|((x[1]&0x80000000)>>31); /* x<<1 */
193 x[1] = (x[1]<<1)|((x[2]&0x80000000)>>31);
194 x[2] = (x[2]<<1)|((x[3]&0x80000000)>>31);
195 x[3] = (x[3]<<1);
196 } else { /* sqrt(1.f * 2**even) = sqrt (1.f) *
197 * 2**(even)/2 */
198 pz->exponent = px->exponent / 2;
199 }
200 s[0]=s[1]=s[2]=s[3]=t[0]=t[1]=t[2]=t[3]=0;
201 q = 0;
202 r = 0x00010000;
203 while(r!=0) { /* compute sqrt[0] */
204 t[0] = s[0]+r;
205 if(t[0]<=x[0]) {
206 s[0] = t[0]+r;
207 x[0] -= t[0];
208 q += r;
209 }
210 x[0] = (x[0]<<1)|((x[1]&0x80000000)>>31); /* x<<1 */
211 x[1] = (x[1]<<1)|((x[2]&0x80000000)>>31);
212 x[2] = (x[2]<<1)|((x[3]&0x80000000)>>31);
213 x[3] = (x[3]<<1);
214 r>>=1;
215 }
216 pz->significand[0] = q;
217 q = 0;
218 r = 0x80000000;
219 while(r!=0) { /* compute sqrt[1] */
220 t[1] = s[1]+r; /* no carry */
221 t[0] = s[0];
222 if(fpu_cmpli(t,x,2)<=0) {
223 c = 0;
224 c = fpu_add3wc(&s[1],t[1],r,c);
225 c = fpu_add3wc(&s[0],t[0],0,c);
226 c = 0;
227 c = fpu_sub3wc(&x[1],x[1],t[1],c);
228 c = fpu_sub3wc(&x[0],x[0],t[0],c);
229 q += r;
230 }
231 x[0] = (x[0]<<1)|((x[1]&0x80000000)>>31); /* x<<1 */
232 x[1] = (x[1]<<1)|((x[2]&0x80000000)>>31);
233 x[2] = (x[2]<<1)|((x[3]&0x80000000)>>31);
234 x[3] = (x[3]<<1);
235 r>>=1;
236 }
237 pz->significand[1] = q;
238 q = 0;
239 r = 0x80000000;
240 while(r!=0) { /* compute sqrt[2] */
241 t[2] = s[2]+r; /* no carry */
242 t[1] = s[1];
243 t[0] = s[0];
244 if(fpu_cmpli(t,x,3)<=0) {
245 c = 0;
246 c = fpu_add3wc(&s[2],t[2],r,c);
247 c = fpu_add3wc(&s[1],t[1],0,c);
248 c = fpu_add3wc(&s[0],t[0],0,c);
249 c = 0;
250 c = fpu_sub3wc(&x[2],x[2],t[2],c);
251 c = fpu_sub3wc(&x[1],x[1],t[1],c);
252 c = fpu_sub3wc(&x[0],x[0],t[0],c);
253 q += r;
254 }
255 x[0] = (x[0]<<1)|((x[1]&0x80000000)>>31); /* x<<1 */
256 x[1] = (x[1]<<1)|((x[2]&0x80000000)>>31);
257 x[2] = (x[2]<<1)|((x[3]&0x80000000)>>31);
258 x[3] = (x[3]<<1);
259 r>>=1;
260 }
261 pz->significand[2] = q;
262 q = 0;
263 r = 0x80000000;
264 while(r!=0) { /* compute sqrt[3] */
265 t[3] = s[3]+r; /* no carry */
266 t[2] = s[2];
267 t[1] = s[1];
268 t[0] = s[0];
269 if(fpu_cmpli(t,x,4)<=0) {
270 c = 0;
271 c = fpu_add3wc(&s[3],t[3],r,c);
272 c = fpu_add3wc(&s[2],t[2],0,c);
273 c = fpu_add3wc(&s[1],t[1],0,c);
274 c = fpu_add3wc(&s[0],t[0],0,c);
275 c = 0;
276 c = fpu_sub3wc(&x[3],x[3],t[3],c);
277 c = fpu_sub3wc(&x[2],x[2],t[2],c);
278 c = fpu_sub3wc(&x[1],x[1],t[1],c);
279 c = fpu_sub3wc(&x[0],x[0],t[0],c);
280 q += r;
281 }
282 x[0] = (x[0]<<1)|((x[1]&0x80000000)>>31); /* x<<1 */
283 x[1] = (x[1]<<1)|((x[2]&0x80000000)>>31);
284 x[2] = (x[2]<<1)|((x[3]&0x80000000)>>31);
285 x[3] = (x[3]<<1);
286 r>>=1;
287 }
288 pz->significand[3] = q;
289 if((x[0]|x[1]|x[2]|x[3])==0) {
290 pz->sticky = pz->rounded = 0;
291 } else {
292 pz->sticky = 1;
293 if(fpu_cmpli(s,x,4)<0) pz->rounded=1; else pz->rounded = 0;
294 }
295 }
296