今回は「多項式 (polynomials)」の計算について取り上げます。NumPy には多項式を取り扱う便利な関数が用意されているので、私たちがプログラムを作る必要はありませんが、Python (NumPy) のお勉強ということで、四則演算を行う簡単なプログラムを作ってみましょう。
n を非負の整数とすると、一変数の多項式は次のように表すことができます。
cnxn + cn-1xn-1 + ... + c2x2 + c1x + c0
一変数の多項式は係数 cn, ..., c0 によって決定できるので、係数を格納した配列で多項式を表すことができます。今回は NumPy に合わせて、[cn, ..., c0] のように係数を格納することにします。すると、多項式の足し算 (引き算) は、対応する係数の足し算 (引き算) になります。掛け算と割り算は筆算と同じ方法でプログラムすることにします。
足し算と引き算は簡単です。プログラムは次のようになります。
リスト : 多項式の足し算と引き算
# 正規化
def normalize(xs):
xs = np.trim_zeros(xs, 'f')
if not len(xs): xs = np.array([0.])
return xs
# 足し算
def polyadd(xs, ys):
xn = len(xs)
yn = len(ys)
if xn >= yn:
zs = xs.copy()
zs[-yn:] += ys
else:
zs = ys.copy()
zs[-xn:] += xs
return normalize(zs)
# 引き算
def polysub(xs, ys):
return polyadd(xs, -ys)
関数名は NumPy に合わせて polyadd と polysub としました。引数 xs, ys の桁数 (長さ) を変数 xn, yn にセットします。桁数が多い引数を copy() して変数 zs にセットし、そこに短いほうの引数を加算します。このとき、係数の位置を合わせることに注意してください。もし、計算結果が 0 になったら np.array([0.]) を返します。この処理を関数 normalize() で行っています。
関数 trim_zeros() は配列の両端から連続した 0 を取り除きます。第 2 引数に 'f' を指定すると左端から、'b' を指定すると右端から削除します。trim_zeros() はビューを返すことに注意してください。空の配列になったら np.array([0.]) を返します。引き算を行う polysub(xs, ys) は polyadd(xs, -ys) を呼び出すだけです。
次は掛け算を行う関数 polymul() を作ります。名前は NumPy に合わせました。
リスト : 掛け算
def polymul(xs, ys):
xn = len(xs)
yn = len(ys)
zn = xn + yn - 1
zs = np.full(zn, 0.)
for i in range(yn):
zs[i:i+xn] += ys[i] * xs
return normalize(zs)
引数 xs, ys の桁数を xn, yn とすると、xs, ys の掛け算は xn - 1 + yn - 1 + 1 = xn + yn - 1 = zn 桁になります。初期値 0 で大きさ zn の配列を full() で生成して変数 zs にセットします。あとは順番に ys の係数を xs に掛け算して、その値を zs に足し算していくだけです。
次は割り算を行う関数 polydiv() を作ります。名前は NumPy に合わせました。
リスト : 割り算
def polydiv(xs, ys):
xn = len(xs)
yn = len(ys)
zs = xs.copy()
qs = []
for _ in range(xn - yn + 1):
temp = zs[0] / ys[0]
zs[:yn] -= temp * ys
qs.append(temp)
zs = zs[1:]
if qs == []: qs = [0.]
return np.array(qs), normalize(zs)
xs を zs にコピーし、qs に空リストをセットします。qs が商で、zs が余りになります。for ループの中で、zs[0] / ys[0] の値を temp にセットして、zs から temp * ys を引き算します。これで zs[0] の係数を 0 にすることができます。あとは temp を qs に追加して、zs の桁を一つ削除します。for ループのあと、qs が空リストならば、qs に np.array([0.]) をセットします。最後に qs と zs を返します。
変数 x に値を代入して多項式の値を求める場合、x の累乗をそのまま計算するよりも、多項式を次のように変形したほうが効率的です。
c4x4 +c3x3 + c2x2 + c1x + c0 = (((c4x + c3)x + c2)x + c1)x + c0
これを「ホーナー法 (Horner's method)」といいます。ホーナー法は functools.reduce() を使うと簡単です。
リスト : 多項式の値を求める (ホーナー法)
def polyval(xs, n):
return functools.reduce(lambda x, y: x * n + y, xs)
ラムダ式の x に累積値、y に xs の値が渡されます。最初、xs[0] が x に、xs[1] が y に渡されるので、ホーナー法として正常に動作します。
それでは実際に試してみましょう。次に示す簡単なテストを行ってみました。
リスト : 簡単なテスト
a = np.array([1.,2.,3.,4.])
b = np.array([5.,6.,7.,8.,9])
for x, y in [(a, a), (a, b), (b, a), (b, b)]:
print(x, "+", y, "=", polyadd(x, y))
print(x, "+", y, "=", np.polyadd(x, y))
print(x, "-", y, "=", polysub(x, y))
print(x, "-", y, "=", np.polysub(x, y))
print(x, "*", y, "=", polymul(x, y))
print(x, "*", y, "=", np.polymul(x, y))
p, q = polydiv(x, y)
print(x, "/", y, "=", p, q)
p, q = np.polydiv(x, y)
print(x, "/", y, "=", p, q)
print("polyval")
print(polyval(a, 3))
print(np.polyval(a, 3))
print(polyval(b, 5))
print(np.polyval(b, 5))
[ 1. 2. 3. 4.] + [ 1. 2. 3. 4.] = [ 2. 4. 6. 8.] [ 1. 2. 3. 4.] + [ 1. 2. 3. 4.] = [ 2. 4. 6. 8.] [ 1. 2. 3. 4.] - [ 1. 2. 3. 4.] = [ 0.] [ 1. 2. 3. 4.] - [ 1. 2. 3. 4.] = [ 0. 0. 0. 0.] [ 1. 2. 3. 4.] * [ 1. 2. 3. 4.] = [ 1. 4. 10. 20. 25. 24. 16.] [ 1. 2. 3. 4.] * [ 1. 2. 3. 4.] = [ 1. 4. 10. 20. 25. 24. 16.] [ 1. 2. 3. 4.] / [ 1. 2. 3. 4.] = [ 1.] [ 0.] [ 1. 2. 3. 4.] / [ 1. 2. 3. 4.] = [ 1.] [ 0.] [ 1. 2. 3. 4.] + [ 5. 6. 7. 8. 9.] = [ 5. 7. 9. 11. 13.] [ 1. 2. 3. 4.] + [ 5. 6. 7. 8. 9.] = [ 5. 7. 9. 11. 13.] [ 1. 2. 3. 4.] - [ 5. 6. 7. 8. 9.] = [-5. -5. -5. -5. -5.] [ 1. 2. 3. 4.] - [ 5. 6. 7. 8. 9.] = [-5. -5. -5. -5. -5.] [ 1. 2. 3. 4.] * [ 5. 6. 7. 8. 9.] = [ 5. 16. 34. 60. 70. 70. 59. 36.] [ 1. 2. 3. 4.] * [ 5. 6. 7. 8. 9.] = [ 5. 16. 34. 60. 70. 70. 59. 36.] [ 1. 2. 3. 4.] / [ 5. 6. 7. 8. 9.] = [ 0.] [ 1. 2. 3. 4.] [ 1. 2. 3. 4.] / [ 5. 6. 7. 8. 9.] = [ 0.] [ 1. 2. 3. 4.] [ 5. 6. 7. 8. 9.] + [ 1. 2. 3. 4.] = [ 5. 7. 9. 11. 13.] [ 5. 6. 7. 8. 9.] + [ 1. 2. 3. 4.] = [ 5. 7. 9. 11. 13.] [ 5. 6. 7. 8. 9.] - [ 1. 2. 3. 4.] = [ 5. 5. 5. 5. 5.] [ 5. 6. 7. 8. 9.] - [ 1. 2. 3. 4.] = [ 5. 5. 5. 5. 5.] [ 5. 6. 7. 8. 9.] * [ 1. 2. 3. 4.] = [ 5. 16. 34. 60. 70. 70. 59. 36.] [ 5. 6. 7. 8. 9.] * [ 1. 2. 3. 4.] = [ 5. 16. 34. 60. 70. 70. 59. 36.] [ 5. 6. 7. 8. 9.] / [ 1. 2. 3. 4.] = [ 5. -4.] [ 25.] [ 5. 6. 7. 8. 9.] / [ 1. 2. 3. 4.] = [ 5. -4.] [ 25.] [ 5. 6. 7. 8. 9.] + [ 5. 6. 7. 8. 9.] = [ 10. 12. 14. 16. 18.] [ 5. 6. 7. 8. 9.] + [ 5. 6. 7. 8. 9.] = [ 10. 12. 14. 16. 18.] [ 5. 6. 7. 8. 9.] - [ 5. 6. 7. 8. 9.] = [ 0.] [ 5. 6. 7. 8. 9.] - [ 5. 6. 7. 8. 9.] = [ 0. 0. 0. 0. 0.] [ 5. 6. 7. 8. 9.] * [ 5. 6. 7. 8. 9.] = [ 25. 60. 106. 164. 235. 220. 190. 144. 81.] [ 5. 6. 7. 8. 9.] * [ 5. 6. 7. 8. 9.] = [ 25. 60. 106. 164. 235. 220. 190. 144. 81.] [ 5. 6. 7. 8. 9.] / [ 5. 6. 7. 8. 9.] = [ 1.] [ 0.] [ 5. 6. 7. 8. 9.] / [ 5. 6. 7. 8. 9.] = [ 1.] [ 0.] polyval 58.0 58.0 4099.0 4099.0
正常に動作していますね。
次は有限体 GF(p) で多項式の四則演算を定義してみましょう。ここでは p を素数に限定します。プログラムは次のようになります。
リスト : GF(p) で多項式を計算する (p は素数)
class GF:
def __init__(self, n):
self.num = n
self.inv = make_inv(n)
# 足し算: polyadd(self, xs, ys)
# 引き算: polysub(self, xs, ys)
# 掛け算: polymul(self, xs, ys)
# 割り算: polydiv(self, xs, ys)
クラス GF は有限体 (ガロア体) を表します。__init__() の引数 n は GF(p) の p を表します。n を self.num にセットし、make_inv() で逆数表を生成して self.inv にセットします。make_inv() は N 色ライツアウトの解法 で作成した関数です。足し算と掛け算は、係数を計算したあと剰余を求めるように修正します。引き算と割り算は、負数と逆数を使って、足し算と掛け算に変換してから計算します。詳細は プログラムリスト２ をお読みくださいませ。
それでは実際に試してみましょう。次に示す簡単なテストを行ってみました。
リスト : 簡単なテスト
def test(p, xs, ys):
xs = np.array(xs)
ys = np.array(ys)
gf = GF(p)
for x, y in [(xs, xs), (xs, ys), (ys, xs), (ys, ys)]:
print(x, "+", y, "=", gf.polyadd(x, y))
print(x, "-", y, "=", gf.polysub(x, y))
print(x, "*", y, "=", gf.polymul(x, y))
p, q = gf.polydiv(x, y)
print(x, "/", y, "=", p, q)
test(2, [1,0,1,1], [1,0,1])
test(3, [2,1,0,2], [2,0,1])
test(5, [4,3,2,1], [3,2,1])
[1 0 1 1] + [1 0 1 1] = [0] [1 0 1 1] - [1 0 1 1] = [0] [1 0 1 1] * [1 0 1 1] = [1 0 0 0 1 0 1] [1 0 1 1] / [1 0 1 1] = [1] [0] [1 0 1 1] + [1 0 1] = [1 1 1 0] [1 0 1 1] - [1 0 1] = [1 1 1 0] [1 0 1 1] * [1 0 1] = [1 0 0 1 1 1] [1 0 1 1] / [1 0 1] = [1 0] [1] [1 0 1] + [1 0 1 1] = [1 1 1 0] [1 0 1] - [1 0 1 1] = [1 1 1 0] [1 0 1] * [1 0 1 1] = [1 0 0 1 1 1] [1 0 1] / [1 0 1 1] = [0] [1 0 1] [1 0 1] + [1 0 1] = [0] [1 0 1] - [1 0 1] = [0] [1 0 1] * [1 0 1] = [1 0 0 0 1] [1 0 1] / [1 0 1] = [1] [0] [2 1 0 2] + [2 1 0 2] = [1 2 0 1] [2 1 0 2] - [2 1 0 2] = [0] [2 1 0 2] * [2 1 0 2] = [1 1 1 2 1 0 1] [2 1 0 2] / [2 1 0 2] = [1] [0] [2 1 0 2] + [2 0 1] = [2 0 0 0] [2 1 0 2] - [2 0 1] = [2 2 0 1] [2 1 0 2] * [2 0 1] = [1 2 2 2 0 2] [2 1 0 2] / [2 0 1] = [1 2] [2 0] [2 0 1] + [2 1 0 2] = [2 0 0 0] [2 0 1] - [2 1 0 2] = [1 1 0 2] [2 0 1] * [2 1 0 2] = [1 2 2 2 0 2] [2 0 1] / [2 1 0 2] = [0] [2 0 1] [2 0 1] + [2 0 1] = [1 0 2] [2 0 1] - [2 0 1] = [0] [2 0 1] * [2 0 1] = [1 0 1 0 1] [2 0 1] / [2 0 1] = [1] [0] [4 3 2 1] + [4 3 2 1] = [3 1 4 2] [4 3 2 1] - [4 3 2 1] = [0] [4 3 2 1] * [4 3 2 1] = [1 4 0 0 0 4 1] [4 3 2 1] / [4 3 2 1] = [1] [0] [4 3 2 1] + [3 2 1] = [4 1 4 2] [4 3 2 1] - [3 2 1] = [4 0 0 0] [4 3 2 1] * [3 2 1] = [2 2 1 0 4 1] [4 3 2 1] / [3 2 1] = [3 4] [1 2] [3 2 1] + [4 3 2 1] = [4 1 4 2] [3 2 1] - [4 3 2 1] = [1 0 0 0] [3 2 1] * [4 3 2 1] = [2 2 1 0 4 1] [3 2 1] / [4 3 2 1] = [0] [3 2 1] [3 2 1] + [3 2 1] = [1 4 2] [3 2 1] - [3 2 1] = [0] [3 2 1] * [3 2 1] = [4 2 0 4 1] [3 2 1] / [3 2 1] = [1] [0]
正常に動作しているようです。興味のある方はいろいろ試してみてください。
次は GF(p) 上で多項式を因数分解するメソッド factorization() を作りましょう。n 次の多項式において、n-1 次以下の多項式で割り切れないものを「既約多項式」といいます。既約多項式は因数分解することができません。これを「既約」といい、因数分解できることを「可約」といいます。たとえば、GF(2) 上の 2 次多項式を因数分解すると次のようになります
式 1, 2, 3 は 1 次多項式 x または x + 1 で割り切れるので可約ですが、式 4 は x と x + 1 で割り切れないので既約です。既約多項式は多項式の素数みたいなもので、これを使って素因数分解のように「試し割り」で多項式を因数分解することができます。なお、今回作成するプログラムは素朴な「試し割り」なので、効率は良くありません。ご注意くださいませ。
具体的にいうと、n 次の多項式 xs を因数分解する場合、1 次の多項式 ys から順番に割り算 (xs / ys) していきます。n/2 次の多項式まで試しても割り切れない場合、その式は因数分解することができません。xs / ys の余りが 0 ならば、ys は xs の因子の一つです。可能な限り ys で除算して、余りが 0 でなくなったならば次の多項式で試し割りします。xs の次数は低くなり、ys の次数は高くなるので、そのうちに xs の次数は ys 以下になります。ここで因数分解を終了します。
プログラムは次のようになります。
リスト : 多項式の因数分解
def factorization(self, xs):
ps = []
g = functools.reduce(lambda x, y: math.gcd(x, y), xs)
if g > 1:
xs = xs // g
ps.append((np.array([g]), 1))
m = (len(xs) - 1) // 2 + 1
zs = itertools.product(range(self.num), repeat=m)
for _ in range(self.num): next(zs)
for z in zs:
ys = np.trim_zeros(np.array(z), 'f')
if len(xs) <= len(ys): break
c = 0
while len(xs) > len(ys):
p, q = self.polydiv(xs, ys)
if q[0] == 0:
xs = p
c += 1
else:
break
if c > 0:
if np.all(xs == ys):
c += 1
xs = np.array([1])
ps.append((ys, c))
if len(xs) > 1: ps.append((xs, 1))
return ps
変数 ps のリストに因子を格納します。最初に係数の最大公約数を求め、1 よりも大きければ、その値で xs を割り算します。そして、ps に (np.array([g]), 1) をセットします。次に、除数の最大桁数を求めて変数 m にセットします。たとえば、xs の次数が 5 の場合、2 次多項式まで試し割りすればいいので、最大桁数は (len(xs) - 1) // 2 + 1 = 2 + 1 = 3 になります。xs の次数が 6 の場合は 3 + 1 = 4 になります。
除数の生成は itertools.product() を使うと簡単です。product() はリストとリストの直積を求める関数ですが、repeat を指定すると重複順列を生成することができます。簡単な例を示しましょう。
>>> for xs in itertools.product([1,2],[3,4]): print(xs) ... (1, 3) (1, 4) (2, 3) (2, 4) >>> for z in itertools.product([0,1],repeat=3): print(z) ... (0, 0, 0) (0, 0, 1) (0, 1, 0) (0, 1, 1) (1, 0, 0) (1, 0, 1) (1, 1, 0) (1, 1, 1)
repeat に m を指定すれば、m-1 次の多項式を生成することができます。定数だけの多項式は不要なので、next() を使って先頭から self.num 個分取り除いています。生成した除数を for ループで順番に取り出し、それを NumPy の配列に変換し、trim_zeros() で先頭の 0 を取り除いてから変数 ys にセットします。あとは、素因数分解のように xs を ys で除算していくだけです。
それでは実際に試してみましょう。
>>> gf2 = GF(2) >>> gf2.factorization(np.array([1,1])) [(array([1, 1]), 1)] >>> gf2.factorization(np.array([1,1,1])) [(array([1, 1, 1]), 1)] >>> gf2.factorization(np.array([1,1,1,1])) [(array([1, 1]), 3)] >>> gf2.factorization(np.array([1,1,1,1,1])) [(array([1, 1, 1, 1, 1]), 1)] >>> gf2.factorization(np.array([1,1,1,1,1,1])) [(array([1, 1]), 1), (array([1, 1, 1]), 2)] >>> gf2.factorization(np.array([1,1,1,1,1,1,1])) [(array([1, 0, 1, 1]), 1), (array([1, 1, 0, 1]), 1)] >>> gf2.factorization(np.array([1,1,1,1,1,1,1,1])) [(array([1, 1]), 7)] >>> gf3 = GF(3) >>> gf3.factorization(np.array([2,2,2,2,2,2,2,2])) [(array([2]), 1), (array([1, 1]), 1), (array([1, 0, 1]), 1), (array([1, 1, 2]), 1), (array([1, 2, 2]), 1)] >>> gf3.factorization(np.array([1,1,1,1,1,1,1,1])) [(array([1, 1]), 1), (array([1, 0, 1]), 1), (array([1, 1, 2]), 1), (array([1, 2, 2]), 1)] >>> gf3.factorization(np.array([2,1,2,1,2,1,2])) [(array([2, 1, 2, 1, 2, 1, 2]), 1)] >>> gf3.factorization(np.array([2,1,2,0,2,1,2])) [(array([1, 1]), 2), (array([1, 1, 2]), 1), (array([2, 1, 1]), 1)] >>> gf3.factorization(np.array([1,1,1,1,1,1])) [(array([1, 1]), 3), (array([1, 2]), 2)] >>> gf5 = GF(5) >>> gf5.factorization(np.array([4,4,4,4,4,4,4,4,4])) [(array([4]), 1), (array([1, 1, 1]), 1), (array([1, 0, 0, 1, 0, 0, 1]), 1)] >>> gf5.factorization(np.array([4,2,4,2,4,2,4,2,4])) [(array([2]), 1), (array([1, 0, 4, 1, 2]), 1), (array([2, 1, 4, 0, 1]), 1)] >>> gf5.factorization(np.array([4,3,2,1,4,3,2,1])) [(array([1, 4]), 3), (array([1, 0, 2]), 1), (array([4, 0, 2]), 1)] >>> gf5.factorization(np.array([1,2,3,4,3,2,1])) [(array([1, 1]), 2), (array([1, 2]), 2), (array([1, 3]), 2)]
ところで、参考 URL 2 によると、GF(p) 上で多項式 xpn - x を因数分解すると、因子に n 次の既約多項式 (xn の係数は 1) がすべて現れるそうです。もちろん、xm - 1, (m = pn - 1) を因数分解してもかまいません。以下にプログラムと実行結果を示します。
リスト : xpn - x を因数分解する
def gf_fact1(p, n):
gf = GF(p)
m = p ** n + 1
xs = np.full(m, 0, dtype = np.int_)
xs[0] = 1
xs[-2] = p - 1
return gf.factorization(xs)
for xs in gf_fact1(2, 3):
if len(xs[0]) == 4: print(xs[0])
for xs in gf_fact1(2, 4):
if len(xs[0]) == 5: print(xs[0])
for xs in gf_fact1(3, 2):
if len(xs[0]) == 3: print(xs[0])
for xs in gf_fact1(3, 3):
if len(xs[0]) == 4: print(xs[0])
for xs in gf_fact1(5, 2):
if len(xs[0]) == 3: print(xs[0])
for xs in gf_fact1(5, 3):
if len(xs[0]) == 4: print(xs[0])
[1 0 1 1] [1 1 0 1] [1 0 0 1 1] [1 1 0 0 1] [1 1 1 1 1] [1 0 1] [1 1 2] [1 2 2] [1 0 2 1] [1 0 2 2] [1 1 0 2] [1 1 1 2] [1 1 2 1] [1 2 0 1] [1 2 1 1] [1 2 2 2] [1 0 2] [1 0 3] [1 1 1] [1 1 2] [1 2 3] [1 2 4] [1 3 3] [1 3 4] [1 4 1] [1 4 2] [1 0 1 1] [1 0 1 4] [1 0 2 1] [1 0 2 4] [1 0 3 2] [1 0 3 3] [1 0 4 2] [1 0 4 3] [1 1 0 1] [1 1 0 2] [1 1 1 3] [1 1 1 4] [1 1 3 1] [1 1 3 4] [1 1 4 1] [1 1 4 3] [1 2 0 1] [1 2 0 3] [1 2 1 3] [1 2 1 4] [1 2 2 2] [1 2 2 3] [1 2 4 2] [1 2 4 4] [1 3 0 2] [1 3 0 4] [1 3 1 1] [1 3 1 2] [1 3 2 2] [1 3 2 3] [1 3 4 1] [1 3 4 3] [1 4 0 3] [1 4 0 4] [1 4 1 1] [1 4 1 2] [1 4 3 1] [1 4 3 4] [1 4 4 2] [1 4 4 4]
GF(p) 上で n 次の既約多項式 xs (xn の係数は 1) を考えます。xs が多項式 xm - 1 を割り切るとき、最小の m を「周期」といいます。そして、m が pn - 1 と等しい場合、その既約多項式 xs を「原始多項式 (primitive polynomial)」といいます。n 次の既約多項式であれば、m = pn - 1 の多項式を割り切りますが、m がそれよりも小さい値の多項式を割り切る既約多項式も存在します。
簡単な例を示しましょう。
GF(3) の既約多項式 f(x) = x2 + 1 (x1 - 1) / f(x) = 0, x + 2 (x2 - 1) / f(x) = 1, 1 (x3 - 1) / f(x) = x, 2x + 2 (x4 - 1) / f(x) = x2 + 2, 0 <== 周期は 4 GF(3) の既約多項式 f(x) = x2 + x + 2 (x1 - 1) / f(x) = 0, x + 2 (x2 - 1) / f(x) = 1, 2x (x3 - 1) / f(x) = x + 2, 2x + 1 (x4 - 1) / f(x) = x2 + 2x + 2, 1 (x5 - 1) / f(x) = x3 + 2x2 + 2x, 2x + 2 (x6 - 1) / f(x) = x4 + 2x3 + 2x2 + 2, x + 1 (x7 - 1) / f(x) = x5 + 2x4 + 2x3 + 2x + 1, x (x8 - 1) / f(x) = x6 + 2x5 + 2x4 + 2x2 + x + 1, 0 <== 周期は 8 (原始多項式)
効率はとても悪いですが、これをそのままプログラムするだけで n 次の原始多項式を求めることができます。
リスト : n 次の原始多項式を求める
def primitive_polynomial(p, n):
gf = GF(p)
xs = np.full(p ** n, 0, dtype = np.int_)
xs[0] = 1
xs[-1] = p - 1
for y in gf.factorization(xs):
ys = y[0]
if len(ys) < n + 1: continue
for m in range(n + 1, p ** n - 1):
zs = np.full(m, 0, dtype = np.int_)
zs[0] = 1
zs[-1] = p - 1
_, q = gf.polydiv(zs, ys)
if q[0] == 0: break
else:
yield ys
xk - 1 (k = pn - 1) を因数分解して n 次の既約多項式 ys を求めます。2 番目の for ループで、xm - 1 を ys で割り算します。m が p ** n - 2 まで割り算して、余りが 0 にならなければ原始多項式です。途中で余りが 0 になったら原始多項式ではありません。break でループを脱出します。
それでは実際に試してみましょう。
>>> for xs in primitive_polynomial(2, 2): print(xs) ... [1 1 1] >>> for xs in primitive_polynomial(2, 3): print(xs) ... [1 0 1 1] [1 1 0 1] >>> for xs in primitive_polynomial(2, 4): print(xs) ... [1 0 0 1 1] [1 1 0 0 1] >>> for xs in primitive_polynomial(2, 5): print(xs) ... [1 0 0 1 0 1] [1 0 1 0 0 1] [1 0 1 1 1 1] [1 1 0 1 1 1] [1 1 1 0 1 1] [1 1 1 1 0 1] >>> for xs in primitive_polynomial(2, 6): print(xs) ... [1 0 0 0 0 1 1] [1 0 1 1 0 1 1] [1 1 0 0 0 0 1] [1 1 0 0 1 1 1] [1 1 0 1 1 0 1] [1 1 1 0 0 1 1] >>> for xs in primitive_polynomial(3, 2): print(xs) ... [1 1 2] [1 2 2] >>> for xs in primitive_polynomial(3, 3): print(xs) ... [1 0 2 1] [1 1 2 1] [1 2 0 1] [1 2 1 1] >>> for xs in primitive_polynomial(3, 4): print(xs) ... [1 0 0 1 2] [1 0 0 2 2] [1 1 0 0 2] [1 1 1 2 2] [1 1 2 2 2] [1 2 0 0 2] [1 2 1 1 2] [1 2 2 1 2] >>> for xs in primitive_polynomial(5, 2): print(xs) ... [1 1 2] [1 2 3] [1 3 3] [1 4 2] >>> for xs in primitive_polynomial(5, 3): print(xs) ... [1 0 3 2] [1 0 3 3] [1 0 4 2] [1 0 4 3] [1 1 0 2] [1 1 1 3] [1 1 4 3] [1 2 0 3] [1 2 1 3] [1 2 2 2] [1 2 2 3] [1 2 4 2] [1 3 0 2] [1 3 1 2] [1 3 2 2] [1 3 2 3] [1 3 4 3] [1 4 0 3] [1 4 1 2] [1 4 4 2]
4 の剰余 (0, 1, 2, 3) は GF(4) になりませんが、GF(2) に新しい数を付け加えることで、数が 4 つある有限体を作ることができます。新しい数を付加して体を拡張することを「体の拡大」といい、拡大された体を「拡大体」といいます。たとえば、実数の世界では代数方程式 x2 + 1 = 0 に解はありませんが、虚数 i を導入した複素数の世界では解くことができますね。複素数の世界は実数の拡大体になるわけです。
GF(2) を拡大する場合、GF(2) 上の原始多項式 f(x) を考えます。2 次の原始多項式には x2 + x + 1 がありますが、GF(2) 上で f(x) = 0 となる解はありません。この式の解の一つを a としましょう。a2 + a + 1 = 0 が成り立つので、a2 = a + 1 であることがわかります。次に、a の累乗を計算します。
a1 = a a2 = a + 1 a3 = a * a2 = a * (a + 1) = a + 1 + a = 1 a4 = a1+3 = a a5 = a2+3 = a + 1 a6 = a3+3 = 1
a, a + 1, 1 が順番に現れています。これを多項式で表すと 1, x, x + 1 になり、1 次多項式がすべて現れていることがわかります。このように、原始多項式の解の累乗を計算すると、n 次未満の多項式がすべて現れます。これに 0 を加えた 4 つの要素 {0, 1, x, x + 1} を数として加法と乗法を考えてみましょう。基本的には式を計算して、それを既約多項式 x2 + x + 1 で割った剰余が値になります。下図を見てください。
+ | 0 1 x x+1 * | 0 1 x x+1 ----+----------------- ----+------------------ 0 | 0 1 x x+1 0 | 0 0 0 0 1 | 1 0 x+1 x 1 | 0 1 x x+1 x | x x+1 0 1 x | 0 x x+1 1 x+1 | x+1 x 1 0 x+1 | 0 x+1 1 x
加法の場合は特に簡単で、係数ごとに加算して mod 2 をとる、つまり係数ごとの排他的論理和になります。乗算の場合は累乗の形式に変換して計算すると簡単です。つまり、am * an = an+m = a(n+m) mod 3 になります。a2 = a + 1 の関係を使って式を簡約してもかまいません。
減法の場合、加法表より負数が自分自身であることがわかるので、加法と同様に係数ごとの排他的論理和になります。徐算の場合、乗法表より 1, x, x + 1 の逆数がすべて存在するので、徐算が成り立つことがわかります。つまり、{0, 1, x, x + 1} は有限体になるわけです。多項式の係数をビットで表すと、加法と乗法は次のようになります。
+ | 00 01 10 11 * | 00 01 10 11 ----+----------------- ----+------------------ 00 | 00 01 10 11 00 | 00 00 00 00 01 | 01 00 11 10 01 | 00 01 10 11 10 | 10 11 00 01 10 | 00 10 11 01 11 | 11 10 01 00 11 | 00 11 01 10
p を素数とし、n を自然数とすると、要素数が pn の有限体は必ず存在することが知られています。これを GF(pn) と表記します。GF(2n) はコンピュータとの相性が良いので、エラー検出、乱数の生成、暗号など多くの分野で応用されています。
ご参考までに、加法表と乗法表を作成するプログラムとその実行例を示します。
リスト : 加法表と乗法表の作成
def add_table(p, zs):
gf = GF(p)
yss = [np.array(y) for y in itertools.product(range(p), repeat=len(zs) - 1)]
tbl = []
for xs in yss:
ls = []
for ys in yss:
_, q = gf.polydiv(gf.polyadd(xs, ys), zs)
ls.append(functools.reduce(lambda x, y: x * p + y, q))
tbl.append(ls)
return np.array(tbl)
def mul_table(p, zs):
gf = GF(p)
yss = [np.array(y) for y in itertools.product(range(p), repeat=len(zs) - 1)]
tbl = []
for xs in yss:
ls = []
for ys in yss:
_, q = gf.polydiv(gf.polymul(xs, ys), zs)
ls.append(functools.reduce(lambda x, y: x * p + y, q))
tbl.append(ls)
return np.array(tbl)
>>> add_table(2, np.array([1,1,1]))
array([[0, 1, 2, 3],
[1, 0, 3, 2],
[2, 3, 0, 1],
[3, 2, 1, 0]])
>>> mul_table(2, np.array([1,1,1]))
array([[0, 0, 0, 0],
[0, 1, 2, 3],
[0, 2, 3, 1],
[0, 3, 1, 2]])
>>> add_table(2, np.array([1,0,1,1]))
array([[0, 1, 2, 3, 4, 5, 6, 7],
[1, 0, 3, 2, 5, 4, 7, 6],
[2, 3, 0, 1, 6, 7, 4, 5],
[3, 2, 1, 0, 7, 6, 5, 4],
[4, 5, 6, 7, 0, 1, 2, 3],
[5, 4, 7, 6, 1, 0, 3, 2],
[6, 7, 4, 5, 2, 3, 0, 1],
[7, 6, 5, 4, 3, 2, 1, 0]])
>>> mul_table(2, np.array([1,0,1,1]))
array([[0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 2, 3, 4, 5, 6, 7],
[0, 2, 4, 6, 3, 1, 7, 5],
[0, 3, 6, 5, 7, 4, 1, 2],
[0, 4, 3, 7, 6, 2, 5, 1],
[0, 5, 1, 4, 2, 7, 3, 6],
[0, 6, 7, 1, 5, 3, 2, 4],
[0, 7, 5, 2, 1, 6, 4, 3]])
>>> add_table(2, np.array([1,0,0,1,1]))
array([[ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15],
[ 1, 0, 3, 2, 5, 4, 7, 6, 9, 8, 11, 10, 13, 12, 15, 14],
[ 2, 3, 0, 1, 6, 7, 4, 5, 10, 11, 8, 9, 14, 15, 12, 13],
[ 3, 2, 1, 0, 7, 6, 5, 4, 11, 10, 9, 8, 15, 14, 13, 12],
[ 4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11],
[ 5, 4, 7, 6, 1, 0, 3, 2, 13, 12, 15, 14, 9, 8, 11, 10],
[ 6, 7, 4, 5, 2, 3, 0, 1, 14, 15, 12, 13, 10, 11, 8, 9],
[ 7, 6, 5, 4, 3, 2, 1, 0, 15, 14, 13, 12, 11, 10, 9, 8],
[ 8, 9, 10, 11, 12, 13, 14, 15, 0, 1, 2, 3, 4, 5, 6, 7],
[ 9, 8, 11, 10, 13, 12, 15, 14, 1, 0, 3, 2, 5, 4, 7, 6],
[10, 11, 8, 9, 14, 15, 12, 13, 2, 3, 0, 1, 6, 7, 4, 5],
[11, 10, 9, 8, 15, 14, 13, 12, 3, 2, 1, 0, 7, 6, 5, 4],
[12, 13, 14, 15, 8, 9, 10, 11, 4, 5, 6, 7, 0, 1, 2, 3],
[13, 12, 15, 14, 9, 8, 11, 10, 5, 4, 7, 6, 1, 0, 3, 2],
[14, 15, 12, 13, 10, 11, 8, 9, 6, 7, 4, 5, 2, 3, 0, 1],
[15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]])
>>> mul_table(2, np.array([1,0,0,1,1]))
array([[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15],
[ 0, 2, 4, 6, 8, 10, 12, 14, 3, 1, 7, 5, 11, 9, 15, 13],
[ 0, 3, 6, 5, 12, 15, 10, 9, 11, 8, 13, 14, 7, 4, 1, 2],
[ 0, 4, 8, 12, 3, 7, 11, 15, 6, 2, 14, 10, 5, 1, 13, 9],
[ 0, 5, 10, 15, 7, 2, 13, 8, 14, 11, 4, 1, 9, 12, 3, 6],
[ 0, 6, 12, 10, 11, 13, 7, 1, 5, 3, 9, 15, 14, 8, 2, 4],
[ 0, 7, 14, 9, 15, 8, 1, 6, 13, 10, 3, 4, 2, 5, 12, 11],
[ 0, 8, 3, 11, 6, 14, 5, 13, 12, 4, 15, 7, 10, 2, 9, 1],
[ 0, 9, 1, 8, 2, 11, 3, 10, 4, 13, 5, 12, 6, 15, 7, 14],
[ 0, 10, 7, 13, 14, 4, 9, 3, 15, 5, 8, 2, 1, 11, 6, 12],
[ 0, 11, 5, 14, 10, 1, 15, 4, 7, 12, 2, 9, 13, 6, 8, 3],
[ 0, 12, 11, 7, 5, 9, 14, 2, 10, 6, 1, 13, 15, 3, 4, 8],
[ 0, 13, 9, 4, 1, 12, 8, 5, 2, 15, 11, 6, 3, 14, 10, 7],
[ 0, 14, 15, 1, 13, 3, 2, 12, 9, 7, 6, 8, 4, 10, 11, 5],
[ 0, 15, 13, 2, 9, 6, 4, 11, 1, 14, 12, 3, 8, 7, 5, 10]])
>>> add_table(3, np.array([1,1,2]))
array([[0, 1, 2, 3, 4, 5, 6, 7, 8],
[1, 2, 0, 4, 5, 3, 7, 8, 6],
[2, 0, 1, 5, 3, 4, 8, 6, 7],
[3, 4, 5, 6, 7, 8, 0, 1, 2],
[4, 5, 3, 7, 8, 6, 1, 2, 0],
[5, 3, 4, 8, 6, 7, 2, 0, 1],
[6, 7, 8, 0, 1, 2, 3, 4, 5],
[7, 8, 6, 1, 2, 0, 4, 5, 3],
[8, 6, 7, 2, 0, 1, 5, 3, 4]])
>>> mul_table(3, np.array([1,1,2]))
array([[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 2, 3, 4, 5, 6, 7, 8],
[0, 2, 1, 6, 8, 7, 3, 5, 4],
[0, 3, 6, 7, 1, 4, 5, 8, 2],
[0, 4, 8, 1, 5, 6, 2, 3, 7],
[0, 5, 7, 4, 6, 2, 8, 1, 3],
[0, 6, 3, 5, 2, 8, 7, 4, 1],
[0, 7, 5, 8, 3, 1, 4, 2, 6],
[0, 8, 4, 2, 7, 3, 1, 6, 5]])
#
# poly.py : 多項式の計算
#
# Copyright (C) 2018-2023 Makoto Hiroi
#
import numpy as np
import functools
# 係数
# [cn, cn-1, ..., c1, c]
# 正規化
def normalize(xs):
xs = np.trim_zeros(xs, 'f')
if not len(xs): xs = np.array([0.])
return xs
# 足し算
def polyadd(xs, ys):
xn = len(xs)
yn = len(ys)
if xn >= yn:
zs = xs.copy()
zs[-yn:] += ys
else:
zs = ys.copy()
zs[-xn:] += xs
return normalize(zs)
# 引き算
def polysub(xs, ys):
return polyadd(xs, -ys)
# 掛け算
def polymul(xs, ys):
xn = len(xs)
yn = len(ys)
zn = xn + yn - 1
zs = np.full(zn, 0.)
for i in range(yn):
zs[i:i+xn] += ys[i] * xs
return normalize(zs)
# 割り算
def polydiv(xs, ys):
xn = len(xs)
yn = len(ys)
zs = xs.copy()
qs = []
for _ in range(xn - yn + 1):
temp = zs[0] / ys[0]
zs[:yn] -= temp * ys
qs.append(temp)
zs = zs[1:]
if qs == []: qs = [0.]
return np.array(qs), normalize(zs)
# 多項式の値を求める (ホーナー法)
def polyval(xs, n):
return functools.reduce(lambda x, y: x * n + y, xs)
#
# polygf.py : GF(p) で多項式を計算する (p は素数)
#
# Copyright (C) 2018-2023 Makoto Hiroi
#
import numpy as np
import functools
import itertools
import math
# 拡張ユークリッドの互除法
def ext_euclid(a, b):
xs = a, 1, 0
ys = b, 0, 1
while ys[0] != 0:
q, z = divmod(xs[0], ys[0])
xs, ys = ys, (z, xs[1] - q * ys[1], xs[2] - q * ys[2])
return xs
# 逆数表の生成
def make_inv(n):
inv = [0] * n
for x in range(1, n):
g, y, _ = ext_euclid(x, n)
if g != 1: raise Exception("not prime number!")
inv[x] = (y + n) % n
return inv
class GF:
def __init__(self, n):
self.num = n
self.inv = make_inv(n)
# 正規化
def normalize(self, xs):
xs %= self.num
xs = np.trim_zeros(xs, 'f')
if not len(xs): xs = np.array([0])
return xs
# 足し算
def polyadd(self, xs, ys):
xn = len(xs)
yn = len(ys)
if xn >= yn:
zs = xs.copy()
zs[-yn:] += ys
else:
zs = ys.copy()
zs[-xn:] += xs
return self.normalize(zs)
# 引き算
def polysub(self, xs, ys):
return self.polyadd(xs, self.num - ys)
# 掛け算
def polymul(self, xs, ys):
xn = len(xs)
yn = len(ys)
zn = xn + yn - 1
zs = np.full(zn, 0, dtype=np.int32)
for i in range(yn):
zs[i:i+xn] += ys[i] * xs
return self.normalize(zs)
# 割り算
def polydiv(self, xs, ys):
xn = len(xs)
yn = len(ys)
zs = xs.copy()
qs = []
for _ in range(xn - yn + 1):
temp = (zs[0] * self.inv[ys[0]]) % self.num
zs[:yn] += self.num - temp * ys
qs.append(temp)
zs = zs[1:] % self.num
if qs == []: qs = [0]
return np.array(qs), self.normalize(zs)
# 因数分解
def factorization(self, xs):
ps = []
g = functools.reduce(lambda x, y: math.gcd(x, y), xs)
if g > 1:
xs = xs // g
ps.append((np.array([g]), 1))
m = (len(xs) - 1) // 2 + 1
zs = itertools.product(range(self.num), repeat=m)
for _ in range(self.num): next(zs)
for z in zs:
ys = np.trim_zeros(np.array(z), 'f')
if len(xs) <= len(ys): break
c = 0
while len(xs) > len(ys):
p, q = self.polydiv(xs, ys)
if q[0] == 0:
xs = p
c += 1
else:
break
if c > 0:
if np.all(xs == ys):
c += 1
xs = np.array([1])
ps.append((ys, c))
if len(xs) > 1: ps.append((xs, 1))
return ps