Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Very slow std::pow() for bases very close to 1

I have a numerical code that solves an equation f(x) = 0, in which I have to raise x to a power p. I solve it using a bunch of things, but in the end I have Newton's method. The solution happens to be equal to x = 1, and thus is the cause of my problems. When the iterated solution gets close to 1, say x = 1 + 1e-13, the time necessary to compute std::pow(x, p) grows tremendously, easily by a factor of 100, making my code unusable.

The machine running this thing is AMD64 (Opteron 6172) on CentOS, and the command is simple y = std::pow(x, p);. Similar behavior shows up on all my machines, all x64. As documented here, this is not only my problem (i.e., someone else is pissed too), appears only on x64 and only for x close to 1.0. Similar thing happens to exp.

Solving this problem is vital to me. Does anyone know if there is any way of going around this slowness?

EDIT: John pointed out that this is due to denormals. The question is then, how to fix this? The code is C++, compiled with g++ for use within GNU Octave. It appears that, although I have set CXXFLAGS to include -mtune=native and -ffast-math, that is not helping and code runs just as slowly.

PSEUDO-SOLUTION FOR NOW: To all who care about this problem, the solutions suggested below didn't work for me personally. I really need the usual speed of std::pow(), but without the sluggishness around x = 1. The solution for me personally is to use the following hack:

inline double mpow(double x, double p) __attribute__ ((const));

inline double mpow(double x, double p)
{
    double y(x - 1.0);
    return (std::abs(y) > 1e-4) ? (std::pow(x, p)) : (1.0 + p * y * (1.0 + (p - 1.0) * y * (0.5 + (1.0 / 6.0) * (p - 2.0) * y)));
}

The bound could be changed, but for -40 < p < 40 the error is smaller than about 1e-11, which is good enough. Overhead is minimal from what I found, thus that solves the issue for me.

like image 632
fledgling Cxx user Avatar asked Feb 04 '13 13:02

fledgling Cxx user


People also ask

Is std :: pow slow?

If you are using double precision ( double ), std::pow(x, n) will be slower than the handcrafted equivalent unless you use -ffast-math, in which case, there is absolutely no overhead. The overhead without using the compiler option is quite large, around 2 orders of magnitude, starting from the third power.

Is C++ POW fast?

YES! Very fast if you only need 'y'/'n' as a long/int which allows you to avoid the slow FPU FSCALE function.

Is POW faster than multiplication C++?

C++11 Performance Tip: Update on When to Use std::pow When not compiling with -ffast-math, direct multiplication was significantly faster than std::pow , around two orders of magnitude faster when comparing x * x * x and code:std::pow(x, 3) .

Is POW function accurate?

ANSWER. There is no problem with the pow function. The problem is that the pow function accepts and returns floating-point numbers (which have a maximum precision of 7 decimal digits). The exp and log functions are used by pow for its calculations.


2 Answers

The obvious workaround is to note that in the reals, a ** b == exp(log(a) * b) and use that form instead. You'll need to check that it doesn't adversely affect the accuracy of your results. Edit: as discussed, this also suffers from the slowdown to almost as great a degree.

The problem isn't denormals, at least not directly; trying to compute exp(-2.4980018054066093e-15) suffers the same slowdown, and -2.4980018054066093e-15 is certainly not denormal.

If you don't care about the accuracy of your results, then scaling either the exponend or the exponent should get you outside the slow zone:

sqrt(pow(a, b * 2))
pow(a * 2, b) / pow(2, b)
...

This bug is known to glibc maintainers: http://sourceware.org/bugzilla/show_bug.cgi?id=13932 - if you are looking for a fix as opposed to a workaround you'd want to commission a floating-point math expert with open source experience.

like image 141
ecatmur Avatar answered Sep 26 '22 19:09

ecatmur


64-bit Linux?

Use the pow() code from FreeBSD.

Linux C library (glibc) has horrible worst-case performance for certain inputs.

See: http://entropymine.com/imageworsener/slowpow/

like image 30
Kirn Gill Avatar answered Sep 25 '22 19:09

Kirn Gill