可以理解为上一道题的扩展板..
然后我们就可以YY出这样一个式子
${\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^cd(ijk)=\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c[gcd(i,j)=gcd(i,k)=gcd(j,k)=1]\lfloor\frac{a}{i}\rfloor\lfloor\frac{b}{j}\rfloor\lfloor\frac{c}{k}\rfloor}$
然后我们枚举第一维,排除掉不和第一维互质的数大力反演就OK啦。
这里介绍另一种很神奇(麻烦)方法,当然这个和反演的关系就没那么大了:
首先设$f(k)=\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}[k=i \times j]$
然后得到$ans=\sum\limits_{i=1}^{a\times b} f(i)\sum\limits_{j=1}{d(i \times j)}$
$ans=\sum\limits_{i=1}^{a\times b} f(i)\sum\limits_{j=1}^{c}{d(i \times j)}$
$ans=\sum\limits_{i=1}^{a\times b} f(i)\sum\limits_{j=1}^{c}\sum\limits_{u=1}^{i}\sum\limits_{v=1}^{j}[gcd(u,v)=1]\lfloor\frac{i}{u}\rfloor\lfloor\frac{j}{v}\rfloor$
$ans=\sum\limits_{i=1}^{a\times b}\sum\limits_{j=1}^{c}[gcd(i,j)=1]\sum\limits_{u=1}^{\frac{a\times b}{i}}f(u\times i)\lfloor \frac{c}{j}\rfloor$
不妨设$S(i)=\sum\limits_{u=1}^{\frac{a\times b}{i}}f(u\times i)$
得到:
$ans=\sum\limits_{i=1}^{a\times b}\sum\limits_{j=1}^{c}[gcd(i,j)=1]S(i)\lfloor \frac{c}{j}\rfloor$
$ans=\sum\limits_{i=1}^{a\times b}\sum\limits_{j=1}^{c}\sum\limits_{d|gcd(i,j)}\mu(d)S(i)\lfloor \frac{c}{j}\rfloor$
不妨设$Q(k)=\sum\limits_{i=1}^{k}\lfloor \frac{k}{i} \rfloor$
得到:
$ans=\sum\limits_{i=1}^{c}\mu(i)\sum\limits_{j=1}^{\frac{a\times b}{i}}S(i\times j)\sum\limits_{k=1}^{\frac{c}{i}}Q(k)$
不妨设$P(k)=\sum\limits_{i=1}^{\frac{a\times b}{k}}S(i\times k)$
最后得到$ans=\sum\limits_{i=1}^{c}\mu(i)P(i)Q(i)$
大力预处理即可。
//CF235E //by Cydiater //2017.2.22 #include <iostream> #include <queue> #include <map> #include <cstring> #include <string> #include <cstdlib> #include <cstdio> #include <cmath> #include <ctime> #include <iomanip> #include <algorithm> #include <bitset> #include <set> #include <vector> #include <complex> using namespace std; #define ll int #define up(i,j,n) for(ll i=j;i<=n;i++) #define down(i,j,n) for(ll i=j;i>=n;i--) #define cmax(a,b) a=max(a,b) #define cmin(a,b) a=min(a,b) const ll MAXN=4e6+5; const ll mod=1073741824; inline ll read(){ char ch=getchar();ll x=0,f=1; while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)f=-1;ch=getchar();} while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();} return x*f; } ll A,B,C,AB,prime[MAXN],P[MAXN],S[MAXN],f[MAXN],Q[MAXN],miu[MAXN],cnt; bool vis[MAXN]; namespace solution{ void Prepare(){ A=read();B=read();C=read();AB=A*B; miu[1]=1; up(i,2,C){ if(!vis[i]){prime[++cnt]=i;miu[i]=-1;} up(j,1,cnt){ if(i*prime[j]>C)break; vis[i*prime[j]]=1; if(i%prime[j])miu[i*prime[j]]=-miu[i]; else break; } } up(i,1,A)up(j,1,B)f[i*j]++; up(i,1,AB)for(ll j=i;j<=AB;j+=i) (S[i]+=f[j])%=mod; up(i,1,AB)for(ll j=i;j<=AB;j+=i) (P[i]+=S[j])%=mod; ll pos; up(i,1,C){ for(ll j=1;j<=i;j=pos+1){ pos=i/(i/j); Q[i]+=(pos-j+1)*(i/j); (Q[i]+=mod)%=mod; } } } void Solve(){ ll ans=0; up(i,1,C){ (ans+=miu[i]*P[i]*Q[C/i])%=mod; (ans+=mod)%=mod; } cout<<ans<<endl; } } int main(){ //freopen("input.in","r",stdin); using namespace solution; Prepare(); Solve(); //cout<<"Time has passed:"<<1.0*clock()/1000<<"s!"<<endl; return 0; }