Create a gist now

Instantly share code, notes, and snippets.

@pburka /Quux.java Secret forked from anonymous/Quux.java
Created Mar 25, 2013

What would you like to do?
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.");
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment