logmod

Computes y such that mod(a^y, p^N) == x

Description

Suppose
   p is an odd prime,
   ord(a) = p-1
   GCD(x, p) = 1
   N is an integer >= 2

Then
   y = logmod(x, a, p, N)

    Returns y such that a^y = x (mod p^N)

Example
Let x= vpi(154), a = vpi(7), p = vpi(17), N = vpi(37)
y = logmod(x, a, p, N)
returns
   y = 2088349219044680767324467844670001776975183904

• MATLAB release: MATLAB 7.2 (R2006a)
• Other requirements: The function mfile modinv.m IS NOT required any more.

[John D'Errico]

Nifty. I had to fix a couple of lines to get this to work, but they were mostly my fault. When I first wrote the vpi toolbox, I had named the prime testing tool primality.m. While this made some sense to me at the time, it was a poor choice, as pointed out by Petter. So if you try to use this code with a new release of vpi, you will need to modify the code to use isprime instead of primality, as well as minv instead of modinv.

Regardless, this is great. It works nicely. Pretty quick too.

x = vpi(154);
a = vpi(7);
N = vpi(37);
p = vpi(17);

y = logmod(x,a,p,N)
y =
   2088349219044680767324467844670001776975183904

powermod(a,y,p^N)
ans =
   154

[Steven Gregory]

Thanks for the 5 stars. I wrote this program in Mathematica because I am doing something number theory -ish (for fun) that required it. I prefer Matlab as I think of Mathematica as a really expensive Swiss army knife. I really don't expect anyone else to need it but when I saw your VPI software I just had to try it out.

I don't have a lot of free time and I hadn't noticed that you updated the VPI software. I'm sorry I hadn't because it would have saved me some effort as well as cleaned up my code somewhat.

I promise as soon as the alligators subside a bit, I will update the program.

I have a Microsoft Word 2003 document that describes the algorithm if anyone is interested. I will add that to the upgrade as soon as I figure out how to do it.