Migrate *the rest* of the math blogs

1 files changed, 307 insertions, 0 deletions

diff --git a/blog/conseq.md b/blog/conseq.md
new file mode 100644
index 0000000..eadc970
--- /dev/null
+++ b/blog/conseq.md
@@ -0,0 +1,307 @@
++++
+rss = "SICP subsection 3.5.2 in Python"
+date = Date(2019, 2, 28)
++++
+@def tags = ["sicp", "fun", "python", "calculus"]
+
+# Infinite Sequences: A Case Study in Functional Python
+
+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
+```
+
+@@colbox-blue
+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, this 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