(* 

This Mathematica file defines a function, jIntegrateGaussian[], which computes Gaussian integrals analytically much faster than the normal Integrate[] function.  It operates mostly by simplifying the integral into a standard form and then applying the following formula recursively:

   Int[Exp[-a x^2 + b x + c],{x,\[Dash]\[Infinity],\[Infinity]}] = Sqrt[\[Pi]/a]Exp[b^2/(4a)+c]

It does its best to thread over sums and products correctly.  It assumes all parameters besides \[ImaginaryI] are real, and it assumes parameters are positive (or less than some constant, etc.) if necessary for an integral to converge.  Obviously this means it is much more likely to produce an incorrect output than Integrate[], especially in the form of an answer when no answer is well defined.  This is the sacrafice made for speed.

The first parameter is the integrand, and the second parameter is a variable or list of variables to be integrated over the real line.  To integrate a variable x over only [0,\[Infinity]) use {x,"Half"} in place of x.

Example 1: 
jIntegrateGaussian[x^3Exp[-3x^2+6*x*z-7y^2+x*y],{x,y}]

Example 2:
jIntegrateGaussian[x^3Exp[-3x^2+6*x*z-7y^2+x*y],{{x,"Half"},y}]

*)

jIntegrateGaussian[expression_,var_]:=jIntegrateGaussian[expression,var,False]; (*default is to print only errors, i.e. not verbose*)
jIntegrateGaussian[expression_,variable_,verbose_]:=Module[{explist,thing,lowpow,norm,isFull,var,integralType,ret,diracArgs,relDiracArg,relDiracArgCoeffList,diracReplacement,heavyEnforcement},
(*this function is the magic juice.  Finally does gaussian integrals fast and without a lot of complaining about variables being positive or whatever*)

ret[answer_]:=answer; (*this function is applied to all answers before they are returned.  For a while it was Simplify[], ExpToTrip[] (since TrigToExp is applied at the begining) or FullSimplify[], but this sometimes caused probalems with confusing the application of the formulas.  On the other hand, *not* using it also sometimes caused problems.  this choice would need some smart refinement*)

thing=Expand[TrigToExp[PowerExpand[expression]]]; (*This converts Sin[]'s and Cos[]'s to Exp[]'s so formulas can easily be applied.  This has the unfortunate side effect of exponentially increasing the size of expression as the number of variables in trig functions increases.  I don't know if there's a way to get around it though without adding a bunch of new forumlas*)

(*now we have to interpret "var" which may come in many forms.  they are, for type "Type",

 x
{x}
{x,y}
{x,y,z}
{x,"Type"}
{x,{y,"Type"}}
{{x,"Type"},y,z}
{{x,"Type"},{y,"Type2"},{z,"Type3"}}

and so on
*)
If[Head[variable]==List&&Length[variable]==0, (*blank list of variables, so just return with ret[] applied *)
Return[ret[thing]];
];

If[Head[variable]==List&&Length[variable]==1, (*is either {x} or {{x,"Type"}}.  Only one var.  So remove outer brackets. *)
Return[ret[jIntegrateGaussian[thing,variable[[1]],verbose]]];
];

If[Length[variable]>0&&(!StringQ[variable[[2]]]), (* Isn't either x or {x,"Type"}, and so it's an actual list of multiple variables (since {{x,"Type"}} was already excluded).  Thread over list *)
If[verbose,Print["Multiple integrals!"];];
thing=expression;
For[i=1,i<=Length[variable],i++,
thing=jIntegrateGaussian[thing,variable[[i]],verbose];
];
Return[ret[thing]];
];

(* We now can assume the variable list is a single variable.  If the expression to integrate is a sum, thread over the sum.  This is placed here before the form of the variable is fully pinned down in an old attempt to minimize the number of printed messages.  This might be improvable. *)
If[Head[thing]==Plus,
If[verbose,Print["Distributing over sum ",thing];];
(*Return[ret[Total[jIntegrateGaussian[#,variable,verbose]&/@thing]]];*)
(*Apparently "total" is unnecessary because the head of `thing' is already sum*)
Return[ret[jIntegrateGaussian[#,variable,verbose]&/@thing]];
];

(* We now can assume (1) it's not a sum and (2) there's only a single variable to integrate*)

(*set integral type*)
isFull=Null;
If[Length[variable]==0,(*default "Full" when none is specified*)
integralType="Full";isFull=True;
var=variable;(*Print["fullint"];*)
];
If[Head[variable]==List,(*Double checks whether the variable list is of the form 'x' or '{x,"Full"}' .*)
var=variable[[1]];
integralType=variable[[2]];
If[integralType=="Full",isFull=True,isFull=False]; 
];
If[!MemberQ[{"Full","Half","Theta","Phi"},integralType],
Print["Error! Can't interpret integral type!"];
Return[ret[thing]];
];

If[verbose,
Switch[integralType,
"Full",Print["Integrating ",thing," with respect to " ,var," over {-\[Infinity],\[Infinity]}..."],
"Half",Print["Integrating ",thing," with respect to " ,var," over {0,\[Infinity]}..."],
"Theta",Print["Integrating ",thing," with respect to " ,var," over {-1,1}..."],
"Phi",Print["Integrating ",thing," with respect to " ,var," over {0,2\[Pi]}..."]
(*"Phi" was previously "Theta", but I think this was a typo*)
]
];



explist=ExpCoefficientList[thing,var,verbose];
lowpow=Count[thing,var]+
If[Length[#]==0,0,thing[[#[[1]]]][[2]]]&[Position[thing,Power[var,_],1]]
;
(*this is the power of the variable outside the exponent*)
(*'Power[x,3]' is the actual explicit functional form of (x^3)*)
(*'Position[thing,Power[var,_]]' would return the list of (possibly multi-index) positions of all explicit powers of x*)
(*'Position[thing,Power[var,_],1]' returns the list of appearances of powers of x at the first level *)
(*We then apparently (?) assume there's only a single element of that list, and we use it to extract the powers from 'thing'*)

(*Now for the Dirac delta function handling *)
diracArgs=#[[1]]&/@Cases[thing,DiracDelta[_]];
(*This is a list of arguments for each Dirac delta function *)
If[Count[diracArgs,DiracDelta[_],Infinity]>0,
Print["Error: Nested Dirac delta functions containing ", var];
Return[ret[thing]];
];
If[Length[diracArgs]>0,
If[verbose,Print["Dirac delta function(s) found."]];
relDiracArg=diracArgs[[First/@Position[diracArgs,var]]];
(*Above is a list of arguments for those Dirac delta functions that contain 'var'*)
If[Count[#[[1]]&/@Cases[thing/Product[DiracDelta[relDiracArg[[r]]],{r,1,Length[relDiracArg]}],DiracDelta[_]],var,Infinity]>0,
(*the first argument of Count is the Dirac arguments remaining*)
Print["Error: I can't understand these Dirac delta functions. ",var," is still left over in one."];
Return[ret[thing]];
];
If[Length[relDiracArg]==0,
If[verbose,Print["None of the Dirac delta functions are relevant."]]
];
If[Length[relDiracArg]>0,
  If[verbose,Print["We have found a relevant Dirac delta function."]];
  relDiracArgCoeffList=CoefficientList[relDiracArg[[1]],var];
  If[Count[relDiracArgCoeffList[[1]],var]!=0||Length[relDiracArgCoeffList]>2,
  Print["Error: Diract delta functions contains terms non-linear in ",var];
  Return[ret[thing]];
  ];
  diracReplacement=((relDiracArg[[1]]-relDiracArgCoeffList[[2]]*var)/-relDiracArgCoeffList[[2]]);
  heavyEnforcement=Switch[integralType,
  "Full",1,
  "Half", HeavisideTheta[diracReplacement],
  "Phi", HeavisideTheta[diracReplacement]* HeavisideTheta[2\[Pi]-diracReplacement],
  "Theta", HeavisideTheta[diracReplacement-1]* HeavisideTheta[1-diracReplacement]
  ];
  
  If[verbose,Print["Integrating DiracDelta[",relDiracArg[[1]],"] over ",var," results in each instance of ",var," being replaced with ", diracReplacement," and an overall multiplicative factor of ",1/Abs[relDiracArgCoeffList[[2]]]]];

  Return[ret[(thing*heavyEnforcement)/DiracDelta[relDiracArg[[1]]]*1/Abs[relDiracArgCoeffList[[2]]]/.{var-> diracReplacement}]];
];
(*_,Print["Error: More than one Dirac delta function contains the variable ", var];Return[ret[thing]];*)(* We are now allowing multiple Dirac delta functions to contain the same var, and only using the first one *)

];

norm=Simplify[(thing/var^lowpow)/Exp[Sum[var^j*explist[[j+1]],{j,1,Length[explist]-1}]]];
If[Length[explist]==0,(*no exponent*)
explist={0};
];
explist[[1]]=0;(*this const term in the exponent appears in the norm, so it shouldn't be double counted*)

If[Count[norm,var,Infinity]!=0,(*this checks that we have broken the expression correctly into parts.  If "norm" has the variable in it, we're screwed and have to figure out what went wrong*) 
Print["Error: I do not understand this form"];
Print["expression: ", expression];
Print["thing: ", thing];
Print["explist: ", explist];
Print["lowpow: ", lowpow];
Print["lowpow search: ",Cases[thing,Power[var,_]]];
Print["norm: ", norm];
Print["remaining var positions: ",Position[norm,var][[1]]];
Print["isFull: ", isFull];
Return[ret[thing]];
];



If[verbose,Print["Applying Formula."];];(*apply formula based on integral type*)
Return[ret[Switch[integralType,
"Full",jGaussianFormula[thing,explist,lowpow,norm,var,isFull,verbose],
"Half",jGaussianFormula[thing,explist,lowpow,norm,var,isFull,verbose],
"Theta",jThetaFormula[thing,explist,lowpow,norm,var,verbose],
"Phi",jPhiFormula[thing,explist,lowpow,norm,var,verbose]
]]];

]; (*End jIntegrateGaussian  *)


jThetaFormula[thing_,explist_,lowpow_,norm_,var_,verbose_]:=Module[{},
(* a "theta"-type integral is just integrating a variable like Cv1v2 (cos between vectors v1 and v2) from -1 to 1.  It expects to integrate with respect to Cos[theta] for some theta, and for the only appearance of Cv1v2 to be linear in the exponent*)

If[lowpow!=0,
Print["Error! Variable ",var," has powers outside the exponent in theta-type integral!"];
Return[thing];
];
If[Length[explist]>2,
Print["Error! Variable ",var," has powers higher than linear in exponent in theta-type integral!"];
Return[thing];
];
If[Length[explist]==1,
If[verbose,Print["Warning: variable ",var," does not appear in theta-type integral!"]];
Return[2*norm];
];
If[Length[explist]==2,
Return[norm*(Exp[explist[[2]]]-Exp[-explist[[2]]])/explist[[2]](*2*Sin[explist\[LeftDoubleBracket]2\[RightDoubleBracket]/\[ImaginaryI]]/(explist\[LeftDoubleBracket]2\[RightDoubleBracket]/\[ImaginaryI])*)];
];
];
jPhiFormula[thing_,explist_,lowpow_,norm_,var_,verbose_]:=Module[{},
(* this definition is completely incomplete *)
(* a "phi"-type integral is just integrating a variable like phi from 0 to 2\[Pi].  Right now, it expect no appearance of phi and so just multiplies by 2\[Pi]*)
If[Simplify[norm==thing],Return[thing*2\[Pi]];];

Print["Error! Can't do this non-trivial phi-type integral yet!"];
Print["norm: ",norm];
Print["thing: ",thing];
Return[thing];
];


jGaussianFormula[thing_,explist_,lowpow_,norm_,var_,isFull_,verbose_]:=Module[{n,T(*\[Equal]posR*),S},
(* This is the normal gaussian integral, either from -\[Infinity] to \[Infinity] ("Full", or no specification) or from 0 to \[Infinity] ("Half").  The latter is typically for the radius r=|x| in radial coordinates.  It expects the x^2 term in the differential of radial coordinates to already be included.  It does *not* add it.  *)


If[verbose,Print["Checking for linear dependence."]];
(*First, let's see if this exponent is all linear and we need evaluate this as a Dirac delta function*)
If[Length[explist]==2,
If[isFull,(*OK, the exponent has only a linear term, so we need to do Dirac or throw an error*)
If[NumberQ[explist[[2]]]&&Re[explist]!=0,
Print["Error: Integral over ",var," does not converge. The exponent has no quadratic term, and the linear term has a real part"];
Return[thing];
];
If[lowpow!=0,
Print["Error: Variable ",var," no quadratic term in exponent and a non-zero power outside the exponent"];
Return[thing];
];
(*OK, we're good to make the replacement*)
If[verbose,
Print["Replacing integral over ",var," with Dirac delta function ",2\[Pi]*DiracDelta[Simplify[explist[[2]]/I]]];
If[!NumberQ[explist[[2]]],Print["Note that I can't guarantee that the coefficient is pure imaginary."];];
];
Return[norm*2\[Pi]*DiracDelta[Simplify[explist[[2]]/I]]];
,
(*not Full {\[Infinity],\[Infinity]} here*)
Print[Row[{"Error! Variable ",var," has linear term but no quadratic term in exponent, and we're not doing a full integral over {-\[Infinity],\[Infinity]} so I can't insert a Dirac delta function!"}]];
Return[thing];
]
,
If [verbose,Print["We're not linear."]];
];

Switch[Length[explist],
2,
Print["Error: We shouldn't be here!  The linear term was handled above."];          
,
1,Print[Row[{"Error! Variable ",var," does not appear in exponent!"}]];Return[thing];,
0,Print[Row[{"Error! Variable ",var," does not appear in exponent! And exponent is empty!"}]];Return[thing];,
_?(#>3&),Print[Row[{"Error! Variable ",var," appears as cubic or higher power in exponent!"}]];Return[thing];
];
If[Length[explist]!=3,Print["Error! How did I get here?"];Return[thing];];

T=-explist[[3]];
S=explist[[2]];
n=lowpow;


If[NumberQ[T]&&T<0,
Print[Row[{"Error! Second power of ",var," has positive coefficent!"}]];Return[thing];
];

If[n==0, (*no powers of variable outside exponent, so answer is simple*)
If[verbose,Print["Simple Gaussian formula used"]];
Return[Simplify[norm*If[isFull,
 Sqrt[\[Pi]]/Sqrt[T] E^(S^2/(4 T)),
(E^(S^2/(4 T)) Sqrt[\[Pi]] (1+Erf[S/(2 Sqrt[T])]))/(2 Sqrt[T])
]]];
];

If[!NumberQ[n]&&verbose,Print["Note that I can't check that the power ",n," outside the main exponent is > -1 because it's not an explicit number.  I'll assume it is."];
];

If[(n>-1||!NumberQ[n])&&S==0, (*kinda simple case *)
If[verbose,Print["Kinda simple Gaussian formula used"];];
Return[Simplify[norm*If[isFull,
1/2 (1+(-1)^n) T^(-(1/2)-n/2) Gamma[(1+n)/2],
1/2 T^(-(1/2)-n/2) Gamma[(1+n)/2]
]]];
];
If[n>-1||!NumberQ[n], (*complicated case *)
If[verbose,
Print["Complicated Gaussian formula used"];
If[!NumberQ[n],Print["Note that I can't check that the power ",n," outside the main exponent is > -1 because it's not an explicit number"];
];
];
Return[Simplify[norm*If[isFull,
1/2 T^(-1-n/2) (-(-1+(-1)^n) S Gamma[1+n/2] Hypergeometric1F1[1+n/2,3/2,S^2/(4 T)]+(1+(-1)^n) Sqrt[T] Gamma[(1+n)/2] Hypergeometric1F1[(1+n)/2,1/2,S^2/(4 T)]),
1/2 T^(-1-n/2) (S Gamma[1+n/2] Hypergeometric1F1[1+n/2,3/2,S^2/(4 T)]+Sqrt[T] Gamma[(1+n)/2] Hypergeometric1F1[(1+n)/2,1/2,S^2/(4 T)])
]]];
];

If[n<=-1,
Print["Error: I don't know how to do a Gaussian integral with powers outside the exp below -1"];Return[thing];
];

Print["Error: I don't know how I got here."];
];

Clear[ExpCoefficientList];
ExpCoefficientList[expression_,var_,verbose_]:=Module[{thing,exps,powers,bases,pow,doit},
(*this function takes an expression like x^2Exp[z^7]y^4(\[Pi]/2)Exp[-x^2-4x+q] and extracts the coefficient list for powers of a particular variable var for the part *inside* the exponent. All exponents are taken together as one. *)

doit[ex_,vr_]:=Module[{},
exps=Position[ex,Power[_?((NumberQ[#]||#==E)&),_],1];
bases=((ex[[#]])[[1]]&)/@exps;
powers=((ex[[#]])[[2]]&)/@exps;
If[Count[bases,NumberQ]>0,(*E, i.e. \[ExponentialE], is *not* considered a number*)
If[verbose,Print["Warning: Correcting bases"]]; powers=Table[powers[[j]]*Log[bases[[j]]],{j,1,Length[bases]}];
];
pow=Total[powers];
Return[(If[Length[#]==1&&verbose,Print["Warning: trivial coefficient list"];];#)&[CoefficientList[pow,vr]]];
];

thing=PowerExpand[expression];
If[Head[thing]==Plus,
Print["Error!  Why plus here?"];
];

If[Head[thing]==Times,
Return[doit[thing,var]]];

If[Head[thing]==Power,
(*If[!NumberQ[thing\[LeftDoubleBracket]1\[RightDoubleBracket]]&&thing\[LeftDoubleBracket]1\[RightDoubleBracket]\[NotEqual]E,(*no exponents*)
If[verbose,Print["Trivial coefficient list"]];Return[CoefficientList[0,var]];,
If[thing\[LeftDoubleBracket]1\[RightDoubleBracket]\[NotEqual]E,(*weird base*)
Print["Warning: Correcting bases"];Return[CoefficientList[thing\[LeftDoubleBracket]2\[RightDoubleBracket]*Log[thing\[LeftDoubleBracket]1\[RightDoubleBracket]],var]];,
Return[CoefficientList[thing\[LeftDoubleBracket]2\[RightDoubleBracket],var]];
];
];*)
Return[ExpCoefficientList[thing*"dummy",var,verbose]]; (*to force threading over the Times operator in rare case that expression contains only one multiplicative term*)
];
];