zoukankan      html  css  js  c++  java
  • 莫比乌斯反演,杜教筛

    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
    ",ans);
      }

     ------------------------------------------------------------------------------

    Codechef DEC16 BOUNCE

    求$sum_{i=1}^nsum_{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_gsum_{i=1}^{lfloor frac{n}{g} floor}sum_{j=1}^{lfloor frac{n}{g} floor}[gcd(i,j)=g]cdot i+sum_gsum_{i=1}^{lfloor frac{n}{g} floor}sum_{j=1}^{lfloor frac{n}{g} floor}[gcd(i,j)=g]cdot j+sum_{i=1}^nsum_{j=1}^n-2$

    三项分别考虑,以第一项为例反演。

    $sum_gsum_dmu(d) cdot dsum_{i=1}^{lfloor frac{n}{d*g} floor}sum_{j=1}^{lfloor frac{n}{d*g} floor}i$

    上面式子中的前半部分就是$1astmucdot n$

    证明:$varphi=mu ast n ,1astmu 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}^nvarphi(i)sum_{j=1}^{lfloor frac{n}{i} floor}f(j)=1$

    将i=1的项提出可得

    $sum_{i=1}^nf(i)=1-sum_{i=2}^nvarphi(i)sum_{j=1}^{lfloor frac{n}{i} floor}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
    ");
              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
    ");else{
            ans=0;
            for (int i=0;i<=10000;i++) tab[i]=-1e9,tabinv[i]=-1e9;
            solve(m);
            printf("%lld
    ",ans);    
          }
        }
      } 
  • 相关阅读:
    kmp 算法
    jdk 和 cglib 的动态代理
    RestTemplate工具类
    bat脚本切换多个工程的分支
    字符串的左旋转
    输入一个正数s,打印出所有和为s的连续正数序列(至少含有两个数)。例如输入15,由于1+2+3+4+5=4+5+6=7+8=15,所以结果打印出3个连续序列1~5、4~6和7~8。
    枚举类型在JPA中的使用
    拾遗
    YAML DEMO
    kiali 1.26 anonymous策略修改为token
  • 原文地址:https://www.cnblogs.com/zhujiangning/p/7256081.html
Copyright © 2011-2022 走看看