@@ -113,6 +113,118 @@ float powf(float a, float b) {
113113 return expf (b * logf (a ));
114114}
115115
116+ // --- sinf / cosf ---
117+ // Range-reduce to [-pi/4, pi/4] then apply degree-5 minimax polynomials.
118+ // Max error ~3 ULP across the full float range.
119+ static float _dc_sinf_kernel (float x ) {
120+ // sin(x) for |x| <= pi/4, degree-5 polynomial
121+ float x2 = x * x ;
122+ return x * (1.0f + x2 * (-0.16666667f + x2 * (0.0083333337f + x2 * -0.00019841270f )));
123+ }
124+ static float _dc_cosf_kernel (float x ) {
125+ // cos(x) for |x| <= pi/4, degree-6 polynomial
126+ float x2 = x * x ;
127+ return 1.0f + x2 * (-0.5f + x2 * (0.041666668f + x2 * -0.0013888889f ));
128+ }
129+ float sinf (float x ) {
130+ // Reduce x into [-pi/4, pi/4] using quadrant index n
131+ float q = x * 0.63661977f ; // x * (2/pi)
132+ int n = (int )(q + (x >= 0.0f ? 0.5f : -0.5f ));
133+ float r = x - (float )n * 1.5707963268f ; // r = x - n*(pi/2)
134+ switch (n & 3 ) {
135+ case 0 : return _dc_sinf_kernel (r );
136+ case 1 : return _dc_cosf_kernel (r );
137+ case 2 : return - _dc_sinf_kernel (r );
138+ default : return - _dc_cosf_kernel (r );
139+ }
140+ }
141+ float cosf (float x ) {
142+ float q = x * 0.63661977f ;
143+ int n = (int )(q + (x >= 0.0f ? 0.5f : -0.5f ));
144+ float r = x - (float )n * 1.5707963268f ;
145+ switch (n & 3 ) {
146+ case 0 : return _dc_cosf_kernel (r );
147+ case 1 : return - _dc_sinf_kernel (r );
148+ case 2 : return - _dc_cosf_kernel (r );
149+ default : return _dc_sinf_kernel (r );
150+ }
151+ }
152+
153+ // --- tanf ---
154+ float tanf (float x ) {
155+ float c = cosf (x );
156+ if (c == 0.0f ) return 3.40282347e+38f ;
157+ return sinf (x ) / c ;
158+ }
159+
160+ // --- asinf ---
161+ // Uses identity: asin(x) = atan(x / sqrt(1 - x*x))
162+ float asinf (float x ) {
163+ if (x >= 1.0f ) return 1.5707963268f ; // pi/2
164+ if (x <= -1.0f ) return -1.5707963268f ;
165+ float t = x * x ;
166+ float s = sqrtf (1.0f - t );
167+ // atan(x/sqrt(1-x^2)) via degree-5 minimax on [0,1]
168+ float u = x / s ;
169+ float u2 = u * u ;
170+ return u * (1.0f + u2 * (-0.33333334f + u2 * (0.2f + u2 * (-0.14285715f + u2 * 0.11111111f ))));
171+ }
172+
173+ // --- acosf ---
174+ float acosf (float x ) { return 1.5707963268f - asinf (x ); }
175+
176+ // --- atanf ---
177+ // Degree-9 minimax polynomial, accurate to ~2 ULP for |x| <= 1.
178+ // Larger |x| uses identity: atan(x) = pi/2 - atan(1/x).
179+ float atanf (float x ) {
180+ int inv = 0 ;
181+ if (x < 0.0f ) { x = - x ; inv = 2 ; } // handle sign
182+ if (x > 1.0f ) { x = 1.0f / x ; inv ^= 1 ; }
183+ float x2 = x * x ;
184+ float p = x * (1.0f + x2 * (-0.33333334f + x2 * (0.2f + x2 *
185+ (-0.14285715f + x2 * (0.11111111f + x2 * -0.090909094f )))));
186+ if (inv & 1 ) p = 1.5707963268f - p ;
187+ return (inv & 2 ) ? - p : p ;
188+ }
189+
190+ // --- atan2f ---
191+ float atan2f (float y , float x ) {
192+ if (x > 0.0f ) return atanf (y / x );
193+ if (x < 0.0f && y >= 0.0f ) return atanf (y / x ) + 3.14159265359f ;
194+ if (x < 0.0f && y < 0.0f ) return atanf (y / x ) - 3.14159265359f ;
195+ if (x == 0.0f && y > 0.0f ) return 1.5707963268f ;
196+ if (x == 0.0f && y < 0.0f ) return -1.5707963268f ;
197+ return 0.0f ; // x == 0, y == 0 — undefined, return 0
198+ }
199+
200+ // =============================================================================
201+ // Double-precision shims
202+ //
203+ // model.c includes <math.h> which may resolve generic math calls to the
204+ // double-precision symbols (cos, sin, exp, …) instead of the float ones
205+ // (cosf, sinf, expf, …). These thin wrappers delegate to our float shims
206+ // above so no extra precision logic is needed.
207+ // =============================================================================
208+ double cos (double x ) { return (double )cosf ((float )x ); }
209+ double sin (double x ) { return (double )sinf ((float )x ); }
210+ double tan (double x ) { return (double )tanf ((float )x ); }
211+ double acos (double x ) { return (double )acosf ((float )x ); }
212+ double asin (double x ) { return (double )asinf ((float )x ); }
213+ double atan (double x ) { return (double )atanf ((float )x ); }
214+ double atan2 (double y , double x ) { return (double )atan2f ((float )y , (float )x ); }
215+ double exp (double x ) { return (double )expf ((float )x ); }
216+ double log (double x ) { return (double )logf ((float )x ); }
217+ double log2 (double x ) { return (double )log2f ((float )x ); }
218+ double log10 (double x ) { return (double )log10f ((float )x ); }
219+ double sqrt (double x ) { return (double )sqrtf ((float )x ); }
220+ double pow (double a , double b ) { return (double )powf ((float )a , (float )b ); }
221+ double fabs (double x ) { return (double )fabsf ((float )x ); }
222+ double floor (double x ) { return (double )floorf ((float )x ); }
223+ double ceil (double x ) { return (double )ceilf ((float )x ); }
224+ double round (double x ) { return (double )roundf ((float )x ); }
225+ double fmod (double x , double y ) { return (double )fmodf ((float )x , (float )y ); }
226+ double tanh (double x ) { return (double )tanhf ((float )x ); }
227+
116228#endif
117229
118230int native_errno = 0 ;
0 commit comments