少し前に別のスレッドでsinとcosをやったので
ついでに数学関係の関数をやってみました。
absをSealさんとは別の方法でやってますが特に他意はありません。
absfをモジュール内の色んな所で使っているので
そのおまけでabsもabsfと同じ方法で入れてみただけです。
M_PI（円周率）やdeg2rad、rad2deg（角度変換系）はそのまま使っていますが
これは定数とかマクロなんで、新たに定義してもコピペに近くなりそうなんで
そのまま使わせてもらいました。
ちなみに、当然ですが本家HSPと比較すると若干速度が遅いです。
精度も微妙……と言いたかったんですが、
試したところ場合によっては本家より正確な数値が出ることもありました。
まあほとんどは下１〜２桁（小数点以下15〜16桁程度）に誤差が出てしまうんですが
その辺は見逃してください。<(_ _)>
#module
// abs
#undef abs
#defcfunc abs int p1
if( p1 < 0 ) {
return -p1;
}
return p1;
// absf
#undef absf
#defcfunc absf double p1
if( p1 < 0.0 ) {
return -p1;
}
return p1;
// limit
#undef limit
#defcfunc limit int p1, int p2, int p3
if(p1 < p2):return p2;
if(p1 > p3):return p3;
return p1;
// limitf
#undef limitf
#defcfunc limitf double p1, double p2, double p3
if(p1 < p2):return p2;
if(p1 > p3):return p3;
return p1;
// logf テイラー展開でlogfの計算 0や負数のときの戻り値はHSPに準拠
#undef logf
#defcfunc logf double p1
if( p1 < 0.0 ) {
prm01 = 0xFFF80000;
drc = 0.0;
lpoke drc, 4, prm01;
}
else:if( p1 == 0.0 ) {
prm01 = 0xFFF00000;
drc = 0.0;
lpoke drc, 4, prm01;
}
else:if( p1 == 2.0 ) {
drc = 0.693147180559945309;
}
else:if( p1 < 2.0 ) {
drc = _logf(p1);
}
else {
prm01 = 1.0;
prm02 = 1.0;
repeat -1
prm02 *= 2.0;
prm03 = p1 / prm02;
if( prm03 <= 2.0 ):break;
prm01 += 1.0;
loop
drc = prm01 * 0.693147180559945309 + _logf(prm03);
}
return drc;
// logfのサブルーチン
#defcfunc _logf double p1
prm_01 = p1 - 1.0;
prm_02 = prm_01;
iprm_01 = 2;
repeat -1
prm_02 *= p1 - 1.0;
prm_03 = prm_02 / iprm_01;
if( absf(prm_03) < 0.0000000000000001 ):break;
if( (iprm_01\2) == 0 ) {
prm_01 -= prm_03;
}
else {
prm_01 += prm_03;
}
iprm_01++;
loop
return prm_01;
// expf マクローリン展開でexpfの計算 オーバーフローやアンダーフローの精度は標準関数より若干低い
#undef expf
#defcfunc expf double p1
if( absf(p1) < 0.00000000000000000001 ):return 1.0;
if( p1 < -745.0 ):return 0.0; // アンダーフロー１
if( p1 > 707.0) { // オーバーフロー
prm01 = 0x7FF00000;
drc = 0.0;
lpoke drc, 4, prm01;
return drc;
}
prm_d = 1.0;
prm_e = 1.0;
drc = 1.0;
lflag = 0;
prmlast = 0.0; // オーバーフローの備え
prm00 = 0x7FF00000; //
lpoke prmlast, 4, prm00; //
repeat -1, 1
await 0
prm_d = drc;
prm_e = prm_e * absf(p1) / double(cnt);
if( prmlast == prm_e ) {
if(p1 <= 0.0) {
drc = 1.0 / drc;
}
lflag = 1;
break;
}
drc += prm_e;
if( (absf(drc-prm_d)/absf(prm_d)) < 0.00000000000000000001 ) {
if(p1 <= 0.0) {
drc = 1.0 / drc;
}
lflag = 1;
break;
}
loop
if( lflag == 1 ) {
return drc;
}
else {
return 0.0; // アンダーフロー２
}
// powf expfとlogfを利用して計算
#undef powf
#defcfunc powf double p1, double p2
if( p1 == 0.0 ) { // 底が0の場合
if( p2 == 0.0 ) { // 指数が0ならば1
return 1.0;
}
else:if( p2 < 0.0 ) { // 指数が負数ならばエラー
pdrc = 0.0;
prm01 = 0x7FF00000;
lpoke pdrc, 4, prm01;
return pdrc;
}
else { // その他の場合は0
return 0.0;
}
}
if( p1 < 0.0 ) { // 底が負数の場合
if( p2 == 0.0 ) { // 指数が0ならば1
return 1.0;
}
else:if( (p2\1.0) == 0.0 ) { // 指数が整数ならばループで計算
pdrc = p1;
repeat int(absf(p2))-1
pdrc *= p1;
loop
return pdrc;
}
else { // その他の場合はエラー
pdrc = 0.0;
prm01 = 0xFFF80000;
lpoke pdrc, 4, prm01;
return pdrc;
}
}
if( (p2\1.0) == 0.0 ) { // 指数が整数ならばループで計算
pdrc = p1;
repeat int(absf(p2))-1
pdrc *= p1;
loop
}
else {
pdrc = expf( p2 * logf(p1) );
}
return pdrc;
// sqrt 平方根の計算方法を利用
#undef sqrt
#defcfunc sqrt double p1
if( p1 == 0.0 ):return 0.0; // 0の場合
if( p1 < 0.0 ) { // 負数の場合はエラー
drc = 0.0;
prm01 = 0xFFF80000;
lpoke drc, 4, prm01;
return drc;
}
prm01 = p1 / 2.0;
prm02 = 0.0;
drc = 1.0;
repeat -1
prm02 = drc;
drc = (drc + p1 / drc) / 2.0;
if( drc == prm02 ):break;
loop
return drc;
// sin テイラー展開でsinの計算
#undef sin
#defcfunc sin double p1
prm01 = p1;
prm02 = prm01;
iprm01 = 1;
drc = prm01;
repeat -1
prm02 *= ( -prm01*prm01 / (2.0*iprm01*(2.0*iprm01+1.0)) );
drc += prm02;
iprm01++;
if( absf(prm02) <= 0.00000000000000000001 ):break;
loop
return drc;
// cos テイラー展開でcosの計算
#undef cos
#defcfunc cos double p1
prm01 = p1;
prm02 = 1.0;
iprm01 = 1;
drc = 1.0;
repeat -1
prm02 *= ( -prm01*prm01 / ((2.0*iprm01-1.0)*2.0*iprm01) );
drc += prm02;
iprm01++;
if( absf(prm02) <= 0.00000000000000000001 ):break;
loop
return drc;
// tan 単純にsin/cosを返す
#undef tan
#defcfunc tan double p1
return sin(p1) / cos(p1);
// atan テイラー展開でatanの計算 戻り値の範囲は -パイ〜パイ（度数にすると -180度より大きく180度以下）
#undef atan
#defcfunc atan double p1, double p2
H_PI = M_PI / 2.0;
Q_PI = M_PI / 4.0;
if( p2 == 0.0 ) {
if( p1 > 0.0 ):return H_PI; // Xが0で Yが正数
if( p1 < 0.0 ):return -H_PI; // Xが0で Yが負数
if( p1 == 0.0 ):return 0.0; // 両方0
}
atdrc = 0.0;
root2 = sqrt(2.0);
root2m1 = root2 - 1.0;
root2p1 = root2 + 1.0;
if( p2 < 0.0 ) {
// X値が負数ならば第二象限か第三象限
if( p1 < 0.0 ) {
// Y値も負数ならば第三象限
Quad = 3;
}
else {
// Y値が正数ならば第二象限
Quad = 2;
}
}
else {
// X値が正数ならば第一象限か第四象限
if( p1 < 0.0 ) {
// Y値が負数ならば第四象限
Quad = 4;
}
else {
// Y値が正数ならば第一象限
Quad = 1;
}
}
x = p1 / p2;
if( x < 0.0 ) {
x = -x;
x2 = -p2 / p1;
}
else:if( p1 != 0 ) {
x2 = p2 / p1;
}
if( x != 0.0 ) {
// 精度を保ちつつ高速化するために4つの条件に分割する
if( (0.0 < x) && (x <= root2m1) ) {
// xが 0 〜 √2-1 の範囲の場合
atdrc = _atan( x );
}
else:if( (root2m1 < x) && (x < 1.0) ) {
// xが √2-1 〜 1 の範囲の場合
atdrc = Q_PI - _atan( (1.0-x) / (1.0+x) );
}
else:if( (1.0 <= x) && (x < root2p1) ) {
// xが 1 〜 √2+1 の範囲の場合
atdrc = Q_PI + _atan( (1.0-x2) / (1.0+x2) );
}
else:if( x >= root2p1 ) {
// xが √2+1 以上の場合
atdrc = H_PI - _atan( x2 );
}
}
else {
atdrc = 0.0;
}
switch Quad
case 2
atdrc = M_PI - atdrc;
swbreak;
case 3
atdrc = atdrc - M_PI;
swbreak;
case 4
atdrc = -atdrc;
swbreak;
default
swbreak;
swend
return atdrc;
// atanのサブルーチン
#defcfunc _atan double p1
sdrc = 0.0;
prm02 = 1.0;
prm03 = 0.0;
prm04 = -1.0;
prm05 = p1;
prm06 = p1*p1;
prmlast = 0.0; // オーバーフローの備え
prm00 = 0x7FF00000; //
lpoke prmlast, 4, prm00; //
repeat -1
prm04 = -prm04;
prm03 = prm04 * prm05 / prm02;
sdrc += prm03;
if( absf(prm03) == prmlast ):break; // オーバーフロー
if( absf(prm03) < 0.00000000000000000001 ):break;
prm05 *= prm06;
prm02 += 2.0;
loop
return sdrc;
#global
log = logf(5.0);
mes strf("5の対数:\t\t\t\t\t%.16f", log);
exp = expf(log);
mes strf("5の対数の指数 精度が高いと5.0になる:\t\t%.16f", exp);
pow = powf(2.0, 10.0);
mes strf("2の10乗 整数乗のときはループで掛けてるだけ:\t%.16f", pow);
root3 = sqrt(3.0);
mes strf("\n3の平方根（sqrt(3.0)の場合）:\t\t\t%.16f", root3);
root3_2 = powf(3.0, 0.5);
mes strf("3の平方根（powf(3.0, 0.5)の場合）:\t\t%.16f", root3_2);
sin60 = sin(deg2rad(60.0));
cos60 = cos(deg2rad(60.0));
tan60 = tan(deg2rad(60.0));
mes strf("\nsin60度 = √3 / 2:\t\t\t%.16f", sin60);
mes strf("cos60度 = 1 / 2:\t\t\t%.16f", cos60);
mes strf("tan60度 = √3 / 1:\t\t\t%.16f", tan60);
arctan = atan(root3, 1.0);
mes strf("\natan(√3 , 1):\t\t\t\t%.16f", arctan);
datan = rad2deg(arctan);
mes strf("atan(√3 , 1)を度数に変換:\t\t%.16f", datan);
グラフィック系もいくつかはAPIを多用すれば比較的簡単にできそうなんですが
redrawで更新しないと画面に反映されなそうなんで、ちょっと迷ってます。
atanに変なバグがあったのをコッソリ修正。
HSPの命令でHSPの命令を作るスレ
