[pypy-commit] pypy default: Re-enable this part of the code and fix it.

[arigo]

Author: Armin Rigo <arigo at tunes.org>
Branch: 
Changeset: r54602:ebca5c7c17dd
Date: 2012-04-21 08:04 +0200
http://bitbucket.org/pypy/pypy/changeset/ebca5c7c17dd/

Log:	Re-enable this part of the code and fix it.

diff --git a/pypy/rlib/rbigint.py b/pypy/rlib/rbigint.py
--- a/pypy/rlib/rbigint.py
+++ b/pypy/rlib/rbigint.py
@@ -40,7 +40,7 @@
 # In that case, do 5 bits at a time.  The potential drawback is that
 # a table of 2**5 intermediate results is computed.
 
-## FIVEARY_CUTOFF = 8   disabled for now
+FIVEARY_CUTOFF = 8
 
 
 def _mask_digit(x):
@@ -456,7 +456,7 @@
 
         # python adaptation: moved macros REDUCE(X) and MULT(X, Y, result)
         # into helper function result = _help_mult(x, y, c)
-        if 1:   ## b.numdigits() <= FIVEARY_CUTOFF:
+        if b.numdigits() <= FIVEARY_CUTOFF:
             # Left-to-right binary exponentiation (HAC Algorithm 14.79)
             # http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf
             i = b.numdigits() - 1
@@ -469,30 +469,51 @@
                         z = _help_mult(z, a, c)
                     j >>= 1
                 i -= 1
-##        else:
-##            This code is disabled for now, because it assumes that
-##            SHIFT is a multiple of 5.  It could be fixed but it looks
-##            like it's more troubles than benefits...
-##
-##            # Left-to-right 5-ary exponentiation (HAC Algorithm 14.82)
-##            # This is only useful in the case where c != None.
-##            # z still holds 1L
-##            table = [z] * 32
-##            table[0] = z
-##            for i in range(1, 32):
-##                table[i] = _help_mult(table[i-1], a, c)
-##            i = b.numdigits() - 1
-##            while i >= 0:
-##                bi = b.digit(i)
-##                j = SHIFT - 5
-##                while j >= 0:
-##                    index = (bi >> j) & 0x1f
-##                    for k in range(5):
-##                        z = _help_mult(z, z, c)
-##                    if index:
-##                        z = _help_mult(z, table[index], c)
-##                    j -= 5
-##                i -= 1
+        else:
+            # Left-to-right 5-ary exponentiation (HAC Algorithm 14.82)
+            # This is only useful in the case where c != None.
+            # z still holds 1L
+            table = [z] * 32
+            table[0] = z
+            for i in range(1, 32):
+                table[i] = _help_mult(table[i-1], a, c)
+            i = b.numdigits()
+            # Note that here SHIFT is not a multiple of 5.  The difficulty
+            # is to extract 5 bits at a time from 'b', starting from the
+            # most significant digits, so that at the end of the algorithm
+            # it falls exactly to zero.
+            # m  = max number of bits = i * SHIFT
+            # m+ = m rounded up to the next multiple of 5
+            # j  = (m+) % SHIFT = (m+) - (i * SHIFT)
+            # (computed without doing "i * SHIFT", which might overflow)
+            j = i % 5
+            if j != 0:
+                j = 5 - j
+            if not we_are_translated():
+                assert j == (i*SHIFT+4)//5*5 - i*SHIFT
+            #
+            accum = r_uint(0)
+            while True:
+                j -= 5
+                if j >= 0:
+                    index = (accum >> j) & 0x1f
+                else:
+                    # 'accum' does not have enough digit.
+                    # must get the next digit from 'b' in order to complete
+                    i -= 1
+                    if i < 0:
+                        break    # done
+                    bi = b.udigit(i)
+                    index = ((accum << (-j)) | (bi >> (j+SHIFT))) & 0x1f
+                    accum = bi
+                    j += SHIFT
+                #
+                for k in range(5):
+                    z = _help_mult(z, z, c)
+                if index:
+                    z = _help_mult(z, table[index], c)
+            #
+            assert j == -5
 
         if negativeOutput and z.sign != 0:
             z = z.sub(c)
diff --git a/pypy/rlib/test/test_rbigint.py b/pypy/rlib/test/test_rbigint.py
--- a/pypy/rlib/test/test_rbigint.py
+++ b/pypy/rlib/test/test_rbigint.py
@@ -379,6 +379,18 @@
         for n, expected in [(37, 9), (1291, 931), (67889, 39464)]:
             v = two.pow(t, rbigint.fromint(n))
             assert v.toint() == expected
+        #
+        # more tests, comparing against CPython's answer
+        enabled = sample(range(5*32), 10)
+        for i in range(5*32):
+            t = t.mul(two)      # add one random bit
+            if random() >= 0.5:
+                t = t.add(rbigint.fromint(1))
+            if i not in enabled:
+                continue    # don't take forever
+            n = randint(1, sys.maxint)
+            v = two.pow(t, rbigint.fromint(n))
+            assert v.toint() == pow(2, t.tolong(), n)
 
     def test_pow_lln(self):
         x = 10L