退役在即,前来温习,遇到毒瘤题。
想出了做法,推出了式子,死在了程序实现上。
这种高精度矩阵快速高精度幂是在逗我吗?
已经尽量工程化的代码:
// TODO 下辈子再做吧
/* Header {{{ */
#include <bits/stdc++.h>
using namespace std;
typedef long long readtype;
typedef long long var;
typedef long double let;
readtype read() {
readtype a = 0, c = getchar(), s = 0;
while (!isdigit(c)) s |= c == '-', c = getchar();
while (isdigit(c)) a = a * 10 + c - 48, c = getchar();
return s ? -a : a;
}
#ifdef LOCAL_LOGGER
#define logger(...) fprintf(stderr, __VA_ARGS__)
#else
#define logger(...);
#endif
/* }}} */
class Poly;
class Bignum;
class Matrix;
class Poly {
private:
static const int N = 1e5 + 1;
static const int MOD = 998244353;
static const int G = 3;
// {{{
inline int add(int a, int b) { return (a + b >= MOD) ? (a + b - MOD) : (a + b); }
inline int mul(var a, int b) { return (a * b) % MOD; }
inline int qpow(int x, var y) {
int res = 1;
while (y) {
if (y & 1) res = mul(res, x);
x = mul(x, x);
y >>= 1;
}
return res;
}
inline int inv(int t) { return qpow(t, MOD - 2); }
int i, j;
int limit, inv_limit, rev[N], g[N];
int tmp_a[N], tmp_b[N];
void InitLimit(const int need_len);
void NTT(int *val, const int type);
// }}}
public:
void Mul(Bignum &a, Bignum &b, Bignum &c);
} poly;
class Bignum {
private:
static const int N = 1e5 + 1;
// {{{
int v[N], len;
int &operator[](int k) { return v[k]; }
void format() {
for (int i = 0; i <= len; ++i)
v[i + 1] += v[i] / 10, v[i] %= 10;
while (v[len + 1]) {
len++;
v[len + 1] += v[len] / 10;
v[len] %= 10;
}
}
// }}}
public:
Bignum() {
memset(v, 0, sizeof(v));
len = 0;
}
Bignum(int x) {
v[len = 0] = x;
format();
}
void read() {
len = 0;
int c = getchar();
while (!isdigit(c)) c = getchar();
while (isdigit(c))
v[len++] = c - '0', c = getchar();
reverse(v, v + len);
len--;
}
void print(char c = 0) {
for (int i = len; i >= 0; --i) putchar(v[i] + '0');
if (c) putchar(c);
}
friend Bignum operator + (Bignum &a, Bignum &b);
friend Bignum operator * (Bignum &a, Bignum &b);
friend Bignum operator / (Bignum &a, int b);
friend class Poly;
friend class Matrix;
};
class Matrix {
private:
int v[2][2];
public:
int *operator [] (int k) { return v[k]; }
};
Bignum n, m, res;
int t;
int main() {
#ifndef ONLINE_JUDGE
freopen("3933.in", "r", stdin);
freopen("3933.out", "w", stdout);
#endif
#ifdef LOCAL_LOGGER
freopen("3933.log", "w", stderr);
#endif
n.read(), t = read(), m.read();
return 0;
}
/* ==== Makefile ==== {{{
CompileAndRun:
make Compile
make Run
Compile:
g++ -o 3933 3933.cpp -g -Wall -DLOCAL_LOGGER
CompileUF:
g++ -o 3933 3933.cpp -g -Wall -DLOCAL_LOGGER -fsanitize=undefined
Run:
./3933 < 3933.in > 3933.out
==================
}}} */
void Poly::InitLimit(const int need_len) {
if (limit > need_len) return ;
for (limit = 1; limit <= need_len; limit <<= 1);
inv_limit = inv(limit);
for (i = 0; i < limit; ++i)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (limit >> 1) : 0);
g[0] = 1, g[1] = qpow(G, (MOD - 1) / limit);
for (int i = 2; i < limit; ++i)
g[i] = mul(g[i - 1], g[1]);
}
void Poly::NTT(int *val, const int type) {
for (int i = 0; i < limit; ++i)
if (i < rev[i]) swap(val[i], val[rev[i]]);
for (int mid = 1; mid < limit; mid <<= 1) {
int R = mid << 1, delta = limit / R;
for (int i = 0; i < limit; i += R) {
for (int j = 0; j < mid; ++j) {
int v = mul(g[j * delta], val[i + mid + j]);
val[i + mid + j] = add(val[i + j], MOD - v);
val[i + j] = add(val[i + j], v);
}
}
}
if (type == 1) return ;
reverse(val + 1, val + limit);
for (int i = 0; i < limit; ++i)
val[i] = mul(val[i], inv_limit);
}
void Poly::Mul(Bignum &a, Bignum &b, Bignum &c) {
c.len = a.len + b.len;
InitLimit(c.len);
for (i = 0; i <= a.len; ++i) tmp_a[i] = a[i];
for (i = 0; i <= b.len; ++i) tmp_b[i] = b[i];
NTT(tmp_a, 1), NTT(tmp_b, 1);
for (int i = 0; i < limit; ++i) {
c[i] = mul(tmp_a[i], tmp_b[i]);
tmp_a[i] = 0, tmp_b[i] = 0;
}
NTT(c.v, -1);
}
Bignum operator + (Bignum &a, Bignum &b) {
Bignum c;
c.len = max(a.len, b.len);
for (int i = 0; i <= c.len; ++i) c[i] = a[i] + b[i];
c.format();
return c;
}
Bignum operator * (Bignum &a, Bignum &b) {
Bignum c;
poly.Mul(a, b, c);
c.format();
return c;
}
Bignum operator / (Bignum &a, int b) {
Bignum c;
c.len = a.len;
int v = 0;
for (int i = a.len; i >= 0; --i) {
(v *= 10) += a[i];
c[i] = v / b, v %= b;
}
while (c.len && !c[c.len]) c.len--;
c.format();
return c;
}