simulate raindrops covering a table; see http://math.stackexchange.com/questions/176383

Question176383.java

import java.util.List;
import java.util.ArrayList;
import java.util.Random;

public class Question176383 {
    public final static Centre origin = new Centre (0,0);

    static List<Centre> centres = new ArrayList<Centre> ();

    static enum Ternary {
	FALSE,UNKNOWN,TRUE;

	public static Ternary and (Ternary t1,Ternary t2) {
	    return values () [Math.min (t1.ordinal (),t2.ordinal ())];
	}

	public static Ternary or (Ternary t1,Ternary t2) {
	    return values () [Math.max (t1.ordinal (),t2.ordinal ())];
	}

	public Ternary negate () {
	    return values () [2 - ordinal ()];
	}
    }

    static class Range {
	double min;
	double max;
	
	public Range (double min,double max) {
	    this.min = min;
	    this.max = max;
	}
    }

    static Random random = new Random (1);

    static class Centre {
	double x;
	double y;

	public Centre (double x,double y) {
	    this.x = x;
	    this.y = y;
	}

	public static Centre random () {
	    double x,y;
	    do {
		x = 1 - 2 * random.nextDouble ();
		y = 1 - 2 * random.nextDouble ();
	    }
	    while (x * x + y * y > 1);
	    return new Centre (x,y);
	}

	public String toString () {
	    return "(" + x + ";" + y + ")";
	}
    }

    static class Square {
	Square [] parts; // four sub-squares, only allocated when needed

	double x; // left
	double y; // top
	double d; // side length

	double dry = Double.MAX_VALUE; // least radius at which the square is partially covered by exactly one drop
	double one = Double.MAX_VALUE; // least radius at which the square is partially covered by more than one drop
	double wet = Double.MAX_VALUE; // least radius at which the square is completely covered

	public Square (double x,double y,double d) {
	    this.x = x;
	    this.y = y;
	    this.d = d;

	    Range range = getRange (origin);
	    // treat squares beyond the table as always covered
	    if (range.min > 1)
		wet = 0;
	    // treat border squares as partially covered at any radius
	    else if ((range.min - 1) * (range.max - 1) < 0)
		dry = 0;
	}

	public boolean isCovered (double r) {
	    // successively search deeper until we can decide whether the square is covered
	    for (int depth = 1;;depth++) {
		Ternary isCovered = isCovered (r,depth);
		if (isCovered != Ternary.UNKNOWN)
		    return isCovered == Ternary.TRUE;
	    }
	}

	int lastCentre;

	public Ternary isCovered (double r,int depth) {
	    // update radius bounds with centres added since last call
	    while (lastCentre < centres.size ()) {
		Range range = getRange (centres.get (lastCentre++));
		wet = Math.min (wet,range.max);
		one = Math.min (one,Math.min (wet,Math.max (dry,range.min)));
		dry = Math.min (dry,range.min);
	    }


	    if (r > wet) return Ternary.TRUE;
	    if (r < one) return Ternary.FALSE;
	    if (depth == 0) return Ternary.UNKNOWN;

	    // coverage is unknown and depth is non-zero: descend into sub-squares

	    if (parts == null) {
		parts = new Square [4];
		double d2 = d / 2;
		for (int i = 0,k = 0;i < 2;i++)
		    for (int j = 0;j < 2;j++)
			parts [k++] = new Square (x + i * d2,y + j * d2,d2);
	    }

	    double allWet = 0; // least radius at which all sub-squares are completely covered

	    // if a subsquare returns false, return false immediately
	    // if a subsquare returns unknown, remember not to return true, but keep going in case we can get false from another sub-square to end the search
	    Ternary allCovered = Ternary.TRUE;
	    for (Square part : parts) {
		Ternary partIsCovered = part.isCovered (r,depth - 1);

		if (partIsCovered == Ternary.FALSE)
		    return Ternary.FALSE;
		allCovered = Ternary.and (allCovered,partIsCovered);
		allWet = Math.max (allWet,part.wet);
	    }

	    // no square return false

	    if (allWet > wet)
		throw new Error ("inconsistent bounds");
	    wet = allWet; // this may or may not improve the bound

	    return allCovered; // return true if all sub-squares returned true
	}

	// get the range of radii in which this square would be partially covered by a drop at centre c
	public Range getRange (Centre c) {
	    double min = Double.MAX_VALUE;
	    double max = 0;

	    // loop over the four corners
	    for (int i = 0,k = 0;i < 2;i++)
                for (int j = 0;j < 2;j++,k++) {
		    double dx = x + i * d - c.x;
		    double dy = y + j * d - c.y;
		    double dist = Math.sqrt (dx * dx + dy * dy);
		    max = Math.max (max,dist);
		    min = Math.min (min,dist);
		}

	    // the maximal distance from the centre is always attained at a corner, but the minimal distance might be attained on the border or in the interior

	    // is the centre within the square's x and/or y ranges?
	    boolean inX = c.x > x && c.x < x + d;
	    boolean inY = c.y > y && c.y < y + d;

	    if (inX && inY) // the centre is inside the square
		min = 0;
	    else if (inX) // the minimal distance is attained on a horizontal edge
		min = c.x < x ? x - c.x : c.x - (x + d);
	    else if (inY) // the minimal distance is attained on a vertical edge
		min = c.y < y ? y - c.y : c.y - (y + d);
	    return new Range (min,max);
	}

	public String toString () {
	    return "(" + x + "," + y + ") [" + d + "], partial : " + (dry == 0);
	}
    }

    public static void main (String [] args) {
	double r0 = 1;    // maximal drop radius
	double r1 = .1;  // minimal drop radius
	int nsteps = 100;
	double step = r0 / nsteps;
	double [] counts = new double [nsteps];  // sum of circle counts required
	double [] counts2 = new double [nsteps]; // sum of squares to estimate the variance

	final int ntrials = 100;

	for (int n = 0;n < ntrials;n++) {
	    System.out.println ("trial " + n + " of " + ntrials);
	    int k = 0;
	    double r = r0;
	
	    Square root = new Square (-1,-1,2);
	    centres.clear ();

	    outer:
	    for (;;) {
		// loop while the current drop radius suffices to cover the table
		while (root.isCovered (r)) {
		    counts [k] += centres.size ();
		    counts2 [k] += centres.size () * centres.size ();
		    k++;
		    r -= step;
		    if (r < r1 - step / 2)
			break outer;
		}
		// the drop radius doesn't suffice to cover the table; add another centre
		centres.add (Centre.random ());
	    }
	}
	
	// print the results
	double r = r0;
	for (int k = 0;r > r1 - step / 2;k++,r -= step) {
	    counts [k] /= ntrials;
	    counts2 [k] /= ntrials;
	    System.out.println (r + " : " + counts [k] + " +- " + Math.sqrt ((counts2 [k] - counts [k] * counts [k]) / (ntrials - 1)));
	}
    }
}