题意：

有$g(n)=\sum_{i=0}^nf(i)^2$，其中$f(0)=0,f(1)=1,f(n)=f(n−2)+2f(n−1)(n≥2)$，给出n，y，x，s，求$x^{g(n*y)}mod(s+1)$。

思路：

有递推式：

$$f(n)^2=f(n-2)^2+4f(n-1)^2+4f(n-2)f(n-1)$$ $$f(n)f(n-1)=f(n-2)f(n-1)+2f(n-1)^2$$ $$g(n)=g(n-1)+f(n)^2$$

因此可以构造矩阵：

$$\left[ \begin{matrix} f(n)^2 \\ f(n-1)^2 \\ f(n)f(n-1) \\ g(n) \\ \end{matrix} \right] = \left[ \begin{matrix} 4&1&4&0 \\ 1&0&0&0 \\ 2&0&1&0 \\ 4&1&4&1 \end{matrix} \right] \left[ \begin{matrix} f(n-1)^2 \\ f(n-2)^2 \\ f(n-1)f(n-2) \\ g(n-1) \\ \end{matrix} \right]\\ 可得 \left[ \begin{matrix} f(ny)^2 \\ f(ny-1)^2 \\ f(ny)f(ny-1) \\ g(ny) \\ \end{matrix} \right] = \left[ \begin{matrix} 4&1&4&0 \\ 1&0&0&0 \\ 2&0&1&0 \\ 4&1&4&1 \end{matrix} \right]^{ny-1} \left[ \begin{matrix} f(1)^2 \\ f(0)^2 \\ f(1)f(0) \\ g(1) \\ \end{matrix} \right]$$

但是这个结果是放在指数上的，这里就要用到欧拉降幂公式了。

$$A^BmodC=A^{Bmodφ(C)+φ(C)}modC$$

还有另一种写法：可以推出$g(n)=\frac{f(n)f(n+1)}{2}$，

这里除以2不能直接用逆元处理，因为并没有保证2和s+1互质，所以可以用到公式$\frac{a}{b}modc=\frac{amod(b*c)}{b}$（b和c不一定要互质）。

还有其他公式$\frac{a}{b}modc=a*b^{φ(c)-1}modc$。

代码：

struct mat
{
    long long a[4][4];
}m,r;
ll MOD;
ll euler(ll n)
{
    ll i,res=n;
    for(i=2;i*i<=n;i++)
    {
        if(n%i==0)
        {
            res=res/i*(i-1);
            while(n%i==0)
                n=n/i;
        }
    }
    if(n>1)res=res/n*(n-1);
    return res;
}
mat mult(mat x,mat y)
{
    mat res={0};
    for(int i=0;i<4;i++)
        for(int j=0;j<4;j++)
            for(int k=0;k<4;k++)
                res.a[i][j]=(res.a[i][j]+(x.a[i][k]*y.a[k][j])%MOD)%MOD;
    return res;
}
mat PowerMod(mat x,ll n)
{
    mat ans={0};
    for(int i=0;i<4;i++)
        ans.a[i][i]=1;
    while(n>0)
    {
        if(n%2==1)
            ans=mult(ans,x);
        x=mult(x,x);
        n=n/2;
    }
    return ans;
}
long long powermod(long long a,long long b,long long c)
{
    long long ans=1;
    long long k;
    k=a;
    k=k%c;
    while(b>0)
    {
        if(b&1)
            ans=(ans*k)%c;
        b>>=1;
        k=(k*k)%c;
    }
    return ans;
}
int main()
{
    int t;
    ll n,y,x,s;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lld%lld%lld%lld",&n,&y,&x,&s);
        MOD=euler(s+1);
        m.a[0][0]=4;m.a[0][1]=1;m.a[0][2]=4;m.a[0][3]=0;
        m.a[1][0]=1;m.a[1][1]=0;m.a[1][2]=0;m.a[1][3]=0;
        m.a[2][0]=2;m.a[2][1]=0;m.a[2][2]=1;m.a[2][3]=0;
        m.a[3][0]=4;m.a[3][1]=1;m.a[3][2]=4;m.a[3][3]=1;
        r=PowerMod(m,n*y-1);
        ll temp=(r.a[3][0]+r.a[3][3])%MOD+MOD;
        printf("%lld\n",powermod(x,temp,s+1));
    }
    return 0;
}