1) The expanded radical by itself is not simpler:

    Simplify`SimplifyCount[Sqrt[((1 + x) (x + y))/y]]
    Simplify`SimplifyCount[Sqrt[1 + x] Sqrt[1/y] Sqrt[x + y]]
    (*
      15
      23
    *)

So if such a transformation is tried, it might be rejected out of hand or before it led to a simpler result. And it's possible that such a transformation won't be tried, because it is known that it doesn't simplify the term, even if further down the line it would cancel with something else. The side-alleys that can be pursued in simplification may be infinite and need to be limited in some way.

2) The transformation you seek is PowerExpand. However, because of the issue in 1) above, adding it does not solve your problem:

Simplify[diff[x, y], x > 0 && y > 0, TransformationFunctions -> {Automatic, PowerExpand}]
(*  Sqrt[((1 + x) y)/(x + y)] - Sqrt[((1 + x) (x + y))/y] + x Sqrt[(1 + x)/(x y + y^2)]  *)

3) The trick is to simplify after the expansion, before the PowerExpand transformation is rejected. Both of the following work:

Simplify[PowerExpand@diff[x, y], x > 0 && y > 0]
(*  0  *)

Simplify[diff[x, y], x > 0 && y > 0, TransformationFunctions -> {Automatic, Simplify@*PowerExpand}]
(*  0  *)

4) Here is a custom transformation for distributing the square roots, assuming they are not in the denominator. Note that Sqrt[u] expand to Power[u, 1/2] and 1/Sqrt[u] expands to Power[u, -1/2]. This is important when using pattern-matching, since the second will not match Sqrt[_].

sqrtExpand[e_] := e /. rad : Sqrt[Times[factors__]] :> 
    Thread[rad, Times] /;                             (* Thread[] distributes the Power[..,1/2] over Times *)
      Simplify[And @@ Thread[List @@ factors >= 0]];  (* The condition is that all factors are nonnegative *)
Assuming[x > 0 && y > 0,
 Simplify[diff[x, y], TransformationFunctions -> {Automatic, Simplify@*sqrtExpand}]
 ]
(*  0  *)

5) Note the using of Assuming[x > 0 && y > 0,..] instead of Simplify[.., x > 0 && y > 0]. Simplify does not add its assumptions to $Assumptions; instead it passes them to certain transformations it uses (apparently). However, it does use any assumptions in $Assumptions. Consequently, the simplification in Simplify@*sqrtExpand would not be done with the assumption x > 0 && y > 0 if I used Simplify[.., x > 0 && y > 0]; however, with Assuming, the assumptions will be used in all functions called during the computation that use $Assumptions.

Well, that is about all I can add. I realize it doesn't completely answer your question, and some of the remarks go beyond what you asked; but I hope it gives you some insight that you seek.

You're welcome. Let me know if you have any questions. To inspect the code, consider changing the Module[{...}, line to

Module[{},

with no localized variables. That way, after a test run you can look at the results in each variable. Also change npts = 10 or some other number, so it's comfortable to check testpts.

You might replace

testpts = Transpose[
     RandomReal[#, npts, WorkingPrecision -> precision] & /@ domains
     ];
   AllTrue[testpts, e1 == e2 /. Thread[vars -> #] &]

with

testpts = Transpose[
     Rationalize[RandomReal[#, npts], Max[Abs@#] * 10^-5] & /@ domains
     ];
   AllTrue[testpts, PossibleZeroQ[e1 - e2 /. Thread[vars -> #] &]

Rationalize is used because you need to pass exact values to e1 and e2 for PossibleZeroQ to work. The factor 10^-5 in the second argument (indirectly) limits the size of the denominator (and therefore the numerator, too) in the fractions produced. Big numbers in the fraction can lead to very long computations sometimes. I don't know if 10^-5 is the best choice, but it doesn't seem like it's probably a bad choice.

I suggest it because I think PossibleZeroQ is a little more rigorous than the simple equality check in my first version.