完全是套用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;

}