莫比乌斯反演升级版 HDU 4746

题目意思,给你一个Q,表示有Q组数据,对于每组数据给你三个数n, m, c
设f(i)表示为i的质因数个数,f(1) = 0
你要求出f(gcd(a, b)) <= c的对数
思路:
如果不考虑任何东西,最暴力最暴力的方法就是枚举a枚举b然后一个一个试,当然这么不行的。
我们需要用到Mobius函数优化下运算速度,设B(i) 为 i | gcd(x, y)的对数设A(i) 为gcd(x, y) == i的对数,由莫比乌斯函数可以得到A(i) = sigma(mu(d) * B(d * i))
这样的话速度加快了很多复杂度为O(n^2 * Q),但是还是会超时。
如果直接继续用分块的思想,对mu函数求一个前缀和,再加速加,复杂度O(n*sqrt(n)*Q)超时照样)因为这样,我们必须考虑分块的策略了。
在纸上画画,可以发现:

A(1) = mu(1)*B(1)+mu(2)*B(2)+…+mu(n)*B(n)
A(2) = A(1) = mu(1)*B(1)+mu(2)*B(2)+…+mu(n)*B(n)

A(d) = A(1) = mu(1)*B(1*d)+mu(2)*B(2*d)+…+mu(n)*B(n*d)
我们要计算的就是从1到d的某些满足条件的值。发现B始终只有5e5种状态,那么设
ans = A(1) + A(2) + A(3) + … + A(d) 提取B(1)到B(n)可以得到
ans = F(1)*B(1) + F(2)*B(2) + F(3)*B(3) + …+F(n)*B(n);
我们预处理出F数组,就可以很高效的分块加速了
因为5e5的最多素因子个数为18,因此枚举素因子个数,然后就能预处理完,很不幸,这样预处理也会超时,Orz真给跪。。
于是我学到了一种很神奇的方式,说不清楚,具体见代码吧

//
//  Created by Running Photon on 2016-03-23
//  Copyright (c) 2015 Running Photon. All rights reserved.
//
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <map>
#include <queue>
#include <string>
#include <sstream>
#include <set>
#include <vector>
#include <stack>
#define ALL(x) x.begin(), x.end()
#define INS(x) inserter(x, x,begin())
#define ll long long
#define CLR(x) memset(x, 0, sizeof x)
using namespace std;
const int inf = 0x3f3f3f3f;
const int MOD = 1e9 + 7;
const int maxn = 5e5 + 10;
const int maxv = 5e5 + 1;
const double eps = 1e-9;


int mu[maxn];
int npri[maxn];
ll prime[maxn];
int f[maxn];
int F[20][maxn];
int sum[20][maxn];
void init() {
    CLR(mu); CLR(npri);
    CLR(f);
    mu[1] = 1;
    int tot = 0;
    for(ll i = 2; i < maxv; i++) {
        if(!npri[i]) {
            f[i] = 1;
            mu[i] = -1;
            prime[++tot] = i;
        }
        for(int j = 1; j <= tot && prime[j] * i < maxv; j++) {
            npri[prime[j] * i] = 1;
            if(i % prime[j]) {
                mu[i * prime[j]] = -mu[i];
            }
            else {
                mu[i * prime[j]] = 0;
//                break;
            }
//f函数的处理
            f[i * prime[j]] = f[i] + 1;
        }
    }
//    for(int i = 1; i <= 100; i++) {
//        printf("f[%d] = %d\n", i, f[i]);
//    }
    for(int i = 1; i < maxv; i++) {
        for(int j = i; j < maxv; j += i) {
            F[f[i]][j] += mu[j / i];
        }
    }
    for(int i = 1; i < 19; i++) {
        for(int j = 1; j < maxv; j++) {
            F[i][j] += F[i - 1][j];
        }
    }
    for(int i = 1; i < maxv; i++) {
        for(int j = 0; j < 19; j++) {
            sum[j][i] = sum[j][i - 1] + F[j][i];
        }
    }
}
ll solve(ll a, ll b, ll c) {
    ll ans = 0;
    for(ll i = 1, la = 0; i <= a; i = la + 1) {
        la = min(a / (a / i), b / (b / i));
        ans += (sum[c][la] - sum[c][i - 1]) * (a / i) * (b / i);
    }
    return ans;
}
int main() {
#ifdef LOCAL
    freopen("in.txt", "r", stdin);
//  freopen("out.txt","w",stdout);
#endif
//  ios_base::sync_with_stdio(0);
    int T;
    scanf("%d", &T);
    ll a, b, c;
    init();
    while(T--) {
        scanf("%lld%lld%lld", &a, &b, &c);
        if(c > 18) {
            printf("%lld\n", a * b);
            continue;
        }
        if(a > b) swap(a, b);
        ll ans = solve(a, b, c);
        printf("%lld\n", ans);
    }

    return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值