10^28开始的10^14个整数的素性概率性测试算法实践报告

无心人 · 发表于 2020-12-9 14:06:26

决定了，程序内不做并发了~，同步太麻烦了~
2.5G主频，非并发，46秒~~~

开始时间: 15:24:58
初始化结束: 15:24:59
筛结束: 15:25:02
测试结束: 15:25:45
通过筛的数字: 12090362
可能素数: 1552676

.·.·. · 发表于 2020-12-9 15:14:21

无心人发表于 2020-12-9 14:06
决定了，程序内不做并发了~，同步太麻烦了~
2.5G主频，非并发，46秒~~~

说一下我的程序问题
我把ord除以2的原因是，sieve事实上是乘过2的
sieve增加1的意思是，我们筛的数字大小加2
这样可以节省一部分内存占用。

另外，才发现，BPSW真的好慢好慢
//一次MillerRabin不用BPSW
Finished release [optimized] target(s) in 6.27s
Running `target/release/fermat`
from=10000000000000000000000000001,len=500
10000000000000000000099999979.0,1552676,12.651374079s
//用一次BPSW
Compiling fermat v0.1.0 (/me/fermat)
Finished release [optimized] target(s) in 6.84s
Running `target/release/fermat`
from=10000000000000000000000000001,len=500
10000000000000000000099999979.1552676,0,29.689700734s

无心人 · 发表于 2020-12-9 16:11:26

筛掉素数平方后，待测试素数从1209万降低到1046万，我是10000内素数的筛，3不参与筛，但是筛9的倍数
2.5G处理器，初始化和筛，1秒左右，测试39秒

平方数测试涉及到以下几个
1、3不筛3*(2k+1)，但是筛3的平方9的倍数
2、1093，3511的平方不筛
3、假设筛是s个连续整数，超过sqrt(s/2)的素数，在筛的范围内平方倍数且是奇数的数字可能并不都能筛掉未筛除数字的，是不是做，需要考虑

刚测试，因为我筛是10^8，sqrt(50000000)=7071.1, 实际只筛小于等于7071的平方，和全筛，筛出来的需要做模幂测试的个数是一样的

无心人 · 发表于 2020-12-9 19:48:28

真郁闷，暴力算1000亿，结果回来远程进去发现，服务器重启了~~~
现在缺少个任务分配平台，暴力塞了1000任务过去，压力太大，崩了服务器
不知道有什么平台，能批量均匀的分配任务，并且自己回收输出

.·.·. · 发表于 2020-12-9 22:41:28

本帖最后由 .·.·. 于 2020-12-10 04:43 编辑

无心人发表于 2020-12-9 19:48
真郁闷，暴力算1000亿，结果回来远程进去发现，服务器重启了~~~
现在缺少个任务分配平台，暴力塞了1000任 ...

直接python？
随便找个threadpool之类的东西自己写一写subprocess.popen就好

话说这玩意难道不应该均分成20个任务提交给20个不同的程序吗？
---
update：fermat正式退出历史舞台

#![feature(slice_fill)]
use rug::{Assign, Integer};
const SIEVE_SIZE_DIV_2:i32=100000;
//const SIEVE_SIZE:i32=SIEVE_SIZE_DIV_2<<1;//the true sieve size, never used.
const PRIME_SIZE:usize=include!("final.txt").len();
const PRIME_MAP:[[i32;2];PRIME_SIZE]=include!("final.txt");
fn main() {
let mut cumsum=0u32;
let now=std::time::Instant::now();
let arguments: Vec<String> = std::env::args().collect();
let begin=if let Some(x)=arguments.get(2){x.parse::<u128>().unwrap_or(0)}else{0};
let len=if let Some(x)=arguments.get(1){x.parse::<u128>().unwrap_or(500)}else{500};
let start=10000000000000000000000000000u128+begin*len;
let mut am1=Integer::from(start);
let mut a=Integer::from(start+1);
println!("from={},len={}",a,len);
let mut b=Integer::from(2);
let mut am1d=Integer::from(0);
let mut prime=0u64;
let mut psp=0u64;
let mut sieve=[true;SIEVE_SIZE_DIV_2 as usize];
let mut prime_curr=[0;PRIME_SIZE];
let mut prime_order=[0;PRIME_SIZE];
let [mut total,mut tests]=[0;2];
for i in 0..PRIME_SIZE{
let t=(start%PRIME_MAP[i][0] as u128) as i32;
let next= if t==0{
start+PRIME_MAP[i][0] as u128
}else{
let mut tmp=start-(t as u128)+PRIME_MAP[i][0] as u128;
if tmp%2==0{tmp+=PRIME_MAP[i][0] as u128}
tmp
};
let count=next/PRIME_MAP[i][0] as u128;
prime_curr[i]=((next-start)>>1) as i32;
if PRIME_MAP[i][1]==-1{// for p^2, should never put 1093^2 and 3511^2 into prime_map.
prime_order[i]=-1;// if o is order of 2 mod p , op is the order of 2 mod p^2(except p=1093,3511...), thus, if 2^(qp^2-1)=1 mod p^2, op|qp^2-1 , p|-1, which is not possible.
}else{
prime_order[i]=PRIME_MAP[i][1]-((count>>1)%PRIME_MAP[i][1] as u128) as i32
}
}
for _i in 0..len{
sieve.fill(true);
for (i,&[prime,order]) in PRIME_MAP.iter().enumerate(){
let mut cur=prime_curr[i];
let mut ord=prime_order[i];
while cur<SIEVE_SIZE_DIV_2{
if ord==0{
ord=order
}else{
sieve[cur as usize]=false
}
cur+=prime;
ord-=1
}
prime_curr[i]=cur-SIEVE_SIZE_DIV_2;
prime_order[i]=ord
}
for test in &sieve{
total+=1;
if *test{
tests+=1;
a+=cumsum;
am1+=cumsum;
cumsum=0;
let aa=am1.find_one(0).unwrap().min(6);
am1d.assign(&am1>>aa);
// if {b.assign(CONSTARR[aa as usize]);let _=b.pow_mod_mut(&am1d,&a);b==1u8}{psp+=1;//Fermat
if {b.assign(2u8);//Miller-Rabin
let aa=a.find_one(0).unwrap();
am1d.assign(&am1>>aa);
let _=b.pow_mod_mut(&am1d,&a);
(b==1 || b==am1 )||{
let mut f=false;
for _ in 0..aa{
let _=b.pow_mod_mut(&Integer::from(2),&a);
if b==am1{f=true;break}
}
if b==1{f=false;println!("{} is PSP(2)",a);psp+=1}
f
}
}{
psp+=1;
prime+=1;
//if a.is_probably_prime(5)==rug::integer::IsPrime::No{
// println!("{} is SPSP(2)",a)
//}else{prime+=1}
}
// }//Fermat
}
cumsum+=2;
}
}
println!("last={},prime={},diff={},total={},miller={},cost={:?}",a,prime,psp-prime,total,tests,now.elapsed());
}

复制代码

fermat.7z (92.53 KB, 下载次数: 0)
（去掉is_probable_prime之后得到）
from=10000000000000000000000000001,len=500
last=10000000000000000000099999979,prime=1552676,diff=0,total=50000000,miller=8377755,cost=10.800631943s
共8377755个数字需要继续送miller-rabin进行计算

无心人 · 发表于 2020-12-10 09:08:47

时间降不下来，发现筛出后，标记合数的只做fermat也不行，降不下来时间，不过候选素数从1000+万降低到600+万

Start Time: 09:11:55
init Time: 09:11:55
start number: 10000000000000000000000000000
sieve Time: 09:11:56
test Time: 09:12:36
to test PSP: 4375043, to test primility 6087984
primes: 1552676

4375043里，我怀疑3的倍数占很大比例

#include<gmp.h>
#include<stdbool.h>
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#define NULL_TERMINATED 0
bool fermatTest( const mpz_t n, unsigned long b)
{
mpz_t r, b1, ns1;
mpz_inits(r, ns1, NULL_TERMINATED);
mpz_init_set_ui(b1, b);
mpz_sub_ui(ns1, n, 1);
mpz_powm(r, b1, ns1, n);
bool result = mpz_cmp_ui(r, 1) == 0;
mpz_clears(r, b1, ns1, NULL_TERMINATED);
return result;
}
bool millerrabinTest( const mpz_t n, unsigned long b )
{
mpz_t ns1, t, b1, r;
unsigned long s;
bool result;
mpz_inits(ns1, t, r, NULL_TERMINATED);
mpz_init_set_ui(b1, b);// b1=b
mpz_sub_ui(ns1, n, 1);// ns1=n-1
s = mpz_scan1(ns1, 0);
mpz_tdiv_q_2exp(t, ns1, s);// ns1=2^s*t, t%2=1
mpz_powm(r, b1, t, n);
if (mpz_cmp_ui(r, 1) == 0 || mpz_cmp(r, ns1) == 0)
result = true;
else
{
result = false;
for (unsigned long i = 1; i < s; i ++)
{
mpz_powm_ui(r, r, 2, n);
if (mpz_cmp(r, ns1) == 0)
{
result = true;
break;
}
}
}
mpz_clears(ns1, t, b1, r, NULL_TERMINATED);
return result;
}
const unsigned long step = 100000000;
char *buff;
unsigned long cnt = 0, pcnt = 0, pscnt = 0;
unsigned long small[1227][2];
unsigned long p100[] = { 2, 3, 5, 7, 11,
13, 17, 19, 23, 29,
31, 37, 41, 43, 47,
53, 59, 61, 67, 71,
73, 79, 83, 89, 97};
//p < 10000
unsigned long ord(unsigned long p)
{
unsigned long factor[16];
unsigned long s = 0;
unsigned long n = p - 1;
mpz_t b, tmp, mp;
mpz_init_set_ui(b, 2);
mpz_init_set_ui(mp, p);
mpz_init(tmp);
for (int i = 0; i < 25; i ++)
{
if (n%p100[i] == 0)
do
{
factor[s++] = p100[i];
n /= p100[i];
} while (n%p100[i] == 0);
}
n = p - 1;
for (int i = 0; i < s; i++)
{
mpz_powm_ui(tmp, b, n / factor[i], mp);
if (mpz_cmp_ui(tmp, 1) == 0)
n /= factor[i];
}
if (n%2==1)
n*=2;
mpz_clears(b, tmp, mp, NULL_TERMINATED);
return n;
}
// 实际测试值=base + index*step + [1..step]
void init( mpz_t base, unsigned long index )
{
mpz_t tmp, tmp1;
mpz_set_ui(base, 10);
mpz_pow_ui(base, base, 28);
mpz_init_set_ui(tmp, step);
mpz_mul_ui(tmp, tmp, index);
mpz_add(base, base, tmp);
buff = (char *)malloc(step);
for (long i = 0; i < step; i += 2)
buff[i] = 0;
for (long i = 1; i < step; i += 2)
buff[i] = 3; //3 = 位组合0:素数 1：PSP(2)
unsigned long i = 0;
mpz_init_set_ui(tmp1, 3);
while (i < 1227)
{
mpz_nextprime(tmp, tmp1);
mpz_set(tmp1, tmp);
small[i][0] = mpz_get_ui(tmp);
small[i][1] = ord(small[i][0]);
i++;
}
mpz_clears(tmp, tmp1, NULL_TERMINATED);
}
bool quickTest( mpz_t n )
{
unsigned long b[] = {3, 5, 7, 11, 13};
for (int i = 0; i < sizeof(b)/sizeof(unsigned long); i++)
if (!millerrabinTest(n, b[i]))
return false;
return true;
}
void sieve( mpz_t base, long idx )
{
mpz_t q, r, q1;
unsigned long r1, q2;
mpz_inits(q, r, q1, NULL_TERMINATED);
//置位奇数中3的倍数为非素数
mpz_mod_ui(q1, base, 3);
q2 = mpz_get_ui(q1);
q2 = 3 - q2;
if (q2%2==0)
q2 += 3;
while (q2 < step)
{
buff[q2] &= 2;
q2 += 6;
}
for (long i = 0; i < 1227; i ++)
{
//合数形如pm, 素数p对2的次数是k, 筛掉不满足m=nk+1的，因为只需要考虑m是奇数，所以，对k不是偶数的，已经乘以2了
mpz_fdiv_qr_ui(q, r, base, small[i][0]);
mpz_mod_ui(q1, q, small[i][1]);
q2 = mpz_get_ui(q1);
r1 = mpz_get_ui(r);
r1 = small[i][0] - r1;
q2 = (q2+1)%small[i][1];
if (r1%2==0)
{
r1 += small[i][0];
q2 = (q2+1)%small[i][1];
}
while (r1 < step)
{
buff[r1] &= 2; //非素数
if (q2 != 1)
{
buff[r1] = 0; //非PSP（2）
}
r1 += 2*small[i][0];
q2 = (q2+2)%small[i][1];
}
//筛掉满足p^2的倍数，其中p是素数且不等于1093，3511
unsigned long sq = small[i][0];
if (sq != 1093 && sq != 3511 && sq <= 7071)
{
sq *= sq;
mpz_mod_ui(q1, base, sq);
q2 = mpz_get_ui(q1);//取模
q2 = sq - q2;//加sq，得到第一个奇数倍数，如果q2=0,结果也是对的
if (q2%2==0)
q2 += sq;
while (q2 < step)
{
buff[q2] = 0; //非PSP（2）
q2 += 2*sq;
}
}
}
//筛掉奇数中9的倍数
mpz_mod_ui(q1, base, 9);
q2 = mpz_get_ui(q1);
q2 = 9 - q2;
if (q2%2==0)
q2 += 9;
while (q2 < step)
{
buff[q2] = 0;
q2 += 18;
}
mpz_clears(q, r, q1, NULL_TERMINATED);
}
//位标0：素数 1：PSP(2)
int check(mpz_t base, long i)
{
mpz_t n;
int result;
mpz_init(n);
mpz_add_ui(n, base, i);
if (fermatTest(n, 2))
{
if (millerrabinTest(n, 2))
{
if (quickTest(n))
{
result = 3;
cnt ++;
}
else
result = 2;
}
else
result = 1;
}
else
result = 0;
if (result > 0 && result < 3)
{
gmp_printf("%Zd type: %d\n", n, result);
}
mpz_clear(n);
return result;
}
void test( mpz_t base, long idx )
{
for (long i = 1; i < step; i += 2)
if (buff[i] == 3) //未筛的数字
{
pcnt ++;
buff[i] = check(base, i);
}
else
if (buff[i] == 2) //判定非素数，且可能PSP（2）的数字
{
mpz_t n;
pscnt++;
mpz_init_set_ui(n, i);
mpz_add(n, n, base);
if (fermatTest( n, 2))
gmp_printf("%Zd: PSP(2)\n", n);
mpz_clear(n);
}
}
void clear( void )
{
free(buff);
}
void showTime(char * pref)
{
time_t rawtime;
struct tm * timeinfo;
char buffer [16];
time (&rawtime);
timeinfo = localtime (&rawtime);
strftime (buffer, 16 ,"%X", timeinfo);
printf("%s: %s\n", pref, buffer);
}
int main( int argc, char * argv[] )
{
mpz_t base;
mpz_init(base);
long idx = 0;
char * ptr;
if (argc > 1)
idx = strtol(argv[1], &ptr, 10);
showTime("Start Time");
init( base, idx );
showTime(" init Time");
gmp_printf("start number: %Zd\n", base);
sieve( base, idx );
showTime("sieve Time");
test( base, idx );
showTime(" test Time");
printf("to test PSP: %d, to test primility %d\n", pscnt, pcnt);
printf("primes: %d\n", cnt);
mpz_clear(base);
clear();
return 0;
}

复制代码

.·.·. · 发表于 2020-12-10 18:02:17

无心人发表于 2020-12-10 09:08
时间降不下来，发现筛出后，标记合数的只做fermat也不行，降不下来时间，不过候选素数从1000+万降低到600+ ...

话说你不会把全体小于5*10^7的素数都扔进去筛素数了吧……
可怕……

无心人 · 发表于 2020-12-22 17:23:15

筛到100万内素数，总时间从40秒降低到35秒

无心人 · 发表于 2020-12-23 09:46:59

#include<gmp.h>
#include<stdbool.h>
#include<stdint.h>
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#define NULL_TERMINATED 0
#define primesize 78494
const uint64_t step = 100000000;
uint64_t cnt = 0, pcnt = 0, pscnt = 0;
uint64_t small[primesize][2];
char * buff;
// 实际测试值=base + index*step + [1..step]
int init( mpz_t base, uint64_t index )
{
char line[64];
FILE *ordFile = fopen("ord6raw.txt", "r");
if (ordFile == NULL)
{
printf("file ord6raw not find!\n");
return -1;
}
uint64_t i = 0;
while (fgets(line, 64, ordFile))
{
uint64_t p, o;
sscanf(line, "%lu%lu", &p, &o);
small[i][0] = p;
small[i++][1] = o;
if (i >= primesize)
break;
}
fclose(ordFile);
mpz_t tmp;
mpz_set_ui(base, 10);
mpz_pow_ui(base, base, 28);
mpz_init_set_ui(tmp, step);
mpz_mul_ui(tmp, tmp, index);
mpz_add(base, base, tmp);
mpz_clear(tmp);
buff = (char *)malloc(step);
for (long i = 0; i < step; i += 2)
buff[i] = 0;
for (long i = 1; i < step; i += 2)
buff[i] = 3; //3 = 位组合0:素数 1：PSP(2)
}
bool fermatTest( const mpz_t n, uint64_t b)
{
mpz_t r, b1, ns1;
mpz_inits(r, ns1, NULL_TERMINATED);
mpz_init_set_ui(b1, b);
mpz_sub_ui(ns1, n, 1);
mpz_powm(r, b1, ns1, n);
bool result = mpz_cmp_ui(r, 1) == 0;
mpz_clears(r, b1, ns1, NULL_TERMINATED);
return result;
}
bool millerrabinTest( const mpz_t n, uint64_t b )
{
mpz_t ns1, t, b1, r;
uint64_t s;
bool result;
mpz_inits(ns1, t, r, NULL_TERMINATED);
mpz_init_set_ui(b1, b);// b1=b
mpz_sub_ui(ns1, n, 1);// ns1=n-1
s = mpz_scan1(ns1, 0);
mpz_tdiv_q_2exp(t, ns1, s);// ns1=2^s*t, t%2=1
mpz_powm(r, b1, t, n);
if (mpz_cmp_ui(r, 1) == 0 || mpz_cmp(r, ns1) == 0)
result = true;
else
{
result = false;
for (uint64_t i = 1; i < s; i ++)
{
mpz_powm_ui(r, r, 2, n);
if (mpz_cmp(r, ns1) == 0)
{
result = true;
break;
}
}
}
mpz_clears(ns1, t, b1, r, NULL_TERMINATED);
return result;
}
bool quickTest( mpz_t n )
{
uint64_t b[] = {3, 5, 7, 11, 13};
for (int i = 0; i < sizeof(b)/sizeof(uint64_t); i++)
if (!millerrabinTest(n, b[i]))
return false;
return true;
}
void sievep( mpz_t base, uint64_t p, uint64_t o)
{
mpz_t q, r, q1;
uint64_t r1, q2;
mpz_inits(q, r, q1, NULL_TERMINATED);
//合数形如pm, 素数p对2的次数是o, 筛掉不满足m=no+1的，因为只需要考虑m是奇数，所以，对o不是偶数的，已经乘以2了
mpz_fdiv_qr_ui(q, r, base, p);
mpz_mod_ui(q1, q, o);
q2 = mpz_get_ui(q1);
r1 = mpz_get_ui(r);
r1 = p - r1;
q2 = (q2+1)%o;
if (r1%2==0)
{
r1 += p;
q2 = (q2+1)%o;
}
while (r1 < step)
{
buff[r1] &= 2; //非素数
if (q2 != 1)
{
buff[r1] = 0; //非PSP（2）
}
r1 += 2*p;
q2 = (q2+2)%o;
}
mpz_clears(q, r, q1, NULL_TERMINATED);
}
void sievep2( mpz_t base, uint64_t sp )
{
mpz_t q1;
uint64_t q2;
mpz_init(q1);
mpz_mod_ui(q1, base, sp);
q2 = mpz_get_ui(q1);//取模
q2 = sp - q2;//加sq，得到第一个奇数倍数，如果q2=0,结果也是对的
if (q2%2==0)
q2 += sp;
while (q2 < step)
{
buff[q2] = 0; //非PSP（2）
q2 += 2*sp;
}
}
void sieve( mpz_t base )
{
mpz_t q1;
uint64_t q2;
mpz_init(q1);
//置位奇数中3的倍数为非素数
mpz_mod_ui(q1, base, 3);
q2 = mpz_get_ui(q1);
q2 = 3 - q2;
if (q2%2==0)
q2 += 3;
while (q2 < step)
{
buff[q2] &= 2;
q2 += 6;
}
for (long i = 0; i < primesize; i ++)
sievep( base, small[i][0], small[i][1] );
sievep( base, 1093, 364 );
sievep( base, 3511, 3510 );
//筛掉奇数中9的倍数
sievep2( base, 9 );
//筛掉满足p^2的倍数，其中p是小于sqrt(step)素数且不等于1093，3511
//筛选范围内有大于sqrt(step)的素数的平方的因子的数字很少
for (long i = 0; i < 1227; i ++)
sievep2( base, small[i][0] * small[i][0]);
mpz_clear(q1);
}
//位标0：素数 1：PSP(2)
int check(mpz_t base, long i)
{
mpz_t n;
int result;
mpz_init(n);
mpz_add_ui(n, base, i);
if (fermatTest(n, 2))
{
if (millerrabinTest(n, 2))
{
if (quickTest(n))
{
result = 3;
cnt ++;
}
else
result = 2;
}
else
result = 1;
}
else
result = 0;
if (result > 0 && result < 3)
{
gmp_printf("%Zd type: %d\n", n, result);
}
mpz_clear(n);
return result;
}
void test( mpz_t base, long idx )
{
for (long i = 1; i < step; i += 2)
if (buff[i] == 3) //未筛的数字
{
pcnt ++;
buff[i] = check(base, i);
}
else
if (buff[i] == 2) //判定非素数，且可能PSP（2）的数字
{
mpz_t n;
pscnt++;
mpz_init_set_ui(n, i);
mpz_add(n, n, base);
if (fermatTest( n, 2))
gmp_printf("%Zd: PSP(2)\n", n);
mpz_clear(n);
}
}
void clear( void )
{
free(buff);
}
void showTime(char * pref)
{
time_t rawtime;
struct tm * timeinfo;
char buffer [16];
time (&rawtime);
timeinfo = localtime (&rawtime);
strftime (buffer, 16 ,"%X", timeinfo);
printf("%s: %s\n", pref, buffer);
}
int main( int argc, char * argv[] )
{
mpz_t base;
mpz_init(base);
uint64_t idx = 0;
char * ptr;
if (argc > 1)
idx = strtol(argv[1], &ptr, 10);
showTime("Start Time");
if (init( base, idx ) == -1)
exit(-1);
showTime(" init Time");
gmp_printf("start number: %Zd\n", base);
sieve( base );
showTime("sieve Time");
test( base, idx );
showTime(" test Time");
printf("to test PSP: %lu, to test primility %lu\n", pscnt, pcnt);
printf("primes: %lu\n", cnt);
mpz_clear(base);
clear();
return 1;
}

复制代码

Start Time: 01:44:11
init Time: 01:44:11
start number: 10000000000000000000000000000
sieve Time: 01:44:14
test Time: 01:44:46
to test PSP: 2920776, to test primility 4062724
primes: 1552676

用素数模2的阶筛100万内素数，然后筛掉1万内素数（不包含1093，3511）的平方

无心人 · 发表于 2020-12-23 09:57:47

小素数范围	需要测试的疑似PSP	需要测试的其他数字	时间
10万	3503020	4874685	37
100万	2920776	4062724	35
10000万	2187737	3044791	34

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

[原创] 10^28开始的10^14个整数的素性概率性测试算法实践报告

点评

点评

点评

点评

点评

点评

点评