Dear Sir/Madam,
As mentioned in here, there is a floating point imprecision problem in
the pow function.
Originally,
#include <stdio.h>
#include <math.h>
int main() {
union {
double d;
unsigned long long int i;
} cast;
cast.d = pow(10., -9);
printf("%llx\n", cast.i);
return 0;
}
3e112e0be826d699
the problem should be fixed after modification
3e112e0be826d695
this modification should also pass the special cases
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
double powi_wrong(double d, int y) {
if (y < 0) {
d = 1.0 / d;
y = -y;
}
unsigned int u = (unsigned int) y;
double rslt = ((u & 1) != 0) ? d : 1.0;
u >>= 1;
do {
d *= d;
if ((u & 1) != 0)
rslt *= d;
u >>= 1;
} while (u > 0);
return rslt;
}
double powi_correct(double a, int b) {
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;
}
int main()
{
printf("pow correct:\n%llx\n", powi_correct(10, -9));
printf("pow(10, -9):\n%llx\n", pow(10.0, -9));
printf("\nSpecial cases:\n");
printf("powi (x, +/-0) is 1 for any x (even a zero, quiet NaN, or
infinity)\n");
printf("pow(5., 0):\n%f\n", pow(5.,0));
printf("powi (+1, y) is 1 for any y (even a quiet NaN)\n");
printf("pow(1., 0):\n%f\n", pow(1.,0));
printf("powi (+/-0, y) is +/-oo and signals the divideByZero exception
for y an odd integer < 0\n");
printf("pow(0., -7):\n%f\n", pow(0.,-7));
printf("pow(-0., -7):\n%f\n", pow(-0.,-7));
printf("powi (+/-0, y) is +oo and signals the divideByZero exception
for finite y < 0 and not an odd integer\n");
printf("pow(0., -8):\n%f\n", pow(0., -8));
printf("pow(-0., -8):\n%f\n", pow(-0., -8));
printf("powi (+/-0, y) is +/-0 for finite y > 0 an odd integer\n");
printf("pow(0., 7):\n%f\n", pow(0., 7));
printf("pow(-0., 7):\n%f\n", pow(-0., 7));
printf("powi (+/-0, y) is +0 for finite y > 0 and not an odd
integer\n");
printf("pow(0., 8):\n%f\n", pow(0., 8));
printf("pow(-0., 8):\n%f\n", pow(0., 8));
long a = 10;
long b = -9;
double c = pow(a, b);
printf("10**-9 == 1e-9:\n%d\n", c == 1e-9);
printf("10**-9 is:\n%.17le\n", c);
printf("%llx", pow(10, -9));
}
result:
pow inaccurate:
3e112e0be826d699
pow(10, -9):
3e112e0be826d695
Special cases:
powi (x, +/-0) is 1 for any x (even a zero, quiet NaN, or infinity)
pow(5., 0):
1.000000
powi (+1, y) is 1 for any y (even a quiet NaN)
pow(1., 0):
1.000000
powi (+/-0, y) is +/-oo and signals the divideByZero exception for y an
odd integer < 0
pow(0., -7):
inf
pow(-0., -7):
-inf
powi (+/-0, y) is +oo and signals the divideByZero exception for finite
y < 0 and not an odd integer
pow(0., -8):
inf
pow(-0., -8):
inf
powi (+/-0, y) is +/-0 for finite y > 0 an odd integer
pow(0., 7):
0.000000
pow(-0., 7):
-0.000000
powi (+/-0, y) is +0 for finite y > 0 and not an odd integer
pow(0., 8):
0.000000
pow(-0., 8):
0.000000
10**-9 == 1e-9:
1
10**-9 is:
1.00000000000000006e-09
3e112e0be826d695
Attached is a patch to fix that.
Best,
Gary
From 53482aa782df3405edc952421ce11413063cdb13 Mon Sep 17 00:00:00 2001
From: Gary Wan <[email protected]>
Date: Mon, 6 Jun 2022 11:01:03 +0800
Subject: [PATCH] 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
_______________________________________________
Mingw-w64-public mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
