使用 Python 的 FFT 和 IFFT 编码热方程解的问题（python fft函数）

25-04-26 1

如果您想了解使用Python的FFT和IFFT编码热方程解的问题的相关知识，那么本文是一篇不可错过的文章，我们将对pythonfft函数进行全面详尽的解释，并且为您提供关于AtCoderAGC019E

如果您想了解使用 Python 的 FFT 和 IFFT 编码热方程解的问题的相关知识，那么本文是一篇不可错过的文章，我们将对python fft函数进行全面详尽的解释，并且为您提供关于AtCoder AGC019E Shuffle and Swap (DP、FFT、多项式求逆、多项式快速幂)、BZOJ.4503. 两个串 (FFT/bitset)、CF 827E Rusty String FFT、CF993E Nikita and Order Statistics 【fft】的有价值的信息。

本文目录一览：

使用 Python 的 FFT 和 IFFT 编码热方程解的问题（python fft函数）
AtCoder AGC019E Shuffle and Swap (DP、FFT、多项式求逆、多项式快速幂)
BZOJ.4503. 两个串 (FFT/bitset)
CF 827E Rusty String FFT
CF993E Nikita and Order Statistics 【fft】

使用 Python 的 FFT 和 IFFT 编码热方程解的问题（python fft函数）

如何解决使用 Python 的 FFT 和 IFFT 编码热方程解的问题

我正在尝试使用 Python 的快速傅立叶变换 (FFT) 来实现热方程（一维）的解决方案。目标是一种有效的数值方案，用于说明在半无限域边界处指定的与时间相关的温度。这在传热文献中被称为 Duhamel 问题（参见，例如 Carslaw 和 Jaeger，1959，第 2.5 节）。对于指定为时间函数的边界处的强迫，已知精确解。在本案例中，我寻求对作为定期收集的离散数据给出的边界历史的响应。

对热方程进行傅里叶变换和适当的边界条件 ( ![equation]https://latex.codecogs.com/gif.download?T%280%2Ct%29%20%3D%20f%28t%29%3B%20%5Clim_%7Bx%20%5Crightarrow%20%5Cinfty%7D%20T%28x%2Ct%29%20%3D%200 ) 和初始条件 ( ![equation]https://latex.codecogs.com/gif.download?T%28x%2C0%29%20%3D%200 ) 相对于时间，问题的解由 ![equation]https://latex.codecogs.com/gif.download?%5Ctilde%7BT%7D%20%3D%20%5Ctilde%7Bf%7D%28%5Calpha%29%20%5Cexp%20%5B-%28i%20%5Calpha%20/%20%5Ckappa%29%5E%7B1/2%7D%20x%5D 给出，其中波浪号表示函数的傅立叶变换，alpha 是变换变量（频率），kappa 是热扩散系数。所以，方案是取边界数据的FFT，乘以 ![equation]https://latex.codecogs.com/gif.download?%5Cexp%20%5B%20-%20%28i%20%5Calpha%20/%20%5Ckappa%29%5E%7B1/2%7D%20x%5D ，并反转（使用 IFFT）以获得 T。

我构建了一个基准问题来验证我的编码。特别是，我考虑了强制函数 ![equation]https://latex.codecogs.com/gif.download?f%28t%29%20%3D%20a%20%5Cexp%20%5B%20-%20b%5E2%20%28t-t_0%29%5E2%5D ，其中 a、b 和 t_0 是固定常量。该函数可以代入经典的 Duhamel 解，结果可以通过数值积分来评估。此外，该函数的FT具有简单的封闭形式，可以代入上一段给出的解，并且可以通过直接数值积分对解的FT进行反演。这两个独立评估的结果是相同的，验证了上面给出的 FT 解决方案。

我已尝试为上述强制函数编写此问题，该函数由以固定时间间隔采样的离散数据给出。我调用Python的FFT例程来生成来自数据的强制函数的 FT。然后，我乘以 x 中的衰减指数，最后调用 IFFT 对给定的 x 固定值处的温度历史进行反演。该代码返回的值表现出一些正确的定性行为（例如，峰值温度的延迟到达；温度在后期逐渐下降），但在数量上（例如，温度过高）和质量上（例如，温度在早期不是渐近零）。我显然在对这个问题的编码中做错了什么，但我一直无法确定问题所在。我在衰减指数中指定频率的方式可能会发现一个可能的错误。

如果有人更熟悉 Python 的 FFT 和 IFFT 的使用，我将不胜感激。

我的代码如下：

def FTtrial3():
    import numpy as np
    import matplotlib.pyplot as plt
# nsam is the number of points given for the boundary data
    nsam = 32
    tmin = 0.
    tmax = 600.
    timedata = np.linspace(tmin,tmax,num=nsam)
# a,b,and t0 are fixed constants defining the boundary forcing
    a = 10.0
    b = 0.02
    time0 = 300.0 
    tempdata = a * np.exp(-b*b*(timedata - time0)**2) 
# take the FFT of the boundary data
    FTtemp = np.fft.fft(tempdata) 
# generate an array of the Fourier frequencies
    tint = (tmax-tmin)/(nsam-1)
    freq = np.fft.fftfreq(nsam,d=tint)
# specify the distance from the boundary for evaluation (in m)
    x = 0.1
# specify the thermal diffusivity (in m^2/s)
    D = 1.2E-05
# compute the FT of the solution
    FTsoln = FTtemp * np.exp(-np.sqrt( 1.j * freq / D) * x)
# invert to obtain the final solution
    soln = np.fft.ifft(FTsoln) 
    solnreal = np.real(soln)
# provide results of independent evaluation of Duhamel''s solution for comparison
    pltdhml = np.array( [ [200.0,2.063E-05],[250.0,2.558E-03],[300.0,6.229E-02],[330.0,2.112E-01],[350.0,3.721E-01],[370.0,5.546E-01],[400.0,7.877E-01],[440.0,9.293E-01],[450.0,9.375E-01],[460.0,9.380E-01],[470.0,9.324E-01],[500.0,8.924E-01],[550.0,7.949E-01],[600.0,6.970E-01],[650.0,6.112E-01],[700.0,5.388E-01],[800.0,4.275E-01],[900.0,3.483E-01],[1000.0,2.902E-01] ] )
    tdhml,Tempdhml = pltdhml.T
# generate plot for comparison of FFT solution to Duhamel (exact) solution
    plt.scatter(tdhml,Tempdhml,color = ''green'',marker = ''o'',label = "Duhamel")
    plt.scatter(timedata,solnreal,color = ''red'',marker = ''^'',label = "FT")
    plt.legend()
    plt.show()
FTtrial3()

AtCoder AGC019E Shuffle and Swap (DP、FFT、多项式求逆、多项式快速幂)

题目链接

https://atcoder.jp/contests/agc019/tasks/agc019_e

题解

tourist 的神仙 E 题啊做不来做不来…… 这题我好像想歪了啊 = =……

首先我们可以考虑，什么样的操作序列才是合法的？有用的位置只有两种，一种是两个序列在这个位置上都是 1, 称作 11 型，另一种是一个 0 一个 1, 称作 01 型。设两种位置分别有 $A$ 个和 $2B$ 个。考虑一个操作序列，交换两个 11 型相当于没交换，每个 11 型只会被交换两次，每个 01 型只会被交换一次。这也就是说，如果我们从 shuffle 之后的 $a_i$ 向 $b_i$ 连边，那么形成的图一定是若干个环加上 $m$ 条链，链的开头结尾都是 01 型，中间是 11 型。对于链来说，上面操作的顺序必须固定；对于环来说，上面操作的顺序可以任意。

下面有两种处理方式。

做法一

设 $dp [i][j]$ 表示把 $j$ 个无标号的 11 型放到 $i$ 条链中，可得 DP 式: $dp [i][j]=\sum_{k\ge 0}\frac {dp [i-1][j-k]}{(k+1)!}$, 其中分母的含义是链上 $(k+1)$ 个点顺序固定，最后的答案是 $A!B!(A+B)!\sum^A_{i=0} dp [B][i]$. $(A+B)!$ 表示将边随意排序，$A!B!$ 表示 11 型和 01 型点之间是有标号的。时间复杂度 $O (n^3)$.

但是我们发现这个 DP 就相当于在给多项式 $\sum_{n\ge 0}\frac {1}{(n+1)!} x^n$ 进行幂运算，于是用多项式快速幂加速即可，时间复杂度 $O (n\log^2n)$ 或 $O (n\log n)$.

做法二

有没有聪明一点的做法？有！设 $dp [i][j]$ 表示目前一共放了 $i$ 个 11 型和 $j$ 个 01 型链 (考虑已经放了的元素的标号，但是每次仅仅是往右添加)，我们强行转移！ $dp [i][j]=j^2\times dp [i][j-1]+ij\times dp [i-1][j]$ 前一个式子是要加一个新的 01 型链，选两个 01 型；后一个式子是要选一个链，再把这条链的结尾端点任意扩展一个位置。答案就是 $\sum_{k}{A\choose k}\times f [k][B]\times ((A-k)!)^2\times {A+B\choose A-k}$, 其中 $A+B\choose A-k$ 是选出位置，$A\choose k$ 是选出编号，$((A-k)!)^2$ 是求出组成环的方案数。时间复杂度 $O (n^2)$, 可以通过。

但是我们发现这个 DP 还可以用多项式优化！令 $g [i][j]=\frac {dp [i][j]}{(j!)^2i!}$, 显然有 $g [i][j]=g [i][j-1]+j\times g [i-1][j]$ 然后这个使用 NE Lattice Path 的方式来理解，就是从 $(0,0)$ 走到 $(i,j)$，每往上走一次路径权值乘上横坐标，求所有路径权值和。考虑另一种 DP，枚举在第 $i$ 列走几步，那么发现第 $i$ 列的生成函数就是 $\frac {1}{1-ix}$, 然后答案就是所有列生成函数之积于是可以分治 NTT + 多项式求逆计算，时间复杂度 $O (n\log^2n)$.

代码

做法二 $O (n^2)$

#include<cstdio>
#include<cstdlib>
#include<iostream>
#include<cassert>
#include<cstring>
#define llong long long
using namespace std;

const int N = 2e4;
const int P = 998244353;
const llong INV2 = 499122177;
llong fact[N+3],finv[N+3];

llong quickpow(llong x,llong y)
{
	llong cur = x,ret = 1ll;
	for(int i=0; y; i++)
	{
		if(y&(1ll<<i)) {y-=(1ll<<i); ret = ret*cur%P;}
		cur = cur*cur%P;
	}
	return ret;
}
llong mulinv(llong x) {return quickpow(x,P-2);}
llong comb(llong x,llong y) {return x<0||y<0||x<y ? 0ll : fact[x]*finv[y]%P*finv[x-y]%P;}

int n,a,b;
char s[N+3],t[N+3];
llong dp[2][N+3];

int main()
{
	fact[0] = 1ll; for(int i=1; i<=N; i++) fact[i] = fact[i-1]*i%P;
	finv[N] = quickpow(fact[N],P-2); for(int i=N-1; i>=0; i--) finv[i] = finv[i+1]*(i+1)%P;
	scanf("%s%s",s+1,t+1); n = strlen(s+1);
	for(int i=1; i<=n; i++)
	{
		if(s[i]==''1'' && t[i]==''1'') {a++;}
		else if(s[i]^t[i]) {b++;}
	}
	b>>=1;
	int cur = 0,prv = 1;
	dp[0][0] = 1ll;
	for(int j=1; j<=b; j++)
	{
		cur^=1; prv^=1;
		dp[cur][0] = dp[prv][0]*j*j%P;
		for(int i=1; i<=a; i++)
		{
			dp[cur][i] = (dp[prv][i]*j*j+dp[cur][i-1]*i*j)%P;
//			printf("dp[%d][%d]=%lld\n",i,j,dp[cur][i]);
		}
	}
	llong ans = 0ll;
	for(int k=0; k<=a; k++)
	{
		ans = (ans+dp[cur][k]*comb(a,k)%P*fact[a-k]%P*fact[a-k]%P*comb(a+b,a-k))%P;
//		printf("k%d ans%lld\n",k,ans);
	}
	printf("%lld\n",ans);
	return 0;
}

BZOJ.4503. 两个串 (FFT/bitset)

题目链接

$Description$

给定两个字符串 S 和 T，求 T 在 S 中出现了几次，以及分别在哪些位置出现。T 中可能有 ''?'' 字符，这个字符可以匹配任何字符。 $|S|,|T|\leq 10^5$。

$Solution$

FFT：

https://www.cnblogs.com/cjyyb/p/8798446.html

显然我们可以同 CF528D 一样枚举 $26$ 个字符然后跑 $FFT$。但字符集太大了，复杂度是 $O (26n\log n)$，以 $FFT$ 的常数肯定 GG。（然而 CF 上一道题 $n$ 更大但 $36n\log n$ 的 $FFT$ 还跑的非常轻松，这就是差距么...）

考虑没有通配符时能怎么做（不考虑 $KMP$）。把每个字符 $a\sim z$ 映射到 $1\sim26$，那么 $S_i,T_i$ 匹配了就是 $S_i=T_i$。考虑 $S,T$ 做差。但是每个位置需要有 $|T|$ 个值求和来表示匹配 $T$ 串的情况，就算不匹配正负相加也可能变成 $0$。所以考虑平方，即若 $f (x)=\sum_{i=0}^{|T|-1}(S_{x+i}-T_i)^2=0$，那么 $x$ 位置匹配了 $T$。把平方和拆开，就可以得到两个常数项和一个 $\sum_{i=0}^{|T|-1} S_{x+i} T_i$，把 $T$ 反转后就可以 $FFT$ 了。

如果有通配符呢？考虑如果 $T_i$ 是通配符，怎么让它不产生影响，也就是贡献是 $0$。那么令通配符 $T_i=0$，外面再乘个 $T_i$ 就可以了。即 f(x)=|T|−1∑i=0(Sx+i−Ti)2Ti&=|T|−1∑i=0S2x+iTi−2|T|−1∑i=0Sx+iT2i+|T|−1∑i=0T3i

第三部分是个常数，前两部分可以分别 $FFT$ 求出来。

//22280kb	3076ms(这慢的...)
#include <cmath>
#include <cctype>
#include <cstdio>
#include <cstring>
#include <algorithm>
//#define gc() getchar()
#define MAXIN 1000000
#define gc() (SS==TT&&(TT=(SS=IN)+fread(IN,1,MAXIN,stdin),SS==TT)?EOF:*SS++)
typedef long long LL;
const int N=(1<<18)+5;
const double PI=acos(-1);

int S[N],T[N],rev[N];
char IN[MAXIN],*SS=IN,*TT=IN;
struct Complex
{
	double x,y;
	Complex(double x=0,double y=0):x(x),y(y) {}
	Complex operator +(const Complex &a) {return Complex(x+a.x, y+a.y);}
	Complex operator -(const Complex &a) {return Complex(x-a.x, y-a.y);}
	Complex operator *(const Complex &a) {return Complex(x*a.x-y*a.y, x*a.y+y*a.x);}
}A[N],B[N],C[N],D[N];

void FFT(Complex *a,int lim,int opt)
{
	for(int i=1; i<lim; ++i) if(i<rev[i]) std::swap(a[i],a[rev[i]]);
	for(int i=2; i<=lim; i<<=1)
	{
		int mid=i>>1; Complex Wn(cos(PI/mid),opt*sin(PI/mid));
		for(int j=0; j<lim; j+=i)
		{
			Complex w(1,0),t;
			for(int k=j; k<j+mid; ++k,w=w*Wn)
				a[k+mid]=a[k]-(t=a[k+mid]*w), a[k]=a[k]+t;
		}
	}
	if(opt==-1) for(int i=0; i<lim; ++i) a[i].x/=lim;
}

int main()
{
	static int pos[N];

	register char c;
	int n=0; while(isalpha(c=gc())) S[n++]=c-96;
	int m=0; while(isalpha(c=gc())||c==''?'') T[m++]=c==''?''?0:c-96;
	std::reverse(T,T+m);//!
	for(int i=0; i<n; ++i) A[i]=Complex(S[i]*S[i],0),C[i]=Complex(S[i]<<1,0);
	int sumT=0;
	for(int i=0,t; i<m; ++i) B[i]=Complex(t=T[i],0),D[i]=Complex(t*t,0),sumT+=t*t*t;

	int lim=1,l=-1;
	while(lim<=n+m-2) lim<<=1,++l;//n+m-2就可以了...~~有点想不通...~~n-1次和m-1次咯
	for(int i=1; i<lim; ++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<l);
	FFT(A,lim,1), FFT(B,lim,1), FFT(C,lim,1), FFT(D,lim,1);
	for(int i=0; i<lim; ++i) A[i]=A[i]*B[i]-C[i]*D[i];//可以直接在这里计算然后只IDFT一次啊（反正都是多项式）
	FFT(A,lim,-1);

	int ans=0;
	for(int i=0; i<=n-m; ++i) if(!(int)(A[m+i-1].x+0.5+sumT)) pos[++ans]=i;
	printf("%d\n",ans);
	for(int i=1; i<=ans; printf("%d\n",pos[i++]));

	return 0;
}

bitset：

可以用 bitset 暴力匹配。维护 bitset<N> Ans 表示每个位置能否匹配 $T$ 串。枚举 $T$ 中的字符 $T_j$，维护 $Ans_i$ 是否能够匹配 $T_j$，也就是 $S_{i+j}$ 处是否是 $T_j$。可以用 $26$ 个 bitset<N> Pos[26] 得到每种字符在 $S$ 的哪些位置上出现过（出现过的位置设为 $1$）。然后枚举 $T_j$（$j\in [0,T|$）的时候，$Ans_i$ 每次都与上 $Pos_{T_j}$ 右移 $j$ 位就可以了。（若是通配符就不用再与了）（注意是右移）最后若 $Ans_i$ 还是 $1$，则 $i$ 位置匹配 $T$ 串。

复杂度 $O (\frac {nm}{w})$。

//2744kb	1372ms
#include <cctype>
#include <cstdio>
#include <bitset>
#include <algorithm>
//#define gc() getchar()
#define MAXIN 1000000
#define gc() (SS==TT&&(TT=(SS=IN)+fread(IN,1,MAXIN,stdin),SS==TT)?EOF:*SS++)
const int N=1e5+3;

char S[N],T[N];
std::bitset<N> Ans,pos[26];
char IN[MAXIN],*SS=IN,*TT=IN;

int main()
{
	static int p[N];

	register char c;
	int n=0; while(isalpha(S[n]=gc())) ++n;//会有个换行符啊=_= 所以while的时候n先别加。
	int m=0; while(isalpha(T[m]=gc())||T[m]==''?'') ++m;

	for(int i=0; i<n; ++i) pos[S[i]-97].set(i);
	Ans.set();
	for(int i=0; i<m; ++i) if(T[i]!=''?'') Ans&=pos[T[i]-97]>>i;

	int ans=0;
	for(int i=0; i<=n-m; ++i) if(Ans[i]==1) p[++ans]=i;
	printf("%d\n",ans);
	for(int i=1; i<=ans; printf("%d\n",p[i++]));

	return 0;
}

CF 827E Rusty String FFT

传送门

如果没有碍事的?的话，判定字符串的循环节直接用KMP的失配数组就可以搞定。现在有了碍事的?，我们就需要考虑更通用的算法。

考虑KMP失配数组判定字符串循环节的本质，发现判定$k$是否为字符串的循环节等价于判定字符串在右移$k$位后能否和原字符串匹配（只考虑二者重叠的部分）。

我们不妨先把?直接看成一个可以匹配任何字符的通配符，而解决带通配符的字符串匹配问题的一个算法就是FFT。

（以下默认下标为$0$~$n-1$）

设字符串为$s$，因为字符集只有$2$，不妨直接枚举不能匹配的两种情况。定义$a_i=[s_i=''V'']$，$b_i=[s_i=''K'']$，另一种情况只需要把定义反过来就行了。

再定义一个数组

$$\begin{align}c_i=\sum_{j=0}^{n-1}a_{i+j}b_j\end{align}$$

（另一种情况式子是一样的）

发现$s$右移$i$位后能和原串匹配当且仅当两种情况的$c_i$都为$0$。

而上面的式子把${a_i}$反转后就是

$$\begin{align}c_i=\sum_{j=0}^{n-1}a^R_{n-i-j-1}b_{j}\end{align}=\left(a^R*b\right)_{n-i-1}$$

跑两遍卷积即可。

然而还是会发现一些问题……从样例就能看出来，这里的?其实并不是通配符，因为在判定循环节时一个?在不同的位置必须代表相同的字符（可能有点抽象，参见第一组样例中２为什么可以匹配但不是循环节）。

然而这个问题其实并不棘手。我们都知道如果一个数真的是循环节那么它的所有倍数也一定是循环节，所以对于那些$c_i$为$００$的位置再判断一下它的倍数是否也都满足$c_i=0$就行了。

（完整证明参见官方题解）

这一步的复杂度是$O(n\log n)$的，不会影响到总复杂度。

另外，考虑到问题的特殊性，其实不必跑两遍卷积，任意跑其中一种情况即可，两种情况的$c_i$分别是$\left(a^Rb\right)_{n-i-1}$和$＋\left(a^Rb\right)_{n＋i-1}$（可以从两种情况的对称性的角度理解）。

（我比较懒写的是NTT）

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int maxn=1048600,p=998244353,g=3;
void NTT(int*,int,int);
int qpow(int,int);
char s[maxn];
bool ans[maxn];
int T,n,A[maxn],B[maxn];
int main(){
	scanf("%d",&T);
	while(T--){
		scanf("%d%s",&n,s);
		int N=1;
		while(N<(n<<1))N<<=1;
		memset(A,0,sizeof(int)*N);
		memset(B,0,sizeof(int)*N);
		for(int i=0;i<n;i++){
			if(s[i]==''V'')A[n-i-1]=1;
			else if(s[i]==''K'')B[i]=1;
		}
		NTT(A,N,1);
		NTT(B,N,1);
		for(int i=0;i<N;i++)A[i]=(long long)A[i]*B[i]%p;
		NTT(A,N,-1);
		int cnt=0;
		for(int i=1;i<=n;i++){
			ans[i]=true;
			for(int j=i;j<=n;j+=i)ans[i]&=!(A[n-j-1]+A[n+j-1]);
			cnt+=ans[i];
		}
		printf("%d\n",cnt);
		for(int i=1;i<=n;i++)
			if(ans[i]){
				printf("%d",i);
				if(--cnt)printf(" ");
			}
		printf("\n");
	}
	return 0;
}
void NTT(int *A,int n,int tp){
	for(int i=1,j=0,k;i<n-1;i++){
		k=n;
		do j^=(k>>=1);while(j<k);
		if(i<j)swap(A[i],A[j]);
	}
	for(int k=2;k<=n;k<<=1){
		int wn=qpow(g,tp>0?(p-1)/k:p-1-(p-1)/k);
		for(int i=0;i<n;i+=k){
			int w=1;
			for(int j=0;j<(k>>1);j++,w=(long long)w*wn%p){
				int a=A[i+j],b=(long long)w*A[i+j+(k>>1)]%p;
				A[i+j]=a+b;
				if(A[i+j]>=p)A[i+j]-=p;
				A[i+j+(k>>1)]=a-b;
				if(A[i+j+(k>>1)]<0)A[i+j+(k>>1)]+=p;
			}
		}
	}
	if(tp<0){
		int inv=qpow(n,p-2);
		for(int i=0;i<n;i++)A[i]=(long long)A[i]*inv%p;
	}
}
int qpow(int a,int b){
	int ans=1;
	for(;b;b>>=1,a=(long long)a*a%p)if(b&1)ans=(long long)ans*a%p;
	return ans;
}

CF993E Nikita and Order Statistics 【fft】

#题目链接 CF993E #题解我们记小于$x$的位置为$1$，否则为$0$ 区间由端点决定，转为两点前缀和相减我们统计出每一种前缀和个数，记为$A[i]$表示值为$i$的位置出现的次数那么对于$k > 0$有 $$ans_k = \sum\limits_{x - y = k} A[x]A[y]$$ 令 $$B[x] = A[n - x]$$ 那么有 $$ans_k = \sum\limits_{x + y = n + k} A[x]B[y]$$ 就成了卷积的形式第$n + k$项系数就是$ans_k \qquad k > 0$

对于$k = 0$，可以直接统计，也可以减去卷积中重复的部分首先减去空串的个数$n + 1$，然后再除以$2$【因为当$x$和$y$相等，大小顺序就可以颠倒了】最后求得的就是$k = 0$的答案

复杂度$O(nlogn)$

#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<vector>
#include<queue>
#include<cmath>
#include<map>
#define LL long long int
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define cls(s,v) memset(s,v,sizeof(s))
#define mp(a,b) make_pair<int,int>(a,b)
#define cp pair<int,int>
using namespace std;
const int maxn = 800005,maxm = 100005,INF = 0x3f3f3f3f;
inline int read(){
	int out = 0,flag = 1; char c = getchar();
	while (c < 48 || c > 57){if (c == ''-'') flag = 0; c = getchar();}
	while (c >= 48 && c <= 57){out = (out << 1) + (out << 3) + c - 48; c = getchar();}
	return flag ? out : -out;
}
struct E{
    double a,b;
    E(){}
    E(double x,double y):a(x),b(y) {}
    E(int x,int y):a(x),b(y) {}
    inline E operator =(const int& b){
        this->a = b; this->b = 0;
        return *this;
    }
    inline E operator =(const double& b){
        this->a = b; this->b = 0;
        return *this;
    }
    inline E operator /=(const double& b){
        this->a /= b; this->b /= b;
        return *this;
    }
};
inline E operator *(const E& a,const E& b){
    return E(a.a * b.a - a.b * b.b,a.a * b.b + a.b * b.a);
}
inline E operator *=(E& a,const E& b){
    return a = E(a.a * b.a - a.b * b.b,a.a * b.b + a.b * b.a);
}
inline E operator +(const E& a,const E& b){
    return E(a.a + b.a,a.b + b.b);
}
inline E operator -(const E& a,const E& b){
    return E(a.a - b.a,a.b - b.b);
}
const double pi = acos(-1);
int R[maxn];
void fft(E* a,int n,int f){
    for (int i = 0; i < n; i++) if (i < R[i]) swap(a[i],a[R[i]]);
    for (int i = 1; i < n; i <<= 1){
        E wn(cos(pi / i),f * sin(pi / i));
        for (int j = 0; j < n; j += (i << 1)){
            E w(1,0),x,y;
            for (int k = 0; k < i; k++,w = w * wn){
                x = a[j + k],y = w * a[j + k + i];
                a[j + k] = x + y; a[j + k + i] = x - y;
            }
        }
    }
    if (f == -1) for (int i = 0; i < n; i++) a[i] /= n;
}
E A[maxn],B[maxn];
int cnt[maxn],N,x,a[maxn];
int main(){
	N = read(); x = read();
	cnt[0]++;
	REP(i,N) cnt[a[i] = a[i - 1] + (read() < x ? 1 : 0)]++;
	for (int i = 0; i <= N; i++) A[i] = cnt[i],B[i] = cnt[N - i];
	int n = 1,L = 0;
	while (n <= (N << 1)) n <<= 1,L++;
	for (int i = 1; i < n; i++) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
	fft(A,n,1); fft(B,n,1);
	for (int i = 0; i < n; i++) A[i] *= B[i];
	fft(A,n,-1);
	A[N].a -= N + 1; A[N].a /= 2;
	for (int i = 0; i <= N; i++)
		printf("%.0lf ",floor(A[N + i].a + 0.3));
	return 0;
}

今天关于使用 Python 的 FFT 和 IFFT 编码热方程解的问题和python fft函数的分享就到这里，希望大家有所收获，若想了解更多关于AtCoder AGC019E Shuffle and Swap (DP、FFT、多项式求逆、多项式快速幂)、BZOJ.4503. 两个串 (FFT/bitset)、CF 827E Rusty String FFT、CF993E Nikita and Order Statistics 【fft】等相关知识，可以在本站进行查询。

本文标签：