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_add. 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_add(QUAD, QUAD, INT *); 00070 00071 #pragma weak c_q_add = __c_q_add 00072 #define c_q_add __c_q_add 00073 00074 double fabs(double); 00075 #pragma intrinsic (fabs) 00076 00077 QUAD 00078 c_q_add(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_add: xhi = %08x%08x\n", *(INT32 *)&xhi, *((INT32 *)&xhi + 1)); 00107 printf("c_q_add: xlo = %08x%08x\n", *(INT32 *)&xlo, *((INT32 *)&xlo + 1)); 00108 printf("c_q_add: yhi = %08x%08x\n", *(INT32 *)&yhi, *((INT32 *)&yhi + 1)); 00109 printf("c_q_add: ylo = %08x%08x\n", *(INT32 *)&ylo, *((INT32 *)&ylo + 1)); 00110 #endif 00111 00112 if ( xptxhi < xptyhi ) 00113 { 00114 tmp1 = xhi; 00115 xhi = yhi; 00116 yhi = tmp1; 00117 xptxhi = xptyhi; 00118 } 00119 00120 if ( fabs(xlo) < fabs(ylo) ) 00121 { 00122 tmp2 = xlo; 00123 xlo = ylo; 00124 ylo = tmp2; 00125 } 00126 00127 if ( xptxhi < 0x7fd ) 00128 { 00129 z.hi = xhi + yhi; 00130 z.lo = xhi - z.hi + yhi; 00131 00132 u = xlo + ylo; 00133 uu = xlo - u + ylo; 00134 00135 lo = z.lo + u; 00136 00137 w = z.hi + lo; 00138 ww = z.hi - w + lo; 00139 00140 rem = z.lo - lo + u; 00141 00142 ww += rem + uu; 00143 z.hi = w + ww; 00144 DBL2LL( z.hi, iz ); 00145 z.lo = w - z.hi + ww; 00146 00147 /* if necessary, round z.lo so that the sum of z.hi and z.lo has at most 00148 107 significant bits 00149 */ 00150 00151 /* first, compute a quarter ulp of z */ 00152 00153 iw = (iz >> DMANTWIDTH); 00154 iqulp = (iw & 0x7ff); 00155 iqulp -= 54; 00156 iqulp <<= DMANTWIDTH; 00157 00158 if ( iqulp > 0 ) 00159 { 00160 LL2DBL( iqulp, qulp ); 00161 iw <<= DMANTWIDTH; 00162 00163 /* Note that the size of an ulp changes at a 00164 * power of two. 00165 */ 00166 00167 if ( iw == iz ) 00168 goto fix; 00169 00170 if ( fabs(z.lo) >= qulp ) 00171 { 00172 qulp = 0.0; 00173 } 00174 else if ( z.lo < 0.0 ) 00175 qulp = -qulp; 00176 00177 z.lo += qulp; 00178 z.lo -= qulp; 00179 } 00180 00181 00182 #ifdef QUAD_DEBUG 00183 printf("q_add: z.hi = %08x%08x\n", *(INT32 *)&z.hi, *((INT32 *)&z.hi + 1)); 00184 printf("q_add: z.lo = %08x%08x\n", *(INT32 *)&z.lo, *((INT32 *)&z.lo + 1)); 00185 #endif 00186 00187 return ( z ); 00188 } 00189 else if ( xptxhi == 0x7ff ) 00190 { 00191 z.hi = xhi + yhi; 00192 z.lo = 0.0; 00193 00194 #ifdef QUAD_DEBUG 00195 printf("q_add: z.hi = %08x%08x\n", *(INT32 *)&z.hi, *((INT32 *)&z.hi + 1)); 00196 printf("q_add: z.lo = %08x%08x\n", *(INT32 *)&z.lo, *((INT32 *)&z.lo + 1)); 00197 #endif 00198 00199 return ( z ); 00200 } 00201 else 00202 { 00203 if ( fabs(yhi) < twop914.d ) 00204 { 00205 z.hi = xhi; 00206 z.lo = xlo; 00207 00208 #ifdef QUAD_DEBUG 00209 printf("q_add: z.hi = %08x%08x\n", *(INT32 *)&z.hi, *((INT32 *)&z.hi + 1)); 00210 printf("q_add: z.lo = %08x%08x\n", *(INT32 *)&z.lo, *((INT32 *)&z.lo + 1)); 00211 #endif 00212 00213 return ( z ); 00214 } 00215 00216 /* avoid overflow in intermediate computations by 00217 computing 4.0*(.25*x + .25*y) 00218 */ 00219 00220 xhi *= 0.25; 00221 xlo *= 0.25; 00222 yhi *= 0.25; 00223 ylo *= 0.25; 00224 00225 z.hi = xhi + yhi; 00226 z.lo = xhi - z.hi + yhi; 00227 00228 u = xlo + ylo; 00229 uu = xlo - u + ylo; 00230 00231 lo = z.lo + u; 00232 00233 w = z.hi + lo; 00234 ww = z.hi - w + lo; 00235 00236 rem = z.lo - lo + u; 00237 00238 ww += rem + uu; 00239 z.hi = w + ww; 00240 DBL2LL( z.hi, iz ); 00241 z.lo = w - z.hi + ww; 00242 00243 /* if necessary, round z.lo so that the sum of z.hi and z.lo has at most 00244 107 significant bits 00245 */ 00246 00247 /* first, compute a quarter ulp of z */ 00248 00249 iw = (iz >> DMANTWIDTH); 00250 iqulp = (iw & 0x7ff); 00251 iqulp -= 54; 00252 iqulp <<= DMANTWIDTH; 00253 00254 if ( iqulp > 0 ) 00255 { 00256 LL2DBL( iqulp, qulp ); 00257 iw <<= DMANTWIDTH; 00258 00259 /* Note that the size of an ulp changes at a 00260 * power of two. 00261 */ 00262 00263 if ( iw == iz ) 00264 goto fix2; 00265 00266 if ( fabs(z.lo) >= qulp ) 00267 { 00268 qulp = 0.0; 00269 } 00270 else if ( z.lo < 0.0 ) 00271 qulp = -qulp; 00272 00273 z.lo += qulp; 00274 z.lo -= qulp; 00275 } 00276 00277 z.hi *= 4.0; 00278 00279 if ( fabs(z.hi) == inf.d ) 00280 { 00281 z.lo = 0.0; 00282 return ( z ); 00283 } 00284 00285 z.lo *= 4.0; 00286 00287 return ( z ); 00288 00289 } 00290 00291 fix: 00292 if ( ((z.hi > 0.0) && (z.lo < 0.0)) || ((z.hi < 0.0) && (z.lo > 0.0)) ) 00293 qulp *= 0.5; 00294 00295 if ( fabs(z.lo) >= qulp ) 00296 { 00297 qulp = 0.0; 00298 } 00299 else if ( z.lo < 0.0 ) 00300 qulp = -qulp; 00301 00302 z.lo += qulp; 00303 z.lo -= qulp; 00304 00305 return ( z ); 00306 00307 fix2: 00308 if ( ((z.hi > 0.0) && (z.lo < 0.0)) || ((z.hi < 0.0) && (z.lo > 0.0)) ) 00309 qulp *= 0.5; 00310 00311 if ( fabs(z.lo) >= qulp ) 00312 { 00313 qulp = 0.0; 00314 } 00315 else if ( z.lo < 0.0 ) 00316 qulp = -qulp; 00317 00318 z.lo += qulp; 00319 z.lo -= qulp; 00320 00321 z.hi *= 4.0; 00322 00323 if ( fabs(z.hi) == inf.d ) 00324 { 00325 z.lo = 0.0; 00326 return ( z ); 00327 } 00328 00329 z.lo *= 4.0; 00330 00331 return ( z ); 00332 } 00333