LOTERY - Editorial

xcwgf666 · October 8, 2015, 3:58am

Problem Link

Practice

Contest

Difficulty

Medium-Hard

Pre-requisites

Math, persistent segment tree

Problem

You are given a number triangle, which is denoted as:

L(n, 1) = n
L(n, k) = (\frac{1}{L(n - 1, k - 1)}-\frac{1}{L(n, k-1)})^{-1}

You are given T on-line queries on finding the value of LCM(L(n, 1), L(n, 2), ..., L(n, k)) modulo 10^9+7.

How to get 30 points

The elements of the triangle, denoted by L(n, k) are actually the denominators of the terms of Leibniz Harmonic Triangle. So, the value of L(n, k) equals to (n \cdot \binom{n-1}{k-1} ). So, the answer for the query denoted by integers n and k will be equal to LCM(L(n, 1), L(n, 2), ..., L(n, k)) = n \cdot LCM(\binom{n-1}{0}, \binom{n-1}{1}, ..., \binom{n-1}{k-1}).

Let’s use the following equality: (n+1) \cdot LCM(\binom{n}{0}, \binom{n}{1}, ..., \binom{n}{k})=LCM(n-k+1, n-k+2, ..., n+1). Using it, it’s easy to see that the answer to the query denoted as (n, k) is LCM(n-k+1, ..., n).

So now we need to come up with a way to calculate the LCM of the consecutive numbers.

For doing it, first of all, we will need a way to factorize integers fast. According to the statement, there won’t be any integer greater than 10^6, so we can use the modified sieve that not only determines whether the number prime or not, but also gives its’ greater prime divisor. For doing that, any time we mark a number as a composite one, we can also remember the prime divisor that was the evidence of its’ compositeness.

For better understanding, have a look at the following pseudocode. Modification is highlighted in bold:

    for i = 2 to n :
        isPrime[i] = true, largestPrimeDivisor[i] = 0
    for i = 2 to n :
        if isPrime[i] :
            <b>largestPrimeDivisor[i] = i</b>
            j = i + i
            while j <= n :
                isPrime[j] = false
                <b>largestPrimeDivisor[j] = i</b>
                j = j + i

Now, when we need to get a prime factorization of an integer, we can simply divide it by its’ largest prime divisor, until we get one and write out these divisors.

The pseudocode that does that:

    primeDivisors = [empty list]
    while n > 1 :
        append largestPrimeDivisor[n] to the end of primeDivisors
        n = n / largestPrimeDivisor[n]

After that, primeDivisors[] will contain the factorization of an integer n. It isn’t hard to make another list that will store each prime divisor exactly once along with its’ power in the decomposition. The runtime of this factorization is O(\log n).

Now let’s apply it to the calculation of the Least Common Multiple of the range of numbers. Consider the prime factorization of LCM(n-k+1, n-k+2, ..., n). Obviously it will contain only prime numbers not exceeding n, and for each of these primes, the power in the LCM will be equal to the maximal its’ power in the integers between n-k+1 and n.

So, the solution will looks as follows:

Calculate the prime factorization for each number between n-k+1 and n.
Using it, find the maximal occurrence power for each prime from 2 to n, inclusively.
The answer will be equal to the product of primes from 2 to n, each in its’ maximal occurrence power.

The factorization of a single integer runs in O(\log n), in a single query we will factorize no more than M integers. Then, we will have to calculate the powers of no more than M integers, each power can be calculated with binary exponentation in O(\log n). So, the final complexity will be O(TM \log M). This solution is enough for getting 30 points by solving subtasks #1 and #2, but it is too slow for solving two last subtasks.

How to get 100 points

Let’s maintain the array D_n[] with the following property: for each k \leq n it holds that LCM(k, k+1, ..., n) = D_n[k] \cdot D_n[k+1] \cdot ... \cdot D_n[n].

Let’s have a few examples:

D_1 = [1]
D_2 = [1, 2]
D_3 = [1, 2, 3]
D_4 = [1, 1, 3, 4]

Suppose that D_{n-1} is constructed correctly. How do we get D_n using it?

For the beginning, let’s take D_n = D_{n-1} with D_n[n] = n. Let’s note that at some k s, the required property won’t be held. More precisely, that will be all k s such that there exists some l \geq k such that GCD(l, n)>1 i.e. there is some integer l greater or equal than k such that it shares some prime divisors with n. In this case, some prime divisors will be counted twice in the LCM, giving the wrong result.

Now, how do we get rid from these extra powers?

For each prime number (say, p) let’s store the stack of the positions of such indexes of D_n[], at which the corresponding element of D_n[] is currently divisible by p. For each position, we can also calculate the power of p with which it will appear - this is basically the maximal number of times that the position number can be divided by p. For each prime number (say, p, again), its’ stack will contain the positions in the decreasing order of the power of p that occurs at this position, having the position with the smallest power of p at the top.

So in order to get rid of counting the same factor twice in the array D_n[], we first decompose n into a product of primes and for each prime (let p be the corresponding prime number) we do the following:

First, we process all the positions from the corresponding stack that have the power of p less or equal to the power of p in n (i.e. the number of times that n can be divided by p without remainder). We divide these positions by p in their respective powers and pop them from stack. For example, if we have n=4, we will firstly have D_4=[1, 2, 3, 4]. But the stack of prime number 2 will have the information that at position 2, D_4[2] is divisible by 2. The power of 2 in 2 is obviously 1 and it is less than the power of 2 in the integer 4. So we divide D_4[2] by the power of 2 in 2 and remove the information that D_4[2] is divisible by 2 (because it is no longer divisible by 2). Finally, we get D_4 = [1, 1, 3, 4].
Then, there can be an entry in the corresponding stack with the power of p greater that the power of p in n. In this case, we divide the value of D_n[] at that position by p in the power it occurs in n. This way, the sought invariant will be held.
Finally, we add a new entry to the corresponding stack, containing the position n and the maximal power p in n.

The divisions can be made with the use of modular inverse, because the modulo 10^9+7 is a prime number.

It is very convenient to use persistent segment tree for calculating D_n[]. The segment tree should be capable of multiplying a single element and finding the product of the segment. Anytime, the won’t be more than O(log M) elements in each stack, and each prime decomposition won’t contain more than O(\log M) distinct prime numbers. So, it will take O(M \log^2 M) to precalculate D_1[], D_2[], …, D_M[].

Using the property of the array D_n[] you can find the LCM of the consecutive range of integers with a single query to the segment tree. So it will take O(\log M) to process a single query. So, the final complexity turns out to be O(T \log M + M \log^2 M).

Setter’s Solution

Can be found here

Tester’s Solution

Can be found here

arpitdec5 · October 12, 2015, 5:29pm

Could not understand the 100 point soln .
Can someone explain it properly or another way of solving thsi problem ?

uhateme · October 12, 2015, 7:40pm

Well, There was another approach I used for solving this problem.

We can easily reduce the above question to following :
For each query of n and k, we need LCM(n−k+1,n−k+2,...,n).

Now, consider an example of LCM(4,5,6,7,8). Here, if we find highest powers of all the primes which divide any of the numbers given and take product, we can get LCM.
That is, for 2, we get 8; for 3, we get 3 ( 6 is divisible by 3 but no number is divisible by 9); for 5, we get 5; for 7 we get 7; for rest of the primes, we get 1.
Thus, LCM(4,5,6,7,8) = 8*3*5*7 = 840.
So, we build our solution on the similar basis.

Say we need to find LCM(A,A+1,A+2,....,B-1,B)
BRIEF SOLUTION:
 1. For all the primes <= sqrt(b) : we find highest power of that prime that divides a number which lies in [A,B].
 2. We find all the primes which divide a number in range [A,B]

EXPLANATION:

PRECALCULATIONS:
1. We can get all primes upto N in O(N) (Sieve of Eratosthenes).
2. We can also calculate PRIMEPRODUCT[N], where PRIMEPRODUCT[i] will give product of all the prime numbers <= i (mod 10^9 + 7) ( done in O(N) )
3. We can calculate Multiplicative Modular Inverse of all PRIMEPRODUCT[n] and store in MMI[n] ( done in O(No. of Primes upto N * log(10^9 + 7))

Now for each query (A,B), we do the following steps : 
 1. if( A == B)     return A;
 2. For each prime <= sqrt(B) : Get Highest power of prime that divides a number in [A,B]. We do this via following snippet :

   num = primes[i];
   highestpower = 1;
   while (num <= b) {
	     diva = ceil(a/ num);
	     divb = floor(b/num);
	     if(diva<=divb){
			highestpower = num;
             num *= primes[i];
             }else{
	             break;
             }
    }

	Correctness Proof : 

	say num * x = X;
	for X to lie in range [A,B], A<=X<=B
	=> A <= num * x <= B
	=> A/num <= x <= B/num<
	=> There should exist an integer x, s.t. A/num <=x<=B/num.
	   We are only concerened with integral values >= A/num and <=B/num
	=> ceil(A/num) <=x <= floor(B/num)
	   Thus, if diva <= divb; there exists a multiple of num in range [A,B].

 3. Now we are done with all primes <= sqrt(b). Thus, for rest of the primes, we know their hightest powers can be the numbers themselves (For query [50,100] , all the primes greater than sqrt(100) can't have a square of prime as a divisor of any number as they will exceed b.
 We can Calculate this by finding all primes in ranges [A/1 , B/1] , [A/2,B/2], .... [A/sqrt(b) , B/sqrt(b)] and avoiding any repitition. But this will give TLE as for this, the time complexity for each query can go upto sqrt(N). So, we use some little trick. We keep a two checks, lowcheck and highcheck. lowcheck tells that we have calculated highest powers for all primes <= lowcheck and highcheck tells that we have calculated highest powers for all primes >=highcheck.
 We initialize highcheck by infinity and lowcheck by sqrt(b)
 Now , study the following snippet : 

 int divisor = 1;
 while(highcheck > lowcheck ){
	high = floor(b/divisor);
	low = ceil(a/divisor);
	if (high >= highcheck) {
		high = highcheck-1;
	}
	if (low <= lowcheck) {
		low = lowcheck+1;
		highcheck = -INFINITY;
	}
	product = (PRIMEPRODUCT[high] * MMI[low-1]) mod 10^9 + 7;
        ans = (ans * product) mod 10^9 + 7
	highcheck = min(highcheck, low);
	divisor = max(divisor+1 , floor(b/highcheck));
 }

 Here, low and high gives the range for which we haven't calculated the primes till now. We get product for this by getting product for all primes <= high and dividing it by product of all primes < low.
 Initially we take divisor =1. We update highcheck as min(highcheck, low) because for each range [A/1 , B/1] , [A/2,B/2], .... [A/sqrt(b) , B/sqrt(b)], high of former range is less than high of later range. So, we are sure that for all primes p s.t. high < p < highcheck, contribution will be of none. we calculate from low to high. thus, we have checked all the primes >= low. 

 Now, we only need to see for high < highcheck as for others, we have already calculated. high = floor(b/divisor). Thus, divisor = max(divisor+1 , b/highcheck); (divisor + 1, to make sure it does not check for same divisor again and again).

 Hence, We can get the answer.

COMPLEXITY:
	PRECALCULATIONS: O(N) + O(N) + O(1000*log(10^9 + 7)) [1000 is the number of primes <= 10^5].
	FOR EACH QUERY : O(65*log(N)) [65 is number of primes upto sqrt(N)] + O(sqrt(N)).
THUS, overall complexity = O(N + 1000*log(10^9 + 7) + 65*Q*log(N) + Q*sqrt(N))

SOLUTION LINK

ceilks · October 12, 2015, 11:45pm

I used a similar computation as @uhateme but somewhat simpler in the second part. I used a static decomposition in primes lower than sqrt(N) and higher than sqrt(N). For the “small” primes i used the same method as @uhatme. for “larger” primes I can be shure that each prime occurs in each interval queried at most once.

So I precomputed a prefix-product array § for an array with the value of the prime at the prime position an d 1 otherwise. The product of large primes in an interval [a,b] can then be easily computed by

p[b]*p[a-1]^-1 mod 1000000007

gnsp · October 13, 2015, 1:04am

This is the way I used for factorizing integers from 1 to m. This does the job with almost no extra cost in terms of time than the Sieve of Eratosthenes. It exploits the fact that the largest exponent of a prime P in an integer N is equal to ( 1 + largest exponent of the prime in N/P ).

  0 ## FAST PRIME FACTORIZATION OF ALL INTEGERS UPTO N
  1 def initSeive(m):
  2     global sv
  3     sv = [[True, {}, i] for i in range(m+1)]   # sv[i][0] = is prime (boolean), 
  4                                                # sv[i][1] = factorization representation (dict) 
  5                                                # sv[i][2] = residue of division by primes
  4     sv[1] = [False, {}, 1]
  5     sv[2] = [True, {2:1}, 1]
  6     lim = int(m**0.5)+2
  7     for i in range(2, lim):
  8         if sv[i][0]:
  9             sv[i][1][i]=1
 10             sv[i][2] = 1
 11             for j in range(2*i, m+1, i):
 12                 sv[j][0] = False
 13                 quotient = sv[j][2]//i
 14                 sv[j][1][i] = 1 + (sv[quotient][1][i] if i in sv[quotient][1] else 0)
 15                 sv[j][2] //= i**(sv[j][1][i])
 16 
 17     for i in range(2,m+1):
 18         if sv[i][2] != 1:
 19             sv[i][1][sv[i][2]] = 1
 20             sv[i][2] = 1

rushilpaul · October 13, 2015, 6:45pm

Can someone explain why this would be true?

(n+1)⋅LCM((^n_0),(^n_1),...,(^n_k))=LCM(n−k+1,n−k+2,...,n+1)

footman · October 14, 2015, 3:54pm

How many elements does each D_i[] array have?
I think we also need to keep track of the updated segment trees values but what would be the time complexity? Is it included in the O(Mlog^2 M) of precalculation?

trickytechies · October 15, 2015, 8:14am

How can someone come up with the fact that the F(N,K) of the problem was related to Leibniz Harmonic Triangle?

kaizer · October 15, 2015, 10:11pm

Hi all, I have just a little a supplement.

Sorry or the delay, the proofs of the main equations is provided here
Also, notice that it was not intended you to know about Leibniz Harmonic Triangle(but it’s good if you know). The relation between given recurrent formula and binomial coefficients can be found if you just write down some initial values of function, and guess about the pattern.
Soon after contest started I noticed that time complexity can be also improved. I’d like to share this details with you(here). And now total time complexity is O(MlogM + TlogM).
Also, would like to say, that I didn’t want O(X*sqrt(X)) solution to be accepted(and my such a solution didn’t pass), but it appeared a lot of solutions(seems like at least more that half of all AC > 30) with such a complexity which passed(some of them are very fast, actually).
And the last one, for those how didn’t understand explanation below, let me give you some clarifications(and example of calculating D[]):

D_i[] contains exactly i elements.

stack[prime] contains all indexes of elements of D[] which are divisible by prime in decreasing order.

for each n, i, D_n[i] <= i

Consider array D_11[]: D = {1,1,1,1,1,1,7,4,9,10,11} // 1-based

stack2 = {10, 8}
stack[3] = {9}, stack[5] = {10}, stack[7] = {7}, stack[11] = {11}

Adding 12:
```
  D[12] = 12

  p = 4:

  p /= 2, D[stack[2 ].top] /= 2 (and then remove from stack, because it's equal one)

  p /= 2, D[stack[2 ].top] = D[8] /= 2 -> D[8] = 2

  p = 3:

  p /= 3, D[9] /= 3 -> D[9] = 3
```
Array D_12[] = {1,1,1,1,1,1,7,2,3,5,11,12}

So LCM(8,9,10,11,12) = 2 * 3 * 5 * 11 * 12

arunnsit · December 5, 2015, 4:17pm

@ceilks can you a suggest a way to calculate p[b]*p[a-1]^-1 mod 1000000007 !

vastolorde95 · April 10, 2017, 5:22pm

I think the correct time complexity for the author’s solution is \mathcal{O}\left(N \log N \log M \sum_{p}\frac{1}{p\log p}\right) where p are primes are less than M. But since \sum_{p}\frac{1}{p\log p} \approx 1.55 for M = 10^5, we can say that the complexity is equivalent to \mathcal{O}(N\log N \log M)

rajat1603 · October 12, 2015, 9:34pm

basically the idea is that for all primes <= sqrt(100000) , we can bruteforce them and for primes larger than that , their highest power will be just 1 , so we have to answer a query of the form (product of all distinct primes from l to r) where l = n - k + 1 and r = n , we this we can use a persistent segment tree , for more details about queries regarding distinct numbers , refer to august long challenge → simple queries editorial . here is my implementation CodeChef: Practical coding for everyone

greendragons · October 13, 2015, 12:22am

Even i used kind of similar approach, for below sqrt(n) precomputed and for higer did in a loop.

greendragons · October 13, 2015, 10:23pm

To understand it better, lets move factor or (n+1) inside, so we get
LCM(C(n,0)(n+1), C(n,1)(n+1), C(n,2)(n+1),…,C(n,k)(n+1))
Now expand the terms according to formula C(n,r) = [n*(n-1)*(n-2)…(n-r+1)]/r!.
we get => [(n+1), (n+1)*n, (n+1)n(n-1)/2,…] Observe that for some denom 1.2.3…d, we have terms (n+1)n(n-1)…(n-d+1), i.e d+1 consecutive terms, which means that if we only consider numerator from (n+1).n.(n-1)…(n-d), there will be some number divisible by d, which means we don’t need to disturb the new term i.e (n-d+1), so these new terms (n+1),n,(n-1)… contribute in LCM.

rushilpaul · October 14, 2015, 12:46am

they contribute in LCM on the RHS, but what about the denominators on the LHS? I don’t understand your proof. And after googling a bit, I found this link: number theory - Is there a direct proof of this lcm identity? - Mathematics Stack Exchange

greendragons · October 14, 2015, 1:51am

What i meant is all terms in denom will cancel out by numerator, for that it is sufficient to use only first d terms, no the newly added term need not take part in division, means we can use it for lcm as it differs from all previous added numbers. It can also be shown using Proof by Induction.