ucchasgfrank
ucchasgfrank

Reputation: 11

Efficiently Counting Divisors of a Number in O(n^(1/3)) in c++

I found an efficient algorithm to count divisors of a number in O(n^(1/3)) here: https://codeforces.com/blog/entry/22317

I've implemented the algorithm in the code below. But it doesn't work properly.

#include <bits/stdc++.h>
using namespace std;

#define ll long long
#define pb push_back
#define endl "\n"
#define pi 2*acos(0)

#define debug(s) cout << #s << " = " << s << endl
#define all(v) (v).begin(), (v).end()
#define mem(a, val) memset(a, val, sizeof a)

vector<int> prime(10000000, 1); // array size is 10^7
void sieve(ll n){
    prime[0]=prime[1]=0;
    for(ll i=2; i<=sqrt(n); i++){
        if(prime[i]==1){
            for(ll j=2; i*j<=n; j++) prime[i*j]=0;
        }
    }
}
bool isPrime(ll n){
    if(n<2) return false;
    return prime[n];
}

vector<int> primeTable;
ll numberOfDivisors(ll n){
    for(int i=0; i<1000000; i++){
        if(prime[i]==1) primeTable.pb(i);
    }
    ll ans=1;
    for(int i=0; ; i++){
        if(primeTable[i]*primeTable[i]*primeTable[i]>n) break;
        ll cnt=1;
        while(n%primeTable[i]==0){
            n=n/primeTable[i];
            cnt++;
        }
        ans*=cnt;
    }
    double sqn=sqrt(n);
    if(isPrime(n)) ans*=2;
    else if((sqn-floor(sqn))==0 and isPrime((ll)sqn)) ans*=3;
    else if(n!=1) ans*=4;

    return ans;
}

int main()
{
    ios_base::sync_with_stdio(false);cin.tie(NULL);
    ll n;
    cin >> n;
    ll N=cbrt(n);
    sieve(N);
    cout << numberOfDivisors(n) << endl;

    return 0;
}

i.e for 15 it returns 2(should be 4), for 100 it returns 6(should be 9), for 999 it returns 8(which is correct), but again for 9999 it returns 6(should be 12). I think something is wrong in between line 42 and 45. But I'm unable to find the bug. Please help me find the bug.

Upvotes: 0

Views: 2228

Answers (1)

nightfury1204
nightfury1204

Reputation: 4654

You set initial value of prime is 1 here vector<int> prime(10000000, 1). Then update the value of prime upto n in seive(ll n) function. So for n+1 to 10000000 prime value will remain 1. In your main function, you ran the sieve(N) for ll N=cbrt(n). For this reason your isPrime(n) is incorrect, because n could be greater than N=cbrt(n). Fix that your code is good to go.

From the blog post, i understand that this algorithm is for finding factors of long long integer. They also mention that to check primality use Miller Rabin. The way you check isPrime(n) will not work for large integer.

Also update this part (change 1000000 to cbrt(n)):

ll numberOfDivisors(ll n){
    for(int i=0; i<1000000; i++){
        if(prime[i]==1) primeTable.pb(i);
    }

Upvotes: 1

Related Questions