LibreOJ 3058 「HNOI2019」白兔之舞

为了好看，设 W=w（

用点普及知识可以知道答案实际上就是 anst=L∑m=1[m≡t(modk)](Lm)Wmx,y

然后套个单位根反演 anst=L∑m=1[m≡t(modk)](Lm)Wmx,y=1kk−1∑d=0ω−dtkL∑m=1(Lm)Wmx,yωdmk=1kk−1∑d=0ω−dtk(Wωdk+I)Lx,y

设 cd=(Wωdk+I)Lx,y。
那么原式就变成了 anst=1kk−1∑d=0cdω−dtk

然后你发现这是 DFT 的标准形式。
然后你敲了个 FFT 上去。
然后然后然后你发现 k 不一定是二的次幂。

然后考虑多点求值（大雾）
然后考虑使用 Bluestein's Algorithm。
（没有写进标签是因为个人感觉这个东西在 OI 界大概也就当个 trick 用）
大概就是把 dt 拆成 (d+t)22−d22−t22。

然后会发现好像没有二次剩余耶（
那么就拆成 (d+t2)−(d2)−(t2) 就好辣（

anst=1kk−1∑d=0cdω−dtk=1kk−1∑d=0cdω−(d+t2)+(d2)+(t2)k=ω(t2)kk−1∑d=0cdω(d2)k⋅ω−(d+t2)k

设 fi=ciω(i2)k,gi=ω−(i2)k，那么就是求一个 anst=ω(t2)kk−1∑d=0fdgd+t

这个可以简单地设 g′2k−i−1=gi 解决。 anst=ω(t2)kk−1∑d=0fdg′2k−d−t−1

当然，要先求原根再动手（

代码：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define add(a,b) (a + b >= mod ? a + b - mod : a + b)
#define dec(a,b) (a < b ? a - b + mod : a - b)
using namespace std;
const int N = 3;
const int K = 1 << 16;
int n,k,ki,L,x,y,mod,len;
int G,rt[K + 5],ans[K + 5];
inline int fpow(int a,int b)
{
    int ret = 1;
    for(;b;b >>= 1)
        (b & 1) && (ret = (long long)ret * a % mod),a = (long long)a * a % mod;
    return ret;
}
struct Matrix
{
    int a[N + 5][N + 5];
    inline Matrix()
    {
        memset(a,0,sizeof a);
    }
    inline Matrix(int)
    {
        memset(a,0,sizeof a);
        for(register int i = 1;i <= n;++i)
            a[i][i] = 1;
    }
    inline int *operator[](const int &x)
    {
        return a[x];
    }
    inline const int *operator[](const int &x) const
    {
        return a[x];
    }
    inline Matrix operator*(const int &x) const
    {
        Matrix ret;
        for(register int i = 1;i <= n;++i)
            for(register int j = 1;j <= n;++j)
                ret[i][j] = (long long)a[i][j] * x % mod;
        return ret;
    }
    inline Matrix operator+(const Matrix &o) const
    {
        Matrix ret;
        for(register int i = 1;i <= n;++i)
            for(register int j = 1;j <= n;++j)
                ret[i][j] = (a[i][j] + o[i][j]) % mod;
        return ret;
    }
    inline Matrix operator*(const Matrix &o) const
    {
        Matrix ret;
        for(register int i = 1;i <= n;++i)
            for(register int j = 1;j <= n;++j)
                for(register int k = 1;k <= n;++k)
                    ret[i][j] = (ret[i][j] + (long long)a[i][k] * o[k][j]) % mod;
        return ret;
    }
} w;
inline Matrix fpow(Matrix a,int b)
{
    Matrix ret(1);
    for(;b;b >>= 1)
        (b & 1) && (ret = ret * a,1),a = a * a;
    return ret;
}
namespace Poly
{
    const int N = 1 << 18;
    const double pi = acos(-1);
    const int W = 1 << 15;
    int n;
    int lg2[N + 5],rev[N + 5];
    struct poly
    {
        int a[N + 5];
        inline const int &operator[](int x) const
        {
            return a[x];
        }
        inline int &operator[](int x)
        {
            return a[x];
        }
        inline void clear(int x = 0)
        {
            memset(a + x,0,(N - x + 1) << 2);
        }
    } f,g;
    struct cp
    {
        double a,b;
        inline void operator+=(const cp &o)
        {
            a += o.a,b += o.b;
        }
        inline cp operator+(const cp &o) const
        {
            return (cp){a + o.a,b + o.b};
        }
        inline cp operator-(const cp &o) const
        {
            return (cp){a - o.a,b - o.b};
        }
        inline cp operator*(const cp &o) const
        {
            return (cp){a * o.a - b * o.b,a * o.b + b * o.a};
        }
        inline void operator*=(const cp &o)
        {
            *this = *this * o;
        }
        inline cp operator*(const double &o) const
        {
            return (cp){a * o,b * o};
        }
        inline cp operator~() const
        {
            return (cp){a,-b};
        }
    } rt[N + 5];
    inline void init(int len)
    {
        for(n = 1;n < len;n <<= 1);
        for(register int i = 2;i <= n;++i)
            lg2[i] = lg2[i >> 1] + 1;
        for(register int i = 0;i <= (n >> 1);++i)
            rt[(n >> 1) + i] = (cp){cos(2 * pi * i / n),sin(2 * pi * i / n)};
        for(register int i = (n >> 1) - 1;i;--i)
            rt[i] = rt[i << 1];
    }
    inline void fft(cp *a,int type,int n)
    {
        type == -1 && (reverse(a + 1,a + n),1);
        int lg = lg2[n] - 1;
        for(register int i = 0;i < n;++i)
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg),
            i < rev[i] && (swap(a[i],a[rev[i]]),1);
        for(register int w = 2,m = 1;w <= n;w <<= 1,m <<= 1)
            for(register int i = 0;i < n;i += w)
                for(register int j = 0;j < m;++j)
                {
                    cp t = rt[m | j] * a[i | j | m];
                    a[i | j | m] = a[i | j] - t,a[i | j] += t;
                }
        if(type == -1)
            for(register int i = 0;i < n;++i)
                a[i].a /= n,a[i].b /= n;
    }
    inline void mul(poly &a,const poly &b,int n)
    {
        static cp f[N + 5],g[N + 5],h[N + 5];
        int lim = 1;
        memset(f,0,sizeof f),memset(g,0,sizeof g);
        for(;lim < (n << 1);lim <<= 1);
        for(register int i = 0;i < n;++i)
            f[i] = (cp){a[i] / W,a[i] % W},g[i] = (cp){b[i] / W,b[i] % W};
        fft(f,1,lim),fft(g,1,lim);
        for(register int i = 0;i < lim;++i)
            h[i] = ~f[(lim - i) % lim];
        for(register int i = 0;i < lim;++i)
            f[i] *= g[i],g[i] *= h[i];
        fft(f,-1,lim),fft(g,-1,lim);
        for(register int i = 0;i < lim;++i)
        {
            long long ac = (f[i].a + g[i].a) / 2 + 0.5;
            long long bd = g[i].a - ac + 0.5;
            long long bcad = f[i].b + 0.5;
            a[i] = ((ac % mod * W % mod * W % mod) % mod + (bcad % mod * W % mod) % mod + bd % mod) % mod;
        }
    }
}
int p[N + 5],t,cnt;
int main()
{
    scanf("%d%d%d%d%d%d",&n,&k,&L,&x,&y,&mod);
    len = k << 1,Poly::init(len << 1);
    t = mod - 1;
    for(register int i = 2;i * i <= t;++i)
    {
        if(!(t % i))
            p[++cnt] = i;
        for(;!(t % i);t /= i);
    }
    t > 1 && (p[++cnt] = t);
    for(register int i = 2,flag;i < mod;++i)
    {
        flag = 1;
        for(register int j = 1;j <= cnt;++j)
            if(fpow(i,(mod - 1) / p[j]) == 1)
            {
                flag = 0;
                break;
            }
        if(flag)
        {
            G = i;
            break;
        }
    }
    for(register int i = 0;i < k;++i)
        rt[i] = fpow(G,(mod - 1) / k * i);
    for(register int i = 1;i <= n;++i)
        for(register int j = 1;j <= n;++j)
            scanf("%d",w.a[i] + j);
    for(register int i = 0;i < k;++i)
        Poly::f[i] = (long long)fpow(w * rt[i] + Matrix(1),L)[x][y] * rt[(long long)i * (i - 1) / 2 % k] % mod;
    for(register int i = 0;i < len - 1;++i)
        Poly::g[len - i - 1] = rt[(k - (long long)i * (i - 1) / 2 % k) % k];
    Poly::mul(Poly::f,Poly::g,len);
    ki = fpow(k,mod - 2);
    for(register int i = 0;i < k;++i)
        ans[i] = (long long)Poly::f[len - i - 1] * rt[(long long)i * (i - 1) / 2 % k] % mod * ki % mod;
    for(register int i = 0;i < k;++i)
        printf("%d\n",ans[i]);
}