牛客练习赛#84 F 莫比乌斯反演+杜教筛+技巧+斐波那契数列和gcd的结论+矩阵快速幂

题意：

给出$n,k$，计算
$$\sum _{i_1=1}^n\sum _{i_2=1}^n...\sum _{i_n=1}^{n}gcd(f_{i_1},f_{i_2},...,f_{i_n})$$

Solution：

这题结论有点多，先给出来，后附第二条的证明，第一条难度太大

结论$1$：

$$gcd(f_{i_1},f_{i_2},....,f_{i_n})=f_{gcd(i_1,i_2,...,i_n)}$$

于是即计算
$$\sum _{i_1=1}^n\sum _{i_2=1}^n...\sum _{i_n=1}^n{f}_{gcd(i_1,i_2,...,i_n)}$$
显然枚举因子，这里因子为$gcd$项，即
$$\sum _{d=1}^n{f}_d\sum _{i_1=1}^n\sum _{i_2=1}^n...\sum _{i_n=1}^n[gcd(i_1,i_2,...,i_n)=d]$$
这是个莫比乌斯反演的经典形式了
$$\sum _{d=1}^n{f}_d\sum _{i_1=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{i_2=1}^{\lfloor \frac{n}{d}\rfloor }...\sum _{i_n=1}^{\lfloor \frac{n}{d}\rfloor }[gcd(i_1,i_2,...,i_n)=1]=\sum _{d=1}^n{f}_d\sum _{i_1=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{i_2=1}^{\lfloor \frac{n}{d}\rfloor }...\sum _{i_n=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{t|gcd(i_1,i_2,...,i_n)}\mu (t)$$
再次优先枚举因子，这里为$t$
$$\sum _{d=1}^n{f}_d\sum _{t=1}^{\lfloor \frac{n}{d}\rfloor }\mu (t)\lfloor \frac{n}{dt}{\rfloor }^k$$

到这里，我的思路为枚举$T=dt$，即计算
$$\sum _{T=1}^n\lfloor \frac{n}{T}{\rfloor }^k\sum _{d|T}\mu (d)f(\frac{T}{d})$$
第二个求和形式极其像莫比乌斯反演的变换形式，但很可惜，斐波那契数列不是一个积性函数，否则将是绝杀，只能作罢

不能化简的话，有一个技巧，第二个求和发现上界是一个下取整，把他令成一个函数
$$R(x)=\sum _{i=1}^x\mu (i)\lfloor \frac{x}{i}\rfloor$$
先把他代入原式我们再分析他的妙处，原式即
$$\sum _{d=1}^n{f}_{d}R(\lfloor \frac{n}{d}\rfloor )$$
首先，这个表达式有下取整，可以数论分块；其次，$R(x)$内部也有下取整，又可以数论分块，并且求$R(x)$需要$\mu$的前缀和，这恰好就是杜教筛！

解决了$R(x)$部分，那么求上式也需要$f$的前缀和，如何处理？

结论$2$：

$$\sum _{i=1}^n{f}_i=f_{n+2}-1$$

那么只需要求斐波那契的一个项即可，矩阵快速幂即可轻松做到

结论$2$证明：

利用累加法，设$S(n)$为斐波那契数列的前$n$项和：
$$\begin{gathered} f_3=f_2+f_1 \\ {f}_4=f_3+f_2 \\ ... \\ {f}_n=f_{n-1}+f_{n-2} \end{gathered}$$
对应位置全部相加，得到：
$$S(n)-f_1-f_2=(S(n)-f_1-f_n)+(S(n)-f_{n-1}-f_n)$$
$f_1=f_2=1$，合并得到
$$S(n)=2f_n+f_{n-1}-1$$
递推公式合并右式，就得到了
$$S(n)=f_{n+2}-1$$
同理是可以得到奇数项和，偶数项和的

// #include
#include
#include
#include
#include
#include
#include
#include
#include
using namespace std;

using ll=long long;
const int N=2e5+5,inf=0x3fffffff;
const long long INF=0x3fffffffffffffff,mod=1e9+9;

ll qm(ll x){return (x%mod+mod)%mod;}

struct matrix
{
    ll a[2][2];
    void initial(){memset(a,0,sizeof(a));}
    void build()
    {
        initial();
        for(int i=0;i<2;i++) a[i][i]=1;
    }
    ll* operator[](int i) const {
        return (ll*)&a[i][0];
    }
    matrix operator*(const matrix& x) const
    {
        matrix ret; ret.initial();
        for(int i=0;i<2;i++)
            for(int j=0;j<2;j++)
                for(int k=0;k<2;k++) ret[i][j]=(ret[i][j]+a[i][k]*x[k][j]%mod)%mod;
        return ret;
    }
}ans={{
    {1,1},
    {0,0}
}},tmp={{
    {1,1},
    {1,0}
}};

matrix qpow(matrix a,ll b)
{
    matrix ret,base=a; ret.build();
    while(b)
    {
        if(b&1) ret=ret*base;
        base=base*base;
        b>>=1;
    }
    return ret;
}

ll f(int x)
{
    if(x<=2)
    {
        if(x==0) return 0;
        return 1;
    }
    return qm((ans*qpow(tmp,x-2))[0][0]);
}

ll fpre(int x){return qm(f(x+2)-1);}

ll qpow(ll a,ll b)
{
    ll ret=1,base=a;
    while(b)
    {
        if(b&1) ret=ret*base%mod;
        base=base*base%mod;
        b>>=1;
    }
    return ret;
}

int mu[N],cnt,prime[N];
bool nt[N];
map<int,ll>map1;

void make_prime()
{
    mu[1]=1;
    for(int i=2;i<N;i++)
    {
        if(!nt[i]) prime[++cnt]=i,mu[i]=mod-1;
        for(int j=1;j<=cnt&&i*prime[j]<N;j++)
        {
            nt[i*prime[j]]=true;
            if(i%prime[j]==0) break;
            else mu[i*prime[j]]=qm(1ll*mu[i]*mu[prime[j]]);
        }
    }
    for(int i=1;i<N;i++) mu[i]=qm(1ll*mu[i]+1ll*mu[i-1]);
}

ll calcmu(int x)
{
    if(x<N) return 1ll*mu[x];
    if(map1.count(x)) return map1[x];
    ll ret=1;
    for(int l=2,r;l<=x;l=r+1)
    {   
        r=x/(x/l);
        ret=qm(ret-(r-l+1)*calcmu(x/l));
    }
    return map1[x]=ret;
}

int n,k;

ll R(int x)
{
    ll ret=0;
    for(int l=1,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        ret=qm(ret+qm(calcmu(r)-calcmu(l-1))*qpow(x/l,k));
    }
    return ret;
}

int main()
{
    make_prime();
    cin>>n>>k; ll res=0;
    for(int l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        res=qm(res+qm(fpre(r)-fpre(l-1))*R(n/l));
    }
    cout<<res;
    return 0;  
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139