Pisano period
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 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 1 to 180.
- Related tasks
Factor
<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 ipow(p, k-1) * pisanoPeriod(p)
}
// 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" putStrLn . see 15 . map (`pisanoPrime` 2) . filter isPrime $ [1 .. 15] putStrLn "PisanoPrime(p,1) for prime p lower than 180" putStrLn . see 15 . map (`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.justifyRight 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
powMod
:: Integral a => a -> a -> a -> a
powMod _ _ k
| k < 0 = error "negative power"
powMod m _ _
| 1 == abs m = 0
powMod m p k
| 1 == abs p = mod v m where v | 1 == p || even k = 1 | otherwise = p
powMod m p k = go p k
where to x y = mod (x * y) m go _ 0 = 1 go u 1 = mod u m go u i | even i = to w w | otherwise = to u (to w w) where w = go u (quot i 2)
-- Fermat primality test probablyPrime
:: Integral a => a -> Bool
probablyPrime p
| p < 2 || even p = 2 == p | otherwise = 1 == powMod p 2 (p - 1)
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
| not $ probablyPrime p = False
isPrime p = go primes
where stop = limitDivisor p go (n:_) | stop < n = True go (n:ns) = (0 /= rem p n) && go ns go [] = True
factor
:: Integral a => a -> [(a, a)]
factor n
| n <= 1 = []
factor n = go n primes
where fun x d c | 0 /= rem x d = (x, c) | otherwise = 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
Java
Use efficient algorithm to calculate period.
<lang Java> import java.util.ArrayList; import java.util.Collections; import java.util.HashMap; import java.util.List; import java.util.Map; import java.util.TreeMap;
public class PisanoPeriod {
public static void main(String[] args) { System.out.printf("Print pisano(p^2) for every prime p lower than 15%n"); for ( long i = 2 ; i < 15 ; i++ ) { if ( isPrime(i) ) { long n = i*i; System.out.printf("pisano(%d) = %d%n", n, pisano(n)); } }
System.out.printf("%nPrint pisano(p) for every prime p lower than 180%n"); for ( long n = 2 ; n < 180 ; n++ ) { if ( isPrime(n) ) { System.out.printf("pisano(%d) = %d%n", n, pisano(n)); } }
System.out.printf("%nPrint pisano(n) for every integer from 1 to 180%n"); for ( long n = 1 ; n <= 180 ; n++ ) { System.out.printf("%3d ", pisano(n)); if ( n % 10 == 0 ) { System.out.printf("%n"); } }
} private static final boolean isPrime(long test) { if ( test == 2 ) { return true; } if ( test % 2 == 0 ) { return false; } for ( long i = 3 ; i <= Math.sqrt(test) ; i += 2 ) { if ( test % i == 0 ) { return false; } } return true; }
private static Map<Long,Long> PERIOD_MEMO = new HashMap<>(); static { PERIOD_MEMO.put(2L, 3L); PERIOD_MEMO.put(3L, 8L); PERIOD_MEMO.put(5L, 20L); } // See http://webspace.ship.edu/msrenault/fibonacci/fib.htm private static long pisano(long n) { if ( PERIOD_MEMO.containsKey(n) ) { return PERIOD_MEMO.get(n); } if ( n == 1 ) { return 1; } Map<Long,Long> factors = getFactors(n); // Special cases // pisano(2^k) = 3*n/2 if ( factors.size() == 1 & factors.get(2L) != null && factors.get(2L) > 0 ) { long result = 3 * n / 2; PERIOD_MEMO.put(n, result); return result; } // pisano(5^k) = 4*n if ( factors.size() == 1 & factors.get(5L) != null && factors.get(5L) > 0 ) { long result = 4*n; PERIOD_MEMO.put(n, result); return result; } // pisano(2*5^k) = 6*n if ( factors.size() == 2 & factors.get(2L) != null && factors.get(2L) == 1 && factors.get(5L) != null && factors.get(5L) > 0 ) { long result = 6*n; PERIOD_MEMO.put(n, result); return result; } List<Long> primes = new ArrayList<>(factors.keySet()); long prime = primes.get(0); if ( factors.size() == 1 && factors.get(prime) == 1 ) { List<Long> divisors = new ArrayList<>(); if ( n % 10 == 1 || n % 10 == 9 ) { for ( long divisor : getDivisors(prime-1) ) { if ( divisor % 2 == 0 ) { divisors.add(divisor); } } } else { List<Long> pPlus1Divisors = getDivisors(prime+1); for ( long divisor : getDivisors(2*prime+2) ) { if ( ! pPlus1Divisors.contains(divisor) ) { divisors.add(divisor); } } } Collections.sort(divisors); for ( long divisor : divisors ) { if ( fibModIdentity(divisor, prime) ) { PERIOD_MEMO.put(prime, divisor); return divisor; } } throw new RuntimeException("ERROR 144: Divisor not found."); } long period = (long) Math.pow(prime, factors.get(prime)-1) * pisano(prime); for ( int i = 1 ; i < primes.size() ; i++ ) { prime = primes.get(i); period = lcm(period, (long) Math.pow(prime, factors.get(prime)-1) * pisano(prime)); } PERIOD_MEMO.put(n, period); return period; } // Use Matrix multiplication to compute Fibonacci numbers. private static boolean fibModIdentity(long n, long mod) { long aRes = 0; long bRes = 1; long cRes = 1; long aBase = 0; long bBase = 1; long cBase = 1; while ( n > 0 ) { if ( n % 2 == 1 ) { long temp1 = 0, temp2 = 0, temp3 = 0; if ( aRes > SQRT || aBase > SQRT || bRes > SQRT || bBase > SQRT || cBase > SQRT || cRes > SQRT ) { temp1 = (multiply(aRes, aBase, mod) + multiply(bRes, bBase, mod)) % mod; temp2 = (multiply(aBase, bRes, mod) + multiply(bBase, cRes, mod)) % mod; temp3 = (multiply(bBase, bRes, mod) + multiply(cBase, cRes, mod)) % mod; } else { temp1 = ((aRes*aBase % mod) + (bRes*bBase % mod)) % mod; temp2 = ((aBase*bRes % mod) + (bBase*cRes % mod)) % mod; temp3 = ((bBase*bRes % mod) + (cBase*cRes % mod)) % mod; } aRes = temp1; bRes = temp2; cRes = temp3; } n >>= 1L; long temp1 = 0, temp2 = 0, temp3 = 0; if ( aBase > SQRT || bBase > SQRT || cBase > SQRT ) { temp1 = (multiply(aBase, aBase, mod) + multiply(bBase, bBase, mod)) % mod; temp2 = (multiply(aBase, bBase, mod) + multiply(bBase, cBase, mod)) % mod; temp3 = (multiply(bBase, bBase, mod) + multiply(cBase, cBase, mod)) % mod; } else { temp1 = ((aBase*aBase % mod) + (bBase*bBase % mod)) % mod; temp2 = ((aBase*bBase % mod) + (bBase*cBase % mod)) % mod; temp3 = ((bBase*bBase % mod) + (cBase*cBase % mod)) % mod; } aBase = temp1; bBase = temp2; cBase = temp3; } return aRes % mod == 0 && bRes % mod == 1 && cRes % mod == 1; }
private static final long SQRT = (long) Math.sqrt(Long.MAX_VALUE);
// Result is a*b % mod, without overflow. public static final long multiply(long a, long b, long modulus) { //System.out.println(" multiply : a = " + a + ", b = " + b + ", mod = " + modulus); long x = 0; long y = a % modulus; long t; while ( b > 0 ) { if ( b % 2 == 1 ) { t = x + y; x = (t > modulus ? t-modulus : t); } t = y << 1; y = (t > modulus ? t-modulus : t); b >>= 1; } //System.out.println(" multiply : answer = " + (x % modulus)); return x % modulus; }
private static final List<Long> getDivisors(long number) { List<Long> divisors = new ArrayList<>(); long sqrt = (long) Math.sqrt(number); for ( long i = 1 ; i <= sqrt ; i++ ) { if ( number % i == 0 ) { divisors.add(i); long div = number / i; if ( div != i ) { divisors.add(div); } } } return divisors; }
public static long lcm(long a, long b) { return a*b/gcd(a,b); }
public static long gcd(long a, long b) { if ( b == 0 ) { return a; } return gcd(b, a%b); }
private static final Map<Long,Map<Long,Long>> allFactors = new TreeMap<Long,Map<Long,Long>>(); static { Map<Long,Long> factors = new TreeMap<Long,Long>(); factors.put(2L, 1L); allFactors.put(2L, factors); }
public static Long MAX_ALL_FACTORS = 100000L;
public static final Map<Long,Long> getFactors(Long number) { if ( allFactors.containsKey(number) ) { return allFactors.get(number); } Map<Long,Long> factors = new TreeMap<Long,Long>(); if ( number % 2 == 0 ) { Map<Long,Long> factorsdDivTwo = getFactors(number/2); factors.putAll(factorsdDivTwo); factors.merge(2L, 1L, (v1, v2) -> v1 + v2); if ( number < MAX_ALL_FACTORS ) { allFactors.put(number, factors); } return factors; } boolean prime = true; long sqrt = (long) Math.sqrt(number); for ( long i = 3 ; i <= sqrt ; i += 2 ) { if ( number % i == 0 ) { prime = false; factors.putAll(getFactors(number/i)); factors.merge(i, 1L, (v1, v2) -> v1 + v2); if ( number < MAX_ALL_FACTORS ) { allFactors.put(number, factors); } return factors; } } if ( prime ) { factors.put(number, 1L); if ( number < MAX_ALL_FACTORS ) { allFactors.put(number, factors); } } return factors; }
} </lang>
- Output:
Print pisano(p^2) for every prime p lower than 15 pisano(4) = 6 pisano(9) = 24 pisano(25) = 100 pisano(49) = 112 pisano(121) = 110 pisano(169) = 364 Print pisano(p) for every prime p lower than 180 pisano(2) = 3 pisano(3) = 8 pisano(5) = 20 pisano(7) = 16 pisano(11) = 10 pisano(13) = 28 pisano(17) = 36 pisano(19) = 18 pisano(23) = 48 pisano(29) = 14 pisano(31) = 30 pisano(37) = 76 pisano(41) = 40 pisano(43) = 88 pisano(47) = 32 pisano(53) = 108 pisano(59) = 58 pisano(61) = 60 pisano(67) = 136 pisano(71) = 70 pisano(73) = 148 pisano(79) = 78 pisano(83) = 168 pisano(89) = 44 pisano(97) = 196 pisano(101) = 50 pisano(103) = 208 pisano(107) = 72 pisano(109) = 108 pisano(113) = 76 pisano(127) = 256 pisano(131) = 130 pisano(137) = 276 pisano(139) = 46 pisano(149) = 148 pisano(151) = 50 pisano(157) = 316 pisano(163) = 328 pisano(167) = 336 pisano(173) = 348 pisano(179) = 178 Print pisano(n) for every integer 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
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
<lang perl>use strict; use warnings; use feature 'say'; use ntheory qw(primes factor_exp lcm);
sub pisano_period_pp {
my($a, $b, $n, $k) = (0, 1, $_[0]**$_[1]); while (++$k) { ($a, $b) = ($b, ($a+$b) % $n); return $k if $a == 0 and $b == 1; }
}
sub pisano_period {
(lcm map { pisano_period_pp($$_[0],$$_[1]) } factor_exp($_[0])) or 1;
}
sub display { (sprintf "@{['%5d' x @_]}", @_) =~ s/(.{75})/$1\n/gr }
say "Pisano periods for squares of primes p <= 50:\n", display( map { pisano_period_pp($_, 2) } @{primes(1, 50)} ),
"\nPisano periods for primes p <= 180:\n", display( map { pisano_period_pp($_, 1) } @{primes(1, 180)} ),
"\n\nPisano periods for integers n from 1 to 180:\n", display( map { pisano_period ($_ ) } 1..180 );</lang>
- Output:
Pisano periods for squares of primes p <= 50: 6 24 100 112 110 364 612 342 1104 406 930 2812 1640 3784 1504 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 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
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.
if m=1 then return 1 end if 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=1 to 180 do p180 &= pisano(n) end for printf(1,"pisano(1..180):\n") pp(p180,{pp_IntFmt,"%4d",pp_IntCh,false})</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(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}
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
Raku
(formerly Perl 6)
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 ) {
[lcm] prime-factors($p).Bag.map: { samewith .key, .value }
}
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 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 1 to 180:" say pisano_period.map(1..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 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]
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
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