Why does Mathematica choose branches as it does in this situation?

Question

Consider these integrals:

$Assumptions = {e>0,e<1/2};
term1 = Integrate[x/(1-x) (1/(1-4x+4x^2))^e - 1/(1-x), {x, 0, 1}]
term2 = Integrate[1/2 x(1-x) (1/(1-4x+4x^2))^e, {x, 0, 1}]
term3 = Integrate[x/(1-x) (1/(1-4x+4x^2))^e - 1/(1-x)
                  + 1/2 x(1-x) (1/(1-4x+4x^2))^e, {x, 0, 1}]

The third integrand is the sum of the first and second integrands. So it follows that term3 == term1 + term2 should be true.

I won't write down the full results because they're long and mostly irrelevant, but term3 has a nonzero imaginary part, whereas neither term1 nor term2 does:

term3 - term1 - term2 == 2 Pi I

Mathematically, I get why this is a reasonable result (I think): Mathematica has to choose a branch for the exponential function, and in the first integral it chooses the branch that gives a purely real result, whereas in the last case it chooses a different branch. Why would it do that, though? Is there something about these expressions that prompts Mathematica to choose the branch that it does in each specific case? For future integrals, if I run into something like this again, can I tell which branch the integrator is going to pick? Or is it something highly dependent on the details of the algorithm and thus random as far as I'm concerned?

I'm also curious to know if I could input these expressions differently so that Mathematica will choose the real branch of the exponential in every case, but that's not a priority since I know what the answer is, and it may be a matter for a separate question in any case.

Answer 1

FWIW, in V11.2, you get term3 == term1 + term2 to be true:

Block[{$Assumptions = {e > 0, e < 1/2}},
 {term1 = Integrate[x/(1 - x) (1/(1 - 4 x + 4 x^2))^e - 1/(1 - x), {x, 0, 1}],
  term2 = Integrate[1/2 x (1 - x) (1/(1 - 4 x + 4 x^2))^e, {x, 0, 1}],
  term3 = Integrate[
    x/(1 - x) (1/(1 - 4 x + 4 x^2))^e - 1/(1 - x) + 
     1/2 x (1 - x) (1/(1 - 4 x + 4 x^2))^e, {x, 0, 1}]}
 ]
(*
{1/(e (-1 + 2 e))
   2^(-1 - 2 e) ((1 - 2 e) Hypergeometric2F1[2 e, 2 e, 1 + 2 e, 1/
       2] - 2^(1 + 2 e)
       e (-1 - EulerGamma + 2 e EulerGamma - Log[2] + 
        e Log[4] + (-1 + 2 e) PolyGamma[0, 
          2 e] + (π - 2 e π) Tan[e π])), 1/(
 12 - 32 e + 
  16 e^2), -1/(e (-3 + 2 e) (-1 + 2 e))
    2^(-2 (1 + 
       e)) (2 (3 - 8 e + 4 e^2) Hypergeometric2F1[2 e, 2 e, 1 + 2 e, 
        1/2] + 4^e e (11 - 8 e + 12 EulerGamma - 32 e EulerGamma + 
         16 e^2 EulerGamma - 32 e Log[2] + 16 e^2 Log[2] + 4 Log[8] + 
         4 (3 - 8 e + 4 e^2) PolyGamma[0, 2 e] - 
         4 (3 - 8 e + 4 e^2) π Tan[e π]))}
*)

term3 - term1 - term2 // Simplify
(*  0  *)