首页 > 代码库 > 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 }
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。