hdu 2879 HeHe 数论 积性函数 优化技巧

题目

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=2879

题目来源：学弟问的。

简要题意：$f(n)$为$x^2\equiv x \mod n\quad (x\in [0,n))$的解的个数，
$\qquad \quad \ \ \ \ $求$\prod\limits_{i=1}^N{f(i)}\mod M$

数据范围：$1\leqslant N\leqslant 10^7;\quad 0<M\leqslant 10^9$

题解

解这题很好的一个方法是暴力打表然后观察，然后会发现$f(n)$和$n$的素因数个数有关。

$$\begin{aligned} x^2\equiv x \mod n \Leftrightarrow &\ x(x-1)\equiv 0 \mod n\\ \Leftrightarrow &\ n\mid x(x-1)\\ \Leftrightarrow &\ \gcd(n,x(x-1))=n \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \end{aligned}$$

显然$f(1)=1$

对于素数$p$来说$f(p)=2$因为根据素数定义$\forall x\in [1,p),\gcd(x,p)=1$则

$$\gcd(p,x(x-1))=\left\{ \begin{aligned} &\gcd(p,x-1)\cdot\gcd(p,x)=1~~~~~~~~~~~~x>1\\ &p~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~x=0,1 \end{aligned} \right.\qquad \qquad \qquad \qquad \qquad \qquad$$ 对于 $p^a$即素数的幂次 $f(p^a)=2$，因为 $\gcd(p^a,x)$是 $p^b(b<a)$而 $\gcd(x,x-1)=1$

再去考虑合数的情况，此处要使用积性函数的性质。

考虑互素的两个数$p,q$，则$f(pq)=f(p)\cdot f(q)$。

这里我无法给出证明，好难啊，想了两天了还是没成功，网上的各种证明没有让我信服的。

这个 属于稍微有点道理的，但是其实两个质数的情况不能直接推出两个互质的数的情况。

假设积性于是$f(n)=\prod\limits_{i=1}^kf(p_i^{a_i})=2^k$

设$n$的素因数个数为$\lambda(n)$则$f(n) = 2^{\lambda(n)}$

实现

实际写代码的话用结论直接去做肯定是会超时的。

需要转化为求$n$内所有素数的倍数的个数之和，设$p_1,p_2\cdots p_k$是$N$内的所有质数。

$g(N)=\sum\limits_{i=1}^k{\lfloor N/p_i\rfloor}$
$\prod\limits_{i=1}^N{f(i)}=2^{g(N)}$

这个玩意可以用很多方法去求，我为了改进效率试了很多种。

预处理素数+分块求和+bitset压缩内存

不分块，不用bitset也能过，但是这么搞之后可以提速，减少内存消耗。

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <cstring>
#include <stack>
#include <queue>
#include <string>
#include <vector>
#include <set>
#include <map>
#include <bitset>

#define pb push_back
#define mp make_pair
#define all(x) (x).begin(),(x).end()
#define sz(x) ((int)(x).size())
#define fi first
#define se second
using namespace std;
typedef long long LL;
typedef vector<int> VI;
typedef pair<int,int> PII;
inline LL powmod(LL a,LL b, LL MOD) {LL res=1;a%=MOD;for(;b;b>>=1){if(b&1)res=res*a%MOD;a=a*a%MOD;}return res;}
// head
LL MOD;
const int N = 10000000+500;
const int M = 670000;

bitset<N> vis;
int p[M];
int tot = 0;

void getPrime(int n) {
    for (int i = 2; i < n; i++) {
        if (!vis[i]) p[++tot] = i;
        for (int j = 1; j<=tot && i*p[j]<n; j++) {
            vis[i*p[j]] = true;
            if (i%p[j]==0) break;
        }
    }
}

int solve(int n) {
    int en, val, pw = 0, be = 1;
    while (p[be] <= n) {
        val = n/p[be];
        en = upper_bound(p+1, p+tot, n/val)-p;
        pw += val*(en-be);
        be = en;
    }
    return powmod(2, pw, MOD);
}

int main() {
    getPrime(N);
    int t, n;
    scanf("%d", &t);
    while (t--) {
        scanf("%d%I64d", &n, &MOD);
        printf("%d\n", solve(n));
    }
    return 0;
}

埃氏筛法预处理答案

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <cstring>
#include <stack>
#include <queue>
#include <string>
#include <vector>
#include <set>
#include <map>

#define pb push_back
#define mp make_pair
#define all(x) (x).begin(),(x).end()
#define sz(x) ((int)(x).size())
#define fi first
#define se second
using namespace std;
typedef long long LL;
typedef vector<int> VI;
typedef pair<int,int> PII;
LL powmod(LL a,LL b, LL MOD) {LL res=1;a%=MOD;for(;b;b>>=1){if(b&1)res=res*a%MOD;a=a*a%MOD;}return res;}
// head
LL MOD;
const int N = 10000000+5;
int res[N];
int tot = 0;
void getPrime(int n) {
    for (int i = 2; i <= n; i++) {
        if (!res[i]) {
            for (int j = i; j<=n; j += i) {
                res[j]++;
            }
        }
        res[i] += res[i-1];
    }
}

int main() {
    getPrime(1e7);
    int t, n;
    scanf("%d", &t);
    while (t--) {
        scanf("%d%I64d", &n, &MOD);
        printf("%I64d\n", powmod(2, res[n], MOD));
    }
    return 0;
}

欧拉筛预处理答案

void getPrime(int n) {
    int temp;
    for (int i = 2; i <= n; i++) {
        if (!res[i]) p[++tot] = i, res[i]++;
        for (int j = 1; j<=tot && (temp=i*p[j])<=n; j++) {
            (res[temp] = res[i])++;
            if (i%p[j]==0) {
                res[temp]--;
                break;
            }
        }
        res[i] += res[i-1];
    }
}

死扣加快好不容易到156ms再也上不去

inline void getPrime(int n) {
    int temp, l = n/2;
    for (int i = 2; i <= l; i++) {
        if (!res[i]) p[++tot] = i, res[i]++;
        for (int j = 1; j<=tot && (temp=i*p[j])<=n; j++) {
            (res[temp] = res[i])++;
            if (i%p[j]==0) {
                res[temp]--;
                break;
            }
        }
        res[i] += res[i-1];
    }
    for (int i = l+1; i <= n; i++) {
        if (!res[i]) res[i]++;
        res[i] += res[i-1];
    }
}