2

Legacy code I am to maintain has a function with a variable x of type double. This variable contains a position measured in meters. It is expected to hold a value in the mm range, let's say 0.0001 < x < 0.01.

Then I noticed a comment in the code like this:

// We will take the 6th power of this variable.
// Convert number to mm to get better accuracy.            
x /= 1.0e-3;
double y = pow(x, 6.0);
  • Does this make sense, is the precision of pow() indeed better when the first parameter is of the order of magnitude of 1? Given the way I think floating point numbers work, I would not expect the size of the mantissa to matter, as long as one is not running into the limits.
  • Should one scale down parameters considerable larger than 1 too?
  • Bonus question: does the value of the second parameter (the power) influence precision? I.e. is taking higher powers less accurate than small powers?

The target platform is Linux 3.10.x for x86_64 and we use the GCC compiler version 4.9.3.

I tried a little test program and it does indeed show larger (but still not very significant) errors when using parameters that are much larger or smaller than 1, but this might also be caused by how the inputs and/or results can be represented by floating point numbers?

#include <stdio.h>
#include <math.h>

void printpow(double d)
{
   double pow6 = pow(d, 6.0);
   printf("% 8g ^ 6 = %.20g\n", d, pow6);
}

int main(int argc, char **argv)
{
   
   printpow(7.0e-10);
   printpow(7.0e-6);
   printpow(7.0e-3);
   printpow(7.0);
   printpow(7.0e3);
   printpow(7.0e6);
   printpow(7.0e10);
   return 0;
}

Gives as output:

7e-10 ^ 6 = 1.1764899999999995812e-55
7e-06 ^ 6 = 1.1764899999999999629e-31
0.007 ^ 6 = 1.1764900000000001195e-13
    7 ^ 6 = 117649
 7000 ^ 6 = 1.1764900000000000603e+23
7e+06 ^ 6 = 1.1764899999999999889e+41
7e+10 ^ 6 = 1.1764900000000000914e+65
slingeraap
  • 484
  • 1
  • 5
  • 14
  • 1
    Relevant? [pow() function in C problems](https://stackoverflow.com/a/48368368/4142924) and [pow() function giving wrong answer](https://stackoverflow.com/questions/56522223/pow-function-giving-wrong-answer). Not all `pow()` implementations are as good as all others. – Weather Vane Nov 11 '20 at 16:25
  • 2
    Well, `pow(d, 6.0)` is basically doing `exp(6.0*log(d))`. The `log` operation is presumably going to be more accurate for smaller values, since the derivative increases as `d` approaches zero from above. For large values of `d`, where the derivative is close to zero, multiple values of `d` end up being mapped onto the same `log` value, so presumably some accuracy is lost. The `pow` function could presumably perform a similar magnitude shift internally, but there's no guarantee that it will go to such effort. – Tom Karzes Nov 11 '20 at 16:30
  • Maybe not relevant, but please don't confuse [accuracy with precision](https://manoa.hawaii.edu/exploringourfluidearth/physical/world-ocean/map-distortion/practices-science-precision-vs-accuracy#:~:text=Accuracy%20refers%20to%20how%20close,item%20are%20to%20each%20other.&text=That%20means%20it%20is%20possible,be%20accurate%20without%20being%20precise.). – Adrian Mole Nov 11 '20 at 16:33
  • 5
    The division by `1e-3` is itself going to introduce some error since the denominator cannot be represented exactly (assuming a binary floating point implementation). Perhaps it just happened to cancel out some other inaccuracies for the particular range of numbers the application was interested in. – Ian Abbott Nov 11 '20 at 16:34
  • 2
    @IanAbbott Yeah, the division doesn't really make much sense given that it could easily be replaced by `x *= 1.0e+3;` – Tom Karzes Nov 11 '20 at 16:36
  • What is your C implementation? What system are you running on? In particular, where does your standard math library come from? – Eric Postpischil Nov 11 '20 at 16:43
  • @EricPostpischil: this is running on Linux 3.10, x86, compiled with GCC 4.9.3. – slingeraap Nov 11 '20 at 16:54
  • Maybe the code does something with `y` afterwards that makes a difference to the accuracy of some final result. After all, that is why functions such as `expm1` and `log1p` exist. – Ian Abbott Nov 11 '20 at 17:17
  • 1
    BTW: `double y = pow(x, 6.0);` so how is `y` used afterwards? Is it multiplied by `1E-18` or ... ? – Support Ukraine Nov 11 '20 at 17:18
  • If you try calculating `double pow6 = pow(d, 6.0);`, `double pow6a = pow(d * 1e3, 6.0) / 1e18);`, `double pow6b = pow(d / 1e-3, 6.0) / 1e18;`, `double pow6c = pow(d * 1e3, 6.0) * 1e-18;` and `double pow6d = pow(d / 1e-3, 6.0) * 1e-18;`, and printing the results, you will probably find the accuracy varies a bit, but not in a particularly consistent way. – Ian Abbott Nov 11 '20 at 17:39
  • if your exponents are always integers then write a powint function yourself instead of using the standard pow that work with float exponents – phuclv Nov 12 '20 at 06:20
  • @TomKarzes: Re “The `log` operation is presumably going to be more accurate for smaller values”: There is no reason to make this presumption. `pow` routines do not spontaneously arise in nature; they are engineered. We should expect the person who implemented the `pow` routine analyzed the function and designed accordingly. The Apple `pow` routine, for example, has results accurate within 1 ULP throughout its domain. – Eric Postpischil Nov 12 '20 at 13:07
  • 2
    There is no general answer to this question as posed, because the accuracy depends on the `pow` implementation. You should add the information about the specific `pow` version to the question (not just in a comment). Then somebody might investigate the source code and provide an answer. – Eric Postpischil Nov 12 '20 at 13:09
  • @EricPostpischil In the upper regions of log, the derivative is close to zero. This means that log(x) and log(x+d) will be *the same value* (in IEEE floating point) for small d, even if x and x+d are not. It's a floating point representation issue. So log(x+d) is *losing information* in that range. It has to, no matter how perfect the log implementation is. If you then invert it with exp, you will not get the same result. It's a subtle point that most people are unaware of. – Tom Karzes Nov 13 '20 at 00:27
  • @TomKarzes: `pow` is not required to use the actual `log` and `exp` routines or to limit itself to the hardware floating-point format, and quality implementations do not do so. Internally, a good `pow` implementation will compute whatever logarithm it needs to whatever precision and accuracy it needs. – Eric Postpischil Nov 13 '20 at 00:40
  • @EricPostpischil Of course, and I explicitly addressed that in my original comment, in which I said "The `pow` function could presumably perform a similar magnitude shift internally, but there's no guarantee that it will go to such effort." So yes, a highly accurate `pow` function will avoid this problem. The sweet spot is where `log` and `exp` are both operating with a derivative >= 1, in which case neither will lose accuracy. – Tom Karzes Nov 13 '20 at 00:53
  • @EricPostpischil I think I overstated the problem. It's not just the derivative that affects the accuracy, but also the magnitude. The derivative times delta `x` gives the absolute difference, but with a floating point representation, it's relative differences that matter, i.e. the difference divided by the magnitude. Since the magnitude of `log(x)` grows more slowly than `x`, this helps reduce the inaccuracy. For `log(x)`, it's not `1/x` that matters, but `1/log(x)`. – Tom Karzes Nov 13 '20 at 02:59
  • @TomKarzes: `pow` routines are engineered, not slapped together haphazardly. – Eric Postpischil Nov 13 '20 at 03:09
  • @EricPostpischil Yes, the results of the glibc `pow` function appear to be much more accurate than just evaluating `exp(y * log(x))`, enough so that I would not attempt to second-guess it and would just use it directly. – Tom Karzes Nov 13 '20 at 03:46
  • @slingeraap The comment *might* have been applicable to the *source* platform where this legacy code was first used. Do you know what the original system was on which this code was deployed? Was the code using floating-point computation when first written? – njuffa Nov 13 '20 at 10:00

0 Answers0