Unsigned PMat32.

PMat32.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>
#include <stdbool.h>
#include <assert.h>
#define INDEX(x, y, cols) ((x) * (cols) + (y))

//We use unsigned integers
//Using signed integers is misaligned with our goal (Positive only integers)
//Complete guide:https://leetarxiv.substack.com/p/positive-only-binary-pmat32
typedef uint32_t PMat32;

//Run : clear && gcc PMat.c -lm -o m.o && ./m.o
float RandomFloat(float min, float max)
{
	return min + ((float)rand() / RAND_MAX) * (max - min);
}

void FloatRandomMatrix(int rows, int cols, float min, float max, float *matrix)
{
	for(int i = 0; i < rows * cols; i++)
	{
		matrix[i] = RandomFloat(min, max);
	}
}

void FloatMultiply(int m, int n, int p, float *A, float *B, float *C)
{
	for(int i = 0; i < m; i++)
	{
		for(int j = 0; j < p; j++)
		{
			float sum = 0.0f;
			for(int k = 0; k < n; k++)
			{
				sum += A[i * n + k] * B[k * p + j];
			}
			C[i * p + j] = sum;
		}
    	}
}

void FloatMatrixPrint(int rows, int cols, float *matrix)
{
	for(int i = 0; i < rows; i++)
	{
		for(int j = 0; j < cols; j++)
		{
			printf("%6.2f ", matrix[i * cols + j]);
		}
		printf("\n");
	}
	printf("\n");
}

void TestFloatMatrix()
{
	//A * B = C
	int rowA = 3; 
	int colA = 4; 
	int rowB = colA; 
	int colB = 4; 
	
	float *A = malloc(rowA * colA * sizeof(float));
	float *B = malloc(rowB * colB * sizeof(float));
	float *C = malloc(rowA * colB * sizeof(float));
	
	//Generate random data
	FloatRandomMatrix(rowA, colA, -5.0f, 5.0f, A);
	FloatRandomMatrix(rowB, colB, -5.0f, 5.0f, B);
	
	//Multiply matrices
	FloatMultiply(rowA, colA, colB, A, B, C);
	
	//Print Results
	FloatMatrixPrint(rowA, colA, A);
	FloatMatrixPrint(rowB, colB, B);
	FloatMatrixPrint(rowA, colB, C);
	//Free memory
	free(A);free(B);free(C);
}

PMat32 SignedFloatToPMat32(float f, PMat32 ONE_CONSTANT,PMat32 MAX_CONSTANT)
{
	//Convert float to scaled signed integer with rounding
	int signedInteger = (f * ONE_CONSTANT + (f >= 0 ? 0.5f : -0.5f));

	//Convert signed integer to PMat32
	PMat32 result = (signedInteger < 0) ? (MAX_CONSTANT + signedInteger) : signedInteger;
	return result;
}

float PMat32ToSignedFloat(PMat32 value, PMat32 ONE_CONSTANT, PMat32 MAX_CONSTANT)
{
	//1. Convert PMat32 to signed integer in two's complement
	int64_t signedInteger = (int64_t)value;
	if(value >= MAX_CONSTANT / 2)
	{
		//Negative number: wrap around
		signedInteger = (int64_t)value - MAX_CONSTANT;
	}
	
	// 2. Scale back to float 
	// Compiler may optimize to ldexp
	float result = (float)signedInteger / ONE_CONSTANT;
	return result;
}

PMat32 PMat32_Add(PMat32 a, PMat32 b, PMat32 MAX_CONSTANT)
{
	return (a + b) & (MAX_CONSTANT - 1);  // modulo is equivalent to & (MAX_CONSTANT - 1);
}

PMat32 PMat32_Multiply(PMat32 a, PMat32 b, int fractionBits, PMat32 MAX_CONSTANT) 
{
	const PMat32 HALF = MAX_CONSTANT >> 1;
	const PMat32 MASK = MAX_CONSTANT - 1;
	// Branchless Decode to signed integers
	int32_t sa = (int32_t)((a ^ HALF) - HALF);
	int32_t sb = (int32_t)((b ^ HALF) - HALF);

	// Fixed-point multiplication with rounding
	int64_t prod = (int64_t)sa * (int64_t)sb;
	int32_t result = (int32_t)((prod + MAX_CONSTANT) >> fractionBits);

	// Rewrap signed result into [0, MAX_CONSTANT)
	return result & (MAX_CONSTANT - 1);
}

void TestPMat32()
{
	int integerBits = 8;
	int fractionBits= 20;
	assert(integerBits + fractionBits < 32);assert(integerBits > -1);assert(fractionBits > -1);

	/*Create mandatory constants*/
	PMat32 ONE_CONSTANT = (1 << fractionBits);
	PMat32 MAX_CONSTANT = (1 << (integerBits + fractionBits));

	/*Test0 : Test Float to PMat32 Conversion and Reverse*/
	float f0 = 19.34;
	PMat32 test0 = SignedFloatToPMat32(f0, ONE_CONSTANT, MAX_CONSTANT);
	float  test0Reverse = PMat32ToSignedFloat(test0, ONE_CONSTANT, MAX_CONSTANT);
	printf("MAX: %u\nONE: %u\nTest0:(%u) where (%.3f %.3f) \n",MAX_CONSTANT,ONE_CONSTANT,test0,f0, test0Reverse);
	
	/*Test1 : Test PMat32 Addition*/
	float f1 = 23.17;PMat32 test1_f1 = SignedFloatToPMat32(f1, ONE_CONSTANT, MAX_CONSTANT);
	float f2 = -22.17;PMat32 test1_f2 = SignedFloatToPMat32(f2, ONE_CONSTANT, MAX_CONSTANT);
	float f3 = 5.17;PMat32 test1_f3 = SignedFloatToPMat32(f3, ONE_CONSTANT, MAX_CONSTANT);
	PMat32 test1_f1PlusF2 = PMat32_Add(test1_f1, test1_f2,MAX_CONSTANT);
	
	float f4 = f1 + f2 + f3;
	PMat32 test1 = PMat32_Add(test1_f1PlusF2, test1_f3,MAX_CONSTANT);
	float  test1Reverse = PMat32ToSignedFloat(test1, ONE_CONSTANT, MAX_CONSTANT);
	printf("Addition tests: %.3f %.3f\n",f4, test1Reverse);
	
	/*Test2 : Test PMat32 Multiplication*/
	float f5 = 2.5;
	float f6 = 3.15;
	PMat32 test2_f5 = SignedFloatToPMat32(f5, ONE_CONSTANT, MAX_CONSTANT);
	PMat32 test2_f6 = SignedFloatToPMat32(f6, ONE_CONSTANT, MAX_CONSTANT);

	float f7 = f5 * f6;PMat32 test2_f7 = SignedFloatToPMat32(f7, ONE_CONSTANT, MAX_CONSTANT);
	PMat32 test2 = PMat32_Multiply(test2_f5, test2_f6, fractionBits, MAX_CONSTANT);
	float test2Reverse = PMat32ToSignedFloat(test2, ONE_CONSTANT, MAX_CONSTANT);
	printf("Multiply tests: %.3f %.3f\n",f7, test2Reverse);
	
}

int main()
{
	//TestFloatMatrix();
	TestPMat32();
	return 0;
}