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.

  • 2
    $\begingroup$ Alternatively, would a blinded approach work (where the algorithm takes variable time, but time taken is uncorrelated to the value being inverted)? $\endgroup$ – poncho Mar 1 '16 at 16:14
  • $\begingroup$ For what value of $p$? How large is your $p$? If $p$ is small you can use a table lookup. $\endgroup$ – D.W. Mar 1 '16 at 18:24
  • $\begingroup$ $p = 2^{64} - 59$ as in linked code. $\endgroup$ – yberman Mar 1 '16 at 20:23
  • 1
    $\begingroup$ @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. $\endgroup$ – 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?

  • $\begingroup$ 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. $\endgroup$ – yberman Mar 1 '16 at 15:35
  • $\begingroup$ @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. $\endgroup$ – fgrieu Mar 1 '16 at 15:45
  • $\begingroup$ @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. $\endgroup$ – fgrieu Mar 1 '16 at 16:06
  • 2
    $\begingroup$ @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 $\endgroup$ – 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$.

  • $\begingroup$ 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. $\endgroup$ – yberman Mar 1 '16 at 20:29
  • 1
    $\begingroup$ 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/… $\endgroup$ – Myria Mar 1 '16 at 23:45
  • $\begingroup$ @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. $\endgroup$ – fgrieu Mar 2 '16 at 10:47
  • 1
    $\begingroup$ I think b*x is still signed when int has 64 bits. $\endgroup$ – CodesInChaos Nov 24 '16 at 9:56
  • $\begingroup$ @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]. $\endgroup$ – fgrieu Nov 25 '16 at 8:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.