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

无心人

.·.·. 发表于 2020-11-28 12:21
fermat可以省略
    if gcd(n, b) > 1
        return false

用fermat的原因是想测试下PSP2基2伪素数的概率，另外，打算用frobenius测试做最终判定是不是素数

.·.·.

无心人 发表于 2020-11-28 14:32
用fermat的原因是想测试下PSP2基2伪素数的概率，另外，打算用frobenius测试做最终判定是不是素数

我说的是
fermat那个算法开省略这一个判断：
    if gcd(n, b) > 1
        return false
    end

然而程序排版导致断句出了问题……

无心人

1. 
   
2. #找到最小的d满足b^d = 1 (mod n)
   
3. function ord(n, b)
   
4.     #if !isprime(n)
   
5.     #    return 0
   
6.     #end
   
7.     global ns1 = n - 1
   
8.     f = factor(Vector, ns1)
   
9.     for i in f
   
10.         if powermod(b, div(ns1, i), n) == 1
    
11.             global ns1 = div(ns1, i)
    
12.         end
    
13.     end
    
14.     ns1
    
15. end
    
16. 
    
17. println(ord(14951, 2))
    
18. println(ord(3391, 2))
    

复制代码

这个函数返回以b为底的n的阶，就是满足 b^d=1 mod n的最小的d

.·.·.

刚刚发现
  isprime(x::BigInt, [reps = 25]) -> Bool

  Probabilistic primality test. Returns true if x is prime with high probability (pseudoprime); and false if x is composite (not prime). The false
  positive rate is about 0.25^reps. reps = 25 is considered safe for cryptographic applications (Knuth, Seminumerical Algorithms).

  julia> isprime(big(3))
  true
isprime是一种概率算法
而且

1. @time for i in range(1, step = 2, stop = 9999999)
   
2.            test = b + i
   
3.            if fermatProbablePrimeTest2(test)
   
4.                if millerRabinProbablePrimeTest(test, 2)
   
5.                    if false
   
6.                        println("$test is SPSP(2)")
   
7.                    end
   
8.                else
   
9.                    println("$test is PSP(2)")
   
10.                end
    
11.            end
    
12.        end

复制代码

10.750330 seconds (87.82 M allocations: 1.502 GiB, 9.61% gc time)

如果不使用isprime的话，计算速度会有一定提升

.·.·.

无心人 发表于 2020-11-28 14:32
用fermat的原因是想测试下PSP2基2伪素数的概率，另外，打算用frobenius测试做最终判定是不是素数

忽然发现你的fermat测试算法不纯粹（而且偏慢）
直接

1. function fermatProbablePrimeTest2(n)
   
2.     powermod(2, n-1, n) == 1
   
3. end

复制代码

会快一点，而且可以发现更多的伪素数

.·.·.

试图Miller-Rabin还是不如直接Fermat快。
现在fermat部分应该已经优化到头了。虽然都是调包侠，但明显Rust比Julia调包更快一些
睡了睡了~

1. use rug::{Assign, Integer};
   
2. fn main() {
   
3.     let now=std::time::Instant::now();
   
4.     let mut a=Integer::from(10000000000000000000000000001u128);
   
5.     let mut am1=Integer::from(10000000000000000000000000000u128);
   
6.     let mut b=Integer::from(2);
   
7.     let mut prime=0u64;
   
8.     for _i in 0..5000000{
   
9.         prime+={b.assign(2u8);let _=b.pow_mod_mut(&am1,&a);b==1u8} as u64;
   
10.         a+=2;
    
11.         am1+=2;
    
12.     }
    
13.     println!("{},{},{:?}",a,prime,now.elapsed());
    
14. }

复制代码

（加上isprobableprime的除非改并行否则打死不改版（改并行时候需要注意结果循环外的变量不能被共享））

1. /*use num_bigint::BigUint;
   
2. fn main() {
   
3.     let now=std::time::Instant::now();
   
4.     let mut a=BigUint::from(100000000000000u64);
   
5.     a=&a*&a+BigUint::from(1u8);
   
6.     let b:BigUint=BigUint::from(2u8);
   
7.     //let mut c:BigUint=BigUint::from(0u8);
   
8.     let mut prime=0u64;
   
9.     for _i in 0..500000{
   
10.         if b.modpow(&(&a-BigUint::from(1u8)),&a)==BigUint::from(1u8){prime+=1}
    
11.         a+=&b;
    
12.     }
    
13.     println!("{},{},{:?}",a,prime,now.elapsed());
    
14. }*/
    
15. /*
    
16. use bigint::U256;
    
17. fn fermat(module:U256)->bool{//module should be larger than 1<<64 and less than 1<<128
    
18.     let [mut a,b,_,_]=module.0;
    
19.     a-=1;
    
20.     let mut res=U256([1,0,0,0]);
    
21.     let mut pow=U256([2,0,0,0]);
    
22.     for _ in 0..64{
    
23.         if a&1==1{
    
24.             res=(res*pow)%module;
    
25.             pow=(pow*pow)%module;
    
26.         }else{
    
27.             pow=(pow*pow)%module;
    
28.         }
    
29.         a>>=1
    
30.     }
    
31.     a=b;
    
32.     while a>1{
    
33.         if a&1==1{
    
34.             res=(res*pow)%module;
    
35.             pow=(pow*pow)%module;
    
36.         }else{
    
37.             pow=(pow*pow)%module;
    
38.         }
    
39.         a>>=1
    
40.     }
    
41.     (res*pow)%module==U256([1,0,0,0])
    
42. }
    
43. fn main() {
    
44.     let now=std::time::Instant::now();
    
45.     let mut a=U256::exp10(28)+U256::from(1);
    
46.     let mut prime=0u32;
    
47.     for _ in 0..50000{
    
48.         if fermat(a){
    
49.             prime+=1;
    
50.         }
    
51.         (a.0)[0]+=2;
    
52.         if (a.0)[0]<2{(a.0)[1]+=1}
    
53.     }
    
54.     println!("{},{},{:?}",a,prime,now.elapsed());
    
55. }*/
    
56. /*use ramp::Int;
    
57. fn main() {
    
58.     let now=std::time::Instant::now();
    
59.     let mut a=Int::from(10000000000000000000000000001u128);
    
60.     let b=Int::from(2u8);
    
61.     let mut c=Int::from(0u8);
    
62.     for _i in 0..500000{
    
63.         c+=b.pow_mod(&(&a-1),&a);
    
64.         a=a+2;
    
65.     }
    
66.     println!("{},{},{:?}",a,c,now.elapsed());
    
67. }*/
    
68. /*
    
69. use gmp::mpz::Mpz;
    
70. fn main() {
    
71.     let now=std::time::Instant::now();
    
72.     let mut a=Mpz::from_str_radix("10000000000000000000000000001",10).unwrap();
    
73.     let b=Mpz::from_str_radix("2",10).unwrap();
    
74.     let mut prime=0u64;
    
75.     for _i in 0..5000000{
    
76.         prime+=(b.powm(&(&a-1),&a)-1).is_zero() as u64;
    
77.         a=a+2;
    
78.     }
    
79.     println!("{},{},{:?}",a,prime,now.elapsed());
    
80. }*/
    
81. /*Pure Miller-Rabin
    
82. 
    
83. use rug::{Assign, Integer};
    
84. fn main() {
    
85.     let now=std::time::Instant::now();
    
86.     let mut a=Integer::from(10000000000000000000000000001u128);
    
87.     let mut am1=Integer::from(10000000000000000000000000000u128);
    
88.     let mut b=Integer::from(2);
    
89.     let mut buffer=Integer::from(0);
    
90.     let mut am1d=Integer::from(0);
    
91.     let mut prime=0u64;
    
92.     let mut psp=0u64;
    
93.     for _i in 0..5000000{
    
94.         prime+={b.assign(2u8);//Miller-Rabin
    
95.                 let aa=a.find_one(0).unwrap();
    
96.                 am1d.assign(&am1>>aa);
    
97.                 let _=b.pow_mod_mut(&am1d,&a);
    
98.                 (b==1 || b==am1 )||{
    
99.                     let mut f=false;
    
100.                     for _ in 0..aa{
     
101.                         let _=b.pow_mod_mut(&Integer::from(2),&a);
     
102.                         if b==1{break}else if b==am1{f=true;break}
     
103.                     }
     
104.                     f
     
105.                 }
     
106.             } as u64;
     
107.         a+=2;
     
108.         am1+=2;
     
109.     }
     
110.     println!("{},{},{:?}",a,prime-psp,now.elapsed());
     
111. }
     
112. */
     
113. use rug::{Assign, Integer};
     
114. fn main() {
     
115.     const CONSTARR:[u128;7]=[1u128<<1,1u128<<2,1u128<<4,1u128<<8,1u128<<16,1u128<<32,1u128<<64];
     
116.     let now=std::time::Instant::now();
     
117.     let mut a=Integer::from(10000000000000000000000000001u128);
     
118.     let mut am1=Integer::from(10000000000000000000000000000u128);
     
119.     let mut b=Integer::from(2);
     
120.     let mut am1d=Integer::from(0);
     
121.     let mut prime=0u64;
     
122.     let mut psp=0u64;
     
123.     for _i in 0..50000000{
     
124.         let aa=am1.find_one(0).unwrap().min(6);
     
125.         am1d.assign(&am1>>aa);
     
126.         if {b.assign(CONSTARR[aa as usize]);let _=b.pow_mod_mut(&am1d,&a);b==1u8}{//Fermat
     
127.             psp+=1;
     
128.             if {b.assign(2u8);//Miller-Rabin
     
129.                 let aa=a.find_one(0).unwrap();
     
130.                 am1d.assign(&am1>>aa);
     
131.                 let _=b.pow_mod_mut(&am1d,&a);
     
132.                 (b==1 || b==am1 )||{
     
133.                     let mut f=false;
     
134.                     for _ in 0..aa{
     
135.                         let _=b.pow_mod_mut(&Integer::from(2),&a);
     
136.                         if b==am1{f=true;break}
     
137.                     }
     
138.                     f
     
139.                 }
     
140.             }{
     
141.                 if a.is_probably_prime(25)==rug::integer::IsPrime::No{
     
142.                     println!("{} is SPSP(2)",a)
     
143.                 }else{prime+=1}
     
144.             }else{println!("{} is PSP(2)",a)}
     
145.         }
     
146.         a+=2;
     
147.         am1+=2;
     
148.     }
     
149.     println!("{}.{},{},{:?}",a,prime,prime-psp,now.elapsed());
     
150. }

复制代码

不知道为什么，Rust快到没朋友
（我这几个小时一直在Rust慢得怀疑人生和Rust快得没朋友之间徘徊）

另，这是Rust比Julia快的地方，只要你的Julia版本高于0.6，请务必执行

1. chmod 644 ~/.julia/packages/Primes/eRBAQ/src/Primes.jl
   
2. sed -i 's/(Any, Cint)/(Ref{BigInt}, Cint)/g' ~/.julia/packages/Primes/eRBAQ/src/Primes.jl
   
3. chmod 444 ~/.julia/packages/Primes/eRBAQ/src/Primes.jl

复制代码

我也不知道为什么能从30s直接变成11s
或许这就是人生吧

最终测速结果：
11.642703 seconds (63.51 M allocations: 1.428 GiB, 5.59% gc time)
（只算Fermat  9.297450 seconds (59.88 M allocations: 1.007 GiB, 7.64% gc time)）
Rust
10000000000000000000100000001.1552676,0,82.901862532s
（搜了前10^8个数字一个PSP都没有真的没问题吗？）

uk702

大牛啊！这都能发现。可以向 julia 发个补丁，向他们确认一下。

1. 另，这是Rust比Julia快的地方，只要你的Julia版本高于0.6，请务必执行
   
2. chmod 644 ~/.julia/packages/Primes/eRBAQ/src/Primes.jl
   
3. sed -i 's/(Any, Cint)/(Ref{BigInt}, Cint)/g' ~/.julia/packages/Primes/eRBAQ/src/Primes.jl
   
4. chmod 444 ~/.julia/packages/Primes/eRBAQ/src/Primes.jl
   
5. 复制代码
   

复制代码

mathematica

.·.·. 发表于 2020-11-30 02:54
试图Miller-Rabin还是不如直接Fermat快。
现在fermat部分应该已经优化到头了。虽然都是调包侠，但明显Rust ...

来这边搞搞马青公式：
https://bbs.emath.ac.cn/thread-17590-1-1.html

.·.·.

无心人 发表于 2020-11-29 15:15
这个函数返回以b为底的n的次数，就是满足 b^d=1 mod n的最小的d

信了你的邪……

1. from=10000000000000000050000000001,len=10000000000
   
2. 10000000000000000060000000001.155106638,0,10335.762851883s
   

复制代码

1. from=10000000000000000040000000001,len=10000000000
   
2. 10000000000000000050000000001.155086324,0,10317.813073111s
   

复制代码

1. from=10000000000000000030000000001,len=10000000000
   
2. 10000000000000000040000000001.155097227,0,10329.304771961s
   

复制代码

1. from=10000000000000000020000000001,len=10000000000
   
2. 10000000000000000030000000001.155112727,0,10291.313970292s
   

复制代码

1. from=10000000000000000010000000001,len=10000000000
   
2. 10000000000000000020000000001.155103639,0,10312.597111373s
   

复制代码

1. from=10000000000000000000000000001,len=10000000000
   
2. 10000000000000000010000000001.155111392,0,10325.273760498s
   

复制代码

mathematica

伪素数是无限的，你想怎么办？
你做这个的意义是什么？