Introduction to Quadratic Residue

定义

令$p$为奇质数。

  • 定义1 对于一个整数$n$,若存在整数$x$,使得$x^2 \equiv n \pmod p$,则称$n$是模$p$的二次剩余,否则称$n$是模$p$的二次非剩余
  • 定义2(勒让德符号) $$\left(\frac{n}{p}\right)=\begin{cases}1, & n是模p的二次剩余 \\ -1, & n是模p的二次非剩余 \\ 0, & p \mid n\end{cases}$$

定理

  • 定理1(欧拉判别准则) $$\left(\frac{n}{p}\right) \equiv n^{\frac{p-1}{2}} \pmod p$$
  • 证明
    若$p \mid n$,则$n \equiv 0 \pmod p$,结论显然成立;
    若$n$是模$p$的二次剩余,则存在$x$满足$x^2 \equiv n \pmod p$,于是$n^{\frac{p-1}{2}} \equiv x^{p-1} \pmod p$,由费马小定理,结论成立;
    若$n$是模$p$的二次非剩余,则不存在$x$满足$x^2 \equiv n \pmod p$。由逆元相关知识,对于任意$1 \le i \le p-1$,存在唯一的$1 \le j \le p-1$,使得$i \neq j$且$ij \equiv n \pmod p$。于是$1$到$p-1$这些数可以两两配对,每一对的乘积都与$n$同余,所以$(p-1)! \equiv n^{\frac{p-1}{2}} \pmod p$,由威尔逊定理,$(p-1)! \equiv -1 \pmod p$,结论成立。

令$1 \le x, y \le p-1$。

  • 定理2 $$x+y \equiv 0 \pmod p \Leftrightarrow x^2 \equiv y^2 \pmod p$$
  • 证明
    $\Rightarrow$:$x+y \equiv 0 \pmod p \Rightarrow x \equiv -y \pmod p$,两边平方即得。
    $\Leftarrow$: $x^2 \equiv y^2 \pmod p \Rightarrow x^2-y^2 \equiv (x-y)(x+y) \equiv 0 \pmod p$,显然$x-y \nmid p$,于是有$x+y \mid p \Rightarrow x+y \equiv 0 \pmod p$。
  • 定理3 在$1$到$p-1$中,模$p$的二次剩余和二次非剩余的个数均为$\frac{p-1}{2}$。
  • 证明
    定理2,$1$到$p-1$中有且仅有$\frac{p-1}{2}$个不同的二次剩余,其余$\frac{p-1}{2}$个即为二次非剩余。并且对于每一个二次剩余$n$,都存在两个不同的$x$满足$x^2 \equiv n \pmod p$。

算法

给定整数$n$和奇质数$p$,求满足$x^2 \equiv n \pmod p$的$x$。

  1. 欧拉判别准则判断$n$是否为模$p$的二次剩余,如果不是则返回$-1$表示无解,如果是$0$则返回$0$;
  2. 从$1$到$p-1$中随机选一个整数$a$,使得$a^2-n$为模$p$的二次非剩余(根据定理3,随机的次数不会太多);
  3. 令$\omega \equiv a^2-n \pmod p$,取$x \equiv (a+\sqrt{\omega})^{\frac{p-1}{2}} \pmod p$,返回$x$(这里$\sqrt{\omega}$可以理解成虚数)。
  • 引理1 $$\omega^{\frac{p}{2}} \equiv -\omega^{\frac{1}{2}} \pmod p$$
  • 证明
    欧拉判别准则,显然成立。
  • 引理2 $$(a+b)^p \equiv a^p+b^p \pmod p$$
  • 证明
    $(a+b)^p \equiv \sum_{i=0}^p {p \choose i}a^ib^{p-i}$,由于$p$是奇质数,所以当且仅当$i=0$或$i=p$时${p \choose i}$没有因子$p$,于是成立。
  • 算法证明
    $$\begin{align} x^2 & \equiv (a+\sqrt{\omega})^{p-1} \\ & \equiv (a+\sqrt{\omega})^p(a+\sqrt{\omega}) \\ & \equiv (a^p+\omega^{\frac{p}{2}})(a+\sqrt{\omega}) \\ & \equiv (a-\sqrt{\omega})(a+\sqrt{\omega}) \\ & \equiv a^2-\omega \equiv n \pmod p \end{align}$$
    由于$x^2 \equiv n \pmod p$有且仅有两根,且都不含$\sqrt{\omega}$项,而由算法给出的$x$是方程的一个根,所以它不含$\sqrt{\omega}$项。

代码

const int MOD = 998244353;
 
#define rnd (1ll * rand() * rand() % MOD)

struct Field {
    int a, b, w;
    Field(int _a = 0, int _b = 0, int _w = 0) : a(_a), b(_b), w(_w) {}
    friend Field operator*(const Field &x, const Field &y) {
        return Field((1ll * x.a * y.a + 1ll * x.b * y.b % MOD * x.w) % MOD, (1ll * x.a * y.b + 1ll * x.b * y.a) % MOD, x.w);
    }
} a[N];

int power(int a, int b) {
    int ret = 1;
    for (; b; b >>= 1) {
        if (b & 1) {
            ret = 1ll * ret * a % MOD;
        }
        a = 1ll * a * a % MOD;
    }
    return ret;
}
 
Field power(Field a, int b) {
    Field ret = Field(1, 0, a.w);
    for (; b; b >>= 1) {
        if (b & 1) {
            ret = ret * a;
        }
        a = a * a;
    }
    return ret;
}
 
int quadratic(int t) {
    int l = power(t, MOD >> 1);
    return l > 1 ? l - MOD : l;
}
 
int solve(int t) {
    int l = quadratic(t);
    if (l < 1) {
        return l;
    }
    int a = rnd, w = (MOD + 1ll * a * a - t) % MOD;
    for (; ~quadratic(w); a = rnd, w = (MOD + 1ll * a * a - t) % MOD) ;
    return power(Field(a, 1, w), MOD + 1 >> 1).a;
}

Four Loop

题目描述

Given two sequences $a$ and $b$ of equal length $n$, find the

$$\sum_{x=1}^n \sum_{y=1}^n \sum_{z=1}^n \sum_{w=1}^n (a_x+a_y+a_z+a_w)^{(b_x \oplus b_y \oplus b_z \oplus b_w)}$$

题意概述

给定两个长度为$n$的序列$a$和$b$,求$\sum_{x=1}^n \sum_{y=1}^n \sum_{z=1}^n \sum_{w=1}^n (a_x+a_y+a_z+a_w)^{(b_x \oplus b_y \oplus b_z \oplus b_w)}$。

数据范围:$1 \le n \le 10^5, \; 1 \le a_i \le 500, \; 1 \le b_i \le 500$。

算法分析

$a_x+a_y+a_z+a_w$的取值范围为$[4,2000]$,$b_x \oplus b_y \oplus b_z \oplus b_w$的取值范围为$[0,511]$,可以分别计算每种情况出现的次数再求和。

令$f_{1,i,j}$表示有多少个$x$满足$a_x=i \land b_x=j$。对于$k \ge 2$,令

$$f_{k,i,j}=\sum_{i_1+i_2=i} \sum_{j_1 \oplus j_2=j} f_{k-1,i_1,j_1}f_{1,i_2,j_2}$$

$f_{4,i,j}$即要求的每种情况出现的次数。转移方程在一维上是乘法卷积,另一维上是异或卷积,可以分别用FFT和FWT进行处理。总时间复杂度为$O(a_{\max}b_{\max}\log(a_{\max}b_{\max}))$。

代码

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
 
int const N = 100005, M = 512, MOD = 998244353, G = 3, INV2 = 499122177;
 
int a[N], b[N], rev[M << 2];
 
int power(int a, int b) {
    int ret = 1;
    for (; b; b >>= 1) {
        if (b & 1) {
            ret = 1ll * ret * a % MOD;
        }
        a = 1ll * a * a % MOD;
    }
    return ret;
}
 
int wn[M << 2], A[M << 2], f[M << 2][M];
 
void init(int n) {
    int len = 1, p = 0;
    for (; len < n; len <<= 1, ++p) ;
    for (int i = 1; i < len; ++i) {
        rev[i] = rev[i >> 1] >> 1 | (i & 1) << p - 1;
    }
    wn[0] = 1, wn[1] = power(G, (MOD - 1) / len);
    for (int i = 2; i < len >> 1; ++i) {
        wn[i] = 1ll * wn[i - 1] * wn[1] % MOD;
    }
}
 
void fft(int *a, int len, bool inv = 0) {
    for (int i = 0; i < len; ++i) {
        if (i < rev[i]) {
            std::swap(a[i], a[rev[i]]);
        }
    }
    for (int i = 1; i < len; i <<= 1) {
        for (int j = 0; j < len; j += i << 1) {
            for (int k = 0; k < i; ++k) {
                int x = a[j + k], y = 1ll * wn[len / (i << 1) * k] * a[j + i + k] % MOD;
                a[j + k] = (x + y) % MOD;
                a[j + i + k] = (MOD + x - y) % MOD;
            }
        }
    }
    if (inv) {
        std::reverse(a + 1, a + len);
        int inv = power(len, MOD - 2);
        for (int i = 0; i < len; ++i) {
            a[i] = 1ll * a[i] * inv % MOD;
        }
    }
}
 
void fwt(int *a, int len, bool inv = 0) {
    for (int i = 1; i < len; i <<= 1) {
        for (int j = 0; j < len; j += i << 1) {
            for (int k = 0; k < i; ++k) {
                int x = a[j + k], y = a[j + i + k];
                a[j + k] = (x + y) % MOD;
                a[j + i + k] = (MOD + x - y) % MOD;
                if (inv) {
                    a[j + k] = 1ll * a[j + k] * INV2 % MOD;
                    a[j + i + k] = 1ll * a[j + i + k] * INV2 % MOD;
                }
            }
        }
    }
}
 
int main() {
    int n;
    scanf("%d", &n);
    for (int i = 1; i <= n; ++i) {
        scanf("%d", &a[i]);
    }
    for (int i = 1; i <= n; ++i) {
        scanf("%d", &b[i]);
        f[a[i]][b[i]] = f[a[i]][b[i]] + 1;
    }
    init(M << 2);
    for (int i = 0; i < M << 2; ++i) {
        for (int j = 0; j < M; ++j) {
            A[j] = f[i][j];
        }
        fwt(A, M);
        for (int j = 0; j < M; ++j) {
            f[i][j] = A[j];
        }
    }
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < M << 2; ++j) {
            A[j] = f[j][i];
        }
        fft(A, M << 2);
        for (int j = 0; j < M << 2; ++j) {
            f[j][i] = A[j];
        }
    }
    for (int i = 0; i < M << 2; ++i) {
        for (int j = 0; j < M; ++j) {
            f[i][j] = power(f[i][j], 4);
        }
    }
    for (int i = 0; i < M << 2; ++i) {
        for (int j = 0; j < M; ++j) {
            A[j] = f[i][j];
        }
        fwt(A, M, 1);
        for (int j = 0; j < M; ++j) {
            f[i][j] = A[j];
        }
    }
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < M << 2; ++j) {
            A[j] = f[j][i];
        }
        fft(A, M << 2, 1);
        for (int j = 0; j < M << 2; ++j) {
            f[j][i] = A[j];
        }
    }
    int ans = 0;
    for (int i = 0; i < M << 2; ++i) {
        for (int j = 0; j < M; ++j) {
            if (f[i][j]) {
                ans = (ans + 1ll * f[i][j] * power(i, j)) % MOD;
            }
        }
    }
    printf("%d\n", ans);
    return 0;
}

Bitwise XOR

题目描述

Given a sequence $a$ and an empty sequence $b$, each element in $a$ is added to $b$ with probability $P$ (independently to the other elements). Denote $s=\bigoplus_{i=1}^{|b|} b_i$, where $\oplus$ denotes bitwise xor and $s=0$ if $b$ is empty. Find the expected value of $s^2$.

题意概述

给定一个长度为$n$的序列$a$,其中每个数都有$P$的概率被选中。令被选中的数的异或和为$s$。求$s^2$的期望。

数据范围:$1 \le n \le 10^5, \; 0 \le a_i \lt 10^9+7$。

算法分析

令$P(s)$表示异或和为$s$的概率。把$s$用二进制表示,最低位为第$0$位。

$$\begin{align} \sum_s s^2P(s) &= \sum_s (\sum_{i=0}^{29} [s二进制第i位为1] 2^i)^2P(s) \\ &= \sum_s P(s) \sum_{i=0}^{29} [s二进制第i位为1] 2^i \sum_{j=0}^{29} [s二进制第j位为1] 2^j \\ &= \sum_s P(s) \sum_{i=0}^{29} \sum_{j=0}^{29} [s二进制第i位和第j位都为1] 2^{i+j} \\ &= \sum_{i=0}^{29} \sum_{j=0}^{29} P(s二进制第i位和第j位都为1) 2^{i+j} \end{align}$$

先枚举$i$和$j$,令$f_{k,0/1,0/1}$表示从前$k$个数中选,异或和第$i$位、第$j$为分别为$0/1$、$0/1$的概率,对$f_{n,1,1}2^{i+j}$求和即为答案。

代码

#include <cstdio>
#include <cstring>
#include <algorithm>
 
int const N = 100005, MOD = 1000000007;
 
int a[N], f[N][2][2];
 
int power(int a, int b) {
    int ret = 1;
    for (; b; b >>= 1) {
        if (b & 1) {
            ret = 1ll * ret * a % MOD;
        }
        a = 1ll * a * a % MOD;
    }
    return ret;
}
 
int main() {
    int n, x, y;
    scanf("%d%d%d", &n, &x, &y);
    int p = 1ll * x * power(y, MOD - 2) % MOD, q = (MOD + 1 - p) % MOD;
    for (int i = 1; i <= n; ++i) {
        scanf("%d", &a[i]);
    }
    int ans = 0, ans2 = 0;
    f[0][0][0] = 1;
    for (int i = 0; i < 30; ++i) {
        for (int j = i; j < 30; ++j) {
            for (int k = 1; k <= n; ++k) {
                for (int _i = 0; _i < 2; ++_i) {
                    for (int _j = 0; _j < 2; ++_j) {
                        int __i = _i ^ (a[k] >> i & 1);
                        int __j = _j ^ (a[k] >> j & 1);
                        f[k][_i][_j] = (1ll * f[k - 1][_i][_j] * q + 1ll * f[k - 1][__i][__j] * p) % MOD;
                    }
                }
            }
            ans = (ans + (1ll << i + j + (i != j)) % MOD * f[n][1][1]) % MOD;
        }
    }
    printf("%d\n", ans);
    return 0;
}

Unbalanced Coin

题目描述

You have an unbalanced coin, the state of which (head or tail) after tossing depends on the previous state. If the head side is up now, the probability that the coin is still head up after tossing is $p$, and the probability that the tail side comes up is $1-p$. Similarly, if the tail side is up now, the probability that the coin is still tail up after tossing is $p$, and the probability that the head side comes up is $1-p$.

Assume the coin is head up initially. You will then toss it for $n$ times, during these $n$ times, every time the head side is up, you will get one score. What is the probability that you will get exactly $k$ score for $k=0,1,2,\cdots,n$?

题意概述

有一枚初始正面朝上的硬币。每次投掷时,硬币有$p$的概率保持不变,$1-p$的概率翻面。求$n$次投掷后,正面朝上的总次数分别为$0,1,2,\cdots,n$的概率。

数据范围:$1 \le n \le 10^5$。

算法分析

令$f_{i,0/1,j}$表示$i$次投掷后,硬币正面/反面朝上,且正面朝上的总次数为$j$概率。容易想到转移方程:

$$\begin{align} f_{i,0,j} &= f_{i-1,0,j-1} \times p+f_{i-1,1,j-1} \times (1-p) \\ f_{i,1,j} &= f_{i-1,0,j} \times (1-p)+f_{i-1,1,j} \times p \end{align}$$

这样直接转移的复杂度是$O(n^2)$,考虑如何优化。

令$f_{i,0/1}$的普通型生成函数为$F_{i,0/1}$,即$F_{i,0/1}=\sum_{j=0}^i f_{i,0/1,j}x^j$。则上述转移方程可表示为:

$$\begin{align} F_{i,0} &= F_{i-1,0} \times px+F_{i-1,1} \times (1-p)x \\ F_{i,1} &= F_{i-1,0} \times (1-p)+F_{i-1,1} \times p \end{align}$$

容易想到其矩阵形式:

$$\begin{equation} \left[\begin{matrix} px & (1-p)x \\ 1-p & p\end{matrix}\right] \left[\begin{matrix} F_{i-1,0} \\ F_{i-1,1} \end{matrix}\right]= \left[\begin{matrix} F_{i,0} \\ F_{i,1} \end{matrix}\right]\end{equation}$$

如此便可以用矩阵快速幂来优化了。其中多项式乘法需要用FFT计算。

代码

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
 
int const N = 200005;
double const PI = acos(-1);
 
int rev[N];
 
struct Complex {
    double a, b;
    Complex(double _a = 0, double _b = 0) : a(_a), b(_b) {}
    friend Complex operator+(Complex const &x, Complex const &y) {
        return Complex(x.a + y.a, x.b + y.b);
    }
    friend Complex operator-(Complex const &x, Complex const &y) {
        return Complex(x.a - y.a, x.b - y.b);
    }
    friend Complex operator*(Complex const &x, Complex const &y) {
        return Complex(x.a * y.a - x.b * y.b, x.a * y.b + x.b * y.a);
    }
    friend Complex operator/(Complex const &x, double const &y) {
        return Complex(x.a / y, x.b / y);
    }
} wn[N], A[2][2][N], B[2][2][N], C[2][2][N];
 
int init(int n) {
    int len = 1, p = 0;
    for (; len < n; len <<= 1, ++p) ;
    for (int i = 1; i < len; ++i) {
        rev[i] = rev[i >> 1] >> 1 | (i & 1) << p - 1;
    }
    for (int i = 0; i < len >> 1; ++i) {
        wn[i] = Complex(cos(2 * PI / len * i), sin(2 * PI / len * i));
    }
    return len;
}
 
void fft(Complex *a, int len, bool inv = 0) {
    for (int i = 0; i < len; ++i) {
        if (i < rev[i]) {
            std::swap(a[i], a[rev[i]]);
        }
    }
    for (int i = 1; i < len; i <<= 1) {
        for (int j = 0; j < len; j += i << 1) {
            for (int k = 0; k < i; ++k) {
                Complex x = a[j + k], y = wn[len / (i << 1) * k] * a[j + i + k];
                a[j + k] = x + y;
                a[j + i + k] = x - y;
            }
        }
    }
    if (inv) {
        std::reverse(a + 1, a + len);
        for (int i = 0; i < len; ++i) {
            a[i] = a[i] / len;
        }
    }
}
 
struct Matrix {
    int n[2][2];
    double *a[2][2];
 
    friend Matrix operator*(Matrix const &x, Matrix const &y) {
        int len = init(x.n[0][0] + y.n[0][0] - 1);
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                for (int k = 0; k < len; ++k) {
                    A[i][j][k] = k < x.n[i][j] ? x.a[i][j][k] : 0;
                    B[i][j][k] = k < y.n[i][j] ? y.a[i][j][k] : 0;
                    C[i][j][k] = 0;
                }
                fft(A[i][j], len);
                fft(B[i][j], len);
            }
        }
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                for (int k = 0; k < 2; ++k) {
                    for (int l = 0; l < len; ++l) {
                        C[i][j][l] = C[i][j][l] + A[i][k][l] * B[k][j][l];
                    }
                }
                fft(C[i][j], len, 1);
            }
        }
        Matrix ret;
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                ret.n[i][j] = 0;
                for (int k = 0; k < 2; ++k) {
                    ret.n[i][j] = std::max(ret.n[i][j], x.n[i][k] + y.n[k][j] - 1);
                }
                ret.a[i][j] = new double[ret.n[i][j]];
                for (int k = 0; k < ret.n[i][j]; ++k) {
                    ret.a[i][j][k] = C[i][j][k].a;
                }
            }
        }
        return ret;
    }
} base, ans;
 
int main() {
    int n;
    double p;
    scanf("%d%lf", &n, &p);
    base.n[0][0] = base.n[0][1] = 2;
    base.n[1][0] = base.n[1][1] = 1;
    base.a[0][0] = new double[2];
    base.a[0][0][0] = 0;
    base.a[0][0][1] = p;
    base.a[0][1] = new double[2];
    base.a[0][1][0] = 0;
    base.a[0][1][1] = 1 - p;
    base.a[1][0] = new double;
    base.a[1][0][0] = 1 - p;
    base.a[1][1] = new double;
    base.a[1][1][0] = p;
    ans.n[0][0] = ans.n[1][1] = 1;
    ans.a[0][0] = new double;
    ans.a[0][0][0] = 1;
    ans.a[1][1] = new double;
    ans.a[1][1][0] = 1;
    int b = n;
    for (; b;) {
        if (b & 1) {
            ans = ans * base;
        }
        if (b >>= 1) {
            base = base * base;
        }
    }
    for (int i = 0; i <= n; ++i) {
        printf("%.10f\n", ans.a[0][0][i] + ans.a[1][0][i]);
    }
    return 0;
}

Tsr and Variance

题目描述

Tsr is a cute boy with handsome moustache.

You are given a sequence with length $n$. Tsr wants you to calculate the sum of variance of each successive subsequence. Note: The variance in this problem should’t be divided by length.

Recall $\overline{a}_{l, r}=\frac{1}{r-l+1} \sum_{i=l}^r a_i$. Then you are supposed to calculate $\sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k-\overline{a}_{i, j})^2$.

题意概述

给定一个长度为$n$的序列,求$\sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k-\overline{a}_{i, j})^2$。有$T$组数据。

数据范围:$1 \le T \le 20, \; 1 \le n \le 10000, \; 1 \le a_i \le 10$。

算法分析

$$\begin{align} \sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k-\overline{a}_{i, j})^2 &= \sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k^2-2a_k \overline{a}_{i, j}+\overline{a}_{i, j}^2) \\ &= \sum_{k=1}^n \sum_{i=1}^k \sum_{j=k}^n a_k^2-2\sum_{i=1}^n \sum_{j=i}^n \overline{a}_{i, j} \sum_{k=i}^j a_k+ \sum_{i=1}^n \sum_{j=i}^n (j-i+1) \overline{a}_{i, j}^2 \\ &= \sum_{k=1}^n k(n-k+1)a_k^2-\sum_{i=1}^n \sum_{j=i}^n (j-i+1) \overline{a}_{i, j}^2 \end{align}$$

前一项可以在$O(n)$时间内求出,只需考虑如何求后一项。令$s_i=\sum_{j=1}^i a_j$。

$$\begin{align} \sum_{i=1}^n \sum_{j=i}^n (j-i+1) \overline{a}_{i, j}^2 &= \sum_{j=1}^n \sum_{i=1}^j \frac{(s_j-s_{i-1})^2}{j-i+1} \\ &= \sum_{j=1}^n \sum_{i=1}^j \frac{s_j^2-2s_js_{i-1}+s_{i-1}^2}{j-i+1} \\ &= \sum_{j=1}^n s_j^2 \sum_{i=1}^j \frac{1}{i}-2\sum_{j=1}^n s_j \sum_{i=1}^j \frac{s_{i-1}}{j-(i-1)}+\sum_{j=1}^n \sum_{i=1}^j \frac{s_{i-1}^2}{j-(i-1)} \end{align}$$

第一项也可以在$O(n)$时间内求出,而后面两项都是卷积的形式,可以用FFT在$O(n\log n)$时间内求出。

代码

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

int const N = 10005, T = 400005;
double const PI = acos(-1);

int a[N], rev[T];

struct Point {
    double x, y;
    Point(double _x = 0, double _y = 0) : x(_x), y(_y) {}
    friend Point const operator + (Point const &a, Point const &b) {
        return Point(a.x + b.x, a.y + b.y);
    }
    friend Point const operator - (Point const &a, Point const &b) {
        return Point(a.x - b.x, a.y - b.y);
    }
    friend Point const 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);
    }
    friend Point const operator / (Point const &a, double const &n) {
        return Point(a.x / n, a.y / n);
    }
} wn[T], A[T], B[T], C[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 main() {
    int T;
    scanf("%d", &T);
    for (; T--;) {
        int n;
        scanf("%d", &n);
        double ans = 0;
        for (int i = 1; i <= n; ++i) {
            scanf("%d", &a[i]);
            ans += 1. * a[i] * a[i] * i * (n - i + 1);
            a[i] += a[i - 1];
        }
        double sum = 0;
        for (int i = 1; i <= n; ++i) {
            sum += 1. / i;
            ans -= sum * a[i] * a[i];
        }
        int len = init(n << 1);
        for (int i = 0; i < n; ++i) {
            A[i] = Point(a[i]);
            B[i] = Point(1. * a[i] * a[i]);
            C[i] = Point(1. / (i + 1));
        }
        for (int i = n; i < len; ++i) {
            A[i] = B[i] = C[i] = Point(0);
        }
        fft(A, len);
        fft(B, len);
        fft(C, len);
        for (int i = 0; i < len; ++i) {
            A[i] = A[i] * C[i];
            B[i] = B[i] * C[i];
        }
        fft(A, len, 1);
        fft(B, len, 1);
        for (int i = 0; i < n; ++i) {
            ans += 2 * a[i + 1] * A[i].x - B[i].x;
        }
        printf("%.10f\n", ans);
    }
    return 0;
}

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

Black-White Balls

题目描述

$n$ black and white balls were put into a bag. Petya doesn’t know exactly how many black balls are there among them. He knows, however, that there are $0, 1, \ldots, n$ black balls among all balls in the bag with equal probability.
Petya took $l$ balls from the bag at random, and $l_1$ of them turned out black, while $l_2$ other turned out white ($l_1+l_2=l$). Now he wants to predict how many black balls were there initially in the bag. Of course, if $l \lt n$, he can’t be sure in his prediction, but he wants to predict a segment $[a, b]$, such that the amount $k$ of black balls belongs to it with probability at least $p$.
You are given $n, l_1, l_2$ and $p$, and you must find such $a$ and $b$ that $b-a$ is minimal possible. If there are several such pairs $(a, b)$, choose the one with the smallest $a$.

题意概述

袋子里有$n$个球,每个球可能是黑色或白色。已知黑球个数在$0, 1, \ldots, n$之间等概率分布。现在从中取出$l$个球,其中有$l_1$个黑球和$l_2$个白球($l_1+l_2=l$)。要求找到一个最短的区间$[a, b]$,使得黑球个数在$[a, b]$中的概率至少为${p \over 100}$。若有多个这样的区间,输出$a$最小的。
数据范围:$1 \le n \le 50, \; 0 \le l_1 \le n, \; 0 \le l_2 \le n-l_1, \; 0 \le p \le 100$。

算法分析

令事件$A$为取出的$l$个球中有$l_1$个黑球,事件$B_i$为袋子里有$i$个黑球。根据贝叶斯定理
$$P(B_i|A)={P(B_i)P(A|B_i) \over \sum_{j=0}^n P(B_j)P(A|B_j)}$$
而我们知道
$$P(B_i)={1 \over n+1}, \; P(A|B_i)={{i \choose l_1}{n-i \choose l_2} \over {n \choose l}}$$
因此可以计算出$P(B_i|A)$,接下来只要枚举区间即可。

代码

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

static int const N = 55;
double a[N];

double get_c(int n, int m) {
  if (n < m)
    return 0;
  double ret = 1;
  for (int i = 0; i < m; ++i)
    ret *= 1. * (n - i) / (m - i);
  return ret;
}

int main() {
  int n, l1, l2, p;
  scanf("%d%d%d%d", &n, &l1, &l2, &p);
  double b = 0;
  for (int i = 0; i <= n; ++i)
    b += a[i] = 1. / (n + 1) * get_c(i, l1) * get_c(n - i, l2) / get_c(n, l1 + l2);
  for (int i = 1; i <= n + 1; ++i)
    for (int j = 0; j + i - 1 <= n; ++j) {
      double sum = 0;
      for (int k = j; k <= j + i - 1; ++k)
        sum += a[k];
      if (sum / b * 100 + 1e-12 >= p)
        return printf("%d %d\n", j, j + i - 1), 0;
    }
  return 0;
}

Fermat’s Last Theorem

题目描述

Given a positive integer $n$ and a positive prime number $p$, find $x, y$ and $z$ such that $x^n+y^n \equiv z^n \pmod p$ and $x, y$ and $z$ are nonzero modulo $p$ or report that there’s no such triple.

题意概述

给定一个整数$n$和一个质数$p$,判断方程$x^n+y^n \equiv z^n \pmod p \; (1 \le x, y, z \lt p)$是否存在整数解,若存在则给出一组解。有$T$组数据。
数据范围:$1 \le T \le 1000, \; 3 \le n \le 10^6, \; 2 \le p \le 10^6$。

算法分析

为了方便,以下同余方程模数均为$p$。
求出$p$的原根$g$,设$x=g^{k_1}, \; y=g^{k_2}, \; z=g^{k_3}$,方程变为
$$g^{k_1n}+g^{k_2n} \equiv g^{k_3n}$$
令$G=g^n$,则
$$G^{k_1}+G^{k_2} \equiv G^{k_3}$$
由于$g$是原根,可知$G \not \equiv 0$,因此
$$
\begin{align}
G^{k_1-k_2}+1 &\equiv G^{k_3-k_2} \\
G^{k_1-k_3}+G^{k_2-k_3} &\equiv 1
\end{align}
$$
根据费马小定理,$G^{p-1} \equiv 1$,所以只要从$1$到$p-1$枚举$i$,判断是否存在$j \le i$满足
$$G^i+1 \equiv G^j \lor G^j+1 \equiv G^i \lor G^i+G^j \equiv 1$$
若存在,则找到了一组解。若枚举到$i$时发现存在$j \lt i$满足$G^i \equiv G^j$,说明开始循环,直接输出无解。根据抽屉原理,枚举的次数不会太多。
此题时限较小,不能每组询问都清空数组,可以使用时间戳。

代码

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

static int const N = 1000005;
int a[N], b[N], pri[N];

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

int get_prt(int p) {
  int top = 0, k = p - 1;
  for (int i = 2; i * i <= k; ++i)
    if (!(k % i)) {
      pri[top++] = i;
      for (; !(k % i); k /= i)
        ;
    }
  if (k > 1)
    pri[top++] = k;
  for (int i = 2; i <= p - 1; ++i) {
    bool fg = 1;
    for (int j = 0; fg && j < top; ++j)
      fg &= power(i, (p - 1) / pri[j], p) > 1;
    if (fg)
      return i;
  }
  return 0;
}

int main() {
  int T;
  scanf("%d", &T);
  for (int n, p; T--;) {
    scanf("%d%d", &n, &p);
    int g = get_prt(p), gn = power(g, n, p), x = 0, y = 0, z = 0;
    for (int i = 1, G = gn; i <= p - 1; ++i, G = 1ll * G * gn % p)
      if (b[G] == T + 1)
        break;
      else {
        a[G] = i, b[G] = T + 1;
        if (b[G + 1] == T + 1) {
          x = power(g, i, p), y = 1, z = power(g, a[G + 1], p);
          break;
        }
        if (b[G - 1] == T + 1) {
          x = power(g, a[G - 1], p), y = 1, z = power(g, i, p);
          break;
        }
        if (b[1 - G + p] == T + 1) {
          x = power(g, i, p), y = power(g, a[1 - G + p], p), z = 1;
          break;
        }
      }
    x ? printf("%d %d %d\n", x, y, z) : puts("-1");
  }
  return 0;
}

Divisibility

题目描述

Inspired by Stephen Graham, the King of Berland started to study algorithms on strings. He was working days and nights, having a feeling that the full potential in this area is still to be unlocked. And he was right!
One day, all the sudden, he made a huge breakthrough by discovering the fact that strings can be magically transformed into integer numbers. It was so simple! You just have to map different letters to different digits and be careful enough not to introduce any leading zeroes.
Here is what he wrote in his textbook about the string ‘lalala’:

  • it can be transformed to an $282828$ by mapping ‘l’ to $2$, and ‘a’ to $8$,
  • it can also be transformed to $909090$ by mapping ‘l’ to $9$, and ‘a’ to $0$,
  • acouple of examples of invalid transformations are $050505$ (the resulting number has a leading zero), $333333$ (different letters are mapped to the same digit), $123456$ (no mapping to the original letters at all).

But then things started to become more interesting. Obviously, it was known from very beginning that a single string can potentially be mapped to a variety of different integer numbers. But the King couldn’t even imagine that all numbers produced by the same string pattern might have common properties!
For example, every single number that can be produced from string ‘lalala’ is always divisible by $259$, irrespective of the letter-to-digit mapping you choose. Fascinating!
So the King ended up with the following problem. For any given string, he wanted to come up with an algorithm to calculate the set of its divisors. A number is called a divisor of the given string if all positive integers, that could possibly be produced from the given string, are divisible by it.
As usual, the King desperately wants you to help him, so stop thinking and start acting!

题意概述

给定一个字符集大小不超过$10$的字符串$s$。定义一个数字是合法的,当且仅当它没有前导零,位数等于$|s|$,且它从高到低的第$i$位和第$j$位相等当且仅当$s_i=s_j$。求出所有能整除每一个合法数字的数。有$n$组数据。
数据范围:$1 \le n \le 100, \; 1 \le |s| \le 14$。

算法分析

如果我们求出了所有合法数字的GCD,那么答案就是GCD的所有约数。
为了方便,我们设字符串的下标从$1$开始。
令$g$等于所有合法数字的GCD。首先考虑它的下界。假设字符$i$在字符串中出现的次数为$a_i$,出现的位置是$p_{i, 1}, p_{i, 2}, \ldots, p_{i, a_i}$。定义字符$i$的价值$v_i=\sum_{j=1}^{a_i} 10^{|s|-p_{i, j}}$。显然,所有合法的数字都可以被表示为$\sum k_iv_i$,其中$k_i$表示字符$i$变成的数字。所以,$g$的下界为$(v_a, v_b, \ldots, v_z)$。
考虑$g$的上界。分两种情况讨论。

  • 引理 $\forall a \le b, \; (a, b) \mid b-a$。

当字符集大小小于$10$时,对于每一个$i$,总能找到一个合法的数$x$使得$x+v_i$也合法。这是因为总有一个合法的数包含$k_i$但不包含$k_i+1$。根据引理,$(x, x+v_i) \mid v_i$。由于$x$和$x+v_i$均合法,因此$g \mid (x, x+v_i)$,即$g \mid v_i$。所以,$g$的上界为$(v_a, v_b, \ldots, v_z)$。结合下界,$g=(v_a, v_b, \ldots, v_z)$。
当字符集大小等于$10$时,要将某个$k_i$加一就必须将另一个$k_j$减一。如果对于不同的$i, j$,有$v_i, v_j \gt 0 \land v_i \lt v_j$,那么$g \mid v_j-v_i$。令$S=\{v_j-v_i \mid v_i, v_j \gt 0 \land v_i \lt v_j\}$。对于一个合法的数,它加上/减去$S$中的某些数后仍是一个合法的数;对于$S$中的数,存在一个合法的数加上/减去它之后仍是一个合法的数。即一个合法的数$x$与$S$中的数加减可以得到所有合法的数。令$h$等于$S$中所有数的GCD,那么显然,$g$的下界为$(x, h)$。又因为$g \mid x$且$g$整除$S$中的所有数,所以$g$的上界也是$(x, h)$。即$g=(x, h)$。
求$g$的所有约数可以先分解质因数再DFS。

代码

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

static int const L = 15;
static int const M = 10000005;
char s[L];
int top, cnt, cc, p[M];
long long v[26], dv[L], dc[L], ans[M];

long long get_gcd(long long a, long long b) {
  for (long long c; b; c = a % b, a = b, b = c)
    ;
  return a;
}

void init() {
  top = 0;
  for (int i = 2; i < M; ++i) {
    if (!p[i])
      p[top++] = i;
    for (int j = 0; j < top; ++j) {
      int k = i * p[j];
      if (k >= M)
        break;
      p[k] = 1;
      if (!(i % p[j]))
        break;
    }
  }
}

void dfs(int t, long long v) {
  if (t == cnt)
    return void(ans[cc++] = v);
  for (int i = 0; i <= dc[t]; ++i, v *= dv[t])
    dfs(t + 1, v);
}

int main() {
  int n;
  scanf("%d", &n), init();
  for (int _ = 1; _ <= n; ++_) {
    printf("Case %d:", _), scanf("%s", s);
    int len = strlen(s), st = 0;
    for (int c = 0; c < 26; ++c) {
      v[c] = 0;
      for (int i = 0; i < len; ++i) {
        v[c] *= 10;
        if (s[i] == c + 'a')
          ++v[c];
      }
      st += v[c] > 0;
    }
    long long gcd = 0;
    if (st < 10)
      for (int i = 0; i < 26; ++i)
        gcd = get_gcd(gcd, v[i]);
    else {
      int k = 0;
      for (int i = 0; i < 26; ++i)
        if (v[i])
          gcd += k++ * v[i];
      for (int i = 0; i < 26; ++i)
        for (int j = 0; j < 26; ++j)
          if (v[i] && v[j] && v[i] < v[j])
            gcd = get_gcd(gcd, v[j] - v[i]);
    }
    cnt = cc = 0;
    for (int i = 0; i < top && 1ll * p[i] * p[i] <= gcd; ++i)
      if (!(gcd % p[i])) {
        dv[cnt] = p[i], dc[cnt] = 0;
        for (; !(gcd % p[i]);)
          gcd /= p[i], ++dc[cnt];
        ++cnt;
      }
    if (gcd > 1)
      dv[cnt] = gcd, dc[cnt] = 1, ++cnt;
    dfs(0, 1), std::sort(ans, ans + cc);
    for (int i = 0; i < cc; ++i)
      printf(" %lld", ans[i]);
    puts("");
  }
  return 0;
}

Gena vs Petya

题目描述

Gena and Petya love playing the following game with each other. There are $n$ piles of stones, the $i$-th pile contains $a_i$ stones. The players move in turns, Gena moves first. A player moves by choosing any non-empty pile and taking an arbitrary positive number of stones from it. If the move is impossible (that is, all piles are empty), then the game finishes and the current player is considered a loser.
Gena and Petya are the world famous experts in unusual games. We will assume that they play optimally.
Recently Petya started to notice that Gena wins too often. Petya decided that the problem is the unjust rules as Gena always gets to move first! To even their chances, Petya decided to cheat and take and hide some stones before the game begins. Since Petya does not want Gena to suspect anything, he will take the same number of stones $x$ from each pile. This number $x$ can be an arbitrary non-negative integer, strictly less that the minimum of $a_i$ values.
Your task is to find the number of distinct numbers $x$ such that Petya will win the game.

题意概述

给定$n$个正整数$a_i$,求有多少个非负整数$x$,满足$x$小于所有给定的数,且所有给定的数减去$x$之后的异或和为$0$。
数据范围:$1 \le n \le 2 \times 10^5, \; 1 \le a_i \le 10^{18}$。

算法分析

Nim游戏后手胜利的条件是所有数的异或和为$0$。
枚举$x$显然不可行。可以尝试从低位到高位依次确定。
考虑低位对高位的影响,易知只有当低位退位时才会改变高位。而对于每一位,若此位上有偶数个$0$,则此位可以填$1$,若有偶数个$1$则可以填$0$(这样才能使得异或和为$0$)。
令$f_{i, j}$表示从低到高处理到第$i$位时,有$j$个数字在第$(i+1)$位退位的方案数。根据退位的性质,这$j$个数字一定是所有给定的数按后$i$位从小到大排序后的前$j$个。可以在枚举$i$的同时基数排序,枚举$j$的同时维护第$i$位上$0$和$1$的个数(会受退位影响),若满足条件则转移(注意此位填$1$时还要考虑此位$0$的退位)。
最后需要特判一下$x$不能等于所有给定的数中最小的数。

代码

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

#define int long long

static int const N = 200005;
static int const M = 65;
int a[N], cnt[M], f[M][N], s[2][N];

signed main() {
  int n;
  scanf("%lld", &n);
  for (int i = 0; i < n; ++i) {
    scanf("%lld", a + i);
    for (int j = 0; j < M; ++j)
      cnt[j] += a[i] >> j & 1;
  }
  f[0][0] = 1;
  for (int i = 0; i < M - 1; ++i) {
    int c0 = n - cnt[i], c1 = cnt[i], c = 0;
    for (int j = 0; j <= n; ++j) {
      if (j)
        if (a[j - 1] >> i & 1)
          ++c0, --c1;
        else
          --c0, ++c1, ++c;
      if (!(c0 & 1))
        f[i + 1][c + c0] += f[i][j];
      if (!(c1 & 1))
        f[i + 1][c] += f[i][j];
    }
    *s[0] = *s[1] = 0;
    for (int j = 0; j < n; ++j)
      s[a[j] >> i & 1][++*s[a[j] >> i & 1]] = a[j];
    for (int j = 1; j <= *s[0]; ++j)
      a[j - 1] = s[0][j];
    for (int j = 1; j <= *s[1]; ++j)
      a[*s[0] + j - 1] = s[1][j];
  }
  int sum = 0;
  for (int i = 1; i < n; ++i)
    sum ^= a[i] - a[0];
  printf("%lld\n", f[M - 1][0] - (sum == 0));
  return 0;
}