Skip to content

Instantly share code, notes, and snippets.

@rubik
Last active August 29, 2015 14:06
Show Gist options
  • Save rubik/bbb903734c6221f81ed0 to your computer and use it in GitHub Desktop.
Save rubik/bbb903734c6221f81ed0 to your computer and use it in GitHub Desktop.
module Main
where
import System.Environment (getArgs)
import Data.List (transpose, sort, groupBy)
import Math.Sieve.Factor (factor, sieve)
-- The Thue-Morse sequence
thueMorse :: [Int]
thueMorse = 0 : interleave (map (1-) thueMorse) (tail thueMorse)
where interleave (x:xs) ys = x : interleave ys xs
-- Fractions in the product P. Each fraction is represented as a tuple.
cc :: [[Integer]]
cc = zipWith (
b -> [[2*n+1, 2*n+2], [2*n+2,2*n+1]]!!b) [0..] thueMorse
-- Sort, group with regard to the first element, and map with f
mapGroup :: (Ord a, Ord b) => ([(a, b)] -> c) -> [(a, b)] -> [c]
mapGroup f = map f . groupBy eq . sort
where eq a b = fst a == fst b
-- Split a list of tuples in two lists, according to the second element of
-- each tuple.
takeOut :: (Ord a, Num a) => [[a]] -> [[[a]]]
takeOut xs = takeOut' xs [] []
where takeOut' [] fs qs = [fs,qs]
takeOut' ([b,e]:xs) fs qs
| e > 0 = takeOut' xs ([b,e]:fs) qs
| e < 0 = takeOut' xs fs ([b,-e]:qs)
| otherwise = takeOut' xs fs qs
-- Subtract two list. The result is a list.
-- > sub x y
-- > [{x - y}, {y - x}]
sub :: (Ord a, Num a) => [(a, a)] -> [(a, a)] -> [[[a]]]
sub x y = takeOut . mapGroup collapse $ x ++ y'
where y' = [(b, -e) | (b,e) <- y]
collapse l = [fst $ head l, sum [e | (b,e) <- l]]
-- Compute the product of [(b,e)] after computing b^e.
prod :: (Num c, Integral b) => [[b]] -> c
prod l = fromIntegral . product $ [b^e | [b,e] <- l]
-- Factor the given number with the given sieve limit.
factor' :: (Integral a) => Int -> a -> [(a, a)]
factor' l = factor (sieve l)
-- Factor a list of numbers and join all the results into a single list.
factorize :: (Integral a) => Int -> [a] -> [(a, a)]
factorize l ns = mapGroup getsum $ concatMap (factor' l) ns
where getsum l = (fst $ head l, sum . map snd $ l)
-- Approximate the infinite product P by computing n levels, i.e. using all
-- the numbers up to 2^n.
approx :: Int -> Double
approx n = fromRational $ prod fp' / prod fq'
where limit = 2^(n-1)
[1:ps,qs] = transpose $ take limit cc
[fp,fq] = map (factorize (2*limit)) [ps,qs]
[fp',fq'] = sub fp fq
main :: IO ()
main = do
args <- getArgs
print . approx . read . head $ args
import operator
import functools
import itertools
import collections
def thue_morse():
'''Generate the Thue-Morse sequence.
Example:
>>> import itertools
>>> list(itertools.islice(thue_morse(), 16))
[0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0]
'''
seqs = [[0]]
yield 0
while True:
new = [int(not b) for seq in seqs for b in seq]
yield from new
seqs.append(new)
def cc():
'''Yield fractions in the product P. A fraction is represented as a tuple.
Example:
>>> import itertools
>>> list(itertools.islice(cc(), 8))
[(1, 2), (4, 3), (6, 5), (7, 8), (10, 9), (11, 12), (13, 14), (16, 15)]
'''
for n, v in enumerate(thue_morse()):
if v:
yield (2*n + 2, 2*n + 1)
else:
yield (2*n + 1, 2*n + 2)
def factors(k):
'''Factor the number n.
Example:
>>> factors(12)
(2, 2, 3)
>>> factors(13)
(13,)
'''
def _gen(n):
if n <= 1:
return
prime = next((x for x in range(2, int(n ** 0.5) + 1) if not n % x), n)
yield prime
yield from _gen(n // prime)
return tuple(_gen(k))
def approx(n):
'''Find an approximation of the infinite fraction by computing n levels.
That means using all the positive integers up to 2^n.
Empirical tests suggest that lim_{n -> +oo} approx(n) = 1 / sqrt(2).
This function avoids the direct approach of multiplying all the fractions
and dividing p / q at the end. Instead, it does the following:
1. It computes all the fractions from cc().
2. It computes the factors of all the numerators and those of all the
denominators.
3. It simplifies all the factors that appear both in the numerators and
in the denominators.
4. It multiplies the remaining numerators and the denominator, dividing
the results.
All of this is done in order to get smaller numbers before dividing.
Python's ints can become arbitrarily large, but dividing two numbers with
more than 100000 takes quite some time.
'''
prod = lambda c: functools.reduce(operator.mul,
(b ** e for b, e in c.items()))
fracs = zip(*itertools.islice(cc(), 2**(n - 1)))
fp, fq = map(lambda s: collections.Counter(
i for r in map(factors, s) for i in r), fracs)
return prod(fp - fq) / prod(fq - fp)
def approx2(n):
# Naive and direct approach, it's ~9x slower than approx(n).
p = q = 1
l = 2 ** n
for i, v in enumerate(thue_morse(), 1):
if i > l:
return p / q
if v:
q *= i
else:
p *= i
def convergents():
'''Generate successive convergents to the infinite fraction.'''
# Currently this function is not used anywhere, but can be useful
p = q = 1
for n, v in enumerate(thue_morse()):
if v:
p *= 2*n + 2
q *= 2*n + 1
else:
p *= 2*n + 1
q *= 2*n + 2
yield p, q
if __name__ == '__main__':
import sys
print(approx(int(sys.argv[1])))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment