--<PRE>

type Quad a = (a,a,a,a)

data Quadtree a = Z a | S (Quad (Quadtree a))

mapQ :: (a -> b) -> Quadtree a -> Quadtree b
mapQ f (Z x) = Z(f x)
mapQ f (S(w,x,y,z)) = S(mapQ f w,mapQ f x,mapQ f y,mapQ f z)

--- Now consider

data QT2 a = Z2 a
           | S2 (QT2(Quad a))

-- S2 :: QT2 (Quad a) -> QT2 a
-- S2 :: QT2 (a,a,a,a) -> QT2 a

-- Note how things of type Quadtree can have any shape
ex1 = S(Z 2,Z 5, Z 0,S(Z 1,Z 4,Z 7,Z 1))


-- Including square things
ex2 = S(S(Z  1, Z  2, Z  3, Z  4)
       ,S(Z  5, Z  6, Z  7, Z  8)
       ,S(Z  9, Z 10, Z 11, Z 12)
       ,S(Z 13, Z 14, Z 15, Z 16))

-- But things of type QT2 can only be square       
ex3 = S2(S2(Z2(( 1, 2, 3, 4),
               ( 5, 6, 7, 8),
               ( 9,10,11,12),
               (13,14,15,16))))

ex2a = mapQ (+1) ex2

-- can we write a mapQ2 ?
-- only if we give an explicit type signature

mapQ2 :: (a -> b) -> QT2 a -> QT2 b
mapQ2 f (Z2 x) = Z2(f x)
mapQ2 f (S2 z) = S2(mapQ2 (mapQuad f) z)
  where mapQuad f (a,b,c,d) = (f a, f b, f c, f d)

-- How could you make a 3 by 3 square matrix?

type Tri a = (a,a,a)

-- Tri Trees
data TT a = Z3 a
          | S3 (TT(Tri a))
          
ex4 = S3 (S3 (Z3 ((1,2,3),
                  (4,5,6),
                  (7,8,9))))
  
-- Now what about square matrices for sizes that are not powers of 2?  

-- Fast exponentiation

power b n = fast 1 b n

fast acc b n 
  | n==0   = acc
  | even n = fast acc (b*b) (half n)
  | odd n  = fast (acc*b) (b*b) (half n)

half :: Int -> Int
half x = div x 2

-- lifted to the type level we get vectors of arbitrary size

type Vector a = Vector_ () a
data Vector_ v w 
  = Zero v
  | Even (Vector_ v (w,w))
  | Odd  (Vector_ (v,w) (w,w))
  
create x n = create_ () x n
create_ :: v -> w -> Int -> Vector_ v w
create_ v w n 
  | n==0   = Zero v
  | even n = Even (create_ v (w,w) (half n))
  | odd n  = Odd  (create_ (v,w) (w,w) (half n))

-----------------------------------------------------
-- To create rectangular arrays, we repeat the definition
-- and place vector at the leaf (Zero) node.

type Rect a = Rect_ () a
data Rect_ v w 
  = ZeroR (Vector v)
  | EvenR (Rect_ v (w,w))
  | OddR  (Rect_ (v,w) (w,w))
  
 

---------------------------------------------

newtype Empty a = E ()          -- the empty vector
newtype Id a = I a              -- the vector of size 1
newtype Pair v w a = P(v a,w a) -- vectors of size v + w

type Square a = Square_ Empty Id a
data Square_ v w a
  = ZeroS (v(v a))
  | EvenS (Square_ v (Pair w w) a)
  | OddS  (Square_ (Pair v w) (Pair w w) a)
  
  
f x = (EvenS (OddS (ZeroS x)))
g x y = P(E(),P(I x,I y))

row1 = P (E (),P(I 3,I 4))
row2 = P (E (),P(I 5,I 6))
table = P(E (),P(I row1, I row2))

z = EvenS(OddS (ZeroS table))

-- 4 by 4 array
h w x y z = P(E(),P(P(I w,I x),P(I y,I z)))
r1 = h 1 2 3 4
r2 = h 5 6 7 8
r3 = h 9 10 11 12
r4 = h 13 14 15 16
tab = h r1 r2 r3 r4
dd = EvenS(EvenS(OddS(ZeroS tab)))

------------------------------------
-- indexing functions

subE i (E ()) = error "no index in empty vector"

subI 0 (I x) = x
subI n (I x) = error "only 0 can index vector of size 1"

subP subv subw vsize i (P (v,w))
  | i < vsize = subv i v
  | i >= vsize = subw (i-vsize) w
  
sub (i,j) m = sub_ subE subI 0 1 (i,j) m


sub_ :: (forall b. Int -> v b -> b) ->
        (forall b. Int -> w b -> b) ->
        Int -> Int -> (Int,Int) -> Square_ v w a -> a
         
sub_ subv subw vsize wsize (i,j) x =
  case x of
    ZeroS vv -> subv i (subv j vv)
    EvenS m -> sub_ subv (subP subw subw wsize)
                    vsize (wsize+wsize) (i,j) m
    OddS m -> sub_ (subP subv subw vsize) (subP subw subw wsize)
                   (vsize+wsize) (wsize+wsize) (i,j) m

  
--</PRE>