找回密码
 欢迎注册
查看: 8970|回复: 3

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

[复制链接]
发表于 2021-1-27 20:05:21 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有账号?欢迎注册

×
首先发一个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.     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
复制代码
不知道大家还有什么神奇优化技巧

补充内容 (2021-1-28 13:58):
筛法筛错了……少了两个该筛的数字,比如vec![105];应该是vec![15,105,195];
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2021-1-28 11:13:29 | 显示全部楼层
你的  5835375 和 157555125 分别对应的 m、n 是多少?最好能提供一个标准的输出,比如,输出 10^8 内所有满足 k+1,k+3,k+7,k+9 都为素数的 k,以及这样的 k 的总个数,这样大家也方便核对检查。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2021-1-28 11:22:17 | 显示全部楼层
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2022-2-15 21:17:26 | 显示全部楼层
本帖最后由 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
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-3-29 17:17 , Processed in 0.285754 second(s), 16 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表