Temporary exact arithmetic sweep for open Diophantine equation

.github/workflows/diophantine-research.yml

@@ -0,0 +1,48 @@
+name: Exact fixed-z Diophantine sweep
+
+on:
+  push:
+    branches:
+      - research/diophantine-sage-20260727
+    paths:
+      - 'research/fixed_z_search.sage'
+      - '.github/workflows/diophantine-research.yml'
+  pull_request:
+    paths:
+      - 'research/fixed_z_search.sage'
+      - '.github/workflows/diophantine-research.yml'
+  workflow_dispatch:
+
+concurrency:
+  group: exact-fixed-z-${{ github.event.pull_request.number || github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  fixed-z:
+    runs-on: ubuntu-latest
+    timeout-minutes: 300
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          - {start: 1, end: 50}
+          - {start: 51, end: 100}
+          - {start: 101, end: 150}
+          - {start: 151, end: 200}
+    steps:
+      - uses: actions/checkout@v4
+      - name: Pull SageMath
+        run: docker pull sagemath/sagemath:latest
+      - name: Run exact fixed-z sweep
+        run: |
+          set -o pipefail
+          docker run --rm -v "$PWD:/work" -w /work sagemath/sagemath:latest \
+            sage research/fixed_z_search.sage "${{ matrix.start }}" "${{ matrix.end }}" \
+            | tee "research/results-fixed-z-${{ matrix.start }}-${{ matrix.end }}.jsonl"
+      - name: Upload exact results
+        if: always()
+        uses: actions/upload-artifact@v4
+        with:
+          name: fixed-z-${{ matrix.start }}-${{ matrix.end }}
+          path: research/results-fixed-z-${{ matrix.start }}-${{ matrix.end }}.jsonl
+          if-no-files-found: warn

research/diophantine_search.sage

@@ -0,0 +1,116 @@
+#!/usr/bin/env sage
+from sage.all import *
+import itertools, json, sys, time
+proof.all(True)
+mode = sys.argv[1] if len(sys.argv) > 1 else 'all'
+
+def emit(kind, **kw):
+    print(json.dumps({'kind': kind, **kw}, sort_keys=True), flush=True)
+
+def verify(x,y,z):
+    x,y,z=ZZ(x),ZZ(y),ZZ(z)
+    return z*z+y*y*z+x*x*x-2 == 0
+
+def map_e1(P):
+    if P.is_zero(): return None
+    X,Y=QQ(P[0]),QQ(P[1]); den=X*(X+2)
+    if den==0: return None
+    y=(2*Y+3*X+6)/den
+    if y.denominator()!=1: return None
+    y=ZZ(y); x=1-y
+    w=2+3*y-(X+1)*y*y
+    if QQ(w).denominator()!=1: return None
+    w=ZZ(w)
+    if (w-y*y)%2: return None
+    z=(w-y*y)//2
+    assert verify(x,y,z)
+    return x,y,z
+
+def map_e2(P):
+    if P.is_zero(): return None
+    U,Y=QQ(P[0]),QQ(P[1]); den=(U-232)*(U-528)
+    if den==0: return None
+    u=4*(37*Y+71*U-235068)/den
+    if u.denominator()!=1: return None
+    u=ZZ(u); y=u+8; x=3*u-7
+    r=(380-U)/148
+    w=74+(QQ(71)/37)*u+r*u*u
+    if w.denominator()!=1: return None
+    w=ZZ(w)
+    if (w-y*y)%2: return None
+    z=(w-y*y)//2
+    assert verify(x,y,z)
+    return x,y,z
+
+def lattice(E, gens, bound, mapper, label):
+    rg=range(-bound,bound+1); r=len(gens); start=time.time(); hits={}; tested=0
+    mult=[[n*P for n in rg] for P in gens]
+    cut=r//2
+    left=[]
+    for cs in itertools.product(*([rg]*cut)):
+        R=E(0)
+        for i,c in enumerate(cs): R += mult[i][c+bound]
+        left.append((cs,R))
+    right=[]
+    for cs in itertools.product(*([rg]*(r-cut))):
+        R=E(0)
+        for j,c in enumerate(cs,start=cut): R += mult[j][c+bound]
+        right.append((cs,R))
+    emit('lattice_start',label=label,rank=r,bound=bound,left=len(left),right=len(right),gens=[str(P) for P in gens])
+    for lc,LP in left:
+        for rc,RP in right:
+            cs=lc+rc
+            if all(c==0 for c in cs): continue
+            tested+=1; out=mapper(LP+RP)
+            if out:
+                key=tuple(map(int,out))
+                if key not in hits:
+                    hits[key]=cs
+                    emit('surface_hit',label=label,coefficients=list(cs),x=str(out[0]),y=str(out[1]),z=str(out[2]),x_digits=len(str(abs(out[0]))))
+    emit('lattice_done',label=label,tested=tested,hits=len(hits),seconds=time.time()-start)
+
+def conics():
+    E=EllipticCurve([0,0,0,0,-2]); P=E(3,5)
+    emit('mordell',rank=int(E.rank()),gens=[str(Q) for Q in E.gens()])
+    for n in range(1,31):
+        Q=n*P; k,v=QQ(Q[0]),QQ(Q[1]); w=abs(v); C=-k**6+40*k**3+32
+        Cn=Conic(QQ,[3*k*k*w*w,-4*w**3,-C]); t=time.time()
+        try:
+            ok,ans=Cn.has_rational_point(point=True)
+            emit('conic',multiple=n,soluble=bool(ok),witness=str(ans),seconds=time.time()-t)
+        except Exception as ex:
+            emit('conic_error',multiple=n,error=repr(ex),seconds=time.time()-t)
+
+def e1():
+    E=EllipticCurve([0,0,1,-1,6]); t=time.time(); gens=E.gens()
+    emit('curve',label='E1',rank=int(E.rank()),gens=[str(P) for P in gens],torsion=[str(P) for P in E.torsion_points()],seconds=time.time()-t)
+    try:
+        pts=E.integral_points(both_signs=True)
+        emit('integral_points',label='E1',count=len(pts),points=[str(P) for P in pts])
+        for P in pts:
+            out=map_e1(P)
+            if out: emit('integral_surface',label='E1',point=str(P),x=str(out[0]),y=str(out[1]),z=str(out[2]))
+    except Exception as ex: emit('integral_error',label='E1',error=repr(ex))
+    lattice(E,gens,100,map_e1,'E1')
+
+def e2():
+    E=EllipticCurve([0,0,0,-476688,133008912]); t=time.time()
+    known=[E(132,8508),E(QQ(23796)/121,QQ(9112140)/1331),E(QQ(18369)/49,QQ(-906975)/343),E(276,-4740),E(QQ(311724)/1369,QQ(-305064180)/50653)]
+    try:
+        gens=E.gens(); rank=E.rank()
+        emit('curve',label='E2',rank=int(rank),gens=[str(P) for P in gens],torsion=[str(P) for P in E.torsion_points()],seconds=time.time()-t)
+    except Exception as ex:
+        emit('curve_error',label='E2',error=repr(ex),seconds=time.time()-t); gens=known
+    try:
+        pts=E.integral_points(both_signs=True)
+        emit('integral_points',label='E2',count=len(pts),points=[str(P) for P in pts])
+        for P in pts:
+            out=map_e2(P)
+            if out: emit('integral_surface',label='E2',point=str(P),x=str(out[0]),y=str(out[1]),z=str(out[2]))
+    except Exception as ex: emit('integral_error',label='E2',error=repr(ex))
+    lattice(E,known,10,map_e2,'E2-known')
+    if len(gens)<=6: lattice(E,gens,8 if len(gens)>=5 else 18,map_e2,'E2-basis')
+
+if mode in ('all','conics'): conics()
+if mode in ('all','e1'): e1()
+if mode in ('all','e2'): e2()

research/eplus_search.sage

@@ -0,0 +1,60 @@
+#!/usr/bin/env sage
+from sage.all import *
+import itertools, json, sys, time
+proof.all(True)
+n = Integer(sys.argv[1]) if len(sys.argv)>1 else Integer(1)
+
+def emit(kind, **kw): print(json.dumps({'kind':kind, **kw}, sort_keys=True), flush=True)
+def verify(x,y,z): return z*z+y*y*z+x*x*x-2 == 0
+
+E0=EllipticCurve([0,0,0,0,2]); P=E0(-1,1); Q=n*P
+k,v=QQ(Q[0]),QQ(Q[1])
+s0=3*(k**4+8*k)/(v**2)
+d0=(k**4-16*k)/(v**2)
+r0=k**2/v
+nq=(k**6+40*k**3-32)/(v**3)
+E=EllipticCurve([0,0,0,s0,0])
+emit('family',multiple=int(n),k=str(k),v=str(v),s0=str(s0),d0=str(d0),curve=str(E))
+
+def to_surface(R):
+    if R.is_zero(): return None
+    T,W=QQ(R[0]),QQ(R[1])
+    x=(T*T-d0)/4; y=W/2
+    U=(T**3-3*r0*T*T+s0*T-nq)/8
+    z=-U
+    if any(q.denominator()!=1 for q in (x,y,z)): return None
+    x,y,z=ZZ(x),ZZ(y),ZZ(z)
+    assert verify(x,y,z)
+    return x,y,z
+
+t=time.time()
+try:
+    rank=E.rank(proof=True); gens=E.gens(proof=True); tors=E.torsion_points()
+    emit('curve_data',multiple=int(n),rank=int(rank),gens=[str(R) for R in gens],torsion=[str(R) for R in tors],seconds=time.time()-t)
+except Exception as ex:
+    emit('curve_error',multiple=int(n),error=repr(ex),seconds=time.time()-t); sys.exit(0)
+
+try:
+    pts=E.integral_points(both_signs=True)
+    emit('integral_points',multiple=int(n),count=len(pts),points=[str(R) for R in pts])
+    for R in pts:
+        out=to_surface(R)
+        if out: emit('surface_hit',multiple=int(n),source='integral_points',point=str(R),x=str(out[0]),y=str(out[1]),z=str(out[2]),x_digits=len(str(abs(out[0]))))
+except Exception as ex: emit('integral_error',multiple=int(n),error=repr(ex))
+
+rank=len(gens)
+bound=60 if rank<=2 else (18 if rank==3 else (7 if rank==4 else 3))
+rg=range(-bound,bound+1); mult=[[c*G for c in rg] for G in gens]
+tors=E.torsion_points(); tested=0; hits={}; start=time.time()
+emit('lattice_start',multiple=int(n),rank=rank,bound=bound,total=(2*bound+1)**rank*len(tors))
+for cs in itertools.product(*([rg]*rank)):
+    R=E(0)
+    for i,c in enumerate(cs): R += mult[i][c+bound]
+    for Tors in tors:
+        S=R+Tors; tested+=1; out=to_surface(S)
+        if out:
+            key=tuple(map(int,out))
+            if key not in hits:
+                hits[key]=(cs,str(Tors))
+                emit('surface_hit',multiple=int(n),source='lattice',coefficients=list(cs),torsion=str(Tors),x=str(out[0]),y=str(out[1]),z=str(out[2]),x_digits=len(str(abs(out[0]))))
+emit('lattice_done',multiple=int(n),tested=tested,hits=len(hits),seconds=time.time()-start)

research/fixed_z_search.sage

@@ -0,0 +1,101 @@
+#!/usr/bin/env sage
+from sage.all import *
+import itertools, json, sys, time
+
+proof.all(False)
+start = ZZ(sys.argv[1])
+end = ZZ(sys.argv[2])
+
+
+def emit(kind, **payload):
+    print(json.dumps({"kind": kind, **payload}, sort_keys=True), flush=True)
+
+
+def verify(x, y, z):
+    return z*z + y*y*z + x*x*x - 2 == 0
+
+
+def local_obstruction(d):
+    # Every solution has a companion with negative z, so only d=-z>0 is needed.
+    if d < 0:
+        return -1
+    fd = factor(d)
+    # x^3=2 modulo 4 or 9 is impossible.
+    for p, e in fd:
+        p, e = ZZ(p), ZZ(e)
+        if p == 2 and e >= 2:
+            return 4
+        if p == 3 and e >= 2:
+            return 9
+        if p in (2, 3) or p % 3 == 2:
+            continue
+        if power_mod(2, (p-1)//3, p) != 1:
+            return p
+    return None
+
+
+for d in range(start, end + 1):
+    if d == 0:
+        continue
+    obstruction = local_obstruction(d)
+    if obstruction is not None:
+        emit("local_skip", d=int(d), obstruction=int(obstruction))
+        continue
+
+    t0 = time.time()
+    c = d**3 * (d**2 - 2)
+    E = EllipticCurve([0, 0, 0, 0, c])
+    emit("curve_start", d=int(d), constant=str(c))
+    try:
+        rank = E.rank(proof=False)
+        gens = E.gens(proof=False)
+        emit("curve_data", d=int(d), rank=int(rank), gens=[str(P) for P in gens], seconds=time.time()-t0)
+    except Exception as exc:
+        emit("rank_error", d=int(d), error=repr(exc), seconds=time.time()-t0)
+        continue
+
+    bound = 150 if len(gens) == 1 else (35 if len(gens) == 2 else (10 if len(gens) == 3 else 4))
+    seen = set()
+    try:
+        if gens:
+            ranges = [range(-bound, bound + 1) for _ in gens]
+            for coeffs in itertools.product(*ranges):
+                if all(n == 0 for n in coeffs):
+                    continue
+                P = E(0)
+                for n, G in zip(coeffs, gens):
+                    P += n*G
+                if P.is_zero():
+                    continue
+                X, Y = QQ(P[0]), QQ(P[1])
+                if X.denominator() != 1 or Y.denominator() != 1:
+                    continue
+                X, Y = ZZ(X), ZZ(Y)
+                if X % d or Y % (d*d):
+                    continue
+                x, y, z = X//d, Y//(d*d), -d
+                if verify(x, y, z):
+                    key = (int(x), int(y), int(z))
+                    if key not in seen:
+                        seen.add(key)
+                        emit("surface_hit", source="lattice", d=int(d), coefficients=list(coeffs), x=str(x), y=str(y), z=str(z), x_digits=len(str(abs(x))))
+    except Exception as exc:
+        emit("lattice_error", d=int(d), error=repr(exc))
+
+    try:
+        pts = E.integral_points(both_signs=True)
+        emit("integral_points", d=int(d), count=len(pts))
+        for P in pts:
+            X, Y = ZZ(P[0]), ZZ(P[1])
+            if X % d or Y % (d*d):
+                continue
+            x, y, z = X//d, Y//(d*d), -d
+            assert verify(x, y, z)
+            key = (int(x), int(y), int(z))
+            if key not in seen:
+                seen.add(key)
+                emit("surface_hit", source="integral_points", d=int(d), point=str(P), x=str(x), y=str(y), z=str(z), x_digits=len(str(abs(x))))
+    except Exception as exc:
+        emit("integral_error", d=int(d), error=repr(exc))
+
+    emit("curve_done", d=int(d), hits=len(seen), seconds=time.time()-t0)

research/lambda_curve_search.sage

@@ -0,0 +1,66 @@
+#!/usr/bin/env sage
+from sage.all import *
+import itertools, json, sys, time
+proof.all(True)
+code=sys.argv[1] if len(sys.argv)>1 else '2'
+values={
+ '2':QQ(2),'3':QQ(3),'5_3':QQ(5,3),'34_15':QQ(34,15),
+ '25_41':QQ(25,41),'58_101':QQ(58,101),'87_101':QQ(87,101),
+ '9129_3943':QQ(9129,3943),'4250453_7007642':QQ(4250453,7007642)
+}
+lam=values[code]
+
+def emit(kind,**kw): print(json.dumps({'kind':kind,**kw},sort_keys=True),flush=True)
+def verify(x,y,z): return z*z+y*y*z+x*x*x-2==0
+A=lam+1/lam; B=3*(1/lam-lam); D=71*(lam-1/lam)
+a=A/8; b=B/8; c=-21*A/8; d=D/8
+# X=a*T and W=a*y gives W^2=X^3+b*X^2+a*c*X+a^2*d.
+E=EllipticCurve([0,b,0,a*c,a*a*d])
+emit('family',code=code,lam=str(lam),curve=str(E),ainvs=[str(q) for q in E.ainvs()])
+
+def to_surface(P):
+ if P.is_zero(): return None
+ X,W=QQ(P[0]),QQ(P[1]); T=X/a; y=W/a
+ x=(T*T-17)/4
+ N1=(T**3-3*T**2-21*T+71)/8
+ z=-lam*N1
+ if any(q.denominator()!=1 for q in (x,y,z)): return None
+ x,y,z=ZZ(x),ZZ(y),ZZ(z); assert verify(x,y,z); return x,y,z,T
+
+t=time.time()
+try:
+ rank=E.rank(proof=True); gens=E.gens(proof=True); tors=E.torsion_points()
+ emit('curve_data',code=code,rank=int(rank),gens=[str(P) for P in gens],torsion=[str(P) for P in tors],seconds=time.time()-t)
+except Exception as ex:
+ emit('curve_error',code=code,error=repr(ex),seconds=time.time()-t)
+ try:
+  rank=E.rank(proof=False); gens=E.gens(proof=False); tors=E.torsion_points()
+  emit('curve_data_unproved',code=code,rank=int(rank),gens=[str(P) for P in gens],torsion=[str(P) for P in tors])
+ except Exception as ex2: emit('curve_error_unproved',code=code,error=repr(ex2)); sys.exit(0)
+
+try:
+ Em=E.global_minimal_model(); iso=E.isomorphism_to(Em)
+ emit('minimal',code=code,curve=str(Em),ainvs=[str(q) for q in Em.ainvs()])
+ pts=Em.integral_points(both_signs=True)
+ emit('minimal_integral_points',code=code,count=len(pts),points=[str(P) for P in pts])
+ for Q in pts:
+  P=iso.inverse()(Q); out=to_surface(P)
+  if out: emit('surface_hit',code=code,source='minimal_integral',point=str(P),x=str(out[0]),y=str(out[1]),z=str(out[2]),T=str(out[3]),x_digits=len(str(abs(out[0]))))
+except Exception as ex: emit('integral_error',code=code,error=repr(ex))
+
+# Modest exact Mordell-Weil lattice scan, with every torsion coset.
+r=len(gens); bound=100 if r==1 else (25 if r==2 else (8 if r==3 else 3))
+rg=range(-bound,bound+1); mult=[[n*P for n in rg] for P in gens]; hits={}; tested=0; start=time.time()
+emit('lattice_start',code=code,rank=r,bound=bound,total=(2*bound+1)**r*len(tors))
+try:
+ for cs in itertools.product(*([rg]*r)):
+  R=E(0)
+  for i,n in enumerate(cs): R += mult[i][n+bound]
+  for Tor in tors:
+   tested+=1; out=to_surface(R+Tor)
+   if out:
+    key=(int(out[0]),int(out[1]),int(out[2]))
+    if key not in hits:
+     hits[key]=(cs,str(Tor)); emit('surface_hit',code=code,source='lattice',coefficients=list(cs),torsion=str(Tor),x=str(out[0]),y=str(out[1]),z=str(out[2]),T=str(out[3]),x_digits=len(str(abs(out[0]))))
+ emit('lattice_done',code=code,tested=tested,hits=len(hits),seconds=time.time()-start)
+except Exception as ex: emit('lattice_error',code=code,tested=tested,error=repr(ex),seconds=time.time()-start)