Problem: There used to be unrest in the Galactic Senate. Several thousand solar systems had declared their intentions to leave the Republic. But fear not! Master Heidi was able to successfully select the Jedi Knights that have restored peace in the galaxy. However, she knows that evil never sleeps and a time may come when she will need to pick another group of Jedi Knights. She wants to be sure she has enough options to do so.

There are n Jedi Knights, each of them with a lightsaber of one of m colors. Given a number k, compute the number of differently colored collections of k lightsabers that some k Jedi Knights might have. Jedi Knights with lightsabers of the same color are indistinguishable (it's not the person, it's the lightsaber color that matters!), and their order does not matter; that is, we consider two collections of Jedi Knights to be different if and only if their vectors of counts of lightsabers of each color (like what you were given in the easy and the medium versions) are different. We count all subsets, not only contiguous subsegments of the input sequence. Output the answer modulo 1009.

Input

The first line of the input contains n (1 ≤ n ≤ 2·105), m (1 ≤ m ≤ n) and k (1 ≤ k ≤ n). The second line contains n integers in the range {1, 2, ..., m} representing colors of the lightsabers of subsequent Jedi Knights.

Output

Output one number: the number of differently colored collections of k lightsabers modulo 1009.

Example

Input

4 3 2
1 2 3 2


Output

4

Solution (in CPP):

Solution:
#include <bits/stdc++.h>
using namespace std;
const int maxn = (1 << 19) | 10, lgn = 20;
const int mod = 1009;
inline int pow_mod(int x, int n) {
  int y = 1;
  while (n) {
    if (n & 1) y = (long long)y * x % mod;
    x = (long long)x * x % mod;
    n >>= 1;
  }
  return y;
}
namespace Poly {
struct comp {
  double r, i;
  comp() {}
  comp(double r, double i) : r(r), i(i) {}
} w[lgn][maxn];
comp operator*(const comp &a, const comp &b) {
  return comp(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r);
}
comp operator+(const comp &a, const comp &b) {
  return comp(a.r + b.r, a.i + b.i);
}
comp operator-(const comp &a, const comp &b) {
  return comp(a.r - b.r, a.i - b.i);
}
comp conj(const comp &a) { return comp(a.r, -a.i); }
void init() {
  static const double pi = acos(-1);
  for (int i = 0; i < lgn; ++i) {
    int n = 1 << i;
    if (n >= maxn) break;
    for (int j = 0; j < n; ++j) {
      if (i && (~j & 1))
        w[i][j] = w[i - 1][j >> 1];
      else {
        double ang = 2. * pi * j / (double)n;
        w[i][j] = comp(cos(ang), sin(ang));
      }
    }
  }
  return;
}
inline int toint(double x, double invn) { return llround(x * invn) % mod; }
void dft(int n, comp *a) {
  for (int i = 0, j = 0; i < n; ++i) {
    if (i < j) swap(a[i], a[j]);
    for (int k = n >> 1; (j ^= k) < k; k >>= 1)
      ;
  }
  for (int hl = 1, l = 2, wi = 1; l <= n; hl = l, l <<= 1, ++wi)
    for (int i = 0; i < n; i += l) {
      comp *x = a + i, *y = x + hl, *z = w[wi];
      for (int j = 0; j < hl; ++j, ++x, ++y, ++z) {
        comp t = *y * *z;
        *y = *x - t;
        *x = *x + t;
      }
    }
  return;
}
struct poly {
  vector<int> p;
  poly() { p.clear(); }
  poly(const vector<int> &q) : p(q) {}
  poly(int n, int *a) {
    p.resize(n);
    for (int i = 0; i < n; ++i) p[i] = a[i];
  }
  inline int size() const { return p.size(); }
  inline void resize(int n) {
    p.resize(n);
    return;
  }
};
poly operator*(const poly &A, const poly &B) {
  int n = A.size(), m = B.size();
  if (n < 10 || m < 10 || n + m - 1 < 80) {
    int N = n + m - 1;
    static int c[maxn];
    memset(c, 0, sizeof(c[0]) * N);
    for (int i = 0; i < n; ++i)
      for (int x = A.p[i], j = 0; j < m; ++j)
        c[i + j] = ((long long)x * B.p[j] + c[i + j]) % mod;
    return poly(N, c);
  }
  int N = 1;
  while (N < n + m - 1) N <<= 1;
  static comp a[maxn], b[maxn], da[maxn], db[maxn], dc[maxn], dd[maxn];
  for (int i = 0; i < n; ++i) {
    int x = (A.p[i] + mod) % mod;
    a[i] = comp(x & 32767, x >> 15);
  }
  for (int i = n; i < N; ++i) a[i] = comp(0, 0);
  for (int i = 0; i < m; ++i) {
    int x = (B.p[i] + mod) % mod;
    b[i] = comp(x & 32767, x >> 15);
  }
  for (int i = m; i < N; ++i) b[i] = comp(0, 0);
  dft(N, a);
  dft(N, b);
  for (int i = 0; i < N; ++i) {
    int j = N - i & N - 1;
    static comp DA, DB, DC, DD;
    DA = (a[i] + conj(a[j])) * comp(.5, 0);
    DB = (a[i] - conj(a[j])) * comp(0, -.5);
    DC = (b[i] + conj(b[j])) * comp(.5, 0);
    DD = (b[i] - conj(b[j])) * comp(0, -.5);
    da[j] = DA * DC;
    db[j] = DA * DD;
    dc[j] = DB * DC;
    dd[j] = DB * DD;
  }
  for (int i = 0; i < N; ++i)
    a[i] = da[i] + db[i] * comp(0, 1), b[i] = dc[i] + dd[i] * comp(0, 1);
  dft(N, a);
  dft(N, b);
  double iN = 1.0 / N;
  N = n + m - 1;
  static int c[maxn];
  for (int i = 0; i < N; ++i) {
    static int da, db, dc, dd;
    da = toint(a[i].r, iN), db = toint(a[i].i, iN), dc = toint(b[i].r, iN),
    dd = toint(b[i].i, iN);
    c[i] = (((long long)dd << 30) + ((long long)db + dc << 15) + da) % mod;
  }
  return poly(N, c);
}
poly operator*(const int &a, const poly &B) {
  static int b[maxn];
  int n = B.size();
  for (int i = 0; i < n; ++i) b[i] = (long long)a * B.p[i] % mod;
  return poly(n, b);
}
poly operator-(const poly &A, const poly &B) {
  static int a[maxn];
  int n = A.size(), m = B.size(), N = max(n, m);
  for (int i = 0; i < N; ++i)
    a[i] = ((i < n ? A.p[i] : 0) - (i < m ? B.p[i] : 0)) % mod;
  return poly(N, a);
}
poly operator+(const poly &A, const poly &B) {
  static int a[maxn];
  int n = A.size(), m = B.size(), N = max(n, m);
  for (int i = 0; i < N; ++i)
    a[i] = ((i < n ? A.p[i] : 0) + (i < m ? B.p[i] : 0)) % mod;
  return poly(N, a);
}
poly inv(int n, const poly &A) {
  if (n == 1) {
    return poly(vector<int>{pow_mod(A.p[0], mod - 2)});
  }
  static poly B0, B, tA;
  B0 = inv(n + 1 >> 1, A);
  if (A.size() < n)
    tA = A, tA.resize(n);
  else
    tA.p.clear(), tA.p.insert(tA.p.end(), A.p.begin(), A.p.begin() + n);
  B = 2 * B0 - B0 * B0 * tA;
  B.resize(n);
  return B;
}
poly rev(const poly &A) {
  static poly B;
  B = A;
  reverse(B.p.begin(), B.p.end());
  return B;
}
poly operator/(const poly &A, const poly &B) {
  static poly rA, rB, C, D;
  int n = A.size(), m = B.size();
  rA = rev(A);
  rB = rev(B);
  C = inv(n - m + 1, rB);
  D = C * rA;
  D.resize(n - m + 1);
  return rev(D);
}
poly operator%(const poly &A, const poly &B) {
  static poly D, ret;
  D = A / B;
  ret = A - B * D;
  ret.resize(B.size() - 1);
  return ret;
}
poly Rem(const poly &A, const poly &B, const poly &D) {
  static poly ret;
  ret = A - B * D;
  ret.resize(B.size() - 1);
  return ret;
}
poly dao(const poly &A) {
  static int a[maxn];
  int n = A.size() - 1;
  for (int i = 0; i < n; ++i) a[i] = (long long)A.p[i + 1] * (i + 1) % mod;
  return poly(n, a);
}
poly ji(const poly &A) {
  static int a[maxn];
  int n = A.size();
  for (int i = 1; i <= n; ++i)
    a[i] = (long long)A.p[i - 1] * pow_mod(i, mod - 2) % mod;
  a[0] = 0;
  return poly(n + 1, a);
}
poly ln(int n, const poly &A) {
  static poly B, C;
  C = dao(A);
  B = inv(n, A);
  B = B * C;
  B = ji(B);
  B.resize(n);
  return B;
}
poly exp(int n, const poly &A) {
  if (n == 1) {
    return poly(vector<int>{1});
  }
  static poly B0, B, C, tA;
  B0 = exp(n + 1 >> 1, A);
  if (A.size() < n)
    tA = A, tA.resize(n);
  else
    tA.p.clear(), tA.p.insert(tA.p.end(), A.p.begin(), A.p.begin() + n);
  C = ln(n, B0);
  B = B0 * (poly(vector<int>{1}) - C + tA);
  B.resize(n);
  return B;
}
poly pow_mod(int n, const poly &A, int N) {
  static poly B;
  int na = A.size(), i = 0;
  while (i < na && A.p[i] == 0) ++i;
  if (i == na || (long long)N * i >= n) return poly(vector<int>(n, 0));
  if (i) {
    B.p.clear();
    for (int j = i; j < na; ++j) B.p.push_back(A.p[j]);
    static poly C;
    C = pow_mod(n - N * i, B, N);
    B.p = vector<int>(N * i, 0);
    B.p.insert(B.p.end(), C.p.begin(), C.p.end());
    return B;
  }
  if (A.p[0] != 1) {
    int t = ::pow_mod(A.p[0], mod - 2), s = ::pow_mod(A.p[0], N);
    B.resize(na);
    for (int i = 0; i < na; ++i) B.p[i] = (long long)A.p[i] * t % mod;
    B = pow_mod(n, B, N);
    for (int i = 0; i < n; ++i) B.p[i] = (long long)s * B.p[i] % mod;
    return B;
  }
  B = ln(n, A);
  for (int i = 0; i < n; ++i) B.p[i] = (long long)B.p[i] * N % mod;
  return exp(n, B);
}
}  // namespace Poly
using namespace Poly;
int a[maxn], n, m, k;
vector<poly> v;
int main() {
  scanf("%d%d%d", &n, &m, &k);
  init();
  for (int x, i = 1; i <= n; ++i) scanf("%d", &x), ++a[x];
  sort(a + 1, a + m + 1);
  reverse(a + 1, a + m + 1);
  for (int x = 1; x <= m; ++x)
    if (a[x]) {
      poly p;
      for (int i = 0; i <= a[x]; ++i) p.p.push_back(1);
      v.push_back(p);
      while (v.size() >= 2 && v.back().p.size() >= v[v.size() - 2].size()) {
        auto t = v.back();
        v.pop_back();
        v.back() = v.back() * t;
      }
    }
  while (v.size() >= 2) {
    auto t = v.back();
    v.pop_back();
    v.back() = v.back() * t;
  }
  printf("%d\n", (v[0].p[k] + mod) % mod);
  return 0;
}