HDU5728 PowMod

11 May 2019 | duanyll | Tags: OI 题解 数论 欧拉函数

由于posts.json的BUG文章详情省略,不过这是道好题

题意

给定$n,m,p < 10^7$,保证$n$没有平方因子,定义

\[k=\sum_{i=1}^{m} \varphi (i*n)\ mod\ 1000000007\]

\[ans=k^{k^{k^{k^{\dots^k}}}}\ mod \ p\]

分析

k的计算

欧拉函数的定义:小于$n$的正整数中与$n$互质的数的数目,$\varphi(1)=1$。欧拉函数是积性函数,即满足:

\[\varphi(ab) = \varphi(a)\varphi(b),gcd(a,b)=1\]

以此可以推出它的通式(证明略):

\[\varphi(x) = x \prod_{i=1}^{n}(1-\frac{1}{1-p_i}),p_i\textrm{是x的质因数}\]

利用它的通式可以得到

\[\varphi(kpn) = p\varphi(kn),p\textrm{是x的质因数}\]

故考虑$\varphi(i*n)$时,如果$i$与$n$互质,就可以直接用积性函数的性质来求。如果$i$与$n$不互质:令$p$是$n$的一个质因子,由于$n$没有平方因子, 那么$gcd(n,n/p)=1$,如果$i$没有因子$p$,就用积性函数求;而如果$i$是$p$的倍数,就用上面那条性质把$p$拿出来。

注意当$p$是质数时$\varphi(p)=p-1$

\[\begin{align} \sum_{i=1}^{m} \varphi (i*n) &= \sum_{i=1}^{m/p}\varphi(i*p*n) + \sum_{i=1,i\%p>0}^{m}\varphi(i*\frac{n}{p})*\varphi(p) \nonumber\\ &= \sum_{i=1}^{m/p}\varphi(i*n)*p + \sum_{i=1,i\%p>0}^{m}\varphi(i*\frac{n}{p})*\varphi(p) \nonumber\\ &= \sum_{i=1}^{m/p}\varphi(i*n)*(\varphi(p)+1) + \sum_{i=1,i\%p>0}^{m}\varphi(i*\frac{n}{p})*\varphi(p) \nonumber\\ &= \sum_{i=1}^{m/p}\varphi(i*n) + \sum_{i=1}^{m/p}\varphi(i*n)*\varphi(p) + \sum_{i=1,i\%p>0}^{m}\varphi(i*\frac{n}{p})*\varphi(p) \nonumber\\ \end{align}\]

此时考虑后两项中$i$的意义:中间一项中$p$的倍数$p*i$对应最后一项中$i$恰好取完了1到m的所有值,故可以合并:

\[\sum_{i=1}^{m} \varphi (i*n) = \sum_{i=1}^{m/p}\varphi(i*n) + \varphi(p)\sum_{i=1}^{m}\varphi(i*\frac{n}{p})\]

此时如果把$m/p$与$n/p$看做整体,设$f(n,m) = \sum_{i=1}^{m} \varphi (i*n)$,有

\[f(n,m) = \left\{ \begin{array}{ll} \varphi(p)*f(n/p,m)+f(n,m/p)&n>1\\ sum[m] & n=1 \end{array} \right.\]

计算时每次枚举n的一个质因子。

计算k的无限次幂

欧拉定理:

\[a^{\varphi(n)} \equiv 1 (mod~n),gcd(a,n)=1\]

可以推出指数循环节公式:

\[a^b \equiv a^{b\%\varphi(c)+\varphi(c)} (mod~c),b \geq \varphi(c)\]

当c是质数时,$\varphi(p)=p-1$,根据费马小定理有

\[a^b \equiv a^{b\%(c-1)} (mod~c)\]

不过上面那个式子与本题无关,本题中那个无限次幂每向上计算一次就会取一个欧拉函数,最后模数会越取越小变成1,膜1结果是0就跳出死循环了。

代码

#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;

typedef long long int64;

const int MAXA = 1e7 + 10;
const int AMXA = 1e7;
const int64 MOD = 1000000007;

int prime[MAXA], prime_cnt;
bool isntp[MAXA];
int64 phi[MAXA], sum_phi[MAXA];

void init_phi() {
    phi[1] = 1;
    for (int i = 2; i <= AMXA; i++) {
        if (!isntp[i]) {
            prime[++prime_cnt] = i;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= prime_cnt && i * prime[j] <= AMXA; j++) {
            isntp[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            } else {
                phi[i * prime[j]] = phi[i] * (prime[j] - 1);
            }
        }
    }
}

template <typename T>
T pow_mod(T a, T b, T MOD) {
    T res = 1;
    while (b) {
        if (b & 1) res = res * a % MOD;
        a = a * a % MOD;
        b /= 2;
    }
    return res;
}

int factors[50], fac_cnt;

int64 f(int64 n, int64 m, int facid) {
    if (n == 1) {
        return sum_phi[m];
    }
    if (m == 0) {
        return 0;
    }
    return ((factors[facid] - 1) * f(n / factors[facid], m, facid - 1) % MOD +
            f(n, m / factors[facid], facid) % MOD) %
           MOD;
}

int64 get_pow(int64 k, int64 p) {
    if (p == 1) {
        return 1;
    }
    int64 tmp = pow_mod(k, get_pow(k, phi[p]), p);
    return (tmp == 0) ? p : tmp;
}

int main() {
    init_phi();
    for (int i = 1; i <= AMXA; i++) {
        sum_phi[i] = sum_phi[i - 1] + phi[i];
        sum_phi[i] %= MOD;
    }
    int64 n, m, p;
    while (cin >> n >> m >> p) {
        // fac_cnt = 0;
        // if (isntp[n]) {
        //     for (int i = 1; i <= prime_cnt && prime[i] * prime[i] <= n; i++)
        //     {
        //         if (n % prime[i] == 0) {
        //             factors[++fac_cnt] = prime[i];
        //             if (!isntp[n / prime[i]]) {
        //                 factors[++fac_cnt] = n / prime[i];
        //             }
        //         }
        //     }
        // } else {
        //     factors[++fac_cnt] = n;
        // }
        fac_cnt = 0;
        int64 tmp = n;
        for (int i = 1; i <= prime_cnt; i++) {
            if (!isntp[tmp]) {
                factors[++fac_cnt] = tmp;
                break;
            } else {
                if (tmp % prime[i] == 0) {
                    factors[++fac_cnt] = prime[i];
                    tmp /= prime[i];
                }
            }
        }
        int64 k = f(n, m, fac_cnt);
        cout << get_pow(k, p) << endl;
    }
}