Index: Doc/lib/libfuncs.tex =================================================================== RCS file: /cvsroot/python/python/dist/src/Doc/lib/libfuncs.tex,v retrieving revision 1.166 diff -c -r1.166 libfuncs.tex *** Doc/lib/libfuncs.tex 2 Jul 2004 06:41:04 -0000 1.166 --- Doc/lib/libfuncs.tex 19 Jul 2004 09:21:08 -0000 *************** *** 691,708 **** arguments must have numeric types. With mixed operand types, the coercion rules for binary arithmetic operators apply. For int and long int operands, the result has the same type as the operands ! (after coercion) unless the second argument is negative; in that ! case, all arguments are converted to float and a float result is ! delivered. For example, \code{10**2} returns \code{100}, but ! \code{10**-2} returns \code{0.01}. (This last feature was added in ! Python 2.2. In Python 2.1 and before, if both arguments were of integer ! types and the second argument was negative, an exception was raised.) ! If the second argument is negative, the third argument must be omitted. ! If \var{z} is present, \var{x} and \var{y} must be of integer types, ! and \var{y} must be non-negative. (This restriction was added in ! Python 2.2. In Python 2.1 and before, floating 3-argument \code{pow()} ! returned platform-dependent results depending on floating-point ! rounding accidents.) \end{funcdesc} \begin{funcdesc}{property}{\optional{fget\optional{, fset\optional{, --- 691,712 ---- arguments must have numeric types. With mixed operand types, the coercion rules for binary arithmetic operators apply. For int and long int operands, the result has the same type as the operands ! (after coercion) unless the second argument is negative and the ! third is absent; in that case, all arguments are converted to float ! and a float result is delivered. For example, \code{10**2} returns ! \code{100}, but \code{10**-2} returns \code{0.01}. (This last ! feature was added in Python 2.2. In Python 2.1 and before, if ! the first two arguments were of integer types and the second argument ! was negative, an exception was raised.) ! ! If the second argument is a negative integer type, and the third ! argument is present, then the first and third arguments must be ! integer types, and must be relatively prime (or an exception will ! be raised). The resulting calculation will be performed in the group ! of integers modulo \var{z}. In other words, the multiplicative inverse ! of \var{x} modulo \var{z} will be calculated, and that value will be ! raised to the absolute value of \var{y}, modulo \var{z}. This behavior was ! added in Python 2.4. \end{funcdesc} \begin{funcdesc}{property}{\optional{fget\optional{, fset\optional{, Index: Lib/test/test_long.py =================================================================== RCS file: /cvsroot/python/python/dist/src/Lib/test/test_long.py,v retrieving revision 1.23 diff -c -r1.23 test_long.py *** Lib/test/test_long.py 3 Feb 2003 15:25:01 -0000 1.23 --- Lib/test/test_long.py 19 Jul 2004 09:21:12 -0000 *************** *** 6,12 **** SHIFT = 15 BASE = 2 ** SHIFT MASK = BASE - 1 ! KARATSUBA_CUTOFF = 35 # from longobject.c # Max number of base BASE digits to use in test cases. Doubling # this will more than double the runtime. --- 6,12 ---- SHIFT = 15 BASE = 2 ** SHIFT MASK = BASE - 1 ! KARATSUBA_CUTOFF = 70 # from longobject.c # Max number of base BASE digits to use in test cases. Doubling # this will more than double the runtime. *************** *** 347,364 **** for z in special: if z != 0 : ! if y >= 0: expected = pow(longx, longy, long(z)) got = pow(x, y, z) checkit('pow', x, y, '%', z) ! else: ! try: ! pow(longx, longy, long(z)) ! except TypeError: ! pass ! else: ! raise TestFailed("pow%r should have raised " ! "TypeError" % ((longx, longy, long(z)),)) # ---------------------------------------- tests of long->float overflow --- 347,369 ---- for z in special: if z != 0 : ! valueError = False ! try: expected = pow(longx, longy, long(z)) + except ValueError: + if longy < 0: + valueError = True #longy not coprime with z + else: + raise TestFailed("pow%r should not have raised " + "ValueError" % ((longx, longy, long(z)),)) + + try: got = pow(x, y, z) checkit('pow', x, y, '%', z) ! except ValueError: ! if not valueError: ! raise TestFailed("pow%r should not have raised " ! "ValueError" % ((x, y, z),)) # ---------------------------------------- tests of long->float overflow Index: Lib/test/test_pow.py =================================================================== RCS file: /cvsroot/python/python/dist/src/Lib/test/test_pow.py,v retrieving revision 1.19 diff -c -r1.19 test_pow.py *** Lib/test/test_pow.py 24 May 2003 20:18:24 -0000 1.19 --- Lib/test/test_pow.py 19 Jul 2004 09:21:12 -0000 *************** *** 95,100 **** --- 95,139 ---- pow(long(i),j,k) ) + def test_powinv(self): + basevals = [144, 144 + (6817*7), 144 - (6817*7)] + expvals = [(0,1), (1,144), (6,5410), (-1,2225), (-6,6473)] + modvals = [6817, 0, -6817] + + for mod in modvals: + for exp, result in expvals: + for base in basevals: + if mod != 0: + def f(v1, v2, v3, t): + x = pow(v1, v2, v3) + if mod > 0: + self.assertEquals(x, result) + if mod < 0: + self.assertEquals(x, result + mod) + self.assertEquals(type(x), t) + else: + def f(v1, v2, v3, t): + self.assertRaises(ValueError, pow, v1, v2, v3) + f(base, exp, mod, int) + f(long(base), exp, mod, long) + f(base, long(exp), mod, long) + f(base, exp, long(mod), long) + f(long(base), long(exp), mod, long) + f(long(base), exp, long(mod), long) + f(base, long(exp), long(mod), long) + f(long(base), long(exp), long(mod), long) + + import random + random.seed("abc") + for x in range(10): + b = random.getrandbits(256) + m = random.getrandbits(256) + try: + r = pow(b, -1, m) + self.assertEquals((r * b) % m, 1L) + except ValueError: + pass + def test_bug643260(self): class TestRpow: def __rpow__(self, other): Index: Objects/intobject.c =================================================================== RCS file: /cvsroot/python/python/dist/src/Objects/intobject.c,v retrieving revision 2.111 diff -c -r2.111 intobject.c *** Objects/intobject.c 26 Jun 2004 23:22:57 -0000 2.111 --- Objects/intobject.c 19 Jul 2004 09:21:15 -0000 *************** *** 611,634 **** } } static PyObject * int_pow(PyIntObject *v, PyIntObject *w, PyIntObject *z) { register long iv, iw, iz=0, ix, temp, prev; ! CONVERT_TO_LONG(v, iv); CONVERT_TO_LONG(w, iw); - if (iw < 0) { - if ((PyObject *)z != Py_None) { - PyErr_SetString(PyExc_TypeError, "pow() 2nd argument " - "cannot be negative when 3rd argument specified"); - return NULL; - } - /* Return a float. This works because we know that - this calls float_pow() which converts its - arguments to double. */ - return PyFloat_Type.tp_as_number->nb_power( - (PyObject *)v, (PyObject *)w, (PyObject *)z); - } if ((PyObject *)z != Py_None) { CONVERT_TO_LONG(z, iz); if (iz == 0) { --- 611,633 ---- } } + /* Helper function for int_pow */ + static PyObject* long_to_int(PyObject *longObject) { + PyIntObject *intObject = NULL; + if (longObject == NULL) + return NULL; + intObject = (PyIntObject *)PyLong_Type.tp_as_number-> + nb_int(longObject); + Py_DECREF(longObject); + return (PyObject *)intObject; + } + static PyObject * int_pow(PyIntObject *v, PyIntObject *w, PyIntObject *z) { register long iv, iw, iz=0, ix, temp, prev; ! CONVERT_TO_LONG(v, iv); CONVERT_TO_LONG(w, iw); if ((PyObject *)z != Py_None) { CONVERT_TO_LONG(z, iz); if (iz == 0) { *************** *** 637,642 **** --- 636,652 ---- return NULL; } } + if (iw < 0) { + if ((PyObject *)z != Py_None) + return long_to_int(PyLong_Type.tp_as_number->nb_power( + (PyObject *)v, (PyObject *)w, (PyObject *)z)); + + /* Return a float. This works because we know that + this calls float_pow() which converts its + arguments to double. */ + return PyFloat_Type.tp_as_number->nb_power( + (PyObject *)v, (PyObject *)w, (PyObject *)z); + } /* * XXX: The original exponentiation code stopped looping * when temp hit zero; this code will continue onwards *************** *** 656,665 **** if (ix / temp != prev) { if (err_ovf("integer exponentiation")) return NULL; ! return PyLong_Type.tp_as_number->nb_power( (PyObject *)v, (PyObject *)w, ! (PyObject *)z); } } iw >>= 1; /* Shift exponent down by 1 bit */ --- 666,675 ---- if (ix / temp != prev) { if (err_ovf("integer exponentiation")) return NULL; ! return long_to_int(PyLong_Type.tp_as_number->nb_power( (PyObject *)v, (PyObject *)w, ! (PyObject *)z)); } } iw >>= 1; /* Shift exponent down by 1 bit */ *************** *** 669,676 **** if (prev!=0 && temp/prev!=prev) { if (err_ovf("integer exponentiation")) return NULL; ! return PyLong_Type.tp_as_number->nb_power( ! (PyObject *)v, (PyObject *)w, (PyObject *)z); } if (iz) { /* If we did a multiplication, perform a modulo */ --- 679,686 ---- if (prev!=0 && temp/prev!=prev) { if (err_ovf("integer exponentiation")) return NULL; ! return long_to_int(PyLong_Type.tp_as_number->nb_power( ! (PyObject *)v, (PyObject *)w, (PyObject *)z)); } if (iz) { /* If we did a multiplication, perform a modulo */ *************** *** 685,692 **** ix = mod; break; case DIVMOD_OVERFLOW: ! return PyLong_Type.tp_as_number->nb_power( ! (PyObject *)v, (PyObject *)w, (PyObject *)z); default: return NULL; } --- 695,703 ---- ix = mod; break; case DIVMOD_OVERFLOW: ! return long_to_int(PyLong_Type.tp_as_number->nb_power( ! (PyObject *)v, (PyObject *)w, (PyObject *)z)); ! default: return NULL; } Index: Objects/longobject.c =================================================================== RCS file: /cvsroot/python/python/dist/src/Objects/longobject.c,v retrieving revision 1.161 diff -c -r1.161 longobject.c *** Objects/longobject.c 28 Jun 2003 20:04:25 -0000 1.161 --- Objects/longobject.c 19 Jul 2004 09:21:18 -0000 *************** *** 11,18 **** /* For long multiplication, use the O(N**2) school algorithm unless * both operands contain more than KARATSUBA_CUTOFF digits (this * being an internal Python long digit, in base BASE). */ ! #define KARATSUBA_CUTOFF 35 #define ABS(x) ((x) < 0 ? -(x) : (x)) --- 11,23 ---- /* For long multiplication, use the O(N**2) school algorithm unless * both operands contain more than KARATSUBA_CUTOFF digits (this * being an internal Python long digit, in base BASE). + * + * For modular exponentiation, use the binary left-to-right algorithm + * unless the exponent contains more than FIVEARY_CUTOFF digits. */ ! #define KARATSUBA_CUTOFF 70 ! #define KARATSUBA_SQUARE_CUTOFF (2 * KARATSUBA_CUTOFF) ! #define FIVEARY_CUTOFF 8 #define ABS(x) ((x) < 0 ? -(x) : (x)) *************** *** 1717,1744 **** return NULL; memset(z->ob_digit, 0, z->ob_size * sizeof(digit)); ! for (i = 0; i < size_a; ++i) { ! twodigits carry = 0; ! twodigits f = a->ob_digit[i]; ! int j; ! digit *pz = z->ob_digit + i; ! ! SIGCHECK({ ! Py_DECREF(z); ! return NULL; ! }) ! for (j = 0; j < size_b; ++j) { ! carry += *pz + b->ob_digit[j] * f; ! *pz++ = (digit) (carry & MASK); ! carry >>= SHIFT; ! } ! for (; carry != 0; ++j) { ! assert(i+j < z->ob_size); ! carry += *pz; ! *pz++ = (digit) (carry & MASK); ! carry >>= SHIFT; ! } ! } return long_normalize(z); } --- 1722,1786 ---- return NULL; memset(z->ob_digit, 0, z->ob_size * sizeof(digit)); ! if (a == b) ! /* ! Efficient squaring per HAC, Algorithm 14.16: ! http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf ! Gives slightly less than a 2x speedup when a == b ! */ ! for (i = 0; i < size_a; ++i) { ! twodigits carry = 0; ! twodigits f = a->ob_digit[i]; ! digit *pz = z->ob_digit + (i << 1); ! digit *pb = a->ob_digit + i + 1; ! digit *pbend = a->ob_digit + size_a; ! ! SIGCHECK({ ! Py_DECREF(z); ! return NULL; ! }) ! ! carry = *pz + f * f; ! *pz++ = (digit)(carry & MASK); ! carry >>= SHIFT; ! ! f <<= 1; /*f is now double-precision */ ! while (pb < pbend) { ! carry += *pz + *pb++ * f; ! *pz++ = (digit)(carry & MASK); ! carry >>= SHIFT; ! } ! if (carry) { ! carry += *pz; ! *pz++ = (digit)(carry & MASK); ! carry >>= SHIFT; ! } ! if (carry) ! *pz += (digit)(carry & MASK); ! assert((carry >> SHIFT) == 0); ! } ! else ! for (i = 0; i < size_a; ++i) { ! twodigits carry = 0; ! twodigits f = a->ob_digit[i]; ! digit *pz = z->ob_digit + i; ! digit *pb = b->ob_digit; ! digit *pbend = b->ob_digit + size_b; ! ! SIGCHECK({ ! Py_DECREF(z); ! return NULL; ! }) ! ! while (pb < pbend) { ! carry += *pz + *pb++ * f; ! *pz++ = (digit)(carry & MASK); ! carry >>= SHIFT; ! } ! if (carry) ! *pz += (digit)(carry & MASK); ! assert((carry >> SHIFT) == 0); ! } return long_normalize(z); } *************** *** 1816,1822 **** } /* Use gradeschool math when either number is too small. */ ! if (asize <= KARATSUBA_CUTOFF) { if (asize == 0) return _PyLong_New(0); else --- 1858,1864 ---- } /* Use gradeschool math when either number is too small. */ ! if (asize <= (a == b ? KARATSUBA_SQUARE_CUTOFF : KARATSUBA_CUTOFF)) { if (asize == 0) return _PyLong_New(0); else *************** *** 1837,1843 **** if (kmul_split(a, shift, &ah, &al) < 0) goto fail; assert(ah->ob_size > 0); /* the split isn't degenerate */ ! if (kmul_split(b, shift, &bh, &bl) < 0) goto fail; /* The plan: * 1. Allocate result space (asize + bsize digits: that's always --- 1879,1892 ---- if (kmul_split(a, shift, &ah, &al) < 0) goto fail; assert(ah->ob_size > 0); /* the split isn't degenerate */ ! if (a == b) { ! bh = ah; ! bl = al; ! Py_INCREF(bh); ! Py_INCREF(bl); ! } ! else ! if (kmul_split(b, shift, &bh, &bl) < 0) goto fail; /* The plan: * 1. Allocate result space (asize + bsize digits: that's always *************** *** 1906,1915 **** Py_DECREF(al); ah = al = NULL; ! if ((t2 = x_add(bh, bl)) == NULL) { ! Py_DECREF(t1); ! goto fail; } Py_DECREF(bh); Py_DECREF(bl); bh = bl = NULL; --- 1955,1969 ---- Py_DECREF(al); ah = al = NULL; ! if (a == b) { ! t2 = t1; ! Py_INCREF(t2); } + else + if ((t2 = x_add(bh, bl)) == NULL) { + Py_DECREF(t1); + goto fail; + } Py_DECREF(bh); Py_DECREF(bl); bh = bl = NULL; *************** *** 1987,1993 **** * of slices, each with a->ob_size digits, and multiply the slices by a, * one at a time. This gives k_mul balanced inputs to work with, and is * also cache-friendly (we compute one double-width slice of the result ! * at a time, then move on, never bactracking except for the helpful * single-width slice overlap between successive partial sums). */ static PyLongObject * --- 2041,2047 ---- * of slices, each with a->ob_size digits, and multiply the slices by a, * one at a time. This gives k_mul balanced inputs to work with, and is * also cache-friendly (we compute one double-width slice of the result ! * at a time, then move on, never backtracking except for the helpful * single-width slice overlap between successive partial sums). */ static PyLongObject * *************** *** 2109,2151 **** Py_DECREF(div); div = temp; } ! *pdiv = div; ! *pmod = mod; return 0; } static PyObject * long_div(PyObject *v, PyObject *w) { ! PyLongObject *a, *b, *div, *mod; CONVERT_BINOP(v, w, &a, &b); ! if (l_divmod(a, b, &div, &mod) < 0) { Py_DECREF(a); Py_DECREF(b); return NULL; } Py_DECREF(a); Py_DECREF(b); - Py_DECREF(mod); return (PyObject *)div; } static PyObject * long_classic_div(PyObject *v, PyObject *w) { ! PyLongObject *a, *b, *div, *mod; CONVERT_BINOP(v, w, &a, &b); if (Py_DivisionWarningFlag && PyErr_Warn(PyExc_DeprecationWarning, "classic long division") < 0) div = NULL; ! else if (l_divmod(a, b, &div, &mod) < 0) div = NULL; - else - Py_DECREF(mod); Py_DECREF(a); Py_DECREF(b); --- 2163,2208 ---- Py_DECREF(div); div = temp; } ! if (pdiv != NULL) ! *pdiv = div; ! else ! Py_DECREF(div); ! if (pmod != NULL) ! *pmod = mod; ! else ! Py_DECREF(mod); return 0; } static PyObject * long_div(PyObject *v, PyObject *w) { ! PyLongObject *a, *b, *div; CONVERT_BINOP(v, w, &a, &b); ! if (l_divmod(a, b, &div, NULL) < 0) { Py_DECREF(a); Py_DECREF(b); return NULL; } Py_DECREF(a); Py_DECREF(b); return (PyObject *)div; } static PyObject * long_classic_div(PyObject *v, PyObject *w) { ! PyLongObject *a, *b, *div; CONVERT_BINOP(v, w, &a, &b); if (Py_DivisionWarningFlag && PyErr_Warn(PyExc_DeprecationWarning, "classic long division") < 0) div = NULL; ! else if (l_divmod(a, b, &div, NULL) < 0) div = NULL; Py_DECREF(a); Py_DECREF(b); *************** *** 2197,2214 **** static PyObject * long_mod(PyObject *v, PyObject *w) { ! PyLongObject *a, *b, *div, *mod; CONVERT_BINOP(v, w, &a, &b); ! if (l_divmod(a, b, &div, &mod) < 0) { Py_DECREF(a); Py_DECREF(b); return NULL; } Py_DECREF(a); Py_DECREF(b); - Py_DECREF(div); return (PyObject *)mod; } --- 2254,2270 ---- static PyObject * long_mod(PyObject *v, PyObject *w) { ! PyLongObject *a, *b, *mod; CONVERT_BINOP(v, w, &a, &b); ! if (l_divmod(a, b, NULL, &mod) < 0) { Py_DECREF(a); Py_DECREF(b); return NULL; } Py_DECREF(a); Py_DECREF(b); return (PyObject *)mod; } *************** *** 2239,2357 **** return z; } ! static PyObject * ! long_pow(PyObject *v, PyObject *w, PyObject *x) { ! PyLongObject *a, *b; ! PyObject *c; ! PyLongObject *z, *div, *mod; ! int size_b, i; ! CONVERT_BINOP(v, w, &a, &b); ! if (PyLong_Check(x) || Py_None == x) { ! c = x; ! Py_INCREF(x); ! } ! else if (PyInt_Check(x)) { ! c = PyLong_FromLong(PyInt_AS_LONG(x)); ! } ! else { ! Py_DECREF(a); ! Py_DECREF(b); ! Py_INCREF(Py_NotImplemented); ! return Py_NotImplemented; ! } ! ! if (c != Py_None && ((PyLongObject *)c)->ob_size == 0) { ! PyErr_SetString(PyExc_ValueError, ! "pow() 3rd argument cannot be 0"); ! z = NULL; ! goto error; ! } ! ! size_b = b->ob_size; ! if (size_b < 0) { ! Py_DECREF(a); ! Py_DECREF(b); ! Py_DECREF(c); ! if (x != Py_None) { ! PyErr_SetString(PyExc_TypeError, "pow() 2nd argument " ! "cannot be negative when 3rd argument specified"); ! return NULL; ! } ! /* Return a float. This works because we know that ! this calls float_pow() which converts its ! arguments to double. */ ! return PyFloat_Type.tp_as_number->nb_power(v, w, x); ! } ! z = (PyLongObject *)PyLong_FromLong(1L); ! for (i = 0; i < size_b; ++i) { ! digit bi = b->ob_digit[i]; ! int j; ! for (j = 0; j < SHIFT; ++j) { ! PyLongObject *temp; ! if (bi & 1) { ! temp = (PyLongObject *)long_mul(z, a); ! Py_DECREF(z); ! if (c!=Py_None && temp!=NULL) { ! if (l_divmod(temp,(PyLongObject *)c, ! &div,&mod) < 0) { ! Py_DECREF(temp); ! z = NULL; ! goto error; ! } ! Py_XDECREF(div); ! Py_DECREF(temp); ! temp = mod; ! } ! z = temp; ! if (z == NULL) ! break; ! } ! bi >>= 1; ! if (bi == 0 && i+1 == size_b) ! break; ! temp = (PyLongObject *)long_mul(a, a); ! Py_DECREF(a); ! if (c!=Py_None && temp!=NULL) { ! if (l_divmod(temp, (PyLongObject *)c, &div, ! &mod) < 0) { ! Py_DECREF(temp); ! z = NULL; ! goto error; ! } ! Py_XDECREF(div); ! Py_DECREF(temp); ! temp = mod; ! } ! a = temp; ! if (a == NULL) { ! Py_DECREF(z); ! z = NULL; ! break; ! } ! } ! if (a == NULL || z == NULL) ! break; ! } ! if (c!=Py_None && z!=NULL) { ! if (l_divmod(z, (PyLongObject *)c, &div, &mod) < 0) { ! Py_DECREF(z); ! z = NULL; ! } ! else { ! Py_XDECREF(div); ! Py_DECREF(z); ! z = mod; ! } ! } ! error: ! Py_XDECREF(a); ! Py_DECREF(b); ! Py_DECREF(c); ! return (PyObject *)z; } static PyObject * --- 2295,2734 ---- return z; } ! /* ! Helper for long_pow ! -------------------- ! Use the Extended Euclidean algorithm to calculate the inverse of a modulo b. ! I.e., calculate c such that (a * c) % b == 1. ! ! This is translated from pseudocode in "Practical Cryptography" section ! 11.3.5, Ferguson & Schneier. ! ! Below is the algorithm in Python: ! ! def invMod(a, b): ! c, d = a, b ! uc, ud = 1, 0 ! while c != 0: ! q = d // c ! c, d = d - (q * c), c ! uc, ud = ud - (q * uc), uc ! if d == 1: ! return ud % b ! else: ! raise ValueError("argument has no inverse") ! */ ! PyLongObject* l_invmod(PyLongObject *a, PyLongObject *b) { ! PyLongObject *c = NULL; ! PyLongObject *d = NULL; ! PyLongObject *uc = NULL; ! PyLongObject *ud = NULL; ! PyLongObject *q = NULL; ! PyLongObject *temp = NULL; ! PyLongObject *temp2 = NULL; ! PyLongObject *retval = NULL; ! ! c = (PyLongObject *)_PyLong_Copy(a); ! if (c == NULL) ! goto error; ! d = (PyLongObject *)_PyLong_Copy(b); ! if (d == NULL) ! goto error; ! uc = (PyLongObject *)PyLong_FromLong(1L); ! if (uc == NULL) ! goto error; ! ud = (PyLongObject *)PyLong_FromLong(0L); ! if (ud == NULL) ! goto error; ! ! while (c->ob_size != 0) { ! Py_XDECREF(q); ! if (l_divmod(d, c, &q, NULL) < 0) ! goto error; ! ! Py_XDECREF(temp); ! temp = (PyLongObject *)long_mul(q, c); ! if (temp == NULL) ! goto error; ! temp2 = (PyLongObject *)long_sub(d, temp); ! if (temp2 == NULL) ! goto error; ! Py_DECREF(d); ! d = c; ! c = temp2; ! temp2 = NULL; ! ! Py_XDECREF(temp); ! temp = (PyLongObject *)long_mul(q, uc); ! if (temp == NULL) ! goto error; ! temp2 = (PyLongObject *)long_sub(ud, temp); ! if (temp2 == NULL) ! goto error; ! Py_DECREF(ud); ! ud = uc; ! uc = temp2; ! temp2 = NULL; ! } ! ! if ((d->ob_size == 1) && (d->ob_digit[0] == 1)) { ! if (l_divmod(ud, b, NULL, &retval) < 0) ! goto error; ! } ! else ! PyErr_SetString(PyExc_ValueError, "argument has no inverse"); ! ! error: ! Py_XDECREF(c); ! Py_XDECREF(d); ! Py_XDECREF(uc); ! Py_XDECREF(ud); ! Py_XDECREF(q); ! Py_XDECREF(temp); ! Py_XDECREF(temp2); ! return retval; ! } ! ! /* ! Helper for long_pow ! -------------------- ! Implement Montgomery Reduction per HAC, Algorithm 14.32: ! http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf ! ! INPUT: T, m, m' ! BASE = 2**15 ! n = |m| (# of digits in m) ! R = BASE**n (next power of BASE greater than m) ! OUTPUT: T*(R**-1) % m ! */ ! PyLongObject* mont_reduce(PyLongObject *T, PyLongObject *m, int mPrime) ! { ! int n, i; ! PyLongObject* A = NULL; ! PyLongObject* retval = NULL; ! ! n = m->ob_size; ! ! /* A = T */ ! A = (PyLongObject *)_PyLong_New(2*n+1); ! if (A == NULL) ! return NULL; ! if (T->ob_size > 2*n) ! goto error; ! memset(A->ob_digit, 0, A->ob_size * sizeof(digit)); ! memcpy(A->ob_digit, T->ob_digit, T->ob_size * sizeof(digit)); ! ! for (i=0; i < n; i++) { ! digit *aptr, *mptr, *mendptr; ! twodigits carry = 0; ! int u; ! ! /* u = a[i] * m' mod base */ ! u = (A->ob_digit[i] * mPrime) & MASK; ! ! /* A += u * m * base**i ! Multiply and add in place: ! A[i+j] += (u * m[j]) + carry */ ! aptr = A->ob_digit + i; ! mptr = m->ob_digit; ! mendptr = mptr + n; ! while (mptr < mendptr) { ! carry += (twodigits)( (*aptr) + (u * (*mptr++)) ); ! *aptr++ = (digit) (carry & MASK); ! carry >>= SHIFT; ! } ! if (carry) ! *aptr += (digit)(carry & MASK); ! assert((carry >> SHIFT) == 0); ! } ! ! /* A = A/base**n (i.e. A >>= n) */ ! memmove(A->ob_digit, A->ob_digit + n, (n+1) * sizeof(digit)); ! A->ob_size = n+1; ! long_normalize(A); ! ! /* If A > m then A = A - m */ ! if (long_compare(A, m) < 0) ! return A; ! else ! retval = (PyLongObject *)long_sub(A, m); ! error: ! Py_DECREF(A); ! return retval; ! } ! static PyObject * ! long_pow(PyObject *v, PyObject *w, PyObject *x) ! { ! PyLongObject *a, *b, *c; /* a,b,c = v,w,x */ ! int negativeOutput = 0; /* if x<0 return negative output */ ! PyLongObject *z = NULL; /* accumulated result */ ! int i, j, k; /* counters */ ! PyLongObject *temp = NULL; ! ! /* Montgomery values */ ! int doMontgomery = 0; ! PyLongObject *R = NULL; ! static PyLongObject *base = NULL; ! int cPrime = 0; ! ! /* 5-ary values */ ! int mask, index; ! PyLongObject *table[32] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ ! 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; ! ! /* a, b, c = v, w, x */ ! CONVERT_BINOP(v, w, &a, &b); ! if (PyLong_Check(x)) { ! c = (PyLongObject *)x; ! Py_INCREF(x); ! } ! else if (PyInt_Check(x)) ! c = (PyLongObject *)PyLong_FromLong(PyInt_AS_LONG(x)); ! else if (x == Py_None) ! c = NULL; ! else { ! Py_DECREF(a); ! Py_DECREF(b); ! Py_INCREF(Py_NotImplemented); ! return Py_NotImplemented; ! } ! ! if (c) { ! /* if modulus == 0: ! raise ValueError() */ ! if (c->ob_size == 0) { ! PyErr_SetString(PyExc_ValueError, ! "pow() 3rd argument cannot be 0"); ! goto error; ! } ! ! /* if modulus < 0: ! negativeOutput = True ! modulus = -modulus */ ! if (c->ob_size < 0) { ! negativeOutput = 1; ! temp = (PyLongObject *)_PyLong_Copy(c); ! if (temp == NULL) ! goto error; ! Py_DECREF(c); ! c = temp; ! temp = NULL; ! c->ob_size *= -1; ! } ! ! /* if modulus == 1: ! return 0 ! This check is not necessary when Montgomery Reduction ! is in the code. However, if Montgomery Reduction is ! removed or disabled, we want to be able to set z=1 ! and be sure that z is inside the group modulo c. ! */ ! if ((c->ob_size == 1) && (c->ob_digit[0] == 1)) { ! z = (PyLongObject *)PyLong_FromLong(0L); ! goto error; ! } ! ! /* if base < 0: ! base = base % modulus ! The l_invmod() Euclidean Algorithm requires a positive value. ! Even if we don't go through that code path, having the base ! positive just makes things easier. ! */ ! if (a->ob_size < 0) { ! if (l_divmod(a, c, NULL, &temp) < 0) ! goto error; ! Py_DECREF(a); ! a = temp; ! temp = NULL; ! } ! } ! ! /*if exponent < 0: */ ! if (b->ob_size < 0) { ! /* if modulus: ! base = invMod(base, modulus) ! exponent = -exponent */ ! if (c) { ! temp = l_invmod(a, c); ! if (temp == NULL) { ! goto error; ! } ! Py_DECREF(a); ! a = temp; ! temp = NULL; ! ! temp = (PyLongObject *)_PyLong_Copy(b); ! if (temp == NULL) ! goto error; ! Py_DECREF(b); ! b = temp; ! temp = NULL; ! b->ob_size *= -1; ! } ! else { ! /* else return a float. This works because we know ! that this calls float_pow() which converts its ! arguments to double. */ ! return PyFloat_Type.tp_as_number->nb_power(v, w, x); ! } ! } ! ! /* At this point a, b, and c are guaranteed non-negative UNLESS ! c is NULL, in which case a may be negative. */ ! ! /* If modulus is odd, use Montgomery Reduction */ ! if (c && (c->ob_digit[0] & 1)) ! { ! doMontgomery = 1; ! ! /*First, calculate the Montgomery parameters (R, m') ! (or in our case, (R, c')). ! ! Then set (a, z) to the Montgomery representations ! of (a, 1), in preparation for the exponentiation loop. ! */ ! ! /* R = BASE ** n for smallest n such that R > c, (n == c->ob_size) */ ! R = _PyLong_New(c->ob_size + 1); ! if (R == NULL) ! goto error; ! memset(R->ob_digit, 0, R->ob_size * sizeof(digit)); ! R->ob_digit[R->ob_size-1] = 1; ! ! /* c' = (-pow(c, -1, BASE)) % BASE */ ! if (base == NULL) { ! base = (PyLongObject *)PyLong_FromLong(BASE); ! if (base == NULL) ! goto error; ! } ! temp = (PyLongObject *)l_invmod(c, base); ! if (temp == NULL) ! goto error; ! cPrime = BASE - temp->ob_digit[0]; ! Py_DECREF(temp); ! temp = NULL; ! ! /* a = a * R % c # keep in mind a * R == a << size_c */ ! temp = _PyLong_New(c->ob_size + a->ob_size); ! if (temp == NULL) ! goto error; ! memset(temp->ob_digit, 0, c->ob_size * sizeof(digit)); ! memcpy(temp->ob_digit + c->ob_size, a->ob_digit, ! a->ob_size * sizeof(digit)); ! ! Py_DECREF(a); ! a = NULL; ! if (l_divmod(temp, c, NULL, &a) < 0) ! goto error; ! Py_DECREF(temp); ! temp = NULL; ! ! /* z = 1 * R % c */ ! if (l_divmod(R, c, NULL, &z) < 0) ! goto error; ! } ! else ! z = (PyLongObject *)PyLong_FromLong(1L); ! ! /* Perform a modular reduction, using Montgomery if applicable */ ! #define REDUCE(X)\ ! if (c) {\ ! if (doMontgomery) {\ ! temp = mont_reduce(X, c, cPrime);\ ! if (temp == NULL) {\ ! Py_DECREF(z);\ ! z = NULL;\ ! goto error;\ ! }\ ! }\ ! else {\ ! if (l_divmod(X, c, NULL, &temp) < 0) {\ ! Py_DECREF(z);\ ! z = NULL;\ ! goto error;\ ! }\ ! }\ ! Py_XDECREF(X);\ ! X = temp;\ ! temp = NULL;\ ! } ! ! /* Multiply two values, then reduce the result */ ! #define MULT(X, Y, result) {\ ! temp = (PyLongObject *)long_mul(X, Y);\ ! if (temp == NULL) {\ ! Py_DECREF(z);\ ! z = NULL;\ ! goto error;\ ! }\ ! Py_XDECREF(result);\ ! result = temp;\ ! temp = NULL;\ ! REDUCE(result)\ ! } ! ! if (b->ob_size <= FIVEARY_CUTOFF) ! /* Left-to-right binary exponentiation (HAC Algorithm 14.79) */ ! /* http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf */ ! for (i = b->ob_size-1; i >= 0; --i) { ! digit bi = b->ob_digit[i]; ! ! for (j = 1 << (SHIFT-1); j != 0; j >>= 1) { ! MULT(z, z, z) ! if (bi & j) ! MULT(z, a, z) ! } ! } ! else { ! /* Left-to-right 5-ary exponentiation (HAC Algorithm 14.82) */ ! table[0] = z; ! for (i=1; i < 32; ++i) ! MULT(table[i-1], a, table[i]) ! ! for (i = b->ob_size-1; i >= 0; --i) { ! digit bi = b->ob_digit[i]; ! ! for (j=2; j >= 0; --j) { ! for (k=0; k < 5; ++k) ! MULT(z, z, z) ! mask = 31 << (j * 5); ! index = (bi & mask) >> (j * 5); ! if (index) ! MULT(z, table[index], z) ! } ! } ! } ! ! if (doMontgomery) ! REDUCE(z) ! ! if (negativeOutput && (z->ob_size != 0)) { ! temp = (PyLongObject *)long_sub(z, c); ! if (temp == NULL) { ! Py_DECREF(z); ! z = NULL; ! goto error; ! } ! Py_DECREF(z); ! z = temp; ! temp = NULL; ! } ! ! error: ! Py_XDECREF(a); ! Py_XDECREF(b); ! Py_XDECREF(c); ! Py_XDECREF(temp); ! Py_XDECREF(R); ! if (b->ob_size > FIVEARY_CUTOFF) ! for (i=1; i < 32; ++i) { ! Py_XDECREF(table[i]); ! } ! return (PyObject *)z; } static PyObject *