Provide ia32 intrinsics for bigint add, sub, quotient digit estimation, and

runtime/lib/bigint.dart

 // Copyright (c) 2014, the Dart project authors.  Please see the AUTHORS file
 // for details. All rights reserved. Use of this source code is governed by a
 // BSD-style license that can be found in the LICENSE file.
 
 // Copyright 2009 The Go Authors. All rights reserved.
 // Use of this source code is governed by a BSD-style
 // license that can be found in the LICENSE file.
 
 /*
  * Copyright (c) 2003-2005  Tom Wu
@@ skipping 394 matching lines @@
         r = -r;
       }
     } else if (_neg) {
       r = -1;
     } else {
       r = 1;
     }
     return r;
   }
 
+  // r_digits[0..used] = digits[0..used-1] + a_digits[0..a_used-1].
+  // used >= a_used > 0.
+  static void _add(Uint32List digits, int used,
+                   Uint32List a_digits, int a_used,
+                   Uint32List r_digits) {
+    var c = 0;
+    for (var i = 0; i < a_used; i++) {
+      c += digits[i] + a_digits[i];
+      r_digits[i] = c & DIGIT_MASK;
+      c >>= DIGIT_BITS;
+    }
+    for (var i = a_used; i < used; i++) {
+      c += digits[i];
+      r_digits[i] = c & DIGIT_MASK;
+      c >>= DIGIT_BITS;
+    }
+    r_digits[used] = c;
+  }
+
+  // r_digits[0..used-1] = digits[0..used-1] - a_digits[0..a_used-1].
+  // used >= a_used > 0.
+  static void _sub(Uint32List digits, int used,
+                   Uint32List a_digits, int a_used,
+                   Uint32List r_digits) {
+    var c = 0;
+    for (var i = 0; i < a_used; i++) {
+      c += digits[i] - a_digits[i];
+      r_digits[i] = c & DIGIT_MASK;
+      c >>= DIGIT_BITS;
+    }
+    for (var i = a_used; i < used; i++) {
+      c += digits[i];
+      r_digits[i] = c & DIGIT_MASK;
+      c >>= DIGIT_BITS;
+    }
+  }
+
   // r = abs(this) + abs(a).
   void _absAddTo(_Bigint a, _Bigint r) {
     var used = _used;
     var a_used = a._used;
     if (used < a_used) {
       a._absAddTo(this, r);
       return;
     }
     if (used == 0) {
       // Set r to 0.
       r._neg = false;
       r._used = 0;
       return;
     }
     if (a_used == 0) {
       _copyTo(r);
       return;
     }
     r._ensureLength(used + 1);
-    var digits = _digits;
-    var a_digits = a._digits;
-    var r_digits = r._digits;
-    var c = 0;
-    for (var i = 0; i < a_used; i++) {
-      c += digits[i] + a_digits[i];
-      r_digits[i] = c & DIGIT_MASK;
-      c >>= DIGIT_BITS;
-    }
-    for (var i = a_used; i < used; i++) {
-      c += digits[i];
-      r_digits[i] = c & DIGIT_MASK;
-      c >>= DIGIT_BITS;
-    }
-    r_digits[used] = c;
+    _add(_digits, used, a._digits, a_used, r._digits);
     r._used = used + 1;
     r._clamp();
   }
 
   // r = abs(this) - abs(a), with abs(this) >= abs(a).
   void _absSubTo(_Bigint a, _Bigint r) {
     assert(_absCompareTo(a) >= 0);
     var used = _used;
     if (used == 0) {
       // Set r to 0.
       r._neg = false;
       r._used = 0;
       return;
     }
     var a_used = a._used;
     if (a_used == 0) {
       _copyTo(r);
       return;
     }
     r._ensureLength(used);
-    var digits = _digits;
-    var a_digits = a._digits;
-    var r_digits = r._digits;
-    var c = 0;
-    for (var i = 0; i < a_used; i++) {
-      c += digits[i] - a_digits[i];
-      r_digits[i] = c & DIGIT_MASK;
-      c >>= DIGIT_BITS;
-    }
-    for (var i = a_used; i < used; i++) {
-      c += digits[i];
-      r_digits[i] = c & DIGIT_MASK;
-      c >>= DIGIT_BITS;
-    }
+    _sub(_digits, used, a._digits, a_used, r._digits);
     r._used = used;
     r._clamp();
   }
 
   // r = abs(this) & abs(a).
   void _absAndTo(_Bigint a, _Bigint r) {
     var r_used = (_used < a._used) ? _used : a._used;
     r._ensureLength(r_used);
     var digits = _digits;
     var a_digits = a._digits;
@@ skipping 396 matching lines @@
       _sqrAdd(digits, i, r_digits, used);
     }
     if (r_used > 0) {
       _mulAdd(digits, i, digits, i, r_digits, 2*i, 1);
     }
     r._used = r_used;
     r._neg = false;
     r._clamp();
   }
 
+  // Indices of the arguments of _estQuotientDigit.
+  static const int _YT = 0;  // Index of top digit of divisor y in args array.
+  static const int _QD = 1;  // Index of estimated quotient digit in args array.
+
+  // Estimate args[_QD] = digits[i]:digits[i-1] ~/ args[_YT].
+  static void _estQuotientDigit(Uint32List args, Uint32List digits, int i) {
+    if (digits[i] == args[_YT]) {
+      args[_QD] = DIGIT_MASK;
+    } else {
+      // Chop off one bit, since a Mint cannot hold 2 DIGITs.
+      var qd = ((digits[i] << (DIGIT_BITS - 1)) | (digits[i - 1] >> 1))
+          ~/ (args[_YT] >> 1);
+      if (qd > DIGIT_MASK) {
+        args[_QD] = DIGIT_MASK;
+      } else {
+        args[_QD] = qd;
+      }
+    }
+  }
+
+
   // Truncating division and remainder.
   // If q != null, q = trunc(this / a).
   // If r != null, r = this - a * trunc(this / a).
   void _divRemTo(_Bigint a, _Bigint q, _Bigint r) {
     if (a._used == 0) return;
     if (_used < a._used) {
       if (q != null) {
         // Set q to 0.
         q._neg = false;
         q._used = 0;
@@ skipping 14 matching lines @@
     }
     else {
       a._copyTo(y);
       _copyTo(r);
     }
     // We consider this and a positive. Ignore the copied sign.
     y._neg = false;
     r._neg = false;
     var y_used = y._used;
     var y_digits = y._digits;
-    var y0 = y_digits[y_used - 1];
-    if (y0 == 0) return;
-    var yt = y0 >> 1;  // Chop off one bit, see below. y is normalized: yt != 0.
+    var yt = y_digits[y_used - 1];
+    if (yt == 0) return;
     var i = r._used;
     var j = i - y_used;
     _Bigint t = (q == null) ? new _Bigint() : q;
     y._dlShiftTo(j, t);
     var r_digits = r._digits;
     if (r._compareTo(t) >= 0) {
       r_digits[r._used++] = 1;
       r._subTo(t, r);
     }
     ONE._dlShiftTo(y_used, t);
     t._subTo(y, y);  // Negate y so we can replace sub with _mulAdd later.
     while (y._used < y_used) {
       y_digits[y._used++] = 0;
     }
-    var qd_digit = new Uint32List(1);
+    Uint32List args = new Uint32List(2);

[Cutch, 2014/09/30 17:52:54]

Can you cache this Uint32List somewhere so we don't allocate one each time?

[regis, 2014/09/30 18:35:58]

There is no optimal place to cache this array, and I removed the ugly global one
I used to have. The cost of an allocation is negligible, compared to the cost of
this whole function.

+    args[_YT] = yt;
     while (--j >= 0) {
-      // Estimate quotient digit.
-      // TODO(regis): Move the expensive mint division below to a function that
-      // can be intrinsified using an uint64_t by uint32_t division instruction,
-      // e.g. qd = _estqd(r_digits, --i, y0).
-      if (r_digits[--i] == y0) {
-        qd_digit[0] = DIGIT_MASK;
-      } else {
-        // Chop off one bit, since a Mint cannot hold 2 DIGITs.
-        var qd =
-            ((r_digits[i] << (DIGIT_BITS - 1)) | (r_digits[i - 1] >> 1)) ~/ yt;
-        if (qd > DIGIT_MASK) {
-          qd_digit[0] = DIGIT_MASK;
-        } else {
-          qd_digit[0] = qd;
-        }
-      }
-      _mulAdd(qd_digit, 0, y_digits, 0, r_digits, j, y_used);
-      if (r_digits[i] < qd_digit[0]) {
+      _estQuotientDigit(args, r_digits, --i);
+      _mulAdd(args, _QD, y_digits, 0, r_digits, j, y_used);
+      if (r_digits[i] < args[_QD]) {
         y._dlShiftTo(j, t);
         r._subTo(t, r);
-        while (r_digits[i] < --qd_digit[0]) {
+        while (r_digits[i] < --args[_QD]) {
           r._subTo(t, r);
         }
       }
     }
     if (q != null) {
       r._drShiftTo(y_used, q);
       if (_neg != a._neg) {
         ZERO._subTo(q, q);
       }
     }
@@ skipping 328 matching lines @@
 class _Reduction {
   _Bigint _convert(_Bigint x);
   _Bigint _revert(_Bigint x);
   void _mulTo(_Bigint x, _Bigint y, _Bigint r);
   void _sqrTo(_Bigint x, _Bigint r);
 }
 
 // Montgomery reduction on _Bigint.
 class _Montgomery implements _Reduction {
   _Bigint _m;
-  int _mp;
-  int _mpl;
-  int _mph;
-  int _um;
   int _mused2;
-  Uint32List _u0digit;
+  Uint32List _rho_mu;
+  static const int _RHO = 0;  // Index of rho in _rho_mu array.
+  static const int _MU = 1;  // Index of mu in _rho_mu array.
 
   _Montgomery(m) {
     _m = m._toBigint();
-    _mp = _m._invDigit();
-    _mpl = _mp & _Bigint.DIGIT2_MASK;
-    _mph = _mp >> _Bigint.DIGIT2_BITS;
-    _um = (1 << (_Bigint.DIGIT_BITS - _Bigint.DIGIT2_BITS)) - 1;
     _mused2 = 2*_m._used;
-    _u0digit = new Uint32List(1);
+    _rho_mu = new Uint32List(2);
+    _rho_mu[_RHO] = _m._invDigit();
+  }
+
+  // args[_MU] = args[_RHO]*digits[i] mod DIGIT_BASE.
+  static void _mulMod(Uint32List args, Uint32List digits, int i) {
+    const int MU_MASK = (1 << (_Bigint.DIGIT_BITS - _Bigint.DIGIT2_BITS)) - 1;
+    var rhol = args[_RHO] & _Bigint.DIGIT2_MASK;
+    var rhoh = args[_RHO] >> _Bigint.DIGIT2_BITS;
+    var dh = digits[i] >> _Bigint.DIGIT2_BITS;
+    var dl = digits[i] & _Bigint.DIGIT2_MASK;
+    args[_MU] =
+        (dl*rhol + (((dl*rhoh + dh*rhol) & MU_MASK) << _Bigint.DIGIT2_BITS))
+        & _Bigint.DIGIT_MASK;
   }
 
   // Return x*R mod _m
   _Bigint _convert(_Bigint x) {
     var r = new _Bigint();
     x.abs()._dlShiftTo(_m._used, r);
     r._divRemTo(_m, null, r);
     if (x._neg && !r._neg && r._used > 0) {
       _m._subTo(r, r);
     }
@@ skipping 11 matching lines @@
   // x = x/R mod _m
   void _reduce(_Bigint x) {
     x._ensureLength(_mused2 + 1);
     var x_digits = x._digits;
     while (x._used <= _mused2) {  // Pad x so _mulAdd has enough room later.
       x_digits[x._used++] = 0;
     }
     var m_used = _m._used;
     var m_digits = _m._digits;
     for (var i = 0; i < m_used; ++i) {
-      // Faster way of calculating u0 = x[i]*mp mod DIGIT_BASE.
-      var j = x_digits[i] & _Bigint.DIGIT2_MASK;
-      _u0digit[0] = (j*_mpl + (((j*_mph + (x_digits[i] >> _Bigint.DIGIT2_BITS)
-          *_mpl) & _um) << _Bigint.DIGIT2_BITS)) & _Bigint.DIGIT_MASK;
-      // Use _mulAdd to combine the multiply-shift-add into one call.
-      _Bigint._mulAdd(_u0digit, 0, m_digits, 0, x_digits, i, m_used);
+      _mulMod(_rho_mu, x_digits, i);
+      _Bigint._mulAdd(_rho_mu, _MU, m_digits, 0, x_digits, i, m_used);
     }
     x._clamp();
     x._drShiftTo(m_used, x);
     if (x._compareTo(_m) >= 0) {
       x._subTo(_m, x);
     }
   }
 
   // r = x^2/R mod _m ; x != r
   void _sqrTo(_Bigint x, _Bigint r) {
@@ skipping 39 matching lines @@
   void _sqrTo(_Bigint x, _Bigint r) {
     x._sqrTo(r);
     _reduce(r);
   }
 
   void _mulTo(_Bigint x, _Bigint y, _Bigint r) {
     x._mulTo(y, r);
     _reduce(r);
   }
 }