5
\$\begingroup\$

I was solving the DIVFACT problem from Sphere Online Judge:

Given a number, find the total number of divisors of the factorial of the number.
Since the answer can be very large, print the answer modulo \$10^9+7\$.

Input

The first line contains \$T\$, the number of test cases

\$T\$ lines follows each containing the number \$N\$.

Constraints:

\$ 1 \le T \le 500\$

\$ 0 \le N \le 50000\$

To solve this I took the following approach: I first calculated the primes up to 50000 using a sieve, then for each prime \$p\$ computed the power of \$p\$ in \$N!\$ given by:

$$ \mathrm{power}(p) = \frac{n}{p} + \frac{n}{p^2} + \frac{n}{p^3} + \ldots $$

… until I got a \$p^\mathrm{power}\$ greater than \$N\$.

Finally I computed the answer for each prime \$p\$:

$$ (\mathrm{power}(p_1)+1) \cdot (\mathrm{power}(p_2)+1) \cdot \ldots \cdot (\mathrm{power}(p_k)+1)$$

assuming there are \$k\$ prime numbers less than 50000.

def prime_numbers(limit=50001):
    '''Prime number generator. Yields the series
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
    using Sieve of Eratosthenes.
    '''
    yield 2
    sub_limit = int(limit**0.5)
    flags = [True, True] + [False] * (limit - 2)
    # Step through all the odd numbers
    for i in range(3, limit, 2):
        if flags[i]:
            continue
        yield i
        # Exclude further multiples of the current prime number
        if i <= sub_limit:
            for j in range(i*i, limit, i<<1):
                flags[j] = True

p=list(prime_numbers())
# print p

for i in range(int(raw_input())):
    n=int(raw_input())
    d=[]
    pr=1
    for num in p:
        if num<=n:
            exp=1
            s=0
            while n>=num**exp:
                s+=int(n/(num**exp))
                exp+=1
            pr*=(s+1)
            # d.append(s)
        else:
            break
    print pr%1000000007

but this code and algo is not fast enough to pass the test cases in the given time limit.Is there any way to improve this code to make it faster

\$\endgroup\$
  • \$\begingroup\$ hmm, couldn't you use logarithms to find the power of each prime that is close to the test number? \$\endgroup\$ – Tadhg McDonald-Jensen Jun 2 '16 at 20:03
  • \$\begingroup\$ @TadhgMcDonald-Jensen can you elobrate I am not getting the point \$\endgroup\$ – john smith Jun 2 '16 at 20:08
  • 2
    \$\begingroup\$ nevermind the log, just tested and math.log isn't exactly accurate so it probably wouldn't work.. but the loop while n>=num**exp : s+=int(n/(num**exp)) ; exp+=1 calculates num**exp twice per iteration, it would be more efficient to do power = 1 ; while n>=power: power*=num ; s+=n//power \$\endgroup\$ – Tadhg McDonald-Jensen Jun 2 '16 at 20:34
  • 1
    \$\begingroup\$ @DarthGizka that is awesome this actually reduced time to 0.37 (and this is the fastest solution for this problem ) I also improved the prime sieve and updated this in your post along with your code \$\endgroup\$ – john smith Jun 4 '16 at 16:54
  • 1
    \$\begingroup\$ @John: I've updated my answer with a new way of speeding things up (by a factor of 20 for the SPOJ limit of 50000). You might find this interesting. \$\endgroup\$ – DarthGizka Jun 6 '16 at 17:31
5
\$\begingroup\$

Things should be fine if you perform strength reduction on n / num**exp. That not only makes things faster, it also gets rid of the messiness inherent in Python's exponentiation operator (so that you don't have to find out how it actually works under the hood).

In the kth iteration of the inner loop you have a quotient q_k that is equal to n / p**k (i.e. n / num**exp in your code).

To obtain n / p**(k+1) in the next iteration, all you have to do is q_k / p because that makes it (n / p**k) / p, which is equal to n / p**(k+1) as desired. In other words, instead of computing the power of the prime anew and then performing a division, simply divide the old quotient by the prime again.

Here's the complete loop, coded in a more pedestrian language:

static ushort[] primes_up_to_LIMIT = U16.SmallPrimesUpTo(50000).ToArray();

...

long divisor_count = 1;

foreach (int prime in primes)
{
    if (prime > n)
        break;

    int e = 0;

    for (int q = n; q >= prime; )
        e += q /= prime; 

    divisor_count = (divisor_count * (e + 1)) % 1000000007;
}

return (int)divisor_count;

It isn't very fast - it clocked 0.12 s on SPOJ - but it should show the principle well enough and it can serve as the basis of a more performant solution.

Note: detailed investigation shows that for once, the lackluster performance of my C# submission on SPOJ is not due to usual problems inherent in .Net (e.g. poor input/output performance) but rather to some strange problem with Mono. On my laptop a worst-case input file clocks 13.7 ms under .Net and 136.4 ms under Mono.

Here is an attempt at rending the loop in Python (renaming some things for clarity):

divisor_count = 1

for current_prime in primes:

    if current_prime > n:
        break

    q = n
    e = 0
    while q >= current_prime:
        q /= current_prime
        e += q

    divisor_count *= e + 1

print divisor_count % 1000000007

One thing worth looking into is the underlying data type of the divisor count (pr in the original code). This can grow quickly, far beyond anything representable in a plain integer, and hence it will go big integer internally. This makes the multiplication pr *= s + 1 more expensive, and the statement has a cost multiplier of max(T) * pi(max(N)) = 500 * 5133.

Reducing the division count modulo 1000000007 after every update by changing

divisor_count *= e + 1

to

divisor_count = (divisor_count * (e + 1)) % 1000000007

improved the time for a worst-case input file from 1.92 s to 1.39 s on my box, so it's something worth keeping in mind. For comparison: the original code clocked 2.64 s for this file.

Belated addition: a few days after the fact an idea hit me out of the blue while I was having a smoke, and I'm reporting it here because it can speed up things by a factor of 20 for the SPOJ limit of 50000 (much more for higher inputs). This idea is based on the observation that half of the primes p satisfy

(n / 2) < p <= n

and hence

1 <= (n / p) < 2

which is just a coy way of saying

n / p = 1

Knowing that there are k primes that satisfy this condition (bisect_right() to the rescue) I can compute the contribution for all these primes in one fell swoop as

(1 + 1)**k

This readily extends to other small quotients in the remaining range of primes, and the fun only stops when the primes get so small that their squares appear as divisors (i.e. p <= sqrt(n)).

Here's how you explain this to your computer in the pythonic manner:

def v6b (n):
    """number of divisors of n! (SPOJ DIVFACT)"""
    divisor_count = 1
    hi = bisect.bisect_right(primes, n)
    for m in range(2, int(math.sqrt(n))):
        lo = bisect.bisect_right(primes, n / m, 0, hi)
        k  = hi - lo
        hi = lo
        divisor_count = (divisor_count * m**k) % 1000000007
        if k < 2:
            break
    # the rest is algorithm v5b but without regular modulo trimming since growth is contained
    for current_prime in primes[0:hi]:
        q = n
        e = 0
        while q >= current_prime:
            q /= current_prime
            e += q
        divisor_count *= e + 1
    return divisor_count % 1000000007

Here are the timings for the usual 500 worst-case inputs, first for the earlier version (here called 'v5b') and then for this new version. I went beyond the SPOJ limit of 50000 to give some idea of the growth behaviour.

# timing v5b ...
500 @    500:    0.016  // cksum 230711749546
500 @   1000:    0.016  // cksum 256088767231
500 @   5000:    0.078  // cksum 248659270107
500 @  10000:    0.156  // cksum 255839576411
500 @  50000:    1.281  // cksum 252317442043
500 @ 100000:    3.313  // cksum 246710058695
500 @ 500000:   44.376  // cksum 244700225412

# timing v6b ...
500 @    500:    0.000  // cksum 230711749546
500 @   1000:    0.015  // cksum 256088767231
500 @   5000:    0.016  // cksum 248659270107
500 @  10000:    0.031  // cksum 255839576411
500 @  50000:    0.063  // cksum 252317442043
500 @ 100000:    0.078  // cksum 246710058695
500 @ 500000:    0.266  // cksum 244700225412

Beyond 100000 or thereabouts the code gets faster with a hand-coded powering function (via repeated squaring) but the builtin exponentiation operator is faster for lower limits, and hence on SPOJ. In C++ and C# things were faster when bailing out of the first loop for k < 12 but in Python it worked best with k < 2.

Note: the time for 500 test cases at the SPOJ limit of 50000 is virtually the same as for the C# version running under Mono on my box, which is no mean feat for an interpreted language...

\$\endgroup\$
  • \$\begingroup\$ what is primes_up_to_sqrt_LIMIT are you saying that in my second part of code I should make a list of primes in range(sqrt(input)) \$\endgroup\$ – john smith Jun 3 '16 at 18:52
  • \$\begingroup\$ That's a brainfart - the symbol name is lying through its teeth. I'll correct it. However, the important thing is the strength reduction and the form it gives to the inner loop. The rest is just to provide context. Sorry for the mess... \$\endgroup\$ – DarthGizka Jun 3 '16 at 18:57
  • \$\begingroup\$ can you edit your code and write it in pythonic of simple english way \$\endgroup\$ – john smith Jun 3 '16 at 19:08
  • 1
    \$\begingroup\$ @john: I tried to make the text clearer and changed exponentiation to Python style. If there is more that needs clarification, just ask and I will do my best. \$\endgroup\$ – DarthGizka Jun 3 '16 at 19:17
  • 1
    \$\begingroup\$ @john: I guess it comes down to input/output performance. I recommend that you do the INTEST and INOUTEST tasks (see links in my answer) to see how Python's builtin I/O and conversion functions fare in that regard. In C# I had cases where I got TLE with the builtin stuff, and switching to homegrown I/O and conversion cut the SPOJ time by a factor of 100 without any change in the problem-solving part of the code. Also, try benching the numeric part of your code (without I/O or conversion between strings and numbers). If the time on your computer is less than 100 ms then it is okay already. \$\endgroup\$ – DarthGizka Jun 3 '16 at 20:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.