-- :l "Macintosh HD:Users:ms:MSS:Vorlesungen:FP-1-2003:Hugsprogs:gauss.hs"
-- :l "MSS/Vorlesungen/FP-1-2003/Hugsprogs/gauss.hs"

--  Lšst linare Gleichungssysteme mittels Gaussverfahren
--  Koeffizientenmatrix und Ergebnis in einer Matrix,
--  Resultat ist die Lšsung
--  Keine Ausgabe eines Lšsungsraumes, falls es das gibt: fehlt
--
--Datenstruktur: [([a1 a2 a3 b] -> (x1 x2 x3 x4)

--   neue Module, deshalb Aenderungen

import Ratio 


vectoradd [] []     = []
vectoradd (x:xs) (y:ys) = (x+y): vectoradd xs ys

vectorscmal a [] = []
vectorscmal a (x:xs) = a*x : vectorscmal a xs 

scalarprod [] [] = fromInteger 0
scalarprod (x:xs) (y:ys) = x*y + scalarprod xs ys

nullnachvorne n xs  =  
     let selplus = filter (\s -> s !! (n-1) /= fromInteger 0) xs
         selnull = filter (\s -> s !! (n-1) == fromInteger 0) xs
     in selplus ++ selnull      

gauss:: (Fractional a) => [[a]] -> [a]
gauss eingabe =
   let order = (length (head eingabe))-1
       dreieck =  gaussr order 1 eingabe
       einok = eingabetest order eingabe
       loesung = gaussloesung order dreieck  
    in if einok then loesung
       else error "Eingabeformat falsch"


gaussschritt n matrix = 
   let matrixa = take (n-1) matrix 
       matrixb = drop (n-1) matrix
       matrixbnn = nullnachvorne n matrixb
       matrixnn = matrixa ++ matrixbnn 
       matrixrest = drop (n-1) matrixnn
       mhead = head matrixrest
       mtail = tail matrixrest
       pivot = mhead !! (n-1)
       restneu = map (\zeile -> vectoradd (vectorscmal (-pivot) zeile)
                                          (vectorscmal (zeile !! (n-1)) mhead))
                     mtail
       matrixneu = (take n matrixnn) ++ restneu  
   in if pivot == fromInteger 0 then error "Abhaengige Loesungen: denn Pivot = 0"
      else matrixneu


eingabetest order eingabe =
    let nlen = map length eingabe
        rest = filter (\x -> x /= (order+1)) nlen
    in if rest /= [] then error "falsche Eingabematrix"
       else True
 
gaussschritttest n matrix = 
   let matrixa = take (n-1) matrix 
       matrixb = drop (n-1) matrix
       matrixbnn = nullnachvorne n matrixb
       matrixnn = matrixa ++ matrixbnn 
       matrixrest = drop (n-1) matrixnn
       mhead = head matrixrest
       mtail = tail matrixrest
       pivot = mhead !! (n-1)
       restneu = map (\zeile -> vectoradd (vectorscmal (-pivot) zeile)
                                          (vectorscmal (zeile !! (n-1)) mhead))
                     mtail
       matrixneu = (take n matrixnn) ++ restneu  
   in (matrixrest,mhead,mtail,pivot,restneu,matrixneu)


gaussdelta eingabe =
   let order = (length (head eingabe))-1
       dreieck =  gaussr order 1 eingabe
   in dreieck

gaussloesung max dreieck = 
    let rev = reverse dreieck
        in loesungr max rev []

loesungr  n [] alteloesung = alteloesung
loesungr  n dreieck alteloesung =
     let hddreieck = head dreieck
         tldreieck = tail dreieck
         tail1 = drop (n-1) hddreieck
         hdtail1 = head tail1
         tltail1 = tail tail1
         (restvector,b) = splitlast tltail1
         neueloesung = (b - (scalarprod alteloesung restvector)) / hdtail1
  in  loesungr  (n-1) tldreieck (neueloesung:alteloesung)

testloesungr  n dreieck alteloesung =
     let hddreieck = head dreieck
         tldreieck = tail dreieck
         tail1 = drop (n-1) hddreieck
         hdtail1 = head tail1
         tltail1 = tail tail1
         (restvector,b) = splitlast tltail1
         neueloesung = (b - (scalarprod alteloesung restvector)) / hdtail1
  in  (hddreieck,tldreieck,tail1,hdtail1,tltail1,restvector,b,neueloesung,
       loesungr  (n-1) tldreieck (neueloesung:alteloesung))
          
gaussr max i eingabe = 
    if i == max 
    then eingabe
    else gaussr max (i+1) (gaussschritt  i eingabe)

splitlast [x] = ([],x)
splitlast (x:xs) = 
       let (ys,y) = splitlast xs
       in (x:ys,y)

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

--test1 = gauss [[1,1,6],[1,-1,10]]
-- test1:: (Num a, Fractional a) => [a]
test1 = gauss test1mat
test1mat:: Num a => [[a]]
test1mat  =  [[1,1,6],[1,-1,10]]
test11 = gaussdelta [[1,1,6],[1,-1,10]]
test12 = gaussr 2 1 [[1,1,6],[1,-1,10]]
test13 = gaussschritt 1 [[1,1,6],[1,-1,10]]
testchinesenmatrix:: Fractional a => [[a]]
testchinesenmatrix = [[5,4,3,2,1496],[4,2,6,3,1175],[3,1,7,5,958],[2,3,5,1,861]]
testchinesen :: (Fractional a, Num a) => [a]
testchinesen = gauss testchinesenmatrix
test4matrix = [[1,1,1,0,205],[1,1,0,1,169],[1,0,1,1,158],[0,1,1,1,194]]::[[Ratio Integer]]
--test4matrix = [[1.0,1.0,1.0,0.0,205.0],[1.0,1.0,0.0,1.0,169.0],[1.0,0.0,1.0,1.0,158.0],[0.0,1.0,1.0,1.0,194.0]]

test4 = gauss test4matrix
test41 = gaussschritttest 1 test4matrix
test42 = gaussschritttest 2 [[1, 1, 1, 0, 205], [0, 0, 1, -1, 36], [0, 1, 0, -1, 47], [0, -1, -1, -1, -194]]
test43 = nullnachvorne 2 [[1, 1, 1, 0, 205], [0, 0, 1, -1, 36], [0, 1, 0, -1, 47], [0, -1, -1, -1, -194]] 

test5 = gauss [[1,1,1,1],[1,1,1,1],[1,1,1,1],[1,1,1,1,1]]
test6 = gauss [[1.5, 2.5]]

-- matrix7 =  [[2.0,4.0,3.0,1.0,79.0],[0.0,2.0,1.0,3.0,41.0],[6.0,0.0,5.0,2.0,88.0],[8.0,4.0,3.0,1.0,127.0]]
matrix7 =  [[2,4,3,1,79],[0,2,1,3,41],[6,0,5,2,88],[8,4,3,1,127]]::[[Ratio Integer]]
test7 = gauss matrix7
test71 = gaussr 4 1 matrix7
sc7 = scalarprod [8,10,6,5]
testmasch = gauss [[8.0,8.0,8.0,8.0,1034.0],[3.0,-5.0,0.0,0.0,0.0],[0.0,7.0,-4.0,0.0,0.0],[0.0,2.0,0.0,-9.0,0.0]] 
testmaschrat = gauss [[8,8,8,8,1034],[3,-5,0,0,0],[0,7,-4,0,0],[0,2,0,-9,0]] 
matrixmaschrat  :: Fractional a => [[a]]
matrixmaschrattest =  [[8,8,8,8,1034],[3,-5,0,0,0],[0,7,-4,0,0],[0,2,0,-9,0]] 
matrixmaschrat =  [[8,8,8,8,1034],[3,-5,0,0,0],[0,7,-4,0,0],[0,2,0,-9,0]]:: Fractional a => [[a]]
testmaschratio:: Fractional a => [a]
testmaschratio = gauss matrixmaschrat

test8mat = [[1,1,1],[2,2,2]]::Fractional a => [[a]]

test8 = gauss test8mat