题解 P1587 【[NOI2016]循环之美】

~~不知道为什么脑子抽筋了想用迭代法写这题，然后就把自己写死了~~ ### 前置芝士：杜教筛/洲阁筛 蛤?不会杜教筛？可以出门左转你站膜板区去学习一下 不过可能你学的是递归版的杜教筛常数会大一点，这篇题解会介绍一下如何用迭代法写杜教筛 ### 前置芝士：莫比乌斯反演 为了做这道题我们可能需要一些比较初级的反演技巧和交换求和号的技巧，比如我们应该知道这些结论都是对的 这篇题解当中我们将会用$\epsilon(i,j)$表示$gcd(i,j)==1$这个式子 $$\lfloor \frac{n}{d} \rfloor$$ $$\sum_{i=1}^{n}g(i)\sum_{d|n}f(d)=\sum_{d=1}^{n}f(d)\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}g(id)$$ $$epsilon(i,j)=\sum_{d|i,d|j}\mu(d)$$ $$\epsilon(dt,j)=\epsilon(d,j)\epsilon(t,j)$$ $$\mu(dt)=\mu(d)\mu(t)\epsilon(d,t)$$ 解释一下第3个式子利用了$\mu$是一个积性函数的性质，如果$d,t$不互质的话$\mu$就是0否则就是二者相乘 并且你还需要知道对于一个$n$来讲，$\lfloor \frac{n}{d} \rfloor$恰好有$2\sqrt{n}$中取值，并且对于后$\sqrt{n}$取值来讲，如果他们的值是x，那么$\lfloor \frac{n}{x} \rfloor$的值互不相同这个结论,当然您可能也需要知道什么叫数论分块(就是如何枚举$\lfloor \frac{n}{d} \rfloor$的取值的一种技巧) 如果您碰巧不会很懂上面的芝士可以去写几道莫比乌斯反演然后来看这题 ## 本题题解 首先我们需要通过题面得出答案事实上等于这个神奇的式子 $$\sum_{i=1}^{n}\sum_{j=1}^{m}\epsilon(i,j)\epsilon(j,k)$$ 为什么呢？ 首先$epsilon(i,j)$这个限制条件是比较好理解的，但是为什么第二个限制条件是对的呢？也就是说我们的分母仅仅和k互质就能推出这个东西是一个纯循环的小数了 首先我们不难理解如果$\frac{1}{y}$是一个纯循环的小数，那么$\frac{x}{y}$也会是一个纯循环的小数，具体证明的话就是两个纯循环小数相加还是纯循环小数，所以一个分数是不是纯循环小数看起来仅仅和分母有关 那么我们假设这个纯循环小数的循环节长度为$l$我们可以用错位相减法得到下面的式子 $$\frac{k^l}{y}-\lfloor \frac{k^l}{y} \rfloor=\frac{1}{y}$$ 也就是这个小数向右挪了$l$位之后的小数部分依然和原来一样 那么把等式两边同时乘y可以得到 $$k^l-\lfloor \frac{k^l}{y} \rfloor y=1$$ 这个东西好像就是取余的定义式子吧？让我们把式子写成这样 $$ k^l \equiv 1 \mod y$$ 根据欧拉定理我们可以知道当$k,y$互质的时候我们就存在这样的$l$否则不存在 所以我们证明了一个分数是纯循环的当且仅当它的分母和进制互质 所以答案自然就是上面提到的式子了 $$\sum_{i=1}^{n}\sum_{j=1}^{m}\epsilon(i,j)\epsilon(j,k)$$ 我们发现$\epsilon(j,k)$这个东西相当的不好搞定，因此让我们先不着急展开这个式子,而是把$\epsilon(i,j)$用$\mu$展开 $$\sum_{j=1}^{m}\epsilon(j,k)\sum_{i=1}^{n}\sum_{d|i,d|j}\mu(d)$$ 交换一下求和号可以得到 $$\sum_{j=1}^{m}\epsilon(j,k)\sum_{d|j}\mu(d)\lfloor \frac{n}{d} \rfloor$$ 然后我们继续交换求和号可以得到 $$\sum_{d=1}^{n}\mu(d)\lfloor \frac{n}{d} \rfloor\sum_{j=1}^{\lfloor \frac{m}{d} \rfloor}\epsilon(jd,k)$$ 看起来可以用刚才的结论了，让我们把$\epsilon(jd,k)$拆开 $$\sum_{d=1}^{n}\mu(d)\lfloor \frac{n}{d} \rfloor\sum_{j=1}^{\lfloor \frac{m}{d} \rfloor}\epsilon(j,k) \epsilon(d,k)$$ 然后似乎可以把$\epsilon(d,k)$提到前面去 $$\sum_{d=1}^{n}\mu(d)\epsilon(d,k)\lfloor \frac{n}{d} \rfloor\sum_{j=1}^{\lfloor \frac{m}{d} \rfloor}\epsilon(j,k)$$ 我们发现函数 $$f(n)=\sum_{i=1}^{n}\epsilon(i,k)=\sum_{i=1}^{n}\epsilon(i\%k,k)=\lfloor \frac{n}{k} \rfloor f(k)+f(n\%k)$$ 似乎是可以在暴力预处理前$k$项之后快速计算的 那么我们的式子就变成了 $$\sum_{d=1}^{n}\mu(d)\epsilon(d,k)\lfloor \frac{n}{d} \rfloor f(\lfloor \frac{m}{d} \rfloor)$$ 看起来我们似乎可以对$\lfloor \frac{n}{d} \rfloor,\lfloor \frac{m}{d} \rfloor$的值进行数论分块了，当然前提是我们需要求出这个东西的前缀和 $$f(d)=\mu(d)\epsilon(d,k)$$ 怎么做呢? 我们设$S(n)$表示下面的式子，然后接下来我们需要对$\lfloor \frac{n}{d} \rfloor$和$\lfloor \frac{m}{d} \rfloor$这两个式子的一共$2\sqrt{n}+2\sqrt{m}$可能的取值处理出对应的$S$函数值 $$S(n)=\sum_{i=1}^{n}\mu(d)\epsilon(d,k)$$ 此时让我们用$\mu$把$\epsilon$展开 $$S(n)=\sum_{i=1}^{n}\mu(d)\sum_{T|d,T|k}\mu(T)$$ 交换一下求和号可以得到 $$S(n)=\sum_{T|k}\mu(T)\sum_{d=1}^{\lfloor \frac{n}{T} \rfloor}\mu(dT)$$ 根据刚才的结论我们可以得到 $$S(n)=\sum_{T|k}\mu(T)\sum_{d=1}^{\lfloor \frac{n}{T} \rfloor}\mu(d)\mu(T)\epsilon(d,T)$$ 把$\mu(T)$提出来可以得到 $$S(n)=\sum_{T|k}\mu^2(T)\sum_{d=1}^{\lfloor \frac{n}{T} \rfloor}\mu(d)\epsilon(d,T)$$ 看起来我们似乎得到了一个可以不断递归的式子？ 我们重新更改一下我们需要求的数的定义，设$G(n,k)$表示 $$\sum_{d=1}^{n}\mu(d)\epsilon(d,k)$$ 则我们可以列出递归式 $$G(n,k)=\sum_{T|k}\mu^2(T)G(\lfloor \frac{n}{d} \rfloor,T)$$ 那么我们打一下表会发现$2000$以内一个数的约数个数最大不超过$40$个 所以看起来我们可以暴力递归下去算咯？ 慢着，似乎我们还没有处理边界条件是什么…… $$G(n,1)=\sum_{d=1}^{n}\mu(d),S(n)=G(n,k)$$ 蛤？$\mu$函数的前缀和？直接杜教筛就行了 ### trick:迭代版杜教筛 众所周知杜教筛其实是一种记忆化搜索的技术 假设我们令$S(n)$表示这个式子 $$S(n)=\sum_{d=1}^{n}\mu(d)$$ 那么我们可以得到 $$\sum_{i=1}^{n}\epsilon(i)=\sum_{i=1}^{n}\sum_{d|i}\mu(\frac{i}{d})=\sum_{d=1}^{n}\sum_{T=1}^{\lfloor \frac{n}{d} \rfloor}\mu(T)$$ 也就是说 $$1=\sum_{i=1}^{n}S(\lfloor \frac{n}{d} \rfloor)$$ 稍微整理一下式子就是 $$S(n)=1-\sum_{i=2}^{n}S(\lfloor \frac{n}{d} \rfloor)$$ 那么我们发现这东西其实就是一个dp式子啊 所以作为一个dp来讲我们当然可以递推的实现它咯 不过问题来了，直接dp的话我们似乎开不下那么大的数组，怎么办呢？ 根据$\lfloor \frac{n}{d} \rfloor$的一个性质是对于后$\sqrt{n}$种值,$n$除他们的取值各不相同，因此我们可以将他们的值存储在下标为$\frac{n}{x}$的位置上 然后我们从**小到大**对于每一个$\lfloor \frac{n}{d} \rfloor$的取值大力dp一下就搞定了，当然不要忘记线性筛出来前$n^{0.67}$部分的点值，否则我们的复杂度将会退化到$O(n^{0.75})$ 这个迭代写法也证明了为什么杜教筛复杂度计算的时候只需展开一层，因为每条转移的边只会被访问一次,后面的部分不是什么高阶无穷小而就是0，~~恩，就是这么暴力~~ 同理我们也可以用这个trick去解决一下处理$G$函数的问题，我们把k的约数之间的转移边建出来之后大力dp也可以做到迭代计算了 上代码~ ```C // luogu-judger-enable-o2 #include<cstdio> #include<algorithm> #include<cmath> using namespace std;const int N=45;const int P=5*1e4+10;const int L=3*1e6+10;const int K=2010; typedef long long ll; int zhi[L];bool book[L];int ct;int sum[L];int n;int m;int k; inline int gcd(int a,int b){if(a<b)swap(a,b);while(b){int c=a%b;a=b;b=c;}return a;} struct arr//加了科技的数组 { ll v1[P];ll v2[P];int mx;int lim; inline ll& operator[](const int& x){return (x<=lim)?v1[x]:v2[mx/x];} }dp1[N],gn,dp2[N],gm;ll f[K];int hd; int v[K];int x[K];int al[N];int nu[K];int val[N];int q[P<<1];int tot; inline void add(int u,int V){v[++ct]=V;x[ct]=al[u];al[u]=ct;} inline void sev(int n,arr& dp)//迭代版杜教筛 { tot=0;for(int i=1,r,t;i<=n;i=r+1)t=n/i,q[++tot]=t,r=n/t; for(int i=tot;i>=1;i--) { int t=q[i];if(t<=L-10){dp[t]=sum[t];continue;}ll ret=1; for(int j=2,pr;j<=t;j=pr+1) pr=(t/(t/j)),ret-=(ll)(pr-j+1)*dp[t/j];dp[t]=ret; } } inline void dypr(int n,arr* dp)//迭代版dp { tot=0;for(int i=1,r,t;i<=n;i=r+1)t=n/i,q[++tot]=t,r=n/t; for(int i=2;i<=hd;i++) for(int j=al[i];j;j=x[j]) {int lc=v[j];int dv=val[lc];for(int k=tot;k>=1;k--)dp[i][q[k]]+=dp[lc][q[k]/dv];} }inline ll getf(ll x){return (x/k)*f[k]+f[x%k];} int main() { sum[1]=1; for(int i=2;i<=L-10;i++) { if(book[i]==false){zhi[++ct]=i;sum[i]=-1;} for(int j=1;j<=ct&&i*zhi[j]<=L-10;j++) { book[i*zhi[j]]=true; if(i%zhi[j]==0)sum[i*zhi[j]]=0;else sum[i*zhi[j]]=-sum[i]; } }scanf("%d%d%d",&n,&m,&k); for(int i=1;i<=k;i++)if(k%i==0&&sum[i]!=0)nu[i]=++hd,val[hd]=i; for(int i=1;i<=L-10;i++)sum[i]+=sum[i-1]; for(int i=k;i>=1;i--) if(nu[i])for(int j=i;j<=k;j+=i)if(nu[j])add(nu[j],nu[i]); for(int i=1;i<=hd;i++)dp1[i].mx=n,dp1[i].lim=sqrt(n); for(int i=1;i<=hd;i++)dp2[i].mx=m,dp2[i].lim=sqrt(m); sev(n,dp1[1]);dypr(n,dp1);gn=dp1[hd]; sev(m,dp2[1]);dypr(m,dp2);gm=dp2[hd]; for(int i=1;i<=k;i++)f[i]=f[i-1]+(gcd(k,i)==1);ll ans=0; for(int i=1,r,tp=1;i<=min(n,m);i=r+1)//大力数论分块就行了 { int t1=n/i;int t2=m/i;int v1=n/t1;int v2=m/t2; if(v1<v2)ans+=(ll)t1*getf(t2)*(gn[v1]-((tp)?gn[i-1]:gm[i-1])),tp=1,r=v1; else ans+=(ll)t1*getf(t2)*(gm[v2]-((tp)?gn[i-1]:gm[i-1])),tp=0,r=v2; }printf("%lld\n",ans);return 0;//拜拜程序~ } ```