Should we still keep input validation to check that the input dimensions allow broadcasting? Or are the generic error messages good enough in this case?

OP here again.  Thank you for the suggestions involving the symbolic package.

> The accuracy is limited not by the gammaln() function, but by the addition operation using IEEE 8-byte doubles.  eps() is 2e-16 and if the difference in scale between two arguments is greater than 1e16 then the smaller argument is lost.

The accuracy of betaln(a, b) can be poor even if the addition operation a + b is exact.  For example betaln(2^62, 2^10) returns -32768 for me, but the actual value is about -37935.2477864453572569.



I am interested in the possibility of contributing a version of betaln.m that improves accuracy at some cost to speed and simplicity.  I have a few questions:

1. How likely is this to be accepted?  Assume for the sake of argument that accuracy is at least 40 bits for the vast majority of the input space, and that speed and simplicity are about 10 times worse than the existing implementation.
2. The existing implementation has input validation that requires the arguments to have equal size or at least one of them to be scalar.  Is there any reason to retain this and not fully support broadcasting?
3. Some definitions of the beta function seem to require (the real parts of) both arguments to be positive, but others relax this requirement.  What should betaln return with one or more non-positive arguments?

Any input would be appreciated.  Thanks.

I've never seen input values to nchoosek like 2^128.  I think that is well outside the norm which is probably why this hasn't been noticed before.  I'm not saying your application doesn't call for these numbers, just that overwhelming number of use cases don't call for these values.

Generally nchoosek blows up very quickly such that it isn't conceivable to try the direct approach of

log (nchoosek (n, k))

Something simple like

nchoosek (1e6,4)
warning: nchoosek: possible loss of precision
warning: called from
    nchoosek at line 151 column 7

ans = 4.1666e+22

is already complaining about lost precision.

Back to betaln:

lnb = gammaln (a) + gammaln (b) - gammaln (a + b)

The accuracy is limited not by the gammaln() function, but by the addition operation using IEEE 8-byte doubles.  eps() is 2e-16 and if the difference in scale between two arguments is greater than 1e16 then the smaller argument is lost.

octave:19> format long
octave:20> 1 + 1e-16
ans = 1
octave:21> sprintf ('%.17d\n', 1 + 2e-16)
ans = 1.0000000000000002

I would go back to the base definition of nchoosek as

c = n! / k!(n-k)!

Take the logarithm of both sides and use the properties of logarithms to write

log (c) = log (n!) - log (k!) - log ((n-k)!)

Then using the relationship

gamma (n + 1) = n!

I can write

log (c) = gammaln (n+1) - gammaln (k+1) - gammaln (n-k+1)

A quick anonymous function for this is

bignk = @(n,k) gammaln (n+1) - gammaln (k+1) - gammaln (n-k+1)

This still has a problem in that the expression "N-K+1" can't be calculated without loss of precision with IEEE doubles.  So, I tried switching to variable precision numbers which are part of the symbolic package.  You will need to have python3-sympy installed on your computer.  After that, install and load the symbolic package with

pkg install -forge symbolic
pkg load symbolic

Finally, create some numbers and try it out.

octave:47> n = vpa (2, 50)
n = (sym) 2.0000000000000000000000000000000000000000000000000
octave:48> n = n^128
n = (sym) 340282366920938463463374607431768211456.00000000000
octave:49> k = vpa (2, 50)
k = (sym) 2.0000000000000000000000000000000000000000000000000
octave:50> k = k^32
k = (sym) 4294967296.0000000000000000000000000000000000000000
octave:51> logc = bignk (n,k)
logc = (sym) 290091236578.66962544142734259366989135742187500000

Maybe this is what you are looking for?

I can confirm that at high ratios of a/b the betaln() function always returns zero.  The cause is straightforward.  In betaln.m the output is calculated as

  lnb = gammaln (a) + gammaln (b) - gammaln (a + b);

When a and b are different enough in size the finite precision of IEEE floating point causes a problem and the sum "a + b" is just equal to "a".

However, it's not clear to me that this is a frequent occurrence.  For example, do constraints from real-world problems limit the sizes of a, b, and their ratio?  I would also observe that Matlab exhibits exactly the same behavior as Octave.  Maybe they are unaware of this corner case, or maybe they consider it unlikely.  I also note that Mathematica doesn't even offer a log beta function.  The help discussions simply point to writing such a function in terms of gammaln in exactly the way Octave has done.