inline void XMScalarSinCos(float* pSin, float* pCos, float Value) { assert(pSin); assert(pCos); // Map Value to y in [-pi,pi], x = 2*pi*quotient + remainder. float quotient = XM_1DIV2PI * Value; if (Value >= 0.0f) { quotient = static_cast<float>(static_cast<int>(quotient + 0.5f)); } else { quotient = static_cast<float>(static_cast<int>(quotient - 0.5f)); } float y = Value - XM_2PI * quotient; // Map y to [-pi/2,pi/2] with sin(y) = sin(Value). float sign; if (y > XM_PIDIV2) { y = XM_PI - y; sign = -1.0f; } else if (y < -XM_PIDIV2) { y = -XM_PI - y; sign = -1.0f; } else { sign = +1.0f; } float y2 = y * y; // 11-degree minimax approximation *pSin = (((((-2.3889859e-08f * y2 + 2.7525562e-06f) * y2 - 0.00019840874f) * y2 + 0.0083333310f) * y2 - 0.16666667f) * y2 + 1.0f) * y; // 10-degree minimax approximation float p = ((((-2.6051615e-07f * y2 + 2.4760495e-05f) * y2 - 0.0013888378f) * y2 + 0.041666638f) * y2 - 0.5f) * y2 + 1.0f; *pCos = sign * p; }
La prima parte è abbastanza semplice: la variabile $Value$ (l'angolo passato alla funzione) viene mappata alla variabile $y$ prima nell'intervallo $[-\pi,\pi]$ e poi in $[-\pi /2,\pi /2]$. Tenendo presente che le applicazioni 3D o 2D lavorano con angoli in intervalli $[0,2\pi]$ (per far ruotare un oggetto fino a 360° in un senso) o $[-2\pi,0]$ (per far ruotare oggetti fino a 360° nell'altro senso), si procede in questo modo distinguendo i due casi:
$(1)$ $Value=[0,2\pi]$
$(2)$ $Value=[-2\pi,0] $
$(1)$ $V=Value/2=[0,1] $
$(2)$ $V=Value/2=[-1,0] $
$(1)$ $t_{1}=(V+1/2)=[1/2,3/2]$
$(2)$ $t_{1}=(V-1/2)=[-3/2,-1/2]$
$(1)$ $t_{2}=(t_1*2\pi)=[\pi,3\pi]=[\pi,-\pi]$
$(2)$ $t_{2}=(t_1*2\pi)=[-3\pi,-\pi]=[\pi,-\pi]$
$(1)$ $y=(Value-t_2)=[0,2\pi]-[\pi,-\pi]=[-\pi,\pi]$
$(2)$ $y=(Value-t_2)=[-2\pi,0]-[\pi,-\pi]=[-\pi,\pi]$
Successivamente mappa $y$ da $[-\pi,\pi]$ a $[-\pi/2,\pi/2]$ verificando semplicemente se $\,\pi/2< y <-\pi/2$. In caso affermativo si esegue semplicemente la traslazione ($y=\pm \pi - y$) per portare l'angolo nel primo o nel quarto quadrante del sistema cartesiano in modo che ricada proprio nell'intervallo desiderato. Se questa traslazione avviene il segno del coseno cambia e deve essere ripristinato. Ed ora la parte più criptica:
float y2 = y * y; // 11-degree minimax approximation *pSin = (((((-2.3889859e-08f * y2 + 2.7525562e-06f) * y2 - 0.00019840874f) * y2 + 0.0083333310f) * y2 - 0.16666667f) * y2 + 1.0f) * y; // 10-degree minimax approximation float p = ((((-2.6051615e-07f * y2 + 2.4760495e-05f) * y2 - 0.0013888378f) * y2 + 0.041666638f) * y2 - 0.5f) * y2 + 1.0f; *pCos = sign * p;
Il codice non è altro che l'implementazione dell'approssimazione polinomiale minimax delle funzioni seno e coseno. Si può spiegare la teoria dietro questa implementazione, senza entrare troppo nei dettagli, in questo modo:
sia $f$ la funzione da approssimare (che non deve essere per forza la funzione seno o coseno) e $p'\in P_n$ il polinomio di grado $n$ che meglio approssima $f$ su $[a,b]$ ($P_n$ insieme dei polinomi di grado $n$). Per trovare $p'$ bisogna definire prima cosa si intende per migliore approssimazione di $f$. Si definisca dunque il massimo errore assoluto tra la funzione $f$ ed il polinomio $p$ come:
$||f-p||_{\infty,[a,b]}=\max\limits_{a \leq x \leq b} |f(x)-p(x)|$
che può essere intesa come la distanza massima tra $f$ e $p$. Per migliore approssimazione intendiamo il polinomio $p' \in P_n$ tale che:
$||f-p'||_{\infty,[a,b]}=\min ||f-p||_{\infty,[a,b]}$
e cioè il polinomio che minimizza il massimo errore tra $f$ e $p \in P_n$. E' per questo motivo che si chiama approssimazione polinomiale minimax. Tale polinomio esiste (teorema di approssimazione di Weierstrass, si veda [3]), è unico (si veda [5]) ed ha una particolare caratteristica che permette di calcolarlo. Si può dimostrare che se $p'$ è il polinomio di grado $n$ che meglio approssima $f$ su $[a,b]$ (secondo la definizione fornita sopra), allora esistono $n+2$ valori $x_i$
$a\leq x_0\le x_1\le x_2\le \cdots \le x_{n+1} \leq b$
tali che
$p'(x_i)-f(x_i)=(-1)^i\,[p'(x_0)-f(x_0)] = \pm ||f-p'|||_{\infty,[a,b]}$
In altre parole, esistono $n+2$ valori in cui $p'-f$ assume un valore uguale alla distanza (errore) massima tra $f$ e $p'$ con il solo segno che si alterna (teorema di equioscillazione di Chebyshev, si veda [5]). Alla luce di quanto appena detto, si può trovare $p'$ in questo modo (algoritmo di Remez, si veda [6]):
1. Si parte da un insieme di $n+2$ valori $x_0, x_1,x_2, \cdots , x_{n+1}$ arbitrari in $[a,b]$. A tale scopo, è consigliabile usare la seguente formula:
$x_i = \displaystyle\frac{a+b}{2}+\frac{b-a}{2}\cos\bigg(\frac{i\pi}{n+1}\bigg)$
2. Si imposta il seguente sistema di $n+2$ equazioni in $n+2$ variabili:
$\begin{cases} p_0 + p_1x_0 + p_2x_0^2 + \cdots + p_nx_0^n - f(x_0) = +\epsilon
\\ p_0 + p_1x_1 + p_2x_1^2 + \cdots + p_nx_1^n - f(x_1) = -\epsilon
\\ p_0 + p_1x_2 + p_2x_2^2 + \cdots + p_nx_2^n - f(x_2) = +\epsilon
\\ \cdots
\\ p_0 + p_1x_{n+1} + p_2x_{n+1}^2 + \cdots + p_nx_{n+1}^n - f(x_{n+1}) = (-1)^{n+1}\epsilon
\end{cases}$
Tale sistema riflette quella che è la principale proprietà di $p'$ nel presentare l'errore massimo ($\epsilon$ in questo caso) con segno alternato in esattamente $n+2$ valori nell'intervallo $[a,b]$. Il sistema è sempre risolvibile e la soluzione è unica (infatti la matrice associata è di Vandermonde, che è non singolare, si vedano [4] e [6]). Risolto il sistema, dunque, si ottengono $n+1$ coefficienti del polinomio $P(x)=p_0+p_1x+\cdots + p_nx^n$ di grado $n$ che può essere considerato una prima approssimazione di $p'$.
3. Calcolare gli $n+2$ estremi di $P-f$ e tornare al punto due usando gli estremi trovati, al posto dei vecchi $x_i$, per costruire il nuovo sistema.
Si può dimostrare che questo processo è convergente con $P \longrightarrow p'$ (si veda [6]). L'iterazione tra i punti due e tre ha fine quando si ritiene che la differenza tra gli estremi di $P-f$ sia piccola abbastanza da giustificare il fatto che $P$ sia diventato un buon candidato di $p'$. Si ricordi che $p'-f$ ha $n+2$ valori in cui si ottiene lo stesso errore massimo cambiato di segno. Dunque quando gli estremi di $P-f$ sono abbastanza vicini tra loro si è sicuri di aver trovato un buon candidato per $p'$.
Scrivere un programma per trovare i coefficienti di $P$ significa risolvere il sistema, calcolare gli estremi e tener conto di una miriade di dettagli e casi speciali. La cosa migliore da fare, dunque, è affidarsi ad un programma dedicato ed affidabile tipo Maple o Mathematica. Chi non ha la possibilità di utilizzare uno di questi due software può provare ad usare Sollya, programma gratuito che si trova in quasi tutte le maggiori distribuzioni GNU\Linux.
Cercare, ad esempio, un polinomio approssimante di grado $11$ in $[-\pi /2, \pi /2]$ per la funzione seno con Sollya restituisce i seguenti coefficienti:
p=fpminimax(sin(x),11,[|D...|],[-Pi/2;Pi/2], fixed, absolute);
p;
x * (0.99999999988985210919167911924887448549270629882812 + x^2 * (-0.166666665414392678457033980521373450756072998046875 + x^2 * (8.3333292644590661879533399769570678472518920898437e-3 + x^2 * (-1.9840702862772285897108304197899997234344482421875e-4 + x^2 * (2.7518855645336515181043068878352642059326171875e-6 + x^2 * (-2.379471364388763277020188979804515838623046875e-8))))))
supnorm(p,sin(x),[-Pi/2;Pi/2],absolute,2^(-40));
[1.32971831729691406029715287037401684717090712475884e-11;1.32971831729808563919394632854087727210908846430794e-11]
L'errore massimo è nell'intervallo $[1.329718317296E-11, 1.329718317298E-11]$, non male. Come si può notare il risultato è molto simile a quello del codice presentato in questo articolo ma non perfettamente uguale in quanto nel sistema usato per trovare i coefficienti usati nella funzione XMScalarSinCos sono stati aggiunti ulteriori vincoli. Uno è quello che la derivata di $P$ deve assumere valore $1$ in $x=0$, come del resto è per la funzione seno. L'altro vincolo fissa il valore ad uno dei due estremi dell'intervallo: considerando che la funzione seno è una funzione dispari $(f(-x)=-f(x))$ e che abbiamo fissato il primo vincolo per $x=0$, si sceglie $x=\pi /2$ come altro estremo dove il valore di $P$ (come per quello della funzione seno) deve valere $1$. Si noti che, essendo la funzione seno dispari, basta calcolare i suoi valori in $[0,\pi /2]$ per trovare, di conseguenza, quelli in $[-\pi /2,\pi / 2]$. Ecco perché si è scelto di mappare l'angolo passato alla funzione a $[-\pi /2,\pi /2 ]$. Riassumendo, i due vincoli aggiuntivi sono:
$P'(0)=f'(0)=(\sin(0))'=\cos(0)=1$
$P(\pi /2)=f(\pi /2)=\sin(\pi /2)=1$
e grazie a questi si riescono a calcolare coefficienti simili a quelli usati in XMScalarSinCos.
Con lo stesso procedimento si può cercare un polinomio approssimante di grado $10$ in $[-\pi /2, \pi /2]$ per la funzione coseno:
p=fpminimax(cos(x),10,[|32...|],[-Pi/2;Pi/2], fixed, absolute);
p;
0.99999999976716935634613037109375 + x^2 * (-0.49999999324791133403778076171875 + x^2 * (4.166663554497063159942626953125e-2 + x^2 * (-1.388835720717906951904296875e-3 + x^2 * (2.476014196872711181640625e-5 + x^2 * (-2.6053749024868011474609375e-7)))))
supnorm(p,cos(x),[-Pi/2;Pi/2],absolute,2^(-40));
[3.1967096142368210141690252469760085416491252052751e-10;3.1967096142396375486744909991642160243426747494351e-10]
Riferimenti:
[1] Elementary Functions: Algorithms and Implementation - Muller
[2] GPGPU Programming for Games and Science - Eberly
Nessun commento:
Posta un commento