diff --git a/Eigensystem.hs b/Eigensystem.hs --- a/Eigensystem.hs +++ b/Eigensystem.hs @@ -1,11 +1,10 @@ - ------------------------------------------------------------------------------ --- Haskell module: Eigensystem --- Date: initialized 2001-03-25, last modified 2001-03-25 --- Author: Jan Skibinski, Numeric Quest Inc. --- Location: http://www.numeric-quest.com/haskell/Eigensystem.hs --- See also: http://www.numeric-quest.com/haskell/QuantumVector.html --- See also: http://www.numeric-quest.com/haskell/Orthogonals.html +-- Haskell module: Eigensystem +-- Date: initialized 2001-03-25, last modified 2001-03-25 +-- Author: Jan Skibinski, Numeric Quest Inc. +-- Location: http://www.numeric-quest.com/haskell/Eigensystem.hs +-- See also: http://www.numeric-quest.com/haskell/QuantumVector.html +-- See also: http://www.numeric-quest.com/haskell/Orthogonals.html -- -- Description: -- @@ -30,7 +29,7 @@ -- manipulation on operators is too difficult, while it is trivial -- to obtain the corresponding results via matrices. One example is a -- computation of a Hermitian conjugate of A: --- < ei | A' | ej > = conjugate < ej | A | ej > +-- < ei | A' | ej > = conjugate < ej | A | ej > -- (Here ' stands for a dagger) -- If however the operator A is made from a product or a sum of simpler -- operators, whose Hermitian conjugates are known to us, then the @@ -54,57 +53,57 @@ import Complex import QuantumVector import LinearAlgorithms (triangular, tridiagonal, triangular2) -import List (findIndex) +import List (findIndex) ---------------------------------------------------------------------------- -- Category: Eigensystem for QuantumVector ---------------------------------------------------------------------------- -eigenvalues :: Ord a => Bool -> Int -> [Ket a] -> (Ket a -> Ket a) -> [Scalar] +eigenvalues :: Ord a => Bool -> Int -> [Ket a] -> (Ket a -> Ket a) -> [Scalar] eigenvalues doTri n es a - -- A list of eigenvalues of operator 'a' - -- obtained after 'n' triangularizations - -- of a matrix corresponding to operator 'a' - -- where - -- 'es' is a list of base vectors - -- 'doTri' declares whether or not we - -- want the initial tridiagonalization - -- (applies to Hermitian operators only) - | doTri == True = f b1 - | otherwise = f b + -- A list of eigenvalues of operator 'a' + -- obtained after 'n' triangularizations + -- of a matrix corresponding to operator 'a' + -- where + -- 'es' is a list of base vectors + -- 'doTri' declares whether or not we + -- want the initial tridiagonalization + -- (applies to Hermitian operators only) + | doTri == True = f b1 + | otherwise = f b where - f c = diagonals $ operator es $ triangular n c - diagonals us = [toBra e <> us e | e <- es] - b = matrix es a - b1 = tridiagonal b - - + f c = diagonals $ operator es $ triangular n c + diagonals us = [toBra e <> us e | e <- es] + b = matrix es a + b1 = tridiagonal b + + eigenpairs :: Ord a => Int -> [Ket a] -> (Ket a -> Ket a) -> ([Scalar], [Ket a]) eigenpairs n es a - -- A pair of lists (eigenvalues, eigenvectors) of hermitian - -- operator 'a' obtained after 'n' triangularizations of 'a' - -- where - -- 'es' is a list of base vectors - -- Note: For a moment this applies only to Hermitian operators - -- until we decide what would be the best way to compute eigenvectors - -- of a triangular matrix: the method from module Orthogonal, power - -- iteration, etc. + -- A pair of lists (eigenvalues, eigenvectors) of hermitian + -- operator 'a' obtained after 'n' triangularizations of 'a' + -- where + -- 'es' is a list of base vectors + -- Note: For a moment this applies only to Hermitian operators + -- until we decide what would be the best way to compute eigenvectors + -- of a triangular matrix: the method from module Orthogonal, power + -- iteration, etc. = (ls, xs) where - (t, q) = triangular2 n b - b = matrix es a - ls = [ tk!!k | (tk, k) <- zip t [0..length t - 1] ] - xs = [compose qk es | qk <- q] + (t, q) = triangular2 n b + b = matrix es a + ls = [ tk!!k | (tk, k) <- zip t [0..length t - 1] ] + xs = [compose qk es | qk <- q] adjoint :: Ord a => [Ket a] -> (Ket a -> Ket a) -> (Ket a -> Ket a) adjoint es a - -- A Hermitian conjugate of operator a, - -- (or a-dagger, or adjoint to a) - -- where 'es' is a list of base vectors - = operator es ms + -- A Hermitian conjugate of operator a, + -- (or a-dagger, or adjoint to a) + -- where 'es' is a list of base vectors + = operator es ms where - ms = [[ conjugate (toBra ei <> vj) | vj <- v] | ei <- es] - v = [a ej | ej <- es] + ms = [[ conjugate (toBra ei <> vj) | vj <- v] | ei <- es] + v = [a ej | ej <- es] ---------------------------------------------------------------------------- @@ -112,28 +111,28 @@ ---------------------------------------------------------------------------- operator :: Ord a => [Ket a] -> [[Scalar]] -> Ket a -> Ket a -operator basis ms x - -- Definition of an operator corresponding - -- to a matrix 'ms' given as a list of scalar - -- columns - -- where - -- 'basis' is a complete list of base vectors - -- 'x' is any ket vector from this space - = a >< x +operator bss ms x + -- Definition of an operator corresponding + -- to a matrix 'ms' given as a list of scalar + -- columns + -- where + -- 'bss' (basis) is a complete list of base vectors + -- 'x' is any ket vector from this space + = a >< x where - a u = case (findIndex (u == ) basis) of - Just k -> compose (ms !! k) basis - Nothing -> error "Out of bounds" + a u = case (findIndex (u == ) bss) of + Just k -> compose (ms !! k) bss + Nothing -> error "Out of bounds" matrix :: Ord a => [Ket a] -> (Ket a -> Ket a) -> [[Scalar]] -matrix basis a - -- List of scalar columns representing - -- the operator 'a' in a given 'basis' +matrix bss a + -- List of scalar columns representing + -- the operator 'a' in a given 'basis' = [[ei' <> vj | ei' <- e'] | vj <- v] where - v = [a ej | ej <- basis] - e' = [toBra ei | ei <- basis] + v = [a ej | ej <- bss] + e' = [toBra ei | ei <- bss] ---------------------------------------------------------------------------- -- Category: Test data @@ -142,32 +141,34 @@ matrixA :: [[Scalar]] matrixA - -- Test matrix A represented as list of scalar columns. - = [ - [1, 2, 4, 1, 5] - , [2, 3, 2, 6, 4] - , [4, 2, 5, 2, 3] - , [1, 6, 2, 7, 2] - , [5, 4, 3, 2, 9] - ] + -- Test matrix A represented as list of scalar columns. + = [ + [1, 2, 4, 1, 5] + , [2, 3, 2, 6, 4] + , [4, 2, 5, 2, 3] + , [1, 6, 2, 7, 2] + , [5, 4, 3, 2, 9] + ] -opA = operator basisA matrixA +opA :: Ket Int -> Ket Int +opA = operator basisA matrixA -basisA = map Ket [1..5::Int] -- or: map Ket "abcde", etc. - +basisA :: [Ket Int] +basisA = map Ket [1..5::Int] -- or: map Ket "abcde", etc. + --------------------------------------------------------------------------- -- Copyright: -- --- (C) 2001 Numeric Quest, All rights reserved +-- (C) 2001 Numeric Quest, All rights reserved -- -- Email: jans@numeric-quest.com -- --- http://www.numeric-quest.com +-- http://www.numeric-quest.com -- -- License: -- --- GNU General Public License, GPL +-- GNU General Public License, GPL -- --------------------------------------------------------------------------- - + diff --git a/Fraction.hs b/Fraction.hs --- a/Fraction.hs +++ b/Fraction.hs @@ -1,69 +1,69 @@ -- Module: -- --- Fraction.hs +-- Fraction.hs -- -- Language: -- --- Haskell +-- Haskell -- -- Description: Rational with transcendental functionalities -- -- --- This is a generalized Rational in disguise. Rational, as a type --- synonim, could not be directly made an instance of any new class --- at all. --- But we would like it to be an instance of Transcendental, where --- trigonometry, hyperbolics, logarithms, etc. are defined. --- So here we are tiptoe-ing around, re-defining everything from --- scratch, before designing the transcendental functions -- which --- is the main motivation for this module. +-- This is a generalized Rational in disguise. Rational, as a type +-- synonim, could not be directly made an instance of any new class +-- at all. +-- But we would like it to be an instance of Transcendental, where +-- trigonometry, hyperbolics, logarithms, etc. are defined. +-- So here we are tiptoe-ing around, re-defining everything from +-- scratch, before designing the transcendental functions -- which +-- is the main motivation for this module. -- --- Aside from its ability to compute transcendentals, Fraction --- allows for denominators zero. Unlike Rational, Fraction does --- not produce run-time errors for zero denominators, but use such --- entities as indicators of invalid results -- plus or minus --- infinities. Operations on fractions never fail in principle. +-- Aside from its ability to compute transcendentals, Fraction +-- allows for denominators zero. Unlike Rational, Fraction does +-- not produce run-time errors for zero denominators, but use such +-- entities as indicators of invalid results -- plus or minus +-- infinities. Operations on fractions never fail in principle. -- -- However, some function may compute slowly when both numerators --- and denominators of their arguments are chosen to be huge. --- For example, periodicity relations are utilized with large --- arguments in trigonometric functions to reduce the arguments --- to smaller values and thus improve on the convergence --- of continued fractions. Yet, if pi number is chosen to --- be extremely accurate then the reduced argument would --- become a fraction with huge numerator and denominator --- -- thus slowing down the entire computation of a trigonometric --- function. +-- and denominators of their arguments are chosen to be huge. +-- For example, periodicity relations are utilized with large +-- arguments in trigonometric functions to reduce the arguments +-- to smaller values and thus improve on the convergence +-- of continued fractions. Yet, if pi number is chosen to +-- be extremely accurate then the reduced argument would +-- become a fraction with huge numerator and denominator +-- -- thus slowing down the entire computation of a trigonometric +-- function. -- -- Usage: -- --- When computation speed is not an issue and accuracy is important --- this module replaces some of the functionalities typically handled --- by the floating point numbers: trigonometry, hyperbolics, roots --- and some special functions. All computations, including definitions --- of the basic constants pi and e, can be carried with any desired --- accuracy. One suggested usage is for mathematical servers, where --- safety might be more important than speed. See also the module --- Numerus, which supports mixed arithmetic between Integer, --- Fraction and Cofra (Complex fraction), and returns complex --- legal answers in some cases where Fraction would produce --- infinities: log (-5), sqrt (-1), etc. +-- When computation speed is not an issue and accuracy is important +-- this module replaces some of the functionalities typically handled +-- by the floating point numbers: trigonometry, hyperbolics, roots +-- and some special functions. All computations, including definitions +-- of the basic constants pi and e, can be carried with any desired +-- accuracy. One suggested usage is for mathematical servers, where +-- safety might be more important than speed. See also the module +-- Numerus, which supports mixed arithmetic between Integer, +-- Fraction and Cofra (Complex fraction), and returns complex +-- legal answers in some cases where Fraction would produce +-- infinities: log (-5), sqrt (-1), etc. -- --- +-- -- Required: -- --- Haskell Prelude --- +-- Haskell Prelude +-- -- Author: -- --- Jan Skibinski, Numeric Quest Inc. +-- Jan Skibinski, Numeric Quest Inc. -- -- Date: -- --- 1998.08.16, last modified 2000.05.31 --- +-- 1998.08.16, last modified 2000.05.31 +-- -- See also bottom of the page for description of the format used --- for continued fractions, references, etc. +-- for continued fractions, references, etc. ------------------------------------------------------------------- module Fraction where @@ -73,57 +73,57 @@ infix 7 :-: ------------------------------------------------------------------- --- Category: Basics +-- Category: Basics ------------------------------------------------------------------- - + data Fraction = Integer :-: Integer - deriving (Eq) + deriving (Eq) num, den :: Fraction -> Integer -num (x:-:y) = x -den (x:-:y) = y - -reduce :: Fraction -> Fraction +num (x:-:_) = x +den (_:-:y) = y + +reduce :: Fraction -> Fraction reduce (x:-:0) - | x < 0 = (-1):-:0 - | otherwise = 1:-:0 + | x < 0 = (-1):-:0 + | otherwise = 1:-:0 reduce (x:-:y) = - (u `quot` d) :-: (v `quot` d) + (u `quot` d) :-: (v `quot` d) where d = gcd u v (u,v) | y < 0 = (-x,-y) - | otherwise = (x,y) - + | otherwise = (x,y) + (//) :: Integer -> Integer -> Fraction x // y = reduce (x:-:y) approx :: Fraction -> Fraction -> Fraction -approx eps (x:-:0) = x//0 +approx _ (x:-:0) = x//0 approx eps x = simplest (x-eps) (x+eps) - where - simplest x y - | y < x = simplest y x - | x == y = x - | x > 0 = simplest' (num x) (den x) (num y) (den y) - | y < 0 = - simplest' (-(num y)) (den y) (-(num x)) (den x) + where + simplest y z + | z < y = simplest z y + | y == z = y + | y > 0 = simplest' (num y) (den y) (num z) (den z) + | z < 0 = - simplest' (-(num z)) (den z) (-(num y)) (den y) | otherwise = 0 :-: 1 simplest' n d n' d' -- assumes 0 < n//d < n'//d' | r == 0 = q :-: 1 | q /= q' = (q+1) :-: 1 | otherwise = (q*n''+d'') :-: n'' - where + where (q,r) = quotRem n d (q',r') = quotRem n' d' (n'':-:d'') = simplest' d' r' d r ------------------------------------------------------------------- --- Category: Instantiation of some Prelude classes +-- Category: Instantiation of some Prelude classes ------------------------------------------------------------------- - + instance Read Fraction where - readsPrec p = + readsPrec p = readParen (p > 7) (\r -> [(x//y,u) | (x,s) <- reads r, ("//",t) <- lex s, (y,u) <- reads t ]) @@ -142,13 +142,13 @@ (x:-:y) * (x':-:y') = reduce ((x*x') :-: (y*y')) negate (x:-:y) = negate x :-: y abs (x:-:y) = abs x :-: y - signum (x:-:y) = signum x :-: 1 + signum (x:-:_) = signum x :-: 1 fromInteger n = fromInteger n :-: 1 - + instance Fractional Fraction where (x:-:0) / (x':-:0) = ((signum x * signum x'):-:0) - (x:-:y) / (x':-:0) = (0:-:1) - (x:-:0) / (x':-:y') = (x:-:0) + (_:-:_) / (_:-:0) = (0:-:1) + (x:-:0) / (_:-:_) = (x:-:0) (x:-:y) / (x':-:y') = reduce ((x*y') :-: (y*x')) recip (x:-:y) = if x < 0 then (-y) :-: (-x) else y :-: x fromRational a = x :-: y @@ -156,10 +156,10 @@ x = numerator a y = denominator a - + instance Real Fraction where - toRational (x :-: 0) = toRational (0%1) - -- or shoud we return some huge number instead? + toRational (_ :-: 0) = toRational ((0::Int)%(1::Int)) + -- or shoud we return some huge number instead? toRational (x :-: y) = toRational (x % y) @@ -189,426 +189,426 @@ numericEnumFromTo n m = takeWhile (<= m) (numericEnumFrom n) numericEnumFromThenTo n n' m = takeWhile p (numericEnumFromThen n n') where p | n' >= n = (<= m) - | otherwise = (>= m) + | otherwise = (>= m) ------------------------------------------------------------------ --- Category: Conversion --- from continued fraction to fraction and vice versa, --- from Taylor series to continued fraction. +-- Category: Conversion +-- from continued fraction to fraction and vice versa, +-- from Taylor series to continued fraction. ------------------------------------------------------------------- -type CF = [(Fraction, Fraction)] +type CF = [(Fraction, Fraction)] fromCF :: CF -> Fraction fromCF x = - -- - -- Convert finite continued fraction to fraction - -- evaluating from right to left. This is used - -- mainly for testing in conjunction with "toCF". - -- - foldr g (1//1) x - where - g :: (Fraction, Fraction) -> Fraction -> Fraction - g u v = (fst u) + (snd u)/v + -- + -- Convert finite continued fraction to fraction + -- evaluating from right to left. This is used + -- mainly for testing in conjunction with "toCF". + -- + foldr g (1//1) x + where + g :: (Fraction, Fraction) -> Fraction -> Fraction + g u v = (fst u) + (snd u)/v -toCF :: Fraction -> CF +toCF :: Fraction -> CF toCF (u:-:0) = [(u//0,0//1)] toCF x = - -- - -- Convert fraction to finite continued fraction - -- - toCF' x [] - where - toCF' u lst = + -- + -- Convert fraction to finite continued fraction + -- + toCF' x [] + where + toCF' u lst = case r of 0 -> reverse (((q//1),(0//1)):lst) - _ -> toCF' (b//r) (((q//1),(1//1)):lst) - where - a = num u - b = den u - (q,r) = quotRem a b + _ -> toCF' (b//r) (((q//1),(1//1)):lst) + where + a = num u + b = den u + (q,r) = quotRem a b approxCF :: Fraction -> CF -> Fraction -approxCF eps [] = 0//1 -approxCF eps x - -- - -- Approximate infinite continued fraction x by fraction, - -- evaluating from left to right, and stopping when - -- accuracy eps is achieved, or when a partial numerator - -- is zero -- as it indicates the end of CF. - -- - -- This recursive function relates continued fraction - -- to rational approximation. - -- - | den h == 0 = h - | otherwise = approxCF' eps x 0 1 1 q' p' 1 - where - h = fst (x!!0) - (q', p') = x!!0 - approxCF' eps x v2 v1 u2 u1 a' n - | abs (1 - f1/f) < eps = approx eps f - | a == 0 = approx eps f - | otherwise = approxCF' eps x v1 v u1 u a (n+1) - where - (b, a) = x!!n - u = b*u1 + a'*u2 - v = b*v1 + a'*v2 - f = u/v - f1 = u1/v1 - +approxCF _ [] = 0//1 +approxCF eps x + -- + -- Approximate infinite continued fraction x by fraction, + -- evaluating from left to right, and stopping when + -- accuracy eps is achieved, or when a partial numerator + -- is zero -- as it indicates the end of CF. + -- + -- This recursive function relates continued fraction + -- to rational approximation. + -- + | den h == 0 = h + | otherwise = approxCF' eps x 0 1 1 q' p' 1 + where + h = fst (x!!0) + (q', p') = x!!0 + approxCF' ept y v2 v1 u2 u1 a' n + | abs (1 - f1/f) < ept = approx ept f + | a == 0 = approx ept f + | otherwise = approxCF' ept y v1 v u1 u a (n+1) + where + (b, a) = y!!n + u = b*u1 + a'*u2 + v = b*v1 + a'*v2 + f = u/v + f1 = u1/v1 +fromTaylorToCF :: (Fractional a) => [a] -> a -> [(a, a)] fromTaylorToCF s x = - -- - -- Convert infinite number of terms of Taylor expansion of - -- a function f(x) to an infinite continued fraction, - -- where s = [s0,s1,s2,s3....] is a list of Taylor - -- series coefficients, such that f(x)=s0 + s1*x + s2*x^2.... - -- - -- Require: No Taylor coefficient is zero - -- - zero:one:[higher m | m <- [2..]] - where - zero = (s!!0, s!!1 * x) - one = (1, -s!!2/s!!1 * x) - higher m = (1 + s!!m/s!!(m-1) * x, -s!!(m+1)/s!!m * x) - + -- + -- Convert infinite number of terms of Taylor expansion of + -- a function f(x) to an infinite continued fraction, + -- where s = [s0,s1,s2,s3....] is a list of Taylor + -- series coefficients, such that f(x)=s0 + s1*x + s2*x^2.... + -- + -- Require: No Taylor coefficient is zero + -- + zero:one:[higher m | m <- [2..]] + where + zero = (s!!0, s!!1 * x) + one = (1, -s!!2/s!!1 * x) + higher m = (1 + s!!m/s!!(m-1) * x, -s!!(m+1)/s!!m * x) + fromFraction :: Fraction -> Double fromFraction = fromRational . toRational - + ------------------------------------------------------------------ --- Category: Auxiliaries +-- Category: Auxiliaries ------------------------------------------------------------------ -fac :: Integer -> Integer +fac :: Integer -> Integer fac = product . enumFromTo 1 integerRoot2 :: Integer -> Integer integerRoot2 1 = 1 integerRoot2 x = -- - -- Biggest integer m, such that x - m^2 >= 0, - -- where x is a positive integer + -- Biggest integer m, such that x - m^2 >= 0, + -- where x is a positive integer -- integerRoot2' 0 x (x `div` 2) x where - integerRoot2' lo hi r y - | c > y = integerRoot2' lo r ((r + lo) `div` 2) y - | c == y = r - | otherwise = - if (r+1)^2 > y then - r - else - integerRoot2' r hi ((r + hi) `div` 2) y - where c = r^2 + integerRoot2' lo hi r y + | c > y = integerRoot2' lo r ((r + lo) `div` 2) y + | c == y = r + | otherwise = + if (r+1)^(2::Int) > y then + r + else + integerRoot2' r hi ((r + hi) `div` 2) y + where c = r^(2::Int) ------------------------------------------------------------------ --- Category: Class Transcendental +-- Category: Class Transcendental -- --- This class declares functions for three data types: --- Fraction, Cofraction (complex fraction) and Numerus --- - a generalization of Integer, Fraction and Cofraction. +-- This class declares functions for three data types: +-- Fraction, Cofraction (complex fraction) and Numerus +-- - a generalization of Integer, Fraction and Cofraction. ------------------------------------------------------------------ class Transcendental a where - pi' :: Fraction -> a - tan' :: Fraction -> a -> a - sin' :: Fraction -> a -> a - cos' :: Fraction -> a -> a - atan' :: Fraction -> a -> a - asin' :: Fraction -> a -> a - acos' :: Fraction -> a -> a + pi' :: Fraction -> a + tan' :: Fraction -> a -> a + sin' :: Fraction -> a -> a + cos' :: Fraction -> a -> a + atan' :: Fraction -> a -> a + asin' :: Fraction -> a -> a + acos' :: Fraction -> a -> a sqrt' :: Fraction -> a -> a - root' :: Fraction -> a-> Integer -> a - power' :: Fraction -> a -> a -> a - exp' :: Fraction -> a -> a - tanh' :: Fraction -> a -> a - sinh' :: Fraction -> a -> a - cosh' :: Fraction -> a -> a - atanh' :: Fraction -> a -> a - asinh' :: Fraction -> a -> a - acosh' :: Fraction -> a -> a - log' :: Fraction -> a -> a - decimal :: Integer -> a -> IO () + root' :: Fraction -> a-> Integer -> a + power' :: Fraction -> a -> a -> a + exp' :: Fraction -> a -> a + tanh' :: Fraction -> a -> a + sinh' :: Fraction -> a -> a + cosh' :: Fraction -> a -> a + atanh' :: Fraction -> a -> a + asinh' :: Fraction -> a -> a + acosh' :: Fraction -> a -> a + log' :: Fraction -> a -> a + decimal :: Integer -> a -> IO () ------------------------------------------------------------------- -- Everything below is the instantiation of class Transcendental -- for type Fraction. See also modules Cofra and Numerus. -- --- Category: Constants +-- Category: Constants ------------------------------------------------------------------- instance Transcendental Fraction where - + pi' eps = - -- - -- pi with accuracy eps - -- - -- Based on Ramanujan formula, as described in Ref. 3 - -- Accuracy: extremely good, 10^-19 for one term of continued - -- fraction - -- - (sqrt' eps d) / (approxCF eps (fromTaylorToCF s x)) - where - x = 1//(640320^3)::Fraction - s = [((-1)^k*(fac (6*k))//((fac k)^3*(fac (3*k))))*((a*k+b)//c) | k<-[0..]] + -- + -- pi with accuracy eps + -- + -- Based on Ramanujan formula, as described in Ref. 3 + -- Accuracy: extremely good, 10^-19 for one term of continued + -- fraction + -- + (sqrt' eps d) / (approxCF eps (fromTaylorToCF s x)) + where + x = 1//(640320^(3::Int))::Fraction + s = [((-1)^k*(fac (6*k))//((fac k)^(3::Int)*(fac (3*k))))*((a*k+b)//c) | k<-[0..]] a = 545140134 - b = 13591409 - c = 426880 - d = 10005 - + b = 13591409 + c = 426880 + d = 10005 + --------------------------------------------------------------------- --- Category: Trigonometry +-- Category: Trigonometry --------------------------------------------------------------------- - - tan' eps 0 = 0 - tan' eps (u:-:0) = 1//0 + + tan' _ 0 = 0 + tan' _ (_:-:0) = 1//0 tan' eps x - -- - -- Tangent x computed with accuracy of eps. - -- - -- Trigonometric identities are used first to reduce - -- the value of x to a value from within the range of [-pi/2,pi/2] - -- - | x >= half_pi' = tan' eps (x - ((1+m)//1)*pi) - | x <= -half_pi' = tan' eps (x + ((1+m)//1)*pi) - --- | absx > 1 = 2 * t/(1 - t^2) - | otherwise = approxCF eps (cf x) - where - absx = abs x - t = tan' eps (x/2) - m = floor ((absx - half_pi)/ pi) - pi = pi' eps - half_pi'= 158//100 - half_pi = pi * (1//2) - cf u = ((0//1,1//1):[((2*r + 1)/u, -1) | r <- [0..]]) - - sin' eps 0 = 0 - sin' eps (u:-:0)= 1//0 + -- + -- Tangent x computed with accuracy of eps. + -- + -- Trigonometric identities are used first to reduce + -- the value of x to a value from within the range of [-pi/2,pi/2] + -- + | x >= half_pi' = tan' eps (x - ((1+m)//1)*p) + | x <= -half_pi' = tan' eps (x + ((1+m)//1)*p) + --- | absx > 1 = 2 * t/(1 - t^2) + | otherwise = approxCF eps (cf x) + where + absx = abs x + _ = tan' eps (x/2) + m = floor ((absx - half_pi)/ p) + p = pi' eps + half_pi'= 158//100 + half_pi = p * (1//2) + cf u = ((0//1,1//1):[((2*r + 1)/u, -1) | r <- [0..]]) + + sin' _ 0 = 0 + sin' _ (_:-:0)= 1//0 sin' eps x = 2*t/(1 + t*t) where t = tan' eps (x/2) - cos' eps 0 = 1 - cos' eps (u:-:0)= 1//0 + cos' _ 0 = 1 + cos' _ (_:-:0)= 1//0 cos' eps x = (1 - p)/(1 + p) where - t = tan' eps (x/2) + t = tan' eps (x/2) p = t*t - + atan' eps x - -- - -- Inverse tangent of x with approximation eps - -- - | x == 1//0 = (pi' eps)/2 - | x == (-1//0) = -(pi' eps)/2 - | x == 0 = 0 - | x > 1 = (pi' eps)/2 - atan' eps (1/x) - | x < -1 = -(pi' eps)/2 - atan' eps (1/x) - | otherwise = approxCF eps ((0,x):[((2*m - 1),(m*x)^2) | m<- [1..]]) - - - asin' eps x - -- - -- Inverse sine of x with approximation eps - -- - | x == 0 = 0//1 - | abs x > 1 = 1//0 - | x == 1 = (pi' eps) *(1//2) - | x == -1 = (pi' eps) * ((-1)//2) - | otherwise = atan' eps (x / (sqrt' eps (1 - x^2))) + -- + -- Inverse tangent of x with approximation eps + -- + | x == 1//0 = (pi' eps)/2 + | x == (-1//0) = -(pi' eps)/2 + | x == 0 = 0 + | x > 1 = (pi' eps)/2 - atan' eps (1/x) + | x < -1 = -(pi' eps)/2 - atan' eps (1/x) + | otherwise = approxCF eps ((0,x):[((2*m - 1),(m*x)^(2::Int)) | m<- [1..]]) - - acos' eps x - -- - -- Inverse cosine of x with approximation eps - -- - | x == 0 = (pi' eps)*(1//2) - | abs x > 1 = 1//0 - | x == 1 = 0//1 - | x == -1 = pi' eps - | otherwise = atan' eps ((sqrt' eps (1 - x^2)) / x) - + + asin' eps x + -- + -- Inverse sine of x with approximation eps + -- + | x == 0 = 0//1 + | abs x > 1 = 1//0 + | x == 1 = (pi' eps) *(1//2) + | x == -1 = (pi' eps) * ((-1)//2) + | otherwise = atan' eps (x / (sqrt' eps (1 - x^(2::Int)))) + + + acos' eps x + -- + -- Inverse cosine of x with approximation eps + -- + | x == 0 = (pi' eps)*(1//2) + | abs x > 1 = 1//0 + | x == 1 = 0//1 + | x == -1 = pi' eps + | otherwise = atan' eps ((sqrt' eps (1 - x^(2::Int))) / x) + --------------------------------------------------------------------- --- Category: Roots +-- Category: Roots --------------------------------------------------------------------- - + sqrt' eps x -- - -- Square root of x with approximation eps - -- - -- The CF pattern is: [(m,x-m^2),(2m,x-m^2),(2m,x-m^2)....] - -- where m is the biggest integer such that x-m^2 >= 0 - -- - | x == 1//0 = 1//0 - | x < 0 = 1//0 - | x == 0 = 0 - | x < 1 = 1/(sqrt' eps (1/x)) - | otherwise = approxCF eps ((m,x-m^2):[(2*m,x-m^2) | r<-[0..]]) - where - m = (integerRoot2 (floor x))//1 - - root' eps x k - -- - -- k-th root of positive number x with approximation eps - -- - | x == (1//0) = 1//0 - | x < 0 = 1//0 - | x == 0 = 0 - | k == 0 = 1//0 - | otherwise = exp' eps ((log' eps x) * (1//k)) - + -- Square root of x with approximation eps + -- + -- The CF pattern is: [(m,x-m^2),(2m,x-m^2),(2m,x-m^2)....] + -- where m is the biggest integer such that x-m^2 >= 0 + -- + | x == 1//0 = 1//0 + | x < 0 = 1//0 + | x == 0 = 0 + | x < 1 = 1/(sqrt' eps (1/x)) + | otherwise = approxCF eps ((m,x-m^(2::Int)):[(2*m,x-m^(2::Int)) | _<-[(0::Integer)..]]) + where + m = (integerRoot2 (floor x))//1 + root' eps x k + -- + -- k-th root of positive number x with approximation eps + -- + | x == (1//0) = 1//0 + | x < 0 = 1//0 + | x == 0 = 0 + | k == 0 = 1//0 + | otherwise = exp' eps ((log' eps x) * (1//k)) + + --------------------------------------------------------------------- --- Category: Powers +-- Category: Powers --------------------------------------------------------------------- - power' eps x y - -- - -- x to power of y with approximation eps - -- - | x == (1//0) = 1//0 - | x < 0 = 1//0 - | x == 0 = 0 - | y == 0 = 1 - | y == (1//0) = 1//0 - | y == (-1//0) = 0 - | otherwise = exp' eps (y * (log' eps x)) - + power' eps x y + -- + -- x to power of y with approximation eps + -- + | x == (1//0) = 1//0 + | x < 0 = 1//0 + | x == 0 = 0 + | y == 0 = 1 + | y == (1//0) = 1//0 + | y == (-1//0) = 0 + | otherwise = exp' eps (y * (log' eps x)) + --------------------------------------------------------------------- --- Category: Exponentials and hyperbolics +-- Category: Exponentials and hyperbolics --------------------------------------------------------------------- - exp' eps x - -- - -- Exponent of x with approximation eps - -- - -- Based on Jacobi type continued fraction for exponential, - -- with fractional terms: - -- n == 0 ==> (1,x) - -- n == 1 ==> (1 -x/2, x^2/12) - -- n >= 2 ==> (1, x^2/(16*n^2 - 4)) - -- For x outside [-1,1] apply identity exp(x) = (exp(x/2))^2 - -- - | x == 1//0 = 1//0 - | x == (-1//0) = 0 - | x == 0 = 1 - | x > 1 = (approxCF eps (f (x*(1//p))))^p - | x < (-1) = (approxCF eps (f (x*(1//q))))^q - | otherwise = approxCF eps (f x) - where - p = ceiling x - q = -(floor x) - f y = (1,y):(1-y/2,y^2/12):[(1,y^2/(16*n^2-4)) | n<-[2..]] - - + exp' eps x + -- + -- Exponent of x with approximation eps + -- + -- Based on Jacobi type continued fraction for exponential, + -- with fractional terms: + -- n == 0 ==> (1,x) + -- n == 1 ==> (1 -x/2, x^2/12) + -- n >= 2 ==> (1, x^2/(16*n^2 - 4)) + -- For x outside [-1,1] apply identity exp(x) = (exp(x/2))^2 + -- + | x == 1//0 = 1//0 + | x == (-1//0) = 0 + | x == 0 = 1 + | x > 1 = (approxCF eps (f (x*(1//p))))^p + | x < (-1) = (approxCF eps (f (x*(1//q))))^q + | otherwise = approxCF eps (f x) + where + p = ceiling x + q = -(floor x) + f y = (1,y):(1-y/2,y^(2::Int)/12):[(1,y^(2::Int)/(16*n^(2::Int)-4)) | n<-[2..]] + + cosh' eps x = - -- - -- Hyperbolic cosine with approximation eps - -- - (a + b)*(1//2) - where - a = exp' eps x - b = 1/a + -- + -- Hyperbolic cosine with approximation eps + -- + (a + b)*(1//2) + where + a = exp' eps x + b = 1/a sinh' eps x = - -- - -- Hyperbolic sine with approximation eps - -- - (a - b)*(1//2) - where - a = exp' eps x - b = 1/a + -- + -- Hyperbolic sine with approximation eps + -- + (a - b)*(1//2) + where + a = exp' eps x + b = 1/a tanh' eps x = - -- - -- Hyperbolic tangent with approximation eps - -- - (a - b)/ (a + b) - where - a = exp' eps x - b = 1/a + -- + -- Hyperbolic tangent with approximation eps + -- + (a - b)/ (a + b) + where + a = exp' eps x + b = 1/a - atanh' eps x - -- - -- Inverse hyperbolic tangent with approximation eps - -- - - | x >= 1 = 1//0 - | x <= -1 = -1//0 - | otherwise = (1//2) * (log' eps ((1 + x) / (1 - x))) - - asinh' eps x - -- - -- Inverse hyperbolic sine - -- - | x == 1//0 = 1//0 - | x == -1//0 = -1//0 - | otherwise = log' eps (x + (sqrt' eps (x^2 + 1))) - + atanh' eps x + -- + -- Inverse hyperbolic tangent with approximation eps + -- + + | x >= 1 = 1//0 + | x <= -1 = -1//0 + | otherwise = (1//2) * (log' eps ((1 + x) / (1 - x))) + + asinh' eps x + -- + -- Inverse hyperbolic sine + -- + | x == 1//0 = 1//0 + | x == -1//0 = -1//0 + | otherwise = log' eps (x + (sqrt' eps (x^(2::Int) + 1))) + acosh' eps x - -- - -- Inverse hyperbolic cosine - -- - | x == 1//0 = 1//0 - | x < 1 = 1//0 - | otherwise = log' eps (x + (sqrt' eps (x^2 - 1))) - + -- + -- Inverse hyperbolic cosine + -- + | x == 1//0 = 1//0 + | x < 1 = 1//0 + | otherwise = log' eps (x + (sqrt' eps (x^(2::Int) - 1))) + --------------------------------------------------------------------- --- Category: Logarithms +-- Category: Logarithms --------------------------------------------------------------------- log' eps x - -- - -- Natural logarithm of strictly positive x - -- - -- Based on Stieltjes type continued fraction for log (1+y) - -- (0,y):(1,y/2):[(1,my/(4m+2)),(1,(m+1)y/(4m+2)),.... - -- (m >= 1, two elements per m) - -- Efficient only for x close to one. For larger x we recursively - -- apply the identity log(x) = log(x/2) + log(2) - -- - | x == 1//0 = 1//0 - | x <= 0 = -1//0 - | x < 1 = -log' eps (1/x) - | x == 1 = 0 - | otherwise = - case (scaled (x,0)) of - (1,s) -> (s//1) * approxCF eps (series 1) - (y,0) -> approxCF eps (series (y-1)) - (y,s) -> approxCF eps (series (y-1)) + (s//1)*approxCF eps (series 1) - where + -- + -- Natural logarithm of strictly positive x + -- + -- Based on Stieltjes type continued fraction for log (1+y) + -- (0,y):(1,y/2):[(1,my/(4m+2)),(1,(m+1)y/(4m+2)),.... + -- (m >= 1, two elements per m) + -- Efficient only for x close to one. For larger x we recursively + -- apply the identity log(x) = log(x/2) + log(2) + -- + | x == 1//0 = 1//0 + | x <= 0 = -1//0 + | x < 1 = -log' eps (1/x) + | x == 1 = 0 + | otherwise = + case (scaled (x,0)) of + (1,s) -> (s//1) * approxCF eps (series 1) + (y,0) -> approxCF eps (series (y-1)) + (y,s) -> approxCF eps (series (y-1)) + (s//1)*approxCF eps (series 1) + where series :: Fraction -> CF series u = (0,u):(1,u/2):[(1,u*((m+n)//(4*m + 2)))|m<-[1..],n<-[0,1]] - scaled :: (Fraction,Integer) -> (Fraction, Integer) - scaled (x, n) - | x == 2 = (1,n+1) - | x < 2 = (x, n) - | otherwise = scaled (x*(1//2), n+1) + scaled :: (Fraction,Integer) -> (Fraction, Integer) + scaled (y, n) + | y == 2 = (1,n+1) + | y < 2 = (y, n) + | otherwise = scaled (y*(1//2), n+1) - + --------------------------------------------------------------------- --- Category: IO +-- Category: IO --------------------------------------------------------------------- - decimal n (u:-:0) = putStr (show u++"//0") + decimal _ (u:-:0) = putStr (show u++"//0") decimal n x - -- - -- Print Fraction with an accuracy to n decimal places, - -- or symbols +/- 1//0 for infinities. - | n <= 0 = decimal 1 x - | x < 0 = putStr (g (-v*10) (den x) n ("-"++show (-u) ++".")) - | otherwise = putStr (g (v*10) (den x) n (show u++".")) - where - (u, v) = quotRem (num x) (den x) - g x y 0 str = str - g x y n str = - case (p, q) of - (_,0) -> str ++ show p - (_,_) -> g (q*10) y (n-1) (str ++ show p) - where - (p, q) = quotRem x y - - - + -- + -- Print Fraction with an accuracy to n decimal places, + -- or symbols +/- 1//0 for infinities. + | n <= 0 = decimal 1 x + | x < 0 = putStr (g (-v*10) (den x) n ("-"++show (-u) ++".")) + | otherwise = putStr (g (v*10) (den x) n (show u++".")) + where + (u, v) = quotRem (num x) (den x) + g _ _ 0 str = str + g y z m str = + case (p, q) of + (_,0) -> str ++ show p + (_,_) -> g (q*10) z (m-1) (str ++ show p) + where + (p, q) = quotRem y z + + + --------------------------------------------------------------------------- -- References: -- @@ -618,46 +618,46 @@ -- http://www.mathsoft.com/asolve/constant/cntfrc/cntfrc.html -- 3. "Efficient on-line computation of real functions using exact floating -- point", by Peter John Potts, Imperial College --- http://theory.doc.ic.ac.uk/~pjp/ieee.html +-- http://theory.doc.ic.ac.uk/~pjp/ieee.html -------------------------------------------------------------------------- -------------------------------------------------------------------------- --- The following representation of continued fractions is used: +-- The following representation of continued fractions is used: -- --- Continued fraction: CF representation: --- ================== ==================== --- b0 + a0 +-- Continued fraction: CF representation: +-- ================== ==================== +-- b0 + a0 -- ------- ==> [(b0, a0), (b1, a1), (b2, a2).....] -- b1 + a1 -- ------- -- b2 + ... -- --- where "a's" and "b's" are Fractions. --- --- Many continued fractions could be represented by much simpler form --- [b1,b2,b3,b4..], where all coefficients "a" would have the same value 1 --- and would not need to be explicitely listed; and the coefficients "b" --- could be chosen as integers. --- However, there are some useful continued fractions that are --- given with fraction coefficients: "a", "b" or both. --- A fractional form can always be converted to an integer form, but --- a conversion process is not always simple and such an effort is not --- always worth of the achieved savings in the storage space or the --- computational efficiency. +-- where "a's" and "b's" are Fractions. -- +-- Many continued fractions could be represented by much simpler form +-- [b1,b2,b3,b4..], where all coefficients "a" would have the same value 1 +-- and would not need to be explicitely listed; and the coefficients "b" +-- could be chosen as integers. +-- However, there are some useful continued fractions that are +-- given with fraction coefficients: "a", "b" or both. +-- A fractional form can always be converted to an integer form, but +-- a conversion process is not always simple and such an effort is not +-- always worth of the achieved savings in the storage space or the +-- computational efficiency. +-- ---------------------------------------------------------------------------- -- -- Copyright: -- --- (C) 1998 Numeric Quest, All rights reserved +-- (C) 1998 Numeric Quest, All rights reserved -- -- -- --- http://www.numeric-quest.com +-- http://www.numeric-quest.com -- -- License: -- --- GNU General Public License, GPL --- +-- GNU General Public License, GPL +-- ----------------------------------------------------------------------------- diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/LinearAlgorithms.hs b/LinearAlgorithms.hs --- a/LinearAlgorithms.hs +++ b/LinearAlgorithms.hs @@ -1,10 +1,9 @@ - ------------------------------------------------------------------------------ --- Haskell module: LinearAlgorithms --- Date: initialized 2001-03-25, last modified 2001-04-01 --- Author: Jan Skibinski, Numeric Quest Inc. --- Location: http://www.numeric-quest.com/haskell/LinearAlgorithms.hs --- See also: http://www.numeric-quest.com/haskell/Orthogonals.html +-- Haskell module: LinearAlgorithms +-- Date: initialized 2001-03-25, last modified 2001-04-01 +-- Author: Jan Skibinski, Numeric Quest Inc. +-- Location: http://www.numeric-quest.com/haskell/LinearAlgorithms.hs +-- See also: http://www.numeric-quest.com/haskell/Orthogonals.html -- -- Description: -- This module provides several _selected_ linear algebra algorithms, @@ -15,14 +14,14 @@ -- Here is a list of implemented algorithms: -- --- + triangular A => R where R is upper triangular --- + triangular2 A => (R, Q) such that R = Q' A Q +-- + triangular A => R where R is upper triangular +-- + triangular2 A => (R, Q) such that R = Q' A Q -- --- + tridiagonal H => T where H is Hermitian and T is --- + tridiagonal2 H => (T, Q) tridiagonal, such that T = Q' H Q +-- + tridiagonal H => T where H is Hermitian and T is +-- + tridiagonal2 H => (T, Q) tridiagonal, such that T = Q' H Q -- --- + subsAnnihilator A => Q such that Q A has zeroed subdiagonals --- + reflection x => y where y is a complex reflection of x +-- + subsAnnihilator A => Q such that Q A has zeroed subdiagonals +-- + reflection x => y where y is a complex reflection of x -- -- Other algoritms, such as solution of linear equations are, at this time, -- imported from module Orthogonals. The latter also deals with triangulization, @@ -36,10 +35,10 @@ ---------------------------------------------------------------------------- module LinearAlgorithms ( - triangular, - triangular2, - tridiagonal, - tridiagonal2, + triangular, + triangular2, + tridiagonal, + tridiagonal2, Scalar,) where import Complex @@ -50,198 +49,198 @@ ---------------------------------------------------------------------------- -- Category: Iterative triangularization -- --- triangular A => R where R is upper triangular --- triangular2 A => (R, Q) such that R = Q' A Q +-- triangular A => R where R is upper triangular +-- triangular2 A => (R, Q) such that R = Q' A Q ---------------------------------------------------------------------------- mult :: [[Scalar]] -> [[Scalar]] -> [[Scalar]] a `mult` b - -- A matrix-product of matrices 'a' and 'b' - -- C = A B - -- where all matrices are represented as lists - -- of scalar columns - = matrix_matrix' (transposed a) b + -- A matrix-product of matrices 'a' and 'b' + -- C = A B + -- where all matrices are represented as lists + -- of scalar columns + = matrix_matrix' (transposed a) b triangular :: Int -> [[Scalar]] -> [[Scalar]] triangular n a - -- A (hopefully) triangular matrix R = Q' A Q obtained by - -- 'n' similarity transformations S(k) of matrix A: - -- Q = S1 S2 S3 .... + -- A (hopefully) triangular matrix R = Q' A Q obtained by + -- 'n' similarity transformations S(k) of matrix A: + -- Q = S1 S2 S3 .... -- -- If matrix A is Hermitian then the result is close -- to a diagonal matrix for sufficiently large n. - | n == 0 = a + | n == 0 = a | otherwise = triangular (n - 1) a1 where - a1 = (q' `mult` a ) `mult` q - q' = subsAnnihilator 0 a - q = adjoint q' - + a1 = (q' `mult` a ) `mult` q + q' = subsAnnihilator 0 a + q = adjoint q' + triangular2 :: Int -> [[Scalar]] -> ([[Scalar]], [[Scalar]]) triangular2 n a - -- A pair of matrices (R, Q) obtained by 'n' - -- similarity transformations, where R = Q' A Q - -- is a (hopefully) triangular matrix, or diagonal - -- if A is Hermitian. The transformation matrix Q - -- is required for computation of eigenvectors - -- of A. + -- A pair of matrices (R, Q) obtained by 'n' + -- similarity transformations, where R = Q' A Q + -- is a (hopefully) triangular matrix, or diagonal + -- if A is Hermitian. The transformation matrix Q + -- is required for computation of eigenvectors + -- of A. = triangular2' n a (unit_matrix n) where - triangular2' n a p - | n == 0 = (a, p) - | otherwise = triangular2' (n - 1) a1 p1 - where - a1 = (q' `mult` a ) `mult` q - p1 = p `mult` q - q' = subsAnnihilator 0 a - q = adjoint q' - + triangular2' o b p + | o == 0 = (b, p) + | otherwise = triangular2' (o - 1) b1 p1 + where + b1 = (q' `mult` b ) `mult` q + p1 = p `mult` q + q' = subsAnnihilator 0 b + q = adjoint q' + ---------------------------------------------------------------------------- -- Category: Tridiagonalization of a Hermitian matrix -- --- + tridiagonal H -> T where H is Hermitian and T is tridiagonal --- + tridiagonal2 H -> (T, Q) such that T = Q' H Q +-- + tridiagonal H -> T where H is Hermitian and T is tridiagonal +-- + tridiagonal2 H -> (T, Q) such that T = Q' H Q ---------------------------------------------------------------------------- tridiagonal :: [[Scalar]] -> [[Scalar]] tridiagonal h - -- A tridiagonal matrix T = Q' H Q, obtained from Hermitian - -- matrix H by a finite number of elementary similarity - -- transformations (Householder reductions). - | n < 3 = h - | otherwise = f (tail es) h 1 + -- A tridiagonal matrix T = Q' H Q, obtained from Hermitian + -- matrix H by a finite number of elementary similarity + -- transformations (Householder reductions). + | n < 3 = h + | otherwise = f (tail es) h 1 where - n = length h - es = unit_matrix n - - f bs a k - | length bs == 1 = a - | otherwise = f (tail bs) a1 (k+1) - where - a1 = (q' `mult` a) `mult` q - q' = [r e | e <- es] - q = adjoint q' - r = reflection u (head bs) - u = replicate k 0 ++ drop k (a!!(k-1)) + n = length h + es = unit_matrix n + f bs a k + | length bs == 1 = a + | otherwise = f (tail bs) a1 (k+1) + where + a1 = (q' `mult` a) `mult` q + q' = [r e | e <- es] + q = adjoint q' + r = reflection u (head bs) + u = replicate k 0 ++ drop k (a!!(k-1)) + tridiagonal2 :: [[Scalar]] -> ([[Scalar]], [[Scalar]]) tridiagonal2 h - -- A pair (T, Q) of matrices, obtained from - -- similarity transformation of Hermitian matrix H - -- where T = Q' H Q is a tridiagonal matrix and Q is unitary - -- transformation made of a finite product of - -- elementary Householder reductions. - | n < 3 = (h, es) - | otherwise = f (tail es) h es 1 + -- A pair (T, Q) of matrices, obtained from + -- similarity transformation of Hermitian matrix H + -- where T = Q' H Q is a tridiagonal matrix and Q is unitary + -- transformation made of a finite product of + -- elementary Householder reductions. + | n < 3 = (h, es) + | otherwise = f (tail es) h es 1 where - n = length h - es = unit_matrix n - - f bs a p k - | length bs == 1 = (a, p) - | otherwise = f (tail bs) a1 p1 (k+1) - where - a1 = (q' `mult` a) `mult` q - q' = [r e | e <- es] - q = adjoint q' - p1 = p `mult` q - r = reflection u (head bs) - u = replicate k 0 ++ drop k (a!!(k-1)) + n = length h + es = unit_matrix n + f bs a p k + | length bs == 1 = (a, p) + | otherwise = f (tail bs) a1 p1 (k+1) + where + a1 = (q' `mult` a) `mult` q + q' = [r e | e <- es] + q = adjoint q' + p1 = p `mult` q + r = reflection u (head bs) + u = replicate k 0 ++ drop k (a!!(k-1)) + ---------------------------------------------------------------------------- -- Category: Elementary unitary transformations -- --- + subsAnnihilator A => Q such that Q A has zeroed subdiagonals --- + reflection x => y where y is a complex reflection of x +-- + subsAnnihilator A => Q such that Q A has zeroed subdiagonals +-- + reflection x => y where y is a complex reflection of x ---------------------------------------------------------------------------- subsAnnihilator :: Int -> [[Scalar]] -> [[Scalar]] subsAnnihilator k a - -- A unitary matrix Q' transforming any n x n - -- matrix A to an upper matrix B, which has - -- zero values below its 'k'-th subdiagonal - -- (annihilates all subdiagonals below k-th) - -- B = Q' A - -- where - -- 'a' is a list of columns of matrix A + -- A unitary matrix Q' transforming any n x n + -- matrix A to an upper matrix B, which has + -- zero values below its 'k'-th subdiagonal + -- (annihilates all subdiagonals below k-th) + -- B = Q' A + -- where + -- 'a' is a list of columns of matrix A -- - -- If k=0 then B is an upper triangular matrix, - -- if k=1 then B is an upper Hessenberg matrix. - -- The transformation Q is built from n - k - 1 - -- elementary Householder transformations of - -- the first n-k-1 columns of iteratively transformed - -- matrix A. - | n < 2 + k = es - | otherwise = f (drop k es) a1 es k + -- If k=0 then B is an upper triangular matrix, + -- if k=1 then B is an upper Hessenberg matrix. + -- The transformation Q is built from n - k - 1 + -- elementary Householder transformations of + -- the first n-k-1 columns of iteratively transformed + -- matrix A. + | n < 2 + k = es + | otherwise = f (drop k es) a1 es k where - n = length a - es = unit_matrix n - a1 = take (n - 1 - k) a + n = length a + es = unit_matrix n + a1 = take (n - 1 - k) a - f bs a p k - | length bs == 1 = p - | otherwise = f (tail bs) a1 p1 (k+1) - where - a1 = [r v |v <- tail a] - p1 = q' `mult` p - q' = [r e | e <- es] - r = reflection u (head bs) - u = replicate k 0 ++ drop k (head a) + f bs b p l + | length bs == 1 = p + | otherwise = f (tail bs) b1 p1 (l+1) + where + b1 = [r v |v <- tail b] + p1 = q' `mult` p + q' = [r e | e <- es] + r = reflection u (head bs) + u = replicate k 0 ++ drop l (head b) reflection :: [Scalar] -> [Scalar] -> [Scalar] -> [Scalar] reflection a e x - -- A vector resulting from unitary complex - -- Householder-like transformation of vector 'x'. + -- A vector resulting from unitary complex + -- Householder-like transformation of vector 'x'. -- - -- The operator of such transformation is defined - -- by mapping vector 'a' to a multiple 'p' of vector 'e' - -- U |a > = p | e > - -- where scalar 'p' is chosen to guarantee unitarity - -- < a | a > = < p e | p e>. + -- The operator of such transformation is defined + -- by mapping vector 'a' to a multiple 'p' of vector 'e' + -- U |a > = p | e > + -- where scalar 'p' is chosen to guarantee unitarity + -- < a | a > = < p e | p e>. -- - -- This transformation is not generally Hermitian, because - -- the scalar 'p' might become complex - unless - -- < a | e > = < e | a >, - -- which is the case when both vectors are real, and - -- when this transformation becomes a simple Hermitian - -- reflection operation. - -- See reference [1] for details. + -- This transformation is not generally Hermitian, because + -- the scalar 'p' might become complex - unless + -- < a | e > = < e | a >, + -- which is the case when both vectors are real, and + -- when this transformation becomes a simple Hermitian + -- reflection operation. + -- See reference [1] for details. -- | d == 0 = x | otherwise = [xk - z * yk |(xk, yk) <- zip x y] where - z = s * bra_ket y x - s = 2/h :+ (-2 * g)/h - h = 1 + g^2 - g = imagPart a_b / d - d = a_a - realPart a_b - y = normalized [ak - bk |(ak, bk) <- zip a b] - p = a_a / (realPart (bra_ket e e)) - b = map ((sqrt p :+ 0) * ) e - a_a = realPart (bra_ket a a) - a_b = bra_ket a b + z = s * bra_ket y x + s = 2/h :+ (-2 * g)/h + h = 1 + g^(2::Int) + g = imagPart a_b / d + d = a_a - realPart a_b + y = normalized [ak - bk |(ak, bk) <- zip a b] + p = a_a / (realPart (bra_ket e e)) + b = map ((sqrt p :+ 0) * ) e + a_a = realPart (bra_ket a a) + a_b = bra_ket a b ---------------------------------------------------------------------------- -- Category: Test data -- ---------------------------------------------------------------------------- -matrixA :: [[Scalar]] -matrixA - -- Test matrix A represented as list of scalar columns. - = [ - [1, 2, 4, 1, 5] - , [2, 3, 2, 6, 4] - , [4, 2, 5, 2, 3] - , [1, 6, 2, 7, 2] - , [5, 4, 3, 2, 9] - ] +-- matrixA :: [[Scalar]] +-- matrixA +-- -- Test matrix A represented as list of scalar columns. +-- = [ +-- [1, 2, 4, 1, 5] +-- , [2, 3, 2, 6, 4] +-- , [4, 2, 5, 2, 3] +-- , [1, 6, 2, 7, 2] +-- , [5, 4, 3, 2, 9] +-- ] ---------------------------------------------------------------------------- -- Module documentation @@ -267,7 +266,7 @@ -- The Schur theorem states that there exists a unitary matrix Q such -- that any nonsingular matrix A can be transformed to an upper triangular -- matrix R via similarity transformation --- R = Q' A Q +-- R = Q' A Q -- which preserves the eigenvalues. Here Q' stands for a Hermitian -- conjugate of Q (adjoint, or Q-dagger). @@ -275,7 +274,7 @@ -- elements, finding such transformation solves the first part of -- the eigenproblem. The second part, finding the eigenvectors of A, -- is trivial since they can be computed from eigenvectors of R: --- | x(A) > = Q | x(R) > +-- | x(A) > = Q | x(R) > -- -- In particular, when matrix A is Hermitian, then the matrix R -- becomes diagonal, and the eigenvectors of R are its normalized @@ -291,7 +290,7 @@ -- The methods known in literature as QR factorization (decomposition) -- methods iteratively compose such unitary matrix Q from a series of -- elementary unitary transformations, Q(1), Q(2)..: --- Q = Q(1) Q(2) Q(3) ... +-- Q = Q(1) Q(2) Q(3) ... -- The most popular method of finding those elementary unitary -- transformations relies on a reflection transformation, so selected as -- to zero out all components of the matrix below its main diagonal. Our @@ -312,8 +311,8 @@ -- if you are interested in eigenvalues only. But when you need -- the eigenvectors as well, then use the function 'triangular2', -- which also delivers the transformation Q, as shown below: --- triangular A => R where R is upper triangular --- triangular2 A => (R, Q) such that R = Q' A Q +-- triangular A => R where R is upper triangular +-- triangular2 A => (R, Q) such that R = Q' A Q -- -- Tridiagonalization of Hermitian matrices: -- ----------------------------------------- @@ -339,8 +338,8 @@ -- -- Again, we are offering here two versions of the tridiagonalization: -- --- + tridiagonal H -> T where H is Hermitian and T is tridiagonal --- + tridiagonal2 H -> (T, Q) such that T = Q' H Q +-- + tridiagonal H -> T where H is Hermitian and T is tridiagonal +-- + tridiagonal2 H -> (T, Q) such that T = Q' H Q -- -- Elementary transformations: -- --------------------------- @@ -353,27 +352,27 @@ -- all matrix components below its first subdiagonal - leading to a so-called -- upper Hessenberg matrix. -- --- + subsAnnihilator A => Q such that Q A has zeroed subdiagonals --- + reflection x => y where y is a complex reflection of x +-- + subsAnnihilator A => Q such that Q A has zeroed subdiagonals +-- + reflection x => y where y is a complex reflection of x -- ---------------------------------------------------------------------------- -- References: --- [1] Xiaobai Sun, On Elementary Unitary and Phi-unitary transformations, --- Duke University, Department Of Computer Science, 1995, --- http://citeseer.nj.nec.com/340881.html +-- [1] Xiaobai Sun, On Elementary Unitary and Phi-unitary transformations, +-- Duke University, Department Of Computer Science, 1995, +-- http://citeseer.nj.nec.com/340881.html --------------------------------------------------------------------------- -- -- Copyright: -- --- (C) 2001 Numeric Quest, All rights reserved +-- (C) 2001 Numeric Quest, All rights reserved -- -- Email: jans@numeric-quest.com -- --- http://www.numeric-quest.com +-- http://www.numeric-quest.com -- -- License: -- --- GNU General Public License, GPL +-- GNU General Public License, GPL -- --------------------------------------------------------------------------- diff --git a/Orthogonals.lhs b/Orthogonals.lhs --- a/Orthogonals.lhs +++ b/Orthogonals.lhs @@ -3,1857 +3,1864 @@ - Indexless linear algebra algorithms - - - -
    -
    -

    - *** -

    -

    - Indexless linear algebra algorithms -

    - -
    - Orthogonalization, linear equations, eigenvalues and eigenvectors -
    - Literate Haskell module Orthogonals.lhs -     -

    - Jan Skibinski, - Numeric Quest Inc., Huntsville, Ontario, Canada -

    - 1998.09.19, last modified 1998.12.28 -

    -
    -

    - It has been argued that the functional paradigm offers more - support for scientific computing than the traditional imperative - programming, such as greater similarity of functional implementation - to mathematical specification of a problem. However, efficiency - of scientific algorithms implemented in Haskell is very low compared - to efficiencies of C or Fortran implementations - notwithstanding - the exceptional descriptive power of Haskell. -

    - It has been also argued that tradition and inertia are partially - responsible for this sore state and that many functional algorithms - are direct translations of their imperative counterparts. -
    - Arrays - with their indexing schemes and destructive updating - are basic tools of imperative programming. But pure functional - languages, which prohibit variable reassignments, cannot compete - with imperative languages by using the same tools and following - similar reasoning and patterns - unless the functional arrays - themselves are designed with performance in mind. This is - a case with Clean, where efficiency of one kind of their arrays - -- strict unboxed array, approaches efficiency of C. -
    - But this has not been done for Haskell arrays yet. They are - lazy, boxed and use auxilliary association lists (index, value) - for initialization -- the latter being mostly responsible for - low efficiency of those algorithms that create many interim - arrays. -
    - It appears, that -- as long as indexing scheme is not used - for lookups and updates -- Haskell lists are more efficient - than arrays -- at least at the currents state of Haskell. -

    - With this in mind, we are attempting to demonstrate here - that the indexing traps can be successfully avoided. - This module implements afresh several typical problems from linear - algebra. Standard Haskell lists are employed instead of arrays - and not a single algorithm ever uses indices for lookups - or updates. -

    - We do not claim high efficiency of these algorithms; consider - them exploratory. However, we do claim that the clarity of - these algorithms is significantly better than of those functionally - similar algorithms that employ indexing schemes. -

    - Two major algorithms have been invented and implemented in Haskell: - one for solving systems of linear equations and one for finding - eigenvalues and eigenvectors of almost any type of a square matrix. - This includes symmetric, hermitian, general complex or nonsymmetric - matrices with real eigenvalues. -

    - Amazingly, both methods are based on the same factorization, akin - to QR method, but not exactly the same as the standard QR one. - A simple trick allows to extend this method to nonsymmetric real - matrices with complex eigenvalues and thus one method applies to - all types of matrices. - It follows that the eigenvalue/eigenvector problem can be consistently - treated all across the board. In addition, no administrative - (housekeeping) boring trivia is required here and that helps to - clearly explain the mechanisms employed. - - -

    -

    -

    - - Contents - -

    -

      -

    • - Notation -

    • - Scalar products and vector normalization -
        -
      • - bra_ket, scalar product -
      • - sum_product, a cousin of bra_ket -
      • - norm, vector norm -
      • - normalized, vector normalized to one -
      -

    • - Transposition and adjoining of matrices -
        -
      • - transposed, transposed matrix -
      • - adjoint, transposed and conjugated matrix -
      -

    • - Products involving matrices -
        -
      • - matrix_matrix, product of two matrices as list of rows -
      • - matrix_matrix', product of two matrices as list of columns -
      • - triangle_matrix', upper triangular matrix times square matrix -
      • - matrix_ket, matrix times ket vector -
      • - bra_matrix, bra vector times matrix -
      • - bra_matrix_ket, matrix multiplied on both sides by vectors -
      • - scalar_matrix, scalar times matrix -
      -

    • - Orthogonalization process -
        -
      • - orthogonals, set of orthogonal vectors -
      • - gram_schmidt, vector perpendicular to a hyperplane -
      - -

    • - Solutions of linear equations by orthogonalization -
        -
      • - one_ket_triangle, triangularization of one vector equation -
      • - one_ket_solution, solution for one unknown vector -
      • - many_kets_triangle, triangularization of several vector equations -
      • - many_kets_solution, solution for several unknown vectors -
      -

    • - Matrix inversion -
        -
      • - inverse, inverse of a matrix -
      -

    • - QR factorization of matrices provided by "many_kets_triangle" -
        -
      • - factors_QR, QR alike factorization of matrices -
      • - determinant, computation of the determinant based on the QR factorization -
      -

    • - Similarity transformations and eigenvalues -
        -
      • - similar_to, matrix obtained by similarity transformation -
      • - iterated_eigenvalues, list of approximations of eigenvalues -
      • - eigenvalues, final approximation of eigenvalues -
      -

    • - Preconditioning of real nonsymmetric matrices -
        -
      • - add_to_diagonal, simple preconditioning method -
      -

    • - Examples of iterated eigenvalues -
        -
      • - Symmetric real matrix -
      • - Hermitian complex matrix -
      • - General complex matrix -
      • - Nonsymmetric real matrix with real eigenvalues -
      • - Nonsymmetric real matrix with complex eigenvalues -
      -

    • - Eigenvectors for distinct eigenvalues -
        -
      • - eigenkets, eigenvectors for distinct eigenvalues -
      -

    • - Eigenvectors for degenerated eigenvalues -
        -
      • - eigenket', eigenvector based on a trial vector -
      - -

    • - Auxiliary functions -
        -
      • - unit_matrix, a unit matrix with 1's on a diagonal -
      • - unit_vector, a vector with one non-zero componenet -
      • - diagonals, vector made of a matrix diagonal -
      -
    - -

    -

    -

    - - Notation - -

    - What follows is written in Dirac's notation, as used - in Quantum Mechanics. Matrices are represented by capital - letters, while vectors come in two varieties: -

      -

    • - Bra vector x, written < x |, is represented by one-row matrix -

    • Ket vector y, written | y >, is represented by one-column matrix -
    -

    - Bra vectors can be obtained from ket vectors by transposition - and conjugation of their components. Conjugation is only - important for complex vectors. -

    - Scalar product of two vectors | x > and | y > is written - as -

    -	< x | y >
-
    - which looks like a bracket and is sometimes called a "bra_ket". - This justifies "bra" and "ket" names introduced by Dirac. There - is a good reason for conjugating the components of "bra-vector": - the scalar product of -
    -	< x | x >
-
    - should be a square of the norm of the vector "x", and that - means that it should be represented by a real number, or complex - number but with its imaginary part equal to zero. -

    -

    -

    -

    -
-> module Orthogonals where
-> import Complex
-> import Ratio
-> import qualified List
-
-
    - - Scalar product and vector normalization - -

    - The scalar product "bra_ket" is a basis of many algorithms - presented here. - - -

    -
-> bra_ket :: (Scalar a, Num a) => [a] -> [a] -> a
-> bra_ket u v =
->       --
->       -- Scalar product of two vectors u and v,
->       -- or < u | v > in Dirac's notation.
->       -- This is equally valid for both: real and complex vectors.
->       --
->       sum_product u (map coupled v)
-
-
    - - Notice the call to function "coupled" in the above implementation - of scalar product. This function conjugates its argument - if it is complex, otherwise does not change it. It is defined - in the class Scalar - specifically designed for this purpose - mainly. -
    - This class also defines a norm of a vector that might be used - by some algorithms. So far we have been able to avoid this. -
    -
-> class Scalar a where
->     coupled    :: a->a
->     norm       :: [a] -> a
->     almostZero :: a -> Bool
->     scaled     :: [a] -> [a]
-
-> instance Scalar Double where
->     coupled x    = x
->     norm u       = sqrt (bra_ket u u)
->     almostZero x = (abs x) < 1.0e-8
->     scaled       = scaled'
-
-> instance Scalar Float where
->    coupled x    = x
->    norm u       = sqrt (bra_ket u u)
->    almostZero x = (abs x) < 1.0e-8
->    scaled       = scaled'
-
-> instance (Integral a) => Scalar (Ratio a) where
->     coupled x    = x
->     -- norm u    = fromDouble ((sqrt (bra_ket u u))::Double)
->     -- Intended hack to silently convert to and from Double.
->     -- But I do not know how to declare it properly.
->     --
->     -- Our type Fraction, when used instead of Ratio a, has its own
->     -- definition of sqrt. No hack would be needed here.
->     almostZero x = abs x < 1e-8
->     scaled       = scaled'
-
-> instance (RealFloat a) => Scalar (Complex a) where
->     coupled (x:+y) = x:+(-y)
->     norm u         = sqrt (realPart (bra_ket u u)) :+ 0
->     almostZero z   = (realPart (abs z)) < 1.0e-8
->     scaled u       = [(x/m):+(y/m) | x:+y <- u]
->        where m = maximum [max (abs x) (abs y) | x:+y <- u]
-
-> norm1 :: (Num a) => [a] -> a
-> norm1 = sum . map abs
-
-> norminf :: (Num a, Ord a) => [a] -> a
-> norminf = maximum . map abs
-
-> matnorm1 :: (Num a, Ord a) => [[a]] -> a
-> matnorm1 = matnorminf . transposed
-
-> matnorminf :: (Num a, Ord a) => [[a]] -> a
-> matnorminf = maximum . map norm1
-
-
-
    - - But we also need a slightly different definition of - scalar product that will appear in multiplication of matrices - by vectors (or vice versa): a straightforward accumulated product - of two lists, where no complex conjugation takes place. - We will call it a 'sum_product". -
    -
-> sum_product :: Num a => [a] -> [a] -> a
-> sum_product u v =
->       --
->       -- Similar to scalar product but without
->       -- conjugations of | u > components
->       -- Used in matrix-vector or vector-matrix products
->       --
->       sum (zipWith (*) u v)
-
-
    - Some algorithms might need vectors normalized to one, although - we'll try to avoid the normalizations due to its high cost - or its inapplicability to rational numbers. Instead, we wiil - scale vectors by their maximal components. -
    -
-> normalized :: (Scalar a, Fractional a) => [a] -> [a]
-> normalized u =
->       [uk/n | uk <- u]
->       where
->           n = norm u
-
-> scaled' u =
->       [uk/um | uk <- u]
->       where
->           um = maximum [abs uk| uk <- u]
-
-
    -
    -

    - - Transposition and adjoining of matrices - -

    - Matrices are represented here by lists of lists. - Function "transposed" converts from row-wise to column-wise - representation, or vice versa. -

    - When transposition is combined with complex conjugation - the resulting matrix is called "adjoint". -

    - A square matrix is called symmetric if it is equal to its transpose -

    -	A = AT
-
    - It is called Hermitian, or self-adjoint, if it equals to - its adjoint -
    -	A = A+
-
-> transposed :: [[a]] -> [[a]]
-> transposed a
->     | null (head a) = []
->     | otherwise = ([head mi| mi <- a])
->                   :transposed ([tail mi| mi <- a])
-
-> adjoint :: Scalar a => [[a]] -> [[a]]
-> adjoint a
->     | null (head a) = []
->     | otherwise = ([coupled (head mi)| mi <- a])
->                   :adjoint ([tail mi| mi <- a])
-
-
    -

    -

    -

    - - Linear combination and sum of two matrices - -

    - One can form a linear combination of two matrices, such - as -

    -	C = alpha A + beta B
-	where
-	    alpha and beta are scalars
-
    - The most generic form of any combination, not necessary - linear, of components of two matrices is given by "matrix_zipWith" - function below, which accepts a function "f" describing such - combination. For the linear combination with two scalars - the function "f" could be defined as: -
    -	f alpha beta a b = alpha*a + beta*b
-
    - For a straightforward addition of two matrices this auxiliary - function is simply "(+)". -
    -
-> matrix_zipWith f a b =
->     --
->     -- Matrix made of a combination
->     -- of matrices a and b - as specified by f
->     --
->     [zipWith f ak bk | (ak,bk) <- zip a b]
-
-> add_matrices a b = matrix_zipWith (+)
-
-
    - -

    -

    -

    - - Products involving matrices - -

    - Variety of products involving matrices can be defined. - Our Haskell implementation is based on lists of lists - and therefore is open to interpretation: sublists - can either represent the rows or the columns of a matrix. -

    - The following definitions are somehow arbitrary, since - one can choose alternative interpretations of lists - representing matrices. -

    - - C = A B - -

    - Inner product of two matrices A B can be expressed quite simply, - providing that matrix A is represented by a list of rows - and B - by a list of columns. Function "matrix_matrix" - answers list of rows, while "matrix_matrix'" - list - of columns. -

    - Major algorithms of this module make use of "triangle_matrix'", - which calculates a product of upper triangular matrix - with square matrix and returns a rectangular list of columns. - -
    -
-> matrix_matrix :: Num a => [[a]] -> [[a]] -> [[a]]
-> matrix_matrix a b
-> --
-> -- A matrix being an inner product
-> -- of matrices A and B, where
-> --     A is represented by a list of rows a
-> --     B is represented by a list of columns b
-> --     result is represented by list of rows
-> -- Require: length of a is equal of length of b
-> -- Require: all sublists are of equal length
->
->       | null a = []
->       | otherwise = ([sum_product (head a) bi | bi <- b])
->                  : matrix_matrix (tail a) b
-
-> matrix_matrix' a b
->       --
->       -- Similar to "matrix_matrix"
->       -- but the result is represented by
->       -- a list of columns
->       --
->       | null b = []
->       | otherwise = ([sum_product ai (head b) | ai <- a])
->                    : matrix_matrix' a (tail b)
-
-
-> triangle_matrix' :: Num a => [[a]] -> [[a]] -> [[a]]
-> triangle_matrix' r q =
->       --
->       -- List of columns of of a product of
->       -- upper triangular matrix R and square
->       -- matrix Q
->       -- where
->       --     r is a list of rows of R
->       --     q is a list of columns of A
->       --
->       [f r qk | qk <- q]
->       where
->           f t u
->               | null t = []
->               | otherwise = (sum_product (head t) u)
->                             : (f (tail t) (tail u))
-
-
-
-
    - - | u > = A | v > - -

    - Product of a matrix and a ket-vector is another - ket vector. The following implementation assumes - that list "a" represents rows of matrix A. -

    -
-> matrix_ket :: Num a => [[a]] -> [a] -> [a]
-> matrix_ket a v = [sum_product ai v| ai <- a]
-
-
    - - < u | = < v | A - -

    - Bra-vector multiplied by a matrix produces - another bra-vector. The implementation below - assumes that list "a" represents columns - of matrix A. It is also assumed that vector - "v" is given in its standard "ket" representation, - therefore the definition below uses "bra_ket" - instead of "sum_product". -

    -
-> bra_matrix :: (Scalar a, Num a) => [a] -> [[a]] -> [a]
-> bra_matrix v a = [bra_ket v ai | ai <- a]
-
-
    - - alpha = < u | A | v > - -

    - This kind of product results in a scalar and is often - used to define elements of a new matrix, such as -

    -	B[i,j] = < ei | A | ej >
-
    - The implementation below assumes that list "a" represents - rows of matrix A. -
    -
-> bra_matrix_ket :: (Scalar a, Num a) => [a] -> [[a]] -> [a] -> a
-> bra_matrix_ket u a v =
->     bra_ket u (matrix_ket a v)
-
-
    - - B = alpha A - -

    - Below is a function which multiplies matrix by a scalar: -

    -
-> scalar_matrix :: Num a => a -> [[a]] -> [[a]]
-> scalar_matrix alpha a =
->       [[alpha*aij| aij <- ai] | ai<-a]
-
-
    -

    -

    -

    - - Orthogonalization process - -

    - - Gram-Schmidt orthogonalization procedure is used here - for calculation of sets of mutually orthogonal vectors. -

    - Function "orthogonals" computes a set of mutually orthogonal - vectors - all orthogonal to a given vector. Such set plus - the input vector form a basis of the vector space. Another - words, they are the base vectors, although we cannot call them - unit vectors since we do not normalize them for two reasons: -
      -
    • - None of the algorithms presented here needs this -- quite - costly -- normalization. -
    • - Some algorithms can be used either with doubles or with - rationals. The neat output of the latter is sometimes desirable - for pedagogical or accuracy reasons. But normalization requires "sqrt" - function, which is not defined for rational numbers. We could - use our module Fraction instead, where "sqrt" is defined, - but we'll leave it for a future revision of this module. -
    -

    - Function "gram_schmidt" computes one vector - orthogonal - to an incomplete set of orthogonal vectors, which form a hyperplane - in the vector space. Another words, "gram_schmidt" vector is - perpendicular to such a hyperlane. - - -

    -
-> orthogonals :: (Scalar a, Fractional a) => [a] -> [[a]]
-> orthogonals x =
->       --
->       -- List of (n-1) linearly independent vectors,
->       -- (mutually orthogonal) and orthogonal to the
->       -- vector x, but not normalized,
->       -- where
->       --     n is a length of x.
->       --
->       orth [x] size (next (-1))
->       where
->           orth a n m
->               | n == 1        = drop 1 (reverse a)
->               | otherwise     = orth ((gram_schmidt a u ):a) (n-1) (next m)
->               where
->                   u = unit_vector m size
->           size = length x
->           next i = if (i+1) == k then (i+2) else (i+1)
->           k = length (takeWhile (== 0) x) -- first non-zero component of x
-
-> gram_schmidt :: (Scalar a, Fractional a) => [[a]] -> [a] -> [a]
-> gram_schmidt a u =
->       --
->       -- Projection of vector | u > on some direction
->       -- orthogonal to the hyperplane spanned by the list 'a'
->       -- of mutually orthogonal (linearly independent)
->       -- vectors.
->       --
->       gram_schmidt' a u u
->       where
->           gram_schmidt' a u v
->               | null a       = v
->               | all (== 0) e = gram_schmidt' (tail a) u v
->               | otherwise    = gram_schmidt' (tail a) u v'
->               where
->                   v' = vectorCombination v (-(bra_ket e u)/(bra_ket e e)) e
->                   e  = head a
->           vectorCombination x c y
->               | null x = []
->               | null y = []
->               | otherwise = (head x + c * (head y))
->                             : (vectorCombination (tail x) c (tail y))
-
-
    -

    -

    -

    - - Solutions of linear equations by orthogonalization - -

    - A matrix equation for unknown vector | x > -

    -	A | x > = | b >
-
    - can be rewritten as -
    -	x1 | 1 > + x2 | 2 > + x3 | 3 > + ... + xn | n > = | b >     (7.1)
-	where
-		| 1 >, | 2 >... represent columns of the matrix A
-
    - For any n-dimensional vector, such as "1", there exist - n-1 linearly independent vectors "ck" that are orthogonal to "1"; - that is, each satisfies the relation: -
    -	< ck | 1 > = 0, for k = 1...m, where m = n - 1
-
    - If we could find all such vectors, then we could multiply - the equation (7.1) by each of them, and end up with m = n-1 - following equations -
    -	< c1 | 2 > x2 + < c1 | 3 > x3 + ... < c1 | n > xn = < c1 | b >
-	< c2 | 2 > x2 + < c2 | 3 > x3 + ... < c2 | n > xn = < c2 | b >
-	.......
-	< cm | 2 > x2 + < cm | 3 > x3 + ... < cm | n > xn = < cm | b >
-
    - But the above is nothing more than a new matrix equation -
    -	A' | x' > = | b' >
-	or
-
-	x2 | 2'> + x3 | 3'> .... + xn | n'> = | b'>
-	where
-	    primed vectors | 2' >, etc. are the columns of the new
-	    matrix A'.
-
    - with the problem dimension reduced by one. - -
    - Taking as an example a four-dimensional problem and writing - down the successive transformations of the original equation - we will end up with the following triangular pattern made of - four vector equations: - -
    -	x1 | 1 > + x2 | 2 > + x3 | 3 >  + x4 | 4 >   = | b >
-		   x2 | 2'> + x3 | 3'>  + x4 | 4'>   = | b'>
-			      x3 | 3''> + x4 | 4''>  = | b''>
-					  x4 | 4'''> = | b'''>
-
    - But if we premultiply each vector equation by a non-zero vector - of our choice, say < 1 | , < 2' |, < 3'' |, and < 4''' | - chosen - correspondingly for equations 1, 2, 3 and 4, then the above - system of vector equations will be converted to much simpler - system of scalar equations. The result is - shown below in matrix representation: - -
    -	| p11  p12   p13   p14 | | x1 | = | q1 |
-	| 0    p22   p23   p24 | | x2 | = | q2 |
-	| 0    0     p33   p34 | | x3 | = | q3 |
-	| 0    0     0     p44 | | x4 | = | q4 |
-
    - In effect, we have triangularized our original matrix A. - Below is a function that does that for any problem size: -
    -
-> one_ket_triangle :: (Scalar a, Fractional a) => [[a]] -> [a] -> [([a],a)]
-> one_ket_triangle a b
->     --
->     -- List of pairs: (p, q) representing
->     -- rows of triangular matrix P and of vector | q >
->     -- in the equation P | x > = | q >, which
->     -- has been obtained by linear transformation
->     -- of the original equation A | x > = | b >
->     --
->     | null a = []
->     | otherwise = (p,q):(one_ket_triangle a' b')
->     where
->         p    = [bra_ket u ak | ak <- a]
->         q    = bra_ket u b
->         a'   = [[bra_ket ck ai | ck <- orth] | ai <- v]
->         b'   = [ bra_ket ck b  | ck <- orth]
->         orth = orthogonals u
->         u    = head a
->         v    = tail a
-
-
    - The triangular system of equations can be easily solved by - successive substitutions - starting with the last equation. - -
    -
-> one_ket_solution :: (Fractional a, Scalar a) => [[a]] -> [a] -> [a]
-> one_ket_solution a b =
->     --
->     -- List representing vector |x>, which is
->     -- a solution of the matrix equation
->     --     A |x> = |b>
->     -- where
->     --     a is a list of columns of matrix A
->     --     b is a list representing vector |b>
->     --
->     solve' (unzip (reverse (one_ket_triangle a b))) []
->     where
->         solve' (a, b) xs
->             | null a  = xs
->             | otherwise = solve' ((tail a), (tail b)) (x:xs)
->             where
->                 x = (head b - (sum_product (tail u) xs))/(head u)
->                 u = head a
-
-
    - The triangularization procedure can be easily extended - to a list of several ket-vectors | b > on the right hand - side of the original equation A | x > = | b > -- instead - of just one: -
    -
-> many_kets_triangle :: (Scalar a, Fractional a) => [[a]] -> [[a]] -> [([a],[a])]
-> many_kets_triangle a b
->     --
->     -- List of pairs: (p, q) representing
->     -- rows of triangular matrix P and of rectangular matrix Q
->     -- in the equation P X = Q, which
->     -- has been obtained by linear transformation
->     -- of the original equation A X = B
->     -- where
->     --     a is a list of columns of matrix A
->     --     b is a list of columns of matrix B
->     --
->     | null a = []
->     | otherwise = (p,q):(many_kets_triangle a' b')
->     where
->         p    = [bra_ket u ak   | ak <- a]
->         q    = [bra_ket u bk   | bk <- b]
->         a'   = [[bra_ket ck ai | ck <- orth] | ai <- v]
->         b'   = [[bra_ket ck bi | ck <- orth] | bi <- b]
->         orth = orthogonals u
->         u    = head a
->         v    = tail a
-
-
    - Similarly, function 'one_ket_solution' can be generalized - to function 'many_kets_solution' that handles cases with - several ket-vectors on the right hand side. -
    -
-> many_kets_solution a b =
->     --
->     -- List of columns of matrix X, which is
->     -- a solution of the matrix equation
->     --     A X = B
->     -- where
->     --     a is a list of columns of matrix A
->     --     b is a list of columns of matrix B
->     --
->     solve' p q emptyLists
->     where
->         (p, q) = unzip (reverse (many_kets_triangle a b))
->         emptyLists = [[] | k <- [1..(length (head q))]]
->         solve' a' b' x
->             | null a'  = x
->             | otherwise = solve' (tail a') (tail b')
->                                 [(f vk xk):xk  | (xk, vk) <- (zip x v)]
->             where
->                 f vk xk = (vk - (sum_product (tail u) xk))/(head u)
->                 u = head a'
->                 v = head b'
-
-
-
    -

    -

    -

    - - Matrix inversion - -

    - Function 'many_kets_solution' can be used to compute - inverse of matrix A by specializing matrix B to a unit - matrix I: -

    -
-	A X = I
-
    - It follows that matrix X is an inverse of A; that is X = A-1. -
    -
-> inverse :: (Scalar a, Fractional a) => [[a]] -> [[a]]
-> inverse a = many_kets_solution a (unit_matrix (length a))
->       --
->       -- List of columns of inverse of matrix A
->       -- where
->       --     a is list of columns of A
-
-
    -

    -

    -

    - - QR factorization of matrices - -

    - The process described above and implemented by - 'many_kets_triangle' function transforms the equation -

    -	A X = B
-
    - into another equation for the same matrix X -
    -	R X = S
-
    - where R is an upper triangular matrix. All operations - performed on matrices A and B during this process are linear, - and therefore we should be able to find a square matrix Q - that describes the entire process in one step. Indeed, assuming - that matrix A can be decomposed as a product of unknown matrix Q - and triangular matrix R and that Q-1 is an inverse of matrix Q - we can reach the last equation by following these steps: -
    -	A X       = B
-	(Q R) X   = B
-	Q-1 Q R X = Q-1 B
-	R X       = S
-
    - It follows that during this process a given matrix B - transforms to matrix S, as delivered by 'many_kets_triangle': -
    -	S = Q-1 B
-
    - from which the inverse of Q can be found: -
    -	Q-1 = S B-1
-
    - Having a freedom of choice of the right hand side matrix B - we can choose the unit matrix I in place of B, and therefore - simplify the definition of Q-1: -
    -	Q-1 = S,  if B is unit matrix
-
    - It follows that any non-singular matrix A can be decomposed - as a product of a matrix Q and a triangular matrix R - -
    -	A = Q R
-
    - where matrices Q-1 and R are delivered by "many_kets_triangle" - as a result of triangularization process of equation: -
    -	A X = I
-
    - The function below extracts a pair of matrices Q and R - from the answer provided by "many_kets_triangle". - During this process it inverts matrix Q-1 to Q. - This factorization will be used by a sequence of similarity - transformations to be defined in the next section. - -
    -
-> factors_QR :: (Fractional a, Scalar a) => [[a]] -> ([[a]],[[a]])
-> factors_QR a =
->       --
->       -- A pair of matrices (Q, R), such that
->       -- A = Q R
->       -- where
->       --     R is upper triangular matrix in row representation
->       --     (without redundant zeros)
->       --     Q is a transformation matrix in column representation
->       --     A is square matrix given as columns
->       --
->       (inverse (transposed q1),r)
->       where
->           (r, q1) = unzip (many_kets_triangle a (unit_matrix (length a)))
-
-
    - -

    -

    -

    - - Computation of the determinant - - - - -

    -
-> determinant :: (Fractional a, Scalar a) => [[a]] -> a
-> determinant a =
->    let (q,r) = factors_QR a
->    -- matrix Q is not normed so we have to respect the norms of its rows
->    in  product (map norm q) * product (map head r)
-
-
    - -Naive division-free computation of the determinant by expanding the first column. -It consumes n! multiplications. - -
    -
-> determinantNaive :: (Num a) => [[a]] -> a
-> determinantNaive [] = 1
-> determinantNaive m  =
->    sum (alternate
->       (zipWith (*) (map head m)
->           (map determinantNaive (removeEach (map tail m)))))
-
-
    - -Compute the determinant with about n^4 multiplications -without division according to the clow decomposition algorithm -of Mahajan and Vinay, and Berkowitz -as presented by Günter Rote: - -Division-Free Algorithms for the Determinant and the Pfaffian: -Algebraic and Combinatorial Approaches. - -
    -
-> determinantClow :: (Num a) => [[a]] -> a
-> determinantClow [] = 1
-> determinantClow m =
->    let lm = length m
->    in  parityFlip lm (last (newClow m
->           (nest (lm-1) (longerClow m)
->               (take lm (iterate (0:) [1])))))
-
-
    - -Compute the weights of all clow sequences -where the last clow is closed and a new one is started. - -
    -
-> newClow :: (Num a) => [[a]] -> [[a]] -> [a]
-> newClow a c =
->    scanl (-) 0
->          (sumVec (zipWith (zipWith (*)) (List.transpose a) c))
-
-
    - -Compute the weights of all clow sequences -where the last (open) clow is extended by a new arc. - -
    -
-> extendClow :: (Num a) => [[a]] -> [[a]] -> [[a]]
-> extendClow a c =
->    map (\ai -> sumVec (zipWith scaleVec ai c)) a
-
-
    - -Given the matrix of all weights of clows of length l -compute the weight matrix for all clows of length (l+1). -Take the result of 'newClow' as diagonal -and the result of 'extendClow' as lower triangle -of the weight matrix. - -
    -
-> longerClow :: (Num a) => [[a]] -> [[a]] -> [[a]]
-> longerClow a c =
->    let diagonal = newClow a c
->        triangle = extendClow a c
->    in  zipWith3 (\i t d -> take i t ++ [d]) [0 ..] triangle diagonal
-
-
    - -Auxiliary functions for the clow determinant. - -
    -
-> {- | Compositional power of a function,
->      i.e. apply the function n times to a value. -}
-> nest :: Int -> (a -> a) -> a -> a
-> nest 0 _ x = x
-> nest n f x = f (nest (n-1) f x)
->
-> {- successively select elements from xs and remove one in each result list -}
-> removeEach :: [a] -> [[a]]
-> removeEach xs =
->    zipWith (++) (List.inits xs) (tail (List.tails xs))
->
-> alternate :: (Num a) => [a] -> [a]
-> alternate = zipWith id (cycle [id, negate])
->
-> parityFlip :: Num a => Int -> a -> a
-> parityFlip n x = if even n then x else -x
->
-> {-| Weight a list of numbers by a scalar. -}
-> scaleVec :: (Num a) => a -> [a] -> [a]
-> scaleVec k = map (k*)
->
-> {-| Add corresponding numbers of two lists. -}
-> {- don't use zipWith because it clips to the shorter list -}
-> addVec :: (Num a) => [a] -> [a] -> [a]
-> addVec x [] = x
-> addVec [] y = y
-> addVec (x:xs) (y:ys) = x+y : addVec xs ys
->
-> {-| Add some lists. -}
-> sumVec :: (Num a) => [[a]] -> [a]
-> sumVec = foldl addVec []
-
-
    - - - -

    -

    -

    - - Similarity transformations and eigenvalues - -

    - Two n-square matrices A and B are called similar if there - exists a non-singular matrix S such that: -

    -	B = S-1 A S
-
    - - It can be proven that: -
      -
    • - Any two similar matrices have the same eigenvalues -
    • - Every n-square matrix A is similar to a triangular matrix - whose diagonal elements are the eigenvalues of A. -
    -

    - If matrix A can be transformed to a triangular or a diagonal - matrix Ak by a sequence of similarity transformations then - the eigenvalues of matrix A are the diagonal elements of Ak. - -

    - - Let's construct the sequence of matrices similar to A -

    -	A, A1, A2, A3...
-
    - by the following iterations - each of which factorizes a matrix - by applying the function 'factors_QR' and then forms a product - of the factors taken in the reverse order: -
    -	A = Q R          = Q (R Q) Q-1    = Q A1 Q-1        =
-	  = Q (Q1 R1) Q-1 = Q Q1 (R1 Q1) Q1-1 Q-1 = Q Q1 A2 Q1-1 Q-1 =
-	  = Q Q1 (Q2 R2) Q1-1 Q-1 = ...
-
-
    - We are hoping that after some number of iterations some matrix - Ak would become triangular and therefore its diagonal - elements could serve as eigenvalues of matrix A. As long as - a matrix has real eigenvalues only, this method should work well. - This applies to symmetric and hermitian matrices. It appears - that general complex matrices -- hermitian or not -- can also - be handled this way. Even more, this method also works for some - nonsymmetric real matrices, which have real eigenvalues only. -
    - The only type of matrices that cannot be treated by this algorithm - are real nonsymmetric matrices, whose some eigenvalues are complex. - There is no operation in the process that converts real elements - to complex ones, which could find their way into diagonal - positions of a triangular matrix. But a simple preconditioning - of a matrix -- described in the next section -- replaces - a real matrix by a complex one, whose eigenvalues are related - to the eigenvalues of the matrix being replaced. And this allows - us to apply the same method all across the board. -
    - It is worth noting that a process known in literature as QR - factorization is not uniquely defined and different algorithms - are employed for this. The algorithms using QR factorization - apply only to symmetric or hermitian matrices, and Q matrix - must be either orthogonal or unitary. -
    - But our transformation matrix Q is not orthogonal nor unitary, - although its first row is orthogonal to all other rows. In fact, - this factorization is only similar to QR factorization. We just - keep the same name to help identify a category of the methods - to which it belongs. -
    - The same factorization is used for tackling two major problems: - solving the systems of linear equations and finding the eigenvalues - of matrices. -
    - Below is the function 'similar_to', which makes a new matrix that is - similar to a given matrix by applying our similarity transformation. -
    - Function 'iterated_eigenvalues' applies this transformation n - times - storing diagonals of each new matrix as approximations of - eigenvalues. -
    - Function 'eigenvalues' follows the same process but reports the last - approximation only. -
    -
-
-> similar_to :: (Fractional a, Scalar a) => [[a]] -> [[a]]
-> similar_to a =
->       --
->       -- List of columns of matrix A1 similar to A
->       -- obtained by factoring A as Q R and then
->       -- forming the product A1 = R Q = (inverse Q) A Q
->       -- where
->       --     a is list of columns of A
->       --
->       triangle_matrix' r q
->       where
->           (q,r) = factors_QR a
-
-> iterated_eigenvalues a n
->       --
->       -- List of vectors representing
->       -- successive approximations of
->       -- eigenvalues of matrix A
->       -- where
->       --     a is a list of columns of A
->       --     n is a number of requested iterations
->       --
->       | n == 0 = []
->       | otherwise = (diagonals a)
->                     : iterated_eigenvalues (similar_to a) (n-1)
-
-> eigenvalues a n
->       --
->       -- Eigenvalues of matrix A
->       -- obtained by n similarity iterations
->       -- where
->       --     a are the columns of A
->       --
->       | n == 0    = diagonals a
->       | otherwise = eigenvalues (similar_to a) (n-1)
-
-
    -

    -

    -

    - - Preconditioning of real nonsymmetric matrices - -

    - As mentioned above, our QR-like factorization method works - well with almost all kind of matrices, but with the exception - of a class of real nonsymmetric matrices that have - complex eigenvalues. -

    - There is no mechanism in that method that would be able to - produce complex eigenvalues out of the real components of - this type of nonsymmetric matrices. Simple trivial replacement - of real components of a matrix by its complex counterparts - does not work because zero-valued imaginary components do - not contribute in any way to production of nontrivial - imaginary components during the factorization process. -
    - What we need is a trick that replaces real nonsymmetric matrix - by a nontrivial complex matrix in such a way that the results - of such replacements could be undone when the series of - similarity transformations finally produced the expected - effect in a form of a triangular matrix. -
    - The practical solution is surprisingly simple: - it's suffice to add any complex number, such as "i", to the - main diagonal of a matrix, and when triangularization is done - -- subtract it back from computed eigenvalues. - The explanation follows. -

    - Consider the eigenproblem for real and nonsymmetric matrix A. -

    -	A | x > = a | x >
-
    - Let us now define a new complex matrix B, such that: -
    -	B = A + alpha I
-	where
-	    I is a unit matrix and alpha is a complex scalar
-
    - It is obvious that matrices A and B commute; that is: -
    -	A B = B A
-
    - It can be proven that if two matrices commute then they - have the same eigenvectors. Therefore we can use vector - | x > of matrix A as an eigenvector of B: -
    -	B | x > = b | x >
-	B | x > = A | x > + alpha I | x >
-		= a | x > + alpha | x >
-		= (a + alpha) | x >
-
    - It follows that eigenvalues of B are related to the eigenvalues - of A by: -
    -	b = a + alpha
-
    - After eigenvalues of complex matrix B have been succesfully - computed, all what remains is to subtract "alpha" from them - all to obtain eigenvalues of A. And nothing has to be done - to eigenvectors of B - they are the same for A as well. - Simple and elegant! -

    - Below is an auxiliary function that adds a scalar to the - diagonal of a matrix: - -

    -
-> add_to_diagonal :: Num a => a -> [[a]] -> [[a]]
-> add_to_diagonal alpha a =
->       --
->       -- Add constant alpha to diagonal of matrix A
->       --
->       [f ai ni | (ai,ni) <- zip a [0..(length a -1)]]
->       where
->           f b k = p++[head q + alpha]++(tail q)
->               where
->                   (p,q) = splitAt k b
->
-
-
-
    -

    -

    -

    - - Examples of iterated eigenvalues - -

    - - - Here is an example of a symmetric real matrix with results - of application of function 'iterated_eigenvalues'. -

    -	| 7  -2  1 |
-	|-2  10 -2 |
-	| 1  -2  7 |
-
-	 [[7.0,     10.0,    7.0],
-	  [8.66667, 9.05752, 6.27582],
-	  [10.7928, 7.11006, 6.09718],
-	  [11.5513, 6.40499, 6.04367],
-	  [11.7889, 6.18968, 6.02142],
-	  [11.8943, 6.09506, 6.01068],
-	  [11.9468, 6.04788, 6.00534],
-	  [11.9733, 6.02405, 6.00267],
-	  [11.9866, 6.01206, 6.00134],
-	  [11.9933, 6.00604, 6.00067],
-	  [11.9966, 6.00302, 6.00034],
-	  [11.9983, 6.00151, 6.00017],
-	  [11.9992, 6.00076, 6.00008],
-	  [11.9996, 6.00038, 6.00004],
-	  [11.9998, 6.00019, 6.00002],
-	  [11.9999, 6.00010, 6.00001],
-	  [11.9999, 6.00005, 6.00001]]
-
-	  The true eigenvalues are:
-	  12, 6, 6
-
-
    - Here is an example of a hermitian matrix. (Eigenvalues of hermitian - matrices are real.) The algorithm works well and converges fast. -
    -	| 2   0     i|
-	[ 0   1   0  |
-	[ -i  0   2  |
-
-	[[2.8     :+ 0.0, 1.0 :+ 0.0, 1.2     :+ 0.0],
-	 [2.93979 :+ 0.0, 1.0 :+ 0.0, 1.06021 :+ 0.0],
-	 [2.97972 :+ 0.0, 1.0 :+ 0.0, 1.02028 :+ 0.0],
-	 [2.9932  :+ 0.0, 1.0 :+ 0.0, 1.0068  :+ 0.0],
-	 [2.99773 :+ 0.0, 1.0 :+ 0.0, 1.00227 :+ 0.0],
-	 [2.99924 :+ 0.0, 1.0 :+ 0.0, 1.00076 :+ 0.0],
-	 [2.99975 :+ 0.0, 1.0 :+ 0.0, 1.00025 :+ 0.0],
-	 [2.99992 :+ 0.0, 1.0 :+ 0.0, 1.00008 :+ 0.0],
-	 [2.99997 :+ 0.0, 1.0 :+ 0.0, 1.00003 :+ 0.0],
-	 [2.99999 :+ 0.0, 1.0 :+ 0.0, 1.00001 :+ 0.0],
-	 [3.0     :+ 0.0, 1.0 :+ 0.0, 1.0     :+ 0.0],
-	 [3.0     :+ 0.0, 1.0 :+ 0.0, 1.0     :+ 0.0],
-	 [3.0     :+ 0.0, 1.0 :+ 0.0, 1.0     :+ 0.0]]
-
-
    - Here is another example: this is a complex matrix and it is not - even hermitian. Yet, the algorithm still works, although its - fluctuates around true values. -
    -	| 2-i   0      i |
-	| 0     1+i  0   |
-	|   i   0    2-i |
-
-	[[2.0     :+ (-1.33333), 1.0 :+ 1.0, 2.0     :+ (-0.666667)],
-	 [1.89245 :+ (-1.57849), 1.0 :+ 1.0, 2.10755 :+ (-0.421509)],
-	 [1.81892 :+ (-1.80271), 1.0 :+ 1.0, 2.18108 :+ (-0.197289)],
-	 [1.84565 :+ (-1.99036), 1.0 :+ 1.0, 2.15435 :+ (-0.00964242)],
-	 [1.93958 :+ (-2.07773), 1.0 :+ 1.0, 2.06042 :+ 0.0777281],
-	 [2.0173  :+ (-2.06818), 1.0 :+ 1.0, 1.9827  :+ 0.0681793],
-	 [2.04357 :+ (-2.02437), 1.0 :+ 1.0, 1.95643 :+ 0.0243654],
-	 [2.03375 :+ (-1.99072), 1.0 :+ 1.0, 1.96625 :+ (-0.00928429)],
-	 [2.01245 :+ (-1.97875), 1.0 :+ 1.0, 1.98755 :+ (-0.0212528)],
-	 [1.99575 :+ (-1.98307), 1.0 :+ 1.0, 2.00425 :+ (-0.0169263)],
-	 [1.98938 :+ (-1.99359), 1.0 :+ 1.0, 2.01062 :+ (-0.00640583)],
-	 [1.99145 :+ (-2.00213), 1.0 :+ 1.0, 2.00855 :+ 0.00212504],
-	 [1.9968  :+ (-2.00535), 1.0 :+ 1.0, 2.0032  :+ 0.00535265],
-	 [2.00108 :+ (-2.00427), 1.0 :+ 1.0, 1.99892 :+ 0.0042723],
-	 [2.00268 :+ (-2.00159), 1.0 :+ 1.0, 1.99732 :+ 0.00158978],
-	 [2.00213 :+ (-1.99946), 1.0 :+ 1.0, 1.99787 :+ (-0.000541867)],
-	 [2.00079 :+ (-1.99866), 1.0 :+ 1.0, 1.9992  :+ (-0.00133514)],
-	 [1.99973 :+ (-1.99893), 1.0 :+ 1.0, 2.00027 :+ (-0.00106525)],
-	 [1.99933 :+ (-1.9996) , 1.0 :+ 1.0, 2.00067 :+ (-0.000397997)],
-	 [1.99947 :+ (-2.00013), 1.0 :+ 1.0, 2.00053 :+ 0.000134972]]
-
-	 The true eigenvalues are
-	 2 - 2i, 1 + i, 2
-
    - Some nonsymmetric real matrices have all real eigenvalues and - our algorithm still works for such cases. Here is one - such an example, which traditionally would have to be treated - by one of the Lanczos-like algorithms, specifically designed - for nonsymmetric real matrices. Evaluation of -
    - - iterated_eigenvalues [[2,1,1],[-2,1,3],[3,1,-1::Double]] 20 - -
    - gives the following results -
    -	[[3.0,     -0.70818,-0.291815],
-	 [3.06743, -3.41538, 2.34795],
-	 [3.02238, -1.60013, 0.577753],
-	 [3.00746, -2.25793, 1.25047],
-	 [3.00248, -1.88764, 0.885154],
-	 [3.00083, -2.06025, 1.05943],
-	 [3.00028, -1.97098, 0.970702],
-	 [3.00009, -2.0148,  1.01471],
-	 [3.00003, -1.99268, 0.992648],
-	 [3.00001, -2.00368, 1.00367],
-	 [3.0,     -1.99817, 0.998161],
-	 [3.0,     -2.00092, 1.00092],
-	 [3.0,     -1.99954, 0.99954],
-	 [3.0,     -2.00023, 1.00023],
-	 [3.0,     -1.99989, 0.999885],
-	 [3.0,     -2.00006, 1.00006],
-	 [3.0,     -1.99997, 0.999971],
-	 [3.0,     -2.00001, 1.00001],
-	 [3.0,     -1.99999, 0.999993],
-	 [3.0,     -2.0,     1.0]]
-
-	 The true eigenvalues are:
-	 3, -2, 1
-
    - Finally, here is a case of a nonsymmetric real matrix with - complex eigenvalues: -
    -	| 2 -3 |
-	| 1  0 |
-
    - The direct application of "iterated_eigenvalues" would - fail to produce expected eigenvalues: -
    -	1 + i sqrt(2) and 1 - i sqrt (2)
-
    - But if we first precondition the matrix by adding "i" to its diagonal: -
    -	| 2+i  -3|
-	| 1     i|
-
    - and then compute its iterated eigenvalues: -
    - - iterated_eigenvalues [[2:+1,1],[-3,0:+1]] 20 - -
    - then the method will succeed. Here are the results: -
    -
-	[[1.0     :+ 1.66667, 1.0     :+   0.333333 ],
-	[0.600936 :+ 2.34977, 1.39906 :+ (-0.349766)],
-	[0.998528 :+ 2.59355, 1.00147 :+ (-0.593555)],
-	[1.06991  :+ 2.413,   0.93009 :+ (-0.412998)],
-	[1.00021  :+ 2.38554, 0.99979 :+ (-0.385543)],
-	[0.988004 :+ 2.41407, 1.012   :+ (-0.414074)],
-	[0.999963 :+ 2.41919, 1.00004 :+ (-0.419191)],
-	[1.00206  :+ 2.41423, 0.99794 :+ (-0.414227)],
-	[1.00001  :+ 2.41336, 0.99999 :+ (-0.413361)],
-	[0.999647 :+ 2.41421, 1.00035 :+ (-0.414211)],
-	[0.999999 :+ 2.41436, 1.0     :+ (-0.41436) ],
-	[1.00006  :+ 2.41421, 0.99993 :+ (-0.414214)],
-	[1.0      :+ 2.41419, 1.0     :+ (-0.414188)],
-	[0.99999  :+ 2.41421, 1.00001 :+ (-0.414213)],
-	[1.0      :+ 2.41422, 1.0     :+ (-0.414218)],
-	[1.0      :+ 2.41421, 0.99999 :+ (-0.414213)],
-	[1.0      :+ 2.41421, 1.0     :+ (-0.414212)],
-	[1.0      :+ 2.41421, 1.0     :+ (-0.414213)],
-	[1.0      :+ 2.41421, 1.0     :+ (-0.414213)],
-	[1.0      :+ 2.41421, 1.0     :+ (-0.414213)]]
-
    - After subtracting "i" from the last result, we will get - what is expected. - -

    -

    -

    - - Eigenvectors for distinct eigenvalues - -

    - Assuming that eigenvalues of matrix A are already found - we may now attempt to find the corresponding aigenvectors - by solving the following homogeneous equation -

    -	(A - a I) | x > = 0
-
    - for each eigenvalue "a". The matrix -
    -	B = A - a I
-
    - is by definition singular, but in most cases it can be - triangularized by the familiar "factors_QR" procedure. -
    -	B | x > = Q R | x > = 0
-
    - It follows that the unknown eigenvector | x > is one of - the solutions of the homogeneous equation: - -
    -	R | x > = 0
-
    - where R is a singular, upper triangular matrix with at least one - zero on its diagonal. -
    - If | x > is a solution we seek, so is its scaled version - alpha | x >. Therefore we have some freedom of scaling choice. - Since this is a homogeneous equation, one of the components - of | x > can be freely chosen, while the remaining components - will depend on that choice. - - To solve the above, we will be working from the bottom up of - the matrix equation, as illustrated in the example below: -
    -	| 0     1     1     3     | | x1 |
-	| 0     1     1     2     | | x2 |      /\
-	| 0     0     2     4     | | x3 | = 0  ||
-	| 0     0     0     0     | | x4 |      ||
-
    - Recall that the diagonal elements of any triangular matrix - are its eigenvalues. - Our example matrix has three distinct eigenvalues: - 0, 1, 2. The eigenvalue 0 has degree of degeneration two. - Presence of degenerated eigenvalues complicates - the solution process. The complication arises when we have to - make our decision about how to solve the trivial scalar equations - with zero coefficients, such as -
    -	0 * x4 = 0
-
    - resulting from multiplication of the bottom row by vector | x >. - Here we have two choices: "x4" could be set to 0, or to any - nonzero number 1, say. By always choosing the "0" option - we might end up with the all-zero trivial vector -- which is - obviously not what we want. Persistent choice of the "1" option, - might lead to a conflict between some of the equations, such as - the equations one and four in our example. -

    - So the strategy is as follows. -

    - If there is at least one zero on the diagonal, find the topmost - row with zero on the diagonal and choose for it the solution "1". - Diagonal zeros in other rows would force the solution "0". - If the diagonal element is not zero than simply solve - an arithmetic equation that arises from the substitutions of - previously computed components of the eigenvector. Since certain - inaccuracies acumulate during QR factorization, set to zero all - very small elements of matrix R. -

    - By applying this strategy to our example we'll end up with the - eigenvector -

    -	< x | = [1, 0, 0, 0]
-
    - -

    - If the degree of degeneration of an eigenvalue of A is 1 then the - corresponding eigenvector is unique -- subject to scaling. - Otherwise an eigenvector found by this method is one of many - possible solutions, and any linear combination of such solutions - is also an eigenvector. This method is not able to find more than one - solution for degenerated eigenvalues. An alternative method, which - handles degenerated cases, will be described in the next section. -

    - The function below calculates eigenvectors corresponding to - distinct selected eigenvalues of any square matrix A, provided - that the singular matrix B = A - a I can still be factorized as Q R, - where R is an upper triangular matrix. - -

    -
-> eigenkets a u
->       --
->       -- List of eigenkets of a square matrix A
->       -- where
->       --     a is a list of columns of A
->       --     u is a list of eigenvalues of A
->       --     (This list does not need to be complete)
->       --
->       | null u        = []
->       | not (null x') = x':(eigenkets a (tail u))
->       | otherwise     = (eigenket_UT (reverse b) d []):(eigenkets a (tail u))
->       where
->           a'  = add_to_diagonal (-(head u)) a
->           x'  = unit_ket a' 0 (length a')
->           b   = snd (factors_QR a')
->           d   = discriminant [head bk | bk <- b] 1
->           discriminant u n
->               | null u    = []
->               | otherwise = x : (discriminant (tail u) m)
->               where
->                   (x, m)
->                       | (head u) == 0     = (n, 0)
->                       | otherwise         = (n, n)
->           eigenket_UT b d xs
->               | null b   = xs
->               | otherwise = eigenket_UT (tail b) (tail d) (x:xs)
->               where
->                   x = solve_row (head b) (head d) xs
->
->           solve_row u n x
->               | almostZero p = n
->               | otherwise    = q/p
->               where
->                   p = head u
->                   q
->                       | null x = 0
->                       | otherwise = -(sum_product (tail u) x)
->
->           unit_ket a' m n
->               | null a'              = []
->               | all (== 0) (head a') = unit_vector m n
->               | otherwise            = unit_ket (tail a') (m+1) n
-
-
    -

    -

    -

    - - Eigenvectors for degenerated eigenvalues - -

    - Few facts: -

      -
    • - Eigenvectors of a general matrix A, which does not have any - special symmetry, are not generally orthogonal. However, they - are orthogonal, or can be made orthogonal, to another set of - vectors that are eigenvectors of adjoint matrix A+; - that is the matrix obtained by complex conjugation and transposition - of matrix A. -
    • - Eigenvectors corresponding to nondegenerated eigenvalues of - hermitian or symmetric matrix are orthogonal. -
    • - Eigenvectors corresponding to degenerated eigenvalues are - in - general - neither orthogonal among themselves, nor orthogonal - to the remaining eigenvectors corresponding to other - eigenvalues. But since any linear combination of such degenerated - eigenvectors is also an eigenvector, we can orthogonalize - them by Gram-Schmidt orthogonalization procedure. -
    - Many practical applications deal solely with hermitian - or symmetric matrices, and for such cases the orthogonalization - is not only possible, but also desired for variety of reasons. -
    - But the method presented in the previous section is not able - to find more than one eigenvector corresponding to a degenerated - eigenvalue. For example, the symmetric matrix -
    -	    |  7  -2   1 |
-	A = | -2  10  -2 |
-	    |  1  -2   7 |
-
    - has two distinct eigenvalues: 12 and 6 -- the latter - being degenerated with degree of two. Two corresponding - eigenvectors are: -
    -	< x1 | = [1, -2, 1] -- for 12
-	< x2 | = [1,  1, 1] -- for 6
-
    - It happens that those vectors are orthogonal, but this is - just an accidental result. However, we are missing a third - distinct eigenvector. To find it we need another method. - One possibility is presented below and the explanation - follows. -
    -
-> eigenket' a alpha eps x' =
->       --
->       -- Eigenket of matrix A corresponding to eigenvalue alpha
->       -- where
->       --     a is a list of columns of matrix A
->       --     eps is a trial inaccuracy factor
->       --         artificially introduced to cope
->       --         with singularities of A - alpha I.
->       --         One might try eps = 0, 0.00001, 0.001, etc.
->       --     x' is a trial eigenvector
->       --
->       scaled [xk' - dk | (xk', dk) <- zip x' d]
->       where
->           b = add_to_diagonal (-alpha*(1+eps)) a
->           d = one_ket_solution b y
->           y = matrix_ket (transposed b) x'
-
-
    - Let us assume a trial vector | x' >, such that -
    -	| x' > = | x > + | d >
-	where
-	    | x > is an eigenvector we seek
-	    | d > is an error of our estimation of | x >
-
    - We first form a matrix B, such that: -
    -	B = A - alpha I
-
    - and multiply it by the trial vector | x' >, which - results in a vector | y > -
    -	B | x' > = |y >
-
    - On another hand: -
    -	B | x' > = B | x > + B | d > = B | d >
-	because
-	    B | x > = A | x > - alpha | x > = 0
-
    - Comparing both equations we end up with: -
    -	B | d > = | y >
-
    - that is: with the system of linear equations for unknown error | d >. - Finally, we subtract error | d > from our trial vector | x' > - to obtain the true eigenvector | x >. -

    - But there is some problem with this approach: matrix B is - by definition singular, and as such, it might be difficult - to handle. One of the two processes might fail, and their failures - relate to division by zero that might happen during either the - QR factorization, or the solution of the triangular system of equations. -

    - But if we do not insist that matrix B should be exactly singular, - but almost singular: -

    -	B = A - alpha (1 + eps) I
-
    - then this method might succeed. However, the resulting eigenvector - will be the approximation only, and we would have to experiment - a bit with different values of "eps" to extrapolate the true - eigenvector. -

    - The trial vector | x' > can be chosen randomly, although some - choices would still lead to singularity problems. Aside from - this, this method is quite versatile, because: -

      -
    • - Any random vector | x' > leads to the same eigenvector - for nondegenerated eigenvalues, -
    • - Different random vectors | x' >, chosen for degenerated - eigenvalues, produce -- in most cases -- distinct eigenvectors. - And this is what we want. If we need it, we can the always - orthogonalize those eigenvectors either internally (always - possible) or externally as well (possible only for hermitian - or symmetric matrices). -
    - It might be instructive to compute the eigenvectors for - the examples used in demonstration of computation of eigenvalues. - We'll leave to the reader, since this module is already too obese. -

    -

    -

    - - Auxiliary functions - -

    - The functions below are used in the main algorithms of - this module. But they can be also used for testing. For example, - the easiest way to test the usage of resources is to use easily - definable unit matrices and unit vectors, as in: - -

    -	one_ket_solution (unit_matrix n::[[Double]])
-			 (unit_vector 0 n::[Double])
-	where n = 20, etc.
-
-
-> unit_matrix :: Num a => Int -> [[a]]
-> unit_matrix m =
->       --
->       -- Unit square matrix of with dimensions m x m
->       --
->       [g 0 k | k <- [0..(m-1)]]
->       where
->       g i k
->           | i == m    = []
->           | i == k    = 1:(g (i+1) k)
->           | otherwise = 0:(g (i+1) k)
->
-
-> unit_vector :: Num a => Int -> Int -> [a]
-> unit_vector i m =
->       --
->       -- Unit vector of length m
->       -- with 1 at position i, zero otherwise
->       [g i k| k <- [0..(m-1)]]
->       where
->           g i k
->               | i == k    = 1
->               | otherwise = 0
-
-> diagonals :: [[a]] -> [a]
-> diagonals a =
->       --
->       -- Vector made of diagonal components
->       -- of square matrix a
->       --
->       diagonals' a 0
->       where
->           diagonals' a n
->               | null a = []
->               | otherwise = (head (drop n (head a)))
->                             :(diagonals' (tail a) (n+1))
-
-
-
    - -
    ------------------------------------------------------------------------------
---
--- Copyright:
---
---      (C) 1998 Numeric Quest Inc., All rights reserved
---
--- Email:
---
---      jans@numeric-quest.com
---
--- License:
---
---      GNU General Public License, GPL
---
------------------------------------------------------------------------------
-
    -
- - -