プログラミング雑感  新JavaScript入門  JavaScript,Neo-Generation  掲示板  表紙
6. 無限級数クラス(1)  8. 指数型無限級数クラス 
プログラミング雑感
Written 4/21/02
7. 無限級数クラス(2)
■四則演算
四則演算は普通に行えます。
プログラムとしては、乗除算より加減算の方が精度の処理の問題もあって、 少し複雑になっています。 定数の加減はちょっと特殊で、 定数は精度が無限にあるということにしているため別コードになっています。
■fractionクラス
これでだいたいの演算が行えるわけですが、 係数は有理数であるときに真に有用になります。 そこで 分数クラスを適当にでっち上げました。 こんな感じに使います。 
  fraction  a, b, c;
  a = fraction(1, 2);
  b = 2;
  c = a / b;
  cout << c << endl;    //1/4
64ビット整数が使えれば、それを分母分子にしています。 また、オーバーフローする場合は例外を投げるようにしています。
■tanx を求める
これで早くも当初の目的を達成することができます。 
  #include "fraction.h"
  #include "taylor.h"
  
  typedef Taylor<fraction>    series;
  
  void main(void) {
      series  a, b, c;
      a = series("sin");
      b = series("cos");
      c = a / b;
      cout << c << endl;
  }
次のような結果となりました。 
  x+1/3x^3+2/15x^5+17/315x^7+62/2835x^9+1382/155925x^11+21844/6081075x^13
これは確かに手で計算できないですね。
割と短い項で打ち切られていますが、 これはオーバーフローをするとその前の項で打ち切られるからです。
■微分・積分
微分・積分は differential ・ integral というメソッドで普通にできます。
積分定数は考えません。 また、x-1 の項があれば、その前で打ち切られます。 
  series  x("x");
  cout << x.integral() << endl; //x2/2
  cout << (x + x * x).differential() << endl;   //1 + 2x
  series  f("sin");
  cout << f.differential().differential() + f << endl;    //d2f/dx2 + f = 0
■逆関数
割り算も手計算ではたいへんでしたが、 逆関数を求めるというのはさらに手間がかかります。 inverse メソッドがこれを簡単に実現します。 ただし、アルゴリズムが簡単になるように、 最小次数は1という制限をつけています。
次の例では、sin-1x を求めています。 
  series  f("sin");
  series  g = f.inverse();
  cout << g << endl;
結果は次のようになりました。 
  x+1/6x^3+3/40x^5+5/112x^7+35/1152x^9+63/2816x^11+231/13312x^13
もっともこれは簡単に計算できて、 
  x + … + 2nCn/(4n(2n+1))x2n+1 + …
となります。こういうのを高校時代によくやったのを思い出します。
また、 
  (tan-1x)' = 1/(1+x2)
なので、ここから tanx を求めることができます。 
  series  x("x");
  series  f = 1 / (1 + x * x);
  series  g = f.integral();     //tan-1x
  series  h = g.inverse();      //tanx
  cout << h << endl;
この結果は次のようです。 
  x+1/3x^3+2/15x^5+17/315x^7+62/2835x^9+1382/155925x^11+21844/6081075x^13
  +929569/638512875x^15+6404582/10854718875x^17
単純に割り算をしたときより多くの項を得ることができました。
■べき乗
べき乗は整数だけでなく分数乗も行えるようにしています。
分数乗のときは (1 + x)a のTaylor展開の公式を用いています。 すなわち、この場合は、最初の項が1に限られています。
例えば、さきほどの sin の逆関数ですが、 
  (sin-1x)' = (1-x2)-1/2
を利用するともっと多くの項を求めることができます。 
  series  x("x");
  fraction    a(-1, 2);
  cout << ((1 - x * x) ^ a).integral() << endl;
次のようになりました。 
  x+1/6x^3+3/40x^5+5/112x^7+35/1152x^9+63/2816x^11+231/13312x^13+143/10240x^15
  +6435/557056x^17+12155/1245184x^19+46189/5505024x^21+88179/12058624x^23
  +676039/104857600x^25+1300075/226492416x^27+5014575/973078528x^29
  +9694845/2080374784x^31
次にもっと多くの項を求める方法を示します。
■関数の合成
関数の合成を行うメンバ関数 synthesize も用意しました。 
  f.synthesize(g)   //f(g(x))
これを用いて sin の逆関数を求めます。 
  series  x("x");
  fraction    b(-1, 2);
  series  f = (1 - x) ^ b;
  f = f.synthesize(x * x);
  cout << f.integral() << endl;
次のようになりました。 
  x+1/6x^3+3/40x^5+5/112x^7+35/1152x^9+63/2816x^11+231/13312x^13+143/10240x^15
  +6435/557056x^17+12155/1245184x^19+46189/5505024x^21+88179/12058624x^23
  +676039/104857600x^25+1300075/226492416x^27+5014575/973078528x^29
  +9694845/2080374784x^31+100180065/23622320128x^33+116680311/30064771072x^35
  +2268783825/635655159808x^37+1472719325/446676598784x^39
  +34461632205/11269994184704x^41+67282234305/23639499997184x^43
  +17534158031/6597069766656x^45+514589420475/206708186021888x^47
  +8061900920775/3448068464705536x^49+5267108601573/2392537302040576x^51
  +61989816618513/29836347531329536x^53+121683714103007/61924494876344320x^55
  +956086325095055/513410357520236544x^57
  +1879204156221315/1062849512059437056x^59
  +7391536347803839/4395513236313604096x^61
  +2077805148460987/1297036692682702848x^63
最後のほうは1行に1項になってしまいましたね。
■tanx
おととい、tanx の計算をしていたら、気が付いてしまいました。 ベルヌイ数が使えますね。
まず、tan を指数で表します（今でも高校で習うかな？）。 
  tanx = (eix - e-ix)/(eix + e-ix)/i
ベルヌイ数 Bn は次のように定義されます（符号で複数の定義があるらしい）。 
  xex/(ex - 1) = B0 + … + Bnxn/n! + …
これと 
  1/(ex + 1) = 1/(ex - 1) - 2/(e2x - 1)
を利用すると簡単に計算できて、 
  tanx = x + … + (-1)n+122n(22n - 1)B2nx2n-1/(2n)! + …
これ、あってますか？
ただし、ベルヌイ数を計算するのもけっこうたいへんです。 定義に沿って、 
  series  x("x");
  series  e("exp");
  series  f = x * e / (e - 1);
  int     b = 1;     //階乗
  for(int i = 0; i < (int)f.size(); i++) {
      if(i >= 2)
          b *= i;
      cout << "B" << i << " = " << f(i) * b << endl;
  }
計算すると、 
  B0 = 1
  B1 = 1/2
  B2 = 1/6
  B3 = 0
  B4 = -1/30
  B5 = 0
  B6 = 1/42
  B7 = 0
  B8 = -1/30
  B9 = 0
  B10 = 5/66
これを利用して、 
  series  x("x");
  series  e("exp");
  series  f = x * e / (e - 1);
  int     c = 1;
  series  tan = 0;
  for(int i = 1; i < (int)f.size(); i += 2) {
      c *= 4;
      if(i % 4 == 1)
          tan += f(i+1) * c * (c - 1) * (x ^ i);
      else
          tan -= f(i+1) * c * (c - 1) * (x ^ i);
  }
  cout << tan << endl;
と tanx を計算すると、 
  x-1/3x^3+2/15x^5-17/315x^7+62/2835x^9
となりました。項数が少ないですね。
気が向いたらもうちょっとなんとかしたいと思います。
first, prev, next, exit