Ricerca chiavi cifrario di Hill

keyHill.c

// Coded by Pietro Squilla
// Semplifica l'uso del cifrario di Hill con carta e penna ricercando matrici con det = +-1;
// osservazione: se la matrice ha determinante 1 o -1 la sua inversa conterrà numeri interi
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

// grandezza alfabeto -> modulo base = MASSIMO + 1
#define MINIMO 0
#define MASSIMO 28

// parametro per la decomposizione LUP
#define TOL 0.001

double randReal(double min, double max);
double **createMatrix(long long int rows, long long int columns);
void printMatrix(long long int rows, long long int columns, double **M);
void freeMatrix(long long int rows, double **M);
double **transposedMatrix(long long int rows, long long int columns, double **M);
double **multiplyMatrices(long long int rows1, long long int columns1, long long int rows2, long long int columns2, double **M1, double **M2);
int LUPDecompose(double **A, int N, double Tol, int *P);
void LUPSolve(double **A, int *P, double *b, int N, double *x);
void LUPInvert(double **A, int *P, int N, double **IA);
double LUPDeterminant(double **A, int *P, int N);

int main(){
    long long int n, i, j;
    double determinante = 0.;
    srand48(time(NULL));
    
    printf("Inserisci la dimensione della matrice: ");
    scanf("%lld",&n);
    
    double **matrice = createMatrix(n,n);    // matrice su cui fare i calcoli
    double **matricetmp = createMatrix(n,n); // matrice backup da printare
    double **inversa = createMatrix(n,n);
    int *P = calloc(n+1,sizeof(int));
    
    fprintf(stderr,"Calcolo matrice in corso...\n");
    while(determinante != 1. && determinante != -1.){
        for(i=0; i<n; i++){
            for(j=0; j<n; j++){
                matrice[i][j] = (int)(randReal(MINIMO,MASSIMO));
                matricetmp[i][j] = matrice[i][j];
            }
        }
        
        LUPDecompose(matrice,n,TOL,P);
        determinante = LUPDeterminant(matrice,P,n);
    }
    fprintf(stderr,"\n");
    
    LUPInvert(matrice,P,n,inversa);
    
    // print matrice
    fprintf(stderr,"Matrice con det = %d:",(int)determinante);
    printMatrix(n,n,matricetmp);
    
    // print inversa
    fprintf(stderr,"\nMatrice inversa:");
    printMatrix(n,n,inversa);
    
    // print matrice congruente in modulo
    fprintf(stderr,"\nMatrice inversa congrente mod %d:",MASSIMO+1);
    for(i=0; i<n; i++){
        for(j=0; j<n; j++){
            inversa[i][j] = (int)(inversa[i][j]) % (MASSIMO + 1);
            
            if(inversa[i][j] < 0){
                inversa[i][j] = (MASSIMO + 1) + inversa[i][j];
            }
        }
    }
    printMatrix(n,n,inversa);
    
    // libero la memoria
    freeMatrix(n,matrice);
    freeMatrix(n,matricetmp);
    freeMatrix(n,inversa);
    free(P);
    
    return 0;
}


double randReal(double min, double max){
    double range = (max - min); 
    double denom = RAND_MAX / range;
    
    return (min + (lrand48() / denom));
}

double **createMatrix(long long int rows, long long int columns){
    long long int i;
    
    double **M = (double **) malloc(rows * sizeof(double *));
    if(M == NULL){
        perror("\nError");
        printf("\n");
        exit(2);
    }
    
    for(i=0; i<rows; i++){
        M[i] = (double *) malloc(columns * sizeof(double));
        if(M == NULL){
            perror("\nError");
            printf("\n");
            exit(3);
        }
    }
    
    return M;
}

void printMatrix(long long int rows, long long int columns, double **M){
    long long int i, j;
    
    printf("\n");
    
    for(i=0; i<rows; i++){
        for(j=0; j<columns; j++){
            printf("%d\t",(int)M[i][j]);
        }
        printf("\n");
    }

    return ;
}

void freeMatrix(long long int rows, double **M){

    while(--rows > -1){
        free(M[rows]);
    }
    
    free(M);
    
    return ;
}

double **transposedMatrix(long long int rows, long long int columns, double **M){
    long long int i, j;
    double **transposedMatrix = createMatrix(rows,columns);

    for(i=0; i<rows; i++){
        for(j=0; j<columns; j++){
            transposedMatrix[j][i] = M[i][j];
        }
    }
    
    return transposedMatrix;
}

double **multiplyMatrices(long long int rows1, long long int columns1, long long int rows2, long long int columns2, double **M1, double **M2){
    long long int i, j, k, sum=0;
    
    if(columns1 != rows2){
        printf("\nThe inserted matrix can not be multiplied!\n");
        exit(1);
    }
    
    double **product = createMatrix(columns1,rows2);

    for(i=0; i<rows1; i++){
        for(j=0; j<columns2; j++){
            for(k=0; k<rows2; k++){
                sum = sum + M1[i][k] * M2[k][j];
            }
            product[i][j] = sum;
            sum = 0;
        }
    }
    
    return product;
}

/*  INPUT: A - array of pointers to rows of a square matrix having dimension N
*          Tol - small tolerance number to detect failure when the matrix is near degenerate
*
*   OUTPUT: Matrix A is changed, it contains both matrices L-E and U as A=(L-E)+U such that P*A=L*U.
*           The permutation matrix is not stored as a matrix, but in an integer vector P of size N+1 
*           containing column indexes where the permutation matrix has "1". The last element P[N]=S+N, 
*           where S is the number of row exchanges needed for determinant computation, det(P)=(-1)^S    
*/
int LUPDecompose(double **A, int N, double Tol, int *P){
    int i, j, k, imax;
    double maxA, *ptr, absA;
    
    for(i=0; i<=N; i++){
        P[i] = i; // Unit permutation matrix, P[N] initialized with N
    }
    
    for(i=0; i<N; i++){
        maxA = 0.0;
        imax = i;
        
        for(k=i; k<N; k++){
            if((absA = fabs(A[k][i])) > maxA){ 
                maxA = absA;
                imax = k;
            }
        }
        
        if(maxA < Tol)
            return 0; // failure, matrix is degenerate
        
        if(imax != i){
            // pivoting P
            j = P[i];
            P[i] = P[imax];
            P[imax] = j;

            // pivoting rows of A
            ptr = A[i];
            A[i] = A[imax];
            A[imax] = ptr;

            // counting pivots starting from N (for determinant)
            P[N]++;
        }
        
        for(j=i+1; j<N; j++){
            A[j][i] /= A[i][i];
            
            for(k = i + 1; k < N; k++){
                A[j][k] -= A[j][i] * A[i][k];
            }
        }
    }
    
    return 1;  // decomposition done
}

/*  INPUT:  A, P filled in LUPDecompose; b - rhs vector; N - dimension
*   OUTPUT: x - solution vector of A*x=b
*/
void LUPSolve(double **A, int *P, double *b, int N, double *x){
    int i, k;
    
    for(i=0; i<N; i++){
        x[i] = b[P[i]];
        
        for(k=0; k<i; k++){
            x[i] -= A[i][k] * x[k];
        }
    }
    
    for(i=N-1; i>=0; i--){
        for(k=i + 1; k<N; k++){
            x[i] -= A[i][k] * x[k];
        }
        
        x[i] = x[i] / A[i][i];
    }
}

/*  INPUT:  A, P filled in LUPDecompose; N - dimension
*   OUTPUT: IA is the inverse of the initial matrix
*/
void LUPInvert(double **A, int *P, int N, double **IA){
    int i, j, k;
    
    for(j=0; j<N; j++){
        for(i=0; i<N; i++){
            if(P[i] == j){
                IA[i][j] = 1.0;
            }else{
                IA[i][j] = 0.0;
            }
            
            for(k=0; k<i; k++){
                IA[i][j] -= A[i][k] * IA[k][j];
            }
        }
        
        for(i=N-1; i>=0; i--){
            for(k=i+1; k<N; k++){
                IA[i][j] -= A[i][k] * IA[k][j];
            }
            IA[i][j] = IA[i][j] / A[i][i];
        }
    }
}

/*  INPUT:  A, P filled in LUPDecompose; N - dimension. 
*   OUTPUT: Function returns the determinant of the initial matrix
*/
double LUPDeterminant(double **A, int *P, int N){
    int i;
    double det = A[0][0];
    
    for(i=1; i<N; i++){
        det *= A[i][i];
    }
    
    if((P[N] - N) % 2 == 0){
        return det;
    }else{
        return -det;
    }
    
    // unpredicted failure
    exit(1);
}