Quux.java

package square;

import static java.lang.Double.MAX_EXPONENT;
import static java.lang.Double.MIN_EXPONENT;
import static java.lang.Double.POSITIVE_INFINITY;
import static java.lang.Double.doubleToRawLongBits;
import static java.lang.Double.longBitsToDouble;
import static java.lang.Math.getExponent;
import static java.lang.Math.rint;
import static square.Quux.Guava.DoubleUtils.SIGNIFICAND_BITS;
import static square.Quux.Guava.DoubleUtils.getSignificand;
import static square.Quux.Guava.DoubleUtils.isFinite;

import java.math.BigInteger;
import java.security.SecureRandom;

import square.Quux.Guava.BigIntegerMath;

public final class Quux {
	static final class NFinder {
		private final BigInteger n = new BigInteger(SquareRoot.BITS,
				new SecureRandom()).pow(2);

		private static final int CALIBRATE_LOOPS = 100000;
		private static final int WARMUP_LOOPS = 200000;
		private static final int TEST_LOOPS = 100000;

		// at least twice as fast as slow
		private static final double FAST_RATIO = 0.50;

		// 75% of the slowest seen in calibration is still slow
		private static final double SLOW_RATIO = 0.75;

		// warm up the JIT
		public void warmUp() {
			System.out.println("Warming up JIT");
			final SecureRandom random = new SecureRandom();
			for (int i = 0; i < WARMUP_LOOPS; i++) {
				BigInteger root = new BigInteger(SquareRoot.BITS, random);
				calibrateAnswer(root);
				SquareRoot.answer(root);
			}
		}

		private long calibrateTime() {
			long min = Long.MAX_VALUE;
			for (int i = 0; i < CALIBRATE_LOOPS; i++) {
				final BigInteger root = BigInteger.ONE
						.shiftLeft(n.bitLength() - 2);
				long time = System.nanoTime();
				calibrateAnswer(root);
				time = System.nanoTime() - time;
				if (time <= 0) {
					throw new Error("Calibration failed. Answered in " + time + " ns.");
				}
				min = Math.min(min, time);

			}
			return min;
		}

		private void calibrateAnswer(BigInteger root) {
			if (n.divide(root).equals(root)) {
				// The goal is to reach this line
				System.out.println("Square root!");
			}
		}

		public BigInteger findN() {
			System.out.println("Calibrating");
			final long minSlowTime = calibrateTime();
			
			System.out.println("Solving");

			BigInteger lowerBound = BigInteger.ZERO;
			final int maxBit = (SquareRoot.BITS + 1) * 2;
			for (int i = 0; i <= maxBit; i++) {
				System.out.println("Solving bit " + (maxBit - i));
				final BigInteger bit = BigInteger.ONE.shiftLeft(maxBit - i);
				final BigInteger candidate = lowerBound.add(bit);
				if (isSlow(candidate, minSlowTime)) {
					lowerBound = candidate;
				}
			}

			return lowerBound.add(BigInteger.ONE);
		}

		private static boolean isSlow(BigInteger candidate, long minSlowTime) {
			for (;;) {
				long min = Long.MAX_VALUE;
				for (int i = 0; i < TEST_LOOPS; i++) {
					long time = System.nanoTime();
					SquareRoot.answer(candidate);
					time = System.nanoTime() - time;
					if (time < minSlowTime * FAST_RATIO) {
						System.out.println("" + time + " < " + (minSlowTime)
								+ " * " + FAST_RATIO);
						return false;
					}
					min = Math.min(time, min);
				}

				if (min >= minSlowTime * SLOW_RATIO) {
					System.out.println("" + min + " >= " + minSlowTime + " * "
							+ SLOW_RATIO);
					return true;
				}

				System.out.println("Result was inconclusive. min=" + min
						+ "; minSlowTime=" + minSlowTime);
			}
		}
	}

	/*
	 * The classes contained in Guava are all derived directly from the Apache
	 * licensed classes of the same names in the Google Guava project.
	 */
	static final class Guava {
		private Guava() {
		}

		/**
		 * A class for arithmetic on values of type {@code BigInteger}.
		 * 
		 * <p>
		 * The implementations of many methods in this class are based on
		 * material from Henry S. Warren, Jr.'s <i>Hacker's Delight</i>,
		 * (Addison Wesley, 2002).
		 * 
		 * <p>
		 * Similar functionality for {@code int} and for {@code long} can be
		 * found in {@link IntMath} and {@link LongMath} respectively.
		 * 
		 * @author Louis Wasserman
		 * @since 11.0
		 */
		public final static class BigIntegerMath {
			/**
			 * Returns the base-2 logarithm of {@code x}, rounded according to
			 * the specified rounding mode.
			 */
			@SuppressWarnings("fallthrough")
			// TODO(kevinb): remove after this warning is disabled globally
			private static int log2Floor(BigInteger x) {
				int logFloor = x.bitLength() - 1;
				return logFloor;
			}

			/*
			 * The maximum number of bits in a square root for which we'll
			 * precompute an explicit half power of two. This can be any value,
			 * but higher values incur more class load time and linearly
			 * increasing memory consumption.
			 */
			static final int SQRT2_PRECOMPUTE_THRESHOLD = 256;

			public static BigInteger sqrtFloor(BigInteger x) {
				/*
				 * Adapted from Hacker's Delight, Figure 11-1.
				 * 
				 * Using DoubleUtils.bigToDouble, getting a double approximation
				 * of x is extremely fast, and then we can get a double
				 * approximation of the square root. Then, we iteratively
				 * improve this guess with an application of Newton's method,
				 * which sets guess := (guess + (x / guess)) / 2. This iteration
				 * has the following two properties:
				 * 
				 * a) every iteration (except potentially the first) has guess
				 * >= floor(sqrt(x)). This is because guess' is the arithmetic
				 * mean of guess and x / guess, sqrt(x) is the geometric mean,
				 * and the arithmetic mean is always higher than the geometric
				 * mean.
				 * 
				 * b) this iteration converges to floor(sqrt(x)). In fact, the
				 * number of correct digits doubles with each iteration, so this
				 * algorithm takes O(log(digits)) iterations.
				 * 
				 * We start out with a double-precision approximation, which may
				 * be higher or lower than the true value. Therefore, we perform
				 * at least one Newton iteration to get a guess that's
				 * definitely >= floor(sqrt(x)), and then continue the iteration
				 * until we reach a fixed point.
				 */
				BigInteger sqrt0;
				int log2 = log2Floor(x);
				if (log2 < Double.MAX_EXPONENT) {
					sqrt0 = sqrtApproxWithDoubles(x);
				} else {
					int shift = (log2 - DoubleUtils.SIGNIFICAND_BITS) & ~1; // even!
					/*
					 * We have that x / 2^shift < 2^54. Our initial
					 * approximation to sqrtFloor(x) will be 2^(shift/2) *
					 * sqrtApproxWithDoubles(x / 2^shift).
					 */
					sqrt0 = sqrtApproxWithDoubles(x.shiftRight(shift))
							.shiftLeft(shift >> 1);
				}
				BigInteger sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1);
				if (sqrt0.equals(sqrt1)) {
					return sqrt0;
				}
				do {
					sqrt0 = sqrt1;
					sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1);
				} while (sqrt1.compareTo(sqrt0) < 0);
				return sqrt0;
			}

			private static BigInteger sqrtApproxWithDoubles(BigInteger x) {
				return DoubleMath.roundToBigIntegerHalfEven(Math
						.sqrt(DoubleUtils.bigToDouble(x)));
			}

			private BigIntegerMath() {
			}
		}

		/**
		 * A class for arithmetic on doubles that is not covered by
		 * {@link java.lang.Math}.
		 * 
		 * @author Louis Wasserman
		 * @since 11.0
		 */
		static final class DoubleMath {
			/*
			 * This method returns a value y such that rounding y DOWN (towards
			 * zero) gives the same result as rounding x according to the
			 * specified mode.
			 */
			private static double roundIntermediateHalfEven(double x) {
				if (!isFinite(x)) {
					throw new ArithmeticException("input is infinite or NaN");
				}
				return rint(x);
			}

			private static final double MIN_LONG_AS_DOUBLE = -0x1p63;
			/*
			 * We cannot store Long.MAX_VALUE as a double without losing
			 * precision. Instead, we store Long.MAX_VALUE + 1 ==
			 * -Long.MIN_VALUE, and then offset all comparisons by 1.
			 */
			private static final double MAX_LONG_AS_DOUBLE_PLUS_ONE = 0x1p63;

			/**
			 * Returns the {@code BigInteger} value that is equal to {@code x}
			 * rounded with the specified rounding mode, if possible.
			 */
			static BigInteger roundToBigIntegerHalfEven(double x) {
				x = roundIntermediateHalfEven(x);
				if (MIN_LONG_AS_DOUBLE - x < 1.0
						& x < MAX_LONG_AS_DOUBLE_PLUS_ONE) {
					return BigInteger.valueOf((long) x);
				}
				int exponent = getExponent(x);
				long significand = getSignificand(x);
				BigInteger result = BigInteger.valueOf(significand).shiftLeft(
						exponent - SIGNIFICAND_BITS);
				return (x < 0) ? result.negate() : result;
			}

			private DoubleMath() {
			}
		}

		/**
		 * Utilities for {@code double} primitives.
		 * 
		 * @author Louis Wasserman
		 */
		final static class DoubleUtils {
			private DoubleUtils() {
			}

			// The mask for the significand, according to the {@link
			// Double#doubleToRawLongBits(double)} spec.
			static final long SIGNIFICAND_MASK = 0x000fffffffffffffL;

			// The mask for the exponent, according to the {@link
			// Double#doubleToRawLongBits(double)} spec.
			static final long EXPONENT_MASK = 0x7ff0000000000000L;

			// The mask for the sign, according to the {@link
			// Double#doubleToRawLongBits(double)} spec.
			static final long SIGN_MASK = 0x8000000000000000L;

			static final int SIGNIFICAND_BITS = 52;

			static final int EXPONENT_BIAS = 1023;

			/**
			 * The implicit 1 bit that is omitted in significands of normal
			 * doubles.
			 */
			static final long IMPLICIT_BIT = SIGNIFICAND_MASK + 1;

			static long getSignificand(double d) {
				if (!isFinite(d))
					throw new IllegalArgumentException("not a normal value");
				int exponent = getExponent(d);
				long bits = doubleToRawLongBits(d);
				bits &= SIGNIFICAND_MASK;
				return (exponent == MIN_EXPONENT - 1) ? bits << 1 : bits
						| IMPLICIT_BIT;
			}

			static boolean isFinite(double d) {
				return getExponent(d) <= MAX_EXPONENT;
			}

			static double bigToDouble(BigInteger x) {
				// This is an extremely fast implementation of
				// BigInteger.doubleValue(). JDK patch pending.
				BigInteger absX = x.abs();
				int exponent = absX.bitLength() - 1;
				// exponent == floor(log2(abs(x)))
				if (exponent < Long.SIZE - 1) {
					return x.longValue();
				} else if (exponent > MAX_EXPONENT) {
					return x.signum() * POSITIVE_INFINITY;
				}

				/*
				 * We need the top SIGNIFICAND_BITS + 1 bits, including the
				 * "implicit" one bit. To make rounding easier, we pick out the
				 * top SIGNIFICAND_BITS + 2 bits, so we have one to help us
				 * round up or down. twiceSignifFloor will contain the top
				 * SIGNIFICAND_BITS + 2 bits, and signifFloor the top
				 * SIGNIFICAND_BITS + 1.
				 * 
				 * It helps to consider the real number signif = absX *
				 * 2^(SIGNIFICAND_BITS - exponent).
				 */
				int shift = exponent - SIGNIFICAND_BITS - 1;
				long twiceSignifFloor = absX.shiftRight(shift).longValue();
				long signifFloor = twiceSignifFloor >> 1;
				signifFloor &= SIGNIFICAND_MASK; // remove the implied bit

				/*
				 * We round up if either the fractional part of signif is
				 * strictly greater than 0.5 (which is true if the 0.5 bit is
				 * set and any lower bit is set), or if the fractional part of
				 * signif is >= 0.5 and signifFloor is odd (which is true if
				 * both the 0.5 bit and the 1 bit are set).
				 */
				boolean increment = (twiceSignifFloor & 1) != 0
						&& ((signifFloor & 1) != 0 || absX.getLowestSetBit() < shift);
				long signifRounded = increment ? signifFloor + 1 : signifFloor;
				long bits = (long) ((exponent + EXPONENT_BIAS)) << SIGNIFICAND_BITS;
				bits += signifRounded;
				/*
				 * If signifRounded == 2^53, we'd need to set all of the
				 * significand bits to zero and add 1 to the exponent. This is
				 * exactly the behavior we get from just adding signifRounded to
				 * bits directly. If the exponent is MAX_DOUBLE_EXPONENT, we
				 * round up (correctly) to Double.POSITIVE_INFINITY.
				 */
				bits |= x.signum() & SIGN_MASK;
				return longBitsToDouble(bits);
			}
		}
	}

	public static void main(String[] args) {
		new NFinder().warmUp();
		for (;;) {
			final BigInteger n = new NFinder().findN();
			System.out.println("Identified N: " + n);
			final BigInteger root = BigIntegerMath.sqrtFloor(n);
			if (root.pow(2).equals(n)) {
				System.out.println("High confidence solution.");
				SquareRoot.answer(root);
				System.exit(0);
			}
			System.out.println("Invalid solution; trying again.");
		}
	}
}