任何更快的算法来计算除数的数量

Question

F系列定义为

F(0) = 1

F(1) = 1

F(i) = i * F(i - 1) * F(i - 2)  for i > 1

任务是找出 F(i) 的不同除数的个数

这个问题来自Timus。我尝试了以下 Python，但它肯定会超出时间限制。这种暴力破解方法不适用于大输入，因为它也会导致整数溢出。

#!/usr/bin/env python

from math import sqrt
n = int(raw_input())

def f(n):
    global arr
    if n == 0:
        return 1
    if n == 1:
        return 1
    a = 1
    b = 1
    for i in xrange(2, n + 1):
        k = i * a * b
        a = b
        b = k
    return b

x = f(n)
cnt = 0
for i in xrange(1, int(sqrt(x)) + 1):
    if x % i == 0:
        if x / i == i:
            cnt += 1
        else:
            cnt += 2
print cnt

有什么优化吗？

编辑 我已经尝试了这个建议，并重写了解决方案：（不是直接存储 F(n) 值，而是一个因素列表）

#!/usr/bin/env python

#from math import sqrt

T = 10000
primes = range(T)

primes[0] = False
primes[1] = False
primes[2] = True
primes[3] = True

for i in xrange(T):
    if primes[i]:
        j = i + i
        while j < T:
            primes[j] = False
            j += i

p = []
for i in xrange(T):
    if primes[i]:
        p.append(i)

n = int(raw_input())

def f(n):
    global p
    if n == 1:
        return 1
    a = dict()
    b = dict()
    for i in xrange(2, n + 1):
        c = a.copy()
        for y in b.iterkeys():
            if c.has_key(y):
                c[y] += b[y]
            else:
                c[y] = b[y]
        k = i
        for y in p:
            d = 0
            if k % y == 0:
                while k % y == 0:
                    k /= y
                    d += 1
                if c.has_key(y):
                    c[y] += d
                else:
                    c[y] = d
                if k < y: break
        a = b
        b = c
    k = 1
    for i in b.iterkeys():
        k = k * (b[i] + 1) % (1000000007)
    return k

print f(n)

它仍然给出了 TL5，速度不够快，但这解决了值 F(n) 的溢出问题。

Answer 1

首先看到这个wikipedia article on the divisor function。简而言之，如果你有一个数字并且你知道它的质因数，你就可以很容易地计算出除数的数量（让 SO 做 TeX 数学）：

$n = \prod_{i=1}^r p_i^{a_i}$

$\sigma_x(n) = \prod_{i=1}^{r} \frac{p_{i}^{(a_{i}+1)x}-1}{p_{i}^x-1}$

反正就是一个简单的函数。

现在，为了解决您的问题，不要将 F(n) 保留为数字本身，而是将其保留为一组素因子和指数大小。然后计算F(n) 的函数只需将F(n-1) 和F(n-2) 的两个集合，将两个集合中相同素因子的指数相加（假设不存在的素因子为零），并另外添加素因子和指数的集合数字i 的大小。这意味着您需要另一个简单的¹ 函数来找到i 的素数。

以这种方式计算F(n)，您只需将上述公式（取自维基百科）应用于集合，这就是您的价值。另请注意，F(n) 可以很快变得非常大。此解决方案还避免使用大数库（因为没有素因数或其指数可能超过 40 亿²）。

---

¹ 当然这对于任意大的i 来说并不是那么简单，否则我们现在不会有任何形式的安全性，但是对于您的应用程序来说它应该足够简单。

² 可能。如果你碰巧找到了一个简单的公式来回答你的问题给定任何n，那么在测试用例中也可以使用大的ns，这个算法很可能会超出时间限制。

Answer 2

这是一个有趣的问题。

F(n) 增长极快。由于F(n) <= F(n+1) 对所有n，我们有

F(n+2) > F(n)²

对于所有n，因此

F(n) > 2^(2^(n/2-1))

为n > 2。这个粗略的估计已经表明，除了最小的n 之外，不能存储这些数字。由此F(100) 需要超过(2^49) 位的存储空间，而128 GB 仅是2^40 位。实际上，F(100) 的素因数分解是

*Fiborial> fiborials !! 100
[(2,464855623252387472061),(3,184754360086075580988),(5,56806012190322167100)
,(7,20444417903078359662),(11,2894612619136622614),(13,1102203323977318975)
,(17,160545601976374531),(19,61312348893415199),(23,8944533909832252),(29,498454445374078)
,(31,190392553955142),(37,10610210054141),(41,1548008760101),(43,591286730489)
,(47,86267571285),(53,4807526976),(59,267914296),(61,102334155),(67,5702887),(71,832040)
,(73,317811),(79,17711),(83,2584),(89,144),(97,3)]

这将需要大约9.6 * 10^20（大约2^70）位 - 其中不到一半是尾随零，但即使将数字存储为带有有效数字和指数的浮点数也不会带来所需的存储空间足够低。

因此，可以考虑素因数分解，而不是存储数字本身。这也允许更容易计算除数的数量，因为

              k                    k
divisors(n) = ∏ (e_i + 1)  if  n = ∏ p_i^e_i
             i=1                  i=1

现在，让我们稍微研究一下F(n) 的质因式分解。我们从

引理：素数p除以F(n)当且仅当p <= n。

这很容易通过归纳证明：F(0) = F(1) = 1 不能被任何素数整除，也没有素数<= 1。

现在假设n > 1 和

A(k) = The prime factors of F(k) are exactly the primes <= k

适用于k < n。那么，由于

F(n) = n * F(n-1) * F(n-2)

F(n) 的素因数集是n、F(n-1) 和F(n-2) 的素因数集的并集。

根据归纳假设，F(k)的素因数集是

P(k) = { p | 1 < p <= k, p prime }

为k < n。现在，如果n 是合数，则n 的所有质因数都比n 小，因此F(n) 的质因数集是P(n-1)，但由于n 不是质数，P(n) = P(n-1) .另一方面，如果n 是质数，则F(n) 的质因数集是

P(n-1) ∪ {n} = P(n)

这样，让我们​​看看一次跟踪F(n) 的素因数分解并更新每个n 的列表/字典需要做多少工作（我忽略了找到n 的因数分解的问题，这对于涉及的小n 来说并不需要很长时间）。

素数p 的条目首先出现在n = p 中，然后为每一个n 更新，总共为F(N) 创建/更新N - p + 1 次。因此有

   ∑   (N + 1 - p) = π(N)*(N+1) -  ∑ p ≈ N²/(2*log N)
p <= N                          p <= N

总共更新。对于N = 10^6，大约3.6 * 10^10 更新，这远远超出了允许的时间（0.5 秒）内所能完成的。

所以我们需要一种不同的方法。让我们单独看一个素数p，然后在F(n) 中追随p 的指数。

让v_p(k) 成为p 在k 的素数分解中的指数。然后我们有

v_p(F(n)) = v_p(n) + v_p(F(n-1)) + v_p(F(n-2))

我们知道v_p(F(k)) = 0 对应于k < p。所以（假设p 不是太小而无法理解发生了什么）：

v_p(F(n))     = v_p(n) +   v_p(F(n-1))  +   v_p(F(n-2))
v_p(F(p))     =   1    +       0        +       0        =  1
v_p(F(p+1))   =   0    +       1        +       0        =  1
v_p(F(p+2))   =   0    +       1        +       1        =  2
v_p(F(p+3))   =   0    +       2        +       1        =  3
v_p(F(p+4))   =   0    +       3        +       2        =  5
v_p(F(p+5))   =   0    +       5        +       3        =  8

所以我们得到了指数的斐波那契数，v_p(F(p+k)) = Fib(k+1) - 有一段时间，因为后来的倍数 p 注入了更多的力量 p，

v_p(F(2*p-1)) =   0    +     Fib(p-1)   +     Fib(p-2)   =     Fib(p)
v_p(F(2*p))   =   1    +     Fib(p)     +     Fib(p-1)   = 1 + Fib(p+1)
v_p(F(2*p+1)) =   0    + (1 + Fib(p+1)) +     Fib(p)     = 1 + Fib(p+2)
v_p(F(2*p+2)) =   0    + (1 + Fib(p+2)) + (1 + Fib(p+1)) = 2 + Fib(p+3)
v_p(F(2*p+3)) =   0    + (2 + Fib(p+3)) + (1 + Fib(p+2)) = 3 + Fib(p+4)

但来自2*p 的额外权力也遵循漂亮的斐波那契模式，我们有v_p(F(2*p+k)) = Fib(p+k+1) + Fib(k+1) 对应0 <= k < p。

对于p 的进一步倍数，我们在指数中得到另一个斐波那契和数，所以

            n/p
v_p(F(n)) =  ∑ Fib(n + 1 - k*p)
            k=1

-- 直到n >= p²，因为p² 的倍数对指数贡献了2，相应的和必须乘以2； p³ 的倍数，3 等。

还可以拆分p 的高次幂的倍数的贡献，因此由于它是p 的倍数，因此可以得到一个斐波那契和数，一个是p² 的倍数，一个是是p³ 等的倍数，产生

              n/p                    n/p²                    n/p³
v_p(F(n)) =    ∑ Fib(n + 1 - k*p)  +  ∑ Fib(n + 1 - k*p²)  +  ∑ Fib(n + 1 - k*p³)  + ...
              k=1                    k=1                     k=1

现在，特别是对于较小的素数，这些和有很多项，以这种方式计算它们会很慢。幸运的是，对于指数为算术级数的斐波那契数之和有一个封闭公式，0 < a <= s

 m
 ∑ Fib(a + k*s) = (Fib(a + (m+1)*s) - (-1)^s * Fib(a + m*s) - (-1)^a * Fib(s - a) - Fib(a)) / D(s)
k=0

在哪里

D(s) = Luc(s) - 1 - (-1)^s

而Luc(k) 是k-th 卢卡斯号码，Luc(k) = Fib(k+1) + Fib(k-1)。

就我们的目的而言，我们只需要以10^9 + 7 为模的斐波那契数，然后除法必须用与D(s) 的模逆的乘法来代替。

使用这些事实，可以在 n <= 10^6 的允许时间内计算 F(n) 模 10^9+7 的除数（在我的旧 32 位机器上大约 0.06 秒），尽管使用 Python，在测试中机器，可能需要进一步优化。