Zenn
🐱

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

2020/12/31に公開
Julia
math
代数学
abstractalgebra
julialang
tech

経緯

Juliaを触りはじめて1ヶ月程。
とある有限群を調べたときに（参照：【初めてのJuliaプログラミング】とある有限群を具体的に求めてみた。）、Juliaで代数学の道具を揃えてみたいと思いました。
AbstractAlgebra.jlが便利そうなので、使い方を備忘録として少しずつ記録していこうと思います。

まずはここから見ると良さそう

英語が苦手な上、プログラミングも初心者すぎて、最初はAbstractAlgebra.jlのトップページを眺めながら途方に暮れていました。
しかし、以下のページを参照することで、徐々に「なんか取っ掛かりが見えてきた…」となってきました。私のように「このドキュメントなんもわからん」な方は、まず最初に読んでみると良いかもしれません。

まずコレを見てみよう！→Constructing mathematical objects in AbstractAlgebra.jl

AbstractAlgebra.jl

AbstractAlgebra.jl
有限体、多項式環、代数体、剰余環など、便利グッズが揃っているようです。苦手な英語を読み解きつつ、簡単そうなところから、じゃんじゃん書いてみます。

Quick Startしてみよう

整数（環）

詳しくはInteger Ring

using AbstractAlgebra

# Rを整数（環）とする
R = ZZ

# 整数（環）における5をaとする
a = R(5)

# 整数（環）の加法単位元
e1 = zero(R)

# 整数（環）の乗法単位元
e2 = one(R)

printlnで出力してみます。

println(R)
println(a)
println(e1)
println(e2)

出力結果は以下です。

Integers
5
0
1

\mathbb{F}_p

詳しくはFinite fields

\mathbb{F}_{13}の場合を書いていきます。

using AbstractAlgebra

# FをZ/13Zとする
F = GF(13)

# Fにおける20をaとする
a = F(20)

# Fにおける39をbとする
b = F(39)

# Fにおける27をcとする
c = F(27)

# Fの加法単位元
e1 = zero(F)

# Fの乗法単位元
e2 = one(F)

# bは加法単位元かどうかチェック
iszero(b)

# cは乗法単位元かどうかチェック
isone(c)

# aの乗法逆元
ai = inv(a)

# cの乗法逆元
ci = inv(c)

# Fの位数
n = order(F)

printlnで出力してみます。

println(F)
println("a=20：",a)
println("b=39：",b)
println("c=27：",c)
println("加法単位元：",e1)
println("乗法単位元：",e2)
println("aの乗法逆元：",ai)
println("cの乗法逆元：",ci)
println("Fの位数：",n)
println("aは加法単位元？：",iszero(a))
println("bは加法単位元？：",iszero(b))
println("aは乗法単位元？：",isone(a))
println("cは乗法単位元？：",isone(c))

出力結果は以下です。

Finite field F_13
a=20：7
b=39：0
c=27：1
加法単位元：0
乗法単位元：1
aの乗法逆元：2
cの乗法逆元：1
Fの位数：13
aは加法単位元？：false
bは加法単位元？：true
aは乗法単位元？：false
cは乗法単位元？：true

\mathbb{Z}/n\mathbb{Z}

詳しくはRing Interface

\mathbb{Z}/12\mathbb{Z}の場合を書いていきます。

using AbstractAlgebra

# RをZ/12Zとする
R = ResidueRing(ZZ, 12)

# Rにおける20をaとする
a = R(20)

# Rにおける36をbとする
b = R(36)

# Rにおける49をbとする
c = R(49)

# Rの加法単位元
e1 = zero(R)

# Rの乗法単位元
e2 = one(R)

# bは加法単位元かどうかチェック
iszero(b)

# cは乗法単位元かどうかチェック
isone(c)

printlnで出力してみます。

println(R)
println("a=20：",a)
println("b=36：",b)
println("c=49：",c)
println("加法単位元：",e1)
println("乗法単位元：",e2)
println("aは加法単位元？：",iszero(a))
println("bは加法単位元？：",iszero(b))
println("aは乗法単位元？：",isone(a))
println("cは乗法単位元？：",isone(c))

出力結果は以下です。

Residue ring of Integers modulo 12
a=20：8
b=36：0
c=49：1
加法単位元：0
乗法単位元：1
aは加法単位元？：false
bは加法単位元？：true
aは乗法単位元？：false
cは乗法単位元？：true

R[X]/(f(X))

\mathbb{Q}[X]/(X^2+1)の場合を書いていきます。

using AbstractAlgebra

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

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

# Sの元aをx^4+3とする
a = S([3,0,0,0,1])

# Sの元bをx^3+3x^2+2x+1とする
b = S([1,2,3,1])

# Sの元cを-(x+2)/5とする
c = S([-2//5,-1//5])

# Tにおけるaをarとする
ar = T(a)

# Tにおけるbをbrとする
br = T(b)

# Tにおけるcをcrとする
cr = T(c)

# dr=br*cr
dr = br*cr

# brの乗法逆元
bri = inv(br)

# crの乗法逆元
cri = inv(cr)

printlnで出力してみます。

println("S：",S)
println("T：",T)
println("a：",a)
println("b：",b)
println("c：",c)
println("T(a)：",ar)
println("T(b)：",br)
println("T(c)：",cr)
println("T(d)*T(c)：",dr)
println("T(b)の乗法逆元：",bri)
println("T(c)の乗法逆元：",cri)

出力結果は以下です。

S：Univariate Polynomial Ring in x over Rationals
T：Residue ring of Univariate Polynomial Ring in x over Rationals modulo x^2 + 1
a：x^4 + 3
b：x^3 + 3*x^2 + 2*x + 1
c：-1//5*x - 2//5
T(a)：4
T(b)：x - 2
T(c)：-1//5*x - 2//5
T(d)*T(c)：1
T(b)の乗法逆元：-1//5*x - 2//5
T(c)の乗法逆元：x - 2

フェルマーの小定理

AbstractAlgebra.jlを使いながらフェルマーの小定理の具体例を書いてみます。

using AbstractAlgebra

#引数：整数a,素数p
#戻り値：F_pでのa^{p-1}
function fermatsSmallpower(a,p)

    # F = F_pとする
    F = GF(p)

    # F_pにおけるaをarとする
    ar = F(a)

    return ar^(p-1)

end

p=13として具体例を計算して、出力してみます。

p = 13

for i in 0:25
    println(i,"^(p-1)=",fermatsSmallpower(i,p))
end

出力結果は以下です。

0^(p-1)=0
1^(p-1)=1
2^(p-1)=1
3^(p-1)=1
4^(p-1)=1
5^(p-1)=1
6^(p-1)=1
7^(p-1)=1
8^(p-1)=1
9^(p-1)=1
10^(p-1)=1
11^(p-1)=1
12^(p-1)=1
13^(p-1)=0
14^(p-1)=1
15^(p-1)=1
17^(p-1)=1
18^(p-1)=1
19^(p-1)=1
20^(p-1)=1
21^(p-1)=1
22^(p-1)=1
23^(p-1)=1
24^(p-1)=1
25^(p-1)=1

平方剰余

AbstractAlgebra.jlを使いながら平方剰余かどうかをチェックするコードを書いてみます。

using AbstractAlgebra

#引数：整数a,素数p
#戻り値：aがpを法として平方剰余ならば1,平方非剰余ならば-1を返す
function quadraticResidue(a,p)
    
    F = GF(p)

    ar = F(a)

    for i in 0:(p-1)
        ir = F(i)

        if ir^2 == ar
            return 1
        end	
    end

    return -1

end

2<p<100,a=-1の場合を計算して、結果を出力してみます（第1補充法則の具体例の計算です）。
前半がp \equiv 1 \pmod{4}の場合、後半がp \equiv 3 \pmod{4}の場合です。

using Primes

onePrimes = Int[]
threePrimes = Int[] 

for p in Primes.primes(3,100)
    if  mod(p,4) == 1
        push!(onePrimes,p)
    else mod(p,4) == 3
        push!(threePrimes,p)
    end 
end


for p in onePrimes
    quadCheck = quadraticResidue(-1,p)
    println(p,"：",quadCheck)
end

for p in threePrimes
    quadCheck = quadraticResidue(-1,p)
    println(p,"：",quadCheck)
end

出力結果は以下です。
-1がpを法として平方剰余ならば1となり、-1がpを法として平方非剰余ならば-1となっています。

5：1
13：1
17：1
29：1
37：1
41：1
53：1
61：1
73：1
89：1
97：1
3：-1
7：-1
11：-1
19：-1
23：-1
31：-1
43：-1
47：-1
59：-1
67：-1
71：-1
79：-1
83：-1

感想

今回扱った例（フェルマーの小定理、平方剰余）に関しては、AbstractAlgebra.jlを使わずともサクサク書けちゃうと思います。
しかし、AbstractAlgebra.jlを使うことで、頭の中にある\mathbb{F}_pをそのままコードに落とし込める感じになり、ストレス少な目でプログラミングできました。
また、pで割って余りを求めて…というのを書かなくて良いので楽でした。

Discussion