注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Memento

这是一间记忆仓库~~

 
 
 

日志

 
 

SPOJ 422 Transposing is Even More Fun  

2012-03-15 20:14:01|  分类: 组合数学 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
        给你一个2^a*2^b的矩阵,这个矩阵在内存中的存储方式是一行一行存储,要你求它的转置矩阵,每次操作都只能交换内存中的两个单元,求最少要多少步能完成任务。

       根据转置矩阵的定义,a[i][j]=a[j][i]。我们把这个式子用矩阵系数对应的在内存中的地址来表示。原矩阵系数a[i][j]在内存中的的地址为2^b*i+j(0<=i<2^a,0<=j<2^b),转置矩阵的系数a[i][j]在内存中的地址应该为2^a*i+j(0<=i<2^b,0<=j<2^a)。a[i][j]=a[j][i],所以我们是把地址为2^b*j+i的值赋给2^a*i+j。仔细观察这个式子,我们发现其实地址的变化是有规律的,因为0<=i<2^b,0<=j<2^a,所以i,j分别可以写成位数为b和位数为a的二进制形式。那么2^b*j+i---->2^a*i+j相当与把地址向右移了a位。由于地址数是有限的,所以每个地址经过若干次右移操作后都会变回自身,也就是说存在循环节。假设一个大小为n的循环节,我们可知用n-1次操作就可以把其全部安置正确。那么假设总共有k个循环节,那么我们需要的最小操作数=2^(a+b)-k。

       怎么样求k呢?仔细一想,我们发现k其实就是在右移这个操作所生成的置换群下,不同构的地址数。这就转化成了Polya求方案数的问题了。抽象出来就是一串有a+b个珠子的项链,我们用两种颜色对其进行染色,定义一次转换操作为逆时针移a个珠子的旋转,如果两种方案经过若干次转换后相同,那么就酸痛一种,求不同的方案数。我们首先想到最基本的项链染色计数问题就是a=1的时候的特殊情况,现在就是做了一个推广。设g=gcd(a,a+b),首先我们可以得到置换群的个数|G|=(a+b)/g,对于每次位i*g个位置的置换群其不动染色的个数C(f)=2^gcd(g*i,(a+b)/g)。所以总方案数k=1/|G|*sum{2^gcd(g*i,(a+b)/g),1<=i<=n/g}。这个公式我们可以用类似于POJ 2154 Color的方法来优化。所以k也在很快的时间就求出来了。

#include<cstdio>
#include<cstring>
#include<algorithm>
#define N 1009
#define mod 1000003
using namespace std;
int a,b;
bool notprime[ N ];
int prime[ N ],pp;
int fac[ 50 ][ 3 ],fp;
int nn[ 32 ];
void init_prime(){
for( int i=2;i<N;i++){
if(!notprime[i]) prime[ pp++ ]=i;
for(int j=0;j<pp&&i*prime[ j ]<N;j++){
notprime[ i*prime[ j ] ]=true;
if(i%prime[ j ]==0) break;
}
}
}
void fact(int n){
fp=0;
for(int i=0;i<pp&&n>1;i++)
if(n%prime[ i ]==0){
for(fac[fp][0]=prime[i],fac[fp][1]=0,fac[fp][2]=1;n%prime[i]==0;n/=prime[i],fac[fp][1]++,fac[fp][2]*=prime[i]);
fp++;
}
if(n>1) fac[fp][0]=n,fac[fp][1]=1,fac[fp++][2]=n;
}
int gcd( int a,int b ){
while( b ) b^=a^=b^=a%=b;
return a;
}
void ex_gcd(long long &x,long long &y,long long a,long long b){
if(!b) x=1,y=0;
else{
ex_gcd(x,y,b,a%b);
long long t=y;
y=x-a/b*y,x=t;
}
}
int pow(int a){
int ret=1;
for(int i=0;a;a>>=1,i++)
if(a&1)
ret=1ll*ret*nn[i]%mod;
return ret;
}
void DFS(int &res,int g,int phi,int d){
if(d==fp){
res=(1ll*phi*pow(g)+res)%mod;
return;
}
for(int i=0,p=fac[d][2];i<=fac[d][1];i++,p/=fac[d][0])
DFS(res,fac[d][2]/p*g,phi*(p-p/fac[d][0]),d+1);
}
int cal(){
if(!a||!b) return 0;
int G=gcd( a,a+b );
fact((a+b)/G);
int ret=0;
DFS(ret,G,1,0);
long long x,y;
ex_gcd(x,y,(a+b)/G,mod);
x%=mod;
ret=((pow(a+b)-x*ret)%mod+mod)%mod;
return ret;
}
int main(){
int cas;
init_prime();
nn[0]=2;
for(int i=1;i<32;i++) nn[ i ]=1ll*nn[ i-1 ]*nn[ i-1 ]%mod;
scanf("%d",&cas);
while(cas--){
scanf("%d%d",&a,&b);
printf("%d\n",cal());
}
}



  评论这张
 
阅读(561)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017