Fast primality checking, factorizing and totient calculation in Competitive Programming*

*This post only works for n \leq10^{18} (which I guess is enough since it’s seldom to have problem which requires BigInteger to save the input)

I am going to explain how to solve phi-phi-phi which simultaneously covers in how to do primality checking and totient calculation. The problem is asking to calculate totient from a value. The value of variable k itself can be ignored since the number of repetition in calculating the totient again and again is not very large until reach \phi(1) .

First, we are going to do sieve until 10^{6} first. Since doing sieve is quite common already, you can check the reference in Sieve Methods : Prime, Divisor, Euler Phi etc. . After doing the sieve, we are going to do calculate phi in k times or until it reaches 1.

For a value n, we are going to factorize the n using primes we have from sieve. The number of primes ranging from 1 \leq p \leq 10^{6} is 78498. After we factorize the value using the primes, there are 3 cases for the leftover value (let it be r ) (if the value is not equal to 1):

  1. r = p^2 , where p is a prime, 10^6 \leq p \leq 10^9.
  2. r = p_1 * p_2 , where p_1 and p_2 are primes, 10^6 \leq p_1, p_2 \leq 10^9 and p_1 \neq p_2
  3. r = p , where p is a prime, 10^6 \leq p \leq 10^{18}

It is easy to check for the first test case, we can take  \lfloor \sqrt r \rfloor and check if the square is equal to r. O(log(r)) .

To check for the third case, we can use Miller–Rabin primality test. Using the deterministic variants which requires the use of first 12 primes, we can check whether it’s prime or not in O(k * log^3(r)) , in this case k = 12 .

The second test case requires us to find the value of p1 and p2, we can use Pollard’s rho algorithm. The problem is only in calculating the function g(x) . Since g(x) use a square function, where x < n \leq 10^{18} , it can cause overflow. A simple multiplication modular such as (x * x) % n won't be good enough since the value of x * x has pass the limit of 64-bit integer. The technique to overcome overflow is using Russian Multiplication. The complexity of Russian Multiplication is O(log(min(a, b)) . Russian Multiplication also can be applied in calculating power function which results in O(log(b)) time and O(1) space. Overall: O(n^{1/4})

The function for calculating totient is: \phi(n) = \prod_{i=1}^{k} (p_i^{{\alpha}_i} - p_i^{{\alpha}_i - 1}) , where n = \prod_{i=1}^{k} p_i^{{\alpha}_i} or using Inclusion–exclusion principle. (I use the second one)

The following code is my submission:

#include <cstdio>
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <string>
#include <stack>
#include <queue>
#include <deque>
#include <vector>
#include <map>
#include <set>
#include <utility>
#include <algorithm>
#include <cmath>
#include <climits>
#ifdef DEBUG
    #include <ctime>
#endif
using namespace std;

// template

// abbreviations

typedef unsigned long long ull;
typedef long long ll;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef vector<ll> vl;
typedef vector<vl> vvl;
typedef vector<string> vs;
typedef pair<int, int> ii;
typedef vector<ii> vii;
typedef map<int, int> mii;
#define a first
#define b second
#define que queue
#define pque priority_queue
#define stk stack
#define pub push_back
#define pob pop_back
#define puf push_front
#define pof pop_front
#define pu push
#define po pop
#define mp make_pair
#define it iterator
#define sz(var) ((int) var.size())
#define rep(it, n) for(int it = 0; it < n; ++it) #define dep(it, n) for(int it = n - 1; it >= 0; --it)
#define rep1(it, n) for(int it = 1; it <= n; ++it) #define dep1(it, n) for(int it = n; it > 0; --it)
#define loop(it, from, to) for(int it = (from); it <= (to); ++it)
#define iter(it, cont) for(__typeof((cont).begin()) it = (cont).begin(); it != (cont).end(); ++it)
#define riter(it, cont) for(__typeof((cont).rbegin()) it = (cont).rbegin(); it != (cont).rend(); ++it)
#define all(cont) (cont).begin(), (cont).end()
#define rng(cont, n) cont, cont + n
#define memclr(var) memset(var, 0, sizeof(var))

const int INF = INT_MAX;
const int NINF = INT_MIN;
const ll INF_LL = LLONG_MAX;
const ll NINF_LL = LLONG_MIN;
const double PI = acos(-1.0);
const int MOD = 1e9 + 7;

inline ll pos_m(ll a, ll c = MOD) { while (a < 0) { a += c; } return a; } inline ll add_m(ll a, ll b, ll c = MOD) { return (a + b) % c; } inline ll mul_m(ll a, ll b, ll c = MOD) { ll ret = 0; while (b) { if (b & 1) ret = add_m(ret, a, c); a = add_m(a, a, c); b >>= 1; } return ret; }
inline ll sub_m(ll a, ll b, ll c = MOD) { return pos_m((a - b) % c, c); }
inline ll pow_mod(ll a, ll b, ll c = MOD) { ll ret = 1; while (b) { if (b & 1) ret = mul_m(a, ret, c); a = mul_m(a, a, c); b >>= 1; } return ret; }

#ifdef DEBUG
    #define debug(fmt, args...) printf("Line %d, in %s\t: " fmt, __LINE__, __FUNCTION__, ##args)
    #define rep_rt() printf("[Run time: %.3fs]\n", ((double) clock()) / CLOCKS_PER_SEC)
#else
    #define debug(...)
#endif

// end of template

#define MAXSQRT3A (int) (1e6)
int is_prime[MAXSQRT3A + 1];
vi primes;

bool is_prime_fun(ll n) {
    int s = 0;
    ll d = n - 1;
    while ((d % 2) == 0) {
        s++;
        d >>= 1;
    }
    for (int i = 0; i < 12; ++i) {
        int prime = primes[i];
        ll fcond = pow_mod(prime, d, n);
        if (fcond == 1)
            continue;
        bool iscomp = true;
        for (int r = 0; r < s; ++r) {
            if (pow_mod(fcond, 1LL << r, n) == n - 1) {
                iscomp = false;
                break;
            }
        }
        if (iscomp)
            return false;
    }
    return true;
}

ll pr_g(ll x, ll n) {
    return add_m(pow_mod(x, 2, n), 1, n);
}

ll pr_get_prime_fact(ll val) {
    for (int i = 0; i < 12; ++i) {
        int prime = primes[i];
        ll x = prime, y = prime, d = 1;
        while (d == 1) {
            x = pr_g(x, val);
            y = pr_g(pr_g(y, val), val);
            d = __gcd(abs(x - y), val);
        }
        if (d != val)
            return d;
    }
    return -1;
}

ll tot_func(vl divs, ll lim) {
    ll ret = 0;
    for (int bm = 0; bm < (1 << sz(divs)); ++bm) {
        ll tot_div = 1;
        for (int i = 0; i < sz(divs); ++i) {
            if (bm & (1 << i))
                tot_div *= divs[i];
        }
        ll sans = lim / tot_div;
        if (__builtin_popcount(bm) % 2 == 0) {
            ret += sans;
        } else {
            ret -= sans;
        }
    }
    return ret;
}

ll sqrt_ll(ll val) {
    ll lb = 0, ub = (int) (1e9);
    while (lb < ub) {         ll m = (lb + ub) >> 1;
        ll sqm = m * m;
        if (sqm < val) {
            lb = m + 1;
        } else {
            ub = m;
        }
    }
    return lb;
}

int main() {
#ifdef DEBUG
    freopen("phi-phi-phi.in", "r", stdin);
#endif

    rep1(val, MAXSQRT3A) {
        is_prime[val] = true;
    }
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; i * i <= MAXSQRT3A; ++i) {
        if (is_prime[i]) {
            for (int j = i; i * j <= MAXSQRT3A; ++j)
                is_prime[i * j] = false;
        }
    }
    rep1(val, MAXSQRT3A) {
        if (is_prime[val])
            primes.pub(val);
    }
    debug("%d\n", sz(primes));

    ll n, k;
    scanf("%lld %lld", &n, &k);
    while (k--) {
        ll lim = n;

        vl divs;
        iter(pp, primes) {
            int p = *pp;
            if (n % p == 0) {
                while (n % p == 0)
                    n /= p;
                divs.pub(p);
                debug("%d\n", p);
            }
        }
        if (n != 1) { // means n is a prime, sq prime or p1 * p2
            ll sqp = sqrt_ll(n);
            if (sqp * sqp == n)
                divs.pub(sqp);
            else if (is_prime_fun(n)) {
                divs.pub(n);
            } else {
                ll p1 = pr_get_prime_fact(n);
                ll p2 = n / p1;
                divs.pub(p1);
                divs.pub(p2);
            }
        }
        n = tot_func(divs, lim);
        if (n == 1)
            break;
    }
    printf("%lld\n", n);

#ifdef DEBUG
    rep_rt();
#endif
    return 0;
}
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s