Untangling with Continued Fractions: part 5

> {-# LANGUAGE MultiParamTypeClasses,FunctionalDependencies,FlexibleInstances,GeneralizedNewtypeDeriving #-}

> module Main where

> import qualified Data.Map as M
> import Control.Monad
> import Prelude
> import Ratio hiding (numerator,denominator)

> infixl 5 .+
> infixl 6 .*

So far I've described a way to use monad notation to represent knots, links and tangles, and I've provided a way to implement a monad so that the notation gives rise to the bracket polynomial when applied to a knot. What happens when we do the same with a tangle?

A tangle corresponds to something with this type:

> type Tangle = (Bool,Bool) -> KnotMonad (Bool,Bool)

Here are two examples of tangles, ie. expressions with this type:

> infinity (i,j) = return (i,j)
> zero (i,j) = do
>   cup (i,j)
>   cap ()

These correspond to these tangles:

And the names of those expressions are the corresponding extended rational values.

Now here's a surprising but nontrivial fact that holds when T is a rational tangle: if we set A²=i (where i²=-1) then the value of the monadic expression corresponding to the tangle is of the form a<[∞]>+b<[0]> where the [r] are the rational tangles above, and <T> is the bracket polynomial (evaluated at A²=i). The rational associated to the original tangle is then given by -ia/b. I believe this discovery is, like many other things I've talked about in this blog, due to Conway.

We now have the makings of a complete algorithm. Using the polynomial monad defined in the previous installment we evaluate our tangle. We then compute the coefficients a and b, which are polynomials in A. Now compute a and b at A²=-i, and then compute -ia/b and use the operations z→z+1, z→z-1 and z→-1/z to reduce this rational to zero, interpreting these operations as twist, antitwist and rotation respectively (as I described here).

There's a shortcut I'm going to take but I don't have a proof it's guaranteed to work so it's risky. Instead of computing polynomials in A, doing polynomial division, and then evaluating at A²=-i, I'm going to do everything from the first step with the assumption that A²=-1. (The catch is that if a and b are polynomials in A, then a and b may have a common divisor that cancels when computing a/b. If we substitute A²=-1 from the beginning we get division of zero by zero. I don't know for sure if this can't happen. It's not a big deal if it does though, I just need to write a polynomial division routine.)

Now if A²=-1 then we could choose A=(1/2)^1/2+(1/2)^1/2i. That suggests using real arithmetic but we need to extract results that are rational numbers. So instead introduce the field K=Q((1/2)^1/2), ie. the field of numbers of the form u+v(1/2)^1/2. We can implement this as:

> data K = K { re2::Rational, im2::Rational } deriving (Eq,Show)

> instance Num K where
>   K a b + K a' b' = K (a+a') (b+b')
>   K a b * K a' b' = K (a*a'+2*b*b') (a*b'+a'*b)
>   negate (K a b) = K (negate a) (negate b)
>   abs _ = error ""
>   signum _ = error ""
>   fromInteger i = K (fromInteger i) 0

> instance Fractional K where
>   recip (K a b) = let r = recip (a*a-2*b) in K (r*a) (-r*b)
>   fromRational x = K x 0

We actually want to work with the field L=K(i), ie. complex numbers built from elements of K. We can implement this as:

> type L = Complex K

(Notice how it's not at all hard to build up algebraic structures like this in Haskell without a symbolic algebra package. For more powerful code in this vein, see David Amos's code).

Now we have our field, our monad is that of vector spaces over L:

> type KnotMonad = V L

Let's define some useful constants. The square root of 1/2, A, B=1/A and i:

> hr2 = K 0 (1%2)
> a = hr2 :+ hr2
> b = hr2 :+ (-hr2)
> i = 0 :+ 1

Here's the knot monad code I've presented before. (This version uses -|()|- in a way that may appear a bit useless. I'm just making explicit when we're working with a 1-dimensional vector space generated by -|()|-, which is much the same thing as the base field.)

> cup :: (Bool,Bool) -> KnotMonad ()
> cup (u,v) = case (u,v) of
>   (False,True) -> (-i * b) .* return ()
>   (True,False) -> (i * a) .* return ()
>   otherwise -> vzero

> cap :: () -> KnotMonad (Bool,Bool)
> cap () = (-i * b) .* return (False,True) .+ (i * a) .* return (True,False)

> over :: Tangle
> over (u,v) = a .* do
>   () <- cup (u,v)
>   cap ()
>   .+
>   b .* return (u,v)

> under :: Tangle
> under (u,v) = b .* do
>   () <- cup (u,v)
>   cap ()
>   .+
>   a .* return (u,v)

Here's where we work out the coefficients a and b. We treat a tangle as the function that it is and look at what happens when we apply it to two different values in (Bool,Bool). Applying to value' to infinity and zero shows that they do in fact compute a and b correctly:

> value' p =
>   let alpha = coefficient (False,False) $ p (False,False)
>       beta  = coefficient (True,False)  $ p (False,True)
>   in (alpha,beta)

From a and b we extract -ia/b as a rational:

> value :: Tangle -> Rational
> value p =
>   let (a,b) = value' p
>   in re2 $ realPart $ -i*a/b

And now we reduce to 0. This isn't the best algorithm for reducing rationals to zero. But it's what I wrote in my first attempt and I'm sticking with it for consistency with part one:

> steps :: Rational -> [String]
> steps 0 = []
> steps q | q<= -1 = "twist":steps (q+1)
>         | -1<q && q<1 = "rotate":steps (-1/q)
>         | q>=1 = "antitwist":steps (q-1)

> untangle t = steps (value t)

Here's the example from the first installment:

> example :: Tangle
> example (a,b) = do
>   (c,d) <- over (a,b)
>   (e,f) <- cap ()
>   (g,h) <- over (c,e)
>   (i,j) <- over (f,d)
>   (m,n) <- cap ()
>   (k,l) <- cap ()
>   (q,r) <- over (h,k)
>   (s,y) <- over (l,i)
>   (o,p) <- over (n,g)
>   (t,u) <- under (p,q)
>   (v,w) <- under (r,s)
>   (x,z) <- over (y,j)
>   cup (o,t)
>   cup (u,v)
>   cup (w,x)
>   return (m,z)

I'll leave you to type untangle example in ghci and refer back to my original video.

To find out more, and to read about the details I've left out, try the many papers by Kauffman, including this one.

---

That was hard! I didn't anticipate how much stuff I'd have to introduce to just sketch out how this works. It's one of those cases where the comments need to be about 30 or 40 times larger than the code and done properly, there's enough material to fill a good few chapters of a book. So I apologise if I skimped on some details along the way, and left out some needed explanations. But I hope that I've at least motivated you to take a peek at knot theory and appreciate how nicely tangles and monads play together.

And remember this is only code to illustrate some mathematics - it's inefficient and not how I'd really untangle big rational tangles! (For example, just using the monad associativity laws you can often rearrange the monad expressions for knots into equivalent forms that execute much quicker.)

---

The rest is similar to code I've used in previous installments:

> swap (x,y) = (y,x)

> class Num k => VectorSpace k v | v -> k where
>     vzero :: v
>     (.+) :: v -> v -> v
>     (.*) :: k -> v -> v
>     (.-) :: v -> v -> v
>     v1 .- v2 = v1 .+ ((-1).*v2)

> data V k a = V { unV :: [(k,a)] } deriving (Show)

> reduce x = filter ((/=0) . fst) $ fmap swap $ M.toList $ M.fromListWith (+) $ fmap swap $ x

> instance Num k => Functor (V k) where
>     fmap f (V as) = V $ map (\(k,a) -> (k,f a)) as

> instance Num k => Monad (V k) where
>     return a = V [(1,a)]
>     x >>= f = join (fmap f x)
>         where join x = V $ concat $ fmap (uncurry scale) $ unV $ fmap unV x
>               scale k1 as = map (\(k2,a) -> (k1*k2,a)) as

> instance Num r => MonadPlus (V r) where
>     mzero = V []
>     mplus (V x) (V y) = V (x++y)

> instance (Num k,Ord a) => VectorSpace k (V k a) where
>     vzero = V []
>     V x .+ V y = V (x ++ y)
>     (.*) k = (>>= (\a -> V [(k,a)]))

> coefficient b (V bs) = maybe 0 id (lookup b (map swap (reduce bs)))

> data Complex a = (:+) { realPart :: a, imagPart :: a } deriving (Eq,Show)

> instance Num a => Num (Complex a) where
>   (a :+ b) + (a' :+ b') = (a+a') :+ (b+b')
>   (a :+ b) * (a' :+ b') = (a*a'-b*b') :+ (a*b'+a'*b)
>   negate (a :+ b) = (-a) :+ (-b)
>   fromInteger n = fromInteger n :+ 0
>   abs (a :+ b) = undefined
>   signum (a :+ b) = undefined

> instance Fractional a => Fractional (Complex a) where
>   recip (a :+ b) = let r = recip (a*a+b*b) in ((a*r) :+ (-b*r))
>   fromRational q = fromRational q :+ 0