Open64 (mfef90, whirl2f, and IR tools)
TAG: version-openad; SVN changeset: 916
|
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 * 00040 * ======================================================================= 00041 * ======================================================================= 00042 */ 00043 00044 00045 #include "defs.h" 00046 #include "quad.h" 00047 00048 /* note that this routine must be kept in sync with the corresponding libc 00049 * routine, q_sub. 00050 */ 00051 00052 typedef union 00053 { 00054 struct 00055 { 00056 UINT32 hi; 00057 UINT32 lo; 00058 } word; 00059 00060 double d; 00061 } du; 00062 00063 static const du twop914 = 00064 {0x79100000, 0x00000000}; 00065 00066 static const du inf = 00067 {0x7ff00000, 0x00000000}; 00068 00069 extern QUAD c_q_sub(QUAD, QUAD, INT *); 00070 00071 #pragma weak c_q_sub = __c_q_sub 00072 #define c_q_sub __c_q_sub 00073 00074 double fabs(double); 00075 #pragma intrinsic (fabs) 00076 00077 QUAD 00078 c_q_sub(QUAD x, QUAD y, INT *p_err ) 00079 { 00080 double xhi, xlo, yhi, ylo; 00081 INT32 ixhi, iyhi; 00082 INT32 xptxhi, xptyhi; 00083 INT64 iz, iw, iqulp; 00084 double w, ww; 00085 double u, uu; 00086 double qulp; 00087 QUAD z; 00088 double tmp1, tmp2, lo, rem; 00089 00090 /* adapted from T. J. Dekker's add2 subroutine */ 00091 00092 *p_err = 0; 00093 00094 xhi = x.hi; xlo = x.lo; 00095 yhi = y.hi; ylo = y.lo; 00096 00097 iyhi = *(INT32 *)&yhi; 00098 xptyhi = (iyhi >> 20); 00099 xptyhi &= 0x7ff; 00100 00101 ixhi = *(INT32 *)&xhi; 00102 xptxhi = (ixhi >> 20); 00103 xptxhi &= 0x7ff; 00104 00105 #ifdef QUAD_DEBUG 00106 printf("c_q_sub: xhi = %08x%08x\n", *(INT32 *)&xhi, *((INT32 *)&xhi + 1)); 00107 printf("c_q_sub: xlo = %08x%08x\n", *(INT32 *)&xlo, *((INT32 *)&xlo + 1)); 00108 printf("c_q_sub: yhi = %08x%08x\n", *(INT32 *)&yhi, *((INT32 *)&yhi + 1)); 00109 printf("c_q_sub: ylo = %08x%08x\n", *(INT32 *)&ylo, *((INT32 *)&ylo + 1)); 00110 #endif 00111 00112 yhi = -yhi; 00113 ylo = -ylo; 00114 00115 if ( xptxhi < xptyhi ) 00116 { 00117 tmp1 = xhi; 00118 xhi = yhi; 00119 yhi = tmp1; 00120 xptxhi = xptyhi; 00121 } 00122 00123 if ( fabs(xlo) < fabs(ylo) ) 00124 { 00125 tmp2 = xlo; 00126 xlo = ylo; 00127 ylo = tmp2; 00128 } 00129 00130 if ( xptxhi < 0x7fd ) 00131 { 00132 z.hi = xhi + yhi; 00133 z.lo = xhi - z.hi + yhi; 00134 00135 u = xlo + ylo; 00136 uu = xlo - u + ylo; 00137 00138 lo = z.lo + u; 00139 00140 w = z.hi + lo; 00141 ww = z.hi - w + lo; 00142 00143 rem = z.lo - lo + u; 00144 00145 ww += rem + uu; 00146 z.hi = w + ww; 00147 DBL2LL( z.hi, iz ); 00148 z.lo = w - z.hi + ww; 00149 00150 /* if necessary, round z.lo so that the sum of z.hi and z.lo has at most 00151 107 significant bits 00152 */ 00153 00154 /* first, compute a quarter ulp of z */ 00155 00156 iw = (iz >> DMANTWIDTH); 00157 iqulp = (iw & 0x7ff); 00158 iqulp -= 54; 00159 iqulp <<= DMANTWIDTH; 00160 00161 if ( iqulp > 0 ) 00162 { 00163 LL2DBL( iqulp, qulp ); 00164 iw <<= DMANTWIDTH; 00165 00166 /* Note that the size of an ulp changes at a 00167 * power of two. 00168 */ 00169 00170 if ( iw == iz ) 00171 goto fix; 00172 00173 if ( fabs(z.lo) >= qulp ) 00174 { 00175 qulp = 0.0; 00176 } 00177 else if ( z.lo < 0.0 ) 00178 qulp = -qulp; 00179 00180 z.lo += qulp; 00181 z.lo -= qulp; 00182 } 00183 00184 00185 #ifdef QUAD_DEBUG 00186 printf("q_add: z.hi = %08x%08x\n", *(INT32 *)&z.hi, *((INT32 *)&z.hi + 1)); 00187 printf("q_add: z.lo = %08x%08x\n", *(INT32 *)&z.lo, *((INT32 *)&z.lo + 1)); 00188 #endif 00189 00190 return ( z ); 00191 } 00192 else if ( xptxhi == 0x7ff ) 00193 { 00194 z.hi = xhi + yhi; 00195 z.lo = 0.0; 00196 00197 #ifdef QUAD_DEBUG 00198 printf("q_add: z.hi = %08x%08x\n", *(INT32 *)&z.hi, *((INT32 *)&z.hi + 1)); 00199 printf("q_add: z.lo = %08x%08x\n", *(INT32 *)&z.lo, *((INT32 *)&z.lo + 1)); 00200 #endif 00201 00202 return ( z ); 00203 } 00204 else 00205 { 00206 if ( fabs(yhi) < twop914.d ) 00207 { 00208 z.hi = xhi; 00209 z.lo = xlo; 00210 00211 #ifdef QUAD_DEBUG 00212 printf("q_add: z.hi = %08x%08x\n", *(INT32 *)&z.hi, *((INT32 *)&z.hi + 1)); 00213 printf("q_add: z.lo = %08x%08x\n", *(INT32 *)&z.lo, *((INT32 *)&z.lo + 1)); 00214 #endif 00215 00216 return ( z ); 00217 } 00218 00219 /* avoid overflow in intermediate computations by 00220 computing 4.0*(.25*x + .25*y) 00221 */ 00222 00223 xhi *= 0.25; 00224 xlo *= 0.25; 00225 yhi *= 0.25; 00226 ylo *= 0.25; 00227 00228 z.hi = xhi + yhi; 00229 z.lo = xhi - z.hi + yhi; 00230 00231 u = xlo + ylo; 00232 uu = xlo - u + ylo; 00233 00234 lo = z.lo + u; 00235 00236 w = z.hi + lo; 00237 ww = z.hi - w + lo; 00238 00239 rem = z.lo - lo + u; 00240 00241 ww += rem + uu; 00242 z.hi = w + ww; 00243 DBL2LL( z.hi, iz ); 00244 z.lo = w - z.hi + ww; 00245 00246 /* if necessary, round z.lo so that the sum of z.hi and z.lo has at most 00247 107 significant bits 00248 */ 00249 00250 /* first, compute a quarter ulp of z */ 00251 00252 iw = (iz >> DMANTWIDTH); 00253 iqulp = (iw & 0x7ff); 00254 iqulp -= 54; 00255 iqulp <<= DMANTWIDTH; 00256 00257 if ( iqulp > 0 ) 00258 { 00259 LL2DBL( iqulp, qulp ); 00260 iw <<= DMANTWIDTH; 00261 00262 /* Note that the size of an ulp changes at a 00263 * power of two. 00264 */ 00265 00266 if ( iw == iz ) 00267 goto fix2; 00268 00269 if ( fabs(z.lo) >= qulp ) 00270 { 00271 qulp = 0.0; 00272 } 00273 else if ( z.lo < 0.0 ) 00274 qulp = -qulp; 00275 00276 z.lo += qulp; 00277 z.lo -= qulp; 00278 } 00279 00280 z.hi *= 4.0; 00281 00282 if ( fabs(z.hi) == inf.d ) 00283 { 00284 z.lo = 0.0; 00285 return ( z ); 00286 } 00287 00288 z.lo *= 4.0; 00289 00290 return ( z ); 00291 00292 } 00293 00294 fix: 00295 if ( ((z.hi > 0.0) && (z.lo < 0.0)) || ((z.hi < 0.0) && (z.lo > 0.0)) ) 00296 qulp *= 0.5; 00297 00298 if ( fabs(z.lo) >= qulp ) 00299 { 00300 qulp = 0.0; 00301 } 00302 else if ( z.lo < 0.0 ) 00303 qulp = -qulp; 00304 00305 z.lo += qulp; 00306 z.lo -= qulp; 00307 00308 return ( z ); 00309 00310 fix2: 00311 if ( ((z.hi > 0.0) && (z.lo < 0.0)) || ((z.hi < 0.0) && (z.lo > 0.0)) ) 00312 qulp *= 0.5; 00313 00314 if ( fabs(z.lo) >= qulp ) 00315 { 00316 qulp = 0.0; 00317 } 00318 else if ( z.lo < 0.0 ) 00319 qulp = -qulp; 00320 00321 z.lo += qulp; 00322 z.lo -= qulp; 00323 00324 z.hi *= 4.0; 00325 00326 if ( fabs(z.hi) == inf.d ) 00327 { 00328 z.lo = 0.0; 00329 return ( z ); 00330 } 00331 00332 z.lo *= 4.0; 00333 00334 return ( z ); 00335 } 00336