【SPOJ-GCDEX】GCD Extreme(欧拉函数)

题目:

SPOJ-GCDEX (洛谷 Remote Judge)

分析:

求:

\[\sum_{i=1}^{n}\sum_{j=i+1}^{n}gcd(i,j)\]

这道题给同届新生讲过,由于种种原因只讲了 \(O(n)\) 预处理欧拉函数 \(O(n)\) 查询的暴力做法,顺带提了一句 “这题能根号查询” 被教练嘴了 QAQ 。以及小恐龙给我说有 \(O(n\log n)\) 预处理 \(O(1)\) 查询的另一种写法。

重点是前几天某学长讲课讲这道题,才知道有 \(O(n)\) 预处理 \(O(1)\) 查询的神仙做法,并且据说因为 SPOJ 上时限 0.237s ,复杂度比这个高的都过不去(只交了最后一种,所以我也不知道是真的假的……)。下面分别介绍这三种做法(在无特殊说明的情况下,本文中所有除法均为向下取整):

方法一

为了方便计算,先给和式中添上 \(\sum_{i=1}^{n}gcd(i,i)=\frac{n(n+1)}{2}\),再换一下枚举顺序,得到:

\[\sum_{i=1}^{n}\sum_{j=1}^{i}gcd(i,j)-\frac{n(n+1)}{2}\]

然后先枚举 \(gcd(i,j)\) 的值,并把 \(i\) 变成 \(i\cdot d\)\(j\) 变成 \(j\cdot d\)

\[\sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{i}{d}}[gcd(i,j)=1]-\frac{n(n+1)}{2}\]

(方括号表示其中式子成立时值为 \(1\) ,否则为 \(0\)

根据欧拉函数的定义( \(\varphi(x)\) 表示不超过 \(x\) 的与 \(x\) 互质的正整数的个数),得到:

\[\sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d}}\varphi(\frac{i}{d})-\frac{n(n+1)}{2}=\sum_{d=1}^{n}d\cdot S_\varphi(\frac{n}{d})-\frac{n(n+1)}{2}\]

(其中 \(S_\varphi\) 表示 \(\varphi\) 的前缀和,即 \(S_\varphi(x)=\sum_{i=1}^{x}\varphi(x)\)

用线性筛处理出欧拉函数前缀和后数论分块即可。预处理 \(O(n)\) ,查询 \(O(\sqrt{n})\)

方法二

设:

\[F(i)=\sum_{j=1}^{i}gcd(i,j)\]

则答案为:

\[\sum_{i=1}^{n}F(i)-\frac{n(n+1)}{2}\]

考虑如何计算 \(F(i)\) 。同样先枚举 \(gcd(i,j)\) 的值:

\[F(i)=\sum_{d|i}d\sum_{j=1}^{\frac{i}{d}}[gcd(i,j)=1]=\sum_{d|i}d\varphi(\frac{i}{d})\]

于是先线性筛 \(\varphi\) ,然后枚举 \(d\) ,每个 \(d\) 向所有 \(F(i\cdot d) (0<i\leq \frac{n}{d})\) 贡献 \(d\cdot \varphi(i)\) 。最后处理 \(F\) 函数的前缀和。预处理时间复杂度为 \(n\) 乘上 \(\sum_{i=1}^{n}\frac{1}{i}\),约为 \(O(n\log n)\) 。可以 \(O(1)\) 查询。

方法三:

接着方法二:

\[F(n)=\sum_{d|n}d\varphi(\frac{n}{d})=\sum_{i\cdot d=n}d\cdot \varphi(i)\]

则答案为:

\[\sum_{i=1}^{n}F(i)-\frac{n(n+1)}{2}\]

考虑如果已经求出了 \(F(n)\) ,则有(其中 \(p\) 是质数且与 \(n\) 互质):

\[F(np)=pF(n)\]

对于和式的每一项,这个 \(p\) 要么乘到 \(d\) 上,要么乘到 \(\varphi(i)\) 上,而欧拉函数是积性函数所以 \(\varphi(ip)=p\cdot \varphi(i)\)

这个性质很好,如果能快速计算 \(n=p^k\) (\(p\) 是质数)的情况就能愉快地线性筛了。(为什么?请看 【知识总结】线性筛_杜教筛_Min25筛 的线性筛部分)

考虑当 \(n=p^k\)\(F(n)=\sum_{i=0}^{k}p^{k-i}\varphi(p^i)\) 。根据欧拉函数的公式:

\[F(n)=n+\sum_{i=1}^{k}p^{k-i}p^i\cdot \frac{p-1}{p}=kp^{k-1}(p-1)+p^k\]

(注意当 \(i=0\)\(p^0=1\) 中没有因子 \(p\) ,要特殊处理)

所以当 \(k>1\)\(F(p^{k-1})=(k-1)p^{k-2}(p-1)+p^{k-1}\) ,给它乘上 \(p\) 再加上 \(p^{k-1}(p-1)\) 就是 \(F(p^k)\)\(k=1\) (即 \(n\) 是质数)的情况直接手算一下,等于 \(2n-1\)

于是可以线性筛出 \(F\) 函数然后处理前缀和, \(O(n)\) 预处理 \(O(1)\) 查询。

代码:

(方法三)

#include <cstdio>
#include <cctype>
#include <cstring>
#include <algorithm>
using namespace std;

namespace zyt
{
    template<typename T>
    inline bool read(T &x)
    {
        char c;
        bool f = false;
        x = 0;
        do
            c = getchar();
        while (c != EOF && c != '-' && !isdigit(c));
        if (c == EOF)
            return false;
        if (c == '-')
            f = true, c = getchar();
        do
            x = x * 10 + c - '0', c = getchar();
        while (isdigit(c));
        if (f)
            x = -x;
        return true;
    }
    template<typename T>
    inline void write(T x)
    {
        static char buf[20];
        char *pos = buf;
        if (x < 0)
            putchar('-'), x = -x;
        do
            *pos++ = x % 10 + '0';
        while (x /= 10);
        while (pos > buf)
            putchar(*--pos);
    }
    typedef long long ll;
    const int N = 1e6 + 10;
    ll f[N];
    int prime[N], last[N], cnt;
    bool mark[N];
    void init()
    {
        f[1] = 1;
        for (int i = 2; i < N; i++)
        {
            if (!mark[i])
                prime[cnt++] = last[i] = i, f[i] = i * 2 - 1;
            for (int j = 0; j < cnt && (ll)i * prime[j] < N; j++)
            {
                int k = i * prime[j];
                mark[k] = true;
                if (i % prime[j] == 0)
                {
                    last[k] = last[i] * prime[j];
                    if (last[k] == k)
                        f[k] = f[i] * prime[j] + (ll)i * (prime[j] - 1);
                    else
                        f[k] = f[i / last[i]] * f[last[k]];
                }
                else
                {
                    last[k] = prime[j];
                    f[k] = f[i] * f[prime[j]];
                }
            }
        }
        for (int i = 1; i < N; i++)
            f[i] += f[i - 1];
    }
    int work()
    {
        int n;
        init();
        while (read(n) && n)
            write(f[n] - (ll)(n + 1) * n / 2), putchar('\n');
        return 0;
    }
}
int main()
{
    return zyt::work();
}

转载于:https://www.cnblogs.com/zyt1253679098/p/10615665.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值