Integral Points on Elliptic Curve: n(n+1)(2n+1)/6=m^2, Rational Points and Integral Points on Elliptic Curve: y^2=x^3-36x

[2003.04.12]n(n+1)(2n+1)/6=m^2の整点, y^2=x^3-36xの有理点と整点

c1-c8,楕円対数の計算プログラム(pari/gp)

■問題: 1²からn²までの平方数の和が平方数m²になるものを全て求めよ。---- 参考文献[5]のExample 1(The Square Pyramid Problem)による。

不定方程式
     1²+2²+3²+...+n² = m² -------- (1)
の正整数解(m,n)を求める。
つまり、楕円曲線
     C: n(n+1)(2n+1)/6 = m²
の整点(m,n)を求める。
     6m² = 2n³+3n²+n

     X = 12n+6,
     Y = 72m

とすると、楕円曲線
     E: Y² = X³-36X ------- (2)
となる。
m,nが正整数なら、X,Yも正整数になるので、楕円曲線Eの整点を求めれば十分である。
最初に、楕円曲線Eの有理点を求める。

■楕円曲線Eの判別式 Δ,j-不変量 j,conductor Nは、それぞれ、
     Δ = 2985984 = 2¹²・3⁶
     j = 1728 = 2⁶・3³
     N = 576 = 2⁶・3²
である。
[pari/gpによる計算]

gp>  e=ellinit([0,0,0,-36,0])
time = 71 ms.
%1 = [0, 0, 0, -36, 0, 0, -72, 0, -1296, 1728, 0, 2985984, 1728, [6.000000000000000000000000000, 0.E-28, -6.000000000000000000000000000]~, 1.070450514037615599317155525, 1.070450514037615599317155525*I, -1.467416107700331192416547425, -4.402248323100993577249642277*I, 1.145864303003395471248349296]
gp>  e.disc
time = 0 ms.
%2 = 2985984
gp>  factor([e.disc])
time = 0 ms.
%3 = [[2, 12; 3, 6]]
gp>  e.j
time = 0 ms.
%4 = 1728
gp>  factor([e.j])
time = 0 ms.
%5 = [[2, 6; 3, 3]]
gp>  ellglobalred(e)
time = 9 ms.
%6 = [576, [1, 0, 0, 0], 16]
gp>  factor([576])
time = 0 ms.
%7 = [[2, 6; 3, 2]]

■楕円曲線Eは、ねじれ点群Z/2Z×Z/2Zを持つ。
楕円曲線Eは、位数2のねじれ点 T₁ = (0,0), T₂ = (6,0), T₁+T₂ = (-6,0)を持ち、Eのねじれ点群は、
E(Q)_tors = {(0,0), (6,0), (-6,0), O}
である。
[pari/gpによる計算]

gp>  elltors(e,1)
time = 110 ms.
%8 = [4, [2, 2], [[0, 0], [6, 0]]]
gp>  elladd(e,[0,0],[6,0])
time = 0 ms.
%9 = [-6, 0]

■楕円曲線EのModell-Weil群のrankは、1であることが分かる。
Cremonaのmwrank3によって、このようにMordell-Weil群の生成元P₁(-3,9)が見つかる。よって、
     E(Q) = Z×Z/2Z×Z/2Z
である。
これは、6が合同数(congruent number)であることを示している。

■楕円曲線Eの全ての整点を、[4]の方法によって求める。
u=1,v=0とする。
以下の形式
     E': y² = f(x) ----- (3)
     f(x) = x^3-36x = x(x+6)(x-6) -------- (4)
を得る。
3次方程式 f(x) = 0の最大の実数根をγとする。 ------ (*)

f(x)=x³-36x=0の根は、３個の実数(有理数)であり、
     γ'' = -6,
     γ' = 0,
     γ = 6
である。
任意のP ∈ E(Q)に対して、有理整数m₁,m₀(ただし、m₀=0,1), ねじれ点 T ∈ E(Q)_torsが存在して、
     P = m₁P₁+m₀T ------ (5)
となる。

(X,Y)がEの整点であるとする。
f(x) > 0 iff (x > γ)または(γ'' < x < γ')

γ'' < x < γ'なる整点(x,y)は、
     (-6, 0), (-3, ±9), (-2, ±8), (0, 0)
の6個である。

よって、x > γの範囲で、Eの整点(x,y)を求める。

gp>  poldisc(x^3-36*x)
time = 0 ms.
%11 = 186624
gp>  rr=polroots(x^3-36*x)
time = 9 ms.
%12 = [-6.000000000000000000000000000 + 0.E-28*I, 0.E-28 + 0.E-28*I, 6.000000000000000000000000000 + 0.E-28*I]~
gp> for(x=-6,0,y2=x^3-36*x;y=floor(sqrt(y2));if(y^2==y2,print([x,y])))[-6, 0]
[-3, 9]
[-2, 8]
[0, 0]
time = 19 ms.

■不等式1
P ∈ E(Q)が(5)で表現されるとすると、
h^(P) >= c₁|m₁| ------- (7)
である。

[pari/gpによる計算]

gp>  read("c1.gp")
time = 14 ms.
gp>  t1=[0,0];t2=[6,0];t3=[-6,0];p1=[-3,9];
time = 0 ms.
gp> H=[ellbil(e,p1,p1)/2]
time = 17 ms.
%20 = [0.4443129374198096198970936673]
gp> ww=polroots(dd);
time = 3 ms.
%21 = [0.4443129374198096198970936673 + 0.E-28*I]~
gp> c1=real(ww[1]);
time = 0 ms.

■不等式2
c₂ = 2*max{|γ|,|γ'|,|γ''|}
とする。任意のx >= c₂に対して、
|∫_x^∞(dt/sqrt{f(t)})| <= 4*sqrt(2)|x|^-1/2 --------- (8)
となる。

[pari/gpによる計算]

gp>  c2=2*max3(abs(real(rr[1])),abs(real(rr[2])),abs(real(rr[3])))
time = 0 ms.
%24 = 12.00000000000000000000000000

■不等式3
u=1,v=0とする。
X₀を正の有理整数で、X₀ > vとする。
このとき、任意のP ∈ E(Q), X(P) >= X₀に対して、
x(P) > 0,
h^(P)-(1/2)*log x(P) <= c₃ ------ (9)
となる。

[pari/gpによる計算]

gp>  c0=log(1)
time = 0 ms.
%25 = 0.E-28
gp>  c3=cc3(c0,e)
time = 2 ms.
%26 = 3.280253577620972887380898170

■不等式1,2,3より、
     |φ(P)|=|(1/ω)∫_x(P)^∞(dt/sqrt{f(t)})| <= (4*sqrt(2)/ω)|x(P)|^-1 <= (4*sqrt(2)/ω)exp(c₃-c₁M²) ----------- (10)
を得る。

■主不等式
P ∈ E₀(R) iff x(P) >= γ iff X(P) >= u²γ+v なので、不等式2,3を満たすために、
     X₀ = floor(max{c₂,u²γ+v,v})+1
とする。

     M² < c₃c₁^-1-c₁^-1log(4*sqrt(2))+c₄c₁^-1(log M+c₇)(log log M+c₈)^r+2 ------- (11)
となる。

[pai/gpによる計算]

gp>  X0=floor(max3(c2,real(rr[3]),0))+1
time = 0 ms.
%32 = 13
gp> factor([e.j])
time = 0 ms.
%33 = [[2, 6; 3, 3]]
gp> he=log(2^6*3^3)
time = 0 ms.
%34 = 7.454719949364000930689128439
gp> A0=AA(e,t1,he)
time = 2 ms.
%35 = 7.454719949364000930689128439
gp>  A1=AA(e,p1,he)
time = 15 ms.
%36 = 7.454719949364000930689128439
gp>  eb=eeb1(e,p1,A0,A1)
time = 3 ms.
%37 = 2.258437048866518060526713459
gp>  c4=cc4(1,eb,A0*A1)
time = 0 ms.
%38 = 4936902317.050147879454625699
gp>  c5=cc5(eb)
time = 0 ms.
%39 = 0.8146730027404535252170846649
gp>  c6=cc6(c5,he)
time = 0 ms.
%40 = 8.269392952104454455906213104
gp>  c7=cc7(c6,1,2)
time = 0 ms.
%41 = 9.056290132664399765323445225
gp>  c8=cc8(c6,1,2)
time = 0 ms.
%42 = 8.553206117125289535766211339
gp> g(m)=m^2-(c3/c1+log(4*sqrt(2))/c1+(c4/c1)*(log(m)+c7)*(log(log(m))+c8)^2)
time = 0 ms.
gp> gdash(m)=2*m-(c4/c1)*((1/m)*(log(log(m))+c8)^2+(log(m)+c7)*2*(log(log(m))+c8)*(1/log(m))*(1/m))
time = 0 ms.
gp> nn(x)=x-g(x)/gdash(x)
time = 0 ms.
gp> x=10^100;for(i=1,320,x=nn(x);print(x))
5.000000000000000000000000000 E99
2.500000000000000000000000000 E99
1.250000000000000000000000000 E99
6.250000000000000000000000000 E98
3.125000000000000000000000000 E98
......
19830718.23629249828764692002
10831790.44853310658556568901
7051695.986205113153231359330
6009174.458720468743436283556
5914535.899365930348371151625
5913738.914333713026544268291
5913738.857774221961236328325
5913738.857774221676385089380
5913738.857774221676385089380
5913738.857774221676385089380
5913738.857774221676385089380
5913738.857774221676385089380
5913738.857774221676385089380
time = 2,091 ms.

主不等式(11)より、
     M < 5913738.857774221676385089380
であることが分かる。
しかし、ここで求めたMの上限は大き過ぎるので、LLL-algorithmによって、上限を下げる。

■不等式(10),(11)を単純に、
     |φ(P)| < K₁exp(-K₂M²), M < K₃ ------- (12)
と書くことができる。
ここで、
     K₁ = (4*sqrt(2)/ω)*exp(c₃),
     K₂ = c₁,
     K₃ = 5913738.857774221676385089380
である。

[pari/gpによる計算]

gp>  K3=5913738.857774221676385089380
time = 0 ms.
%48 = 5913738.857774221676385089379
gp>  K1=(4*sqrt(2)/2*real(e.omega[1]))*exp(c3)
time = 0 ms.
%49 = 80.48364124922876671085960767
gp>  K2=c1
time = 0 ms.
%50 = 0.4443129374198096198970936673
gp>  KK0(1,2,K3)
time = 0 ms.
%51 = 559556916447180.0971638365428

K₀=10²⁰とする。 A=[1, 0 ; [K₀*φ(P₁)] K₀] にLLL-algorithmを適用して、reduced basis {b₁,b₂}を求めると、
     b₁ = [ -531939361, 294952455]^t
となる。このとき、
     ||b₁|| >= 2^(1/2)*2*K₃*sqrt(2) = 4*K₃ ----- (16)
ならば、
     M² <= K₂^-1(log(K₀K₁) - log(2^-3*||b₁||²/2²-K₃²)-K₃) ------ (17)
が成立する。

[pari/gpによる計算]

gp>  default(realprecision,30)
   realprecision = 38 significant digits (30 digits displayed)
time = 0 ms.
gp>  ec=ellinit([0,0,0,-36,0])
time = 5 ms.
%52 = [0, 0, 0, -36, 0, 0, -72, 0, -1296, 1728, 0, 2985984, 1728, [6.00000000000000000000000000000, 0.E-38, -6.00000000000000000000000000000]~, 1.07045051403761559931715552584, 1.07045051403761559931715552584*I, -1.46741610770033119241654742581, -4.40224832310099357724964227743*I, 1.14586430300339547124834929663]
gp>  K0=10^20
time = 0 ms.
%53 = 100000000000000000000
gp>  a1=floor(K0*phi(ec,p1,K0))
time = 3 ms.
gp>  aaa=[1,0;a1,K0]
time = 0 ms.
%55 = 
[1 0]

[55448510981686876881 100000000000000000000]

gp>  bbb=qflll(aaa,1)
time = 14 ms.
%56 = 
[-531939361 -10252797340]

[294952455 5685023459]

gp>  b1n=nr(bbb)
time = 0 ms.
%57 = 608240441.348333353321588719356
gp>  b1n-2^(1/2)*2*K3*sqrt(2)
time = 0 ms.
%58 = 584585485.917236466616048361837
gp>  M2=(1/K2)*(log(K0*K1)-log(sqrt(b1n^2/2^3-K3^2)-K3))
time = -9 ms.
%59 = 70.40458927575631660219646764
gp>  sqrt(M2)
time = 0 ms.
%60 = 8.390744262325977755207277353

よって、
M <= 8
が得られた。

次に、K₃=8, K₀=10⁶とする。 A=[1, 0 ; [K₀*φ(P₁)] K₀] に、再度LLL-algorithmを適用する。
[pari/gpによる計算]

gp>  K3=8
time = 0 ms.
%61 = 8
gp>  KK0(1,2,K3)
time = 0 ms.
%62 = 1024.00000000000000000000000000
gp> KK0(1,2,K3))
time = 0 ms.
%65 = 1024.00000000000000000000000000
gp>  K0=10^6
time = 0 ms.
%66 = 1000000
gp>  a1=floor(K0*phi(ec,p1,K0))
time = 1 ms.
%67 = 554485
gp>  aaa=[1,0;a1,K0]
time = 0 ms.
%68 = 
[1 0]

[554485 1000000]

gp>  bbb=qflll(aaa,1)
time = 0 ms.
%69 = 
[-312 -1349]

[173 748]

gp>  b1n=nr(bbb)
time = 0 ms.
%70 = 356.753416241526690723635366210
gp>  b1n-2^(1/2)*2*K3*sqrt(2)
time = 0 ms.
%71 = 324.753416241526690723635366210
gp>  M2=(1/K2)*(log(K0*K1)-log(sqrt(b1n^2/2^3-K3^2)-K3))
time = 0 ms.
%72 = 30.23526460870034559191255606
gp>  sqrt(M2)
time = 0 ms.
%73 = 5.498660255798710973121817002

よって、
M <= 5
が得られた。

■|m₁| <= 5, m₀=0,1, T ∈ E(Q)_torsについて、P=m₁P₁+m₀Tが整点かどうかを確認する。

[pari/gpによる計算]

gp>  check1(5,e,p1)
[-3, -9]
[-3, 9]
There are 2 integral points.
time = 6 ms.
gp>  check1t(5,e,p1,t1)
[12, -36]
[0, 0]
[12, 36]
There are 3 integral points.
time = 5 ms.
gp>  check1t(5,e,p1,t2)
[294, -5040]
[-2, 8]
[6, 0]
[-2, -8]
[294, 5040]
There are 5 integral points.
time = 4 ms.
gp>  check1t(5,e,p1,t3)
[18, 72]
[-6, 0]
[18, -72]
There are 3 integral points.
time = 3 ms.

よって、楕円曲線Eの整点(X,Y)は、
     (±6, 0), (-3, ±9), (-2, ±8), (0, 0), (12,±36), (18, ±72), (294,±5040)
の13個に限る。

■最後に、楕円曲線Cの整点(m,n)を求める。
Xは6の倍数,Yは72の倍数でなければならないので、楕円曲線Cの整点(m,n)は、
     (0, 0), (0, -1), (±1, 1), (±70, 24)
の6個に限る。

よって、(1)の解で、m,nが正整数になるものは、
     (1, 1), (70, 24)
の2個に限る。

[2004.08,22追記]
■SIMATH-4.6(simcalc)を使って、楕円曲線Eの整点を求める。
simcalcによると、楕円曲線EのMordell-Weil群E(Q)のrankは1であり、その生成元は(-3,9)である。
楕円曲線Eの整点を(Pと-Pを同一視して)求めると、
     (-6,0), (-3,9), (-2,8), (0,0), (6,0), (12,36), (18,72), (294,5040)
の8個となる。
この中で(-6,0),(0,0),(6,0)は位数2のねじれ点であり、それ以外の5点のそれぞれを(-1)倍した点も整点であるので、楕円曲線Eは合計3+5*2=13個の整点を持つ。

[simcalcによる計算]

> E=EC(0,0,0,-36,0)
                   simcalc in free(): warning: junk pointer, too high to make sense.
         E = EC(-36, 0)
> basismwg(E)
  basis :  PT(-3, 9, 1)
simcalc in free(): warning: junk pointer, too high to make sense.
         @ = 1
> faintp(E)
  all nontrivial integral points modulo negation :
  PT(-6, 0, 1) = PT(-6, 0, 1) + 0*PT(-3, 9, 1)
  PT(0, 0, 1) = PT(0, 0, 1) + 0*PT(-3, 9, 1)
  PT(6, 0, 1) = PT(6, 0, 1) + 0*PT(-3, 9, 1)
  PT(-3, 9, 1) = PT(0, 1, 0) + PT(-3, 9, 1)
  PT(18, 72, 1) = PT(-6, 0, 1) - PT(-3, 9, 1)
  PT(12, 36, 1) = PT(0, 0, 1) + PT(-3, 9, 1)
  PT(-2, 8, 1) = PT(6, 0, 1) - PT(-3, 9, 1)
  PT(294, 5040, 1) = PT(6, 0, 1) + 2*PT(-3, 9, 1)
 
simcalc in free(): warning: junk pointer, too high to make sense.
         @ = PT(-6, 0, 1)

[2007.09.17追記]
■不定方程式(1)の別の解法が文献[6](p.424-427)に記載されている。
mが偶数の場合、不定方程式y^2=8x^4+1に帰着している。
mが奇数の場合、数列{M_n},{G_n}
α=2+sqrt{3}, β=2-sqrt{3},
M_n=(αⁿ+βⁿ)/2, G_n=(αⁿ-βⁿ)/(α-β)
および平方剰余、２次整数環Z[sqrt{3}]の基本単数を使って、解決している。

[参考文献]

[1]Joseph H. Silverman, John Tate(著), 足立恒雄, 木田雅成, 小松啓一, 田谷久雄(訳), "楕円曲線論入門", シュプリンガー・フェアラーク東京, 1995, p127, ISBN4-431-70683-6, {3900円}.
[2]Joseph H. Silverman, "The Arithmetic of Elliptic Curves", GTM 106, Springer-Verlag New York Inc., 1986, p59-60, ISBN0-387-96203-4.
[3]Nigel P. Smart, "The Algorithmic Resolution of Diophantine Equations", LMSST 41, Cambridge University Press, 1998, ISBN0-521-64633-2.
[4]R. J. Stroeker, N. Tzanakis, "Solving elliptic diophantine equations by estimating linear forms in elliptic logarithms", Acta. Math., 3, p209-220, 1994.
[5]John Cremona, "The Arithmetic of Elliptic Curves", University of Exeter.
[6]Henri Cohen, "Number Theory Volume I: Tools and Diophantine Equations", Springer Science+Business Media, LLC, 2007, p424-427, ISBN978-0-387-49922-2.

Last Update: 2007.09.17

H.Nakao