BZOJ3557: [Ctsc2014]随机数
BZOJ上数据范围没写清楚，实际数据范围如这个表格所示。
把$M_n$和$x$看成$\mathbb{F}_2$下的$m$维的列向量，则有转移矩阵 $$A = \begin{pmatrix} 0 & 0 & \cdots & 0 & -x_0\\ 1 & 0 & \cdots & 0 & -x_1\\ 0 & 1 & \cdots & 0 & -x_2\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & -x_{m-1}\\ \end{pmatrix}$$ 该矩阵为多项式 $$p(z) = z^m+\sum_{i=0}^{m-1}x_iz^i$$ 的友阵。因此，若把$M_n$看成$F$中的元素，则有 $$M_n(z) = z^nM_0(z)$$ 其中$F = \mathbb{F}_2[z]/(p(z))$为模$p(z)$意义下的$\mathbb{F}_2$上的多项式构成的环。当然，也可以直接通过观察得出这个结论。注意到$n < m$时$F$中的$z^n$与$\mathbb{F}_2[z]$中的相同，于是可以以复杂度 $$O\bigl(m\log m\max(1,\log k-\log m)\bigr)$$ 解决第一种情况。
所有$2^m-1$个非零元素都能被生成，是数列在这些元素中等概率取值的必要条件。该条件成立时有 $$\begin{aligned} 2^m-1 &= |\mathopen{<}z\mathclose{>}M_0(z)| = |\mathopen{<}z\mathclose{>}|\\ &\le |F^\times| \le |F \setminus \{0\}| = 2^m-1 \end{aligned}$$ 其中$F^\times$为$F$中所有有乘法逆元的元素构成的乘法群。由此得$\mathopen{<}z\mathclose{>} = F^\times = F \setminus \{0\}$，那么$F$为大小为$2^m$的有限域，且$\forall f(z) \in F^\times$，有$(f(z))^{2^m} = f(z)$，故 $$\begin{aligned} M_k(z) &= (z^{k \cdot 2^l})^{2^{m-l}}M_0(z)\\ &= \Bigl(M_{k \cdot 2^l}(z)\bigl(M_0(z)\bigr)^{-1}\Bigr)^{2^{m-l}}M_0(z)\\ &= \bigl(M_{k \cdot 2^l}(z)\bigr)^{2^{m-l}}\Bigl(\bigl(M_0(z)\bigr)^{2^{m-l}}\Bigr)^{2^l-1}\\ \end{aligned}$$ 从而可以以复杂度$O(m^2\log m)$解决第二种情况。
由上述推导可知，$x$是好的等价于$F$是一个域，且$z$是$F^\times$的生成元。$F$是一个域当且仅当$p(z)$是$\mathbb{F}_2$上的既约多项式。由此，或许可以用某些方法生成第二种情况的数据。
这题挺好写的，就一多项式除法就没了，不知道为啥没几个人写。
#include<algorithm>
#include<vector>
#include<cstdio>
#define RAN(a)a.begin(),a.end()
using namespace std;
typedef unsigned long long u64;
typedef unsigned u32;
namespace num{
const u32 p=1811939329;
const u32 g=13;
inline u32 imod(u32 a){
return a<p?a:a-p;
}
inline u32 ipow(u32 a,u32 n){
u32 s=1;
for(;n;n>>=1){
if(n&1)
s=(u64)s*a%p;
a=(u64)a*a%p;
}
return s;
}
inline u32 inv(u32 a){
return ipow(a,p-2);
}
}
using namespace num;
class poly_base{
public:
vector<u32>a;
poly_base(){}
explicit poly_base(int n):a(n){}
void fft(int n,bool f){
a.resize(n);
if(n<=1)
return;
for(int i=0,j=0;i<n;++i){
if(i<j)
std::swap(a[i],a[j]);
int k=n>>1;
while((j^=k)<k)
k>>=1;
}
vector<u32>w(n/2);
w[0]=1;
for(int i=1;i<n;i<<=1){
for(int j=i/2-1;~j;--j)
w[j<<1]=w[j];
int m=(p-1)/i/2;
u64 s=ipow(g,f?p-1-m:m);
for(int j=1;j<i;j+=2)
w[j]=s*w[j-1]%p;
for(int j=0;j<n;j+=i<<1){
u32*b=&a[0]+j,*c=b+i;
for(int k=0;k<i;++k){
u32 v=(u64)w[k]*c[k]%p;
c[k]=imod(b[k]+p-v);
b[k]=imod(b[k]+v);
}
}
}
}
void dft(int n){
fft(n,0);
}
void idft(){
int n=a.size();
fft(n,1);
u64 f=inv(n);
for(int i=0;i<n;++i)
a[i]=f*a[i]%p;
}
};
class poly:public poly_base{
public:
poly(){}
explicit poly(int n):poly_base(n){}
u32&operator[](int i){
return a[i];
}
const u32&operator[](int i)const{
return a[i];
}
int size()const{
return a.size();
}
void swap(poly&b){
a.swap(b.a);
}
static int len(int n){
while(n^n&-n)
n+=n&-n;
return n;
}
void fix(){
while(size()&&!a.back())
a.pop_back();
}
void mul(poly b){
fix();
b.fix();
int n=len(size()+b.size()-1);
dft(n);
b.dft(n);
for(int i=0;i<n;++i)
a[i]=(u64)a[i]*b[i]%p;
idft();
for(int i=0;i<n;++i)
a[i]&=1;
fix();
}
void sqr(){
fix();
int n=len(2*size()-1);
dft(n);
for(int i=0;i<n;++i)
a[i]=(u64)a[i]*a[i]%p;
idft();
for(int i=0;i<n;++i)
a[i]&=1;
fix();
}
void mod(int n){
a.resize(n);
}
void inv(int n){
int m=len(n);
mod(m);
vector<u32>b(1,1);
a.swap(b);
for(int i=2;i<=m;i<<=1){
poly c(i);
for(int j=0;j<i;++j)
c[j]=b[j];
sqr();
mul(c);
mod(i);
}
mod(n);
}
void div(poly b){
fix();
b.fix();
int n=size()-b.size()+1;
if(n<=0){
a.clear();
return;
}
reverse(RAN(a));
mod(n);
reverse(RAN(b.a));
b.inv(n);
mul(b);
mod(n);
reverse(RAN(a));
fix();
}
void mod(poly b){
fix();
b.fix();
int m=b.size();
if(size()>=m){
poly c=*this;
div(b);
mul(b);
mod(m-1);
for(int i=0;i<m-1;++i)
a[i]^=c[i];
fix();
}
}
};
struct ano{
operator u64(){
u64 x=0;
signed c=getchar();
while(c<48)
c=getchar();
while(c>47)
x=x*10+c-48,c=getchar();
return x;
}
}buf;
poly basis(int n){
poly a(n+1);
a[n]=1;
return a;
}
int main(){
int n=buf;
poly p(n+1);
p[n]=1;
for(int i=0;i<n;++i)
p[i]=buf;
poly s(n);
for(int i=0;i<n;++i)
s[i]=buf;
if(buf==0){
u64 k=buf;
if(k<n){
s.mul(basis(k));
s.mod(p);
}else{
int i=0;
while((k&(2<<i)-1)<n)
++i;
s.mul(basis(k&(1<<i)-1));
s.mod(p);
poly a=basis(1<<i);
k>>=i;
for(;k;k>>=1){
if(k&1){
s.mul(a);
s.mod(p);
}
if(k>1){
a.sqr();
a.mod(p);
}
}
}
}else{
int l=buf%n;
poly a(n);
for(int i=0;i<n;++i)
a[i]=buf;
s.swap(a);
for(int i=0;i<n-l;++i){
s.sqr();
s.mod(p);
a.sqr();
a.mod(p);
}
for(int i=0;i<l;++i){
s.mul(a);
s.mod(p);
if(i<l-1){
a.sqr();
a.mod(p);
}
}
}
s.mod(n);
for(int i=0;i<n;++i)
printf("%d",s[i]);
puts("");
}
Rewrited @ 姫野星奏生誕祭2019
Edited @ 2019-06-28
