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>  e=ellinit([0,0,0,-9,81])
time = 200 ms.
%32 = [0, 0, 0, -9, 81, 0, -18, 324, -81, 432, -69984, -2787696, -6912/239, [-5.015099644971482909244449343, 2.507549822485741454622224671 - 3.140607956549944875841943233*I, 2.507549822485741454622224671 + 3.140607956549944875841943233*I]~, 2.124874750200380372884939643, -1.062437375100190186442469821 + 0.5555883588141383881784168922*I, -2.657786728371786154340090028 - 4.94599652 E-29*I, 1.328893364185893077170045014 - 2.173411877405180254606148302*I, 1.180555675149431606555725242]
gp>  e.disc
time = 0 ms.
%33 = -2787696
gp>  e.j
time = 0 ms.
%34 = -6912/239
gp>  factor([e.disc])
time = 4 ms.
%35 = [[-1, 1; 2, 4; 3, 6; 239, 1]]
gp>  factor([e.j])
time = 1 ms.
%36 = [[-1, 1; 2, 8; 3, 3; 239, -1]]
gp>  ellglobalred(e)
time = 10 ms.
%37 = [8604, [1, 0, 0, 0], 12]
gp>  factor([8604])
time = 2 ms.
%38 = [[2, 2; 3, 2; 239, 1]]

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

gp>  elltors(e,1)
time = 99 ms.
%39 = [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>  poldisc(x^3-9*x+81)
time = 3 ms.
%42 = -174231
gp>  rr=polroots(x^3-9*x+81)
time = 37 ms.
%43 = [-5.015099644971482909244449343 + 0.E-28*I, 2.507549822485741454622224671 - 3.140607956549944875841943233*I, 2.507549822485741454622224671 + 3.140607956549944875841943233*I]~

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

[pari/gpによる計算]

gp>  read("c1.gp")
time = 71 ms.
gp>  p1=[-3,9];p2=[-15/4,63/8];t1=[0];
time = 0 ms.
gp>  H=HH2(e,p1,p2)
time = 119 ms.
%45 = 
[0.2303858627991179521527112701 -0.4714879686742258760678927537]

[-0.4714879686742258760678927537 1.168856325034384324776937512]

gp>  dd=matdet(H-x*matid(2))
time = 1 ms.
%46 = x^2 - 1.399242187833502276929648782*x + 0.04698706832670508411012207090
gp>  ww=polroots(dd)
time = 5 ms.
%47 = [0.03442743284761073836633356289 + 0.E-28*I, 1.364814754985891538563315220 + 0.E-28*I]~
gp>  c1=min(real(ww[1]),real(ww[2]))
time = 0 ms.
%48 = 0.03442743284761073836633356289

■不等式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]))
time = 2 ms.
%49 = 10.03019928994296581848889868

■不等式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.
%50 = 0.E-28
gp>  c3=cc3(c0,e)
time = 4 ms.
%51 = 2.933679987341000232672282109

■不等式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
time = 0 ms.
%53 = 11
gp>  factor([e.j])
time = 0 ms.
%54 = [[-1, 1; 2, 8; 3, 3; 239, -1]]
gp>  he=log(3^4)
time = 0 ms.
%55 = 4.394449154672438765580980947
gp>  A0=AA(e,t1,he)
time = 9 ms.
%56 = 24.03034836538438580719490188
gp>  A1=AA(e,p1,he)
time = 46 ms.
%57 = 16.23088450550429511709565855
gp>  A2=AA(e,p2,he)
time = 12 ms.
%58 = 17.91380682882059825700047169
gp>  eb=eeb2(e,p1,p2,A0,A1,A2)
time = 25 ms.
%59 = 1.044517009167470794935226451
gp>  c4=cc4(2,eb,A0*A1*A2)
time = 1 ms.
%60 = 752708443513102777102360.1015
gp>  c5=cc5(eb)
time = 0 ms.
%61 = 0.04355458638880544089850357123
gp>  c6=cc6(c5,he)
time = 0 ms.
%62 = 4.438003741061244206479484518
gp>  c7=cc7(c6,2,1)
time = 0 ms.
%63 = 5.162400921621189515896716640
gp>  c8=cc8(c6,2,1)
time = 0 ms.
%64 = 4.699274796068189233099483930
gp>  g(m)=m^2-(c3/c1+log(4*sqrt(2))/c1+(c4/c1)*(log(m)+c7)*(log(log(m))+c8)^4)
time = 0 ms.
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))
time = 0 ms.
gp>  nn(x)=x-g(x)/gdash(x)
time = 0 ms.
gp> x=10^100;for(i=1,300,x=nn(x);print(x))
5.000000000000000000000000000 E99
2.500000000000000000000000000 E99
1.250000000000000000000000000 E99
6.250000000000000000000000000 E98
3.125000000000000000000000000 E98
......
10438628997671819.24127327466
5425290989442511.425174653576
3100079031022099.857772674053
2217418671155232.180866448482
2037050898896241.139909234971
2028773086561823.431563607510
2028755541812728.211489983689
2028755541733904.307242603909
2028755541733904.307241012875
2028755541733904.307241012875
2028755541733904.307241012875
2028755541733904.307241012875
2028755541733904.307241012875
2028755541733904.307241012875
2028755541733904.307241012875
2028755541733904.307241012875
2028755541733904.307241012875
2028755541733904.307241012875
2028755541733904.307241012875
2028755541733904.307241012875
2028755541733904.307241012875
2028755541733904.307241012875
time = 1,794 ms.

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

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

[pari/gpによる計算]

gp>  K3=2028755541733904.307241012875
time = 0 ms.
%88 = 2028755541733904.307241012874
gp>  K1=(4*sqrt(2)/2*real(e.omega[1]))*exp(c3)
time = 0 ms.
%89 = 112.9690195375335745566609874
gp>  K2=c1
time = 0 ms.
%90 = 0.03442743284761073836633356289
gp>  KK0(2,1,K3)
time = 0 ms.
%91 = 9.817615517244816998237896737 E47

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)
   realprecision = 105 significant digits (100 digits displayed)
time = 1 ms.
gp>  ec=ellinit([0,0,0,-9,81])
time = 44 ms.
gp>  ec=ellinit([0,0,0,-9,81])
time = 25 ms.
%92 = [0, 0, 0, -9, 81, 0, -18, 324, -81, 432, -69984, -2787696, -6912/239, [-5.015099644971482909244449343658671686184794481613745076410707315942575390065379758620913582655434452, 2.507549822485741454622224671829335843092397240806872538205353657971287695032689879310456791327717226 - 3.140607956549944875841943233345214600579718592213725291743167276300071492813005074986055244174766996*I, 2.507549822485741454622224671829335843092397240806872538205353657971287695032689879310456791327717226 + 3.140607956549944875841943233345214600579718592213725291743167276300071492813005074986055244174766996*I]~, 2.124874750200380372884939643156123616440540356483296417341217233069119137307271379765721917510006923, -1.062437375100190186442469821578061808220270178241648208670608616534559568653635689882860958755003461 + 0.5555883588141383881784168922377660448522397244853580209118437459567988879100580546221255814316234064*I, -2.657786728371786154340090028377598051271174069473138193813456983867443612567487070680428973381680650 - 2.13572298 E-106*I, 1.328893364185893077170045014188799025635587034736569096906728491933721806283743535340214486690840325 - 2.173411877405180254606148302638903408320782385574084597869511969282979500335111100968639369590639633*I, 1.180555675149431606555725242362851636769110659681188864284553420424952894954491270544788607797738453]
gp>  K0=10^50
time = 0 ms.
%93 = 100000000000000000000000000000000000000000000000000
gp>  a1=floor(K0*phi(ec,p1,K0)+0.5)
time = 18 ms.
%94 = 82184716303503834794016770144123154758030557458348
gp>  a2=floor(K0*phi(ec,p2,K0)+0.5)
time = 23 ms.
%95 = 86340370905260401297200875459673620931474624375992
gp>  aaa=[1,0,0;0,1,0;a1,a2,K0]
time = 0 ms.
%96 = 
[1 0 0]

[0 1 0]

[82184716303503834794016770144123154758030557458348 86340370905260401297200875459673620931474624375992 100000000000000000000000000000000000000000000000000]

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

[15067757935118734 16055651220972700 72398916443544413]

[-21879585984404525 -7906991468798613 50288325982457660]

gp>  b1n=nr(bbb)
time = 1 ms.
%98 = 28674693276221069.59680740602258382132409643548931726340696146147229509822633161874758473642291353283
gp>  b1n-2^(2/2)*1*K3*sqrt(2^2+2)
time = 0 ms.
%99 = 18735861496037612.794649731390257587191
gp>  M2=(1/K2)*(log(K0*K1)-log(sqrt(b1n^2/2^2-2*K3^2)-2*K3))
time = 1 ms.
%100 = 2411.333560557942531683322308
gp>  sqrt(M2)
time = 0 ms.
%101 = 49.10533128447401617919912225

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

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

gp>  K3=49
time = 0 ms.
%105 = 49
gp>  KK0(2,1,K3)
time = 0 ms.
%106 = 13832640.89993750976359098915704352398592328424728527602911893669653639539747883146675462027737573531
gp>  K0=10^9
time = 0 ms.
%107 = 1000000000
gp>  a1=floor(K0*phi(ec,p1,K0)+0.5)
time = 6 ms.
%108 = 821847163
gp>  a2=floor(K0*phi(ec,p2,K0)+0.5)
time = 6 ms.
%109 = 863403709
gp>  aaa=[1,0,0;0,1,0;a1,a2,K0]
time = 1 ms.
%110 = 
[1 0 0]

[0 1 0]

[821847163 863403709 1000000000]

gp>  bbb=qflll(aaa,1)
time = 1 ms.
%111 = 
[-53 191 1040]

[728 -703 364]

[-585 450 -1169]

gp>  b1n=nr(bbb)
time = 1 ms.
%112 = 935.4239680487131926887012504142257893490035455570359795227231434571669689936694446935469573335018239
gp>  b1n-2^(2/2)*1*K3*sqrt(2^2+2)
time = 0 ms.
%113 = 695.3739732559617390653674110930484329363406924526823069363192718665728520028525720926407328892470809
gp>  M2=(1/K2)*(log(K0*K1)-log(sqrt(b1n^2/2^2-2*K3^2)-2*K3))
time = 1 ms.
%114 = 567.9109398713675544802039195
gp>  sqrt(M2)
time = 0 ms.
%115 = 23.83088206238635126027117032

よって、
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 = 7,952 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)

[参考文献]

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

H.Nakao