BZOJ2194: 快速傅立叶之二(NTT,卷积)

bzoj2194 飞快傅立叶之二

2194: 急迅傅立叶之二

Time Limit: 10 Sec  Memory
Limit: 259 MB

2179: FFT飞速傅立叶

Time Limit: 10 Sec  Memory
Limit: 259 MB

难题轮廓

给定 class=”math inline”>\(S(n,m)\)表示第二类Sterling数,定义函数 class=”math inline”>\(f(n)\)
\[f(n) =
\sum_{i=0}^n\sum_{j=0}^iS(i,j)*2^j*(j!)\]
给定正整数\(n,(n\leq
10^5)\),求\(f(n)\)

Time Limit:10 SecMemory Limit:259 MB
Submit:1776Solved:1055
[Submit][Status][Discuss]

2194: 急迅傅立叶之二

Time Limit:10 SecMemory Limit:259 MB
Submit:829Solved:474
[Submit][Status][Discuss]

Description

请计算C[k]=sigma(a[i]*b[i-k]) 在那之中 k < = i < n ,並且有 n
< = 10 ^ 5。 a,b中的元素均为小于等于100的非负整数。

 

Description

交付四个n位10进制整数x和y,你须要总计x*y。

题解

笔者们都精通第二类Sterling数的递推公式为
\[S(i,j) = S(i-1,j-1) + j*S(i-1,j),(1
\leq j \leq i-1)\]
且有境界\(S(i,i) = 1(0 \leq i),S(i,0) = 0(1
\leq i)\)
其次类Sterling数\(S(i,j)\)的意义是把\(i\)个成分划分成\(j\)个冬天的汇聚的方案
借使允许空集结的留存的话,方案即为\(m^n\)
大家使用容斥原理,枚举至少有微微空集结空会集,那么有
\[S(n,m) =
\frac{1}{m!}\sum_{k=1}^{m}C_m^k(m-k)^n(-1)^k\]
设\(g(n) =
\sum_{i=0}^nS(n,i)2^i(i!)\)
那正是说大家将\(S(n,m)\)代入\(g(n)\)化简得
\[g(n) =
\sum_{m=0}^n2^m(m!)\sum_{k=0}^m\frac{(-1)^k}{k!}\frac{(m-k)^n}{(m-k)!}\]
那么将\(g(n)\)带入答案表明式中,有
\[ans =
\sum_{n=0}^x\sum_{m=0}^n2^m(m!)\sum_{k=0}^m\frac{(-1)^k}{k!}\frac{(m-k)^n}{(m-k)!}\]
那时候大家发现每一遍最外层的\(n ->
(n+1)\)时,都约等于在其间的\(\frac{(m-k)^n}{(m-k)!}\)一项上又增加了八个\(\frac{(m-k)^{n+1}}{(m-k)!}\)
因此大家把这一项做等比数列求和.
设\(g(x) = \frac{x^{n+1} –
x}{(x-1)(x!)}\)
那就是说上式变成了
\[ans =
\sum_{m=0}^n2^m(m!)\sum_{k=0}^m\frac{(-1)^k}{k!}g(m-k)\]
于是大家在\(\sum_{k=0}^m\frac{(-1)^k}{k!}g(m-k)\)进行FFT总括卷积
与上述同类就只剩下了八个sigma式,for循环一边就能够.
复杂度\(O(nlogn)\)
图片 1

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
template<typename T>inline void read(T &x){
    x=0;char ch;bool flag = false;
    while(ch=getchar(),ch<'!');if(ch == '-') ch=getchar(),flag = true;
    while(x=10*x+ch-'0',ch=getchar(),ch>'!');if(flag) x=-x;
}
const int maxn = 600010;
const int mod = 998244353;
const int pri_rt = 3;
int w[maxn];
inline int qpow(int x,int p){
    int ret = 1;
    for(;p;p>>=1,x=1LL*x*x%mod) if(p&1) ret=1LL*ret*x % mod;
    return ret;
}
inline void FNT(int *x,int n,int p){
    for(int i=0,t=0;i<n;++i){
        if(i > t) swap(x[i],x[t]);
        for(int j=n>>1;(t^=j)<j;j>>=1);
    }
    for(int m=2;m<=n;m<<=1){
        int k = m>>1;
        int wn = qpow(pri_rt,p == 1 ? (mod-1)/m : (mod-1) - (mod-1)/m);
        for(int i=1;i<k;++i) w[i] = 1LL*w[i-1]*wn % mod;
        w[0] = 1;
        for(int i=0;i<n;i+=m){
            for(int j=0;j<k;++j){
                int u = 1LL*x[i+j+k]*w[j] % mod;
                x[i+j+k] = x[i+j] - u;
                if(x[i+j+k] < 0) x[i+j+k] += mod;
                x[i+j] += u;
                if(x[i+j] >= mod) x[i+j] -= mod;
            }
        }
    }
    if(p == -1){
        int inv = qpow(n,mod-2);
        for(int i=0;i<n;++i) x[i] = 1LL*x[i]*inv % mod;
    }
}
int fac[maxn],inv[maxn];
inline void init(int n){
    fac[0] = 1;
    for(int i=1;i<=n;++i) fac[i] = 1LL*fac[i-1]*i % mod;
    inv[n] = qpow(fac[n],mod-2);
    for(int i = n-1;i>=0;--i) inv[i] = 1LL*inv[i+1]*(i+1) % mod;
}
int A[maxn],B[maxn];
int main(){
    int n;read(n);
    int len;for(len=1;len <= (n+1);len<<=1);len<<=1;
    init(n);
    for(int i=0;i<=n;++i){
        if(i&1) A[i] = -inv[i] + mod;
        else A[i] = inv[i];
    }
    for(int i=2;i<=n;++i){
        B[i] = qpow(i,n+1) - i + mod;
        if(B[i] < 0) B[i] += mod;
        B[i] = (1LL*B[i]*qpow(i-1,mod-2)%mod*inv[i]) % mod;
    }B[1] = n;
    FNT(A,len,1);FNT(B,len,1);
    for(int i=0;i<len;++i) A[i] = 1LL*A[i]*B[i] % mod;
    FNT(A,len,-1);
    int ans = 1;
    for(int i=1,f2=2;i<=n;++i){
        ans = (ans + 1LL*A[i]*f2%mod*fac[i]) % mod;
        f2 = (f2<<1) % mod;
    }printf("%d\n",ans);
    getchar();getchar();
    return 0;
}

Description

请计算C[k]=sigma(a[i]*b[i-k]) 当中 k < = i < n ,而且有 n
< = 10 ^ 5。 a,b中的成分均为小于等于100的非负整数。

Description

请计算C[k]=sigma(a[i]*b[i-k]) 当中 k < = i < n ,况且有 n
< = 10 ^ 5。 a,b中的成分均为小于等于100的非负整数。

Input

率先行一个卡尺头N,接下去N行,第i+2..i+N-1行,每行三个数,依次表示a[i],b[i]
(0 < = i < N)。

Input

率先行贰个正整数n。 第二行描述叁个位数为n的正整数x。
第三行描述七个位数为n的正整数y。

Input

率先行二个整数N,接下去N行,第i+2..i+N-1行,每行五个数,依次表示a[i],b[i]
(0 < = i < N)。

Input

第一行二个整数N,接下去N行,第i+2..i+N-1行,每行八个数,依次表示a[i],b[i]
(0 < = i < N)。

Output

出口N行,每行一个整数,第i行输出C[i-1]。

Output

输出一行,即x*y的结果。

Output

出口N行,每行贰个整数,第i行输出C[i-1]。

Output

输出N行,每行贰个整数,第i行输出C[i-1]。

Sample Input

5
3 1
2 4
1 1
2 4
1 4

Sample Input

1
3
4

Sample Input

5
3 1
2 4
1 1
2 4
1 4

Sample Input

5
3 1
2 4
1 1
2 4
1 4

Sample Output

24
12
10
6
1

Sample Output

12

多少范围:
n<=60000

Sample Output

24
12
10
6
1

Sample Output

24
12
10
6
1

一看题的名字就清楚那题什么形式了。

那不是卷积的样式,所以我们要把五个数组中的多少个扭转,那样就改成卷积的花样了。然后就能够用FFT化解了。

形如c[k]=∑(a[i]*b[i-k])的样式得以将贰个数组翻转,再用FFT。

#include
#include
#include
#include
#include
#include
#define F(i,j,n) for(int i=j;i<=n;i++)
#define D(i,j,n) for(int i=j;i>=n;i--)
#define ll long long
#define maxn 300005
using namespace std;
int n,m,len,rev[maxn];
struct cp
{
    double x,y;
    inline cp operator +(cp b){return (cp){x+b.x,y+b.y};}
    inline cp operator -(cp b){return (cp){x-b.x,y-b.y};}
    inline cp operator *(cp b){return (cp){x*b.x-y*b.y,x*b.y+y*b.x};}
}a[maxn],b[maxn],c[maxn];
const double PI=acos(-1.0);
inline int read()
{
 int x=0,f=1;char ch=getchar();
 while (ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}
 while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
 return x*f;
}
inline void fft(cp *x,int n,int flag)
{
 F(i,0,n-1) if (rev[i]>i) swap(x[i],x[rev[i]]);
 for(int m=2;m<=n;m<<=1)
 {
  cp wn=(cp){cos(2.0*PI/m*flag),sin(2.0*PI/m*flag)};
  for(int i=0;i>1;
   F(j,0,mid-1)
   {
    cp u=x[i+j],v=x[i+j+mid]*w;
    x[i+j]=u+v;x[i+j+mid]=u-v;
    w=w*wn;
   }
  }
 }
    if(flag==-1) F(i,0,n-1) x[i].x/=n;
}
int main()
{
 n=read();int nn=n;
 F(i,0,n-1) a[i].x=read(),b[n-1-i].x=read();
 n=2*n-1;m=1;
 while (m>=1;}
  rev[i]=ret;
 }
 fft(a,n,1);fft(b,n,1);
 F(i,0,n-1) c[i]=a[i]*b[i];
 fft(c,n,-1);
 F(i,0,nn-1) printf("%lld\n",(ll)(c[nn-1+i].x+0.5));
 return 0;
}

飞快傅立叶之二 2194: 火速傅立叶之二 Time
Limit:10 SecMemory Limit:259 MB Submit:829Solved:474
[Submit][Status][Discuss] Description
请计算C[k]=sigma(a[i]*…

HINT

 

HINT

#include<map>
#include<cmath>
#include<queue>
#include<cstdio>
#include<complex>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define cp complex<double>
#define ll long long
#define PI acos(-1)
#define N 200010
int n,m,c[N],L=-1,r[N];
cp a[N],b[N];
char s1[N],s2[N];
void FFT(cp *x,int f)
{
    for(int i=0;i<n;i++) if(i<r[i]) swap(x[i],x[r[i]]);
    for(int i=1;i<n;i<<=1)
    {
        cp wn(cos(PI/i),f*sin(PI/i));
        for(int j=0;j<n;j+=(i<<1))
        {
            cp w(1,0),X,Y;
            for(int k=0;k<i;k++,w*=wn)
            {
                X=x[j+k];Y=w*x[j+k+i];
                x[j+k]=X+Y;x[j+k+i]=X-Y;
            }
        }
    }
}
int main()
{
    scanf("%d%s%s",&n,s1,s2);n--;
    for(int i=0;i<=n;i++) a[i]=s1[n-i]-'0';
    for(int i=0;i<=n;i++) b[i]=s2[n-i]-'0';
    m=n<<1;for(n=1;n<=m;n<<=1) L++;
    for(int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<L);
    FFT(a,1);FFT(b,1);
    for(int i=0;i<=n;i++) a[i]*=b[i];
    FFT(a,-1);
    for(int i=0;i<=m;i++) c[i]=(int)(a[i].real()/n+0.1);
    for(int i=0;i<=m;i++)
    {
        if(c[i]>9)
        {
            c[i+1]+=c[i]/10;
            c[i]%=10;if(i==m)m++;
        }
    }
    for(int i=m;i>=0;i--) printf("%d",c[i]);
    return 0;
}

 

HINT

Source

#include<map>
#include<cmath>
#include<queue>
#include<cstdio>
#include<complex>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define cp complex<double>
#define inf 1000000007
#define ll long long
#define PI acos(-1)
#define N 400010
inline int rd()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}
cp a[N],b[N];
int c[N],n,m,L=-1,r[N];
void FFT(cp *x,int f)
{
    for(int i=0;i<n;i++) if(i<r[i]) swap(x[i],x[r[i]]);
    for(int i=1;i<n;i<<=1)
    {
        cp wn(cos(PI/i),f*sin(PI/i));
        for(int j=0;j<n;j+=(i<<1))
        {
            cp w(1,0),X,Y;
            for(int k=0;k<i;k++,w*=wn)
            {
                X=x[j+k];Y=w*x[i+j+k];
                x[j+k]=X+Y;x[i+j+k]=X-Y;
            }
        }
    }
}
int main()
{
    n=rd()-1;
    for(int i=0;i<=n;i++){a[i]=rd();b[n-i]=rd();}
    m=n<<1;for(n=1;n<=m;n<<=1) L++;
    for(int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<L);
    FFT(a,1);FFT(b,1);
    for(int i=0;i<n;i++) a[i]*=b[i];
    FFT(a,-1);
    for(int i=m/2;i<=m;i++) printf("%d\n",(int)(a[i].real()/n+0.1));
    return 0;
}

 

Source

难题中给的公式不佳搞

作者们依照套路,将$B$翻转一下

*$$C = \sum_0^n a_i \ b_{n – 1 – i

  • k}$$**

那时候背后的姿态就只与$k$有关了

设$$D(n – 1 + k) = \sum_0^n a_i *
b_{n – 1 – i + k}$$

直接NTT

#include<cstdio>#define swap x ^= y, y ^= x, x ^= y#define LL long long using namespace std;const int MAXN = 3 * 1e5 + 10;inline int read() {    char c = getchar(); int x = 0, f = 1;    while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}    while(c >= '0' && c <= '9') x = x * 10 + c - '0',c = getchar();    return x * f;}const int P = 998244353, g = 3, gi = 332748118;int N;int LL a[MAXN], b[MAXN], r[MAXN];LL fastpow(LL a, int p, int mod) {    LL base = 1;    while {        if(p & 1) base = (base * a) % mod;        a =  % mod; p >>= 1;    }    return base % mod;}LL NTT(LL *A, int type, int N, int mod) {    for(int i = 0; i < N; i++)         if(i < r[i]) swap(A[i], A[r[i]]);    for(int mid = 1; mid < N; mid <<= 1) {        LL W = fastpow( (type == 1) ? g : gi, (P - 1) / (mid << 1), mod );         for(int j = 0; j < N; j += (mid << 1)) {            int w = 1;             for(int k = 0; k < mid; k++, w =  % P) {                LL x = A[j + k] % P, y = w * A[j + k + mid] % P;                A[j + k] =  % P;                A[j + k + mid] = (x - y + P) % P;            }         }    }    if(type == -1) {        LL inv = fastpow(N, mod - 2, mod);        for(int i = 0; i < N; i++)            A[i] = (A[i] * inv) % mod;    }}int main() {    #ifdef WIN32    freopen("a.in","r",stdin);    #endif    N = read();    for(int i = 0; i < N; i++)         a[i] = read(), b[N - i] = read();    int limit = 1, L = 0;    while(limit <= N + N) limit <<=1, L++;    for(int i = 0; i < limit; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));    NTT(a, 1, limit, P); NTT(b, 1, limit, P);    for(int i = 0; i <  limit; i++) a[i] = (a[i] * b[i]) % P;    NTT(a, -1, limit, P);    for(int i = 0; i < N * 2; i++)        printf("%d\n",a[i] % P);    return 0;} 

相关文章