Probabilities of a slightly altered dice

In art class a few months ago, I made a dice out of clay. This dice isn’t quite your typical dice. Its structure is tampered with ever so slightly, so that just \frac{1}{64} of the material is removed.

Let us model this crude, imperfect dice as a hypothetical, perfect dice. Suppose this dice is represented as a cube, composed of 64 smaller ‘cubelets’ (imagine the cubelets as 1x1x1 and the dice as 4x4x4). Furthermore, we imagine four ‘layers’ of cubelets, such that each layer is a 4×4 square.

In this dice, we remove one of the 16 cubelets from the second layer, like this:

(keep in mind that this is still a solid cube)

Assuming that the cube is made out of uniform density material (and that the missing cubelet has no mass), what is the probability of landing on each of the six sides when the cube is thrown?

A mathematical interpretation

I will solve a problem which may or may not be equivalent to throwing a dice. Instead of throwing the dice, we place the dice on a table with a random orientation, and let go. We then observe which side the dice lands on.

This appears to be a legit and identical interpretation, but there is a flaw. With an uneven dice, it is not guaranteed that each orientation would appear with uniform likelihood. Nevertheless, for our purposes we ignore this and we assume that this is an identical problem.

Which side will it land on? The center of gravity of the dice must be directly on one of the six sides; this side is the side that the dice will fall onto.

Consider a sphere of some radius, centered at the cube’s center of gravity and entirely contained in the cube:

Lines are drawn from the center of gravity to the eight vertices of the cube, essentially dividing the sphere into six pyramid-like parts.

The chance of landing on a side of the cube is proportional to the side’s surface area when ‘projected’ onto the surface of the sphere. This is because the side can be divided into infinitely many small ‘pyramids’, each one of which corresponds to the dice landing on that side when left to drop at that specific angle.

Such a portion of the surface area is called a solid angle, which is measured in steradians (sr). Whereas the entire circle is denoted as 2 \pi radians, the entire sphere is denoted by 4 \pi sr. Likewise, one hemisphere would be 2 \pi sr, and so on.

There is an efficient formula for calculating the solid angle of a tetrahedral projection onto a sphere, given by Oosterom and Strackee. If \vec{a}, \vec{b}, \vec{c} are vectors of three vertices of a tetrahedron (in relationship to the fourth), then the subtended solid angle is given by:

\tan(\frac{1}{2} \Omega) = \frac{|\vec{a} \vec{b} \vec{c}|}{abc + c(\vec{a} \cdot \vec{b}) + b(\vec{a} \cdot \vec{c}) + a(\vec{b} \cdot \vec{c})}

Here a, b, c are magnitudes of the vectors, the dot represents a dot product, and |\vec{a} \vec{b} \vec{c}| is the determinant of the scalar product of the vectors.

The magnitude of the 3-D euclidean distance:

a = \sqrt{a_x^2 + a_y^2 + a_z^2}

The dot product of two vectors is:

\vec{a} \cdot \vec{b} = a_x b_x + a_y b_y + a_z b_x


|\vec{a} \vec{b} \vec{c}| = |a_x b_y c_z + b_x c_y a_z + c_x a_y b_z - a_x b_z c_y - b_x c_z a_y - c_x a_z b_y|

Calculating the probability

Let us position the cube on a 3D coordinate, such that the corners are the origin and (4,4,4). The empty cube is the cube with center (2.5,2.5,2.5).

To find the center of gravity, we just take the averages of the centroids of all the smaller cubes (the missing cube contributing 0 to the sum). This turns out to be (k,k,k) with k= \frac{251}{126} or about 1.9921.

We have k = \frac{251}{126}; now we let k' = 4-k. To facilitate vector manipulations, we will now transpose the cube such that the center of gravity is at (0,0,0):

Because of symmetry, it is obvious that three of the six sides should have one equal probability, and the other three sides to have a different, but again equal, probability. Thus it is not strictly necessary to go through the calculations for all six sides; we only need to do them for two of the sides.

A face forms a solid angle for a square pyramid, not a tetrahedron (for which we have a formula for). To solve this, we note that a square pyramid can be divided into two tetrahedrons, and the solid angle for it is simply the sum of solid angles for the two tetrahedrons.

Now by doing these calculations, we get a probability of 16.740% of landing on each of the three heavier sides, and 16.594% on each of the three lighter sides (sides closer to the missing cubelet).

By comparison, the probability of landing on any side on a fair dice is 16.667%. As only \frac{1}{64} of the dice’s mass is removed, the resulting difference is slightly less than 1%.

Project Euler 299: Three similar triangles

Being the last Project Euler problem before the summer break, Problem 299 is quite an interesting problem. Solving it involves both geometry and number theory.

Points A, B, C, D are represented on a coordinate plane as (a,0), (b,0), (0,c), and (0,d) respectively, all of which have integer coordinates.

P is a point on AC with integer coordinates such that triangles DCP, DBP, and PAB are similar.

It can be shown that in order for the triangles to be similar, a=c.

For b+d < 100,000,000, how many triples (a,b,d) exist such that point P exists?

Initial observations

Before we can start coding, there is some geometry work to be done.

It is given that it must be necessary for a=c. Why is this? Suppose that a \neq c. Then \angle OAC \neq \angle OCA. Next, \angle DCP \neq \angle PAB. It is obvious that triangles DCP and PAB cannot be similar if \angle DCP \neq \angle PAB.

Since \angle COA is a right angle and \triangle COA is isosceles, it follows that \angle OCA = \angle OAC = 45^\circ. So \angle DCP = \angle PAB = 135^\circ and also \angle DPB = 135^\circ.

Working backwards, we know that triangles DCP, DPB, and PAB are all similar, Then \angle CDP = \angle PDB and \angle DBP = \angle PBA; lines DP and PB are angular bisectors of \triangle DOB, thus P is the incenter of \triangle DOB.

So the distance from P to the three sides OB, OD, and DB are equal. This also means P can be represented as (i,i) since its x and y coordinates are equal.

We should note at this point that there is an additional case, where \angle CDP = \angle DBP and \angle ABP = \angle BDP. Then \angle DPC = \angle BDP, so lines AC and BD are parallel. However, in this case P is no longer the incenter.

We shall consider the two as separate cases, and refer to them as the incenter and parallel case respectively.

Incenter case

We first consider the incenter case. Note that in this case, a is uniquely determined by b and d. For any pair (b,d), there is only one possible AC passing through the incenter.

We need to find pairs of (b,d) such that there exists a point (i,i) where the distance from the point (i,i) to BD is i (and that i is integral).

Line BD can be expressed by the equation

y = -\frac{d}{b}x + d ,

or in standard form,

dx + by - bd = 0 .

Recall the distance from point to line formula giving the distance between a point (m,n) to line Ax + By + C = 0:

d = \frac{|Am+Bn+C|}{\sqrt{A^2 + B^2}} .

By substitution, we have

i = \frac{|di+bi-bd|}{\sqrt{b^2 + d^2}} .


i^2 (b^2 + d^2) = (di+bi-bd)^2 .

It is necessary and sufficient for b^2+d^2 to be a perfect square, as then i will be uniquely determined and will be an integer.

Thus the incenter case reduces to finding all pairs (b,d) for b+d < 100,000,000 where b^2 + d^2 is a perfect square.

Parallel case

Now we consider the case when AC || BD.

Let X be the circumcenter of \triangle DPB:

Notice that here there are actually two possible places for P. We can ignore this fact for now.

Let T be a point (not shown) such that DPBC is cyclic. Then \angle DTB = 45^\circ because \angle DPB = 135^\circ, therefore \angle DXB = 90^\circ.

As a consequence of the parallelism of AC and BD, b and d must be equal. Since \angle X is a right angle, it follows that the coordinates of X is (b,b), and that DOBX is a square.

The circle around X is centered at (b,b) and has a radius of b, thus its equation is

(x-b)^2 + (y-b)^2 = b^2 .

The line AC has an equation of

y=c-x .

By substitution into the circle equation:

(x-b)^2 + (c-x-b)^2 = b^2 ,


2x^2 + (-2c)x + (b^2+c^2-2bc)=0 ;

Applying the quadratic formula and dividing by 2 gives

x = \frac{c \pm \sqrt{c^2 - 2(b-c)^2}}{2} .

Here it is sufficient for c^2 - 2(b-c)^2 to be a perfect square, as then x will be an integer.

We prove this by using a parity argument: if c is odd, then c^2 is odd as well, and the expression inside the radical is odd; Supposing that it is a perfect square, the square root of that is odd, and when added to c, makes an even number. A similar argument applies if c is even.

We can substitute f for b-c giving the perfect square

c^2 - 2f^2.

If we let q^2 be the perfect square, we get

q^2 + 2f^2 = c^2.

Essentially the problem reduces down to finding integral solutions to the above equation, with the limit set to

2(c+f) < 100,000,000.

Writing the code

We are now ready to write a program to enumerate integer pairs to satisfy our equations.

We will start with the incenter case, which is somewhat more basic and easier to deal with.

Recall that we have derived this equation (replacing the variables with x, y, z):

x^2 + y^2 = z^2 ,

with the limit being on the sum of x+y. Obviously, this is a pythagorean triple. Enumerating pythagorean triples can be done very efficiently.

If m and n are coprime integers with an odd sum, and with m<n, then the primitive pythagorean triples can be parameterized by the formulas:

x = n^2 - m^2

y = 2mn

z = n^2 + m^2

Just by using this formula and little else, primitive pythagorean triples triples can be enumerated very quickly.

It is not very difficult to enumerate the non-primitive triples, either. Suppose we have generated the triple (3,4,5). To count the number of similar triples with x+y < 100000, sum the x and y values of the triple, which is in this case 7. Then divide the limit by 7, which in this case is \frac{100,000}{7}, so we have 14285 such triangles.

A nearly identical approach can be used to enumerate pairs for the parallel case. Here we have

x^2 + 2y^2 = z^2 ,

which is nearly identical to the previous equation and can be parameterized as:

x = n^2 - 2m^2

y = 2mn

z = n^2 + 2m^2

This case is slightly different from the previous case, in the sense that we no longer require x to be positive, so we do not need the restriction that m<n. Additionally, we need to divide by y+z, instead of x+y as we did before.

This is easy to implement in Java:

public class Main{
  final static long L = 100000000;

  static long gcd(long m, long n) {
    if (m < n) { long t = m; m = n; n = t; }
    long r = m % n;
    if (r == 0) return n;
    else return gcd(n, r); }

  static long incenterCase(){
    long count = 0;
    for(long n = 1; n < L/2; n++)
      for(long m = 1; m < n; m++){
        if((m+n) % 2 == 0) continue;
        if(gcd(m,n)!=1) continue;
        long b = n*n - m*m;
        long d = 2*n*m;
        long sum = b+d;
        if(sum >= L) break;
        if(b == d) count += L/sum;
        else count += 2*(L/sum); }
    return count; }

  static long parallelCase(){
    long count = 0;
    for(long n = 1; n < L; n+=2)
      for(long m = 1; m < L; m++){
        if(gcd(m,n)!=1) continue;
        long g = 2*n*m;
        long a = n*n + 2*m*m;
        long b = g+a;
        if(b > L/2) break;
        count += (L-1)/(2*b); }
    return count; }

  public static void main(String[] args) {
    System.out.println(incenterCase() + parallelCase()); }

This code generates the correct answer in about 25 seconds on my machine.