UVA 11426 - GCD - Extreme (II) 莫比乌斯反演

最新推荐文章于 2019-08-29 16:39:45 发布

ZJM__

最新推荐文章于 2019-08-29 16:39:45 发布

阅读量255

点赞数

分类专栏：莫比乌斯反演文章标签：莫比乌斯反演数论线性筛积性函数

本文链接：https://blog.csdn.net/Loi_Seavot/article/details/90245128

版权

莫比乌斯反演专栏收录该内容

3 篇文章 0 订阅

订阅专栏

坚持写博客！

题目：求
$\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}gcd(i,j)$
把结果写成表格的形式容易发现：
$\frac{\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)-\sum_{i=1}^{n}i}{2}$
那么问题转变成了求：
$\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)$
太眼熟了，直接化简：
$\sum_{d=1}^{n}d\sum_{i=1}^{n/d}\mu (d)(\frac{n}{i*d})^2$
利用换元，令T=id
$\sum_{T=1}^{n}(\frac{n}{T})^2\sum_{d|T}d*\mu (T/d)$
注意数据范围：1e7，不能用之前的nlnn预处理，查了好多资料，基本理解了求狄利克雷卷积后的积性函数的方法：
$\sum_{d|T}d*\mu (T/d)$

只需要知道该函数在T为素数和T为p^k时为何值即可。

当T为素数，显然H(T)=p-1；

当T为p^k时，H（T）= p^k - p^k-1

这里解释一下其他情况:

当i,p互素时，H(T) = H(i)*H§;积性函数性质

当T为p^(k-1)*q * p时 H(T)=H(q)*H(p^k);积性函数性质

综上，我们处理一个low数组表示数i包含最小的素因子及其次数，比如low[12] = 4，因为12最小素因子为2，其次数为2，所以low[12] = 2^2 = 4;

然后当i%p = 0时用low[i] == i？来判断i是不是p^k ，如果是，带入推算得到的式子，如果不是就直接套性质即可。

参考：https://www.cnblogs.com/zwfymqz/p/9337898.html
代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<queue>
#include<stack>
#include<map>
#include<ctime>
#define up(i,x,y) for(int i = x;i <= y;i ++)
#define down(i,x,y) for(int i = x;i >= y;i --)
#define mem(a,b) memset((a),(b),sizeof(a))
#define mod(x) ((x)%MOD)
#define lson p<<1
#define rson p<<1|1
using namespace std;
typedef long long ll;
const int SIZE = 4000010;
const int INF = 2147483640;
const double eps = 1e-8;

inline void RD(int &x)
{
    x = 0;  char c; c = getchar();
    bool flag = 0;
    if(c == '-')    flag = 1;
    while(c < '0' || c > '9')   {if(c == '-')   {flag = 1;} c = getchar();}
    while(c >= '0' && c <= '9') x = (x << 1) + (x << 3) + c - '0',c = getchar();
}

int primes,prime[SIZE],mu[SIZE];
bool vis[SIZE];
ll H[SIZE],low[SIZE];//low[i]表示p1^a1,H是积性函数
void Init()
{
//    mu[1]=1;
    primes=0;
    H[1] = 1;
    for (ll i=2;i<=4000001;i++)
    {
        if (!vis[i])//素数
        {
            prime[primes++]=i;
        //    mu[i]=-1;
            H[i] = i-1;//素数情况
            low[i] = i;
        }  
        for (ll j=0;j<primes && i*prime[j]<=4000001;j++)
        {  
            vis[i*prime[j]]=1;  
            if(i%prime[j])//i pj互素
            {
                H[i*prime[j]] = H[i] * H[prime[j]];
                low[i*prime[j]] = prime[j];
            //    mu[i*prime[j]]=-mu[i];
            }
            else
            {
                low[i*prime[j]] = (low[i] * prime[j]);
                if(low[i] == i)
                {
                //    if(i*prime[j] == 4) printf("faq\n");
                    H[i*prime[j]] = i*(prime[j]-1);//特殊判断,p^k情况
                }
                else    H[i*prime[j]] = H[i / low[i]] * H[prime[j]*low[i]];//
            //    mu[i*prime[j]]=0;
                break;
            }
        }
    }
    for(int i = 1;i <= 4000001;i ++)
    {
    //    printf("i: %d %lld\n",i,H[i]);
        H[i] = H[i-1] + H[i];//前缀和
    }

}
ll n;
void solve()
{
    ll last,ans = 0;
    for(ll i = 1;i <= n;i = last+1)
    {
        last = n/(n/i);
        ans += (n/i)*(n/i)*(H[last] - H[i-1]);
    }
    ans -= (n)*(n+1)/2;
    ans /= 2;
    printf("%lld\n",ans);
}
int main(int argc, char const *argv[])
{
    Init();
    while(scanf("%lld",&n) && n)
    {
        solve();
    }
    return 0;
}