FriCAS Code

This text explains how stmult routine in 'sttaylor.spad' works.  I think
the best explanation is to give short reconstruction how the routine was
written.

1) Let us look first at general pattern of stream routines.  Good
example is the simplified version of 'map':

    map(f,x)== delay
      empty? x => empty()
      concat(f frst x, map(f,rst x))

Note that the result is produced by the 'delay' function.  'delay(f)'
takes function 'f' as argument and produces from it a stream, such
that evaluation of the stream causes f to be called.  The name 'delay'
is related to the fact that 'delay' effectively delays evaluation
of its argument to the moment when first value from the stream
is requested.  When stream is evaluated?  Evaluation is triggered
by call to 'empty?'.  In other words, before trying to access
first element of the stream we need to call 'empty?' to cause
evaluation, after the call we know if the element is there (and
if it is present we can use it).  So above the first line calls
'empty?' which either give us true, in which case we simply return
empty stream or we may use the first element of input stream to
produce first element of the result result.  Now, once we have
first element of the result in hand we can produce the result by
concatenating first element of the result with the rest of the
result.  In case of 'map' the result is obtained by recursive
application of 'map' to the rest of input stream.

Important remark: stream are produced and should be consumed
in lazy way.  'delay' makes sure that no computation takes
place before first element of the result is requested.  However
for single request map should produce only one element and to
produce it should consume only one element form input stream.
Using 'concat' makes sure that caller will get the computed
element and that request for further elements will go to 'map'
applied to the rest of stream.

Let me also say that the code above looks like unbounded
recursion, but in fact there is _no_ recursion at runtime:
what looks like recursive call to 'map' is done _after_
call to outer map has finished.

2) Having pattern for general operation, we now apply it in a
concrete situation.  For example, we would like to multiply
element 'n' of the stream by 'n+1'.  To do this we need
to know what 'n' is, so we need to pass 'n' between invocations
of function.  So we need an auxiliary function and a wrapper:

   stmuln_aux(n0, x) == delay
       empty? x => empty()
       n := n0 + 1
       firstval := n*frst x
       concat(firstval, stmuln_aux(n, rst x))

   stmuln(x) == stmuln_aux(0, x)

Note: due to stupidity of Spad compiler I can not simply increment
'n0', but must introduce extra variable 'n'.

3) Back to multiplication.  For multiplication formula is:

   (f*g)(n) = \sum_{i=0}^{n} f(i)*g(n-i)

This is more complicated then previous case because I need terms
f(0),...,f(n) and g(0),...,g(n) to produce the result.
Naive implementation (assuming infinite streams) may look like:

   stmult(n, x, y) == delay
       valn := sum(coefficient(x, i)*coefficient(y, n - i), i=0..n)
       concat(valn, stmult(n + 1, x, y))

   mult(x, y) == stmult(0, x, y)

However, 'coefficient' for streams has to step trough stream from
the start up to requested element, so it is slow if 'n' gets large.
It is easy to see that we can step trough 'x' as we sum saving
half of the operations:

   stmult(n, x, y) == delay
       valn : A := 0
       xp := x
       for i in 0..n repeat
           empty? xp => break
           valn := frst xp * coefficient(y, n - i)
           xp := rst xp
       concat(valn, stmult(n + 1, x, y))

We still use 'coefficient' on 'y'.  Since we need to go trough 'y'
backwards, it is a bit harder to get rid of.  Easy way it to
have a list which stores the first 'n' elements of 'y' in backward
order:

   stmult(n, x, y0, ll0) == delay
       yval : A
       y := y0
       if empty? y then
           yval := 0
       else
           yval := frst y
           y := rst y
       ll := cons(yval, ll0)
       xp := x
       llp := ll
       valn : A := 0
       for i in 0..n repeat
           empty? xp => break
           valn := valn + frst xp * first(llp)
           llp := rest(llp)
           xp := rst xp
       concat(valn, stmult(n + 1, x, y, ll))

This may be (I did not test it) correct routine to multiply infinite
power series, but it does not conform to expectation of rest of FriCAS.
Namely, if we are given to finite series it will produce infinite
series -- after finite number of nonzero terms we will get infinite
stream of zeros.  There is also unneeded eagerness: if first element
if 'y' is zero we do not need element 'n' of 'x' to produce first element
of the result.  Also, if 'y' is finite there is no need to append
zeros.  First, let us handle leading zeros of 'y'.  If first
element of 'y' is zero we get zero as first element of the result and
rest of result is 'x' times rest of 'y':

     ...
     yval := frst y
     n = 0 and yval = 0 => concat(0, stmult(n, x, rst y, ll))

Note: we need to test if 'n' is 0 to know that we really have first
element and not some zero in the middle.

Next, if we get to the end of 'y' instead of extending 'll' by zero
we may drop first element of 'x'.  We need to be a little
careful if 'x' is also empty (possibly after dropping several
elements).  And we need to properly adjust 'n'.  Adjusting 'n' is
easier if we start from 'n = -1' instead of 'n = 0'.  Also, if at start
'y' is empty then the result is 0 (empty):

       if empty? y then
            n < 0 => return empty()
            empty? x => return empty()
            x := rst x

If in the loop we discover that 'x' is empty we can also
return 0.  Finally, collecting the changes above we get:

   stmult(n0, x0, y0, ll0) == delay
       x := x0
       y := y0
       ll := ll0
       n = n0
       if empty? y then
           n < 0 => return empty()
           empty? x => return empty()
           x := rst x
       else
           yval := frst y
           n < 0 and yval = 0 => concat(0, stmult(n, x, rst y, ll))
           y := rst y
           ll := cons(yval, ll)
           n := n + 1
       xp := x
       llp := ll
       valn : A := 0
       for i in 0..n repeat
           empty? xp =>
               i = 0 => return empty()
               break
           valn := valn + frst xp * first(llp)
           llp := rest(llp)
           xp := rst xp
       concat(valn, stmult(n + 1, x, y, ll))

This is almost the routine I want to use.  However, there is one
extra problem: if series has unevaluated tail printing routines
print 'O(..)' term even if after further evaluation it turns out
that in fact the tail is empty.  To avoid such artifacts printing
routines check if tail is explicitly empty.  Note that testing
if tail is empty causes evaluation and printing routines are not
allowed to evaluate extra terms.  But 'explicitlyEmpty?' does not
cause evaluation, even if 'explicitlyEmpty?' return false
'empty?' may return 'true'.  So, to make printing routines
happy we need to make sure that the tail is explicitly empty.
This is done adding after loop test:

         ....
      explicitlyEmpty? rst x and
        explicitlyEmpty? y => concat(res, empty())

Note that if 'x' is empty we call return inside loop, so after
loop 'x' is always nonempty and calling 'rst' is safe.  Note that
if 'x' has only one element and 'y' is empty then in the next iteration
we will throw out the only element of 'x' and get empty result.
So the test allows us to produce earlier the empty result.
OTOH if that test fails in principle we can get more elements,
so must use general form.

Martin also asked why the new version uses less space than the
old one and why the old one used so much space.

First, note that stream data is stored in nested anonymous
functions.  The new version creates single function for
each term of the series, so when we want n terms n such
functions are created.  Each function takes fixed amount
of space for variables (the space is actually taken by "frozen"
values of variables from outside which the function uses).
The variables may reference large data structure (for
example x and y may be big), but fortunately such large
structures may be shared.  In particular parts of 'x' and 'y' are
shared so that at each step we add (compute) only one extra element
of 'x' and of 'y'.  The list 'll' has in general 'n + 1' elements, but
we share 'n' elements with previous iteration.  So at each step we
allocate fixed amount of memory.  More precisely, due to calls to
'empty?' element 'n' of 'x' and element 'n' of 'y' are computed
and stored.  'stmult' stores element 'n' of 'x*y' plus fixed storage
(of order 10-15 machine words).  In other words, except for
original arguments and the final result memory usage is linear in n.

The old version was:

   (x:ST A) * (y:ST A) == delay
     empty? y => zro()
     empty? x => zro()
     concat(frst x * frst y,frst x * rst y + rst x * y)

It is easy to see that this version allocates stream node for
each term in the product, that is of order n^2 nodes.  Theoretically
most of those nodes may be garbage collected.  I did not check
if the FriCAS code for some reason (for example due to keeping
unnecessary extra reference to the data) did not allow collection,
or if this was problem with SBCL garbage collector, but the fact
was that the space usage grow quadratically...

The new version allocates much less memory, so there is no
chance for garbage collector (or FriCAS programmers) to screw up.