LOJ6686 Stupid GCD（数论，欧拉函数，杜教筛）
做题重心转移到 LOJ 了。
至于为什么，如果你知道“……”的密码，就去看吧。
LOJ 上用户自创题大多数都不可做，今天看到个可做题（而且还是个水题），就来做了一发。
明显枚举立方根。（以下令 $m=\lfloor\sqrt[3]{n}\rfloor$）
$$\sum\limits_{i=1}^m\sum\limits_{j=i^3}^{\min(n,(i+1)^3-1)}\gcd(i,j)$$
由于 $i=m$ 比较特殊，我们把它拎出来：（其实就是把 $\min$ 拆开）
$$\sum\limits_{i=1}^{m-1}\sum\limits_{j=i^3}^{(i+1)^3-1}\gcd(i,j)+\sum_{j=m^3}^n\gcd(m,j)$$
现在考虑 $\sum\limits_{i=1}^n\gcd(x,i)$ 怎么求。
$$\sum\limits_{i=1}^n\sum\limits_{d|(x,i)}\varphi(d)$$
$$\sum\limits_{d|x}\varphi(d)\lfloor\frac{n}{d}\rfloor$$
那么套回去：
$$\sum\limits_{i=1}^{m-1}\sum\limits_{d|i}\varphi(d)(\lfloor\frac{(i+1)^3-1}{d}\rfloor-\lfloor\frac{i^3-1}{d}\rfloor)+\sum\limits_{d|m}\varphi(d)(\lfloor\frac{n}{d}\rfloor-\lfloor\frac{m^3-1}{d}\rfloor)$$
后面那部分可以直接枚举因数，暴力求欧拉函数。我这里预处理了前 $10^7\approx m^{\frac{2}{3}}$ 个欧拉函数，时间复杂度是远远小于 $O(m^{\frac{2}{3}})$ 的。
接下来继续推前面的：
$$\sum\limits_{i=1}^{m-1}\sum\limits_{d|i}\varphi(d)(\lfloor\frac{(i+1)^3-1}{d}\rfloor-\lfloor\frac{i^3-1}{d}\rfloor)$$
$$\sum\limits_{d=1}^{m-1}\varphi(d)\sum\limits_{i=1}^{\lfloor\frac{m-1}{d}\rfloor}(\lfloor\frac{(id+1)^3-1}{d}\rfloor-\lfloor\frac{(id)^3-1}{d}\rfloor)$$
$$\sum\limits_{d=1}^{m-1}\varphi(d)\sum\limits_{i=1}^{\lfloor\frac{m-1}{d}\rfloor}(\lfloor\frac{(id)^3+3(id)^2+3(id)}{d}\rfloor-\lfloor\frac{(id)^3-1}{d}\rfloor)$$
$$\sum\limits_{d=1}^{m-1}\varphi(d)\sum\limits_{i=1}^{\lfloor\frac{m-1}{d}\rfloor}((i^3d^2+3i^2d+3i)-(\lfloor i^3d^2-\frac{1}{d}\rfloor))$$
$$\sum\limits_{d=1}^{m-1}\varphi(d)\sum\limits_{i=1}^{\lfloor\frac{m-1}{d}\rfloor}((i^3d^2+3i^2d+3i)-(i^3d^2-1))$$
$$\sum\limits_{d=1}^{m-1}\varphi(d)\sum\limits_{i=1}^{\lfloor\frac{m-1}{d}\rfloor}(3i^2d+3i+1)$$
$$\sum\limits_{d=1}^{m-1}d\varphi(d)\sum\limits_{i=1}^{\lfloor\frac{m-1}{d}\rfloor}3i^2+\sum\limits_{d=1}^{m-1}\varphi(d)\sum\limits_{i=1}^{\lfloor\frac{m-1}{d}\rfloor}(3i+1)$$
分块杜教筛即可。
时间复杂度：由于每次分块杜教筛都跑得过 $10^{10}$，我就姑且当它是 $O(m^{\frac{2}{3}})$ 的。所以时间复杂度 $O(n^{\frac{2}{9}})$。
upd: 已经会证复杂度是 $O(m^{\frac{2}{3}})$ 了。如果用 map 是可以做到这个复杂度（忽略那个 $\log$ 啦），但如果用其他什么奇怪的记忆化方法复杂度可能退化到 $O(m^{\frac{5}{6}})$。
代码很丑。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128 lll;
const int maxn=,mod=,inv6=;
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define ROF(i,a,b) for(int i=(a);i>=(b);i--)
#define MEM(x,v) memset(x,v,sizeof(x))
inline lll read(){
char ch=getchar();lll x=,f=;
while(ch<'' || ch>'') f|=ch=='-',ch=getchar();
while(ch>='' && ch<='') x=x*+ch-'',ch=getchar();
return f?-x:x;
}
lll n,nn;
int pr[maxn],pl,phi[maxn],pre1[maxn],pre2[maxn],ans;
bool vis[maxn];
map<ll,int> p1,p2;
ll cb,cb1;
void init(int upr){
phi[]=;
FOR(i,,upr){
if(!vis[i]) phi[pr[++pl]=i]=i-;
FOR(j,,pl){
if(i*pr[j]>upr) break;
vis[i*pr[j]]=true;
if(i%pr[j]) phi[i*pr[j]]=phi[i]*(pr[j]-);
else{
phi[i*pr[j]]=phi[i]*pr[j];break;
}
}
}
pre1[]=pre2[]=;
FOR(i,,upr){
pre1[i]=(pre1[i-]+phi[i])%mod;
pre2[i]=(pre2[i-]+1ll*phi[i]*i)%mod;
}
}
int s1(ll n){
int nn=n%mod;
return 1ll*nn*(nn+)/%mod;
}
int s1(ll l,ll r){return (s1(r)-s1(l-)+mod)%mod;}
int s2(ll n){
int nn=n%mod;
return 1ll*nn*(nn+)%mod*(*nn+)%mod*inv6%mod;
}
int s2(ll l,ll r){return (s2(r)-s2(l-)+mod)%mod;}
int S1(ll n){
if(n<=) return pre1[n];
if(p1.count(n)) return p1[n];
int nn=n%mod,s=1ll*nn*(nn+)/%mod;
for(ll l=,r;l<=n;l=r+){
r=n/(n/l);
s=(s-1ll*(r-l+)%mod*S1(n/l)%mod+mod)%mod;
}
return p1[n]=s;
}
int S2(ll n){
if(n<=) return pre2[n];
if(p2.count(n)) return p2[n];
int nn=n%mod,s=1ll*nn*(nn+)%mod*(*nn+)%mod*inv6%mod;
for(ll l=,r;l<=n;l=r+){
r=n/(n/l);
s=(s-1ll*s1(l,r)*S2(n/l)%mod+mod)%mod;
}
return p2[n]=s;
}
int calc_phi(ll n){
if(n<=) return phi[n];
ll s=n;
for(ll i=;i*i<=n;i++) if(n%i==){
s=s/i*(i-);
while(n%i==) n/=i;
}
if(n>) s=s/n*(n-);
return s%mod;
}
int calc(ll x){
int s=;
for(ll i=;i*i<=x;i++) if(x%i==){
s=(s+1ll*calc_phi(i)*((n/i)%mod-(nn/i)%mod+mod)%mod)%mod;
if(i*i!=x) s=(s+1ll*calc_phi(x/i)*((n/(x/i))%mod-(nn/(x/i))%mod+mod)%mod)%mod;
}
return s;
}
int main(){
n=read();
cb=pow(n,1.0/);cb1=cb-;
nn=(lll)cb*cb*cb-;
init(min(cb,10000000ll));
for(ll l=,r;l<=cb1;l=r+){
r=cb1/(cb1/l);
ans=(ans+3ll*s2(,cb1/l)*(S2(r)-S2(l-)+mod)%mod+(3ll*s1(,cb1/l)+cb1/l)%mod*(S1(r)-S1(l-)+mod)%mod)%mod;
}
printf("%d\n",(ans+calc(cb))%mod);
}
upd：
被出题人针对着卡，结果真卡掉了……
其实是开立方根出精度问题了，改一改就好了。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128 lll;
const int maxn=,mod=,inv6=;
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define ROF(i,a,b) for(int i=(a);i>=(b);i--)
#define MEM(x,v) memset(x,v,sizeof(x))
inline lll read(){
char ch=getchar();lll x=,f=;
while(ch<'' || ch>'') f|=ch=='-',ch=getchar();
while(ch>='' && ch<='') x=x*+ch-'',ch=getchar();
return f?-x:x;
}
lll n,nn;
int pr[maxn],pl,phi[maxn],pre1[maxn],pre2[maxn],ans;
bool vis[maxn];
map<ll,int> p1,p2;
ll cb,cb1;
void init(int upr){
phi[]=;
FOR(i,,upr){
if(!vis[i]) phi[pr[++pl]=i]=i-;
FOR(j,,pl){
if(i*pr[j]>upr) break;
vis[i*pr[j]]=true;
if(i%pr[j]) phi[i*pr[j]]=phi[i]*(pr[j]-);
else{
phi[i*pr[j]]=phi[i]*pr[j];break;
}
}
}
pre1[]=pre2[]=;
FOR(i,,upr){
pre1[i]=(pre1[i-]+phi[i])%mod;
pre2[i]=(pre2[i-]+1ll*phi[i]*i)%mod;
}
}
int s1(ll n){
int nn=n%mod;
return 1ll*nn*(nn+)/%mod;
}
int s1(ll l,ll r){return (s1(r)-s1(l-)+mod)%mod;}
int s2(ll n){
int nn=n%mod;
return 1ll*nn*(nn+)%mod*(*nn+)%mod*inv6%mod;
}
int s2(ll l,ll r){return (s2(r)-s2(l-)+mod)%mod;}
int S1(ll n){
if(n<=) return pre1[n];
if(p1.count(n)) return p1[n];
int nn=n%mod,s=1ll*nn*(nn+)/%mod;
for(ll l=,r;l<=n;l=r+){
r=n/(n/l);
s=(s-1ll*(r-l+)%mod*S1(n/l)%mod+mod)%mod;
}
return p1[n]=s;
}
int S2(ll n){
if(n<=) return pre2[n];
if(p2.count(n)) return p2[n];
int nn=n%mod,s=1ll*nn*(nn+)%mod*(*nn+)%mod*inv6%mod;
for(ll l=,r;l<=n;l=r+){
r=n/(n/l);
s=(s-1ll*s1(l,r)*S2(n/l)%mod+mod)%mod;
}
return p2[n]=s;
}
int calc_phi(ll n){
if(n<=) return phi[n];
ll s=n;
for(ll i=;i*i<=n;i++) if(n%i==){
s=s/i*(i-);
while(n%i==) n/=i;
}
if(n>) s=s/n*(n-);
return s%mod;
}
int calc(ll x){
int s=;
for(ll i=;i*i<=x;i++) if(x%i==){
s=(s+1ll*calc_phi(i)*((n/i)%mod-(nn/i)%mod+mod)%mod)%mod;
if(i*i!=x) s=(s+1ll*calc_phi(x/i)*((n/(x/i))%mod-(nn/(x/i))%mod+mod)%mod)%mod;
}
return s;
}
int main(){
n=read();
cb=pow(n,1.0/);
while((lll)cb*cb*cb<=n) cb++;
while((lll)cb*cb*cb>n) cb--;
cb1=cb-;
nn=(lll)cb*cb*cb-;
init(min(cb,10000000ll));
for(ll l=,r;l<=cb1;l=r+){
r=cb1/(cb1/l);
ans=(ans+3ll*s2(,cb1/l)*(S2(r)-S2(l-)+mod)%mod+(3ll*s1(,cb1/l)+cb1/l)%mod*(S1(r)-S1(l-)+mod)%mod)%mod;
}
printf("%d\n",(ans+calc(cb))%mod);
}
最新文章
- JDBC数据库连接池技术
- LUA+resty 搭建验证码服务器
- [宽度优先搜索] HDU 1372 Knight Moves
- [iOS]ios archives 出现的是other items而不是iOS Apps的解决方案
- UNIX/Linux网络编程基础：图解TCP/IP协议栈
- careercup-递归和动态规划 9.5
- ASP.NET MVC学习系列 WebAPI初探
- cordova热更新
- (*p)++ 与 *p++ 与 ++*p 拨开一团迷雾
- Longest Common Substring（最长公共子序列）
- JavaScript设计模式之策略模式
- VS + QT 出现 LNK2001 无法解析的外部符号 QMetaObject 的问题
- CODEVS 2455 繁忙的都市 SCOI2005（洛谷 P2330）
- redis（四）
- Hadoop Yarn源码 - day2
- 跟我学SharePoint 2013视频培训课程——网站导航及页面元素（2）
- json数据在前端（javascript）和后端（php）转换
- Java非递归的方式获取目录中所有文件（包括目录）
- 27-拓扑排序-poj1094
- java初学4