I. Integrals of the form (∫_((Subscript[z, 1]))^((Subscript[z, 2])))((P(z))/(Q(
                z))) dz

In order to take definite integrals of rational functions in the complex plane, first we decompose the rational function into partial fractions so that the antiderivative F(z) will be represented as a sum of logarithms and integer powers. Thus F(z) will have a regular and simple branch cut structure consisting of half-lines (-∞ + i (Subscript[y, 0]), 
        (Subscript[x, 0]) + i (Subscript[y, 0])) for each logarithm in the antiderivative. Suppose F(z) has a discontinuity at the point (Subscript[z, 0]) of the segment of integration [(Subscript[z, 1]), (Subscript[z, 2])] and is continuous on [(Subscript[z, 1]), (Subscript[z, 0])) and ((Subscript[z, 0]), z2] (continuous from one side at the endpoints). Define

F(
          (Subscript[z, 0]) - 0) = (lim_(z -> (Subscript[z, 0]) - 0)) F(z),  
      F((Subscript[z, 0]) + 0) = (lim_(z -> (Subscript[z, 0]) + 0)) F(z).

Note that all the limits here and below are directional and are taken along the integration path so that (Subscript[z, 0]) ± ϵ lies on the segment [(Subscript[z, 1]), (Subscript[z, 2])]; in other words,

\(TraditionalForm\`\(lim\+\(ϵ -> 0\)\) \(F(
            z\_0 + ϵ)\) ≡ \(lim\+\(ϵ\^′ -> 0\)\) \(F(
            z\_0 + \(ϵ\^′\)(
                z\_2 - z\_1))\), \ ϵ\^′ ∈ ℝ\).

If both limits exist we can make F(z) continuous on [(Subscript[z, 1]), (Subscript[z, 0])] by letting F((Subscript[z, 0])) = F((Subscript[z, 0]) - 0) and continuous on [(Subscript[z, 0]), (Subscript[z, 2])] by letting F((Subscript[z, 0])) = F((Subscript[z, 0]) + 0). We obtain

(∫_((Subscript[z, 1]))^((Subscript[z, 2])))f(
              z) dz = (∫_((Subscript[z, 1]))^((Subscript[z, 0])))f(
                  z) dz + (∫_((Subscript[z, 0]))^((Subscript[z, 2])))f(
                  z) dz=

=(F((Subscript[z, 0]) - 0) - 
                F((Subscript[z, 1]))) + (F((Subscript[z, 2])) - F((Subscript[z, 0]) + 0))=

=F((Subscript[z, 2])) - 
          F((Subscript[z, 1])) - (F((Subscript[z, 0]) + 0) - F((Subscript[z, 0]) - 0)).

This is a generalization of the first fundamental theorem of calculus: we have to subtract the sum of the values of all the 'jumps' of the antiderivative from F((Subscript[z, 2])) - F((Subscript[z, 1])) to obtain the value of the integral. To find the principal value we need to replace F((Subscript[z, 0]) + 0) - F((Subscript[z, 0]) - 0) with (lim_(ϵ -> +0))(
        F((Subscript[z, 0]) + ϵ) - F((Subscript[z, 0]) - ϵ)).

For the terms F(z) = (1/(((z - (Subscript[z, 0])))^p))of the antiderivative we only have to check if (Subscript[z, 0]) belongs to [(Subscript[z, 1]), (Subscript[z, 2])]. If it does the integral doesn't exist in the ordinary sense; if we're looking for the principal value, then for even powers p we can consider the antiderivative to be continuous at (Subscript[z, 0]): if p is even then

(lim_(ϵ -> +0))(
          F((Subscript[z, 0]) + ϵ) - F((Subscript[z, 0]) - ϵ)) = 0.

For F(z) = ln(z - (Subscript[z, 0])) the problem is to find the intersection of the segment [(Subscript[z, 1]), (Subscript[z, 2])] with the half-line (-∞ + i Im (Subscript[z, 0]), 
        Re (Subscript[z, 0]) + i Im (Subscript[z, 0])). If (Subscript[z, 1]) and (Subscript[z, 2]) lie below and above the branch cut respectively and [(Subscript[z, 1]), (Subscript[z, 2])] intersects the branch cut at a point other than (Subscript[z, 0]), then the value of the discontinuity will be 2  π i because the argument of (z - (Subscript[z, 0])) changes by 2  π. If the point of intersection is (Subscript[z, 0]) then

\(TraditionalForm\`\(lim\+\(ϵ -> \(+0\)\)\)(
          F(z\_0 + ϵ) - 
            F(z\_0 - ϵ)) = \(\(\(lim\+\(ϵ -> \(+0\)\)\)(ln(ϵ) - ln(\(-ϵ\)))\)\(=\)\)\)

\(TraditionalForm\`\(\(=\)\(\(lim\+\(ϵ -> \(+0\)\)\) \(i(
              arg\ ϵ - 
                arg(\(-ϵ\)))\) = π\ i\)\)\).

Care should be taken to consider all special cases such as one of the points (Subscript[z, 1]), (Subscript[z, 2]) lying on the branch cut: we have to remember that the branch cut is 'pasted' to the upper half-plane and calculate the discontinuities of F(z) accordingly.

Finally, we need to find the decomposition of a given function f(z) into partial fractions. If

f(z) = ((P(z))/(Q(z))),  
      Q(z) = (∏_(i = 1)^n)(Subscript[(((z - (Subscript[z, i])))^α), i]),  
      (Subscript[z, i]) ≠ (Subscript[z, j]),  deg P < deg Q,

then partial fraction decomposition is

\(TraditionalForm\`\(P(z)\)\/\(Q(z)\) = ∑\+\(i = 1\)\%n\(∑\+\(j = 1\)\%\(α\_i\)A\_\(i\ j\)\/\((z - z\_i)\)\^j\)\)

and each (Subscript[A, i j]) is given by

\(TraditionalForm\`A\_\(i\ j\) = \(1\/\(\((α\_i - 
                    j)\)!\)\) \(lim\+\(z -> 
                z\_i\)\) \(d\^\(α\_i - j\)\/d\ 
                z\^\(α\_i - j\)\) \((\(\(P(z)\)\/\(Q(
                    z)\)\) \((z - z\_i)\)\^α\_i)\)\).

For small values of (Subscript[α, i]):

(Subscript[α, i]) = 1 :  
            (Subscript[A, i 1]) = ((P((Subscript[z, i])))/((Q^′)((Subscript[z, i]))));

(Subscript[α, i])=2 :  (Subscript[A, i 1])=(2/((Q^″)((Subscript[z, i]))))((P^′)((Subscript[z, i]))-P((Subscript[z, i]))(∑_({{j = 1}, {j ≠ i}})^n)((Subscript[α, j])/((Subscript[z, i]) - (Subscript[z, j])))), \(A\_\(i\ 2\) = \(2 \( P(z\_i)\)\)\/\(\(Q\^″\)(z\_i)\)\).

In practice this is not as simple because we have to deal with multiple roots which may, for example, be represented as \(TraditionalForm\`\(\(|\)\(a\)\(|\)\)\) and a under assumption a ≥ 0. Unfortunately, there is no guarantee that simplifications will transform all equal expressions to the same form, so we have to explicitly check if any two roots are equal. To simplify this task, we first sort the roots into lexicographic order: (Subscript[z, 1]) ≲ (Subscript[z, 2]) iff

Re (Subscript[z, 1]) < Re (Subscript[z, 2]) ∨ 
      Re (Subscript[z, 1]) == Re (Subscript[z, 2]) ∧ 
        Im (Subscript[z, 1]) ≤ Im (Subscript[z, 2])

so that equal roots will be grouped. Also we use FactorSquareFree to pull out multiple factors since simply taking list of roots for a polynomial with machine-precision coefficients can lead to two instances of a multiple root being numerically unequal.

Clear[int]
Options[int] = {PrincipalValue -> False, 
          Assumptions -> {}};

    
int[ffz_, {z_Symbol, z1_, z2_}, opt___?OptionQ] := 
      Module[
{init, main, jump, eq, ls, lse, sgn, pvF, 
          assum, ans}, 

SetAttributes[
          init, HoldAll]; 
init[fz_, pz_, qz_, a0_, 
            Lroot_] := 
          Module[
{ordQ, qfact, Lqzfact, 
              Lfactroot}, 
fz = 
              Together[fz]; 
{pz, qz} = {Numerator[fz], 
                Denominator[
                  fz]}; 
If[! PolynomialQ[pz, 
                    z] || ! PolynomialQ[qz, 
                    z], 
Throw[$Failed] ]; 
fz = 
              PolynomialQuotient[pz, qz, z]; 
pz = 
              PolynomialRemainder[pz, qz, z]; 
a0 = 
              Coefficient[qz, z, 
                Exponent[qz, z]]; 
Lqzfact = 
              Rest[FactorSquareFreeList[
                  qz]]; 
Lroot = {}; 
Do[
qfact = 
                Function @@ {z, 
                    Lqzfact[[i, 
                      1]]}; 
Lfactroot = 
                Array[Root[qfact, #] &, 
                  Exponent[qfact[z], z]]; 
Lroot = {Lroot, 
                  Table[Lfactroot, {Lqzfact[[i, 
                        2]]} ] }, 
{i, 
                Length[Lqzfact]}
]; 
ordQ = ls[Re[#1], Re[#2]] || 
                  eq[Re[#1], Re[#2]] && 
                    lse[Im[#1], 
                      Im[#2]] &; 
Lroot = {#[[1]], 
                    Length[#]} & /@ 
Split[
                  Sort[Flatten[Lroot], ordQ], 
                  eq]
]; 

main[] := 
          Module[
{fz = ffz, pz, qz, qz2, Lroot, 
              Lroot2, 
ad, jmp, a0, z0, deg, 
              A}, 
init[fz, pz, qz, a0, 
              Lroot]; 
ad = 
              Integrate[fz, z]; 
jmp = 
              0; 
Do[
{z0, deg} = 
                Lroot[[
                  i]]; 
Do[
If[
                  deg == 1, 
A = 
                    pz/D[qz, z] /. z -> z0, 
Lroot2 = 
                    Delete[Lroot, i]; 
qz2 = 
                    Times @@ ((z - #[[1]])^#[[2]] &) /@ 
                        Lroot2; 
A = 
                    D[pz/qz2, {z, deg - j}]/(a0 (deg - j)!) /. 
                      z -> z0
]; 
If[
                  Catch[eq[A, 0]], Continue[] ]; 
If[
                  j == 1, 
ad += 
                    A Log[z - 
                          z0], 
ad += -A/((j - 
                              1) (z - z0)^(j - 
                                1)) ]; 
jmp += 
                  A jump[z0, j], 
{j, 
                  deg}
], 
{i, 
                Length[Lroot]}
]; 
(ad /. 
                  z -> z2) - (ad /. z -> z1) - 
              jmp
]; 

jump[z0_, deg_] := 
          Module[
{x0, y0, x1, y1, x2, y2, 
              x}, 
{x0, y0} = {Re[z0], 
                Im[z0]}; 
{x1, y1} = {Re[z1], 
                Im[z1]}; 
{x2, y2} = {Re[z2], 
                Im[z2]}; 
If[
              eq[z1, z0] || 
                eq[z2, z0], 
Throw[$Failed] ]; 
If[
              deg ≠ 1, 
If[
                pvF && OddQ[deg] || ! eq[
                      Abs[z0 - z1] + Abs[z2 - z0] - Abs[z2 - z1], 
                      0], 
Return[
                  0] ]; 
Throw[$Failed]
]; 
If[
              eq[y1, y0] && eq[y2, y0], 
If[
                ls[x1, x0] && ls[x2, x0] || 
                  ls[x0, x1] && ls[x0, x2], 
Return[
                  0] ]; 
If[pvF, 
                Return[π I sgn[
                      x1 - x2] ] ]; 
Throw[$Failed]
]; 
If[
              ls[y1, y0] && ls[y2, y0] || 
                lse[y0, y1] && lse[y0, y2], 
Return[
                0] ]; 
x = (y0 - 
                      y1) (x2 - x1)/(y2 - y1) + 
                x1; 
Switch[
              sgn[x - x0], 
1, 
              0, 
-1, 
              2  π I sgn[y2 - y1], 
0, 
              If[pvF, π I sgn[y2 - y1], 
                Throw[$Failed] ]
]
]; 

eq = 
          Module[{ans}, 
ans = 
                Simplify[#1 - #2 == 0, 
                  assum]; 
If[! MatchQ[ans, 
                    True | False], 
                Throw[$Failed] ]; 
ans] &; 
ls = sgn[#1 - #2] < 0 &; 
lse = 
          sgn[#1 - #2] ≤ 0 &; 
sgn = 
          Module[{ans}, 
ans = 
                Simplify[Sign[#], 
                  assum]; 
If[! MatchQ[
                    ans, -1 | 0 | 1], 
                Throw[$Failed] ]; 
ans] &; 

{pvF, 
            assum} = {PrincipalValue, 
              Assumptions} /. 
Flatten[{opt, 
                Options[int]}]; 
ans /; (ans = 
                Catch[main[]]) =!= $Failed
]

The algorithm works for rational functions with approximate as well as exact coefficients.

int[(((x^5) - (x^3) + 4  x)/((((x^2) - 1))^3)), {x, -1 - I, -1 + 
            I}, PrincipalValue -> True] // 
      ComplexExpand
int[(((x^5) - (x^3) + 4  x)/((((x^2) - 1. ))^3)), {x, -1 - 
          I, -1 + I}, PrincipalValue -> True]

We can check the result numerically.

Module[
{ϵ = (10^(-6)), 
        fx = (((x^5) - (x^3) + 4  x)/((((x^2) - 1))^3))}, 
NIntegrate[fx, {x, -1 - I, -1 - I ϵ}, 
          MaxRecursion -> 15, 
          WorkingPrecision -> 30] + 
NIntegrate[
          fx, {x, -1 + I ϵ, -1 + I}, 
          MaxRecursion -> 15, 
          WorkingPrecision -> 30]
]

Also the algorithm can handle some simple assumptions.

int[(1/((x^2) - 1)), {x, a - I a, -a + I a}, 
      Assumptions -> a > 0]

int[(1/((x^2) - 1)), {x, 0, a}, Assumptions -> a > 1, 
      PrincipalValue -> True]

II. Generalization

To be able to handle a wider class of functions, first we perform some transformations on the antiderivative in order to express it in terms of a smaller subset of functions for which the branch cuts structure is known. For example, EllipticE[z] can be converted to Hypergeometric2F1. Branch cuts are described as combinations of a function parametrizing the branch cut as a curve in the complex plane and the domain of the parameter values. A branch cut discontinuity running from -∞ to 0 can be represented as {# &, # < 0 &}, meaning {z | z = a ∧ a < 0}.

Also the algorithm can handle EllipticE[ϕ, z] and EllipticE[z, m] as functions of z (for real values of parameters ϕ and m only). For EllipticE[ϕ, z] the branch cut structure is

\(TraditionalForm\`\(\(|\)\(ϕ\)\(|\)\(\(\(≥\)\(π\/2\)\) : \ 
          z ∈ \((1, ∞)\)\)\)\);

\(TraditionalForm\`\(\(|\)\(ϕ\)\(|\)\(\(\(<\)\(π\/2\)\) : \ 
          z ∈ \((\(csc\^2\) ϕ, ∞)\)\)\)\).

For EllipticE[z, m]:

m < 0 :  x=0∧\[VerticalSeparator]y\[VerticalSeparator]>arcsh  (1/(Sqrt[-m]));

0 < m < 1 :  x=(π/2) + π k∧\[VerticalSeparator]y\[VerticalSeparator]>arcch  (1/(Sqrt[m]));

1 < m :  x=(π/2) + π k∨y=0∧arcsin  (1/(Sqrt[m])) + π k < 
                x < -arcsin (1/(Sqrt[m])) + π(k + 1),

where z = x + i y,  
      k ∈ ℕ.

If the antiderivative contains an expression of the form F(g(z)) where F is discontinuous, then parametrizing the path of integration [(Subscript[z, 1]), (Subscript[z, 2])] as z = ((Subscript[z, 2]) - (Subscript[z, 1])) t + (Subscript[z, 1]) and the branch cut of F as z = B(a) we obtain g(z(t)) = B(a) at the point of intersection. Equating the real and imaginary parts of both expressions we find the values of t and a and determine the points of discontinuity of F(g(z)). For EllipticE[z, m] we need to perform an additional step of finding all the branch cuts that can cross the integration path. For instance, in the case 0 < m < 1,  Re (Subscript[z, 1]) < Re (Subscript[z, 2]) this task can be reduced to finding all the values k such that Re (Subscript[z, 1]) < (π/2) + π k < Re (Subscript[z, 2]).

To determine the direction in which we cross the branch cut, we check the value of g(z(t + ϵ)) for some small positive ϵ (ie further along the intergation path); this method can fail in some situations. To find the left and right limits of the antiderivative at the point (Subscript[z, 0]), we use the following method: we start with F((Subscript[z, 0]) - 0) = F((Subscript[z, 0]) + 0) = F(z) and for each subexpression (Subscript[F, i])(z) having a discontinuity at (Subscript[z, 0]) we replace (Subscript[F, i])(z) with (Subscript[F, i])((Subscript[z, 0]) - 0) in F((Subscript[z, 0]) - 0) and replace (Subscript[F, i])(z) with (Subscript[F, i])((Subscript[z, 0]) + 0) in F((Subscript[z, 0]) + 0). All other subexpressions are simply evaluated at (Subscript[z, 0]). This is not guaranteed to work in general but still it will give correct results in such cases as F(z) being a sum or a product of (Subscript[F, i])(z) even if two or more of (Subscript[F, i]) have a discontinuity at the same point (Subscript[z, 0]).

Also the implementation given here doesn't handle special cases such as one of the points (Subscript[z, 1]), (Subscript[z, 2]) lying on the branch cut, or integration path going through a branch point. This implementation handles integrals for which the antiderivative can be expressed in terms of Power, Log, CosIntegral, CoshIntegral, ChebyshevT, ChebyshevU, LegendreP, EllipticE, Hypergeometric2F1, HypergeometricPFQ, inverse trigonometric and inverse hyperbolic functions (those are converted to logarithms).

Clear[int]
int[fz_, {z_Symbol, z1_?NumericQ, z2_?NumericQ}] := 
      Module[
{finddcf, finddct, finddcz, finddc, finddce, 
          findbc, findconst, 
Lsing = {}, ad, ad0, jmp, 
          fdcz, type, f, farg, param, 
          Lt0}, 

finddcf[$ad_] := 
          Module[
{ad = $ad, fdcz}, 
ad = 
              ad /. 
{EllipticE[x_] -> π/2 
                         Hypergeometric2F1[-1/2, 1/2, 1, 
                          x], 
ChebyshevT[n_, x_] -> 
                      Hypergeometric2F1[-n, n, 
                        1/2, (1 - x)/2], 
ChebyshevU[
                        n_, x_] -> 
                      Hypergeometric2F1[-n, n + 2, 
                        3/2, (1 - x)/2], 
LegendreP[n_, 
                        x_] -> 
                      Hypergeometric2F1[-n, n + 1, 
                        1, (1 - x)/
                          2]} // 
TrigToExp; 
fdcz = 
              Cases[ad, 
s : 
                        expr_^p_Rational | 
(Log | 
                                CosIntegral | CoshIntegral)[
                            expr_] | 
EllipticE[_, expr_] | 
                          EllipticE[
                            expr_, _] | 
HoldPattern[
                            Hypergeometric2F1[___, 
                              expr_]] | 
(HoldPattern[
                                HypergeometricPFQ[L1_, L2_, 
                                  expr_]] /; 
Length[
                                  L1] == 
                                Length[L2] + 
                                  1) /; 
! FreeQ[expr, 
                          z] :> 
{Head[
                        s], 
Function @@ {z, 
                          s}, 
Function @@ {z, 
                          expr}, 
{}}, 
{0, ∞}] // Union; 
fdcz = 
              Replace[fdcz, 
{{EllipticE, 
                      f : Function[z, EllipticE[ϕ_, farg_]], 
                      g : Function[z, 
                          farg_], {}} -> 
{EllipticE1, 
                      f, g, {ϕ}}, 
{EllipticE, 
                      f : Function[z, EllipticE[farg_, m_]], 
                      g : Function[z, 
                          farg_], {}} -> 
{EllipticE2, 
                      f, g, {m}}}, 
{1}]; 
{ad, 
              fdcz}
]; 

finddct[type_, f_, farg_, param_] := 
          Module[
{ans = {}, solQ, g, Lsol, Lsol2, Lcond, 
              Lx, bcparam, bccond, 
rez, imz, rez0, imz0, 
              z0, t, a}, 
solQ[Lsol_] := 
              Head[Lsol] =!= Solve && 
                Lsol =!= {}; 
g = (z2 - z1) # + 
                  z1 &; 
{bcparam, bccond} = 
              findbc[type, param]; 
rez = 
              Re[z2 - z1] t + Re[z1]; 
imz = 
              Im[z2 - z1] t + 
                Im[z1]; 
Do[
Lsol = 
                Solve[farg[
                      z] == bcparam[[i]][a], z]; 
If[! solQ[Lsol], 
                Return[{}] ]; 
Lsol = 
                Union[Lsol]; 
Do[
z0 = 
                  z /. Lsol[[j]]; 
{rez0, imz0} = 
                  Simplify[
ComplexExpand[{Re[z0], Im[z0]}, 
                      TargetFunctions -> {Re, 
                          Im}], 
bccond[[i]][a] && 
                      C[_] ∈ Integers]; 
If[
                  FreeQ[{rez0, imz0}, 
                    ArcTan[_, 0]], 
Lsol2 = 
                    Solve[{rez == rez0, imz == imz0}, {t, 
                        a}], 
Lx = 
                    Cases[{rez0, imz0}, 
                        ArcTan[x_, 0] -> x, ∞] // 
                      Union; 
Lcond = 
                    Distribute[{# < 0, # > 0} & /@ Lx, 
                      List]; 
Lsol2 = Cases[
Solve[
Simplify[{rez == rez0, 
                                  imz == 
                                    imz0}, #] /. 
Abs[
                                  x_] :> 
                                FunctionExpand[
                                  Abs[x], #], 
{t, 
                              a}], 
x_ /; 
                            Im[a /. x] == 0 && 
                              And @@ (# /. x)] & /@ 
                      Lcond; 
Lsol2 = 
                    Flatten[Lsol2, 
                      1]
]; 
If[! solQ[
                      Lsol2], 
                  Continue[] ]; 
Lsol2 = {t, a} /. 
                    Union[Lsol2]; 
Lsol2 = 
                  findconst[Lsol2, type, 
                    bcparam[[i]], 
                    bccond[[i]], 
                    param]; 
Lsol2 = 
                  Cases[Lsol2, {t_, a_, 
                          c_} /; 
Im[t] == 
                            Im[a] == 0 && 
                          0 < t < 
                            1 && 
(bccond[[i]][a] /. 
                              C[1] -> c) && 
Chop[
                              N[farg[g[
                                        t]] - bcparam[[i]][a] /. C[1] -> c]] == 0 -> 
                      t]; 
ans = {ans, 
                    Lsol2}, 
{j, 
                  Length[Lsol]}
], 
{i, 
                Length[bcparam]}
]; 
ans // 
              Flatten
]; 

finddcz[type_, f_, farg_, Lt0_] := 
          Module[
{g, ans = {}, sgn, 
              t, ϵ = 
                10^-6}, 
g = (z2 - z1) # + 
                  z1 &; 
Do[
t = 
                Lt0[[i]]; 
If[
                type === EllipticE2, 
If[
                  Im[farg[g[t]]] == 0, 
sgn = 
                    Sign[Im[farg[g[t + ϵ]]] - 
                        Im[farg[
                            g[t - ϵ]]]], 
sgn = 
                    Sign[Re[farg[g[t + ϵ]]] - 
                        Re[farg[
                            g[t - ϵ]]]]
], 
sgn = 
                  Sign[Im[
                      f[g[t + ϵ]]]]
]; 
AppendTo[ans, {g[t], sgn}], 
{i, 
                Length[Lt0]}
]; 
ans
]; 

finddc[$Lsing_] := 
          Module[
{ordQ, Lsing = $Lsing, ans = 0, up, lo, 
              type, f, farg, param, z0, sgn}, 
ordQ = 
              Re[#1] < Re[#2] || 
                  Re[#1] == Re[#2] && 
                    Im[#1] ≤ Im[#2] &; 
Lsing = 
              Sort[Lsing, 
                ordQ[#1[[5]], #2[[5]] ] &]; 
Lsing = 
              Split[Lsing, #1[[5]] == #2[[5]] &]; 
Do[
up = lo = 
                  ad; 
Do[
{type, f, 
                    farg, param, z0, sgn} = 
                  Lsing[[i, 
                    j]]; 
If[
                  type === EllipticE2, 
{up, lo} = 
                    finddce[f, farg, param, z0, sgn, up, 
                      lo], 
up = 
                    up /. f[z] -> 
                        Re[f[z0]] + 
                          I sgn Abs[Im[f[z0]]]; 
lo = 
                    lo /. f[z] -> 
                        Re[f[z0]] - 
                          I sgn Abs[
                              Im[f[z0]]]
], 
{j, 
                  Length[Lsing[[i]] ] }
]; 
{up, lo} = {up, lo} /. 
                  z -> z0; 
ans += 
                up - lo, 
{i, 
                Length[Lsing]}
]; 
ans
]; 

finddce[f_, 
            farg_, param_, z0_, $sgn_, up_, lo_] := 
          Module[
{sgn = $sgn, Earg, newup, newlo, m, k, 
              mid, rs, is, Em}, 
Earg = 
              farg[z0]; 
m = 
              param[[1]]; 
k = IntegerPart[Earg/π]; 
{rs, is} = 
              Sign[{Re[Earg], Im[Earg]}]; 
Em = 
              EllipticE[
                m]; 
Which[
Im[
                  Earg] == 0, 
If[
                Abs[Re[Earg]] < π/
                      2 + π Abs[k], 
sgn *= 
                  rs; 
mid = (2  k + rs) Em - 
                    I rs Im[Em]; 
{newup, 
                    newlo} = {2  mid - f[z0], 
                    f[z0]}, 
sgn *= -rs; 
mid = (2  k + rs) Em + 
                    I rs Im[Em]; 
{newup, 
                    newlo} = {2  mid - f[z0], 
                    f[z0]}
], 
m < 
                
                0, 
sgn *= -is; 
{newup, newlo} = {-Conjugate[f[z0]], 
                  f[z0]}, 
0 < m < 
                1, 
sgn *= 
                is; 
{newup, 
                  newlo} = {(rs + is) Em + 
                    f[z0], (4  k + rs - is) Em - 
                    Conjugate[
                      f[z0]]}, 
True, 
sgn *= 
                is; 
{newup, 
                  newlo} = {(rs + is) Em + 
                    f[z0], 
(4  k + rs - is) 
                      Re[Em] + (rs - is) I Im[Em] - 
                    Conjugate[
                      f[z0]]}
]; 
{newup, 
                newlo} = ExpandAll[{newup, newlo}]; 
If[
              sgn == 1, 
{up /. f[z] -> newup, 
                lo /. f[z] -> newlo}, 
{up /. 
                  f[z] -> newlo, 
                lo /. f[z] -> 
                    newup}]
]; 

findbc[type_, param_] := 
          Module[
{ϕ, m}, 
Switch[
              type, 
Hypergeometric2F1 | 
                HypergeometricPFQ, 
{{# &}, {# > 
                      1 &}}, 
EllipticE1, 
ϕ = 
                param[[1]]; 
{{# &}, {If[
                    Abs[ϕ] ≥ π/2, # > 1 &, # > 
                        Csc[ϕ]^2 &]}}, 
EllipticE2, 
m = 
                param[[1]]; 
Which[
m < 
                  0, {{I  # &}, {Abs[#] > 
                        ArcSinh[1/Sqrt[-m]] &}}, 
0 < 
                  m < 1, {{π/2 + π C[1] + I  # &}, {Abs[#] > 
                        ArcCosh[1/Sqrt[m]] &}}, 
m > 
                  1, {{π/2 + π C[1] + 
                        I  # &, # &}, 
{# ∈ 
                        Reals &, 
                    ArcSin[1/Sqrt[m]] + π C[1] < # < -ArcSin[
                              1/Sqrt[m]] + π (C[1] + 
                                1) &}}
], 
_, {{# &}, {# < 
                      0 &}}
]
]; 

findconst[Lsol_, type_, bcparam_, 
            bccond_, param_] := 
          Module[
{ans = {}, m, Lsolc, rng, c, t, tt, 
              a}, 
If[
              type =!= EllipticE2, 
Return[Append[#, 
                      C[1]] & /@ 
                  Lsol] ]; 
Do[
{t, a} = 
                Lsol[[
                  i]]; 
If[! FreeQ[
                    bccond, C[1]], 
m = 
                  param[[1]]; 
c = 
                  Floor[(a - 
                          ArcSin[
                            1/Sqrt[m]])/π]; 
AppendTo[
                  ans, {t, a, C[1]} /. 
                    C[1] -> c], 
If[! FreeQ[t, 
                      C[1]], 
c = 
                    C[1] /. Solve[tt == t, 
                          C[1]][[1]]; 
rng = 
                    Sort[c /. {{tt -> 0}, {tt -> 
                              1}}]; 
Lsolc = {t, a, 
                        C[1]} /. 
Table[{C[1] -> 
                            i}, {i, 
                          Ceiling[
                            rng[[1]] ], 
                          Floor[rng[[2]] ] } ]; 
ans = 
                    Join[ans, Lsolc], 
AppendTo[
                    ans, {t, a, 
                      C[1]}]
]
], 
{i, 
                Length[Lsol]}
]; 
ans
]; 

ad0 = 
          Integrate[fz, 
            z]; 
If[! FreeQ[ad0, Integrate], 
          Return[$Failed] ]; 
{ad, fdcz} = 
          finddcf[ad0]; 
Do[
{type, f, 
              farg, param} = 
            fdcz[[i]]; 
Lt0 = 
            finddct[type, f, farg, 
              param]; 
AppendTo[
                Lsing, {type, f, farg, param, 
                  Sequence @@ #}] & /@ 
finddcz[type, f, 
              farg, Lt0], 
{i, 
            Length[fdcz]}
]; 
jmp = 
          finddc[Lsing]; 
(ad0 /. 
              z -> z2) - (ad0 /. z -> z1) - 
          jmp
]

Here is a quite complicated example where the integrand has three discontinuities along the integration path.

int[(Sqrt[z - (Sqrt[z])]), {z, -(3/4) - 2  I, -(3/4) + 2  I}] // 
      N
NIntegrate[(Sqrt[z - (Sqrt[z])]), {z, -(3/4) - 2  I, -(3/4) + 
          2  I}]

Here integral comes out in terms of CosIntegral.

int[
        Cos[1 - 2  t + I (t - 1)]/(1 - 2  t + I (t - 1)), {t, 0, 
          2}] // N
NIntegrate[
      Cos[1 - 2  t + I (t - 1)]/(1 - 2  t + I (t - 1)), {t, 0, 
        2}]

Here integral is expressed in terms of EllipticE. We have to use arbitrary-precision arithmetics to calculate the value of EllipticE on the branch cut correctly.

int[(Sqrt[1 + 2 Sinh[z]^2]), {z, -(3/4), (5/4) - I π}] // 
        N[#, 20] & // N
NIntegrate[(Sqrt[1 + 2 Sinh[z]^2]), {z, -(3/4), 
        (5/4) - I π}]

Copyright © 2002, Maxim Rytin