CF1103E Radix sum【DFT,扩域】
题目描述：给出一个长度为 \(n\) 的序列 \(a_i\)，对于每一个 \(p\in [0,n)\cap \Z\)，求有多少个长为 \(n\) 的整数序列 \(i_j\) 满足：\(i_j\in [1,n]\)，且在十进制不进位加法下 \(\sum_{j=1}^na_{i_j}=p\)。答案对 \(2^{58}\) 取模。
数据范围：\(n\le 10^5,0\le x_i<10^5\)
实际上就相当于对于一个幂级数 \(F\)，在 \(5\) 维长度为 \(10\) 的循环卷积意义下，求 \(F^n\)。
于是考虑如何做 dft 和 idft。这个跟平常的 ntt 差别实际上在于：
- 模 \(2^{58}\) 意义下除以 \(10^5\)？
- \(\omega_{10}\) 怎么办？
首先看第一个问题，我们把模数改为 \(2^{64}\)，也就是 ULL 自然溢出。\(5^5\) 就可以算逆元，\(2^5\) 显然就可以直接除掉。
关于第二个问题，众所周知 \(\omega_{10}\) 是无理数，这个模数还这么丑，于是就只能暴力扩域。设 \(x=\omega_5\)，\(\omega_{10}=-x^3\)，则可以将用到的数都表示为 \(x\) 的多项式。
大家都知道余数定理：设 \(F,G\) 为两个多项式，\(G(a)=0\)，则 \((F\bmod G)(a)=F(a)\)。我们取 \(G\) 为五次分圆多项式 \(1+x+x^2+x^3+x^4\)。于是在系数模 \(2^{64}\) ，整体模 \(G\) 的意义下，得到的三次多项式将 \(\omega_5\) 代入就得到答案。
剩下的问题在于，如何求出最后的答案。我们知道答案都为整数，然而我们得到的是 \(a_0+a_1\omega_5+a_2\omega_5^2+a_3\omega_5^3\)。因为它是整数，所以可以设它为 \(a_0-y\)，则 \(\omega_5\) 为 \(y+a_1x+a_2x^2+a_3x^3=0\) 的根。根据分圆多项式在有理数域上的不可约性，\(G\) 是被 \((x-\omega_5)\) 整除的，次数最小的，首1的有理多项式，所以 \(G|(y+a_1x+a_2x^2+a_3x^3)\)，所以 \(a_1=a_2=a_3=y=0\)，答案即为 \(a_0\)。
#include<bits/stdc++.h>
#define Rint register int
using namespace std;
typedef unsigned long long LL;
const int N = 100000, po10[5] = {1, 10, 100, 1000, 10000};
const LL inv = 6723469279985657373ull, mod = (1ull << 58) - 1;
template<typename T>
inline void read(T &x){
int ch = getchar(); x = 0; bool f = false;
for(;ch < '0' || ch > '9';ch = getchar()) f |= ch == '-';
for(;ch >= '0' && ch <= '9';ch = getchar()) x = x * 10 + ch - '0';
if(f) x = -x;
}
struct comp {
LL a[4];
comp(){memset(a, 0, sizeof a);}
comp operator = (const comp &o){memcpy(a, o.a, sizeof a); return *this;}
comp operator * (LL val) const {
comp res;
for(Rint i = 0;i < 4;++ i) res.a[i] = a[i] * val;
return res;
}
comp operator + (const comp &o) const {
comp res;
for(Rint i = 0;i < 4;++ i) res.a[i] = a[i] + o.a[i];
return res;
}
comp operator - (const comp &o) const {
comp res;
for(Rint i = 0;i < 4;++ i) res.a[i] = a[i] - o.a[i];
return res;
}
comp operator * (const comp &o) const {
comp res; static LL x[7]; memset(x, 0, sizeof x);
for(Rint i = 0;i < 4;++ i)
for(Rint j = 0;j < 4;++ j) x[i + j] += a[i] * o.a[j];
for(Rint i = 2;~i;-- i){
for(Rint j = 0;j < 4;++ j) x[i + j] -= x[i + 4];
x[i + 4] = 0;
}
res.a[0] = x[0]; res.a[1] = x[1]; res.a[2] = x[2]; res.a[3] = x[3];
return res;
}
} w10[11], A[N];
int n;
inline comp ksm(comp a, int b){
comp res; res.a[0] = 1;
for(;b;b >>= 1, a = a * a) if(b & 1) res = res * a;
return res;
}
void dft(comp *A){
static comp B[10];
for(Rint d = 0;d < 5;++ d)
for(Rint i = 0;i < N;++ i) if(!((i / po10[d]) % 10)){
memset(B, 0, sizeof B);
for(Rint j = 0;j < 10;++ j)
for(Rint k = 0;k < 10;++ k) B[j] = B[j] + A[i + k * po10[d]] * w10[j * k % 10];
for(Rint j = 0;j < 10;++ j) A[i + j * po10[d]] = B[j];
}
}
void idft(comp *A){
static comp B[10];
for(Rint d = 0;d < 5;++ d)
for(Rint i = 0;i < N;++ i) if(!((i / po10[d]) % 10)){
memset(B, 0, sizeof B);
for(Rint j = 0;j < 10;++ j)
for(Rint k = 0;k < 10;++ k) B[j] = B[j] + A[i + k * po10[d]] * w10[10 - j * k % 10];
for(Rint j = 0;j < 10;++ j) A[i + j * po10[d]] = B[j];
}
for(Rint i = 0;i < N;++ i) A[i] = A[i] * inv;
}
int main(){
read(n); w10[0].a[0] = 1; w10[1].a[3] = -1;
for(Rint i = 2;i <= 10;++ i) w10[i] = w10[i - 1] * w10[1];
for(Rint i = 1, x;i <= n;++ i){read(x); ++ A[x].a[0];}
dft(A); for(Rint i = 0;i < N;++ i) A[i] = ksm(A[i], n); idft(A);
for(Rint i = 0;i < n;++ i) cout << ((A[i].a[0] >> 5) & mod) << endl;
}
