-
Notifications
You must be signed in to change notification settings - Fork 6
/
853 Pisano Periods 1 -- v2.pl
70 lines (48 loc) · 1.21 KB
/
853 Pisano Periods 1 -- v2.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#!/usr/bin/perl
# Pisano Periods 1
# https://projecteuler.net/problem=853
# Runtime: 0.345s
use 5.036;
use List::Util qw(first);
use ntheory qw(:all);
sub fibonacci ($n) {
lucasu(1, -1, $n);
}
sub fibmod ($d, $n) {
(lucas_sequence($n, 1, -1, $d))[0];
}
sub divisors_le ($n, $k) {
my @d = (1);
my @pp = grep { $_->[0] <= $k } factor_exp($n);
foreach my $pp (@pp) {
my ($p, $e) = @$pp;
my @t;
my $r = 1;
for my $i (1 .. $e) {
$r *= $p;
foreach my $u (@d) {
push(@t, $u * $r) if ($u * $r <= $k);
}
}
push @d, @t;
}
return @d;
}
sub pisano_period_pp ($p, $k = 1) {
$p**($k - 1) * first { fibmod($_, $p) == 0 } divisors($p - kronecker($p, 5));
}
sub my_pisano_period ($n) {
return 0 if ($n <= 0);
return 1 if ($n == 1);
my $d = lcm(map { pisano_period_pp($_->[0], $_->[1]) } factor_exp($n));
foreach my $k (0 .. 2) {
my $t = $d << $k;
if (fibmod($t + 1, $n) == 1) {
return $t;
}
}
die "Error for n = $n";
}
my $n = 120;
my $limit = 1e9;
say vecsum(grep { my_pisano_period($_) == $n } divisors_le(fibonacci($n), $limit));