Finding the total divisors of a number up to 10^12 as fast as possible?

Question

I need to calculate the total number of divisors of a number N (without caring about what the values of the divisors are) and do so within 40-80 operations for all such numbers N. How can I do it? This is not a homework question. I tried out Pollard's Rho algorithm but even that turned out too slow for my purposes. Here is my code in python. How can I improve its performance, if possible?

def is_prime(n):    
    if n < 2:
        return False
    ps = [2,3,5,7,11,13,17,19,23,29,31,37,41,
         43,47,53,59,61,67,71,73,79,83,89,97]
    def is_spsp(n, a):
        d, s = n-1, 0
        while d%2 == 0:
            d /= 2; s += 1
        t = pow(int(a),int(d),int(n))
        if t == 1:
            return True
        while s > 0:
            if t == n-1:
                return True
            t = (t*t) % n
            s -= 1
        return False
    if n in ps: return True
    for p in ps:
        if not is_spsp(n,p):
            return False
    return True

def gcd(a,b):
        while b: a, b = b, a%b
        return abs(a)

def rho_factors(n, limit=100):
    def gcd(a,b):
        while b: a, b = b, a%b
        return abs(a)
    def rho_factor(n, c, limit):
        f = lambda x:    (x*x+c) % n
        t, h, d = 2, 2, 1
        while d == 1:
            if limit == 0:
                raise OverflowError('limit exceeded')
            t = f(t); h = f(f(h)); d = gcd(t-h, n)
        if d == n:
            return rho_factor(n, c+1, limit)
        if is_prime(d):
            return d
        return rho_factor(d, c+1, limit)
    if -1 <= n <= 1: return [n]
    if n < -1: return [-1] + rho_factors(-n, limit)
    fs = []
    while n % 2 == 0:
        n = n // 2; fs = fs + [2]
    if n == 1: return fs
    while not is_prime(n):
        f = rho_factor(n, 1, limit)
        n = int(n / f)
        fs = fs + [f]
    return sorted(fs + [n])

def divs(n):
    if(n==1):
        return 1
    ndiv=1
    f=rho_factors(n)
    l=len(f)
    #print(f)
    c=1
    for x in range(1,l):
        #print(f[x])
        if(f[x]==f[x-1]):
            c=c+1
        else:
            ndiv=ndiv*(c+1)
            c=1
       # print ("C",c,"ndiv",ndiv)
    ndiv=ndiv*(c+1)
    return ndiv

Answer 1

First, do you mean find the total number of divisors, the number of primes in its factorization, or the number of distinct prime divisors? For instance, 12 = 2 * 2 * 3 has 6 divisors (1,2,3,4,6,12), 3 primes in its factorization (2,2,3), and 2 distinct prime divisors (2,3). Do you want 6, 3, or 2 as your result? I'm going to assume you want the second of these for the rest of this, but I don't think anything materially changes if you were interested in one of the others...

Second, you're going to have to fully factor your number. There is no known shortcut that can find the number of prime factors without finding the factors themselves. (With the notable exception that you can test whether the number of factors is ==1 or >=2 quickly.)

10^12 is not that big. You only need to test divisors up to the square root of the number, which is at most 10^6. Say a divide takes 20 cycles on a modern CPU at 2GHz, that's only 10 milliseconds to test a million divisors.

#include <stdio.h>
#include <stdlib.h>
int main(int argc, char *argv[]) {
  long long n = atoll(argv[1]);
  for (int i = 2; i < 1000000; i++) {
    while (n % i == 0) { printf("%d\n", i); n /= i; }
  }
  if (n > 1) printf("%lld\n", n);
}

Takes 23 milliseconds on my machine. Wonder where that other 13 milliseconds went?

Python is about 10x slower, as this code still takes only 0.23 seconds on my machine:

import sys
n = int(sys.argv[1])
for i in xrange(2, 1000000):
  while n%i==0: print i; n/=i
if n>1: print n

How fast do you want it?

Answer 2

I remember there is a solution based on sum of digits in a number and other features
As an example, 3411 is divisible by 9, because 3+4+1+1 = 9, sum of digits is divisible by 9, than a number is also divisible by 9. with other numbers rules are similar.

Answer 3

I remember solving this problem before on SPOJ, but I don't remember the exact method I used (it'd be great if you provide the problem ID). Did you try the naive method here? It has a complexity of O(sqrt n), which is about O(10 ^ 6) in the worst case. The modulo operator might be a bit slow, but it's worth giving it a try. Here's how it should look when done in C++:

int cntdiv = 0;
for(int i = 2; i * i <= x; i ++) {
    if(x % i == 0) {
        cntdiv += 1 + (i * i != x);
    }
}
//cntdiv is now your count