三角関数・対数関数・指数関数の高速計算法
三角関数・対数関数・指数関数を高速に計算するコードを作成しました。いずれも，関数の変数の変域を小さい範囲に収めることにより，テーラー展開
の項の数を減らして高速化したものです。
三角関数
2倍角の公式を用いて三角関数を高速に計算する方法を考えました。c言語のdouble型で効果があります。 解説 cos.pdf
以下は，c言語でcosθを求める関数のコードで，math.hの三角関数より高速です。角の範囲が -2π≦θ≦2π
ならば，誤差の問題はありません。
また，若干遅くなりますが、コメントアウトした行をコードに組み込めば、任意の角についても求められます。
sinθを求めるには，piを円周率として，cos_fast(x-pi/2) または cos_fast(pi/2-x)
により求めるのがよいでしょう。
//高速三角関数
#define PAI 3.14159265358979323846
double cos_fast(double x)
{
static double
waru[5]={1.0/(3*4*5*6*7*8*9*10),-1.0/(3*4*5*6*7*8),1.0/(3*4*5*6),-1.0/(3*4),1.0};
double y,*p=waru;
int i,q;
//q=(int)(x/(2.0*PAI));
//x=x-q*(2.0*PAI);
x=x/32.0;
x=x*x;
y=-1.0/(3*4*5*6*7*8*9*10*11*12);
do{
y=y*x+(*p);
p++;
}while (p<waru+5);
y=y*x;
for (i=0;i<5;i++) y=y*(4.0-y);
return 1.0-y/2.0;
}
double sin_fast(double kaku) {return cos_fast(kaku-PAI/2);}
double tan_fast(double kaku) {return
cos_fast(kaku-PAI/2)/cos_fast(kaku);}
#undef PAI
対数関数
以下は，c言語で log(x)
を求める関数のコードで，関数値のテーブル引きを用いてテーラー展開の項の数を減らしました。math.hの対数関数より高速です。x=1の近傍での相対
誤差を少なくする工夫をしました。速度とメモリはトレードオフであり，テーブルのデータ数を増やせばさらに高速化できます。
実装では，まず，double型を指数部と仮数部に分離し，仮数部 1≦w＜2 に対して logw
の値を求めます。
logx のマクローリン展開はh=(x-z)/(x+z)として
log(x/z)=2（h+h3/3+h5/5+h7/7+h9/9
+…）
となることを用いています。
//高速対数関数
#define LOG2 0.693147180559945309417 //log(2)
double log_fast(double x)
{
static double table[17]={
.0
, // log( 16 /16)
.0606246218164348425806 , // log( 17 /16)
.1177830356563834545387 , // log( 18 /16)
.17185025692665922234 , // log( 19 /16)
.2231435513142097557662 , // log( 20 /16)
.2719337154836417588316 , // log( 21 /16)
.3184537311185346158102 , // log( 22 /16)
.3629054936893684531378 , // log( 23 /16)
.405465108108164381978 , // log( 24 /16)
.4462871026284195115325 , // log( 25 /16)
.4855078157817008078017 , // log( 26 /16)
.5232481437645478365168 , // log( 27 /16)
.5596157879354226862708 , // log( 28 /16)
.5947071077466927895143 , // log( 29 /16)
.6286086594223741377443 , // log( 30 /16)
.6613984822453650082602 , // log( 31 /16)
.6931471805599453094172 , // log( 32 /16)
};
unsigned long long w,kasuu16;
int q;
double y,h,z;
w=*(unsigned long long*)&x;
q=(((int)(w>>47)&0x1F)+1)>>1;
kasuu16=(w &
0xFFFFFFFFFFFFFULL)^0x4030000000000000ULL;//仮数*16
16<=kasuu16<32
h=*(double*)&kasuu16;
z=(double)(q+16);
h=(h-z)/(h+z);
z=h*h;
y=(2.0/9)*z+2.0/7;
y=y*z+2.0/5;
y=y*z+2.0/3;
y=y*z+2.0;
y=y*h;
return ((int)(w>>52)-1023)*LOG2+table[q]+y;
}
#undef LOG2
指数関数
以下は，c言語で exp(x)
を求める関数のコードで，関数値のテーブル引きを用いてテーラー展開の項の数を減らしました。math.hの指数関数より高速です。
//高速指数関数
#define LOG2 0.693147180559945309417 //log(2)
double exp_fast(double x)
{
static double
waru[7]={1.0/(2*3*4*5*6),1.0/(2*3*4*5),1.0/(2*3*4),1.0/(2*3),1.0/2,1.0,1.0};
static unsigned long long table[16]={
0x059B0D3158574ull, // 2^( 1
/32)-1
0x11301D0125B51ull, // 2^( 3
/32)-1
0x1D4873168B9AAull, // 2^( 5
/32)-1
0x29E9DF51FDEE1ull, // 2^( 7
/32)-1
0x371A7373AA9CBull, // 2^( 9
/32)-1
0x44E086061892Dull, // 2^(
11 /32)-1
0x5342B569D4F82ull, // 2^(
13 /32)-1
0x6247EB03A5585ull, // 2^(
15 /32)-1
0x71F75E8EC5F74ull, // 2^(
17 /32)-1
0x82589994CCE13ull, // 2^(
19 /32)-1
0x93737B0CDC5E5ull, // 2^(
21 /32)-1
0xA5503B23E255Dull, // 2^(
23 /32)-1
0xB7F76F2FB5E47ull, // 2^(
25 /32)-1
0xCB720DCEF9069ull, // 2^(
27 /32)-1
0xDFC97337B9B5Full, // 2^(
29 /32)-1
0xF50765B6E4540ull, // 2^(
31 /32)-1
};
double y=1.0/(2*3*4*5*6*7),*p=waru,z,r;
int q;
unsigned long long w;
z=x*(16.0/LOG2);
q=(int)z-(x<0);
r=x-((q<<1)+1)*(LOG2/32.0);
w=(unsigned long long)(1023+(q>>4))<<52
^ table[q & 0xF];
z=*(double*)&w;
do{
y=y*r+(*p);
p++;
}while (p<waru+7);
return y*z;
}
#undef LOG2
余弦・正弦関数
以下は，c言語でcosθ，sinθを同時に求める関数のコードで，速度優先で作りました。角の範囲が -π≦θ≦π
ならば，誤差の問題はありません。
//高速余弦正弦関数
void sincos(double kaku,double *x,double *y)
{
static double
waru[8]={-1.0/(3*4*5*6*7*8),-1.0/(2*3*4*5*6*7),1.0/(3*4*5*6),1.0/(2*3*4*5),-1.0/(3*4),-1.0/(2*3),1.0,1.0};
double
c=1.0/(3*4*5*6*7*8*9*10),s=1.0/(2*3*4*5*6*7*8*9),*p=waru,z;
int i;
kaku=kaku/32.0;
z=kaku*kaku;
do{
c=c*z+(*p);
p++;
s=s*z+(*p);
p++;
}while (p<waru+8);
c=c*z;
s=s*kaku;
for (i=0;i<5;i++){
s=s*(2.0-c);
c=c*(4.0-c);
}
*x=1.0-c/2.0;
*y=s;
}
ボックス・ミュラー法への応用 (Box–Muller transform)
seiki.c
ボックス・ミュラー法は，[0，1) の実数の一様乱数から，正規分布（ガウス分布）に従う擬似乱数を生成するのに応用される方法です。 解説 逆関数法とボックス・ミュラー法
高速な三角関数と対数関数のルーチンを用いることにより，math.hの組み込み関数を利用するときと比べて2倍ほど高速に正
規分布の乱数を発生できます。Box–Muller法のcos2πt，sin2πt
の部分を三角関数を使わずに済ますPolar法も組み込んでありますから，速度を比較してみて下さい。
乱数は平均が0，標準偏差が1の正規分布に従
いますので，平均m，標準偏差sの正
規分布の乱数を得るには，s*ransuu_seiki()+m として下さい。
なお，カイ２乗分布，指数分布，ポアソン分布，二項分布も組み込んであります。