Open64: osprey/common/util/c_q

00001 /*
00002 
00003   Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
00004 
00005   This program is free software; you can redistribute it and/or modify it
00006   under the terms of version 2 of the GNU General Public License as
00007   published by the Free Software Foundation.
00008 
00009   This program is distributed in the hope that it would be useful, but
00010   WITHOUT ANY WARRANTY; without even the implied warranty of
00011   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
00012 
00013   Further, this software is distributed without any warranty that it is
00014   free of the rightful claim of any third person regarding infringement 
00015   or the like.  Any license provided herein, whether implied or 
00016   otherwise, applies only to this software file.  Patent licenses, if 
00017   any, provided herein do not apply to combinations of this program with 
00018   other software, or any other product whatsoever.  
00019 
00020   You should have received a copy of the GNU General Public License along
00021   with this program; if not, write the Free Software Foundation, Inc., 59
00022   Temple Place - Suite 330, Boston MA 02111-1307, USA.
00023 
00024   Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
00025   Mountain View, CA 94043, or:
00026 
00027   http://www.sgi.com
00028 
00029   For further information regarding this notice, see:
00030 
00031   http://oss.sgi.com/projects/GenInfo/NoticeExplan
00032 
00033 */
00034 
00035 
00036 /* =======================================================================
00037  * =======================================================================
00038  *
00039  *  Module: c_q_div.c
00040  *  $Revision$
00041  *  $Date$
00042  *  $Author$
00043  *  $Source$
00044  *
00045  * =======================================================================
00046  * =======================================================================
00047  */
00048 
00049 
00050 #include "defs.h"
00051 #include "quad.h"
00052 
00053 /* note that this routine must be kept in sync with the corresponding libc
00054  * routine, q_div.
00055  */
00056 
00057 typedef union
00058 {
00059   struct
00060   {
00061     UINT32 hi;
00062     UINT32 lo;
00063   } word;
00064 
00065   double  d;
00066 } du;
00067 
00068 
00069 /* const1 is 1.0 + 2^(53 - 53/2), i.e. 1.0 + 2^27 */
00070 
00071 static  const du  const1 =
00072 {0x41a00000,  0x02000000};
00073 
00074 static const du   twopm968 =
00075 {0x03700000,  0x00000000};
00076 
00077 static const du   twopm54 =
00078 {0x3c900000,  0x00000000};
00079 
00080 static const du   twop52 =
00081 {0x43300000,  0x00000000};
00082 
00083 static const du   inf =
00084 {0x7ff00000,  0x00000000};
00085 
00086 extern  QUAD  c_q_div(QUAD, QUAD, INT *);
00087 
00088 #if defined(BUILD_OS_DARWIN)
00089 /* Can't use "pragma weak" to create aliases in Mach-O */
00090 QUAD c_q_div(QUAD x, QUAD y, INT *p_err );
00091 QUAD __c_q_div(QUAD x, QUAD y, INT *p_err ) { return c_q_div(x, y, p_err); }
00092 #else /* defined(BUILD_OS_DARWIN) */
00093 #pragma weak c_q_div = __c_q_div
00094 #define c_q_div __c_q_div
00095 #endif /* defined(BUILD_OS_DARWIN) */
00096 
00097 double  fabs(double);
00098 #pragma intrinsic (fabs)
00099 
00100 #define EXPBIAS 0x3ff
00101 
00102   /* computes the quotient of two long doubles */
00103 
00104 QUAD
00105 c_q_div(QUAD x, QUAD y, INT *p_err )
00106 {
00107 double  xhi, xlo, yhi, ylo;
00108 INT32 n;
00109 INT64 ixhi, iyhi;
00110 INT64 xptxhi, xptyhi;
00111 INT64 ix, iy, iz;
00112 double  c, cc, w, ww;
00113 double  quarterulp;
00114 double  xfactor, yfactor;
00115 double  p, hc, tc;
00116 double  hyhi, tyhi;
00117 double  z, zz;
00118 QUAD  u;
00119 QUAD  result;
00120 
00121   /* adapted from T. J. Dekker's div2 subroutine */
00122 
00123 #ifdef QUAD_DEBUG
00124   printf("c_q_div: xhi = %08x%08x\n", *(INT32 *)&xhi, *((INT32 *)&xhi + 1));
00125   printf("c_q_div: xlo = %08x%08x\n", *(INT32 *)&xlo, *((INT32 *)&xlo + 1));
00126   printf("c_q_div: yhi = %08x%08x\n", *(INT32 *)&yhi, *((INT32 *)&yhi + 1));
00127   printf("c_q_div: ylo = %08x%08x\n", *(INT32 *)&ylo, *((INT32 *)&ylo + 1));
00128 #endif
00129 
00130   *p_err = 0;
00131 
00132   xhi = x.hi; xlo = x.lo;
00133   yhi = y.hi; ylo = y.lo;
00134 
00135   /* extract exponents of y and x for some quick screening */
00136 
00137   DBL2LL(yhi, iyhi);
00138   xptyhi = (iyhi >> DMANTWIDTH);
00139   xptyhi &= 0x7ff;
00140 
00141   DBL2LL(xhi, ixhi);
00142   xptxhi = (ixhi >> DMANTWIDTH);
00143   xptxhi &= 0x7ff;
00144 
00145   /* Avoid underflows and overflows in forming the products
00146      x*const1.d, y*const1.d, x*y, and hx*hy by scaling if
00147      necessary.  x and y are also be scaled if tx*ty is
00148      a denormal.
00149   */
00150 
00151   if ( (0x30d < xptxhi) && (xptxhi < 0x4f1) && 
00152        (0x30d < xptyhi) && (xptyhi < 0x4f1) 
00153      )
00154   {
00155     /* normal case */
00156 
00157     c = xhi/yhi;
00158 /*
00159     u.ld = __qprod(c, yhi);
00160 */
00161     p = c*const1.d;
00162   
00163     hc = (c - p) + p;
00164     tc = c - hc;
00165   
00166     p = yhi*const1.d;
00167   
00168     hyhi = (yhi - p) + p;
00169     tyhi = yhi - hyhi;
00170   
00171     u.hi = c*yhi;
00172     u.lo = (((hc*hyhi - u.hi) + hc*tyhi) + hyhi*tc) + tc*tyhi;
00173 
00174     cc = ((xhi - u.hi - u.lo + xlo) - c*ylo)/yhi;
00175 
00176     z = c + cc;
00177     zz = (c - z) + cc;
00178 
00179 #ifdef QUAD_DEBUG
00180   printf("c_q_div: z = %08x%08x\n", *(INT32 *)&z, *((INT32 *)&z + 1));
00181   printf("c_q_div: zz = %08x%08x\n", *(INT32 *)&zz, *((INT32 *)&zz + 1));
00182 #endif
00183 
00184     /* if necessary, round zz so that the sum of z and zz has
00185        at most 107 significant bits
00186     */
00187 
00188     if ( fabs(z) >= twopm968.d )
00189     {
00190       /* determine true exponent of z + zz as a 107 bit number */
00191 
00192       DBL2LL(z, iz);
00193       iz >>= DMANTWIDTH;
00194       iz <<= DMANTWIDTH;
00195       LL2DBL(iz, w); /* w = dtwofloor(z) */
00196 
00197       if ( (z == w) && ((z > 0.0 && zz < 0.0) ||
00198            (z < 0.0 && zz > 0.0))
00199          )
00200         w *= 0.5;
00201 
00202       /* round zz if it's less than 1/4 ulp of w */
00203 
00204       quarterulp = twopm54.d*fabs(w);
00205 
00206       if ( fabs(zz) < quarterulp )
00207       {
00208         if ( zz >= 0.0 )
00209         {
00210           zz = (quarterulp + zz) - quarterulp;
00211         }
00212         else
00213         {
00214           zz = quarterulp + (zz - quarterulp);
00215         }
00216 
00217         w = z + zz;
00218         ww = (z - w) + zz;
00219 
00220         z = w;
00221         zz = ww;
00222       }
00223     }
00224 
00225     result.hi = z;
00226     result.lo = zz;
00227 
00228 #ifdef QUAD_DEBUG
00229   printf("c_q_div: result.hi = %08x%08x\n", *(INT32 *)&result.hi, *((INT32 *)&result.hi + 1));
00230   printf("c_q_div: result.lo = %08x%08x\n", *(INT32 *)&result.lo, *((INT32 *)&result.lo + 1));
00231 #endif
00232 
00233     return ( result );
00234   }
00235 
00236   if ( (xptxhi < 0x7ff) && (xptyhi < 0x7ff) )
00237   {
00238     if ( (xhi == 0.0) || (yhi == 0.0) )
00239     {
00240       result.hi = xhi/yhi;
00241       result.lo = 0.0;
00242 
00243 #ifdef QUAD_DEBUG
00244   printf("c_q_div: result.hi = %08x%08x\n", *(INT32 *)&result.hi, *((INT32 *)&result.hi + 1));
00245   printf("c_q_div: result.lo = %08x%08x\n", *(INT32 *)&result.lo, *((INT32 *)&result.lo + 1));
00246 #endif
00247 
00248       return ( result );
00249     }
00250 
00251     xfactor = 1.0;
00252     yfactor = 1.0;
00253 
00254     /* Scale x and y appropriately and use previous algorithm.
00255        Then unscale result.
00256     */
00257 
00258     if ( xptxhi <= 0x30d )
00259     {
00260       xhi *= twop52.d;  /* first, make sure x is normal */
00261       xlo *= twop52.d;
00262 
00263       /* now scale x so that its exponent is 0x30e */
00264 
00265       DBL2LL(xhi, ix);
00266       ix >>= DMANTWIDTH;
00267       n = (ix & 0x7ff);
00268 
00269       ix = (0x30e - n) + EXPBIAS;
00270       ix <<= DMANTWIDTH;
00271       LL2DBL(ix, c);
00272       xhi *= c;
00273       xlo *= c;
00274       ix = (n - 52 - 0x30e) + EXPBIAS;
00275       ix <<= DMANTWIDTH;
00276       LL2DBL(ix, xfactor);
00277     }
00278 
00279     if ( xptyhi <= 0x30d )
00280     {
00281       yhi *= twop52.d;  /* first, make sure y is normal */
00282       ylo *= twop52.d;
00283 
00284       /* now scale y so that its exponent is 0x30e */
00285 
00286       DBL2LL(yhi, iy);
00287       iy >>= DMANTWIDTH;
00288       n = (iy & 0x7ff);
00289 
00290       iy = (0x30e - n) + EXPBIAS;
00291       iy <<= DMANTWIDTH;
00292       LL2DBL(iy, c);
00293       yhi *= c;
00294       ylo *= c;
00295       iy = (0x30e - n + 52) + EXPBIAS;
00296       iy <<= DMANTWIDTH;
00297       LL2DBL(iy, yfactor);
00298     }
00299 
00300     if ( xptxhi >= 0x4f1 )
00301     {
00302       /* scale x so that its exponent is 0x4f1 */
00303 
00304       DBL2LL(xhi, ix);
00305       ix >>= DMANTWIDTH;
00306       n = (ix & 0x7ff);
00307 
00308       ix = (0x4f1 - n) + EXPBIAS;
00309       ix <<= DMANTWIDTH;
00310       LL2DBL(ix, c);
00311       xhi *= c;
00312       xlo *= c;
00313       ix = (n - 0x4f1) + EXPBIAS;
00314       ix <<= DMANTWIDTH;
00315       LL2DBL(ix, xfactor);
00316     }
00317 
00318     if ( xptyhi >= 0x4f1 )
00319     {
00320       /* scale y so that its exponent is 0x4f1 */
00321 
00322       DBL2LL(yhi, iy);
00323       iy >>= DMANTWIDTH;
00324       n = (iy & 0x7ff);
00325 
00326       iy = (0x4f1 - n) + EXPBIAS;
00327       iy <<= DMANTWIDTH;
00328       LL2DBL(iy, c);
00329       yhi *= c;
00330       ylo *= c;
00331       iy = (0x4f1 - n) + EXPBIAS;
00332       iy <<= DMANTWIDTH;
00333       LL2DBL(iy, yfactor);
00334     }
00335 
00336     c = xhi/yhi;
00337 
00338 /*
00339     u.ld = __qprod(c, yhi);
00340 */
00341     p = c*const1.d;
00342   
00343     hc = (c - p) + p;
00344     tc = c - hc;
00345   
00346     p = yhi*const1.d;
00347   
00348     hyhi = (yhi - p) + p;
00349     tyhi = yhi - hyhi;
00350   
00351     u.hi = c*yhi;
00352     u.lo = (((hc*hyhi - u.hi) + hc*tyhi) + hyhi*tc) + tc*tyhi;
00353 
00354     cc = ((xhi - u.hi - u.lo + xlo) - c*ylo)/yhi;
00355 
00356     z = c + cc;
00357     zz = (c - z) + cc;
00358 
00359 #ifdef QUAD_DEBUG
00360   printf("c_q_div: z = %08x%08x\n", *(INT32 *)&z, *((INT32 *)&z + 1));
00361   printf("c_q_div: zz = %08x%08x\n", *(INT32 *)&zz, *((INT32 *)&zz + 1));
00362 #endif
00363 
00364     /* if necessary, round zz so that the sum of z and zz has
00365        at most 107 significant bits
00366     */
00367 
00368     if ( fabs(z) >= twopm968.d )
00369     {
00370       /* determine true exponent of z + zz as a 107 bit number */
00371 
00372       DBL2LL(z, iz);
00373       iz >>= DMANTWIDTH;
00374       iz <<= DMANTWIDTH;
00375       LL2DBL(iz, w); /* w = dtwofloor(z) */
00376 
00377       if ( (z == w) && ((z > 0.0 && zz < 0.0) ||
00378            (z < 0.0 && zz > 0.0))
00379          )
00380         w *= 0.5;
00381 
00382       /* round zz if it's less than 1/4 ulp of w */
00383 
00384       quarterulp = twopm54.d*fabs(w);
00385 
00386       if ( fabs(zz) < quarterulp )
00387       {
00388         if ( zz >= 0.0 )
00389         {
00390           zz = (quarterulp + zz) - quarterulp;
00391         }
00392         else
00393         {
00394           zz = quarterulp + (zz - quarterulp);
00395         }
00396 
00397         w = z + zz;
00398         ww = (z - w) + zz;
00399 
00400         z = w;
00401         zz = ww;
00402       }
00403     }
00404 
00405     if ( ((xfactor <= 1.0) && (1.0 <= yfactor)) ||
00406          ((yfactor <= 1.0) && (1.0 <= xfactor))
00407        )
00408     {
00409       z = z*(xfactor*yfactor);
00410     }
00411     else
00412     {
00413       z *= xfactor;
00414       z *= yfactor;
00415     }
00416 
00417     if ( (z == 0.0) || (fabs(z) == inf.d) )
00418     {
00419       result.hi = z;
00420       result.lo = 0.0;
00421 
00422 #ifdef QUAD_DEBUG
00423   printf("c_q_div: result.hi = %08x%08x\n", *(INT32 *)&result.hi, *((INT32 *)&result.hi + 1));
00424   printf("c_q_div: result.lo = %08x%08x\n", *(INT32 *)&result.lo, *((INT32 *)&result.lo + 1));
00425 #endif
00426 
00427       return ( result );
00428     }
00429 
00430     if ( ((xfactor <= 1.0) && (1.0 <= yfactor)) ||
00431          ((yfactor <= 1.0) && (1.0 <= xfactor))
00432        )
00433     {
00434       zz = zz*(xfactor*yfactor);
00435     }
00436     else
00437     {
00438       zz *= xfactor;
00439       zz *= yfactor;
00440     }
00441 
00442     result.hi = z + zz;
00443     result.lo = (z - result.hi) + zz;
00444 
00445 #ifdef QUAD_DEBUG
00446   printf("c_q_div: result.hi = %08x%08x\n", *(INT32 *)&result.hi, *((INT32 *)&result.hi + 1));
00447   printf("c_q_div: result.lo = %08x%08x\n", *(INT32 *)&result.lo, *((INT32 *)&result.lo + 1));
00448 #endif
00449 
00450     return ( result );
00451   }
00452 
00453   result.hi = xhi/yhi;
00454   result.lo = 0.0;
00455 
00456 #ifdef QUAD_DEBUG
00457   printf("c_q_div: result.hi = %08x%08x\n", *(INT32 *)&result.hi, *((INT32 *)&result.hi + 1));
00458   printf("c_q_div: result.lo = %08x%08x\n", *(INT32 *)&result.lo, *((INT32 *)&result.lo + 1));
00459 #endif
00460 
00461   return ( result );
00462 }
00463
osprey/common/util/c_q_div.c
