Nothing says fun like formal power series!

A formal power series is a (usually) infinite polynomial in x. For example

*1 + x + x*^{2}* + x*^{3}* + x*^{4}* + …*

This is an expression, not a number. If we give a value to x, we may get a number, or the evaluation may run away (diverge) on us (if |x|≥1).

They’re called *formal* power series because we don’t normally try to evaluate them, we just manipulate them formally and symbolically.

For example, if we multiply the above series by *1-x*, we get … *1*. So *1-x* is the multiplicative inverse of that series, even though the series diverges for most values of *x*. These are purely formal manipulations.

So where’s the fun? Well recently my colleagues and I produced a Lucid interpreter. It’s written in Python and handles the full language described in the Lucid book, plus it supports an extra space dimension. I’ll tell you about it in another post and make it available through github.

Anyway pyLucid plus formal power series spells fun. Clearly a formal power series is determined by the infinite sequence of its coefficients, and such a sequence can be represented in a straightforward way as a space vector. Thus the power series 1 is *1 sby 0*, the series *x* is *0 sby 1 sby 0*, the series *x*^{2} is *0 sby 0 sby 1 sby 0 *and the power series first given is simply *1*.

More precisely 1 sby 0 is the vector

1, 0, 0, 0, …

which represents

1 + 0x + 0x^2+ 0x^3 + …

0 sby 1 sby 0 is the vector

0, 1, 0, 0, 0, …

which represents the series

0 + 1x + 0x^2+ 0x^3 + …

0 sby 0 sby 1 sby 0 represents

0 + 0x + 1x^2+ 0x^3 + …

and 1 is the vector

1, 1, 1, 1, …

which represents the original series above.

What about operations on power series? Addition is just that: if pyLucid vectors P and Q represent two power series, P+Q represents their sum, since addition is coefficientwise.Multiplication is more complex.

(I’m going using ASCII notation instead of proper sub and super scripting. I tried but the cockamamy new wordpress editor ate my html).

Anyway we can break down P to p0+xP’ where P’ is p1+p2x+p3x^2+… . Similarly Q is q0+xQ’ and multiplying them we get

p*0q0 + (p0q1+p1q0)x + P’Q’x^2*

This defines multiplication of P and Q in terms of multiplication of P’ and Q’. This is recursion but P’ and Q’ are in no sense simpler than P and Q. However we can use it to write pyLucid code because it produces two coefficients before it recurses. It doesn’t deadlock.

The pyLucid definition of product of two series represented as space vectors is

*pprod(p,q) = r*

where

p0 = init p; q0 = init q;

pp = succ p; qp = succ q;

r0 = p0*q0; r = (r0 sby (p0*qp + q0*pp)) + (0 sby 0 sby pprod(pp,qp));

end;

(Dangnabbit, there doesn’t seem to be a way to indent a block in Gutenberg, WordPress’ newfangled editor.)

We can test *pprod* by multiplying 1 (our original power series) by itself, and we get

*1 + 2x + 3x^2 + 4x^3 + …*

which is correct. In other words, *pprod(1,1)* is

*1, 2, 3, 4, …*

Division is even trickier, I’ll spare you the explanation, the code is

*pdiv(q,w) = t*

where

q0 = init q;

w0 = init w;

r = q0/w0;

v = succ(q – r*w);

t = r sby pdiv(v,w);

end;

We can test it by setting *one = 1 sby 0* and calculating *pdiv(one,1)* which gives coefficients

1, -1, 0, 0, 0, …

which is correct since the result is *1-x.*

Formal power series can be integrated and differentiated (formally) and that’s where things get interesting. The derivative of

*p0 + p1x + p2x^2 + p3x^3 + …*

is

*p1 + 2p2x + 3p3x^2 + …*

and this is easy to code:

*pderiv(g) = h*

where

gp = succ g;

k = 1 sby k+1;

h = gp*k;

end;

Integration is even simpler, the integral of

p*0 + p1x + p2x^2 + p3x^3 + …*

is

*c + p0x + p1x^2/2 + p2x^3/3 + …*

(c is the constant of integration.)

The code is

*pinteg(c,s) = d where*

i = 1 sby i+1;

d = c sby s/i;

end;

Simple enough … but now let’s think about the power series corresponding to the exponential function e^x. This function is its own derivative and therefore also its own integral. More precisely, e^x is one plus the integral of e^x. In code, we get

*ex = pinteg(1,ex)*

and this works! Even though it’s recursive, we get all the coefficients:

1.00000 1.00000 0.50000 0.16667 0.04167 0.00833 0.00139 0.00020 …

Let’s compute e! By evaluating ex at x=1. Evaluating doesn’t always work but sometimes we can get approximations. The code

*peval(p,a) = v*

where

v = init p +(0 sby a * peval(succ p, a));

end;

gives us the sequence of partial sums, which sometimes converges. Sure enough, *peval(ex,1)* yields

1.00000 2.00000 2.50000 2.66667 2.70833 2.71667 2.71806 2.71825 2.71828 …

In the same way, sin(x) and cos(x) are each others derivatives/integrals and the pyLucid code

*sinx = pinteg(0,cosx)*

cosx = pinteg(1,sinx)

works. For example, the coefficients of *sinx* are

0.00000 1.00000 0.00000 -0.16667 0.00000 0.00833 0.00000 -0.00020 …

Now for the fireworks! Arctan(x) is really interesting and its derivative is 1/(1+x^2). We can integrate this in code:

*x2 = 0 sby 0 sby 1 sby 0; one = 1 sby 0;*

atanx = pinteg(0,pdiv(one,one+x2));

This does the job, giving the coefficients of arctan(x) as

0.00000 1.00000 0.00000 -0.33333 0.00000 0.20000 0.00000 -0.14286 …

Now it so happens that the arctan of 1/2 plus the arctan of 1/3 is pi/4. So the expression

*4*(peval(atanx,1/2)+peval(atanx,1/3))*

is what we want, and sure enough we get

0.00000 3.33333 3.33333 3.11728 3.11728 3.14558 3.14558 3.14085 3.14085 3.14174 3.14174 3.14156 3.14156 3.14160 3.14160 3.14159 …

The amazing thing about all this is that pi emerges from a relatively simple program (given in full below) in which only small integers appear.

Enough fun for now. I’ll leave you with the complete program and a promise to tell you more about the interpreter and to make it publicly available.

4*(peval(atanx,1/2)+peval(atanx,1/3))

where

x2 = 0 sby 0 sby 1 sby 0;

atanx = pinteg(0,pdiv(one,one+x2));

one = 1 sby 0;

peval(p,a) = v
where
v = init p +(0 sby a * peval(succ p, a));
end;
pdiv(q,w) = t
where
q0 = init q;
w0 = init w;
r = q0/w0;
v = succ(q - r*w);
t = r sby pdiv(v,w);
end;
pinteg(c,s) = d where
i = 1 sby i+1;
d = c sby s/i;
end;
columns = 16;
rows = 1;
numformat = '%7.5f';

end

(Sorry this is the best I can do with the persnickety editor which even screws up a simple paste. Land o’ Goshen)