Skip to content

Instantly share code, notes, and snippets.

@MurageKibicho

MurageKibicho/FMPZ_IndexCalculus.c

Last active October 18, 2025 14:01
Show Gist options
  • Save MurageKibicho/92b7cbad071f234dd01e0e84db704629 to your computer and use it in GitHub Desktop.
Save MurageKibicho/92b7cbad071f234dd01e0e84db704629 to your computer and use it in GitHub Desktop.
Download ZIP
Big Integer Index Calculus till RREF using FLINT and LibGMP in C. Needs FLINT 3.4 or greater for fmpz_mod_mat
Raw
FMPZ_IndexCalculus.c
////Complete walkthrough: https://leetarxiv.substack.com/p/row-reduction-over-finite-fields
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <stdbool.h>
#include <time.h>
#include <math.h>
#include <gmp.h>
#include <flint/flint.h>
#include <flint/fmpz.h>
#include <flint/fmpz_factor.h>
#include <flint/fmpz_mod_mat.h>
#include <flint/fmpz_mod.h>
#define STB_DS_IMPLEMENTATION
#include "stb_ds.h"
//clear && gcc FMPZ_IndexCalculus.c -lm -lflint -lgmp -lmpfr -o m.o && ./m.o
struct system_struct
{
fmpz_t prime;
fmpz_t primeMinusOne;
fmpz_t primitiveRoot;
fmpz_t temp0;
fmpz_t temp1;
flint_rand_t randomState;
fmpz_factor_t factorization;
fmpz_mod_ctx_t moduloContext;
fmpz_mod_mat_t rrefSolver;
fmpz_mod_mat_t relations;
fmpz_t **crtSolutions;
fmpz_t *crtDividends;
fmpz_t *crtInverses;
fmpz_t *combinedCRT;
int smoothnessBound;
int factorBaseLength;
int extraRelations;
int *factorBase;
int *primitiveRoots;
slong *rank;
};
typedef struct system_struct System[1];
double OdlyzkoBound(int N)
{
double logBase2 = (double) N;
double c = 1 / (4 * log(2));
c = sqrt(c);
double B = c * sqrt(logBase2 * log(logBase2));
B = pow(2, B);
return B;
}
bool FindGenerator(fmpz_t testGenerator, fmpz_t prime, fmpz_t primeMinusOne, fmpz_factor_t factorization)
{
bool result = false;
fmpz_t testResult, exponent;
fmpz_init(testResult);
fmpz_init(exponent);
for(slong i = 0; i < factorization->num; i++)
{
fmpz_divexact(exponent, primeMinusOne, &factorization->p[i]);
fmpz_powm(testResult, testGenerator, exponent, prime);
if(fmpz_cmp_ui(testResult, 1) == 0)
{
result = false;
fmpz_clear(exponent);
fmpz_clear(testResult);
return result;
}
}
fmpz_clear(exponent);
fmpz_clear(testResult);
result = true;
return result;
}
void System_CreateFactorBase(System system)
{
mpz_t prime,current;
mpz_inits(prime,current,NULL);
int primeHolder = 0;
int count = 0;
mpz_set_ui(prime,1);
if(system->factorBaseLength == 0)
{
while(primeHolder < system->smoothnessBound)
{
mpz_nextprime(prime, prime);
primeHolder = mpz_get_ui(prime);
arrput(system->factorBase, primeHolder);
system->factorBaseLength += 1;
}
system->smoothnessBound = system->factorBase[system->factorBaseLength - 1];
system->factorBaseLength = arrlen(system->factorBase);
}
else
{
while(count < system->factorBaseLength)
{
mpz_nextprime(prime, prime);
primeHolder = mpz_get_ui(prime);
system->smoothnessBound = primeHolder;
arrput(system->factorBase, primeHolder);
count += 1;
}
}
assert(arrlen(system->factorBase) == system->factorBaseLength);
mpz_clears(prime,current,NULL);
}
void System_InitFromInteger(System system, int prime, int smoothnessBound, int primitiveRootCount, int factorBaseLength, int extraRelations)
{
assert(prime > 0);
assert(primitiveRootCount > 0);
assert(extraRelations > 0);
system->smoothnessBound = (int)ceil(OdlyzkoBound(log2(prime)));
if(smoothnessBound > 0){system->smoothnessBound = smoothnessBound;}
system->factorBaseLength = factorBaseLength;
system->extraRelations = extraRelations;
system->factorBase = NULL;
system->primitiveRoots = NULL;
fmpz_init(system->prime);
fmpz_init(system->primeMinusOne);
fmpz_init(system->primitiveRoot);
fmpz_init(system->temp0);
fmpz_init(system->temp1);
flint_rand_init(system->randomState);
fmpz_factor_init(system->factorization);
fmpz_set_ui(system->prime, prime);
fmpz_sub_ui(system->primeMinusOne,system->prime, 1);
fmpz_factor(system->factorization, system->primeMinusOne);
System_CreateFactorBase(system);
fmpz_mod_ctx_init(system->moduloContext, system->primeMinusOne);
fmpz_mod_mat_init(system->relations, arrlen(system->factorBase) + system->extraRelations, arrlen(system->factorBase) + 1, system->moduloContext);
fmpz_mod_mat_init(system->rrefSolver, arrlen(system->factorBase) + system->extraRelations, arrlen(system->factorBase) + 1, system->moduloContext);
//Find primitive roots
bool isPrimitiveRoot = false;
for(int i = 3; i < 10000; i++)
{
fmpz_set_ui(system->primitiveRoot, i);
isPrimitiveRoot = FindGenerator(system->primitiveRoot, system->prime, system->primeMinusOne, system->factorization);
if(isPrimitiveRoot)
{
arrput(system->primitiveRoots, i);
primitiveRootCount -= 1;
if(primitiveRootCount == 0){break;}
}
}
//Create rank and crt holders
system->rank = calloc(system->factorization->num, sizeof(slong));
system->crtSolutions = malloc(system->factorization->num * sizeof(fmpz_t *));
for(slong i = 0; i < system->factorization->num; i++)
{
system->crtSolutions[i] = malloc(arrlen(system->factorBase) * sizeof(fmpz_t));
for(int j = 0; j < arrlen(system->factorBase); j++)
{
fmpz_init(system->crtSolutions[i][j]);
}
}
system->combinedCRT = malloc(arrlen(system->factorBase) * sizeof(fmpz_t));
for(int j = 0; j < arrlen(system->factorBase); j++)
{
fmpz_init(system->combinedCRT[j]);
}
system->crtDividends = malloc(system->factorization->num * sizeof(fmpz_t));
system->crtInverses = malloc(system->factorization->num * sizeof(fmpz_t));
for(int i = 0; i < system->factorization->num; i++)
{
fmpz_init(system->crtDividends[i]);
fmpz_init(system->crtInverses[i]);
fmpz_divexact(system->crtDividends[i], system->primeMinusOne, &system->factorization->p[i]);
fmpz_invmod(system->crtInverses[i], system->crtDividends[i],&system->factorization->p[i]);
}
}
void CollectDiscreteLogRelations(System system, int primitiveRootIndex)
{
int trialsMax = 10000;
int trials = trialsMax;
int found = 0;
assert(primitiveRootIndex > -1);
assert(primitiveRootIndex < arrlen(system->primitiveRoots));
fmpz_set_ui(system->primitiveRoot, system->primitiveRoots[primitiveRootIndex]);
fmpz_mod_ctx_set_modulus(system->moduloContext, system->primeMinusOne);
while(found < arrlen(system->factorBase) + system->extraRelations && trialsMax > -1)
{
bool foundSmooth = false;
while(!foundSmooth && trials > -1)
{
fmpz_randm(system->temp1, system->randomState, system->primeMinusOne);
fmpz_mod_set_fmpz(fmpz_mod_mat_entry(system->relations, found, arrlen(system->factorBase)), system->temp1, system->moduloContext);
fmpz_powm(system->temp0, system->primitiveRoot, system->temp1, system->prime);
for(size_t j = 0; j < arrlen(system->factorBase); j++)
{
int power = 0;
while(fmpz_mod_ui(system->temp1,system->temp0, system->factorBase[j]) == 0)
{
fmpz_set_ui(system->temp1, system->factorBase[j]);
fmpz_divexact(system->temp0, system->temp0, system->temp1);
power += 1;
}
fmpz_mod_set_ui(fmpz_mod_mat_entry(system->relations, found, j), power, system->moduloContext);
}
if(fmpz_cmp_ui(system->temp0, 1) == 0)
{
foundSmooth = true;
found += 1;
break;
}
trials -= 1;
}
trialsMax -= 1;
}
if(found != arrlen(system->factorBase) + system->extraRelations )
{
assert(found == arrlen(system->factorBase) + system->extraRelations);
}
}
void System_FindRREF(System system)
{
for(slong i = 0; i < system->factorization->num; i++)
{
fmpz_mod_ctx_set_modulus(system->moduloContext, system->primeMinusOne);
fmpz_mod_mat_set_fmpz_mat(system->rrefSolver, system->relations, system->moduloContext);
fmpz_mod_ctx_set_modulus(system->moduloContext, &system->factorization->p[i]);
_fmpz_mod_mat_reduce(system->rrefSolver, system->moduloContext);
system->rank[i] = fmpz_mod_mat_rref(system->rrefSolver, system->rrefSolver, system->moduloContext);
if(system->rank[i] != arrlen(system->factorBase))
{
printf("Maybe increase no. of relations(Rank: %lu, FactorBase: %lu)\n", system->rank[i], arrlen(system->factorBase));
assert(system->rank[i] == arrlen(system->factorBase));
}
else
{
for(size_t j = 0; j < arrlen(system->factorBase); j++)
{
fmpz_set(system->crtSolutions[i][j], fmpz_mod_mat_entry(system->rrefSolver, j, arrlen(system->factorBase)));
}
}
}
//Combine solutions without lifting (assumes all powers are 1)
for(size_t i = 0; i < arrlen(system->factorBase); i++)
{
fmpz_set_ui(system->combinedCRT[i], 0);
for(slong j = 0; j < system->factorization->num; j++)
{
fmpz_mul(system->temp0, system->crtInverses[j], system->crtDividends[j]);
fmpz_mul(system->temp0, system->temp0, system->crtSolutions[j][i]);
fmpz_add(system->combinedCRT[i], system->combinedCRT[i], system->temp0);
}
fmpz_mod(system->combinedCRT[i], system->combinedCRT[i], system->primeMinusOne);
}
}
void System_Clear(System system)
{
fmpz_clear(system->prime);
fmpz_clear(system->primeMinusOne);
fmpz_clear(system->primitiveRoot);
fmpz_clear(system->temp0);
fmpz_clear(system->temp1);
flint_rand_clear(system->randomState);
fmpz_factor_clear(system->factorization);
fmpz_mod_mat_clear(system->relations, system->moduloContext);
fmpz_mod_mat_clear(system->rrefSolver, system->moduloContext);
fmpz_mod_ctx_clear(system->moduloContext);
for(int j = 0; j < arrlen(system->factorBase); j++)
{
fmpz_clear(system->combinedCRT[j]);
}
free(system->combinedCRT);
arrfree(system->factorBase);
arrfree(system->primitiveRoots);
free(system->rank);
for(slong i = 0; i < system->factorization->num; i++)
{
for(int j = 0; j < arrlen(system->factorBase); j++)
{
fmpz_clear(system->crtSolutions[i][j]);
}
free(system->crtSolutions[i]);
}
free(system->crtSolutions);
for(int i = 0; i < system->factorization->num; i++)
{
fmpz_clear(system->crtDividends[i]);
fmpz_clear(system->crtInverses[i]);
}
free(system->crtDividends);
free(system->crtInverses);
}
void System_PrintIntArray(int length, int *array)
{
for(int i = 0; i < length; i++)
{
printf("%3d, ", array[i]);
}
printf("\n");
}
void PrintSystem(System system)
{
printf("Prime: ");fmpz_print(system->prime);printf("\n");
printf("Prime-1: ");fmpz_print(system->primeMinusOne);printf("\n");
printf("Factorization: ");
for(slong i = 0; i < system->factorization->num; i++)
{
printf("(");
fmpz_print(&system->factorization->p[i]); printf(" ^ ");
fmpz_print(&system->factorization->exp[i]); printf(")");
}
printf("\n");
printf("Primitive Roots:[");
for(size_t i = 0; i < arrlen(system->primitiveRoots); i++)
{
printf("%3d, ",system->primitiveRoots[i]);
}
printf("]\n");
printf("SmoothnessBound: %d\n",system->smoothnessBound);
printf("Factor Base (length %d): [", system->factorBaseLength);
for(size_t i = 0; i < arrlen(system->factorBase); i++)
{
printf("%3d, ",system->factorBase[i]);
}
printf("]\n");
printf("Current Primitive root: ");fmpz_print(system->primitiveRoot);printf("\n");
printf("Relations\n");
for(size_t i = 0; i < arrlen(system->factorBase) + system->extraRelations; i++)
{
printf("[");
for(size_t j = 0; j < arrlen(system->factorBase) + 1; j++)
{
fmpz_print(fmpz_mod_mat_entry(system->relations, i, j));printf(", ");
}
printf("],\n");
}
printf("CRT Solutions\n");
for(slong i = 0; i < system->factorization->num; i++)
{
printf("(");fmpz_print(&system->factorization->p[i]); printf(" ^ ");fmpz_print(&system->factorization->exp[i]); printf(") | ");
for(size_t j = 0; j < arrlen(system->factorBase); j++)
{
fmpz_print(system->crtSolutions[i][j]);printf(", ");
}
printf("\n");
}
printf("Final solution\n");
for(size_t i = 0; i < arrlen(system->factorBase); i++)
{
fmpz_print(system->combinedCRT[i]);printf(", ");
}
printf("\n");
}
void TestSystem()
{
System chap3;
int prime = 100003;
int smoothnessBound = 100;
int primitiveRootCount = 20;
int factorBaseLength = 10;
int extraRelations = 25;
int primitiveRootIndex = 5;
System_InitFromInteger(chap3, prime, smoothnessBound, primitiveRootCount,factorBaseLength, extraRelations);
CollectDiscreteLogRelations(chap3, primitiveRootIndex);
System_FindRREF(chap3);
PrintSystem(chap3);
System_Clear(chap3);
}
int main()
{
TestSystem();
flint_cleanup();
return 0;
}
Raw
Sage_Python_RREF.py
# Complete SageMath solution for RREF modulo composite numbers
# Now works for any composite modulus
def hensel_lift_matrix(M_p, p, k):
"""
Lift a matrix from mod p to mod p^k using proper Hensel lifting
"""
rows, cols = M_p.dimensions()
# Convert to integers for lifting
M_int = matrix(ZZ, [[int(x) for x in row] for row in M_p.rows()])
# For RREF, we need to be careful about pivots
# This is a simplified approach - in practice, you'd solve the system properly
M_pk = copy(M_int)
# Reduce modulo p^k
M_pk_mod = M_pk.change_ring(IntegerModRing(p^k))
return M_pk_mod
def modular_rref_composite(M, modulus):
"""
Compute RREF modulo a composite number by factoring and CRT
"""
print(f"Computing RREF modulo {modulus}")
# Factor the modulus
factors = factor(modulus)
print(f"Modulus factors: {list(factors)}")
moduli_solutions = {}
for p, k in factors:
p_k = p^k
print(f"\n=== Solving modulo {p}^{k} = {p_k} ===")
if k == 1:
# Prime modulus - use GF(p)
try:
M_p = Matrix(GF(p), [[int(x) for x in row] for row in M.rows()])
rref_p = M_p.rref()
moduli_solutions[p_k] = rref_p
print(f"RREF mod {p}:")
print(rref_p)
print(f"Rank mod {p}: {rref_p.rank()}")
except Exception as e:
print(f"Error solving mod {p}: {e}")
# Fallback: use echelon form
M_p = Matrix(IntegerModRing(p), [[int(x) for x in row] for row in M.rows()])
rref_p = M_p.echelon_form()
moduli_solutions[p_k] = rref_p
print(f"Echelon form mod {p} (fallback):")
print(rref_p)
else:
# Prime power modulus
try:
# Try direct echelon form for prime powers
M_pk = Matrix(IntegerModRing(p_k), [[int(x) for x in row] for row in M.rows()])
rref_pk = M_pk.echelon_form()
moduli_solutions[p_k] = rref_pk
print(f"Echelon form mod {p}^{k}:")
print(rref_pk)
except Exception as e:
print(f"Error solving mod {p}^{k}: {e}")
# Fallback: use Hensel lifting from mod p
M_p = Matrix(GF(p), [[int(x) for x in row] for row in M.rows()])
rref_p = M_p.rref()
print(f"RREF mod {p} (base for lifting):")
print(rref_p)
# Lift to mod p^k
rref_pk = hensel_lift_matrix(rref_p, p, k)
moduli_solutions[p_k] = rref_pk
print(f"Lifted RREF mod {p}^{k}:")
print(rref_pk)
# Combine using CRT
if len(moduli_solutions) == 0:
print("No solutions found for any modulus!")
return None
print(f"\n=== Combining {len(moduli_solutions)} solutions using CRT ===")
final_rref = combine_crt_solutions(moduli_solutions, modulus)
return final_rref
def combine_crt_solutions(solutions_dict, modulus):
"""
Combine solutions from different moduli using CRT
"""
moduli = list(solutions_dict.keys())
if len(moduli) == 1:
return solutions_dict[moduli[0]]
# For matrix CRT, we combine element-wise
sample_sol = solutions_dict[moduli[0]]
rows, cols = sample_sol.dimensions()
result = matrix(ZZ, rows, cols)
for i in range(rows):
for j in range(cols):
residues = []
mods = []
for m in moduli:
sol_m = solutions_dict[m]
residues.append(int(sol_m[i,j]))
mods.append(m)
# Use CRT to combine
try:
combined_val = crt(residues, mods)
result[i,j] = combined_val
except Exception as e:
print(f"CRT failed for element ({i},{j}): {e}")
result[i,j] = 0 # Fallback
# Reduce modulo the original modulus
result_mod = result.change_ring(IntegerModRing(modulus))
return result_mod
# Your test code
modulus = 100002
matrix_rows = [
[7, 3, 0, 0, 0, 1, 0, 0, 0, 0, 7781, ],
[8, 1, 0, 1, 1, 0, 0, 0, 0, 0, 3246, ],
[1, 2, 2, 0, 2, 0, 0, 0, 0, 0, 7497, ],
[4, 2, 0, 0, 0, 1, 0, 0, 0, 0, 10491, ],
[1, 0, 0, 1, 2, 0, 0, 0, 0, 0, 86958, ],
[1, 0, 5, 0, 0, 0, 0, 0, 0, 0, 67290, ],
[0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 33792, ],
[10, 0, 1, 0, 0, 0, 0, 1, 0, 0, 94936, ],
[4, 2, 1, 0, 1, 0, 0, 0, 0, 0, 21039, ],
[0, 2, 0, 1, 1, 0, 0, 1, 0, 0, 90010, ],
[0, 0, 3, 1, 0, 0, 0, 0, 0, 0, 61328, ],
[0, 0, 0, 2, 1, 0, 1, 0, 0, 0, 63056, ],
[2, 2, 0, 1, 0, 0, 1, 0, 0, 0, 91681, ],
[2, 0, 0, 1, 0, 0, 1, 0, 0, 1, 13288, ],
[2, 0, 1, 1, 0, 0, 0, 0, 0, 0, 10412, ],
[0, 3, 0, 0, 0, 0, 2, 0, 0, 0, 16161, ],
[1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 76914, ],
[7, 1, 0, 0, 1, 0, 1, 0, 0, 0, 90440, ],
[3, 2, 3, 0, 0, 0, 0, 0, 0, 0, 70954, ],
[1, 0, 2, 2, 0, 0, 0, 0, 0, 0, 823, ],
[0, 3, 0, 0, 1, 0, 0, 0, 0, 1, 49058, ],
[0, 2, 3, 0, 0, 0, 0, 0, 1, 0, 68043, ],
[3, 2, 0, 0, 0, 0, 0, 1, 1, 0, 21144, ],
[0, 2, 0, 2, 0, 0, 0, 0, 0, 1, 60675, ],
[1, 1, 1, 4, 0, 0, 0, 0, 0, 0, 42563, ],
[0, 1, 0, 0, 0, 0, 0, 0, 2, 1, 87658, ],
[1, 1, 2, 2, 0, 1, 0, 0, 0, 0, 74025, ],
[3, 1, 1, 0, 0, 1, 0, 1, 0, 0, 88135, ],
[2, 0, 0, 0, 1, 1, 0, 0, 0, 0, 16913, ],
[6, 2, 1, 0, 0, 1, 0, 0, 0, 0, 55950, ],
[2, 0, 0, 0, 1, 0, 1, 0, 1, 0, 30238, ],
[1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 34517, ],
[3, 2, 1, 0, 1, 0, 0, 1, 0, 0, 10513, ],
[0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 11081, ],
[2, 1, 2, 0, 0, 0, 1, 1, 0, 0, 9826, ],
]
print("Original augmented matrix:")
M_original = matrix(ZZ, matrix_rows)
print(f"Dimensions: {M_original.dimensions()}")
print(f"Modulus: {modulus} = {factor(modulus)}")
# First, let's check the rank over rationals
print(f"\n=== QUICK RANK CHECK ===")
coefficient_matrix = M_original[:, :8]
constants = M_original[:, 8]
print(f"Rank of coefficient matrix (QQ): {coefficient_matrix.rank()}")
print(f"Rank of augmented matrix (QQ): {M_original.rank()}")
# Now compute RREF
final_rref = modular_rref_composite(M_original, modulus)
if final_rref is not None:
print(f"\n=== FINAL RESULT ===")
print(f"RREF modulo {modulus}:")
print(final_rref)
# Analysis
rows, cols = final_rref.dimensions()
variables = cols - 1
print(f"\n=== SOLUTION ANALYSIS ===")
print(f"Matrix dimensions: {rows} x {cols}")
print(f"Variables: {variables}")
# Check consistency
consistent = True
pivot_cols = []
current_row = 0
for i in range(rows):
all_zero = True
for j in range(variables):
if final_rref[i,j] != 0:
all_zero = False
if current_row == i and final_rref[i,j] == 1:
pivot_cols.append(j)
current_row += 1
break
if all_zero and final_rref[i, variables] != 0:
print(f"Inconsistent equation in row {i}: 0 = {final_rref[i, variables]}")
consistent = False
if consistent:
print("System is consistent")
print(f"Pivot columns: {pivot_cols}")
print(f"Free variables: {[j for j in range(variables) if j not in pivot_cols]}")
print(f"Rank: {len(pivot_cols)}")
else:
print("System is inconsistent")
else:
print("Failed to compute RREF")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment