【AbstractAlgebra.jl】Juliaで代数学をやってみたいんじゃ② ～剰余環を使って√2やiをつくろう～

前回記事：【AbstractAlgebra.jl】Juliaで代数学をやってみたいんじゃ①

Introduction

Juliaで代数学をやってみるための便利パッケージAbstractAlgebra.jl
に慣れるために、遊んでみようと思います。今回は、剰余環を使って遊びます。

私なりに確認はしていますが、素人プログラミング＆素人数学のため、本記事は誤りを含む可能性もあります。適宜、公式ドキュメント等、信頼性の高い文献の参照をお願いいたします。

今回の記事の内容についてはTaketo Sanoさんが先駆者です。
Taketo Sanoさんの記事は、数学の解説もすごく充実していているので、必見です！Swiftで数をつくっています。

[Taketo Sanoさんの記事]
Swiftで代数学入門 〜 7. 代数拡大で数を作ろう！
√2をつくる

私のこの記事では、数学の解説は少ないので、ぜひ上記の記事を参照してみてください。

$\sqrt{2}$をつくろう

Juliaで$\sqrt{2}$を書いてみると

julia> sqrt(2)
1.4142135623730951

ですが、これは

julia> sqrt(2)*sqrt(2)
2.0000000000000004

2回掛けても、きっかり$2$にはなってくれないようです。

せっかくなので$\sqrt{2}$をつくってみます。

まず、$\mathbb{Q}[X]/(X^2-2)$をJuliaの世界に出現させます。

using AbstractAlgebra

# SをQ[x]とする
S, x = PolynomialRing(QQ, "x")

# TをQ[x]/(x^2-2)とする
T = ResidueRing(S, x^2-2)

$\mathbb{Q}[X]/(X^2-2)$の元$a,b$を

$$a=X+(X^2-2),b=-X+(X^2-2)$$

とします。この$a,b$もJuliaの世界に出現させてみます。

# fをQ[x]の元 xとする
f = S([0,1])

# gをQ[x]の元 -xとする
g = S([0,-1])

# Q[x]/(x^2-2)におけるfをaとする
a = T(f)

# Q[x]/(x^2-2)におけるgをbとする
b = T(g)

$a^2,b^2$をprintlnで出力してみます。

# a^2を出力
println("a^2=",a^2)

# b^2を出力
println("b^2=",b^2)

出力結果は以下です。

a^2=2
b^2=2

2乗して2になりました！

この項では

$$\mathbb{Q}(\sqrt{2}) \simeq \mathbb{Q}[X]/(X^2-2)$$

であることを使いながら、$a+b\sqrt{2}(a,b \in \mathbb{Q})$と$a+bX+(X^2-2)(a,b \in \mathbb{Q})$を対応させて考えています。
つまり、$\sqrt{2}$と$X+(X^2-2)$が対応し、$-\sqrt{2}$と$-X+(X^2-2)$が対応していると考えて、プログラムを書いています。

そのため、通常想定される実数の$\sqrt{2}$とそっくりそのまま同じものというわけではありません。

この先の項についても、同様の考え方をしています。

$i$をつくろう

Juliaでは複素数が以下のように使えます。

julia> complex(0,1)
0 + 1im

julia> complex(0,1)*complex(0,1)
-1 + 0im

でも、せっかく剰余環が使えるので$i$をつくってみます。

まず、$\mathbb{Q}[X]/(X^2+1)$をJuliaの世界に出現させます。

using AbstractAlgebra

# SをQ[x]とする
S, x = PolynomialRing(QQ, "x")

# TをQ[x]/(x^2+1)とする
T = ResidueRing(S, x^2+1)

$\mathbb{Q}[X]/(X^2+1)$の元$a,b$を

$$a=X+(X^2+1),b=-X+(X^2+1)$$

とします。この$a,b$もJuliaの世界に出現させてみます。

# fをQ[x]の元 xとする
f = S([0,1])

# gをQ[x]の元 -xとする
g = S([0,-1])

# Q[x]/(x^2+1)におけるfをaとする
a = T(f)

# Q[x]/(x^2+1)におけるgをbとする
b = T(g)

$a^2,b^2$をprintlnで出力してみます。

# a^2を出力
println("a^2=",a^2)

# b^2を出力
println("b^2=",b^2)

出力結果は以下です。

a^2=-1
b^2=-1

2乗して-1になりました！

$1$の原始$3$乗根$ω$をつくろう

今度は$1$の原始$3$乗根（よく$ω$で表現されている数）をつくってみます。

まず、$\mathbb{Q}[X]/(X^2+X+1)$をJuliaの世界に出現させます。

using AbstractAlgebra

# SをQ[x]とする
S, x = PolynomialRing(QQ, "x")

# TをQ[x]/(x^2+x+1)とする
T = ResidueRing(S, x^2+x+1)

$\mathbb{Q}[X]/(X^2+X+1)$の元$a,b$を

$$a=X+(X^2+X+1),b=-X-1+(X^2+X+1)$$

とします。この$a,b$もJuliaの世界に出現させてみます。

# fをQ[x]の元 xとする
f = S([0,1])

# gをQ[x]の元 -x-1とする
g = S([-1,-1])

# Q[x]/(x^2+x+1)におけるfをaとする
a = T(f)

# Q[x]/(x^2+x+1)におけるgをbとする
b = T(g)

「$a,b$が$1$（乗法単位元）かどうか」と「$a,b$は同じものかどうか」をチェックし、$a^3,b^3$をprintlnで出力してみます。

println("a=1?：",isone(a))
println("b=1?：",isone(b))
println("a=b?：",a==b)
println("a^3=",a^3)
println("b^3=",b^3)

出力結果は以下です。

a=1?：false
b=1?：false
a=b?：false
a^3=1
b^3=1

3乗して1になりました！

$1$の原始$p$乗根をつくろう

最後です。

私は円分体という代数体が好きなので、$1$の原始$p$乗根（$p$は奇素数）をいっぱいつくって、喜びたいと思います。

まずは円分多項式$X^{p-1}+X^{p-2}+...+X^2+X+1$を生成するための関数を書いてみます。

# 以下2つを返す関数。但しpは奇素数を想定
# - 円分多項式 X^(p-1}+X^(p-2)+...+X+1の係数を格納した配列
# - minN以上maxN以下の素数pを格納した配列
# minN=3 maxN=10の場合、以下を返す。
# (Array[[1, 1, 1, 0, 0, 0, 0, 0], [1, 1, 1, 1, 1, 0, 0, 0], [1, 1, 1, 1, 1, 1, 1]], [3, 5, 7])
# 引数:正の整数minN,maxN
# 戻り値：circlePoly（配列、円分多項式の係数を格納）、primesArray（対象となる素数）
function circlePolyGenerate(minN,maxN)

    circlePoly = Array[]
    coe = Int[]

    i = 1

    primesArray = Primes.primes(minN,maxN)
    maxP = primesArray[end]

    for p in primesArray

        for i in 1:(p)
            push!(coe,1)
            i = i+1
        end
        if p < maxP
            while size(coe)[1]<=maxP
                push!(coe,0)
            end
        end

        push!(circlePoly,coe)

        coe = Int[]
        i = 1

    end

    return (circlePoly,primesArray)

end

以下を$f(X)$とします（$p$は奇素数）。

$$f(X)=X^{p-1}+X^{p-2}+...+X^2+X+1$$

$\mathbb{Q}[X]/(f(X))$の元$a$を以下とします。

$$a = X+(f(X))$$

Juliaの世界に$\mathbb{Q}[X]/(f(X))$と$a,a^2,a^3\cdots a^{p-1}$を出現させてみます。

そして「$a^i (1\leqq i \leqq p-1)$が$1$かどうか」をチェックし、「$a^i (1\leqq i \leqq p-1)$を$p$乗した結果」をprintlnで出力してみます。

今回は$5 \leqq p \leqq 10$とします。

using AbstractAlgebra

minN = 5
maxN = 13

(circlePoly,primesArray) = circlePolyGenerate(minN,maxN)  

n = size(primesArray)[1]

for i in 1:n

    f = circlePoly[i]
    p = primesArray[i]

    # SをQ[x]とする
    S, x = PolynomialRing(QQ, "x")

    # TをQ[x]/(f)とする
    T = ResidueRing(S, S(f))

    # Q(x)におけるxをgとする
    g = S([0,1])

    # Q[x]/(f)におけるgをaとする
    a = T(g)

    println("p=",p)
    println("円分多項式:",S(f))
    
    for j in 1:(p-1)

        b = a^j

        println("a^",j,":",b)
        println("a^",j,"=1?:",isone(b))
        println("a^",j,"のp乗=",b^p,"\n")
    end

    println("\n")

end

出力結果は以下です。

p=5
円分多項式:x^4 + x^3 + x^2 + x + 1
a^1:x
a^1=1?:false
a^1のp乗=1  

a^2:x^2     
a^2=1?:false
a^2のp乗=1  

a^3:x^3     
a^3=1?:false
a^3のp乗=1  

a^4:-x^3 - x^2 - x - 1
a^4=1?:false
a^4のp乗=1



p=7
円分多項式:x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
a^1:x
a^1=1?:false
a^1のp乗=1

a^2:x^2
a^2=1?:false
a^2のp乗=1

a^3:x^3
a^3=1?:false
a^3のp乗=1

a^4:x^4
a^4=1?:false
a^4のp乗=1

a^5:x^5
a^5=1?:false
a^5のp乗=1

a^6:-x^5 - x^4 - x^3 - x^2 - x - 1
a^6=1?:false
a^6のp乗=1



p=11
円分多項式:x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
a^1:x
a^1=1?:false
a^1のp乗=1

a^2:x^2
a^2=1?:false
a^2のp乗=1

a^3:x^3
a^3=1?:false
a^3のp乗=1

a^4:x^4
a^4=1?:false
a^4のp乗=1

a^5:x^5
a^5=1?:false
a^5のp乗=1

a^6:x^6
a^6=1?:false
a^6のp乗=1

a^7:x^7
a^7=1?:false
a^7のp乗=1

a^8:x^8
a^8=1?:false
a^8のp乗=1

a^9:x^9
a^9=1?:false
a^9のp乗=1

a^10:-x^9 - x^8 - x^7 - x^6 - x^5 - x^4 - x^3 - x^2 - x - 1
a^10=1?:false
a^10のp乗=1



p=13
円分多項式:x^12 + x^11 + x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
a^1:x
a^1=1?:false
a^1のp乗=1

a^2:x^2
a^2=1?:false
a^2のp乗=1

a^3:x^3
a^3=1?:false
a^3のp乗=1

a^4:x^4
a^4=1?:false
a^4のp乗=1

a^5:x^5
a^5=1?:false
a^5のp乗=1

a^6:x^6
a^6=1?:false
a^6のp乗=1

a^7:x^7
a^7=1?:false
a^7のp乗=1

a^8:x^8
a^8=1?:false
a^8のp乗=1

a^9:x^9
a^9=1?:false
a^9のp乗=1

a^10:x^10
a^10=1?:false
a^10のp乗=1

a^11:x^11
a^11=1?:false
a^11のp乗=1

a^12:-x^11 - x^10 - x^9 - x^8 - x^7 - x^6 - x^5 - x^4 - x^3 - x^2 - x - 1
a^12=1?:false
a^12のp乗=1

$p$乗して1になりました！

感想

「数学の世界がコンピュータの中に来たぞー！」という気持ちになり、楽しく遊べました！
数学書の演習問題等と連携させつつ、もっと色々書いてみたいな…と思っています。