プログラミング雑感  新JavaScript入門  JavaScript,Neo-Generation  掲示板  表紙
7. 無限級数クラス(2)  9. 因数分解(1) 
プログラミング雑感
Written 5/05/02
8. 指数型無限級数クラス
前回、最後にtanx の級数展開がベルヌイ数によって表されることに言及しベルヌイ数にを求めたところ、 B10までしか求まりませんでした。
ベルヌイ数とは、 
  xex/(ex - 1) = sum(Bnxn/n!, n=0,…,∞)
と表せるのでした。ここでsumというのは、 
  sum(ai i=p,…,q) = ap + ap+1 + … + aq
を表すことにします。Σ記号を使うと書きにくいのでこのようにします。
次のような指数型無限級数を考えると、 高次（？）のベルヌイ数を求めることができます。 
  f(x) = sum(anxn/n!, n=0,…,∞)
前のとは各項に、/n! がついているところが違います。
■指数型とは
なぜ指数型というと、ex は、 
  ex = sum(xn/n!, n=0,…,∞)
だから、 
  an = 1
と非常に簡単になるからです。
前のは ex を表そうとするとすぐにオーバーフローしてしまいましたが、 この an をデータとすればそういうことは起こらなくなります。
■ファイル
完成品が次のファイルです。 
  taylorE.h
Eはexponentialのeです。 前のよりデータ構造が簡単で、バグが少ないはずです。
■データ構造
プロトタイプは次のようにしました。 
  class TaylorE : public vector<fraction> {
  public:
      TaylorE() { resize(100); }
      TaylorE(fraction a, int l);
      TaylorE(int a, int l);
      TaylorE(string func, int maxlength);
      TaylorE(const TaylorE& t) {
          (vector<fraction>&)(*this) = t;
      }
      TaylorE& operator =(const TaylorE& t) {
          (vector<fraction>&)(*this) = t;
          return *this;
      }
      TaylorE& operator +=(const TaylorE& a);
      TaylorE& operator -=(const TaylorE& a);
      TaylorE& operator *=(const TaylorE& a);
      TaylorE& operator /=(const TaylorE& a);
      TaylorE& operator ^=(int a);
      TaylorE& operator ^=(fraction a);
      TaylorE& operator +=(fraction a);
      TaylorE& operator -=(fraction a);
      TaylorE& operator *=(fraction a);
      TaylorE& operator /=(fraction a);
      TaylorE& operator +=(int a) { operator +=((fraction)a); return *this; }
      TaylorE& operator -=(int a) { operator -=((fraction)a); return *this; }
      TaylorE differential() const;
      TaylorE integral() const;
      TaylorE inverse() const;
      int getMin() const;
  };
今回はテンプレートは使わずに、vector<fraction> の派生クラスとしています。 また、最初の項は定数項を表すことにしました。 そして、前後の0の項も0の値を持ちます。
すなわち、vector<fraction> が 
  { 0, 0, 1, 2, 0, 0 }
という配列なら、 
  x2/2! + 2x3/3!
を表します。
このような表現のため、パブリックな派生クラスとしました。
最後に、getMin() は最小次数を返します。上の例なら 2 を返します。
■四則演算
加減算は簡単になりますが、乗除算はやや難しくなります。 
  f(x) = sum(aixi/i!, i=0,…,l)
  g(x) = sum(bixi/i!, i=0,…,l)
  h(x) = sum(zixi/i!, i=0,…,l)
としたときに、 
  h = f * g
の各項の係数は、 
  zi/i! = sum(aj/j!bi-j/(i-j)!, j=0,…,i)
すなわち、 
  zi = sum(iCjajbi-j, j=0,…,i)
となってcombinationの計算をしなければなりません。 これがオーバーフローの原因になります。 
  comb.cpp
  -------------------------------------------------------------------
  #include "fraction.h"
  
  void main(void) {
      for(int i = 1; i < 100; i++) {
              fraction    c = 1;
          for(int j = 0; j <= i; j++) {
              try {
                  cout << "<SUB>" << i << "</SUB>C<SUB>"
                          << j << "</SUB> = " << c << endl;
                  c *= fraction(i - j, j + 1);
              }
              catch(int e) {
                  cout << "<SUB>" << i << "</SUB>C<SUB>"
                          << (j + 1) << "</SUB> : overflow." << endl;
                  return;
              }
          }
      }
  }
このプログラムは次々にiCjを求めます。 その結果が次です。 
  1C0   = 1
  ...
  67C0  = 1
  67C1  = 67
  67C2  = 2211
  ...
  67C29 = 7886597962249166160
  67C30 : overflow.
これを避けるために、各項の係数の式を、 combinationの対称性を考慮して次のようにします。 
  zi = sum(iCj(ajbi-j + ai-jbj), j=0,…,i/2-1) + iCi/2ai/2bi/2  (i : even)
  zi = sum(iCj(ajbi-j + ai-jbj), j=0,…,(i-1)/2)             (i : odd)
これを次のような漸化式を用いて計算すると、 必要のないオーバーフローを避けられることが多くなります。 
  pi/2 = ai/2bi/2                         (i : even)
  p(i-1)/2 = a(i-1)/2b(i+1)/2 + a(i+1)/2b(i-1)/2  (i : odd)
  pj = (i-j)/(j+1)pj+1 + ajbi-j + ai-jbj
  zi = p0
除算もこれと同じ計算を用います。基本的な式は、 
  h = f / g
  f = h * g
  ai = sum(iCjzjbi-j, j=0,…,i)
  zi = (ai - sum(iCjzjbi-j, j=0,…,i-1)) / b0
最後の式が漸化式になってます。
ただし、これはb0が0でないときで、 実際には最小次数が除式の方が大きいときに例外を投げるなどの処理が必要です。
■ベルヌイ数
これでベルヌイ数を求めることができます。 
  bernoulli.cpp
  -------------------------------------------------------------------
  #include "taylorE.h"
  
  void main(void) {
      TaylorE x("x");
      TaylorE e("exp");
      TaylorE f = x * e / (e - 1);
      for(int i = 0; i < (int)f.size(); i++)
          cout << "B<SUB>" << i << "</SUB> = " << f[i] << endl;
  }
結果は次のようです。 
  B0  = 1       B12 = -691/2730         B26 = 8553103/6
  B1  = 1/2     B14 = 7/6               B28 = -23749461029/870
  B2  = 1/6     B16 = -3617/510         B30 = 8615841276005/14322
  B4  = -1/30   B18 = 43867/798         B32 = -7709321041217/510
  B6  = 1/42    B20 = -174611/330       B34 = 2577687858367/6
  B8  = -1/30   B22 = 854513/138
  B10 = 5/66    B24 = -236364091/2730
値が0のものは省いてます。
■微分・積分
微分は単なるシフトになります。 
  f(x)  = sum(anxn/n!, n=0,…,∞)
  f'(x) = sum(nanxn-1/n!, n=1,…,∞)
        = sum(anxn-1/(n-1)!, n=1,…,∞)
        = sum(an+1xn/n!, n=0,…,∞)
積分も同様です。
■逆関数
逆関数も基本式を示しましょう。
ここで最小次数は1です。 
  x = f(g(x))
  f(x) = sum(zixi/i!, i=1,…,l)
  g(x) = sum(aixi/i!, i=1,…,l)
  hi(x) = fi(x) = sum(bijxj/j!, j=1,…,l)
  b1j = aj
  bij = sum(jCkbi-1kaj-k, k=1,…,j)
これを最初の式に代入して、 
  1 = a1z1
  0 = sum(bjizj/j!, j=1,…,i) (i>=2)
  a1 = 1/z1
  ai = -sum(bjizj/j!, j=2,…,i)/z1 (i>=2)
最後の2式が漸化式となっています。
最後の式も前と同様に高い項から計算していきます。 
  aii = biizi
  aji = aj+1i/(j+1)+bjizj
  ai = -a2i/2
あまりいい例を思いつかないですが、 tanx を求めておきましょう。 
  TaylorE x("x");
  TaylorE f = 1 / (1 + x * x);
  TaylorE g = f.integral();     //tan-1x
  TaylorE h = g.inverse();      //tanx
  cout << h << endl;
この結果は次のようです。 
  x/1!+2x^3/3!+16x^5/5!+272x^7/7!+7936x^9/9!+353792x^11/11!+22368256x^13/13!
当然ながらあまり高次まで求まりません。
■べき乗
べき乗はまともに微分してTaylor展開の公式を使って求めます。 
  f(x) = ga(x)
  f(n)(x) = hn(x)ga-n(x)
  hn+1(x) = h'n(x) + (a - n)hn(x)g'(x)
最後の漸化式を計算します。 
  g(0) = 1
を前提としているので、 
  f(n)(0) = hn(0)
となります。 
  TaylorE e = TaylorE("exp", 100);
  cout << (e ^ fraction(1, 2)) << endl;
と e1/2 を計算すると、 
  1+1/2x/1!+1/4x^2/2!+1/8x^3/3!+1/16x^4/4!+1/32x^5/5!+1/64x^6/6!+1/128x^7/7!
  +1/256x^8/8!+1/512x^9/9!+1/1024x^10/10!+1/2048x^11/11!+1/4096x^12/12!
  +1/8192x^13/13!+1/16384x^14/14!+1/32768x^15/15!+1/65536x^16/16!
  +1/131072x^17/17!+1/262144x^18/18!+1/524288x^19/19!+1/1048576x^20/20!
  +1/2097152x^21/21!+1/4194304x^22/22!+1/8388608x^23/23!+1/16777216x^24/24!
と、これもあまり高次まで求まりませんね。
■tanx
最後に tanx を求めておきましょう。 
  tan.cpp
  -------------------------------------------------------------------
  #include "taylorE.h"
  
  void main(void) {
      cout << TaylorE("sin") / TaylorE("cos") << endl;
  }
この結果は、 
  x/1!+2x^3/3!+16x^5/5!+272x^7/7!+7936x^9/9!+353792x^11/11!+22368256x^13/13!
  +1903757312x^15/15!+209865342976x^17/17!+29088885112832x^19/19!
  +4951498053124096x^21/21!+1015423886506852352x^23/23!
となりました。前回よりは多くの項が求まりましたが、 ベルヌイ数を用いたほうがより多くの項が表すことができます。 
  tan2.cpp
  -------------------------------------------------------------------
  #include "TaylorE.h"
  
  inline fraction abs(const fraction& a) {
    return a.getBunsi() < 0 ? -a : a;
  }
  
  void main(void) {
      TaylorE x("x");
      TaylorE e("exp");
      TaylorE f = x * e / (e - 1);
      cout << "x";
      for(int i = 2; i <= (int)f.size() / 2; i++)
          cout << "+" << "2<SUP>" << i * 2 << "</SUP>(2<SUP>"
                  << i * 2 << "</SUP>-1)" << abs(f[i*2])
                  << "x<SUP>" << i * 2 - 1 << "</SUP>/" << i * 2 << "!";
      cout << endl;
  }
この結果は、 
  x+24(24-1)1/30x3/4!+26(26-1)1/42x5/6!+28(28-1)1/30x7/8!+210(210-1)5/66x9/10!
  +212(212-1)691/2730x11/12!+214(214-1)7/6x13/14!+216(216-1)3617/510x15/16!
  +218(218-1)43867/798x17/18!+220(220-1)174611/330x19/20!+222(222-1)854513/138x21/22!
  +224(224-1)236364091/2730x23/24!+226(226-1)8553103/6x25/26!
  +228(228-1)-2020342747/870x27/28!+230(230-1)8615841276005/14322x29/30!
  +232(232-1)-145255103/510x31/32!+234(234-1)2577687858367/6x33/34!
  +236(236-1)18660653648619265/6831740x35/36!
となりました。
以上ですが、楽しんでいただけたでしょうか。
first, prev, next, exit