Uniform Chemistry
Contest: SWERC 2022-2023 - Online Mirror (Unrated, ICPC Rules, Teams Preferred)
Difficulty rating: 3200
Tags: dp, math

In a parallel universe there are $$$n$$$ chemical elements, numbered from $$$1$$$ to $$$n$$$. The element number $$$n$$$ has not been discovered so far, and its discovery would be a pinnacle of research and would bring the person who does it eternal fame and the so-called SWERC prize.

There are $$$m$$$ independent researchers, numbered from $$$1$$$ to $$$m$$$, that are trying to discover it. Currently, the $$$i$$$-th researcher has a sample of the element $$$s_i$$$. Every year, each researcher independently does one fusion experiment. In a fusion experiment, if the researcher currently has a sample of element $$$a$$$, they produce a sample of an element $$$b$$$ that is chosen uniformly at random between $$$a+1$$$ and $$$n$$$, and they lose the sample of element $$$a$$$. The elements discovered by different researchers or in different years are completely independent.

The first researcher to discover element $$$n$$$ will get the SWERC prize. If several researchers discover the element in the same year, they all get the prize. For each $$$i = 1, \, 2, \, \dots, \, m$$$, you need to compute the probability that the $$$i$$$-th researcher wins the prize.

Input format:
The first line contains two integers $$$n$$$ and $$$m$$$ ($$$2 \le n \le 10^{18}$$$, $$$1 \le m \le 100$$$) — the number of elements and the number of researchers.

The second line contains $$$m$$$ integers $$$s_1, \, s_2, \, \dots, \, s_m$$$ ($$$1 \le s_i < n$$$) — the elements that the researchers currently have.

Output format:
Print $$$m$$$ floating-point numbers. The $$$i$$$-th number should be the probability that the $$$i$$$-th researcher wins the SWERC prize. Your answer is accepted if each number differs from the correct number by at most $$$10^{-8}$$$.

Constraints:
Time limit: 2.0 s
Memory limit: 256.0 MB

Note:
In the first sample, all researchers will discover element $$$2$$$ in the first year and win the SWERC prize.

In the second sample, the last researcher will definitely discover element $$$3$$$ in the first year and win the SWERC prize. The first two researchers have a $$$50\%$$$ chance of discovering element $$$2$$$ and a $$$50\%$$$ chance of discovering element $$$3$$$, and only element $$$3$$$ will bring them the prize.

In the third sample, each researcher has an independent $$$50\%$$$ chance of discovering element $$$3$$$ in the first year, in which case they definitely win the SWERC prize. Additionally, if they all discover element $$$2$$$ in the first year, which is a $$$12.5\%$$$ chance, then they will all discover element $$$3$$$ in the second year and all win the prize.

Example input:
2 3
1 1 1

Example output:
1.0 1.0 1.0

Example input:
3 3
1 1 2

Example output:
0.5 0.5 1.0

Editorial:
Small constraints
The small version of this problem, which was used in the onsite competition, had constraints $$$n \le 100$$$, $$$m \le 10$$$.
The standard approach for this type of problem would be to use dynamic programming to compute the probability that each researcher wins the prize for each possible state of the process — the set of elements that each researcher has. However, even for the small constraints the number of such sets can be as high as $$$100^{10}$$$, which is clearly too big for us to process them all one by one.
Therefore the key idea is to use the independence of the processes that each researcher follows. Suppose we compute for each researcher $$$i$$$ and each year $$$j$$$ the probability $$$p_{ij}$$$ that the $$$i$$$-th researcher discovers element $$$n$$$ in the $$$j$$$-th year. In order for them to win the SWERC prize when doing so, the other researchers must have not discovered element $$$n$$$ yet. The probability that researcher $$$k$$$ has not yet discovered element $$$n$$$ by the $$$j$$$-th year can be computed as $$$1-\sum_{t<j}p_{kt}$$$, and because of the independence of different researchers, the probability that all other researchers have not yet discovered element $$$n$$$ is equal to $$$\prod_{k\ne i}(1-\sum_{t<j}p_{kt})$$$, and the probability that the $$$i$$$-th researcher wins the SWERC prize is equal to
$$$$$$\sum_j p_{ij}\prod_{k\ne i}\left(1-\sum_{t<j}p_{kt}\right)$$$$$$
Now we just need to compute the probability $$$p_{ij}$$$ that the $$$i$$$-th researcher discovers element $$$n$$$ in the $$$j$$$-th year for all $$$i$$$ and $$$j$$$. This can be done using dynamic programming that computes $$$q_{ij}$$$: the probability that a researcher that has element $$$n-i$$$ in the beginning (in other words, is $$$i$$$ elements away from the goal) discovers element $$$n$$$ in the $$$j$$$-th year. The probabilities that we want are then found as $$$p_{ij}=q_{n-s_i,j}$$$.
From the fusion experiment definition we directly obtain:
$$$$$$q_{ij}=\frac{1}{i}\sum_{0\le k<i} q_{k,j-1}$$$$$$
We can then apply this formula in increasing order of $$$i$$$, or in increasing order of $$$j$$$, to find all those probabilities, starting from $$$q_{00}=1$$$, $$$q_{i0}=0$$$ for $$$i > 0$$$, and then apply the formula that combines those values into our answer.
The dynamic programming has $$$O(n^2)$$$ states, and each state is processed in $$$O(n)$$$, so its running time is $$$O(n^3)$$$. The final formula is computed in $$$O(n^2 m)$$$ for each of the $$$m$$$ researchers, so the overall running time is $$$O(n^2(n+m^2))$$$, which is fast enough for $$$n \le 100$$$, $$$m \le 10$$$.
Large constraints
The above solution can be optimized to run in $$$O(n(n+m))$$$ using relatively standard techniques, but that still won't be enough to solve the large version that was given in the mirror round: $$$n \le 10^{18}$$$, $$$m \le 100$$$. Here we need to come up with more radical improvements.
The first big improvement is to notice that the distance to $$$n$$$ decreases twice every year on average, and therefore the typical number of years needed to reach $$$n$$$ is likely to be proportional to $$$\log n$$$. Moreover, for each $$$\varepsilon$$$ there likely exists a constant $$$C$$$ such that the number of years is below $$$C \log n$$$ with probability $$$1-\varepsilon$$$.
To prove this formally, we can notice that each year the distance to $$$n$$$ decreases at least twice with probability of at least $$$\frac{1}{2}$$$, so after $$$k$$$ years the distance to $$$n$$$ decreases at least twice with probability at least $$$1-\frac{1}{2^k}$$$. We can choose $$$k$$$ such that $$$\frac{1}{2^k}<\frac{\varepsilon}{log_2 n}$$$, and then the number of years will not exceed $$$k \log_2 n$$$ with probability at least $$$1-\varepsilon$$$. This is a very simple estimate, using Bernstein inequalities one can prove much more.
But we don't really need a formal proof during the round, instead the intuition about the distance decreasing twice on average can be coupled with an experiment that measures how quickly the numbers $$$q_{ij}$$$ described above become very small. In practice, it turns out that we can always assume $$$j \le 100$$$ and get the required precision.
This observation improves the running time to $$$O((n+m)\log n)$$$, which is still not good enough. We will describe two independent ways to proceed further.
The exact approach. Let us examine the recurrence for $$$q_{ij}$$$ more closely. If we apply it repeatedly, we end up having the following sum:
$$$$$$ q_{ij} = \sum_{0<a_1<\dots<a_j=i}\prod_{k=1}^{j}\frac{1}{a_k} $$$$$$
This prompts us to consider the following polynomial:
$$$$$$ f(x) = \frac{1}{i}\prod_{k=1}^{i-1}\left(1+\frac{1}{k}x\right) $$$$$$
The coefficient next to $$$x^{j-1}$$$ in this polynomial is equal to $$$q_{ij}$$$, because to get it we can use any of the ways to pick $$$j-1$$$ terms of the product from which we take the part with $$$x$$$ and then we get a product of the reciprocals of their indices, which is precisely what the sum for $$$q_{ij}$$$ has.
So now we need to compute the first $$$O(\log n)$$$ coefficients of this polynomial quickly. Let us take its natural logarithm (in the formal power series sense), and substitute the Taylor series for $$$\log(1+y)$$$:
$$$$$$ \log f(x) = -\log i + \sum_{k=1}^{i-1}\log\left(1+\frac{1}{k}x\right) = -\log i - \sum_{k=1}^{i-1} \sum_{t=1}^{\infty} \frac{(-1)^t}{tk^t}x^t = -\log i - \sum_{t=1}^{\infty} \frac{(-1)^t}{t}x^t \sum_{k=1}^{i-1} \frac{1}{k^t} $$$$$$
Since the first $$$O(\log n)$$$ coefficients of $$$f(x)$$$ can be found by exponentation of the formal power series and depend only on the first $$$O(\log n)$$$ coefficients of $$$\log f(x)$$$, all we need now is a way to find the sums of inverse powers $$$\sum_{k=1}^{i-1} \frac{1}{k^t}$$$ for $$$t$$$ up to $$$C \log n$$$.
This can be done using the Euler-Maclaurin formula, or using the following trick which is likely equivalent to it in some sense. We denote $$$S_t(a,b)=\sum_{k=a}^{b}k^{-t}$$$, and notice that using the Taylor series for $$$\log(1+y)$$$ we get:
$$$$$$ S_1(a,b)=\sum_{k=a}^{b}\log\left(1+\frac{1}{k}\right)+\sum_{k=a}^{b}\left(\frac{1}{k}-\log\left(1+\frac{1}{k}\right)\right)=\log(b+1) - \log a + \sum_{k=a}^{b}\sum_{p=2}^{\infty}\frac{(-1)^p}{pk^p}= $$$$$$ $$$$$$ =\log(b+1) - \log a+\sum_{p=2}^{\infty}\frac{(-1)^p}{p} S_p(a,b) $$$$$$
Similarly for $$$t>1$$$ using the Taylor series for $$$(1+y)^{-(t-1)}$$$ we get:
$$$$$$ S_t(a,b)=\sum_{k=a}^{b}\frac{1}{t-1}\left(\frac{1}{k^{t-1}}-\frac{1}{(k+1)^{t-1}}\right)+\sum_{k=a}^{b}\left(\frac{1}{k^t}-\frac{1}{t-1}\left(\frac{1}{k^{t-1}}-\frac{1}{(k+1)^{t-1}}\right)\right)= $$$$$$ $$$$$$ =\frac{1}{t-1}\left(\frac{1}{a^{t-1}} - \frac{1}{(b+1)^{t-1}}\right)+\sum_{k=a}^{b}\left(\frac{1}{k^t}-\frac{1}{(t-1)k^{t-1}}\left(1-\frac{1}{(1+\frac{1}{k})^{t-1}}\right)\right)= $$$$$$ $$$$$$ =\frac{1}{t-1}\left(\frac{1}{a^{t-1}} - \frac{1}{(b+1)^{t-1}}\right)+\sum_{k=a}^{b} \frac{1}{(t-1)k^{t-1}}\sum_{p=2}^{\infty}(-1)^p {t-2+p \choose p} \frac{1}{k^p}= $$$$$$ $$$$$$ =\frac{1}{t-1}\left(\frac{1}{a^{t-1}} - \frac{1}{(b+1)^{t-1}}\right)+ \frac{1}{t-1}\sum_{p=t+1}^{\infty}(-1)^{p-t+1} {p-1 \choose p-t+1} S_p(a,b) $$$$$$
These formulas express $$$S_t(a,b)$$$ using $$$S_p(a,b)$$$ where $$$p>t$$$. Since for sufficiently large $$$t$$$ $$$S_t(a,b)$$$ can be computed with good enough precision directly because $$$\frac{1}{k^t}$$$ decreases very quickly, this allows to compute $$$S_t(a,b)$$$ fast for all $$$t$$$ up to $$$C \log n$$$, and therefore to obtain the polynomial we need and to solve the problem.
The approximate approach. Alternatively, we can approximate our problem with its continuous version: suppose the result of the fusion experiment was not an integer, but a uniformly distributed real number, and we were interested in the probability to get a number less than 1 after exactly $$$j$$$ years if we start with $$$i$$$.
This yields the following formulas:
$$$$$$ q'_{i1}=\frac{1}{i} $$$$$$
$$$$$$ q'_{ij}=\frac{1}{i}\int_{1}^{i} q'_{k,j-1} dk $$$$$$
The integral is actually not very hard to compute, and we get
$$$$$$ q'_{ij}=\frac{\log^{j-1} i}{i(j-1)!} $$$$$$
This approach does not give good enough approximation by itself, but we can improve on it in the following way: suppose we want to find some number $$$q_{ij}$$$. Let us choose some boundary $$$b$$$, for example $$$b=10^6$$$, and for all values of $$$i < b$$$ we can just compute $$$q_{ij}$$$ using the original recurrence. When $$$i$$$ is larger than or equal to $$$b$$$, we will separate the process into two stages: while $$$i$$$ is larger than or equal to $$$b$$$, and when it becomes smaller, and we will use the continuous approximation only for the first stage.
Getting from $$$i$$$ to $$$b$$$ in the continuous setting is the same as getting from $$$i/b$$$ to $$$1$$$, so we can use the above formula for $$$q'$$$ for this.
To glue the two stages together, we can notice that whenever $$$i$$$ becomes less than $$$b$$$, because of the uniform choice all values of $$$i$$$ between 0 and $$$b-1$$$ are equally likely, and therefore this situation is the same as if we had $$$i=b$$$ on the previous step. Therefore we can use the following formula to compute our approximate answer:
$$$$$$ q_{ij} \approx \sum_{k=1}^{j} q'_{\frac{i}{b},k} q_{b,j-k+1} $$$$$$
This (barely) gives good enough precision for $$$b=10^6$$$.
Precision issues. The solutions to this problem, especially the exact approach involving formal power series exponentiation and Taylor expansions, can easily suffer from floating-point precision issues. We are very confident that our reference solutions are correct because we have several different approaches, implemented in different languages, using different floating-point types, all agree on the answers.
To achieve this, we had to make some changes that we would like to mention to make it easier to follow in our footsteps. In the exact solution, we used the above formulas to compute $$$S_t(a,b)$$$ only for sufficiently large values of $$$a$$$ (at least $$$a \ge 2$$$), and added the first few terms directly, as that allowed for better convergence. We used those formulas only for small values of $$$t$$$ ($$$t \le 5$$$, or even just $$$t \le 2$$$), as otherwise the large binomial coefficients have amplified the precision errors too much. When computing the exponential of a formal power series using an $$$O(n\log n)$$$ approach (which was not necessary in this problem), we found that using the identity $$$exp(f(x))=(exp(\frac{f(x)}{8}))^8$$$ can improve the precision dramatically.
In addition, we noticed that if we use the above formulas to compute $$$S_t(a,b)$$$ only for $$$a \ge 10^6$$$, then we need them only for $$$t=1$$$ and $$$t=2$$$, and even for those we need only a couple of terms, making this part of the solution really short:
$$$$$$ S_1(a,b) \approx \log(b+1)-\log a+\frac{1}{2}S_2(a,b) $$$$$$ $$$$$$ S_2(a,b) \approx \frac{1}{a} - \frac{1}{b+1} $$$$$$