Prime time programming

The purpose of computing is insight, not numbers.” – Richard Hamming

My interest in pursuing this series is not in the numbers, nor the Mathematics, but in the process of constructing reliable solutions. I’m also learning Scala, and want to understand how its hybrid OO/FP nature may offer fresh perspectives to my pursuit of programming.

The largest factor

The third problem from Project Euler asks for the largest prime factor of a number N = 600,851,475,143. One approach would be produce the prime factors of N and return the max of that set. An obvious way to do that is to generate primes and try dividing each into N, but:

  1. Generating (or testing for) primes requires more code than generating divisors of N.

    We can find all factors of N by testing divisors 2 … n (where n = √ N). The sieve of Eratosthenes yields all primes at most n with effort of O(n log n), compared with O(n) for simply dividing N by all of the candidate divisors.

  2. Factoring just one single number doesn’t appear to provide a return on the investment of building the collection of primes.

    Had the problem called for factoring multiple numbers of comparable size, then the cost of constructing a set of primes could have been amortized over multiple uses. I’ll revisit this later.

  3. The smallest factor of N must be prime.

    The smallest factor of N can’t have any factors smaller than itself (else it wouldn’t be the smallest). If it doesn’t have any factors smaller than itself, it’s prime.

We can base a simple design on those ideas:

  • Removing all occurrences of the smallest factor of N leaves a quotient. Unless that quotient is 1, its largest prime factor is the largest prime factor of N. If the quotient is 1, the last factor removed is the largest prime factor of N.
  • We find factors by checking increasingly larger candidate divisors, beginning with 2. We don’t need to restart the divisors for a new quotient, because all smaller divisors have already been eliminated.
  • If the square of the current candidate divisor exceeds the current quotient then that quotient itself is prime (because all of its other possible divisors have already been eliminated).

First version

It’s more trouble to write out those bullet points in English than it is to write the code:

  def lastFactor(t: Long) = {
    def loop(q: Long, f: Long, d: Long): Long =
      if (q == 1)
      else if (q < d * d)
      else if (q % d == 0)
        loop(q / d, d, d)
        loop(q, f, d + 1L)
    loop(t, 1L, 2L)

As an aside, let me mention an issue of style. I’ve noticed that published functional programs tend to be written using short (some would say “cryptic”) variable names. Here’s the same definition, using words instead of initials:

  def lastFactor(target: Long) = {
    def loop(quotient: Long, factor: Long, divisor: Long): Long =
      if (quotient == 1)
      else if (quotient < divisor * divisor)
      else if (quotient % divisor == 0)
        loop(quotient / divisor, divisor, divisor)
        loop(quotient, factor, divisor + 1)
    loop(target, 1, 2)

My speculation is that the functional style tends to emphasize the structure of the function itself; longer names seem to put the emphasis on the variables. A “meatier” style makes it harder to see the “bones”. As I experiment with various approaches to a problem, I find myself using shorter names in an attempt to highlight the structural differences. I’ll stay with that style here (with apologies to anyone who finds it harder to read).

Back on the main topic, notice that the inner function loop deals with two different issues in the two tail-recursive calls:

  1. loop(q / d, d, d) ensures that all occurrences of the current trial divisor are factored out (retaining the current divisor as the new largest factor), and
  2. loop(q, f, d + 1L) proceeds to the next trial divisor (with the current quotient and factor).

As I mentioned in an earlier post, Scala compiles direct tail recursion into tight, fast bytecode that doesn’t grow the activation stack. On my laptop, it finds the largest prime factor of 600,851,475,143 in about 66µs. (If you don’t want to waste a few microseconds running the code yourself, the answer is 6,857.)

By way of contrast, an equivalent function using a pre-computed list of primes requires only about 15µs. However, producing that list (all primes from 2 to 775,147) via a simple sieve takes about 61ms (yes, that’s 61,000µs). Clearly we’d need to use that list quite a few times to recover the investment.

A factor saved…

Even if I didn’t know that future Project Euler problems will revisit factoring (my experiments are ahead of my writing), I would prefer not to throw information away. How much would the solution change if we decided to keep all factors instead of just the largest? Not much; just a few adjustements to the inner loop function:

  • Parameter f: Long (the most recent factor) becomes fs: List[Long] (a list of all factors).
  • The initial value for fs is Nil (the empty list).
  • Everywhere f gets a new value, fs gets a value consed onto the front.

The changes are highlighted in the new function below:

  def allFactors(t: Long) = {
    def loop(q: Long, fs: List[Long], d: Long): List[Long] =
      if (q == 1)
      else if (q < d * d)
        q :: fs
      else if (q % d == 0)
        loop(q / d, d :: fs, d)
        loop(q, fs, d + 1)
    loop(t, Nil, 2)

Running that version with the specified value of N returns the complete list of factors…

scala> allFactors(600851475143L)
res19: List[Long] = List(6857, 1471, 839, 71)

…and does so almost for free, taking about 67µs instead of 66µs on my laptop. (I’ve only run a few million tests, so that small a difference is almost statistical noise.) That’s not really surprising, as the only extra overhead is a call to :: as each factor is saved. Even if our target had been a power of two (producing the most factors) there would have been only about 40 factors on the list.

…is a factor folded

“Cut-and-paste” is regarded as bad practice in OO. The DRY principle says that if we find ourselves repeating the same code, or pattern of code, we should extract a reusable form of the common logic. The same principle applies in FP.

Notice that lastFactor and allFactors have exactly the same higher-level structure; the differences are summarized in the following table, which also suggests a general case for some unknown type T.

  lastFactor allFactors general case
Long Long Long
Long List[Long] T
1 Nil ?
LongLongLong LongList[Long]List[Long] LongTT
(a, b) → a
(f, fs) → f :: fs

So, just for fun here is a generic Scala function that processes the prime factors of an argument, using a caller-suppllied "zero" and composition function, with the differences highlighted:

  def foldFactors[T](t: Long)(z: T)(f: (Long, T) => T) = {
    def loop(q: Long, acc: T, d: Long): T =
      if (q == 1)
      else if (q < d * d)
        f(q, acc)
      else if (q % d == 0)
        loop(q / d, f(d, acc), d)
        loop(q, acc, d + 1)
    loop(t, z, 2)

We can get the "last" factor and the complete list of factors by providing the appropriate zero and composition function for each.

  scala> foldFactors[Long](n)(1)((a,b) => a)       
  res32: Long = 6857
  scala> foldFactors[List[Long]](n)(Nil)(_ :: _)
  res33: List[Long] = List(6857, 1471, 839, 71)

To be realistic, allFactors is really our best solution. It returns the factors as List[Long]], which means that we get foldLeft for free, along with everything else that List provides. But this was an irresistible opportunity to point out that higher-order functions and type parameterization address the DRY principle that we already know from OO.

To be continued

I think that's enough to digest for now. There's another punch line that resulted from exploring this problem, but it deserves its own write-up. If you want a head-start on it, think about the fact that the result of allFactors is a bit inconvenient for some purposes, especially when the argument has many factors, as in:

  scala> allFactors(86400)
  res35: List[Long] = List(5, 5, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2)

Contrast that with the "normal" Mathematical representation of 52 × 33 × 27 for a hint. That's where I found myself reaching for the hybrid nature of Scala.

I noticed that Basildon Coder had a post on this same problem. I held off reading it until after working through my own approach, so that I could look at the similarities and differences in two independently-created solutions. Perhaps you'll find that interesting as well.

Trackbacks are closed, but you can post a comment.


  • Russ  On 2008-04-21 at 14:14

    That really is a first-class analysis, well played. I particularly liked the use of a table to identify a general case for the foldFactors function. For the problems I’ve tackled so far the code has been fairly disposable (with the possible exception of a (quite limited) prime sieve), but you’ve dug out a couple of useful nuggets here I think. Can’t wait for the next one!

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: