JZPKIL

#include <cstdio>
#include <ctime>
#include <algorithm>

#define int64 long long
#ifdef WIN32
#define fmt64 "%I64d"
#else
#define fmt64 "%lld"
#endif
#define M 3010
#define P 1000000007LL

using namespace std;

double sum;
bool np[M + 5]; int pm[M / 2], ptot, Total;
int64 fac[M + 5], inv[M + 5], inf[M + 5], B[M + 5];
int T, x, y; int64 n;
int64 fa[100]; int id[100], tot;

int64 mul(int64 x, int64 y, int64 z = P) {
    return (x * y - (int64)(x / (long double)z * y + 1e-3) * z + z) % z;
}

int64 fpm(int64 a, int64 b) {
    int64 res = 1; a %= P;
    for (; b; b >>= 1, a = a * a % P) {
        if (b & 1) res = res * a % P;
    }
    return res;
}

int64 fpm(int64 a, int64 b, int64 c) {
    int64 res = 1;
    for (; b; b >>= 1, a = mul(a, a, c)) {
        if (b & 1) res = mul(res, a, c);
    }
    return res;
}

int64 gcd(int64 a, int64 b) {
    for (int64 c; b; c = a, a = b, b = c % b);
    return a;
}

void prime() {
    int i, j;
    for (i = 2; i <= M; ++i) {
        if (!np[i]) {
            pm[ptot++] = i;
            inv[i] = fpm(i, P - 2);
        }
        for (j = 0; j < ptot && i * pm[j] <= M; ++j) {
            np[i * pm[j]] = 1;
            inv[i * pm[j]] = inv[i] * inv[pm[j]] % P;
            if (i % pm[j] == 0) break;
        }
    }
}

void factorial() {
    fac[0] = 1;
    for (int i = 1; i <= M; ++i) {
        fac[i] = fac[i - 1] * i % P;
    }
    inf[M] = fpm(fac[M], P - 2);
    for (int i = M; i; --i) {
        inf[i - 1] = inf[i] * i % P;
    }
}

bool prime(int64 n) {
    for (int i = 0; i < 12; ++i) {
        int64 a = pm[i], p = n - 1;
        p /= p & -p;
        for (a = fpm(a, p, n); p < n; p <<= 1, a = mul(a, a, n)) {
            if (a == 1 || a == -1) break;
        }
        if (p >= n) return 0;
    }
    return 1;
}

void rho(int64 n) {
    int i; int64 x, y, d;
    if (prime(n)) {
        for (i = 0; i < tot; ++i)
            if (fa[i] == n) { ++id[i]; return; }
        fa[tot] = n, id[tot++] = 1;
        return;
    }
    for (; ; ) {
        x = pm[++Total % 13], y = x, d = 1;
        for (i = 1; d == 1; ++i) {
            x = mul(x, x, n) + 13;
            d = gcd(abs(x - y), n);
            if (i == (i & -i)) y = x;
        }
        if (d == n) continue;
        return rho(d), rho(n / d);
    }
}

void divisor(int64 n) {
    int i; tot = 0;
    for (i = 0; i < ptot && pm[i] * pm[i] <= n; ++i)
        if (n % pm[i] == 0) {
            fa[tot] = pm[i], id[tot++] = 1;
            for (n /= pm[i]; n % pm[i] == 0; ++id[tot - 1], n /= pm[i]);
        }
    if (n > 1 && n <= pm[ptot - 1]) {
        fa[tot] = n, id[tot++] = 1;
    } else if (n > 1) rho(n);
}

void bernoulli() {
    int i, j;
    B[0] = 1;
    for (i = 1; i < M; ++i) {
        for (j = 0; j < i; ++j)
            B[i] += fac[i] * inf[j] % P * inf[i + 1 - j] % P * B[j] % P, B[i] %= P;
        B[i] = (P - B[i]) % P;
    }
    B[1] = (P - B[1]) % P;
}

int64 bernoulli(int64 n, int m) {
    int64 res = 0, base = n %= P;
    for (int i = m; ~i; --i) {
        res += fac[m] * inf[i] % P * inf[m + 1 - i] % P * B[i] % P * base % P, res %= P;
        base = base * n % P;
    }
    return res;
}

int64 h2(int64 p, int k, int z) {
    int xx = x;
    if (xx < z) swap(xx, z);
    int64 base = fpm(p, k * z), temp = fpm(p, xx - z);
    int64 res = base;
    for (int i = 0; i < k; ++i) {
        base *= temp, base %= P;
        res += base, res %= P;
    }
    return res;
}

int64 calc() {
    int i, j; int64 res = 0;
    divisor(n);
    for (i = 1; i <= y + 1; ++i) {
        if (!B[y + 1 - i]) continue;
        int64 t = fac[y] * inf[y + 1 - i] % P * inf[i] % P * B[y + 1 - i] % P;
        for (j = 0; j < tot; ++j) {
            t *= (h2(fa[j], id[j], i) - h2(fa[j], id[j] - 1, i) * fpm(fa[j], y) % P + P) % P, t %= P;
        }
        res += t, res %= P;
    }
    res *= fpm(n, y), res %= P;
    return res;
}

int64 mul0(int64 a, int64 b, int64 c) {
    int64 res = 0;
    for (; b; b >>= 1, a <<= 1, a %= c)
        if (b & 1) res += a, res %= c;
    return res;
}

int main() {
    freopen("jzpkil.in", "r", stdin);
    freopen("jzpkil.out", "w", stdout);

    prime();
    factorial();
    bernoulli();
    for (scanf("%d", &T); T--; ) {
        scanf(fmt64"%d%d", &n, &x, &y);
        if (x == y) {
            printf(fmt64"\n", bernoulli(n, x) * fpm(n, x) % P);
        } else {
            printf(fmt64"\n", calc());
        }
    }
}



  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值