mirror of
https://https.git.savannah.gnu.org/git/gnulib.git
synced 2026-04-28 06:33:36 +00:00
361 lines
12 KiB
C
361 lines
12 KiB
C
/* Remainder.
|
|
Copyright (C) 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/>. */
|
|
|
|
#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. */
|
|
#if defined _MSC_VER && !defined __clang__
|
|
# 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;
|
|
}
|
|
}
|