(Free42) not a bug - COMB

[Werner]

.. but an easy fix for a tiny 'problem'.
COMB obviously computes the result as

n/1*(n-1)/2*(n-2)/3...*(n-p)/p

But it will give an Out of Range error for eg C(20408,10204) while the result is 1.47e6141, well within range. That can easily be fixed computing COMB as

1E-5*n/1*(n-1)/2*(n-2)/3...*(n-p)/p*1E5

This will allow up to COMB(20420,10210), and 10210 is the largest value p for which C(n,p) doesn't overflow but C(n,p)*p does (the last value in the chain calculation before the final division by p), so 1E-5 is sufficiently small. 1E-4 is not ;-)

Cheers, Werner

[Paul Dale]

This looks like it would only address the specific example. Other cases would still overflow.

The 34S uses logΓ with a forced rounding to integer: \( ^nC_r = e^{logΓ(n) - logΓ(r) - logΓ(n - r)} \)

This risks cancellation of digits by subtraction.


Pauli

[Albert Chan]

We can shift the argument: \( \binom{n}{r} = {n \over n-r}\; \binom{n-1}{r} \)

→ C(20408, 10204) = 2 × C(20407, 10204)

[Werner]

(07-13-2020 12:03 PM)Paul Dale Wrote:  This looks like it would only address the specific example. Other cases would still overflow.

The 34S uses logΓ with a forced rounding to integer: \( ^nC_r = e^{logΓ(n) - logΓ(r) - logΓ(n - r)} \)

This risks cancellation of digits by subtraction.


Pauli

No, it doesn't. As long as the final result is less than 1e6145, it will return it - otherwise it overflows, as it should. There are hundreds of cases like this for which the current COMB overflows. It's nitpicking, perhaps - but to simulate an HP, you have to pay attention to detail :-)
Werner

[Werner]

(07-13-2020 12:10 PM)Albert Chan Wrote:  We can shift the argument: \( \binom{n}{r} = {n \over n-r}\; \binom{n-1}{r} \)

→ C(20408, 10204) = 2 × C(20407, 10204)

Doesn't always work: COMB(20421,10167)=COMB(20420,10167)*(20421/10167)
But Free42's COMB can't compute COMB(20420,10167) either.

Werner

[Thomas Okken]

That looks like a good fix. I'll apply it in the next release.

Unless any more serious bugs turn up, the next release won't be very soon. I think I'll hold off until Free42+ is ready, so I can release it and a regular Free42 update at the same time.

(It won't actually be called Free42+ but I'll just use that as the working title. Work begins this week...)

[Werner]

(07-13-2020 12:22 PM)Werner Wrote:  There are hundreds of cases like this for which the current COMB overflows. It's nitpicking, perhaps - but to simulate an HP, you have to pay attention to detail :-)
Werner

Not so nitpicking, take COMB(n,p) with p= 1229, just an example.
The built-in COMB can compute C(n,1229) up till n=45116057, returning Out of Range for larger n.
My RPN routine (in another thread) goes up to n=45377959. Quite a difference.

Thanks for the fix, Thomas, and of course it isn't urgent.

Cheers, Werner

[Thomas Okken]

I'm coding the fix so that it starts the product with 1/s instead of 1, and then multiplies by s at the end, where s is 10^(1+floor(log10(p))), or 2^(1+floor(log2(p))) in the binary version. That fixes all the cases mentioned in this thread.

https://github.com/thomasokken/free42/co...487ac559e6

[Werner]

(07-14-2020 03:38 AM)Thomas Okken Wrote:  I'm coding the fix so that it starts the product with 1/s instead of 1, and then multiplies by s at the end, where s is 10^(1+floor(log10(p))), or 2^(1+floor(log2(p))) in the binary version. That fixes all the cases mentioned in this thread.

https://github.com/thomasokken/free42/co...487ac559e6

For Free42 Decimal, the largest value of C(n,p)*p that overflows while C(n,p) does not is for n=20421,p=10167, and that is also the only case that requires the constant s to be 1E5 i.o. 1E4. 1E5 will thus work for all cases in Free42 Decimal.
Likewise, s=2^10 will work for all cases in Free42 Binary. Anything larger will also work, of course.

Cheers, Werner