Skip to content

Commit 7f52241

Browse files
MitchellThompkinsMitchellThompkins
committed
T(1e-01) and T(1e-03) are double literals promoted to T, making them
slightly larger than the corresponding long double values 0.1L and 0.001L (which was prevoiusly used). This caused find_exponent to return the wrong exponent for those exact boundary values, e.g. log(0.1L) returned 2*ln(0.1). This actually effected all calls to gcem::log; it just only manifests when log lands exactly on 0.1 or 0.001 in long double. add test case for failing log case add 8 epsilon of error add tests for find_exponent make epsilon ULP choice configurable and add mantissa.hpp tests no reason to use T(8) wording updates fix fall-through for 10.0L and 100.0L
1 parent 2ff0096 commit 7f52241

6 files changed

Lines changed: 158 additions & 24 deletions

File tree

‎include/gcem_incl/find_exponent.hpp‎

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -34,15 +34,21 @@ llint_t
3434
find_exponent(const T x, const llint_t exponent)
3535
noexcept
3636
{
37+
// Subtract one epsilon from T(1) so that values whose x*10^n product lands
38+
// within one epsilon below 1.0 (e.g. 0.001L in 80-bit, which rounds to
39+
// ~0.5 eps below 1.0 after scaling by 1000) are not pushed into the wrong
40+
// decade. One epsilon is the unique correct choice: it must be > 0.5 eps
41+
// to catch 0.001L, and < 1.5 eps to avoid misclassifying the next
42+
// representable value below the boundary.
3743
return( // < 1
38-
x < T(1e-03) ? \
44+
x * T(1000) < T(1) - GCLIM<T>::epsilon() ? \
3945
find_exponent(x * T(1e+04), exponent - llint_t(4)) :
40-
x < T(1e-01) ? \
46+
x * T(10) < T(1) - GCLIM<T>::epsilon() ? \
4147
find_exponent(x * T(1e+02), exponent - llint_t(2)) :
42-
x < T(1) ? \
48+
x < T(1) - GCLIM<T>::epsilon() ? \
4349
find_exponent(x * T(10), exponent - llint_t(1)) :
4450
// > 10
45-
x > T(10) ? \
51+
x >= T(10) ? \
4652
find_exponent(x / T(10), exponent + llint_t(1)) :
4753
x > T(1e+02) ? \
4854
find_exponent(x / T(1e+02), exponent + llint_t(2)) :
@@ -52,6 +58,7 @@ noexcept
5258
exponent );
5359
}
5460

61+
5562
}
5663

5764
#endif

‎include/gcem_incl/mantissa.hpp‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -34,9 +34,13 @@ T
3434
mantissa(const T x)
3535
noexcept
3636
{
37-
return( x < T(1) ? \
38-
mantissa(x * T(10)) :
39-
x > T(10) ? \
37+
// Mirror the one-epsilon tolerance used in find_exponent so that
38+
// (mantissa(x), find_exponent(x)) stays a consistent decomposition of x.
39+
return( x < T(1) - GCLIM<T>::epsilon() ? \
40+
mantissa(x * T(10)) :
41+
x < T(1) ? \
42+
T(1) :
43+
x >= T(10) ? \
4044
mantissa(x / T(10)) :
4145
// else
4246
x );

‎tests/Makefile‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -101,6 +101,9 @@ fabsl:
101101
factorial:
102102
$(GCEM_MAKE_CALL)
103103

104+
find_exponent:
105+
$(GCEM_MAKE_CALL)
106+
104107
fmod:
105108
$(GCEM_MAKE_CALL)
106109

@@ -140,6 +143,9 @@ log_binomial_coef:
140143
log:
141144
$(GCEM_MAKE_CALL)
142145

146+
mantissa:
147+
$(GCEM_MAKE_CALL)
148+
143149
log1p:
144150
$(GCEM_MAKE_CALL)
145151

‎tests/find_exponent.cpp‎

Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,59 @@
1+
/*################################################################################
2+
##
3+
## Copyright (C) 2016-2026 Keith O'Hara
4+
##
5+
## This file is part of the GCE-Math C++ library.
6+
##
7+
## Licensed under the Apache License, Version 2.0 (the "License");
8+
## you may not use this file except in compliance with the License.
9+
## You may obtain a copy of the License at
10+
##
11+
## http://www.apache.org/licenses/LICENSE-2.0
12+
##
13+
## Unless required by applicable law or agreed to in writing, software
14+
## distributed under the License is distributed on an "AS IS" BASIS,
15+
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16+
## See the License for the specific language governing permissions and
17+
## limitations under the License.
18+
##
19+
################################################################################*/
20+
21+
#define TEST_PRINT_PRECISION_1 6
22+
#define TEST_PRINT_PRECISION_2 18
23+
24+
#include "gcem_tests.hpp"
25+
26+
int main()
27+
{
28+
print_begin("find_exponent");
29+
30+
// long double: positive exponents
31+
32+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, 0.0L, 1.0L, 0LL);
33+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, 0.0L, 5.0L, 0LL);
34+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, 1.0L, 10.0L, 0LL);
35+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, 1.0L, 11.0L, 0LL);
36+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, 1.0L, 99.9L, 0LL);
37+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, 2.0L, 100.0L, 0LL);
38+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, 2.0L, 200.0L, 0LL);
39+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, 2.0L, 123.45L, 0LL);
40+
41+
// long double: negative exponents, including the boundary cases (0.1L and
42+
// 0.001L round toward zero in 80-bit, placing them just below their nominal
43+
// decade boundary)
44+
45+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, -1.0L, 0.5L, 0LL);
46+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, -1.0L, 0.1L, 0LL);
47+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, -2.0L, 0.01L, 0LL);
48+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, -3.0L, 0.001L, 0LL);
49+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, -4.0L, 2.0e-4L, 0LL);
50+
51+
// double: verify the same boundary values work for double
52+
53+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, -1.0L, 0.1, 0LL);
54+
GCEM_TEST_EXPECTED_VAL(gcem::internal::find_exponent, -3.0L, 0.001, 0LL);
55+
56+
print_final("find_exponent");
57+
58+
return 0;
59+
}

‎tests/log.cpp‎

Lines changed: 19 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/*################################################################################
22
##
3-
## Copyright (C) 2016-2024 Keith O'Hara
3+
## Copyright (C) 2016-2026 Keith O'Hara
44
##
55
## This file is part of the GCE-Math C++ library.
66
##
@@ -31,23 +31,25 @@ int main()
3131

3232
//
3333

34-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, 0.5L);
35-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, 0.00199900000000000208L);
36-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, 1.0L);
37-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, 1.5L);
38-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, 41.5L);
39-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, 123456789.5L);
40-
41-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, 0.0L);
42-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, -1.0L);
43-
44-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, 1e-500L);
45-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, std::numeric_limits<long double>::min());
46-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, std::numeric_limits<double>::max());
34+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, 0.001L);
35+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, 0.1L);
36+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, 0.5L);
37+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, 0.00199900000000000208L);
38+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, 1.0L);
39+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, 1.5L);
40+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, 41.5L);
41+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, 123456789.5L);
42+
43+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, 0.0L);
44+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, -1.0L);
45+
46+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, 1e-500L);
47+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, std::numeric_limits<long double>::min());
48+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, std::numeric_limits<double>::max());
4749

48-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, -std::numeric_limits<long double>::infinity());
49-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, std::numeric_limits<long double>::infinity());
50-
GCEM_TEST_COMPARE_VALS(gcem::log,std::log, std::numeric_limits<long double>::quiet_NaN());
50+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, -std::numeric_limits<long double>::infinity());
51+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, std::numeric_limits<long double>::infinity());
52+
GCEM_TEST_COMPARE_VALS(gcem::log, std::log, std::numeric_limits<long double>::quiet_NaN());
5153

5254
//
5355

‎tests/mantissa.cpp‎

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,56 @@
1+
/*################################################################################
2+
##
3+
## Copyright (C) 2016-2026 Keith O'Hara
4+
##
5+
## This file is part of the GCE-Math C++ library.
6+
##
7+
## Licensed under the Apache License, Version 2.0 (the "License");
8+
## you may not use this file except in compliance with the License.
9+
## You may obtain a copy of the License at
10+
##
11+
## http://www.apache.org/licenses/LICENSE-2.0
12+
##
13+
## Unless required by applicable law or agreed to in writing, software
14+
## distributed under the License is distributed on an "AS IS" BASIS,
15+
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16+
## See the License for the specific language governing permissions and
17+
## limitations under the License.
18+
##
19+
################################################################################*/
20+
21+
#define TEST_PRINT_PRECISION_1 6
22+
#define TEST_PRINT_PRECISION_2 18
23+
24+
#include "gcem_tests.hpp"
25+
26+
int main()
27+
{
28+
print_begin("mantissa");
29+
30+
// long double: values whose mantissa is unambiguous
31+
32+
GCEM_TEST_EXPECTED_VAL(gcem::internal::mantissa, 1.0L, 1.0L);
33+
GCEM_TEST_EXPECTED_VAL(gcem::internal::mantissa, 5.0L, 5.0L);
34+
GCEM_TEST_EXPECTED_VAL(gcem::internal::mantissa, 9.99L, 9.99L);
35+
GCEM_TEST_EXPECTED_VAL(gcem::internal::mantissa, 2.0L, 20.0L);
36+
GCEM_TEST_EXPECTED_VAL(gcem::internal::mantissa, 1.2345L, 123.45L);
37+
GCEM_TEST_EXPECTED_VAL(gcem::internal::mantissa, 5.0L, 0.5L);
38+
GCEM_TEST_EXPECTED_VAL(gcem::internal::mantissa, 2.0L, 0.02L);
39+
40+
// long double: boundary cases withing the ULP-tolerant snap. 0.1L and
41+
// 0.001L round toward zero in 80-bit, so repeated scaling by 10 lands just
42+
// below 1.0 rather than exactly on it. mantissa() must snap these to 1.0
43+
// to stay consistent with find_exponent().
44+
45+
GCEM_TEST_EXPECTED_VAL(gcem::internal::mantissa, 1.0L, 0.1L);
46+
GCEM_TEST_EXPECTED_VAL(gcem::internal::mantissa, 1.0L, 0.001L);
47+
48+
// double: verify the same boundary values work for double
49+
50+
GCEM_TEST_EXPECTED_VAL(gcem::internal::mantissa, 1.0, 0.1);
51+
GCEM_TEST_EXPECTED_VAL(gcem::internal::mantissa, 1.0, 0.001);
52+
53+
print_final("mantissa");
54+
55+
return 0;
56+
}

0 commit comments

Comments
 (0)