Version 1 (modified by 11 years ago) (diff) | ,
---|

#### Table of Contents

# 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:

- when the function is not continuous (for instance, a step function)
- when the function is not differentiable (for instance, the
*x+abs(x)*function) - sometimes when the function is not smooth
- sometimes even when the function is not analytic

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:

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:

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)*:

Dividing by *y* up and down gives:

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:

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

## Conclusion

You should now be able to fix lower-order coefficients in the minimax polynomial for possible performance improvements.

Please report any trouble you may have had with this document to sam@hocevar.net. You may then carry on to the next section: additional tips.