diff --git a/examples/MM0916.hs b/examples/MM0916.hs
new file mode 100644
--- /dev/null
+++ b/examples/MM0916.hs
@@ -0,0 +1,110 @@
+-- |  http://www2.stetson.edu/~efriedma/mathmagic/0916.html
+-- On an N×N chessboard, when we place Q queens, 
+-- what is the maximum number of squares 
+-- that can be attacked exactly A times?
+
+{-# language LambdaCase #-}
+
+import Prelude hiding ((&&),(||),not,and,or)
+import qualified Prelude as P
+import OBDD 
+import OBDD.Linopt
+
+import Data.Ix (inRange)
+import qualified Data.Array as A
+import Control.Monad ( guard )
+import System.Environment ( getArgs )
+import Data.List (sort)
+import qualified Data.Map.Strict as M
+
+main = getArgs >>= \ case
+  [] -> run 6 2 0
+  [n,q,a] -> run (read n) (read q) (read a)  
+
+run n q a = putStrLn $ form n q a
+        $ linopt ( board n q a ) 
+        $ M.fromList 
+        $ zip ((\ p -> Var p Attacked) <$> positions n) (repeat 1)
+        ++ zip ((\ p -> Var p Queen) <$> positions n) (repeat 0)
+
+type Position = (Int,Int)
+
+positions :: Int -> [ Position ]
+positions n = (,) <$> [1..n] <*> [1..n]
+
+data Type =  Attacked | Queen deriving (Eq, Ord, Show)
+data Var = Var !Position !Type deriving (Eq, Ord, Show)
+
+type Bit = OBDD Var
+
+queen p = variable $ Var p Queen
+attacked p = variable $ Var p Attacked
+
+
+header n q a w = unwords 
+   [ "n =", show n
+   , "q =", show q
+   , "a =", show a
+   , "m =", show w 
+   ] 
+
+form n q a (Just (w,m)) = unlines $ header n q a w  : do
+  row <- [1..n]
+  return $ do
+    col <- [1..n]
+    let c = if m M.! Var (row,col) Queen then 'Q' 
+            else if m M.! Var (row,col) Attacked then '+'
+            else '.'
+    [ c, ' ' ]
+
+for = flip map
+
+board :: Int -> Int -> Int -> Bit
+board n q a = let r = ray n  in and 
+    $ (  exactly q $ queen <$> positions n )
+    : ( for ( positions n) $ \ p -> 
+      (not $ queen p) || (not $ attacked p) )
+    ++ ( for (positions n) $ \ p -> 
+     implies (attacked p) $ exactly a $ for directions $ \ d ->
+       r A.! (d,p) )
+    
+
+-- | ray n ! (d,p) == looking in direction d from p,
+--  there is (at least one) queen (which might be on p)
+ray n = 
+    let bounds = (((-1,-1),(1,1)),((1,1),(n,n)))
+        result =  A.array bounds $ do
+          (d,p) <- A.range bounds
+          let q = shift d p
+          return ( (d,p)
+              , queen p
+                  || if onboard n q then result A.! (d,q) else false
+              )
+    in  result
+  
+
+directions = filter (/= (0,0)) 
+   $ (,) <$> [ -1 .. 1 ] <*> [ -1 .. 1 ]
+
+onboard n (x,y) = inRange (1,n) x P.&& inRange (1,n) y
+shift (dx,dy) (x,y) = (x+dx,y+dy)
+
+exactly :: Int -> [Bit] -> Bit
+exactly k xs = 
+  if k <= 8 P.&& length xs <= 8 
+  then exactly_direct k xs
+  else exactly_rectangle k xs
+
+exactly_rectangle n xs = last $ 
+  foldl ( \ cs x -> zipWith ( \ a b -> ite x a b )
+                    (false : cs) cs 
+        ) (true : replicate n false) xs
+
+exactly_direct k xs = atmost k xs && atleast k xs
+
+atmost k xs = not $ atleast (k+1) xs
+atleast k xs = or $ for (select k xs) and
+
+select 0 xs = return []
+select k [] = []
+select k (x:xs) = select k xs ++ ( (x:) <$> select (k-1) xs )
diff --git a/obdd.cabal b/obdd.cabal
--- a/obdd.cabal
+++ b/obdd.cabal
@@ -1,10 +1,31 @@
 Name:                obdd
-Version:             0.5.0
+Version:             0.6.0
 Cabal-Version:       >= 1.8
 Build-type: Simple
 Synopsis:            Ordered Reduced Binary Decision Diagrams
-Description:         Construct, combine and query OBDDs;
-                     an efficient representation for formulas in propositional logic
+Description:
+  Construct, combine and query OBDDs;
+  an efficient representation for formulas in propositional logic.
+  .
+  This is mostly educational.
+  The BDDs do not share nodes and this might introduce inefficiencies.
+  .
+  An important (for me, in teaching) feature is
+  that I can immediately draw the BDD to an X11 window (via graphviz).
+  For example, to show the effect of different variable orderings,
+  try this in ghci:
+  .
+  > import qualified Prelude as P
+  > import OBDD
+  > let f [] = false; f (x:y:zs) = x && y || f zs
+  > display P.$ f P.$ P.map variable [1,2,3,4,5,6]
+  > display P.$ f P.$ P.map variable [1,4,2,5,3,6]
+
+  If you want better performance,
+  use <http://vlsi.colorado.edu/%7Efabio/CUDD/ CUDD>
+  <https://hackage.haskell.org/package/cudd Haskell bindings>,
+  see <https://gitlab.imn.htwk-leipzig.de/waldmann/min-comp-sort this example>.
+
 category:	     Logic
 License:             GPL
 License-file:        LICENSE
@@ -55,3 +76,9 @@
     Build-Depends: base, containers, obdd
     Ghc-Options: -rtsopts    
 
+test-suite obdd-mm0916
+    Hs-Source-Dirs : examples
+    Type: exitcode-stdio-1.0
+    Main-Is: MM0916.hs
+    Build-Depends: base, containers, obdd, array
+    Ghc-Options: -rtsopts
diff --git a/src/OBDD.hs b/src/OBDD.hs
--- a/src/OBDD.hs
+++ b/src/OBDD.hs
@@ -1,13 +1,9 @@
 -- | Reduced ordered binary decision diagrams,
 -- pure Haskell implementation.
--- (c) Johannes Waldmann, 2008
+-- (c) Johannes Waldmann, 2008 - 2016
 --
 -- This module is intended to be imported qualified
 -- because it overloads some Prelude names.
---
--- For a similar, but much more elaborate project, see
--- <http://www.informatik.uni-kiel.de/~mh/lehre/diplomarbeiten/christiansen.pdf>
--- but I'm not sure where that source code would be available.
 
 module OBDD 
 
diff --git a/src/OBDD/Data.hs b/src/OBDD/Data.hs
--- a/src/OBDD/Data.hs
+++ b/src/OBDD/Data.hs
@@ -1,6 +1,7 @@
 {-# language GeneralizedNewtypeDeriving #-}
 {-# language RecursiveDo #-}
 {-# language FlexibleContexts #-}
+{-# language TupleSections #-}
 
 -- | implementation of reduced ordered binary decision diagrams.
 
@@ -15,6 +16,7 @@
 , number_of_models
 , some_model, all_models
 , fold, foldM
+, full_fold, full_foldM
 , toDot, display
 -- * for internal use
 , Node (..)
@@ -43,6 +45,7 @@
 
 import Data.Set ( Set )
 import qualified Data.Set as S
+import Data.Bool (bool)
 
 import Control.Arrow ( (&&&) )
 import Control.Monad.State.Strict
@@ -51,6 +54,7 @@
 import Control.Monad.Fix
 import Control.Monad ( forM, guard, void )
 import qualified Control.Monad ( foldM )
+import Data.Functor.Identity
 import System.Process
 import Data.List (isPrefixOf, isSuffixOf)
 
@@ -76,23 +80,26 @@
                -- (unary will be simulated by binary)
             }
 
+
+-- | Apply function in each node, bottom-up.
+-- return the value in the root node.
+-- Will cache intermediate results.
+-- You might think that 
+-- @count_models = fold (\b -> if b then 1 else 0) (\v l r -> l + r)@ 
+-- but that's not true because a path might omit variables.
+-- Use @full_fold@ to fold over interpolated nodes as well.
 fold :: Ord v 
      => ( Bool -> a )
      -> ( v -> a -> a -> a )
      -> OBDD v -> a
-fold leaf branch o =
-    let f = leaf False ; t = leaf True
-        m0 = M.fromList 
-           [(icore_false,f), (icore_true,t)]
-        m = foldl ( \ m (i,n) -> 
-            let val = case n of
-                    Branch v l r -> 
-                        branch v (m M.! l) (m M.! r) 
-            in M.insert i val m
-          ) m0 $ IM.toAscList $ core o
-    in  m M.! top o
-
+fold leaf branch o = runIdentity 
+   $ foldM ( return . leaf )
+           ( \ v l r -> return $ branch v l r )
+           o
 
+-- | Run action in each node, bottum-up.
+-- return the value in the root node.
+-- Will cache intermediate results.
 foldM :: (Monad m, Ord v)
      => ( Bool -> m a )
      -> ( v -> a -> a -> m a )
@@ -109,7 +116,49 @@
           ) m0 $ IM.toAscList $ core o
     return $ m M.! top o
 
+-- | Apply function in each node, bottom-up.
+-- Also apply to interpolated nodes: when a link
+-- from a node to a child skips some variables:
+-- for each skipped variable, we run the @branch@ function
+-- on an interpolated node that contains this missing variable,
+-- and identical children.
+-- With this function, @number_of_models@
+-- can be implemented as 
+-- @full_fold vars (bool 0 1) ( const (+) )@.
+-- And it actually is, see the source.
+full_fold :: Ord v 
+     => Set v
+     -> ( Bool -> a )
+     -> ( v -> a -> a -> a )
+     -> OBDD v -> a
+full_fold vars leaf branch o = runIdentity 
+   $ full_foldM vars 
+           ( return . leaf )
+           ( \ v l r -> return $ branch v l r )
+           o
 
+full_foldM :: (Monad m, Ord v)
+     => Set v 
+     -> ( Bool -> m a )
+     -> ( v -> a -> a -> m a )
+     -> OBDD v -> m a
+full_foldM vars leaf branch o = do
+    let vs = S.toAscList vars
+        low = head vs
+        m = M.fromList $ zip vs $ tail vs
+        up v = M.lookup v m
+        interpolate now goal x | now == goal = return x
+        interpolate (Just now) goal x = 
+            branch now x x >>= interpolate (up now) goal
+    (a,res) <- foldM 
+          ( \ b -> (Just low ,) <$> leaf b )
+          ( \ v (p,l) (q,r) -> do
+                l' <- interpolate p (Just v) l
+                r' <- interpolate q (Just v) r
+                (up v,) <$> branch v l' r'
+          ) o
+    interpolate a Nothing res  
+
 icore_false = 0 :: Index  
 icore_true = 1 :: Index  
               
@@ -120,26 +169,8 @@
 -- all variables that were used to construct it, since some  nodes may have been removed
 -- because they had identical children.
 number_of_models :: Ord v => Set v -> OBDD v ->  Integer
-number_of_models vs o = 
-    let fun o vs = do
-            m <- get
-            case access o of
-                   Leaf c -> case c of
-                        False -> return 0
-                        True -> return $ 2 ^ length vs
-                   Branch v l r -> do
-                       let ( pre, _ : post ) = span (/= v) vs
-                       case M.lookup ( top o ) m of
-                          Just x -> return $ ( 2 ^ length pre ) * x
-                          Nothing -> do
-                             xl <- fun l post
-                             xr <- fun r post
-                             let xlr = xl + xr
-                             m <- get
-                             put $! M.insert ( top o ) xlr m
-                             return $ ( 2 ^ length pre ) * xlr
-    in evalState ( fun o $ reverse $ S.toAscList vs ) M.empty
-    
+number_of_models vars o = 
+  full_fold vars (bool 0 1) ( const (+) ) o
 
 empty :: OBDD v
 empty = OBDD 
@@ -197,20 +228,11 @@
 -- | list of all models (WARNING not using 
 -- variables that had been deleted)
 all_models :: Ord v => OBDD v -> [ Map v Bool ]
-all_models s = case access s of
-    Leaf True -> return  M.empty
-    Leaf False -> [ ]
-    Branch v l r -> do
-        let nonempty_children = do
-                 ( p, t ) <- [ (False, l), (True, r) ]        
-                 guard $ case access t of
-                      Leaf False -> False
-                      _ -> True
-                 return ( p, t )
-        (p, t) <- nonempty_children
-        m <- all_models t
-        return $ M.insert v p m 
-        
+all_models = 
+  fold ( bool [] [ M.empty ] )
+       ( \ v l r -> (M.insert v False <$> l)
+                 ++ (M.insert v True  <$> r) )
+
 select_one :: [a] -> IO a
 select_one xs | not ( Prelude.null xs ) = do
     i <- System.Random.randomRIO ( 0, length xs - 1 )
diff --git a/src/OBDD/Linopt.hs b/src/OBDD/Linopt.hs
--- a/src/OBDD/Linopt.hs
+++ b/src/OBDD/Linopt.hs
@@ -1,40 +1,39 @@
-module OBDD.Linopt where
+module OBDD.Linopt (linopt) where
 
-import OBDD (OBDD, fold)
+import OBDD (OBDD, full_fold)
 
 import qualified Data.Map.Strict as M
-
-type Item v w = (w, [(v,Bool)])
+import Data.Bool (bool)
 
 -- | solve the constrained linear optimisation problem:
 -- returns an assignment that is a model of the BDD
 -- and maximises the sum of weights of variables.
+-- The set of keys of the weight map *must* be the
+-- full set of variables.
 linopt :: ( Ord v , Num w, Ord w ) 
        => OBDD v 
        -> M.Map v w 
        -> Maybe (w, M.Map v Bool)
-linopt d m = ( \(w,kvs) -> (w,M.fromList kvs) ) <$>
-   fold ( \ leaf  -> if leaf then Just (0, []) else Nothing )
-       ( \ v ml mr -> case (ml,mr) of
-          (Nothing, Just r) -> Just $   add m v $ fill m v r
-          (Just l, Nothing) -> Just $ noadd m v $ fill m v l
+linopt d m = full_fold (M.keysSet m) 
+   ( bool Nothing ( Just (0, M.empty) ))
+   ( \ v ml mr -> case (ml,mr) of
+          (Just l, Nothing) -> Just $ noadd m v l
+          (Nothing, Just r) -> Just $   add m v r
           (Just l,  Just r) -> Just $
-                  let l' = noadd m v $ fill m v l 
-                      r' =   add m v $ fill m v r
+                  let l' = noadd m v l
+                      r' =   add m v r
                   in  if fst l' >= fst r' then l' else r'
+          -- the following *can* happen for
+          -- interpolated nodes directly above False:
+          (Nothing, Nothing) -> Nothing
        )
        d
 
-fill :: (Ord v, Num w) => M.Map v w -> v -> Item v w -> Item v w
-fill m v (w, xs) = 
-    let vs = (case xs of
-               [] -> id
-               (u,_):_ -> takeWhile (\(k,v) -> k > u) ) 
-           $ dropWhile (\(k,_) -> k >= v) 
-           $ M.toDescList m
-    in  foldr (add m) (w, xs) $ map fst vs
-
-noadd, add :: (Ord v, Num w) => M.Map v w -> v -> Item v w -> Item v w
-noadd m v (w,xs) = (w          , (v,False) : xs)
-add   m v (w,xs) = (w + m M.! v, (v, True) : xs)
+type Item v w = (w, M.Map v Bool)
 
+noadd, add :: (Ord v, Num w) 
+           => M.Map v w -> v -> Item v w -> Item v w
+noadd m v (w, b) = 
+  (w                          , M.insert v False b)
+add   m v (w, b) = 
+  (w + M.findWithDefault 0 v m, M.insert v True  b)
diff --git a/src/OBDD/Operation.hs b/src/OBDD/Operation.hs
--- a/src/OBDD/Operation.hs
+++ b/src/OBDD/Operation.hs
@@ -4,11 +4,12 @@
 module OBDD.Operation 
 
 ( (&&), (||), not, and, or
-, bool, implies, equiv, xor
+, ite, bool, implies, equiv, xor
 , unary, binary
 , instantiate
 , exists, exists_many
 , fold, foldM
+, full_fold, full_foldM
 )
 
 where
@@ -40,6 +41,9 @@
 
 bool :: Ord v => OBDD v -> OBDD v -> OBDD v -> OBDD v
 bool f t p = (f && not p) || (t && p)
+
+ite :: Ord v => OBDD v -> OBDD v -> OBDD v -> OBDD v
+ite i t e = bool e t i
 
 equiv :: Ord v => OBDD v -> OBDD v -> OBDD v
 equiv = symmetric_binary (==)