Mophues HDU - 4746 (莫比乌斯反演)

Mophues

\[ Time Limit: 10000 ms\quad Memory Limit: 262144 kB \]

题意

求出满足 \(gcd\left(a,b\right) = k\),其中\(1\leq a\leq n,1\leq b \leq m\)\(k\) 的因子数 \(\leq P\)

思路

\(g\left(x\right)\) 表示 \(gcd\left(a, b\right) | x\) 的对数
\(f\left(x\right)\) 表示 \(gcd\left(a, b\right) = x\) 的对数
根据莫比乌斯反演有
\[ f\left(n\right) = \sum_{n|d} g\left(d\right)\\ g\left(n\right) = \sum_{n|d} \mu\left(\frac{d}{n}\right) f\left(d\right) \\ \]
根据题意
\[ f\left(x\right) = \lfloor\frac{n}{x}\rfloor \lfloor\frac{m}{x}\rfloor \\ \]
那么可以得到
\[ \begin{aligned} ans &= \sum_{k\in 条件} \sum_{k|d} \mu\left(\frac{d}{k}\right) \lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor\\ &= \sum_{d}\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor \sum_{k|d}\mu\left(\frac{d}{k}\right) \left(k \in 条件\right) \end{aligned} \]
\(n \leq 10^{5}\),其中最多的因子数不会超过 \(19\)
\(c[d][p]\) 表示 \(\sum_{k|d}\mu\left(\frac{d}{k}\right)\),其中 \(k\) 的因子数\(\leq p\)。通过 \(nlogn\) 暴力打表出 \(p<=19\) 的情况,当 \(p>19\) 时,所有的对数都满足条件。
由于 \(\lfloor\frac{n}{d}\rfloor \lfloor\frac{m}{d}\rfloor\) 的存在,需要整除分块,所以对于求出来的 \(c[d][p]\) 还需要统计前缀和。

/***************************************************************
    > File Name    : a.cpp
    > Author       : Jiaaaaaaaqi
    > Created Time : 2019年07月17日 星期三 14时38分33秒
 ***************************************************************/

#include <map>
#include <set>
#include <list>
#include <ctime>
#include <cmath>
#include <stack>
#include <queue>
#include <cfloat>
#include <string>
#include <vector>
#include <cstdio>
#include <bitset>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define  lowbit(x)  x & (-x)
#define  mes(a, b)  memset(a, b, sizeof a)
#define  fi         first
#define  se         second
#define  pii        pair<int, int>

typedef unsigned long long int ull;
typedef long long int ll;
const int    maxn = 5e5 + 10;
const int    maxm = 1e5 + 10;
const ll     mod  = 1e9 + 7;
const ll     INF  = 1e18 + 100;
const int    inf  = 0x3f3f3f3f;
const double pi   = acos(-1.0);
const double eps  = 1e-8;
using namespace std;

int n, m, P;
int cas, tol, T;

int pri[maxn], mob[maxn], cnt[maxn];
ll c[maxn][20];
bool ispri[maxn];

void handle() {
    mes(pri, 0), mes(ispri, 1);
    tol = 0;
    int mx = 5e5;
    mob[1] = 1;
    cnt[1] = 0;
    for(int i=2; i<=mx; i++) {
        if(ispri[i]) {
            pri[++tol] = i;
            mob[i] = -1;
            cnt[i] = 1;
        }
        for(int j=1; j<=tol&&i*pri[j]<=mx; j++) {
            ispri[i*pri[j]] = false;
            cnt[i*pri[j]] = cnt[i]+1;
            if(i%pri[j] == 0) {
                mob[i*pri[j]] = 0;
                break;
            } else {
                mob[i*pri[j]] = -mob[i];
            }
        }
    }
    mes(c, 0);
    for(int i=1; i<=mx; i++) {
        for(int j=i; j<=mx; j+=i) {
            c[j][cnt[i]] += mob[j/i];
        }
    }
    for(int i=1; i<=mx; i++) {
        for(int j=1; j<=19; j++) {
            c[i][j] += c[i][j-1];
        }
    }
    for(int i=1; i<=mx; i++) {
        for(int j=0; j<=19; j++) {
            c[i][j] += c[i-1][j];
        }
    }
}

int main() {
    handle();
    scanf("%d", &T);
    while(T--) {
        scanf("%d%d%d", &n, &m, &P);
        if(n > m)   swap(n, m);
        if(P > 19) {
            printf("%lld\n", 1ll*n*m);
            continue;
        }
        ll ans = 0;
        int l = 0, r;
        for(int i=1; i<=n; i++) {
            r = min(n/(n/i), m/(m/i));
            ans += 1ll*(n/i)*(m/i)*(c[r][P]-c[l][P]);
            i = r;
            l = r;
        }
        printf("%lld\n", ans);
    }
    return 0;
}

转载于:https://www.cnblogs.com/Jiaaaaaaaqi/p/11201772.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值