分析:
硬核数据范围欺诈。。。数据开大了可还行。。。
考场上想了三个小时都以为是骗分。。。
其实还是蛮简单的。就是有些坑
二维的平移,不一定平移(i,j)时,每个i*j的矩阵都必须一模一样。
考场上被这个毒了好久。。。
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
#define SF scanf
#define PF printf
#define MAXN 1030
#define MOD 1000000007
using namespace std;
typedef long long ll;
struct Matrix{
ll a[MAXN][MAXN];
}e,f,bg;
void mul(Matrix &x,Matrix &y,int n){
static Matrix tmp;
memset(tmp.a,0,sizeof tmp.a);
for(int i=0;i<(1<<n);i++)
for(int j=0;j<(1<<n);j++)
for(int k=0;k<(1<<n);k++)
(tmp.a[j][k]+=x.a[j][i]*y.a[i][k]%MOD)%=MOD;
for(int i=0;i<(1<<n);i++)
for(int j=0;j<(1<<n);j++)
x.a[i][j]=tmp.a[i][j];
}
ll fsp(ll x,int y){
ll res=1;
while(y){
if(y&1)
res=res*x%MOD;
x=x*x%MOD;
y>>=1;
}
return res;
}
ll phi(ll x){
ll res=1;
ll x1=x;
for(int i=2;1ll*i*i<=x;i++){
if(x%i==0){
res*=(i-1);
x/=i;
while(x%i==0)
x/=i,res*=i;
}
}
if(x!=1)
res*=(x-1);
// PF("[%lld %lld]\n",x1,res);
return res;
}
int tran(int x,int tot,int n){
return ((x>>n)|(x<<(tot-n)))&((1<<tot)-1);
}
ll query(int tot,int m,int n){
memset(bg.a,0,sizeof bg.a);
int m1=m;
static int blo[20];
int cnt=0;
for(int i=1;i<tot;i++)
blo[++cnt]=(1<<i)|(1<<(i-1));
if(tot>1)
blo[++cnt]=1|(1<<(tot-1));
static int t[MAXN];
memset(t,0,sizeof t);
t[0]=1;
for(int i=0;i<(1<<tot);i++){
int pos=i&(-i);
for(int j=1;j<=cnt;j++)
if((i&blo[j])==blo[j]&&(blo[j]&pos))
t[i]+=t[i^blo[j]];
}
for(int i=0;i<(1<<tot);i++)
for(int j=0;j<(1<<tot);j++){
if(i&j)
continue;
bg.a[i][j]+=t[((1<<tot)-1)^i^j];
}
for(int i=0;i<(1<<tot);i++)
for(int j=0;j<(1<<tot);j++)
f.a[i][j]=e.a[i][j]=bg.a[i][j];
m--;
while(m){
if(m&1)
mul(f,e,tot);
mul(e,e,tot);
m>>=1;
}
ll ans=0;
for(int i=0;i<(1<<tot);i++)
ans=(ans+f.a[i][tran(i,tot,n)])%MOD;
// PF("[%d %d %d %lld]\n",tot,n,m1,ans);
return ans;
}
int gcd(int x,int y){
if(y==0)
return x;
return gcd(y,x%y);
}
ll calc(int n,int m,int i,int j){
int kd=gcd(n/i,m/j)*i;
return 1ll*query(kd,j,i)*phi(n/i)%MOD*phi(m/j)%MOD;
}
ll solve(int n,int m){
ll ans=0;
for(int i=1;i<=n;i++)
for(int j=1;1ll*j*j<=m;j++)
if(n%i==0&&m%j==0){
ans=(ans+calc(n,m,i,j))%MOD;
if(j*j!=m)
ans=(ans+calc(n,m,i,m/j))%MOD;
}
ans=ans*fsp(1ll*n*m%MOD,MOD-2)%MOD;
return ans;
}
int main(){
freopen("anziehung.in","r",stdin);
freopen("anziehung.out","w",stdout);
int n,m;
SF("%d%d",&n,&m);
PF("%lld",solve(n,m));
}