解説

各素因数ごとに 一般化Lucasの定理 で$\binom{N}{R}$ $mod$ $P_i ^ {E_i}$を計算し、中国剰余定理で解を復元します。
一般化Lucasの定理の詳細は参考文献をご覧ください。
計算量は前計算$O(M)$時間・空間、クエリ$O(logNlogM + logM)$時間になります。
内訳は以下の通りです。

• 前計算では、$P_i$の倍数を除いた階乗$x!_{P_i}$ $mod$ $P_i ^ {E_i}$を$P_i ^ {E_i}$まで計算すればよいです。
  計算量は時間・空間ともに$O(ΣP_i ^ {E_i}) \leq O(ΠP_i ^ {E_i}) = O(M)$です。
• クエリでは$N,R,N-R$を$P_i$進数表記に変換し、その後連続した$E_i$桁ごとに二項係数の積を計算します。
  時間計算量は前者が各$O(log_{P_i} N)$、後者は実装方針によりますが$O(log_{P_i} N)$や$O(E_i × log_{P_i} N)$になります。
  ただし$ΣE_i \leq logM$なので、どちらの実装であっても$O(ΣE_i × log_{P_i} N) \leq O(logNlogM)$時間で抑えられます。
• 中国剰余定理の復元は$O(Πlog P_i ^ {E_i}) = O(logM)$です。

なお本問の制約範囲外ですが、抽象化の際には$N,R = 0$や$N < R$、$M = 1$のケースに注意してください。

参考文献

一般化Lucasの定理と東大数学2021 - mathlog
https://mathlog.info/articles/1851

解答例(Python)

解答例では、二項係数が与えられるたびに階乗を計算しています。

一般化Lucasの定理 実装例(Python)

$O(√M + M)$時間・$O(M)$空間で階乗を前計算します。
二項係数クエリは1個あたり$O(logNlogM)$で計算します。

この実装例を用いた解答例を示します。