Euler Totient / φ Function

Euler Totient / Phi Function φ(n) counts the number of positive integers less than or equal to n which are relatively prime to n, i.e. do not have any common divisor with n except 1.

Formula for Euler Phi function:

Here, the product is over the distinct prime numbers which divide n. Now, you can just factorize n and calculate φ(n) pretty easily. But, will that be efficient for a task such as you are asked to find φ(n) over a range of integers?

If you look closely to the formula, you will see that, we multiply n with (p-1)/p for each prime p that divides n. Now recall what do we do when we run Sieve of Eratosthenes for marking primes / non-primes. On the outer loop of the sieve, we determine if a number is prime, then in inner loop, instead of setting flags, if we can keep multiplying the number with (p-1)/p where p is the current prime number from outer loop, at the end of the iterations, we can actually generate φ(n) for each n over the range we are asked for. Here is a source code example that does exactly the same thing. Just take a look and try to understand what are the loops doing here, and how we are performing calculation and storing results.

#include <cstdio>

const int MAX = 1000001;
int phi[MAX];

void euler_phi() {
    phi[1] = 1;
    for(int i=2; i<MAX; i++) {
        if(!phi[i]) {
            phi[i] = i-1;
            for(int j=(i<<1); j<MAX; j+=i) {
                if(!phi[j]) phi[j] = j;
                phi[j] = phi[j]/i*(i-1);
            }
        }
    }
}

int main() {
    euler_phi();
    for(int t=1; t<MAX; t++) printf("%d\n", phi[t]);
    return 0;
}

A bit of explanation on what we are doing here: Initially the phi[] array is set to 0 (as it is declared global). We know that, phi[1] = 1 and phi[n] = n-1 when n is a prime number. So, similar to sieve algorithm, first we check if current number is prime in the outer loop, if phi[i] = 0, it means i is prime. So, we update it with i-1 accordingly. Now, for all the multiples of i greater than i, which starts from 2*i, calling it j in inner loop, we need to multiply phi[j] by (i-1) / i. Here, we first check if phi[j] is 0, i.e. visiting it for the first time, in this case we set it with j, as I said earlier that, for φ(n) we will multiply n with (p-1)/p where p are the primes that divide n. Also, one thing to note here: a * b / c and a / c * b are not always same for integer calculation unless you can assure that c divides a. In this case it does, why? cause this is basically a prime factoring algorithm, and c = i here. As we are traversing i's multiples, it is guaranteed that a=phi[j] can be divided by c=i and instead of a * b / c format, we will always use a / c * b in these types of situations, because it will help preventing overflow many times.

Now, think about the optimizations we could apply here, and try applying them, like discarding even numbers, starting inner loop from squares, increment inner loop twice the prime number each time, won't work here. Because, we need to go through every numbers in inner loop, as we are trying to find Totient function for every n in the range 1 to MAX. Test the code on smaller range, and try to check if it is doing this correctly, play around with it.

Divisor Function

Prime Factorization and Divisor Function σ_x(n)

σ_x(n) is defined as the sum of x^th powers of the distinct positive divisors of n. The function can be expressed as:

Here, r = ω(n), which is the number of distinct positive prime factors of n. p_i for i = 1 to r, are the prime factors and a_i is the maximum power of p_iby which n is divisible.

Clearly, this function can be used for various problems, for example, when x = 0, it simply means the number of distinct positive divisors of n, if x = 1, it is the sum of distinct positive divisors of n, for x = 2, its the sum of squares of positive divisors of n and so on.

For programming contest practice, there are a few problems that requires sum of divisors, or number of divisors, which can be calculated by simply calculating the prime factors with their count, and then evaluating the function shown above with appropriate value on x. Also, the form of function definition can be changed when you set a value for x. After putting x = 0, we get this:

So this is easy, you just need to find frequency of each prime, and then multiply each (a_i + 1) for all primes, you get the number of divisors of n.

For sum of divisor, the idea is similar, here x = 1. The following code shows how to find sum of distinct positive divisors of n:

#include <cstdio>
#include <cmath>
using namespace std;

#define sq(x) ((x)*(x))
#define i64 unsigned long long
#define MAX 784
#define LMT 28

unsigned flag[MAX/64];
unsigned primes[5761460], total;

#define chkC(n) (flag[n>>6]&(1<<((n>>1)&31)))
#define setC(n) (flag[n>>6]|=(1<<((n>>1)&31)))

/*
Regular sieve of eratosthenes, bitwise implementation
*/
void sieve()
{
    unsigned i, j, k;
    flag[0]|=0;
    for(i=3;i<LMT;i+=2)
        if(!chkC(i))
            for(j=i*i,k=i<<1;j<MAX;j+=k)
                setC(j);
    primes[(j=0)++] = 2;
    for(i=3;i<MAX;i+=2)
        if(!chkC(i))
            primes[j++] = i;
    total = j;
}

/*
finds n^p in log(p) time
*/
i64 power(unsigned n, unsigned p)
{
    i64 x=1, y=n;
    while(p > 0)
    {
        if(p&1) x *= y;
        y *= y;
        p >>= 1;
    }
    return x;
}

/*
calculates the sigma(1) function, we don't need to find prime frequencies.
*/
inline void update(i64 &sigma1, i64 n, unsigned p)
{
    if(p==1) sigma1 *= (n+1);
    else sigma1 *= ((power(n,p+1)-1)/(n-1));
}

/*
Factorization function, we do not need to store the primes here,
instead, whenever a prime is found, you update corresponding prime
and frequency with sigma 1
*/
void factor(i64 n, i64 &sigma1)
{
    unsigned i, v;
    i64 t;
    for(i=0, t=primes[i]; i<total && t*t <= n; t = primes[++i])
    {
        if(n % t == 0)
        {
            v = 0;
            while(n % t == 0)
            {
                v++;
                n /= t;
            }
            update(sigma1, primes[i], v);
        }
    }
    if(n>1) update(sigma1, n, 1);
}

/*
Our beloved main function
*/
int main()
{
    int t, x;
    i64 n, sigma1;
    sieve();
    scanf("%d", &t);
    for(x=1; x<=t; x++)
    {
        scanf("%llu", &n);
        factor(n, sigma1);
        printf("%llu\n",sigma1);
    }
    return 0;
}

We just need the values of σ, so here we will not store the prime factors. Some similar problem would be to find the number of odd divisors of n, or testing if a number is square free, i.e. no perfect square divides n, going to leave these as an exercise for readers.

A closer look to σ function from wikipedia.

Hacker Cup 2013 Round 1

Facebook Hacker Cup 2013 Round 1

Facebook hacker cup 2013: Round 1 has been finished, and the results are out as well. The problem set will be found here.

Problem A: Card Game 20 pts problem. If you take a little bit time observing the results, you can easily see what is really happening here. Any number X can be maximum in a set if all the other numbers are less or equal to X, so if we have at least K-1 such numbers, we will be able to make K sized subset where X is the maximum. If we can make M such subsets, X will be added to result M times.

Now, say, for X, there are P numbers which are less or equal to X and P >= K-1. So, we can choose K-1 numbers from P numbers in C(P, K-1) ways, everyone knows that. So, first sort all the numbers, and start from the Kth number (it will be pointless to take previous items, as they will not be maximum for any other subset). So, the Kth number will be added C(K-1, K-1) times, (K+1)th number will be added C(K, K-1) times, next one for C(K+1, K-1) times and so on while the last number will be added C(N-1, K-1) times, consider N as total number of integers.

So the solution is simple (Note, you will need modular inverse to properly count combinations, and the overall complexity is O(n lg n) ):

typedef __int64 i64;

const i64 MOD = 1000000007;
const int MAX = 10009;

i64 a[MAX];
i64 fact[MAX];

class Euclid {
public:
    i64 x, y, d;
    Euclid() {}
    Euclid(i64 _x, i64 _y, i64 _d) : x(_x), y(_y), d(_d) {}
};

Euclid egcd(i64 a, i64 b) {
    if(!b) return Euclid(1, 0, a);
    Euclid r = egcd(b, a % b);
    return Euclid(r.y, r.x - a / b * r.y, r.d);
}

i64 modinv(i64 a, i64 n) {
    Euclid t = egcd(a, n); // here n is always prime, no check needed
    i64 r = t.x % n;
    return (r < 0 ? r + n : r);
}

int main() {
    //freopen("card_game.txt", "r", stdin);
    //freopen("card_game.out", "w", stdout);
    int test, cs, n, k, i;
    i64 sum, num, deno;
    for(fact[0] = i = 1; i < MAX; i++) fact[i] = (fact[i-1] * i) % MOD;
    scanf("%d", &test);
    for(cs = 1; cs <= test; cs++) {
        scanf("%d %d", &n, &k);
        for(i = 0; i < n; i++) scanf("%I64d", &a[i]);
        sort(a, a + n);
        for(sum = 0, num = deno = fact[k-1], i = k-1; i < n;) {
            sum = (sum + (((a[i] * num) % MOD) * modinv(deno, MOD)) % MOD) % MOD;
            i++; num = (num * i) % MOD; deno = (deno * (i-k+1)) % MOD;
        }
        printf("Case #%d: %I64d\n", cs, sum);
    }
    return 0;
}

Problem B: Security 35 pts. I got a cross here :( Will talk about it after describing the solution. Basically two forms of the same string is given, and each one has M segments, all the segments have equal length. As we know, string 1 is actually the real string while string 2 is the permuted string (only whole segments are permuted). So we can run a maximum bipartite matching among the segments of two strings. If they do not match, the the answer is obviously possible. While comparing, if the current character of both string is equal, or any of them is a ? mark, then we also consider them to be equal.

The tricky part of this solution is, we can have multiple matches, and in those case we need to choose the match which will give us lexicographically smallest string. (I failed in this part). It can be achieved in many ways, like putting weights in matching and taking minimum weighted match, or just trying bipartite matching cutting every edge once and see if we get a better match. We then need to copy the letters from string 2 segments to string 1 segments where there is a ? in string 1. Finally, on the remaining ?, just put 'a'. Not posting the code, because it was not fully correct, but the idea is ok (confirmed by others who were able to solve it correctly).

Problem C: Dead Pixels 45 pts. Looks tricky at a first glance. Here, notable thing is, each problem has only one query, i.e. there is not dynamic update + query. So, we don't really need to do anything like segment tree here. Although, segment tree can be used to solve this one as well.

The solution is actually simple, we just sort the points by X axis, and run sliding window of width equal to picture width. We maintain a set of Y values in current window. When we skip a X value, we remove (if available) corresponding y values from the set and update number of possible positions, and then we add the y values (if available) of new X that we've just touched by our window and update accordingly. Here is my solution, pretty simple I should say.

const double EPS = 1e-9;
const int INF = 0x7f7f7f7f;

const int MAX = 1000032;

pii p[MAX];
int W, H, P, Q, N, X, Y, a, b, c, d;
map< int, int > M;
map< int, int >::iterator it;
set< int > S;
set< int > :: iterator Curr, Prev, Next;

void generate() {
    p[0].ff = X, p[0].ss = Y;
    for(int i = 1; i < N; i++) {
        p[i].ff = (p[i - 1].ff * a + p[i - 1].ss * b + 1) % W;
        p[i].ss = (p[i - 1].ff * c + p[i - 1].ss * d + 1) % H;
    }
}

int lowerbound(pii *a, int st, int en, int val) {
    int idx, cnt, stp;
    cnt = en - st;
    while(cnt > 0) {
        stp = cnt >> 1; idx = st + stp;
        if(a[idx].ff < val) st = ++idx, cnt -= stp+1;
        else cnt = stp;
    }
    return st;
}

int upperbound(pii *a, int st, int en, int val) {
    int idx, cnt, stp;
    cnt = en - st;
    while(cnt > 0) {
        stp = cnt >> 1; idx = st + stp;
        if(a[idx].ff <= val) st = ++idx, cnt -= stp+1;
        else cnt = stp;
    }
    return st;
}

inline int span(int lt, int rt) {
    return (rt-Q)-(lt+1)+1;
}

void insert(int cv, int &zeros) {
    int lt, rt, gap;
    it = M.find(cv);
    if(it == M.end()) {
        M.insert(pii(cv, 1));
        S.insert(cv);
        Next = Prev = Curr = S.find(cv);
        Prev--; Next++;
        lt = (Curr == S.begin()) ? -1 : *Prev;
        rt = (Next == S.end()) ? H : *Next;
        gap = span(lt, rt); if(gap > 0) zeros -= gap;
        gap = span(lt, cv); if(gap > 0) zeros += gap;
        gap = span(cv, rt); if(gap > 0) zeros += gap;
    }
    else it->ss++;
}

void remove(int cv, int &zeros) {
    int lt, rt, gap;
    it = M.find(cv); if(it == M.end()) return;
    if(it->ss == 1) {
        Next = Prev = Curr = S.find(cv);
        Prev--; Next++;
        lt = (Curr == S.begin()) ? -1 : *Prev;
        rt = (Next == S.end()) ? H : *Next;
        gap = span(lt, cv); if(gap > 0) zeros -= gap;
        gap = span(cv, rt); if(gap > 0) zeros -= gap;
        gap = span(lt, rt); if(gap > 0) zeros += gap;
        S.erase(Curr);
        M.erase(it);
    }
    else it->ss--;
}


int main() {
    //READ("dead_pixels.txt");
    //WRITE("dead_pixels.out");
    int test, cs, i, j, zeros, tmp, total, st, en;
    scanf("%d", &test);
    for(cs = 1; cs <= test; cs++) {
        scanf("%d %d %d %d %d %d %d %d %d %d %d", &W, &H, &P, &Q, &N, &X, &Y, &a, &b, &c, &d);
        generate();
        sort(p, p + N);
        M.clear(); S.clear();
        zeros = H - Q + 1; total = 0;
        tmp = upperbound(p, 0, N, P-1);
        for(i = 0; i < tmp; i++) {
            insert(p[i].ss, zeros);
        }
        for(i = 0; i <= W - P; i++) {
            total += zeros;
            st = lowerbound(p, 0, N, i);
            en = upperbound(p, 0, N, i);
            for(j = st; j < en; j++) {
                remove(p[j].ss, zeros);
            }
            st = lowerbound(p, 0, N, i+P);
            en = upperbound(p, 0, N, i+P);
            for(j = st; j < en; j++) {
                insert(p[j].ss, zeros);
            }
        }
        printf("Case #%d: %d\n", cs, total);
    }
    return 0;
}

However, could not make it to top 500 because of the time penalty and the mistake in problem B. Overall, problems were nice. And again, 45 pts one was easier than the 35 pts.