题解 P3803 【【模板】多项式乘法(FFT)】
小菜鸟
2019-11-05 20:26:48
这题...我可能是唯一一个用三种语言各写一遍的人
~~然而Pascal卡不过去~~
并且我应该是第一个用C写FFT通过的(
所以这篇题解的重点不在数学,而在如何用C愉快地食用FFT(
~~既然前面有P党卡过去发题解那这个应该也行吧~~
---
先重复一下做法:把两个多项式变换成点值表示,直接$O(n)$乘起来,然后变回系数表示就是了。。。
---
然后进正题
由于C++有面向对象、函数重载等一系列特性,复数运算很容易以优美的方式呈现
然而C既不能重载函数,也难以面向对象...
然而在研究`std::complex`的时候,我们发现它的头文件`include`了一个`complex.h`库...
于是上网学习,发现C也有复数库(其实甚至不需要库,在C99中复数已经是内建类型了)
于是就极度舒适(
---
FFT需要的基本用法:
```c
#include<complex.h>//为了方便而包含,实际可以省略,但声明变量需将complex改为_Complex
double complex a;//声明基于double的复数变量
creal(a)//取a的实部
cimag(a)//取a的虚部
a=1+2*_Complex_I//_Complex_I即是单位虚数i
//另外不担心重名可以用I代替_Complex_I
//然后不考虑可移植性的话,在GCC下提供1.0fi,1.0i,1.0li等写法
//同时从C++14开始,类似的1.0if,1.0i,1.0il是标准的虚数字面量
```
然后大量的C数学函数均有针对复数特化的泛型函数,在名字前写前缀`c`就是了
于是有了神奇的复数,我们就可以愉快地搞FFT了(
---
然后是不开IO优化卡不过去的代码
```cpp
#include<stdio.h>
#include<ctype.h>
#include<math.h>
#include<complex.h>
char gc(void)
{
static char buf[1<<16],*p1=buf,*p2=buf;
if(p1==p2)
{
p2=(p1=buf)+fread(buf,1,1<<16,stdin);
if(p1==p2)return EOF;
}
return *p1++;
}
// #define gc getchar
void read(int *x)
{
*x=0;
char c=gc();
while(!isdigit(c))c=gc();
while(isdigit(c))
{
*x=*x*10+c-48;
c=gc();
}
}
void write(int x)
{
if(x>=10)write(x/10);
putchar(x%10+48);
}
#define N (1<<21)
#define PI 3.141592653589793238462643383
int n,m;
double complex a[N],b[N];
int rev[N];
void getrev(int n)
{
for(int i=0;i<n;++i)
{
rev[i]=(rev[i>>1]>>1)|(i&1?n>>1:0);
}
}
void fft(double complex *a,int n,int dft)
{
for(int i=0;i<n;++i)
{
if(i<rev[i])
{
double complex temp=a[i];
a[i]=a[rev[i]];
a[rev[i]]=temp;
}
}
for(int step=1;step<n;step<<=1)
{
const double complex wn=cos(PI/step)+sin(dft*PI/step)*_Complex_I;
for(int j=0;j<n;j+=step<<1)
{
double complex wnk=1;
for(int k=j;k<j+step;++k)
{
double complex x=a[k];
double complex y=a[k+step]*wnk;
a[k]=x+y;
a[k+step]=x-y;
wnk*=wn;
}
}
}
if(dft==-1)
{
double inv=1./n;
for(int i=0;i<n;++i)a[i]*=inv;
}
}
int main(void)
{
read(&n),read(&m);
for(int i=0;i<=n;++i)
{
int x;
read(&x);
a[i]=x;
}
for(int i=0;i<=m;++i)
{
int x;
read(&x);
b[i]=x;
}
int len=1;
while(len<=n+m+1)len<<=1;
getrev(len);
fft(a,len,1);
fft(b,len,1);
for(int i=0;i<len;++i)a[i]*=b[i];
fft(a,len,-1);
for(int i=0;i<n+m+1;++i)
{
write((int)(creal(a[i])+0.5));
putchar(' ');
}
}
```
~~awsl为什么这题突然减小时限卡常数啊~~
警告,若编译器有`__STDC_NO_COMPLEX__`宏,则不支持复数(gcc可以放心
最后,虽然我是个C++&面向对象的死忠粉,但还是要说一句:
## C牛逼!!!
另外为什么C的语法高亮萎掉了啊求修复(被迫写`cpp`QwQ
---
### 2020.10.14更新
初赛前颓废学了一下Rust
考虑到FFT代码相对简单于是就来练练手
---
简要介绍Rust,一个安全性和效率都相当不错的新兴语言
由Mozilla的大佬创造维护
定位和C++类似,零开销,多范式
但具有更多安全的、流行的特性
基本上可以把UB扼杀在编译期
~~而且拥有比GNU STL清晰一些的标准库~~
---
这篇题解本来就不是讲思路而是讲代码的...
所以下面直接给出代码实现(PS下面有效率对比)
```rust
const P:i64=998244353;
const g:i64=3;
const gi:i64=332748118;
fn power(mut a:i64,mut b:i64)->i64
{
let mut res:i64=1;
while b!=0
{
if (b&1)!=0
{
res=res*a%P;
}
a=a*a%P;
b>>=1;
}
return res;
}
fn getrev(rev:&mut Vec<usize>,len:usize)
//Rust对于所有权的控制很严格,只能用可变引用来改变外部的东西
//Rust的引用与C++的指针类似
{
rev.clear();
rev.push(0);
for i in 1..len
{
rev.push((rev[i>>1]>>1)|(if (i&1)==0{0}else{len>>1}));
}
}
fn ntt(arr:&mut Vec<i64>,n:usize,dft:bool,rev:&Vec<usize>)
{
for i in 0..n
{
if i<rev[i]
{
arr.swap(i,rev[i]);//在Rust中,一个对象不能同时存在两个可变引用
//所以直接std::mem::swap(i,rev[i])是unsafe的,需要用Vec自己的方法
}
}
let mut step=1;
while step<n
{
let wn:i64=power(if dft{g}else{gi},(P-1)/(step as i64*2));
let mut j=0;
while(j<n)
{
let mut wnk=1;
for k in j..j+step
{
let x=arr[k];
let y=arr[k+step]*wnk%P;
arr[k]=(x+y)%P;
arr[k+step]=(x-y+P)%P;
wnk=wnk*wn%P;
}
j+=step*2;
}
step*=2;
}
if !dft
{
let inv=power(n as i64,P-2);
for i in 0..n
{
arr[i]=arr[i]*inv%P;
}
}
}
extern "C"
{
fn getchar()->u8;
}
fn read()->i64//rust没有原生的比较简洁的快读,不得不学习巨佬cosmicAC的
{
let mut res=0;
let mut ch;
loop
{
unsafe
{
ch=getchar();
}
if ch>=48 && ch<=57
{
break;
}
}
loop
{
res=res*10+(ch as i64)-48;
unsafe
{
ch=getchar();
}
if ch<48 || ch>57
{
break;
}
}
return res;
}
fn main()
{
let n:usize=read() as usize;//Rust在safe mode里不能定义全局变量
let m:usize=read() as usize;
let n=n+1;
let m=m+1;
let mut a:Vec<i64>=Vec::new();
let mut b:Vec<i64>=Vec::new();
let mut len=1;
while len<=n+m
{
len*=2;
}
for i in 0..n
{
let x=read();
a.push(x);
}
while a.len()<len
{
a.push(0);
}
for i in 0..m
{
let x=read();
b.push(x);
}
while b.len()<len
{
b.push(0);
}
let mut rev:Vec<usize>=Vec::new();
for i in 0..len
{
rev.push(0);
}
getrev(&mut rev,len);
ntt(&mut a,len,true,&rev);
ntt(&mut b,len,true,&rev);
let mut res:Vec<i64>=Vec::new();
for i in 0..len
{
res.push(a[i]*b[i]%P);
}
ntt(&mut res,len,false,&rev);
for i in 0..n+m-1
{
print!("{} ",res[i]);
}
}
```
C++在使用普通快读的情况下是1.4s(有没有`fread`区别不大)
Rust是1.8s...
看来GCC开发这么多年毕竟是比rustc成熟的233