Changes between Version 1 and Version 2 of research/trig


Ignore:
Timestamp:
Oct 10, 2011, 1:35:03 AM (9 years ago)
Author:
sam
Comment:

minimax for sin(x)

Legend:

Unmodified
Added
Removed
Modified
  • research/trig

    v1 v2  
    22
    33Our research notes about implementation of fast trigonometric functions.
     4
     5== Minimax polynomial for sin(x) ==
     6
     7=== Absolute error ===
     8
     9Suppose we want to approximate sin(x) on [-π/2; π/2] with a polynomial P(x) such that the '''absolute''' error is never more than E:
     10
     11{{{
     12#!latex
     13\[\big\vert\sin(x) - P(x)\big\vert \le E \qquad \forall x \in \bigg[-\frac{\pi}{2}, \frac{\pi}{2}\bigg]\]
     14}}}
     15
     16We know sin(x) is an odd function, so instead we look for a polynomial Q(x) such that P(x) = xQ(x²):
     17
     18{{{
     19#!latex
     20\[\big\vert\sin(x) - xQ(x^2)\big\vert \le E \qquad \forall x \in \bigg[-\frac{\pi}{2}, \frac{\pi}{2}\bigg]\]
     21}}}
     22
     23Substitute y for x² and reduce the range to positive values:
     24
     25{{{
     26#!latex
     27\[\big\lvert\sin(\sqrt{y}) - \sqrt{y}Q(y)\big\rvert \le E \qquad \forall y \in \bigg[0, \frac{\pi^2}{4}\bigg]\]
     28}}}
     29
     30Divide through by √y:
     31
     32{{{
     33#!latex
     34\[\bigg\lvert\frac{\sin(\sqrt{y})}{\sqrt{y}} - Q(y)\bigg\rvert \le \frac{E}{|\sqrt{y}|} \qquad \forall y \in \bigg[0, \frac{\pi^2}{4}\bigg]\]
     35}}}
     36
     37If we want to force the asymptotic behaviour in x=0, we substitute Q(y) with 1+yR(y):
     38
     39{{{
     40#!latex
     41\[\bigg\lvert\frac{\sin(\sqrt{y})}{\sqrt{y}} - 1 - yR(y)\bigg\rvert \le \frac{E}{|\sqrt{y}|} \qquad \forall y \in \bigg[0, \frac{\pi^2}{4}\bigg]\]
     42}}}
     43
     44Divide through by y:
     45
     46{{{
     47#!latex
     48\[\bigg\lvert\frac{\sin(\sqrt{y})-\sqrt{y}}{y\sqrt{y}} - R(y)\bigg\rvert \le \frac{E}{|y\sqrt{y}|} \qquad \forall y \in \bigg[0, \frac{\pi^2}{4}\bigg]\]
     49}}}
     50
     51We then use the following code:
     52
     53{{{
     54#!cpp
     55
     56static real myfun(real const &x)
     57{
     58    real y = sqrt(x);
     59    return (sin(y) - y) / (x * y);
     60}
     61
     62static real myerr(real const &x)
     63{
     64    real y = sqrt(x);
     65    return re(x * y);
     66}
     67
     68RemezSolver<6> solver;
     69solver.Run(real::R_1 >> 400, real::R_PI_2 * real::R_PI_2, myfun, myerr, 15);
     70}}}
     71
     72These are the resulting R, Q and P:
     73
     74{{{
     75#!latex
     76\begin{eqnarray*}
     77R(y) & = & a_0 + a_1 y + a_2 y^2 + a_3 y^3 + a_4 y^4 + a_5 y^5 + a_6 y^6 \\
     78Q(y) & = & 1 + a_0 y + a_1 y^2 + a_2 y^3 + a_3 y^4 + a_4 y^5 + a_5 y^6 + a_6 y^7 \\
     79P(x) & = & x + a_0 x^3 + a_1 x^5 + a_2 x^7 + a_3 x^9 + a_4 x^{11} + a_5 x^{13} + a_6 x^{15} \\
     80a_0 & = & -1.66666666666658080941942898789420724.10^{-1} \\
     81a_1 & = & +8.33333333326271609442503773834687308.10^{-3} \\
     82a_2 & = & -1.98412698200591143928364634696492885.10^{-4} \\
     83a_3 & = & +2.75573160733868922065738227278330896.10^{-6} \\
     84a_4 & = & -2.50518513021429359590028300127165228.10^{-8} \\
     85a_5 & = & +1.60472959182597740337401201006549498.10^{-10} \\
     86a_6 & = & -7.36458957326227991327065122848667046.10^{-13} \\
     87\end{eqnarray*}
     88}}}