Problem: Find the number of sequences of N non-negative integers A_1, A_2, ..., A_N that satisfy the following conditions:

* L \leq A_1 + A_2 + ... + A_N \leq R
* When the N elements are sorted in non-increasing order, the M-th and (M+1)-th elements are equal.



Since the answer can be enormous, print it modulo 10^9+7.

Constraints

* All values in input are integers.
* 1 \leq M < N \leq 3 \times 10^5
* 1 \leq L \leq R \leq 3 \times 10^5

Input

Input is given from Standard Input in the following format:


N M L R


Output

Print the number of sequences of N non-negative integers, modulo 10^9+7.

Examples

Input

4 2 3 7


Output

105


Input

2 1 4 8


Output

3


Input

141592 6535 89793 238462


Output

933832916

Solution (in CPP):

Solution:
#include <bits/stdc++.h>

using namespace std;


#ifdef zxc

#include <sys/resource.h>
#include "debug.h"
#include "profile.h"
#define debug(...) cerr << "[" << #__VA_ARGS__ << "]:", debug_out(__VA_ARGS__)
#define profile(X) LOG_DURATION(X)
#else
#define debug(...) 42
#define profile(...) 42
#endif

using ll = int64_t;
using ld = long double;
const ld EPS = 1e-9;
//const ll MOD = 924844033;
const ld PI = 3.141592653589793;
const int maxn = 600001;
const ll MOD = 1e9 + 7;
using cd = complex<double>;

void fast_io() {
    ios_base::sync_with_stdio(false);
    cin.tie(0);
}


ll add(ll a, ll b) {
    a += b;
    if (a >= MOD) {
        a -= MOD;
    }
    return a;
}

ll mul(ll a, ll b) {
    return a * b % MOD;
}

ll pw(ll a, ll b) {
    ll res = 1;
    for (; b; b >>= 1) {
        if (b & 1) {
            res = mul(res, a);
        }
        a = mul(a, a);
    }
    return res;
}

ll fac[maxn];
ll inv[maxn];

void precalc(int n) {
    fac[0] = 1;
    for (int i = 1; i < n; ++i) {
        fac[i] = mul(fac[i - 1], i);
    }
    inv[n - 1] = pw(fac[n - 1], MOD - 2);
    for (int i = n - 2; i >= 0; --i) {
        inv[i] = mul(inv[i + 1], i + 1);
    }
}

ll C(ll n, ll k) {
    if (k < 0 || k > n) {
        return 0;
    }
    return fac[n] * inv[n - k] % MOD * inv[k] % MOD;
}

inline ll distribute(ll n, ll sum) {
    return C(sum + n, n);
}

int close_lg(int n) { return n > 1 ? 32 - __builtin_clz(n - 1) : 0; }

vector<cd> make_complex_roots(int n) {
    vector<cd> w(n);
    for (int i = 0; i < n; ++i) {
        ld ang = 2 * PI * i / n;
        w[i] = cd(cos(ang), sin(ang));
    }
    return w;
}

struct Fourier { // substitute T with cd or ll
    Fourier(vector<cd> roots) : LG(close_lg(roots.size())),
                                sz(roots.size()), tail(sz), w(move(roots)) {
        for (int i = 1; i < tail.size(); ++i) {
            tail[i] = tail[i >> 1] >> 1 | ((i & 1) << (LG - 1));
        }
    }

    void fft(vector<cd> &a, bool inv) const { // a should be a power of two
        int n = a.size();
        int k = close_lg(n);
        assert(1 << k == n);
        for (int i = 0; i < n; ++i) {
            int cur_tail = tail[i] >> (LG - k);
            if (cur_tail < i) {
                swap(a[i], a[cur_tail]);
            }
        }
        for (int step = 0; step < k; ++step) { // merge of segments of length (1<<step) each
            int half = 1 << step;
            int len = half << 1;
            int go = 1 << (LG - step - 1);
            for (int block = 0; block < n; block += len) {
                int pw_ind = 0;
                for (int j = block; j < block + half; ++j) {
                    cd fi = a[j];
                    cd se = a[j + half] * w[pw_ind];
                    a[j] = fi + se;
                    a[j + half] = fi - se;
                    pw_ind += go;
                }
            }
        }
        if (inv) {
            reverse(a.begin() + 1, a.end());
            for (auto &el:a) {
                el /= n;
            }
        }
    }

    int LG, sz;
    vector<int> tail;
    vector<cd> w;
};

vector<cd> mult(vector<cd> fi, vector<cd> se, const Fourier &fft) {
    int e = 1 << close_lg(fi.size() + se.size());
    fi.resize(e);
    se.resize(e);
    for(int i=0; i<e; ++i){
        fi[i] = cd{fi[i].real(), se[i].real()};
    }
    fft.fft(fi, 0);
    for (int i = 0; i < e; ++i)
        fi[i] *= fi[i];
    fft.fft(fi, 1);
    for(int i=0; i<e; ++i){
        fi[i] = cd{fi[i].imag()/2,0};
    }
    return fi;
}

vector<ll> mult_mod(vector<ll> fi, vector<ll> se, const Fourier &fft) {
    profile("mult mod");
    vector<cd> fi_arr[2];
    vector<cd> se_arr[2];
    int shifts[2] = {0, 15};
    {
        for (int i = 0; i < 2; ++i) {
            fi_arr[i].resize(fi.size());
            for (int e = 0; e < fi.size(); ++e) {
                fi_arr[i][e] = (fi[e] >> shifts[i]) & ((1 << 15) - 1);
            }
        }
    }
    {
        for (int i = 0; i < 2; ++i) {
            se_arr[i].resize(se.size());
            for (int e = 0; e < se.size(); ++e) {
                se_arr[i][e] = (se[e] >> shifts[i]) & ((1 << 15) - 1);
            }
        }
    }
    vector<ll> res(fi.size() + se.size());
    for (int i = 0; i < 2; ++i) {
        for (int j = 0; j < 2; ++j) {
            auto multed = mult(fi_arr[i], se_arr[j], fft);
            multed.resize(fi.size() + se.size());
            ll coef = 1ll << (15 * (i + j));
            coef %= MOD;
            {
                for (int e = 0; e < res.size(); ++e) {
                    res[e] += mul(coef, ll(round(multed[e].real())) % MOD);
                    res[e] %= MOD;
                }
            }
        }
    }
    return res;
}

vector<ll> generate_vec(int n, int m, int mx) {
    vector<ll> basic_poly(mx + 1);
    for (int i = 0; i + m <= mx && i <= n - m; ++i) {
        basic_poly[i + m] = C(n - m, i);
        if (i & 1) {
            basic_poly[i + m] = (MOD - basic_poly[i + m]) % MOD;
        }
    }
    vector<ll> res(mx + 1);
    for (int k = 1; k <= mx; ++k) {
        for (int e = 1; e <= mx && e * k <= mx; ++e) {
            res[e * k] = add(res[e * k], basic_poly[e]);
        }
    }
    return res;
}

vector<ll> generate_vec2(int n, int mx) {
    vector<ll> res(mx + 1);
    for (int i = 0; i <= mx; ++i) {
        res[i] = C(n + i - 1, n - 1);
    }
    return res;
}

signed main() {
#ifdef zxc
    struct rlimit rl;
    const rlim_t kStackSize = 512L * 1024L * 1024L;
    assert(!getrlimit(RLIMIT_STACK, &rl));
    rl.rlim_cur = kStackSize;
    assert(!setrlimit(RLIMIT_STACK, &rl));
    freopen("../input.txt", "r", stdin);
//    freopen("../kolya_output.txt", "w", stdout);
#else
    //    mt19937 mt(std::chrono::high_resolution_clock::now().time_since_epoch().count());
    //            freopen("", "r", stdin);
    //        freopen("hack.out", "w", stdout);
#endif

    auto solve = [&](int _) {
        int n, m, l, r;
        cin >> n >> m >> l >> r;
        ll all = add(distribute(n, r), MOD - distribute(n, l - 1));
        Fourier fft(make_complex_roots(1 << 20));
        auto fi = generate_vec(n, m, r);
        auto se = generate_vec2(n, r);
        vector<ll> res = mult_mod(fi, se, fft);
        res.resize(r + 1);
        vector<ll> rres(r + 1);
        for (int i = 1; i <= r; ++i) {
            rres[i] = res[i];
            if (i >= m) {
                rres[i] = add(rres[i], MOD - res[i - m]);
            }
        }
        ll ans = accumulate(rres.begin() + l, rres.begin() + r + 1, 0, add);
        cout << add(all, MOD - ans * C(n, m) % MOD);
    };


    precalc(maxn);

    fast_io();
    cout.precision(9);
    cout << fixed;
    auto start = chrono::steady_clock::now();

    for (int i = 1; i <= 1; ++i) solve(i);
    auto end = chrono::steady_clock::now();

    debug("time", chrono::duration_cast<chrono::milliseconds>(end - start).count() / 1000.0);

}