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_sqrt.c
00040  *  $Revision$
00041  *  $Date$
00042  *  $Author$
00043  *  $Source$
00044  *
00045  * =======================================================================
00046  * =======================================================================
00047  */
00048 
00049 
00050 #ifndef _LP64
00051 #include <sgidefs.h>
00052 #endif /* _LP64 */
00053 #include <math.h>
00054 #include "defs.h"
00055 #include "quad.h"
00056 
00057 /* quad square root */
00058 
00059 /* note that this routine must be kept in sync with the corresponding libm
00060  * routine, qsqrt.
00061  */
00062 
00063 typedef union
00064 {
00065   struct
00066   {
00067     UINT32 hi;
00068     UINT32 lo;
00069   } word;
00070 
00071   double  d;
00072 } du;
00073 
00074 
00075 #ifndef __MATH_H__
00076 double  sqrt(double);
00077 
00078 #if _MIPS_ISA != _MIPS_ISA_MIPS1
00079 #pragma intrinsic (sqrt)
00080 #endif
00081 #endif /* __MATH_H__ */
00082 
00083 static  double  __dtwofloor(double);
00084 
00085 /* const1 is 1.0 + 2^(53 - 53/2), i.e. 1.0 + 2^27 */
00086 
00087 static  const du  const1 =
00088 {0x41a00000,  0x02000000};
00089 
00090 static const du   twopm53 =
00091 {0x3ca00000,  0x00000000};
00092 
00093 static const du   twopm54 =
00094 {0x3c900000,  0x00000000};
00095 
00096 static const du   twopm6 =
00097 {0x3f900000,  0x00000000};
00098 
00099 static const du   twop3 =
00100 {0x40200000,  0x00000000};
00101 
00102 static const du   twop108 =
00103 {0x46b00000,  0x00000000};
00104 
00105 #if defined(BUILD_OS_DARWIN)
00106 /* Can't use "pragma weak" to create aliases in Mach-O */
00107 QUAD c_q_sqrt(QUAD x, INT *p_err );
00108 QUAD __c_q_sqrt(QUAD x, INT *p_err ) { return c_q_sqrt(x, p_err); }
00109 #else /* defined(BUILD_OS_DARWIN) */
00110 #pragma weak c_q_sqrt = __c_q_sqrt
00111 #define c_q_sqrt __c_q_sqrt
00112 #endif /* defined(BUILD_OS_DARWIN) */
00113 
00114 QUAD
00115 c_q_sqrt(QUAD x, INT *p_err )
00116 {
00117 INT64 ix, xpt;
00118 INT32 n;
00119 double  xfactor;
00120 double  w;
00121 double  quarterulp;
00122 double  p;
00123 double  hc, tc;
00124 QUAD  c;
00125 QUAD  u, z;
00126 
00127   /* adapted from T. J. Dekker's sqrt2 subroutine */
00128 
00129   *p_err = 0;
00130 
00131   /* extract exponent of x for some quick screening */
00132 
00133   DBL2LL(x.hi, ix);
00134   xpt = (ix >> DMANTWIDTH);
00135   xpt &= 0x7ff;
00136 
00137   if ( ix < 0 )
00138   {
00139     z.hi = sqrt(x.hi);
00140     z.lo = 0.0;
00141 
00142     return ( z );
00143   }
00144 
00145   /* Avoid underflows and overflows in forming the products
00146      c.hi*const1.d, c.hi*c.hi, and hc*hc by scaling if
00147      necessary.  x should also be scaled if tc*tc is
00148      a denormal.
00149   */
00150 
00151   if ( (0x06d < xpt) && (xpt < 0x7fb) )
00152   {
00153     /* normal case */
00154 
00155     c.hi = sqrt(x.hi);
00156 
00157     /*u.ld = _prodl(c.hi, c.hi);*/
00158 
00159     p = c.hi*const1.d;
00160   
00161     hc = (c.hi - p) + p;
00162     tc = c.hi - hc;
00163   
00164     u.hi = c.hi*c.hi;
00165     u.lo = ((((hc*hc - u.hi) + hc*tc) + hc*tc) + tc*tc);
00166 
00167     c.lo = 0.5*(((x.hi - u.hi) - u.lo) + x.lo)/c.hi;
00168 
00169     /* if necessary, round c.lo so that the sum of c.hi and c.lo
00170        has at most 107 significant bits
00171     */
00172 
00173     /* Compute a quarter ulp of c.hi */
00174 
00175     DBL2LL(c.hi, ix);
00176     ix >>= DMANTWIDTH;
00177     ix <<= DMANTWIDTH;
00178     LL2DBL(ix, w); /* w = dtwofloor(c.hi) */
00179 
00180     /* Note that the size of an ulp changes at a
00181      * power of two.
00182      */
00183 
00184     if ( (c.hi == w) && (c.lo < 0.0) )
00185       w *= 0.5;
00186 
00187     /* round c.lo if it's less than 1/4 ulp of w */
00188 
00189     quarterulp = twopm54.d*fabs(w);
00190 
00191     if ( fabs(c.lo) < quarterulp )
00192     {
00193       if ( c.lo >= 0.0 )
00194       {
00195         c.lo = (quarterulp + c.lo) - quarterulp;
00196       }
00197       else
00198       {
00199         c.lo = quarterulp + (c.lo - quarterulp);
00200       }
00201 
00202     }
00203 
00204     z.hi = c.hi + c.lo;
00205     z.lo = (c.hi - z.hi) + c.lo;
00206 
00207     return ( z );
00208   }
00209 
00210   if ( xpt < 0x7ff )
00211   {
00212     if ( x.hi == 0.0 )
00213     {
00214       z.lo = 0.0;
00215       z.hi = x.hi;
00216 
00217       return ( z );
00218     }
00219 
00220     if ( xpt <= 0x06d )
00221     {
00222       x.hi *= twop108.d;
00223       x.lo *= twop108.d;
00224       xfactor = twopm54.d;
00225     }
00226     else
00227     {
00228       x.hi *= twopm6.d;
00229       x.lo *= twopm6.d;
00230       xfactor = twop3.d;
00231     }
00232 
00233     c.hi = sqrt(x.hi);
00234 
00235     /*u.ld = _prodl(c.hi, c.hi);*/
00236   
00237     p = c.hi*const1.d;
00238   
00239     hc = (c.hi - p) + p;
00240     tc = c.hi - hc;
00241   
00242     u.hi = c.hi*c.hi;
00243     u.lo = ((((hc*hc - u.hi) + hc*tc) + hc*tc) + tc*tc);
00244 
00245     c.lo = 0.5*(((x.hi - u.hi) - u.lo) + x.lo)/c.hi;
00246 
00247     /* if necessary, round c.lo so that the sum of c.hi and c.lo
00248        has at most 107 significant bits
00249     */
00250 
00251     /* Compute a quarter ulp of c.hi */
00252 
00253     DBL2LL(c.hi, ix);
00254     ix >>= DMANTWIDTH;
00255     ix <<= DMANTWIDTH;
00256     LL2DBL(ix, w); /* w = dtwofloor(c.hi) */
00257 
00258     /* Note that the size of an ulp changes at a
00259      * power of two.
00260      */
00261 
00262     if ( (c.hi == w) && (c.lo < 0.0) )
00263       w *= 0.5;
00264 
00265     /* round c.lo if it's less than 1/4 ulp of w */
00266 
00267     quarterulp = twopm54.d*fabs(w);
00268 
00269     if ( fabs(c.lo) < quarterulp )
00270     {
00271       if ( c.lo >= 0.0 )
00272       {
00273         c.lo = (quarterulp + c.lo) - quarterulp;
00274       }
00275       else
00276       {
00277         c.lo = quarterulp + (c.lo - quarterulp);
00278       }
00279 
00280     }
00281 
00282     c.hi *= xfactor;
00283     c.lo *= xfactor;
00284 
00285     z.hi = c.hi + c.lo;
00286     z.lo = (c.hi - z.hi) + c.lo;
00287 
00288     return ( z );
00289   }
00290 
00291   z.lo = 0.0;
00292   z.hi = sqrt(x.hi);
00293 
00294   return ( z );
00295 }
00296
osprey/common/util/c_q_sqrt.c
