Pisano period: Difference between revisions
(→{{header|Go}}: Replaced 'isPrime' function with a more efficient one.) |
(Added Sidef) |
||
Line 321: | Line 321: | ||
[3, 8, 6, 20, 24, 16, 12, 24, 60, 10, 24, 28, 48, 40, 24, 36, 24, 18, 60, 16, 30, 48, 24, 100, 84, 72, 48, 14, 120, 30, 48, 40, 36, 80, 24, 76, 18, 56, 60, 40, 48, 88, 30, 120, 48, 32, 24, 112, 300, 72, 84, 108, 72, 20, 48, 72, 42, 58, 120, 60, 30, 48, 96, 140, 120, 136, 36, 48, 240, 70, 24, 148, 228, 200, 18, 80, 168, 78, 120, 216, 120, 168, 48, 180, 264, 56, 60, 44, 120, 112, 48, 120, 96, 180, 48, 196, 336, 120, 300, 50, 72, 208, 84, 80, 108, 72, 72, 108, 60, 152, 48, 76, 72, 240, 42, 168, 174, 144, 120, 110, 60, 40, 30, 500, 48, 256, 192, 88, 420, 130, 120, 144, 408, 360, 36, 276, 48, 46, 240, 32, 210, 140, 24, 140, 444, 112, 228, 148, 600, 50, 36, 72, 240, 60, 168, 316, 78, 216, 240, 48, 216, 328, 120, 40, 168, 336, 48, 364, 180, 72, 264, 348, 168, 400, 120, 232, 132, 178, 120] |
[3, 8, 6, 20, 24, 16, 12, 24, 60, 10, 24, 28, 48, 40, 24, 36, 24, 18, 60, 16, 30, 48, 24, 100, 84, 72, 48, 14, 120, 30, 48, 40, 36, 80, 24, 76, 18, 56, 60, 40, 48, 88, 30, 120, 48, 32, 24, 112, 300, 72, 84, 108, 72, 20, 48, 72, 42, 58, 120, 60, 30, 48, 96, 140, 120, 136, 36, 48, 240, 70, 24, 148, 228, 200, 18, 80, 168, 78, 120, 216, 120, 168, 48, 180, 264, 56, 60, 44, 120, 112, 48, 120, 96, 180, 48, 196, 336, 120, 300, 50, 72, 208, 84, 80, 108, 72, 72, 108, 60, 152, 48, 76, 72, 240, 42, 168, 174, 144, 120, 110, 60, 40, 30, 500, 48, 256, 192, 88, 420, 130, 120, 144, 408, 360, 36, 276, 48, 46, 240, 32, 210, 140, 24, 140, 444, 112, 228, 148, 600, 50, 36, 72, 240, 60, 168, 316, 78, 216, 240, 48, 216, 328, 120, 40, 168, 336, 48, 364, 180, 72, 264, 348, 168, 400, 120, 232, 132, 178, 120] |
||
4.208398 seconds (4.23 k allocations: 249.266 KiB) |
4.208398 seconds (4.23 k allocations: 249.266 KiB) |
||
</pre> |
|||
=={{header|Sidef}}== |
|||
<lang ruby>func pisano_period_pp(p,k) is cached { |
|||
assert(k.is_pos, "k = #{k} must be positive") |
|||
assert(p.is_prime, "p = #{p} must be prime") |
|||
var (a, b, n) = (0, 1, p**k) |
|||
for k in (1..Inf) { |
|||
(a, b) = (b, (a+b) % n) |
|||
if ([a,b] == [0,1]) { |
|||
return k |
|||
} |
|||
} |
|||
} |
|||
func pisano_period(n) { |
|||
n.factor_map {|p,k| pisano_period_pp(p, k) }.lcm |
|||
} |
|||
say "Pisano periods for squares of primes p <= 15:" |
|||
say 15.primes.map {|p| pisano_period_pp(p, 2) } |
|||
say "\nPisano periods for primes p <= 180:" |
|||
say 180.primes.map {|p| pisano_period_pp(p, 1) } |
|||
say "\nPisano periods for integers n from 2 to 180:" |
|||
say pisano_period.map(2..180)</lang> |
|||
{{out}} |
|||
<pre> |
|||
Pisano periods for squares of primes p <= 15: |
|||
[6, 24, 100, 112, 110, 364] |
|||
Pisano periods for primes p <= 180: |
|||
[3, 8, 20, 16, 10, 28, 36, 18, 48, 14, 30, 76, 40, 88, 32, 108, 58, 60, 136, 70, 148, 78, 168, 44, 196, 50, 208, 72, 108, 76, 256, 130, 276, 46, 148, 50, 316, 328, 336, 348, 178] |
|||
Pisano periods for integers n from 2 to 180: |
|||
[3, 8, 6, 20, 24, 16, 12, 24, 60, 10, 24, 28, 48, 40, 24, 36, 24, 18, 60, 16, 30, 48, 24, 100, 84, 72, 48, 14, 120, 30, 48, 40, 36, 80, 24, 76, 18, 56, 60, 40, 48, 88, 30, 120, 48, 32, 24, 112, 300, 72, 84, 108, 72, 20, 48, 72, 42, 58, 120, 60, 30, 48, 96, 140, 120, 136, 36, 48, 240, 70, 24, 148, 228, 200, 18, 80, 168, 78, 120, 216, 120, 168, 48, 180, 264, 56, 60, 44, 120, 112, 48, 120, 96, 180, 48, 196, 336, 120, 300, 50, 72, 208, 84, 80, 108, 72, 72, 108, 60, 152, 48, 76, 72, 240, 42, 168, 174, 144, 120, 110, 60, 40, 30, 500, 48, 256, 192, 88, 420, 130, 120, 144, 408, 360, 36, 276, 48, 46, 240, 32, 210, 140, 24, 140, 444, 112, 228, 148, 600, 50, 36, 72, 240, 60, 168, 316, 78, 216, 240, 48, 216, 328, 120, 40, 168, 336, 48, 364, 180, 72, 264, 348, 168, 400, 120, 232, 132, 178, 120] |
|||
</pre> |
|||
By assuming that '''Wall-Sun-Sun primes''' do not exist, we can compute the Pisano period more efficiently, as illustrated below on Fermat numbers '''F_n = 2^(2^n) + 1''': |
|||
<lang ruby>func pisano_period_pp(p, k=1) { |
|||
(p - kronecker(5, p)).divisors.first_by {|d| fibmod(d, p) == 0 } * p**(k-1) |
|||
} |
|||
func pisano_period(n) { |
|||
return 0 if (n <= 0) |
|||
return 1 if (n == 1) |
|||
var d = n.factor_map {|p,k| pisano_period_pp(p, k) }.lcm |
|||
3.times {|k| |
|||
var t = d<<k |
|||
if ((fibmod(t, n) == 0) && (fibmod(t+1, n) == 1)) { |
|||
return t |
|||
} |
|||
} |
|||
} |
|||
for k in (1..8) { |
|||
say ("Pisano(F_#{k}) = ", pisano_period(2**(2**k) + 1)) |
|||
}</lang> |
|||
{{out}} |
|||
<pre> |
|||
Pisano(F_1) = 20 |
|||
Pisano(F_2) = 36 |
|||
Pisano(F_3) = 516 |
|||
Pisano(F_4) = 14564 |
|||
Pisano(F_5) = 2144133760 |
|||
Pisano(F_6) = 4611702838532647040 |
|||
Pisano(F_7) = 28356863910078205764000346543980814080 |
|||
Pisano(F_8) = 3859736307910542962840356678888855900560939475751238269689837480239178278912 |
|||
</pre> |
</pre> |
Revision as of 13:27, 3 March 2020
The Fibonacci sequence taken modulo 2 is a periodic sequence of period 3 : 0, 1, 1, 0, 1, 1, ...
For any integer n, the Fibonacci sequence taken modulo n is periodic and the period is called the Pisano period.
Let call pisano, the Pisano period (pisano(2) = 3). If m and n are coprime, pisano(m*n) = lcm(pisano(m),pisano(n)).
Therefore calculate the Pisano period of an integer m is accomplished by calculating the Pisano periods of the prime powers in the prime decomposition of m.
- Task
Write 2 functions: pisanoPrime(p,k) and pisano(m).
pisanoPrime(p,k) return the Pisano period of pk where p is prime and k is a positive integer.
pisano(m) uses pisanoPrime to return the Pisano period of m where m is a positive integer.
Print pisanoPrime(p,2) for every prime lower than 15.
Print pisanoPrime(p,1) for every prime lower than 180.
Print pisano(m) for every integer from 2 to 180.
- Related tasks
Go
<lang go>package main
import "fmt"
func gcd(a, b uint) uint {
if b == 0 { return a } return gcd(b, a%b)
}
func lcm(a, b uint) uint {
return a / gcd(a, b) * b
}
func ipow(x, p uint) uint {
prod := uint(1) for p > 0 { if p&1 != 0 { prod *= x } p >>= 1 x *= x } return prod
}
// Gets the prime decomposition of n. func getPrimes(n uint) []uint {
var primes []uint for i := uint(2); i <= n; i++ { div := n / i mod := n % i for mod == 0 { primes = append(primes, i) n = div div = n / i mod = n % i } } return primes
}
// OK for 'small' numbers. func isPrime(n uint) bool {
switch { case n < 2: return false case n%2 == 0: return n == 2 case n%3 == 0: return n == 3 default: d := uint(5) for d*d <= n { if n%d == 0 { return false } d += 2 if n%d == 0 { return false } d += 4 } return true }
}
// Calculates the Pisano period of 'm' from first principles. func pisanoPeriod(m uint) uint {
var p, c uint = 0, 1 for i := uint(0); i < m*m; i++ { p, c = c, (p+c)%m if p == 0 && c == 1 { return i + 1 } } return 1
}
// Calculates the Pisano period of p^k where 'p' is prime and 'k' is a positive intege. func pisanoPrime(p uint, k uint) uint {
if !isPrime(p) || k == 0 { return 0 // can't do this one } return pisanoPeriod(ipow(p, k))
}
// Calculates the Pisano period of 'm' using pisanoPrime. func pisano(m uint) uint {
primes := getPrimes(m) primePowers := make(map[uint]uint) for _, p := range primes { primePowers[p]++ } var pps []uint for k, v := range primePowers { pps = append(pps, pisanoPrime(k, v)) } if len(pps) == 1 { return pps[0] } f := pps[0] for i := 1; i < len(pps); i++ { f = lcm(f, pps[i]) } return f
}
func main() {
for p := uint(2); p < 15; p++ { pp := pisanoPrime(p, 2) if pp > 0 { fmt.Printf("pisanoPrime(%2d: 2) = %d\n", p, pp) } } fmt.Println() for p := uint(2); p < 180; p++ { pp := pisanoPrime(p, 1) if pp > 0 { fmt.Printf("pisanoPrime(%3d: 1) = %d\n", p, pp) } } fmt.Println() fmt.Println("pisano(n) for integers 'n' from 2 to 180 are:") for n := uint(2); n <= 180; n++ { fmt.Printf("%3d ", pisano(n)) if n != 2 && (n-1)%15 == 0 { fmt.Println() } } fmt.Println()
}</lang>
- Output:
pisanoPrime( 2: 2) = 6 pisanoPrime( 3: 2) = 24 pisanoPrime( 5: 2) = 100 pisanoPrime( 7: 2) = 112 pisanoPrime(11: 2) = 110 pisanoPrime(13: 2) = 364 pisanoPrime( 2: 1) = 3 pisanoPrime( 3: 1) = 8 pisanoPrime( 5: 1) = 20 pisanoPrime( 7: 1) = 16 pisanoPrime( 11: 1) = 10 pisanoPrime( 13: 1) = 28 pisanoPrime( 17: 1) = 36 pisanoPrime( 19: 1) = 18 pisanoPrime( 23: 1) = 48 pisanoPrime( 29: 1) = 14 pisanoPrime( 31: 1) = 30 pisanoPrime( 37: 1) = 76 pisanoPrime( 41: 1) = 40 pisanoPrime( 43: 1) = 88 pisanoPrime( 47: 1) = 32 pisanoPrime( 53: 1) = 108 pisanoPrime( 59: 1) = 58 pisanoPrime( 61: 1) = 60 pisanoPrime( 67: 1) = 136 pisanoPrime( 71: 1) = 70 pisanoPrime( 73: 1) = 148 pisanoPrime( 79: 1) = 78 pisanoPrime( 83: 1) = 168 pisanoPrime( 89: 1) = 44 pisanoPrime( 97: 1) = 196 pisanoPrime(101: 1) = 50 pisanoPrime(103: 1) = 208 pisanoPrime(107: 1) = 72 pisanoPrime(109: 1) = 108 pisanoPrime(113: 1) = 76 pisanoPrime(127: 1) = 256 pisanoPrime(131: 1) = 130 pisanoPrime(137: 1) = 276 pisanoPrime(139: 1) = 46 pisanoPrime(149: 1) = 148 pisanoPrime(151: 1) = 50 pisanoPrime(157: 1) = 316 pisanoPrime(163: 1) = 328 pisanoPrime(167: 1) = 336 pisanoPrime(173: 1) = 348 pisanoPrime(179: 1) = 178 pisano(n) for integers 'n' from 2 to 180 are: 3 8 6 20 24 16 12 24 60 10 24 28 48 40 24 36 24 18 60 16 30 48 24 100 84 72 48 14 120 30 48 40 36 80 24 76 18 56 60 40 48 88 30 120 48 32 24 112 300 72 84 108 72 20 48 72 42 58 120 60 30 48 96 140 120 136 36 48 240 70 24 148 228 200 18 80 168 78 120 216 120 168 48 180 264 56 60 44 120 112 48 120 96 180 48 196 336 120 300 50 72 208 84 80 108 72 72 108 60 152 48 76 72 240 42 168 174 144 120 110 60 40 30 500 48 256 192 88 420 130 120 144 408 360 36 276 48 46 240 32 210 140 24 140 444 112 228 148 600 50 36 72 240 60 168 316 78 216 240 48 216 328 120 40 168 336 48 364 180 72 264 348 168 400 120 232 132 178 120
Julia
The factorization method gives about a 3X speedup over just looking for the 0, 1 that marks the repeating sequence. <lang julia>using Primes
function pisano(p)
p < 2 && return 1 lastn, n = 0, 1 for i in 0:p^2 lastn, n = n, (lastn + n) % p if lastn == 0 && n == 1 return i + 1 end end return 1
end
pisanoprime(p, k) = (@assert(isprime(p)); pisano(p^k)) primedecomp(n) = [p for p in collect(factor(n))] pisanotask(n) = mapreduce(p -> pisanoprime(p[1], p[2]), lcm, primedecomp(n), init=1)
for i in 1:15
if isprime(i) println("pisanoPrime($i, 2) = ", pisanoprime(i, 2)) end
end
for i in 1:180
if isprime(i) println("pisanoPrime($i, 1) = ", pisanoprime(i, 1)) end
end
println("\nPisano(n) for n from 2 to 180:\n", [pisano(i) for i in 2:180])
@time [pisano(i) for i in 9999500:10000000]
println("\nPisano(n) using pisanoPrime for n from 2 to 180:\n", [pisanotask(i) for i in 2:180])
@time [pisanotask(i) for i in 9999500:10000000]</lang>
- Output:
pisanoPrime(2, 2) = 6 pisanoPrime(3, 2) = 24 pisanoPrime(5, 2) = 100 pisanoPrime(7, 2) = 112 pisanoPrime(11, 2) = 110 pisanoPrime(13, 2) = 364 pisanoPrime(2, 1) = 3 pisanoPrime(3, 1) = 8 pisanoPrime(5, 1) = 20 pisanoPrime(7, 1) = 16 pisanoPrime(11, 1) = 10 pisanoPrime(13, 1) = 28 pisanoPrime(17, 1) = 36 pisanoPrime(19, 1) = 18 pisanoPrime(23, 1) = 48 pisanoPrime(29, 1) = 14 pisanoPrime(31, 1) = 30 pisanoPrime(37, 1) = 76 pisanoPrime(41, 1) = 40 pisanoPrime(43, 1) = 88 pisanoPrime(47, 1) = 32 pisanoPrime(53, 1) = 108 pisanoPrime(59, 1) = 58 pisanoPrime(61, 1) = 60 pisanoPrime(67, 1) = 136 pisanoPrime(71, 1) = 70 pisanoPrime(73, 1) = 148 pisanoPrime(79, 1) = 78 pisanoPrime(83, 1) = 168 pisanoPrime(89, 1) = 44 pisanoPrime(97, 1) = 196 pisanoPrime(101, 1) = 50 pisanoPrime(103, 1) = 208 pisanoPrime(107, 1) = 72 pisanoPrime(109, 1) = 108 pisanoPrime(113, 1) = 76 pisanoPrime(127, 1) = 256 pisanoPrime(131, 1) = 130 pisanoPrime(137, 1) = 276 pisanoPrime(139, 1) = 46 pisanoPrime(149, 1) = 148 pisanoPrime(151, 1) = 50 pisanoPrime(157, 1) = 316 pisanoPrime(163, 1) = 328 pisanoPrime(167, 1) = 336 pisanoPrime(173, 1) = 348 pisanoPrime(179, 1) = 178 Pisano(n) for n from 2 to 180: [3, 8, 6, 20, 24, 16, 12, 24, 60, 10, 24, 28, 48, 40, 24, 36, 24, 18, 60, 16, 30, 48, 24, 100, 84, 72, 48, 14, 120, 30, 48, 40, 36, 80, 24, 76, 18, 56, 60, 40, 48, 88, 30, 120, 48, 32, 24, 112, 300, 72, 84, 108, 72, 20, 48, 72, 42, 58, 120, 60, 30, 48, 96, 140, 120, 136, 36, 48, 240, 70, 24, 148, 228, 200, 18, 80, 168, 78, 120, 216, 120, 168, 48, 180, 264, 56, 60, 44, 120, 112, 48, 120, 96, 180, 48, 196, 336, 120, 300, 50, 72, 208, 84, 80, 108, 72, 72, 108, 60, 152, 48, 76, 72, 240, 42, 168, 174, 144, 120, 110, 60, 40, 30, 500, 48, 256, 192, 88, 420, 130, 120, 144, 408, 360, 36, 276, 48, 46, 240, 32, 210, 140, 24, 140, 444, 112, 228, 148, 600, 50, 36, 72, 240, 60, 168, 316, 78, 216, 240, 48, 216, 328, 120, 40, 168, 336, 48, 364, 180, 72, 264, 348, 168, 400, 120, 232, 132, 178, 120] 13.821726 seconds (7 allocations: 4.281 KiB) Pisano(n) using pisanoPrime for n from 2 to 180: [3, 8, 6, 20, 24, 16, 12, 24, 60, 10, 24, 28, 48, 40, 24, 36, 24, 18, 60, 16, 30, 48, 24, 100, 84, 72, 48, 14, 120, 30, 48, 40, 36, 80, 24, 76, 18, 56, 60, 40, 48, 88, 30, 120, 48, 32, 24, 112, 300, 72, 84, 108, 72, 20, 48, 72, 42, 58, 120, 60, 30, 48, 96, 140, 120, 136, 36, 48, 240, 70, 24, 148, 228, 200, 18, 80, 168, 78, 120, 216, 120, 168, 48, 180, 264, 56, 60, 44, 120, 112, 48, 120, 96, 180, 48, 196, 336, 120, 300, 50, 72, 208, 84, 80, 108, 72, 72, 108, 60, 152, 48, 76, 72, 240, 42, 168, 174, 144, 120, 110, 60, 40, 30, 500, 48, 256, 192, 88, 420, 130, 120, 144, 408, 360, 36, 276, 48, 46, 240, 32, 210, 140, 24, 140, 444, 112, 228, 148, 600, 50, 36, 72, 240, 60, 168, 316, 78, 216, 240, 48, 216, 328, 120, 40, 168, 336, 48, 364, 180, 72, 264, 348, 168, 400, 120, 232, 132, 178, 120] 4.208398 seconds (4.23 k allocations: 249.266 KiB)
Sidef
<lang ruby>func pisano_period_pp(p,k) is cached {
assert(k.is_pos, "k = #{k} must be positive") assert(p.is_prime, "p = #{p} must be prime")
var (a, b, n) = (0, 1, p**k)
for k in (1..Inf) { (a, b) = (b, (a+b) % n)
if ([a,b] == [0,1]) { return k } }
}
func pisano_period(n) {
n.factor_map {|p,k| pisano_period_pp(p, k) }.lcm
}
say "Pisano periods for squares of primes p <= 15:" say 15.primes.map {|p| pisano_period_pp(p, 2) }
say "\nPisano periods for primes p <= 180:" say 180.primes.map {|p| pisano_period_pp(p, 1) }
say "\nPisano periods for integers n from 2 to 180:" say pisano_period.map(2..180)</lang>
- Output:
Pisano periods for squares of primes p <= 15: [6, 24, 100, 112, 110, 364] Pisano periods for primes p <= 180: [3, 8, 20, 16, 10, 28, 36, 18, 48, 14, 30, 76, 40, 88, 32, 108, 58, 60, 136, 70, 148, 78, 168, 44, 196, 50, 208, 72, 108, 76, 256, 130, 276, 46, 148, 50, 316, 328, 336, 348, 178] Pisano periods for integers n from 2 to 180: [3, 8, 6, 20, 24, 16, 12, 24, 60, 10, 24, 28, 48, 40, 24, 36, 24, 18, 60, 16, 30, 48, 24, 100, 84, 72, 48, 14, 120, 30, 48, 40, 36, 80, 24, 76, 18, 56, 60, 40, 48, 88, 30, 120, 48, 32, 24, 112, 300, 72, 84, 108, 72, 20, 48, 72, 42, 58, 120, 60, 30, 48, 96, 140, 120, 136, 36, 48, 240, 70, 24, 148, 228, 200, 18, 80, 168, 78, 120, 216, 120, 168, 48, 180, 264, 56, 60, 44, 120, 112, 48, 120, 96, 180, 48, 196, 336, 120, 300, 50, 72, 208, 84, 80, 108, 72, 72, 108, 60, 152, 48, 76, 72, 240, 42, 168, 174, 144, 120, 110, 60, 40, 30, 500, 48, 256, 192, 88, 420, 130, 120, 144, 408, 360, 36, 276, 48, 46, 240, 32, 210, 140, 24, 140, 444, 112, 228, 148, 600, 50, 36, 72, 240, 60, 168, 316, 78, 216, 240, 48, 216, 328, 120, 40, 168, 336, 48, 364, 180, 72, 264, 348, 168, 400, 120, 232, 132, 178, 120]
By assuming that Wall-Sun-Sun primes do not exist, we can compute the Pisano period more efficiently, as illustrated below on Fermat numbers F_n = 2^(2^n) + 1: <lang ruby>func pisano_period_pp(p, k=1) {
(p - kronecker(5, p)).divisors.first_by {|d| fibmod(d, p) == 0 } * p**(k-1)
}
func pisano_period(n) {
return 0 if (n <= 0) return 1 if (n == 1)
var d = n.factor_map {|p,k| pisano_period_pp(p, k) }.lcm
3.times {|k| var t = d<<k if ((fibmod(t, n) == 0) && (fibmod(t+1, n) == 1)) { return t } }
}
for k in (1..8) {
say ("Pisano(F_#{k}) = ", pisano_period(2**(2**k) + 1))
}</lang>
- Output:
Pisano(F_1) = 20 Pisano(F_2) = 36 Pisano(F_3) = 516 Pisano(F_4) = 14564 Pisano(F_5) = 2144133760 Pisano(F_6) = 4611702838532647040 Pisano(F_7) = 28356863910078205764000346543980814080 Pisano(F_8) = 3859736307910542962840356678888855900560939475751238269689837480239178278912