SPOJ422：Transposing is Even More Fun（Polya）

传送门

题意：
给你一个2a∗2b 的矩阵，在内存中的存放方式是先存第一行的，再存第二行
的……每行也是从左到右存放。现在你想求它的转置矩阵（也是一样的储存方
式），但是只能用交换操作，问需要交换多少步。

题解：
直接贴论文了：

为了描述方便起见，先看一个例子：假设 a=5， b=3。那么考虑元素(12,1)，用二进制表示是(01010,001)，那么它原来的地址就是 01010 001，而新的地址就是 001 01010。可见这个转置操作其实就是把每个元素的地址循环向右移动了b位。
考虑地址的循环节如果有k个，那么答案就是2a+b−k。

问题转化为求地址的循环节，显然，两个地址若能通过右移b的某倍数步，那么他们处于一个循环。
相当于是求在右移b的倍数步这个置换群中找染色方案。
成功转化为Polya裸题，但是n√枚举约数会TLE。
考虑线筛的时候记录最小质因数，这样找出一个大数的所有质因数接近O(1).直接DFS找约数就好了。

#include<iostream>#include<cstdio>#include<vector>using namespace std;struct IO{    streambuf *ib,*ob;    inline void init(){        ios::sync_with_stdio(false);        cin.tie(NULL);cout.tie(NULL);        ib=cin.rdbuf();ob=cout.rdbuf();    }    inline int read(){        static char ch;static int i,f;        ch=ib->sbumpc();i=0,f=1;        while(!isdigit(ch)){            if(ch==-1)return false;            if(ch=='-')f=-1;            ch=ib->sbumpc();        }        while(isdigit(ch)){            i=(i+(i<<2)<<1)+ch-'0';            ch=ib->sbumpc();        }        return ((f>0)?i:-i);    }    inline void W(int x){        static int buf[50];        if(!x){ob->sputc('0');return;}        if(x<0){ob->sputc('-');x=-x;}        while(x){buf[++buf[0]]=x%10;x/=10;}        while(buf[0]){ob->sputc(buf[buf[0]--]+'0');}    }}io;const int Mod=1e6+3;inline int gcd(int x,int y){return (y?(gcd(y,x%y)):x);}const int Maxn=1e6;int ans,Primecnt,prime[Maxn+50],phi[Maxn+50],IsNotPrime[Maxn+50],_pow[Maxn],inv[Maxn],pr[Maxn];inline void seive(){    phi[1]=1;inv[1]=1;_pow[1]=2;pr[1]=1;    for(int i=2;i<=Maxn;i++){        _pow[i]=(_pow[i-1]<<1)%Mod;        inv[i]=(-1ll*inv[Mod%i]*(Mod/i))%Mod;        if(!IsNotPrime[i]){prime[++Primecnt]=i;phi[i]=i-1;pr[i]=i;}        for(int j=1;j<=Primecnt;j++){            int k=i*prime[j];            if(k>Maxn)break;            IsNotPrime[k]=1;            if(!(i%prime[j])){phi[k]=phi[i]*prime[j];pr[k]=pr[i];break;}            phi[k]=phi[i]*(prime[j]-1);pr[k]=prime[j];        }    }}int a,b,lim,g,factory[Maxn],cnt[Maxn],tail;inline void fac(int x){    tail=0;int last=0;    while(pr[x]!=1){        if(pr[x]==last)cnt[tail]++;        else cnt[++tail]=1,factory[tail]=pr[x],last=pr[x];        x/=pr[x];    }}inline void dfs(int pos,int d){    if(pos>tail){(ans+=1ll*phi[lim/d]*_pow[d*g]%Mod)%=Mod;return;}    dfs(pos+1,d);    for(int i=1;i<=cnt[pos];i++)d*=factory[pos],dfs(pos+1,d);}int main(){    io.init();seive();    int T=io.read();    while(T--){        ans=0;        a=io.read(),b=io.read();        if(a==0||b==0){io.W(0);io.ob->sputc('\n');continue;}        g=gcd(a+b,b),lim=(a+b)/g;        fac(lim);dfs(1,1);        ans=1ll*ans*inv[lim]%Mod;        ans=(_pow[a+b]-ans+Mod)%Mod;        io.W(ans);io.ob->sputc('\n');    }}