公式解説の \(s_n, t_n, u_n\) の間に成り立つ関係式は、約数・倍数ゼータ変換の形になっています。そのため、\(O(N + M\log\log M)\) 時間で計算できます。計算方法は noshi91 さんの記事 などが参考になります。

約数・倍数ゼータ変換は和で計算することが多いですが、一般には可換モノイドが載ります。\(u_n\) の計算は可換モノイド \((\mathbb{Z}, \max)\) 上での約数ゼータ変換として行えます。