[libc++] Bugfix to std::binomial_distribution<int>

[atmnpatel]

Summary

The current implementation of binomial_distribution is not guaranteed to converge for certain extreme configurations of the engine and distribution. This is due to a mistake in the implementation of the algorithm from the given reference paper. The algorithm in the paper is guaranteed to terminate but has redundant statements. The current implementation simplified away the redundancy into a while loop, but it excludes the return condition of the case where a good sample cannot be returned for the particular sample being used from the uniform distribution, which is what causes the infinite loop. This change guarantees termination by recognizing that a good sample cannot be returned and returning 0 after breaking the loop. This is also in contrast to the paper because the return value as specified in the paper violates basic checks in at least a subset of the extreme cases where the current implementation fails to terminate. This default return value of 0 is satisfactory for the extreme case known so far.

Since this is only meant to affect extreme cases where the algorithm does not terminate anyways, the behavior is expected to remain exactly the same for all non-extreme cases that have been terminating so far.

Fixes: https://bugs.llvm.org/show_bug.cgi?id=44847

[ldionne]

Is there any way to add a test for this? Perhaps something based on the reproducer in https://llvm.org/PR44847?

[atmnpatel]

Yep, sorry I forgot to add it before. It's a bit of a wordy and entirely based on the bug report.

[ldionne]

Thanks for your contribution!

What email address do you want this committed with? In the future, please consider using arc diff to upload patches, since that automatically handles attribution.

Also, I'll admit I am thorn apart on this change. One the one hand I can confirm it fixes the issue reported in PR44847, however I am not sure of what it does to the algorithm itself. Indeed, our tests for <random> are not exactly great, and I don't understand or have access to the algorithm that's implemented in binomial_distribution. So I can't say for sure that it doesn't have other side effects such as making binomial_distribution not binomial anymore, or similar. This is not a great place to be in.

Perhaps you can walk me through your thought process for deciding when the algorithm should break? @EricWF @mclow.lists have either of you fixed bugs like this in <random>? How do you do it?

[atmnpatel]

Sorry, I was transitioning between machines while updating the commit and I couldn't figure out arc diff - I'd like it to be committed under [e-mail].

I have access to the paper and can confirm that although the implementation is a convoluted implementation of the algorithm in the paper, it does indeed follow the algorithm as best as it should - the most direct implementation of the algorithm is at least partially incorrect (I implemented it and it doesn't pass the most basic tests we have as far as I can tell). The original algorithm consists of a few for-loops (with return conditions) as well as a catch-all return condition, which as I understand it, has been merged rather cleverly into a single while loop while maintaining its behavior (the original algorithm as written is almost certainly erroneous probably due to typographical mistakes). This edit should cause the while loop to terminate when it has looped through an empty loop - which should correspond to the case in the original algorithm where it terminates when it is known that a good sample cannot be constructed given the initial sample from uniform distribution, and returns a fairly default value.

If we assume that the current implementation is in fact a functional implementation of the binomial_distribution, this edit will not change that fact, and will only return 0 in the cases where the current implementation will not terminate. Also I'd like to be really upfront about this - It returns 0 as opposed to following the paper in the catch-all return condition because if the paper is followed, the mean is nowhere near it should be, and this is particular evident in the case of the bug when the probability is set to be extremely small. If instead we return 0, we do end up with a better approximation to the binomial distribution in terms of mean, kurtosis, etc. Since this bug is really hard to reproduce, I can't speak to how generically this will remain true but it is empirically true given what we know now. I did find newer papers with authors that I can contact in case of confusion on the same topic of sampling from a binomial distribution, but given that we have an implementation that is known to work, I wasn't sure if it was appropriate to change the sampling method entirely while claiming to fix a bug that is clearly caused by a simple case of missing a termination guarantee. That is, I would be open to implementing a newer method for binomial_distribution if this seems like an unsatisfactory resolution.

[EricWF]

Could you please upload this diff with more context?
(https://llvm.org/docs/Phabricator.html#requesting-a-review-via-the-web-interface)

[atmnpatel]

Better? Or did I misunderstand the request? The updated patch paragraph from the Code Reviews page seems to be primarily for changes in the content of the patch rather than changes in the commit message/context provided.

[ldionne]

In D74997#1918498, @atmnpatel wrote:

Sorry, I was transitioning between machines while updating the commit and I couldn't figure out arc diff - I'd like it to be committed under [e-mail].

Ok, will do.

[...]

If we assume that the current implementation is in fact a functional implementation of the binomial_distribution, this edit will not change that fact, and will only return 0 in the cases where the current implementation will not terminate. Also I'd like to be really upfront about this - It returns 0 as opposed to following the paper in the catch-all return condition because if the paper is followed, the mean is nowhere near it should be, and this is particular evident in the case of the bug when the probability is set to be extremely small. [...]

Ok, reading your explanation, I believe I follow. Thanks for explaining. I'm fine with the change.

In D74997#1918581, @EricWF wrote:

Could you please upload this diff with more context?

@atmnpatel What Eric meant here is to upload a patch with some Diff context. For example, see https://reviews.llvm.org/D76288 where you can expand/collapse the code surrounding the changes, and where it says instead Context not available. We prefer reviewing patches with Diff context cause it's easier, but this time around it's okay.

Thanks a lot for your contribution, I'm going to commit this for you now under [e-mail].

libcxx/include/random

binomial_distribution<_IntType>::operator()(_URNG& __g, const param_type& __pr)
     if (__u < 0)
         return __pr.__r0_;
     double __pu = __pr.__pr_;
     double __pd = __pu;
     result_type __ru = __pr.__r0_;
     result_type __rd = __ru;
     while (true)
     {
+        bool __break = true;
         if (__rd >= 1)
         {
             __pd *= __rd / (__pr.__odds_ratio_ * (__pr.__t_ - __rd + 1));
             __u -= __pd;
+            __break = false;
             if (__u < 0)
                 return __rd - 1;
         }
         if ( __rd != 0 )
             --__rd;
         ++__ru;
         if (__ru <= __pr.__t_)
         {
             __pu *= (__pr.__t_ - __ru + 1) * __pr.__odds_ratio_ / __ru;
             __u -= __pu;
+            __break = false;
             if (__u < 0)
                 return __ru;
         }
+        if (__break)
+            return 0;
     }
 }
 
 template <class _CharT, class _Traits, class _IntType>
 basic_ostream<_CharT, _Traits>&
 operator<<(basic_ostream<_CharT, _Traits>& __os,
            const binomial_distribution<_IntType>& __x)
 {