Use GMP & BCMath native functions for powerMod()

src/BigInteger.php

@@ -468,33 +468,9 @@ public function powerMod($exp, $mod) : BigInteger
             throw DivisionByZeroException::divisionByZero();
         }
 
-        // Special case: the algorithm below fails with 0 power 0 mod 1 (returns 1 instead of 0)
-        if ($this->value === '0' && $exp->value === '0' && $mod->value === '1') {
-            return BigInteger::zero();
-        }
-
-        // Special case: the algorithm below fails with power 0 mod 1 (returns 1 instead of 0)
-        if ($exp->value === '0' && $mod->value === '1') {
-            return BigInteger::zero();
-        }
-
-        $x = $this;
+        $result = Calculator::get()->powmod($this->value, $exp->value, $mod->value);
 
-        $res = BigInteger::one();
-
-        // numbers are positive, so we can use remainder() instead of mod(), this is less expensive right now
-        $x = $x->remainder($mod);
-
-        while ($exp->isPositive()) {
-            if ($exp->remainder(2)->isEqualTo(1)) {
-                $res = $res->multipliedBy($x)->remainder($mod);
-            }
-
-            $exp = $exp->quotient(2);
-            $x = $x->multipliedBy($x)->remainder($mod);
-        }
-
-        return $res;
+        return new BigInteger($result);
     }
 
     /**

src/Internal/Calculator.php

@@ -237,13 +237,24 @@ abstract public function divQR(string $a, string $b) : array;
     /**
      * Exponentiates a number.
      *
-     * @param string $a The base.
+     * @param string $a The base number.
      * @param int    $e The exponent, validated as an integer between 0 and MAX_POWER.
      *
      * @return string The power.
      */
     abstract public function pow(string $a, int $e) : string;
 
+    /**
+     * Raises a number into power with modulo.
+     *
+     * @param string $base The base number; must be positive or zero.
+     * @param string $exp  The exponent; must be positive or zero.
+     * @param string $mod  The modulo; must be strictly positive.
+     *
+     * @return string The power.
+     */
+    abstract function powmod(string $base, string $exp, string $mod) : string;
+
     /**
      * Returns the greatest common divisor of the two numbers.
      *

src/Internal/Calculator/BcMathCalculator.php

@@ -72,6 +72,14 @@ public function pow(string $a, int $e) : string
         return \bcpow($a, (string) $e, 0);
     }
 
+    /**
+     * {@inheritdoc}
+     */
+    public function powmod(string $base, string $exp, string $mod) : string
+    {
+        return \bcpowmod($base, $exp, $mod, 0);
+    }
+
     /**
      * {@inheritDoc}
      */

src/Internal/Calculator/GmpCalculator.php

@@ -74,6 +74,14 @@ public function pow(string $a, int $e) : string
         return \gmp_strval(\gmp_pow($a, $e));
     }
 
+    /**
+     * {@inheritdoc}
+     */
+    public function powmod(string $base, string $exp, string $mod) : string
+    {
+        return \gmp_strval(\gmp_powm($base, $exp, $mod));
+    }
+
     /**
      * {@inheritdoc}
      */

src/Internal/Calculator/NativeCalculator.php

@@ -228,6 +228,40 @@ public function pow(string $a, int $e) : string
         return $result;
     }
 
+    /**
+     * {@inheritdoc}
+     */
+    public function powmod(string $base, string $exp, string $mod) : string
+    {
+        // special case: the algorithm below fails with 0 power 0 mod 1 (returns 1 instead of 0)
+        if ($base === '0' && $exp === '0' && $mod === '1') {
+            return '0';
+        }
+
+        // special case: the algorithm below fails with power 0 mod 1 (returns 1 instead of 0)
+        if ($exp === '0' && $mod === '1') {
+            return '0';
+        }
+
+        $x = $base;
+
+        $res = '1';
+
+        // numbers are positive, so we can use remainder instead of modulo
+        $x = $this->divR($x, $mod);
+
+        while ($exp !== '0') {
+            if (in_array($exp[-1], ['1', '3', '5', '7', '9'])) { // odd
+                $res = $this->divR($this->mul($res, $x), $mod);
+            }
+
+            $exp = $this->divQ($exp, '2');
+            $x = $this->divR($this->mul($x, $x), $mod);
+        }
+
+        return $res;
+    }
+
     /**
      * Adapted from https://cp-algorithms.com/num_methods/roots_newton.html
      *