当前位置:网站首页>FFT模板
FFT模板
2022-08-10 07:56:00 【ThXe】
#include<bits/stdc++.h>
typedef long long ll;
typedef unsigned long long ull;
const int N = 5000007;
const double PI = acos(-1);
inline int read()
{
register int x = 0, f = 1;
register 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;
}
struct Complex
{
double x, y;
Complex(double x = 0, double y = 0) : x(x), y(y) {
}
}a[N], b[N];
Complex operator * (Complex J, Complex Q) {
//模长相乘,幅度相加
return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);
}
Complex operator - (Complex J, Complex Q) {
return Complex(J.x - Q.x, J.y - Q.y);
}
Complex operator + (Complex J, Complex Q) {
return Complex(J.x + Q.x, J.y + Q.y);
}
struct fft {
int n, m;int res, ans[N];
int limit = 1;
int L;//二进制的位数
int R[N];
void FFT(Complex* A, int type)
{
for (int i = 0; i < limit; ++i)
if (i < R[i])
swap(A[i], A[R[i]]);
//i小于R[i]时才交换,防止同一个元素交换两次,回到它原来的位置。
//从底层往上合并
for (int mid = 1; mid < limit; mid <<= 1) {
//待合并区间长度的一半,最开始是两个长度为1的序列合并,mid = 1;
Complex wn(cos(PI / mid), type * sin(PI / mid));//单位根w_n^1;
for (int len = mid << 1, pos = 0; pos < limit; pos += len) {
//len是区间的长度,pos是当前的位置,也就是合并到了哪一位
Complex w(1, 0);//幂,一直乘,得到平方,三次方...
for (int k = 0; k < mid; ++k, w = w * wn) {
//只扫左半部分,蝴蝶变换得到右半部分的答案,w 为 w_n^k
Complex x = A[pos + k];//左半部分
Complex y = w * A[pos + mid + k];//右半部分
A[pos + k] = x + y;//左边加
A[pos + mid + k] = x - y;//右边减
}
}
}
if (type == 1) return;
for (int i = 0; i <= limit; ++i)
A[i].x /= limit;
//最后要除以limit也就是补成了2的整数幂的那个N,将点值转换为系数
//(前面推过了点值与系数之间相除是N)
}
void init()
{
n = read(), m = read();
//读入多项式的每一项,保存在复数的实部
for (int i = 0; i <= n; ++i)
a[i].x = read();
for (int i = 0; i <= m; ++i)
b[i].x = read();
while (limit <= n + m)
limit <<= 1, L++;
//也可以写成:limit = 1 << int(log2(n + m) + 1);
// 补成2的整次幂,也就是N
for (int i = 0; i < limit; ++i)
R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
}
void output()
{
FFT(a, 1);//FFT 把a的系数表示转化为点值表示
FFT(b, 1);//FFT 把b的系数表示转化为点值表示
//计算两个系数表示法的多项式相乘后的点值表示
for (int i = 0; i <= limit; ++i)
a[i] = a[i] * b[i];
//对应项相乘,O(n)得到点值表示的多项式的解C,利用逆变换完成插值得到答案C的点值表示
FFT(a, -1);
for (int i = 0; i <= n + m; ++i)
//这里的 x 和 y 是 double 的 hhh
printf("%d ", (int)(a[i].x + 0.5));//注意要+0.5,否则精度会有问题
}
}FFT;
int main()
{
FFT.init();
FFT.output();
}
边栏推荐
猜你喜欢
【MySQL】SQL语句
添加spark的相关依赖和打包插件(第六弹)
PLSQL学习第四天
协同工具满足70%-90%的工作需求,成为企业香饽饽
ATH10 sensor reads temperature and humidity
MySQL索引事务
Add spark related dependencies and packaging plugins (sixth bullet)
Uni-app develops WeChat applet using local images as background images
解决win10win7win8系统找不到指定的模块,注册不了大漠插件的问题
快速输入当前日期与时间
随机推荐
本地生活商家如何通过短视频赛道,提升销量曝光量?
神经网络的三种训练方法,神经网络训练全过程
Rust learning: 6.1_Slices of composite types
navicat for mysql 连接时报错:1251-Client does not support authentication protocol requested by server
pytest之parametrize参数化
Synchronization lock synchronized traces the source
90. (cesium house) cesium height monitoring events
Quickly enter the current date and time
30条实用MySQL优化法则
Rust学习:6.3_复合类型之元组
DGIOT三千万电表集抄压测
Power function Exponential function Logarithmic function
34. Talk about why you want to split the database?What methods are there?
神经网络样本太少怎么办,神经网络训练样本太少
.NET-8. My Thought Notes
VS2013-debug assembly code-generate asm file-structure memory layout-function parameter stack-calling convention
数据库公共字段自动填充
[In-depth study of 4G/5G/6G topic-56]: L3 signaling control-5-radio bearer management
【Rust指南】使用Cargo工具高效创建Rust项目 | 理解Rust特别的输入输出语句
winget包管理器