egcd(_, 0, 1, 0) :- !.
egcd(A, B, X, Y) :-
    divmod(A, B, Q, R),
    egcd(B, R, S, X),
    Y is S - Q*X.

modinv(A, B, N) :-
    egcd(A, B, X, Y),
    A*X + B*Y =:= 1,
    N is X mod B.