欧拉工程第66题：Diophantine equation

题目链接

脑补知识：佩尔方差

上面说的貌似很明白，最小的i，对应最小的解

然而我理解成，一个循环的解了，然后就是搞不对，后来，仔细看+手工推导发现了问题。i从0开始变量，知道第一个满足等式的解就是最小解。

问题转化为求根号n的连分数问题，分子就是x，分母就是y

要求的分子，分母，问题又转化为：根号n的连分数表示，对，求出其连分数表示就OK了

先求出a的序列是什么？

第64题，就是求a的序列的。

a求出来了，要求出分子分母的表达式。

第65题，就是已经知道了a的序列，求分子，当然也可以求分母的

分子，分母求出来了，在验证：X*X-D*Y*Y=1时候就是最小解

问题真是一环套一环的。

Python程序：

import time as time 

start = time.time()

def getD(N):
    x_max ,y_max= 0,0
    D = 0
    x,y = 0,0
    for S in range(2,N+1):
        x,y = resolve(S)
        if x>x_max:
            x_max ,y_max= x,y
            D = S
    return D,x_max,y_max

def resolve(S):
    m = 0
    d = 1
    a0 = int(S**0.5)
    if a0*a0 == S :return -1,-1;
    a= a0
    li = [a]
    x,y = 1,1
    while x*x-S*y*y!=1:
        m = d*a - m
        d = (S - m*m)/d
        a = int((a0 + m)/d)
        li.append(a)
        x = getX(li)
        y = getY(li)
#     print li
    return x,y;

def getX(li):
    x0 = 1
    x1 = li[0]
    li = li[1:]
    for l in li:
        x = l * x1 + x0
        x0 = x1
        x1 = x
    return x 

def getY(li):
    y0 = 0
    y1 = 1
    li = li[1:]
    for l in li:
        y = l * y1 + y0
        y0 = y1
        y1 = y
    return y 

if __name__ == ‘__main__‘:

    start = time.time()
    N = 1000
    D ,x_max,y_max= getD(N)
    print "running time={0}seconds,D={1},x_max={2},y_max={3}".format(time.time()-start,D,x_max,y_max)

求的是最小解X的最大值时候的D，答案是661

然而：

x_max=16421658242965910275055840472270471049

y_max=638728478116949861246791167518480580

这个值好大的

附一：python程序：

from math import sqrt
from time import time

def prefect_sqrt(n):
    return int(sqrt(n))**2 == n 

def floor_root(n):
    return int(sqrt(n))

def chakravals(n):
    x_max = 0
    for d in range(2,n+1):
        if not prefect_sqrt(d):
            p1 = floor_root(d)
            q1 = 1
            m1 = p1**2 - d
#             print -p1,(-p1) % abs(m1),(-p1) % abs(m1)
            if (-p1) % abs(m1) ==0:
                x1 = abs(m1)
            else:
                x1 = (-p1) % abs(m1)

            while m1!=1:
                p0 = p1
                q0 = q1
                m0 = m1
                x0 = x1

                p1 = (p0 * x0 +d *q0)/abs(m0)
                q1 = (p0 + x0)/abs(m0)

                m1 = (x0**2 -d)/m0 

                if (-x0)%abs(m1) ==0:
                    x1 = abs(x0)
                else:
                    x1 = (-x0)%abs(m1)
            if p1>x_max:
                x_max = p1
                d_max = d
                print "d= %04d  x = %d"%(d_max,x_max)
    print 

if __name__==‘__main__‘:
    start = time()
    chakravals(1000)
    end = time()
    print "time elapse=%f"%(end - start)

Java程序：

这个跑的好慢的

package project61;
import java.math.*;  

public class P66
{
    public static final int precision = 500;
    public static void main( String args[] )
    {
        BigInteger Max = new BigInteger("0");
        int ans = 0;  

outer:  for (int D_i = 2; D_i <= 1000; D_i++)
        {
            BigDecimal D = new BigDecimal(D_i);
            BigDecimal SD = calculation.Sqrt(D);  

            BigDecimal SD_i = SD.setScale(0, BigDecimal.ROUND_FLOOR);
            if (SD_i.multiply(SD_i).equals(D))
                continue;  

            int a[] = calculation.toContinuedFraction(SD, 100);  

            for (int i = 1; i < 100; i++)
            {
                Fraction temp = new Fraction(a[i],1);  

                for (int j = i - 1; j >= 0; j--)
                    temp = Fraction.Compute(a[j], temp);  

                BigInteger y_2 = temp.denominator.multiply(temp.denominator);
                BigInteger x_2 = temp.numerator.multiply(temp.numerator);
                BigInteger result = x_2.subtract(y_2.multiply(D.toBigIntegerExact())).subtract(BigInteger.ONE);  

                if (result.equals(BigInteger.ZERO))
                {
                    if (temp.numerator.compareTo(Max) > 0)
                    {
                        Max = temp.numerator;
                        ans = D_i;
                    }  

                    continue outer;
                }  

            }
            System.out.print("Warning!\n");  

        }
        System.out.print(ans+"\n");
    }
}  

class Fraction
{
    public BigInteger numerator;
    public BigInteger denominator;  

    private Fraction()
    {  

    }
    public Fraction (int numerator, int denominator)
    {
        this.numerator = BigInteger.valueOf(numerator);
        this.denominator = BigInteger.valueOf(denominator);
    }
    public static Fraction Compute(int p1, Fraction p2)
    {
        Fraction ans = new Fraction();
        ans.numerator = p2.denominator.add(p2.numerator.multiply(BigInteger.valueOf(p1)));
        ans.denominator = p2.numerator.add(BigInteger.ZERO);
        return ans;
    }
}  

class calculation
{
    private static final BigDecimal N0 = new BigDecimal(0);
    private static final BigDecimal N1 = new BigDecimal(1);
    private static final BigDecimal N2 = new BigDecimal(2);  

    public static BigDecimal Sqrt(BigDecimal In)
    {
        BigDecimal N = new BigDecimal(1);
        while(true)
        {
            BigDecimal NN = N.multiply(N);
            NN = NN.add(In);
            NN = NN.divide(N2);
            NN = NN.divide(N, P66.precision, BigDecimal.ROUND_FLOOR);  

            if (NN.equals(N))
                break;  

            N = NN;
        }  

        return N;
    }
    public static int[] toContinuedFraction(BigDecimal In, int l)
    {
        int ans[] = new int[l];  

        BigDecimal temp = In.add(N0);
        for (int i = 0; i < l; i++)
        {
            ans[i] = Integer.valueOf(temp.setScale(0, BigDecimal.ROUND_FLOOR).toString()).intValue();
            temp = temp.subtract(temp.setScale(0, BigDecimal.ROUND_FLOOR));
            temp = N1.divide(temp, P66.precision, BigDecimal.ROUND_FLOOR);
        }
        return ans;
    }
}

最后两个程序在网上复制过来的