Dear Sir/Madam,
I am writing to inform you the problem of pow() floating point
imprecision and am trying to give solution to tackle this problem.
As mentioned in here [1], there is a floating point imprecision problem
in the pow function.
Attached is the patch files to fix this problem and the following is the
explanation of the problem and suggested solution.
Reason for Causing pow() Floating Point Imprecision
The main reason for the imprecision is located in line 126 to 130 in
"mingw-w64/mingw-w64-crt/math/powi.def.h" [2].
if (y < 0)
{
d = __FLT_CST(1.0) / d;
y = -y;
}
When calculating 'x ^ y', the program will inverse the value of 'x' and
negate the value of 'y' every time the program see a reciprocal (y < 0).
This is equivalent as:
x ^ (-y) , given y > 0
= (1 / x) ^ y
This is correct mathematically. However, the code in line 126 to 130
will give a imprecision value in computer memory because of limited size
of data type.
Consider the following example:
3. ^ -5
In mingw [2] current algorithm, it will perform
x = 1 / 3
y = 5
thus, (1 / 3) ^ 5
Again, this is correct mathematically. However, the value of (1 / 3) is
in fact 0.333333...
Computer cannot represent the value of '1 / 3' precisely in computer
memory under IEEE 754 because of limited size of data type.
This function yield a imprecise value in the very first step, and the
imprecision problem will mostly be enlarged in the future calculation.
We can observe the following result in different environment and
compiler.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main()
{
union {
double d;
unsigned long long int i;
} cast;
double a = 3.;
int b = -5;
printf("pow(%.2f, %d)\n", a, b);
cast.d = pow(a, b);
printf("%llx\n", cast.i);
printf("%.17e\n\n", cast.d);
return 0;
}
mingw pow [2]:
pow(3.00, -5)
3f70db20a88f4695
4.11522633744855915e-03
llvm pow [3]:
pow(3.00, -5)
3f70db20a88f4696
4.11522633744856002e-03
pow()
under NixOS 22.05 (Quokka),
gcc (GCC) 11.3.0
Copyright (C) 2021 Free Software Foundation, Inc.:
pow(3.00, -5)
3f70db20a88f4696
4.11522633744856002e-03
I believe the reason for yielding an different result is because of line
126 to 130 in "mingw-w64/mingw-w64-crt/math/powi.def.h".
And the reason for causing imprecision mentioned in this post [1] is
also because of the same reason. '10' can be represented precisely under
IEEE 754 but '0.1' cannot.
However, llvm pow() [3] can give a result align with the linux pow
although llvm [3] and mingw [2] are using the same algorithm which is
exponential by squaring.
Therefore, I suggest mingw to follow the algorithm of llvm [3] which
gives a more accurate result.
It is believe that the current alogrithm in mingw powi() would like to
preserve the value when the result of 'x^y' is near the the subnormal.
Instead of performing '1 / (x^y)' to give a more precise result, mingw
current algorithm would like to perform '(1 / x)^y' to give a more
comprehensive result because '1 / x' is less than 1 and will not
overflow.
However, even the value of 'x' is '0<x<1' (cannot overflow), the current
algorithm [2] of mingw will still operate 'x <- 1 / x' and 'y <- (-y)'.
This is totally unnecessary and will give a more imprecise value after
this operatoin (line 126 to 130). Since floating 'x' cannot be
represented precisely in computer already, and '1 / x' will give a more
imprecise value most of the time.
Again, I suggest mingw to follow the algorithm of llvm [3] which gives a
more accurate result.
Potential problem when result near subnormal
The potential problem of llvm [3] is that llvm [3] algorithm will result
'0' when the result of 'x^y' is near subnormal.
The maximum value that a `double` may represent is `0x1p+1023`, while
the minimum positive (subnormal) value is `0x1p-1074`
powi(2, -1025) = 0x1p-1025
And llvm [3] gives:
powi(2, -1025) = 1 / powi(2, 1025)
= 1 / +infinity
= 0
Suggested Solution
The problem of the current mingw pow is that it will do reciprocal
whenever the program see a 'y < 0'.
Even '0 < x < 1' and will not cause an overflow result for
'abs(x)^abs(y)', the current mingw pow [2] will still do the reciprocal
in the very beginning.
In fact, it is unnecessary and producing imprecision.
Thus, my suggestion is to add an overflow detection code in the
algorithm base on llvm algorithm [3].
The program should execute 'x <- 1 / x' and 'y -< -y' only if the
program detect 'y < 0' and detect an overflow result of 'abs(x) ^
abs(y)'.
The overflow detection code comes from the following fact:
x^y > m , given x > 0, y > 0 and m = maximum value of
specific data type.
log(x^y) > log(m)
y * log(x) > log(m)
After adding the overflow detection code, the algorithm should give
appropriate result in all range of specific data type.
Under the following situation (2.^-1025):
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main()
{
union {
double d;
unsigned long long int i;
} cast;
double a = 2.;
int b = -1025;
printf("pow(%.2f, %d)\n", a, b);
cast.d = pow(a, b);
printf("%llx\n", cast.i);
printf("%.17e\n\n", cast.d);
return 0;
}
my suggested solution:
pow(2.00, -1025)
2000000000000
2.78134232313400173e-309
mingw [2] pow():
pow(2.00, -1025)
2000000000000
2.78134232313400173e-309
llvm [3] pow():
pow(2.00, -1025)
0
0.00000000000000000e+00
pow()
under NixOS 22.05 (Quokka),
gcc (GCC) 11.3.0
Copyright (C) 2021 Free Software Foundation, Inc.:
pow(2.00, -1025)
2000000000000
2.78134232313400173e-309
Also, a more accurate result under some common case (10.^-9):
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main()
{
union {
double d;
unsigned long long int i;
} cast;
double a = 10.;
int b = -9;
printf("pow(%.2f, %d)\n", a, b);
cast.d = pow(a, b);
printf("%llx\n", cast.i);
printf("%.17e\n\n", cast.d);
return 0;
}
my suggested solution:
pow(10.00, -9)
3e112e0be826d695
1.00000000000000006e-09
mingw [2] pow():
pow(10.00, -9)
3e112e0be826d699
1.00000000000000089e-09
llvm [3] pow():
pow(10.00, -9)
3e112e0be826d695
1.00000000000000006e-09
pow()
under NixOS 22.05 (Quokka),
gcc (GCC) 11.3.0
Copyright (C) 2021 Free Software Foundation, Inc.:
pow(10.00, -9)
3e112e0be826d695
1.00000000000000006e-09
The suggested solution should solve the problem of both imprecision and
result near subnormal.
Hope the recommendation can be accepted.
Best,
Gary
Links:
------
[1] https://github.com/msys2/MINGW-packages/issues/11733
[2]
https://github.com/mingw-w64/mingw-w64/blob/0b02ea1d7dbf8afbaced7122a8a2fa715e0db0cf/mingw-w64-crt/math/powi.def.h#L126
[3]
https://github.com/llvm/llvm-project/blob/1483fb33b314ae02ec96b735c5ff15ce595322bf/compiler-rt/lib/builtins/powidf2.cFrom 53482aa782df3405edc952421ce11413063cdb13 Mon Sep 17 00:00:00 2001
From: Gary Wan <[email protected]>
Date: Mon, 6 Jun 2022 11:01:03 +0800
Subject: [PATCH 1/5] fix pow() floating point imprecision
---
mingw-w64-crt/math/powi.def.h | 39 ++++++++++++-----------------------
1 file changed, 13 insertions(+), 26 deletions(-)
diff --git a/mingw-w64-crt/math/powi.def.h b/mingw-w64-crt/math/powi.def.h
index 1997727fd..249fb64e4 100644
--- a/mingw-w64-crt/math/powi.def.h
+++ b/mingw-w64-crt/math/powi.def.h
@@ -121,33 +121,20 @@ __FLT_ABI(__powi) (__FLT_TYPE x, int y)
return (odd_y && signbit(x) ? -__FLT_HUGE_VAL : __FLT_HUGE_VAL);
}
- d = __FLT_ABI(fabs) (x);
- if (y < 0)
- {
- d = __FLT_CST(1.0) / d;
- y = -y;
+ double a = x;
+ int b = y;
+
+ const int recip = b < 0;
+ double r = 1;
+ while (1) {
+ if (b & 1)
+ r *= a;
+ b /= 2;
+ if (b == 0)
+ break;
+ a *= a;
}
+ return recip ? 1 / r : r;
- if (!y)
- rslt = __FLT_CST(1.0);
- else if (y == 1)
- rslt = d;
- else
- {
- unsigned int u = (unsigned int) y;
- rslt = ((u & 1) != 0) ? d : __FLT_CST(1.0);
- u >>= 1;
- do
- {
- d *= d;
- if ((u & 1) != 0)
- rslt *= d;
- u >>= 1;
- }
- while (u > 0);
- }
- if (signbit (x) && odd_y)
- rslt = -rslt;
- return rslt;
}
--
2.36.0
From 8dfff0e31313ba2db40b881e2956a264aea96509 Mon Sep 17 00:00:00 2001
From: Gary Wan <[email protected]>
Date: Wed, 8 Jun 2022 11:08:52 +0800
Subject: [PATCH 2/5] fix pow() floating point imprecision, add credit
---
mingw-w64-crt/math/powi.def.h | 27 ++++++++++++++++++---------
1 file changed, 18 insertions(+), 9 deletions(-)
diff --git a/mingw-w64-crt/math/powi.def.h b/mingw-w64-crt/math/powi.def.h
index 249fb64e4..90f515420 100644
--- a/mingw-w64-crt/math/powi.def.h
+++ b/mingw-w64-crt/math/powi.def.h
@@ -122,18 +122,27 @@ __FLT_ABI(__powi) (__FLT_TYPE x, int y)
}
- double a = x;
- int b = y;
-
- const int recip = b < 0;
+ // credit: llvm/llvm-project
+ // This project uses source code from the files llvm-project/compiler-rt/lib/builtins/powidf2.c
+ // from https://github.com/llvm/llvm-project/blob/1483fb33b314ae02ec96b735c5ff15ce595322bf/compiler-rt/lib/builtins/powidf2.c
+ // Copyright (c) 2009-2015 by the contributors listed in CREDITS.TXT
+ // CREDITS.TXT: https://github.com/llvm/llvm-project/blob/1483fb33b314ae02ec96b735c5ff15ce595322bf/compiler-rt/CREDITS.TXT
+ // licensed under the Apache 2.0 license.
+ // Followed by the whole Apache 2.0 license text.
+ // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+ // See https://llvm.org/LICENSE.txt for license information.
+ // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+
+
+ const int recip = y < 0;
double r = 1;
while (1) {
- if (b & 1)
- r *= a;
- b /= 2;
- if (b == 0)
+ if (y & 1)
+ r *= x;
+ y /= 2;
+ if (y == 0)
break;
- a *= a;
+ x *= x;
}
return recip ? 1 / r : r;
--
2.36.0
From b67781c5325f335e45de7e8969c210411b3e8a3d Mon Sep 17 00:00:00 2001
From: Gary Wan <[email protected]>
Date: Fri, 10 Jun 2022 11:54:55 +0800
Subject: [PATCH 3/5] deal wtih situation near subnormal
---
mingw-w64-crt/math/powi.def.h | 33 ++++++++++++++++++++++++++++-----
1 file changed, 28 insertions(+), 5 deletions(-)
diff --git a/mingw-w64-crt/math/powi.def.h b/mingw-w64-crt/math/powi.def.h
index 90f515420..e4caf1292 100644
--- a/mingw-w64-crt/math/powi.def.h
+++ b/mingw-w64-crt/math/powi.def.h
@@ -136,14 +136,37 @@ __FLT_ABI(__powi) (__FLT_TYPE x, int y)
const int recip = y < 0;
double r = 1;
+ double a_ = a;
+ int b_ = b;
+
+
while (1) {
- if (y & 1)
- r *= x;
- y /= 2;
- if (y == 0)
+ if (b & 1)
+ r *= a;
+ b /= 2;
+ if (b == 0)
break;
- x *= x;
+ a *= a;
+ }
+
+ if (recip && fpclassify(r) == FP_INFINITE)
+ {
+ printf("in\n");
+ a = 1 / a_;
+ b = abs(b_);
+ r = 1;
+
+ while (1) {
+ if (b & 1)
+ r *= a;
+ b /= 2;
+ if (b == 0)
+ break;
+ a *= a;
+ }
+ return r;
}
+
return recip ? 1 / r : r;
}
--
2.36.0
From d07487dea08a14a03cedac84ef53187acf97525d Mon Sep 17 00:00:00 2001
From: Gary Wan <[email protected]>
Date: Fri, 10 Jun 2022 16:55:21 +0800
Subject: [PATCH 4/5] fix pow() imprecision and deal with result near subnormal
---
mingw-w64-crt/math/powi.def.h | 41 +++++++++++++----------------------
1 file changed, 15 insertions(+), 26 deletions(-)
diff --git a/mingw-w64-crt/math/powi.def.h b/mingw-w64-crt/math/powi.def.h
index e4caf1292..355de7ebc 100644
--- a/mingw-w64-crt/math/powi.def.h
+++ b/mingw-w64-crt/math/powi.def.h
@@ -76,7 +76,7 @@ __FLT_ABI(__powi) (__FLT_TYPE x, int y)
{
int x_class = fpclassify (x);
int odd_y = y & 1;
- __FLT_TYPE d, rslt;
+ __FLT_TYPE rslt;
if (y == 0 || x == __FLT_CST(1.0))
return __FLT_CST(1.0);
@@ -135,38 +135,27 @@ __FLT_ABI(__powi) (__FLT_TYPE x, int y)
const int recip = y < 0;
- double r = 1;
- double a_ = a;
- int b_ = b;
+ __FLT_TYPE r = 1;
+ int overflow = abs(y) * log(__FLT_ABI(fabs) (x)) > (sizeof(__FLT_TYPE) * 8 - 1) * __FLT_LOGE2;
+ if (recip && overflow)
+ {
+ x = 1 / x;
+ y = -y;
+ }
while (1) {
- if (b & 1)
- r *= a;
- b /= 2;
- if (b == 0)
+ if (y & 1)
+ r *= x;
+ y /= 2;
+ if (y == 0)
break;
- a *= a;
+ x *= x;
}
- if (recip && fpclassify(r) == FP_INFINITE)
- {
- printf("in\n");
- a = 1 / a_;
- b = abs(b_);
- r = 1;
-
- while (1) {
- if (b & 1)
- r *= a;
- b /= 2;
- if (b == 0)
- break;
- a *= a;
- }
+ if (recip && overflow)
return r;
- }
-
+
return recip ? 1 / r : r;
}
--
2.36.0
From 60f2d2bb52f8154bbd1b3596a35980864866cbc9 Mon Sep 17 00:00:00 2001
From: Gary Wan <[email protected]>
Date: Fri, 10 Jun 2022 17:35:01 +0800
Subject: [PATCH 5/5] modify return code
---
mingw-w64-crt/math/powi.def.h | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
diff --git a/mingw-w64-crt/math/powi.def.h b/mingw-w64-crt/math/powi.def.h
index 355de7ebc..c98d107e8 100644
--- a/mingw-w64-crt/math/powi.def.h
+++ b/mingw-w64-crt/math/powi.def.h
@@ -153,9 +153,9 @@ __FLT_ABI(__powi) (__FLT_TYPE x, int y)
x *= x;
}
- if (recip && overflow)
- return r;
-
- return recip ? 1 / r : r;
-
+ if(recip) {
+ return overflow ? r : 1/r;
+ } else {
+ return r;
+}
}
--
2.36.0
_______________________________________________
Mingw-w64-public mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
