完全是套用polya模版……
#include<iostream>
#include<stdio.h>
#include<algorithm>
#include<iomanip>
#include<cmath>
#include<string>
using namespace std;
int prime[30005],m,mod=1000000007;
bool f[30005];
__int64 extend(__int64 a,__int64 b,__int64 &x,__int64 &y)
{
__int64 d;
if(b==0)
{
x=1;y=0;
return a;
}
else
{
d=extend(b,a%b,x,y);
__int64 t=x;
x=y;
y=t-a/b*y;
}
return d;
}
void init()
{
int i,j;
m=0;
for(i=2;i<=30000;i++)
if(f[i]==0)
{
prime[m++]=i;
for(j=i*i;j<=30000;j+=i)
f[j]=1;
}
}
__int64 euler(__int64 n)
{
__int64 ans=1,i;
for(i=0;i<m&&prime[i]*prime[i]<=n;i++)
{
if(n%prime[i]==0)
{
ans*=prime[i]-1;
n/=prime[i];
while(n%prime[i]==0)
{
ans*=prime[i];
n/=prime[i];
}
}
}
if(n>1) ans*=n-1;
return ans%mod;
}
__int64 pows(__int64 a,__int64 b)
{
__int64 ans=1;
a=a%mod;
while(b)
{
if(b&1) ans=(ans*a)%mod;
b>>=1;
a=(a*a)%mod;
}
return ans;
}
int main()
{
init();
int t,i,j,n,s,k=0;
__int64 ans;
cin>>t;
while(t--)
{
cin>>n>>s;
ans=0;
for(i=1;i<=s;i++)
{
if(s%i==0)
ans=(ans+pows(n,i)*euler(s/i))%mod;
}
if(s&1)
{
ans=(ans+pows(n,(s+1)/2)*s)%mod;
}
else
{
ans=(ans+pows(n,s/2)*(s/2))%mod;
ans=(ans+pows(n,s/2+1)*(s/2))%mod;
}
__int64 x,y;
extend(s<<1,mod,x,y);
x=(x%mod+mod)%mod;
ans=(ans*x)%mod;
printf("Case #%d: %d
",++k,ans);
}
return 0;
}