Options

Fun With Recurrences

The following puzzles were inspired by John's Puzzle 104 and Puzzle 105 in Lecture 36.

Puzzle MD 1:

From the New York Times (2013)

There are 10 steps in front of you that you are about to climb. At any point, you can either take one step up, or you can jump two steps up. In how many unique ways can you climb the 10 steps?

Puzzle MD 2: Suppose there were \(10^7\) steps in front of you. Let \(N\) be the number of unique ways you can go up these steps. What are the last 10 digits of \(N\)?

(\(\star\)) What about the last 10 digits after \(10^{100}\) steps?

Puzzle MD 3: Suppose you could either go up the stairs 3 different ways: using your left foot to go up one, your right foot to go up one, or hop up two at a time. There's still \(10^7\) steps in front of you, what are the last 10 digits of the count now?

(\(\star\)) What about the last 10 digits after \(10^{100}\) steps?

Puzzle MD 4: Now suppose you can go either 1 step, 2 steps, or 3 steps at a time. Once again, what are the last 10 digits of the count for \(10^7\) steps?

(\(\star\)) What about the last 10 digits after \(10^{100}\) steps?

Puzzle MD 5: Finally, assume you can go up either 1 step, 2 steps, or two different ways of going 5 steps, or three different ways of going 9 steps. What are the last 10 digits for the count for \(10^7\) steps?

(\(\star\)) What about the last 10 digits after \(10^{100}\) steps?

Comments

  • 1.
    edited June 4

    Puzzle MD 1.

    For one step there's 1 way.

    For two steps there's 2 ways.

    For three steps there's 3 ways.

    For five steps there's 5 ways.

    ...

    For \(n\) steps there's \(f_{n-1}\) where \(f_{n-1}\) is the \(n-1\)st [Fibonacci number].

    Here's a script that outputs the answer (I'm using nix-shell to run Haskell):

    #! /usr/bin/env nix-shell
    #! nix-shell -p "haskell.packages.ghc822.ghc" -i "cabal exec -- runghc"
    
    module Main where
    
    fib :: Int -> Integer
    fib n | n < 0     = 0
          | n == 0    = 1
          | otherwise = fib (n - 1) + fib (n - 2)
    
    main :: IO ()
    main = (putStrLn . show . fib) 10
    

    This outputs 89, which is the answer.

    Comment Source:**Puzzle MD 1**. For one step there's 1 way. For two steps there's 2 ways. For three steps there's 3 ways. For five steps there's 5 ways. ... For \\(n\\) steps there's \\(f_{n-1}\\) where \\(f_{n-1}\\) is the \\(n-1\\)st [Fibonacci number]. Here's a script that outputs the answer (I'm using [`nix-shell`](https://nixos.org/) to run Haskell): <pre> #! /usr/bin/env nix-shell #! nix-shell -p "haskell.packages.ghc822.ghc" -i "cabal exec -- runghc" module Main where fib :: Int -> Integer fib n | n &lt; 0 = 0 | n == 0 = 1 | otherwise = fib (n - 1) + fib (n - 2) main :: IO () main = (putStrLn . show . fib) 10 </pre> This outputs *89*, which is the answer.
  • 2.
    edited June 4

    Puzzle MD 2.

    This problem is a little trickier, since it requires an efficient means of calculating Fibonacci numbers.

    However, in Exercise 1.19 of SICP § 1.2.4 an \(\mathcal{O}(\mathrm{ln}(n))\) algorithm is provided.

    Below I've translated the answer of this problem into Haskell and used it for my puzzle.

    #! /usr/bin/env nix-shell
    #! nix-shell -p "haskell.packages.ghc822.ghc" -i "cabal exec -- runghc"
    
    module Main where
    
    fib :: Int -> Integer
    fib n = fibIter 1 0 0 1 n
      where
        fibIter a b p q n 
          | n <= 0    = a
          | even n    = fibIter a b (p*p + q*q) (2*p*q + q*q) (n `div` 2)
          | otherwise = fibIter (b*q + a*q + a*p) (b*p + a*q) p q (n - 1)
    
    main :: IO ()
    main = putStrLn . show $ fib (10^7) `mod` (10^10)
    

    This program prints out the answer which is 7120937501

    Comment Source:**Puzzle MD 2.** This problem is a little trickier, since it requires an efficient means of calculating Fibonacci numbers. However, in Exercise 1.19 of [SICP § 1.2.4](https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-11.html#%_sec_1.2.4) an \\(\mathcal{O}(\mathrm{ln}(n))\\) algorithm is provided. Below I've translated the answer of this problem into Haskell and used it for my puzzle. <pre> #! /usr/bin/env nix-shell #! nix-shell -p "haskell.packages.ghc822.ghc" -i "cabal exec -- runghc" module Main where fib :: Int -> Integer fib n = fibIter 1 0 0 1 n where fibIter a b p q n | n &lt;= 0 = a | even n = fibIter a b (p*p + q*q) (2*p*q + q*q) (n `div` 2) | otherwise = fibIter (b*q + a*q + a*p) (b*p + a*q) p q (n - 1) main :: IO () main = putStrLn . show $ fib (10^7) `mod` (10^10) </pre> This program prints out the answer which is *7120937501*
  • 3.

    thanks for sharing it ..following

    Comment Source:thanks for sharing it ..following
  • 4.

    Thanks Pierre!

    My last post had a typo, I just fixed it.

    So, while Haskell is my favorite language, I know not everyone is into it.

    Since you are reading along, would you prefer if I wrote in another language?

    I think I can manage Scala, Python, Clojure, Java, JavaScript, C, C++, or Mathematica.

    Comment Source:Thanks Pierre! My last post had a typo, I just fixed it. So, while Haskell is my favorite language, I know not everyone is into it. Since you are reading along, would you prefer if I wrote in another language? I think I can manage Scala, Python, Clojure, Java, JavaScript, C, C++, or Mathematica.
  • 5.
    edited June 4

    Puzzle MD 2: Since \(\lambda x . x \textrm{ mod } n\) gives a homomorphism from \(\mathbb{Z}\) to \(\mathbb{Z}_n\), we're allowed to just ignore all but the last 10 digits on each step of recursion with \(n = 10^{10}\).

    -- normal_fibs = 0 : 1 : zipWith (+) fibs (tail fibs)
    fibs = 0 : 1 : zipWith (\x y -> mod (x + y) (10 ^ 10)) fibs (tail fibs)
    fib n = fibs!!(n + 1)
    
    main :: IO ()
    main = (putStrLn . show . fib) (10 ^ 7)
    

    The online repl was too slow so I ran it in python and got 7120937501.

    Puzzle MD 3: Same idea, but now there are 2 ways to get to step n from step n-1.

    pells = 0 : 1 : zipWith (\x y -> mod (x + (2 * y)) (10 ^ 10)) pells (tail pells)
    

    Which gives 4570126209.

    Puzzle MD 4:

    tribs = 0 : 0 : 1 : zipWith3 (\x y z -> mod (x + y + z) (10 ^ 10)) tribs (tail tribs) (tail (tail tribs))
    trib n = tribs!!(n + 2)
    

    Which gives 8389777176.

    All of these probably have much faster solutions of the form round(\(c_0*c_1^n\)), but I could use the haskell practice. Thanks for the puzzles!

    Comment Source:**Puzzle MD 2**: Since \\(\lambda x . x \textrm{ mod } n\\) gives a homomorphism from \\(\mathbb{Z}\\) to \\(\mathbb{Z}_n\\), we're allowed to just ignore all but the last 10 digits on each step of recursion with \\(n = 10^{10}\\). -- normal_fibs = 0 : 1 : zipWith (+) fibs (tail fibs) fibs = 0 : 1 : zipWith (\x y -> mod (x + y) (10 ^ 10)) fibs (tail fibs) fib n = fibs!!(n + 1) main :: IO () main = (putStrLn . show . fib) (10 ^ 7) The online repl was too slow so I ran it in python and got 7120937501. **Puzzle MD 3**: Same idea, but now there are 2 ways to get to step n from step n-1. pells = 0 : 1 : zipWith (\x y -> mod (x + (2 * y)) (10 ^ 10)) pells (tail pells) Which gives 4570126209. **Puzzle MD 4**: tribs = 0 : 0 : 1 : zipWith3 (\x y z -> mod (x + y + z) (10 ^ 10)) tribs (tail tribs) (tail (tail tribs)) trib n = tribs!!(n + 2) Which gives 8389777176. All of these probably have much faster solutions of the form round(\\(c_0*c_1^n\\)), but I could use the haskell practice. Thanks for the puzzles!
  • 6.
    edited June 4

    Hi Matthew the closest you get about matlab .m scripts (python?) even better to me..best

    Comment Source:Hi Matthew the closest you get about matlab .m scripts (python?) even better to me..best
  • 7.

    Hey Alex!

    Great job.

    Based on your answers I made some challenge variations of these problems, where the number of steps is \(10^{100}\).

    All of these probably have much faster solutions of the form round(\(c_0*c_1^n\)), but I could use the haskell practice. Thanks for the puzzles!

    Yeah, there's a solution like this. While it is cute, it's not numerically stable. There is also a similar solution that can be derived taking the \(z\)-transform.

    Comment Source:Hey Alex! Great job. Based on your answers I made some challenge variations of these problems, where the number of steps is \\(10^{100}\\). > All of these probably have much faster solutions of the form round(\\(c_0*c_1^n\\)), but I could use the haskell practice. Thanks for the puzzles! Yeah, there's a solution like this. While it is cute, it's not numerically stable. There is also a similar solution that can be derived taking the \\(z\\)-transform.
  • 8.

    Hi Matthew the closest you get about matlab .m scripts (python?) even better to me..best

    I will try to post python solutions to MD 1, MD 2 and the new MD 2 \((\star)\) later today.

    Comment Source:> Hi Matthew the closest you get about matlab .m scripts (python?) even better to me..best I will try to post python solutions to **MD 1**, **MD 2** and the new **MD 2 \\((\star)\\)** later today.
  • 9.
    edited June 5

    Puzzle MD 2 (\(\star\)): The homomorphism from \(\mathbb{Z}\) to \(\mathbb{Z}_n\) still works when we throw in multiplication, so we can still still put mods in everywhere:

    fib :: Integer -> Integer
    fib n = fibIter 1 0 0 1 n
      where
        fibIter a b p q n 
          | n <= 0    = a
          | even n    = fibIter a b ((p*p+q*q)`mod`(10^10)) ((2*p*q+q*q)`mod`(10^10)) (n `div` 2)
          | otherwise = fibIter ((b*q+a*q+a*p)`mod`(10^10)) ((b*p+a*q)`mod`(10^10)) p q (n - 1)
    
    main :: IO ()
    main = putStrLn . show $ fib (10^100)
    

    Which gives 2460937501.

    Comment Source:**Puzzle MD 2 (\\(\star\\))**: The homomorphism from \\(\mathbb{Z}\\) to \\(\mathbb{Z}_n\\) still works when we throw in multiplication, so we can still still put mods in everywhere: fib :: Integer -> Integer fib n = fibIter 1 0 0 1 n where fibIter a b p q n | n <= 0 = a | even n = fibIter a b ((p*p+q*q)`mod`(10^10)) ((2*p*q+q*q)`mod`(10^10)) (n `div` 2) | otherwise = fibIter ((b*q+a*q+a*p)`mod`(10^10)) ((b*p+a*q)`mod`(10^10)) p q (n - 1) main :: IO () main = putStrLn . show $ fib (10^100) Which gives 2460937501.
  • 10.

    You are a machine, Alex!

    Great work.

    Comment Source:You are a machine, Alex! Great work.
  • 11.
    edited June 5

    Puzzle MD 3 (\(\star\)): Same idea, but now \(T_{pq}\) brings \(a,b\) to \(bq+2aq+ap, bp+aq\) as described in that excercise Matthew shared, and the math checks out:

    fib n = fibIter 1 0 0 1 n
      where
        fibIter a b p q n 
          | n <= 0    = a
          | even n    = fibIter a b ((p*p+q*q)`mod`(10^10)) ((2*p*q+2*q*q)`mod`(10^10)) (n `div` 2)
          | otherwise = fibIter ((b*q+2*a*q+a*p)`mod`(10^10)) ((b*p+a*q)`mod`(10^10)) p q (n - 1)
    

    Which gives 9020126209.

    Comment Source:**Puzzle MD 3 (\\(\star\\))**: Same idea, but now \\(T_{pq}\\) brings \\(a,b\\) to \\(bq+2aq+ap, bp+aq\\) as described in that excercise Matthew shared, and the math checks out: fib n = fibIter 1 0 0 1 n where fibIter a b p q n | n <= 0 = a | even n = fibIter a b ((p*p+q*q)`mod`(10^10)) ((2*p*q+2*q*q)`mod`(10^10)) (n `div` 2) | otherwise = fibIter ((b*q+2*a*q+a*p)`mod`(10^10)) ((b*p+a*q)`mod`(10^10)) p q (n - 1) Which gives 9020126209.
  • 12.
    edited June 5

    Hey Pierre,

    Here is a solution to Puzzle MD 1 in Python 3:

    #! /usr/bin/env nix-shell
    #! nix-shell -p python3 -i python3
    
    def fib(n):
        a = 1
        b = 0
        for _ in range(n):
            temp = a
            a = a + b
            b = temp
        return a
    
    if __name__ == "__main__":
        print(fib(10))
    

    This outputs 89 like the Haskell program in #1.

    The problem with this program is that it's slow.

    To run the loop for _ in range(n): ... it needs to do \(\mathcal{O}(n)\) steps!

    It can manage \(10^7\) steps okay. Here my solution to Puzzle MD 2:

    #! /usr/bin/env nix-shell
    #! nix-shell -p python3 -i python3
    
    def fib_mod(n, modulus):
        a = 1
        b = 0
        for _ in range(n):
            temp = a
            a = (a + b) % modulus
            b = temp
        return a
    
    if __name__ == "__main__":
        print(fib_mod(10**7, 10**10))
    

    This takes my computer about 2 seconds to get the answer 7120937501

    But to do MD 2 \((\star)\), it would need to do \(10^{100}\) steps. At the rate my computer is going, this would take around ≈ \(5×10^{79}\) times longer than age of the universe to compute the answer to that puzzle.

    Comment Source:Hey Pierre, Here is a solution to **Puzzle MD 1** in Python 3: <pre> #! /usr/bin/env nix-shell #! nix-shell -p python3 -i python3 def fib(n): a = 1 b = 0 for _ in range(n): temp = a a = a + b b = temp return a if __name__ == "__main__": print(fib(10)) </pre> This outputs *89* like the Haskell program in [#1](https://forum.azimuthproject.org/discussion/comment/18920/#Comment_18920). The problem with this program is that it's *slow*. To run the loop `for _ in range(n): ...` it needs to do \\(\mathcal{O}(n)\\) steps! It can manage \\(10^7\\) steps *okay*. Here my solution to **Puzzle MD 2**: <pre> #! /usr/bin/env nix-shell #! nix-shell -p python3 -i python3 def fib_mod(n, modulus): a = 1 b = 0 for _ in range(n): temp = a a = (a + b) % modulus b = temp return a if __name__ == "__main__": print(fib_mod(10**7, 10**10)) </pre> This takes my computer about 2 seconds to get the answer *7120937501* But to do **MD 2 \\((\star)\\)**, it would need to do \\(10^{100}\\) steps. At the rate my computer is going, this would take around ≈ \\(5×10^{79}\\) times longer than age of the universe to compute the answer to that puzzle.
  • 13.

    Here's a more "pythonic" version using pattern matching:

    def fib_mod(n, modulus):
      a, b = 1, 0
      for _ in range(n):
        a, b = (a + b) % modulus, a
      return a
    
    Comment Source:Here's a more "pythonic" version using pattern matching: def fib_mod(n, modulus): a, b = 1, 0 for _ in range(n): a, b = (a + b) % modulus, a return a
  • 14.
    edited June 5

    Hehe, I started writing the code idiomatically like that if you can believe it.

    But Pierre said that he prefers the code to be as close to matlab as possible. I was trying to write the code in a clear and simple way, albeit a little clumsily.

    Comment Source:Hehe, I started writing the code idiomatically like that if you can believe it. But Pierre said that he prefers the code to be as close to matlab as possible. I was trying to write the code in a clear and simple way, albeit a little clumsily.
  • 15.

    Thank you!

    Comment Source:Thank you!
  • 16.

    Puzzle MD 4 (\(\star\)): This makes more sense to me in Mathematica with matrix multiplication, but it doesn't agree with my haskell solution. M is a matrix and v is the vector {a, b, c}:

    trib[n0_,mod_]:=Block[{n=n0,M,v},
      M={{1,1,1},{1,0,0},{0,1,0}};
      v={1,0,0};
      While[n>0,
        If[EvenQ[n],
          M=Mod[M.M,mod];n=n/2,
          v=Mod[M.v,mod];n=n-1];
      ];
      v[[1]]
    ];
    
    trib[10^7,10^10] (* 7912206465 *)
    trib[10^100,10^10] (* 6198029313 *)
    
    Comment Source:**Puzzle MD 4 (\\(\star\\))**: This makes more sense to me in Mathematica with matrix multiplication, but it doesn't agree with my haskell solution. M is a matrix and v is the vector {a, b, c}: trib[n0_,mod_]:=Block[{n=n0,M,v}, M={{1,1,1},{1,0,0},{0,1,0}}; v={1,0,0}; While[n>0, If[EvenQ[n], M=Mod[M.M,mod];n=n/2, v=Mod[M.v,mod];n=n-1]; ]; v[[1]] ]; trib[10^7,10^10] (* 7912206465 *) trib[10^100,10^10] (* 6198029313 *)
  • 17.

    Hey Alex,

    Let me write this up in Haskell and see if I can corroborate your answers...

    Comment Source:Hey Alex, Let me write this up in Haskell and see if I can corroborate your answers...
  • 18.
    edited June 5

    Hey Alex,

    Here is my solution (using Haskell):

    #! /usr/bin/env nix-shell
    #! nix-shell -p "haskell.packages.ghc822.ghcWithPackages (pkgs: [pkgs.matrix])" -i "cabal exec -- runghc"
    
    module Main where
    
    import Control.Applicative ((<$>))
    import qualified Data.Matrix as Matrix
    
    dot :: Num a => Matrix.Matrix a -> Matrix.Matrix a -> Matrix.Matrix a
    dot = Matrix.multStd
    
    matrixExp :: Integral a => Matrix.Matrix a -> Integer -> a -> Matrix.Matrix a
    matrixExp m n modulus
      | n <= 0    = Matrix.identity (Matrix.nrows m)
      | odd n     = (`mod` modulus) <$> m `dot` (matrixExp m (n-1) modulus)
      | otherwise = (`mod` modulus) <$> mHalfN `dot` mHalfN
        where mHalfN = matrixExp m (n `div` 2) modulus
    
    trib :: Integer -> Integer -> Integer
    trib n modulus = Matrix.getElem 1 1 $ (matrixExp u n modulus) `dot` unitImpulse
      where
        u           = Matrix.fromLists [ [1, 1, 1]
                                       , [1, 0, 0]
                                       , [0, 1, 0] ]
        unitImpulse = Matrix.fromLists [ [1] 
                                       , [0]
                                       , [0] ]
    
    main :: IO ()
    main = do 
        putStrLn $ "MD 4: " ++ show (trib (10^7) (10^10))
        putStrLn $ "MD 4 (*): " ++ show (trib (10^100) (10^10))
    

    This code agrees with 7912206465 and 6198029313


    The \(U\) operator in the code above is the same as the one John Baez mentions in Lecture 36, comment #83.

    In digital signal processing, operators like \(U\) can sometimes be thought of as digital filters.

    Consider the special case where \(U\) characterizes linear difference equations such as the Fibbonaccis, the Tribonnacis, or the Pell numbers. When we evaluate this recurrence, we are looking at \(U\)'s impulse response. Specifically, impulse is characterized by a vector:

    $$ \delta = \begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix} $$ This is referred to as unit impulse in the digital signal processing literature and in my code.

    Now let \(v = U^n \delta\). Then \(v_{1,1}\) is the \(n\)th term in the recurrence!

    I know John Baez is an expert on applying category theory to circuit design. Hopefully he can stop by and provide a few insights on how this connects to analogue filters.

    Comment Source:Hey Alex, Here is my solution (using Haskell): <pre> #! /usr/bin/env nix-shell #! nix-shell -p "haskell.packages.ghc822.ghcWithPackages (pkgs: [pkgs.matrix])" -i "cabal exec -- runghc" module Main where import Control.Applicative ((<$>)) import qualified Data.Matrix as Matrix dot :: Num a => Matrix.Matrix a -> Matrix.Matrix a -> Matrix.Matrix a dot = Matrix.multStd matrixExp :: Integral a => Matrix.Matrix a -> Integer -> a -> Matrix.Matrix a matrixExp m n modulus | n <= 0 = Matrix.identity (Matrix.nrows m) | odd n = (`mod` modulus) <$> m `dot` (matrixExp m (n-1) modulus) | otherwise = (`mod` modulus) <$> mHalfN `dot` mHalfN where mHalfN = matrixExp m (n `div` 2) modulus trib :: Integer -> Integer -> Integer trib n modulus = Matrix.getElem 1 1 $ (matrixExp u n modulus) `dot` unitImpulse where u = Matrix.fromLists [ [1, 1, 1] , [1, 0, 0] , [0, 1, 0] ] unitImpulse = Matrix.fromLists [ [1] , [0] , [0] ] main :: IO () main = do putStrLn $ "MD 4: " ++ show (trib (10^7) (10^10)) putStrLn $ "MD 4 (*): " ++ show (trib (10^100) (10^10)) </pre> This code agrees with *7912206465* and *6198029313* -------------------------------------------- The \\(U\\) operator in the code above is the same as the one John Baez mentions in [Lecture 36, comment #83](https://forum.azimuthproject.org/discussion/comment/19032/#Comment_19032). In digital signal processing, operators like \\(U\\) can sometimes be thought of as [*digital filters*](https://en.wikipedia.org/wiki/Digital_filter). Consider the special case where \\(U\\) characterizes [linear difference equations](https://en.wikipedia.org/wiki/Linear_difference_equation) such as the Fibbonaccis, the Tribonnacis, or the Pell numbers. When we evaluate this recurrence, we are looking at \\(U\\)'s [impulse response](https://en.wikipedia.org/wiki/Digital_filter#Impulse_response). Specifically, impulse is characterized by a vector: \[ \delta = \begin{bmatrix} 1 \\\\ 0 \\\\ \vdots \\\\ 0 \end{bmatrix} \] This is referred to as *unit impulse* in the digital signal processing literature and in my code. Now let \\(v = U^n \delta\\). Then \\(v_{1,1}\\) is the \\(n\\)th term in the recurrence! I know John Baez is an expert on applying category theory to circuit design. Hopefully he can stop by and provide a few insights on how this connects to analogue filters.
  • 19.
    edited June 9

    This realm of recurrences, discrete dinamics and Z-transform is very likable to stroll by and have productive reflections (and connected with generating functions). The case of continuous time and anolog filter desing is very much related. One focused on the Laplace transform instead of the Z one, but much algebra is preserved mutatis mutandis. I can't say I remember a lot from my college days, but from passive RLC circuits (with capacitors and inductances added to resistors and current and potential suources) one could easily formulate blocks with transfer functions in the Laplace domain that got composed by multiplication. There was also a formula for signal feeback. The global transfer function was a rational one, much connected with the ordinary diffential equations describing the behaviour. Poles and zeros were instrumental in the description. Baez and Jason Erbele studied the categorical angle here, but he could say better!

    Comment Source:This realm of recurrences, discrete dinamics and Z-transform is very likable to stroll by and have productive reflections (and connected with generating functions). The case of continuous time and anolog filter desing is very much related. One focused on the Laplace transform instead of the Z one, but much algebra is preserved mutatis mutandis. I can't say I remember a lot from my college days, but from passive RLC circuits (with capacitors and inductances added to resistors and current and potential suources) one could easily formulate [blocks](https://en.wikipedia.org/wiki/Laplace_transform#s-domain_equivalent_circuits_and_impedances) with transfer functions in the Laplace domain that got composed by multiplication. There was also a formula for signal feeback. The global transfer function was a rational one, much connected with the ordinary diffential equations describing the behaviour. Poles and zeros were instrumental in the description. Baez and Jason Erbele studied the categorical angle [here](https://arxiv.org/abs/1405.6881), but he could say better!
  • 20.
    edited June 13

    I thought I would take the chance to provide the general solution to all of these puzzles.

    If anyone has any questions I can motivate the solution here.

    #! /usr/bin/env nix-shell
    #! nix-shell -p "haskell.packages.ghc822.ghcWithPackages (pkgs: [pkgs.matrix])" -i "cabal exec -- runghc"
    
    module Main where
    
    import Control.Applicative ((<$>))
    import qualified Data.Matrix as Matrix
    import Data.Matrix ((<->))
    
    dot :: Num a => Matrix.Matrix a -> Matrix.Matrix a -> Matrix.Matrix a
    dot = Matrix.multStd
    
    matrixExp :: Integral a => Matrix.Matrix a -> Integer -> a -> Matrix.Matrix a
    matrixExp m n modulus
      | n <= 0    = Matrix.identity (Matrix.nrows m)
      | odd n     = (`mod` modulus) <$> m `dot` (matrixExp m (n-1) modulus)
      | otherwise = (`mod` modulus) <$> mHalfN `dot` mHalfN
        where mHalfN = matrixExp m (n `div` 2) modulus
    
    series :: [Integer] -> Integer -> Integer -> Integer
    series coeffs n modulus = Matrix.getElem 1 1 $ (matrixExp u n modulus)
      where
        m = length coeffs
        u = Matrix.fromLists [coeffs]
            <->
            (Matrix.extendTo 0 (m - 1) l $ Matrix.identity (m - 1))
    
    fib :: Integer -> Integer -> Integer
    fib = series [1,1]
    
    pell :: Integer -> Integer -> Integer
    pell = series [2,1]
    
    trib :: Integer -> Integer -> Integer
    trib = series [1,1,1]
    
    extreme :: Integer -> Integer -> Integer
    extreme = series [1,1,0,0,5,0,0,0,3]
    
    main :: IO ()
    main = do 
      putStrLn $ "MD 1:     " ++ show (fib 10 (10^10))
      putStrLn $ "MD 2:     " ++ show (fib (10^7) (10^10))
      putStrLn $ "MD 2 (*): " ++ show (fib (10^100) (10^10))
      putStrLn $ "MD 3:     " ++ show (pell (10^7) (10^10))
      putStrLn $ "MD 3 (*): " ++ show (pell (10^100) (10^10))
      putStrLn $ "MD 4:     " ++ show (trib (10^7) (10^10))
      putStrLn $ "MD 4 (*): " ++ show (trib (10^100) (10^10))
      putStrLn $ "MD 5:     " ++ show (extreme (10^7) (10^10))
      putStrLn $ "MD 5 (*): " ++ show (extreme (10^100) (10^10))
    

    This outputs:

    MD 1:     89
    MD 2:     7120937501
    MD 2 (*): 2460937501
    MD 3:     4570126209
    MD 3 (*): 9020126209
    MD 4:     7912206465
    MD 4 (*): 6198029313
    MD 5:     6777244134
    MD 5 (*): 4404226404
    
    Comment Source:I thought I would take the chance to provide the general solution to all of these puzzles. If anyone has any questions I can motivate the solution here. <pre> #! /usr/bin/env nix-shell #! nix-shell -p "haskell.packages.ghc822.ghcWithPackages (pkgs: [pkgs.matrix])" -i "cabal exec -- runghc" module Main where import Control.Applicative ((<$>)) import qualified Data.Matrix as Matrix import Data.Matrix ((<->)) dot :: Num a => Matrix.Matrix a -> Matrix.Matrix a -> Matrix.Matrix a dot = Matrix.multStd matrixExp :: Integral a => Matrix.Matrix a -> Integer -> a -> Matrix.Matrix a matrixExp m n modulus | n <= 0 = Matrix.identity (Matrix.nrows m) | odd n = (`mod` modulus) <$> m `dot` (matrixExp m (n-1) modulus) | otherwise = (`mod` modulus) <$> mHalfN `dot` mHalfN where mHalfN = matrixExp m (n `div` 2) modulus series :: [Integer] -> Integer -> Integer -> Integer series coeffs n modulus = Matrix.getElem 1 1 $ (matrixExp u n modulus) where m = length coeffs u = Matrix.fromLists [coeffs] <-> (Matrix.extendTo 0 (m - 1) l $ Matrix.identity (m - 1)) fib :: Integer -> Integer -> Integer fib = series [1,1] pell :: Integer -> Integer -> Integer pell = series [2,1] trib :: Integer -> Integer -> Integer trib = series [1,1,1] extreme :: Integer -> Integer -> Integer extreme = series [1,1,0,0,5,0,0,0,3] main :: IO () main = do putStrLn $ "MD 1: " ++ show (fib 10 (10^10)) putStrLn $ "MD 2: " ++ show (fib (10^7) (10^10)) putStrLn $ "MD 2 (*): " ++ show (fib (10^100) (10^10)) putStrLn $ "MD 3: " ++ show (pell (10^7) (10^10)) putStrLn $ "MD 3 (*): " ++ show (pell (10^100) (10^10)) putStrLn $ "MD 4: " ++ show (trib (10^7) (10^10)) putStrLn $ "MD 4 (*): " ++ show (trib (10^100) (10^10)) putStrLn $ "MD 5: " ++ show (extreme (10^7) (10^10)) putStrLn $ "MD 5 (*): " ++ show (extreme (10^100) (10^10)) </pre> This outputs: <pre> MD 1: 89 MD 2: 7120937501 MD 2 (*): 2460937501 MD 3: 4570126209 MD 3 (*): 9020126209 MD 4: 7912206465 MD 4 (*): 6198029313 MD 5: 6777244134 MD 5 (*): 4404226404 </pre>
  • 21.

    i don't really program or know category theory but as an aside this is similar to problems you get in statistical mechanics, and also 'econophysics'--application of stat mech to econ---exchanges among individuals . you get very different results depending on how you do the combinatorics, and the dynamic approaches are very different from static ones which assume you are already at equilibrium --dynamic approaches ask how get to the optimum. I'm unclear on this but sometimes it may be np-complete or noncomputable.

    Comment Source:i don't really program or know category theory but as an aside this is similar to problems you get in statistical mechanics, and also 'econophysics'--application of stat mech to econ---exchanges among individuals . you get very different results depending on how you do the combinatorics, and the dynamic approaches are very different from static ones which assume you are already at equilibrium --dynamic approaches ask how get to the optimum. I'm unclear on this but sometimes it may be np-complete or noncomputable.
  • 22.

    I didn't actually use any category theory here - this was inspired by some graph theory John presented in Lecture 36 when describing free categories.

    Comment Source:I didn't actually use any category theory here - this was inspired by some graph theory John presented in [Lecture 36](https://forum.azimuthproject.org/discussion/2204/lecture-36-categories-from-graphs) when describing free categories.
Sign In or Register to comment.