尤拉計劃748:倒丟蕃多方程

對於丟蕃多方程的倒寫
1/x² +1/y²=13/z²
設(x,y,z)是一組正整數解，其中gcd(x,y,z)=1,為原始解，
設S(n)是滿足1<=x,y,z<=n，x<=y的原始解之和。
對於n=100,有2組原始解，各是(2,3,6)和(5,90,18)，所以s(100)=124,
已知S(1000)=1470,S(100000)=2340084
求S(10¹⁶)的最後9位數

測試程式碼，在pypy3上算到10⁴需要3秒，算到10⁵需要600秒。

import math
import time

def f(n,m):
 s=0
 for x in range(2,n//3):
  for y in range(x+1,n+1,2):
   gxy=math.gcd(x,y)
   if (x*x+y*y)%m==0:
    zz=m*(x*x*y*y)/(x*x+y*y)
    if zz % 1==0:
     z=int(pow(zz,0.5))
     if z<n and z*z==zz and math.gcd(gxy,z)==1:s+=x+y+z;print(x,y,z);break
return s

t=time.time();print(f(10000,13));print(time.time()-t)

2 3 6
...
2997 3922 8586
70677
3.3921940326690674

把對y的迴圈改為對w=x*y/z的迴圈，快了一倍。 觀察w的平方根，都是質數或者質數的積(25)。

def f2(n,m):
 s=0
 for x in range(2,n//3): 
  w0=int(x*pow(2/m,0.5))+1
  if w0%2==0:w0+=1
  for w in range(w0,int(pow((x*x+n*n)/m,0.5)+1),2):
   y0=pow(m*w*w-x*x,0.5)
   if y0%1 ==0:
    y=int(y0)
    gxy=math.gcd(x,y)
    z1=y*x/w
    if z1 % 1==0:
     z=int(z1)
     if z<n and math.gcd(gxy,z)==1:s+=x+y+z;print(x,y,z,int(w**0.5));break
 return s

t=time.time();print(f2(10000,13));print(time.time()-t)

2 3 6 1
5 90 18 5
30 85 102 5
102 1037 366 17
117 598 414 13
493 918 1566 17
522 2987 1854 29
667 2958 2346 29
675 2150 2322 25
1258 4773 4386 37
1450 1725 4002 25
2173 5658 7314 41
2997 3922 8586 37
70677
1.4536311626434326

如果w不迴圈，而是隻羅列質數和質數之積，有了數量級的提高。算到10⁵縮短到1秒多。

def prime(n):
 sieve = [True] * (n//2+1)
 for number in range(3, int(n ** 0.5) + 1,2):
  if sieve[(number-1)//2]:
   for index in range(number * number, n, number*2):
    sieve[(index-1)//2] = False
 return [x for x in range(1,n,2) if sieve[(x-1)//2]]

def f3(n,m):
 s=0
 maxp=int((n/(m**0.5))**0.5)
 p=prime(maxp)
 pp=[i*j for i in p for j in p if i<=j and i*j<=maxp]
 for x in range(2,n//3):
  w0=int(x*pow(2/m,0.5))+1
  if w0%2==0:w0+=1
  for w in [j*j for j in pp]:
   if not(w>=w0 and w<int(pow((x*x+n*n)/m,0.5)+1)):continue
   y0=pow(m*w*w-x*x,0.5)
   if y0%1 ==0:
    y=int(y0)
    gxy=math.gcd(x,y)
    z1=y*x/w
    if z1 % 1==0:
     z=int(z1)
     if z<n and math.gcd(gxy,z)==1:s+=x+y+z;print(x,y,z,w**0.5);break
 return s

t=time.time();print(f3(10000,13));print(time.time()-t)
70677
0.09000492095947266
t=time.time();print(f3(100000,13));print(time.time()-t)
2052798
1.3252480030059814

這與題目公佈的正確值有差距，經過檢查，遺漏了w=125的結果。在函式f4中加上，將來隨著n的增加，還有可以加上更多的質數之積。 考慮到x，y都是w的倍數，所以x的迴圈可以放到w的迴圈內部，又有數量級的提高，可以在1秒鐘內計算10⁶，而上一種演算法需要42秒（質數表不全）。

def f4(n,m):
 s=0
 maxp=int((n/(m**0.5))**0.5)
 p=prime(maxp)
 pp=[i*j for i in p for j in p if i<=j and i*j<=maxp]
 ppp=[i*j*k for i in p for j in p for k in p if i>1 and i<=j and j<=k and i*j*k<=maxp]
 pppp=[i*j*i*j for i in p for j in p if i>1 and i<=j and i*j*i*j<=maxp]
 pp+=ppp+pppp
 for w in [j*j for j in pp]:
  for k in range(1,n//int(w**0.5)+1):
   x=k*int(w**0.5)
   w0=int(x*pow(2/m,0.5))+1
   if not(w>=w0 and w<int(pow((x*x+n*n)/m,0.5)+1)):break
   y0=pow(m*w*w-x*x,0.5)
   if y0%1 ==0:
    y=int(y0)
    gxy=math.gcd(x,y)
    z1=y*x/w
    if z1 % 1==0:
     z=int(z1)
     if z<n and math.gcd(gxy,z)==1:s+=x+y+z;print(x,y,z,int(w**0.5))
 return s

t=time.time();print(f4(100000,13));print(time.time()-t)
2340084
0.1716001033782959
t=time.time();print(f4(10**6,13));print(time.time()-t)
85126107
0.48520898818969727
t=time.time();print(f3(10**6,13));print(time.time()-t)
77588781
42.492141008377075

觀察n很大的數，以及把m=13換成其他的平方和，似乎w中不存在3 7 11,而存在1 5 13 17，這些數共同的特徵是可以表示成4*n+1的形式。

把程式碼中p=prime(maxp)改為p=[i for i in prime(maxp) if i%4==1]，結果完全不變，效能又有提高。

import random as rd
import time
import math

def ran(n):
    if n==1:return 1
    if (n%4!=1):return 0
    while True:
        a=rd.randint(1,10**6)%(n-1)+1 #1~ n-1
        b=pow(a,(n-1)//4,n);
        if(b*b%n==n-1):
          if b<=(n-1)/2 :return b
          else: return n-b

def solve(p):
    global c1
    c1+=1
    if(p==2): return [1,1];
    if(p%4!=1): return mp[-1,-1];
    A=ran(p); B=1
    tmp=1;
    M=(tmp*A*A+tmp*B*B)//p
    while M>1:
        u=(A%M+M)%M;v=(B%M+M)%M
        if(2*u>M): u-=M
        if(2*v>M): v-=M
        ta=A
        A=(tmp*u*A+tmp*v*B)//M
        B=(tmp*v*ta-tmp*u*B)//M
        M=(tmp*u*u+tmp*v*v)//M
    return [min(abs(A),abs(B)),max(abs(A),abs(B))]

def e(l,s2):
 l1=[]
 for s1 in l:
  u,v=s1[0],s1[1]
  a,b=s2[0],s2[1]
  l1.append([abs(u*a+v*b),abs(v*a-u*b)])
  l1.append([abs(u*a-v*b),abs(v*a+u*b)])
  l2=[]
 for x in l1:
  if x not in l2 and 0 not in x:
   l2.append(x)
 return l2 

def test2(sol,pl,m):
 sl=[sol[(p-1)//4] for p in pl]
 sl+=sl
 sl+=sl
 l1=[solve(m)]
 for i in range(0,len(sl)):
  l1=e(l1,sl[i])
 l2=[]
 for x in l1:
  if x not in l2 and 0 not in x:
   l2.append(x)
 return l2

def f7(n,m):
 s=0
 maxp=int((n/(m**0.5))**0.5)
 p=[i for i in prime(maxp) if i%4==1]

 #sol=[solve(px) for px in p]
 sol=[1 for i in range(maxp//4+1)]
 for px in p:sol[(px-1)//4]=solve(px)
 pp=[(i*j,[i,j]) for i in p for j in p if i<=j and i*j<=maxp]
 ppp=[(i*j*k,[i,j,k]) for i in p for j in p for k in p if i>1 and i<=j and j<=k and i*j*k<=maxp]
 pppp=[(i*j*i*j,[i,j,i,j]) for i in p for j in p if i>1 and i<=j and i*j*i*j<=maxp]
 pp+=ppp+pppp
 for w,lp in [(j[0]*j[0],j[1]) for j in pp]:
  r=test2(sol,lp,m)
  for i in r:
   x=i[0]
   y=i[1]
   gxy=math.gcd(x,y)
   z1=y*x/w
   if z1 % 1==0:
    z=int(z1)   
    if z<n and math.gcd(gxy,z)==1:s+=x+y+z;#print(x,y,z,int(w**0.5))
 return s

c1=0
t=time.time();print(f7(10**7,13));print(time.time()-t,c1)

上述是採用費馬降階法求整數化為平方和的程式碼