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

# [擂台] 一起竞速吗？——计算区间[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 的总个数，这样大家也方便核对检查。

 本帖最后由 tprime 于 2022-2-15 21:23 编辑 我的实现 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

 您需要登录后才可以回帖 登录 | 欢迎注册 本版积分规则 回帖并转播 回帖后跳转到最后一页

GMT+8, 2023-9-23 23:26 , Processed in 0.067982 second(s), 16 queries .

Powered by Discuz! X3.5

© 2001-2023 Discuz! Team.