about summary refs log tree commit diff homepage
path: root/blog/conseq.md
blob: 92c553228808fb1e813dd0c8529886253fb4e4d6 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
+++
title = "Infinite Sequences: A Case Study in Functional Python"
rss = "SICP subsection 3.5.2 in Python"
date = Date(2019, 2, 28)
tags = ["fun", "math", "python", "sicp"]
+++

# {{title}}

In this article, we will only consider sequences defined by a function
whose domain is a subset of the set of all integers.  Such sequences will be
*visualized*, i.e. we will try to evaluate the first few (thousand) elements,
using functional programming paradigm, where functions are more similar
to the ones in math (in contrast to imperative style with side effects
confusing to inexperenced coders).  The idea is taken from [subsection 3.5.2
of SICP] and adapted to Python, which, compare to Scheme, is significantly
more popular: Python is pre-installed on almost every modern Unix-like system,
namely macOS, GNU/Linux and the \*BSDs; and even at MIT, the new 6.01 in Python
has recently replaced the legendary 6.001 (SICP).

One notable advantage of using Python is its huge **standard** library.
For example the *identity sequence* (sequence defined by the identity function)
can be imported directly from ``itertools``:

```python
>>> from itertools import count
>>> positive_integers = count(start=1)
>>> next(positive_integers)
1
>>> next(positive_integers)
2
>>> for _ in range(4): next(positive_integers)
... 
3
4
5
6
```

To open a Python emulator, simply lauch your terminal and run `python`.
If that is somehow still too struggling, navigate to [the interactive shell]
on Python.org.

*Let's get it started* with somethings everyone hates: recursively defined
sequences, e.g. the famous Fibonacci ($F_n = F_{n-1} + F_{n-2}$,
$F_1 = 1$ and $F_0 = 0$).  Since [Python does not support] [tail recursion],
it's generally **not** a good idea to define anything recursively (which is,
ironically, the only trivial *functional* solution in this case)
but since we will only evaluate the first few terms
(use the **Tab** key to indent the line when needed):

```python
>>> def fibonacci(n, a=0, b=1):
...     # To avoid making the code look complicated,
...     # n < 0 is not handled here.
...     return a if n == 0 else fibonacci(n - 1, b, a + b)
... 
>>> fibo_seq = (fibonacci(n) for n in count(start=0))
>>> for _ in range(7): next(fibo_seq)
... 
0
1
1
2
3
5
8
```

!!! note "Note"

    The `fibo_seq` above is just to demonstrate how `itertools.count`
    can be use to create an infinite sequence defined by a function.
    For better performance, the following should be used instead:

    ```python
    def fibonacci_sequence(a=0, b=1):
        yield a
        yield from fibonacci_sequence(b, a+b)
    ```

It is noticable that the elements having been iterated through (using `next`)
will disappear forever in the void (oh no!), but that is the cost we are
willing to pay to save some memory, especially when we need to evaluate a
member of (arbitrarily) large index to estimate the sequence's limit.
One case in point is estimating a definite integral using [left Riemann sum].

```python
def integral(f, a, b):
    def left_riemann_sum(n):
        dx = (b-a) / n
        def x(i): return a + i*dx
        return sum(f(x(i)) for i in range(n)) * dx
    return left_riemann_sum
```

The function `integral(f, a, b)` as defined above returns a function taking
$n$ as an argument.  As $n\to\infty$, its result approaches
$\int_a^b f(x)\mathrm d x$.  For example, we are going to estimate
$\pi$ as the area of a semicircle whose radius is $\sqrt 2$:

```python
>>> from math import sqrt
>>> def semicircle(x): return sqrt(abs(2 - x*x))
... 
>>> pi = integral(semicircle, -sqrt(2), sqrt(2))
>>> pi_seq = (pi(n) for n in count(start=2))
>>> for _ in range(3): next(pi_seq)
... 
2.000000029802323
2.514157464087051
2.7320508224700384
```

Whilst the first few aren't quite close, at index around 1000,
the result is somewhat acceptable:

```
3.1414873191059525
3.1414874770617427
3.1414876346231577
```

Since we are comfortable with sequence of sums, let's move on to sums of
a sequence, which are called series.  For estimation, again, we are going to
make use of infinite sequences of partial sums, which are implemented as
`itertools.accumulate` by thoughtful Python developers.  [Geometric] and
[p-series] can be defined as follow:

```python
from itertools import accumulate as partial_sums

def geometric_series(r, a=1):
    return partial_sums(a*r**n for n in count(0))

def p_series(p):
    return partial_sums(1 / n**p for n in count(1))
```

We can then use these to determine whether a series is convergent or divergent.
For instance, one can easily verify that the $p$-series with $p = 2$
converges to $\pi^2 / 6 \approx 1.6449340668482264$ via

```python
>>> s = p_series(p=2)
>>> for _ in range(11): next(s)
... 
1.0
1.25
1.3611111111111112
1.4236111111111112
1.4636111111111112
1.4913888888888889
1.511797052154195
1.527422052154195
1.5397677311665408
1.5497677311665408
1.558032193976458
```

We can observe that it takes quite a lot of steps to get the precision we would
generally expect ($s_{11}$ is only precise to the first decimal place;
second decimal places: $s_{101}$; third: $s_{2304}$).
Luckily, many techniques for series acceleration are available.
[Shanks transformation] for instance, can be implemented as follow:

```python
from itertools import islice, tee

def shanks(seq):
    return map(lambda x, y, z: (x*z - y*y) / (x + z - y*2),
               *(islice(t, i, None) for i, t in enumerate(tee(seq, 3))))
```

In the code above, `lambda x, y, z: (x*z - y*y) / (x + z - y*2)` denotes
the anonymous function $(x, y, z) \mapsto \frac{xz - y^2}{x + z - 2y}$
and `map` is a higher order function applying that function to
respective elements of subsequences starting from index 1, 2 and 3 of `seq`.
On Python 2, one should import `imap` from `itertools` to get the same
[lazy] behavior of `map` on Python 3.

```python
>>> s = shanks(p_series(2))
>>> for _ in range(10): next(s)
... 
1.4500000000000002
1.503968253968257
1.53472222222223
1.5545202020202133
1.5683119658120213
1.57846371882088
1.5862455815659202
1.5923993101138652
1.5973867787856946
1.6015104548459742
```

The result was quite satisfying, yet we can do one step futher
by continuously applying the transformation to the sequence:

```python
>>> def compose(transform, seq):
... 	yield next(seq)
... 	yield from compose(transform, transform(seq))
... 
>>> s = compose(shanks, p_series(2))
>>> for _ in range(10): next(s)
... 
1.0
1.503968253968257
1.5999812811165188
1.6284732442271674
1.6384666832276524
1.642311342667821
1.6425249569252578
1.640277484549416
1.6415443295058203
1.642038043478661
```

Shanks transformation works on every sequence (not just sequences of
partial sums).  Back to previous example of using left Riemann sum
to compute definite integral:

```python
>>> pi_seq = compose(shanks, map(pi, count(2)))
>>> for _ in range(10): next(pi_seq)
... 
2.000000029802323
2.978391111182236
3.105916845397819
3.1323116570377185
3.1389379264270736
3.140788413965646
3.140921512857936
3.1400282163913436
3.1400874774021816
3.1407097229603256
>>> next(islice(pi_seq, 300, None))
3.1415061302492413
```

Now having series defined, let's see if we can learn anything
about power series. Sequence of partial sums of power series
$\sum c_n (x - a)^n$ can be defined as

```python
from operator import mul

def power_series(c, start=0, a=0):
    return lambda x: partial_sums(map(mul, c, (x**n for n in count(start))))
```

We can use this to compute functions that can be written as
[Taylor series][]:

```python
from math import factorial
def exp(x):
    return power_series(1/factorial(n) for n in count(0))(x)

def cos(x):
    c = ((1 - n%2) * (1 - n%4) / factorial(n) for n in count(0))
    return power_series(c)(x)

def sin(x):
    c = (n%2 * (2 - n%4) / factorial(n) for n in count(1))
    return power_series(c, start=1)(x)
```

Amazing!  Let's test 'em!

```python
>>> e = compose(shanks, exp(1)) # this should converges to 2.718281828459045
>>> for _ in range(4): next(e)
... 
1.0
2.749999999999996
2.718276515152136
2.718281825486623
```

Impressive, huh? For sine and cosine, series acceleration is not even necessary:

```python
>>> from math import pi as PI
>>> s = sin(PI/6)
>>> for _ in range(5): next(s)
... 
0.5235987755982988
0.5235987755982988
0.49967417939436376
0.49967417939436376
0.5000021325887924
>>> next(islice(cos(PI/3), 8, None))
0.500000433432915
```

[subsection 3.5.2 of SICP]: https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-24.html#%_sec_3.5.2
[the interactive shell]: https://www.python.org/shell
[Python does not support]: http://neopythonic.blogspot.com/2009/04/final-words-on-tail-calls.html
[tail recursion]: https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-11.html#call_footnote_Temp_48
[left Riemann sum]: https://en.wikipedia.org/wiki/Riemann_sum#Left_Riemann_sum
[Geometric]: https://en.wikipedia.org/wiki/Geometric_series
[p-series]: https://math.oregonstate.edu/home/programs/undergrad/CalculusQuestStudyGuides/SandS/SeriesTests/p-series.html
[Shanks transformation]: https://en.wikipedia.org/wiki/Shanks_transformation
[lazy]: https://en.wikipedia.org/wiki/Lazy_evaluation
[Taylor series]: https://en.wikipedia.org/wiki/Taylor_series