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, \angle CAK = 2 \angle KAB and \angle ABC = 2 \angle KAB. Let a = CB, b = AC, c = AB, d = AK, and x = KB. 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 AK 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 \triangle CAK is similar to \triangle ABC. Being an exterior angle, CKA = 3 \theta, and also \angle CAB = 3 \theta. Both of the triangles have an angle of measure 2 \theta and another angle of measure 3 \theta, thus they are similar.

From the similar triangles ratio

\frac{b}{d} = \frac{a}{c}

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

d = \frac{bc}{a}

Similarly, the side CK can be written as a-x. Then we have the ratio

\frac{a}{b} = \frac{b}{a-x}

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

x = \frac{a^2 - b^2}{a}

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 KB = BF and \triangle KBF is an isosceles triangle.

Again, since the exterior angle here is 2 \theta, both \angle BKF = \angle BFK = \theta. Also with this construction, \triangle AKF \sim \triangle KBF, and is also isosceles, hence d = AK = KF.

With this construction, we can write the ratio

\frac{x}{d} = \frac{d}{c+x}

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

(\frac{a^2-b^2}{a})(c+\frac{a^2-b^2}{a}) = (\frac{bc}{a})^2

Simplifying this gives

(a^2 - b^2) (a^2 - b^2 + ac) = b^2 c^2

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

a^4+a^3 c-2 a^2 b^2-a b^2 c+b^4-b^2 c^2 = 0

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

c = \frac{a^3 - ab^2 \pm \sqrt{(ab^2-a^3)^2 + 4b^2(b^4 + a^4 - 2a^2 b^2)}}{2b^2}

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:

D = (ab^2-a^3)^2 + 4b^2(b^4 + a^4 - 2a^2 b^2)

Expanding this gives

D = a^6 - 7a^2 b^4 + 2 a^4 b^2 + 4b^6

Fortunately, this expression has an interesting factorization:

D = (a^2+4b^2) (a+b)^2 (a-b)^2

Or we can also write

\sqrt{D} = (a+b) (a-b) \sqrt{a^2 + 4b^2}

We’ve simplified this problem to finding values where a^2 + 4b^2 is a perfect square, that is:

a^2 + (2b)^2 = k^2

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).

Notes on Catalan’s Conjecture

Alright so I’ve been a bit busy lately with school and all the homework and such, leaving me with much less time to study and write about mathematics than during the summer. But after a three week break I’m now blogging again. Yay. Anyways.

Catalan’s conjecture is an interesting conjecture in number theory. In this equation

x^a - y^b = 1

The only solution with x, y, a, b > 1 is 3^2 - 2^3 = 1. That is to say, there are no two other powers with a difference of 1.

This was conjectured in 1844, but only recently proven in 2002 by Mihăilescu, so it’s known also as Mihăilescu’s theorem. Mihăilescu’s proof of the theorem uses some very sophisticated techniques, and is beyond the scope of this blog post.

It turns out that Catalan’s conjecture is useful for solving some interesting number theory problems. We can sometimes reduce a problem down to a very specific instance of Catalan’s conjecture, then apply the proof, essentially overkilling the problem. Here are some examples.

Problem 1

Find all triples of positive integers (a,b,c) where

a^2 + 2^{b+1} = 3^c


We look at powers modulo 4. It can be seen that 4 | 2^{b+1} when b is positive, thus 2^{b+1} \equiv 0 \mod 4 and a^2 \equiv 3^c \mod 4.

Next note that 3 \equiv -1 \mod 4 so 3^c \equiv (-1)^c \mod 4. As 3^c is always odd and 2^{b+1} is always even, a^2 must be odd hence a must be odd. Taking the modulo 4, a^2 \equiv 1 \mod 4 (if a \equiv 1 or a \equiv 3). Thus (-1)^c \equiv 1 \mod 4, implying that c is even.

Then c can be written as 2m for some m. The equation is then a^2 + 2^{b+1} = 3^{2m}, and with differences of squares,

2^{b+1} = 3^{2m} - a^2 = (3^m+a) (3^m-a) .

The product of 3^m+a and 3^m-a is a power of 2, so both must themselves be powers of 2. Then there exists x and y where x,y \in \mathbb{N}_0 and x<y such that

3^m - a = 2^x

3^m + a = 2^y

x+y = b+1

Adding the first two yields

2 \cdot 3^m = 2^x + 2^y

Putting x=0 causes a parity error. If x \geq 2 then the RHS is divisible by 4 while the LHS clearly isn’t. This makes x=1, giving

3^m = 1 + 2^{y-1} .

This is a case of Catalan’s conjecture. The only solutions to (m,y) are (1,2) and (2,4). Now we can just substitute the solutions for (m,y) to get solutions for (a,b,c). They are (1,2,2) and (7,4,4).

Problem 1a

Here is a very similar problem: find all integer solutions to

3^a + 4^b = c^2

The solution is left to the reader.

Hint: rewrite as a difference of squares and subtract

Problem 2

Find all triples (x,y,n) of non-negative integers such that

x^n = y^4 + 4


Obviously n cannot be 0. We first look at n=1. Then, x = y^4 + 4, and any triple (y^4+4, y, 1) satisfies the equation.

Next we consider when n=2. Then x^2 = y^4 + 4, or as a difference of squares, (x+y^2) (x-y^2) = 4. The solution then is x=2 and y=0, or (2,0,2).

Now we try n>2. Suppose that y is even. Then in the equation x^n = y^4 + 4, 4 | y^4 and 4 | RHS so 4 | LHS, thus x is even. Then x^n = 2^n k^n for some k, and it is clear that n \leq 2, which is a contridiction.

So y must be odd. The RHS, y^4 + 4, can be written as y^4 + 4y^2 + 4 - 4y^2, or (y+2+2y)(y+2-2y). This gives us

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

We claim that for odd y, y^2 + 2y + 2 \perp y^2-2y+2. Let k be the GCD of y^2+2y+2 and y^2-2y+2 so that k | y^2+2y+2 and k | y^2-2y+2. Then k|4y (the difference). As y^2+2y+2 is odd, and k divides it, k must be odd so k|y. But y^2+2y+2 = y(y+2)+2 so k|y and k|y(y+2)+2, thus k=1.

Hence for some integers u and v where uv=x,

u^n = y^2 + 2y + 2

v^n = y^2 - 2y + 2

Alternatively, u^n = (y+1)^2 + 1 and v^n = (y-1)^2 + 1.

From Catalan’s conjecture, the equation x^a - y^2 = 1 with a>2 has no solutions other than 1^n-0^2=1. But we have the simultaneous equations u^n - (y+1)^2 = 1 and v^n - (y-1)^2 = 1. This is clearly impossible.

Thus the only solutions are (y^4+4,y,1) and (2,0,2).

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 ABC, with all angles smaller than 120^\circ. Point T is inside the triangle such that p+q+r is at its minimum (T can always be determined uniquely).

We define a Torricelli triangle as one where a, b, c, p, q, r are all integers.

Generate all Torricelli triangles with p + q + r \leq 120000.

The Fermat point

In triangle ABC, with P being the point such that PA+PB+PC is minimized, P is called the Fermat point of triangle ABC.

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

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

We now show that P must be at the intersection of AD and the circumcircle of \triangle BCD.

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

BP \times CD + PC \times BD > PD \times BC

Equality only occurs when PCDB is cyclic. Since BCD is equilateral, we have BC = CD = BD so we can divide it out to get

PB + PC > PD

Adding PA to both sides:

PA + PB + PC > PA + PD

Now if P did lie on the circumcircle of \triangle BCD, we would have PA + PB + PC = PA + PD, so P must lie on the circumcircle of \triangle BCD in order to be optimal.

So we know P is on the circle. But now we assume that P is not on AD. As we had before, PA + PB + PC = PA + PD.

Next if P is not on AD, then AP + PD > AD, or by substitution,

PA + PB + PC > AD

Of course if P were actually on AD, then PA + PB + PC = AD, and the sum PA+PB+PC would be optimal.

This proves the optimal place for P to be on the intersection of AD and the circumcircle of \triangle BCD.

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

Further, we see that quadrilaterals BPCD, APCE, and FAPB are all cyclic. As \angle BDC = 60^\circ, \angle BPC = 120^\circ, and similarly \angle APC = 120^\circ and \angle APB = 120^\circ.

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 120^\circ, we can apply the cosine law (\cos 120^\circ = -\frac{1}{2}):

a^2 = q^2 + r^2 - 2qr \cos 120

a^2 = q^2 + r^2 +qr

A similar formula applies to sides b and c. We call (x,y) a pair if x^2 + xy + y^2 is a perfect square.

We have found a Torricelli triangle if for three integers p, q, r, all of (p,q), (p,r), (q,r) are all pairs.

This leaves us with an outline of an algorithm:

  1. Generate all pairs (x,y) under the limit (with x+y < 120000 and x^2 + xy + y^2 being a square)
  2. Sort the list of pairs and index them (to be explained in step 3)
  3. For each pair (a,b) in the list, search through the list to check if there exists some c where (a,c) is a pair and (b,c) is a pair. We index the list to drastically improve searching time.
  4. If an (a,b,c) triple is found and a+b+c \leq 120000, then mark a+b+c 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 120000^2 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 a^2 + b^2 = c^2 can be generated by a = m^2 - n^2, b = 2mn, c = m^2 + n^2 where m and n are coprime and are of opposite parity.

Can a similar parameterization be found for the equation a^2 + ab + b^2 equation?

Letting x = \frac{a}{c}, and y = \frac{b}{c} and dividing by c^2, we get

x^2 + xy + y^2 = 1

which is the equation of an ellipse.

Originally we wanted to find integer solutions to the equation a^2 + ab + b^2 = c^2. This is equivalent to finding all rational points on the ellipse x^2 + xy + y^2 = 1:

It is easy to see why: if x and y are rational points such that x^2 + xy + y^2 = 1 then x = \frac{a}{c} and y = \frac{b}{c} where a, b, c are positive integers, and we arrive back at the original form of the equation.

Also in order for a, b, c to be positive, we will only consider points on the ellipse in the first quadrant, with both x and y 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 (0,-1) to be the first point of the line. Then the equation of the line is y = tx-1, where t is valid only if it is positive and greater than 1.

Substituting into the equation of the ellipse:

x^2 +tx^2 - x + t^2 x^2 - 2tx + 1 = 1


x (t^2+t+1) = 2t + 1

x = \frac{2t+1}{t^2 + t + 1}

Now evaluating for y:

y = t (\frac{2t+1}{t^2 + t + 1})-1 = \frac{t^2-1}{t^2+t+1}

Notice now that for x and y to be rational, we just need the slope t to be rational. So we can write t as \frac{m}{n} where m and n are positive integers and m>n.

Simplifying the expressions for x and y:

x = \frac{2 \frac{m}{n} + 1}{(\frac{m}{n})^2 + \frac{m}{n} + 1} = \frac{2mn+n^2}{m^2 + mn + n^2}

y = \frac{(\frac{m}{n})^2 - 1}{(\frac{m}{n})^2 + \frac{m}{n} + 1} = \frac{m^2 - n^2}{m^2 + mn + n^2}

Since we defined x = \frac{a}{c} and y = \frac{b}{c}, we have a valid triple if a = 2mn+n^2, b = m^2 - n^2, and c = m^2 + mn + n^2.

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

Still this isn’t quite sufficient. Consider what would happen if m \equiv n (\mod 3). We have a = 2mn+n^2 which can be written as n(n-m+3m) which is divisible by 3. Next b = (m+n) (m-n) which is divisible by 3. Finally, c = m^2 + mn + n^2 which can be written as (m-n)^2 + 3mn, again divisible by 3! So if m \equiv n \mod 3, then the resulting triple is not primitive!

But it turns out that we can ignore all cases where m \equiv n \mod 3. Given a triple (2mn+n^2, m^2-n^2, m^2+mn+n^2), with m \equiv n \mod 3, it is possible to find an identical parameterization giving (\frac{2mn+n^2}{3}, \frac{m^2-n^2}{3}, \frac{m^2 + mn + n^2}{3}).

As m \equiv n mod 3, we have m+2n \equiv 0 \mod 3 and also m-n \equiv 0 \mod 3. Thus we can have u and v where u and v are positive integers:

u = \frac{m+2n}{3}

v = \frac{m-n}{3}

Multiplying by 3, 3u = m+2n, and 3v = m-n. Combining the two and substituting we get

m = u+2v

n = u-v

If we substitute u and v for m and n in the triple (2mn+n^2, m^2-n^2, m^2+mn+n^2), we get:

2mn + n^2 = 3 (u^2-v^2)

m^2 - n^2 = 3(2uv + v^2)

m^2 + mn + n^2 = 3 (u^2 + uv + v^2)

So this proves that when m \equiv n \mod 3, 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 n operations instead of n^2 operations. So we can loop up to \sqrt{120000}, 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++){

  // Sort pairs list

  // 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];

      pair<int64,int64> next = pairs[ia];
      if(next.first != a) break;

      pair<int64,int64> next = pairs[ib];
      if(next.first != b) break;

    // Find common objects between va and vb
    for(int ja=0; ja<va.size(); ja++){
      for(int jb=0; jb<vb.size(); 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;

  return 0;

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.

On some number-theoretic properties of right triangles (Project Euler 218)

The two hundred and eighteenth problem of Project Euler is quite interesting, but different. It resembles more closely a mathematics olympiad problem than a programming challenge. Its answer is somewhat surprising, too.

Original problem

The problem can be stated as follows:

A right triangle is considered perfect if (1): it is primitive (GCD of the three sides is 1) and (2): its hypotenuse is a perfect square.

A perfect triangle is considered superperfect if its area is divisible by the perfect numbers 6 and 28.

How many perfect triangles are there with hypotenuse below 10^{16} that are not also superperfect?


It turns out that every perfect triangle is also superperfect, meaning that for any limit there are no perfect triangles that are not superperfect.

Looking on the forums, it seems that a large group of solvers counted the triangles for a smaller limit, like 10^{9} or 10^{12} , and getting 0, assumed the answer applied for 10^{16} .

In this article I will attempt to prove, mathematically, that it is indeed impossible for a perfect triangle not to be superperfect.


Let’s rephrase this problem a bit: Prove that if a primitive right triangle has a hypotenuse that is a perfect square then its area must be a multiple of 6 and 28.

If the area is a multiple of 6 and 28, then it is a multiple of \mathrm{LCM}(6,28) = 84 . If we let p, q, and c be the sides of the right triangle (with c as the hypotenuse), then the area is \frac{pq}{2} .

Since 84 | \frac{pq}{2} , it follows that 168 | pq . As c is a perfect square, we write c as r^2 and since \mathrm{GCD}(p,q,c)=1 , it also follows that \mathrm{GCD}(p,q,r)=1 . This is what we shall now prove.

Lemma 1

For positive integers p, q, r where p^2 + q^2 = r^4 and \mathrm{GCD}(p,q,r)=1 , then 168|pq .

From the Euclid theorem of Pythagorean triples, the sides of a primitive right triangle with sides a, b, and c can be represented as:

a = u^2 - v^2

b = 2uv

c = u^2 + v^2

.. where u>v and u and v are of opposite parity.

Thus by applying the theorem to our triple:

p = u^2 - v^2

q = 2uv

r^2 = u^2 + v^2

Notice here that the third equation here itself represents a pythagorean triple. We then apply the same formula again, this time using m and n for the integers:

u = m^2 - n^2

v = 2mn

r = m^2 + n^2

Substituting for p, q, r:

\begin{array}{rcl} p &=& u^2 - v^2 \\ &=& (m^2-n^2) - (2mn)^2 \\ &=& m^4 - 2m^2n^2 + n^4 - 4m^2n^2 \\&=& m^4 + n^4 - 6m^2 n^2 \\ \\ q &=& 2uv \\ &=& 2(m^2-n^2)(2mn) \\ &=& 4mn(m^2-n^2) \\ \\ r &=& m^2 + n^2 \end{array}


pq = 4mn(m^2-n^2)(m^4+n^4-6m^2n^2)

.. which we will prove to be divisible by 168. Since 168=8*3*7, the expression must be proved to be divisible by 8, by 3, and by 7.

Lemma 2

8 | 4mn(m^2-n^2)(m^4+n^4-6m^2n^2)

Proof of Lemma 2

The proof is almost trivial. In order for the triple to primitive, m and n must be of opposite parity, meaning mn is even.

Because 2|mn and 4 is already a factor of the expression, it follows that 8|pq .

Lemma 3

3 | mn(m^2-n^2)

Proof of Lemma 3

Rewrite the expression as

mn(m+n)(m-n) .

If 3|m or 3|n , then 3|mn .

If m \equiv n \mod 3 , then m-n \equiv 0 \mod 3 .

The only other scenario is when m \not \equiv n \mod 3 , then either m \equiv 1 \mod 3 and n \equiv 2 \mod 3 or vice versa. Either way, m+n \equiv 0 \mod 3 .

Lemma 4

If 7 \nmid m , then m^2 \equiv 1,2,4 \mod 7 .

Proof of Lemma 4

We construct a table. The first column of this table is m \mod 7 while the second column is m^2 \mod 7:

1 1
2 4
3 2
4 2
5 4
6 1

The only possible values are 1, 2, and 4.

Lemma 5

If 7 \nmid m^2 - n^2 , then either m^2 \equiv 2n^2 \mod 7 or n^2 \equiv 2m^2 \mod 7 .

Proof of Lemma 5

Because 7 \nmid m^2 - n^2 , m^2 \not\equiv n^2 \mod 7 . Then, not counting reflective cases, there are three cases we need to consider:

Case 1: m^2 \equiv 1, n^2 \equiv 2 . Then n^2 \equiv 2m^2 \mod 7 .

Case 2: m^2 \equiv 1, n^2 \equiv 4 . Then m^2 \equiv 2n^2 \mod 7 .

Case 3: m^2 \equiv 2, n^2 \equiv 4 . Then n^2 \equiv 2m^2 \mod 7 .

One of these (or their reflection) apply for whatever value of m and n.

Lemma 6

7 | m^4 + n^4 - 6m^2n^2

Proof of Lemma 6

Without loss of generality, rewrite the expression as a congruence mod 7, in terms of m. Since m^2 \equiv 2n^2 \mod 7,

\begin{array}{rcl} m^4 + n^4 - 6m^2n^2 &=& m^4 + (2m^2)^2 - 6m^2(2m^2) \\ &=& m^4 + 4m^4 - 12m^4 \\ &=& -7m^4 \end{array}

The result follows.


Socks in a Drawer

I’ve began studying some probability and statistics for fun, and I’m using a book called Fifty Challenging Problems in Probability by Mosteller:

Printing math textbooks is a great way to use up the school’s printing credits at the end of the year.

The first problem is about socks in a drawer. At first glance it looks really, deceptively simple. But when I start doing the problem I find that it’s a bit trickier than I expected.

Solutions to the problems are provided at the back of the book. For this problem, I find that the book’s solution is really ugly and unintuitive. Thus I’m going to give an alternate solution here.

Problem (my version)

A drawer contains red and black socks. I randomly take a sock out of the drawer. I then randomly take another sock out of the drawer.

The probability that both are red is exactly \frac{1}{2}. Determine how many socks of each color there are in the drawer.


Let r be the number of red socks in the drawer, and b be the number of black socks. Then the total number of socks in the drawer is r+b.

The probability that a sock you randomly take is red would be


Now the probability that the next sock you take would also be red is


This is because there is one less red sock in the drawer, and at the same time there is one less sock (of any color) in the drawer.

The probability that both socks will be red is the two probabilities multiplied:

P = \frac{1}{2} = \frac{r}{r+b} \cdot \frac{r-1}{r+b-1}

Now we’re solving the equation to find possible values of r and b. There are two variables, and one equation. However, both variables have to be positive integers. This type of equation is called a diophantine equation.

We begin by simplifying the equation into a quadratic.

\begin{array}{rcl} \frac{r(r-1)}{(r+b)(r+b-1)} &=& \frac{1}{2} \\(r+b)(r+b-1)&=& 2r(r-1) \\ r^2+2rb+b^2-r-b&=&2r^2-2r \\ r^2+(-2b-1)r+(b-b^2)&=&0 \end{array}

Now using the quadratic formula to solve for r:

r = \frac{(2b+1) + \sqrt{(-2b-1)^2-4(b-b^2)}}{2}

It’s not necessary to consider the negative root of the quadratic. Simplifying the quadratic further, we get

r = \frac{(2b+1)+ \sqrt{8b^2+1}}{2}

Here we can make a substitution. Let k = \sqrt{8b^2+1} . Then,

r = \frac{(2b+1)+ k}{2}

In this equation r will be an integer as long as b is an integer and k is an odd integer.

We are now finding integer solutions to the equation

k = \sqrt{8b^2+1}

..where k is also odd.

Simplifying, we get:

\begin{array}{rcl} k^2 &=& 8b^2+1 \\ k^2-8b^2&=&1 \end{array}

At this point we can drop the condition that k is odd, because if k is even then left hand side is even, but the right hand side is 1. As a result, k cannot be even if it’s a solution to this equation.

We’ve basically converted this problem into a form of Pell’s equation, which generally has this form:

x^2 - Dy^2 = 1

..where D is a constant and x and y are unknown integers. Our instance of Pell’s equation is the case where D=8.

Fortunately for us, the Pell’s equation is well studied. We can use existing results in number theory to help us with this problem.

By trial and error, it can be found that the smallest nontrivial solution, or the fundamental solution to the equation is (k,b) = (3,1). Once the fundamental solution to a Pell type equation has been found, all other solution pairs can be generated. And there are infinitely many such pairs.

Generating solution pairs can be done with this recurrence relation:

\begin{array}{rcl} k_{i+1} &=& k_1 k_i + D b_1 b_i \\ b_{i+1} &=& k_1 b_i + b_1 k_i \end{array}

It’s simple to write a computer program to generate these pairs. For instance, the first ten pairs of (k,b) are:


All that is left is to substitute values of b to get integer values of r. The first ten pairs of (r,b) are:


We are done; we can use the same recurrence relation to generate arbitrarily large sock combinations.

Random Math Problems (3)

School finished yesterday. After one exam next week (for math), I’m done grade 10. I’ll probably have more time to write blog posts in the summer.

Problem 1

(from Geometry Revisited)

Prove that in this triangle (a being m+n), that

a(p^2+mn) = b^2m + c^2n

This seems like a really weird thing to prove, but apparently it’s known as Stewart’s theorem and is useful in some places.


Recall the Law of Cosines:

\cos C = \frac{a^2+b^2-c^2}{2ab}

We can apply the Law of Cosines twice to \angle AXB and \angle AXC:

\cos \angle AXB = \frac{p^2+m^2-c^2}{2mp}

\cos \angle AXC = \frac{p^2+n^2-b^2}{2np}

Adding the two together,

\cos \angle AXB + \cos \angle AXC = \frac{p^2+m^2-c^2}{2mp} + \frac{p^2+n^2-b^2}{2np}

Notice that because \cos(\theta) = \cos(180-\theta), the sum of cosines of supplementary angles is 0. Therefore,

\begin{array}{rcl} 0 &=& \frac{p^2+m^2-c^2}{2mp} + \frac{p^2 + n^2 - b^2}{2np} \\ &=& n(p^2+m^2-c^2) + m(p^2 + n^2 - b^2) \\ &=& np^2 + nm^2 - nc^2 + mp^2 + mn^2 - mb^2 \end{array}

Moving the two negative terms to the left hand side and factoring,

\begin{array}{rcl} b^2m + c^2n &=& np^2 + nm^2 + mp^2 + mn^2 \\ &=& n(p^2) + m(p^2) + n(mn) + m(mn) \\ &=& (m+n)(p^2+mn) \end{array}

The result follows:

b^2m + c^2n = a(p^2+mn)

Problem 2

(from the 2010 UVM contest)

How many ways are there to arrange the letters of the word NOODLES, so that the consonants (NDLS) appear in alphabetical order (although not necessarily consecutively)?


The four consonants must appear in the order of D, L, N, S. The three vowels, OOE, must appear somewhere in the five spaces between the consonants:

Thus the answer to the problem is the same as the number of ways to pack OOE into five boxes, a box being able to hold an arbitrary amount of letters, and permutations of letters within a box being counted as distinct.

This is half of the number of ways to pack ABC into five boxes, since there are two O’s in OOE. We will now calculate how many ways there are to pack ABC into five boxes.

We can split up this problem into several cases:

Case 1: All three letters ABC are in the same box. For any box there are 3! = 6 ways to permute ABC, while there are five boxes. The total for this case is 6*5 = 30.

Case 2: The three letters all go in three different boxes. This is simply P(5,3) which is 60.

Case 3: Two of the letters go in a box, while the third letter goes in a different box. This is a little bit tricky, so we break it in a few parts.

First of all we choose the two boxes of the five we use. The number of ways to do this is \binom{5}{2} which is 10. Secondly, if the two boxes are labelled Box 1 and Box 2, we can either put two letters in Box 1 and one in Box 2, or the other way around, giving an additional two ways.

Finally there are 3! or 6 ways to permute the three letters. The final count for case 3 is 10 * 2 * 6 or 120.

The number of ways to pack ABC into five boxes is the sum of the three cases, or 210. The answer to the problem is half of that, 105.

Problem 3

(again from the same UVM contest)

Point P is inside rectangle ABCD. If the distance from P to three corners of the rectangle are 5, 10, and 14, find the distance from P to the forth corner.


We first draw perpendicular lines from P to the four sides of the rectangle and label them a, b, c, and d, like this:

From pythagorean, we find that

b^2 + c^2 = 5^2

b^2 + d^2 = 10^2

a^2 + d^2 = 14^2

As we are looking for the length of DP, we can manipulate the three equations to get a value for a^2 + c^2:

(b^2 + c^2) + (a^2 + d^2)-(b^2+d^2) = a^2+c^2

Finally we can now substitute in real values to find a^2 + c^2:

a^2+c^2 = 5^2 + 14^2 - 10^2 = 121

The line DP can be found by taking a square root, it is 11.

General formula

When three of the four diagonals of a rectangle are known, the forth one can be calculated. Suppose we know m, n, and o. The formula for x is given by:

x = \sqrt{m^2 + n^2 - o^2}

Problem 4

This problem is from the Number Theory chapter of The Art and Craft of Problem Solving.

Prove that

\frac{n^3 + 2n}{n^4 + 3n^2 + 1}

is in lowest terms for all positive integers n.


Another way of saying this problem is that

\textrm{GCD}(n^3 + 2n,n^4 + 3n^2 + 1) = 1

for all positive integers n.

Then there does not exist any integer c, where c > 1, c | n^3+2n and c | n^4+3n^2+1. Supposing c exists, and it divides n^3+2n, we can prove by contradiction that c cannot simultaneously divide n^4+3n^2+1.

If c divides n^3+2n or n(n^2+2), then either c | n or c | n^2+2.

We can rewrite n^4 + 3n^2 + 1:

n^3 + 3n^2 + 1 = n^2(n^2+3) + 1

This way it can be seen that n \perp n^2(n^2+3) , as two consecutive integers are relatively prime.

To show that n^2+2 is also relatively prime, we rewrite it this way:

n^3 + 3n^2 + 1 = (n^2+1)(n^2+2)-1

This shows that, n^2+2 \perp (n^2+1)(n^2+2)-1.

As both n and n^2+2 are relatively prime to n^4 + 3n^2 + 1, it follows that c \perp n^4 + 3n^2 + 1, which is what needs to be shown.

Digit sum divisibility rule proofs in bases other than 10

The divisibility test for 3 and 9 is fairly well known. To find out of any given integer is divisible by 9, add up the digits of the number; if the result is divisible by 9 then the number is divisible by 9.

What’s a bit less well known is that in any base b, the same trick applies for b-1 and all of its divisors.

For example, if we are doing arithmetic in hexadecimal (base 16), the divisibility rule we use for base 10 applies not for 3 and 9, but instead for 3, 5, and 15.

Suppose we were to confirm that \textrm{831f3f}_{16} is divisible by \textrm{f}_{16} . Adding up the digits:

\begin{array}{l} 8_{16} + 3_{16} + 1_{16} + \textrm{f}_{16} + 3_{16} + \textrm{f}_{16} \\ = \textrm{2d}_{16} \end{array}

And since \textrm{2d}_{16} is divisible by \textrm{f}_{16} , so is the original number.


Let n be an integer, b be the base, and m be the number of digits in n. Then we can represent n:

n = d_{m-1} b^{m-1} + d_{m-2} b^{m-2} + \cdots + d_1 b + d_0

We can rearrange this:

\begin{array}{rl} n =& [d_{m-1} (b^{m-1}-1) + \cdots + d_1 (b-1)] \\ &+ (d_{m-1} + d_{m-2} + \cdots + d_1 + d_0) \end{array}

Each of the expressions b-1 , b^2-1 , up to b^{m-1}-1 are divisible by b-1 . I’ve explored the proof for this in an earlier blog post.

Thus the entire first group of the above expression is divisible by b-1. What remains is the second group, which happens to be the sum of digits in n.

The sum of digits in n is divisible by b-1 if and only if n is divisible by b-1, which is what we wanted to prove.

For a factor of b-1, this also works as the entire first expression is divisible by any factor of b-1.

The Sieve of Sundaram

The Sieve of Eratosthenes is probably the best known algorithm for generating primes. Together with wheel factorization and other optimization options, it can generate primes very quickly.

But a lesser well known algorithm for sieving primes is the Sieve of Sundaram. This algorithm was discovered in 1934 by Sundaram; like the sieve of Eratosthenes it finds all prime numbers up to a certain integer.

The algorithm

A simplified version of the algorithm, using N as the limit to which we want to find primes to:

m =: Floor of N/2
L =: List of numbers from 1 to m
For every solution (i,j) to i + j + 2ij < m:
    Remove i + j + 2ij from L

For each k remaining in L:
    2k + 1 is prime.

In practice we can find solutions to i + j + 2ij < m by using two nested for loops:

For i in 0 to m:
    For j in i to m:
        L[i + j + 2ij] =: False

Here i is always less than j, because the two are interchangeable and filtering it twice would be a waste.

We don’t actually need to loop j from 0 to m. From the inequality i + j + 2ij < m, we can solve for j: j < \frac{m-i}{2i+1}. The new algorithm:

m =: Floor of N/2
L =: Boolean array of length m
Fill L with true

For i in 0 to m:
    For j in i to (m-i)/(2i+1):
        L[i + j + 2ij] =: False

For each k remaining in L:
    2k + 1 is prime.

Why this algorithm works

In the algorithm, 2k+1 is prime where k can be written as i + j + 2ij where i and j are integers. We can rewrite this:

\begin{array}{l} 2(i+j+2ij) + 1 \\ = 2i + 2j + 4ij + 1 \\ = (2i+1) (2j+1) \end{array}

Both 2i+1 and 2j+1 are odd numbers, and any number that can be written as the product of two odd numbers are composite.

Of the odd numbers, those that cannot be written as the product of two odd numbers are prime. We’ve filtered everything that can be written as (2i+1)(2j+1) so we are left with the odd prime numbers.

This algorithm only gets the odd prime numbers, but fortunately there is only one even prime number, 2.

Benchmarks with the Sieve of Eratosthenes

Here’s an implementation of the Sieve of Sundaram:

#include <stdio.h>
#include <stdlib.h>

typedef unsigned long long ll;
int main() {
 ll n = 100000000LL;
 ll m = n/2;
 char *sieve = malloc(m);

 for(ll i=0; i<m; i++)
  sieve[i] = 1;

 for(ll i=1; i<m; i++)
  for(ll j=i; j<=(m-i)/(2*i+1); j++)
   sieve[i+j+2*i*j] = 0;

 ll s=1;
 for(ll i=1; i<m; i++)
  if(sieve[i]) s++;

 printf("%llu", s);
 return 0;

This code counts the number of primes below 100 million, which should be 5761455. The above code runs in 9.725 seconds.

Here’s an alternative, an implementation of the more standard Sieve of Eratosthenes:

#include <stdio.h>
#include <stdlib.h>

typedef unsigned long long ll;
int main(){
 ll lim = 100000000LL;
 char *sieve = malloc(lim);

 for(int i=0; i<lim; i++)
  sieve[i] = 1;

 int s=0;
 for(int i=2; i<lim; i++){
   for(int j=2; j<=lim/i; j++)
    sieve[i*j] = 0;

 printf("%d", s);
 return 0;

I was surprised to find that the Sieve of Eratosthenes actually ran faster. It completed in 7.289 seconds.

I expected the Sieve of Sundaram to be faster because according to Wikipedia this algorithm uses O(n \log(n)) operations, while the Sieve of Eratosthenes uses O(n \log(n) \log \log(n)).

Challenge of the Week 02/23/2010

The Challenge of the Week is a weekly math challenge site run by the University of Washington. It has a collection of interesting problems. In each week, people are to submit solutions. At the end of the week, the solutions are revealed.

I think I’ve managed to solve this week’s challenge. I’ll reproduce the problem here:

Each of n numbers, x_1, x_2, \cdots x_n is selected from the set \{-1, 1\}. Prove that if

x_1 x_2 x_3 x_4 + x_2 x_3 x_4 x_5 + \cdots x_n x_1 x_2 x_3 = 0

then n is a multiple of 4.

My solution

I immediately noticed a few things about this problem.

The list of x_n values wraps around, so there really isn’t a start and end to the list. I thought of it as a circle:

A valid group of four is any element and the three elements after it. So for a circle of n elements, there are n groups of four:

So here x_5, circled in red, is the first element of the group. The next three elements circled in blue are the rest of the group. So this group would be x_5 x_6 x_1 x_2.

Because each of x_n is either 1 or -1, the product of the group of four must be either 1 or -1.

For it to equal zero, the number of groups being 1 (positive) and the number of groups being -1 (negative) must be the same. As well as an individual element being either positive or negative, a group of four is also either positive or negative.

This immediately eliminates all circles of size n where n is odd, because if there are n odd groups, there can be no way to split the groups up into equal amounts of positive and negative groups.

Now here’s my approach of proving that n cannot be of the form 4k + 2, or in other words when n is even but not divisible by 4.

Let’s consider the circle as initially all positive:

In the above image, a plus sign beside an element refers to the sign of the group of four starting on that element.

The 16 in the middle refers to the value of the whole circle, or the whole expression. We need that to be zero.

Now let’s make one of the elements negative:

Doing this makes all four groups containing that particular element negative. So we have four less positives, and four more negatives. Instead of 16, the total sum is 8.

Do it again and the sum is 0, which is exactly what we need:

Reading from the topmost clockwise, this list would be [1,1,1,-1,1,1,1,1,1,1,-1,1,1,1,1,1].

Changing one element from positive to negative one does not always negate the four groups containing it. Instead, if flips the signs of the four groups.

This is simply because -1 * -1 = 1.

I’ll make an example from the previous diagram:

Here, negating one element actually increased the overall value! It made several negative groups back into positives.

There are four general cases:

What it means is, by changing an element into -1, we flip its four containing groups. Doing so will take away 8 from the total, take away 4 from the total, or add 4 to the total.

Note that there is no case where (----) becomes (++++) because that would mean the head element is already negative, and cannot be made negative again.

Remember that we start with the total being n, and we want the total to be 0. Each step, you can change the total by -8, -4, or 4, each of which is a multiple of 4.

If the total is initially a multiple of 4, it will still be a multiple of 4 after any iteration.

If it’s not already a multiple of 4, nothing can make it a multiple of 4, and there is no way to make it zero by negating elements.