1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#ifdef LOC
#include "debuglib.hpp"
#else
#define _USE_MATH_DEFINES
#include <bits/stdc++.h>
#define deb(...)
#define DBP(...)
#endif
using namespace std;
using   ll         = long long;
using   vi         = vector<int>;
using   pii        = pair<int, int>;
#define pb           push_back
#define mp           make_pair
#define x            first
#define y            second
#define rep(i, b, e) for (int i = (b); i < (e); i++)
#define each(a, x)   for (auto& a : (x))
#define all(x)       (x).begin(), (x).end()
#define sz(x)        int((x).size())

constexpr int MOD = int(1e9) + 7;

ll modPow(ll a, ll e, ll m) {
	assert(e >= 0);
	ll t = 1 % m;
	while (e) {
		if (e % 2) t = t*a % m;
		e /= 2; a = a*a % m;
	}
	return t;
}

struct Zp {
	ll x;
	Zp() : x(0) {}
	Zp(ll a) : x(a%MOD) { if (x < 0) x += MOD; }
	#define OP(c,d) Zp& operator c##=(Zp r) { x = x d; return *this; } Zp operator c(Zp r) const { Zp t = *this; return t c##= r; }
	OP(+, +r.x - MOD*(x+r.x >= MOD));
	OP(-, -r.x + MOD*(x-r.x < 0));
	OP(*, *r.x % MOD);
	OP(/, *r.inv().x % MOD);
	Zp operator-() const { return Zp()-*this; }
	Zp inv() const { return pow(MOD-2); }
	Zp pow(ll e) const{ return modPow(x,e,MOD); }
	void print() { cerr << x; }
};

// (1) sum i=0..n: x^i * y^(n-i)
// (2) sum i=0..n: x^i * y^(n-i) * (i+1)
// (3) x^n
tuple<Zp, Zp, Zp> getFunnySums(Zp x, Zp y, ll n) {
	assert(n >= 0);
	if (n == 0) return {1, 1, 1};
	if (n % 2 == 0) {
		auto [a, b, xn] = getFunnySums(x, y, n-1);
		xn *= x;
		return {a*y + xn, b*y + xn*(n+1), xn};
	}
	auto [a, b, xn] = getFunnySums(x*x, y*y, (n-1)/2);
	return {a*(x+y), b*y*2-a*y + b*x*2, xn*x};
}

int main() {
	cin.sync_with_stdio(0); cin.tie(0);
	int n, k, m;
	cin >> n >> k >> m;

	Zp kInv = Zp(k).inv();

	vector<Zp> pEqual(m+1); // P[sum = i]
	vector<Zp> eEqual(m+1); // E[length | sum = i] * P[sum = i]
	{
		Zp ps = 1, es = 0;
		pEqual[0] = 1;
		rep(i, 1, m) {
			pEqual[i] = ps * kInv;
			eEqual[i] = es * kInv + pEqual[i];
			ps += pEqual[i];
			es += eEqual[i];
			if (i-k >= 0) {
				ps -= pEqual[i-k];
				es -= eEqual[i-k];
			}
		}
	}

	vector<Zp> pNice(m+1); // P[sum >= i and without last < i]
	vector<Zp> eNice(m+1); // E[length | sum >= i and without last < i] * P[...]
	{
		Zp pa = 0, pb = 0, ea = 0, eb = 0;
		rep(i, 0, m) {
			pNice[i] = pa*i + pb;
			eNice[i] = ea*i + eb;
			int t = min(i+k+1, m);
			Zp q = kInv * pEqual[i];
			pa -= q; pb += q*t;
			q = kInv * (eEqual[i] + pEqual[i]);
			ea -= q; eb += q*t;
			if (i-k-1 >= 0) {
				q = kInv * pEqual[i-k-1];
				pa += q; pb -= q*i;
				q = kInv * (eEqual[i-k-1] + pEqual[i-k-1]);
				ea += q; eb -= q*i;
			}
		}
		pNice[0] = 1;
	}

	Zp ans = 1;

	rep(t, max(m-k, 0), m) {
		Zp acc = eEqual[t];
		if (n > 1) {
			auto [a, b, c] = getFunnySums(pNice[t+1], pNice[t], n-2);
			acc *= a*pNice[t] + c*pNice[t+1];
			acc += b * pEqual[t] * eNice[t+1];
			acc += (a*n-b) * pEqual[t] * eNice[t];
		}
		ans += acc * (k-m+t+1) * kInv;
	}

	cout << ans.x << '\n';
}