３．逆三角関数を造る　Ark Tangent & PI

三角関数の次は、逆三角関数です。
三角関数は角度から座標を求め、逆三角関数は座標から角度を求める用途が多い
と思います。
前者の結果は (cosθ, sinθ)、後者の結果は atan(y/x) といった形になるので、
三角関数は cos, sin を造りましたが、逆三角関数は atan を造ってみます。
必要なら cos, sin から tan を計算出来たように、atan から asin, acos を
計算することも容易です。

tan では微分に体力が必要でした。今回は少し楽をします。

tan の逆関数を f(x) とすると、
f(x)=atan(x)　つまり、x=tan(f(x))　です。
逆関数の微分係数は、元の関数の微分係数の逆数ですから、


df(x)/dx = ( dx/df(x) )-1
　　　　 = ( tan'(f(x)) )-1
　　　　 = ( 1/cos2(f(x)) )-1
　　　　 = ( (cos2(f(x))+sin2(f(x)))/cos2(f(x)) )-1
　　　　 = ( 1+tan2(f(x)) )-1
　　　　 = 1/(1+x2)


この式の形は、なんだか見覚えがあります。等比級数ですね。
  

    

      　　　　 Σ
       ∞
      ak
      = 1+a+a2+a3+a4+・・　　　 　　[1]　
    

    

       k=0
    

    
　

    

      　　　　aΣ
       ∞
      ak
      = a+a2+a3+a4+a5+・・　　　 a*[1]　
    

    

       k=0
    

    
　

    

      　(1-a)Σ
       ∞
      ak
      = 1-a∞　　　　　　　　　[1]-a*[1]　
    

    

       k=0
    

    

      
      = 1　　　( |a|<1 )
    

    
　

    

      　　　　Σ
       ∞
      ak
      = 1/(1-a)
    

    

       k=0
    

  


a を -x2 に置きかえると、さきほどの df(x)/dx の形になります。


Σ(-x2)k=1/(1+x2)　　( |-x2|<1 )
つまり、
df(x)/dx=Σ(-x2)k
        =1-x2+x4-x6+x8-x10+・・・
        
積分して、
f(x)=x-x3/3+x5/5-x7/7+x9/9-x11/11+・・・+Const
　　=Σ(-1)kx(k*2+1)/(k*2+1) + Const
f(0)=0だから、Const=0

テーラー展開の定義式と微妙に違いますが、これもテーラー展開です。
微分すると階乗が相殺されたり、偶数項が0になったわけです。


面倒な tan の高階微分をせずに済みました。

なんだかインチキくさいですね。でも、インチキではないです。たぶん。
atan の級数展開は、後述の π に関連して相当に興味を惹く対象だったようで、
この方法は１７世紀には(なんと江戸時代の和算家にも)認識されていたそうです。

さて、この便利な形が使えるのは、|-x2|<1 の範囲に限られるのでした。
これは、-1<x<1 ですから、-π/4<f(x)<π/4 です。
これでは、tan の周期(π)の半分しかありません。もうひとひねり必要です。
tan の全周期をカバーするために、f(x) を π/2 ずらしてみます。


tan(f(x)-π/2) = sin(f(x)-π/2) / cos(f(x)-π/2)
　　　　　　　 = cos(f(x)) / (-sin(f(x)))
　　　　　　　 = -1/tan(f(x))
　　　　　　　 = -1/x
両辺のatanをとると、
f(x)-π/2 = atan(-1/x)
f(x) = atan(-1/x)+π/2


うまくいきました。
必要な範囲は、|-x2|>=1 ですから、ゼロ除算は避けられます。

余談ですが、この式からπを求める事が出来ます。

tan(π/4) = 1
π/4 = atan(1) 
 　 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 +・・・
π　 = 4*(1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 +・・・)

実際に計算機にかける場合は、符号周期でまとめて、

π　 = 4*( (1-1/3) + (1/5-1/7) + (1/9-1/11) +・・・ )
 　 = 4*( (3-1)/(1*3) + (7-5)/(5*7) + (11-9)/(9*11) +・・・ )
 　 = 8*( 1/(1*3) + 1/(5*7) + 1/(9*11) +・・・ )
とするとスピードが上がります。

πの計算例




計算項数110100100010000




級数値43.03.133.1403.1414




　たかがπの近似に、ずいぶんな計算量が必要ですね。
　それもそのはずで、n項目は4/(2*n-1)を足すか引くかしているので、1000項め近辺でも
　4/(2000-1)=0.002程度の振幅があり、小数点以下2桁目がかろうじて安定している程度です。

　ところが、真値を挟んで行ったり来たりしているので、1000項めと1001項めの平均をとると、
　唐突に精度が上がります。




計算項数110100100010000




級数値+次項/23.33.133.141543.14159213.141592650




　この場合の振幅は、4/(2*n-1)-4/(2*(n+1)-1) = 8/(4*n2-1) です。
　nがある程度大きい場合の近似誤差は、前者は 1/n、後者は 1/(n2) ですから、
　まさに次元が違います。

　πを求めるには、これよりも速い方法や明快な方法がたくさんあるようです。
　至近のトピックでは、1兆1429億531万8634桁目から314159265358が出現したなんてのも。
　(10の12乗桁の中で、特定の12桁の数が出現する確率を論ずるのは野暮というものです）
　また、求めたπにはこんな楽しみかたも。

　π関連の話は非常に面白いのですが、ここは整数のサイトなので、このへんで。


余談終了。

それでは実装です。


　ここより先の内容は現在制作中です。



３．平方根を造る　Newton-Raphson Method　〜　ニュートンラプソン法　〜


unsigned int ISqrt(unsigned int a)
{
    unsigned int x,y;
    if(a<=1) return a;
    y=0xc0000000;       //a桁判定初期値
    x=0x0000c000;       //sqrt(a)桁判定初期値
    while((y & a)==0){  //sqrt(a)<=x<sqrt(a)*2
        x>>=1;
        y>>=2;
    }
    x=(x+a/x)>>1;
    x=(x+a/x)>>1;
//  x=(x+a/x)>>1;
    return x;
}






４．平方根と三角関数の使用例　その１．　〜　ブツシュミ(物理シミュレーション)　〜


空気のようなスクリーンセーバ
　msaver.zip(12KB) for Win9X/2K/XP
もどる　次へ