ProjectEuler Problem 66

gistfile1.java

import java.math.BigInteger;
import java.util.Arrays;
import java.util.HashMap;
import java.util.Map;

// @see http://www37.atwiki.jp/uwicoder/pages/1639.html
// x^2-Dy^2=1の特殊解を求める。
// (s^2-D)t^2+(2s)t+1=0, s=0
// (-D)t^2+1=0
// t=1/√Dの連分数展開からconvergentsをつくって-D(num)^2+(den)^2=1になる最初の(num,den)をとってくればよい。
public class E66 {
	public static void solve()
	{
		BigInteger max = BigInteger.ZERO;
		int maxd = -1;
		for(int D = 1;D <= 1000;D++){
			int sq = (int)Math.sqrt(D);
			if(sq * sq == D)continue;
			long[] res = new long[200];
			int pos = makeCF(0, 1, D, D, res);
			
			BigInteger num = BigInteger.ONE;
			BigInteger nump = BigInteger.ZERO;
			BigInteger den = BigInteger.ZERO;
			BigInteger denp = BigInteger.ONE;
			
			for(int i = 0;;i++){
				if(i > 0 && res[i] == 0)i = pos;
				BigInteger numn = num.multiply(b(res[i])).add(nump);
				BigInteger denn = den.multiply(b(res[i])).add(denp);
				nump = num; num = numn;
				denp = den; den = denn;
				if(i > 0 && b(-D).multiply(num).multiply(num)
						.add(den.multiply(den)).equals(BigInteger.ONE)){
					break;
				}
			}
			//x=den, y=num
//			tr(D, num, den, den.multiply(den).subtract(b(D).multiply(num.multiply(num))));
			if(den.compareTo(max) > 0){
				max = den;
				maxd = D;
			}
		}
		tr(max, maxd); // [16421658242965910275055840472270471049, 661]
	}
	
	static BigInteger b(long n){ return BigInteger.valueOf(n); }
	
	public static int makeCF(long a, long b, long c, long d, long[] ret)
	{
		Map<Long, Integer> cache = new HashMap<Long, Integer>();
		for(int pos = 0;;pos++){
			long code = (long)c*123456789+(long)b*987654321+a;
			if(cache.containsKey(code)){
				return cache.get(code);
			}
			cache.put(code, pos);
			long p = (int)Math.floor((a+b*Math.sqrt(d))/c);
			long den = (a-c*p)*(a-c*p)-b*b*d;
			if(den == 0)return 0;
			long numa = c*(a-c*p);
			long numb = -c*b;
			long g = gcd(gcd(Math.abs(numa), Math.abs(numb)), Math.abs(den));
			numa /= g;
			numb /= g;
			den /= g;
			ret[pos] = p;
			a = numa; b = numb; c = den;
		}
	}
	
	public static long gcd(long a, long b) {
		if(a == 0)return b;
		if(b == 0)return a;
		while (b > 0){
			long c = a;
			a = b;
			b = c % b;
		}
		return a;
	}
	
	public static void main(String[] args) {
		long s = System.currentTimeMillis();
		solve();
		long g = System.currentTimeMillis();
		tr(g-s+"ms");
	}
	
	static void tr(Object... o) { System.out.println(Arrays.deepToString(o)); }
}