Re: [問題] 請問該如何寫成副程式

作者: PhysiAndMath (老師說要愛數學)   2015-03-29 22:19:43
※ 引述《syxuan ()》之銘言:
: 遇到一個問題
: (x-a1)(x-a2)(x-a3)(x-a4)...(x-an)
: 要找出方程式的某個次方的係數
: 下面是只有四項要找三項的迴圈
: for(i[1] = 3; i[1] <= 4; i[1]++) {
: for(i[2] = 2; i[2] <= (i[1]-1); i[2]++) {
: for(i[3] = 1; i[3] <= (i[2]-1); i[3]++){
: sum = sum + a[i[1]]*a[i[2]]*a[i[3]];
: printf("i1=%d, i2=%d, i3=%d, sum=%d\n", i[1], i[2], i[3], sum);
: }
: }
: }
: 不知道要怎麼用副程式的方式寫成可以有n項取m次方的係數
臨時想到的
不知道有沒有更好的方法。
Pn(x)=(x-a1)(x-a2)(x-a3)(x-a4)...(x-an)
第m次方的係數等於
1/(2πi)∮Pn(z)/z^(m+1) dz
我開玩笑的
以下認真
基本idea是這樣的
最高次數是m+n
要找出第n次項係數
如果能夠有效率的窮舉會乘出x的n次方的可能性
問題基本上應該算是解決了
如果只乘了n次x,表示總共選出了m個ai
考慮以下窮舉的演算法
(以下所有計數都從1開始)
1。
一開始有m個指標p,p[i]=&a[i]
2。
回傳所有指標所指的組合的乘積
3。
如果p[m]不指向a[n+m],則p[m]=p[m]+1
然後回(2。)
否則就檢查p[m-1]
4。
如果被檢查的p[j]+1 == p[j+1]而且j==1就停下來
如果j!=1但p[j]+1 == p[j+1]就檢查p[j-1] 回到(4。)
否則
p[j]=p[j]+1;
for(int i=j+1;i<=m;i++)
p[i]=p[i-1]+1;
回到(2。)
5。
把2。回傳的結果加總起來就是結果
觀察(2。)的回傳結果
第一次
a[1]…a[m-1]a[m]
第二次
a[1]…a[m-1]a[m+1]

第n次
a[1]…a[m-1]a[m+n]
第n+1次
a[1]…a[m-2]a[m]a[m+1]

第2n-1次
a[1]…a[m-2]a[m]a[n]
現在定義b[m][i]= Σa[j] ,j=i to n+m ,i= m to n+m
則前n次可以寫成a[1]…a[m-1]b[m][m]
繼續定義b[j-1][i] = Σa[k]b[j][k+1] , k=i to n+j-1 ,i=j-1 to n+j-1
發現a[1]…a[j-2]b[j-1][j-1]可以算完前面j-2個指標不動的情況下所有的組合
得到b[1][1]就是所要的答案
提供參考
作者: abcdcba (wade)   2015-03-30 01:00:00
電腦不會複變哈哈
作者: firejox (Tangent)   2015-03-30 17:57:00
用蒙地卡羅阿wwww
作者: PhysiAndMath (老師說要愛數學)   2015-03-30 23:30:00
痾…真要做也可以用gaussian quadrature

Links booklink

Contact Us: admin [ a t ] ucptt.com