最近正好月考,然而并不想去月考,于是来颓 blog 吧。
 
给你一个 $n-1$ 次多项式 $F(x)$ 和一个 $m-1$ 次多项式 $G(x)$,让你求 $(F\times G)(x)$ 的各项系数。$n,m\leqslant 10^{6}$。
我们不妨将 $F(x)$ 的 $i$ 次项系数记为 $F[i]$
不难发现直接暴力算是 $O(n^{2})$ 的,因此我们需要优化。
不过为了优化,我们得先扯远点。
点值表达 如果我们选取 $n$ 个点 $(x_{0},y_{0}),(x_{1},y_{1}),\cdots,(x_{n-1},y_{n-1})$,并且其中 $x_{i}$ 两两不同,我们就可以唯一地确定出一个 $n-1$ 次多项式。
就比如说
瞎写的(
我们可以列出如下的三元一次方程组:
解得
注意到 $(F\times G)(x)=F(x)G(x)$(废话),我们只要知道了 $F(x)$ 和 $G(x)$ 的点值表达,就可以 $O(n)$ 的计算出 $(F\times G)(x)$ 的点值表达了。因为 $(F\times G)(x)$ 的次数是 $n+m-2$,我们选出前 $n+m-1$ 个自然数即可。
但是这还不够。注意到,如果我们选取的点的横坐标如果很普通,我们首先需要 $O(n^{2})$ 转成点值表达,乘完之后又要转回系数表达,效率甚至不如暴力(
于是我们还需要优化,于是我们还需要再扯远一点。
复数 对,你没看错,扯到复数了。
基本概念 根据初中学习的数学知识,我们知道有些一元二次方程是没有实数根的。就比如说
我们知道,它的判别式是 $\Delta=b^{2}-4ac=-4\lt 0$,因此它没有实数根。
于是我们就定义了虚数单位 $i$,并规定 $i^{2}=-1$。形如 $x+yi$ 这样的数被称为复数。它的模长被定义为它到原点的距离,即 $\sqrt{x^{2}+y^{2}}$;幅角被定义为与 $x$ 轴正半轴的夹角,即 $\operatorname{tan}^{-1}\cfrac{y}{x}$。
因为 $i$ 不是实数,它不能被画在我们现有的数轴上。那我们就再拿来一条数轴,将两条数轴垂直放置,垂足为原点。或者说,我们可以将这理解成平面直角坐标系,$x$ 就是实轴,$y$ 轴就是虚轴,$x+yi$ 就对应了点 $(x,y)$。
上图展示了 $4+i$ 和 $3+4i$。
复数的运算其实没什么出乎意料的:
我们来单独考虑一下复数相乘的几何意义。
假设我们有两个复数 $c_{1}$,$c_{2}$,它们的模长分别是 $r_{1}$,$r_{2}$,幅角分别是 $\alpha_{1}$,$\alpha_{2}$。不难发现我们有
总结成一句话,就是「模长相乘,幅角相加」。
单位圆与单位根 单位圆就是指半径为 $1$ 的圆。不过一般我们都是把它画在原点的。
考虑这么一个方程
它的所有复数根。
因为复数相乘意味着模长相乘,如果一个复数的 $n$ 次方等于 $1$,它自身的模长也应该是 $1$。如果它的幅角是 $a$,我们应该有
不难想象出我们有 $n$ 个这样的复数,它们的幅角通式是 $\cfrac{2k\pi}{n}(k\in[0,n)\cap\mathbb{Z})$。我们称其中幅角等于 $\cfrac{2\pi}{n}$ 的复数,即 $\cos\cfrac{2\pi}{n}+i\sin\cfrac{2\pi}{n}$ 为 $n$ 次单位根 $\omega_{n}$,我们就可以把这 $n$ 个复数表示为 $\omega_{n}^{k}(k\in[0,n)\cap\mathbb{Z})$。
它有如下的一些性质
因为它们的模长都是 $1$,幅角相等就相等了。
然后就没了,读者自证不难(
因为 $\omega_{n}^{\frac{n}{2}}$ 的幅角是 $\cfrac{2\times\frac{n}{2}\pi}{n}=\pi$,不难发现它就是 $-1$ 了。
运行过程 那么问题来了,我们刚才扯了这么些,有什么用呢?
就是说,有个叫「让 · 巴普蒂斯 · 约瑟夫 · 傅里叶」(Jean Baptiste Joseph Fourier)的神仙有一天大开脑洞,掏出 $n$ 次单位根想要求点值表达。
首先,我们通过在高次补 $0$ 的方式,将这个多项式的项数(也就是次数 $+1$)补到 $2$ 的非负整数次幂。
然后一巴掌把这个多项式拍成两半,按奇偶性拼成两个多项式:
于是我们有
现在我们假设 $0\leqslant k<\cfrac{n}{2}$,将 $x=\omega_{n}^{k}$ 和 $x=\omega_{n}^{k+\frac{n}{2}}$ 代入
我们注意到,如果说我们求出了 $F_{1}(x)$ 和 $F_{2}(x)$ 的点值表达,我们就能够 $O(n)$ 地求出 $F(x)$ 的点值表达了。至此,我们不难想到分治,边界条件就是 $n=1$,这时什么也不用做直接返回即可。
简单地贴一下代码:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 void  work (complex * F,int  n)  {     if (n==1 )          return ;      complex  F1[n>>1 ],F2[n>>1 ];  for (int  i=0 ;i<n;i+=2 ){  	 F1[i]=F[2 *i];  F2[i]=F[2 *i+1 ];      }  work(F1,n>>1 );  work(F2,n>>1 );  complex unit(cos(2*Pi/n),sin(2*Pi/n)),tmp=1;  for (int  i=0 ;i<(n>>1 );++i,tmp=tmp*unit){  	 complex  w=F2[i+(n>>1 )]*tmp;  F[i]=F1[i]+w;  F[i+(n>>1 )]=F1[i]-w;      } } 
 
不过,以上代码是我现写出来的,保证其不正确性不保证其正确性。
不过,这样递归运算效率还是太低了。我们来考虑一下,递归到底层后,F 数组的每一个位置上实际存的是几次项:
1 2 3 4 0 1 2 3 4 5 6 7 0 2 4 6|1 3 5 7 0 4|2 6|1 5|3 7 0|4|2|6|1|5|3|7 
 
写成二进制看看:
位置(十进制) 
位置(二进制) 
次数(十进制) 
次数(二进制) 
 
 
$0$ 
$000$ 
$0$ 
$000$ 
 
$1$ 
$001$ 
$4$ 
$100$ 
 
$2$ 
$010$ 
$2$ 
$010$ 
 
$3$ 
$011$ 
$6$ 
$110$ 
 
$4$ 
$100$ 
$1$ 
$001$ 
 
$5$ 
$101$ 
$5$ 
$101$ 
 
$6$ 
$110$ 
$3$ 
$011$ 
 
$7$ 
$111$ 
$7$ 
$111$ 
 
 
 
注意到两边的数字的二进制是镜像的。也就是说,我们只要把位置编号的二进制位的最后 $\operatorname{log}n$ 位左右镜像过来,就可以得到递归到最底层后,这个位置上的系数所对应的项的次数了。
我们可以 $O(n)$ 地处理处每个数的镜像:
1 2 for (int  i=0 ;i<n;++i)     rev[i]=(rev[i>>1 ]>>1 )|((i&1 )?(n>>1 ):0 ); 
 
中间有一个位或运算符,我们以它为分界线,将上面这行代码分成左右两部分。左边就是 $i$ 这个数除了最后一位以外的所有位的镜像;右边特判了一下 $i$ 的最后一位是否为 $1$,如果是的话,就在最高位补一个 $1$。
接下来的步骤我不是很能解释得清楚,因为我也是背的板子。总之这个东西写出来差不多长这样:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 void  work (complex * F,int  n)  {     for (int  i=0 ;i<n;++i)          if (i<rev[i])              swap(F[i],F[rev[i]]);      for (int  len=2 ,p=1 ;len<=n;len<<=1 ,p<<=1 ){          complex  unit (cos (Pi/p),sin (Pi/p))  ;          for (int  i=0 ;i<n;i+=len){              complex  cur (1 ,0 )  ;              for (int  j=i;j<i+p;++j){                  complex  tmp=F[j+p]*cur;                  F[j+p]=F[j]-tmp;                  F[j]=F[j]+tmp;                  cur=cur*unit;              }          }      } } 
 
那么现在还差最后一步,将点值表达转回系数表达。我们将 $F(\omega_{n}^{0}),F(\omega_{n}^{1}),\cdots,F(\omega_{n}^{n-1})$ 分别记为 $y_{0},y_{1},\cdots,y_{n-1}$,我们求点值的过程可以用矩阵表达成这样:
我们可以试图寻找左边那个矩阵的逆矩阵。直接摆结论的话,它差不多长这样:
现在我们要证明这两个矩阵乘起来是单位矩阵。令第一个矩阵为 $A$,第二个矩阵为 $B$,$A\times B=C$。不难发现 $A[i][j]=\omega_{n}^{ij}$,$B[i][j]=\cfrac{1}{n}\omega_{n}^{-ij}$,我们有
若 $i=j$,不难发现此时 $C[i][j]=1$。
否则,即 $i\neq j$,设 $i-j=l$,我们有
综上所述,矩阵 $C$ 是单位矩阵,因此矩阵 $B$ 是矩阵 $A$ 的逆矩阵。
因此我们有
注意到这个过程和我们之前将系数转化为点值表达的过程极为相似。这就是在启示我们,如果说我们把之前用的 $\omega_{n}^{0},\omega_{n}^{1},\cdots,\omega_{n}^{n-1}$ 换成 $\omega_{n}^{0},\omega_{n}^{-1},\cdots,\omega_{n}^{-n+1}$,然后对着点值表达来一遍 FFT,然后再除以 $n$,就得到了原多项式的系数了。
注意到这两个过程只有用到的单位根不一样,我们可以将上面那段代码的
1 complex  unit (cos (Pi/p),sin (Pi/p))  ;
 
改成
1 complex  unit (cos (Pi/p),tp*sin (Pi/p))  ;
 
然后调用时再传一个参数 tp 进去。tp=1 表示是系数转点值,tp=-1 表示是点值转系数。
完整代码如下:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 #include <cstdio>  #include <cmath>  #define  maxn 2097152 const  double  Pi=acos (-1.0 );int  n,m;int  rev[maxn+1 ];struct  complex {     double  r,c;      complex (double  _r=0 ,double  _c=0 ){          r=_r;          c=_c;      } }; complex  F[maxn+1 ],G[maxn+1 ];inline  complex  operator +(complex & a,complex & b){     return  complex (a.r+b.r,a.c+b.c); } inline  complex  operator -(complex & a,complex & b){     return  complex (a.r-b.r,a.c-b.c); } inline  complex  operator *(complex & a,complex & b){     return  complex (a.r*b.r-a.c*b.c,a.r*b.c+a.c*b.r); } template <typename  _tp>inline  void  swap (_tp& x,_tp& y)  {    _tp tmp=x;      x=y;      y=tmp; } inline  void  FAQ (complex  F[],int  tp)  {     for (register  int  i=0 ;i<n;++i)          if (i<rev[i])              swap(F[i],F[rev[i]]);      for (register  int  len=2 ,p=1 ;len<=n;len<<=1 ,p<<=1 ){          register  complex  unit (cos (Pi/p),tp*sin (Pi/p))  ;          for (register  int  i=0 ;i<n;i+=len){              register  complex  cur (1 ,0 )  ;              for (register  int  j=i;j<i+p;++j){                  register  complex  tmp=F[j+p]*cur;                  F[j+p]=F[j]-tmp;                  F[j]=F[j]+tmp;                  cur=cur*unit;              }          }      } } int  main ()  {     scanf ("%d%d" ,&n,&m);      ++n;      for (register  int  i=0 ;i<n;++i)          scanf ("%d" ,&F[i].r);      ++m;      for (register  int  i=0 ;i<m;++i)          scanf ("%d" ,&G[i].r);      for (m+=n-1 ,n=1 ;n<m;n<<=1 );      for (register  int  i=0 ;i<n;++i)          rev[i]=(rev[i>>1 ]>>1 )|((i&1 )?(n>>1 ):0 );      FAQ(F,1 );      FAQ(G,1 );      for (register  int  i=0 ;i<n;++i)          F[i]=F[i]*G[i];      FAQ(F,-1 );      for (register  int  i=0 ;i<m;++i)          printf ("%0.0lf " ,F[i].r/n);      return  0 ; } 
 
不要问我 cltstream 哪去了(
然后您就可以切掉这道板子题  了。
其实上面的代码是我从另一道题里复制过来然后现改的,说不定会改出错(
最好还是自己写吧(
什么?想知道是哪道题?往下看(
例题 这是一道裸的卷积题。
考虑两个长度为 $n$ 的数组 $F$ 和 $G$,现在我们想求一个数组 $H$,它满足
在本页面往上翻,翻到这个式子:
发现了吗?这两个过程其实是一样的。
于是,我们如下构造两个多项式
然后直接一波 FFT 套上去,输出次数最低的 $n$ 项的系数就好。
那么这题呢?首先我们把 $j\lt i$ 和 $j\gt i$ 分开计算。
先考虑 $j\lt i$,令 $F[i]=q_{i}$,$G[i]=\begin{cases}&\cfrac{1}{i^{2}}\;\;&(i\gt 0)\\&0&(i=0)\end{cases}$,那么
直接套板子就行。
对于 $j>i$ 的情况,我们将数组 $F$ 左右翻转,然后继续套板子就行(
具体还是看代码吧:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 #include \lt cstdio> #include <cmath>  #define  maxn 2097152 const  double  Pi=acos (-1.0 );int  n,m;int  rev[maxn+1 ];struct  complex {     double  r,c;      complex (double  _r=0 ,double  _c=0 ){          r=_r;          c=_c;      } }; complex  F1[maxn+1 ],F2[maxn+1 ],G[maxn+1 ];inline  complex  operator +(complex & a,complex & b){     return  complex (a.r+b.r,a.c+b.c); } inline  complex  operator -(complex & a,complex & b){     return  complex (a.r-b.r,a.c-b.c); } inline  complex  operator *(complex & a,complex & b){     return  complex (a.r*b.r-a.c*b.c,a.r*b.c+a.c*b.r); } template <typename  _tp>inline  void  swap (_tp& x,_tp& y)  {    _tp tmp=x;      x=y;      y=tmp; } inline  void  FAQ (complex  F[],int  tp)  {     for (register  int  i=0 ;i<n;++i)          if (i<rev[i])              swap(F[i],F[rev[i]]);      for (register  int  len=2 ,p=1 ;len<=n;len<<=1 ,p<<=1 ){          register  complex  unit (cos (Pi/p),tp*sin (Pi/p))  ;          for (register  int  i=0 ;i<n;i+=len){              register  complex  cur (1 ,0 )  ;              for (register  int  j=i;j<i+p;++j){                  register  complex  tmp=F[j+p]*cur;                  F[j+p]=F[j]-tmp;                  F[j]=F[j]+tmp;                  cur=cur*unit;              }          }      } } int  main ()  {     scanf ("%d" ,&n);      for (register  int  i=0 ;i<n;++i){          scanf ("%lf" ,&F1[i].r);          F2[n-i-1 ].r=F1[i].r;      }      for (register  int  i=1 ;i<n;++i)          G[i]=1.0 /i/i;      for (m=n,n=1 ;n<m;n<<=1 );      n<<=1 ;      for (register  int  i=0 ;i<n;++i)          rev[i]=(rev[i>>1 ]>>1 )|((i&1 )?(n>>1 ):0 );      FAQ(G,1 );      FAQ(F1,1 );      for (register  int  i=0 ;i<n;++i)          F1[i]=F1[i]*G[i];      FAQ(F1,-1 );      FAQ(F2,1 );      for (register  int  i=0 ;i<n;++i)          F2[i]=F2[i]*G[i];      FAQ(F2,-1 );      for (register  int  i=0 ;i<m;++i)          printf ("%lf\n" ,(F1[i].r-F2[m-i-1 ].r)/n);      return  0 ; } 
 
注意到朴素的 FFT 使用的是单位复根。然而它们有一个十分大的缺陷,就是必须要用 double 存。这会带来精度上的误差,一个直接的结果就是,对于只有整数参与的多项式乘法,跑完 FFT 却会出现小数。
这就启示我们,能不能用其他的什么东西替换掉单位复根。
设有两个互质的正整数 $a$ 和 $p$,最小的使得 $a^{k}\equiv 1\pmod{p}$ 的 $k$ 被称为 $a$ 模 $p$ 的阶,记作 $\delta_{p}(a)$。
如果说 $\delta_{p}(a)=\varphi(p)$,我们就称 $a$ 是模 $p$ 的一个原根。
现在我们找一个质数 $p=an+1$,其中 $a$ 是一个正整数,$n$ 是 $2$ 的非负整数次幂。然后我们找到它的原根 $g$,并定义 $\omega_{n}=g^{a}$。让我们来看看单位复根有的性质现在的原根有没有:
$\omega_{n}^{0},\omega_{n}^{1},\cdots,\omega_{n}^{n-1}$ 互不相同。这是为了保证点值表达的合法。
虽然我不会证,不过我们的确有 $g^{0},g^{1},\cdots,g^{n-1}$ 在模 $p$ 意义下互不相同,$a$ 次幂自然也一样。
Updated on 2019-03-17
似乎这么证明有些问题?
原根的这个性质实际上是 $g^{0},g^{1},\cdots,g^{p-2}$ 在模 $p$ 意义下互不相同。因此如果 $(n-1)a<p-1$,那么上面的结论就是正确的。
以上。
 
$\omega_{2n}^{2k}=\omega_{n}^{k}$。这是为了让我们可以分治。
根据定义,$\omega_{2n}=g^{\frac{a}{2}}$,就是将现在的 $p$ 进一步拆成 $\cfrac{a}{2}\cdot2n+1$。不难发现
因而原根有上述性质。
$\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}$,或者说 $\omega_{n}^{\frac{n}{2}}=-1$。这同样是为了让我们可以分治。
因为 $p=an+1$,根据费马小定理,我们有
因而 $\omega_{n}^{\frac{n}{2}}\equiv\pm 1\pmod{p}$。又因为 $\omega_{n}^{0}=1$,而 $\omega_{n}^{\frac{n}{2}}\not\equiv\omega_{n}^{0}\pmod{p}$,我们就得到 $\omega_{n}^{\frac{n}{2}}\equiv -1\pmod{p}$。
若 $k\neq 0$,$\sum_{i=0}^{n-1}(\omega_{n}^{k})^{i}=0$。这是为了实现逆变换。不过这个很明显,就留作习题吧。
以上,我们成功地用原根取代了单位复根。一般情况下,我们会取 $p=998244353=7\times 17\times 2^{23}+1$,它的原根是 $3$。
需要注意的是,我们还有一个可以优化的小细节。注意到 $\omega_{n}=g^{a}=g^{\frac{p-1}{n}}$,我们可以将这些值预处理出来,就不用每次都跑快速幂了。对于 $998244353$ 来说,我们需要预处理 $\omega_{2},\omega_{4},\cdots,\omega_{2^{23}}$。注意到
我们一遍快速幂算出 $\omega_{2^{23}}$ 然后倒着推出剩下的即可。
关于 $\omega_{n}^{-k}$,注意到它就是 $\omega_{n}^{k}$ 在模 $998244353$ 意义下的逆元。那么我们把上面两个式子中的 $3$ 换成它在模 $998244353$ 意义下的逆元即可,这个数是 $332748118$。
然后我们把 FFT 板子里的单位复根全部换成原根,运算换成模意义下的就行了。
代码:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 #include <cstdio>  #define  re register #define  maxn 2097152 #define  mod 998244353 #define  swap(a,b) a^=b,b^=a,a^=b namespace  cltstream{     #define  size 1048576       char  cltin[size+1 ],*ih=cltin,*it=cltin;      inline  char  gc ()  {          #ifdef  ONLINE_JUDGE               if (ih==it){                  it=(ih=cltin)+fread(cltin,1 ,size,stdin );                  if (ih==it)                      return  EOF;              }              return  *ih++;          #else                return  getchar();          #endif        }      char  cltout[size+1 ],*oh=cltout,*ot=cltout+size;      inline  void  pc (char  c)  {          if (oh==ot){              fwrite(cltout,1 ,size,stdout );              oh=cltout;          }         *oh++=c;      }      #define  clop() fwrite(cltstream::cltout,1,cltstream::oh-cltstream::cltout,stdout)       #undef  size       template  <typename  _tp>      inline  void  read (_tp& x)  {          int  sn=1 ;          char  c=gc();          for (;c!=45 &&(c<48 ||c>57 )&&c!=EOF;c=gc());          if (c==45 &&c!=EOF)              sn=-1 ,c=gc();          for (x=0 ;c>=48 &&c<=57 &&c!=EOF;x=(x<<3 )+(x<<1 )+(c^48 ),c=gc());          x*=sn;      }      template  <typename  _tp>      inline  void  write (_tp x,char  text=-1 )  {          if (x<0 )              pc(45 ),x=-x;          if (!x)              pc(48 );          else {              int  digit[22 ];              for (digit[0 ]=0 ;x;digit[++digit[0 ]]=x%10 ,x/=10 );              for (;digit[0 ];pc(digit[digit[0 ]--]^48 ));          }          if (text>=0 )              pc(text);      } } int  n,m;int  unit[2 ][24 ],rev[maxn+1 ],F[maxn+1 ],G[maxn+1 ];inline  int  cltpow (re int  x,re int  y)  {     re int  res=1 ;      for (;y;){          if (y&1 )              res=1L L*res*x%mod;          x=1L L*x*x%mod;          y>>=1 ;      }      return  res; } inline  void  NTT (int * F,re int  tp)  {     for (re int  i=0 ;i<n;++i)          if (i<rev[i])              swap(F[i],F[rev[i]]);      for (re int  k=1 ,p=1 ;p<n;++k,p<<=1 )          for (re int  i=0 ;i<n;i+=p<<1 )              for (re int  j=i,tmp=1 ;j<i+p;++j,tmp=1L L*tmp*unit[tp][k]%mod){                  re int  x=F[j],y=1L L*F[j+p]*tmp%mod;                  F[j]=(x+y)%mod;                  F[j+p]=(x-y+mod)%mod;              } } int  main ()  {     unit[0 ][23 ]=cltpow(3 ,119 );      unit[1 ][23 ]=cltpow(332748118 ,119 );      for (re int  i=0 ;i<2 ;++i)          for (re int  j=22 ;j>=0 ;--j)              unit[i][j]=1L L*unit[i][j+1 ]*unit[i][j+1 ]%mod;      cltstream::read(n);      cltstream::read(m);      ++n;      for (re int  i=0 ;i<n;++i)          cltstream::read(F[i]);      ++m;      for (re int  i=0 ;i<m;++i)          cltstream::read(G[i]);      for (m+=n-1 ,n=1 ;n<m;n<<=1 );      for (re int  i=0 ;i<n;++i)          rev[i]=(rev[i>>1 ]>>1 )|((i&1 )?(n>>1 ):0 );      NTT(F,0 );      NTT(G,0 );      for (re int  i=0 ;i<n;++i)          F[i]=1L L*F[i]*G[i]%mod;      NTT(F,1 );      n=cltpow(n,mod-2 );      for (re int  i=0 ;i<m;++i)          cltstream::write(1L L*F[i]*n%mod,i<m-1 ?32 :-1 );      clop();      return  0 ; } 
 
以上。