尤拉計劃386題(反鏈的最大長度) 論壇python解法學習筆記

這是瑞士的Lucy_Hedgehog 在pe論壇386題第一個帖子中給出的，也是目前我見過的當中最快的。 原作者只給出了主要片斷，並說利用了非標準庫，我東拼西湊，終於讓它執行起來。沒有達到作者說的速度，也許我找的模組不夠高效。

# This is just a code snippet.
# It needs non-standard python libraries.
# 4/15/2013: Added a compact version of pi_table.
# This is using Python 3.x

def pi_table(n):
    '''Returns a table containing pi(m) for all integer m
       in V = {n//1, n//2, n//3, n//4, ... , 1}'''
    r = int(n**0.5)
    while r*r  <= n: r += 1
    V = [n//i for i in range(1, n//r+1)] + list(range(r-1,0,-1))
    # C[i] == number or primes or products of primes > q up to i
    C = {i:i-1 for i in V}
    for q in range(2,r):
        if C[q] > C[q-1]:  # q is prime
            pi = C[q-1]  # number of primes smaller than q
            q2 = q*q
            for v in V:
                if v < q2: break
                C[v] -= C[v//q] - pi
    return C

def factorization_count(n):
    '''Returns a table where the indices ares
       tuples of exponents and the values are the
       number of times these exponents occur in
       the prime factorization of the integers 1 .. n'''
    def rec(idx, exp, r):
        w = pi[r] - idx
        if w <= 0: return
        E = tuple(sorted(exp + (1,)))
        T[E] += w
        for i in range(idx, len(P)):
            p = P[i]
            if p*p > r: break
            rec(i+1, E, r // p)
            q = p*p
            k = 2
            while q <= r:
                F = tuple(sorted(exp + (k,)))
                T[F] += 1
                if q * p <= r:
                    rec(i+1, F, r // q)
                k += 1
                q *= p            
    pi = pi_table(n)
    m = max(100, int(n**0.5)+5)
    P = all_primes(m)
    T = defaultdict(int)
    T[()] = 1
    rec(0,(),n)
    return T

def P386(n):
    x = Polynomial([0,1])
    res = 0
    T = factorization_count(n)
    for E in T:
      t = Polynomial([1])
      for e in E:
        t *= Polynomial([1]*(e+1))
      w = max(t.coef)  #lt modifed from t.coeffs
      res += w * T[E]
    return res

def all_primes(n):  ##copy from Stack Overflow 
    numbers = set(range(n, 1, -1))
    primes = []
    while numbers:
        p = numbers.pop()
        primes.append(p)
        numbers.difference_update(set(range(p*2, n+1, p)))
    return primes

def timeo1(x):
    start = time.clock( )
    print (P386(x))
    stend = time.clock( )
    thetime = stend-start
    return thetime


--執行時間
>>> print(timeo1(1000000000))
5844921982.0
9.915745406882706 （原帖說是6秒）

>>> print(timeo1(100000000))
528755790.0
2.2057292688625694

如果沒有安裝numpy，下載並安裝pip install "numpy-1.12.0b1-cp35-none-win32.whl" 需要的import語句有

import numpy as numpy
from numpy.polynomial import Polynomial as Polynomial
from collections import defaultdict
import time