题目66：考察丢番图方程x^2 − Dy^2 = 1

欧拉计划

Diophantine equation

Consider quadratic Diophantine equations of the form:



For example, when D=13, the minimal solution in x is

It can be assumed that there are no solutions in positive integers when D is square.

By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:



Hence, by considering minimal solutions in x for D ≤ 7, the largest x is obtained when D=5.

Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is obtained.

题目：

考虑如下形式的二次丢番图方程：



例如当 D=13 时, x 的最小解是

可以认为当 D 时平方数时方程无正整数解。

通过寻找当 D = {2, 3, 5, 6, 7} 时 x 的最小解，我们得到：



因此对于 D ≤ 7, x 的最小解的最大值在 D=5 时取到。

找出 D ≤ 1000 中使得 x 的最小值取到最大的 D 的值。

jerryxjr1220

最笨最暴力的解法：

1. #coding:utf-8
   
2. maxx = 1
   
3. for D in range(2,1001):
   
4.         print (D)
   
5.         if (D**0.5) % 1 != 0:
   
6.                 y = 1
   
7.                 while True:
   
8.                         x = (D * y * y + 1) ** 0.5
   
9.                         if x % 1 == 0:
   
10.                                 if x > maxx:
    
11.                                         maxx = int(x)
    
12.                                 break
    
13.                         else:
    
14.                                 y += 1
    
15. print (maxx)

复制代码

2178548422
[Finished in 2751.0s]

jerryxjr1220

2178548422 * 2178548422 - 778 * 78104745 * 78104745 = 1

芒果加黄桃

1. # encoding:utf-8
   
2. # x**2 - D*y**2 = 1
   
3. from time import time
   
4. from math import sqrt
   
5. def euler066():
   
6.     max_x, r_D, r_y = 0, 0, 0
   
7.     for D in range(2, 1001):
   
8.         if int(sqrt(D)) == sqrt(D):
   
9.             continue
   
10.         y = 0
    
11.         flag = False
    
12.         while not flag:
    
13.             y += 1
    
14.             x = sqrt(D * (y ** 2) + 1)
    
15.             if x == int(x):
    
16.                 print(D, int(x), y)
    
17.                 y = 1
    
18.                 if int(x) > max_x:
    
19.                     max_x, r_D, r_y = int(x), D, y
    
20.                 flag = True
    
21.                 break
    
22.     print(max_x, r_D, r_y)
    
23. if __name__ == '__main__':
    
24.     start = time()
    
25.     euler066()
    
26.     print('cost %.6f sec' % (time() - start))
    

复制代码

跑到D=500多等不了了....

gunjang

采用连分数法解Pell方程
答案：16421658242965910275055840472270471049， 661

1. # encoding:utf-8
   
2. '''
   
3. Pell方程
   
4. X^2 - DY^2=1
   
5. 运用连分数理论
   
6. '''
   
7. from math import sqrt
   
8. 
   
9. #sqrt(n) = a0+1/(a1+1/(a2+1/....+1/an))
   
10. #the sequence is repeating, a1,a2,...ax,a1,a2,...ax
    
11. #√13=[3;(1,1,1,1,6,1,1,1,1,6,1,1,1,1,6,...)]
    
12. #√13=[3;(1,1,1,1,6)], period=5
    
13. def getContinuedFractionsOfSqrt(n):
    
14.         a0 = int(sqrt(n))
    
15.         if n == a0*a0:
    
16.                 return a0, []
    
17. 
    
18.         fractions = []
    
19.         c = 1
    
20.         b = a0
    
21.         fdict = {}
    
22.         while True:
    
23.                 #(sqrt(n)-b)/c ^ -1 ==> a1 + (sqrt(n)-b1)/c1
    
24.                 c1 = (n - b*b)//c
    
25.                 tb = b
    
26.                 #now (sqrt(n)+tb) / c1 == a1 + (sqrt(n)-b1)/c1
    
27.                 # a1c1-b1==tb, b1 <= a0 ===> a1 <= (a0+tb)/c1 => a1 = (a0+tb)//c1
    
28.                 a1 = (a0+tb)//c1
    
29.                 b1 = a1*c1-tb
    
30.                 if fdict.get((a1,b1,c1), 0) > 0:
    
31.                         break
    
32.                 fdict[(a1,b1,c1)] = 1
    
33.                 fractions.append(a1)
    
34.                 c = c1
    
35.                 b = b1
    
36. 
    
37.         return a0, fractions
    
38. 
    
39. def calculteContinuedFractions(a0, fractions, fcount):
    
40.         #return franction a/b (a,b)
    
41.         i = (fcount-1) % len(fractions)
    
42.         a,b = 1, fractions[i]
    
43.         for cnt in range(fcount-1):
    
44.                 #print(a, b)
    
45.                 i = (i-1) % len(fractions)
    
46.                 an1 = fractions[i] #A(n-1)
    
47.                 #1/(an1 + (a/b)) ===> b / (an1*b+a)
    
48.                 a,b = b, an1*b+a
    
49.         #add A0, a0+(a/b) ==> (a0*b+a)/b
    
50.         a += a0*b
    
51.         return a,b       
    
52. 
    
53.        
    
54. def solvePell(D): #may not min X,Y
    
55.         a0, fractions = getContinuedFractionsOfSqrt(D)
    
56.         m = len(fractions)
    
57.         if m == 0:
    
58.                 return 0,0
    
59.         #print(a0, fractions)
    
60. 
    
61.         if m%2 == 0: #even
    
62.         #x0=P(mn-1),y0=Q(mn-1), n=1,2,3...
    
63.         #n=1
    
64.                 x, y = calculteContinuedFractions(a0, fractions, m-1)
    
65.         else: #odd
    
66.         #x0=P(2mn-1),y0=Q(2mn-1), n=0,1,2,3...
    
67.         #n=1
    
68.                 x, y = calculteContinuedFractions(a0, fractions, 2*m-1)       
    
69.         return x,y
    
70. 
    
71. #print(solvePell(13)) #649, 180
    
72. #print(solvePell(778)) #2178548422, 78104745
    
73. 
    
74. maxx = 0
    
75. saveD = 0
    
76. for D in range(2, 1000+1):
    
77.         x,y = solvePell(D)
    
78.         if x > maxx:
    
79.                 maxx = x
    
80.                 saveD = D
    
81. print(maxx, saveD) #16421658242965910275055840472270471049 661

复制代码

grf1973

def findx(n):   
    m = int(n**0.5)
    if  m * m == n :
        return 0
    ha,ka = 1, 0
    d,a,b,c = m, n, -m,1
    hb = d; kb = 1
    if hb * hb - n * kb * kb == 1:
       return hb
    while True:
        nC = a - b * b
        nC = nC / c
        nd = int((a**0.5 - b) / nC)
        nb = -b - nd * nC
        c = nC; d = nd; b = nb
        hc = d*hb+ha
        kc = d * kb + ka
        ha,ka,hb,kb = hb,kb,hc,kc        
        if hc*hc==kc*kc*n+1:
            return hc         
    return 0


ans = 0; tmp = 0
for i in range(1,1001):
    z = findx(i)
    if z > tmp:
        tmp=z
        ans=i

print(ans, tmp)

661 16421658242965910275055840472270471049
[Finished in 0.1s]

guosl

方程x^2-Dy^2=1的解可以通过将D的平方根展开成连分数来得到。时间可以控制到0.01秒以下。

1. #include <iostream>
   
2. #include <vector>
   
3. #include <cmath>
   
4. #include <algorithm>
   
5. #include <mpirxx.h>
   
6. using namespace std;
   
7. 
   
8. int gcd(int x, int y)
   
9. {
   
10.   if (y == 0)
    
11.     return x;
    
12.   int z = x % y;
    
13.   while (z != 0)
    
14.   {
    
15.     x = y;
    
16.     y = z;
    
17.     z = x % y;
    
18.   }
    
19.   return y;
    
20. }
    
21. 
    
22. int gcd3(int x, int y, int z)
    
23. {
    
24.   int a = gcd(x, y);
    
25.   a = gcd(a, z);
    
26.   return a;
    
27. }
    
28. 
    
29. struct stTriples
    
30. {
    
31.   int a, b, c;
    
32.   bool operator==(const stTriples &t) const
    
33.   {
    
34.     return (a == t.a && b == t.b && c == t.c);
    
35.   }
    
36. };
    
37. 
    
38. void get_triples(int nD, const stTriples &tr1, int &n, stTriples &tr2)
    
39. {
    
40.   int a1 = tr1.a * tr1.c, b1 = -tr1.b * tr1.c, c1 = tr1.a * tr1.a * nD - tr1.b * tr1.b;
    
41.   n = int((sqrt(nD) * double(a1) + double(b1)) / double(c1));
    
42.   b1 -= n * c1;
    
43.   int k = gcd3(abs(a1), abs(b1), abs(c1));
    
44.   a1 /= k;
    
45.   b1 /= k;
    
46.   c1 /= k;
    
47.   tr2.a = a1;
    
48.   tr2.b = b1;
    
49.   tr2.c = c1;
    
50. }
    
51. 
    
52. int get_cycle(int nD, vector<int>& vI)
    
53. {
    
54.   int p = 0;
    
55.   int q = int(sqrt((double)nD));
    
56.   vI.push_back(q);
    
57.   if (q * q == nD)
    
58.     return 0;
    
59.   vector<stTriples> vTriples;
    
60.   stTriples tr1, tr2;
    
61.   tr1.a = 1;
    
62.   tr1.b = -q;
    
63.   tr1.c = 1;
    
64.   int n;
    
65.   while (true)
    
66.   {
    
67.     get_triples(nD, tr1, n, tr2);
    
68.     auto itr = find(vTriples.begin(), vTriples.end(), tr2);
    
69.     if (itr != vTriples.end())
    
70.     {
    
71.       p = int(vTriples.end() - itr);
    
72.       break;
    
73.     }
    
74.     vI.push_back(n);
    
75.     vTriples.push_back(tr2);
    
76.     tr1 = tr2;
    
77.   }
    
78.   return p;
    
79. }
    
80. 
    
81. int main(void)
    
82. {
    
83.   mpz_class xMax = 0;
    
84.   int nD = 0;
    
85.   for (int k = 2; k < 1000; ++k)
    
86.   {
    
87.     vector<int> vI;
    
88.     int n = get_cycle(k, vI);
    
89.     if (n == 0)
    
90.       continue;
    
91.     int m = vI.size();
    
92.     int s = m - n;
    
93.     int kk = n - 1;
    
94.     if ((n & 1) == 1)
    
95.     {
    
96.       kk = 2 * n - 1;
    
97.       for (int i = s; i < m; ++i)
    
98.         vI.push_back(vI[i]);
    
99.     }
    
100.     mpz_class p[3], q[3];
     
101.     p[0] = vI[0], p[1] = vI[0] * vI[1] + 1;
     
102.     q[0] = 1, q[1] = vI[1];
     
103.     for (int i = 2; i <= kk; ++i)
     
104.     {
     
105.       p[2] = vI[i] * p[1] + p[0];
     
106.       q[2] = vI[i] * q[1] + q[0];
     
107.       p[0] = p[1];
     
108.       p[1] = p[2];
     
109.       q[0] = q[1];
     
110.       q[1] = q[2];
     
111.     }
     
112.     if (p[1] > xMax)
     
113.     {
     
114.       xMax = p[1];
     
115.       nD = k;
     
116.     }
     
117.   }
     
118.   cout << nD << endl;
     
119.   return 0;
     
120. }
     

复制代码