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