//ICos(a)=163*cos(a*PI/512)
int ICos(int a)
{
int ratio=1;
unsigned int b,c;
///////////////////
// 1.引数の整形 //
///////////////////
if(a<0){ //引数が負なら引数反転
a=-a;
}
a%=1024; //引数を0〜2PIに
// a&=0x3ffff
if(a>=512){ //引数>=PIなら引数をPIシフトし、戻り値反転
a-=512;
ratio=-1;
}
if(a==256) return 0;
if(a>256){ //引数>=PI/2なら引数を鏡像にし、戻り値反転
a=512-a;
ratio=-ratio;
}
///////////////////////
// 2.守備範囲の判定 //
///////////////////////
//およそπ/4以上の領域では、sinで計算した方が高精度
// if(a>104) return ratio*ISin(256-a); //閾値(104)は実計算で決定
if(a>104) return ratio*ISinR(256-a); //閾値(104)は実計算で決定
/////////////////////////////
// 3.マクローリン式の計算 //
/////////////////////////////
b=a*a;
/*
// aは8bit以下となった
// ４項６次までテイラー展開。(演算子^はべき乗を示す。a^b=pow(a,b))
// 163 -163/2!*(a/163)^2 +163/4!*(a/163)^4 -163/6!*(a/163)^6
// =163 -a^2/2/163 +a^4/24/163^3 -a^6/720/163^5
// =163 -a^2/2/163 +a^4/24/163^3 -(a^4/163^3)*(a^2/163)/720/163
// アンダーフロー抑制//分子を変数範囲いっぱいに拡大
// ={163^4*2+a^4/24*2-[a^2*163^2 +(a^4/163^2)*(a^2)/720*2]}/163/163/163/2
// ={1411823522+a^4/12-[a^2*26569 +(a^4/26569)*(a^2)/360]}/8661494 //31bit/24bit
// unsigned32で作業 a^4+αまでオーバーフロー回避
// より高速/高精度が要求される場合、引数のレンジ別に使用する項数(次数)を使い分ける
c=(1411823522+b*b/12-(b*26569 +(b*b/26569)/360*b))/8661494; //４項６次
c=(1411823522+b*b/12-b*26569)/8661494; //３項４次
c=(1367343104-b*(25731-13*b/161))>>23; //３項４次分母を強引に２のべき乗化
c=(1411823522-b*26569)/8661494; //２項２次
*/
c=(1367343104-b*25731)>>23; //２項２次分母を強引に２のべき乗化
ratio*=c; //cにratioを掛けて終了
return ratio;
}
//ISin(a)=163*sin(a*PI/512)
int ISin(int a)
{
int ratio=1;
unsigned int b,c,d;
///////////////////
// 1.引数の整形 //
///////////////////
if(a<0){ //引数が負なら引数と戻り値反転
a=-a;
ratio=-1;
}
a%=1024; //引数を0〜2PIに
// a&=0x3ffff
if(a>=512){ //引数>=PIなら引数をPIシフトし、戻り値反転
a-=512;
ratio=-ratio;
}
if(a>256){ //引数>=PI/2なら引数を鏡像に
a=512-a;
}
///////////////////////
// 2.守備範囲の判定 //
///////////////////////
//およそπ/4以上の領域では、cosで計算した方が高精度
// if(a>151) return ratio*ICos(256-a);
if(a>151) return ratio*ICosR(256-a);
/////////////////////////////
// 3.マクローリン式の計算 //
/////////////////////////////
b=a*a;
d=a;
/*
// aは8bit以下となったた
// ３項５次までテイラー展開。(演算子^はべき乗を示す。a^b=pow(a,b))
// a -163/3!*(a/163)^3 +163/5!*(a/163)^5
// =a -a^3/6/163^2 +a^5/120/163^5
// アンダーフロー抑制//分子を変数範囲いっぱいに拡大
// =[a*163^3*2 -a^3/3*163 +a^5/60/163^2]/163^3/2
// =[a*163^3*2*3 -a^3*163 +(a^4*3/163^2)*a/60]/163^3/2/3
// =[a*25984482 -a^3*163 +(a^4*3/26569)*a/60]/25984482 //32bit/25bit
// unsigned32で作業 a^4+αまでオーバーフロー回避
// より高速/高精度が要求される場合、引数のレンジ別に使用する項数(次数)を使い分ける
c=(d*25984482 )/25984482; //１項１次
c=(d*25984482 -d*b*163)/25984482; //２項３次
c=(d*25984482 +(b*b*3/26569)*d/60 -d*b*163)/25984482; //３項５次
c=(d*16777216-b*d*105)>>24; //２項３次分母を強引に２のべき乗化
*/
c=(d*(16777216-b*105))>>24; //２項３次分母を強引に２のべき乗化
ratio*=c; //cにratioを掛けて終了
return ratio;
}
//引数チェック省略版。ISin()以外から呼ぶのは危険。
int ICosR(int a)
{
unsigned int b,c;
b=a*a;
c=(1367343104-b*25731)>>23; //２項２次分母を強引に２のべき乗化
return (int)c;
}
//引数チェック省略版。ICos()以外から呼ぶのは危険。
int ISinR(int a)
{
unsigned int b,c,d;
b=a*a;
d=a;
c=(d*(16777216-b*105))>>24; //２項３次分母を強引に２のべき乗化
return (int)c;
}
いかがでしょうか。
精度的には、最大で ±1 ( ±1/±163 ＝ 0.6% ) の誤差が発生します。
次数を上げれば、速度と引換えに精度は向上します。
ここで処理を任された側の関数は、すでに ０〜Π/２ の範囲に集約され
た引数を受取ります。
そこで、引数の整形を省略した関数( int ICosR(int), ISinR(int) )を
別途準備し、若干の速度アップを図りました。関数の数が増えて煩雑に思
われる場合は、ICos(int), ISin(int)に任せても良いです。
半径163pxの円弧を、純正の実数関数(白)とICos,ISin関数(黒)
で描いたものです。π/4の前後で1pxのずれが発生しています。
速度は、実数( double/64bit )の sin, cos に比べ、６ 倍くらい出ます。
手元の Pen4-2.26GHz DualCH-DDR266 では、５８ 百万回/sec. 程度です。
データ長が半分である事を差し引いてもまずまずの結果と思います。
更なるスピードをお望みの場合は・・・テーブルですね。
トリッキーな内容はありませんが、注意すべき点を ２ つ
まず、% (剰余)演算子の扱いについて、第一 and/or 第二オペランドが負数
の場合は不定義のはずですが、MicroSoft VisualC++ 6.0は独自に定義して
値を返します。両オペランドの絶対値の剰余に、第一オペランドの符号が
付いた値が返ります。例えば、(-5)%3 == -(5%3) == -2 です。
mod(a,b)= a<0 ? -(abs(a)%abs(b)) : abs(a)%abs(b); なようです。
同じ事を Excel2000 で試してみると、(-5)%3 == 1 です。
これは、-5 = 3*(-2) + 1 との解釈です。
感覚的には Excel の方がしっくりくるし、今回作った関数の中でも負数の
反転が不用になって好都合なのですが。
次に、int の演算順序には注意してください。
2*3/4 == 1 ですが、2/4*3 == 0 です。
そんな事はわかっておる、と怒られそうですが、私はよく間違えます。
おまけ：正しく接する
