UOJ#272. 【清华集训2016】石家庄的工人阶级队伍比较坚强

传送门

设运算 \(op1,op2\)，一个表示三进制不进位的加法，一个表示不退位的减法
设 \(cnt1[x],cnt2[x]\) 分别表示 \(x\) 转成三进制后 \(1/2\) 的个数
那么
\(f_{i,x}=\sum f_{i-1,y}b_{cnt1[x~op2~y],cnt2[x~op2~y]}\)
设 \(B_{x,y}=b_{cnt1[x~op2~y],cnt2[x~op2~y]}\)
那么可以发现 \(B_{x,y}=B_{x~op2~y,0}\)
那么我们要求的就是 \(f\) 与 \(B\) 的第一行的 \(t\) 次卷积的卷积
其中下标运算为 \(op1\)
那么我们求出 \(f\) 和 \(B\) 的"点值表达"，快速幂之后变换回去即可
下标运算可以看成是每一位的模 \(3\) 的循环卷积，用三次单位根 \(FWT\)，每一层手动做一遍长度为 \(3\) 的 \(DFT\)
由于题目中 \(p\) 的性质，可以得到 \(3\perp p\)，所以 \(3\) 有逆元
注意到 \(\omega_3^2+\omega_3+1=0\)
把所有数字用 \(a+b\omega\) 表示，重定义运算即可

# include <bits/stdc++.h>
using namespace std;
typedef long long ll;

namespace IO {
    const int maxn(1 << 21 | 1);

    char ibuf[maxn], obuf[maxn], *iS, *iT, *oS = obuf, *oT = obuf + maxn - 1, c, st[66];
    int tp, f;

    inline char Getc() {
        return iS == iT ? (iT = (iS = ibuf) + fread(ibuf, 1, maxn, stdin), (iS == iT ? EOF : *iS++)) : *iS++;
    }

    template <class Int> inline void In(Int &x) {
        for (f = 1, c = Getc(); c < '0' || c > '9'; c = Getc()) f = c == '-' ? -1 : 1;
        for (x = 0; c >= '0' && c <= '9'; c = Getc()) x = (x << 1) + (x << 3) + (c ^ 48);
        x *= f;
    }

    inline void Flush() {
        fwrite(obuf, 1, oS - obuf, stdout);
        oS = obuf;
    }

    inline void Putc(char c) {
        *oS++ = c;
        if (oS == oT) Flush();
    }

    template <class Int> void Out(Int x) {
        if (!x) Putc('0');
        if (x < 0) Putc('-'), x = -x;
        while (x) st[++tp] = x % 10 + '0', x /= 10;
        while (tp) Putc(st[tp--]);
    }
}

using IO :: In;
using IO :: Out;
using IO :: Putc;
using IO :: Flush;

const int maxn(531441);

int mod, m, t, n, bin[20], b[20][20], cnt1[maxn], cnt2[maxn], phi, inv3;

inline int Pow(ll x, int y) {
    ll ret = 1;
    for (; y; y >>= 1, x = x * x % mod)
        if (y & 1) ret = ret * x % mod;
    return ret;
}

inline void Inc(int &x, int y) {
    x = x + y >= mod ? x + y - mod : x + y;
}

inline void Dec(int &x, int y) {
    x = x - y < 0 ? x - y + mod : x - y;
}

inline int Add(int x, int y) {
    return x + y >= mod ? x + y - mod : x + y;
}

inline int Sub(int x, int y) {
    return x - y < 0 ? x - y + mod : x - y;
}

struct Complex {
    int a, b;

    inline Complex(int _a = 0, int _b = 0) {
        a = _a, b = _b;
    }

	inline Complex W1() {
		return Complex(Sub(0, b), Sub(a, b));
	}

	inline Complex W2() {
		return Complex(Sub(b, a), Sub(0, a));
	}

    inline Complex operator +(Complex y) const {
        return Complex(Add(a, y.a), Add(b, y.b));
    }

    inline Complex operator -(Complex y) const {
        return Complex(Sub(a, y.a), Sub(b, y.b));
    }

    inline Complex operator *(Complex y) const {
        return Complex(Sub((ll)a * y.a % mod, (ll)b * y.b % mod), Sub(Add((ll)a * y.b % mod, (ll)b * y.a % mod), (ll)b * y.b % mod));
    }

    inline Complex operator *(int y) const {
        return Complex((ll)a * y % mod, (ll)b * y % mod);
    }
} coef[maxn], f[maxn], tmp[3];

inline Complex PowComplex(Complex x, int y) {
	Complex ret = Complex(1, 0);
	for (; y; y >>= 1, x = x * x) if (y & 1) ret = ret * x;
	return ret;
}

inline void DFWT(Complex *p, int opt) {
	int i, j, k, t;
	for (i = 1; i < n; i *= 3)
		for (j = 0, t = i * 3; j < n; j += t)
			for (k = 0; k < i; ++k) {
				tmp[0] = p[j + k], tmp[1] = p[j + k + i], tmp[2] = p[j + k + i + i];
				p[j + k] = tmp[0] + tmp[1] + tmp[2];
				p[j + k + i] = tmp[0] + tmp[1].W1() + tmp[2].W2();
				p[j + k + i + i] = tmp[0] + tmp[1].W2() + tmp[2].W1();
				if (opt == -1) {
					swap(p[j + k + i], p[j + k + i + i]);
					p[j + k] = p[j + k] * inv3;
					p[j + k + i] = p[j + k + i] * inv3;
					p[j + k + i + i] = p[j + k + i + i] * inv3;
				}
			}
}

int main() {
	freopen("b.in", "r", stdin);
	int i, j, x;
	In(m), In(t), In(mod);
	for (i = bin[0] = 1; i < 20; ++i) bin[i] = bin[i - 1] * 3;
	n = bin[m];
	if (mod == 1) {
		for (i = 0; i < n; ++i) Putc('0'), Putc('\n');
		return Flush(), 0;
	}
	x = phi = mod;
	for (i = 2; i * i <= x; ++i)
		if (x % i == 0) {
			phi -= phi / i;
			while (x % i == 0) x /= i;
		}
	if (x > 1) phi -= phi / x;
	inv3 = Pow(3, phi - 1);
	for (i = 0; i < n; ++i) {
		cnt1[i] = cnt1[i / 3] + (i % 3 == 1);
		cnt2[i] = cnt2[i / 3] + (i % 3 == 2);
	}
	for (i = 0; i < n; ++i) In(f[i].a);
	for (i = 0; i <= m; ++i)
		for (j = 0; j <= m - i; ++j) In(b[i][j]);
	for (i = 0; i < n; ++i) coef[i].a = b[cnt1[i]][cnt2[i]];
	DFWT(coef, 1), DFWT(f, 1);
	for (i = 0; i < n; ++i) f[i] = f[i] * PowComplex(coef[i], t);
	DFWT(f, -1);
	for (i = 0; i < n; ++i) Out(f[i].a), Putc('\n');
	return Flush(), 0;
}