Pisano period

From Rosetta Code
Revision as of 01:30, 5 March 2020 by Blek (talk | contribs)
Pisano period is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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 length of the periodic cycle is referred to as the Pisano period.

Prime numbers are straightforward; the Pisano period of a prime number p is simply: pisano(p). The Pisano period of a composite number c may be found in different ways. It may be calculated directly: pisano(c), which works, but may be time consuming to find, especially for larger integers, or, it may be calculated by finding the least common multiple of the Pisano periods of each composite component.

E.G. Given a Pisano period function: pisano(x), and a least common multiple function lcm(x, y):

   pisano(m × n) is equivalent to lcm(pisano(m), pisano(n)) where m and n are coprime

A formulae to calculate the pisano period for integer powers k of prime numbers p is:

   pisano(pk) == p(k-1)pisano(p)

The equation is conjectured, no exceptions have been seen.

If a positive integer i is split into its prime factors then the second and first equations above can be applied to generate the pisano period.

Task

Write 2 functions: pisanoPrime(p,k) and pisano(m).

pisanoPrime(p,k) should return the Pisano period of pk where p is prime and k is a positive integer.

pisano(m) should use pisanoPrime to return the Pisano period of m where m is a prime 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 1 to 180.

Related tasks


Factor

Works with: Factor version 0.99 2020-01-23

<lang factor>USING: formatting fry grouping io kernel math math.functions math.primes math.primes.factors math.ranges sequences ;

pisano-period ( m -- n )
   [ 0 1 ] dip [ sq <iota> ] [ ] bi
   '[ drop tuck + _ mod 2dup [ zero? ] [ 1 = ] bi* and ]
   find 3nip [ 1 + ] [ 1 ] if* ;
pisano-prime ( p k -- n )
   over prime? [ "p must be prime." throw ] unless
   ^ pisano-period ;
pisano ( m -- n )
   group-factors [ first2 pisano-prime ] [ lcm ] map-reduce ;
show-pisano ( upto m -- )
   [ primes-upto ] dip
   [ 2dup pisano-prime "%d %d pisano-prime = %d\n" printf ]
   curry each nl ;

15 2 show-pisano 180 1 show-pisano

"n pisano for integers 'n' from 2 to 180:" print 2 180 [a,b] [ pisano ] map 15 group [ [ "%3d " printf ] each nl ] each</lang>

Output:
2 2 pisano-prime = 6
3 2 pisano-prime = 24
5 2 pisano-prime = 100
7 2 pisano-prime = 112
11 2 pisano-prime = 110
13 2 pisano-prime = 364

2 1 pisano-prime = 3
3 1 pisano-prime = 8
5 1 pisano-prime = 20
7 1 pisano-prime = 16
11 1 pisano-prime = 10
13 1 pisano-prime = 28
17 1 pisano-prime = 36
19 1 pisano-prime = 18
23 1 pisano-prime = 48
29 1 pisano-prime = 14
31 1 pisano-prime = 30
37 1 pisano-prime = 76
41 1 pisano-prime = 40
43 1 pisano-prime = 88
47 1 pisano-prime = 32
53 1 pisano-prime = 108
59 1 pisano-prime = 58
61 1 pisano-prime = 60
67 1 pisano-prime = 136
71 1 pisano-prime = 70
73 1 pisano-prime = 148
79 1 pisano-prime = 78
83 1 pisano-prime = 168
89 1 pisano-prime = 44
97 1 pisano-prime = 196
101 1 pisano-prime = 50
103 1 pisano-prime = 208
107 1 pisano-prime = 72
109 1 pisano-prime = 108
113 1 pisano-prime = 76
127 1 pisano-prime = 256
131 1 pisano-prime = 130
137 1 pisano-prime = 276
139 1 pisano-prime = 46
149 1 pisano-prime = 148
151 1 pisano-prime = 50
157 1 pisano-prime = 316
163 1 pisano-prime = 328
167 1 pisano-prime = 336
173 1 pisano-prime = 348
179 1 pisano-prime = 178

n pisano 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 

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 integer. 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) == 0 {
       return 1
   }
   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 1 to 180 are:")
   for n := uint(1); n <= 180; n++ {
       fmt.Printf("%3d ", pisano(n))
       if n != 1 && n%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 1 to 180 are:
  1   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 

Haskell

<lang Haskell>import qualified Data.Text as T

main = do

  putStrLn $ "PisanoPrime(p,2) for prime p lower than 15"
  print.map (flip pisanoPrime 2).filter isPrime $ [1..15]
  putStrLn $ "PisanoPrime(p,1) for prime p lower than 180"
  print.map (flip pisanoPrime 1).filter isPrime $ [1..180]
  let ns = [1..180] :: [Int]
  let xs = map pisanoPeriod ns
  let ys = map pisano ns
  let zs = map pisanoConjecture ns
  putStrLn $ "Pisano(m) for m from 1 to 180"
  putStrLn.see 15 $ map pisano [1..180]
  putStrLn $ "map pisanoPeriod [1..180] == map pisano [1..180] = " ++ (show $ xs == ys)
  putStrLn $ "map pisanoPeriod [1..180] == map pisanoConjecture [1..180] = " ++ (show $ ys == zs)

bagOf :: Int -> [a] -> a bagOf _ [] = [] bagOf n xs = let (us,vs) = splitAt n xs in us : bagOf n vs

see :: Show a => Int -> [a] -> String see n = unlines.map unwords.bagOf n.map (T.unpack.T.justifyLeft 3 ' '.T.pack.show)

fibMod :: Integral a => a -> [a] fibMod 1 = repeat 0 fibMod n = fib

   where fib = 0 : 1 : zipWith (\x y -> rem (x+y) n) fib (tail fib)

pisanoPeriod :: Integral a => a -> a pisanoPeriod m | m <= 0 = 0 pisanoPeriod 1 = 1 pisanoPeriod m = go 1 (tail $ fibMod m)

   where
   go t (0:1:_) = t
   go t (_:xs)  = go (succ t) xs

primes :: Integral a => [a] primes = 2:3:5:7:[p | p <- [11,13..], isPrime p]

limitDivisor :: Integral a => a -> a limitDivisor = floor.(+0.05).sqrt.fromIntegral

isPrime :: Integral a => a -> Bool isPrime p | p < 2 = False isPrime p = go primes

   where
   stop = limitDivisor p
   go (n:_) | stop < n = True
   go (n:ns) = if 0 == rem p n then False else go ns
   go [] = True

factor :: Integral a => a -> [(a,a)] factor n | n <= 1 = [] factor n = if null ans then [(n,1)] else ans

   where
   ans = go n primes
   fun x d c = if 0 /= rem x d then (x,c) else fun (quot x d) d (succ c)
   go 1 _  = []
   go _ [] = []
   go x (d:ds) | 0 /= rem x d = go x $ dropWhile ((0 /=).rem x) ds
   go x (d:ds) = let (u,c) = fun (quot x d) d 1 in (d,c):go u ds

pisanoPrime :: Integral a => a -> a -> a pisanoPrime p k | p <= 0 || k < 0 = 0 pisanoPrime p k = pisanoPeriod $ p^k

pisano :: Integral a => a -> a pisano m | m < 1 = 0 pisano 1 = 1 pisano m = foldl1 lcm.map (uncurry pisanoPrime) $ factor m

pisanoConjecture :: Integral a => a -> a pisanoConjecture m | m < 1 = 0 pisanoConjecture 1 = 1 pisanoConjecture m = foldl1 lcm.map (uncurry pisanoPrime') $ factor m

   where pisanoPrime' p k = (p^(k-1))*(pisanoPeriod p)</lang>
Output:
PisanoPrime(p,2) for prime p lower than 15
[6,24,100,112,110,364]
PisanoPrime(p,1) for prime p lower than 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(m) for m from 1 to 180
1   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

map pisanoPeriod [1..180] == map pisano [1..180] = True
map pisanoPeriod [1..180] == map pisanoConjecture [1..180] = True

Julia

<lang julia>using Primes

const pisanos = Dict{Int, Int}() function pisano(p)

   p < 2 && return 1
   (i = get(pisanos, p, 0)) > 0 && return i
   lastn, n = 0, 1
   for i in 1:p^2
       lastn, n = n, (lastn + n) % p
       if lastn == 0 && n == 1
           pisanos[p] = i
           return i
       end
   end
   return 1

end

pisanoprime(p, k) = (@assert(isprime(p)); p^(k-1) * pisano(p)) pisanotask(n) = mapreduce(p -> pisanoprime(p[1], p[2]), lcm, collect(factor(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]) println("\nPisano(n) using pisanoPrime for n from 2 to 180:\n", [pisanotask(i) for i in 2:180])

</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]

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]

Perl 6

Works with: Rakudo version 2020.02

Didn't bother making two differently named routines, just made a multi that will auto dispatch to the correct candidate. <lang perl6>use Prime::Factor;

constant @fib := 1,1,*+*…*;

my %cache;

multi pisano-period (Int $p where *.is-prime, Int $k where * > 0 = 1 ) {

   return %cache{"$p|$k"} if %cache{"$p|$k"};
   my $fibmod = @fib.map: * % $p**$k;
   %cache{"$p|$k"} = (1..*).first: { !$fibmod[$_-1] and ($fibmod[$_] == 1) }

}

multi pisano-period (Int $p where * > 0, Int $k where * > 0 = 1 ) {

   [lcm] prime-factors($p).squish.map: { samewith $_, $k }

}


put "Pisano period (p, 2) for primes less than 50"; put (map { pisano-period($_, 2) }, ^50 .grep: *.is-prime )».fmt('%4d');

put "\nPisano period (p, 1) for primes less than 180"; .put for (map { pisano-period($_, 1) }, ^180 .grep: *.is-prime )».fmt('%4d').batch(15);

put "\nPisano period (p, 1) for integers 1 to 180"; .put for (1..180).map( { pisano-period($_) } )».fmt('%4d').batch(15);</lang>

Output:
Pisano period (p, 2) for primes less than 50
   6   24  100  112  110  364  612  342 1104  406  930 2812 1640 3784 1504

Pisano period (p, 1) for primes less than 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 period (p, 1) for integers 1 to 180
   1    3    8    3   20   24   16    3    8   60   10   24   28   48   40
   3   36   24   18   60   16   30   48   24   20   84    8   48   14  120
  30    3   40   36   80   24   76   18   56   60   40   48   88   30   40
  48   32   24   16   60   72   84  108   24   20   48   72   42   58  120
  60   30   16    3  140  120  136   36   48  240   70   24  148  228   40
  18   80  168   78   60    8  120  168   48  180  264   56   30   44  120
 112   48  120   96  180   24  196   48   40   60   50   72  208   84   80
 108   72   24  108   60  152   48   76   72  240   42   56  174  144  120
  10   60   40   30   20   48  256    3   88  420  130  120  144  408   40
  36  276   48   46  240   32  210  140   24  140  444   16  228  148  120
  50   18   72  240   60  168  316   78  216   60   48   24  328  120   40
 168  336   48   28  180   72  264  348  168   80   30  232  132  178  120

Phix

<lang Phix>function pisano_period(integer m) -- Calculates the Pisano period of 'm' from first principles. (copied from Go)

   integer p = 0, c = 1
   for i=0 to m*m-1 do
       {p, c} = {c, mod(p+c,m)}
       if p == 0 and c == 1 then
           return i + 1
       end if
   end for
   return 1

end function

function pisanoPrime(integer p, k)

   if not is_prime(p) or k=0 then ?9/0 end if
   return power(p,k-1)*pisano_period(p)

end function

function pisano(integer m) -- Calculates the Pisano period of 'm' using pisanoPrime.

   sequence s = prime_factors(m, true, get_maxprime(m))&0,
            pps = {}
   integer k = 1, p = s[1]
   for i=2 to length(s) do
       integer n = s[i]
       if n!=p then
           pps = append(pps,pisanoPrime(p,k))
           {k,p} = {1,n}
       else
           k += 1
       end if
   end for
   return lcm(pps)

end function

procedure p(integer k, lim) -- test harness

   printf(1,"pisanoPrimes")
   integer pdx = 1, c = 0
   while true do
       integer p = get_prime(pdx)
       if p>=lim then exit end if
       c += 1
       if c=7 then puts(1,"\n            ") c = 1
       elsif pdx>1 then puts(1,", ") end if
       printf(1,"(%3d,%d)=%3d",{p,k,pisanoPrime(p,k)})
       pdx += 1
   end while
   printf(1,"\n")

end procedure p(2,15) p(1,180)

sequence p180 = {} for n=2 to 180 do p180 &= pisano(n) end for printf(1,"pisano(2..180):\n") pp(p180,{pp_IntFmt,"%4d",pp_IntCh,false})

printf(1,"\nTiming comparison:\n") for i=1 to 2 do

   atom t0 = time()
   string fn = {"pisano_period","pisano"}[i],
          via = {"",", via pisanoPrime()"}[i]  
   integer f = routine_id(fn)
   for n=2 to 10000 do
       {} = f(n)
   end for
   t0 = time()-t0
   printf(1,"%s(2..10000)%s:%s\n",{fn,via,elapsed(t0)})

end for</lang>

Output:
pisanoPrimes(  2,2)=  6, (  3,2)= 24, (  5,2)=100, (  7,2)=112, ( 11,2)=110, ( 13,2)=364
pisanoPrimes(  2,1)=  3, (  3,1)=  8, (  5,1)= 20, (  7,1)= 16, ( 11,1)= 10, ( 13,1)= 28
            ( 17,1)= 36, ( 19,1)= 18, ( 23,1)= 48, ( 29,1)= 14, ( 31,1)= 30, ( 37,1)= 76
            ( 41,1)= 40, ( 43,1)= 88, ( 47,1)= 32, ( 53,1)=108, ( 59,1)= 58, ( 61,1)= 60
            ( 67,1)=136, ( 71,1)= 70, ( 73,1)=148, ( 79,1)= 78, ( 83,1)=168, ( 89,1)= 44
            ( 97,1)=196, (101,1)= 50, (103,1)=208, (107,1)= 72, (109,1)=108, (113,1)= 76
            (127,1)=256, (131,1)=130, (137,1)=276, (139,1)= 46, (149,1)=148, (151,1)= 50
            (157,1)=316, (163,1)=328, (167,1)=336, (173,1)=348, (179,1)=178
pisano(2..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}

Timing comparison:
pisano_period(2..10000):8.9s
pisano(2..10000), via pisanoPrime():3.4s

Python

Uses the SymPy library.

<lang python>from sympy import isprime, lcm, factorint, primerange from functools import reduce


def pisano1(m):

   "Simple definition"
   if m < 2:
       return 1
   lastn, n = 0, 1
   for i in range(m ** 2):
       lastn, n = n, (lastn + n) % m
       if lastn == 0 and n == 1:
           return i + 1
   return 1

def pisanoprime(p, k):

   "Use conjecture π(p ** k) == p ** (k − 1) * π(p) for prime p and int k > 1"
   assert isprime(p) and k > 0
   return p ** (k - 1) * pisano1(p)

def pisano_mult(m, n):

   "pisano(m*n) where m and n assumed coprime integers"
   return lcm(pisano1(m), pisano1(n))

def pisano2(m):

   "Uses prime factorization of m"
   return reduce(lcm, (pisanoprime(prime, mult)
                       for prime, mult in factorint(m).items()), 1)


if __name__ == '__main__':

   for n in range(1, 181):
       assert pisano1(n) == pisano2(n), "Wall-Sun-Sun prime exists??!!"
   print("\nPisano period (p, 2) for primes less than 50\n ",
         [pisanoprime(prime, 2) for prime in primerange(1, 50)])
   print("\nPisano period (p, 1) for primes less than 180\n ",
         [pisanoprime(prime, 1) for prime in primerange(1, 180)])
   print("\nPisano period (p) for integers 1 to 180")
   for i in range(1, 181):
       print(" %3d" % pisano2(i), end="" if i % 10 else "\n")</lang>
Output:
Pisano period (p, 2) for primes less than 50
  [6, 24, 100, 112, 110, 364, 612, 342, 1104, 406, 930, 2812, 1640, 3784, 1504]

Pisano period (p, 1) for primes less than 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 period (p) for integers 1 to 180
   1   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

REXX

<lang rexx>/*REXX pgm calculates pisano period for a range of N, and pisanoPrime(N,m) [for primes]*/ numeric digits 500 /*ensure enough decimal digits for Fib.*/ parse arg lim1 lim2 lim3 . /*obtain optional arguments from the CL*/ if lim1== | lim1=="," then lim1= 15 - 1 /*Not specified? Then use the default.*/ if lim2== | lim2=="," then lim2= 180 - 1 /* " " " " " " */ if lim3== | lim3=="," then lim3= 180 /* " " " " " " */ call fib /*generate some Fibonacci numbers. */

         do i=1  for max(lim1, lim2, lim3);   call pisano(i)      /*find some pisano #s*/
         end   /*i*/
   do p=1  for lim1;  if \isPrime(p)  then iterate  /*Not prime?  Then skip this number*/
   say '   pisanoPrime('right(p, length(lim1))", 2) = "right(pisanoPrime(p, 2), 5)
   end   /*pp*/

say

   do p=1  for lim2;  if \isPrime(p)  then iterate  /*Not prime?  Then skip this number*/
   say '   pisanoPrime('right(p, length(lim2))", 1) = "right(pisanoPrime(p, 1), 5)
   end   /*pp*/

say say center(' pisano numbers for 1──►'lim3" ", 20*4 - 1, "═") /*display a title. */ $=

   do j=1  for lim3;      $= $ right(@.j, 3)    /*append pisano number to the  $  list.*/
   if j//20==0  then do;  say substr($, 2);  $= /*only display twenty numbers to a line*/
                     end
   end
                          say substr($, 2)      /*possible display any residuals──►term*/

exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ isPrime: parse arg #; if #<2 then return 0; if wordpos(#,'2 3 5 7 11')\==0 then return 1

        if #//2==0  then return 0;  do k=3  by 2  while k*k<=#; if #//k==0  then return 0
                                    end  /*k*/
        return 1                                /*primality test for smallish primes.  */

/*──────────────────────────────────────────────────────────────────────────────────────*/ fib: procedure expose fib.; parse arg x; fib.=.; if x== then x= 1000

    fib.0= 0;  fib.1= 1;                        if fib.x\==.  then return fib.x
                 do j=2  for x-1;    _= j-1;    __= j-2;           fib.j=  fib.__ + fib._
                 end   /*j*/;    return fib.j

/*──────────────────────────────────────────────────────────────────────────────────────*/ pisano: procedure expose @. fib.; parse arg m; if m==1 then do; @.m=1; return 1; end

               do j=1;  _= j+1;                 if fib.j//m==0 & fib._//m==1  then leave
               end   /*j*/
       @.m= j;                   return j

/*──────────────────────────────────────────────────────────────────────────────────────*/ pisanoPrime: procedure expose @. fib.; parse arg m,n; return pisano(m**n)</lang>

output   when using the default inputs:

(Shown at   3/4   size.)

   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 numbers for 1──►180 ══════════════════════════
  1   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

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)
   1..Inf -> first_by {
       (a, b) = (b, (a+b) % n)
       (a == 0) && (b == 1)
   }

}

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

zkl

Library: GMP

GNU Multiple Precision Arithmetic Library for prime testing

<lang zkl>var [const] BI=Import("zklBigNum"); // libGMP

fcn pisanoPeriod(p){

  if(p<2) return(0);
  lastn,n,t := 0,1,0;
  foreach i in ([0..p*p]){
     t,n,lastn = n, (lastn + n) % p, t;
     if(lastn==0 and n==1) return(i + 1);
  }
  1

} fcn pisanoPrime(p,k){

  _assert_(BI(p).probablyPrime(), "%s is not a prime number".fmt(p));
  pisanoPeriod(p.pow(k))

}</lang> <lang zkl>println("Pisano period (p, 2) for primes less than 50:"); [1..50].pump(List,BI,"probablyPrime",Void.Filter, pisanoPrime.fp1(2))

  .concat(" ","   ").println();

println("Pisano period (p, 1) for primes less than 180:"); [1..180].pump(List,BI,"probablyPrime",Void.Filter, pisanoPrime.fp1(1))

 .pump(Void,T(Void.Read,14,False),fcn{ vm.arglist.apply("%4d".fmt).concat().println() });</lang>
Output:
Pisano period (p, 2) for primes less than 50:
   6 24 100 112 110 364 612 342 1104 406 930 2812 1640 3784 1504
Pisano period (p, 1) for primes less than 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

<lang zkl>fcn pisano(m){

  primeFactors(m).pump(Dictionary().incV)  //18 --> (2,3,3) --> ("2":1, "3":2)
    .reduce(fcn(z,[(k,v])){ lcm(z,pisanoPrime(k.toInt(),v)) },1)

}

fcn lcm(a,b){ a / a.gcd(b) * b } fcn primeFactors(n){ // Return a list of prime factors of n

  acc:=fcn(n,k,acc,maxD){  // k is 2,3,5,7,9,... not optimum
     if(n==1 or k>maxD) acc.close();
     else{

q,r:=n.divr(k); // divr-->(quotient,remainder) if(r==0) return(self.fcn(q,k,acc.write(k),q.toFloat().sqrt())); return(self.fcn(n,k+1+k.isOdd,acc,maxD)) # both are tail recursion

     }
  }(n,2,Sink(List),n.toFloat().sqrt());
  m:=acc.reduce('*,1);      // mulitply factors
  if(n!=m) acc.append(n/m); // opps, missed last factor
  else acc;

}</lang> <lang zkl>println("Pisano(m) for integers 1 to 180:"); [1..180].pump(List, pisano, "%4d".fmt)

 .pump(Void,T(Void.Read,14,False),fcn{ vm.arglist.concat().println() });</lang>
Output:
Pisano(m) for integers 1 to 180:
   1   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