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

Memento

这是一间记忆仓库~~

 
 
 

日志

 
 

HDU 3441 Rotation  

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

  下载LOFTER 我的照片书  |
        这题思路还是很清晰的,两次Polya,第一次求出B*B的方块的方案数,然后再求K个方块和中间的小方块连在一起的方案数,外围的方块可以等效为珠子,可用颜色数为第一次Polya求出的方案数。有一个比较蛋疼的地方就是B*B*K+1=A*A,A^2的范围为10^18,所以不能直接质因数分解。把等式变为B*B*K=(A-1)*(A+1),我们分解A-1和A+1的质因数,然后就可以DFS枚举B并确定K,在对B,K分别求方案数。这样就解决了这个问题。中间过程有好几个,所以代码量不小,而且要相互之间组织好,不然很容易出bug。最后效率还不错,居然又登顶了- -

HDU 3441 Rotation - endlesscount - Memento
 

#include<cstdio>
#include<iostream>
#include<cstring>
#include<vector>
#include<algorithm>
#define mod 1000000007
#define N 32000
using namespace std;
bool notprime[N];
int prime[N],pp;
long long nn[ 64 ];
long long fac[ 50 ][ 3 ],fp;
long long ans;
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 a ){
int i,temp;
vector<int> fa;
for( i=0,temp=a-1;i<pp&&temp>1;i++ )
for( ;temp%prime[ i ]==0;temp/=prime[ i ],fa.push_back( prime[ i ] ) );
if( temp>1 ) fa.push_back( temp );
for( i=0,temp=a+1;i<pp&&temp>1;i++ )
for( ;temp%prime[ i ]==0;temp/=prime[ i ],fa.push_back( prime[ i ] ) );
if( temp>1 ) fa.push_back( temp );
sort( fa.begin( ),fa.end( ) );
fp=0;
memset( fac,0,sizeof( fac ) );
fac[ fp ][ 0 ]=fa[ 0 ],fac[ fp ][ 1 ]=1,fac[ fp++ ][ 2 ]=fa[ 0 ];
for( int i=1;i<fa.size( );i++ ){
if( fa[ i ]==fa[ i-1 ] ){
fac[ fp-1 ][ 1 ]++;
fac[ fp-1 ][ 2 ]*=fa[ i ];
}
else{
fac[ fp ][ 0 ]=fac[ fp ][ 2 ]=fa[ i ];
fac[ fp++ ][ 1 ]=1;
}
}
}
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;
}
}
long long pow( long long b ){
long long ret=1;
for( int i=0;b;b>>=1,i++ )
if( b&1 )
ret=( ret*nn[ i ] )%mod;
return ret;
}
void DFS( long long g,long long phi,int d ){
if( d==fp ){
phi%=mod;
ans=( ans+phi*pow( g ) )%mod;
return;
}
for( long long i=0,p=fac[d][2];i<=fac[d][1];i++,p/=fac[d][0] )
DFS( fac[d][2]/p*g,phi*(p-p/fac[d][0]),d+1 );
}
void calc( long long b,long long c ){
long long x,y;
ex_gcd( x,y,4,mod );
x=( x%mod+mod )%mod;
nn[ 0 ]=c;
for( int i=1;i<64;i++ ) nn[ i ]=( nn[ i-1 ]*nn[ i-1 ] )%mod;
c=( pow( b*b )+2*pow( b*b/4+( b%2?1:0 ) )+pow( b*b/2+( b%2?1:0 ) ) )%mod;
c=( c*x )%mod;
nn[ 0 ]=c;
for( int i=1;i<64;i++ ) nn[ i ]=( nn[ i-1 ]*nn[ i-1 ] )%mod;
}
long long cal( long long a,long long b,long long k,long long c ){
calc( b,c );
ans=0;
DFS( 1,1,0 );
long long x,y;
ex_gcd( x,y,k,mod );
x=( x%mod+mod )%mod;
ans=( ans*x )%mod;
ans=( ans*c )%mod;
return ans;
}
void dfs( long long &res,long long a,long long b,long long c,int d ){
if( d==fp ){
res=( res+cal( a,b,( a*a-1 )/b/b,c ) )%mod;
}
else{
for( long long i=0,p=1;2*i<=fac[ d ][ 1 ];i++,p*=fac[ d ][ 0 ] ){
fac[ d ][ 1 ]-=2*i;
fac[ d ][ 2 ]/=p*p;
dfs( res,a,b*p,c,d+1 );
fac[ d ][ 1 ]+=2*i;
fac[ d ][ 2 ]*=p*p;
}
}
}
int main(){
// freopen( "in.in","r",stdin );
init_prime();
int a,c,T;
cin>>T;
for(int cas=1;cas<=T;cas++){
printf( "Case %d: ",cas );
cin>>a>>c;
if( a==1 ) cout<<c<<endl;
else{
fact( a );
long long res=0;
dfs( res,a,1,c,0 );
cout<<res<<endl;
}
}
}


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

历史上的今天

评论

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

页脚

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