| /* Remainder. |
| Copyright (C) 2011-2020 Free Software Foundation, Inc. |
| |
| This program is free software: you can redistribute it and/or modify |
| it under the terms of the GNU General Public License as published by |
| the Free Software Foundation; either version 3 of the License, or |
| (at your option) any later version. |
| |
| This program 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 General Public License for more details. |
| |
| You should have received a copy of the GNU General Public License |
| along with this program. If not, see <https://www.gnu.org/licenses/>. */ |
| |
| #if ! defined USE_LONG_DOUBLE |
| # include <config.h> |
| #endif |
| |
| /* Specification. */ |
| #include <math.h> |
| |
| #include <float.h> |
| #include <stdlib.h> |
| |
| #ifdef USE_LONG_DOUBLE |
| # define FMOD fmodl |
| # define DOUBLE long double |
| # define MANT_DIG LDBL_MANT_DIG |
| # define L_(literal) literal##L |
| # define FREXP frexpl |
| # define LDEXP ldexpl |
| # define FABS fabsl |
| # define TRUNC truncl |
| # define ISNAN isnanl |
| #else |
| # define FMOD fmod |
| # define DOUBLE double |
| # define MANT_DIG DBL_MANT_DIG |
| # define L_(literal) literal |
| # define FREXP frexp |
| # define LDEXP ldexp |
| # define FABS fabs |
| # define TRUNC trunc |
| # define ISNAN isnand |
| #endif |
| |
| /* MSVC with option -fp:strict refuses to compile constant initializers that |
| contain floating-point operations. Pacify this compiler. */ |
| #ifdef _MSC_VER |
| # pragma fenv_access (off) |
| #endif |
| |
| #undef NAN |
| #if defined _MSC_VER |
| static DOUBLE zero; |
| # define NAN (zero / zero) |
| #else |
| # define NAN (L_(0.0) / L_(0.0)) |
| #endif |
| |
| /* To avoid rounding errors during the computation of x - q * y, |
| there are three possibilities: |
| - Use fma (- q, y, x). |
| - Have q be a single bit at a time, and compute x - q * y |
| through a subtraction. |
| - Have q be at most MANT_DIG/2 bits at a time, and compute |
| x - q * y by splitting y into two halves: |
| y = y1 * 2^(MANT_DIG/2) + y0 |
| x - q * y = (x - q * y1 * 2^MANT_DIG/2) - q * y0. |
| The latter is probably the most efficient. */ |
| |
| /* Number of bits in a limb. */ |
| #define LIMB_BITS (MANT_DIG / 2) |
| |
| /* 2^LIMB_BITS. */ |
| static const DOUBLE TWO_LIMB_BITS = |
| /* Assume LIMB_BITS <= 3 * 31. |
| Use the identity |
| n = floor(n/3) + floor((n+1)/3) + floor((n+2)/3). */ |
| (DOUBLE) (1U << (LIMB_BITS / 3)) |
| * (DOUBLE) (1U << ((LIMB_BITS + 1) / 3)) |
| * (DOUBLE) (1U << ((LIMB_BITS + 2) / 3)); |
| |
| /* 2^-LIMB_BITS. */ |
| static const DOUBLE TWO_LIMB_BITS_INVERSE = |
| /* Assume LIMB_BITS <= 3 * 31. |
| Use the identity |
| n = floor(n/3) + floor((n+1)/3) + floor((n+2)/3). */ |
| L_(1.0) / ((DOUBLE) (1U << (LIMB_BITS / 3)) |
| * (DOUBLE) (1U << ((LIMB_BITS + 1) / 3)) |
| * (DOUBLE) (1U << ((LIMB_BITS + 2) / 3))); |
| |
| DOUBLE |
| FMOD (DOUBLE x, DOUBLE y) |
| { |
| if (isfinite (x) && isfinite (y) && y != L_(0.0)) |
| { |
| if (x == L_(0.0)) |
| /* Return x, regardless of the sign of y. */ |
| return x; |
| |
| { |
| int negate = ((!signbit (x)) ^ (!signbit (y))); |
| |
| /* Take the absolute value of x and y. */ |
| x = FABS (x); |
| y = FABS (y); |
| |
| /* Trivial case that requires no computation. */ |
| if (x < y) |
| return (negate ? - x : x); |
| |
| { |
| int yexp; |
| DOUBLE ym; |
| DOUBLE y1; |
| DOUBLE y0; |
| int k; |
| DOUBLE x2; |
| DOUBLE x1; |
| DOUBLE x0; |
| |
| /* Write y = 2^yexp * (y1 * 2^-LIMB_BITS + y0 * 2^-(2*LIMB_BITS)) |
| where y1 is an integer, 2^(LIMB_BITS-1) <= y1 < 2^LIMB_BITS, |
| y1 has at most LIMB_BITS bits, |
| 0 <= y0 < 2^LIMB_BITS, |
| y0 has at most (MANT_DIG + 1) / 2 bits. */ |
| ym = FREXP (y, &yexp); |
| ym = ym * TWO_LIMB_BITS; |
| y1 = TRUNC (ym); |
| y0 = (ym - y1) * TWO_LIMB_BITS; |
| |
| /* Write |
| x = 2^(yexp+(k-3)*LIMB_BITS) |
| * (x2 * 2^(2*LIMB_BITS) + x1 * 2^LIMB_BITS + x0) |
| where x2, x1, x0 are each integers >= 0, < 2^LIMB_BITS. */ |
| { |
| int xexp; |
| DOUBLE xm = FREXP (x, &xexp); |
| /* Since we know x >= y, we know xexp >= yexp. */ |
| xexp -= yexp; |
| /* Compute k = ceil(xexp / LIMB_BITS). */ |
| k = (xexp + LIMB_BITS - 1) / LIMB_BITS; |
| /* Note: (k - 1) * LIMB_BITS + 1 <= xexp <= k * LIMB_BITS. */ |
| /* Note: 0.5 <= xm < 1.0. */ |
| xm = LDEXP (xm, xexp - (k - 1) * LIMB_BITS); |
| /* Note: Now xm < 2^(xexp - (k - 1) * LIMB_BITS) <= 2^LIMB_BITS |
| and xm >= 0.5 * 2^(xexp - (k - 1) * LIMB_BITS) >= 1.0 |
| and xm has at most MANT_DIG <= 2*LIMB_BITS+1 bits. */ |
| x2 = TRUNC (xm); |
| x1 = (xm - x2) * TWO_LIMB_BITS; |
| /* Split off x0 from x1 later. */ |
| } |
| |
| /* Test whether [x2,x1,0] >= 2^LIMB_BITS * [y1,y0]. */ |
| if (x2 > y1 || (x2 == y1 && x1 >= y0)) |
| { |
| /* Subtract 2^LIMB_BITS * [y1,y0] from [x2,x1,0]. */ |
| x2 -= y1; |
| x1 -= y0; |
| if (x1 < L_(0.0)) |
| { |
| if (!(x2 >= L_(1.0))) |
| abort (); |
| x2 -= L_(1.0); |
| x1 += TWO_LIMB_BITS; |
| } |
| } |
| |
| /* Split off x0 from x1. */ |
| { |
| DOUBLE x1int = TRUNC (x1); |
| x0 = TRUNC ((x1 - x1int) * TWO_LIMB_BITS); |
| x1 = x1int; |
| } |
| |
| for (; k > 0; k--) |
| { |
| /* Multiprecision division of the limb sequence [x2,x1,x0] |
| by [y1,y0]. */ |
| /* Here [x2,x1,x0] < 2^LIMB_BITS * [y1,y0]. */ |
| /* The first guess takes into account only [x2,x1] and [y1]. |
| |
| By Knuth's theorem, we know that |
| q* = min (floor ([x2,x1] / [y1]), 2^LIMB_BITS - 1) |
| and |
| q = floor ([x2,x1,x0] / [y1,y0]) |
| are not far away: |
| q* - 2 <= q <= q* + 1. |
| |
| Proof: |
| a) q* * y1 <= floor ([x2,x1] / [y1]) * y1 <= [x2,x1]. |
| Hence |
| [x2,x1,x0] - q* * [y1,y0] |
| = 2^LIMB_BITS * ([x2,x1] - q* * [y1]) + x0 - q* * y0 |
| >= x0 - q* * y0 |
| >= - q* * y0 |
| > - 2^(2*LIMB_BITS) |
| >= - 2 * [y1,y0] |
| So |
| [x2,x1,x0] > (q* - 2) * [y1,y0]. |
| b) If q* = floor ([x2,x1] / [y1]), then |
| [x2,x1] < (q* + 1) * y1 |
| Hence |
| [x2,x1,x0] - q* * [y1,y0] |
| = 2^LIMB_BITS * ([x2,x1] - q* * [y1]) + x0 - q* * y0 |
| <= 2^LIMB_BITS * (y1 - 1) + x0 - q* * y0 |
| <= 2^LIMB_BITS * (2^LIMB_BITS-2) + (2^LIMB_BITS-1) - 0 |
| < 2^(2*LIMB_BITS) |
| <= 2 * [y1,y0] |
| So |
| [x2,x1,x0] < (q* + 2) * [y1,y0]. |
| and so |
| q < q* + 2 |
| which implies |
| q <= q* + 1. |
| In the other case, q* = 2^LIMB_BITS - 1. Then trivially |
| q < 2^LIMB_BITS = q* + 1. |
| |
| We know that floor ([x2,x1] / [y1]) >= 2^LIMB_BITS if and |
| only if x2 >= y1. */ |
| DOUBLE q = |
| (x2 >= y1 |
| ? TWO_LIMB_BITS - L_(1.0) |
| : TRUNC ((x2 * TWO_LIMB_BITS + x1) / y1)); |
| if (q > L_(0.0)) |
| { |
| /* Compute |
| [x2,x1,x0] - q* * [y1,y0] |
| = 2^LIMB_BITS * ([x2,x1] - q* * [y1]) + x0 - q* * y0. */ |
| DOUBLE q_y1 = q * y1; /* exact, at most 2*LIMB_BITS bits */ |
| DOUBLE q_y1_1 = TRUNC (q_y1 * TWO_LIMB_BITS_INVERSE); |
| DOUBLE q_y1_0 = q_y1 - q_y1_1 * TWO_LIMB_BITS; |
| DOUBLE q_y0 = q * y0; /* exact, at most MANT_DIG bits */ |
| DOUBLE q_y0_1 = TRUNC (q_y0 * TWO_LIMB_BITS_INVERSE); |
| DOUBLE q_y0_0 = q_y0 - q_y0_1 * TWO_LIMB_BITS; |
| x2 -= q_y1_1; |
| x1 -= q_y1_0; |
| x1 -= q_y0_1; |
| x0 -= q_y0_0; |
| /* Move negative carry from x0 to x1 and from x1 to x2. */ |
| if (x0 < L_(0.0)) |
| { |
| x0 += TWO_LIMB_BITS; |
| x1 -= L_(1.0); |
| } |
| if (x1 < L_(0.0)) |
| { |
| x1 += TWO_LIMB_BITS; |
| x2 -= L_(1.0); |
| if (x1 < L_(0.0)) /* not sure this can happen */ |
| { |
| x1 += TWO_LIMB_BITS; |
| x2 -= L_(1.0); |
| } |
| } |
| if (x2 < L_(0.0)) |
| { |
| /* Reduce q by 1. */ |
| x1 += y1; |
| x0 += y0; |
| /* Move overflow from x0 to x1 and from x1 to x0. */ |
| if (x0 >= TWO_LIMB_BITS) |
| { |
| x0 -= TWO_LIMB_BITS; |
| x1 += L_(1.0); |
| } |
| if (x1 >= TWO_LIMB_BITS) |
| { |
| x1 -= TWO_LIMB_BITS; |
| x2 += L_(1.0); |
| } |
| if (x2 < L_(0.0)) |
| { |
| /* Reduce q by 1 again. */ |
| x1 += y1; |
| x0 += y0; |
| /* Move overflow from x0 to x1 and from x1 to x0. */ |
| if (x0 >= TWO_LIMB_BITS) |
| { |
| x0 -= TWO_LIMB_BITS; |
| x1 += L_(1.0); |
| } |
| if (x1 >= TWO_LIMB_BITS) |
| { |
| x1 -= TWO_LIMB_BITS; |
| x2 += L_(1.0); |
| } |
| if (x2 < L_(0.0)) |
| /* Shouldn't happen, because we proved that |
| q >= q* - 2. */ |
| abort (); |
| } |
| } |
| } |
| if (x2 > L_(0.0) |
| || x1 > y1 |
| || (x1 == y1 && x0 >= y0)) |
| { |
| /* Increase q by 1. */ |
| x1 -= y1; |
| x0 -= y0; |
| /* Move negative carry from x0 to x1 and from x1 to x2. */ |
| if (x0 < L_(0.0)) |
| { |
| x0 += TWO_LIMB_BITS; |
| x1 -= L_(1.0); |
| } |
| if (x1 < L_(0.0)) |
| { |
| x1 += TWO_LIMB_BITS; |
| x2 -= L_(1.0); |
| } |
| if (x2 < L_(0.0)) |
| abort (); |
| if (x2 > L_(0.0) |
| || x1 > y1 |
| || (x1 == y1 && x0 >= y0)) |
| /* Shouldn't happen, because we proved that |
| q <= q* + 1. */ |
| abort (); |
| } |
| /* Here [x2,x1,x0] < [y1,y0]. */ |
| /* Next round. */ |
| x2 = x1; |
| #if (MANT_DIG + 1) / 2 > LIMB_BITS /* y0 can have a fractional bit */ |
| x1 = TRUNC (x0); |
| x0 = (x0 - x1) * TWO_LIMB_BITS; |
| #else |
| x1 = x0; |
| x0 = L_(0.0); |
| #endif |
| /* Here [x2,x1,x0] < 2^LIMB_BITS * [y1,y0]. */ |
| } |
| /* Here k = 0. |
| The result is |
| 2^(yexp-3*LIMB_BITS) |
| * (x2 * 2^(2*LIMB_BITS) + x1 * 2^LIMB_BITS + x0). */ |
| { |
| DOUBLE r = |
| LDEXP ((x2 * TWO_LIMB_BITS + x1) * TWO_LIMB_BITS + x0, |
| yexp - 3 * LIMB_BITS); |
| return (negate ? - r : r); |
| } |
| } |
| } |
| } |
| else |
| { |
| if (ISNAN (x) || ISNAN (y)) |
| return x + y; /* NaN */ |
| else if (isinf (y)) |
| return x; |
| else |
| /* x infinite or y zero */ |
| return NAN; |
| } |
| } |
