Project Euler 58: Investigate the number of primes that lie on the diagonals of the spiral grid

Project Euler 58: Investigate the number of primes that lie on the diagonals of the spiral grid

Written by on 27 August 2011

Topics: Project Euler

In Problem 58 of Project Euler we are revisiting the spiral we worked on in Problem 28. I have since then learned that it is called an Ulam spiral after the Polish born American mathematician Stanislaw Ulam. The problem reads:

Starting with 1 and spiralling anticlockwise in the following way, a square spiral with side length 7 is formed.

37 36 35 34 33 32 31
38 17 16 15 14 13 30
39 18  5  4  3 12 29
40 19  6  1  2 11 28
41 20  7  8  9 10 27
42 21 22 23 24 25 26
43 44 45 46 47 48 49

It is interesting to note that the odd squares lie along the bottom right diagonal, but what is more interesting is that 8 out of the 13 numbers lying along both diagonals are prime; that is, a ratio of 8/13  62%.

If one complete new layer is wrapped around the spiral above, a square spiral with side length 9 will be formed. If this process is continued, what is the side length of the square spiral for which the ratio of primes along both diagonals first falls below 10%?

In problem 28 we sought and found a closed form solution to finding the sum of the diagonals. In this problem we need to investigate each corner of the spiral. So we need to generate each corner an check it.  I have found no other way than brute force to create the corners and check them. I found no pattern in the primes except that they never occurs in the lower right diagonal. The numbers there are always perfect squares due to the nature of the spiral.

Generating the corner numbers

For now let us assume that we have a method called isPrime, which returns true if the number is a prime, and false otherwise. Such that we can focus on the main loop. It is really simple though, so lets jump right into it.

int noOfPrimes = 3;
int sl = 2;
int c = 9;

while( ((double)noOfPrimes)/(2*sl+1) > 0.10){
    sl += 2;
    for(int i = 0; i < 3; i++){
        c += sl;
        if(isPrime(c)) noOfPrimes++;
    }
    c+= sl;
}

As we can see there are two loops. An inner loop that traverses the corners of the current loop in the spiral and checks the primality of each corner and an outer loop that increases the side length as the spiral grows.
The length between two corners is the side length minus 1.

As we can also note the code just not check the lower right diagonal as we know it will never be a prime.

Checking primality of the corners

I quickly discovered that this operation is the time consuming one. I started out with an implementation that used trial division which I copied from Problem 51. The code does the job but the time was not impressive.

The sidelength of the spiral when the ratio falls below 10% is 26241
Solution took 258 ms

I made a trial division based presieving primes to create a significant speed up to the algorithm, but I decided to look into the Rabin-Miller primality test, which is a probability based method, but it can rather easily be made deterministic numbers up to 3 * 1014. This upper bound for determinism should be more than enough for this purpose.  I wont go into details with how it works here, since I sure I can’t give a satisfying account of that yet. It might come later.

However, with my implementation of it (which you can see in the source code) the result becomes.

The sidelength of the spiral when the ratio falls below 10% is 26241
Solution took 30 ms

Wrapping up Problem 58

I don’t think the problem in question was particularly difficult to solve. However, we need an efficient primality check. For this purpose I made an implementation of the Rabin-Miller primality check.

You can as always dive into the source code. If you have questions or comments to the solution or the code, or maybe another solution method. Please let me know and let us have a discussion about it.

Blog image

The header image for this post is kindly shared under the creative commons license by PBoothe. It shows the Ulam spiral with the primes marked as dots. The colour of the dots are chosen randomly according to the creator.

It is a bit difficult to see it on the section I have cut out, but it is very clear on the original image – There are certain lines forming of the primes, these can be represented by quadratic polynomials.

9 Comments For This Post I'd Love to Hear Yours!

  1. SuprDewd says:

    Thanks for the Miller-Rabin primality test. I have my own implementation of it based on Java’s BigInteger.isProbablePrime() method, and it’s not performing nearly as good as your implementation.

  2. Kristian says:

    Hi SuprDewd

    I saw your Rabin-Miller implementation and I think that the main difference is that I implemented the modular power method while you use the one that comes with the BigInteger class.

    The BigInteger class definitely has performance issues compared to the integral types like int and longs. However, in cases where you would normally use Rabin-Miller for testing primality – such as cryptology – the numbers are so large that you would need the BigInteger class. So I don’t think there is anything wrong with the implementation you made.

    /Kristian

  3. SuprDewd says:

    I changed your code to work with a BigInteger, and it only took 250ms (vs. your 70ms).
    I think there’s something more than just the BigInteger.ModPow making my code slower. And I also think there are (non-performance) issues with my code.

  4. Kristian says:

    hmmm okay. You might even try to replace my implementation with the native BigInteger.ModPow and you might get a further speed up. What the differences are then I don’t know.

    /Kristian

  5. Muhammad al-Khwarizmi says:

    I found this terrifically useful for my Python solution, which runs in 2 seconds after paring away tons of needless complexity:

    from __future__ import division
    import copy
    from my_euler import *
    
    def proportion_of_primes(total_primes, n):
        return total_primes / ((n - 1) * 4 + 1)
    
    n = 2
    last_corner = 9
    total_primes = 3
    
    while proportion_of_primes(total_primes, n) > 0.1:
        n += 1
    
        for i in xrange(4):
            last_corner = last_corner + (n - 1) * 2
    
            if is_prime(last_corner):
                total_primes += 1
    
    # side length is desired
    print 2 * n - 1
  6. Kristian says:

    Glad you found it useful. Was it the way to generate the corner numbers, or the Rabin-Miller part you found useful?

  7. Shobhit says:

    Thanks for the Miller-Rabin primality test.
    I would look into it. I am not sure if I would be able to follow it completely.
    This one was pretty simple. The last challenging one was Poker Hands. I am looking for a tough one now.

  8. jiang zhiqiang says:
    using System;
    using System.Collections;
    
    class ProjectEuler058
    {
      const  int N = 30001;
      static BitArray primes;
    
      static void PrimeSieve(int n)
      {
        primes = new BitArray(n + 1, true);
        primes[0] = primes[1] = false;
        int m = (int)Math.Sqrt(n);
        for (int i = 2; i &lt;= m; i++)
          if (primes[i])
            for (int k = i &lt;&lt; 1; k &lt;= n; k += i)
              primes[k] = false;
      }
    
      static int SpiralCalc()
      {
        for (int n = 0, i = 1, k = 2; ; i += k, k += 2)
        {
          if (primes[i += k]) n++;
          if (primes[i += k]) n++;
          if (primes[i += k]) n++;
          if (n * 5 &lt;= k) return k + 1;
        }
      }
    
      static void Main()
      {
        PrimeSieve(N * N);
        Console.WriteLine(SpiralCalc());
      }
    }
    
  9. Juan says:

    Great blog… I use to read you and dreamshire to get a better insight in Project Euler problems and at the time I read this problem I knew it was necessary to use a “primality test”, but after reading and get Fermat Primality Test and later Miller-Rabin Primality Test, I failed to test it manually… I had to see it in code, and it really helped… worked like a charm.

    Submitting the answer was a little bit bittersweet, but I prefer to think that I learnt something instead of the taste of cheating hehehe.

    Thanks!

Leave a Comment Here's Your Chance to Be Heard!

You can use these html tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

You can use short tags like: [code language="language"][/code]

The code tag supports the following languages: as3, bash, coldfusion, c, cpp, csharp, css, delphi, diff,erlang, groovy, javascript, java, javafx, perl, php, plain, powershell, python, ruby, scala, sql, tex, vb, xml

You can use [latex]your formula[/latex] if you want to format something using latex

This site uses cookies.