(*^ DELETE EVERYTHING ABOVE THIS LINE.
    DELETE EVERYTHING BELOW THE COMMENT LINE AT THE END OF THIS FILE.
    THEN SAVE THE FILE IN TEXT ONLY (ASCII) FORM.

Instructions for using this package, called HFT, are contained in the
LaTeX document "Mathematica and Harmonic Function Theory", which can
be obtained by e-mail from axler@math.msu.edu.

Load this file into Mathematica using the << command.

HFT

Initialization

Print

*)
If[ $VersionNumber < 2, 
Print["This package will not work correctly with versions"];
Print["of Mathematica earlier than Mathematica 2.0."] ]

Print[" "]

Print["HFT; Version 4.03, 1 October 1996"]

Print[" "]

Print["This Mathematica package is designed to"]

Print["accompany the book Harmonic Function Theory,"]

Print["by Sheldon Axler, Paul Bourdon, and Wade Ramey,"]

Print["published by Springer-Verlag."]

Print[" "]

Print["This package is copyrighted by the authors but"]

Print["is distributed free of charge by electronic mail."]

Print[" "]

Print["Comments, suggestions, bug reports, requests for"]

Print["the most recent version, and requests for the"]

Print["documentation should be sent by electronic mail to:"]

Print["axler@math.msu.edu"]

Print[" "]

Print["***********************************"]

Print["*  The computer is unpacking the  *"]

Print["*  HFT package.  This may take a  *"]

Print["*  minute or so.                  *"]

BeginPackage["HFT`"]

Unprotect[ Gradient, Jacobian ]

Remove[ Gradient, Jacobian ]
(*

Usage

*)
AnnularDirichlet::usage = 
	"AnnularDirichlet[ f, g, r, R, x ] gives the function of x " <>
	"harmonic on an annular region that equals f on the inner " <>
	"sphere of radius r and equals g on the outer sphere of radius R."

AntiLaplacian::usage = "AntiLaplacian[f, x] gives a polynomial " <>
	"whose Laplacian equals f.  " <>
	"Here f must be a polynomial function of x."

Ball::usage = 
	"Ball is one of the values that may be assigned to the option " <>
	"Orthonormal."

BasisH::usage = 
	"BasisH[m, x] gives a basis for the space of harmonic " <>
	"polynomials homogeneous of degree m in the variable x."

BergmanKernel::usage = "BergmanKernel[x, y] gives the Bergman " <>
	"reproducing kernel for the unit ball."

BergmanKernelH::usage = "BergmanKernel[z, w] gives the Bergman " <>
	"reproducing kernel for the upper half-space."

BergmanProjection::usage = "BergmanProjection[ f, x] gives the " <>
	"orthogonal projection of a polynomial f, as a function of x, " <>
	"onto the Bergman space of harmonic functions on the ball."

BiDirichlet::usage = "BiDirichlet[f , x] solves the " <>
	"BiDirichlet problem with boundary data f, "<>
	"as a function of x."

Delta::usage = "Delta[j] is the vector that equals 1 in the j-th " <>
	"coordinate and 0 in the other coordinates."

Dimension::usage = "Dimension[x] is the dimension of the vector x."

DimensionH::usage = 
	"DimensionH[m, n] gives the vector space dimension of the " <>
	"space of spherical harmonics of degree m in n variables."

Dirichlet::usage = "Dirichlet[f, x] gives the harmonic function " <>
	"that equals f on the unit sphere.  Here f must be a " <>
	"polynomial function of x."

Divergence::usage = "Divergence[f, x] gives the divergence of f " <>
	"with respect to x."

ExpandNorm::usage = 
	"ExpandNorm[f] gives f with all terms of the form |x + y| " <>
	"replaced by Sqrt[ |x|^2 + |y|^2 + 2 x.y ]."

ExteriorDirichlet::usage =
	"ExteriorDirichlet[f , x] solves the " <>
	"exterior Dirichlet problem with boundary data f, " <>
	"as a function of x."

ExteriorNeumann::usage =
	"ExteriorNeumann[f , x] solves the " <>
	"exterior Neumann problem with boundary data f, " <>
	"as a function of x."

Gradient::usage = "Gradient[f, x] gives the gradient of f with " <>
	"respect to x."

HarmonicConjugate::usage = "HarmonicConjugate[u, {x,y}] gives " <>
	"the harmonic conjugate of u with respect to the variables x, y."

HarmonicDecomposition::usage = "HarmonicDecomposition[u, x] " <>
	"gives the harmonic decomposition of u with respect to " <>
	"the variable x."

HilbertSchmidt::usage = "HilbertSchmidt[A] gives the " <>
	"HilbertSchmidt norm of a matrix A."

Homogeneous::usage = "Homogeneous[f, d, x] gives the homogenous " <>
	"part of f of degree d, with respect to x."

IntegrateBall::usage = "IntegrateBall[f, x] gives the integral " <>
	"of f, as a function of x, over the unit ball."

IntegrateSphere::usage = "IntegrateSphere[f, x] gives the " <>
	"integral of f, as a function of x, over the unit sphere " <>
	"with respect to normalized surface area measure."

Inversion::usage = "Inversion[x] gives the inversion of x."

Jacobian::usage = "Jacobian[f, x] gives the Jacobian of f with " <>
	"respect to x."

Kelvin::usage = "Kelvin[u, x] gives the Kelvin transform of u, " <>
	"thought of as a function of x."

KelvinM::usage = "KelvinM is the modified Kelvin transform, " <>
	"as defined in Chapter 7 of Harmonic Function Theory."

Laplacian::usage = "Laplacian[f, x] gives the Laplacian of f " <>
	"with respect to x."

Multiple::usage = 
	"Multiple is an option for AntiLaplacian.  " <>
	"The default value is None.  The value Norm^2 produces the " <>
	"unique antiLaplacian that is a polynomial multiple of " <>
	"Norm[x]^2, where x is the variable."

Neumann::usage = "Neumann[f, g, x] solves the " <>
	"Neumann problem of finding a function of x whose " <>
	"outward normal derivative (on the unit sphere) is f " <>
	"and whose Laplacian is g."

Norm::usage = "Norm[x] gives the Euclidean norm of x."

NormalD::usage = 
	"NormalD[f, z] gives the outward normal derivative of f, " <>
 	"as a function of z, with respect to the unit ball."

Orthonormal::usage = 
	"Orthonormal is an option for BasisH.  " <>
	"The default value is None.  The value Ball produces a " <>
	"basis that is orthonormal with respect to volume measure " <>
	"on the ball.  The value Sphere produces a " <>
	"basis that is orthonormal with respect to surface area " <>
	"measure on the sphere."

Partial::usage = "Partial[f, x[j]] gives the partial derivative " <>
	"of f with respect to x[j]."

Phi::usage = "Phi[z] is the modified inversion defined in " <>
	"Chapter 7 of Harmonic Function Theory."

PoissonKernel::usage =
	"PoissonKernel[x, z] gives the Poisson kernel for the unit ball."

PoissonKernelH::usage = "PoissonKernelH[x, y, t] gives the " <>
	"Poisson kernel for the upper half-space."

S::usage = "S is the south pole."

Schwarz::usage = "Schwarz[x] gives the maximum of |u[x]|, where " <>
	"u ranges over all harmonic functions on the unit ball with " <>
	"u[0] = 0 and |u| < 1."

SetDimension::usage = "SetDimension[x, n] sets the Dimension of x to n."

Singularity::usage = 
	"Singularity is an option for AntiLaplacian.  " <>
	"The default value is None."

Sphere::usage =
	"Sphere is one of the values that may be assigned to the option " <>
	"Orthonormal."

Subscripts::usage = "TurnOff[ Subscripts ] will force coordinates to " <>
	"print in the form x[j]."

SurfaceArea::usage = "SurfaceArea[n] gives the surface area of " <>
	"the unit ball in n dimensional real Euclidean space."

Taylor::usage = "Taylor[f, d, x] gives the Taylor expansion of f " <>
	"to degree d, with respect to x."

Togetherness::usage = "TurnOff[ Togetherness ] turns off the feature " <>
	"that simplifies some output.  TurnOn[ Togetherness ] turns " <>
	"this feature back on."

Tracer::usage = "Tracer is the same as the built-in Mathematica " <>
    "function Trace, used for debugging."

TurnOff::usage = "TurnOff[F] turns feature F off."

TurnOn::usage = "TurnOn[F] turns feature F on."

Volume::usage = "Volume[n] gives the volume of the unit ball in " <>
	"n-dimensional real Euclidean space."

ZonalHarmonic::usage = "ZonalHarmonic[m, x, z] gives the " <>
	"zonal harmonic of degree m with pole x."
(*

Protection

*)
Begin["`private`"]

Unprotect[ Gradient, Im, IntegerQ, Jacobian, Limit, Transpose ]

protectedWords = Unprotect[ Abs, Dot, 
	IdentityMatrix, Integer, Power, Times ]
(*

Togetherness

*)
postPower = (# //. { (a_ b_)^m_?IntegerQ :> a^m b^m, 
	(1/a_)^m_?IntegerQ :> 1/a^m,
	Rational[1, a_]^m_?IntegerQ :> 1/a^m
	} )& ;

SetAttributes[ postTogether, HoldAll ]

postTogether[Apart[f_]] := Apart[f]

postTogether[Together[f_]] := Together[f]

postTogether[Simplify[f_]] := Simplify[f]

postTogether[Expand[f_]] := Expand[f]

postTogether[ Set[a_, b_] ] := ( a = postTogether[b] )

shorter[f_] := With[ {g = Together[f]},
	If[ StringLength[ ToString[ Format[ g, InputForm ] ] ] < 
	StringLength[ ToString[ Format[ f, InputForm ] ] ], g, f ] ]

shorter/: f_Symbol[ a___, Literal[shorter][x_], b___ ] :=
	f[a, x, b] /; MemberQ[ Attributes[f], HoldAll ] ||
	MemberQ[ Attributes[f], HoldFirst ] ||
	MemberQ[ Attributes[f], HoldRest ]

postTogether[f_] := Map[shorter, f, {0, -2}] //.
	Literal[shorter][j_] -> j

usersPost = If[ ToString[$Post] === "$Post", Identity, $Post ]

Togetherness::Off = "Togetherness has been turned off.  " <>
	"Use the command  TurnOn[Togetherness]  to turn it back on."

TurnOn[Togetherness] ^:= ($Post = Composition[ postTogether, postPower,
	usersPost ] ; togetherness = True ; )

TurnOff[Togetherness] ^:= ($Post = Composition[ postPower,
	usersPost ] ; togetherness = False ; )

TurnOn[Togetherness]
(*

Calculus on Rn

IdentityMatrix

*)
identityMatrix = IdentityMatrix ;

Remove[ IdentityMatrix ]

HFT`IdentityMatrix::usage =
	"IdentityMatrix[n] is the n-by-n identity matrix."

HFT`IdentityMatrix[ a_?NumberQ ] := identityMatrix[a]

matrix[IdentityMatrix[_]] ^= True ;

IdentityMatrix/: a_ . IdentityMatrix[_] := a

IdentityMatrix/: IdentityMatrix[_] . v_ := v
(*

Transpose

*)
transpose = Transpose ;

Remove[ Transpose ]

HFT`Transpose::usage =
	"Transpose[A] is the transpose of a matrix A."

HFT`Transpose[ a_?MatrixQ ] := transpose[a]

matrix[ Transpose[_] ] = True

Transpose[ Transpose[ x_ ] ] := x

Transpose[{c_ x_?vector}] := c Transpose[{x}]

Transpose[ c_ x_?matrix ] := c Transpose[x]

Transpose[ IdentityMatrix[n_] ] := IdentityMatrix[n]

Transpose[{0}] := {0}
(*

Trace

*)
SetAttributes[ Tracer, HoldAll ]

Tracer = Trace

Attributes[ Tracer ] = Attributes[ Trace ]

options = Options[Trace]

Unprotect[ Trace ]

Remove[ Trace ]

Options[ Tracer ] = options

HFT`Trace::usage =
	"Trace[A] is the trace of a matrix A."

Trace [ x_ ] := Sum[ x[[j,j]], 
	{j, Length[x]} ] /; MatrixQ[x] && Length[x] == Length[ x[[1]] ]

Trace[ IdentityMatrix[n_] ] := n

Trace[ Transpose[A_] ] := Trace[A]

Trace[ c_ A_ ] := c Trace[A] /; matrix[A]
(*

Norm

*)
Norm[ {x_?vector, y_} ] := ( unmakeVector[y] ;
	Sqrt[ Norm[x]^2 + y^2 ] )

Norm[ x_ ] := 0 /; ( makeVector[x]; False )

Norm[x_List] := Sqrt[ x.x ]

Norm[0] = 0 ;

Norm[c_?NumberQ x_] := ( makeVector[x]; Abs[c] Norm[x] )

Norm[c_ x_?vector] := ( Im[c] = 0; Abs[c] Norm[x] )

Format[ Norm[x_] ] := SequenceForm[ "|", x, "|"]
(*

ExpandNorm

*)
ExpandNorm[f_] := 
	f //. Norm[g_ + h_] :> Sqrt[Norm[g]^2 + Norm[h]^2 + 2 g . h]
(*

Abs

*)
Abs[Norm[x_]^p_.] := ( makeVector[x]; Norm[x]^p )

Abs/: Abs[ x_[j_] ] ^ n_?EvenQ := x[j]^n /; vector[x]

Abs/: Abs[c_]^p_?EvenQ := c^p /; Im[c] == 0

Abs/: Abs[c_^p_?EvenQ] := c^p /; Im[c] == 0
(*

vector

*)
vector[ Literal[Plus[a__]] ] := If[ TrueQ[ Or @@ Thread[ vector[{a}] ] ],
	makeVector[ Plus[a] ] ; True ]

vector[ Literal[Times[a__]] ] := Or @@ Thread[ vector[{a}] ]

vector[A_?matrix . x_?vector] := True

vector[x_?vector . Y_?matrix] := True

vector[ a_[ ".", _ ] ] := True

vector[ a_[ _, "." ] ] := True

0[j_] := 0

vector[ Partial[_][b_][x_] ] := True /; vector[ b[x] ]

vector[ Gradient[_][_] ] = True
(*

makeVector

*)
vectors = {} ;

vectorValuedFunctions = {} ;

makeVector[x_Symbol] := If[ MemberQ[ vectors, x ], , If[ useSubscripts,
	Format[ x[j_] ] := Subscripted[ x[j] ] ] ;
	vectors = Append[ vectors, x ]; 
	vector[x] ^= True ;  ]

makeVector[c_?NumberQ x_] := makeVector[x]

makeVector[ Literal[ Plus[a__] ] ] := makeVector /@ {a}

makeVector[f_Symbol[_]] := If[ MemberQ[ vectorValuedFunctions, f ], ,
If[ useSubscripts, Format[ f[x_][j_] ] := Subscripted[ f[j] ][x] ] ;
vectorValuedFunctions = Append[ vectorValuedFunctions, f ]; 
vector[f[_]] ^= True ; vector[ Derivative[_][f][_] ] = True ; 
vector[ Partial[_][f][_] ] = True ; ] /; f =!= List

makeVector[ Partial[_][u_][_] ] := makeVector[ u[_] ]

Attributes[unmakeVector] = {Listable}

unmakeVector[x_Symbol] := If[ MemberQ[ vectors, x ],
	If[ useSubscripts, Format[ x[j$_] ] = . ] ;
	vectors = Delete[ vectors, Position[ vectors, x ][[1,1]] ] ;
	vector[x] ^= False ; ]
(*

matrix

*)
matrix[ Literal[Times[a__]] ] := Or @@ Thread[ matrix[{a}] ]

matrix[ Literal[Plus[a__]] ] :=  Or @@ Thread[ matrix[{a}] ]

matrix[a_?matrix . x_?matrix] := True

matrix[ Jacobian[_][_] ] = True

matrix[ {_} ] = True

0[_, _] := 0

matrix[ Transpose[x_] ] := matrix[x]
(*

makeMatrix

*)
matrices = {} ;

makeMatrix[x_Symbol] := If[ TrueQ[matrix[x]], , If[ useSubscripts, 
		Format[ x[j_, k_] ] := Subscripted[ x[j, k] ] ] ;
	matrices = Append[ matrices, x ]; 
	matrix[x] ^= True ; vector[x] ^= False ]
(*

Subscripts

*)
TurnOff[ Subscripts ] ^:= ( If[ useSubscripts, 
(Format[ #[j$_] ] = .)& /@ vectors ; 
(Format[ #[x$_][j$_] ] = .)& /@ vectorValuedFunctions ; 
(Format[ #[j$_, k$_] ] = .)& /@ matrices ;
Format[ Partial[{j__}] ] = . ] ; 
Format[ HilbertSchmidt[A_] ] = . ;
useSubscripts = False ; )

TurnOn[ Subscripts ] ^:= If[ TrueQ[ useSubscripts ], ,
(Format[ #[j$_] ] := Subscripted[ #[j$] ])& /@ vectors ; 
(Format[ #[x$_][j$_] ] := 
	Subscripted[ #[j$] ][x$])& /@ vectorValuedFunctions ; 
(Format[ #[j$_, k$_] ] := Subscripted[ #[j$, k$] ])& /@ matrices ;
Format[ Partial[{j__}] ] := HoldForm[ Subscripted[ Partial[j] ] ] ;
Format[ HilbertSchmidt[A_] ] := HoldForm[
	SequenceForm[ "|", A, "|", Subscript[2] ] ] ;
useSubscripts = True ; ]

TurnOn[ Subscripts ]
(*

Dot

*)
Attributes[ Dot ] = {OneIdentity}

b_ . 0 = 0 ;

0 . b_ = 0 ;

Literal[ Plus[x__] ] . z_ := Plus @@ ( (# . z)& /@ {x} )

a_ . Literal[ Plus[x__] ] := Plus @@ ( (a . #)& /@ {x} )

(c_ x_) . y_ := c (x.y) /; vector[x] || matrix[x] || NumberQ[c]

{c_ x_} . y_ := c ({x}.y) /; vector[x] || matrix[x] || NumberQ[c]

a_ . (c_ y_) := c (a.y) /; vector[y] || matrix[y] || NumberQ[c]

a_ . {c_ y_} := c (a.{y}) /; vector[y] || matrix[y] || NumberQ[c]

x_ . x_ := Norm[x]^2 /; matrix[x] =!= True && Head[x] =!= List 

Delta[j_] . Delta[k_] := Delta[j, k]

Delta[j_] . x_Symbol := ( makeVector[x]; x[j] ) /;
	matrix[x] =!= True

a_Symbol . Delta[j_] := ( makeVector[a]; a[j] ) /;
	matrix[a] =!= True

Delta[j_] . f_Symbol[x_] := ( makeVector[f[_]]; f[x][j] ) /;
	matrix[f[x]] =!= True /; f =!= List

b_Symbol[x_] . Delta[j_] := ( makeVector[b[_]]; b[x][j] ) /;
	matrix[b[x]] =!= True /; b =!= List

Delta[j_] . Gradient[u_][x_] := Partial[{j}][u][x]

Jacobian[b_][x_] . Delta[j_] := Partial[{j}][b][x]

Delta[j_] . Jacobian[b_][x_] := Gradient[ b . Delta[j] ][x]

Delta[j_] . Y_?matrix := Y[ j, "." ]

A_?matrix . Delta[j_] := A[ ".", j ]

Delta[k_] . Y[ j_, "." ] := Y[j, k]

A[ j_, "." ] . Delta[k_] := A[j, k]

Delta[j_] . Y[ ".", k_ ] := Y[j, k]

A[ ".", k_ ] . Delta[j_] := A[j, k]

Delta[j_] . Delta[k_] := Delta[j, k]

a_?VectorQ . Delta[j_] := a[[j]]

Delta[j_] . v_?VectorQ := v[[j]]

A_?MatrixQ . Delta[j_] := Transpose[A][[j]]

Delta[j_] . Y_?MatrixQ := A[[j]]

x_ . A_?identityMatrixMultiple := A[[1,1]] x /; primitiveVector[x]

A_?identityMatrixMultiple . x_ := A[[1,1]] x /; primitiveVector[x]

identityMatrixMultiple[A_] := MatrixQ[A] && ( Length[A] > 1 ) &&
	Length[A] === Length[ A[[1]] ] &&
	( And @@ ( (A[[#,#]] === A[[1,1]] &) /@ Range[ Length[A] ] ) ) &&
	( And @@ Flatten[ Table[ A[[m,n]] === 0 
		|| m === n, {m, Length[A]}, {n, Length[A]} ] ] )

a_List . x_Symbol := ( makeVector[x]; 
	a . Array[ x, Length[a] ] ) /; matrix[x] =!= True

x_Symbol . y_List := ( makeVector[x];
    Array[ x, Length[y] ] . y ) /; matrix[x] =!= True

a_List . B_Symbol := 
	a . Array[ B, Dimension[B] ] /; matrix[B]

B_Symbol . c_List := 
	Array[ B, Dimension[B] ] . c /; matrix[B]

x_?primitiveVector . y_?primitiveVector :=
	y . x /; ( If[ matrix[y] =!= True, makeVector[y] ];
	If[ matrix[x] =!= True, makeVector[x] ];
	Order[ ToString[y],ToString[x] ] == 1 &&
	matrix[y] =!= True && matrix[x] =!= True )

primitiveVector[x_Symbol] := True /; Context[x] =!= "System`"

primitiveVector[ Gradient[_][_] ] = True

primitiveVector[ Delta[_] ] = True

Format[ Power[w_ . z_, p_] ] := SequenceForm["(", w . z, ")"] ^ p

Format[ a_Dot . b_Dot ] := SequenceForm[ "(", a, ").(", b, ")" ]

Format[ a_Dot . b_ ]  := SequenceForm[
	"(", a, ")", ".", b]

Format[ a_ . b_Dot ] := SequenceForm[
	a, ".(", b, ")" ]

Format[ x_ . y_ ] := SequenceForm[x, ".", y]
(*

degree

*)
degree[f_, z_] := 
	degree1[ Expand[f /. Norm[z] -> z[normReplacement] ], z ]

degree1[ 0, _ ] := -Infinity

degree1[ Literal[Plus[a__]], z_ ] := Max[ Thread[ degree1[{a}, z] ] ]

degree1[ f_, z_ ] := Exponent[ f /. z[_] -> z, z ]
(*

monomialQ

*)
   monomialQ[f_, x_Symbol] := PolynomialQ[f, 
   Cases[ Variables[x], x[_] ]] && 
   FreeQ[f, Plus] && (Head[f] === Times && !Apply[Or, (FreeQ[f[[#1]],
   x] & ) /@ Range[Length[f]]] || Head[f] === Power || Head[f] ===  x)
(*

polynomialDecomposition

*)
polynomialDecomposition[ 0, x_ ] := Null

polynomialDecomposition[ p_, x_ ] := With[ {p1 = 
	Expand[ p /. Norm[x] -> x[normReplacement] ]},
	With[ {t = Map[ {#, degree[#,x]}&, 
	If[ Head[p1] === Plus, Apply[ List, p1 ], {p1} ] ]},  
	With[ {s = Table[ If[ MemberQ[t, {_,j}], 
	{Apply[ Plus, Transpose[ Cases[ t, {_,j}] ][[1]] ], j},
	Null ]
	, {j, 0, Max[Transpose[t][[2]]]} ]}, Delete[ s, Position[s, Null] ] ]
	] ] /. x[normReplacement] -> Norm[x]
(*

Delta

*)
vector[Delta[_]] ^:= True

Delta[j_, k_] := 0 /; j != k

Delta[j_, k_] := 1 /; j == k

Delta[j_][k_] := Delta[j, k]

Norm[Delta[_]] ^= 1 ;
(*

Partial

*)
Partial[f_, {x_[_], 0}] := ( makeVector[x]; f )

Partial[f_, {y_Symbol, 0}] := ( unmakeVector[y]; f )

Partial[f_, {z_, n_Integer}] :=
	Partial[ postTogether[ Partial[f, {z, n-1}] ], z ] /; Positive[n]

Partial[f_, z1_, z__] := Partial[ Partial[f, z1], z ]

Partial[f_, x_[_]] := ( makeVector[x]; 0 ) /; FreeQ[f, x]

Partial[f_, y_Symbol] := ( unmakeVector[y]; 0 /; FreeQ[f, y] )

Partial[x_[j_], x_[k_]] := ( makeVector[x]; Delta[j,k] )

Partial[y_, y_] := 1

Partial[x_, x_[j_]] := ( makeVector[x]; Delta[j] )

Partial[ Literal[Plus[a__]], y_ ] := Plus @@
	Thread[ Partial[{a}, y] ]

Partial[f_ g_, y_] := f Partial[g, y] + Partial[f, y] g

Partial[f_ . g_, y_] := ( If[ TrueQ[ matrix[f] ], , makeVector[f] ] ;
	If[ TrueQ[ matrix[g] ], , makeVector[g] ] ;
	f . Partial[g, y] + Partial[f, y] . g )

Partial[ f_[v_List], t_ ] := Plus @@ 
	((Partial[{#}][f][v] Partial[v[[#]], t]&) /@
	Range[ Length[v] ] )

Partial[f_^p_, y_] := p f^(p-1) Partial[f, y] +
							f^p Log[f] Partial[p, y]

Partial[f_[g_], x_[j_]] := ( makeVector[x];
		If[ TrueQ[ vector[g] ],
		    If[ TrueQ[ vector[f[g]] ],
		   	  Jacobian[f][g] . Partial[g, x[j]],
		   	  Gradient[f][g] . Partial[g, x[j]]],
		    Derivative[1][f][g] Partial[g, x[j]] ] )

Partial[f_[g_], y_Symbol] := ( unmakeVector[y] ;
		If[ TrueQ[ vector[g] ],
		    If[ TrueQ[ vector[f[g]] ],
		   	  Jacobian[f][g] . Partial[g, y],
		   	  Gradient[f][g] . Partial[g, y]],
		    Derivative[1][f][g] Partial[g, y] ] )

Partial[j_List][Partial[k_List][f_]] := Partial[Sort[Join[k, j]]][f]
(*

NormalD

*)
NormalD[f_, z_Symbol] := 
	Gradient[ f, z ] . z /. Norm[z] -> 1

NormalD[f_, v_List] := ( unmakeVector[v] ;
	Together[Gradient[ f, v ] . v] /. Norm[v]^2 -> 1 )
(*

Gradient

*)
Gradient[f_, x_Symbol] := ( makeVector[x]; 0 ) /; FreeQ[f, x]

Gradient[x_[j_], x_Symbol] := ( makeVector[x]; Delta[j] )

Gradient[ Literal[Plus[a__]], x_Symbol ] := Plus @@
	Thread[ Gradient[{a}, x] ]

Gradient[f_ g_, x_Symbol] := f Gradient[g, x] + Gradient[f, x] g

Gradient[f_ . g_, x_Symbol] := f . Jacobian[g, x] + 
					g . Jacobian[f, x]

Gradient[f_^p_, x_Symbol] := p f^(p-1) Gradient[f, x] +
							f^p Log[f] Gradient[p, x]

Gradient[f_[y_], x_Symbol] := ( makeVector[x];
			If[TrueQ[vector[y]],
				Gradient[f][y] . Jacobian[y,x],
				Derivative[1][f][y] Gradient[y, x] ] )

Gradient[Norm] := #/Norm[#] &

Gradient[Dimension] := 0 &

Gradient[f_, v_List] := ( unmakeVector[v] ;
	Partial[ f, # ] & /@ v )
(*

Divergence

*)
Divergence[f_, x_Symbol] := ( makeVector[x]; 0 ) /; FreeQ[f, x]

Divergence[x_, x_Symbol] := ( makeVector[x]; Dimension[x] )

Divergence[ Literal[Plus[a__]], x_Symbol ] := Plus @@
	Thread[ Divergence[{a}, x] ]

Divergence[ f_ g_, x_Symbol] :=
	f Divergence[g, x] + Gradient[f, x].g   /;
	( makeVector[x]; vector[g] )

Divergence[ A_?matrix . f_, x_Symbol ] := 
	Trace[ A . Jacobian[f, x] ] /; FreeQ[A, x]

Divergence[ f_ . T_?matrix, x_Symbol ] :=
	Trace[ Transpose[T] . Jacobian[f, x] ] /; FreeQ[T, x]

Divergence[ u_List, v_List ] := ( unmakeVector[v] ;
	Plus @@ Thread[ partial[ u,v ] ] /. partial -> Partial )
(*

Jacobian

*)
Jacobian[f_, x_Symbol] := ( makeVector[f]; makeVector[x]; 
	0 ) /; FreeQ[f, x]

Jacobian[x_, x_Symbol] := ( makeVector[x]; IdentityMatrix[ Dimension[x] ] )

Jacobian[ Literal[Plus[a__]], x_Symbol ] := Plus @@
	Thread[ Jacobian[{a}, x] ]

Jacobian[f_ g_?vector, x_Symbol] := f Jacobian[g, x] + 
	Transpose[ {Gradient[f, x]} ] . {g}

Jacobian[ f_Symbol[x_], x_Symbol ] := ( makeVector[x];
	makeVector[ f[_] ] ; Jacobian[f][x] ) /; f =!= List

Jacobian[ f_[y_] , x_Symbol ] := If[ TrueQ[ vector[y] ],
	Jacobian[f][y] . Jacobian[y, x],
	Transpose[{Derivative[1][f][y]}] . {Gradient[y, x]} ]

Jacobian[ A_?matrix . f_ , x_Symbol ] := ( makeVector[f];
	A . Jacobian[f, x] ) /; FreeQ[ A, x ]

Jacobian[ f_ . T_?matrix , x_Symbol ] := ( makeVector[f];
	Transpose[T] . Jacobian[f, x] ) /; FreeQ[ T, x ]

Jacobian[ f_List, v_List ] := ( unmakeVector[v] ;
	Table[ Partial[ f[[m]], v[[n]] ], 
	{m, Length[f]}, {n, Length[v]} ] )
(*

HilbertSchmidt

*)
HilbertSchmidt[A_?MatrixQ] := Plus @@ (Plus @@ (A^2))
(*

Laplacian

*)
Laplacian[f_, v_List] := ( unmakeVector[v] ;
	Plus @@ ( Partial[f, {#, 2}] & /@ v ) )

Laplacian[f_, x_Symbol] := ( makeVector[x]; 0 ) /; FreeQ[f, x]

Laplacian[x_[j_], x_Symbol] := ( makeVector[x]; 0 )

Laplacian[ Literal[Plus[a__]], x_Symbol ] := Plus @@
	Thread[ Laplacian[{a}, x] ]

Laplacian[f_ g_, x_Symbol] := f Laplacian[g, x] +
							Laplacian[f, x] g +
							2 Gradient[f, x] . Gradient[g, x]

Laplacian[f_^p_, x_Symbol] := p f^(p-1) Laplacian[f, x] +
		p (p-1) f^(p-2) Norm[Gradient[f, x]]^2 +
		2 p f^(p-1) Log[f] Gradient[f, x] . Gradient[p, x] +
		f^p Log[f] Laplacian[p, x] +
		2 f^(p-1) Gradient[f, x] . Gradient[p, x] +
		f^p Log[f]^2 Norm[Gradient[p,x]]^2

Laplacian[f_[a_. x_ + b_.], x_Symbol] := ( makeVector[x];
	If[ b == 0, , makeVector[b] ] ;
	a^2 Laplacian[f][a x + b] )/; FreeQ[a, x] && FreeQ[b, x]

Laplacian[a_ . x_, x_Symbol] := ( makeVector[x]; 0) /;
	FreeQ[a, x]

Laplacian[x_ . y_, x_Symbol] := ( makeVector[x]; 0) /;
	FreeQ[y, x]

Laplacian[ (A_ . x_) . y_, x_Symbol ] := (makeVector[x]; 0) /;
	FreeQ[A,x] && FreeQ[y,x]

Laplacian[ (x_ . Y_) . z_, x_Symbol ] := (makeVector[x]; 0) /;
	FreeQ[Y,x] && FreeQ[z,x]

Laplacian[ a_ . (B_ . x_), x_Symbol ] := (makeVector[x]; 0) /;
	FreeQ[a,x] && FreeQ[B,x]

Laplacian[ a_ . (x_ . Y_), x_Symbol ] := (makeVector[x]; 0) /;
	FreeQ[a,x] && FreeQ[Y,x]

Laplacian[ (A_ . x_) . x_, x_Symbol ] := (makeVector[x];
	2 Trace[A]) /; FreeQ[A,x]

Laplacian[ (x_ . A_) . x_, x_Symbol ] := (makeVector[x];
	2 Trace[A]) /; FreeQ[A,x]

Laplacian[ x_ . (A_ . x_), x_Symbol ] := (makeVector[x];
	2 Trace[A]) /; FreeQ[A,x]

Laplacian[ x_ . (x_ . A_), x_Symbol ] := (makeVector[x];
	2 Trace[A]) /; FreeQ[A,x]

Laplacian[ b_ . Gradient[u_][b_], b_Symbol ] := (makeVector[b] ;
	2 Laplacian[ u[b], b ] + b . Gradient[ Laplacian[ u[b], b ], b ] )

Laplacian[ Gradient[u_][x_] . x_, x_Symbol ] := (makeVector[x] ;
	2 Laplacian[ u[x], x ] + x . Gradient[ Laplacian[ u[x], x ], x ] )

Laplacian[f_[g_], x_Symbol] := 
		   Derivative[2][f][g] Norm[Gradient[g, x]]^2 +
		      Derivative[1][f][g] Laplacian[g, x] /;
		   ( makeVector[x]; vector[g] =!= True )   

Laplacian[Norm] := (Dimension[#] - 1)/Norm[#] &

Laplacian[ Norm[A_ . x_], x_ ] := ( makeVector[x];
	( Norm[A.x]^2 HilbertSchmidt[A]^2 -
	Norm[A.x.A]^2 ) / Norm[A.x]^3 ) /; FreeQ[A, x ]

Laplacian[ Norm[x_ . Z_], x_ ] := ( makeVector[x];
	( Norm[x.Z]^2 HilbertSchmidt[Z]^2 -
	Norm[x.Z.Transpose[Z]]^2 ) / Norm[x.Z]^3 ) /; FreeQ[Z, x ]

Laplacian[Dimension] := 0 &

Laplacian[f_, x_Symbol, y_Symbol] := ( unmakeVector[y];
	Laplacian[f, x] + Partial[f, {y,2}] )

(Laplacian^n_)[ f_, x___ ] := (Laplacian^(n-1))[ Laplacian[ f, x ], x ] /;
	Head[n] == Integer && Positive[n]

Laplacian[ (Laplacian^n_.)[ f_ ] ] := (Laplacian^(n+1))[f]

(Laplacian^m_.)[ (Laplacian^n_.)[ u_ ] ] := (Laplacian^(m+n))[u]
(*

Dimension and SetDimension

*)
Dimension[a_ x_] := Dimension[x] /; vector[x] || matrix[x]

Dimension[a_?NumberQ x_] := ( If [matrix[x] =!= True, makeVector[x] ]; 
	Dimension[x] )

Dimension[ a_. S + x_ ] := Dimension[x]

Dimension[x_ + w_] := If[ StringLength[ToString[Dimension[x]]] < 
	StringLength[ToString[Dimension[w]]],
	Dimension[x], Dimension[w]]

Dimension[v_?VectorQ] := Length[v]

Dimension[A_?MatrixQ] := { Length[A], Length[ A[[1]] ] }

IntegerQ[ Dimension[_] ] ^:= True 

SetDimension[ x_Symbol, {m_, n_} ] := ( Dimension[x] = {m, n }; 
	makeMatrix[x] ; IntegerQ[m] = True ; IntegerQ[n] = True ; )

SetDimension[ x_Symbol, n_ ] := ( Dimension[x] = n ; makeVector[x] ;
	IntegerQ[n] = True ; )

SetDimension[ f_Symbol[_], n_ ] := ( Dimension[ f[_] ] = n ;
	makeVector[ f[_] ] ; IntegerQ[n] = True ; ) /; f =!= List
(*

Homogeneous

*)
Homogeneous[f_, d_, x_Symbol, c_Symbol] := ( makeVector[x];
	makeVector[c];
	SetDimension[ c, Dimension[x] ];
	Homogeneous[ f, d, Union[ Cases[ Level[ f, -1 ], x[_] ] ], 
	( Union[ Cases[ Level[ f, -1 ], x[_] ] ] /.
	x -> c ) ] /.
	Array[c, Length[Union[ Cases[ Level[ f, -1 ], x[_] ] ]]] -> c
	) /; Not[ 
	MemberQ[ Level[ f, {-1} ], x ] ] && Head[d] == Integer

Homogeneous[f_, d_, x_Symbol, c_Symbol] :=
	( makeVector[x];
	makeVector[c];
	SetDimension[ c, Dimension[x] ] ;
	Homogeneous[ f, d, x, Array[c, Dimension[x] ] ] /.
	Array[c, Dimension[x] ] -> c) /;
	MemberQ[ Level[ f, {-1} ], x ] && 
	Head[Dimension[x]] == Integer && Head[d] == Integer

Homogeneous[f_, d_, x_Symbol, c_List] := ( makeVector[x] ;
	unmakeVector[c] ;
	Module[ {z}, 
	Homogeneous[ f , d, Array[x, Length[c]], c ] /.
	x[j_] -> z[j] /. x -> c /. z[j_] -> x[j] ] ) /;
	Head[d] == Integer

Homogeneous[f_, d_, v_List, c_Symbol] := ( makeVector[c];
	unmakeVector[v] ;
	SetDimension[ c, Length[v] ] ;
	Homogeneous[ f, d, v, Array[c, Dimension[c]] ] ) /;
	Head[d] == Integer

Homogeneous[f_, d_, v_List, c_List] := 
  ( unmakeVector[v] ; unmakeVector[c] ;
  If[ togetherness && Not[Norm[c] === 0], 
  		TurnOff[ Togetherness ] ; Message[Togetherness::Off] ] ;
  With[{mI = multiIndices[d,Length[v]]},
	Plus @@ ( ( ( Apply[Partial, 
		Prepend[Transpose[{v, mI[[#]]}], f]]
		/. Thread[v -> c] )
		(Times @@ ((v - c)^mI[[#]])) /
		multiFactorial[mI[[#]]] ) & /@
		Range[ Length[mI] ] ) ] ) /; Head[d] == Integer

Homogeneous[f_, d_, v_List] := Homogeneous[ f, d, v,
	Table[0, {Length[v]}] ] /; Head[d] == Integer

Homogeneous[f_, d_, x_Symbol] := ( Homogeneous[ f, d, 
	Union[ Cases[ Level[ f, -1 ], x[_] ] ] ] /. Table[0,
	{ Length[Union[ Cases[ Level[ f, -1 ], x[_] ] ]] } ] -> 0 ) /; Not[ 
	MemberQ[ Level[ f, {-1} ], x ] ] && Head[d] == Integer

Homogeneous[f_, d_, x_Symbol] := 
	( Homogeneous[ f, d, x, Table[0, {Dimension[x]}] ] 
	/. Table[0,{ Dimension[x]} ] -> 0 ) /;
	MemberQ[ Level[ f, {-1} ], x ] &&
	Head[Dimension[x]] == Integer && Head[d] == Integer

multiFactorial[a_List] := Times @@ (a!)

multiIndices[k_, 1] := {{k}}

multiIndices[0, n_] := {Table[0, {n}]}

multiIndices[k_Integer, n_Integer] := multiIndices[k, n] =
	Apply[Join,
		Table[Map[Prepend[#,k-j]&, multiIndices[j,n-1]],
			{j,0,k}
		] ]
(*

Taylor

*)
Taylor[f_, d_, vc___] := 
	Plus @@ ( Homogeneous[f, #, vc] & /@ Range[0, d] ) /;
	Head[d] == Integer
(*

Volume

*)
Volume[n_Integer] := Pi^(n/2) / Gamma[n/2 + 1] /; Positive[n]
(*

SurfaceArea

*)
SurfaceArea[n_Integer] := 2 Pi^(n/2) / Gamma[n/2] /; Positive[n]
(*

IntegrateBall

*)
IntegrateBall[f_ , v_List] := ( unmakeVector[v] ;
	Module[{z}, SetDimension[ z, Length[v] ] ;
	IntegrateBall[ f /. 
	Thread[v -> Array[z, Length[v]] ], z ] ] )

IntegrateBall[f_, x_Symbol] := ( makeVector[x];
	f Volume[ Dimension[x] ] ) /; FreeQ[f, x]

IntegrateBall[f_, x_Symbol] := Integrate[f /. 
	{Norm[x]^2 -> x^2, Norm[x] -> Abs[x], x[1] -> x},
	{x, -1, 1}] /; Dimension[x] == 1

IntegrateBall[f_, x_Symbol] := ( makeVector[x];
	Module[ {r}, 
	( If[  Head[ Dimension[x]] === Integer, ,
	Limit[ c_. r^(Dimension[x] + m_.), r -> 0 ] := 0 
	/; ( (m > 0 || m == 0) && FreeQ[c, r] ) ] ;
	Dimension[x] Volume[ Dimension[x] ] *
	Integrate[ r^(Dimension[x]-1) *
	IntegrateSphere[f /. {Norm[x] -> r, x[j_] -> r x[j],
	x -> r x}, x ], {r, 0, 1} ] ) ] )
(*

IntegrateSphere

*)
IntegrateSphere[f_ , v_List] := ( unmakeVector[v];
	Module[{z}, SetDimension[ z, Length[v] ] ;
	IntegrateSphere[ f /. 
	Thread[v -> Array[z, Length[v]] ], z ] ] )

IntegrateSphere[f_, x_Symbol] := ( makeVector[x];
	integrateSphere1[ Expand[f
	/. { Norm[x] -> 1,
	x . z_ :> Array[ x, Dimension[x] ] . Array[ z, Dimension[x] ],
	z_ . x :> Array[ z, Dimension[x] ] . Array[ x, Dimension[x] ] }
	], x] )

integrateSphere1[ Literal[Plus[a__]], x_ ] := Plus @@
	Thread[ integrateSphere1[{a}, x] ]

integrateSphere1[f_, x_] := f /; FreeQ[f, x]

integrateSphere1[c_. x_[_], x_] := 0

integrateSphere1[c_ f_, x_] := 
	c integrateSphere1[f, x] /; FreeQ[c, x]

integrateSphere1[x_[_]^j_, x_] := integrateSphere3[ {j}, Dimension[x] ]

integrateSphere1[f_, x_] := integrateSphere2[ 
	f , x ] /; monomialSpecialQ[f,x]

integrateSphere2[f_, x_] := integrateSphere3[
	(List @@ f) /. x[_]^b_ -> b,
	Dimension[x] ]

integrateSphere3[b_, n_] := 0 /; Or @@ OddQ[b]

integrateSphere3[b_, n_] := 
	( Times @@ ((b - 1)!!) ) /
	Times @@ Range[ n, n + (Plus @@ b)  - 2, 2 ]

monomialSpecialQ[f_, x_Symbol] := 
	Head[f] === Times && 
		 And @@ ( MatchQ[ f[[#]], x[_]^_ ] & /@ Range[ Length[f] ] )
(*

Boundary Value Problems

Dirichlet

*)
Dirichlet[p_, x_Symbol] := 
	Plus @@ (#[[1]]& /@
	HarmonicDecomposition[Limit[ p /. Norm[x] -> r, r -> 1], x])

Dirichlet[f_, g_, z_Symbol] :=
	( makeVector[z];
	With[ { t = AntiLaplacian[g, z ]},
	t + Dirichlet[f - t, z] ] )

Dirichlet[f_, v_List] := explicit[ Dirichlet, f, v ]

Dirichlet[f_, g_, v_List] := explicit[ Dirichlet, f, g, v ]
(*

HarmonicDecomposition

*)
d[0,0] = 1

lap[ f_, v_ ] := Plus @@ (D[f, {#,2}]& /@ v)

HarmonicDecomposition[ 0, x_ ] := {}

HarmonicDecomposition[ p_, x_Symbol ] := ( makeVector[x];
	With[ {deg=degree[p /. Norm[x] -> norm,x],
	n = Dimension[x]},
	powers = {}; hd[_] := 0;
	With[ {partition = 
	Partition[ Flatten[ Expand[(
	m = #[[2]];
	k = Floor[m/2];
	If[m > 1, d[1,1] = 1/ (2 (2m+n-4) )];
	v = Select[ Variables[#[[1]]], Not[FreeQ[#,x]]& ];
	laplac[0] = #[[1]];
	laplac[1] = lap[ #[[1]], v ];
	Do[ d[j,j] = d[j-1,j-1] (2m + n - 2j) /
	( 2j (2m + n - 4j + 2)(2m + n - 4j) );
	laplac[j] = lap[ laplac[j-1], v ], {j, 2, k} ];
	Do[ d[i,j] = d[i,j-1]/ ( 2(i-j) ( 2(m-i-j-1)+n ) ),
	{i, 0, k - 1}, {j, i+1, k} ];
	Table[ {
	Sum[ d[i,j] Norm[x]^(2(j-i)) laplac[j], {j,i,k}]
	, 2i}, {i,0,k}] )& /@
	polynomialDecomposition[ p /. Norm[x] -> norm, x ]]
	//. {a_. + norm^e_. b_., j_} -> {{a,j}, {b, j+e}} ], 2 ]},
	powers = Union[ powers, Transpose[partition][[2]] ];
	(hd[#[[2]]] = hd[#[[2]]] + #[[1]])& /@
	partition ];
	aa = Select[ {hd[#],#}& /@ powers, (#[[1]] =!= 0)& ];
	Clear[hd]; aa ] )

HarmonicDecomposition[f_, v_List] := explicit[
	HarmonicDecomposition, f, v ]
(*

AntiLaplacian

*)
Options[ AntiLaplacian ] = { Multiple -> None, Singularity -> None }

AntiLaplacian[0, x_, options___?OptionQ] := 0

AntiLaplacian[p_, x_Symbol, options___?OptionQ] := ( makeVector[x]; 
	If[ ((Multiple  /. {options} /. Options[AntiLaplacian]) ===
	Norm^2) || Not[ FreeQ[p, Norm[x]] ],
	antiLaplacianN[ Expand[p], x, options ],
	antiLaplacian[ Expand[p], x ] ] )

antiLaplacianN[Literal[Plus[a__]], z_, options___?OptionQ ] := Plus @@
	Thread[ Unevaluated[ antiLaplacianN[{a}, z, options] ] ]

antiLaplacianN[ c_ g_, z_, options___?OptionQ ] := 
	c antiLaplacianN[g, z, options] /; FreeQ[c, z]

antiLaplacianN[ g_, z_, options___?OptionQ ] :=
	antiLap1[ If[ Head[g] === Times,
		Select[ g, FreeQ[#, Norm[z]]& ],
		If[ FreeQ[g, Norm[z]], g , 1] ],
	          If[ Head[g] === Times,
	    Select[ g, Not[FreeQ[#, Norm[z]]]& ],
	    If[ FreeQ[g, Norm[z]], 1 , g] ],
	 z, options ]

antiLap1[ p_, f_, z_, options___?OptionQ ] :=
	Plus @@ ( antiLap2[ #[[1]], #[[2]], f, z, options ]& /@
	HarmonicDecomposition[p, z] )

antiLap2[ p_, j_, f_, z_, options___?OptionQ ] :=
	Plus @@ ( antiLap3[ #[[1]], #[[2]], j, f, z, options ]& /@
	polynomialDecomposition[p, z] )

antiLap3[ p_, m_, j_, f_, z_, options___?OptionQ ] :=
	With[ {n = Dimension[z]},
	p *
	If[ ((Singularity  /. {options} /. Options[AntiLaplacian])
		=!= None),
		Integrate[ t^(1 - 2m - n) *
	    Integrate[t^(-1 + 2m + n + j) f /. Norm[z] -> t, t],
	    t ],
	With[ {int =
	Integrate[s^(-1 + 2m + n + j) f /. Norm[z] -> s, {s,0,t}]},
	Integrate[ t^(1 - 2m - n) *
	If[ FreeQ[int,Indeterminate] && FreeQ[int,Infinity] &&
		FreeQ[int,Limit] &&
		((Singularity  /. {options} /. Options[AntiLaplacian]) ===
	    None), int,
	Message[ AntiLaplacian::singularity ];
	Integrate[t^(-1 + 2m + n + j) f /. Norm[z] -> t, t] ],
	t ] ] ]  /. t -> Norm[z] ]

AntiLaplacian::singularity =
	"Possible singularity at 0.  HFT.m is switching to another method."

AntiLaplacian[f_ , v_List, options___?OptionQ] := ( unmakeVector[v];
	Module[{z}, SetDimension[z, Length[v] ] ;
	AntiLaplacian[ (f /. Norm[v]^2 -> Norm[z]^2) /. 
	Thread[v -> Array[z, Length[v]] ], z, options ]
	/. Prepend[Thread[Array[z, Length[v]] -> v ], 
	Norm[z] -> Norm[v] ] ] )

antiLaplacian[Literal[Plus[a__]], z_ ] := Plus @@
	Thread[ antiLaplacian[{a}, z] ]

antiLaplacian[g_, z_] := g z[1]^2/2 /; FreeQ[g, z]

antiLaplacian[g_, z_] := antiLaplacian3[
	Transpose[ Sort[ 
	({Exponent[ g, z[#] ], #} & /@ Map[ #[[1]]&, Variables[ g ] ] )
		] ], g, z ] /; 
	monomialQ[g,z]

antiLaplacian[c_ g_, z_] := 
	c antiLaplacian[g, z] /; FreeQ[c, z]

antiLaplacian3[{e_, p_}, g_, z_] := z[ p[[-1]] ]^(e[[-1]]+2) /
			( (e[[-1]]+2) (e[[-1]]+1) ) /; Length[e] == 1

antiLaplacian3[{e_, p_}, g_, z_] := 
	g z[ Last[p] ]^2 / ( (Last[e]+1)(Last[e]+2) ) -
	Plus @@
	( ( (e[[#]]-1) e[[#]]
	antiLaplacian[ z[ Last[p] ]^2 g / z[ p[[#]] ]^2, z ] /
	( (Last[e]+1)(Last[e]+2) ) ) & /@
	Range[ Position[Positive[e], True][[1,1]], Length[e] - 1 ] ) /;
	Length[e] > 1
(*

ExteriorDirichlet

*)
ExteriorDirichlet[f_, z_] := Factor[ 
	Kelvin[ Dirichlet[f, z], z ] ]

ExteriorDirichlet[f_, g_, z_Symbol] := Factor[
	With[ { t = AntiLaplacian[g, z, Singularity -> 0 ]},
	t + ExteriorDirichlet[f - t, z] ] ]

ExteriorDirichlet[f_, v_List] :=
	explicit[ ExteriorDirichlet, f, v ]

ExteriorDirichlet[f_, g_, v_List] :=
	explicit[ ExteriorDirichlet, f, g, v ]
(*

AnnularDirichlet

*)
annular[ 0, r_, s_, z_ ] = 0

annular[ f_, r_, s_, z_ ] := Plus @@ 
	( If[ #[[2]] === #[[3]],
	(Log[Norm[z]] - Log[s]) r^#[[2]] #[[1]] / 
	(Log[r] - Log[s]),
	(Norm[z]^(2#[[2]]-2#[[3]]) -
	s^(2#[[2]]-2#[[3]])) r^#[[2]] #[[1]] / 
	(r^(2#[[2]]-2#[[3]]) -
	s^(2#[[2]]-2#[[3]]))]& /@ 
	( Flatten[ append[ 
	HarmonicDecomposition[#[[1]], z], #[[2]]]& /@
	polynomialDecomposition[ f, z], 1 ] ) 
	) /; Dimension[z] == 2

annular[ f_, r_, s_, z_ ] := Plus @@ 
	( ( (Norm[z]^(2-Dimension[z]+2#[[2]]-2#[[3]]) -
	s^(2-Dimension[z]+2#[[2]]-2#[[3]])) r^#[[2]] #[[1]] / 
	(r^(2-Dimension[z]+2#[[2]]-2#[[3]]) -
	s^(2-Dimension[z]+2#[[2]]-2#[[3]])))& /@ 
	( Flatten[ append[ HarmonicDecomposition[#[[1]], z], #[[2]]]& /@
	polynomialDecomposition[ f, z], 1 ] ) ) 

append[ b_, d_ ] := Append[#,d]& /@ b

AnnularDirichlet[ f_, g_, r_, s_, z_Symbol ] := 
	annular[ f, r, s, z ] + annular[ g, s, r, z ]

AnnularDirichlet[ f_, g_, r0_, r1_, v_List ] := ( unmakeVector[v] ;
	Module[{x}, SetDimension[x, Length[v] ] ;
	AnnularDirichlet[ f /. Thread[v -> Array[x, Length[v]] ], 
	g /. Thread[v -> Array[x, Length[v]] ], r0, r1, x ]
	/. Thread[Array[x, Length[v]] -> v ] /. Norm[x] -> Norm[v] ] )

AnnularDirichlet[ f0_, f1_, g_, r0_, r1_, z_Symbol ] :=
	With[ { t = AntiLaplacian[g, z]},
	t + AnnularDirichlet[f0 - t /. Norm[z] -> r0,
	f1 - t/. Norm[z] -> r1, r0, r1, z] ]

AnnularDirichlet[f0_, f1_, g_, r0_, r1_, v_List] := (
	unmakeVector[v];
	Module[{x}, SetDimension[x, Length[v] ] ;
	AnnularDirichlet[
	f0 /. Thread[v -> Array[x, Length[v]] ], 
	f1 /. Thread[v -> Array[x, Length[v]] ],
	r0, r1,
	g /. Thread[v -> Array[x, Length[v]] ], x ]
	/. Thread[Array[x, Length[v]] -> v ] /. Norm[x] -> Norm[v] ] )
(*

Neumann

*)
Neumann[f_, v_List] := explicit[ Neumann, f, v ]

Neumann[f_, g_, v_List] := explicit[ Neumann, f, g, v ]

Neumann[0, x_Symbol] := 0

Neumann[f_, x_Symbol] := Apply[ Plus,
	Map[ If[ #[[2]] === 0, 0, #[[1]] / #[[2]] ]&, 
	polynomialDecomposition[ Dirichlet[ f, x ], x ] ] ]

Neumann[f_, g_, x_Symbol] := With[ {v = AntiLaplacian[g, x]},
	v + Neumann[ f - NormalD[v,x], x ] - Limit[ 
	v /. {x[_] -> 0, Norm[x] -> t}, t -> 0, Direction -> -1 ] ]
(*

ExteriorNeumann

*)
ExteriorNeumann[f_, v_List] := explicit[ ExteriorNeumann, f, v ]

ExteriorNeumann[f_, g_, v_List] := explicit[
	ExteriorNeumann, f, g, v ]

ExteriorNeumann[0, x_Symbol] := 0

ExteriorNeumann[f_, x_Symbol] := Apply[ Plus,
	Map[ Norm[x]^(2-Dimension[x]-2#[[2]])
	#[[1]] / (Dimension[x]+#[[2]]-2)&, 
	polynomialDecomposition[ Dirichlet[ f, x ], x ] ]
	] /; Dimension[x] =!= 2

ExteriorNeumann[f_, x_Symbol] := Apply[ Plus,
	Map[ Norm[x]^(-2#[[2]])
	#[[1]] If[#[[2]]===0, 0, 1/#[[2]]]&, 
	polynomialDecomposition[ Dirichlet[ f, x ], x ] ] 
	] /; Dimension[x] === 2

ExteriorNeumann[f_, g_, x_Symbol] := With[
	{v = AntiLaplacian[g, x, Singularity -> 0]},
	v + ExteriorNeumann[ f + NormalD[v,x], x ] ]
(*

BiDirichlet

*)
BiDirichlet[f_, v_List] := explicit[ BiDirichlet, f, v ]

BiDirichlet[f_, g_, v_List] := explicit[ BiDirichlet, f, g, v ]

BiDirichlet[ f_, z_Symbol ] := Together[ 
	With[ {p = Dirichlet[ f, z ]},
	p + (1 - Norm[z]^2) z . Gradient[ p, z] / 2 ] ]

BiDirichlet[ f_, g_, z_Symbol ] := 
	normalIntegral[ Dirichlet[g,z], z ] + BiDirichlet[ f, z ]

normalIntegral[p_, z_] :=
	( Norm[z]^2 - 1) p / 2
(*

Spherical Harmonics

ZonalHarmonic

*)
ZonalHarmonic[ 0, x_, z_ ] := 1

ZonalHarmonic[ m_Integer, x_Symbol, z_Symbol ] := 
	( makeVector[x] ; makeVector[z] ;
	Plus @@ 
	( (-1)^# (Dimension[x + z] + 2m - 2)
	( Times @@ 
	Range[ Dimension[x + z], Dimension[x + z] + 2m - 2# -4, 2 ] )
	(x . z)^(m - 2#) Norm[x]^(2#) Norm[z]^(2#)
	/ (#! (m - 2#)! 2^#) & /@ Range[0, m/2] ) )

ZonalHarmonic[ m_Integer, x_List, z_List ] := 
	( unmakeVector[x] ; unmakeVector[z] ;
	Expand[ 
	Plus @@ 
	( (-1)^# (Length[x] + 2m - 2)
	( Times @@ Range[ Length[x], Length[x] + 2m - 2# -4, 2 ] )
	(x . z)^(m - 2#) Norm[x]^(2#) Norm[z]^(2#)
	/ (#! (m - 2#)! 2^#) & /@ Range[0, m/2] ) ] )
(*

DimensionH

*)
DimensionH[0, 2] = 1 ;

DimensionH[m_, n_] := Binomial[n + m - 1, n - 1] - 
	Binomial[n + m - 3, n - 1]
(*

BasisH

*)
Options[ BasisH ] = { Orthonormal -> None } ;

BasisH[ 0, x_, options___?OptionQ ] := With[
	{orthonormal = Orthonormal  /. {options} /. Options[BasisH]},
	If[ orthonormal === None, {1},
	Which[ orthonormal === Sphere, {1},
	  orthonormal === Ball, {1/Sqrt[ Volume[ Dimension[x] ] ]} ] ] ]

BasisH[m_Integer, x_Symbol, options___?OptionQ ] := postTogether[ With[
	{ orthonormal = Orthonormal  /. {options} /. Options[BasisH] },
	  Module[ {xx}, With[ {ttt = Expand[
	Expand[ ( Kelvin[#, xx]& /@
	( (Partial @@ Prepend[
	Transpose[ {Array[xx, 2], #} ],
	Log[Norm[xx]]]) & /@
	Join[ Prepend[#, 0]& /@ multiIndices[m, 1],
		Prepend[#, 1]& /@ multiIndices[m-1, 1] ] ) ) ] ]
	/. xx -> x } ,
	If[ orthonormal === None, 
	  noLeadingMinus /@
	  (Transpose[ FactorTermsList /@ ttt ][[2]] ),
	  orthonormalize[ ttt, (Which[
	  orthonormal === Sphere, IntegrateSphere,
	  orthonormal === Ball, IntegrateBall
	  ]
	  [#1 #2, x ] &) ]
	   ] ] ] ] ] /; Dimension[x] == 2

BasisH[m_Integer, x_Symbol, options___?OptionQ ] :=
	postTogether[ With[
	{ orthonormal = Orthonormal  /. {options} /. Options[BasisH] },
	  Module[ {xx}, With[ {ttt = Expand[
	Expand[ ( Kelvin[#, xx]& /@
	( (Partial @@ Prepend[
	Transpose[ {Array[xx, Dimension[x]], #} ],
	Norm[xx]^(2 - Dimension[x])]) & /@
	Join[ Prepend[#, 0]& /@ multiIndices[m, Dimension[x]-1],
		Prepend[#, 1]& /@ multiIndices[m-1, Dimension[x]-1] ] ) ) ] ] 
		/. xx -> x },
	If[ orthonormal === None, 
	  noLeadingMinus /@
	  (Transpose[ FactorTermsList /@ ttt ][[2]] ),
	  orthonormalize[ ttt, (Which[
	  orthonormal === Sphere, IntegrateSphere,
	  orthonormal === Ball, IntegrateBall
	  ] [#1 #2, x ] &) ] ] ] ] ] ] /;
	Head[Dimension[x]] == Integer

BasisH[m_Integer, v_List, options___?OptionQ ] := ( unmakeVector[v] ;
  Module[ {z}, SetDimension[ z, Length[v] ] ;
  Expand[ BasisH[ m, z, options ] /. Thread[Array[z, Length[v]] -> v ] ]
    /. Norm[z] -> Norm[v]] )

noLeadingMinus[x_] := If[ Head[x] === Plus,
		 	If[Head[ x[[1]] ] === Times,
				If[ x[[1, 1]] < 0, Expand[- x], x, x],
				If[ x[[1]] < 0, Expand[- x], x, x] ],
			x ]		
(*

orthonormalize

*)
orthonormalize[ b_List, innerProduct_ ] := Module[ {e, j, t},
  e[1] = b[[1]] / Sqrt[ innerProduct[b[[1]],b[[1]]] ];
  Do[
  t = b[[j]] -
   ( Plus @@ ( innerProduct[ b[[j]], e[#] ] e[#] & /@ Range[j-1] ) ) ;
  e[j] = t / Sqrt[ innerProduct[t,t] ],
  {j, 2, Length[b]} ];
  e[#]& /@ Range[Length[b]] ]
(*

The Kelvin Transform

Inversion

*)
Inversion[Infinity] := 0

Inversion[x_] := ( If[ Head[x] == Symbol, makeVector[x] ];
	If[Norm[x] === 0, Infinity, x/Norm[x]^2] )

Inversion[x_, y_] := ( If[ Head[x] == Symbol, makeVector[x] ];
	If[Norm[x] === 0 && y === 0, Infinity, 
	( x - y S ) / ( Norm[x]^2 + y^2)] )
(*

Kelvin

*)
Kelvin[ Literal[Plus[u__]], x_ ] := Plus @@
	Thread[ Kelvin[{u}, x] ]

Kelvin[u_, x_List] := ( unmakeVector[x] ;
	Norm[x]^(2 - Length[x]) * (u /. Thread[x -> Inversion[x]]) )

Kelvin[u_, x_Symbol] := ( makeVector[x] ;
	Norm[x]^(2 - Dimension[x])
	(u /. {x -> Inversion[x], x[j_] :> x[j]/Norm[x]^2}) )
(*

Phi

*)
Phi[ {x_?vector, y_} ] := Phi[x, y]

Phi[z_List] := ( unmakeVector[z];
	If[ Last[z] === -1 && Norm[ Drop[z, -1] ] === 0, Infinity,
	2 Inversion[ z + Append[ Table[0, {Length[z]-1}], 1 ] ] -
	Append[ Table[0, {Length[z]-1}], 1 ] ] )

Phi[z_] := ( If[ Not[ z === 0 || z === Infinity || MatchQ[ z, _. S ] ],
	SetDimension[S, Dimension[z] ], makeVector[z] ] ;
	2 Inversion[z - S] + S )

Phi[ 0, -1 ] = Infinity ;

Phi[x_, y_] := ( If[ x =!= 0, 
	SetDimension[S, Dimension[x] + 1]; makeVector[x];
	unmakeVector[y] ] ; 
	{ 2 x / ( Norm[x]^2 + (y + 1)^2 ), 
	2 (y + 1) / ( Norm[x]^2 + (y + 1)^2 ) - 1 } )
(*

S

*)
Norm[S] ^= 1 ;

S . z_Symbol := -z[Dimension[z]] /; matrix[z] =!= True

a_Symbol . S := -a[Dimension[a]] /; matrix[z] =!= True

makeVector[S] ;

S[j_] := -1 /; j == Dimension[S]

S[j_] := 0 /; j < Dimension[S]
(*

KelvinM

*)
KelvinM[ Literal[Plus[u__]], x_ ] := Plus @@
	Thread[ KelvinM[{u}, x] ]

KelvinM[u_, x_List] := ( unmakeVector[x] ;
	2^( (Length[x]-2)/2 ) *
	(Norm[x]^2 + 2 Last[x] + 1) ^ (2 - Length[x]) *
	(u /. Thread[x -> Phi[x]]) )

KelvinM[u_, z_Symbol] := ( makeVector[z] ;
	2^( (Dimension[z]-2)/2 ) *
	Norm[ z - S ]^(2 - Dimension[z]) *
	(u /. {z -> Phi[z], z[j_] :> Phi[z].Delta[j]}) )

KelvinM[u_, x_Symbol, y_Symbol] := ( makeVector[x] ;
	unmakeVector[y] ;
	2^( (Dimension[x]-1)/2 ) *
	(Norm[x]^2 + (y+1)^2)^((1 - Dimension[x])/2) *
	( u /. {x -> 2 x / ( Norm[x]^2 + (y + 1)^2 ),
	y -> 2 (y + 1) / ( Norm[x]^2 + (y + 1)^2 ) - 1,
	x[j_] :> 2 x[j] / ( Norm[x]^2 + (y + 1)^2 ) } ) )
(*

Kernels

PoissonKernel

*)
PoissonKernel[ x_, z_ ] := ( makeVector[x]; makeVector[z]; 
  (1 - Norm[x]^2 Norm[z]^2)/
  (1 - 2 x.z + Norm[x]^2 Norm[z]^2)^(Dimension[x+z]/2) )
(*

PoissonKernelH

*)
PoissonKernelH[z_, w_] := ( makeVector[z]; makeVector[w];
	2/( Dimension[z] Volume[ Dimension[z] ] ) *
	( z.Delta[ Dimension[z] ] + w.Delta[ Dimension[z] ] ) /
	( Norm[z]^2 + Norm[w]^2 - 2 z.w + 
	4 z.Delta[ Dimension[z] ] w.Delta[ Dimension[z] ] )^
	( Dimension[z]/2 ) )

PoissonKernelH[ x_, y_, t_, u_] := ( makeVector[x] ; makeVector[t] ;
	unmakeVector[y] ; unmakeVector[u] ;
	2/( (Dimension[x]+1) Volume[Dimension[x]+1] ) *
	(y + u) / ( Norm[x-t]^2 + (y+u)^2 )^( (Dimension[x]+1)/2 ) )
(*

BergmanKernel

*)
BergmanKernel[x_, y_] := ( makeVector[x]; makeVector[y];
	With[ {n = Dimension[x + y]},
	( n + (8 x.y - 2 n - 4) Norm[x]^2 Norm[y]^2 +
	(n-4) Norm[x]^4 Norm[y]^4 ) /
	( n Volume[n] (1 - 2 x.y + Norm[x]^2 Norm[y]^2)^(1 + n/2) ) ] )
(*

BergmanKernelH

*)
BergmanKernelH[z_, w_] := ( makeVector[z]; makeVector[w];
	With[ {n = Dimension[z + w]},
	4 ( n ( z.Delta[n]^2 + w.Delta[n]^2 ) - Norm[z]^2 - Norm[w]^2 +
	2 z.w + (2n - 4) z.Delta[n] w.Delta[n] ) / (
	n Volume[n] (Norm[z]^2 + Norm[w]^2 - 2 z.w +
	4 z.Delta[n] w.Delta[n])^(1 + n/2) )
    ] )

BergmanKernelH[ x_, y_, t_, u_ ] := ( makeVector[x] ; makeVector[t] ;
	unmakeVector[t]; unmakeVector[u] ;
	With[ {n = Dimension[x] + 1},
	4 ( (n-1)(y + u)^2 - Norm[x - t]^2 ) / (
	n Volume[n] (Norm[x - t]^2 + (y + u)^2)^(1 + n/2)) ] )
(*

Miscellaneous

HarmonicConjugate

*)
HarmonicConjugate[ u_, x_Symbol, y_Symbol ] := (
	unmakeVector[x] ; unmakeVector[y] ;
	Integrate[ D[ u, x ], {y, 0, y} ] -
	Integrate[ ( D[u, y] /. y -> 0 ), {x, 0, x} ] )
(*

BergmanProjection

*)
BergmanProjection[0, x_] := 0

BergmanProjection[u_, x_Symbol] :=
	( makeVector[x] ;
	Apply[ Plus, Map[ Apply[Plus, #]& , Map[ Map[ 
	Function[ j, ( (Dimension[x] + 2 #[[2]]) #[[1]] / 
	(#[[2]] + j + Dimension[x]) )& ][#[[2]]], #[[1]] ]&, 
	Map[ {polynomialDecomposition[ #[[1]] ,x ], #[[2]]}&, 
	Map[ {Dirichlet[ #[[1]], x], #[[2]]}& , 
	polynomialDecomposition[u, x] ] ] ] ] ] )

BergmanProjection[u_ , v_List] := explicit[ BergmanProjection, u, v ]
(*

Schwarz

*)
Schwarz[0] = 0 ;

Schwarz[x_] := ( If[ Head[x] == List, unmakeVector[x] ] ;
	Module[ {t, u},
	With[ {s = Integrate[
		Cos[u]^( (Dimension[x]-2) ) /
		(1 + 1/t^2  - 2 Sin[u]/t)^(Dimension[x]/2),
		{u, 0, Pi/2} ] },
	Together[
	( (1 - 1/t^2) (Dimension[x] - 1) Volume[Dimension[x]-1]/
	(Dimension[x] Volume[ Dimension[x] ])
	( s - (s /. t -> -t) ) )
		 /. t -> 1/Norm[x] ] ] ] /.
		 Literal[Plus[t_]] :>
		 Factor[ Select[Plus[t], FreeQ[#, ArcTan]& ] ] +
		 Factor[ Select[Plus[t], Not[ FreeQ[#, ArcTan] ] & ] ] /.
			(-1 + t_) (1 + t_) -> - (1 - t^2) ) /;
			EvenQ[ Dimension[x] ]

Schwarz[x_] := ( If[ Head[x] == List, unmakeVector[x] ] ;
	Apart[ ( Module[ {t, u},
	With[ {s = Integrate[
		Cos[u]^( (Dimension[x]-2) ) /
		(1 + t^2  - 2 Sin[u] t)^(Dimension[x]/2),
		{u, 0, Pi/2} ] },
	Together[
	( (1 - t^2) (Dimension[x] - 1) Volume[Dimension[x]-1]/
	(Dimension[x] Volume[ Dimension[x] ])
	( s - (s /. t -> -t) ) ) /. t -> Norm[x] /.
		 {((1 + Norm[x])^2)^Rational[k_,2] :> (1 + Norm[x])^k, 
	((-1 + Norm[x])^2)^Rational[k_,2] :> (1 - Norm[x])^k,
	((-1 - Norm[x])^2)^Rational[k_,2] :> (1 + Norm[x])^k } 
	] ] ] /.
	 Literal[Plus[t_]] :>
		 Factor[ Select[Plus[t], FreeQ[#, Power[_, 1/2]]& ] ] +
		 Factor[ Select[Plus[t], Not[ FreeQ[#, Power[_, 1/2]] ] & ] ] /.
		 (-1 + t_) (1 + t_) -> - (1 - t^2)   ) ] ) /;
			OddQ[ Dimension[x] ]
(*

IntegerQ

*)
IntegerQ[ m_?IntegerQ + n_?IntegerQ ] = True ;

IntegerQ[ m_?IntegerQ n_?IntegerQ ] = True ;
(*

Times

*)
Format[ Times[-1, x_?vector] ] := SequenceForm["-",x]

Format[ Literal[Times[a__]]   /; 
  Or @@ Thread[ vector[{a}] ] ]  :=
  With[ {t = Position[ Thread[ vector[{a}] ], True ][[1,1]]},
  SequenceForm[ If[ Head[Times @@ Delete[ {a}, t ]] === Plus, 
  SequenceForm[ "(", First[Delete[ {a}, t ]], ")"],
  Times @@ Delete[ {a}, t ] ], " ", 
  If[Head[{a}[[t]]] === Plus, SequenceForm["(",{a}[[t]],")"], {a}[[t]] ]
  ] ]

Format[ Times[-1, x_?matrix] ] := SequenceForm["-",x]

Format[ Literal[Times[a__]]   /; 
  Or @@ Thread[ matrix[{a}] ] ]  :=
  With[ {t = Position[ Thread[ matrix[{a}] ], True ][[1,1]]},
  SequenceForm[ If[ Head[Times @@ Delete[ {a}, t ]] === Plus, 
  SequenceForm[ "(", First[Delete[ {a}, t ]], ")"],
  Times @@ Delete[ {a}, t ] ], " ", 
  If[Head[{a}[[t]]] === Plus, SequenceForm["(",{a}[[t]],")"], {a}[[t]] ]
  ] ]
(*

explicit

*)
explicit[ function_, f_, v_ ] := ( unmakeVector[v];
	Module[{z}, SetDimension[z, Length[v] ] ;
	function[ f /. 
	Thread[v -> Array[z, Length[v]] ], z ]
	/. Prepend[Thread[Array[z, Length[v]] -> v ], 
	Norm[z] -> Norm[v] ] ] )

explicit[function_, f_, g_, v_List] := ( unmakeVector[v];
	Module[{z}, SetDimension[z, Length[v] ] ;
	function[
	f /. Thread[v -> Array[z, Length[v]] ], 
	g /. Thread[v -> Array[z, Length[v]] ],
	z ]
	/. Prepend[Thread[Array[z, Length[v]] -> v ], 
	Norm[z] -> Norm[v] ] ] )
(*

Closing

*)
Protect[ Release[protectedWords] ] ;

End[ ]

EndPackage[ ]

Print["***********************************"]

Print[" "]

Print["You can now use the functions in this package."]
(*

DELETE EVERYTHING ABOVE THE COMMENT LINE AT THE BEGINNING OF THIS FILE.
SAVE THE FILE IN TEXT ONLY (ASCII) FORM AFTER
DELETING EVERYTHING BELOW THIS LINE. ^*)