(より詳しい解説は後日ブログで公開します)

非負整数 \(u, v, w\) によって \(BC u + CA v + AB w\) と書ける非負整数の集合を \(I\) とすると，\(S = \displaystyle\prod_p \sum_{i\in I} p^{-2i}\) となります (\(p\) は素数をわたる．Euler 積)．

\(I\) の母関数 \(F(x) = \sum_{i\in I} x^i\) を考えます．さらに，\(i = BC u + CA v + AB w\) という表し方を数える母関数 \(G(x) = \dfrac{1}{(1-x^{BC}) (1 - x^{CA}) (1 - x^{AB})}\) も考えます．\(i \bmod {ABC}\) を決めると \(u \bmod A, v \bmod B, w \bmod C\) が決まるので，表し方の自由度は \(BC (u \bmod A) + CA (v \bmod B) + AB (w \bmod C)\) に対してどこに \(ABC\) を足していくかというぶんだけあります．よって \(F(x), G(x)\) を適当なところから \(ABC\) 項おきに見るとそれぞれ \(0, \ldots, 0, 1, 1, 1, 1, \ldots\) と \(0, \ldots, 0, 1, 3, 6, 10, \ldots\) のようになるので，\(G(x) = F(x) \dfrac{1}{(1 - x^{ABC})^2}\) がわかり，\(F(x) = \dfrac{(1 - x^{ABC})^2}{(1-x^{BC}) (1 - x^{CA}) (1 - x^{AB})}\) です．

\(F(x)\) の形から，Riemann ゼータ関数を用いて \(S = \dfrac{\zeta(2BC) \zeta(2CA) \zeta(2AB)}{\zeta(2ABC)^2}\) と計算できます．これであとは Bernoulli 数 \(B_n = [x^n/n!] \dfrac{x}{\mathrm{e}^x-1}\) が計算できればよいです．

\(B_n = [x^n/n!] \dfrac{\log(1+(\mathrm{e}^x-1))}{(\mathrm{e}^x-1)} = [x^n/n!] \displaystyle\sum_{k\ge0} \dfrac{(-1)^k}{k+1} (\mathrm{e}^x-1)^k = [x^n/n!] \displaystyle\sum_{k=0}^n \dfrac{(-1)^k}{k+1} (\mathrm{e}^x-1)^k\) です．これを \(\mathrm{e}^x\) の多項式として展開したときの係数列が得られればよいです．\(\dfrac{\log(1+y)}{y}\) の満たす微分方程式を用いると漸化式が得られます．詳しくは EntropyIncreaser さんの記事を参照してください．\(0^n, \ldots, n^n\) の計算に線型時間の素数篩を用いることで， \(O(n)\) 時間で \(B_n\) を求められます．