记$S(n)$为整数$n$的所有约数所构成的集合。

若$S(n)$的某个子集A中只包含一个元素，或者A中任意一个元素均不能整除其它元素，则称A为$S(n)$的反链。

例如$S(30)=\{1, 2, 3, 5, 6, 10, 15, 30\}$，$\{2, 5, 6\}$不是$S(30)$的反链，$\{2, 3, 5\}$是$S(30)$的反链。

记$N(n)$为$S(n)$的反链的最大长度，求$\sum_{n=1}^{10^8}N(n)$。

本题难度:

解答

根据OEIS A096825中的公式，若 $$n=p_1^{e_1}\ldots p_m^{e_m},$$ 是n的质因数分解，那么最大反链长度就是多项式$\prod_{i=1}^m(1+x+\ldots+x^{e_i})$中$x^k$前的系数，其中$k=\lfloor(e_1+\ldots+e_m)/2\rfloor$。

遍历n，依次计算n的质因数分解，再递归求出系数即得结果$528755790$。

注：因需cache装饰器作记忆化搜索，以下代码为Python 3，且代码中打印了进度信息。

from functools import *

target=10000
d=[0]*target
n=2
while n < target:
    for i in range(n*n,target,n):
        d[i]=d[i]+1
    n=n+1
    while n < target and d[n]>0:
        n=n+1

primeList=[k for k in range(2,target) if d[k]==0]

def decompose(n):
    z=[]
    for p in primeList:
        if p*p>n:
            break
        k=0
        while n%p==0: 
            n//=p
            k+=1
        if k: 
            z.append(k)
    if n>1:
        z.append(1)
    return tuple(z)

@cache
def cal(exponents,n):
    if n < 0: return 0
    if n==0: return 1
    if sum(exponents) < n: return 0
    if len(exponents)==1: return 1
    return sum(cal(tuple(exponents[1:]),n-i)for i in range(exponents[0]+1))

s=1
tick=10**6
for i in range(2,10**8+1):
    exp=decompose(i)
    s+=cal(exp,sum(exp)//2)
    if i%tick==0:
        print(i//tick,"percent completed")
        
print(s)