Primes

题目描述

This is an interactive problem.
For two positive integers $x, y$, we define $\pi(x, y)$ to be the number of distinct primes that divide both $x$ and $y$. For example $\pi(2, 3) = 0, \; \pi(8, 16) = 1$ and $\pi(30, 105) = 2$.
For two positive integers $a, b$, where $a \le b$, we define $S(a, b)$ to be the sum of values $\pi(x, y)$ over all pairs of integers $(x, y)$ satisfying $a \le x \lt y \le b$.
Your task is to compute the values $S(a, b)$ for many query pairs $(a, b)$. To make your task more challenging, all the queries have to be answered online.

题意概述

给定两个整数$a, b$,定义$\pi(x, y)$为所有整除$x$和$y$的质数的个数,$S(a, b)$为所有满足$a \le x \lt y \le b$的$\pi(x, y)$的和,求$S(a, b)$。有$q$组询问,强制在线。

数据范围:$1 \le q \le 5 \times 10^4, \; 1 \le a \le b \le 10^6$。

算法分析

分别考虑每个质数对答案的贡献。若在区间$[a, b]$中有$k$个数能被质数$p$整除,则$p$对答案的贡献为${k \choose 2}$。

首先筛出所有质数。但是每次询问直接枚举质数会超时。考虑对于一个整数$n$,$\lfloor {n \over i} \rfloor$只有$O(\sqrt{n})$种不同的取值。因此可以对$\lfloor {a \over i} \rfloor$和$\lfloor {b \over i} \rfloor$进行分段枚举。

代码

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

int const N = 1000005;

int tp, prime[N], rec[N], vis[N], sum[N];

void init() {
    for (int i = 2; i < N; ++i) {
        if (!vis[i]) {
            prime[tp++] = i;
        }
        for (int j = 0; j < tp && i * prime[j] < N; ++j) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                break;
            }
        }
    }
    for (int i = 2; i < N; ++i)
        sum[i] = sum[i - 1] + !vis[i];
}

int main() {
    init();
    int q;
    scanf("%d", &q);
    for (; q--;) {
        int a, b;
        scanf("%d%d", &a, &b);
        long long ans = 0;
        for (int i = 0; i < tp && prime[i] < b;) {
            int cnt = b / prime[i] - (a - 1) / prime[i];
            int nxt = 1e9;
            if (b / prime[i]) nxt = std::min(nxt, b / (b / prime[i]));
            if ((a - 1) / prime[i]) nxt = std::min(nxt, (a - 1) / ((a - 1) / prime[i]));
            ans += 1ll * cnt * (cnt - 1) / 2 * (sum[nxt] - i);
            i = sum[nxt];
        }
        printf("%lld\n", ans);
        fflush(stdout);
    }
    return 0;
}

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

Counting Divisors (Square)

题目描述

Let $\sigma_0(n)$ be the number of positive divisors of $n$.
For example, $\sigma_0(1)=1, \sigma_0(2)=2$ and $\sigma_0(6)=4$.
Let $S_2(n)=\sum_{i=1}^n \sigma_0(i^2)$.
Given $N$, find $S_2(N)$.

题意概述

令$\sigma_0(n)$表示$n$的正因数个数,$S_2(n)=\sum_{i=1}^n \sigma_0(i^2)$。给定$N$,求$S_2(N)$。
数据范围:$1 \le N \le 10^{12}$。

算法分析

为了方便,我们令$f(n)=\sigma_0(n^2)$,$\omega(n)$表示$n$的不同质因子个数。
设$n=\prod_{k=1}^{\omega(n)} p_k^{a_k}$,其中$p_k$表示不同质数。
先来推$f(n)$。对于$n$的每个因数$d$,将$d$的某些质因子$p_k$的次数加上$a_k$得到$d’$,有$d’ \nmid n, d’ \mid n^2$。这样的方案数有$(2^{\omega(d)}-1)$种,再算上$d$本身,可得
$$
\begin{align}
f(n)=\sum_{d \mid n} 2^{\omega(d)}
\end{align}
$$
化简$2^{\omega(d)}$。它等价于$d$的每个质因数要么选一个,要么不选。那么它就等于$d$的因数中square-free的数的个数。因此
$$
\begin{align}
2^{\omega(d)}=\sum_{i \mid d} \mu^2(i)
\end{align}
$$
回到原式
$$
\begin{align}
S_2(n)=\sum_{i=1}^n \sum_{j \mid i} \sum_{k \mid j} \mu^2(k)
\end{align}
$$
还是毫无头绪。我们继续化简$f(n)$。考虑它的卷积形式$f=\mu^2 \ast 1 \ast 1=\mu^2 \ast \sigma_0$。所以
$$
\begin{align}
S_2(n)&=\sum_{i=1}^n \sum_{j \mid i} \mu^2(j) \sigma_0\left({i \over j}\right) \\
&=\sum_{ij \le n} \mu^2(i) \sigma_0(j)
\end{align}
$$
那么只要用反演求出$\mu^2$的前缀和、用分块求出$\sigma_0$的前缀和,就可以利用分块的技巧了。具体如下
$$
\begin{align}
\sum_{i=1}^n \mu^2(i)&=\sum_{i=1}^n [i\text{最大平方因子}=1] \\
&=\sum_{i=1}^n \sum_{j^2 \mid i} \mu(j) \\
&=\sum_{i=1}^{\lfloor \sqrt{n} \rfloor} \mu(i) \left\lfloor {n \over i^2} \right\rfloor \\
\sum_{i=1}^n \sigma_0(i)&=\sum_{i=1}^n \left\lfloor {n \over i} \right\rfloor
\end{align}
$$

代码

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

void read(long long &t) {
  char c; while ((c = getchar()) < '0' || c > '9') ; t = c - '0';
  while ((c = getchar()) >= '0' && c <= '9') (t *= 10) += c - '0';
}

int N;
char mui[100000000];
int d[100000000], sum[100000000], prime[10000000];
bool vis[100000000];
long long n, sigma[100000000];

void init() {
  int top = 0; mui[1] = 1, sigma[1] = 1;
  for (int i = 2; i < N; ++ i) {
    if (! vis[i]) prime[top ++] = i, mui[i] = -1, sigma[i] = 2, d[i] = 1;
    for (int j = 0; j < top; ++ j) {
      int k = i * prime[j]; if (k >= N) break;
      vis[k] = true;
      if (i % prime[j]) mui[k] = -mui[i], sigma[k] = sigma[i] << 1, d[k] = 1;
      else { mui[k] = 0, sigma[k] = sigma[i] / (d[i] + 1) * (d[i] + 2), d[k] = d[i] + 1; break; }
    }
  }
  for (int i = 1; i < N; ++ i) sum[i] = sum[i - 1] + mui[i] * mui[i], sigma[i] += sigma[i - 1];
}

inline long long get_mui(long long n) {
  if (n < N) return sum[n];
  int q = sqrt(n); long long ret = 0;
  for (register int i = 1; i <= q; ++ i) if (mui[i]) ret += mui[i] * (n / i / i);
  return ret;
}

inline long long get_sigma(long long n) {
  if (n < N) return sigma[n];
  long long ret = 0;
  for (register long long i = 1, j; i <= n; i = j + 1)
    j = n / (n / i), ret += (j - i + 1) * (n / i);
  return ret;
}

long long calc(long long n) {
  int q = sqrt(n); long long ret = 0;
  for (register int i = 1; i <= q; ++ i) if (mui[i]) ret += get_sigma(n / i);
  for (register long long i = q + 1, j, cur, last = sum[q]; i <= n; i = j + 1, last = cur)
    j = n / (n / i), cur = get_mui(j), ret += (cur - last) * get_sigma(n / i);
  return ret;
}

signed main() {
  long long T; read(T);
  if (T > 800) N = 20000; else N = 100000000;
  init();
  while (T --) read(n), printf("%lld\n", calc(n));
  return 0;
}

Nearly Prime Numbers

题目概述

Nearly prime number is an integer positive number for which it is possible to find such primes $P_1$ and $P_2$ that given number is equal to $P_1 \times P_2$. There is given a sequence on $N$ integer positive numbers, you are to write a program that prints “Yes” if given number is nearly prime and “No” otherwise.

题意概述

给定$N$个数$a_i$,分别判断它们是否仅由两个质数相乘得到。
数据范围:$1 \le N \le 10, \; 1 \le a_i \le 10^9$。

算法分析

可以用Miller-Rabin算法快速判断一个数是否是质数。
Miller-Rabin基于以下几个事实:

  • 费马小定理:对于素数$p$和任意整数$a$,有$a^p \equiv a \pmod p$。反之,满足$a^p \equiv a \pmod p$,$p$也几乎一定是素数。
  • 伪素数:如果$n$是一个正整数,且存在和$n$互素的正整数$a$满足$a^{n-1} \equiv 1 \pmod n$,那么$n$是基于$a$的伪素数。如果一个数是伪素数,那么它几乎肯定是素数。
  • 二次探测定理:如果$p$是奇素数,则$x^2 \equiv 1 \pmod p$的解为$x \equiv \pm 1 \pmod p$。

代码

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

#define int long long
#define random(x) ((1ll * rand() << 30 ^ 1ll * rand() << 15 ^ rand()) % (x) + 1)

using namespace std;

struct IOManager {
  template <typename T> inline bool read(T &x) { char c = '\n'; bool flag = false; x = 0; while (~ c && ! isdigit(c = getchar()) && c != '-') ; c == '-' && (flag = true, c = getchar()); if (! ~ c) return false; while (isdigit(c)) x = (x << 3) + (x << 1) + c - '0', c = getchar(); return (flag && (x = -x), true); }
  inline bool read(char &c) { c = '\n'; while (~ c && ! (isprint(c = getchar()) && c != ' ')) ; return ~ c; }
  inline int read(char s[]) { char c = '\n'; int len = 0; while (~ c && ! (isprint(c = getchar()) && c != ' ')) ; if (! ~ c) return 0; while (isprint(c) && c != ' ') s[len ++] = c, c = getchar(); return (s[len] = '\0', len); }
  template <typename T> inline IOManager operator > (T &x) { read(x); return *this; }
  template <typename T> inline void write(T x) { if (x < 0) putchar('-'), x = -x; if (x > 9) write(x / 10); putchar(x % 10 + '0'); }
  inline void write(char c) { putchar(c); }
  inline void write(char s[]) { int pos = 0; while (s[pos] != '\0') putchar(s[pos ++]); }
  template <typename T> inline IOManager operator < (T x) { write(x); return *this; }
} io;

struct Solver {
private:
  int n;
  void input() { io > n; }
  int multiple(int a, int b, int mod) {
    int ret = 0; a %= mod;
    while (b) { b & 1 && ((ret += a) %= mod), (a <<= 1) %= mod, b >>= 1; }
    return ret;
  }
  int power(int a, int b, int mod) {
    int ret = 1; a %= mod;
    while (b) { b & 1 && ((ret *= a) %= mod), (a *= a) %= mod, b >>= 1; }
    return ret;
  }
  bool is_prime(int n) {
    if (n == 2) return true;
    if (n < 2 || ! (n & 1)) return false;
    int c = 0, p = n - 1, cnt = 30;
    while (! (p & 1)) ++ c, p >>= 1;
    while (cnt --) {
      int a = random(n - 1), x = power(a, p, n);
      for (int i = 0; i < c; ++ i) {
        int y = multiple(x, x, n);
        if (y == 1 && x != 1 && x != n - 1) return false; x = y;
      }
      if (x != 1) return false;
    }
    return true;
  }
  void process() {
    int a; io > a; int q = sqrt(a);
    for (int i = 2; i <= q; ++ i)
      if (! (a % i) && is_prime(i) && is_prime(a / i)) { io < (char *) "Yes\n"; return; }
    io < (char *) "No\n";
  }

public:
  void solve() { input(); while (n --) process(); }
} solver;

signed main() {
  solver.solve();
  return 0;
}

Arpa and a List of Numbers

题目描述

Arpa has found a list containing $n$ numbers. He calls a list bad if and only if it is not empty and gcd of numbers in the list is $1$.
Arpa can perform two types of operations:

  • Choose a number and delete it with cost $x$.
  • Choose a number and increase it by $1$ with cost $y$.

Arpa can apply these operations to as many numbers as he wishes, and he is allowed to apply the second operation arbitrarily many times on the same number.
Help Arpa to find the minimum possible cost to make the list good.

题意概述

给定一个长度为$n$的序列,第$i$个数字为$a_i$。有两种操作:删除一个数,消耗$x$;将一个数加$1$,消耗$y$。求使得序列中所有数的最大公约数大于$1$的最少消耗。
数据范围:$1 \le n \le 5 \times 10^5, \; 1 \le x, y \le 10^9, \; 1 \le a_i \le 10^6$。

算法分析

枚举$[2, 10^6]$中的素数$p$作为最大公约数(合数作为最大公约数的消耗一定不比质数少)。对于序列中的每个数,要么将它删去,要么将它变成不小于它且最小的$p$的倍数。设$d_i=\left\lceil {a_i \over p} \right\rceil \times p -a_i$。显然,若$x \le d_i \times y$,那么将它删去更优,否则将它变成$p$的倍数更优。也就是说,对于$d_i \ge {x \over y}$的数,我们将它删去,否则将它变成$p$的倍数。预处理出前缀和,根据调和级数,可知时间复杂度为$O(a_i\log a_i)$。

代码

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

using namespace std;

void minify(long long &a, long long b) { b < a && (a = b); }

struct IOManager {
  template <typename T>
  inline bool read(T &x) {
    char c; bool flag = false; x = 0;
    while (~ c && ! isdigit(c = getchar()) && c != '-') ;
    c == '-' && (flag = true, c = getchar());
    if (! ~ c) return false;
    while (isdigit(c)) x = (x << 3) + (x << 1) + c - '0', c = getchar();
    return (flag && (x = -x), true);
  }
  inline bool read(char &c) {
    c = '\n';
    while (~ c && ! (isprint(c = getchar()) && c != ' ')) ;
    return ~ c;
  }
  inline int read(char s[]) {
    char c; int len = 0;
    while (~ c && ! (isprint(c = getchar()) && c != ' ')) ;
    if (! ~ c) return 0;
    while (isprint(c) && c != ' ') s[len ++] = c, c = getchar();
    return (s[len] = '\0', len);
  }
  template <typename T>
  inline IOManager operator > (T &x) {
    read(x); return *this;
  }
  template <typename T>
  inline void write(T x) {
    x < 0 && (putchar('-'), x = -x);
    x > 9 && (write(x / 10), true);
    putchar(x % 10 + '0');
  }
  inline void write(char c) {
    putchar(c);
  }
  inline void write(char s[]) {
    int pos = 0;
    while (s[pos] != '\0') putchar(s[pos ++]);
  }
  template <typename T>
  inline IOManager operator < (T x) {
    write(x); return *this;
  }
} io;

struct Solver {
private:
  static const int N = 1000000;
  static const int M = 1000000;
  int n, x, y, a[N + 1], cnt[N + 1 << 1];
  long long sum[N + 1 << 1];
  int top, prime[N];
  bool vis[N + 1 << 1];
  void input() {
    io > n > x > y;
    for (int i = 1; i <= n; ++ i) io > a[i], ++ cnt[a[i]], sum[a[i]] += a[i];
  }
  void init() {
    for (int i = 1; i <= M << 1; ++ i) cnt[i] += cnt[i - 1], sum[i] += sum[i - 1];
    for (int i = 2; i <= M; ++ i) {
      if (! vis[i]) prime[++ top] = i;
      for (int j = 1; j <= top && i * prime[j] <= M; ++ j) {
        vis[i * prime[j]] = true;
        if (! (i % prime[j])) break;
      }
    }
  }
  void process() {
    int p = x / y;
    long long ans = 1000000000000000000ll;;
    for (int i = 1; i <= top; ++ i) {
      long long tmp = 0;
      int delta = prime[i] - p - 1;
      for (int j = 0; j <= N; j += prime[i]) {
        if (delta < 0) {
          tmp += (((long long) cnt[j + prime[i]] - cnt[j]) * (j + prime[i]) - (sum[j + prime[i]] - sum[j])) * y;
        } else {
          tmp += ((long long) cnt[j + delta] - cnt[j]) * x + (((long long) cnt[j + prime[i]] - cnt[j + delta]) * (j + prime[i]) - (sum[j + prime[i]] - sum[j + delta])) * y;
        }
      }
      minify(ans, tmp);
    }
    io < ans < '\n';
  }

public:
  void solve() {
    input(), init(), process();
  }
} solver;

int main() {
  solver.solve();
  return 0;
}

Vasya’s Function

题目描述

Vasya is studying number theory. He has denoted a function $f(a, b)$ such that:

  • $f(a, 0)=0$;
  • $f(a, b)=1+f(a, b-(a, b))$, where $(a, b)$ is the greatest common divisor of $a$ and $b$.

Vasya has two numbers $x$ and $y$, and he wants to calculate $f(x, y)$. He tried to do it by himself, but found out that calculating this function the way he wants to do that might take very long time. So he decided to ask you to implement a program that will calculate this function swiftly.

题意概述

函数$f(a, b)$的定义如下:
$$
\begin{align}
f(a, 0)&=0 \\
f(a, b)&=1+f(a, b-(a, b))
\end{align}
$$
给定$x, y$,求$f(x, y)$。
数据范围:$1 \le x, y \le 10^{12}$。

算法分析

有个显然的结论:$f(a, b)=f(ka, kb) \; (k \gt 0)$。所以只需要考虑$(a, b)=1$的情况。
设$(x, y)=1$,由于$x$始终不变,因此先将$x$分解质因数。质因子最多只有$12$个,因为$2 \times 3 \times 5 \times 7 \times \cdots \times 31 \times 37 \gt 10^{12}$。找出最大的$s$满足$s \lt y$且$s$是$x$某个质因子的倍数,将答案加上$y-s$,并令$y=s$。将$x, y$同时除以$(x, y)$,问题转化成了一个更小的子问题,可以迭代处理。

代码

import std.stdio;
import std.math;

static immutable int N = 13;

uint top;
long x, y, p, q, ans;
long[N] prime, cnt;

long gcd(long a, long b) {
  return b ? gcd(b, a % b) : a;
}

int main() {
  readf(" %d %d", &x, &y);
  p = x, q = cast(long) sqrt(real(x));
  for (int i = 2; i <= q; ++ i) {
    while (! (p % i)) {
      p /= i;
      if (prime[top] == i) ++ cnt[top];
      else prime[++ top] = i, cnt[top] = 1;
    }
  }
  if (p > 1) prime[++ top] = p, cnt[top] = 1;
  while (y) {
    long max = 0;
    for (int i = 1; i <= top; ++ i)
      if (cnt[i])
        if (y / prime[i] * prime[i] > max)
          max = y / prime[i] * prime[i];
    ans += y - max, y = max;
    long g = gcd(x, y);
    x /= g, y /= g;
    for (int i = 1; i <= top; ++ i)
      while (! (g % prime[i]))
        g /= prime[i], -- cnt[i];
  }
  writeln(ans);
  return 0;
}

Notes on Module

逆元

对于模意义下的除法,$a/x \equiv b \pmod m \Leftrightarrow a \times x^{-1} \equiv b \pmod m$,则称$x^{-1}$为$x$在模$m$意义下的逆元。
关于逆元的求法可以参考Notes on Multiplicative Inverse
考虑线性同余方程$ax \equiv b \pmod m$,有
$$(a/(a, m))x \equiv b/(a, m) \pmod {m/(a, m)}$$
显然,当$(a, m) \not \mid b$时,原方程无解;否则,有
$$x \equiv (a/(a, m))^{-1} \times (b/(a, m)) \pmod {m/(a, m)}$$
因此原方程的解为
$$x \equiv (a/(a, m))^{-1} \times (b/(a, m)) + (m/(a, m)) \times k \pmod m \; (0 \le k \lt (a, m))$$
原方程可能有多解,也可能无解。


费马小定理

当$p$是素数时,对任意整数$x$都有$x^p \equiv x \pmod p$。若$p \not \mid x$,则
$$
\begin{align}
x^{p-1} &\equiv 1 \pmod p \\
x^{-1} &\equiv x^{p-2} \pmod p
\end{align}
$$
当$m$不是素数且$x$与$m$互质时,有
$$x^{\varphi(m)} \equiv 1 \pmod m$$
其中$\varphi(m)$为欧拉函数,其值等于不大于$m$的正整数中与$m$互质的数的个数。
利用线性筛,每次发现质数$p$时把它倍数的欧拉函数乘上$(p-1)/p$,就可以求出$1$到$n$的欧拉函数。


线性同余方程组

给定由若干个形如$a_ix \equiv b_i \pmod {m_i}$的方程组成的方程组,求它的解。如果有解,则一定有无穷多解,并可以表示为$x \equiv b \pmod m$的形式,因此问题转化为求解$b$和$m$。如果我们能求解方程组
$$
\left\{
\begin{align}
x &\equiv b_i \pmod {m_i} \\
ax &\equiv b_{i+1} \pmod {m_{i+1}}
\end{align}
\right.
$$
那么只要对方程逐个求解即可得到答案。
因为$x \equiv b_i \pmod {m_i}$,所以$x=tm_i+b_i$,代入第二个式子,得到
$$a(tm_i+b_i) \equiv b_{i+1} \pmod {m_{i+1}}$$
移项后得到
$$atm_i \equiv b_{i+1}-ab_i \pmod {m_{i+1}}$$
这是一个一次同余方程。当$(am_i, m_{i+1}) \not \mid b_{i+1}-ab_i$时原方程组无解。


中国剩余定理

如果同余方程组满足
$$\forall i \in [1, n], \; a_i=1 \land \forall i, j \in [1, n] \land i \neq j, \; (m_i, m_j)=1$$
那么答案的形式一定是
$$x \equiv b \pmod {\prod m_i}$$
反之,对于一个合数$n$,如果$n=ab \land (a, b)=1$,那么$x \bmod n$的值确定了之后,$x \bmod a$和$x \bmod b$的值也就确定了。因此,以$n$为模数来考虑和以$a$、$b$为模数来考虑是等价的。这个定理叫做中国剩余定理。
所以,对于模合数的情况只要考虑模若干$p^k$($p$为素数)的情况就可以了。如果$n$是一个square-free的数,那么问题就转化成模素数的情况,从而变得容易求解。
若$n$不是square-free的数,可以将每个数字$x$用一个二元组$(a, b)$表示,代表$x=ap_i^b$。有以下运算规则
$$
\begin{align}
(a, b)\times(c, d)&=(ac, b+d) \\
(a, b)/(c, d)&=(ac^{-1}, b-d)
\end{align}
$$
关于扩展中国剩余定理可参考Strange Way to Express Integers


$n!$

在计数问题中经常用到$n!$。因此有必要了解$n!$在模$p$意义下的一些性质。
假设$p$是素数,$n!=ap^e \; (a \not \mid p)$,求$a \bmod p$和$e$,$e$是$n!$能够迭代整除$p$的次数。显然$$e=n/p+n/p^2+n/p^3+ \cdots$$
因为$n/d$等于不超过$n$的能被$d$整除的正整数个数。
接下来计算$a \bmod p$。可以发现,不能被$p$整除的项在模$p$意义下呈周期性,它们的积为
$$((p-1)!)^{n/p} \times (n \bmod p)!$$
事实上,根据威尔逊定理,$(p-1)! \equiv -1 \pmod p$。因为除了$1$和$p-1$之外的其余项都可以和各自的逆元相乘得到$1$。
接下来计算可以被$p$整除的项的积。可以被$p$整除的项为$p, 2p, 3p, \ldots, (n/p)p$,将$p$消去后得到$1, 2, 3, \ldots, n/p$。因此问题的范围就由$n$缩小到了$n/p$。


${n \choose k} \bmod p$

当$p$不太大时,根据卢卡斯定理,令
$$n=\sum n_ip^i \; (0 \le n_i \lt p), \; k=\sum k_ip^i \; (0 \le k_i \lt p)$$

$${n \choose k} \equiv \prod {n_i \choose k_i} \pmod p$$
还有另一种思路。首先,把${n \choose k}$写成阶乘的形式
$${n \choose k}={n! \over k!(n-k)!}$$

$$n!=a_1p^{e_1}, \; k!=a_2p^{e_2}, \; (n-k)!=a_3p^{e_3}$$
从中可以看出,当$e_1 \gt e_2+e_3$时,${n \choose k}$可以被$p$整除;当$e_1=e_2+e_3$时,${n \choose k}$无法被$p$整除,此时${n \choose k} \equiv a_1(a_2a_3)^{-1} \pmod p$。
关于扩展卢卡斯定理可参考Ceizenpok’s Formula

Notes on Math Theory

Lucas

$n=(a_ka_{k-1}a_{k-2}\ldots a_0)_p, \; m=(b_kb_{k-1}b_{k-2}\ldots b_0)_p, \; {n \choose m} \equiv \prod_{i=0}^k {a_i \choose b_i} \pmod p$.
$n!$中$p$的次数$=\sum_{i=1}^{\infty} [{n \over p^i}]={n-S_p(n) \over p-1}$
${n \choose m}$中$p$的次数$p^\alpha \mid \mid {n \choose m} \Leftrightarrow p^\alpha \mid {n \choose m} \land p^{\alpha+1} \not \mid {n \choose m}$.
${n \choose m} \rightarrow p \not \mid {n \choose m} \Leftrightarrow {n-S_p(n) \over p-1}={m-S_p(m) \over p-1}+{n-m-S_(n-m) \over p-1} \Leftrightarrow m+(n-m)$在$p$进制下无进位。
$(1+x)^n \rightarrow x^m$系数${n \choose m}$.
$(1+x)^{p^k} \equiv 1+x^{p^k} \pmod p$.
$(1+x)^{a_kp^k+ \cdots +a_1p+a_0}=(1+x)^{a_kp^k} \cdots (1+x)^{a_1p}(1+x)^{a_0}$.
$\equiv (1+x^{p^k})^{a_k} \cdots (1+x^p)^{a_1}(1+x)^{a_0} \pmod p$.
$n!!=1 \times 3 \times 5 \times \cdots \times n$.


证明${n \choose k} (0 \le k \le n)$为奇数的个数是$2$的幂次。


Euler

$x^{\varphi(n)} \equiv 1 \pmod n, (x, n)=1$.
$x^n \equiv x^{n-\varphi(n)} \pmod n$.
$\bmod n$缩系.
$a_1 \times a_2 \times \cdots \times a_{\varphi(n)} \equiv a_1x \times a_2x \times \cdots \times a_{\varphi(n)}x \pmod n$.
$n$为素数,$x^{p-1} \equiv 1 \pmod p, p \not \mid x$.
$x^p \equiv x \pmod p$.
$\sum_{d \mid n} \varphi(d)=n$.
$\varphi(p^{\alpha})=p^{\alpha}-p^{\alpha-1}$.
$n=p_1^{\alpha_1} \cdots p_k^{\alpha_k}$.
$\sum_{\beta_i \le \alpha_i} \varphi(p_1^{\beta_1} \cdots p_k^{\beta_k})=\sum_{\beta_i \le \alpha_i} \prod_{i=1}^k \varphi(p_i^{\beta_i})=\prod_{i=1}^k \sum_{\beta_i \le \alpha_i} \varphi(p_i^{\beta_i})=\prod p_i^{\alpha_i}$.


求满足$n \mid 2^n-1$的$n$的个数。
$n=1$成立。
$n>1$时,设$p$为$n$最小的素因子。
满足$p \mid 2^{d_0}-1$的最小的$d_0$,$\forall p \mid 2^n-1, \; d_0 \mid n$.
$\begin{cases} p \mid 2^n-1\\ p \mid 2^{p-1}-1\\ p \mid 2^d-1 \end{cases} \Rightarrow \begin{cases} d \mid n\\ d \mid p-1 \end{cases} \Rightarrow d=1 \Rightarrow p \mid 1 \Rightarrow n=1$.
$ax+by=(a, b)$.
$nx+(p-1)y=(n, p-1)$.
$p \mid 2^{(n, p-1)}-1$.
∵$(n, p-1)=1$.
∴$p \mid 1$.


$1, x, x^2, \ldots, x^{\varphi(n)-1}$构成$\bmod n$的缩系,则称$x$为$n$的原根,记作$g$。
$p$是素数时,${g, g^2, \ldots, g^{p-1}} \equiv {1, 2, \ldots, p-1} \pmod p$.
$n$的原根有$\varphi(\varphi(n))$个。


给定$k$,求$\sum_{i=1}^{p-1} i^k \pmod p$。
${1, 2, \ldots, p-1} \Rightarrow {a, 2a, \ldots, (p-1)a}$.
$\sum i^k \equiv \sum (a_i)^k \Rightarrow (a^k-1) \sum i^k \equiv 0 \pmod p$.


二项同余式

$x^k \equiv n \pmod p$无解或有$(k, p-1)$个解。
$f(x) \equiv 0 \pmod p$.
$x^k \equiv 1 \pmod 1$有$d$个解,$d=(k, p-1)=ak+b(p-1)$。
$\begin{cases} x^{p-1} \equiv 1 \pmod p\\ x^k \equiv 1 \pmod 1 \end{cases} \Rightarrow x^{(k, p-1)} \equiv 1 \pmod p$至多$d$个解。
${x^{p-1}-1 \over x^d-1}=x^{p-1-d}+x^{p-1-2d}+ \cdots +x^d+1 \equiv 1 \pmod p$至多$p-1-d$个解。
$x^k \equiv 1 \pmod 1$恰有$d$个解。
$x \rightarrow p$的缩系,$x^k$有${p-1 \over (k, p-1)}$个不同的取值。
$x^{p-1} \equiv 1 \pmod p$且$x^d \not \equiv 1 \pmod p, d \lt p-1$,$x$是原根。
令$\gamma(n)$表示次数为$n$的数的个数。$p-1=q_1^{\alpha_1} \cdots q_k^{\alpha_k}, \gamma(q^{\alpha})=q^{\alpha}-q^{\alpha-1}$,且$\gamma$是积性函数。
$x^{q^{\alpha}} \equiv 1 \pmod p$.
$x^{q^{\alpha-1}} \not \equiv 1 \Leftrightarrow \forall d \lt q^{alpha}, \; x^d \not \equiv 1 \pmod p$.
$x^{(q^{\alpha}, d)} \equiv 1 \Leftrightarrow x^{q^{\beta}} \equiv 1 \pmod p$.
$\gamma(l_1l_2)=\gamma(l_1)\gamma(l_2)$.
$h_1$次数是$l_1$,$h_2$次数是$l_2 \Rightarrow h_1h_2$的次数是$l_1l_2$。
$x$次数是$l_1l_2 \Rightarrow x_1$次数是$l_1$,$x_2$次数是$l_2$。


Prime

$\lim \sum_{p=1}^{\infty} (1-{1 \over p})=0$.


$f(x) \equiv 0 \pmod {p^{\alpha}}$的解数$=f(x) \equiv 0 \pmod p$的解数$\Leftarrow \begin{cases} f(x) \equiv 0\\ f'(x) \equiv 0 \end{cases} \pmod p$无解。
$x^2+1 \equiv 0 \pmod {p^{\alpha}}, p \ge 3$无解或恰有$2$解。


勾股数的通式

$a^2+b^2=c^2 (a, b, c \gt 0)$.
$(a, b, c)=1$本原。
$a^2+b^2=c^2 \Rightarrow b^2=c^2-a^2=(c-a)(c+a)$.
$(c-a, c+a)=(2c, c-a)=(2, c-a)=1$.
$2 \not \mid c-a$.
设$b$为偶数,$a$为奇数,则$c$为奇数。
$\begin{cases} {c-a \over 2}=m^2\\ {c+a \over 2}=n^2 \end{cases}$.
$\begin{cases} b=2mn\\ a=n^2-m^2\\ c=n^2+m^2 \end{cases}$.


多项式

拉格朗日差值公式:$f(x)=\sum_{i=1}^n \prod_{j \neq i} {(x-x_j) \over (x_i-x_j)}f(x_i)$.
牛顿差值公式:$f(x)=a_0+a_1(x-x_1)+a_2(x-x_1)(x-x_2)+ \cdots +a_n \prod_{i=1}^n (x-x_i)$.


有理根 Eisenstein

$f \in Z[x], f(x)=a_nx^n+ \cdots +a_1x+a_0$.
$\pm {p \over q} \Rightarrow p \mid a_n, q \mid a_0$.
$f(x)=c_nx^n+c_{n-1}x^{n-1}+ \cdots +c_1x+c_0$可约$\Leftarrow p \mid c_0, p \mid c_1, \ldots, p \mid c_{n-1}, p \not \mid c_n, p^2 \not \mid c_0$.
$f(x)=c_nx^n+c_{n-1}x^{n-1}+ \cdots +c_1x+c_0$可约$\Leftarrow p \mid (c_0, c_1, \ldots, c_{n-1}), p \not \mid c_n, p^2 \not \mid c_0$.
$f(x)=x^{p-1}+x^{p-2}+ \cdots +x+1={x^p-1 \over {x-1}}$.
令$y=x-1$,$f(x)=g(y)={(y+1)^p-1 \over y}=y^{p-1}+{p \choose 1}y^{p-2}+ \cdots +{p \choose p-1}$.


升幂定理

令$V_p(n)$表示$n$中$p$的次数。
$5^{210000}-2^{210000}$里有多少$3$的幂次。$2$次。
$(1+p)^k-1$里有多少$p$的幂次。$1+V_p(k)$次。
$5^{33333}+1$里有多少$3$的幂次。$2$次。
$p^{\alpha} \mid \mid a-b, (a, p)=1, (b, p)=1, p \ge 3$.
$V_p(a^m-b^m)=\alpha+V_p(m)$.
设$p^{\alpha} \mid \mid a^n-b^n$.
$p^{\alpha} \mid \mid a^{pn}-b^{pn}$.
$p^{\alpha} \mid \mid a^{kn}-b^{kn}, (k, p)=1$.


二次剩余

$ax^2+bx+c \equiv 0 \pmod p$.
$x^2 \equiv a \pmod p$.
勒让德符号:$({a \over p})=\begin{cases} 0 & p \mid a\\ 1 & a是p的二次剩余\\ -1 & a不是p的二次剩余 \end{cases}$.
${p-1 \over (2, p-1)}={p-1 \over 2}$.
$1^2, 2^2, 3^2, \ldots, ({p-1 \over 2})^2$.
二次剩余数量$=$非二次剩余数量$={p-1 \over 2}$。
Euler判别法:$n^{p-1 \over 2} \equiv ({n \over p}) \pmod p, (n, p)=1$.
Gauss Lama:$n, 2n, 3n, \ldots, {p-1 \over 2}n \pmod p$的最小正剩余中大于${p-1 \over 2}$的个数为$m$,则$({n \over p})=(-1)^m$。
$({2 \over p})=(-1)^{p^2-1 \over 8}=(-1)^{[{p \over 2}]-[{p \over 4}]}$.
二次互反:$({q \over p})({p \over q})=(-1)^{{p-1 \over 2}{q-1 \over 2}}$.


Summary

$n=(a_ka_{k-1}a_{k-2}\ldots a_0)_p, m=(b_kb_{k-1}b_{k-2}\ldots b_0)_p, {n \choose m} \equiv \prod_{i=0}^k {a_i \choose b_i} \pmod p$.

$x^{\varphi(n)} \equiv 1 \pmod n, (x, n)=1$.
$x^n \equiv x^{n-\varphi(n)} \pmod n$.
$n$为素数,$x^{p-1} \equiv 1 \pmod p, p \not \mid x$.
$x^p \equiv x \pmod p$.

$x^k \equiv n \pmod p$无解或有$(k, p-1)$个解。
$x^k \equiv 1 \pmod 1$有$d$个解,$d=(k, p-1)=ak+b(p-1)$。
$x \Rightarrow p$的缩系,$x^k$有${p-1 \over (k, p-1)}$个不同的取值。
$\bmod p$的原根有$\varphi(p-1)$个。

$\lim \sum_{p=1}^{\infty} (1-{1 \over p})=0$.

$f(x) \equiv 0 \pmod {p^{\alpha}}$的解数$=f(x) \equiv 0 \pmod p$的解数$\Leftarrow \begin{cases} f(x) \equiv 0\\ f'(x) \equiv 0 \end{cases} \pmod p$无解。

$\begin{cases} b=2mn\\ a=n^2-m^2\\ c=n^2+m^2 \end{cases}$.

拉格朗日差值公式$f(x)=\sum_{i=1}^n \prod_{j \neq i} {(x-x_j) \over (x_i-x_j)}f(x_i)$.
牛顿差值公式$f(x)=a_0+a_1(x-x_1)+a_2(x-x_1)(x-x_2)+ \cdots +a_n \prod_{i=1}^n (x-x_i)$.

$f(x)=c_nx^n+c_{n-1}x^{n-1}+ \cdots +c_1x+c_0$可约$\Leftarrow p \mid c_0, p \mid c_1, \ldots, p \mid c_{n-1}, p \not \mid c_n, p^2 \not \mid c_0$.
$f(x)=c_nx^n+c_{n-1}x^{n-1}+ \cdots +c_1x+c_0$可约$\Leftarrow p \mid (c_0, c_1, \ldots, c_{n-1}), p \not \mid c_n, p^2 \not \mid c_0$.

$p^{\alpha} \mid \mid a-b, (a, p)=1, (b, p)=1, p \ge 3$.
$V_p(a^m-b^m)=\alpha+V_p(m)$.

Fox and Dinner

题目描述

Fox Ciel is participating in a party in Prime Kingdom. There are $n$ foxes there (include Fox Ciel). The $i$-th fox is $a_i$ years old.
They will have dinner around some round tables. You want to distribute foxes such that:

  • each fox is sitting at some table,
  • each table has at least 3 foxes sitting around it,
  • the sum of ages of any two adjacent foxes around each table should be a prime number.

If $k$ foxes $f_1, f_2, \ldots, f_k$ are sitting around table in clockwise order, then for $1 \le i \le k-1$: $f_i$ and $f_{i+1}$ are adjacent, and $f_1$ and $f_k$ are also adjacent.
If it is possible to distribute the foxes in the desired manner, find out a way to do that.

题意概述

给定$n$个数$a_i$,你需要将它们分成若干个圈,使得每个圈中至少有$3$个数,且任意两个相邻的数的和为质数。
数据范围:$3 \le n \le 200, \; 2 \le a_i \le 10000$。

算法分析

因为每个数都大于$1$,所以任意两个相邻的数必定是一奇一偶。对于任意一个奇数$a_i$和任意一个偶数$a_j$,如果$a_i+a_j$是质数,那么就由$a_i$向$a_j$连一条流量为$1$的边。由于每个数都必须与另外两个数相连,因此可以由一个假设的源点向每个奇数的点连一条流量为$2$的边,由每个偶数的点向一个假设的汇点连一条流量为$2$的边。找出这张图的最大流,如果源点连出的边和连向汇点的边全部满流则说明有解。最后遍历找出所有环即可。

代码

#include <iostream>
#include <cstring>
#include <cmath>
#include <vector>
using namespace std;
struct edge {
    int v, f, nxt;
} e[50001];
int n, nume = 1, top, src, sink, a[202], h[202], g[202], dist[202], que[202];
bool vis[202];
vector<int> ans[201];
void add_edge(int u, int v, int f) {
    e[++nume].v = v, e[nume].f = f, e[nume].nxt = h[u], h[u] = nume;
    e[++nume].v = u, e[nume].f = 0, e[nume].nxt = h[v], h[v] = nume;
}
void bfs() {
    memset(dist, 0, sizeof(dist));
    que[1] = src, dist[src] = 1;
    int s = 0, t = 1;
    while (s < t) {
        int u = que[++s];
        for (int i = h[u]; i; i = e[i].nxt) {
            if (e[i].f && !dist[e[i].v]) {
                que[++t] = e[i].v, dist[e[i].v] = dist[u] + 1;
            }
        }
    }
}
int dfs(int u, int delta) {
    if (u == sink) return delta;
    int ret = 0;
    for (int i = g[u]; i; i = e[i].nxt) {
        if (e[i].f && dist[e[i].v] == dist[u] + 1) {
            int dd = dfs(e[i].v, min(e[i].f, delta));
            e[i].f -= dd, e[i ^ 1].f += dd;
            if (e[i].f) g[u] = i;
            delta -= dd, ret += dd;
        }
    }
    return ret;
}
int max_flow() {
    int ret = 0;
    while (1) {
        bfs();
        if (!dist[sink]) return ret;
        for (int i = 0; i <= n + 1; ++i) g[i] = h[i];
        ret += dfs(src, 1e9);
    }
}
bool is_prime(int t) {
    int k = (int) sqrt(t);
    for (int i = 2; i <= k; ++i) if (!(t % i)) return false;
    return true;
}
void find(int u, int t) {
    vis[u] = true, ans[t].push_back(u);
    for (int i = h[u]; i; i = e[i].nxt) {
        if (e[i].v && e[i].v != n + 1 && e[i].f == !(a[u] & 1) && !vis[e[i].v]) find(e[i].v, t);
    }
}
int main()
{
    cin >> n;
    src = 0, sink = n + 1;
    for (int i = 1; i <= n; ++i) {
        cin >> a[i];
        if (a[i] & 1) add_edge(src, i, 2);
        else add_edge(i, sink, 2);
    }
    for (int i = 1; i <= n; ++i) {
        for (int j = 1; j <= n; ++j) {
            if (a[i] & 1 && !(a[j] & 1) && is_prime(a[i] + a[j])) add_edge(i, j, 1);
        }
    }
    max_flow();
    for (int i = h[0]; i; i = e[i].nxt) {
        if (e[i].f) {
            cout << "Impossible" << endl;
            return 0;
        }
    }
    for (int i = h[n + 1]; i; i = e[i].nxt) {
        if (e[i].f < 2) {
            cout << "Impossible" << endl;
            return 0;
        }
    }
    for (int i = 1; i <= n; ++i) if (!vis[i]) top++, find(i, top);
    cout << top << endl;
    for (int i = 1; i <= top; ++i) {
        cout << ans[i].size() << ' ';
        for (vector<int>::iterator iter = ans[i].begin(); iter != ans[i].end(); ++iter) {
            cout << *iter << ' ';
        }
        cout << endl;
    }
    return 0;
}