Identity with nested sum taken over divisors of $\gcd$'s For computational reasons,
I want to show that the following holds true:
Let $n_1,n_2,N\in \mathbb{N}$. One has
$$\Large
\sum_{a\mid \gcd(n_1,N)}\sum_{b\mid \gcd(n_2,\frac{n_1N}{a^2})} ab
=\sum_{g\mid \gcd(n_1,n_2)}\sum_{d\mid \gcd(\frac{n_1n_2}{g^2},N)}gd.
$$
It is true numerically, but I do not know how to prove it mathematically
Any comment ?
Thanks. 
PS:  Note that all sums in the above equation runs over positive divisors only.
 A: barto's comment to the OP's question suggest exchanging the names of $n_1$ and $N$ : the question becomes
$$ \sum_{a|(n_1, N)}\sum_{b|(n_2, n_1N/a^2)}ab \overset{?}{=} \sum_{g|(n_2, N)}\sum_{d|(n_1, n_2N/g^2)}gd .$$
In other words, the left-hand side must be unchanged if we swap $n_1$ and $n_2$. In terms of the generating series, the identity you want to show is therefore equivalent to proving that for each $N\geq 1$, the function
$$ (1) \qquad F(s, w) = \sum_{n_1, n_2\geq 1}n_1^{-s} n_2^{-w} \sum_{a|(n_1, N)}\sum_{b|(n_2, n_1N/a^2)} ab $$
satisfies $F(s, w) = F(w, s)$. That's by unicity of the coefficients of a Dirichlet series, granted that it has finite abscissa of absolute convergence. To check this, note that the coefficients of $F(s, w)$ are positive, and
$$ F(s, w) \leq \sum_{n_1, n_2\geq 1}n_1^{-s}n_2^{-w} \sum_{a|n_1}a\sum_{b|n_2}b \leq \sum_{n_1, n_2\geq 1}d(n_1)n_1^{1-s}d(n_2)n_2^{1-w} = \zeta(s-1)^2\zeta(w-1)^2 $$
where $d(n)$ stands for the divisor counting function. If $s$ and $w$ are both $>2$, the left hand side is finite so this implies absolute convergence of $F(s, w)$, and you are guaranteed unicity of coefficients.
The standard procedure to evaluate the LHS of (1) to switch summation, and use the fact that $(a, b)|c$ is equivalent to $a|c$ and $b|c$.

Push the summations over $n_1$ and $n_2$ to the end. You have
$$ F(s, w) = \sum_{\substack{a, b\geq 1\\ a|N}}\sum_{\substack{n_1, n_2\geq 1\\ a|n_1, b|n_2 \\ a^2b | n_1N}} abn_1^{-s}n_2^{-w}. $$
Rewrite $n_1 = an_1'$ and $n_2=bn_2'$, then
$$ F(s, w) = \sum_{\substack{a, b\geq 1\\ a|N}}\sum_{\substack{n_1', n_2'\geq 1\\ ab | n_1'N}} a^{1-s}b^{1-w}(n_1')^{-s}(n_2')^{-w}. $$
Note that
$$ ab|n_1'N \Leftrightarrow \frac{ab}{(ab, N)}\Big|n_1'\frac{N}{(ab, N)} $$
so that actually $ab/(ab, N) | n_1'$. Rewriting $n_1' = (ab/(ab, N))n_1''$, you get
$$ F(s, w) = \sum_{\substack{a, b\geq 1\\ a|N}}a^{1-2s}b^{1-w-s}(ab, N)^s \sum_{\substack{n_1'', n_2'\geq 1}} (n_1'')^{-s}(n_2')^{-w} = H(s,w)\zeta(s)\zeta(w) $$
where
$$ H(s, w) = \sum_{\substack{a, b\geq 1\\ a|N}}a^{1-2s}b^{1-w-s}(ab, N)^s $$
and we want to show that $H(s, w)=H(w, s)$. Since $a|N$, then $(ab, N)=a(b, N/a)$ and so
$$ H(s, w) = \sum_{\substack{a, b\geq 1\\ a|N}}a^{1-s}b^{1-w-s}(b, N/a)^s $$
Let $u=(b, N/a)$, and write $b=ub'$, then
$$ H(s, w) = \sum_{\substack{a, b', u\geq 1\\au|N \\(b', N/(au))=1}}u^{1-w}a^{1-s}(b')^{1-w-s} .$$
Upon writing $d=au$, this is rewritten as
$$ H(s, w) = \sum_{\substack{d|N, b\geq 1 \\ (b, N/d)=1}} \Big(\sum_{d=au}a^{1-s}u^{1-w}\Big)b^{1-w-s} $$
whose symmetry is obvious, because $a$ and $u$ play the same role.
You can probably reverse-engineer the changes of variables to deduce a magic 1-1 correspondance between $(a, b)$ and $(d, g)$ in your original sum and obfuscate the rest.
A: $\newcommand{\lcm}[0]{\mathrm{lcm}}$Allow me to rewrite this as
$$\large
\sum_{a \mid \gcd(x, z)}
\ 
\sum_{b \mid \gcd(y,\frac{x z}{a^2})} a b
=
\sum_{c \mid \gcd(x, y)}
\ 
\sum_{d\ \mid \gcd(z, \frac{x y}{c^2})} c d.
$$ 
I will provide a proof (barring mistakes) in the special case when $\gcd(x, y) = 1$. I have some hope to reduce the general case to this one.
So in this case we have $\lcm(x, y) = x y$, and thus
$$
\gcd(z, x y) = \gcd(z, x) \gcd(z, y),
$$
where the two factors are coprime.
One sees thus that it is enough to prove that
$$
\gcd(y, z) = \gcd(y, \frac{x z}{a^{2}})
$$
for each $a \mid \gcd(x, z)$.
Clearly if $t \mid \gcd(y, \dfrac{x z}{a^{2}})$, then $t \mid x z$, and as $t \mid y$, we have $\gcd(t, x) = 1$, so that $t \mid z$. It follows that $t \mid \gcd(y, z)$, or
$$
\gcd(y, \frac{x z}{a^{2}}) \mid \gcd(y, z).
$$
Conversely, suppose $t \mid \gcd(y, z)$. Since $a \mid x$, we will have $\gcd(t, a) = 1$. Now
$$
t \mid x z = a^{2} \cdot \frac{x z}{a^{2}},
$$
and thus $t \mid \dfrac{x z}{a^{2}}$. This shows that
$$
\gcd(y, z) \mid \gcd(y, \frac{x z}{a^{2}}).
$$
A: I do not have a proof, but I will try to classify this problem in terms of linear operators. Maybe someone later can prove the claimed identity using this idea.
Let $f:\mathbb{N}\to \mathbb{R}$ be a multiplicative arithmetic function. 
Define a linear operator $T(n)$ on $f$ as follows:
$$
(T(n)f)(N):=\sum_{a\mid \gcd(n,N)}f(\frac{Nn}{a}).
$$
It is clear that $(T(1)f)(N)=f(N)$.
For two positive integers $n_1,n_2$ one has
$$
\big(T(n_2)\circ T(n_1)f\big)(N)= \sum_{a\mid \gcd(n_1,N)}\sum_{b\mid \gcd(n_2,\frac{Nn}{a^2})}f(\frac{Nn_1n_2}{ab}).
$$
On the other hand, one has
$$
\Big(\sum_{g\mid \gcd(n_1,n_2)}(T(\frac{n_1n_2}{g^2})f)(N)\Big)=
\sum_{g\mid \gcd(n_1,n_2)}\sum_{d\mid \gcd(\frac{n_1n_2}{g^2},N)}f(\frac{n_1n_2}{g d}).
$$
Comparing the last two equations with the original identity in the first post, we can now say that the proof of this identity is equivalent to proving that 
$$
T(n_2)\circ T(n_1)=\sum_{g\mid \gcd(n_1,n_2)}T(\frac{n_1n_2}{g^2}).
$$
Remark
Some functions like $\sigma$ divisor function, and $\tau$ function satisfy similar identities.
A: Dirichlet convolution looks potentially good here, but I found a simpler approach.  Let $n_3 = N$, so the problem reads more symmetrically:
$$ \Large
\sum_{a\mid \gcd(n_1,n_3)}\sum_{b\mid \gcd(n_2,\frac{n_1 n_3}{a^2})} ab
=\sum_{g\mid \gcd(n_1,n_2)}\sum_{d\mid \gcd(\frac{n_1n_2}{g^2},n_3)}gd $$

For some reason this looks like a Venn diagram problem or like Associtivity.


*

*$ (A \cap B) \cap C = A \cap (B \cap C)$ or

*$ (A+B)+C = A+(B+C)$

*$ (f \ast g)\ast h = f \ast (g \ast h)$ (this is dirichlet convolution)

*An appropriate identity using mutual information

What is being intersected or convolved here?  How about the following definitions which I figured out after the computation below:


*

*Let $A$ be the set of divisors of $a$, $B$ be the set of divisors of $b$

*Let $N_k$ be the set of divisors of $n_k$, with $k = 1, 2, 3$


Looking at the problem we have to compare these two rather complicated sets.
$$ \left\{ab: a\big| n_1,\, a\big|n_3,\, b\big|n_2,\, b\Big| \frac{n_1}{a}\cdot\frac{n_3}{a} \right\} \text{ or }
\left\{gd: g\big| n_1,\, g\big|n_2,\, d\big|n_3,\, d\Big| \frac{n_1}{g}\cdot\frac{n_2}{g} \right\}$$
In terms of set theory we have that $A \subseteq N_1 \cap N_3$ and the following set theory computation:
\begin{eqnarray*}
B &\subseteq & N_2 \cap \bigg( (N_1 \backslash A) \cup (N_3 \backslash A) \bigg) 
\\
&=& N_2 \cap \bigg( (N_1 \cup N_3 ) \backslash A) \bigg) \\
&=& \bigg(N_2 \cap (N_1 \cup N_3) \bigg)\backslash A
\end{eqnarray*}
the set of divisors of $ab$ are going to the the union of divisors of $A\cup B$, yet
$$ A \cup B \subseteq (N_1 \cap N_3) \cup (N_2 \cap N_3) \cup(N_1 \cap N_2) $$
Note these unions and intersections have multiplicity so we are dealing with multisets.  union and intersection are carefully defined in Wikipedia.
the interpretation here is that $ab$ divides any factors that appear in at least two the three numbers $n_1, n_2, n_3$. 
Conversely, if we start with a set $C \subset  (N_1 \cap N_3) \cup (N_2 \cap N_3) \cup(N_1 \cap N_2) $ can be split into 
$$ C = A \cup B  \text{ with } A \subseteq N_1 \cap N_3 \text{ and } A \cap B = \varnothing $$
Alternatively we could have split $C = G \cup D$ with $G \subset N_1 \cap N_2$\dots in fact for any function the product
$$ \Large
\sum_{a\mid \gcd(n_1,n_3)}\sum_{b\mid \gcd(n_2,\frac{n_1 n_3}{a^2})} F(ab)
=\sum_{g\mid \gcd(n_1,n_2)}\sum_{d\mid \gcd(\frac{n_1n_2}{g^2},n_3)} F(gd) $$

Older stuff
It should be enough to say $ab$ divides $M= \frac{\gcd(n_1, n_2)\gcd(n_2, n_3)\gcd(n_1, n_3)}{\gcd(n_1, n_2, n_3)^2}$ and specify a divisor $b|\gcd(n_2, M)$.  I found this formula with a Venn diagram.  It should be correct, and it should be clear that $g,d$ gives the same multiset.
A: 
Note: This answer is a supplement to the nice observation of @Elsa and a note to OPs comment regarding a proof based on induction.

We consider arithmetical functions which are complex-valued functions $f:\mathbb{N}\to \mathbb{C}$ defined on the set of positive integers.
An arithmetical function is called multiplicative if  $f(n)\ne 0$  for at least one integer $n$ and
$$f(mn)=f(m)f(n)\qquad\qquad   m,n\in \mathbb{N} \text{  with }(m,n)=1$$
An arithmetical function  is  called completely  multiplicative  if $f(n)\ne 0$  for at least one integer $n$ and
$$f(mn)=f(m)f(n)\qquad\qquad   m,n\in \mathbb{N}$$
Recall the (Dirichlet) convolution $f\star g$ of $f$ and $g$
$$(f\star g)(n)=\sum_{d|n}f(d)g\left(\frac{n}{d}\right)\qquad\qquad n\in \mathbb{N}$$

We now introduce according to the answer of Elsa a linear operator $T_n$ on $f$ foreach $n\in \mathbb{N}$ defined pointwise
$$(T_n(f))(N) := \sum_{a|(n,N)}f\left(\frac{Nn}{a}\right)\qquad\qquad N\in\mathbb{N}$$

According to Elsas elaboration OPs expression is equivalent with the following identity

\begin{align*}
T_{n_1}\circ T_{n_2}=\sum_{g|(n_1,n_2)}T_{\frac{n_1n_2}{g^2}}\tag{1}
\end{align*}

Now observe that $\sigma_k$ the sum of the $k$-th powers of (positive) divisors of $n$
$\sigma_k(n):=\sum_{d|n}\sigma^{k}$ fulfills according to Introduction to Arithmetical Functions by Paul J. McCarthy following relation
\begin{align*}
\sigma_k(n_1)\star \sigma_k(n_2)=\sum_{g|(n_1,n_2)}\sigma_k\left(\frac{n_1n_2}{g^2}\right)\cdot g^k\tag{2}
\end{align*}
which indicates the strong relationship with the expression (1) of the linear operator $T_n$.

The expression (2) is a so-called Busche-Ramanujan Identity.
It is fulfilled by multiplicative arithmetical functions which are the convolution of completely multiplicative functions. These functions are called specially multiplicative functions.

The following is valid for specially multiplicative functions (Theorem 1.12 by McCarthy):

If $f$ is a multiplicative function then the following statements are equivalent:
  
  
*
  
*$f$ is a convolution of two completely multiplicative functions
  
*There is a multiplicative function $F$ such that for all $m$ and $n$
$$f(mn)=\sum_{d|(m,n)}f\left(\frac{m}{d}\right)f\left(\frac{n}{d}\right)F(d)$$
  
  
*
  
*There is a completely multiplicative function $B$ such that for all $m$ and $n$
  
  
  $$f(m)f(n)=\sum_{d|(m,n)}f\left(\frac{mn}{d^2}\right)B(d)$$
  
  
*
  
*For all primes $p$ and all $\alpha \geq 1$,
  
  
  $$f\left(p^{\alpha+1}\right)=f\left(p\right)f\left(p^{\alpha}\right)
+f\left(p^{\alpha-1}\right)\left(f(p^2)-f(p)^2\right)$$

Note: Since OPs expression has been already successfully answered we know that (1) is true. Therefore we may conclude that the linear operator $T_n$ fulfills similar relations according to the theorem above.
