/*
  Author:  Pate Williams (c) 1997

  Exercise III.2.2
  "Find the inverses of the following matrices
  mod N. Write the entries in the inverse matrix
  as nonnegative integers less than N.
  (a) | 1 3 | mod 5    (b) | 1 3 | mod 29
      | 4 3 |              | 4 3 |
  (c) | 15 17 | mod 26 (d) | 40  0 | mod 841
      |  4  9 |            | 0  21 |
  (e) | 197  62 | mod 841
      | 603 271 |". -Neal Koblitz-
  See "A Course in Number Theory and Cryptography"
  by Neal Koblitz second edition page 77.
*/

#include <malloc.h>
#include <stdio.h>

long Extended_Euclidean(long b, long n)
{
  long b0 = b, n0 = n, t = 1, t0 = 0, temp, q, r;

  q = n0 / b0;
  r = n0 - q * b0;
  while (r > 0) {
    temp = t0 - q * t;
    if (temp >= 0) temp = temp % n;
    else temp = n - (- temp % n);
    t0 = t;
    t = temp;
    n0 = b0;
    b0 = r;
    q = n0 / b0;
    r = n0 - q * b0;
  }
  if (b0 != 1) return 0;
  else return t % n;
}

void inverse(long N, long **A)
{
  long d = (A[0][0] * A[1][1] - A[0][1] * A[1][0]) % N;
  long B[2][2], C[2][2], i, j;

  if (d < 0) d += N;
  d = Extended_Euclidean(d, N);
  B[0][0] = (d * A[1][1]) % N;
  B[1][1] = (d * A[0][0]) % N;
  B[0][1] = (- d * A[0][1]) % N;
  B[1][0] = (- d * A[1][0]) % N;
  C[0][0] = (A[0][0] * B[0][0] + A[0][1] * B[1][0]) % N;
  C[0][1] = (A[0][0] * B[0][1] + A[0][1] * B[1][1]) % N;
  C[1][0] = (A[1][0] * B[0][0] + A[1][1] * B[1][0]) % N;
  C[1][1] = (A[1][0] * B[0][1] + A[1][1] * B[1][1]) % N;
  for (i = 0; i < 2; i++) {
    for (j = 0; j < 2; j++)
      printf("%3ld ", A[i][j]);
    for (j = 0; j < 2; j++) {
      if (B[i][j] < 0) B[i][j] += N;
      printf("%3ld ", B[i][j]);
    }
    for (j = 0; j < 2; j++) {
      if (C[i][j] < 0) C[i][j] += N;
      printf("%3ld ", C[i][j]);
    }
    printf("\n");
  }
  printf("\n");
}

int main(void)
{
  long A[5][2][2] = {{{1, 3}, {4, 3}},
                     {{1, 3}, {4, 3}},
                     {{15, 17}, {4, 9}},
                     {{40, 0}, {0, 21}},
                     {{197, 62}, {603, 271}}};
  long N[5] = {5, 29, 26, 841, 841}, **B, i, j, k;

  B = calloc(2, sizeof(long *));
  for (i = 0; i < 2; i++)
    B[i] = calloc(2, sizeof(long));
  for (i = 0; i < 5; i++) {
    for (j = 0; j < 2; j++)
      for (k = 0; k < 2; k++)
        B[j][k] = A[i][j][k];
    inverse(N[i], B);
  }
  return 0;
}