What does ReduceFreeVariables really do? And can we rely on it?

Very frequently, in answers on this site, the undocumented function ReduceFreeVariables 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:

ReduceFreeVariables[expr1]
(*  {a, m, n, rho, u, v}  *)

ReduceFreeVariables[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 ReduceFreeVariables do?
2. Is it reliable?
• I think I know, but I don't think I can explain it. Certainly not tonight. Hint/intuition: Start at the top level. If the head is Plus or Times, consider the arguments as a list of elements: if there are no interdependencies treating the elements as variables, then you're done; however if one of the terms is a and another has a in one of its arguments, the refine and repeat until no more changes. (I don't know how to explain "refine," but I suspect x[0,a^2,u,v,m] is replaced by a^2, u, v, m. BTW, this is all guesswork, anyway.) Math fns. like Sin[] might get special treatment. Jan 2, 2020 at 5:49
• Yes, math functions get special treatment: ReduceFreeVariables[Sin[x] + Cos[x]] vs. ReduceFreeVariables[f[x] + g[x]] Jan 2, 2020 at 5:53
• Variables[Level[expr1, {-1}]] and Variables[Level[expr2, {-1}]] give the same outputs. Jan 2, 2020 at 16:03
• But Variables[Level[x[1] + x[2], {-1}]] does give the desired output for indexed variables. (Both Variables[x[1] + x[2]] and ReduceFreeVariables[x[1] + x[2]] do.) Jan 3, 2020 at 1:29

Here's an explanation by example (update below). It's a function that seems to do what ReduceFreeVariables 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_] := InternalDependsOnQ[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 InternalDependsOnQ. 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]
ReduceFreeVariables[expr1]
(*
{a, m, n, rho, u, v}
{a, m, n, rho, u, v}
*)

getFreeVariables[expr2]
ReduceFreeVariables[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]
ReduceFreeVariables[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]
ReduceFreeVariables[exprtmp2]
(*
{True, w, y, f[Sin[x]], g[z]}
ReduceFreeVariables[True - 3 Cos[2 + y] + f[Sin[x]] + g[z] + h[w, y]]
*)

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

exprtmp2 = exprtmp + And;
getFreeVariables[exprtmp2]
ReduceFreeVariables[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 ReduceFreeVariables.

Example:

ReduceFreeVariables[x^y < 4]
ReduceFreeVariables[x^y < 4, "Algebraic"]
ReduceFreeVariables[{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.