题目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

最笨最暴力的解法：

#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
[Finished in 2751.0s]

jerryxjr1220

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多等不了了....

gunjang

采用连分数法解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=[3;(1,1,1,1,6,1,1,1,1,6,1,1,1,1,6,...)]
#√13=[3;(1,1,1,1,6)], 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[i]
        for cnt in range(fcount-1):
                #print(a, b)
                i = (i-1) % len(fractions)
                an1 = fractions[i] #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

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秒以下。

#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[i]);
    }
    mpz_class p[3], q[3];
    p[0] = vI[0], p[1] = vI[0] * vI[1] + 1;
    q[0] = 1, q[1] = vI[1];
    for (int i = 2; i <= kk; ++i)
    {
      p[2] = vI[i] * p[1] + p[0];
      q[2] = vI[i] * q[1] + q[0];
      p[0] = p[1];
      p[1] = p[2];
      q[0] = q[1];
      q[1] = q[2];
    }
    if (p[1] > xMax)
    {
      xMax = p[1];
      nD = k;
    }
  }
  cout << nD << endl;
  return 0;
}