题目66：考察丢番图方程x^2 − Dy^2 = 1
Diophantine equationConsider 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 的值。
最笨最暴力的解法：
#coding:utf-8
maxx = 1
for D in range(2,1001):
print (D)
if (D**0.5) % 1 != 0:
y = 1
while True:
x = (D * y * y + 1) ** 0.5
if x % 1 == 0:
if x > maxx:
maxx = int(x)
break
else:
y += 1
print (maxx)
2178548422
2178548422 * 2178548422 - 778 * 78104745 * 78104745 = 1 # encoding:utf-8
# x**2 - D*y**2 = 1
from time import time
from math import sqrt
def euler066():
max_x, r_D, r_y = 0, 0, 0
for D in range(2, 1001):
if int(sqrt(D)) == sqrt(D):
continue
y = 0
flag = False
while not flag:
y += 1
x = sqrt(D * (y ** 2) + 1)
if x == int(x):
print(D, int(x), y)
y = 1
if int(x) > max_x:
max_x, r_D, r_y = int(x), D, y
flag = True
break
print(max_x, r_D, r_y)
if __name__ == '__main__':
start = time()
euler066()
print('cost %.6f sec' % (time() - start))
跑到D=500多等不了了.... 采用连分数法解Pell方程
答案：16421658242965910275055840472270471049， 661
# encoding:utf-8
'''
Pell方程
X^2 - DY^2=1
运用连分数理论
'''
from math import sqrt
#sqrt(n) = a0+1/(a1+1/(a2+1/....+1/an))
#the sequence is repeating, a1,a2,...ax,a1,a2,...ax
#√13=
#√13=, period=5
def getContinuedFractionsOfSqrt(n):
a0 = int(sqrt(n))
if n == a0*a0:
return a0, []
fractions = []
c = 1
b = a0
fdict = {}
while True:
#(sqrt(n)-b)/c ^ -1 ==> a1 + (sqrt(n)-b1)/c1
c1 = (n - b*b)//c
tb = b
#now (sqrt(n)+tb) / c1 == a1 + (sqrt(n)-b1)/c1
# a1c1-b1==tb, b1 <= a0 ===> a1 <= (a0+tb)/c1 => a1 = (a0+tb)//c1
a1 = (a0+tb)//c1
b1 = a1*c1-tb
if fdict.get((a1,b1,c1), 0) > 0:
break
fdict[(a1,b1,c1)] = 1
fractions.append(a1)
c = c1
b = b1
return a0, fractions
def calculteContinuedFractions(a0, fractions, fcount):
#return franction a/b (a,b)
i = (fcount-1) % len(fractions)
a,b = 1, fractions
for cnt in range(fcount-1):
#print(a, b)
i = (i-1) % len(fractions)
an1 = fractions #A(n-1)
#1/(an1 + (a/b)) ===> b / (an1*b+a)
a,b = b, an1*b+a
#add A0, a0+(a/b) ==> (a0*b+a)/b
a += a0*b
return a,b
def solvePell(D): #may not min X,Y
a0, fractions = getContinuedFractionsOfSqrt(D)
m = len(fractions)
if m == 0:
return 0,0
#print(a0, fractions)
if m%2 == 0: #even
#x0=P(mn-1),y0=Q(mn-1), n=1,2,3...
#n=1
x, y = calculteContinuedFractions(a0, fractions, m-1)
else: #odd
#x0=P(2mn-1),y0=Q(2mn-1), n=0,1,2,3...
#n=1
x, y = calculteContinuedFractions(a0, fractions, 2*m-1)
return x,y
#print(solvePell(13)) #649, 180
#print(solvePell(778)) #2178548422, 78104745
maxx = 0
saveD = 0
for D in range(2, 1000+1):
x,y = solvePell(D)
if x > maxx:
maxx = x
saveD = D
print(maxx, saveD) #16421658242965910275055840472270471049 661 def findx(n):
m = int(n**0.5)
ifm * 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
本帖最后由 guosl 于 2021-3-21 11:58 编辑
方程x^2-Dy^2=1的解可以通过将D的平方根展开成连分数来得到。时间可以控制到0.01秒以下。
#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <mpirxx.h>
using namespace std;
int gcd(int x, int y)
{
if (y == 0)
return x;
int z = x % y;
while (z != 0)
{
x = y;
y = z;
z = x % y;
}
return y;
}
int gcd3(int x, int y, int z)
{
int a = gcd(x, y);
a = gcd(a, z);
return a;
}
struct stTriples
{
int a, b, c;
bool operator==(const stTriples &t) const
{
return (a == t.a && b == t.b && c == t.c);
}
};
void get_triples(int nD, const stTriples &tr1, int &n, stTriples &tr2)
{
int a1 = tr1.a * tr1.c, b1 = -tr1.b * tr1.c, c1 = tr1.a * tr1.a * nD - tr1.b * tr1.b;
n = int((sqrt(nD) * double(a1) + double(b1)) / double(c1));
b1 -= n * c1;
int k = gcd3(abs(a1), abs(b1), abs(c1));
a1 /= k;
b1 /= k;
c1 /= k;
tr2.a = a1;
tr2.b = b1;
tr2.c = c1;
}
int get_cycle(int nD, vector<int>& vI)
{
int p = 0;
int q = int(sqrt((double)nD));
vI.push_back(q);
if (q * q == nD)
return 0;
vector<stTriples> vTriples;
stTriples tr1, tr2;
tr1.a = 1;
tr1.b = -q;
tr1.c = 1;
int n;
while (true)
{
get_triples(nD, tr1, n, tr2);
auto itr = find(vTriples.begin(), vTriples.end(), tr2);
if (itr != vTriples.end())
{
p = int(vTriples.end() - itr);
break;
}
vI.push_back(n);
vTriples.push_back(tr2);
tr1 = tr2;
}
return p;
}
int main(void)
{
mpz_class xMax = 0;
int nD = 0;
for (int k = 2; k < 1000; ++k)
{
vector<int> vI;
int n = get_cycle(k, vI);
if (n == 0)
continue;
int m = vI.size();
int s = m - n;
int kk = n - 1;
if ((n & 1) == 1)
{
kk = 2 * n - 1;
for (int i = s; i < m; ++i)
vI.push_back(vI);
}
mpz_class p, q;
p = vI, p = vI * vI + 1;
q = 1, q = vI;
for (int i = 2; i <= kk; ++i)
{
p = vI * p + p;
q = vI * q + q;
p = p;
p = p;
q = q;
q = q;
}
if (p > xMax)
{
xMax = p;
nD = k;
}
}
cout << nD << endl;
return 0;
}
页:
[1]