I think I got it, although I don't know if there is a simpler way to do this:
f(xy) gcd(f(x), y) = f(x) f(y) and f(yx) gcd(f(y), x) = f(y) f(x) imply gcd(f(x),y) = gcd(f(y),x) for all x,y.
Two important properties:
[1] If we let f(x)=n, then n=gcd(f(x),n)=gcd(f(n),x) implying n divides x (notation: n|x). That is, f(x)|x for all x.
[2] Note that gcd(f(x),y)|y and gcd(f(y),x)|x. If gcd(x,y)=1, it follows that gcd(f(x),y)=gcd(f(y),x)=1. Then f(xy) gcd(f(x), y) = f(x) f(y) becomes f(xy) = f(x) f(y). Therefore f(x) must be a

multiplicative function.
Now we check specific values. First, by property [1], f(1) must be 1.
Then we check f(p) for prime p. By property [1], f(p) is either 1 or p:
Case 1: f(p)=1. Then x=y=p gives f(p

^{2}) gcd(f(p),p) = f(p)f(p) ---> f(p

^{2})=1. Likewise, x=p

^{2}, y=p gives f(p

^{3})=1, and so on. By induction, f(p

^{n})=1.
Case 2: f(p)=p. Then x=y=p gives f(p

^{2}) gcd(f(p),p) = f(p)f(p) ---> f(p

^{2})=p. Likewise, x=p

^{2}, y=p gives f(p

^{3})=p, and so on. By induction, f(p

^{n})=p.
Thus, for each prime p, we have a choice whether to assign f(p

^{n})=1 for all n≥1, or f(p

^{n})=p for all n≥1; that is, whether to include p as a prime factor in the output or not. Together with with the multiplicative function property [2], we can summarize possible f as follows: Any function satisfying f(xy) gcd(f(x), y) = f(x) f(y) must be of the form f

_{Q}(x), where Q is a (possibly infinite) subset of the primes, and f

_{Q}(x) is the product of all prime factors of x that are in Q.
----
There is one other thing: The above reasoning does not actually verify that f=f

_{Q}(x) satisfies the equation. To verify it, we need one more thing:
Consider a function f=f

_{Q}(x). To check that f

_{Q}(xy) gcd(f

_{Q}(x), y) = f

_{Q}(x) f

_{Q}(y), we will look at each prime p:
Case 1: p is not in Q: None of the factors in the equation will contain the factor p then.
Case 2: p is in Q but divides neither x nor y: None of the factors in the equation will contain the factor p then.
Case 3: p is in Q and divides x, but not y: The factors f

_{Q}(xy) and f

_{Q}(x) will each contain exactly one factor of p, and only these factors.
Case 4: p is in Q and divides y, but not x: The factors f

_{Q}(xy) and f

_{Q}(y) will each contain exactly one factor of p, and only these factors.
Case 5: p is in Q and divides both x and y: The factors f

_{Q}(xy), gcd(f

_{Q}(x), y), f

_{Q}(x) and f

_{Q}(y) will each contain exactly one factor of p.
In all cases, the factors of p are balanced on both sides. Doing this over all primes p shows that the equation holds.