We will ignore a few corner cases, such as $$$i=j$$$ or $$$j=k$$$, and other corner cases involving $$$0$$$. The contributions of those cases are relatively simple or classical, so we omit them here; Consequently, some summation bounds below are not written with full formal tightness.
We actually want to compute
$$$ \sum\limits_{i=0}^m\sum\limits_{j=0}^m\sum\limits_{k=0}^m[\gcd(i\oplus j,j\oplus k)|n] $$$As usual, enumerate the value of the $$$\gcd$$$. Thus
$$$ \begin{aligned} =&\sum\limits_{d|n} \sum\limits_{i=0}^m\sum\limits_{j=0}^m\sum\limits_{k=0}^m[\gcd(i\oplus j,j\oplus k)=d]\\ =&\sum\limits_{d|n}\sum\limits_{p=i\oplus j}\sum\limits_{q=j\oplus k}[\gcd(p,q)=d]c(p,q) \end{aligned} $$$where $$$p=i\oplus j$$$, $$$q=j\oplus k$$$, and $$$c(p,q)$$$ is the number of tuples $$$(i,j,k)$$$ with those values of $$$p$$$ and $$$q$$$.
Apply Möbius inversion to the condition $$$\gcd(p,q)=d$$$:
$$$ \begin{aligned} =&\sum\limits_{d|n} \sum\limits_{d|p}\sum\limits_{d|q}[\gcd(\frac{p}{d},\frac{q}{d})=1]c(p,q)\\ =&\sum\limits_{d|n} \sum\limits_{d|p}\sum\limits_{d|q}\sum\limits_{d'|\frac{p}{d},d'|\frac{q}{d}}\mu(d')c(p,q)\\ \end{aligned} $$$where $$$\mu$$$ is the Möbius function.
We now change the summation order and pull $$$d'$$$ outward.
$$$ \begin{aligned} =&\sum\limits_{d|n}\sum\limits_{d'}^{\frac{2m}{d}}\mu(d')\sum\limits_{p}\sum\limits_{q}[dd'|p\land dd'|q]c(p,q)\\ \end{aligned} $$$Writing $$$T=dd'$$$ we obtain
$$$ \begin{aligned} =&\sum\limits_{T=1}^{2m}\sum\limits_{d|\gcd(n,T)}\mu(\frac{T}{d})\sum\limits_{p}\sum\limits_{q}[T|p\land T|q]c(p,q)\\ \end{aligned} $$$The outer sum over $$$T$$$ runs up to $$$2m$$$ because $$$p,q\in[0,2m]$$$.
The inner quantity depends on $$$T$$$ only. Let $$$f(T)$$$ be the number of tuples $$$(i,j,k)$$$ such that $$$T$$$ is a divisor of both $$$i\oplus j$$$ and $$$j\oplus k$$$.
Then the entire sum becomes
$$$ \begin{aligned} =&\sum\limits_{T=1}^{2m}\sum\limits_{d|\gcd(n,T)}\mu(\frac{T}{d})f(T)\\ \end{aligned} $$$Due to the series sum of the Harmonic series, we can enumerate $$$T$$$ and $$$d$$$ in $$$\mathcal O(m\ln m)$$$. Hence, the main task is to compute $$$f(T)$$$ efficiently for all $$$T$$$.
Fix $$$T$$$. We enumerate $$$j$$$ and count the number of pairs $$$(i,k)$$$ such that $$$i\oplus j$$$ and $$$j\oplus k$$$ are multiples of $$$T$$$.
Equivalently, we count the number of pairs $$$(d_1,d_2)$$$ under the enumerated $$$j$$$ such that
$$$ i\oplus j=d_1T, j\oplus k=d_2T $$$with the additional constraints $$$d_1T\oplus j\le m$$$ and $$$d_2T\oplus j\le m$$$ (because $$$(i,k)$$$ must lie in $$$[0,m]$$$).
Equivalently, if we can enumerate $$$d_1,d_2$$$ to get $$$d_1T$$$ and $$$d_2T$$$, we only need to find the number of $$$j$$$ that satisfy both constraints. If we enumerate $$$d_1$$$ to find the set of $$$j$$$ and record it, we can enumerate $$$d_2$$$ similarly, and what we need to calculate is the size of the intersection of the two sets.
Consider enumerating $$$d_1$$$ to get $$$d_1T$$$, due to bitmask theory, the $$$j$$$ such that $$$d_1T\oplus j\le m$$$ will be distributed in $$$\mathcal O(\log m)$$$ intervals.
By using Segment Tree or Fenwick Tree naively, we will solve the problem in $$$\mathcal O(m\ln m\log^2 m)$$$, which is too slow to get accepted in the problem.
By LiHn's observation, actually, we can apply a Trie instead of a Segment Tree to record the intervals, because each interval can be formed by a subtree in the Trie. This will allow us to avoid incurring a logarithmic cost to find intervals in data structures every time.
By Leftist_G's observation, there are only two types of intervals, but one type of intervals has $$$\mathcal O(\log)$$$, and another has $$$\mathcal O(1)$$$. We can avoid using Segment Tree or Fenwick Tree in the former.
By liuhangxin's observation, we can avoid using any data structure. We can naively save the intervals, sort them, and then apply two pointers to calculate the intersection. Note that if we implement the method naively, the complexity will turn to $$$\mathcal O(m\ln m\log m\log\log m)$$$ (although it may be enough to get Accepted). He used a designed order and implmented carefully to ensure his complexity.
All solutions above can solve the problem in $$$\mathcal O(m\ln m\log m)$$$.
Note that some not-so-good solutions may apply __int128 to avoid integer overflow. Don't forget to deal with corner cases.
