環のパズル
1.「奇数＝平方和」問題
このワークシートはMath by Codeの一部です。
２以上の整数を高々２つの平方数の和に分解する問題、
「２平方数和」問題に取り組んでみよう。
もとの整数が奇素数の場合はどうだろう？
「２より大きい素数ｐがガウス合成元なら、ｐ≡1(mod 4)でm2+n2と分解できる。」
（理由）
ｐは非単元のa,bを使ってp=abと分解でき、形式的にpの2乗をpとｐの共役元の積とみなすこともできる。
p2=pp'=ab(ab)'=aa' bb'=N(a)N(b)となるので、p=N(a)=N(b)に限るね。
a=m+niとおくとN(a)=m2+n2だから、素数ｐが平方数の和になる。
・m,nともに偶数ならp=N(a)≡0+0=0(mod 4)となり、素数pが４の倍数になり、あり得ない。
・mが偶数、nが奇数ならp=N(a)≡0+1=1(mod 4)となる。
（m,nともに奇数ならp≡1+1=2(mod 4)だから、p=2だけ。）
（例）
2=12+12=(1+i)(1-i)
5=22+12=(2+i)(2-i)
13=32+22=(3+2i)(3-2i)
17=42+12=(4+i)(4-i)
29=52+22=(5+2i)(5-2i)
では、整数が合成数の場合はどうだろうか。
整数が合成数ならば、素因数分解できるね。
その素因数それぞれがガウス整数のノルムとしてみなすことができる。
すると、合成数の素因数分解がノルムの積分解になる。
しかも、ノルムの積は１つのガウス整数のノルムになる。つまり、
「ノルム積は合成できる。N(a)N(b)=N(ab)」
（例）
25=5*5=N(2+i)N(2+i)=N(3+4i)=32+42
65=5*13=N(2-i)N(3+2i) から、(2-i)(3+2i)=8+i 。だから、 65=N(8+i)=82+12 5*13=N(1+2i)N(3-2i)から、(1+2i)(3-2i)=7+4i 。だから、65=N(7+4i)=72+42
85=5*17=N(2-i)N(4+i)から、(2-i)(4+i)=9-2i。だから、85=N(9-2i)=92+22
5*17=N(1+2i)N(3-2i)から、(1+2i)(4-i)=6+7i 。だから、85=N(6+7i)=72+62
145=5*29=N(1+2i)N(5-2i)から、(1+2i)(5-2i)=9+8i。だから、145=N(9+8i)=92+82
5*29=N(2+i)N(5-2i)から、(2+i)(5-2i)=12+i 。だから、145=N(12+i)=122+12こうしてみると、1(mod 4)型の素数p,qの積に分解できれば、p,qはN(a),N(b)と置き換えられる。
だから、pq=N(a)N(b)=N(ab)=m2+n2
ガウス整数どうしのおためしの積は次のアプリを使ってみよう。
（一般化）
イメージ的な語り方をしてみる。
ガウス整数のノルム写像Nは、
（複素数である）ガウス整数m+ni→（ふつうの整数である）平方数の和m2+n2
という2世界のパイプ役だ。a=m+ni→通常の整数。
重要なのは、ノルムの積は積のノルム、N(a)N(b)=N(ab)という２世界での積の結果同一性だ。
つまり、Nは準同型写像だ。積を先にしてもあとでNをしてもOKということだね。
素数をすべてガウスノルムを通してガウスの世界に移動すると、
その積が１つのガウス整数にまとめられる。１つのガウス整数はノルムを使ってふつうの整数世界に
m2+n2として返せるのだ。
ガウス世界に移動できない素因数qはそのまま残るから、
素因数が平方数であれば、M2＋N2にできる。
数学的な記述としては、
２より大きい奇数ｐは
・ｐが1(mod 4)型の素数ならｐはガウス合成数だから、ガウスのノルム=m2+n2と表現可能。
・ｐが合成数なら、ｐを構成する素因数qが1(mod 4)型ならガウスのノルムN(a)の置き換えができて、
ノルムの積は積のノルムだから、１つのガウスのノルムに置き換えられるから、m2+n2と表現可能。
・ｐが合成数なら、ｐを構成する素因数ｑが3(mod 4)型ならqはガウス素元で分解きない。ただし、素因数が偶数q2rになっているなら、q2r(m2+n2)=(qrm)2+(qrn)2となり、M2＋N2にできる。
ただし、ガウス整数a=m+niのノルムN(a)は、原点中心の点aの半径２乗だから、aの単元倍,1,-1,i,-i倍しても
変わらない。つまり、N(a)からaを検索すると４通りの可能性がある。
だから、N(a)N(b)...をN(ab...)にするとき、ノルムは同じでも、別のm2+n2の表現になりうる。
つまり、
2 つの整数 x,yを用いて n=x2+y2
⟺ nを素因数分解したときの
4k+3型の素数の指数が全て偶数
2.整数を高々２つの平方和に直す手順
(例）
2025を2平方数和で表そう。
・2025=3*3*3*3*5*5
・3*3*3*3=92
・5*5=N(2+i)N(2-i)=N(3+4i)=32+42
・2025=92(32+42)=(9*3)2+(9*4)2=272+362
この手順を任意の整数でもできるようしたい。
A・素因数分解をする。
B・ｘ＝3(mod 4)型の素因数ｘの指数がすべて偶数かどうかの判定をする。
C・素因数ｘの偶数乗をｙ2の形にするために、ｘと素因数からｙを求める。
D・素因数の型がｘ＝1(mod 4)=N(a)となるガウス整数a=m+niを求める。
E・N(a)N(b)=N(ab)=m2+n2の形にして、ｙ2をかけた、my2+ny2を返す。
＜A。素因数分解＞
質問：1000000以下の数を素因数分解するコードはどうやって作りますか。
大きな方針は、整数Nを1000以下の素数で小さい順に割り続けることです。それだけです。
（このように、車輪から作ることをやりたくない人は、using Primesを先頭にかきましょう。
ただし、Primesをインストールしていないとエラー表示がでますので、
その表示を呼んで、インストールしておきましょう。）
くわしく書くとこうなります。
Nの最初の素因数を見つけるためにNの約数リストを求める関数divisors(N)を作ります。
素因数分解する関数がfactors(N)です。
・最初に1000以下の素数リストprimeを作ります。
・Nの約数リストdivsの2番目が、Nの最小の素因数factorです。
・factorIdはfactorが素数リストprimeの何番目かを登録したものです。
そのあとは繰り返しになります。
NがfactorId番目のprimeで割り切れるなら、primesに登録し、新規ならptypeにも登録します。
Nをprimeで割った商をNとして、factorId番号をふやしていけば、最後はNは1になります。
繰り返しが終わると、
primesには重複ありの素因子が順に入り、
ptypeには重複なしの素因子が入ります。
最後に素因子ごとの指数を求めるために
primesの中の素数の繰り返しをptypeの素因子ごとに集計しましょう。それがresです。
[IN]julia
#===============
function divisors(N)
return filter(m -> N % m ==0, 1:N)
end
function factors(N)
prime =filter(n->length(filter( m -> n % m==0,1:n))==2, 2:1000) #1000以下の素数
prims=[] #重複ありの素因子
ptype=[] #重複なしの素因子
divs = divisors(N)
preP = 0
if N!=1
factor = divs[2] #最小素因子
factorId = findfirst(x ->x ==factor,prime) #最小素因子のインデックス
while N >= prime[factorId]
if N % prime[factorId] == 0
curP = prime[factorId]
push!(prims, curP)
N = N÷curP
if curP != preP
push!(ptype,curP)
preP = curP
end
else
factorId +=1
end
end
res = [length([x for x in prims if item == x]) for item in ptype] #素因数分解の指数
return prims, ptype, res
else
return false
end
end
pArray,pType,idx=factors(2025)
[OUT]
#===========================
(Any[3, 3, 3, 3, 5, 5], Any[3, 5], [4, 2])
geogebraのCASにはFactorsコマンドあります。
これに数値を入れるだけで、素数の種類と種類ごとの指数が行列形式で表示されます。
CASで素因数分解
＜B.平方和に表現できるかどうかを判定する。＞
質問：素因数分解された整数が２平方数和に直せるかどうかの判定をするコードはどうやってかいたらよいでしょうか。
素数の種別リストpListと、種別ごとの指数リストpIndexを渡します。
種別リストを順にみて、もし3(mod 4)型があれば、それの指数が奇数1(mod 4)型なら
フラグを１を立てます。それ以外ではフラグは０です。
１つでもフラグの１があれば、平方和に表現できないと判定します。
まず、日本語をそのまま手続き的に翻訳して、
関数pJudgeを作りましょう。
[IN]julia
#======================
function pJudge(pList,pIndex)
res = 0
for pos in 1:length(pList)
if pList[pos] % 4 ==3
if pIndex[pos] % 2 == 0
res +=0
else
res +=1
end
end
end
return res
end
println(pJudge([3, 5], [5, 2]))
println(pJudge([3, 5], [4, 2]))
#=============
[OUT]
1
0
すごく単純なことをしているわりにかさばっているので、コンパクトな宣言型に
書き換えてみよう。状態変化ではなく、リストを使ってリストを作るという関数型をイメージします。
remは、素数は４で割った剰余にして、指数は２で割った剰余にして、Zipでペアにしたものです。
res01は、(3(mod 4) , 1(mod 2))のときだけ１、それ以外０にしてフラグの合計を出したものです。
行数は減りましたが、わかりやすくはないですね。
[IN]julia
#======================
function pJudge(pList,pIndex)
rem = [(v[1] % 4, v[2] % 2) for v in zip(pList,pIndex)]
res01= sum([(if (x[1]==3 && x[2]== 1) 1 else 0 end) for x in rem])
return res01
end
println(pJudge([3, 5], [5, 2]))
println(pJudge([3, 5], [4, 2]))
#=============
[OUT]
1
0
geogebraのfactorsの結果に絡めるために宣言型でやってみましょう。
CASではなく数式でもFactors関数が使えます。関数が返した行列は２次元配列faの値の受け取り方は
fa(1)が１つめの素数の情報で、fa(1,1)が素数、fa(1,2)が指数です。
そこで、mod(fa(1,1), 4)==3 && mod(fa(1,2),2)==2 がtrueになると、
２平方和ができないフラグを立てます。それ以外をfalseにして合計すれば、
全体としての判定に使えますね。
素数リストのサイズに関係なく支えようにするために、行列の行サイズlength(fa)を得て、
それをsequenceにすることで、１番目の素数から最後まで調べられますよ。
２平方和に分解できるかの判定
＜C.素因数の偶数乗を２乗形にする＞
質問：素因数xpyqzrのp,q,rが偶数のとき、xpyqzr=a2にするにはどうしたら良いでしょうか。
pが偶数のとき、p=2nとかける。ｘp=x2n=(xn)2
ということは、素因数の指数を２で割ったものをかければよいですね。
素因数分解のデータから、計算が必要な素数だけ取り出して指数を半分にした累乗の積を
求めるのがpsqr関数です。
たとえば、
zipで、素数と指数をペアにしたタプルを作り、そのリストをpzipとします。
pzipのうち、素数が3(mod 4)で、指数が0(mod 2)のペアだけ抜き出し、pfiltとします。
その素数の（指数の２分の１）乗をそれぞれ求めてマッピングしたリストにしてから、
リストの要素をすべてかけて数値に次元下げ、reduceしたものがresです。
pfiltにひっかからない素数と指数のペアはこの計算から除外されますね。
[IN]Julia
#====================
function psqr(pList,pIndex)
pzip = [v for v in zip(pList,pIndex)]
pfilt= filter(x -> x[1] % 4==3 && x[2] % 2== 0 ,pzip)
res= reduce(*,map(x-> x[1]^(x[2]÷2),pfilt))
return res
end
println(psqr([3, 5], [5, 2]))
println(psqr([3, 5], [4, 2]))
println(psqr([3, 5, 7], [4, 2, 2]))
println(psqr([3, 5, 7], [4, 90, 4]))
#=======================
[OUT]
1
9
63
441
＜D・素因数の型がp＝1(mod 4)=N(a)となるガウス整数a＞
質問：素数p=1(mod 4)がN(a)となるガウス整数aを求めるコードはどうやってかけますか。
2=12+12=(1+i)(1-i)
5=22+12=(2+i)(2-i)
13=32+22=(3+2i)(3-2i)
17=42+12=(4+i)(4-i)
29=52+22=(5+2i)(5-2i)
p=m2+n2と仮定する。m=nとするとm=sqrt(p/2)となります。
m＜nとするなら、mを１からfoor(sqrt(p/2))の範囲で探索すれば、p=m2+n2となるnが見つかりますね。
[In]Julia
#============
function sum2square(p)
for m in 1:floor(sqrt(p/2))
np = p - m^2
n = floor(sqrt(np))
if n^2 == np
return Int(m), Int(n)
end
end
end
oneMod4 = [5, 13, 17, 29, 37, 41]
print([sum2square(x) for x in oneMod4])
oneMod4 = [5, 5, 13, 29]
print([sum2square(x) for x in oneMod4])
#===========
[OUT]
[(1, 2), (2, 3), (1, 4), (2, 5), (1, 6), (4, 5)][(1, 2), (1, 2), (2, 3), (2, 5)]
a=m+niに対して、ガウスの単元+1,-1,i,-iの４通り倍すると、pをガウス素元分解できますが、
N(a)はすべて、pのままですね。
＜E・N(a)N(b)=N(ab)=m2+n2の形にする＞
質問：複数のガウス素元の積から、そのノルムのもとになる実部と虚部を返すコードはどうかきますか。
実部raelと虚部imagから複素数complex(real,imag)を構成して積を計算しましょう。
ガウス素元の実部と虚部の組のリストprimeListを受取り、複素数リストcListにして返します。
リストのすべての要素をただただ、かけ算した答えをreduce(*,cList)によって出しましょう。
たとえば、
85=5*17=N(2-i)N(4+i)から、(2-i)(4+i)=9-2i。だから、85=N(9-2i)=92+22
のように、9,-2がわかればそれがノルムの計算に使えるので、m,nにあたる実部、虚部だけ返します。
だだし、ガウス素元をノルムで渡すのではなく、ガウス素元にして渡す必要がありますね。
たとえば、5というノルムから、それがノルムになるガウス素元をsum2suareで(1,2)のように
実部と虚部にして渡す必要がありますで、注意してください。
[IN]Julia
#=====================
function NormMulti(primeList)
cList=[complex(x[1],x[2]) for x in primeList]
ret = reduce(*,cList)
return real(ret), imag(ret)
end
oneMod4 = [5, 13, 17, 29, 37, 41]
NormMulti([sum2square(x) for x in oneMod4])
[OUT]
#==========
(5656, 4077)
質問：2より大きな整数n=x^2+y^2となる、x,yを求める手順のコードをまとめましょう。
[IN]Julia
# 2より大きい整数を高々２個の平方数和に直す。
#A・素因数分解をする。
function divisors(N)
return filter(m -> N % m ==0, 1:N)
end
function factors(N)
prime =filter(n->length(filter( m -> n % m==0,1:n))==2, 2:1000) #1000以下の素数
prims =Int[] #重複ありの素因子
ptype=Int[] #重複なしの素因子
divs = divisors(N)
preP = 0
if N!=1
factor = divs[2] #最小素因子
factorId = findfirst(x ->x ==factor,prime) #最小素因子のインデックス
while N >= prime[factorId]
if N % prime[factorId] == 0
curP = prime[factorId]
push!(prims, curP)
N = N÷curP
if curP != preP
push!(ptype,curP)
preP = curP
end
else
factorId +=1
end
end
res = [length([x for x in prims if item == x]) for item in ptype] #素因数分解の指数
return prims, ptype, res
else
return false
end
end
#B・ｘ＝3(mod 4)型の素因数ｘの指数がすべて偶数かどうかの判定をする。
function pJudge(pList,pIndex)
rem = [(v[1] % 4, v[2] % 2) for v in zip(pList,pIndex)]
res01= sum([(if (x[1]==3 && x[2]== 1) 1 else 0 end) for x in rem])
return res01
end
#C・素因数ｘの偶数乗をｙ2の形にするために、ｘと指数からｙを求める。
function psqr(pList,pIndex)
pzip = [v for v in zip(pList,pIndex)]
pfilt= filter(x -> x[1] % 4==3 && x[2] % 2== 0 ,pzip)
res= reduce(*,map(x-> x[1]^(x[2]÷2),pfilt))
return res
end
#D・素因数の型がｘ＝1(mod 4)=N(a)となるガウス整数aを求める。
function sum2square(p)
if p==2
return 1,1
else
for m in 1:floor(sqrt(p/2))
np = p - m^2
n = floor(sqrt(np))
if n^2 == np
return Int(m), Int(n)
end
end
end
end
#E・N(a)N(b)=N(ab)=m2+n2の形にする
function NormMulti(primeList)
cList=[complex(x[1],x[2]) for x in primeList]
ret = reduce(*,cList)
return Int(real(ret)), Int(imag(ret))
end
#AからEまでをつないで、入力の受取りから画面表示までをまとめて実行する。
function DoSum2Sq(N)
Num = factors(N)
factrs = Num[1]
ptype = Num[2]
pindex = Num[3]
if pJudge(ptype,pindex) ==0
print("$(N)は分解できます。")
y = psqr(ptype,pindex)
NumGuessP=filter(x->x % 4!=3 ,factrs)
re,im = NormMulti([sum2square(x) for x in NumGuessP])
m = abs(re*y)
n = abs(im*y)
println(" $N => $factrs ==> $m^2 + $n^2 = $(m^2+n^2)")
return true
else
print("$(N)は２平方和の和に分解できません。")
Num = factors(N)
factrs = Num[1]
println(" $N => $factrs")
return false
end
end
for item in 2:20
DoSum2Sq(item)
end
DoSum2Sq(90)
for item in 0:10
DoSum2Sq(2000 + 5*item)
end
#==================
[OUT]
2は分解できます。 2 => [2] ==> 1^2 + 1^2 = 2
3は２平方和の和に分解できません。 3 => [3]
4は分解できます。 4 => [2, 2] ==> 0^2 + 2^2 = 4
5は分解できます。 5 => [5] ==> 1^2 + 2^2 = 5
6は２平方和の和に分解できません。 6 => [2, 3]
7は２平方和の和に分解できません。 7 => [7]
8は分解できます。 8 => [2, 2, 2] ==> 2^2 + 2^2 = 8
9は分解できます。 9 => [3, 3] ==> 3^2 + 0^2 = 9
10は分解できます。 10 => [2, 5] ==> 1^2 + 3^2 = 10
11は２平方和の和に分解できません。 11 => [11]
12は２平方和の和に分解できません。 12 => [2, 2, 3]
13は分解できます。 13 => [13] ==> 2^2 + 3^2 = 13
14は２平方和の和に分解できません。 14 => [2, 7]
15は２平方和の和に分解できません。 15 => [3, 5]
16は分解できます。 16 => [2, 2, 2, 2] ==> 4^2 + 0^2 = 16
17は分解できます。 17 => [17] ==> 1^2 + 4^2 = 17
18は分解できます。 18 => [2, 3, 3] ==> 3^2 + 3^2 = 18
19は２平方和の和に分解できません。 19 => [19]
20は分解できます。 20 => [2, 2, 5] ==> 4^2 + 2^2 = 20
90は分解できます。 90 => [2, 3, 3, 5] ==> 3^2 + 9^2 = 90
2000は分解できます。 2000 => [2, 2, 2, 2, 5, 5, 5] ==> 44^2 + 8^2 = 2000
2005は分解できます。 2005 => [5, 401] ==> 39^2 + 22^2 = 2005
2010は２平方和の和に分解できません。 2010 => [2, 3, 5, 67]
2015は２平方和の和に分解できません。 2015 => [5, 13, 31]
2020は分解できます。 2020 => [2, 2, 5, 101] ==> 24^2 + 38^2 = 2020
2025は分解できます。 2025 => [3, 3, 3, 3, 5, 5] ==> 27^2 + 36^2 = 2025
2030は２平方和の和に分解できません。 2030 => [2, 5, 7, 29]
2035は２平方和の和に分解できません。 2035 => [5, 11, 37]
2040は２平方和の和に分解できません。 2040 => [2, 2, 2, 3, 5, 17]
2045は分解できます。 2045 => [5, 409] ==> 37^2 + 26^2 = 2045
2050は分解できます。 2050 => [2, 5, 5, 41] ==> 33^2 + 31^2 = 2050