Open64: osprey/arith/ieee_fadd.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 /*
00042  *  Portable IEEE 754 floating-point addition and subtraction
00043  */
00044 
00045 #include "arith.internal.h"
00046 
00047 
00048 int
00049 ar_ifadd32 (AR_IEEE_32 *x,
00050       const AR_IEEE_32 *a,
00051       const AR_IEEE_32 *b,
00052       int roundmode)
00053 {
00054 
00055   int res = AR_STAT_OK;
00056   unsigned int lbits, y_lbits, rbits;
00057   unsigned int shift, diffsign, carry, x_expo;
00058   AR_IEEE_32 y;
00059 
00060   /* Ensure that the first argument has the largest exponent,
00061    * for simplicity.
00062    */
00063   if (b->expo > a->expo)
00064     y = *a, *x = *b;
00065   else
00066     y = *b, *x = *a;
00067   x_expo = x->expo;
00068 
00069   /* If either x or y is a NaN, it's the result. */
00070   if (IS_IEEE32_NaN(x)) {
00071     return res | AR_STAT_UNDEFINED;
00072   }
00073   if (IS_IEEE32_NaN(&y)) {
00074     *x = y;
00075     return res | AR_STAT_UNDEFINED;
00076   }
00077 
00078   /* Test for infinities */
00079   if (x_expo > AR_IEEE32_MAX_EXPO)
00080     if (y.expo > AR_IEEE32_MAX_EXPO)
00081       if (x->sign == y.sign) {
00082         /* infinities of same sign - return infinity */
00083         if (x->sign)
00084           res |= AR_STAT_NEGATIVE;
00085         return res | AR_STAT_OVERFLOW;
00086       } else {
00087         /* infinities of different signs - quiet NaN */
00088         QNaNIEEE32 (x);
00089         return res | AR_STAT_UNDEFINED;
00090       }
00091     else {
00092       /* infinity + finite quantity - return infinity */
00093       if (x->sign)
00094         res |= AR_STAT_NEGATIVE;
00095       return res | AR_STAT_OVERFLOW;
00096     }
00097   else if (y.expo > AR_IEEE32_MAX_EXPO) {
00098     *x = y;
00099     if (x->sign)
00100       res |= AR_STAT_NEGATIVE;
00101     return res | AR_STAT_OVERFLOW;
00102   }
00103 
00104   /* Test for denorms (they have zero exponents) to determine the
00105    * values of the implicit normalization bits.
00106    */
00107   if (ar_state_register.ar_denorms_trap &&
00108     ((!x_expo && IS_IEEE32_NZ_COEFF(x)) ||
00109      (!y.expo && IS_IEEE32_NZ_COEFF(&y)))) {
00110     /* operand is a denorm and denorms cause a trap */
00111     x->expo = AR_IEEE32_MAX_EXPO + 1;
00112     return res | AR_STAT_UNDEFINED;
00113   }
00114   lbits = !!x_expo;
00115   y_lbits = !!y.expo;
00116 
00117   /* Shift the coefficient of the second argument down (right). We
00118    * do this in parts so as not to overflow registers.
00119    */
00120   rbits = 0;
00121   shift = (x_expo - lbits) - (y.expo - y_lbits);
00122   if (shift > AR_IEEE32_COEFF_BITS + AR_IEEE32_ROUND_BITS) {
00123     rbits = !!(y.coeff0 | y.coeff1);
00124     y.sign = y.expo = y_lbits = 0;
00125     y.coeff0 = y.coeff1 = 0;
00126   } else
00127     for (; shift; shift--) {
00128       /* Sticky bit shifting */
00129       rbits = (rbits & 1) | (rbits >> 1) |
00130         ((y.coeff1 & 1) << (AR_IEEE32_ROUND_BITS-1));
00131       SHRIGHTIEEE32 (y);
00132       y.coeff0 |= y_lbits << (AR_IEEE32_C0_BITS - 1);
00133       y_lbits = 0;
00134     }
00135 
00136   /* If signs differ, complement the first argument; this is equivalent
00137    * to negating and then subtracting one, in a two's complement sense.
00138    */
00139   if (diffsign = (x->sign ^ y.sign)) {
00140     lbits ^= MASKR (2);
00141     NOTIEEE32 (*x);
00142   }
00143 
00144   /* Compute sum of coefficients */
00145   if (diffsign) {
00146     rbits += MASKR (AR_IEEE32_ROUND_BITS);
00147     carry = rbits >> AR_IEEE32_ROUND_BITS;
00148     rbits &= MASKR (AR_IEEE32_ROUND_BITS);
00149   } else
00150     carry = 0;
00151   ADDIEEE32 (*x, carry, *x, y);
00152   lbits = (carry + lbits + y_lbits) & MASKR (2);
00153 
00154   /* Opposite signs */
00155   if (diffsign)
00156     if (lbits & 2) {
00157       /* Complement */
00158       lbits ^= MASKR (2);
00159       NOTIEEE32 (*x);
00160       rbits ^= MASKR (AR_IEEE32_ROUND_BITS);
00161     } else {
00162       /* End-around carry */
00163       x->sign ^= 1;
00164       carry = 1 + rbits;
00165       rbits = carry & MASKR (AR_IEEE32_ROUND_BITS);
00166       carry >>= AR_IEEE32_ROUND_BITS;
00167       INCIEEE32 (*x, carry);
00168       lbits = (carry + lbits) & MASKR (2);
00169     }
00170 
00171   /* Shift rightward if carry, or if undenormalized */
00172   if (lbits & 2 || lbits && !x_expo) {
00173     rbits = (rbits >> 1) | rbits & 1 |
00174       ((x->coeff1 & 1) << (AR_IEEE32_ROUND_BITS - 1));
00175     SHRIGHTIEEE32 (*x);
00176     x->coeff0 |= lbits << (AR_IEEE32_C0_BITS - 1);
00177     lbits >>= 1;
00178     if (x_expo)
00179       x_expo++;
00180     else
00181       x_expo = 2;
00182   }
00183 
00184   /* Check for a zero result, as a special case */
00185   if (!(lbits | x->coeff0 | x->coeff1 | rbits)) {
00186     x->sign &= !diffsign; /* -0 + -0 == -0 */
00187     x->expo = 0;
00188     if (x->sign)
00189       res |= AR_STAT_NEGATIVE;
00190     return res | AR_STAT_ZERO;
00191   }
00192 
00193   return ar_i32norm (x_expo, lbits, rbits, x, roundmode);
00194 }
00195 
00196 
00197 int
00198 ar_ifsub32 (AR_IEEE_32 *x,
00199       const AR_IEEE_32 *a,
00200       const AR_IEEE_32 *b,
00201       int roundmode)
00202 {
00203 
00204   AR_IEEE_32 nb;
00205 
00206   /* If either a or b is a NaN, it's the result. */
00207   if (IS_IEEE32_NaN(a)) {
00208     *x = *a;
00209     return AR_STAT_UNDEFINED;
00210   }
00211   if (IS_IEEE32_NaN(b)) {
00212     *x = *b;
00213     return AR_STAT_UNDEFINED;
00214   }
00215 
00216   nb = *b;
00217   nb.sign ^= 1;
00218   return ar_ifadd32 (x, a, &nb, roundmode);
00219 }
00220 
00221 
00222 int
00223 ar_ifadd64 (AR_IEEE_64 *x,
00224       const AR_IEEE_64 *a,
00225       const AR_IEEE_64 *b,
00226       int roundmode)
00227 {
00228 
00229   int res = AR_STAT_OK;
00230   unsigned int lbits, y_lbits, rbits;
00231   unsigned int shift, diffsign, carry, x_expo;
00232   AR_IEEE_64 y;
00233 
00234   /* Ensure that the first argument has the largest exponent,
00235    * for simplicity.
00236    */
00237   if (b->expo > a->expo)
00238     y = *a, *x = *b;
00239   else
00240     y = *b, *x = *a;
00241   x_expo = x->expo;
00242 
00243   /* If either x or y is a NaN, it's the result. */
00244   if (IS_IEEE64_NaN(x)) {
00245     return res | AR_STAT_UNDEFINED;
00246   }
00247   if (IS_IEEE64_NaN(&y)) {
00248     *x = y;
00249     return res | AR_STAT_UNDEFINED;
00250   }
00251 
00252   /* Test for infinities */
00253   if (x_expo > AR_IEEE64_MAX_EXPO)
00254     if (y.expo > AR_IEEE64_MAX_EXPO)
00255       if (x->sign == y.sign) {
00256         /* infinities of same sign - return infinity */
00257         if (x->sign)
00258           res |= AR_STAT_NEGATIVE;
00259         return res | AR_STAT_OVERFLOW;
00260       } else {
00261         /* infinities of different signs - quiet NaN */
00262         QNaNIEEE64 (x);
00263         return res | AR_STAT_UNDEFINED;
00264       }
00265     else {
00266       /* infinity + finite quantity - return infinity */
00267       if (x->sign)
00268         res |= AR_STAT_NEGATIVE;
00269       return res | AR_STAT_OVERFLOW;
00270     }
00271   else if (y.expo > AR_IEEE64_MAX_EXPO) {
00272     *x = y;
00273     if (x->sign)
00274       res |= AR_STAT_NEGATIVE;
00275     return res | AR_STAT_OVERFLOW;
00276   }
00277 
00278   /* Test for denorms (they have zero exponents) to determine the
00279    * values of the implicit normalization bits.
00280    */
00281   if (ar_state_register.ar_denorms_trap &&
00282     ((!x_expo && IS_IEEE64_NZ_COEFF(x)) ||
00283      (!y.expo && IS_IEEE64_NZ_COEFF(&y)))) {
00284     /* operand is a denorm and denorms cause a trap */
00285     x->expo = AR_IEEE64_MAX_EXPO + 1;
00286     return res | AR_STAT_UNDEFINED;
00287   }
00288   lbits = !!x_expo;
00289   y_lbits = !!y.expo;
00290 
00291   /* Shift the coefficient of the second argument down (right). We
00292    * do this in parts so as not to overflow registers.
00293    */
00294   rbits = 0;
00295   shift = (x_expo - lbits) - (y.expo - y_lbits);
00296   if (shift > AR_IEEE64_COEFF_BITS + AR_IEEE64_ROUND_BITS) {
00297     rbits = !!(y.coeff0 | y.coeff1 | y.coeff2 | y.coeff3);
00298     y.sign = y.expo = y_lbits = 0;
00299     y.coeff0 = y.coeff1 = y.coeff2 = y.coeff3 = 0;
00300   } else
00301     for (; shift; shift--) {
00302       /* Sticky bit shifting */
00303       rbits = (rbits & 1) | (rbits >> 1) |
00304         ((y.coeff3 & 1) << (AR_IEEE64_ROUND_BITS-1));
00305       SHRIGHTIEEE64 (y);
00306       y.coeff0 |= y_lbits << (AR_IEEE64_C0_BITS - 1);
00307       y_lbits = 0;
00308     }
00309 
00310   /* If signs differ, complement the first argument; this is equivalent
00311    * to negating and then subtracting one, in a two's complement sense.
00312    */
00313   if (diffsign = (x->sign ^ y.sign)) {
00314     lbits ^= MASKR (2);
00315     NOTIEEE64 (*x);
00316   }
00317 
00318   /* Compute sum of coefficients */
00319   if (diffsign) {
00320     rbits += MASKR (AR_IEEE64_ROUND_BITS);
00321     carry = rbits >> AR_IEEE64_ROUND_BITS;
00322     rbits &= MASKR (AR_IEEE64_ROUND_BITS);
00323   } else
00324     carry = 0;
00325   ADDIEEE64 (*x, carry, *x, y);
00326   lbits = (carry + lbits + y_lbits) & MASKR (2);
00327 
00328   /* Opposite signs */
00329   if (diffsign)
00330     if (lbits & 2) {
00331       /* Complement */
00332       lbits ^= MASKR (2);
00333       NOTIEEE64 (*x);
00334       rbits ^= MASKR (AR_IEEE64_ROUND_BITS);
00335     } else {
00336       /* End-around carry */
00337       x->sign ^= 1;
00338       carry = 1 + rbits;
00339       rbits = carry & MASKR (AR_IEEE64_ROUND_BITS);
00340       carry >>= AR_IEEE64_ROUND_BITS;
00341       INCIEEE64 (*x, carry);
00342       lbits = (carry + lbits) & MASKR (2);
00343     }
00344 
00345   /* Shift rightward if carry, or if undenormalized */
00346   if (lbits & 2 || lbits && !x_expo) {
00347     rbits = (rbits >> 1) | rbits & 1 |
00348       ((x->coeff3 & 1) << (AR_IEEE64_ROUND_BITS - 1));
00349     SHRIGHTIEEE64 (*x);
00350     x->coeff0 |= lbits << (AR_IEEE64_C0_BITS - 1);
00351     lbits >>= 1;
00352     if (x_expo)
00353       x_expo++;
00354     else
00355       x_expo = 2;
00356   }
00357 
00358   /* Check for a zero result, as a special case */
00359   if (!(lbits | x->coeff0|x->coeff1|x->coeff2|x->coeff3 | rbits)) {
00360     x->sign &= !diffsign; /* -0 + -0 == -0 */
00361     x->expo = 0;
00362     if (x->sign)
00363       res |= AR_STAT_NEGATIVE;
00364     return res | AR_STAT_ZERO;
00365   }
00366 
00367   return ar_i64norm (x_expo, lbits, rbits, x, roundmode);
00368 }
00369 
00370 
00371 int
00372 ar_ifsub64 (AR_IEEE_64 *x,
00373       const AR_IEEE_64 *a,
00374       const AR_IEEE_64 *b,
00375       int roundmode)
00376 {
00377 
00378   AR_IEEE_64 nb;
00379 
00380   /* If either a or b is a NaN, it's the result. */
00381   if (IS_IEEE64_NaN(a)) {
00382     *x = *a;
00383     return AR_STAT_UNDEFINED;
00384   }
00385   if (IS_IEEE64_NaN(b)) {
00386     *x = *b;
00387     return AR_STAT_UNDEFINED;
00388   }
00389 
00390   nb = *b;
00391   nb.sign ^= 1;
00392   return ar_ifadd64 (x, a, &nb, roundmode);
00393 }
00394 
00395 #ifdef __mips
00396 /* Use native addition for MIPS */
00397 int
00398 ar_ifadd128(AR_IEEE_128 *x,
00399       const AR_IEEE_128 *a,
00400       const AR_IEEE_128 *b,
00401       int roundmode)
00402 {
00403    AR_TYPE ty = AR_Float_IEEE_NR_128;
00404    *(long double *)x = *(long double *)a + *(long double *)b;
00405    return AR_status((const AR_DATA *) x, &ty);
00406 }
00407 
00408 int
00409 ar_ifsub128( AR_IEEE_128 *x,
00410       const  AR_IEEE_128 *a,
00411       const  AR_IEEE_128 *b,
00412       int roundmode)
00413 {
00414    AR_TYPE ty = AR_Float_IEEE_NR_128;
00415    *(long double *)x = *(long double *)a - *(long double *)b;
00416    return AR_status((const AR_DATA *) x, &ty);
00417 }
00418 
00419 #else
00420 
00421 int
00422 ar_ifadd128(AR_IEEE_128 *x,
00423       const AR_IEEE_128 *a,
00424       const AR_IEEE_128 *b,
00425       int roundmode)
00426 {
00427 
00428   int res = AR_STAT_OK;
00429   unsigned int lbits, y_lbits, rbits;
00430   unsigned int shift, diffsign, carry, x_expo;
00431   AR_IEEE_128 y;
00432 
00433   /*
00434    * Use native arithmetic for MIPS.
00435    */
00436   if (HOST_IS_MIPS) {
00437     AR_TYPE ty = AR_Float_IEEE_NR_128;
00438 
00439     *(long double *)x = *(long double *)a + *(long double *)b;
00440     return AR_status((AR_DATA *) x, &ty);
00441   }
00442 
00443   /* Ensure that the first argument has the largest exponent,
00444    * for simplicity.
00445    */
00446   if (b->expo > a->expo)
00447     y = *a, *x = *b;
00448   else
00449     y = *b, *x = *a;
00450   x_expo = x->expo;
00451 
00452   /* If either x or y is a NaN, it's the result. */
00453   if (IS_IEEE128_NaN(x)) {
00454     return res | AR_STAT_UNDEFINED;
00455   }
00456   if (IS_IEEE128_NaN(&y)) {
00457     *x = y;
00458     return res | AR_STAT_UNDEFINED;
00459   }
00460 
00461   /* Test for infinities */
00462   if (x_expo > AR_IEEE128_MAX_EXPO)
00463     if (y.expo > AR_IEEE128_MAX_EXPO)
00464       if (x->sign == y.sign) {
00465         /* infinities of same sign - return infinity */
00466         if (x->sign)
00467           res |= AR_STAT_NEGATIVE;
00468         return res | AR_STAT_OVERFLOW;
00469       } else {
00470         /* infinities of different signs - quiet NaN */
00471         QNaNIEEE128 (x);
00472         return res | AR_STAT_UNDEFINED;
00473       }
00474     else {
00475       /* infinity + finite quantity - return infinity */
00476       if (x->sign)
00477         res |= AR_STAT_NEGATIVE;
00478       return res | AR_STAT_OVERFLOW;
00479     }
00480   else if (y.expo > AR_IEEE128_MAX_EXPO) {
00481     *x = y;
00482     if (x->sign)
00483       res |= AR_STAT_NEGATIVE;
00484     return res | AR_STAT_OVERFLOW;
00485   }
00486 
00487   /* Test for denorms (they have zero exponents) to determine the
00488    * values of the implicit normalization bits.
00489    */
00490   if (ar_state_register.ar_denorms_trap &&
00491     ((!x_expo && IS_IEEE128_NZ_COEFF(x)) ||
00492      (!y.expo && IS_IEEE128_NZ_COEFF(&y)))) {
00493     /* operand is a denorm and denorms cause a trap */
00494     x->expo = AR_IEEE128_MAX_EXPO + 1;
00495     return res | AR_STAT_UNDEFINED;
00496   }
00497   lbits = !!x_expo;
00498   y_lbits = !!y.expo;
00499 
00500   /* Shift the coefficient of the second argument down (right). We
00501    * do this in parts so as not to overflow registers.
00502    */
00503   rbits = 0;
00504   shift = (x_expo - lbits) - (y.expo - y_lbits);
00505   if (shift > AR_IEEE128_COEFF_BITS + AR_IEEE128_ROUND_BITS) {
00506     rbits = !!(y.coeff0 | y.coeff1 | y.coeff2 |
00507            y.coeff3 | y.coeff4 | y.coeff5 | y.coeff6);
00508     y.sign = y.expo = y_lbits = 0;
00509     y.coeff0 = y.coeff1 = y.coeff2 =
00510     y.coeff3 = y.coeff4 = y.coeff5 = y.coeff6 = 0;
00511   } else
00512     for (; shift; shift--) {
00513       /* Sticky bit shifting */
00514       rbits = (rbits & 1) | (rbits >> 1) |
00515         ((y.coeff6 & 1) << (AR_IEEE128_ROUND_BITS-1));
00516       SHRIGHTIEEE128 (y);
00517       y.coeff0 |= y_lbits << (AR_IEEE128_C0_BITS - 1);
00518       y_lbits = 0;
00519     }
00520 
00521   /* If signs differ, complement the first argument; this is equivalent
00522    * to negating and then subtracting one, in a two's complement sense.
00523    */
00524   if (diffsign = (x->sign ^ y.sign)) {
00525     lbits ^= MASKR (2);
00526     NOTIEEE128 (*x);
00527   }
00528 
00529   /* Compute sum of coefficients */
00530   if (diffsign) {
00531     rbits += MASKR (AR_IEEE128_ROUND_BITS);
00532     carry = rbits >> AR_IEEE128_ROUND_BITS;
00533     rbits &= MASKR (AR_IEEE128_ROUND_BITS);
00534   } else
00535     carry = 0;
00536   ADDIEEE128 (*x, carry, *x, y);
00537   lbits = (carry + lbits + y_lbits) & MASKR (2);
00538 
00539   /* Opposite signs */
00540   if (diffsign)
00541     if (lbits & 2) {
00542       /* Complement */
00543       lbits ^= MASKR (2);
00544       NOTIEEE128 (*x);
00545       rbits ^= MASKR (AR_IEEE128_ROUND_BITS);
00546     } else {
00547       /* End-around carry */
00548       x->sign ^= 1;
00549       carry = 1 + rbits;
00550       rbits = carry & MASKR (AR_IEEE128_ROUND_BITS);
00551       carry >>= AR_IEEE128_ROUND_BITS;
00552       INCIEEE128 (*x, carry);
00553       lbits = (carry + lbits) & MASKR (2);
00554     }
00555 
00556   /* Shift rightward if carry, or if undenormalized */
00557   if (lbits & 2 || lbits && !x_expo) {
00558     rbits = (rbits >> 1) | rbits & 1 |
00559       ((x->coeff6 & 1) << (AR_IEEE128_ROUND_BITS - 1));
00560     SHRIGHTIEEE128 (*x);
00561     x->coeff0 |= lbits << (AR_IEEE128_C0_BITS - 1);
00562     lbits >>= 1;
00563     if (x_expo)
00564       x_expo++;
00565     else
00566       x_expo = 2;
00567   }
00568 
00569   /* Check for a zero result, as a special case */
00570   if (!(lbits | rbits | x->coeff0 | x->coeff1 | x->coeff2 |
00571               x->coeff3 | x->coeff4 | x->coeff5 | x->coeff6)) {
00572     x->sign &= !diffsign; /* -0 + -0 == -0 */
00573     x->expo = 0;
00574     if (x->sign)
00575       res |= AR_STAT_NEGATIVE;
00576     return res | AR_STAT_ZERO;
00577   }
00578 
00579   return ar_i128norm (x_expo, lbits, rbits, x, roundmode);
00580 }
00581 
00582 
00583 int
00584 ar_ifsub128(AR_IEEE_128 *x,
00585       const AR_IEEE_128 *a,
00586       const AR_IEEE_128 *b,
00587       int roundmode)
00588 {
00589 
00590   AR_IEEE_128 nb;
00591 
00592   /* If either a or b is a NaN, it's the result. */
00593   if (IS_IEEE128_NaN(a)) {
00594     *x = *a;
00595     return AR_STAT_UNDEFINED;
00596   }
00597   if (IS_IEEE128_NaN(b)) {
00598     *x = *b;
00599     return AR_STAT_UNDEFINED;
00600   }
00601 
00602   nb = *b;
00603   nb.sign ^= 1;
00604   return ar_ifadd128 (x, a, &nb, roundmode);
00605 }
00606 #endif /* __mips */
00607 
00608 
00609 static char USMID [] = "\n%Z%%M%  %I% %G% %U%\n";
00610 static char rcsid [] = "$Id: ieee_fadd.c,v 1.1.1.1 2005/10/21 19:00:00 marcel Exp $";