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 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 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