SPOJ GUANGGUN

Problem link: SPOJ Problem Set (classical): 9952. 111…1 Squared

Interesting and fun problem. This is one of those problems that you can sense a pattern from miles away. Just do some bruteforce and print the patterns. You can use pen and paper too, but Python is surely a better option when you have it. All I did is write some bruteforce in my python interpreter, and the pattern is instantly visible, here is the output of f(n)² for n = 1 to 50:

 1 121 12321 1234321 123454321 12345654321 1234567654321 123456787654321 12345678987654321 1234567900987654321 123456790120987654321 12345679012320987654321 1234567901234320987654321 123456790123454320987654321 12345679012345654320987654321 1234567901234567654320987654321 123456790123456787654320987654321 12345679012345678987654320987654321 1234567901234567900987654320987654321 123456790123456790120987654320987654321 12345679012345679012320987654320987654321 1234567901234567901234320987654320987654321 123456790123456790123454320987654320987654321 12345679012345679012345654320987654320987654321 1234567901234567901234567654320987654320987654321 123456790123456790123456787654320987654320987654321 12345679012345679012345678987654320987654320987654321 1234567901234567901234567900987654320987654320987654321 123456790123456790123456790120987654320987654320987654321 12345679012345679012345679012320987654320987654320987654321 1234567901234567901234567901234320987654320987654320987654321 123456790123456790123456790123454320987654320987654320987654321 12345679012345679012345679012345654320987654320987654320987654321 1234567901234567901234567901234567654320987654320987654320987654321 123456790123456790123456790123456787654320987654320987654320987654321 12345679012345679012345679012345678987654320987654320987654320987654321 1234567901234567901234567901234567900987654320987654320987654320987654321 123456790123456790123456790123456790120987654320987654320987654320987654321 12345679012345679012345679012345679012320987654320987654320987654320987654321 1234567901234567901234567901234567901234320987654320987654320987654320987654321 123456790123456790123456790123456790123454320987654320987654320987654320987654321 12345679012345679012345679012345679012345654320987654320987654320987654320987654321 1234567901234567901234567901234567901234567654320987654320987654320987654320987654321 123456790123456790123456790123456790123456787654320987654320987654320987654320987654321 12345679012345679012345679012345679012345678987654320987654320987654320987654320987654321 1234567901234567901234567901234567901234567900987654320987654320987654320987654320987654321 123456790123456790123456790123456790123456790120987654320987654320987654320987654320987654321 12345679012345679012345679012345679012345679012320987654320987654320987654320987654320987654321 1234567901234567901234567901234567901234567901234320987654320987654320987654320987654320987654321 123456790123456790123456790123456790123456790123454320987654320987654320987654320987654320987654321 

Beautiful, isn't it?

SPOJ POLCONST

Problem link: SPOJ Problem Set (classical): 17707. Constructible Regular Polygons

As stated in the problem description, you are required to utilize this page. As there are only 5 Fermat Primes, you can easily do a pre-processing up to 10⁶ in constant time (2⁵×5) and answer each query in O(1).

SPOJ FINDPRM

SPOJ Problem Set (Classical) 5970. Finding Primes

Fun problem! First it wants you to run a sieve from 2 to N, and then do the similar operations starting from N to 2, except, this time, you have to mark the factors of N and then find the next largest N1 not yet marked, then mark factors of N1 and so on... Finally, you have to tell how many numbers will be unmarked by both algorithms up to N.

First observation is, you only need to consider those numbers which are primes, so, we can run a sieve up to 10⁷, which will allow us to test whether a number is prime or not in O(1) time.

Now, for each N, we can pre-calculate the result. Just note that, say if you know ans[n-1] and want to determine ans[n] from it, this inequality will hold |ans[n] - ans[n-1]| <= 1. Because the current N will be added as an unmarked, or there may be a number which you added in the past, will be removed by one of the factors of this N. Or, if N is not a prime, then it will be already removed by the first algorithm i.e. sieve.

If N is a prime number, then no N1 has removed this N already, so in this case the ans[n] = ans[n-1] + 1.

If N is an even number, and if N/2 was a prime, you have added that with your result, but this should be removed at this stage. So in such a case, ans[n] = ans[n-1].

Why don't you need to consider other factors like N/3, N/5 ... ? 2 is the only even prime factor and it can produce even numbers by multiplying other primes with it. For other primes, p = 3, 5 and so on, if N/p is prime, then N cannot be even. So the basic idea is as follows:

 for( i = 2; i <= N; i++ ) {     ans[i] = ans[i - 1];     if( i is even and i/2 is prime ) ans[i] = ans[i] - 1;     if( i is prime ) ans[i] = ans[i] + 1; } 

I had absolutely no clue at the first glance. Really nice problem.

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.

Pollard's Rho in Java

This is a Pollard's Rho implementation in java. Not very fast, but works for uva online judge. The reason behind using java is the default support of the BigInteger class.

import java.math.BigInteger;
import java.security.SecureRandom;
import java.io.*;
import java.util.*;

public class PollardRho {

    private final static BigInteger ZERO = new BigInteger("0");
    private final static BigInteger ONE  = new BigInteger("1");
    private final static BigInteger TWO  = new BigInteger("2");
    private final static SecureRandom random = new SecureRandom();

    static Vector<BigInteger> v = new Vector<BigInteger>();

    public static BigInteger rho(BigInteger N) {

        BigInteger divisor;
        BigInteger c  = new BigInteger(N.bitLength(), random);
        BigInteger x  = new BigInteger(N.bitLength(), random);
        BigInteger xx = x;

        if (N.mod(TWO).compareTo(ZERO) == 0) return TWO;

        do {
            x  =  x.multiply(x).mod(N).add(c).mod(N);
            xx = xx.multiply(xx).mod(N).add(c).mod(N);
            xx = xx.multiply(xx).mod(N).add(c).mod(N);
            divisor = x.subtract(xx).gcd(N);
        } while((divisor.compareTo(ONE)) == 0);

        return divisor;
    }

    public static void factor(BigInteger N) {

        if (N.compareTo(ONE) == 0) return;

        if (N.isProbablePrime(20)) {
            v.add(N);
            return;
        }

        BigInteger divisor = rho(N);
        factor(divisor);
        factor(N.divide(divisor));
    }

    public static void main(String[] args) throws Exception {

        String string = "";
        InputStreamReader input = new InputStreamReader(System.in);
        BufferedReader reader = new BufferedReader(input);

        while(null != (string = reader.readLine())) {
            BigInteger N = new BigInteger(string);
            v.clear();
            factor(N);
            Collections.sort(v);
            for(int i = 0; i < v.size(); i++) System.out.println(v.get(i));
            System.out.println();
        }
    }
}

I have seen this piece of code years ago somewhere in the internet, but can't remember exactly where. So, I will be glad if anyone can comment/mail me the original source. However, for spoj, I need to write a much much better version of this in C++ :( Exam sucks life...

Matrix Exponentiation

Introduction:

Don't be confused with the title, this article has nothing to do with "how to calculate an exponent of a given matrix", rather it will discuss on how to use this technique to solve a specific class of problems.

Sometimes we face some problems, where, we can easily derive a recursive relation (mostly suitable for dynamic programming approach), but the given constraints make us about to cry, there comes the matrix exponentiation idea. The situation can be made more clear with the following example:

Let, a problem says: find f(n) : n'th Fibonacci number. When n is sufficiently small, we can solve the problem with a simple recursion, f(n) = f(n-1) + f(n-2), or, we can follow dynamic programming approach to avoid the calculation of same function over and over again. But, what will you do if the problem says, given 0 < n < 1000000000, find f(n) % 999983 ? No doubt dynamic programming will fail!

We'll develop the idea on how and why these types of problems could be solved by matrix exponentiation later, first lets see how matrix exponentiation can help is to represent recurrence relations.

Prerequisite:

Before continuing, you need to know:

• Given two matrices, how to find their product, or given the product matrix of two matrices, and one of them, how to find the other matrix.
• Given a matrix of size d x d, how to find its n'th power in O( d³log(n) ).

Patterns:

What we need first, is a recursive relation, and what we want to do, is to find a matrix M which can lead us to the desired state from a set of already known states. Let, we know k states of a given recurrence relation, and want to find the (k+1)th state. Let M be a k x k matrix, and we build a matrix A:[k x 1] matrix from the known states of the recurrence relation, now we want to get a matrix B:[k x 1] which will represent the set of next states, i.e. M x A = B, as shown below:

| f(n) | | f(n+1) | | f(n-1) | | f(n) | M x | f(n-2) | = | f(n-1) | | ...... | | ...... | | f(n-k) | |f(n-k+1)|

So, if we can design M accordingly, job's done!, the matrix will then be used to represent the recurrence relation.

• Type 1:

  Lets start by the simplest one, f(n) = f(n-1) + f(n-2).
  So, f(n+1) = f(n) + f(n-1)
  Let, we know, f(n) and f(n-1); we want to get f(n+1)
  From the above situation, matrix A and B can be formed as shown below:
  

  
  
  
  
  

  Matrix A	Matrix B
  | f(n) | | f(n-1) |	| f(n+1) | | f(n) |

  
  [Note: matrix A will be always designed in such a way that, every state on which f(n+1) depends, will be present]
  So, now, we need to design a 2x2 matrix M such that, it satisfies M x A = B as stated above.
  The first element of B is f(n+1) which is actually f(n) + f(n-1). To get this, from matrix A, we need, 1 f(n) and 1 f(n-1). So, the 1st row of M will be [1 1].
  

  | 1 1 | x | f(n) | = | f(n+1) | | ----- | | f(n-1) | | ------ |

  
  [Note: ----- means, we are not concerned about this value]
  Similarly, 2nd item of B is f(n) which we can get by simply taking 1 f(n) from A. So, the 2nd row of M is [1 0].
  

  | ----- | x | f(n) | = | ------ | | 1 0 | | f(n-1) | | f(n) |

  
  [I hope you know how a matrix multiplication is done and how the values ar assigned!]
  Thus we get the desired 2 x 2 matrix M:
  

  | 1 1 | x | f(n) | = | f(n+1) | | 1 0 | | f(n-1) | | f(n) |

  
  

  If you are confused about how the above matrix is calculated, you might try doing it this way:

  
  

  We know, the multiplication of an n x n matrix M with an n x 1 matrix A will generate an n x 1 matrix B, i.e. M x A = B. The k'th element in the product matrix B is the product of k'th row of the n x n matrix M with the n x 1 matrix A in the left side.
  In the above example, the 1st element in B is f(n+1) = f(n) + f(n-1). So, it's the product of 1st row of matrix M and matrix B. Let, the first row of M is [x y]. So, according to matrix multiplication,
  

  
   x * f(n) + y * f(n-1) = f(n+1) = f(n) + f(n-1)
  ⇒ x = 1, y = 1
  

  
  Thus we can find the first row of matrix M is [1 1]. Similarly, let, the 2nd row of matrix M is [x y], and according to matrix multiplication:
  

  
   x * f(n) + y * f(n-1) = f(n)
  ⇒ x = 1, y = 0
  

  
  Thus we get the second row of M is [1 0].
  

• Type 2:

  Now, we make it a bit complex: find f(n) = a * f(n-1) + b * f(n-2), where a, b are some constants.
  This tells us, f(n+1) = a * f(n) + b * f(n-1).
  By this far, this should be clear that the dimension of the matrices will be equal to the number of dependencies, i.e. in this particular example, again 2. So, for A and B, we can build two matrices of size 2 x 1:
  

  
  
  
  
  

  Matrix A	Matrix B
  | f(n) | | f(n-1) |	| f(n+1) | | f(n) |

  
  Now for f(n+1) = a * f(n) + b * f(n-1), we need [a b] in the first row of objective matrix M instead of [1 1] from the previous example. Because, now we need a of f(n)'s and b of f(n-1)'s.
  

  | a b | x | f(n) | = | f(n+1) | | ----- | | f(n-1) | | ------ |

  
  And, for the 2nd item in B i.e. f(n), we already have that in matrix A, so we just take that, which leads, the 2nd row of the matrix M will be [1 0] as the previous one.
  

  | ----- | x | f(n) | = | ------ | | 1 0 | | f(n-1) | | f(n) |

  
  So, this time we get:
  

  | a b | x | f(n) | = | f(n+1) | | 1 0 | | f(n-1) | | f(n) |

  
  Pretty simple as the previous one...



• Type 3:

  We've already grown much older, now lets face a bit complex relation: find f(n) = a * f(n-1) + c * f(n-3).
  Ooops! a few minutes ago, all we saw were contiguous states, but here, the state f(n-2) is missing. Now? what to do?

  
  

  Actually, this is not a problem anymore, we can convert the relation as follows: f(n) = a * f(n-1) + 0 * f(n-2) + c * f(n-3), deducing f(n+1) = a * f(n) + 0 * f(n-1) + c * f(n-2). Now, we see that, this is actually a form described in Type 2. So, here, the objective matrix M will be 3 x 3, and the elements are:
  

  | a 0 c | | f(n) | | f(n+1) | | 1 0 0 | x | f(n-1) | = | f(n) | | 0 1 0 | | f(n-2) | | f(n-1) |

  
  These are calculated in the same way as Type 2. [Note, if you find it difficult, try on pen and paper!]
  



• Type 4:

  Life is getting complex as hell, and Mr. problem now asks you to find f(n) = f(n-1) + f(n-2) + c where c is any constant.

  
  

  Now, this is a new one and all we have seen in past, after the multiplication, each state in A transforms to its next state in B.
  f(n) = f(n-1) + f(n-2) + c
  f(n+1) = f(n) + f(n-1) + c
  f(n+2) = f(n+1) + f(n) + c
  ................................. so on
  So, normally we can't get it through the previous fashions. But, how about now we add c as a state?
  

  | f(n) | | f(n+1) | M x | f(n-1) | = | f(n) | | c | | c |

  
  Now, its not much hard to design M according to the previous fashion. Here it is done, but don't forget to verify on yours:
  

  | 1 1 1 | | f(n) | | f(n+1) | | 1 0 0 | x | f(n-1) | = | f(n) | | 0 0 1 | | c | | c |




• Type 5:

  Lets put it altogether: find matrix suitable for f(n) = a * f(n-1) + c * f(n-3) + d * f(n-4) + e.
  I would leave it as an exercise to reader. The final matrix is given here, check if it matches with your solution. Also find matrix A and B.
  

  | a 0 c d 1 | | 1 0 0 0 0 | | 0 1 0 0 0 | | 0 0 1 0 0 | | 0 0 0 0 1 |

  
  [Note: you may take a look back to Type 3 and 4]
  
• 
  

  Type 6:

  
  

  
  Sometimes, a recurrence is given like this:
  

  
  f(n) = if n is odd, f(n-1) else, f(n-2)
  In short:
  f(n) = (n&1) * f(n-1) + (!(n&1)) * f(n-2)
  

  
  Here, we can just split the functions in the basis of odd even and keep 2 different matrix for both of them and calculate separately. Actually, there might appear many different patterns, but these are the basic patterns.
  



• Type 7:

  Sometimes we may need to maintain more than one recurrence, where they are interrelated. For example, let a recurrence relation be:
  g(n) = 2g(n-1) + 2g(n-2) + f(n), where, f(n) = 2f(n-1) + 2f(n-2). Here, recurrence g(n) is dependent upon f(n) and the can be calculated in the same matrix but of increased dimensions. Lets design the matrices A, B then we'll try to find matrix M.

  
  
  
  

  Matrix A	Matrix B
  | g(n) | | g(n-1) | | f(n+1) | | f(n) |	| g(n+1) | | g(n) | | f(n+2) | | f(n+1) |

  
  Here, g(n+1) = 2g(n) + 2g(n-1) + f(n+1) and f(n+2) = 2f(n+1) + 2f(n).
  Now, using the above process, we can generate the objective matrix M as follows:
  

  | 2 2 1 0 | | 1 0 0 0 | | 0 0 2 2 | | 0 0 1 0 |

So, these are the basic categories of recurrence relations which are used to be solved by this simple technique.

Analysis:

Now that we have seen how matrix multiplication can be used to maintain recurrence relations, we are back to out first question, how this helps us in solving recurrences on a huge range.

Recall the recurrence f(n) = f(n-1) + f(n-2).
We already know that:

M x | f(n) | = | f(n+1) | | f(n-1) | | f(n) |

............(1)

How about we multiply M multiple times? Like this:

M x M x | f(n) | = | f(n+1) | | f(n-1) | | f(n) |

Replacing from (1):

M x M x | f(n) | = M x | f(n+1) | = | f(n+2) | | f(n-1) | | f(n) | | f(n+1) |

So, we finally get:

M^2 x | f(n) | = | f(n+2) | | f(n-1) | | f(n+1) |

Similarly we can show:

M^3 x | f(n) | = | f(n+3) | | f(n-1) | | f(n+2) |

M^4 x | f(n) | = | f(n+4) | | f(n-1) | | f(n+3) |

............................... ............................... ...............................

M^k x | f(n) | = | f(n+k) | | f(n-1) | |f(n+k-1)|

Thus we can get any state f(n) by simply raising the power of objective matrix M to n-1 in O( d³log(n) ), where d is the dimension of square matrix M. So, even if n = 1000000000, still this can be calculated pretty easily as long as d³ is sufficiently small.

Related problems:

UVa 10229 : Modular Fibonacci
UVa 10870 : Recurrences
UVa 11651 : Krypton Number System
UVa 10754 : Fantastic Sequence
UVa 11551 : Experienced Endeavour

Regards, Zobayer Hasan.