+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: Not HP Calculators (/forum-7.html)
+--- Forum: Not quite HP Calculators - but related (/forum-8.html)
+--- Thread: (Free42) not a bug - COMB (/thread-15338.html)
(Free42) not a bug - COMB - Werner - 07-13-2020 11:36 AM
.. 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
RE: (Free42) not a bug - COMB - Paul Dale - 07-13-2020 12:03 PM
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
RE: (Free42) not a bug - COMB - Albert Chan - 07-13-2020 12:10 PM
We can shift the argument: \( \binom{n}{r} = {n \over n-r}\; \binom{n-1}{r} \)
→ C(20408, 10204) = 2 × C(20407, 10204)
RE: (Free42) not a bug - COMB - Werner - 07-13-2020 12:22 PM
(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
RE: (Free42) not a bug - COMB - Werner - 07-13-2020 12:28 PM
(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} \)Doesn't always work: COMB(20421,10167)=COMB(20420,10167)*(20421/10167)
→ C(20408, 10204) = 2 × C(20407, 10204)
But Free42's COMB can't compute COMB(20420,10167) either.
Werner
RE: (Free42) not a bug - COMB - Thomas Okken - 07-13-2020 12:58 PM
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...)
RE: (Free42) not a bug - COMB - Werner - 07-13-2020 01:16 PM
(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 :-)Not so nitpicking, take COMB(n,p) with p= 1229, just an example.
Werner
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
RE: (Free42) not a bug - COMB - Thomas Okken - 07-14-2020 03:38 AM
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/commit/50675effbaed4ce885dce5a91aed98487ac559e6
RE: (Free42) not a bug - COMB - Werner - 07-14-2020 06:57 AM
(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/commit/50675effbaed4ce885dce5a91aed98487ac559e6
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