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

August 25, 2010

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

Simplifying:

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++){
        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;
}

Project Euler 299: Three similar triangles

July 29, 2010

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

Simplifying:

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 ,

Or,

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)

June 20, 2010

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?

Answer

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.

Proof

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}

Then,

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.

Q.E.D.


Project Euler 285

April 3, 2010

Problem 285 of Project Euler was released yesterday; it requires some geometry and probability sense.

The problem goes like this:

Albert chooses two real numbers a and b randomly between 0 and 1, and a positive integer k.

He plays a game with the following formula:

k' = \sqrt{(ka+1)^2 + (kb+1)^2}

The result k’ is rounded to the nearest integer. If k’ rounded is equal to k, then Albert gets k points. Otherwise, he gets nothing.

Find the expected value of the total score if Albert plays this game with k = 1, 2, \cdots , 10^5.

For example (from the problem statement), let k=6, a=0.2, and b=0.85.

Here the expression \sqrt{(ka+1)^2 + (kb+1)^2} gives the number 6.484. When rounded to the nearest integer, it becomes 6, which is equal to k. Therefore, Albert gets 6 points.

Yup – another expected value problem. But the concept of this problem is simple:

Hosted by imgur.com

For any value of k, there are certain values of a and b where Albert will get points. Otherwise, he does not get points. Here the gray blob is the area where he does get points, and the white area is where he does not.

The red square is the sample space of the random variables a and b, since 0<a<1 and 0<b<1.

The probability of Albert getting the points is the fraction of the square that’s gray.

The area of the red square is always 1. The area of the gray area varies with k. Thus the idea is to calculate the area of the gray area for each value of k from 1 to 10^5.

Now we can get to the details of the problem.

Notice that if round(k') = k, then k-0.5 < k' < k+0.5 (the difference between less-equal and less is irrelevant in this problem):

Hosted by imgur.com

In a similar way, the equation of the larger circle is (a+\frac{1}{k})^2 + (b+\frac{1}{k})^2 < (\frac{k+0.5}{k})^2.

These are equations of circles. A circle has the equation (x-a)^2 + (y-b)^2 = r^2 where (a,b) is the center and r is the radius.

Thus both the smaller and the larger circles have a center at (-\frac{1}{k},-\frac{1}{k}). The radii of the circles are different: the radius of the smaller circle is \frac{k-0.5}{k} and the radius of the larger circle \frac{k+0.5}{k}.

Plotting the inequalities:

Hosted by imgur.com

We can see now that this is a calculus problem! Indeed, the shaded area is the area of the larger circle inside the boxarea of the smaller circle inside the box.

Initially I tried integrating the circles directly, but I ran into problems while integrating and I was unable to correctly integrate the equation. Instead, I integrated an equivalent problem:

Hosted by imgur.com

It’s obvious that the areas do not change with this transformation.

So now for the smaller circle the equation is a^2 + (b+\frac{1}{k})^2 > (\frac{k-0.5}{k})^2. Solving for b, we get:

b > \sqrt{(\frac{k-0.5}{k})^2-a^2} - \frac{1}{k}

.. And this equation for the larger circle:

b < \sqrt{(\frac{k+0.5}{k})^2-a^2} - \frac{1}{k}

The variable of integration here is a, and k is only a constant. Therefore we can replace the radius with r and get this equation to integrate:

b = \sqrt{r^2-a^2} - \frac{1}{k}

With the help of an integration table, this can be integrated fairly easily.

Hosted by imgur.com

There’s one additional problem – computing the integrating ‘limits’ of the two circles. Let’s call them L_S and L_L for small and large respectively.

Both of them can be computed using the pythagorean theorem.

Hosted by imgur.com

The formula for L_S (it should be obvious by now what L_L should be):

L_S = \sqrt{(\frac{k-0.5}{k})^2-\frac{1}{k^2}}

Now we can get this really big formula for computing the area for any integer k:

\begin{array}{l} \int_{\frac{1}{k}}^{\sqrt{(\frac{k+0.5}{k})^2-\frac{1}{k^2}}} (\sqrt{(\frac{k+0.5}{k})^2-a^2}-\frac{1}{k}) da \\ - \int_{\frac{1}{k}}^{\sqrt{(\frac{k-0.5}{k})^2-\frac{1}{k^2}}} (\sqrt{(\frac{k-0.5}{k})^2-a^2}-\frac{1}{k}) da \end{array}

Also, fortunately for me, it is impossible for the upper bound to cross the circle at the top or right, I mean, we never have something like this -

Hosted by imgur.com

Anyways, here’s a Haskell implementation of the algorithm:

integral a r k = 0.5 * (a * sqrt (r^2 - a^2) + r^2 * atan (a/sqrt (r^2-a^2))) - a/k

area k = let area_1 1 = 0
             area_1 k = integral (sqrt (((k-0.5)/k)^2 - (1/k)^2)) ((k-0.5)/k) k - integral (1/k) ((k-0.5)/k) k
             area_2 k = integral (sqrt (((k+0.5)/k)^2 - (1/k)^2)) ((k+0.5)/k) k - integral (1/k) ((k+0.5)/k) k
         in  area_2 k - area_1 k

main = print . sum . zipWith (*) [1..] $ map area [1..10^5]

The code runs in about 3 seconds.

There seems to be an edge case where k=1 (probably an asymptote) so I have to specially make it be zero.

This problem was released 10 pm last night and I was tired, lol, but I still managed to solve it and get into the top 20.

Hosted by imgur.com


Project Euler 280 (Revisited)

March 27, 2010

A while ago I wrote this blog post on this problem. The previous code was able to produce the answer in about 20 minutes, using several gigabytes of RAM.

You may want to read the previous blog post first.

There is a simpler and more efficient solution (used by inamori, zwuupeape, and many others).

The new approach

We do not need to know the exact value of the solution, only an approximate. This approximation has to be accurate only to six digits.

Any any step, there is a certain probability that the ant will reach the final state in that step.

Let’s call this function P(s): the probability that the ant will reach the final state on step s:

Hosted by imgur.com

At very small values of s, the probability is zero because it takes a certain number of steps to get to the final position. 44 is the smallest number of steps to get from the initial to the final position, and the probability increases from there.

But it only increases until a certain point. After that point, the probability decreases. After all, if we run this for thousands of steps, how likely is the ant to still not have reached the final position?

However, the graph continues forever to the right. Even after thousands of steps, there is still a chance, however small, that the final position hasn’t been reached.

Let’s define Q(s) as the weighted average of the probabilities of each value of P(s) from 0 to s:

Hosted by imgur.com

The expected value and the answer for this problem would simply be the weighted average for the entire graph.

Since s must be an integer greater than zero, we don’t have to do any integration with curves. A bit misleading because I’ve drawn P(s) and Q(s) as curves, when they are not.

But because Q(s) cannot be expressed algebraically, we cannot evaluate it as it approaches infinity. Instead, we pick a point that gives us enough precision:

Hosted by imgur.com

Either that, or we keep calculating the weighted average while increasing s until we have enough precision.

Calculating the probabilities

In order to calculate a weighted average Q(s), we need to be able to calculate P(s); while P(s) is the probability function for the final state, we’ll need to know the probabilities for the in-between states as well.

The probabilities of the states on step s depend most upon the probabilities of the previous step, s-1.

Hosted by imgur.com

So if on step s-1 there’s a 0.6 probability that we’re on state X, then on step s we should have a 0.3 probability of being on either Y or Z as there’s an equal chance of moving to either state.

And if multiple states can go to our state, we can get the probability by adding up the probabilities:

Hosted by imgur.com

Coding a solution

At first we start on the initial state, with a probability of 1 (certain chance) of being on it.

On each step we generate a new set of probabilities using rules I described above, and so on.

A lot of code is stolen from my previous work on this problem. The State class is essentially unchanged in this version.

Here’s the code:

import java.util.*;

public class Main{

	public static void main(String[] args){

		final int ITERATIONS = 2501;

		// Initialize probability map.
		// The initial step is 0, containing the initial position with a probability of 1.
		Map<State,Double> prevProbs = new LinkedHashMap<State,Double>();
		prevProbs.put(State.INIT_STATE, 1.0);

		// Distributions for the final state.
		double[] finalValues = new double[ITERATIONS+1];
		Arrays.fill(finalValues, 0);

		// Keep the time.
		long startTime = System.currentTimeMillis();

		// Loop counter is actually the step number of the next step.
		for(int step=1; step<ITERATIONS; step++){

			// Make new probabilities.
			Map<State,Double> nextProbs = new LinkedHashMap<State,Double>();

			// Loop through old probabilities and add
			Set<State> stateSet = prevProbs.keySet();

			for(State state : stateSet){

				// Probability of being on this state right now
				double stateProb = prevProbs.get(state);

				// List of next states
				List<State> nextStates = nextStates(state);

				// Probability factor for each of the next states
				double multiplier = 1.0 / nextStates.size();

				// Added value to each of the next states
				double added = stateProb * multiplier;

				// Now go and update all the state probabilities.
				for(State nextState : nextStates){

					// If we've got the final state, handle it separately.
					if(nextState.equals(State.FINAL_STATE)){
						finalValues[step] += added;
						continue;
					}

					// Could be the first time.
					if(!nextProbs.containsKey(nextState))
						nextProbs.put(nextState, added);

					// Might not be the first time.
					else {
						// Old value, we need to add to it instead of replacing it.
						double previous = nextProbs.get(nextState);
						nextProbs.put(nextState, previous + added);
					}
				}
			}

			// Our new probabilities are the old probabilities for the next loop.
			prevProbs = nextProbs;

			System.out.println(step + " " + weightedAverage(finalValues));
		}

		System.out.println((System.currentTimeMillis() - startTime) + "ms");
	}

	// Calculates a weighted average of a double array with the indices as
	// values and the values in the array as weights.
	static double weightedAverage(double[] array){

		double sumValues = 0;
		for(int i=0; i<array.length; i++)
			sumValues += i * array[i];

		return sumValues;
	}

	// Returns a list of possible continuation states from the current one.
	static List<State> nextStates(State state){

		// Current and changing position of the ant
		int antX = state.ant % 5;
		int antY = state.ant / 5;

		// Whether it can go into each of the four directions (N,S,E,W respectively).
		boolean[] possibleDirs = new boolean[4];
		Arrays.fill(possibleDirs, true);

		// Take out some directions if it's in the corner.
		if(antY == 0) possibleDirs[0] = false; // Can't go north
		if(antY == 4) possibleDirs[1] = false; // Can't go south
		if(antX == 4) possibleDirs[2] = false; // Can't go east
		if(antX == 0) possibleDirs[3] = false; // Can't go west

		// Construct a list of continuations.
		List<State> nextStates = new ArrayList<State>();

		// Loop through the four directions.
		for(int i=0; i<4; i++){

			// Cannot go this direction.
			if( !(possibleDirs[i])) continue;

			int newAntX = antX;
			int newAntY = antY;

			// Modify direction.
			switch(i){
				case 0: newAntY--; break;
				case 1: newAntY++; break;
				case 2: newAntX++; break;
				case 3: newAntX--; break;
			}

			// Start constructing new state object.
			int oldAnt = state.ant; // old ant position
			int newAnt = 5*newAntY + newAntX;
			int[] board = state.board.clone();
			boolean carrying = state.carrying;

			// Carrying a seed. Notice that a square can contain
			// two seeds at once (but not more); seeds are indistinguishable
			// so we just need to keep track of the number of seeds
			// on each square.
			if(carrying){
				board[oldAnt] --;
				board[newAnt] ++;
			}

			// Drop off the seed.
			if(newAntY == 0 && board[newAnt] == 1 && carrying)
				carrying = false;

			// Pick up a new seed.
			if(newAntY == 4 && board[newAnt] == 1 && !carrying)
				carrying = true;

			// Treat the five done positions the same.
			if(donePosition(board))
				nextStates.add(State.FINAL_STATE);

			// Add it normally.
			else nextStates.add(new State(board, newAnt, carrying));
		}

		return nextStates;
	}

	// Is the board in the done position?
	static boolean donePosition(int[] b){
		return b[0]==1 && b[1]==1 && b[2]==1 && b[3]==1 && b[4]==1;
	}
}

class State{

	static final State INIT_STATE = new State(
	new int[]{
		0,0,0,0,0,
		0,0,0,0,0,
		0,0,0,0,0,
		0,0,0,0,0,
		1,1,1,1,1
	}, 12, false);

	// Consider all final states the same; there is
	// no ant position.
	static final State FINAL_STATE = new State(
	new int[]{
		1,1,1,1,1,
		0,0,0,0,0,
		0,0,0,0,0,
		0,0,0,0,0,
		0,0,0,0,0
	}, -1, true);

	// 25 board.
	int[] board;

	int ant;
	boolean carrying;

	State(int[] board, int ant, boolean carrying){
		this.board = board;
		this.ant = ant;
		this.carrying = carrying;
	}

	State(State s){
		this(s.board, s.ant, s.carrying);
	}

	public boolean equals(Object o){
		State s = (State) o;
		return Arrays.equals(s.board, board) && s.ant == ant && s.carrying == carrying;
	}

	public int hashCode(){
		return Arrays.hashCode(board) + ant;
	}

	// For debugging mostly.
	public String toString(){
		StringBuilder ret = new StringBuilder("\n");
		for(int i=0; i<5; i++){
			for(int j=0; j<5; j++){
				int pos = 5*i + j;
				if(ant == pos) ret.append("#");
				else ret.append(board[pos] >= 1? "+" : "-");
			}
			ret.append("\n");
		}
		return ret.toString();
	}
}

As we run this program, the answer becomes more and more accurate. We only need 2300 or 2500 steps to get the six digits of accuracy we need.

Running the program:

Hosted by imgur.com

This takes just over a minute to run — not bad at all.


Project Euler 280

March 6, 2010

Project Euler 280 is an interesting problem involving probability and combinatorics:

There is a 5×5 grid. An ant starts in the center square of the grid and walks randomly. Each step, the ant moves to an adjacent square (but not off the grid).

In each of the five bottom-most squares, there is a seed. When the ant reaches such a square, he picks up the seed (if he isn’t already carrying one).

When the ant reaches a square in the top-most row while carrying a seed, he drops off the seed (if there isn’t already a seed there).

This ‘game’ ends when all five seeds are in the five squares of the top-most row.

The problem asks for the expected (average) number of turns it would take to get from the initial to the ending position. It requires six digits of precision.

The Monto Carlo Approach

Perhaps the easiest way to tackle this problem is with a Monto Carlo simulation. This uses a computer to actually run the game many times, using random number generators.

This is my straight-forward Monto Carlo implementation in Java:

import java.util.*;

public class Main{
	public static void main(String[] args){
		Random r = new Random();
		int stepsum = 0;
		int tries = 0;
		while(true){
			stepsum += new Simulation(r).simulate();
			tries ++;
			System.out.println((double) stepsum / (double) tries);
		}
	}
}

class Simulation{

	static final int[] init = {
		0,0,0,0,0,
		0,0,0,0,0,
		0,0,0,0,0,
		0,0,0,0,0,
		1,1,1,1,1
	};
	State state;
	Random rand;
	int steps;

	Simulation(Random rand){
		state = new State(init, 12, false);
		this.rand = rand;
		steps = 0;
	}

	int simulate(){
		while(!done())
			step();
		return steps;
	}

	boolean done(){
		int[] b = state.board;
		return b[0]==1 && b[1]==1 && b[2]==1 && b[3]==1 && b[4]==1;
	}

	boolean step(){
		steps ++;

		int antX = state.ant % 5;
		int antY = state.ant / 5;
		int dir;

		while(true){
			// 0:N 1:S 2:E 3:W
			dir = rand.nextInt(4);
			if(antY==0 && dir==0) continue;
			if(antY==4 && dir==1) continue;
			if(antX==4 && dir==2) continue;
			if(antX==0 && dir==3) continue;
			break;
		}

		switch(dir){
			case 0: antY--; break;
			case 1: antY++; break;
			case 2: antX++; break;
			case 3: antX--; break;
		}

		int oldAnt = state.ant;
		state.ant = 5*antY + antX;

		if(state.carrying){
			state.board[oldAnt]--;
			state.board[state.ant]++;
		}

		if(antY == 0 && state.board[state.ant] == 1 && state.carrying){
			// drop off
			state.carrying = false;
		}

		if(antY == 4 && state.board[state.ant] == 1 && !state.carrying){
			// pick up
			state.carrying = true;
		}

		return true;
	}
}

class State{
	// 25 board.
	int[] board;

	int ant;
	boolean carrying;

	State(int[] board, int ant, boolean carrying){
		this.board = board.clone();
		this.ant = ant;
		this.carrying = carrying;
	}

	State(State s){
		this(s.board, s.ant, s.carrying);
	}

	public boolean equals(Object o){
		State s = (State) o;
		return Arrays.equals(s.board, board) && s.ant == ant;
	}

	public int hashCode(){
		return Arrays.hashCode(board) + ant;
	}
}

Running this program will produce something close to the final result. But this will not nearly be accurate enough to actually solve the problem: running this for a few hours would only produce two or three digits of precision, while six is required.

It’s not exactly impossible, as on the forums Rodinio claims to have run a similar Monto Carlo simulation on a computer cluster, giving enough precision to start making guesses for the last digits. To get so many digits on a Monto Carlo simulation requires an astronomical amount of tries.

Obviously this is not the most efficient approach.

Introducing the Markov Chain

A more efficient solution can be implemented using Markov Chains.

A Markov Chain solves the problem of discrete random processes. What is a Markov Chain, and how is it used?

Let’s say we have a finite group of states, S = \{ s_1, s_2, \cdots , s_n \}. Make a matrix M, of size n * n.

On each step (or move) of the process, the probability of going from state s_i to s_j is M_{ij}.

The process starts on an initial state. From each state there is a certain probability of going to every other state.

Notice that in a Markov Chain, what happens next does not depend on what has happened before. The only data is the current state you’re on.

I’ll explain all this with a very simple example:

On a 1×3 game board, Bob (the happy face) starts at square 1. Each step, he randomly moves to an adjacent square with equal probability. The game ends when he reaches square 3.

Of course the only time he has a ‘choice’ in which square to move to is when he’s on square 2 (half-half chance of moving to 1 and 3). If he’s already on square 1, he has a 100% chance of moving to square 2, and if he’s on square 3, well, the game is over.

We can represent this game like this:

Another way of showing the above would be something like this:

This is the matrix, of course:

\left[ \begin{array}{ccc} 0 & 1 & 0 \\ \frac{1}{2} & 0 & \frac{1}{2} \\ 0 & 0 & 1 \end{array} \right]

This type of matrix is called a Stochastic matrix. What’s special about this matrix is that each row of the matrix adds up to 1.

A state that can only return to itself is called an absorbing state. If in the matrix, M_{ii} = 1 for any value of i, then state i is an absorbing state.

Here, state 3 is an absorbing state.

The next step is very important: determining the time until absorption.

From our matrix, take out all of the absorbing rows, and all of the left columns to make it a square matrix:

It seems that by doing this we are losing information. But we’re actually not. We are treating each absorbing state as the same; and since each row originally added up to 1, the probability of going to the absorbing state is simply 1 minus the sum of the rest of the row.

This form is considered the canonical form of the stochastic matrix. Our matrix in canonical form looks like this:

\left[ \begin{array}{cc} 0 & 1 \\ \frac{1}{2} & 0 \end{array} \right]

Next we subtract it from the identity matrix.

\left[ \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right] - \left[ \begin{array}{cc} 0 & 1 \\ \frac{1}{2} & 0 \end{array} \right] = \left[ \begin{array}{cc} 1 & -1 \\ -\frac{1}{2} & 1 \end{array} \right]

Next we invert this matrix and get:

\left[ \begin{array}{cc} 1 & -1 \\ -\frac{1}{2} & 1 \end{array} \right]^{-1} = \left[ \begin{array}{cc} 2 & 2 \\ 1 & 2 \end{array} \right]

The expected absorption time from state s_n is simply the sum of row n in the above matrix.

So because we started on state 1, the expected length of the game is 2+2 or 4.

All of this seems like magic, because I haven’t given any proof for this. If you really want to prove that this method works, you can get a textbook about Markov chains.

“Markov Chains” by J.R.Norris (Cambridge 1997) is a decent textbook if you want to know more about this. I would recommend this as a less dense introduction to Markov chains.

An attempt using Markov Chains

Now that we know some basics of Markov chains, we need to apply it to solve the problem.

If you attach a Set or something to the Monto Carlo simulation program, you would find that there are 10270 distinct states for this problem. If we count the five final states as the same state, then we have 10266 states.

This is exactly the same as what we’ve just done. Except instead of a 3*3 matrices, we’re now working with matrices with tens of thousands of rows and columns. Woo hoo.

You can probably see here that this is going to be a problem. A 10265*10265 (final state not counted) matrix contains over 100 million elements. Most of these elements will be zero. If we store all of this data in double types, that’s 800mb of memory used, just for the matrix.

It’s easy to underestimate the size of this problem. But there’s a few optimizations we can make.

We could save some time and space by subtracting from the identity matrix at the same time as filling up the matrix.

Inverting a 10265*10265 matrix is no small task either. Instead there’s a simple optimization to avoid inverting such a huge matrix:

Remember that we want to find the inverse of the matrix, and get the sum of all the elements in the first row. Let M be our matrix and M' be the inverse.

What’s important here is that finding the sum of all the rows of a matrix is the same as multiplying by a column of ones. This produces a one-columned matrix, each cell representing the sum of the corresponding row.

If C is the column of ones, and S is the sum matrix (what we want), then M' * C = S. We can also write that as M * S = C.

See, this way we don’t have to compute the inverse. We just have to solve the above equation for S.

To avoid implementing the standard matrix solve function, I’ll be using the JAMA library for matrix manipulation.

Finally, the source code for my implementation of what I just discussed:

import java.util.*;
import Jama.*;

public class Main{

	// For a state, maps to all possible states from it.
	// All next-states have the same probability.
	static Map<State,List<State>> states = new LinkedHashMap<State,List<State>>();

	public static void main(String[] args){

		// Keep the time.
		long startTime = System.currentTimeMillis();

		// Initialize the state map.
		addAllStates(State.INIT_STATE);

		// To construct the matrix, we need to map all the states to a position
		// in the matrix (an integer). Because our map is already ordered, we
		// use its position in the map as the position in the matrix. In order
		// to avoid doing a linear search through the keyset to find its position,
		// we cache its position in a map.

		Set<State> keySet = states.keySet();
		Map<State,Integer> positionMap = new HashMap<State,Integer>();

		// Set up position map.
		Iterator<State> iterator = keySet.iterator();
		int pos = 0;
		while(iterator.hasNext()){
			positionMap.put(iterator.next(), pos);
			pos++;
		}

		// Finally we can begin constructing our matrix.
		// This takes up about 900mb of memory.
		double[][] matrix = new double[states.size()][states.size()];

		// Fill up the matrix.
		Collection<List<State>> values = states.values();
		Iterator<List<State>> rowIterator = values.iterator();

		for(int row=0; row<matrix.length; row++){

			// Fill up one row.
			List<State> thisRow = rowIterator.next();

			// What to fill with?
			double fill = 1.0 / thisRow.size();

			// An optimization: why not do the subtraction step here
			// instead of subtracting from identity matrix?
			matrix[row][row] = 1;

			for(State st : thisRow){
				// The final state doesn't count.
				if(st.equals(State.FINAL_STATE)) continue;

				int col = positionMap.get(st);
				matrix[row][col] = -fill;
			}
		}

		// Finishing up.
		Matrix bigMatrix = new Matrix(matrix);
		Matrix onesColumn = new Matrix(states.size(), 1, 1);
		Matrix sums = bigMatrix.solve(onesColumn);

		System.out.println(sums.get(0,0));
		System.out.println(System.currentTimeMillis() - startTime + "ms");

	}

	// Returns a list of possible continuation states from the current one.
	static List<State> nextStates(State state){

		// Current and changing position of the ant
		int antX = state.ant % 5;
		int antY = state.ant / 5;

		// Whether it can go into each of the four directions (N,S,E,W respectively).
		boolean[] possibleDirs = new boolean[4];
		Arrays.fill(possibleDirs, true);

		// Take out some directions if it's in the corner.
		if(antY == 0) possibleDirs[0] = false; // Can't go north
		if(antY == 4) possibleDirs[1] = false; // Can't go south
		if(antX == 4) possibleDirs[2] = false; // Can't go east
		if(antX == 0) possibleDirs[3] = false; // Can't go west

		// Construct a list of continuations.
		List<State> nextStates = new ArrayList<State>();

		// Loop through the four directions.
		for(int i=0; i<4; i++){

			// Cannot go this direction.
			if( !(possibleDirs[i])) continue;

			int newAntX = antX;
			int newAntY = antY;

			// Modify direction.
			switch(i){
				case 0: newAntY--; break;
				case 1: newAntY++; break;
				case 2: newAntX++; break;
				case 3: newAntX--; break;
			}

			// Start constructing new state object.
			int oldAnt = state.ant; // old ant position
			int newAnt = 5*newAntY + newAntX;
			int[] board = state.board.clone();
			boolean carrying = state.carrying;

			// Carrying a seed. Notice that a square can contain
			// two seeds at once (but not more); seeds are indistinguishable
			// so we just need to keep track of the number of seeds
			// on each square.
			if(carrying){
				board[oldAnt] --;
				board[newAnt] ++;
			}

			// Drop off the seed.
			if(newAntY == 0 && board[newAnt] == 1 && carrying)
				carrying = false;

			// Pick up a new seed.
			if(newAntY == 4 && board[newAnt] == 1 && !carrying)
				carrying = true;

			// Treat the five done positions the same.
			if(donePosition(board))
				nextStates.add(State.FINAL_STATE);

			// Add it normally.
			else nextStates.add(new State(board, newAnt, carrying));

		}

		return nextStates;
	}

	// Recursively add all continuation states.
	static void addAllStates(State state){

		// Don't add if it's already added.
		if(states.containsKey(state)) return;

		List<State> nexts = nextStates(state);
		states.put(state, nexts);

		// Recurse (but not if we've reached the final state).
		for(State next : nexts)
			if( !(next.equals(State.FINAL_STATE)))
				addAllStates(next);

	}

	// Is the board in the done position?
	static boolean donePosition(int[] b){
		return b[0]==1 && b[1]==1 && b[2]==1 && b[3]==1 && b[4]==1;
	}
}

class State{

	static final State INIT_STATE = new State(
	new int[]{
		0,0,0,0,0,
		0,0,0,0,0,
		0,0,0,0,0,
		0,0,0,0,0,
		1,1,1,1,1
	}, 12, false);

	// Consider all final states the same; there is
	// no ant position.
	static final State FINAL_STATE = new State(
	new int[]{
		1,1,1,1,1,
		0,0,0,0,0,
		0,0,0,0,0,
		0,0,0,0,0,
		0,0,0,0,0
	}, -1, true);

	// 25 board.
	int[] board;

	int ant;
	boolean carrying;

	State(int[] board, int ant, boolean carrying){
		this.board = board;
		this.ant = ant;
		this.carrying = carrying;
	}

	State(State s){
		this(s.board, s.ant, s.carrying);
	}

	public boolean equals(Object o){
		State s = (State) o;
		return Arrays.equals(s.board, board) && s.ant == ant && s.carrying == carrying;
	}

	public int hashCode(){
		return Arrays.hashCode(board) + ant;
	}

	// For debugging mostly.
	public String toString(){
		StringBuilder ret = new StringBuilder("\n");
		for(int i=0; i<5; i++){
			for(int j=0; j<5; j++){
				int pos = 5*i + j;
				if(ant == pos) ret.append("#");
				else ret.append(board[pos] >= 1? "+" : "-");
			}
			ret.append("\n");
		}
		return ret.toString();
	}
}

Running it:

So my program takes about 19 minutes to run, using up over 2GB of ram. The construction of the matrix only takes about 2 seconds, and the rest of the time is used on solving the huge matrix.

JAMA is probably not the fastest matrix library; the running time might be cut down a bit if we use libraries designed for sparse matrices.

But to get this to run in under a minute, a completely different approach would be needed. For now, though, I’m pretty happy with getting it in 19 minutes.


Projecteuler-solutions and Polymath

February 28, 2010

A couple of months ago (May 2009), I created a project on Google Code to share Project Euler solutions.

At first I started with only a hundred or so answers, but I found more, and people began contributing them, and by July or August I had the answer to every Project Euler problem.

But problems were still coming out, and regular posters like jeneshicc and inamori were less willing to publish their solutions.

So in September, I set up a forum to discuss solutions. This quickly turned into a polymath project.

A polymath project is a large collaborative project, where many people work together on a math problem. The idea was perhaps first inspired by a blog post by Gowers.

There were previous successful attempts at polymath projects. One I found particularly interesting was this one by Terry Tao, which is about Q6 of the International Math Olympiad. This “mini-polymath” project is similar to Projecteuler-solutions in the way that both are solving already-solved problems with an unknown (but existent) solution.

I feel that this project (specifically its forums) provides a good example of a successful polymath project. From September to present, the community has solved 24 Project Euler problems. Some easier problems were mostly of individual effort (just posting the answer), while more difficult problems were of a group effort. A particularly hard problem (257 I think) took the community over a week.

Even though the project has existed for months, I still get emails and forum posts suggesting that I take the whole thing down.

Experiences with Projecteuler-solutions

Overall, I think the project is very successful. Right now, there are 200 registered users (albeit a lot of them spammers), and over 900 posts.

There are perhaps 10 or 15 serial contributers: people who come up with good insights for multiple problems.

Then there are many more people who occasionally come up with something useful. Of course, any new information is welcome in such a group effort.

There have been a couple of trolls. Although they do not contribute to problems, they provoke some entertaining discussions. As they do not really get in the way of solving the problems, I don’t ban them.

More recently, there has been a lot of spammers, linking to porn sites and similar crap. They were a bit harder to deal with, since I hosted the forums on forumer, which does not allow me to install better Captcha’s to stop the bots. What I eventually did was to limit thread creation to a specific group, which stopped the spammers.

Here are some impressions I got from being a moderator:

  1. The group effort is very powerful. No matter how hard the problem is, I’m confident that it will be solved eventually. Any individual will give up after some amount of time, but it’s much harder for the entire group to give up. Until it is solved, someone is always trying the problem.
  2. Much of the work is done offline. On one hand, we have a lot of low-quality observations. On the other hand, we have posts that make giant leaps of progress, but are done by one person. Indeed, for easier problems one person can come up with the entire solution. But for harder problems, multiple of these ‘giant leaps’ are required, sometimes by different people building up from the work done by others.
  3. The forum is probably the best format for this kind of project. The other polymath projects used the wiki format. I think forums are better because ideas are ordered chronologically making it easier to view progress, whereas in a wiki the ideas are organized by topic and many mini-discussions arise, while less is being discussed about the problem as a whole.
  4. The community usually does fine without any moderation. Members can create topics, work on the problem, and come up with solutions on their own. Even trolls are dealt with by the community, and I’m rarely forced to ban them.

Perhaps Forumer is not the best service to host such a project. I am not able to install most plugins, such as support for code syntax highlighting, for Latex, and for custom Captcha’s. But since my parents won’t pay for my own server hosting, I think Forumer is a good free service.


Follow

Get every new post delivered to your Inbox.

Join 69 other followers