魈凯KBS • 1个月前
typedef unsigned char byte; typedef unsigned int word; typedef unsigned long long ull; typedef long long ll; const word mod=1004535809,size=21; //NTT模数,结果不超过(1<<size) char io[(1<<size)+1];//输入输出 word a[1<<size],b[1<<size];//多项式/高精 word realid[1<<size],root[1<<size],inverse[1<<size]; //迭代后系数所在的位置,单位根及其逆元 word i,id,floor; ll num1,num2; char *top; //循环变量 inline ll pow(ll a,ll b){//快速幂 register ll ans=1; for(;b;b>>=1){
if(b&1) (ans*=a)%=mod;
(a*=a)%=mod;
}
return ans;
} inline void loading(){//预处理迭代后位置,单位根及其逆元
root[0]=inverse[0]=1;
num1=pow(3,mod>>size);
num2=pow(num1,mod-2);
for(i=1;i<1<<size;i++){
root[i]=num1*root[i-1]%mod;
inverse[i]=num2*inverse[i-1]%mod;
for(id=i,floor=0;floor<size;floor++,id>>=1)
realid[i]=realid[i]<<1|(id&1);
}
} inline void read(){ //直接用fread,然后指针前移读入数据,效率较高,写法较简单 //因fread特性,DEV-CPP下不能以数字结尾(不停读入最后一个字符),但linux下可以
top=io+(1<<size);
fread(io+1,1,1<<size,stdin);
while('0'>*top||*top>'9') top--;
for(i=0;'0'<=*top&&*top<='9';top--,i++)
a[realid[i]]=*top-'0';//直接放到迭代后的位置
while('0'>*top||*top>'9') top--;
for(i=0;'0'<=*top&&*top<='9';top--,i++)
b[realid[i]]=*top-'0';//b同理
} inline void DFT(){ //非迭代版DFT
for(floor=0;floor<size;floor++)
for(i=0;i<1<<size;i+=(1<<(floor+1)))
for(id=0;id<1<<floor;id++){
num1=a[i+id];//蝴蝶变换
num2=a[i+id+(1<<floor)];
(num2*=root[id<<size-floor-1])%=mod;
a[i+id]=(num1+num2)%mod;//放回原位
a[i+id+(1<<floor)]=(num1+mod-num2)%mod;
num1=b[i+id];//b同理
num2=b[i+id+(1<<floor)];
(num2*=root[id<<size-floor-1])%=mod;
b[i+id]=(num1+num2)%mod;
b[i+id+(1<<floor)]=(num1+mod-num2)%mod;
}
} inline void IDFT(){
for(floor=0;floor<21;floor++)//与DFT相同
for(i=0;i<0x200000;i+=1<<floor+1)
for(id=0;id<1<<floor;id++){
num1=root[i+id];
num2=root[i+id+(1<<floor)];
(num2*=inverse[id<<20-floor])%=mod;//乘上单位根逆元
root[i+id]=(num1+num2)%mod;
root[i+id+(1<<floor)]=(num1+mod-num2)%mod;
}
} inline void write(){//fwrite输出,同输入
num2=pow(1<<size,mod-2);
top=io+(1<<size);
num1=0;
for(i=0;i<1<<size;i++){
num1+=num2*root[i]%mod;//最后要乘上n的逆元
top--;
*top=num1%10+'0';
num1/=10;
}
while(*top=='0'&&top!=io+(1<<size)) top++;
if(top==io+(1<<size)) putchar('0');
else fwrite(top,1,io+(1<<size)-top,stdout);
} int main(){
loading();
read();
DFT();
for(i=0;i<1<<size;i++)
root[realid[i]]=(ll)(a[i])*b[i]%mod;//点值相乘
IDFT();
write();
return 0;
}
Comments: