Mathematica PowerMod inverse and mpz_powm in C

Question

I have implemented an algorithm in Mathematica that uses PowerMod to find a modular inverse. I now need to implement this algorithm in C, and I've decided to use gmp and its function mpz_powm, which apparently does the same thing. The problem is, I'm not getting the same values. For instance, when run on Mathematica it gives:

PowerMod[30030, -1, 43] = 35

while mpz_pwm gives 16. And:

PowerMod[30030, -1, 71] = 8

whereas mpz_pwm gives 46. (I hope these are enough examples.) I have also tried various modular inverse algorithms I've managed to find online and they give different values too. I suspect, though, that Mathematica is right. Does anyone know what is happening here?

Isn't this really a question about C programming, not the software *Mathematica*? I am not sure this is the right site for the question. — Verbeia, Aug 13 '14 at 04:13
This question appears to be off-topic because it does not concern the technical computing software Mathematica by Wolfram Research. — , Aug 13 '14 at 05:27

Daniel Lichtblau · Answer 1 · 2014-08-13T15:10:57.553

There is something that seems relevant at this link.

In short, it appears that mpz_pwm interprets -1 as an unsigned short 2^16-1. One of the responses in the thread at this link is from a GMP developer, to the effect that this is a bug in the documentation.

As for the 2^16-1, I get that by some experimenting.

PowerMod[30030,2^16-1,43]                                               

(* Out[5]= 16 *)

So far so good. It falls apart when the third argument is 71. But PowerMod[30030, -1, 71] is not in fact 8. That holds when the third argument is 79. And we still do not get PowerMod[30030, 2^16-1,79] equal to 46: that holds when the third argument is 73. So I'm guessing that these were the arguments you had actually supplied to powerMod and mpz_pwm respectively.

If GMP does not have a function to do directly what you need, you can get the modular inverse from the extended GCD. I'll show the code one would use in Mathematica but I'm sure it can be adapted for the purpose at hand.

myModularInverse[n_,p_] := With[{inv=ExtendedGCD[n,p][[2,1]]},
  If[inv>0, inv, inv+p]]

Quick examples:

myModularInverse[30030, 43]                                            

(* Out[28]= 35 *)

myModularInverse[30030, 79]                                            

(* Out[29]= 8 *)

Edit: Some more digging indicates that GMP has a function mpz_invert that appears to do what is wanted here. See the GMP number theory functions. Also there is a gmp_invert function which I guess is at a higher level in the gmp hierarchy.

score 1 · Answer 2 · answered Aug 13 '14 at 16:21

The native GMP mpz_powm does not directly handle negative exponents. You need to calculate the inverse using mpz_invert.

I implemented the behavior you want in gmpy2 (a Python wrapper for GMP). My code is at:

https://code.google.com/p/gmpy/source/browse/trunk/src/gmpy2_pow.c

Note on the code: mpz_inoc and mpz_cloc are replacments for mpz_init and mpz_clear that cache objects behind the scenes.

Mathematica PowerMod inverse and mpz_powm in C

2 Answers2