对于一个k进制下的分数
nm((n,m)=1)
n
m
(
(
n
,
m
)
=
1
)
,求他的小数部分的过程就是不断的执行n=n*k%m这样一个过程
若他是一个循环小数,设他小数点后第一位为x,当
x≠0
x
≠
0
时,有
x≡xkl(mod m)
x
≡
x
k
l
(
m
o
d
m
)
,于是
kl≡1,(k,m)=1
k
l
≡
1
,
(
k
,
m
)
=
1
,
x=0
x
=
0
时此时
m=1
m
=
1
,也有
(k,m)=1
(
k
,
m
)
=
1
于是我们要求的其实是这么个柿子
为了方便下文中的 xy x y 都代表 ⌊xy⌋ ⌊ x y ⌋
∑mj=1[(j,k)=1]∑ni=1∑d|i,d|jμ(d) ∑ j = 1 m [ ( j , k ) = 1 ] ∑ i = 1 n ∑ d | i , d | j μ ( d )
∑dμ(d)nd∑j[(dj,k)=1] ∑ d μ ( d ) n d ∑ j [ ( d j , k ) = 1 ]
∑dμ(d)nd[(d,k)=1]∑j[(j,k)=1] ∑ d μ ( d ) n d [ ( d , k ) = 1 ] ∑ j [ ( j , k ) = 1 ]
记
g(n,k)=∑ni=1(i,k)=1
g
(
n
,
k
)
=
∑
i
=
1
n
(
i
,
k
)
=
1
,因为
(a,b)=(a%b,b)
(
a
,
b
)
=
(
a
%
b
,
b
)
所以有
g(n,k)=nkg(k,k)+g(n%k,k)
g
(
n
,
k
)
=
n
k
g
(
k
,
k
)
+
g
(
n
%
k
,
k
)
,于是
g(n,k)
g
(
n
,
k
)
可以在
O(k)
O
(
k
)
的预处理后
O(1)
O
(
1
)
求
∑d[(d,k)=1]μ(d)ndg(md,k) ∑ d [ ( d , k ) = 1 ] μ ( d ) n d g ( m d , k )
我们可以对 nd,md n d , m d 分块,那么我们需要快速求出前面那个东西的前缀和
记 f(n,k)=∑ni=1[(i,k)=1]μ(i) f ( n , k ) = ∑ i = 1 n [ ( i , k ) = 1 ] μ ( i ) ,不妨设 k k 的每个质因子次数不超过1,当时,记 k=pq( (p,q)=1 ) k = p q ( ( p , q ) = 1 )
f(n,k)=f(n,q)−∑ni=1[(i,q)=1,p|i]μ(i) f ( n , k ) = f ( n , q ) − ∑ i = 1 n [ ( i , q ) = 1 , p | i ] μ ( i )
=f(n,q)−∑n/pi=1[(ip,q)=1]μ(ip) = f ( n , q ) − ∑ i = 1 n / p [ ( i p , q ) = 1 ] μ ( i p )
因为 μ(ip) μ ( i p ) 当 p2|ip p 2 | i p 即 p|i p | i 时=0没有贡献,所以
f(n,k)=f(n,q)−μ(p)∑n/pi=1[(i,p)=1,(i,q)=1]μ(i) f ( n , k ) = f ( n , q ) − μ ( p ) ∑ i = 1 n / p [ ( i , p ) = 1 , ( i , q ) = 1 ] μ ( i )
=f(n,q)+∑n/pi=1[(i,k)=1]μ(i) = f ( n , q ) + ∑ i = 1 n / p [ ( i , k ) = 1 ] μ ( i )
=f(n,q)+f(n/p,k) = f ( n , q ) + f ( n / p , k )
我们可以递归求
f
f
的值,用一个哈希表记忆化,因为每次递归下去两个值有一个会
/p
/
p
,所以总的状态数不会太多(感性/理性分析一下会发现确实不会很多)
边界的话
k=1
k
=
1
时是个杜教筛,建议预处理一下
n<=1000
n
<=
1000
时候的
f(n,k)
f
(
n
,
k
)
这样会快很多
code:
#include<set>
#include<map>
#include<deque>
#include<queue>
#include<stack>
#include<cmath>
#include<ctime>
#include<bitset>
#include<string>
#include<vector>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<climits>
#include<complex>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;
int gcd(int a,int b){return !a?b:gcd(b%a,a);}
const int maxn = 5100000;
const int hashmod = 2333333;
int K;
int p[maxn>>2],pri,mu[maxn],smu[maxn],mp[maxn];
int g[2005],ff[2005][2005];
void pre()
{
smu[1]=mu[1]=1;
for(int i=2;i<maxn;i++)
{
if(!mp[i]) p[++pri]=i,mu[i]=-1,mp[i]=i;
for(int j=1,k=p[j]*i;k<maxn;j++,k=p[j]*i)
{
mp[k]=p[j];
if(i%p[j]==0) { mu[k]=0;break; }
mu[k]=-mu[i];
}
smu[i]=smu[i-1]+mu[i];
}
int tk=1;
while(K>1)
{
int pi=mp[K]; tk*=pi;
while(K%pi==0) K/=pi;
}K=tk;
for(int i=1;i<=K;i++) g[i]=g[i-1]+(gcd(i,K)==1);
for(int k=1;k<=K;k++) if(K%k==0)
{
for(int i=1;i<2005;i++)
ff[k][i]=ff[k][i-1]+(gcd(k,i)==1?mu[i]:0);
}
}
struct edge{int y1,y2,c,nex;};
struct HashTable
{
edge a[hashmod<<2]; int len,fir[hashmod];
int Hash(int x){return x%hashmod;}
void ins(int n,int k,int c)
{
int x=Hash(n+k);
if(!fir[x]) fir[x]=x,a[x]=(edge){n,k,c,0};
else a[++len]=(edge){n,k,c,fir[x]},fir[x]=len;
}
int query(int n,int kk,int &ok)
{
int x=Hash(n+kk);
ok=1;
for(int k=fir[x];k;k=a[k].nex)
if(a[k].y1==n&&a[k].y2==kk) return a[k].c;
ok=0;
return -1;
}
}hf,hu;
int calg(int n,int k){return g[k]*(n/k)+g[n%k];}
int calu(int n)
{
if(n<maxn) return smu[n];
int ok,c=hu.query(n,0,ok);
if(ok==0)
{
c=1;
for(int i=2,r;i<=n;i=r+1)
{
r=n/(n/i);
c=c-calu(n/i)*(r-i+1);
}
hu.ins(n,0,c);
}
return c;
}
int calf(int n,int k)
{
if(k==1) return calu(n);
if(n<2005) return ff[k][n];
int ok,c=hf.query(n,k,ok);
if(ok==0)
{
int pi=mp[k];
c=calf(n,k/pi)+calf(n/pi,k);
hf.ins(n,k,c);
}
return c;
}
int n,m;
int main()
{
//freopen("tmp.in","r",stdin);
//freopen("tmp.out","w",stdout);
hf.len=hu.len=hashmod-1;
scanf("%d%d%d",&n,&m,&K); pre();
int u=min(n,m); ll ans=0;
for(int i=1,r;i<=u;i=r+1)
{
r=min(n/(n/i),m/(m/i));
ans+=(ll)calg(m/i,K)*(n/i)*(calf(r,K)-calf(i-1,K));
}
printf("%lld\n",ans);
return 0;
}