洛谷P3312 [SDOI2014]数表 积性函数+反演+树状数组

分析

我们设 s d ( i ) sd(i) sd(i) i i i的所有约数之和
如果不考虑a对答案的影响:

a n s = ∑ i = 1 n ∑ j = 1 m s d ( g c d ( i , j ) ) ans=\sum_{i=1}^{n}\sum_{j=1}^{m}sd(gcd(i,j)) ans=i=1nj=1msd(gcd(i,j))

正常的思路,我们枚举 g c d ( i , j ) = d gcd(i,j)=d gcd(i,j)=d

a n s = ∑ d = 1 m i n ( n , m ) s d ( d ) ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = = d ] ans=\sum_{d=1}^{min(n,m)}sd(d)\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d] ans=d=1min(n,m)sd(d)i=1nj=1m[gcd(i,j)==d]

i , j i,j i,j提取出个 d d d

a n s = ∑ d = 1 m i n ( n , m ) s d ( d ) ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ [ g c d ( i , j ) = = 1 ] ans=\sum_{d=1}^{min(n,m)}sd(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,j)==1] ans=d=1min(n,m)sd(d)i=1dnj=1dm[gcd(i,j)==1]

我们可以使用反演函数替换掉后面那个 [ g c d ( i , j ) = = 1 ] [gcd(i,j)==1] [gcd(i,j)==1]

a n s = ∑ d = 1 m i n ( n , m ) s d ( d ) ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ ∑ t ∣ g c d ( i , j ) μ ( t ) ans=\sum_{d=1}^{min(n,m)}sd(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{t \mid gcd(i,j)}\mu(t) ans=d=1min(n,m)sd(d)i=1dnj=1dmtgcd(i,j)μ(t)

i , j i,j i,j提取出个 t t t

a n s = ∑ d = 1 m i n ( n , m ) s d ( d ) ∑ t = 1 m i n ( ⌊ n d ⌋ , ⌊ m d ⌋ ) ⌊ n t d ⌋ ⌊ m t d ⌋ ans=\sum_{d=1}^{min(n,m)}sd(d)\sum_{t=1}^{min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\lfloor\frac{n}{td}\rfloor\lfloor\frac{m}{td}\rfloor ans=d=1min(n,m)sd(d)t=1min(dn,dm)tdntdm

我们令 T = t d T=td T=td

a n s = ∑ T = 1 m i n ( n , m ) ⌊ n T ⌋ ⌊ n T ⌋ ∑ d ∣ T μ ( T d ) ⋅ s d ( d ) ans=\sum_{T=1}^{min(n,m)}\lfloor\frac{n}{T}\rfloor\lfloor\frac{n}{T}\rfloor\sum_{d \mid T}\mu(\frac{T}{d}) \cdot sd(d) ans=T=1min(n,m)TnTndTμ(dT)sd(d)

我们发现 ∑ d ∣ T μ ( T d ) ⋅ s d ( d ) \sum_{d \mid T}\mu(\frac{T}{d}) \cdot sd(d) dTμ(dT)sd(d)是积性函数,可以用线性筛求出

现在,我们有a的限制:

a n s = ∑ T = 1 m i n ( n , m ) ⌊ n T ⌋ ⌊ n T ⌋ ∑ d ∣ T μ ( T d ) ⋅ s d ( d ) ( d < = a ) ans=\sum_{T=1}^{min(n,m)}\lfloor\frac{n}{T}\rfloor\lfloor\frac{n}{T}\rfloor\sum_{d \mid T}\mu(\frac{T}{d}) \cdot sd(d)\quad(d<=a) ans=T=1min(n,m)TnTndTμ(dT)sd(d)(d<=a)

这样的话,我们只能离线的对所有查询进行处理。首先我们根据a的大小对所以查询排序,a每变大一点,我们都可以插入符合条件的 s d sd sd,这个我们可以用树状数组维护
对于一个符合条件的d,我们可以给d的所有倍数插入相应的值

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
template <typename T>
void out(T x) { cout << x << endl; }
ll fast_pow(ll a, ll b, ll p) {ll c = 1; while(b) { if(b & 1) c = c * a % p; a = a * a % p; b >>= 1;} return c;}
ll exgcd(ll a, ll b, ll &x, ll &y) { if(!b) {x = 1; y = 0; return a; } ll gcd = exgcd(b, a % b, y, x); y-= a / b * x; return gcd; }
const ll N = 1e5 + 100;
const ll mod = 1ll << 31;
struct node
{
    ll n, m, a, id;
    bool operator < (const node b) const
    {
        return a < b.a;
    }
}q[N], dd[N];
ll prime[N], mu[N], tot;
ll sd[N], sp[N];
bool mark[N];
void get_mu_sd()
{
    mark[0] = mark[1] = true;
    tot = 0;
    mu[1] = sd[1] = sp[1] = 1;
    for(ll i = 2; i < N; i ++)
    {
        if(!mark[i])
        {
            prime[tot ++] = i;
            mu[i] = -1;
            sd[i] = sp[i] = i + 1;
        }
        for(ll j = 0; j < tot && i * prime[j] < N; j ++)
        {
            mark[i * prime[j]] = true;
            if(i % prime[j] == 0)
            {
                sd[i * prime[j]] = sd[i] / sp[i] * (sp[i] * prime[j] + 1);
                sp[i * prime[j]] = sp[i] * prime[j] + 1;
                break;
            }
            sd[i * prime[j]] = (prime[j] + 1) * sd[i];
            sp[i * prime[j]] = prime[j] + 1;
            mu[i * prime[j]] = -mu[i];
        }
    }
    for(ll i = 1; i < N; i ++)
        dd[i].id = i, dd[i].a = sd[i];
    sort(dd + 1, dd + N);
}
ll tree[N], ans[N];
ll lowbit(ll x)
{
    return x & (-x);
}
void add(ll pos, ll val)
{
    for(ll i = pos; i < N; i += lowbit(i))
        tree[i] = tree[i] + val;
}
ll query(ll pos)
{
    ll num = 0;
    for(ll i = pos; i; i -= lowbit(i))
        num = num + tree[i];
    return num;
}
void add(ll x)
{
    for(ll i = 1; i * x < N; i ++)
        add(i * x, mu[i] * sd[x]);
}
ll cal(ll n, ll m)
{
    ll ans = 0;
    if(n > m)
        swap(n, m);
    for(ll i = 1, j = 1; i <= n; i = j + 1)
    {
        j = min(n / (n / i), m / (m / i));
        ans = (ans + (n / i) * (m / i) % mod * (query(j) - query(i - 1)) % mod + mod) % mod;
    }
    return ans;
}
int main()
{
    ios::sync_with_stdio(false);
    get_mu_sd();
    memset(tree, 0, sizeof(tree));
    ll t;
    cin >> t;
    for(ll i = 1; i <= t; i ++)
    {
        cin >> q[i].n >> q[i].m >> q[i].a;
        q[i].id = i;
    }
    sort(q + 1, q + 1 + t);
    ll cnt = 1;
    for(ll i = 1; i <= t; i ++)
    {
        while(dd[cnt].a <= q[i].a && cnt < N)
            add(dd[cnt ++].id);
        ans[q[i].id] = cal(q[i].n, q[i].m);
    }
    for(ll i = 1; i <= t; i ++)
        cout << ans[i] << endl;
}

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值