(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* ::Input::Initialization:: *)
If[$VersionNumber<7,Print["This package may not work correctly with versions of \!\(\*
StyleBox[\"Mathematica\",\nFontSlant->\"Italic\"]\) earlier than \!\(\*
StyleBox[\"Mathematica\",\nFontSlant->\"Italic\"]\) 7. Versions of this package that will work correctly with earlier versions of \!\(\*
StyleBox[\"Mathematica\",\nFontSlant->\"Italic\"]\) are available from Sheldon Axler, who can be contacted at axler@sfsu.edu."]]


(* ::Input::Initialization:: *)
Print["HFT12.m, version 12.03, 20 December 2020; for use with \!\(\*
StyleBox[\"Mathematica\",\nFontSlant->\"Italic\"]\) versions 7 through 12 (and probably later versions of \!\(\*
StyleBox[\"Mathematica\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\")\",\nFontSlant->\"Italic\"]\)."]


(* ::Input::Initialization:: *)
Print["The HFT12.m \!\(\*
StyleBox[\"Mathematica\",\nFontSlant->\"Italic\"]\) package is designed for computing with harmonic functions."]


(* ::Input::Initialization:: *)
Print["Documentation for the use of this package and information about the algorithms used in it is available in the document titled \!\(\*
StyleBox[\"Computing\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"with\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"Harmonic\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"Functions\",\nFontSlant->\"Italic\"]\), which is available in both nb and pdf formats."]


(* ::Input::Initialization:: *)
Print["The most recent version of this HFT12.m package and its documentation \!\(\*
StyleBox[\"Computing\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"with\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"Harmonic\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"Functions\",\nFontSlant->\"Italic\"]\) are available at http://www.axler.net."]


(* ::Input::Initialization:: *)
Print["For addional information about harmonic functions, see the book \!\(\*
StyleBox[\"Harmonic\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"Function\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"Theory\",\nFontSlant->\"Italic\"]\) (second edition), by Sheldon Axler, Paul Bourdon, and Wade Ramey, published by Springer."]


(* ::Input::Initialization:: *)
Print["This package is copyrighted by Sheldon Axler but is distributed without charge."]


(* ::Input::Initialization:: *)
Print["Comments, suggestions, and bug reports should be sent to axler@sfsu.edu."]


(* ::Input::Initialization:: *)
BeginPackage["HFT`"]


(* ::Input::Initialization:: *)
\[CapitalDelta]::"usage"="\!\(\*SubscriptBox[\(\[CapitalDelta]\), \(x\)]\)[f] gives the Laplacian of f with respect to x."


(* ::Input::Initialization:: *)
annulus::"usage"="annulus is an option for region in dirichlet."


(* ::Input::Initialization:: *)
antiLaplacian::"usage"="antiLaplacian[u, x] gives a function whose Laplacian equals u. Here u must be a polynomial function of \!\(\*FormBox[\(x\),
TraditionalForm]\) or a sum of terms, each of which is a polynomial in x times a function of \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]."


(* ::Input::Initialization:: *)
basisH::"usage"="basisH[m, x] gives a basis for the space of harmonic polynomials homogeneous of degree m in the variable x."


(* ::Input::Initialization:: *)
bergmanKernel::"usage"="bergmanKernel[x, y] gives the Bergman reproducing kernel for the unit ball."


(* ::Input::Initialization:: *)
bergmanKernelH::"usage"="bergmanKernelH[z, w] gives the Bergman reproducing kernel for the upper half-space."


(* ::Input::Initialization:: *)
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."


(* ::Input::Initialization:: *)
biDirichlet::"usage"="biDirichlet[p , x] solves the biDirichlet problem with boundary data p, as a function of x."


(* ::Input::Initialization:: *)
\[Delta]::"usage"="\!\(\*SubscriptBox[\(\[Delta]\), \(j\)]\) is the vector that equals 1 in the j-th coordinate and 0 in the other coordinates."


(* ::Input::Initialization:: *)
dimension::"usage"="dimension[x] is the dimension of x, which might be a vector or a matrix."


(* ::Input::Initialization:: *)
dimHarmonic::"usage"="dimHarmonic[m, n] gives the vector space dimension of the space of harmonic polynomials on \!\(\*SuperscriptBox[\(R\), \(n\)]\) homogeneous of degree m."


(* ::Input::Initialization:: *)
dirichlet::"usage"="dirichlet[p, x] gives the harmonic function that equals p on the unit sphere. Here p must be a polynomial function of x."


(* ::Input::Initialization:: *)
divergence::"usage"="divergence[f, x] gives the divergence of f with respect to x."


(* ::Input::Initialization:: *)
expandNorm::"usage"="expandNorm[f] gives f with all terms of the form \[LeftDoubleBracketingBar]x + y\[RightDoubleBracketingBar] replaced by \!\(\*SqrtBox[\(\*SuperscriptBox[\(\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\), \(2\)]\\\  + \\\ \*SuperscriptBox[\(\[LeftDoubleBracketingBar]y\[RightDoubleBracketingBar]\), \(2\)]\\\  + \\\ 2\\\ x . y\)]\)."


(* ::Input::Initialization:: *)
exteriorNeumann::"usage"="exteriorNeumann[p , x] solves the exterior Neumann problem with boundary data p, as a function of x."


(* ::Input::Initialization:: *)
exteriorSphere::"usage"="exteriorSphere is an option for dirichlet, specifying that the region should be the exterior of the unit sphere."


(* ::Input::Initialization:: *)
gradient::"usage"="gradient[f, x] gives the gradient of f with respect to x."


(* ::Input::Initialization:: *)
harmonicConjugate::"usage"="harmonicConjugate[u, x, y] gives the harmonic conjugate of u with respect to the variables x, y."


(* ::Input::Initialization:: *)
harmonicDecomposition::"usage"="harmonicDecomposition[p, x] gives the harmonic decomposition of p with respect to the variable x."


(* ::Input::Initialization:: *)
hilbertSchmidt::"usage"="hilbertSchmidt[A] gives the Hilbert-Schmidt norm of a matrix A."


(* ::Input::Initialization:: *)
homogeneous::"usage"="homogeneous[f, d, x] gives the homogenous part of f of degree d, with respect to x."


(* ::Input::Initialization:: *)
If[ $VersionNumber< 10.2, Hyperplane::"usage"="Hyperplane[n, c] represents the hyperplane with normal n given by the points x that satisfy n.x = c."]


(* ::Input::Initialization:: *)
integrateBall::"usage"="integrateBall[f, x] gives the integral of f, as a function of x, over the unit ball."


(* ::Input::Initialization:: *)
integrateEllipsoidArea::"usage"="integrateEllipsoidArea[f, b, c, d, x] gives the integral of f/\[LeftDoubleBracketingBar]\[Del]q\[RightDoubleBracketingBar] as a function of x with repsect to surface area measure over the ellipsoid {x: q(x) = 0}, where q(x) = b.\!\(\*SuperscriptBox[\(x\), \(2\)]\)+ c.x + d. Here b and c must be lists whose length equals dimension[x] and f must be a polynomial function of x. If c is the 0 vector, then c can be omitted. If d = -1, then d can be omitted."


(* ::Input::Initialization:: *)
integrateEllipsoidVolume::"usage"="integrateEllipsoidVolume[f, b, c, d, x] gives the integral of f as a function of x with repsect to volume measure over the region {x: q(x) < 0}, where q(x) = b.\!\(\*SuperscriptBox[\(x\), \(2\)]\)+ c.x + d. Here b and c must be lists whose length equals dimension[x] and f must be a polynomial function of x. If c is the 0 vector, then c can be omitted. If d = -1, then d can be omitted."


(* ::Input::Initialization:: *)
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."


(* ::Input::Initialization:: *)
jacobian::"usage"="jacobian[f, x] gives the Jacobian of f with respect to x."


(* ::Input::Initialization:: *)
kelvin::"usage"="kelvin[u, x] gives the Kelvin transform of u, thought of as a function of x."


(* ::Input::Initialization:: *)
kelvinH::"usage"="kelvinH is the modified Kelvin transform, as defined in Chapter 7 of Harmonic Function Theory."


(* ::Input::Initialization:: *)
laplacian::"usage"="laplacian equals \[CapitalDelta], the Laplacian."


(* ::Input::Initialization:: *)
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 \[LeftDoubleBracketingBar]x\!\(\*SuperscriptBox[\(\[RightDoubleBracketingBar]\), \(2\)]\), where x is the variable."


(* ::Input::Initialization:: *)
neumann::"usage"="neumann[p, x] solves the Neumann problem of finding a harmonic function of x whose outward normal derivative (on the unit sphere) is p."


(* ::Input::Initialization:: *)
norm::"usage"="norm[x] gives the norm of a vector x."


(* ::Input::Initialization:: *)
normalD::"usage"="normalD[f, z] gives the outward normal derivative of f, as a function of z, with respect to the unit ball."


(* ::Input::Initialization:: *)
partial::"usage" = "partial[f, \!\(\*SubscriptBox[\(x\), \(j\)]\)] gives the partial derivative of f with respect to \!\(\*SubscriptBox[\(x\), \(j\)]\)."


(* ::Input::Initialization:: *)
\[CapitalPhi]::"usage"="\[CapitalPhi][z] is the modified reflection defined in Chapter 7 of Harmonic Function Theory."


(* ::Input::Initialization:: *)
poissonKernel::"usage"="poissonKernel[x, z] gives the Poisson kernel for the unit ball."


(* ::Input::Initialization:: *)
poissonKernelH::"usage"="poissonKernelH[x, y, t] gives the Poisson kernel for the upper half-space."


(* ::Input::Initialization:: *)
quadratic::"usage"="quadratic is an option for multiple in antiLaplacian and for region in dirichlet."


(* ::Input::Initialization:: *)
reflection::"usage"="reflection[x] gives the reflection of x in the unit sphere."


(* ::Input::Initialization:: *)
region::"usage"="region is an option for dirichlet."


(* ::Input::Initialization:: *)
southPole::"usage"="southPole is the south pole."


(* ::Input::Initialization:: *)
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."


(* ::Input::Initialization:: *)
setDimension::"usage"="setDimension[x, n] sets the dimension of x to n."


(* ::Input::Initialization:: *)
singularity::"usage"="singularity is an option for antiLaplacian. The default value is None."


(* ::Input::Initialization:: *)
surfaceArea::"usage"="surfaceArea[n] gives the surface area of the unit ball in n-dimensional real Euclidean space."


(* ::Input::Initialization:: *)
taylor::"usage"="taylor[f, d, x] gives the Taylor expansion of f to degree d, with respect to x."


(* ::Input::Initialization:: *)
togetherness::"usage"="turnOff[ togetherness ] turns off the feature that simplifies some output. turnOn[ togetherness ] turns this feature back on."


(* ::Input::Initialization:: *)
trace::"usage"="trace[A] gives the trace of the matrix A."


(* ::Input::Initialization:: *)
turnOff::"usage"="turnOff[F] turns feature F off."


(* ::Input::Initialization:: *)
turnOn::"usage"="turnOn[F] turns feature F on."


(* ::Input::Initialization:: *)
vectorLast::"usage" =
"turnOff[ vectorLast] turns off the feature that displays the vector or matrix last in the product of scalar or scalars and a vector or matrix. \[IndentingNewLine]turnOn[ vectorLast ] turns this feature back on."


(* ::Input::Initialization:: *)
volume::"usage"="volume[n] gives the volume of the unit ball in n-dimensional real Euclidean space."


(* ::Input::Initialization:: *)
zeroToZero::"usage"="turnOn[ zeroToZero ] turns on the feature that changes 0. to 0 via $Post.  turnOff[ zeroToZero ] turns this feature back off."


(* ::Input::Initialization:: *)
zonalHarmonic::"usage"="zonalHarmonic[m, x, z] gives the zonal harmonic of degree m with pole x."


(* ::Input::Initialization:: *)
Begin["`private`"]


(* ::Input::Initialization:: *)
Unprotect[IntegerQ,Limit,Transpose];


(* ::Input::Initialization:: *)
protectedWords=Unprotect[Abs,Dot,IdentityMatrix, Integer, Power, Times]


(* ::Input::Initialization:: *)
SetAttributes[postTogether,HoldAll]


(* ::Input::Initialization:: *)
postTogether[Apart[f_]]:=Apart[f]


(* ::Input::Initialization:: *)
postTogether[Together[f_]]:=Together[f]


(* ::Input::Initialization:: *)
postTogether[Simplify[f_]]:=Simplify[f]


(* ::Input::Initialization:: *)
postTogether[Expand[f_]]:=Expand[f]


(* ::Input::Initialization:: *)
postTogether[a_=b_]:=a=postTogether[b]


(* ::Input::Initialization:: *)
shorter[f_] := With[ {g = Together[f]},
	If[ StringLength[ ToString[ Format[ g, InputForm ] ] ] < 
	StringLength[ ToString[ Format[ f, InputForm ] ] ], g, f ] ]


(* ::Input::Initialization:: *)
shorter/:f_Symbol[a___,HoldPattern[shorter][x_],b___]:=f[a,x,b]/;MemberQ[Attributes[f],HoldAll]||MemberQ[Attributes[f],HoldFirst]||MemberQ[Attributes[f],HoldRest]


(* ::Input::Initialization:: *)
postTogether[f_]:=Map[shorter,f,{0,-2}]//.HoldPattern[shorter][j_]->j


(* ::Input::Initialization:: *)
zeroToZero[f_] := (f /. 0. -> 0);


(* ::Input::Initialization:: *)
usersPost=If[ToString[$Post]==="$Post",Identity,$Post];


(* ::Input::Initialization:: *)
togetherness::Off="togetherness has been turned off.  "<>"The command  turnOn[togetherness]  will turn it back on."


(* ::Input::Initialization:: *)
togetherness::On="togetherness has been turned on.  "<>"The command  turnOff[togetherness]  will turn it back off."


(* ::Input::Initialization:: *)
togetherness::AlreadyOn="togetherness is already turned on."


(* ::Input::Initialization:: *)
togetherness::AlreadyOff="togetherness is already turned off."


(* ::Input::Initialization:: *)
zeroToZero::Off="zeroToZero has been turned off.  "<>"The command  turnOn[zeroToZero]  will turn it back on."


(* ::Input::Initialization:: *)
zeroToZero::On="zeroToZero has been turned on.  "<>"The command  turnOff[zeroToZero]  will turn it back off."


(* ::Input::Initialization:: *)
zeroToZero::AlreadyOn="zeroToZero is already turned on."


(* ::Input::Initialization:: *)
zeroToZero::AlreadyOff="zeroToZero is already turned off."


(* ::Input::Initialization:: *)
vectorLast::Off="vectorLast has been turned off.  "<>"The command  turnOn[vectorLast]  will turn it back on."


(* ::Input::Initialization:: *)
vectorLast::On="vectorLast has been turned on.  "<>"The command  turnOff[vectorLast]  will turn it back off."


(* ::Input::Initialization:: *)
vectorLast::AlreadyOn="vectorLast is already turned on."


(* ::Input::Initialization:: *)
vectorLast::AlreadyOff="vectorLast is already turned off."


(* ::Input::Initialization:: *)
turnOn[togetherness]^:= If[ togetherness2 , Message[togetherness::AlreadyOn],$Post=If[zerotozero,Composition[zeroToZero, postTogether, usersPost],Composition[postTogether,usersPost]];togetherness2=True; Message[togetherness::On]]


(* ::Input::Initialization:: *)
turnOn[zeroToZero]^:= If[ zerotozero , Message[zeroToZero::AlreadyOn],$Post=If[togetherness2,Composition[zeroToZero, postTogether, usersPost],Composition[zeroToZero,usersPost]];zerotozero=True; Message[zeroToZero::On]]


(* ::Input::Initialization:: *)
turnOn[vectorLast] ^:= If[vectorMatricesLast, Message[vectorLast::AlreadyOn],protectTimes = Unprotect[Times];
timesAdjust1; timesAdjust2;
Protect[Evaluate[protectTimes]]; vectorMatricesLast = True; Message[vectorLast::On]]


(* ::Input::Initialization:: *)
$Post=Composition[postTogether, usersPost];


(* ::Input::Initialization:: *)
togetherness2=True; zerotozero = False; vectorMatricesLast = True;


(* ::Input::Initialization:: *)
turnOff[zeroToZero]^:= If[ zerotozero ,$Post=If[togetherness2,Composition[postTogether,usersPost], usersPost];zerotozero=False;Message[zeroToZero::Off], 
Message[zeroToZero::AlreadyOff]]


(* ::Input::Initialization:: *)
turnOff[togetherness]^:= If[ togetherness2 ,$Post=If[zerotozero,Composition[zeroToZero,usersPost], usersPost];togetherness2=False;Message[togetherness::Off], 
Message[togetherness::AlreadyOff]]


(* ::Input::Initialization:: *)
turnOff[vectorLast] ^:= If[vectorMatricesLast,
 \!\(\*
TagBox[
RowBox[{
FormBox[
RowBox[{
TagBox[
RowBox[{
FormBox[
RowBox[{
RowBox[{
RowBox[{"protectTimes", " ", "=", " ", 
RowBox[{"Unprotect", "[", "Times", "]"}]}], ";"}], " ", "\[IndentingNewLine]", 
RowBox[{
RowBox[{
TagBox[
FormBox[
RowBox[{"HoldPattern", "[", 
RowBox[{"Times", "[", "a__", "]"}], "]"}],
StandardForm],
Format], " ", "=", " ", "."}], ";"}], " "}],
StandardForm], " "}],
Format], " "}],
StandardForm], "\[IndentingNewLine]", " ", 
RowBox[{
RowBox[{
TagBox[
FormBox[
RowBox[{"HoldPattern", "[", 
RowBox[{"Times", "[", 
RowBox[{
RowBox[{"-", "1"}], ",", "a__"}], "]"}], " ", "]"}],
StandardForm],
Format], "=", " ", "."}], ";"}], "\[IndentingNewLine]"}],
Format]\)Protect[Evaluate[protectTimes]];vectorMatricesLast = False; Message[vectorLast::Off],Message[vectorLast::AlreadyOff]]


(* ::Input::Initialization:: *)
matrix[IdentityMatrix[_]]^=True;


(* ::Input::Initialization:: *)
IdentityMatrix/:a_ . IdentityMatrix[_]:=a


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


(* ::Input::Initialization:: *)
Off[Transpose::"nmtx"]


(* ::Input::Initialization:: *)
matrix[Transpose[_]]=True


(* ::Input::Initialization:: *)
Transpose[Transpose[x_]]:=x


(* ::Input::Initialization:: *)
Transpose[{c_ x_?vector}]:=c Transpose[{x}]


(* ::Input::Initialization:: *)
Transpose[c_ x_?matrix]:=c Transpose[x]


(* ::Input::Initialization:: *)
Transpose[IdentityMatrix[n_]]:=IdentityMatrix[n]


(* ::Input::Initialization:: *)
Transpose[{0}]:={0}


(* ::Input::Initialization:: *)
trace[IdentityMatrix[n_]]:=n


(* ::Input::Initialization:: *)
trace[Transpose[A_]]:=trace[A]


(* ::Input::Initialization:: *)
trace[c_ A_]:=c trace[A]/;NumericQ[c]


(* ::Input::Initialization:: *)
trace[c_ A_]:=c trace[A]/;matrix[A]


(* ::Input::Initialization:: *)
trace[A_]:=Tr[A] /; MatrixQ[A]


(* ::Input::Initialization:: *)
Format[norm[x_], StandardForm] := DoubleBracketingBar[x]


(* ::Input::Initialization:: *)
MakeExpression[RowBox[{"\[LeftDoubleBracketingBar]", x_, "\[RightDoubleBracketingBar]"}],StandardForm]  := MakeExpression[RowBox[{"norm", "[", x, "]"}], StandardForm]


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


(* ::Input::Initialization:: *)
norm[0] :=0


(* ::Input::Initialization:: *)
norm[x_List]:= Sqrt[Plus @@ (x^2)]


(* ::Input::Initialization:: *)
norm[c_ x_?vector]:=(realize[c];Abs[c] norm[x])


(* ::Input::Initialization:: *)
norm[x_]:=0/;(makeVector[x];False)


(* ::Input::Initialization:: *)
norm[c_?NumericQ x_]:=Abs[c] norm[x]


(* ::Input::Initialization:: *)
real = {};


(* ::Input::Initialization:: *)
realize[c_] := (real = real\[Union]{c})


(* ::Input::Initialization:: *)
expandNorm[f_]:=f//.norm[g_+h_]:>Sqrt[norm[g]^2+norm[h]^2+2 g . h]


(* ::Input::Initialization:: *)
Abs[norm[x_]^p_.]:=norm[x]^p


(* ::Input::Initialization:: *)
Abs/:Abs[Subscript[x_, j_]]^p_?EvenQ:=Subscript[x, j]^p/;vector[x]


(* ::Input::Initialization:: *)
Abs/:Abs[c_]^p_?EvenQ:=c^p/;MemberQ[real, c] ||MemberQ[real,1/c]||MemberQ[real,c^(p/2)]||MemberQ[real,1/c^(p/2)]==0


(* ::Input::Initialization:: *)
vector[HoldPattern[+a__]]:=If[TrueQ[Or@@Thread[vector[{a}]]],makeVector[+a];True]


(* ::Input::Initialization:: *)
vector[HoldPattern[Times[a__]]]:=Or@@Thread[vector[{a}]]


(* ::Input::Initialization:: *)
vector[A_?matrix . x_?vector]:=True


(* ::Input::Initialization:: *)
vector[x_?vector . Y_?matrix]:=True


(* ::Input::Initialization:: *)
vector[a_[".",_]]:=True


(* ::Input::Initialization:: *)
vector[a_[_,"."]]:=True


(* ::Input::Initialization:: *)
Subscript[0, _]=0


(* ::Input::Initialization:: *)
vector[Subscript[partial, _][b_][x_]]:=True/;vector[b[x]]


(* ::Input::Initialization:: *)
vector[gradient[_][_]]=True


(* ::Input::Initialization:: *)
vectors={};


(* ::Input::Initialization:: *)
vectorValuedFunctions={};


(* ::Input::Initialization:: *)
makeVector[x_Symbol]:=If[MemberQ[vectors,x],Null,x[j_]=Subscript[x, j];AppendTo[vectors,x];vector[x]^=True;]


(* ::Input::Initialization:: *)
makeVector[c_?NumericQ x_]:=makeVector[x]


(* ::Input::Initialization:: *)
makeVector[HoldPattern[+a__]]:=makeVector/@{a}


(* ::Input::Initialization:: *)
makeVector[f_Symbol[_]]:=If[MemberQ[vectorValuedFunctions,f],Null,f[x_][j_]=Subscript[f, j][x];
		AppendTo[vectorValuedFunctions,f];vector[f[_]]^=True;vector[(f^(_))[_]]=True;vector[Subscript[partial, _][f][_]]=True;]/;f=!=List


(* ::Input::Initialization:: *)
makeVector[Subscript[partial, _][u_][_]]:=makeVector[u[_]]


(* ::Input::Initialization:: *)
Attributes[unmakeVector]={Listable};


(* ::Input::Initialization:: *)
unmakeVector[x_Symbol]:= (Unprotect[Abs];Abs/:Abs[x]^p_?EvenQ:=x^p;If[MemberQ[vectors,x],(x[j$_])=.;vectors=Delete[vectors,Position[vectors,x][[1,1]]];vector[x]^=False;])


(* ::Input::Initialization:: *)
matrix[HoldPattern[Times[a__]]]:=Or@@Thread[matrix[{a}]]


(* ::Input::Initialization:: *)
matrix[HoldPattern[+a__]]:=Or@@Thread[matrix[{a}]]


(* ::Input::Initialization:: *)
matrix[a_?matrix . x_?matrix]:=True


(* ::Input::Initialization:: *)
matrix[jacobian[_][_]]=True;


(* ::Input::Initialization:: *)
matrix[{_}]=True;


(* ::Input::Initialization:: *)
Subscript[0, _,_]=0


(* ::Input::Initialization:: *)
matrix[Transpose[x_]]:=matrix[x]


(* ::Input::Initialization:: *)
matrices={};


(* ::Input::Initialization:: *)
makeMatrix[x_Symbol]:=If[TrueQ[matrix[x]],Null,x[j_,k_]=Subscript[x, j,k];AppendTo[matrices,x];matrix[x]^=True;vector[x]^=False]


(* ::Input::Initialization:: *)
Attributes[Dot]={OneIdentity}


(* ::Input::Initialization:: *)
b_ . 0 = 0;


(* ::Input::Initialization:: *)
0 . b_ = 0;


(* ::Input::Initialization:: *)
HoldPattern[+x__] . z_:=Plus@@(#1 . z&)/@{x}


(* ::Input::Initialization:: *)
a_ . HoldPattern[+x__]:=Plus@@(a . #1&)/@{x}


(* ::Input::Initialization:: *)
(c_ x_) . y_:=c x . y/;vector[x]||matrix[x]||NumericQ[c]


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


(* ::Input::Initialization:: *)
a_ . (c_ y_):=c a . y/;vector[y]||matrix[y]||NumericQ[c]


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


(* ::Input::Initialization:: *)
x_ . x_:=norm[x]^2/;vector[x]


(* ::Input::Initialization:: *)
f_[x_] . f_[x_] := norm[f[x]]^2 /; Head[f]==Symbol &&vector[f[x]]


(* ::Input::Initialization:: *)
Derivative[s_][f_][x_] . Derivative[s_][f_][x_] := norm[Derivative[s][f][x] ]^2/; Head[f] == Symbol &&vector[f[x]]


(* ::Input::Initialization:: *)
Subscript[partial, j_][f_][x_] . Subscript[partial, j_][f_][x_] := norm[ Subscript[partial, j][f][x]]^2/; Head[f] == Symbol &&vector[f[x]]


(* ::Input::Initialization:: *)
Subscript[\[Delta], j_] . Subscript[\[Delta], k_]:=Subscript[\[Delta], j,k]


(* ::Input::Initialization:: *)
Subscript[\[Delta], j_] . x_Symbol:=(makeVector[x];Subscript[x, j])/;matrix[x]=!=True


(* ::Input::Initialization:: *)
a_Symbol . Subscript[\[Delta], j_]:=(makeVector[a];Subscript[a, j])/;matrix[a]=!=True


(* ::Input::Initialization:: *)
Subscript[\[Delta], j_] . f_Symbol[x_]:=(makeVector[f[_]];f[x][j])/;matrix[f[x]]=!=True/;f=!=List


(* ::Input::Initialization:: *)
b_Symbol[x_] . Subscript[\[Delta], j_]:=(makeVector[b[_]];b[x][j])/;matrix[b[x]]=!=True/;b=!=List


(* ::Input::Initialization:: *)
Subscript[\[Delta], j_] . gradient[u_][x_]:=\!\(\*SubscriptBox[\(partial\), \({j}\)]\)[u][x]


(* ::Input::Initialization:: *)
gradient[u_][x_] . Subscript[\[Delta], j_]:=\!\(\*SubscriptBox[\(partial\), \({j}\)]\)[u][x]


(* ::Input::Initialization:: *)
jacobian[b_][x_] . Subscript[\[Delta], j_]:=\!\(\*SubscriptBox[\(partial\), \({j}\)]\)[b][x]


(* ::Input::Initialization:: *)
Subscript[\[Delta], j_] . jacobian[b_][x_]:=gradient[b . Subscript[\[Delta], j]][x]


(* ::Input::Initialization:: *)
Subscript[\[Delta], j_] . Y_?matrix:=Subscript[Y, j,"."]


(* ::Input::Initialization:: *)
A_?matrix . Subscript[\[Delta], j_]:=\!\(\*SubscriptBox[\(A\), \("\<.\>", j\)]\)


(* ::Input::Initialization:: *)
Subscript[\[Delta], k_] . Subscript[Y, j_,"."]:=Subscript[Y, j,k]


(* ::Input::Initialization:: *)
Subscript[A, j_,"."] . Subscript[\[Delta], k_]:=Subscript[A, j,k]


(* ::Input::Initialization:: *)
Subscript[\[Delta], j_] . \!\(\*SubscriptBox[\(Y\), \("\<.\>", k_\)]\):=Subscript[Y, j,k]


(* ::Input::Initialization:: *)
\!\(\*SubscriptBox[\(A\), \("\<.\>", k_\)]\) . Subscript[\[Delta], j_]:=Subscript[A, j,k]


(* ::Input::Initialization:: *)
a_?VectorQ . Subscript[\[Delta], j_]:=a[[j]]


(* ::Input::Initialization:: *)
Subscript[\[Delta], j_] . v_?VectorQ:=v[[j]]


(* ::Input::Initialization:: *)
A_?MatrixQ . Subscript[\[Delta], j_]:=Transpose[A][[j]]


(* ::Input::Initialization:: *)
Subscript[\[Delta], j_] . Y_?MatrixQ:=A[[j]]


(* ::Input::Initialization:: *)
x_ . A_?identityMatrixMultiple:=A[[1,1]] x/;primitiveVector[x]


(* ::Input::Initialization:: *)
A_?identityMatrixMultiple . x_:=A[[1,1]] x/;primitiveVector[x]


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


(* ::Input::Initialization:: *)
v_ . Transpose[{x_Symbol}] . {x_Symbol} := (v . x)x


(* ::Input::Initialization:: *)
a_List . x_Symbol:=(makeVector[x];a . Table[Subscript[x, j],{j,Length[a]}])/;matrix[x]=!=True


(* ::Input::Initialization:: *)
x_Symbol . y_List:=(makeVector[x];Table[Subscript[x, j],{j,Length[y]}] . y)/;matrix[x]=!=True


(* ::Input::Initialization:: *)
a_List . B_Symbol:=a . Table[ Subscript[B, j,k], {j, dimension[B][[1]]}, {k,  dimension[B][[2]]}]/;matrix[B]


(* ::Input::Initialization:: *)
B_Symbol . c_List:=Table[ Subscript[B, j,k], {j, dimension[B][[1]]}, {k,  dimension[B][[2]]}] . c/;matrix[B]


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


(* ::Input::Initialization:: *)
primitiveVector[x_Symbol]:=True/;Context[x]=!="System`"


(* ::Input::Initialization:: *)
primitiveVector[gradient[_][_]]=True


(* ::Input::Initialization:: *)
primitiveVector[\[Delta][_]]=True


(* ::Input::Initialization:: *)
degree[f_,z_]:=degree1[Expand[f/.norm[z]->Subscript[z, normReplacement]],z]


(* ::Input::Initialization:: *)
degree[f_,z_]:=degree1[Expand[f/.{norm[z]->Subscript[z, normReplacement], z . _ ->Subscript[z, normReplacement2], _ . z -> Subscript[z, normReplacement3]}] ,z]


(* ::Input::Initialization:: *)
degree1[0,_]:=-\[Infinity]


(* ::Input::Initialization:: *)
degree1[HoldPattern[+a__],z_]:=Max[Thread[degree1[{a},z]]]


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


(* ::Input::Initialization:: *)
uninomialQ[f_, x_Symbol] := MatchQ[f,Subscript[x, _]] || MatchQ[f, Subscript[x, _]^_?((IntegerQ[#1] &&Positive[#1])&)]


(* ::Input::Initialization:: *)
monomialQ[f_, x_Symbol] := uninomialQ[f,x] || (Head[f] === Times && (And@@(uninomialQ[f[[#1]],x]&)/@Range[Length[f]]))


(* ::Input::Initialization:: *)
polynomialDecomposition[0,x_Symbol]:=Null


(* ::Input::Initialization:: *)
polynomialDecomposition[p_,x_Symbol]:=With[{p1=Expand[p/.norm[x]->Subscript[x, normReplacement]]},With[{t=({#1,degree[#1,x]}&)/@If[Head[p1]===Plus,List@@p1,{p1}]},With[{s=Table[If[MemberQ[t,{_,j}],{Plus@@Transpose[Cases[t,{_,j}]][[1]],j},Null],{j,0,Max[Transpose[t][[2]]]}]},Delete[s,Position[s,Null]]]]]/.Subscript[x, normReplacement]->norm[x]


(* ::Input::Initialization:: *)
exponentsMod2[t_,x_] :=Drop[Mod[Cases[t  (Times @@ Array[x, dimension[x]+1])^2, Power[Subscript[x,_],_]][[All,2]],2], -1]


(* ::Input::Initialization:: *)
block[p_, x_, ex_] := Select[p, (exponentsMod2[#,x] == ex) &]


(* ::Input::Initialization:: *)
blocks[{p_, m_}, x_] := {{p, m, ConstantArray[0, dimension[x]]}} /; m == 0


(* ::Input::Initialization:: *)
blocks[{p_, m_}, x_] := {{p, m, exponentsMod2[p,x]}} /;Not[ Head[p] === Plus]


(* ::Input::Initialization:: *)
blocks[{p_, m_}, x_] := Module[{q=p, b = {}, ex, bq},While[Head[q]===Plus,ex = exponentsMod2[q[[1]], x]; bq = block[q, x, ex]; AppendTo[b, {bq,m,ex}]; q = q - bq];If[q===0,b, Append[b, {q, m, exponentsMod2[q,x]}]]]


(* ::Input::Initialization:: *)
blockPolynomialDecomposition[p_, x_] :=
 Flatten[blocks[#, x]& /@ polynomialDecomposition[p,x],1]


(* ::Input::Initialization:: *)
vector[Subscript[\[Delta], _]]^:=True


(* ::Input::Initialization:: *)
Subscript[\[Delta], j_,k_]:=0/;j!=k


(* ::Input::Initialization:: *)
Subscript[\[Delta], j_,k_]:=1/;j==k


(* ::Input::Initialization:: *)
Subscript[\[Delta], j_][k_]:=Subscript[\[Delta], j,k]


(* ::Input::Initialization:: *)
norm[Subscript[\[Delta], _]]^=1;


(* ::Input::Initialization:: *)
partial[f_,{Subscript[x_, _],0}]:=(makeVector[x];f)


(* ::Input::Initialization:: *)
partial[f_,{y_Symbol,0}]:=(unmakeVector[y];f)


(* ::Input::Initialization:: *)
partial[f_,{z_,n_Integer}]:=Module[{g = f}, Do[g = Together[partial[g,z]],{n}];g]/;Positive[n]


(* ::Input::Initialization:: *)
partial[f_,z1_,z__]:=Module[{g =partial[f,z1]}, Do[g = partial[g, {z}[[j]]], {j, Length[{z}]}]; g]


(* ::Input::Initialization:: *)
partial[f_,Subscript[x_, _]]:=(makeVector[x];0)/;FreeQ[f,x]


(* ::Input::Initialization:: *)
partial[f_,y_Symbol]:=(unmakeVector[y];0/;FreeQ[f,y])


(* ::Input::Initialization:: *)
partial[Subscript[x_, j_],Subscript[x_, k_]]:=(makeVector[x];Subscript[\[Delta], j,k])


(* ::Input::Initialization:: *)
partial[y_,y_]:=1


(* ::Input::Initialization:: *)
partial[x_,Subscript[x_, j_]]:=(makeVector[x];Subscript[\[Delta], j])


(* ::Input::Initialization:: *)
partial[HoldPattern[+a__],y_]:=Plus@@Thread[partial[{a},y]]


(* ::Input::Initialization:: *)
partial[f_ g_,y_]:=f partial[g,y]+partial[f,y] g


(* ::Input::Initialization:: *)
partial[f_ . g_,y_]:=(If[TrueQ[matrix[f]],Null,makeVector[f]];If[TrueQ[matrix[g]],Null,makeVector[g]];f . partial[g,y]+partial[f,y] . g)


(* ::Input::Initialization:: *)
partial[f_[v_List],t_]:=Plus@@(\!\(\*SubscriptBox[\(partial\), \({#1}\)]\)[f][v] * partial[v[[#1]],t]&)/@Range[Length[v]]


(* ::Input::Initialization:: *)
partial[f_^p_,y_]:=p f^(p-1) partial[f,y]+f^p Log[f] partial[p,y]


(* ::Input::Initialization:: *)
partial[f_[g_],Subscript[x_, j_]]:=(makeVector[x];If[TrueQ[vector[g]],If[TrueQ[vector[f[g]]],jacobian[f][g] . partial[g,Subscript[x, j]],gradient[f][g] . partial[g,Subscript[x, j]]],Derivative[1][f][g] partial[g,Subscript[x, j]]])


(* ::Input::Initialization:: *)
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]])


(* ::Input::Initialization:: *)
Subscript[partial, j_List][Subscript[partial, k_List][f_]]:=Subscript[partial, Sort[Join[k,j]]][f]


(* ::Input::Initialization:: *)
normalD[f_,z_Symbol]:=gradient[f,z] . z/.norm[z]->1


(* ::Input::Initialization:: *)
normalD[f_,v_List]:=(unmakeVector[v];Together[gradient[f,v] . v]/.norm[v]^2->1)


(* ::Input::Initialization:: *)
normalD[f_, q_, z_Symbol] := gradient[f, z] . gradient[q, z]/Sqrt[gradient[q, z] . gradient[q, z]]


(* ::Input::Initialization:: *)
normalD[f_,q_,v_List]:=(unmakeVector[v];Together[gradient[f, v] . gradient[q, v]/Sqrt[gradient[q, v] . gradient[q, v]]])


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


(* ::Input::Initialization:: *)
gradient[Subscript[x_, j_],x_Symbol]:=(makeVector[x];Subscript[\[Delta], j])


(* ::Input::Initialization:: *)
gradient[HoldPattern[+a__],x_Symbol]:=Plus@@Thread[gradient[{a},x]]


(* ::Input::Initialization:: *)
gradient[f_ g_,x_Symbol]:=f gradient[g,x]+gradient[f,x] g


(* ::Input::Initialization:: *)
gradient[f_ . g_,x_Symbol]:=f . jacobian[g,x]+g . jacobian[f,x]


(* ::Input::Initialization:: *)
gradient[f_^p_,x_Symbol]:=p f^(p-1) gradient[f,x]+f^p Log[f] gradient[p,x]


(* ::Input::Initialization:: *)
gradient[f_[y_],x_Symbol]:=(makeVector[x];If[TrueQ[vector[y]],gradient[f][y] . jacobian[y,x],Derivative[1][f][y] gradient[y,x]])


(* ::Input::Initialization:: *)
gradient[ Gamma[t_, f_], x_Symbol]:= -E^-f f^(t-1) gradient[f, x] /; FreeQ[t, x]


(* ::Input::Initialization:: *)
gradient[norm]:=#1/norm[#1]&


(* ::Input::Initialization:: *)
gradient[dimension]:=0&


(* ::Input::Initialization:: *)
gradient[f_,v_List]:=(unmakeVector[v];(partial[f,#1]&)/@v)


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


(* ::Input::Initialization:: *)
divergence[x_,x_Symbol]:=(makeVector[x];dimension[x])


(* ::Input::Initialization:: *)
divergence[HoldPattern[+a__],x_Symbol]:=Plus@@Thread[divergence[{a},x]]


(* ::Input::Initialization:: *)
divergence[f_ g_,x_Symbol]:=f divergence[g,x]+gradient[f,x] . g/;(makeVector[x];vector[g])


(* ::Input::Initialization:: *)
divergence[A_?matrix . f_,x_Symbol]:=trace[A . jacobian[f,x]]/;FreeQ[A,x]


(* ::Input::Initialization:: *)
divergence[f_ . T_?matrix,x_Symbol]:=trace[Transpose[T] . jacobian[f,x]]/;FreeQ[T,x]


(* ::Input::Initialization:: *)
divergence[u_List,v_List]:=(unmakeVector[v];Plus@@Thread[partial2[u,v]]/.partial2->partial)


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


(* ::Input::Initialization:: *)
jacobian[x_,x_Symbol]:=(makeVector[x];IdentityMatrix[dimension[x]])


(* ::Input::Initialization:: *)
jacobian[HoldPattern[+a__],x_Symbol]:=Plus@@Thread[jacobian[{a},x]]


(* ::Input::Initialization:: *)
jacobian[f_ g_?vector,x_Symbol]:=f jacobian[g,x]+Transpose[{gradient[f,x]}] . {g}


(* ::Input::Initialization:: *)
jacobian[f_Symbol[x_],x_Symbol]:=(makeVector[x];makeVector[f[_]];jacobian[f][x])/;f=!=List


(* ::Input::Initialization:: *)
jacobian[f_[y_],x_Symbol]:=If[TrueQ[vector[y]],jacobian[f][y] . jacobian[y,x],Transpose[{Derivative[1][f][y]}] . {gradient[y,x]}]


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


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


(* ::Input::Initialization:: *)
dimension[jacobian[f_][x_]]:= {dimension[f[_]],dimension[x] }


(* ::Input::Initialization:: *)
jacobian[f_List,v_List]:=(unmakeVector[v];Table[partial[f[[m]],v[[n]]],{m,Length[f]},{n,Length[v]}])


(* ::Input::Initialization:: *)
hilbertSchmidt[A_?MatrixQ]:=Plus@@Plus@@(A^2)


(* ::Input::Initialization:: *)
hilbertSchmidt[IdentityMatrix[n_]]:=n


(* ::Input::Initialization:: *)
hilbertSchmidt[ c_?NumericQ A_] := Abs[c] hilbertSchmidt[A]


(* ::Input::Initialization:: *)
hilbertSchmidt[ c_ A_] := Abs[c] hilbertSchmidt[A] /; matrix[A]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], v_List][f_]:=(unmakeVector[v];Plus@@(partial[f,{#1,2}]&)/@v)


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][f_]:=(makeVector[x];0)/;FreeQ[f,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][Subscript[x_, j_]]:=(makeVector[x];0)


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][HoldPattern[+a__]]:=Plus@@Thread[Subscript[\[CapitalDelta], x][{a}]]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][f_ g_]:=f Subscript[\[CapitalDelta], x][g]+Subscript[\[CapitalDelta], x][f] g+2 gradient[f,x] . gradient[g,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][f_^p_]:=p f^(p-1) Subscript[\[CapitalDelta], x][f]+p (p-1) f^(p-2) gradient[f,x] . gradient[f,x] +2 p f^(p-1) Log[f] gradient[f,x] . gradient[p,x]+f^p Log[f] Subscript[\[CapitalDelta], x][p]+2 f^(p-1) gradient[f,x] . gradient[p,x]+f^p Log[f]^2 gradient[p,x] . gradient[p,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][f_[a_. x_+b_.]]:=(makeVector[x];If[b==0,Null,makeVector[b]];a^2 \[CapitalDelta][f][a x+b])/;FreeQ[a,x]&&FreeQ[b,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][a_ . x_]:=(makeVector[x];0)/;FreeQ[a,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][x_ . y_]:=(makeVector[x];0)/;FreeQ[y,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][(A_ . x_) . y_]:=(makeVector[x];0)/;FreeQ[A,x]&&FreeQ[y,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][(x_ . Y_) . z_]:=(makeVector[x];0)/;FreeQ[Y,x]&&FreeQ[z,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][a_ . (B_ . x_)]:=(makeVector[x];0)/;FreeQ[a,x]&&FreeQ[B,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][a_ . (x_ . Y_)]:=(makeVector[x];0)/;FreeQ[a,x]&&FreeQ[Y,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][(A_ . x_) . x_]:=(makeVector[x];2 trace[A])/;FreeQ[A,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][(x_ . A_) . x_]:=(makeVector[x];2 trace[A])/;FreeQ[A,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][x_ . (A_ . x_)]:=(makeVector[x];2 trace[A])/;FreeQ[A,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][x_ . (x_ . A_)]:=(makeVector[x];2 trace[A])/;FreeQ[A,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], b_Symbol][b_ . gradient[u_][b_]]:=(makeVector[b];2 Subscript[\[CapitalDelta], b][u[b]]+b . gradient[Subscript[\[CapitalDelta], b][u[b]],b])


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][gradient[u_][x_] . x_]:=(makeVector[x];2 Subscript[\[CapitalDelta], x][u[x]]+x . gradient[Subscript[\[CapitalDelta], x][u[x]],x])


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][f_[g_]]:=Derivative[2][f][g]  gradient[g,x] . gradient[g,x] + Derivative[1][f][g] Subscript[\[CapitalDelta], x][g]/;(makeVector[x];vector[g]=!=True)


(* ::Input::Initialization:: *)
\[CapitalDelta][norm]:=(dimension[#1]-1)/norm[#1]&


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][norm[A_ . x_]]:=(makeVector[x];(norm[A . x]^2 hilbertSchmidt[A]^2-norm[A . x . A]^2)/norm[A . x]^3)/;FreeQ[A,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][norm[x_ . Z_]]:=(makeVector[x];(norm[x . Z]^2 hilbertSchmidt[Z]^2-norm[x . Z . Transpose[Z]]^2)/norm[x . Z]^3)/;FreeQ[Z,x]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][Gamma[t_, f_]] :=E^-f f^(t-1) norm[ gradient[f, x] ]^2-E^-f f^(t-1) Subscript[\[CapitalDelta], x][f]  - (t-1)E^-f f^(t-2) norm[ gradient[f, x] ]^2 /; FreeQ[t, x]


(* ::Input::Initialization:: *)
\[CapitalDelta][dimension]:=0&


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol,y_Symbol][f_]:=(unmakeVector[y];Subscript[\[CapitalDelta], x][f]+partial[f,{y,2}])


(* ::Input::Initialization:: *)
Unprotect[Power]; (Subscript[\[CapitalDelta], x__]^n_)[f_]:=(Subscript[\[CapitalDelta], x]^(n-1))[Subscript[\[CapitalDelta], x][f]]/;Head[n]==Integer&&Positive[n]


(* ::Input::Initialization:: *)
Subscript[\[CapitalDelta], x_Symbol][u_Symbol[norm[x_]^m_ x_]] := (makeVector[x];norm[x]^(m-2) (norm[x]^(2+m) \[CapitalDelta][u][norm[x]^m x] + m(m + dimension[x]) x . gradient[u][norm[x]^m x]+ m(m+2)norm[x]^m (jacobian[gradient[u]][norm[x]^m x] . x) . x))


(* ::Input::Initialization:: *)
laplacian = \[CapitalDelta]


(* ::Input::Initialization:: *)
dimension[a_ x_]:=dimension[x]/;vector[x]||matrix[x]


(* ::Input::Initialization:: *)
dimension[a_?NumericQ x_]:=(If[matrix[x]=!=True,makeVector[x]];dimension[x])


(* ::Input::Initialization:: *)
dimension[a_. southPole+x_]:=dimension[x]


(* ::Input::Initialization:: *)
dimension[x_+w_]:=If[StringLength[ToString[dimension[x]]]<StringLength[ToString[dimension[w]]],dimension[x],dimension[w]]


(* ::Input::Initialization:: *)
dimension[v_?VectorQ]:=Length[v]


(* ::Input::Initialization:: *)
dimension[A_?MatrixQ]:={Length[A],Length[A[[1]]]}


(* ::Input::Initialization:: *)
IntegerQ[dimension[_]]^:=True


(* ::Input::Initialization:: *)
setDimension::Matrix = "`1` will be considered to be a matrix of size `2`-by-`3`."


(* ::Input::Initialization:: *)
setDimension[x_Symbol,{m_,n_}]:=(dimension[x]={m,n};makeMatrix[x];If[Not[IntegerQ[m]],Unprotect[IntegerQ];IntegerQ[m]=True];If[Not[IntegerQ[n]],Unprotect[IntegerQ];IntegerQ[n]=True];
			With[{messageOptions = Options[$FrontEnd, MessageOptions][[1]]}, SetOptions[$FrontEnd, MessageOptions->{Global`WarningAction->{}}];Message[setDimension::Matrix, x, m, n]; SetOptions[$FrontEnd, messageOptions];])


(* ::Input::Initialization:: *)
setDimension::Vector = "`1` will be considered to be a vector in `2`-dimensional real Euclidean space."


(* ::Input::Initialization:: *)
setDimension2[x_Symbol,n_]:=(dimension[x]=n;makeVector[x];If[Not[IntegerQ[n]],Unprotect[IntegerQ];IntegerQ[n]=True];)


(* ::Input::Initialization:: *)
setDimension[x_Symbol,n_]:=(setDimension2[x, n];With[{messageOptions = Options[$FrontEnd, MessageOptions][[1]]}, SetOptions[$FrontEnd, MessageOptions->{Global`WarningAction->{}}];Message[setDimension::Vector, x, n]; SetOptions[$FrontEnd, messageOptions];])


(* ::Input::Initialization:: *)
setDimension::VectorValuedFunction = "`1` will be considered to be a function taking values in `2`-dimensional Euclidean space."


(* ::Input::Initialization:: *)
setDimension[f_Symbol[_],n_]:=(dimension[f[_]]=n;makeVector[f[_]];If[Not[IntegerQ[n]],Unprotect[IntegerQ];IntegerQ[n]=True];
		With[{messageOptions = Options[$FrontEnd, MessageOptions][[1]]}, SetOptions[$FrontEnd, MessageOptions->{Global`WarningAction->{}}];Message[setDimension::VectorValuedFunction, f, n]; SetOptions[$FrontEnd, messageOptions];])/;f=!=List


(* ::Input::Initialization:: *)
homogeneous[f_,d_,x_Symbol,c_Symbol]:=(makeVector[x];makeVector[c];setDimension2[c,dimension[x]];homogeneous[f,d,Union[Cases[Level[f,-1],Subscript[x, _]]],Union[Cases[Level[f,-1],Subscript[x, _]]]/.x->c]/.Table[ Subscript[c, j], {j, Length[Union[Cases[Level[f,-1],Subscript[x, _]]]]}]->c)/;!(MemberQ[Level[f,{-1}],x])&&Head[d]==Integer


(* ::Input::Initialization:: *)
homogeneous[f_,d_,x_Symbol,c_Symbol]:=(makeVector[x];makeVector[c];setDimension2[c,dimension[x]];homogeneous[f,d,x,Table[Subscript[c, j], {j, dimension[x]}]]/.Table[Subscript[c, j], {j, dimension[x]}]->c)/;MemberQ[Level[f,{-1}],x]&&Head[dimension[x]]==Integer&&Head[d]==Integer


(* ::Input::Initialization:: *)
homogeneous[f_,d_,x_Symbol,c_List]:=(makeVector[x];unmakeVector[c];Module[{z},homogeneous[f,d,Table[Subscript[x, j], {j,Length[c]}],c]/.Subscript[x, j_]->Subscript[z, j]/.x->c/.Subscript[z, j_]->Subscript[x, j]])/;Head[d]==Integer


(* ::Input::Initialization:: *)
homogeneous[f_,d_,v_List,c_Symbol]:=(makeVector[c];unmakeVector[v];setDimension2[c,Length[v]];homogeneous[f,d,v,Table[Subscript[c, j], {j, dimension[c]}]])/;Head[d]==Integer


(* ::Input::Initialization:: *)
homogeneous[f_,d_,v_List,c_List]:=(unmakeVector[v];unmakeVector[c];If[togetherness2&&!(norm[c]===0),turnOff[togetherness]];With[{mI=multiIndices[d,Length[v]]},Plus@@(((partial@@Prepend[Transpose[{v,mI[[#1]]}],f]/.Thread[v->c]) Times@@((v-c)^mI[[#1]]))/multiFactorial[mI[[#1]]]&)/@Range[Length[mI]]])/;Head[d]==Integer


(* ::Input::Initialization:: *)
homogeneous[f_,d_,v_List]:=homogeneous[f,d,v,Table[0,{Length[v]}]]/;Head[d]==Integer


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


(* ::Input::Initialization:: *)
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


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


(* ::Input::Initialization:: *)
multiIndices[k_,1]:={{k}}


(* ::Input::Initialization:: *)
multiIndices[0,n_]:={Table[0,{n}]}


(* ::Input::Initialization:: *)
multiIndices[k_,n_]:=multiIndices[k,n]=Join@@Table[(Prepend[#1,k-j]&)/@multiIndices[j,n-1],{j,0,k}]


(* ::Input::Initialization:: *)
multiIndices[k_,1, {t_}]:=If[Mod[k-t,2]==0,{{k}},{}]


(* ::Input::Initialization:: *)
multiIndices[0,n_, tt_]:=If[And @@ EvenQ[tt],{Table[0,{n}]},{}]


(* ::Input::Initialization:: *)
multiIndices[k_,n_, tt_]:=multiIndices[k,n, tt]=Join@@Table[(Prepend[#1,k-j]&)/@multiIndices[j,n-1, Delete[tt,1]],{j,Mod[k-tt[[1]],2],k, 2}]


(* ::Input::Initialization:: *)
taylor[f_,d_,vc___]:=Plus@@(homogeneous[f,#1,vc]&)/@Range[0,d]/;Head[d]==Integer


(* ::Input::Initialization:: *)
volume[n_Integer]:=\[Pi]^(n/2)/Gamma[n/2+1]/;Positive[n]


(* ::Input::Initialization:: *)
surfaceArea[n_Integer]:=(2 \[Pi]^(n/2))/Gamma[n/2]/;Positive[n]


(* ::Input::Initialization:: *)
integrateBall[f_,v_List]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];integrateBall[f/.Thread[v->Table[Subscript[z, j], {j, Length[v]}]],z]])


(* ::Input::Initialization:: *)
integrateBall[f_,x_Symbol]:=(makeVector[x];
		f volume[dimension[x]])/;FreeQ[f,x]


(* ::Input::Initialization:: *)
integrateBall[f_,x_Symbol]:=\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(-1\), \(1\)]\(\((f /. {
\*SuperscriptBox[\(norm[x]\), \(2\)] -> 
\*SuperscriptBox[\(x\), \(2\)], norm[x] -> Abs[x], 
\*SubscriptBox[\(x\), \(1\)] -> x})\) \[DifferentialD]x\)\)/;dimension[x]==1


(* ::Input::Initialization:: *)
integrateBall[f_,x_Symbol]:=(makeVector[x];Module[{r},If[Head[dimension[x]]===Integer,Null,Unprotect[Limit]; Limit[c_. r^(dimension[x]+m_.),r->0]:=0/;(m>0||m==0)&&FreeQ[c,r]];dimension[x] volume[dimension[x]] Integrate[Expand[r^(dimension[x]-1) integrateSphere[f/.{norm[x]->r,Subscript[x, j_]->r Subscript[x, j],x->r x},x]],{r,0,1}, Assumptions->dimension[x]>=1]])


(* ::Input::Initialization:: *)
integrateSphere[f_,v_List]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];integrateSphere[f/.Thread[v->Table[ Subscript[z, j], {j,Length[v]}]],z]])


(* ::Input::Initialization:: *)
integrateSphere[f_,x_Symbol]:=(makeVector[x];integrateSphere1[Expand[f/.{norm[x]->1,x . z_:>Table[Subscript[x, j],{j, dimension[x]}] . Table[Subscript[z, j],{j, dimension[x]}],z_ . x:>Table[Subscript[z, j],{j, dimension[x]}] . Table[Subscript[x, j],{j, dimension[x]}]}],x])


(* ::Input::Initialization:: *)
integrateSphere1[HoldPattern[+a__],x_]:=Plus@@Thread[integrateSphere1[{a},x]]


(* ::Input::Initialization:: *)
integrateSphere1[f_,x_]:=f/;FreeQ[f,x]


(* ::Input::Initialization:: *)
integrateSphere1[c_. Subscript[x_, _],x_]:=0


(* ::Input::Initialization:: *)
integrateSphere1[c_ f_,x_]:=c integrateSphere1[f,x]/;FreeQ[c,x]


(* ::Input::Initialization:: *)
integrateSphere1[Subscript[x_, _]^j_,x_]:=integrateSphere3[{j},dimension[x]]


(* ::Input::Initialization:: *)
integrateSphere1[f_,x_]:=integrateSphere2[f,x]/;monomialSpecialQ[f,x]


(* ::Input::Initialization:: *)
integrateSphere2[f_,x_]:=integrateSphere3[List@@f/.Subscript[x, _]^b_->b,dimension[x]]


(* ::Input::Initialization:: *)
integrateSphere3[b_,n_]:=0/;Or@@OddQ[b]


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


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


(* ::Input::Initialization:: *)
integrateEllipsoidArea[f_,b__,v_List]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];integrateEllipsoidArea[f/.Thread[v->Table[ Subscript[z, j], {j,Length[v]}]],b,z]])


(* ::Input::Initialization:: *)
integrateEllipsoidArea[f_,b_,c_,d_,x_Symbol]:= integrateEllipsoidArea[ f /. Table[ Subscript[x, j]->Subscript[x, j] -c[[j]]/(2b[[j]]),{j,dimension[x]}], b, d - Sum[ c[[j]]^2/b[[j]], {j,dimension[x]}]/4, x]


(* ::Input::Initialization:: *)
integrateEllipsoidArea[f_,b_,c_List,x_Symbol]:= integrateEllipsoidArea[f, b, c, -1, x]


(* ::Input::Initialization:: *)
integrateEllipsoidArea[f_,b_,d_,x_Symbol]:= -d integrateEllipsoidArea[f, -b/d, x]


(* ::Input::Initialization:: *)
integrateEllipsoidArea[f_,b_,x_Symbol]:=(makeVector[x];integrateEllipsoid1[Expand[f/.{norm[Sqrt[b] Table[Subscript[x, j],{j, dimension[x]}]]->1,x . z_:>Table[Subscript[x, j],{j, dimension[x]}] . Table[Subscript[z, j],{j, dimension[x]}],z_ . x:>Table[Subscript[z, j],{j, dimension[x]}] . Table[Subscript[x, j],{j, dimension[x]}],norm[x]^k_-> Sum[Subscript[x, j]^2,{j,1, dimension[x]}]^(k/2)}],b,x])


(* ::Input::Initialization:: *)
integrateEllipsoid1[HoldPattern[+a__],b_,x_]:=Plus@@Table[ integrateEllipsoid1[{a}[[j]],b, x],{j, 1,Length[{a}]}]


(* ::Input::Initialization:: *)
integrateEllipsoid1[f_, b_,x_]:= f dimension[x] volume[dimension[x]]/(2 Sqrt[Product[b[[j]],{j, dimension[x]}]])/;FreeQ[f,x]


(* ::Input::Initialization:: *)
integrateEllipsoid1[f_,b_,x_]:=0/;Or@@OddQ[Table[Exponent[f,Subscript[x, j ]],{j,dimension[x]}]]


(* ::Input::Initialization:: *)
integrateEllipsoid1[c_ f_,b_,x_]:=c integrateEllipsoid1[f,b, x]/;FreeQ[c,x]


(* ::Input::Initialization:: *)
integrateEllipsoid1[f_,b_,x_]:=dimension[x] volume[dimension[x]] Product[ (Exponent[f,Subscript[x, j]] - 1)!!, {j, dimension[x]}]/(2 Sqrt[Product[ b[[j]]^(Exponent[f,Subscript[x, j]]+1), {j, dimension[x]}]]Product[ k, {k, dimension[x], dimension[x] + Sum[Exponent[f,Subscript[x, j]],{j,dimension[x]}]-2, 2}])


(* ::Input::Initialization:: *)
integrateEllipsoidVolume[f_,b__,v_List]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];integrateEllipsoidVolume[f/.Thread[v->Table[ Subscript[z, j], {j,Length[v]}]],b,z]])


(* ::Input::Initialization:: *)
integrateEllipsoidVolume[f_,b_,c_,d_,x_Symbol]:= integrateEllipsoidVolume[ f /. Table[ Subscript[x, j]->Subscript[x, j] -c[[j]]/(2b[[j]]),{j,dimension[x]}], b, d - Sum[ c[[j]]^2/b[[j]], {j,dimension[x]}]/4, x]


(* ::Input::Initialization:: *)
integrateEllipsoidVolume[f_,b_,c_List,x_Symbol]:= integrateEllipsoidVolume[f, b, c, -1, x]


(* ::Input::Initialization:: *)
integrateEllipsoidVolume[f_,b_,d_,x_Symbol]:= integrateEllipsoidVolume[f, -b/d, x]


(* ::Input::Initialization:: *)
integrateEllipsoidVolume[f_,b_,x_Symbol]:=(makeVector[x];integrateEllipsoid1V[Expand[f/.{x . z_:>Table[Subscript[x, j],{j, dimension[x]}] . Table[Subscript[z, j],{j, dimension[x]}],z_ . x:>Table[Subscript[z, j],{j, dimension[x]}] . Table[Subscript[x, j],{j, dimension[x]}],norm[x]^k_-> Sum[Subscript[x, j]^2,{j,1, dimension[x]}]^(k/2)}],b,x])


(* ::Input::Initialization:: *)
integrateEllipsoid1V[f_, b_, x_] := integrateBall[f /. Table[Subscript[x, j]->Subscript[x, j]/Sqrt[b[[j]]], {j, dimension[x]}],x]/Sqrt[Times @@ b]


(* ::Input::Initialization:: *)
Options[dirichlet]={region->Sphere, \[CapitalDelta]->0}


(* ::Input::Initialization:: *)
dirichlet[p_, x_Symbol, lastCoordinate_, options___] := dirHalfSpace[p /. lastCoordinate->0,x,lastCoordinate, options]/; And[ Head[lastCoordinate] =!=  Rule, Head[lastCoordinate] =!= List]


(* ::Input::Initialization:: *)
dirHalfSpace[p_,x_,lastCoordinate_] := Module[{d = If[degree[p,x]==-\[Infinity],0,degree[p,x]/2],pL=p,dirHalf = p,k}, Do[ pL =Expand[(-1)/((2k-1)2k) Subscript[\[CapitalDelta], x][pL]];dirHalf = dirHalf + lastCoordinate^(2k) pL, {k,d}]; dirHalf]
(* The Expand command above greatly speeds up this calculation. *)


(* ::Input::Initialization:: *)
dirHalfSpace[p_,x_,lastCoordinate_, options_]:= Module[ {t = antiLaplacian[(\[CapitalDelta] /. {options} /. Options[dirichlet]),x, singularity -> 0]}, t + dirHalfSpace[p-t,x,lastCoordinate]]


(* ::Input::Initialization:: *)
dirichlet[p_, x_Symbol,options__] := With[ {t = antiLaplacian[(\[CapitalDelta] /. {options} /. Options[dirichlet]),x, singularity -> 0]},Factor[t+dirichlet[p-t,x, Sequence @@ DeleteCases[{options}, \[CapitalDelta]->_]]]]/; And[(region /. {options} /. Options[dirichlet]) === exteriorSphere, MemberQ[{options}, \[CapitalDelta]-> _]]


(* ::Input::Initialization:: *)
dirichlet[p_, x_Symbol,options__] := With[ {t = antiLaplacian[(\[CapitalDelta] /. {options} /. Options[dirichlet]),x]}, t+dirichlet[p-t,x, Sequence @@ DeleteCases[{options}, \[CapitalDelta]->_]]] /; MemberQ[{options}, \[CapitalDelta]-> _]


(* ::Input::Initialization:: *)
dirichlet[p_, x_Symbol,options___] := With[ {region2 = (region /. {options} /. Options[dirichlet])}, Which[region2 === Sphere, dirichlet1[p,x],
Head[region2]===quadratic, p + antiLaplacian[-Subscript[\[CapitalDelta], x][p], x, multiple-> region2],
region2 === exteriorSphere, exteriorDirichlet[p,x],
region2 === HalfSpace, dirichlet[p /. norm[x]^m_-> Sum[Subscript[x, j]^2,{j,1, dimension[x]}]^(m/2), x, Subscript[x, dimension[x] ]],
Head[region2] === annulus, If[Head[p]===List, annularDirichlet[p[[1]],p[[2]],region2[[1]], region2[[2]],x], annularDirichlet[p,p,region2[[1]], region2[[2]],x]]]]


(* ::Input::Initialization:: *)
dirichlet1[p_,x_]:=Plus@@(#1[[1]]&)/@harmonicDecomposition[p/.norm[x]->1,x] /; PolynomialQ[p /. norm[x]->1, Cases[p, Subscript[x, _],\[Infinity]]]


(* ::Input::Initialization:: *)
dirichlet[f_,v_List, opt___]:=explicit[dirichlet,f,v, opt] 


(* ::Input::Initialization:: *)
exteriorDirichlet[f_,z_]:=Factor[kelvin[dirichlet[f,z],z]]


(* ::Input::Initialization:: *)
annular[0,r_,s_,z_]=0


(* ::Input::Initialization:: *)
annular[f_,r_,s_,z_]:=Plus@@(If[#1[[2]]===#1[[3]],((Log[norm[z]]-Log[s]) r^#1[[2]] #1[[1]])/(Log[r]-Log[s]),((norm[z]^(2 #1[[2]]-2 #1[[3]])-s^(2 #1[[2]]-2 #1[[3]])) r^#1[[2]] #1[[1]])/(r^(2 #1[[2]]-2 #1[[3]])-s^(2 #1[[2]]-2 #1[[3]]))]&)/@Flatten[(append[harmonicDecomposition[#1[[1]],z],#1[[2]]]&)/@polynomialDecomposition[f,z],1]/;dimension[z]==2


(* ::Input::Initialization:: *)
annular[f_,r_,s_,z_]:=Plus@@(((norm[z]^(2-dimension[z]+2 #1[[2]]-2 #1[[3]])-s^(2-dimension[z]+2 #1[[2]]-2 #1[[3]])) r^#1[[2]] #1[[1]])/(r^(2-dimension[z]+2 #1[[2]]-2 #1[[3]])-s^(2-dimension[z]+2 #1[[2]]-2 #1[[3]]))&)/@Flatten[(append[harmonicDecomposition[#1[[1]],z],#1[[2]]]&)/@polynomialDecomposition[f,z],1]


(* ::Input::Initialization:: *)
append[b_,d_]:=(Append[#1,d]&)/@b


(* ::Input::Initialization:: *)
annularDirichlet[f_,g_,r_,s_,z_Symbol]:=annular[f,r,s,z]+annular[g,s,r,z]


(* ::Input::Initialization:: *)
d[0,0]=1


(* ::Input::Initialization:: *)
lap[f_,v_]:=Plus@@(\!\(
\*SubscriptBox[\(\[PartialD]\), \({#1, 2}\)]f\)&)/@v


(* ::Input::Initialization:: *)
harmonicDecomposition[0,x_]:={}


(* ::Input::Initialization:: *)
harmonicDecomposition[p_,x_Symbol]:=(makeVector[x];With[{deg=degree[p/.norm[x]->normtemp,x],n=dimension[x]},powers={};hd[_]:=0;With[{partition=Partition[Flatten[Expand[((m=#1[[2]];k=Floor[m/2];If[m>1,d[1,1]=1/(2 (2 m+n-4))];v=Select[Variables[#1[[1]]],!(FreeQ[#1,x])&];laplac[0]=#1[[1]];laplac[1]=lap[#1[[1]],v];Do[d[j,j]=(d[j-1,j-1] (2 m+n-2 j))/(2 j (2 m+n-4 j+2) (2 m+n-4 j));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[{\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(j = i\), \(k\)]\(d[i, j]\ 
\*SuperscriptBox[\(norm[x]\), \(2\ \((j - i)\)\)]\ laplac[j]\)\),2 i},{i,0,k}])&)/@polynomialDecomposition[p/.norm[x]->normtemp,x]]//.{a_.+normtemp^e_. b_.,j_}->{{a,j},{b,j+e}}],2]},powers=powers\[Union]Transpose[partition][[2]];((hd[#1[[2]]]=hd[#1[[2]]]+#1[[1]])&)/@partition];aa=Select[({hd[#1],#1}&)/@powers,#1[[1]]=!=0&];Clear[hd];aa])


(* ::Input::Initialization:: *)
harmonicDecomposition[f_,v_List]:=explicit[harmonicDecomposition,f,v]


(* ::Input::Initialization:: *)
Options[antiLaplacian]={multiple->None,singularity->None}


(* ::Input::Initialization:: *)
antiLaplacian[0,x__,OptionsPattern[]]:=0


(* ::Input::Initialization:: *)
antiLaplacian[0.,x__,OptionsPattern[]]:=0.


(* ::Input::Initialization:: *)
antiLaplacian[p_,x_Symbol,OptionsPattern[]]:=(makeVector[x];
If[Head[OptionValue[multiple]]===quadratic,
antiLaplacianQ[p, x, OptionValue[multiple]],
If[OptionValue[multiple]=== norm^2||!(FreeQ[p,norm[x]]),antiLaplacianN[Expand[p],x,OptionValue[singularity]],antiLaplacian1[Expand[p],x]]])


(* ::Input::Initialization:: *)
antiLaplacian[p_,x_Symbol,lastCoordinate_,OptionsPattern[]]:= Module[{t =  p, a =antiLaplacian[p, x, singularity->OptionValue[singularity]],s, k},s = a; Do[t = Subscript[\[CapitalDelta], x,lastCoordinate][a] - t; a = antiLaplacian[ t, x]; s =s +(-1)^k a,{k,1,degree[p,lastCoordinate]/2}];s]/; PolynomialQ[p,lastCoordinate]


(* ::Input::Initialization:: *)
antiLaplacian[f_,v_List,OptionsPattern[]]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];antiLaplacian[f/.norm[v]^2->norm[z]^2/.Thread[v->Table[Subscript[z, j], {j, Length[v]}]],z,{multiple-> OptionValue[multiple],singularity->OptionValue[singularity]}]/.Prepend[Thread[Table[Subscript[z, j], {j, Length[v]}]->v],norm[z]->norm[v]]])


(* ::Input::Initialization:: *)
antiLaplacianQ[p_, x_, q_ ] := 
Module[{k,j,pp,ph,m},pp = (p /. norm[x]^m_?EvenQ-> Sum[Subscript[x, k]^2,{k, dimension[x]}]^(m/2));ph = blockPolynomialDecomposition[pp, x];
Which[ Length[q] ==1, antiLaplacianQ[pp,x,quadratic[q[[1]],-1]],
Length[q]==2 , If[Head[q[[2]]]===List,antiLaplacianQ[pp,x,quadratic[q[[1]],q[[2]],-1]],(q[[1]] . Array[x,Length[q[[1]]]]^2+  q[[2]])Sum[antiLaplacianQ4[ph[[j]], x, List @@ q], {j,1, Length[ph]}]],
Length[q]==3,(q[[1]] . Array[x,Length[q[[1]]]]^2+ q[[2]] . Array[x,Length[q[[2]]]]+ q[[3]])Sum[antiLaplacianQ3[ph[[j]], x, List @@ q], {j,1, Length[ph]}]]]


(* ::Input::Initialization:: *)
antiLaplacianQ4[p_, x_, q_ ] := antiLaplacianQ2[p, x, q[[1]]] /; q[[2]] ===0


(* ::Input::Initialization:: *)
antiLaplacianQ4[p_, x_, q_ ] := Module[{f},With[{m = p[[2]], q2 = q[[1]], q0 = q[[2]], p3= p[[3]]},
f[m+1] = 0;
f[m] = antiLaplacianQ2[p, x, q2];
Do[f[m-j] =  antiLaplacianQ2[{-Subscript[\[CapitalDelta], x][q0 f[m-j+2]],m-j,p3},x,q2],{j,1,m} ];
Sum[f[j], {j,0, m}]
]]


(* ::Input::Initialization:: *)
antiLaplacianQ3[p_, x_, q_ ] := antiLaplacianQ4[p, x, {q[[1]],q[[3]]}] /;  q[[2]] === Array[0,dimension[x]]


(* ::Input::Initialization:: *)
antiLaplacianQ3[p_, x_, q_ ] := Module[{f},With[{m = p[[2]], q2 = q[[1]], q1 = q[[2]] . Array[x,Length[q[[2]]]], q0 = q[[3]]},
f[m+1] = 0;
f[m] = antiLaplacianQ2[p, x, q2];
Do[f[m-j] =  antiLaplacianQ2a[-Subscript[\[CapitalDelta], x][q1 f[m-j+1] + q0 f[m-j+2]],x,q2],{j,1,m} ];
Sum[f[j], {j,0, m}]
]]


(* ::Input::Initialization:: *)
antiLaplacianQ2a[0, x_, q_] := 0


(* ::Input::Initialization:: *)
antiLaplacianQ2a[p_, x_, q_ ] := 
With[{ph = blockPolynomialDecomposition[p, x]},Sum[antiLaplacianQ2[ph[[j]], x, q], {j,1, Length[ph]}]
]


(* ::Input::Initialization:: *)
antiLaplacianQ2[p_, x_, q_] :=  p[[1]]/(2 Plus @@ q) /; p[[2]] == 0


(* ::Input::Initialization:: *)
antiLaplacianQ2[p_, x_, q_] :=  With[{p1=p[[1]], twoA = 2(Plus @@ q), n = dimension[x], multiIndicesmn =multiIndices[p[[2]],dimension[x], p[[3]]]},Module[{d,s},((Plus @@ Map[Subscript[d, #] (Times @@ (Array[x,n]^#))/multiFactorial[#]&, multiIndicesmn]) /. Solve[Map[(Expand[(twoA + 4  Plus @@ (# q)) Subscript[d, #] + ((Plus @@ (# (# - 1) q Table[Subscript[d, ReplacePart[#, #[[i]]-2,i]],{i, n}]))/. Subscript[d, s_ ]:> Sum[Subscript[d, ReplacePart[s, s[[j]]+2,j]],{j,n}])] == multiFactorial[#] Coefficient[p1,Times @@ (Array[x,n]^#)])&, multiIndicesmn], Map[Subscript[d, #]&, multiIndicesmn]])[[1]]]]


(* ::Input::Initialization:: *)
antiLaplacianN[HoldPattern[+a__],z_,singularityValue_]:=Plus@@Thread[Unevaluated[antiLaplacianN[{a},z,singularityValue]]]


(* ::Input::Initialization:: *)
antiLaplacianN[c_ g_,z_,singularityValue_]:=c antiLaplacianN[g,z,singularityValue]/;FreeQ[c,z]


(* ::Input::Initialization:: *)
antiLaplacianN[g_,z_,singularityValue_]:=antiLap1[If[Head[g]===Times,Select[g,FreeQ[#1,norm[z]]&],If[FreeQ[g,norm[z]],g,1]],If[Head[g]===Times,Select[g,!(FreeQ[#1,norm[z]])&],If[FreeQ[g,norm[z]],1,g]],z,singularityValue]


(* ::Input::Initialization:: *)
antiLap1[p_,f_,z_,singularityValue_]:=Plus@@(antiLap2[#1[[1]],#1[[2]],f,z,singularityValue]&)/@harmonicDecomposition[p,z]


(* ::Input::Initialization:: *)
antiLap2[p_,j_,f_,z_,singularityValue_]:=Plus@@(antiLap3[#1[[1]],#1[[2]],j,f,z,singularityValue]&)/@polynomialDecomposition[p,z]


(* ::Input::Initialization:: *)
antiLap3[p_,m_,j_,f_,z_,singularityValue_]:=With[{n=dimension[z]},p If[singularityValue=!=None,\[Integral]t^(1-2 m-n) (\[Integral](t^(-1+2 m+n+j) f/.norm[z]->t)\[DifferentialD]t)\[DifferentialD]t,With[{int=\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(t\)]\(\((
\*SuperscriptBox[\(s\), \(\(-1\) + 2\ m + n + j\)]\ f /. norm[z] -> s)\) \[DifferentialD]s\)\)},\[Integral]t^(1-2 m-n) If[FreeQ[int,Indeterminate]&&FreeQ[int,\[Infinity]]&&FreeQ[int,ConditionalExpression] && FreeQ[int,Limit]&&FreeQ[int,If]&&FreeQ[int,Integrate]&&singularityValue===None,int,\[Integral](t^(-1+2 m+n+j) f/.norm[z]->t)\[DifferentialD]t]\[DifferentialD]t]]/.t->norm[z]]


(* ::Input::Initialization:: *)
antiLaplacian1[HoldPattern[+a__],z_]:=Plus@@Thread[antiLaplacian1[{a},z]]


(* ::Input::Initialization:: *)
antiLaplacian1[g_,z_]:=1/2 g Subscript[z, 1]^2/;FreeQ[g,z]


(* ::Input::Initialization:: *)
antiLaplacian1[g_,z_]:=antiLaplacian3[Transpose[Sort[({Exponent[g,Subscript[z, #1]],#1}&)/@(#1[[2]]&)/@Variables[g]]],g,z]/;monomialQ[g,z]


(* ::Input::Initialization:: *)
antiLaplacian1[c_ g_,z_]:=c antiLaplacian1[g,z]/;FreeQ[c,z]


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


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


(* ::Input::Initialization:: *)
neumann[f_,v_List]:=explicit[neumann,f,v]


(* ::Input::Initialization:: *)
neumann[f_,g_,v_List]:=explicit[neumann,f,g,v]


(* ::Input::Initialization:: *)
neumann[0,x_Symbol]:=0


(* ::Input::Initialization:: *)
neumann[f_,x_Symbol]:=Plus@@(If[#1[[2]]===0,0,#1[[1]]/#1[[2]]]&)/@polynomialDecomposition[dirichlet[f,x],x]


(* ::Input::Initialization:: *)
neumann[f_,g_,x_Symbol]:=With[{v=antiLaplacian[g,x]},v+neumann[f-normalD[v,x],x]-Limit[v/.{Subscript[x, _]->0,norm[x]->t},t->0,Direction->-1]]


(* ::Input::Initialization:: *)
antiLaplacianQNa[0,x_Symbol, q_ ] := 0


(* ::Input::Initialization:: *)
antiLaplacianQNa[p_, x_Symbol, q_ ] := 
With[{ph = blockPolynomialDecomposition[p, x]},Sum[antiLaplacianQN[ph[[j]], x, q], {j,1, Length[ph]}]
]


(* ::Input::Initialization:: *)
makePD[l_,x_] :=
MapIndexed[{Subscript[x, First[#2]],#1}&,l]


(* ::Input::Initialization:: *)
antiLaplacianQN[p_, x_Symbol, q_] :=  With[{p1=p[[1]], n = dimension[x], multiIndicesmn =multiIndices[p[[2]],dimension[x], p[[3]]]},Module[{c,d,r}, 
c := (q / (2(q . #)+4q) Subscript[d, #])& /@ multiIndicesmn;
r:=Laplacian[Table[Times @@ (Array[x,n]^ReplacePart[#,#[[i]]+2,i]),{i,n}]& /@ multiIndicesmn,Array[x,dimension[x]]];
((Plus @@ Map[Subscript[d, #] (Times @@ (Array[x,n]^#))&, multiIndicesmn]) /.
Solve[Map[  Plus @@ Flatten[c  Apply[D,Prepend[ makePD[#,x],r]]]
==(*multiFactorial[#] Coefficient[p1,Times @@ (Array[x,n]^#)]&,multiIndicesmn]*)Apply[D,Prepend[ makePD[#,x],p1]]&,multiIndicesmn],Map[Subscript[d, #]&, multiIndicesmn]
])[[1]]]]


(* ::Input::Initialization:: *)
neumannInvMap[p_,x_Symbol,q_] :=
Plus @@(If[PossibleZeroQ[2q . #[[1]]],#[[2]],(#[[2]]/(2q . #[[1]]))](Times @@ (Array[x,dimension[x]]^#[[1]]))& /@CoefficientRules[p,Array[x,dimension[x]]])


(* ::Input::Initialization:: *)
neumann2[p_,x_Symbol,q_] :=
With[{ph = Association[#[[2]]->#[[1]]&/@polynomialDecomposition[p,x]],n=dimension[x], qx = q . Array[x,dimension[x]]^2-1},
Module[{f,m,pp},
m = Max[Keys[ph]]-2;
f[m+2] = 0;
f[m+1] = 0;
f[-1] = 0;
Do[
pp = Simplify[Subscript[\[CapitalDelta], x][neumannInvMap[Lookup[ph,m+2-j,0]+f[m+2-j],x,q]]];
f[m-j] = antiLaplacianQNa[pp,x,q];
,{j,0,m} 
];
Sum[f[j], {j,0, m}]
]]


(* ::Input::Initialization:: *)
neumann1[p_, x_, q_] := neumannInvMap[p,x,q] - neumannInvMap[(q . Array[x, dimension[x]]^2 -1)neumann2[p,x,q],x,q]


(* ::Input::Initialization:: *)
neumann[f_,v_List, options__]:=explicit[neumann,f,v, options]


(* ::Input::Initialization:: *)
neumann[p_, x_Symbol, options__] := Module[{quad =(region /. options), neu}, 
Which[
Length[quad] == 1,neumann1[p, x, quad[[1]]],
Length[quad] ==2,If[Head[quad[[2]]]===List,neumann[p, x, region-> quadratic[quad[[1]], quad[[2]], -1]],neumann1[p, x, -quad[[1]]/quad[[2]]]/(-quad[[2]])],
Length[quad] == 3,  
neu = neumann[p /. Table[ Subscript[x, j]->Subscript[x, j] -quad[[2]][[j]]/(2quad[[1]][[j]]),{j,dimension[x]}], x, region->quadratic[quad[[1]],Simplify[quad[[3]]-Sum[quad[[2]][[j]]^2/quad[[1]][[j]],{j,dimension[x]}]/4]]] /. Table[ Subscript[x, j]->Subscript[x, j] +quad[[2]][[j]]/(2quad[[1]][[j]]),{j,dimension[x]}]; neu -( neu/. Table[Subscript[x, j]-> 0, {j, dimension[x]}])
]]


(* ::Input::Initialization:: *)
neumann[f_, g_, x_Symbol, options__] := Module[{v=antiLaplacian[g,x],quad =(region /. options),q},
Which[Length[quad] == 1, q = quad[[1]] . Array[x, dimension[x]]^2,
Length[quad] == 2, If[Head[quad[[2]]]===List, q = quad[[1]] . Array[x, dimension[x]]^2+ quad[[2]] . Array[x, dimension[x]],  q = quad[[1]] . Array[x, dimension[x]]^2],
Length[quad]==3, q = quad[[1]] . Array[x, dimension[x]]^2+ quad[[2]] . Array[x, dimension[x]]
];
v+neumann[f-gradient[v,x] . gradient[q,x],x, options]]


(* ::Input::Initialization:: *)
exteriorNeumann[f_,v_List]:=explicit[exteriorNeumann,f,v]


(* ::Input::Initialization:: *)
exteriorNeumann[f_,g_,v_List]:=explicit[exteriorNeumann,f,g,v]


(* ::Input::Initialization:: *)
exteriorNeumann[0,x_Symbol]:=0


(* ::Input::Initialization:: *)
exteriorNeumann[f_,x_Symbol]:=Plus@@((norm[x]^(2-dimension[x]-2 #1[[2]]) #1[[1]])/(dimension[x]+#1[[2]]-2)&)/@polynomialDecomposition[dirichlet[f,x],x]/;dimension[x]=!=2


(* ::Input::Initialization:: *)
exteriorNeumann[f_,x_Symbol]:=Plus@@(norm[x]^(-2 #1[[2]]) #1[[1]] If[#1[[2]]===0,0,1/#1[[2]]]&)/@polynomialDecomposition[dirichlet[f,x],x]/;dimension[x]===2


(* ::Input::Initialization:: *)
exteriorNeumann[f_,g_,x_Symbol]:=With[{v=antiLaplacian[g,x,singularity->0]},v+exteriorNeumann[f+normalD[v,x],x]]


(* ::Input::Initialization:: *)
biDirichlet[f_,v_List]:=explicit[biDirichlet,f,v]


(* ::Input::Initialization:: *)
biDirichlet[f_,g_,v_List]:=explicit[biDirichlet,f,g,v]


(* ::Input::Initialization:: *)
biDirichlet[f_,z_Symbol]:=Together[With[{p=dirichlet[f,z]},p+1/2 (1-norm[z]^2) z . gradient[p,z]]]


(* ::Input::Initialization:: *)
biDirichlet[f_,g_,z_Symbol]:=normalIntegral[dirichlet[g,z],z]+biDirichlet[f,z]


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


(* ::Input::Initialization:: *)
zonalHarmonic[0,x_,z_]:=1


(* ::Input::Initialization:: *)
zonalHarmonic[m_Integer,x_Symbol,z_Symbol]:=(makeVector[x];makeVector[z];Plus@@(1/(#1! (m-2 #1)! 2^#1) ((-1)^#1 (dimension[x+z]+2 m-2) Times@@Range[dimension[x+z],dimension[x+z]+2 m-2 #1-4,2] (x . z)^(m-2 #1) norm[x]^(2 #1) norm[z]^(2 #1))&)/@Range[0,m/2])


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


(* ::Input::Initialization:: *)
dimHarmonic[0,2]=1;


(* ::Input::Initialization:: *)
dimHarmonic[m_,n_]:=Binomial[n+m-1,n-1]-Binomial[n+m-3,n-1]


(* ::Input::Initialization:: *)
ball[f_,g_, x_]:= integrateBall[ f g, x]


(* ::Input::Initialization:: *)
sphere[f_,g_, x_]:= integrateSphere[ f g, x]


(* ::Input::Initialization:: *)
basisH[0,x_]:={1}


(* ::Input::Initialization:: *)
basisH[m_Integer,x_Symbol]:=Together[Module[{xx},With[{ttt=Expand[Expand[(kelvin[#1,xx]&)/@(partial@@Prepend[Transpose[{{Subscript[xx, 1],Subscript[xx, 2]},#1}],Log[norm[xx]]]&)/@Join[(Prepend[#1,0]&)/@multiIndices[m,1],(Prepend[#1,1]&)/@multiIndices[m-1,1]]]]/.xx->x},noLeadingMinus/@Transpose[FactorTermsList/@ttt][[2]]]]]/;dimension[x]==2


(* ::Input::Initialization:: *)
basisH[m_Integer,x_Symbol]:=Together[Module[{xx},With[{ttt=Expand[Expand[(kelvin[#1,xx]&)/@(partial@@Prepend[Transpose[{Table[ Subscript[xx, j], {j, dimension[x]}],#1}],norm[xx]^(2-dimension[x])]&)/@Join[(Prepend[#1,0]&)/@multiIndices[m,dimension[x]-1],(Prepend[#1,1]&)/@multiIndices[m-1,dimension[x]-1]]]]/.xx->x},noLeadingMinus/@Transpose[FactorTermsList/@ttt][[2]]]]]/;Head[dimension[x]]==Integer


(* ::Input::Initialization:: *)
basisH[m_Integer,x_Symbol, ip_]:= (inprod = ip; If [ip === Ball, inprod = ball]; If[ip=== Sphere, inprod = sphere];orthonormalize[ basisH[m, x], inprod[#1, #2, x] & ])/; Head[dimension[x]] == Integer


(* ::Input::Initialization:: *)
basisH[m_Integer,v_List]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];Expand[basisH[m,z]/.Thread[Table[ Subscript[z, j], {j,Length[v]}]->v]]/.norm[z]->norm[v]])


(* ::Input::Initialization:: *)
basisH[m_Integer,v_List, ip_]:= (inprod = ip; If [ip === Ball, inprod = ball]; If[ip=== Sphere, inprod = sphere];orthonormalize[ basisH[m, v], inprod[#1, #2, v] & ])


(* ::Input::Initialization:: *)
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]


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


(* ::Input::Initialization:: *)
reflection[\[Infinity]]:=0


(* ::Input::Initialization:: *)
reflection[x_]:=(If[Head[x]==Symbol,makeVector[x]];If[norm[x]===0,\[Infinity],x/norm[x]^2])


(* ::Input::Initialization:: *)
reflection[x_List, Sphere[c_, r_]]:= reflection[x, Sphere[Table[c . Subscript[\[Delta], j],{j,1,Length[x]}],r]] /; Head[c] =!= List


(* ::Input::Initialization:: *)
reflection[x_, Sphere[c_List, r_]]:= reflection[Table[x . Subscript[\[Delta], j],{j,1,Length[c]}], Sphere[c, r]] /; Head[x] =!= List


(* ::Input::Initialization:: *)
reflection[x_, Sphere[c_, r_]]:= c + r^2 reflection[x-c]


(* ::Input::Initialization:: *)
reflection[x_, Hyperplane[b_, c_]]:= x - (2((x . b)-c)b)/norm[b]^2


(* ::Input::Initialization:: *)
reflection[x_List, Hyperplane[b_, c_]]:= reflection[x, Hyperplane[Table[b . Subscript[\[Delta], j],{j,1,Length[x]}],c]] /; Head[b] =!= List


(* ::Input::Initialization:: *)
reflection[x_, Hyperplane[b_List,c_]]:= reflection[Table[x . Subscript[\[Delta], j],{j,1,Length[b]}], Hyperplane[b, c]] /; Head[x] =!= List


(* ::Input::Initialization:: *)
kelvin[HoldPattern[+u__],x_]:=Plus@@Thread[kelvin[{u},x]]


(* ::Input::Initialization:: *)
kelvin[u_,x_List]:=(unmakeVector[x];PowerExpand[Together[norm[x]^(2-Length[x]) (u/.Thread[x->reflection[x]])],Plus @@ (x^2)])


(* ::Input::Initialization:: *)
kelvin[u_,x_Symbol]:=(makeVector[x];norm[x]^(2-dimension[x]) (u/.{x->reflection[x],Subscript[x, j_]:>Subscript[x, j]/norm[x]^2}))


(* ::Input::Initialization:: *)
\[CapitalPhi][z_List]:=(unmakeVector[z];If[Last[z]===-1&&norm[Drop[z,-1]]===0,\[Infinity],2 reflection[z+Append[Table[0,{Length[z]-1}],1]]-Append[Table[0,{Length[z]-1}],1]])


(* ::Input::Initialization:: *)
\[CapitalPhi][z_]:=(If[!(z===0||z===\[Infinity]||MatchQ[z,_. southPole]),makeVector[z];If[dimension[southPole]=!= dimension[z],setDimension[southPole,dimension[z]]]];2 reflection[z-southPole]+southPole)


(* ::Input::Initialization:: *)
\[CapitalPhi][0,-1]=\[Infinity];


(* ::Input::Initialization:: *)
\[CapitalPhi][x_,y_]:= (makeVector[x];unmakeVector[y];{(2 x)/(norm[x]^2+(y+1)^2),(2 (y+1))/(norm[x]^2+(y+1)^2)-1})


(* ::Input::Initialization:: *)
norm[southPole]^=1;


(* ::Input::Initialization:: *)
southPole . z_Symbol:=-Subscript[z, dimension[z]]/;matrix[z]=!=True


(* ::Input::Initialization:: *)
a_Symbol . southPole:=-Subscript[a, dimension[a]]/;matrix[z]=!=True


(* ::Input::Initialization:: *)
makeVector[southPole];


(* ::Input::Initialization:: *)
N[Subscript[southPole, j_Real]] := Subscript[southPole, Round[j]]/; Round[j] == j;


(* ::Input::Initialization:: *)
Subscript[southPole, j_]:=-1/;j==dimension[southPole]


(* ::Input::Initialization:: *)
Subscript[southPole, j_]:=0/;j<dimension[southPole]


(* ::Input::Initialization:: *)
kelvinH[HoldPattern[+u__],x_]:=Plus@@Thread[kelvinH[{u},x]]


(* ::Input::Initialization:: *)
kelvinH[u_,x_List]:=(unmakeVector[x];2^(1/2 (Length[x]-2)) (norm[x]^2+2 Last[x]+1)^(2-Length[x]) (u/.Thread[x->\[CapitalPhi][x]]))


(* ::Input::Initialization:: *)
kelvinH[u_,z_Symbol]:=(makeVector[z];2^(1/2 (dimension[z]-2)) norm[z-southPole]^(2-dimension[z]) (u/.{z->\[CapitalPhi][z],Subscript[z, j_]:>\[CapitalPhi][z] . Subscript[\[Delta], j]}))


(* ::Input::Initialization:: *)
kelvinH[u_,x_Symbol,y_Symbol]:=(makeVector[x];unmakeVector[y];2^(1/2 (dimension[x]-1)) (norm[x]^2+(y+1)^2)^(1/2 (1-dimension[x])) (u/.{x->(2 x)/(norm[x]^2+(y+1)^2),y->(2 (y+1))/(norm[x]^2+(y+1)^2)-1,Subscript[x, j_]:>(2 Subscript[x, j])/(norm[x]^2+(y+1)^2)}))


(* ::Input::Initialization:: *)
poissonKernel[x_,z_]:=(makeVector[x];makeVector[z];(1-norm[x]^2 norm[z]^2)/(1-2 x . z+norm[x]^2 norm[z]^2)^(1/2 dimension[x+z]))


(* ::Input::Initialization:: *)
poissonKernelH[z_,w_]:=(makeVector[z];makeVector[w];(2 (z . Subscript[\[Delta], dimension[z]]+w . Subscript[\[Delta], dimension[z]]))/((dimension[z] volume[dimension[z]]) (norm[z]^2+norm[w]^2-2 z . w+4 z . Subscript[\[Delta], dimension[z]] w . Subscript[\[Delta], dimension[z]])^(dimension[z]/2)))


(* ::Input::Initialization:: *)
poissonKernelH[x_,y_,t_,u_]:=(makeVector[x];makeVector[t];unmakeVector[y];unmakeVector[u];(2 (y+u))/(((dimension[x]+1) volume[dimension[x]+1]) (norm[x-t]^2+(y+u)^2)^(1/2 (dimension[x]+1))))


(* ::Input::Initialization:: *)
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))])


(* ::Input::Initialization:: *)
bergmanKernelH[z_,w_]:=(makeVector[z];makeVector[w];With[{n=dimension[z+w]},(4 (n ((z . Subscript[\[Delta], n])^2+(w . Subscript[\[Delta], n])^2)-norm[z]^2-norm[w]^2+2 z . w+(2 n-4) z . Subscript[\[Delta], n] w . Subscript[\[Delta], n]))/(n volume[n] (norm[z]^2+norm[w]^2-2 z . w+4 z . Subscript[\[Delta], n]w . Subscript[\[Delta], n])^(1+n/2))])


(* ::Input::Initialization:: *)
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))])


(* ::Input::Initialization:: *)
harmonicConjugate[u_,x_Symbol,y_Symbol]:=(unmakeVector[x];unmakeVector[y];\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(y\)]\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]u \[DifferentialD]y\)\)-\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(x\)]\(\((
\*SubscriptBox[\(\[PartialD]\), \(y\)]u /. y -> 0)\) \[DifferentialD]x\)\))


(* ::Input::Initialization:: *)
bergmanProjection[0,x_]:=0


(* ::Input::Initialization:: *)
bergmanProjection[u_,x_Symbol]:=(makeVector[x];Plus@@(Plus@@#1&)/@(Function[j,((dimension[x]+2 #1[[2]]) #1[[1]])/(#1[[2]]+j+dimension[x])&][#1[[2]]]/@#1[[1]]&)/@({polynomialDecomposition[#1[[1]],x],#1[[2]]}&)/@({dirichlet[#1[[1]],x],#1[[2]]}&)/@polynomialDecomposition[u,x])


(* ::Input::Initialization:: *)
bergmanProjection[u_,v_List]:=explicit[bergmanProjection,u,v]


(* ::Input::Initialization:: *)
schwarz[0]=0;


(* ::Input::Initialization:: *)
schwarz[x_]:=PowerExpand[Simplify[(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},GenerateConditions->False]},Together[((1-1/t^2) (dimension[x]-1) volume[dimension[x]-1] (s-(s/.t->-t)))/(dimension[x] volume[dimension[x]])/.t->1/norm[x]]]]/.HoldPattern[+t_]:>Factor[Select[+t,FreeQ[#1,ArcTan]&]]+Factor[Select[+t, !(FreeQ[#1,ArcTan])&]]/.(-1+t_) (1+t_)->-(1-t^2))]/. {(-1+norm[x])^p_?EvenQ:>(1-norm[x])^p, Abs[1 + norm[x]]->1 + norm[x], Abs[-1 + norm[x]]->1 - norm[x]}]


(* ::Input::Initialization:: *)
Unprotect[IntegerQ];IntegerQ[m_?IntegerQ+n_?IntegerQ]=True;


(* ::Input::Initialization:: *)
IntegerQ[m_?IntegerQ n_?IntegerQ]=True;


(* ::Input::Initialization:: *)
vectorOrMatrix[c_] := vector[c] || matrix[c]


(* ::Input::Initialization::Closed:: *)
Unprotect[Times];


(* ::Input::Initialization:: *)
timesAdjust1 := \!\(\*
TagBox[
FormBox[
RowBox[{"HoldPattern", "[", 
RowBox[{"Times", "[", "a__", "]"}], "]"}],
StandardForm],
Format]\):=(With[{t=Position[Table[vectorOrMatrix[{a}[[j]]],{j, Length[{a}]}],True][[1,1]]},
Row[{If[Head[Times@@Delete[{a},t]]===Plus,
Row[{"(",First[Delete[{a},t]],")"}],Times@@Delete[{a},t]],
If[Head[{a}[[t]]]===Plus,Row[{"(",{a}[[t]],")"}],{a}[[t]]]}, "\!\(\*
StyleBox[\" \",\nFontSize->14]\)" ]])/; ( Or@@Table[vectorOrMatrix[{a}[[j]]],{j, Length[{a}]-1}])&&Not[{a}[[1]] === -1]
(* The space inside the quotes is in 14 point size. *)


(* ::Input::Initialization:: *)
timesAdjust2 := \!\(\*
TagBox[
FormBox[
RowBox[{"HoldPattern", "[", 
RowBox[{"Times", "[", 
RowBox[{
RowBox[{"-", "1"}], ",", "a__"}], "]"}], " ", "]"}],
StandardForm],
Format]\):=-(With[{t=Position[Table[vectorOrMatrix[{a}[[j]]],{j, Length[{a}]}],True][[1,1]]},Row[{If[Head[Times@@Delete[{a},t]]===Plus,Row[{"(",First[Delete[{a},t]],")"}],Times@@Delete[{a},t]],If[Head[{a}[[t]]]===Plus,Row[{"(",{a}[[t]],")"}],{a}[[t]]]}, "\!\(\*
StyleBox[\" \",\nFontSize->14]\)"]]) /; And[Or@@Table[vectorOrMatrix[{a}[[j]]],{j, Length[{a}]-1}], Length[{a}] > 1]
(* The space inside the quotes is in 14 point size. *)
(* A separate instruction is needed for expressions of the form
-expression because Mathematica stores such items as Times[-1, expression]. *)


(* ::Input::Initialization:: *)
timesAdjust1


(* ::Input::Initialization:: *)
timesAdjust2


(* ::Input::Initialization:: *)
explicit[function_,f_,v_List, opt___]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];function[f/.Thread[v->Table[ Subscript[z, j], {j,Length[v]}]],z, Sequence @@ ({opt} /.Thread[v->Table[ Subscript[z, j], {j,Length[v]}]])]/.Append[Thread[Table[ Subscript[z, j], {j,Length[v]}]->v],z->v]])


(* ::Input::Initialization:: *)
explicit[function_,f_,g_,v_List]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];function[f/.Thread[v->Table[ Subscript[z, j], {j,Length[v]}]],g/.Thread[v->Table[ Subscript[z, j], {j,Length[v]}]],z]/.Append[Thread[Table[ Subscript[z, j], {j,Length[v]}]->v],z->v]])


(* ::Input::Initialization:: *)
Protect[Evaluate[protectedWords]];


(* ::Input::Initialization:: *)
End[]


(* ::Input::Initialization:: *)
*)
EndPackage[ ]
(*
EndPackage[]


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