It seems to me that the OP almost contains the answer: the gcd of $f$ and $f^\prime$ (let's assume characteristic zero) contains all the irreducible factors of $f$ which appear with exponent greater than $1.$ This is surely enough to figure out if the polynomial is a square.
EDIT to answer the revised version of the question:
Write down $$\sum_{i=0}^n a_i x^i = (\sum_{j=0}^{n/2} b_j x^j)^2.$$
This gives a collection of $n+1$ quadratic equations in $3n/2 + 2$ variables. You now eliminate the $b_j$ to get the variety of perfect squares. Needless to say, this is not algorithmically very pleasant (the degree is going to be exponential in $n$), but you can use successive resultants or Grobner bases to do it for small degrees, and you might see a pattern.
Another Edit
If you have Mathematica, you can perform the above-mentioned experiments with the program below:
genpoly[deg_, name_, var_] := Sum[name[i] var^i, {i, 0, deg}]
quadraticeq[deg_, name1_, name2_, var_]:= Eliminate[MapThread[Equal, {CoefficientList[genpoly[2deg, name1, var], var],CoefficientList[Expand[genpoly[deg, name2, var]^2], var]}], Table[name2[i], {i, 0, deg}]]
(for example, to see what the variety is describing quadratic polynomials which are squares, you do:
quadraticeq[1, a, b, x]
a and b are dummy variables, a[0], ..., a[2 deg] are the variables you care about. For quadratic polynomials you get (no surprise):
4 a[0] a[2]==a[1]^2
While for quartic polynomials you get:
a[0] a[3]^2==a[1]^2 a[4]&&-4 a[0] a[1] a[2]+8 a[0]^2 a[3]==-a[1]^3&&8 a[0] a[3] a[4]==a[1] (-a[3]^2+4 a[2] a[4])&&16 a[0] a[4]^2==-a[2] a[3]^2+4 a[2]^2 a[4]-2 a[1] a[3] a[4]&&8 a[1] a[4]^2==a[3] (-a[3]^2+4 a[2] a[4])&&a[0] (-4 a[2]^2+2 a[1] a[3])+16 a[0]^2 a[4]==-a[1]^2 a[2]&&a[0] (-4 a[2] a[3]+8 a[1] a[4])==-a[1]^2 a[3]
which is a little more painful.