Kyoya and Train

题目描述

Kyoya Ootori wants to take the train to get to school. There are $n$ train stations and $m$ one-way train lines going between various stations. Kyoya is currently at train station $1$, and the school is at station $n$. To take a train, he must pay for a ticket, and the train also takes a certain amount of time. However, the trains are not perfect and take random amounts of time to arrive at their destination. If Kyoya arrives at school strictly after $t$ time units, he will have to pay a fine of $x$.
Each train line is described by a ticket price, and a probability distribution on the time the train takes. More formally, train line $i$ has ticket cost $c_i$, and a probability distribution $p_{i, k}$ which denotes the probability that this train will take $k$ time units for all $1 \le k \le t$. Amounts of time that each of the trains used by Kyouya takes are mutually independent random values (moreover, if Kyoya travels along the same train more than once, it is possible for the train to take different amounts of time and those amounts are also independent one from another).
Kyoya wants to get to school by spending the least amount of money in expectation (for the ticket price plus possible fine for being late). Of course, Kyoya has an optimal plan for how to get to school, and every time he arrives at a train station, he may recalculate his plan based on how much time he has remaining. What is the expected cost that Kyoya will pay to get to school if he moves optimally?

题意概述

有$n$座城市和$m$条铁路,第$i$条铁路从$a_i$号城市出发,到达$b_i$号城市,票价为$c_i$。某人要在$t$个单位时间内从$1$号城市到$n$号城市,如果超过时间就会被罚款$x$。已知每次搭乘第$i$条铁路有$p_{i, k}$的概率需要$k \; (1 \le k \le t)$个单位时间。求在采取最优策略的情况下总花费的期望。
数据范围:$2 \le n \le 50, \; 1 \le m \le 100, \; 1 \le t \le 20000, \; 0 \le x \le 10^6$。

算法分析

首先考虑最暴力的DP。令$f_{i, j}$表示当前时刻为$j$,在最优策略下从$i$号城市到$n$号城市的总花费的期望。那么
$$
f_{i, j}=
\begin{cases}
\min(c_e+\sum_{k=1}^t p_{e, k} \cdot f_{b_e, j+k} \mid e \in [1, m] \land a_e=i), & i \lt n \land j \le t \\
d_{i, n}+x , & i \lt n \land j \gt t \\
0, & i=n \land j \le t \\
x, & i=n \land j \gt t
\end{cases}
$$
第二部分中的$d_{i, n}$表示在只考虑铁路票价的情况下从$i$号城市到$n$号城市的最小花费(因为已经超时了,所以不如走总票价最少的路)。
这样做的复杂度是$O(mt^2)$的,需要对它进行优化。令
$$S_{e, j}=\sum_{k=1}^t p_{e, k} \cdot f_{b_e, j+k}$$
转移方程的第一部分变为
$$f_{i, j}=\min(c_e+S_{e, j} \mid e \in [1, m] \land a_e=i)$$
对于$S$,它类似于卷积,可以将其中一部分翻转后用FFT求值;对于$f$,可以用CDQ分治,统计较大的$j$对较小的$j$的贡献。
具体来说,假设我们要计算$l \le j \le r$的$f_{i, j}$,令$mid=\lfloor {l+r \over 2} \rfloor$,那么先计算$mid \lt j \le r$的$f_{i, j}$,用这些来更新$l \le j \le mid$的$S_{e, j}$($m$条边分别用FFT更新,复杂度$O(m(r-l)\log (r-l))$),最后计算$l \le j \le mid$的$f_{i, j}$。当$l=r$时$f$可以直接由$S$得到。总复杂度降至$O(mt\log^2t)$。
实现时细节很多,要注意DP边界条件和FFT下标变换。

代码

/*
 * Your nature demands love and your happiness depends on it.
 */

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>

static int const N = 55;
static int const M = 105;
static int const T = 100000;
static double const PI = acos(-1);
int rev[T];

class Point {
public:
  double x, y;
  Point(double _x = 0, double _y = 0) : x(_x), y(_y) {}
  std::pair<double, double> pair() {
    return std::make_pair(x, y);
  }
  friend Point operator+(Point const &a, Point const &b) {
    return Point(a.x + b.x, a.y + b.y);
  }
  friend Point operator-(Point const &a, Point const &b) {
    return Point(a.x - b.x, a.y - b.y);
  }
  friend Point operator*(Point const &a, Point const &b) {
    return Point(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
  }
  Point operator/(double const &n) {
    return Point(x / n, y / n);
  }
} wn[T], A[T], B[T];

int init(int n) {
  int m = n, l = 0;
  for (n = 1; n <= m; n <<= 1, ++l)
    ;
  for (int i = 1; i < n; ++i)
    rev[i] = rev[i >> 1] >> 1 | (i & 1) << l - 1;
  for (int i = 0; i < n >> 1; ++i)
    wn[i] = Point(cos(2 * PI / n * i), sin(2 * PI / n * i));
  return n;
}

void fft(Point *a, int n, int inv = 0) {
  for (int i = 0; i < n; ++i)
    if (i < rev[i])
      std::swap(a[i], a[rev[i]]);
  for (int i = 1; i < n; i <<= 1)
    for (int j = 0; j < n; j += i << 1)
      for (int k = 0; k < i; ++k) {
        Point x = a[j + k], y = wn[n / (i << 1) * k] * a[j + k + i];
        a[j + k] = x + y, a[j + k + i] = x - y;
      }
  if (inv) {
    std::reverse(a + 1, a + n);
    for (int i = 0; i < n; ++i)
      a[i] = a[i] / n;
  }
}

int n, m, t, x, mp[N][N];
double p[M][T], s[M][T], f[N][T];
struct Line {
  int a, b, c;
} li[M];

void update(int l, int r) {
  int mid = l + r >> 1, len = init(r - l + r - mid - 2);
  for (int i = 1; i <= m; ++i) {
    for (int j = 0; j < len; ++j)
      A[j] = B[j] = 0;
    for (int j = mid + 1; j <= r; ++j)
      A[j - mid - 1] = f[li[i].b][r - j + mid + 1];
    for (int j = 1; j <= r - l; ++j)
      B[j - 1] = p[i][j];
    fft(A, len), fft(B, len);
    for (int j = 0; j < len; ++j)
      A[j] = A[j] * B[j];
    fft(A, len, 1);
    for (int j = l; j <= mid; ++j)
      s[i][j] += A[r - j - 1].x;
  }
}

void solve(int l, int r) {
  if (l == r) {
    for (int i = 1; i <= m; ++i)
      f[li[i].a][l] = std::min(f[li[i].a][l], s[i][l] + li[i].c);
    return;
  }
  int mid = l + r >> 1;
  solve(mid + 1, r), update(l, r), solve(l, mid);
}

int main() {
  scanf("%d%d%d%d", &n, &m, &t, &x);
  memset(mp, 0x3f, sizeof mp);
  for (int i = 1; i <= m; ++i) {
    scanf("%d%d%d", &li[i].a, &li[i].b, &li[i].c);
    mp[li[i].a][li[i].b] = std::min(mp[li[i].a][li[i].b], li[i].c);
    for (int j = 1; j <= t; ++j)
      scanf("%lf", &p[i][j]), p[i][j] /= 100000;
  }
  for (int i = 1; i <= n; ++i)
    for (int j = 1; j <= n; ++j)
      for (int k = 1; k <= n; ++k)
        mp[j][k] = std::min(mp[j][k], mp[j][i] + mp[i][k]);
  for (int i = 1; i <= n; ++i)
    for (int j = 0; j <= t << 1; ++j)
      if (i == n)
        if (j <= t)
          f[i][j] = 0;
        else
          f[i][j] = x;
      else
        if (j <= t)
          f[i][j] = 1e9;
        else
          f[i][j] = x + mp[i][n];
  update(0, t << 1), solve(0, t), printf("%.8lf\n", f[1][0]);
  return 0;
}

Optimized Input and Output in C++

Background

C++标准库<iostream>中的cin, cout相比<cstdio>中的scanf, printf要方便许多,但速度也较为缓慢。在输入输出的数据类型不是很特殊时,可以使用以下代码来代替cin, cout(要包含<cctype>, <cstdio>,不要包含<iostream>)。

Code

class InputStream {
private:
  static int const BUFFER_SIZE = 10000000;
  char *buffer;

public:
  InputStream() {
    buffer = new char[BUFFER_SIZE]();
    fread(buffer, sizeof(char), BUFFER_SIZE, stdin);
  }
  InputStream &operator>>(int &n) {
    for (; *buffer < '0' || *buffer > '9'; ++buffer)
      ;
    for (n = 0; *buffer >= '0' && *buffer <= '9';
         (n *= 10) += *(buffer++) - '0')
      ;
    return *this;
  }
  InputStream &operator>>(char &c) {
    for (; !isprint(*buffer) || *buffer == ' '; ++buffer)
      ;
    c = *(buffer++);
    return *this;
  }
  InputStream &operator>>(char *s) {
    char *t = s;
    for (; !isprint(*buffer) || *buffer == ' '; ++buffer)
      ;
    for (; isprint(*buffer) && *buffer != ' '; *(t++) = *(buffer++))
      ;
    *t = '\0';
    return *this;
  }
} cin;

class OutputStream {
private:
  static int const BUFFER_SIZE = 10000000;
  char *buffer, *cur;

public:
  OutputStream() { buffer = cur = new char[BUFFER_SIZE](); }
  ~OutputStream() { fwrite(buffer, sizeof(char), cur - buffer, stdout); }
  OutputStream &operator<<(int const &n) {
    for (sprintf(cur, "%d", n); *cur; ++cur)
      ;
    return *this;
  }
  OutputStream &operator<<(char const &c) {
    *(cur++) = c;
    return *this;
  }
  OutputStream &operator<<(char *s) {
    for (char *t = s; *t; *(cur++) = *(t++))
      ;
    return *this;
  }
} cout;

static char const endl = '\n';

Introduction to KMP

Background

在字符串理论中,经常要求一个字符串$T$在另一个字符串$S$中所有出现的位置。这个过程就是字符串匹配。我们称$S$为文本,$T$为模式。下面介绍一种时间复杂度为线性的算法KMP,可以解决此类问题。

Introduction

令$S_i$表示字符串$S$的第$i$个字符,$S_{i, j}$表示字符串$S$的第$i$到第$j$个字符构成的字符串。为了方便讲解,字符串下标从$1$开始。
在KMP中有一个重要的数组,我们称它为$nxt$。对于一个字符串$S$,$nxt_i=\max(x \mid x \lt i \land S_{1, x}=S_{i-x+1, i})$。特别的,令$nxt_1=0$。
考虑如何求出模式$T$的$nxt$数组。最暴力的暴力时间复杂度为$O(n^2)$,但这实际上浪费了很多有用的信息。其实可以从$nxt_{i-1}$推到$nxt_i$。
我们来形象化这个过程

    v
acbacabacbacb
           ^

假设我们已经求出$nxt_{12}=5, \; nxt_5=2$,当前两个指针分别指向$S_{12}$和$S_{nxt_{12}}=S_5$,求$nxt_{13}$。在这里,两个指向$S_p, S_q$的指针表示$S_{p-q+1, p}=S_{1, q}$。
比较两个指针的后一位。易知,若$S_{13}=S_6$,则$nxt_{13}=nxt_{12}+1=6$,可以将两个指针都向后移一位。
否则,我们已经知道$nxt_{12}=5$,即$S_{1, 5}=S_{8, 12}$,而$nxt_5=2$,即$S_{1, 2}=S_{4, 5}=S_{11, 12}$。这意味着我们可以把指向$S_{nxt_{12}}$的指针重新指向$S_{nxt_{nxt_{12}}}=S_2$,如下

 v
acbacbaacbacb
           ^

此时仍然满足$S_{12-2+1, 12}=S_{1, 2}$。再次比较两个指针的后一位,发现$S_{13}=S_3$,于是$nxt_{13}=nxt_{2}+1=3$,将两个指针向后移一位。当然,若$S_{13} \neq S_3$则需要将指向$S_{nxt_5}$的指针再指向$S_{nxt_{nxt_5}}$,重复上述过程,直到指针指向$S_0$,那么$nxt_{13}=0$。
代码并不复杂(该代码中的字符串下标从$0$开始)

void kmp(char c[], int len) {
  nxt[0] = -1;
  for (int i = 1; i < len; ++i) {
    for (nxt[i] = nxt[i - 1]; ~nxt[i] && c[nxt[i] + 1] != c[i];
         nxt[i] = nxt[nxt[i]])
      ;
    if (c[nxt[i] + 1] == c[i])
      ++nxt[i];
  }
}

在得到模式$T$的$nxt$数组后,我们就可以用类似的方法求出$T$在文本$S$中所有出现的位置了。思路也是维护两个指针,一个在$S$上,指向$S_p$,另一个在$T$上,指向$T_q$,同样要满足$S_{p-q+1, p}=T_{1, q}$。显然,当$q=|T|$时,就找到了一个出现的位置。

Application

KMP算法除了可以用于字符串匹配,还能求出一个循环串的最小循环节。循环串是由至少两个相同的字符串拼接起来得到的字符串。
具体来讲,设$i=|S|$,$nxt_i \gt 0 \land i-nxt_i \mid i \Rightarrow S_{1, i}$是循环串,它的最小循环节是$S_{1, i-nxt_i}$。
要证明以上结论,只需要证明:

  • 引理1 所有非循环串均不满足上述条件。
    证明 显然,在满足上述条件的情况下,$S$的确是循环串,$S_{1, i-nxt_i}$就是$S$的一个循环节。

  • 引理2 对于两个字符串$A, B$,若$AB=BA$,则$AB$是循环串,$A_{1, (|A|, |B|)}=B_{1, (|A|, |B|)}$是一个循环节。其中$(a, b)$表示$a, b$的最大公约数。
    证明 设$|A| \le |B|$,由$AB=BA$可知,$A=B_{1, |A|}=B_{|A|+1, 2|A|}=\cdots$,那么当$(|A|) \mid (|B|)$时,$A$是$B$的一个循环节(或$A=B$),因此$AB$是循环串。
    当$(|A|) \nmid (|B|)$时,令$k=\left\lfloor {|B| \over |A|} \right\rfloor$。
    考虑字符串$S=AB, \; T=S_{(k-1)|A|+1, |S|}$,可知$T_{1, |A|}=T_{|T|-|A|+1, |T|}=A, \; T_{1, |T|-|A|}=T_{|A|+1, |T|}=B_{(k-1)|A|+1, |B|}$。令$A’=B_{(k-1)|A|+1, |B|}, \; B’=A$,那么我们又得到了$A’B’=B’A’$的形式,此时$|A’|=|B| \bmod |A|, \; |B’|=|A|$。因此,最后一定能得到长为$(|A|, |B|)$的字符串,它是$A, B$的一个循环节,也是$AB$的一个循环节。

  • 引理3 若字符串$S$是循环串,则$S$的最小循环节是$S_{1, i-nxt_i}$。
    证明 令$x$为最小循环节的长度,只要证明不存在$y$满足$y \lt x$且$S_{1, i-y}=S_{y+1, i}$。
    假设存在这样的$y$,那么$S_{y+1, y+x}=S_{1, x}=S_{x+1, 2x}$。令$A=S_{y+1, x}, \; B=S_{x+1, y+x}$,可以发现$A$是$S_{1, x}$的后缀、$S_{y+1, y+x}$的前缀,$B$是$S_{x+1, 2x}$的前缀、$S_{y+1, y+x}$的后缀,即$S_{1, x}=AB=BA$,根据引理2,$S_{1, x}$是循环串,这与$x$是最小循环节的长度矛盾。因此这样的$y$不存在,而因为$S_{1, x}$是一个循环节,所以$S_{1, i-x}=S_{x+1, i}$,根据KMP的原理,$nxt_i=i-x$,即$S$的最小循环节是$S_{1, i-nxt_i}$。

由此命题得证。

Picnic Cows

题目描述

It’s summer vocation now. After tedious milking, cows are tired and wish to take a holiday. So Farmer Carolina considers having a picnic beside the river. But there is a problem, not all the cows consider it’s a good idea! Some cows like to swim in West Lake, some prefer to have a dinner in Shangri-la, and others want to do something different. But in order to manage expediently, Carolina coerces all cows to have a picnic!
Farmer Carolina takes her $N$ cows to the destination, but she finds every cow’s degree of interest in this activity is so different that they all loss their interests. So she has to group them to different teams to make sure that every cow can go to a satisfied team. Considering about the security, she demands that there must be no less than $T$ cows in every team. As every cow has its own interest degree of this picnic, we measure this interest degree’s unit as “Moo~“. Cows in the same team should reduce their Moo~ to the one who has the lowest Moo~ in this team – It’s not a democratical action! So Carolina wishes to minimize the TOTAL reduced Moo~s and groups $N$ cows into several teams.
For example, Carolina has $7$ cows to picnic and their Moo~ are ‘$8 \; 5 \; 6 \; 2 \; 1 \; 7 \; 6$’ and at least $3$ cows in every team. So the best solution is that cow No. $2, 4, 5$ in a team (reduce $((2-1)+(5-1))$ Moo~) and cow No. $1, 3, 6, 7$ in a team (reduce $((7-6)+(8-6))$ Moo~), the answer is $8$.

题意概述

有$N$头牛,每头牛有一个快乐值。要将牛分成若干组,使得每组至少有$T$头牛。分在一组中的牛的快乐值都会变成那组中最小的快乐值。求分组前后总快乐值之差的最小值。
数据范围:$1 \lt T \le N \le 4 \times 10^5$。

算法分析

将牛的快乐值$a_i$从小到大排序。根据贪心策略,一定存在最优分组方案是将排序后的牛分成若干连续段。
令$f_i$表示排序后前$i$头牛分组前后总快乐值之差的最小值,$s_i=\sum_{j=1}^i a_j$,$b_i=a_{i+1}$。可以得到转移方程
$$f_i=\min(f_j+(s_i-s_j)-(i-j)b_j \mid i-j \ge T)$$
假设$k \lt j \lt i$,且决策$j$比决策$k$更优,那么
$$
\begin{align}
f_j+(s_i-s_j)-(i-j)b_j &\le f_k+(s_i-s_k)-(i-k)b_k \\
(f_j-s_j+j \cdot b_j)-(f_k-s_k+k \cdot b_k) &\le i(b_j-b_k)
\end{align}
$$
$i$单调递增,因此可以斜率优化。但要注意两点:当$i-j \lt T$时$f_j$不能转移到$f_i$;当$i \lt T$时$f_i$没有意义。

代码

/*
 * You attempt things that you do not even plan because of your extreme
 * stupidity.
 */

#include <algorithm>
#include <cstdio>
#include <cstring>

#define int long long

static int const N = 400005;
int a[N], f[N], s[N], b[N], que[N];

int dety(int i, int j) {
  return f[j] - s[j] + b[j] * j - (f[i] - s[i] + b[i] * i);
}

int detx(int i, int j) { return b[j] - b[i]; }

signed main() {
  for (int n, k; ~scanf("%lld%lld", &n, &k);) {
    for (int i = 1; i <= n; ++i)
      scanf("%lld", &a[i]);
    std::sort(a + 1, a + n + 1);
    for (int i = 1; i <= n; ++i)
      s[i] = s[i - 1] + a[i], b[i - 1] = a[i];
    int qb = 0, qe = 1;
    que[0] = 0;
    for (int i = 1; i <= n; ++i) {
      for (; qb + 1 < qe &&
             dety(que[qb], que[qb + 1]) <= i * detx(que[qb], que[qb + 1]);
           ++qb)
        ;
      f[i] = f[que[qb]] + s[i] - s[que[qb]] - b[que[qb]] * (i - que[qb]);
      if (i - k + 1 >= k) {
        for (;
             qb + 1 < qe &&
             dety(que[qe - 2], que[qe - 1]) * detx(que[qe - 1], i - k + 1) >=
                 dety(que[qe - 1], i - k + 1) * detx(que[qe - 2], que[qe - 1]);
             --qe)
          ;
        que[qe++] = i - k + 1;
      }
    }
    printf("%lld\n", f[n]);
  }
  return 0;
}

Print Article

题目描述

Zero has an old printer that doesn’t work well sometimes. As it is antique, he still like to use it to print articles. But it is too old to work for a long time and it will certainly wear and tear, so Zero use a cost to evaluate this degree.
One day Zero want to print an article which has $N$ words, and each word $i$ has a cost $C_i$ to be printed. Also, Zero know that print $k$ words in one line will cost
$$\left(\sum_{i=1}^k C_i\right)^2+M$$
$M$ is a const number.
Now Zero want to know the minimum cost in order to arrange the article perfectly.

题意概述

一篇文章中有$N$个单词,输入第$i$个单词的代价为$C_i$。将连续$k$个单词排在一行的总代价为$\left(\sum_{i=1}^k C_i\right)^2+M$。求输入这篇文章的最小代价。
数据范围:$0 \le N \le 5 \times 10^5, \; 0 \le M \le 1000$。

算法分析

令$f_i$表示输入前$i$个单词的最小代价,$s_i=\sum_{j=1}^i C_j$。考虑第$i$个单词的所在行,可以得到转移方程
$$f_i=\min(f_j+(s_i-s_j)^2+M)$$
假设$k \lt j \lt i$,且$j$这个决策更优,即
$$f_j+(s_i-s_j)^2+M \le f_k+(s_i-s_k)^2+M$$
化简得
$$
\begin{align}
f_j-2s_is_j+s_j^2 &\le f_k-2s_is_k+s_k^2 \\
f_j+s_j^2-(f_k+s_k^2) &\le 2s_i(s_j-s_k) \\
{f_j+s_j^2-(f_k+s_k^2) \over 2s_j-2s_k} &\le s_i
\end{align}
$$
也就是说,在满足这个不等式时,对于状态$i$,决策$j$比决策$k$更优。
我们知道$s_i$单调不递减,即如果对于状态$i_0$,$k \lt j \lt i$且决策$j$更优,那么对于状态$i_1 \gt i_0$,决策$j$也一定更优,因此我们可以抛弃决策$k$。
对于每个$i$,假设平面中有点$(2s_i, f_i+s_i^2)$。对于不等式左侧这个东西,把它看成平面直角坐标系里直线的斜率。由于要最小化$f_i$,因此需要找到一个$k$使得不存在$k \lt j \lt i$满足不等式。可以发现我们只要用单调队列维护前$(i-1)$个点的下凸壳(斜率单调递增)。这样时间复杂度就是$O(N)$的了。

代码

#include <algorithm>
#include <cstdio>
#include <cstring>

static int const N = 500005;
int a[N], s[N], f[N], que[N];

int detx(int i, int j) { return s[j] - s[i] << 1; }

int dety(int i, int j) { return f[j] + s[j] * s[j] - f[i] - s[i] * s[i]; }

int main() {
  for (int n, m; ~scanf("%d%d", &n, &m);) {
    for (int i = 1; i <= n; ++i)
      scanf("%d", &a[i]);
    for (int i = 1; i <= n; ++i)
      s[i] = s[i - 1] + a[i];
    int qb = 0, qe = 1;
    for (int i = 1; i <= n; ++i) {
      for (; qb + 1 < qe &&
             dety(que[qb], que[qb + 1]) <= s[i] * detx(que[qb], que[qb + 1]);
           ++qb)
        ;
      f[i] = f[que[qb]] + (s[i] - s[que[qb]]) * (s[i] - s[que[qb]]) + m;
      for (; qb + 1 < qe &&
             dety(que[qe - 2], que[qe - 1]) * detx(que[qe - 1], i) >=
                 dety(que[qe - 1], i) * detx(que[qe - 2], que[qe - 1]);
           --qe)
        ;
      que[qe++] = i;
    }
    printf("%d\n", f[n]);
  }
  return 0;
}

Ellipse

题目描述

Math is important!! Many students failed in 2+2’s mathematical test, so let’s AC this problem to mourn for our lost youth..
Look at this sample picture:
A Ellipse
An ellipse in a plane centered on point $O$. The $L, R$ lines will be vertical through the $X$-axis. The problem is calculating the blue intersection area. But calculating the intersection area is dull, so I have to turn to you, a talent of programmer. Your task is telling me the result of calculation. (defined $\pi=3.14159265$, the area of an ellipse $A=\pi ab$)

题意概述

给定$a, b, l, r$,求椭圆${x^2 \over a^2}+{y^2 \over b^2}=1$在直线$x=l$与$x=r$之间的面积。
数据范围:$-a \le l \le r \le a$。

算法分析

这是一个不规则图形,无法直接计算,但我们可以用积分来求面积。设$f(n)$表示直线$x=n$与椭圆的两个交点之间的距离。
根据辛普森公式
$$
\int_a^b f(x) \, {\rm d}x \approx {b-a \over 6} \left(f(a)+4f\left({a+b \over 2}\right)+f(b)\right)
$$
但这样的精度并不是很高。考虑二分,令
$$g(l, r)=\int_l^r f(x) \, {\rm d}x, \; h(l, r)={r-l \over 6} \left(f(l)+4f\left({l+r \over 2}\right)+f(r)\right)$$
当$h(l, r)$与$h\left(l, {l+r \over 2}\right)+h\left({l+r \over 2}, r\right)$的差足够小时,令
$$g(l, r)=h(l, r)$$
否则
$$g(l, r)=g\left(l, {l+r \over 2}\right)+g\left({l+r \over 2}, r\right)$$
这样就可以达到精度要求。

代码

/*
 * Today is what happened to yesterday.
 */

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>

static double const EPS = 1e-8;
int a, b;

double get_y(double x) { return 2 * b * sqrt((1 - x * x / a / a)); }

double get_s(double l, double r) {
  return (get_y(l) + 4 * get_y((l + r) / 2) + get_y(r)) * (r - l) / 6;
}

double calc(double l, double r) {
  double mid = (l + r) / 2;
  if (fabs(get_s(l, r) - get_s(l, mid) - get_s(mid, r)) < EPS)
    return get_s(l, r);
  return calc(l, mid) + calc(mid, r);
}

int main() {
  int T;
  scanf("%d", &T);
  for (; T--;) {
    int l, r;
    scanf("%d%d%d%d", &a, &b, &l, &r);
    printf("%.3lf\n", calc(l, r));
  }
  return 0;
}

Clever Y

题目描述

Little Y finds there is a very interesting formula in mathematics:
$$
X^Y \bmod Z=K
$$
Given $X, Y, Z$, we all know how to figure out $K$ fast. However, given $X, Z, K$, could you figure out $Y$ fast?

题意概述

给定$X, Z, K$,求满足$X^Y \bmod Z=K$的$Y$的最小值。若不存在$0 \le Y \lt Z$满足条件,则输出无解。
数据范围:$0 \le X, Z, K \le 10^9$。

算法分析

显然$K \ge Z$时无解。在其他情况下,条件可以写成
$$
X^Y \equiv K \pmod Z
$$
特殊处理两种情况:$X, K$中至少有一个为$0$时,$Y=1$或无解;$K=1$时,$Y=0$。
先考虑$(X, Z)=1$的情况。令$M=\lceil \sqrt{Z} \rceil, \; Y=aM-b \; (0 \lt a, b \le M)$,那么
$$
\begin{align}
X^{aM-b} &\equiv K \pmod Z \\
X^{aM} &\equiv K \cdot X^b \pmod Z
\end{align}
$$
可以在$O(M)$的时间内预处理出$K \cdot X^b$,存在哈希表中,接着$O(M)$枚举$a$,判断哈希表中是否存在$X^{aM}$,若存在,则找到解$Y=aM-b$。由于要求最小值,因此哈希表中有重复时应记录较大的$b$。
下面考虑$(X, Z) \gt 1$的情况。令$D=(X, Z)$,易知若$K \nmid D$则无解。设$X=Dx, \; K=Dk, \; Z=Dz$,据同余的性质可得
$$
\begin{align}
(Dx)^Y &\equiv Dk \pmod{Dz} \\
x \cdot (Dx)^{Y-1} &\equiv k \pmod z
\end{align}
$$
因为$D=(X, Z)$,所以$(x, z)=1$,即$x$在模$z$意义下有逆元,因此
$$
(Dx)^{Y-1} \equiv k \cdot x^{-1} \pmod z
$$
令$X’=Dx=X, \; Y’=Y-1, \; K’=k \cdot x^{-1}, \; Z’=z$,则
$$
(X’)^{Y’} \equiv K’ \pmod{Z’}
$$
这和初始时的方程具有一样的形式。若$(X’, Z’) \gt 1$,则重复以上过程,否则可以参照$(X, Z)=1$的情况。

代码

/*
 * You will be singled out for promotion in your work.
 */

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>

class HashTable {
private:
  static int const MOD = 3214567;
  int numr, h[MOD], id[MOD], val[MOD], nxt[MOD];

public:
  void clear() { memset(h, numr = 0, sizeof h); }

int count(int const &n) {
    for (int i = h[n % MOD]; i; i = i[nxt])
      if (i[id] == n)
        return 1;
    return 0;
  }

int &operator[](int const &n) {
    for (int i = h[n % MOD]; i; i = i[nxt])
      if (i[id] == n)
        return i[val];
    ++numr, numr[id] = n, numr[nxt] = h[n % MOD], h[n % MOD] = numr;
    return numr[val];
  }
} mp;

int ex_gcd(int a, int b, int &x, int &y) {
  if (!b) {
    x = 1, y = 0;
    return a;
  }
  int ret = ex_gcd(b, a % b, y, x);
  y -= a / b * x;
  return ret;
}

int get_gcd(int a, int b) {
  int x, y;
  return ex_gcd(a, b, x, y);
}

int get_inv(int a, int b) {
  int x, y;
  return ex_gcd(a, b, x, y), (x + b) % b;
}

int bsgs(int a, int b, int c) {
  a %= c;
  if (a == 0)
    return b == 0 ? 1 : -1;
  if (b == 0)
    return -1;
  if (b == 1)
    return 0;
  int k = 0;
  for (int t; (t = get_gcd(a, c)) > 1; ++k) {
    if (b % t)
      return -1;
    c /= t, b /= t, b = 1ll * b * get_inv(a / t, c) % c;
  }
  mp.clear();
  int d = a, m = sqrt(c) + 1;
  for (int i = 1; i <= m; ++i, d = 1ll * d * a % c)
    mp[1ll * d * b % c] = i;
  d = 1ll * d * get_inv(a, c) % c;
  int e = d;
  for (int i = 1; i <= m; ++i, e = 1ll * e * d % c)
    if (mp.count(e))
      return i * m - mp[e] + k;
  return -1;
}

int main() {
  for (int a, b, c; scanf("%d%d%d", &a, &c, &b), a || b || c;) {
    int ans = bsgs(a, b, c);
    if (~ans)
      printf("%d\n", ans);
    else
      puts("No Solution");
  }
  return 0;
}

Mathjax in WordPress

Background

Mathjax允许我们在Blog中插入漂亮的数学公式。WordPress里有许多Mathjax插件,但大多都不尽如人意,总是缺少这样那样的功能。下面介绍一种不使用插件的方法,可以近乎完美地解决问题。

Steps

  • 点击仪表盘左侧“外观”选项卡中的“编辑”选项,进入“编辑主题”;
  • 在右侧的“主题文件”中,选择“主题页眉”(header.php);
  • 找到其中的<head>标签,在下面插入如下代码;

<script type="text/x-mathjax-config">
  MathJax.Hub.Config({
    extensions: ["tex2jax.js"],
    jax: ["input/TeX", "output/HTML-CSS"],
    tex2jax: {
      inlineMath: [ ['\$','\$'], ["\\(","\\)"] ],
      displayMath: [ ['\$\$','\$\$'], ["\\[","\\]"] ],
      processEscapes: true
    },
    "HTML-CSS": { fonts: ["TeX"] }
  });
</script>

<script type="text/javascript" src="https://cdn.bootcss.com/mathjax/2.7.4/latest.js?config=TeX-MML-AM_CHTML"></script>

  • 点击更新文件,完成后回到站点,查看Mathjax效果;
  • 根据自己的需要修改配置(以上代码适用于大部分情况)。

Known Issues

  • 每次更换主题后需要重复以上步骤。

Ceizenpok’s Formula

题目描述

Dr. Ceizenp’ok from planet i1c5l became famous across the whole Universe thanks to his recent discovery – the Ceizenpok’s formula. This formula has only three arguments: $n, k$ and $m$, and its value is a number of $k$-combinations of a set of $n$ modulo $m$.
While the whole Universe is trying to guess what the formula is useful for, we need to automate its calculation.

题意概述

求${n \choose k} \bmod m$。
数据范围:$1 \le n \le 10^{18}, \; 0 \le k \le n, \; 2 \le m \le 10^6$。

算法分析


$$
x={n \choose k}, \; m=p_1^{a_1}p_2^{a_2} \cdots p_q^{a_q}
$$
可以列出方程组
$$
\left\{
\begin{array}{c}
x \equiv r_1 \pmod{p_1^{a_1}} \\
x \equiv r_2 \pmod{p_2^{a_2}} \\
\cdots \\
x \equiv r_q \pmod{p_q^{a_q}}
\end{array}
\right.
$$
由于模数两两互质,所以该方程组在模$m$意义下有唯一解。
考虑如何求$r_i$。实际上,我们要求的就是${n \choose k} \bmod p_i^{a_i}$。我们知道
$$
{n \choose k}={n! \over k!(n-k)!}
$$
那么只要求出$n! \bmod p_i^{a_i}, \; k! \bmod p_i^{a_i}, \; (n-k)! \bmod p_i^{a_i}$的值,就可以用逆元求出${n \choose k} \bmod p_i^{a_i}$。
对于如何求$n! \bmod p_i^{a_i}$,令
$$
f(n)=n! \bmod p_i^{a_i}
$$
由$x \equiv x+p_i^{a_i} \pmod{p_i^{a_i}}$,可得
$$
f(n)=\left(f\left(\left\lfloor{n \over p_i}\right\rfloor\right)\cdot p_i^{\lfloor n/p_i \rfloor}\cdot\left(\prod_{i \in [1, p_i^{a_i}], \; p_i\not\mid i}i\right)^{\lfloor n/p_i^{a_i} \rfloor}\cdot\prod_{i \in [1, n \bmod p_i^{a_i}], \; p_i\not\mid i}i\right)\bmod p_i^{a_i}
$$
但是$k!, \; (n-k)!$在模$p_i^{a_i}$意义下可能不存在逆元,因此需要将$n!, \; k!, \; (n-k)!$中的$p_i$因子提取出来,求出逆元后再乘回去。
这样就得到了所有$r_i$。用中国剩余定理求解方程组即可。

代码

/*
 * Your lucky number has been disconnected.
 */

#include <algorithm>
#include <cstdio>
#include <cstring>

#define int long long

int power(int a, int b, int m) {
  int ret = 1;
  for (; b; b >>= 1)
    b & 1 && ((ret *= a) %= m), (a *= a) %= m;
  return ret;
}

void ex_gcd(int a, int b, int &x, int &y) {
  if (!b) {
    x = 1, y = 0;
    return;
  }
  ex_gcd(b, a % b, y, x), y -= a / b * x;
}

int get_inv(int a, int b) {
  int x, y;
  ex_gcd(a, b, x, y);
  return (x + b) % b;
}

int get_fac(int n, int p, int k) {
  int m = power(p, k, 1e9), ret = 1;
  for (; n; n /= p) {
    if (n / m) {
      int rec = 1;
      for (int i = 2; i < m; ++i)
        if (i % p)
          (rec *= i) %= m;
      (ret *= power(rec, n / m, m)) %= m;
    }
    for (int i = n % m; i > 1; --i)
      if (i % p)
        (ret *= i) %= m;
  }
  return ret;
}

int calc(int N, int K, int M, int p, int k) {
  int a = get_fac(N, p, k), b = get_fac(K, p, k), c = get_fac(N - K, p, k),
      cnt = 0;
  for (int i = N; i; i /= p)
    cnt += i / p;
  for (int i = K; i; i /= p)
    cnt -= i / p;
  for (int i = N - K; i; i /= p)
    cnt -= i / p;
  int m = power(p, k, 1e9),
      ret = a * get_inv(b, m) % m * get_inv(c, m) % m * power(p, cnt, m) % m;
  return ret * (M / m) % M * get_inv(M / m, m) % M;
}

signed main() {
  int N, K, M, ans = 0;
  scanf("%lld%lld%lld", &N, &K, &M);
  for (int i = 2, t = M; t > 1; ++i)
    if (!(t % i)) {
      int k = 0;
      for (; !(t % i); ++k, t /= i)
        ;
      (ans += calc(N, K, M, i, k)) %= M;
    }
  printf("%lld\n", ans);
  return 0;
}

Strange Way to Express Integers

题目描述

Elina is reading a book written by Rujia Liu, which introduces a strange way to express non-negative integers. The way is described as following:

  • Choose $k$ different positive integers $a_1, a_2, \ldots, a_k$. For some non-negative $m$, divide it by every $a_i$ to find the remainder $r_i$. If $a_1, a_2, \ldots, a_k$ are properly chosen, $m$ can be determined, then the pairs $(a_i, r_i)$ can be used to express $m$.

“It is easy to calculate the pairs from $m$,” said Elina. “But how can I find $m$ from the pairs?”
Since Elina is new to programming, this problem is too difficult for her. Can you help her?

题意概述

给定同余方程组
$$
\left\{
\begin{array}{c}
x \equiv r_1 \pmod{a_1} \\
x \equiv r_2 \pmod{a_2} \\
\cdots \\
x \equiv r_k \pmod{a_k}
\end{array}
\right.
$$
求$x$的最小非负整数解。
数据范围:所有输入输出均可以表示为64位整型。

算法分析

对于两个方程
$$
\begin{align}
x &\equiv r_1 \pmod{a_1} \\
x &\equiv r_2 \pmod{a_2}
\end{align}
$$
考虑将它们合并。易知
$$
x=r_1+k_1a_1=r_2+k_2a_2
$$
移项得
$$
k_1a_1=r_2-r_1+k_2a_2
$$
由贝祖定理可知,这个方程有整数解的充要条件为$(a_1, a_2) \mid r_2-r_1$。方程两边同时除以$(a_1, a_2)$
$$
\begin{align}
k_1{a_1 \over (a_1, a_2)}&={r_2-r_1 \over (a_1, a_2)}+k_2{a_2 \over (a_1, a_2)} \\
k_1{a_1 \over (a_1, a_2)} &\equiv {r_2-r_1 \over (a_1, a_2)} \pmod{{a_2 \over (a_1, a_2)}} \\
k_1 &\equiv \left( {a_1 \over (a_1, a_2)} \right)^{-1} \cdot {r_2-r_1 \over (a_1, a_2)} \pmod{{a_2 \over (a_1, a_2)}}
\end{align}
$$
将其代回$x=r_1+k_1a_1$,得
$$
x \equiv r_1+k_1a_1 \pmod{{a_1a_2 \over (a_1, a_2)}}
$$
至此已成功将两个方程合并为一个相同形式的方程。

代码

/*
 * Your life would be very empty if you had nothing to regret.
 */

#include <algorithm>
#include <cstdio>
#include <cstring>

#define int long long

int ex_gcd(int a, int b, int &x, int &y) {
  if (!b) {
    x = 1, y = 0;
    return a;
  }
  int ret = ex_gcd(b, a % b, y, x);
  y -= a / b * x;
  return ret;
}

int get_gcd(int a, int b) {
  int x, y;
  return ex_gcd(a, b, x, y);
}

int get_inv(int a, int b) {
  int x, y;
  return ex_gcd(a, b, x, y), (x += b) %= b;
}

signed main() {
  for (int k, a, r; ~scanf("%lld%lld%lld", &k, &a, &r);) {
    int flag = 0;
    for (int ai, ri; --k;) {
      scanf("%lld%lld", &ai, &ri);
      int gcd = get_gcd(a, ai);
      if ((r - ri) % gcd)
        flag = 1;
      r += get_inv(a / gcd, ai / gcd) * ((ri - r) / gcd) % (ai / gcd) * a;
      (a /= gcd) *= ai, ((r %= a) += a) %= a;
    }
    if (flag)
      puts("-1");
    else
      printf("%lld\n", r);
  }
  return 0;
}