#138. 類歐幾裡得算法

以下除法均為向下取整，定義

(

)

∑

(

)

f(a, b, c, n, k_1, k_2) = \sum\limits_{x = 0} ^{n} x ^{k_1} \left(\frac{a \times x + b}{c}\right) ^ {k_2}

f(a,b,c,n,k1,k2)=x=0∑nxk1(ca×x+b)k2。

∑

(

)

(

≤

)

\sum_{x = 0} ^{n} x ^{k_1} \left(\frac{a \times x + b}{c}\right) ^ {k_2}, (1 \le n, a, c \le 10 ^ 9, 0 \le b \le 10 ^ 9, 0 \le k_1 + k_2 \le 10)\\

x=0∑nxk1(ca×x+b)k2,(1≤n,a,c≤109,0≤b≤109,0≤k1+k2≤10)

≥

′

a \ge c, a' = a\ \%\ c, k = \frac{a}{c}

a≥c,a′=a % c,k=ca：

′

−

′

−

\sum_{x = 0} ^{n} x ^{k_1} \left(kx + \frac{a' \times x + b}{c} \right) ^ {k_2}\\ \sum_{x = 0} ^{n} x ^{k_1}\sum_{i = 0} ^{k_2} C_{k_2} ^{i} (kx) ^ i \left(\frac{a' \times x + b}{c} \right) ^{k_2 - i}\\ \sum_{i = 0} ^{k_2} \left(C_{k_2} ^{i} k ^{i}\right) \sum_{x = 0} ^{n} x ^{i + k_1} \left( \frac{a' \times x + b}{c} \right) ^{k_2 - i}\\ f(a, b, c, n, k_1, k_2) = \sum_{i = 0} ^{k_2} C_{k_2} ^{i} k ^{i} f(a', b, c, n, i + k_1, k_2 - i)\\

x=0∑nxk1(kx+ca′×x+b)k2x=0∑nxk1i=0∑k2Ck2i(kx)i(ca′×x+b)k2−ii=0∑k2(Ck2iki)x=0∑nxi+k1(ca′×x+b)k2−if(a,b,c,n,k1,k2)=i=0∑k2Ck2ikif(a′,b,c,n,i+k1,k2−i)

b \ge c, b' = b \ \%\ c, k = \frac{b}{c}

b≥c,b′=b % c,k=cb：

\sum_{x = 0} ^{n} x ^{k_1} \left(\frac{a \times x + b'}{c} + k \right) ^{k_2}\\ \sum_{x = 0} ^{n} x ^{k_1} \sum_{i = 0} ^{k_2} C_{k_2} ^{i} k ^{i} \left( \frac{a \times x + b'}{c} \right) ^{k_2 - i}\\ \sum_{i = 0} ^{k_2} C_{k_2} ^{i} k^{i} \sum_{x = 0} ^{n} x ^{k_1} \left(\frac{a \times x + b'}{c} \right) ^{k_2 - i}\\ f(a, b, c, n, k_1, k_2) = \sum_{i = 0} ^{k_2} C_{k_2} ^{i} k ^{i} f(a, b', c, n, k_1, k_2 - i)\\

x=0∑nxk1(ca×x+b′+k)k2x=0∑nxk1i=0∑k2Ck2iki(ca×x+b′)k2−ii=0∑k2Ck2ikix=0∑nxk1(ca×x+b′)k2−if(a,b,c,n,k1,k2)=i=0∑k2Ck2ikif(a,b′,c,n,k1,k2−i)

a < c , \ b < c, m = \frac{a \times n + b}{c}

a<c, b<c,m=ca×n+b：

−

[

≥

]

−

n ^{k} = \sum_{i = 1} ^{n} i ^{k} - (i - 1) ^{k}\\ \sum_{x = 0} ^{n} x ^{k_1} \left(\frac{a \times x + b}{c}\right) ^ {k_2}\\ \sum_{x = 0} ^{n} x ^{k_1} \sum_{i = 1} ^{m} \left(i ^{k_2} - (i - 1) ^{k_2}\right)[\frac{a \times x + b}{c} \ge i]\\ \sum_{i = 0} ^{m - 1} \left((i + 1) ^{k_2} - i ^{k_2} \right) \sum_{x = 0} ^{n} x^{k_1} [\frac{a \times x + b}{c} \ge i + 1]\\ \sum_{i = 0} ^{m - 1} \left((i + 1) ^{k_2} - i ^{k_2} \right) \sum_{x = 0} ^{n} x^{k_1} - \sum_{i = 0} ^{m - 1} \left((i + 1) ^{k_2} - i ^{k_2} \right) \sum_{x = 0} ^{\frac{c \times i + c - b - 1}{a}} x ^{k_1}\\

nk=i=1∑nik−(i−1)kx=0∑nxk1(ca×x+b)k2x=0∑nxk1i=1∑m(ik2−(i−1)k2)[ca×x+b≥i]i=0∑m−1((i+1)k2−ik2)x=0∑nxk1[ca×x+b≥i+1]i=0∑m−1((i+1)k2−ik2)x=0∑nxk1−i=0∑m−1((i+1)k2−ik2)x=0∑ac×i+c−b−1xk1

前項可以通過插值得到，我們考慮後項求和，

不難得到

−

(i + 1) ^{k_2} - i ^{k_2}

(i+1)k2−ik2，是一個

k_2 - 1

k2−1次的多項式，我們可以通過插值預處理出其系數，設其為

A(x)

A(x)。

對于

\sum\limits_{i = 0} ^{\frac{c \times i + c - b - 1}{a}} x ^{k_1}

i=0∑ac×i+c−b−1xk1，也是一個多項式，最高次幂為

k_1 + 1

k1+1，設其為

B(x)

B(x)。

−

\sum_{i = 0} ^{m - 1} \left((i + 1) ^{k_2} - i ^{k_2} \right) \sum_{x = 0} ^{n} x^{k_1} - \sum_{x = 0} ^{m - 1} \sum_{i = 0} ^{k_2 - 1} A_i x ^{i} \sum_{j = 0} ^{k_1 + 1} B_j \left( \frac{c \times x + c - b - 1}{a} \right) ^{j}\\ \sum_{i = 0} ^{m - 1} \left((i + 1) ^{k_2} - i ^{k_2} \right) \sum_{x = 0} ^{n} x^{k_1} - \sum_{i = 0} ^{k_2 - 1} A_i \sum_{j = 0} ^{k_1 + 1} B_j \sum_{x = 0} ^{m - 1} x ^{i} \left( \frac{c \times x + c - b - 1}{a} \right) ^{j}\\ f(a, b, c, n, k_1, k_2) = \sum_{i = 0} ^{m - 1} \left((i + 1) ^{k_2} - i ^{k_2} \right) \sum_{x = 0} ^{n} x^{k_1} - \sum_{i = 0} ^{k_2 - 1} \sum_{j = 0} ^{k_1 + 1} A_i \times B_j \times f(a', b', c', m - 1, i, j)\\

i=0∑m−1((i+1)k2−ik2)x=0∑nxk1−x=0∑m−1i=0∑k2−1Aixij=0∑k1+1Bj(ac×x+c−b−1)ji=0∑m−1((i+1)k2−ik2)x=0∑nxk1−i=0∑k2−1Aij=0∑k1+1Bjx=0∑m−1xi(ac×x+c−b−1)jf(a,b,c,n,k1,k2)=i=0∑m−1((i+1)k2−ik2)x=0∑nxk1−i=0∑k2−1j=0∑k1+1Ai×Bj×f(a′,b′,c′,m−1,i,j)

隻要預處理各種插值系數，即可在

log

⁡

O(T \times k ^4 \log n)

O(T×k4logn)的複雜度内求解。

#include<bits/stdc++.h>
using namespace std;
const int MAXN = 15;
const int P = 1e9 + 7;
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
template <typename T> void chkmax(T &x, T y) {x = max(x, y); }
template <typename T> void chkmin(T &x, T y) {x = min(x, y); } 
template <typename T> void read(T &x) {
x = 0; int f = 1;
char c = getchar();
for (; !isdigit(c); c = getchar()) if (c == '-') f = -f;
for (; isdigit(c); c = getchar()) x = x * 10 + c - '0';
x *= f;
}
template <typename T> void write(T x) {
if (x < 0) x = -x, putchar('-');
if (x > 9) write(x / 10);
putchar(x % 10 + '0');
}
template <typename T> void writeln(T x) {
write(x);
puts("");
}
struct info {
int a[MAXN][MAXN];
};
int sum[MAXN][MAXN];
int binom[MAXN][MAXN];
int power(int x, int y) {
if (y == 0) return 1;
int tmp = power(x, y / 2);
if (y % 2 == 0) return 1ll * tmp * tmp % P;
else return 1ll * tmp * tmp % P * x % P;
}
void update(int &x, int y) {
x += y;
if (x >= P) x -= P;
}
info func(int n, int a, int b, int c) {
assert(n >= 0 && a >= 0 && b >= 0 && c >= 0);
info ans;
memset(ans.a, 0, sizeof(ans.a));
if (a == 0 || (1ll * a * n + b) / c == 0) {
for (int k1 = 0; k1 <= 10; k1++) {
int mul = 0, now = 1;
for (int i = 0; i <= k1 + 1; i++) {
update(mul, 1ll * now * sum[k1][i] % P);
now = 1ll * now * n % P;
      }
int base = (1ll * a * n + b) / c % P; now = 1;
for (int k2 = 0; k1 + k2 <= 10; k2++) {
ans.a[k1][k2] = 1ll * now * mul % P;
now = 1ll * now * base % P;
      }
    }
return ans;
  }
if (a >= c) {
info tmp = func(n, a % c, b, c);
for (int k1 = 0; k1 <= 10; k1++)
for (int k2 = 0; k1 + k2 <= 10; k2++) {
int now = 1, base = a / c;
for (int i = 0; i <= k2; i++) {
update(ans.a[k1][k2], 1ll * binom[k2][i] * now % P * tmp.a[k1 + i][k2 - i] % P);
now = 1ll * now * base % P;
      }
    }
return ans;
  }
if (b >= c) {
info tmp = func(n, a, b % c, c);
for (int k1 = 0; k1 <= 10; k1++)
for (int k2 = 0; k1 + k2 <= 10; k2++) {
int now = 1, base = b / c;
for (int i = 0; i <= k2; i++) {
update(ans.a[k1][k2], 1ll * binom[k2][i] * now % P * tmp.a[k1][k2 - i] % P);
now = 1ll * now * base % P;
      }
    }
return ans;
  }
int m = (1ll * a * n + b) / c;
info tmp = func(m - 1, c, c - b - 1, a);
for (int k1 = 0; k1 <= 10; k1++) {
int all = 0, now = 1;
for (int i = 0; i <= k1 + 1; i++) {
update(all, 1ll * now * sum[k1][i] % P);
now = 1ll * now * n % P;
    }
for (int k2 = 0; k1 + k2 <= 10; k2++) {
ans.a[k1][k2] = 1ll * power(m, k2) * all % P;
for (int i = 0; i <= k2 - 1; i++)
for (int j = 0; j <= k1 + 1; j++) {
int coef = 1ll * binom[k2][i] * sum[k1][j] % P;
update(ans.a[k1][k2], P - 1ll * coef * tmp.a[i][j] % P);
      }
    }
  }
return ans;
}
void Lagrange(int n, int *x, int *y, int *a) {
static int p[MAXN], q[MAXN];
memset(p, 0, sizeof(p)); p[0] = 1;
for (int i = 1; i <= n; i++) {
for (int j = i - 1; j >= 0; j--) {
p[j + 1] = (p[j + 1] + p[j]) % P;
p[j] = (P - 1ll * p[j] * x[i] % P) % P;
    }
  }
for (int i = 1; i <= n; i++) {
memset(q, 0, sizeof(q));
for (int j = n - 1; j >= 0; j--)
q[j] = (p[j + 1] + 1ll * q[j + 1] * x[i]) % P;
int now = 1;
for (int j = 1; j <= n; j++)
if (j != i) now = 1ll * now * (x[i] - x[j]) % P;
now = power((P + now) % P, P - 2);
for (int j = 0; j <= n; j++)
q[j] = 1ll * q[j] * now % P;
for (int j = 0; j <= n; j++)
a[j] = (a[j] + 1ll * q[j] * y[i]) % P;
  }
}
int main() {
for (int i = 0; i <= 10; i++) {
static int pos[MAXN], now[MAXN];
now[0] = power(0, i);
for (int j = 1; j <= i + 2; j++) {
now[j] = (now[j - 1] + power(j, i)) % P;
pos[j] = j;
    }
Lagrange(i + 2, pos, now, sum[i]);
  }
for (int i = 0; i <= 10; i++) {
binom[i][0] = 1;
for (int j = 1; j <= i; j++)
binom[i][j] = binom[i - 1][j - 1] + binom[i - 1][j];
  }
int T; read(T);
while (T--) {
int n, a, b, c, k1, k2;
read(n), read(a), read(b), read(c), read(k1), read(k2);
writeln(func(n, a, b, c).a[k1][k2]);
  }
return 0;
}

#138. 類歐幾裡得算法