Minimum quadrilateral inscribed in a square

A problem that I’ve seen lately reduces to the following problem:

We have a square, and we put a point on each side of the square. Then we connect the four points to create a quadrilateral. How can we make this quadrilateral have the smallest possible perimeter?

Intuitively, you may believe that this natural, obvious configuration should produce the least perimeter:

Attempt with Calculus

How can we prove that this indeed gives us the smallest possible perimeter?

A first attempt might be to give variables to the side lengths, and somehow find the minimum perimeter using algebra and calculus tools. So there are four independent points — let’s parameterize them with four variables, and assume the side length of the square is 1:

Then we want to minimize this expression:

At this point, it isn’t clear how to proceed — there doesn’t seem to be any way to minimize this expression of four variables.

Proof by Net

We’ll have to try something different. It’s hard to make sense of anything when there are four independent variables. Instead, if we expand things out a bit, things start to become more manageable:

What we did was reflect the square three times, and each time the square is reflected, the inscribed quadrilateral goes with it. By taking only the relevant parts of the quadrilateral, we get the green path.

Now we might have a solution. If we had a different green path, can we reverse the steps and get the original quadrilateral back? Basically, the following requirements have to be met:

• The path has to cross all three of the internal lines BC, BA, and DA.
• The path’s position on the bottom-most line, DC must be the same when reflected onto the top-most line DC.

With these requirements in mind, the shortest green path that satisfies these requirements is a straight line connecting a point on the bottom left to its reflected point on the top right:

Our intuition at the start was well-founded.

Now notice that this isn’t the only possible shortest path. If we move the entire green line to the left or right, we get a different path of the same length!

For instance, the degenerate ‘quadrilateral’ formed by connecting two opposite corners has the same perimeter as the one we get by connecting the midpoints. Neat, huh?

Understanding Harmonic Conjugates (sort of)

For many people (for me at least), the Harmonic Conjugate is a difficult concept to understand. I didn’t really get it the first time I saw it, at Mathcamp. Let’s take the definition of the harmonic conjugate:

AB and CD are harmonic conjugates if this equation holds:

If you’re like me, you’re thinking along the lines of “But why? Why is this defined this way? Why would we spend so much time proving things about this weird concept? What’s the point, what’s the use?”

Even now, I can’t really give you an intuitive explanation of why this equality is so important. On the other hand, I could certainly come up with a few problems in which the concept of the harmonic conjugate turns to be useful.

Apollonius and Fleeing Ships

Apollonius’s problem was this: you are in control of a ship (point A on diagram), and you are in pursuit of another ship (point B). The other ship is fleeing in a straight line in some direction:

Your speed is (obviously) faster than the speed of the other ship: say they’re going at 30 km/h and you’re going at 50 km/h. Additionally, your ship is required to go in a straight line.

In which direction should you set off in order to intercept the fleeing ship?

Solution with Harmonic Conjugates

The first step of the solution is to construct harmonic conjugates CD so that their ratio is the ratio of your speed to the other ship’s speed (we’ll prove later that this is actually possible; assume we can do this for now):

Next, draw a circle with diameter CD:

There is a point where the ray from B (their ship) intersects this circle. Now go to this point immediately, in a straight line: the ship will be there.

The Proof

In order to prove that this works, we’ll need to take a step back and look at how we constructed the points C and D. The solution turns out to be evident directly from the construction of the harmonic conjugates.

Again, let’s assume our desired ratio is 5/3. Starting with the points A and B, the first step is constructing some point P so that:

This is fairly easy to do. Draw a circle of radius 5 around A, and draw a circle of radius 3 around B — the intersection P of these two circles forms the correct ratio. (if the circles don’t intersect, just scale everything up and try again)

Next, dropping the internal and external angle bisectors of the new triangle gives the harmonic conjugates C and D:

Why angle bisectors? From the angle bisector theorems (which I won’t prove here):

Combining the two proves that C and D are indeed harmonic conjugates to AB.

As a corollary, notice that because of angle bisecting, the angle CPD is always a right angle — hence, the locus of all points P forms a circle with diameter CD.

Returning to the ship problem, since each point P is defined as a point so that , it follows that when both ships travel to such a point P, they will meet at the same time.

Calculating the Plane Angles of Platonic Solids

What is the angle between two planes of a tetrahedron?

The angle between any two edges of the tetrahedron is , so it’s easy to (falsely) conclude that the angle between two faces must also be .

But this isn’t the case: the plane angle is defined as the angle between two lines on the two planes that are both perpendicular to the edge. None of the edges in a tetrahedron is perpendicular to any other, so the answer of is invalid.

We can try to compute the angle using regular Euclidean solid geometry, but things tend to get messy. A different way to approach the problem is by using vector geometry: using vector methods we can easily calculate the plane angle of the tetrahedron (as well as the icosahedron and the dodecahedron).

Assume symmetry. We represent three concurrent edges of the polyhedron as three vectors beginning at the same point: , , and ; let and be the angles between the vectors (by symmetry we’re assuming that the two alpha’s are equal):
(in case my poor drawing skills does not make this apparent, vectors a and b form one face of the polyhedron and c and b form another face)

For simplicity, let’s also say the lengths of each of the three vectors is 1.

We want to compute the angle between the plane formed by and , and the plane formed by and . Hence let and be the perpendiculars to and respectively each ending at the same point as their respective vectors:

For any two vectors, the dot product is defined with being the angle between the vectors. Given that the lengths of the vectors are all 1, we have:

Also, by vector addition:

Hence and .

We want to find the angle between and — call this angle . Then

The dot product of vectors x and y is simply:

Additionally . Hence the cosine of the angle is:

We can now use this newly derived formula to calculate plane angles! For example…

Tetrahedron

In the tetrahedron, both and are 60:
So  and .

Icosahedron

In the icosahedron, but :
The top ‘cap’ is a regular pentagon, which has a vertex angle of 108; each side of the pentagon constitutes a side of an equilateral triangle. Since , , and .

Dodecahedron

Computing angles for the dodecahedron works a bit differently from the tetrahedron and icosahedron. Instead of using existing edges as vectors, we construct an equilateral triangle by connecting three vertices:

So (since it’s part of a regular pentagon) and . Then and .

Varignon’s theorem proved in one line with vectors

Today I was reading some math book when the author mentions Varignon’s theorem, and gives a proof. The proof was not very long, but it was somewhat confusing. On Wikipedia, several more short proofs were given, but they were all more confusing than need be.

I remembered seeing the theorem proven using vector geometry before, but I couldn’t find the text (nor any other page / book that proves it this way) –

[image shamelessly taken from Wikipedia]

Varignon’s theorem states that in any quadrilateral, if we join the midpoints of the sides, then we get a parallelogram.

In the diagram, it suffices to prove that vector HG is equal to vector EF — vectors must have both the same orientation and length to be equal. This works since any method that proves HG = EF can also prove HE = GF. The proof goes as follows –

And we’re done. (the last step is due to symmetry of HG and EF)

Problem 10 of the 2011 Euclid Math Contest (remix)

I wrote the Euclid contest this year two days ago, on tuesday. There were 10 problems, and the tenth problem was a nice geometry problem. Three subproblems with a nice neat triangle on the right, with the subproblems getting progressively harder and harder. As I proceeded with the problem, I couldn’t help but imagine it as a Project Euler problem — instead of finding one such integer triangle, it would ask you to find the sum of the perimeters of all such integer triangles with perimeter less than 10^12 or some large number.

A modified problem

In the above diagram, and . Let , , , , and . Write an algorithm to find triangles satisfying these conditions where a, b, c are all integers.

Similar triangles

It is difficult to try to find integer triangles with such strange requirements as these. It seems that the line is completely unnecessary, but if we take it out, there doesn’t seem to be any way to relate the angle ratios to integer side lengths.

We can prove that is similar to . Being an exterior angle, , and also . Both of the triangles have an angle of measure and another angle of measure , thus they are similar.

From the similar triangles ratio

We can write d in terms of the three sides of the triangle:

Similarly, the side can be written as . Then we have the ratio

Solving for x allows us to express it in terms of the three sides of the triangle, again:

Constructing another similar triangle

Our goal here is to relate the lengths a, b, c with a simple equation, which then the problem turns into a number theory problem. Since we can write the lengths d and x in terms of a, b, and c, we can also relate any of a, b, c, d, x in an equation.

Again, there doesn’t seem to be a way to relate all of the variables together, in a way that any solution implies the original angle ratio required — unless we make a construction.

Here we extend AB to F, so that and is an isosceles triangle.

Again, since the exterior angle here is , both . Also with this construction, , and is also isosceles, hence .

With this construction, we can write the ratio

Perfect! Cross multiplying and then substituting in the original sides of the triangles gives

Simplifying this gives

Number theory magic

Now that we have an equation relating the three side lengths — it’s easy to verify that any three integers satisfying the triangle inequality gives a triangle with the required conditions — we can use number theory to try to generate integers that fit the equation.

If we expand the equation, we get

It makes sense to solve for c, which can be done using just the quadratic formula:

We need the discriminant D — the expression inside the square root — to be a perfect square, where we are allowed to have integer values for a and b. If we can get D to be a perfect square, then c will turn out to be a rational number. Then multiplying all three variables by a constant gives integer values for all three.

So we defined D:

Expanding this gives

Fortunately, this expression has an interesting factorization:

Or we can also write

We’ve simplified this problem to finding values where is a perfect square, that is:

This is just finding Pythagorean triples where one of the two sides are even! For instance, in the triple (3,4,5), we have a=3 and b=2. However, substituting a=3, b=2 into the quadratic formula gives c=5. This is almost a solution, only that the sides have to satisfy the triangle inequality (two sides have to add up to more than the third side).

The next Pythagorean triple (5,12,13) gives a=5 and b=6. Substituting this in gives c=11/9, which does satisfy the triangle inequality. Multiplying everything by 9 gives a=45, b=54, c=11 as the smallest working combination.

With this method, it is possible to quickly find arbitrarily many such triples, using Pythagorean triples as a starting point (which can be generated quickly with known methods).

Rotating a Hyperbola

The equation of a hyperbola gives it some very interesting properties under a series of stretches. Consider the hyperbola :

If we stretch the hyperbola horizontally by a factor of 2 about the x axis — replacing in the equation with , we have the equation .

Now, if we compress the resulting hyperbola vertically by a factor of 2 — replacing in the equation with and simplifying, we get the original hyperbola, . So by a horizontal stretch followed by a vertical compression, we get the original hyperbola.

However, consider a point on the hyperbola , say . Stretching the hyperbola horizontally by a factor of 2, this point gets transformed to . Even though each point has a different image point, however, the graph as a whole remains identical.

More generally, if we stretch by a factor of instead of 2, where is any positive number, we still obtain the original hyperbola after the two stretches. The point gets transformed to the point . We call this combinations of transformations a hyperbolic rotation. Intuitively, the points on the hyperbola are ‘rotated’ downwards.

There are some interesting things about the hyperbolic rotation. For one, by choosing the right for the transformation, it is possible to transform any point on the hyperbola into any other point on the same arm of the hyperbola (by allowing to be smaller than 1). Also, since the hyperbolic rotation is composed of two simple stretches, parallel lines as well as ratios of a line segments are preserved. The area of a figure is also preserved since one stretch reduces the area, while the other multiplies the area by the same amount.

With these facts, we can prove things about the hyperbola itself:

Theorem: Let P be a point on the hyperbola . Let be the line passing through P tangent to the hyperbola. Let X be the intersection of and the x-axis, and Y be the intersection of with the y-axis. Prove that P is the midpoint of XY.

The proof is very simple with the hyperbolic rotation. Consider when . The theorem obviously holds here, because of symmetry. But by applying a hyperbolic rotation, we can transform the point into any other point on the hyperbola. Since ratios between line segments don’t change with a hyperbolic rotation, P is always the midpoint, and we are done.

Calculus magic

You have probably seen some version of this curve at some time or another:

The blue lines are drawn by putting evenly spaced points in the interval for both the horizontal and vertical axis and connecting them.

As approaches infinity, the question is, what is the area occupied by the blue lines? The area is clearly finite since the entire shape is bounded by a 1×1 square, but how much of that square is occupied?

Intuitively it might seem that the blue curve is a circle, however that intuition is wrong. In the above diagram, the red curve is the actual circle, and clearly the blue curve is outside the circle. But fortunately the area can be found with some calculus.

If we choose some point on the x-axis, then obviously that point is connected to point on the y-axis. The slope of this line is and the equation of this line is . More conveniently, we can this expression as , intuitively the line formed when we choose as our x-axis point.

Let be the bounding curve for the blue lines. To find , we want to choose the value , given , so that is maximized. To do this, we take the partial derivative :

The optimal value of is the one for which , so solving for we get . Then we have

Integrating from 0 to 1 gives us the answer:

Notes on Mercator’s Projection

The following work is my math essay assignment (an essay on any topic related to the Math 30 curriculum). I found it interesting enough that I’ll post it here on my math blog. Yay.

Introduction

How do you make a world map? We find world maps everywhere — but the making of world maps is more complicated than most people realize. The earth is round, and a map is flat. How should one represent the surface of the earth, a sphere, on a flat sheet of paper?

The problem of displaying the surface of a sphere on a map has concerned cartographers, or map makers, for centuries. Over the years, cartographers have come up with a multitude of map projections, or systematic methods of transferring the surface of the earth onto flat paper.

There are many types of map projections, but one of the most common types is the cylindrical projection. Imagine a semi-transparent, hollow globe with a light inside of it. Now wrap a piece of blank map paper around it, like a cylinder:

Let the earth be represented as a unit sphere. The longitude (east-west location) shall be denoted as , and the latitude (north-south location) as . is the center of the earth, and a point on the equator.

Now a light source is placed at point ; a point on the surface of the globe will be projected in a straight line to point on the cylinder such that points , , and are collinear. Using some trigonometry, we have ; since , the length of is .

After repeating this procedure for every point on the sphere, the cylinder is unrolled into a flat sheet. The result is the following map:

Here, if a point on the earth is known (we know its latitude and longitude), we know its position on the map. The position coordinates of the map, can be described in terms of and in the following equations:

This projection has several major drawbacks: most notably, as the latitude increases toward , the slope of becomes increasingly steep, and the -coordinates of the point of the map increases towards infinity. At high latitudes, the map is very distorted.

The map is not equal-area — in an equal-area map, shapes of equal area on the sphere have equal area on the projection. Looking at the map we created, this is clearly not the case. Neither is the map conformal, or shape preserving. In a conformal map, the angle between any two lines must be the same as their counterparts on the map.

Mercator’s Projection

One of the most well known cylindrical map projections is Mercator’s projection. Invented by Belgian cartographer Gerardus Mercator in 1569, Mercator addresses some of the problems in older cylindrical projections. Unlike our projection, Mercator’s projection is conformal, and it has been used for navigation for centuries:

However, the equations for Mercator’s projection are rather complicated. The position on Mercator’s projection, when given the latitude and longitude, is:

This is a strange and seemingly random formula. But when we use this formula, we get a conformal map. How does this work?

Deriving Mercator’s Projection

Before deriving the formula for Mercator’s Projection, we take a step back and ask ourselves — what was the projection for in the first place?

Unlike most maps today which are used to decorate walls, Mercator’s map was more than that — it was used for nautical navigation. Suppose you were to travel by ship from Vancouver to Honolulu, equipped with nothing more than a compass. Which direction should you travel?

Now if you have a copy of Mercator’s map, this becomes simple. You draw a straight line from Vancouver to Honolulu, and with a protractor you find that the direction is south to the parallel. With a compass, it would be easy for you to maintain this constant heading for the entire duration of the trip.

(Footnote — The straight line on Mercator’s map is not straight on the Earth — it is actually a curved line, called a rhumb line. That is, after taking an initial bearing, one proceeds along the same bearing (relative to the north pole) for the duration of the trip. In the example, you would have to adjust your heading intermittently so that you are always heading at the direction south to the parallel.)

This technique is possible because it relies on one crucial property of Mercator’s map: that it is conformal. It would fail if we tried the same technique on the cylindrical projection given in the introduction.

However, notice that on Mercator’s map, all parallels have the same length: a parallel is simply a line across the width of the entire map. But on Earth, a sphere, parallels are longer when they are closer to the equator and shorter when near the poles:

Let the Earth be a unit sphere with center . is a point on the surface of the Earth and is its corresponding point on the equator. Let be a point such that is parallel to and is perpendicular to . As , the circumference of the Earth at the equator is . But the circumference at latitude is only .

It is easy to see why: as in the diagram, it follows that angles and are equal. Next, , therefore . The result follows.

On the map, the length of the equator and the parallel at latitude are equal. But on Earth, the parallel at latitude is smaller than the equator by a factor of ; in the projection any line at latitude ends up being stretched horizontally by a factor of or .

Finally, in order to satisfy the property of conformality, any part of the map stretched horizontally by some factor, say , must be stretched vertically by the same factor . Doing so preserves the angles of that part of the map:

Imagine a very small piece of land at latitude . When it is represented on the map, it is stretched horizontally by . Hence to preserve conformity, it must also be stretched vertically by .

But this can only work if the piece of land has no area. But as all pieces of land, no matter how small, has some area, the latitude of the piece of land is not , but goes from to for some small number . Similarly, the space this piece of land occupies on the map is not just , but goes from for some small number . Now as approaches 0, the ratio of to approaches :

Or equivalently,

Solving for and integrating both sides gives:

If we let be the denominator, , it turns out that the numerator is the derivative of , or . Thus we can substitute into the integral:

We defined the map so that the equator lies on the x-axis, therefore when , . Substituting and solving for , we find that . Also, since for , we can remove the absolute value brackets. Therefore,

This is what is desired.

An Alternative Map Projection

Mercator’s Projection is merely a conformal map, not a perfect map. It has several drawbacks — the scale varies greatly from place to place: as the latitude increases, the map gets stretched more and more. This severely distorts larger figures, especially areas near the poles. At a latitude greater than about , the map is practically unusable.

The Gall-Peters Projection is an alternative cylindrical projection:

When a map is used in navigation, it is important that angles are preserved in the map. But when a map is used to present statistical data, it’s often more important to preserve areas in the map: otherwise one may be deceived into thinking one country is bigger than another when they are actually the same size.

The Gall-Peters projection does this, trading conformity for equal-area — two countries of the same area will have the same area on the map, no matter where they are on the map. This map projection is produced by the following equations:

A discussion on why this set of equations produces an equal-area map is beyond the scope of this essay.

Conclusion

It is a tricky matter to accurately depict the surface of the earth on a flat sheet of paper. The method of cylindrical projection produces several considerably different maps.

The first map we considered uses only simple trigonometry, but the map produced is not very useful. The second map we investigate, Mercator’s Projection, preserves angles on the map. The last map we consider, the Gall-Peters projection, preserves areas on the map.

Each of these map projections are both similar and very different. They are similar in that they are all produced in some way by wrapping a cylinder around a spherical Earth. But the mathematics needed to produce them differ.

Fermat points and parameterizing the 120 degree integer triangle (Project Euler 143)

Problem 143 of Project Euler is an interesting geometry problem, and one notorious for having one of the lowest solvers to time released ratios of all the problems on the site: under 600 solvers in over 3 years as time of writing.

So what is the problem?

The problem is quite simple. We have a triangle , with all angles smaller than . Point is inside the triangle such that is at its minimum ( can always be determined uniquely).

We define a Torricelli triangle as one where , , , , , are all integers.

Generate all Torricelli triangles with .

The Fermat point

In triangle , with being the point such that is minimized, is called the Fermat point of triangle .

We can use Ptolemy’s theorem to locate the Fermat point in the triangle.

Let be a point such that is equilateral, then construct the circumcircle of :

We now show that must be at the intersection of and the circumcircle of .

Assume that does not lie on the circumcircle of . Then, according to Ptolemy’s theorem:

Equality only occurs when is cyclic. Since is equilateral, we have so we can divide it out to get

Adding to both sides:

Now if did lie on the circumcircle of , we would have , so must lie on the circumcircle of in order to be optimal.

So we know is on the circle. But now we assume that is not on . As we had before, .

Next if is not on , then , or by substitution,

Of course if were actually on , then , and the sum would be optimal.

This proves the optimal place for to be on the intersection of and the circumcircle of .

If we repeat this for the other three sides, we notice that is the intersection of the three circumcircles, and also the intersection of , , and :

Further, we see that quadrilaterals , , and are all cyclic. As , , and similarly and .

Designing an algorithm

We have enough information now to start working on an algorithm. Let us come back to the previous diagram:

So knowing that the central angles are all , we can apply the cosine law ():

A similar formula applies to sides and . We call a pair if is a perfect square.

We have found a Torricelli triangle if for three integers , , , all of , , are all pairs.

This leaves us with an outline of an algorithm:

1. Generate all pairs under the limit (with and being a square)
2. Sort the list of pairs and index them (to be explained in step 3)
3. For each pair in the list, search through the list to check if there exists some where is a pair and is a pair. We index the list to drastically improve searching time.
4. If an triple is found and , then mark as found. This is easily implemented as a bit array of size 120000 with all bits initially 0 and flipped to 1 when a sum is encountered.
5. Add up the sums and output

Step 1 takes up most of the time in this algorithm. To do this by ‘brute force’, we need an outer loop from 1 up to 120000, then another inner loop again from 1 to 120000 (actually a bit less), which is essentially operations, which is in the billions. An implementation of this algorithm would take about 3 minutes.

Parameterizing the a^2 + ab + b^2 = c^2 equation

It is well known that all primitive triples of the form can be generated by , , where and are coprime and are of opposite parity.

Can a similar parameterization be found for the equation equation?

Letting , and and dividing by , we get

which is the equation of an ellipse.

Originally we wanted to find integer solutions to the equation . This is equivalent to finding all rational points on the ellipse :

It is easy to see why: if and are rational points such that then and where , , are positive integers, and we arrive back at the original form of the equation.

Also in order for , , to be positive, we will only consider points on the ellipse in the first quadrant, with both and being positive. We do this by choosing a point on the ellipse, and from there combine the equations of the ellipse and the line.

Let us choose the point to be the first point of the line. Then the equation of the line is , where is valid only if it is positive and greater than 1.

Substituting into the equation of the ellipse:

Simplifying:

Now evaluating for :

Notice now that for and to be rational, we just need the slope to be rational. So we can write as where and are positive integers and .

Simplifying the expressions for and :

Since we defined and , we have a valid triple if , , and .

If and are not coprime, then neither will , , and be coprime as they will share a common divisor. So in order for to be primitive we will need and to be coprime.

Still this isn’t quite sufficient. Consider what would happen if . We have which can be written as which is divisible by 3. Next which is divisible by 3. Finally, which can be written as , again divisible by 3! So if , then the resulting triple is not primitive!

But it turns out that we can ignore all cases where . Given a triple , with , it is possible to find an identical parameterization giving .

As , we have and also . Thus we can have and where and are positive integers:

Multiplying by 3, , and . Combining the two and substituting we get

If we substitute and for and in the triple , we get:

So this proves that when , the case is nonprimitive and has already been covered by a smaller (possibly primitive) case. We are done.

Implementation in C++

Now using the new parameterization scheme, we can generate all pairs in about operations instead of operations. So we can loop up to , which cuts down the time down to about 2 seconds:

#include <cstdio>
#include <vector>
#include <algorithm>
#include <utility>
using namespace std;

// Limit and square root of limit
#define L 120000
#define LSQ 347
typedef long long int int64;

int64 gcd(int64 a, int64 b)
{
  int64 c = a % b;
  while(c != 0)
  {
    a = b;
    b = c;
    c = a % b;
  }
  return b;
}

int main()
{
  // Store pairs in here
  vector< pair<int64,int64> > pairs;

  // Use the parameterization
  for(int64 u=1; u<LSQ; u++){
    for(int64 v=1; v<u; v++){
      if(gcd(u,v) != 1) continue;
      if((u-v) % 3 == 0) continue;
      int64 a = 2*u*v + v*v;
      int64 b = u*u - v*v;
      if(a+b > L) break;
      // From coprime pairs make composite pairs
      for(int k=1; k*(a+b)<L; k++){
        pairs.push_back(make_pair(k*a,k*b));
        pairs.push_back(make_pair(k*b,k*a));
      }
    }
  }

  // Sort pairs list
  sort(pairs.begin(),pairs.end());

  // Create index
  int index[L];
  for(int i=0; i<L; i++) index[i] = -1;
  for(int i=0; i<pairs.size(); i++){
    if(index[pairs[i].first] == -1)
      index[pairs[i].first] = i;
  }

  // Which sums have been reached?
  bool sums[L];
  for(int i=0; i<L; i++) sums[i] = false;

  // Iterate through all pairs
  for(int i=0; i<pairs.size(); i++){
    int64 a = pairs[i].first;
    int64 b = pairs[i].second;

    // Construct vectors for indices
    vector<int64> va;
    vector<int64> vb;

    // Fetch indices
    int ia = index[a];
    int ib = index[b];

    while(ia<pairs.size()){
      pair<int64,int64> next = pairs[ia];
      if(next.first != a) break;
      va.push_back(next.second);
      ia++;
    }

    while(ib<pairs.size()){
      pair<int64,int64> next = pairs[ib];
      if(next.first != b) break;
      vb.push_back(next.second);
      ib++;
    }

    // Find common objects between va and vb
    for(int ja=0; ja<va.size(); ja++){
      for(int jb=0; jb<vb.size(); jb++){
        if(va[ja]==vb[jb]){
          // Potential c found
          int64 c = va[ja];
          if(a+b+c<L) sums[a+b+c] = true;
        }
      }
    }
  }

  // Tally up sums
  int64 s = 0;
  for(int i=0; i<L; i++)
    if(sums[i]) s+=i;

  printf("%lld\n",s);
  return 0;
}

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 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 radians, the entire sphere is denoted by sr. Likewise, one hemisphere would be 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 , , are vectors of three vertices of a tetrahedron (in relationship to the fourth), then the subtended solid angle is given by:

Here , , are magnitudes of the vectors, the dot represents a dot product, and is the determinant of the scalar product of the vectors.

The magnitude of the 3-D euclidean distance:

The dot product of two vectors is:

Finally,

Calculating the probability

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

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 with or about 1.9921.

We have ; now we let . To facilitate vector manipulations, we will now transpose the cube such that the center of gravity is at :

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 of the dice’s mass is removed, the resulting difference is slightly less than 1%.