One way function with cryptography

find_irreducible_and_primitive_polynomial.py

# Step 1: Install SageMath (if you haven't already)
# Follow the installation instructions at: https://www.sagemath.org/download.html

# Step 2: Import necessary libraries
from sage.all import *

# Step 3: Define the finite field and polynomial ring
n = 129
F = GF(2)  # Define the finite field with 2 elements
R = PolynomialRing(F, 'x')  # Define the polynomial ring over the finite field

# Step 4: Write a function to generate polynomials in reverse order
def generate_polynomials_reverse_order(n):
    for i in reversed(range(2**n)):
        yield R(Integer(i).bits())

# Step 6: Iterate through the polynomials and find the primitive ones
for poly in generate_polynomials_reverse_order(n):
    if poly.is_irreducible() and poly.is_primitive():
        print("Found primitive polynomial:", poly)
        break

####################################################################################################

def is_irreducible_primitive(polynomial):
    R.<x> = PolynomialRing(ZZ)
    sage_poly = R(polynomial)
    
    if not sage_poly.is_irreducible():
        return False

    if not sage_poly.is_primitive():
        return False

    return True

# Calculate the polynomial coefficients for 2^129 - 25
number = 2**129 - 25
binary_representation = bin(number)[2:]  # Convert the number to binary and remove the '0b' prefix
polynomial = [int(bit) for bit in binary_representation]  # Convert each binary digit to an integer and create a list

# Test if the polynomial is irreducible and primitive
result = is_irreducible_primitive(polynomial)
print(f"The polynomial {polynomial} is irreducible and primitive: {result}")

One-way-function.go

/*
package main

import (
    "fmt"
    "math/big"
    "math/rand"
    "time"
)

type Matrix [][]*big.Int

func main() {
    // Define dimensions, probability, and irreducible polynomials
    m, n := 32, 32
    p := 0.01
    pX := new(big.Int)
    //2 ** 129 - 25
    pX.SetString("680564733841876926926749214863536422887", 10) // Replace with your desired irreducible polynomial
    gX := new(big.Int)
    //2 ** 64
    gX.SetString("18446744073709551615", 10) // Replace with your desired irreducible polynomial

    // Initialize matrices A and B with random values
    A := randomMatrix(m, n)
    B := randomMatrix(m, n)

    // Apply the one-way function
    D := oneWayFunction(A, B, p, pX, gX)

    // Print the result
    fmt.Println("Matrix D:")
    printMatrix(D)
}

func randomMatrix(m, n int) Matrix {
    rand.Seed(time.Now().UnixNano())
    M := make(Matrix, m)
    for i := range M {
        M[i] = make([]*big.Int, n)
        for j := range M[i] {
            M[i][j] = big.NewInt(rand.Int63())
        }
    }
    return M
}

func printMatrix(M Matrix) {
    for i := range M {
        for j := range M[i] {
            fmt.Printf("%v\t", M[i][j])
        }
        fmt.Println()
    }
}

func oneWayFunction(A, B Matrix, p float64, pX, gX *big.Int) Matrix {
    APrime := pseudoRandomZeroReplacement(A, p)
    BPrime := pseudoRandomZeroReplacement(B, p)
    APrime = extension(APrime, pX)
    BPrime = extension(BPrime, pX)
    CPrime := hadamardProduct(APrime, BPrime)
    DPrime := moduloReduction(CPrime, gX)
    D := shrinking(DPrime)
    return D
}

func pseudoRandomZeroReplacement(M Matrix, p float64) Matrix {
    rand.Seed(time.Now().UnixNano())
    MPrime := make(Matrix, len(M))
    for i := range M {
        MPrime[i] = make([]*big.Int, len(M[i]))
        for j := range M[i] {
            if rand.Float64() < p {
                MPrime[i][j] = big.NewInt(0)
            } else {
                MPrime[i][j] = new(big.Int).Set(M[i][j])
            }
        }
    }
    return MPrime
}

func extension(M Matrix, pX *big.Int) Matrix {
    MPrime := make(Matrix, len(M))
    for i := range M {
        MPrime[i] = make([]*big.Int, len(M[i]))
        for j := range M[i] {
            MPrime[i][j] = new(big.Int).Mod(M[i][j], pX)
        }
    }
    return MPrime
}

func hadamardProduct(A, B Matrix) Matrix {
    C := make(Matrix, len(A))
    for i := range A {
        C[i] = make([]*big.Int, len(A[i]))
        for j := range A[i] {
            C[i][j] = new(big.Int).Mul(A[i][j], B[i][j])
        }
    }
    return C
}

func moduloReduction(M Matrix, gX *big.Int) Matrix {
    MPrime := make(Matrix, len(M))
    for i := range M {
        MPrime[i] = make([]*big.Int, len(M[i]))
        for j := range M[i] {
            MPrime[i][j] = new(big.Int).Mod(M[i][j], gX)
        }
    }
    return MPrime
}

func shrinking(M Matrix) Matrix {
    MPrime := make(Matrix, len(M))
    twoTo64 := big.NewInt(0).Exp(big.NewInt(2), big.NewInt(64), nil)
    for i := range M {
        MPrime[i] = make([]*big.Int, len(M[i]))
        for j := range M[i] {
            MPrime[i][j] = new(big.Int).Mod(M[i][j], twoTo64)
        }
    }
    return MPrime
}
*/

package main

import (
	cryptoRandom "crypto/rand"
	"fmt"
	"math/big"
	mathRandom "math/rand"
)

type Polynomial struct {
	Coefficients []int
}

func (p *Polynomial) Degree() int {
	return len(p.Coefficients) - 1
}

func (p *Polynomial) Eval(x *big.Int) *big.Int {
	result := big.NewInt(0)
	temp := big.NewInt(1)
	for _, c := range p.Coefficients {
		term := new(big.Int).Mul(temp, big.NewInt(int64(c)))
		result.Add(result, term)
		temp.Mul(temp, x)
	}

	return result
}

type GF struct {
	p *Polynomial
}

func NewGF(p *Polynomial) *GF {
	return &GF{p: p}
}

func (gf *GF) FromInteger(x *big.Int) *big.Int {
	return gf.p.Eval(x)
}

func main() {
	m, n := 64, 64
	probability := 0.01

	// Step 1: Create matrices A and B with given dimensions m and n, and perform Pseudo-Random Zero Replacement
	A := createMatrix(m, n, probability)
	B := createMatrix(m, n, probability)

	// Step 2: Number fields extension
	polynomial := &Polynomial{Coefficients: []int{128, 127, 126, 125, 124, 123, 122, 121, 120, 119, 118, 117, 116, 115, 114, 113, 112, 111, 110, 109, 108, 107, 106, 105, 104, 103, 102, 101, 100, 99, 98, 97, 96, 95, 94, 93, 92, 91, 90, 89, 88, 87, 86, 85, 84, 83, 82, 81, 80, 79, 78, 77, 76, 75, 74, 73, 72, 71, 70, 69, 68, 67, 66, 65, 64, 63, 62, 61, 60, 59, 58, 57, 56, 55, 54, 53, 52, 51, 50, 49, 48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 8, 7, 6, 5, 3, 1}}
	K := NewGF(polynomial)
	A_ext := applyMap(A, K.FromInteger)
	B_ext := applyMap(B, K.FromInteger)

	// Step 3: Hadamard Product
	C_ext := hadamardProduct(A_ext, B_ext)
  
	// Step 4: Modulo Reduction
	mod := new(big.Int).Exp(big.NewInt(2), big.NewInt(64), nil)
	D := applyMap(C_ext, func(x *big.Int) *big.Int { return new(big.Int).Mod(x, mod) })

	// Print Matrix D
	fmt.Println("\nMatrix:")
	printMatrix(D)
}

func createMatrix(m, n int, probability float64) [][]*big.Int {
	matrix := make([][]*big.Int, m)
	for i := 0; i < m; i++ {
		matrix[i] = make([]*big.Int, n)
		for j := 0; j < n; j++ {
			if mathRandom.Float64() > probability {
				matrix[i][j], _ = cryptoRandom.Int(cryptoRandom.Reader, new(big.Int).Exp(big.NewInt(2), big.NewInt(64), nil))
			} else {
				matrix[i][j] = big.NewInt(0)
			}
		}
	}
	return matrix
}

func applyMap(matrix [][]*big.Int, f func(*big.Int) *big.Int) [][]*big.Int {
	m, n := len(matrix), len(matrix[0])
	result := make([][]*big.Int, m)
	for i := 0; i < m; i++ {
		result[i] = make([]*big.Int, n)
		for j := 0; j < n; j++ {
			result[i][j] = f(matrix[i][j])
		}
	}
	return result
}

func hadamardProduct(A, B [][]*big.Int) [][]*big.Int {
	m, n := len(A), len(A[0])
	C := make([][]*big.Int, m)
	for i := 0; i < m; i++ {
		C[i] = make([]*big.Int, n)
		for j := 0; j < n; j++ {
			C[i][j] = new(big.Int).Mul(A[i][j], B[i][j])
		}
	}
	return C
}

func printMatrix(matrix [][]*big.Int) {
	for _, row := range matrix {
		for _, cell := range row {
			fmt.Printf("%d\t", cell)
		}
		fmt.Println()
	}
}

One-way-function.tex

\begin{aligned}
&x = 2 \\
&p(x) = x^{128} + x^{127} + x^{126} + x^{125} + x^{124} + x^{123} + x^{122} + x^{121} + x^{120} + x^{119} + x^{118} + x^{117} + x^{116} + x^{115} + x^{114} + x^{113} + x^{112} + x^{111} + x^{110} + x^{109} + x^{108} + x^{107} + x^{106} + x^{105} + x^{104} + x^{103} + x^{102} + x^{101} + x^{100} \\
&+ x^{99} + x^{98} + x^{97} + x^{96} + x^{95} + x^{94} + x^{93} + x^{92} + x^{91} + x^{90} + x^{89} + x^{88} + x^{87} + x^{86} + x^{85} + x^{84} + x^{83} + x^{82} + x^{81} + x^{80} \\
&+ x^{79} + x^{78} + x^{77} + x^{76} + x^{75} + x^{74} + x^{73} + x^{72} + x^{71} + x^{70} + x^{69} + x^{68} + x^{67} + x^{66} + x^{65} + x^{64} + x^{63} + x^{62} + x^{61} + x^{60} \\
&+ x^{59} + x^{58} + x^{57} + x^{56} + x^{55} + x^{54} + x^{53} + x^{52} + x^{51} + x^{50} + x^{49} + x^{48} + x^{47} + x^{46} + x^{45} + x^{44} + x^{43} + x^{42} + x^{41} + x^{40} \\
&+ x^{39} + x^{38} + x^{37} + x^{36} + x^{35} + x^{34} + x^{33} + x^{32} + x^{31} + x^{30} + x^{29} + x^{28} + x^{27} + x^{26} + x^{25} + x^{24} + x^{23} + x^{22} + x^{21} + x^{20} \\
&+ x^{19} + x^{18} + x^{17} + x^{16} + x^{15} + x^{14} + x^{13} + x^{12} + x^{11} + x^{10} + x^8 + x^7 + x^6 + x^5 + x^3 + 1
\end{aligned}

\newline

\begin{aligned}
&\mathbf{probability = 0.01}\\
&\text{Given two matrices } A \text{ and } B \text{ of dimensions } m \times n, \text{ with elements } a_{ij}, b_{ij} \in \mathbb{F}_{2^{64}}, \\
&\text{we can construct a one-way function using the following steps:} \\
\\
    &\textbf{Step 1: Pseudo-Random Zero Replacement.} \\
    &\text{Replace certain elements of } A \text{ and } B \text{ pseudo-randomly with 0 according to a probability } p: \\
    &A_{ij} = \left\{ \begin{array}{ll}
        0 & \text{with probability } p \\
        a_{ij} & \text{otherwise}
        \end{array} \right., \\
    &B_{ij} = \left\{ \begin{array}{ll}
        0 & \text{with probability } p \\
        b_{ij} & \text{otherwise}
        \end{array} \right. \\
\\
    &\textbf{Step 2: Number fields extension.} \\
    &\text{Extend the elements of } A'' \text{ and } B'' \text{ to a larger polynomial field using the given (binary) irreducible primitive polynomial } p(x): \\
    &A'_{ij} = \text{CopyElement}(A_{ij}, p(x)), \quad B'_{ij} = \text{CopyElement}(B_{ij}, p(x)). \\
\\
    &\textbf{Step 3: Hadamard Product.} \\
    &\text{Compute the Hadamard product } C' \text{ of the extended matrices } A' \text{ and } B': \\
    &C_{ij} = A'_{ij} \odot B'_{ij}.\\
    &C_{ij} \in \mathbb{F}_{2^{129}} \\
\\
    &\textbf{Step 4: Modulo Reduction.} \\
    &\text{Obtain a temporary result } D' \text{ by taking the modulo of each element of } C \text{ with } 2^{64}: \\
    &D_{ij} = C_{ij} \mod 2^{64}, \\
\\
\end{aligned}

Onw-way-function.cpp

#if 1

#include <iostream>
#include <vector>
#include <random>
#include <NTL/GF2EX.h>
#include <NTL/ZZ_pEX.h>
#include <NTL/ZZ_p.h>
#include <NTL/matrix.h>
#include <NTL/mat_ZZ_p.h>

/*
	https://libntl.org/doc/tour-struct.html
	ZZ: big integers
	ZZ_p: big integers modulo p
	zz_p: integers mod "single precision" p
	GF2: integers mod 2
	ZZX: univariate polynomials over ZZ
	ZZ_pX: univariate polynomials over ZZ_p
	zz_pX: univariate polynomials over zz_p
	GF2X: polynomials over GF2
	ZZ_pE: ring/field extension over ZZ_p
	zz_pE: ring/field extension over zz_p
	GF2E: ring/field extension over GF2
	ZZ_pEX: univariate polynomials over ZZ_pE
	zz_pEX: univariate polynomials over zz_pE
	GF2EX: univariate polynomials over GF2E
*/

NTL::Mat<NTL::ZZ> CreateMatrixWithRandomNumber(std::uint64_t m, std::uint64_t n, double probability) {
	std::random_device rd;
	std::mt19937 gen(rd());
	std::uniform_real_distribution<double> dist(0.0, 1.0);
	std::uniform_int_distribution<uint64_t> dist_int(0, std::numeric_limits<uint64_t>::max());

	NTL::Mat<NTL::ZZ> matrix;
	matrix.SetDims(m, n);
	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			//1.Generates pseudorandom number probabilities with normalization to a uniform distribution (0~1).
			//2.According to the generated regularized probability, if this probability is greater than the set probability, then the matrix is assigned elements that are not zero, otherwise this probability is less than the set probability, then the matrix is assigned elements from a pseudo-randomly generated function or a uniformly distributed randomly generated function, which is generated from an infinite domain of discrete elements of the integers, regularized to a finite domain GF(2^64)
			matrix[i][j] = dist(gen) > probability ? NTL::ZZ(dist_int(gen)) : NTL::ZZ(0);
		}
	}

	return matrix;
}

NTL::Mat<NTL::ZZ_p> CreateMatrixWithZeroAndModulo(std::uint64_t m, std::uint64_t n, const NTL::ZZ& modulo_number) {
	NTL::Mat<NTL::ZZ_p> matrix;
	matrix.SetDims(m, n);
	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			matrix[i][j] = 0;
			matrix[i][j].init(modulo_number);
		}
	}

	return matrix;
}

NTL::ZZ GF2XToZZ(const NTL::GF2X& poly) {
	NTL::ZZ result(0);
	long degree = NTL::deg(poly);
	for (long i = degree; i >= 0; --i) {
		// Get the coefficient at position i and convert it to an integer (0 or 1)
		int coef = NTL::rep(NTL::coeff(poly, i));
		// Shift the result left by 1 bit and add the coefficient
		result = (result << 1) | coef;
	}
	if(degree > 64)
		std::cout << "Note: that this is a number exceeding 64-bits from matrix value, please be cautious." << std::endl;
	return result;
}

NTL::GF2X ZZToGF2X(const NTL::ZZ& value) {
	NTL::GF2X poly;
	long num_bits = NTL::NumBits(value);
	for (long i = 0; i < num_bits; ++i) {
		// Get the bit at position i
		int bit_value = NTL::bit(value, i);
		// Set the coefficient at position i
		NTL::SetCoeff(poly, i, bit_value);
	}
	if(num_bits > 64)
		std::cout << "Note: that this is a number exceeding 64bits from matrix value, please be cautious." << std::endl;
	return poly;
}

NTL::GF2EX ZZToGF2EX(const NTL::ZZ& value) {
	NTL::GF2EX poly;
	long number_bits = NTL::NumBits(value);
	for (long i = 0; i < number_bits; ++i) {
		// Get the bit at position i
		int bit_value = NTL::bit(value, i);
		// Set the coefficient at position i
		NTL::SetCoeff(poly, i, bit_value);
	}
	return poly;
}

NTL::ZZ GF2EXToZZ(const NTL::GF2EX& poly) {
	NTL::ZZ result(0);
	long degree = NTL::deg(poly);
	for (long i = degree; i >= 0; --i) {
		// Get the coefficient at position i and convert it to a GF2X object
		NTL::GF2E coef = NTL::coeff(poly, i);
		NTL::GF2X coef_value = NTL::rep(coef);
		// Convert the GF2X object to an integer (0 or 1)
		int bit_value = NTL::IsOne(coef_value);
		// Shift the result left by 1 bit and add the coefficient
		result = (result << 1) | bit_value;
	}
	return result;
}

NTL::ZZ_p GF2X_to_ZZ_p(const NTL::GF2X& gf2x) {
	NTL::ZZ number(0);
	long degree = NTL::deg(gf2x);
	for (long i = degree; i >= 0; --i) {
		// Get the coefficient at position i and convert it to an integer (0 or 1)
		int coeff = NTL::rep(NTL::coeff(gf2x, i));
		// Shift the result left by 1 bit and add the coefficient
		number = (number << 1) | coeff;
	}
	return NTL::to_ZZ_p(number);
}

NTL::GF2X ZZ_p_to_GF2X(const NTL::ZZ_p& zz_p) {
	NTL::GF2X gf2x;
	NTL::ZZ num = NTL::rep(zz_p);
	long number_bits = NTL::NumBits(num);
	for (long i = 0; i < number_bits; ++i) {
		if (NTL::bit(num, i)) {
			NTL::SetCoeff(gf2x, i);
		}
	}
	return gf2x;
}

int main() {
	//Number 680564733841876926926749214863536422887 is Formatted irreducible and primitive polynomial in GF(2)
	//x = 2, p(x) = x^128 + x^127 + x^126 + x^125 + x^124 + x^123 + x^122 + x^121 + x^120 + x^119 + x^118 + x^117 + x^116 + x^115 + x^114 + x^113 + x^112 + x^111 + x^110 + x^109 + x^108 + x^107 + x^106 + x^105 + x^104 + x^103 + x^102 + x^101 + x^100 + x^99 + x^98 + x^97 + x^96 + x^95 + x^94 + x^93 + x^92 + x^91 + x^90 + x^89 + x^88 + x^87 + x^86 + x^85 + x^84 + x^83 + x^82 + x^81 + x^80 + x^79 + x^78 + x^77 + x^76 + x^75 + x^74 + x^73 + x^72 + x^71 + x^70 + x^69 + x^68 + x^67 + x^66 + x^65 + x^64 + x^63 + x^62 + x^61 + x^60 + x^59 + x^58 + x^57 + x^56 + x^55 + x^54 + x^53 + x^52 + x^51 + x^50 + x^49 + x^48 + x^47 + x^46 + x^45 + x^44 + x^43 + x^42 + x^41 + x^40 + x^39 + x^38 + x^37 + x^36 + x^35 + x^34 + x^33 + x^32 + x^31 + x^30 + x^29 + x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^8 + x^7 + x^6 + x^5 + x^3 + 1

	//This operation affects all ZZ_p classes
	NTL::ZZ_p::init(NTL::conv<NTL::ZZ>("680564733841876926926749214863536422887"));

	std::uint64_t m = 64, n = 64;
	double probability = 0.01;

	// Step 1: Create matrices A and B with given dimensions m and n, and perform Pseudo-Random Zero Replacement
	NTL::Mat<NTL::ZZ> A = CreateMatrixWithRandomNumber(m, n, probability);
	NTL::Mat<NTL::ZZ> B = CreateMatrixWithRandomNumber(m, n, probability);

	// Step 2: Number fields extension
	NTL::Mat<NTL::ZZ_p> A_ext = NTL::Mat<NTL::ZZ_p>(NTL::INIT_SIZE, n, m);
	NTL::Mat<NTL::ZZ_p>	B_ext = NTL::Mat<NTL::ZZ_p>(NTL::INIT_SIZE, n, m);
	// Convert elements of A and B to GF(p(x))
	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			NTL::ZZ_p ZZp_A(0);
			NTL::ZZ_p ZZp_B(0);
			
			ZZp_A = NTL::to_ZZ_p(A[i][j]);
			ZZp_B = NTL::to_ZZ_p(B[i][j]);

			A_ext[i][j] = ZZp_A;
			B_ext[i][j] = ZZp_B;

			//std::cout << "A_ext: " << A_ext[i][j] << std::endl;
			//std::cout << "B_ext: " << B_ext[i][j] << std::endl;
		}
	}

	// Step 3: Hadamard Product
	NTL::Mat<NTL::ZZ_p> C_ext = NTL::Mat<NTL::ZZ_p>(NTL::INIT_SIZE, n, m);
	std::cout << "\nBig Number Matrix:" << std::endl;
	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			NTL::ZZ_p ZZp_C(0);
			ZZp_C = A_ext[i][j] * B_ext[i][j];

			C_ext[i][j] = ZZp_C;
			std::cout << "C_ext: " << C_ext[i][j] << std::endl;
		}
	}

	// Step 4: Modulo Reduction
	NTL::Mat<NTL::ZZ_p> D;
	D.SetDims(m, n);
	NTL::GF2X big_number_poly;
	NTL::GF2X modulo_poly;
	for (int i = 0; i < 64; ++i) {
		//18446744073709551616 in GF(2^64)
		NTL::SetCoeff(modulo_poly, i);
	}
	std::cout << "Modulo Polynomial: " << modulo_poly << std::endl;

	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			big_number_poly = ZZ_p_to_GF2X(C_ext[i][j]);
			D[i][j] = GF2X_to_ZZ_p(big_number_poly % modulo_poly);
		}
	}

	// Print Matrix
	std::uint64_t BitIndex = 0;
	std::cout << "\nNumber Matrix: " << std::endl;
	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			std::cout << D[i][j] << ",";
		}
		std::cout << '\n';
	}

	return 0;
}

#else

#include <iostream>
#include <vector>
#include <random>
#include <NTL/GF2EX.h>
#include <NTL/ZZ_pEX.h>
#include <NTL/ZZ_p.h>
#include <NTL/matrix.h>
#include <NTL/mat_ZZ_p.h>

/*
	https://libntl.org/doc/tour-struct.html
	ZZ: big integers
	ZZ_p: big integers modulo p
	zz_p: integers mod "single precision" p
	GF2: integers mod 2
	ZZX: univariate polynomials over ZZ
	ZZ_pX: univariate polynomials over ZZ_p
	zz_pX: univariate polynomials over zz_p
	GF2X: polynomials over GF2
	ZZ_pE: ring/field extension over ZZ_p
	zz_pE: ring/field extension over zz_p
	GF2E: ring/field extension over GF2
	ZZ_pEX: univariate polynomials over ZZ_pE
	zz_pEX: univariate polynomials over zz_pE
	GF2EX: univariate polynomials over GF2E
*/

NTL::Mat<NTL::ZZ> CreateMatrixWithRandomNumber(std::uint64_t m, std::uint64_t n, double probability) {
	std::random_device rd;
	std::mt19937 gen(rd());
	std::uniform_real_distribution<double> dist(0.0, 1.0);
	std::uniform_int_distribution<uint64_t> dist_int(0, std::numeric_limits<uint64_t>::max());

	NTL::Mat<NTL::ZZ> matrix;
	matrix.SetDims(m, n);
	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			matrix[i][j] = dist(gen) > probability ? NTL::ZZ(dist_int(gen)) : NTL::ZZ(0);
		}
	}

	return matrix;
}

NTL::Mat<NTL::ZZ_p> CreateMatrixWithZeroAndModulo(std::uint64_t m, std::uint64_t n, const NTL::ZZ& modulo_number) {
	NTL::Mat<NTL::ZZ_p> matrix;
	matrix.SetDims(m, n);
	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			matrix[i][j] = 0;
			matrix[i][j].init(modulo_number);
		}
	}

	return matrix;
}

NTL::ZZ GF2XToZZ(const NTL::GF2X& poly) {
	NTL::ZZ result(0);
	long degree = NTL::deg(poly);
	for (long i = degree; i >= 0; --i) {
		// Get the coefficient at position i and convert it to an integer (0 or 1)
		int coef = NTL::rep(NTL::coeff(poly, i));
		// Shift the result left by 1 bit and add the coefficient
		result = (result << 1) | coef;
	}
	if(degree > 64)
		std::cout << "Note: that this is a number exceeding 64-bits from matrix value, please be cautious." << std::endl;
	return result;
}

NTL::GF2X ZZToGF2X(const NTL::ZZ& value) {
	NTL::GF2X poly;
	long num_bits = NTL::NumBits(value);
	for (long i = 0; i < num_bits; ++i) {
		// Get the bit at position i
		int bit_value = NTL::bit(value, i);
		// Set the coefficient at position i
		NTL::SetCoeff(poly, i, bit_value);
	}
	if(num_bits > 64)
		std::cout << "Note: that this is a number exceeding 64bits from matrix value, please be cautious." << std::endl;
	return poly;
}

NTL::GF2EX ZZToGF2EX(const NTL::ZZ& value) {
	NTL::GF2EX poly;
	long num_bits = NTL::NumBits(value);
	for (long i = 0; i < num_bits; ++i) {
		// Get the bit at position i
		int bit_value = NTL::bit(value, i);
		// Set the coefficient at position i
		NTL::SetCoeff(poly, i, bit_value);
	}
	return poly;
}

NTL::ZZ GF2EXToZZ(const NTL::GF2EX& poly) {
	NTL::ZZ result(0);
	long degree = NTL::deg(poly);
	for (long i = degree; i >= 0; --i) {
		// Get the coefficient at position i and convert it to a GF2X object
		NTL::GF2E coef = NTL::coeff(poly, i);
		NTL::GF2X coef_value = NTL::rep(coef);
		// Convert the GF2X object to an integer (0 or 1)
		int bit_value = NTL::IsOne(coef_value);
		// Shift the result left by 1 bit and add the coefficient
		result = (result << 1) | bit_value;
	}
	return result;
}

int main() {
	std::uint64_t m = 64, n = 64;
	double probability = 0.01;

	// Step 1: Create matrices A and B with given dimensions m and n, and perform Pseudo-Random Zero Replacement
	NTL::Mat<NTL::ZZ> A = CreateMatrixWithRandomNumber(m, n, probability);
	NTL::Mat<NTL::ZZ> B = CreateMatrixWithRandomNumber(m, n, probability);

	// Step 2: Number fields extension
	NTL::GF2X primitive_poly;
	NTL::SetCoeff(primitive_poly, 128);
	for (int i = 1; i <= 127; ++i) {
		if (i != 9 && i != 4 && i != 2 && i != 1) {
			NTL::SetCoeff(primitive_poly, i);
		}
	}
	NTL::SetCoeff(primitive_poly, 0);
	NTL::GF2E::init(primitive_poly);

	NTL::Mat<NTL::GF2X> A_ext, B_ext;
	// Resize A_ext and B_ext to have the same dimensions as A and B
	A_ext.SetDims(m, n);
	B_ext.SetDims(m, n);

	// Convert elements of A and B to GF2E
	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			A_ext[i][j] = ZZToGF2X(A[i][j]);
			B_ext[i][j] = ZZToGF2X(B[i][j]);

			//std::cout << "A_ext" << A_ext[i][j] << std::endl;
			//std::cout << "B_ext" << B_ext[i][j] << std::endl;
		}
	}

	// Step 3: Hadamard Product
	NTL::Mat<NTL::GF2X> C_ext;
	C_ext.SetDims(m, n);
	std::cout << "\nBig Number Matrix:" << std::endl;
	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			C_ext[i][j] = A_ext[i][j] * B_ext[i][j];
			std::cout << "C_ext" << C_ext[i][j] << std::endl;
		}
	}

	// Step 4: Modulo Reduction
	NTL::Mat<NTL::GF2X> D;
	D.SetDims(m, n);

	NTL::GF2X modulo_poly;
	for (int i = 0; i < 64; ++i) {
		NTL::SetCoeff(modulo_poly, i);
	}
	std::cout << "Modulo Polynomial :" << modulo_poly << std::endl;
	
	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			D[i][j] = C_ext[i][j] % modulo_poly;
		}
	}

	// Print Matrix
	std::uint64_t BitIndex = 0;
	std::cout << "\nNumber Matrix:" << std::endl;
	for (int i = 0; i < m; ++i) {
		for (int j = 0; j < n; ++j) {
			NTL::ZZ big_number_value = GF2XToZZ(D[i][j]);
			std::cout << big_number_value << ",";
		}
		std::cout << '\n';
	}

	return 0;
}

#endif

Onw-way-function.java

import java.math.BigInteger;
import java.util.Random;

public class BigNumberMatrix {
    public static BigInteger[][] createMatrixWithRandomNumber(int m, int n, double probability) {
        BigInteger[][] matrix = new BigInteger[m][n];
        Random random = new Random();

        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                if (random.nextDouble() > probability) {
                    matrix[i][j] = new BigInteger(64, random);
                } else {
                    matrix[i][j] = BigInteger.ZERO;
                }
            }
        }

        return matrix;
    }

    public static void main(String[] args) {
        int m = 128;
        int n = 128;
        double probability = 0.01;

        BigInteger[][] A = createMatrixWithRandomNumber(m, n, probability);
        BigInteger[][] B = createMatrixWithRandomNumber(m, n, probability);

        BigInteger[][] C = new BigInteger[m][n];
        BigInteger bigModuloNumber = new BigInteger("680564733841876926926749214863536422887");
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                C[i][j] = A[i][j].multiply(B[i][j]);
                C[i][j] = C[i][j].mod(bigModuloNumber);
            }
        }

        BigInteger moduloNumber = new BigInteger("18446744073709551616");
        BigInteger[][] D = new BigInteger[m][n];

        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                D[i][j] = C[i][j].mod(moduloNumber);
            }
        }

        System.out.println("\nNumber Matrix: ");
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                System.out.print(D[i][j] + ",");
            }
            System.out.println();
        }
    }
}

Onw-way-function.py

```
# Step 1: Install SageMath (if you haven't already)
# Follow the installation instructions at: https://www.sagemath.org/download.html

# Step 2: Import necessary libraries
from sage.all import *

# Step 3: Define the finite field and polynomial ring
n = 129
F = GF(2)  # Define the finite field with 2 elements
R = PolynomialRing(F, 'x')  # Define the polynomial ring over the finite field

# Step 4: Write a function to generate polynomials in reverse order
def generate_polynomials_reverse_order(n):
    for i in reversed(range(2**n)):
        yield R(Integer(i).bits())

# Step 6: Iterate through the polynomials and find the primitive ones
for poly in generate_polynomials_reverse_order(n):
    if poly.is_irreducible() and poly.is_primitive():
        print("Found primitive polynomial:", poly)
        break
```
#Found primitive polynomial: x^128 + x^127 + x^126 + x^125 + x^124 + x^123 + x^122 + x^121 + x^120 + x^119 + x^118 + x^117 + x^116 + x^115 + x^114 + x^113 + x^112 + x^111 + x^110 + x^109 + x^108 + x^107 + x^106 + x^105 + x^104 + x^103 + x^102 + x^101 + x^100 + x^99 + x^98 + x^97 + x^96 + x^95 + x^94 + x^93 + x^92 + x^91 + x^90 + x^89 + x^88 + x^87 + x^86 + x^85 + x^84 + x^83 + x^82 + x^81 + x^80 + x^79 + x^78 + x^77 + x^76 + x^75 + x^74 + x^73 + x^72 + x^71 + x^70 + x^69 + x^68 + x^67 + x^66 + x^65 + x^64 + x^63 + x^62 + x^61 + x^60 + x^59 + x^58 + x^57 + x^56 + x^55 + x^54 + x^53 + x^52 + x^51 + x^50 + x^49 + x^48 + x^47 + x^46 + x^45 + x^44 + x^43 + x^42 + x^41 + x^40 + x^39 + x^38 + x^37 + x^36 + x^35 + x^34 + x^33 + x^32 + x^31 + x^30 + x^29 + x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^8 + x^7 + x^6 + x^5 + x^3 + 1


from sage.all import *
import random

m = 64
n = 64
probability = 0.01

# Step 1: Create matrices A and B with given dimensions m and n, and perform Pseudo-Random Zero Replacement
def random_matrix(m, n):
    M = Matrix(ZZ, m, n)
    for i in range(m):
        for j in range(n):
            M[i, j] = random.randint(0, 2**64 - 1)
    return M

def pseudo_random_zero_replacement(M, p):
    M_prime = deepcopy(M)
    for i in range(M.nrows()):
        for j in range(M.ncols()):
            if random.random() < p:
                M_prime[i, j] = 0
    return M_prime

def polynomial_to_number(polynomial):
    bits = [int(coef) for coef in polynomial.list()]
    bits_reversed = bits[::-1]  # Reverse the list since the polynomial is represented in reverse
    binary_str = ''.join(map(str, bits_reversed))
    number = int(binary_str, 2)
    return number

A = random_matrix(m, n)
B = random_matrix(m, n)
A_prime = pseudo_random_zero_replacement(A, probability)
B_prime = pseudo_random_zero_replacement(B, probability)

# Step 2: Number fields extension
R = PolynomialRing(GF(2), 'x')
primitive_poly = R('x^128 + x^127 + x^126 + x^125 + x^124 + x^123 + x^122 + x^121 + x^120 + x^119 + x^118 + x^117 + x^116 + x^115 + x^114 + x^113 + x^112 + x^111 + x^110 + x^109 + x^108 + x^107 + x^106 + x^105 + x^104 + x^103 + x^102 + x^101 + x^100 + x^99 + x^98 + x^97 + x^96 + x^95 + x^94 + x^93 + x^92 + x^91 + x^90 + x^89 + x^88 + x^87 + x^86 + x^85 + x^84 + x^83 + x^82 + x^81 + x^80 + x^79 + x^78 + x^77 + x^76 + x^75 + x^74 + x^73 + x^72 + x^71 + x^70 + x^69 + x^68 + x^67 + x^66 + x^65 + x^64 + x^63 + x^62 + x^61 + x^60 + x^59 + x^58 + x^57 + x^56 + x^55 + x^54 + x^53 + x^52 + x^51 + x^50 + x^49 + x^48 + x^47 + x^46 + x^45 + x^44 + x^43 + x^42 + x^41 + x^40 + x^39 + x^38 + x^37 + x^36 + x^35 + x^34 + x^33 + x^32 + x^31 + x^30 + x^29 + x^28 + x^27 + x^26 + x^25 + x^24 + x^23 + x^22 + x^21 + x^20 + x^19 + x^18 + x^17 + x^16 + x^15 + x^14 + x^13 + x^12 + x^11 + x^10 + x^8 + x^7 + x^6 + x^5 + x^3 + 1')
print("primitive_poly\n", primitive_poly)
primitive_poly_number = polynomial_to_number(primitive_poly)
print("primitive_poly_number\n", primitive_poly_number)
K.<x> = GF(2**128, name='x', impl='ntl', modulus=primitive_poly)
A_ext = A_prime.apply_map(lambda a_ij: K.from_integer(a_ij))
B_ext = B_prime.apply_map(lambda b_ij: K.from_integer(b_ij))

# Step 3: Hadamard Product
C_ext = matrix(K, m, n, lambda i, j: A_ext[i, j] * B_ext[i, j])

# Step 4: Modulo Reduction
D = C_ext.apply_map(lambda c_ij: Integer(c_ij.to_integer()) % 2**64)

# Print Matrix D
print("Matrix:")
print(D)