hdu 5728 （公式推导+指数循环节）

PowMod

Time Limit: 3000/1500 MS (Java/Others) Memory Limit: 262144/262144 K (Java/Others)
Total Submission(s): 1221 Accepted Submission(s): 428

Problem Description
Declare:
k=∑（m，i=1）φ(i∗n) mod 1000000007

n is a square-free number.

φ is the Euler’s totient function.

find:
ans=k^(k^(k^(k…k))) mod p

There are infinite number of k

Input
Multiple test cases(test cases ≤100), one line per case.

Each line contains three integers, n,m and p.

1≤n,m,p≤107

Output
For each case, output a single line with one integer, ans.

Sample Input
1 2 6
1 100 9

Sample Output
4
7

这题比较难求的是k的值，只要知道k的值这题就是一个指数循环节定理的裸题了。
观察这个式子显然直接求k是不行的，发现n的因子中是没有重复的数的，这意味着当p整除n且i%p！=0时 phi(n/p*i*p)可以写成phi(p)*phi(i*n/p) 而当i%p==0时，也就是说i含有质因子p这个式子可以写成
phi(i*n)=phi((i* (n/p)*p)=p*phi((i/p)*n) 令i/p =j 所以可以写成p*phi(j*n)
有了这两个式子，则原式= ∑（1 <= i <= m && i%p!=0）phi(p)φ(i∗n/p) + ∑（1 <= j <= m /p && j%p ==0）p*phi(j*n) 发现p=phi(p)+1 于是两个式子能和并在一起即 ∑（1 <= i <= m ）phi(p)φ(i∗n/p) +∑（1 <= i <= m / p ）φ(i∗n) 就可以递归的算出k的值了，后面一部分则是应用指数循环节定理A^B%C = A^(B%phi(C)+phi(C));所以依旧可以递归的求出。

#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <string>
#include <string.h>
#include <set>
#include <map>
#include <queue>
#include <math.h>
#include <stack>
#include <deque>
#include <map>
using namespace std;
const int maxn = 1e7+100;
const int mod = 1e9+7;
long long  cnt,n,m,pp;
long long phi[maxn],p[maxn],isp[maxn],rec[maxn],sum[maxn];
void init()  //欧拉筛法 O(n)
{
    int n=10000005;
    memset(isp,false,sizeof(isp));
    phi[1]=1;
    for(int i=2;i<n;i++)
    {
        if(!isp[i])
        {
            p[isp[0]++]=i;
            phi[i]=i-1;
        }
        for(int j=0;j<isp[0] && i*p[j]<n;j++)
        {
            isp[i*p[j]]=true;
            if(i%p[j]==0)
            {
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
            else phi[i*p[j]]=phi[i]*(p[j]-1);
        }
    }
    sum[0]=0;
    for(int i=1;i<maxn;i++)
        sum[i]=(sum[i-1]+phi[i])%mod;
}
long long  fk(long long  m,long long  n,int deep)
{
    if (m==0) return 0;
    if (n==1) return sum[m];
    return (phi[rec[deep]]*fk(m,n/rec[deep],deep+1)%mod+fk(m/rec[deep],n,deep))%mod;
}
long long f(long long  x,long long k,long long md)
{
    long long res=1;
    while(k){
        if(k&1) res=res*x%md;
        k>>=1;
        x=x*x%md;
    }
    return res;
}
long long solve(long long  k,long long md)
{
    if (md==1) return 1;
    return f((long long)k,solve(k,phi[md])+phi[md],md);
}
int main(){
     init();
     while (scanf("%I64d%I64d%I64d",&n,&m,&pp)!=EOF)
     {
         cnt=0;
         long long  tmp=n;
         for (int i=0;i<isp[0];i++)
         {
             if (tmp%p[i]==0)
             {
                 rec[cnt++]=p[i];
                 tmp/=p[i];
             }
             if (tmp==1) break;
         }
         if (tmp!=1) rec[cnt++]=tmp;
         long long  k=fk(m,n,0)%mod;
         printf("%I64d\n",solve(k,pp)%pp);
     }
}