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

.·.·.

首先发一个Pari/GP的代码

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)]))

复制代码

过一会儿写筛法&&多线程
初步计数的结果[n,小于2^n的满足k+1,k+3,k+7,k+9均为素数的k]如下：

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
首先放一个头文件,miller-rabin.rs，用witness进行素性判定

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

复制代码

程序：

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. 
   
10. real        0m19.608s
    
11. user        3m26.773s
    
12. 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. 
    
11. real        12m9.250s
    
12. user        126m49.848s
    
13. sys        0m4.234s

复制代码

不知道大家还有什么神奇优化技巧

补充内容 (2021-1-28 13:58):
筛法筛错了……少了两个该筛的数字，比如vec![105];应该是vec![15,105,195];

uk702

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

mathe

https://bbs.emath.ac.cn/thread-182-1-1.html

tprime

我的实现 https://github.com/ktprime/ktprime/blob/master/Ktprime.cpp
github 名称就是k（1-9）生素数,（素数，孪生，3生---9生素数或自定义类型， 最大范围到10^16）
10^12以内大概要几秒

运行如下：
AMD Ryzen 7 5800H with Radeon Graphics         L1/L2/L3 cache = 32/512/2057 kb
Prime quadruplets
: p, p + 2, p + 6, p + 8
Twin primes
thread(3) 2.09%, ~= 3.45 sec, Twin primes ~= 224150355
thread(4) 36.57%, ~= 3.16 sec, Twin primes ~= 224377560, err ~= 10.1363%%
thread(3) 71.04%, ~= 3.01 sec, Twin primes ~= 224343405, err ~= 822129631577663.8750%%
PI2(100000000000) = 224376048 (2.96 sec)

>> k41
Prime quadruplets
: p, p + 2, p + 6, p + 8

>> e12
thread(5) 1.26%, ~= 6.89 sec, Prime quadruplets ~= 8427510
thread(7) 22.10%, ~= 6.84 sec, Prime quadruplets ~= 8402940, err ~= 21888724040326892.0000%%
thread(5) 42.94%, ~= 6.86 sec, Prime quadruplets ~= 8417682, err ~= 17.5439%%
thread(7) 63.78%, ~= 6.38 sec, Prime quadruplets ~= 8390655, err ~= 21914280052049396.0000%%
thread(6) 84.62%, ~= 5.94 sec, Prime quadruplets ~= 8395569, err ~= 5.8565%%
PI4(1000000000000) = 8398278 (5.81 sec)

>> k61
Prime sextuplets
: p, p + 4, p + 6, p + 10, p + 12, p + 16

>> e13
thread(1) 4.03%, ~= 11.63 sec, Prime sextuplets ~= 298760
thread(2) 70.52%, ~= 11.81 sec, Prime sextuplets ~= 303380, err ~= 154.6392%%
PI6(10000000000000) = 303828 (12.19 sec)


>> h
        [P: Print time use]
        [D: Debug log]
        [S: Save result to file]
        [R: Runtime check pattern]
        [A: Save/Continue last task]
        [K: Calculate of Prime/Twin[Cousin]/Ktuplet Prime k(0 - 8)]
        [M: Monitor progress m(0 - 30)]
        [H: Help]
        [F: Factorial of wheel prime factor f(7 - 29)]
        [T: Threads number t(2 - 64)]
        [C: Cpu L1/L2 data cache size (L1:16-128, L2:128-1024)]

        [B: Benchmark (start) (end) (gptk)]
        [U: Unit test (n 1 - 10000) (gptk 0 - 2)]
        [N: Number of patterns (start)
        [I: List base pow index (powbase) (start) (end)]
        [L: List multi gptk (start) (end/count) (step)]


        All command/config as follow:
        B, B e9 e10 123
        C31, C128000
        K41, K52, KK26, KK268 T2-32
        H, Hk, A, D, S, R
        U, U 1000 012, U 1000+2 2
        N 120000*1000
        I 2 10 20
        L 2e9-100 1000 10
        L e9-100 2e9*2 1e8+2