fermat_factorization_method_2.pl 1.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869
  1. #!/usr/bin/perl
  2. # Daniel "Trizen" Șuteu
  3. # Date: 13 September 2017
  4. # https://github.com/trizen
  5. # Fermat's factorization method.
  6. # Theorem:
  7. # If the absolute difference between the prime factors of a
  8. # semiprime `n` is known, then `n` can be factored in polynomial time.
  9. # Based on the following quadratic equation:
  10. # x^2 + (a - b)*x - a*b = 0
  11. #
  12. # which has the solutions:
  13. # x₁ = -a
  14. # x₂ = +b
  15. # See also:
  16. # https://en.wikipedia.org/wiki/Fermat%27s_factorization_method
  17. use 5.020;
  18. use strict;
  19. use warnings;
  20. use experimental qw(signatures);
  21. use ntheory qw(vecprod sqrtint is_prime is_square valuation);
  22. sub fermat_factorization ($n) {
  23. # Check for primes and negative numbers
  24. return () if ($n <= 1);
  25. return ($n) if is_prime($n);
  26. # Check for divisibility by 2
  27. if (!($n & 1)) {
  28. my $v = valuation($n, 2);
  29. return ((2) x $v, __SUB__->($n >> $v));
  30. }
  31. my $p = sqrtint($n);
  32. my $q = $p * $p - $n;
  33. until (is_square($q)) {
  34. $q += 2 * $p++ + 1;
  35. }
  36. my $s = sqrtint($q);
  37. my ($x1, $x2) = (
  38. ($p + $s),
  39. ($p - $s),
  40. );
  41. return sort { $a <=> $b } (
  42. __SUB__->($x1),
  43. __SUB__->($x2)
  44. );
  45. }
  46. foreach my $n (160587846247027, 5040, 65127835124, 6469693230) {
  47. my @f = fermat_factorization($n);
  48. say join(' * ', @f), " = $n";
  49. die 'error' if vecprod(@f) != $n;
  50. }