A few years ago, Pleng showed me a really cool trick you can do with lazy infinite lists in Haskell.

Kublai: Oh no, not Haskell! 😱

Don’t worry. There will be no mandatory Haskell in this post (although it is a really nice language which you should check out!). Instead, today Erik asked me if that trick also works using infinite Python generators, and turns out it does!

First, as a warm-up, let’s define an infinite generator for the positive integers.

Kublai: Oh, oh, I know how to do this! That’s easy!

def ints():
    cnt = 1
    while True:
        yield cnt
        cnt += 1

Nice! That works. However, we’re going to define it in a slightly different way. Recursively! The commented line is Haskell and below that is the translation into Python.

# ints = 1 : map (+1) ints

def ints():
    yield 1
    yield from map(lambda x: x + 1, ints())

Kublai: WTF? How does that work?

It’s actually pretty simple. The first positive integer is 1, obviously. The remaining integers after 1 are just the positive integers but with 1 added to each one. As simple as that.

We’ll need a function for printing out the first few elements of that generator. (Assuming you already imported the correct modules.)

# take n f

def take(n, f):
    return list(itertools.islice(f, n))

And now let’s test it with ints().

>>> take(10, ints())
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

Yay! It works!

Next, let’s define a function to integrate a Taylor series, where the coefficients of the Taylor series will be expressed as an infinite generator. Recall that the integral of $x^n$ is $\frac{1}{n+1}x^{n+1}+C$. That means our integral is some leading constant plus all the terms of the original Taylor series shifted over by one and divided by $n+1$.

# integrate f c = c : zipWith (/) f ints

def integrate(f, c):
    yield c
    yield from map(operator.truediv, f, ints())

Now time for some magic. Recall that the derivative of $e^x$ is itself and $e^x = \int_0^x e^y dy + 1$ so the integration constant is 1. So, let’s use these properties to define $e^x$!

# expSeries = integrate expSeries 1

def expSeries():
    yield from integrate(expSeries(), 1)

Let’s print it out.

>>> take(10, expSeries())
[1, 1.0, 0.5, 0.16666666666666666, 0.041666666666666664, 0.008333333333333333, 0.001388888888888889, 0.0001984126984126984, 2.48015873015873e-05, 2.7557319223985893e-06]

Whoa! That’s the Taylor series for $e^x$!

Kublai: Wait, we never even told Python what $e^x$ is, other than that it’s equal to its own integral! What kind of crazy magic is this?

We can also evaluate $e^x$ using that Taylor series.

# evalAt n f x = foldr (\a acc -> a + acc * x) 0 (take n f)

def evalAt(n, f, x):
    return functools.reduce(lambda acc, a: a + acc * x, reversed(take(n, f)))

Let’s compare it to the one in the Python standard library.

>>> evalAt(100, expSeries(), 2)
7.38905609893065
>>> math.exp(2)
7.38905609893065

Pretty close!

Here’s the punchline. Recall that $\sin x = \int_0^x \cos y dy$ and $\cos x = -\int_0^x \sin y dy + 1$. Let’s encode that in Python!

# sine = integrate cosine 0
# cosine = map negate (integrate sine -1)

def sine():
    yield from integrate(cosine(), 0)

def cosine():
    yield from map(operator.neg, integrate(sine(), -1))

Now let’s evaluate sine:

>>> evalAt(100, sine(), 2)
0.9092974268256817
>>> math.sin(2)
0.9092974268256817

Kublai: 🤯

Pleng noticed that my Python code can be pretty slow since Python generators aren’t memoized unlike Haskell lists (and also hits the Python recursion limit).

Fortunately, someone on StackOverflow has already solved this problem using a decorator, so here’s the fix:

def memoize(f):
    cache = {}
    @functools.wraps(f)
    def wrapper(*args, **kwargs):
        k = args, frozenset(kwargs.items())
        it = cache[k] if k in cache else f(*args, **kwargs)
        cache[k], result = itertools.tee(it)
        return result
    return wrapper

@memoize # Add this line
def ints():
    yield 1
    yield from map(lambda x: x + 1, ints())

@memoize # Add this line
def expSeries():
    yield from integrate(expSeries(), 1)

print(take(expSeries(), 1000))

And now it’s fast!

Also, Isaac noted that this trick also works using Python’s fractions module to get the answer as a rational number instead of a float.