zoukankan      html  css  js  c++  java
  • [快速数论变换 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
    ", (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__ ("	mull %%ebx
    	divl %%ecx
    " :"=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
    ", ans);
    	}
    	return 0;
    }
    

      

      

    给时光以生命,而不是给生命以时光。
  • 相关阅读:
    Python 资源大全中文版
    Life is short.,You need Python
    哪些 Python 库让你相见恨晚?
    中国裁判文书网全网最新爬虫分析
    关于pycharm导入其他项目时出现找不到python无法运行的问题
    禅道项目管理软件配置及使用教程
    curl
    fusionpbx 中文 汉化
    kafka operation
    golang包管理工具——glide
  • 原文地址:https://www.cnblogs.com/Candyouth/p/5437789.html
Copyright © 2011-2022 走看看