2014.10.14

前へ
次へ
ホームページトップへ戻る

MYCPU80でCP/Mを!
超巨大基板の8080互換HCMOS・CPUでCP/Mを走らせてしまおうという、なんとも狂気なプロジェクトです!


[第52回]


●チェビシェフの逆展開式を求める(2)

前回はチェビシェフの11次の逆展開式を手計算で求めました。
今回はBorland C++を使ってプログラムで求めてみます。

[第45回]ではチェビシェフの多項式をC++を使って求めました。
そのときのプログラムの後ろに逆展開式を求めるプログラムを追加しました。

/// chebychev polynomials
/// 14/10/1 10/2 10/11 10/12
//
#include <stdio.h>
//
void main()
{
// i---> Ti(x)
// j---> x^j
// x[i][j]---> 'keisu' of x
        int x[21][21];
        int i,j;
        for (i=0;i<=20;i++){
                for(j=0;j<=20;j++){
                        x[i][j]=0;
                        }
                }
        x[0][0]=1;// T0
        x[1][1]=1;// T1
// keisan
        for(i=2;i<=20;i++){
                for(j=i;j>=2;j--){
                        x[i][j]=x[i-1][j-1]*2-x[i-2][j];
                        }
                x[i][1]=x[i-1][0]*2-x[i-2][1];//1次の項
                x[i][0]=-x[i-2][0];//定数項
                }
// print
        printf("T0(x)=1\n");
        printf("T1(x)=x\n");
        for(i=2;i<=20;i++){
                printf("T%d(x)=",i);//左辺
                for(j=20;j>1;j--){
                        if(x[i][j]>1&&i==j)printf("%dx^%d",x[i][j],j);//第1項(符号を省く)
                        else if(x[i][j]>1)printf("+%dx^%d",x[i][j],j);//符号
                        else if(x[i][j]<-1)printf("%dx^%d",x[i][j],j);//符号
                        }
                if(x[i][1]>1)printf("+%dx",x[i][j]);//1次の項&符号
                else if(x[i][1]<-1)printf("%dx",x[i][j]);//1次の項&符号
                if(x[i][0]==0)printf("\n");//定数項なし
                else if(x[i][0]==1)printf("+1\n");//定数項
                else if(x[i][0]==-1)printf("-1\n");//定数項
                }
//
// x^n=f(Tn)
//
        int jj;
        float a,b;
        float T[21][21];
        for (i=0;i<=20;i++){
                for(j=0;j<=20;j++){
                        T[i][j]=0;
                        }
                }
        T[0][0]=1;// x0
        T[1][1]=1;// x1
        for(i=2;i<=20;i++){
                a=x[i][i];//除数
                T[i][i]=1/a;//第1項の係数
                for(j=i-1;j>=0;j--){
                        if(x[i][j]!=0){
                                b=-x[i][j]/a;//第n項の最初の係数
                                for(jj=j;jj>=0;jj--){
                                        if(x[i][jj]!=0)T[i][jj]=T[i][jj]+T[j][jj]*b;
                                        }
                                }
                        }
                }
// print
        int n,m;
        printf("1=T0\n");
        printf("x=T1\n");
        for(i=2;i<=20;i++){
                printf("x^%d=",i);//左辺
                n=1/T[i][i];
                printf("1/%d(T%d",n,i);//第1項
                for(j=i-1;j>=0;j--){
                                if(T[i][j]>0){
                                        b=T[i][j]/T[i][i];
                                        m=b;
                                        if(m==1)printf("+T%d",j);
                                        else printf("+%dT%d",m,j);
                                        }
                                }
                printf(")\n");
                }
return;
}
//

// x^n=f(Tn)
から後ろが今回追加した部分です。
今回もTnの係数を格納するために2次元の配列Tを使いますが、今回は分数の計算をしますからfloat型の配列を使います。
計算では分数は小数として扱われます。

式の係数を求めるための計算は
for(i=2;i<=20;i++){
の{}の中のたった7行のみです。

なんでこうなるの?
と疑問に思われるかも知れませんが、前回説明しました手計算での方法をそのまま2次から順にやっていくと、こうなります。
こういうのはパズルみたいなもんです。

i,jに実際の値をあてはめて、上の式の配列がどこの係数値を示しているか、またどこの係数としてどういう値が入るか紙とエンピツで検算してみてください。
「おお、そういうことをやっているのか!」と納得されることと思います。
実際に紙に書いてみればわかりますので、解説はいたしません。
どうぞ皆様ご自身にてお確かめください。

計算は小数で行なわれますが、結果をそのまま小数で表示しますと、ピンときません。
実際の結果は分数に直すことができる値になります。
そこで結果の表示は分数に変換して表示するようにしました。

下が実行結果です。
前半は[第45回]と同じチェビシェフの多項式です。
後半はそれを元に計算して求めた逆展開式です。

T0(x)=1
T1(x)=x
T2(x)=2x^2-1
T3(x)=4x^3-3x
T4(x)=8x^4-8x^2+1
T5(x)=16x^5-20x^3+5x
T6(x)=32x^6-48x^4+18x^2-1
T7(x)=64x^7-112x^5+56x^3-7x
T8(x)=128x^8-256x^6+160x^4-32x^2+1
T9(x)=256x^9-576x^7+432x^5-120x^3+9x
T10(x)=512x^10-1280x^8+1120x^6-400x^4+50x^2-1
T11(x)=1024x^11-2816x^9+2816x^7-1232x^5+220x^3-11x
T12(x)=2048x^12-6144x^10+6912x^8-3584x^6+840x^4-72x^2+1
T13(x)=4096x^13-13312x^11+16640x^9-9984x^7+2912x^5-364x^3+13x
T14(x)=8192x^14-28672x^12+39424x^10-26880x^8+9408x^6-1568x^4+98x^2-1
T15(x)=16384x^15-61440x^13+92160x^11-70400x^9+28800x^7-6048x^5+560x^3-15x
T16(x)=32768x^16-131072x^14+212992x^12-180224x^10+84480x^8-21504x^6+2688x^4-128x^2+1
T17(x)=65536x^17-278528x^15+487424x^13-452608x^11+239360x^9-71808x^7+11424x^5-816x^3+17x
T18(x)=131072x^18-589824x^16+1105920x^14-1118208x^12+658944x^10-228096x^8+44352x^6-4320x^4+162x^2-1
T19(x)=262144x^19-1245184x^17+2490368x^15-2723840x^13+1770496x^11-695552x^9+160512x^7-20064x^5+1140x^3-19x
T20(x)=524288x^20-2621440x^18+5570560x^16-6553600x^14+4659200x^12-2050048x^10+549120x^8-84480x^6+6600x^4-200x^2+1
1=T0
x=T1
x^2=1/2(T2+T0)
x^3=1/4(T3+3T1)
x^4=1/8(T4+4T2+3T0)
x^5=1/16(T5+5T3+10T1)
x^6=1/32(T6+6T4+15T2+10T0)
x^7=1/64(T7+7T5+21T3+35T1)
x^8=1/128(T8+8T6+28T4+56T2+35T0)
x^9=1/256(T9+9T7+36T5+84T3+126T1)
x^10=1/512(T10+10T8+45T6+120T4+210T2+126T0)
x^11=1/1024(T11+11T9+55T7+165T5+330T3+462T1)
x^12=1/2048(T12+12T10+66T8+220T6+495T4+792T2+462T0)
x^13=1/4096(T13+13T11+78T9+286T7+715T5+1287T3+1716T1)
x^14=1/8192(T14+14T12+91T10+364T8+1001T6+2002T4+3003T2+1716T0)
x^15=1/16384(T15+15T13+105T11+455T9+1365T7+3003T5+5005T3+6435T1)
x^16=1/32768(T16+16T14+120T12+560T10+1820T8+4368T6+8008T4+11440T2+6435T0)
x^17=1/65536(T17+17T15+136T13+680T11+2380T9+6188T7+12376T5+19448T3+24310T1)
x^18=1/131072(T18+18T16+153T14+816T12+3060T10+8568T8+18564T6+31824T4+43758T2+24310T0)
x^19=1/262144(T19+19T17+171T15+969T13+3876T11+11628T9+27132T7+50388T5+75582T3+92378T1)
x^20=1/524288(T20+20T18+190T16+1140T14+4845T12+15504T10+38760T8+77520T6+125970T4+167960T2+92378T0)

の式までは「数値計算」(コロナ社)に載っているものと同じであることを確認しました。
またx11は前回手計算で求めたものと同じになっています。
ですので多分x10も正しく求められているものと推測できます。
問題はそれ以降の式が果たして正しいのかどうなのかということです。

11まで正しく求められていることからすると、プログラムのアルゴリズムはこれで間違いないようです。
気になりますのは途中の計算を小数で行なっている点です。
最後のほうになるとかなり大きな数になっていますから、ひょっとすると途中の計算の過程で誤差が発生しているのではないかとちょっと気になります。
途中で発生しているかもしれない誤差を考慮しないで、m=b; なんて簡単に実数を整数にしてしまっていますから、ここでひょっとすると端数が切り捨てられているかも知れないと少し心配になりました。

しかしこんな大きな値の計算を手計算で検算するのはとてもたまりません。
第一それじゃあプログラムで求めた意味がないじゃあありませんか。
なんとか余り面倒な計算をしないで、上で得られた結果が正しいかどうかを確かめる方法はないものか。
しばらく数字をながめていましたら。
おお。
できるじゃありませんか。

皆様もしばらくの間、上の数列をじっくりと眺めてみてください。
何かにお気付きになりませんでしょうか?

MYCPU80でCP/Mを![第52回]
2014.10.14upload

前へ
次へ
ホームページトップへ戻る