[Polya计数] SPOJ 422 TRANSP2 - Transposing is Even More Fun

以下来自JYT的ppt

n个点的环 m种颜色 本质不同的染色为

$$\sum_{d|n} {\phi ({n \over d })*f_d}$$

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

inline char nc(){
  static char buf[100000],*p1=buf,*p2=buf;
  return p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++;
}
inline void read(int &x){
  char c=nc(),b=1;
  for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
  for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

const int N=1000005;
const int P=1000003;

int num,prime[N];
int d[N],phi[N];
int _pow[N],inv[N];

inline void Pre(int n){
  _pow[0]=1; for (int i=1;i<=n;i++) _pow[i]=(_pow[i-1]<<1)%P;
  inv[1]=1; for (int i=2;i<=n;i++) inv[i]=(ll)(P-P/i)*inv[P%i]%P;
  phi[1]=1;
  for (int i=2;i<=n;i++){
    if (!d[i]) prime[++num]=i,d[i]=i,phi[i]=i-1;
    for (int j=1;j<=num && (ll)prime[j]*i<=n;j++){
      d[i*prime[j]]=prime[j];
      if (i%prime[j]==0){
    phi[i*prime[j]]=phi[i]*prime[j];
    break;
      }else
    phi[i*prime[j]]=phi[i]*(prime[j]-1);
    }
  }
}
int n,m;
int p[N],q[N],cnt;
ll Ans=0;

inline void dfs(int t,int cur){
  if (t==cnt+1){
    Ans+=(ll)phi[n/cur]*_pow[m*cur]%P;
    return;
  }
  for (int i=0;i<=q[t];i++){
    dfs(t+1,cur);
    cur*=p[t];
  }
}

inline ll Solve(){
  cnt=0;
  int last=0,t=n;
  while (t>1){
    if (d[t]!=last) p[++cnt]=d[t],q[cnt]=0,last=d[t];
    q[cnt]++;
    t/=d[t];
  }
  Ans=0;
  dfs(1,1); Ans%=P;
  return Ans*inv[n]%P;
}

int main(){
  int Q,g,a,b;
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  read(Q); Pre(1e6);
  while (Q--){
    read(a),read(b); if (!a || !b) { printf("0\n"); continue; }
    g=__gcd(a+b,b);
    n=(a+b)/g; m=g;
    printf("%lld\n",(_pow[a+b]+P-Solve())%P);
  }
  return 0;
}