Day 1
T1 数字表格
题目大意
· 求(prodlimits_{i=1}^nprodlimits_{j=1}^mFibonacci(gcd(i,j)))，(Tleq1000)，(n,mleq10^6)
思路
· 一言不合化式子（不失一般性地假设(n<m)）$$egin{aligned}
ans&=prod_{i=1}nprod_{j=1}mf((i,j))
&=prod_{d=1}nf(d){sumlimits_{i=1}nsumlimits_{j=1}m[(i,j)=d]}
&=prod_{d=1}nf(d){sumlimits_{k=1}^{frac nd}mu(k)leftlfloorfrac n{dk}
ight
floorleftlfloorfrac m{dk}
ight
floor}
end{aligned}$$
· 到这里直接(Oleft(Tn^{frac 34}
ight))暴力计算答案，常数压的好就可以通过此题。
· 更神奇的：$$prod_{d=1}nf(d){sumlimits_{k=1}^{frac nd}mu(k)leftlfloorfrac n{dk}
ight
floorleftlfloorfrac m{dk}
ight
floor}=prod_{T=1}nleft(prod_{d|T}f(d){muleft(frac Td
ight)}
ight)^{leftlfloorfrac nT
ight
floorleftlfloorfrac mT
ight
floor}$$
· 我们设(g(n)=prodlimits_{d|n}f(d)^{muleft(frac nd
ight)})
· 很类似莫比乌斯反演，(f(n)=prodlimits_{d|n}g(n))，这个结论没什么意义，不过可以从另一个角度来做这道题。
· 如何预处理(g(i))呢？(g(i))并不是积性函数，但我们可以枚举所有斐波那契数，然后再枚举它们的倍数来计算对(g(i))的贡献，时间复杂度(O(nln n))。
· 预处理完成后(Oleft(left(sqrt n+sqrt m
ight)log n
ight))回答单次询问。
· 时间复杂度(Oleft(nln n+Tleft(sqrt n+sqrt m
ight)log n
ight))。
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N = 1000003;
const int p = 1000000007;
bool notp[N];
int mu[N], prime[N], num = 0, g[N], nig[N];
int ipow(int a, int b) {
int r = 1, w = a;
while (b) {
if (b & 1) r = 1ll * r * w % p;
w = 1ll * w * w % p;
b >>= 1;
}
return r;
}
void Euler_shai() {
mu[1] = 1;
for (int i = 2; i <= 1000000; ++i) {
if (!notp[i]) {
prime[++num] = i;
mu[i] = -1;
}
for (int j = 1, cj; j <= num && (cj = prime[j] * i) <= 1000000; ++j) {
notp[cj] = true;
if (i % prime[j] == 0) {
mu[cj] = 0;
break;
}
mu[cj] = -mu[i];
}
}
for (int i = 0; i <= 1000000; ++i) g[i] = nig[i] = 1;
}
int Fib[N], niFib[N];
int solve(int n, int m) {
int ret = 1;
for (int ri, T = 1; T <= n; ++T) {
ri = min(n / (n / T), m / (m / T));
ret = 1ll * ret * ipow(1ll * g[ri] * nig[T - 1] % p, 1ll * (n / T) * (m / T) % (p - 1)) % p;
T = ri;
}
return ret;
}
int main() {
freopen("product.in", "r", stdin);
freopen("product.out", "w", stdout);
Euler_shai();
Fib[0] = 0; Fib[1] = 1; niFib[1] = niFib[0] = 1;
for (int i = 2; i <= 1000000; ++i) {
Fib[i] = (Fib[i - 1] + Fib[i - 2]) % p;
niFib[i] = ipow(Fib[i], p - 2);
}
for (int i = 1; i <= 1000000; ++i)
for (int j = i; j <= 1000000; j += i)
if (mu[j / i] == 1) {
g[j] = 1ll * g[j] * Fib[i] % p;
nig[j] = 1ll * nig[j] * niFib[i] % p;
} else if (mu[j / i] == -1) {
g[j] = 1ll * g[j] * niFib[i] % p;
nig[j] = 1ll * nig[j] * Fib[i] % p;
}
for (int i = 2; i <= 1000000; ++i) {
g[i] = 1ll * g[i] * g[i - 1] % p;
nig[i] = 1ll * nig[i] * nig[i - 1] % p;
}
int T, n, m; scanf("%d", &T);
while (T--) {
scanf("%d%d", &n, &m);
if (n > m) n ^= m ^= n ^= m;
printf("%d
", solve(n, m));
}
return 0;
}
T2 树点涂色
题目大意
· 根为1的有根树上初始每个节点都有不同颜色，一条路径的权值为这条路径上点的不同颜色数。
· 修改：把x到根这条链上所有点都涂上一种之前从没出现过的颜色。
· 查询：链的权值和子树内一点到根的路径的权值的最大值。
· (nleq10^5,mleq10^5)
思路
· 把具有相同颜色的链看成重链，那么重链上相邻的点只可能是父子关系，不可能是兄弟关系。
· 线段树维护dfs序，每个节点维护这个区间内所有点到根路径的权值的最大值，链的权值可以两个端点和lca差分一下求出，子树内一点到根的权值的最大值直接在线段树上查询子树的dfs序区间内的最大值。
· 修改操作类似lct的access操作，每次修改x点，相当于把x到根的这条路径变成重链，同时要相应的断掉以前的一些重链。在重链和轻链变换的时候顺便在线段树上修改即可。
· 时间复杂度(O(mlog^2n))。
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int N = 100003;
int n, deep[N], pos[N];
void check_max(int &a, int b) {if (b > a) a = b;}
namespace SegmentTree {
int ma[N << 2], lazy[N << 2];
void BuildTree(int rt, int l, int r) {
if (l == r) {ma[rt] = deep[pos[l]]; return;}
int mid = (l + r) >> 1;
BuildTree(rt << 1, l, mid);
BuildTree(rt << 1 | 1, mid + 1, r);
ma[rt] = max(ma[rt << 1], ma[rt << 1 | 1]);
}
void pushdown(int rt) {
if (lazy[rt]) {
int num = lazy[rt];
lazy[rt << 1] += num;
ma[rt << 1] += num;
lazy[rt << 1 | 1] += num;
ma[rt << 1 | 1] += num;
lazy[rt] = 0;
}
}
void add(int rt, int l, int r, int L, int R, int num) {
if (L <= l && r <= R) {lazy[rt] += num; ma[rt] += num; return;}
int mid = (l + r) >> 1;
pushdown(rt);
if (L <= mid) add(rt << 1, l, mid, L, R, num);
if (R > mid) add(rt << 1 | 1, mid + 1, r, L, R, num);
ma[rt] = max(ma[rt << 1], ma[rt << 1 | 1]);
}
int query(int rt, int l, int r, int L, int R) {
if (L <= l && r <= R) return ma[rt];
int mid = (l + r) >> 1, s = 0;
pushdown(rt);
if (L <= mid) check_max(s, query(rt << 1, l, mid, L, R));
if (R > mid) check_max(s, query(rt << 1 | 1, mid + 1, r, L, R));
return s;
}
}
struct Enode {int nxt, to;} E[N << 1];
int cnt = 0, point[N], st[N], top;
void ins(int u, int v) {E[++cnt] = (Enode) {point[u], v}; point[u] = cnt;}
int m, L[N], R[N], dfsID = 0, sz[N], fath[N];
bool used[N];
void dfs() {
st[top = 1] = 1; deep[1] = 1;
while (top) {
int u = st[top];
if (!used[u]) {
L[u] = ++dfsID; pos[dfsID] = u;
for (int i = point[u]; i; i = E[i].nxt)
if (E[i].to != fath[u]) {
st[++top] = E[i].to;
fath[E[i].to] = u;
deep[E[i].to] = deep[u] + 1;
}
used[u] = true;
} else {
for (int i = point[u]; i; i = E[i].nxt)
if (E[i].to != fath[u])
sz[u] += sz[E[i].to];
R[u] = L[u] + sz[u];
++sz[u];
--top;
}
}
}
namespace LCT {
struct node *null;
struct node {
node *ch[2], *fa; int id, mn;
void setc(node *r, int c) {ch[c] = r; if (r != null) r->fa = this;}
bool pl() {return fa->ch[1] == this;}
bool check() {return fa == null || (fa->ch[0] != this && fa->ch[1] != this);}
void count() {
mn = id;
if (ch[0] != null) if (L[ch[0]->mn] < L[mn]) mn = ch[0]->mn;
if (ch[1] != null) if (L[ch[1]->mn] < L[mn]) mn = ch[1]->mn;
}
} pool[N];
void init() {
null = &pool[0]; null->ch[0] = null->ch[1] = null->fa = null;
node *t;
for (int i = 1; i <= n; ++i) {
t = pool + i; t->id = t->mn = i;
t->ch[0] = t->ch[1] = null;
t->fa = &pool[fath[i]];
}
}
void rotate(node *r) {
node *f = r->fa;
int c = r->pl();
if (f->check()) r->fa = f->fa;
else f->fa->setc(r, f->pl());
f->setc(r->ch[c ^ 1], c);
r->setc(f, c ^ 1);
f->count();
}
void splay(node *r) {
for (; !r->check(); rotate(r))
if (!r->fa->check()) rotate(r->pl() == r->fa->pl() ? r->fa : r);
r->count();
}
void access(node *r) {
node *y = null;
while (r != null) {
splay(r);
if (r->ch[1] != null) SegmentTree::add(1, 1, n, L[r->ch[1]->mn], R[r->ch[1]->mn], 1);
r->ch[1] = y; r->count();
if (y != null) SegmentTree::add(1, 1, n, L[y->mn], R[y->mn], -1);
y = r; r = r->fa;
}
}
}
int F[N][18];
int work(int x, int y) {
if (deep[x] < deep[y]) x ^= y ^= x ^= y;
int xx = x, yy = y;
for (int i = 17; i >= 0; --i)
if (deep[F[x][i]] >= deep[y])
x = F[x][i];
if (x == y) return SegmentTree::query(1, 1, n, L[xx], L[xx]) - SegmentTree::query(1, 1, n, L[x], L[x]) + 1;
for (int i = 17; i >= 0; --i)
if (F[x][i] != F[y][i])
x = F[x][i], y = F[y][i];
int numx, numy, numz;
int ret = SegmentTree::query(1, 1, n, L[xx], L[xx]) - (numx = SegmentTree::query(1, 1, n, L[x], L[x])) + 1
+ SegmentTree::query(1, 1, n, L[yy], L[yy]) - (numy = SegmentTree::query(1, 1, n, L[y], L[y])) + 1;
numz = SegmentTree::query(1, 1, n, L[fath[x]], L[fath[x]]);
if (numz != numx && numz != numy) ++ret;
return ret;
}
int main() {
freopen("paint.in", "r", stdin);
freopen("paint.out", "w", stdout);
scanf("%d%d", &n, &m);
int u, v;
for (int i = 1; i < n; ++i) {
scanf("%d%d", &u, &v);
ins(u, v); ins(v, u);
}
dfs();
SegmentTree::BuildTree(1, 1, n);
LCT::init();
for (int i = 1; i <= n; ++i) F[i][0] = fath[i];
for (int j = 1; j <= 17; ++j)
for (int i = 1; i <= n; ++i)
F[i][j] = F[F[i][j - 1]][j - 1];
int x, y, op;
while (m--) {
scanf("%d%d", &op, &x);
if (op == 2) scanf("%d", &y);
switch (op) {
case 1:
LCT::access(&LCT::pool[x]);
break;
case 2:
printf("%d
", work(x, y));
break;
case 3:
printf("%d
", SegmentTree::query(1, 1, n, L[x], R[x]));
break;
}
// for (int i = 1; i <= n; ++i) printf("%d ", SegmentTree::query(1, 1, n, L[i], L[i])); puts("");
}
return 0;
}
T3 序列计数
题目大意
· 长度为(n)的序列，满足序列中的数都是不大于(m)的正整数，(n)个数的和是(p)的倍数，而且这(n)个数中至少有一个是质数。
· 求满足条件的序列的个数。
· (1leq nleq 10^9,1leq m leq 2 imes 10^7,1leq p leq 100)
思路
· 容斥一下，在满足前两个条件的情况下，统计(n)个数无质数限制的序列个数再减去(n)个数都不是质数的序列个数。
· 设(f(i,j))表示长度为(i)，和模(p)为(j)的序列个数。
· 提前预处理模(p)为(j)有多少个数，即(f(1,j))。欧拉筛求出(m)以内的质数。
· (f(i,j)=sumlimits_{k,l}[(k+l)mod p=j]f(i-1,k) imes f(1,l))
· dp可以用倍增优化，时间复杂度(Oleft(m+p^2log n
ight))。
· 当然矩乘也是可以的。
#include<cstdio>
#include<bitset>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N = 20000003;
const int mo = 20170408;
int f[103], g[103], F[103], n, m, p;
int nump = 0, prime[1270608];
bitset <N> notp;
void Euler_shai() {
notp[1] = 1;
for (int i = 2; i <= m; ++i) {
if (!notp[i]) prime[++nump] = i;
for (int j = 1; j <= nump && prime[j] * i <= m; ++j) {
notp[prime[j] * i] = 1;
if (i % prime[j] == 0) break;
}
}
}
int main() {
freopen("count.in", "r", stdin);
freopen("count.out", "w", stdout);
scanf("%d%d%d", &n, &m, &p);
int num = n, tot = 0;
while (num) {num >>= 1; ++tot;}
for (int i = 1; i <= m; ++i) ++f[i % p];
for (int i = 0; i < p; ++i) F[i] = (f[i] %= mo);
for (int i = tot - 2; i >= 0; --i) {
memset(g, 0, p << 2);
for (int j = 0; j < p; ++j)
for (int k = 0; k < p; ++k)
(g[(j + k) % p] += 1ll * f[j] * f[k] % mo) %= mo;
memcpy(f, g, p << 2);
if ((n >> i) & 1) {
memset(g, 0, p << 2);
for (int j = 0; j < p; ++j)
for (int k = 0; k < p; ++k)
(g[(j + k) % p] += 1ll * f[j] * F[k] % mo) %= mo;
memcpy(f, g, p << 2);
}
}
int ans = f[0];
Euler_shai();
memset(f, 0, p << 2); memset(F, 0, p << 2);
for (int i = 1; i <= m; ++i) if (notp[i]) ++f[i % p];
for (int i = 0; i < p; ++i) F[i] = (f[i] %= mo);
for (int i = tot - 2; i >= 0; --i) {
memset(g, 0, p << 2);
for (int j = 0; j < p; ++j)
for (int k = 0; k < p; ++k)
(g[(j + k) % p] += 1ll * f[j] * f[k] % mo) %= mo;
memcpy(f, g, p << 2);
if ((n >> i) & 1) {
memset(g, 0, p << 2);
for (int j = 0; j < p; ++j)
for (int k = 0; k < p; ++k)
(g[(j + k) % p] += 1ll * f[j] * F[k] % mo) %= mo;
memcpy(f, g, p << 2);
}
}
printf("%d
", (ans - f[0] + mo) % mo);
return 0;
}
Day 2
T1 新生舞会
题目大意
· (n)个男生和(n)个女生组成舞伴，第(i)个男生和第(j)个女生组成舞伴会产生(a_{ij})的喜悦程度和(b_{ij})的不协调程度。
· 有(n)个喜悦程度(a_1',a_2',dots a_n')和不协调程度(b_1',b_2',dots b_n')，最大化(C=frac{a_1'+a_2'+dots+a_n'}{b_1'+b_2'+dots+b_n'})
· (1leq nleq 100,1leq a_{ij},b_{ij}leq 10^4)，保留6位小数。
思路
· 分数规划，把分母乘到左边， 左边再移向到右边，二分(C)的值，通过移项后的式子的最大值来判断(C)偏大还是偏小。
· 移项后的式子的最大值可以用最大费用最大流来计算。
· 用数组存边比结构体快一倍！！！
· 把实数乘(10^7)变成整数可以减小常数。
· 时间复杂度(Oleft(n^2sqrt nlog 10^{12}
ight))。
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int in() {
int k = 0; char c = getchar();
for (; c < '0' || c > '9'; c = getchar());
for (; c >= '0' && c <= '9'; c = getchar())
k = k * 10 + (c ^ 48);
return k;
}
const int N = 103;
int b[N][N], n, S, T;
ll a[N][N], w[N * N * 3];
int nxt[N * N * 3], to[N * N * 3], c[N * N * 3], cnt, point[N << 1];
void ins(int u, int v, int C, ll W) {
++cnt; nxt[cnt] = point[u]; to[cnt] = v; c[cnt] = C; w[cnt] = W; point[u] = cnt;
++cnt; nxt[cnt] = point[v]; to[cnt] = u; c[cnt] = 0; w[cnt] = -W; point[v] = cnt;
}
bool inq[N << 1];
ll dis[N << 1];
int qu[100000], pre[N << 1];
bool spfa() {
memset(dis, 0xc0, (T + 1) << 3);
int p = 0, q = 1, v;
ll t;
qu[1] = S; dis[S] = 0; inq[S] = true;
while (p != q) {
++p; if (p == 100000) p = 0;
int u = qu[p]; inq[u] = false;
for (int i = point[u]; i; i = nxt[i])
if (c[i] && (t = dis[u] + w[i]) > dis[v = to[i]]) {
pre[v] = i;
dis[v] = t;
if (!inq[v]) {
++q; if (q == 100000) q = 0;
qu[q] = v;
inq[v] = true;
}
}
}
return dis[T] != dis[0];
}
ll mcmf(ll C) {
cnt = 1; memset(point, 0, (T + 1) << 2);
for (int i = 1; i <= n; ++i) {
ins(S, i, 1, 0), ins(i + n, T, 1, 0);
ll *ta = a[i]; int *tb = b[i];
for (int j = 1; j <= n; ++j) {
++ta; ++tb;
ins(i, j + n, 1, *ta - C * *tb);
}
}
ll f = 0;
while (spfa()) {
int u = T, e;
for (u = T, e = pre[u]; u != S; e = pre[u = to[e ^ 1]]) --c[e], ++c[e ^ 1];
f += dis[T];
}
return f;
}
int main() {
freopen("ball.in", "r", stdin);
freopen("ball.out", "w", stdout);
n = in();
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= n; ++j)
a[i][j] = in(), a[i][j] *= 10000000ll;
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= n; ++j)
b[i][j] = in();
S = (n << 1) + 1; T = S + 1;
ll left, right, s1 = 0, s2 = 0;
for (int i = 1; i <= n; ++i)
s1 += a[i][i], s2 += b[i][i];
left = s1 / s2 + 1;
s1 = 0; s2 = 0;
ll mx;
for (int i = 1; i <= n; ++i) {
mx = 0;
for (int j = 1; j <= n; ++j)
if (a[i][j] > mx)
mx = a[i][j];
s1 += mx;
mx = 0x7fffffff;
for (int j = 1; j <= n; ++j)
if (b[i][j] < mx)
mx = b[i][j];
s2 += mx;
}
right = s1 / s2;
ll mid;
while (left < right) {
mid = (left + right + 1) >> 1;
if (mcmf(mid) >= 0) left = mid;
else right = mid - 1;
}
if (left % 10 >= 5) left += 10; left /= 10;
printf("%I64d.%06I64d
", left / 1000000, left % 1000000);
return 0;
}
T2 硬币游戏
题目大意
· (n)个同学每个同学猜一个长度为(m)的HT序列（H表示正面朝上，T表示反面朝上）
· 不断地扔硬币，组成一个HT序列。当某个同学的猜的序列突然出现在扔出来的硬币序列中，则停止扔硬币并且这个同学获胜。
· 问每个同学获胜的概率是多少。
· (1leq n,mleq 300)
思路
· 这道题我想了一整天，最后问了zyf2000才知道怎么做。
· 某个同学猜的序列在硬币序列中出现的位置肯定是最后。
· 首先设(f(i,j))表示长度为(i)的硬币序列，只有第(i)个位置和(m)个序列中某一个序列匹配，并且这个匹配的序列是第(j)个序列的概率。
· 那么第(i)个同学获胜的概率就是(sumlimits_{jgeq 1}f(j,i))
· 设(N)表示没有扔出来和同学猜的序列匹配的硬币序列的概率，这个(N)有点抽象，更具体的：(N=sumlimits_{jgeq 1}left(1-sumlimits_{k=1}^nf(j,k)
ight))
· 然后假设A同学猜的序列是TTH，B同学猜的序列是HTT，那么可以构造一个以TTH结尾的硬币序列，并且1到序列长度-3的所有位置都不是能匹配A或B的位置，所以匹配位置只可能是序列长度-2
，序列长度-1或序列长度（即最后一个位置）
· 这样就有一个等式了：P(NTTH)=P(A)+P(BTH)+P(BH)
· 很明显吧，如果一个序列的前缀是另一个序列的后缀，那么就会产生贡献。
· 所有同学获胜的概率和为1，(n+1)个变元，(n+1)个等式，高斯消元即可。
· 时间复杂度(O(n^3))。
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int N = 703;
int cnt;
double a[N][N], f[N];
void gauss() {
for (int i = 1; i <= cnt; ++i) {
int tmp = i;
for (int j = i + 1; j <= cnt; ++j)
if (fabs(a[j][i]) > fabs(a[tmp][i]))
tmp = j;
if (tmp != i)
for (int j = i; j <= cnt + 1; ++j)
swap(a[i][j], a[tmp][j]);
for (int j = i + 1; j <= cnt; ++j)
if (a[j][i] != 0) {
double k = a[j][i] / a[i][i];
for (int t = i; t <= cnt + 1; ++t)
a[j][t] -= a[i][t] * k;
}
}
for (int i = cnt; i >= 1; --i) {
for (int j = i + 1; j <= cnt; ++j)
a[i][cnt + 1] -= a[i][j] * f[j];
f[i] = a[i][cnt + 1] / a[i][i];
}
}
int end[N];
double fac[N];
int n, m, Acnt = 1, ch[N * N][2], fail[N * N], qu[N * N], deep[N * N], to[N * N];
char s[N];
void insert(int id) {
int tmp = 1;
for (int i = 1; i <= m; ++i) {
int v = s[i] == 'H' ? 0 : 1;
if (ch[tmp][v] == 0) ch[tmp][v] = ++Acnt;
tmp = ch[tmp][v];
}
end[id] = tmp; to[tmp] = id;
}
int dfsID = 0, L[N * N], R[N * N], dfs_a[N * N];
void dfs(int x) {
L[x] = ++dfsID;
if (to[x])
dfs_a[dfsID] = to[x];
else {
if (ch[x][0]) dfs(ch[x][0]);
if (ch[x][1]) dfs(ch[x][1]);
}
R[x] = dfsID;
}
void BuildFail() {
int p = 0, q = 1; qu[1] = 1;
while (p != q) {
int u = qu[++p];
for (int i = 0; i <= 1; ++i)
if (ch[u][i]) {
int p = fail[u];
while (p && ch[p][i] == 0) p = fail[p];
fail[ch[u][i]] = p ? ch[p][i] : 1;
deep[ch[u][i]] = deep[u] + 1;
qu[++q] = ch[u][i];
}
}
dfs(1);
}
void Jump(int id) {
int p = fail[end[id]];
while (p != 1) {
for (int i = L[p]; i <= R[p]; ++i)
if (dfs_a[i])
a[dfs_a[i]][id] += fac[m - deep[p]];
p = fail[p];
}
}
int main() {
freopen("game.in", "r", stdin);
freopen("game.out", "w", stdout);
scanf("%d%d", &n, &m);
for (int i = 1; i <= n; ++i) {
scanf("%s", s + 1);
insert(i);
}
fac[0] = 1;
for (int i = 1; i <= m; ++i) fac[i] = fac[i - 1] / 2;
BuildFail();
cnt = n + 1;
for (int i = 1; i <= n; ++i) {
a[i][i] = 1; a[i][cnt] = -1;
Jump(i);
}
for (int i = 1; i <= n; ++i) a[cnt][i] = 1;
a[cnt][cnt + 1] = 1;
gauss();
for (int i = 1; i <= n; ++i) printf("%.10lf
", f[i]);
return 0;
}
T3 相关分析
题目大意
· (n)个二元组((x,y))构成一个序列。
· 修改([L,R],S,T)：区间内所有二元组的(x)都加上(S)，(y)都加上(T)。
· 修改([L,R],S,T)：区间内所有二元组的(x)都变为(S+i)，(y)都变为(T+i)，(i)为二元组的下标。
· 查询([L,R])：计算(frac{sumlimits_{i=L}^R(x_i-ar x)(y_i-ar y)}{sumlimits_{i=L}^R(x_i-ar x)^2})
· (1leq n,mleq10^5,1leq Lleq Rleq n,0leq|S|,|T|leq10^5,0leq|x_i|,|y_i|leq10^5)
思路
· 线段树维护4个信息:(sum x_i,sum y_i,sum x_i^2,sum x_iy_i)
· 把要计算的式子全拆开，所需要的信息就是上面的4个信息。
· 区间加：(sum (x_i+S)^2=sum(x_i^2+2Sx_i+S^2),sum (x_i+S)(y_i+T)=sum(x_iy_i+Tx_i+Sy_i+ST))
· 区间覆盖：(sum x_i^2)可以用平方和公式(frac {n(n+1)(2n+1)}6)计算（直接预处理平方和可以减小常数），(sum x_iy_i)同理。
· 维护一个区间加标记和区间覆盖标记就可以了。
· 时间复杂度(O(mlog n))。
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int in() {
int k = 0, fh = 1; char c = getchar();
for (; c < '0' || c > '9'; c = getchar())
if (c == '-') fh = -1;
for (; c >= '0' && c <= '9'; c = getchar())
k = k * 10 + (c ^ 48);
return k * fh;
}
const int N = 100003;
int n, m, len;
double x[N], y[N], xl, yl;
struct PP {
double sumx, sumy, sumxx, sumxy;
PP (double _sumx = 0, double _sumy = 0, double _sumxx = 0, double _sumxy = 0)
: sumx(_sumx), sumy(_sumy), sumxx(_sumxx), sumxy(_sumxy) {}
PP operator + (const PP &A) const {
return PP(sumx + A.sumx, sumy + A.sumy, sumxx + A.sumxx, sumxy + A.sumxy);
}
};
namespace ST {
double sum_x[N << 2], sum_y[N << 2], sum_xy[N << 2], sum_xx[N << 2], lazy_x[N << 2], lazy_y[N << 2];
double cover_x[N << 2], cover_y[N << 2];
bool cover_mark[N << 2];
ll getsum(int l, int r) {
return (1ll * r * (r + 1) * ((r << 1) + 1) - 1ll * (l - 1) * l * (((l - 1) << 1) + 1)) / 6;
}
void coverit(int rt, int l, int r, const double &S, const double &T) {
lazy_x[rt] = lazy_y[rt] = 0;
cover_x[rt] = S; cover_y[rt] = T; cover_mark[rt] = true;
len = r - l + 1;
sum_x[rt] = (S + l + S + r) * len / 2;
sum_y[rt] = (T + l + T + r) * len / 2;
sum_xx[rt] = getsum(S + l, S + r);
sum_xy[rt] = S * T * len + (S + T) * (1ll * (l + r) * len / 2) + getsum(l, r);
}
double Slen, Tlen;
void update(int rt, int l, int r, const double &S, const double &T) {
if (cover_mark[rt]) {
if (l != r) {
int mid = (l + r) >> 1;
coverit(rt << 1, l, mid, cover_x[rt], cover_y[rt]);
coverit(rt << 1 | 1, mid + 1, r, cover_x[rt], cover_y[rt]);
}
cover_mark[rt] = false;
}
lazy_x[rt] += S; lazy_y[rt] += T;
len = r - l + 1; Slen = S * len; Tlen = T * len;
sum_xx[rt] += S * 2 * sum_x[rt] + S * Slen;
sum_xy[rt] += sum_x[rt] * T + sum_y[rt] * S + S * Tlen;
sum_x[rt] += Slen;
sum_y[rt] += Tlen;
}
void pushdown(int rt, int l, int r) {
if (lazy_x[rt] || lazy_y[rt]) {
int mid = (l + r) >> 1;
update(rt << 1, l, mid, lazy_x[rt], lazy_y[rt]);
update(rt << 1 | 1, mid + 1, r, lazy_x[rt], lazy_y[rt]);
lazy_x[rt] = lazy_y[rt] = 0;
}
if (cover_mark[rt]) {
int mid = (l + r) >> 1;
coverit(rt << 1, l, mid, cover_x[rt], cover_y[rt]);
coverit(rt << 1 | 1, mid + 1,r, cover_x[rt], cover_y[rt]);
cover_mark[rt] = false;
}
}
void pushup(int rt) {
sum_x[rt] = sum_x[rt << 1] + sum_x[rt << 1 | 1];
sum_y[rt] = sum_y[rt << 1] + sum_y[rt << 1 | 1];
sum_xx[rt] = sum_xx[rt << 1] + sum_xx[rt << 1 | 1];
sum_xy[rt] = sum_xy[rt << 1] + sum_xy[rt << 1 | 1];
}
void build(int rt, int l, int r) {
if (l == r) {
xl = x[l]; yl = y[l];
sum_x[rt] = xl; sum_y[rt] = yl;
sum_xx[rt] = xl * xl;
sum_xy[rt] = xl * yl;
cover_mark[rt] = false;
return;
}
int mid = (l + r) >> 1;
build(rt << 1, l, mid);
build(rt << 1 | 1, mid + 1, r);
pushup(rt);
}
void add(int rt, int l, int r, int L, int R, const double &S, const double &T) {
if (L <= l && r <= R) {
update(rt, l, r, S, T);
return;
}
int mid = (l + r) >> 1;
pushdown(rt, l, r);
if (L <= mid) add(rt << 1, l, mid, L, R, S, T);
if (R > mid) add(rt << 1 | 1, mid + 1, r, L, R, S, T);
pushup(rt);
}
void cover(int rt, int l, int r, int L, int R, const double &S, const double &T) {
if (L <= l && r <= R) {
coverit(rt, l, r, S, T);
return;
}
int mid = (l + r) >> 1;
pushdown(rt, l, r);
if (L <= mid) cover(rt << 1, l, mid, L, R, S, T);
if (R > mid) cover(rt << 1 | 1, mid + 1, r, L, R, S, T);
pushup(rt);
}
PP getans(int rt, int l, int r, int L, int R) {
if (L <= l && r <= R) return PP(sum_x[rt], sum_y[rt], sum_xx[rt], sum_xy[rt]);
PP tl, tr;
int mid = (l + r) >> 1;
pushdown(rt, l, r);
if (L <= mid) tl = getans(rt << 1, l, mid, L, R);
if (R > mid) tr = getans(rt << 1 | 1, mid + 1, r, L, R);
return tl + tr;
}
}
int main() {
freopen("relative.in", "r", stdin);
freopen("relative.out", "w", stdout);
n = in(); m = in(); //scanf("%d%d", &n, &m);
for (int i = 1; i <= n; ++i) scanf("%lf", x + i);
for (int i = 1; i <= n; ++i) scanf("%lf", y + i);
ST::build(1, 1, n);
int L, R, op; double S, T, xb, yb, fz, fm;
PP r;
while (m--) {
op = in(); //scanf("%d", &op);
if (op == 1) {
L = in(); R = in();//scanf("%d%d", &L, &R);
r = ST::getans(1, 1, n, L, R);
xb = r.sumx / (R - L + 1);
yb = r.sumy / (R - L + 1);
fz = r.sumxy - r.sumy * xb - r.sumx * yb + xb * yb * (R - L + 1);
fm = r.sumxx - r.sumx * 2 * xb + xb * xb * (R - L + 1);
printf("%.10lf
", fz / fm);
} else if (op == 2) {
L = in(); R = in();
scanf("%lf%lf", &S, &T);
ST::add(1, 1, n, L, R, S, T);
} else {
L = in(); R = in();
scanf("%lf%lf", &S, &T);
ST::cover(1, 1, n, L, R, S, T);
}
/*
for (int i = 1; i <= n; ++i) {
PP r = ST::getans(1, 1, n, i, i);
printf("(%.3lf, %.3lf) ", r.sumx, r.sumy);
}
puts("");*/
}
return 0;
}