float-funcs
float 演算関数の検討
SR-プロセッサに fpu を搭載したが，newlib/math の初等関数を用いると実行時間がかなり大きい。マウスのコントロールなどの用途では精度をある程度犠牲にしても高速な演算をしたいことがある。
そこで、いくつかの関数を独自に検討、実装し、精度と速度の兼ね合いを取れるようにした。
具体的には 平方根、三角関数、指数関数、対数関数を実装した。libm の関数名とかぶらないように、例えば sqrtf() --> fsqrt() のように fで始まる関数名とした。
また、簡単のため値域のチェックなど一切行っていない。
平方根 fsqrt(x)
- アルゴリズム
平方根はニュートン法で求めるが，平方根をそのままニュートン法の漸化式にすると除算が必要となる。そのため逆数平方根 \( \frac{1}{\sqrt(x)} \) をニュートン法で求めたのちに \( \sqrt{x} = x \frac{1}{\sqrt{x}} \) とする。
また、ニュートン法の場合，初期値が真値に近いほど少ない繰り返しで精度が上がるため、初期値の求め方を工夫する。
- 漸化式
\[ f(x) = \frac{1}{x^2} - a, f(x) = 0 の根は \frac{1}{\sqrt{a}}\\ ニュートン法の漸化式：\\ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\\ x_{n+1} = \frac{1}{2} x_n ( 3 - {a}{x_n}^2) \]
- 初期値
ieee754 の 32bit float のバイナリ表現をうまく利用して、整数演算で \( \frac{1}{\sqrt{a}} \) の近似値を求める。
- code
typedef union {float f; unsigned u;} fu_t; float invsqrt(float a) { unsigned ix = (0xbe800000 - ((fu_t)a).u) >> 1; float x = ((fu_t)ix).f; // Xn+1 = Xn * (3 - A * Xn^2) / 2 x = x * (3.0f - a * x * x) / 2; x = x * (3.0f - a * x * x) / 2; x = x * (3.0f - a * x * x) / 2; return x; } float fsqrt(float x) { return x * invsqrt(x); }
三角関数 fsin(x) fcos(x) fsincos(x)
- アルゴリズム
- 周期関数なので、\( -\pi \le x \lt \pi \) の範囲の計算を行う
\[ n = \lfloor \frac{2 x}{\pi} \rfloor\\ b = x - \frac{n \pi}{2} \small{0 \le b \lt \frac{\pi}{2}}\\ \cos{x} = \sin{(x + \frac{\pi}{2})} \sin{(-x)} = -\sin{x}\\ を利用して、\sin{b}, \cos{b} から -\pi \le x \lt \pi の範囲に \]
- \( \sin{b}, \cos{b} \) を Maclaurin 展開で近似
少ない項数で近似誤差を最小にするような係数の補正を行った。\[ \sin{b} = b - \frac{b^3}{3!} + \frac{b^5}{5!} - \frac{b^7}{7!} + \cdots\\ \cos{b} = 1 - \frac{b^2}{2!} + \frac{b^4}{4!} - \frac{b^6}{6!} + \cdots\\ \]
- code
#define Cs1 (float)(1.0/(2*3) - 38e-6) #define Cs2 (float)(1.0/(2*3*4*5) - 180e-6) #define Cc1 (float)(1.0/(2) + 4e-6) #define Cc2 (float)(1.0/(2*3*4) + 15e-6) #define Cc3 (float)(1.0/(2*3*4*5*6)) void fsincos(float x, float *s, float *c) { int sgn = x < 0.0f; x = sgn ? -x : x; int ix = ((int)(x * (float)(4/M_PI)) + 1) >> 1; x -= ix * (float)(M_PI/2); float xx = x*x; float ss, cc; ss = x * (1.0f - xx * (Cs1 - xx * Cs2)); cc = (1.0f - xx * (Cc1 - xx * (Cc2 - xx * Cc3))); if(sgn){ switch(ix & 0x3){ case 0: *s = -ss; *c = cc; break; case 1: *s = -cc; *c = -ss; break; case 2: *s = ss; *c = -cc; break; case 3: *s = cc; *c = ss; break; } }else{ switch(ix & 0x3){ case 0: *s = ss; *c = cc; break; case 1: *s = cc; *c = -ss; break; case 2: *s = -ss; *c = -cc; break; case 3: *s = -cc; *c = ss; break; } } } float fsin(float x) { int sgn = x < 0.0f ? -1 : 1; x *= sgn; int ix = ((int)(x * (float)(4/M_PI)) + 1) >> 1; x -= ix * (float)(M_PI/2); float xx = x*x; sgn = (ix & 0x2) ? -sgn : sgn; switch(ix & 0x1){ case 0: return (x * (1.0f - xx * (Cs1 - xx * Cs2)))*sgn; case 1: return (1.0f - xx * (Cc1 - xx * (Cc2 - xx * Cc3)))*sgn; } } float fcos(float x) { x = x < 0.0f ? -x : x; int ix = ((int)(x * (float)(4/M_PI)) + 1) >> 1; x -= ix * (float)(M_PI/2); float xx = x*x; switch(ix & 0x3){ case 0: return (1.0f - xx * (Cc1 - xx * (Cc2 - xx * Cc3))); case 1: return -(x * (1.0f - xx * (Cs1 - xx * Cs2))); case 2: return -(1.0f - xx * (Cc1 - xx * (Cc2 - xx * Cc3))); case 3: return (x * (1.0f - xx * (Cs1 - xx * Cs2))); } }
指数関数 fexp(x)
- アルゴリズム
- 算法変換
\( \mathrm{e}^{x} \) を２のべき乗と\( \pm \frac{1}{2}\log{2} \) の範囲の指数関数\( \mathrm{e}^{b} \) の積に変換する。\[ \mathrm{e}^{x} = 2^n \mathrm{e}^{b} \small{-\frac{1}{2}\log{2} \le b \le \frac{1}{2}\log{2}}\\ n = \lfloor \frac{x}{\log{2}} \rfloor\\ b = x \bmod \log{2} \] - \( \mathrm{e}^{b} \) を Maclaurin 展開で計算する。
\[ \mathrm{e}^{x} = 1 + \frac{x}{1!} + \frac{x^2}{2!} + \cdots + \frac{x^n}{n!} + \cdots \]
- \( 2^n \mathrm{e}^{b} \) を計算
- code
#include <stdio.h> #include <math.h> typedef union {float f; unsigned u;} fu_t; #define LN2 ((float)M_LN2) #define INVLN2 ((float)(1.0/M_LN2)) float fexp(float x) { float a, b, expb; int n; //const float cf5[] = {9.9999985252e-01,4.9999200435e-01,1.6667133594e-01,4.1890954711e-02,8.3186421629e-03,}; const float cf6[] = {1.0000000000e+00,4.9999996902e-01,1.6666648685e-01,4.1669474588e-02,8.3571913396e-03,1.3624404424e-03,}; // 1 1/2! 1/3! 1/4! 1/5! 1/6! // range reduction -0.5ln2 < b < 0.5ln2 a = x * INVLN2 + (x < 0.0f ? -0.5f : 0.5f); n = (int)a; b = x - n * LN2; // calc exp(b) by Maclaurin (+ coeff adjust) // expb = 1.0f+b*(cf5[0]+b*(cf5[1]+b*(cf5[2]+b*(cf5[3]+b*(cf5[4]))))); expb = 1.0f+b*(cf6[0]+b*(cf6[1]+b*(cf6[2]+b*(cf6[3]+b*(cf6[4]+b*(cf6[5])))))); // exp(x) = 2^n * exp(b) add n to ieee754 exponent direct unsigned iy = (((fu_t)expb).u + (n << 23)) & 0x7fffffff; // multiply 2^n return ((fu_t)iy).f; }
対数関数 flog(x)
- アルゴリズム
- 算法変換
\( n \) 及び \( 1+f \) は ieee754 のバイナリ表現から整数演算で直接得ることが出来る。\[ a = 2^n (1+f) \small{0 \le f \lt 1}\\ \log{a} = \log_2{2^n} \log{2} + \log{(1+f)}\\ = n \log{2} + \log{(1+f)} \]
- \( \log{(1+f)} \) をMaclaurin展開で計算する。 \( f \) の範囲は \( \small{0 \le f \lt 1} \)
収束の高速化のため、変数変換(除算が一回必要)。
少ない項数で誤差を最小にするよう、非線形最適化で係数を補正した。\[ x = 1 + f\\ y = \frac{(x-1)}{(x+1)}\\ \frac{\log{x}}{2} = y + \frac{y^3}{3} + \cdots + \frac{y^{2n-1}}{2n-1} + \cdots \] - \( \log{a} = n \log{2} + \log{(1+f)} \) を計算
- code
#include <stdio.h> #include <math.h> typedef union {float f; unsigned u;} fu_t; #define LN2 ((float)M_LN2) float flog(float a) { int n = (int)( (((fu_t)a).u << 1) - 0x7f000000 ) >> 24; // get ieee754 8bit exponent (2^n) unsigned ix = (((fu_t)a).u & 0x7fffff) | 0x3f800000; // get 23bit fraction (1+frac) float x = ((fu_t)ix).f; float y = (x - 1.0f) / (x + 1.0f); // 除算1回 float f, yy = y * y; const float cf5[] = {2.0000000764e+00,6.6665553119e-01,4.0044387116e-01,2.7786767296e-01,2.8657712787e-01,}; // const float cf6[] = {2.0000000000e+00,6.6666561224e-01,4.0004315684e-01,2.8480352473e-01,2.3043217603e-01,1.7447693180e-01,}; f = y *(cf5[0] + yy*(cf5[1] + yy*(cf5[2] + yy*(cf5[3] + yy*cf5[4])))); // ほぼfloat 限界 // f = y *(cf6[0] + yy*(cf6[1] + yy*(cf6[2] + yy*(cf6[3] + yy*(cf6[4] + yy*cf6[5]))))); // float 限界 return n * LN2 + f; }
