C言語で学ぶ大学の数学

C言語で学ぶ大学の数学 「Laurent多項式の加法」

この記事では、Laurent多項式の加法をC言語で実行する方法を、前回紹介した方法 A と方法 B それぞれの場合に説明します。

１．方法 Aで定義した１変数Laurent多項式の加法

まずはLaurent多項式を方法Aで定義した場合の加法を説明します。

例として

a(x) = x^(-3) - x^(-1) + 2 + x

b(x) = x^(-2) + x^(-1)

に対し、c(x)= a(x)+b(x)を計算します。

方法 Aをおさらいしておくと、まずa(x)とb(x)の最低次数と最高次数を格納する変数を用意しました。

int ldeg_a, ldeg_b, deg_a, deg_b;

そして次数を代入します。ただし最高次数、最低次数ともに０以上とし、ldeg_a ,ldeg_b は符号を除去します。

deg_a = 1;

ldeg_a = 3;

deg_b = 0;

ldeg_b = 2;

次に c(x) の次数を考えます。c(x) を直接計算すると

c(x) = x^(-3) + x^(-2) + 2 + x

となります。このとき、c(x)の最高次数は deg_a と deg_b の大きいほう、最低次数は ldeg_a と ldeg_b の大きいほうになっていると分かります。

一般に c(x) の次数は次のように決定できます。

int deg_c, ldeg_c;

if(deg_a > deg_b)　deg_c = deg_a;

else　deg_c = deg_b;

if(ldeg_a > ldeg_b)　ldeg_c = ldeg_a;

else　ldeg_c = ldeg_b;

次に a(x), b(x) ,c(x)の係数にあたる配列変数を定義します。

float a[deg_a +ldeg_a +1];

float b[deg_b +ldeg_b +1];

float c[deg_c +ldeg_c +1];

そして a(x), b(x) の係数を代入します。

for(i=0; i<=deg_a +ldeg_a; i++){

　printf("coefficient of x^%d of a(x)=" , i -ldeg_a);

　scanf("%f",&a[i]);

}

for(i=0; i<=deg_b +ldeg_b; i++){

　printf("coefficient of x^%d of b(x)=" , i -ldeg_b);

　scanf("%f",&b[i]);

}

などを使うと便利です。

（配列変数に係数をまとめて代入することもできます。その場合は「配列変数の定義」と「代入」を同時に行います。例えば a(x) を正規化した多項式 a(x) = 1 - x^2 + 2x^3 + x^4 を代入するには

float a[] = {1, 0, -1, 2, 1};

とするか、

float a[5] = {1, 0, -1, 2, 1};

とすればできます。しかし、変数 deg_a , ldeg_a を使って

int deg_a = 1;

int ldeg_a = 3;

float a[deg_a +ldeg_a +1] = {1, 0, -1, 2, 1};

とするとエラーが発生します。詳しくは分かりませんが、配列変数にまとめて代入するやり方では[]の中身に「main 関数の中で定義した変数」は使えないようです。

ですが main 関数の外に

#define deg_a 1

#define ldeg_a 3

と定義した場合は

float a[deg_a +ldeg_a +1] = {1, 0, -1, 2, 1};

とすれば正しく代入することができます。

次数の上限を固定するやり方であれば #define を使ってもいいですが、今回はややこしいので for 文を使うか

float a[] = {1, 0, -1, 2, 1};

float b[] = {1, 1, 0};

とすることとします。）

ここから加法を計算します。

問題となるのが、最低次数のぶんa[], b[], c[]の引数がずれていて、一般には

c[n] = a[n] + b[n]

のような関係が成り立たないことです。

また、後でソースコードを見た時、引数のずれへの対処が正しいかどうか簡単にチェックできた方が都合が良いです。

そこで、次数がマイナスの部分と、次数が０以上の部分を別々に計算します。

まずは次数が０以上の部分から加法を計算します。

あらかじめint型の変数 n を用意しておきます。

for(n=0; n<=deg_c; n++){

　c[ldeg_c +n] = 0;

　if(n<=deg_a)　c[ldeg_c +n] += a[ldeg_a +n];

　if(n<=deg_b)　c[ldeg_c +n] += b[ldeg_b +n];

}

ポイントは引数がいずれも「ldeg_a +n」のようになっていることです。

これは x^n の係数が a[ldeg_a +n] に格納されていることに起因します。

つまりnを 0 からdeg_a まで動かせば、a(x)の０以上の項を次数が低いものから順番に見ていることになります。

ldeg_aなどで引数がずれていますが、if 文の使い方はおおむね通常の１変数多項式の加法と同じです。

では次に、次数がマイナスの項を計算します。

for(n=1; n<=ldeg_c; n++){

　c[ldeg_c -n] = 0;

　if(n<=ldeg_a)　c[ldeg_c -n] += a[ldeg_a -n];

　if(n<=ldeg_b)　c[ldeg_c -n] += b[ldeg_b -n];

}

ポイントは、次数が-1の項からスタートして、最低次数の項を最後に計算するようにしていることです。

ifの中身の「n<=ldeg_a」は「ldeg_a -n≧0」と同値です。つまりa[]の引数が０未満のときにa[]を加えています。

このようにしてLaurent多項式の加法は計算できます。

なお、これらを一つにまとめると最低次数から計算することができます。

具体的には次のようになります。

for(n=-ldeg_c; n<=deg_c; n++){

　c[ldeg_c +n] = 0;

　if(n>=-ldeg_a && n<=deg_a) c[ldeg_c +n] += a[ldeg_a +n];

　if(n>=-ldeg_b && n<=deg_b) c[ldeg_c +n] += b[ldeg_b +n];

}

２．方法 B で定義した１変数Laurent多項式の加法

今度は、方法 B で Laurent多項式を定義した場合の加法について説明します。

例として方法 A と同じ多項式

a(x) = x^(-3) - x^(-1) + 2 + x

b(x) = x^(-2) + x^(-1)

に対し、c(x)= a(x)+b(x)を計算します。

まずは a(x), b(x) の次数をそのまま格納します。

int deg_a = 1;

int ldeg_a = -3;

int deg_b = -1;

int ldeg_b = -2;

次に c(x) の次数ですが、最高次数は deg_a とdeg_b のうち大きいほう、最低次数は ldeg_a とldeg_b のうち小さいほうになります。

int deg_c, ldeg_c;

if(deg_a > deg_b)　deg_c = deg_a;

else　deg_c = deg_b;

if(ldeg_a < ldeg_b)　ldeg_c = ldeg_a;

else　ldeg_c = ldeg_b;

そして a(x), b(x) の係数を代入します。正規化するとそれぞれ 1 - x^2 + 2x^3 + x^4,

1 + x になるので

float a[] = {1, 0, -1, 2, 1};

float b[] = {1, 1};

とします。b(x) は方法 Aのときと異なっていることに注意します。

それから c(x) の係数を格納するための配列変数を用意します。

float c[deg_c - ldeg_c +1];

ldeg_c の係数がマイナスになっていますが、deg_c ≧ ldeg_c という不等号が常に成り立つので特に問題ありません。

では、ここから加法を計算します。

ポイントはやはり a(x) と b(x) から同じ次数の項を取ってくることです。

a(x) の「x^n の係数」が正規化により「x^(n - ldeg_a)の係数」になっている、つまり「a[n - ldeg_a]」に格納されています。なのでだいたいの n に対し

c[n - ldeg_c] = a[n - ldeg_a] + b[n - ldeg_b]

という関係式が成り立っています。

これに引数に対する条件を反映させて計算すると

for(n=ldeg_c; n<=deg_c; n++){

　c[n - ldeg_c] = 0;

　if(n>=ldeg_a && n<=deg_a) c[n - ldeg_c] += a[n - ldeg_a];

　if(n>=ldeg_b && n<=deg_b) c[n - ldeg_c] += b[n - ldeg_b];

}

となります。よく見ると方法 A の最後のまとまりで ldeg_x を -ldeg_x (x は a, b, c)に置き換えたものになっています。

次の記事では Laurent多項式の乗法を計算する方法を説明します。