Squinting at Power Series …

                                 Squinting at Power Series

                                      M. Douglas McIlroy
                                   AT&T Bell Laboratories
                                 Murray Hill, New Jersey 07974


                                         ABSTRACT

       Data streams are an ideal vehicle for handling power series. Stream
       implementations can be read off directly from simple recursive equations
       that define operations such as multiplication, substitution, exponentiation,
       and reversion of series. The bookkeeping that bedevils these algorithms
       when they are expressed in traditional languages is completely hidden
       when they are expressed in stream terms. Communicating processes are
       the key to the simplicity of the algorithms. Working versions are
       presented in newsqueak, the language of Pike's ``squint'' system; their
       effectiveness depends critically on the stream protocol.


Series and streams
      Power series are natural objects for stream processing. Pertinent computations are
neatly describable by recursive equations. CSP (communicating sequential process)
languages 1, 2 are good vehicles for implementation. This paper develops the algorithmic
ideas and reduces them to practice in a working CSP language, 3 not coincidentally
illustrating the utility of concurrent programming notions in untangling the logic of some
intricate computations on sequences.
       This elegant approach to power series computation, first demonstrated but never
published by Kahn and MacQueen as an application of their data stream system, is still
little known. 4 Power series are represented as streams of coefficients, the exponents
being given implicitly by ordinal position. (This is an infinite analogue of the familiar
representation of a polynomial as an array of coefficients.) Thus the power series for the
exponential function
                                           1      1       1
            ex =    0 x n / n!
                   n=
                                 = 1 + x + _ x2 + _ x3 + _ x4 + . . .
                                           2
                                            _
                                                  6
                                                   _      __
                                                         24
                                                 -2-


is represented as a stream of rationals:

                    1 1 __
                     _ _ 1
              1, 1, _ , _ , _ , . . .
                    2 6 24

      Throughout the paper power series will be denoted by symbols such as F, with the
functional form F(x) being used when the variable must be made explicit. Subscripts
denote individual coefficients:

                          
              F(x) =     0 F i x i .
                        i=


The corresponding stream will be referred to by the same symbol in program font: F. In
programs the names of power series variables contain capital letters; scalars are all lower
case. When new power series are calculated from old, the input series will always be
called F and G.
      In this section all programs are written in pseudocode; later they will be translated
into the language ``newsqueak'' to be run in the ``squint'' system.* The methods are
strictly formal; analytic interpretation of results depends on the convergence of the power
series in question. Thus assertions such as  i x i = 1/( 1 - x) are understood to hold
only where the series converge.
Addition. The sum, S = F + G, is easy. Since S i = F i + G i , all we need is a looping
process, which gets pairs of coefficients from the two inputs and puts their sum on the
output S. The following pseudocode suffices.

      # calculate power series S = F + G
      loop forever
              let f = get(F),
                   g = get(G)
              put(f+g, S)

Multiplication by a term. Simple termwise operations on single power series, such as
multiplication by a constant, integration, or differentiation are equally easy to program
and are left to the reader's imagination. Multiplication of a power series by x to yield
P = x F involves a one-element delay:
* The name of the language reveals descent from a mousier ancestor intended for programming a terminal
with multiple asynchronous windows. The name ``squint'' suggests ``squeak interpreter.''
                                             -3-


     # calculate power series P = x*F
     put(0, P)
     loop forever
             put(get(F), P)

Multiplication. Calculation of the product, P = F G, is more challenging. By equating
coefficients in

                                              
                 P n x n =   F i x i   G j x j ,
             n =0          i = 0     j = 0    

we obtain the familiar convolution for terms of the product,
                     n
             Pn =     Fi Gn -i ,
                    i =0
                                                                                         (1)


a most unpromising formula to program directly. To calculate the nth coefficient we
must have stored up the first n terms of both inputs, thus sacrificing most benefits of the
stream representation. It will be best to look in another direction, guided by the adage:
When iteration gets sticky, try recursion.
     We treat streams recursively, as we do lists, by handling the leading term and
recurring on the tail. This means rewriting a power series as a first term plus x times
another power series:
                          _
              F = F 0 + x F.
         _
The tail F is a whole power series beginning with a constant term; it is ripe for recursion.
     Writing the multiplication P = F G recursively, we obtain
                    _                     __         _         _ _
                                                                 _
            P 0 + x P = F 0 G 0 + x (F 0 G + G 0 F ) + x 2 F G .

Equate coefficients to find the first term of the product,

             P0 = F0 G0 ,

leaving for the tail,
             _        _
                      _      _     _ _
                                     _
             P = F 0 G + G 0 F + x F G,                                                  (2)

a sum of a constant F 0 times a power series, another constant G 0 times another power
                                                                    _ _
                                                                      _
series, and x times a third power series, which is itself a product F G. We already know
                                            -4-


how to multiply power series by a constant and how to add power series. We also know
how to multiply power series by x; recursion does the rest.
     In the following pseudocode, f*G denotes an auxiliary stream for the product f G,
and so on. Each auxiliary stream will be computed by another process. Thus f G will be
computed by a process that multiplies a power series by a constant. When two auxiliary
streams share an input stream, both see every element of the input.

     # calculate power series P = F*G
     let f = get(F),
                                                                           _
                                                                           _
         g = get(G)                       # F and G now "contain" Fbar and G
     put(f*g, P)                          # F0 G0
     let fG = f*G,
         gF = g*F,
         xFG = x*(F*G)       # here is the recursion
     loop forever
             put(get(fG) + get(gF) + get(xFG), P)

The convolution (1) has been hidden completely; all the necessary bookkeeping and
storage management is embodied in the topology and contents of the auxiliary streams.
Streams fan out to feed multiple consumers. For instance G enters into both f*G and
F*G. The topology ramifies as the recursion progresses; more auxiliary streams appear
at every stage, making the recursive picture of Figure 1. To write the program, though,
we need not think about the topology, for the program is a direct image of the
mathematics (2), and the mathematics is simple.
Substitution. Substitution, or composition, of power series, defined by S = F(G), may
be done similarly, using multiplication as a subprocess. The recursive formulation is
                     _                      _ _
                                             _             _
                                                           _
             S 0 + x S = F 0 + (G 0 + x G ) F (G 0 + x G ).

Equate coefficients to find the constant term,
                                              _
            S 0 = F 0 + G 0 . ( first term of F (G) ).

Induction shows that in general S 0 is an infinite sum.
                    
            S0 =     F i G0 .
                            i
                   i =0


To keep the problem finite, we further stipulate that G 0 = 0. Now
                                          -5-



                  F                                                 G




                 ×G 0                                             ×F 0


                                               _ __
                                               F ×G

                                          ×x




                                           +




                                               F×G



    Figure 1. Data flow for multiplication. The outer box represents the process for the
    tail of F G, equation (2). The similar inner box does not receive the first terms of F
    or G, and its output does not enter into the first two terms of the product F G. The
    newsqueak implementation uses a separate process for each box, and one or more
    processes for each stream splitting.

    Viewed as a three-dimensional scene receding into the distance, the picture shows a
    simple repeating structure connected on top by F and G buses. The buses grow,
    extending one level further with the receipt of each input term. A daisy chain of
    output connections leads from back to front.


                       _           _ _
                                   _
           S = S 0 + x S = F 0 + x G F (G) ,                                           (3)

from which pseudocode can be read off.
                                                     -6-


     # calculate power series S = F(G)
     put(get(F), S)          # F0
                                                                     _
     let FG = F(G)                                  # the recursion, F (G)
     let g = get(G)                                 # discard the first term (must be 0)
                                                      _ _
                                                      _
     let GF = G*FG                                  # G F (G)
     loop forever
             put(get(GF), S)

A subtlety: the G in let FG = F(G) stands for the whole power series, while the G in
                                          __
let GF = G*FG stands only for the tail G. In the same way, the meaning of F
changes between get(F) and let FG = F(G).
Exponentiation. We wish to compute X = e F . From the recursive representation,
                                _
                      F0
              X = e        e   xF
                                    ,                                                    (4)

we see that

                      F0
              X = e        . ( 1 + terms in x).

To make the constant term of X rational, we further require that F 0 = 0. Hence X 0 = 1.
     Unfortunately equation (4) does not lead to an algorithm. Formal substitution yields
                             _
              X(x) = 1 . X(x X ).                                                        (5)
                                                                   _
Since the first term of X(x) depends on the first term of X(x X ), we are stuck in an
infinite regress. By induction we may conclude that the first term is an infinite product of
ones, but what are the higher terms? We can ``cheat,'' and substitute F in the known
power series for the exponential. It is more satisfying, though, to build from nothing and
let the program ``discover'' the exponential for itself.
     A neat trick suffices. The identity,

                   d
              X = _ eF = eF F = X F ,
                  __
                  dx

may be read as a differential equation for X, with the formal solution,

                           x
              X(x) =        X(t)        F  (t) dt + c,
                           0
                                            -7-


where the constant of integration c is the constant term X 0 = 1. Now we can sketch the
program.

     # calculate X = exp(F), where F 0 = 0
     let D = deriv(F),
         P = X*D,
         I = integ(P, 1)
     loop forever
             put(get(I), X)

Here the stream operator deriv differentiates and the operator integ integrates with a
given constant of integration. The calculation involves a data stream loop: X enters as an
input into the calculation of X. Unlike the recursion in equation (5), this self-reference is
benign. Deadlock is avoided because the integral operator produces the first term of its
output (the constant of integration) before it needs any input. This result gets fed onto
stream X to help calculate the next term, and so on.
      The method of calculating power series by integration feeding on itself has a
distinguished lineage. It may be recognized as an instance of the classical Picard method
of successive approximations for solving y  = f (x,y). 5 It was first used in stream
context by Kahn and MacQueen. Abelson and Sussman give the degenerate case of
exp (x). 6
     Evidently the network of streams for calculating a product or an exponential entails
careful scheduling. Certain lazy-evaluation systems, such as Miranda, handle such
scheduling invisibly. 7 For example, in terms of lazily-evaluated infinite lists, we may
write the following functional pseudocode for power-series addition and multiplication.

     let add(F, G) = cons(head(F)*head(G), add(tail(F), tail(G)))

     let cmul(c, F) = cons(c*head(F), cmul(c, tail(F)))

     let mul(F,G) =
         let f = head(F),
             g = head(G),
             Fbar = tail(F),
             Gbar = tail(G)
         in cons(f*g, add(cmul(f, Gbar),
                          add(cmul(g, Fbar),
                              cons(0, mul(Fbar, Gbar)))))
                                           -8-


These functions follow from the same recursive equations as did the previous
pseudocode, with head and cons playing the role of get and put.

     Lazily-evaluated languages are, alas, not widely disseminated. Moreover, they do
not inherently offer control over the exploitation of parallelism. An alternative exists in
CSP languages. These languages require some infrastructure to be built to handle the
scheduling. Once that is done, it becomes possible to write more efficient (though more
verbose) code in CSP style. We shall return to this comparison at the end of the paper.
In the meantime we proceed to an implementation in a channel-based CSP language.



Translation to newsqueak

      Pike's newsqueak language supports asynchronous processes communicating over
typed channels. Channels and functions are first-class citizens. Values of both kinds
may be assigned to variables, passed as arguments, and returned by functions. Processes
and channels may be created at will. There is just one kind of interprocess
communication--rendezvous of two processes on a channel. Consequently the queuing
and fanout that are necessary for power series computations must be programmed
explicitly. A fair choice is made among competing rendezvous. When not waiting for
rendezvous, all processes progress, at dithered rates, a few atomic operations per
interpreter cycle.

      Most of newsqueak will be familiar to experienced programmers. The
communication semantics descends mainly from CSP, 1 the syntax from C, and the type
system from Pascal and ML. 8 Unusual features will be described as they are needed. For
details see the accompanying paper and the manual. 3, 9

      At base we need rationals, declared as numerator-denominator pairs, and power
series, declared as channels that carry rationals,

     # rationals and power series - basic declarations
     type rat: struct of { num: int; den: int; };


     type ps: chan of rat;

along with some obvious support routines, declared as follows. The literal definitions of
these functions, all one-liners, have been left out, simply to avoid discussing the finicky,
but straightforward, matter of data constructors.
                                                 -9-


     ratmk: prog(i:int, j:int) of rat;                                    #   construct a rational i/j
     ratadd: prog(r:rat, s:rat) of rat;                                   #   return sum of rationals
     ratsub: prog(r:rat, s:rat) of rat;                                   #   subtract
     ratmul: prog(r:rat, s:rat) of rat;                                   #   multiply
     ratprint: prog(r:rat);                                               #   print a rational


     psmk: prog() of ps;                                                  # construct a power series

     A program to print (endlessly) a power series is easy. To write it, though, we need
some peculiar newsqueak syntax. Beware, in newsqueak := is not a simple assignment.
It declares and initializes a new variable of a type inferred from the initializing
expression. Assignment is represented by = as in C. The prefix operator