def mod_inv(a, p):
    """Return x with x * a == 1 (mod p)  for integers a, p.

    p >= 1 must hold.
    Raise ZerodivsionError if inversion is not possible
    """
    assert p >= 1
    a1, a2, x1, x2 = p, a % p, 0, 1
    # Invert a modulo p with the extended Euclidean algorithm
    while a2 > 0:
        # assert x1 * a % p == a1 % p and  x2 * a % p == a2 % p
        q, r = divmod(a1, a2)
        a1, a2, x1, x2 = a2, r, x2, x1 - q * x2
    if a1 != 1:
        raise ZeroDivisionError("Cannot invert %d mod %d" % (a, p))
    # assert x1 * a % p == 1
    return x1 % p