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

7 min read読了の目安(約6700字

前回記事:【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,b1(乗法単位元)かどうか」と「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になりました!

感想

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