BZOJ

Solution

$$k^x\equiv 1 \pmod b$$，有解的条件是 $$\gcd(k,b)=1$$

$ans=F(n,m,k)=\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=1][\gcd(j,k)=1]$

$F(n,m,k)=\sum_{j=1}^m[(j,k)=1]\sum_{i=1}^n[(i,j)=1]$

$=\sum_{d|k}\mu(d)\sum_{i=1}^n\sum_{j=1}^{m/d}[(i,jd)=1]$

$=\sum_{d|k}\mu(d)\sum_{i=1}^n[(i,d)=1]\sum_{j=1}^{m/d}[(i,j)=1]=\sum_{d|k}\mu(d)F(m/d,n,d)$

$F(n,m,1)=\sum_{i=1}^n\sum_{j=1}^m[(i,j)=1]=\sum_{d=1}^{\min(n,m)}\mu(d)\frac n d\frac m d$

$$\mu$$ 的前缀和以及 $$F$$ 分别进行记忆化即可，可能需要注意卡常。

Code

#include <algorithm>
#include <cstdio>
using namespace std;
typedef long long ll;
const int N=1000000,M=500000,P=999917,INF=1e9;
template <typename Tp> inline int getmin(Tp &x,Tp y){return y<x?x=y,1:0;}
template <typename Tp> inline int getmax(Tp &x,Tp y){return y>x?x=y,1:0;}
template <typename Tp> inline void read(Tp &x)
{
x=0;int f=0;char ch=getchar();
while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
if(ch=='-') f=1,ch=getchar();
while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
if(f) x=-x;
}
struct data{
int n,m,k;
bool operator == (const data &b)const{return n==b.n&&m==b.m&&k==b.k;}
int operator % (const int &b)const{return (n*1000000009ll+m*879121319ll+k)%b;}
};
struct Hasha{
int& operator [] (int y)
{
int x=y%P;
ins(x,y);return w[p];
}
}ma;
struct Hashb{
ll& operator [] (data y)
{
int x=y%P;
ins(x,y);return w[p];
}
}mb;
int n,m,k,tot,pri[N+10],mu[N+10],vis[N+10];
int min(int x,int y){return x<y?x:y;}
void init()
{
mu[1]=1;
for(int i=2;i<=N;i++)
{
if(!vis[i]) pri[++tot]=i,mu[i]=-1;
for(int j=1;j<=tot&&i*pri[j]<=N;j++)
{
vis[i*pri[j]]=1;
if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
mu[i*pri[j]]=-mu[i];
}
}
for(int i=2;i<=N;i++) mu[i]+=mu[i-1];
}
int S(int n)
{
if(n<=N) return mu[n];
int &res=ma[n];
if(res!=-INF) return res;
res=1;
for(int i=2,j;i<=n;i=j+1)
{
j=n/(n/i);
res-=(j-i+1)*S(n/i);
}
return res;
}
ll F(int n,int m,int k)
{
if(n==0||m==0) return 0ll;
if(k==1)
{
if(n>m) swap(n,m);
data hs=(data){n,m,k};ll &res=mb[hs];
if(~res) return res;
res=0ll;
for(int i=1,j;i<=n;i=j+1)
{
j=min(n/(n/i),m/(m/i));
res+=(ll)(S(j)-S(i-1))*(n/i)*(m/i);
}
return res;
}
data hs=(data){n,m,k};ll &res=mb[hs];
if(~res) return res;
int t;res=0ll;
for(int i=1;i*i<=k;i++)
if(k%i==0)
{
t=k/i;res+=(mu[i]-mu[i-1])*F(m/i,n,i);
if(i*i!=k) res+=(mu[t]-mu[t-1])*F(m/t,n,t);
}
return res;
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt","r",stdin);
#endif
init();