From 91141bbbb6612e0ef977db2cf2f5f369e26e0f3a Mon Sep 17 00:00:00 2001
From: Viktor <vic99999@yandex.ru>
Date: Sat, 10 Apr 2021 10:14:16 +0200
Subject: [PATCH] =?UTF-8?q?Miller=E2=80=93Rabin=20primality=20test=20(#213?=
 =?UTF-8?q?9)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

* Miller–Rabin primality test

* add tests for some big numbers

* Update isPrime.js

* Update isPrime.js

Co-authored-by: Jos de Jong <wjosdejong@gmail.com>
---
 src/function/utils/isPrime.js                 | 55 ++++++++++++++++++-
 .../unit-tests/function/utils/isPrime.test.js |  5 ++
 2 files changed, 57 insertions(+), 3 deletions(-)

diff --git a/src/function/utils/isPrime.js b/src/function/utils/isPrime.js
index bffe2ca102..f03bdddbf0 100644
--- a/src/function/utils/isPrime.js
+++ b/src/function/utils/isPrime.js
@@ -58,13 +58,62 @@ export const createIsPrime = /* #__PURE__ */ factory(name, dependencies, ({ type
       }
       if (n.lte(3)) return n.gt(1)
       if (n.mod(2).eq(0) || n.mod(3).eq(0)) return false
+      if (n.lt(Math.pow(2, 32))) {
+        const x = n.toNumber()
+        for (let i = 5; i * i <= x; i += 6) {
+          if (x % i === 0 || x % (i + 2) === 0) {
+            return false
+          }
+        }
+        return true
+      }
 
-      for (let i = 5; n.gte(i * i); i += 6) {
-        if (n.mod(i).eq(0) || n.mod(i + 2).eq(0)) {
-          return false
+      function modPow (base, exponent, modulus) {
+        // exponent can be huge, use non-recursive variant
+        let accumulator = 1
+        while (!exponent.eq(0)) {
+          if (exponent.mod(2).eq(0)) {
+            exponent = exponent.div(2)
+            base = base.mul(base).mod(modulus)
+          } else {
+            exponent = exponent.sub(1)
+            accumulator = base.mul(accumulator).mod(modulus)
+          }
         }
+        return accumulator
       }
 
+      // https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Deterministic_variants
+      const Decimal = n.constructor.clone({ precision: n.toFixed(0).length * 2 })
+      n = new Decimal(n)
+      let r = 0
+      let d = n.sub(1)
+      while (d.mod(2).eq(0)) {
+        d = d.div(2)
+        r += 1
+      }
+      let bases = null
+      // https://en.wikipedia.org/wiki/Miller–Rabin_primality_test#Testing_against_small_sets_of_bases
+      if (n.lt('3317044064679887385961981')) {
+        bases = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41].filter(x => x < n)
+      } else {
+        const max = Math.min(n.toNumber() - 2, Math.floor(2 * Math.pow(n.toFixed(0).length * Math.log(10), 2)))
+        bases = []
+        for (let i = 2; i <= max; i += 1) {
+          bases.push(max)
+        }
+      }
+      for (let i = 0; i < bases.length; i += 1) {
+        const a = bases[i]
+        const adn = modPow(n.sub(n).add(a), d, n)
+        if (!adn.eq(1)) {
+          for (let i = 0, x = adn; !x.eq(n.sub(1)); i += 1, x = x.mul(x).mod(n)) {
+            if (i === r - 1) {
+              return false
+            }
+          }
+        }
+      }
       return true
     },
 
diff --git a/test/unit-tests/function/utils/isPrime.test.js b/test/unit-tests/function/utils/isPrime.test.js
index 004aac088a..02d60e8417 100644
--- a/test/unit-tests/function/utils/isPrime.test.js
+++ b/test/unit-tests/function/utils/isPrime.test.js
@@ -44,4 +44,9 @@ describe('isPrime', function () {
     assert.throws(function () { isPrime(new Date()) }, /TypeError: Unexpected type of argument/)
     assert.throws(function () { isPrime({}) }, /TypeError: Unexpected type of argument/)
   })
+
+  it('should work fast for huge values', function () {
+    assert.strictEqual(isPrime(bignumber('2305843009213693951')), true)
+    assert.strictEqual(isPrime(bignumber('230584300921369395')), false)
+  })
 })