Python for Everyone!
Oregon Curriculum Network
There's a real difference between verifying a suspected identity, in this case involving arctangents and phi ($\phi$), and providing a proof.
My impression is some of the Ramanujan convergent series have been verified, which is easier to do through brute force computation, rather than formally proved. Am I right?
Consider this monster for example:
$$\frac{1}{\pi} = \frac{2\sqrt{2}}{9801} \sum^\infty_{k=0} \frac{(4k)!(1103+26390k)}{(k!)^4 396^{4k}}$$Has someone provided a formal proof?
In the absence of proof, there's still computational verification. As I posted to edu-sig:
Without having all the higher level math it might take to actually prove some identity, and while starting to learn of those identities nonetheless, extended precision would seem especially relevant.
"Seriously, I can generate Pi with that summation series?" Lets check.
Just writing the program is a great workout.
The full post is linked below, in the For Further Reading section.
Here's a pattern David Koski was discovering, using a calculator:
from math import atan as arctan
Ø = (1 + rt2(5))/2
arctan( Ø ** 1) - arctan( Ø ** -1) == arctan(1/2)
arctan( Ø **-1) - arctan( Ø ** -3) == arctan(1/3)
arctan( Ø **-3) - arctan( Ø ** -5) == arctan(1/7)
arctan( Ø **-5) - arctan( Ø ** -7) == arctan(1/18)
arctan( Ø **-7) - arctan( Ø ** -9) == arctan(1/47)
arctan( Ø **-9) - arctan( Ø **-11) == arctan(1/123)
. . .The three dots at the end mean "and so on", implying there's some rule. Do you see it? The exponents of Ø go negative from 1, skipping two, staying odd. The ==, remember, means "is equal too". These are assertions, which gets me thinking of unittest.
Remember the conventional picture of what arctan is doing. Picture a right triangle on a unit circle and some angle at the origin, defined by the radius (movable) and the bottom edge (fixed). The angle's measure is the arctan of opposite to adjacent edges.
Given the limitations of floating point numbers, we'll switch to gmpy2, which, unlike the Decimal type, supports arctan. I've used the Decimal type to verify the Ramanujan monster (above) but it has fewer methods, so lets exit the Standard Library in favor of a 3rd party solution.
We may enlist the unittest.TestCase method assertAlmostEqual to verify the pattern so far. We're in control of how many decimal places to check. I go with 30 in the code cell below.
I'm organizing my tests to run inside of the verify function. That's slightly unusual, but why not practice? Here's a link to the documentation.
I'm also making use of two versions of arctan, one of which takes two arguments, y, x and computes arctan(y, x). Only this latter implementation is strong enough to take us out 30 places, and then over 90 places in the next exercise.
In [1]:
import gmpy2
from gmpy2 import atan2 as arctan, atan as at
from unittest import TestCase, TextTestRunner, TestSuite
# arctan(y, x) = arctan(y/x)
# at(x) = atan(x)
gmpy2.get_context().precision=300
Ø = (1 + gmpy2.sqrt(5))/2
class TestArcTan(TestCase):
def test_claim(self):
self.assertAlmostEqual(at(Ø) - arctan(1, Ø), arctan(1,2), places=30)
self.assertAlmostEqual(arctan( 1, Ø) - arctan(1, Ø** 3), arctan(1,3), places=30)
self.assertAlmostEqual(arctan( 1, Ø**3) - arctan(1, Ø** 5), arctan(1,7), places=30)
self.assertAlmostEqual(arctan( 1, Ø**5) - arctan(1, Ø** 7), arctan(1,18), places=30)
self.assertAlmostEqual(arctan( 1, Ø**7) - arctan(1, Ø** 9), arctan(1,47), places=30)
self.assertAlmostEqual(arctan( 1, Ø**9) - arctan(1, Ø**11), arctan(1,123), places=30)
def verify():
"""
how true are these expressions?
using a function to stuff an empty TestSuite, taking from
an existing class, is not the usual pattern
"""
def suite():
suite = TestSuite()
suite.addTest(TestArcTan('test_claim'))
return suite
runner = TextTestRunner()
runner.run(suite())
verify()
We've looked at the left side of the equation, but what about the right side? What is this sequence of fractions 1/2, 1/3, 1/7, 1/18... in the resulting term. What are the next identities? How might we generate this sequence by following some rule?
A useful resource in such cases is the On-Line Encyclopedia of Integer Sequences, which comes through for us when we search on just those denominators. We find this is a "bisection of the Lucas Numbers". Koski confirmed this was the sequence he was using.
We learn from this web page that this sequence may be generated as the sum of consecutive odd-position Fibonacci numbers. Here's my first attempt at a corresponding generator:
In [2]:
def Fibo(a=1, b=0):
# 1, 0, 1, 1, 2, 3, 5, 8, 13, 21...
while True:
yield a
a, b = b, a + b
def A005248():
"""
Skip half the Fibs, sum consecutive terms from
the other half:
1 + 1, 1 + 2, 2 + 5, 5 + 13...
"""
fibo_0 = Fibo()
fibo_1 = Fibo()
# advance 2 positions
next(fibo_1); next(fibo_1) #1, 0
while True:
yield next(fibo_0) + next(fibo_1) # 1 + 1, 1 + 2, 2 + 5...
next(fibo_0) # skip over next Fib
next(fibo_1) # skip over next Fib
a005248 = A005248()
print([next(a005248) for _ in range(20)])
Now we have the ability to build successive arctan expressions and verify them. Rather then depend on unittest again, lets simply compute the floating point difference between the two sides of the equation and print it out.
If it stays extremely small, we will consider the expression verified to some degree of precision.
In [3]:
def build_expr():
a005248 = A005248()
exp1 = 1
exp2 = -1
template = "(at(Ø**{}) - at(Ø**{})) - arctan(1, {})"
while True:
yield template.format(exp1, exp2, next(a005248))
exp1 -= 2
exp2 -= 2
expr_gen = build_expr()
for _ in range(20):
expr = next(expr_gen)
print("Expr: {}".format(expr))
print("Diff: {0:4.20e}".format(eval(expr)))
Looking good!
Notice the differences come back in exponential notation and end with large negative exponents.
A related sequence identified by Koski is:
arctan(Ø) == arctan(1/2) + arctan(1/3) + arctan(1/7) + arctan(1/18) . . .Again, we're using half the Lucas numbers, a so-called bisection. Might we verify?
Lets generalize with a class named Lucas built around the Fibonaccis. The goal is to make it easy to initialize however we like. Maybe we only want ever other Lucas number, beginning with some specific start value.
Start where you like, and skip as you like.
In [4]:
from functools import partial
gmpy2.get_context().precision=1000
Ø = (1 + gmpy2.sqrt(5))/2 # recomputed because more precise than ever before
class Lucas:
def __init__(self, start=0, skip=0):
lucas = partial(Fibo, 2, 1)
self.seq = lucas()
self.skip = skip
if start > 0: # only used once
for _ in range(start):
next(self.seq)
def __next__(self):
val = next(self.seq)
if self.skip > 0:
for _ in range(self.skip):
next(self.seq)
return val
def __iter__(self):
return self
s = Lucas(skip=1)
terms = [next(s) for _ in range(500)]
result = sum(map(lambda x: gmpy2.atan2(1,x), terms))
print("{0:200.198f}".format(result))
print("{0:200.198f}".format(gmpy2.atan(Ø)))
They match to a satisfying degree, wouldn't you say?
For a final act, lets verify an identity grabbed from MathWorld, that lets us generate Fibonacci and Lucus numbers with any skip distance between the first two.
We'll always say:
F(1, offset) == 1
but then
F(2, offset) == 1 + offset;
F(3, offset) == F(2, offset) + offset
F(4, offset) == F(3) + F(2) # offset is fixed
F(5, offset) == F(4) + F(3)
...and so on.
In [5]:
def F(n, offset):
num_term0 = (offset + gmpy2.sqrt(offset ** 2 + 4))**n
num_term1 = (offset - gmpy2.sqrt(offset ** 2 + 4))**n
denom = 2**n * gmpy2.sqrt(offset ** 2 + 4)
return gmpy2.round_away((num_term0 - num_term1)/denom)
genexp = (F(gmpy2.mpfr(n), 1) for n in range(1, 20))
for num in genexp:
print(int(num), end=" ")
print()
genexp = (F(gmpy2.mpfr(n), 2) for n in range(1, 20))
for num in genexp:
print(int(num), end=" ")
print()
I see we've picked up the Fibonacci and lambda numbers ( A000129).
For Further Reading: