See http://math.stackexchange.com/a/926138/2582

tm.hs

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

tm.py

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])))