On Removing Singularities from Rational Functions

Introduction
Suppose we have the function

> f x = 1/(x+x^2) - 1/(x+2*x^2)

Some basic algebraic manipulation shows that in the limit as x→0, f(x)→1. But we can't simply compute f 0 because this computation involves division by zero at intermediate stages. How can we automate the process of computing the limit without implementing symbolic algebra?

I've already described one way to remove singularities from a function. But that approach is very limited in its applicability.

This article is about a variation on the approach to formal power series that nicely showcases some advantages of lazy lists. It will allow us to form Laurent series of functions so we can keep track of the singularities.

The usual Haskell approach to power series allows you to examine the coefficients of any term in the power series of the functions you can form. These series can't be used, however, to evaluate the function. Doing so requires summing an infinite series, but we can't do so reliably because no matter how many terms in a power series we add, we can never be sure that there aren't more large terms further downstream that we haven't reached yet. And if we want to perform computations completely over the rationals, say, we don't want to be dealing with infinite sums.

I'd like to look at a way of working with power series that allows us to perform exact computations making it possible to answer questions like "what is the sum of all the terms in this power series starting with the x^n term?" By extending to Laurent series, and implementing the ability to selectively sum over just the terms with non-negative powers, we can compute functions like f above at 0 and simply skip over the troublesome http://en.wikipedia.org/wiki/Pole_%28complex_analysis%29.

Power Series
When I previously discussed power series I used code that worked with the coefficients of the power series. This time we want to work with values of the function so it makes sense to store, not the coefficients a_i but the terms themselves, a_ixⁱ. So instead of a list of coefficients, Num a => [a] we need a representation that looks a little like:

> data Power a = Power (a -> [a])

where we pass x in as an argument to the function contained in a Power. But we also want to allow Laurent series so we need to also store an offset to say which (possibly negative) term our series starts with:

> data Laurent a = Laurent (a -> (Int, [a]))

But this fails us for at least two reasons:

1. We have the individual terms, but to evaluate the function requires summing all of the terms in an infinite list.
2. If we have a Laurent series, then we need to store values of a_ixⁱ for x=0 and i<0. We'll end up with division by zero errors.

Partial Sum Series
So here's what we'll do instead. Suppose our power series is Σ_i=n^∞a_ixⁱ. We'll store the terms s_j=Σ_i=j^∞a_ix^i-j. Our type will look like:

> data Partial a = Partial (a -> (Int, [a]))

> instance Eq (Partial a)
> instance Show (Partial a)

It's straightforward to add two functions in this form. We just add them term by term after first aligning them so that the xⁱ term in one is lined up with the xⁱ term in the other:

> instance Num a => Num (Partial a) where

>  Partial f + Partial g = Partial $ \x ->
>     let (m, xs) = f x
>         (n, ys) = g x
>         pad 0 _ ys = ys
>         pad n x ys = let z:zs = pad (n-1) x ys
>                      in x*z : z : zs
>         l = min m n
>     in (l, zipWith (+) (pad (m-l) x xs) (pad (n-l) x ys))

Notice the slight subtlety in the alignment routine pad. By the definition above, the jth term has a factor of x^j built into it. So we need to multiply by x each time we pad our list on the left.

Now we need to multiply series. We know from ordinary power series that we need some sort of convolution. But it looks like for this case we have an extra complication. We appear to need to difference our representation to get back the original terms, convolve, and then resum. Amazingly, we don't need to do this at all. We can convolve 'in place' so to speak.

Here's what an ordinary convolution looks like when we want to multiply the sequence of terms (a_i) by (b_i):

In this example, the blue diagonal corresponds to the terms that are summed to get the 4th term in the result.

However, we wish to work with partial sums s_j=Σ_i=j^∞a_ix^i-j and t_j=Σ_i=j^∞t_ix^i-j, constructing the partial sums of the convolution of a and b from s and t. The partial sums of the convolution can be derived from the partial sums by tweaking the convolution so it looks like this:

The blue terms work just like before and need to be summed. But we also need to subtract off the red terms, weighted by a factor of x. That's it! (I'll leave that as an exercise to prove. The inclusion-exclusion principle helps.)

The neat thing is that the red terms for each sum are a subset of the blue terms needed for the next element. We don't need to perform two separate sums. We can share much of the computation between the red and blue terms. All we need to do is write an ordinary convolution routine that additionally returns not just the blue terms, but a pair containing the blue sum and the red sum.

>  Partial f * Partial g = Partial $ \x ->
>     let (m, xs) = f x
>         (n, ys) = g x
>         (outer, inner) = convolve xs ys
>         f' a b = a-x*b -- (the subtraction I mentioned above)
>     in (m+n, zipWith f' outer inner)
>  fromInteger n = let n' = fromInteger n in Partial $ \_ -> (0, n' : repeat 0)

>  negate (Partial f) = Partial $ \x -> let (m, xs) = f x
>                                 in (m, map negate xs)
>  signum = error "signum not implemented"
>  abs    = error "signum not implemented"

This is an ordinary convolution routine tweaked to return the partial sum inner.

> convolve (a0:ar@(a1:as)) ~(b0:br@(b1:bs)) =
>  let (inner, _) = convolve ar br
>      ab = map (a0 *) bs
>      ba = map (* b0) as
>  in (a0*b0 : a0*b1+a1*b0
>            : zipWith3 (\a b c -> a+b+c) inner ab ba, 0 : inner)

The code is very similar to the usual power series multiplication routine. We can also use the same method described by McIlroy to divide our series.

As our series are a munged up version of the usual power series it's pretty surprising that it's possible to divide with so little code:

> instance Fractional a => Fractional (Partial a) where
>   fromRational n = let n' = fromRational n in Partial $ \_ -> (0, n' : repeat 0)
>   recip (Partial f) = Partial $ \x ->
>      let nibble (n, x:xs) | x==0      = nibble (n+1, xs)
>                           | otherwise = (n, (x:xs))
>          (n, xs) = nibble (f x)
>      in (-n, rconvolve x xs)

In effect, rconvolve solves the equation convolve a b==1:

> rconvolve x (a0:ar@(a1:as)) =
>   let (outer, inner) = convolve ar result
>       f a b = x*b-a
>       r = -1/f a0 a1
>       result = recip a0 : (map (r *) $ zipWith f outer inner)
>   in result

Note one ugly quirk of this code. I need to 'nibble' off leading zeroes from the series. This requires our underlying type a to have computable equality. (In principle we can work around this using parallel or).

That's it. We can now write a function to compute the positive part of a rational function. (By positive part, I mean all of the terms using non-negative powers of x.)

> pos f z = let Partial g = f $ Partial $ \x -> (1, 1 : repeat 0)
>               (n, xs) = g z
>           in if n>0
>               then z^n*head xs
>               else xs!!(-n)

Here are some examples:

> test1 = let f x = (1+2*x)/(3-4*x*x)
>         in pos (\x -> 1/(f x-f 0)/x) (0::Rational)

> test2 = pos (\x -> 1/(1+4*x+3*x^2+x^3) - 1/(1+x)) (1::Rational)

The original example I started with:

> test3 = pos (\x -> 1/(x+x^2) - 1/(x+2*x^2)) (0::Rational)

No division by zero anywhere!

Conclusions
The code works. But it does have limitations. As written it only supports rational functions. It's not hard to extend to square roots. (Try writing the code - it makes a nice exercise.) Unfortunately, any implementation of square root will (I think) require a division by x. This means that you'll be able to compute the positive part away from zero, but not at zero.

This method can't be extended fully to transcendental functions. But it is possible to add partial support for them. In fact, So with a little work we can still compute the positive part of functions like 1/sqrt(cos x-1) away from x==0. But applying cos to an arbitrary rational function may need more complex methods. I encourage you to experiment.

Note that this code makes good use of laziness. If your function has no singularities then you might find it performs no computations beyond what is required to compute the ordinary numerical value.