What does Reduce`FreeVariables really do? And can we rely on it?

Question

Very frequently, in answers on this site, the undocumented function Reduce`FreeVariables comes up in the context of generating a list of variables inside an expression. It is especially recommended in answers and comments by reputable members Carl Woll (here, here, here, ...) and J.M. (here, here, here, here, ...).

Based on the name FreeVariables and how it has been used on this site, I was under the impression that it only generates a List of Symbols appearing in its input expression. I went so far as to document it as such in my post of undocumented functions (under 'Some more stuff').

However, it is not the case that it generates a List of Symbols. Indeed, it can output more complicated expressions. And, it even returns drastically different results with very little modification to the input. Consider the following very similar expressions:

expr1 = f[g[d[rho]],g[d[s[rho, m, n]]]]*(-x[0,a^2,u,v,m]-2*a*x[1,a^2,u,v,m]);

expr2 = f[g[d[rho]],g[d[s[rho, m, n]]]]*(-x[0,a^2,u,v,m]-2(*a*)*x[1,a^2,u,v,m]);

In the expr2 a factor of a is missing in front of the x[...] function. Yet, the output of ReduceFreeVariables` is quite different for the two examples:

Reduce`FreeVariables[expr1]
(*  {a, m, n, rho, u, v}  *)

Reduce`FreeVariables[expr2]
(*  {f[g[d[rho]],g[d[s[rho, m, n]]]], x[0,a^2,u,v,m], x[1,a^2,u,v,m]}  *)

Ordinarily, this behavior should be no cause for alarm. After all, it is an undocumented function. However, given the level of recommendation by various users on this site for this function, I would like to ask the questions (which I direct especially at the SE members who use it):

1. What exactly does Reduce`FreeVariables do?
2. Is it reliable?

Answer 1

Here's an explanation by example (update below). It's a function that seems to do what Reduce`FreeVariables does. Please let me know if you break it.

(* check dependence of x on y
 * by construction x does not depend on x *)
dependsOnQ[x_][y_] := dependsOnQ[x, y];
dependsOnQ[x_, x_] := False; (* b/c we apply Or in the next def *)
dependsOnQ[x_, y_List] := Or @@ dependsOnQ[x] /@ y;
dependsOnQ[x_, y_] := Internal`DependsOnQ[x, y];
(* True if f is a NumericFunction *)
numericFunctionQ[f_Symbol] := MemberQ[Attributes@f, NumericFunction];
numericFunctionQ[f_] := False;
(* get variables (Variables[] ignores powers) *)
myVariables[f_Symbol?numericFunctionQ[x__]] := myVariables[{x}];
myVariables[list_List] := Union[myVariables /@ list // Flatten];
myVariables[x_?NumericQ] := Sequence[];
myVariables[other_] := other;
(* iterative step: check current Variables[expr] to see 
 * if they are free with respect to each other;
 * replace dependent ones by their arguments   *)
replaceTopLevelUnfreeVariables[expr_] := With[{vars = myVariables[expr]},
   Union@DeleteCases[Replace[vars,
      {e_ /; dependsOnQ[e, vars] :> Sequence @@ e}, 1
      ], x_?NumericQ]
   ];
(* Iterate replaceTopLevelUnfreeVariables until all are free *)
getFreeVariables[expr_] := 
  FixedPoint[replaceTopLevelUnfreeVariables, {expr}];

Update: What getFreeVariables[expr] does is start from the top level with a list {expr} and descends until all expressions in the list do not have the a NumericFunction head and are not dependent on any of the other expressions, as defined by Internal`DependsOnQ. If the head is a NumericFunction, it is replaced with its arguments. (There is some management so that the list of current variables is kept flat and free of duplicates.) If the head is not a NumericFunction but is of the form e = f[x1, x2,...], then it is checked to see if e is dependent on the other current variables in the list. If it is dependent on some of them, it is replaced by its arguments x1, x2,..., with numeric arguments discarded. The list of current variables is then reevaluated. The process repeats until the list stops changing. (myVariables is similar to Variables, except that Variables does not reduce all numeric functions to their arguments.)

Test examples:

getFreeVariables[expr1]
Reduce`FreeVariables[expr1]
(*
  {a, m, n, rho, u, v}
  {a, m, n, rho, u, v}
*)

getFreeVariables[expr2]
Reduce`FreeVariables[expr2]
(*
  {f[d[rho]^3, d[s[rho, m, n]]^3], x[0, a^2, u, v, m], 
   x[1, a^2, u, v, m]}
  {f[d[rho]^3, d[s[rho, m, n]]^3], x[0, a^2, u, v, m], 
   x[1, a^2, u, v, m]}
*)

exprtmp = f@Sin[x] - 3 Cos[y + 2] + g[z] + h[w, y];
getFreeVariables[exprtmp]
Reduce`FreeVariables[exprtmp]
(*
  {w, y, f[Sin[x]], g[z]}
  {w, y, f[Sin[x]], g[z]}
*)

Adding True or False breaks it, but adding some other system symbols does not. I'm not sure why the boolean constants are special cases, and other special cases might exist. One just has to add such special cases to getFreeVariables to fix it. It would be genuinely interesting if both gave a result and the results were different.

exprtmp2 = exprtmp + True;
getFreeVariables[exprtmp2]
Reduce`FreeVariables[exprtmp2]
(*
  {True, w, y, f[Sin[x]], g[z]}
  Reduce`FreeVariables[True - 3 Cos[2 + y] + f[Sin[x]] + g[z] + h[w, y]]
*)

exprtmp2 = exprtmp + Cos;
getFreeVariables[exprtmp2]
Reduce`FreeVariables[exprtmp2]
(*
  {Cos, w, y, f[Sin[x]], g[z]}
  {Cos, w, y, f[Sin[x]], g[z]}
*)

exprtmp2 = exprtmp + And;
getFreeVariables[exprtmp2]
Reduce`FreeVariables[exprtmp2]
(*
  {And, w, y, f[Sin[x]], g[z]}
  {And, w, y, f[Sin[x]], g[z]}
*)

Update:

@rogerl in this answer revealed an optional second argument to Reduce`FreeVariables.

Example:

Reduce`FreeVariables[x^y < 4]
Reduce`FreeVariables[x^y < 4, "Algebraic"]
Reduce`FreeVariables[{x^y < 4, x}, "Algebraic"]
(*
  {x, y}
  {x^y}
  {x, y}
*)

It seems that when "Algebraic" is specified, transcendental functions are not reduced unless they depend on another free variable.