blob: 1c3a70be4c00ba76c25ac0c72305563c49fdea91 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
#include "ultra64.h"
#include "math.h"
#include "macros.h"
#pragma weak sinf = __sinf
// Coefficients of a degree 9 polynomial approximation of sine. It is not the Maclaurin polynamial, but some as-yet
// undetermined more uniform approximation.
static const du P[] = {
{ 1.0 },
{ -0.16666659550427756 },
{ 0.008333066246082155 },
{ -0.0001980960290193795 },
{ 0.000002605780637968037 },
};
static const du rpi = { 1 / M_PI }; // "reciprocal of pi"
// pihi + pilo is the closest double to pi, this representation allows more precise calculations since pi itself is not
// an exact float
static const du pihi = { 3.1415926218032837 };
static const du pilo = { 3.178650954705639E-8 };
static const fu zero = { 0x00000000 };
/**
* 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 the polynomial approximation to sin(x)
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| < 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;
polyApprox = ((P[4].d * xSq + P[3].d) * xSq + P[2].d) * xSq + P[1].d;
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 usable)
if (xpt < 310) {
dx = x;
dn = dx * rpi.d;
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;
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)
// 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;
}
return zero.f;
}
|
