hdu 5730 Shell Necklace——多项式求逆+拆系数FFT

时间：2022-11-13 15:30:58

题目：http://acm.hdu.edu.cn/showproblem.php?pid=5730

可以用分治FFT。但自己只写了多项式求逆。

和COGS2259几乎很像。设A(x)，指数是长度，系数是方案。 \( A(x)^{k} \) 的 m 次项系数表示 k 个连续段组成长度为 m 的序列的方案数。

\( B(x)=1+F(x)+F^{2}(x)+F^{3}(x)+... \)

\( B(x) = \frac{1}{1-F(x)} \)（通过计算B(x)的逆来看出这个式子）

然后多项式求逆就行了。

注意模数 \( 313=2^{3}*3*13 \) ，原根是10，但那个 23 太小了！如果 len 大于3的话就会除出小数，所以不能直接用NTT！

那么就用FFT。FFT不能中途取模，所以最大的值是 312×312×10000=9734400000，会让FFT的精度变得很低。所以用拆系数FFT。

#include<iostream>

#include<cstdio>

#include<cstring>

#include<algorithm>

#include<cmath>

#define db double

#define ll long long

using namespace std;

const int N=1e5+,M=(<<)+,mod=;

const db pi=acos(-);

int n,a[M],b[M],tp[M],len,r[M],base;

struct cpl{db x,y;}A[M],B[M],Ta[M],Tb[M],Tc[M],Td[M],Ini,I;

cpl operator+ (cpl a,cpl b){return (cpl){a.x+b.x,a.y+b.y};}

cpl operator- (cpl a,cpl b){return (cpl){a.x-b.x,a.y-b.y};}

cpl operator* (cpl a,cpl b){return (cpl){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}

cpl cnj(cpl a){return (cpl){a.x,-a.y};}

int rdn()

{

  int ret=;bool fx=;char ch=getchar();

  while(ch>''||ch<''){if(ch=='-')fx=;ch=getchar();}

  while(ch>=''&&ch<='') ret=ret*+ch-'',ch=getchar();

  return fx?ret:-ret;

}

void upd(int &x){x>=mod?x-=mod:;}

int pw(int x,int k)

{int ret=;while(k){if(k&)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=;}return ret;}

void fft(cpl *a,bool fx)

{

  for(int i=;i<len;i++)

    if(i<r[i])swap(a[i],a[r[i]]);

  for(int R=;R<=len;R<<=)

    {

      int m=R>>;

      cpl wn=(cpl){ cos(pi/m),fx?-sin(pi/m):sin(pi/m) };

      for(int i=;i<len;i+=R)

    {

      cpl w=I;

      for(int j=;j<m;j++,w=w*wn)

        {

          cpl x=a[i+j], y=w*a[i+m+j];

          a[i+j]=x+y;  a[i+m+j]=x-y;

        }

    }

    }

  if(!fx)return;

  for(int i=;i<len;i++)a[i].x/=len,a[i].y/=len;

}

void mtt(int n,int *a,int *b,int *c)

{

  for(len=;len<n<<;len<<=);

  for(int i=;i<len;i++)r[i]=(r[i>>]>>)+((i&)?len>>:);

  for(int i=;i<n;i++)A[i]=(cpl){ a[i]/base,a[i]%base }; for(int i=n;i<len;i++)A[i]=Ini;

  for(int i=;i<n;i++)B[i]=(cpl){ b[i]/base,b[i]%base }; for(int i=n;i<len;i++)B[i]=Ini;

  fft(A,);  fft(B,);

  cpl ta,tb,tc,td;

  A[len]=A[]; B[len]=B[];

  for(int i=,j=len;i<len;i++,j--)

    {

      ta=(A[i]+cnj(A[j]))*(cpl){0.5,};

      tb=(A[i]-cnj(A[j]))*(cpl){,-0.5};

      tc=(B[i]+cnj(B[j]))*(cpl){0.5,};

      td=(B[i]-cnj(B[j]))*(cpl){,-0.5};

      Ta[i]=ta*tc; Tb[i]=ta*td; Tc[i]=tb*tc; Td[i]=tb*td;

    }

  A[len]=B[len]=Ini;

  for(int i=;i<len;i++)A[i]=Ta[i]+Tb[i]*(cpl){,};

  for(int i=;i<len;i++)B[i]=Tc[i]+Td[i]*(cpl){,};

  fft(A,);  fft(B,);

  for(int i=,Da,Db,Dc,Dd;i<n;i++)

    {

      Da=(ll)(A[i].x+0.5)%mod;  Db=(ll)(A[i].y+0.5)%mod;

      Dc=(ll)(B[i].x+0.5)%mod;  Dd=(ll)(B[i].y+0.5)%mod;

      c[i]=(Da*base*base+(Db+Dc)*base+Dd)%mod+mod; upd(c[i]);

    }

}

void getinv(int n,int *a,int *b)

{

  if(n==){b[]=pw(a[],mod-);return;}

  getinv(n+>>,a,b);

  mtt(n,a,b,tp);

  mtt(n,tp,b,tp);

  for(int i=;i<n;i++)b[i]=((b[i]<<)-tp[i])%mod+mod,upd(b[i]);

}

int main()

{

  base=sqrt(mod); I.x=;

  while()

    {

      memset(a,,sizeof a);memset(b,,sizeof b);

      n=rdn(); if(!n)return ;

      for(int i=;i<=n;i++)a[i]=rdn();

      for(int i=;i<=n;i++)a[i]=mod-a[i]%mod,upd(a[i]);

      a[]++;

      getinv(n+,a,b);

      printf("%d\n",b[n]);

    }

  return ;

}

相关文章