diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,27 @@
+Copyright (c) Henning Thielemann 2018
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. Neither the name of the author nor the names of his contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
diff --git a/src/Graph.hs b/src/Graph.hs
new file mode 100644
--- /dev/null
+++ b/src/Graph.hs
@@ -0,0 +1,93 @@
+module Main where
+
+import Text.Printf (printf)
+
+import qualified Data.Map as Map
+import Data.Map (Map)
+
+
+type Node = (Int, Int)
+
+pattern :: Node -> [Bool]
+pattern (a,b) = map (flip elem [a,b]) [0..5]
+
+ident, label :: Node -> String
+ident (a,b) =
+   map (\x -> if x then 'A' else '_') $ pattern (a,b)
+
+label (a,b) =
+   map (\x -> if x then ':' else '.') $ pattern (a,b)
+
+move :: [((Node, Node), String)]
+move =
+   (((0,0), (0,0)), "loop above") :
+   [(((a,b),(a-1,b-1)), "") | a<-[1..5], b<-[a..5]]
+
+roll :: [((Node, Node), String)]
+roll =
+   [(((0,1),(0,0)), "")]
+   ++
+   [(((0,1),(0,1)), "loop above")]
+   ++
+   [(((0,1),(0,a)), "out=-80,in=170") | a<-[2..5]]
+   ++
+   [((from,to), Map.findWithDefault "" (from,to) arcs) |
+       b<-[2..5], let b' = pred b, a<-[0..5],
+       let from = (0,b); to = (min a b', max a b')]
+
+arcs :: Map (Node, Node) String
+arcs =
+   Map.fromList $
+      (((0,2),(1,3)), "out=-10,in=-170") :
+      (((0,3),(2,3)), "out=45,in=-80") :
+      (((0,3),(2,5)), "out=-20,in=-160") :
+      (((0,4),(3,4)), "out=45,in=-80") :
+      (((0,5),(3,4)), "out=90,in=-80") :
+      (((0,5),(4,5)), "out=45,in=-80") :
+      []
+
+edge1 ::  ((Node, Node), String) -> String
+edge1 (((a,b),(c,d)), option) =
+   printf "(%s) edge [%s,dashed] (%s)" (ident (a,b)) option (ident (c,d))
+
+edge2 :: ((Node, Node), String) -> String
+edge2 (((a,b),(c,d)), option) =
+   printf "(%s) edge [%s] (%s)" (ident (a,b)) option (ident (c,d))
+
+pos :: Node -> String
+pos (a,b) =
+   printf "\\node[state%s] (%s) [at={(%d/2+%d*2,1.5*(%d-%d))}] {%s};"
+      (if a==0 && b/=0 then ",rolling" else "")
+      (ident (a,b)) a b a b (label (a,b))
+
+posCircularSector :: Node -> String
+posCircularSector p@(0,0) =
+   printf "\\node[state] (%s) [at={(0,0)}] {%s};" (ident p) (label p)
+posCircularSector (a,b) =
+   let af = fromIntegral a
+       bf = fromIntegral b
+       cos90, sin90 :: Float -> Float
+       cos90 x = cos (x/180*pi)
+       sin90 x = sin (x/180*pi)
+   in printf "\\node[state] (%s) [at={(%f,%f)}] {%s};"
+         (ident (a,b))
+         (2*bf*cos90(60/bf*(af-bf)))
+         (2*bf*sin90(60/bf*(af-bf)))
+         (label (a,b))
+
+allNodes :: [Node]
+allNodes = [(a,b) | a<-[0..5],b<-[a..5]]
+
+graph :: String
+graph =
+   unlines $
+      "\\begin{tikzpicture}[auto,rolling/.style=double]" :
+      map pos allNodes ++
+      "\\path[->]" :
+      map edge1 move ++ map edge2 roll ++
+      ";" :
+      "\\end{tikzpicture}" :
+      []
+
+main :: IO ()
+main = putStr graph
diff --git a/src/Matrix.hs b/src/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/src/Matrix.hs
@@ -0,0 +1,268 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Main where
+
+import qualified Numeric.LAPACK.Orthogonal as Ortho
+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Layout as Layout
+import qualified Numeric.LAPACK.Matrix.Square as Square
+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
+import qualified Numeric.LAPACK.Matrix as Matrix
+import qualified Numeric.LAPACK.Vector as Vector
+import qualified Numeric.LAPACK.Shape as ExtShape
+import qualified Numeric.Netlib.Class as Class
+import Numeric.LAPACK.Matrix
+         (ShapeInt,
+          (#-#), (#*#), (#\##), (#*##), (##*#), (#*|), (#*\), (\*#), (|||))
+
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Boxed as BoxedArray
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Shape ((::+))
+import Data.Array.Comfort.Storable ((!))
+
+import qualified Data.Complex as Complex
+import Data.Complex (Complex)
+
+import Text.Printf (printf)
+
+import qualified Data.List.Key as Key
+import qualified Data.List as List
+import Data.Semigroup ((<>))
+import Data.Tuple.HT (swap)
+
+
+type ShapeInt2 = ShapeInt ::+ ShapeInt
+
+data Score = P1 | P2 | P3 | P4 | P5 | P6 deriving (Eq, Ord, Enum, Bounded, Show)
+type ScoreShape = Shape.Enumeration Score
+type TriShape = Shape.UpperTriangular ScoreShape
+type TriShapeInt = ExtShape.IntIndexed TriShape
+type Transition = Matrix.Square TriShape
+type Transformation = Matrix.Square TriShapeInt
+type TallMatrix = Matrix.Tall TriShape ShapeInt
+type Matrix = Matrix.General TriShape TriShape
+type Vector = Vector.Vector TriShape
+type Eigenvalues = Vector.Vector TriShapeInt
+
+
+triShape :: TriShape
+triShape = Shape.upperTriangular Shape.Enumeration
+
+transition :: Transition Double
+transition =
+   Square.fromFull $ Matrix.scale (1/6) $
+   Matrix.fromRowArray triShape $ BoxedArray.fromList triShape $
+   map (Vector.fromList triShape) $
+      [6,1,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :
+      [0,1,1,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0] :
+      [0,1,0,1,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0,0,0] :
+      [0,1,0,0,1,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0,0] :
+      [0,1,0,0,0,1,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0] :
+      [0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :
+      [0,0,1,0,0,0,0,0,0,0,0,6,0,0,0,0,0,0,0,0,0] :
+      [0,0,1,1,0,0,0,0,0,0,0,0,6,0,0,0,0,0,0,0,0] :
+      [0,0,1,0,1,0,0,0,0,0,0,0,0,6,0,0,0,0,0,0,0] :
+      [0,0,1,0,0,1,0,0,0,0,0,0,0,0,6,0,0,0,0,0,0] :
+      [0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :
+      [0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,6,0,0,0,0,0] :
+      [0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,6,0,0,0,0] :
+      [0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,6,0,0,0] :
+      [0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :
+      [0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,0] :
+      [0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0] :
+      [0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :
+      [0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6] :
+      [0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :
+      [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :
+      []
+
+-- cf. Graph.move, roll
+transitionAlt :: Transition Double
+transitionAlt =
+   Matrix.scale (1/6) $ Matrix.transpose $ ArrMatrix.fromVector $
+   Array.fromAssociations 0 (Layout.square MatrixShape.RowMajor triShape) $
+      ((((P1,P1),(P1,P1)), 6) :
+       [(((a,b),(pred a, pred b)), 6) | a<-[P2 .. P6], b<-[a .. P6]] ++
+       [((from,to), 1) |
+         b<-[P2 .. P6], let b' = pred b, a<-[P1 .. P6],
+         let from = (P1,b); to = (min a b', max a b')])
+
+initial1 :: Vector Double
+initial1 =
+   Vector.scale (1/6) $
+   Vector.fromList triShape [1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
+
+initial2 :: Vector Double
+initial2 =
+   Vector.scale (1/36) $
+   Vector.fromList triShape [1,2,2,2,2,2,1,2,2,2,2,1,2,2,2,1,2,2,1,2,1]
+
+eigenSystem ::
+   (a ~ Complex Double) =>
+   (Transformation a, Eigenvalues a, Transformation a)
+eigenSystem =
+   let (vr,d,vlAdj) =
+         Square.eigensystem $ Square.mapSize ExtShape.IntIndexed transition
+       scal = Square.takeDiagonal $ vlAdj <> vr
+   in (vr, d, Vector.recip scal \*# vlAdj)
+
+similarityTransformation ::
+   (Class.Floating a, Vector.RealOf a ~ ar, Class.Real ar) =>
+   (Transformation a, Eigenvalues a, Transformation a) -> Transition a
+similarityTransformation (trafo, diag, trafoInv) =
+   Square.mapSize ExtShape.deconsIntIndexed $ trafo #*\ diag ##*# trafoInv
+
+reconstruct :: Transition (Complex Double)
+reconstruct = similarityTransformation eigenSystem
+
+reconstructionError :: (Complex Double, ((Score,Score), (Score,Score)))
+reconstructionError =
+   swap $ Vector.argAbsMaximum $ ArrMatrix.toVector $
+   Matrix.sub (Matrix.fromReal transition) reconstruct
+
+
+unwrapEigenvector :: Vector.Vector TriShapeInt a -> Vector a
+unwrapEigenvector = Array.mapShape ExtShape.deconsIntIndexed
+
+eigenvectors :: [(Complex Double, Vector (Complex Double))]
+eigenvectors =
+   let (m,v,_minv) = eigenSystem
+   in  Key.sort (negate . Complex.magnitude . fst) $
+       zip (Vector.toList v) (map unwrapEigenvector $ Matrix.toColumns m)
+
+tallMatrixFromColumns :: [Vector Double] -> TallMatrix Double
+tallMatrixFromColumns =
+   either id (error "matrix not tall") . Matrix.caseTallWide .
+   Matrix.fromColumns (Matrix.height transition)
+
+dominantEigenmatrix :: TallMatrix Double
+dominantEigenmatrix =
+   tallMatrixFromColumns $
+   map (Array.map Complex.realPart . snd) $ take 2 eigenvectors
+
+{-
+This is like the first step of the Schur decomposition.
+
+https://cs.stackexchange.com/questions/33068/partially-diagonalizing-a-matrix
+https://math.stackexchange.com/questions/3967036/partial-eigendecomposition-of-a-positive-semi-definite-matrix
+-}
+partialEigenSystem ::
+   (Matrix.LiberalSquare TriShape ShapeInt2 Double,
+    Matrix.LiberalSquare TriShape ShapeInt2 Double)
+partialEigenSystem =
+   let (ml,mr) = dominantEigenmatrixPairs
+       scal = Square.takeDiagonal $ Square.fromFull $ Matrix.adjoint ml #*# mr
+       mlScal = Matrix.generalizeTall $ ml #*\ Vector.recip scal
+       ortho = Matrix.generalizeTall $ Ortho.complement ml
+   in (Square.liberalFromFull $
+         mlScal  |||  ortho #-# mlScal #*# (Matrix.adjoint mr #*# ortho),
+       Square.liberalFromFull $ Matrix.generalizeTall mr ||| ortho)
+
+partiallyDiagonalized :: Matrix.LiberalSquare ShapeInt2 ShapeInt2 Double
+partiallyDiagonalized =
+   case fromInteger 1 :: Int of
+      0 -> let mr = snd partialEigenSystem in mr #\## transition #*## mr
+      _ ->
+         let (ml,mr) = partialEigenSystem
+         in Matrix.adjoint ml #*## transition #*## mr
+
+{-
+U are the given eigenvectors.
+Uorth is the basis for their orthogonal complement.
+
+U _|_ Uorth
+M*U _|_ M*Uorth - does not apply
+
+Thus, (U | Uorth) cannot diagonalize the transition matrix.
+-}
+
+
+filterEig :: (Double -> Bool) -> Transition (Complex Double)
+filterEig p =
+   let (trafo, diag, trafoInv) = eigenSystem
+   in similarityTransformation
+         (trafo,
+          Array.map (\x -> if p (Complex.magnitude x) then 1 else 0) diag,
+          trafoInv)
+
+powerCoeff :: Complex Double
+powerCoeff =
+   (filterEig (\x->0.99<x && x<1.01) #*| Vector.fromReal initial2) ! (P1,P1)
+
+
+outerProduct ::
+   (Eq height, Eq width, Shape.C height, Shape.C width, Class.Floating a) =>
+   Vector.Vector height a ->
+   Vector.Vector width a ->
+   Matrix.General height width a
+outerProduct = Matrix.outer MatrixShape.RowMajor
+
+rank2matrix :: Matrix.General ShapeInt ShapeInt Double
+rank2matrix =
+   Matrix.add
+      (Vector.autoFromList [38,37,20,80,38]
+       `outerProduct`
+       Vector.autoFromList [97,30,19,49,37])
+      (Vector.autoFromList [7,43,89,64,20]
+       `outerProduct`
+       Vector.autoFromList [4,96,44,49,80])
+
+eigenvectorPairs ::
+   [(Complex Double, (Vector (Complex Double), Vector (Complex Double)))]
+eigenvectorPairs =
+   Key.sort (negate . Complex.magnitude . fst) $
+   (\(vr,d,vlAdj) ->
+      zip (Vector.toList d) $
+      zip (map unwrapEigenvector $ Matrix.toRows vlAdj)
+          (map unwrapEigenvector $ Matrix.toColumns vr)) $
+   Square.eigensystem $ Square.mapSize ExtShape.IntIndexed transition
+
+dominantEigenmatrixPairs :: (TallMatrix Double, TallMatrix Double)
+dominantEigenmatrixPairs =
+   (tallMatrixFromColumns $ map (Array.map Complex.realPart . fst . snd) $
+      take 2 eigenvectorPairs,
+    tallMatrixFromColumns $ map (Array.map Complex.realPart . snd . snd) $
+      take 2 eigenvectorPairs)
+
+rankReduced :: Matrix (Complex Double)
+rankReduced =
+   foldl1 Matrix.add $
+   map
+      (\(l,(vl,vr)) ->
+         Matrix.scale (l / Vector.inner vl vr) $ outerProduct vr vl)
+      (take 2 eigenvectorPairs)
+
+powerParameters :: Vector Double -> [(Complex Double, Complex Double)]
+powerParameters initial =
+   let initialC = Vector.fromReal initial
+   in map
+         (\(l,(vl,vr)) ->
+            (l, (vr!(P1,P1)) * Vector.inner vl initialC / Vector.inner vl vr))
+         eigenvectorPairs
+
+
+probabilities :: Vector Double -> [Double]
+probabilities = map (! (P1,P1)) . iterate (transition #*|)
+
+probabilitiesApprox1 :: [Double]
+probabilitiesApprox1 =
+   map (1-) $ iterate (0.9083578427694454*) 0.9237265161044081
+
+probabilitiesApprox2 :: [Double]
+probabilitiesApprox2 =
+   map (1-) $ iterate (0.9083578427694454*) 1.025929108036849
+
+formatProbabilities :: [String]
+formatProbabilities =
+   List.zipWith5
+      (printf "%3d & %1.5f & %1.5f & %1.5f & %1.5f \\\\")
+      [0::Int ..]
+      (probabilities initial1)
+      probabilitiesApprox1
+      (probabilities initial2)
+      probabilitiesApprox2
+
+
+main :: IO ()
+main = mapM_ putStrLn $ take 101 formatProbabilities
diff --git a/wuerfelschlange.cabal b/wuerfelschlange.cabal
new file mode 100644
--- /dev/null
+++ b/wuerfelschlange.cabal
@@ -0,0 +1,60 @@
+Name:             wuerfelschlange
+Version:          0.1
+License:          BSD3
+License-File:     LICENSE
+Author:           Henning Thielemann <haskell@henning-thielemann.de>
+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>
+Homepage:         http://code.henning-thielemann.de/wuerfelschlange/
+Category:         Math
+Synopsis:         Code for the dice chain problem
+Description:
+  Code for the dice chain problem
+  <https://de.wikipedia.org/wiki/Würfelschlange>.
+  It is similar to the Kruskal Count.
+Tested-With:      GHC==7.4.2, GHC==7.8.4, GHC==8.6.5
+Cabal-Version:    1.14
+Build-Type:       Simple
+
+Flag buildGraph
+  description: Build graph generator
+  default:     True
+
+Flag buildMatrix
+  description: Build matrix computations
+  default:     True
+
+Source-Repository this
+  Tag:         0.1
+  Type:        darcs
+  Location:    http://code.henning-thielemann.de/wuerfelschlange/
+
+Source-Repository head
+  Type:        darcs
+  Location:    http://code.henning-thielemann.de/wuerfelschlange/
+
+Executable wuerfelschlange-graph
+  If flag(buildGraph)
+    Build-Depends:
+      containers >=0.4 && <0.7,
+      base >=4.5 && <5
+  Else
+    Buildable: False
+
+  Main-Is:          src/Graph.hs
+  GHC-Options:      -Wall
+  Default-Language: Haskell98
+
+Executable wuerfelschlange-matrix
+  If flag(buildMatrix)
+    Build-Depends:
+      lapack >=0.4 && <0.6,
+      netlib-ffi >=0.1.1 && <0.2,
+      comfort-array >=0.5 && <0.6,
+      utility-ht >=0.0.13 && <0.1,
+      base >=4.5 && <5
+  Else
+    Buildable: False
+
+  Main-Is:          src/Matrix.hs
+  GHC-Options:      -Wall
+  Default-Language: Haskell98