\(\text{MAX}=5\times10^5\) とします。

\(x=0,1,\ldots,\text{MAX}\) に対する答えをまとめて計算することを考えます。

まず、ある閾値 \(X\) を用意し、以下のように \(2\) 種類の解法を用意します。

[1] \(B_i\le X\) のとき

\(d[a][b]\) を「\(A_i=a,B_i=b\) を満たす \(i\) の個数」としたカウンタを用意し、それを計算します。その後、各 \(b\) に対し \(d\) の情報を元に、\(B_i=b\) を満たす \(i\) に対する \((x+A_i)\bmod B_i\) の値の総和をまとめて計算します。

具体的には、\(\displaystyle c_b=\sum_{a=0}^{b-1}d[a][b]\) とし、\(\displaystyle g_b(r)=\sum_{a=0}^{b-1}d[a][b]\bigl((r+a)\bmod b\bigr)\) とおきます。すると、\(\displaystyle g_b(0)=\sum_{a=0}^{b-1}a\,d[a][b]\) であり、\(\displaystyle g_b(r+1)=g_b(r)+c_b-b\times d[b-1-r][b]\ (0\le r<b-1)\) が成り立つため、\(g_b(0),g_b(1),\ldots,g_b(b-1)\) を \(O(b)\) 時間で計算することができます。

各 \(x\) に対する寄与は \(g_b(x\bmod b)\) なので、これを各 \(x\) の答えに加えます。

各 \(b\) に対する計算量は \(O(\text{MAX}+b)\) であるため、合計して \(O(X\times\text{MAX}+X^2)\) となります。

[2] \(B_i>X\) のとき

\(\displaystyle (x+A_i)\bmod B_i=x+A_i-B_i\left\lfloor\frac{x+A_i}{B_i}\right\rfloor\) を用います。

\(x+A_i\) の部分は初期値と傾きとしてまとめて扱うことができます。また、\(\displaystyle \left\lfloor\frac{x+A_i}{B_i}\right\rfloor\) の値が変わる箇所は \(\displaystyle O\left(\frac{\text{MAX}}{B_i}\right)\) 箇所であるため、変化する箇所を記録し、まとめて imos 法で答えを計算することができます。

変化する箇所を記録する計算量は各 \(i\) に対して \(\displaystyle O\left(\frac{\text{MAX}}{B_i}\right)\) 時間となり、まとめて \(\displaystyle O\left(N\times\frac{\text{MAX}}{X}\right)\) 時間と評価することができます。

最後に差分配列の累積和を計算する \(O(\text{MAX})\) 時間を合わせると、計算量は \(\displaystyle O\left(\text{MAX}+N\times\frac{\text{MAX}}{X}\right)\) となります。

[1], [2] を合わせ、入力・出力の計算量も含めることで、全体の計算量は \(\displaystyle O\left(N+Q+X\times\text{MAX}+X^2+N\times\frac{\text{MAX}}{X}\right)\) となります。

\(\displaystyle X=\left\lceil\sqrt{\text{MAX}}\right\rceil\) とすることで、この計算量は \(\displaystyle O\left((N+\text{MAX})\sqrt{\text{MAX}}+Q\right)\) と評価することができ、この解法は高速に動作します。

以上を適切に実装することでこの問題に正答することができます。

実装例(Python3)