# COUNTIT - EDITORIAL

**URL:** <https://discuss.codechef.com/t/countit-editorial/27153>\
**Category:** editorial\
**Tags:** editorial, junechallenge, june19\
**Created:** [June 7, 2019, 8:05pm UTC](https://discuss.codechef.com/t/countit-editorial/27153 "2019-06-07T20:05:52Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![teja349](https://discuss.codechef.com/user_avatar/discuss.codechef.com/teja349/32/3653_2.png) [@teja349](https://discuss.codechef.com/u/teja349)\
**Post date:** [June 7, 2019, 8:05pm UTC](https://discuss.codechef.com/t/countit-editorial/27153/1 "2019-06-07T20:05:52Z")

</div>

# PROBLEM LINK:

[Practice](https://www.codechef.com/problems/COUNTIT)

[Contest: Division 1](https://www.codechef.com/JUNE19A/problems/COUNTIT)

[Contest: Division 2](https://www.codechef.com/JUNE19B/problems/COUNTIT)

**Setter:** [Andrey Filimonov](https://www.codechef.com/users/filyan)

**Tester:** [Radoslav Dimitrov](https://www.codechef.com/users/radoslav192)

**Editorialist:** [Teja Vardhan Reddy](https://www.codechef.com/users/teja349)

# DIFFICULTY:

Medium - Hard

# PREREQUISITES:

Lagrange Interpolation, combinatrics, counting

# PROBLEM:

Consider all matrices with N rows (numbered 1 through N) and M columns (numbered 1 through M) containing only integers between 1 and K (inclusive). For each such matrix A, let’s form a sequence L\_1, L\_2, \ldots, L\_{N+M}:

For each i (1 \le i \le N), L\_i is the maximum of all elements in the i-th row of A. Lets call this N elements as sequence B.  
For each i (1 \le i \le M), L\_{N+i} is the maximum of all elements in the i-th column of A.Lets call these M elements as sequence C.

Find the number of different sequences $$ formed this way.

# EXPLANATION

Let us try to analyse a bit what kind of sequences can be obtained.  
**Claim 1:** Maximum in B and C must be equal.  
**Proof:** Proof by contradiction.  
Lets say B\_{max} \gt C\_{max}  
Let B\_{max} occur at (i,j), then the C\_j \geq B\_{max}. And since C\_{max}\>=C\_{j} which implies C\_{max}\>=B\_{max} which is a contradiction.

So now maximum in B and C is equal.

**Claim:** Any array B and C with same maximum can be obtained.  
**Proof:** Lets say maximum in B and C is x. Now lets say x occurs in B at position p and x occurs in C at position q. We will provide a construction such that rest can be anything from 1 to x in both B and C.

A[p][i]=C\_i for each i (1 \le i \le M)  
A[i][q]=B\_i for each i (1 \le i \le N)  
Fill rest of the values as 1.

Now lets see how ith row looks

i \ne p, 1,1,1,...q-1 times , B\_i,1,1,1....,1. Hence, maximum will be B\_i.

i = p, C\_1,C\_2,.....C\_{q-1},C\_q,C\_{q+1},.....C\_m. Hence , maximum will be C\_{max} = C\_q = B\_p.

Similarly we can see along columns that array C is also satisfied.

Now to count the number of pairs of arrays B and C. We will show how to count for a fixed maximum x and later we can take summation across x from 1 to k.

Number of n length arrays with each element belonging to [1,x] = x^n (since each element has x choices).

Now try to figure out why above is not the answer we are searching.

In the above case, we are counting array which have elements belonging to [1,x] but all elements belong to [1,x-1] which makes maximum \< x.

So we need to subtract them to get the required answer = x^n - (x-1)^n. This is a polynomial in x of size n.

Similarly, number of C sequences will be x^m-(x-1)^m. This is a polynomial in x of size m.

Now, we want \sum\_{x=1}^{k} (x^n - (x-1)^n)\* (x^m - (x-1)^m). This is a polynomial in x of n+m.

We can prove that \sum\_{i=1}^k i^n is a polynomial of degree of n+1 in k. Proof attached [here](https://math.stackexchange.com/questions/18983/why-is-sum-limits-k-1n-km-a-polynomial-with-degree-m1-in-n)

Hence, \sum\_{x=1}^{k} (x^n - (x-1)^n)\* (x^m - (x-1)^m) (Lets call this S(k)). S(k) will be a polynomial of degree n+m+1 using above two conclusions.

We can calculate answers for S(1),S(2),.....S(n+m+2) by iterating over k and calculating S(k) from S(k-1).

Now if (k\lt n+m+3), we can answer it here itself. So we assume k \gt n+m+2.  
Now we know n+m+2 distinct value pairs satisfied by S. Since S is a n+m+1 degree polynomial, using these n+m+2 pairs , we can find the exact polynomial S using [Lagrange interpolation](https://en.wikipedia.org/wiki/Lagrange_polynomial).

Lets say we have n+m+2 pairs (x\_1,y\_1),.....(x\_{n+m+2},y\_{n+m+2}).

S = \sum\limits\_{i=1}^{n+m+2} y\_i\* l\_i.  
l\_i = \prod\limits\_{(1\<=j\<=n+m+2,j \ne i)}\frac{x-x\_j}{x\_i-x\_j}.

Now we need to substitute x=k and x\_i = i and respective values of y\_i which we calculated above.

Let p = \prod\_{i=1}^{n+m+2}(k-i).

Now l\_i = \frac{p/(k-i)}{\prod\limits\_{(1\<=j\<=n+m+2,j \ne i)}{(i-j)}}

Now lets try to simplify, \prod\limits\_{(1\<=j\<=n+m+2,j \ne i)}{(i-j)}   
 = (i-1)\*(i-2)\*...2\*1\*(-1)\*(-2)\*(i-(n+m+2))  
= (i-1)! \* (n+m+2-i)! \*(-1)^{n+m+2-i}.

So we precompute p and factorial values, then we can find each l\_i in O(log(MOD)) time because we will need to find a few inverse modulos.

And substituting l\_i in S gives the required S(k).

And thats a wrap!!

# TIME COMPLEXITY

Finding first n+m+2 values takes O((n+m)\*log(n+m)) because we need fast exponentation for each term.

Finding p takes O(n) time.

precomputing factorials take O(n) time.

Computing each l\_i takes O(log(MOD)) time as mentioned above.

Finally combining all and getting S(k) takes O(n+m) time.

So total complexity = O((n+m)\*log(MOD)).

# SOLUTIONS:

> **Setter's Solution**
>
> ```
> #define y1 askjdkasldjlkasd
> #include <bits/stdc++.h>
> #undef y1
>  
> using namespace std;
>  
> #define pb push_back
> #define mp make_pair
> #define fi(a, b) for(int i=a; i<=b; i++)
> #define fj(a, b) for(int j=a; j<=b; j++)
> #define fo(a, b) for(int o=a; o<=b; o++)
> #define fdi(a, b) for(int i=a; i>=b; i--)
> #define fdj(a, b) for(int j=a; j>=b; j--)
> #define fdo(a, b) for(int o=a; o>=b; o--)
> #define sz(x) (int)x.size()
>  
> typedef long long ll;
> typedef long double ld;
> typedef vector<int> vi;
> typedef pair<int, int> pii;
> typedef vector<pii> vpii;
> typedef vector<ll> vl;
> typedef pair<ll, ll> pll;
> typedef vector<pll> vpll;
> typedef vector<ll> vll;
>  
> #ifdef LOCAL
> #define err(...) fprintf(stderr, __VA_ARGS__ )
> #else
> #define err(...) while (0)
> #endif
>  
> double START_TIME;
>  
> void exit() {
> #ifdef LOCAL    
> cerr << "TIME: " << setprecision(5) << fixed << (clock() - START_TIME) / CLOCKS_PER_SEC << endl;
> #endif  
> exit(0);
> }
>  
> template<typename A, typename B>
> ostream& operator<<(ostream& os, pair<A, B> p) {
> os << "(" << p.first << ", " << p.second << ")";
> return os;
> }
>  
> template<typename T>
> ostream& operator<<(ostream& os, vector<T> v) {
> fi(0, sz(v) - 1) {
> os << v[i] << " ";
> }
> return os;
> }
>  
> template<typename T>
> ostream& operator<<(ostream& os, set<T> t) {
> for (auto z : t) {
> os << z << " ";
> }
> return os;
> }
>  
> template<typename T1, typename T2>
> ostream& operator<<(ostream& os, map<T1, T2> t) {
> cerr << endl;
> for (auto z : t) {
> os << "\t" << z.first << " -> " << z.second << endl;
> }
> return os;
> }
>  
> #ifdef LOCAL
> #define dbg(x) {cerr << __LINE__ << "\t" << #x << ": " << x << endl;}
> #define dbg0(x, n) {cerr << __LINE__ << "\t" << #x << ": "; for (int ABC = 0; ABC < n; ABC++) cerr << x[ABC] << ' '; cerr << endl;}
> #else
> #define dbg(x) while(0){}
> #define dbg0(x, n) while(0){}
> #endif
>  
> #ifdef LOCAL
> #define ass(x) if (!(x)) { cerr << __LINE__ << "\tassertion failed: " << #x << endl, abort(); }
> #else
> #define ass(x) assert(x)
> #endif
>  
> ///////////////////////////////////////////////////
>  
> const int MOD = 1e9 + 7;
> const int MAX = 2e5 + 41;
>  
> int add(int a, int b, int MOD) {
> a += b;
> if (a >= MOD) a -= MOD;
> return a;
> }
>  
> int sub(int a, int b, int MOD) {
> a -= b;
> if (a < 0) a += MOD;
> return a;
> }
>  
> int mul(int a, int b, int MOD) {
> return (ll) a * b % MOD;
> }
>  
> int bp(int x, int d, int MOD) {
> int res = 1;
> while (d) {
> if (d & 1) res = mul(res, x, MOD);
> x = mul(x, x, MOD);
> d >>= 1;
> }
> return res;
> }
>  
> int inv(int x, int MOD) {
> return bp(x, MOD - 2, MOD);
> }
>  
> namespace FFT {
>  
> const int MAGIC = 200;
>  
> const int ROOT_PW = (1 << 20);
>  
> int MODS[3] = {985661441, 976224257, 975175681};
> int ROOTS[3] = {717, 315, 1335};
> int IROOTS[3] = {92105044, 951431260, 590949158}; 
>  
> void fft (vi & a, int root, int iroot, int root_pw, int mod, bool invert) {
> int n = (int) a.size();
>  
> for (int i=1, j=0; i<n; ++i) {
> int bit = n >> 1;
> for (; j>=bit; bit>>=1)
> j -= bit;
> j += bit;
> if (i < j)
> swap (a[i], a[j]);
> }
>  
> for (int len=2; len<=n; len<<=1) {
> int wlen = invert ? iroot : root;
> for (int i=len; i<root_pw; i<<=1)
> wlen = int (wlen * 1ll * wlen % mod);
> for (int i=0; i<n; i+=len) {
> int w = 1;
> for (int j=0; j<len/2; ++j) {
> int u = a[i+j], v = int (a[i+j+len/2] * 1ll * w % mod);
> a[i+j] = u+v < mod ? u+v : u+v-mod;
> a[i+j+len/2] = u-v >= 0 ? u-v : u-v+mod;
> w = int (w * 1ll * wlen % mod);
> }
> }
> }
> if (invert) {
> int nrev = inv(n, mod);
> fi(0, n - 1) {
> a[i] = int (a[i] * 1ll * nrev % mod);
> }
> }
> }
>  
> void multiply (const vi & a, const vi & b, vi & res, int root, int iroot, int root_pw, int mod, bool square = false) {
> vi fa (a.begin(), a.end()), fb (b.begin(), b.end());
> size_t n = 1;
> while (n < max (a.size(), b.size())) n <<= 1;
> n <<= 1;
> fa.resize (n), fb.resize (n);
>  
> fi(0, sz(fa) - 1) {
> fa[i] %= mod;
> }   
> fi(0, sz(fb) - 1) {
> fb[i] %= mod;
> }
>  
> fft (fa, root, iroot, root_pw, mod, false);
> if (!square) {  
> fft (fb, root, iroot, root_pw, mod, false);
> } else {
> fb = fa;
> }
> fi(0, (int) n - 1) {
> fa[i] = mul(fa[i], fb[i], mod);
> }
> fft (fa, root, iroot, root_pw, mod, true);
>  
> res.resize (n);
> fi(0, (int) n - 1) {
> res[i] = fa[i];
> }
> }
>     
> const int INVMOD0INMODS1 = inv(MODS[0], MODS[1]);
> const int INV2 = inv(mul(MODS[0], MODS[1], MODS[2]), MODS[2]);
>  
> int crt(const vi &rems, int MOD) {
> int k = (rems[1] - rems[0]) % MODS[1];
> if (k < 0) k += MODS[1];
> k = mul(k, INVMOD0INMODS1, MODS[1]);
> ll x = (ll) k * MODS[0] + rems[0];
>  
> k = (int) ( (rems[2] - x) % MODS[2]);
> if (k < 0) k += MODS[2];
> k = mul(k, INV2, MODS[2]);
>  
> int res = 0;
> res = add(res, (int) (x % MOD), MOD);
>  
> int tmp = mul(MODS[0], MODS[1], MOD);
> tmp = mul(tmp, k, MOD);
> res = add(res, tmp, MOD);
>  
> return res;
> }   
>  
> vi multiplybrute(const vi &a, const vi &b, int MOD) {       
> vi res(sz(a) + sz(b) - 1, 0);
> fi(0, sz(a) - 1) {
> if (a[i]) {
> fj(0, sz(b) - 1) {
> res[i + j] = add(res[i + j], mul(a[i], b[j], MOD), MOD);
> }
> }
> }
> return res;
> }
>  
> vi multiply(const vi &a, const vi &b, int MOD, bool square = false) {           
> if ( (ll) sz(a) * sz(b) <= MAGIC) {
> return multiplybrute(a, b, MOD);
> }
> vector<vi> resp;
> fi(0, 2) {  
> vi c;
> multiply(a, b, c, ROOTS[i], IROOTS[i], ROOT_PW, MODS[i], square);
> resp.pb(c);
> }
> vi res;
> fi(0, sz(resp[0]) - 1) {
> vi rems;
> fo(0, 2) {
> rems.pb(resp[o][i]);                
> }
> int g = crt(rems, MOD);
> res.pb(g);
> }
> return res;
> }
> };
>  
> vi getinversepolynomial(vi vec) {//return inverse polynomial
> vi v = vec;
> vi res(1, inv(vec[0], MOD));
> int n = sz(vec);
> for (int i = 1; true; i++) {
>         
> int len = (1 << i);
>         
> vi a(len, 0);       
> fj(0, len / 2 - 1) {
> a[j] = mul(res[j], 2, MOD);
> }
>         
> vi res2 = FFT::multiply(res, res, MOD, true);
> if (sz(res2) > len) res2.resize(len);
>  
> vi c;
> c.insert(c.begin(), v.begin(), v.begin() + min(len, sz(v)));
> vi b = FFT::multiply(res2, c, MOD);
>  
> res.resize(len);
> fj(0, len - 1) {
> res[j] = sub(a[j], b[j], MOD);
> }
>  
> if (len > n) break;
> }
>  
> return res;
> } 
>  
> int f[MAX];//factorials
> int invf[MAX];//inverse factorials
> int n, m, k;
>  
> int getc(int n, int k) {//return combinations
> if (n == k) return 1;
> if (k > n) return 0;
> return mul(f[n], mul(invf[k], invf[n - k], MOD), MOD);
> }
>  
> void init() {
> f[0] = 1;
> fi(1, MAX - 1) f[i] = mul(f[i - 1], i, MOD);
> fi(0, MAX - 1) invf[i] = inv(f[i], MOD);
> }
>  
> int b[MAX];//b[i] = i-th Bernoulli number of second kind
>  
> void findbernoulli() {//calculates Bernoulli numbers of second kind
> vi ve;
> fi(0, n + m + 1) {
> ve.pb(invf[i + 1]);
> }   
> ve = getinversepolynomial(ve);      
> ve.resize(n + m + 2);   
> fi(0, n + m + 1) {
> b[i] = mul(ve[i], f[i], MOD);       
> }   
> b[1] = inv(2, MOD);         
> }
>  
> int degsum[MAX];//degsum[i] = sum_{j = 1}{k} j ^ i
>  
> void findfauhalber() {//calculates Fauhalbers sums, save them in degsum
> vi a = vi(n + m + 2, 0);
> fi(0, n + m + 1) {
> int v1 = ::b[i];        
> int v2 = invf[i];
> int v3 = inv(bp(k, i, MOD), MOD);       
> int v = mul(v1, mul(v2, v3, MOD), MOD);     
> a[i] = v;
> }       
> vi b = vi(n + m + 2, 0);
> fi(0, n + m + 1) {
> b[i] = invf[i];
> }
> vi c = FFT::multiply(a, b, MOD);
> fi(0, n + m + 1) {
> c[i] = mul(c[i], f[i], MOD);        
> c[i] = mul(c[i], bp(k, i, MOD), MOD);
> c[i] = sub(c[i], ::b[i], MOD);
> c[i] = mul(c[i], inv(i, MOD), MOD);
> }   
> fi(0, n + m) {
> degsum[i] = c[i + 1];
> }   
> }
>  
> void precalc() {
> // n = 100 * 1000;
> // m = 100 * 1000;
> findbernoulli();        
> findfauhalber();            
> }
>  
> void solve() {
> int ans = 0;
> {//sum_{i = 1}{k} i ^ (n + m) + (i - 1) ^ (n + m)
> int ab = mul(2, degsum[n + m], MOD);
> ab = sub(ab, bp(k, n + m, MOD), MOD);
> ans = add(ans, ab, MOD);
> }
> {//sum_{i = 1}{k} (i - 1) ^ n * i ^ m       
> fj(0, m) {
> int sign = (j % 2 == 1 ? MOD - 1 : 1);
> int v1 = getc(m, j);
> int v2 = degsum[n + m - j];
> int v = mul(v1, mul(v2, sign, MOD), MOD);
> ans = sub(ans, v, MOD);
> }
> }
> {//sum_{i = 1}{k} (i - 1) ^ m * i ^ n       
> fj(0, n) {
> int sign = (j % 2 == 1 ? MOD - 1 : 1);
> int v1 = getc(n, j);
> int v2 = degsum[m + n - j];
> int v = mul(v1, mul(v2, sign, MOD), MOD);
> ans = sub(ans, v, MOD);
> }
> }
> printf("%d\n", ans);
> }
>  
> int main() {
> #ifdef LOCAL
> freopen("input.txt", "r", stdin);
> freopen("output.txt", "w", stdout);
> START_TIME = (double)clock();
> #endif
>  
> init();
> int t;
> scanf("%d", &t);    
> while (t--) {
> scanf("%d %d %d", &n, &m, &k);  
> precalc();
> solve();
> }   
>  
> exit();
> return 0;
> }
> 
> ```

[details = “Tester’s Solution”]  
#include \<bits/stdc++.h\>  
#define endl ‘\n’

```
//#pragma GCC optimize ("O3")
//#pragma GCC target ("sse4")
 
#define SZ(x) ((int)x.size())
#define ALL(V) V.begin(), V.end()
#define L_B lower_bound
#define U_B upper_bound
#define pb push_back
 
using namespace std;
template<class T, class T2> inline int chkmax(T &x, const T2 &y) { return x < y ? x = y, 1 : 0; }
template<class T, class T2> inline int chkmin(T &x, const T2 &y) { return x > y ? x = y, 1 : 0; }
const int MAXN = (1 << 20);
const int mod = (int)1e9 + 7; 
 
int pw(int x, int p)
{
    int r = 1;
    while(p)
    {
        if(p & 1) r = r * 1ll * x % mod;
        x = x * 1ll * x % mod;
        p >>= 1;
    }
 
    return r;
}
 
int inv(int x) { return pw(x, mod - 2); } 
 
int n, m, x;
 
void read()
{
    cin >> n >> m >> x;
}
 
int y[MAXN];
int fact[MAXN];
 
void solve()
{
    y[0] = 0;
    for(int i = 1; i <= n + m; i++)
        y[i] = (y[i - 1] + (pw(i, n) - pw(i - 1, n) + mod) * 1ll * (pw(i, m) - pw(i - 1, m) + mod)) % mod;
 
    if(x <= n + m)
    {
        cout << y[x] << endl;
        return;
    }
 
    // F(k) = SUM y[i] * L(i, k)
    // L(i, k) = PRODUCT (k - j) / (i - j)
 
    int product_up = 1;
    for(int i = 0; i <= n + m; i++)
        product_up = product_up * 1ll * (x - i) % mod; 
 
    fact[0] = 1;
    for(int i = 1; i <= n + m; i++)
        fact[i] = fact[i - 1] * 1ll * i % mod;
 
    int ans = 0;
    for(int i = 0; i <= n + m; i++)
    {
        int v = product_up * 1ll * inv(x - i) % mod;
        v = v * 1ll * inv(fact[i]) % mod;
        v = v * 1ll * inv(fact[n + m - i]) % mod;
        if((n + m - i) & 1) 
            v = (mod - v) % mod;
 
        v = v * 1ll * y[i] % mod;
        ans = (ans + v) % mod;
    }
 
    cout << ans << endl;
}
 
int main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
 
    int T;
    cin >> T;
    while(T--)
    {
        read();
        solve();
    }
 
    return 0;
}

```

[/details]

[details = “Editorialist’s Solution”]  
//teja349  
#include \<bits/stdc++.h\>  
#include   
#include   
#include   
#include   
#include   
#include   
#include   
#include   
#include   
#include   
#include   
#include   
#include   
#include \<ext/pb\_ds/assoc\_container.hpp\>  
#include \<ext/pb\_ds/tree\_policy.hpp\>  
//setbase - cout \<\< setbase (16); cout \<\< 100 \<\< endl; Prints 64  
//setfill - cout \<\< setfill (‘x’) \<\< setw (5); cout \<\< 77 \<\< endl; prints xxx77  
//setprecision - cout \<\< setprecision (14) \<\< f \<\< endl; Prints x.xxxx  
[//cout.precision](https://cout.precision)(x) cout\<\<fixed\<\<val; // prints x digits after decimal in val

```
using namespace std;
using namespace __gnu_pbds;
 
#define f(i,a,b) for(i=a;i<b;i++)
#define rep(i,n) f(i,0,n)
#define fd(i,a,b) for(i=a;i>=b;i--)
#define pb push_back
#define mp make_pair
#define vi vector< int >
#define vl vector< ll >
#define ss second
#define ff first
#define ll long long
#define pii pair< int,int >
#define pll pair< ll,ll >
#define sz(a) a.size()
#define inf (1000*1000*1000+5)
#define all(a) a.begin(),a.end()
#define tri pair<int,pii>
#define vii vector<pii>
#define vll vector<pll>
#define viii vector<tri>
#define mod (1000*1000*1000+7)
#define pqueue priority_queue< int >
#define pdqueue priority_queue< int,vi ,greater< int > >
#define flush fflush(stdout) 
#define primeDEN 727999983
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
 
// find_by_order() // order_of_key
typedef tree<
int,
null_type,
less<int>,
rb_tree_tag,
tree_order_statistics_node_update>
ordered_set;
#define int ll
int powi(int a,int b){
    int ans=1;
    while(b){
        if(b%2){
            ans*=a;
            ans%=mod;
        }
        a*=a;
        a%=mod;
        b/=2;
    }
    return ans;
}
int getinv(int i){
    return powi(i,mod-2);
}
int fact[412345],getsum[412345],pn[412345],pm[412345];
 
main(){
    std::ios::sync_with_stdio(false); cin.tie(NULL);
    int t;
    cin>>t;
   
    int i;
    fact[0]=1;
    f(i,1,412345){
        fact[i]=fact[i-1]*i;
        fact[i]%=mod;
    }
    //cout<<"das"<<endl;
    while(t--){
        int n,m;
        cin>>n>>m;
         int k;
        cin>>k;
        rep(i,2*(n+m)+10){
            pn[i]=powi(i,n);
            pm[i]=powi(i,m);
            //cout<<pn[i]<<endl;
        }
        int val;
        getsum[0]=0;    
        f(i,1,2*(n+m)+10){
            val=pn[i]-pn[i-1];
            val*=pm[i]-pm[i-1];
            val%=mod;
            val+=mod;
            val%=mod;
            getsum[i]=getsum[i-1]+val;
            getsum[i]%=mod;
        }
        if(k<2*(n+m)+10){
            cout<<getsum[k]<<endl;
            continue;
        }
        ll gg=1;
        f(i,1,2*(n+m)+10){
            gg*=(k-i);
            gg%=mod;
        }
        ll wow=0;
        f(i,1,2*(n+m)+10){
            val=fact[i-1];
            val*=fact[2*(n+m)+9-i];
            val%=mod;
            if((2*(n+m)+9-i)%2){
                val*=-1;
                val%=mod;
                val+=mod;
                val%=mod;
            }
            val=getinv(val);
            val*=gg;
            val%=mod;
            val*=getinv(k-i);
            val%=mod;
            val*=getsum[i];
            val%=mod;
            wow+=val;
            wow%=mod;
        }
        cout<<wow<<endl;
 
    }
    return 0;   
}

```

[/details]

Feel free to Share your approach, If it differs. Suggestions are always welcomed. 🙂

---

<div class="post-metadata">

**Author:** ![tejas\_919](https://discuss.codechef.com/user_avatar/discuss.codechef.com/tejas_919/32/3569_2.png) [@tejas\_919](https://discuss.codechef.com/u/tejas_919)\
**Post date:** [June 19, 2019, 7:14am UTC](https://discuss.codechef.com/t/countit-editorial/27153/2 "2019-06-19T07:14:29Z")

</div>

For the people who are new to Lagrange interpolation method, [this question](https://codeforces.com/problemset/problem/622/F) can be a good upsolve

---

<div class="post-metadata">

**Author:** ![islingr](https://www.codechef.com/sites/all/themes/abessive/images/user_default_thumb.jpg) [@islingr](https://discuss.codechef.com/u/islingr)\
**Post date:** [June 19, 2019, 12:27pm UTC](https://discuss.codechef.com/t/countit-editorial/27153/3 "2019-06-19T12:27:32Z")

</div>

Here’s another solution that I think _should_ work but I gave up trying to implement it since I figured it out around 1:00 AM the night before and I had never implemented FFT.

After you get the relation \sum\_{x = 1}^{k} (x^m - (x - 1)^m) (x^n - (x - 1)^n) rewrite it as such

\sum\_{x = 1}^{k} (x^m - (x - 1)^m) (x^n - (x - 1)^n) = \sum\_{x = 0}^{k - 1} ((x + 1)^m - x^m) ((x + 1)^n - x^n)\\ = \sum\_{x = 0}^{k - 1} (x + 1)^{m + n} + x^{m + n} - x^m(x + 1)^n - (x+1)^mx^n\\ = 2 S\_{m + n} + k^{m + n} - \sum\_{x = 0}^{k - 1} \left(x^m(x + 1)^n + (x+1)^mx^n \right)\\

where S\_{r} = 1^r + 2^r + \dots + (k - 1)^r.  
Now, looking at x^m(x + 1)^n

\sum\_{x = 0}^{k - 1} x^m(x + 1)^n = \sum\_{x = 0}^{k - 1} x^m\sum\_{j = 0}^n \binom{n}{j} x^j = \sum\_{j = 0}^n \binom{n}{j} \sum\_{x = 0}^{k - 1} x^{m + j} = \sum\_{j = 0}^n \binom{n}{j} S\_{m + j}

Thus, the final answer would be

2 S\_{m + n} + k^{m + n} - \sum\_{j = 0}^n \binom{n}{j} S\_{m + j} - \sum\_{j = 0}^m \binom{m}{j} S\_{n + j}\\

So, once you find S\_r for all r \le m + n, you can simply find the answer in O(m + n) time since the factorials of binomial coefficient can be precomputed.

If we let

f(x) = \sum\_{r = 0}^{\infty} S\_r \frac{x^r}{r!} = \sum\_{r = 0}^{\infty} \sum\_{j = 0}^{k - 1} j^r \frac{x^r}{r!} = \sum\_{j = 0}^{k - 1} \sum\_{r = 0}^{\infty} \frac{(jx)^r}{r!} = \sum\_{j = 0}^{k - 1} e^{jx} = \frac{e^{kx} - e^x}{e^x - 1}.

Thus, if we find the first m + n coefficients in the power series of f(x) we’d be done.

The numerator is easy to deal with since it’s simply

e^{kx} - e^x = \sum\_{j = 0}^{\infty} (k^j - 1) \frac{x^j}{j!}

For the denominator, note that for each power series there exists it’s inverse, ie, for each g(x) = \sum\_{j = 0}^{\infty} a\_j x^j there exists an h(x) = \sum\_{j = 0}^{\infty} b\_jx^j such that h(x) g(x) = 1.  
Therefore, we need to find the inverse power series of

e^x - 1 = \sum\_{j = 1}^{\infty} \frac{x^j}{j!}.

However since we only need to find the first m + n coefficients of f(x) we can truncate the sequence at m + n, ie, find the power series

\frac{1}{\sum\_{j = 1}^{m + n} \frac{x^j}{j!}}.

after which we can simply multiply the two sequences to get our required S.

Finding the inverse is a classic application of Newton-Ralphson’s method, for a polynomial f let g\_j be the polynomial with coefficients such that the first 2^j coefficients for g\_jf is 1. Then we get the relation

h(g) = f - g^{-1} \implies g\_{n + 1} = g\_n - \frac{h(g\_n)}{h'(g\_n)} = g\_n - \frac{f - g\_n^{-1}}{g\_n^{-2}} = g\_n(2 - g\_nf) 

this is precisely the method described in CF user adamant’s article [here](https://codeforces.com/blog/entry/55572). There is also an alternate proof described in [CF 438E’s editorial](https://codeforces.com/blog/entry/12513).

The most important thing that has not been described here is to multiply two polynomials efficiently. This is known as FFT and you can read more about it in [adamant’s article](https://codeforces.com/blog/entry/55572) linked above and the much-revered [CLRS](https://g.co/kgs/N3Bap4).  
Note that we need to multiply the polynomials modulo 10^9+7 for which reading the section “Multiplication of arbitrary modulus” on page 6 on adamant’s article will be helpful.  
Along with this [cp-algorithms.com](https://cp-algorithms.com/algebra/fft.html) (translation of [e-maxx.ru](http://e-maxx.ru)) has an article on FFT.

Also, while searching around I also found [SERSUM - Editorial](https://discuss.codechef.com/t/sersum-editorial-reupload/18702) which was a problem on Long May 2018 Div 1 and the part about computing L in the editorial is exactly the problem about computing S tackled there.  
And [PFRUIT - Editorial](https://discuss.codechef.com/t/pfruit-editorial/19423) too.  
[admant also has a blog post on computing S](https://codeforces.com/blog/entry/60756).

The final time complexity of the solution would be O((m + n) \log(m + n)) however, note that you can simply precompute the coefficients of \frac{e^x}{e^x - 1} till \max(m + n) and then multiply k^j to the j-th coefficient to get the sequence \frac{e^{kx}}{e^x - 1} from here, all you need to do is subtract the two series which would give a time complexity of O(a \log a + (m + n) T) where T is the number of test cases and a = \max(m + n) which is 2\times 10^5 here.

I tried to implement my solution at like 7:00 AM but I gave up when I realised when something was wrong with my attempt to find \frac{x}{e^x - 1} since I started getting waaayyy off answers if I tried to find more than ~30 coefficients of the series. [Here’s my failed try if you want to see.](https://www.codechef.com/viewsolution/24836864)

**tl;dr:** GenFun are cool, kids.

---

<div class="post-metadata">

**Author:** ![tusg25](https://discuss.codechef.com/user_avatar/discuss.codechef.com/tusg25/32/1532_2.png) [@tusg25](https://discuss.codechef.com/u/tusg25)\
**Post date:** [June 19, 2019, 2:09pm UTC](https://discuss.codechef.com/t/countit-editorial/27153/4 "2019-06-19T14:09:58Z")

</div>

> [@islingr](#):
>
> , you can simply find the answer in O(m+n)O(m + n) time since the factorials of binomial coefficient can be p

Although I could see problem was based on Lagrange Interpolation the moment I got it converted to the formula but was a bit confused seeing jtnydv25 solution having of execution time of around 4.5 sec on the problem. It took me days to figure out this approach so I find this better than the original problem.

---

<div class="post-metadata">

**Author:** ![maverick16](https://discuss.codechef.com/user_avatar/discuss.codechef.com/maverick16/32/12193_2.png) [@maverick16](https://discuss.codechef.com/u/maverick16)\
**Post date:** [June 20, 2019, 9:26am UTC](https://discuss.codechef.com/t/countit-editorial/27153/5 "2019-06-20T09:26:52Z")

</div>

What are some good resources out there so that one can improve in solving math based challenges like this?

---

<div class="post-metadata">

**Author:** ![savaliya\_vivek](https://discuss.codechef.com/user_avatar/discuss.codechef.com/savaliya_vivek/32/4976_2.png) [@savaliya\_vivek](https://discuss.codechef.com/u/savaliya_vivek)\
**Post date:** [June 21, 2019, 11:51am UTC](https://discuss.codechef.com/t/countit-editorial/27153/6 "2019-06-21T11:51:02Z")

</div>

can anyone explain…  
how get (x^n -(x-1)^n) \*(x^m -(x-1)^m) this equation??

---

<div class="post-metadata">

**Author:** ![nagpaljatin141](https://discuss.codechef.com/user_avatar/discuss.codechef.com/nagpaljatin141/32/28700_2.png) [@nagpaljatin141](https://discuss.codechef.com/u/nagpaljatin141)\
**Post date:** [June 21, 2019, 2:33pm UTC](https://discuss.codechef.com/t/countit-editorial/27153/7 "2019-06-21T14:33:45Z")

</div>

Since max(Row) = max(Column) (explained above),  
For K = x, x^n is all possible combinations(with max 0 to x-1) &  
For K = x-1, (x-1)^n is all possible combinations(with max 0 to x-2)  
Therefore, For max equals x-1, (x^n-(x-1)^n) is total possible combinations 🙂 @savaliya_vivek

---

<div class="post-metadata">

**Author:** ![rananjay23](https://discuss.codechef.com/user_avatar/discuss.codechef.com/rananjay23/32/28830_2.png) [@rananjay23](https://discuss.codechef.com/u/rananjay23)\
**Post date:** [June 25, 2019, 9:29am UTC](https://discuss.codechef.com/t/countit-editorial/27153/9 "2019-06-25T09:29:39Z")

</div>

> [@teja349](#):
>
> We can calculate answers for S(1),S(2),…S(n+m+2)S(1),S(2),…S(n+m+2) by iterating over kk and calculating S(k)S(k) from S(k−1)S(k-1).

Can some one explain the paragraphs to me after the above line ?

---

<div class="post-metadata">

**Author:** ![anon76330767](https://www.codechef.com/sites/all/themes/abessive/images/user_default_thumb.jpg) [@anon76330767](https://discuss.codechef.com/u/anon76330767)\
**Post date:** [June 28, 2019, 11:20am UTC](https://discuss.codechef.com/t/countit-editorial/27153/10 "2019-06-28T11:20:25Z")

</div>

these too  
i. [TopCoder Statistics - Problem Statement](http://www.topcoder.com/stat?c=problem_statement&pm=10239)  
ii. [TopCoder Statistics - Problem Statement](http://www.topcoder.com/stat?c=problem_statement&pm=8725)