Merge pull request #1 from bitwiseshiftleft/master

Update

1 files changed, 268 insertions, 66 deletions

diff --git a/core/bn.js b/core/bn.js
index 672fd08..01dd349 100644
--- a/core/bn.js
+++ b/core/bn.js
@@ -1,4 +1,5 @@
 /**
+ * @constructor
  * Constructs a new bignum from another bignum, a number or a hex string.
  */
 sjcl.bn = function(it) {
@@ -9,7 +10,7 @@ sjcl.bn.prototype = {
   radix: 24,
   maxMul: 8,
   _class: sjcl.bn,
-  
+
   copy: function() {
     return new this._class(this);
   },
@@ -18,17 +19,17 @@ sjcl.bn.prototype = {
    * Initializes this with it, either as a bn, a number, or a hex string.
    */
   initWith: function(it) {
-    var i=0, k, n, l;
+    var i=0, k;
     switch(typeof it) {
     case "object":
       this.limbs = it.limbs.slice(0);
       break;
-      
+
     case "number":
       this.limbs = [it];
       this.normalize();
       break;
-      
+
     case "string":
       it = it.replace(/^0x/, '');
       this.limbs = [];
@@ -59,14 +60,14 @@ sjcl.bn.prototype = {
     }
     return (difference === 0);
   },
-  
+
   /**
    * Get the i'th limb of this, zero if i is too large.
    */
   getLimb: function(i) {
     return (i >= this.limbs.length) ? 0 : this.limbs[i];
   },
-  
+
   /**
    * Constant time comparison function.
    * Returns 1 if this >= that, or zero otherwise.
@@ -83,7 +84,7 @@ sjcl.bn.prototype = {
     }
     return (greater | ~less) >>> 31;
   },
-  
+
   /**
    * Convert to a hex string.
    */
@@ -99,7 +100,7 @@ sjcl.bn.prototype = {
     }
     return "0x"+out;
   },
-  
+
   /** this += that.  Does not normalize. */
   addM: function(that) {
     if (typeof(that) !== "object") { that = new this._class(that); }
@@ -112,7 +113,7 @@ sjcl.bn.prototype = {
     }
     return this;
   },
-  
+
   /** this *= 2.  Requires normalized; ends up normalized. */
   doubleM: function() {
     var i, carry=0, tmp, r=this.radix, m=this.radixMask, l=this.limbs;
@@ -127,7 +128,7 @@ sjcl.bn.prototype = {
     }
     return this;
   },
-  
+
   /** this /= 2, rounded down.  Requires normalized; ends up normalized. */
   halveM: function() {
     var i, carry=0, tmp, r=this.radix, l=this.limbs;
@@ -154,14 +155,21 @@ sjcl.bn.prototype = {
     }
     return this;
   },
-  
+
   mod: function(that) {
+    var neg = !this.greaterEquals(new sjcl.bn(0));
+
     that = new sjcl.bn(that).normalize(); // copy before we begin
     var out = new sjcl.bn(this).normalize(), ci=0;
-    
+
+    if (neg) out = (new sjcl.bn(0)).subM(out).normalize();
+
     for (; out.greaterEquals(that); ci++) {
       that.doubleM();
     }
+
+    if (neg) out = that.sub(out).normalize();
+
     for (; ci > 0; ci--) {
       that.halveM();
       if (out.greaterEquals(that)) {
@@ -170,15 +178,15 @@ sjcl.bn.prototype = {
     }
     return out.trim();
   },
-  
+
   /** return inverse mod prime p.  p must be odd. Binary extended Euclidean algorithm mod p. */
   inverseMod: function(p) {
     var a = new sjcl.bn(1), b = new sjcl.bn(0), x = new sjcl.bn(this), y = new sjcl.bn(p), tmp, i, nz=1;
-    
+
     if (!(p.limbs[0] & 1)) {
       throw (new sjcl.exception.invalid("inverseMod: p must be odd"));
     }
-    
+
     // invariant: y is odd
     do {
       if (x.limbs[0] & 1) {
@@ -189,13 +197,13 @@ sjcl.bn.prototype = {
         }
         x.subM(y);
         x.normalize();
-        
+
         if (!a.greaterEquals(b)) {
           a.addM(p);
         }
         a.subM(b);
       }
-      
+
       // cut everything in half
       x.halveM();
       if (a.limbs[0] & 1) {
@@ -203,20 +211,20 @@ sjcl.bn.prototype = {
       }
       a.normalize();
       a.halveM();
-      
+
       // check for termination: x ?= 0
       for (i=nz=0; i<x.limbs.length; i++) {
         nz |= x.limbs[i];
       }
     } while(nz);
-    
+
     if (!y.equals(1)) {
       throw (new sjcl.exception.invalid("inverseMod: p and x must be relatively prime"));
     }
-    
+
     return b;
   },
-  
+
   /** this + that.  Does not normalize. */
   add: function(that) {
     return this.copy().addM(that);
@@ -226,7 +234,7 @@ sjcl.bn.prototype = {
   sub: function(that) {
     return this.copy().subM(that);
   },
-  
+
   /** this * that.  Normalizes and reduces. */
   mul: function(that) {
     if (typeof(that) === "number") { that = new this._class(that); }
@@ -240,7 +248,7 @@ sjcl.bn.prototype = {
       for (j=0; j<bl; j++) {
         c[i+j] += ai * b[j];
       }
-     
+
       if (!--ii) {
         ii = this.maxMul;
         out.cnormalize();
@@ -256,22 +264,18 @@ sjcl.bn.prototype = {
 
   /** this ^ n.  Uses square-and-multiply.  Normalizes and reduces. */
   power: function(l) {
-    if (typeof(l) === "number") {
-      l = [l];
-    } else if (l.limbs !== undefined) {
-      l = l.normalize().limbs;
-    }
+    l = new sjcl.bn(l).normalize().trim().limbs;
     var i, j, out = new this._class(1), pow = this;
 
     for (i=0; i<l.length; i++) {
       for (j=0; j<this.radix; j++) {
-        if (l[i] & (1<<j)) {
-          out = out.mul(pow);
-        }
+        if (l[i] & (1<<j)) { out = out.mul(pow); }
+        if (i == (l.length - 1) && l[i]>>(j + 1) == 0) { break; }
+
         pow = pow.square();
       }
     }
-    
+
     return out;
   },
 
@@ -282,14 +286,184 @@ sjcl.bn.prototype = {
 
   /** this ^ x mod N */
   powermod: function(x, N) {
-    var result = new sjcl.bn(1), a = new sjcl.bn(this), k = new sjcl.bn(x);
-    while (true) {
-      if (k.limbs[0] & 1) { result = result.mulmod(a, N); }
-      k.halveM();
-      if (k.equals(0)) { break; }
-      a = a.mulmod(a, N);
-    }
-    return result.normalize().reduce();
+    x = new sjcl.bn(x);
+    N = new sjcl.bn(N);
+
+    // Jump to montpowermod if possible.
+    if ((N.limbs[0] & 1) == 1) {
+      var montOut = this.montpowermod(x, N);
+
+      if (montOut != false) { return montOut; } // else go to slow powermod
+    }
+
+    var i, j, l = x.normalize().trim().limbs, out = new this._class(1), pow = this;
+
+    for (i=0; i<l.length; i++) {
+      for (j=0; j<this.radix; j++) {
+        if (l[i] & (1<<j)) { out = out.mulmod(pow, N); }
+        if (i == (l.length - 1) && l[i]>>(j + 1) == 0) { break; }
+
+        pow = pow.mulmod(pow, N);
+      }
+    }
+
+    return out;
+  },
+
+  /** this ^ x mod N with Montomery reduction */
+  montpowermod: function(x, N) {
+    x = new sjcl.bn(x).normalize().trim();
+    N = new sjcl.bn(N);
+
+    var i, j,
+      radix = this.radix,
+      out = new this._class(1),
+      pow = this.copy();
+
+    // Generate R as a cap of N.
+    var R, s, wind, bitsize = x.bitLength();
+
+    R = new sjcl.bn({
+      limbs: N.copy().normalize().trim().limbs.map(function() { return 0; })
+    });
+
+    for (s = this.radix; s > 0; s--) {
+      if (((N.limbs[N.limbs.length - 1] >> s) & 1) == 1) {
+        R.limbs[R.limbs.length - 1] = 1 << s;
+        break;
+      }
+    }
+
+    // Calculate window size as a function of the exponent's size.
+    if (bitsize == 0) {
+      return this;
+    } else if (bitsize < 18)  {
+      wind = 1;
+    } else if (bitsize < 48)  {
+      wind = 3;
+    } else if (bitsize < 144) {
+      wind = 4;
+    } else if (bitsize < 768) {
+      wind = 5;
+    } else {
+      wind = 6;
+    }
+
+    // Find R' and N' such that R * R' - N * N' = 1.
+    var RR = R.copy(), NN = N.copy(), RP = new sjcl.bn(1), NP = new sjcl.bn(0), RT = R.copy();
+
+    while (RT.greaterEquals(1)) {
+      RT.halveM();
+
+      if ((RP.limbs[0] & 1) == 0) {
+        RP.halveM();
+        NP.halveM();
+      } else {
+        RP.addM(NN);
+        RP.halveM();
+
+        NP.halveM();
+        NP.addM(RR);
+      }
+    }
+
+    RP = RP.normalize();
+    NP = NP.normalize();
+
+    RR.doubleM()
+    var R2 = RR.mulmod(RR, N);
+
+    // Check whether the invariant holds.
+    // If it doesn't, we can't use Montgomery reduction on this modulus.
+    if (!RR.mul(RP).sub(N.mul(NP)).equals(1)) {
+      return false;
+    }
+
+    var montIn = function(c) { return montMul(c, R2) },
+    montMul = function(a, b) {
+      // Standard Montgomery reduction
+      var k, carry, ab, right, abBar, mask = (1 << (s + 1)) - 1;
+
+      ab = a.mul(b);
+
+      right = ab.mul(NP);
+      right.limbs = right.limbs.slice(0, R.limbs.length);
+
+      if (right.limbs.length == R.limbs.length) {
+        right.limbs[R.limbs.length - 1] &= mask;
+      }
+
+      right = right.mul(N);
+
+      abBar = ab.add(right).normalize().trim();
+      abBar.limbs = abBar.limbs.slice(R.limbs.length - 1);
+
+      // Division.  Equivelent to calling *.halveM() s times.
+      for (k=0; k < abBar.limbs.length; k++) {
+        if (k > 0) {
+          abBar.limbs[k - 1] |= (abBar.limbs[k] & mask) << (radix - s - 1)
+        }
+
+        abBar.limbs[k] = abBar.limbs[k] >> (s + 1)
+      }
+
+      if (abBar.greaterEquals(N)) {
+        abBar.subM(N)
+      }
+
+      return abBar;
+    },
+    montOut = function(c) { return montMul(c, 1); };
+
+    pow = montIn(pow);
+    out = montIn(out);
+
+    // Sliding-Window Exponentiation (HAC 14.85)
+    var h, precomp = {}, cap = (1 << (wind - 1)) - 1;
+
+    precomp[1] = pow.copy();
+    precomp[2] = montMul(pow, pow);
+
+    for (h=1; h<=cap; h++) {
+      precomp[(2 * h) + 1] = montMul(precomp[(2 * h) - 1], precomp[2]);
+    }
+
+    var getBit = function(exp, i) { // Gets ith bit of exp.
+      var off = i % exp.radix;
+
+      return (exp.limbs[Math.floor(i / exp.radix)] & (1 << off)) >> off;
+    }
+
+    for (i = x.bitLength() - 1; i >= 0; ) {
+      if (getBit(x, i) == 0) {
+        // If the next bit is zero:
+        //   Square, move forward one bit.
+        out = montMul(out, out)
+        i = i - 1;
+      } else {
+        // If the next bit is one:
+        //   Find the longest sequence of bits after this one, less than `wind`
+        //   bits long, that ends with a 1.  Convert the sequence into an
+        //   integer and look up the pre-computed value to add.
+        var l = i - wind + 1;
+
+        while (getBit(x, l) == 0) {
+          l++;
+        }
+
+        var indx = 0;
+        for (j = l; j <= i; j++) {
+          indx += getBit(x, j) << (j - l)
+          out = montMul(out, out);
+        }
+
+        out = montMul(out, precomp[indx])
+
+        i = l - 1;
+      }
+    }
+
+    return montOut(out);
   },
 
   trim: function() {
@@ -300,7 +474,7 @@ sjcl.bn.prototype = {
     l.push(p);
     return this;
   },
-  
+
   /** Reduce mod a modulus.  Stubbed for subclassing. */
   reduce: function() {
     return this;
@@ -310,7 +484,7 @@ sjcl.bn.prototype = {
   fullReduce: function() {
     return this.normalize();
   },
-  
+
   /** Propagate carries. */
   normalize: function() {
     var carry=0, i, pv = this.placeVal, ipv = this.ipv, l, m, limbs = this.limbs, ll = limbs.length, mask = this.radixMask;
@@ -320,7 +494,10 @@ sjcl.bn.prototype = {
       carry = (l-m)*ipv;
     }
     if (carry === -1) {
-      limbs[i-1] -= this.placeVal;
+      limbs[i-1] -= pv;
+    }
+    while (limbs.length > 0 && limbs[limbs.length-1] === 0) {
+      limbs.pop();
     }
     return this;
   },
@@ -336,35 +513,39 @@ sjcl.bn.prototype = {
     limbs[i] += carry;
     return this;
   },
-  
+
   /** Serialize to a bit array */
   toBits: function(len) {
     this.fullReduce();
-    len = len || this.exponent || this.limbs.length * this.radix;
+    len = len || this.exponent || this.bitLength();
     var i = Math.floor((len-1)/24), w=sjcl.bitArray, e = (len + 7 & -8) % this.radix || this.radix,
         out = [w.partial(e, this.getLimb(i))];
     for (i--; i >= 0; i--) {
-      out = w.concat(out, [w.partial(this.radix, this.getLimb(i))]);
+      out = w.concat(out, [w.partial(Math.min(this.radix,len), this.getLimb(i))]);
+      len -= this.radix;
     }
     return out;
   },
-  
+
   /** Return the length in bits, rounded up to the nearest byte. */
   bitLength: function() {
     this.fullReduce();
     var out = this.radix * (this.limbs.length - 1),
         b = this.limbs[this.limbs.length - 1];
-    for (; b; b >>= 1) {
+    for (; b; b >>>= 1) {
       out ++;
     }
     return out+7 & -8;
   }
 };
 
+/** @memberOf sjcl.bn
+* @this { sjcl.bn }
+*/
 sjcl.bn.fromBits = function(bits) {
   var Class = this, out = new Class(), words=[], w=sjcl.bitArray, t = this.prototype,
       l = Math.min(this.bitLength || 0x100000000, w.bitLength(bits)), e = l % t.radix || t.radix;
-  
+
   words[0] = w.extract(bits, 0, e);
   for (; e < l; e += t.radix) {
     words.unshift(w.extract(bits, e, t.radix));
@@ -384,6 +565,9 @@ sjcl.bn.prototype.radixMask = (1 << sjcl.bn.prototype.radix) - 1;
  * i.e. a prime of the form 2^e + sum(a * 2^b),where the sum is negative and sparse.
  */
 sjcl.bn.pseudoMersennePrime = function(exponent, coeff) {
+  /** @constructor
+  * @private
+  */
   function p(it) {
     this.initWith(it);
     /*if (this.limbs[this.modOffset]) {
@@ -401,7 +585,7 @@ sjcl.bn.pseudoMersennePrime = function(exponent, coeff) {
   ppr.fullOffset = [];
   ppr.fullFactor = [];
   ppr.modulus = p.modulus = new sjcl.bn(Math.pow(2,exponent));
-  
+
   ppr.fullMask = 0|-Math.pow(2, exponent % ppr.radix);
 
   for (i=0; i<coeff.length; i++) {
@@ -415,9 +599,12 @@ sjcl.bn.pseudoMersennePrime = function(exponent, coeff) {
   ppr._class = p;
   ppr.modulus.cnormalize();
 
-  /** Approximate reduction mod p.  May leave a number which is negative or slightly larger than p. */
+  /** Approximate reduction mod p.  May leave a number which is negative or slightly larger than p.
+   * @memberof sjcl.bn
+   * @this { sjcl.bn }
+   */
   ppr.reduce = function() {
-    var i, k, l, mo = this.modOffset, limbs = this.limbs, aff, off = this.offset, ol = this.offset.length, fac = this.factor, ll;
+    var i, k, l, mo = this.modOffset, limbs = this.limbs, off = this.offset, ol = this.offset.length, fac = this.factor, ll;
 
     i = this.minOffset;
     while (limbs.length > mo) {
@@ -426,7 +613,7 @@ sjcl.bn.pseudoMersennePrime = function(exponent, coeff) {
       for (k=0; k<ol; k++) {
         limbs[ll+off[k]] -= fac[k] * l;
       }
-      
+
       i--;
       if (!i) {
         limbs.push(0);
@@ -438,7 +625,10 @@ sjcl.bn.pseudoMersennePrime = function(exponent, coeff) {
 
     return this;
   };
-  
+
+  /** @memberof sjcl.bn
+  * @this { sjcl.bn }
+  */
   ppr._strongReduce = (ppr.fullMask === -1) ? ppr.reduce : function() {
     var limbs = this.limbs, i = limbs.length - 1, k, l;
     this.reduce();
@@ -452,11 +642,14 @@ sjcl.bn.pseudoMersennePrime = function(exponent, coeff) {
     }
   };
 
-  /** mostly constant-time, very expensive full reduction. */
+  /** mostly constant-time, very expensive full reduction.
+   * @memberof sjcl.bn
+   * @this { sjcl.bn }
+   */
   ppr.fullReduce = function() {
     var greater, i;
     // massively above the modulus, may be negative
-    
+
     this._strongReduce();
     // less than twice the modulus, may be negative
 
@@ -464,7 +657,7 @@ sjcl.bn.pseudoMersennePrime = function(exponent, coeff) {
     this.addM(this.modulus);
     this.normalize();
     // probably 2-3x the modulus
-    
+
     this._strongReduce();
     // less than the power of 2.  still may be more than
     // the modulus
@@ -473,7 +666,7 @@ sjcl.bn.pseudoMersennePrime = function(exponent, coeff) {
     for (i=this.limbs.length; i<this.modOffset; i++) {
       this.limbs[i] = 0;
     }
-    
+
     // constant-time subtract modulus
     greater = this.greaterEquals(this.modulus);
     for (i=0; i<this.limbs.length; i++) {
@@ -484,6 +677,10 @@ sjcl.bn.pseudoMersennePrime = function(exponent, coeff) {
     return this;
   };
 
+
+  /** @memberof sjcl.bn
+  * @this { sjcl.bn }
+  */
   ppr.inverse = function() {
     return (this.power(this.modulus.sub(2)));
   };
@@ -494,18 +691,24 @@ sjcl.bn.pseudoMersennePrime = function(exponent, coeff) {
 };
 
 // a small Mersenne prime
+var sbp = sjcl.bn.pseudoMersennePrime;
 sjcl.bn.prime = {
-  p127: sjcl.bn.pseudoMersennePrime(127, [[0,-1]]),
+  p127: sbp(127, [[0,-1]]),
 
   // Bernstein's prime for Curve25519
-  p25519: sjcl.bn.pseudoMersennePrime(255, [[0,-19]]),
+  p25519: sbp(255, [[0,-19]]),
+
+  // Koblitz primes
+  p192k: sbp(192, [[32,-1],[12,-1],[8,-1],[7,-1],[6,-1],[3,-1],[0,-1]]),
+  p224k: sbp(224, [[32,-1],[12,-1],[11,-1],[9,-1],[7,-1],[4,-1],[1,-1],[0,-1]]),
+  p256k: sbp(256, [[32,-1],[9,-1],[8,-1],[7,-1],[6,-1],[4,-1],[0,-1]]),
 
   // NIST primes
-  p192: sjcl.bn.pseudoMersennePrime(192, [[0,-1],[64,-1]]),
-  p224: sjcl.bn.pseudoMersennePrime(224, [[0,1],[96,-1]]),
-  p256: sjcl.bn.pseudoMersennePrime(256, [[0,-1],[96,1],[192,1],[224,-1]]),
-  p384: sjcl.bn.pseudoMersennePrime(384, [[0,-1],[32,1],[96,-1],[128,-1]]),
-  p521: sjcl.bn.pseudoMersennePrime(521, [[0,-1]])
+  p192: sbp(192, [[0,-1],[64,-1]]),
+  p224: sbp(224, [[0,1],[96,-1]]),
+  p256: sbp(256, [[0,-1],[96,1],[192,1],[224,-1]]),
+  p384: sbp(384, [[0,-1],[32,1],[96,-1],[128,-1]]),
+  p521: sbp(521, [[0,-1]])
 };
 
 sjcl.bn.random = function(modulus, paranoia) {
@@ -531,4 +734,3 @@ sjcl.bn.random = function(modulus, paranoia) {
     }
   }
 };
-