首页 > 代码库 > hnu Dirichlet's Theorem

hnu Dirichlet's Theorem

 

 

  1 /*  2 求ax+b  x属于区间[L,R];范围内素数的个数。  3 a*R+b<=10^12 ;  R-L+1<=10^6  4   5 枚举,超时。  6 1.如果GCD(a,b)>1 那么a+b 2*a+b ..都会是合数。此时只有判断b是否为素数。  7   8 2.如果GCD(a,b)=1 那么就可以列式子  9   ax+b %p = 0  其中p为素数。 如果满足,那么ax+b就是合数。 10   筛选整个区间即可。 11   12 由于最大的数字是10^12,只能被sqrt(10^12)的素数整除,所以筛选10^6内的素数。 13 由于区间L,R 10^6,开一个bool hash[ ]。 14  15 ax+b % py = 0  ===>  ax+py = -b; 16  17 根据扩展欧几里得求出 最小的x,此时的x可以为0.  while(x<0) x+p; 18  19 求出最小的x,关键还要求出第一个满足在区间[ L ,R ]里的数字。 20  21 temp = L%p; 22 x = x - temp; 23  24 while(a*(L+x)+b<=p) { //关于此处等号,是一个问题 既然a*x+b 是合数,怎么会=p,加了也不会错。 25  x = x + p; 26 } 27  28 这样的L+x就是区间[L ,R]里的第一个满足的数字。 29 而且x可以为0,刚好用hash的时候,直接对x进行哈希。 30  31 while(x<(R-L+1)){//不能等于,从0 --R-L  有 R-L+1个了。 32 hash[x] = false; 33 x = x+p; 34 } 35  36 3.最后求出结果。扫一遍哈希。 37  38  39 需要注意的是,由于a*x+b <=2的情况,所以对x==0 || x<=1 进行特判。 40 */ 41 #include<iostream> 42 #include<stdio.h> 43 #include<cstring> 44 #include<cstdlib> 45 #include<math.h> 46 using namespace std; 47 typedef long long LL; 48  49 const int maxn = 1e6+1; 50 int prime[maxn],len = 0; 51 bool s [maxn]; 52 bool hash1[maxn]; 53 void init() 54 { 55     int i,j; 56     memset(s,true,sizeof(s)); 57     for(i=2;i<maxn;i++) 58     { 59         if(s[i]==false) continue; 60         prime[++len] = i; 61         for(j=i*2;j<maxn;j=j+i) 62             s[j]=false; 63     } 64     s[0] = s[1] = false; 65 } 66 bool isprime(LL n) 67 { 68     LL i,ans; 69     if(n<maxn) return s[n]; 70     ans = (LL)sqrt(n*1.0); 71  72     for(i=1; i<=len && prime[i]<=ans; i++) 73     { 74         if(n%prime[i]==0) return false; 75     } 76     return true; 77 } 78 LL Ex_GCD(LL a,LL b,LL &x,LL& y) 79 { 80     if(b==0) 81     { 82         x=1; 83         y=0; 84         return a; 85     } 86     LL g=Ex_GCD(b,a%b,x,y); 87     LL hxl=x-(a/b)*y; 88     x=y; 89     y=hxl; 90     return g; 91 } 92 int main() 93 { 94     LL a,b,L,U,x,y; 95     LL i,p; 96     int t = 0; 97     init(); 98     while(scanf("%I64d",&a)>0) 99     {100         if(a==0)break;101         scanf("%I64d%I64d%I64d",&b,&L,&U);102         LL g = Ex_GCD(a,b,x,y);103         if(g>1)104         {105             if(L==0 && isprime(b))106                 printf("Case %d: 1\n",++t);107             else printf("Case %d: 0\n",++t);108         }109         else if(g==1)/** gcd(a,b) == 1   **/110         {111             memset(hash1,true,sizeof(hash1));112             if(L==0)113                 hash1[0] = isprime(b);114             if(L<=1)115                 hash1[1-L] = isprime(a+b);116             LL length = U-L+1;117             LL MAX = a*U+b;118             for(i=1; i<=len; i++)119             {120                 p = prime[i];121                 if(a%p==0)continue;122                 if(p*p>MAX)break;;123 124                 g = Ex_GCD(a,p,x,y);// ax+py = -b;125                 x = (x*-b) % p;126                 while(x<0) x=x+p;127 128                 LL temp = L%p;129                 x = x - temp;130                 while(x<0) x=x+p;131 132                 while(a*(x+L)+b<=p)133                 {134                     x = x+p;135                 }136                 while(x<length)137                 {138                     hash1[x]=false;139                     x=x+p;140                 }141             }142             LL hxl = 0;143             for(i=0; i<length; i++) if(hash1[i]==true) hxl++;144             printf("Case %d: %I64d\n",++t,hxl);145         }146     }147     return 0;148 }