modular_fibonacci_cassini.pl 1.3 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455
  1. #!/usr/bin/perl
  2. # An efficient algorithm for computing the nth-Fibonacci number (mod m).
  3. # Algorithm from:
  4. # https://metacpan.org/source/KRYDE/Math-NumSeq-72/lib/Math/NumSeq/Fibonacci.pm
  5. # See also:
  6. # https://en.wikipedia.org/wiki/Fibonacci_number
  7. use 5.020;
  8. use warnings;
  9. use experimental qw(signatures);
  10. use Math::GMPz;
  11. use Math::Prime::Util::GMP qw(consecutive_integer_lcm gcd);
  12. sub fibmod ($n, $m) {
  13. $n = Math::GMPz->new("$n");
  14. $m = Math::GMPz->new("$m");
  15. my ($f, $g, $a) = (0, 1, -2);
  16. foreach my $bit (split(//, substr(Math::GMPz::Rmpz_get_str($n, 2), 1))) {
  17. ($g *= $g) %= $m;
  18. ($f *= $f) %= $m;
  19. my $t = ($g << 2) - $f + $a;
  20. $f += $g;
  21. if ($bit) {
  22. ($f, $g, $a) = ($t - $f, $t, -2);
  23. }
  24. else {
  25. ($g, $a) = ($t - $f, 2);
  26. }
  27. }
  28. return ($g % $m);
  29. }
  30. sub fibonacci_factorization ($n, $B = 10000) {
  31. my $k = consecutive_integer_lcm($B); # lcm(1..B)
  32. my $F = fibmod($k, $n); # Fibonacci(k) (mod n)
  33. return gcd($F, $n);
  34. }
  35. say fibonacci_factorization("121095274043", 700); #=> 470783 (p+1 is 700-smooth)
  36. say fibonacci_factorization("544812320889004864776853", 3000); #=> 333732865481 (p-1 is 3000-smooth)