-- This program solves the first laboratory exercise on our introductory
-- Java course: see 
-- http://www.cs.chalmers.se/ComputingScience/Education/Courses/d1pt/d1ptb/lab1/index.html

import List
import Random
import GraphicsUtils

-- We represent vectors as lists, rather than as a new type. This is more
-- general: we can represent vectors of any dimension. It also simplifies
-- the programming. It is slightly less efficient for 2D vectors, but who
-- cares?

type Vector = [Double]

vector2 x y = [x,y]
getX [x,y] = x
getY [x,y] = y

v1 .* v2 = sum (zipWith (*) v1 v2)		-- scalar product
v1 .+ v2 = zipWith (+) v1 v2			-- vector addition

type Matrix = [Vector]

matrix2 a b c d = [vector2 a b, vector2 c d]
rows m = m
cols m = transpose m

det [] = 1
det rs = foldr (-) 0 [a * det m | (a:_, m) <- zip rs (minors (map tail rs))]

minors (x:xs) = xs:map (x:) (minors xs)
minors [] = []

m #* v = map (.* v) m				-- matrix-vector multiplication

-- Define a symbolic expression type for testing matrix operations.
-- This isn't part of the lab exercise, but it's fun.
-- Example: 
-- Main> det (matrix2 (var"a") (var"b") (var"c") (var"d"))
-- a*d - b*c  

data Polynomial = Poly [(Int,[String])] deriving Eq

instance Show Polynomial where
  show (Poly []) = "0"
  show (Poly terms) = foldr1 join (map showTerm terms)
    where join s1 ('-':s2) = s1 ++ " - " ++ 
				if take 2 s2=="1*" then drop 2 s2 else s2
	  join s1 s2 = s1 ++ " + " ++ s2
	  showTerm (k,[]) = show k
	  showTerm (1,vars) = showVars vars
	  showTerm (k,vars) = show k ++ "*" ++ showVars vars
	  showVars = foldr1 (\v v' -> v++"*"++v')

instance Num Polynomial where
  fromInt 0 = Poly []
  fromInt n = Poly [(n,[])]
  negate (Poly terms) = Poly [(-k,vs) | (k,vs) <- terms]
  Poly ts + Poly ts' = Poly (merge ts ts')
    where merge [] ts' = ts'
	  merge ts [] = ts
	  merge ((k,vs):ts) ((k',vs'):ts')
	    | vs < vs' = (k,vs):merge ts ((k',vs'):ts')
	    | vs == vs' = if k+k'==0 then merge ts ts' 
			  else (k+k',vs):merge ts ts'
	    | vs > vs' = (k',vs'):merge ((k,vs):ts) ts'
  Poly ts * Poly ts' = Poly (sortBy (\(_,vs) (_,vs') -> compare vs vs')
				[(k*k',sort(vs++vs')) 
				| (k,vs) <- ts, (k',vs') <- ts'])

var s = Poly [(1,[s])]

-- Now for the picture: we define transformations corresponding to the
-- four parts of the image.

blue, red, green, yellow :: (Matrix, Vector)
blue = (matrix2 0.2 (-0.25) 0.25 0.2, vector2 0 64)
red = (matrix2 0.64 0.04 (-0.04) 0.84, vector2 0 64)
yellow = (matrix2 (-0.19) 0.29 0.29 0.19, vector2 0 18)
green = (matrix2 0 0 0 0.17, vector2 0 0)

chooseTransformation p =
  alts !! length (takeWhile (<p) (scanl1 (+) probs))
  where alts = [blue, red, yellow, green]
	probs= [0.1,  0.75,0.1,    0.05]

--transform :: (Matrix,Vector) -> Vector -> Vector
transform (m,b) v = m #* v .+ b

-- convert a random integer to a probability.

probability n = fromInt (n `mod` 1000) / 999

-- the list of points to fill (avoid storing the entire list)

points () = generate (map (chooseTransformation.probability) 
		          (randoms (StdGen 0 1)))
		     (vector2 0 0)
  where generate (t:ts) v = v : generate ts (transform t v)

-- it only remains to draw the points.

pointGraphic v = withColor Green$ polygon [(x,y),(x+1,y),(x+1,y+1),(x,y+1)]
  where (x,y) = (round (getX v) + 200, round (getY v))

main = runGraphics$ do
  w <- openWindow "Fractal" (400,400)
  sequence_ [drawInWindow w (pointGraphic p) | p <- take 10000 (points ())]
  getKey w
  closeWindow w