(* :Title: Gaussian quadrature *)

(* :Author: Mark Fisher *)

(* :Context: GaussianQuadrature` *)

(* :Package Version:
	1.0 (February 1999)
	1.1 (April 1999) Use $MinPrecision to get precision up to WorkingPrecision
	1.2 (June 1999)  Apply N after computing exact solutions for abscissas
						and weights
*)

(* :Mathematica Version: 3.0 *)

(* :Summary:
This package provides two functions: GaussianQuadraturePairs and
GaussianProductRule.

GaussianQuadraturePairs produces a list of n abscissa/weight pairs {x,
w} for n-point Gaussian quadrature. Four types of orthogonal
polynomials are supported: HermiteH, LaguerreL, LegendreP, and
ChebyshevT. (The Standard Package NumericalMath`GaussianQuadrature`
computes abscissa/weight pairs only for the LegendreP case.)
WorkingPrecision can be specified (including Infinity). The option
ScalingParameters can be used to rescale the interval and/or weight
function, depending on the type.

In order to guarantee that all digits are correct up to the
WorkingPrecision specified, the abscissas and weights are first
computed exactly and then computed numerically with a precision of
Max[WorkingPrecision, $MachinePrecision + 1]. If WorkingPrecision is
$MachinePrecision or less, N is applied to the result to reduce the
precision to $MachinePrecision.

GaussianProductRule takes a sequence of pairs (produced by
GaussianQuadraturePairs) and returns a tensor array of abscissa/weight
tuples that can be used for multidimensional Gaussian quadrature.
*)

(* :Sources:
Kenneth L. Judd (1998), Numerical Methods in Economics, MIT Press, chapter 7.
Abramowitz and Stegun (1964), Handbook of Mathematical Tables, chapter 25.
*)

(* :Examples:
p1 = GaussianQuadraturePairs[10, Type -> HermiteH,
	ScalingParameters -> {0, 1}]
p2 = GaussianQuadraturePairs[7, Type -> LegendreP,
	ScalingParameters -> {-3, 10}]
prod = GaussianProductRule[p1, p2]
Needs["Graphics`Graphics3D`"]
ListSurfacePlot3D[prod, Axes -> True, BoxRatios -> {1, 1, .5},
	PlotRange -> All]
*)

BeginPackage["GaussianQuadrature`"]

GaussianQuadrature::usage = "GaussianQuadrature is a package for doing
Gaussian quadrature."

GaussianQuadraturePairs::usage = "GaussianQuadraturePairs[n] returns n
pairs of abscissas and weights suitable for n-point Gaussian quadrature.
The default settings for the options are Type -> HermiteH, WorkingPrecision
-> $MachinePrecision, and ScalingParameters -> Automatic. Other valid types
are LaguerreL, LegendreP, and ChebyshevT. WorkingPrecision can be set to
\[Infinity]. Setting WorkingPrecsion less than $MachinePrecision has no
effect."

Type::usage = "Type is an option that specifies which type of
orthogonal polynomials to use for Gaussian quadrature. Valid types are
HermiteH (the default), LaguerreL, LegendreP, and ChebyshevT."

ScalingParameters::usage = "ScalingParameters is an options that
specifies how to rescale the abscissas and the weights. The default is
ScalingParameters -> Automatic, which does no rescaling. To affect
rescaling, the proper usage is ScalingParameters -> {a, b}. The
interpretation of the scaling parameters depends on the type of
orthogonal polynomials used. If the type is either LegendreP or
ChebyshevT, then a and b specify the minimum and maximum of a finite
interval (a < b). If the type is HermiteH, then a and b specify the
mean and the standard deviation of a Gaussian distribution. If the
type is LaguerreL, then a and b specify the coefficient of x in the
exponential weighting function and the lower bound of the integral, as
in \!\(\[Integral]\_b\%\[Infinity]\( \[ExponentialE]\^\(\(-a\)\ x\)\)
f[x] \[DifferentialD]x\), where a \[NotEqual] 0."

GaussianProductRule::usage = "GaussianProductRule[pair1, pair2, ...]
takes a sequence of lists of abscissa/weight pairs (produced by
GaussianQuadraturePairs) and returns a tensor array of abscissa/weight
tuples that can be used for multidimensional Gaussian quadrature. By
setting the option Flatten -> True, GaussianProductRule will return a
list of tuples instead."

Begin["`Private`"]

Options[GaussianQuadraturePairs] = {
	Type -> HermiteH,
	WorkingPrecision :> $MachinePrecision,
	ScalingParameters -> Automatic}

GaussianQuadrature::badtype = "Unknown type: ``."

GaussianQuadrature::badprec = "WorkingPrecision should be either
a positive integer or \[Infinity]."

GaussianQuadrature::badscale = "ScalingParameters should be either
Automatic or a pair of numbers."

GaussianQuadraturePairs[n_Integer?Positive, opts___?OptionQ] :=
	Module[{type, wp, sp, a, b, A, B, C, wfun, x, abscissas, weights},
	{type, wp, sp} = {Type, WorkingPrecision, ScalingParameters} /.
		{opts} /. Options[GaussianQuadraturePairs];
	Switch[wp, _Integer?Positive | Infinity, {},
			_, Message[GaussianQuadrature::badprec]; Return[$Failed]];
	Switch[sp,
		Automatic, {A, B, C} = {1, 0, 1},
		{_, _},
			{a, b} = SetPrecision[sp, Infinity];
			Switch[type,
				HermiteH, {A, B, C} = {Sqrt[2] b, a, 1/Sqrt[Pi]},
				LegendreP | ChebyshevT, {A, B, C} = {0, a, 0} + (b-a)/2,
				LaguerreL, {A, B, C} = {1/a, b, Exp[-a b]/a},
				_,	Message[GaussianQuadrature::badtype, type];
					Return[$Failed]],
		_, Message[GaussianQuadrature::badscale]; Return[$Failed]];
	wfun = Function @ Evaluate @
		Switch[type,
			HermiteH,   2^(n-1) n! Sqrt[\[Pi]]/(n HermiteH[n - 1, #])^2,
			LaguerreL,  #/((n+1)^2 LaguerreL[n + 1, #]^2),
			LegendreP,  -2/((n+1) LegendreP[n + 1, #] D[LegendreP[n, #], #]),
			ChebyshevT, Pi/n
		];
	abscissas = x /. Solve[type[n, x] == 0];
	weights = wfun /@ abscissas;
	{abscissas, weights} =
		N[{abscissas, weights}, Max[wp, $MachinePrecision + 1]];
	If[wp < Infinity,
		{abscissas, weights} = {abscissas, weights} /. Complex[x_, _] -> x];
	If[wp <= $MachinePrecision,
		{abscissas, weights} = N[{abscissas, weights}]];
	Sort @ Transpose[{A abscissas + B, C weights}]
	]

Options[GaussianProductRule] = {Flatten -> False}

GaussianProductRule[pairs:{{_, _} ..} .., opts___?OptionQ] :=
	Module[{len, weights, abscissas, tuples},
	len = Length[{pairs}] - 1;
	(
	weights   = Outer[Times, ##]& @@ Map[Last,  {pairs}, {2}];
	abscissas = Outer[List,  ##]& @@ Map[First, {pairs}, {2}];
	tuples = MapThread[Append, {abscissas, weights}, Length[{pairs}]];
	If[TrueQ[Flatten /. {opts} /. Options[TensorFromPairs]],
		Flatten[tuples, len],
		tuples]
	) /; len > 0
	]

End[]
EndPackage[]