Assuming you already found the number of ways to distribute N elements… You should have some array A1, A2, A3, A4, …, AN with number of ways to distribute 1, 2, 3, 4, …, N elements respectively. So at this moment you only need to find a fast way to compute the multiplication of a range of this array.

A hard approach for that could be using some data structure like range-trees, or any equivalent technique like that (i.e. operate with sqrt intervals). But even this could be a quite slow for that number of cases.

So, an alternative way is using modular inverse, like in that solution which you linked. It uses cumulative multiplication (just the same as cumulative sum) to know in O(1) time, the multiplication between the first K <= N elements in the sequence. For example for K = 8 this array return the value of the multiplication between the first eight elements (A1*A2*A3*A4*A5*A6*A7*A8) … but, modulated by M. This is (A1*A2*A3*A4*A5*A6*A7*A8) % M.

If you want only the multiplication between A5 and A8, you need to divide previous result by (A1*A2*A3), this is (A1*A2*A3*A4*A5*A6*A7*A8) / (A1*A2*A3) or using M, ((A1*A2*A3*A4*A5*A6*A7*A8) % M) / ((A1*A2*A3) % M) using the same array A. This division will lead you to some wrong solutions due modulated values by M… for that reason you could use inverse modular and now rewrite the formula in this way: ((A1*A2*A3*A4*A5*A6*A7*A8) % M) / ((A1*A2*A3) % M) = ((A1*A2*A3*A4*A5*A6*A7*A8) % M) * Modular_Inverse((A1*A2*A3) % M)…

Greetings…