【COGS2187】帕秋莉的超级多项式

19 篇文章 0 订阅
7 篇文章 0 订阅

Description

求:


Solution

直接模拟即可,具体操作可参考之前的博客。

多项式幂次还没有写过,其实很简单, F k = e k ln ⁡ F F^k=e^{k\ln F} Fk=eklnF即可,注意两点:

  • 如果 F F F有常数项不是 1 1 1或者没有常数项要特殊处理
  • 泰勒展开后的 k k k并没有出现在指数上,所以直接对模数取模而不是模数-1。

Code

/************************************************
 * Au: Hany01
 * Date: Jul 30th, 3018
 * Prob: 帕秋莉的超级多项式
 * Email: hany01@foxmail.com
 * Inst: Yali High School
************************************************/

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;
#define rep(i, j) for (register int i = 0, i##_end_ = (j); i < i##_end_; ++ i)
#define For(i, j, k) for (register int i = (j), i##_end_ = (k); i <= i##_end_; ++ i)
#define Fordown(i, j, k) for (register int i = (j), i##_end_ = (k); i >= i##_end_; -- i)
#define Set(a, b) memset(a, b, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define x first
#define y second
#define pb(a) push_back(a)
#define mp(a, b) make_pair(a, b)
#define SZ(a) ((int)(a).size())
#define INF (0x3f3f3f3f)
#define INF1 (2139062143)
#define debug(...) fprintf(stderr, __VA_ARGS__)
#define y1 wozenmezhemecaia

template <typename T> inline bool chkmax(T &a, T b) { return a < b ? a = b, 1 : 0; }
template <typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }

inline int read() {
	static int _, __; static char c_;
    for (_ = 0, __ = 1, c_ = getchar(); c_ < '0' || c_ > '9'; c_ = getchar()) if (c_ == '-') __ = -1;
    for ( ; c_ >= '0' && c_ <= '9'; c_ = getchar()) _ = (_ << 1) + (_ << 3) + (c_ ^ 48);
    return _ * __;
}

const int MOD = 998244353, g0 = 3, maxn = 1e5 + 5;

int rev[maxn << 2], powg[maxn << 2], ipowg[maxn << 2], ig0, invn;

inline int ad(int x, int y) { if ((x += y) >= MOD) return x - MOD; return x; }

inline LL Pow(LL a, LL b) {
	static LL Ans;
	for (Ans = 1; b; b >>= 1, (a *= a) %= MOD) if (b & 1) (Ans *= a) %= MOD;
	return Ans;
}

inline void Init(int N) {
	ig0 = Pow(g0, MOD - 2);
	for (register int i = 1; i <= N; i <<= 1)
		powg[i] = Pow(g0, (MOD - 1) / i), ipowg[i] = Pow(ig0, (MOD - 1) / i);
}

inline void pre(int n) {
	static int N, cnt;
	for (N = 1, cnt = 0; N < n; N <<= 1, ++ cnt);
	rep(i, n) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (cnt - 1));
	invn = Pow(n, MOD - 2);
}

inline void NTT(int* a, int n, int ty) {
	rep(i, n) if (i < rev[i]) swap(a[i], a[rev[i]]);
	for (register int i = 2, p = 1; i <= n; p = i, i <<= 1) {
		register int w0 = ty > 0 ? powg[i] : ipowg[i];
		for (register int j = 0; j < n; j += i) {
			register int w = 1;
			rep(k, p) {
				register int x = a[j + k], y = (LL)a[j + k + p] * w % MOD;
				a[j + k] = ad(x, y), a[j + k + p] = ad(x, MOD - y);
				w = (LL)w * w0 % MOD;
			}
		}
	}
	if (ty < 1) rep(i, n) a[i] = (LL)a[i] * invn % MOD;
}

inline void Mult(int* a, int* b, int* c, int n) {
	static int A[maxn << 2], B[maxn << 2];
	rep(i, n) A[i] = a[i], B[i] = b[i];
	pre(n), NTT(A, n, 1), NTT(B, n, 1);
	rep(i, n) A[i] = (LL)A[i] * B[i] % MOD;
	NTT(A, n, -1);
	rep(i, n) c[i] = A[i];
}

namespace Inv {
	static int A[maxn << 2], B[maxn << 2], a[maxn << 2];
	void Inv_(int* b, int n) {
		if (n == 1) { b[0] = Pow(a[0], MOD - 2); return; }
		Inv_(b, n >> 1);
		rep(i, n) A[i] = a[i], B[i] = b[i];
		For(i, n, n << 1) A[i] = B[i] = 0;
		pre(n << 1), NTT(A, n << 1, 1), NTT(B, n << 1, 1);
		rep(i, n << 1) A[i] = (LL)A[i] * B[i] % MOD * B[i] % MOD;
		NTT(A, n << 1, -1);
		rep(i, n) b[i] = ad(ad(b[i], b[i]), MOD - A[i]);
	}
	inline void Inv(int* x, int* y, int n) {
		rep(i, n) a[i] = x[i], y[i] = 0;
		Inv_(y, n);
	}
}

inline void Int(int* a, int* b, int n) {
	Fordown(i, n - 1, 0) b[i + 1] = Pow(i + 1, MOD - 2) * a[i] % MOD;//downto!!!
	b[0] = 0;
}

inline void Der(int *a, int* b, int n) {
	rep(i, n - 1) b[i] = (LL)a[i + 1] * (i + 1) % MOD;
}

inline void Ln(int* a, int* b, int n) {
	static int A[maxn << 2], B[maxn << 2];
	rep(i, n) A[i] = B[i] = a[i];
	For(i, n, n << 1) A[i] = B[i] = 0;
	Inv:: Inv(A, A, n), Der(B, B, n << 1);
	Mult(A, B, A, n << 1), Int(A, A, n << 1);
	rep(i, n) b[i] = A[i];
}

namespace Exp {
	static int a[maxn << 2], A[maxn << 2];
	void Exp_(int* b, int n) {
		if (n == 1) { b[0] = 1; return; }
		Exp_(b, n >> 1), Ln(b, A, n);
		rep(i, n) A[i] = ad(a[i], MOD - A[i]);
		For(i, n, n << 1) A[i] = 0, b[i] = 0;
		++ A[0], Mult(b, A, b, n << 1);
	}
	void Exp(int* x, int* y, int n) {
		rep(i, n) a[i] = x[i], y[i] = 0;
		Exp_(y, n);
	}
}

namespace Sqrt
{
    static int A[maxn << 2], B[maxn << 2], a[maxn << 2];
    void Sqrt_(int* b, int n) {
        if (n == 1) { b[0] = sqrt(a[0]); return; }
        Sqrt_(b, n >> 1);
        rep(i, n) A[i] = b[i];
        For(i, n, n << 1) A[i] = 0;
        pre(n << 1), NTT(A, n << 1, 1);
        rep(i, n << 1) A[i] = (LL)A[i] * A[i] % MOD;
        NTT(A, n << 1, -1);
        rep(i, n) A[i] = ad(A[i], a[i]), B[i] = ad(b[i], b[i]);
        Inv:: Inv(B, B, n);
        pre(n << 1), NTT(A, n << 1, 1), NTT(B, n << 1, 1);
        rep(i, n << 1) A[i] = (LL)A[i] * B[i] % MOD;
        NTT(A, n << 1, -1); rep(i, n) b[i] = A[i];
    }
    void Sqrt(int* x, int* y, int n) {
        rep(i, n) a[i] = x[i], y[i] = 0;
        Sqrt_(y, n);
    }
}

inline void PPow(int* a, int k, int* b, int n) {
	Ln(a, b, n);
	rep(i, n) b[i] = (LL)b[i] * k % MOD;
	Exp:: Exp(b, b, n);
}

int main()
{
	freopen("polynomial.in", "r", stdin);
	freopen("polynomial.out", "w", stdout);

	static int n, k, f[maxn << 2], N;

	for (n = read(), k = read(), N = 1; N < n; N <<= 1);
	Init(N << 1);
	rep(i, n) f[i] = read();
	Sqrt:: Sqrt(f, f, N);
	Inv:: Inv(f, f, N);
	Int(f, f, N);
	Exp:: Exp(f, f, N);
	Inv:: Inv(f, f, N), ++ f[0];
	Ln(f, f, N), ++ f[0];
	PPow(f, k, f, N);
	Der(f, f, N);
	rep(i, n - 1) printf("%d ", f[i]);
	puts("0");

    return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值