C Function to Identify Conic Sections by Matrix Form

ConicIdentify.c

//Code Explained : https://leetarxiv.substack.com/p/conic-sections-matrix-theory
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>
#include <stdbool.h>
#include <assert.h>
#define COMPARE_EPSILON 1e-6f
#define INDEX(x, y, cols) ((x) * (cols) + (y))
typedef enum matrix_type_enum MatrixType;
enum matrix_type_enum{MATRIX_INT,MATRIX_FLOAT,MATRIX_DOUBLE, MATRIX_TYPE_LENGTH};

//Run: clear && gcc ConicIdentify.c -lm -o m.o && ./m.o
double Determinant2x2(MatrixType type, void *matrix)
{
	switch (type)
	{
		case MATRIX_INT:
		{
			int *m = (int*)matrix;
			return (m[0] * m[3] - m[1] * m[2]);  
		}
		case MATRIX_FLOAT:
		{
			float *m = (float*)matrix;
			return (m[0] * m[3] - m[1] * m[2]);
		}
		case MATRIX_DOUBLE:
		{
			double *m = (double*)matrix;
			return (m[0] * m[3] - m[1] * m[2]);
		}
		default:
		return NAN;  // Not a Number (invalid type)
	}
}
double Determinant3x3(MatrixType type, void *matrix)
{
	switch(type)
	{
		case MATRIX_INT:
		{
			int *m = (int*)matrix;
			return m[0]*(m[4]*m[8] - m[5]*m[7]) - m[1]*(m[3]*m[8] - m[5]*m[6]) + m[2]*(m[3]*m[7] - m[4]*m[6]);
		}
		case MATRIX_FLOAT:
		{
			float *m = (float*)matrix;
			return m[0]*(m[4]*m[8] - m[5]*m[7]) - m[1]*(m[3]*m[8] - m[5]*m[6]) + m[2]*(m[3]*m[7] - m[4]*m[6]);
		}
		case MATRIX_DOUBLE:
		{
			double *m = (double*)matrix;
			return m[0]*(m[4]*m[8] - m[5]*m[7]) - m[1]*(m[3]*m[8] - m[5]*m[6]) + m[2]*(m[3]*m[7] - m[4]*m[6]);
		}
		default:
			return NAN;
	}	
}

bool IsSymmetric(int rows, int cols, MatrixType type, void *matrix)
{
	if(rows != cols)
	{
		return false;
	}
	for(int i = 0; i < rows; i++)
	{
		for(int j = 0; j < cols; j++)
		{
			int index_ij = INDEX(i, j, cols);
			int index_ji = INDEX(j, i, cols);	
			switch(type)
			{
				case MATRIX_INT:
				{
					int value_ij = ((int*)matrix)[index_ij];
					int value_ji = ((int*)matrix)[index_ji];
					if(value_ij != value_ji){return false;}
					break;
				}
				case MATRIX_FLOAT:
				{
					float value_ij = ((float*)matrix)[index_ij];
					float value_ji = ((float*)matrix)[index_ji];
					if(fabs(value_ij - value_ji) > COMPARE_EPSILON){return false;}
					break;
				}
				case MATRIX_DOUBLE:
				{
					double value_ij = ((double*)matrix)[index_ij];
					double value_ji = ((double*)matrix)[index_ji];
					if(fabs(value_ij - value_ji) > COMPARE_EPSILON){return false;}
					break;
				}
				default:
				{
					printf("%d\n", type);
					assert(type < MATRIX_TYPE_LENGTH);
					assert(type > -1);
					return false;
				}
			}
		}	
	}
	return true;
}

void IdentifyConicSection(int rows, int cols, MatrixType type, void *matrix)
{
	//Step 1: Ensure this is a 3*3 matrix
	assert(rows == 3 && cols == 3);
	
	//Step 2: Extract Q = [A B; B C] from the conic matrix:
	//[ A   B   F ]
	//[ B   C   G ]
	//[ F   G   H ]
	void *Q = NULL;int Q_rows = 2, Q_cols = 2;
	switch(type)
	{
		case MATRIX_INT:
		{
			int Q_int[4];
			int *m = (int *)matrix;
			Q_int[0] = m[0]; Q_int[1] = m[1]; // First row: A, B
			Q_int[2] = m[3]; Q_int[3] = m[4]; // Second row: B, C
			Q = Q_int;
			break;
		}
		case MATRIX_FLOAT:
		{
			float Q_float[4];
			float *m = (float *)matrix;
			Q_float[0] = m[0]; Q_float[1] = m[1];
			Q_float[2] = m[3]; Q_float[3] = m[4];
			Q = Q_float;
			break;
		}
		case MATRIX_DOUBLE:
		{
			double Q_double[4];
			double *m = (double *)matrix;
			Q_double[0] = m[0]; Q_double[1] = m[1];
			Q_double[2] = m[3]; Q_double[3] = m[4];
			Q = Q_double;
			break;
		}
		default:
			fprintf(stderr, "Invalid matrix type\n");
			return;
	}

	//Step 2: Ensure Q is symmetric
	assert(true == IsSymmetric(Q_rows, Q_cols, type, Q));
	//Step 3: Ensure M is symmetric
	assert(true == IsSymmetric(rows, cols, type, matrix));
	//Step 4: Find the determinant of Q
	double determinantQ = Determinant2x2(type, Q);
	//Step 5: Find the determinant of M
	double determinantM = Determinant3x3(type, matrix);
	//Print out type
	if(determinantM == 0.0f)
	{
		printf("Degenerate conic: ");
		if(determinantQ == 0)
		{
			printf("Parallel lines\n");		
		}	
		else
		{
			printf("Intersecting lines\n");			
		}
	}
	else
	{
		printf("Non-Degenerate conic: ");	
		if(determinantQ > 0)
		{
			printf("Ellipse\n");		
		}
		else if(determinantQ == 0)
		{
			printf("Parabola\n");		
		}
		else
		{
			printf("Hyperbola\n");		
		}
	}
}

void TestConicSection()
{
	//General form
	//Ax^2 + Bxy + Cy^2 + Fx + Gy + H = 0
	
	/*Ellipse*/
	//1x^2 + 0xy + 1y^2 + 4x + 6y + -23 = 0
	double conic0[] =  
	{
		1.0, 0.0, 4.0,  // A=1.0, B=0.0, F=4.0
		0.0, 1.0, 0.0,  // B=0.0, C=1.0, G=6.0
		4.0, 0.0, -23.0 // F=4.0, G=6.0, H=-23.0
	};
	IdentifyConicSection(3, 3, MATRIX_DOUBLE, conic0);
	
	/*Hyperbola*/
	//3x^2 + -10xy + 3y^2 + 14x + -2y + 3 = 0
	double conic1[] =  
	{
		1.0  , -10.0, 14.0,  // A=3.0,   B=-10.0,F=14.0
		-10.0, 3.0  , -2.0,  // B=-10.0, C=3.0,  G=-2.0
		14.0 , -2.0 , 3.0    // F=14.0,  G=-2.0, H= 3.0
	};
	IdentifyConicSection(3, 3, MATRIX_DOUBLE, conic1);
	
}
int main()
{
        TestConicSection();
        return 0;
}