Potential Bug in mod()

AUTHORS

@@ -229,5 +229,6 @@ Michael Greminger <[email protected]>
 Kiku <[email protected]>
 MaybePixem <[email protected]>
 Aly Khaled <[email protected]>
+Praise Nnamonu <[email protected]>
 
 # Generated by tools/update-authors.js

src/function/arithmetic/gcd.js

@@ -1,14 +1,18 @@
+import { isInteger, nearlyEqual } from '../../utils/number.js'
+import { nearlyEqual as bigNearlyEqual } from '../../utils/bignumber/nearlyEqual.js'
 import { factory } from '../../utils/factory.js'
 import { createMatAlgo01xDSid } from '../../type/matrix/utils/matAlgo01xDSid.js'
 import { createMatAlgo04xSidSid } from '../../type/matrix/utils/matAlgo04xSidSid.js'
 import { createMatAlgo10xSids } from '../../type/matrix/utils/matAlgo10xSids.js'
 import { createMatrixAlgorithmSuite } from '../../type/matrix/utils/matrixAlgorithmSuite.js'
-import { gcdNumber } from '../../plain/number/index.js'
+// import { gcdNumber } from '../../plain/number/index.js'
 import { ArgumentsError } from '../../error/ArgumentsError.js'
 
 const name = 'gcd'
 const dependencies = [
   'typed',
+  'config',
+  'round',
   'matrix',
   'equalScalar',
   'BigNumber',
@@ -23,7 +27,7 @@ function is1d (array) {
   return !array.some(element => Array.isArray(element))
 }
 
-export const createGcd = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, BigNumber, DenseMatrix, concat }) => {
+export const createGcd = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, config, round, equalScalar, BigNumber, DenseMatrix, concat }) => {
   const matAlgo01xDSid = createMatAlgo01xDSid({ typed })
   const matAlgo04xSidSid = createMatAlgo04xSidSid({ typed, equalScalar })
   const matAlgo10xSids = createMatAlgo10xSids({ typed, DenseMatrix })
@@ -57,7 +61,7 @@ export const createGcd = /* #__PURE__ */ factory(name, dependencies, ({ typed, m
   return typed(
     name,
     {
-      'number, number': gcdNumber,
+      'number, number': _gcdNumber,
       'BigNumber, BigNumber': _gcdBigNumber,
       'Fraction, Fraction': (x, y) => x.gcd(y)
     },
@@ -89,6 +93,28 @@ export const createGcd = /* #__PURE__ */ factory(name, dependencies, ({ typed, m
     }
   )
 
+  /**
+ * Calculate gcd for numbers
+ * @param {number} a
+ * @param {number} b
+ * @returns {number} Returns the greatest common denominator of a and b
+ * @private
+ */
+  function _gcdNumber (a, b) {
+    if (!isInteger(a) || !isInteger(b)) {
+      throw new Error('Parameters in function gcd must be integer numbers')
+    }
+
+    // https://en.wikipedia.org/wiki/Euclidean_algorithm
+    let r
+    while (b !== 0) {
+      r = _modNumber(a, b)
+      a = b
+      b = r
+    }
+    return (a < 0) ? -a : a
+  }
+
   /**
    * Calculate gcd for BigNumbers
    * @param {BigNumber} a
@@ -104,10 +130,56 @@ export const createGcd = /* #__PURE__ */ factory(name, dependencies, ({ typed, m
     // https://en.wikipedia.org/wiki/Euclidean_algorithm
     const zero = new BigNumber(0)
     while (!b.isZero()) {
-      const r = a.mod(b)
+      const r = _modBigNumber(a, b)
       a = b
       b = r
     }
     return a.lt(zero) ? a.neg() : a
   }
+
+  /**
+   * Calculate the mod of two numbers
+   * @param {number} x
+   * @param {number} y
+   * @returns {number} res
+   * @private
+   */
+  function _modNumber (x, y) {
+    if (y === 0) {
+      return x
+    }
+    // then y < 0 or y > 0
+    const div = x / Math.abs(y)
+    if (nearlyEqual(div, round(div), config.epsilon)) {
+      const result = x - y * round(div)
+      return nearlyEqual(result, round(result), config.epsilon)
+        ? round(result)
+        : result
+    } else {
+      return x - Math.abs(y) * Math.floor(x / Math.abs(y))
+    }
+  }
+
+  /**
+   * Calculate the mod of two BigNumbers
+   * @param {number} x
+   * @param {number} y
+   * @returns {number} res
+   * @private
+   */
+  function _modBigNumber (x, y) {
+    if (y === 0) {
+      return x
+    }
+    // then y < 0 or y > 0
+    const div = x.div(y.abs())
+    if (bigNearlyEqual(div, round(div), config.epsilon)) {
+      const result = x.sub(y.mul(round(div)))
+      return bigNearlyEqual(result, round(result), config.epsilon)
+        ? round(result)
+        : result
+    } else {
+      return x.sub(y.mul(div.floor()))
+    }
+  }
 })

src/function/arithmetic/mod.js

@@ -1,22 +1,26 @@
+import { nearlyEqual as bigNearlyEqual } from '../../utils/bignumber/nearlyEqual.js'
+import { nearlyEqual } from '../../utils/number.js'
 import { factory } from '../../utils/factory.js'
 import { createMatAlgo02xDS0 } from '../../type/matrix/utils/matAlgo02xDS0.js'
 import { createMatAlgo03xDSf } from '../../type/matrix/utils/matAlgo03xDSf.js'
 import { createMatAlgo05xSfSf } from '../../type/matrix/utils/matAlgo05xSfSf.js'
 import { createMatAlgo11xS0s } from '../../type/matrix/utils/matAlgo11xS0s.js'
 import { createMatAlgo12xSfs } from '../../type/matrix/utils/matAlgo12xSfs.js'
-import { modNumber } from '../../plain/number/index.js'
+// import { modNumber } from '../../plain/number/index.js'
 import { createMatrixAlgorithmSuite } from '../../type/matrix/utils/matrixAlgorithmSuite.js'
 
 const name = 'mod'
 const dependencies = [
   'typed',
+  'config',
+  'round',
   'matrix',
   'equalScalar',
   'DenseMatrix',
   'concat'
 ]
 
-export const createMod = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, DenseMatrix, concat }) => {
+export const createMod = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, round, matrix, equalScalar, DenseMatrix, concat }) => {
   const matAlgo02xDS0 = createMatAlgo02xDS0({ typed, equalScalar })
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo05xSfSf = createMatAlgo05xSfSf({ typed, equalScalar })
@@ -62,13 +66,23 @@ export const createMod = /* #__PURE__ */ factory(name, dependencies, ({ typed, m
   return typed(
     name,
     {
-      'number, number': modNumber,
+      'number, number': _modNumber,
 
       'BigNumber, BigNumber': function (x, y) {
         if (y.isNeg()) {
           throw new Error('Cannot calculate mod for a negative divisor')
+        } else if (y.isZero()) {
+          return x
+        }
+        const div = x.div(y)
+        if (bigNearlyEqual(div, round(div), config.epsilon)) {
+          const result = x.sub(y.mul(round(div)))
+          return bigNearlyEqual(result, round(result), config.epsilon)
+            ? round(result)
+            : result
+        } else {
+          return x.sub(y.mul(div.floor()))
         }
-        return y.isZero() ? x : x.mod(y)
       },
 
       'Fraction, Fraction': function (x, y) {
@@ -87,4 +101,35 @@ export const createMod = /* #__PURE__ */ factory(name, dependencies, ({ typed, m
       sS: matAlgo12xSfs
     })
   )
+
+  /**
+ * Calculate the modulus of two numbers
+ * @param {number} x
+ * @param {number} y
+ * @returns {number} res
+ * @private
+ */
+  function _modNumber (x, y) {
+    if (y > 0) {
+      // We don't use JavaScript's % operator here as this doesn't work
+      // correctly for x < 0 and x === 0
+      // see https://en.wikipedia.org/wiki/Modulo_operation
+
+      // To ensure precision with float approximation 
+      const div = x / y
+      if (nearlyEqual(div, round(div), config.epsilon)) {
+        const result = x - y * round(div)
+        return nearlyEqual(result, round(Math.abs(result)), config.epsilon)
+          ? round(Math.abs(result))
+          : result
+      } else { // then Math.floor is precise
+        return x - y * Math.floor(x / y)
+      }
+    } else if (y === 0) {
+      return x
+    } else { // y < 0
+      // TODO: implement mod for a negative divisor
+      throw new Error('Cannot calculate mod for a negative divisor')
+    }
+  }
 })

test/unit-tests/function/arithmetic/mod.test.js

@@ -34,6 +34,12 @@ describe('mod', function () {
     approx.equal(mod(-5, 3), 1)
   })
 
+  it('should handle precise approximation of float approximation', function () {
+    assert.strictEqual(mod(0.1, 0.01), 0)
+    assert.strictEqual(mod(0.15, 0.05), 0)
+    assert.strictEqual(mod(1.23456789, 0.00000000001), 0)
+  })
+
   it('should throw an error if the divisor is negative', function () {
     assert.throws(function () { mod(10, -4) })
   })