(*:Name: NumericalMath`CauchyPrincipalValue` *)

(*:Title: Numerical Evaluation of Cauchy Principal Values *)

(*:Author: Jerry B. Keiper *)

(*:Summary:
This package introces a function for numerically evaluating Cauchy
principal values of integrals that may not be Riemann integrable.
The built-in function NIntegrate is used in a symmetric way to get
cancelation.
*)

(*:Context: NumericalMath`CauchyPrincipalValue` *)

(*:Package Version: Mathematica 2.0 *)

(* :Copyright: Copyright 1990-1999,  Wolfram Research, Inc.
*)

(* :History:
	Version 1.2 by Jerry B. Keiper, January 1990.
	Revised by Jerry B. Keiper, November 1990.
*)

(*:Keywords: Cauchy principal value, integration, divergent integral *)

(*:Source:
	E. T. Whittaker & G. N. Watson, A Course of Modern Analysis,
	Cambridge University Press, (4th ed.), Section 4.5.
*)

(*:Mathematica Version: 2.0 *)

(*:Limitation:
	Since CauchyPrincipalValue is a numerical function and it
	relies on symmetry near the singular points, loss of digits
	due to cancelation can be a problem.  The user must know
	and specify EXACTLY where the singularities are, or the
	result can be grossly in error.  Again this is because of
	the extreme sensitivity to perturbations.  
*)

BeginPackage["NumericalMath`CauchyPrincipalValue`"]

CauchyPrincipalValue::usage = 
"CauchyPrincipalValue[f[x], {x, a, {b, eps}, c}, options] gives the Cauchy
principal value of NIntegrate[f[x], {x, a, c}, options] with a singularity
at x == b.  This is accomplished as NIntegrate[f[x], {x, a, b-eps, b, b+eps, c},
options], where the portion NIntegrate[f[x], {x, b-eps, b, b+eps}, options]
is evaluated as NIntegrate[f[b+t] + f[b-t], {t, 0, eps}, options].  If eps
does not appear, a default value will be used.  More generally,
CauchyPrincipalValue[f[x], {x, x0, ..., xi, ..., {xj}, ..., xn}] gives
the integral along the path {x0, ..., xi, ..., xj, ..., xn} using the
Cauchy principle value at those abscissas enclosed in {}, any of which
may include a user-specified eps."

Begin["NumericalMath`CauchyPrincipalValue`Private`"]

CauchyPrincipalValue[f_, range_List, options___] :=
	Module[{e}, e /; (e = CPV0[f, range, options]) =!= Fail]

CPV0[f_, range_, options___] :=
    Module[{a, b, eps, epsa, epsc, i, ilim, sum=0, temp, wprec},
	ilim = Length[range]-1;
	wprec = WorkingPrecision /. {options} /. Options[NIntegrate];
	For[ i = 2, i <= ilim, i++,
	    If[ i == 2, 
		a = range[[2]];
		If[ Head[a] === List, Return[Fail]],
	      (* else *)
		a = c
	    ];
	    b = range[[i+1]];
	    If[ Head[b] =!= List, 
		c = b;
		If[ a =!= c,
		    temp = NIntegrate[f, Evaluate[{range[[1]],a,c}], options];
		    If[ !NumberQ[temp], Return[Fail]];
		    sum += temp;
		];
		Continue[]  (* next i in For[ ] loop *)
	    ];
	    c = range[[i+2]];
	    If[ Head[c] === List,
		If[ i == ilim, Return[Fail]];
		c = (c[[1]] + b[[1]])/2
	    ];
	    Which[
		Length[b] > 2,
		    Return[Fail],
		Length[b] == 2,
		    eps = b[[2]];
		    b = b[[1]],
		True,
		    b = b[[1]];
		    eps = 1;
		    epsa = N[b-a, wprec]/2;
		    epsc = N[c-b, wprec]/2;
		    If[Abs[epsa] < 1, eps = epsa];
		    If[Abs[epsc] < Abs[eps], eps = epsc]
	    ];
	    temp = CPV1[f,{range[[1]],a,{b,eps},c},options];
	    If[ !NumberQ[temp], Return[Fail]];
	    sum += temp
	];
	sum
    ]

arg[z_?(#!=0&)]:=Arg[z]

arg[z_?(#==0&)]:=0


CPV1[f_, {x_, a_, {b_, eps_}, c_}, options___] :=
    Module[{sum, t, fp, fm},
	sum = Abs[N[arg[b - a]] - N[arg[b - c]]];
	If[ sum > 1, sum = 2 N[Pi] - sum];
	If[ sum < .00000000001,
		NIntegrate[f, {x, a, c}, options],
	  (* else *)
		sum = NIntegrate[f, {x, a, b-eps}, options];
		sum += NIntegrate[f, {x, b+eps, c}, options];
		fp = f /. x -> b+t;
		fm = f /. x -> b-t;
		sum + NIntegrate[fp + fm, {t, 0, eps}, options]
	]
    ]

End[] (* NumericalMath`CauchyPrincipalValue`Private` *)

Protect[CauchyPrincipalValue];

EndPackage[] (* NumericalMath`CauchyPrincipalValue` *)


