How can I separate a separable function

I would take a logarithmic derivative with respect to one variable - it should be then independent of the other one, then integrate it back over the first variable and exponentiate. The second function is found by plain division. Here is the code:

ClearAll[getGX];
getGX[expr_, xvar_, yvar_] :=
  With[{dlogg = D[Log[expr], xvar] // FullSimplify},
     Exp[Integrate[dlogg, xvar]] /; FreeQ[dlogg, yvar]];

Clear[getHY];
getHY[expr_, xvar_, yvar_] := FullSimplify[(#/getGX[#, xvar, yvar]) &[expr]]

A test function:

ftest[x_, y_] := (x^2 + 1)*y^3 *Exp[-x - y] 

Now,

getGX[ftest[x,y],x,y]

(* E^-x (1+x^2)  *)

getHY[ftest[x,y],x,y]

(* E^-y y^3 *)

The integration constant ambiguity translates into an ambiguity of how you split the function, since this operation is only defined up to a multiplicative constant factor by which you can multiply one function, and divide the other one.

What I'd do:

gaussian[x_, y_] := 1/(2 π σ^2) Exp[-((x^2 + y^2)/(2 σ^2))];
f[x_, y_] = D[gaussian[x, y], x, y]

Exp[Select[Expand[PowerExpand[Log[Together[f[x, y]]]]], #]] & /@
    {FreeQ[#, x | y] &, ! FreeQ[#, x] &, ! FreeQ[#, y] &}
{1/(2 Pi σ^6), E^(-(x^2/(2 σ^2))) x, E^(-(y^2/(2 σ^2))) y}

The snippet separates out the constant factor, the factors with x, and the factors with y.

More examples:

f = -Pi Cos[x]^2 Sin[y]^3/E;
Exp[Select[Expand[PowerExpand[Log[Together[f]]]], #]] & /@
    {FreeQ[#, x | y] &, ! FreeQ[#, x] &, ! FreeQ[#, y] &}
{-Pi/E, Cos[x]^2, Sin[y]^3}

We see that negative constant factors are reproduced.

f = w[x] z[y];
Exp[Select[Expand[PowerExpand[Log[Together[f]]]], #]] & /@
    {FreeQ[#, x | y] &, ! FreeQ[#, x] &, ! FreeQ[#, y] &}
{1, w[x], z[y]}

The implicit constant factor of 1 is detected.

While the two existing answers are good, I've always approached this sort of problem with FactorList and there should be an answer showing how to use this old Mathematica function in this way. Much time having passed since this question was asked, I'll use current tools to process the list of factors.

I decided to separate the task into two parts, which makes the solution a little more complicated than necessary. OTOH, I've sometimes wanted to know the factors, including ones that cannot put in the form $g(x)h(y)$, separated by type.

ClearAll[sepFac, sepParts, dependsOn];
sepParts[expr_, vars_List] := With[{
      factors = GroupBy[FactorList[expr], 
        Function[fac, 
         Pick[dependsOn @@ vars, Internal`DependsOnQ[fac, #] & /@ vars]]
        ],
      apply = Function[{h, lev}, Apply[h, #, lev] &]},
   apply[Times, {1}]@ apply[Power, {2}]@ factors
   ];
sepFac::nonsep = "The factor `` could not be separated.";
sepFac[expr_, vars_List] := With[{parts = sepParts[expr, vars]},
   Message[sepFac::nonsep, #] & /@ 
    Values@KeySelect[parts, Length[#] >= 2 &];
   (MapAt[dependsOn[] # &, (* multiply first by constant factor *)
     dependsOn /@ vars,
     1] /. parts
     ) /; Max[Length /@ Keys[parts]] <= 1
   ];

Examples:

gaussian[x_, y_] := 1/(2 π σ^2) Exp[-((x^2 + y^2)/(2 σ^2))]
sepFac[D[gaussian[x, y], x, y], {x, y}]

(*
  {(E^(-(x^2/(2 σ^2))) x)/(2 π σ^6), 
   E^(-(y^2/(2 σ^2))) y}
*)

sepFac[(x + y) x y, {x, y}]

sepFac::nonsep :  The factor x+y could not be separated.

(*  sepFac[x y (x + y), {x, y}]  *)

sepParts[(x + y) x y, {x, y}]

(*
  <|dependsOn[] -> 1, dependsOn[x] -> x, dependsOn[y] -> y, 
   dependsOn[x, y] -> x + y|>
*)