Version 2 (modified by sam, 11 years ago) (diff)

--

Remez tutorial 5/5: additional tips

This part of the tutorial is constantly updated. Check its history to know what is new.

When not to use Remez?

There are cases when you should expect the Remez algorithm to potentially perform badly:

There are cases where you should not try to use the Remez algorithm at all:

  • when the range is not finite, eg. [0,+∞]
  • when the function to approximate has an infinite derivative at a point contained in or near the approximation range, eg. sqrt(x) or cbrt(x) in 0

What if I want to use Remez anyway?

If you need to approximate a function f(x) over [a,+∞] and for some reason you want to use the Remez exchange algorithm, you can still through a change of variable: y = 1/x. The function to approximate becomes f(1/y) and the new range is [0,1/a] (see “changing variables” for how to deal with 1/x in 0. Your minimax polynomial will use 1/x as its variable; please be aware that computing 1/x at runtime may be expensive.

When you want arbitrary runtime precision

There are In the previous section, we saw how to fix lower-order parameters took advantage of the symmetry of sin(x) to build the following minimax expression:

\[\max_{y \in [0, \pi^2/4]}{\dfrac{\bigg\lvert\dfrac{\sin(\sqrt{y})}{\sqrt{y}} - Q(y)\bigg\rvert}{\dfrac{1}{|\sqrt{y}|}}} = E\]

Leading to the following coefficient, amongst others:

const double a1 = 9.999999766327004083636926544075452569800e-1;

This is an interesting value, because it is very close to 1. Many CPUs can load the value 1 very quickly, which can be a potential runtime gain.

The brutal way

Now we may wonder: what would be the cost of directly setting a1 = 1 here? Let’s see the error value:

No image "bad-error.png" attached to doc/maths/remez/tutorial-additional-tips

Duh. Pretty bad, actually. Maximum error is about 10 times worse.

The clever way

The clever way involves some more maths. Instead of looking for polynomial Q(y) and setting Q(0) = 1 manually, we search instead for R(y) where Q(y) = 1 + yR(y):

\[\max_{y \in [0, \pi^2/4]}{\dfrac{\bigg\lvert\dfrac{\sin(\sqrt{y})}{\sqrt{y}} - 1 - yR(y)\bigg\rvert}{\dfrac{1}{|\sqrt{y}|}}} = E\]

Dividing by y up and down gives:

\[\max_{y \in [0, \pi^2/4]}{\dfrac{\bigg\lvert\dfrac{\sin(\sqrt{y})-\sqrt{y}}{y\sqrt{y}} - R(y)\bigg\rvert}{\dfrac{1}{|y\sqrt{y}|}}} = E\]

Once again, we get a form suitable for the Remez algorithm.

Source code

#include "lol/math/real.h"
#include "lol/math/remez.h"

using lol::real;
using lol::RemezSolver;

real f(real const &y)
{
    real sqrty = sqrt(y);
    return (sin(sqrty) - sqrty) / (y * sqrty);
}

real g(real const &y)
{
    return re(y * sqrt(y));
}

int main(int argc, char **argv)
{
    RemezSolver<3, real> solver;
    solver.Run("1e-1000", real::R_PI_2 * real::R_PI_2, f, g, 30);
    return 0;
}

Only f and g changed here, as well as the polynomial degree. The rest is the same as in the previous section.

Compilation and execution

Build and run the above code:

make
./remez

After all the iterations the output should be as follows:

Final error: 4.613764035583927492630791542066631050746e-9
Polynomial estimate:
x**0*-1.666665708792613295297026208726412532085e-1
+x**1*8.333017134192478312472752663154642556843e-3
+x**2*-1.980660618990173078599247774582918018345e-4
+x**3*2.600038299369332084157603784882057661080e-6

We can therefore write the corresponding C++ function:

double fastsin2(double x)
{
    const double a3 = -1.666665708792613295297026208726412532085e-1
    const double a5 = 8.333017134192478312472752663154642556843e-3
    const double a7 = -1.980660618990173078599247774582918018345e-4
    const double a9 = 2.600038299369332084157603784882057661080e-6

    return x + x*x*x * (a3 + x*x * (a5 + x*x * (a7 + x*x * a9))));
}

Note that because of our change of variables, the polynomial coefficients are now a3, a5, a7

Analysing the results

Let’s see the new error curve:

No image "better-error.png" attached to doc/maths/remez/tutorial-additional-tips

Excellent! The loss of precision is clearly not as bad as before.

Conclusion

Please report any trouble you may have had with this document to sam@hocevar.net. You may then return to the Remez documentation.