BZOJ3944 杜教筛

(被owenowl要求加上的内容)flx 教我杜教筛,虽然他没做这道题。

传送门,戳这里

Problem

Description

题目

Input

一共T+1行
第1行为数据组数T(T<=10)
第2~T+1行每行一个 非负整数N,代表一组询问

Output

一共T行,每行两个用空格分隔的数ans1,ans2

Sample

Input
6
1
2
8
13
30
2333
Output
1 1
2 0
22 -2
58 -3
278 -3
1655470 2

Solution

杜教筛裸题。

下面是我学杜教筛的一点经验,和大家分享一下。
如果您已经精通杜教筛,请跳过下面的内容,直接看注意事项。

如果您还不会莫比乌斯反演,请先学习反演再看本文。否则,可能一头雾水。

一般求Phi函数和mu函数都用线性筛,可以在O(n)的时间内求出1~n的所有phi值和mu值。下面是线性筛的代码,主要思想就是筛质数,然后根据phi和mu的积性函数性质筛后面的数。

#include<bits/stdc++.h>
const int MAXN = 5000000;
int prime[MAXN+10];
ll mu[MAXN+10],phi[MAXN+10],cnt;
bool isprime[MAXN+10];
void init()
{
    mu[1]=1,phi[1]=1;
    for(register int i=2;i<=MAXN;i++)
    {
        if(!isprime[i]) prime[++cnt]=i,phi[i]=i-1,mu[i]=-1;
        for(register int j=1;j<=cnt;j++)
        {
            if(prime[j]*i>MAXN) break;
            isprime[i*prime[j]]=1;
            if((i%prime[j])==0) {phi[i*prime[j]]=phi[i]*prime[j],mu[i*prime[j]]=0;break;}
            else phi[i*prime[j]]=phi[i]*phi[prime[j]],mu[i*prime[j]]=-mu[i];
        }
    }
    for(register int i=1;i<=MAXN;i++) phi[i]+=phi[i-1];
    for(register int i=1;i<=MAXN;i++) mu[i]+=mu[i-1];
    return ;
}

Markdown炸了。。没法打数学公式。。

但是对于n很大求前缀和的情况(比如本题),O(n)的复杂度就不够了。这个时候就需要一种新的方法————杜教筛。

设h = f ⊗ g,设f的前缀和为S,而h的前缀和易求,那么可以用杜教筛。经过推导,可以得到:
杜教筛公式

首先考虑对μ函数的前缀和(欧拉函数同理):
构造一个g(i)=1!那么,有:

M(x)为mu[x]的前缀和:

求phi函数的时候,如何构造g(x)留给读者自行思考。

杜教筛关键在于构造g函数与f函数进行卷积,再通过杜教筛的公式递归求值。(可以是另一个杜教筛)
有几个注意事项:

1、首先,在用杜教筛之前,先用线性筛把开始较小的一些数(1~n^(2/3))的phi和mu给筛出来,再调用杜教筛,这样可以优化一下复杂度。

2、在用杜教筛枚举i求phi(x/i)时,分段枚举,对于一段x/i相同的i的区间只枚举一次,再乘以这一段的大小即可。否则可能会T掉

3、对于线性筛筛出来的值,要用数组存,不要用map,否则会多一个log,很可能T掉;

Code bzoj 3944

附上AC代码,前几十行写了读入优化,可以跳过。
本题AC需要靠运气,或者您去学洲阁筛吧。

/*********************************************
    Problem: 3944
    User: dxymaster
    Language: C++
    Result: Accepted
    Time:11032 ms
    Memory:103976 kb
**********************************************/

#pragma comment(linker, "/STACK:1024000000,1024000000")//防止爆栈
#include<bits/stdc++.h>
using namespace std;
#ifdef __attribute__
#define fast __attribute__((optimize(3), target("no-ieee-fp,arch=corei7")))
#else
#define fast
#endif
#define inline fast inline
namespace io
{
    // #define USE_FREAD_IO
#ifdef USE_FREAD_IO
#define BUF_SIZE 5000010
    char buf[BUF_SIZE]; int cur = BUF_SIZE; FILE *in = stdin;
    inline char gnc()
    {
        if (cur == BUF_SIZE) { fread(buf, BUF_SIZE, 1, in); cur = 0; }
        return buf[cur++];
    }
#else
    inline char gnc()
    {
        return (char)getchar();
    }
#endif
    template<typename T>
    inline void gi(T &dx)
    {
        dx = 0;
        int yc = gnc();
        bool nega = false;
        while (yc < '0' || yc > '9') { nega = (yc == '-' ? true : nega); yc = gnc(); }
        while (yc >= '0' && yc <= '9') { dx = (T)(dx * 10 + yc - '0'); yc = gnc(); }
        if (nega) { dx = -dx; }
    }
    void gc(char &a)
    {
        do a = gnc(); while (!isgraph(a));
    }
    void gss(char *c)
    {
        *c = gnc();
        while (!isgraph(*c)) *c = gnc();
        while (isgraph(*c)) *++c = gnc();
        *c = 0;
    }
#if __cplusplus >= 201103l || (defined(_MSVC_LANG) && _MSVC_LANG >= 201103l)
    template<typename T, typename ...Args> inline void gi(T &a, Args &...args)
    { gi(a); gi(args...); }
#else
    template<typename t1, typename t2> inline void gi(t1 &a, t2 &b) { gi(a); gi(b); }
    template<typename t1, typename t2, typename t3> inline void gi(t1 &a, t2 &b, t3 &c) { gi(a); gi(b); gi(c); }
    template<typename t1, typename t2, typename t3, typename t4> inline void gi(t1 &a, t2 &b, t3 &c, t4 &d) { gi(a); gi(b); gi(c); gi(d); }
    template<typename t1, typename t2, typename t3, typename t4, typename t5> inline void gi(t1 &a, t2 &b, t3 &c, t4 &d, t5 &e) { gi(a); gi(b); gi(c); gi(d); gi(e); }
#endif
}
using namespace io;
typedef long long ll;
int T,n;
map <ll,ll>m,p;
const int MAXN = 5000000;
int prime[MAXN+10];
ll mu[MAXN+10],phi[MAXN+10],cnt;
bool isprime[MAXN+10];
inline ll get_mu(ll x)
{
    if(x<MAXN) return mu[x];
    if(m[x])return m[x];
    ll res=1;
    for(register ll i=2,j=0;i<=x;i=j+1)
    {
        j=x/(x/i);
        res-=get_mu(x/i)*((ll)j-i+1);
        if(j==x)break;
    }
    m[x]=res;
    return res;
}
inline ll get_phi(ll x)
{   
        if(x<MAXN) return phi[x];
        if(p[x]) return p[x];
        ll res=0;
        if(x&1) res=(x+1)/2*x;
        else res=(x/2)*(x+1);
        for(register ll i=2,j=0;i<=x;i=j+1)
        {
            j=x/(x/i);
            res-=get_phi(x/i)*((ll)j-i+1);
            if(j==x)break;
        }
        p[x]=res;
        return res; 
}
inline void init()
{
    mu[1]=1,phi[1]=1;
    for(register int i=2;i<=MAXN;i++)
    {
        if(!isprime[i]) prime[++cnt]=i,phi[i]=i-1,mu[i]=-1;
        for(register int j=1;j<=cnt;j++)
        {
            if(prime[j]*i>MAXN) break;
            isprime[i*prime[j]]=1;
            if((i%prime[j])==0) {phi[i*prime[j]]=phi[i]*prime[j],mu[i*prime[j]]=0;break;}
            else phi[i*prime[j]]=phi[i]*phi[prime[j]],mu[i*prime[j]]=-mu[i];
        }
    }
    for(register int i=1;i<=MAXN;i++) phi[i]+=phi[i-1];
    for(register int i=1;i<=MAXN;i++) mu[i]+=mu[i-1];
    return ;
}
int main()
{
    init();
    gi(T);
    while(T--)
    {
        gi(n);
        printf("%lld %lld\n",get_phi(n),get_mu(n));
    }
    return 0;
}