| /* Fused multiply-add. | |
| Copyright (C) 2007, 2009, 2011-2026 Free Software Foundation, Inc. | |
| This file is free software: you can redistribute it and/or modify | |
| it under the terms of the GNU Lesser General Public License as | |
| published by the Free Software Foundation, either version 3 of the | |
| License, or (at your option) any later version. | |
| This file is distributed in the hope that it will be useful, | |
| but WITHOUT ANY WARRANTY; without even the implied warranty of | |
| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
| GNU Lesser General Public License for more details. | |
| You should have received a copy of the GNU Lesser General Public License | |
| along with this program. If not, see <https://www.gnu.org/licenses/>. */ | |
| /* Written by Bruno Haible <bruno@clisp.org>, 2011. */ | |
| #if ! defined USE_LONG_DOUBLE | |
| # include <config.h> | |
| #endif | |
| /* Specification. */ | |
| #include <math.h> | |
| #include <limits.h> | |
| #include <stdlib.h> | |
| #if HAVE_FEGETROUND | |
| # include <fenv.h> | |
| #endif | |
| #include "float+.h" | |
| #include "integer_length.h" | |
| #ifdef USE_LONG_DOUBLE | |
| # define FUNC fmal | |
| # define DOUBLE long double | |
| # define FREXP frexpl | |
| # define LDEXP ldexpl | |
| # define MIN_EXP LDBL_MIN_EXP | |
| # define MANT_BIT LDBL_MANT_BIT | |
| # define L_(literal) literal##L | |
| #elif ! defined USE_FLOAT | |
| # define FUNC fma | |
| # define DOUBLE double | |
| # define FREXP frexp | |
| # define LDEXP ldexp | |
| # define MIN_EXP DBL_MIN_EXP | |
| # define MANT_BIT DBL_MANT_BIT | |
| # define L_(literal) literal | |
| #else /* defined USE_FLOAT */ | |
| # define FUNC fmaf | |
| # define DOUBLE float | |
| # define FREXP frexpf | |
| # define LDEXP ldexpf | |
| # define MIN_EXP FLT_MIN_EXP | |
| # define MANT_BIT FLT_MANT_BIT | |
| # define L_(literal) literal##f | |
| #endif | |
| #undef MAX | |
| #define MAX(a,b) ((a) > (b) ? (a) : (b)) | |
| #undef MIN | |
| #define MIN(a,b) ((a) < (b) ? (a) : (b)) | |
| /* MSVC with option -fp:strict refuses to compile constant initializers that | |
| contain floating-point operations. Pacify this compiler. */ | |
| #if defined _MSC_VER && !defined __clang__ | |
| # pragma fenv_access (off) | |
| #endif | |
| /* Work around GCC 4.2.1 bug on OpenBSD 5.1/SPARC64. */ | |
| #if defined __GNUC__ && defined __sparc | |
| # define VOLATILE volatile | |
| #else | |
| # define VOLATILE | |
| #endif | |
| /* It is possible to write an implementation of fused multiply-add with | |
| floating-point operations alone. See | |
| Sylvie Boldo, Guillaume Melquiond: | |
| Emulation of FMA and correctly-rounded sums: proved algorithms using | |
| rounding to odd. | |
| <https://www.lri.fr/~melquion/doc/08-tc.pdf> | |
| But is it complicated. | |
| Here we take the simpler (and probably slower) approach of doing | |
| multi-precision arithmetic. */ | |
| /* We use the naming conventions of GNU gmp, but vastly simpler (and slower) | |
| algorithms. */ | |
| typedef unsigned int mp_limb_t; | |
| #define GMP_LIMB_BITS 32 | |
| static_assert (sizeof (mp_limb_t) * CHAR_BIT == GMP_LIMB_BITS); | |
| typedef unsigned long long mp_twolimb_t; | |
| #define GMP_TWOLIMB_BITS 64 | |
| static_assert (sizeof (mp_twolimb_t) * CHAR_BIT == GMP_TWOLIMB_BITS); | |
| /* Number of limbs needed for a single DOUBLE. */ | |
| #define NLIMBS1 ((MANT_BIT + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS) | |
| /* Number of limbs needed for the accumulator. */ | |
| #define NLIMBS3 (3 * NLIMBS1 + 1) | |
| /* Assuming 0.5 <= x < 1.0: | |
| Convert the mantissa (x * 2^DBL_MANT_BIT) to a sequence of limbs. */ | |
| static void | |
| decode (DOUBLE x, mp_limb_t limbs[NLIMBS1]) | |
| { | |
| /* I'm not sure whether it's safe to cast a 'double' value between | |
| 2^31 and 2^32 to 'unsigned int', therefore play safe and cast only | |
| 'double' values between 0 and 2^31 (to 'unsigned int' or 'int', | |
| doesn't matter). | |
| So, we split the MANT_BIT bits of x into a number of chunks of | |
| at most 31 bits. */ | |
| enum { chunk_count = (MANT_BIT - 1) / 31 + 1 }; | |
| /* Variables used for storing the bits limb after limb. */ | |
| mp_limb_t *p = limbs + NLIMBS1 - 1; | |
| mp_limb_t accu = 0; | |
| unsigned int bits_needed = MANT_BIT - (NLIMBS1 - 1) * GMP_LIMB_BITS; | |
| /* The bits bits_needed-1...0 need to be ORed into the accu. | |
| 1 <= bits_needed <= GMP_LIMB_BITS. */ | |
| /* Unroll the first 4 loop rounds. */ | |
| if (chunk_count > 0) | |
| { | |
| /* Here we still have MANT_BIT-0*31 bits to extract from x. */ | |
| enum { chunk_bits = MIN (31, MANT_BIT - 0 * 31) }; /* > 0, <= 31 */ | |
| x *= (mp_limb_t) 1 << chunk_bits; | |
| mp_limb_t d = (int) x; /* 0 <= d < 2^chunk_bits. */ | |
| x -= d; | |
| if (!(x >= L_(0.0) && x < L_(1.0))) | |
| abort (); | |
| if (bits_needed < chunk_bits) | |
| { | |
| /* store bits_needed bits */ | |
| accu |= d >> (chunk_bits - bits_needed); | |
| *p = accu; | |
| if (p == limbs) | |
| goto done; | |
| p--; | |
| /* hold (chunk_bits - bits_needed) bits */ | |
| accu = d << (GMP_LIMB_BITS - (chunk_bits - bits_needed)); | |
| bits_needed = GMP_LIMB_BITS - (chunk_bits - bits_needed); | |
| } | |
| else | |
| { | |
| /* store chunk_bits bits */ | |
| accu |= d << (bits_needed - chunk_bits); | |
| bits_needed -= chunk_bits; | |
| if (bits_needed == 0) | |
| { | |
| *p = accu; | |
| if (p == limbs) | |
| goto done; | |
| p--; | |
| accu = 0; | |
| bits_needed = GMP_LIMB_BITS; | |
| } | |
| } | |
| } | |
| if (chunk_count > 1) | |
| { | |
| /* Here we still have MANT_BIT-1*31 bits to extract from x. */ | |
| enum { chunk_bits = MIN (31, MAX (MANT_BIT - 1 * 31, 0)) }; /* > 0, <= 31 */ | |
| x *= (mp_limb_t) 1 << chunk_bits; | |
| mp_limb_t d = (int) x; /* 0 <= d < 2^chunk_bits. */ | |
| x -= d; | |
| if (!(x >= L_(0.0) && x < L_(1.0))) | |
| abort (); | |
| if (bits_needed < chunk_bits) | |
| { | |
| /* store bits_needed bits */ | |
| accu |= d >> (chunk_bits - bits_needed); | |
| *p = accu; | |
| if (p == limbs) | |
| goto done; | |
| p--; | |
| /* hold (chunk_bits - bits_needed) bits */ | |
| accu = d << (GMP_LIMB_BITS - (chunk_bits - bits_needed)); | |
| bits_needed = GMP_LIMB_BITS - (chunk_bits - bits_needed); | |
| } | |
| else | |
| { | |
| /* store chunk_bits bits */ | |
| accu |= d << (bits_needed - chunk_bits); | |
| bits_needed -= chunk_bits; | |
| if (bits_needed == 0) | |
| { | |
| *p = accu; | |
| if (p == limbs) | |
| goto done; | |
| p--; | |
| accu = 0; | |
| bits_needed = GMP_LIMB_BITS; | |
| } | |
| } | |
| } | |
| if (chunk_count > 2) | |
| { | |
| /* Here we still have MANT_BIT-2*31 bits to extract from x. */ | |
| enum { chunk_bits = MIN (31, MAX (MANT_BIT - 2 * 31, 0)) }; /* > 0, <= 31 */ | |
| x *= (mp_limb_t) 1 << chunk_bits; | |
| mp_limb_t d = (int) x; /* 0 <= d < 2^chunk_bits. */ | |
| x -= d; | |
| if (!(x >= L_(0.0) && x < L_(1.0))) | |
| abort (); | |
| if (bits_needed < chunk_bits) | |
| { | |
| /* store bits_needed bits */ | |
| accu |= d >> (chunk_bits - bits_needed); | |
| *p = accu; | |
| if (p == limbs) | |
| goto done; | |
| p--; | |
| /* hold (chunk_bits - bits_needed) bits */ | |
| accu = d << (GMP_LIMB_BITS - (chunk_bits - bits_needed)); | |
| bits_needed = GMP_LIMB_BITS - (chunk_bits - bits_needed); | |
| } | |
| else | |
| { | |
| /* store chunk_bits bits */ | |
| accu |= d << (bits_needed - chunk_bits); | |
| bits_needed -= chunk_bits; | |
| if (bits_needed == 0) | |
| { | |
| *p = accu; | |
| if (p == limbs) | |
| goto done; | |
| p--; | |
| accu = 0; | |
| bits_needed = GMP_LIMB_BITS; | |
| } | |
| } | |
| } | |
| if (chunk_count > 3) | |
| { | |
| /* Here we still have MANT_BIT-3*31 bits to extract from x. */ | |
| enum { chunk_bits = MIN (31, MAX (MANT_BIT - 3 * 31, 0)) }; /* > 0, <= 31 */ | |
| x *= (mp_limb_t) 1 << chunk_bits; | |
| mp_limb_t d = (int) x; /* 0 <= d < 2^chunk_bits. */ | |
| x -= d; | |
| if (!(x >= L_(0.0) && x < L_(1.0))) | |
| abort (); | |
| if (bits_needed < chunk_bits) | |
| { | |
| /* store bits_needed bits */ | |
| accu |= d >> (chunk_bits - bits_needed); | |
| *p = accu; | |
| if (p == limbs) | |
| goto done; | |
| p--; | |
| /* hold (chunk_bits - bits_needed) bits */ | |
| accu = d << (GMP_LIMB_BITS - (chunk_bits - bits_needed)); | |
| bits_needed = GMP_LIMB_BITS - (chunk_bits - bits_needed); | |
| } | |
| else | |
| { | |
| /* store chunk_bits bits */ | |
| accu |= d << (bits_needed - chunk_bits); | |
| bits_needed -= chunk_bits; | |
| if (bits_needed == 0) | |
| { | |
| *p = accu; | |
| if (p == limbs) | |
| goto done; | |
| p--; | |
| accu = 0; | |
| bits_needed = GMP_LIMB_BITS; | |
| } | |
| } | |
| } | |
| if (chunk_count > 4) | |
| { | |
| /* Here we still have MANT_BIT-4*31 bits to extract from x. */ | |
| /* Generic loop. */ | |
| for (size_t k = 4; k < chunk_count; k++) | |
| { | |
| size_t chunk_bits = MIN (31, MANT_BIT - k * 31); /* > 0, <= 31 */ | |
| x *= (mp_limb_t) 1 << chunk_bits; | |
| mp_limb_t d = (int) x; /* 0 <= d < 2^chunk_bits. */ | |
| x -= d; | |
| if (!(x >= L_(0.0) && x < L_(1.0))) | |
| abort (); | |
| if (bits_needed < chunk_bits) | |
| { | |
| /* store bits_needed bits */ | |
| accu |= d >> (chunk_bits - bits_needed); | |
| *p = accu; | |
| if (p == limbs) | |
| goto done; | |
| p--; | |
| /* hold (chunk_bits - bits_needed) bits */ | |
| accu = d << (GMP_LIMB_BITS - (chunk_bits - bits_needed)); | |
| bits_needed = GMP_LIMB_BITS - (chunk_bits - bits_needed); | |
| } | |
| else | |
| { | |
| /* store chunk_bits bits */ | |
| accu |= d << (bits_needed - chunk_bits); | |
| bits_needed -= chunk_bits; | |
| if (bits_needed == 0) | |
| { | |
| *p = accu; | |
| if (p == limbs) | |
| goto done; | |
| p--; | |
| accu = 0; | |
| bits_needed = GMP_LIMB_BITS; | |
| } | |
| } | |
| } | |
| } | |
| /* We shouldn't get here. */ | |
| abort (); | |
| done: ; | |
| #ifndef USE_LONG_DOUBLE /* On FreeBSD 6.1/x86, 'long double' numbers sometimes | |
| have excess precision. */ | |
| if (!(x == L_(0.0))) | |
| abort (); | |
| #endif | |
| } | |
| /* Multiply two sequences of limbs. */ | |
| static void | |
| multiply (mp_limb_t xlimbs[NLIMBS1], mp_limb_t ylimbs[NLIMBS1], | |
| mp_limb_t prod_limbs[2 * NLIMBS1]) | |
| { | |
| enum { len1 = NLIMBS1 }; | |
| enum { len2 = NLIMBS1 }; | |
| for (size_t k = len2; k > 0; ) | |
| prod_limbs[--k] = 0; | |
| for (size_t i = 0; i < len1; i++) | |
| { | |
| mp_limb_t digit1 = xlimbs[i]; | |
| mp_twolimb_t carry = 0; | |
| for (size_t j = 0; j < len2; j++) | |
| { | |
| mp_limb_t digit2 = ylimbs[j]; | |
| carry += (mp_twolimb_t) digit1 * (mp_twolimb_t) digit2; | |
| carry += prod_limbs[i + j]; | |
| prod_limbs[i + j] = (mp_limb_t) carry; | |
| carry = carry >> GMP_LIMB_BITS; | |
| } | |
| prod_limbs[i + len2] = (mp_limb_t) carry; | |
| } | |
| } | |
| DOUBLE | |
| FUNC (DOUBLE x, DOUBLE y, DOUBLE z) | |
| { | |
| if (isfinite (x) && isfinite (y)) | |
| { | |
| if (isfinite (z)) | |
| { | |
| /* x, y, z are all finite. */ | |
| if (x == L_(0.0) || y == L_(0.0)) | |
| return z; | |
| if (z == L_(0.0)) | |
| return x * y; | |
| /* x, y, z are all non-zero. | |
| The result is x * y + z. */ | |
| { | |
| int e; /* exponent of x * y + z */ | |
| int sign; | |
| mp_limb_t sum[NLIMBS3]; | |
| size_t sum_len; | |
| { | |
| int xys; /* sign of x * y */ | |
| int zs; /* sign of z */ | |
| int xye; /* sum of exponents of x and y */ | |
| int ze; /* exponent of z */ | |
| mp_limb_t summand1[NLIMBS3]; | |
| size_t summand1_len; | |
| mp_limb_t summand2[NLIMBS3]; | |
| size_t summand2_len; | |
| { | |
| mp_limb_t zlimbs[NLIMBS1]; | |
| mp_limb_t xylimbs[2 * NLIMBS1]; | |
| { | |
| DOUBLE xn; /* normalized part of x */ | |
| DOUBLE yn; /* normalized part of y */ | |
| DOUBLE zn; /* normalized part of z */ | |
| int xe; /* exponent of x */ | |
| int ye; /* exponent of y */ | |
| mp_limb_t xlimbs[NLIMBS1]; | |
| mp_limb_t ylimbs[NLIMBS1]; | |
| xys = 0; | |
| xn = x; | |
| if (x < 0) | |
| { | |
| xys = 1; | |
| xn = - x; | |
| } | |
| yn = y; | |
| if (y < 0) | |
| { | |
| xys = 1 - xys; | |
| yn = - y; | |
| } | |
| zs = 0; | |
| zn = z; | |
| if (z < 0) | |
| { | |
| zs = 1; | |
| zn = - z; | |
| } | |
| /* xn, yn, zn are all positive. | |
| The result is (-1)^xys * xn * yn + (-1)^zs * zn. */ | |
| xn = FREXP (xn, &xe); | |
| yn = FREXP (yn, &ye); | |
| zn = FREXP (zn, &ze); | |
| xye = xe + ye; | |
| /* xn, yn, zn are all < 1.0 and >= 0.5. | |
| The result is | |
| (-1)^xys * 2^xye * xn * yn + (-1)^zs * 2^ze * zn. */ | |
| if (xye < ze - MANT_BIT) | |
| { | |
| /* 2^xye * xn * yn < 2^xye <= 2^(ze-MANT_BIT-1) */ | |
| return z; | |
| } | |
| if (xye - 2 * MANT_BIT > ze) | |
| { | |
| /* 2^ze * zn < 2^ze <= 2^(xye-2*MANT_BIT-1). | |
| We cannot simply do | |
| return x * y; | |
| because it would round differently: A round-to-even | |
| in the multiplication can be a round-up or round-down | |
| here, due to z. So replace z with a value that doesn't | |
| require the use of long bignums but that rounds the | |
| same way. */ | |
| zn = L_(0.5); | |
| ze = xye - 2 * MANT_BIT - 1; | |
| } | |
| /* Convert mantissas of xn, yn, zn to limb sequences: | |
| xlimbs = 2^MANT_BIT * xn | |
| ylimbs = 2^MANT_BIT * yn | |
| zlimbs = 2^MANT_BIT * zn */ | |
| decode (xn, xlimbs); | |
| decode (yn, ylimbs); | |
| decode (zn, zlimbs); | |
| /* Multiply the mantissas of xn and yn: | |
| xylimbs = xlimbs * ylimbs */ | |
| multiply (xlimbs, ylimbs, xylimbs); | |
| } | |
| /* The result is | |
| (-1)^xys * 2^(xye-2*MANT_BIT) * xylimbs | |
| + (-1)^zs * 2^(ze-MANT_BIT) * zlimbs. | |
| Compute | |
| e = min (xye-2*MANT_BIT, ze-MANT_BIT) | |
| then | |
| summand1 = 2^(xye-2*MANT_BIT-e) * xylimbs | |
| summand2 = 2^(ze-MANT_BIT-e) * zlimbs */ | |
| e = MIN (xye - 2 * MANT_BIT, ze - MANT_BIT); | |
| if (e == xye - 2 * MANT_BIT) | |
| { | |
| /* Simply copy the limbs of xylimbs. */ | |
| for (size_t i = 0; i < 2 * NLIMBS1; i++) | |
| summand1[i] = xylimbs[i]; | |
| summand1_len = 2 * NLIMBS1; | |
| } | |
| else | |
| { | |
| size_t ediff = xye - 2 * MANT_BIT - e; | |
| /* Left shift the limbs of xylimbs by ediff bits. */ | |
| size_t ldiff = ediff / GMP_LIMB_BITS; | |
| size_t shift = ediff % GMP_LIMB_BITS; | |
| for (size_t i = 0; i < ldiff; i++) | |
| summand1[i] = 0; | |
| if (shift > 0) | |
| { | |
| mp_limb_t carry = 0; | |
| for (size_t i = 0; i < 2 * NLIMBS1; i++) | |
| { | |
| summand1[ldiff + i] = (xylimbs[i] << shift) | carry; | |
| carry = xylimbs[i] >> (GMP_LIMB_BITS - shift); | |
| } | |
| summand1[ldiff + 2 * NLIMBS1] = carry; | |
| summand1_len = ldiff + 2 * NLIMBS1 + 1; | |
| } | |
| else | |
| { | |
| for (size_t i = 0; i < 2 * NLIMBS1; i++) | |
| summand1[ldiff + i] = xylimbs[i]; | |
| summand1_len = ldiff + 2 * NLIMBS1; | |
| } | |
| /* Estimation of needed array size: | |
| ediff = (xye - 2 * MANT_BIT) - (ze - MANT_BIT) <= MANT_BIT + 1 | |
| therefore | |
| summand1_len | |
| = (ediff + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS + 2 * NLIMBS1 | |
| <= (MANT_BIT + GMP_LIMB_BITS) / GMP_LIMB_BITS + 2 * NLIMBS1 | |
| <= 3 * NLIMBS1 + 1 | |
| = NLIMBS3 */ | |
| if (!(summand1_len <= NLIMBS3)) | |
| abort (); | |
| } | |
| if (e == ze - MANT_BIT) | |
| { | |
| /* Simply copy the limbs of zlimbs. */ | |
| for (size_t i = 0; i < NLIMBS1; i++) | |
| summand2[i] = zlimbs[i]; | |
| summand2_len = NLIMBS1; | |
| } | |
| else | |
| { | |
| size_t ediff = ze - MANT_BIT - e; | |
| /* Left shift the limbs of zlimbs by ediff bits. */ | |
| size_t ldiff = ediff / GMP_LIMB_BITS; | |
| size_t shift = ediff % GMP_LIMB_BITS; | |
| for (size_t i = 0; i < ldiff; i++) | |
| summand2[i] = 0; | |
| if (shift > 0) | |
| { | |
| mp_limb_t carry = 0; | |
| for (size_t i = 0; i < NLIMBS1; i++) | |
| { | |
| summand2[ldiff + i] = (zlimbs[i] << shift) | carry; | |
| carry = zlimbs[i] >> (GMP_LIMB_BITS - shift); | |
| } | |
| summand2[ldiff + NLIMBS1] = carry; | |
| summand2_len = ldiff + NLIMBS1 + 1; | |
| } | |
| else | |
| { | |
| for (size_t i = 0; i < NLIMBS1; i++) | |
| summand2[ldiff + i] = zlimbs[i]; | |
| summand2_len = ldiff + NLIMBS1; | |
| } | |
| /* Estimation of needed array size: | |
| ediff = (ze - MANT_BIT) - (xye - 2 * MANT_BIT) <= 2 * MANT_BIT | |
| therefore | |
| summand2_len | |
| = (ediff + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS + NLIMBS1 | |
| <= (2 * MANT_BIT + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS + NLIMBS1 | |
| <= 3 * NLIMBS1 + 1 | |
| = NLIMBS3 */ | |
| if (!(summand2_len <= NLIMBS3)) | |
| abort (); | |
| } | |
| } | |
| /* The result is | |
| (-1)^xys * 2^e * summand1 + (-1)^zs * 2^e * summand2. */ | |
| if (xys == zs) | |
| { | |
| /* Perform an addition. */ | |
| size_t i; | |
| mp_limb_t carry; | |
| sign = xys; | |
| carry = 0; | |
| for (i = 0; i < MIN (summand1_len, summand2_len); i++) | |
| { | |
| mp_limb_t digit1 = summand1[i]; | |
| mp_limb_t digit2 = summand2[i]; | |
| sum[i] = digit1 + digit2 + carry; | |
| carry = (carry | |
| ? digit1 >= (mp_limb_t)-1 - digit2 | |
| : digit1 > (mp_limb_t)-1 - digit2); | |
| } | |
| if (summand1_len > summand2_len) | |
| for (; i < summand1_len; i++) | |
| { | |
| mp_limb_t digit1 = summand1[i]; | |
| sum[i] = carry + digit1; | |
| carry = carry && digit1 == (mp_limb_t)-1; | |
| } | |
| else | |
| for (; i < summand2_len; i++) | |
| { | |
| mp_limb_t digit2 = summand2[i]; | |
| sum[i] = carry + digit2; | |
| carry = carry && digit2 == (mp_limb_t)-1; | |
| } | |
| if (carry) | |
| sum[i++] = carry; | |
| sum_len = i; | |
| } | |
| else | |
| { | |
| /* Perform a subtraction. */ | |
| /* Compare summand1 and summand2 by magnitude. */ | |
| while (summand1[summand1_len - 1] == 0) | |
| summand1_len--; | |
| while (summand2[summand2_len - 1] == 0) | |
| summand2_len--; | |
| if (summand1_len > summand2_len) | |
| sign = xys; | |
| else if (summand1_len < summand2_len) | |
| sign = zs; | |
| else | |
| { | |
| size_t i = summand1_len; | |
| for (;;) | |
| { | |
| --i; | |
| if (summand1[i] > summand2[i]) | |
| { | |
| sign = xys; | |
| break; | |
| } | |
| if (summand1[i] < summand2[i]) | |
| { | |
| sign = zs; | |
| break; | |
| } | |
| if (i == 0) | |
| /* summand1 and summand2 are equal. */ | |
| return L_(0.0); | |
| } | |
| } | |
| if (sign == xys) | |
| { | |
| /* Compute summand1 - summand2. */ | |
| size_t i; | |
| mp_limb_t carry; | |
| carry = 0; | |
| for (i = 0; i < summand2_len; i++) | |
| { | |
| mp_limb_t digit1 = summand1[i]; | |
| mp_limb_t digit2 = summand2[i]; | |
| sum[i] = digit1 - digit2 - carry; | |
| carry = (carry ? digit1 <= digit2 : digit1 < digit2); | |
| } | |
| for (; i < summand1_len; i++) | |
| { | |
| mp_limb_t digit1 = summand1[i]; | |
| sum[i] = digit1 - carry; | |
| carry = carry && digit1 == 0; | |
| } | |
| if (carry) | |
| abort (); | |
| sum_len = summand1_len; | |
| } | |
| else | |
| { | |
| /* Compute summand2 - summand1. */ | |
| size_t i; | |
| mp_limb_t carry; | |
| carry = 0; | |
| for (i = 0; i < summand1_len; i++) | |
| { | |
| mp_limb_t digit1 = summand1[i]; | |
| mp_limb_t digit2 = summand2[i]; | |
| sum[i] = digit2 - digit1 - carry; | |
| carry = (carry ? digit2 <= digit1 : digit2 < digit1); | |
| } | |
| for (; i < summand2_len; i++) | |
| { | |
| mp_limb_t digit2 = summand2[i]; | |
| sum[i] = digit2 - carry; | |
| carry = carry && digit2 == 0; | |
| } | |
| if (carry) | |
| abort (); | |
| sum_len = summand2_len; | |
| } | |
| } | |
| } | |
| /* The result is | |
| (-1)^sign * 2^e * sum. */ | |
| /* Now perform the rounding to MANT_BIT mantissa bits. */ | |
| while (sum[sum_len - 1] == 0) | |
| sum_len--; | |
| /* Here we know that the most significant limb, sum[sum_len - 1], is | |
| non-zero. */ | |
| { | |
| /* How many bits the sum has. */ | |
| unsigned int sum_bits = | |
| integer_length (sum[sum_len - 1]) + (sum_len - 1) * GMP_LIMB_BITS; | |
| /* How many bits to keep when rounding. */ | |
| unsigned int keep_bits; | |
| /* How many bits to round off. */ | |
| unsigned int roundoff_bits; | |
| if (e + (int) sum_bits >= MIN_EXP) | |
| /* 2^e * sum >= 2^(MIN_EXP-1). | |
| result will be a normalized number. */ | |
| keep_bits = MANT_BIT; | |
| else if (e + (int) sum_bits >= MIN_EXP - MANT_BIT) | |
| /* 2^e * sum >= 2^(MIN_EXP-MANT_BIT-1). | |
| result will be a denormalized number or rounded to zero. */ | |
| keep_bits = e + (int) sum_bits - (MIN_EXP + MANT_BIT); | |
| else | |
| /* 2^e * sum < 2^(MIN_EXP-MANT_BIT-1). Round to zero. */ | |
| return L_(0.0); | |
| /* Note: 0 <= keep_bits <= MANT_BIT. */ | |
| if (sum_bits <= keep_bits) | |
| { | |
| /* Nothing to do. */ | |
| roundoff_bits = 0; | |
| keep_bits = sum_bits; | |
| } | |
| else | |
| { | |
| int round_up; | |
| roundoff_bits = sum_bits - keep_bits; /* > 0, <= sum_bits */ | |
| { | |
| #if HAVE_FEGETROUND && defined FE_TOWARDZERO | |
| /* Cf. <https://pubs.opengroup.org/onlinepubs/9699919799/functions/fegetround.html> */ | |
| int rounding_mode = fegetround (); | |
| if (rounding_mode == FE_TOWARDZERO) | |
| round_up = 0; | |
| # if defined FE_DOWNWARD /* not defined on sh4 */ | |
| else if (rounding_mode == FE_DOWNWARD) | |
| round_up = sign; | |
| # endif | |
| # if defined FE_UPWARD /* not defined on sh4 */ | |
| else if (rounding_mode == FE_UPWARD) | |
| round_up = !sign; | |
| # endif | |
| #else | |
| /* Cf. <https://pubs.opengroup.org/onlinepubs/9699919799/basedefs/float.h.html> */ | |
| int rounding_mode = FLT_ROUNDS; | |
| if (rounding_mode == 0) /* toward zero */ | |
| round_up = 0; | |
| else if (rounding_mode == 3) /* toward negative infinity */ | |
| round_up = sign; | |
| else if (rounding_mode == 2) /* toward positive infinity */ | |
| round_up = !sign; | |
| #endif | |
| else | |
| { | |
| /* Round to nearest. */ | |
| round_up = 0; | |
| /* Test bit (roundoff_bits-1). */ | |
| if ((sum[(roundoff_bits - 1) / GMP_LIMB_BITS] | |
| >> ((roundoff_bits - 1) % GMP_LIMB_BITS)) & 1) | |
| { | |
| /* Test bits roundoff_bits-1 .. 0. */ | |
| bool halfway = | |
| ((sum[(roundoff_bits - 1) / GMP_LIMB_BITS] | |
| & (((mp_limb_t) 1 << ((roundoff_bits - 1) % GMP_LIMB_BITS)) - 1)) | |
| == 0); | |
| if (halfway) | |
| for (int i = (roundoff_bits - 1) / GMP_LIMB_BITS - 1; i >= 0; i--) | |
| if (sum[i] != 0) | |
| { | |
| halfway = false; | |
| break; | |
| } | |
| if (halfway) | |
| /* Round to even. Test bit roundoff_bits. */ | |
| round_up = ((sum[roundoff_bits / GMP_LIMB_BITS] | |
| >> (roundoff_bits % GMP_LIMB_BITS)) & 1); | |
| else | |
| /* Round up. */ | |
| round_up = 1; | |
| } | |
| } | |
| } | |
| /* Perform the rounding. */ | |
| { | |
| size_t i = roundoff_bits / GMP_LIMB_BITS; | |
| { | |
| size_t j = i; | |
| while (j > 0) | |
| sum[--j] = 0; | |
| } | |
| if (round_up) | |
| { | |
| /* Round up. */ | |
| sum[i] = | |
| (sum[i] | |
| | (((mp_limb_t) 1 << (roundoff_bits % GMP_LIMB_BITS)) - 1)) | |
| + 1; | |
| if (sum[i] == 0) | |
| { | |
| /* Propagate carry. */ | |
| while (i < sum_len - 1) | |
| { | |
| ++i; | |
| sum[i] += 1; | |
| if (sum[i] != 0) | |
| break; | |
| } | |
| } | |
| /* sum[i] is the most significant limb that was | |
| incremented. */ | |
| if (i == sum_len - 1 && (sum[i] & (sum[i] - 1)) == 0) | |
| { | |
| /* Through the carry, one more bit is needed. */ | |
| if (sum[i] != 0) | |
| sum_bits += 1; | |
| else | |
| { | |
| /* Instead of requiring one more limb of memory, | |
| perform a shift by one bit, and adjust the | |
| exponent. */ | |
| sum[i] = (mp_limb_t) 1 << (GMP_LIMB_BITS - 1); | |
| e += 1; | |
| } | |
| /* The bit sequence has the form 1000...000. */ | |
| keep_bits = 1; | |
| } | |
| } | |
| else | |
| { | |
| /* Round down. */ | |
| sum[i] &= ((mp_limb_t) -1 << (roundoff_bits % GMP_LIMB_BITS)); | |
| if (i == sum_len - 1 && sum[i] == 0) | |
| /* The entire sum has become zero. */ | |
| return L_(0.0); | |
| } | |
| } | |
| } | |
| /* The result is | |
| (-1)^sign * 2^e * sum | |
| and here we know that | |
| 2^(sum_bits-1) <= sum < 2^sum_bits, | |
| and sum is a multiple of 2^(sum_bits-keep_bits), where | |
| 0 < keep_bits <= MANT_BIT and keep_bits <= sum_bits. | |
| (If keep_bits was initially 0, the rounding either returned zero | |
| or produced a bit sequence of the form 1000...000, setting | |
| keep_bits to 1.) */ | |
| { | |
| /* Split the keep_bits bits into chunks of at most 32 bits. */ | |
| unsigned int chunk_count = (keep_bits - 1) / GMP_LIMB_BITS + 1; | |
| /* 1 <= chunk_count <= ceil (sum_bits / GMP_LIMB_BITS) = sum_len. */ | |
| static const DOUBLE chunk_multiplier = /* 2^-GMP_LIMB_BITS */ | |
| L_(1.0) / ((DOUBLE) (1 << (GMP_LIMB_BITS / 2)) | |
| * (DOUBLE) (1 << ((GMP_LIMB_BITS + 1) / 2))); | |
| unsigned int shift = sum_bits % GMP_LIMB_BITS; | |
| DOUBLE fsum; | |
| if (MANT_BIT <= GMP_LIMB_BITS) | |
| { | |
| /* Since keep_bits <= MANT_BIT <= GMP_LIMB_BITS, | |
| chunk_count is 1. No need for a loop. */ | |
| if (shift == 0) | |
| fsum = (DOUBLE) sum[sum_len - 1]; | |
| else | |
| fsum = (DOUBLE) | |
| ((sum[sum_len - 1] << (GMP_LIMB_BITS - shift)) | |
| | (sum_len >= 2 ? sum[sum_len - 2] >> shift : 0)); | |
| } | |
| else | |
| { | |
| int k; | |
| k = chunk_count - 1; | |
| if (shift == 0) | |
| { | |
| /* First loop round. */ | |
| fsum = (DOUBLE) sum[sum_len - k - 1]; | |
| /* General loop. */ | |
| while (--k >= 0) | |
| { | |
| fsum *= chunk_multiplier; | |
| fsum += (DOUBLE) sum[sum_len - k - 1]; | |
| } | |
| } | |
| else | |
| { | |
| /* First loop round. */ | |
| { | |
| VOLATILE mp_limb_t chunk = | |
| (sum[sum_len - k - 1] << (GMP_LIMB_BITS - shift)) | |
| | (sum_len >= k + 2 ? sum[sum_len - k - 2] >> shift : 0); | |
| fsum = (DOUBLE) chunk; | |
| } | |
| /* General loop. */ | |
| while (--k >= 0) | |
| { | |
| fsum *= chunk_multiplier; | |
| { | |
| VOLATILE mp_limb_t chunk = | |
| (sum[sum_len - k - 1] << (GMP_LIMB_BITS - shift)) | |
| | (sum[sum_len - k - 2] >> shift); | |
| fsum += (DOUBLE) chunk; | |
| } | |
| } | |
| } | |
| } | |
| fsum = LDEXP (fsum, e + (int) sum_bits - GMP_LIMB_BITS); | |
| return (sign ? - fsum : fsum); | |
| } | |
| } | |
| } | |
| } | |
| else | |
| return z; | |
| } | |
| else | |
| return (x * y) + z; | |
| } |