BigInt tweaks and fixes

jmat_bigint.js

@@ -520,6 +520,7 @@ Jmat.BigInt.rshift = function(a, b) {
   }
 
   if(a.minus) result = result.neg();
+  B.stripInPlace_(result.a);
   return result;
 };
 Jmat.BigInt.prototype.rshift = function(b) {
@@ -895,9 +896,8 @@ Jmat.BigInt.prototype.comparer = function(b) {
   var sign = this.getSign();
   b = Math.abs(b);
 
-  var n = Jmat.BigInt.getNumDigits(this);
   var r = 0;
-  for(var i = 0; i < n && r <= b; i++) {
+  for(var i = 0; i < this.a.length && r <= b; i++) {
     r *= this.radix;
     r += this.a[i];
   }
@@ -1264,12 +1264,25 @@ Jmat.BigInt.primeCache_ = [];
 // Returns prime factors in array.
 // Will not factorize if the problem is too difficult (only uses simple checks), last element of result is 0 then to indicate the error.
 // Result may be probabilistic since a probabilistic prime test is used.
+// For inputs smaller than 2, includes the non-composite factors -1, 0 and 1 in the result
 Jmat.BigInt.factorize = function(a) {
   var B = Jmat.BigInt;
 
-  // avoid infinite loops
-  if(a.eqr(0)) return [B(0)];
-  if(a.eqr(1)) return [B(1)];
+  var result = [];
+  if(a.minus) {
+    a = B.neg(a);
+    result.push(B(-1));
+  }
+  if(a.lter(2)) {
+    if(result.length == 0 || !a.eqr(1)) result.push(a); // return [0] if x is 0, [1] if x is 1. Also avoids infinite loops.
+    return result;
+  }
+
+  if (a.ltr(Jmat.Real.BIGGESTJSINT)) {
+    var r = Jmat.Real.factorize(a);
+    for(var i = 0; i < r.length; i++) result.push(B(r[i]));
+    return result;
+  }
 
   // returns a factor, or a itself if end reached (a is prime), or 0 if undetermined because the problem is too hard
   var f = function(a) {
@@ -1303,7 +1316,6 @@ Jmat.BigInt.factorize = function(a) {
     return B(0); // Not found
   }
 
-  var result = [];
   for(;;) {
     var b = f(a);
     if(b.eqr(0)) {
@@ -1852,10 +1864,11 @@ Jmat.BigInt.isPrimeSimple = function(n) {
 };
 
 //do rounds of Miller-Rabin primality test
-//base = the potential witnesses to try as an array, e.g. [2, 3]. The members of the array may be either regular JS number or BigInt
+//opt_base = the potential witnesses to try as an array, e.g. [2, 3]. The members of the array may be either regular JS number or BigInt. Optional parameter, if not given, chosen automatically (random)
 //requires that n is big enough, at least 3. First use Jmat.BigInt.isPrimeSimple, and only use this function if that returns -1.
-Jmat.BigInt.isPrimeMillerRabin = function(n, base) {
+Jmat.BigInt.isPrimeMillerRabin = function(n, opt_base) {
   var B = Jmat.BigInt;
+  var base = opt_base || B.chooseMillerRabinBase_(n);
 
   // choose s and odd d such that n = 2^s * d
   var d = n.divr(2);
@@ -1914,11 +1927,11 @@ Jmat.BigInt.chooseMillerRabinBase_ = function(n) {
 //Is *probably* prime at least.
 //If you wish to control amount of miller rabin rounds and bases, use isPrimeMillerRabin directly with own bases (first use isPrimeSimple).
 Jmat.BigInt.isPrime = function(n) {
+  if(n.ltr(Jmat.Real.BIGGESTJSPRIME)) return Jmat.Real.isPrime(n.toInt());
   var init = Jmat.BigInt.isPrimeSimple(n);
   if(init != -1) return !!init;
 
-  var base = Jmat.BigInt.chooseMillerRabinBase_(n);
-  return Jmat.BigInt.isPrimeMillerRabin(n, base);
+  return Jmat.BigInt.isPrimeMillerRabin(n);
 };
 
 /*
@@ -1976,6 +1989,7 @@ Jmat.BigInt.getNumDigits = function(x) {
 };
 
 //strips array a (removing leading zeroes unless it is zero), modifying the object
+//NOTE: do not give a BigInt, a must be an array
 Jmat.BigInt.stripInPlace_ = function(a) {
   while(a[0] == 0 && a.length > 1) a.shift(); //JS array shift function (a is array, not BigInt)
 };
@@ -1994,7 +2008,7 @@ Jmat.BigInt.copystrip_ = function(a) {
 };
 
 //makes copy if needed object needs to be altered, or original if it doesn't have to be stripped
-//a is array
+//NOTE: do not give a BigInt, a must be an array
 Jmat.BigInt.maybecopystrip_ = function(a, allowmodify) {
   if(a.length <= 1 || a[0] != 0) return a;
   if(!allowmodify) a = Jmat.BigInt.copyarray_(a);

jmat_matrix.js

@@ -2778,8 +2778,10 @@ Jmat.Matrix.rref = function(a) {
     var p = pivots[k];
     // make corresponding elements of row above zero using row operations
     for(var y = k - 1; y >= 0; y--) {
-      submul(a, p + 1, k, a.e[y][p], y);
-      a.e[y][p] = C(0); // make extra-sure it's 0, avoid numerical imprecision
+      if(!(a.e[y][p].eqr(0))) {
+        submul(a, p + 1, k, a.e[y][p], y);
+        a.e[y][p] = C(0); // make extra-sure it's 0, avoid numerical imprecision
+      }
     }
   }
 

jmat_real.js

@@ -70,6 +70,7 @@ Jmat.Real.SQRTPI = Math.sqrt(Math.PI); // gamma(0.5)
 Jmat.Real.EM = 0.57721566490153286060; // Euler-Mascheroni constant
 Jmat.Real.APERY = 1.2020569; // Apery's constant, zeta(3)
 Jmat.Real.BIGGESTJSINT = 9007199254740992; // largest number that JS (float64) can represent as integer: 2^53, 0x20000000000000, 9007199254740992
+Jmat.Real.BIGGESTJSPRIME = 9007199254740881; // largest prime number that JS (float64) can represent as integer, that is, the biggest prime smaller than Jmat.Real.BIGGESTJSINT.
 
 ////////////////////////////////////////////////////////////////////////////////
 // Categories
@@ -455,12 +456,19 @@ Jmat.Real.smallestPrimeFactor = function(x) {
   return x;
 };
 
-//factorize: returns prime factors as array of real integers, sorted from smallest to largest. x must be real non-negative integer.
+//factorize: returns prime factors as array of real integers, sorted from smallest to largest. x must be integer.
 Jmat.Real.factorize = function(x) {
   if(x > Jmat.Real.BIGGESTJSINT) return undefined; //too large for the floating point's integer precision, will cause crash
-  if(x <= 2) return [x]; // return [0] if x is 0, [1] if x is 1
   var x = Math.round(x);
   var result = [];
+  if(x < 0) {
+    x = -x;
+    result.push(-1);
+  }
+  if(x <= 2) {
+    if(result.length == 0 || x != 1) result.push(x); // return [0] if x is 0, [1] if x is 1
+    return result;
+  }
   for(;;) {
     if(x < 1) break;
     var y = Jmat.Real.smallestPrimeFactor(x);
@@ -657,6 +665,7 @@ Jmat.Real.pascal_triangle = function(n, p) {
 //greatest common divisor
 Jmat.Real.gcd = function(x, y) {
   if(!Jmat.Real.isInt(x) || !Jmat.Real.isInt(y)) return NaN; //prevents infinite loop if both x and y are NaN. Also, reals are not supported here.
+  if(Math.abs(x) > Jmat.Real.BIGGESTJSINT || Math.abs(y) > Jmat.Real.BIGGESTJSINT) return NaN; // does not work above JS integer precision
  //Euclid's algorithm
  for(;;) {
    if(y == 0) return Math.abs(x); //if x or y are negative, the result is still positive by the definition