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

mathematica 发表于 2020-11-30 15:40:30

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

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

无心人 发表于 2020-11-30 17:04:13

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

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

.·.·. 发表于 2020-11-30 17:06:46

离成功优化powermod大概只差一步了fn pow_mod(mut base:u128,mut exp:u128,module:u128)->u128{
    // assert!(module < 1<<94 && module&1==1 && exp>0)
    let mut acc = 1;
    //let m64=(module>>64) as u64;
    let m126=(1<<126)%module;
    let m157=(m126<<31)%module;
    let mulmod=|mut x,mut y|{
      let =;
      x=lx*hy+ly*hx;
      y=lx*ly+((x&9223372036854775807)<<63);x>>=63;//y=low_mid, x=high_mid<<126
      //lx*ly+y+x<<126+hx*hy<<126
      x+=hx*hy;
      //lx*ly+y+x<<126
      hx=x&2147483647;
      x>>=31;
      //(y+hx*m126+x*m157)%module
      y+=hx*m126+x*m157;
      //y-(((y>>64)as u64)/m64) as u128*module
      y%module
      //lx.max(ly) < 1<<63 ==> lx*ly < 1<<126
      // y < 9223372036854775807<<63 < 1<<126
      //hx <2147483647*module< 1<<31+94 = 1<<126
      // x <(module*module>>(63*2+31))*module < 1<<(94*2-63*2-31+94)=1<<125
      //Thus, lx*ly+y+hx*m126+x*m157 < 1<<128, safe to add them together and executing % ops.
    };
    while exp > 1 {
      if (exp & 1) == 1 {
            acc = mulmod(acc,base);
            base = mulmod(base,base);
      }else{
            base = mulmod(base,base);
      }
      exp >>= 1;
    }
    mulmod(acc,base)
}
fn main(){
    const CONSTARR:=;
    let now=std::time::Instant::now();
    let arguments: Vec<String> = std::env::args().collect();
    let begin=arguments.get(2).and_then(|x|Some(x.parse::<u128>().unwrap_or(0))).unwrap_or(0);
    let len=arguments.get(1).and_then(|x|Some(x.parse::<u128>().unwrap_or(10000000))).unwrap_or(10000000);
    let mut am1=10000000000000000000000000000u128+begin*len;
    let mut a=10000000000000000000000000001u128+begin*len;
    println!("from={},len={}",a,len);
    let mut b=2;
    let mut am1d=0;
    let mut prime=0u64;
    let mut psp=0u64;
    for _i in 0..len>>1{
      let aa=am1.trailing_zeros().min(6);
      am1d=&am1>>aa;
      if {b=CONSTARR; pow_mod(b,am1d,a)==1u128}{//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,psp-prime,now.elapsed());
}这个程序是gmp的1/2，但如果拿掉取模步，速度是gmp的两倍多。
另外，或许应该直接写C调GMP，GMP写了Miller-Rabin测试，不知速度如何，但小改一下应该比Fermat快。
毕竟Miller-Rabin可以少算好几次乘法。

.·.·. 发表于 2020-11-30 18:13:22

本帖最后由 .·.·. 于 2020-11-30 19:06 编辑

无心人 发表于 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

#不使用Primes，直接手写一个isprime
isprime(x::BigInt) = ccall((:__gmpz_probab_prime_p,:libgmp),Cint, (Ref{BigInt}, Cint), x, 25) > 0
const b = big(10)^28

function fermatProbablePrimeTest2(n)
    a = powermod(2, n-1, n)
    return a == 1
end

function millerRabinProbablePrimeTest( n, b )
    k = trailing_zeros(n-1)
    m = (n-1) >> k
    a = powermod(b, m, n)
    if (a == 1) || (a == (n-1))
      return true
    else
      for i in 1:k-1
            a = powermod(a, 2, n)
            if a == n - 1
                return true
            end
      end
    end
    return false
end
@time begin
    Threads.@threads for i in range(1, step = 2, stop = 9999999)
      test = b + i
      if fermatProbablePrimeTest2(test)
            if millerRabinProbablePrimeTest(test, 2)
                if !isprime(test)
                  println("$test is SPSP(2)")
                end
            else
                println("$test is PSP(2)")
            end
      end
    end
end

无心人 发表于 2020-11-30 21:59:19

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

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

.·.·. 发表于 2020-11-30 22:40:16

本帖最后由 .·.·. 于 2020-12-1 00:56 编辑

无心人 发表于 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
这一行可以带来非常可怕的性能提升

无心人 发表于 2020-12-1 10:16:49

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

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

1、设定初始值为base，则s位置代表的整数是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 = 0
                start=start+2*p
                j=j+2
               

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

无心人 发表于 2020-12-1 14:43:10

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

.·.·. 发表于 2020-12-1 18:54:38

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

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

无心人 发表于 2020-12-1 19:00:06

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

我感觉模乘不需要除法

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