# PFRUIT - Editorial

**URL:** <https://discuss.codechef.com/t/pfruit-editorial/19423>\
**Category:** editorial\
**Tags:** editorial, hard, july18\
**Created:** [July 20, 2018, 7:05am UTC](https://discuss.codechef.com/t/pfruit-editorial/19423 "2018-07-20T07:05:51Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![melfice](https://discuss.codechef.com/user_avatar/discuss.codechef.com/melfice/32/287830_2.png) [@melfice](https://discuss.codechef.com/u/melfice)\
**Post date:** [July 20, 2018, 7:05am UTC](https://discuss.codechef.com/t/pfruit-editorial/19423/1 "2018-07-20T07:05:51Z")

</div>

### PROBLEM LINK:

[Div 1](http://www.codechef.com/JULY18A/problems/PFRUIT)  
[Div 2](http://www.codechef.com/JULY18B/problems/PFRUIT)  
[Practice](http://www.codechef.com/problems/PFRUIT)

**Author:** [Full name](http://www.codechef.com/users/author_nick)  
**Tester:** [Full name](http://www.codechef.com/users/tester_nick)  
**Editorialist:** [Oleksandr Kulkov](http://www.codechef.com/users/melfice)

### DIFFICULTY:

HARD

### PREREQUISITES:

Moment-generating functions, FFT, computing functions on polynomials

### PROBLEM:

There are k chefs and m workers. For each worker we know integers a and b such that worker can pick from a to b fruits. There are n groups of workers, each group having same a and b. Each worker uniformly choose some number from his [a,b] and picks that amount of fruit.

If the total number of fruits picked during the day is S, then goodness of the day is F(S)=S^k. You have to find expected value of F modulo 998244353. Additionally we know that n\cdot k \leq 10^5.

### QUICK EXPLANATION:

Bring all stuff you know about combinatorics and probability theory, add some stuff you don’t know, blend it in mixer until you get the solution.

### EXPLANATION:

p=998244353 is prime number and p-1=2^{23} \cdot 7 \cdot 17, so there is root of unity of order 2^{23} modulo p and we can use this mod for fft which is strong hint here. Let’s try to get some convenient look for the answer.

Let’s say that each group is described by triple \{a\_i,b\_i,c\_i\} where c\_i is number of workers in the group. Let’s find [moment-generating function](https://en.wikipedia.org/wiki/Moment-generating_function) of random variable S:

M\_S(t) = \mathbb E[e^{tS}]=\sum\limits\_{k=0}^\infty \dfrac{t^k\mathbb E(S^k)}{k!}

We will only need to get coefficient near t^k in expansion of M\_S(t) and multiply it by k! to find an answer.

We should note that amounts of fruits picked by each worker are independent random variables, thus moment-generating function of S=\sum\limits\_{i=1}^m s\_i is product of moment-generating functions for each separate worker. Let’s find it for worker with given numbers a and b.

M\_{s\_i}(t)=\mathbb E[e^{ts\_i}]=\dfrac{1}{b-a+1}\sum\limits\_{k=a}^b e^{tk}=\dfrac{1}{b-a+1}\cdot \dfrac{e^{ta}-e^{t(b+1)}}{1-e^t}

Now we can say that for whole moment-generating function holds:

M\_S(t) = \dfrac{1}{\prod\_{i=1}^n(b\_i-a\_i+1)^{c\_i}}\prod\limits\_{i=1}^n\left[\dfrac{e^{ta\_i}-e^{t(b\_i+1)}}{1-e^t}\right]^{c\_i}

Coefficient on the left is constant and can be calculated in advance. As for product on the right, we can efficiently calculate it by using polynomial \log and \exp, that is:

\prod\limits\_{i=1}^n\left[\dfrac{e^{ta\_i}-e^{t(b\_i+1)}}{1-e^t}\right]^{c\_i} = \exp\left(\sum\limits\_{i=1}^n c\_i \left[ta\_i+\log(1-e^{t(b\_i-a\_i+1)})-\log(1-e^t)\right]\right)

For each summand you may explicitly calculate first k coefficients in expansions of 1-e^{t(b\_i-a\_i+1)} and 1-e^t using formula e^x = \sum\limits\_{k=0}^\infty \dfrac{x^k}{k!}. After this you may perform \log and \exp as described in [this](http://codeforces.com/blog/entry/55572) and [this](http://codeforces.com/blog/entry/56422) articles. This will give you O(nk \log nk) solution. There are few things you should admit before using \log straight ahead. Since coefficient near x^0 must be evaluated manually, you would like it to be equal 0 because you work in modulo field. Argument of your \log operator must have coefficient 1 near x^0. To obtain this you should note following:

1. Coefficient near x^0 is initially 0 for both numenator and denominator so you should divide both of them by x.
2. Coefficient near x^0 of numenator after this won’t be equal to 1 so you should divide numenator by this value and multiply final answer by this value raised to power c\_i.

# ALTERNATIVE SOLUTION 1:

There is also simpler O(nk \log m \log k) solution not messing with \log and \exp but only counting fractions in multipliers using formula for polynomial inverses, which can be found in same articles as for \log and \exp and binary exponentiation to power c\_i, though you need good constant to pass with such solution.

# ALTERNATIVE SOLUTION 2:

This is the one I wrote at first. After that I realized that there is simpler approach which is at main part of editorial now.

To get expected value for S^k we should find M\_S^{(k)}(0) since:

M\_S(t) = \sum\limits\_{n=0}^\infty \dfrac{t^n\mathbb E(S^n)}{n!} \implies M\_S^{(k)}(t) = \sum\limits\_{n=0}^{\infty} \dfrac{t^n\mathbb E(S^{k+n})}{n!} \implies M\_S^{(k)}(0) = \mathbb E(S^k)

1. Now let’s denote G\_i(t) = (e^{ta\_i}+\dots+e^{tb\_i})^{c\_i} and introduce auxiliary array (jump to 2. to understand why do we need them):

g\_{ij} = G\_i^{(j)}(0)=\left[\left(e^{ta\_i}+\dots+e^{tb\_i}\right)^{c\_i}\right]^{(j)}(0) = c\_i(c\_i-1)\dots(c\_i-j+1)(a\_i^j+\dots+b\_i^j)^{c\_i-j}

This will help us to calculate derivatives. Main problem in calculating these values is sum on large segment of powers, i.e. a\_i^j+\dots+b\_i^j. We will start with following telescopic identity (see [this](https://en.wikipedia.org/wiki/Faulhaber%27s_formula) for other ways of calculating these sums, like Faulhaber formula):

(n+1)^{k+1}-1 = \sum\limits\_{m=1}^n[(m+1)^{k+1}-m^{k+1}]=\sum\limits\_{p=0}^{k}\dbinom{k+1}{p} (1^p+\dots+n^p)

Let’s denote S\_p=1^p+\dots+n^p and rewrite it in the following more convenient way:

\sum\limits\_{p=0}^{k} \dfrac{S\_p}{p!} \cdot \dfrac{1}{(k-p+1)!}= \dfrac{(n+1)^{k+1}-1}{(k+1)!}=C\_{k}

If we now denote A\_p = \dfrac{S\_p}{p!} and B\_{k-p}=\dfrac{1}{(k-p+1)!} we can say that:

\sum\limits\_{p=0}^k A\_p B\_{k-p} = C\_{k}

If we multiply both parts by x^k we will see that sequence C\_k corresponds to multiplications of A\_k and B\_{k} as polynomials. Thus to recover A\_k we have to find polynomial inverse for B\_k and multiply C\_k by it. Finding inverse for polynomial is routine procedure which can be done in O(n \log n) for first n coefficients using following formula:

Q\_{k+1}=Q\_k(2-PQ\_k)

Here P is polynomial for which we find iverse and Q\_k is polynomial having first \sim 2^k coefficients equal to ones of 1/P(x) series. You may check [this](http://codeforces.com/blog/entry/55572) or [this](http://codeforces.com/blog/entry/56422) resources to know how to calculate inverse series. This allows you to find A\_k in O(n\log n) which in turns provides you with necessary S\_p for all p from 0 to k and allows you to calculate all g\_{ij} in O(nk \log k). Now that we have g\_{ij} what did we calculate them for?  
2. Answer is to find M\_S^{(k)}(0). It will have constant coefficient of \frac{1}{\prod\_{i=1}^n(b\_i-a\_i+1)^{c\_i}} and the rest is k-th derivative of

\prod\limits\_{i=1}^n \left(\sum\limits\_{k=a\_i}^{b\_i}e^{tk}\right)^{c\_i}=\prod\limits\_{i=1}^n G\_i(t)

in point 0. By [general Leibniz rule](https://en.wikipedia.org/wiki/General_Leibniz_rule#More_than_two_factors) it is equal to the following sum:

\sum\limits\_{k\_1+\dots+k\_n=k}\dbinom{k}{k\_1,\dots,k\_n}G\_1^{(k\_1)}(0)\times\dots\times G\_n^{(k\_n)}(0)=

=k!\sum\limits\_{k\_1+\dots+k\_n=k} \dfrac{g\_{1k\_1}}{k\_1!} \times \dots \times \dfrac{g\_{nk\_n}}{k\_n!}

Here after k! goes coefficient near x^k in product of n polynomials:

P\_i(x) = \sum\limits\_{j=0}^k \dfrac{g\_{ij}}{j!} x^j

This product can be calculated in O(nk \log n \log k) with divide and conquer technique. Bringing all mentioned earlier stuff together finally allows us to solve the problem for good. You can also solve problem in O(nk \log nk) if you calculate \log of each single polynomial, then sum them up and then calculate \exp of the result. This is however much more complicated to implement. Please, refer to [this](http://codeforces.com/blog/entry/55572) or [this](http://codeforces.com/blog/entry/56422) articles mentioned before to learn about \log and \exp over polynomials.

# NOTE ON POWER SUMS:

In alternative solution 2 you have to calculate S\_p=1^p+2^p+\dots+n^p. There is more direct way to calculate them if you consider exponential generating function:

S(x) = \sum\limits\_{p=0}^\infty \dfrac{S\_p}{p!}x^p=\sum\limits\_{k=1}^n \sum\limits\_{p=0}^\infty \dfrac{k^px^p}{p!}=\sum\limits\_{k=1}^n e^{kx}=\dfrac{e^x-e^{x(n+1)}}{1-e^x}

The fraction on the right is actually what we got in alternative solution 2 in closed form and you once again will need to reduce numenator and denominator by x and find inverse series of denominator, which still be O(k \log k).

### AUTHOR’S AND TESTER’S SOLUTIONS:

Author’s solution can be found [here](http://www.codechef.com/download/Solutions/JULY18/setter/PFRUIT.cpp).

Tester’s solution can be found [here](http://www.codechef.com/download/Solutions/JULY18/tester/PFRUIT.cpp).

---

<div class="post-metadata">

**Author:** ![ohweonfire](https://discuss.codechef.com/user_avatar/discuss.codechef.com/ohweonfire/32/3653_2.png) [@ohweonfire](https://discuss.codechef.com/u/ohweonfire)\
**Post date:** [July 23, 2018, 7:32am UTC](https://discuss.codechef.com/t/pfruit-editorial/19423/2 "2018-07-23T07:32:59Z")

</div>

Additionally, the first solution is related to [Bernoulli Numbers](https://en.wikipedia.org/wiki/Bernoulli_Numbers) , and the alternative solution is related to [Stirling numbers of the second kind](https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind). If you find the formulas above tough to read or understand, take a look at these numbers. Much of the formulas in the solution is actually equivalent to the identities of them.

Author’s (alternative) solution in a simpler way:

Calculate \sum\binom S0,\sum\binom S1,\dots,\sum\binom Sk instead of \sum S^k (Divided by \prod(b\_i-a\_i+1)^{c\_i} in the end to get expectation). To write S^k in sums of \sum\binom S0,\sum\binom S1,\dots,\sum\binom Sk, we use [Stirling numbers of the second kind](https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind).

Define P(x)=\sum\_{i=0}^kx^i\sum\binom Si, Q(a,b,x)= \sum\_{i=0}^kx^i\sum\_{j=a}^b\binom ji. We can see if we add a new worker (a,b), new P'(x)=P(x)\cdot Q(a,b,x). That means we only need to get \prod Q^{c\_i}(a\_i,b\_i,x). Every single Q can be calculated directly because \sum\_{j=a}^b\binom ji=\binom {b+1}{i+1}-\binom a{i+1}.

Finally, we can use FFT and divide and conquer to get O(nk\log n\log k) solution. As described in the editorial, we can also use \ln and \exp to get faster but harder to implement O(nk\log nk) solution.

Most contestants’ (first) solution in a simpler way:

Define polynomial P\_2(x)=\sum\_{i=0}^k\frac{x^i}{i!}\sum S^i and Q\_2(x)=\sum\_{i=0}^k\frac{x^i}{i!}\sum\_{j=a}^bj^i. To get Q efficiently, we use [Bernoulli Numbers](https://en.wikipedia.org/wiki/Bernoulli_Numbers). Also, we can see if we add a new worker (a,b), new P\_2'(x)=P\_2(x)\cdot Q\_2(a,b,x). And the rest is exactly the same as the first solution.

To calculate Stirling numbers efficiently, we use the explicit formula of the numbers and FFT. See author’s code.

To calculate Bernoulli numbers, we use polynomial inverses, because the exponential generating function of Bernoulli numbers is \frac x{e^x-1}.

Things described above may come differently as in the editorial, but they lead to the same formulas.

---

<div class="post-metadata">

**Author:** ![vijju123](https://discuss.codechef.com/user_avatar/discuss.codechef.com/vijju123/32/195207_2.png) [@vijju123](https://discuss.codechef.com/u/vijju123)\
**Post date:** [July 26, 2018, 2:14am UTC](https://discuss.codechef.com/t/pfruit-editorial/19423/3 "2018-07-26T02:14:43Z")

</div>

```auto
QUICK EXPLANATION:
Bring all stuff you know about combinatorics and probability theory, add some stuff you don't know, blend it in mixer until you get the solution.

```

```auto
*scrolls down to see the mathematical formulas*

```

…  
…

```auto
*dies from Cancermathics*

```