单位根反演

发表于 2018-08-13 | 分类于学习笔记 |

引子

假设现在我们要求一个式子

$\sum _ {i=0}^{n} {n \choose 2i}$

不妨构造二项式

$(1+x)^n=\sum _ {i=0}^{n} {n \choose i} x^i$

我们现在将$x=1$,$x=-1$ 带入，则有

$\begin {align} 2^n &= \sum _ {i=0}^{n} {n \choose i} \\ 0 &= \sum _ {i=0}^{n} (-1)^i {n \choose i} \end {align}$

两式相加，则有

$2^n=2*\sum _ {i=0}^{\lfloor \frac{n}{2} \rfloor} {n \choose 2i}$

问题已经解决，$Ans=2^{2*n-1}$。

那如果我们现在想计算

$\sum _ {i=0} ^{n} {n \choose ki}$

该怎么办？

简介

单位根反演就是做以上工作的东西。

简而言之，如果能快速求出某个多项式的点值，那我们就可能能求出$x$的特定倍数幂的系数和。

形式化的，我们可以求

$\sum _ {i=0}^{\lfloor \frac{n}{k} \rfloor} [x^{ik}] f(x)$

知识储备

单位根

根据$FFT$的前置知识，我们可以知道，对于一个$n$，若有$x^n-1=0$，则称$x$是$n$次单位根，记为$\omega _ n$。

单位根在复平面上的分布将单位圆均分为$n$份，在第一象限的最靠近实轴的单位根记为$\omega n^1$，一般也可记为$\omega {n}$。其余的按逆时针排序，分别为$\omega _ n^{2 \cdots n-1}$

所以我们不难发现如下几个性质：

$ \omega n^i +\omega n^{i+n/2} (2|n)$
$\omega n^i=\omega n^{i \% n}$
$\omega {kn}^{ki}=\omega {n}^{i}$

由于性质3的存在，我们可以定义本原单位根。

本原单位根的定义是：若$\omega,n$满足$\omega^n-1=0$，且$n$是最小满足条件的，则称$\omega$是$n$的本原单位根。

我们不难发现，对于单位根$\omega {n}^{i}$，若$\gcd(n,i)==1$，则$\omega n^i$一定是n的本原单位根。

否则，设$d=\gcd(i,n)$，则$\omega {n}^{i}=\omega {n’d}^{i’d}=\omega {n’}^{i’}$，$(\omega {n’}^{i’})^{n’}-1=0$。并不是$n$的本原单位根。

性质：

$\forall k,[n \mid k]=\frac{1}{n}\sum _ {i=0}^{n-1} \omega _ {n}^{ik}$

证明：

若$n \mid k$，则有

$Ans=\frac{1}{n} \sum _ {i=0}^{n-1} \omega _ {n}^{ink'}=\frac{1}{n}\sum _ {i=0}^{n-1} 1=1$

否则，这玩意是个等比数列求和，根据高中知识，我们有

$Ans=\frac{1}{n} \frac{\omega _ {n}^{0}-\omega _ {n}^{nk}}{1-\omega _ {n}^{k}}=0$

得证。

应用

现在我们考虑一开始提出的那个问题，求

$\sum _ {i=0} ^{n} {n \choose ki}$

这个就相当于在求(设$N=nk$)

$f(x)=\sum _ {i=0}^{N} {N \choose i} x^i$

的次数为$k$的倍数的项的系数和。

我们不妨设

$f(x)=\sum _ {i=0}^{N} a _ i x^i<$

那我们只需要求

$\begin {align} Ans &= \sum _ {i=0}^{N} ([x^i]f(x)) [k \mid i] \\ &= \sum _ {i=0}^{N} ([x^i]f(x)) \frac{1}{k} \sum _ {j=0}^{k-1} \omega _ {k}^{ij} \\ &= \sum _ {i=0}^{N} ([x^i]f(x)) \frac{1}{k}\sum _ {j=0}^{k-1} (\omega _ {k}^{j})^i \\ &= \sum _ {i=0}^{N} a _ i \sum _ {j=0}^{k-1} (\omega _ {k}^{j})^i\\ &= \sum _ {j=0}^{k-1} \sum _ {i=0}^{N} (\omega _ {k}^{j})^i\\ &= \frac{1}{k}\sum _ {j=0}^{k-1} f(\omega _ {k}^{j}) \end{align}$

即可。

同时注意到

$f(x)=(x+1)^n$

例题

[LOJ 6485] LJJ 学二项式定理

Link

题意

给定$n,s,a 0,a 1,a 2,a 3$，求

$[\sum _ {i=0}^{n} ({n \choose i} s^i a _ {i \bmod4}) ] \bmod 998244353$

题解

我们不妨不关注$a$数组部分，先构造二项式

$\begin {align} f(x) &= \sum _ {i=0}^{n} {n \choose i} s^i x^i *1^{n-i} \\ &= (sx+1)^n \end {align}$

现在只要分别求出$f(x)$除$4$余$0 \cdots 3$的系数之和即可，直接套用之前所讲的方法即可。

需要注意，在模$P$意义下，设$g$为$P$的原根，则$\omega _ {n}^{1}=g^{\frac{P-1}{n}}$。

第二种构造方法：

$\begin {align} f(x) &= \sum _ {i=0}^{n} {n \choose n-i} s^{n-i} x^i \\ f(x) &=(s+x)^n \end {align}$

注意，这里使用了$n-i$替换了$i$，所以要处理一下$a$数组的下标。

常系数齐次线性递推学习笔记

发表于 2018-06-28 | 分类于学习笔记 |

简介

什么是常系数齐次线性递推？

设有数列${a 0,a 1,a 2 \cdots }$满足递推关系$a n=\sum\limits {i=1}^{k}a {n-i}*f _ i$，则称该数列满足k阶齐次线性递推关系。

矩阵乘法

这个很普及的东西还是提一下吧。

假设我们有一个满足k阶齐次线性递推关系的数列${a 0,a 1,a 2 \cdots}$，它满足的齐次线性递推关系为$a n=\sum\limits {i=1}^{k}a {n-i}*f i$，现在我们要求$a n$。

假设我们现在维护着一个列矩阵，它的行数为k。

$A= \begin{bmatrix} a _ {n} \\ a _ {n-1} \\ \cdots \\ a _ {n-k+2} \\ a _ {n-k+1} \end{bmatrix}$

我们考虑如何让这个矩阵的所有元素都向前递推一格，不难设计出k行k列的转移矩阵。

$M= \begin{bmatrix} f _ 1 &f _ 2 &f _ 3 &f _ 4 & \cdots &f _ {k-2} &f _ {k-1} \\ 1 &0 &0 &0 & \cdots &0 &0 \\ 0 &1 &0 &0 & \cdots &0 &0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots\\ 0 &0 &0 &0 & \cdots &1 &0 \end{bmatrix}$

我们可以得到

$\begin{bmatrix} f _ 1 &f _ 2 &f _ 3 &f _ 4 & \cdots &f _ {k-2} &f _ {k-1} \\ 1 &0 &0 &0 & \cdots &0 &0 \\ 0 &1 &0 &0 & \cdots &0 &0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots\\ 0 &0 &0 &0 & \cdots &1 &0 \end{bmatrix} \times \begin{bmatrix} a _ {n-1} \\ a _ {n-2} \\ \cdots \\ a _ {n-k+1} \\ a _ {n-k} \end{bmatrix} = \begin{bmatrix} a _ {n} \\ a _ {n-1} \\ \cdots \\ a _ {n-k+2} \\ a _ {n-k+1} \end{bmatrix}$

现在我们要求$a _ n$，根据矩阵乘法的结合律，易得到：

$\left( \begin{bmatrix} f _ 1 &f _ 2 &f _ 3 &f _ 4 & \cdots &f _ {k-2} &f _ {k-1} \\ 1 &0 &0 &0 & \cdots &0 &0 \\ 0 &1 &0 &0 & \cdots &0 &0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots\\ 0 &0 &0 &0 & \cdots &1 &0 \end{bmatrix} \right) ^n \times \begin{bmatrix} a _ {k-1} \\ a _ {k-2} \\ \cdots \\ a _ {1} \\ a _ {0} \end{bmatrix} \times \begin{bmatrix} a _ {n+k-1} \\ a _ {n+k-2} \\ \cdots \\ a _ {n+1} \\ a _ {n} \end{bmatrix}$

所以我们只需要算出$M^n \times A$，然后取最后一个数就可以了。

可以使用矩阵快速幂，复杂度$\mathcal{O}(k^3log _ 2n)$

特征多项式

特征值，特征向量:

若有常数$\lambda$,向量$\vec v$，满足$\lambda\vec v=A\vec v$ ，则称向量$\vec v$为矩阵$A$的一组特征向量，$\lambda$为矩阵$A$的一组特征值。

秩为k的矩阵有k组线性不相关的特征向量。

特征多项式

对关系式进行变换：$(\lambda I-A)\vec v=0$。

这个等式有解的充要条件是$det(\lambda I-A)=0$,大致可以看做向量被拍扁之类的。

可以得到一个$k$次多项式，这个多项式的值等于$0$时把这个方程称为矩阵$A$的特征方程。这个多项式称为矩阵$A$的特征多项式。

特征多项式记为$f(x)=det(\lambda I-A)$。

其中，$det()$为行列式函数，$I$为单位矩阵。

$k$个解自然就是$k$个特征值。（并不需要解出来，怎么处理后面会说）

跟据算数基本定理，最终的多项式有$k$个解，可以写作$f(x)=\prod i(\lambda i-x)$。

Cayley-Hamilton定理

根据Cayley-Hamilton定理，可知$f(A)=O$，$O$为0矩阵。

下面给出一个简单证明：

$f(A)=\prod\limits {i}(\lambda {i} I - A)$

考虑将$f(A)$得到的矩阵分别乘特征向量，如果最后都得到了$0$矩阵，因为这几个特征向量线性不相关，那么可证明$f(A)$乘以任意向量都会得到$0$矩阵，从而$f(A)$为$0$矩阵。

现在问题转换为证明$f(A)$乘任意特征向量都会得到$0$矩阵。

先证明：$(\lambda i I-A) \times (\lambda j I -A) = (\lambda j I -A) \times (\lambda i I -A)$

展开即可，不再赘述。

现在考虑任意一个特征向量 $\vec v _ i$

$f(A) \times \vec v _ {i} = \prod _ {j!=i} \times (\lambda _ {j} I - A) \times (\lambda _ i I -A) \times \vec v _ {i} \\ \text{由特征向量和特征值的定义可知，} (\lambda _ i I -A) \times \vec v _ {i}=0。 \\ \therefore f(A) \times \vec v _ {i} =O$

证毕。

优化递推

设转移矩阵为$M$。

根据矩阵快速幂那套理论，我们只要求出来$M^n$就可以了。

我们考虑$M$的特征多项式$f(x)$，这是一个$k$次多项式。

我们将$M$带入，就会得到一个关于$M$的$k$次多项式。

我们可以将$M^n$分解为$M^n=f(M) \times g(M) +R(M)$。

由于$f(M)=O$，所以$M^n=R(M)$，$R(M)$是一个次数不超过$k-1$的多项式。

也就是说，我们只要求出来$M^n \% f(M)$ 就万事大吉了。

但是我们怎么求呢？

我们考虑快速幂的过程（就是倍增）。

假设我们现在已知 $g(M)=M^{2^i} \% f(M)$，现在要求$h(M)= M^{2^{i+1}} \% f(M)$。

一个直接的想法就是令$H(M)=g(M) \times g(M)$。

但是这样做，$H(M)$ 的次数是$2k-2$的。

那我们考虑原本的递推关系，$a n=\sum\limits {i=1}^{k}a {n-i}*f i$。

不难得到$M^n=\sum\limits {i=1} ^{k} M^{n-i} \times f {i}$。

所以我们可以用这个式子将多余出来的系数都向前压一位。

这样我们就得到了一个$\mathcal{O}(k^2\log _ 2 n)$的做法。

那有没有优化的余地呢？
我们从倍增的过程入手，可以发现$H(M)=g(M) \times g(M)$的过程可以由FFT加速至$\mathcal{O}(klog _ 2k)$。

现在只要解决压系数的过程就行了。

我们只要让$h(M) =H(M) \% f(M)$就行了。

等等，$f(M)$怎么求？

根据定义，$f(x)=det(xI - M)$，得到

$f(x) = |x I - M| = \begin{bmatrix} x- a _ 1 & -a _ 2 & -a _ 3 & \cdots & -a _ {k - 2} & -a _ {k - 1} & -a _ k \\ -1 & x & 0 & \cdots & 0 & 0 & 0 \\ 0 & -1 & x & \cdots & 0 & 0 & 0 \\ 0 & 0 & -1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -1 & x & 0 \\ 0 & 0& 0 & \cdots & 0 & -1 & x \end{bmatrix}$

我们对其进行第一列展开，得到

$f(x) = (x - a _ 1)M _ {11} + (-a _ 2)M _ {12} + \cdots + (-a _ k)M _ {1n} = x ^ k - a _ 1 x ^ {k - 1} - a _ 2x ^ {k - 2} - \cdots - a _ k$

$M _ {i,j}$代表$M$的算术余子式。

只要直接上多项式取模就行了。

最后的总复杂度$\mathcal{O}(k \log 2 k \log 2 n)$

我们手动模拟一下多项式取模的过程，其实可以发现我们手动向前压的过程就是在暴力取模。

附上BZOJ4161的$\mathcal{O}(k^2\log _ 2 n)$代码

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<vector>
#include<set>
#include<map>
#include<queue>
#include<stack>
#include<cassert>
typedef long long ll;
typedef unsigned long long ull;
using namespace std;

const int P=1000000007;
const int MAXN=4010;

int n,k,ans;
int f[MAXN],h[MAXN];

struct Matrix{
    int a[MAXN];
    Matrix (){memset(a,0,sizeof a);}
    int& operator [] (const int &i) {return a[i];}
    int operator [] (const int &i) const {return a[i];}
    inline Matrix operator * (const Matrix &rhs) const 
    {
        Matrix ret;
        for(int i=0;i<k;i++)
            for(int j=0;j<k;j++)
                (ret[i+j]+=1ll*a[i]*rhs[j]%P)%=P;
        for(int i=2*k-2;i>=k;ret[i--]=0) // notice ret[i--]=0
            for(int j=1;j<=k;j++)
                (ret[i-j]+=1ll*ret[i]*f[j]%P)%=P;
        return ret;
    }
}res;

Matrix ksm(Matrix a,int b) // 严格来说，不是快速幂
{
    Matrix ret;
    ret[0]=1;
    for(;b;a=a*a,b>>=1) if(b&1) ret=ret*a;
    return ret;
}

int main()
{
    scanf("%d%d",&n,&k);
    for(int i=1;i<=k;i++) scanf("%d",&f[i]),f[i]=f[i]>0?f[i]:f[i]+P;
    for(int i=0;i<k;i++) scanf("%d",&h[i]),h[i]=h[i]>0?h[i]:h[i]+P;
    if(n<k) printf("%d\n",h[n]);
    res[1]=1;ans=0;
    res=ksm(res,n);
    for(int i=0;i<k;i++)  ans=(ans+1ll*res[i]*h[i]%P)%P;
    printf("%d\n",ans);
    #ifdef LOCAL
        system("pause");
    #endif
}

杜教筛初步

发表于 2018-06-23 | 分类于学习笔记 |

之前学杜教筛的时候感觉没有学的很明白，现在重新学习一遍。

杜教筛简介

何为杜教筛？
说直白一点就是在低于线性的时间内求出某个积性函数的前缀和。

假设现在我们有一个积性函数$f(i)$，设$S(i)=\sum\limits_{j=1}^{i} f(i)$。

我们要求$S(n)$，假设我们找到了另一个任意数论函数$g(i)$ 。

那么可以很轻松地得到如下的推导：

$\begin{align} \sum_{i=1}^n (f*g) (i) &= \sum _{i=1} ^{n} \sum _{d|i} g(d) f( \frac {i} {d } ) \\ &=\sum _{d=1} ^{n} g(d) \sum _{d|i} ^{n} f ( \frac {i} {d} ) \\ &=\sum_{d=1} ^{n} g(d) \sum _{i=1} ^{ \lfloor \frac {n} {d} \rfloor } f(i) \\ &=\sum _{d=1} ^{n} g(d) S( \lfloor \frac {n} {d} \rfloor) \end{align}$

然后我们做进一步推导：

$\begin {align} &\begin {split} \sum _{i=1} ^{n} (f*g) (i) &= \sum _{d=1} ^{n} g(d) S( \lfloor \frac {n} {d} \rfloor) \\ &= g(1)S(n) + \sum _{d=2} ^{n} g(d) S( \lfloor \frac {n} {d} \rfloor ) \\ \end {split} \\ &\therefore g(1)S(n) = \sum _{i=1} ^{n} (f*g) (i) -\sum _{d=2} ^{n} g(d) S( \lfloor \frac {n} {d} \rfloor ) \end {align}$

也就是说，如果我们能够找到一个积性函数$g$，并且我们能够快速计算$(f*g)$的前缀和，$g$的前缀和，那么我们的问题就可以用我们最后得到的那个式子进行计算，具体来说就是可以用数论分块，然后最后一项进行递归处理。

例题

BZOJ 3944

Link

题意

求$\sum\limits{i=1}^{n} \varphi(i)$和$\sum\limits{i=1}^{n} \mu(i)$。

$n\leqslant 2^{31}-1$

题解

我们首先考虑$\sum\limits_{i=1}^{n} \varphi(i)$。

按照刚才那套理论，我们需要找到一个函数$g$，并且$(g*\varphi)$和$g$的前缀和均能够容易求出。

我们不难想到，$\varphi *I=id$，并且$id$和$1$的函数均可以很容易地求出。

所以不妨令$g$函数就是恒等函数$I$。

那么我们就有：

$\begin {align} S(n) &= \sum _{i=1} ^{n} i - \sum_{i=2} ^{n} S( \lfloor \frac {n} {i} \rfloor) \end {align}$

然后就可以前半部分用公式，后半部分数论分块+递归。

我们再来考虑$\sum\limits_{i=1}^{n} \mu(i)$。

按照刚才那套理论，我们需要找到一个函数$g$，并且$(g*\mu)$和$g$的前缀和均能够容易求出。

我们不难想到$\mu *1=\varepsilon$。

那么我们不妨令$g$为恒等函数$I$。

于是我们带入得到：

$S(n) = 1-\sum _{i=2} ^{n} S( \lfloor \frac {n} {i} \rfloor)$

直接数论分块+递归即可。

洛谷P3768 简单的数学题

Link

题意

求$\sum\limits{i=1} ^{n}\sum\limits\limits {j=1} ^{n} ij*gcd(i,j)$

$n \leqslant 10^{10}$

题解

我们首先对原式进行反演。

$\begin {align} &\begin {split} Ans &= \sum _{i=1} ^{n} \sum _{j=1} ^{n} ij*gcd(i,j) \\ &= \sum _{i=1} ^{n} \sum _{j=1} ^{n} ij \sum _{d|i,d|j} \varphi(d) \\ &= \sum _{d=1} ^{n} \varphi(d)\sum _{i|d} ^{n} \sum _{j|d} ^{n} ij \\ &= \sum _{d=1} ^{n} \varphi(d) d^2 \sum_{i=1} ^{\lfloor \frac {n} {d} \rfloor} \sum _{j=1} ^{ \lfloor \frac {n} {d} \rfloor} ij\\ &= \sum _{d=1} ^{n} \varphi(d) d^2 (\sum_{i=1} ^{\lfloor \frac {n} {d} \rfloor} i)^2 \\ &= \sum _{d=1} ^{n} \varphi(d) d^2 \sum_{i=1} ^{\lfloor \frac {n} {d} \rfloor} i ^3 \end {split} \\ \end {align}$

现在的问题就是快速求出$\varphi(d)d^2$的前缀和。

我们设$f(i)=\varphi(i)*i^2$，$S(i)=\sum\limits_{j=1}^{i} f(j)$。

现在我们要找到函数$g$，使得我们能够快速计算$(f*g)$的前缀和，$g$的前缀和。

也就是说，我们需要能够快速求以下这个式子：

$\begin {align} &\begin {split} Ans &= \sum _{i=1} ^{n} \sum _{d|i} f(d)*g(\frac{i}{d}) \\ &= \sum _{i=1} ^{n} \sum _{d|i} \varphi(d)d^2*g(\frac{i}{d}) \end{split} \\ &\text{我们不妨令}g(i)=i^2 \\ &\begin {split} Ans &=\sum _{i=1} ^{n} \sum _{d|i} \varphi(d) i^2 \\ &=\sum _{i=1} ^{n}i^3 \end {split} \end{align}$

所以我们只需要使$g(i)=i^2$，然后继续使用刚才的刚才的方法就可以完成此题。

[51NOD 1594] Gcd and Phi

发表于 2018-05-17 | 分类于题解 |

题目大意：给定n，求$F(n) = \sum{i=1}^{n} \sum{j=1}^{n} \varphi( \gcd(\varphi(i),\varphi(j)))$ ，多组数据。

$n \leqslant 10^5 ,T \leqslant 5$

$F(n) = \sum_{i=1}^{n} \sum_{j=1}^{n} \varphi( \gcd(\varphi(i),\varphi(j)))$

不妨设$p[x]=\sum\limits_{i=1}^{n} [\varphi(i)==x]$

不妨枚举$\varphi$的取值，易作如下变形

$\begin{align} F(n) &= \sum_{i=1}^{n} \sum_{j=1}^{n} \varphi( \gcd(\varphi(i),\varphi(j)) ) \\ F(n) &= \sum_{i=1}^{n} \sum_{j=1}^{n} \varphi(\gcd(i,j)) \times p[i] \times p[j] \\ F(n) &= \sum_{i=1}^{n} \sum_{j=1}^{n} p[i] \times p[j] \times \sum_{d|i,d|j}\varphi(d) \times \varepsilon(\gcd(\frac{i}{d},\frac{j}{d})) \\ F(n) &=\sum_{d=1}^{n} \sum_{i=1}^{\lfloor{\frac{n}{d}}\rfloor} \sum_{j=1}^{\lfloor\frac{n}{d}\rfloor} \varphi(d) \times p[i d] \times p[j d] \times \varepsilon(\gcd(i,j)) \\ F(n) &=\sum_{d=1}^{n} \sum_{i=1}^{\lfloor{\frac{n}{d}}\rfloor} \sum_{j=1}^{\lfloor\frac{n}{d}\rfloor} \varphi(d) \times p[i d] \times p[j d] \times \sum_{d'|i,d'|j} \mu(d') \\ F(n) &=\sum_{d=1}^{n} \sum_{d'=1} ^{\lfloor\frac{n}{d}\rfloor}\sum_{i=1}^{\lfloor{\frac{n}{dd'}}\rfloor} \sum_{j=1}^{\lfloor\frac{n}{dd'}\rfloor} \mu(d') \times \varphi(d) \times p[i d d'] \times p[j d d'] \\ F(n) &=\sum_{d=1}^{n} \varphi(d) \sum_{d'=1} ^{\lfloor\frac{n}{d}\rfloor} \mu(d') \sum_{i=1}^{\lfloor{\frac{n}{dd'}}\rfloor} \sum_{j=1}^{\lfloor\frac{n}{dd'}\rfloor} p[i d d'] \times p[j d d'] \end{align}$

不妨设$G(x)=\sum\limits{i=1}^{\lfloor \frac{n}{x} \rfloor } \sum\limits{j=1}^{\lfloor \frac{n}{x} \rfloor} p[ix] \times p[jx]$

显然，

$\begin{align} G(x) &= \sum\limits_{i=1}^{\lfloor \frac{n}{x} \rfloor } \sum\limits_{j=1}^{\lfloor \frac{n}{x} \rfloor} p[ix] \times p[jx] \\ G(x) &= \sum\limits_{i=1}^{\lfloor \frac{n}{x} \rfloor } p[ix] \sum\limits_{j=1}^{\lfloor \frac{n}{x} \rfloor} p[jx] \\ G(x) &=(\sum_{i=1}^{\lfloor \frac{n}{x}\rfloor} p[ix])^2 \end{align}$

这个东西可以$\mathcal{O}(n \ log_2 n)$求。

我们看看$F(n)$变成了什么

$\begin {align} F(n) &=\sum_{d=1}^{n} \varphi(d) \sum_{d'=1} ^{\lfloor\frac{n}{d}\rfloor} \mu(d') \sum_{i=1}^{\lfloor{\frac{n}{dd'}}\rfloor} \sum_{j=1}^{\lfloor\frac{n}{dd'}\rfloor} p[i d d'] \times p[j d d']\\ F(n) &= \sum_{d=1}^{n} \varphi(d) \sum_{d'=1} ^{\lfloor\frac{n}{d}\rfloor} \mu(d') \times G(dd') \end {align}$

这个东西也可以$\mathcal{O}(n \log_2 n)$求。

好，做完了。

代码如下

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<vector>
#include<set>
#include<map>
#include<queue>
#include<stack>
typedef long long ll;
typedef unsigned long long ull;
using namespace std;
const int MAXN=2E6+10;
int phi[MAXN],mu[MAXN],n,pri[MAXN],tot,p[MAXN],T;
ll f[MAXN],ans;
bool vis[MAXN];
void solve()
{
    scanf("%d",&n);tot=0;
    memset(vis,0,sizeof vis);
    memset(f,0,sizeof f);
    memset(p,0,sizeof p);
    phi[1]=mu[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i]) pri[++tot]=i,mu[i]=-1,phi[i]=i-1;
        for(int j=1;j<=tot&&pri[j]*i<=n;j++)
        {
            vis[pri[j]*i]=1;
            if(i%pri[j]==0) {mu[i*pri[j]]=0;phi[i*pri[j]]=phi[i]*pri[j];break;}
            else mu[i*pri[j]]=-mu[i],phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
    }
    for(int i=1;i<=n;i++) ++p[phi[i]];
    ans=0;
    for(int i=1;i<=n;i++)
        for(int j=i;j<=n;j+=i)
            f[i]+=p[j];
    for(int i=1;i<=n;i++) f[i]=f[i]*f[i];
    for(int i=1;i<=n;i++)
        if(mu[i])
            for(int j=1;i*j<=n;j++)
                ans+=mu[i]*f[i*j]*phi[j];
    printf("%lld\n",ans);
}
int main()
{
    scanf("%d",&T);
    while(T--) solve();
}

APIO2018 游记

发表于 2018-05-13 | 分类于游记 |

这篇游记是回石家庄的火车上写的，~~所以可能有口胡~~

Day0 2018.5.10

从石家庄坐火车去北京，还是一样的火车站，还是一样的地铁站。

~~话说前一天竟然将电源落在机房了，不爽~~

到宾馆的时候大概是11点左右，但是那时候房间都没有腾出来，不爽。

所以闭~~（包）~~队先带我们去吃了个饭~~（话说还真有酱爆鸡丁这种操作）~~。

吃完，回酒店，竟然还没好，差评。

没办法，只好先去试机。

学校好远。。。

先敲了个后缀自动机，略微调了一会，还好。

又写了写后缀数组和$NTT$，基本上一遍过，感觉$OK$。

话说ztb和czy都调了不短的时间，23333。

~~临走时顺便白嫖了点水~~

房间终于好了，不愧是五星级酒店，爽到

Day1 2018.5.11

~~上午一开始先来了张刺激人心的图片~~，如下

哇，竟然有秦岳，害怕

首先先讲了讲折纸与信息技术，感觉很懵逼。

接下来是秦岳的图片拼接算法和~~（游戏）~~人工智能。

感觉挺有意思的。

~~话说首师大附中的午饭还不错~~

下午首先讲的是二分，竟然有HEOI2016原题，感觉听懂了不少。

讲课的时候竟然还能开弹幕，学到了。~~（性感AP，在线IO）~~

然后是hht的玄学一般图匹配。

晚上在宾馆摸着。

明天考试好虚啊。

Day2 2018.5.12

今天考试。

考成mengbier。

没了。

Day 3 2018.5.13

今天是北大的人来讲课啊。。。

上午讲了些什么杜教筛，洲阁筛和$Min_25$筛和比较高端的数论知识。

下午去80中，讲了些后缀自动机。~~（话说你画我猜还真是好玩）~~

晚上坐火车回石家庄。

2018.6.22 @ 高铁

Zhang_RQ

In solitude, where we are least alone.

27 日志

8 分类

32 标签