超高次の多項式の原始関数を求めたいのですが

質問

f(x) = (n-x)(n-1-x)(n-2-x).....(n-m-x)
n: 大きな自然数(例えば1000000など)
m: n>mの大きな自然数(例えば100000など)

という多項式 f(x)の原始関数を高速に求めるアルゴリズムを考えて
います．

f(x)を具体的に展開してから原始関数を求めれば簡単だと思い，上記
の式を展開するプログラムを書いたのですが，組み合わせの計算を
する必要が生じて，mの値が大きな時に高速に計算できませんでした．

原始関数を直接導出しようと，いろいろ場合分けして考えてみたので
すが挫折しました．

アドバイス頂けませんでしょうか？

よろしくお願いします．

回答No.3

stomachman · ベストアンサー

10万次多項式だなんて、すんごいですね。プログラムを書いたと仰るんだから、本当に計算する必要があるのでしょうけれど、一体どういう応用で出て来るんでしょうか。fの原始関数をどう使うんでしょうか。興味があります。というのは…
　ガンマ関数の近似式を使わず、数式処理としての不定積分を考えている、という状況ですから、多項式の係数を浮動小数点表示で丸めて表したんではおそらく意味がない。きっと整数で保持するのでしょう。たとえばxに何か整数値を与えてf(x)を計算するとしますと、因子(n-k-x)の値がだいたい６桁ぐらいだとして、それを10万個掛け算するので、答はおおざっぱに言って100万桁ぐらい。64bit(=10進で約20桁)の5万倍です。答を入れるメモリだけで500Kbyteぐらい必要になる。また「高速で」と仰る意味がよく分からないけれど、64bitの足し算に1ナノ秒掛かるパソコンだと５万倍長の足し算ひとつが50マイクロ秒程度、これを10万回やったら数秒掛かる。掛け算はFFTなどを使って旨く工夫しても足し算の何十倍かぐらい掛かるんで、１回1ミリ秒。10万回の掛け算なら1分は掛かりそうです。パソコンだとこれだけ時間を掛けて、ようやくひとつのxについてのf(x)の値が計算できるわけです。もし10万倍速いスパコンが使えたとすると、1ミリ秒でf(x)が出る（理想的なパフォーマンスで使えた場合の話ですが）。
　一方、f(x)を展開した場合には、係数が10万個発生します。そのひとつひとつがざっと100万桁。つまり、単に多項式の係数を全部記憶するだけでおよそ百ギガバイト級のメモリが要る。そして、係数１個を計算するのに必要な掛け算の回数が千回程度というすばらしく旨い方法があったとしても、100万桁×10万個の係数が算出されねばならないんだから、パソコンなら何日がかり。ソレ以下にはならんでしょう。だから、どのみちスパコンがなくちゃ始まらない規模の話です。

　幾つか方策を考えてみましたが、どうも結局、単に展開してから積分するのが速そうだと思います。で、展開も最も愚直にやるのが一番良さそうです。すなわち多項式
f[k]=(n-x)(n-1-x)…(n-k-x)
を展開したものを
f[k]=Σ{j=0～k}A[k,j](x^j)
とするとき、
f[k+1]=((n-k+1)-x)f[k]
=(n-k+1)A[k,0]+Σ{j=1～k}((n-k+1)A[k,j]-A[k,j-1])(x^j)-A[k,k](x^(k+1))
従って、
A[1,0]=n
A[1,1]=-1
A[k+1,0]=(n-k+1)A[k,0]
A[k+1,j]=(n-k+1)A[k,j]-A[k,j-1]　(1≦j≦k)
A[k+1,k+1] = -A[k,k]
という漸化式によって
f[m]=Σ{j=0～m}A[m,j](x^j)
を作る。
　A[k+1,j]を計算するのに必要なのはA[k,j], A[k,j-1]だけですから、逆に言えばA[k,j]とA[k,j+1]の計算が終わったらA[k-1,j]は忘れて構わない。つまりメモリの同じ場所に上書きして行けるから、100万桁×10万個ぐらいの係数を記憶できれば良い。従って、百ギガバイトぐらいのメモリで足りそうです。
　この漸化式でA[m,j](1≦j≦m)を得るまでに、掛け算がだいたい10万×10万回のオーダーで発生しますが、乗数の一方はほんの数桁しかない小さい数に限られています。だから、足し算の何倍か程度の時間で計算できるのではないか。だとすると、自家用のスパコンをお持ちなら何時間、パソコンなら何年かで答が出る、という程度の規模に収まるでしょう。

　いやいや、そんな高精度は必要ない、という話なのかも知れません。ですが、展開して扱おうとすると、零点がほぼ整数値のところに来るようにするためには、係数の精度が極めて高くなくてはならない。
　一方、原始関数の用途によっては、うんと簡単な近似法が適用できるかもしれません。たとえば、x～n+m/2の近辺だけを見れば、f(x)もF(x)もsin関数とほとんど見分けがつかないでしょう。

補足

回答どうもありがとうございました．リプライ遅くなりまして申し訳ございませんでした．

ご提示頂いた漸化式，展開毎に10万回×10万回程度の掛け算が発生する点など非常にありがたいアドバイスを頂きました．連立方程式に持ち込んだ場合と比べてメモリーの消費量が少なくて済みそうというのは確かに大きな利点です．

詳細を検討させて頂こうと思います．

どうもありがとうございました．

回答No.2

Tacosan

現状どう計算しているのかわかりません (「組合せの計算」が必要なところがちょっと思いつかない) しやってみないとわかりませんが, ベタに展開するという方針も悪くないような気がします. これなら, 理論上は m の線形オーダーになるんじゃないかな.

補足

お返事ありがとうございました．Reply遅くなりまして申し訳ございません.

展開した場合，mの線形オーダーに収まるかもしれないという点，もしよろしければご教授頂けませんでしょうか？

(前後しますが)No.3で回答頂いた方から漸化式が作れるとのアドバイスを頂き，漸化式を用いて前向きに(xの次数を増やす方向で)展開していった場合，展開する回数はmステップで済むということが理解できました．今回扱っている多項式の形状の特徴なのですが，項の数が多い(m個？)ので，その計算にも各展開毎にm回計算が必要になるのではというアドバイスも頂きました．

もし多項式の係数が疎な場合(展開した際x^iの個数がmに比べて少ない場合)にはmの線形オーダーで収まるのでは？と理解しています．

もし今回のような密な多項式の場合でもmの線形オーダーで収まるアイディアご存知でしたらアドバイス頂ければ幸いです．

ご回答ありがとうございました．

回答No.1

noname#221368

まず「思いつき」である事を、お断りしておきます。次に、高次多項式を実際に扱われているので、桁落ち対策として、多倍長処理などはやっている、は前提です。

　f(x)の次数はわかっているので（L次とします）、原始関数は(L＋1)次です。(L＋1)次の多項式を未定定数を含む形で仮定し、その微分（式aとします）を出すのは一瞬だと思います。
　f(x)の零点もわかっているので、a(零点)＝0とし、これを未定定数に関する連立方程式と考えます。このとき条件数より、未定定数の数が一個多くなりますが、それは「定数項＝n(n-1)・・・(n-m)」で解決できます。
　そうすると後は、上記連立方程式をいかに高速に解くかです。メモリとスピードから、共役勾配法なんて手はあると思います。

　結果は、やってみないと？、ですが・・・。

補足

アドバイスどうもありがとうございます．

なるほど...確かにそうですね．ご指摘通り，連立方程式に持ち込めそう
です．

＃f(x)の式の形の良い面(?)を見てあげるべきでした

大きな連立方程式を計算なのですが，シンプルに手計算のように変数を
1つづつ消していったとしても，1つの変数を消すために必要な主な計算
は，連立方程式の式のかず分(＝変数のかず回)程度の代入になり，それ
を変数が１つになるまで繰り返すとしても，全体でｆ(x)の次数(m-1で
しょうか)の2乗程度のオーダーで済みそうな気がします．

私が組み合わせを計算して式を展開していた時から比べると劇的
な進歩です！

連立方程式の高速解法もご指摘いただいてありがとうございます．
早速勉強してみます．

どうもありがとうございました．