Multiplicative inverse mod 2^w, without Euclid or Bezout.

multiplicative-inverse.lisp

(defun multiplicative-inverse (x bitwidth)
  "Iteratively finds the multiplicative inverse of x mod 2^bitwidth.

   The induction step uses the fact that, for a given inverse x^-1
   with (x * x^-1) mod 2^w == 1, (x * x^-1) mod 2^{w - 1} == 1 as
   well.  We start with the trivial inverse mod 2, and increase the
   bitwidth iteratively: given inv_{w - 1} = x^-1 mod 2^{w - 1},
   x^-1 mod 2^w is either inv_{w - 1}, or inv_{w - 1} + 2^{w - 1}."
  (assert (oddp x))
  (let ((inverse 1))
    ;; For any odd x, (x * 1) mod 2 == 1, so
    ;; we begin the search at mod (1 << 2) == mod 4.
    (loop for i from 2 upto bitwidth
       do (let ((mask (1- (ash 1 i))))
            (setf inverse (logior inverse
                                  (if (= 1 (logand (* inverse x) mask))
                                      0
                                      (ash 1 (1- i)))))
            (assert (= 1 (logand (* inverse x) mask))))
       finally
         (assert (= 1 (mod (* inverse x) (ash 1 bitwidth))))
         (return inverse))))