diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2015, Henning Thielemann
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * 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.
+
+    * Neither the name of Henning Thielemann nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS 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 COPYRIGHT
+OWNER 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/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#!/usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/linear-circuit.cabal b/linear-circuit.cabal
new file mode 100644
--- /dev/null
+++ b/linear-circuit.cabal
@@ -0,0 +1,56 @@
+Name:                linear-circuit
+Version:             0.0
+Synopsis:            Compute resistance of linear electrical circuits
+Description:
+  Compute resistance of linear electrical circuits.
+  .
+  For examples see test directory.
+Homepage:            http://hub.darcs.net/thielema/linear-circuit
+License:             BSD3
+License-File:        LICENSE
+Author:              Henning Thielemann
+Maintainer:          haskell@henning-thielemann.de
+Category:            Math
+Build-Type:          Simple
+Cabal-Version:       >=1.10
+
+Source-Repository this
+  Tag:         0.0
+  Type:        darcs
+  Location:    http://hub.darcs.net/thielema/linear-circuit
+
+Source-Repository head
+  Type:        darcs
+  Location:    http://hub.darcs.net/thielema/linear-circuit
+
+Library
+  Exposed-Modules:
+    Math.LinearCircuit
+  Build-Depends:
+    comfort-graph >=0.0 && <0.1,
+    hmatrix >=0.16 && <0.17,
+    containers >=0.4 && <0.6,
+    utility-ht >=0.0.11 && <0.1,
+    base >=4.5 && <5
+  Hs-Source-Dirs:      src
+  Default-Language:    Haskell2010
+  GHC-Options:         -Wall
+
+Test-Suite test-linear-circuit
+  Type:                exitcode-stdio-1.0
+  Hs-Source-Dirs:      test
+  Main-is:             Main.hs
+  Other-Modules:
+    ResistorCube
+    Tree
+  Build-Depends:
+    linear-circuit,
+    QuickCheck >=2 && <3,
+    comfort-graph,
+    non-empty >0.2 && <0.4,
+    transformers >=0.4 && <0.5,
+    containers,
+    utility-ht,
+    base
+  Default-Language:    Haskell2010
+  GHC-Options:         -Wall
diff --git a/src/Math/LinearCircuit.hs b/src/Math/LinearCircuit.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/LinearCircuit.hs
@@ -0,0 +1,96 @@
+module Math.LinearCircuit (resistance) where
+
+import qualified Data.Graph.Comfort as Graph
+import Data.Graph.Comfort (Graph)
+
+import qualified Numeric.Container as NC
+import qualified Numeric.LinearAlgebra.HMatrix as HMatrix
+import qualified Data.Packed.Matrix as Matrix
+import qualified Data.Packed.Vector as Vector
+import Numeric.LinearAlgebra.HMatrix (Field, (<\>))
+import Data.Packed.Matrix (Matrix)
+import Data.Packed.Vector (Vector)
+
+import qualified Data.Map as Map
+import qualified Data.List as List
+import Data.Monoid (mconcat)
+
+import Control.Functor.HT (outerProduct)
+import Data.Bool.HT (if')
+
+
+voltageMatrix ::
+   (Graph.Edge edge, Ord node, Field a) =>
+   Graph edge node a nodeLabel -> Matrix a
+voltageMatrix gr =
+   Matrix.fromLists $
+   outerProduct
+      (\e n ->
+         if' (Graph.from e == n) 1 $
+         if' (Graph.to   e == n) (-1) $
+         0)
+      (Graph.edges gr)
+      (Graph.nodes gr)
+
+{- |
+It is almost currentMatrix = trans voltageMatrix,
+except that a row is deleted in currentMatrix.
+-}
+currentMatrix ::
+   (Graph.Edge edge, Ord node, Field a) =>
+   Graph edge node a nodeLabel -> node -> node -> Matrix a
+currentMatrix gr _src dst =
+   Matrix.fromLists $
+   outerProduct
+      (\n e ->
+         if' (Graph.from e == n) 1 $
+         if' (Graph.to   e == n) (-1) $
+         0)
+      (List.delete dst $ Graph.nodes gr)
+      (Graph.edges gr)
+
+resistanceMatrix ::
+   (Graph.Edge edge, Ord node, Field a) =>
+   Graph edge node a nodeLabel -> Matrix a
+resistanceMatrix gr =
+   HMatrix.diag $ Vector.fromList $
+   Map.elems $ Graph.edgeLabels gr
+
+fullMatrix ::
+   (Graph.Edge edge, Ord node, Field a) =>
+   Graph edge node a nodeLabel -> node -> node -> Matrix a
+fullMatrix gr src dst =
+   let currents = currentMatrix gr src dst
+       voltages = voltageMatrix gr
+   in  Matrix.fromBlocks
+          [[NC.konst 0 (1, Matrix.cols currents),
+               Matrix.asRow $ Vector.fromList $
+               map (\n -> if n==src then 1 else 0) $ Graph.nodes gr],
+           [resistanceMatrix gr, voltages],
+           [currents, NC.konst 0 (Matrix.rows currents, Matrix.cols voltages)]]
+
+rhs ::
+   (Graph.Edge edge, Ord node, Field a) =>
+   Graph edge node a nodeLabel -> node -> node -> Vector a
+rhs gr src dst =
+   mconcat
+      [NC.konst 0 1,
+       NC.konst 0 (length (Graph.edges gr)),
+       Vector.fromList $
+       map (\n -> if n==src then 1 else 0) $
+       List.delete dst $ Graph.nodes gr]
+
+
+solution ::
+   (Graph.Edge edge, Ord node, Field a) =>
+   Graph edge node a nodeLabel -> node -> node -> Vector a
+solution gr src dst =
+   fullMatrix gr src dst <\> rhs gr src dst
+
+resistance ::
+   (Graph.Edge edge, Ord node, Field a) =>
+   Graph edge node a nodeLabel -> node -> node -> a
+resistance gr src dst =
+   solution gr src dst
+   `NC.atIndex`
+   (length (Graph.edges gr) + length (takeWhile (dst/=) $ Graph.nodes gr))
diff --git a/test/Main.hs b/test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Main.hs
@@ -0,0 +1,22 @@
+module Main where
+
+import qualified ResistorCube
+import qualified Tree
+
+import qualified Test.QuickCheck as QC
+
+
+approx :: Double -> Double -> Bool
+approx x y = abs (x-y) < 1e-8
+
+test :: (QC.Testable prop) => String -> prop -> IO ()
+test msg prop =
+   putStr (msg ++ ": ") >> QC.quickCheck prop
+
+main :: IO ()
+main = do
+   test "resistor cube" (approx ResistorCube.resistance (5/6))
+   test "resistor tree"
+      (\x -> approx (Tree.treeResistance x) (Tree.graphResistance x))
+   test "orientation of resistors"
+      (uncurry approx . Tree.flippedResistances)
diff --git a/test/ResistorCube.hs b/test/ResistorCube.hs
new file mode 100644
--- /dev/null
+++ b/test/ResistorCube.hs
@@ -0,0 +1,41 @@
+{- |
+Consider a cube of resistors of equal resistance.
+What is the overall resistance from one corner to the opposite one?
+-}
+module ResistorCube where
+
+import qualified Math.LinearCircuit as LinearCircuit
+
+import qualified Data.Graph.Comfort as Graph
+import Data.Graph.Comfort (Graph)
+
+import Control.Applicative (liftA2, liftA3)
+
+
+data Coord = C0 | C1 deriving (Eq, Ord, Show, Enum, Bounded)
+data Corner = Corner Coord Coord Coord deriving (Eq, Ord, Show)
+
+
+allCoords :: [Coord]
+allCoords = [minBound .. maxBound]
+
+dimEdges ::
+   (Coord -> Coord -> Coord -> Corner) ->
+   [(Graph.UndirEdge Corner, Double)]
+dimEdges corner =
+   liftA2
+      (\a b -> (Graph.undirEdge (corner C0 a b) (corner C1 a b), 1))
+      allCoords allCoords
+
+graph :: Graph Graph.UndirEdge Corner Double ()
+graph =
+   Graph.fromList
+      (map (flip (,) ()) $ liftA3 Corner allCoords allCoords allCoords)
+      (dimEdges (\x y z -> Corner x y z) ++
+       dimEdges (\y z x -> Corner x y z) ++
+       dimEdges (\z x y -> Corner x y z))
+
+resistance :: Double
+resistance =
+   LinearCircuit.resistance graph
+      (Corner C0 C0 C0) (Corner C1 C1 C1)
diff --git a/test/Tree.hs b/test/Tree.hs
new file mode 100644
--- /dev/null
+++ b/test/Tree.hs
@@ -0,0 +1,184 @@
+{- |
+Arrange resistors according to a tree of parallel and serial compositions.
+Compare resistance of trees with the general graph resistance computation.
+-}
+module Tree where
+
+import qualified Math.LinearCircuit as LinearCircuit
+
+import qualified Test.QuickCheck as QC
+
+import qualified Data.Graph.Comfort as Graph
+import Data.Graph.Comfort (Graph)
+
+import qualified Control.Monad.Trans.Class as MT
+import qualified Control.Monad.Trans.State as MS
+import Control.Monad (liftM, liftM2, replicateM)
+
+import qualified Data.Map as Map; import Data.Map (Map)
+
+import qualified Data.NonEmpty.Class as NonEmptyC
+import qualified Data.NonEmpty as NonEmpty
+import qualified Data.Foldable as Fold
+import qualified Data.List.Match as Match
+import Data.Functor.Classes (Eq1, Ord1, Show1, eq1, compare1, showsPrec1)
+import Data.Monoid (mappend)
+import Data.Ord.HT (comparing)
+import Data.Eq.HT (equating)
+
+
+data T =
+     Resistance Double
+   | Serial (NonEmpty.T [] T)
+   | Parallel (NonEmpty.T [] T)
+   deriving (Show)
+
+instance QC.Arbitrary T where
+   arbitrary =
+      let res = liftM Resistance $ QC.choose (0,1)
+          go 0 = res
+          go size =
+            let subTree n =
+                  let x = QC.resize (div size n) QC.arbitrary
+                  in  liftM2 NonEmpty.cons x (replicateM (n-1) x)
+            in  QC.frequency $
+                  (3, res) :
+                  (1, liftM Serial   (QC.choose (1,size) >>= subTree)) :
+                  (1, liftM Parallel (QC.choose (1,size) >>= subTree)) :
+                  []
+      in  QC.sized go
+   shrink tree =
+      case tree of
+         Resistance res ->
+            let simpleRess = [0,1]
+            in  if elem res simpleRess
+                  then []
+                  else map Resistance simpleRess
+         Parallel xs -> NonEmpty.flatten xs ++ map Parallel (QC.shrink xs)
+         Serial xs   -> NonEmpty.flatten xs ++ map Serial   (QC.shrink xs)
+
+
+parallel2 :: Double -> Double -> Double
+parallel2 0 0 = 0
+parallel2 x y = x*y / (x+y)
+
+treeResistance :: T -> Double
+treeResistance x =
+   case x of
+      Resistance res -> res
+      Serial xs -> Fold.foldl1 (+) $ fmap treeResistance xs
+      Parallel xs -> Fold.foldl1 parallel2 $ fmap treeResistance xs
+
+
+newtype EdgeId = EdgeId Int
+   deriving (Eq, Ord, Show)
+
+instance Enum EdgeId where
+   fromEnum (EdgeId n) = n
+   toEnum = EdgeId
+
+newEdgeId :: (Monad m) => MS.StateT EdgeId m EdgeId
+newEdgeId = do
+   n <- MS.get
+   MS.put $ succ n
+   return n
+
+data Edge a =
+   Edge {
+      edgeId :: EdgeId,
+      edgeFrom, edgeTo :: a
+   }
+   deriving (Show)
+
+instance Eq (Edge a) where (==) = equating edgeId
+instance Ord (Edge a) where compare = comparing edgeId
+
+instance Eq1 Edge where eq1 = (==)
+instance Ord1 Edge where compare1 = compare
+instance Show1 Edge where showsPrec1 = showsPrec
+
+instance Fold.Foldable Edge where
+   foldMap f (Edge _ x y) = mappend (f x) (f y)
+
+instance Graph.Edge Edge where
+   from (Edge _ n _) = n
+   to (Edge _ _ n) = n
+
+instance Graph.Reverse Edge where
+   reverseEdge (Edge n from to) = Edge n to from
+
+
+newtype Node = Node Int
+   deriving (Eq, Ord, Show)
+
+instance Enum Node where
+   fromEnum (Node n) = n
+   toEnum = Node
+
+newNode :: (Monad m) => MS.StateT Node m Node
+newNode = do
+   n <- MS.get
+   MS.put $ succ n
+   return n
+
+edgesFromTree ::
+   T -> (Node, Node) ->
+   MS.StateT EdgeId (MS.State Node) (Map (Edge Node) Double)
+edgesFromTree tree (from, to) =
+   case tree of
+      Resistance res -> do
+         e <- newEdgeId
+         return $ Map.singleton (Edge e from to) res
+      Serial xs -> do
+         ns <- sequence $ Match.replicate (NonEmpty.tail xs) $ MT.lift newNode
+         fmap Map.unions $ sequence $
+            NonEmpty.flatten $
+            NonEmptyC.zipWith edgesFromTree xs $
+            NonEmpty.mapAdjacent (,) $
+            NonEmpty.cons from $ NonEmpty.snoc ns to
+      Parallel xs -> do
+         fmap Map.unions $ mapM (flip edgesFromTree (from,to)) $
+            NonEmpty.flatten xs
+
+graphFromTree :: T -> (Graph Edge Node Double (), (Node, Node))
+graphFromTree tree =
+   let ((edgeMap, globalEnds), lastNode) =
+         flip MS.runState (Node 0) $ flip MS.evalStateT (EdgeId 0) $ do
+            ends <- MT.lift $ liftM2 (,) newNode newNode
+            edges <- edgesFromTree tree ends
+            return (edges, ends)
+   in  (Graph.fromMap
+          (Map.fromList $ map (flip (,) ()) [Node 0 .. pred lastNode])
+          edgeMap,
+        globalEnds)
+
+graphResistance :: T -> Double
+graphResistance =
+   uncurry (uncurry . LinearCircuit.resistance) . graphFromTree
+
+
+
+data
+   FlippedGraph =
+      FlippedGraph (Graph Edge Node (Double, Bool) ()) (Node, Node)
+   deriving (Show)
+
+instance QC.Arbitrary FlippedGraph where
+   arbitrary = do
+      (graph, ends) <- fmap graphFromTree QC.arbitrary
+      flpGraph <-
+         Graph.traverseEdge (\res -> liftM ((,) res) QC.arbitrary) graph
+      return $ FlippedGraph flpGraph ends
+
+flippedResistances :: FlippedGraph -> (Double, Double)
+flippedResistances (FlippedGraph graph ends) =
+   let flippedGraph =
+          Graph.fromMap
+             (Graph.nodeLabels graph)
+             (Map.fromList $
+              map
+                 (\(e, (res, flp)) ->
+                    (if flp then Graph.reverseEdge e else e, res)) $
+              Map.toList $ Graph.edgeLabels graph)
+   in  (uncurry (LinearCircuit.resistance (Graph.mapEdge fst graph)) ends,
+        uncurry (LinearCircuit.resistance flippedGraph) ends)