diff --git a/Eigensystem.hs b/Eigensystem.hs new file mode 100644 --- /dev/null +++ b/Eigensystem.hs @@ -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 +-- +--------------------------------------------------------------------------- + + diff --git a/EigensystemNum.hs b/EigensystemNum.hs new file mode 100644 --- /dev/null +++ b/EigensystemNum.hs @@ -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" diff --git a/Fraction.hs b/Fraction.hs new file mode 100644 --- /dev/null +++ b/Fraction.hs @@ -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 +-- +-- +-- +-- http://www.numeric-quest.com +-- +-- License: +-- +-- GNU General Public License, GPL +-- +----------------------------------------------------------------------------- diff --git a/LinearAlgorithms.hs b/LinearAlgorithms.hs new file mode 100644 --- /dev/null +++ b/LinearAlgorithms.hs @@ -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 +-- +--------------------------------------------------------------------------- + diff --git a/Makefile b/Makefile new file mode 100644 --- /dev/null +++ b/Makefile @@ -0,0 +1,5 @@ + +html: Orthogonals.html QuantumVector.html Tensor.html + +%.html: %.lhs + ln -s $< $@ diff --git a/Orthogonals.lhs b/Orthogonals.lhs new file mode 100644 --- /dev/null +++ b/Orthogonals.lhs @@ -0,0 +1,1869 @@ + + + + + + Indexless linear algebra algorithms + + + + + + + + diff --git a/QuantumVector.lhs b/QuantumVector.lhs new file mode 100644 --- /dev/null +++ b/QuantumVector.lhs @@ -0,0 +1,1245 @@ + + + + + + Quantum vector + + + +
+

+ *** +

+

+ Quantum vector +

+

+ + + Jan Skibinski, + Numeric Quest Inc., Huntsville, Ontario, Canada +
+ Literate Haskell module QuantumVector.lhs +

+ Initialized: 2000-05-31, last modified: 2000-06-10 + +

+ +
+ +

+

+

+ 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. +

+ 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]. +

+ 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. + +

+ 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 + Momenta -- a non-trivial example designed to compute Clebsch-Gordan coefficients + of a transformation from one basis of angular momenta to another. +

+ 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, x <> y, + and a closure operation, a >< x, which can be considered + an application of a quantum operator a to a vector x. + At this stage our formalism applies only to discrete cases, but + we hope to generalize it on true Hilbert space as well. + +

+

+

+ + Contents + +

+ +

+

+

+ + 1. Infix operators + +

+ 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. + + +

+
+> 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
+
+
+

+

+

+ + 2. Vector space + +

+ Definition. A set V of elements x ,y ,z ,...is called a vector + (or linear) space over a complex field C if +

+ 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. +

+ 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. +

+ 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. +

+
+> 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
+
+
+

+

+

+ + 3. Ket vector space + +

+ We submit that the following datatype and accompanying + operations define a complex vector space, which we will call + the ket vector space. +

+
+> 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
+
+
+ + A tensor product of two ket spaces is also a ket space. +
+
+> (*>) :: (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
+
+
+
+ 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 Ket (3::Int) creates a base vector |3>, + annotated by Int. + Similarly, Ket (2::Int,1::Int), creates a base vector + |(2,1)> annotated by a tuple of Ints. Those two + vectors belong to two different bases. +

+ The eight examples below illustrate the eight defining equations + of the vector space, given in section 1. All of them evaluate + to True. +

+
+	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
+
+ The ket expressions can be pretty printed, as shown below. +
+	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>
+
+ 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. +
+
+
+> 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)
+
+
+
+ But those auxilliary functions refer to vectors from the + conjugated space bra, which we shall now define below. +

+

+

+ + 4. Bra vector space + +

+ 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. +

+ 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. +

+
+> 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
+
+
+
+ A tensor product of two bra spaces is also a bra space. +
+
+> (<*) :: (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
+
+
+ + 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. +
+
+> 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)
+
+
+
+

+

+

+ + 5. Bra and Ket spaces as inner product spaces + +

+ + 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 +

+	9:  < x | y >          = < y | x >*
+	10: < a x | b y >      = a* b < x | y >
+	11: < z | a x + b y >  = a < z | x > + b < z, y >
+	    where
+		a, b, c are the complex scalars
+
+ 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: +
+
+
+
+> (<>) :: 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
+>
+
+
+ The expressions below illustrate the definitions 9-11. + They are all true. +
+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
+
+

+

+

+ + 6. Linear operators + +

+ + Linear operators, or simply operators, are functions from vector + in representation a a to vector in representation b + +

+	a :: Ket a -> Ket b
+
+ although quite often the operations are performed + on the same representation. The linear operators A are defined by +
+	A (c1 | x > + c2 | y > ) = c1 A | x > + c2 A | y >
+
+ +

+ 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. +

+ 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 +

+	A = B
+
+ 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. + +

+

+

+ + 6.1. Operator notation + +

+ The notation +

+	| y > = A | x >
+
+ 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: +
+	y = a >< x
+	   where
+		a :: Ket Int -> Ket (Int, Int)
+		a = ......
+
+ 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. +

+ The above can be also written as +

+	| y > = | A x >
+
+ where the right hand side is just a defining label saying that the + resulting vector has been produced by operator A acting on | x >. +

+ Linear operators can also act on the bra vectors +

+	< y | = < x | A
+		<---
+
+ 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: +
+	< y | = < A x |
+
+ which says that bra y is obtained from ket y, + where | y > = | A x >, as before. In Haskell we will write + it as +
+	y = toBra $ a >< x
+
+
+ +

+

+

+ + 6.2. Renaming the representation + +

+ One simple example of an operator is label "new" + which renames a vector representation by adding extra label + "new" in the basis vectors Ket a. 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. +

+
+	label :: (Ord a, Ord b) => b -> Ket a -> Ket (b, a)
+	label i (Ket a) = Ket (i, a)
+	label i x       = (label i) >< x
+
+
+

+

+

+ + 6.3. Closure formula, or identity operator + +

+ 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: +

+	| x > = ck | k >   (sum over k)
+
+ To recover the component ck we form + the inner product +
+	    ck = < k | x >
+
+ Putting it back to the previous equation: +
+	| x > = < k | x > | k >      (sum over k)
+	      = | k > < k | x >
+	      = Id | x >
+	where
+	    Id = | k > < k |        (sum over k)
+
+ 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: +
+	id (Ket 1 +> 10 |> Ket 2) ==> | 1 > + 10 | 2 >
+
+ but Haskell's "id" does not know anything about representations; + it just gives us back the same vector | x > in our original + representation. +

+ 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 closure function and + coresponding operator (><) that implement + the closure formula for a given operator. +

+
+> 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
+
+
+
+

+

+

+ + 6.4. Changing the representation + +

+ The silly label function found in the comment of the + section 6.1 uses in fact the closure relation. But we could + define is simpler than that: +

+
+> label i (Ket x) = Ket (i, x)
+
+
+ and then apply a closure to a vector x, as in: +
+	closure (label 0) (Ket 2 +> 7 |> Ket 3)
+		==> 1.0 |(0,2)> + 7.0 |(0,3)>
+
+ Somewhat more realistic example involves "rotation" of + the old basis with simulaneous base renaming: +
+
+> 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"
+
+
+ The example function + rot 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. +
+	rot >< (Ket 1 +> 7 |> Ket 2) ==>
+		5.65685 |(1,1)> + -4.24264 |(1,2)>
+
+ +

+

+

+ + 6.5. Implementation of the operator equation A | x > = | y > + +

+ 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. +

+ We will be using linear operators to evaluate equations + like this: +

+	| y > = A | x >
+
+ 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 +
+	Ket a -> Ket b
+
+ 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 +
+	| y > = A | k > < k | x >            (sum over k)
+
+ This will be implemented in Haskell as: +
+	y = a >< x
+
+ 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 +
+	y = a' x
+
+ 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! + + +

+

+

+ + 6.6. Inverse operator + +

+ An operator B = A-1 that inverses the + equation +

+	| y > = A | x >
+	  y   = a >< x -- where a :: Ket a -> Ket b
+
+ into +
+	| x > = B | y >
+	  x   = b >< y -- where b :: Ket b -> Ket a
+
+ is called the inverse operator. +

+ For example, the inverse operator to the operator label i + is: +

+
+> label' :: (Ord a, Ord b) => Ket (a, b) -> Ket b
+> label' (Ket (i, x)) = Ket x
+
+
+ It is easy to check that applying the operator A and its inverse + A-1 in succession to any ket | x > one should + obtain the same vector | x > again, as in: + +
+	A-1 A | x > = | x >
+
+	-- Haskell example
+	label' >< (label 0 >< x) == x
+	   where
+		x = Ket 1 +> 10 |> Ket 7
+	==> True
+
+ Once again, notice the omnipresent closure operator in Haskell + implementation. Tempting as it might be to implement the + above example as +
+	-- Do not do it in Haskell!!!
+	(label' . label 0) >< x == x
+	    where
+	       x = Ket 1 +> 10 |> Ket 7
+	==> True
+
+ this is not a recommended way. Although this example would work, + but a similar example for rotation operations would + fail in a spectacular way. The correct way is to insert the + closure operator between two rotations: +
+	rot' >< (rot >< x) == x
+	    where
+		x = Ket 1 +> 10 |> Ket 2
+	==> True
+
+ where the inverse operator rot' is defined below: + +
+
+> 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"
+
+
+

+

+

+ + 6.7. Matrix representation of an operator + +

+

+ The scalar products +

+	< k | A l' > = < k | A | l' >
+
+ such that | k > and | l' > are the base vectors + (in general belonging to two different bases), form a transformation + matrix Akl'. +

+ In Haskell this matrix is formed as +

+	k <> a >< l'
+	    where
+	       k  = ... :: Bra b
+	       l' = ... :: Ket a
+	       a  = ... :: Ket a -> Ket b
+
+ +

+

+

+ + 6.8. Adjoint operator + + +

+ 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 + classical adjoint operator. + + +

+ With every linear operator A we can associate an adjoint + operator B = A+, also known as Hermitian conjugate + operator, such that equality of the two scalar + products +

+	< A+ u | x > = < u | A x >
+
+ holds for every vector | u > and | x >. + In Haskell notation the above can be written as: +
+	(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
+
+
+ For example, the operator rot' is adjoint + to operator rot +
+	(toBra (rot' >< u) <> x) == (toBra u <> rot >< x)
+	    where
+		x = Ket 1 +> 10 |> Ket 2
+		u = Ket (1,1) +> 4 |> Ket (1,2)
+	==> True
+
+
+ It can be shown that +
+	(A+)+ = A
+
+ Matrix A+ is conjugate transposed to A, as + proven below + +
+	= A+kl'
+	= < k | A+ | l' >
+	= < k | A+ l' >
+	= < A+ l' | k >*
+	= < l' | A | k >*
+	= A*l'k
+
+ + +

+

+

+ + 6.9. Unitary operator + +

+ 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. +

+	< A x | A x > = < x | x >
+
+
+ The example of this is rotation transformation, which indeed + preserves the norm of any vector x, as shown in this Haskell + example +
+	(toBra u <> u) == (toBra x <> x)
+	    where
+		u = rot >< x
+		x = Ket 1 +> 10 |> Ket 2
+
+	==> True
+
+

+ Inverse and adjoint operators of unitary operators are equal +

+	A-1 = A+
+
+ which indeed is true for our example operator rot. +

+ Computation of the adjont operators A+ 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. + + +

+

+

+ + 6.10. Hermitian operator + +

+ A Hermitian operator is a self adjoint operator; that is +

+	< A u | x > = < u | A x >
+
+ Another words: A+ = A. +

+ Notice however, that this relation holds only for the + vectors in the same representation, since in general + the operators + A and A+ have distinct signatures, unless + types a, b are the same: +

+	a  :: Ket a -> Ket b -- operator A
+	a' :: Ket b -> Ket a -- operator A+
+
+ Elements of hermitian matrices must therefore satisfy: +
+	 Aij = (Aji)*
+
+ In particular, their diagonal elements must be real. +

+ Our example operator rot 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: +

+	(toBra (a >< u) <> x) == (toBra u <> a >< x)
+	where
+	    a v = 4 |> v
+
+	    x = Ket 1 +> Ket 2
+	    u = Ket 2
+
+	==> True
+
+ Here is a short quote from [3]. +
+ Why do we care whether an operator is Hermitian? + It's because of a few theorems: + +
    +
  1. + The eigenvalues of Hermitian operators are always real. +
  2. + The expectation values of Hermitian operators are always real. +
  3. + The eigenvectors of Hermitian operators span the Hilbert space. +
  4. + The eigenvectors of Hermitian operators belonging to distinct eigenvalues are orthogonal. +
+ 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! +
+

+

+

+ + 7. Showing kets and bras + +

+ Lastly, here are show functions for pretty printing of Dirac + vectors. +

+
+> 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 ") "
+
+
+

+

+

+ + 8. Data Tuple for tensor products + +

+ 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: +

+	Ket 1 *> Ket (2, 1) * Ket '+' ==> | 1; (2,1); '+' >
+
+ See module Momenta for practical example of tensor products + of vector spaces. +
+
+> 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
+
+
+

+

+

+ + 9. References + +

+

+

+

+

+ + 10. Copyright and license + + +

+--
+-- Copyright:
+--
+--      (C) 2000 Numeric Quest, All rights reserved
+--
+--      Email: jans@numeric-quest.com
+--
+--      http://www.numeric-quest.com
+--
+-- License:
+--
+--      GNU General Public License, GPL
+--
+
+
+
+ + + + diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -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. diff --git a/Roots.hs b/Roots.hs new file mode 100644 --- /dev/null +++ b/Roots.hs @@ -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 +-- +----------------------------------------------------------------------------- diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#! /usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/Tensor.lhs b/Tensor.lhs new file mode 100644 --- /dev/null +++ b/Tensor.lhs @@ -0,0 +1,974 @@ + + + + + + N-dimensional tensors + + + + + + + + diff --git a/numeric-quest.cabal b/numeric-quest.cabal new file mode 100644 --- /dev/null +++ b/numeric-quest.cabal @@ -0,0 +1,25 @@ +Name: numeric-quest +Version: 0.1 +License: GPL +Author: Jan Skibinski +Maintainer: Henning Thielemann +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