# Constant time multiplicative inverse within a word

I was playing with an algorithm which at one step, calculated $f(x) = x^{-1} \mod p$ for $0 < x < p = 2^{64}-59$ (note $p$ is a prime). I used Knuth's Vol 2 Algorithm X algorithm for calculating inverse modulo a prime, using the Extended GCD, but I was wondering if there was a way for making it run in constant time, within a word?

The trick I found online for producing constant time inverses used Montgomery Multiplication, and therefore wouldn't fit within a 64-bit word.

My objective is to write a 64-bit cipher which (is probably insecure, definitely slow, and) uses inversion for confusion and diffusion.

Edit: Specified $p = 2^{64}-59$, a prime.

Edit2: Specified the algorithm attributed to Knuth with link to page which has attribution.

• Alternatively, would a blinded approach work (where the algorithm takes variable time, but time taken is uncorrelated to the value being inverted)? – poncho Mar 1 '16 at 16:14
• For what value of $p$? How large is your $p$? If $p$ is small you can use a table lookup. – D.W. Mar 1 '16 at 18:24
• $p = 2^{64} - 59$ as in linked code. – yberman Mar 1 '16 at 20:23
• @poncho, wrt a "blinded approach" the only think I could think of would to insert random multiples then strip them off? Calculate $(xr)^{-1}$ for random $r$, then multiply by $r$. I need to think about how to pull that off. – yberman Mar 1 '16 at 20:40

You can compute $f(x) = x^{\phi(p)-1} \bmod p$ in constant time, using $O(\log p)$ constant time modular multiplies. If $p$ is prime, this reduces to $f(x) = x^{p-2} \bmod p$

BTW: calling it 'Knuth's' algorithm is, in general, not very helpful; Knuth gives hundreds of different algorithms in 'The Art of Computer Programming'. I assume you mean the Extended Euclidean algorithm?

• That is correct. His Extended Euclidean algorithm. I am not sure how to take the an exponent if $2^{32} < p <2^{64}$, because products can overflow? That is why I didn't not use a Fermat's little theorem type based approach. – yberman Mar 1 '16 at 15:35
• @Yosef: simple, use Knuth's algorithm! TAOCP Volume 2, 4.6.3, Algorithm A (it's over 3400 years old), reducing modulo $p$ at each step. Or any of the following ones. – fgrieu Mar 1 '16 at 15:45
• @poncho: the algorithm I point to is the one to implement exponentiation (by right-to-left exponent scanning); so it is constant time (if modular multiplication and reduction is), and combined with your technique yields constant-time modular inverse. Designating this "Knuth's algorithm" was meant as an illustration of your note. – fgrieu Mar 1 '16 at 16:06
• @fgrieu: oops, sorry about that. In any case, his complaint was that it was inconvenient to do a modular multiply, given that his language doesn't have support for intermediate values > 64 bits – poncho Mar 1 '16 at 16:09

This answers the original question in the particular case of $p=2^{64}$. It quickly computes the modular inverse of any odd 64-bit integer, hopefully on any compilers with 64-bit support conforming to C99 or later. It is constant-time if basic operations are (which is often a worry for multiplication).

#include <stdint.h>                     // for uint64_t and uint32_t

// given odd a, compute x such that a*x = 1 over 64 bits.
uint64_t invmod(uint64_t a) {
uint32_t x = (((a+2u)&4u)<<1)+a;    // low  4 bits of inverse
x =    (2u-1u*a*x)*x;               // low  8 bits of inverse
x =    (2u-1u*a*x)*x;               // low 16 bits of inverse
x =    (2u-1u*a*x)*x;               // low 32 bits of inverse
return (2u-1u*a*x)*x;               //     64 bits of inverse
}


Justification: $a\cdot x\equiv1\pmod{2^k}\implies a\cdot(2-a\cdot x)\cdot x\equiv1\pmod{2^{2\cdot k}}$

Note: the (revised) code uses 2u and 1u* so that all quantities involved have an unsigned type. In C, unsigned quantities yield a well-defined behavior: the results are reduced modulo $2^k$ for some $k$, such that the unsigned type represents integers $x$ with $0\le x<2^k$.

• Thank you. I should have specified that $p = 2^{64} - 59$, the largest prime less than $2^{64}$. I failed to specify $p$ was prime. I thought it was implicit for the domain specification $0 < x < p$, but in retrospect, it is not. I also thought there would be a generic solution for all $p$, which makes the problem harder than a specific $p$. That being said, this is really slick. – yberman Mar 1 '16 at 20:29
• This code has undefined behavior on a platform that has a 64-bit int (most don't). Putting the letter u after each 0x2 solves this problem because the promotion to at least unsigned int will propagate to the multiplications, due to the left-to-right evaluation. See stackoverflow.com/questions/24795651/… – Myria Mar 1 '16 at 23:45
• @Myria: the code you commented works at least on a platform with 64-bit int using 2-complement that does not raise exception on integer underflow. But indeed the suffix ushall be used here. I originally posted 0x2in hope that was more recognizable than than 2u in my code, but I incorrectly remembered the rules governing signedness of hex constants without suffix (contrary to decimal ones, they can be unsigned for some intervals of value, where I wrongly remembered they allways are unsigned). Thanks for correcting me. – fgrieu Mar 2 '16 at 10:47
• I think b*x is still signed when int has 64 bits. – CodesInChaos Nov 24 '16 at 9:56
• @CodesInChaos: That was a good catch! For int of 64 bits, and band x of type uint32_t, indeed b*x is a (signed) int, potentially overflows, and is not formally specified, including in its low-order 32 bits. A generic solution to that issue of the C language is 1u*b*x, which really forces the use of unsigned arithmetic. I now have simplified the code, and hopefully it is OK. The 1u*a thing seems necessary for an hypothetical 80 or 128-bit int [re-posted with edit following revision of the answer]. – fgrieu Nov 25 '16 at 8:14