machine-prime 1.5.7

use machine_prime::{
    is_prime,strong_fermat_128,
    one_mont_128,two_mont_128,
    to_mont_128,lucas_128,
    mul_inv2_128,PRIME_TABLE_128
    };

/*

  Demonstration of an almost surely correct deterministic primality test for n < 2^128, that is also faster than ECPP
  
  Argument:  All semiprimes of the form (2x+1)(4x+1) are eliminated by the fermat witnesses. 
  There are probably no Carmichael numbers of the form p,q,r 3 mod 4 considering that they are much rarer than 
  the semiprimes. Likewise we expect no other composites to pass the fermat component, the most probable candidate
  is of the form (4x+1)(12x+1), with a probability of existing less than 0.065. 
  Additionally there are no known pseudoprimes to 2 and the lucas test under 2^64. The few pseudoprimes to both a
  lucas test and a strong fermat test, occur at a rate of approximately 1/42000 per fermat pseudoprime. Therefore
  in the worst case we can estimate a counterexample exists with a probability of less than 1.5E-6.  
  
  
  
  Note that running this example requires using the "internal","lucas" and "wide" features
  
  Approximate run time 21.5 Fermat tests
  
  This is probably superseded by the QFT feature which provides a test with complexity of about 7.5 fermat tests
  while having a much stronger argument for infallibility than the BPSW.  
*/

fn strong_test(x: u128) -> bool{
   // We cannot use this algorithm for n < 2^64 due to a Montgomery transform "error" 
   // when using 64-bit inputs for a 128-bit transform. Machine-prime's 64-bit is_prime is also faster 
   // and  proven to be correct
    if x < 1u128<<64{
       return is_prime(x as u64);
    }
    // X must be odd to ensure that the inverse over 2^128 can be calculated
    if x&1==0{
       return false;
    }
    
    let mut idx = 0usize;
   // Let's perform trial division for the average case
    while idx < 256{
    
      let prod = x.wrapping_mul(PRIME_TABLE_128[idx]);
      
      if prod <= PRIME_TABLE_128[idx+1]{
         return prod==1;  
      }
      idx+=2;
    }
   // Initialise necessary values
   
    let inv = mul_inv2_128(x);

    let tzc = x.wrapping_sub(1).trailing_zeros();
    let one = one_mont_128(x);
    let oneinv = x.wrapping_sub(one);
    let two = two_mont_128(one, x);

   // Witness using 2. Nearly all composites will be eliminated by this
   // You can remove this if it is being used as a complement to the is_prime_128 library function
    if !strong_fermat_128(x, tzc, two, one, oneinv, inv) {
        return false;
    }

   // Fermat component alone should be sufficient
   // The last 5 witnesses were computed from 59000 pseudoprimes to the prime witnesses [2;47]
   for i in [3,5,7,11,13,17,19,23,29,31,37,41,43,47,511,659,679,8129,70157]{
       let montwitness = to_mont_128(i,x);
       // reuse the prior values from the witness-2 check
       if !strong_fermat_128(x,tzc,montwitness,one,oneinv,inv){
          return false;
       }
   }
   // It is possible that there exist wieferich primes to the witnesses above
   // meaning that a perfect square with no nonquadratic residues may exist, resulting in the Lucas/QFT parameter
   // searches running infinitely.
   // However, this is almost certainly false. No wieferich primes to 2 exist between 2^64 and 2^92, let alone to 
   // 16 other primes in addition.If they do then we can simply check if x == sqrt(x)^2.  
   
   // Add a Lucas sequence test 
    lucas_128(x,one,two,inv)
   // Or conversely a QFT
    // qft(x,one,two,oneinv,inv)
}

fn main(){
  for i in 0..100_000{
    if strong_test(u128::MAX-i){
       println!("2^128-{} is definitely prime",i+1);
    }
  }
}