assert(false)

[H. S. Teoh]

On Fri, Jun 21, 2013 at 12:36:13AM +0100, Joseph Rushton Wakeling wrote:
[...]
>     private size_t diceImpl(Rng, Range)(ref Rng rng, Range proportions)
>     if (isForwardRange!Range && isNumeric!(ElementType!Range) && isForwardRange!Rng)
>     {
>         double sum = reduce!("(assert(b >= 0), a + b)")(0.0, proportions.save);
>         enforce(sum > 0, "Proportions in a dice cannot sum to zero");
>         immutable point = uniform(0.0, sum, rng);
>         assert(point < sum);
>         auto mass = 0.0;
> 
>         size_t i = 0;
>         foreach (e; proportions)
>         {
>             mass += e;
>             if (point < mass) return i;
>             i++;
>         }
>         // this point should not be reached
>         assert(false);
>     }
[...]

On Fri, Jun 21, 2013 at 09:45:55AM -0700, Ali Çehreli wrote:
> On 06/21/2013 02:56 AM, Joseph Rushton Wakeling wrote:
> 
> > I did also have to think about it for a few seconds.  You have
> > a sum of values from an array; you generate a
> > uniformly-distributed random number point in [0.0, sum); you
> > sequentially sum values from the same array, and return when
> > that partial sum ("mass") exceeds the value of point.  Since
> > point < sum you must do this eventually
> 
> I spent a considerable amount of time reading that function as well. :)
> 
> I was afraid whether two separate floating point sums would give the
> same result. After all, the initial sum is a 'double' but the type
> of the elements of 'proportions' may be different. (I don't know the
> answer.)
[...]

Hmm. This is a very good point. I actually found two potential problems
with this code, due to the peculiarities of floating-point arithmetic:

If 'proportions' contains an NaN, then sum will also be nan, and the
enforce line will fail (and with a misleading message).

If 'proportions' contains +inf, then sum will also be +inf, so the
enforce will pass. Now what does uniform() return if sum is +inf? I
tested this, and found that uniform(0.0, float.infinity) consistently
returns double.infinity, and uniform(0.0, real.infinity) consistently
returns real.infinity. Both will cause the next assert() will fail. (Had
this assert not been there, the foreach could terminate without
returning, e.g. if the last element of 'proportions' is +inf, then
point=+inf and point < mass is never true.)

There's also a potential pitfall if 'proportions' is a range of floats,
with some adjacent elements that differ too much in magnitude such that
summing them as floats will give a different answer from summing them as
doubles. Fortunately, reduce is called with 0.0 as seed value, which is
a double, so this problem doesn't occur here (each element is promoted
to double before summation, so no loss of accuracy occurs).

The reverse problem of a range of reals on x86, where sum may lose
precision from very small real values that don't fit in a double,
doesn't cause any problem either, because mass is also a double, and any
precision loss is duplicated in both the reduce and the foreach, so
there's no possibility of the loop terminating without returning.

So, long story short, the assert(false) actually will never be reached.
But I can't see the compiler perform this level of analysis (involving
floating-point analysis on hypothetical template arguments!) to reach
that conclusion! So yeah, the assert(false) should be left there. :)

(Don't you just love floating-point? They make life so much more
interesting than if we'd been able to use mathematically-exact reals!
:-P)


T

-- 
Written on the window of a clothing store: No shirt, no shoes, no service.