1 /* expf.c: Computes e**x of a 32-bit float as outlined in [1]
3 Copyright (C) 2001, 2002 Jesus Calvino-Fraga, jesusc@ieee.org
5 This library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
10 This library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
15 You should have received a copy of the GNU Lesser General Public
16 License along with this library; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
19 /* [1] William James Cody and W. M. Waite. _Software manual for the
20 elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */
22 /* Version 1.0 - Initial release */
32 #define SDCC_FLOAT_LIB
35 // TODO: share with other temps
37 static data unsigned char expf_y[4];
38 static data unsigned char n;
46 mov _sign_bit, c // remember sign
47 clr acc.7 // and make input positive
50 rlc a // extract exponent
53 // input is a very small number, so e^x is 1.000000
59 // TODO: check exponent for very small values, and return zero
62 add a, #0xE8 // is x >= 0.69314718
70 jnc expf_no_range_reduction
72 mov (_expf_y + 0), dpl // keep copy of x in "exp_y"
73 mov (_expf_y + 1), dph
84 lcall ___fsmul // x * 1.442695041 = x / ln(2)
89 lcall ___fs2uchar // n = int(x * 1.442695041)
94 ljmp fs_return_inf // exponent overflow
100 lcall expf_scale_and_add
101 mov (_expf_y + 0), dpl
102 mov (_expf_y + 1), dph
109 lcall expf_scale_and_add
110 expf_no_range_reduction:
113 // Compute e^x using the cordic algorithm. This works over an
114 // input range of 0 to 0.69314712. Can be extended to work from
115 // 0 to 1.0 if the results are normalized, but over the range
116 // we use, the result is always from 1.0 to 1.99999988 (fixed
121 rlc a // extract exponent to acc
123 mov r1, dpl // mantissa to r4/r3/r2/r1
128 // first, align the input into a 32 bit long
129 cjne a, #121, exp_lshift
133 // exp_a is greater than 121, so left shift
138 // exp_a is less than 121, so right shift
142 exp_noshift: // r4/r3/r2/r1 = x
144 mov (_expf_y + 0), a // y = 1.0;
147 mov (_expf_y + 3), #0x20
148 mov dptr, #__fs_natural_log_table
156 movc a, @a+dptr // r7/r6/r5/b = table[i]
169 subb a, b // compare x to table[i]
176 jc exp_cordic_skip // if table[i] < x
180 mov r1, a // x -= table[i]
191 mov r5, (_expf_y + 1) // r7/r6/r5/b = y >> i
192 mov r6, (_expf_y + 2)
193 mov r7, (_expf_y + 3)
195 lcall __fs_cordic_rshift_r765_unsigned
201 mov (_expf_y + 1), a // y += (y >> i)
210 cjne r0, #27, exp_cordic_loop
211 mov r4, (_expf_y + 3)
212 mov r3, (_expf_y + 2)
213 mov r2, (_expf_y + 1)
214 mov r1, (_expf_y + 0)
215 lcall fs_normalize_a // end of cordic
218 add a, _n // ldexpf(x, n)
220 lcall fs_round_and_return
222 jnb _sign_bit, expf_done
230 lcall ___fsdiv // 1.0 / x
237 #pragma less_pedantic
242 static void dummy1(void) _naked
252 lcall ___uchar2fs // turn n into float
253 lcall ___fsmul // n * scale_factor
262 mov dpl, (_expf_y + 0)
263 mov dph, (_expf_y + 1)
266 lcall ___fsadd // x += (n * scale_factor)
275 static void dummy(void) _naked
297 djnz r0, fs_lshift_loop
306 #else // not MATH_ASM_MCS51
309 #define P0 0.2499999995E+0
310 #define P1 0.4160288626E-2
311 #define Q0 0.5000000000E+0
312 #define Q1 0.4998717877E-1
314 #define P(z) ((P1*z)+P0)
315 #define Q(z) ((Q1*z)+Q0)
317 #define C1 0.693359375
318 #define C2 -2.1219444005469058277e-4
320 #define BIGX 88.72283911 /* ln(XMAX) */
321 #define EXPEPS 1.0E-7 /* exp(1.0E-7)=0.0000001 */
322 #define K1 1.4426950409 /* 1/ln(2) */
324 float expf(const float x)
327 float xn, g, r, z, y;
335 if(y<EXPEPS) return 1.0;