Any faster algorithm to compute the number of divisors

Question

The F series is defined as

F(0) = 1

F(1) = 1

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

The task is to find the number of different divisors for F(i)

This question is from Timus . I tried the following Python but it surely gives a time limit exceeded. This bruteforce approach will not work for a large input since it will cause integer overflow as well.

#!/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

Any optimization?

EDIT I have tried the suggestion, and rewrite the solution: (not storing the F(n) value directly, but a list of factors)

#!/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)

And it still gives TL5, not faster enough, but this solves the problem of overflow for value F(n).

Answer 1

First see this wikipedia article on the divisor function. In short, if you have a number and you know its prime factors, you can easily calculate the number of divisors (get SO to do TeX math):

$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}$

Anyway, it's a simple function.

Now, to solve your problem, instead of keeping F(n) as the number itself, keep it as a set of prime factors and exponent sizes. Then the function that calculates F(n) simply takes the two sets for F(n-1) and F(n-2), sums the exponents of the same prime factors in both sets (assuming zero for nonexistent ones) and additionally adds the set of prime factors and exponent sizes for the number i. This means that you need another simple¹ function to find the prime factors of i.

Computing F(n) this way, you just need to apply the above formula (taken from Wikipedia) to the set and there's your value. Note also that F(n) can quickly get very large. This solution also avoids usage of big-num libraries (since no prime factor nor its exponent is likely to go beyond 4 billion²).

---

¹ Of course this is not so simple for arbitrarily large i, otherwise we wouldn't have any form of security right now, but for your application it should be simple enough.

² Well it might. If you happen to figure out a simple formula answering your question given any n, then large ns would also be possible in the test case, for which this algorithm is likely going to give a time limit exceeded.

Answer 2

That is a fun problem.

The F(n) grow extremely fast. Since F(n) <= F(n+1) for all n, we have

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

for all n, and thus

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

for n > 2. That crude estimate already shows that one cannot store these numbers for any but the smallest n. By that F(100) requires more than (2^49) bits of storage, and 128 GB are only 2^40 bits. Actually, the prime factorisation of F(100) is

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

and that would require about 9.6 * 10^20 (roughly 2^70) bits - a little less than half of them are trailing zeros, but even storing the numbers à la floating point numbers with a significand and an exponent doesn't bring the required storage down far enough.

So instead of storing the numbers themselves, one can consider the prime factorisation. That also allows an easier computation of the number of divisors, since

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

Now, let us investigate the prime factorisations of the F(n) a little. We begin with the

Lemma: A prime p divides F(n) if and only if p <= n.

That is easily proved by induction: F(0) = F(1) = 1 is not divisible by any prime, and there are no primes <= 1.

Now suppose that n > 1 and

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

holds for k < n. Then, since

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

the set prime factors of F(n) is the union of the sets of prime factors of n, F(n-1) and F(n-2).

By the induction hypothesis, the set of prime factors of F(k) is

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

for k < n. Now, if n is composite, all prime factors of n are samller than n, hence the set of prime factors of F(n) is P(n-1), but since n is not prime, P(n) = P(n-1). If, on the other hand, n is prime, the set of prime factors of F(n) is

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

With that, let us see how much work it is to track the prime factorisation of F(n) at once, and update the list/dictionary for each n (I ignore the problem of finding the factorisation of n, that doesn't take long for the small n involved).

The entry for the prime p appears first for n = p, and is then updated for each further n, altogether it is created/updated N - p + 1 times for F(N). Thus there are

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

updates in total. For N = 10^6, about 3.6 * 10^10 updates, that is way more than can be done in the allowed time (0.5 seconds).

So we need a different approach. Let us look at one prime p alone, and follow the exponent of p in the F(n).

Let v_p(k) be the exponent of p in the prime factorisation of k. Then we have

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

and we know that v_p(F(k)) = 0 for k < p. So (assuming p is not too small to understand what goes on):

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

So we get Fibonacci numbers for the exponents, v_p(F(p+k)) = Fib(k+1) - for a while, since later multiples of p inject further powers of 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)

but the additional powers from 2*p also follow a nice Fibonacci pattern, and we have v_p(F(2*p+k)) = Fib(p+k+1) + Fib(k+1) for 0 <= k < p.

For further multiples of p, we get another Fibonacci summand in the exponent, so

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

-- until n >= p², because multiples of p² contribute two to the exponent, and the corresponding summand would have to be multiplied by 2; for multiples of p³, by 3 etc.

One can also split the contributions of multiples of higher powers of p, so one would get one Fibonacci summand due to it being a multiple of p, one for it being a multiple of p², one for being a multiple of p³ etc, that yields

              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

Now, in particular for the smaller primes, these sums have a lot of terms, and computing them that way would be slow. Fortunately, there is a closed formula for sums of Fibonacci numbers whose indices are an arithmetic progression, for 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

where

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

and Luc(k) is the k-th Lucas number, Luc(k) = Fib(k+1) + Fib(k-1).

For our purposes, we only need the Fibonacci numbers modulo 10^9 + 7, then the division must be replaced by a multiplication with the modular inverse of D(s).

Using these facts, the number of divisors of F(n) modulo 10^9+7 can be computed in the allowed time for n <= 10^6 (about 0.06 seconds on my old 32-bit box), although with Python, on the testing machines, further optimisations might be necessary.