[杜教筛] Codechef January Challenge 2018 #SQRGOOD Simplify the Square Root

版权声明:本文为博主原创文章,未经博主允许随意转载。 https://blog.csdn.net/u014609452/article/details/79106226

二分转化为μ2的前缀和。
然后转化为O(n13)的运算,但是需要预处理μ的前缀和,大力杜教筛求和。
然后感谢阿爷教我把二分改成了迭代,小范围内一个一个挪,用rho求μ(n),然后就能卡进去了。
复杂度似乎是萎的吧。

#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<vector>
#include<map>
#include<algorithm>
#define pb push_back
using namespace std;
typedef long long ll;

namespace Rho{
  int prime[9]={2,3,5,7,11,13,17,19,23};  
  unsigned long long RR;  
  inline long long R(long long Mo){return (RR+=4179340454199820289ll)%Mo;}  
  inline long long mul(long long x,long long y,long long Mo){  
    long long tmp=(x*y-(long long)((long double)x/Mo*y+0.5)*Mo);  
    return tmp<0?tmp+Mo:tmp;  
  }  
  inline long long Pow(long long x,long long k,long long Mo){  
    long long ans=1;  
    for(;k;k>>=1,x=mul(x,x,Mo))if(k&1)ans=mul(ans,x,Mo);  
    return ans;  
  }  
  inline bool MR(long long n){  
    if(n<=1)return false;  
    for(int i=0;i<9;i++)if(n==prime[i])return true;  
    long long d=n-1;int tmp=0;  
    while((d&1)==0)d>>=1,tmp++;  
    for(int i=0;i<9;i++){  
      long long x=Pow(prime[i],d,n),p=x;  
      for(int j=1;j<=tmp;j++){  
    x=mul(x,x,n);  
    if((x==1)&&(p!=1)&&(p!=n-1))return false;  
    p=x;  
      }  
      if(x!=1)return false;  
    }  
    return true;  
  }  
  inline long long f(long long x,long long c,long long Mo){return (mul(x,x,Mo)+c)%Mo;}  
  inline long long gcd(long long x,long long y){return x==0?y:gcd(y%x,x);}  
  inline long long check(long long c,long long n){  
    long long x=R(n),y=f(x,c,n),p=n;  
    while((x!=y)&&((p==n)||(p==1))){  
      if(x>y)p=gcd(n,x-y);  
      else p=gcd(n,y-x);  
      x=f(x,c,n);y=f(f(y,c,n),c,n);  
    }  
    return p;  
  }

  vector<ll> v;

  inline void rho(long long n){  
    if(n<=1)return;  
    if(MR(n)){v.pb(n);return;}  
    while(true){  
      long long tmp=check(R(n-1)+1,n);  
      if(tmp!=n && tmp!=1){rho(tmp),rho(n/tmp);return;}  
    }  
  }

  inline bool nmu2(ll n){
    v.clear(); rho(n);
    sort(v.begin(),v.end());
    for (int i=0;i<v.size();i++) if (v[i]==v[i+1]) return 1;
    return 0;
  }
}

const int maxn=20000000;
const int N=maxn+5;

int prime[7000000],num;
int mu[N],mu2[N]; 

const int P=10000007;

inline void Pre(int n){
  mu[1]=1; int *vst=mu2;
  for (int i=2;i<=n;i++){
    if (!vst[i]) prime[++num]=i,mu[i]=-1; ll t;
    for (int j=1;j<=num && (t=(ll)i*prime[j])<=n;j++){
      vst[t]=1;
      if (i%prime[j]==0){
    mu[t]=0; break;
      }else
    mu[t]=-mu[i];
    }
  }
  for (int i=1;i<=n;i++) mu2[i]=((bool)mu[i])+mu2[i-1],mu[i]=mu[i-1]+mu[i];
}

map<int,int> Map;  

inline int Sum(int n){  
  if (n<=maxn) return mu[n];  
  if (Map.find(n)!=Map.end()) return Map[n];  
  int tem=1; int l,r;  
  for (l=2;l*l<=n;l++) tem-=Sum(n/l);  
  for (int t=n/l;l<=n;l=r+1,t--){  
    r=n/t;  
    tem-=(r-l+1)*Sum(t);  
  }  
  return Map[n]=tem;  
}

inline ll S(ll n){
  if (n<=maxn) return n-mu2[n];
  int x=sqrt(n); ll ret=0,cur,lst=0;
  for (int i=1,j;i<=x;i=j+1){
    ll t=n/i/i; j=sqrt(n/t);
    ret+=((cur=Sum(j))-lst)*t;
    lst=cur;
  }
  return n-ret;
}

const double PI=acos(-1.0);
const double _K=1.0/(1-6.0/PI/PI);
const double K=2.550546096730430440286486962;

ll n;

inline ll Mod(ll x,ll f){
  if (f>=n){
    while (f>=n){
      f-=Rho::nmu2(x);
      if (f<n) return x;
      x--;
    }
  }else{
    while (f<n){
      f+=Rho::nmu2(x+1);
      if (f>=n) return x+1; 
      x++;
    }
  }
}

int main(){
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  Pre(maxn);
  scanf("%lld",&(n));
  ll L,R,MID,f1,f2;
  L=1; R=max(1e10,2.550546098*n);
  f1=S(L); f2=S(R);
  while (L+1<R){
    if (f2-n<=1000){
      R=Mod(R,f2);
      break;
    }
    if (n-f1<=1000){
      R=Mod(L,f1);
      break;
    }
    MID=L+(long double)(n-f1)/(f2-f1)*(R-L);
    ll f=S(MID);
    if (f==n){
      R=Mod(MID,f); break;
    }
    if (f>n)
      f2=f,R=MID;
    else
      f1=f,L=MID;
  }
  printf("%lld\n",R);
  return 0;
}
阅读更多

没有更多推荐了,返回首页