Prime # counting algos, each better or more pythonic then previous

Moebius.py

#!/usr/bin/env python3
# coding: utf-8
# various algos counting prime numbers < N
# each better or more pythonic then previous

import itertools
import numpy as np
import datetime
from operator import mul
from functools import reduce

GAMMA = 0.57721566490153286061
li2 = 1.045163780117492784844588889194  # Li(x) = li(x) - li(2)


def factors(n):
    # with list
    return [p for p in ambi_sieve(n + 1) if n % p == 0]


def factors2(n):
    # with generator
    return (p for p in ambi_sieve(n + 1) if n % p == 0)  # generator


def factors3(n):
    # with yield
    for p in ambi_sieve(n + 1):
        if n % p == 0:
            yield p  # yield


def ambi_sieve(n):
    # DEAD link http://tommih.blogspot.com/2009/04/fast-prime-number-generator.html
    s = np.arange(3, n, 2)
    for m in range(3, int(n ** 0.5) + 1, 2):
        if s[(m - 3) // 2]:
            s[(m * m - 3) // 2::m] = 0
    return np.r_[2, s[s > 0]]


def all_factors(prime_dict):
    series = [[p**e for e in range(maxe + 1)]
              for p, maxe in prime_dict.items()]
    for multipliers in itertools.product(*series):
        yield reduce(mul, multipliers)


def divisors(factors):
    # Generates all divisors, unordered, from the prime factorization.
    ps = sorted(set(factors))
    omega = len(ps)

    def rec_gen(n=0):
        if n == omega:
            yield 1
        else:
            pows = [1]
            for j in range(factors.count(ps[n])):
                pows += [pows[-1] * ps[n]]
            for q in rec_gen(n + 1):
                for p in pows:
                    yield p * q

    for p in rec_gen():
        yield p


def m1(n, p):
    return -1 if (n / p) % p else 0


def m2(n, f):
    # zero if duplicated factor, -1 otherwise
    return 0 if (n / f) % f == 0 else - 1


def m(n, f):
    # lambda from above in explicit way
    if (n / f) % f == 0:
        return 0
    else:
        return -1


def mu(n):
    # Works. But difficult to understand code for n == 1
    return reduce(mul, [m1(n, p) for p in factors3(n)], 1)


def mobius(n):
    # My definition of Mobius.
    if n == 1:
        return 1
    else:
        # zero if duplicated factor, -1 otherwise
        # multiply all. if either is duplicated = 0
        return reduce(mul, [m2(n, k) for k in factors(n)], 1)


def li(x):
    # programmingpraxis.com/2011/07/29/approximating-pi
    # mathworld.wolfram.com/LogarithmicIntegral.html
    return GAMMA + np.log(np.log(x)) + sum(pow(np.log(x), k) /
                                           (reduce(mul, range(1, k + 1)) * k) for k in range(1, 101))


def Li(x):
    return li(x) - li2


def R(x):
    # mathworld.wolfram.com/RiemannPrimeCountingFunction.html
    # suboptimal. don't bother with k where mobius(k)==0
    return sum(mobius(k) / k * li(pow(x, 1.0 / k)) for k in range(1, 101))


def Ri(x):
    # mathworld.wolfram.com/RiemannPrimeCountingFunction.html
    # with generator computing Mobius(k) twice but ignoring k where Mobius == 0
    return sum(mobius(k) / k * li(pow(x, 1.0 / k)) for k in range(1, 101) if mobius(k) != 0)


def Ri2(x):
    # mathworld.wolfram.com/RiemannPrimeCountingFunction.html
    # with tuple ignoring k where Mobius == 0
    mob = tuple(mobius(k) for k in range(1, 101))
    return sum(mob[k - 1] / k * li(pow(x, 1.0 / k)) for k in range(1, 101) if mob[k - 1] != 0)


def Ri3(x):
    # mathworld.wolfram.com/RiemannPrimeCountingFunction.html
    # with generator ignoring k where Mobius == 0
    mob = ((k, mobius(k)) for k in range(1, 101))
    return sum(mu / k * li(pow(x, 1.0 / k)) for (k, mu) in mob if mu != 0)


# tuple of Mobius(k) computed outside R(x) calls
moby = tuple(mobius(k) for k in range(1, 101))


def Ri4(x):
    # mathworld.wolfram.com/RiemannPrimeCountingFunction.html
    # with tuple ignoring k where Mobius == 0
    return sum(moby[k - 1] / k * li(pow(x, 1.0 / k)) for k in range(1, 101) if moby[k - 1] != 0)


N = 1e100
# print( ambi_sieve(N), factors(N), divisors(factors(N)), mobius3(N))
tc = datetime.datetime.now()
print("ambi sieve", ambi_sieve(1e6), datetime.datetime.now() - tc)
# print divisors([2, 3, 5])

# print sorted(all_factors({2:3, 3:2, 5:1}))

# print [m(N, k) for k in factors(N)], mobius3(N)

# print [mobius3(n) for n in range(100)]

# print [mu(n) for n in range(N)]

# mob = [mobius3(k) for k in range(1, N)]
# print mob, len(mob)

tc = datetime.datetime.now()
print("Ri4", Ri4(N), datetime.datetime.now() - tc)
tc = datetime.datetime.now()
print("Ri2", Ri2(N), datetime.datetime.now() - tc)
tc = datetime.datetime.now()
print("Ri", Ri(N), datetime.datetime.now() - tc)
tc = datetime.datetime.now()
print("R", R(N), datetime.datetime.now() - tc)

# print reduce(mul, [m(n, k) for k in factors(n)], 1)

# print [mobius2(i) for i in range(1000)]

# print [factors(i) for i in range(1000)]

# for N in range(10, 10000, 1000):
#   print R(N), li(N), Li(N), N/np.log(N)

# print factors(N), factors2(N), factors3(N)