[快速数论变换 NTT]

先粘一个模板。这是求高精度乘法的

#include <bits/stdc++.h>
#define maxn 1010
using namespace std;

char s[maxn];

typedef long long ll;
ll A[maxn], B[maxn];

const int md = 998244353, G = 3;

int n;

ll power_mod(ll a, ll b){
	ll ret = 1;
	while(b > 0){
		if(b & 1)ret = ret * a % md;
		b >>= 1;
		a = a * a % md;
	}return ret;
}

void NTT(ll A[], int n, int type){
	for(int i = 0, j = 0; i < n; i ++){
		if(i < j)swap(A[i], A[j]);
		for(int t = n >> 1; (j ^= t) < t; t >>= 1);
	}

	for(int k = 2; k <= n; k <<= 1){
		ll wn = type > 0 ? power_mod(G, (md - 1) / k) : power_mod(G, (md - 1) - (md - 1) / k);
		for(int i = 0; i < n; i += k){
			ll w = 1;
			for(int j = 0; j < k >> 1; j ++){
				ll T = w * A[i+j+(k>>1)] % md;
				A[i+j+(k>>1)] = (A[i+j] - T) % md;
				A[i+j] = (A[i+j] + T) % md;
				w = w * wn % md;
			}
		}
	}

	if(type < 0){
		ll inv = power_mod(n, md - 2);
		for(int i = 0; i < n; i ++)
		    A[i] = A[i] * inv % md;
	}
}

int main(){
	freopen("mul.in", "r", stdin);
	freopen("mul.out", "w", stdout);
	scanf("%s", s);
	int n1 = strlen(s);
	for(int i = n1 - 1; ~i; i --)
		A[n1-i-1] = s[i] ^ 48;
    scanf("%s", s);
	int n2 = strlen(s);
	for(int i = n2 - 1; ~i; i --)
		B[n2-i-1] = s[i] ^ 48;
	for(n = 1; n <= n1 + n2; n <<= 1);
	NTT(A, n, 1), NTT(B, n, 1);
	for(int i = 0; i < n; i ++)
	    A[i] = A[i] * B[i] % md;
	NTT(A, n, -1);
	for(int i = 0; i < n; i ++){
		A[i] = (A[i] + md) % md;
		A[i+1] += A[i] / 10, A[i] %= 10;
	}
	while(n && A[n] == 0) n --;
	for(int i = n; i >= 0; i --)
	    printf("%lld", A[i]);
	return 0;
}

[SDOI 2015]序列统计

取指标变成加法问题,构造生成函数F(x),F(x)的n次方就是答案。

注意指标要mod M-1

还有注意0的问题,如果x!=0就不用管数据中的0,否则需要两种情况求和(然而数据中x并不等于0我就没写)

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#define maxn 50010
using namespace std;
typedef long long ll;
const int md = 1004535809, G = 3;

int N, M, X, S, a[maxn], Log[maxn], n;
ll F[maxn], H[maxn], ret[maxn];

int p[maxn], primes, vis[maxn];

ll power_mod(ll a, ll b, ll md){
	ll ret = 1;
	while(b > 0){
		if(b & 1) ret = ret * a % md;
		b >>= 1;
		a = a * a % md;
	}return ret;
}

void pre_prime(){
	int n = 20000;
	for(int i = 2; i <= n; i ++){
		if(!vis[i])p[++ primes] = i;
		for(int j = 1; j <= primes; j ++){
			if(i * p[j] > n)break;
			vis[i * p[j]] = true;
			if(i % p[j] == 0)break;
		}
	}
}

int getG(int M){
	pre_prime();
	int ret = 2;
	for(; ; ret ++){
		bool flag = true;
		for(int i = 1; i <= primes; i ++){
			if(M-1 < p[i])break;
			if((M-1) % p[i] == 0){
				if(power_mod(ret, (M - 1) / p[i], M) == 1){
					flag = false;
					break;
				}
			}
		}
		if(flag)return ret;
	}
}

void NTT(ll A[], int n, int type){
	for(int i = 0, j = 0; i < n; i ++){
		if(i > j)swap(A[i], A[j]);
		for(int t = n >> 1; (j ^= t) < t; t >>= 1);
	}
	
	for(int k = 2; k <= n; k <<= 1){
		ll wn = type > 0 ? power_mod(G, (md-1)/k, md) : power_mod(G, md-1-(md-1)/k, md);
		for(int i = 0; i < n; i += k){
			ll w = 1;
			for(int j = 0; j < k >> 1; j ++){
				ll T = w * A[i+j+(k>>1)];
				A[i+j+(k>>1)] = (A[i+j] - T + md) % md;
				A[i+j] = (A[i+j] + T) % md;
				w = w * wn % md;
			}
		}
	}
	
	if(type < 0){
		ll inv = power_mod(n, md - 2, md);
		for(int i = 0; i < n; i ++)
			A[i] = A[i] * inv % md;
	}
}

void CF(ll A[]){
	for(int i = 0; i < n; i ++)H[i] = F[i];
	NTT(A, n, 1), NTT(H, n, 1);
	for(int i = 0; i < n; i ++)
	    A[i] = A[i] * H[i] % md;
	NTT(A, n, -1);
	for(int i = (M - 1); i < n; i ++)
		(A[i % (M - 1)] += A[i]) %= md, A[i] = 0;
}

void power(int b){
	ret[0] = 1;
	while(b){
		if(b & 1)CF(ret);
		b >>= 1;
		CF(F);
	}
}

int main(){
	scanf("%d%d%d%d", &N, &M, &X, &S);
	for(int i = 1; i <= S; i ++)
		scanf("%d", &a[i]);

	for(n = 1; n <= M; n <<= 1); n <<= 1;
	
	ll nw = 1, Gn = getG(M); 
	for(int i = 0; i < M-1; i ++)
		Log[nw] = i, nw = nw * Gn % M;

	for(int i = 1; i <= S; i ++)
	    if(a[i])F[Log[a[i]]] = 1;

	power(N);
	printf("%lld\n", (ret[Log[X]] + md) % md);
	return 0;
}

[COGS 2287]疯狂的机器人

这里有Catalan数的链接

#include <bits/stdc++.h>
#define maxn 300010
using namespace std;
typedef long long ll;
int n;
const int md = 998244353, G = 3;

ll power_mod(ll a, ll b = md - 2){
	ll ret = 1;
	while(b > 0){
		if(b & 1)ret = ret * a % md;
		b >>= 1;
		a = a * a % md;
	}return ret;
}

ll fac[maxn], inv[maxn], f[maxn];

void NTT(ll* A, int n, int type){
	for(int i = 0, j = 0; i < n; i ++){
		if(i > j)swap(A[i], A[j]);
		for(int t = n >> 1; (j ^= t) < t; t >>= 1);
	}
	
	for(int k = 2; k <= n; k <<= 1){
		ll wn = power_mod(G, type > 0 ? (md-1) / k : md-1 - (md-1) / k);
		for(int i = 0; i < n; i += k){
			ll w = 1;
			for(int j = 0; j < k >> 1; j ++){
				ll T = w * A[i+j+(k>>1)] % md;
				A[i+j+(k>>1)] = (A[i+j] - T) % md;
				A[i+j] = (A[i+j] + T) % md;
				w = w * wn % md;
			}
		}
	}
	
	if(type == -1){
		ll inv = power_mod(n);
		for(int i = 0; i < n; i ++)
			A[i] = A[i] * inv % md;
	}
}

int main(){
	freopen("crazy_robot.in", "r", stdin);
	freopen("crazy_robot.out", "w", stdout);
	scanf("%d", &n);
	fac[0] = inv[0] = 1;
	for(int i = 1; i <= 2 * n; i ++)
	    fac[i] = fac[i-1] * i % md;
	inv[2 * n] = power_mod(fac[2 * n]);
	for(int i = 2 * n - 1; i >= 1; i --)
		inv[i] = inv[i+1] * (i+1) % md;
	for(int i = 0; i <= n; i ++){
		if(i & 1){f[i] = 0;continue;}
		f[i] = fac[i] * inv[i / 2] % md * inv[i / 2 + 1] % md * inv[i] % md;
	}
	
	int _N = 1;
	for(; _N <= n + n; _N <<= 1);
	NTT(f, _N, 1);
	for(int i = 0; i < _N; i ++)
	    f[i] = f[i] * f[i] % md;
	NTT(f, _N, -1);
	ll ans = 0;
	for(int i = 0; i <= n; i ++)
		(ans += inv[n - i] * f[i]) %= md;
	cout << (ans * fac[n] % md + md) % md << endl;
	return 0;
}

[NOI十连测3007]卡常大法好。

#include <cstring>
#include <cstdio>
#define maxn 150010

const int md = 998244353, G = 3;

int R;
#define fastcall __attribute__((optimize("-Os")))
#define IL __inline__ __attribute__((always_inline))
fastcall IL int mul_mod(int a, int b){
    __asm__ __volatile__ ("\tmull %%ebx\n\tdivl %%ecx\n" :"=d"(R):"a"(a),"b"(b),"c"(md));
    return R;
}

int n, N;

int f[maxn], g[maxn], h[maxn], fac[maxn], inv[maxn], pw[maxn], tmp[maxn];

fastcall IL int power_mod(int a, long long b = md - 2){
	int ret = 1;
	while(b > 0){
		if(b & 1)ret = mul_mod(ret, a);
		b >>= 1;
		a = mul_mod(a, a);
	}return ret;
}

fastcall IL void swap(int& a, int& b){
	a ^= b ^= a ^= b;
}

fastcall IL void NTT(int* A, int n, int type){
	for(int i = 0, j = 0; i < n; i ++){
		if(i > j)swap(A[i], A[j]);
		for(int t = n >> 1; (j ^= t) < t; t >>= 1);
	}

	for(int k = 2; k <= n; k <<= 1){
		int wn = power_mod(G, type > 0 ? (md-1)/k : (md-1)-(md-1)/k);
		for(int i = 0, m = k >> 1; i < n; i += k){
			int w = 1;
			for(int j = 0; j < k >> 1; j ++){
				int T = mul_mod(w, A[i+j+m]);
				A[i+j+m] = (A[i+j] - T + md) % md;
				A[i+j] = (A[i+j] + T) % md;
				w = mul_mod(w, wn);
			}
		}
	}

	if(type < 0){
		int inv = power_mod(n);
		for(int i = 0; i < n; i ++)
		    A[i] = mul_mod(A[i], inv);
	}
}

fastcall void Getinv(int* a, int* b, int n){
	if(n == 1){b[0] = power_mod(a[0]); return;}
	Getinv(a, b, n >> 1);
	memcpy(tmp, a, sizeof(a[0]) * n);
	NTT(tmp, n << 1, 1);
	NTT(b, n << 1, 1);
	for(int i = 0; i < n << 1; i ++)
	    b[i] = mul_mod(b[i], (2ll - mul_mod(b[i], tmp[i]) + md) % md);
	NTT(b, n << 1, -1);
	memset(b + n, 0, sizeof(a[0]) * n);
}

int F[30010][16];
int K;

fastcall void cdq(int l, int r){
	if(l == r)return;
	int mid = l + r >> 1, len = r - l + 1;
	cdq(l, mid);
	for(n = 1; n <= len; n <<= 1);
	for(int i = 0; i < n; i ++)f[i] = h[i] = 0;
	for(int i = l; i <= mid; i ++)
	    f[i-l] = mul_mod(F[i][K-1] + F[i][K], inv[i]);
    h[0] = 0;
	for(int i = 1; i < n; i ++)
	    h[i] = mul_mod(g[i], inv[i-1]);
	if(n<=32) {
		for(int i = 0; i < n; i ++)tmp[i] = 0;
		for(int i = 0 ; i < n ; ++ i) {
			for(int j = 0 ; j <= i ; ++ j) {
				tmp[i] += mul_mod(h[j], f[i-j]);
				if(tmp[i]>=md) tmp[i] -= md;
			}
		}
		for(int i = 0 ; i < n ; ++ i) h[i] = tmp[i];
	} else {
		NTT(f, n, 1), NTT(h, n, 1);
		for(int i = 0; i < n; i ++)
	    	h[i] = mul_mod(h[i], f[i]);
		NTT(h, n, -1);
	}
	for(int i = mid + 1; i <= r; i ++){
	    F[i][K] += mul_mod(h[i-l], fac[i-1]);
	    if(F[i][K] >= md)F[i][K] -= md;
	}
	cdq(mid + 1, r);
}

int main(){
	N = 30000;
	for(n = 1; n <= N; n <<= 1);
	fac[0] = inv[0] = 1;
	for(int i = 1; i <= n; i ++)fac[i] = mul_mod(fac[i - 1], i);
	inv[n] = power_mod(fac[n]);
	for(int i = n - 1; i >= 1; i --)inv[i] = mul_mod(inv[i + 1], i + 1);
	for(long long i = 0; i <= n; i ++)pw[i] = power_mod(2, i * (i - 1) / 2);
	for(int i = 1; i < n; i ++)h[i] = mul_mod(pw[i], inv[i - 1]);
	for(int i = 0; i < n; i ++)f[i] = mul_mod(pw[i], inv[i]);
	Getinv(f, g, n);
	n <<= 1;
	NTT(g, n, 1);
	NTT(h, n, 1);
	for(int i = 0; i < n; i ++)
	    g[i] = mul_mod(g[i], h[i]);
	NTT(g, n, -1);
	for(int i = 1; i <= N; i ++)
	    g[i] = mul_mod(g[i], fac[i - 1]);
 
	for(int i = 0; i <= N; i ++)F[i][0] = pw[i];
	for(int i = 1; i <= 15; i ++)K = i, cdq(0, N);

	static int S[20][20];
	S[0][0] = 1;
	for(int i = 1; i <= 15; i ++)
	    for(int j = 1; j <= i; j ++){
	        S[i][j] = mul_mod(S[i-1][j], j) + S[i-1][j-1];
	        if(S[i][j] >= md)S[i][j] -= md;
		}

	int test, n, m;
	scanf("%d", &test);
	while(test --){
		scanf("%d%d", &n, &m);
		int ans = 0;
		for(int i = 0; i <= m; i ++){
			ans += mul_mod(mul_mod(S[m][i], F[n][i]), fac[i]);
			if(ans >= md)ans -= md;
		}
		printf("%d\n", ans);
	}
	return 0;
}

  

  

posted @ 2016-04-27 10:08  _Horizon  阅读(965)  评论(0编辑  收藏  举报