16
$\begingroup$

I'm familiar with the idea of basic tail recursion elimination, where functions that return the direct result of a call to themselves can be rewritten as iterative loops.

foo(...):
    # ...
    return foo(...)

I also understand that, as a special case, the function can still be rewritten if the recursive call is wrapped in a call to cons.

foo(...):
    # ...
    return (..., foo(...))

What property of cons allows this? What functions other than cons can wrap a recursive tail call without destroying our ability to rewrite it iteratively?

GCC (but not Clang) is able to optimize this example of "tail recursion modulo multiplication," but it's unclear what mechanism allows it to discover this or how it makes its transformations.

pow(x, n):
    if n == 0: return 1
    else if n == 1: return x
    else: return x * pow(x, n-1)
$\endgroup$
3
  • 1
    $\begingroup$ In your Godbolt compiler explorer link, your function has if(n==0) return 0; (not return 1 like in your question). x^0 = 1, so that's a bug. Not that it matters for the rest of the question, though; the iterative asm checks for that special case first. But strangely, the iterative implementation introduces a multiply of 1 * x that wasn't present in the source, even if we make a float version. gcc.godbolt.org/z/eqwine (and gcc only succeeds with -ffast-math.) $\endgroup$ Commented Mar 16, 2019 at 10:01
  • $\begingroup$ @PeterCordes Good catch. The return 0 has been fixed. The multiplication by 1 is interesting. I'm not sure what to make of it. $\endgroup$
    – Maxpm
    Commented Mar 16, 2019 at 18:06
  • $\begingroup$ I think it's a side-effect of the way GCC transforms when turning it into a loop. Clearly gcc has some missed optimizations here, e.g. missing it for float without -ffast-math, even though it's the same value being multiplied every time. (Except for the 1.0f` which might be the sticking point?) $\endgroup$ Commented Mar 16, 2019 at 18:11

2 Answers 2

15
$\begingroup$

While GCC likely uses ad-hoc rules, you can derive them in the following way. I'll use pow to illustrate since you're foo is so vaguely defined. Also, foo might best be understood as an instance of last-call optimization with respect to single-assignment variables as the language Oz has and as discussed in Concepts, Techniques, and Models of Computer Programming. The benefit of using single-assignment variables is it allows remaining within a declarative programming paradigm. Essentially, you can have each field of the struct foo returns represented by single-assignment variables that are then passed to foo as additional arguments. foo then becomes a tail-recursive void returning function. No particular cleverness is needed for this.

Returning to pow, first, transform into continuation passing style. pow becomes:

pow(x, n):
    return pow2(x, n, x => x)

pow2(x, n, k):
    if n == 0: return k(1)
    else if n == 1: return k(x)
    else: return pow2(x, n-1, y => k(x*y))

All calls are tail calls now. However, the control stack has been moved into the captured environments in the closures representing the continuations.

Next, defunctionalize the continuations. Since there is only one recursive call, the resulting data structure representing the defunctionalized continuations is a list. We get:

pow(x, n):
    return pow2(x, n, Nil)

pow2(x, n, k):
    if n == 0: return applyPow(k, 1)
    else if n == 1: return applyPow(k, x)
    else: return pow2(x, n-1, Cons(x, k))

applyPow(k, acc):
    match k with:
        case Nil: return acc
        case Cons(x, k):
            return applyPow(k, x*acc)

What applyPow(k, acc) does is take a list, i.e. free monoid, like k=Cons(x, Cons(x, Cons(x, Nil))) and make it into x*(x*(x*acc)). But since * is associative and generally forms a monoid with unit 1, we can reassociate this into ((x*x)*x)*acc, and, for simplicity, tack a 1 on to start, producing (((1*x)*x)*x)*acc. The key thing is that we can actually partially compute the result even before we have acc. That means instead of passing k around as a list which is essentially some incomplete "syntax" that we'll "interpret" at the end, we can "interpret" it as we go. The upshot is that we can replace Nil with the unit of the monoid, 1 in this case, and Cons with the operation of the monoid, *, and now k represents the "running product". applyPow(k, acc) then becomes just k*acc which we can inline back into pow2 and simplify producing:

pow(x, n):
    return pow2(x, n, 1)

pow2(x, n, k):
    if n == 0: return k
    else if n == 1: return k*x
    else: return pow2(x, n-1, k*x)

A tail-recursive, accumulator-passing style version of the original pow.

Of course, I'm not saying GCC does all this reasoning at compile-time. I don't know what logic GCC uses. My point is simply having done this reasoning once, it's relatively easy to recognize the pattern and immediately translate the original source code into this final form. However, the CPS transform and defunctionalization transform are completely general and mechanical. From there fusion, deforestation, or supercompilation techniques could be used to attempt to eliminate the reified continuations. The speculative transformations could be thrown away if it is not possible to eliminate all the allocation of the reified continuations. I suspect, though, that that would be too expensive to do all the time, in full generality, hence more ad-hoc approaches.

If you want to get ridiculous, you can check out the paper Recycling Continuations which also uses CPS and representations of continuations as data, but does something similar to but different to tail-recursion-modulo-cons. This describes how you might produce pointer-reversing algorithms by transformation.

This pattern of CPS transforming and defunctionalizing is a quite powerful tool for understanding, and is used to good effect in a series of papers I list here.

$\endgroup$
3
  • $\begingroup$ The technique GCC uses in place of the Continuation-Passing Style you show here is, I believe, Static Single Assignment Form. $\endgroup$
    – Davislor
    Commented Dec 30, 2017 at 8:46
  • $\begingroup$ @Davislor While related to CPS, SSA does not affect the control flow of a procedure nor does it reify the stack (or otherwise introduce data structures that would need to be dynamically allocated). As related to SSA, CPS "does too much" which is why Administrative Normal Form (ANF) is a closer fit to SSA. So GCC uses SSA, but SSA doesn't lead to the control stack being viewable as a manipulable data structure. $\endgroup$ Commented Dec 30, 2017 at 9:55
  • $\begingroup$ Right. I was responding to, “I'm not saying GCC does all this reasoning at compile-time. I don't know what logic GCC uses.” My answer, similarly, was showing that the transformation is theoretically justified, not saying that it’s the implementation method any given compiler uses. (Although, as you know, many compilers do transform a program into CPS during optimization.) $\endgroup$
    – Davislor
    Commented Dec 30, 2017 at 17:11
9
$\begingroup$

I’m going to beat around the bush for a while, but there is a point.

Semigroups

The answer is, the associative property of the binary reduction operation.

That’s pretty abstract, but multiplication is a good example. If x, y and z are some natural numbers (or integers, or rational numbers, or real numbers, or complex numbers, or N×N matrices, or any of a whole bunch more things), then x×y is the same kind of number as both x and y. We started with two numbers, so it’s a binary operation, and got one, so we reduced the count of numbers we had by one, making this a reduction operation. And (x×yz is always the same as x×(y×z), which is the associative property.

(If you already know all this, you can skip to the next section.)

A few more things you often see in computer science that work the same way:

  • adding any of those kinds of numbers instead of multiplying
  • concatenating strings ("a"+"b"+"c" is "abc" whether you start with "ab"+"c" or "a"+"bc")
  • Splicing two lists together. [a]++[b]++[c] is similarly [a,b,c] either from back to front or front to back.
  • cons on a head and tail, if you think of the head as a singleton list. That’s just concatenating two lists.
  • taking the union or the intersection of sets
  • Boolean and, Boolean or
  • bitwise &, | and ^
  • composition of functions: (fg)∘h x = f∘(gh) x = f(g(h(x)))
  • maximum and minimum
  • addition modulo p

Some things that don’t:

  • subtraction, because 1-(1-2) ≠ (1-1)-2
  • xy = tan(x+y), because tan(π/4 + π/4) is undefined
  • multiplication over the negative numbers, because -1 × -1 is not a negative number
  • division of integers, which has all three problems!
  • logical not, because it has only one operand, not two
  • int print2(int x, int y) { return printf( "%d %d\n", x, y ); }, as print2( print2(x,y), z ); and print2( x, print2(y,z) ); have different output.

It’s a useful enough concept that we named it. A set with an operation that has these properties is a semigroup. So, the real numbers under multiplication are a semigroup. And your question turns out to be one of the ways this kind of abstraction becomes useful in the real world. Semigroup operations can all be optimized the way you’re asking about.

Try This At Home

As far as I know, this technique was first described in 1974, in Daniel Friedman and David Wise’s paper, “Folding Stylized Recursions into Iterations”, although they assumed a few more properties than it turns out they needed to.

Haskell is a great language to illustrate this in, because it has the Semigroup typeclass in its standard library. It calls the operation of a generic Semigroup the operator <>. Since lists and strings are instances of Semigroup, their instances define <> as the concatenation operator ++, for example. And with the right import, [a] <> [b] is an alias for [a] ++ [b], which is [a,b].

But, what about numbers? We just saw that numeric types are semigroups under either addition or multiplication! So which one gets to be <> for a Double? Well, either one! Haskell defines the types Product Double, where (<>) = (*) (that is the actual definition in Haskell), and also Sum Double, where (<>) = (+).

One wrinkle is that you used the fact that 1 is the multiplicative identity. A semigroup with an identity is called a monoid and is defined in the Haskell package Data.Monoid, which calls the generic identity element of a typeclass mempty. Sum, Product and list each has an identity element (0, 1 and [], respectively), so they are instances of Monoid as well as Semigroup. (Not to be confused with a monad, so just forget I even brought those up.)

That’s enough information to translate your algorithm into a Haskell function using monoids:

module StylizedRec (pow) where

import Data.Monoid as DM

pow :: Monoid a => a -> Word -> a
{- Applies the monoidal operation of the type of x, whatever that is, by
 - itself n times.  This is already in Haskell as Data.Monoid.mtimes, but
 - let’s write it out as an example.
 -}
pow _ 0 = mempty -- Special case: Return the nullary product.
pow x 1 = x      -- The base case.
pow x n = x <> (pow x (n-1)) -- The recursive case.

Importantly, note that this is tail recursion modulo semigroup: every case is either a value, a tail-recursive call, or the semigroup product of both. Also, this example happened to use mempty for one of the cases, but if we hadn’t needed that, we could have done it with the more general typeclass Semigroup.

Let’s load this program up in GHCI and see how it works:

*StylizedRec> getProduct $ pow 2 4
16
*StylizedRec> getProduct $ pow 7 2
49

Remember how we declared pow for a generic Monoid, whose type we called a? We gave GHCI enough information to deduce that the type a here is Product Integer, which is an instance of Monoid whose <> operation is integer multiplication. So pow 2 4 expands recursively to 2<>2<>2<>2, which is 2*2*2*2 or 16. So far, so good.

But our function uses only generic monoid operations. Previously, I said that there is another instance of Monoid called Sum, whose <> operation is +. Can we try that?

*StylizedRec> getSum $ pow 2 4
8
*StylizedRec> getSum $ pow 7 2
14

The same expansion now gives us 2+2+2+2 instead of 2*2*2*2. Multiplication is to addition as exponentiation is to multiplication!

But I gave one other example of a Haskell monoid: lists, whose operation is concatenation.

*StylizedRec> pow [2] 4
[2,2,2,2]
*StylizedRec> pow [7] 2
[7,7]

Writing [2] tells the compiler that this is a list, <> on lists is ++, so [2]++[2]++[2]++[2] is [2,2,2,2].

Finally, an Algorithm (Two, in Fact)

By simply replacing x with [x], you convert the generic algorithm that uses recursion modulo a semigroup into one that creates a list. Which list? The list of elements the algorithm applies <> to. Because we used only semigroup operations that lists have too, the resulting list will be isomorphic to the original computation. And since the original operation was associative, we can equally well evaluate the elements from back to front or from front to back.

If your algorithm ever reaches a base case and terminates, the list will be non-empty. Since the terminal case returned something, that will be the final element of the list, so it will have at least one element.

How do you apply a binary reduction operation to every element of a list in order? That’s right, a fold. So you can substitute [x] for x, get a list of elements to reduce by <>, and then either right-fold or left-fold the list:

*StylizedRec> getProduct $ foldr1 (<>) $ pow [Product 2] 4
16
*StylizedRec> import Data.List
*StylizedRec Data.List> getProduct $ foldl1' (<>) $ pow [Product 2] 4
16

The version with foldr1 actually exists in the standard library, as sconcat for Semigroup and mconcat for Monoid. It does a lazy right fold on the list. That is, it expands [Product 2,Product 2,Product 2,Product 2] to 2<>(2<>(2<>(2))).

This is not efficient in this case because you can’t do anything with the individual terms until you generate all of them. (At one point I had a discussion here about when to use right folds and when to use strict left folds, but it went too far afield.)

The version with foldl1' is a strictly-evaluated left fold. That is to say, a tail-recursive function with a strict accumulator. This evaluates to (((2)<>2)<>2)<>2, calculated immediately and not later when it’s needed. (At least, there are no delays within the fold itself: the list being folded is generated here by another function that might contain lazy evaluation.) So, the fold calculates (4<>2)<>2, then immediately calculates8<>2, then 16. This is why we needed the operation to be associative: we just changed the grouping of the parentheses!

The strict left fold is the equivalent of what GCC is doing. The leftmost number in the previous example is the accumulator, in this case a running product. At each step, it’s multiplied by the next number in the list. Another way to express it that is: you iterate over the values to be multiplied, keeping the running product in an accumulator, and on each iteration, you multiply the accumulator by the next value. That is, it’s a while loop in disguise.

It can sometimes be made just as efficient. The compiler might be able to optimize away the list data structure in memory. In theory, it has enough information at compile time to figure out it should do so here: [x] is a singleton, so [x]<>xs is the same as cons x xs. Each iteration of the function might be able to re-use the same stack frame and update the parameters in place.

Either a right fold or a strict left fold could be more appropriate, in a particular case, so know which one you want. There are also some things only a right fold can do (such as generate interactive output without waiting for all the input, and operate on an infinite list). Here, though, we’re reducing a sequence of operations to a simple value, so a strict left fold is what we want.

So, as you can see, it is possible to automatically optimize tail-recursion modulo any semigroup (one example of which is any of the usual numeric types under multiplication) to either a lazy right fold or a strict left fold, in one line of Haskell.

Generalizing Further

The two arguments of the binary operation don’t have to be the same type, so long as the initial value is the same type as your result. (You can of course always flip the arguments to match the order of the kind of fold you’re doing, left or right.) So you might repeatedly add patches to a file to get an updated file, or starting with an initial value of 1.0, divide by integers to accumulate a floating-point result. Or prepend elements to the empty list to get a list.

Another type of generalization is to apply the folds not to lists but to other Foldable data structures. Often, an immutable linear linked list is not the data structure you want for a given algorithm. One issue I did not get into above is that it’s a lot more efficient to add elements to the front of a list than to the back, and when the operation is not commutative, applying x on the left and the right of the operation aren’t the same. So you would need to use another structure, such as a pair of lists or binary tree, to represent an algorithm that could apply x on the right of <> as well as to the left.

Also note that the associative property allows you to regroup the operations in other useful ways, such as divide-and-conquer:

times :: Monoid a => a -> Word -> a
times _ 0 = mempty
times x 1 = x
times x n | even n    = y <> y
          | otherwise = x <> y <> y
  where y = times x (n `quot` 2)

Or automatic parallelism, where each thread reduces a subrange to a value that is then combined with the others.

$\endgroup$
1
  • 2
    $\begingroup$ We can do an experiment to test that associativity is the key to GCC's ability to do this optimization: a pow(float x, unsigned n) version gcc.godbolt.org/z/eqwine only optimizes with -ffast-math, (which implies -fassociative-math. Strict floating point is of course not associative because different temporaries = different rounding). The introduces a 1.0f * x that wasn't present in the C abstract machine (but which will always give an identical result). Then n-1 multiplications like do{res*=x;}while(--n!=1) are the same as the recursive, so this is a missed optimization. $\endgroup$ Commented Mar 16, 2019 at 10:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.