math.stirlingS2(n, i) and hence math.bellNumbers(n) rounding errors and incorrect values

[brubsby]

I wrote the following line of js to reproduce this page from oeis:

console.log([...Array(27).keys()].reduce((prev, n)=>prev += "\n" + n.toString().padStart(2)+math.bellNumbers(n).toString().padStart(22), ""))

but the output has floating point errors and incorrect values

 0                     1
 1                     1
 2                     2
 3                     5
 4                    15
 5                    52
 6                   203
 7                   877
 8                  4140
 9                 21147
10                115975
11                678570
12               4213597
13              27644437
14             190899322
15            1382958545
16           10480142147
17     82864869804.00002
18     682076806158.9999
19         5832742205057
20        51724158235372
21    474869816156751.06
22      4506715738447324
23     44152005855084340
24    445958869294805250
25   4638590332230000000
26  49631246523618755000

attempting to work around this default behavior by passing in n as a bignumber also doesn't work:

console.log([...Array(27).keys()].reduce((prev, n)=>prev += "\n" + n.toString().padStart(2)+math.bellNumbers(math.bignumber(n)).toString().padStart(22), ""))
Uncaught TypeError: Cannot implicitly convert a number with >15 significant digits to BigNumber (value: 51090942171709440000). Use function bignumber(x) to convert to BigNumber.
    at convert (math.js:40:116329)
    at Array.C.t (math.js:40:28692)
    at Function.r (math.js:40:30983)
    at Function.ae (math.js:40:31507)
    at divideScalar (math.js:40:31921)
    at Function.number | BigNumber, number | BigNumber (math.js:40:555793)
    at stirlingS2 (math.js:40:31734)
    at Function.number | BigNumber (math.js:40:556157)
    at Object.bellNumbers (math.js:40:31680)
    at <anonymous>:1:98

both of which seem like incorrect behaviors

[gwhitney]

@brubsby thank you so much for the report. Note that JavaScript's "max safe integer" is 2^53-1, so we would only expect correct values in the number version up to B_22, as B_23 is larger than this. Yet as you report there are floating-point errors well before this; in particular, isInteger(bellNumbers(17)) returns false, which is nonsense. So that is definitely a bug that needs tracking down.

And indeed, you should be able to pass a bignumber to bellNumbers() to get greater accuracy, yet that is not working either, so that is also an additional bug.

I think the first place to look is whether the stirling numbers of the second kind are behaving properly or not.

[gwhitney]

Ah, indeed: stirlingS2(14, 13) returns 91.00000000048178, and isInteger returns false on that. So that should be tracked down first.

[gwhitney]

Ah, yes, the current implementation for Stirling numbers uses the sum formula, which suffers from producing very large intermediate values. The recurrence relation is much more practical for accurate computation. I will file a PR switching the implementation as soon as possible. That will not necessarily cure all evils here, but it will be a good start.

[brubsby]

Thank you! I will admit this is my first time using the library, but it seems like the function should return at least B_23 and up as a bignumber by default (after the bugs are fixed), regardless of the input parameter's type. I know this is kind of against the general rule for the library, but asking for Bell number 23 seems like explicit enough instruction to use a higher precision bignumber.

Something about asking a math library for an exact integer from a sequence and getting a different number, when it could've given you the correct one, seems wrong.

Although there would eventually come a point where the sequence grows too large for bignumbers as well, so I guess it could go either way. But I'm not sure if there would be any negative consequences to returning bignumbers when the user wasn't expecting it.

[gwhitney]

OK, so I've just submitted the PR which should allow you to create the table you want, up to B_22 for number input, and I actually don't know how large it will manage with BigNumber input, I only tried up to B_50 and it didn't break a sweat.

should return at least B_23 and up as a bignumber by default (after the bugs are fixed), regardless of the input parameter's type

As you point out, that is very much against the design of mathjs, so I don't see that as a change that will happen. Remember, one of the main goals of mathjs is full modularity as to the selection of types that are available, and indeed it by default provides not just the main module but also one that has support only for the built-in JavaScript number type. So pretty much bellNumber(n) when n is an ordinary JavaScript number is obliged to return a number. That's why mathjs also attempts to educate its clients about the challenges the number type faces.

A personal project of mine will be to extend mathjs with (possibly optionally loaded, through an extension module) full native bigint support throughout all relevant functions, at which point you will, in principle, be able to compute arbitrarily large Bell numbers.

[josdejong]

yes indeed, the general idea is number in ->number out, BigNumber in -> BigNumber out.

BigInt support will be very welcome in mathjs itself too 😄 . We have to be careful with breaking IE11 support right now though.