It is a surprisingly concrete and useful application of abstract nonsense.
One approach of transforming mutually primitive recursive equations to first-order tail-recursive computations (i.e., for-loops by another name), forcefully advocated by Danvy, is the continuation-passing style (CPS) transformation followed by defunctionalization. A creative intermediate rearrangement step might be needed, however. As was pointed by Gibbons, that intermediate rearrangement often depends on associativity, which is not generally given. There is also a dissatisfaction from first introducing lots of functions and closures (in the CPS transformation step), and trying to eliminate them later.
On the other hand, Oliveira has described a direct and general method of turning mutually primitive recursive equations over natural numbers to a for-loop, with no assumptions about the operations -- not even associativity. However, the explanation is embedded in an extensive category-theoretical analysis of the for-loop combinator, and is hard to understand and easy to dismiss. Phrases like `mutual catamorphisms' and `exponential adjunction' tend to turn off practical programmers.
This page presents the Oliveira's method in simpler terms, and on simple but interesting examples. The method becomes easy, and even somewhat intuitive, to apply -- without background in abstract algebra. (The background, however, helps one appreciate that the method is quite more general than presented.)
The second event was reading Oliveira's note on the under-appreciated for-loop, the day before Gibbons' talk. The bulk of the note is the category-theoretical analysis of the for-loop combinator and showing it is left-dual to itself -- or, in general, that Peano recursion is adjoint to tail-recursion, via contravariant exponential adjunction. Yet, Sec 2 has caught my eye: it was about calculating for-loops from recursive equations. Stood out was the phrase ``In a sense, the mutual recursion law (7) gives us a hint on how many global variables are needed and how they `are born' in computer programs, out of the maths definitions themselves.'' Not often one sees recursion laws and global mutable variables mentioned in the same sentence.
Listening to Gibbons' talk and the subtractorial challenge reminded me of Oliveira's note. Would Oliveira's method apply to subtractorial? Would the lack of associativity be a stumbling block? After the talk I jotted down the calculation. The method applied, rather easily; associativity was not needed. The resulting code seemed to work. Yet I did not believe it. It just seemed too magical.
The next morning I tried to derive the efficient subtractorial code in an obvious, and obviously believable style. The key word is deriving: obtaining the for-loop code by obviously correct steps from the intuitively adequate specification. It turns out easily possible. In fact, any factorial-like function -- applying a binary operation to a progression of integers -- can be converted to an accumulator-passing style -- without assuming anything about the operation. The method is simple and does not require CPS transformation or defunctionalization, which I always felt to be a distraction.
I then compared thus derived for-loop code to what I calculated by Oliveira's method the evening before. To my astonishment, the results were identical.
<Programming>'23, Tokyo, Japan, March 15, 2023. Vol 6.
José N. Oliveira.
A Note on the Under-Appreciated For-Loop
Technical Report TR-HASLab:01:2020, HASLab/U.Minho and INESC TEC, 2020
<http://www.di.uminho.pt/~jno/ps/haslabtr202010.pdf>
J.N. Oliveira. Program Design by Calculation, 2019
Draft of textbook in preparation, current version: October 2019
Informatics Department, University of Minho
<http://www.di.uminho.pt/~jno/ps/pdbc.pdf>
recursion from iteration The inverse transformation
The method is simple and does not require the CPS transformation or defunctionalization. It is, however, particular to the factorial-like functions. See the following sections for a more general method.
Subtractorial is like the factorial with subtraction instead of multiplication (and the base case being 0 rather than 1):
let rec subt n = if n = 0 then 0 else n - subt (n-1)
subt 5;;
- : int = 3
subt 9;;
- : int = 5
To make sure that our solution does not depend on a particular property of
subtraction (not that there are that many), we generalize to an
arbitrary operation op and arbitrary initial value z:
let rec gfac (op:int->'a->'a) (z:'a) n : 'a =
if n = 0 then z else op n (gfac op z (n-1))
so that gfac ( * ) 1 is the familial factorial and gfac (-) 0 is
the subtractorial. In general, the second argument of op (and its
result) do not have to be integers.
Gratifyingly, the generalization helps us to see that
gfac is right fold. We use the OCaml's formulation
val List.fold_right : ('a -> 'b -> 'b) -> 'a list -> 'b -> 'b
with the properties
List.fold_right op [] z ≡ z
List.fold_right op (h::t) z ≡ op h (List.fold_right op t z)
In fact, the two properties uniquely define the right fold: the so called
universality property. Our gfac has similar properties:
gfac op z 0 ≡ z
gfac op z (n+1) ≡ op (n+1) (gfac op z n)
This similarity demonstrates -- in fact, proves, by appeal to
universality -- that gfac is representable as the right fold:
let gfac op z n = List.fold_right op (List.init n (fun i -> n - i)) zThe
List.init expression in parentheses returns the list [n;n-1;...1]:
List.init 5 (fun i -> 5 - i);;
- : int list = [5; 4; 3; 2; 1]
It seems like a step in a wrong direction: right fold is not
tail-recursive, and so does not obviously correspond to a loop.
Furthermore, we introduced the list [n;n-1;...1] as an intermediate
data structure, which takes time and memory to construct.
However, the solution is all but one step away. Right fold (on a finite list) can always be converted to the left fold, on the reverse list. Left fold, in the OCaml formulation
List.fold_left : ('a -> 'b -> 'a) -> 'a -> 'b list -> 'a
is defined by the recursive equations
List.fold_left op z [] ≡ z
List.fold_left op z (h::t) ≡ List.fold_left op t (op z h)
and is clearly tail-recursive.
The connection between the right and left folds is clear
from the picture (for a sample list of three elements):
List.fold_right ⊕ [a;b;c] z ≡ (a ⊕ (b ⊕ (c ⊕ z))) ≡ ((z ⊖ c) ⊖ b) ⊖ a
List.fold_left ⊖ z [a;b;c] ≡ (((z ⊖ a) ⊖ b) ⊖ c)
where ⊕ is an arbitrary binary operation, written in infix,
and ⊖ is the same operation with the arguments flipped.
We stress that nothing is assumed about the properties of ⊕.
We thus obtain
let gfac op z n = List.fold_left (Fun.flip op) z (List.init n (fun i -> i+1))where
List.init n (fun i -> i+1) is the list [1;...n].
That's it. Left fold, especially over the list [1;...n] is
the accumulating for-loop in disguise: the loop index is the list
element. Thus
let gfac op z n = for_loop 1 n op zwhere the for-loop, in the functional style, is
let for_loop (from:int) (to_:int) (op:int->'a->'a) (z:'a) : 'a =
let rec loop acc i = if i > to_ then acc else loop (op i acc) (i+1)
in loop z from
The for-loop is more familiar in the imperative style:
let for_loop (from:int) (to_:int) (op:int->'a->'a) (z:'a) : 'a =
let acc = ref z in
for i = from to to_ do acc := op i !acc done; !acc
Coming back to subtractorial specifically, it can be computed as
let subt n =
let acc = ref 0 in
for i = 1 to n do acc := i - !acc done; !acc
subt 5;;
- : int = 3
subt 9;;
- : int = 5
The Oliveira's method is expressed in one, rather cryptic equation, called the mutual recursion law:
⎧ f ⋅ in = h ⋅ F <f,g>
<f,g> = ⦇ <h,k> ⦈ ⇔ ⎨
⎩ g ⋅ in = k ⋅ F <f,g>
Here is what it amounts to, in a more familiar notation and applied to
natural numbers. (The mutual recursion law is more general: it applies
to any inductive data type.)
Suppose functions f : int -> a and g : int -> b
(for some types a and b) are defined by the
following system of recursive equations:
f 0 = zf
f (n+1) = s (f n) (g n)
g 0 = zg
g (n+1) = t (f n) (g n)
where zf : a and zg : b are some values, and
s : a -> b -> a and t : a -> b -> b are some functions.
Then f n and g n can be computed with the for-loop:
(f n, g n) = for_loop 0 (n-1) (fun _ (x,y) -> (s x y, t x y)) (zf,zg)where, to remind,
let for_loop (from:int) (to_:int) (op:int->'a->'a) (z:'a) : 'a =
let acc = ref z in
for i = from to to_ do acc := op i !acc done; !acc
In other words, f n and g n are computed by the for-loop with the
state zf and zg updated by the pair of functions s and t.
The next two sections apply this method on concrete examples. The rest of this section, which may be skipped, explains the notation of the mutual recursion law and shows how the recursive equations came out of it.
First, the notation.
The product of two functions f:a->c and g:b->d (where a,b,c
and d are some types) is notated f*g : a*b -> c*d. It is the
function defined as
(f*g) (x,y) = (f x, g y)where
(x,y) is the notation for pairs (as in OCaml).
We may also tuple functions: if h1: a->c and h2: a->d,
(which have to have the same argument type), the tupling
<h1,h2> : a -> c*d is the function defined as
<h1,h2> x = (h1 x, h2 x)We write
a+b for a disjoint union of the types a and b (the
OCaml notation is the (a,b) either data type, with the constructors
Left and Right). The same notation is extended to
functions (such as f and g above):
f+g : a+b -> c+d is the function
(f+g) (Left x) = Left (f x)
(f+g) (Right y) = Right (g y)
If h1' : a -> c and h2' : b -> c (the result types must be the
same), their cotupling [h1',h2'] : a+b -> c is defined as
[h1',h2'] (Left x) = h1' x
[h1',h2'] (Right y) = h2' y
The co-tupling hence corresponds to the case analysis.
If we are interested in for-loops over natural numbers, then the remaining notation in the mutual recursion law is
in = [const 0, succ]
F X = 1 + X
It is easier to understand what F X does if we substitute X with
some function f. Then, purely formally,
F f is 1 + f, where 1 is a constant function. In other words,
F f is the function
(1+f) (Left x) = c
(1+f) (Right y) = f y
where c is the constant associated with 1.
We are now ready to tackle the mutual recursion law. Let's start with
f ⋅ in = h ⋅ F <f,g>Suppose
f : int -> a. The left-hand--side means:
(f ⋅ in) = f ⋅ [const 0 n, succ n] = [f ⋅ const 0, f ⋅ succ] = [const (f 0), f ⋅ succ]since the functional composition distributes into tupling, as easy to see. For the right-hand side, we have
h ⋅ F <f,g> = h ⋅ (1 + <f,g>)Now,
(1 + <f,g>) produces a sum, so h must be a function that
accepts sum, that is, do a case analysis.
Any such function can be represented as a cotupling
h = [s',s] of some functions s' and s.
Continuing for the right-hand side:
[s',s] ⋅ F <f,g> = [s',s] ⋅ (1 + <f,g>) = [s'⋅1, s⋅<f,g>] = [const (s' c), s⋅<f,g>]The left- and the right-hand--sides must be equal, from which follows:
f 0 = s' c
f ⋅ succ = s⋅<f,g>
Let's call the result s' c as zf. For the second line,
we apply n to both sides. The result is:
f 0 = zf
f (n+1) = s (f n, g n)
Which is the recursive equation for f we saw earlier.
The equation g ⋅ in = k ⋅ F <f,g> is understood analogously.
Once we figured out what h = [s',s] (and, analogously,
k = [t',t]) are, the part
<f,g> = ⦇ <h,k> ⦈of the mutual recursion law says that
f and g can be computed as
as a fold whose state is a pair: the values of f n and g n at the
current step n. The initial state is determined by
s' and t' and updated by s and t. The fold on natural numbers
is the for-loop (as Oliveira's note explains at the beginning).
gfact, which is, recall, defined
as
let rec gfac (op:int->'a->'a) (z:'a) n : 'a =
if n = 0 then z else op n (gfac op z (n-1))
Notating gfac op z n as f n, the above OCaml code can be re-written
as the following equations:
f 0 = z
f n = op n (f (n-1))
which can be adjusted to match the form of the mutual-recursion law:
f 0 = z
f (n+1) = op (n+1) (f n)
The right-hand should contain only functions
applied to n, but we have n+1. Therefore, we have to introduce another
function g n = n+1. It may also be defined by the
similar recursive equations:
f 0 = z
f (n+1) = op (g n) (f n)
g 0 = 1
g (n+1) = 1 + g n
The equations match the mutual-recursion law; the function
s is fun x y -> op y x and t is the successor. We immediately obtain
let gfac op z n =
let (fn, gn) = for_loop 0 (n-1) (fun _ (x,y) -> (op y x, y+1)) (z,1)
in fn
or, to expand
let gfac op z n =
let fn = ref z in
let gn = ref 1 in
for i = 0 to (n-1) do
fn := op !gn !fn;
gn := !gn + 1
done; !fn
One notices that the current value of gn is i+1 where i is the
loop index. The code can then be simplified to
let gfac op z n =
let fn = ref z in
for i = 1 to n do
fn := op i !fn
done; !fn
which is what we have derived in Any factorial-like function, as a for-loop.
cos x = SUM_{i=0} (-1)^i / (2i)! x^(2i)
Suppose x is fixed; let's f n be an approximation of cos x:
the partial sum from i=0 through i=n. It is immediate that
f 0 = 1
f (n+1) = (f n) + (-1)^(n+1) / (2n+2)! x^(2n+2)
The second factor on the second line looks like a complicated
function of n; call it g. We see
g 0 = (-1)^(1) / 2! x^2 = -x^2/2
g (n+1) = (-1)^(n+2) / (2n+4)! x^(2n+4)
= - (-1)^(n+1)/((2n+2)! (2n+3) (2n+4)) x^(2n+2) x^2
= (g n) (- x^2/((2n+3)(2n+4)))
= (g n) (- x^2/(h n))
where we have introduced h n = (2n+3)(2n+4). It is a rather simple
function and does not have to be computed inductively. For the sake
of example, we nevertheless do it:
h 0 = 12
h (n+1) = (2n+3+2)(2n+4+2) = (h n) + 8n + 18
Finally, we introduce k n = 8n + 18 and calculate
k 0 = 18
k (n+1) = k n + 8
Collecting all the equations and applying the mutual-recursion law gives:
let cosn x n =
let x2 = -. x *. x in
let (fn, gn, hn, kn) =
for_loop 0 (n-1)
(fun _ (f,g,h,k) -> (f+.g, g*.x2/. float_of_int h, h+k, k+8))
(1., x2 /. 2., 12, 18)
in fn
or, imperatively
let cosn x n =
let x2 = -. x *. x in
let fn = ref 1. and gn = ref (x2 /. 2.) and hn = ref 12 and kn = ref 18 in
for i = 0 to n-1 do
fn := !fn +. !gn;
gn := !gn *. x2 /. float_of_int !hn;
hn := !hn + !kn;
kn := !kn + 8
done; !fn
cosn (Float.pi /. 3.) 5;;
- : float = 0.499999996390943224
Each iteration includes the constant amount of work:
floating point addition, multiplication, division, and a couple of
integer additions. If one has not seen the derivation,
the constants like 12 and 18 in the code appear strange, coming
out of the blue.
Since the loop index is left unused, we can convert the loop to a while loop
and iterate until the difference from the previous cos x estimate
becomes small enough:
let cosn eps x =
let x2 = -. x *. x in
let fn = ref 1. and gn = ref (x2 /. 2.) and hn = ref 12 and kn = ref 18 in
while Float.abs !gn > eps do
fn := !fn +. !gn;
gn := !gn *. x2 /. float_of_int !hn;
hn := !hn + !kn;
kn := !kn + 8
done; !fn
cosn 1e-10 (Float.pi /. 3.);;
- : float = 0.500000000021777691