账号 自动登录 找回密码 密码 欢迎注册
 搜索

# [擂台] 一起竞速吗？——计算区间[m,n]中所有使得k+1,k+3,k+7,k+9均为素数的k的个数

### 马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

×

1. c=0;for(i=1,oo,k=(2*i-1)*105;if(isprime(k-4),if(isprime(k-2),if(isprime(k+2),if(isprime(k+4),if(k>2^32,break);c+=1;print((k-5)/10))))))

1. for(t=1,32,c=1/*notice that 11,13,17,19 was not calculated*/;for(i=1,oo,k=(2*i-1)*105;if(isprime(k-4),if(isprime(k-2),if(isprime(k+2),if(isprime(k+4),if(k>2^t,break);c+=1;/*print((k-5)/10)*/)))));print([t,c-(t<=3)]))

1. 19:57:48> for(t=1,32,c=1/*notice that 11,13,17,19 was not calculated*/;for(i=1,oo,k=(2*i-1)*105;if(isprime(k-4),if(isprime(k-2),if(isprime(k+2),if(isprime(k+4),if(k>2^t,break);c+=1;/*print((k-5)/10)*/)))));print([t,c-(t<=3)]))
2. [1, 0]
3. [2, 0]
4. [3, 0]
5. [4, 1]
6. [5, 1]
7. [6, 1]
8. [7, 2]
9. [8, 2]
10. [9, 2]
11. [10, 2]
12. [11, 2]
13. [12, 4]
14. [13, 4]
15. [14, 6]
16. [15, 8]
17. [16, 11]
18. [17, 13]
19. [18, 18]
20. [19, 31]
21. [20, 53]
22. [21, 109]
23. [22, 160]
24. [23, 266]
25. [24, 424]
26. [25, 709]
27. [26, 1227]
28. [27, 2039]
29. [28, 3421]
30. [29, 5826]
31. [30, 10130]
32. [31, 17511]
33. [32, 30633]
34. cpu time = 26,537 ms, real time = 26,572 ms.

update: (naive)Rust

1. fn mod_exp(mut x: u64, mut d: u64, n: u64) -> u64 {
2.     let mut ret: u64 = 1;
3.     if n<4294967296{
4.         while d != 0 {
5.             if d & 1 == 1 {
6.                 ret = (ret*x)%n
7.             }
8.             d >>= 1;
9.             x = (x*x)%n;
10.         }
11.     }else{
12.         while d != 0 {
13.             if d & 1 == 1 {
14.                 ret = ((ret as u128*x as u128)%n as u128) as u64
15.             }
16.             d >>= 1;
17.             x = ((x as u128*x as u128)%n as u128) as u64;
18.         }
19.     }
20.     ret
21. }
22. pub fn miller_rabin(n: u64) -> bool {
23.     // we have a strict upper bound, so we can just use the witness
24.     // table of Pomerance, Selfridge & Wagstaff and Jeaschke to be as
25.     // efficient as possible, without having to fall back to
26.     // randomness.
27.     const WITNESSES: &[(u64, &[u64])] =
28.         &[(2_046, &[2]),
29.           (1_373_652, &[2, 3]),
30.           (9_080_190, &[31, 73]),
31.           (25_326_000, &[2, 3, 5]),
32.           (4_759_123_140, &[2, 7, 61]),
33.           (1_112_004_669_632, &[2, 13, 23, 1662803]),
34.           (2_152_302_898_746, &[2, 3, 5, 7, 11]),
35.           (3_474_749_660_382, &[2, 3, 5, 7, 11, 13]),
36.           (341_550_071_728_320, &[2, 3, 5, 7, 11, 13, 17]),
37.           (0xFFFF_FFFF_FFFF_FFFF, &[2, 3, 5, 7, 11, 13, 17, 19, 23])
38.          ];

39.     let mut d = n - 1;
40.     let mut s = 0;
41.     while d & 1 == 0 { d /= 2; s += 1 }

42.     let witnesses =
43.         WITNESSES.iter().find(|&&(hi, _)| hi >= n)
44.             .map(|&(_, wtnss)| wtnss).unwrap();
45.     'next_witness: for &a in witnesses.iter() {
46.         let mut power = mod_exp(a, d, n);
47.         assert!(power < n);
48.         if power == 1 || power == n - 1 { continue 'next_witness }
49.         if n<4294967296{
50.             for _r in 0..s {
51.                 power = (power*power)%n;
52.                 if power == 1 { return false }
53.                 if power == n - 1 {
54.                     continue 'next_witness
55.                 }
56.             }
57.         }else{
58.             for _r in 0..s {
59.                 power = ((power as u128*power as u128)%n as u128) as u64;
60.                 if power == 1 { return false }
61.                 if power == n - 1 {
62.                     continue 'next_witness
63.                 }
64.             }
65.         }
66.         return false
67.     }

68.     true
69. }

1. include!{"miller-rabin.rs"}
2. fn main() {
3.     let mut i=0u64.wrapping_sub(1);
4.     let mut n=1;
5.     let mut ctr=1;
6.     let mut limit=1<<n;
7.     loop{
8.         i+=2;
9.         let j=i*105;
10.         if j>limit {
11.             println!("{:?}",[n,ctr-(n<=3) as i32]);
12.             limit<<=1;
13.             n+=1;
14.             if n>34{break}else{continue}
15.         }
16.         if miller_rabin(j-4) && miller_rabin(j-2) && miller_rabin(j+2) && miller_rabin(j+4) {
17.             ctr+=1
18.         }
19.     }
20. }

update2: naive Rust的变种（向Rayon以及筛法靠拢）
1. include!{"miller-rabin.rs"}
2. //use rayon::prelude::*;
3. fn main() {
4.     let mut prepare=vec![105];
5.     let mut prepare_len=210;
6.     println!("total:{}",1+(0..(1u64<<32)/prepare_len).into_iter().fold(0,|s,i|s+{
7.         let start=prepare_len*i;
8.         prepare.iter().map(|&x|x+start).fold(0,|s,j|{
9.             if miller_rabin(j-4) && miller_rabin(j-2) && miller_rabin(j+2) && miller_rabin(j+4) {
10.                 s+1
11.             }else{
12.                 s
13.             }
14.         })
15.     }))
16. }

Rayon+筛法，为了保证计算效率，程序计数最好停止在筛的边界（而非之前随便写的那个边界）
1. include!{"miller-rabin.rs"}
2. extern crate rayon;
3. use rayon::prelude::*;
4. fn not_divide(d:u64,n:u64)->bool{
5.     (n/d)*d!=n
6. }
7. fn sieve(v:&mut Vec<u64>,len:&mut u64,prime:u64){
8.     let mut vv=v.par_iter().copied().fold(||Vec::new(),|mut s,x|{s.append(&mut (0..prime).map(|p|*len*p+x).filter(|&x|not_divide(prime,x-4)&&#172;_divide(prime,x-2)&&#172;_divide(prime,x+2)&&#172;_divide(prime,x+4)).collect());s})
9.                                      .reduce(||Vec::new(),|mut s,mut x|{s.append(&mut x);s});
10.     *len*=prime;
11.     vv.shrink_to_fit();
12.     vv.sort_unstable();
13.     println!("{:?}",vv.len());
14.     *v=vv
15. }
16. fn main() {
17.     let end = 12;
18.     let mut prepare=vec![105];
19.     let mut prepare_len=210;
20.     sieve(&mut prepare,&mut prepare_len,11);
21.     sieve(&mut prepare,&mut prepare_len,13);
22.     sieve(&mut prepare,&mut prepare_len,17);
23.     sieve(&mut prepare,&mut prepare_len,19);
24.     sieve(&mut prepare,&mut prepare_len,23);
25.     sieve(&mut prepare,&mut prepare_len,29);
26.     sieve(&mut prepare,&mut prepare_len,31);
27.     println!("total:{:?}(from 0 to {})",1+(0u64..end).into_par_iter().fold(||0,|s,i|s+{
28.         let start=prepare_len*i;
29.         prepare.iter().map(|&x|x+start).fold(0,|s,j|{
30.             if miller_rabin(j-4) && miller_rabin(j-2) && miller_rabin(j+2) && miller_rabin(j+4) {
31.                 s+1
32.             }else{
33.                 s
34.             }
35.         })
36.     }).reduce(||0,|a,b|a+b),end*prepare_len)
37. }

1. $rustc main.rs -C panic=abort -C opt-level=3 -C target-cpu=native -C codegen-units=1 -C lto -L ../target/release/deps -o main && time ./main 2. 7 3. 63 4. 819 5. 12285 6. 233415 7. 5835375 8. total:326240(from 0 to 77636318760) 9. real 0m19.608s 10. user 3m26.773s 11. sys 0m0.110s 复制代码 以及 1.$ rustc main.rs -C panic=abort -C opt-level=3 -C target-cpu=native -C codegen-units=1 -C lto -L ../target/release/deps -o main && time ./main
2. 7
3. 63
4. 819
5. 12285
6. 233415
7. 5835375
8. 157555125
9. total:5907339(from 0 to 2406725881560)

10. real        12m9.250s
11. user        126m49.848s
12. sys        0m4.234s

 你的  5835375 和 157555125 分别对应的 m、n 是多少？最好能提供一个标准的输出，比如，输出 10^8 内所有满足 k+1,k+3,k+7,k+9 都为素数的 k，以及这样的 k 的总个数，这样大家也方便核对检查。