sys_math, cosf OK and documented, sinf, coss, sins documented (#289)

* OK, data imported, document

* spec

* Match cosf, document sinf and cosf

* Some more documentation, switch fu back

* Format

* Little more documentation

* Document coss too

* Move macros

* Review, add some .s

* Fix headers, review

include/convert.h

@@ -1,6 +1,8 @@
 #ifndef _ULTRA64_CONVERT_H_
 #define _ULTRA64_CONVERT_H_
 
+#include "stdint.h"
+
 #define OS_CLOCK_RATE           62500000LL
 #define OS_CPU_COUNTER          (OS_CLOCK_RATE*3/4)
 

include/functions.h

@@ -3269,9 +3269,9 @@ void Check_ClearRGBA16(s16* buffer);
 // void Check_DrawRegionLockErrorMessage(void);
 // void Check_ExpansionPak(void);
 // void Check_RegionIsSupported(void);
-f32 func_80179300(f32 param_1);
-f32 func_80179400(s32 param_1);
-f32 pow_int(f32 x, s32 pow);
+f32 func_80179300(f32 n);
+f32 func_80179400(s32 n);
+f32 pow_int(f32 base, s32 exp);
 f32 sin_rad(f32 rad);
 f32 cos_rad(f32 rad);
 f32 Rand_ZeroFloat(f32 scale);

include/macros.h

@@ -1,8 +1,9 @@
 #ifndef _MACROS_H_
 #define _MACROS_H_
 
-#include "convert.h"
 #include "stdint.h"
+#include "convert.h"
+#include "z64.h"
 
 #define ARRAY_COUNT(arr) (s32)(sizeof(arr) / sizeof(arr[0]))
 #define ARRAY_COUNTU(arr) (u32)(sizeof(arr) / sizeof(arr[0]))
@@ -116,6 +117,8 @@ extern GraphicsContext* __gfxCtx;
 #define CLAMP_MAX(x, max) ((x) > (max) ? (max) : (x))
 #define CLAMP_MIN(x, min) ((x) < (min) ? (min) : (x))
 
+#define ROUND(x) (s32)(((x) >= 0.0) ? ((x) + 0.5) : ((x) - 0.5))
+
 #define SWAP(type, a, b)  \
     {                     \
         type _temp = (a); \

include/math.h

@@ -12,12 +12,11 @@
 #define DEGTORAD(x) ((x) * (M_PI / 180.0f))
 
 typedef union {
+    f64 d;
     struct {
         u32 hi;
         u32 lo;
     } word;
-
-    f64 d;
 } du;
 
 typedef union {

include/variables.h

@@ -1670,7 +1670,7 @@ extern PadMgr* padmgrContext;
 // extern UNK_TYPE4 controllerInputsCaptured;
 // extern UNK_TYPE4 D_801D1538;
 extern UNK_PTR D_801D1540;
-extern f32 D_801D1570[13];
+// extern f32 sFactorialTbl[13];
 extern Vec3f D_801D15B0;
 extern Vec3s D_801D15BC;
 extern RSPMatrix D_801D1DE0;

spec

@@ -633,7 +633,6 @@ beginseg
     include "build/src/code/sys_math.o"
     include "build/src/code/sys_math3d.o"
     include "build/data/code/sys_math3d.bss.o"
-    include "build/data/code/code_801D1570.data.o"
     include "build/data/code/code_801D15B0.data.o"
     include "build/src/code/sys_math_atan.o"
     include "build/src/code/sys_matrix.o"

src/code/sys_math.c

@@ -1,15 +1,93 @@
+/**
+ * Maths library: two factorials, integer power, wrappers for libultra's sins and coss (the main ones used), and some
+ * random functions moved from OoT's z_actor
+ */
 #include "global.h"
 
-#pragma GLOBAL_ASM("asm/non_matchings/code/sys_math/func_80179300.s")
+static f32 sFactorialTbl[] = { 1.0f,    1.0f,     2.0f,      6.0f,       24.0f,       120.0f,      720.0f,
+                               5040.0f, 40320.0f, 362880.0f, 3628800.0f, 39916800.0f, 479001600.0f };
 
-#pragma GLOBAL_ASM("asm/non_matchings/code/sys_math/func_80179400.s")
+// Rename to Math_FactorialF
+/**
+ * Takes a float, returns the factorial of it(s trunctation), iteratively
+ * Unused
+ */
+f32 func_80179300(f32 n) {
+    f32 ret = 1.0f;
+    s32 i;
 
-#pragma GLOBAL_ASM("asm/non_matchings/code/sys_math/pow_int.s")
+    //! @bug No check for negative argument. Will return 1.0f if the argument truncates to a negative int.
+    for (i = n; i > 1; i--) {
+        ret *= i;
+    }
+    return ret;
+}
 
-#pragma GLOBAL_ASM("asm/non_matchings/code/sys_math/sin_rad.s")
+// Rename to Math_Factorial
+/**
+ * Takes an int and returns its factorial as a float. Uses the lookup table for n <= 12.
+ * Unused
+ */
+f32 func_80179400(s32 n) {
+    f32 ret;
+    s32 i;
 
-#pragma GLOBAL_ASM("asm/non_matchings/code/sys_math/cos_rad.s")
+    //! @bug No check for negative argument. Will read the array out-of-bounds if the argument is negative.
+    //! (The OoT version does an unsigned check instead, which will return sFactorialTbl[12] for a negative argument.)
+    if (n > 12) {
+        ret = sFactorialTbl[12];
+        for (i = 13; i <= n; i++) {
+            ret *= i;
+        }
+    } else {
+        ret = sFactorialTbl[n];
+    }
+    return ret;
+}
 
-#pragma GLOBAL_ASM("asm/non_matchings/code/sys_math/Rand_ZeroFloat.s")
+// Rename to Math_PowF
+/**
+ * Returns base^{exp} if exp is nonnegative and 1.0f otherwise.
+ * Unused
+ */
+f32 pow_int(f32 base, s32 exp) {
+    f32 ret = 1.0f;
 
-#pragma GLOBAL_ASM("asm/non_matchings/code/sys_math/randPlusMinusPoint5Scaled.s")
+    while (exp > 0) {
+        exp--;
+        ret *= base;
+    }
+
+    return ret;
+}
+
+// Rename to Math_SinF
+/**
+ * Takes an angle in radians and returns the sine.
+ */
+f32 sin_rad(f32 rad) {
+    return sins(RADF_TO_BINANG(rad)) * SHT_MINV;
+}
+
+// Rename to Math_CosF
+/**
+ * Takes an angle in radians and returns the cosine.
+ */
+f32 cos_rad(f32 rad) {
+    return coss(RADF_TO_BINANG(rad)) * SHT_MINV;
+}
+
+/**
+ * Returns a pseudo-random floating-point number between 0.0f and scale. Originally in z_actor in OoT.
+ */
+f32 Rand_ZeroFloat(f32 scale) {
+    return Rand_ZeroOne() * scale;
+}
+
+// Rename to Rand_CenteredFloat
+/**
+ * Returns a pseudo-random floating-point number between (- scale / 2) and (scale / 2). Originally in z_actor in OoT.
+ */
+f32 randPlusMinusPoint5Scaled(f32 scale) {
+    return Rand_Centered() * scale;
+}

src/libultra/gu/cosf.c

@@ -1,3 +1,70 @@
-#include "global.h"
+#include "math.h"
+#include "macros.h"
 
-#pragma GLOBAL_ASM("asm/non_matchings/boot/cosf/__cosf.s")
+// A slightly tweaked form of the coefficients of the Maclaurin series of sine up to x^9
+// [https://mathworld.wolfram.com/MaclaurinSeries.html].
+// The commented versions do not match.
+static const du P[] = {
+    { 1.0 },                     // 1
+    { -0.16666659550427756 },    // -1/3! = 1/6
+    { 0.008333066246082155 },    // 1/5! = 1/120
+    { -0.0001980960290193795 },  // -1/7! = -1/5040
+    { 0.000002605780637968037 }, // 1/9! = 1/362880
+};
+
+static const du rpi = { 1 / 3.14159265358979323846 }; // 1/M_PI, "reciprocal of pi"
+
+static const du pihi = { 3.1415926218032837 };
+
+static const du pilo = { 3.178650954705639E-8 }; // pihi + pilo is the closest double to pi
+
+static const fu zero = { 0x00000000 };
+
+/**
+ * Computes the cosine of a float, returning a float. It essentially computes sin(x+pi/2) by the same method as __sinf,
+ * without as many size checks.
+ */
+f32 __cosf(f32 x) {
+    f32 absx;
+    f64 dx;         // x promoted to double
+    f64 xSq;        // square of dx
+    f64 polyApprox; // Most of Maclaurin polynomial of sin(x) of degree 9
+    f64 dn;         // n promoted to double
+    s32 n;          // number of multiples of pi away from the first half-period
+    f64 result;
+    s32 ix = *(s32*)&x;   // Type-pun x into an s32, i.e. its IEEE-754 hex representation
+    s32 xpt = (ix >> 22); // Obtain the exponent of x (actually 2 * exponent + 127)
+
+    xpt &= 0x1FF; // Remove the sign bit
+
+    // |x| < 2^{28} (beyond this range, floats are too sparse to make the trig functions useable)
+    if (xpt < 310) {
+        absx = (x > 0) ? x : -x;
+        dx = absx;
+
+        dn = dx * rpi.d + 0.5;
+        n = ROUND(dn);
+        dn = n;
+
+        dn -= 0.5;
+
+        dx -= dn * pihi.d;
+        dx -= dn * pilo.d;
+
+        xSq = SQ(dx);
+        polyApprox = ((P[4].d * xSq + P[3].d) * xSq + P[2].d) * xSq + P[1].d;
+        result = dx + (dx * xSq) * polyApprox; // Actual Maclaurin polynomial for sin(x)
+
+        if (n % 2 == 0) {
+            return (f32)result;
+        }
+        return -(f32)result;
+    }
+
+    // if x is NaN
+    if (x != x) {
+        return __libm_qnan_f;
+    }
+
+    return zero.f;
+}

src/libultra/gu/coss.c

@@ -1,5 +1,8 @@
 #include "global.h"
 
-short coss(unsigned short x) {
+/**
+ * Compute the sine of a hex angle and return a short, using the formula cos(x) = sin(x+pi).
+ */
+s16 coss(u16 x) {
     return sins(x + 0x4000);
 }

src/libultra/gu/sinf.c

@@ -1,77 +1,79 @@
 #include "math.h"
+#include "macros.h"
 
+// A slightly tweaked form of the coefficients of the Maclaurin series of sine up to x^9
+// [https://mathworld.wolfram.com/MaclaurinSeries.html].
+// The commented versions do not match.
 static const du P[] = {
-    { 0x3FF00000, 0x00000000 }, { 0xBFC55554, 0xBC83656D }, { 0x3F8110ED, 0x3804C2A0 },
-    { 0xBF29F6FF, 0xEEA56814 }, { 0x3EC5DBDF, 0x0E314BFE },
+    { 1.0 },                     // 1
+    { -0.16666659550427756 },    // -1/3! = 1/6
+    { 0.008333066246082155 },    // 1/5! = 1/120
+    { -0.0001980960290193795 },  // -1/7! = -1/5040
+    { 0.000002605780637968037 }, // 1/9! = 1/362880
 };
 
-static const du rpi = { 0x3FD45F30, 0x6DC9C883 };
+static const du rpi = { 1 / 3.14159265358979323846 }; // 1/M_PI, "reciprocal of pi"
 
-static const du pihi = { 0x400921FB, 0x50000000 };
+static const du pihi = { 3.1415926218032837 };
 
-static const du pilo = { 0x3E6110B4, 0x611A6263 };
+static const du pilo = { 3.178650954705639E-8 }; // pihi + pilo is the closest double to pi
 
 static const fu zero = { 0x00000000 };
 
-extern float __libm_qnan_f;
+/**
+ * Returns the sine of a float as a float, using the Maclaurin series and shifting.
+ */
+f32 __sinf(f32 x) {
+    f64 dx;         // x promoted to double
+    f64 xSq;        // square of dx
+    f64 polyApprox; // Most of Maclaurin polynomial of sin(x) of degree 9
+    f64 dn;         // n promoted to double
+    s32 n;          // number of multiples of pi away from the first half-period
+    f64 result;
+    s32 ix = *(s32*)&x;   // Type-pun x into an s32, i.e. its IEEE-754 hex representation
+    s32 xpt = (ix >> 22); // Obtain the exponent of x (actually 2 * exponent + 127)
 
-float __sinf(float x) {
-    double dx;  // double x
-    double xsq; // x squared
-    double poly;
-    double dn;
-    int n;
-    double result;
-    int ix; // int x
-    int xpt;
-
-    ix = *(int*)&x;
-    xpt = (ix >> 22) & 0x1FF;
+    xpt &= 0x1FF; // Remove the sign bit
 
+    // |x| < 1
     if (xpt < 255) {
         dx = x;
+
+        // |x| > 2^{-12}: for x smaller in magnitude than this, sin(x) - x is too small for a float to register the
+        // error
         if (xpt >= 230) {
-            xsq = dx * dx;
+            xSq = SQ(dx);
+            polyApprox = ((P[4].d * xSq + P[3].d) * xSq + P[2].d) * xSq + P[1].d;
 
-            poly = (((((P[4].d * xsq) + P[3].d) * xsq) + P[2].d) * xsq) + P[1].d;
-
-            result = ((dx * xsq) * poly) + dx;
-
-            return result;
-        } else {
-            return x;
+            result = dx + (dx * xSq) * polyApprox;
+            return (f32)result;
         }
+        return x;
     }
 
+    // |x| < 2^{28} (beyond this range, floats are too sparse to make the trig functions useable)
     if (xpt < 310) {
         dx = x;
-
         dn = dx * rpi.d;
-
-        if (dn >= 0) {
-            n = dn + 0.5;
-        } else {
-            n = dn - 0.5;
-        }
-
+        n = ROUND(dn);
         dn = n;
 
+        // Bring x to the first half-period where the Maclaurin polynomial can be used
         dx -= dn * pihi.d;
         dx -= dn * pilo.d;
 
-        xsq = dx * dx;
+        xSq = SQ(dx);
+        polyApprox = ((P[4].d * xSq + P[3].d) * xSq + P[2].d) * xSq + P[1].d;
+        result = dx + (dx * xSq) * polyApprox; // Actual Maclaurin polynomial for sin(x)
 
-        poly = (((((P[4].d * xsq) + P[3].d) * xsq) + P[2].d) * xsq) + P[1].d;
-
-        result = ((dx * xsq) * poly) + dx;
-
-        if ((n & 0x1) == 0) {
-            return result;
-        } else {
-            return -(float)result;
+        // add a minus sign if x is an odd number of multiples of pi away from the first half-period
+        if (n % 2 == 0) {
+            return (f32)result;
         }
+        return -(f32)result;
     }
 
+    // if x is NaN
     if (x != x) {
         return __libm_qnan_f;
     }

src/libultra/gu/sins.c

@@ -77,8 +77,12 @@ static s16 sintable[0x400] = {
     0x7FFE, 0x7FFF,
 };
 
-short sins(unsigned short x) {
-    short val;
+/**
+ * Compute the sine of a hex angle and return a short, using a lookup table for the first quadrant, extrapolating to the
+ * others using the equivalents of sin(pi-x) = sin(x) and sin(pi+x) = -sin(x).
+ */
+s16 sins(u16 x) {
+    s16 val;
 
     x >>= 4;
     if ((x & 0x400) != 0) {

tools/disasm/files.txt

@@ -561,7 +561,7 @@
             0x801D1520 : "graph",
             0x801D1530 : "padmgr",
             0x801D1540 : "speed_meter",
-            0x801D1570 : "",
+            0x801D1570 : "sys_math",
             0x801D15B0 : "",
             0x801D15D0 : "sys_math_atan",
             0x801D1DE0 : "sys_matrix",