OJ-Problems-Source/.ACM-Templates/.Not classified/数学数论.cpp
2016-08-13 23:35:41 +08:00

635 lines
22 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

//快速幂
ll powMod(ll a, ll b, ll m) {
ll r = 1;
for (a %= m; b; b >>= 1) { if (b & 1) { r = r * a % m; } a = a * a % m; }
return r;
}
//素数筛
//Eratosthenes O(nloglogn)
const int N = 10000000; //~110ms
bitset<N> isprime;
void getPrime() {
isprime.set(); isprime[0] = isprime[1] = false;
for (int i = 2; i < N; i++) {
if (isprime[i]) {
for (ll j = (ll)i * i; j < N; j += i) { isprime[j] = false; }
}
}
}
//Euler O(n) prime[0]为个数
const int N = 10000000; //~110ms
int prime[N]; //3711111 for [2, 10^9)
void getPrime() {
for (int i = 2; i < N; i++) {
if (!prime[i]) { prime[++prime[0]] = i; }
for (int j = 1; j <= prime[0] && prime[j] * i < N; j++) {
prime[prime[j] * i] = 1;
if (i % prime[j] == 0) { break; }
}
}
}
//Euler O(n)
const int N = 10000000; //~95ms
int prime[N >> 3]; bitset<N> isprime;
void getPrime() {
isprime.set(); isprime[0] = isprime[1] = false;
for (int i = 2; i < N; i++) {
if (isprime[i]) { prime[++prime[0]] = i; }
for (int j = 1; j <= prime[0] && prime[j] * i < N; j++) {
isprime[prime[j] * i] = false;
if (i % prime[j] == 0) { break; }
}
}
}
//[a, b]区间内素数个数
bitset<N> isprime, isprimesmall;
ll segPrime(ll a, ll b) {
ll ret = 0; isprime.set(); isprimesmall.set();
for (int i = 2; (ll)i * i <= b; i++) {
if (isprimesmall[i]) {
for (ll j = (ll)i * i; (ll)j * j <= b; j += i) { isprimesmall[j] = false; }
for (ll j = max(2ll, (a + i - 1) / i)) * i; j <= b; j += i) { isprime[j - a] = false; }
}
}
for (ll i = 0; i <= r - l; i++) { ret += isprime[i]; }
return ret;
}
//分解质因数
ll factor[100], facCnt;
void getFactors(ll x) {
facCnt = 0;
for (int i = 2, xx = sqrt(x + 0.5); i <= xx; i++) {
if (x % i == 0) {
factor[facCnt++] = i;
while (x % i == 0) { x /= i; }
}
}
if (x != 1) { factor[facCnt++] = x; }
}
//分解质因数及个数 预处理素数表
ll factor[100][2], facCnt;
void getFactors(ll x) {
facCnt = 0;
for (int i = 1; prime[i] <= x / prime[i]; i++) {
factor[facCnt][1] = 0;
if (x % prime[i] == 0) {
factor[facCnt][0] = prime[i];
while (x % prime[i] == 0) { factor[facCnt][1]++; x /= prime[i]; }
facCnt++;
}
}
if (x != 1) { factor[facCnt][0] = x; factor[facCnt++][1] = 1; }
}
//Miller-Rabin素性测试 素数返回true 错误(伪素数)概率为1/4^Times
const int Times = 7;
int WIT[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022}; //7 bases for n < 2^64
ll mulMod(ll a, ll b, ll m) {
ll r = 0;
for (a %= m, b %= m; b; b >>= 1) { if (b & 1) { r = (r + a) % m; } a = (a << 1) % m; }
return r;
}
ll powMod(ll a, ll b, ll m) {
ll r = 1;
for (a %= m; b; b >>= 1) { if (b & 1) { r = mulMod(r, a, m); } a = mulMod(a, a, m); }
return r;
}
bool Miller_Rabin(ll n) {
if (n == 2) { return true; }
if (n < 2 || (n & 1) == 0) { return false; }
ll m = n - 1; int k = 0;
while ((m & 1) == 0) { k++; m >>= 1; }
for (int i = 0; i < Times; i++) {
ll a = WIT[i], x = powMod(a, m, n), y = 0;
//ll a = rand() % (n - 1) + 1;
for (int j = 0; j < k; j++, x = y) {
y = mulMod(x, x, n);
if (y == 1 && x != 1 && x != n - 1) { return false; }
}
if (y != 1) { return false; }
}
return true;
}
//pollard rho质因素分解
//对n进行素因子分解, 存入factor, k设置为107左右即可
ll factor[100], facCnt; //质因素分解结果(无序)
ll pollard_rho(ll x, ll c) {
ll i = 1, k = 2, x0 = rand() % (x - 1) + 1, y = x0;
while (true) {
x0 = (mulMod(x0, x0, x) + c) % x;
ll d = llabs(__gcd(y - x0, x));
if (d != 1 && d != x) { return d; }
if (y == x0) { return x; }
if (++i == k) { y = x0; k <<= 1; }
}
}
void findfac(ll n, int k = 107) {
if (n == 1) { return; }
if (Miller_Rabin(n)) { factor[facCnt++] = n; return; }
ll p = n; int c = k;
while (p >= n) { p = pollard_rho(p, c--); } //k值变化, 防止死循环
findfac(p, k); findfac(n / p, k);
}
//求单个数的欧拉函数
ll eular(ll n) {
ll ret = 1;
while ((n & 1) == 0) { n >>= 1; ret <<= 1; }
for (ll i = 3; i * i <= n; i += 2) {
if (n % i == 0) { n /= i; ret *= i - 1; while (n % i == 0) { n /= i; t *= i; } }
}
return n > 1 ? ret * (n - 1) : ret;
}
//欧拉函数筛 O(nloglogn)
const int N = 10000000; //~400ms
int phi[N] = {0, 1};
void getPhi() {
for (int i = 2; i < N; i++) {
if (!phi[i]) {
for (int j = i; j < N; j += i) {
if (!phi[j]) { phi[j] = j; } phi[j] -= phi[j] / i;
}
}
}
}
//素数 + 欧拉函数筛 O(n)
const int N = 10000000; //~150ms
int prime[N >> 3], phi[N] = {0, 1}; bitset<N> isprime;
void getPrimePhi() {
isprime.set(); isprime[0] = isprime[1] = false;
for (int i = 2; i < N; i++) {
if (isprime[i]) { prime[++prime[0]] = i; phi[i] = i - 1; }
for (int j = 1, k; j <= prime[0] && prime[j] * i < N; j++) {
isprime[k = prime[j] * i] = false;
if (i % prime[j] == 0) { phi[k] = phi[i] * prime[j]; break; }
phi[k] = phi[i] * (prime[j] - 1);
}
}
}
//素数 + 莫比乌斯函数筛 O(n)
const int N = 10000000; //150ms
int prime[N >> 3], miu[N] = {0, 1}; bitset<N> isprime;
void getPrimeMiu() {
isprime.set(); isprime[0] = isprime[1] = false;
for (int i = 2; i < N; i++) {
if (isprime[i]) { prime[++prime[0]] = i; miu[i] = -1; }
for (int j = 1, k; j <= prime[0] && prime[j] * i < N; j++) {
isprime[k = prime[j] * i] = false;
if (i % prime[j] == 0) { miu[k] = 0; break; }
miu[k] = -miu[i];
}
}
}
//素数 + 欧拉函数 + 莫比乌斯函数筛 O(n)
const int N = 10000000; //~230ms
int prime[N >> 3], phi[N] = {0, 1}, miu[N] = {0, 1}; bitset<N> isnprime;
void getPrimePhiMiu() {
for (int i = 2; i < N; i++) {
if (!isnprime[i]) { prime[++prime[0]] = i; phi[i] = i - 1; miu[i] = -1; }
for (int j = 1, k; j <= prime[0] && prime[j] * i < N; j++) {
isnprime[k = prime[j] * i] = true;
if (i % prime[j] == 0) { phi[k] = phi[i] * prime[j]; miu[k] = 0; break; }
phi[k] = phi[i] * (prime[j] - 1); miu[k] = -miu[i];
}
}
}
//素数 + 约数个数筛 O(n)
const int N = 10000000; //~200ms
bitset<N> isprime;
int prime[N >> 3], faccnt[N] = {0, 1}, d[N]; //d[i]表示i的最小质因子的幂次
void getPrimeFaccnt() {
isprime.set(); isprime[0] = isprime[1] = false;
for (int i = 2; i < N; i++) {
if (isprime[i]) { prime[++prime[0]] = i; faccnt[i] = 2; d[i] = 1; }
for (int j = 1, k; j <= prime[0] && prime[j] * i < N; j++) {
isprime[k = prime[j] * i] = false;
if (i % prime[j] == 0) {
faccnt[k] = faccnt[i] / (d[i] + 1) * (d[i] + 2); d[k] = d[i] + 1; break;
}
faccnt[k] = faccnt[i] << 1; d[k] = 1;
}
}
}
//A^B的约数之和为:
//sum = [1+p1+p1^2+...+p1^(a1*B)]*...*[1+pn+pn^2+...+pn^(an*B)].
//等比数列求和 1+a+a^2+...+a^b
ll sumPow(ll a, ll b, ll m) {
ll r = 1; a %= m;
for (ll t = 1; b; b >>= 1) {
if (b & 1) { r = (r * a + t) % m; }
t = t * (a + 1) % m;
a = a * a % m;
}
return r;
}
//求逆元(ax = 1(mod m)的x值)
//扩展欧几里得(求ax + by = gcd(a, b)的解), 求出的x为a对b的模逆元
ll exgcd(ll a, ll b, ll &x, ll &y) {
if (b == 0) { x = 1; y = 0; return a; }
ll d = exgcd(b, a % b, y, x); y -= a / b * x; return d;
}
//解不定方程ax + by = c 求得的只是其中一组解
//对于不定整数方程ax + by = c, 若c mod gcd(a, b) = 0, 则该方程存在整数解, 否则不存在整数解
//在找到ax + by = gcd(a, b)的一组解x0, y0后可得到ax + by = c的一组解x1 = x0 * (c / gcd(a, b)), y1 = y0 * (c / gcd(a,b))
//ax + by = c的其他整数解满足
//x = x1 + b / gcd(a, b) * t, y = y1 - a / gcd(a, b) * t(其中t为任意整数)
//求ax + by = c的一组解
bool linear_equation(int a, int b, int c, int &x, int &y) {
int d = exgcd(a, b, x, y);
if (c % d) { return false; }
int k = c / d; x *= k; y *= k; return true;
}
//求ax = b (mod p)循环节内的所有解
ll ans[N], cnt;
bool linear_equation(ll a, ll b, ll p) {
ll x, y, d = exgcd(a, p, x, y); cnt = 0;
if (b % d) { return false; }
x = (x % p + p) % p;
ans[++cnt] = x * (b / d) % (p / d);
for (int i = 1; i < d; i++) { ans[++cnt] = (ans[1] + i * p / d) % n; }
return true;
}
//线性预处理逆元
ll Inv[N] = {1, 1};
void getInv(int m) {
for (ll i = 2; i < m; i++) { Inv[i] = (m - m / i) * Inv[m % i] % m; }
}
//扩展欧几里得求逆元
ll modReverse(ll a, ll m) {
ll x, y, d = exgcd(a, m, x, y);
if (d == 1) { return (x % m + m) % m; } else { return -1; }
}
//费马小定理, m为素数, a与m互质
ll inv(ll a, ll m) { return powMod(a, m - 2, m); }
//只能求0 < a < m的情况,a和m互质
ll inv(ll a, ll m) {
if (a == 1) { return 1; }
return inv(m % a, m) * (m - m / a) % m;
}
//中国剩余定理 求模线性方程组x = a[i] (mod m[i]) m[i]可以不互质
//[1, n]内解的个数为(n - x) / m1 + (x != 0)
bool merge(ll a1, ll m1, ll a2, ll m2, ll &a3, ll &m3) {
ll d = __gcd(m1, m2), c = a2 - a1;
if (c % d != 0) { return false; }
c = (c % m2 + m2) % m2 / d; m1 /= d; m2 /= d;
c = c * inv(m1, m2) % m2 * m1 * d + a1;
m3 = m1 * m2 * d; a3 = (c % m3 + m3) % m3;
return true;
}
ll CRT(ll a[], ll m[], int k) {
ll a1 = a[0], m1 = m[0];
for (int i = 1; i < k; i++) {
ll a2 = a[i], m2 = m[i], m3, a3;
if (!merge(a1, m1, a2, m2, a3, m3)) { return -1; }
a1 = a3; m1 = m3;
}
return (a1 % m1 + m1) % m1;
}
//模线性方程组 需扩展欧几里得
//求解方程ax ≡ b (mod m) 相当于求解方程ax + my = b (x, y为整数)
int m[10], a[10]; //模数为m, 余数为a, X % m = a
bool solve(int &m0, int &a0, int m, int a) {
ll y, x, d = exgcd(m0, m, x, y);
if (abs(a - a0) % d) { return false; }
x *= (a - a0) / d; x %= m / d;
a0 = (x * m0 + a0); m0 *= m / d; a0 %= m0;
if (a0 < 0) { a0 += m0; }
return true;
}
//无解返回false, 有解返回true
//解的形式最后为a0 + m0 * t (0 <= a0 < m0)
bool MLES(int &m0, int &a0, int n) { //解为X = a0 + m0 * k
bool flag = true; m0 = 1; a0 = 0;
for (int i = 0; i < n; i++) {
if (!solve(m0, a0, m[i], a[i])) { flag = false; break; }
}
return flag;
}
//求原根
ll fac[N];
ll getRoot(ll n) {
int cnt = 0;
for (ll i = 2; i * i < n - 1; i++) { if ((n - 1) % i == 0) { fac[cnt++] = i; fac[cnt++] = (n - 1) / i; } }
for (int i = 2, j;; i++) {
for (j = 0; j < cnt; j++) { if (powMod(i, fac[j], n) == 1) { break; } }
if (j == cnt) { return i; }
}
}
//线性基
//异或线性基
//若要查询第k小子集异或和, 则把k写成二进制, 对于是1的第i位, 把从低位到高位第i个不为0的数异或进答案
//若要判断是否有非空子集的异或和为0, 如果不存在自由基, 那么说明只有空集的异或值为0, 需要高斯消元来判断
struct XORBase {
int a[64];
void clear() { memset(a, 0, sizeof(a)); }
void ins(ll x) {
for (int i = 62; i >= 0; i--) {
if (x & (1 << i)) {
if (a[i]) { x ^= a[i]; }
else { a[i] = x; }
break;
}
}
}
//查询最大子集异或和
void query() {
ll ret = 0;
for (int i = 62; i >= 0; i--) { ret = max(ret, ret ^ a[i]); }
return ret;
}
};
//实数线性基
//ins返回要插入的数是否可以被之前的数线性表示出来, 返回true表示不能, false表示可以
int m;
struct Base {
double a[N][N]; bool v[N];
void clear() { memset(a, 0, sizeof(a)); memset(v, 0, sizeof(v)); }
bool ins(double *x) {
for (int i = 0; i < m; i++) {
if (fabs(x[i]) > 1e-6) {
if (v[i]) {
double t = x[i] / a[i][i];
for (int j = 0; j < m; j++) { x[j] -= t * a[i][j]; }
} else {
v[i] = 1;
for (int j = 0; j < m; j++) { a[i][j] = x[j]; }
return true;
}
}
}
return false;
}
};
//离散对数 大步小步算法 Baby-Step Giant-Step
//BSGS(a, b, p): 求ax = b (mod p)的最小非负整数解, 若无解则返回 -1
//rev(a, p): 扩展欧几里得求逆元
//powMod(base, pow, mod): 快速幂
//mulMod(a, b, mod): 快速乘(这里用快速乘是为了避免爆long long int, 实际有时可以不用)
unordered_map<ll, ll> Hash;
ll BSGS(ll a, ll b, ll p) {
if (b >= p) { return -1; }
a %= p, b %= p;
if (!a && !b) { return 1; } //a和b都是p的倍数的话, 就相当于0^x = 0 (mod p)了, 那么最小非负整数解就是1
if (!a) { return -1; } //如果a是p的倍数但是b不是, 就相当于0^x = t (mod p), t > 0, 无解
Hash.clear();
ll m = ceil(sqrt(p)), tmp = 1 % p; //tmp = a^j
for (ll j = 0; j < m; j++) { //预处理出a^j mod p的值
Hash[tmp] = j; tmp = mulMod(tmp, a, p);
}
tmp = rev(powMod(a, m, p), p); //tmp = a^(-m)
for (ll i = 0; i < m; i++) {
if (Hash.find(b) != Hash.end()) { return i * m + Hash[b]; }
b = mulMod(b, tmp, p);
}
return -1;
}
//高斯消元 求线性方程组的解 O(n^3)
//有equ个方程, var个变元, 增广矩阵行数为equ, 列数为var + 1, 下标从0开始
int a[N][N], x[N]; //增广矩阵, 解集
int freex[N], freenum;//自由变元 (多解枚举自由变元可以使用)
//返回值为-1表示无解, 为0是唯一解, 否则返回自由变元个数
int Gauss(int equ, int var) {
int mxrow, col, k; freenum = 0;
for (k = 0, col = 0; k < equ && col < var; k++, col++) {
mxrow = k;
for (int i = k + 1; i < equ; i++) { if (abs(a[i][col]) > abs(a[mxrow][col])) { mxrow = i; } }
if (a[mxrow][col] == 0) { k--; freex[freenum++] = col; continue; } //自由变元
if (mxrow != k) { for (int j = col; j <= var; j++) { swap(a[k][j], a[mxrow][j]); } }
for (int i = k + 1; i < equ; i++) {
if (a[i][col]) {
int x = abs(a[i][col]), y = abs(a[k][col]), lcm = x / __gcd(x, y) * y, tx = lcm / x, ty = lcm / y;
if (a[i][col] * a[k][col] < 0) { ty = -ty; }
for (int j = col; j <= var; j++) {
a[i][j] = a[i][j] * tx - a[k][j] * ty; //a[i][j] = (a[i][j] % M + M) % M;
}
}
}
}
for (int i = k; i < equ; i++) { if (a[i][col]) { return -1; } } //无解
if (k < var) { return var - k; } //自由变元个数
for (int i = var - 1; i >= 0; i--) { //唯一解,回代
for (int j = i + 1; j < var; j++) {
if (a[i][j]) { a[i][var] -= a[i][j] * x[j]; /*a[i][var] = (a[i][var] % M + M) % M;*/ }
}
//while (a[i][var] % a[i][i]) { a[i][var] += M; }
x[i] = a[i][var] / a[i][i]; //x[i] = (a[i][var] * inv(a[i][i], M)) % M;
}
return 0;
}
//高斯消元 (浮点数)
const double eps = 1e-9;
const int N = 205;
double a[N][N], x[N]; //方程的左边的矩阵和等式右边的值, 求解之后x存的就是结果
int equ, var; //方程数和未知数个数
//返回0表示无解, 1表示有解
int Gauss() {
int i, j, k, col, mxrow;
for (k = 0, col = 0; k < equ && col < var; k++, col++) {
mxrow = k;
for (i = k + 1; i < equ; i++) {
if (fabs(a[i][col]) > fabs(a[mxrow][col])) { mxrow = i; }
}
if (fabs(a[mxrow][col]) < eps) { return 0; }
if (k != mxrow) {
for (j = col; j < var; j++) { swap(a[k][j], a[mxrow][j]); }
swap(x[k], x[mxrow]);
}
x[k] /= a[k][col];
for (j = col + 1; j < var; j++) { a[k][j] /= a[k][col]; }
a[k][col] = 1;
for (i = 0; i < equ; i++) {
if (i != k) {
x[i] -= x[k] * a[i][k];
for (j = col + 1; j < var; j++) { a[i][j] -= a[k][j] * a[i][col]; }
a[i][col] = 0;
}
}
}
return 1;
}
//自适应simpson积分
//给定一个函数f(x), 求[a, b]区间内f(x)到x轴所形成区域的面积
double simpson(double l, double r) {return (f(l) + f(r) + 4 * f((l + r) / 2.0)) * (r - l) / 6.0;}
double rsimpson(double l, double r) {
double mid = (l + r) / 2.0;
if (fabs(simpson(l, r) - simpson(l, mid) - simpson(mid, r)) < EPS) {
return simpson(l, mid) + simpson(mid, r);
}
return rsimpson(l, mid) + rsimpson(mid, r);
}
//FFT O(nlogn)
//以下n必须为2的幂, op为1时是求DFT, op为-1时为求IDFT
typedef complex<double> comp;
const double PI = acos(-1.0);
void fft(comp a[], int n, int op) {
for (int i = 1, j = 0; i < n - 1; i++) {
for (int s = n; j ^= s >>= 1, ~j & s;);
if (i < j) { swap(a[i], a[j]); }
}
for (int i = 1; i < n; i <<= 1) {
comp wn(cos(PI / i), op * sin(PI / i));
for (int j = 0; j < n; j += i << 1) {
comp w(1, 0);
for (int k = 0; k < i; k++, w *= wn) {
comp x = a[j + k], y = w * a[i + j + k];
a[j + k] = x + y; a[i + j + k] = x - y;
}
}
}
if (op == -1) { for (int i = 0; i < n; i++) { a[i] = comp(a[i].real() / n, a[i].imag()); } }
}
//求高精度乘法 HDU 1402
comp a[N], b[N];
char str1[N / 2], str2[N / 2];
int sum[N];
int main() {
while (~scanf("%s%s", str1, str2)) {
memset(a, 0, sizeof(a)); memset(b, 0, sizeof(b));
int n = strlen(str1), m = strlen(str2), len = 1;
while (len < n * 2 || len < m * 2) { len <<= 1; }
for (int i = 0; i < n; i++) { a[i] = comp(str1[n - 1 - i] - '0', 0); }
for (int i = 0; i < m; i++) { b[i] = comp(str2[m - 1 - i] - '0', 0); }
fft(a, len, 1); fft(b, len, 1);
for (int i = 0; i < len; i++) { a[i] *= b[i]; }
fft(a, len, -1);
for (int i = 0; i < len; i++) { sum[i] = (int)(a[i].real() + 0.5); }
for (int i = 0; i < len; i++) { sum[i + 1] += sum[i] / 10; sum[i] %= 10; }
len = n + m - 1;
while (sum[len] <= 0 && len > 0) { len--; }
for (int i = len; i >= 0; i--) { putchar(sum[i] + '0'); } puts("");
}
}
//NTT O(nlogn)
//998244353 = 119 * 2^23 + 1, 原根为3; 1004535809 = 479 * 2^21 + 1, 原根为3
//786433 = 3 * 2^18 + 1, 原根为10; 880803841 = 105 * 2^23 + 1, 原根为26
//P是素数且N必须是P - 1的因子
//op为1时是求FNT, op为-1时为求IFNT
const int P = 998244353, G = 3, N = 262144, K = 17;
ll g[N + 5], ng[N + 5], Inv[N + 5] = {1, 1};
void initG() {
g[K] = powMod(G, (P - 1) / N, P); ng[K] = powMod(g[K], P - 2, P);
for (int i = K - 1; i >= 0; i--) { g[i] = g[i + 1] * g[i + 1] % P; ng[i] = ng[i + 1] * ng[i + 1] % P; }
for (ll i = 2; i <= N; i++) { Inv[i] = (P - P / i) * Inv[P % i] % P; }
}
void ntt(ll a[], int n, int op) {
for (int i = 1, j = 0; i < n - 1; i++) {
for (int s = n; j ^= s >>= 1, ~j & s;);
if (i < j) { swap(a[i], a[j]); }
}
for (int d = 0; (1 << d) < n; d++) {
int m = 1 << d; ll w0 = op == 1 ? g[d] : ng[d];
for (int i = 0; i < n; i += m << 1) {
for (int j = 0, w = 1; j < m; j++, w = w * w0 % P) {
ll &x = a[i + j + m], &y = a[i + j], t = w * x % P;
x = y - t; y = y + t;
if (x < 0) { x += P; } if (y >= P) { y -= P; }
}
}
}
if (op == -1) { for (int i = 0; i < n; i++) { a[i] = a[i] * Inv[n] % P; } }
}
//多项式求逆元 O(nlogn) 即求B(x)满足A(X) * B(x) = 1 (mod x^n), deg(B) <= deg(A)
void polyInv(ll a[], ll b[], int n) {
if (n == 1) { b[0] = powMod(a[0], P - 2, P); return; }
polyInv(a, b, n >> 1);
int k = 1;
while (k < n << 1) { k <<= 1; }
for (int i = 0; i < n; i++) { tmp[i] = a[i]; }
for (int i = n; i < k; i++) { tmp[i] = b[i] = 0; }
ntt(tmp, k, 1); ntt(b, k, 1);
for (int i = 0; i < k; i++) {
b[i] = b[i] * (2 - tmp[i] * b[i] % P) % P;
if (b[i] < 0) { b[i] += P; }
}
ntt(b, k, -1);
for (int i = n; i < k; i++) { b[i] = 0; }
}
//多项式除法 O(nlogn) 即求D(x)和R(x)满足A(x) = D(x) * B(x) + R(x), deg(D) <= deg(A) - deg(B), deg(R) < deg(B)
ll a0[N], b0[N];
void polyDiv(ll a[], int n, ll b[], int m, ll d[], ll r[]) {
int k = 1, t = n - m + 1;
while (k < t << 1) { k <<= 1; }
for (int i = 0; i < k; i++) { a0[i] = b0[i] = 0; }
for (int i = 0; i < m; i++) { a0[i] = b[m - i - 1]; }
polyInv(a0, b0, t);
for (int i = t; i < k; i++) { b0[i] = 0; }
for (int i = 0; i < t; i++) { a0[i] = a[n - i - 1]; }
for (int i = t; i < k; i++) { a0[i] = 0; }
ntt(b0, k, 1); ntt(a0, k, 1);
for (int i = 0; i < k; i++) { a0[i] = a0[i] * b0[i] % P; }
ntt(a0, k, -1);
reverse(a0, a0 + t);
for (int i = 0; i < t; i++) { d[i] = a0[i]; }
for (k = 1; k < n << 1; k <<= 1);
for (int i = t; i < k; i++) { a0[i] = 0; }
for (int i = 0; i < m; i++) { b0[i] = b[i]; }
for (int i = m; i < k; i++) { b0[i] = 0; }
ntt(a0, k, 1); ntt(b0, k, 1);
for (int i = 0; i < k; i++) { a0[i] = a0[i] * b0[i] % P; }
ntt(a0, k, -1);
for (int i = 0; i < m; i++) { r[i] = (a[i] - a0[i]) % P; }
for (int i = m; i < k; i++) { r[i] = 0; }
}
//多项式求对数函数 O(nlogn) a[0] = 1
void polyLn(ll a[], ll b[], int n) {
polyInv(a, tmp2, n);
int k = 1;
while (k < n << 1) { k <<= 1; }
for (int i = 0; i < n - 1; i++) { b[i] = a[i + 1] * (i + 1) % P; }
for (int i = n - 1; i < k; i++) { b[i] = 0; }
ntt(b, k, 1); ntt(tmp2, k, 1);
for (int i = 0; i < k; i++) { b[i] = b[i] * tmp2[i] % P; }
ntt(b, k, -1);
for (int i = n - 1; i >= 0; i--) { b[i] = b[i - 1] * Inv[i] % P; } b[0] = 0;
}
//多项式求指数函数 O(nlogn) a[0] = 0
void polyExp(ll a[], ll b[], int n) {
if (n == 1) { b[0] = 1; return; }
polyExp(a, b, n >> 1); polyLn(b, tmp, n);
int k = 1;
while (k < n << 1) { k <<= 1; }
for (int i = 0; i < n; i++) { tmp[i] -= a[i]; if (tmp[i] < 0) { tmp[i] += P; } }
if (++tmp[0] == P) { tmp[0] = 0; }
for (int i = n; i < k; i++) { tmp[i] = b[i] = 0; }
ntt(tmp, k, 1); ntt(b, k, 1);
for (int i = 0; i < k; i++) { b[i] = b[i] * tmp[i] % P; }
ntt(b, k, -1);
for (int i = n; i < k; i++) { b[i] = 0; }
}
//多项式求平方根 O(nlogn) a[0] = 1
void polySqrt(ll a[], ll b[], int n) {
if (n == 1) { b[0] = 1; return; }
polySqrt(a, b, n >> 1); polyInv(b, tmp2, n);
int k = 1;
while (k < n << 1) { k <<= 1; }
for (int i = 0; i < n; i++) { tmp[i] = a[i]; }
for (int i = n; i < k; i++) { tmp[i] = b[i] = 0; }
ntt(tmp, k, 1); ntt(b, k, 1); ntt(tmp2, k, 1);
for (int i = 0; i < k; i++) { b[i] = (b[i] * b[i] + tmp[i]) % P * Inv[2] % P * tmp2[i] % P; }
ntt(b, k, -1);
for (int i = n; i < k; i++) { b[i] = 0; }
}
//快速沃尔什变换 即求C[i] = sum{j ? k = i}(A[j] * B[k]) ?为任一二元逻辑位运算
void fwt(ll a[], int n) {
for (int d = 1; d < n; d <<= 1) {
for (int i = 0, m = d << 1; i < n; i += m) {
for (int j = 0; j < d; j++) {
ll x = a[i + j], y = a[i + j + d];
//xor: a[i + j] = x + y; a[i + j + d] = x - y;
//and: a[i + j] = x + y;
//or: a[i + j + d] = x + y;
}
}
}
}
void ufwt(ll a[], int n) {
for (int d = 1; d < n; d <<= 1) {
for (int i = 0, m = d << 1; i < n; i += m) {
for (int j = 0; j < d; j++) {
ll x = a[i + j], y = a[i + j + d];
//xor: a[i + j] = (x + y) >> 1; a[i + j + d] = (x - y) >> 1;
//and: a[i + j] = x - y;
//or: a[i + j + d] = y - x;
}
}
}
}