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

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

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

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

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

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

     X = 6n+3,
     Y = 18m+9

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

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

gp> read("c12.gp");
time = 1 ms.
gp> e=ellinit([0,0,0,-9,81])
%2 = [0, 0, 0, -9, 81, 0, -18, 324, -81, 432, -69984, -2787696, -6912/239, Vecsmall([1]), [Vecsmall([128, -1])], [0, 0, 0, 0, 0, 0, 0, 0]]
gp> e.disc
%3 = -2787696
gp> e.j
%4 = -6912/239
gp> factor(e.disc)
%5 =
[ -1 1]

[  2 4]

[  3 6]

[239 1]

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

[  2  8]

[  3  3]

[239 -1]

gp> ellglobalred(e)
%7 = [8604, [1, 0, 0, 0], 12, [2, 2; 3, 2; 239, 1], [[2, 4, 0, 3], [2, -1, 0, 4], [1, 5, 0, 1]]]
gp> factor(8604)
%8 =
[  2 2]

[  3 2]

[239 1]

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

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

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

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

f(x)=x³-9x+81=0の根は、１個の実数と２個の虚数であり、
     γ'' ≒ 2.507549822485741454622224671 + 3.140607956549944875841943233*sqrt(-1),
     γ' ≒ 2.507549822485741454622224671 - 3.140607956549944875841943233*sqrt(-1),
     γ ≒ -5.015099644971482909244449343
である。
任意の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-9*x+81
%10 = (x)->x^3-9*x+81
gp> poldisc(f(x))
%11 = -174231
gp> factor(f(x))
%12 =
[x^3 - 9*x + 81 1]

gp> rr=polroots(f(x))
time = 1 ms.
%13 = [-5.0150996449714829092444493436586716862 + 0.E-38*I, 2.5075498224857414546222246718293358431 - 3.1406079565499448758419432333452146006*I, 2.5075498224857414546222246718293358431 + 3.1406079565499448758419432333452146006*I]~

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

[pari/gpによる計算]

gp> p1=[-3,9];p2=[-15/4,63/8];t1=[0];
gp> H=HH2(e,p1,p2)
time = 1 ms.
%15 =
[0.23038586279911795215271127028625550689 -0.47148796867422587606789275355968469678]

[-0.47148796867422587606789275355968469678 1.1688563250343843247769375129758392515]

gp> dd=matdet(H-x*matid(2))
%16 = x^2 - 1.3992421878335022769296487832620947584*x + 0.046987068326705084110122071226770769744
gp> ww=polroots(dd)
%17 = [0.034427432847610738366333563126332942180 + 0.E-38*I, 1.3648147549858915385633152201357618162 + 0.E-38*I]~
gp> c1=min(real(ww[1]),real(ww[2]))
time = -6 ms.
%18 = 0.034427432847610738366333563126332942180

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

■不等式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)
%24 = 0.E-38
gp> c3=cc3(c0,e)
time = 1 ms.
%25 = 2.9336799873410002326722821098791591306

■不等式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[1]),0))+1
%26 = 11
gp> factor(e.j)
%27 =
[ -1  1]

[  2  8]

[  3  3]

[239 -1]

gp> he=log(3^4)
%28 = 4.3944491546724387655809809476901028186
gp> om=omega_pe(-9,81)
%29 = [2.1248747502003803728849396431561236164 + 1.1111767176282767763568337844755320897*I, 2.1248747502003803728849396431561236164 - 1.1111767176282767763568337844755320897*I]
gp> om1=om[1];om2=om[2];
gp> A0=AA(e,t1,he,om)
%31 = 7.6504437409494757532862117077458296334
gp> A1=AA(e,p1,he,om)
time = 2 ms.
%32 = 5.1673603264978260935791665462900104523
gp> A2=AA(e,p2,he,om)
time = 20 ms.
%33 = 5.7031454245393464396320591018308838002
gp> eb=eeb2(e,p1,p2,A0,A1,A2,om)
time = 1 ms.
%34 = 2.2194678189347779213019302896307597563
gp> c4=cc4(2,eb,A0*A1*A2)
%35 = 35270815627116.902386853901347012301379
gp> c5=cc5(eb)
%36 = 0.79726744594591780901099344226782543171
gp> c6=cc6(c5,he)
%37 = 5.1917166006183565745919743899579282503
gp> c7=cc7(c6,2,1)
%38 = 5.9161137811783018840092065114161048183
gp> c8=cc8(c6,2,1)
%39 = 5.4529876556253016012119738015282555326
gp> g(m)=m^2-(c3/c1+log(4*sqrt(2))/c1+(c4/c1)*(log(m)+c7)*(log(log(m))+c8)^4)
%40 = (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))
%41 = (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)
%42 = (x)->x-g(x)/gdash(x)
(14:42) 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
........
157494321735.68450817136879773325574003
79324056017.950095821305276353445965056
40768713938.607838450139185793440062565
22468237835.727304367766749960740235348
14919426438.597522463583712560910698475
12948645773.861841697701385325473020779
12791265638.902647584851748385297708358
12790244312.880317667387492803127039021
12790244269.846286432859826947477742300
12790244269.846286356457164531706556910
12790244269.846286356457164531706556669
12790244269.846286356457164531706556669
12790244269.846286356457164531706556669
12790244269.846286356457164531706556669
12790244269.846286356457164531706556669
12790244269.846286356457164531706556669
.........
12790244269.846286356457164531706556669
12790244269.846286356457164531706556669
12790244269.846286356457164531706556669
time = 29 ms.

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

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

[pari/gpによる計算]

gp> K3=12790244269.846286356457164531706556669
%45 = 12790244269.846286356457164531706556669
gp> K1=(4*sqrt(2)/2*real(om1))*exp(c3)
%46 = 112.96901953753357455666098751398240262
gp> K2=c1
%47 = 0.034427432847610738366333563126332942180
gp> KK0(2,1,K3)
%48 = 246010350013335753701143962935344.02989

K₀=10⁵⁰とする。 A=[1 0 0 ; 0 1 0 ; [K₀*φ(P₁)] [K₀*φ(P₂)] K₀] にLLL-algorithmを適用して、reduced basis {b₁,b₂,b₃}を求めると、
     b₁ = [ 10792794932040232, 15067757935118734, -21879585984404525 ]^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
%50 = 100000000000000000000000000000000000000000000000000
gp> a1=floor(K0*phi(e,p1,K0)+0.5)
time = 1 ms.
%51 = 82184716303503834794016770144123154758030557458348
gp> a2=floor(K0*phi(e,p2,K0)+0.5)
time = -7 ms.
%52 = 86340370905260401297200875459673620931474624375992
gp> aaa=[1,0,0;0,1,0;a1,a2,K0]
%53 =
[1 0 0]

[0 1 0]

[82184716303503834794016770144123154758030557458348 86340370905260401297200875459673620931474624375992 100000000000000000000000000000000000000000000000000]

gp> bbb=qflll(aaa,1)
time = 1 ms.
%54 =
[ 10792794932040232 -7246502287178037 -137249143204007733]

[ 15067757935118734 16055651220972700   72398916443544413]

[-21879585984404525 -7906991468798613   50288325982457660]

gp> b1n=nr(bbb)
%55 = 28674693276221069.59680740602258382132409643548931726340696146147229509822633161874758473642291353283
gp> b1n-2^(2/2)*1*K3*sqrt(2^2+2)
%56 = 28674630617076776.23721234667597424157635541810607505925326
gp> M2=(1/K2)*(log(K0*K1)-log(sqrt(b1n^2/2^2-2*K3^2)-2*K3))
%57 = 2400.83907030590258649360967366164851129
gp> sqrt(M2)
%58 = 48.9983578327468295827381262351281865664

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

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

gp> default(realprecision,35)
gp> K3=49
%65 = 49
gp> KK0(2,1,K3)
%66 = 13832640.899937509763590989157043524
gp> K0=10^9
%67 = 1000000000
gp> a1=floor(K0*phi(e,p1,K0)+0.5)
time = 1 ms.
%68 = 821847163
gp> a2=floor(K0*phi(e,p2,K0)+0.5)
time = 2 ms.
%69 = 863403709
gp> aaa=[1,0,0;0,1,0;a1,a2,K0]
%70 =
[        1         0          0]

[        0         1          0]

[821847163 863403709 1000000000]

gp> bbb=qflll(aaa,1)
%71 =
[ -53  191  1040]

[ 728 -703   364]

[-585  450 -1169]

gp> b1n=nr(bbb)
%72 = 935.42396804871319268870125041422579
gp> b1n-2^(2/2)*1*K3*sqrt(2^2+2)
%73 = 695.37397325596173906536741109304843
gp> M2=(1/K2)*(log(K0*K1)-log(sqrt(b1n^2/2^2-2*K3^2)-2*K3))
%74 = 567.91093987136755448020391562584246
gp> sqrt(M2)
%75 = 23.830882062386351260271170241535840

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

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

[pari/gpによる計算]

gp> check2(23,e,p1,p2)
[5112, 365499]
[1099, -36433]
[33, 189]
[-5, 1]
[39, -243]
[24, -117]
[0, 9]
[3, -9]
[7, 19]
[9, 27]
[-3, -9]
[513, 11619]
[513, -11619]
[-3, 9]
[9, -27]
[7, -19]
[3, 9]
[0, -9]
[24, 117]
[39, 243]
[-5, -1]
[33, -189]
[1099, 36433]
[5112, -365499]
There are 24 integral points.
time = 298 ms.

よって、楕円曲線Eの整点(X,Y)は、
     (-5,±1), (-3,±9), (0,±9), (3,±9), (7,±19), (9,±27),
     (24,±117), (33,±189), (39,±243), (513,±11619), (1099,±36433), (5112,±365499)
の24個に限る。

■最後に、楕円曲線Cの整点(m,n)を求める。
m,nを整数とすると、X=6n+3は3の倍数,Y=18m+9は9の倍数であり、どちらも奇数でなければならないので、
     (X,Y) = (3, ±9), (9, ±27), (33, ±189), (39, ±243), (513, ±11619)
のみが条件を満たす。よって、楕円曲線Cの整点(m,n)は、
     (0, 0), (-1, 0), (-2, 1), (1,1), (-11, 5), (10, 5), (-14, 6), (13, 6), (-646, 85), (645, 85)
の10個に限る。

よって、(1)の解で、m,nが正整数になるものは、
     (1, 1), (10, 5), (13, 6), (645, 85)
の4個に限る。
よって、三角数m(m+1)/2で、1²からn²までの平方数の和に等しいものは、
     1, 55, 91, 208335
の4個に限る。

[2004.08.22追記]
■SIMATH-4.6(simcalc)によって、楕円曲線Eの整点を求める。
simcalcによると、楕円曲線EのMordell-Weil群E(Q)のrankは2であり、その生成元は(3,9),(-3,9)である。
また、楕円曲線Eの整点を(Pと-Pを同一視して)求めると、
     (-5,1), (-3,9), (0,9), (3,9), (7,19), (9,27), (24,117), (33,189), (39,243), (513,11619),
     (1099,36433), (5112,365499)
の12個となる。
これらの点のそれぞれを(-1)倍した点も整点なので、楕円曲線Eは合計2*12=24個の整点を持つ。

[simcalcによる計算]

> E=EC(0,0,0,-9,81)
                   simcalc in free(): warning: junk pointer, too high to make sense.
         E = EC(-9, 81)
> basismwg(E)
  basis :  PT(3, 9, 1)
           PT(-3, 9, 1)
simcalc in free(): warning: junk pointer, too high to make sense.
         @ = 2
> faintp(E)
  all nontrivial integral points modulo negation :
  PT(3, 9, 1) = PT(0, 1, 0) + PT(3, 9, 1) + 0*PT(-3, 9, 1)
  PT(-5, 1, 1) = PT(0, 1, 0) - 2*PT(3, 9, 1) + 0*PT(-3, 9, 1)
  PT(1099, 36433, 1) = PT(0, 1, 0) + 4*PT(3, 9, 1) + 0*PT(-3, 9, 1)
  PT(24, 117, 1) = PT(0, 1, 0) + 2*PT(3, 9, 1) - PT(-3, 9, 1)
  PT(9, 27, 1) = PT(0, 1, 0) - PT(3, 9, 1) + PT(-3, 9, 1)
  PT(-3, 9, 1) = PT(0, 1, 0) + 0*PT(3, 9, 1) + PT(-3, 9, 1)
  PT(0, 9, 1) = PT(0, 1, 0) - PT(3, 9, 1) - PT(-3, 9, 1)
  PT(33, 189, 1) = PT(0, 1, 0) - 2*PT(3, 9, 1) - PT(-3, 9, 1)
  PT(7, 19, 1) = PT(0, 1, 0) + 0*PT(3, 9, 1) - 2*PT(-3, 9, 1)
  PT(39, 243, 1) = PT(0, 1, 0) + PT(3, 9, 1) + 2*PT(-3, 9, 1)
  PT(513, 11619, 1) = PT(0, 1, 0) + PT(3, 9, 1) - 3*PT(-3, 9, 1)
  PT(5112, 365499, 1) = PT(0, 1, 0) - 3*PT(3, 9, 1) - 3*PT(-3, 9, 1)
 
simcalc in free(): warning: junk pointer, too high to make sense.
         @ = PT(3, 9, 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