我们
TM怎么又要上文化课。。我 哔哔哔哔哔哔
有 \(T\) 组数据，求
\[ \prod_{i = 1}^{n} \prod_{j = 1}^{m} fib[\gcd(i, j)] \]
$ 1 \leq n, m \leq 10 ^ 6, 1 \leq T \leq 1000 $
又是一道莫比乌斯反演套路题。。
\[ \begin{aligned} &\prod_{i = 1}^{n} \prod_{j = 1}^{m} fib[\gcd(i, j)]\&=\prod_{d = 1}^{n} fib[d]^{\sum_{i = 1}^n\sum_{j = 1}^n [(i, j) = d]} \end{aligned} \]
写在幂次那里不好看，我们单独提出来。。。
\[ \begin{aligned} & \sum_{i = 1}^n\sum_{j = 1}^n [(i, j) = d]\&= \sum_{d | x} \mu(\frac x d) \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{x} \rfloor\&= \sum_{i = 1} ^ {\lfloor \frac{n}{d} \rfloor} \mu(i) \lfloor \frac{n}{id} \rfloor \lfloor \frac{m}{id} \rfloor \end{aligned} \]
我们令 \(T = id\) 那么原式就变成
\[ \begin{aligned} &\prod_{d = 1}^{n} fib[d]^{\sum_{i = 1} ^ {\lfloor \frac{n}{d} \rfloor} \mu(i) \lfloor \frac{n}{id} \rfloor \lfloor \frac{m}{id} \rfloor}\&= \prod_{T = 1}^n (\prod_{d | T} fib[d]^{\mu(\frac{T}{d} )})^{\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor} \end{aligned} \]
预处理中间那个括号的内容，然后外面套整除分块即可。
其实不预处理，套两层整除分块 \(\mathcal LOJ\) 的机子可以卡过 QAQ
复杂度是 \(\mathcal O(n \log n + T \sqrt n \log n)\) 的。
#include <bits/stdc++.h>
#define For(i, l, r) for(register int i = (l), _end_ = (int)(r); i <= _end_; ++i)
#define Fordown(i, r, l) for(register int i = (r), _end_ = (int)(l); i >= _end_; --i)
#define Set(a, v) memset(a, v, sizeof(a))
using namespace std;
bool chkmin(int &a, int b) {return b < a ? a = b, 1 : 0;}
bool chkmax(int &a, int b) {return b > a ? a = b, 1 : 0;}
inline int read() {
int x = 0, fh = 1; char ch = getchar();
for (; !isdigit(ch); ch = getchar() ) if (ch == '-') fh = -1;
for (; isdigit(ch); ch = getchar() ) x = (x<<1) + (x<<3) + (ch ^ '0');
return x * fh;
}
typedef long long ll;
int n, m;
const int N = 1e6 + 1e3;
ll mu[N], prime[N / 4], fib[N], ifib[N], cnt, g[N], ig[N];
bitset<N> is_prime;
const int Mod = 1e9 + 7, Pmod = Mod - 1;
ll fpm(ll x, ll power) {
ll res = 1;
for (; power; power >>= 1, (x *= x) %= Mod)
if (power & 1) (res *= x) %= Mod;
return res;
}
void Init(int maxn) {
is_prime.set();
is_prime[0] = is_prime[1] = false;
int res; g[0] = ig[0] = g[1] = ifib[1] = mu[1] = fib[1] = 1;
For (i, 2, maxn) {
if (is_prime[i]) {
prime[++cnt] = i;
mu[i] = -1;
}
For (j, 1, cnt) {
res = prime[j] * i;
if (res > maxn) break ;
is_prime[res] = false;
if (i % prime[j]) mu[res] = -mu[i];
else {mu[res] = 0; break ;}
}
g[i] = 1;
fib[i] = fib[i - 1] + fib[i - 2];
if (fib[i] >= Mod) fib[i] -= Mod;
ifib[i] = fpm(fib[i], Mod - 2);
}
For (i, 1, maxn) {
for (register int j = i, t = 1; j <= maxn; j += i, ++ t)
if (mu[t] == 1) g[j] *= fib[i], g[j] = g[j] >= Mod ? g[j] % Mod : g[j];
else if (mu[t] == -1) g[j] *= ifib[i], g[j] = g[j] >= Mod ? g[j] % Mod : g[j];
(g[i] *= g[i - 1]) %= Mod;
ig[i] = fpm(g[i], Mod - 2);
}
}
ll Calc() {
cerr << "Calc: " << endl;
ll res = 1; ll Nextd, n_, m_;
if (n > m) swap(n, m);
For(d, 1, n) {
n_ = n / d; m_ = m / d;
Nextd = min(n / n_, m / m_);
(res *= fpm(g[Nextd] * ig[d - 1] % Mod, n_ * m_ % Pmod) ) %= Mod;
d = Nextd;
}
return res;
}
void File() {
#ifdef zjp_shadow
freopen ("2000.in", "r", stdin);
freopen ("2000.out", "w", stdout);
#endif
}
int main() {
File();
Init((int)1e6);
int cases = read();
while (cases --) {
n = read(); m = read();
printf ("%lld\n", Calc() );
}
return 0;
}看我之前的 题解 就好啦。
似乎说是一个套路题？染上没有的颜色就是 access 的过程，然后套上线段树维护就好啦。
求有多少个长度为 $ n $ 的序列，至少有一个数是质数，序列中的数都是不超过 $ m $ 的正整数，而且这 $ n $ 个数的和是 $ p $ 的倍数。答案对于 \(20170408\) 取模。
$ 1 \leq n \leq 10 ^ 9, 1 \leq m \leq 2 \times 10 ^ 7, 1 \leq p \leq 100 $。
看到数据就知道是思博矩幂了，可以做到 \(\mathcal O(m + p^3 \log n)\) 。
其实这是个循环矩阵，无限自乘还是循环的。我们只需要维护其中一行就知道其它所有行了，那每次拿个向量搞搞得出一行就行啦。
这样可以做到 \(\mathcal O(m + p^2 \log n)\) 。
其实那个过程本质上是循环卷积，但是和 \(FFT\) 关于 \(2^k\) 循环似乎不太契合。
那么这个多项式快速幂每次只能做到 \(\mathcal O(m + p \log p \log n)\) 。
其实可以利用 ln 变成加法，然后 exp 回去即可，做到 \(\mathcal O(m + p \log p)\) 。
但由于模数的鬼畜，不能用 NTT 利用 MTT 实现就好啦。。。
显然我的做法是最暴力的矩幂啦 qwq
#include <bits/stdc++.h>
#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl
using namespace std;
template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }
inline int read() {
int x(0), sgn(1); char ch(getchar());
for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
return x * sgn;
}
void File() {
#ifdef zjp_shadow
freopen ("2002.in", "r", stdin);
freopen ("2002.out", "w", stdout);
#endif
}
const int N = 2e7 + 1e3, M = 110, Mod = 20170408;
int n, m, p;
bitset<N> is_prime; int prime[N], plen;
void Linear_Sieve(int maxn) {
is_prime.set();
is_prime[0] = is_prime[1] = false;
For (i, 2, maxn) {
if (is_prime[i])
prime[++ plen] = i;
For (j, 1, plen) {
int res = prime[j] * i;
if (res > maxn) break;
is_prime[res] = false;
if (!(i % prime[j])) break;
}
}
}
int ptot[M], tot[M];
struct Mat {
int a[M][M];
Mat() { Set(a, 0); }
inline void Unit() {
Rep (i, p) a[i][i] = 1;
}
inline Mat friend operator * (const Mat &lhs, const Mat &rhs) {
Mat res;
Rep (i, p) Rep (k, p) Rep (j, p)
res.a[i][j] = (res.a[i][j] + 1ll * lhs.a[i][k] * rhs.a[k][j]) % Mod;
return res;
}
};
inline Mat fpm(Mat x, int power) {
Mat res; res.Unit();
for (; power; power >>= 1, x = x * x)
if (power & 1) res = res * x;
return res;
}
Mat f, trans;
int main () {
File();
n = read(); m = read(); p = read();
Linear_Sieve(m);
For (i, 1, m) {
++ tot[i % p];
if (is_prime[i]) ++ ptot[i % p];
}
Rep (i, p) Rep (j, p)
trans.a[(i + j) % p][i] += tot[j];
f.a[0][0] = 1; f = fpm(trans, n) * f;
int ans = f.a[0][0];
Rep (i, p) Rep (j, p)
trans.a[(i + j) % p][i] -= ptot[j];
Set(f.a, 0); f.a[0][0] = 1; f = fpm(trans, n) * f;
ans = (ans + Mod - f.a[0][0]) % Mod;
printf ("%d\n", ans);
return 0;
}左右各有 \(n\) 个点的二分图， \(a_{i, j}, b_{i, j}\) 为 \(i, j\) 匹配能得到的两个权值。
一个方案中有 $ n $ 对匹配，假设每组匹配的权值 \(a_{i, j}\) 分别是 $ a‘_1, a‘_2, \ldots, a‘_n $， \(b_{i, j}\) 分别是 $ b‘_1, b‘_2, \ldots, b‘_n $。令
求一个完美匹配，使得
\[ C = \frac{a'_1 + a'_2 + \cdots + a'_n}{b'_1 + b'_2 + \cdots + b'_n} \]
尽量大。
这 round1 模板题是真多。。。
直接分数规划，二分 \(C\) ，边权就变成 \(a_{i, j} - C \times b_{i, j}\) 。
那么就是找一个完美匹配使得匹配的边权和不小于 \(0\) 。
然后用 KM 或者费用流跑最大权匹配就行了。
#include <bits/stdc++.h>
#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl
using namespace std;
template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }
inline int read() {
int x(0), sgn(1); char ch(getchar());
for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
return x * sgn;
}
void File() {
#ifdef zjp_shadow
freopen ("2003.in", "r", stdin);
freopen ("2003.out", "w", stdout);
#endif
}
const int inf = 0x3f3f3f3f;
template<int Maxn, int Maxm>
struct MCMF {
int Head[Maxn], Next[Maxm], to[Maxm], cap[Maxm], e; double val[Maxm];
void Init() {
e = 1; Set(Head, 0);
}
inline void add_edge(int u, int v, double cost, int flow) {
to[++ e] = v; Next[e] = Head[u]; Head[u] = e;
val[e] = cost; cap[e] = flow;
}
inline void Add(int u, int v, double cost, int flow) {
add_edge(u, v, cost, flow); add_edge(v, u, -cost, 0);
}
int S, T;
int pre[Maxn]; double dis[Maxn]; bitset<Maxn> inq;
bool Spfa() {
queue<int> Q;
For (i, 1, T) dis[i] = - 1e8;
Set(pre, 0); inq.reset();
Q.push(S); dis[S] = 1;
while (!Q.empty()) {
int u = Q.front(); Q.pop(); inq[u] = false;
for (int i = Head[u], v = to[i]; i; v = to[i = Next[i]]) {
if (cap[i] && chkmax(dis[v], dis[u] + val[i])) {
pre[v] = i; if (!inq[v]) inq[v] = true, Q.push(v);
}
}
}
return pre[T];
}
double Run() {
int sum_flow = 0; double sum_cost = 0;
while (Spfa()) {
int flow = inf;
for (int u = T; u != S; u = to[pre[u] ^ 1]) chkmin(flow, cap[pre[u]]);
for (int u = T; u != S; u = to[pre[u] ^ 1]) {
cap[pre[u]] -= flow;
cap[pre[u] ^ 1] += flow;
sum_cost += val[pre[u]] * flow;
}
sum_flow += flow;
}
return sum_cost;
}
};
const int N = 110;
const double eps = 1e-8;
MCMF<N * 2, N * N * 2> T;
int n; double A[N][N], B[N][N];
#define In(x) (x)
#define Out(x) ((x) + n)
inline bool check(double k) {
T.Init();
T.T = (T.S = Out(n) + 1) + 1;
For (i, 1, n) {
T.Add(T.S, In(i), 0, 1);
T.Add(Out(i), T.T, 0, 1);
}
For (i, 1, n) For (j, 1, n)
T.Add(In(i), Out(j), A[i][j] - k * B[i][j], 1);
return T.Run() > 0;
}
int main () {
File();
n = read();
For (i, 1, n) For (j, 1, n) A[i][j] = read();
For (i, 1, n) For (j, 1, n) B[i][j] = read();
double l = 0.0, r = 1e8;
while (r - l > eps) {
double mid = (l + r) / 2.0;
if (check(mid)) l = mid; else r = mid;
}
printf ("%.6lf\n", l);
return 0;
}选出 $ n $ 个同学， 每个同学猜一个长度为 $ m $ 的序列，当某一个同学猜的序列在硬币序列中出现时，就不再扔硬币了，并且这个同学胜利。为了保证只有一个同学胜利，同学们猜的 $ n $ 个序列两两不同。
最后求出每个同学胜利的概率。
$ 1 \leq n, m \leq 300 $
六道题里唯一没有那么思博的题。。
暴力 ac自动机 + 高斯消元 是 \(\mathcal O((nm)^3)\) 的，只有 \(40\) 分。
考虑优化，显然我们最后要求的只有 \(n\) 个变量的值，其他变量都是毫无意义的。
那么我们只需要求出这些变量之间的相关方程即可。
我们多设一个 \(f_s\) 为不匹配到其中任意一个同学的概率，那么我们在当前串后面添加一个同学 \(i\) 的拥有的字符串 \(S_i\) 就能对应到它的概率。
这个似乎利用了一个结论，如果不匹配到 \(n\) 个结束节点，停到任意一个点的概率似乎是相同的。
但是这样会把 \(S_i\) 有个前缀和 \(S_j\) 一个后缀是一样的概率多算上去，我们减掉就行了，这个出现的概率是 \(2^{-l}\) ，\(l\) 为这个前缀的长度。
利用 KMP 实现就可以做到 \(\mathcal O(n^2(n + m))\) 啦。
我懒，用的暴力匹配。。也过啦。。
#include <bits/stdc++.h>
#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl
using namespace std;
template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }
inline int read() {
int x(0), sgn(1); char ch(getchar());
for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
return x * sgn;
}
void File() {
#ifndef ONLINE_JUDGE
freopen ("2004.in", "r", stdin);
freopen ("2004.out", "w", stdout);
#endif
}
const int N = 310;
double G[N][N];
void Gauss(int n) {
For (i, 1, n) {
For (j, i + 1, n) if (fabs(G[j][i]) > fabs(G[i][i])) swap(G[i], G[j]);
if (!G[i][i]) continue;
Fordown (j, n + 1, i) G[i][j] /= G[i][i];
For (j, i + 1, n) Fordown (k, n + 1, i) G[j][k] -= G[j][i] * G[i][k];
}
Fordown (i, n, 1) {
For (j, 1, i - 1)
G[j][n + 1] -= G[j][i] * G[i][n + 1], G[j][i] = 0;
}
}
int n, m; double inv[N];
char str[N][N];
int main() {
File();
n = read(); m = read();
For (i, 1, n) scanf ("%s", str[i] + 1);
inv[0] = 1.0;
For (i, 1, m) inv[i] = inv[i - 1] / 2.0;
For (i, 1, n) {
G[i][i] = 1.0; G[i][n + 1] = - 1.0;
For (j, 1, n) {
For (p, 2, m) {
int len = 1;
while (p + len - 1 <= m && str[i][len] == str[j][p + len - 1]) ++ len;
if (p + len == m + 2) G[i][j] += inv[p - 1];
}
}
}
For (i, 1, n) G[n + 1][i] = 1.0; G[n + 1][n + 2] = 1.0;
Gauss(n + 1);
For (i, 1, n)
printf ("%.10lf\n", G[i][n + 2]);
return 0;
}一个序列，每个元素是一个点 \((x_i, y_i)\) 。
有 \(m\) 个操作，共有三种：
\(\texttt{1 L R}:\) 求 \([L, R]\) 的回归直线的斜率。
\(\texttt{2 L R S T}:\) 把 \(i \in [L, R]\) \(x_i\) 加上 \(S\) ，\(y_i\) 加上 \(T\) 。
\(\texttt{2 L R S T}:\) 把 \(i \in [L, R]\) \(x_i\) 变成 \(S + i\) ，\(y_i\) 变成 \(T + i\) 。
$ 1 \leq n, m \leq 10 ^ 5 $
昨天才讲的回归直线。。今天就用上啦QAQ
\[ \hat{a} = \frac{\sum\limits_{i = L} ^ R (x_i - \bar{x})(y_i - \bar{y})}{\sum\limits_{i = L} ^ R (x_i - \bar{x}) ^ 2} \]
这个数学书上常常用另外一种化简的方式。。
令 \(n = R - L + 1\) 即区间长度。
\[ \hat{a} = \frac{(\sum\limits_{i = L} ^ R x_i y_i) - n\bar{x}\bar{y}}{(\sum\limits_{i = L} ^ R x_i^2) - n \bar{x}^2} \]
然后维护区间 \(x, y, xy, x^2\) 的和就行了。
#include <bits/stdc++.h>
#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl
using namespace std;
template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }
inline int read() {
int x(0), sgn(1); char ch(getchar());
for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
return x * sgn;
}
void File() {
#ifdef zjp_shadow
freopen ("2005.in", "r", stdin);
freopen ("2005.out", "w", stdout);
#endif
}
const int N = 1e5 + 1e3;
using ll = double;
ll pre1(ll x) {
return x * (x + 1) / 2;
}
ll pre2(ll x) {
return x * (x + 1) * (2 * x + 1) / 6;
}
#define ls o << 1
#define rs o << 1 | 1
#define mid ((l + r) >> 1)
#define lson ls, l, mid
#define rson rs, mid + 1, r
template<int Maxn>
struct Segment_Tree {
ll sumxy[Maxn], sumx[Maxn], sumy[Maxn], sumxx[Maxn];
int tag[Maxn], S[Maxn], T[Maxn];
inline void modify(int o, int l, int r, int opt, int sv, int tv) {
int len = r - l + 1;
chkmax(tag[o], opt);
if (opt == 2) {
sumxy[o] += sv * sumy[o] + tv * sumx[o] + 1ll * len * sv * tv;
sumxx[o] += 2 * sv * sumx[o] + 1ll * len * sv * sv;
sumx[o] += 1ll * sv * len;
sumy[o] += 1ll * tv * len;
S[o] += sv; T[o] += tv;
} else {
sv += l; tv += l;
sumxy[o] = 1ll * len * sv * tv + pre1(len - 1) * (sv + tv) + pre2(len - 1);
sumxx[o] = 1ll * len * sv * sv + 2 * pre1(len - 1) * sv + pre2(len - 1);
sumx[o] = 1ll * len * sv + pre1(len - 1);
sumy[o] = 1ll * len * tv + pre1(len - 1);
S[o] = sv - l; T[o] = tv - l;
}
}
inline void push_down(int o, int l, int r) {
if (tag[o]) {
modify(lson, tag[o], S[o], T[o]);
modify(rson, tag[o], S[o], T[o]);
tag[o] = S[o] = T[o] = 0;
}
}
inline void push_up(int o) {
sumxy[o] = sumxy[ls] + sumxy[rs];
sumxx[o] = sumxx[ls] + sumxx[rs];
sumx[o] = sumx[ls] + sumx[rs];
sumy[o] = sumy[ls] + sumy[rs];
}
void Update(int o, int l, int r, int opt, int ul, int ur, int sv, int tv) {
if (ul <= l && r <= ur) {
modify(o, l, r, opt, sv, tv); return;
}
push_down(o, l, r);
if (ul <= mid) Update(lson, opt, ul, ur, sv, tv);
if (ur > mid) Update(rson, opt, ul, ur, sv, tv);
push_up(o);
}
ll Query(int o, int l, int r, int ql, int qr, ll *sum) {
if (ql <= l && r <= qr) return sum[o];
push_down(o, l, r); ll res = 0;
if (ql <= mid) res += Query(lson, ql, qr, sum);
if (qr > mid) res += Query(rson, ql, qr, sum);
push_up(o); return res;
}
void Build(int o, int l, int r, int *x, int *y) {
if (l == r) {
sumxy[o] = 1ll * x[l] * y[r];
sumxx[o] = 1ll * x[l] * x[r];
sumx[o] = x[l]; sumy[o] = y[r];
return;
}
Build(lson, x, y); Build(rson, x, y); push_up(o);
}
};
Segment_Tree<N << 2> T;
int n, m, x[N], y[N];
#define root 1, 1, n
int main () {
File();
n = read(); m = read();
For (i, 1, n) x[i] = read();
For (i, 1, n) y[i] = read();
T.Build(root, x, y);
while (m --) {
int opt = read(), l = read(), r = read();
if (opt == 1) {
ll sumxy = T.Query(root, l, r, T.sumxy),
sumxx = T.Query(root, l, r, T.sumxx),
sumx = T.Query(root, l, r, T.sumx),
sumy = T.Query(root, l, r, T.sumy), len = r - l + 1;
double averx = 1.0 * sumx / len, avery = 1.0 * sumy / len,
k = 1.0 * (sumxy - len * averx * avery) / (sumxx - len * averx * averx);
printf ("%.6lf\n", k);
} else {
int sv = read(), tv = read();
T.Update(root, opt, l, r, sv, tv);
}
}
return 0;
}原文：https://www.cnblogs.com/zjp-shadow/p/10409375.html