The Music of Streams on a strict instrument




The cleverest copyist is the one whose music is performed with the most ease without the performer guessing why. -- J.J.Rousseau, Dictionary of Music.
The appeal to music is not just a flourish or a veiled reference. It is the point: how pleasant is a strict language to play the classical themes from analysis: infinite power series.

As the instrument to to deal with infinite structures, lazy evaluation -- evaluating an expression only when needed and to the needed extent, memoizing the result -- has the undeniable appeal. Just look at the example, chosen to show off Haskell at

    primes = filterPrime [2..]
      where filterPrime (p:xs) = p : filterPrime [x | x <- xs, x `mod` p /= 0]
Haskell does not hold the monopoly on lazy evaluation: almost any strict language supports it, just not by default. For example, the straightforward implementation of primes in OCaml is as follows:
    type 'a ll = 'a lv Lazy.t and  'a lv = Cons of 'a * 'a ll | Nil
    let rec filter p : 'a ll -> 'a ll = function
      | lazy Nil -> lazy Nil
      | lazy (Cons (x,t)) ->
          let t = lazy (Lazy.force @@ filter p t) in 
          if p x then lazy (Cons (x,t)) else t
    let rec iota x : int ll = lazy (Cons (x,iota (x + 1)))
    let primes = 
      let rec filter_prime = function
        | lazy (Cons (p,xs)) -> 
            lazy (Cons (p,filter_prime @@ filter (fun x -> x mod p <> 0) xs))
      in filter_prime @@ iota 2
Haskell however looks quite more attractive in comparison. It is not just the clutter of lazy and force that spoils. It is the line lazy (Lazy.force @@ filter p t). The simple filter p t in its place would have type-checked as well. Guess how that code would have worked though. Thus not only the programmers have to use lazy and force; they have to know where to use them.

One may wonder if the burden of correctly placing lazy annotations is the inherent defect of strict languages, or one can do better. And what a better example to test this than the most elegant and perhaps the most sophisticated use of lazy evaluation: power streams. Such was the challenge posed by Kim-Ee Yeoh: write Doug McIlroy's power serious one-liners without thunks and lazy annotations. Doug McIlroy started his ``Music of streams'' paper with the J.J.Rousseau's quote. Uncannily, the quote talks about our challenge: play power streams without the performer guessing their implementation -- without even being aware if he is working in a strict or lazy language.

M.D. McIlroy: Power serious: power series in ten one-liners. July 2007

M.D. McIlroy: The music of streams
Information Processing Letters 77 (2001) 189-195
<> [9K]
The complete OCaml code for the article, with more examples and tests

The discussion thread on Haskell-Cafe, 25 December 2015 and January 2016


Power series: Haskell v. OCaml

Our subject is hence the infinite power series a0 + a1*x + a2*x^2 + ... + an*x^n + ..., or, in Horner rule, a0 + x*(a1 + x*(a2 + ...)). It is represented as the non-ending stream of its coefficients [a0, a1, ... an, ...]. The series is always infinite, although it may have a zero tail. McIlroy later describes the extension to represent polynomials finitely, which is out of scope for this article.

The table below contrasts McIlroy's famous one-liners for power stream manipulation, written in Haskell and OCaml. The Haskell code comes verbatim from McIlroy's web page, which also explains the underlying mathematics. In his convention, a series variable has suffix s, or t when it is a tail. The OCaml code is written with the library to be introduced later; its functions are qualified with the module name I. We show the inferred signatures in the comments. When comparing code one has to keep in mind that Haskell and OCaml differ in many respects, not just the evaluation strategy. OCaml is in general a little bit more verbose and ungainly at places: its syntax had to be compatible with the original (heavy) Caml, which in turn inherited the syntax from ML of late 70's, limited by the parsing technology of the day. Recursive definitions have to be marked explicitly. OCaml also does not have overloading: That is why we see not just + but also +. (for float addition) and +% (the addition of infinite series).

Coerce scalar to series
    series f = f : repeat 0
    let series f = I.cons f I.repeat 0.
    (* val series : float -> float I.t *)
Promote integer constants
    fromInteger c = series(fromInteger c)
    let int x = series @@ float_of_int x
    (* val int : int -> float I.t *)
    negate (f:ft) = -f : -ft
    let rec negate fs = I.decon fs @@ fun f ft -> I.cons (-. f) negate ft
... faster
    negate = map negate
    let negate = (fun x -> -. x)
    (f:ft) + (g:gt) = f+g : ft+gt
    let rec (+%) fs gs = I.decon2 fs gs @@ fun f ft g gt -> I.cons (f +. g) (fun (ft,gt) -> ft +% gt) (ft,gt)
... faster
    (+) = zipWith (+)
    let (+%) = I.zip_with (+.)
    (f:ft) * gs@(g:gt) = f*g : ft*gs + series(f)*gt
    let rec ( *% ) fs gs = I.decon2 fs gs @@ fun f ft g gt -> I.cons (f *. g) (fun gt -> ft *% gs +% series f *% gt) gt
    -- Tying-the-knot trick
    (f:ft) / (g:gt) = qs where qs = f/g : series(1/g)*(ft-qs*gt)
    let ( /% ) fs gs = I.decon2 fs gs @@ fun f ft g gt -> I.fix @@
              I.cons (f /. g) (fun qs -> series (1. /. g) *% (ft -% qs *% gt))
    (f:ft) # gs@(0:gt) = f : gt*(ft#gs)
    (* Since # is reserved in OCaml, we use %% *)
    let rec ( %% ) fs gs = I.decon2 fs gs @@ fun f ft 0. -> I.cons f (fun gt -> gt *% (ft %% gs))
Reversion (compositional inverse)
    -- Tying-the-knot trick
    revert (0:ft) = rs where rs = 0 : 1/(ft#rs)
    let revert fs = I.decon fs @@ fun 0. ft -> I.fix @@ I.cons 0. @@ fun rs -> int 1 /% (ft %% rs)
    (* val revert : float I.t -> float I.t *)
    -- integral from 0 to x
    int fs = 0 : zipWith (/) fs [1..]
    let integ = I.cons 0. (fun fs -> I.zip_with (/.) fs (iota 1.))
    (* val integ : float I.t -> float I.t *)
    -- type (Num a,Enum a)=>[a]->[a]
    diff (_:ft) = zipWith (*) ft [1..]
    let diff fs = I.decon fs @@ fun _ ft -> I.zip_with ( *. ) ft (iota 1.)
    (* val diff : float I.t -> float I.t *)
    tans = revert(int(1/(1:0:1)))
    let tans = revert @@ integ (int 1 /% from_list [1.;0.;1.])
    (* val tans : float I.t *)
Sine and cosine
    -- (Implicit) Mutual recursion
    sins = int coss
    coss = 1 - int sins
    let (sins,coss) = I.fix2 (fun (s,c) -> (integ c, int 1 -% integ s))
Final test
    test4 = tans - sins/coss
    let test4 = tans -% sins /% coss
The series for tan x is in terms of the series for its functional inverse, arctan x = INT dx/(1+x^2). The final test (Music of Streams paper, Example 4) exercises all facilities, and is demanding. It seems to run a bit faster in the OCaml byte-code interpreter than in GHCi. As promised, there are no thunks or lazy annotations.


DSL for infinite streams

The OCaml power-series one-liners rely on the library for infinite streams of the following interface:
    module type INFSER = sig
      type 'a t
      val decon  : 'a t -> ('a -> 'a t -> 'w t) -> 'w t  (* deconstructor *)
      (* Another deconstructor, although technically unnecessary, but convenient *)
      val decon2 : 'a t -> 'b t -> ('a -> 'a t -> 'b -> 'b t -> 'w t) -> 'w t
      val cons   : 'a -> ('b -> 'a t) -> 'b -> 'a t      (* constructor *)
      val repeat   : 'a -> 'a t                          (* same elems *)
      val map      : ('a -> 'b) -> ('a t -> 'b t)
      val zip_with : ('a -> 'b -> 'c) -> ('a t -> 'b t -> 'c t) 
      val take : int -> 'a t -> 'a list
      val fix  : ('a t -> 'a t) -> 'a t
      val fix2 : ('a t * 'b t -> 'a t * 'b t) -> 'a t * 'b t
Using the interface is straightforward, although its design is not. The case in point is cons, with an unexpected type. Another surprise is the deconstructor decon, whose result must be a power series. A function that deconstructs and consumes an infinite series must itself be producing an infinite series. INFSER lets us re-implement the earlier filter as:
    let rec filter : ('a -> bool) -> 'a I.t -> 'b I.t = fun p l ->
      I.decon l (fun x t -> 
        let t = filter p t in 
        if p x then I.cons x (fun x -> x) t else t)
Now we can write filter p t for the recursive invocation, without ill-effects.

Again, INFSER lets us write power serious without the clutter of lazy, force, thunks and other suspensions that had to be correctly placed by the user. It almost seems that the programmer does not have to be aware of the default evaluation strategy. Yet the recursion gives the strictness away. Although we can write recursive definitions just like in Haskell (see, e.g., the code for addition and multiplication), and although the explicit fixpoint for tying the knot may occur with either strategy (cf. mfix in Haskell) and should probably be encouraged given its subtlety, recursive series definitions seem to be possible only in a lazy language. As Doug McIlroy pointed out in the discussion thread, the exponential series can be defined in Haskell simply as

    exps = 1 + integral exps
The similar definition in OCaml let rec exps = int 1 +% integ exps is not accepted: ``This kind of expression is not allowed as right-hand side of let rec''. One is forced to use the explicit fixpoint:
    let exps = I.fix @@ fun exps -> int 1 +% integ exps
Again one wonders how deep is this problem and if something can be done, in a hypothetical strict language -- or even in OCaml. Recall, a recursive definition let rec f = e may be treated as a syntax sugar for let f = fix (fun f -> e). Such a `macro-expansion' is valid in any language, strict or lazy. If only we could sneak-in our own fix -- based on the inferred type of the expression or in some other way... OCaml attributes do give us `the other way'. With the appropriate PPX extension, one can write the desired
    let[@stream I] rec exps = int 1 +% integ exps
and having it expand to the earlier I.fix code. Again, such an expansion is valid regardless of the evaluation strategy. Thus elegant recursive power-series definitions can in principle be written in a strict language, or even in OCaml with some extensions.

Speaking of elegance, I can't help but quote the paragraph from Doug McIlroy's web page:

``Extensions to handle polynomials make a practical package, doubled in size, not as pretty, but much faster and capable of feats like pascal. To see the dramatic speedup, try a bigger test like take 20 tans.''
Indeed, ``not as pretty, but much faster''.



In the upshot, well-chosen combinators are far more important than strict/lazy evaluation differences. A language expressive enough to support mini-languages (embedded DSLs) lets us deal with infinite series without much trouble or clutter. Strict evaluation is not an obstacle: if we can define a DSL, we can make it follow any evaluation strategy we want. Tom Ellis summarized it well: ``we can have access to laziness within strict languages with no less elegance than we have access to IO in pure languages.''

I thank Tom Ellis for the productive discussion and Kim-Ee Yeoh for posing such an interesting and well-defined problem.