论坛 › 难题征解 › 10^28开始的10^14个整数的素性概率性测试算法实践报告

无心人 发表于 2020-11-28 14:32:43

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


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

.·.·. 发表于 2020-11-28 14:38:57

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

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

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

无心人 发表于 2020-11-29 15:15:32

#找到最小的d满足b^d = 1 (mod n)
function ord(n, b)
    #if !isprime(n)
    #    return 0
    #end
    global ns1 = n - 1
    f = factor(Vector, ns1)
    for i in f
      if powermod(b, div(ns1, i), n) == 1
            global ns1 = div(ns1, i)
      end
    end
    ns1
end

println(ord(14951, 2))
println(ord(3391, 2))


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

.·.·. 发表于 2020-11-29 18:48:40

本帖最后由 .·.·. 于 2020-11-29 23:07 编辑

刚刚发现
isprime(x::BigInt, ) -> 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是一种概率算法
而且

@time for i in range(1, step = 2, stop = 9999999)
         test = b + i
         if fermatProbablePrimeTest2(test)
               if millerRabinProbablePrimeTest(test, 2)
                   if false
                     println("$test is SPSP(2)")
                   end
               else
                   println("$test is PSP(2)")
               end
         end
       end

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

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

.·.·. 发表于 2020-11-29 22:41:04

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

忽然发现你的fermat测试算法不纯粹（而且偏慢）
直接function fermatProbablePrimeTest2(n)
    powermod(2, n-1, n) == 1
end
会快一点，而且可以发现更多的伪素数

.·.·. 发表于 2020-11-30 02:54:34

本帖最后由 .·.·. 于 2020-11-30 04:20 编辑

试图Miller-Rabin还是不如直接Fermat快。
现在fermat部分应该已经优化到头了。虽然都是调包侠，但明显Rust比Julia调包更快一些
睡了睡了~
use rug::{Assign, Integer};
fn main() {
    let now=std::time::Instant::now();
    let mut a=Integer::from(10000000000000000000000000001u128);
    let mut am1=Integer::from(10000000000000000000000000000u128);
    let mut b=Integer::from(2);
    let mut prime=0u64;
    for _i in 0..5000000{
      prime+={b.assign(2u8);let _=b.pow_mod_mut(&am1,&a);b==1u8} as u64;
      a+=2;
      am1+=2;
    }
    println!("{},{},{:?}",a,prime,now.elapsed());
}（加上isprobableprime的除非改并行否则打死不改版（改并行时候需要注意结果循环外的变量不能被共享））
/*use num_bigint::BigUint;
fn main() {
    let now=std::time::Instant::now();
    let mut a=BigUint::from(100000000000000u64);
    a=&a*&a+BigUint::from(1u8);
    let b:BigUint=BigUint::from(2u8);
    //let mut c:BigUint=BigUint::from(0u8);
    let mut prime=0u64;
    for _i in 0..500000{
      if b.modpow(&(&a-BigUint::from(1u8)),&a)==BigUint::from(1u8){prime+=1}
      a+=&b;
    }
    println!("{},{},{:?}",a,prime,now.elapsed());
}*/
/*
use bigint::U256;
fn fermat(module:U256)->bool{//module should be larger than 1<<64 and less than 1<<128
    let =module.0;
    a-=1;
    let mut res=U256();
    let mut pow=U256();
    for _ in 0..64{
      if a&1==1{
            res=(res*pow)%module;
            pow=(pow*pow)%module;
      }else{
            pow=(pow*pow)%module;
      }
      a>>=1
    }
    a=b;
    while a>1{
      if a&1==1{
            res=(res*pow)%module;
            pow=(pow*pow)%module;
      }else{
            pow=(pow*pow)%module;
      }
      a>>=1
    }
    (res*pow)%module==U256()
}
fn main() {
    let now=std::time::Instant::now();
    let mut a=U256::exp10(28)+U256::from(1);
    let mut prime=0u32;
    for _ in 0..50000{
      if fermat(a){
            prime+=1;
      }
      (a.0)+=2;
      if (a.0)<2{(a.0)+=1}
    }
    println!("{},{},{:?}",a,prime,now.elapsed());
}*/
/*use ramp::Int;
fn main() {
    let now=std::time::Instant::now();
    let mut a=Int::from(10000000000000000000000000001u128);
    let b=Int::from(2u8);
    let mut c=Int::from(0u8);
    for _i in 0..500000{
      c+=b.pow_mod(&(&a-1),&a);
      a=a+2;
    }
    println!("{},{},{:?}",a,c,now.elapsed());
}*/
/*
use gmp::mpz::Mpz;
fn main() {
    let now=std::time::Instant::now();
    let mut a=Mpz::from_str_radix("10000000000000000000000000001",10).unwrap();
    let b=Mpz::from_str_radix("2",10).unwrap();
    let mut prime=0u64;
    for _i in 0..5000000{
      prime+=(b.powm(&(&a-1),&a)-1).is_zero() as u64;
      a=a+2;
    }
    println!("{},{},{:?}",a,prime,now.elapsed());
}*/
/*Pure Miller-Rabin

use rug::{Assign, Integer};
fn main() {
    let now=std::time::Instant::now();
    let mut a=Integer::from(10000000000000000000000000001u128);
    let mut am1=Integer::from(10000000000000000000000000000u128);
    let mut b=Integer::from(2);
    let mut buffer=Integer::from(0);
    let mut am1d=Integer::from(0);
    let mut prime=0u64;
    let mut psp=0u64;
    for _i in 0..5000000{
      prime+={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==1{break}else if b==am1{f=true;break}
                  }
                  f
                }
            } as u64;
      a+=2;
      am1+=2;
    }
    println!("{},{},{:?}",a,prime-psp,now.elapsed());
}
*/
use rug::{Assign, Integer};
fn main() {
    const CONSTARR:=;
    let now=std::time::Instant::now();
    let mut a=Integer::from(10000000000000000000000000001u128);
    let mut am1=Integer::from(10000000000000000000000000000u128);
    let mut b=Integer::from(2);
    let mut am1d=Integer::from(0);
    let mut prime=0u64;
    let mut psp=0u64;
    for _i in 0..50000000{
      let aa=am1.find_one(0).unwrap().min(6);
      am1d.assign(&am1>>aa);
      if {b.assign(CONSTARR);let _=b.pow_mod_mut(&am1d,&a);b==1u8}{//Fermat
            psp+=1;
            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}
                  }
                  f
                }
            }{
                if a.is_probably_prime(25)==rug::integer::IsPrime::No{
                  println!("{} is SPSP(2)",a)
                }else{prime+=1}
            }else{println!("{} is PSP(2)",a)}
      }
      a+=2;
      am1+=2;
    }
    println!("{}.{},{},{:?}",a,prime,prime-psp,now.elapsed());
}不知道为什么，Rust快到没朋友
（我这几个小时一直在Rust慢得怀疑人生和Rust快得没朋友之间徘徊）

另，这是Rust比Julia快的地方，只要你的Julia版本高于0.6，请务必执行
chmod 644 ~/.julia/packages/Primes/eRBAQ/src/Primes.jl
sed -i 's/(Any, Cint)/(Ref{BigInt}, Cint)/g' ~/.julia/packages/Primes/eRBAQ/src/Primes.jl
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)
（只算Fermat9.297450 seconds (59.88 M allocations: 1.007 GiB, 7.64% gc time)）
Rust
10000000000000000000100000001.1552676,0,82.901862532s
（搜了前10^8个数字一个PSP都没有真的没问题吗？）

uk702 发表于 2020-11-30 07:39:44

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

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

mathematica 发表于 2020-11-30 09:07:19

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

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

.·.·. 发表于 2020-11-30 15:17:26

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

信了你的邪……
from=10000000000000000050000000001,len=10000000000
10000000000000000060000000001.155106638,0,10335.762851883s
from=10000000000000000040000000001,len=10000000000
10000000000000000050000000001.155086324,0,10317.813073111s
from=10000000000000000030000000001,len=10000000000
10000000000000000040000000001.155097227,0,10329.304771961s
from=10000000000000000020000000001,len=10000000000
10000000000000000030000000001.155112727,0,10291.313970292s
from=10000000000000000010000000001,len=10000000000
10000000000000000020000000001.155103639,0,10312.597111373s
from=10000000000000000000000000001,len=10000000000
10000000000000000010000000001.155111392,0,10325.273760498s

mathematica 发表于 2020-11-30 15:39:52

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

查看完整版本: 10^28开始的10^14个整数的素性概率性测试算法实践报告