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

litong · 发表于 2023-9-25 12:12:16

首先发一个Pari/GP的代码

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

复制代码

如果只准备计数的话，可以用这个

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]如下：

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)]))[1, 0][2, 0][3, 0][4, 1][5, 1][6, 1][7, 2][8, 2][9, 2][10, 2][11, 2][12, 4][13, 4][14, 6][15, 8][16, 11][17, 13][18, 18][19, 31][20, 53][21, 109][22, 160][23, 266][24, 424][25, 709][26, 1227][27, 2039][28, 3421][29, 5826][30, 10130][31, 17511][32, 30633]cpu time = 26,537 ms, real time = 26,572 ms.

复制代码

update: (naive)Rust
首先放一个头文件,miller-rabin.rs，用witness进行素性判定

fn mod_exp(mut x: u64, mut d: u64, n: u64) -> u64 { let mut ret: u64 = 1; if n<4294967296{ while d != 0 { if d & 1 == 1 { ret = (ret*x)%n } d >>= 1; x = (x*x)%n; } }else{ while d != 0 { if d & 1 == 1 { ret = ((ret as u128*x as u128)%n as u128) as u64 } d >>= 1; x = ((x as u128*x as u128)%n as u128) as u64; } } ret}pub fn miller_rabin(n: u64) -> bool { // we have a strict upper bound, so we can just use the witness // table of Pomerance, Selfridge & Wagstaff and Jeaschke to be as // efficient as possible, without having to fall back to // randomness. const WITNESSES: &[(u64, &[u64])] = &[(2_046, &[2]), (1_373_652, &[2, 3]), (9_080_190, &[31, 73]), (25_326_000, &[2, 3, 5]), (4_759_123_140, &[2, 7, 61]), (1_112_004_669_632, &[2, 13, 23, 1662803]), (2_152_302_898_746, &[2, 3, 5, 7, 11]), (3_474_749_660_382, &[2, 3, 5, 7, 11, 13]), (341_550_071_728_320, &[2, 3, 5, 7, 11, 13, 17]), (0xFFFF_FFFF_FFFF_FFFF, &[2, 3, 5, 7, 11, 13, 17, 19, 23]) ]; let mut d = n - 1; let mut s = 0; while d & 1 == 0 { d /= 2; s += 1 } let witnesses = WITNESSES.iter().find(|&&(hi, _)| hi >= n) .map(|&(_, wtnss)| wtnss).unwrap(); 'next_witness: for &a in witnesses.iter() { let mut power = mod_exp(a, d, n); assert!(power < n); if power == 1 || power == n - 1 { continue 'next_witness } if n<4294967296{ for _r in 0..s { power = (power*power)%n; if power == 1 { return false } if power == n - 1 { continue 'next_witness } } }else{ for _r in 0..s { power = ((power as u128*power as u128)%n as u128) as u64; if power == 1 { return false } if power == n - 1 { continue 'next_witness } } } return false } true}

复制代码

程序：

include!{"miller-rabin.rs"}fn main() { let mut i=0u64.wrapping_sub(1); let mut n=1; let mut ctr=1; let mut limit=1<<n; loop{ i+=2; let j=i*105; if j>limit { println!("{:?}",[n,ctr-(n<=3) as i32]); limit<<=1; n+=1; if n>34{break}else{continue} } if miller_rabin(j-4) && miller_rabin(j-2) && miller_rabin(j+2) && miller_rabin(j+4) { ctr+=1 } }}

复制代码

update2: naive Rust的变种（向Rayon以及筛法靠拢）

include!{"miller-rabin.rs"}//use rayon::prelude::*;fn main() { let mut prepare=vec![105]; let mut prepare_len=210; println!("total:{}",1+(0..(1u64<<32)/prepare_len).into_iter().fold(0,|s,i|s+{ let start=prepare_len*i; prepare.iter().map(|&x|x+start).fold(0,|s,j|{ if miller_rabin(j-4) && miller_rabin(j-2) && miller_rabin(j+2) && miller_rabin(j+4) { s+1 }else{ s } }) }))}

复制代码

Rayon+筛法，为了保证计算效率，程序计数最好停止在筛的边界（而非之前随便写的那个边界）

include!{"miller-rabin.rs"}extern crate rayon;use rayon::prelude::*;fn not_divide(d:u64,n:u64)->bool{ (n/d)*d!=n}fn sieve(v:&mut Vec<u64>,len:&mut u64,prime:u64){ 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)&¬_divide(prime,x-2)&¬_divide(prime,x+2)&¬_divide(prime,x+4)).collect());s}) .reduce(||Vec::new(),|mut s,mut x|{s.append(&mut x);s}); *len*=prime; vv.shrink_to_fit(); vv.sort_unstable(); println!("{:?}",vv.len()); *v=vv}fn main() { let end = 12; let mut prepare=vec![105]; let mut prepare_len=210; sieve(&mut prepare,&mut prepare_len,11); sieve(&mut prepare,&mut prepare_len,13); sieve(&mut prepare,&mut prepare_len,17); sieve(&mut prepare,&mut prepare_len,19); sieve(&mut prepare,&mut prepare_len,23); sieve(&mut prepare,&mut prepare_len,29); sieve(&mut prepare,&mut prepare_len,31); println!("total:{:?}(from 0 to {})",1+(0u64..end).into_par_iter().fold(||0,|s,i|s+{ let start=prepare_len*i; prepare.iter().map(|&x|x+start).fold(0,|s,j|{ if miller_rabin(j-4) && miller_rabin(j-2) && miller_rabin(j+2) && miller_rabin(j+4) { s+1 }else{ s } }) }).reduce(||0,|a,b|a+b),end*prepare_len)}

复制代码

算得

$ 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 ./main763819122852334155835375total:326240(from 0 to 77636318760)real 0m19.608suser 3m26.773ssys 0m0.110s

复制代码

以及

$ 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 ./main763819122852334155835375157555125total:5907339(from 0 to 2406725881560)real 12m9.250suser 126m49.848ssys 0m4.234s

复制代码

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

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

ogemama · 发表于 2023-9-25 12:12:33

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

oabipalehov · 发表于 2023-9-25 12:12:48

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

erawaliktevi · 发表于 2023-9-25 12:13:14

我的实现 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

		自动登录	找回密码
密码			立即注册