Skip to content

Instantly share code, notes, and snippets.

@joriki
Last active January 11, 2024 08:52
Show Gist options
  • Save joriki/614f0355b5da6f27cf9efc0d56d30624 to your computer and use it in GitHub Desktop.
Save joriki/614f0355b5da6f27cf9efc0d56d30624 to your computer and use it in GitHub Desktop.
Estimate the probability that the triangle formed by three uniformly random points on the sphere contains its circumcentre; see https://math.stackexchange.com/questions/4842652.
public class Question4840933 {
final static long ntrials = 2000000;
public static void main(String [] args) {
double [] [] x = new double [3] [3];
double [] c = new double [3];
long count = 0;
for (long n = 0;n < ntrials;n++) {
// the first point can be fixed arbitrarily
x [0] [0] = 0;
x [0] [1] = 0;
x [0] [2] = 1;
// the azimuth of the second point can be fixed arbitrarily
double cos1 = 1 - 2 * Math.random();
double sin1 = Math.sqrt (1 - cos1 * cos1);
x [1] [0] = 0;
x [1] [1] = sin1;
x [1] [2] = cos1;
// the third point has to be in general position
double cos2 = 1 - 2 * Math.random();
double sin2 = Math.sqrt (1 - cos2 * cos2);
double phi = 2 * Math.PI * Math.random();
double cos = Math.cos(phi);
double sin = Math.sin(phi);
x [2] [0] = sin2 * cos;
x [2] [1] = sin2 * sin;
x [2] [2] = cos2;
// form the differences from the first point
x [1] [2]--;
x [2] [2]--;
// form the cross product of the differences
for (int i = 0;i < 3;i++)
c [i] = x [1] [(i + 1) % 3] * x [2] [(i + 2) % 3] - x [2] [(i + 1) % 3] * x [1] [(i + 2) % 3];
normalize (c);
// restore the original points
x [1] [2]++;
x [2] [2]++;
// project the points along the cross product
for (int j = 0;j < 3;j++) {
double dot = dot (x [j],c);
for (int i = 0;i < 3;i++)
x [j] [i] -= dot * c [i];
normalize (x [j]);
}
// check whether one of the angles is reflex
double angle = 0;
for (int j = 0;j < 3;j++)
angle += Math.acos(dot (x [j],x [(j + 1) % 3]));
if (angle < 2 * Math.PI - 1e-11)
count++;
}
System.out.println(count / (double) ntrials);
}
static void normalize (double [] x) {
double norm = 0;
for (int i = 0;i < 3;i++)
norm += x [i] * x [i];
norm = 1 / Math.sqrt(norm);
for (int i = 0;i < 3;i++)
x [i] *= norm;
}
static double dot (double [] x,double [] y) {
double dot = 0;
for (int i = 0;i < 3;i++)
dot += x [i] * y [i];
return dot;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment