当前位置:首页 » Project Euler » 详细页

Diophantine equation 题号:66 难度: 25 中英对照

Consider quadratic Diophantine equations of the form:

x2 – Dy2 = 1

For example, when D=13, the minimal solution in x is 6492 – 13×1802 = 1.

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:

32 – 2×22 = 1
22 – 3×12 = 1
92 – 5×42 = 1
52 – 6×22 = 1
82 – 7×32 = 1

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.

Solution

此问题实际上属于['佩尔方程'](http://blog.sina.com.cn/s/blog_5d06e2390100ll92.html)问题。解决的思路我们可以先求出 $1000​$ 以内,每个 $D​$ 值对应的满足方程$x^2-Dy^2=1​$ 的 $x​$ 的最小值。再比较这些 $x​$ 值,从而找到 $x​$ 的最大值,及其对应的 $D​$ 值。(注意,当 $D​$ 是一个完全平方数时,是不存在正整数解的,跳过这些 $D​$ 值)。 对于一个给定的非完全平方数 $D$ ,由[64题](http://www.maths.date/cn/pe/solve/64) ,可知,$D$ 可以分解为连分数的形式 $D=[a_0;a_1,a_2,…a_{r-1},a_r]$,周期为 $r$ 。对于这个 $D$,对应佩尔方程的$x^2-Dy^2=1$的 $x$ 取最小值的解,可以使用如下方法计算,方法的具体说明请参考['佩尔方程'](http://blog.sina.com.cn/s/blog_5d06e2390100ll92.html): $$ p_0=a_0, \quad p_1=a_1p_0+1,\quad ...,\quad p_{r-1}=a_{r-1}p_{r-2}+p_{r-3} $$ $$ q_0=1, \quad q_1=a_1,\quad q_2=a_2q_1+q_0...,\quad q_{r-1}=a_{r-1}q_{r-2}+q_{r-3} $$ $$ MinX=p_{r-1}\quad 当r为偶数 $$ $$ MinX=p_{r-1}^2+Dq_{r-1}^2\quad 当r为奇数 $$

Code

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public final class p66 {
	
	public static void main(String[] args) {
		long start = System.nanoTime();
		int result = run();
        long end = System.nanoTime();
		System.out.println( result );
		System.out.println (( end - start )/1000000 +"ms");
	}
	
	
	//对每个D值计算最小解x,再找出其中最大的x,及其对应的D值
	public static int run() {
		int minN = -1;
		BigInteger maxX = BigInteger.ZERO;
		for (int n = 2; n <= 1000; n++) {
			if (isSquare(n))
				continue;
			BigInteger x = smallestSolutionX(n);
			if (x.compareTo(maxX) > 0) {
				minN = n;
				maxX = x;
			}
		}
		return minN;
	}
	
	
	//返回最小解的x>0,使得存在某个值y,满足 x^2 - n y^2 = 1.
	//n不能是个完全平方数
	private static BigInteger smallestSolutionX(int n) {
		List<BigInteger>[] contFrac = sqrtToContinuedFraction(n);

		List<BigInteger> temp = new ArrayList<>();
		temp.addAll(contFrac[0]);
		temp.addAll(contFrac[1].subList(0, contFrac[1].size() - 1));
		Fraction val = new Fraction(temp.get(temp.size() - 1), BigInteger.ONE);
		
		for (int i = temp.size() - 2; i >= 0; i--){
			val = new Fraction(val.denominator, val.numerator).add(new Fraction(temp.get(i), BigInteger.ONE));
		}
		
		if (contFrac[1].size() % 2 == 0)
			return val.numerator;
		else
			return val.numerator.pow(2).add(val.denominator.pow(2).multiply(BigInteger.valueOf(n)));
	}
	
	
	// 返回sqrt(n)的周期连续的分数项
	//result[0]是最小的非周期的前缀,result[1]是最小的周期尾部 例如sqrt(23)=[4;(1,3,1,8)],result[0]=4,result[1]=[1,3,1,8]
	@SuppressWarnings("unchecked")
	private static List<BigInteger>[] sqrtToContinuedFraction(int n) {
		List<BigInteger> terms = new ArrayList<>();
		Map<QuadraticSurd,Integer> seen = new HashMap<>();
		QuadraticSurd val = new QuadraticSurd(BigInteger.ZERO, BigInteger.ONE, BigInteger.ONE, BigInteger.valueOf(n));
		do {
			seen.put(val, seen.size());
			BigInteger flr = val.floor(); //对应的是a值
			terms.add(flr);
			val = val.subtract(new QuadraticSurd(flr, BigInteger.ZERO, BigInteger.ONE, val.d)).reciprocal();
		} while (!seen.containsKey(val));
		return new List[]{terms.subList(0, seen.get(val)), terms.subList(seen.get(val), terms.size())};
	}
	
	public static boolean isSquare(int x) {
		if (x < 0)
			return false;
		int y = (int) Math.sqrt(x);
		return y * y == x;
	}
	
	// 表示为(a + b * sqrt(d)) / c。 d不能为完全平方数
	private static class QuadraticSurd {
		
		public final BigInteger a, b, c, d;
		
		
		public QuadraticSurd(BigInteger a, BigInteger b, BigInteger c, BigInteger d) {
			if (c.signum() == 0)
				throw new IllegalArgumentException();
			
			// Simplify
			if (c.signum() == -1) {
				a = a.negate();
				b = b.negate();
				c = c.negate();
			}
			BigInteger gcd = a.gcd(b).gcd(c);
			if (!gcd.equals(BigInteger.ONE)) {
				a = a.divide(gcd);
				b = b.divide(gcd);
				c = c.divide(gcd);
			}
			
			this.a = a;
			this.b = b;
			this.c = c;
			this.d = d;
		}
		
		
		public QuadraticSurd subtract(QuadraticSurd other) {
			if (!d.equals(other.d))
				throw new IllegalArgumentException();
			return new QuadraticSurd(a.multiply(other.c).subtract(other.a.multiply(c)), b.multiply(other.c).subtract(other.b.multiply(c)), c.multiply(other.c), d);
		}
		
		
		public QuadraticSurd reciprocal() {
			return new QuadraticSurd(a.multiply(c).negate(), b.multiply(c), b.multiply(b).multiply(d).subtract(a.multiply(a)), d);
		}
		
		
		public BigInteger floor() {
			BigInteger temp = sqrt(b.multiply(b).multiply(d));
			if (b.signum() == -1)
				temp = temp.add(BigInteger.ONE).negate();
			temp = temp.add(a);
			if (temp.signum() == -1)
				temp = temp.subtract(c.subtract(BigInteger.ONE));
			return temp.divide(c);
		}
		
		public static BigInteger sqrt(BigInteger x) {
			if (x.signum() == -1)
				throw new IllegalArgumentException("Square root of negative number");
			BigInteger y = BigInteger.ZERO;
			for (int i = (x.bitLength() - 1) / 2; i >= 0; i--) {
				y = y.setBit(i);
				if (y.multiply(y).compareTo(x) > 0)
					y = y.clearBit(i);
			}
			return y;
		}
		
		public boolean equals(Object obj) {
			if (!(obj instanceof QuadraticSurd))
				return false;
			else {
				QuadraticSurd other = (QuadraticSurd)obj;
				return a.equals(other.a) && b.equals(other.b) && c.equals(other.c) && d.equals(other.d);
			}
		}
		
		
		public int hashCode() {
			return a.hashCode() + b.hashCode() + c.hashCode() + d.hashCode();
		}
		
		
		public String toString() {
			return String.format("(%d + %d*sqrt(%d)) / %d", a, b, d, c);
		}
		
	}
	
	//分数的类和相关计算
	final static class Fraction implements Comparable<Fraction> {
		
		public final BigInteger numerator;    // 分子一定与分母互质
		public final BigInteger denominator;  
		
		
		public Fraction(BigInteger numer, BigInteger denom) {
			if (denom.signum() == 0)
				throw new ArithmeticException("Division by zero");
			
			//化为标准形式
			if (denom.signum() == -1) {
				numer = numer.negate();
				denom = denom.negate();
			}
			BigInteger gcd = numer.gcd(denom);
			if (!gcd.equals(BigInteger.ONE)) {
				numer = numer.divide(gcd);
				denom = denom.divide(gcd);
			}
			
			numerator = numer;
			denominator = denom;
		}
		
		
		public Fraction add(Fraction other) {
			return new Fraction(numerator.multiply(other.denominator).add(other.numerator.multiply(denominator)), denominator.multiply(other.denominator));
		}
		
		
		public Fraction subtract(Fraction other) {
			return new Fraction(numerator.multiply(other.denominator).subtract(other.numerator.multiply(denominator)), denominator.multiply(other.denominator));
		}
		
		
		public Fraction multiply(Fraction other) {
			return new Fraction(numerator.multiply(other.numerator), denominator.multiply(other.denominator));
		}
		
		
		public Fraction divide(Fraction other) {
			return new Fraction(numerator.multiply(other.denominator), denominator.multiply(other.numerator));
		}
		
		
		public boolean equals(Object obj) {
			if (!(obj instanceof Fraction))
				return false;
			Fraction other = (Fraction)obj;
			return numerator.equals(other.numerator) && denominator.equals(other.denominator);
		}
		
		
		public int compareTo(Fraction other) {
			return numerator.multiply(other.denominator).compareTo(other.numerator.multiply(denominator));
		}
		
		
		public int hashCode() {
			return numerator.hashCode() + denominator.hashCode();
		}
		
		
		public String toString() {
			return numerator + "/" + denominator;
		}
		
	}

}
661
133ms