Open64: osprey/arith/ieee_fdiv.c Source File

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