@@ -7,6 +7,112 @@ void *memcpy(void *dst, const void *src, size_t n) {
77void * memset (void * s , int c , size_t n ) {
88 return mp_fun_table .memset_ (s , c , n );
99}
10+
11+ // =============================================================================
12+ // libm shims — safe replacements for arm-none-eabi Newlib math functions.
13+ //
14+ // libm.a is NOT linked (LIBM_PATH := in Makefile) because mpy_ld cannot
15+ // resolve libm's internal Newlib dependencies (__errno, _REENT, kernel
16+ // helpers). Calling an unresolved libm symbol causes a silent HardFault.
17+ //
18+ // These shims use only integer ops and IEEE 754 bit manipulation — no
19+ // external dependencies. They cover all single-precision functions
20+ // commonly used by DEEPCRAFT / TensorFlow Lite inference models.
21+ //
22+ // If a future model needs a function not listed here, the build will fail
23+ // with a clear "undefined symbol" link error.
24+ // =============================================================================
25+
26+ // Helper: bit-cast float - unsigned int without UB
27+ static inline unsigned int _f2u (float f ) { union { float f ; unsigned int u ; } v ; v .f = f ; return v .u ; }
28+ static inline float _u2f (unsigned int u ) { union { float f ; unsigned int u ; } v ; v .u = u ; return v .f ; }
29+
30+ // --- expf ---
31+ // Degree-5 minimax polynomial + IEEE 754 exponent scaling.
32+ // Max error ~1.5 ULP for |x| <= 88.
33+ float expf (float x ) {
34+ if (x > 88.0f ) return 3.40282347e+38f ;
35+ if (x < -88.0f ) return 0.0f ;
36+ int n = (int )(x * 1.4426950409f + (x >= 0.0f ? 0.5f : -0.5f ));
37+ float r = x - (float )n * 0.6931471806f ;
38+ float p = 1.0f + r * (1.0f + r * (0.5f + r * (0.16666667f +
39+ r * (0.041666667f + r * 0.0083333333f ))));
40+ return p * _u2f ((unsigned int )((n + 127 ) << 23 ));
41+ }
42+
43+ // --- logf ---
44+ // Natural log via exponent extraction + polynomial on [sqrt(0.5), sqrt(2)].
45+ float logf (float x ) {
46+ if (x <= 0.0f ) return -3.40282347e+38f ;
47+ int e = (int )((_f2u (x ) >> 23 ) & 0xFF ) - 127 ;
48+ float m = _u2f ((_f2u (x ) & 0x007FFFFFu ) | 0x3F000000u ); // m in [0.5, 1)
49+ if (m < 0.70710678f ) { m *= 2.0f ; e -= 1 ; }
50+ float f = (m - 1.0f ) / (m + 1.0f ), f2 = f * f ;
51+ float p = f * (2.0f + f2 * (0.66666667f + f2 * (0.4f + f2 *
52+ (0.28571429f + f2 * 0.22222222f ))));
53+ return p + (float )e * 0.6931471806f ;
54+ }
55+
56+ // --- log2f ---
57+ float log2f (float x ) { return logf (x ) * 1.4426950409f ; }
58+
59+ // --- log10f ---
60+ float log10f (float x ) { return logf (x ) * 0.4342944819f ; }
61+
62+ // --- tanhf ---
63+ // tanh(x) = (e^2x - 1)/(e^2x + 1); clamped for |x| > 9 (result = ±1).
64+ float tanhf (float x ) {
65+ if (x > 9.0f ) return 1.0f ;
66+ if (x < -9.0f ) return -1.0f ;
67+ float e2x = expf (2.0f * x );
68+ return (e2x - 1.0f ) / (e2x + 1.0f );
69+ }
70+
71+ // --- sqrtf ---
72+ // Three Newton-Raphson iterations on an IEEE 754 seed estimate (~6 correct bits).
73+ // On Cortex-M4F the compiler may emit vsqrt.f32 directly and never call this.
74+ float sqrtf (float x ) {
75+ if (x <= 0.0f ) return 0.0f ;
76+ float r = _u2f (0x1fbb4000u + (_f2u (x ) >> 1 )); // seed: ~half the exponent
77+ r = 0.5f * (r + x / r );
78+ r = 0.5f * (r + x / r );
79+ r = 0.5f * (r + x / r );
80+ return r ;
81+ }
82+
83+ // --- fabsf --- (often compiled to vabs.f32 inline, but shim for safety)
84+ float fabsf (float x ) { return _u2f (_f2u (x ) & 0x7FFFFFFFu ); }
85+
86+ // --- floorf ---
87+ float floorf (float x ) { int i = (int )x ; return (float )(i - (x < (float )i )); }
88+
89+ // --- ceilf ---
90+ float ceilf (float x ) { int i = (int )x ; return (float )(i + (x > (float )i )); }
91+
92+ // --- roundf ---
93+ float roundf (float x ) { return (x >= 0.0f ) ? floorf (x + 0.5f ) : ceilf (x - 0.5f ); }
94+
95+ // --- fmodf ---
96+ float fmodf (float x , float y ) {
97+ if (y == 0.0f ) return 0.0f ;
98+ float q = x / y ;
99+ int n = (int )q ;
100+ return x - (float )n * y ;
101+ }
102+
103+ // --- powf ---
104+ // a^b = exp(b * ln(a)); handles integer-exponent cases first for accuracy.
105+ float powf (float a , float b ) {
106+ if (b == 0.0f ) return 1.0f ;
107+ if (a == 0.0f ) return 0.0f ;
108+ if (a < 0.0f ) {
109+ int ib = (int )b ;
110+ if ((float )ib != b ) return 0.0f ; // non-integer exponent of negative base
111+ return (ib & 1 ) ? - expf (b * logf (- a )) : expf (b * logf (- a ));
112+ }
113+ return expf (b * logf (a ));
114+ }
115+
10116#endif
11117
12118int native_errno = 0 ;
@@ -37,7 +143,7 @@ static mp_obj_t dc_make_new(const mp_obj_type_t *type, size_t n_args, size_t n_k
37143mp_obj_t mpy_init (mp_obj_fun_bc_t * self , size_t n_args , size_t n_kw , mp_obj_t * args ) {
38144 // This must be first, it sets up the globals dict and other things
39145 MP_DYNRUNTIME_INIT_ENTRY
40-
146+
41147 // Populate type
42148 dc_type .base .type = (void * )& mp_type_type ;
43149 dc_type .flags = MP_TYPE_FLAG_NONE ;
@@ -57,4 +163,4 @@ mp_obj_t mpy_init(mp_obj_fun_bc_t *self, size_t n_args, size_t n_kw, mp_obj_t *a
57163 // This must be last, it restores the globals dict
58164 MP_DYNRUNTIME_INIT_EXIT
59165
60- }
166+ }
0 commit comments