diff --git a/include/convert.h b/include/convert.h index 31bdb034f7..456e4b1982 100644 --- a/include/convert.h +++ b/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) diff --git a/include/functions.h b/include/functions.h index 0948b502e9..aba291a288 100644 --- a/include/functions.h +++ b/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); diff --git a/include/macros.h b/include/macros.h index 36af1534bc..e0e6c8e7c5 100644 --- a/include/macros.h +++ b/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); \ diff --git a/include/math.h b/include/math.h index 316bc57102..20fa352faa 100644 --- a/include/math.h +++ b/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 { diff --git a/include/variables.h b/include/variables.h index f2fd0f2b4e..c6fb559823 100644 --- a/include/variables.h +++ b/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; diff --git a/spec b/spec index 0754aca9c9..6b028302d3 100644 --- a/spec +++ b/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" diff --git a/src/code/sys_math.c b/src/code/sys_math.c index 4142cd5982..9b274c7e24 100644 --- a/src/code/sys_math.c +++ b/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; +} diff --git a/src/libultra/gu/cosf.c b/src/libultra/gu/cosf.c index c5af1c9076..72d154d700 100644 --- a/src/libultra/gu/cosf.c +++ b/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; +} diff --git a/src/libultra/gu/coss.c b/src/libultra/gu/coss.c index 466c36505a..ebed284053 100644 --- a/src/libultra/gu/coss.c +++ b/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); } diff --git a/src/libultra/gu/sinf.c b/src/libultra/gu/sinf.c index 01923221c0..5f3455ec54 100644 --- a/src/libultra/gu/sinf.c +++ b/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; } diff --git a/src/libultra/gu/sins.c b/src/libultra/gu/sins.c index 6dd5f23fcf..69b3b7e391 100644 --- a/src/libultra/gu/sins.c +++ b/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) { diff --git a/tools/disasm/files.txt b/tools/disasm/files.txt index c67e0b7a7c..eb4be6eb39 100644 --- a/tools/disasm/files.txt +++ b/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",