Sqr関数高速化案
Sqr関数高速化案
関数sqr(function.sbp内)の高速化案です。
現関数Sqrの不具合:Sqr(x)のxが1から離れれば離れるほど計算時間がかかる。
改良案1(Sqrt関数)
改良案2(Sqrt4関数) 改良案1プラス除算を可能な限り削除。
(何故か改良案1のほうが速い。)
高速化の効果:
x^0.5 , pow(x,0.5) , Sqr(x) , Sqrt(x) , Sqrt4(x)について、
xの値を 1.001から10^300まで変えて比較した結果を示します。
使用CPU:Pentium4 1.9Ghz
計算時間は、1回当りのマイクロ秒。
x x^0.5 pow(x,0.5) Sqr(x) Sqrt(x) Sqrt4(x)
1.001 6.3 7.0 1.6 1.8 2.1
2 7.4 7.3 2.2 2.2 2.6
10^3 18.1 18.1 3.4 1.9 2.1
10^6 31.4 31.4 4.7 2.0 2.2
10^60 316 314 26.5 2.1 2.5
10^300 2100 2100 127 2.1 2.5
(注)字下げを行なうため、全角スペースを使っています。
Function Sqrt(number As Double) As Double
Dim s,last As Double
Dim i,j,jj,k as Word
If number=0 Then
Sqrt=0
Exit Function
ElseIf number>0 Then
s=number
i=Varptr(s)+6
jj=GetWord(i)
j=jj\32
k=jj mod 32
j=(j+511)*16+k
SetWord(i,j)
Do
last=s
s=(number/s+s)/2
Loop While s<>last
Else
'error
print
print "Error:Sqr argument is illegal"
Sqrt=-1
Exit Function
End If
Sqrt=s
End Function
Function Sqrt4(number As Double) As Double
'4次の漸化式により1/√Xを求め、X*1/√Xを計算することで√Xを求める。
Dim s,h As Double
Dim i,j,jj,k as Word
const eps=1e-4
If number=0 Then
Sqrt4=0
Exit Function
ElseIf number>0 Then
s=1/number
i=Varptr(s)+6
jj=GetWord(i)
j=jj\32
k=jj mod 32
j=(j+511)*16+k
SetWord(i,j)
Do
h=1-number*s*s
s=s*(1+h*(8+h*(6+5*h))/16)
Loop While Abs(h)>eps
Else
'error
print
print "Error:Sqr argument is illegal"
Sqrt4=-1
Exit Function
End If
Sqrt4=s*number
End Function
現関数Sqrの不具合:Sqr(x)のxが1から離れれば離れるほど計算時間がかかる。
改良案1(Sqrt関数)
改良案2(Sqrt4関数) 改良案1プラス除算を可能な限り削除。
(何故か改良案1のほうが速い。)
高速化の効果:
x^0.5 , pow(x,0.5) , Sqr(x) , Sqrt(x) , Sqrt4(x)について、
xの値を 1.001から10^300まで変えて比較した結果を示します。
使用CPU:Pentium4 1.9Ghz
計算時間は、1回当りのマイクロ秒。
x x^0.5 pow(x,0.5) Sqr(x) Sqrt(x) Sqrt4(x)
1.001 6.3 7.0 1.6 1.8 2.1
2 7.4 7.3 2.2 2.2 2.6
10^3 18.1 18.1 3.4 1.9 2.1
10^6 31.4 31.4 4.7 2.0 2.2
10^60 316 314 26.5 2.1 2.5
10^300 2100 2100 127 2.1 2.5
(注)字下げを行なうため、全角スペースを使っています。
Function Sqrt(number As Double) As Double
Dim s,last As Double
Dim i,j,jj,k as Word
If number=0 Then
Sqrt=0
Exit Function
ElseIf number>0 Then
s=number
i=Varptr(s)+6
jj=GetWord(i)
j=jj\32
k=jj mod 32
j=(j+511)*16+k
SetWord(i,j)
Do
last=s
s=(number/s+s)/2
Loop While s<>last
Else
'error
print "Error:Sqr argument is illegal"
Sqrt=-1
Exit Function
End If
Sqrt=s
End Function
Function Sqrt4(number As Double) As Double
'4次の漸化式により1/√Xを求め、X*1/√Xを計算することで√Xを求める。
Dim s,h As Double
Dim i,j,jj,k as Word
const eps=1e-4
If number=0 Then
Sqrt4=0
Exit Function
ElseIf number>0 Then
s=1/number
i=Varptr(s)+6
jj=GetWord(i)
j=jj\32
k=jj mod 32
j=(j+511)*16+k
SetWord(i,j)
Do
h=1-number*s*s
s=s*(1+h*(8+h*(6+5*h))/16)
Loop While Abs(h)>eps
Else
'error
print "Error:Sqr argument is illegal"
Sqrt4=-1
Exit Function
End If
Sqrt4=s*number
End Function
河川屋さん、アドバイスありがとうございます。
こちらのほうでも、速度検証をしてみました。結果は、ダントツでSqrt、Sqrt4が勝っていました。次回のバージョンアップでは、Sqrの演算内容をSqrtのものに変更させていただこうと思います。
s=number
i=Varptr(s)+6
jj=GetWord(i)
j=jj\32
k=jj mod 32
j=(j+511)*16+k
SetWord(i,j)
それはそうと、既存のSqr関数と河川屋さんのSqrt関数を見比べると、上記のようなコードが際立って見えますね。大きい桁を考慮して、Double型変数の指数部を直接的に(しかも整数演算のみで)計算するってのがポイントですね。
> (注)字下げを行なうため、全角スペースを使っています。
このフォーラムでは、半角スペース、タブ文字共に正常に表示できるように改良してありますので、そちらのほう、バシバシつかってもらって大丈夫です。
イグトランスさん:
> そういえばインラインアセンブラの搭載という話はどうなったのでしょうか。
> CPUの浮動小数点命令などが使えればより速くなりそうですが。
> あるいはコンパイラ組み込みの関数にするなんてのも。
>
> #最近ABへの要望ばっかりですいません。>山本様
インラインアセンブラというと、かなり大作業になってしまうので、スケジュールを調整しつつ、検討していきたいといったところでしょうか。
どっちみちVer5.0で64ビットコンパイラに対応するので、機械語コードの一新と同時に考えてみたいと思います。
こちらのほうでも、速度検証をしてみました。結果は、ダントツでSqrt、Sqrt4が勝っていました。次回のバージョンアップでは、Sqrの演算内容をSqrtのものに変更させていただこうと思います。
s=number
i=Varptr(s)+6
jj=GetWord(i)
j=jj\32
k=jj mod 32
j=(j+511)*16+k
SetWord(i,j)
それはそうと、既存のSqr関数と河川屋さんのSqrt関数を見比べると、上記のようなコードが際立って見えますね。大きい桁を考慮して、Double型変数の指数部を直接的に(しかも整数演算のみで)計算するってのがポイントですね。
> (注)字下げを行なうため、全角スペースを使っています。
このフォーラムでは、半角スペース、タブ文字共に正常に表示できるように改良してありますので、そちらのほう、バシバシつかってもらって大丈夫です。
イグトランスさん:
> そういえばインラインアセンブラの搭載という話はどうなったのでしょうか。
> CPUの浮動小数点命令などが使えればより速くなりそうですが。
> あるいはコンパイラ組み込みの関数にするなんてのも。
>
> #最近ABへの要望ばっかりですいません。>山本様
インラインアセンブラというと、かなり大作業になってしまうので、スケジュールを調整しつつ、検討していきたいといったところでしょうか。
どっちみちVer5.0で64ビットコンパイラに対応するので、機械語コードの一新と同時に考えてみたいと思います。
最後に編集したユーザー 山本 [ 2005年6月02日(木) 12:05 ], 累計 5 回
> > そういえばインラインアセンブラの搭載という話はどうなったのでしょうか。
確かに、インラインアセンブラはより高度な処理を可能にします。
しかし、果たしてABにインラインアセンブラが必要でしょうか。
アセンブラは、ドライバやOS、コンパイラなどを開発するごく限られた一部の人間以外は殆ど使用しません。(一部画像処理分野でも需要があるらしい。)
今時分アセンブラでプログラムを組むことも少なくなってきています。というか、無くなりました。
さらに、仮にAB上でアセンブラが必要になったとして、解決が困難かというとそうでもないと思えます。
必要であればMASMやNASMといったアセンブラを簡単に入手できますし、C++でインラインアセンブラを使用してDLLを作成しAB側でインポートする、といったことも可能ですから。
僕は、ABにインラインアセンブラを搭載する必要性は強くないと考えます。
> > CPUの浮動小数点命令などが使えればより速くなりそうですが。
現時点でのABの標準関数はすべてソースコード形式で配布されているため、コンパイル毎に機械語に翻訳されているというのが現状です。
如何に優れたアルゴリズムを利用しても、これでは・・・
FPUをダイレクトに叩ければ最速の処理速度が期待できますが、ソースコードレベルではFPUを扱うことは出来ません。そこで、
> そういえばインラインアセンブラの搭載という話はどうなったのでしょうか。
> あるいはコンパイラ組み込みの関数にするなんてのも。
こういった意見に行き着くわけですね。組み込み関数なら機械語レベルでの最適化も期待できますし、悪いことは無いわけです。
> どっちみちVer5.0で64ビットコンパイラに対応するので、機械語コードの一新と同時に考えてみたいと思います。
そのことですが、前から質問しようと思っていた事項があります。
この際ですので、こちらで質問させていただきます。
今現在64ビットの命令規格は様々なものがあります。
例えば、Itanium IA-64, AMD64, EM64T(AMD64との互換)などが主力でしょうか。
近い将来、少なくとも2年以内に64ビットの波が一般ユーザーの元にも押し寄せると僕は踏んでいますが、ABはどの規格に対応するのでしょうか?
まさか全部ってことはないでしょうw
それとも、未決定事項ですか?
確かに、インラインアセンブラはより高度な処理を可能にします。
しかし、果たしてABにインラインアセンブラが必要でしょうか。
アセンブラは、ドライバやOS、コンパイラなどを開発するごく限られた一部の人間以外は殆ど使用しません。(一部画像処理分野でも需要があるらしい。)
今時分アセンブラでプログラムを組むことも少なくなってきています。というか、無くなりました。
さらに、仮にAB上でアセンブラが必要になったとして、解決が困難かというとそうでもないと思えます。
必要であればMASMやNASMといったアセンブラを簡単に入手できますし、C++でインラインアセンブラを使用してDLLを作成しAB側でインポートする、といったことも可能ですから。
僕は、ABにインラインアセンブラを搭載する必要性は強くないと考えます。
> > CPUの浮動小数点命令などが使えればより速くなりそうですが。
現時点でのABの標準関数はすべてソースコード形式で配布されているため、コンパイル毎に機械語に翻訳されているというのが現状です。
如何に優れたアルゴリズムを利用しても、これでは・・・
FPUをダイレクトに叩ければ最速の処理速度が期待できますが、ソースコードレベルではFPUを扱うことは出来ません。そこで、
> そういえばインラインアセンブラの搭載という話はどうなったのでしょうか。
> あるいはコンパイラ組み込みの関数にするなんてのも。
こういった意見に行き着くわけですね。組み込み関数なら機械語レベルでの最適化も期待できますし、悪いことは無いわけです。
> どっちみちVer5.0で64ビットコンパイラに対応するので、機械語コードの一新と同時に考えてみたいと思います。
そのことですが、前から質問しようと思っていた事項があります。
この際ですので、こちらで質問させていただきます。
今現在64ビットの命令規格は様々なものがあります。
例えば、Itanium IA-64, AMD64, EM64T(AMD64との互換)などが主力でしょうか。
近い将来、少なくとも2年以内に64ビットの波が一般ユーザーの元にも押し寄せると僕は踏んでいますが、ABはどの規格に対応するのでしょうか?
まさか全部ってことはないでしょうw
それとも、未決定事項ですか?
当初の話題『Sqr関数高速化案』から大きくかけ離れた話題になってしまいましたね。別スレッドたてた方がいいのかな...(--;
> 64bitコンパイラになった場合、32bitCPUでの動作はどういうものになるのでしょうか?
> 1:高速になる
> 2:低速になる
> 3:変わらない
> 4:動作しない
お答えしておきます。答えは、「すべての可能性が否定できない」です。
まず、純粋な32ビットCPU&OSのペアであり、なおかつコンパイラが64ビットコードを生成した場合、「4.動作しない」が答えとなります。ま、これは説明不要、当然ですよね。
次に、本来64ビットを意識して生成されたEXEファイルだが、32ビットコードとの互換性を意識していた場合、「2.低速になる」ことは避けられないでしょう。
そして、完全に32ビットCPUを意識して生成された機械語コードであれば、32ビットと64ビット、どちらでも動作速度に大差はないと思われます。つまり、「3.変わらない」が答えとなります。
最後に、残された選択肢「1.高速になる」ですが、一目見て「んなあほな」と言いたくなります。しかし、大いにあり得ます。
少し込み入った話になりますが、かつてIntel社は現在のIA-32アーキテクチャの後継としてIA-64アーキテクチャを提唱しました。
しかしこのIA-64、あろうことかIA-32との互換性はありません。
取り敢えずシミュレーションによりIA-32コードも実行できます。ただし、CPU内部でシミュレートするためIA-32上で実行する方が高速になると言われています。
これが、この選択肢の可能性の所以です。
> 64bitコンパイラになった場合、32bitCPUでの動作はどういうものになるのでしょうか?
> 1:高速になる
> 2:低速になる
> 3:変わらない
> 4:動作しない
お答えしておきます。答えは、「すべての可能性が否定できない」です。
まず、純粋な32ビットCPU&OSのペアであり、なおかつコンパイラが64ビットコードを生成した場合、「4.動作しない」が答えとなります。ま、これは説明不要、当然ですよね。
次に、本来64ビットを意識して生成されたEXEファイルだが、32ビットコードとの互換性を意識していた場合、「2.低速になる」ことは避けられないでしょう。
そして、完全に32ビットCPUを意識して生成された機械語コードであれば、32ビットと64ビット、どちらでも動作速度に大差はないと思われます。つまり、「3.変わらない」が答えとなります。
最後に、残された選択肢「1.高速になる」ですが、一目見て「んなあほな」と言いたくなります。しかし、大いにあり得ます。
少し込み入った話になりますが、かつてIntel社は現在のIA-32アーキテクチャの後継としてIA-64アーキテクチャを提唱しました。
しかしこのIA-64、あろうことかIA-32との互換性はありません。
取り敢えずシミュレーションによりIA-32コードも実行できます。ただし、CPU内部でシミュレートするためIA-32上で実行する方が高速になると言われています。
これが、この選択肢の可能性の所以です。
自己レス。
1.プログラムの誤り。
指数を切り出す際に、桁指定を誤っていた。
(だからといって計算結果が間違っているわけではないが。)
2.毒を食らわば
√を求める範囲を√1~√2になるよう調整し、
第1近似を2次式で計算(誤差は最悪で1%くらい)したもの。
ただし、速度は前とあまり変わっていない。
実数演算しても指数部は全く動かないため、整数演算に置き換えれば
もっと速くなるはず。
3.指数、対数
こっちは、黙っていても訂正してもらえるものだと思っていたけど.....
exp(x)は、Abs(x)が相当大きい(X=100とか。)場合でないと効果が現れないから
直す意味は無いかもしれませんが、
Logは、Sqr以上に差が出ます。
なお、非正規数(約10^-308~10^-324)については、正しく計算しません。
尤も、以前のlogでは、この範囲の値を与えると永久ループに陥ります。
いずれにせよ、新旧ルーチンとも、IEEE754規格外の動作となっています。
Function Sqrt(number As Double) As Double
Dim s,x,last As Double
Dim i,j,j0,jj,k as Word
If number=0 Then
Sqrt=0
Exit Function
ElseIf number>0 Then
x=number
i=Varptr(x)+6
jj=GetWord(i)
j=jj\16
k=jj mod 16
SetWord(i,16*(1024-1)+k)
if x<>1 then
'第一近似値を求める
's=0.4125*x+0.5991 '1次式の場合
s=x*(-0.0705*x+0.6241)+0.4474 '2次式の場合
Do
last=s
s=(x/s+s)/2
Loop While last<>s
Else
s=1
End If
Else
'error
Exit Function
End If
i=Varptr(s)+6
k=GetWord(i) mod 16
j0=((j-1)\2+512)*16+k
SetWord(i,j0)
if j mod 2 =0 then s=s*_System_SQRT2
Sqrt=s
End Function
Function Exp_E(x As Double) As Double
Dim i,j,n As Long, k As Integer
'Dim E,Y As Double
Dim w,x2 As Double
If x>=0 Then
k=Int(x/_System_LOG2+0.5)
Else
k=Int(x/_System_LOG2-0.5)
End If
x=x-k*_System_LOG2
x2=x*x
w=x2/22
i=18
While i>=6
w=x2/(w+i)
i=i-4
Wend
Y=(2+w+x)/(2+w-x)
'e^x (-0.5<x<+0.5) を求める別案(こちらのほうが遅いため不採用)
'E=1
'Y=1
'n=1
'Do
'E=E*x/n
'Y=Y+E
'n=n+1
'Loop while Abs(E)>1e-16
i=Varptr(Y)+6
j=GetWord(i)+k*16
SetWord(i,j)
Exp_E=Y
End Function
Function Log_E(x As Double) As Double
Dim i,j,jj As Long, k As Long
Dim s As Double
If x<=0 Then
'error
Exit Function
End If
i=Varptr(x)+6
jj=GetWord(i)
k=jj\16-1023
j=jj mod 16
SetWord(i,16*1023+j)
x=x-1
s=0
i=9
While i>=1
s=i*x/(2+i*x/(2*i+1+s))
i=i-1
Wend
Log_E=_System_LOG2*k+x/(1+s)
End Function
1.プログラムの誤り。
指数を切り出す際に、桁指定を誤っていた。
(だからといって計算結果が間違っているわけではないが。)
2.毒を食らわば
√を求める範囲を√1~√2になるよう調整し、
第1近似を2次式で計算(誤差は最悪で1%くらい)したもの。
ただし、速度は前とあまり変わっていない。
実数演算しても指数部は全く動かないため、整数演算に置き換えれば
もっと速くなるはず。
3.指数、対数
こっちは、黙っていても訂正してもらえるものだと思っていたけど.....
exp(x)は、Abs(x)が相当大きい(X=100とか。)場合でないと効果が現れないから
直す意味は無いかもしれませんが、
Logは、Sqr以上に差が出ます。
なお、非正規数(約10^-308~10^-324)については、正しく計算しません。
尤も、以前のlogでは、この範囲の値を与えると永久ループに陥ります。
いずれにせよ、新旧ルーチンとも、IEEE754規格外の動作となっています。
Function Sqrt(number As Double) As Double
Dim s,x,last As Double
Dim i,j,j0,jj,k as Word
If number=0 Then
Sqrt=0
Exit Function
ElseIf number>0 Then
x=number
i=Varptr(x)+6
jj=GetWord(i)
j=jj\16
k=jj mod 16
SetWord(i,16*(1024-1)+k)
if x<>1 then
'第一近似値を求める
's=0.4125*x+0.5991 '1次式の場合
s=x*(-0.0705*x+0.6241)+0.4474 '2次式の場合
Do
last=s
s=(x/s+s)/2
Loop While last<>s
Else
s=1
End If
Else
'error
Exit Function
End If
i=Varptr(s)+6
k=GetWord(i) mod 16
j0=((j-1)\2+512)*16+k
SetWord(i,j0)
if j mod 2 =0 then s=s*_System_SQRT2
Sqrt=s
End Function
Function Exp_E(x As Double) As Double
Dim i,j,n As Long, k As Integer
'Dim E,Y As Double
Dim w,x2 As Double
If x>=0 Then
k=Int(x/_System_LOG2+0.5)
Else
k=Int(x/_System_LOG2-0.5)
End If
x=x-k*_System_LOG2
x2=x*x
w=x2/22
i=18
While i>=6
w=x2/(w+i)
i=i-4
Wend
Y=(2+w+x)/(2+w-x)
'e^x (-0.5<x<+0.5) を求める別案(こちらのほうが遅いため不採用)
'E=1
'Y=1
'n=1
'Do
'E=E*x/n
'Y=Y+E
'n=n+1
'Loop while Abs(E)>1e-16
i=Varptr(Y)+6
j=GetWord(i)+k*16
SetWord(i,j)
Exp_E=Y
End Function
Function Log_E(x As Double) As Double
Dim i,j,jj As Long, k As Long
Dim s As Double
If x<=0 Then
'error
Exit Function
End If
i=Varptr(x)+6
jj=GetWord(i)
k=jj\16-1023
j=jj mod 16
SetWord(i,16*1023+j)
x=x-1
s=0
i=9
While i>=1
s=i*x/(2+i*x/(2*i+1+s))
i=i-1
Wend
Log_E=_System_LOG2*k+x/(1+s)
End Function