莫比乌斯反演，杜教筛
BZOJ4176
#include <cstdio> #include <map> #define LL long long using namespace std; map <LL,LL> mpa; map <LL,LL> mpb; const LL mo=1e9+7; int b[10000001],phi[10000001],ss[10000001]; int miu[10000001],summiu[10000001],cnt,n; void euler(int n=1e7){ b[1]=1;miu[1]=1;summiu[1]=1; for (LL i=2;i<=n;i++){ if (b[i]==0) {ss[++cnt]=i;miu[i]=-1;} summiu[i]=summiu[i-1]+miu[i]; for (LL j=1;j<=cnt;j++){ if (i*ss[j]>n) break; b[i*ss[j]]=1; if (i%ss[j]==0) {miu[i*ss[j]]=0;break;} else miu[i*ss[j]]=miu[i]*(-1); } } } LL getsummiu(LL num){ if (num<=1e7) return(summiu[num]); if (mpb[num]!=0) return(mpa[num]);//可以不使用map，可能出现的值为n/i的根号n个值 LL res=1;//筛phi时为n*(n+1)/2 for (LL i=2,last;i<=num;i=last+1){ last=min((num/(num/i)),num); res-=(last-i+1)*getsummiu(num/i)%mo; res%=mo; } mpb[num]=1;mpa[num]=res; return(res); } LL f(LL n){ LL ret=0; for (int i=1,last;i<=n;i=last+1){ last=min(n/(n/i),n); ret+=(last-i+1)*(n/i)%mo; ret%=mo; } return(ret*ret%mo); } int main(){ scanf("%d",&n); euler(); LL ans=0; for (int d=1,last;d<=n;d=last+1){ last=min(n/(n/d),n); ans+=(getsummiu(last)-getsummiu(d-1))*f(n/d)%mo; ans%=mo; } ans%=mo;ans+=mo;ans%=mo; printf("%lld\n",ans); }
------------------------------------------------------------------------------
Codechef DEC16 BOUNCE
求$\sum_{i=1}^n\sum_{j=1}^n[j>\frac{lob1}{lob2}]\cdot[ j<\frac{upb1}{upb2}]\cdot(\frac{lcm(i,j)}{i}+\frac{lcm(i,j)}{j}-2)$
先不考虑两条直线带来的贡献，考虑化简
$\sum_g\sum_{i=1}^{\lfloor \frac{n}{g} \rfloor}\sum_{j=1}^{\lfloor \frac{n}{g} \rfloor}[gcd(i,j)=g]\cdot i+\sum_g\sum_{i=1}^{\lfloor \frac{n}{g} \rfloor}\sum_{j=1}^{\lfloor \frac{n}{g} \rfloor}[gcd(i,j)=g]\cdot j+\sum_{i=1}^n\sum_{j=1}^n-2$
三项分别考虑，以第一项为例反演。
$\sum_g\sum_d\mu(d) \cdot d\sum_{i=1}^{\lfloor \frac{n}{d*g}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{d*g}\rfloor}i$
上面式子中的前半部分就是$1\ast\mu\cdot n$
证明:$\varphi=\mu \ast n ,1\ast\mu \cdot n \ast \mu \ast n=\mu \cdot n \ast n \ast 1 \ast \mu=\varepsilon$
我们称该函数为f，则$\varphi \ast f=\varepsilon$
考虑杜教筛，要求一个函数的前缀和，可以考虑构造另一个函数，使得我们能求出另一个函数的前缀和与两个函数狄利克雷卷积的前缀和。
这里我们找到了$\varphi$
$\sum_{i=1}^n\varphi(i)\sum_{j=1}^{\lfloor \frac{n}{i} \rfloor}f(j)=1$
将i=1的项提出可得
$\sum_{i=1}^nf(i)=1-\sum_{i=2}^n\varphi(i)\sum_{j=1}^{\lfloor \frac{n}{i} \rfloor}f(j)$
递归求解即可。
对于给出的限制可以用类欧几里得算法求解。
#include <bits/stdc++.h> #define LL long long using namespace std; const int lim=1.5e7; const LL mo=1e9+7; const LL inv2=(mo+1)/2; const LL inv6=(mo+1)/6; struct data{ LL f,g,h; }; LL upb1,upb2,lob1,lob2,ans; int T,n,a[2000001],bt[lim+1],miu[lim+1],ss[lim+1],cnt; LL tab[200001],sumphi[lim+1],tabinv[200001],suminvphi[lim+1],phi[lim+1],m; LL conv[lim+1]; char st[2000001]; void updupb(LL x,LL y){ if (upb1*y>x*upb2) upb1=x,upb2=y; } void updlob(LL x,LL y){ if (lob1*y<x*lob2) lob1=x,lob2=y; } data likegcd(LL a,LL b,LL c,LL n){ data ret;ret.f=ret.g=ret.h=0; LL N=n%mo,N_N1=N*(N+1)%mo,N_N1_2N1=N%mo*(N+1)%mo*(2*N+1)%mo; if (a==0) return(ret); if (a>=c||b>=c){ LL N_N1_inv2=N_N1*inv2%mo,N_N1_2N1_inv6=N_N1_2N1*inv6%mo; data t=likegcd(a%c,b%c,c,n); ret.f+=t.f+a/c*N_N1_inv2%mo+b/c*(N+1)%mo;ret.f%=mo; ret.g+=t.g+a/c*N_N1_2N1_inv6%mo+(b/c)*N_N1_inv2%mo;ret.g%=mo; ret.h+=(a/c)*(a/c)%mo*N_N1_2N1_inv6%mo+ (b/c)*(b/c)%mo*(N+1)%mo+ (a/c)*(b/c)%mo*N_N1%mo+ t.h+ 2*(a/c)%mo*t.g%mo+ 2*(b/c)%mo*t.f%mo; ret.h%=mo; return(ret); }else{ LL m=(a*n+b)/c,M=m%mo; data t=likegcd(c,c-b-1,a,m-1); ret.f=N*M%mo-t.f;ret.f%=mo; ret.g=(N_N1%mo*M%mo-t.f-t.h)%mo*inv2%mo;ret.g%=mo; ret.h=N*M%mo*(M+1)%mo-2*t.g-2*t.f-ret.f;ret.h%=mo; return(ret); } } LL calc(LL n){ LL N=n%mo; if (lob1!=0){ LL sumi=0,sumj=0; data t=likegcd(upb1,0,upb2,n); sumi+=t.g; sumi-=(upb2+n/upb2*upb2)%mo*((n/upb2)%mo)%mo*inv2%mo; t=likegcd(lob1,0,lob2,n); sumi-=t.g; t=likegcd(lob2,0,lob1,n*lob1/lob2); sumj+=t.g; LL tupb=(n*lob1/lob2)/lob1; sumj-=(lob1+tupb%mo*lob1)%mo*((tupb)%mo)%mo*inv2%mo; t=likegcd(upb2,0,upb1,n*upb1/upb2); sumj-=t.g; LL ran=n*upb1/upb2-n*lob1/lob2,rr=n*upb1/upb2,rl=n*lob1/lob2+1; sumj+=N*((rl+rr)%mo)%mo*(ran%mo)%mo*inv2%mo; return((sumi+sumj)%mo); }else{ LL sumi=0,sumj=0; data t=likegcd(upb1,0,upb2,n); sumi+=t.g; sumi-=(upb2+n/upb2*upb2)%mo*((n/upb2)%mo)%mo*inv2%mo; t=likegcd(upb2,0,upb1,n*upb1/upb2); sumj-=t.g; LL ran=n*upb1/upb2,rr=n*upb1/upb2,rl=1; sumj+=N*((rl+rr)%mo)%mo*(ran%mo)%mo*inv2%mo; return((sumi+sumj)%mo); } } LL calcnod(LL n){ LL N=n%mo; if (lob1!=0){ LL nod=0; data t=likegcd(upb1,0,upb2,n); nod+=t.f;nod%=mo; nod-=n/upb2; t=likegcd(lob1,0,lob2,n); nod-=t.f;nod%=mo; return(nod); }else{ LL nod=0; data t=likegcd(upb1,0,upb2,n); nod+=t.f;nod%=mo; nod-=n/upb2;nod%=mo; return(nod); } } LL getphi(LL n,LL po){ if (po<=lim) return(sumphi[po]); if (tab[n/po]!=-1e9) return(tab[n/po]); LL ret=po%mo*((po+1)%mo)%mo*inv2%mo; for (LL i=2,nxt;i<=po;i=nxt+1){ nxt=po/(po/i); ret-=(nxt-i+1)%mo*getphi(n,po/i)%mo; ret%=mo; } tab[n/po]=ret; return(ret); } LL getinvphi(LL n,LL po){ if (po<=lim) return(suminvphi[po]); if (tabinv[n/po]!=-1e9) return(tabinv[n/po]); LL ret=1; for (LL i=2,nxt;i<=po;i=nxt+1){ nxt=po/(po/i); LL sm=getphi(n,nxt)-getphi(n,i-1); ret-=sm%mo*getinvphi(n,po/i)%mo; ret%=mo; } tabinv[n/po]=ret; return(ret); } void solve(LL n){ for (LL i=1,nxt;i<=n;i=nxt+1){ nxt=n/(n/i); LL sm=getinvphi(n,nxt)-getinvphi(n,i-1); ans+=sm*calc(n/i)%mo; ans%=mo; } LL t=calcnod(n); ans-=2*t%mo; ans%=mo;ans+=mo;ans%=mo; } void euler(){ bt[1]=1;miu[1]=1;phi[1]=1;conv[1]=1; for (int i=1;i<=lim;i++){ if (!bt[i]) {ss[++cnt]=i;miu[i]=-1;phi[i]=i-1;conv[i]=1-i;} for (int j=1;j<=cnt&&ss[j]*i<=lim;j++){ bt[ss[j]*i]=1; if (i%ss[j]==0){ miu[i*ss[j]]=0; phi[i*ss[j]]=phi[i]*ss[j]; conv[i*ss[j]]=conv[i]; break; }else miu[i*ss[j]]=(-1)*miu[i],phi[i*ss[j]]=(ss[j]-1)*phi[i],conv[i*ss[j]]=conv[i]*conv[ss[j]]; } } } void init(){ euler(); for (int i=1;i<=lim;i++) sumphi[i]=sumphi[i-1]+phi[i],sumphi[i]%=mo, suminvphi[i]=suminvphi[i-1]+conv[i],suminvphi[i]%=mo; } int gcd(int x,int y){ if (x%y==0) return(y);else return(gcd(y,x%y)); } int main(){ init(); scanf("%d",&T); while (T--){ scanf("%lld%d",&m,&n); scanf("%s",&st); int lstx=0,lsty=0,flag=1; for (int i=1;i<=n;i++) if (st[i-1]=='R'||st[i-1]=='L'){ a[i]=0; if (st[i-1]=='R'&&lstx) flag=0; if (st[i-1]=='L'&&!lstx) flag=0; lstx^=1; }else{ a[i]=1; if (st[i-1]=='T'&&lsty) flag=0; if (st[i-1]=='D'&&!lsty) flag=0; lsty^=1; } if (!flag){ printf("0\n"); continue; } if (a[1]==0) for (int i=1;i<=n;i++) a[i]^=1; int cnt0=0,cnt1=0; upb1=1e9,upb2=1,lob1=0,lob2=1; for (int i=1;i<=n;i++) if (a[i]){ updupb(cnt0+1,cnt1+1); cnt1++; }else{ updlob(cnt0+1,cnt1+1); cnt0++; } int gc=gcd(upb1,upb2);upb1/=gc;upb2/=gc; gc=gcd(lob1,lob2);lob1/=gc;lob2/=gc; if (upb1*lob2<=lob1*upb2) printf("0\n");else{ ans=0; for (int i=0;i<=10000;i++) tab[i]=-1e9,tabinv[i]=-1e9; solve(m); printf("%lld\n",ans); } } }
