Index calculus to solve matrix in reduced row echelon form

Solve.c

//Full guide: https://leetarxiv.substack.com/p/row-reduction-over-finite-fields
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

// Function to compute modular inverse using Fermat's Little Theorem
long modinv(long a, long p)
{
    long result = 1;
    long power = p - 2;
    a = a % p;
    
    while (power > 0)
    {
        if (power % 2 == 1)
            result = (result * a) % p;
        a = (a * a) % p;
        power /= 2;
    }
    return result;
}

// Function to print a matrix
void printMatrix(long **matrix, int numRows, int numCols)
{
    for (int i = 0; i < numRows; i++)
    {
        for (int j = 0; j < numCols; j++)
        {
            printf("%ld ", matrix[i][j]);
        }
        printf("\n");
    }
}

// Function to perform RREF modulo p
void rref_mod(long **matrix, int numRows, int numCols, long p)
{
    int i = 0, j = 0;
    
    while (i < numRows && j < numCols)
    {
        // Find the pivot row
        if (matrix[i][j] == 0)
        {
            for (int k = i + 1; k < numRows; k++)
            {
                if (matrix[k][j] != 0)
                {
                    // Swap rows
                    long *temp = matrix[i];
                    matrix[i] = matrix[k];
                    matrix[k] = temp;
                    break;
                }
            }
        }

        if (matrix[i][j] == 0)
        {
            j++;
            continue;
        }

        // Normalize the pivot row
        long pivot = matrix[i][j];
        long inv_pivot = modinv(pivot, p);
        for (int col = 0; col < numCols; col++)
        {
            matrix[i][col] = (matrix[i][col] * inv_pivot) % p;
        }

        // Eliminate other rows
        for (int k = 0; k < numRows; k++)
        {
            if (k != i && matrix[k][j] != 0) 
            {
                long factor = matrix[k][j];
                for (int col = 0; col < numCols; col++)
                {
                    matrix[k][col] = (matrix[k][col] - factor * matrix[i][col]) % p;
                    // Ensure positive result
                    if (matrix[k][col] < 0) matrix[k][col] += p;
                }
            }
        }

        i++;
        j++;
    }
}

int main() {
    // Initialize matrix
    int numRows = 5;
    int numCols = 6;
    long p = 7043;
    
    // Allocate memory for matrix
    long **A = (long **)malloc(numRows * sizeof(long *));
    for (int i = 0; i < numRows; i++) {
        A[i] = (long *)malloc(numCols * sizeof(long));
    }
    
    // Initialize matrix values
    long values[5][6] = {
        {3, 0, 1, 2, 0, 346},
        {6, 0, 0, 1, 1, 171},
        {0, 3, 1, 0, 1, 153},
        {0, 4, 1, 0, 0, 442},
        {3, 5, 1, 0, 0, 458}
    };
    
    for (int i = 0; i < numRows; i++) {
        for (int j = 0; j < numCols; j++) {
            A[i][j] = values[i][j];
        }
    }
    
    printf("Original matrix:\n");
    printMatrix(A, numRows, numCols);
    
    // Perform RREF
    rref_mod(A, numRows, numCols, p);
    
    printf("\nRREF matrix (mod %ld):\n", p);
    printMatrix(A, numRows, numCols);
    
    // Free memory
    for (int i = 0; i < numRows; i++) {
        free(A[i]);
    }
    free(A);
    
    return 0;
}

Solve.py

#!/usr/bin/python

import time, random

def modinv(a, p):
    # Fermat's Little Theorem
    return pow(a, p-2, p)
    
def printMatrix(m):
    for row in m:
        print(str(row))



def rref_mod(matrix, p):
    if not matrix: return

    numRows = len(matrix)
    numCols = len(matrix[0])

    i, j = 0, 0
    while i < numRows and j < numCols:
        # Find the pivot row (swap if necessary)
        if matrix[i][j] == 0:
            for k in range(i+1, numRows):
                if matrix[k][j] != 0:
                    matrix[i], matrix[k] = matrix[k], matrix[i]
                    break

        if matrix[i][j] == 0:
            j += 1
            continue

        # Normalize the pivot row (using modular inverse)
        pivot = matrix[i][j]
        inv_pivot = modinv(pivot, p)
        matrix[i] = [(x * inv_pivot) % p for x in matrix[i]]

        # Eliminate other rows
        for k in range(numRows):
            if k != i and matrix[k][j] != 0:
                factor = matrix[k][j]
                matrix[k] = [(y - factor * x) % p 
                            for x, y in zip(matrix[i], matrix[k])]

        i += 1
        j += 1

    return matrix


A = [
    [3, 0, 1, 2, 0, 346],
    [6, 0, 0, 1, 1, 171],
    [0, 3, 1, 0, 1, 153],
    [0, 4, 1, 0, 0, 442],
    [3, 5, 1, 0, 0, 458]
]

p = 7043  # Prime modulus
A_rref = rref_mod(A, p)
printMatrix(A_rref)