hdu5728PowMod

题目：

给定 n,m,p

先得到   k=∑mi=1φ(i∗n) mod 1000000007

k=∑mi=1φ(i∗n) mod 1000000007

其中n为非平方数

再计算ans=kkkk...k mod p

这里有无穷个k

思路：

1。求k

欧拉函数是非完全积性函数，φ(a*b)=φ(a)*φ(b)，当gcd(a,b)=1;

φ(i*n)=φ(a)*φ(b)     a*b=i*n     gcd(a,b)=1;

遍历出一个n的约数

显然第一次遍历到的是一个素数p，并且是n的最小素因子

利用   φ(a*b)=φ(a)*φ(b)，当gcd(a,b)=1; p和（n/p） 互素        ∑mi=1φ(i∗n)= φ(p）* ∑mi=1φ(i∗n/p)

似乎不太对， p和（n/p） 互素  ，但是 p和i不一定互素

当p和i不互素时 ，此时的i=kp，并且p是一个质数

φ(i∗n)=p*φ(i∗n/p)

前面的φ(p）* ∑mi=1φ(i∗n/p)已经 加上了 φ(p）个 [即p-1个]φ(i∗n/p)

也即  ∑mi=1φ(i∗n)= φ(p）* ∑mi=1φ(i∗n/p) +∑  φ(i∗n/p)

在∑  φ(i∗n/p)中，i=kp  ，∑  φ(i∗n/p)=∑  φ(k∗n)    ,k=1,2,.....,m/p

即最后

∑(i=1,m)φ(i*n) = φ(pi) * ∑(i=1,m)φ(i*n/pi) + ∑(i=1,m/pi)φ(i*n) ;

这样k便可以递归的求出来了

2。求ans

我们知道欧拉定理得到的 指数循环节： A^x = A^(x % φ(C) + φ(C)) (mod C)  (x >= φ(C))

ans=k^k....(mod c) = k^(k^k... % φ(C) + φ(C)) (mod C)  (k^k...... >= φ(C))

而其中的k^k^k^k... % φ(C)

又可以为 ans'=k^k..... (modφ(C) )= k^(k^k... %φ(φ(C)) + φ(φ(C))) (modφ(C))  (k^k...... >= φ(C))

往下递归下去

φφφφφφφφφφφ （C）会收敛到1

任何数mod1都为 0可以结束递归返回快速幂的答案

#include<bits/stdc++.h>
using namespace std;
//const int M =10000010,N=1000000007 ;
const int M = 1e7 + 5 ;
const int N = 1e9 + 7 ;
int prime[M],euler[M],s[M];
int res;	
int n,m,p;
void init()
{
	res=0;
	euler[1]=1;
	for(int i=2;i<M;i++)
	{
		if(!euler[i])
		{
		 prime[res++]=i;
		 euler[i]=i-1;	
		}
		for(int j=0;j<res&&i*prime[j]<M;j++)
		{
			if(i%prime[j]==0)
            {
                euler[prime[j]*i]=euler[i]*prime[j];
                break;
            }
            euler[prime[j]*i]=euler[i]*(prime[j]-1);
		}
		
	}
	for(int i=1;i<M;i++)
	s[i]=(euler[i]+s[i-1])%N;
}

long long qu(long long a,long long b,long long p)
{
	long long  ans=1;
	a=a%p;
	while(b)
	{
		if(b&1) ans=(ans*a)%p;
		b>>=1;
		a=(a*a)%p;
	}
	return ans%p;
}
long long f(long long k,long long p)
{   
	if(p==2) return 0;
	return qu(k,f(k,euler[p])+euler[p],p);
}
long long ans(long long n,long long m)
{
    if( m < 1 ) return 0 ;
    if( m == 1 ) return euler[n];
    if( n == 1 ) return s[m] ;
    if( euler[n] == n-1 )  return ( ans(1,m) * (n-1) % N+ ans(n,m/n) ) % N ;
     for( int i = 2 ; i*i <= n ; i++ )
    {
        if( n%i ) continue ;
        return ((i-1)* ans(n/i,m) % N+ ans(n,m/i) ) % N ;
    } 
}
int main()
{
	init();
	while(~scanf("%d%d%d",&n,&m,&p))
	 cout<< f(ans(n,m),p)<<endl;	
	return 0;
}