多项式全家桶学习笔记（基础）

好的一点一点啃吧

Todo List
  • 任意模数多项式乘法逆
  • 多项式开根（加强版）
  • 多项式快速幂 (加强版)

目录

零.预处理逆元

原来其实没写这个。。。但是多项式反三角卡常必须要用，那就写一下。

记

\[\operatorname{inv}(i)=i^{-1}\pmod P \]

有：

\[\operatorname{inv}(i)=-\left[\frac{P}{i}\right]\operatorname{inv}(P\ \bmod i) \]

Proof

显然

\[P-\left[\frac{P}{i}\right]i\equiv -\left[\frac{P}{i}\right]i\pmod P \]

所以

\[\operatorname{inv}(i)=\frac{1}{i}=-\frac{\left[\frac{P}{i}\right]}{P-\left[\frac{P}{i}\right]i}=-\frac{\left[\frac{P}{i}\right]}{P\ \bmod i}=-\left[\frac{P}{i}\right]\operatorname{inv}(P\ \bmod i) \]

所以可以递推，复杂度\(O(n)\)。

Code 
ll prep(ll n){
    ll lmt=1;
    inv[0]=inv[1]=1;
    while(lmt<n)lmt<<=1;
    for(ll i=2;i<lmt;i++)inv[i]=(P-P/i)*inv[P%i]%P;
    return lmt;
}

一.模数为NTT模数的部分

1. 多项式乘法（FFT）

orz myy！！！

模板

给定序列 \(a_i,b_i\)，下标从 \(0\) 开始，求

\[c_i=\sum_{k=0}^{i}a_kb_{i-k} \]

显然 \(|c|=|a|+|b|-1\)。

首先，我们可以找一个 \(n=2^k\)，使得 \(n>c\)，至于为什么要求 \(n=2^k\) 后面再说。

于是：

\[c_r=\sum_{p,q}[(p+q)\bmod n=r]a_pb_q \]

记 \(\omega\) 为 \(n\) 次本原单位根，则根据单位根反演：

\[[n\bmod k=0]=\frac{1}{n}\sum_{i=0}^{n-1}\omega^{ik} \]

所以

\[\begin{aligned} &[(p+q)\bmod n=r] \cr =&[(p+q-r)\bmod n=0] \cr =&\frac{1}{n}\sum_{k=0}^{n-1}\omega^{k(p+q-r)}\cr =&\frac{1}{n}\sum_{k=0}^{n-1}\omega^{-kr}\omega^{kp}\omega^{kq} \end{aligned} \]

所以

\[\begin{aligned} c_r&=\sum_{p,q}[(p+q)\bmod n=r]a_pb_q\\ &=\sum_{p,q}\frac{1}{n}\sum_{k=0}^{n-1}\omega^{-kr}\omega^{kp}\omega^{kq}a_pb_q\\ &=\frac{1}{n}\sum_{k=0}^{n-1}\omega^{-rk}\sum_{p,q}\omega^{pk}a_p\omega^{qk}a_q\\ &=\frac{1}{n}\sum_{k=0}^{n-1}\omega^{-rk}\sum_p\omega^{pk}a_p\sum_q\omega^{qk}a_q \end{aligned} \]

于是我们发现了这两个式子：

\[g_m=\sum_{k=0}^{n-1}\omega^{mk}f_k \]

\[f_m=\sum_{k=0}^{n-1}\omega^{-mk}g_k \]

根据单位根反演，这两式等价。

根据上面求 \(c\) 的式子，我们发现只要对于每一种关系式，知道右边即可快速求得左边，就可以解决这个问题。

另一方面，我们发现第一个式子等价于

\[g_m=F(\omega^m) \]

其中 \(F\) 为数列 \(f\) 的 OGF。

所以这个操作相当于，给出一个多项式的系数表示，对于所有 \(k\) 求出这个多项式在 \(\omega^k\) 的值，也就是在系数表示和点值表示之间转化。

我们称这样的一个过程为离散傅里叶变换(\(\text{DFT}\))，反过来就是 \(\text{IDFT}\)。

下面我们介绍一种可以在 \(O(n\log n)\) 时间内完成这个操作的算法，叫做快速傅里叶变换 \(\text{FFT}\)。

下面我们只考虑 \(\text{DFT}\)，\(\text{IDFT}\) 同理。

设 \(\omega_n\) 表示 \(n\) 次本原单位根，则有：

\[\omega_{2n}^{2m}=\omega_n^m \]

\[\omega_{2n}^{m}=\omega_{2n}^{m+n} \]

设 \(A_0(x)\) 表示 \(A(x)\) 偶次项的和，\(A_1(x)\) 表示 \(A(x)\) 奇次项的和，有：

\[\begin{aligned} A(\omega_n^m)&=A_0((\omega_n^m)^2)+\omega_n^mA_1((\omega_n^m)^2)\\ &=A_0(\omega_{\frac{n}{2}}^m)+\omega_n^mA_1(\omega_{\frac{n}{2}}^m) \end{aligned} \]

\[\begin{aligned} A(\omega_n^{m+\frac{n}{2}}))&=A_0((\omega_n^m)^2)+\omega_n^{m+\frac{n}{2}}A_1((\omega_n^m)^2)\\ &=A_0(\omega_{\frac{n}{2}}^m)-\omega_n^mA_1(\omega_{\frac{n}{2}}^m) \end{aligned} \]

我们将上述两式称为蝴蝶操作。

于是可以考虑迭代/递归，考虑到递归常数太大，使用迭代法进行处理。

Code 
int lmt,rev[N];
struct C{
    double x,y;
    C(double a=0,double b=0){x=a,y=b;}
};
C operator+(C a,C b){return C(a.x+b.x,a.y+b.y);}
C operator-(C a,C b){return C(a.x-b.x,a.y-b.y);}
C operator*(C a,C b){return C(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
inline void init(int n){
    lmt=1;int t=0;
    while(lmt<n)lmt<<=1,t++;
    for(int i=1;i<lmt;i++)rev[i]=(rev[i>>1]>>1)|(i&1)<<(t-1);
}
inline void FFT(vector<C>&A,int lmt,int tp){
    for(int i=0;i<lmt;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
    for(int mid=1;mid<lmt;mid<<=1){
        C Wn(cos(Pi/mid),tp*sin(Pi/mid));
        for(int R=mid<<1,j=0;j<lmt;j+=R){
            C w(1,0);
            for(int k=0;k<mid;k++,w=w*Wn){
                C x=A[j+k],y=w*A[j+k+mid];
                A[j+k]=x+y,A[j+k+mid]=x-y;
            }
        }
    }
} 
vec mul(vec f,vec g,int len){
    init(len);
    f.resize(lmt),g.resize(lmt);
    vector<C>a(lmt),b(lmt);
    for(int i=0;i<lmt;i++)a[i].x=f[i],b[i].x=g[i],a[i].y=b[i].y=0;
    FFT(a,lmt,1);FFT(b,lmt,1);
    for(int i=0;i<lmt;i++)a[i]=a[i]*b[i];
        FFT(a,lmt,-1);
        vec ret(lmt);
        for(int i=0;i<lmt;i++)ret[i]=lround(a[i].x/lmt);
        return ret;
    }

时间复杂度 \(O(n\log n)\)。

2. 多项式乘法（NTT）

模板同上

注意到 FFT 要用复数的原因是因为他有一个单位根，而我们注意到由于有了模数，所以可以取它的一个原根 \(g\)，并且用 \(g\) 来代替单位根。

具体来说，我们有

\[e^{\frac{i\pi}{mid}}=g^{\frac{P-1}{2mid}} \]

不过这玩意有个限制，就是要求 \(2mid|P-1\) ，否则上式显然不成立。因此，我们一般来说要求 \(V_2(P-1)\ge 20\) ，此处 \(V_2(x)\) 表示 \(x\) 的因子中 \(2\) 的次幂。

通常我们把满足这个要求的模数称为 NTT 模数，典型例子有 \(998244353=1+2^{23}\times 7\times 17\)（\(g=3\)），\(1004535809=1+2^{21}\times 479\)（\(g=3\)）等。

Code 
inline int qpow(int a,int k){
	int ret=1;
	while(k){
		if(k&1)ret=1LL*ret*a%P;
		a=1LL*a*a%P;
		k>>=1;
	}
	return ret%P;
}
inline void init(int n){
	lmt=1;int t=0;
	while(lmt<n)lmt<<=1,t++;
	for(int i=1;i<lmt;i++)rev[i]=(rev[i>>1]>>1)|(i&1)<<(t-1);
}
inline void NTT(vec&A,int lmt,int tp){
	if(A.size()<lmt)A.resize(lmt);
        for(int i=0;i<lmt;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
        for(int m=1,t=0;m<lmt;m<<=1,t++)
            for(int j=0,Wn=pw[t+1];j<lmt;j+=m<<1)
                for(int k=0,w=1,x,y;k<m;k++,w=1LL*w*Wn%P)
                    x=A[j+k],y=1LL*w*A[j+k+m]%P,A[j+k]=(x+y)%P,A[j+k+m]=(x-y+P)%P;
        if(tp==1)return;
        reverse(A.begin()+1,A.begin()+lmt);
        for(int i=0,inv=qpow(lmt,P-2);i<=lmt;i++)A[i]=1LL*A[i]*inv%P;
} 
vec mul(vec f,vec g,int len){
	init(len);
	f.resize(lmt),g.resize(lmt);
	NTT(f,lmt,1);NTT(g,lmt,1);
	for(int i=0;i<lmt;i++)f[i]=1LL*f[i]*g[i]%P;
	NTT(f,lmt,-1);
	return f;
}

时间复杂度还是 \(O(n\log n)\)。

3. 多项式乘法逆

模板

首先设我们已经找到了一个多项式 \(G_1(x)\)，使得

\[G_1(x)F(x)\equiv 1\pmod {x^{\lceil\frac{n}{2}\rceil}} \]

则将其与原式相减即可得到

\[G(x)-G_1(x)\equiv 0\pmod {x^{\lceil\frac{n}{2}\rceil}} \]

平方：

\[G(x)^2-2G(x)G_1(x)+G_1(x)^2\equiv 0\pmod {x^n} \]

两边同乘 \(F(x)\)，根据 \(F(x)G(x)\equiv 1\pmod {x^n}\) 可得：

\[G(x)-2G_1(x)+G_1^2(x)F(x)\equiv 0\pmod {x^n} \]

移项则有

\[G(x)\equiv 2G_1(x)-G_1^2(x)F(x)\pmod {x^n} \]

over。

时间复杂度 \(T(n)=T(\frac{n}{2})+n\log n\)，根据主定理可知 \(T(n)=n\log n\)。

Code 
void getinv(vec&f,vec&g,int len){
	if(f.size()<(len<<2))f.resize(len<<2);
	if(g.size()<(len<<2))g.resize(len<<2);
	if(len==1){g[0]=qpow(f[0],P-2);return;}
	getinv(f,g,len+1>>1);
	init(len<<1);
	vec c(lmt);
	copy(f.begin(),f.begin()+len,c.begin());
	NTT(c,lmt,1),NTT(g,lmt,1);
	for(int i=0;i<lmt;i++)g[i]=1LL*(2LL-1LL*g[i]*c[i]%P+P)%P*g[i]%P;
	NTT(g,lmt,-1);
	for(int i=len;i<lmt;i++)g[i]=0;
}

4. 多项式除法

模板

占位。

4.5 多项式求导，积分

并不想说什么。

Code 
vec getdev(vec f,int len){
	vec g(len-1);
	for(int i=1;i<len;i++)g[i-1]=1LL*i*f[i]%P;
	return g;
}
vec getinvdev(vec f,int len){
	vec g(len+1);
	for(int i=1;i<len;i++)g[i]=1LL*f[i-1]*inv[i]%P;
	return g;
}

复杂度 \(O(n)\)。

5. 多项式 ln

模板

\[G(x)\equiv \ln F(x)\pmod{x^n} \]

两边同时求导

\[G'(x)\equiv \frac{F'(x)}{F(x)}\pmod{x^{n-1}} \]

积回去

\[G(x)\equiv \int \frac{F'(x)}{F(x)}dx\pmod{x^n} \]

over。

复杂度 \(O(n\log n)\).

Code 
vec getln(vec f,int len){
	vec a=getdev(f,len),b;
	getinv(f,b,len);
	return getinvdev(mul(a,b,len<<1),len);
}

5.5 多项式牛顿迭代

若

\[G(F(x))\equiv 0\pmod {x^n} \]

\[G(F_0(x))\equiv 0\pmod {x^{\lceil\frac{n}{2}\rceil}} \]

有

\[F\equiv F_0-\frac{G(F_0)}{G'(F_0)}\pmod{x^n} \]

Proof

根据泰勒展开：

\[G(F)=G(F_0)+G'(F_0)(F-F_0)+G''(F_0)(F-F_0)^2+\cdots \]

因为

\[F\equiv F_0 \pmod {x^{\lceil\frac{n}{2}\rceil}} \]

所以

\[F^2\equiv F_0^2\pmod {x^n} \]

所以

\[0\equiv G(F)\equiv G(F_0)+G'(F_0)(F-F_0) \pmod {x^n} \]

所以

\[F\equiv F_0-\frac{G(F_0)}{G'(F_0)}\pmod{x^n} \]

6 多项式exp

模板

\[G(x)\equiv e^{F(x)}\pmod {x^n} \]

\[\ln G-F\equiv 0 \pmod {x^n} \]

设

\[y=\ln G-F \]

则

\[dy=d(\ln G-F) \]

\[dy=\frac{dG}{G} \]

\[\frac{dy}{dG}=\frac{1}{G} \]

设

\[\ln G_0-F\equiv 0 \pmod {x^{\lceil\frac{n}{2}\rceil}} \]

则根据牛顿迭代公式：

\[G\equiv G_0-\frac{\ln G_0-F}{\frac{1}{G_0}}\equiv G_0(1-\ln G_0+F)\pmod{x^n} \]

于是递归即可。

复杂度 \(T(n)=T(\frac{n}{2})+O(n\log n)\)，根据主定理可知 \(T(n)=O(n\log n)\)。

Code 
void getexp(vec f,vec&g,int len){
	if(g.size()<(len<<2))g.resize(len<<2);
	if(len==1){g[0]=1;return;}
	getexp(f,g,len+1>>1);
	init(len<<1);
	vec d=getln(g,len),e(f.begin(),f.begin()+len);
        NTT(g,lmt,1),NTT(d,lmt,1),NTT(e,lmt,1);
   	for(int i=0;i<lmt;i++)g[i]=1LL*(1-d[i]+e[i]+P)*g[i]%P;
    	NTT(g,lmt,-1);
    	for(int i=len;i<lmt;i++)g[i]=0;
}

7 多项式快速幂

模板

\[G\equiv F^k\pmod {x^n} \]

\[\ln G\equiv k\ln F\pmod {x^n} \]

\[G\equiv e^{k\ln F}\pmod {x^n} \]

over。

复杂度 \(O(n\log n)\)。

Code 
vec getpow(vec f,int len,int k){
	vec g(len),h=getln(f,len);
	for(int i=0;i<len;i++)h[i]=1LL*h[i]*k%P;
	getexp(h,g,len);
	return g;
}

8.多项式开根

模板

\[G\equiv \sqrt F\pmod {x^n} \]

\[\ln G\equiv \frac{F}{2}\pmod{x^n} \]

\[G\equiv e^\frac{F}{2}\pmod{x^n} \]

over。

复杂度 \(O(n\log n)\)。

Code 
vec getsqrt(vec f,int len){
	vec h=getln(f,len),g;
	for(int i=0;i<len;i++)h[i]=1LL*h[i]*inv[2]%P;
	getexp(h,g,len);
	return g;
}

9.多项式三角函数

模板

我们有

\[\begin{cases} e^{ix}=\cos x+i\sin x\\ e^{-ix}=\cos x-i\sin x \end{cases} \]

解得

\[\begin{cases} \cos x=\frac{e^{ix}+e^{-ix}}{2}\\ \sin x=\frac{e^{ix}-e^{-ix}}{2i} \end{cases} \]

又

\[\left(\frac{-1}{P}\right)=(-1)^\frac{P-1}{2}=1 \]

所以存在 \(i_0\in \mathbb Z_P\)，使得 \(i_0^2\equiv -1\pmod P\)。

比如当 \(P=998244353\) 时，\(i=86583718\)。

于是

\[\begin{cases} \cos F=\frac{e^{i_0F}+e^{-i_0F}}{2}\\ \sin F=\frac{e^{i_0F}-e^{-i_0F}}{2i_0} \end{cases} \]

over。

总复杂度 \(O(n\log n)\)。

Q：如果 \(P\) 不是固定的怎么做？

A：暴力枚举 \(i\) 即可，当然也可以用 Cipolla 算法。

Code 
vec sin(vec f,int len){
	vec x(len),X,Y,g(len);
	for(int i=0;i<len;i++)x[i]=1LL*img*f[i]%P;
	getexp(x,X,len),getinv(X,Y,len);
	for(int i=0,inv=qpow(img<<1,P-2);i<len;i++)g[i]=1LL*(X[i]-Y[i]+P)%P*inv%P;
	return g;
}
vec cos(vec f,int len){
	vec x(len),X,Y,g(len);
	for(int i=0;i<len;i++)x[i]=1LL*img*f[i]%P;
	getexp(x,X,len),getinv(X,Y,len);
	for(int i=0;i<len;i++)g[i]=1LL*(X[i]+Y[i])%P*inv[2]%P;
	return g;
}

10.多项式反三角函数

例题

和前面多项式\(\ln\)差不多。。。

首先推一发 \(\arcsin\)：

\[\begin{aligned} y&=\arcsin x\\ \sin y&=x\\ \text d(\sin y)&=\text dx\\ \text dy\times \cos y&=\text dx\\ \frac{dy}{dx}&=\frac{1}{\cos y}\\ \frac{dy}{dx}&=\frac{1}{\sqrt{1-\sin y^2}}\\ \frac{dy}{dx}&=\frac{1}{\sqrt{1-x^2}} \end{aligned} \]

返回原题

\[\begin{aligned} G&\equiv \arcsin F\pmod {x^n}\\ G'&\equiv \frac{F'}{\sqrt{1-F^2}}\pmod {x^{n-1}}\\ G&\equiv \int \frac{F'}{\sqrt{1-F^2}}dx\pmod {x^n} \end{aligned} \]

\(y=\arctan x\) 同理可得：

\[G\equiv \int \frac{F'}{1+F^2}dx\pmod {x^n} \]

over。

总复杂度依然是\(O(n\log n)\)。

Code 
vec arcsin(vec f,int len){
	vec x=getdev(f,len);
	init(len<<1);
	vec y(lmt);
	NTT(f,lmt,1);
	for(int i=0;i<lmt;i++)y[i]=(1+P-1LL*f[i]*f[i]%P)%P;
	NTT(y,lmt,-1);
	for(int i=len;i<lmt;i++)y[i]=0;
	vec z=getsqrt(y,len);
	y.clear();
	getinv(z,y,len);
	NTT(x,lmt,1),NTT(y,lmt,1);
	for(int i=0;i<lmt;i++)x[i]=1LL*x[i]*y[i]%P;
	NTT(x,lmt,-1);
	return getinvdev(x,len);
}
vec arctan(vec f,int len){
	vec x=getdev(f,len);
	init(len<<1);
	vec y(lmt),z;
	NTT(f,lmt,1);
	for(int i=0;i<lmt;i++)y[i]=(1+1LL*f[i]*f[i]%P)%P;
	NTT(y,lmt,-1);
	for(int i=len;i<lmt;i++)y[i]=0;
	getinv(y,z,len);
	NTT(x,lmt,1),NTT(z,lmt,1);
	for(int i=0;i<lmt;i++)x[i]=1LL*x[i]*z[i]%P;
	NTT(x,lmt,-1);
	return getinvdev(x,len);
}

二.模数不为NTT模数的部分

1. 多项式乘法（MTT）

模板

orz myy！！！

1.1 9 次NTT的MTT

考虑到值域最大可以达到 \(10^{23}\)，我们可以取 \(3\) 个NTT模数（如 \(998244353\)，\(1004535809\)，\(167772161\)），用 CRT 合并即可。

共计使用 \(9\) 次NTT。

1.2 7 次FFT的MTT

取 \(M_0=[\sqrt M]\)，设 \(F(x)=M_0A(x)+B(x)\)，\(G(x)=M_0C(x)+D(x)\)，所以 \(F(x)G(x)=M_0^2A(x)C(x)+M_0(A(x)D(x)+B(x)C(x))+B(x)D(x)\)。

注意到中间系数为 \(M_0\) 的两项我们可以直接合并，因此可以省去一次 IDFT，共计 \(7\) 次 DFT。

1.3 DFT 二合一

这里二合一的前置要求是：两个数列都是实数数列。

我们考虑两个实多项式\(A(x),B(x)\)，希望能一次DFT同时求出这两个的DFT。

令

\[P(x)=A(x)+iB(x) \]

\[Q(x)=A(x)-iB(x) \]

为了方便起见，记 \(F_p[k]=P(\omega^k),F_q[k]=Q(\omega^k)\)，并定义 \(X=\frac{2\pi jk}{2L}\)。

所以

\[\begin{aligned} F_p[k]&=A(\omega_{2L}^k)+iB(\omega_{2L}^k)\\ &=\sum_{j=0}^{2L-1}A_j\omega_{2L}^{jk}+iB_j\omega_{2L}^{jk}\\ &=\sum_{j=0}^{2L-1}(A_j+iB_j)(\cos X+i\sin X) \end{aligned} \]

\[\begin{aligned} F_q[k]&=A(\omega_{2L}^k)-iB(\omega_{2L}^k)\\ &=\sum_{j=0}^{2L-1}A_j\omega_{2L}^{jk}-iB_j\omega_{2L}^{jk}\\ &=\sum_{j=0}^{2L-1}(A_j-iB_j)(\cos X+i\sin X)\\ &=\sum_{j=0}^{2L-1}(A_j\cos X+B_j\sin X)+i(A_j\sin X-B_j\cos X)\\ &=\operatorname{conj}\left(\sum_{j=0}^{2L-1}(A_j\cos X+B_j\sin X)-i(A_j\sin X-B_j\cos X)\right)\\ &=\operatorname{conj}\left(\sum_{j=0}^{2L-1}(A_j+iB_j)(\cos(-X)+i\sin(-X))\right)\\ &=\operatorname{conj}\left(\sum_{j=0}^{2L-1}(A_j+iB_j)\omega_{2L}^{-jk}\right)\\ &=\operatorname{conj}\left(\sum_{j=0}^{2L-1}(A_j+iB_j)\omega_{2L}^{(2L-k)j}\right)\\ &=\operatorname{conj}(F_p[2L-k]) \end{aligned} \]

所以只要一次 \(\text{DFT}\) 即可同时算出 \(F_p,F_q\)，从而

\[\text{DFT}(A[k])=\frac{F_p[k]+F_q[k]}{2} \]

\[\text{DFT}(B[k])=i\frac{F_p[k]-F_q[k]}{2} \]

于是一次DFT即可。

1.4 IDFT 二合一

这里二合一的要求是：IDFT的结果必须是实数数列。

考虑前面两个式子，我们知道 \(\text{DFT}(A[k]),\text{DFT}(B[k])\) 之后，就可以得到 \(F_p,F_q\)，就可以求出 \(P(x)\)，于是就可以求出 \(A(x),B(x)\)。

于是一次 IDFT 即可。

1.5 4 次 DFT 的 MTT

考虑上面的7次DFT的方法，总共用了4次DFT与3次IDFT。

于是4次DFT两两合并压缩为2次，3次IDFT两两合并压缩为2次。

Q：为什么不能再压缩了？

A：再压缩就不符合压缩的条件了。

总共使用4次\(\text{DFT}\)。

Code 
vec MTT(vec f,vec g,int n,int m,int P){
	init(n+m);
	const int lim=(1<<15)-1;
	vector<C>a(lmt),b(lmt),c(lmt),d(lmt);
	vec ans(lmt); 
	for(int i=0;i<n;i++)a[i]=C(f[i]&lim,f[i]>>15);
	for(int i=0;i<m;i++)b[i]=C(g[i]&lim,g[i]>>15);
	FFT(a,lmt,1),FFT(b,lmt,1);
	for(int i=0;i<lmt;i++){
		int t=(lmt-i)&(lmt-1);
		c[i]=C((a[i].x+a[t].x)*0.5,(a[i].y-a[t].y)*0.5)*b[i];
		d[i]=C((a[i].y+a[t].y)*0.5,(a[t].x-a[i].x)*0.5)*b[i];
	}
	FFT(c,lmt,-1),FFT(d,lmt,-1);
	for(int i=0;i<lmt;i++)c[i]=c[i]/lmt,d[i]=d[i]/lmt;
	for(int i=0;i<lmt;i++){
		long long p=c[i].x+0.5,o=c[i].y+0.5,x=d[i].x+0.5,u=d[i].y+0.5;
		ans[i]=1LL*(p%P+((o+x)%P<<15)+(u%P<<30))%P;
	}
	return ans;
}
posted @ 2020-08-02 21:18  happydef  阅读(875)  评论(5)  编辑  收藏  举报