級数展開する上で、やはり２項定理をスキップする訳にいかないと考え直して、それについて触れることにした。
浅くしか理解していないが、とりあえず分かったつもりのレベルまで書いてみたい。
関連のrubyプログラムを作ったので、その解説も記す。分量的にはこの方が大きくなるかもしれない。
プログラム本体と関連のサンプルプログラムをbinomer.zipという圧縮ファイルに同梱した。
なお、2015/11/12の公開以降にメソッドの追加と修正を行った。それについては２項定理とそのプログラミング処理・追加編を参照されたい。
また、2015/12/16の公開後、級数に具体的な値を代入して計算するための calcメソッドを改良した。引数としてBigDecimalオブジェクトを与えると、計算結果をBigDecimalで返すようにした。Floatよりも精度の高い値を得ることができる。
２項定理の基本部分は簡単だ。
(a+x)^n を展開するとどうなるかをみる。まず、nが2, 3, 4の時にどうなるかを示す。
(a+x)^2 = aa + 2ax + xx → 係数の合計1+2+1=2^2=4
(a+x)^3 = aaa + 3aax + 3axx + xxx → 係数の合計1+3+3+1=2^3=8
(a+x)^4 = aaaa + 4aaax + 6aaxx + 4axxx + xxxx → 係数の合計1+4+6+4+1=2^4=16
４乗の場合について '+' で区切られた各項をみると、箱を４個用意して、箱にa, xのどちらかを入れる形になっている。
４個全部にaを入れるやり方は１通り、aを３個・xを１個入れるやり方が４通り、aを２個・xを２個入れるやり方は６通り、aを１個・xを３個入れるやり方が４通り、４個全部にxを入れるやり方は１通りということになる。
４個の箱全部にaを入れた状態からスタートして、xをどのように入れ込むか、どのように増やしていくかの試みといえる。
何通りであるかの数は、４個から何も選ばないケース、一つだけ選ぶケース、二つ選ぶケース、三つ選ぶケース、四つ選ぶケースの組み合わせの数になる。こうした組み合わせの数を２項係数あるいはncrとかnckというようだ。
ncrは、２乗の場合なら 1, 2, 1 となり、３乗だと 1, 3, 3, 1 また、４乗の場合は 1, 4, 6, 4, 1 のように上って下りる形になっている。
このような数を横に並べ、それを上から順に積み上げたものをパスカルの三角形というらしい。
一番上には１を二つ横に配置し、２段目が 1, 2, 1、３段目は 1, 3, 3, 1 などと配置する。
上下がちょっとずつずれるので、真下とか真上には数がなく、上をみても下をみても、斜めの位置に数が置かれている。
各段の左端と右端の１は別として、ある位置の数は、その斜め左上と斜め右上の数を足したものになっている。この規則性を利用すれば、ncrの算出が楽になる。
(a+x)^n を一般化して書くと、下のような式になる。
(a+x)^n =
a^n +
n/1 * a^(n-1) * x +
n*(n-1) / (2*1) * a^(n-2) * x^2 +
n*(n-1)*(n-2) / (3*2*1) * a^(n-3) * x^3 +
n*(n-1)*(n-2)*(n-3) / (4*3*2*1) * a^(n-4) * x^4 + ……
n*(n-1)*……*(n-r+1) / (r*……*2*1) * a^(n-r) * x^r + ……
n! / n! * x^n
最後の n! は、nの階乗の意味である。「n! / n!」は１なので、最後の項は x^n とだけ書いてもいいのだが、式の展開の流れからn!の階乗を記した。
途中に出てくる (r*……*2*1) は、rの階乗(r!)のことである。
上は (a+x)^n の展開だが、(a-x)^n の場合はどうなるだろうか。
いうまでもなく、x^r の項のrが奇数ならマイナス、偶数だとプラスの符号をつけることになる。
機械的な処理を念頭に置くなら、ncrを計算する時の分母の r! について、全要素をマイナスにしてやればいい。
つまり (r*……*2*1) → (-r*……*-2*-1) とする。階乗してから最後に -1 を掛ける訳ではないので注意。
結局、(a-x)^n の場合のncrを取り出すと下のようになる。
n/-1
n*(n-1) / (-2*-1)
n*(n-1)*(n-2) / (-3*-2*-1)
n*(n-1)*(n-2)*(n-3) / (-4*-3*-2*-1)
n*(n-1)*……*(n-r+1) / (-r*……*-2*-1)
(a-x)^n は、三角関数にかかわる sqrt(1-x^2) = (1-x^2)^(1/2) の展開に関連して用いることになる。
この辺から私には「？」が出るようになる。
(a+x)^n のnが自然数なら簡単だ。
n=0 の時も、どんな数も０乗すれば１になると割り切ってしまえばそれで済む。
(a+x)^-1 は、(a+x)^1 と掛け合わせた時に１になるのだから 1/(a+x) である。それくらいは分かる。
ただ、(a+x)^-1 を展開するのに、２項定理をそのまま用いることができる、というのが感覚的にどうにも腑に落ちない。
(a+x)^n のnが負の整数の場合も２項定理が成立することの証明はともかくとして、具体的な展開を行う際に留意する点について書いてみたい。
nが自然数、たとえば４の場合でいうと、ncrを計算するのに r: 1, 2, 3, 4 までで終わりになる。ncrの分母の部分は、その時々のrの階乗である。
具体的には 1!=1, 2!=2, 3!=6, 4!=24 となる。
ncrの分子の方は、4, 4*(4-1)=12, 4*(4-1)*(4-2)=24, 4*(4-1)*(4-2)*(4-3)=24 である。
ここで r が４で終わりにならず、更に１ずつ増え続けるケースを考えてみよう。
r=5 の場合、分母は 5!=120 である。
一方、分子の方は 4*(4-1)*(4-2)*(4-3)*(4-4) = 0 となる。
n=4, r=5 の時に分子が０になり、これ以降、rがどんどん増え続けても分子は０のままである。したがって、ncrは０だ。
一般に r = n+1 になった段階で、ncrが０になり、その後、rがどんな値になっても、その項は無視してかまわないことになる。
さて、n=-1 の場合はどうだろうか。rが１からスタートして１ずつ増えるとすれば、永遠に r = n+1 になることはない。つまり、(a+x)^n のnが負の整数の時は、無限級数になることを予感させる。
実際、1/(1+x) というのは 1 - x + x^2 - x^3 + x^4 …… という無限級数になる。この辺の事情については双曲線と対数そして積分で触れたが、２項定理による展開でも同じ結果を得る。
ともあれ、(a+x)^-1 の展開式をみてみよう。
(a+x)^-1 = a^-1 + -a^-2x + a^-3x^2 + -a^-4x^3 + a^-5x^4 + -a^-6x^5 + ……
上に出てくる a^-1 というのは 1/a のことであり、a^-2 は 1/a^2、a^-3 は 1/a^3 のことである。
aが１より大きければ、この部分は次第に小さくなる。しかし、aが１より小さい場合は、どんどん大きくなっていく。
xの方については、言うまでもなくその逆だ。xが１より大きければどんどん大きくなり、１より小さければどんどん小さくなる。
このように書くと、一見、対象性がちゃんと備わっているようにみえるが、なんだか奇妙だ。
たとえば、a=2, x=0.5 の場合にどうなるだろうか。
(2+0.5)^-1 = 1/2.5 = 4/10 = 0.4 である。
これは、aとxを入れ替えても同じはずだ。しかし、２項定理によって得られた級数を用いて計算すると、aとxを入れ替えると結果が大きく違ってくる。
級数を無限にたどるのは無理なので、第10項まで計算すると下のようになる。
a=2, x=0.5: 0.40000009536743164
a=0.5, x=2: 1677722.0
１番目の方は、第10項でやめずに無限に続けていけば0.4になりそうだ。
しかし、２番目の方は、無限に計算を続けてもどうにもならない。
足し算の交換法則という基本的なルールが保たれないものを「妥当である」とみなしていいのだろうか。
いろいろな解説を読んでも「それでいいのだ」というのが結論のようだ。だとすると、「まあ、いいか」という気になってくる。我ながら情けなくはあるが。
展開式が実際に機能するのは、xの絶対値がaの絶対値より小さい場合のようである。このような制約条件がなぜ生まれてくるのか、ちゃんとした根拠・証明があるのだろうが、私にはよく分からないので先に進むことにする。
ちなみに、指数nが自然数であれば、aとxを入れ替えても同じ結果になる。具体的な数値を当てはめるまでもなく、展開式をみれば明らかだ。対象性がきちんと保たれている。
そして、展開式が無限級数になることもない。ncrも、上ったら下りるという対象性が保持されている。
無限級数になった途端、素朴な対象性がどこかへいってしまう。ひたすら上り続ける(あるいは下り続ける)無限の一方向性がなせるイタズラだろうか。
指数nが負の整数の場合だけでなく、有理数の時も２項定理が使えるという。ニュートンは、このことを知っており、実際に用いたらしい。
実は、指数nは、有理数にとどまらず実数・複素数にも拡張することができるという。私にはとても追いつけないので、そこまでは触れないことにする。
指数nを1/2にするというのは、平方根を求めることに相当する。1/3なら立方根(３乗根)だ。２項定理によって、それらを算出する無限級数を得ることができることになる。なんだか魔法の杖を授けてもらった気分だ。ある意味、ペテンにかけられている気分といってもいい。
指数nを有理数にする場合も、２項定理の適用方法は同じだ。これまで述べてきた計算方法をそのまま使えばいい。
指数nが1/2の場合、級数の各項のncrを計算するとき、相変わらず r は 1, 2, 3, 4, …… と自然数の値をとる。したがって、ncrの分母は r!(rの階乗)である。
ncrの分子の方は、1/2*(1/2-1)*(1/2-2)*……*(1/2-r+1) となる。
参考まで (1+x)^(1/2) = sqrt(1+x) を２項定理で展開した式を掲げてみる。
(1+x)^(1/2) = 1 + 1/2x - 1/8x^2 + 1/16x^3 - 5/128x^4 + 7/256x^5 + ……
上の展開式は、平方根算出のための級数展開と近似値計算のところで掲げた展開式と同じである。
そこでは、２項定理は用いず平方根の展開式を「a + bx + cx^2 + dx^3 + ex^4 + ……」と仮定し、その上で a, b, c, d, e …… を求めた。
理解力のある人なら、「結局、両者は同じものなのだ」と見抜くのだろうが、私には不思議な感覚が残ったままである。
２項定理を扱うのに便利なように、rubyというプログラミング言語を用いて Binomer というクラスをつくった。
展開するのは、常に (a+x)^n であるという前提に立つ。
サンプルプログラムを掲げながら、関連の解説を加えてみる。
指数nが２の場合に、その展開式を出力するプログラムは下のとおり。
-------- bnr01.rb ここから
require "binomer"
bnr = Binomer.new(2)
puts bnr.cax(:string)
-------- bnr01.rb ここまで
上の３行を実行すると、下の文字が出力される。
a^2 + 2ax + x^2
Binomerクラスのオブジェクト(上ではbnr)を生成するとき、(a+x)^n の指数nを指定する。
こうして生成されたbnrについては、それを展開式用の文字に変換するメソッド cax などいくつかのメソッドが用意されている。
cax は、係数(coefficient), a, x の三つを示しているつもり。
これに対して cx というメソッドもある。こちらは a=1 という仮定に立って、係数とxだけを得る時に使う。
cax, cx の両方とも、各項を a^2 とか 2ax などの文字にしたものを配列に記録して返す。配列の各要素は、それぞれの「項」を文字にしたもの。
cax(:string), cx(:string) のように引数として :string を与えると、配列でなく文字列を返す。展開式として表示する場合は、この方が便利。
引数を :string5 とすれば、第５番目までの「項」が出力される。:string8 なら８番目まで。
指数nが-1の場合のサンプルを掲げてみよう。第５項まで出力するサンプル。
-------- bnr02.rb ここから
require "binomer"
bnr = Binomer.new(-1, 5)
puts bnr.cx(:string)
-------- bnr02.rb ここまで
上のプログラムの出力結果は下のとおり。
1 - x + x^2 - x^3 + x^4
上のプログラムで「Binomer.new(-1, 5)」として第２引数まで指定している。
第１引数は ncrのnであり、第２引数は ncrのrである。
第２引数のrを省略すると、nが自然数ならrがnと同じ値になり、そうでなければ r=10 となる。
最初のプログラム bnr01.rb では「Binomer.new(2)」としてrを省略した。なので、rはnと同じ２とみなされる。
展開式の a, x に実際に値を代入するプログラムも掲げてみよう。calc() というメソッドを使う。
sqrt(1/2) を求めて、その逆数をとって sqrt(2) つまり２の平方根を得る。
sqrt(1/2) = (1 - 1/2)^(1/2) を２項定理で展開して計算する。
本来は無限級数だが、第50項まで計算する。
-------- bnr03.rb ここから
require "binomer"
bnr = Binomer.new(Rational(1,2), 50) # (a+x)^(1/2)
val = bnr.calc(1, -0.5) # val = (1 - 1/2)^(1/2)
puts 1.0 / val
-------- bnr03.rb ここまで
Rational(1,2) は 1/2 の意味。0.5 と書いても同じ結果になる。
Binomerは、ncr, aの指数, xの指数にかんして、Rational(有理数型)として管理している。0.5を与えても内部で Rational(1,2) に変換される。
Binomerは、複素数には対応しないので注意。
実際に計算を行う calc() メソッドは、浮動小数点数(floatまたはBigDecimal)の値を返す。
与えられた引数がBigDecimalであればBigDecimalを返し、それ以外の場合はFloatを返す。
引数は、基本的に二つ与えるが、一つだけ指定した時は、それがxに代入される。そして、aには１が代入される。
なので、bnr03.rb に出てくる「calc(1, -0.5)」は、第１引数を省略して「calc(-0.5)」と書いてもよい。
参考まで bnr03.rb をBigDecimal用に書き換えたものを掲げてみると下のとおり。
-------- bnr03b.rb ここから
require "binomer"
bnr = Binomer.new(Rational(1,2), 400)
dlen = 100
rval = bnr.calc(BigDecimal.new("-0.5", dlen))
bin_sqrt2 = rval.power(-1, dlen)
sqrt2 = BigDecimal.new("2").sqrt(dlen)
print Bnm::dcmp(bin_sqrt2, sqrt2)
-------- bnr03b.rb ここまで
上に出てくる dlen は、数値の最大有効桁数(整数部と小数部を合わせた桁数)。BigDecimalの計算で使っている。
この dlen の指定を省略するとデフォルトの値になる。デフォルト値はパソコンやＯＳによって違うかもしれないが、40くらいではないかと思われる。
上のプログラムでは、２項定理で計算した√２のほか、BigDecimalのsqrtメソッドで計算した√２も算出して、くい違いがないかどうか比較している。
Bnm::dcmp() は、その比較のためのメソッド。比較の結果を文字列で返す。
２項定理での展開を300項までにすると、95桁のところでくい違いが出る。そのため400項までとした。
アークサインは、三角関数サインの逆関数である。ここでは、これを arcsin と表記する。
sin(30度)=1/2 なので、arcsin(1/2)=30度である。そして、30度を６倍すれば180度になる。180度は、ラジアンで表現すると円周率 3.14 だ。
arcsinの展開式が分かれば、そのxに 1/2 つまり 0.5 を代入して６倍することによって、円周率が算出される。それを Binomer を利用して試してみる。
説明は省くが、sqrt(1-x^2) の展開式にcircleという名前をつけたとすると、arcsinは下の式で表現される。
arcsin(x) = 2 * ∫(circle) - x * circle
circleは、(1-x^2)^(1/2) なので、２項定理を用いて展開式を得ることができる。
これらのことを利用してプログラムしすると、下のようになる。まだ触れていなかったメソッドがあれこれ出てくるが、メソッドの役割は容易に想像できると思う。
-------- bnr04.rb ここから
require "binomer"
circle = Binomer.new(Rational(1,2), 50, -2)
area1 = circle.copy.integral.inner("c *= 2")
area2 = circle.copy.inner("x += 1")
arcsin = area1 - area2
puts arcsin.calc(0.5) * 6.0
-------- bnr04.rb ここまで
circle つまり (1-x^2)^(1/2) を得る際に、「Binomer.new(Rational(1,2), 50, -2)」と三つの引数を与えている。
第１引数は、(a+x)^n のn(つまりncrのn)、第２引数は、第何項まで展開するかの値(つまりncrのr)である。今回は第50項まで展開する。
第３引数は (a+x)^n のxにかかわる指定で、これを -1 にすれば (a-x)^n の展開になり、-2 にすれば (a-x^2)^n の展開になる。2 だと (a+x^2)^n である。
指定しなければ 1 とみなされる。つまり、(a+x)^n となる。
circle.copyは、circleのコピーを生成して返す。戻り値は、もちろん Binomerオブジェクト。
copyには dupという別名も割り当てている。
integralは、展開式(Binomerオブジェクト)のxについて、各項積分を施す。戻り値は、積分処理された新たなBinomerオブジェクト。自分自身に変更を加える訳ではない。
circleに対して直接integralを適用しても、circleの中身が変更されることはない。なので、circle.copy.integral のようにcopyを適用しなくてもいいのだが、念のためcopyを用いてみた。
なお、integral! を適用すると、自分自身に変更を加える。自分自身が各項積分を施したものになる。
integralには inte という別名もある。
また、differential(別名 diff)という微分処理を施すメソッドも用意されている。
innerメソッドについて説明する前に、Binomerオブジェクトの内部情報に触れておこう。
オブジェクトの中には @bary という名前の配列がある。これには各項の「係数(c)」 「aの指数(a)」 「xの指数(x)」の三つから構成される配列が複数記録されている。複数というのは、「項」の数だけという意味である。
第１項の係数(c)は @bary[0][0] に入っており、aの指数は @bary[0][1]、xの指数は @bary[0][2] に記録されている。
第２項は @bary[1]、第３項は @bary[2] という具合に記録されている。
innerメソッドは、各項の c, a, x を一律に変更する時に用いる。
「inner(“c *= 2”)」とすると、全部の「項」のc(係数)の値が２倍になる。
「inner(“x += 1”)」とした場合は、全部の「項」のxの指数の値が１だけプラスされる。
innerメソッドは、処理を施した新たなBinomerオブジェクトを返す。自分自身に変更を加えたい場合は inner! を適用する。
innerの引数(文字列)は、その中に出てくる変数名 c, a, x が該当の要素に置換され、その上で ruby の eval() で実行される。
変数 c, a, x が複数回出てきても、また、混在してもかまわない。三つとも有理数型(Rational)なので、c.denominator(分母), c.numerator(分子)のような記述が可能。
rubyの文法に従っていれば大丈夫だと思うが、複雑なものだとちゃんと実行されないかもしれない。こった処理は、後述の to_a() で配列を取り出し、処理を加えてから from_a() でそれを反映させるのが無難。
'-', '+' は、メソッドというより演算子のイメージだが、rubyではメソッドとして定義できるので、ここでは「メソッド」と記す。
bnr04.rb には「arcsin = area1 - area2」という式が出てくる。area1もarea2もBinomerオブジェクトだ。
この時の '-' メソッドは、xの指数が同じである「項」について、前者(area1)の係数から公社(area2)の係数を引き算して、新たなオブジェクト(arcsin)を生成する。
xの指数が共通しない「項」は、新たなオブジェクト(arcsin)にそのまま引き継がれる。そのままといっても、後者のarea2については、その係数を０から引き算した値(つまり -1 を掛け算した結果)を引き継ぐ。
xの指数が共通するか否かを問わず、aの指数は、前者(area1)のものがそのまま引き継がれる。なので、展開式に具体的な値を代入して計算する場合、aに１でない値を代入しても意味のない値が算出されるので注意されたい。
'+' メソッドの方は、'-' が引き算するのに対して、係数の足し算を行って新たなBinomerオブジェクトを生成する。
なお、'+', '-' の後にBinomerオブジェクトでなく定数がきた場合は、定数項についての足し算または引き算になる。定数項がない場合は、指定された定数が新たに定数項となる。
掛け算の '*'、割り算の '/' などについては２項定理とそのプログラミング処理・追加編で触れる。
蛇足だが、arcsinは、次の１行でその展開式を得ることもできる。
arcsin = Binomer.new(-Rational(1,2), 50, -2).integral
以下、サンプルプログラムは掲げないが、まだ言及していないメソッドや仕様について記す。
変数bnrにBinomerオブジェクトが代入されている場合、bnr.to_a は、内部情報の @bary のコピー(配列)を返す。
戻り値の配列には、係数、aの指数、xの指数が複数個(展開式の「項」の個数分)記録されている。以下では便宜上、このような配列をcax配列ということにする。
なお、to_a の a は、array(配列の意味)の頭文字のつもり。
変数bnrにBinomerオブジェクトが代入されており、aryにはcax配列が代入されている場合、bnr.from_a(ary) とすれば、aryを内部情報として取り込む。
従来の @bary が破棄されて、新たに ary のコピーが @bary に代入される。
from_a の戻り値は、更新されたBinomerオブジェクト自身。cax配列を返す訳ではないので注意。一応 from_a! という別名も割り当てている。仕様は同じ。
Binomer.new(……) としてBinomerオブジェクトを生成する場合、実は内部で expansionを呼び出している。expansionは展開の意味。
expansionの引数の与え方は、当然ながら Binomer.new() の場合と同じである。
expansionの戻り値は、Binomerオブジェクト自身。
一応 expansion! という別名も割り当てているが、仕様は全く同じ。
オブジェクトを生成する時は何も引数を指定せず、後から expansion を呼び出すことができる。たとえば下のように書ける。
bnr = Binomer.new
bnr.expansion(Rational(1,2), 50, -2)
expansionの第１引数としてcax配列を与えた時は、第２・第３引数があったとしても無視し、そのcax配列のコピーを @bary にセットする。つまり、from_a と同じ挙動をとる。
この仕様は、Binomer.new(……) の場合も同じ。
変数bnrにBinomerオブジェクトが代入されている場合、bnr.coefficient は、cax配列から係数だけを抜き出して、配列の形で返す。
その戻り値には、展開式の「項」の個数分の係数が記録されている。
Binomerクラスに関する説明は以上でおわり。
今回のBinomerでは複素数を扱えない。rubyには複素数を扱うためのComplexクラスが用意されているので、プログラム的に複素数まで拡張するのは難しくないかもしれない。
ただ、私が２項定理に複素数を適用することについてまだ理解できていないため有理数までとした。
〜 「２項定理とそのプログラミング処理」おわり 〜
Copyright (C) T. Yoshiizumi, 2015 All rights reserved.