Problem at hand is finding the total number of divisors of a given number efficiently. The divisor of a number is a number less than or equal to the given number and is perfectly divisible by it. For example : divisors of 4 -> 1,2,4; divisors of 12 -> 1,2,3,4,6,12; As we are trying to find the number of divisors of a given number:
Examples :
Input : 12 Output : 6 Input : 15 Output : 4
Now we already know the conventional method we use to find the number of divisors of a given number. If you are not familiar with it, I would strongly recommend to please go through that first.Here is the link The link mentioned above provides the list of all the divisors but it can be easily modified to find the total number of divisors. The complexity of this conventional algorithm in O (sqrt(n)) where n is the number of which we are trying to find the number of divisors. This article will provide a better alternative to that by the help of an algorithm which will find the number of divisors in time complexity of O(n^(1/3)). The algorithm will provide the result up to 10^18.
Concept :
We are going to use a mathematical concept here which states that if N = a^x * b^y * c^z than total number of divisors are (x+1) * (y+1) * (z+1). If you want to read more about it, you can look here So, Now we will loop over all the primes in the range of [2, N^(1/3)] and try to reduce N to its prime factorization. We will split our number N to X*Y where X contains only prime factors in the range [2, N^(1/3)] and Y will have all the prime factors which are greater than N^(1/3). Here, we will individually calculate the no of divisors of X and Y and as we multiply them, it will provide us with the answer for N. Proof for above statement – We can see that gcd (X, Y) will be 1 because we are dividing them in the form that X and Y do not have any common factor.Let the count of divisors of N be denoted by F(N), then it is easy to prove that this function will be multiplicative in nature i.e. F(X*Y) = F(X)*F(Y). Continuing with the concept, we will calculate F(X) with the help of naive approach that is by the trial division to prime factorise X and then counting the divisors. What looks tricky is calculate the value of F(Y) but if we really think about it, that’s not really tough to calculate, we just need to check some conditions about what the value of Y can be. To be specific, there are only 4 cases: Y is a prime number, then F(Y) = 2; Y is a product of 2 prime numbers, then F(Y) = 4; Y is square of a prime number, then F(Y) = 3; Y is 1, then F(Y) = 1; Y can not be a product of 3 or more than 3 prime numbers. We will prove this by the help of contradiction. Say, Y = P * Q * R. We already know that all the factors of Y will be greater than N^(1/3), so the product of P, Q and R will become greater than N but we also know that Y will surely be less than N. These 2 statements contradict each other. Similarly, it can be done for 4 or more factors. It proves that Y cannot have more than 2 factors.
Pseudo Code :
primes = array containing primes till 10^6
N = input()
ans = 1
for all i in primes :
if i*i*i > N:
break
count = 1
while N divisible by i:
N = N/i
count = count + 1;
ans = ans * count
if N is prime:
ans = ans * 2
else if N is square of a prime:
ans = ans * 3
else if N != 1:
ans = ans * 4
Pre-Requisites:
To implement this above psuedo code, there first need to calculate all the primes upto the range of 10^6. This can be done by the help of sieve of eratosthenes. Other than that, we will also need to use Miller-Rabin primality check algorithm which is a probabilistic algorithm for checking whether a number is prime or not. If you like, you can also use Fermat Primality Test, although the former is preferred normally. I have implemented the below algorithm using Miller-Rabin test, if you would like to know more about any of them, you can look them up here for reference.
Implementation in c++ :
#include "iostream"
#include "cmath"
#include "cstdlib"
using std::cout;
using std::endl;
const int N = 1e6+5;
int isPrime[N+1];
long long int prime[N+1];
int size = 0;
long long mulmod(long long a,long long b,long long c);
long long modulo(long long a,long long b,long long c){
long long x=1,y=a; // long long is taken to avoid overflow of intermediate results
while(b > 0){
if(b%2 == 1){
x=mulmod(x,y,c);
}
y = mulmod(y,y,c); // squaring the base
b /= 2;
}
return x%c;
}
long long mulmod(long long a,long long b,long long c){
long long x = 0,y=a%c;
while(b > 0){
if(b%2 == 1){
x = (x+y)%c;
}
y = (y*2)%c;
b /= 2;
}
return x%c;
}
bool Miller_rabin(long long p,int iteration){
if(p<2){
return false;
}
if(p!=2 && p%2==0){
return false;
}
long long s=p-1;
while(s%2==0){
s/=2;
}
for(int i=0;i<iteration;i++){
long long a=rand()%(p-1)+1,temp=s;
long long mod=modulo(a,temp,p);
while(temp!=p-1 && mod!=1 && mod!=p-1){
mod=mulmod(mod,mod,p);
temp *= 2;
}
if(mod!=p-1 && temp%2==0){
return false;
}
}
return true;
}
void sieve_of_eratosthenes(){
isPrime[0]=isPrime[1]=1;
for (int i=2;i*i<=N;i++){
if (isPrime[i])
continue;
for (int j=i*i;j<=N;j+=i){
isPrime[j]=1;
}
}
for (int i=2;i<=N;i++){
if (isPrime[i]==0){
prime[size++] = i;
}
}
}
int divisors(long long int num){
int k = 20; // Number of iterations of Miller-Rabin primality test;
int ans=1; // as 1 is a divisor for all the integers.
int count=0;
int i=0;
while(i < size && prime[i]*prime[i]*prime[i]<=num)
{
count = 1;
while((num % prime[i]) == 0){
num /= prime[i];
count++;
}
ans *= count;
i++;
}
// Starting the condition to calculate F(Y).
// We can see that value of num is equal to Y only.
// check if num is prime number
if (Miller_rabin(num,k)){
ans *= 2;
}
else if (num == 1){
ans *= 1;
}
else if ((long long)sqrt(num)*(long long)sqrt(num) == num && Miller_rabin(sqrt(num),k)){
ans *= 3;
}
else if (num != 1){
ans *= 4;
}
return ans;
}
int main(){
sieve_of_eratosthenes();
cout << "Number of divisors of 6 are " << divisors(6) << endl;
cout << "Number of divisors of 12 are " << divisors(12) << endl;
cout << "Number of divisors of 13 are " << divisors(13) << endl;
cout << "Number of divisors of 28 are " << divisors(28) << endl;
cout << "Number of divisors of 67 are " << divisors(67) << endl;
cout << "Number of divisors of 456 are " << divisors(456) << endl;
return 0;
}
Output :
Number of divisors of 6 are 4 Number of divisors of 12 are 6 Number of divisors of 13 are 2 Number of divisors of 28 are 6 Number of divisors of 67 are 2 Number of divisors of 456 are 16
PS :
This trick was directly used at one of the ACM ICPC regionals in the year 2015. See Problem F | Codeforces Gym . For practice purposes, you can also try it out on a problem which can be solved in sqrt(N). It will work like a charm.
