(* ::Package:: *)

(* :Title:  SF1626 *)
(* :Context: SF1626 *)
(* :Name: SF1626 *)
(* :Author: Goeran Andersson *)
(* :Email: goeran@kth.se *)  
(* :Summary: 
	This package provides some support methods for multivariable calculus. 
	The package is developed for educational purposes and is free ware. 
*)
(* :Package Version: 2011-03-18 *)
(* :Mathematica Version: 8.0 *)


BeginPackage["SF1626`"];


(* prepare for reload *)
Unprotect[
	Taylor,QuadraticFormQ,SymmetricMatrix,CompleteSquare,
	NegativeDefiniteMatrixQ,IndefiniteMatrixQ,SemidefiniteMatrixQ,PositiveSemidefiniteMatrixQ,NegativeSemidefiniteMatrixQ
];
ClearAll[
	Taylor,QuadraticFormQ,SymmetricMatrix,CompleteSquare,
	NegativeDefiniteMatrixQ,IndefiniteMatrixQ,SemidefiniteMatrixQ,PositiveSemidefiniteMatrixQ,NegativeSemidefiniteMatrixQ
];


(* Help info *)
CompleteSquare::usage="CompleteSquare[p, x] completes the quadratic polynomial and CompleteSquare[p, {x1, x2, ...}] completes a multivariable quadratic polynomial.";
SymmetricMatrix::usage="SymmetricMatrix[q, {x1, x2, ...}] returns the symmetric coefficient matrix A so that x.A.x = q.";
Taylor::usage="Taylor[f[x, y, ...], {x, y, ...}, {x0, y0, ...}, n] returns a truncated Taylor polynomial of degree n.";
QuadraticFormQ::usage="QuadraticFormQ[q, {x1, x2, ...}] returns True if q[x1, x2, ...] is a quadratic form, and False otherwise.";
NegativeDefiniteMatrixQ::usage="NegativeDefiniteMatrixQ[m] tests whether the Hermitian matrix m is a negative definite matrix.";
PositiveSemidefiniteMatrixQ::usage="PositiveSemidefiniteMatrixQ[m] tests whether Hermitian matrix m is a positive semidefinite matrix.";
NegativeSemidefiniteMatrixQ::usage="NegativeSemidefiniteMatrixQ[m] tests whether Hermitian matrix m is a negative semidefinite matrix.";
SemidefiniteMatrixQ::usage="SemidefiniteMatrixQ[m] tests whether Hermitian matrix m is a semidefinite matrix.";
IndefiniteMatrixQ::usage="IndefiniteMatrixQ[m] tests whether Hermitian matrix m is a indefinite matrix.";


(* package code *)
Begin["`Private`"];

(* CompleteSquare *)
CompleteSquare[q_List,x:Except[_List]]:=
Module[{a,b,c,\[FormalC]},
	If[Exponent[Last[q],x]==2,
		{c,b,a}=CoefficientList[Last[q],x];
		Join[Most[q],Simplify[{a (x+b/(2 a))^2,\[FormalC]-b^2/(4 a)}]]/.\[FormalC]->c,
		q
	]
];
CompleteSquare[q:Except[_List],x:Except[_List]]:=If[Exponent[q,x]!=2,q, Plus@@CompleteSquare[{q},x]]/;PolynomialQ[q,x];
CompleteSquare[q:Except[_List],{x:Except[_List]}]:=CompleteSquare[q,x];
CompleteSquare[q:Except[_List],x:{Except[_List], Except[_List]..}]:=Plus@@Fold[CompleteSquare,{q},Most[x]]/;QuadraticFormQ[q,x];
CompleteSquare[q:Except[_List],x:{Except[_List], Except[_List]..}]:=Plus@@Fold[CompleteSquare,{q},x]/;PolynomialQ[q,x];
CompleteSquare[q_,___]:=q;

(* SymmetricMatrix *)
SymmetricMatrix[q_,x:{_Symbol..} | {_?(Head[Head[#1]] === Symbol & )..}]:=
Module[{A,a,s,z},
	A=Array[a,{Length[x],Length[x]}];
	z=Thread[Flatten[A]->0];
	s=SolveAlways[x.A.x==q,x];
	A=A/. s[[1]]/. z;
	1/2 (A+Transpose[A])
]/;QuadraticFormQ[q,x];

(* Taylor *)
Taylor[f_,x:{_Symbol..} | {_?(Head[Head[#1]] === Symbol & )..},x0:{__},n_Integer?NonNegative]:=
Module[{d,h,t,i},
	d[i_]:=\!\(
\*SubscriptBox[\(\[PartialD]\), \({x, i}\)]f\)/. Thread[x->x0];
	h[0]={};h[i_]:=h[i]=Join[h[i-1],{x-x0}];
	t[i_]:=Expand[Dot@@Join[{d[i]},h[i]]];
	\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 0\), \(n\)]
\*FractionBox[\(t[i]\), \(i!\)]\)
]/;Length[x]==Length[x0];

(* QuadraticFormQ *)
QuadraticFormQ[q_,x:{_Symbol..} | {_?(Head[Head[#1]] === Symbol & )..}]:=PolynomialQ[q,x]&&And@@Thread[Apply[Plus,First/@CoefficientRules[q,x],{1}]==2];
QuadraticFormQ[_,_]:=False;

(* ...DefiniteMatrixQ *)
(* PositiveDefiniteMatrixQ is included in Mathematica *)
NegativeDefiniteMatrixQ[m_?HermitianMatrixQ]:=PositiveDefiniteMatrixQ[-m];
PositiveSemidefiniteMatrixQ[m_?HermitianMatrixQ]:=Length[NullSpace[m]]>0&&And@@NonNegative[Chop[Eigenvalues[m]]];
NegativeSemidefiniteMatrixQ[m_?HermitianMatrixQ]:=Length[NullSpace[m]]>0&&And@@NonPositive[Chop[Eigenvalues[m]]];
SemidefiniteMatrixQ[m_?HermitianMatrixQ]:=PositiveSemidefiniteMatrixQ[m]||NegativeSemidefiniteMatrixQ[m]
IndefiniteMatrixQ[m_?HermitianMatrixQ]:=!Or@@PositiveDefiniteMatrixQ/@{m,-m}&&!SemidefiniteMatrixQ[m];

End[];


(* protection *)
Protect[
	CompleteSquare,SymmetricMatrix,Taylor,QuadraticFormQ,
	NegativeDefiniteMatrixQ,PositiveSemidefiniteMatrixQ,NegativeSemidefiniteMatrixQ,SemidefiniteMatrixQ,IndefiniteMatrixQ
];

EndPackage[];
