1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
|
-- |
-- Module: Math.NumberTheory.Recurrences.Linear
-- Copyright: (c) 2011 Daniel Fischer
-- Licence: MIT
-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>
--
-- Efficient calculation of linear recurrent sequences, including Fibonacci and Lucas sequences.
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE PostfixOperators #-}
module Math.NumberTheory.Recurrences.Linear
( factorial
, factorialFactors
, fibonacci
, fibonacciPair
, lucas
, lucasPair
, generalLucas
) where
import Data.Bits
import Data.List.Infinite (Infinite(..), (...))
import qualified Data.List.Infinite as Inf
import Numeric.Natural
import Math.NumberTheory.Primes
-- | Infinite zero-based table of factorials.
--
-- >>> take 5 factorial
-- [1,1,2,6,24]
--
-- The time-and-space behaviour of 'factorial' is similar to described in
-- "Math.NumberTheory.Recurrences.Bilinear#memory".
factorial :: (Num a, Enum a) => Infinite a
factorial = Inf.scanl (*) 1 (1...)
{-# SPECIALIZE factorial :: Infinite Int #-}
{-# SPECIALIZE factorial :: Infinite Word #-}
{-# SPECIALIZE factorial :: Infinite Integer #-}
{-# SPECIALIZE factorial :: Infinite Natural #-}
-- | Prime factors of a factorial.
--
-- > factorialFactors n == factorise (factorial !! n)
--
-- >>> factorialFactors 10
-- [(Prime 2,8),(Prime 3,4),(Prime 5,2),(Prime 7,1)]
factorialFactors :: Word -> [(Prime Word, Word)]
factorialFactors n
| n < 2
= []
| otherwise
= map (\p -> (p, mult (unPrime p))) [minBound .. precPrime n]
where
mult :: Word -> Word
mult p = go np np
where
np = n `quot` p
go !acc !x
| x >= p = let xp = x `quot` p in go (acc + xp) xp
| otherwise = acc
-- | @'fibonacci' k@ calculates the @k@-th Fibonacci number in
-- /O/(@log (abs k)@) steps. The index may be negative. This
-- is efficient for calculating single Fibonacci numbers (with
-- large index), but for computing many Fibonacci numbers in
-- close proximity, it is better to use the simple addition
-- formula starting from an appropriate pair of successive
-- Fibonacci numbers.
fibonacci :: Num a => Int -> a
fibonacci = fst . fibonacciPair
{-# SPECIALIZE fibonacci :: Int -> Int #-}
{-# SPECIALIZE fibonacci :: Int -> Word #-}
{-# SPECIALIZE fibonacci :: Int -> Integer #-}
{-# SPECIALIZE fibonacci :: Int -> Natural #-}
-- | @'fibonacciPair' k@ returns the pair @(F(k), F(k+1))@ of the @k@-th
-- Fibonacci number and its successor, thus it can be used to calculate
-- the Fibonacci numbers from some index on without needing to compute
-- the previous. The pair is efficiently calculated
-- in /O/(@log (abs k)@) steps. The index may be negative.
fibonacciPair :: Num a => Int -> (a, a)
fibonacciPair n
| n < 0 = let (f,g) = fibonacciPair (-(n+1)) in if testBit n 0 then (g, -f) else (-g, f)
| n == 0 = (0, 1)
| otherwise = look (finiteBitSize (0 :: Word) - 2)
where
look k
| testBit n k = go (k-1) 0 1
| otherwise = look (k-1)
go k g f
| k < 0 = (f, f+g)
| testBit n k = go (k-1) (f*(f+shiftL1 g)) ((f+g)*shiftL1 f + g*g)
| otherwise = go (k-1) (f*f+g*g) (f*(f+shiftL1 g))
{-# SPECIALIZE fibonacciPair :: Int -> (Int, Int) #-}
{-# SPECIALIZE fibonacciPair :: Int -> (Word, Word) #-}
{-# SPECIALIZE fibonacciPair :: Int -> (Integer, Integer) #-}
{-# SPECIALIZE fibonacciPair :: Int -> (Natural, Natural) #-}
-- | @'lucas' k@ computes the @k@-th Lucas number. Very similar
-- to @'fibonacci'@.
lucas :: Num a => Int -> a
lucas = fst . lucasPair
{-# SPECIALIZE lucas :: Int -> Int #-}
{-# SPECIALIZE lucas :: Int -> Word #-}
{-# SPECIALIZE lucas :: Int -> Integer #-}
{-# SPECIALIZE lucas :: Int -> Natural #-}
-- | @'lucasPair' k@ computes the pair @(L(k), L(k+1))@ of the @k@-th
-- Lucas number and its successor. Very similar to @'fibonacciPair'@.
lucasPair :: Num a => Int -> (a, a)
lucasPair n
| n < 0 = let (f,g) = lucasPair (-(n+1)) in if testBit n 0 then (-g, f) else (g, -f)
| n == 0 = (2, 1)
| otherwise = look (finiteBitSize (0 :: Word) - 2)
where
look k
| testBit n k = go (k-1) 0 1
| otherwise = look (k-1)
go k g f
| k < 0 = (shiftL1 g + f,g+3*f)
| otherwise = go (k-1) g' f'
where
(f',g')
| testBit n k = (shiftL1 (f*(f+g)) + g*g,f*(shiftL1 g + f))
| otherwise = (f*(shiftL1 g + f),f*f+g*g)
{-# SPECIALIZE lucasPair :: Int -> (Int, Int) #-}
{-# SPECIALIZE lucasPair :: Int -> (Word, Word) #-}
{-# SPECIALIZE lucasPair :: Int -> (Integer, Integer) #-}
{-# SPECIALIZE lucasPair :: Int -> (Natural, Natural) #-}
-- | @'generalLucas' p q k@ calculates the quadruple @(U(k), U(k+1), V(k), V(k+1))@
-- where @U(i)@ is the Lucas sequence of the first kind and @V(i)@ the Lucas
-- sequence of the second kind for the parameters @p@ and @q@, where @p^2-4q /= 0@.
-- Both sequences satisfy the recurrence relation @A(j+2) = p*A(j+1) - q*A(j)@,
-- the starting values are @U(0) = 0, U(1) = 1@ and @V(0) = 2, V(1) = p@.
-- The Fibonacci numbers form the Lucas sequence of the first kind for the
-- parameters @p = 1, q = -1@ and the Lucas numbers form the Lucas sequence of
-- the second kind for these parameters.
-- Here, the index must be non-negative, since the terms of the sequence for
-- negative indices are in general not integers.
generalLucas :: Num a => a -> a -> Int -> (a, a, a, a)
generalLucas p q k
| k < 0 = error "generalLucas: negative index"
| k == 0 = (0,1,2,p)
| otherwise = look (finiteBitSize (0 :: Word) - 2)
where
look i
| testBit k i = go (i-1) 1 p p q
| otherwise = look (i-1)
go i un un1 vn qn
| i < 0 = (un, un1, vn, p*un1 - shiftL1 (q*un))
| testBit k i = go (i-1) (un1*vn-qn) ((p*un1-q*un)*vn - p*qn) ((p*un1 - (2*q)*un)*vn - p*qn) (qn*qn*q)
| otherwise = go (i-1) (un*vn) (un1*vn-qn) (vn*vn - 2*qn) (qn*qn)
{-# SPECIALIZE generalLucas :: Int -> Int -> Int -> (Int, Int, Int, Int) #-}
{-# SPECIALIZE generalLucas :: Word -> Word -> Int -> (Word, Word, Word, Word) #-}
{-# SPECIALIZE generalLucas :: Integer -> Integer -> Int -> (Integer, Integer, Integer, Integer) #-}
{-# SPECIALIZE generalLucas :: Natural -> Natural -> Int -> (Natural, Natural, Natural, Natural) #-}
shiftL1 :: Num a => a -> a
shiftL1 n = n + n
|
