[[TOC]] = Trigonometric functions = Our research notes about implementation of fast trigonometric functions. == Other implementations == * OpenBSD [http://www.openbsd.org/cgi-bin/cvsweb/src/lib/libm/src/k_sin.c?rev=1.3;content-type=text%2Fplain sin(double)]. * OpenBSD [http://www.openbsd.org/cgi-bin/cvsweb/src/lib/libm/src/k_sinf.c?rev=1.3;content-type=text%2Fplain sinf(float)]. == Evaluating minimax polynomials == === Relative error on float ULPs === Let’s check how much an error on the last bit affects float variables, using relative error: {{{ #!cpp for (float f = 1e-25f; f < 50.0f; f *= 1.001f) { union { float f; uint32_t x; } u = { f }, v = u, w = u; ++v.x; --w.x; printf("%e: %e %e\n", u.f, (v.f - u.f) / u.f, (u.f - w.f) / u.f); } }}} Plotting the results gives the answer: between 5.96e-8 and 1.19e-7. === Relative error on double precision ULPs === How about double? {{{ #!cpp for (double d = 1e-25; d < 50.0; d *= 1.001) { union { double d; uint64_t x; } u = { d }, v = u, w = u; ++v.x; --w.x; printf("%e: %e %e\n", u.d, (v.d - u.d) / u.d, (u.d - w.d) / u.d); } }}} Result: between 1.11e-16 and 2.22e-16. == Minimax polynomial for sin(x) == === Absolute error === Suppose we want to approximate sin(x) on [-π/2; π/2] with a polynomial P(x) such that the '''absolute''' error is never more than E: {{{ #!latex $\max_{x \in [-\pi/2, \pi/2]}{\big\vert\sin(x) - P(x)\big\vert} = E$ }}} We know sin(x) is an odd function, so instead we look for a polynomial Q(x) such that P(x) = xQ(x²), and we reduce the range to positive values: {{{ #!latex $\max_{x \in [0, \pi/2]}{\big\vert\sin(x) - xQ(x^2)\big\vert} = E$ }}} Substitute y for x²: {{{ #!latex $\max_{y \in [0, \pi^2/4]}{\big\lvert\sin(\sqrt{y}) - \sqrt{y}Q(y)\big\rvert} = E$ }}} Divide through by √y: {{{ #!latex $\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$ }}} If we want to force the asymptotic behaviour at x=0, we substitute Q(y) with 1+yR(y): {{{ #!latex $\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$ }}} Divide through by y: {{{ #!latex $\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$ }}} We then use the following code: {{{ #!cpp static real myfun(real const &y) { real x = sqrt(y); return (sin(x) - x) / (x * y); } static real myerr(real const &y) { return re(y * sqrt(y)); } RemezSolver<6> solver; solver.Run(real::R_1 >> 400, real::R_PI_2 * real::R_PI_2, myfun, myerr, 40); }}} These are the resulting R, Q and P: {{{ #!latex \begin{eqnarray*} R(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 \\ Q(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 \\ P(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} \\ \end{eqnarray*} }}} With the following coefficients: {{{ #!cpp a0 = -1.666666666666580938362041558393413542600e-1; a1 = +8.333333333262715528278347301093116699226e-3; a2 = -1.984126982005911547055378498482154233331e-4; a3 = +2.755731607338689059680115311170244593349e-6; a4 = -2.505185130214293461336864464029272975945e-8; a5 = +1.604729591825977276746033955401272849354e-10; a6 = -7.364589573262279656101943192883347174447e-13; E = 1.098969630370672683831702893969063712485e-16; }}} === Relative error === Searching for '''relative error''' instead: {{{ #!latex $\max_{x \in [-\pi/2, \pi/2]}{\dfrac{\big\vert\sin(x) - P(x)\big\vert}{|\sin(x)|}} = E$ }}} Using the same method as for absolute error, we get: {{{ #!latex $\max_{y \in [0, \pi^2/4]}{\dfrac{\bigg\lvert\dfrac{\sin(\sqrt{y})-\sqrt{y}}{y\sqrt{y}} - R(y)\bigg\rvert}{\bigg\lvert\dfrac{\sin(y)}{y\sqrt{y}}\bigg\rvert}} = E$ }}} {{{ #!cpp static real myfun(real const &y) { real x = sqrt(y); return (sin(x) - x) / (x * y); } static real myerr(real const &y) { real x = sqrt(y); return sin(x) / (x * y); } RemezSolver<6> solver; solver.Run(real::R_1 >> 400, real::R_PI_2 * real::R_PI_2, myfun, myerr, 40); }}} {{{ #!cpp a0 = -1.666666666666618655330686951220600109093e-1; a1 = +8.333333333285509269011390197274270193963e-3; a2 = -1.984126982504390943378441670654599796999e-4; a3 = +2.755731659890484005079622148869365590255e-6; a4 = -2.505188017067512158000464673183918592305e-8; a5 = +1.604809231079007402834515514014751626424e-10; a6 = -7.373308642081174610234470417310065337523e-13; E = 1.536616934979294924070319645582179175464e-16; }}} === Another example === Trying to reproduce [http://www.insomniacgames.com/remez-exchange-algorithm/ Mike Acton]'s experimental result: {{{ #!latex $\max_{x \in [-1, 1]}{\dfrac{\big\vert\sin(\frac\pi{2}x) - P(x)\big\vert}{|\sin(\frac\pi{2}x)|}} = E$ }}} Using P(x) = xQ(x²): {{{ #!latex $\max_{y \in [0, 1]}{\dfrac{\bigg\lvert\dfrac{\sin(\frac\pi{2}\sqrt{y})}{\sqrt{y}} - Q(y)\bigg\rvert}{\bigg\lvert\dfrac{\sin(\frac\pi{2}\sqrt{y})}{\sqrt{y}}\bigg\rvert}} = E$ }}} {{{ #!cpp static real myfun(real const &y) { real x = sqrt(y); return sin(real::R_PI_2 * x) / x; } RemezSolver<4> solver; solver.Run(real::R_1 >> 400, real::R_1, myfun, myfun, 40); }}} {{{ #!cpp a0 = +1.570796318452974170937444182514099057668; a1 = -6.459637106518004316895285560913086847441e-1; a2 = +7.968967893119537554369879233657335437244e-2; a3 = -4.673766402124326801818077428151491285886e-3; a4 = +1.514849803879821510688050931907848092814e-4; E = 5.310632770140865101487747414605249755889e-9; }}}