fix(combinatorics): stirlingS2 returns plain number instead of BigNumber for S(n,n) with n>0

AUTHORS

@@ -287,6 +287,7 @@ Dheemanth D <165369664+Dheemanth07@users.noreply.github.com>
 Adriano Fantini <adrfantini@users.noreply.github.com>
 Chibuike Anselm <chibuikeanselm@Chibuikes-MacBook-Pro.local>
 Jos de Jong <email@josdejong.com>
+JSap0914 <JSap0914@users.noreply.github.com>
 rrbharath <ryanrbharath@gmail.com>
 Lê Vạn Bảo Trọng <32262183+chuanghiduoc@users.noreply.github.com>
 Hugohong258 <ugfifkgfufyu@gmail.com>

HISTORY.md

@@ -2,6 +2,8 @@
 
 # unpublished changes since 15.2.0
 
+- Feat: support bigint in stirlingS2(n,k) and ensure that result is
+  always of correct type (#3663). Thanks @JSap0914 and @gwhitney.
 - Docs: #3401 add pages for String and Boolean types, clarify that concat
   works on strings, and correct obsolete string manipulation examples (#3745)
   Thanks @Hugohong258.

src/function/combinatorics/stirlingS2.js

@@ -1,5 +1,4 @@
 import { factory } from '../../utils/factory.js'
-import { isNumber } from '../../utils/is.js'
 
 const name = 'stirlingS2'
 const dependencies = [
@@ -15,7 +14,9 @@ const dependencies = [
   'isInteger',
   'number',
   '?bignumber',
-  'larger'
+  '?bigint',
+  'larger',
+  'typeOf'
 ]
 
 export const createStirlingS2 = /* #__PURE__ */ factory(name, dependencies, (
@@ -32,73 +33,95 @@ export const createStirlingS2 = /* #__PURE__ */ factory(name, dependencies, (
     isInteger,
     number,
     bignumber,
-    larger
+    bigint,
+    larger,
+    typeOf
   }
 ) => {
-  const smallCache = []
-  const bigCache = []
+  const caches = { number: [], bignumber: [], bigint: [] }
   /**
-   * The Stirling numbers of the second kind, counts the number of ways to partition
-   * a set of n labelled objects into k nonempty unlabelled subsets.
+   * The Stirling numbers of the second kind, counts the number of ways to
+   * partition a set of n labelled objects into k nonempty unlabelled subsets.
    * stirlingS2 only takes integer arguments.
    * The following condition must be enforced: k <= n.
    *
    *  If n = k or k = 1 <= n, then s(n,k) = 1
    *  If k = 0 < n, then s(n,k) = 0
    *
-   * Note that if either n or k is supplied as a BigNumber, the result will be
-   * as well.
+   * Note that if either n or k is supplied as a bigint, the result will be as
+   * well. If not, and if either is a BigNumber, the result will be a BigNumber.
+   * Otherwise, the result will be a plain JavaScript number.
    *
    * Syntax:
    *
    *   math.stirlingS2(n, k)
    *
    * Examples:
    *
-   *    math.stirlingS2(5, 3) //returns 25
+   *    math.stirlingS2(5, 3) // returns 25
+   *    math.stirlingS2(math.bignumber(8), 3) // returns BigNumber 966
+   *    math.stirlingS2(18n, 3) // returns 64439010n
    *
    * See also:
    *
    *    bellNumbers
    *
-   * @param {Number | BigNumber} n    Total number of objects in the set
-   * @param {Number | BigNumber} k    Number of objects in the subset
-   * @return {Number | BigNumber}     S(n,k)
+   * History:
+   *
+   *    v15.3  added bigint support
+   *    v1.7   created
+   *
+   * @param {Number | BigNumber | bigint} n    Total number of objects in the set
+   * @param {Number | BigNumber | bigint} k    Number of subsets
+   * @return {Number | BigNumber}     S2(n,k)
    */
   return typed(name, {
-    'number | BigNumber, number | BigNumber': function (n, k) {
+    'number | BigNumber | bigint, number | BigNumber | bigint': function (n, k) {
       if (!isInteger(n) || isNegative(n) || !isInteger(k) || isNegative(k)) {
-        throw new TypeError('Non-negative integer value expected in function stirlingS2')
+        throw new TypeError(
+          'Non-negative integer value expected in function stirlingS2')
       } else if (larger(k, n)) {
-        throw new TypeError('k must be less than or equal to n in function stirlingS2')
+        throw new TypeError(
+          'k must be less than or equal to n in function stirlingS2')
+      }
+
+      let cache = caches.number
+      let make = number
+      if (typeOf(n) === 'bigint' || typeOf(k) === 'bigint') {
+        cache = caches.bigint
+        make = bigint
+      } else if (typeOf(n) === 'BigNumber' || typeOf(k) === 'BigNumber') {
+        cache = caches.bignumber
+        make = bignumber
       }
+      // Since it is impractical to compute stirling numbers for indices
+      // outside the range of plain number, safe now to convert:
+      n = number(n)
+      k = number(k)
 
-      const big = !(isNumber(n) && isNumber(k))
-      const cache = big ? bigCache : smallCache
-      const make = big ? bignumber : number
-      const nn = number(n)
-      const nk = number(k)
       /* See if we already have the value: */
-      if (cache[nn] && cache[nn].length > nk) {
-        return cache[nn][nk]
+      if (cache[n] && cache[n].length > k) {
+        return cache[n][k]
       }
       /* Fill the cache */
-      for (let m = 0; m <= nn; ++m) {
-        if (!cache[m]) {
-          cache[m] = [m === 0 ? make(1) : make(0)]
-        }
-        if (m === 0) continue
-        const row = cache[m]
-        const prev = cache[m - 1]
-        for (let i = row.length; i <= m && i <= nk; ++i) {
-          if (i === m) {
-            row[i] = 1
-          } else {
-            row[i] = addScalar(multiplyScalar(make(i), prev[i]), prev[i - 1])
+      if (n >= cache.length || k >= cache[n].length) {
+        for (let m = 0; m <= n; ++m) {
+          if (!cache[m]) {
+            cache[m] = [m === 0 ? make(1) : make(0)]
+          }
+          if (m === 0) continue
+          const row = cache[m]
+          const prev = cache[m - 1]
+          for (let i = row.length; i <= m && i <= k; ++i) {
+            if (i === m) {
+              row[i] = make(1)
+            } else {
+              row[i] = addScalar(multiplyScalar(make(i), prev[i]), prev[i - 1])
+            }
           }
         }
       }
-      return cache[nn][nk]
+      return cache[n][k]
     }
   })
 })

test/node-tests/doc.test.js

@@ -377,16 +377,14 @@ const knownUndocumented = new Set([
 const knownNoHistory = new Set([
   'config', 'typed', 'derivative', 'leafCount', 'lsolve', 'lsolveAll', 'lup',
   'lusolve', 'lyap', 'qr', 'rationalize', 'resolve', 'schur', 'simplify',
-  'simplifyConstant', 'simplifyCore', 'slu', 'sylvester', 'symbolicEqual',
-  'usolve', 'usolveAll', 'abs', 'add', 'cbrt', 'ceil', 'cube', 'divide',
-  'dotDivide', 'dotMultiply', 'dotPow', 'exp', 'expm', 'expm1', 'fix', 'floor',
-  'gcd', 'hypot', 'lcm', 'log10', 'log1p', 'log2', 'mod', 'multiply',
-  'norm', 'nthRoot', 'nthRoots', 'pow', 'round', 'sign', 'sqrt', 'sqrtm',
-  'square', 'subtract', 'unaryMinus', 'unaryPlus', 'xgcd', 'bitAnd', 'bitNot',
-  'bitOr', 'bitXor', 'leftShift', 'rightArithShift', 'rightLogShift',
-  'bellNumbers', 'catalan', 'composition', 'stirlingS2', 'arg', 'conj', 'im',
+  'simplifyConstant', 'slu', 'sylvester', 'symbolicEqual',
+  'usolve', 'usolveAll', 'abs', 'cube', 'divide', 'expm',
+  'hypot', 'log10', 'log1p', 'log2', 'multiply',
+  'norm', 'pow', 'round', 'sign', 'sqrt', 'sqrtm',
+  'square', 'unaryMinus', 'unaryPlus', 'bitNot',
+  'bellNumbers', 'catalan', 'composition', 'arg', 'conj', 'im',
   're', 'compile', 'evaluate', 'help', 'parser', 'distance', 'intersect',
-  'and', 'not', 'nullish', 'or', 'xor', 'column', 'concat', 'count', 'cross',
+  'not', 'nullish', 'column', 'count', 'cross',
   'ctranspose', 'det', 'diag', 'diff', 'dot', 'eigs', 'fft', 'filter',
   'flatten', 'forEach', 'getMatrixDataType', 'identity', 'ifft', 'inv', 'kron',
   'mapSlices', 'matrixFromColumns', 'matrixFromFunction', 'matrixFromRows',
@@ -395,16 +393,15 @@ const knownNoHistory = new Set([
   'transpose', 'zeros', 'solveODE', 'bernoulli', 'combinations',
   'combinationsWithRep', 'factorial', 'gamma', 'kldivergence', 'lgamma',
   'multinomial', 'permutations', 'pickRandom', 'random', 'randomInt',
-  'compareNatural', 'compareText', 'deepEqual', 'equal', 'equalText', 'larger',
-  'largerEq', 'smaller', 'smallerEq', 'unequal', 'setCartesian',
-  'setDifference', 'setDistinct', 'setIntersect', 'setIsSubset',
+  'compareText', 'equalText',
+  'setCartesian', 'setDifference', 'setDistinct', 'setIntersect', 'setIsSubset',
   'setMultiplicity', 'setPowerset', 'setSize', 'setSymDifference', 'setUnion',
   'freqz', 'zpk2tf', 'erf', 'zeta', 'corr', 'cumsum', 'mad', 'max', 'mean',
   'median', 'min', 'mode', 'prod', 'quantileSeq', 'std', 'sum', 'variance',
   'acos', 'acosh', 'acot', 'acoth', 'acsc', 'acsch', 'asec', 'asech', 'asin',
-  'asinh', 'atan', 'atan2', 'atanh', 'cos', 'cosh', 'cot', 'coth', 'csc',
-  'csch', 'sec', 'sech', 'sin', 'sinh', 'tan', 'tanh', 'to', 'toBest', 'bin',
-  'clone', 'hasNumericValue', 'hex', 'isBounded', 'isFinite', 'isInteger',
+  'asinh', 'atan', 'atanh', 'cos', 'cosh', 'cot', 'coth', 'csc',
+  'csch', 'sec', 'sech', 'sin', 'sinh', 'tan', 'tanh', 'toBest', 'bin',
+  'clone', 'hasNumericValue', 'hex', 'isBounded', 'isInteger',
   'isNaN', 'isNegative', 'isNumeric', 'isPositive', 'isPrime', 'isZero',
   'oct', 'print', 'typeOf', 'bignumber', 'chain', 'number', 'splitUnit'
 ])

test/unit-tests/function/combinatorics/stirlingS2.test.js

@@ -3,6 +3,8 @@ import math from '../../../../src/defaultInstance.js'
 const stirlingS2 = math.stirlingS2
 
 describe('stirlingS2', function () {
+  const bn = math.bignumber
+
   it('should calculate the number of ways to partition a set of n objects into k non-empty subsets', function () {
     assert.strictEqual(stirlingS2(5, 3), 25)
     assert.strictEqual(stirlingS2(0, 0), 1)
@@ -15,7 +17,6 @@ describe('stirlingS2', function () {
   })
 
   it('should calculate the stirlingS2 of n items taken k at a time with BigNumbers', function () {
-    const bn = math.bignumber
     assert.deepStrictEqual(stirlingS2(bn(7), bn(5)), bn(140))
     assert.deepStrictEqual(stirlingS2(bn(8), bn(6)), bn(266))
     assert.deepStrictEqual(stirlingS2(bn(14), 13), bn(91))
@@ -26,6 +27,21 @@ describe('stirlingS2', function () {
     assert.deepStrictEqual(stirlingS2(bn(23), 9), bn('12320068811796900'))
     assert.deepStrictEqual(stirlingS2(bn(50), 14),
       bn('16132809270066494376125322988035691981158490930'))
+    // S(n,n) = 1 must return BigNumber when both inputs are BigNumber (n > 0)
+    // S(n,n) = 1 must return BigNumber when both inputs are BigNumber (n > 0)
+    assert.deepStrictEqual(stirlingS2(bn(1), bn(1)), bn(1))
+    assert.deepStrictEqual(stirlingS2(bn(3), bn(3)), bn(1))
+    assert.deepStrictEqual(stirlingS2(bn(7), bn(7)), bn(1))
+    // mixed: one BigNumber input still produces BigNumber when n===k
+    assert.deepStrictEqual(stirlingS2(bn(4), 4), bn(1))
+  })
+
+  it('should calculate stirlingS2 with bigints', function () {
+    assert.strictEqual(stirlingS2(5, 0n), 0n)
+    assert.strictEqual(stirlingS2(7n, 4), 350n)
+    assert.strictEqual(stirlingS2(30n, 15n), 12879868072770626040000n)
+    assert.strictEqual(stirlingS2(bn(8), 5n), 1050n)
+    assert.strictEqual(stirlingS2(19n, bn(19)), 1n)
   })
 
   it('should not work with non-integer and negative input', function () {