【BZOJ】3561: DZY Loves Math VI

题意

\(\sum_{i=1}^{n} \sum_{j=1}^{m} lcm(i, j)^{gcd(i, j)}\)\(n, m<=500000\)

分析

很显然要死推莫比乌斯

题解

\(n \le m\)

\[ \begin{aligned} ans & = \sum_{i=1}^{n} \sum_{j=1}^{m} lcm(i, j)^{gcd(i, j)} \\ & = \sum_{i=1}^{n} \sum_{j=1}^{m} (\frac{ij}{gcd(i, j)})^{gcd(i, j)} \\ & = \sum_{d=1}^{n} \sum_{i=1}^{a} \sum_{j=1}^{b} \left( \frac{ijdd}{d} \right)^{d} \sum_{k|(i, j)} \mu(k) \ \ \left( a=\left \lfloor \frac{n}{d} \right \rfloor, b=\left \lfloor \frac{m}{d} \right \rfloor \right) \\ & = \sum_{d=1}^{n} d^d \sum_{k=1}^{a} \mu(k) \sum_{k|i}^{a} i^d \sum_{k|j}^{b} j^d \\ & = \sum_{d=1}^{n} d^d \sum_{k=1}^{a} \mu(k) k^{2d} \sum_{i=1}^{\left \lfloor \frac{a}{k} \right \rfloor} i^d \sum_{j=1}^{\left \lfloor \frac{b}{k} \right \rfloor} j^d \\ & = \sum_{d=1}^{n} d^d \sum_{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor} \mu(k) k^{2d} \sum_{i=1}^{\left \lfloor \frac{n}{kd} \right \rfloor} i^d \sum_{j=1}^{\left \lfloor \frac{m}{kd} \right \rfloor} j^d \\ \end{aligned} \]

于是我们对于每一个\(d\),暴力维护一下\(\mu(k) k^{2d}\),暴力维护一下\(\displaystyle \sum_{i=1}^{\left \lfloor \frac{m}{kd} \right \rfloor} j^d\),总复杂度\(O(nlogn)\)

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mo=1000000007, N=500005;
int mu[N], p[N], pcnt, np[N], c[N], C[N], b[N];
int ipow(int a, int b) {
    int x=1;
    for(; b; b>>=1, a=(ll)a*a%mo) if(b&1) x=(ll)x*a%mo;
    return x;
}
void init(int n) {
    mu[1]=1;
    for(int i=2; i<=n; ++i) {
        if(!np[i]) {
            p[pcnt++]=i;
            mu[i]=-1;
        }
        for(int j=0; j<pcnt; ++j) {
            int t=p[j]*i;
            if(t>n) break;
            np[t]=1;
            if(i%p[j]==0) {
                mu[t]=0;
                break;
            }
            mu[t]=-mu[i];
        }
    }
}
int main() {
    int n, m, ans=0;
    scanf("%d%d", &n, &m);
    if(n>m) {
        swap(n, m);
    }
    init(n);
    for(int i=1; i<=m; ++i) {
        c[i]=1;
    }
    for(int d=1; d<=n; ++d) {
        int A=ipow(d, d);
        int nn=n/d, mm=m/d;
        for(int k=1; k<=mm; ++k) {
            c[k]=(ll)c[k]*k%mo;
            C[k]=C[k-1]+c[k];
            if(C[k]>=mo) {
                C[k]-=mo;
            }
        }
        int temp=0;
        for(int k=1; k<=nn; ++k) if(mu[k]) {
            temp+=(ll)c[k]*c[k]%mo*C[nn/k]%mo*C[mm/k]%mo*mu[k];
            if(temp>=mo) {
                temp-=mo;
            }
            if(temp<0) {
                temp+=mo;
            }
        }
        ans+=(ll)A*temp%mo;
        if(ans>=mo) {
            ans-=mo;
        }
    }
    printf("%d\n", ans);
    return 0;
}

转载于:https://www.cnblogs.com/iwtwiioi/p/4986030.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值