/*
	Author:	James Pate Williams (c) 2001

	Exercise 1.3 from _Cryprography Theory and Practice_
	by Douglas R. Stinson p. 41. Originally, I had this
	problem incorrect. Martin Manscher of the Technical
	University of Denmark suggested the following
	solution.
*/

#include <math.h>
#include <stdio.h>

void extendedEuclid(long a, long b, long *x, long *y, long *d) {
	/* calculates a * *x + b * *y = gcd(a, b) = *d */
	long q, r, x1, x2, y1, y2;

	if (b == 0) {
		*d = a, *x = 1, *y = 0;
		return;
	}
	x2 = 1, x1 = 0, y2 = 0, y1 = 1;
	while (b > 0) {
		q = a / b, r = a - q * b;
		*x = x2 - q * x1, *y = y2 - q * y1;
		a = b, b = r;
		x2 = x1, x1 = *x, y2 = y1, y1 = *y;
	}
	*d = a, *x = x2, *y = y2;
}

long inverse(long a, long n) {
	/* computes the inverse of a modulo n */

	long d, x, y;
	
	extendedEuclid(a, n, &x, &y, &d);
	if (d == 1)
		return x;
	return 0;
}

int main(void) {
	int As11, As12, As21, As22, Asqr, detA, i, j, k, l, n = 26;
	int count1 = 0, count2 = 0, count3 = 0, count4 = 0;
	int count5 = 0, count6 = 0, count7 = 0, count8 = 0;
	long a11, m = 26;

	/*  
	   A = | i  j |
	       | k  l |

           det(A) = (i * l - j * k) mod 26

	   We want to find A such that det(A) = -+ 1 and A * A = 1
	*/
	for (i = 0; i < n; i++) {
		for (j = 0; j < n; j++) {
			for (k = 0; k < n; k++) {
				for (l = 0; l < n; l++) {
					detA = (i * l - j * k) % n;
					if (detA < 0)
						detA += 26;
					As11 = (i * i + j * k) % n;
					As12 = (i * j + j * l) % n;
					As21 = (k * i + l * k) % n;
					As22 = (k * j + l * l) % n;
					Asqr = As11 == 1 && As12 == 0 && As21 == 0 && As22 == 1;
					if (Asqr) {
						count1++;
						if (detA == n - 1) {
							count2++;
							a11 = (1 - i * i) % m;
							if (a11 < 0)
								a11 += m;
							if (inverse(a11, m) != 0)
								count3++;
							else if (a11 == 0)
								count4++;
							else if (a11 == 13)
								count5++;
							else if (a11 % 2 == 0)
								count6++;
						}
						else if (detA == 1)
							count7++;
					}
					else
						count8++;
				}
			}
		}
	}
	printf("number of involutary matrices with det(A) = + 1: %3d\n", count7);
	printf("number of involutary matrices with det(A) = - 1: %3d\n", count2);
	printf("det(A) = - 1, (1 - a[1][1] ^ 2) %% 26 in Z_26   : %3d\n", count3);
	printf("det(A) = - 1, (1 - a[1][1] ^ 2) %% 26 = 0       : %3d\n", count4);
	printf("det(A) = - 1, (1 - a[1][1] ^ 2) %% 26 = 13      : %3d\n", count5);
	printf("det(A) = - 1, (1 - a[1][1] ^ 2) %% 26 even      : %3d\n", count6);
	printf("total number of involutary matrices            : %3d\n", count1);
	printf("total number of matrices tested                : %3d\n", count1 + count8);
	printf("26 ** 4 = 26 ^ 4 = %f\n", pow(n, 4));
	return 0;
}