usr/src/lib/libbc/libc/gen/common/_Qfdiv.c

   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 }