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

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

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

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

×

首先发一个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 ./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![105];应该是vec![15,105,195];

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

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

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

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