C言語で学ぶ大学の数学

C言語で学ぶ大学の数学 「形式べき級数の逆元」

この記事では、定数項が０でない多項式の「乗法逆元」のような形式べき級数を導入し、その係数を次数が小さいほうから逐次的に計算する方法について解説します。

１．形式べき級数とは

多項式とは、環 R を係数とする単項式の有限個の和のことです。

それに対し、形式べき級数とは環 R を係数とする単項式の無限和のことです。

例えば、x^n の係数が全て 1 であるような形式べき級数

1 + x + x^2 + x^3 + …

などがあります。

形式べき級数のなす環を形式べき級数環といいます。形式べき級数は多項式を拡張した概念なので、その和は次数ごとに計算し、積はa(x)とb(x)それぞれの単項式から次数の和がn になる組み合わせを全て選んできて掛けたものを足し合わせます。

２．形式べき級数の逆元

形式べき級数 a(x) = a_0 + a_1*x + a_2*x^2 + a_3*x^3 + … の逆元

b(x) = b_0 + b_1*x + b_2*x^2 + b_3*x^3 + …

が存在するとき、b(x)の係数を計算することを考えます。

a(x)*b(x) = 1 の係数を０次の項から比較すると、

a_0*b_0 = 1

a_0*b_1 + a_1*b_0 = 0

a_0*b_2 + a_1*b_1 + a_2*b_0 = 0

a_0*b_3 + a_1*b_2 + a_2*b_1 + a_3*b_0 = 0

と続いていき、これを変形すると

b_0 = 1/a_0

b_1 = -a_1*b_0/a_0

b_2 = -(a_1*b_1 + a_2*b_0)/a_0

b_3 = -(a_1*b_2 + a_2*b_1 + a_3*b_0)/a_0

一般に

b_n = -(a_1*b_(n-1) + a_2*b_(n-2) + … + a_(n-1)*b_1 + a_n*b_0)/a_0

となり、{b_n}をn=0から逐次的に計算することができます。見ての通り、a(x)の逆元が存在するためにはa_0≠0でなければなりません。

上の計算を計算機で実行する際に注意することは、計算機に入力した a(x)の項の次数より大きい次数の b(x)の係数は計算できないという点です。

例えばa(x)の100次の項までを使ってb(x)を計算する場合、正確に計算できるのはb(x)の100次までの項で、b(x)の101次の項はa_101 の値を格納していないので計算できません。

しかしa(x) が多項式の場合、つまりあるNより大きい任意の自然数 n について a(x)の x^n の係数が 0 である場合には、b(x) の係数は次数がいくつであろうと原理的には計算することができます。

いずれにせよ、a(x)の有限個の項の係数しか計算機に代入することはできないため、以下では a(x) を単なる多項式とみてb(x) の項の計算をすることとします。

（注：a(x)の係数が次数 n を使った初等関数で表示できたり、n について周期関数になっていたり、漸化式で定義されるような場合は無限項の係数を代入しているのと同じことになります。そのような場合はユーザ関数を適切に定義すればb(x)の係数を任意次数まで計算できますが、ややこしいので今回は触れません）

３．C言語で多項式の逆元を計算する

前節で説明した b_n を次数の小さい項の係数を使って表示する方法を利用して、C言語で b_n を計算する方法について説明します。

ここでは a(x) の次数を５次以下とします。例として

a(x)=2 + x + x^3 + x^4 -4*x^5

の逆元を5次まで計算することとします。

#include<stdio.h>

#include<math.h>

#define DEG_A 5

void calcterm(float a[], float b[], int n){

　int i;

　b[n] = 0;

　for(i=1; i<=n; i++)

　　if(i<=DEG_A) b[n] += a[i]*b[n-i];

　b[n] = -b[n]/a[0];

}

int main(){

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

　int deg_b = 5;

　float b[deg_b +1];

　int n;

　b[0] = 1/a[0];

　for(n=1;n<=deg_b;n++) calcterm(a,b,n);

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

　　if(fabsf(b[n])>0.001){

　　　printf("%.3lfx^%d + ", b[n], n);

　　}

　}

　printf("\n");

　return 0;

}

#define DEG_A 5 は「DEG_A」という定数を定義していて、その値は５と定めています。

calcterm 関数はb[n]を計算する関数です。nは3番目の変数で指定します。返り値が void なので計算を終えるとb[n] に結果を格納します。

main 関数の中では、まずa(x)の係数をa[]に格納し、計算したい逆元の次数をdeg_b に格納します。

各b[n]の計算について、b[0]だけは calcterm 関数で計算できないので別でa[0]の逆数を代入しています。

n次まで全て計算し終えてから次のfor 文のところで計算結果を表示しています。ここで fabsf はfloat型の変数の絶対値を返す関数で、これを使うために math.h をインクルードしています。

printf 関数内の「%.3lf」というのは「float型またはdouble型の変数を小数点以下3桁まで表示する」ときに使います。表示したい桁数を変えるときは異なる数字にします。

a(x)が形式べき級数の一部ではなく単なる多項式である場合は、deg_b に代入する値だけを変えることで任意の次数まで逆元を計算することができます。

ちなみに計算結果が正しいかを確かめるには、b[n]を全て計算し終えた後で

float check;

int i;

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

　check = 0;

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

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

　}

　if(fabsf(check)>0.001) printf("%.3lfx^%d + ", check, n);

}

のようにしてa(x)とb(x)の積を計算し、これが１になることを調べれば良いです。

４．（おまけ）形式べき級数の無限和の解釈

無限和といえば数列や関数列の和を思い浮かべると思います。

例えば数列の和

1/2 + 1/2^2 + 1/2^3 + 1/2^4 + … + 1/2^n + …

はユークリッド空間の位相で 1 に収束します。これは「第n項までの部分和 S_n からなる数列{S_n}が１に収束する」という意味です。

また関数項級数

1 + x + x^2/2! + x^3/3! + x^4/4! + … + x^n/n! + …

は一様ノルムの定める位相で指数関数 exp(x) に収束します。これは「第n項までの部分和 f_n(x)からなる関数列{f_n} がexp(x)に収束する」という意味です。

同じように、形式べき級数も「第n項までの部分和をa_n(x)とすると、多項式の列{a_n(x)}が形式べき級数 a(x)に収束する」と考えたくなります。

実際にそのような解釈を与える形式べき級数環の位相があります。

Rを可換環、R[[x]]を形式べき級数環としたとき、x^nR[[x]]を０または最低次数がn以上の形式べき級数全体の集合とします。つまり

x^nR[[x]] = {x^n*f(x) | f(x)∈R[[x]]}

です。このとき、0∈R[[x]] の基本近傍系が{x^nR[[x]]}_(n=0,1,2,…)であり、f(x)の基底近傍系がそれをf(x)だけ平行移動したもの{f(x) + x^nR[[x]]}_(n=0,1,2,…)であるようなR[[x]]の位相が定義できます。ただし

f(x) + x^nR[[x]] = {f(x) + x^n*g(x) | g(x)∈R[[x]]}

です。

この位相に関して{a_n(x)}が形式べき級数a(x)に収束することは次のように確認できます。

Uをa(x)の任意の開近傍としたとき、mを十分大きく取ると a(x) の基本近傍系からUに含まれるa(x) + x^mR[[x]]が得られます。このとき、n≧mとするとa_n(x)は m 次以下の項がa(x)と等しくなります。つまりn≧mのときa_n(x)はa(x) + x^mR[[x]]に含まれます。これはa_n(x)がa(x)に収束していることを意味します。

（証明のイメージは、a(x)のどんなに小さい近傍を取って来ても、多項式の列{a_n(x)}の n が十分大きい項は全てその近傍に含まれている、という感じです）

このような位相を、R=ℂのときは h進位相（h-adic topology）といいます。（h なのは変数をプランク定数を意識した h にしてℂ[[h]]と表すことが多いからで、ℂ上の代数に h を加えて非可換な代数を構成する「量子化」と関係があります。）

またR=ℤ/pℤのときは、R[[x]]はp進整数環とみなせます。