view release on metacpan or  search on metacpan

lib/Math/Prime/Util/RandomPrimes.pm  view on Meta::CPAN

  my($status,$prime,$prime_seed,$prime_gen_counter,$cert)
     = _ST_Random_prime($k, $seed);
  croak "Shawe-Taylor random prime failure" unless $status;
  croak "Shawe-Taylor random prime failure: prime $prime failed certificate verification!" unless verify_prime($cert);

  return ($prime,$cert);
}

sub _seed_plus_one {
    my($s) = @_;
    for (my $i = length($s)-1; $i >= 0; $i--) {
        vec($s, $i, 8)++;
        last unless vec($s, $i, 8) == 0;
    }
    return $s;
}

sub _ST_Random_prime {  # From FIPS 186-4
  my($k, $input_seed) = @_;
  croak "Shawe-Taylor random prime must have length >= 2"  if $k < 2;
  $k = int("$k");

  croak "Shawe-Taylor random prime, invalid input seed"
     unless defined $input_seed && length($input_seed) >= 32;

  if (!defined $Digest::SHA::VERSION) {
    eval { require Digest::SHA;
           my $version = $Digest::SHA::VERSION;
           $version =~ s/[^\d.]//g;
           $version >= 4.00; }
      or do { croak "Must have Digest::SHA 4.00 or later"; };
  }

  my $k2 = tobigint(powint(2,$k-1));    # $k2 is a bigint

  if ($k < 33) {
    my $seed = $input_seed;
    my $prime_gen_counter = 0;
    my $kmask    = 0xFFFFFFFF >> (32-$k);    # Does the mod operation
    my $kstencil = (1 << ($k-1)) | 1;        # Sets high and low bits
    while (1) {
      my $seedp1 = _seed_plus_one($seed);
      my $cvec = Digest::SHA::sha256($seed) ^ Digest::SHA::sha256($seedp1);
      # my $c = Math::BigInt->from_hex('0x' . unpack("H*", $cvec));
      # $c = $k2 + ($c % $k2);
      # $c = (2 * ($c >> 1)) + 1;
      my($c) = unpack("N*", substr($cvec,-4,4));
      $c = ($c & $kmask) | $kstencil;
      $prime_gen_counter++;
      $seed = _seed_plus_one($seedp1);
      my ($isp, $cert) = is_provable_prime_with_cert($c);
      return (1,$c,$seed,$prime_gen_counter,$cert) if $isp;
      return (0,0,0,0) if $prime_gen_counter > 10000 + 16*$k;
    }
  }
  my($status,$c0,$seed,$prime_gen_counter,$cert)
     = _ST_Random_prime( (($k+1)>>1)+1, $input_seed);
  return (0,0,0,0) unless $status;
  $cert = cmpint($c0,"18446744073709551615") <= 0
         ?  ""  :  _strip_proof_header($cert);
  my $iterations = int(($k + 255) / 256) - 1;  # SHA256 generates 256 bits
  my $old_counter = $prime_gen_counter;
  my $c02 = lshiftint($c0);  # $c02 = 2*$c0
  my $xstr = '';
  for my $i (0 .. $iterations) {
    $xstr = Digest::SHA::sha256_hex($seed) . $xstr;
    $seed = _seed_plus_one($seed);
  }
  my $x = fromdigits($xstr,16);
  $x = $k2 + ($x % $k2);
  my $t = cdivint($x, $c02);
  _make_big_gcds() if $_big_gcd_use < 0;
  while (1) {
    my $c = add1int(mulint($t,$c02));
    if ($c > 2*$k2) {
      $t = cdivint($k2, $c02);
      $c = add1int(mulint($t,$c02));
    }
    $prime_gen_counter++;

    # Don't do the Pocklington check unless the candidate looks prime
    my $looks_prime = 0;
    if (MPU_USE_GMP) {
      # MPU::GMP::is_prob_prime has fast tests built in.
      $looks_prime = Math::Prime::Util::GMP::is_prob_prime($c);
    } else {
      # No GMP, so first do trial divisions, then a SPSP test.
      $looks_prime = gcd($c, 111546435) == 1;
      if ($looks_prime && $_big_gcd_use && $c > $_big_gcd_top) {
        $looks_prime = gcd($c, $_big_gcd[0]) == 1 &&
                       gcd($c, $_big_gcd[1]) == 1 &&
                       gcd($c, $_big_gcd[2]) == 1 &&
                       gcd($c, $_big_gcd[3]) == 1;
      }
      $looks_prime = 0 if $looks_prime && !is_strong_pseudoprime($c, 3);
    }

    if ($looks_prime) {
      # We could use a in (2,3,5,7,11,13), but pedantically use FIPS 186-4.
      my $astr = '';
      for my $i (0 .. $iterations) {
        $astr = Digest::SHA::sha256_hex($seed) . $astr;
        $seed = _seed_plus_one($seed);
      }
      my $a = fromdigits($astr,16);
      $a = addint(modint($a,$c-3),2);
      my $z = powmod($a, lshiftint($t), $c);
      if (gcd($z-1,$c) == 1 && powmod($z, $c0, $c) == 1) {
        croak "Shawe-Taylor random prime failure at ($k): $c not prime"
          unless is_prob_prime($c);
        $cert = "[MPU - Primality Certificate]\nVersion 1.0\n\n" .
                 "Proof for:\nN $c\n\n" .
                 "Type Pocklington\nN $c\nQ $c0\nA $a\n" .
                 $cert;
        return (1, $c, $seed, $prime_gen_counter, $cert);
      }
    } else {
      # Update seed "as if" we performed the Pocklington check from FIPS 186-4
      for my $i (0 .. $iterations) {
        $seed = _seed_plus_one($seed);
      }
    }
    return (0,0,0,0) if $prime_gen_counter > 10000 + 16*$k + $old_counter;
    $t++;
  }
}

sub random_safe_prime {
  my $bits = int("$_[0]");
  croak "random_safe_prime, bits must be >= 3" unless $bits >= 3;
  return (5,7)[urandomb(1)] if $bits == 3;
  return 11 if $bits == 4;
  return 23 if $bits == 5;
  return (47,59)[urandomb(1)] if $bits == 6;
  return (83,107)[urandomb(1)] if $bits == 7;

  # Without GMP (e.g. Calc), this can be significantly faster.
  # With GMP, they are about the same.
  return _random_safe_prime_large($bits) if $bits > 35;

  my($p,$q);
  while (1) {
    $q = Math::Prime::Util::random_nbit_prime($bits-1);
    my $qm = modint($q, 1155);  # a nice native int
    next if ($qm % 3) != 2
         || ($qm % 5) == 2 || ($qm % 7) == 3 || ($qm % 11) == 5;
    $p = mulint(2, $q) + 1;
    # This is sufficient, but we'll do the full test including pre-tests.
    #last if is_pseudoprime($p,2);  # p is prime if q is prime
    last if is_prob_prime($p);
  }
  return $p;
}

sub _random_safe_prime_large {
  my $bits = shift;
  croak "Not enough bits for large random_safe_prime" if $bits <= 35;

  # Set first and last two bits
  my $base = addint(lshiftint(1, $bits-1),3);
  # Fill in lower portion with random bits, leaving 32 upper
  $base = addint($base, lshiftint(urandomb($bits - 35), 2));

  while (1) {
    my($p,$q,$qmod,$pmod);
    # 1. generate random nbit p
    $p = lshiftint(urandomb(32), $bits-33);
    $p = addint($base, $p);

    # 2. p = 2q+1  =>  q = p>>1
    $q = rshiftint($p);

    # 3. Force q mod 6 = 5
    $qmod = modint(add1int($q),6);
    if ($qmod > 0) {
      $q = subint($q,$qmod);
      $q = addint($q,12) if 1+logint($q,2) != $bits-1;
      $p = add1int(lshiftint($q));
    }