1 files changed, 99 insertions, 55 deletions

diff --git a/libs/libc/math.c b/libs/libc/math.c
index 606281c..4d31e04 100644
--- a/libs/libc/math.c
+++ b/libs/libc/math.c
@@ -2,104 +2,148 @@
 
 #include <math.h>
 
-f32 powf(f32 base, f32 exp)
+f64 mceil(f64 x)
 {
-	return (f32)pow(base, exp);
+	if (x == 0.0)
+		return x;
+
+	f64 out;
+	__asm__ volatile("frndint\n" : "=t"(out) : "0"(x));
+	if (out < x)
+		return out + 1.0;
+	return out;
 }
 
-f64 pow(f64 base, f64 exp)
+f64 mfloor(f64 x)
 {
+	if (x == 0.0)
+		return x;
+
 	f64 out;
-	__asm__ volatile("fyl2x;"
-			 "fld %%st;"
-			 "frndint;"
-			 "fsub %%st,%%st(1);"
-			 "fxch;"
-			 "fchs;"
-			 "f2xm1;"
-			 "fld1;"
-			 "faddp;"
-			 "fxch;"
-			 "fld1;"
-			 "fscale;"
-			 "fstp %%st(1);"
-			 "fmulp;"
-			 : "=t"(out)
-			 : "0"(base), "u"(exp)
-			 : "st(1)");
+	__asm__ volatile("frndint\n" : "=t"(out) : "0"(x));
+	if (out > x)
+		return out - 1.0;
 	return out;
 }
 
-// TODO: More efficient sqrt?
-
-f32 sqrtf(f32 num)
+f64 mexp(f64 exp)
 {
-	return powf(num, .5);
+	f64 out;
+	__asm__ volatile("fldl2e\n"
+			 "fmulp\n"
+			 "fld1\n"
+			 "fld %%st(1)\n"
+			 "fprem\n"
+			 "f2xm1\n"
+			 "faddp\n"
+			 "fscale\n"
+			 "fstp %%st(1)"
+			 : "=t"(out)
+			 : "0"(exp));
+	return out;
 }
 
-f64 sqrt(f64 num)
+f64 mexp2(f64 exp)
 {
-	return pow(num, .5);
+	f64 out;
+	__asm__ volatile("fld1\n"
+			 "fld %%st(1)\n"
+			 "fprem\n"
+			 "f2xm1\n"
+			 "faddp\n"
+			 "fscale\n"
+			 "fstp %%st(1)"
+			 : "=t"(out)
+			 : "0"(exp));
+	return out;
 }
 
-/**
- * Interpolations
- */
+f64 mlog(f64 x)
+{
+	f64 out;
+	__asm__ volatile("fldln2\n"
+			 "fld %%st(1)\n"
+			 "fyl2x\n"
+			 "fstp %%st(1)"
+			 : "=t"(out)
+			 : "0"(x));
+	return out;
+}
 
-f32 lerpf(f32 from, f32 to, f32 trans)
+f64 mlog2(f64 x)
 {
-	return from + (to - from) * trans;
+	f64 out;
+	__asm__ volatile("fld1\n"
+			 "fld %%st(1)\n"
+			 "fyl2x\n"
+			 "fstp %%st(1)"
+			 : "=t"(out)
+			 : "0"(x));
+	return out;
 }
 
-f64 lerp(f64 from, f64 to, f64 trans)
+f64 mpow(f64 base, f64 exp)
 {
-	return from + (to - from) * trans;
+	if (exp == 0)
+		return 1;
+	if (exp == 1)
+		return base;
+	if (base == 0)
+		return 0;
+	if (exp == (f64)((s32)exp)) {
+		f64 out = base;
+		for (u32 i = 0; i < FABS(exp) - 1; i++)
+			out *= base;
+		if (exp < 0)
+			out = 1.0 / out;
+		return out;
+	}
+	return mexp2(exp * mlog2(base));
 }
 
-f32 blerpf(f32 a, f32 b, f32 c, f32 d, f32 transx, f32 transy)
+// TODO: More efficient sqrt?
+
+f64 msqrt(f64 num)
 {
-	return lerpf(lerpf(a, b, transx), lerpf(c, d, transx), transy);
+	return mpow(num, .5);
 }
 
-f64 blerp(f64 a, f64 b, f64 c, f64 d, f64 transx, f64 transy)
+f64 mcubic(f64 x, f64 a, f64 b, f64 c, f64 d)
 {
-	return lerp(lerp(a, b, transx), lerp(c, d, transx), transy);
+	return a * mpow(x, 3) + b * mpow(x, 2) + c * x + d;
 }
 
 /**
- * Trigonometric functions
+ * Interpolations
  */
 
-f32 sinf(f32 angle)
+f64 mlerp(f64 from, f64 to, f64 trans)
 {
-	f32 ret = 0.0;
-	__asm__ volatile("fsin" : "=t"(ret) : "0"(angle));
-	return ret;
+	return from + (to - from) * trans;
 }
 
-f64 sin(f64 angle)
+f64 mblerp(f64 a, f64 b, f64 c, f64 d, f64 transx, f64 transy)
 {
-	f64 ret = 0.0;
-	__asm__ volatile("fsin" : "=t"(ret) : "0"(angle));
-	return ret;
+	return mlerp(mlerp(a, b, transx), mlerp(c, d, transx), transy);
 }
 
-f32 cosf(f32 angle)
-{
-	return sinf(angle + (f32)M_PI_2);
-}
+/**
+ * Trigonometric functions
+ */
 
-f64 cos(f64 angle)
+f64 msin(f64 angle)
 {
-	return sin(angle + (f64)M_PI_2);
+	f64 ret = 0.0;
+	__asm__ volatile("fsin" : "=t"(ret) : "0"(angle));
+	return ret;
 }
 
-f32 tanf(f32 angle)
+f64 mcos(f64 angle)
 {
-	return (f32)tan(angle);
+	return msin(angle + (f64)M_PI_2);
 }
 
-f64 tan(f64 angle)
+f64 mtan(f64 angle)
 {
 	f64 ret = 0.0, one;
 	__asm__ volatile("fptan" : "=t"(one), "=u"(ret) : "0"(angle));