projecteuler No.66 Diophantine equation

题目链接：

Problem 66 Diophantine equation

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

通过人数：9417

题目分析：

这道题实质上需要面临的是求佩尔方程的最小正整数解。至于接下来的寻找使得最小正整数解最大的方程系数就是比较容易的了。

而关于佩尔方程的解法，这篇来自网络的文章较详细地介绍了连分数解法。

于是就有了具体的解题方法了，因为涉及高精度运算，所以选用java语言。

解题过程：

零、首先先贴一下完整代码

【详细分析在后面】

package test;

import java.math.*;


public class Main 
{
	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;
	}
	// return p1 + 1/p2
}


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, Main.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, Main.precision, BigDecimal.ROUND_FLOOR);
		}
		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, Main.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, Main.precision, BigDecimal.ROUND_FLOOR);
		}
		return ans;
	}
}

这部分具体的算法在之前的一道题：projecteuler No.64 Odd period square roots中较详细地讲了，在此就不多说了。

这两部分封装在了这个类中。第一个函数以Main中定义的常量precision为精度使用牛顿迭代法开平方，第二个函数对大于1的无理数In求前l位的连分展开并返回。

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);

这段为主函数对这个类的调用，将sqrt(D_i)的连分展开前100项存在a中。

二、分数类

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;
	}
	// return p1 + 1/p2
}

这个类是一个非常简单的基于大整数的分数类。Compute函数返回p1+1/p2，用于连分数化简为普通分数。

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

这段是主函数对这个类的调用，根据连分展开a[]，化简其前i项并存入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;
}

由于佩尔方程最小正解为连分展开中某项的分子分母，于是将其带入尝试。

若result为0，则为最小解。判断并更新当前Max值并保存对应的D值。

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

这句话对在连分数前100项没有解的方程报警。

最后没有报警，输出答案为661.

分析：这道题由于涉及精度比较高，1s内跑不完，但1min内可以跑完，所以时间上还是没有什么问题的。

用python的话好像对高精度处理会方便不少...可惜我不会....最近开学了，也没有时间自学了...

以上只是我做题时的解法。

如果有更好的解法、更好的思路，欢迎评论讨论~O(∩_∩)O~