Divisor Sum Implementations

#	User	Rating
1	tourist	3985
2	jiangly	3814
3	jqdai0815	3682
4	Benq	3529
5	orzdevinwang	3526
6	ksun48	3517
7	Radewoosh	3410
8	hos.lyric	3399
9	ecnerwala	3392
9	Um_nik	3392

#	User	Contrib.
1	cry	169
2	maomao90	162
2	Um_nik	162
4	atcoder_official	161
5	djm03178	158
6	-is-this-fft-	157
7	adamant	155
8	awoo	154
8	Dominater069	154
10	luogu_official	150

I found these implementations ^^

Sorry, I am not sure if there are some corner cases or I calculated wrong complexity. If there are, please correct me.

I am wondering if there are some other implementations ^^ (you can comment below and I will add to the post)

Implementation 0: Trivial

/// Implementation: Trivial
/// Time  Complexity: O(n)
/// Space Complexity: O(1)

int main()
{
    int n;
    cin >> n;

    ll sum = 0;
    /// n = sigma(x | n % x == 0)
    for (int i = 1; i <= n; ++i) /// O(n)
        if (n % i == 0) /// If (i) is a divisor
            sum += i;
    
    cout << sum;
    return 0;
}

Implementation 1: Optimized Trivial

/// Implementation: Optimized Trivial
/// Time  Complexity: O(√n)
/// Space Complexity: O(1)

typedef long long ll;
int main()
{
    ll n;
    cin >> n;
    int sqrtn = sqrt(n);
    
    ll sum = 0;
    /// if (i) is a divisor then (n) / (i) is also a divisor
    for (ll i = 1; i <= sqrtn; ++i) /// O(√n)
    {
        if (n % i != 0) continue; /// If (i) not a divisor
        sum += i;
        ll t = n / i;
        if (i != t) sum += t;     /// If (t) is a unique divisor
    }
    
    cout << sum;
    return 0;
}

Implementation 2: Naive Factorization

/// Implementation: Naive | Multiplicative Formula
/// Time  Complexity: O(n ^ 2)
/// Space Complexity: O(1)

typedef long long ll;
ll pw(int x, int n) /// O(n)
{
    ll res = 1;
    for (int i = 1; i <= n; ++i) res *= x;
    return res;
}

bool isPrime(int n) /// O(n)
{
    int cnt = 0;
    for (int i = 1; i <= n; ++i) cnt += (n % i == 0);
    return (cnt == 2);
}

int main() /// O(n ^ 2) ?
{
    int n;
    cin >> n;

    ll sum = 1;
    /// n = p1 ^ f1 * p2 ^ 2 * ... * pk ^ fk
    /// sum = ∏(divisor prime p){ p ^ (f + 1) - 1) / (p - 1) }
    for (int i = 1; i <= n; ++i) /// O(n * f(n)) ≤ O(n + divisor_count ^ 2) ≤ O(n ^ 2)
    {
        if (n % i != 0) continue;
        if (isPrime(i) == false) continue; /// O(i)
        
        int p = i, f = 0;
        do f++, n /= p; while (n % p == 0);  /// O(log_p(n))
        sum *= (pw(p, f + 1) - 1) / (p - 1); /// O(f)
    }
    
    cout << sum;
    return 0;
}

Implementation 3: Factorization

/// Implementation: Optimized | Multiplicative Formula
/// Time  Complexity: O(n√n)
/// Space Complexity: O(1)

typedef long long ll;
ll pw(ll x, ll n) /// O(log n)
{
    ll res = 1;
    for (; n > 0; x = x * x, n >>= 1)
        if (n & 1) res = res * x;
    
    return res;
}

bool isPrime(ll n) /// O(sqrt n)
{ 
    if (n < 2) return false; 
    if (n < 4) return true; 
    if (n % 2 == 0 || n % 3 == 0) return false; 
  
    int sqrtn = sqrt(n);
    for (ll i = 5; i <= sqrtn; i += 6) 
        if ((n % i == 0) || (n % (i + 2) == 0)) 
            return false; 
  
    return true; 
} 

int main() /// O(n√n)
{
    ll n;
    cin >> n;

    ll sum = 1;
    /// n = p1 ^ f1 * p2 ^ 2 * ... * pk ^ fk
    /// sum = ∏(divisor prime p){ p ^ (f + 1) - 1) / (p - 1) }
    for (ll i = 1; i <= n; ++i) /// O(n * f(n)) ≤ O(n + sqrt(divisor_count)) ≤ O(n * sqrt(n))
    {
        if (n % i != 0) continue;
        if (isPrime(i) == false) continue; /// O(sqrt i)  

        ll p = i; int f = 0;
        do f++, n /= p; while (n % p == 0);  /// O(log_p(n))
        sum *= (pw(p, f + 1) - 1) / (p - 1); /// O(log f)
    }
    
    cout << sum;
    return 0;
}

Implementation 4: Miller-rabin

/// Implementation: Miller-rabin Primality Test | Multiplicative Formula
/// Time  Complexity: O(n * log^5(n))
/// Space Complexity: O(1)

typedef long long ll;
ll bigmodMul(ll a, ll b, ll m = MOD) /// O(log n)
{
    ll res = 0;
    for (a %= m, b %= m; b > 0; a <<= 1, b >>= 1) {
        if (a >= m) a -= m;
        if (b & 1) { res += a; if (res >= m) res -= m; }
    }

    return res;
}
ll bigmodPow(ll x, ll n, ll m = MOD) /// O(log^2(n))
{
    ll res = 1;
    for (x %= m; n > 0; x = bigmodMul(x, x, m), n >>= 1)
        if (n & 1) res = bigmodMul(res, x, m);
    
    return res;
}

ll pw(ll x, ll n) /// O(log n)
{
    ll res = 1;
    for (; n > 0; x = x * x, n >>= 1)
        if (n & 1) res = res * x;
    
    return res;
}

bool isPrime(ll p) /// O(log^5(p)) - O(log p) rounds performed - O(log^2(p)) primality tests - O(log^2(p)) calculation
{
    if (p < 2) return false;
    if (p < 4) return true;
    if (p % 2 == 0 || p % 3 == 0) return false;
 
    ll q = p - 1;
    int k = 0;
    while ((q & 1) == 0) /// O(log p)
        q >>= 1, k++;
 
    ll a = rand() % (p - 4) + 2;
    ll t = bigmodPow(a, q, p); /// O(log^2(q))
 
    bool ok = (t == 1) || (t == p-1);
    for (int i = 1; i <= k && !ok; i++) { /// O(log k * log t)
        t = bigmodMul(t, t, p); /// O(log t)
        ok = (t == p - 1);
    }
 
    return ok;
}

int main() /// O(n * log5(n))
{
    ll n;
    cin >> n;
    ll sum = 1;
    /// n = p1 ^ f1 * p2 ^ 2 * ... * pk ^ fk
    /// sum = ∏(divisor prime p){ p ^ (f + 1) - 1) / (p - 1) }
    for (ll i = 1; i <= n; ++i) /// O(n * f(n)) ≤ O(divisor_count * log^5(n)) ≤ O(n * log^5(n))
    {
        if (n % i != 0) continue;
        if (isPrime(i) == false) continue; /// O(log^5(i))

        ll p = i; int f = 0;
        do f++, n /= p; while (n % p == 0);  /// O(log_p(n))
        sum *= (pw(p, f + 1) - 1) / (p - 1); /// O(log f)
    }
    
    cout << sum;
    return 0;
}

Implementation 5: Sieve + Factorization

/// Implementation: Sieve | Factorization | Multiplicative Formula
/// Time  Complexity: O(n) precalculation + O(log n) query
/// Space Complexity: O(n) sieve + O(1) query

typedef long long ll;
ll pw(ll x, ll n) /// O(log n)
{
    ll res = 1;
    for (; n > 0; x = x * x, n >>= 1)
        if (n & 1) res = res * x;
    
    return res;
}

vector<int> prime; /// List of prime numbers
vector<int> lpf;   /// lpf[x] = Lowest_Prime_Factor of x

void sieve(int lim = LIM) /// O(n)
{
    prime.assign(1, 2);
    lpf.assign(lim + 1, 2);

    for (int i = 3; i <= lim; i += 2)
    {
        if (lpf[i] == 2) prime.pb(lpf[i] = i);
        for (int j = 0; j < sz(lpf) && prime[j] <= lpf[i] && prime[j] * i <= lim; ++j)
            lpf[prime[j] * i] = prime[j];
    }
}

int main()
{
    int n;
    cin >> n;
    sieve(n); /// O(n) for sieving

    ll sum = 1;
    /// n = p1 ^ f1 * p2 ^ 2 * ... * pk ^ fk
    /// sum = ∏(divisor prime p){ p ^ (f + 1) - 1) / (p - 1) }
    for (int p, f; n > 1; ) /// O(log n)
    {
        for (p = lpf[n], f = 0; (n > 1) && (p == lpf[n]); n /= p) f++; /// O(log n)
        sum *= (pw(p, f + 1) - 1) / (p - 1);                           /// O(log n)
    }
    
    cout << sum;
    return 0;
}

Implementation 6: Sieve + Miller-rabin + Pollard-rho

/// Implementation: Sieve | Factorization | Pollard-rho | Miller-rabin | Multiplicative Formula
/// Time  Complexity: O(√n) precal + O(√(√(n)) polylog n) query
/// Space Complexity: O(√n) sieve + O(log n / log(log n)) query

const int MOD = 1e9 + 7;
typedef long long ll;
ll bigmodMul(ll a, ll b, ll m = MOD) /// O(log n)
{
    ll res = 0;
    for (a %= m, b %= m; b > 0; a <<= 1, b >>= 1) {
        if (a >= m) a -= m;
        if (b & 1) { res += a; if (res >= m) res -= m; }
    }

    return res;
}
ll bigmodPow(ll x, ll n, ll m = MOD) /// O(log^2(n))
{
    ll res = 1;
    for (x %= m; n > 0; x = bigmodMul(x, x, m), n >>= 1)
        if (n & 1) res = bigmodMul(res, x, m);
    
    return res;
}

ll pw(ll x, ll n) /// O(log n)
{
    ll res = 1;
    for (; n > 0; x = x * x, n >>= 1)
        if (n & 1) res = res * x;
    
    return res;
}

bool isPrime(ll p) /// O(log^5(p)) - O(log p) rounds performed - O(log^2(p)) primality tests - O(log^2(p)) calculation
{
    if (p < 2) return false;
    if (p < 4) return true;
    if (p % 2 == 0 || p % 3 == 0) return false;
 
    ll q = p - 1;
    int k = 0;
    while ((q & 1) == 0) /// O(log p)
        q >>= 1, k++;
 
    ll a = rand() % (p - 4) + 2;
    ll t = bigmodPow(a, q, p); /// O(log^2(q))
 
    bool ok = (t == 1) || (t == p-1);
    for (int i = 1; i <= k && !ok; i++) { /// O(log k * log t)
        t = bigmodMul(t, t, p); /// O(log t)
        ok = (t == p - 1);
    }
 
    return ok;
}

map<int, int> M;
int sqrtn;
ll sum = 0; 

vector<int> prime; /// List of prime numbers
vector<int> lpf;   /// lpf[x] = Lowest_Prime_Factor of x

void sieve(int lim = LIM) /// O(n)
{
    prime.assign(1, 2);
    lpf.assign(lim + 1, 2); /// O(n)

    for (int i = 3; i <= lim; i += 2) /// O(n)
    {
        if (lpf[i] == 2) prime.pb(lpf[i] = i);
        for (int j = 0; j < sz(lpf) && prime[j] <= lpf[i] && prime[j] * i <= lim; ++j)
            lpf[prime[j] * i] = prime[j];
    }
}

ll rho(ll n, ll c) {  /// O(√(√(n)) * polylog n)
    ll x = 2, y = 2, i = 2, k = 2, d;
    while (true) {
        x = (bigmodMul(x, x, n) + c); /// O(log n)
        if (x >= n)	x -= n;
        d = __gcd(abs(x - y), n);     /// O(log n)
        if (d > 1) return d;
        if (i++ == k) y = x, k <<= 1;
    }
    return n;
}
 
void big_fact(ll n) { /// O(√(√(n)) * polylog n)
    if (n == 1) /// O(1)
        return ;
 
    if (n < 1e+9) { /// O(√(n) / ln(√(n)) * log n / log log n)
        for (int p : prime) /// O(pair<int, int>(√(n)) * log(M.size)) = O(√(n) / ln(√(n)) * log n / log(log n))
        {
            if (p > sqrtn) break;
            while (n % p == 0) M[p]++, n /= p; /// O(log_p(n) * log(M.size)) = O(log_p(n) * log(n) / log(log n))
        }

        if (n > 1) M[n]++; /// O(log(M.size)) = O(log n / log(log n))
        return ;
    }
 
    if (isPrime(n))            /// O(log^5(n))
        return M[n]++, void(); /// O(log(M.size)) = O(log n / log(log n))
 
    ll d = n;
    for (int i = 2; d == n; i++) /// O(√(√(n)) * polylog n)
        d = rho(n, i);           /// O(√(√(n)) * polylog n)
 
    big_fact(    d);
    big_fact(n / d);
}

int main() /// O(√n + √(√(n)) * polylog n)
{
    ll n;
    cin >> n;

    sqrtn = sqrt(n);   /// O(√(n))   
    sieve(sqrtn);      /// O(√(n))
    
    srand(time(NULL));
    big_fact(n);       /// O(√(√(n)) * polylog n)
    
    ll sum = 1;
    /// n = ∏(e ∈ M){ e.first ^ e.second }
    /// sum = ∏(divisor prime e.first){ e.first ^ (e.second + 1) - 1) / (e.first - 1) }
    for (pair<int, int> e : M)     /// O(log n / log(log n))
        sum *= (pw(e.first, e.second + 1) - 1) / (e.first - 1); /// O(log n)

    cout << sum;    
    return 0;
}

Implementation 7: Divisor Sum Sieve

vector<int> prime;
vector<int> ds; /// divisor sum
void sieve(int lim) /// Precalculation: O(n log n)
{
    prime.clear();
    if (n < 2) return ;

    ds.assign(n + 1, 1);
    for (int i = 2; i <= lim; ++i)
        for (int j = i; j <= lim; j += i)
            ds[j] += i;
    
    return ; /// If you want to get primes, remove this line
    prime.assign(1, 2);
    for (int i = 3; i <= n; i += 2)
        if (ds[i] == i + 1)
            prime.push_back(i);
}

int main()
{
    int n;
    cin >> n;
    sieve(n);      /// O(n) for sieving
    cout << ds[n]; /// O(1) query
    return 0;
}

Compilation:

Implementation	Main Algorithm	Time Complexity	Space Complexity	Coding Space	Coding Complex
0. Trivial	Implementation	O(n)	O(1)	Short	Simple
1. Trivial	Implementation	O(√n)	O(1)	Short	Simple
2. Factorization + Math	Math	O(n ^ 2) ?	O(1)	Normal	Normal
3. Factorization + Math	Factorization	O(n√n)	O(1)	Normal	Normal
4. Factorization + Math	Miller-rabin	O(n * log^5(n))	O(1)	Long	Normal
5. Factorization + Math	Sieve	O(n) + O(log n)	O(n)	Normal	Normal
6. Factorization + Math	Pollard-rho	O(√n) + O(√(√(n)) polylog n) query	O(√n) + O(log n / log(log n)) query	Very Long	Complex

Planning:

Add relative blogs
Add some more implementations
Add tutorial & reasoning for each implementation

Comments (2)

Show archived | Write comment?

feodorv

5 years ago, # |

You can add one simple sieve for divisor summation for all numbers from interval $$$[1 \dots N]$$$. For example, this is the solution for SPOJ problem DIVSUM:

Code

#include <stdio.h>

#define MAXN 500000

int sum[MAXN+1];

void sieve( void )
{
  int i, j;
  for( i = 1; i <= MAXN; ++i)
    for( j = i+i; j <= MAXN; j += i) sum[j] += i;
}

int main( void )
{
  int t;
  sieve();
  scanf( "%d", &t);
  while( t-- > 0 )
  {
    int n;
    scanf( "%d", &n);
    printf( "%d\n", sum[n]);
  }
  return 0;
}

Also the main method for summation of the divisors of the single (not very large) number $$$N$$$ is prime factorization with ready prime list till $$$\sqrt{N}$$$. So you can easily solve SPOJ problem DIVSUM2:

#include <stdio.h>

#if defined( _WIN32 )
typedef __int64 az_int64_t;
typedef unsigned __int64 az_uint64_t;
#define I64(x) x ## I64
#define F64 "I64"
#else
typedef long long az_int64_t;
typedef unsigned long long az_uint64_t;
#define I64(x) x ## ll
#define F64 "ll"
#endif

#define LIMIT 100000000
#define HAS(_n) (pbits[(_n)>>5] & (1u << ((_n) & 0x1f)))
#define SET(_n) (pbits[(_n)>>5] |= (1u << ((_n) & 0x1f)))
unsigned int pbits[LIMIT/32+1];
int pc, pr[LIMIT/10];

void sieve( void )
{
  int i, j;
  pc = 0;
  for( i = 2; i <= LIMIT; ++i)
  {
    if( !HAS( i ) ) pr[pc++] = i;
    for( j = 0; j < pc && i * pr[j] <= LIMIT; ++j)
    {
      int num = i * pr[j];
      SET( num );
      if( i % pr[j] == 0 ) break;
    }
  }
}

az_int64_t solution( az_int64_t n )
{
  az_int64_t o = n, ans = 1;
  int i;
  for( i = 0; i < pc && (az_int64_t) pr[i] * pr[i] <= n; ++i)
  {
    int count = 0;
    az_int64_t sum = 1, pw = 1;
    while( n % pr[i] == 0 )
    {
      pw *= pr[i];
      sum += pw;
      count++;
      n /= pr[i];
    }
    if( count > 0 ) ans *= sum;
  }
  if( n > 1 ) ans *= n + 1;
  return ans - o;
}

int main( void )
{
  int t;
  sieve();
  scanf( "%d", &t);
  while( t-- > 0 )
  {
    az_int64_t n;
    scanf( "%" F64 "d", &n);
    printf( "%" F64 "d\n", solution( n ));
  }
  return 0;
}

PS Please, pay attention that the SPOJ problems in question are about summation of proper divisors.

→ Reply

SPyofgame

5 years ago, # ^ |

← Rev. 4 →

Thank you for the solutions ^^

SPyofgame's blog

I found these implementations ^^

Compilation:

Planning: