Nothing says fun like formal power series!

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

1 + x + x² + x³ + x⁴ + …

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).

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² 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

p0q0 + (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

p0 + 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)