SPOJ 4491PGCD

这是GCD八题中最难的一题,自己写的是各种超时,只好参考网上大牛的解题报告,研究了一上午,终于给搞懂了.....

题目要求:gcd(i,j)=素数的对数,其中1<=i,j<=10000000

对于求解gcd(i,j)=k的对数利用莫比乌斯反演我们已经可以求解,所以最单纯的想法便是枚举素数然后逐个求解,这样肯定会TLE的。

采取了网上一大牛的思路：

我们要求解的结果为 Sigma(p, Sigma(p|d, mu(d/p)*f(d)))，其中f(d)为d|gcd(i,j)的个数,显然为(n/d)*(m/d)个.

我们对以上式子进行转化可以得到等价的表示方法为: Sigma(d, f(d)*Sigma(p, mu(d/p))),如果我们能够求解出Sigma(p,mu(d/p))的话,我们便可以在O(n)的时间内求解出上述的结果.

另g(x)=Sigma(p,mu(x/p)), g(x)是具有一下性质的:

考虑质数p'，g(p'x) = sigma(p | p'x, μ(p'x/p))。

<1>当x % p' == 0，那么考虑sigma中的变量p的所有取值，它和g(x)中p是相同的。而μ(x)这个函数，如果x有平方因子的话就等于0，因此当p != p'时μ(p'x/p) = 0，

当p == p'是，μ(p'x/p) = μ(x)。所以g(p'x) = μ(x)。

<2>当x % p' != 0，同样考虑p，会发现它的取值只比g(x)中的p多出一个p'。同理按照p是否等于p'讨论，可以得到g(p'x) = -g(x) + μ(x)。

这样的话我们便可以在求解mu的时候进行g[x]的求解.

        /*
        LL ans = 0;
        FOR(i, 1, a)
        {
            ans += 1LL*(a/i)*(b/i)*g[i];
        }
        cout << ans << endl;
        */

利用上述的方法便可以求解,但是SPOJ太狠了,上面的代码都TLE.后来经大牛的提醒我们对于上述的表达式是可以分段进行处理的即的一定的区间范围内(n/d)和(m/d)是不变的.

分段处理:

<1>首先需要处理一段区间内的g(x)的值,因此我们在预处理的时候需要处理g(x)的前缀和

<2>对于当前d求解最大的d'使得floor(n/d)=floor(n/d'), d'*floor(n/d)<=n,因此最大的d'为 floor(n/floor(n/d)),这样不断地取n和m的最小值,分段成一个一个的区间进行求解

#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cassert>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <deque>
#include <list>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <algorithm>
#include <functional>
using namespace std;

//----------------------------------------------------------
#define CL(a,b)       memset(a,b,sizeof(a))
#define CLQ(q)        while(!q.empty())q.pop();
#define FOR(i,a,b)    for(int i=a;i<=b;++i)
#define FD(i,a,b)     for(int i=a;i>=b;--i)
#define FOS(i,a,b,c)  for(int i=a;i<=b;i+=c)
#define FS(i,a)       for(int i=0;a[i];++i)
#define REP(i,n)      for(int i=0;i<n;++i)
#define PR2(a,n,m)    for(int i=0;i<n;++i){for(int j=0;j<m;++j)cout<<a[i][j]<<" ";puts("");}
#define max(a,b)      ((a)>(b)?(a):(b))
#define min(a,b)      ((a)<(b)?(a):(b))
#define checkMax(a,b) {if(a<b)a=b;}
#define checkMin(a,b) {if(a>b)a=b;}
#define READ(a)       freopen(a,"r",stdin)
#define WRITE(a)      freopen(a,"w",stdout)
#define sqr(x)        ((x)*(x))
#define inf           0x3f3f3f3f
#define INF           0x3f3f3f3f3f3f3f3fLL
#define eps           1e-10
typedef long long LL;
const double pi  = acos(-1.0);
const double hpi = asin(1.0);
//-----------------------------------------------------------

const int MAXN = 10000010;
int check[MAXN];
int prime[MAXN];
int mu[MAXN];
int g[MAXN];
int sum[MAXN];

void Moblus()
{
    CL(check, 0);
    mu[1] = 1;
    int tot = 0;
    for(int i = 2; i <= 10000000; ++i)
    {
        if(!check[i])
        {
            g[i] = 1;     /
            prime[tot++] = i;
            mu[i] = -1;
        }
        for(int j = 0; j < tot; ++j)
        {
            if(i * prime[j] > 10000000) break;
            check[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                mu[i * prime[j]] = 0;
                g[i * prime[j]] = mu[i];
                break;
            }
            else
            {
                g[i*prime[j]] = mu[i] - g[i];
                mu[i * prime[j]] = -mu[i];
            }
        }
    }
    for(int i = 1; i <= 10000000; ++i)
        sum[i] = sum[i-1] + g[i];
}

int main()
{
    //READ("aa.in"); WRITE("bb.out");

    Moblus();
    int T, a, b;
    cin >> T;
    while(T--)
    {
        cin >> a >> b;
        if(a > b) swap(a, b);
        /*
        LL ans = 0;
        FOR(i, 1, a)
        {
            ans += 1LL*(a/i)*(b/i)*g[i];
        }
        cout << ans << endl;
        */
        LL ans = 0;
        for(int i = 1, last; i <= a; i = last + 1)
        {
            last = min(a/(a/i), b/(b/i));
            ans += 1LL*(a/i)*(b/i)*(sum[last]-sum[i-1]);
        }
        cout << ans << endl;

    }
    return 0;
}