Exploring SML

Last summer, one of the books which silently judged me from my shelf, collecting dust, was Elements of ML Programming by Jeffrey Ullman. I avoided reading a lot of CS stuff then—mainly because I was programming at work and didn't feel like doing more of the same in my free time. But during the fall semester I was back to a life of stress and pain carefree leisure. For fifteen weeks, I had an appetite for CS material again since most of my classes were mathematics (sorry Math Department, I still love you too).

And then I went back to work (and school!) spring of this year. This post has taken longer than I thought.

So, I finally cracked it open to get my programming fix, and oh boy, is it a good one! You can consider this post half book review and half language review.

You Should Own This Book

The Lord of the Rings, The Phantom Tollbooth, and Elements of ML Programming all have one feature in common. Cool maps. Ullman one-upped Tolkein and Norton Juster, too; Elements' map is right on the front cover.

There are nine chapters, and they cover everything you might need to get started with the language. It doesn't assume prior experience with functional programming, but it is written for someone who already knows a common and frequently-taught language like, uh... Pascal. (It was written in 1994, OK?)

I was impressed by how much useful material there is that's not SML-specific. There are tons of worked examples of complexity analysis for various functions, and discussions on cool algorithms like the ones we'll see below.

What's SML?

Standard Meta-Language was created by Robin Milner at the University of Edinburgh. When a language infers the datatype of a value, it's probably using a Hindley-Milner type system (like SML's) to do it. SML is strongly typed, which means you can determine the types of all expressions without running the program. It's also functional, so it usually computes without side effects. There's lots of recursion. Like, lots.

Here's a simple example:

fun gcd (m, 0) = m
  | gcd (m, n) = gcd (n, m mod n)

It's Euclid's GCD algorithm in just two lines. How cool is that?

The fun keyword means this is a fun language to use a function, of course, but what happens after would've been strange to a Pascal programmer. The SML environment matches the pattern of the input. If the input is a tuple of some integer m and 0, then the GCD is m. Any other input matches the second line, (m, n), and the function recurses on that. Also notice that the keyword mod is used for modulo, not % or rem. (Actually rem is still a keyword, but it does something different.)

Here's another little example, which checks whether a number is prime.

fun divisible (_, 0) = false
  | divisible (_, 1) = false
  | divisible (n, i) = divisible (n, i - 1) orelse n mod i = 0

fun isPrime 2 = true
  | isPrime n = not (divisible (n, ceil(Math.sqrt(real n))))

There's a helper function, divisible, which uses something in its pattern-matching we haven't seen. The _ is a wildcard, and matches any pattern. It also indicates we're throwing away the value in that case. If it wasn't weird enough already, SML uses orelse for logical OR, and andalso for logical AND.

That code isn't pretty (my first SML code so don't judge). I arrived at a cleaner implementation with some thought.

fun divisible (p, q, l) =
    q <= l andalso (p mod q = 0 orelse divisible (p, q + 1, l))

fun isPrime 2 = true
  | isPrime n = n > 1 andalso not (divisible (n, 2, Real.floor (Math.sqrt (Real.fromInt n))))

Now all the recursion is wrapped up in the divisible function; the primality check never calls itself.

Fibonacci Sequence

We can't talk about recursion without first talking about recursion without mentioning the Fibonacci sequence. You can guess how the recursive implementation looks.

fun fib 0 = 0
  | fib 1 = 1
  | fib n = fib (n - 1) + fib (n - 2)

You can also guess this implementation makes cold molasses look like Usain Bolt. Here's the (more verbose) memoized implementation:

(* Memoized version *)
fun mfib n =
  let
    open IntRedBlackMap
    
    fun mfib' (0, memo) = (0, memo)
      | mfib' (1, memo) = (1, memo)
      | mfib' (n, memo) =
        case find (memo, n) of
          SOME v => (v, memo)  (* Found it; return v and memo *)
        | NONE =>
          let
            val (v1, memo1) = mfib' (n - 1, memo)
            val (v2, memo2) = mfib' (n - 2, memo1)
            val result = v1 + v2
            val new_memo = insert (memo2, n, result)
          in
            (result, new_memo)
          end
  in
    #1 (mfib' (n, empty)) (* Start with an empty map *)
  end

Syntactically, this looks a lot different. Right at the top you see a let block, which lets you write out variables in some expression, e.g. let <vals> in <expr> end. Then comes an open statement, which basically brings the data structure IntRedBlackMap and its functions into scope.

Then we have a function declaration inside a function declaration. We could even rewrite the primality test above to declare the divisible function this way. Then it wouldn't be visible outside that internal scope. You can also see where Rust's Option<T> type came from in that case statement!

The purely recursive version had to recompute all the previous Fibonacci numbers on every call. Here the mfib' function is caching all previously-computed results for reuse, making this implementation lots faster. Finally the original mfib function returns the expression #1 (mfib' (...)), which grabs the first element of the tuple returned by mfib'.

This is a little more verbose than it has to be, but it shows off the idea of tuples, pattern matching, and a bit of the standard library. Next, I'll bring in a code example from the book itself.

Karatsuba Multiplication

We haven't dealt with lists in SML yet, even though lists are probably the most important data structures in functional languages. I liked Case Study 3.6, which builds up to polynomial multiplication via the Karatsuba algorithm. We can represent polynomials as lists of coefficients, and SML has terrific support for list operations.

One of the basic operations which we'll use is polynomial addition. For some polynomials \[c(x)=c_0+c_1x+c_2x^2+\cdots+c_nx^n \quad\mathrm{and}\quad d(x)=d_0+d_1x+d_2x^2+\cdots+d_nx^n\] the sum is just \[(c + d)(x)=(c_0+d_0)+(c_1+d_1)x+(c_2+d_2)x^2+\cdots +(c_n+d_n)x^n.\] If our polynomials are two lists of coefficients, we just add the coefficients element-wise. Here's the recursive implementation:

(* polynomial addition of P and Q *)
fun padd (P, nil) = P
  | padd (nil, Q) = Q
  | padd ((p:real)::ps, q::qs) = (p + q)::padd (ps, qs)

The keyword nil matches the empty list, as in LISP. Line 3 has the recursive case, where it matches the lists to the patterns p::ps and q::qs. The :: is the "cons" operator, or list concatenation. In a pattern it splits the head of the list off and stores it in p or q, and puts the tail in ps or qs.

Another basic operation is scalar multiplication (three guesses on how that works). The code is similar:

(* scalar multiplication of P by q *)
fun smult (nil, _) = nil
  | smult ((p:real)::ps, q) = (p * q)::smult(ps, q)

These two operations are enough to give us a complete polynomial multiplication algorithm in a two-line function.

(* polynomial multiplication of P and Q *)
fun pmult (_, nil) = nil
  | pmult (P, q::qs) = padd(smult(P, q), 0.0::pmult (P, qs))

It's not immediately obvious (to me!) that this does what I just claimed. Ullman proves the algorithm works with induction. The following comes almost verbatim from the book:

Basis: If the second polynomial is empty, the result is empty.
Induction: If the second polynomial \(Q\) can be written as \(q+Sx\), then \[PQ=Pq+PSx.\] The product \(Pq\) is a scalar multiplication. \(PS\) is a recursive application of this algorithm with a smaller second argument. The only thing that really needs explaining is the subexpression 0.0::pmult(P, qs), which multiplies by \(x\) (i.e. shifts right one by inserting zero). ∎

He then walks through the analysis to show that the running time is in \(\mathcal{O}(n^2)\) assuming the two polynomials have length \(n\) (as you'd probably expect it to be!)

Multiplying Faster
As it turns out, there is a faster algorithm in \(\mathcal{O}(n^{\mathrm{lg}\,3})\) (where \(\mathrm{lg}(x)\) is log base 2; that's Donald Knuth's notation). Ullman points out that you can do multiplication in \(n\,\mathrm{lg}\,n\) time with an FFT, but we're not going to write an FFT routine in SML.

We can split each polynomial into two smaller polynomials. For polynomials \(P\) and \(Q\), splitting at the \(s\)th term we'd get \[P=T+x^sU\quad\mathrm{and}\quad Q=V+x^sW.\] So the product becomes \[PQ=TV+x^s(TW+UV)+x^{2s}UW.\] This makes life easier. Even though we now have four polynomial multiplication problems (\(TV\), \(TW\), \(UV\), & \(UW\)), they're all half-sized, and we just need some extra linear-time shifts and additions.

But multiplication is still a quadratic-time operation. The trick that the Karatsuba-Ofman algorithm uses is to eliminate one multiplication by rewriting the middle part \(TW+UV\) as \[TW+UV=(T+U)(V+W)-TV-UW.\] Notice the products \(TV\) and \(UW\) we're subtracting were already computed anyway, so we only need the one multiplication \((T+U)(V+W)\) instead of the two before. Now if \(n\) is large enough, our algorithm will beat out the quadratic-time one, because this function calls itself only three times, where it would have called itself four times without the trick. Ullman walks through the complexity analysis. Buy the book if you want to read the proof.

Now, the code:

(* length computes degree + 1 of a polynomial *)
fun length nil = 0
  | length (_::ps) = 1 + length ps

(* psub computes P - Q *)
fun psub (P, Q) = padd (P, smult (Q, ~1.0) )

(* bestSplit computes the optimal size for the
   low-order half of polynomials with lengths m &
   n. It's min{m, n} if one is less than half the
   other, otherwise, it's half the larger. *)
fun bestSplit (m, n) =
    if 2 * n <= m then n
    else if 2 * m <= n then m
    else if n <= m then m div 2
    else n div 2

(* shift computes x^n times a polynomial *)
fun shift (P, 0) = P
  | shift (P, n) = 0.0::shift (P, n - 1)

(* carve computes two polynomials where the first
   is the n low-order terms of P, and the second
   are the remaining terms divided by x^n. *)
fun carve (P, 0) = (nil, P)
  | carve (nil, _) = (nil, nil)
  | carve (p::ps, n) =
    let
      val (qs, rs) = carve (ps, n - 1)
    in
      (p::qs, rs)
    end

(* komult computes the product of two polynomials
   using the Karatsuba-Ofman algorithm which runs
   in O(n^lg 3) time. *)
fun komult (_, nil) = nil
  | komult (nil, _) = nil
  | komult (P, [q]) = smult (P, q)
  | komult ([p], Q) = smult (Q, p)
  | komult (P, Q) =
    let
      val n = length P
      val m = length Q
      val s = bestSplit (n, m)
      val (T, U) = carve (P, s)
      val (V, W) = carve (Q, s)
      val TV = komult (T, V)
      val UW = komult (U, W)
      val TUVW = komult (padd (T, U), padd (V, W))
      val middle = psub (psub (TUVW, TV), UW)
    in
      padd (padd (TV, shift (middle, s)), shift (UW, 2 * s))
    end

At this point, there's no new syntax in this code we haven't seen before, besides the fact that you have to use '~' for negating a real number. I'm not going to explain the entire program. One of the exercises is to write this function...

(* Generate a polynomial with all coefficients 1.0 for test purposes. *)
fun genPoly 0 = nil
  | genPoly n = 1.0 :: genPoly (n - 1)

...and see at what point komult is faster than pmult. It's visibly noticable at \(n=1000\) on my machine.

Back to the Present

I used SML a couple days ago for a simple problem in my number theory class. The question was to find the least "rectangular number" which is also "triangular." You can make a rectangle with a rectangular number of dots, and you can make a triangle with a triangular number of dots. For example, \(6\). In fact, \(6\) is the answer (spoiler).

. . .
. . .

and

  .
 . .
. . .

I wanted to do a search with a computer to see what other numbers like this exist. Let's think about the definitions.

Mathematically, a rectangular number is the product of two consecutive positive integers, so like \(n=a(a+1)\), and a triangular number is the sum of the first however-many integers, like \(6=1+2+3\) above, or in general, \(n=\frac{b(b+1)}{2}\). You get a rectangular-triangular number when these are equal to each other, so we can set the equations equal and solve for \(b\) in terms of \(a\). \[ \begin{aligned} a(a+1) &= \frac{1}{2}b(b+1) \\ b^2 + b - 2a(a+1) &= 0 \end{aligned} \] And solving for \(b\) via the quadratic formula, we get \[ \begin{aligned} b &= \frac{-1 \pm \sqrt{1 - 4(1)(-2a(a+1))}}{2} \\ &= \frac{-1 \pm \sqrt{1 + 8a(a+1)}}{2}. \end{aligned} \] You can see that \(b\) is an integer only when the discriminant \(1 + 8a(a+1)\) is a perfect square. Since it's of the form \(8k + 1\) it must be an odd perfect square. So if we search for all \(a\) such that \(8a^2 + 8a + 1\) is a perfect square, then we will find a pair \((a, b)\) satisfying our original equation.

If we want to search at large integer values, we should probably use arbitrary-precision integer types. SML calls these IntInf, and they're in the standard library (standard basis in ML-speak). You'll also notice I nested a function inside a function inside a function. That was fun. :)

(* Search for integer solutions to the equation a(a+1) = 0.5 * b(b+1). *)
fun searchSolns (maxA : IntInf.int) =
    (* To search efficiently, find `a` such that
     *   8 * a * (a + 1) + 1
     * is an odd perfect square, in which case `b` is
     *   (sqrt(a) - 1) / 2
     * and (a, b) are solutions to
     *   a(a+1) = 0.5 * b(b+1).
     *)
    let
        (* Find the integer floored square root of n via the Babylonian method. *)
        fun isqrt (n : IntInf.int) : IntInf.int =
            if n < 0 then raise Fail "Negative inputs forbidden!"
            else if n = 0 then 0
            else
                let
                    fun loop x =
                        let
                            val next = (x + n div x) div 2
                        in
                            if next >= x then x
                            else loop next
                        end
                in
                    loop (IntInf.max (1, n div 2))
                end

        (* For some a, check whether it can form a triangular-rectangular number with some b. 
         * Then cons that pair with acc until a > maxA.
         *)
        fun check (a, acc) =
            if a > maxA then
                List.rev acc
            else
                let
                    val target =  8 * a * (a + 1) + 1
                    val s = isqrt target
                in
                    if s * s = target then
                        let
                            val b = (s - 1) div 2
                        in
                            check (a + 1, (a, b, a * (a + 1)) :: acc)
                        end
                    else
                        check (a + 1, acc)
                end
    in
        check (1, [])
    end

And there we go. Here's a run of the program:

ethan@hostname $ sml search.sml
Standard ML of New Jersey [Version 110.99.9; 64-bit; November 4, 2025]
[opening pell.sml]
[autoloading]
[library $SMLNJ-BASIS/basis.cm is stable]
[library $SMLNJ-BASIS/(basis.cm):basis-common.cm is stable]
[autoloading done]
val searchSolns = fn : IntInf.int -> (?.intinf * ?.intinf * ?.intinf) list
- (* Loading done, run the function. *)
- searchSolns 1000000;
val it =
  [(2,3,6),(14,20,210),(84,119,7140),(492,696,242556),(2870,4059,8239770),
   (16730,23660,279909630),(97512,137903,9508687656),
   (568344,803760,323015470680)] : (?.intinf * ?.intinf * ?.intinf) list
- 

The output is a list of triples \((a, b, n)\) where \(n\) is the rectangular-triangular number. This sequence of \(n\) is A029549 in the OEIS. Some remarks from that page:

Triangular numbers that are twice other triangular numbers. - Don N. Page
Triangular numbers that are also pronic numbers. These will be shown to have a Pythagorean connection in a paper in preparation. - Stuart M. Ellerstein (ellerstein(AT)aol.com), Mar 09 2002
In other words, triangular numbers which are products of two consecutive numbers. E.g., a(2) = 210: 210 is a triangular number which is the product of two consecutive numbers: 14 * 15. - Shyam Sunder Gupta, Oct 26 2002

So that's interesting!

Conclusion

The Language Review

I really enjoy SML now, and I'll keep using it for things in my mathematical career. It's now my favorite functional language (I don't think Rust strictly counts as functional) and the recursive way of doing things makes proofs really easy; just do induction. There's so much more functionality in the language than I covered in this post. I'll be looking for more excuses to use it across the board.

The Book Review

Jeffrey Ullman is a great author, that much is obvious. Now I want more of his books. I got this one from Dr. Stephen Rainwater (my academic hero, thanks Dr. Rainwater!), since he retired last year. I'm not sure what else to say. The writing is clear and concise, the typography is gorgeous, the examples are deep, well-thought-out, and would be useful no matter what language you want to work in.