- 给定两个正整数 (d,p)。
- 有一个数组 (a_{1sim 5000})，其中 (a_1=x,a_2=y)（(x,y) 未知），(a_{3sim 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)。
- (2le dle 10)，(3le ple10^9+9) 且 (p) 为质数
乘积转化
(xy) 有一个经典的转化方式就是 (xy=frac{(x+y)^2-x^2-y^2}{2})。
由于已有加法运算，我们只需要实现乘常数运算、减法运算、平方运算即可。
基本运算
先考虑实现乘常数运算 * i v k 表示令 (a_k=(a_i imes v) operatorname{mod} p)（(v) 为某个常数）。这就是经典的快速乘，过程中只涉及加法操作。
然后考虑实现减法运算 - i j k 表示令 (a_k=(a_i-a_j) operatorname{mod} p)。就相当于令 (a_k=(a_i+a_j imes (p-1)) operatorname{mod} p)。
注意前面乘常数运算的部分有个小问题，就是若 (v=0) 我们需要令 (a_k=0)，因此得先预处理出一个位置值为 (0)。这只要任选一个位置自减即可，自减过程虽然会用到乘常数运算，但此时常数为 (p-1) 肯定不为 (0)，所以不会出问题。
平方运算
仅通过线性运算难以得出平方级别的东西，所以我们必须考虑利用 (d) 次幂运算。
考虑如果利用二项式定理将 ((x+i)^d) 拆开，其中包含常数项以及 (x) 的 (1sim 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_isum_{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_{0sim 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
",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
",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
",M(D(Sqr(A(1,2)),A(Sqr(1),Sqr(2))),X+1>>1)),0;//xy=((x+y)^2-x^2-y^2)/2
}