From 30614c758c35229b412d0daccd906ef7f087b95a Mon Sep 17 00:00:00 2001
From: Anonimus <anonimus.srb2k@gmail.com>
Date: Wed, 1 Oct 2025 10:09:31 -0400
Subject: [PATCH] 64-bit FixedMul and FixedSqrt, integer sqrt functions

---
 src/m_fixed.c | 44 ++++++++++++++++++++++++++++++++++++++++++++
 src/m_fixed.h | 23 +++++++++++++++++++++++
 2 files changed, 67 insertions(+)

diff --git a/src/m_fixed.c b/src/m_fixed.c
index 5e5e589fb..c5b4ab03f 100644
--- a/src/m_fixed.c
+++ b/src/m_fixed.c
@@ -61,6 +61,45 @@ fixed_t FixedSqrt(fixed_t x)
 #endif
 }
 
+// Shitty 64-bit duplicate so certain distancing edgecases in k_odds die
+INT64 FixedSqrt64(INT64 x)
+{
+#ifdef HAVE_SQRT
+	const float fx = FixedToFloat64(x);
+	float fr;
+#ifdef HAVE_SQRTF
+	fr = sqrtf(fx);
+#else
+	fr = (float)sqrt(fx);
+#endif
+	return FloatToFixed64(fr);
+#else
+	// The neglected art of Fixed Point arithmetic
+	// Jetro Lauha
+	// Seminar Presentation
+	// Assembly 2006, 3rd- 6th August 2006
+	// (Revised: September 13, 2006)
+	// URL: http://jet.ro/files/The_neglected_art_of_Fixed_Point_arithmetic_20060913.pdf
+	register UINT64 root, remHi, remLo, testDiv, count;
+	root = 0;         /* Clear root */
+	remHi = 0;        /* Clear high part of partial remainder */
+	remLo = x;        /* Get argument into low part of partial remainder */
+	count = (31 + (FRACBITS >> 1));    /* Load loop counter */
+	do
+	{
+		remHi = (remHi << 2) | (remLo >> 62); remLo <<= 2;  /* get 2 bits of arg */
+		root <<= 1;   /* Get ready for the next bit in the root */
+		testDiv = (root << 1) + 1;    /* Test radical */
+		if (remHi >= testDiv)
+		{
+			remHi -= testDiv;
+			root += 1;
+		}
+	} while (count-- != 0);
+	return root;
+#endif
+}
+
 fixed_t FixedHypot(fixed_t x, fixed_t y)
 {
 	// Moved the code from R_PointToDist2 to here,
@@ -1122,6 +1161,11 @@ FUNCMATH FUNCINLINE static inline fixed_t FixedMulC(fixed_t a, fixed_t b)
 	return (fixed_t)((((INT64)a * b)) / FRACUNIT);
 }
 
+FUNCMATH FUNCINLINE static inline INT64 FixedMul64C(fixed_t a, fixed_t b)
+{
+	return (INt64)((((INT64)a * b)) / FRACUNIT);
+}
+
 FUNCMATH FUNCINLINE static inline fixed_t FixedDivC2(fixed_t a, fixed_t b)
 {
 	INT64 ret;
diff --git a/src/m_fixed.h b/src/m_fixed.h
index b7722bd15..2080c396f 100644
--- a/src/m_fixed.h
+++ b/src/m_fixed.h
@@ -55,11 +55,21 @@ FUNCMATH FUNCINLINE static ATTRINLINE float FixedToFloat(fixed_t x)
 	return x / (float)FRACUNIT;
 }
 
+FUNCMATH FUNCINLINE static ATTRINLINE float FixedToFloat64(INT64 x)
+{
+	return (float)x / (float)FRACUNIT;
+}
+
 FUNCMATH FUNCINLINE static ATTRINLINE fixed_t FloatToFixed(float f)
 {
 	return (fixed_t)(f * FRACUNIT);
 }
 
+FUNCMATH FUNCINLINE static ATTRINLINE INT64 FloatToFixed64(float f)
+{
+	return (INT64)(f * FRACUNIT);
+}
+
 // for backwards compat
 #define FIXED_TO_FLOAT(x) FixedToFloat(x) // (((float)(x)) / ((float)FRACUNIT))
 #define FLOAT_TO_FIXED(f) FloatToFixed(f) // (fixed_t)((f) * ((float)FRACUNIT))
@@ -79,6 +89,13 @@ FUNCMATH FUNCINLINE static ATTRINLINE fixed_t FixedMul(fixed_t a, fixed_t b)
 	return (fixed_t)(((UINT64)((INT64)a * b)) >> FRACBITS);
 }
 
+FUNCMATH FUNCINLINE static ATTRINLINE INT64 FixedMul64(INT64 a, INT64 b)
+{
+	// Need to cast to unsigned before shifting to avoid undefined behaviour
+	// for negative integers
+	return (INT64)(((UINT64)(a * b)) >> FRACBITS);
+}
+
 /**	\brief	The FixedDiv2 function
 
 	\param	a	fixed_t number
@@ -145,6 +162,12 @@ FUNCMATH FUNCINLINE static ATTRINLINE fixed_t FixedPercentageFloat(float fpercen
 
 */
 FUNCMATH fixed_t FixedSqrt(fixed_t x);
+FUNCMATH INT64 FixedSqrt64(INT64 x);
+
+// Theoretically, we can feed this any integer value and it'll output an
+// accurate square root * 256.
+#define IntSqrt(x) (INT32)(FixedSqrt(x) / 256)
+#define IntSqrt64(x) (INT64)(FixedSqrt64(x) / 256)
 
 /**	\brief	The FixedHypot function