s = Cos[x] Sin[x]; (*The function to match*)
s0 = {-Sin[x] Sqrt[1 - Sin[x]^2],
Sin[x] Sqrt[1 - Sin[x]^2]};(*the cndidates*)
s1 = FullSimplify[
Reduce[# == s, x, Reals] & /@ s0] ; (*the validity domains*)
hh[x_] :=
Piecewise[
Transpose@{s0,
ReplaceAll[#, ((xx : Element[u__, __]) && z__) ->
Hold@Resolve[Exists[u, xx, z], Reals]] & /@ s1}];
Plot[{ReleaseHold@hh[x], Sin[x] Cos[x]}, {x, 0, 2 Pi},
PlotStyle -> {{Thick, Dashed, Red}, Blue}] // Quiet