Re: nchoosek warning, checks and tests

[Francesco Potortì]

>>>>   if (n == 1)
>>>> -    if (k > v/2)
>>>> -      k = v - k;
>>>> +    k = min (k, v-k);          # improve precision at next step
>>>> +    A = round (exp (sum (log (v-k+1:v)) - sum (log (2:k))));
>>>> +    if (A > 1/(exp (2*k*eps) - 1))
>>>> +      warning ("nchoosek", "nchoosek: possible loss of precision");
>>>>     endif
>>>> -    A = round (prod ((v-k+1:v)./(1:k)));
>>
>>>Why did you replace
>>>round (prod ((v-k+1:v)./(1:k)))
>>>by
>>>round (exp (sum (log (v-k+1:v)) - sum (log (2:k))))
>>>?

>> When comparing nchoosek to bincoeff, it turns out that the first one is
>> thought for integer and not very big numbers, while the second is
>> tailored for fractional and big numbers.  In fact, bincoeff uses the
>> gamma function, which allows real arguments, and is faster for big
>> numbers, because nchoosek creates an array to be summed or multiplied.
>>
>> So I think it is important for nchoosek to emit a warning when the
>> result is not exact (Matlab does the same).  However, as far as I can
>> tell, the error bounds for the two expressions above (the prod and the
>> sum of logs) is the same, so using one or the other is the same, as long
>> as we want to have an exact result.  On the other hand, the sum of logs
>> avoids overflow in some cases, such as v=300, k=140.
>
>No, the existing expression doesn't overflow for these numbers. 

I stand corrected.  Sorry for not being scrupolous enough, and thank you
for spotting this error of mine.

>In fact it simply cannot overflow if the result doesn't overflow.

Yes, I see.

>> Important: this assumes that the log and exp functions produce a
>> relative error near eps, which I am not sure is true in Octave's
>> implementation.  If this assumption is wrong, then the prod
>> implementation is preferable.
>
>Whether true or not, the exp-log implementation calculates 2*k
>logarithms and one exponential, which is certainly slower, not
>speaking about more room for errors to accumulate. On the other hand,
>my expression does k divisions and multiplications.

Yes, your implementation is preferable.

>> Some more considerations:
>>
>> rethinking about it, the above expression
>>  if (A > 1/(exp (2*k*eps) - 1))
>> is equivalent with the much simpler
>>  if (A*2*k*eps > 1)
>> moreover, this is not entirely accurate and should be instead:
>>  if (A*2*k*eps > 0.5)
>>
>> Note that Matlab is not so picky about rounding errors, and as long as
>> its help page is accurate, gives a warning only when (A > 1e15).

>Actually, I don't exactly get the reason for the warning, because just
>about every computation in Octave (and Matlab) is inexact. That's just
>a property of floating-point numbers. Surely we don't warn at "sin(1)"
>that the result is only accurate to 15-16 decimal places...

Yes.  However, nchoosek is a function that works with integer arguments.
Ideally nchoosek should do integer arithmetic, in my opinion.  However,
we can resort to floating point to get more significant digits.
Apparently, this is the rationale behind Matlab's warning, and I think
it is reasonable.  On the other hand bincoeff is the right way to go in
other cases.

>but if it's in Matlab, then why not. Maybe it's for naive users who
>could expect the result to be computed exactly.

I wish to correct my latest changeset as you pointed out (use prod
rather than sum of logs) and as I discovered while writing (correct the
test that triggers the warning).  Should I create one more changeset to
be chained with the one I just posted?

I would create an alternative changeset with respect to the current
version in svannah, but I don't understand enough of Mercurial to do
that.  Please let me know what to do.

-- 
Francesco Potortì (ricercatore)        Voice: +39 050 315 3058 (op.2111)
ISTI - Area della ricerca CNR          Fax:   +39 050 315 2040
via G. Moruzzi 1, I-56124 Pisa         Email: address@hidden
(entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/