pari-gpで初等数論

2008/03/13

富山大学理工学研究部(理学) 木村巌

iwao@sci.u-toyama.ac.jp

この文書では,pari-gpを用いて,初等数論(整数の剰余から、素体の平方剰余まで)の計算を具体的に行う.証明や詳しい解説については,例えば,小野孝「数論序説」,裳華房,を見よ.



行頭に gp> とあるのが,gpのプロンプトである.\\ 以下行末までコメントである.

整数の商と余り

a, bが整数の時,b = qa + r, (0≦r<a)と一通りに書ける.有理整数環ZEuclid性.

gpqを実際に求めるには,b\a, rを求めるにはb%a.

gp> a = 3; b = 5; \\ 5 = 1*3 + 2
gp> b \ a
1
gp> b % a
2

ZEuclid整域なので,単項イデアル整域(PID),更に素元分解一意整域(UDF).

与えられた整数nが素数であることを証明するのは,isprime(n). 桁数が増えると(1,000桁程度)大変遅い.与えられた整数nが強擬素数であることを示すのがispseudoprime(n).

gp> isprime(257)
1

n番目の素数を与えるのがprime(n). 始めのn個の素数からなるベクトルを返す

のがprimes(n). n以下の素数の個数を返すのがprimepi(n). primepi()素朴な実装なので,nは現在の最大素数primelimitを越えることは出来ない.この値はdefault(primelimit)で見ることが出来,default(primelimit, 16777216)のように再設定できる.

gp> prime(25)
97
gp> primes(25)
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
gp> primepi(100)
25

素数についてのループをするために,forprime()が用意されている.

gp> forprime(p=2, 10^5, if(Mod(2, p^2)^(p-1) == Mod(1, p^2), print(p))

これはWieferich素数1093, 3511を表示する.

自然数nの素因子分解を行うのがfactorint(n). ECM, MPQS, Pollard-Brent Rho, Shanks's SQUFOFなどを組み合わせて行う.整数や多項式などの因子分解はfactor(n)でできる.

gp> factorint(2^32 + 1) \\ 2^32+1 = 641*6700417
[641 1]

[6700417 1]

最大公約数はgcd(a, b), 最小公倍数はlcm(a, b).

gp> gcd(a, b)
1
gp> lcm(a, b)
21

gcdを自分で書くと次のようになるだろう:

gp> mygcd(a, b) = {local(c); if(a>b,c=b;b=a;a=c);
	if(a==1 || b%a == 0, return(a), mygcd(b%a, a))}

lcmは次のようになる:

gp> mylcm(a, b) = a*b/mygcd(a,b)

整数s, tが存在して,as + bt = gcd(a, b) が成り立つ.s, tを求めるには,bezout(a,b)をつかう.Bezoutは人名で,ベズーと読む.

gp> bezout(a, b)
[2, -1, 1] \\ s = 2, t = -1, gcd = 1, つまり 2*3 + (-1)*5 = 1.

有理整数環の剰余環

有理整数環Zを整数nを法として類別した剰余環Z/nZの元はMod(a, n)で表される.加減乗算は通常通りできる.anが互いに素なら,逆元も計算できる:

gp> Mod(1,5) + Mod(2, 5)
Mod(3, 5)
gp> Mod(2, 5)*Mod(3, 5)
Mod(1, 5)
gp> Mod(1,5)/Mod(2, 5)
Mod(3, 5)
gp> Mod(2, 5)^4 \\ Fermatの小定理から,Mod(1, 5)になるはず
Mod(1, 5)

中国式剰余定理

x = a (mod m), x = b (mod n)なるxを,m, nが互いに素な時に求めるのが

chinese(Mod(a, m), Mod(b, n)):

gp> chinese(Mod(2,3),Mod(3, 5))
Mod(8, 15)

自分で書くと,次の用になるだろう:

gp> mychinese(am,bn)={local(a,b,m,n,bez,i); 
    a=component(am,2); b=component(bn,2); 
    m=component(am,1); n=component(bn,1);
    bez=bezout(m,n);
    if(bez[3]!=1,error("moduli are not mutually prime"));
    return(Mod(b*bez[1]*m+a*bez[2]*n, m*n))}

am=Mod(a, m)のとき,component(am, 2)aを,component(am, 1)mを与える.

条件が3つ以上あるときは,入れ子にする:

gp> chinese(chinese(Mod(2,3),Mod(3, 5)), Mod(5,7))
Mod(68, 105)

剰余環の乗法群の構造

nの剰余環Z/nZの乗法群(Z/nZ)^*の構造は,znstar(n)で知ることができる.

返り値は3成分のベクトル[s, t, u]で,sは乗法群(Z/nZ)^*の位数,よって

Euler関数φ(n)の値と等しい(これはeulerphi(n)でも求まる).tは乗法群

(Z/nZ)^*を巡回成分で書いたときの,各巡回成分の位数からなるベクトル.

uは各巡回成分の生成元からなるベクトル.



gp> znstar(15) \\ (Z/15Z)^* は,
[8, [4, 2], [Mod(8, 15), Mod(11, 15)]] \\ 位数8, Z/4ZZ/2Zに同型で,
それぞれの巡回成分の生成元は8 mod 15, 11 mod 15.

nが素数pなら,Z/pp個の元からなる体である.

よって,Z/pの乗法群(Mod(0, p)以外の元全体)は巡回群である.

その生成元を,mod pの原始根というのだった.原始根を求める関数はznprimroot(n).

gp> znprimroot(5)
Mod(2, 5)
znprimroot(n)n4もしくは奇素数の冪に対して結果を返す(この時のみ(Z/nZ)^*は巡回群).

znprimroot(p)を自分で書くと,次のようになるだろう.ただし,次の例は簡単に書くことを主眼にしており,効率は全く度外視している.

gp> myznprimroot(p)={local(g);
	for(g=2,p-1,if(length(Set(vector(p-1,i,Mod(g,p)^i))) == p-1,
	return(Mod(g,p))))}

znprimroot(n)の返り値をgとすると,Mod(a, n)は全てgの巾で表される:

Mod(a,n) = g^j for some j. jを求めるのが(つまり,(Z/nZ)^*で離散対数問題を解くのが)znlog().

gp> g = znprimroot(97)
Mod(5, 97)
gp> znlog(Mod(60,97), g)
43
gp> g^43
Mod(60,97)

位数を求めるのはznorder(Mod(a,n)). これらも,n4もしくは奇素数の冪に対して結果を返す.

平方剰余

整数aが法pについて平方剰余とは,x^2 = a (mod p)となる整数xが存在するこ

と.Legendre記号(x/p)は,pxを割り切るとき9, 整数aが法pについて平方剰

余の時1, それ以外の時-1を返す.乗法的に拡張して,Kronecker記号が定義さ

れる.gpではkronecker(a, b).

gp> kronecker(-1, 11)
-1
gp> kronecker(2, 17)
1

p素数として,p*を,p4を法として13かに応じて,p, -pと定義する.

2次体Q(√p*)が円のp分体に含まれることを数値的に検証してみよう.Gaussg(p)=Σ(n/-p)ζ^n, ζ1の原始p乗根,を

gp> g(p) ={local(c);c=polcyclo(p);sum(n=1,p-1,kronecker(n,-p)*Mod(x,c)^n)}

と定義する.polcyclo(p)p次円分多項式を与える.

g(p)の平方g(p)^2が,p4を法として13かに応じて,p, -pになる(つまりg(p)^2=p*, よって√p* = ±g(p)∈p分体)ことを見るのだった.

今の場合,1の原始p乗根ζMod(x, c), c = polcyclo(p)である.

gp> forprime(p=3,100,print(p, ", ",component(g(p)^2,2)))

3, -3

5, 5

7, -7

11, -11

13, 13

17, 17

19, -19

23, -23

29, 29

となっている.

pari/gpの頁に戻る