by 河川屋 » 2005年6月09日(木) 21:39
自己レス。
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