51nod 1575 Gcd and Lcm

题目:https://www.51nod.com/onlineJudge/questionCode.html#!problemId=1575

题意:
T 组数据,每组数据给出一个正整数 n ,求 ni=1ij=1ik=1lcm(gcd(i,j),gcd(i,k))mod232
T10,n109

题解:
观察到 jk 是无关的、可轮换的,不妨尝试将 jk 化简成相同的形式。
注意到 ji[gcd(i,j)=d]=jdid[gcd(id,jd)=1]=φ(id) ,所以我们可以对原式进行初步化简:

i=1nj=1ik=1ilcm(gcd(i,j),gcd(i,k))=i=1nj=1ik=1id1d2[gcd(i,j)=d1][gcd(i,k)=d2]lcm(d1,d2)=i=1nd1|id2|ilcm(d1,d2)j=1i[gcd(i,j)=d1]k=1i[gcd(i,k)=d2]=i=1nd1|id2|ilcm(d1,d2)φ(id1)φ(id2)

为了消除 d1d2 相互的限制,我们需要将 lcm(d1,d2) 拆开。
对于一般情况,不难得到 lcm(n,m)=d[gcd(nd,md)=1]ndmdd=d|n,d|md|nd,d|mdμ(d)d2nddmddd ,所以我们可以尝试进行代换:

i=1nd1|id2|ilcm(d1,d2)φ(id1)φ(id2)=i=1nd1|id2|id|d1,d|d2d|d1d,d|d2dμ(d)d2d1ddd2dddφ(id1dddd)φ(id2dddd)=i=1nd|id|iddμ(d)d2d1dd|iddd1ddφ(iddd1dd)d2dd|iddd2ddφ(iddd2dd)=i=1nd|id|iddμ(d)d2(d′′dd|iddd′′ddφ(iddd′′dd))2

在上面的最后一步,终于将 jk 对应的 d1d2 表示成类似的形式,并且上式可以表示成:

i=1nd|id|iddμ(d)d2(d′′dd|iddd′′ddφ(iddd′′dd))2=i=1nxyz=ixμ(y)y2(d|zdφ(zd))2=i=1n(idid2μ(idφ)2)(i)

其中 表示狄利克雷卷积。

因此所求是一个积性函数的前缀和,不妨设 f=idid2μ(idφ)2 ,可以发现 f(pk)=pk1((2k+1)φ(pk+1)+1)(k>0) ,或者是 f(pk)=p2f(pk1)+φ(pk)(2pk1)(k>0) ,可以尝试利用筛法来计算答案。

这里首先给出 f(pk) 的复杂推导,令 g=idφ,h=idid2μ ,不难得到 g(pk)=(k+1)pkkpk1(k>0)h(1)=1,h(pk)=pk(1p)(k>1) ,则有

f(pk)=pi|pkg2(pi)h(pki)=g2(pk)+i=0k1((i+1)piipi1)2)pki(1p)=g2(pk)+pki=0k1(i+1)2pi+1+(i+1)(3i+1)pii(3i+2)pi1+i2pi2=g2(pk)+pk(i=0ki2pi+i=0k1(i+1)(3i+1)pii=0k2((i+1)(3i+5)pi)+i=0k3(i+2)2pi)

由于 i2+(i+3)(3i+1)(i+1)(3i+5)+(i+2)2=0 ,故只需考虑 i=1,k2,k1,k 时的系数。现在考虑上和式外面的 g2(pk) ,则有
f(pk)[pk1]=1f(pk)[p2k2]=k2(k2)2+(k1)(3k5)(k1)(3k1)=0f(pk)[p2k1]=2k(k+1)(k1)2+k(3k2)=2k1f(pk)[p2k]=(k+1)2k2=2k+1

f(pk)=(2k+1)(p2kp2k1)+pk1=pk1((2k+1)φ(pk+1)+1) ,但是这样做非常麻烦,不妨回到原式来做如下推导
f(pk)=i=1pkj=1pklcm(gcd(i,pk),gcd(j,pk))=i=0kj=0kpmax(i,j)φ(pki)φ(pkj)=i=0kpiφ(pki)(2j=0i1φ(pkj)+φ(pki))=pk(2pk1)+i=0k1φ(pk)(2pk(pki+pki1))=(2k+2)p2k2kp2k1pki=0k1(p2kip2ki2)=(2k+1)(p2kp2k1)+pk1

相对简单许多。

由于本题模数固定,且 N=max(n) 不是很大,所以我采取了分段打表的做法。比如要计算 f(x)(x[L,R]) 的值,可以先筛出不超过 R 的质数,再枚举每个质数去计算对 f(x)(x[L,R]) 的贡献,这样做的复杂度是 O(Rx=LΩ(x)+RlnR)=O((RL+1)loglogR+RlnR) ,选取一个恰当的阈值 K ,分别计算 [1,K],[K+1,2K],,[N1KK+1,N] 的值,在本地打表装入代码,在处理询问时只需要处理一个长度不超过 K2 的区间,这样的复杂度是 O(T(KloglogN+NlnN)) 的,本地打表的复杂度和所有数字的质因子个数之和正相关。

显然 K 越小越好,但这样会增加代码的长度,需要具体问题具体分析。对于本题,所存储的值不超过 mod=232 ,用十进制表示需要 10 个字节,用十六进制表示需要 8 个字节,用八十五进制表示需要 5 个字节(ASCII 可见字符集可以支持九十三进制),用 p 进制表示需要 logpmod 个字节。
假设可以用 M 个字节存表,则 logpmodNKM ,取 K=logpmodNM 即可。

代码:(略去打表相关的代码)

#include <cmath>
#include <stdio.h>
#include <algorithm>
const int maxn = 31623, delta = 100000, maxd = 10001, maxl = delta + 1;
int tot, pr[maxn], d[maxn], tf[maxd];
int solve(int L, int R) // [L, R]
{
    if(L > R)
        return 0;
    static int _val[maxl], _rem[maxl];
    int *val = _val - L, *rem = _rem - L, lim = (int)ceil(sqrt(R));
    for(int i = L; i <= R; ++i)
    {
        val[i] = 1;
        rem[i] = i;
    }
    for(int i = 0; i < tot && pr[i] <= lim; ++i)
        for(int j = R / pr[i] * pr[i]; j >= L; j -= pr[i])
        {
            rem[j] /= pr[i];
            int cnt = 1, pk = pr[i];
            for( ; rem[j] >= maxn && rem[j] % pr[i] == 0; rem[j] /= pr[i], ++cnt, pk *= pr[i], val[j] *= pr[i]);
            if(rem[j] < maxn)
                for( ; d[rem[j]] == pr[i]; rem[j] /= pr[i], ++cnt, pk *= pr[i], val[j] *= pr[i]);
            val[j] *= (cnt << 1 | 1) * (pr[i] - 1) * pk + 1;
        }
    int ret = 0;
    for(int i = L; i <= R; ++i)
    {
        if(rem[i] > 1)
            val[i] *= 3 * (rem[i] - 1) * rem[i] + 1;
        ret += val[i];
    }
    return ret;
}
inline int solve(int n) // [1, n]
{
    int idx = n / delta, low = idx * delta, upp = low + delta;
    return n - low <= upp - n ? tf[idx] + solve(low + 1, n) : tf[idx + 1] - solve(n + 1, upp);
}
int main()
{
    for(int i = 2; i < maxn; ++i)
    {
        if(!d[i])
            pr[tot++] = d[i] = i;
        for(int j = 0, k; (k = i * pr[j]) < maxn; ++j)
        {
            d[k] = pr[j];
            if(d[i] == pr[j])
                break;
        }
    }
    // parse table data from string to array, satisfied that tf[i] = solve(1, i * delta).
    int t;
    scanf("%d", &t);
    while(t--)
    {
        int n;
        scanf("%d", &n);
        printf("%u\n", solve(n));
    }
    return 0;
}
©️2020 CSDN 皮肤主题: 技术黑板 设计师: CSDN官方博客 返回首页
实付0元
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值