function [ mInv, isInverse ] = matrixModularInverse( M, verbose)
%MATRIXMODULARINVERSE Computes the inverse of the matrix, mod 26.
% This function computes the inverse of the matrix, mod 26. This is useful
% when trying to find the decryption matrix for a Hill cipher.
% We denote the inverse as M^(-1)
% 
% The algorithm involves using Matlab's built-in inverse function 
% (which uses decimals), factoring out the determinant, 
% inverting it, and multiplying it back in. 

% Used as a tutorial for CSSE/MA479
% Matt Boutell, Spring 2009

% As an example, let's compute M^(-1) (mod 26)
% where M = |4 5| = |a b|
%           |1 2|   |c d|

% Simple: use Matlab's inverse function, right? 
% Only problem is that gives decimals (with no mod)
format rat % for viewing debug output

[rows, cols] = size(M);
if rows ~= cols
    fprintf('NO INVERSE! M is not square');
    isInverse = 0;
    return;
end

mInv = inv(M);
if verbose
    mInv
end


% In our example, this is: 
% |    2/3   -5/3 |
% |   -1/3    4/3 |
%
% However, we know that M^(-1) can always be written as 1/det(M) * (a matrix
% of integers, like the 2x2 case:
% M^(-1) = 1/(ad-bc) * | d -b |
%                      |-c  a |
%
% In our example, det = (4)(2) - (1)(5) = 3
% and we could factor out the 1/3 and rewrite M^(-1) as: 
% M^(-1) = 1/3 * | 2 -5 |
%                |-1  4 |

determinant = round(det(M));
if verbose
    determinant
end

if (gcd(determinant, 26) == 1)
    isInverse = 1;
else
    fprintf('NO INVERSE! gcd(det, 26)= %d\n', gcd(determinant,26))
    isInverse = 0;
end

% Factoring out 1/mInv to get integers. In the 2x2 case, we could use 
% the determinant formula, but that doesn't help with larger matrices.
mInvIntegers = mInv * determinant;
if verbose
    mInvIntegers
end


% Find determinant^(-1) mod 26 and multiply it back in.
% Note this only works when gcd(d, 26) == 1
% One nice thing is that gcd returns x and y from the extended 
% Euclidean algorithm. x is d^(-1) (mod 26)

% In our example, 3^(-1) = 9, so 
% M^(-1) = 9 * | 2 -5 | (mod 26) = | 18  7 |
%              |-1  4 |            | 17 10 |


[g, detInv, y]= gcd(determinant, 26); 
if verbose
    detInv
end

mInv = detInv * mInvIntegers;
if verbose
    mInv
end

mInv = mod(mInv, 26);

% Rounding on the next line fixes inaccuracies in doubles. 
% We still need to convert back to double once
% we round, so that we can multiply by another matrix.
mInv = double(round(mInv)); 
if verbose
    mInv
end

% Check by multiplying m * mInv; it should be the identity.
diffs = mod(M*mInv, 26) ~= eye(rows); 
if sum(sum(diffs)) ~= 0
    fprintf('Error! M * mInv is not I');
end
end
