51nod1237(EES解法，省空间)

这里用的是最近比较流行的EES方法解决的
不了解的可以看下这篇文章：
https://blog.csdn.net/Feynman1999/article/details/82874491
这种方法只要能找出f(pe)f(p^e)f(pe)的表达式，并且当e+1e+1e+1时可以O(1)O(1)O(1)维护，就可以用EES筛法去做
对于本题而言，考虑函数S(n)=∑i=1ngcd(i,n)S(n)=\sum_{i=1}^ngcd(i,n)S(n)=∑i=1ngcd(i,n)，可知若能快速求出其和函数的值记做res，则2∗res−(1+2+3+...+n)2*res-(1+2+3+...+n)2∗res−(1+2+3+...+n)即为ans
如何求res呢？
可知函数S(n)S(n)S(n)是积性函数，实际上S(n)=∑d∣Nφ(Nd)∗dS(n)=\sum_{d|N}\varphi(\frac{N}{d})*dS(n)=∑d∣Nφ(dN)∗d，可以看到他是欧拉函数和单位函数的狄利克雷卷积。那我们列出S(p)=2p−1S(p)=2p-1S(p)=2p−1,S(pe+1)=p∗S(Pe)+φ(pe+1)S(p^{e+1})=p*S(P^e)+\varphi(p^{e+1})S(pe+1)=p∗S(Pe)+φ(pe+1) ，那只要维护一个欧拉函数即可，φ(p)=p−1,φ(pe+1)=p∗φ(pe)\varphi(p)=p-1,\varphi(p^{e+1})=p*\varphi(p^{e})φ(p)=p−1,φ(pe+1)=p∗φ(pe)
维护p1,p0p^1,p^0p1,p0的前缀即可（看了EES才知道这里前缀是什么意思）
没有特别要注意的地方

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=1e9+7;
const int ni6=166666668;
const int ni2=500000004;
ll n,M;
vector<int> pre[2],hou[2],primes,pri_ni;inline int add(const int x, const int v) {return x + v >= mod ? x + v - mod : x + v;
}
inline int dec(const int x, const int v) {return x - v < 0 ? x - v + mod : x - v;
}//这里res是n/枚举的数
int dfs(ll res, int last, ll f){//最大质因子是prime[last-1] 但将1放在外面值显然一样//cout<<n<<" "<<res<<endl;int t=dec((res > M ? hou[1][n/res] : pre[1][res]),pre[1][primes[last]-1]);//cout<<t<<endl;int ret= (ll)t*f%mod;//这里需修改for(int i=last;i<(int) primes.size();++i){int p = primes[i];if((ll)p*p > res) break;ll phi=p-1;ll ftp=(2*p-1);for(ll q=p,nres=res,nf=f*ftp%mod;q*p<=res;q*=p){//nf需修改ret = add(ret,dfs(nres/=p,i+1,nf));//枚举更大的数phi = phi*p%mod;ftp = add(p*ftp%mod,phi);nf = f*ftp%mod;//继续枚举当前素数的指数ret =add(ret,nf);//指数大于1时，记上贡献}}return ret;
}inline int ff(ll x){x%=mod;return x*(x+1)%mod*ni2%mod;
}inline int fff(ll x){x%=mod;return x*(x+1)%mod*(2*x+1)%mod*ni6%mod;
}int solve(ll n){M=sqrt(n);for(int i=0;i<2;++i){pre[i].clear();pre[i].resize(M+1);hou[i].clear();hou[i].resize(M+1);}primes.clear();primes.reserve(M+1);pri_ni.clear();pri_ni.reserve(M+1);for(int i=1;i<=M;++i){pre[0][i]=i-1;hou[0][i]=(n/i-1)%mod;pre[1][i]=dec(ff(i),1);;hou[1][i]=dec(ff(n/i),1);}for(int p=2;p<=M;++p){if(pre[0][p]==pre[0][p-1]) continue;primes.push_back(p);const ll q=(ll)p*p,m=n/p;const int pnt0=pre[0][p-1],pnt1=pre[1][p-1];const int mid=M/p;const int End=min((ll)M,n/q);for(int i=1;i<=mid;++i){hou[0][i]=dec(hou[0][i],dec(hou[0][i*p],pnt0));hou[1][i]=dec(hou[1][i],dec(hou[1][i*p],pnt1)*(ll)p%mod);//hou[2][i]=dec(hou[2][i],dec(hou[2][i*p],pnt2)*q%mod);}for(int i=mid+1;i<=End;++i){hou[0][i]=dec(hou[0][i],dec(pre[0][m/i],pnt0));hou[1][i]=dec(hou[1][i],dec(pre[1][m/i],pnt1)*(ll)p%mod);//hou[2][i]=dec(hou[2][i],dec(pre[2][m/i],pnt2)*q%mod);}for(int i=M;i>=q;--i){pre[0][i]=dec(pre[0][i],dec(pre[0][i/p],pnt0));pre[1][i]=dec(pre[1][i],dec(pre[1][i/p],pnt1)*(ll)p%mod);//pre[2][i]=dec(pre[2][i],dec(pre[2][i/p],pnt2)*q%mod);}}primes.push_back(M+1);for (int i = 1; i <= M; i++) {pre[1][i] = dec(2*pre[1][i]%mod, pre[0][i]);//2*p-1hou[1][i] = dec(2*hou[1][i]%mod, hou[0][i]);}return n>1 ? add(dfs(n,0,1),1) : 1;
}int main()
{//freopen("in.txt","r",stdin);ios::sync_with_stdio(false);cin>>n;cout<<((2LL*solve(n)-ff(n))%mod+mod)%mod<<endl;return 0;
}

另外附上杜教筛的代码：

#include<bits/stdc++.h>
#define pb push_back
#define mp make_pair
#define bit(a,b) ((a>>b)&1) //from 0
#define all(x) (x).begin(),(x).end()
const int ni2=500000004;
using namespace std;
typedef long long ll;
int debug_num=0;const int mod=1e9+7;
const int maxn=1e7+10;
bool valid[maxn];
int ans[maxn/10];
ll phi[maxn];
const int maxm=(1LL<<37)/maxn;
ll help1[maxm];
bool vis[maxm];
ll tot,up,m;inline ll add(const ll x, const ll v) {//return x + v >= mod ? x + v - mod : x + v;return (x+v)%mod;
}inline ll dec(const ll x, const ll v) {//return x - v < 0 ? x - v + mod : x - v;return (x-v)%mod;
}void get_prime(int n)
{memset(valid,true,sizeof(valid));tot=0;phi[1]=1;for(int i=2;i<=n;++i){if(valid[i]){ans[++tot]=i;phi[i]=i-1;}for(int j=1;j<=tot && ans[j]*i<=n;++j){int tp=ans[j]*i;valid[tp]=false;if(i%ans[j]==0){phi[tp]=phi[i]*ans[j];break;}else{phi[tp]=phi[i]*(ans[j]-1);}}}for(int i=1;i<=n;++i){phi[i]=add(phi[i],phi[i-1]);}
}ll get_phi(ll n)
{return (n<=up) ? phi[n] : help1[m/n] ;
}inline  ff(ll x){x%=mod;return x*(x+1)%mod*ni2%mod;
}void solve1(ll n)
{int t;if(n<=up || vis[t=m/n]) return ;vis[t]=true;help1[t]=ff(n);for(ll l=2,r;l<=n;l=r+1){r=n/(n/l);solve1(n/r);help1[t]=dec(help1[t]%mod,(r-l+1)%mod*get_phi(n/r)%mod);}
}ll solve2(ll n)
{ll ans=0;for(ll l=1,r;l<=n;l=r+1){r=n/(n/l);ll tp=n/r;if(tp<=up){ans=add(ans,(ff(r)-ff(l-1))*phi[tp]%mod);}else{ans=add(ans,(ff(r)-ff(l-1))*help1[m/tp]%mod);}}return ans;
}int main()
{//freopen("in.txt","r",stdin);up=maxn-10;get_prime(up);//cout<<clock()<<endl;ll n;cin>>n;m=n;if(m>up){memset(vis,false,sizeof(vis));solve1(n);}cout<<((2*solve2(n)%mod-ff(n))%mod+mod)%mod<<endl;return 0;
}

51nod1237(EES解法，省空间)相关推荐

最大公共子序列、子串、可重叠重复子串
最长公共子序列寻找两个给定序列的子序列,该子序列在两个序列中以相同的顺序出现,但是不必要是连续的举例:X=ABCBDAB,Y=BDCABA.序列 BCA是X和Y的一个公共子序列,但不是X和Y的最长 ...
面向对象的软件工程应用浅研
来源:http://www.biyelww.com/ [摘要]随着面向对象研究的不断深入,面向对象技术的应用越来越广泛,面向对象的思想被应用到许多不同的领域.在介绍软件工程方法的基础上分析了面向对象的 ...
华为机试HJ32：密码截取
作者:翟天保Steven 版权声明:著作权归作者所有,商业转载请联系作者获得授权,非商业转载请注明出处题目描述: 该题目是一道密码加密题,密码混合在复杂字符串中,是一个对称子字符串,比如12321A ...
关于CS找实习的不完整经验
转:http://bbs.byr.cn/#!article/GoAbroad/309095 逛飞版好久了,好像很少看到版里有说关于来美之后找实习的事情.现在大多数14fall的同学应该搞定了出国前的手 ...
谷歌、苹果、亚马逊等大厂技术面试真题
技术面试题是许多顶尖科技公司面试的主要内容,其中一些难题会令许多面试者望而却步,但其实这些题是有合理的解决方法的. 7.1 准备事项多数求职者只是通读一遍问题和解法,囫囵吞枣.这好比试图单凭看问题和 ...
【Leetcode】面试题 01.08. 零矩阵（Zero Matrix）
一.题目编写一种算法,若M × N矩阵中某个元素为0,则将其所在的行与列清零. 二.示例三.解法解法-辅助空间 #时间:O(n^2) #空间:O(n^2) class Solution:def ...
Day16-20 - Python语言进阶
数据结构和算法算法:解决问题的方法和步骤评价算法的好坏:渐近时间复杂度和渐近空间复杂度. 渐近时间复杂度的大O标记: - 常量时间复杂度 - 布隆过滤器 / 哈希存储 - 对数时间复杂度 - ...
背包问题详解-动态规划
目录 01背包问题暴力法空间优化背包问题的所有情况完全背包问题状态转移方程分析空间优化方法二:记录第i件物品的数据k 方法三:转换成01背包多重背包方法一:记录第i件物品的数据k 方 ...
重庆赛区ACM热身赛-8529. Cake
[问题描述] 小 W 和小 R 同月同日生,今天是他们的生日~ 但是只有一个生日蛋糕,切成了 n 块(每块是角度为 ai 的扇形). 现在他们两人要拿走连续的若干块蛋糕(最终没有蛋糕剩余).他们想知道 ...

51nod1237(EES解法，省空间)

代码

51nod1237(EES解法，省空间)相关推荐

最新文章

热门文章