## 题目描述

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?

## 算法分析

$$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}$$

$$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)$$

## 代码

/*
*/

#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);
return 0;
}


## 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]();
}
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

    v
acbacabacbacb
^


 v
acbacbaacbacb
^


void kmp(char c[], int len) {
nxt = -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];
}
}


## Application

KMP算法除了可以用于字符串匹配，还能求出一个循环串的最小循环节。循环串是由至少两个相同的字符串拼接起来得到的字符串。

• 引理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}$。

## 题目描述

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$.

## 算法分析

$$f_i=\min(f_j+(s_i-s_j)-(i-j)b_j \mid i-j \ge T)$$

\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;
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;
}


## 题目描述

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.

## 算法分析

$$f_i=\min(f_j+(s_i-s_j)^2+M)$$

$$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}

## 代码

#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;
}


## 题目描述

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: 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$)

## 算法分析

$$\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)$$

$$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;
}


## 题目描述

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^Y \equiv K \pmod Z$$

\begin{align} X^{aM-b} &\equiv K \pmod Z \\ X^{aM} &\equiv K \cdot X^b \pmod Z \end{align}

\begin{align} (Dx)^Y &\equiv Dk \pmod{Dz} \\ x \cdot (Dx)^{Y-1} &\equiv k \pmod z \end{align}

$$(Dx)^{Y-1} \equiv k \cdot x^{-1} \pmod z$$

$$(X’)^{Y’} \equiv K’ \pmod{Z’}$$

## 代码

/*
* 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;
}


## Background

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

## Steps

• 点击仪表盘左侧“外观”选项卡中的“编辑”选项，进入“编辑主题”；
• 找到其中的<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

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

## 题目描述

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.

## 算法分析

$$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.$$

$${n \choose k}={n! \over k!(n-k)!}$$

$$f(n)=n! \bmod 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}$$

## 代码

/*
* 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;
}


## 题目描述

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.$$

## 算法分析

\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$$

\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 \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;
}