▼
Posts: 56
Threads: 10
Joined: May 2006
Hi,
I do not know any calculator that calculates the binomial coefficients C(n,k) when n is not an integer. Of course k must be an integer >=0, but n can be any real or complex number. To be precise
C(z,k) = z(z-1)(z-2)...(z-k+1)/k!
is the coefficient of x^k in the binomial expansion of (1+x)^z (for |x|<1, of course).
For the sake of motivation, one may ask what what are the second, fourth, sixth, eight-order, etc. asymptotic expansions of the Lorentz factor
1/sqrt(1-v^2/c^2) = 1 + ...
from relativity. The coefficients of (v/c)^2, (v/c)^4, (v/c)^6, (v/c)^8, etc., are simply
-C(-1/2,1) = 0.5
C(-1/2,2) = 0.375
-C(-1/2,3) = 0.3125
C(-1/2,4) = 0.2734375
As far as I could check neither the 34S nor even the HP50G implement C(n,k) other than for integers with n>=k>=0 (returning Domain Errors otherwise), and yet no essential change in the algorithm should be needed to extend the domain (if anything, less error-checking would be necessary).
Besides the obvious application to finding Taylor coefficients, there are combinatorial ones. For instance, whereas C(n,k) is the number of ways of choosing a k-element subset of an n-element set; we have that abs(C(-n,k)) is the number of ways of choosing a k-element multiset from an n-set; in other words, a k-element "subset" chosen from an n-element set, but allowing repeated choices.
[In reality the above can be reduced to the usual binomial coefficients via C(-n,k) = (-1)^k C(n+k-1,k), but it's hard to beat the beauty and simplicity of the formula abs(C(-n,k)). However, note that this only applies to negative integer values of n and not in general to other real values, much less complex values of n.]
Any hope that this modification, both useful and likely trivial to implement, could make its way to a future ROM? Even the extension of C(x,y) to all real x (even if not complex) would be useful, and both welcome and unique among all calculators.
Cheers,
Eduardo
▼
Posts: 4,587
Threads: 105
Joined: Jul 2005
Buenas dias Eduardo,
Good news first: C y,x works for complex numbers as well, and it also computes C 0.5,1 .
The further extension to negative arguments is not implemented, however. We'll discuss that.
d:-)
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
Negative integral arguments fail due to the poles in the gamma function. Negative non-integral arguments should work.
- Pauli
▼
Posts: 4,587
Threads: 105
Joined: Jul 2005
Quote:
Negative non-integral arguments should work.
But the WP 34S doesn't know this and fails ;-) Try, for example, -0.5 ENTER 1 Cy,x ...
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
I was mistaken. For the real case, the code checks for either argument being negative and fails early.
- Pauli
Posts: 56
Threads: 10
Joined: May 2006
Quote:
Negative integral arguments fail due to the poles in the gamma function. Negative non-integral arguments should work.
- Pauli
Negative integers can be reduced to positive ones, per
C(-n,k) = (-1)^k * C(n+k-1,k)
as I had preemptively remarked.
Is C(x,y) implemented using the gamma function in the 34s? Of course we have the formula
C(x,y) = Gamma(x+1)/(Gamma(y+1)*Gamma(x-y+1))
but my instincts say direct application of this formula would lead to large rounding errors in some cases and overflow in others.
Eduardo
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
These functions are implemented using LnGamma. This avoids the overflow issues, although there is still some risk of loss of accuracy I think.
The reflection formula you mention could be used. Is there similar for permutations?
- Pauli
▼
Posts: 56
Threads: 10
Joined: May 2006
Quote:
These functions are implemented using LnGamma. This avoids the overflow issues, although there is still some risk of loss of accuracy I think.
The reflection formula you mention could be used. Is there similar for permutations?
- Pauli
Pauli,
I'm certain that some accuracy is lost using LnGamma to compute C(x,y). This must definitely be the case when x is close to a negative integer. I can swear there must be algorithms to give a better approximation in such cases, but I'm not an expert. If there is a chance to possibly re-implementing C(x,y), I can try looking in the specialized literature.
Your question about the permutations is absolutely pertinent. For any real or complex z and k an integer >=0, define
P(z,k) = z(z-1)(z-2)...(z-k+1)
Then what I asked at the top of the thread was to extend the definition of C(z,k) as P(z,k)/k! .
The above definition of P(z,k) is known in discrete math/computer science under the quirky name of "z to the k-th falling power". Using a notation introduced, I believe, by Donald Knuth, it is denoted by z^k with the exponent k *underlined*. There is a similar definition of z to the k-th rising power, z^k with the exponent k *overlined* (looks like k is complex-conjugated) given by
z(z+1)(z+2)...(z+k-1)
but I will presently focus only on the falling powers, that is P(z,k) as defined above.
In discrete mathematics, falling powers are very useful. P(z,k) are polynomials in z of degree k, and they are the natural polynomials to use in discrete interpolation. The "finite-differences formula", the discrete analog of Tayor's formula, reads
f(z) = c0 + c1*P(z,1) + c2*P(z,2)/2! + c3*P(z,3)/3! + ... +cn*P(z,n)/n!
= c0 + c1*C(z,1) + c2*C(z,2) + c3*C(z,3) + ... +cn*C(z,n)
The above formula gives the only polynomial of degree n with given differences
c0 = f(0)
c1 = Delta_1(f)(0) = f(1)-f(0)
c2 = Delta_2(f)(0)/2! = f(2)-2f(1)+f(0)
...
(I won't write the definition of the higher difference quotients to avoid too much clutter.)
The difference quotients depend only on the values f(0), f(1), ..., f(n), so the finite-differences formula is basically a more explicit and cleaner form of the Lagrange polynomial interpolation formula in the special (but very common and useful case) when the values of the interpolated function are known at n+1 *equally spaced* points. It's easy to write a user program to compute the finite differences of a given function. In conjunction with P(z,x), this would solve many polynomial interpolation problems easily. (Heck, I do volunteer to write a quick-n-dirty library to compute the finite differences if there's interest.)
Summing up the above, it makes perfect sense to extend both P(z,k) and C(z,k) to all complex z and integer k>=0. Of course using the Gamma function one may extend the above even to arbitrary real or complex k; this is almost the same as asking for the calculation of Euler's beta (which the 34S computes already).
Finally, we have the following general reflection formula (directly from the definition):
P(-z,k) = (-1)^k P(z+k-1,k)
and the same identity holds for C(x,y) instead of P(x,y).
Eduardo
▼
Posts: 4,587
Threads: 105
Joined: Jul 2005
Buenas tardes Eduardo,
Quote:
Summing up the above, it makes perfect sense to extend both P(z,k) and C(z,k) to all complex z and integer k>=0.
I doubt it. You wrote in the OP:
Quote:
Any hope that this modification, both useful and likely trivial to implement, could make its way to a future ROM?
The more we go into details the less trivial it becomes obviously. And the use isn't overwhelming either.
Quote:
Of course using the Gamma function one may extend the above even to arbitrary real or complex k; this is almost the same as asking for the calculation of Euler's beta (which the 34S computes already).
This argument sounds pretty "financial" to me. The WP 34S already contains a lot of ... ummh ... off-standard functions - we have to draw a line once. Feel free to program whatever you like and as long as you like - memory is large enough. But as soon as individual requests generate extra overhead - I simply doubt the benefits balancing the costs. Sorry.
d:-/
▼
Posts: 56
Threads: 10
Joined: May 2006
Walter,
Quote:
This argument sounds pretty "financial" to me. The WP 34S already contains a lot of ... ummh ... off-standard functions - we have to draw a line once. Feel free to program whatever you like and as long as you like - memory is large enough. But as soon as individual requests generate extra overhead - I simply doubt the benefits balancing the costs. Sorry.
Did you read "finite math" as "finance math" in my earlier post? I don't understand where "financial" came from above.
I don't know whether the 34s implements C(n,k) using the Gamma function or else computing a k-fold product, but I definitely had the latter in mind writing my original post. In such case, the exact same algorithm should work, provided only that the routine is not exited (it is not a genuine Domain Error to provide non-integer n).
Even if the 34s uses Gamma to compute C(n,k), the exact same algorithm would work *except*, in this case, for negative integers n (but otherwise for any real or complex n), although precision would be lost as soon as n-k+1 is close to a negative integer. Here, of course, I'm assuming Gamma defined for all complex numbers except its poles.
What I really should do is look at the code for C(n,k) and P(n,k) and see just what change(s) would be needed. I really can see that the tiniest of changes would extend the domain. I only blame myself for not making the usefulness of these enhancements appealing enough with the explicit examples of applicability I provided (binomial expansion of fractional and complex powers, finite-differences method of polynomial interpolation). As a mathematician, I always end up writing too many details that tend to lengthen posts...
Eduardo
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
The code for C(n,k) and P(n,k) lives in two places -- one for the real implementation and another for the complex.
The real version start from that link -- combinations first and permutations below. The complex flavours live elsewhere. Changing the latter will be much easier since they are keystroke programs. The C code for the real versions would be more difficult to change & space is much more of a concern here but it certainly would be possible.
As I've mentioned before, the 34S uses log gamma to compute both. Adding the reflection formula would avoid the precision issues for values approaching negative integers since they wouldn't occur.
- Pauli
▼
Posts: 56
Threads: 10
Joined: May 2006
Thanks for the links. I will take a look.
Eduardo
Quote:
The code for C(n,k) and P(n,k) lives in two places -- one for the real implementation and another for the complex.
The real version start from that link -- combinations first and permutations below. The complex flavours live elsewhere. Changing the latter will be much easier since they are keystroke programs. The C code for the real versions would be more difficult to change & space is much more of a concern here but it certainly would be possible.
- Pauli
Posts: 3,229
Threads: 42
Joined: Jul 2006
Quote:
There is a similar definition of z to the k-th rising power, z^k with the exponent k *overlined* (looks like k is complex-conjugated) given by
z(z+1)(z+2)...(z+k-1)
but I will presently focus only on the falling powers, that is P(z,k) as defined above.
My plans for extra factorial functions got vetoed :-( I wrote code for n!! and !n (double factorial and sub-factorial) functions. There are half a dozen or more extra factorial like functions that could have been included if we had space. Rising and falling factorials would have been part of the suite.
- Pauli
▼
Posts: 3,283
Threads: 104
Joined: Jul 2005
It looks like the implementation can easily done in keystroke programs and might become part of the library. No need to load the suite of function for those who don't need it. The only drawback: The self made functions will have trouble to preserve the stack the same way internal code does.
Posts: 735
Threads: 34
Joined: May 2007
Quote:
As far as I could check neither the 34S nor even the HP50G implement C(n,k) other than for integers with n>=k>=0 (returning Domain Errors otherwise)
For both HP-42S and HP-15C we can use simple programs using the formula:
HP-42S
00 { 24-Byte Prgm }
01 LBL "COMB"
02 X<>Y
03 1
04 +
05 GAMMA
06 LASTX
07 RCL- ST Z
08 GAMMA
09 /
10 X<>Y
11 N!
12 /
13 END
HP-15C
001 - 42,21,13 LBL C
002 - 36 ENTER
003 - 36 ENTER
004 - 42 0 x!
005 - 43 33 R^
006 - 42 0 x!
007 - 43 36 LSTx
008 - 43 33 R^
009 - 30 -
010 - 42 0 x!
011 - 10 /
012 - 34 x<>y
013 - 10 /
014 - 43 32 RTN
Kind regards
Thomas
PS: This "trick" works as well for the
HP 20b.
Edited: 28 Jan 2013, 3:59 a.m.
▼
Posts: 4,587
Threads: 105
Joined: Jul 2005
Correct - you'll get an early overflow, however, when Gamma(x+1) overflows. Not really relevant in Eduardo's problem, but for COMBs with large arguments. And we've got just one code.
d:-/
|
