Open64: osprey/arith/ieee_fmul.c Source File

00001 /*
00002  * Copyright 2004 PathScale, Inc.  All Rights Reserved.
00003  */
00004 
00005 /*
00006 
00007   Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
00008 
00009   This program is free software; you can redistribute it and/or modify it
00010   under the terms of version 2 of the GNU General Public License as
00011   published by the Free Software Foundation.
00012 
00013   This program is distributed in the hope that it would be useful, but
00014   WITHOUT ANY WARRANTY; without even the implied warranty of
00015   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
00016 
00017   Further, this software is distributed without any warranty that it is
00018   free of the rightful claim of any third person regarding infringement 
00019   or the like.  Any license provided herein, whether implied or 
00020   otherwise, applies only to this software file.  Patent licenses, if 
00021   any, provided herein do not apply to combinations of this program with 
00022   other software, or any other product whatsoever.  
00023 
00024   You should have received a copy of the GNU General Public License along
00025   with this program; if not, write the Free Software Foundation, Inc., 59
00026   Temple Place - Suite 330, Boston MA 02111-1307, USA.
00027 
00028   Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
00029   Mountain View, CA 94043, or:
00030 
00031   http://www.sgi.com
00032 
00033   For further information regarding this notice, see:
00034 
00035   http://oss.sgi.com/projects/GenInfo/NoticeExplan
00036 
00037 */
00038 
00039 
00040 /*
00041  *  Portable IEEE 754 floating-point multiplication evaluator
00042  */
00043 
00044 #include "arith.internal.h"
00045 
00046 
00047 
00048 int
00049 ar_ifmul32 (AR_IEEE_32 *x,
00050       const AR_IEEE_32 *a,
00051       const AR_IEEE_32 *b,
00052       int roundmode)
00053 {
00054 
00055   int i, res = AR_STAT_OK;
00056   unsigned long x_lbits, z_lbits, rbits, carry;
00057   signed int x_expo;
00058   AR_IEEE_32 x2, y, y2, z;
00059 
00060   /* If either a or b is a NaN, it's the result. */
00061   if (IS_IEEE32_NaN(a)) {
00062     *x = *a;
00063     return res | AR_STAT_UNDEFINED;
00064   }
00065   if (IS_IEEE32_NaN(b)) {
00066     *x = *b;
00067     return res | AR_STAT_UNDEFINED;
00068   }
00069 
00070   /* Test for infinities */
00071   if (b->expo > AR_IEEE32_MAX_EXPO)
00072     if (a->expo == 0 && !IS_IEEE32_NZ_COEFF(a)){
00073       /* zero * inf = quiet NaN */
00074       QNaNIEEE32 (x);
00075       return res | AR_STAT_UNDEFINED;
00076     } else {
00077       /* anything * infinity = infinity */
00078       if (i = a->sign ^ b->sign)
00079         res |= AR_STAT_NEGATIVE;
00080       *x = *b;
00081       x->sign = i;
00082       return res | AR_STAT_OVERFLOW;
00083     }
00084   if (a->expo > AR_IEEE32_MAX_EXPO)
00085     if (b->expo == 0 && !IS_IEEE32_NZ_COEFF(b)){
00086       /* infinity * zero = quiet NaN */
00087       QNaNIEEE32 (x);
00088       return res | AR_STAT_UNDEFINED;
00089     } else {
00090       /* infinity * anything = infinity */
00091       if (i = a->sign ^ b->sign)
00092         res |= AR_STAT_NEGATIVE;
00093       *x = *a;
00094       x->sign = i;
00095       return res | AR_STAT_OVERFLOW;
00096     }
00097 
00098   /* Test for denorms (they have zero exponents) to determine the
00099    * values of the implicit normalization bits; make them explicit.
00100    */
00101   if (ar_state_register.ar_denorms_trap &&
00102     ((!a->expo && IS_IEEE32_NZ_COEFF(a)) ||
00103      (!b->expo && IS_IEEE32_NZ_COEFF(b)))) {
00104     /* operand is a denorm and denorms cause a trap */
00105     x->expo = AR_IEEE32_MAX_EXPO + 1;
00106     return res | AR_STAT_UNDEFINED;
00107   }
00108   ZEROIEEE32 (y);
00109   y.coeff1 = !!a->expo;
00110   y2 = *a;
00111   z = *b;
00112   z_lbits = !!b->expo;
00113   x_expo = a->expo + b->expo + !a->expo + !b->expo - AR_IEEE32_EXPO_BIAS;
00114   if (x_expo <= 0)
00115     x_expo--;
00116   i = a->sign ^ b->sign;
00117 
00118   /* Sum the pyramid */
00119   if (z.coeff1 & 1) {
00120     x2 = *a;
00121     *x = y;
00122   }
00123   else {
00124     ZEROIEEE32 (*x);
00125     ZEROIEEE32 (x2);
00126   }
00127   x->sign = i;
00128   x_lbits = z_lbits & y.coeff1;
00129   SHLEFTIEEE32_2 (y, y2);
00130   SHRIGHTIEEE32 (z);
00131   z.coeff0 |= z_lbits << (AR_IEEE32_C0_BITS - 1);
00132   for (i = 1; i < AR_IEEE32_COEFF_BITS + 1; i++) {
00133 
00134     if (z.coeff1 & 1) {
00135       carry = 0;
00136       ADDIEEE32 (x2, carry, x2, y2);
00137       ADDIEEE32 (*x, carry, *x, y);
00138       x_lbits += carry;
00139     }
00140 
00141     SHLEFTIEEE32_2 (y, y2);
00142     SHRIGHTIEEE32 (z);
00143   }
00144 
00145   /* Extract rounding bits */
00146   rbits = x2.coeff0 >> (AR_IEEE32_C0_BITS - AR_IEEE32_ROUND_BITS);
00147   if (x2.coeff0 & MASKR (AR_IEEE32_C0_BITS - AR_IEEE32_ROUND_BITS) |
00148     x2.coeff1)
00149     rbits |= 1; /* sticky bit */
00150 
00151   /* Normalize and round */
00152   return ar_i32norm (x_expo, x_lbits, rbits, x, roundmode);
00153 }
00154 
00155 
00156 
00157 int
00158 ar_ifmul64 (AR_IEEE_64 *x,
00159       const AR_IEEE_64 *a,
00160       const AR_IEEE_64 *b,
00161       int roundmode)
00162 {
00163 
00164   int i, res = AR_STAT_OK;
00165   unsigned long x_lbits, z_lbits, rbits, carry;
00166   signed int x_expo;
00167   AR_IEEE_64 x2, y, y2, z;
00168 
00169   /* If either a or b is a NaN, it's the result. */
00170   if (IS_IEEE64_NaN(a)) {
00171     *x = *a;
00172     return res | AR_STAT_UNDEFINED;
00173   }
00174   if (IS_IEEE64_NaN(b)) {
00175     *x = *b;
00176     return res | AR_STAT_UNDEFINED;
00177   }
00178 
00179   /* Test for infinities */
00180   if (b->expo > AR_IEEE64_MAX_EXPO)
00181     if (a->expo == 0 && !IS_IEEE64_NZ_COEFF(a)) {
00182       /* zero * inf = quiet NaN */
00183       QNaNIEEE64 (x);
00184       return res | AR_STAT_UNDEFINED;
00185     } else {
00186       /* anything * infinity = infinity */
00187       if (i = a->sign ^ b->sign)
00188         res |= AR_STAT_NEGATIVE;
00189       *x = *b;
00190       x->sign = i;
00191       return res | AR_STAT_OVERFLOW;
00192     }
00193   if (a->expo > AR_IEEE64_MAX_EXPO)
00194     if (b->expo == 0 && !IS_IEEE64_NZ_COEFF(b)) {
00195       /* infinity * zero = quiet NaN */
00196       QNaNIEEE64 (x);
00197       return res | AR_STAT_UNDEFINED;
00198     } else {
00199       /* infinity * anything = infinity */
00200       if (i = a->sign ^ b->sign)
00201         res |= AR_STAT_NEGATIVE;
00202       *x = *a;
00203       x->sign = i;
00204       return res | AR_STAT_OVERFLOW;
00205     }
00206 
00207   /* Test for denorms (they have zero exponents) to determine the
00208    * values of the implicit normalization bits; make them explicit.
00209    */
00210   if (ar_state_register.ar_denorms_trap &&
00211     ((!a->expo && IS_IEEE64_NZ_COEFF(a)) ||
00212      (!b->expo && IS_IEEE64_NZ_COEFF(b)))) {
00213     /* operand is a denorm and denorms cause a trap */
00214     x->expo = AR_IEEE64_MAX_EXPO + 1;
00215     return res | AR_STAT_UNDEFINED;
00216   }
00217   ZEROIEEE64 (y);
00218   y.coeff3 = !!a->expo;
00219   y2 = *a;
00220   z = *b;
00221   z_lbits = !!b->expo;
00222   x_expo = a->expo + b->expo + !a->expo + !b->expo - AR_IEEE64_EXPO_BIAS;
00223   if (x_expo <= 0)
00224     x_expo--;
00225   i = a->sign ^ b->sign;
00226 
00227   /* Sum the pyramid */
00228   if (z.coeff3 & 1) {
00229     x2 = *a;
00230     *x = y;
00231   }
00232   else {
00233     ZEROIEEE64 (*x);
00234     ZEROIEEE64 (x2);
00235   }
00236   x->sign = i;
00237   x_lbits = z_lbits & y.coeff3;
00238   SHLEFTIEEE64_2 (y, y2);
00239   SHRIGHTIEEE64 (z);
00240   z.coeff0 |= z_lbits << (AR_IEEE64_C0_BITS - 1);
00241   for (i = 1; i < AR_IEEE64_COEFF_BITS + 1; i++) {
00242 
00243     if (z.coeff3 & 1) {
00244       carry = 0;
00245       ADDIEEE64 (x2, carry, x2, y2);
00246       ADDIEEE64 (*x, carry, *x, y);
00247       x_lbits += carry;
00248     }
00249 
00250     SHLEFTIEEE64_2 (y, y2);
00251     SHRIGHTIEEE64 (z);
00252   }
00253 
00254   /* Extract rounding bits */
00255   rbits = x2.coeff0 >> (AR_IEEE64_C0_BITS - AR_IEEE64_ROUND_BITS);
00256   if (x2.coeff0 & MASKR (AR_IEEE64_C0_BITS - AR_IEEE64_ROUND_BITS) |
00257     x2.coeff1 | x2.coeff2 | x2.coeff3)
00258     rbits |= 1; /* sticky bit */
00259 
00260   /* Normalize and round */
00261   return ar_i64norm (x_expo, x_lbits, rbits, x, roundmode);
00262 }
00263 
00264 #ifdef __mips
00265 
00266 int
00267 ar_ifmul128(AR_IEEE_128 *x,
00268       const AR_IEEE_128 *a,
00269       const AR_IEEE_128 *b,
00270       int roundmode)
00271 {
00272    AR_TYPE ty = AR_Float_IEEE_NR_128;
00273    *(long double *)x = *(long double *)a * *(long double *)b;
00274    return AR_status((const AR_DATA *) x, &ty);
00275 }
00276 
00277 #else
00278 
00279 int
00280 ar_ifmul128(AR_IEEE_128 *x,
00281       const AR_IEEE_128 *a,
00282       const AR_IEEE_128 *b,
00283       int roundmode)
00284 {
00285 
00286   int i, res = AR_STAT_OK;
00287   unsigned long x_lbits, z_lbits, rbits, carry;
00288   signed int x_expo;
00289   AR_IEEE_128 x2, y, y2, z;
00290 
00291   /*
00292    * Use native arithmetic for MIPS.
00293    */
00294   if (HOST_IS_MIPS) {
00295     AR_TYPE ty = AR_Float_IEEE_NR_128;
00296 
00297     *(long double *)x = *(long double *)a * *(long double *)b;
00298     return AR_status((AR_DATA *) x, &ty);
00299   }
00300 
00301   /* If either a or b is a NaN, it's the result. */
00302   if (IS_IEEE128_NaN(a)) {
00303     *x = *a;
00304     return res | AR_STAT_UNDEFINED;
00305   }
00306   if (IS_IEEE128_NaN(b)) {
00307     *x = *b;
00308     return res | AR_STAT_UNDEFINED;
00309   }
00310 
00311   /* Test for infinities */
00312   if (b->expo > AR_IEEE128_MAX_EXPO)
00313     if (a->expo == 0 && !IS_IEEE128_NZ_COEFF(a)) {
00314       /* zero * inf = quiet NaN */
00315       QNaNIEEE128 (x);
00316       return res | AR_STAT_UNDEFINED;
00317     } else {
00318       /* anything * infinity = infinity */
00319       if (i = a->sign ^ b->sign)
00320         res |= AR_STAT_NEGATIVE;
00321       *x = *b;
00322       x->sign = i;
00323       return res | AR_STAT_OVERFLOW;
00324     }
00325   if (a->expo > AR_IEEE128_MAX_EXPO)
00326     if (b->expo == 0 && !IS_IEEE128_NZ_COEFF(b)) {
00327       /* infinity * zero = quietNaN */
00328       QNaNIEEE128 (x);
00329       return res | AR_STAT_UNDEFINED;
00330     } else {
00331       /* infinity * anything = infinity */
00332       if (i = a->sign ^ b->sign)
00333         res |= AR_STAT_NEGATIVE;
00334       *x = *a;
00335       x->sign = i;
00336       return res | AR_STAT_OVERFLOW;
00337     }
00338 
00339   /* Test for denorms (they have zero exponents) to determine the
00340    * values of the implicit normalization bits; make them explicit.
00341    */
00342   if (ar_state_register.ar_denorms_trap &&
00343     ((!a->expo && IS_IEEE128_NZ_COEFF(a)) ||
00344      (!b->expo && IS_IEEE128_NZ_COEFF(b)))) {
00345     /* operand is a denorm and denorms cause a trap */
00346     x->expo = AR_IEEE128_MAX_EXPO + 1;
00347     return res | AR_STAT_UNDEFINED;
00348   }
00349   ZEROIEEE128 (y);
00350   y.coeff6 = !!a->expo;
00351   y2 = *a;
00352   z = *b;
00353   z_lbits = !!b->expo;
00354   x_expo = a->expo + b->expo + !a->expo + !b->expo - AR_IEEE128_EXPO_BIAS;
00355   if (x_expo <= 0)
00356     x_expo--;
00357   i = a->sign ^ b->sign;
00358 
00359   /* Sum the pyramid */
00360   if (z.coeff6 & 1) {
00361     x2 = *a;
00362     *x = y;
00363   }
00364   else {
00365     ZEROIEEE128 (*x);
00366     ZEROIEEE128 (x2);
00367   }
00368   x->sign = i;
00369   x_lbits = z_lbits & y.coeff6;
00370   SHLEFTIEEE128_2 (y, y2);
00371   SHRIGHTIEEE128 (z);
00372   z.coeff0 |= z_lbits << (AR_IEEE128_C0_BITS - 1);
00373   for (i = 1; i < AR_IEEE128_COEFF_BITS + 1; i++) {
00374 
00375     if (z.coeff6 & 1) {
00376       carry = 0;
00377       ADDIEEE128 (x2, carry, x2, y2);
00378       ADDIEEE128 (*x, carry, *x, y);
00379       x_lbits += carry;
00380     }
00381 
00382     SHLEFTIEEE128_2 (y, y2);
00383     SHRIGHTIEEE128 (z);
00384   }
00385 
00386   /* Extract rounding bits */
00387   rbits = x2.coeff0 >> (AR_IEEE128_C0_BITS - AR_IEEE128_ROUND_BITS);
00388   if (x2.coeff0 & MASKR (AR_IEEE128_C0_BITS - AR_IEEE128_ROUND_BITS) |
00389     x2.coeff1 | x2.coeff2 | x2.coeff3 | x2.coeff3 | x2.coeff5 | x2.coeff6)
00390     rbits |= 1; /* sticky bit */
00391 
00392   /* Normalize and round */
00393   return ar_i128norm (x_expo, x_lbits, rbits, x, roundmode);
00394 }
00395 #endif
00396 
00397 
00398 static char USMID [] = "\n%Z%%M%  %I% %G% %U%\n";
00399 static char rcsid [] = "$Id: ieee_fmul.c,v 1.1.1.1 2005/10/21 19:00:00 marcel Exp $";