>I think you have not given any evidence that solve is changing sqrt(a*b) to sqrt(a) * sqrt(b).
I think that i have such evidence. Let's consider system of the equations:
solve([a*b=(-ma^2/2+Ea*Eb)^2,a=Ea^2-ma^2,b=Eb^2-mb^2],[Eb,a,b]);
(1) a*b=(-ma^2/2+Ea*Eb)^2
(2) a=Ea^2-ma^2
(3) b=Eb^2-mb^2
(sorry, but i can not give more simple example)
To solve this system one need to substitute (2) and (3) into (1). It gives quadratic equation. From this equation one can easy get Eb. It has term sqrt((Ea^2-ma^2)*(ma^2-4*mb^2)) (maxima gives right answer for Eb). Now, to get b one need to substitute Eb into (3). It is easy to see that b should have the term sqrt((Ea^2-ma^2)*(ma^2-4*mb^2)). But maxima says that there is the term sqrt(ma - Ea)*sqrt(ma + Ea)*sqrt(2 mb - ma)*sqrt(2 mb + ma) So, i think that maxima divides sqrt((Ea^2-ma^2)*(ma^2-4*mb^2)) into sqrt(ma - Ea)*sqrt(ma + Ea)*sqrt(2 mb - ma)*sqrt(2 mb + ma)