関数の速度(3)

前々回、expがsin/cosより速いのが気になった。
また、logがatanより速いのも気になった。

exp(x) 35ns
sin(x) 45ns
cos(x) 45ns
log(x) 43ns
atan(x) 87ns

数学的には、expとsin/cosは同類である。

 \sin{x} = \frac{e^{ix} - e^{-ix}}{2i}
 \cos{x} = \frac{e^{ix} + e^{-ix}}{x}

logとatanもテーラー展開すれば似たようなものである。

 \log{(1 + x)} = x - \frac{x^2}{2} + \frac{x^3}{3} - \cdots
 \mbox{atan}{x} = x - \frac{x^3}{3} + \frac{x^5}{5} - \cdots

それでも全然速度が違うのは不思議だ。
テーラー展開より連分数のほうが速いらしいので、
実際にはそちらを使うのだろう。
ここに、logの連分数展開が載っている。


http://www.tac.tsukuba.ac.jp/~yaoki/conf.pdf

 \frac{\log(1 + x)}{x} = \frac{1}{1+\frac{\frac{x}{2}}{1+\frac{\frac{x}{2\times 3}}{1+\frac{\frac{2x}{2\times 3}}{1+\frac{\frac{2x}{2\times 5}}{1+\cdots}}}}}

これらを素直にコードにすると、こうなる。


// Taylor展開
double log2(double x) {
static const int N_log3 = 10;
x -= 1;
double y = x / N_log3;
for(int i = N_log3 - 1; i > 0; i--)
y = x * (1.0 / i - y);
return y;
}


// 連分数展開
double log3(double x) {
static const int N_log2 = 5;
x -= 1;
double y = 1 + (N_log2 + 1) / 2 * x / (4 * (N_log2 / 2) + 2);
for(int i = N_log2 - 1; i > 0; i--)
y = 1 + ((i + 1) / 2) * x / ((4 * (i / 2) + 2) * y);
y = x / y;
return y;
}

そのPDFには、テーラー展開10項までと連分数展開5項まででは、x = 1.1なら精度は同等だが、それ以上なら連分数展開のほうがよいとしてある。しかし、連分数展開は複雑で、テーラー展開41nsに対し、88nsほどかかる。
ただし、連分数の分母・分子に同じ多項式を掛けていき、(多項式)/(多項式)の形にする、すなわち、

 \frac{60x+60x^2+11x^3}{60+90x+36x^2+3x^3}

の形にしておくと、次のようなコードとなって、


double log4(double x) {
x -= 1;
return ((11 * x + 60) * x + 60) * x
/ (((3 * x + 36) * x + 90) * x + 60);
}

13ns程度となる。
連分数展開の精度がよい理由は、この分数をテーラー展開してみると分かる。

x-1/2x^2+1/3x^3-1/4x^4+1/5x^5-1/6x^6+57/400x^7-99/800x^8
+2603/24000x^9-1529/16000x^10+13513/160000x^11
-71779/960000x^12+212029/3200000x^13-375997/6400000x^14
+10004611/192000000x^15-1183341/25600000x^16
+52492101/1280000000x^17-279436751/7680000000x^18
+826445593/25600000000x^19-1466573977/51200000000x^20
+39038155127/1536000000000x^21-23092158589/1024000000000x^22
+204895498097/10240000000000x^23+...

6項まで正しい。これは当然で、元の分数で分母に4つ、分子に3つパラメータがあって、分数だから-1で6つパラメータがあることになるためである。しかし、第7項は57/400で1/7=57/399に非常に近い。それ以下も近いため、精度がよくなると思われる。連分数が精度がよいというのは、有理関数の精度が高いということだろう。


ただし、このときxが[1, 2]で誤差が最大2.5e-5程度となる。
精度をどの程度にするかによるが、1e-7程度とするとあと3項必要で、

 \frac{1890x+6510x^2+6720x^3+2020x^4+48x^5}{945+4200x+6300x^2+3600x^3+600x^4}

となる。
これをテーラー展開すると、

x-1/2x^2+1/3x^3-1/4x^4+1/5x^5-1/6x^6+1/7x^7-1/8x^8+1/9x^9
-3175/31752x^10+265/2916x^11-428195/5143824x^12
+6216415/81015228x^13-14807875/208324872x^14
+433798525/6562233468x^15+...

これで、18ns程度となる。