Patch: 3-argument pow() function

Andrew KUCHLING (fnord@cs.mcgill.ca)
3 Aug 1994 14:44:36 GMT

-----BEGIN PGP SIGNED MESSAGE-----

The following patch adds an optional third argument to the pow() builtin
function. pow(x,y) still computes x to the power of y, as before. But
pow(x,y,z) will now compute x to the power of y, modulo z. pow(x,y,z) is
*much* faster than pow(x,y) % z, because the modulo operation can be
performed after every intermediate operation, thus avoiding the need to
generate a very large result for pow(x,y).

It's mostly intended for use with long integers. For floating point
numbers, it just calls (fmod(pow(x,y), z). For ordinary integers, it should
provide a theoretical speedup over the existing exponentiation code, but you
can't raise integers to very large powers anyway, because the result won't
fit in an integer.

I wrote this because I'm working on a PGP module for Python. As an example,
here's an RSA encryption/decryption using keys previously generated by PGP:
(Sorry for the long lines.)

develop@quartz% python
>>> n= 34303374782200890257026213439194813524709428368411760108466060114374302355213380102758494564333732314126567943631081L
>>> e= 17L
>>> d= 1793640511487628248733396781134369334625329587890810986062971297361096118883133442448242634828621732589125755914257L
>>> def encrypt(data): # RSA encryption
... return pow(data, e, n)
...
>>> def decrypt(data): # RSA decryption
... return pow(data, d, n)
...
>>> v=encrypt(19721015)
>>> print v
26569145839974975265330455446932665188282825500235559542097848414057871021969872844013264870800518470658012072241940L
>>> decrypt(v)
19721015L
>>> decrypt(encrypt(123456789L))
123456789L
>>> encrypt(decrypt(9876543210L))
9876543210L

If you try this, you'll find decryption is still rather slow, due to the
repeated squarings of x. I'm still working on further optimizing the long
integer routines.

If you're interested in the modification, simply apply the patch and
recompile. The patch has been tried on Python 1.0.3, although it should
work on 1.0.2 as well, and it may work on even older versions. If you have
any problems, or (horrors!) find a bug or mathematical error, please let me
know. If no one reports any problems, after a while I'll ask GvR to make
the 3-argument pow() an official part of the language.

Andrew Kuchling
fnord@binkley.cs.mcgill.ca

diff -C2 ../origpy/Include//object.h Include//object.h
*** ../origpy/Include//object.h Tue Jan 11 11:22:02 1994
- --- Include//object.h Mon Aug 1 22:44:04 1994
***************
*** 136,139 ****
- --- 136,140 ----
typedef object * (*unaryfunc) PROTO((object *));
typedef object * (*binaryfunc) PROTO((object *, object *));
+ typedef object * (*trinaryfunc) PROTO((object *, object *, object *));
typedef int (*inquiry) PROTO((object *));
typedef int (*coercion) PROTO((object **, object **));
***************
*** 151,155 ****
binaryfunc nb_remainder;
binaryfunc nb_divmod;
! binaryfunc nb_power;
unaryfunc nb_negative;
unaryfunc nb_positive;
- --- 152,156 ----
binaryfunc nb_remainder;
binaryfunc nb_divmod;
! trinaryfunc nb_power;
unaryfunc nb_negative;
unaryfunc nb_positive;
diff -C2 ../origpy/Objects//floatobject.c Objects//floatobject.c
*** ../origpy/Objects//floatobject.c Tue Jan 4 22:04:56 1994
- --- Objects//floatobject.c Mon Aug 1 23:16:48 1994
***************
*** 262,272 ****

static object *
! float_pow(v, w)
floatobject *v;
floatobject *w;
{
double iv, iw, ix;
iv = v->ob_fval;
iw = w->ob_fval;
/* Sort out special cases here instead of relying on pow() */
if (iw == 0.0)
- --- 262,278 ----

static object *
! float_pow(v, w, z)
floatobject *v;
floatobject *w;
+ floatobject *z;
{
double iv, iw, ix;
iv = v->ob_fval;
iw = w->ob_fval;
+ /* XXX Doesn't handle overflows if z!=None yet; it may never do so :(
+ * The z parameter is really only going to be useful for integers and
+ * long integers. Maybe something clever with logarithms could be done.
+ * [AMK]
+ */
/* Sort out special cases here instead of relying on pow() */
if (iw == 0.0)
***************
*** 274,278 ****
if (iv == 0.0) {
if (iw < 0.0) {
! err_setstr(ValueError, "0.0 to the negative power");
return NULL;
}
- --- 280,284 ----
if (iv == 0.0) {
if (iw < 0.0) {
! err_setstr(ValueError, "0.0 to a negative power");
return NULL;
}
***************
*** 291,294 ****
- --- 297,301 ----
return NULL;
}
+ if (z!=None) ix=fmod(ix, z->ob_fval); /* XXX To Be Rewritten */
return newfloatobject(ix);
}
diff -C2 ../origpy/Objects//intobject.c Objects//intobject.c
*** ../origpy/Objects//intobject.c Tue Jan 4 22:05:01 1994
- --- Objects//intobject.c Mon Aug 1 23:15:00 1994
***************
*** 331,339 ****

static object *
! int_pow(v, w)
intobject *v;
intobject *w;
{
! register long iv, iw, ix;
iv = v->ob_ival;
iw = w->ob_ival;
- --- 331,341 ----

static object *
! int_pow(v, w, z)
intobject *v;
intobject *w;
+ intobject *z;
{
! register long iv, iw, iz, ix, temp, prev;
! char zset=0;
iv = v->ob_ival;
iw = w->ob_ival;
***************
*** 342,356 ****
return NULL;
}
ix = 1;
! while (--iw >= 0) {
! long prev = ix;
! ix = ix * iv;
! if (iv == 0)
! break; /* 0 to some power -- avoid ix / 0 */
! if (ix / iv != prev)
! return err_ovf("integer pow()");
}
return newintobject(ix);
! }

static object *
- --- 344,379 ----
return NULL;
}
+ if (z!=None) {
+ iz=z->ob_ival;
+ zset=1;
+ }
+ /*
+ * XXX: GvR's original exponentiation code stopped looping
+ * when temp hit zero; this code will continue onwards
+ * unnecessarily, but at least it won't cause any errors.
+ * Hopefully the speed improvement from the fast exponentiation
+ * will compensate for the slight inefficiency.
+ */
+ temp=iv;
ix = 1;
! while (iw > 0) {
! prev=ix; /* Save value for overflow check */
! if (iw & 1) {
! ix=ix*temp;
! if (temp == 0)
! break; /* Avoid ix / 0 */
! if (ix / temp != prev)
! return err_ovf("integer pow()");
! }
! prev=temp;
! temp *= temp; /* Square the value of temp */
! if(temp/prev!=prev)
! return err_ovf("integer pow()");
! if (zset) { /* If we did a multiplication, perform a modulo */
! ix = ix % iz;
! temp = temp % iz;
! }
! iw=iw>>1; /* Shift exponent down by 1 bit */
}
return newintobject(ix);
! }
!

static object *
***************
*** 542,546 ****
(binaryfunc)int_mod, /*nb_remainder*/
(binaryfunc)int_divmod, /*nb_divmod*/
! (binaryfunc)int_pow, /*nb_power*/
(unaryfunc)int_neg, /*nb_negative*/
(unaryfunc)int_pos, /*nb_positive*/
- --- 565,569 ----
(binaryfunc)int_mod, /*nb_remainder*/
(binaryfunc)int_divmod, /*nb_divmod*/
! (trinaryfunc)int_pow, /*nb_power*/
(unaryfunc)int_neg, /*nb_negative*/
(unaryfunc)int_pos, /*nb_positive*/
diff -C2 ../origpy/Objects//longobject.c Objects//longobject.c
*** ../origpy/Objects//longobject.c Tue Feb 15 16:01:36 1994
- --- Objects//longobject.c Tue Aug 2 22:45:54 1994
***************
*** 584,588 ****
static object *long_mod PROTO((longobject *, longobject *));
static object *long_divmod PROTO((longobject *, longobject *));
! static object *long_pow PROTO((longobject *, longobject *));
static object *long_neg PROTO((longobject *));
static object *long_pos PROTO((longobject *));
- --- 584,588 ----
static object *long_mod PROTO((longobject *, longobject *));
static object *long_divmod PROTO((longobject *, longobject *));
! static object *long_pow PROTO((longobject *, longobject *, longobject *));
static object *long_neg PROTO((longobject *));
static object *long_pos PROTO((longobject *));
***************
*** 953,961 ****

static object *
! long_pow(a, b)
longobject *a;
longobject *b;
{
! longobject *z;
int size_b, i;

- --- 953,962 ----

static object *
! long_pow(a, b, c)
longobject *a;
longobject *b;
+ longobject *c;
{
! longobject *z, *div, *mod;
int size_b, i;

***************
*** 965,983 ****
return NULL;
}
- -
z = (longobject *)newlongobject(1L);
- -
INCREF(a);
for (i = 0; i < size_b; ++i) {
digit bi = b->ob_digit[i];
int j;
!
for (j = 0; j < SHIFT; ++j) {
longobject *temp;
!
if (bi & 1) {
temp = (longobject *)long_mul(z, a);
DECREF(z);
! z = temp;
if (z == NULL)
break;
- --- 966,988 ----
return NULL;
}
z = (longobject *)newlongobject(1L);
INCREF(a);
for (i = 0; i < size_b; ++i) {
digit bi = b->ob_digit[i];
int j;
!
for (j = 0; j < SHIFT; ++j) {
longobject *temp;
!
if (bi & 1) {
temp = (longobject *)long_mul(z, a);
DECREF(z);
! if ((object*)c!=None && temp!=NULL) {
! l_divmod(temp, c, &div, &mod);
! XDECREF(div);
! DECREF(temp);
! temp=mod;
! }
! z=temp;
if (z == NULL)
break;
***************
*** 988,991 ****
- --- 993,1002 ----
temp = (longobject *)long_mul(a, a);
DECREF(a);
+ if ((object*)c!=None && temp!=NULL) {
+ l_divmod(temp, c, &div, &mod);
+ XDECREF(div);
+ DECREF(temp);
+ temp=mod;
+ }
a = temp;
if (a == NULL) {
***************
*** 1352,1355 ****
- --- 1363,1367 ----
#define UF (unaryfunc)
#define BF (binaryfunc)
+ #define TF (trinaryfunc)
#define IF (inquiry)

***************
*** 1361,1365 ****
BF long_mod, /*nb_remainder*/
BF long_divmod, /*nb_divmod*/
! BF long_pow, /*nb_power*/
UF long_neg, /*nb_negative*/
UF long_pos, /*tp_positive*/
- --- 1373,1377 ----
BF long_mod, /*nb_remainder*/
BF long_divmod, /*nb_divmod*/
! TF long_pow, /*nb_power*/
UF long_neg, /*nb_negative*/
UF long_pos, /*tp_positive*/
diff -C2 ../origpy/Python//bltinmodule.c Python//bltinmodule.c
*** ../origpy/Python//bltinmodule.c Fri Jul 15 18:28:39 1994
- --- Python//bltinmodule.c Tue Aug 2 22:00:16 1994
***************
*** 821,839 ****
object *args;
{
! object *v, *w, *x;
! if (!getargs(args, "(OO)", &v, &w))
! return NULL;
! if (v->ob_type->tp_as_number == NULL ||
w->ob_type->tp_as_number == NULL) {
err_setstr(TypeError, "pow() requires numeric arguments");
return NULL;
}
! if (coerce(&v, &w) != 0)
return NULL;
! x = (*v->ob_type->tp_as_number->nb_power)(v, w);
! DECREF(v);
! DECREF(w);
return x;
}

static object *
- --- 821,853 ----
object *args;
{
! object *v, *w, *z, *x;
! z=None;
! if (!getargs(args, "(OO)", &v, &w)) {
! err_clear();
! if (!getargs(args, "(OOO)", &v, &w, &z)) {
! return NULL;
! }
! }
! if (v->ob_type->tp_as_number == NULL ||
! (z!=None && z->ob_type->tp_as_number == NULL) ||
w->ob_type->tp_as_number == NULL) {
err_setstr(TypeError, "pow() requires numeric arguments");
return NULL;
}
! if (coerce(&v, &w) != 0)
return NULL;
! if (z!=None) { /* XXX Are all these coercions needed? [AMK] */
! if (coerce(&w, &z) != 0)
! return NULL;
! if (coerce(&v, &z) != 0)
! return NULL;
! }
! x = (*v->ob_type->tp_as_number->nb_power)(v, w, z);
! DECREF(v);
! DECREF(w);
! if (z!=None) {DECREF(w); DECREF(v); DECREF(z); DECREF(z);}
return x;
}
+

static object *

-----BEGIN PGP SIGNATURE-----
Version: 2.6

iQBUAgUBLj+s+kEKhMyDwxU1AQGiOQH2JPwrScknHbMPnxsgD2k5/FgCKEszGAPN
W8LhnOukkjuZzTq0iaivG3wcTI5gp5X9v+pOIQdD1UndQsimHsAB
=E8b8
-----END PGP SIGNATURE-----