HDU 3923 polya+乘法逆元+扩展欧几里德

#include<iostream>
#include <stdlib.h>
#include <algorithm>
#include <stdio.h>
#include<vector>
using namespace std;
int prime[8084];
int isprime[10001];
int num,n,p;
const int Mod=1000000007;
void getprime()
{
num=0;
for(int i=2;i<=10000;i++)if(!isprime[i])
{
prime[num++]=i;
for(int j=1;j*i<=10000;j++)
{
isprime[i*j]=1;
}
}
}
long long euler(int x)
{
long long res=x;
for(int i=0;prime[i]*prime[i]<=x&&i<num;i++)
{
if(x%prime[i]==0)
{
res=res/prime[i]*(prime[i]-1);
while(x%prime[i]==0)
{
x/=prime[i];
}
}
}
if(x>1) res=res/x*(x-1);
return res;
}
long long fun(int a,int b)
{
long long ret=1;
long long aa=a;
while(b>0)
{
if(b&1)ret=(ret*aa)%1000000007;
aa=(aa)*(aa)%1000000007;
b>>=1;
}
return ret%1000000007;
}
__int64 Extend_euclid(__int64 a,__int64 b,__int64 &x,__int64 &y)
{
__int64 d=0,t=0;
if(b==0)
{
x=1;
y=0;
return a;
}
else
{
d=Extend_euclid(b,a%b,x,y);
t=x;
x=y;
y=t-a/b*y;
}
return d;
}
__int64 Bignum_Div(__int64 a,__int64 b,__int64 mod)
{
__int64 x=0,y=0;
Extend_euclid(b,mod,x,y);//扩展gcd求出(2*n)*x-mod*k=1的x
__int64 ans= a*x%mod;//ans/(2*n)%p = ((ans*x)/(2*n*x))%p=ans*x%mod; (p=1000000007,显然gcd(2*n,p)=1)
while(ans<0)
ans+=mod;
return ans;
}

int main()
{
getprime();
int c,Case;
scanf("%d",&Case);
for(int k=1;k<=Case;k++)
{
scanf("%d%d",&c,&n);
long long ans = 0;
for (int i = 1; i <= n; i++)if (n % i == 0)
{
ans = (ans+fun(c, i)* euler(n / i))%1000000007;
}
if (n & 1) ans =(ans+ n*(fun(c, n / 2 + 1)))%1000000007;
else ans = (ans+n/2*(fun(c, n / 2)+ fun(c, n / 2 + 1)))%1000000007;
ans=Bignum_Div(ans,2*n,1000000007);//不能直接ans/(2*n)，因为ans是mod%1000000007的结果
printf("Case #%d: %lld\n",k,ans);

}
return 0;
}