853 Pisano Periods 1 -- v2.pl 1.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071
  1. #!/usr/bin/perl
  2. # Pisano Periods 1
  3. # https://projecteuler.net/problem=853
  4. # Runtime: 0.345s
  5. use 5.036;
  6. use List::Util qw(first);
  7. use ntheory qw(:all);
  8. sub fibonacci ($n) {
  9. lucasu(1, -1, $n);
  10. }
  11. sub fibmod ($d, $n) {
  12. (lucas_sequence($n, 1, -1, $d))[0];
  13. }
  14. sub divisors_le ($n, $k) {
  15. my @d = (1);
  16. my @pp = grep { $_->[0] <= $k } factor_exp($n);
  17. foreach my $pp (@pp) {
  18. my ($p, $e) = @$pp;
  19. my @t;
  20. my $r = 1;
  21. for my $i (1 .. $e) {
  22. $r *= $p;
  23. foreach my $u (@d) {
  24. push(@t, $u * $r) if ($u * $r <= $k);
  25. }
  26. }
  27. push @d, @t;
  28. }
  29. return @d;
  30. }
  31. sub pisano_period_pp ($p, $k = 1) {
  32. $p**($k - 1) * first { fibmod($_, $p) == 0 } divisors($p - kronecker($p, 5));
  33. }
  34. sub my_pisano_period ($n) {
  35. return 0 if ($n <= 0);
  36. return 1 if ($n == 1);
  37. my $d = lcm(map { pisano_period_pp($_->[0], $_->[1]) } factor_exp($n));
  38. foreach my $k (0 .. 2) {
  39. my $t = $d << $k;
  40. if (fibmod($t + 1, $n) == 1) {
  41. return $t;
  42. }
  43. }
  44. die "Error for n = $n";
  45. }
  46. my $n = 120;
  47. my $limit = 1e9;
  48. say vecsum(grep { my_pisano_period($_) == $n } divisors_le(fibonacci($n), $limit));