Implementace funkce sinus jako ochutnávka libm

Sinus jako C funkce

Součástí specifikace programovacího jazyka C je i popis standardní knihovny funkcí. Kromě každodenních funkcí pro práci se stringy, soubory, atd je součástí specifikace i práce s reálnými čísly (IEEE float) a sada několika základních reálných funkcí, typicky poskytovaných v knihovně libm.

První specifikace C89 připouští jen argumenty a hodnoty typu double, pozdější specifikace C99 zavádí i verze pro float a long double. Například pro funkci sinus tak existují sin(3), sinf(3) a sinl(3). Pozdější standardy už v matematické knihovně nic podstatného nemění.

Zatímco pro práci se strojovými čísly jsou specifikovány i detaily aritmetických operací, způsob výpočtu dalších matematických funkcí specifikován není; záleží tedy na implementaci.

Jako příklad prozkoumáme implementaci funkce sinus v systému OpenBSD, pro jednoduchost sinf(), tedy single-precision verzi. Ve FreeBSD, NetBSD a několika dalších systémech včetně Androidu je implementace velmi podobná, ne-li stejná: všechny vycházejí z matematické knihovny někdejšího operačního systému Solaris od Sun Microsystems. Na sinusoidě se od roku 1993 pravda mnoho nezměnilo, zato na procesorech a jejich práci s reálnými čísly ano.

Naším záměrem je ukázat aplikace prostředků matematické analýzy při práci s reálnou funkcí, především Taylorův polynom, odhad jeho chyby, a pojem tečny založený na derivaci. Nevyhneme se ale několika vysloveně technickým odbočkám k počítání se strojovými čísly, které má předem omezenou přesnost a je tedy samo zatíženo chybou.

Aproximace

Sinus je spojitá reálná funkce, takže žádná její strojová implementace z podstaty nemůže být "správně": místo reálných čísel máme v počítači k dispozici jen konečnou diskrétní podmnožinu racionálních čísel, totiž strojových čísel (float). Nabízí se otázka, co se tedy "opravdu" počítá. Zdrojů chyb je několik:

• použitá reprezentace reálných čísel má omezenou přesnost
• tato přesnost nemá všude na reálné přímce stejný absolutní dopad
• místo sinusoidy počítáme její aproximaci Taylorovým polynomem
• ani to ne, Taylorův polynom aproximujeme jeho tečnou

Připomeňme obecný tvar Taylorova polynomu. Má-li reálná funkce $f$ na nějakém okolí bodu $a$ konečnou $(n+1)$-ní derivaci, potom její Taylorův polynom stupně $n$ v bodě $a$ je $$ T_n(x) = \sum_{k=0}^n \frac{f^{(k)}(a)}{k!}(x-a)^k, $$ kde $f^{(k)}$ je $k$-tá derivace funkce $f$, speciálně $f^{(0)}$ je funkce $f$ sama. Krajním případem je Taylorův polynom stupně jedna, což je právě rovnice tečny v daném bodě: $$ t(x) = f(a) + f'(a) \cdot (x-a). $$

V případě funkce sinus není s existencí derivace problém: má konečnou derivaci všech řádů všude na reálné přímce (je to hladká funkce). Taylorův rozvoj funkce sinus na okolí nuly je potom $$ \sin(0) + \frac{\cos(0)}{1}x + \frac{-\sin(0)}{2!}x^2 + \frac{-\cos(0)}{3!}x^3 + \cdots $$ neboli $$ x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \cdots $$ Všechny sudé členy díky $\sin(0)$ vypadnou, sinus je lichá funkce.

Podle Taylorovy věty je chyba takové polynomiální aproximace v každém bodě $x$ z uvažovaného okolí rovna $$ R_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}(x-a)^{n+1}, $$ pro nějaké číslo $\xi$ ležící mezi $a$ a $x$.

Implementace, kterou popíšeme níže, používá Taylorův polynom stupně 13 pro sinus v okolí nuly, totiž pro body z intervalu $(-\frac{\pi}{4}, +\frac{\pi}{4})$. Pro každý bod $x$ z toho intervalu potom existuje nějaké $\xi < |x|$ tak, že chyba použité aproximace je podle výše uvedeného vzorce rovna $$ R_{13}(x) = \frac{-\sin(\xi)}{14!}x^{14} $$ a tedy není v absolutní hodnotě větší než $\frac{1}{14!} \cdot (\frac{\pi}{4})^{14} \approx 0.000000000000389807317$.

Implementace

Implementace samotné funkce sinf() je jednoduchá, projdeme ji celou včetně komentářů.

/* s_sinf.c -- float version of s_sin.c.
 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
 */

Jak jsme uvedli výše, původní specifikace matematické knihovny jazyka C počítala jen s čísly typu double. Toto je verze pro (single) float.

/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */

V úvodu jsme zmínili, že tato implementace vychází z původního kódu Solaris OS. Systém Solaris je nyní i oficiálně mrtvý, ale knihovna je již několik let udržována jako Freely Distributable Math Library (FDLIBM) a její různé reinkarnace. Tentýž základ používá FreeBSD, NetBSD, Android a Julia (OpenLibm).

Linuxové systémy typicky používají implementaci z GNU libc, jejíž matematická knihovna dnes vychází z IBM Accurate Mathematical Library. Touto implementací se zabývat nebudeme, případným masochistům mezi čtenáři v tom ale nic nebrání.

Na okamžik si vypůjčíme ještě komentář z (double) funkce sin(), který je ve float verzi smazán, ale osvětluje, co se bude dít:

/* sin(x)
* Return sine function of x.
*
* kernel function:
*      __kernel_sin            ... sine function on [-pi/4,pi/4]
*      __kernel_cos            ... cose function on [-pi/4,pi/4]
*      __ieee754_rem_pio2      ... argument reduction routine
*
* Method.
*      Let S,C and T denote the sin, cos and tan respectively on
*      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
*      in [-pi/4 , +pi/4], and let n = k mod 4.
*      We have
*
*          n        sin(x)      cos(x)        tan(x)
*     ----------------------------------------------------------
*          0          S           C             T
*          1          C          -S            -1/T
*          2         -S          -C             T
*          3         -C           S            -1/T
*     ----------------------------------------------------------
*
* Special cases:
*      Let trig be any of sin, cos, or tan.
*      trig(+-INF)  is NaN, with signals;
*      trig(NaN)    is that NaN;
*
* Accuracy:
*      TRIG(x) returns trig(x) nearly rounded
*/

Předně se využije toho, že sinus je lichá periodická funkce, a argument se redukuje do intervalu $\left(-\frac{\pi}{4}, +\frac{\pi}{4}\right)$. To je naše okolí bodu $a=0$ z Taylorovy věty. Na reálné přímce to znamená posunout dané číslo o vhodný celočíselný násobek $\frac{\pi}{2}$, ale my nejsme na reálné přímce. Různě posunuté "kopie" intervalu $\left(-\frac{\pi}{4}, +\frac{\pi}{4}\right)$ jsou ve strojové reprezentaci různě dlouhé, mají různý počet prvků, a tyto prvky jsou v nich rozmístěny různým způsobem. Redukce argumentu do $\left(-\frac{\pi}{4}, +\frac{\pi}{4}\right)$ tedy není úplně triviální, ba dokonce je poměrně složitá. Děje se pomocí neveřejné funkce \_\_ieee754_rem_pio2f(), která převede dané číslo $x$ do tvaru $x = y_0 + y_1 - n \frac{\pi}{2}$, kde $y_0 + y_1$ je z intervalu $\left(-\frac{\pi}{4}, +\frac{\pi}{4}\right)$, a vrátí příslušné $n$, neboli informaci o tom, v které posunuté kopii intervalu $\left(-\frac{\pi}{4}, +\frac{\pi}{4}\right)$ se původní číslo $x$ nachází. Čísla $y_0$ a $y_1$ se nazývají head a tail čísla $x$; číslo $y_0$ je nejlepší možný representant čísla $x$ v intervalu $\left(-\frac{\pi}{4}, +\frac{\pi}{4}\right)$ a $y_1$ opravuje chybu, které se $y_0$ dopouští, protože rozložení strojových čísel v obou intervalech je různé.

Vnitřnosti této redukční funkce zkoumat nebudeme; nepřekvapí nás ale, že její implementace stojí na předpočítaných násobcích a zlomcích čísla $\frac{\pi}{2}$ (v dané přesnosti), a na intimní znalosti rozmístění strojových čísel v reálné přímce. Níže se k ní vrátíme na několika příkladech, nyní pokračujeme ve čtení s_sinf.c.

#include "math.h"
#include "math_private.h"

float
sinf(float x)
{
    float y[2],z=0.0;
    int32_t n, ix;

    GET_FLOAT_WORD(ix,x);

    /* |x| ~< pi/4 */
    ix &= 0x7fffffff;
    if(ix <= 0x3f490fdb) return __kernel_sinf(x,z,0);

    /* sin(Inf or NaN) is NaN */
    else if (ix>=0x7f800000) return x-x;

    /* argument reduction needed */
    else {
        n = __ieee754_rem_pio2f(x,y);
        switch(n&3) {
        case 0: return  __kernel_sinf(y[0],y[1],1);
        case 1: return  __kernel_cosf(y[0],y[1]);
        case 2: return -__kernel_sinf(y[0],y[1],1);
        default:
            return -__kernel_cosf(y[0],y[1]);
        }
    }
}

Funkce sinus je periodická s periodou $2\pi$. Při redukci argumentu tedy pomocí n % 4 potřebného posunu o $n$-násobek $\frac{\pi}{2}$ zjistíme, jestli máme po redukci počítat funkci $\sin(x)$ nebo $\cos(x)$ nebo $-\sin(x)$ nebo $-\cos(x)$. (Detail: n % 4 je v binární reprezentaci totéž co n & 3.) Nejde tu o první čtyři derivace (totiž nultou až třetí derivaci) funkce $\sin(x)$, ale o posun argumentu o násobek $n \frac{\pi}{2}$ pro $n = 0,1,2,3$ a periodicky dále; stačí si prohlednout grafy těchto čtyř funkcí.

Pro argument z intervalu $\left( -\frac{\pi}{4}, +\frac{\pi}{4} \right)$ pak budeme počítat hodnotu sinu pomocí funkce \_\_kernel_sinf respektive \_\_kernel_cosf. Pro příliš velké nebo neplatné hodnoty argumentu hned končíme. Ostatní argumenty nejprve redukujeme, a pak počítáme příslušnou funkci.

Strojová čísla jako integery

První technická odbočka: test na $|x| < \frac{\pi}{4}$ se děje jinak, než bychom čekali, totiž místo porovnání čísel x < M_PI_4 se bere float word čísla $x$ (makro je definováno v math_private.h), tedy hledíme na 32 bitů, které v paměti representují hodnotu čísla $x$, jako na 32-bitové přirozené číslo (unsigned int). Zamaskováním sign bitu pomocí 0x7fffffff zároveň ignorujeme znaménko. Magická konstanta 0x3f490fdb je pak právě hodnota čísla M_PI_4, což snadno ověříme:

#include <stdio.h>
#include <math.h>

typedef union
{
    float f;
    u_int32_t u;
} ufloat;

#define GET_FLOAT_WORD(_u,_f) { \
    ufloat uf;      \
    uf.f = (_f);        \
    (_u) = uf.u;        \
}

int
main()
{
    int32_t i;
    float f = M_PI_4;
    GET_FLOAT_WORD(i, f);
    printf("float %.8e, uint %#x\n", f, i);
    return 0;
}

Podobných tanců s převodem mezi strojovým (float) číslem a jeho integerovou podobou je knihovna libm plná. Smyslem těchto obratů je (nebo alespoň v době psaní bylo) co nejvíce omezit použití floating-point aritmetiky a samotné výpočty provádět v integerové aritmetice. Ta byla na tehdejších procesorech rychlejší a přesnější; to už dnes není pravda, ale kód tak zásadní jako knihovna matematických funkcí má zjevně dlouhou setrvačnost.

Podobně test na NaN a INF ("not a number", resp. "infinity") se místo nyní standardních funkcí isnan() a isinf() děje testováním float wordu na >= 0x7f800000, což opět vychází z binární reprezentace čísel (nebo v tomto případě "ne-čísel"). Pro takové argumenty vracíme x-x, což je v obou případech NaN; přitom od specifikace C99 existuje standardní funkce nan().

Odvážný čtenář může vybrané partie knihovny libm přepsat do čisté floating-point aritmetiky a otestovat na všech platformách, které libm používají, co všechno se změnilo.

Vnitřnosti

Jádrem implementace je potom interní funkce \_\_kernel_sinf(). Vypůjčíme si opět nejdříve komentář z double verze, který popisuje, co přesně budeme počítat:

/* __kernel_sin(x, y, iy)
 * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
 * Input x is assumed to be bounded by ~pi/4 in magnitude.
 * Input y is the tail of x.
 * Input iy indicates whether y is 0. (if iy=0, assume y to be 0).
 *
 * Algorithm
 *  1. Since sin(-x) = -sin(x), we need only to consider positive x.
 *  2. if x < 2^-27 (hx<0x3e400000), return x with inexact if x!=0.
 *  3. sin(x) is approximated by a polynomial of degree 13 on
 *     [0,pi/4]
 *                   3            13
 *      sin(x) ~ x + S1*x + ... + S6*x
 *     where
 *  
 *  |sin(x)         2     4     6     8     10     12  |     -58
 *  |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
 *  |  x                               |
 *
 *  4. sin(x+y) = sin(x) + sin'(x')*y
 *          ~ sin(x) + (1-x*x/2)*y
 *     For better accuracy, let
 *           3      2      2      2      2
 *      r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
 *     then                   3    2
 *      sin(x) = x + (S1*x + (x *(r-y/2)+y))
 */

Počítáme již s argumentem redukovaným do intervalu $\left( -\frac{\pi}{4}, +\frac{\pi}{4} \right)$, navíc bez újmy na obecnosti kladným. Pro argument dostatečně blízký nule rovnou vracíme $x$ jako hodnotu $\sin(x)$; to je nepřesné, s výjimkou případu $\sin(0) = 0$, ale přijatelné, díky známé limitě $$ \lim_{x \to 0}{\frac{\sin(x)}{x}} = 1 . $$

Pro ostatní argumenty se konečně bude počítat Taylorův polynom, a to stupně 13. K výpočtu takového polynomu potřebujeme znát sedm koeficientů stojících u lichých mocnin $x^1, x^3, \dots, x^{13}$ (viz výše). První z nich je roven jedné, ostatní jsou předpočítány:

#include "math.h"
#include "math_private.h"

static const float
half =  5.0000000000e-01,/* 0x3f000000 */
S1  = -1.6666667163e-01, /* 0xbe2aaaab */
S2  =  8.3333337680e-03, /* 0x3c088889 */
S3  = -1.9841270114e-04, /* 0xb9500d01 */
S4  =  2.7557314297e-06, /* 0x3638ef1b */
S5  = -2.5050759689e-08, /* 0xb2d72f34 */
S6  =  1.5896910177e-10; /* 0x2f2ec9d3 */

Nabízí se otázka, proč počítáme právě polynom stupně $13$, neboť pro funkci sinus je aproximace Taylorovým polynomem na okolí nuly s každým dalším stupněm přesnější. Necháme na čtenáři, aby současnou implementaci pozměnil na Taylorův polynom stupně $15$ a změřil rozdíly. Lze očekávat, že zlepšení mezi $T_{13}$ a $T_{15}$ se na intervalu $\left( -\frac{\pi}{4}, +\frac{\pi}{4} \right)$ při single float reprezentaci už neprojeví.

float
__kernel_sinf(float x, float y, int iy)
{
    float z,r,v;
    int32_t ix;
    GET_FLOAT_WORD(ix,x);
    ix &= 0x7fffffff;           /* high word of x */
    if(ix<0x32000000)           /* |x| < 2**-27 */
       {if((int)x==0) return x;}        /* generate inexact */
    z   =  x*x;
    v   =  z*x;
    r   =  S2+z*(S3+z*(S4+z*(S5+z*S6)));
    if(iy==0) return x+v*(S1+z*r);
    else      return x-((z*(half*y-v*r)-y)-v*S1);
}

Pro zmatení nepřítele se nyní head a tail redukovaného argumentu jmenují $x$ a $y$. Pomocí Taylorova polynomu spočteme přibližnou hodnotu $\sin(x)$, a tuto případně ještě upravíme: je-li y != 0 (tato informace se z nějakého důvodu předává indikátorem iy), má být "skutečný" argument v intervalu okolo nuly o kousek jinde, totiž v bodě $x+y$. Hodnotu v tomto posunutém bodě aproximujeme tečnou; přitom derivací funkce $\sin$ je funkce $\cos$, jejíž hodnotu (tj. směrnici tečny) opět aproximujeme, totiž jejím Taylorovým polynomem $1 - \frac{x^2}{2!}$ stupně dva. Všimněme si také, jakým výrazem se hodnota polynomu vyčísluje. Místo naivního násobení až do příslušné mocniny a následného sčítání takových členů se inteligentní úpravou minimalizuje počet potřebných sčítání a násobení.

Redukce

Druhá technická odbočka: vrátíme se ještě na chvíli k redukci předloženého argumentu do intervalu $\left( -\frac{\pi}{4}, +\frac{\pi}{4} \right)$. Nebudeme zkoumat samotnou redukční proceduru, podíváme se jen na výsledky pro několik význačných hodnot. K tomu používáme následující kód:

$ cat reduction.c
#include <stdio.h>
#include <math.h>
#include "math_private.h"

struct point {
    const char *name;
    float val;
} points[] = {
    { "zero",        0.0              },
    { "pi/8",        0.5 * M_PI_4     },
    { "pi/4",        M_PI_4           },
    { "pi/4",        M_PI_4           },
    { "pi/2",        M_PI_2           },
    { "3pi/4",       3.0 * M_PI_4     },
    { "pi",          M_PI             },
    { "3pi/2",       3.0 * M_PI_2     },
    { "2pi",         2.0 * M_PI       },
    { "3pi",         3.0 * M_PI       },
    { "4pi",         4.0 * M_PI       },
    { "10^6 pi/4",   1.0e6 * M_PI_4   },
    { "10^6 pi/2",   1.0e6 * M_PI_2   },
    { "10^7 pi/4",   1.0e7 * M_PI_4   },
    { "10^7 pi/2",   1.0e7 * M_PI_2   },
    { "10^100 pi/4", 1.0e100 * M_PI_4 },
    { "10^100 pi/2", 1.0e100 * M_PI_2 },
    { "10^120 pi/4", 1.0e120 * M_PI_4 },
    { "10^120 pi/2", 1.0e120 * M_PI_2 },
    { NULL, 0 }
};

void
print(struct point *p)
{
    float y[2];
    int32_t k;

    k = __ieee754_rem_pio2f(p->val, y);
    printf("%-11s % .8e (k = %d), head % .8e, tail % .8e\n",
        p->name, p->val, k, y[0], y[1]);
}

int
main()
{
    struct point *p = points;
    while (p && p->name)
        print(p++);
    return 0;
}

$ cc -Wall -pedantic -c e_rem_pio2f.c
$ cc -Wall -pedantic -c k_rem_pio2f.c
$ cc -Wall -pedantic -lm -o reduction reduction.c e_rem_pio2f.o k_rem_pio2f.o

Na číslech z intervalu $\left( -\frac{\pi}{4}, +\frac{\pi}{4} \right)$ se neděje nic: head je číslo samo, tail je nula.

zero         0.00000000e+00 (k = 0), head  0.00000000e+00, tail  0.00000000e+00
pi/8         3.92699093e-01 (k = 0), head  3.92699093e-01, tail  0.00000000e+00

Hned v sousedních intervalech už se dopouštíme chyby: číslo $\frac{\pi}{2}$ hraje v intervalu $\left( \frac{\pi}{4}, \frac{3\pi}{4} \right)$ stejnou roli, jako číslo $0$ v intervalu $\left( -\frac{\pi}{4}, +\frac{\pi}{4} \right)$, přitom redukcí $\frac{\pi}{2}$ není čistá nula. Podobně pro další; například $\sin(4\pi)$ už dokonce nebude nula.

pi/4         7.85398185e-01 (k = 1), head -7.85398126e-01, tail -1.58934199e-08
pi/2         1.57079637e+00 (k = 1), head  4.37113883e-08, tail  1.72084569e-15
3pi/4        2.35619450e+00 (k = 2), head -7.85398185e-01, tail  2.78178049e-08
pi           3.14159274e+00 (k = 2), head  8.74227766e-08, tail  3.44169138e-15
3pi/2        4.71238899e+00 (k = 3), head  1.19248806e-08, tail -1.83697028e-16
2pi          6.28318548e+00 (k = 4), head  1.74845553e-07, tail  6.88338275e-15
3pi          9.42477798e+00 (k = 6), head  2.38497613e-08, tail -3.67394056e-16
4pi          1.25663710e+01 (k = 8), head  3.49691106e-07, tail  1.37667655e-14

Čím dále je argument od nuly, tím je absolutní chyba větší, jak strojová čísla postupně "řídnou" (relativní chyba zůstává podobná). Například počítání sinu v $10^6 \frac{\pi}{4}$ znamená sinus přibližně v $0.0241$:

10^6 pi/4    7.85398188e+05 (k = 0), head  2.41025519e-02, tail -1.72213133e-10
10^6 pi/2    1.57079638e+06 (k = 0), head  4.82051037e-02, tail -3.44426265e-10
10^7 pi/4    7.85398150e+06 (k = 0), head -1.33974478e-01, tail -5.44742162e-09
10^7 pi/2    1.57079630e+07 (k = 7), head  1.30284739e+00, tail -2.48039100e-08

Pro skutečně velké hodnoty už je chyba nepřijatelná. To je zřejmě smyslem následující poznámky v manuálu sin: A large magnitude argument may yield a result with little or no significance.

10^100 pi/4  inf (k = 0), head -nan, tail -nan
10^100 pi/2  inf (k = 0), head -nan, tail -nan
10^120 pi/4  inf (k = 0), head -nan, tail -nan
10^120 pi/2  inf (k = 0), head -nan, tail -nan

Vidíme tedy, že už samotná reprezentace reálných čísel strojovými čísly přináší při výpočtu funkce sinus problémy, které nesouvisí s použitou aproximací.

Závěr

Implementace sinusoidy ve standardní knihovně libm se děje pomocí Taylorova polynomu. Tato implementace je náchylná k různým druhům chyb, obzvláště pro velké argumenty; jedním z nich je samotný způsob reprezentace reálných čísel strojovými čísly v počítači.

Aproximace Taylorovým polynomem není jediná možnost, jak sinus počítat. Přesnějších výsledků (nejen pro sinus, ale i pro mnoho jiných funkcí) lze dosáhnou s použitím Čebyševových polynomů.

Některé moderní procesory mají i pro funkce jako sinus samostatnou instrukci, jako například fsin v instrukční sadě x86/x64. To umožňuje implementovat sinus v asembleru. Bohužel, použitá hardwarová implementace může být někdy zoufale nepřesná.

Případné zájemce odkazujeme na knihu Computer Approximations.

---

Tento článek byl vznikl v rámci řešení projektu Institucionální podpora Českého vysokého učení technického v Praze (CZ.02.2.69/0.0/0.0/16_015/0002382).

---