Jesse Tilly at Memphis JUG

This month’s Memphis Java Users’ Group meeting featured Jesse Tilly of IBM Rational Software, who spoke to us on static analysis. He will be doing a more product-intensive session, “What is IBM® Rational® Software Analyzer® Telling Me?”, at the upcoming IBM Rational Software Conference. (Don’t be misled by all those “circle-R”s; I just linked to the title from the conference web site.)

For our meeting, Jesse left the branding iron at home. He began with an overview of the history and benefits of static analysis. The major portion of the presentation offered a practical approach to analysis as part of a development project, including a detailed how-to on interpreting and using analysis results. Jesse finished with return to history, drawing unexpected parallels with the analysis of Enigma traffic at Bletchley Park during WWII—the background for Allen Turing’s later theoretical work that led to the computers we program today.

Because Jesse had an early flight, our regular door-prize drawing followed his presentation. In our lightning talk segment, Matt Stine introduced Morph AppSpace, I presented on “Structured Functional Programming” (pdf here), and Walter Heger gave a quick look at jGears.

---

Recommended reading:

Encryption and cryptanalysis are deeply entwined with computing, whether in history(Codebreakers: The Inside Story of Bletchley Park) or in imagination (Cryptonomicon).

Two highly-respected tools for static analysis in Java are FindBugs and PMD; both web sites offer excellent documentation and other reference material.

BuilderBuilder: The Model in Haskell

This post describes the first model in Haskell for the BuilderBuilder task. We will develop the model incrementally until we have rough parity with the Java version.

I’m experimenting with ways to distinguish user input from system output in transcripts of interactive sessions. This time I’m trying color, using a medium blue for output. I will appreciate feedback on whether that works for you.

Step one: defining and using a type

The simplest possible Haskell version of our model for a Java field is:

data JField = JField String String

However, the PMNOPML (“Pay Me Now Or Pay Me Later”) principle says that we’ll regret it if we stop there. In fact, later comes quickly.

We can create an instance of JField in a source file:

field1 = JField "name1" "Type1"

To do the same in a ghci session, prefix each definition with let, as in:

*Main> let field1 = JField "name1" "Type1"

Step two: showing the data

Trying to look at the instance yields a fairly opaque error message.

*Main> field1

<interactive>:1:0:
    No instance for (Show JField)
      arising from a use of `print' at <interactive>:1:0-5
    Possible fix: add an instance declaration for (Show JField)
    In a stmt of a 'do' expression: print it

Remember that in Java the default definition of toString() returns something like com.localhost.builderbuilder.JFieldDTO@53f67e; that’s also obscure at first glance. Haskell just goes a bit further, complaining that we haven’t defined how to show a JField instance. We can ask for a default implementation by adding deriving Show to a data type definition:

data JField = JField String String deriving Show

After loading that change, we get back a string that resembles the field’s defining expression:

*Main> field1
JField "name1" "Type1"

Step three: referential transparency

Our first model represented a Java class by its package, class name, and enclosed fields. The Haskell equivalent is:

data JClass = JClass String String [JField] deriving Show

The square brackets mean “list of …”, so a JClass takes two strings and a list of JField values. I’ll say more about lists in a moment, but first let’s deal with referential transparency.

We can build a class incrementally:

field1 = JField "name1" "Type1"
field2 = JField "name2" "Type2"
class1 = JClass "com.sample.foo" "TestClass" [field1, field2]

or all at once:

class1 = JClass "com.sample.foo"
                "TestClass"
                [   JField "name1" "Type1" ,
                    JField "name2" "Type2"
                ]

and get the same result:

*Main> class1
JClass "com.sample.foo" "TestClass" [JField "name1" "Type1",JField "name2" "Type2"]

As mentioned previously, only one of those definitions of class1 can go in our program. To Haskell, name = expression is a permanent commitment. From that point forward, we can use name and expression interchangeably, because they are expected to mean the same thing. That expectation would break if we were allowed to give name another meaning later (in the same scope).

Consequently, we can define a class using previously defined fields, or we can just write everything in one definition, nesting the literal fields inside the class definition. As we’ll see later, this also has implications for how we write functions; a “pure” function and its definition are also interchangeable.

Step four: lists

The array is the most fundamental multiple-valued data structure in Java; the list plays a corresponding role in Haskell. In fact, lists are so important that there are a few syntactical short-cuts for dealing with lists.

• Type notation: If t is any Haskell type, then [t] represents a list of values of that type.
• Empty lists: Square brackets with no content, written as [], indicate a list of length zero.
• Literal lists: Square brackets, enclosing a comma-separated sequence of values of the same type, represent a literal list.
• Constructing lists: The : operator constructs a new list from its left argument (a single value) and right argument (a list of the same type).

For example, ["my","dog","has","fleas"] is a literal value that has type [String] and contains four strings. "my":["dog","has","fleas"] and "my":"dog":"has":"fleas":[] are equivalent expressions that compute the list instead of stating it as a literal value.

By representing the fields in a class with a list, we achieve two benefits:

• The number of fields can vary from class to class.
• The order of the fields is significant.

Step five: types and records

Given a JField, how do we get its name? Or its type? We can define functions:

fieldName (JField n _) = n
fieldType (JField _ t) = t

and do the same for the JClass data:

package   (JClass p _ _ ) = p
className (JClass _ n _ ) = n
fields    (JClass _ _ fs) = fs

but all that typing seems tiresome.

Before solving that problem, let’s note two other limitations of our current implementation:

• Definitions using multiple String values leave us with the burden of remembering the meaning of each strings.
• The derived show method leaves us with a similar problem; it doesn’t help distinguish values of the same type.

If you suspect that I’m going to pull another rabbit out of Haskell’s hat, you’re right. In fact, two rabbits.

Type declarations

We can make our code more readable by defining synonyms that help us remember why we’re using a particular type. By adding these definitions:

type Name     = String
type JavaType = String
type Package  = String

we can rewrite our data definitions to be more informative:

data JField = JField Name JavaType deriving Show
data JClass = JClass Package Name [JField] deriving Show

Record syntax

The second rabbit is a technique to get Haskell to do even more work for us. We represent each component of a data type as a name with an explicit type—all in curly braces, separated by commas:

data JField = JField {
     fieldName :: Name ,
     fieldType :: JavaType
} deriving Show

data JClass = JClass {
     package   :: Package ,
     className :: Name ,
     fields    :: [JField]
} deriving Show

When we use this syntax, Haskell creates the accessor functions automagically, and enables a more explicit and flexible notation to create values. All of these definitions:

field1  = JField "name1" "Type1"
field1a = JField {fieldName = "name1", fieldType = "Type1"}
field1b = JField {fieldType = "Type1", fieldName = "name1"}

produce equivalent results:

*Main> field1
JField {fieldName = "name1", fieldType = "Type1"}
*Main> field1a
JField {fieldName = "name1", fieldType = "Type1"}
*Main> field1b
JField {fieldName = "name1", fieldType = "Type1"}

Step last: that’s it!

We have covered quite a bit of ground! The complete source code for the model appears at the end of this post. With both Java and Haskell behind us, we have most of the basic ideas we’ll need for the Erlang and Scala versions.

---

Recommended reading:

Real World Haskell, which is also available
for the Amazon Kindle
or on-line at the book’s web site. I really can’t say enough good things about this book.

---

The current BuilderBuilder mode in Haskell, along with sample data, is this:

-- BuilderBuilder.hs

-- data declarations

type Name     = String
type JavaType = String
type Package  = String

data JField = JField {
     fieldName :: Name ,
     fieldType :: JavaType
} deriving Show

data JClass = JClass {
     package   :: Package ,
     className :: Name ,
     fields    :: [JField]
} deriving Show

-- sample data for demonstration and testing

field1  = JField "name1" "Type1"
field1a = JField {fieldName = "name1", fieldType = "Type1"}
field1b = JField {fieldType = "Type1", fieldName = "name1"}

field2 = JField "name2" "Type2"

class1 = JClass "com.sample.foo" "TestClass" [field1, field2]

studentDto = JClass {
    package   = "edu.bogusu.registration" ,
    className = "StudentDTO" ,
    fields    = [
        JField {
            fieldName = "id" ,
            fieldType = "String"
        },
        JField {
            fieldName = "firstName" ,
            fieldType = "String"
        },
        JField {
            fieldName = "lastName" ,
            fieldType = "String"
        },
        JField {
            fieldName = "hoursEarned" ,
            fieldType = "int"
        },
        JField {
            fieldName = "gpa" ,
            fieldType = "float"
        }
    ]
}

---

Updated 2009-05-09 to correct formatting and add category.

BuilderBuilder: Haskell Preliminaries

The next step in the BuilderBuilder project is to develop a model in Haskell that is analogous to the Java model in the previous post. This post will introduce just enough Haskell to get started; the next post will get into the BuilderBuilder model.

Environment:

I’m using GHC 6.10.1, obtained from the Haskell web site. There are a variety of platform-specific binaries; I used the classic configure/make/install process on OSX. (For Java programmers, make is what we used instead of ant back in the Jurassic era.) Consult the Haskell Implementations page for details on obtaining Haskell for your preferred platform.

The complete development environment consists of two windows: one running a text editor, and the other running ghci, the interactive Haskell shell that comes with GHC.

Haskell introduction:

Use your text editor to create a file named bb1.hs with this content:

-- bb1.hs

-- simplest possible data declarations

data JField = JField String String

-- sample data for demonstration and testing

field1 = JField "id" "String"

-- sample function

helloField :: JField -> String
helloField (JField n t) = "Hello, " ++ n ++ ", of type " ++ t

Then run ghci as follows, where user input is underlined:

your-prompt-here$ ghci
GHCi, version 6.10.1: http://www.haskell.org/ghc/    for help
Loading package ghc-prim ... linking ... done.
Loading package integer ... linking ... done.
Loading package base ... linking ... done.
Prelude> :l bb1.hs
[1 of 1] Compiling Main             ( bb1.hs, interpreted )
Ok, modules loaded: Main.
*Main> helloField field1
"Hello, id, of type String"
*Main>

We started ghci, told it to load our source file (the :l … line), and then invoked the helloField function on the sample field. Now let’s examine the Haskell features used in that code. The lines beginning with double-hyphens are comments, and will be ignored in the description.

Defining data types

Because Haskell emphasizes functions, it’s no surprise that the syntax for defining data types is very lightweight. The Java BuilderBuilder model represents a field with two strings, one for the name and one for the type. The simplest possible Haskell equivalent is:

data JField = JField String String

This defines a data type named JField. It has a constructor (also named JField) that takes two strings, distinguished only by the order in which they are written.

Defining values

The next line of code defines an instance of this type:

field1 = JField "id" "String"

The equal sign means “is defined as“. That statement defines field1 as the instance of JField constructed on the right-hand side. It is not declaring and initializing a mutable variable. Within the current scope, attempting to redefine field1 will produce an error. (More about scope later.)

Defining functions

Finally, we have a simple function that converts a JField to a String.

helloField :: JField -> String
helloField (JField n t) = "Hello, " ++ n ++ ", of type " ++ t

Everything in Haskell has a type, including functions. The double colon means “is of type“, so the type of helloField is function from JField to String.

The value of applying helloField to a JField containing strings n and t is defined by the expression on the right-hand side. Haskell regards strings as lists of characters; the ++ operator concatenates lists of any type. The names n and t are only meaningful within that definition, similar to the local variables in this Java fragment:

public static String helloField(IJField f) {
    String n = f.getName();
    String t = f.getType();
    return "Hello, " + n + ", of type" + t;
}

Type inference

Java requires that we explicitly declare the local variables as type String. But in Haskell, because JField is specified to have two String values, the compiler can infer the types of n and t In fact, the entire first line of helloField is not necessary. The defining equation in the second line explicitly uses a JField on the left and constructs a String on the right. Therefore, the compiler can infer JField -> String as the type of the function. Haskell’s type inference allows us to write very compact code without giving up strong, static typing.

To see that in action, add the following line to the end of your bb1.hs file:

hiField (JField n _) = "Hi, " ++ n

(The underscore is a wild card, showing the presence of a second value but indicating that we don’t need it in this function.)

Reloading bb1.hs in ghci allows us to see type inference at work.

*Main> :l bb1.hs
[1 of 1] Compiling Main             ( bb1.hs, interpreted )
Ok, modules loaded: Main.
*Main> hiField field1
"Hi, id"
*Main> :type hiField
hiField :: JField -> [Char]

As we’ll see later in this series, Scala brings type inference to the JVM environment. Coming from the dynamic language side, the Diamondback Ruby research project is adding type inference to Ruby. So perhaps type inference is (finally) an idea whose time has come.

We’ll pick up more Haskell details along the way, but we have enough to start defining our first BuilderBuilder model. That will be the subject of the next post.

---

Updated 2009-05-09 to fix formatting.

BuilderBuilder: The Agenda

The “BuilderBuilder” posts will be my exploration of one question: “What does Functional Programming mean to me as a practicing OO software developer?” Despite contemporary buzz about how FP will save the world from bugs and single-threaded code, I don’t want to read (or write) another clever Fibonacci number generator, or another theoretical explanation of category theory and monads. I want to see what FP is like when it sits down the aisle from me in work clothes.

To pursue that goal, I’m going to explore a programming problem that meets these criteria:

1. Its purpose would be recognizable to a working programmer.
2. It is large enough that a clever one-liner won’t suffice.
3. It isn’t (consciously) chosen to cater to either OO or FP.
4. It will require me to try things outside my daily comfort zone.

That last point means that I’ll make mistakes and do stupid things. I hope the record, warts and all, helps someone else (including me, later on) avoid them.

The problem I’ve chosen is a specific code generation task, which I’ll describe in the next post. There is a wide range of opinions on source code generation, but source code is one thing that any programmer understands, regardless of other application domain specialization.

Because my goal is exploration, and I’m starting with unequal familiarity with the languages I’ll use, the code almost certainly won’t be idiomatic or make best use of all available libraries. In fact, I’ll probably do more than necessary “by hand” to avoid bias toward more familiar tools. I may be able to address some of that along the way, based on accumulated experience (or helpful feedback, hint, hint… .

OK, enough meta-talk. Next up, the task.

Scala and Programming 2.0

I commented elsewhere on how the “Architecture of Participation” idea may be percolating into the field of programming languages. I am especially interested in seeing whether the adoption of Scala provides evidence of this phenomenon.

Scala is a strongly, statically typed language implemented on the JVM—all characteristics that raise eyebrows (if not noses) in some circles. However, Scala’s ultralight approach to syntax is very much in line with the current taste for highly flexible notation and internal DSLs as a tool of expression.

Type inference is a very attractive compiler feature. And it’s great to get performance improvements “for free” every time the JVM team makes HotSpot smarter about JIT compilation, method in-lining, etc. But every time I revisit the ideas in Martin Odersky’s Scala talk at Javapolis 2007, I’m impressed with the design of Scala as a notation that invites participation.

Time will tell.

Project Euler 5: Multiplicity

The fifth problem from Project Euler asks for “the smallest number that is evenly divisible by all of the numbers from 1 to 20″. The key, for me, was to paraphrase that as “the least common multiple of 1 to 20″.

The least common multiple (lcm) of two natural numbers is their product divided by their greatest common divisor (gcd, a well-known lab rat). So, the stock definitions of lcm and gcd, along with one quick fold, are all we need.

  def divisibleByAll(lo: Int, hi: Int) = {
    def gcd(a: Int, b: Int) = gcd1(a max b, a min b)
    def gcd1(a: Int, b: Int): Int = if (b == 0) a else gcd(b, a % b)
    def lcm(a: Int, b: Int) = a * (b / gcd(a, b))
    (lo to hi).foldLeft(1)(lcm(_, _))
  }

There is one subtlety in this code; gcd is only there to insure that the (natural) arguments to gcd1 are correctly ordered (larger first). That is only needed once; if a >= b, then a % b cannot be larger than b, so the recursion in gcd1 is guaranteed to maintain the larger-first argument ordering.

Finally, this design encounters integer overflow at 1..24, changing the inner functions to work with Long gets us up to 1..42 before long overflow kicks in.

It should be obvious that the answers for 1..21 and 1..22 are the same as for 1..20, but I had to do a double-take while scanning for the limits. Duh.

Three X eerhT

The fourth task from Project Euler is to find the largest palindrome that is a product of two three-digit numbers. (I’m assuming decimal. )

Name no one man

The first sub-task is to determine whether the decimal representation of an integer is palindromic. A direct approach yields this:

  def reverseInt(n: Int) = {
    def rev(a: Int, b: Int): Int =
      if (a == 0)
        b
      else
        rev(a / 10, b * 10 + a % 10)
    rev(n, 0)
  }

  def isPalindrome(n: Int) = {
    reverseInt(n) == n
  }

I’m not approaching this by way of string conversion, for three reasons:

• Nothing in the problem statement involves strings;
• Expressing the solution in terms of string manipulation doesn’t make the solution significantly easier to understand; and
• Doing so would be slower.

Now on to the more interesting part.

Seek no monkees

The problem can be solved by searching a two-dimensional space, where each dimension ranges over the three-digit integers. A moment’s thought provides the following additional facts we can use:

• Multiplication is symmetric, so there’s no need to examine both products of two different numbers (a * b and b * a).
• The Scala expression 100 to 999 provides the three-digit integers.
• We can ignore 100 as a candidate, because any multiple of 100 must end in a zero, which isn’t the first digit of any multi-digit decimal number in standard notation.
• There’s a solution, because 101 × 101 = 10201, which is palindromic. (It’s nice to know that the search function doesn’t need to deal with an empty solution space.)

Going for the obvious, we can construct the products of non-redundant pairs of three-digit numbers, filter the products for palindromes, and capture the maximum of the palindromes:

  def maxOfPalindromes() = {
    (
      for {
        a <- 101 to 999
        b  a max b)
  }

It would be easy to parameterize this for the range to be searched (either by bounds or by number of digits), but I want to focus on the search performance.

Top spot

The version above takes about 152 ms on my laptop, which seems entirely too much time for this small a task. The following diagram shows the search space, as a context for the rest of the post.

The names in the diagram describe the state of the search at a given instant:

• lo and hi — the lower and upper bounds (inclusive) of the range being searched;
• a and b — the current position in the search space;
• c — the value at that position (the product of a and b); and
• x — the largest palindrome found thus far.

In his landmark book, A Discipline of Programming (published 32 years ago!), Edsger Dijkstra stated the “Linear Search Theorem”. I’ll paraphrase it informally for our purposes as:

When searching for the maximum value having some property, examine the values in decreasing order; the first satisfactory value found will be the maximum.

Assuming that we apply this to both a and b, the green area has been examined and the blue area has not. We can express the next state of the search in a decision table (each row’s condition assumes that conditions on all previous rows failed):

Condition	Action	State changes
		next a	next b	next x
a < lo	Search complete x is the result	n/a
b < a	Current row searched, move to next blue row	a - 1	hi	x
c <= x	c and rest of this row eliminated, move to next blue row	a - 1	hi	x
isPalindrome(c)	c is new max palindrome (rest of this row is smaller values)	a - 1	hi	c
otherwise	c eliminated, continue searching current row	a	b - 1	x

The decision table translates directly to a tail-recursive function. I’ll use a bit of scope management to avoid computing c unless it will be used.

  def maxPalindromeIn(lo: Int, hi: Int) = {
    def loop(a: Int, b: Int, x: Int): Int = {
      if (a < lo)
        x
      else if (b < a)
        loop(a - 1, hi, x)
      else {
        val c = a * b
        if (c <= x)
          loop(a - 1, hi, x)
        else if (isPalindrome(c))
          loop(a - 1, hi, c)
        else
          loop(a , b - 1, x)
      }
    }
    loop(hi, hi, 0)
  }

  def maxPalindrome3x3() = maxPalindromeIn(101, 999)

This runs in about 445µs, a speedup factor of well over 300.

I prefer π

Let me round this post off with a few observations.

Performance is not the most important issue. I probably need to state that clearly, having quoted run time so often in this and recent Project Euler posts. I believe that speed, small memory footprint, and any other measure of economy, should be subordinate and subsequent to good design and correctness.

Those who forget history are condemned to repeat it. Over twenty years ago, when the dreaded “G” word was still occasionally heard in polite programming circles, a misguided letter to the editor of the Communications of the ACM attempted to claim that goto was necessary to effective programming and offered a small problem in defense of that claim. Dijkstra’s response addressed the errors and approach of the claim and its “solution”; his demonstration used the idea of the two-dimensional application of the Linear Search theorem.

There are important ideas that transcend programming model. The ideas in A Discipline of Programming are all presented in an imperative style, yet the clarity of design makes the concepts important to FP as well. That book is not an easy read (I worked rather than read through it when I got my copy in 1976) but it rewards the investment richly. David Gries gave a tutorial-style presentation of many of the same core ideas in The Science of Programming. I cannot possibly recommend them highly enough.

Compare and contrast. It’s useful to see how others approach a problem, after you’ve solved it yourself.

• The Haskell web site has Haskell solutions to the Project Euler problems; another site offers Python-based solutions of the first several. However, these sites have very little explanation of the solution or design strategy (my focus).
• You can register at Project Euler to track your own progress, and read the comments posted by others when you’ve completed a problem. (Comments are closed for the earlier problems.)
• Russ Gray (Basildon Coder) is working through the problems in a variety of languages and posting on his experiences. His discussion of the first problem is what tipped me over the edge into this series.
• Dianne Marsh (a friend from Java Posse Roundup) is also doing Project Euler problems in Scala. Apparently she’s roped in collegues from SRT: Bill Wagner in C#, Darrell Hawley in Python, Marina Fedner in Ruby, and Anne Marsan in Matlab.
• I’m also aware of a partial set of solutions in F# by Chris Smith

My RSS reader is subscribed to the fellow solvers above; many of them have other interesting things to say beyond Project Euler!

---

Recommended reading:

• A Discipline of Programming (Prentice-Hall Series in Automatic Computation)
• The Science of Programming (Monographs in Computer Science)

A composite function

The previous post began working with problem 3 from Project Euler. The allFactors function produced a list of the prime factors of its argument (all occurrences of all prime factors), such as:

  scala> val nums = allFactors(96000)
  nums: List[Long] = List(5, 5, 5, 3, 2, 2, 2, 2, 2, 2, 2, 2)

It would be nice to have a representation that’s a bit closer to the conventional 5³ × 3¹ × 2⁸, which makes clear the exponent for each prime factor.

In thinking about solutions to that issue, we get a prime use case for Scala’s composite nature as a hybrid OO/FP language. In addition, we encounter another heuristic for avoiding errors (functional or otherwise).

A separate piece

One easy solution is to write a separate function that collapses a list of values into a list of value–count pairs, for which we can define a type alias as:

  type PrimePower = (Long, Int)

The partition method on List is perfect for our current need; given a condition (represented as a boolean-valued function), it splits the list into the prefix that satisfies the condition and the remainder.

  scala> nums.partition(_ == 5)
  res82: (List[Long], List[Long]) = (List(5, 5, 5),List(3, 2, 2, 2, 2, 2, 2, 2, 2))

Of course, the length of the prefix is the exponent, so the function almost writes itself:

  def groupCount(fs: List[Long]): List[PrimePower] = fs match {
    case Nil => Nil
    case p :: _ =>
      val (ps, rest) = fs.partition(_ == p)
      (p, ps.length) :: groupCount(rest)
  }

If factor list is empty, the result is an empty list. On the other hand, if fs begins with some prime p, we split the prefix containing all ps from the rest of the list, and return a list containing the PrimePower for p consed onto the groupCounted rest of the original list.

  scala> groupCount(nums)
  res86: List[()] = List((5,3), (3,1), (2,8))

The advantages of writing this as a separate function are ease of testing and availability for re-use. On the other hand, it would be interesting to see the consequences of directly producing the prime powers.

One list to bind them

For reference, here’s the function from last time, with the inner alternatives labeled for reference and highlighting on the parts we need to consider:

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

The result will be List[PrimePower], which is the new type of the second argument and result of loop.

Case 2 deals with a quotient which we know to be prime, because all possible divisors have been exhausted. Given that q occurs exactly once, we replace q :: fs with (q, 1) :: fs as the result.

Case 3 requires the most thought. It only accounts for a single occurrence of a divisor; q / d removes one factor of d from q and d :: fs saves that occurrence in the list being accumulated.

Help me, Rhonda!

We could use a helper function to handle d, based on whether it is another occurrence of the previous divisor or a new divisor not previously seen. Because that distinction is based on data in the list, it seems natural to drive the logic by matching on the list:

  def countFactor(f: Long, fs: List[PrimePower]) = fs match {
    case (d, c) :: rest if f == d => (d, c + 1) :: rest
    case _                        => (f, 1) :: fs
  }

To test countFactor, we need to iterate it over a sequence of values for f with fs initially Nil. This looks like a job for a fold!

  def testCountFactor(ls: List[Long]) =
    ls.foldRight(Nil: List[PrimePower])(countFactor(_, _))

In the interest of space, I won’t include all the test scenarios: Nil, a list of length 1, all values equal, all values different, etc. As the old saying goes, checking those is “left as an exercise for the reader”.

True confessions

Substituting a call to countFactor into case 3, along with the other changes already discussed, gives us this…

  def allFactorsCounted(t: Long) = {
    def loop(q: Long, fs: List[PrimePower], d: Long): List[PrimePower] =
      if (q == 1)
        fs                                    // case 1
      else if (q < d * d)
        (q, 1) :: fs                          // case 2
      else if (q % d == 0)
        loop(q / d, countFactor(d, fs), d)    // case 3
      else
        loop(q, fs, d + 1)                    // case 4
    loop(t, Nil, 2)
  }

…which behaves as follows…

  scala> val nums2 = allFactorsCounted(96000)
  nums2: List[()] = List((5,1), (5,2), (3,1), (2,8))

…or, I should say, “misbehaves as follows”! What went wrong?

The root cause of this error is the special-case code in case 2, which is doing two things at once: handling a prime factor (in its own way!) and terminating the recursion. Put another way, the design breaks the symmetry between cases 2 and 3, both of which identify a prime factor and do something with it. We can make that symmetry more obvious by rewriting case 2 as in the following:

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

And now each case in the inner loop has exactly one responsibility:

1. Terminates the computation because all factors have been found;
2. Removes q as a known factor;
3. Removes d as a known factor;
4. Moves past an unsuccessful divisor.

Mixing concerns in a single part of the code creates these risks:

• The code can become harder to understand;
• Symmetries can be hidden or broken;
• Change driven by (only) one concern can have unintended consequences.

For me, the moral of this error is this:

Each part of the code (especially alternatives within a single function) should “Do one thing, and do it well!”

An object lesson

In the OO world, primitive obsession is regarded as a code smell. Notice that q and fs together represent the current state of our factorization; we have the invariant that q times the product of fs equals the original number to be factored! Failing to maintain that invariant would likely break the code. The clue to watch for is that the (tail-)recursive calls of the inner loop function either alter both of those parameters or neither of them. This is telling us that they are related aspects of a single concept.

OO gives us a technique for managing related data: put them in an object as private fields, expose the minimum methods required by the object’s client, and ensure that those methods maintain the invariants. We can refactor the quotient and factor list into the following class:

  class Factorization(val n: Long, val fs: List[PrimePower]) {
    def isComplete = n == 1
    def move(d: Long): Factorization = {
      def remove(i: Int, q: Long, d: Long): Factorization =
        if (q % d == 0)
          remove(i + 1, q / d, d)
        else
          new Factorization(q, (d, i) :: fs)
      if (n % d == 0)
        remove(0, n, d)
      else
        this
    }
  }

Each instance of Factorization represents some stage of the factorization of a number by n (the non-factored part) and fs (the factors previously found). An instance can tell whether it is a complete factorization (all factors have been extracted). Finally, an instance can return a Factorization resulting from move-ing all occurrences of a given divisor from the non-factored part to the factor list.

Notice that a Factorization is immutable; the result of a move will be the same instance if the divisor isn’t a factor, or a new instance resulting from the extraction of that divisor.

With this class defined, our factorization function is even more conceptually streamlined:

  def factorize(n: Long): List[PrimePower] = {
    def loop(f: Factorization, d: Long): Factorization = {
      if (f.isComplete)
        f
      else if (f.n < d * d)
        loop(f.move(f.n), d)
      else
        loop(f.move(d), d + 1)
    }
    loop(new Factorization(n, Nil), 2).fs
  }

If the factorization is complete, then its factor list is the desired result. Otherwise, ask the factorization to remove a candidate divisor and continue with the result of that request. The inner loop function is only concerned with managing the sequence of divisors until the factorization is complete.

So the second moral of this post is this:

Whenever part of a function’s parameter list represents some concept, consider making that explicit by using a class for the concept.

Doing so will allow that concept to be implemented and tested independently, and will improve the clarity of the remaining code. Although the example in this post is a very small programming task, I think it makes a nice demonstration of the value of Scala’s hybrid nature.

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)
        f
      else if (q < d * d)
        q
      else if (q % d == 0)
        loop(q / d, d, d)
      else
        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)
        factor
      else if (quotient < divisor * divisor)
        quotient
      else if (quotient % divisor == 0)
        loop(quotient / divisor, divisor, divisor)
      else
        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)
        fs
      else if (q < d * d)
        q :: fs
      else if (q % d == 0)
        loop(q / d, d :: fs, d)
      else
        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
factor type	Long	Long	Long
accumulator type	Long	List[Long]	T
"zero" value	1	Nil	?
composition type	Long → Long → Long	Long → List[Long] → List[Long]	Long → T → T
composition function	replacement (a, b) → a	cons (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)
        acc
      else if (q < d * d)
        f(q, acc)
      else if (q % d == 0)
        loop(q / d, f(d, acc), d)
      else
        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 5² × 3³ × 2⁷ 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.

Go with the flow

The second problem from Project Euler asks for the sum of the even Fibonacci numbers which are at most 4,000,000.

An aside on the numbers

The problem’s description of the Fibonacci sequence is flawed (or at least non-standard); the problem statement begins the sequence with 1 and 2, showing the first 10 terms as:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89

It is common to see the definition given as:

fib₁ = 1, fib₂ = 1, fib_n = fib_n-1 + fib_n-2

I’m persuaded by Dijkstra’s line of reasoning, and prefer to use the definition…

fib₀ = 0, fib₁ = 1, fib_n = fib_n-1 + fib_n-2

…not only for the benefit of a zero origin, but also because it invites the mind to consider completing the story. The equations…

fib_n = fib_n-1 + fib_n-2

…and…

fib_n – fib_n-1 = fib_n-2

…are clearly equivalent, but the second one may help us realize that we can extend the conventional list to be infinite in both directions…

n	…	-7	-6	-5	-4	-3	-2	-1	0	1	2	3	4	5	6	7	…
fibn	…	13	-8	5	-3	2	-1	1	0	1	1	2	3	5	8	13	…

…making fib a total function over the integers. I’ll put off discussing that interesting option until another time, and for the rest of this post use only natural numbers.

Basic solutions

The naive recursive function for the n^th Fibonacci number…

  def fibR(n: Int): Int = {
    if (n == 0)
      0
    else if (n == 1)
      1
    else
      fibR(n - 1) + fibR(n - 2)
  }

…requires O(exp n) time to evaluate. We could go for a tail-recursive (iterative) version…

  def fibI(n: Int): Int = {
    def fib0(i: Int, a: Int, b: Int): Int =
      if (i == n) a else fib0(i + 1, b, a + b)
    fib0(0, 0, 1)
  }

…which executes in O(n) time. That’s an improvement, but still leaves us dealing with the Fibonacci sequence one number at a time.

Gently down the stream

The problem calls for a sum across multiple values from the Fibonacci sequence, and both of the above functions calculate “upwards” to reach the specified value. There’s no point in going through the sequence more than once (regardless of how efficiently we can do that), so let’s think of the sequence as a Stream[Int].

Element-wise addition of streams can be done as:

  def add(s1: Stream[Int], s2: Stream[Int]): Stream[Int] =
       Stream.cons(s1.head + s2.head, plus(s1.tail, s2.tail))

A stream of the Fibonacci numbers is then:

  val fibS: Stream[Int] = Stream.cons(0, Stream.cons(1, add(fibS, fibS tail)))

And now we’re back to familiar territory from the first problem. We simply filter for the even values, apply the limit, and sum the results:

  scala> fibS.filter(_ % 2 == 0).takeWhile(_ <= 4000000).foldLeft(0)(_ + _)
  res17: Int = 4613732

We could stop with that solution, but there are a couple of additional techniques worth a mention.

“Go back, Jack, and do it again”

First, we could replace the stream with a simple Iterator[Int] that holds a pair of consecutive Fibonacci numbers:

  class Fibs extends Iterator[Int] {
    var state = (0, 1)
    def hasNext = true
    def next = {
      val (a, b) = state
      state = (b, a + b)
      a
    }
  }

Thanks to the Iterator trait, when we implement hasNext (is there another value?) and next (return the next value), we get everything needed to complete the solution:

  scala> new Fibs().filter(_ % 2 == 0).takeWhile(_ <= 4000000).foldLeft(0)(_ + _)
  res20: Int = 4613732

However, this still generates three times as many values as needed, because only every third Fibonacci number is even:

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144…

How could we generate only the even Fibonacci numbers? Notice this progression in the numbers above:

fib₀ = 0
fib₃ = 2
fib₆ = 8
fib₉ = 34
fib₁₂ = 144
…

So we’re interested in producing the sequence fib_3n rather than fib_n if we want to obtain the even values. In our Fib class above, a and b are consecutive values, fib_n and fib_n+1. From the original definition, we can apply just a little algebra:

fib_n+2 = b + a
fib_n+3 = fib_n+2 + fib_n+1 = (b + a) + b = 2 × b + a
fib_n+4 = fib_n+3 + fib_n+2 = (2 × b + a) + (b + a) = 3 × b + 2 × a

If n is a multiple of 3, then (n + 3) is the next multiple of 3. So the above formulas for fib_n+3 and fib_n+4 give us a way to write an iterator that produces only even Fibonacci numbers:

  class EvenFibs extends Iterator[Int] {
    var state = (0, 1)
    def hasNext = true
    def next = {
      val (a, b) = state
      state = (2 * b + a, 3 * b + 2 * a)
      a
    }
  }

This new iterator eliminates the need for the filter, reducing the solution to:

  scala> new EvenFibs().takeWhile(_ <= 4000000).foldLeft(0)(_ + _)
  res30: Int = 4613732

This calculates exactly the values of interest, then uses standard methods to provide the sum. These last three versions retain the ability to re-use the parts that generate the Fibonacci numbers (or even Fibonacci numbers) for purposes other than the sum.

For sake of completeness, here’s the completely optimized version…

  def sumEvenFibs(n: Int) = {
    def iter(a: Int, b: Int, s: Int): Int =
      if (a > n) s else iter(2 * b + a, 3 * b + 2 * a, a + s)
    iter(0, 1, 0)
  }

…which buys performance by giving up that reusability.