1 /* logf.c: Computes the natural log 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 */
31 #define SDCC_FLOAT_LIB
34 // TODO: share with other temps
35 static data unsigned char logf_tmp[4];
42 // extract the two input, placing it into:
43 // sign exponent mantiassa
44 // ---- -------- ---------
45 // x: sign_a exp_a r4/r3/r2
49 jnb sign_a, logf_zero_check
50 // TODO: set errno to EDOM (negative numbers not allowed)
55 // TODO: set errno to ERANGE (zero not allowed)
66 mov (_logf_tmp + 0), a // y = 0
67 mov (_logf_tmp + 1), a
68 mov (_logf_tmp + 2), a
69 mov (_logf_tmp + 3), a
70 mov dptr, #__fs_natural_log_table
73 mov ar7, r4 // r7/r6/r5 = x >> i
78 lcall __fs_cordic_rshift_r765_unsigned
79 mov a, r1 // check if x + (x >> i) > 1.0
89 mov a, r1 // x = x + (x >> i)
101 clr a // y = y + log_table[i]
103 add a, (_logf_tmp + 0)
104 mov (_logf_tmp + 0), a
107 addc a, (_logf_tmp + 1)
108 mov (_logf_tmp + 1), a
111 addc a, (_logf_tmp + 2)
112 mov (_logf_tmp + 2), a
115 addc a, (_logf_tmp + 3)
116 mov (_logf_tmp + 3), a
123 cjne r0, #30, logf_cordic_loop
124 // at this point, _logf_tmp has the natural log of the positive
125 // normalized fractional part (0.5 to 1.0 -> 0.693 to 0.0)
126 mov r4, (_logf_tmp + 3)
127 mov r3, (_logf_tmp + 2)
128 mov r2, (_logf_tmp + 1)
129 mov r1, (_logf_tmp + 0)
134 cjne a, #126, logf_exponent
135 // if the input exponent was 126, then we don't need to add
136 // anything and we can just return the with we have already
138 // TODO: which of these gives best accuracy???
139 ljmp fs_zerocheck_return
140 //ljmp fs_round_and_return
143 // the input exponent was greater than 126
149 // the input exponent was less than 126
154 // r0 has abs(exp - 126)
156 // put the log of faction into b, so we can use a to compute
157 // the offset to be added to it or subtracted from it
159 // multiply r0 by log(2), or r0 * 0xB17218
181 // turn r0 * log(2) into a proper float
184 // now just add log(fractional) +/- log(2) * abs(exp - 126)
185 ljmp fsadd_direct_entry
187 #pragma less_pedantic
191 #else // not MATH_ASM_MCS51
196 /*Constants for 24 bits or less (8 decimal digits)*/
197 #define A0 -0.5527074855E+0
198 #define B0 -0.6632718214E+1
202 #define C0 0.70710678118654752440
203 #define C1 0.693359375 /*355.0/512.0*/
204 #define C2 -2.121944400546905827679E-4
206 float logf(const float x) _FLOAT_FUNC_REENTRANT
208 #if defined(SDCC_mcs51) && defined(SDCC_MODEL_SMALL) \
209 && !defined(SDCC_NOOVERLAY)
213 float f, z, w, znum, zden, xn;
236 Rz=z+z*(w*A(w)/B(w));
238 return ((xn*C2+Rz)+xn*C1);