[WP34S] A suggestion about C(x,y)



#3

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


#4

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:-)


#5

Negative integral arguments fail due to the poles in the gamma function. Negative non-integral arguments should work.

- Pauli


#6

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 ...

#7

I was mistaken. For the real case, the code checks for either argument being negative and fails early.


- Pauli

#8

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


#9

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


#10

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


#11

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:-/


#12

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


#13

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


#14

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


#15

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


#16

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.

#17

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.


#18

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:-/


Possibly Related Threads…
Thread Author Replies Views Last Post
  Prime Connectivity Kit Suggestion toml_12953 1 1,679 12-06-2013, 10:41 PM
Last Post: Michael de Estrada
  HP Prime suggestion to avoid Numeric/Symbolic confusion Chris Pem10 4 1,918 11-19-2013, 05:49 AM
Last Post: bluesun08
  HP Prime - Revision Suggestion - Setting the Clock Bill Triplett 5 2,249 11-15-2013, 12:36 AM
Last Post: Joe Horn
  HP Prime - Cross product suggestion bluesun08 13 4,759 11-08-2013, 01:49 AM
Last Post: Patrice
  HP Prime Emulator/Conn. kit suggestion Han 2 1,538 09-27-2013, 11:23 AM
Last Post: Han
  WP-34s feature suggestion: "Follow jump" shortcut Andrew Nikitin 3 1,555 06-12-2013, 01:42 AM
Last Post: Walter B
  [WP34S] WP34S firmware on the AT91SAM7L-STK dev kit? jerome ibanes 1 1,248 10-04-2012, 04:59 PM
Last Post: Paul Dale
  wp34s Emulator suggestion Damir 25 6,670 06-07-2012, 05:02 PM
Last Post: Chris Tvergard
  [wp34s] Incomplete Gamma on the wp34s Les Wright 18 5,348 12-06-2011, 11:07 AM
Last Post: Namir
  [wp34s] Romberg Integration on the wp34s Les Wright 2 1,519 12-04-2011, 10:49 PM
Last Post: Les Wright

Forum Jump: