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

.·.·. 发表于 2021-1-27 20:05:21

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

首先发一个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())过一会儿写筛法&&多线程
初步计数的结果如下：
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())
































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, &)] =
      &[(2_046, &),
          (1_373_652, &),
          (9_080_190, &),
          (25_326_000, &),
          (4_759_123_140, &),
          (1_112_004_669_632, &),
          (2_152_302_898_746, &),
          (3_474_749_660_382, &),
          (341_550_071_728_320, &),
          (0xFFFF_FFFF_FFFF_FFFF, &)
         ];

    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!("{:?}",);
            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!;
    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!;
    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 ./main
7
63
819
12285
233415
5835375
total:326240(from 0 to 77636318760)

real      0m19.608s
user      3m26.773s
sys      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 ./main
7
63
819
12285
233415
5835375
157555125
total:5907339(from 0 to 2406725881560)

real        12m9.250s
user        126m49.848s
sys        0m4.234s不知道大家还有什么神奇优化技巧

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

uk702 发表于 2021-1-28 11:13:29

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

mathe 发表于 2021-1-28 11:22:17

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

tprime 发表于 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
      
      
      
      
      
      /Ktuplet Prime k(0 - 8)]
      
      
      
      
      

      
      
      [N: Number of patterns (start)
      
      


      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

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