関数の速度（3） - inamori’s diary

前々回、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程度となる。