@Carl Love  A and B having polynomial entries, A.B = Id ==> A invertible, that's all.
BTW, your potential objection is also valit for the provided answer, if A = <p,0;0,1/p>, where det(A) gives 1. But if we disect such a simple problem, everything becomes complicated (more or less artificially).

Actually, you solved the problem! For a)  it's enough to compute  simplify(A . B) , obtaining the unit matrix, without any assumption on p. (Of course it's very useful to know what a singular matrix is.)

@fhn2022 A smooth discontinuous function? Analytic? In mathematics such things do not exist
[f analytic ==> f smooth ==> f continuous].

@fhn2022 Why are you so sure that the function is continuous? You have many complex square roots! Denominators too!

@fhn2022 You may need to increase Digits but then even a longer time will be needed.
If the problem is important for you, try to optimize/simplify the computations.

I just want to mention that a similar problem is here:

https://www.mapleprimes.com/questions/222071-Drawing-An-Apollonian-Gasket#answer239140

@fhn2022 I did not change anything in your code. The first output is just a check with the value of T2 for the those values of alf, y and x. You may delete it.
Only the method of integration in int was changed to method=_CubaCuhre. As you see, it works for alf > 0.277.

@ecterrab  Just want to notice that for an analytic f,  d/dz |f(z)|  exists in a nbd of a point  iff  f is constant.

@jschatzman A standard way to check symbolic results is to compute numerical values for several numerical values of the parameters (like evalf above). If the values agree (sometimes working with high precision if possible), there are big chances to have a correct symbolic result.

@jschatzman You are wrong. The integral for f=+-1/(2*Pi) is Pi/2 and must be considered in the sense of Cauchy principal value.
You can see this using a numerical evaluation:

evalf(Int(sin(x)*exp(-I*x)/x, x = -infinity .. infinity));

                          1.570796327+Float(undefined)*I

@Axel Vogt Actually, for a correct result, Heaviside(0) must be 1/2.

Note however that for f = +- 1/(2Pi), the integral must be considered in the sense of Cauchy principal value.

@Preben Alsholm I do not trust much in solve/identity for such complicated equations.
For example, sol:=solve(identity(EQ,X),[x0,y0,A]);
gives several solutions, but not the correct one: [x0 = 0, y0 = 0, A = 0]

@ecterrab  As mentioned in the definition, `diff/abs`  implements  d/dx |f(x+iy)|,  (for f analytic), not  d/dz |f(z)|. The latter does not exist because it's not a Wirtinger derivative. So,  obviously  d/dx conj(x+iy) = 1.
[for f analytic,  d/dx f = d/dz f].
`diff/abs`  is not suited for a general usage; it works for real variables (as in the provided example). A general solution will have to implement the rest of the nonanalytic functions.

@Carl Love I mean the example without x=y=0. The result of Maple (a+b+x)*(a+x+y)  is different from your simulation.

I don't think that automatic simplification can explain this.

By the way, your use of %%%f  is very ingenious!