1 /* Test of remainder*() function family.
2 Copyright (C) 2012 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
18 my_ldexp (DOUBLE x, int d)
32 const DOUBLE TWO_MANT_DIG =
33 /* Assume MANT_DIG <= 5 * 31.
35 n = floor(n/5) + floor((n+1)/5) + ... + floor((n+4)/5). */
36 (DOUBLE) (1U << ((MANT_DIG - 1) / 5))
37 * (DOUBLE) (1U << ((MANT_DIG - 1 + 1) / 5))
38 * (DOUBLE) (1U << ((MANT_DIG - 1 + 2) / 5))
39 * (DOUBLE) (1U << ((MANT_DIG - 1 + 3) / 5))
40 * (DOUBLE) (1U << ((MANT_DIG - 1 + 4) / 5));
42 /* Randomized tests. */
43 for (i = 0; i < SIZEOF (RANDOM) / 5; i++)
44 for (j = 0; j < SIZEOF (RANDOM) / 5; j++)
46 DOUBLE x = L_(16.0) * RANDOM[i]; /* 0.0 <= x <= 16.0 */
47 DOUBLE y = RANDOM[j]; /* 0.0 <= y < 1.0 */
50 DOUBLE z = REMAINDER (x, y);
51 ASSERT (z >= - L_(0.5) * y);
52 ASSERT (z <= L_(0.5) * y);
53 z -= x - (int) ((L_(2.0) * x + y) / (L_(2.0) * y)) * y;
54 ASSERT (/* The common case. */
55 (z > - L_(2.0) * L_(16.0) / TWO_MANT_DIG
56 && z < L_(2.0) * L_(16.0) / TWO_MANT_DIG)
57 || /* rounding error: 2x+y / 2y computed too large */
58 (z > y - L_(2.0) * L_(16.0) / TWO_MANT_DIG
59 && z < y + L_(2.0) * L_(16.0) / TWO_MANT_DIG)
60 || /* rounding error: 2x+y / 2y computed too small */
61 (z > - y - L_(2.0) * L_(16.0) / TWO_MANT_DIG
62 && z < - y + L_(2.0) * L_(16.0) / TWO_MANT_DIG));
66 for (i = 0; i < SIZEOF (RANDOM) / 5; i++)
67 for (j = 0; j < SIZEOF (RANDOM) / 5; j++)
69 DOUBLE x = L_(1.0e9) * RANDOM[i]; /* 0.0 <= x <= 10^9 */
70 DOUBLE y = RANDOM[j]; /* 0.0 <= y < 1.0 */
73 DOUBLE z = REMAINDER (x, y);
75 ASSERT (z >= - L_(0.5) * y);
76 ASSERT (z <= L_(0.5) * y);
78 /* Determine the quotient 2x+y / 2y in two steps, because it
80 int q1 = (int) (x / y / L_(65536.0));
81 int q2 = (int) ((L_(2.0) * (x - q1 * L_(65536.0) * y) + y)
83 DOUBLE q = (DOUBLE) q1 * L_(65536.0) + (DOUBLE) q2;
86 /* The absolute error of z can be up to 1e9/2^MANT_DIG.
87 The absolute error of r can also be up to 1e9/2^MANT_DIG.
88 Therefore the error of z - r can be twice as large. */
90 ASSERT (/* The common case. */
91 (z > - L_(2.0) * L_(1.0e9) / TWO_MANT_DIG
92 && z < L_(2.0) * L_(1.0e9) / TWO_MANT_DIG)
93 || /* rounding error: 2x+y / 2y computed too large */
94 (z > y - L_(2.0) * L_(1.0e9) / TWO_MANT_DIG
95 && z < y + L_(2.0) * L_(1.0e9) / TWO_MANT_DIG)
96 || /* rounding error: 2x+y / 2y computed too small */
97 (z > - y - L_(2.0) * L_(1.0e9) / TWO_MANT_DIG
98 && z < - y + L_(2.0) * L_(1.0e9) / TWO_MANT_DIG));
103 int large_exp = (MAX_EXP - 1 < 1000 ? MAX_EXP - 1 : 1000);
104 DOUBLE large = my_ldexp (L_(1.0), large_exp); /* = 2^large_exp */
105 for (i = 0; i < SIZEOF (RANDOM) / 10; i++)
106 for (j = 0; j < SIZEOF (RANDOM) / 10; j++)
108 DOUBLE x = large * RANDOM[i]; /* 0.0 <= x <= 2^large_exp */
109 DOUBLE y = RANDOM[j]; /* 0.0 <= y < 1.0 */
112 DOUBLE z = REMAINDER (x, y);
113 /* Regardless how large the rounding errors are, the result
114 must be >= -y/2, <= y/2. */
115 ASSERT (z >= - L_(0.5) * y);
116 ASSERT (z <= L_(0.5) * y);