Process two 32-bit digits as one 64-bit digit in all bigint intrinsics on x64.

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 1020 matching lines @@
       _lShiftTo(nsh, r);
     }
     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;
+    assert((y_used & 1) == 0);
     var y_digits = y._digits;
     Uint32List args = new Uint32List(4);
     args[_YT_LO] = y_digits[y_used - 2];
     args[_YT] = y_digits[y_used - 1];
+    var r_digits = r._digits;
     var i = r._used;
+    if ((i & 1) == 1) {
+      // For 64-bit processing, make sure r has an even number of digits.
+      r_digits[i++] = 0;
+    }
     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_digits[r._used] = 0;  // Set leading zero for 64-bit processing.
       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;
     }
     y_digits[y._used] = 0;  // Set leading zero for 64-bit processing.
-    while (--j >= 0) {
+    while (j > 0) {
       var d0 = _estQuotientDigit(args, r_digits, --i);
+      j -= d0;
       var d1 = _mulAdd(args, _QD, y_digits, 0, r_digits, j, y_used);
       // _estQuotientDigit and _mulAdd must agree on the number of digits to
       // process.
       assert(d0 == d1);
       if (d0 == 1) {
         if (r_digits[i] < args[_QD]) {
           y._dlShiftTo(j, t);
           r._subTo(t, r);
           while (r_digits[i] < --args[_QD]) {
             r._subTo(t, r);
           }
         }
       } else {
-        // TODO(regis): Is this necessary, since intrinsified estimation is more
-        // accurate?
-        if ((r_digits[i] < args[_QD_HI]) ||
-            ((r_digits[i] == args[_QD_HI]) && (r_digits[i-1] < args[_QD]))) {
+        assert(d0 == 2);
+        assert(r_digits[i] <= args[_QD_HI]);
+        if ((r_digits[i] < args[_QD_HI]) || (r_digits[i-1] < args[_QD])) {
           y._dlShiftTo(j, t);
           r._subTo(t, r);
-          assert(args[_QD] > 0);
-          while (r_digits[i] < --args[_QD]) {
+          if (args[_QD] == 0) {
+            --args[_QD_HI];
+          }
+          --args[_QD];
+          assert(r_digits[i] <= args[_QD_HI]);
+          while ((r_digits[i] < args[_QD_HI]) || (r_digits[i-1] < args[_QD])) {
             r._subTo(t, r);
+            if (args[_QD] == 0) {
+              --args[_QD_HI];
+            }
+            --args[_QD];
+            assert(r_digits[i] <= args[_QD_HI]);
           }
         }
+        --i;
       }
     }
     if (q != null) {
       r._drShiftTo(y_used, q);
       if (_neg != a._neg && q._used > 0) {
         q._neg = !q._neg;
       }
     }
     r._used = y_used;
     r._clamp();
@@ skipping 351 matching lines @@
   static const int _X_HI = 1;  // Index of high 32-bits of x (64-bit only).
   static const int _RHO = 2;  // Index of rho.
   static const int _RHO_HI = 3;  // Index of high 32-bits of rho (64-bit only).
   static const int _MU = 4;  // Index of mu.
   static const int _MU_HI = 5;  // Index of high 32-bits of mu (64-bit only).
 
   _Montgomery(m) {
     _m = m._toBigint();
     _mused2 = 2*_m._used;
     _args = new Uint32List(6);
+    // Determine if we can process digit pairs by calling an intrinsic.
+    _digits_per_step = _mulMod(_args, _args, 0);
     _args[_X] = _m._digits[0];
-    _args[_X_HI] = _m._digits[1];
-    _digits_per_step = _invDigit(_args);
+    if (_digits_per_step == 1) {
+      _invDigit(_args);
+    } else {
+      assert(_digits_per_step == 2);
+      _args[_X_HI] = _m._digits[1];
+      _invDigitPair(_args);
+    }
   }
 
-  // Return -1/this % DIGIT_BASE, useful for Montgomery reduction.
-  //
+  // Calculates -1/x % DIGIT_BASE, x is 32-bit digit.
   //         xy == 1 (mod m)
   //         xy =  1+km
   //   xy(2-xy) = (1+km)(1-km)
   // x(y(2-xy)) = 1-k^2 m^2
   // x(y(2-xy)) == 1 (mod m^2)
   // if y is 1/x mod m, then y(2-xy) is 1/x mod m^2
   // Should reduce x and y(2-xy) by m^2 at each step to keep size bounded.
   //
-  // TODO(regis): Intrinsify this method and provide a 64-bit inverse digit pair
-  // on 64-bit platforms.
-  //
   // Operation:
   //   args[_RHO] = 1/args[_X] mod DIGIT_BASE.
-  //   return 1.
-  // Note: intrinsics on 64-bit platform process a digit pair:
-  //   args[_RHO.._RHO_HI] = 1/args[_X.._X_HI] mod DIGIT_BASE^2.
-  //   return 2.
-  static int _invDigit(Uint32List args) {
+  static void _invDigit(Uint32List args) {
     var x = args[_X];
     var y = x & 3;    // y == 1/x mod 2^2
     y = (y*(2 - (x & 0xf)*y)) & 0xf;  // y == 1/x mod 2^4
     y = (y*(2 - (x & 0xff)*y)) & 0xff;  // y == 1/x mod 2^8
     y = (y*(2 - (((x & 0xffff)*y) & 0xffff))) & 0xffff; // y == 1/x mod 2^16
     // Last step - calculate inverse mod DIGIT_BASE directly;
     // Assumes 16 < DIGIT_BITS <= 32 and assumes ability to handle 48-bit ints.
     y = (y*(2 - x*y % _Bigint.DIGIT_BASE)) % _Bigint.DIGIT_BASE;
     // y == 1/x mod DIGIT_BASE
-    // We really want the negative inverse, and - DIGIT_BASE < y < DIGIT_BASE.
-    args[_RHO] = (y > 0) ? _Bigint.DIGIT_BASE - y : -y;
-    return 1;
+    y = -y;  // We really want the negative inverse.
+    args[_RHO] = y & _Bigint.DIGIT_MASK;
   }
 
+
+  // Calculates -1/x % DIGIT_BASE^2, x is a pair of 32-bit digits.
+  // Operation:
+  //   args[_RHO.._RHO_HI] = 1/args[_X.._X_HI] mod DIGIT_BASE^2.
+  static void _invDigitPair(Uint32List args) {
+    var xl = args[_X];  // Lower 32-bit digit of x.
+    var y = xl & 3;    // y == 1/x mod 2^2
+    y = (y*(2 - (xl & 0xf)*y)) & 0xf;  // y == 1/x mod 2^4
+    y = (y*(2 - (xl & 0xff)*y)) & 0xff;  // y == 1/x mod 2^8
+    y = (y*(2 - (((xl & 0xffff)*y) & 0xffff))) & 0xffff;  // y == 1/x mod 2^16
+    y = (y*(2 - ((xl*y) & 0xffffffff))) & 0xffffffff;  // y == 1/x mod 2^32
+    var x = (args[_X_HI] << _Bigint.DIGIT_BITS) | xl;
+    y = (y*(2 - ((x*y) & 0xffffffffffffffff))) & 0xffffffffffffffff;
+    // y == 1/x mod DIGIT_BASE^2
+    y = -y;  // We really want the negative inverse.
+    args[_RHO] = y & _Bigint.DIGIT_MASK;
+    args[_RHO_HI] = (y >> _Bigint.DIGIT_BITS) & _Bigint.DIGIT_MASK;
+  }
+
+
   // Operation:
   //   args[_MU] = args[_RHO]*digits[i] mod DIGIT_BASE.
   //   return 1.
   // Note: intrinsics on 64-bit platform may process a digit pair:
   //   args[_MU.._MU_HI] = args[_RHO.._RHO_HI]*digits[i..i+1] mod DIGIT_BASE^2.
   //   return 2.
   static int _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;
@@ skipping 93 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);
   }
 }