【CF1060H】Sophisticated Device（构造）

题目链接

• 给定两个正整数 \(d,p\)。
• 有一个数组 \(a_{1\sim 5000}\)，其中 \(a_1=x,a_2=y\)（\(x,y\) 未知），\(a_{3\sim 5000}=1\)。
• 你可以给出三种指令：+ i j k 表示令 \(a_k=(a_i+a_j)\ \operatorname{mod}\ p\)；^ i k 表示令 \(a_k=a_i^d\ \operatorname{mod}\ p\)；f i 表示输出 \(a_i\)。
• 要求执行给出的指令后，对于任意 \(x,y\)，输出的数都恰好是 \(xy\ \operatorname{mod}\ p\)。
• \(2\le d\le 10\)，\(3\le p\le10^9+9\) 且 \(p\) 为质数

乘积转化

\(xy\) 有一个经典的转化方式就是 \(xy=\frac{(x+y)^2-x^2-y^2}{2}\)。

由于已有加法运算，我们只需要实现乘常数运算、减法运算、平方运算即可。

基本运算

先考虑实现乘常数运算 * i v k 表示令 \(a_k=(a_i\times v)\ \operatorname{mod}\ p\)（\(v\) 为某个常数）。这就是经典的快速乘，过程中只涉及加法操作。

然后考虑实现减法运算 - i j k 表示令 \(a_k=(a_i-a_j)\ \operatorname{mod}\ p\)。就相当于令 \(a_k=(a_i+a_j\times (p-1))\ \operatorname{mod}\ p\)。

注意前面乘常数运算的部分有个小问题，就是若 \(v=0\) 我们需要令 \(a_k=0\)，因此得先预处理出一个位置值为 \(0\)。这只要任选一个位置自减即可，自减过程虽然会用到乘常数运算，但此时常数为 \(p-1\) 肯定不为 \(0\)，所以不会出问题。

平方运算

仅通过线性运算难以得出平方级别的东西，所以我们必须考虑利用 \(d\) 次幂运算。

考虑如果利用二项式定理将 \((x+i)^d\) 拆开，其中包含常数项以及 \(x\) 的 \(1\sim d\) 次项。

如果我们将 \(x^d,(x+1)^d,\cdots(x+d)^d\) 配上系数 \(v\) 相加，就能得到方程：

\[x^2=\sum_{i=0}^dv_i(x+i)^d=\sum_{i=0}^dv_i\sum_{j=0}^dC_d^jx^ji^{d-j}=\sum_{j=0}^d(\sum_{i=0}^dC_d^ji^{d-j}v_i)x^j \]

分别考虑每种次数的项，得到 \(d+1\) 个方程，其中第 \(j\) 个方程可以表示为：

\[\sum_{i=0}^dC_d^ji^{d-j}\cdot v_i=[j=2] \]

\(d+1\) 个方程，\(d+1\) 个未知数，高斯消元即可解出 \(v_{0\sim d}\)。

然后就可以通过代入 \(\sum_{i=0}^dv_i(x+i)^d\) 求出 \(x^2\)，实现平方运算。

代码：\(O(d^3)\)

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Rg register
#define RI Rg int
#define Cn const
#define CI Cn int&
#define I inline
#define W while
#define N 10
using namespace std;
int o,d,X,tot=2,C[N+5][N+5];
I int A(CI x,CI y) {return ++tot,printf("+ %d %d %d\n",x,y,tot),tot;}//加法运算
I int M(CI x,RI v) {RI u=x,t=o?o:0;W(v) v&1&&(t=t?A(t,u):u),u=A(u,u),v>>=1;return t;}//乘常数运算
I int D(CI x,CI y) {RI z=M(y,X-1);return A(x,z);}//减法运算
namespace G
{
	I int QP(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
	int a[N+5][N+5],v[N+5];I void Gauss()//高斯消元
	{
		RI i,j,k,t;for(i=0;i<=d;++i)
		{
			if(!a[i][i]) {for(j=i;j<=d&&!a[j][i];++j);for(swap(v[i],v[j]),k=i;k<=d;++k) swap(a[i][k],a[j][k]);}
			for(j=i+1;j<=d;++j) for(t=1LL*(X-a[j][i])*QP(a[i][i],X-2)%X,v[j]=(v[j]+1LL*t*v[i])%X,k=i;k<=d;++k) a[j][k]=(a[j][k]+1LL*t*a[i][k])%X;
		}
		for(i=d;~i;--i) for(v[i]=1LL*v[i]*QP(a[i][i],X-2)%X,j=i-1;~j;--j) v[j]=(v[j]-1LL*v[i]*a[j][i]%X+X)%X;
	}
}
I int Sqr(CI x) {RI t;for(RI i=0,p,s;i<=d;++i) p=i?A(p,5000):x,++tot,printf("^ %d %d\n",p,tot),s=M(tot,G::v[i]),t=i?A(t,s):s;return t;}//代入求平方
int main()
{
	RI i,j;scanf("%d%d",&d,&X),o=D(1,1);//令a[o]=0
	for(C[0][0]=i=1;i<=d;++i) for(C[i][0]=j=1;j<=i;++j) C[i][j]=(C[i-1][j-1]+C[i-1][j])%X;//预处理组合数
	RI p;for(i=0;i<=d;++i) for(j=0;j<=d;++j) G::a[j][i]=1LL*C[d][j]*G::QP(i,d-j)%X;G::v[2]=1,G::Gauss();//列出方程
	return printf("f %d\n",M(D(Sqr(A(1,2)),A(Sqr(1),Sqr(2))),X+1>>1)),0;//xy=((x+y)^2-x^2-y^2)/2
}