C言語で学ぶ大学の数学

C言語で学ぶ大学の数学 準備編「多項式に多項式を代入する」

この記事ではC言語で以下の計算をする方法を説明します。

・１変数多項式に１変数多項式を代入する

・次数が決められた多項式にユーザ関数を使って多項式を代入する

１．１変数多項式に１変数多項式を代入する

前回の記事では、多項式に数値を代入する方法を解説しました。その場合は多項式を「関数」とみなし、数値を代入することで関数のとる値を計算していたことになります。

今回は多項式に多項式を代入します。これは多項式を「関数」とみなし、「関数の合成」を計算していることになります。

具体例として、

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

b(x) = -1 + x^2

という多項式について、a(b(x)) を求めたいと思います。ちなみに正しい結果は

a(b(x)) = 2 + x^2 - 3x^4 + x^6

です。

まずは上の二つの多項式を配列変数に格納します。係数は今回は int 型にします。

int deg_a = 3;

int deg_b = 2;

int a[deg_a +1], b[deg_b +1];

a[0] = 1;

a[1] = -2;

a[2] = 0;

a[3] = 1;

b[0] = -1;

b[1] = 0;

b[2] = 1;

ここから計算に入ります。計算の順番は色々ありますが、前回の計算の順番に倣って「a[i]*(b(x))^i 」を i=0 から順に足していくこととします。

前回は「x^i」をpowx という変数に格納していました。今回は「(b(x))^i」をpowb[]という配列変数に格納します。i の最大値は a(x)の次数３なので、(b(x))^i の次数の最大値は６です。つまりpowb[] は６次多項式を想定して

int powb[deg_a *deg_b +1];

と定義します。

また(b(x))^0 = 1 に基づいて初期値を代入します。

for(i=0; i<=deg_a*deg_b; i++){

　if(i==0) powb[i] = 1;

　else powb[i] = 0;

}

多項式に数値を代入するプログラムのときは、a(x) の次数が上がるごとに powx に x を掛けていました。

同じように今回は a(x) の次数が上がるごとに「 powb[]の表す多項式」に 「b(x)」を掛けていきます。その部分のやり方は次のように計算します。

for(n=0; n<=deg_a*deg_b; n++){

　c[n] = 0;

　for(j=0; j<=deg_a*deg_b; j++){

　　if(n-j>=0 && n-j<=deg_b)　c[n] += powb[j]*b[n-j];

　}

}

for(n=0; n<=deg_a*deg_b; n++){

　powb[n] = c[n];

}

ポイントは、一度 c[] という変数に計算結果を代入して、全て計算し終えてから powb[] に代入していることです。そうしないと計算途中で powb[] の値が変わってしまい正確な値が得られません。

上で説明したことを組み合わせて、次のようなソースコードを書きます。

#include<stdio.h>

int main(){

　int deg_a = 3;

　int a[deg_a +1];

　a[0] = 1;

　a[1] = -2;

　a[2] = 0;

　a[3] = 1;

　int deg_b = 2;

　int b[deg_b +1];

　b[0] = -1;

　b[1] = 0;

　b[2] = 1;

　int i, j, n;

　int result[deg_a *deg_b +1], powp[deg_a *deg_b +1], c[deg_a *deg_b +1];

　for(i=0; i<=deg_a *deg_b; i++){

　　result[i] = 0;　　……(※１)

　　if(i==0)　powp[i] = 1;

　　else　powp[i] = 0;

　}

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

　　for(n=0; n<=deg_a *deg_b; n++)　result[n] += a[i]*powp[n];　　……(※２)

　　for(n=0; n<=deg_a *deg_b; n++){

　　　c[n] = 0;

　　　for(j=0; j<=deg_a *deg_b; j++){

　　　　if(n-j>=0 && n-j<=deg_b) c[n] += powp[j]*b[n-j];

　　　}

　　}

　　for(n=0; n<=deg_a *deg_b; n++){

　　　powp[n] = c[n];

　　}

　}

　for(i=0; i<=deg_a *deg_b; i++){

　　if(result[i] != 0){

　　　printf("%dx^%d ", result[i], i);　　……(※３)

　　　if(i < deg_a *deg_b)　printf("+ ");

　　}

　}

　printf("\n");

　return 0;

}

※１　result[] の初期化がどこかで必要なので、powp[] に初期値を代入するときに一緒に初期化しています。

※２　多項式に数値を代入する計算のときの「a[i]*x^i を足す」部分に対応しています。

※３　計算結果を表示させています。result[i] != 0 という条件は int 型で使用するものなので、係数を float 型にする場合は abs(result[i])>1e-5 などのようにします。また「+」の出力を次の行に分けていますが、どう出力するかは好みの問題なので自由にしてください。

２．次数が決められた多項式にユーザ関数を使って多項式を代入する

上で行った計算を、今度はユーザ関数を使って実行する方法を説明します。

配列変数を多用する都合上、ユーザ関数の返り値は void とします。

返り値 void は、通常は返り値が存在しない場合に使用します。例えば

void print_poly(float a[6+1]){

　int i;

　for(i=0; i<=6; i++){

　　if(abs(a[i]) > 1e-5){

　　　printf("%fx^%d ", a[i], i);

　　　if(i < 6)　printf("+ ");

　　}

　}

　printf("\n");

}

というユーザ関数は、float 型の変数を係数とする６次多項式を printf で表示する関数で、返り値（return の後の数値のこと）はありません。

一方で、配列変数の数値をまとめて変更するときなどにも使用されます。

例えば、「多項式 a(x)と b(x) の積を計算し、その値を a(x) の係数 a[] に代入する」ユーザ関数は

void prod_poly(float a[6+1], float b[6+1]){

　int n, i;

　float c[6+1];

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

　　c[n] = 0;

　　for(i=0; i<=6; i++){

　　　if(n-i>=0 && n-i<=6) c[n] += a[i]*b[n-i];

　　}

　}

　for(n=0; n<=6; n++)　a[n] = c[n];

}

のように定義されます。

さて、多項式に多項式を代入する計算ですが、上の二つのユーザ関数ともう一つのユーザ関数で計算できます。次のようにユーザ関数を定義します。

void subst_poly(float a[6 +1], float b[6 +1], float result[6 +1]){

　int i, n;

　float powb[6 +1];

　for(i=0; i<=6; i++){//result,powpの初期化

　　result[i] = 0;

　　if(i==0)　powb[i] = 1;

　　else powb[i] = 0;

　}

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

　　for(i=0; i<=6; i++){

　　　result[i] += a[n]*powb[i];

　　}

　prod_poly(powp, b);　//ここでユーザ関数 prod_poly を使用

　}

}

subst_poly の中身で prod_poly を使っているため、ソースコードには subst_poly より上に prod_poly を書く必要があります。

計算の面倒な部分を subst_poly に押し付けたので、main 関数は簡単になります。

#include<stdio.h>

#include<stdlib.h>

void prod_poly(float a[6+1], float b[6+1]){

　（省略）

}

void subst_poly(float a[6 +1], float b[6 +1], float result[6 +1]){

　（省略）

}

void print_poly(float a[6+1]){

　（省略）

}

int main(){

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

　float b[6 +1] = {-1, 0, 1, 0, 0, 0, 0};

　float result[6+1];

　subst_poly(a, b, result); //ここで a(b(x))を計算

　print_poly(result); //ここでa(b(x))を表示

　return 0;

}

気を付けることとしては、a[] と b[] を６次多項式とみなして格納しているところです。というのも、ユーザ関数の引数がいずれも６次多項式として定義しているからです。配列変数の型が全て float なのにも注意します。

また、ユーザ関数を使う時は、引数の[]を省略します。

一般の次数の多項式の場合は、６としているところをa(x) の次数と b(x) の次数の積に取ります。また subst_poly の中の３は a(x) の次数に置き換えます。