2

I have created a function to compute large modular exponents. I am aware that this function is built into the python language. My function is incorrect for numbers with larger than 17 digits, and I can't figure out why. Any help is greatly appreciated.

from random import randint

def modpow(x,n,m):
  if n == 0:
    return 1
  elif n == 1:
    return x
  elif n%2 == 0:
    return modpow(x*(x%m),n/2,m)%m
  elif n%2 == 1:
    return (x *  modpow(x*(x%m),(n-1)/2,m)%m )%m

for i in range(5,32):
  x = ''
  n = ''
  m = ''
  for j in range(i):
    x += str(randint(0,9))
    n += str(randint(0,9))
    m += str(randint(0,9))
  x = int(x)
  n = int(n)
  m = int(m)
  if pow(x,n,m) != modpow(x,n,m):
    print(i,x,n,m)

Sample output:

17 38508450670424585 67111951647554134 59005802612594983
18 24027200512104766 205942669690724726 816654795945860553
...

I've run this several times and it always starts failing at i = 17, I'm not quite sure why.

Matt Pennington
  • 560
  • 1
  • 8
  • 21
  • 4
    If you are using Python 3, the `/` operator returns a floating point value. You probably want to use `//` to return an integer result with both Python 2 and 3. – casevh Apr 05 '14 at 22:19
  • Just ran the code with `from random import *` and `m = int(m)` and it seems to do well. EDIT: Python 2.7 – etna Apr 05 '14 at 22:25
  • I've made the suggested edits above. I am using Py3. – Matt Pennington Apr 05 '14 at 22:33
  • Replacing the / with // worked, thanks casevh! I must have been losing some precision when it converted n/2 to a double. – Matt Pennington Apr 05 '14 at 22:34
  • ``int`` in python is 'unlimited' precision. ``double`` is only 64bit, so yeah that was probably your problem. – aruisdante Apr 05 '14 at 22:35

1 Answers1

0

My guess is that Python's pow() works on doubles under the hood. The log base 2 of the first value that fails (38508450670424585) is about 55, but a double has only 53 bits of precision, so an integer that is larger than 2^53 can't necessarily be represented exactly. My guess is that if you check the log base 2 of the last number that succeeds (i.e. for x=16), it will be less than 53.

Emmet
  • 6,192
  • 26
  • 39