[洛谷P2158]仪仗队

难度：提高+/省选-

题目描述

        作为体育委员，C君负责这次运动会仪仗队的训练。仪仗队是由学生组成的N * N的方阵，为了保证队伍在行进中整齐划一，C君会跟在仪仗队的左后方，根据其视线所及的学生人数来判断队伍是否整齐(如下图)。 现在，C君希望你告诉他队伍整齐时能看到的学生人数。

输入输出格式

输入格式:

        共一个数N(1≤N≤40000)。

输出格式:

        共一个数，即C君应看到的学生人数。

输入输出样例

Sample input

4

Sample output

9

题解

        以左下方的观察点为原点建系，发现可以观察到的人总在正上方和正右方以及右上方一个n-1阶的方阵中。
        以图中n=6为例，考察右上方5阶方阵中的观察情况。
        (1，5)        (2，5)        (3，5)        (4，5)        (5，5)
        (1，4)        (2，4)        (3，4)        (4，4)        (5，4)
        (1，3)        (2，3)        (3，3)        (4，3)        (5，3)
        (1，2)        (2，2)        (3，2)        (4，2)        (5，2)
        (1，1)        (2，1)        (3，1)        (4，1)        (5，1)
        其中红色标出的是观察不到的人的坐标。发现可以被观察到等价于其横纵坐标互素。
        那么若令f(x)表示n阶方阵中横纵坐标互素的点的个数，则答案即为f(N)+2。注意N=1时特判为0。
        下面探讨f(x)的求法，很自然想到递推。
        (1，5)        (2，5)        (3，5)        (4，5)        (5，5)
        (1，4)        (2，4)        (3，4)        (4，4)        (5，4)
        (1，3)        (2，3)        (3，3)        (4，3)        (5，3)
        (1，2)        (2，2)        (3，2)        (4，2)        (5，2)
        (1，1)        (2，1)        (3，1)        (4，1)        (5，1)
        如上图，在n=4(绿色部分)的基础上递推到n=5的情况，容易发现(5，5)一定会被剔除，然后(5，5)的正左方和正下方一定是成对剔除的。显然，如果令φ(x)表示小于x且与x互素的正整数的个数，则有以下递推式:

$$f(x)=f(x-1)+2*φ(x)$$

        现在问题的关键在于求解φ(x)。
        φ(x)在数论上是有定义的，称为欧拉函数(Euler’s totient function)。现在介绍欧拉函数的线性筛法。
        为了理解这个筛法，要了解以下三个基本原理:

• 若p为素数，则φ(p)=p-1
• 若p为素数且p†k，则φ(kp)=(p-1)φ(k)
• 若p为素数且p|k，则φ(kp)=pφ(k)

        证明见题后注解。
        用这三个原理结合欧拉筛素数法，可得到下面筛法:

int euler[40001] = {0};//记录欧拉函数值
int mark[40001] = {0};//记录是不是素数
int prim[40001], tot = 0;//储存当前已知的素数,tot 记录个数
for (int i = 2; i <= 40000; i++) {
    if (!mark[i]) {
    	prim[++tot] = i;
    	euler[i] = i - 1;//判断为素数,直接给欧拉函数赋值
    }
    for (int j = 1; j <= tot; j++) {
    	if (i * prim[j] > 40000)break;
    	mark[i * prim[j]] = 1;//标记这个数一定不是素数
    	if (i % prim[j] == 0) {
                euler[i * prim[j]] = euler[i] * prim[j];
                break;
            }
    	else euler[i * prim[j]] = euler[i] * (prim[j] - 1);
    }
}

        这个筛法的可行性除了与欧拉筛法有关外，还利用了一个事实:每一个合数都可以写成k*p的形式并且p<=k。prim数组中存放了所有不大于i的质数，由上面的原理可知在计算φ(i*prim[j])时φ(i)一定已经被求出了。
        这个筛法很重要，它是很多数论方面题目的重要算法。
        利用这个算法再结合递推式，本题目迎刃而解。

#include<iostream>

using namespace std;
int euler[40001] = {0};
int mark[40001] = {0};
int prim[40001], tot = 0;
int A[40001] = {0};

int main() {
    for (int i = 2; i <= 40000; i++) {
        if (!mark[i]) {
            prim[++tot] = i;
            euler[i] = i - 1;
        }
        for (int j = 1; j <= tot; j++) {
            if (i * prim[j] > 40000)break;
            mark[i * prim[j]] = 1;
            if (i % prim[j] == 0) {
                euler[i * prim[j]] = euler[i] * prim[j];
                break;
            } else euler[i * prim[j]] = euler[i] * (prim[j] - 1);
        }
    }
    A[1] = 1;
    for (int i = 2; i <= 40000; i++)A[i] = A[i - 1] + 2 * euler[i];
    int n;
    cin >> n;
    if (n == 1)cout << 0;
    else cout << A[n - 1] + 2;
    return 0;
}

注解

        这里给出欧拉函数性质后两个的证明。
        首先欧拉函数是积性函数，即在(a，b)=1时有φ(ab)=φ(a)φ(b)。
        那么在p†k时，由于p是素数，必有(p，k)=1，故φ(pk)=φ(p)φ(k)=(p-1)φ(k)。在p|k时，不妨设k=qp^s(p†q)，那么φ(pk)=φ(p^s+1q)=φ(p^s+1)φ(q)。这里φ(p^s+1)为p^s+1-p^s，即p^s(p-1)，因此φ(pk)=(p-1)p^sφ(q)。并且φ(k)=φ(q)φ(p^s)=φ(q)p^s-1(p-1)。代入即得φ(pk)=pφ(k)。

【扩展】欧拉定理

        提到欧拉函数不得不提欧拉定理，它的内容如下：

【欧拉定理】若a与m互质，则有$a^{\phi(m)}\equiv 1(mod m)$。
【推论】若a与m互质，则有$a^b\equiv a^{b\%\phi(m)}(mod m)$。

        这个定理可以用来实现超级快速幂（指数极大，需要高精度）,在a与m不互质的情况下，有欧拉函数扩展形式：

【扩展欧拉定理】对于a、b与m，在$b\geq \phi(m)$时，有$a^b\equiv a^{b\% \phi(m)+\phi(m)}(mod m)$。

        这样在b巨大，但m很小时，可以将b减小到m的欧拉函数级别。
        另一个扩展：直接求欧拉函数。由上文可知，欧拉函数可以通过线性筛筛出来，但是对于很大的数去求欧拉函数，并且只求一个的情况，线性筛比较低效。这里可以直接用欧拉函数公式在$O(\sqrt n)$下计算出。
        如果m的所有互异质因子为$p_1,p_2,\cdots,p_s$，则有：

$$\phi(m)=m\prod_{i=1}^s\frac {p_i-1} {p_i}$$

        可以通过下面的代码计算：

inline int getEuler(int x) {
    euler = 1;
    for (int i = 2; 1ll * i * i <= x; i++) {
        if (x % i == 0) {
            x /= i;
            euler *= (i - 1);
            while (x % i == 0)x /= i, euler *= i;
        }
    }
    if (x > 1)euler *= x - 1;
    return euler;
};

【扩展】线性筛

        线性筛可以在线性时间内求解积性函数，这里的欧拉函数就是积性函数。利用线性筛，我们不仅可以筛出欧拉函数还可以筛出其它积性函数。

约数个数

        规定$d(x)$表示x的约数个数，比如$d(6)=4$，是因为6含有1、2、3、6四个因子，$d(x)$是积性函数。1~n的约数个数可以在线性时间内筛出。为了保证能够筛出，这里需要引入一个辅助数组$num[x]$，其表示x除去所有最小质因子后的数，这样可以得到下面筛法代码：

#include<bits/stdc++.h>

#define N 1000000
using namespace std;
int d[N + 1], prim[N + 1], tot, num[N + 1];
bool vis[N + 1];

int main() {
    d[1] = 1;
    for (int i = 2; i <= N; i++) {
        if (!vis[i])prim[++tot] = i, d[i] = 2, num[i] = 1;
        for (int j = 1; j <= tot; j++) {
            if (1ll * i * prim[j] > N)break;
            vis[i * prim[j]] = true;
            if (i % prim[j])d[i * prim[j]] = d[i] * d[prim[j]], num[i * prim[j]] = i;//积性函数
            else {
                d[i * prim[j]] = d[num[i * prim[j]] = num[i]] + d[i];//这一步是关键
                break;
            }
        }
    }
    for (int i = 1; i <= 100; i++)cout << d[i] << endl;//只输出前100验证
    return 0;
}

        显然在$i\%prim[j]\not =0$时，两者互质，直接利用积性函数性质即可。在两者不互质时，利用$num$数组求解。期间需要更新$num$数组。
        这里不易理解的是$d[i * prim[j]] = d[num[i * prim[j]] = num[i]] + d[i]$这一句。这里$i * prim[j]$的最小质因子显然为$prim[j]$，故其因子组成应该来自两部分，一部分是$i$本身就有的，另一部分是$i$除去最小质因子后的数（即$num[i]$）的所有因子与$i * prim[j]$全部最小质因子的乘积。从这个角度不难理解算法正确性。

约数和

        规定$d(x)$为x的所有因子之和，比如$d(6)=12$。可以证明$d(x)$也为积性函数，可以线性筛出。
        由约数个数中$num$数组的启发，我们应该保留$num$数组，但这样仍然求不出来，故需要引入另一个数组$minn[x]$，表示x全部最小质因子的乘积。和约数个数类似，代码改一下就好：

#include<bits/stdc++.h>

#define N 10000
using namespace std;
int d[N + 1], prim[N + 1], tot, num[N + 1], minn[N + 1];
bool vis[N + 1];

int main() {
    d[1] = 1;
    for (int i = 2; i <= N; i++) {
        if (!vis[i])prim[++tot] = i, d[i] = i + 1, num[i] = 1, minn[i] = i;
        for (int j = 1; j <= tot; j++) {
            if (1ll * i * prim[j] > N)break;
            vis[i * prim[j]] = true;
            if (i % prim[j])d[i * prim[j]] = d[i] * d[prim[j]], num[i * prim[j]] = i, minn[i * prim[j]] = prim[j];
            else {
                d[i * prim[j]] = d[i] + d[num[i * prim[j]] = num[i]] * (minn[i * prim[j]] = minn[i] * prim[j]);
                break;
            }
        }
    }
    for (int i = 1; i <= 100; i++)cout << d[i] << endl;
    return 0;
}