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

mathematica

无心人 发表于 2020-11-26 16:18
用这种方式，测试的每10^7需要34.87s
10^14需要4035天，还需要优化~~~

你用的是什么软件？郭先强的软件？

无心人

100亿20并发测试开始，明天看结果
1000万4并发10秒

1479秒100亿，20并发，无结果~

.·.·.

离成功优化powermod大概只差一步了

1. fn pow_mod(mut base:u128,mut exp:u128,module:u128)->u128{
   
2.     // assert!(module < 1<<94 && module&1==1 && exp>0)
   
3.     let mut acc = 1;
   
4.     //let m64=(module>>64) as u64;
   
5.     let m126=(1<<126)%module;
   
6.     let m157=(m126<<31)%module;
   
7.     let mulmod=|mut x,mut y|{
   
8.         let [lx,ly,mut hx,hy]=[x&9223372036854775807,y&9223372036854775807,x>>63,y>>63];
   
9.         x=lx*hy+ly*hx;
   
10.         y=lx*ly+((x&9223372036854775807)<<63);x>>=63;//y=low_mid, x=high_mid<<126
    
11.         //lx*ly+y+x<<126+hx*hy<<126
    
12.         x+=hx*hy;
    
13.         //lx*ly+y+x<<126
    
14.         hx=x&2147483647;
    
15.         x>>=31;
    
16.         //(y+hx*m126+x*m157)%module
    
17.         y+=hx*m126+x*m157;
    
18.         //y-(((y>>64)as u64)/m64) as u128*module
    
19.         y%module
    
20.         //lx.max(ly) < 1<<63 ==> lx*ly < 1<<126
    
21.         // y < 9223372036854775807<<63 < 1<<126
    
22.         //hx <2147483647*module< 1<<31+94 = 1<<126
    
23.         // x <(module*module>>(63*2+31))*module < 1<<(94*2-63*2-31+94)=1<<125
    
24.         //Thus, lx*ly+y+hx*m126+x*m157 < 1<<128, safe to add them together and executing % ops.
    
25.     };
    
26.     while exp > 1 {
    
27.         if (exp & 1) == 1 {
    
28.             acc = mulmod(acc,base);
    
29.             base = mulmod(base,base);
    
30.         }else{
    
31.             base = mulmod(base,base);
    
32.         }
    
33.         exp >>= 1;
    
34.     }
    
35.     mulmod(acc,base)
    
36. }
    
37. fn main(){
    
38.     const CONSTARR:[u128;7]=[1u128<<1,1u128<<2,1u128<<4,1u128<<8,1u128<<16,1u128<<32,1u128<<64];
    
39.     let now=std::time::Instant::now();
    
40.     let arguments: Vec<String> = std::env::args().collect();
    
41.     let begin=arguments.get(2).and_then(|x|Some(x.parse::<u128>().unwrap_or(0))).unwrap_or(0);
    
42.     let len=arguments.get(1).and_then(|x|Some(x.parse::<u128>().unwrap_or(10000000))).unwrap_or(10000000);
    
43.     let mut am1=10000000000000000000000000000u128+begin*len;
    
44.     let mut a=10000000000000000000000000001u128+begin*len;
    
45.     println!("from={},len={}",a,len);
    
46.     let mut b=2;
    
47.     let mut am1d=0;
    
48.     let mut prime=0u64;
    
49.     let mut psp=0u64;
    
50.     for _i in 0..len>>1{
    
51.         let aa=am1.trailing_zeros().min(6);
    
52.         am1d=&am1>>aa;
    
53.         if {b=CONSTARR[aa as usize]; pow_mod(b,am1d,a)==1u128}{//Fermat
    
54.             psp+=1;/*
    
55.             if {b.assign(2u8);//Miller-Rabin
    
56.                 let aa=a.find_one(0).unwrap();
    
57.                 am1d.assign(&am1>>aa);
    
58.                 let _=b.pow_mod_mut(&am1d,&a);
    
59.                 (b==1 || b==am1 )||{
    
60.                     let mut f=false;
    
61.                     for _ in 0..aa{
    
62.                         let _=b.pow_mod_mut(&Integer::from(2),&a);
    
63.                         if b==am1{f=true;break}
    
64.                     }
    
65.                     f
    
66.                 }
    
67.             }{
    
68.                 if a.is_probably_prime(25)==rug::integer::IsPrime::No{
    
69.                     println!("{} is SPSP(2)",a)
    
70.                 }else{prime+=1}
    
71.             }else{println!("{} is PSP(2)",a)}*/
    
72.         }
    
73.         a+=2;
    
74.         am1+=2;
    
75.     }
    
76.     println!("{}.{},{},{:?}",a,prime,psp-prime,now.elapsed());
    
77. }

复制代码

这个程序是gmp的1/2，但如果拿掉取模步，速度是gmp的两倍多。
另外，或许应该直接写C调GMP，GMP写了Miller-Rabin测试，不知速度如何，但小改一下应该比Fermat快。
毕竟Miller-Rabin可以少算好几次乘法。

.·.·.

无心人 发表于 2020-11-30 17:04
100亿20并发测试开始，明天看结果
1000万4并发10秒

？？
把Primes.jl改一改吧
改完，1000万不并发也就11秒

刚刚试了一下
如果一定不想改Primes.jl的话
这样启动Julia（linux，windows最好先用set改一下环境变量）
JULIA_NUM_THREADS=6 julia

1. #不使用Primes，直接手写一个isprime
   
2. isprime(x::BigInt) = ccall((:__gmpz_probab_prime_p,:libgmp),Cint, (Ref{BigInt}, Cint), x, 25) > 0
   
3. const b = big(10)^28
   
4. 
   
5. function fermatProbablePrimeTest2(n)
   
6.     a = powermod(2, n-1, n)
   
7.     return a == 1
   
8. end
   
9. 
   
10. function millerRabinProbablePrimeTest( n, b )
    
11.     k = trailing_zeros(n-1)
    
12.     m = (n-1) >> k
    
13.     a = powermod(b, m, n)
    
14.     if (a == 1) || (a == (n-1))
    
15.         return true
    
16.     else
    
17.         for i in 1:k-1
    
18.             a = powermod(a, 2, n)
    
19.             if a == n - 1
    
20.                 return true
    
21.             end
    
22.         end
    
23.     end
    
24.     return false
    
25. end
    
26. @time begin
    
27.     Threads.@threads for i in range(1, step = 2, stop = 9999999)
    
28.         test = b + i
    
29.         if fermatProbablePrimeTest2(test)
    
30.             if millerRabinProbablePrimeTest(test, 2)
    
31.                 if !isprime(test)
    
32.                     println("$test is SPSP(2)")
    
33.                 end
    
34.             else
    
35.                 println("$test is PSP(2)")
    
36.             end
    
37.         end
    
38.     end
    
39. end
    

复制代码

无心人

我刚思考了下，用frobenius扫同样的10^10的时间太长了，不太合适（现在还没出结果~~~，看来frobenius测试实在不适合批量测试~~~）
还是用2，3，5，7，11，13六个MR测试代替吧
这种理论上，有1/4^6的会漏掉，大概是0.0244%，实际上很小
下面可能测试下，六次MR测试失败的概率，是不是会低于亿分之一~~~~

因为我们是统计，所以，对于很小概率的错误，是可以允许的

.·.·.

无心人 发表于 2020-11-30 21:59
我刚思考了下，用frobenius扫同样的10^10的时间太长了，不太合适（现在还没出结果~~~，看来frobenius测试实 ...

……你试试去掉isprime扫前10^7个数字
如果用isprime(x::BigInt) = ccall((:__gmpz_probab_prime_p,:libgmp),Cint, (Ref{BigInt}, Cint), x, 25) > 0
开销也就额外增加1/10，而不是10倍
整个程序我已经给了
难道有BUG不能执行吗？

我昨天扫了600亿，用了大约10300秒，一个Fermat伪素数都没看见（用Rust）
换Julia，6线程大概13000秒？
isprime(x::BigInt) = ccall((:__gmpz_probab_prime_p,:libgmp),Cint, (Ref{BigInt}, Cint), x, 25) > 0
这一行可以带来非常可怕的性能提升

无心人

我说下我的筛法思路
初始化一个表small，表里是10000内大于3的素数和模2的阶

建立数组b[]，长度是s
for i=0 to (s-1)/2
        b[2*i]=0
        b[2*i+1]=1

1、设定初始值为base，则s[i]位置代表的整数是base+i
for (p, d) in small
        if d是奇数
                d = 2*d
        j = base/p
        start=j*p
        if start是偶数
                start += p
                j=j+1
        stopJ=(base+s-1)/p
        while j<=stopJ
                i = start - base
                if (j%d!=1)
                  s[i] = 0
                start=start+2*p
                j=j+2
               

列表里所有的p处理完后
2、遍历s[i]
  for i in 0..s-1
        if s[i] == 1
          s[i] = FermatTest(base+i, 2)
       
3、位标1的，对2做MR测试，没通过的置0，输出PSP
4、位标1的，对3，5，7，11，13做MR测试，没通过的输出SPSP

无心人

有没有人写段汇编的128位模乘啊

.·.·.

无心人 发表于 2020-12-1 14:43
有没有人写段汇编的128位模乘啊

我不知道该怎么写128bit的模
但如果不做这次取模，fermat的速度可以乘以2
这大概表示，fermat的优化已经到头了，就算汇编未必比gmp快多少

无心人

.·.·. 发表于 2020-12-1 18:54
我不知道该怎么写128bit的模
但如果不做这次取模，fermat的速度可以乘以2
这大概表示，fermat的优化已 ...

我感觉模乘不需要除法