(* ::Package:: *)

BeginPackage["Landen`"];


Print["Landen integration package by Dante Manna, Luis Medina, Victor H. Moll, and Armin Straub
	accompanying the paper \"A new numerical method for the integration of rational functions\"
	-- Tulane University -- Version 1.0.0 (2009/02/13)"];


NLandenIntegrate::usage=
	"NLandenIntegrate[f] numerically computes the definite integral of the rational function f using Landen transformations.
	The following options are available:
		MethodOrder - order of the Landen transformations (defaults to 2)
		WorkingPrecision - the precision used for the Landen iterations (defaults to MachinePrecision)
		PrecisionGoal - the precision that the method tries to achieve for the integral (defaults to half the WorkingPrecision)
		MaxIterations - maximal number of Landen iterations to be performed before giving up (defaults to $RecursionLimit)
		StepMonitor - a function that is called after each Landen iteration with the current rational function as a parameter (as a list of coefficients suitable for RatFromCoeffs)
	Usage:
		NLandenIntegrate[1/(1+2x+3x^2)]
		NLandenIntegrate[1/(1+2x+3x^2),PrecisionGoal\[Rule]40]
		NLandenIntegrate[1/(1+2x+3x^2),MethodOrder\[Rule]3,StepMonitor\[Rule](Print[RatFromCoeffs[#,x]]&),MaxIterations\[Rule]3,WorkingPrecision\[Rule]Infinity]";
GenerateLandenTransforms::usage=
	"GenerateLandenTransforms[p] generates the Landen transformations of degree up to p.
	The following options are available:
		MethodOrder - order of the Landen transformations (defaults to 2)
		FunctionName - the Landen transformations will be available under this name (defaults to LandenStep)
		FileName - If specified the generated Landen transformations (including existing ones of the same name) will be saved to this file
	Usage:
		GenerateLandenTransforms[12]
		GenerateLandenTransforms[8,MethodOrder->3]
		GenerateLandenTransforms[FileName\[Rule]\"LandenTransforms.m\"]";
GenerateLandenTransform::usage=
	"GenerateLandenTransform[p,m] generates a Landen transformation for rational functions B/A where
		p is the degree of A (B is of degree p-2), and
		m is the order of the transformation.
	Usage:
		f=GenerateLandenTransform[3,2];
		f[{{b0,b1},{a0,a1,a2,a3}}]";
RatFromCoeffs::usage=
	"RatFromCoeffs[{Ln,Ld},x] returns the rational function f=A/B in x such that
		Ln are the coefficients of A, and
		Ld are the coefficients of B
	each list starting with the highest order coefficient.
	Usage:
		RatFromCoeffs[{{1},{3,2,1}},x]";
CoeffsFromRat::usage=
	"CoeffsFromRat[f] takes a rational function f=A/B and returns a list {Ln, Ld} where
		Ln are the coefficients of A, and
		Ld are the coefficients of B
	each list starting with the highest order coefficient.
	The following options are available:
		Variable - to specify which variable f is to be considered a rational function in (defaults to Automatic)
		LandenFormat - if set to True the lists Ln,Ld will be left-padded with zeroes so that Length[Ln]+2==Length[Ld] (defaults to False)
	Usage:
		CoeffsFromRat[1/(1+2x+3x^2)]
		CoeffsFromRat[1/(1+a x+3x^2),Variable\[Rule]x]
		CoeffsFromRat[1/(1+2x+3x^2+x^3),LandenFormat\[Rule]True]";


Begin["Private`"];


(* ::Subsection:: *)
(*Landen Integration*)


Options[NLandenIntegrate]={
	WorkingPrecision->Automatic,
	PrecisionGoal->Automatic,
	MaxIterations:>$RecursionLimit,
	MethodOrder->2,
	StepMonitor->(None&)};
NLandenIntegrate::largenumerator="The integral does not exist because the numerator is too large.";
NLandenIntegrate::lowprec="The precision dropped below the precision goal of `1`.";
NLandenIntegrate::noconvergence="No convergence observed after `1` iterations.";
NLandenIntegrate::infiniteprecision="You may want to set PrecisionGoal to something less than Infinity or decrease MaxIterations.";
NLandenIntegrate::nolandentransform="The Landen transformation of degree `1` and order `2` is not available. You may use GenerateLandenTransforms[`1`,MethodOrder->`2`] to generate it.";
NLandenIntegrate[rat_,OptionsPattern[]]:=Module[{p,B,A,B1,A1,AC,c,c1,T,prec,precg},
	(* Set up precision parameters. *)
	prec=OptionValue[WorkingPrecision];
	precg=OptionValue[PrecisionGoal];
	If[precg===Automatic,
		If[prec===Automatic,prec=MachinePrecision];precg=prec/2,
		If[prec===Automatic,prec=2precg]];
	If[precg===Infinity&&OptionValue[MaxIterations]===$RecursionLimit,
		Message[NLandenIntegrate::infiniteprecision]];
	(* Convert rational function to list of coefficients. *)
	If[MatchQ[rat,{_,_}],
		{B,A}=rat,
		{B,A}=CoeffsFromRat[rat];
		If[Head[A]===Symbol,Return[$Failed]]];
	(* Numerator has to have degree at least two less. *)
	If[Length[B]+2>Length[A],
		Message[NLandenIntegrate::largenumerator];
		Return[Indeterminate]];
	B=Join[Table[0,{Length[A]-Length[B]-2}],B];
	(* Set up the candidates that the denominator may converge to. *)
	p=Length[A]-1;
	If[EvenQ[p]&&FreeQ[A,_Complex],
		AC=Module[{x},{Reverse[CoefficientList[(1+x^2)^(p/2),x]]}],
		AC=Module[{x},Table[Reverse[CoefficientList[(x-I)^q (x+I)^(p-q),x]],{q,0,p}]]];
	(* Use working precision. *)
	{B,A}=SetPrecision[{B,A},prec];
	{B1,A1,c1}={B,A,1};
	(* Pick the Landen transformation to use. *)
	T=Global`LandenStep[#,OptionValue[MethodOrder]]&;
	Catch[Do[
		{B,A,c}={B1,A1,c1};
		With[{BA1=T[{B,A}]},
			If[!ListQ[BA1],
				Message[NLandenIntegrate::nolandentransform,p,OptionValue[MethodOrder]];
				Throw[$Failed]];
			{B1,A1}=BA1];
		(* Normalize so that the denominator has leading coefficient 1. *)
		If[A1[[1]]!=0,
			{B1,A1}={B1,A1}/A1[[1]];
			c1=B1[[1]]/A1[[1]],
			c1=Indeterminate];
		(* Tell the outside world about our progress. *)
		OptionValue[StepMonitor][{B1,A1}];
		(* Celebrate if error estimates drop below precision goal. *)
		If[RelError[c1,c]<=10^(-precg)&&Min[Table[RelErrorList[A1,C],{C,AC}]]<=10^(-precg),
			Throw[c1 Pi]];
		(* Or go home if current precision is hopeless. *)
		If[Precision[A1[[1]]]<precg&&precg!=Infinity,
			Message[NLandenIntegrate::lowprec,precg];
			Throw[Indeterminate]];,
		{OptionValue[MaxIterations]}];
	(* Oh no. No convergence observed. *)
	Message[NLandenIntegrate::noconvergence,OptionValue[MaxIterations]];
	Indeterminate]
];


(* ::Subsection:: *)
(*Helper Functions*)


RelError[xn_,x_]:=Abs[(xn-x)/If[x!=0,x,1]]
RelErrorList[xn_,x_]:=Max[Table[RelError[xn[[k]],x[[k]]],{k,1,Length[x]}]]


PolyFromCoeffs[coeffs_List,x_]:=
	Total[coeffs Reverse[x^Range[0,Length[coeffs]-1]]];
RatFromCoeffs[{numcoeffs_List,dencoeffs_List},x_]:=
	PolyFromCoeffs[numcoeffs,x]/PolyFromCoeffs[dencoeffs,x];


Options[CoeffsFromRat]={
	Variable->Automatic,
	LandenFormat->False};
CoeffsFromRat::notrat="`1` is not a rational function in one variable.";
CoeffsFromRat[f_,OptionsPattern[]]:=Module[{vars,x,n,d,Ln,Ld,p},
	(* Determine variable. *)
	If[OptionValue[Variable]===Automatic,
		vars=Variables[f];
		If[Length[vars]==1&&Head[First[vars]]===Symbol,
			x=First[vars],
			Message[CoeffsFromRat::notrat,f];
			Return[$Failed]],
		x=OptionValue[Variable]];
	(* Get coefficients. *)
	n=Numerator[f];
	d=Denominator[f];
	If[PolynomialQ[n,x]&&PolynomialQ[d,x],
		{Ln,Ld}=Reverse/@{CoefficientList[n,x],CoefficientList[d,x]};
		If[OptionValue[LandenFormat]===True,
			p=Max[Length[Ld],Length[Ln]+2];
			Ln=Join[Table[0,{p-2-Length[Ln]}],Ln];
			Ld=Join[Table[0,{p-Length[Ld]}],Ld];
		];
		Return[{Ln,Ld}];
	];
	Message[CoeffsFromRat::notrat,f];
	$Failed
];


(* ::Subsection:: *)
(*Landen Transform Generator*)


GenerateLandenTransform[p_,m_]:=Module[{P,Q,A,B,A1,B1,F,E1,Z,C1,T,M1,M2,z,x,a,b,
	s=m*p-2,\[Nu]=Floor[p/2],\[Lambda]=Floor[(m*p-2)/2]},
	(* R=P/Q such that R (cot (\[Phi]))==cot (m \[Phi])*)
	P=Sum[(-1)^i*Binomial[m,2*i]*z^(m-2*i),{i,0,Floor[m/2]}];
	Q=Sum[(-1)^i*Binomial[m,2*i+1]*z^(m-1-2*i),{i,0,Floor[(m-1)/2]}];
	(* B/A is the rational function to be Landen transformed *)
	A=Sum[a[i] z^(p-i),{i,0,p}];
	B=Sum[b[i] z^(p-2-i),{i,0,p-2}];
	(* Denominator of the Landen transform *)
	A1=Resultant[PolynomialRemainder[A,P-x*Q,z],P-x*Q,z,Method->"SylvesterMatrix"];
	E1=Sum[Coefficient[A1,x,p-i]*(P)^(p-i)*(Q)^i,{i,0,p}];
	C1=Table[Coefficient[B*PolynomialQuotient[E1,A,z],z,s-i],{i,0,s}];
	T[x_,a_,b_]:=Sum[(-1)^(a-x+j)*Binomial[a,x-j]*Binomial[b,j],{j,0,x}];
	M1[j_,\[Alpha]_,\[Beta]_,\[Gamma]_]:=(-1)^(j+\[Alpha]-\[Beta])*C1[[2*j+1]]*2^(2*(\[Alpha]-\[Beta]))*\[Alpha]*(1/(2*\[Alpha]-\[Beta]))*Binomial[2*\[Alpha]-\[Beta],\[Beta]]
		*Binomial[\[Nu]-\[Alpha]-1+\[Beta],\[Gamma]]*(T[\[Lambda]+\[Alpha]*m,2*j,s-2*j]+T[\[Lambda]-\[Alpha]*m,2*j,s-2*j]);
	M2[j_,\[Alpha]_,\[Beta]_,\[Gamma]_]:=(-1)^(j+\[Beta])*C1[[2*j+2]]*2^(2*\[Beta]+1)*Binomial[\[Alpha]+\[Beta],2*\[Beta]+1]*Binomial[\[Nu]-2-\[Beta],\[Gamma]]
		*(T[\[Lambda]+\[Alpha]*m,2*j+1,s-2*j-1]-T[\[Lambda]-\[Alpha]*m,2*j+1,s-2*j-1]);
	(* Numerator of the Landen transform *)
	B1=1/2^s(Sum[Binomial[\[Nu]-1,\[Gamma]]Sum[(-1)^j C1[[2 j+1]] T[\[Lambda],2 j,s-2 j],{j,0,\[Lambda]}]x^(2\[Gamma]),{\[Gamma],0,\[Nu]-1}]
		+Sum[Sum[Sum[Sum[M1[j,\[Alpha],\[Beta],\[Gamma]],{\[Beta],0,\[Alpha]}],{\[Alpha],1,\[Nu]-1-\[Gamma]}],{j,0,\[Lambda]}]x^(2\[Gamma]),{\[Gamma],0,\[Nu]-2}]
		+Sum[Sum[Sum[Sum[M1[j,\[Alpha],\[Beta],\[Gamma]],{\[Beta],\[Alpha]-\[Nu]+\[Gamma]+1,\[Alpha]}],{\[Alpha],\[Nu]-\[Gamma],\[Nu]-1}],{j,0,\[Lambda]}]x^(2\[Gamma]),{\[Gamma],1,\[Nu]-1}]
		+Sum[Sum[Sum[Sum[M2[j,\[Alpha],\[Beta],\[Gamma]],{\[Beta],0,\[Alpha]-1}],{\[Alpha],1,\[Nu]-1-\[Gamma]}],{j,0,\[Lambda]-1}]x^(2\[Gamma]+1),{\[Gamma],0,\[Nu]-2}]
		+Sum[Sum[Sum[Sum[M2[j,\[Alpha],\[Beta],\[Gamma]],{\[Beta],0,\[Alpha]-1}],{\[Alpha],\[Nu]-\[Gamma],\[Nu]-1}],{j,0,\[Lambda]-1}]x^(2\[Gamma]+1),{\[Gamma],1,\[Nu]-2}]);
	(* Extract coefficients of the Landen transform F=J/H *)
	F={Table[Coefficient[B1,x,i],{i,p-2,0,-1}],Table[Coefficient[A1,x,i],{i,p,0,-1}]};
	F=Expand[F];
	Function[F/.{a[n_]:>#[[2,n+1]],b[n_]:>#[[1,n+1]]}]
];


DefineLandenTransform[p_Integer,m_Integer,name_,T_]:=Module[{DefineLanden,varnames,varpatterns},
	DefineLanden[var_,val_]:=(name[var,m]:=val);
	varnames={Table[ToExpression["b"<>ToString[i]],{i,0,p-2}],
		Table[ToExpression["a"<>ToString[i]],{i,0,p}]};
	varpatterns={Table[ToExpression["b"<>ToString[i]<>"_"],{i,0,p-2}],
		Table[ToExpression["a"<>ToString[i]<>"_"],{i,0,p}]};
	DefineLanden[varpatterns,T[varnames]]];


Options[GenerateLandenTransforms]={
	FileName->None,
	FunctionName->Global`LandenStep,
	MethodOrder->2};
GenerateLandenTransforms[p0_Integer:0,OptionsPattern[]]:=Module[{p=p0,m,file,name,T,T0},
	m=OptionValue[MethodOrder];
	file=OptionValue[FileName];
	name=OptionValue[FunctionName];
	(* Generate the Landen transform of degree p and order m,
		and use it to establish the Landen transforms of degree less than p, too. *)
	If[p>1,
		If[OddQ[p],p++];
		T=GenerateLandenTransform[p,m];
		Do[
			T0=Function[ba,Drop[#,p-q]&/@T[Join[Table[0,{p-q}],#]&/@ba]];
			DefineLandenTransform[q,m,name,T0],
			{q,2,p-1}];
		DefineLandenTransform[p,m,name,T]];
	If[file=!=None,
		Put[file];
		Save[file,Evaluate[name]]];
];


End[];


EndPackage[];