[SWERC2015]Saint John Festival

第一次做像点样子的计算几何吧……

很显然要包含在某个三角形里,就要被包含在大点的凸包里。所以求出大点组成的凸包,然后对于所有小点判断它是否在那个凸包里即可。

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <deque>
#include <cmath>
#include <vector>
using R = double;
const R eps = 1e-8;
int sign(R x) {
  if(fabs(x) < eps) {
    return 0;
  } else {
    if(x < 0) return -1;
    else return 1;
  }
}
struct Point {
  R x, y;
  Point(R qx = 0, R qy = 0) {
    x = qx; y = qy;
  }
};
typedef Point Vector;
Vector operator +(const Vector &a, const Vector &b) {
  return Vector(a.x + b.x, a.y + b.y);
}
Vector operator -(const Point &a, const Point &b) {
  return Vector(a.x - b.x, a.y - b.y);
}
R dot(const Vector &a, const Vector &b) {
  return a.x * b.x + a.y * b.y;
}
R times(const Vector &a, const Vector &b) {
  return a.x * b.y - a.y * b.x;
}
bool cmp(const Point &a, const Point &b) {
  if(sign(a.x - b.x) == 0) {
    return a.y < b.y;
  } else {
    return a.x < b.x;
  }
}
const int maxn = 10005;
Point P[maxn]; int L;
std::vector<Point> bot, top;
void andrew() {
  // bot.resize(L); top.resize(L);
  std::sort(P + 1, P + 1 + L, cmp);
  for(int i = 1; i <= L; i ++) {
    if(i != 1 && sign(P[i].x - P[i - 1].x) == 0) continue;
    while(bot.size() > 1 && sign(times(P[i] - bot.back(), bot.back() - bot[bot.size() - 2])) >= 0) {
      bot.pop_back();
    }
    bot.push_back(P[i]);
  }
  for(int i = L; i >= 1; i --) {
    if(i != L && sign(P[i].x - P[i + 1].x) == 0) continue;
    while(top.size() > 1 && sign(times(P[i] - top.back(), top.back() - top[top.size() - 2])) >= 0) {
      top.pop_back();
    }
    top.push_back(P[i]);
  }
  std::reverse(top.begin(), top.end());
#ifdef LOCAL
  printf("top :");
  for(auto p : top) {
    printf(" (%.3f, %.3f)", p.x, p.y);
  }
  puts("");
  printf("bot :");
  for(auto p : bot) {
    printf(" (%.3f, %.3f)", p.x, p.y);
  }
  puts("");
#endif
}
R get_p(Point l, Point r, R x) {
  R delta = x - l.x;
  R k = (r.y - l.y) / (r.x - l.x);
#ifdef LOCAL
  printf("k between (%.2lf, %.2lf) and (%.2lf, %.2lf) : %.5lf\n", l.x, l.y, r.x, r.y, k);
#endif
  return l.y + k * delta;
}
bool cmp2(const Point &a, const Point &b) {
  return sign(a.x - b.x) == -1;
}
bool query(const Point &p) {
  if(sign(top.front().x - p.x) == 1 || sign(p.x - top.back().x) == 1) return false;
  auto lp = (std::lower_bound(bot.begin(), bot.end(), p, cmp2));
  auto rp = (std::lower_bound(top.begin(), top.end(), p, cmp2));
  R l;
  if(sign((*lp).x - p.x) == 0) {
    l = (*lp).y;
  } else {
    auto lp2 = lp - 1;
    l = get_p(*lp2, *lp, p.x);
  }
  R r;
  if(sign((*rp).x - p.x) == 0) {
    r = (*rp).y;
  } else {
    auto rp2 = rp - 1;
    r = get_p(*rp2, *rp, p.x);
  }
#ifdef LOCAL
  printf("bd of (%.2lf, %.2lf) : [%.2lf, %.2lf]\n", p.x, p.y, l, r);
#endif
  return (sign(l - p.y) <= 0 && sign(p.y - r) <= 0);
}

int main() {
  scanf("%d", &L);
  for(int i = 1; i <= L; i ++) {
    scanf("%lf%lf", &P[i].x, &P[i].y);
  }
  andrew();
  int S; scanf("%d", &S);
  int cnt = 0;
  while(S --) {
    R x, y; scanf("%lf%lf", &x, &y);
    Point p(x, y);
    if(query(p)) cnt ++;
  }
  printf("%d\n", cnt);
  return 0;
}

[BZOJ 1070][SCOI2007]修车

啊啊啊我为什么要去颓费用流啊……

这个题可以保证每个人的车最后都会被修,所以说也就是让你最小化总等待时间辣!

直接考虑给一个技术人员建时间轴是麻烦而不可取的。我们考虑一点,对于每一个技术人员,若总共修了\(t\)辆车,那么第\(i\)个来找他修车的人会影响后面人的等待时间,换言之我们可以认为第\(i\)个来修车的人给答案做出了\(c\cdot (t - i + 1)\)(这里用\(c\)表示他自己的修车时间,注意他也会影响他自己)的贡献,似乎可以苟最小费用最大流呢。

但问题在于我们并不能提前知道\(t\)是多少,并且这样的话费用不是递增的,一个人可能会优先考虑到后面的位置去修车。那么我们考虑倒过来,去考虑倒数第\(i\)个来修车的人的贡献,显然这个贡献是\(c\cdot i\)。然后接下来就肥肠简单了,,,

代码:

/**************************************************************
	Problem: 1070
	User: danihao123
	Language: C++
	Result: Accepted
	Time:576 ms
	Memory:13928 kb
****************************************************************/

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <utility>
#include <algorithm>
#include <queue>
#include <cassert>
typedef long long ll;
int m, n;
const int maxno = 100005;
const int maxe = 400005;
int first[maxno];
int next[maxe], from[maxe], to[maxe], flow[maxe], cap[maxe];
ll cost[maxe];
int gcnt = 0;
void add_edge(int u, int v, int f, ll c) {
  gcnt ++;
  next[gcnt] = first[u];
  first[u] = gcnt;
  from[gcnt] = u; to[gcnt] = v;
  cap[gcnt] = f; flow[gcnt] = 0;
  cost[gcnt] = c;
}
int rev(int i) {
  return ((i - 1) ^ 1) + 1;
}
void ins_edge(int u, int v, int f, ll c = 0LL) {
#ifdef LOCAL
  printf("Inserting Edge (%d, %d, %d, %lld)\n", u, v, f, c);
#endif
  add_edge(u, v, f, c);
  add_edge(v, u, 0, -c);
}

const ll LINF = 0x7f7f7f7f7f7f7f7fLL;
int tot;
bool spfa(int s, int t, int &f, ll &c) {
  static ll d[maxno];
  static bool inq[maxno];
  static int a[maxno], p[maxno];
  std::fill(d, d + tot + 1, LINF);
  std::fill(inq, inq + tot + 1, false);
  std::fill(a, a + tot + 1, 0);
  std::fill(p, p + tot + 1, 0);
  d[s] = 0;
  std::queue<int> Q; Q.push(s); inq[s] = true;
  a[s] = 0x7fffffff;
  while(!Q.empty()) {
    int u = Q.front(); Q.pop();
    inq[u] = false;
    for(int i = first[u]; i; i = next[i]) {
      if(cap[i] > flow[i]) {
        int v = to[i];
        if(d[v] > d[u] + cost[i]) {
          d[v] = d[u] + cost[i];
          a[v] = std::min(a[u], cap[i] - flow[i]); p[v] = i;
          if(!inq[v]) {
            Q.push(v); inq[v] = true;
          }
        }
      }
    }
  }
  if(a[t] == 0) return false;
  f += a[t];
  c += (ll(a[t])) * d[t];
  for(int e = p[t]; e; e = p[from[e]]) {
    flow[e] += a[t];
    flow[rev(e)] -= a[t];
  }
  return true;
}
void mcmf(int s, int t, int &f, ll &c) {
  while(spfa(s, t, f, c));
}

int pos[105][105][2];
int main() {
  scanf("%d%d", &m, &n);
  int s = 0, t = 1; tot = 1;
  for(int i = 1; i <= n; i ++) {
    ins_edge(s, ++ tot, 1);
  }
  for(int i = 1; i <= m; i ++) {
    for(int j = 1; j <= n; j ++) {
      pos[i][j][0] = ++ tot;
      pos[i][j][1] = ++ tot;
      ins_edge(tot - 1, tot, 1);
      ins_edge(tot, t, 1);
    }
  }
  for(int i = 1; i <= n; i ++) {
    for(int j = 1; j <= m; j ++) {
      int T; scanf("%d", &T);
      for(int k = 1; k <= n; k ++) {
        int id = pos[j][k][0]; int id2 = pos[j][k][1];
        ins_edge(i + 1, id, 1, k * T);
      }
    }
  }
  int f = 0; ll c = 0;
  mcmf(s, t, f, c);
#ifdef LOCAL
  printf("%d\n", f);
#endif
  printf("%.2lf\n", (double(c)) / (double(n)));
  return 0;
}

[BZOJ 1449][JSOI2009]球队收益

二次函数费用流的入门题……我真太弱了……

可以给比赛、队伍都建点,然后\(S\)向每个比赛连容量为1的边,每个比赛向队伍连容量为\(1\)的边,来表示赢家是谁。这样一来一次比赛只有一个队伍赢了。

考虑怎么处理那个二次函数费用。费用函数为\(f(x, y) = C x^2 + D y^2\),由于一个队伍的总比赛次数是已知的,所以\(y\)不需要,不妨假设一个队伍有\(t\)场比赛,则费用函数为\(f(x) = C x^2 + D(t - x)^2\)。

对这个函数做差分:\(\Delta f(x) = f(x + 1) - f(x) = 2Cx - 2D(t - x) + C + D\),然后这个差分是单调不降的。因此我们对所有差分都从队伍向\(T\)连边(费用为这一点的差分),这样一来的话会优先选择费用小的边,保证这个队伍的费用一直合法。

代码:

/**************************************************************
	Problem: 1449
	User: danihao123
	Language: C++
	Result: Accepted
	Time:732 ms
	Memory:1448 kb
****************************************************************/

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <queue>
#include <cassert>
typedef long long ll;
const int maxn = 5005;
const int maxm = 1005;
const int maxno = maxn + maxm + 5;
const int maxe = (1000 * 2 + 1000 * 3) * 2 + 50;
int first[maxno];
int next[maxe], from[maxe], to[maxe], flow[maxe], cap[maxe];
ll cost[maxe];
int gcnt = 0;
void add_edge(int u, int v, int f, ll c) {
  gcnt ++;
  next[gcnt] = first[u];
  first[u] = gcnt;
  from[gcnt] = u; to[gcnt] = v;
  cap[gcnt] = f; flow[gcnt] = 0;
  cost[gcnt] = c;
}
int rev(int i) {
  return ((i - 1) ^ 1) + 1;
}
void ins_edge(int u, int v, int f, ll c) {
  add_edge(u, v, f, c);
  add_edge(v, u, 0, -c);
}

int n, m;
const ll LINF = 0x7f7f7f7f7f7f7f7fLL;
bool spfa(int s, int t, int &f, ll &c) {
  static ll d[maxno];
  static bool inq[maxno];
  static int a[maxno], p[maxno];
  std::fill(d, d + n + m + 2, LINF);
  std::fill(inq, inq + n + m + 2, false);
  std::fill(a, a + n + m + 2, 0);
  std::fill(p, p + n + m + 2, 0);
  d[s] = 0;
  std::queue<int> Q; Q.push(s); inq[s] = true;
  a[s] = 0x7fffffff;
  while(!Q.empty()) {
    int u = Q.front(); Q.pop();
    inq[u] = false;
    for(int i = first[u]; i; i = next[i]) {
      if(cap[i] > flow[i]) {
        int v = to[i];
        if(d[v] > d[u] + cost[i]) {
          d[v] = d[u] + cost[i];
          a[v] = std::min(a[u], cap[i] - flow[i]); p[v] = i;
          if(!inq[v]) {
            Q.push(v); inq[v] = true;
          }
        }
      }
    }
  }
  if(a[t] == 0) return false;
  f += a[t];
  c += (ll(a[t])) * d[t];
  for(int e = p[t]; e; e = p[from[e]]) {
    flow[e] += a[t];
    flow[rev(e)] -= a[t];
  }
  return true;
}
void mcmf(int s, int t, int &f, ll &c) {
  while(spfa(s, t, f, c));
}

ll win[maxn], lose[maxn], C[maxn], D[maxn];
ll tot[maxn];
int main() {
  scanf("%d%d", &n, &m);
  ll ans = 0;
  for(int i = 1; i <= n; i ++) {
    scanf("%lld%lld%lld%lld", &win[i], &lose[i], &C[i], &D[i]);
    tot[i] = win[i] + lose[i];
  }
  int s = 0, t = n + m + 1;
  for(int i = 1; i <= m; i ++) {
    int a, b; scanf("%d%d", &a, &b);
    ins_edge(0, i + n, 1, 0);
    ins_edge(i + n, a, 1, 0);
    ins_edge(i + n, b, 1, 0);
    tot[a] ++; tot[b] ++;
  }
  for(int i = 1; i <= n; i ++) {
    ans += C[i] * win[i] * win[i] + D[i] * (tot[i] - win[i]) * (tot[i] - win[i]);
    for(ll j = win[i]; j <= (tot[i] - lose[i] - 1LL); j ++) {
      ins_edge(i, t, 1, 2LL * C[i] * j - 2LL * D[i] * (tot[i] - j) + C[i] + D[i]);
    }
  }
  int f = 0; mcmf(s, t, f, ans);
  printf("%lld\n", ans);
  return 0;
}

[LibreOJ 2146][SHOI2017]寿司餐厅

建了半天图……搞出来了几个神奇的最小割……

然后发现这TM不就是最简单的最大权闭合子图吗……

首先和正负点权普通的最大权闭合子图一样,然后任意\(d_{i, i}\)去依赖点\(i\)。对于一个\(d_{i, j}(i < j)\),我们要保证那一天\([i, j]\)全部被吃了,所以说它需要依赖\(d_{i + 1, j}\)和\(d_{i, j - 1}\)。

每种寿司的费用也不难处理,我们把\(mx^2\)和\(cx\)分开考虑。\(mx^2\)单独建点,很显然代号为所有\(x\)都要依赖它。对于\(cx\),可以视为所有点有一个\(x\)的费用。

然后van了……

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <cassert>
#include <algorithm>
#include <utility>
#include <vector>
#include <queue>
const int maxn = 5000005;
const int maxm = 5000005;
int first[maxn];
int next[maxm], to[maxm], cap[maxm], flow[maxm];
int gcnt = 0;
void add_edge(int u, int v, int f) {
  gcnt ++;
  next[gcnt] = first[u];
  first[u] = gcnt;
  to[gcnt] = v;
  cap[gcnt] = f;
  flow[gcnt] = 0;
}
int rev(int i) {
  if(1 & i) {
    return i + 1;
  } else {
    return i - 1;
  }
}
void ins_edge(int u, int v, int f) {
  add_edge(u, v, f); add_edge(v, u, 0);
}

int s, t;
int d[maxn];
bool bfs() {
#ifdef LOCAL
  // puts("A new round :");
#endif
  static bool vis[maxn];
  memset(vis, 0, sizeof(vis));
  d[s] = 1; vis[s] = true;
  std::queue<int> Q; Q.push(s);
  while(!Q.empty()) {
    int u = Q.front(); Q.pop();
    for(int i = first[u]; i; i = next[i]) {
      int v = to[i];
      if(cap[i] > flow[i] && !vis[v]) {
        d[v] = d[u] + 1; vis[v] = true;
#ifdef LOCAL
        // printf("d[%d] : %d\n", v, d[v]);
#endif
        Q.push(v);
      }
    }
  }
  return vis[t];
}
int cur[maxn];
int dfs(int u, int a) {
#ifdef LOCAL
  // printf("State (%d, %d)\n", u, a);
#endif
  if(a == 0 || u == t) return a;
  int fl = 0;
  for(int &i = cur[u]; i; i = next[i]) {
    int v = to[i]; int f;
    if(d[v] == d[u] + 1 && (f = dfs(v, std::min(cap[i] - flow[i], a))) > 0) {
      fl += f; a -= f;
      flow[i] += f; flow[rev(i)] -= f;
      if(a == 0) break;
    }
  }
  if(a > 0) d[u] = -1;
  return fl;
}
int tot;
int dinic() {
  int ans = 0;
  while(bfs()) {
    for(int i = 0; i <= tot; i ++) cur[i] = first[i];
    ans += dfs(s, 0x7fffffff);
  }
  return ans;
}

int ma_p[105], ma_m[1005];
int ma_d[105][105];
int main() {
  tot = 1;
  int n, m; scanf("%d%d", &n, &m);
  int ans = 0;
  s = 0, t = 1;
  for(int i = 1; i <= n; i ++) {
    int a; scanf("%d", &a);
    ma_p[i] = ++ tot;
#ifdef LOCAL
    printf("p[%d] : %d\n", i, tot);
#endif
    if(!ma_m[a]) {
      ma_m[a] = ++ tot;
#ifdef LOCAL
      printf("m[%d] : %d\n", a, tot);
#endif
      ins_edge(tot, 1, m * a * a);
    }
    ins_edge(ma_p[i], 1, a);
    ins_edge(ma_p[i], ma_m[a], 0x7fffffff);
  }
  for(int i = 1; i <= n; i ++) {
    ma_d[i][i] = ++ tot;
    for(int j = i + 1; j <= n; j ++) {
      ma_d[i][j] = ++ tot;
#ifdef LOCAL
      printf("ma_d[%d][%d] : %d\n", i, j, tot);
#endif
    }
  }
  for(int i = 1; i <= n; i ++) {
    int d_i; scanf("%d", &d_i);
    if(d_i >= 0) {
      ans += d_i;
      ins_edge(0, ma_d[i][i], d_i);
    } else {
      ins_edge(ma_d[i][i], 1, -d_i);
    }
    ins_edge(ma_d[i][i], ma_p[i], 0x7fffffff);
    for(int j = i + 1; j <= n; j ++) {
      scanf("%d", &d_i);
      int np = ma_d[i][j];
      if(d_i >= 0) {
        ans += d_i;
        ins_edge(0, np, d_i);
        ins_edge(np, ma_d[i][j - 1], 0x7fffffff);
        ins_edge(np, ma_d[i + 1][j], 0x7fffffff);
      } else {
        ins_edge(np, 1, -d_i);
        ins_edge(np, ma_d[i][j - 1], 0x7fffffff);
        ins_edge(np, ma_d[i + 1][j], 0x7fffffff);
      }
    }
  }
  ans -= dinic();
#ifdef LOCAL
  for(int i = 0; i <= tot; i ++) {
    for(int j = first[i]; j; j = next[j]) {
      if(!(j & 1)) continue;
      int v = to[j];
      if(cap[j] == flow[j]) {
        printf("Edge (%d, %d, %d) deleted.\n", i, v, cap[j]);
      }
    }
  }
#endif
  printf("%d\n", ans);
  return 0;
}

[LibreOJ 2142][SHOI2017]相逢是问候

f**********ck我终于肝掉这道题了……在这种辣鸡题上浪费了大量时间

首先一堆\(c\)套起来的幂让我们想到了上帝与集合的正确用法那道题(这题题解我屯到现在没写,还是算了吧就写这题得了)……那让我们试试欧拉降幂公式?

观察欧拉降幂公式

\[a^x\equiv a^{x\mod \varphi(m) + \varphi(m)}\pmod{m}(a\geq m)\]

在多次降幂之后,最后肯定存在一个膜数为1,最后在对这一个地方进行操作的话,虽然\(A[i]\)的位置会被移到更高的次幂上,但是开始出现膜数1的地方一定只能取到1了,所以再对这些地方操作是不合理的。

每次取\(\varphi(x)\)都会使数至少下降一半,所以一个数最多被操作\(\log_2 p\)次。

然后我就跪了(虽然在BZOJ上水过去了),后来手肝gprof发现瓶颈在快速幂(逃

然后这TM怎么优化……膜了一发题解发现正解非常的KBS……

对于一个幂,他的指数肯定可以表示为\(2^{16} a + b\)且\(a\)最大的形式……然后我们可能用到的膜数是很少的,对于每种膜数大力预处理可能的\(b\)的答案和\(a\)的答案,然后求一次幂就\(O(1)\)了……

代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <cmath>
#include <algorithm>
#include <utility>
#include <climits>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/hash_policy.hpp>
const int maxn = 50005;
const int maxno = maxn << 2;
const int siz = 65536;
typedef long long ll;
inline ll phi(ll x) {
  ll bd = sqrt(x + 0.5);
  ll ans = x;
  for(ll i = 2; i <= bd; i ++) {
    if(x % i == 0LL) {
      ans = ans / i * (i - 1LL);
      while(x % i == 0LL) {
        x /= i;
      }
    }
  }
  if(x > 1LL) ans = ans / x * (x - 1LL);
  return ans;
}
int n, m;
ll p, c;
ll f[100]; int fn, bd;
ll pow1[70][siz], pow2[70][siz];
inline void process() {
  f[0] = p; fn = 1LL;
  for(int i = 1; ; i ++) {
    f[i] = phi(f[i - 1]);
    if(f[i] == 1LL) {
      fn = i;
      break;
    }
  }
  ll tmp = 1LL; bd = 0;
  if(c != 1LL) {
    while(tmp * c < p) {
      bd ++; tmp *= c;
    }
  } else {
    bd = 0x7fffffff;
  }
  f[69] = LLONG_MAX;
  for(int i = 0; i <= fn; i ++) {
    ll c1 = c % f[i]; ll c2 = c1;
    int t = 16; while(t --) c2 = (c2 * c2) % f[i];
    pow1[i][0] = 1LL; if(i == fn) pow1[i][0] = 0;
    for(int j = 1; j < siz; j ++) pow1[i][j] = (pow1[i][j - 1] * c1) % f[i];
    pow2[i][0] = 1LL; if(i == fn) pow2[i][0] = 0;
    for(int j = 1; j < siz; j ++) pow2[i][j] = (pow2[i][j - 1] * c2) % f[i];
  }
  int i = 69;
  for(int g = 0; g <= 0; g ++) {
    ll c1 = c % f[i]; ll c2 = c1;
    int t = 16; while(t --) c2 = (c2 * c2) % f[i];
    pow1[i][0] = 1LL;
    for(int j = 1; j < siz; j ++) pow1[i][j] = (pow1[i][j - 1] * c1) % f[i];
    pow2[i][0] = 1LL;
    for(int j = 1; j < siz; j ++) pow2[i][j] = (pow2[i][j - 1] * c2) % f[i];
  }
}
inline ll pow_mod(ll a, ll b, ll p) {
  return (pow1[p][b & (siz - 1)] * pow2[p][b >> 16]) % f[p];
}
ll A[maxn];
int be_done[maxn]; ll sumv[maxno];
bool break_down[maxno];
inline void maintain(int o) {
  int lc = o << 1, rc = o << 1 | 1;
  sumv[o] = (sumv[lc] + sumv[rc]);
  if(sumv[o] >= p) sumv[o] -= p;
  break_down[o] = (break_down[lc] && break_down[rc]);
}

inline ll get_ans(int i) {
  register ll up;
  for(int j = be_done[i]; j >= 0; j --) {
    if(j == be_done[i]) {
      up = A[i];
    } else {
      if(up > bd) {
        up = pow_mod(c % f[j], up, j) + f[j];
      } else {
        up = pow_mod(c, up, 69);
      }
    }
  }
  return up;
}
void build_tree(int o, int L, int R) {
  if(L == R) {
    be_done[L] = 0;
    break_down[o] = 0;
    sumv[o] = A[L];
  } else {
    int M = (L + R) / 2;
    build_tree(o << 1, L, M);
    build_tree(o << 1 | 1, M + 1, R);
    maintain(o);
  }
}
int ql, qr;
void modify(int o, int L, int R) {
  if(L == R) {
    be_done[L] ++;
    if(c == 1LL) {
      sumv[o] = 1LL;
      break_down[o] = true;
    } else {
      sumv[o] = get_ans(L) % p;
      if((A[L] != 0LL && be_done[L] == fn) || (A[L] == 0LL && be_done[L] == fn + 1)) break_down[o] = true;
    }
  } else {
    int M = (L + R) / 2;
    if(ql <= M && !break_down[o << 1]) {
      modify(o << 1, L, M);
    }
    if(qr > M && !break_down[o << 1 | 1]) {
      modify(o << 1 | 1, M + 1, R);
    }
    maintain(o);
  }
}
ll query(int o, int L, int R) {
  if(ql <= L && R <= qr) {
    return sumv[o];
  } else {
    int M = (L + R) / 2;
    ll ans = 0;
    if(ql <= M) {
      ans = (ans + query(o << 1, L, M));
      if(ans >= p) ans -= p;
    }
    if(qr > M) {
      ans = (ans + query(o << 1 | 1, M + 1, R));
      if(ans >= p) ans -= p;
    }
    return ans;
  }
}

int main() {
  // freopen("17.in", "r", stdin);
  // freopen("dummy", "w", stdout);
  scanf("%d%d%lld%lld", &n, &m, &p, &c);
  process();
  for(int i = 1; i <= n; i ++) scanf("%lld", &A[i]);
  build_tree(1, 1, n);
  while(m --) {
    int op; scanf("%d%d%d", &op, &ql, &qr);
    if(op == 0) {
      modify(1, 1, n);
    } else {
      printf("%lld\n", query(1, 1, n));
      // fflush(stdout);
    }
  }
  return 0;
}

[BZOJ 1857][SCOI2010]传送带

三分套三分入门题……

策略肯定是从A走到AB上一点,然后再走到CD上的一个点,再向D走。

很显然答案函数是一个关于那两个点下凸的东西(不会证?GeoGebra之类的东西画一下就好啦!还不如像我这样口胡),所以我们可以先对第一维三分,然后套上对第二维的三分……

代码:

/**************************************************************
	Problem: 1857
	User: danihao123
	Language: C++
	Result: Accepted
	Time:244 ms
	Memory:820 kb
****************************************************************/

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <cmath>
typedef double R;
const R eps = 1e-6;
int sign(R x) {
  if(fabs(x) < eps) {
    return 0;
  } else {
    if(x < 0) return -1;
    else return 1;
  }
}
struct Point {
  R x, y;
  Point(R qx = 0, R qy = 0) {
    x = qx; y = qy;
  }
};
typedef Point Vector;
Vector operator +(const Vector &a, const Vector &b) {
  return Vector(a.x + b.x, a.y + b.y);
}
Vector operator -(const Point &a, const Point &b) {
  return Vector(a.x - b.x, a.y - b.y);
}
Vector operator *(R x, const Vector &v) {
  return Point(v.x * x, v.y * x);
}
Vector operator *(const Vector &v, R x) {
  return Point(v.x * x, v.y * x);
}
R dot(const Vector &a, const Vector &b) {
  return a.x * b.x + a.y * b.y;
}
R times(const Vector &a, const Vector &b) {
  return a.x * b.y - a.y * b.x;
}
R dist(const Point &a, const Point &b) {
  return sqrt(dot(a - b, a - b));
}
bool cmp(const Point &a, const Point &b) {
  if(sign(a.x - b.x) == 0) {
    return a.y < b.y;
  } else {
    return a.x < b.x;
  }
}
Point A, B, C, D;
R p, q, r;
Vector D_AB, D_DC;
R f(const Point &AB, const Point &CD) {
  return (dist(AB, A) / p + dist(CD, D) / q + dist(AB, CD) / r);
}
R F(Point AB) {
  R L = 0, R = 1;
  int T = 200;
  while(T --) {
    double M1 = L + (R - L) / 3;
    double M2 = R - (R - L) / 3;
    Point P1 = D + M1 * D_DC;
    Point P2 = D + M2 * D_DC;
    double f1 = f(AB, P1), f2 = f(AB, P2);
    if(f1 < f2) {
      R = M2;
    } else {
      L = M1;
    }
  }
  return f(AB, D + L * D_DC);
}
R solve() {
  R L = 0, R = 1;
  int T = 200;
  while(T --) {
    double M1 = L + (R - L) / 3;
    double M2 = R - (R - L) / 3;
    Point P1 = A + M1 * D_AB;
    Point P2 = A + M2 * D_AB;
    double F1 = F(P1), F2 = F(P2);
    if(F1 < F2) {
      R = M2;
    } else {
      L = M1;
    }
  }
  return F(A + L * D_AB);
}

int main() {
  scanf("%lf%lf%lf%lf", &A.x, &A.y, &B.x, &B.y);
  scanf("%lf%lf%lf%lf", &C.x, &C.y, &D.x, &D.y);
  scanf("%lf%lf%lf", &p, &q, &r);
  D_AB = B - A; D_DC = C - D;
  printf("%.2lf\n", solve());
  return 0;
}

[BZOJ 3944]Sum

杜教筛模板题……(然而常数卡卡卡

我还是讲一讲基础的东西吧……

先考虑求欧拉函数的前缀和,我们先找一个它的卷积(这里用\(f=\varphi\ast 1\),且设\(g = 1\)),发现有(用\(S\)表示原函数的前缀和):

\[\sum_{i = 1}^n f(i) = \sum_{i = 1}^n g(i)\times S(\lfloor\frac{n}{i}\rfloor) = \frac{n(n + 1)}{2}\]

第二个和式的第一项是\(g(1)\times S(n) = S(n)\),也就是我们需要的答案。这样的话我们把第一项分割出来,后面的项大力求,然后用\(\frac{n(n + 1)}{2}\)减去那些项就得到了第一项。并且注意到\(\lfloor\frac{n}{i}\rfloor\)的取值种类是\(\sqrt{n}\)级别的,又可以大力优化一波。然后考虑预处理\(O(n^{\frac{2}{3}})\)级别的前缀和(时间和空间都允许),再优化一波。然后考虑把没预处理的状态扔一个哈希表里(反正看起来不会很多),又优化一波……

这样复杂度就是\(O(n^{\frac{2}{3}})\)的了……(又不会证了

求莫比乌斯函数的前缀和的方法类似(并且更简单一点)……

代码(又被卡成苟了):

/**************************************************************
	Problem: 3944
	User: danihao123
	Language: C++
	Result: Accepted
	Time:11292 ms
	Memory:127740 kb
****************************************************************/

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/hash_policy.hpp>
typedef long long ll;
const int N = 3500000;
const int maxn = N + 5;
ll phi[maxn], phi_S[maxn];
ll mu[maxn], mu_S[maxn];
void sieve() {
  static bool vis[maxn];
  static int prm[maxn]; int cnt = 0;
  phi[1] = 1; mu[1] = 1;
  for(int i = 2; i <= N; i ++) {
    if(!vis[i]) {
      prm[cnt ++] = i;
      mu[i] = -1; phi[i] = i - 1;
    }
    for(int j = 0; j < cnt && prm[j] * (ll(i)) <= (ll(N)); j ++) {
      int v = prm[j] * i;
      vis[v] = true;
      if(i % prm[j] == 0) {
        mu[v] = 0; phi[v] = phi[i] * prm[j];
        break;
      } else {
        mu[v] = mu[i] * -1;
        phi[v] = phi[i] * phi[prm[j]];
      }
    }
  }
  for(int i = 1; i <= N; i ++) {
    phi_S[i] = phi_S[i - 1] + phi[i];
    mu_S[i] = mu_S[i - 1] + mu[i];
  }
}
__gnu_pbds::gp_hash_table<ll, ll> mu_d;
ll mu_sum(ll n) {
  if(n <= (ll(N))) return mu_S[n];
  if(mu_d.find(n) != mu_d.end()) {
    return mu_d[n];
  }
  ll ans = 1;
  for(ll i = 2; i <= n;) {
    ll nx = n / (n / i);
    ans -= (nx - i + 1LL) * mu_sum(n / i);
    i = nx + 1LL;
  }
  mu_d[n] = ans;
  return ans;
}
ll pre_sum(ll n) {
  return (n * (n + 1LL) / 2LL);
}
__gnu_pbds::gp_hash_table<ll, ll> phi_d;
ll phi_sum(ll n) {
  if(n <= (ll(N))) return phi_S[n];
  if(phi_d.find(n) != phi_d.end()) {
    return phi_d[n];
  }
  ll ans = pre_sum(n);
  for(ll i = 2; i <= n;) {
    ll nx = n / (n / i);
    ans -= (nx - i + 1LL) * phi_sum(n / i);
    i = nx + 1LL;
  }
  phi_d[n] = ans;
  return ans;
}

int main() {
  sieve();
  int T; scanf("%d", &T);
  while(T --) {
    int n; scanf("%d", &n);
    printf("%lld %lld\n", phi_sum(n), mu_sum(n));
  }
  return 0;
}

[BZOJ 4816][SDOI2017]数字表格

当年在考场上一点反演都不会啊啊啊啊啊啊

这题虽然牵扯到了对积的反演,但其实也不是很难。先让我们大力推导一波(逃

考虑枚举柿子中的\(f[(i, j)]\)中的最大公约数(设为\(k\)),然后式子变成了(这里不妨偷个懒,钦定\(n\leq m\)):

\[\prod_{k = 1}^n f[k]^{\sum_{i = 1}^n \sum_{j = 1}^m\:[(i, j) = k]}\]

然后发现上面那个指数就是Problem b的式子,显然可化为\(\sum_{d = 1}^n \mu(d) \lfloor\frac{n}{kd}\rfloor \lfloor\frac{m}{kd}\rfloor\)。然后似乎化简没有余地了,但是根据直觉,我们可以钦定\(T = kd\),然后枚举\(T\)。

然后柿子变成了:

\[\prod_{T = 1}^n (\prod_{k\mid T} f[k]^{\mu(\tfrac{T}{k})})^{\lfloor\tfrac{n}{T}\rfloor \lfloor\tfrac{m}{T}\rfloor}\]

柿子中的\(\prod_{k\mid T} f[k]^{\mu(\tfrac{T}{k})}\)是一个类似狄利克雷卷积的东西(取完对数就是了),可以枚举约数预处理。并且这个东西很不错的一点就是上面的指数是莫比乌斯函数,可以取到的值只有那三种,不需要快速幂。

然后剩下的部分就和通常的反演题一般无二了……

代码(卡线过的蒟蒻……求轻喷……但我在LOJ上跑的还挺快的):

/**************************************************************
    Problem: 4816
    User: danihao123
    Language: C++
    Result: Accepted
    Time:50840 ms
    Memory:33048 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <utility>
#include <cctype>
#include <cassert>
const int N = 1000000;
const int maxn = N + 5;
typedef long long ll;
const ll ha = 1000000007LL;
bool p_flag = false;
ll pow_mod(ll a, ll b) {
  if(!b) return 1LL;
  ll ans = pow_mod(a, b / 2LL);
  ans = (ans * ans) % ha;
  if(1LL & b) ans = (ans * a) % ha;
#ifdef LOCAL
  if(p_flag) printf("%lld^%lld : %lld\n", a, b, ans);
  if(b == ha - 2LL) p_flag = false;
#endif
  return ans;
}
ll inv(ll x) {
  // if(x == 1LL) p_flag = true;
  return pow_mod(x, ha - 2LL);
}
int miu[maxn];
void sieve() {
  static int prm[maxn]; int cnt = 0;
  static bool vis[maxn];
  miu[1] = 1;
  for(int i = 2; i <= N; i ++) {
    if(!vis[i]) {
      miu[i] = -1;
      prm[cnt ++] = i;
    }
    for(int j = 0; j < cnt && prm[j] * i <= N; j ++) {
      int v = prm[j] * i;
      vis[v] = true;
      if(i % prm[j] == 0) {
        miu[v] = 0;
        break;
      } else {
        miu[v] = miu[i] * -1;
      }
    }
  }
}
ll f[maxn], inv_d[maxn], F[maxn];
void process() {
  sieve();
  f[1] = 1LL; inv_d[1] = 1LL;
  for(int i = 2; i <= N; i ++) {
    f[i] = (f[i - 1] + f[i - 2]) % ha;
    inv_d[i] = inv(f[i]);
    assert(f[i] != 0LL); assert(inv_d[i] != 0LL);
  }
  for(int i = 1; i <= N; i ++) F[i] = 1LL;
  for(int i = 1; i <= N; i ++) {
    for(int j = i; j <= N; j += i) {
      int res = j / i;
      if(miu[res] == 1) {
        F[j] = (F[j] * f[i]) % ha;
#ifdef LOCAL
        if(j <= 3) printf("f[%d](%lld) -> F[%d]\n", i, f[i], j);
#endif
        continue;
      }
      if(miu[res] == -1) {
        F[j] = (F[j] * inv_d[i]) % ha;
#ifdef LOCAL
        if(j <= 3) printf("inv_f[%d](%lld) -> F[%d]\n", i, inv_d[i], j);
#endif
      }
    }
  }
  F[0] = 1LL;
  bool flag = true;
  for(int i = 1; i <= N; i ++) {
    F[i] = (F[i] * F[i - 1]) % ha;
    if(flag && F[i] == 0LL) {
      printf("SF[%d] has been 0.\n", i);
      flag = false;
    }
  }
}
 
ll calc(ll n, ll m) {
  if(n > m) std::swap(n, m);
  ll ans = 1LL;
  for(ll i = 1; i <= n;) {
    ll nx = std::min(n / (n / i), m / (m / i));
#ifdef LOCAL
    // printf("nx of %lld : %lld\n", i, nx);
#endif
    ll mulv = pow_mod((F[nx] * inv(F[i - 1])) % ha, (n / i) * (m / i));
    ans = (ans * mulv) % ha;
    i = nx + 1LL;
  }
  return ans;
}
 
int main() {
  process();
  int T; scanf("%d", &T);
  while(T --) {
    int n, m; scanf("%d%d", &n, &m);
    printf("%lld\n", calc(n, m));
  }
  return 0;
}

[BZOJ 5059]前鬼后鬼的守护

这不是可并堆裸题吗……我这种傻逼只会用线段树合并做

显然一个集合的中位数可以用平衡树或者权值线段树求出,然后为了方便合并我用了权值线段树。

然后考虑对于一个递减的序列,选择中位数一定更优,一个递增的序列当然xjb选即可。我们可以用一个栈来维护当前的所有答案区间,然后考虑一个新的点加进来。这个点的答案如果比栈首答案小,那么就违反了单调不降的要求,我们需要把两段合并。然后完了……

代码:

/**************************************************************
	Problem: 5059
	User: danihao123
	Language: C++
	Result: Accepted
	Time:2348 ms
	Memory:250520 kb
****************************************************************/

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <cassert>
#include <algorithm>
#include <utility>
#include <stack>
const int BUFSIZE = 20 * 1024 * 1024;
int n;
struct Node {
  Node *lc, *rc;
  int sumv;
};
Node pool[BUFSIZE];
Node *nil, *cur;
void init_pool() {
  cur = nil = pool;
  nil -> lc = nil -> rc = nil;
  nil -> sumv = 0;
}
Node *alloc_node(int v = 0, Node *lc = nil, Node *rc = nil) {
  cur ++;
  cur -> lc = lc; cur -> rc = rc;
  cur -> sumv = v;
  return cur;
}
void modify(Node* &o, int L, int R, int p, int v) {
  if(o == nil) {
    o = alloc_node(0);
  }
  o -> sumv += v;
  if(L < R) {
    int M = (L + R) / 2;
    if(p <= M) {
      modify(o -> lc, L, M, p, v);
    } else {
      modify(o -> rc, M + 1, R, p, v);
    }
  }
}
Node *merge(Node* &tag, Node* &src) {
  if(src == nil) return tag;
  if(tag == nil) {
    return src;
  }
  tag -> sumv += src -> sumv;
  tag -> lc = merge(tag -> lc, src -> lc);
  tag -> rc = merge(tag -> rc, src -> rc);
  return tag;
}
int kth(Node *o, int L, int R, int k) {
  if(L == R) {
    return L;
  } else {
    int M = (L + R) / 2;
    if(k <= o -> lc -> sumv) {
      return kth(o -> lc, L, M, k);
    } else {
      return kth(o -> rc, M + 1, R, k - o -> lc -> sumv);
    }
  }
}
void gou_tree(Node *o, int L, int R) {
  printf("Tree (%d, %d, %d)\n", L, R, o -> sumv);
  int M = (L + R) / 2;
  if(L < R) {
    gou_tree(o -> lc, L, M);
    gou_tree(o -> rc, M + 1, R);
  }
}
int lsiz;
const int maxn = 500005;
int A[maxn], A2[maxn];
void break_up() {
  std::sort(A2 + 1, A2 + 1 + n);
  lsiz = std::unique(A2 + 1, A2 + 1 + n) - A2 - 1;
  for(int i = 1; i <= n; i ++) {
    A[i] = std::lower_bound(A2 + 1, A2 + 1 + lsiz, A[i]) - A2;
  }
}
struct Interval {
  Node *rt;
  int l, r; int siz, ans;
  Interval() {
    rt = nil;
    siz = ans = 0;
  }
  Interval(int p) {
    rt = nil; modify(rt, 1, lsiz, p, 1);
    siz = 1; ans = p;
  }
};
typedef long long ll;
ll solve() {
  break_up(); init_pool();
  std::stack<Interval> S;
  for(int i = 1; i <= n; i ++) {
    Interval nv(A[i]);
    nv.l = nv.r = i;
    while(!S.empty() && S.top().ans > nv.ans) {
      Interval g = S.top(); S.pop();
#ifdef LOCAL
    printf("ans of [%d, %d] : %d\n", g.l, g.r, A2[g.ans]);
#endif
      nv.l = g.l; nv.siz += g.siz;
      nv.rt = merge(nv.rt, g.rt);
#ifdef LOCAL
      gou_tree(nv.rt, 1, lsiz);
#endif
      nv.ans = kth(nv.rt, 1, lsiz, (nv.siz + 1) / 2);
    }
    S.push(nv);
  }
  ll ans = 0;
  while(!S.empty()) {
    Interval g = S.top(); S.pop();
#ifdef LOCAL
    printf("ans of [%d, %d] : %d\n", g.l, g.r, A2[g.ans]);
#endif
    for(int i = g.l; i <= g.r; i ++) {
      ans += abs(A2[A[i]] - A2[g.ans]);
    }
  }
  return ans;
}

int main() {
  scanf("%d", &n);
  for(int i = 1; i <= n; i ++) {
    scanf("%d", &A[i]);
    A2[i] = A[i];
  }
  printf("%lld\n", solve());
  return 0;
}

[BZOJ 3527][ZJOI2014]力

现在才明白卷积真的是the deep, dark fantasies……(逃

首先约去\(q_i\),得到:

\[E_j = \sum_{i < j}\frac{q_i}{(j - i)^2} - \sum_{j < i}\frac{q_i}{(j - i)^2}\]

注意到如果很容易得到等式前半部分的高效求法,后半部分把序列反过来就能做了。

那么我们会发现,设\(g_x = \frac{1}{x^2}\),然后式子前半部分(姑且称作\(p_j\))可表示为:

\[p_j = \sum_{i < j} q_i g_{j - i}\]

这不就是个卷积吗?FFT一波即可。

代码:

/**************************************************************
    Problem: 3527
    User: danihao123
    Language: C++
    Result: Accepted
    Time:9984 ms
    Memory:28952 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cctype>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <utility>
#include <queue>
#include <cassert>
#include <complex>
typedef long double R;
typedef std::complex<R> C;
const R eps = 1e-3;
int sign(R x) {
  if(fabs(x) < eps) {
    return 0;
  } else {
    if(x < 0) {
      return -1;
    } else {
      return 1;
    }
  }
}
const int maxn = 400005;
const double pi = acos(-1);
int bi, len;
int flip(int x) {
  int ans = 0;
  for(int i = 0; i < bi; i ++) {
    if(x & (1 << i)) {
      ans += 1 << (bi - i - 1);
    }
  }
  return ans;
}
void fft(C *A, double g = 1) {
  for(int i = 0; i < len; i ++) {
    int R = flip(i);
#ifdef LOCAL
    // printf("The flipping of %d is %d\n", i, R);
#endif
    if(i < R) {
      std::swap(A[R], A[i]);
    }
  }
  for(int L = 1; L < len; L <<= 1) {
    C xi_n(cos(pi / (double(L))), sin(g * pi / (double(L))));
    for(int i = 0; i < len; i += L << 1) {
      C xi(1, 0);
      for(int j = i; j < i + L; j ++) {
        C a = A[j], b = A[j + L];
        A[j] = a + xi * b;
        A[j + L] = a - xi * b;
        xi = xi * xi_n;
      }
    }
  }
}
 
int main() {
  int n;
  static C A[maxn], rd[maxn]; static R ans[maxn];
  static R q[maxn];
  scanf("%d", &n);
  bi = 0; len = 1;
  while(len <= 2 * n) {
    len <<= 1;
    bi ++;
  }
  for(int i = 0; i < n; i ++) {
    scanf("%Lf", &q[i]); A[i] = q[i];
  }
  assert(sign(rd[0].real()) == 0);
  for(int i = 1; i < n; i ++) {
    rd[i] = 1.00 / ((R(i)) * (R(i)));
#ifdef LOCAL
    printf("rd[%d] : %.5Lf\n", i, rd[i].real());
#endif
  }
  fft(A); fft(rd);
  for(int i = 0; i < len; i ++) A[i] *= rd[i];
  fft(A, -1);
  for(int i = 0; i < n; i ++) {
    ans[i] += A[i].real() / len;
#ifdef LOCAL
    printf("delta_v of %d : %.5Lf\n", i, A[i].real() / len);
#endif
  }
  std::reverse(q, q + n);
  for(int i = 0; i < len; i ++) A[i] = 0;
  for(int i = 0; i < n; i ++) {
    A[i] = q[i];
  }
  fft(A);
  for(int i = 0; i < len; i ++) A[i] *= rd[i];
  fft(A, -1);
  for(int i = 0; i < n; i ++) {
    ans[i] -= A[n - 1 - i].real() / len;
    printf("%.3Lf\n", ans[i]);
  }
  return 0;
}