完全是套用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\n",++k,ans);
    }
    return 0;
}