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

[2003.05.11]n(n+1)(n+2)/6=m(m+1)/2の整点, y²=x³-144x+1296の有理点と整点

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

■問題Ａ: 平方数m²で、1からn(n+1)/2までの連続する三角数の和に等しくなるものを全て求めよ。
不定方程式
     (1・2)/2+(2・3)/2+(3・4)/2+...+n(n+1)/2 = m² -------- (a1)
の正整数解(m,n)を求める。

つまり、楕円曲線
     C₁: n(n+1)(n+2)/6 = m²
の整点(m,n)を求める。
ここで、
     X = 6n+6,
     Y = 36m

とすると、楕円曲線
     E₁: Y² = X³-36X ------- (a2)
となる。
m,nが正整数なら、X,Yも正整数になるので、楕円曲線E₁の整点を求めれば十分である。
Y^2=X^3-36Xの有理点と整点で既に示したように、楕円曲線E₁の整点は、
     (±6, 0), (-3, ±9), (-2, ±8), (0, 0), (12,±36), (18, ±72), (294,±5040)
の13個に限る。
m,nが整数とすると、X=6n+6は6の倍数、Y=36mは36の倍数なので、楕円曲線C₁の整点(m,n)は、
     (0,-2), (0,-1), (0, 0), (±1, 1), (±2, 2), (±140, 48)
の9個に限る。
よって、(a1)の解で、m,nが正整数になるのは、
     (1, 1), (2, 2), (140, 48)
の3個に限る。
つまり、問題Ａを満たす平方数m²は、
     1, 4, 19600
の3個に限る。

■問題Ｂ: 三角数m(m+1)/2で、1からn(n+1)/2までの連続する三角数の和に等しくなるものを全て求めよ。

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

     X = 12(n+1),
     Y = 36(2m+1)

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

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

gp> read("c12.gp");
time = 14 ms.
gp> e=ellinit([0,0,0,-144,1296])
%2 = [0, 0, 0, -144, 1296, 0, -288, 5184, -20736, 6912, -1119744, -534491136, -110592/179, Vecsmall([1]), [Vecsmall([128, -1])], [0, 0, 0, 0, 0, 0, 0, 0]]
gp> e.disc
%3 = -534491136
gp> e.j
%4 = -110592/179
gp> ellglobalred(e)
%5 = [1611, [2, 0, 0, 4], 4, [3, 2; 179, 1], [[2, -1, 0, 4], [1, 5, 0, 1]]]
gp> factor(e.disc)
%6 =
[ -1  1]

[  2 12]

[  3  6]

[179  1]

gp> factor(e.j)
%7 =
[ -1  1]

[  2 12]

[  3  3]

[179 -1]

gp> factor(1611)
%8 =
[  3 2]

[179 1]

■楕円曲線Eは、自明なねじれ点群{O}を持つ。
[pari/gpによる計算]

gp> elltors(e)
%9 = [1, [], []]

■楕円曲線EのModell-Weil群のrankは、2であることが分かる。
Cremonaのmwrank3によって、このようにMordell-Weil群の生成元P₁(0,36),P₂(48,324)が見つかる。よって、
     E(Q) = Z×Z
である。

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

f(x)=x³-144x+1296=0の根は、１個の実数と２個の虚数であり、
     γ'' ≒ 7.575306764442961497746139877 + 5.306205586503963424151249099*sqrt(-1),
     γ' ≒ 7.575306764442961497746139877 - 5.306205586503963424151249099*sqrt(-1),
     γ ≒ -15.15061352888592299549227975
である。
任意のP ∈ E(Q)に対して、有理整数m₁,m₂が存在して、
     P = m₁P₁+m₂P₂ ------ (5)
となる。

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

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

gp> f(x)=x^3-144*x+1296
%10 = (x)->x^3-144*x+1296
gp> poldisc(f(x))
time = -3 ms.
%11 = -33405696
gp> rr=polroots(f(x))
time = 1 ms.
%12 = [-15.150613528885922995492279755584786450 + 0.E-38*I, 7.5753067644429614977461398777923932248 - 5.3062055865039634241512490993181956734*I, 7.5753067644429614977461398777923932248 + 5.3062055865039634241512490993181956734*I]~
gp> factor(f(x))
%13 =
[x^3 - 144*x + 1296 1]

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

[pari/gpによる計算]

gp> p1=[0,36];p2=[48,324];t1=[0];
gp> H=HH2(e,p1,p2)
time = 1 ms.
%15 =
[0.17026182626464410560603174839512503263 -0.30537866986083332725413814444054485456]

[-0.30537866986083332725413814444054485456 0.74317635666607524617557902250693275584]

gp> dd=matdet(H-x*matid(2))
%16 = x^2 - 0.91343818293071935178161077090205778848*x + 0.033278431716698652728242055268180095359
gp> ww=polroots(dd)
%17 = [0.038014067860143153991285402741076699421 + 0.E-38*I, 0.87542411507057619779032536816098108905 + 0.E-38*I]~
gp> c1=min(real(ww[1]),real(ww[2]))
%18 = 0.038014067860143153991285402741076699421

■不等式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(rr[1]),abs(rr[2]),abs(rr[3]))
%19 = 30.301227057771845990984559511169572899

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

■不等式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)
となる。ただし、r=rank(E)=2である。

[pai/gpによる計算]

gp> X0=floor(max3(c2,real(rr[1]),0))+1
%22 = 31
gp> factor(e.j)
%23 =
[ -1  1]

[  2 12]

[  3  3]

[179 -1]

gp> he=log(3^4)
%24 = 4.3944491546724387655809809476901028186
gp> om=omega_pe(-144,1296)
%25 = [1.4748693591176675578857991565282814667 + 0.65246303538110892367330286105363575772*I, 1.4748693591176675578857991565282814667 - 0.65246303538110892367330286105363575772*I]
gp> om1=om[1];om2=om[2];
gp> A0=AA(e,t1,he,om)
%27 = 8.4912576069323132872335977970055048518
gp> A1=AA(e,p1,he,om)
time = 13 ms.
%28 = 4.3944491546724387655809809476901028186
gp> A2=AA(e,p2,he,om)
time = 1 ms.
%29 = 4.3944491546724387655809809476901028186
gp> eb=eeb2(e,p1,p2,A0,A1,A2,om)
time = 1 ms.
%30 = 2.2194678189347779213019302896307597563
gp> c4=cc4(2,eb,A0*A1*A2)
%31 = 25652316465063.642095534713171808532845
gp> c5=cc5(eb)
%32 = 0.79726744594591780901099344226782543171
gp> c6=cc6(c5,he)
%33 = 5.1917166006183565745919743899579282503
gp> c7=cc7(c6,2,1)
%34 = 5.9161137811783018840092065114161048183
gp> c8=cc8(c6,2,1)
%35 = 5.4529876556253016012119738015282555326
gp> g(m)=m^2-(c3/c1+log(4*sqrt(2))/c1+(c4/c1)*(log(m)+c7)*(log(log(m))+c8)^4)
%36 = (m)->m^2-(c3/c1+log(4*sqrt(2))/c1+(c4/c1)*(log(m)+c7)*(log(log(m))+c8)^4)
gp> gdash(m)=2*m-(c4/c1)*((1/m)*(log(log(m))+c8)^4+(log(m)+c7)*4*(log(log(m))+c8)^3*(1/log(m))*(1/m))
%39 = (m)->2*m-(c4/c1)*((1/m)*(log(log(m))+c8)^4+(log(m)+c7)*4*(log(log(m))+c8)^3*(1/log(m))*(1/m))
gp> nn(x)=x-g(x)/gdash(x)
%40 = (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
......
79058626742.673816081394418885450737863
40260402534.714685123123850606332786532
21517983034.360650654370761529925322976
13281316387.269607356943483559372021396
10662090146.025100334830078663056753389
10325689401.725526464004858126549959428
10319910356.213226705723798705450326608
10319908647.816001153304818782618345206
10319908647.815851852257081785144728478
10319908647.815851852257081784004450193
10319908647.815851852257081784004450193
10319908647.815851852257081784004450193
.....
10319908647.815851852257081784004450193
10319908647.815851852257081784004450193
10319908647.815851852257081784004450193
10319908647.815851852257081784004450193
time = 24 ms.

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

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

[pari/gpによる計算]

gp> K3=10319908647.815851852257081784004450193
%43 = 10319908647.815851852257081784004450193
gp> K1=(4*sqrt(2)/2*real(om1))*exp(c3)
%44 = 156.82293314442302002308727617214338732
gp> K2=c1
%45 = 0.038014067860143153991285402741076699421
gp> KK0(2,1,K3)
%46 = 129224369336101355531324680441856.67159

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

[pari/gpによる計算]

gp> default(realprecision,100)
gp> K0=10^50
%48 = 100000000000000000000000000000000000000000000000000
gp> a1=floor(K0*phi(e,p1,K0)+0.5)
time = 2 ms.
%49 = 69514738651492599912791387887984369050742880155213
gp> a2=floor(K0*phi(e,p2,K0)+0.5)
time = 1 ms.
%50 = 19681232180225008935930294004268900375698080946496
gp> aaa=[1,0,0;0,1,0;a1,a2,K0]
%51 =
[1 0 0]

[0 1 0]

[69514738651492599912791387887984369050742880155213 19681232180225008935930294004268900375698080946496 100000000000000000000000000000000000000000000000000]

gp> bbb=qflll(aaa,1)
%52 =
[-17025378326681629  4330230206849322  59738315282011456]

[  8914936734530043 23321369383585760 -78956458516445463]

[ 10080577850777571 -7600081067290717 -25987329821190254]

gp> b1n=nr(bbb)
%53 = 21701558790854197.17124316440699721681287234023289972418319794981066110912449114633822091537073183113
gp> b1n-2^(2/2)*1*K3*sqrt(2^2+2)
%54 = 21701508233833438.60255067348390579019664057163642631906614
gp> M2=(1/K2)*(log(K0*K1)-log(sqrt(b1n^2/2^2-2*K3^2)-2*K3))
%55 = 2190.27758109719260548425708365447111722
gp> sqrt(M2)
%56 = 46.8004015057263435598032858018276889119

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

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

gp> default(realprecision,35)
gp> K3=46
%58 = 46
gp> KK0(2,1,K3)
%59 = 11444329.612970084321574280449387487
gp> K0=10^9
%60 = 1000000000
gp> a1=floor(K0*phi(e,p1,K0)+0.5)
%61 = 695147386
gp> a2=floor(K0*phi(e,p2,K0)+0.5)
%62 = 196812322
gp> aaa=[1,0,0;0,1,0;a1,a2,K0]
%63 =
[        1         0          0]

[        0         1          0]

[695147386 196812322 1000000000]

gp> bbb=qflll(aaa,1)
%64 =
[-266  101 -1037]

[ -97  756  1610]

[ 204 -219   404]

gp> b1n=nr(bbb)
%65 = 348.97134552854049480859946444525530
gp> b1n-2^(2/2)*1*K3*sqrt(2^2+2)
%66 = 123.61828919248810977444932957231329
gp> M2=(1/K2)*(log(K0*K1)-log(sqrt(b1n^2/2^2-2*K3^2)-2*K3))
%67 = 566.40202691853341102288347404745862
gp> sqrt(M2)
%68 = 23.799202232817246494630643680948922

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

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

gp> K3=23
%69 = 23
gp> KK0(2,1,K3)
%70 = 1430541.2016212605401967850561734359
gp> K0=10^8
%71 = 100000000
gp> a1=floor(K0*phi(e,p1,K0)+0.5)
%72 = 69514739
gp> a2=floor(K0*phi(e,p2,K0)+0.5)
%73 = 19681232
gp> aaa=[1,0,0;0,1,0;a1,a2,K0]
%74 =
[       1        0         0]

[       0        1         0]

[69514739 19681232 100000000]

gp> bbb=qflll(aaa,1)
%75 =
[ 266 -256 -165]

[  97   76  659]

[-204  163  -15]

gp> b1n=nr(bbb)
%76 = 348.97134552854049480859946444525530
gp> b1n-2^(2/2)*1*K3*sqrt(2^2+2)
%77 = 236.29481736051430229152439700878430
gp> M2=(1/K2)*(log(K0*K1)-log(sqrt(b1n^2/2^2-2*K3^2)-2*K3))
%78 = 490.45194184557229972402254323326853
gp> sqrt(M2)
%79 = 22.146149594129727627450923692517537

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

■|m₁|,|m₂| <= 22について、P=m₁P₁+m₂P₂が整点かどうかを確認する。

[pari/gpによる計算]

gp> check2(22,e,p1,p2)
[420, -8604]
[108, 1116]
[-8, 44]
[-15, -9]
[9, -27]
[24, -108]
[252, -3996]
[4953, 348579]
[40, 244]
[12, 36]
[4, 28]
[-12, -36]
[0, -36]
[48, -324]
[48, 324]
[0, 36]
[-12, 36]
[4, -28]
[12, -36]
[40, -244]
[4953, -348579]
[252, 3996]
[24, 108]
[9, 27]
[-15, 9]
[-8, -44]
[108, -1116]
[420, 8604]
There are 28 integral points.
time = 156 ms.

よって、楕円曲線Eの整点(X,Y)は、
     (-15,±9), (-12,±36), (-8,±44), (0,±36), (4,±28), (9,±27), (12,±36),
     (24,±108), (40,±244), (48,±324), (108,±1116), (252,±3996), (420,±8604), (4953,±348579)
の28個に限る。

■最後に、楕円曲線C₂の整点(m,n)を求める。
m,nを整数とすると、X=12(n+1)は12の倍数,Y=36(2m+1)は36の倍数かつ(Y/36)は奇数でなければならないので、
     (X,Y) = (-12,±36), (0,±36), (12,±36), (24,±108), (252, ±3996), (420,±8604)
のみが条件を満たす。よって、楕円曲線C₂の整点(m,n)は、
     (-1,-2), (0,-2), (-1,-1), (0,-1), (-1,0), (0,0),
     (-2,1), (1,1), (-56,20), (55,20), (-120,34), (119,34)
の12個に限る。

よって、(1)の解で、m,nが正整数になるものは、
     (1, 1), (55,20), (119, 34)
の3個に限る。
つまり、問題Ｂを満たす三角数m(m+1)/2は、
     1, 1540, 7140
の3個に限る。

[2004.08.22追記]
■SIMATH-4.6(simcalc)によって、楕円曲線Eの整点を求める。
simcalcによって、楕円曲線EのMordell-Weil群E(Q)のrankは2であり、その生成元は(0,36), (12,36)である。
また、楕円曲線Eの整点を(Pと-Pを同一視して)求めると、
  &nbso; (-15,9), (-12,36), (-8,44), (0,36), (4,28), (9,27), (12,36), (24,108), (40,244), (48,324),
  &nbso; (108,1116), (252,3996), (420,8604), (4953,348579)
の14個となる。これらの点のそれぞれを(-1)倍した点も整点となるので、楕円曲線E葉合計2*14=28個の整点を持つ。

[simcalcによる計算]

> E=EC(0,0,0,-144,1296)
                       simcalc in free(): warning: junk pointer, too high to make sense.
         E = EC(-144, 1296)
> basismwg(E)
  basis :  PT(0, 36, 1)
           PT(12, 36, 1)
simcalc in free(): warning: junk pointer, too high to make sense.
         @ = 2
> faintp(E)
  all nontrivial integral points modulo negation :
  PT(0, 36, 1) = PT(0, 1, 0) + PT(0, 36, 1) + 0*PT(12, 36, 1)
  PT(4, 28, 1) = PT(0, 1, 0) - 2*PT(0, 36, 1) + 0*PT(12, 36, 1)
  PT(252, 3996, 1) = PT(0, 1, 0) + 3*PT(0, 36, 1) + 0*PT(12, 36, 1)
  PT(-15, 9, 1) = PT(0, 1, 0) + 2*PT(0, 36, 1) - PT(12, 36, 1)
  PT(24, 108, 1) = PT(0, 1, 0) + PT(0, 36, 1) - PT(12, 36, 1)
  PT(12, 36, 1) = PT(0, 1, 0) + 0*PT(0, 36, 1) + PT(12, 36, 1)
  PT(-12, 36, 1) = PT(0, 1, 0) - PT(0, 36, 1) - PT(12, 36, 1)
  PT(48, 324, 1) = PT(0, 1, 0) - 2*PT(0, 36, 1) - PT(12, 36, 1)
  PT(108, 1116, 1) = PT(0, 1, 0) - PT(0, 36, 1) + 2*PT(12, 36, 1)
  PT(-8, 44, 1) = PT(0, 1, 0) + 0*PT(0, 36, 1) + 2*PT(12, 36, 1)
  PT(9, 27, 1) = PT(0, 1, 0) - PT(0, 36, 1) - 2*PT(12, 36, 1)
  PT(40, 244, 1) = PT(0, 1, 0) + 2*PT(0, 36, 1) + 2*PT(12, 36, 1)
  PT(420, 8604, 1) = PT(0, 1, 0) - PT(0, 36, 1) - 3*PT(12, 36, 1)
  PT(4953, 348579, 1) = PT(0, 1, 0) + 4*PT(0, 36, 1) + 3*PT(12, 36, 1)
 
simcalc in free(): warning: junk pointer, too high to make sense.
         @ = PT(0, 36, 1)

[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.

Last Update: 2021.03.18

H.Nakao