Mercurial > hg > octave-nkf > gnulib-hg

annotate lib/fmod.c @ 17848:ab58d4870664

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
version-etc: new year * doc/gnulib.texi: * lib/version-etc.c (COPYRIGHT_YEAR): Update copyright date. * all files: Run 'make update-copyright'.
author Paul Eggert <eggert@cs.ucla.edu>
date Thu, 01 Jan 2015 01:38:23 +0000
parents 344018b6e5d7
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
16486
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
1 /* Remainder.
17848
ab58d4870664 version-etc: new year
Paul Eggert <eggert@cs.ucla.edu>
parents: 17587
diff changeset
2 Copyright (C) 2011-2015 Free Software Foundation, Inc.
16486
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
3
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
4 This program is free software: you can redistribute it and/or modify
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
5 it under the terms of the GNU General Public License as published by
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
6 the Free Software Foundation; either version 3 of the License, or
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
7 (at your option) any later version.
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
8
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
9 This program is distributed in the hope that it will be useful,
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
12 GNU General Public License for more details.
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
13
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
14 You should have received a copy of the GNU General Public License
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
16
16565
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
17 #if ! defined USE_LONG_DOUBLE
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
18 # include <config.h>
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
19 #endif
16486
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
20
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
21 /* Specification. */
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
22 #include <math.h>
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
23
16565
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
24 #include <float.h>
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
25 #include <stdlib.h>
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
26
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
27 #ifdef USE_LONG_DOUBLE
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
28 # define FMOD fmodl
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
29 # define DOUBLE long double
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
30 # define MANT_DIG LDBL_MANT_DIG
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
31 # define L_(literal) literal##L
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
32 # define FREXP frexpl
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
33 # define LDEXP ldexpl
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
34 # define FABS fabsl
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
35 # define TRUNC truncl
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
36 # define ISNAN isnanl
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
37 #else
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
38 # define FMOD fmod
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
39 # define DOUBLE double
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
40 # define MANT_DIG DBL_MANT_DIG
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
41 # define L_(literal) literal
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
42 # define FREXP frexp
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
43 # define LDEXP ldexp
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
44 # define FABS fabs
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
45 # define TRUNC trunc
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
46 # define ISNAN isnand
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
47 #endif
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
48
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
49 /* MSVC with option -fp:strict refuses to compile constant initializers that
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
50 contain floating-point operations. Pacify this compiler. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
51 #ifdef _MSC_VER
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
52 # pragma fenv_access (off)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
53 #endif
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
54
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
55 #undef NAN
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
56 #if defined _MSC_VER
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
57 static DOUBLE zero;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
58 # define NAN (zero / zero)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
59 #else
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
60 # define NAN (L_(0.0) / L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
61 #endif
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
62
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
63 /* To avoid rounding errors during the computation of x - q * y,
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
64 there are three possibilities:
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
65 - Use fma (- q, y, x).
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
66 - Have q be a single bit at a time, and compute x - q * y
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
67 through a subtraction.
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
68 - Have q be at most MANT_DIG/2 bits at a time, and compute
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
69 x - q * y by splitting y into two halves:
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
70 y = y1 * 2^(MANT_DIG/2) + y0
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
71 x - q * y = (x - q * y1 * 2^MANT_DIG/2) - q * y0.
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
72 The latter is probably the most efficient. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
73
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
74 /* Number of bits in a limb. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
75 #define LIMB_BITS (MANT_DIG / 2)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
76
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
77 /* 2^LIMB_BITS. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
78 static const DOUBLE TWO_LIMB_BITS =
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
79 /* Assume LIMB_BITS <= 3 * 31.
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
80 Use the identity
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
81 n = floor(n/3) + floor((n+1)/3) + floor((n+2)/3). */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
82 (DOUBLE) (1U << (LIMB_BITS / 3))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
83 * (DOUBLE) (1U << ((LIMB_BITS + 1) / 3))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
84 * (DOUBLE) (1U << ((LIMB_BITS + 2) / 3));
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
85
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
86 /* 2^-LIMB_BITS. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
87 static const DOUBLE TWO_LIMB_BITS_INVERSE =
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
88 /* Assume LIMB_BITS <= 3 * 31.
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
89 Use the identity
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
90 n = floor(n/3) + floor((n+1)/3) + floor((n+2)/3). */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
91 L_(1.0) / ((DOUBLE) (1U << (LIMB_BITS / 3))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
92 * (DOUBLE) (1U << ((LIMB_BITS + 1) / 3))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
93 * (DOUBLE) (1U << ((LIMB_BITS + 2) / 3)));
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
94
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
95 DOUBLE
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
96 FMOD (DOUBLE x, DOUBLE y)
16486
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
97 {
16565
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
98 if (isfinite (x) && isfinite (y) && y != L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
99 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
100 if (x == L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
101 /* Return x, regardless of the sign of y. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
102 return x;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
103
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
104 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
105 int negate = ((!signbit (x)) ^ (!signbit (y)));
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
106
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
107 /* Take the absolute value of x and y. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
108 x = FABS (x);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
109 y = FABS (y);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
110
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
111 /* Trivial case that requires no computation. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
112 if (x < y)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
113 return (negate ? - x : x);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
114
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
115 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
116 int yexp;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
117 DOUBLE ym;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
118 DOUBLE y1;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
119 DOUBLE y0;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
120 int k;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
121 DOUBLE x2;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
122 DOUBLE x1;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
123 DOUBLE x0;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
124
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
125 /* Write y = 2^yexp * (y1 * 2^-LIMB_BITS + y0 * 2^-(2*LIMB_BITS))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
126 where y1 is an integer, 2^(LIMB_BITS-1) <= y1 < 2^LIMB_BITS,
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
127 y1 has at most LIMB_BITS bits,
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
128 0 <= y0 < 2^LIMB_BITS,
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
129 y0 has at most (MANT_DIG + 1) / 2 bits. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
130 ym = FREXP (y, &yexp);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
131 ym = ym * TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
132 y1 = TRUNC (ym);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
133 y0 = (ym - y1) * TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
134
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
135 /* Write
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
136 x = 2^(yexp+(k-3)*LIMB_BITS)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
137 * (x2 * 2^(2*LIMB_BITS) + x1 * 2^LIMB_BITS + x0)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
138 where x2, x1, x0 are each integers >= 0, < 2^LIMB_BITS. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
139 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
140 int xexp;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
141 DOUBLE xm = FREXP (x, &xexp);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
142 /* Since we know x >= y, we know xexp >= yexp. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
143 xexp -= yexp;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
144 /* Compute k = ceil(xexp / LIMB_BITS). */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
145 k = (xexp + LIMB_BITS - 1) / LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
146 /* Note: (k - 1) * LIMB_BITS + 1 <= xexp <= k * LIMB_BITS. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
147 /* Note: 0.5 <= xm < 1.0. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
148 xm = LDEXP (xm, xexp - (k - 1) * LIMB_BITS);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
149 /* Note: Now xm < 2^(xexp - (k - 1) * LIMB_BITS) <= 2^LIMB_BITS
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
150 and xm >= 0.5 * 2^(xexp - (k - 1) * LIMB_BITS) >= 1.0
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
151 and xm has at most MANT_DIG <= 2*LIMB_BITS+1 bits. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
152 x2 = TRUNC (xm);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
153 x1 = (xm - x2) * TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
154 /* Split off x0 from x1 later. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
155 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
156
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
157 /* Test whether [x2,x1,0] >= 2^LIMB_BITS * [y1,y0]. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
158 if (x2 > y1 || (x2 == y1 && x1 >= y0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
159 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
160 /* Subtract 2^LIMB_BITS * [y1,y0] from [x2,x1,0]. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
161 x2 -= y1;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
162 x1 -= y0;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
163 if (x1 < L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
164 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
165 if (!(x2 >= L_(1.0)))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
166 abort ();
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
167 x2 -= L_(1.0);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
168 x1 += TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
169 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
170 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
171
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
172 /* Split off x0 from x1. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
173 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
174 DOUBLE x1int = TRUNC (x1);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
175 x0 = TRUNC ((x1 - x1int) * TWO_LIMB_BITS);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
176 x1 = x1int;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
177 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
178
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
179 for (; k > 0; k--)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
180 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
181 /* Multiprecision division of the limb sequence [x2,x1,x0]
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
182 by [y1,y0]. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
183 /* Here [x2,x1,x0] < 2^LIMB_BITS * [y1,y0]. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
184 /* The first guess takes into account only [x2,x1] and [y1].
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
185
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
186 By Knuth's theorem, we know that
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
187 q* = min (floor ([x2,x1] / [y1]), 2^LIMB_BITS - 1)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
188 and
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
189 q = floor ([x2,x1,x0] / [y1,y0])
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
190 are not far away:
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
191 q* - 2 <= q <= q* + 1.
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
192
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
193 Proof:
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
194 a) q* * y1 <= floor ([x2,x1] / [y1]) * y1 <= [x2,x1].
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
195 Hence
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
196 [x2,x1,x0] - q* * [y1,y0]
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
197 = 2^LIMB_BITS * ([x2,x1] - q* * [y1]) + x0 - q* * y0
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
198 >= x0 - q* * y0
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
199 >= - q* * y0
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
200 > - 2^(2*LIMB_BITS)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
201 >= - 2 * [y1,y0]
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
202 So
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
203 [x2,x1,x0] > (q* - 2) * [y1,y0].
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
204 b) If q* = floor ([x2,x1] / [y1]), then
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
205 [x2,x1] < (q* + 1) * y1
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
206 Hence
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
207 [x2,x1,x0] - q* * [y1,y0]
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
208 = 2^LIMB_BITS * ([x2,x1] - q* * [y1]) + x0 - q* * y0
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
209 <= 2^LIMB_BITS * (y1 - 1) + x0 - q* * y0
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
210 <= 2^LIMB_BITS * (2^LIMB_BITS-2) + (2^LIMB_BITS-1) - 0
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
211 < 2^(2*LIMB_BITS)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
212 <= 2 * [y1,y0]
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
213 So
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
214 [x2,x1,x0] < (q* + 2) * [y1,y0].
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
215 and so
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
216 q < q* + 2
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
217 which implies
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
218 q <= q* + 1.
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
219 In the other case, q* = 2^LIMB_BITS - 1. Then trivially
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
220 q < 2^LIMB_BITS = q* + 1.
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
221
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
222 We know that floor ([x2,x1] / [y1]) >= 2^LIMB_BITS if and
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
223 only if x2 >= y1. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
224 DOUBLE q =
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
225 (x2 >= y1
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
226 ? TWO_LIMB_BITS - L_(1.0)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
227 : TRUNC ((x2 * TWO_LIMB_BITS + x1) / y1));
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
228 if (q > L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
229 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
230 /* Compute
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
231 [x2,x1,x0] - q* * [y1,y0]
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
232 = 2^LIMB_BITS * ([x2,x1] - q* * [y1]) + x0 - q* * y0. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
233 DOUBLE q_y1 = q * y1; /* exact, at most 2*LIMB_BITS bits */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
234 DOUBLE q_y1_1 = TRUNC (q_y1 * TWO_LIMB_BITS_INVERSE);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
235 DOUBLE q_y1_0 = q_y1 - q_y1_1 * TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
236 DOUBLE q_y0 = q * y0; /* exact, at most MANT_DIG bits */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
237 DOUBLE q_y0_1 = TRUNC (q_y0 * TWO_LIMB_BITS_INVERSE);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
238 DOUBLE q_y0_0 = q_y0 - q_y0_1 * TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
239 x2 -= q_y1_1;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
240 x1 -= q_y1_0;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
241 x1 -= q_y0_1;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
242 x0 -= q_y0_0;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
243 /* Move negative carry from x0 to x1 and from x1 to x2. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
244 if (x0 < L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
245 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
246 x0 += TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
247 x1 -= L_(1.0);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
248 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
249 if (x1 < L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
250 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
251 x1 += TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
252 x2 -= L_(1.0);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
253 if (x1 < L_(0.0)) /* not sure this can happen */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
254 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
255 x1 += TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
256 x2 -= L_(1.0);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
257 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
258 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
259 if (x2 < L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
260 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
261 /* Reduce q by 1. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
262 x1 += y1;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
263 x0 += y0;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
264 /* Move overflow from x0 to x1 and from x1 to x0. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
265 if (x0 >= TWO_LIMB_BITS)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
266 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
267 x0 -= TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
268 x1 += L_(1.0);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
269 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
270 if (x1 >= TWO_LIMB_BITS)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
271 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
272 x1 -= TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
273 x2 += L_(1.0);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
274 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
275 if (x2 < L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
276 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
277 /* Reduce q by 1 again. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
278 x1 += y1;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
279 x0 += y0;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
280 /* Move overflow from x0 to x1 and from x1 to x0. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
281 if (x0 >= TWO_LIMB_BITS)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
282 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
283 x0 -= TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
284 x1 += L_(1.0);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
285 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
286 if (x1 >= TWO_LIMB_BITS)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
287 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
288 x1 -= TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
289 x2 += L_(1.0);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
290 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
291 if (x2 < L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
292 /* Shouldn't happen, because we proved that
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
293 q >= q* - 2. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
294 abort ();
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
295 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
296 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
297 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
298 if (x2 > L_(0.0)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
299 || x1 > y1
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
300 || (x1 == y1 && x0 >= y0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
301 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
302 /* Increase q by 1. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
303 x1 -= y1;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
304 x0 -= y0;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
305 /* Move negative carry from x0 to x1 and from x1 to x2. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
306 if (x0 < L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
307 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
308 x0 += TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
309 x1 -= L_(1.0);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
310 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
311 if (x1 < L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
312 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
313 x1 += TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
314 x2 -= L_(1.0);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
315 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
316 if (x2 < L_(0.0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
317 abort ();
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
318 if (x2 > L_(0.0)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
319 || x1 > y1
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
320 || (x1 == y1 && x0 >= y0))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
321 /* Shouldn't happen, because we proved that
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
322 q <= q* + 1. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
323 abort ();
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
324 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
325 /* Here [x2,x1,x0] < [y1,y0]. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
326 /* Next round. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
327 x2 = x1;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
328 #if (MANT_DIG + 1) / 2 > LIMB_BITS /* y0 can have a fractional bit */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
329 x1 = TRUNC (x0);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
330 x0 = (x0 - x1) * TWO_LIMB_BITS;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
331 #else
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
332 x1 = x0;
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
333 x0 = L_(0.0);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
334 #endif
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
335 /* Here [x2,x1,x0] < 2^LIMB_BITS * [y1,y0]. */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
336 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
337 /* Here k = 0.
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
338 The result is
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
339 2^(yexp-3*LIMB_BITS)
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
340 * (x2 * 2^(2*LIMB_BITS) + x1 * 2^LIMB_BITS + x0). */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
341 {
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
342 DOUBLE r =
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
343 LDEXP ((x2 * TWO_LIMB_BITS + x1) * TWO_LIMB_BITS + x0,
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
344 yexp - 3 * LIMB_BITS);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
345 return (negate ? - r : r);
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
346 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
347 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
348 }
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
349 }
16486
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
350 else
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
351 {
16565
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
352 if (ISNAN (x) || ISNAN (y))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
353 return x + y; /* NaN */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
354 else if (isinf (y))
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
355 return x;
16486
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
356 else
16565
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
357 /* x infinite or y zero */
0badbd4d62e4 fmod, fmodl: Fix computation for large quotients x / y.
Bruno Haible <bruno@clisp.org>
parents: 16486
diff changeset
358 return NAN;
16486
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
359 }
8b3ead70232c fmod-ieee: Work around test failures on OSF/1, mingw.
Bruno Haible <bruno@clisp.org>
parents:
diff changeset
360 }