numeric-quest (empty) → 0.1
raw patch · 12 files changed
+5488/−0 lines, 12 filesdep +basedep +haskell98build-type:Customsetup-changed
Dependencies added: base, haskell98
Files
- Eigensystem.hs +173/−0
- EigensystemNum.hs +37/−0
- Fraction.hs +663/−0
- LinearAlgorithms.hs +379/−0
- Makefile +5/−0
- Orthogonals.lhs +1869/−0
- QuantumVector.lhs +1245/−0
- README +5/−0
- Roots.hs +110/−0
- Setup.lhs +3/−0
- Tensor.lhs +974/−0
- numeric-quest.cabal +25/−0
+ Eigensystem.hs view
@@ -0,0 +1,173 @@++------------------------------------------------------------------------------+-- 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:+--+-- This module extends the QuantumVector module by providing functions+-- to calculate eigenvalues and eigenvectors of Hermitian operators.+-- Such toolkit is of primary importance due to pervasiveness of+-- eigenproblems in Quantum Mechanics.+--+-- This module is organized in three layers:+--+-- 1. Interface to module QuantumVector, where all function signatures+-- are expressed in terms of linear operators, Dirac vectors and scalars.+--+-- Here the operators are defined directly via maps from input to+-- output vectors. In many cases it is much easier to define the operators+-- directly rather than to rely on their matrix representation.+--+-- 2. Conversion layer between operators and their matrix representation.+--+-- Sometimes it is more convenient to start with an underlying matrix+-- representation of an operator. There are also cases where a direct+-- 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 >+-- (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+-- direct approach from the upper layer could be easier and perhaps more+-- efficient in some cases.+--+-- 3. Implementation layer is stored in a separate module LinearAlgorithms,+-- where matrices are represented as lists of columns of scalars, and+-- vectors -- as lists of scalars.+--+-- This layer is completely independendent of the other two and can be+-- reused separately for applications other than those caring for the+-- QuantumVector module and its notation. It can also be reimplemented+-- via Haskell arrays, or perhaps by some other means, such as trees+-- of nodes relating square blocks of data to support paralleism.+--+-- See also bottom of the page for references and license.+-----------------------------------------------------------------------------++module Eigensystem (eigenvalues, adjoint) where+import Complex+import QuantumVector+import LinearAlgorithms (triangular, tridiagonal, triangular2)+import List (findIndex) ++----------------------------------------------------------------------------+-- Category: Eigensystem for QuantumVector+----------------------------------------------------------------------------++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+ where+ 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.+ = (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]++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+ where+ ms = [[ conjugate (toBra ei <> vj) | vj <- v] | ei <- es]+ v = [a ej | ej <- es]+++----------------------------------------------------------------------------+-- Category: Conversion from operators to matrices and vice versa+----------------------------------------------------------------------------++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+ where+ a u = case (findIndex (u == ) basis) of+ Just k -> compose (ms !! k) basis+ 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'+ = [[ei' <> vj | ei' <- e'] | vj <- v]+ where+ v = [a ej | ej <- basis]+ e' = [toBra ei | ei <- basis]++----------------------------------------------------------------------------+-- 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]+ ]++opA = operator basisA matrixA++basisA = map Ket [1..5::Int] -- or: map Ket "abcde", etc.+ +---------------------------------------------------------------------------+-- Copyright:+--+-- (C) 2001 Numeric Quest, All rights reserved+--+-- Email: jans@numeric-quest.com+--+-- http://www.numeric-quest.com +--+-- License:+--+-- GNU General Public License, GPL+--+---------------------------------------------------------------------------++
+ EigensystemNum.hs view
@@ -0,0 +1,37 @@+module EigensystemNum where++import Orthogonals+import List++mult :: Num a => [[a]] -> [[a]] -> [[a]]+mult x y = matrix_matrix x (transposed y)++matSqr :: Num a => [[a]] -> [[a]]+matSqr x = mult x x++powerIter :: (Fractional a, Ord a) => [[a]] -> [([[a]],[[a]])]+powerIter x = tail (iterate+ (\(_,z)->let s=normalize (matSqr z) in (s,(mult x s)))+ ([],x)+ )++normalize :: (Fractional a, Ord a) => [[a]] -> [[a]]+normalize x = map (map (/(matnorm1 x))) x++getGrowth :: (Fractional a, Ord a) => ([[a]],[[a]]) -> a+getGrowth (x,y) = uncurry (/) (maximumBy+ (\(_,xc) (_,xa) -> compare (abs xc) (abs xa))+ (concat (zipWith zip y x))+ )++specRadApprox :: (Fractional a, Ord a) => [[a]] -> [a]+specRadApprox = map getGrowth . powerIter++eigenValuesApprox :: (Scalar a, Fractional a) => [[a]] -> [[a]]+eigenValuesApprox = map diagonals . iterate similar_to++limit :: (Num a, Ord a) => a -> [a] -> a+limit tol (x0:x1:xs) = if abs (x1-x0) < tol * abs x0+ then x0+ else limit tol (x1:xs)+limit _ _ = error "Only infinite sequences are allowed"
+ Fraction.hs view
@@ -0,0 +1,663 @@+-- Module:+--+-- Fraction.hs+--+-- Language:+--+-- 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.+--+-- 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.+--+-- 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. +--+-- +-- Required:+--+-- Haskell Prelude +-- +-- Author:+--+-- Jan Skibinski, Numeric Quest Inc.+--+-- Date:+--+-- 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. +-------------------------------------------------------------------++module Fraction where++import Data.Ratio++infix 7 :-:++-------------------------------------------------------------------+-- Category: Basics+-------------------------------------------------------------------+ +data Fraction = Integer :-: Integer+ deriving (Eq)++num, den :: Fraction -> Integer+num (x:-:y) = x+den (x:-:y) = y+ +reduce :: Fraction -> Fraction+reduce (x:-:0)+ | x < 0 = (-1):-:0+ | otherwise = 1:-:0+reduce (x:-:y) =+ (u `quot` d) :-: (v `quot` d)+ where+ d = gcd u v+ (u,v)+ | y < 0 = (-x,-y)+ | otherwise = (x,y) + +(//) :: Integer -> Integer -> Fraction+x // y = reduce (x:-:y)++approx :: Fraction -> Fraction -> Fraction+approx eps (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)+ | 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 + (q,r) = quotRem n d+ (q',r') = quotRem n' d'+ (n'':-:d'') = simplest' d' r' d r++-------------------------------------------------------------------+-- Category: Instantiation of some Prelude classes+-------------------------------------------------------------------+ +instance Read Fraction where+ readsPrec p = + readParen (p > 7) (\r -> [(x//y,u) | (x,s) <- reads r,+ ("//",t) <- lex s,+ (y,u) <- reads t ])++instance Show Fraction where+ showsPrec p (x:-:y)+ | y == 1 = showsPrec p x+ | otherwise = showParen (p > 7) (shows x . showString "/" . shows y)++instance Ord Fraction where+ compare (x:-:y) (x':-:y') = compare (x*y') (x'*y)++instance Num Fraction where+ (x:-:y) + (x':-:y') = reduce ((x*y' + x'*y):-:(y*y'))+ (x:-:y) - (x':-:y') = reduce ((x*y' - x'*y):-:(y*y'))+ (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+ 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) + (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+ where+ 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 (x :-: y) = toRational (x % y)+++instance RealFrac Fraction where+ properFraction (x :-: y) = (fromInteger q, r :-: y)+ where (q,r) = quotRem x y++instance Enum Fraction where++ toEnum = fromIntegral+ fromEnum = truncate -- dubious+ enumFrom = numericEnumFrom+ enumFromTo = numericEnumFromTo+ enumFromThen = numericEnumFromThen+ enumFromThenTo = numericEnumFromThenTo++numericEnumFrom :: Real a => a -> [a]+numericEnumFromThen :: Real a => a -> a -> [a]+numericEnumFromTo :: Real a => a -> a -> [a]+numericEnumFromThenTo :: Real a => a -> a -> a -> [a]+--+-- Prelude does not export these, so here are the copies+++numericEnumFrom n = n : (numericEnumFrom $! (n+1))+numericEnumFromThen n m = iterate ((m-n)+) n+numericEnumFromTo n m = takeWhile (<= m) (numericEnumFrom n)+numericEnumFromThenTo n n' m = takeWhile p (numericEnumFromThen n n')+ where p | n' >= n = (<= m)+ | otherwise = (>= m)++------------------------------------------------------------------+-- Category: Conversion+-- from continued fraction to fraction and vice versa,+-- from Taylor series to continued 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++toCF :: Fraction -> CF+toCF (u:-:0) = [(u//0,0//1)]+toCF x =+ --+ -- 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 +++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+ ++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)+ ++fromFraction :: Fraction -> Double+fromFraction = fromRational . toRational+ +------------------------------------------------------------------+-- Category: Auxiliaries+------------------------------------------------------------------++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+ --+ 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++------------------------------------------------------------------+-- Category: Class Transcendental+--+-- 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+ 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 ()++-------------------------------------------------------------------+-- Everything below is the instantiation of class Transcendental+-- for type Fraction. See also modules Cofra and Numerus.+--+-- 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..]]+ a = 545140134+ b = 13591409+ c = 426880+ d = 10005+ +---------------------------------------------------------------------+-- Category: Trigonometry+---------------------------------------------------------------------+ + tan' eps 0 = 0+ tan' eps (u:-: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+ 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' eps x = (1 - p)/(1 + p)+ where+ 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)))++ + 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)+ +---------------------------------------------------------------------+-- 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))+ ++---------------------------------------------------------------------+-- 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))+ +---------------------------------------------------------------------+-- 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..]]+ + + cosh' eps x =+ --+ -- 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++ tanh' eps x =+ --+ -- 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)))+ + acosh' eps x+ --+ -- Inverse hyperbolic cosine+ --+ | x == 1//0 = 1//0+ | x < 1 = 1//0+ | otherwise = log' eps (x + (sqrt' eps (x^2 - 1)))+ +---------------------------------------------------------------------+-- 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 + 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)++ +---------------------------------------------------------------------+-- Category: IO+---------------------------------------------------------------------+ decimal n (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+ + + +---------------------------------------------------------------------------+-- References:+--+-- 1. Classical Gosper notes on continued fraction arithmetic:+-- http://www.inwap.com/pdp10/hbaker/hakmem/cf.html+-- 2. Pages on numerical constants represented as continued fractions:+-- 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+--------------------------------------------------------------------------++--------------------------------------------------------------------------++-- The following representation of continued fractions is used:+--+-- 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. +--+----------------------------------------------------------------------------+--+-- Copyright:+--+-- (C) 1998 Numeric Quest, All rights reserved+--+-- <jans@numeric-quest.com>+--+-- http://www.numeric-quest.com +--+-- License:+--+-- GNU General Public License, GPL+-- +-----------------------------------------------------------------------------
+ LinearAlgorithms.hs view
@@ -0,0 +1,379 @@++------------------------------------------------------------------------------+-- 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,+-- supporting computation of eigenvalues and eigenvectors of dense+-- matrices of small size. This module is to be utilized by module+-- Eigensystem, which redefines the eigenproblems in terms of+-- linear operators (maps) and abstract Dirac vectors.++-- 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+--+-- + 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+--+-- Other algoritms, such as solution of linear equations are, at this time,+-- imported from module Orthogonals. The latter also deals with triangulization,+-- so you can compare the results from two different approaches:+-- orthogonalization vs. Householder reduction used in this module.+-- In essence the former method is a bit faster but overflows for large+-- number of iterations since, for typing reasons - its algorithms+-- avoid the normalization of vectors.+-- For full documentation of this module, and for references and the license,+-- go to the bottom of the page.+----------------------------------------------------------------------------++module LinearAlgorithms (+ triangular,+ triangular2,+ tridiagonal,+ tridiagonal2,+ Scalar,) where++import Complex+import Orthogonals hiding (Scalar)++type Scalar = Complex Double++----------------------------------------------------------------------------+-- Category: Iterative triangularization+--+-- 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++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 ....+ --+ -- If matrix A is Hermitian then the result is close+ -- to a diagonal matrix for sufficiently large n.+ | n == 0 = a+ | otherwise = triangular (n - 1) a1+ where+ 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.+ = 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'+ ++----------------------------------------------------------------------------+-- 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 :: [[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+ 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))+++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+ 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))+++----------------------------------------------------------------------------+-- 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 :: 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+ --+ -- 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++ 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)+++reflection :: [Scalar] -> [Scalar] -> [Scalar] -> [Scalar]+reflection a e 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>.+ --+ -- 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++----------------------------------------------------------------------------+-- 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]+ ]++----------------------------------------------------------------------------+-- Module documentation+-- ====================++-- Representation of vectors, matrices and scalars:+-- ------------------------------------------------+-- We have chosen to follow the same scheme as used in module Orthogonals:+-- vectors are represented here as lists of scalars, while matrices --+-- as lists of scalar columns (vectors). But while scalars over there are+-- generic and cover a range of types, the scalars of this module are+-- implemented as Complex Double. Although all algorithms here+-- operate on complex matrices and complex vectors, they will work+-- on real matrices without modifications. If however, the performance+-- is a premium it will be a trivial exercise to customize all these+-- algorithms to real domain. Perhaps the most important change should+-- be then made to a true workhorse of this module, the function 'reflection',+-- in order to convert it to a real reflection of a vector in a hyperplane+-- whose normal is another vector.+--+-- Schur triangularization of any matrix:+-- --------------------------------------+-- 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+-- which preserves the eigenvalues. Here Q' stands for a Hermitian+-- conjugate of Q (adjoint, or Q-dagger).++-- Since the eigenvalues of a triangular matrix R are its diagonal+-- 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) >+--+-- In particular, when matrix A is Hermitian, then the matrix R+-- becomes diagonal, and the eigenvectors of R are its normalized+-- columns; that is, the unit vectors. It follows that the eigenvectors+-- of A are then the columns of matrix Q.+-- But when A is not Hermitian one must first find the eigenvectors+-- of a triangular matrix R before applying the above transformation.+-- Fortunately, it is easier to find eigenvectors of a triangular matrix+-- R than those of the square matrix A.+--+-- Implementation of Schur triangularization via series of QR factorizations:+-- --------------------------------------------------------------------------+-- 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) ...+-- 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+-- implementation uses a complex variety of such a 'reflection', described+-- in the reference [1]. The columnar reduction of the lower portion of+-- the matrix to zeros is also known under the name of Householder+-- reduction, or Householder transformation. This is, however, not the+-- only possible choice for elementary transformations; see for example+-- our module Orthogonals, where such transformations are perfomed via+-- Gram-Schmidt orthogonalization procedure instead.+--+-- The iterative functions 'triangular' and 'triangular2' attempt to+-- triangularize any complex matrix A by a series of similarity+-- transformation, known in literature as QR decomposition.+-- Function 'triangular' does not deliver the transformation Q but+-- only a transformed matrix A, which should be close to triangular+-- form after a sufficient number of iterations. Use this function+-- 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+--+-- Tridiagonalization of Hermitian matrices:+-- -----------------------------------------+-- While the above functions are iterative and require a bit of+-- experimentation with a count of iterations to figure out whether+-- the required accuracy has yet been achieved, the tridiagonalization+-- methods transform any matrix A to a tridiagonal form in a finite+-- number of elementary transformations.+--+-- However, our implementation is not generic because it performs+-- tridiagonalization only on Hermitian matrices. It uses the same+-- unitary 'reflection', as the triangularization does.+--+-- Why would you care for such tridiagonalization at all? Many world+-- class algorithms use it as a first step to precondition the original+-- matrix A for faster convergence and for better stability and accuracy.+-- Its cost is small in comparison to the overall cost incurred during+-- the iterative stage. What's more, the triangularization iteration+-- does preserve the shape of tridiagonal matrix at each step - bringing+-- it only closer to the diagonal shape. So the tridiagonalization+-- is a recommended option to be executed before the iterative+-- triangulariation.+--+-- 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+--+-- Elementary transformations:+-- ---------------------------+-- All the above algorithms heavily rely on the function 'reflection'+-- which defines a complex reflection transformation of a vector. One use+-- of this function is to perform a Householder reduction of a column-vector,+-- to zero out all of its components but one. For example, the unitary+-- transformation 'subsAnnihilator 0' annihilates all subdiagonals lying+-- below the main diagonal. Similarly, 'subsAnnihilator 1' would zero out+-- 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+--+----------------------------------------------------------------------------+-- 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 +---------------------------------------------------------------------------+--+-- Copyright:+--+-- (C) 2001 Numeric Quest, All rights reserved+--+-- Email: jans@numeric-quest.com+--+-- http://www.numeric-quest.com +--+-- License:+--+-- GNU General Public License, GPL+--+---------------------------------------------------------------------------+
+ Makefile view
@@ -0,0 +1,5 @@++html: Orthogonals.html QuantumVector.html Tensor.html++%.html: %.lhs+ ln -s $< $@
+ Orthogonals.lhs view
@@ -0,0 +1,1869 @@+<html>+<head>+<BASE HREF="http://www.numeric-quest.com/haskell/Orthogonals.html">++<title>+ Indexless linear algebra algorithms+</title>+</head>+<body>+<ul>+<center>+<h1>+ ***+</h1>+<h1>+ Indexless linear algebra algorithms+</h1>+<b>+<br>+ Orthogonalization, linear equations, eigenvalues and eigenvectors+<br>+ Literate Haskell module <i>Orthogonals.lhs</i>+</b>+<p>+ Jan Skibinski, <a href="http://www.numeric-quest.com/news/">+ Numeric Quest Inc.</a>, Huntsville, Ontario, Canada+<p>+ 1998.09.19, last modified 1998.12.28+</center>+<hr>+<p>+ 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.+<dd>+ 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.+<dd>+ 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.+<dd>+ 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.+<dd>+ 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.+<p>+ 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.+<dd>+ 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.+<p>+ 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.+<dd>+ 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.++</i>+<p>+<hr>+<p>+<b>+ Contents+</b>+<p>+<ul>+<p><li>+ Notation+<p><li>+ Scalar products and vector normalization+<ul>+<li><b>+ bra_ket</b>, scalar product+<li><b>+ sum_product</b>, a cousin of bra_ket+<li><b>+ norm</b>, vector norm+<li><b>+ normalized</b>, vector normalized to one+</ul>+<p><li>+ Transposition and adjoining of matrices+<ul>+<li><b>+ transposed</b>, transposed matrix+<li><b>+ adjoint</b>, transposed and conjugated matrix+</ul>+<p><li>+ Products involving matrices+<ul>+<li><b>+ matrix_matrix</b>, product of two matrices as list of rows+<li><b>+ matrix_matrix'</b>, product of two matrices as list of columns+<li><b>+ triangle_matrix'</b>, upper triangular matrix times square matrix+<li><b>+ matrix_ket</b>, matrix times ket vector+<li><b>+ bra_matrix</b>, bra vector times matrix+<li><b>+ bra_matrix_ket</b>, matrix multiplied on both sides by vectors+<li><b>+ scalar_matrix</b>, scalar times matrix+</ul>+<p><li>+ Orthogonalization process+<ul>+<li><b>+ orthogonals</b>, set of orthogonal vectors+<li><b>+ gram_schmidt</b>, vector perpendicular to a hyperplane+</ul>++<p><li>+ Solutions of linear equations by orthogonalization+<ul>+<li><b>+ one_ket_triangle</b>, triangularization of one vector equation+<li><b>+ one_ket_solution</b>, solution for one unknown vector+<li><b>+ many_kets_triangle</b>, triangularization of several vector equations+<li><b>+ many_kets_solution</b>, solution for several unknown vectors+</ul>+<p><li>+ Matrix inversion+<ul>+<li><b>+ inverse</b>, inverse of a matrix+</ul>+<p><li>+ QR factorization of matrices provided by "many_kets_triangle"+<ul>+<li><b>+ factors_QR</b>, QR alike factorization of matrices+<li><b>+ determinant</b>, computation of the determinant based on the QR factorization+</ul>+<p><li>+ Similarity transformations and eigenvalues+<ul>+<li><b>+ similar_to</b>, matrix obtained by similarity transformation+<li><b>+ iterated_eigenvalues</b>, list of approximations of eigenvalues+<li><b>+ eigenvalues</b>, final approximation of eigenvalues+</ul>+<p><li>+ Preconditioning of real nonsymmetric matrices+<ul>+<li><b>+ add_to_diagonal</b>, simple preconditioning method+</ul>+<p><li>+ Examples of iterated eigenvalues+<ul>+<li>+ Symmetric real matrix+<li>+ Hermitian complex matrix+<li>+ General complex matrix+<li>+ Nonsymmetric real matrix with real eigenvalues+<li>+ Nonsymmetric real matrix with complex eigenvalues+</ul>+<p><li>+ Eigenvectors for distinct eigenvalues+<ul>+<li><b>+ eigenkets</b>, eigenvectors for distinct eigenvalues+</ul>+<p><li>+ Eigenvectors for degenerated eigenvalues+<ul>+<li><b>+ eigenket'</b>, eigenvector based on a trial vector+</ul>++<p><li>+ Auxiliary functions+<ul>+<li><b>+ unit_matrix</b>, a unit matrix with 1's on a diagonal+<li><b>+ unit_vector</b>, a vector with one non-zero componenet+<li><b>+ diagonals</b>, vector made of a matrix diagonal+</ul>+</ul>++<p>+<hr>+<p>+<b>+ Notation+</b>+<p>+ 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:+<ul>+<p><li>+ Bra vector x, written < x |, is represented by one-row matrix+<p><li> Ket vector y, written | y >, is represented by one-column matrix+</ul>+<p>+ Bra vectors can be obtained from ket vectors by transposition+ and conjugation of their components. Conjugation is only+ important for complex vectors.+<p>+ Scalar product of two vectors | x > and | y > is written+ as+<pre>+ < x | y >+</pre>+ 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+<pre>+ < x | x >+</pre>+ 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.+<p>+<hr>+<p>+<pre>++> module Orthogonals where+> import Complex+> import Ratio+> import qualified List++</pre>+<b>+ Scalar product and vector normalization+</b>+<p>+ The scalar product "bra_ket" is a basis of many algorithms+ presented here.+++<pre>++> 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)++</pre>++ 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.+<dd>+ 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.+<pre>++> 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+++</pre>++ 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".+<pre>++> 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)++</pre>+ 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.+<pre>++> 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]++</pre>+<hr>+<p>+<b>+ Transposition and adjoining of matrices+</b>+<p>+ Matrices are represented here by lists of lists.+ Function "transposed" converts from row-wise to column-wise+ representation, or vice versa.+<dd>+ When transposition is combined with complex conjugation+ the resulting matrix is called "adjoint".+<p>+ A square matrix is called symmetric if it is equal to its transpose+<pre>+ A = A<sup>T</sup>+</pre>+ It is called Hermitian, or self-adjoint, if it equals to+ its adjoint+<pre>+ A = A<sup>+</sup>++> 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])++</pre>+<p>+<hr>+<p>+<b>+ Linear combination and sum of two matrices+</b>+<p>+ One can form a linear combination of two matrices, such+ as+<pre>+ C = alpha A + beta B+ where+ alpha and beta are scalars+</pre>+ 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:+<pre>+ f alpha beta a b = alpha*a + beta*b+</pre>+ For a straightforward addition of two matrices this auxiliary+ function is simply "(+)".+<pre>++> 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 (+)++</pre>++<p>+<hr>+<p>+<b>+ Products involving matrices+</b>+<p>+ 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.+<dd>+ The following definitions are somehow arbitrary, since+ one can choose alternative interpretations of lists+ representing matrices.+<p>+<b>+ C = A B+</b>+<p>+ 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.+<dd>+ 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.++<pre>++> 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))++++</pre>+<b>+ | u > = A | v >+</b>+<p>+ Product of a matrix and a ket-vector is another+ ket vector. The following implementation assumes+ that list "a" represents rows of matrix A.+<pre>++> matrix_ket :: Num a => [[a]] -> [a] -> [a]+> matrix_ket a v = [sum_product ai v| ai <- a]++</pre>+<b>+ < u | = < v | A+</b>+<p>+ 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".+<pre>++> bra_matrix :: (Scalar a, Num a) => [a] -> [[a]] -> [a]+> bra_matrix v a = [bra_ket v ai | ai <- a]++</pre>+<b>+ alpha = < u | A | v >+</b>+<p>+ This kind of product results in a scalar and is often+ used to define elements of a new matrix, such as+<pre>+ B[i,j] = < ei | A | ej >+</pre>+ The implementation below assumes that list "a" represents+ rows of matrix A.+<pre>++> bra_matrix_ket :: (Scalar a, Num a) => [a] -> [[a]] -> [a] -> a+> bra_matrix_ket u a v =+> bra_ket u (matrix_ket a v)++</pre>+<b>+ B = alpha A+</b>+<p>+ Below is a function which multiplies matrix by a scalar:+<pre>++> scalar_matrix :: Num a => a -> [[a]] -> [[a]]+> scalar_matrix alpha a =+> [[alpha*aij| aij <- ai] | ai<-a]++</pre>+<p>+<hr>+<p>+<b>+ Orthogonalization process+</b>+<p>++ Gram-Schmidt orthogonalization procedure is used here+ for calculation of sets of mutually orthogonal vectors.+<dd>+ 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:+<ul>+<li>+ None of the algorithms presented here needs this -- quite+ costly -- normalization.+<li>+ 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.+</ul>+<p>+ 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.+++<pre>++> 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))++</pre>+<p>+<hr>+<p>+<b>+ Solutions of linear equations by orthogonalization+</b>+<p>+ A matrix equation for unknown vector | x >+<pre>+ A | x > = | b >+</pre>+ can be rewritten as+<pre>+ x1 | 1 > + x2 | 2 > + x3 | 3 > + ... + xn | n > = | b > (7.1)+ where+ | 1 >, | 2 >... represent columns of the matrix A+</pre>+ 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:+<pre>+ < ck | 1 > = 0, for k = 1...m, where m = n - 1+</pre>+ 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+<pre>+ < 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 >+</pre>+ But the above is nothing more than a new matrix equation+<pre>+ A' | x' > = | b' >+ or++ x2 | 2'> + x3 | 3'> .... + xn | n'> = | b'>+ where+ primed vectors | 2' >, etc. are the columns of the new+ matrix A'.+</pre>+ with the problem dimension reduced by one.++<dd>+ 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:++<pre>+ x1 | 1 > + x2 | 2 > + x3 | 3 > + x4 | 4 > = | b >+ x2 | 2'> + x3 | 3'> + x4 | 4'> = | b'>+ x3 | 3''> + x4 | 4''> = | b''>+ x4 | 4'''> = | b'''>+</pre>+ 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:++<pre>+ | p11 p12 p13 p14 | | x1 | = | q1 |+ | 0 p22 p23 p24 | | x2 | = | q2 |+ | 0 0 p33 p34 | | x3 | = | q3 |+ | 0 0 0 p44 | | x4 | = | q4 |+</pre>+ In effect, we have triangularized our original matrix A.+ Below is a function that does that for any problem size:+<pre>++> 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++</pre>+ The triangular system of equations can be easily solved by+ successive substitutions - starting with the last equation.++<pre>++> 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++</pre>+ 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:+<pre>++> 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++</pre>+ 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.+<pre>++> 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'+++</pre>+<p>+<hr>+<p>+<b>+ Matrix inversion+</b>+<p>+ Function 'many_kets_solution' can be used to compute+ inverse of matrix A by specializing matrix B to a unit+ matrix I:+<pre>++ A X = I+</pre>+ It follows that matrix X is an inverse of A; that is X = A<sup>-1</sup>.+<pre>++> 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++</pre>+<p>+<hr>+<p>+<b>+ QR factorization of matrices+</b>+<p>+ The process described above and implemented by+ 'many_kets_triangle' function transforms the equation+<pre>+ A X = B+</pre>+ into another equation for the same matrix X+<pre>+ R X = S+</pre>+ 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<sup>-1</sup> is an inverse of matrix Q+ we can reach the last equation by following these steps:+<pre>+ A X = B+ (Q R) X = B+ Q<sup>-1</sup> Q R X = Q<sup>-1</sup> B+ R X = S+</pre>+ It follows that during this process a given matrix B+ transforms to matrix S, as delivered by 'many_kets_triangle':+<pre>+ S = Q<sup>-1</sup> B+</pre>+ from which the inverse of Q can be found:+<pre>+ Q<sup>-1</sup> = S B<sup>-1</sup>+</pre>+ 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<sup>-1</sup>:+<pre>+ Q<sup>-1</sup> = S, if B is unit matrix+</pre>+ It follows that any non-singular matrix A can be decomposed+ as a product of a matrix Q and a triangular matrix R++<pre>+ A = Q R+</pre>+ where matrices Q<sup>-1</sup> and R are delivered by "many_kets_triangle"+ as a result of triangularization process of equation:+<pre>+ A X = I+</pre>+ 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<sup>-1</sup> to Q.+ This factorization will be used by a sequence of similarity+ transformations to be defined in the next section.++<pre>++> 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)))++</pre>++<p>+<hr>+<p>+<b>+ Computation of the determinant+</b>++<!-- added by Henning Thielemann -->++<pre>++> 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)++</pre>++Naive division-free computation of the determinant by expanding the first column.+It consumes n! multiplications.++<pre>++> determinantNaive :: (Num a) => [[a]] -> a+> determinantNaive [] = 1+> determinantNaive m =+> sum (alternate+> (zipWith (*) (map head m)+> (map determinantNaive (removeEach (map tail m)))))++</pre>++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:+<a href="http://page.inf.fu-berlin.de/~rote/Papers/pdf/Division-free+algorithms.pdf">+Division-Free Algorithms for the Determinant and the Pfaffian:+Algebraic and Combinatorial Approaches</a>.++<pre>++> 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])))))++</pre>++Compute the weights of all clow sequences+where the last clow is closed and a new one is started.++<pre>++> newClow :: (Num a) => [[a]] -> [[a]] -> [a]+> newClow a c =+> scanl (-) 0+> (sumVec (zipWith (zipWith (*)) (List.transpose a) c))++</pre>++Compute the weights of all clow sequences+where the last (open) clow is extended by a new arc.++<pre>++> extendClow :: (Num a) => [[a]] -> [[a]] -> [[a]]+> extendClow a c =+> map (\ai -> sumVec (zipWith scaleVec ai c)) a++</pre>++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.++<pre>++> 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++</pre>++Auxiliary functions for the clow determinant.++<pre>++> {- | 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 []++</pre>++++<p>+<hr>+<p>+<b>+ Similarity transformations and eigenvalues+</b>+<p>+ Two n-square matrices A and B are called similar if there+ exists a non-singular matrix S such that:+<pre>+ B = S<sup>-1</sup> A S+</pre>++ It can be proven that:+<ul>+<li>+ Any two similar matrices have the same eigenvalues+<li>+ Every n-square matrix A is similar to a triangular matrix+ whose diagonal elements are the eigenvalues of A.+</ul>+<p>+ 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.++<p>++ Let's construct the sequence of matrices similar to A+<pre>+ A, A1, A2, A3...+</pre>+ 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:+<pre>+ A = Q R = Q (R Q) Q<sup>-1</sup> = Q A1 Q<sup>-1</sup> =+ = Q (Q1 R1) Q<sup>-1</sup> = Q Q1 (R1 Q1) Q1<sup>-1</sup> Q<sup>-1</sup> = Q Q1 A2 Q1<sup>-1</sup> Q<sup>-1</sup> =+ = Q Q1 (Q2 R2) Q1<sup>-1</sup> Q<sup>-1</sup> = ...++</pre>+ 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.+<dd>+ 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.+<dd>+ 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.+<dd>+ 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.+<dd>+ The same factorization is used for tackling two major problems:+ solving the systems of linear equations and finding the eigenvalues+ of matrices.+<dd>+ Below is the function 'similar_to', which makes a new matrix that is+ similar to a given matrix by applying our similarity transformation.+<dd>+ Function 'iterated_eigenvalues' applies this transformation n+ times - storing diagonals of each new matrix as approximations of+ eigenvalues.+<dd>+ Function 'eigenvalues' follows the same process but reports the last+ approximation only.+<pre>+++> 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)++</pre>+<p>+<hr>+<p>+<b>+ Preconditioning of real nonsymmetric matrices+</b>+<p>+ 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.+<dd>+ 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.+<dd>+ 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.+<dd>+ 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.+<p>+ Consider the eigenproblem for real and nonsymmetric matrix A.+<pre>+ A | x > = a | x >+</pre>+ Let us now define a new complex matrix B, such that:+<pre>+ B = A + alpha I+ where+ I is a unit matrix and alpha is a complex scalar+</pre>+ It is obvious that matrices A and B commute; that is:+<pre>+ A B = B A+</pre>+ 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:+<pre>+ B | x > = b | x >+ B | x > = A | x > + alpha I | x >+ = a | x > + alpha | x >+ = (a + alpha) | x >+</pre>+ It follows that eigenvalues of B are related to the eigenvalues+ of A by:+<pre>+ b = a + alpha+</pre>+ 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!+<p>+ Below is an auxiliary function that adds a scalar to the+ diagonal of a matrix:++<pre>++> 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+>+++</pre>+<p>+<hr>+<p>+<b>+ Examples of iterated eigenvalues+</b>+<p>+++ Here is an example of a symmetric real matrix with results+ of application of function 'iterated_eigenvalues'.+<pre>+ | 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++</pre>+ Here is an example of a hermitian matrix. (Eigenvalues of hermitian+ matrices are real.) The algorithm works well and converges fast.+<pre>+ | 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]]++</pre>+ 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.+<pre>+ | 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+</pre>+ 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+<br>+<i>+ iterated_eigenvalues [[2,1,1],[-2,1,3],[3,1,-1::Double]] 20+</i>+<br>+ gives the following results+<pre>+ [[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+</pre>+ Finally, here is a case of a nonsymmetric real matrix with+ complex eigenvalues:+<pre>+ | 2 -3 |+ | 1 0 |+</pre>+ The direct application of "iterated_eigenvalues" would+ fail to produce expected eigenvalues:+<pre>+ 1 + i sqrt(2) and 1 - i sqrt (2)+</pre>+ But if we first precondition the matrix by adding "i" to its diagonal:+<pre>+ | 2+i -3|+ | 1 i|+</pre>+ and then compute its iterated eigenvalues:+<br>+<i>+ iterated_eigenvalues [[2:+1,1],[-3,0:+1]] 20+</i>+<br>+ then the method will succeed. Here are the results:+<pre>++ [[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)]]+</pre>+ After subtracting "i" from the last result, we will get+ what is expected.++<p>+<hr>+<p>+<b>+ Eigenvectors for distinct eigenvalues+</b>+<p>+ Assuming that eigenvalues of matrix A are already found+ we may now attempt to find the corresponding aigenvectors+ by solving the following homogeneous equation+<pre>+ (A - a I) | x > = 0+</pre>+ for each eigenvalue "a". The matrix+<pre>+ B = A - a I+</pre>+ is by definition singular, but in most cases it can be+ triangularized by the familiar "factors_QR" procedure.+<pre>+ B | x > = Q R | x > = 0+</pre>+ It follows that the unknown eigenvector | x > is one of+ the solutions of the homogeneous equation:++<pre>+ R | x > = 0+</pre>+ where R is a singular, upper triangular matrix with at least one+ zero on its diagonal.+<dd>+ 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.+</pre>+ To solve the above, we will be working from the bottom up of+ the matrix equation, as illustrated in the example below:+<pre>+ | 0 1 1 3 | | x1 |+ | 0 1 1 2 | | x2 | /\+ | 0 0 2 4 | | x3 | = 0 ||+ | 0 0 0 0 | | x4 | ||+</pre>+ 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+<pre>+ 0 * x4 = 0+</pre>+ 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.+<p>+ So the strategy is as follows.+<p>+ 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.+<p>+ By applying this strategy to our example we'll end up with the+ eigenvector+<pre>+ < x | = [1, 0, 0, 0]+</pre>++<p>+ 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.+<p>+ 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.++<pre>++> 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++</pre>+<p>+<hr>+<p>+<b>+ Eigenvectors for degenerated eigenvalues+</b>+<p>+ Few facts:+<ul>+<li>+ 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<sup>+</sup>;+ that is the matrix obtained by complex conjugation and transposition+ of matrix A.+<li>+ Eigenvectors corresponding to nondegenerated eigenvalues of+ hermitian or symmetric matrix are orthogonal.+<li>+ 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.+</ul>+ 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.+<dd>+ 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+<pre>+ | 7 -2 1 |+ A = | -2 10 -2 |+ | 1 -2 7 |+</pre>+ has two distinct eigenvalues: 12 and 6 -- the latter+ being degenerated with degree of two. Two corresponding+ eigenvectors are:+<pre>+ < x1 | = [1, -2, 1] -- for 12+ < x2 | = [1, 1, 1] -- for 6+</pre>+ 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.+<pre>++> 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'++</pre>+ Let us assume a trial vector | x' >, such that+<pre>+ | x' > = | x > + | d >+ where+ | x > is an eigenvector we seek+ | d > is an error of our estimation of | x >+</pre>+ We first form a matrix B, such that:+<pre>+ B = A - alpha I+</pre>+ and multiply it by the trial vector | x' >, which+ results in a vector | y >+<pre>+ B | x' > = |y >+</pre>+ On another hand:+<pre>+ B | x' > = B | x > + B | d > = B | d >+ because+ B | x > = A | x > - alpha | x > = 0+</pre>+ Comparing both equations we end up with:+<pre>+ B | d > = | y >+</pre>+ 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 >.+<p>+ 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.+<p>+ But if we do not insist that matrix B should be exactly singular,+ but almost singular:+<pre>+ B = A - alpha (1 + eps) I+</pre>+ 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.+<p>+ 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:+<ul>+<li>+ Any random vector | x' > leads to the same eigenvector+ for nondegenerated eigenvalues,+<li>+ 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).+</ul>+ 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.+<p>+<hr>+<p>+<b>+ Auxiliary functions+</b>+<p>+ 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:++<pre>+ 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))+++</pre>++<pre>+-----------------------------------------------------------------------------+--+-- Copyright:+--+-- (C) 1998 Numeric Quest Inc., All rights reserved+--+-- Email:+--+-- jans@numeric-quest.com+--+-- License:+--+-- GNU General Public License, GPL+--+-----------------------------------------------------------------------------+</pre>+</ul>+</body>++<SCRIPT language="Javascript">+<!--++// FILE ARCHIVED ON 20010628005806 AND RETRIEVED FROM THE+// INTERNET ARCHIVE ON 20030626101500.+// JAVASCRIPT APPENDED BY WAYBACK MACHINE, COPYRIGHT INTERNET ARCHIVE.+// ALL OTHER CONTENT MAY ALSO BE PROTECTED BY COPYRIGHT (17 U.S.C.+// SECTION 108(a)(3)).++ var sWayBackCGI = "http://web.archive.org/web/20010628005806/";++ function xLateUrl(aCollection, sProp) {+ var i = 0;+ for(i = 0; i < aCollection.length; i++)+ if (aCollection[i][sProp].indexOf("mailto:") == -1 &&+ aCollection[i][sProp].indexOf("javascript:") == -1)+ aCollection[i][sProp] = sWayBackCGI + aCollection[i][sProp];+ }++ if (document.links) xLateUrl(document.links, "href");+ if (document.images) xLateUrl(document.images, "src");+ if (document.embeds) xLateUrl(document.embeds, "src");++ if (document.body && document.body.background)+ document.body.background = sWayBackCGI + document.body.background;++//-->++</SCRIPT>+</html>
+ QuantumVector.lhs view
@@ -0,0 +1,1245 @@+<html>+<head>+<BASE HREF="http://www.numeric-quest.com/haskell/QuantumVector.html">++<title>+ Quantum vector+</title>+</head>+<body>+<center>+<h1>+ ***+</h1>+<h1>+ Quantum vector+</h1>+<p>+<b>++ Jan Skibinski, <a href=http://www.numeric-quest.com/news/>+ Numeric Quest Inc.</a>, Huntsville, Ontario, Canada+<br>+ Literate Haskell module <em>QuantumVector.lhs</em>+<p>+ Initialized: 2000-05-31, last modified: 2000-06-10+</b>+</center>++<blockquote>+<em>+<p>+<hr>+<p>+ This is our attempt to model the abstract Dirac's formalism+ of Quantum Mechanics in Haskell. Although we have been+ developing quantum mechanical applications and examples for some time [2], the+ machinery used there is tightly coupled to a concrete+ representation of states and observables by complex vectors+ and matrices. implemented mainly as Haskell lazy lists.+<p>+ However, the Dirac's formalism in Hilbert space is much more+ abstract than that, and many problems of Quantum Mechanics can be+ solved without referring to any+ particular matrix representation, but using certain generic properties+ of operators, such as their commutative relations instead.+ Haskell seems to be well suited for such abstract tasks,+ even in its current form that does not support any+ of the abstract notions of computer algebra as yet.+ This has been already recognized by Jerzy Karczmarczuk [1],+ where he proposes a very interesting representation of Hilbert+ space and illustrates it by several powerful examples.+ But the task is not trivial and far from being complete.+ Quantum Mechanics presents many challenges to any formalism+ and only by careful examination of many of its facets+ and alternative approaches, a consistent model of+ Dirac's formalism can be developed for Haskell. Hoping to+ help with solving this problem, we present here a computing+ abstract, which is quite different from that of [1].+<p>+ We recognize a quantum state as an abstract vector | x >,+ which can be represented in one of many possible bases -- similar+ to many alternative representations of a 3D vector in rotated systems+ of coordinates. A choice of a particular basis is controlled+ by a generic type variable, which can be any Haskell object+ -- providing that it supports a notion of equality and ordering.+ A state which is composed of many quantum subsystems, not+ necessarily of the same type, can be represented in a vector+ space considered to be a tensor product of the subspaces.++<p>+ With this abstract notion we proceed with Haskell definition of two+ vector spaces: Ket and its dual Bra. We demonstrate+ that both are properly defined according to the abstract+ mathematical definition of vector spaces. We then introduce inner+ product and show that our Bra and Ket can be indeed+ considered the vector spaces with inner product. Multitude+ of examples is attached in the description. To verify+ the abstract machinery developed here we also provide the basic library+ module <a href="http://www.numeric-quest.com/haskell/Momenta.html">+ Momenta</a> -- a non-trivial example designed to compute Clebsch-Gordan coefficients+ of a transformation from one basis of angular momenta to another.+<p>+ Section 6 is a rehash of known definitions of linear operators+ with the emphasis on both Dirac and Haskell notations and on+ Haskell examples. The formalism developed here centers around+ two operations: a scalar product of two vectors, <b>x <> y</b>,+ and a closure operation, <b>a >< x</b>, which can be considered+ an application of a quantum operator <b>a</b> to a vector <b>x</b>.+ At this stage our formalism applies only to discrete cases, but+ we hope to generalize it on true Hilbert space as well.+</em>+<p>+<hr>+<p>+<b>+ Contents+</b>+<ul>+<li>+ 1. Infix operators+<li>+ 2. Vector space+<li>+ 3. Ket vector space+<li>+ 4. Bra vector space+<li>+ 5. Bra and Ket spaces as inner product spaces+<li>+ 6. Linear operators+<ul>+<li> 6.1. Operator notation+<li>+ 6.2. Renaming the representation+<li>+ 6.3. Closure formula, or identity operator+<li>+ 6.4. Changing the representation+<li>+ 6.5. Implementation of the operator equation A | x > = | y >+<li>+ 6.6. Inverse operator+<li>+ 6.7. Matrix representation of an operator+<li>+ 6.8. Adjoint operator+<li>+ 6.9. Unitary operator+<li>+ 6.10. Hermitian operator+</ul>+<li>+ 7. Showing kets and bras+<li>+ 8. Data Tuple for tensor products+<li>+ 9. References+<li>+ 10. Copyright and license++</ul>++<p>+<hr>+<p>+<b>+ 1. Infix operators+</b>+<p>+ Haskell requires that fixities of infix operators are defined+ at the top of the module. So here they are. They are+ to be explained later.++</b>+<pre>++> module QuantumVector where+> import Complex -- our Scalar is Complex Double+> import Fraction hiding (reduce) -- to bypass enum bug in Ratio and for better pretty printing+> import List (nub)++> infixl 7 *> -- tensor product of two kets+> infixl 7 <* -- tensor product of two bras++> -- scalar-ket multiplication+> infix 6 |>+> -- scalar-bra multiplication+> infix 6 <|+++> infixl 5 +> -- sum of two kets+> infixl 5 <+ -- sum of two bras+++> infix 4 <> -- inner product+> infix 5 >< -- closure++</pre>+<p>+<hr>+<p>+<b>+ 2. Vector space+</b>+<p>+ Definition. A set V of elements x ,y ,z ,...is called a vector+ (or linear) space over a complex field C if+<ul>+<li>+ (a) vector addition + is defined in V such that V is an+ abelian group under addition, with identity element 0+<pre>+ 1: <b>x</b> + <b>y</b> = <b>y</b> + <b>x</b>+ 2: <b>x</b> + (<b>y</b> + <b>z</b>) = (<b>x</b> + <b>y</b>) + <b>z</b>+ 3: <b>0</b> + <b>x</b> = <b>x</b> + <b>0</b>++</pre>+<p>+<li>+ (b) the set is close with respect to scalar multiplication+ and vector addition+<pre>+ 4: a (<b>x</b> + <b>y</b>) = a <b>x</b> + a <b>y</b>+ 5: (a + b) <b>x</b> = a <b>x</b> + b <b>x</b>+ 6: a (b <b>x</b>) = (a b) <b>x</b>+ 7: 1 <b>x</b> = <b>x</b>+ 8: 0 <b>x</b> = <b>0</b>+ where+ a, b, c are complex scalars+</pre>+</ul>+ Definition. The maximum number of linearly independent vectors+ in V or, what is the same thing, the minimum number of linearly+ independent vectors required to span V is the dimension r of+ vector space V.+<p>+ Definition. A set of r linearly independent vectors is called+ a basis of the space. Each vector of the space is then a unique+ linear combination of the vectors of this basis.+<p>+ Based on the above definitions we will define two vector+ spaces: ket space and its dual -- bra space, which, in addition+ to the above properties, will also support+ several common operations -- grouped below in the class+ DiracVector.+<pre>++> class DiracVector a where+> add :: a -> a -> a+> scale :: Scalar -> a -> a+> reduce :: a -> a+> basis :: a -> [a]+> components :: a -> [Scalar]+> compose :: [Scalar] -> [a] -> a+> dimension :: a -> Int+> norm :: a -> Double+> normalize :: a -> a++> dimension x = length (basis x)+>+> normalize x+> | normx == 0 = x+> | otherwise = compose cs (basis x)+> where+> cs = [a*v :+ b*v |a :+ b <- components x]+> v = 1 / normx+> normx = norm x++</pre>+<p>+<hr>+<p>+<b>+ 3. Ket vector space+</b>+<p>+ We submit that the following datatype and accompanying+ operations define a complex vector space, which we will call+ the ket vector space.+<pre>++> type Scalar = Complex Double++> data Ket a =+> KetZero -- zero ket vector+> | Ket a -- base ket vector+> | Scalar :|> Ket a -- scaling ket vectors+> | Ket a :+> Ket a -- spanning ket space++</pre>++ A tensor product of two ket spaces is also a ket space.+<pre>++> (*>) :: (Ord a, Ord b) => Ket a -> Ket b -> Ket (Tuple a b)+> Ket a *> Ket b = Ket (a :* b)+> x *> KetZero = KetZero+> KetZero *> y = KetZero+> x *> y = foldl1 (:+>) [((Bra a <> x) * (Bra b <> y)) :|> Ket (a :* b)+> | Ket a <- basis x, Ket b <- basis y]+++> (|>) :: Ord a => Scalar -> Ket a -> Ket a+> --+> -- Multiplication of ket by scalar+> --+> s |> (x :+> y) = (s |> x) +> (s |> y)+> s |> KetZero = KetZero+> 0 |> x = KetZero+> s |> (s2 :|> x) = (s * s2) |> x+> s |> x = s :|> x+++> (+>) :: Ord a => Ket a -> Ket a -> Ket a+> --+> -- Addition of two kets+> --+> x +> KetZero = x+> KetZero +> x = x+> x +> y = reduce (x :+> y)+++> instance (Eq a, Ord a) => Eq (Ket a) where+> --+> -- Two ket vectors are equal if they have identical+> -- components+> --+> x == y = and [c k x == c k y | k <- basis x]+> where+> c k x = (toBra k) <> x+++</pre>+ The data Ket is parametrized by type variable "a", which can be+ anything that can be compared for equality and ordered: integer,+ tuple, list of integers, etc. For example, the data+ constructor <code>Ket (3::Int)</code> creates a base vector <code>|3></code>,+ annotated by Int.+ Similarly, <code>Ket (2::Int,1::Int)</code>, creates a base vector+ <code>|(2,1)></code> annotated by a tuple of Ints. Those two+ vectors belong to two different bases.+<p>+ The eight examples below illustrate the eight defining equations+ of the vector space, given in section 1. All of them evaluate+ to True.+<pre>++ 1: Ket 2 +> Ket 3 == Ket 3 +> Ket 2+ 2: Ket 1 +> (Ket 2 +> Ket 3) == (Ket 1 +> Ket 2) +> Ket 3+ 3: Ket 1 +> KetZero == KetZero +> Ket 1+ 4: 5 |> (Ket 2 +> Ket 3) == 5 |> Ket 2 +> 5 |> Ket 3+ 5: (5 + 7) |> Ket 2 == 5 |> Ket 2 +> 7 |> Ket 2+ 6: 2 |> (4 |> Ket 2) == 8 |> Ket 2+ 7: 1 |> Ket 2 == Ket 2+ 8: 0 |> Ket 2 == KetZero+</pre>+ The ket expressions can be pretty printed, as shown below.+<pre>+ Ket 2 +> Ket 3 ==> 1.0 |2> + 1.0 |3>+ 5 |> (Ket 2 +> Ket 3) ==> 5.0 |2> + 5.0 |3>+ 2 |> (4 |> Ket 2) ==> 8.0 |2>+</pre>+ In order to support all those identities we also need several+ additional functions for reducing the vector to its canonical form,+ for composing the ket vector, and for extracting the ket+ basis and the ket components -- as shown below.+<pre>+++> reduceKet :: Ord a => Ket a -> Ket a+> reduceKet x+> --+> -- Reduce vector `x' to its canonical form+> --+> = compose cs ks+> where+> ks = basis x+> cs = [toBra k <> x | k <- ks]+++> ketBasis :: Ord a => Ket a -> [Ket a]+> --+> -- Sorted list of unique base vectors of the ket vector+> --+> ketBasis KetZero = []+> ketBasis (Ket k) = [Ket k]+> ketBasis (s :|> x) = [x]+> ketBasis (k1 :+> k2) = nub (ketBasis k1 ++ ketBasis k2)+++> toBra :: Ord a => Ket a -> Bra a+> --+> -- Convert from ket to bra vector+> --+> toBra (Ket k) = Bra k+> toBra (x :+> y) = toBra x :<+ toBra y+> toBra (p :|> x) = (conjugate p) :<| toBra x+++> instance Ord a => DiracVector (Ket a) where+> add = (+>)+> scale = (|>)+> reduce = reduceKet+> basis = ketBasis+> components x = [toBra e <> x | e <- basis x]+> compose xs ks = foldl1 (:+>) [fst z :|> snd z | z <- zip xs ks]+>+> norm KetZero = 0+> norm x = sqrt $ realPart (toBra x <> x)+++</pre>+ But those auxilliary functions refer to vectors from the+ conjugated space bra, which we shall now define below.+<p>+<hr>+<p>+<b>+ 4. Bra vector space+</b>+<p>+ Definition. Let V be the defining n-dimensional complex vector+ space. Associate with the defining n-dimensional complex vector+ space V a conjugate (or dual) n-dimensional vector space+ obtained by complex conjugation of elements x in V.+<p>+ We will call this space the bra space, and the corresponding vectors+ - the bra vectors. Further, we submit that the following datatype and the corresponding+ operations define bra space in Haskell.+<pre>++> data Bra a =+> BraZero -- zero bra vector+> | Bra a -- base bra vector+> | Scalar :<| Bra a -- scaling bra vectors+> | Bra a :<+ Bra a -- spanning bra space+++</pre>+ A tensor product of two bra spaces is also a bra space.+<pre>++> (<*) :: (Ord a, Ord b) => Bra a -> Bra b -> Bra (Tuple a b)+> Bra a <* Bra b = Bra (a :* b)+> x <* BraZero = BraZero+> BraZero <* y = BraZero+> x <* y = foldl1 (:<+) [((x <> Ket a) * (y <> Ket b)) :<| Bra (a :* b)+> | Bra a <- basis x, Bra b <- basis y]++> (<|) :: Ord a => Scalar -> Bra a -> Bra a+> s <| (x :<+ y) = (s <| x) <+ (s <| y)+> s <| BraZero = BraZero+> 0 <| x = BraZero+> s <| (s2 :<| x) = (s * s2) <| x+> s <| x = s :<| x+++> (<+) :: Ord a => Bra a -> Bra a -> Bra a+> --+> -- Sum of two bra vectors+> --+> x <+ BraZero = x+> BraZero <+ x = x+> x <+ y = reduce (x :<+ y)+++> instance (Eq a, Ord a) => Eq (Bra a) where+> --+> -- Two bra vectors are equal if they have+> -- identical components+> --+> --+> x == y = and [c b x == c b y | b <- basis x]+> where+> c b x = x <> toKet b++</pre>++ Similarly to what we have done for ket vectors, we also define several+ additional functions for reducing the bra vector to its canonical form,+ for composing the bra vector, and for extracting the bra+ basis and the bra components -- as shown below.+<pre>++> reduceBra :: Ord a => Bra a -> Bra a+> reduceBra x+> --+> -- Reduce bra vector `x' to its canonical form+> --+> = compose cs bs+> where+> bs = basis x+> cs = [x <> toKet b | b <- bs]+++> braBasis :: Ord a => Bra a -> [Bra a]+> --+> -- List of unique basis of the bra vector+> --+> braBasis BraZero = []+> braBasis (Bra b) = [Bra b]+> braBasis (s :<| x) = [x]+> braBasis (b1 :<+ b2) = nub (braBasis b1 ++ braBasis b2)+++> toKet :: Ord a => Bra a -> Ket a+> --+> -- Convert from bra to ket vector+> --+> toKet (Bra k) = Ket k+> toKet (x :<+ y) = toKet x :+> toKet y+> toKet (p :<| Bra k) = (conjugate p) :|> Ket k+++> instance Ord a => DiracVector (Bra a) where+> add = (<+)+> scale = (<|)+> reduce = reduceBra+> basis = braBasis+> components x = [x <> toKet e | e <- basis x]+> compose xs ks = foldl1 (:<+) [fst z :<| snd z | z <- zip xs ks]+>+> norm BraZero = 0+> norm x = sqrt $ realPart (x <> toKet x)+++</pre>+<p>+<hr>+<p>+<b>+ 5. Bra and Ket spaces as inner product spaces+</b>+<p>++ Definition. A complex vector space V is an inner product space+ if with every pair of elements x ,y from V there is associated+ a unique inner (or scalar) product < x | y > from C, such that+<pre>+ 9: < x | y > = < y | x ><sup>*</sup>+ 10: < a x | b y > = a<sup>*</sup> b < x | y >+ 11: < z | a x + b y > = a < z | x > + b < z, y >+ where+ a, b, c are the complex scalars+</pre>+ We submit that the dual ket and bra spaces are inner product+ spaces, providing that the inner product is defined by the operator+ <> given below:+<pre>++++> (<>) :: Ord a => Bra a -> Ket a -> Scalar+> --+> -- Inner product, or the "bra-ket" product+> --+> BraZero <> x = 0+> x <> KetZero = 0+> Bra i <> Ket j = d i j+> (p :<| x) <> (q :|> y) = p * q * (x <> y)+> (p :<| x) <> y = p * (x <> y)+> x <> (q :|> y) = q * (x <> y)+> x <> (y1 :+> y2) = (x <> y1) + (x <> y2)+> (x1 :<+ x2) <> y = (x1 <> y) + (x2 <> y)+++> d :: Eq a => a -> a -> Scalar+> d i j+> --+> -- Classical Kronecker's delta+> -- for instances of Eq class+> --+> | i == j = 1+> | otherwise = 0+>++</pre>+ The expressions below illustrate the definitions 9-11.+ They are all true.+<pre>+9: (toBra x <> y) == conjugate (toBra y <> x)+10: (toBra (a |> x) <> (b |> y)) == (conjugate a)*b*(toBra x <> y)+11: (toBra z <> (a |> x +> b |> y)) == a*(toBra z <> x) + b*(toBra z <> y)+ where+ x = (2 :+ 3) |> Ket 2+ y = ((1:+2) |> Ket 3) +> Ket 2+ z = Ket 2 +> Ket 3+ a = 2:+1+ b = 1+</pre>+<p>+<hr>+<p>+<b>+ 6. Linear operators+</b>+<p>++ Linear operators, or simply operators, are functions from vector+ in representation a <em>a</em> to vector in representation <em>b</em>++<pre>+ a :: Ket a -> Ket b+</pre>+ although quite often the operations are performed+ on the same representation. The linear operators A are defined by+<pre>+ A (c1 | x > + c2 | y > ) = c1 A | x > + c2 A | y >+</pre>++<p>+ We will describe variety of special types+ of operators, such as inverse, unitary, adjoint and hermitian.+ This is not an accident that the names of those operators+ resemble names from matrix calculus, since+ Dirac vectors and operators can be viewed as matrices.+<p>+ With the exception of variety of examples, no significant+ amount of Haskell code will be added here. This section+ is devoted mainly to documentation; we feel that it is important+ to provide clear definitions of the operators, as seen from+ the Haskell perspective. Being a strongly typed language,+ Haskell might not allow for certain relations often shown+ in traditional matrix calculus, such as+<pre>+ A = B+</pre>+ since the two operators might have in fact two distinct signatures.+ In matrix calculus one only compares tables of unnamed numbers,+ while in our Haskell formalism we compare typed+ entieties.+ For this reason, we will be threading quite+ slowly here, from one definition to another to assure that+ they are correct from the perspective of+ typing rules of Haskell.++<p>+<hr>+<p>+<b>+ 6.1. Operator notation+</b>+<p>+ The notation+<pre>+ | y > = A | x >+</pre>+ is pretty obvious: operator A acting on vector | x > produces+ vector | y >. It is not obvious though whether both vectors+ use the same representation. The Haskell version of the above+ clarifies this point, as in this example:+<pre>+ y = a >< x+ where+ a :: Ket Int -> Ket (Int, Int)+ a = ......+</pre>+ In this case it is seen the two vectors have distinct+ representations. The operator >< will be explained soon+ but for now treat is as an application of an operator+ to a vector, or some kind of a product of the two.+<p>+ The above can be also written as+<pre>+ | y > = | A x >+</pre>+ where the right hand side is just a defining label saying that the+ resulting vector has been produced by operator A acting on | x >.+<p>+ Linear operators can also act on the bra vectors+<pre>+ < y | = < x | A+ <---+</pre>+ providing that they have correct signatures. This postfix notation+ though is a bit awkward, and not supported by Haskell. To avoid+ confusion we will be using the following notation instead:+<pre>+ < y | = < A x |+</pre>+ which says that bra y is obtained from ket y,+ where | y > = | A x >, as before. In Haskell we will write+ it as+<pre>+ y = toBra $ a >< x++</pre>++<p>+<hr>+<p>+<b>+ 6.2. Renaming the representation+</b>+<p>+ One simple example of an operator is <em>label "new"</em>+ which renames a vector representation by adding extra label+ <em>"new"</em> in the basis vectors <em>Ket a</em>. Silly+ as it sounds, this and other similar re-labeling operations+ can be actually quite useful; for example,+ we might wish to distinguish between old and new bases, or+ just to satisfy the Haskell typechecker.+<pre>++ label :: (Ord a, Ord b) => b -> Ket a -> Ket (b, a)+ label i (Ket a) = Ket (i, a)+ label i x = (label i) >< x++</pre>+<p>+<hr>+<p>+<b>+ 6.3. Closure formula, or identity operator+</b>+<p>+ Although the general Dirac formalism often refers to+ abstract vectors | x >, our implementation must+ be more concrete than that -- we always represent the+ abstract vectors in some basis of our choice, as in:+<pre>+ | x > = c<sub>k</sub> | k > (sum over k)+</pre>+ To recover the component c<sub>k</sub> we form+ the inner product+<pre>+ c<sub>k</sub> = < k | x >+</pre>+ Putting it back to the previous equation:+<pre>+ | x > = < k | x > | k > (sum over k)+ = | k > < k | x >+ = Id | x >+ where+ Id = | k > < k | (sum over k)+</pre>+ we can see that the vector | x > has been abstracted away. The formula+ says that vector | x > can be decomposed in any basis+ by applying identity operator Id to it. This is also known+ as a closure formula. Well, Haskell has the "id" function too,+ and we could apply it to any ket, as in:+<pre>+ id (Ket 1 +> 10 |> Ket 2) ==> | 1 > + 10 | 2 >+</pre>+ but Haskell's "id" does not know anything about representations;+ it just gives us back the same vector | x > in our original+ representation.+<p>+ We need something more accurately depicting the closure+ formula | k > < k |, that would allow us to change+ the representation if we wanted to, or leave it alone+ otherwise. Here is the <em>closure</em> function and+ coresponding operator (><) that implement+ the closure formula for a given <em>operator</em>.+<pre>++> closure :: (DiracVector a, DiracVector b) => (a -> b) -> a -> b+> closure operator x =+> compose' (components x) (map operator (basis x))+> where+> compose' xs ks = foldl1 add (zipWith scale xs ks)++> operator >< x = closure operator x+++</pre>+<p>+<hr>+<p>+<b>+ 6.4. Changing the representation+</b>+<p>+ The silly <em>label</em> function found in the comment of the+ section 6.1 uses in fact the closure relation. But we could+ define is simpler than that:+<pre>++> label i (Ket x) = Ket (i, x)++</pre>+ and then apply a closure to a vector x, as in:+<pre>+ closure (label 0) (Ket 2 +> 7 |> Ket 3)+ ==> 1.0 |(0,2)> + 7.0 |(0,3)>+</pre>+ Somewhat more realistic example involves "rotation" of+ the old basis with simulaneous base renaming:+<pre>++> rot :: Ket Int -> Ket (Int, Int)+> rot (Ket 1) = normalize $ Ket (1,1) +> Ket (1,2)+> rot (Ket 2) = normalize $ Ket (1,1) +> (-1) |> Ket (1,2)+> rot (Ket _) = error "exceeded space dimension"++</pre>+ The example function+ <em>rot</em> assumes transformation from+ two-dimensional basis [| 1 >, | 2 >] to another+ two-dimensional basis [| (1,1) >, | (1,2) >] by+ expressing the old basis by the new one. Given this+ transformation we can apply the closure to any vector | x >+ represented in the old basis; as a result we will get+ the same vector | x > but represented in the new+ basis.+<pre>+ rot >< (Ket 1 +> 7 |> Ket 2) ==>+ 5.65685 |(1,1)> + -4.24264 |(1,2)>+</pre>++<p>+<hr>+<p>+<b>+ 6.5. Implementation of the operator equation A | x > = | y >+</b>+<p>+ The Haskell implementation of the closure formula is not just+ a useless simulation of the theoretical closure - it is one of the+ workhorses of the apparatus employed here.+<p>+ We will be using linear operators to evaluate equations+ like this:+<pre>+ | y > = A | x >+</pre>+ The resulting vector | y > can have either the same+ representation as | x > or different - depending on+ the nature of operator A. The most general type of+ A is+<pre>+ Ket a -> Ket b+</pre>+ but more often than not the basis will be the same as before.+ But how we define the operator A itself? The best way is+ to specify how it acts on the base vectors | k >. If we can chose+ as our basis the eigenvectors of A this would be even better,+ because the definition of A would be then extremely simple.+ After inserting the identity | k >< k | between the+ operator A and vector | x > in the above equation one gets+<pre>+ | y > = A | k > < k | x > (sum over k)+</pre>+ This will be implemented in Haskell as:+<pre>+ y = a >< x+</pre>+ The closure formula will take care of the rest and it will+ produce the result | y > . The examples previously given+ do just that. One caveat though: since operator A will+ only be defined for the basis, but not for other vectors,+ skipping the closure formula and coding directly+<pre>+ y = a' x+</pre>+ is not advisable.+ This will certainly fail for vectors other than basis unless+ one makes extra provisions for that. This is what we did+ in module Momenta, before we had the closure support ready.+ Using the closure is safe and this is the way to go!+++<p>+<hr>+<p>+<b>+ 6.6. Inverse operator+</b>+<p>+ An operator B = A<sup>-1</sup> that inverses the+ equation+<pre>+ | y > = A | x >+ y = a >< x -- where a :: Ket a -> Ket b+</pre>+ into+<pre>+ | x > = B | y >+ x = b >< y -- where b :: Ket b -> Ket a+</pre>+ is called the inverse operator.+<p>+ For example, the inverse operator to the operator <em>label i</em>+ is:+<pre>++> label' :: (Ord a, Ord b) => Ket (a, b) -> Ket b+> label' (Ket (i, x)) = Ket x++</pre>+ It is easy to check that applying the operator A and its inverse+ A<sup>-1</sup> in succession to any ket | x > one should+ obtain the same vector | x > again, as in:++<pre>+ A<sup>-1</sup> A | x > = | x >++ -- Haskell example+ label' >< (label 0 >< x) == x+ where+ x = Ket 1 +> 10 |> Ket 7+ ==> True+</pre>+ Once again, notice the omnipresent closure operator in Haskell+ implementation. Tempting as it might be to implement the+ above example as+<pre>+ -- Do not do it in Haskell!!!+ (label' . label 0) >< x == x+ where+ x = Ket 1 +> 10 |> Ket 7+ ==> True+</pre>+ this is not a recommended way. Although this example would work,+ but a similar example for <em>rotation</em> operations would+ fail in a spectacular way. The correct way is to insert the+ closure operator between two rotations:+<pre>+ rot' >< (rot >< x) == x+ where+ x = Ket 1 +> 10 |> Ket 2+ ==> True+</pre>+ where the inverse operator <em>rot'</em> is defined below:++<pre>++> rot' :: Ket (Int, Int) -> Ket (Int)+> rot' (Ket (1,1)) = normalize $ Ket 1 +> Ket 2+> rot' (Ket (1,2)) = normalize $ Ket 1 +> (-1) |> Ket 2+> rot' (Ket (_,_)) = error "exceeded space dimension"++</pre>+<p>+<hr>+<p>+<b>+ 6.7. Matrix representation of an operator+</b>+<p>+<p>+ The scalar products+<pre>+ < k | A l' > = < k | A | l' >+</pre>+ such that | k > and | l' > are the base vectors+ (in general belonging to two different bases), form a transformation+ matrix Akl'.+<p>+ In Haskell this matrix is formed as+<pre>+ k <> a >< l'+ where+ k = ... :: Bra b+ l' = ... :: Ket a+ a = ... :: Ket a -> Ket b+</pre>++<p>+<hr>+<p>+<b>+ 6.8. Adjoint operator+</b>+<font color="teal">+<p>+ Our definition of adjoint operator is different+ than that in theory of determinants. Many books, not necessarily+ quantum mechanical oriented, refer to the latter as <em>+ classical adjoint operator</em>.+</font>++<p>+ With every linear operator A we can associate an adjoint+ operator B = A<sup>+</sup>, also known as Hermitian conjugate+ operator, such that equality of the two scalar+ products+<pre>+ < A<sup>+</sup> u | x > = < u | A x >+</pre>+ holds for every vector | u > and | x >.+ In Haskell notation the above can be written as:+<pre>+ (toBra (b >< u) <> x) == toBra u <> a >< x+ where+ a = ... :: Ket a -> Ket b+ b = ... :: Ket b -> Ket a+ x = ... :: Ket a+ u = ... :: Ket b++</pre>+ For example, the operator <em>rot'</em> is adjoint+ to operator <em>rot</em>+<pre>+ (toBra (rot' >< u) <> x) == (toBra u <> rot >< x)+ where+ x = Ket 1 +> 10 |> Ket 2+ u = Ket (1,1) +> 4 |> Ket (1,2)+ ==> True++</pre>+ It can be shown that+<pre>+ (A<sup>+</sup>)<sup>+</sup> = A+</pre>+ Matrix A<sup>+</sup> is conjugate transposed to A, as+ proven below++<pre>+ = A<sup>+</sup>kl'+ = < k | A<sup>+</sup> | l' >+ = < k | A<sup>+</sup> l' >+ = < A<sup>+</sup> l' | k ><sup>*</sup>+ = < l' | A | k ><sup>*</sup>+ = A<sup>*</sup>l'k+</pre>+++<p>+<hr>+<p>+<b>+ 6.9. Unitary operator+</b>+<p>+ Unitary transformations preserve norms of vectors.+ We say, that the norm of a vector is invariant under unitary+ transformation.+ Operators describing such transformations are called+ unitary operators.+<pre>+ < A x | A x > = < x | x >++</pre>+ The example of this is rotation transformation, which indeed+ preserves the norm of any vector x, as shown in this Haskell+ example+<pre>+ (toBra u <> u) == (toBra x <> x)+ where+ u = rot >< x+ x = Ket 1 +> 10 |> Ket 2++ ==> True+</pre>+<p>+ Inverse and adjoint operators of unitary operators are equal+<pre>+ A<sup>-1</sup> = A<sup>+</sup>+</pre>+ which indeed is true for our example operator <em>rot</em>.+<p>+ Computation of the adjont operators A<sup>+</sup> from A+ is quite easy since the process is rather mechanical, as+ described in the previous section. On the other hand, finding+ inverse operators is not that easy, with the exception of+ some simple cases, such as our example 2D rotation.+ It is therefore important to know whether a given operator+ is unitary, as this would allow us to replace inverse+ operators by adjoint operators.+++<p>+<hr>+<p>+<b>+ 6.10. Hermitian operator+</b>+<p>+ A Hermitian operator is a self adjoint operator; that is+<pre>+ < A u | x > = < u | A x >+</pre>+ Another words: A<sup>+</sup> = A.+<p>+ Notice however, that this relation holds only for the+ vectors in the same representation, since in general+ the operators+ A and A<sup>+</sup> have distinct signatures, unless+ types a, b are the same:+<pre>+ a :: Ket a -> Ket b -- operator A+ a' :: Ket b -> Ket a -- operator A<sup>+</sup>+</pre>+ Elements of hermitian matrices must therefore satisfy:+<pre>+ Aij = (Aji)<sup>*</sup>+</pre>+ In particular, their diagonal elements must be real.+<p>+ Our example operator <em>rot</em> is not hermitian,+ since it describes transformation from one basis+ to another.+ But here is a simple example of a hermitian operator, which+ multiplies any ket by scalar 4. It satisfies our definition:+<pre>+ (toBra (a >< u) <> x) == (toBra u <> a >< x)+ where+ a v = 4 |> v++ x = Ket 1 +> Ket 2+ u = Ket 2++ ==> True+</pre>+ Here is a short quote from [3].+<blockquote>+ Why do we care whether an operator is Hermitian?+ It's because of a few theorems:++<ol>+<li>+ The eigenvalues of Hermitian operators are always real.+<li>+ The expectation values of Hermitian operators are always real.+<li>+ The eigenvectors of Hermitian operators span the Hilbert space.+<li>+ The eigenvectors of Hermitian operators belonging to distinct eigenvalues are orthogonal.+</ol>+ In quantum mechanics, these characteristics are essential if you+ want to represent measurements with operators. Operators must be+ Hermitian so that observables are real. And, you must be able to+ expand in the eigenfunctions - the expansion coefficients+ give you probabilities!+</blockquote>+<p>+<hr>+<p>+<b>+ 7. Showing kets and bras+</b>+<p>+ Lastly, here are show functions for pretty printing of Dirac+ vectors.+<pre>++> instance (Show a, Eq a, Ord a) => Show (Ket a) where+> showsPrec n KetZero = showString "| Zero >"+> showsPrec n (Ket j) = showString "|" . showsPrec n j . showString ">"+> showsPrec n (x :|> k) = showsScalar n x . showsPrec n k+> showsPrec n (j :+> k) = showsPrec n j . showString " + " . showsPrec n k++> instance (Show a, Eq a, Ord a) => Show (Bra a) where+> showsPrec n BraZero = showString "< Zero |"+> showsPrec n (Bra j) = showString "<" . showsPrec n j . showString "|"+> showsPrec n (x :<| k) = showsScalar n x . showsPrec n k+> showsPrec n (j :<+ k) = showsPrec n j . showString " + " . showsPrec n k+++> showsScalar n x@(a :+ b)+> | b == 0 = showsPrec n a . showString " "+> | otherwise = showString "(" .showsPrec n x . showString ") "++</pre>+<p>+<hr>+<p>+<b>+ 8. Data Tuple for tensor products+</b>+<p>+ A state vector of several subsystems is modelled as a ket parametrized+ by a type variable Tuple, which is similar to ordinary () but is+ shown differently. Tensor product of several simple states leads+ to deeply entangled structure, with many parenthesis obstructing+ readability. What we really want is a simple notation for easy+ visualization of products of several states, as in:+<pre>+ Ket 1 *> Ket (2, 1) * Ket '+' ==> | 1; (2,1); '+' >+</pre>+ See module Momenta for practical example of tensor products+ of vector spaces.+<pre>++> data Tuple a b = a :* b+> deriving (Eq, Ord)++> instance (Show a, Show b) => Show (Tuple a b) where+> showsPrec n (a :* b) = showsPrec n a . showString "; " . showsPrec n b++</pre>+<p>+<hr>+<p>+<b>+ 9. References+</b>+<p>+<ul>+<p>+<li>++ [1] Jerzy Karczmarczuk, Scientific computation and functional+ programming, Dept. of Computer Science, University of Caen, France,+ Jan 20, 1999, <a href="http://www.info.unicaen.fr/~karczma/">+ http://www.info.unicaen.fr/~karczma/</a>+<p>+<li>+ [2] Jan Skibinski, Collection of Haskell modules,+ Numeric Quest Inc., <a href="http://www.numeric-quest.com/haskell/">+ http://www.numeric-quest.com/haskell/"</a>+<p>+<li>+ [3] Steven Pollock, University of Colorado,+ <a href="http://www.colorado.edu/physics/phys3220/3220_fa97/notes/notes_table.html">+ Quantum Mechanics, Physics 3220 Fall 97, lecture notes</a>++</ul>+<p>+<hr>+<p>+<b>+ 10. Copyright and license+</b>++<pre>+--+-- Copyright:+--+-- (C) 2000 Numeric Quest, All rights reserved+--+-- Email: jans@numeric-quest.com+--+-- http://www.numeric-quest.com+--+-- License:+--+-- GNU General Public License, GPL+--++</pre>+</blockquote>+</body>++<SCRIPT language="Javascript">+<!--++// FILE ARCHIVED ON 20010421035521 AND RETRIEVED FROM THE+// INTERNET ARCHIVE ON 20030715011358.+// JAVASCRIPT APPENDED BY WAYBACK MACHINE, COPYRIGHT INTERNET ARCHIVE.+// ALL OTHER CONTENT MAY ALSO BE PROTECTED BY COPYRIGHT (17 U.S.C.+// SECTION 108(a)(3)).++ var sWayBackCGI = "http://web.archive.org/web/20010421035521/";++ function xLateUrl(aCollection, sProp) {+ var i = 0;+ for(i = 0; i < aCollection.length; i++)+ if (aCollection[i][sProp].indexOf("mailto:") == -1 &&+ aCollection[i][sProp].indexOf("javascript:") == -1)+ aCollection[i][sProp] = sWayBackCGI + aCollection[i][sProp];+ }++ if (document.links) xLateUrl(document.links, "href");+ if (document.images) xLateUrl(document.images, "src");+ if (document.embeds) xLateUrl(document.embeds, "src");++ if (document.body && document.body.background)+ document.body.background = sWayBackCGI + document.body.background;++//-->++</SCRIPT>+</html>
+ README view
@@ -0,0 +1,5 @@+http://web.archive.org/web/20010520121707/www.numeric-quest.com/haskell/++The Literate Haskell files are actually HTML files.+To make your browser happy, you can start 'make html'+in order to make links *.html links to *.lhs files.
+ Roots.hs view
@@ -0,0 +1,110 @@+module Roots where +import Data.Complex+import Data.List(genericLength)++roots :: RealFloat a => a -> Int -> [Complex a] -> [Complex a]+roots eps count as =+ --+ -- List of complex roots of a polynomial+ -- a0 + a1*x + a2*x^2...+ -- represented by the list as=[a0,a1,a2...]+ -- where+ -- eps is a desired accuracy+ -- count is a maximum count of iterations allowed+ -- Require: list 'as' must have at least two elements+ -- and the last element must not be zero + roots' eps count as []+ where+ roots' eps count as xs + | length as <= 2 = x:xs+ | otherwise = + roots' eps count (deflate x bs [last as]) (x:xs)+ where+ x = laguerre eps count as 0+ bs = drop 1 (reverse (drop 1 as))+ deflate z bs cs+ | bs == [] = cs+ | otherwise = + deflate z (tail bs) (((head bs)+z*(head cs)):cs)+++laguerre :: RealFloat a => a -> Int -> [Complex a] -> Complex a -> Complex a +laguerre eps count as x+ --+ -- One of the roots of the polynomial 'as',+ -- where+ -- eps is a desired accuracy+ -- count is a maximum count of iterations allowed+ -- x is initial guess of the root+ -- This method is due to Laguerre.+ --+ | count <= 0 = x+ | magnitude (x - x') < eps = x'+ | otherwise = laguerre eps (count - 1) as x'+ where+ x' = laguerre2 eps as as' as'' x+ as' = polynomial_derivative as+ as'' = polynomial_derivative as' + laguerre2 eps as as' as'' x+ -- One iteration step+ | magnitude b < eps = x+ | magnitude gp < magnitude gm = + if gm == 0 then x - 1 else x - n/gm+ | otherwise = + if gp == 0 then x - 1 else x - n/gp+ where+ gp = g + delta+ gm = g - delta+ g = d/b+ delta = sqrt ((n-1)*(n*h - g2))+ h = g2 - f/b+ b = polynomial_value as x+ d = polynomial_value as' x+ f = polynomial_value as'' x+ g2 = g^2+ n = genericLength as++polynomial_value :: Num a => [a] -> a -> a+polynomial_value as x =+ --+ -- Value of polynomial a0 + a1 x + a2 x^2 ...+ -- evaluated for 'x',+ -- where 'as' is a list [a0,a1,a2...]+ --+ foldr (u x) 0 as+ where+ u x a b = a + b*x++polynomial_derivative :: Num a => [a] -> [a]+polynomial_derivative as+ --+ -- List of coefficients for derivative of polynomial+ -- a0 + a1 x + a2 x^2 ...+ --+ | as == [] = []+ | otherwise = deriv 1 (drop 1 as) []+ where+ deriv n bs cs+ | bs == [] = reverse2 cs+ | otherwise = deriv (n+1) (tail bs) ((n*(head bs)):cs) + reverse2 cs+ | cs == [] = []+ | otherwise = reverse cs++++-----------------------------------------------------------------------------+--+-- Copyright:+--+-- (C) 1998 Numeric Quest Inc., All rights reserved+--+-- Email:+--+-- jans@numeric-quest.com+--+-- License:+--+-- GNU General Public License, GPL+-- +-----------------------------------------------------------------------------
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ Tensor.lhs view
@@ -0,0 +1,974 @@+<html>+<head>+<BASE HREF="http://www.numeric-quest.com/haskell/Tensor.html">++<title>+ N-dimensional tensors+</title>+</head>+<body>+<ul>+<center>+<h1>+ ***+</h1>+<h1>+ N-dimensional tensors+</h1>+<b>+<br>+ Literate Haskell module <i>Tensor.lhs</i>+</b>+<p>+ Jan Skibinski, <a href="http://www.numeric-quest.com/news/">+ Numeric Quest Inc.</a>, Huntsville, Ontario, Canada+<p>+ 1999.10.08, last modified 1999.10.16++</center>+<p>+<hr>+<p>+<i>+ This is a quick sketch of what might be a basis of a real+ Tensor module. This module has quite a few limitations (listed below).+ I'd like to get some feedback on what should be a better+ way to design it properly. Nevertheless, this module works+ and is able to tackle complex and mundane manipulations+ in the very straightforward way.+<p>+ There are few arbitrary decisions we have taken. For example,+ we consider a scalar to be a tensor of rank 0. This forces us to+ do conversions between true scalars and such tensors, but it also+ saves us a lot of headache related to typing restrictions. This+ is a typical price paid for (too much?) generalization.+<p>+ To get rid of those awful sums appearing in multiplications+ of tensors we do introduce Einstein's summation convention by the way of+ text examples -- followed by the equivalent Haskell examples.+ Hopefully it is clear and be well appreciated for its economy+ of notation, which is standard in the tensor calculus.+<p>+ Datatype <code>Tensor</code> defined here is an instance+ of class <code>Eq</code>, <code>Show</code> and <code>Num</code>.+ That means that one can compare tensors for equality and perform+ basic numerical calculations, such as addition, negation,+ subtraction, multiplication, etc. -- using standard notation+ <code>(==), (/=), (+), (-), (*)</code>. In addition, several+ customized operations, such as <code> (<*>)</code>+ and <code>(<<*>>)</code> are defined for+ variety of inner products.++<p>+ Limitations of this module:+<ul>++<p>+<li>+ Tensor components are Doubles. Why not Fraction, Complex, etc?+ For a moment we will leave this question aside, and+ return to it some time later. But we consider it+ the important question -- which is evident from the attempts of+ such generalization in some of our other modules:+ <a href="http://www.numeric-quest.com/haskell/Orthogonals.html">+ Orthogonals</a> and+ <a href="http://www.numeric-quest.com/haskell/fractions.html">+ Fraction</a>.+<p>+<li>+ We are well aware that the decision to represent tensors+ as nested objects will have significant impact on access+ (and update -- if supported) of such data structure. Linear+ arrays seem to be better suited for such tasks, where all+ indices must be explicitely computed first, but the access+ time is linear. In contrary, the hierarchical data structure+ defined here require very little effort in index computing+ but the access time depends on the depth of the data tree.+<p>+ But speed has not been tested yet, so we really do not know+ how inefficient this module is and all of the above is+ just a pure speculation. Certain operations of this module+ seem to be quite well matched with this tree-like data structure,+ and because of it this design decision might be not so bad+ after all.++<p>+<li>+ The shape of tensors defined here involves two parameters:+ dimension and rank. Rank is associated with the+ depth of the tensor tree and corresponds to a total number+ of indices by which you can access the individual components.+ No limits are imposed on ranks and there are binary operations+ which involve tensors of different ranks.+ Dimension is associated with the breadth of the tree and+ correspond to a number of values each index can take.+ Dimension is fixed via constant <code>dims</code>. At first it might+ seem as a severe limitation, but in fact one should never+ mix tensors with different dimensions. One usually works+ either with three-dimensional tensors (classical mechanics,+ electrodynamics, elasticity, etc.) or the four-dimentional+ tensors (relativity theory).+</ul>+<p>+</i>+<p>+<hr>+<p>+<b>+ Tensor datatype+</b>+<p>+<pre>++> module Tensor where+> import Data.Array(inRange)+> infixl 9 # -- used for tensor indexing+> infixl 9 ## -- used for indices expressed as lists+> infixl 7 <*> -- inner product with one bound+> infixl 7 <<*>> -- inner product with two bounds++</pre>+ Indices will assume values from range (1,dims) (defined below).+<p>++ Tensor can contain a scalar value or a list of tensors.+ This recursively defines tensor of any rank in n-D space.++<pre>++> data Tensor = S Double+> | T [Tensor]+++</pre>+ There is no way we could specify the length of the list+ <code>[Tensor]</code> in the data declaration. Typing is not+ concerned with shapes.+ We could of course use more specific representation of+ this data structure, such as:+<pre>+ data Tensor = S Double | T Tensor Tensor Tensor+</pre>+ but then we would severily limit ourselves to three-dimensional+ tensors.+<p>++ Rank is either 0 (scalars), 1 (vectors), or higher: 2, 3, 4 ...+<pre>++> rank :: Tensor -> Int+> rank t = rank' 0 t where+> rank' n (S a) = n+> rank' n (T xs) = rank' (n+1) (head xs)++</pre>+ Here we define our tensor dimension as constant for this+ module. All binary operations on tensors require the+ same dimensions, so it makes sense to treat dimensions+ as constants. But ranks can be different.+<pre>++> dims :: Int+> dims = 3++</pre>++<p>+<hr>+<p>+<b>+ Showing+</b>+<p>+ Tensors are printed as recursive lists with a word "Tensor"+ prepended++<pre>++> instance Show Tensor where+> showsPrec 0 (S a) = showString "Tensor " . showsPrec 0 a+> showsPrec n (S a) = showsPrec n a++> showsPrec 0 (T xs) = showString "Tensor " . showList' 0 xs+> showsPrec n (T xs) = showList' n xs++> showList' n [] = showString "[]"+> showList' n (x:xs) = showChar '[' . showsPrec (n+1) x . showRem (n+1) xs+> where+> showRem n [] = showChar ']'+> showRem n (x:xs) = showChar ',' . showsPrec n x . showRem n xs++</pre>+<p>+<hr>+<p>+<b>+ Input+</b>+<p>++ Although tensors are printed as structured list+ it is easier to input data via flat lists.+ But make sure that the length of the list is one+ of: dims^0, dims^1, dims^2, dims^3, dims^4, etc.+<p>+ This function is quite inefficient for ranks higher than 4.+ Compare, for example, timings of:+<pre>+ tensor [1..3^6]+ tensor [1..3^3] * tensor [1..3^3]+</pre>+ Although both expressions create tensors of the same rank 6,+ but the execution of the latter is much faster. This is+ because the function <code>tensor</code> spends much+ of its effort on recursively restructuring the flat lists+ into the lists-of-lists-of-lists...+<pre>+++> tensor :: [Double] -> Tensor+> tensor xs+> | size == 1 = S (head xs)+> | q /= 0 = error "Length is not a power of dims"+> | otherwise = T (tlist p xs)+> where+> (p,q) = rnk 1 (quotRem size dims)+> rnk m (1, v) = (m, v)+> rnk m (u, 0) = rnk (m+1) (quotRem u dims)+> rnk m (u, v) = (m, v)+> size = length xs+> group n xs = group' n xs [] where+> group' n xs as+> | length xs == 0 = reverse as+> | length xs < n = reverse (xs:as)+> | otherwise = group' n (drop n xs) ((take n xs):as)+>+> tlist 1 xs = map S xs+> tlist rnk xs = tlist' (rnk-1) (map S xs)+> where+> tlist' 0 zs = zs+> tlist' n zs = tlist' (n-1) (map T (group dims zs))+>++</pre>+<p>+<hr>+<p>+<b>+ Extraction and conversion+</b>+<p>++ Tensor components are also tensors and can be extracted+ via (#) operator++<pre>++> ( # ) :: Tensor -> Int -> Tensor+> a@(S a1) # 1 = S a1+> a@(S a1) # i = error "out of range"+> a@(T xs) # i = xs!!(i-1)++> ( ## ) :: Tensor -> [Int] -> Tensor+> a ## [] = a+> a ## (x:xs) = (a#x) ## xs++</pre>++ Tensors of rank 0 can be converted to scalars; i.e.,+ simple numbers of type Double.+<pre>++> scalar :: Tensor -> Double+> scalar (S a) = a+> scalar (T xs) = error "rank not 0"++</pre>++ Tensors of rank 1 can be converted to vectors; i.e.,+ lists with "dims" components of type Double+<pre>++> vector :: Tensor -> [Double]+> vector (S a) = error "rank not 1"+> vector a@(T xs)+> | rank a /= 1 = error "rank not 1"+> | otherwise = map scalar xs++</pre>+<p>+<hr>+<p>+<b>+ Useful tensors: epsilon and delta+</b>+<p>+ Function "epsilon' i j k" emulates values of the pseudo-tensor Eijk.+ It is valid only for three-dimensional tensors.+ It takes three indices i,j,k from the range (1,3)+ and returns one of the three values:+ 0.0, 1.0, -1.0 -- depending on the rules specified below:+<pre>++> epsilon' :: Int -> Int -> Int -> Double+> epsilon' i j k+> | dims /= 3 = error "not 3-dims"+> | outside (1,3) i j k = error "Not in range"+> | (i == j) || (i == k) || (j == k) = 0+> | otherwise = epsilon1 i j k+> where+> epsilon1 i j k+> | (i == 1) && (j == 2) && (k == 3) = 1+> | (i == 3) && (j == 2) && (k == 1) = -1+> | otherwise = epsilon1 j k i+> outside (p,q) a b c =+> (not $ inRange (p,q) a) ||+> (not $ inRange (p,q) b) ||+> (not $ inRange (p,q) c)++</pre>+ Function "delta' i j" emulates Kronecker's delta:+<pre>++> delta' :: Int -> Int -> Double+> delta' i j+> | i == j = 1+> | otherwise = 0++</pre>++ Delta' and epsilon' can be converted to tensors++<pre>++> delta = tensor [delta' i j | i <- [1..dims], j <- [1..dims]]+> epsilon = tensor [epsilon' i j k | i <- [1..3], j <- [1..3], k <- [1..3]]++</pre>+ The components delta[ij] and epsilon[i,j,k] can be extracted+ and converted to numbers. For example:+<pre>+ scalar (epsilon#1#2#3) = 1+ scalar (epsilon#1#1#3) = 0,+ scalar (epsilon#3#2#1) = -1+</pre>+<p>+<hr>+<p>+<b>+ Dot product+</b>+<p>+ Dot product of two tensors of rank 1 could be defined as+ tensor of rank 0. This is not the most efficient implementation+ but we still want the dot product to be recognised as+ tensor, so we loose on speed here:+<pre>++> dot :: Tensor -> Tensor -> Tensor+> dot a b = S (sum [scalar (a#i) * scalar (b#i) | i <- [1..dims]])++</pre>++<p>+<hr>+<p>+<b>+ Cross product - valid for 3D space only+</b>+<p>+ The cross product of two vectors is another vector:+ C = A x B. The pseudotensor Eijk is used to compute+ such cross product.+<p>+ First, here are numerical components of C, C[i]:+<pre>++> cross' :: Tensor -> Tensor -> Int -> Double+> cross' a b i = sum [(epsilon' i j k)* scalar (a#j) * scalar (b#k)|+> j<-[1..3],k<-[1..3], j/=k]++</pre>+ And here is the full vector C (as tensor of rank 1):++<pre>++> cross :: Tensor -> Tensor -> Tensor+> cross a b = tensor (map (cross' a b) [1..3])++</pre>++ Example:+<pre>+ cross (tensor [1..3]) (tensor [1,8,1]) ==> Tensor [-22.0, 2.0, 6.0]+</pre>++<p>+<hr>+<p>+<b>+ Equality of tensors+</b>+<p>+ Tensor can be admitted to class <code>Eq</code>. We only need to+ define either equality or nonequality operation. We've chosen+ to define the former: two tensors are equal if they have the same+ rank and equal components:+<pre>++> instance Eq Tensor where+> (==) a b+> | ranka /= rank b = False+> | ranka == 0 = scalar a == scalar b+> | otherwise = and [(a#i) == (b#i) | i <- [1..dims]]+> where+> ranka = rank a+>++</pre>+++<p>+<hr>+<p>+<b>+ Tensor as instance of class Num+</b>+<p>+ To admit tensors to class <code>Num</code> we have to+ support all the operations from that class. Here is+ the class Num declaration taken from the Prelude:+<pre>+class (Eq a, Show a) => Num a where+ (+), (-), (*) :: a -> a -> a+ negate :: a -> a+ abs, signum :: a -> a+ fromInteger :: Integer -> a++ -- Minimal complete definition: All, except negate or (-)+ x - y = x + negate y+ negate x = 0 - x+</pre>+ All operations but <code>(*)</code> are straightforward,+ meaningful and easy to implement. The semantics of multiplication+ <code>(*)</code> is, however, not so obvious and it is up to us+ how to define it: as an inner product or as an outer+ product. We have chosen the latter, which means that the+ operation <code>c = a * b</code> produces a new tensor <code>c</code>+ whose rank is a sum of the ranks of tensors being+ multiplied:+<pre>+ rank c = rank a + rank b+</pre>++ Suffice to add that tensor products are generally not+ commutative; that is:+<pre>+ a * b /= b * a++</pre>+ That said, here is the instantiation of <code>Num</code>+ for datatype Tensor:+<pre>++> instance Num Tensor where+> (+) a b+> | ranka /= rank b = error "different ranks"+> | ranka == 0 = S (scalar a + scalar b)+> | otherwise = T [a#i + b#i | i <- [1..dims]]+> where+> ranka = rank a++> negate a@(S a1) = S (negate a1)+> negate a@(T xs) = T (map negate xs)++> abs a@(S a1) = S (abs a1)+> abs a@(T xs) = T (map abs xs)++> signum a@(S a1) = S (signum a1)+> signum a@(T xs) = T (map signum xs)++> fromInteger n = S (fromInteger n)++> (*) a@(S a1) b@(S b1) = S (a1*b1)+> (*) a@(S a1) b@(T xs) = T (map (a*) (take dims xs))+> (*) a@(T xs) b = T (map (*b) (take dims xs))++</pre>+ Having defined the operation <code>(*)</code> as an outer product+ such operation will generally increase the rank of the outcome.+ For example, if <code>a</code> is a tensor of rank 2 (matrix) and+ <code>b</code> is a tensor of rank 1 (vector) then the result is+ a tensor of rank 3:+<pre>+ c = a * b, that is+ c[ijk] = a[ij] b[k]+</pre>+ But this is not what is typically considered a multiplication+ of tensors; we are more often than not interested in the inner+ products, informally described below.++<p>+<hr>+<p>+<b>+ Contraction+</b>+<p>+<p>+ Eistein's indexing convention of tensors is based on+ the distinction between free indices and bound indices.+ Free indices appear in the tensorial expressions, such+ as <code>A[ijkl]</code>, once only and they indicate+ a freedom for substitution of any specific index+ from the range of valid indices. This range is (1,3)+ for 3D tensors. The expression <code>A[ijkl]</code>+ represents in fact one of 3^4 possible components+ of the tensor <code>A</code>.+<p>+ Bound indices, on the other hand, appear in pairs+ (and only in pairs) and they indicate the summation of+ tensor expression over the valid range. For example,+<pre>+ A[kkj] = A[11j] + A[22j] + A[33j]+</pre>+ Note that the index "j" is still free, and that means+ that the above represents three equations for j = 1,2,3.+<p>+ A process of converting of a pair of free indices+ to a pair of bound indices is called contraction. As+ a result a rank of a tensor (or expression involving+ several tensors) is being reduced+ by two.+<p>+ The function <code>contract</code> below accepts a tensor of a+ rank bigger or equal 2 and two integers m,n from the range (1,rank a)+ which indicate positions of the two indices to be used for+ contraction. The result is a tensor with its rank reduced+ by two.++<pre>+++> contract :: Int -> Int -> Tensor -> Tensor+> contract m n a+> | m >= n = error "wrong ordering"+> | outside m n = error "not in range"+> | ranka < 2 = error "cannot contract"+> | ranka == 2 = S (sum [scalar (a#i#i) | i <- [1..dims]])+> | ranka > 2 = tensor [summa m n us a | us <- freeIndices (ranka-2)]+> where+> ranka = rank a+>+> outside p q = (not $ inRange (1,ranka) p)+> ||(not $ inRange (1,ranka) q)+> summa p q xs a = sum [scalar (a##(insert p q xs r)) |+> r <- [1..dims]]++> -- Insert element r at positions m n to the list+> -- of indices xs+> insert m n xs r = us++[r]++ws++[r]++zs+> where+> (us,vs) = splitAt (m-1) xs+> (ws,zs) = splitAt (n - m - 1) vs+>+> freeIndices 1 = [[x] | x <- [1..dims]]+> freeIndices n = [x:y | x <- [1..dims], y <- freeIndices (n-1)]+++</pre>++ Let's take for example tensor <code>delta</code> and contract+ it in its two indices:+<pre>+ delta [kk] = delta[1,1] + delta[2,2] + delta[3,3] = 1 + 1 + 1 = 3+</pre>+ The same can be done in Haskell:+<pre>+ contract 1 2 delta ==> Tensor 3.0+ rank (contract 1 2 delta) ==> 0+</pre>+++<p>+<hr>+<p>+<b>+ Inner product+</b>+<p>+ The inner product of two tensors can be considered+ as two-phase process: first the outer product is+ formed and then a contraction is applied to a selected+ pair of indices. There are countless possibilities+ of defining such inner products, since we can choose+ any pair, or even more than one pair, of indices+ to become bound.+<p>+ How do we usually multiply tensors? Here is one example,+ which is equivalent to matrix-vector multiplication:+<pre>+ C[i] = A[ij] B[j]+</pre>+ Notice two types of indices: index "i" is free since+ it appears only once on both sides of the equation. It means+ that you can freely substitute 1,2 or 3 for "i". So in fact+ we have here three equations:+<pre>+ C[1] = A[1j] B[j]+ C[2] = A[2j] B[j]+ C[3] = A[3j] B[j]+</pre>+ Index "j" is bound - it appears two times on the right hand+ side, but not on the left side. Bound indices signify summation+ from 1 to 3. So the above in fact means:+<pre>+ C[1] = A[11] B[1] + A[12] B[2] + A[13] B[3]+ C[2] = A[21] B[1] + A[22] B[2] + A[23] B[3]+ C[3] = A[31] B[1] + A[32] B[2] + A[33] B[3]+</pre>+ The economy of notation is evident in our first form above.+ How will we do it in Haskell?+<p>+ To obtain the above result we will first form the outer product+ of matrix A and vector B, obtain a tensor of rank 3,+ and then contract it in indices 2 and 3 to obtain a+ the final expected result (inner product):+<pre>+ c = contract 2 3 (a * b)+</pre>+ This approach is quite inefficient storage-wise and+ speed-wise and a direct customized encoding which avoids creating+ outer products is recommended instead.+<p>+ The system of equations+<pre>+ C[i] = A[ij] B[j]+</pre>+ could obviously be represented explicite as:+<pre>+ c i = sum [scalar(a#i#j) * scalar(b#j) | j <- [1..dims]]+ -- valid for i = 1..dims+</pre>+ But when efficiency is not a premium we could still+ take advantage of function <code>contract</code>+ to write clear code that avoids the explicit sums. The+ operator <code> <*></code>, introduced below, allows+ us to write the same function as:+<pre>+ c = a <*> b -- the output is a tensor of rank 1+ c' i = (a <*> b)#i -- the output is a tensor of rank 0+ c'' i = scalar ((a <*> b)#i) -- the output is a number+</pre>++<p>+<hr>+<p>+<b>+ Convenience operators for inner products+</b>+<p>+ Variety of specialized functions for inner products+ could be defined. We will show few examples here+ and introduce specialized convenience operators+ for most common types of inner products. Please+ note that the proposed operators are not standard+ in any way, and we are not trying to suggest that+ they are important. Just treat them as examples.+<p>+ The semantics of operator <code> <*> </code> has+ been chosen to support matrix-vector or vector-matrix+ multiplications. But this operator is more general+ than that, because it also handles products with scalars+ (tensors of rank 0), and generally any products+ of any two tensors with bounds imposed on one pair+ of indices: last index of the first tensor and first+ index of the second tensor.++<pre>++> a <*> b+> | (ranka == 0) || (rankb == 0) = a * b+> | otherwise = contract ranka (ranka + 1) (a * b)+> where+> ranka = rank a+> rankb = rank b++</pre>++ Take for example a classical identity:+<pre>+ A[i] = delta[ij] B[j], where delta is a Kronecker's delta+</pre>+ Here is an example of how we can use it in Haskell:+<pre>+ delta <*> tensor [4,5,6]) ==> Tensor [4.0, 5.0, 6.0]+ (delta <*> tensor [4,5,6])#1 ==> Tensor 4.0+</pre>++ Let's try something more complex, for example a constitutive equation+ relating the stress tensor S[ij] with the deformation tensor G[kl].+ The tensor C[ijkl] is an anisotropic tensor of material constants:+ 81 altogether. In fact, due to all sorts of symmetries this number+ could be reduced to twenty-something for the most complex crystals,+ and to two independent components for the isotropic materials.+ Anyway, the relation is linear and can be written as follows:+<pre>+ S[ij] = C[ijkl] G[kl]+</pre>+ This represents 9 equations (i,j->1,2,3) and expands heavily+ to sums over k and l on the right-hand side.+ We need to impose two bounds in two pairs of indices to+ support above example. Here is another specialized operator+ for inner product with two specificly selected bounds.+<pre>++> a <<*>> b+> | (ranka < 2) || (rankb < 2) = error "rank too small"+> | otherwise = contract (ranka-1) ranka+> (contract ranka (ranka+2) (a * b))+> where+> ranka = rank a+> rankb = rank b++</pre>+ Here is a dummy, but easy to generate example of the above:++<pre>+ tensor [1..81] <<*>> tensor [1..9]++ ==> s = Tensor [[ 285.0, 690.0, 1095.0],+ [1500.0, 1905.0, 2310.0],+ [2715.0, 3120.0, 3525.0]]++ (tensor [1..81] <<*>> tensor [1..9])#1#1 = Tensor 285.0+</pre>+<p>+<hr>+<p>+<b>+ Double cross products+</b>+<p>+ Here is another useful example of tensor multiplication.+ Say you want to compute a cross product of three vectors:+<pre>+ D = C X (A x B )+</pre>+ In index notation this could be expressed as:+<pre>+ D[i] = E[ijk] C[j] E[kpq] A[p] B[q]+</pre>+ This represents three equations for i=1,2,3. All other indices+ j,k,p,q are bound; that is, they appear in pairs on the right+ hand side, indicating four sums. Although you can calculate+ it directly, and this Haskell module can do it easily, we can+ simplify this equation by organizing it differently and+ using this identity:+<pre>+ E[ijk] = E[kij]+</pre>+ (Even permutation of indices does not change a sign of pseudo-tensor+ E.)+<pre>+ D[i] = E[kij] E[kpq] C[j] A[p] B[q]+</pre>+ Now here is another useful identity, which gets rid of the+ bound index "k" (sitting in the first position above):+<pre>+ E[kij] E[kpq] = delta[ip] delta[jq] - delta[iq] delta[jp]+</pre>+ After substitution and using identity <code>delta[ij] G[j] = G[i]</code>+ the <code>C x (A x B)</code> transforms to:+<pre>+ D[i] = C[j] B[j] A[i] - C[j] A[j] B[i]+</pre>+ We still have three scalar equations, but they are less complex:+ there is only one summation (over the "j") on the right hand side.+<p>+ You should easily recognize that <code>C[j] B[j]</code>+ represents the scalar product. Therefore our double cross product+ can be represented as a difference of two vectors:+<pre>+ D = C x (A x B) = (C o B) A - (C o A) B+</pre>++ Now, let us see how this module handles this. Let's take an+ example of three randomly chosen vectors A, B, C. The direct+ method is straightforward, although it involves quite a lot+ of multiplications and summations (which would not be so+ evident if we have not done all those preliminary examinations+ above).+<pre>++> d_standard = cross c (cross a b) where+> a = tensor [1,2,3]+> b = tensor [3,1,8]+> c = tensor [5,2,4]++</pre>+ On the other hand we could encode the equivalent equation:+<pre>+ D = (C o B) A - (C o A) B+</pre>+ as:+<pre>++> d_simpler =+> tensor [n1 * scalar (a#i) - n2 * scalar (b#i) | i <- [1..dims]] where+>+> a = tensor [1,2,3]+> b = tensor [3,1,8]+> c = tensor [5,2,4]+> n1 = scalar (c `dot` b)+> n2 = scalar (c `dot` a)++</pre>++ Both <code>d_standard</code> and <code>d_simpler</code>+ lead to the same result:+<pre>+ ==> Tensor [-14.0, 77.0, -21.0]+</pre>+<p>+<hr>+<p>+<b>+ Vector transformation+</b>+<p>+ A vector can be decomposed in any system of reference. The best+ choice is any orthogonal system of reference, where all base+ unit vectors are mutually perpendicular (orthogonal), since this+ simplifies the computations. The base vectors <code>e[1], e[2], e[3]</code>+ are usually chosen as vectors of length one (we say that they are+ normalized to one), and hence they are called "orthonormal".+ They obey the orthonormality relations for their scalar products:+<pre>+ e[i] o e[j] = delta[ij]+</pre>+ where the Kronecker's "delta" has been defined before.+<p>+ Here is an example of the vector decomposition:+<pre>+ A = A[i] e[i] (summation over "i"!)+</pre>+ The components A[i] of the vector A obviously depend on the choice+ of the base system. The same vector A will have different+ components in two different systems of references:+<pre>+ A'[i] e'[i] = A[i] e[i]+</pre>+ where primes refer to the new system. Now, if we multiply both+ sides of the above equation by a base vector <code>e'[k]</code>,+ using the scalar (dot) product definition, we will get:+<pre>+ A'[i] e'[k] o e'[i] = A[i] e'[k] o e[i]+</pre>+ The new base vectors are mutually orthonormal, so+<pre>+ e'[k] o e'[i] = delta[ki]+</pre>+ and the left hand side will be transformed to:+<pre>+ A'[i] delta[ki] = A'[k]+</pre>+ But the base vectors on the right hand side are taken from+ two different systems, and therefore they are not mutually+ orthonormal. All such nine scalar products form the components of the+ transormation tensor, R:+<pre>+ R[ki] = e'[k] o e[i]+</pre>+ As a result, our original equation can be expressed as+ a new equation defining transformation of the vector A:+<pre>+ A'[k] = R[ki] A[i]+</pre>+ This gives us a rule how to compute new components A'[k] of vector+ A from its old components and transformation tensor R[ki].+<p>+ You might want to run some exercise choosing the old+ system with the base vectors:+<pre>+ e#1=tensor [1,0,0]+ e#2=tensor [0,1,0]+ e#3=tensor [0,0,1],+</pre>+ where "e" can be considered a tensor of rank 2:+<pre>+ e = tensor [1,0,0,+ 0,1,0,+ 0,0,1]+</pre>+ and the new system obtained from the old one by rotation+ around the axis 3 (x3, or z) by an angle "alpha". Some+ trigonometry will be involved to compute the new base+ vectors, e'[i]. The next step is to compute tensor R[ki]+<pre>++ r = tensor [scalar (e'#k `dot` e#i)|k<-[1..dims], i<-[1..dims]]++</pre>+ and finally use operator <code> <*></code> to compute new components+ of vector A:+<pre>+ a' = r <*> a+</pre>+<p>+<hr>+<p>+ Related page on this site:+ <a href="http://www.numeric-quest.com/haskell/index.html">+ Collection of Haskell modules</a>++<pre>+-----------------------------------------------------------------------------+--+-- Copyright:+--+-- (C) 1999 Numeric Quest Inc., All rights reserved+--+-- Email:+--+-- jans@numeric-quest.com+--+-- License:+--+-- GNU General Public License, GPL+--+-----------------------------------------------------------------------------+</pre>+</ul>+</body>++<SCRIPT language="Javascript">+<!--++// FILE ARCHIVED ON 20010630021753 AND RETRIEVED FROM THE+// INTERNET ARCHIVE ON 20030626102305.+// JAVASCRIPT APPENDED BY WAYBACK MACHINE, COPYRIGHT INTERNET ARCHIVE.+// ALL OTHER CONTENT MAY ALSO BE PROTECTED BY COPYRIGHT (17 U.S.C.+// SECTION 108(a)(3)).++ var sWayBackCGI = "http://web.archive.org/web/20010630021753/";++ function xLateUrl(aCollection, sProp) {+ var i = 0;+ for(i = 0; i < aCollection.length; i++)+ if (aCollection[i][sProp].indexOf("mailto:") == -1 &&+ aCollection[i][sProp].indexOf("javascript:") == -1)+ aCollection[i][sProp] = sWayBackCGI + aCollection[i][sProp];+ }++ if (document.links) xLateUrl(document.links, "href");+ if (document.images) xLateUrl(document.images, "src");+ if (document.embeds) xLateUrl(document.embeds, "src");++ if (document.body && document.body.background)+ document.body.background = sWayBackCGI + document.body.background;++//-->++</SCRIPT>+</html>
+ numeric-quest.cabal view
@@ -0,0 +1,25 @@+Name: numeric-quest+Version: 0.1+License: GPL+Author: Jan Skibinski+Maintainer: Henning Thielemann <numeric-quest@henning-thielemann.de>+Homepage: http://www.haskell.org/haskellwiki/Numeric_Quest+Package-URL: http://darcs.haskell.org/numeric-quest/+Category: Math+Build-Depends: base, haskell98+Synopsis: Math and quantum mechanics+Description: List based linear algebra, similtaneous linear equations, eigenvalues and eigenvectors, roots of polynomials, transcendent functions with arbitrary precision implemented by continued fractions, quantum operations, tensors+GHC-Options: -Wall+Hs-source-dirs: .+Exposed-modules:+ Eigensystem+ EigensystemNum+ Fraction+ LinearAlgorithms+ Orthogonals+ QuantumVector+ Roots+ Tensor+Data-Files:+ Makefile+ README