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> read("c12.gp");
time = 1 ms.
gp> e=ellinit([0,0,0,-36,0])
time = 1 ms.
%2 = [0, 0, 0, -36, 0, 0, -72, 0, -1296, 1728, 0, 2985984, 1728, Vecsmall([1]), [Vecsmall([128, 1])], [0, 0, 0, 0, 0, 0, 0, 0]]
gp> e.disc
%3 = 2985984
gp> factor(e.disc)
%4 =
[2 12]

[3  6]

gp> e.j
%5 = 1728
gp> factor(e.j)
%6 =
[2 6]

[3 3]

gp> ellglobalred(e)
time = 1 ms.
%7 = [576, [1, 0, 0, 0], 16, [2, 6; 3, 2], [[6, -6, 0, 4], [2, -1, 0, 4]]]
gp> factor(576)
%8 =
[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)
time = 1 ms.
%9 = [4, [2, 2], [[-6, 0], [0, 0]]]
gp> elladd(e,[0,0],[6,0])
%10 = [-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> f(x)=x^3-36*x
%11 = (x)->x^3-36*x
gp> poldisc(f(x))
%12 = 186624
gp> factor(f(x))
time = 2 ms.
%13 =
[x - 6 1]

[    x 1]

[x + 6 1]

gp> rr=polroots(f(x))
%14 = [-6.0000000000000000000000000000000000000 + 0.E-38*I, 0.E-38 + 0.E-38*I, 6.0000000000000000000000000000000000000 + 0.E-38*I]~
gp> for(x=-6,0,y2=f(x);y=floor(sqrt(y2));if(y^2==y2,print([x,y])))
[-6, 0]
[-3, 9]
[-2, 8]
[0, 0]

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

[pari/gpによる計算]

gp> t1=[0,0];t2=[6,0];t3=[-6,0];p1=[-3,9];
gp> H=[ellbil(e,p1,p1)/2]
%17 = [0.44431293741980961989709366749034464329]
gp> dd=H[1]-x
%18 = -x + 0.44431293741980961989709366749034464329
gp> ww=polroots(dd)
%19 = [0.44431293741980961989709366749034464329 + 0.E-38*I]~
gp> c1=real(ww[1])
%20 = 0.44431293741980961989709366749034464329

■不等式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])))
%21 = 12.000000000000000000000000000000000000

■不等式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)
%22 = 0.E-38
gp> c3=cc3(c0,e)
%23 = 3.2802535776209728873808981706082474146

■不等式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> he=log(2^6*3^3)
%52 = 7.4547199493640009306891284395166365224
gp> om=omega_pe(-36,0)
%53 = [2.1409010280752311986343110516952903019, 2.1409010280752311986343110516952903019*I]
gp> om1=om[1];om2=om[2];
gp> A0=AA(e,t1,he,om)
%55 = 7.4547199493640009306891284395166365224
gp> A1=AA(e,p1,he,om)
%56 = 7.4547199493640009306891284395166365224
gp> eb=eeb1(e,p1,A0,A1,om)
%57 = 2.4175450998807598383942856716569623999
gp> c4=cc4(1,eb,A0*A1)
%58 = 3305035440.8995543505175522560333864667
gp> c5=cc5(eb)
%59 = 0.88275260375796288835450500835286186183
gp> c6=cc6(c5,he)
%60 = 8.3374725531219638190436334478694983842
gp> c7=cc7(c6,1,2)
%61 = 9.1243697336819091284608655693276749522
gp> c8=cc8(c6,1,2)
%62 = 8.6212857181427988989036316825804802311
gp> g(m)=m^2-(c3/c1+log(4*sqrt(2))/c1+(c4/c1)*(log(m)+c7)*(log(log(m))+c8)^3)
%63 = (m)->m^2-(c3/c1+log(4*sqrt(2))/c1+(c4/c1)*(log(m)+c7)*(log(log(m))+c8)^3)
gp> gdash(m)=2*m-(c4/c1)*((1/m)*(log(log(m))+c8)^3+(log(m)+c7)*3*(log(log(m))+c8)^2*(1/log(m))*(1/m))
%64 = (m)->2*m-(c4/c1)*((1/m)*(log(log(m))+c8)^3+(log(m)+c7)*3*(log(log(m))+c8)^2*(1/log(m))*(1/m))
gp> nn(x)=x-g(x)/gdash(x)
%65 = (x)->x-g(x)/gdash(x)
gp> x=10^100;for(i=1,350,x=nn(x);print(x))
5.0000000000000000000000000000000000000 E99
2.5000000000000000000000000000000000000 E99
1.2500000000000000000000000000000000000 E99
6.2500000000000000000000000000000000000 E98
3.1250000000000000000000000000000000000 E98
......
19636380767.240956451219946465053493003
9818200354.7982995153773413799593199393
4909119541.5423859668697844034252276636
2454597348.1100648161520235537209392225
1227371537.4203682950035147127818712680
613826930.35418118630986796975856430887
307186664.61348084417266289452988306317
154121279.04124427829413574632500708852
78077612.143269774966737678601316771225
40979273.688179286984247938650713625750
24075142.139169544721173671241601144282
18020040.794827197606739757424939079998
16960743.653133269441295002083206115828
16925894.490420676212956229431698077985
16925856.612670965990331244442567482048
16925856.612626215265630373290368932671
16925856.612626215265630310826063921853
16925856.612626215265630310826063921853
16925856.612626215265630310826063921853
16925856.612626215265630310826063921853
16925856.612626215265630310826063921853
16925856.612626215265630310826063921853
16925856.612626215265630310826063921853
16925856.612626215265630310826063921853
......
16925856.612626215265630310826063921853
16925856.612626215265630310826063921853
16925856.612626215265630310826063921853
16925856.612626215265630310826063921853
time = 471 ms.

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

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

[pari/gpによる計算]

gp> K3=16925856.612626215265630310826063921853
%67 = 16925856.612626215265630310826063921853
gp> K1=(4*sqrt(2)/2*real(om1))*exp(c3)
time = 1 ms.
%68 = 160.96728249845753342171921537018673012
gp> K2=c1
%69 = 0.44431293741980961989709366749034464329
gp> KK0(1,2,K3)
%70 = 4583753953138921.2501280065799169098881

K₀=10³⁰とする。 A=[1, 0 ; [K₀*φ(P₁)] K₀] にLLL-algorithmを適用して、reduced basis {b₁,b₂}を求めると、
     b₁ = [ -322709562803153, 178937647369860]^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> K0=10^30
%71 = 1000000000000000000000000000000
gp> a1=floor(K0*phi(e,p1,K0))
%72 = 554485109816868768813402553348
gp> aaa=[1,0;a1,K0]
%73 =
[                             1                               0]

[554485109816868768813402553348 1000000000000000000000000000000]

gp> bbb=qflll(aaa,1)
%74 =
[-322709562803153 -1298627613777081]

[ 178937647369860   720069675036403]

gp> b1n=nr(bbb)
%75 = 368998839525089.17845530743753658571751
gp> b1n-2^(1/2)*2*K3*sqrt(2)
%76 = 368998771821662.72795044637501534241325
gp> M2=(1/K2)*(log(K0*K1)-log(sqrt(b1n^2/2^3-K3^2)-K3))
%77 = 93.755226771696866674466224136201237251
gp> sqrt(M2)
%78 = 9.6827282710864537962132568595201856359

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

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

gp> K3=9
%79 = 9
gp> KK0(1,2,K3)
%80 = 1296.0000000000000000000000000000000000
gp> K0=10^6
%81 = 1000000
gp> a1=floor(K0*phi(e,p1,K0))
%82 = 554485
gp> aaa=[1,0;a1,K0]
%83 =
[     1       0]

[554485 1000000]

gp> bbb=qflll(aaa,1)
%84 =
[-312 -1349]

[ 173   748]

gp> b1n=nr(bbb)
%85 = 356.75341624152669072363536621060608140
gp> b1n-2^(1/2)*2*K3*sqrt(2)
%86 = 320.75341624152669072363536621060608140
gp> M2=(1/K2)*(log(K0*K1)-log(sqrt(b1n^2/2^3-K3^2)-K3))
%87 = 31.815783139351204872084383767823555967
gp> sqrt(M2)
%88 = 5.6405481240169562861515350358336472220

よって、
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.
gp> check1t(5,e,p1,t1)
[12, -36]
[0, 0]
[12, 36]
There are 3 integral points.
time = 1 ms.
gp> check1t(5,e,p1,t2)
[294, -5040]
[-2, 8]
[6, 0]
[-2, -8]
[294, 5040]
There are 5 integral points.
gp> check1t(5,e,p1,t3)
[18, 72]
[-6, 0]
[18, -72]
There are 3 integral points.

よって、楕円曲線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個に限る。

[2020.07.23追記]
自明でない整点(m,n)=(70,24)をピラミッド型に24段積み上げた黄金球で表現すると、以下のようになる。

[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}]の基本単数を使って、解決している。

[2021.03.17追記]
"c1-c8,楕円対数の計算プログラム(pari/gp)"の関数eeb1(),eeb2(),eeb3()の誤りを修正し、計算結果を見直し、訂正した。

[参考文献]

[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: 2021.03.18

H.Nakao