/*-------------------------------------------------------------------------
expf.c - Computes e**x of a 32-bit float as outlined in [1]
Copyright (C) 2001, 2002, Jesus Calvino-Fraga, jesusc@ieee.org
This library is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this library; see the file COPYING. If not, write to the
Free Software Foundation, 51 Franklin Street, Fifth Floor, Boston,
MA 02110-1301, USA.
As a special exception, if you link this library with other files,
some of which are compiled with SDCC, to produce an executable,
this library does not by itself cause the resulting executable to
be covered by the GNU General Public License. This exception does
not however invalidate any other reasons why the executable file
might be covered by the GNU General Public License.
-------------------------------------------------------------------------*/
/* [1] William James Cody and W. M. Waite. _Software manual for the
elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */
/* Version 1.0 - Initial release */
#define __SDCC_MATH_LIB
#include <math.h>
#include <errno.h>
#include <stdbool.h>
#ifdef MATH_ASM_MCS51
#define __SDCC_FLOAT_LIB
#include <float.h>
// TODO: share with other temps
static __bit sign_bit;
static __data unsigned char expf_y[4];
static __data unsigned char n;
float expf(float x)
{
x;
__asm
mov c, acc.7
mov _sign_bit, c // remember sign
clr acc.7 // and make input positive
mov r0, a
mov c, b.7
rlc a // extract exponent
add a, #153
jc expf_not_zero
// input is a very small number, so e^x is 1.000000
mov dptr, #0
mov b, #0x80
mov a, #0x3F
ljmp expf_exit
expf_not_zero:
// TODO: check exponent for very small values, and return zero
mov _n, #0
mov a, dpl
add a, #0xE8 // is x >= 0.69314718
mov a, dph
addc a, #0x8D
mov a, b
addc a, #0xCE
mov a, r0
addc a, #0xC0
mov a, r0
jnc expf_no_range_reduction
expf_range_reduction:
mov (_expf_y + 0), dpl // keep copy of x in "exp_y"
mov (_expf_y + 1), dph
mov (_expf_y + 2), b
mov (_expf_y + 3), a
mov r0, #0x3B
push ar0
mov r0, #0xAA
push ar0
mov r0, #0xB8
push ar0
mov r0, #0x3F
push ar0
lcall ___fsmul // x * 1.442695041 = x / ln(2)
dec sp
dec sp
dec sp
dec sp
lcall ___fs2uchar // n = int(x * 1.442695041)
mov a, dpl
mov _n, a
add a, #128
jnc expf_range_ok
lcall fs_return_inf // exponent overflow
ljmp expf_exit
expf_range_ok:
mov r0,#0x00
mov r1,#0x80
mov r2,#0x31
mov r3,#0xBF
lcall expf_scale_and_add
mov (_expf_y + 0), dpl
mov (_expf_y + 1), dph
mov (_expf_y + 2), b
mov (_expf_y + 3), a
mov r0,#0x83
mov r1,#0x80
mov r2,#0x5E
mov r3,#0x39
lcall expf_scale_and_add
expf_no_range_reduction:
// Compute e^x using the cordic algorithm. This works over an
// input range of 0 to 0.69314712. Can be extended to work from
// 0 to 1.0 if the results are normalized, but over the range
// we use, the result is always from 1.0 to 1.99999988 (fixed
// exponent of 127)
expf_cordic_begin:
mov c, b.7
rlc a // extract exponent to acc
setb b.7
mov r1, dpl // mantissa to r4/r3/r2/r1
mov r2, dph
mov r3, b
mov r4, #0
// first, align the input into a 32 bit long
cjne a, #121, exp_lshift
sjmp exp_noshift
exp_lshift:
jc exp_rshift
// exp_a is greater than 121, so left shift
add a, #135
lcall fs_lshift_a
sjmp exp_noshift
exp_rshift:
// exp_a is less than 121, so right shift
cpl a
add a, #122
lcall fs_rshift_a
exp_noshift: // r4/r3/r2/r1 = x
clr a
mov (_expf_y + 0), a // y = 1.0;
mov (_expf_y + 1), a
mov (_expf_y + 2), a
mov (_expf_y + 3), #0x20
mov dptr, #__fs_natural_log_table
mov r0, a // r0 = i
exp_cordic_loop:
clr a
movc a, @a+dptr
mov b, a
inc dptr
clr a
movc a, @a+dptr // r7/r6/r5/b = table[i]
mov r5, a
inc dptr
clr a
movc a, @a+dptr
mov r6, a
inc dptr
clr a
movc a, @a+dptr
mov r7, a
inc dptr
clr c
mov a, r1
subb a, b // compare x to table[i]
mov a, r2
subb a, r5
mov a, r3
subb a, r6
mov a, r4
subb a, r7
jc exp_cordic_skip // if table[i] < x
clr c
mov a, r1
subb a, b
mov r1, a // x -= table[i]
mov a, r2
subb a, r5
mov r2, a
mov a, r3
subb a, r6
mov r3, a
mov a, r4
subb a, r7
mov r4, a
mov b, (_expf_y + 0)
mov r5, (_expf_y + 1) // r7/r6/r5/b = y >> i
mov r6, (_expf_y + 2)
mov r7, (_expf_y + 3)
mov a, r0
lcall __fs_cordic_rshift_r765_unsigned
mov a, (_expf_y + 0)
add a, b
mov (_expf_y + 0), a
mov a, (_expf_y + 1)
addc a, r5
mov (_expf_y + 1), a // y += (y >> i)
mov a, (_expf_y + 2)
addc a, r6
mov (_expf_y + 2), a
mov a, (_expf_y + 3)
addc a, r7
mov (_expf_y + 3), a
exp_cordic_skip:
inc r0
cjne r0, #27, exp_cordic_loop
mov r4, (_expf_y + 3)
mov r3, (_expf_y + 2)
mov r2, (_expf_y + 1)
mov r1, (_expf_y + 0)
mov exp_a, #121
lcall fs_normalize_a // end of cordic
mov a, #127
add a, _n // ldexpf(x, n)
mov exp_a, a
lcall fs_round_and_return
jnb _sign_bit, expf_done
push dpl
push dph
push b
push acc
mov dptr, #0
mov b, #0x80
mov a, #0x3F
lcall ___fsdiv // 1.0 / x
dec sp
dec sp
dec sp
dec sp
expf_done:
clr acc.7 // Result is always positive!
expf_exit:
__endasm;
#pragma less_pedantic
}
static void dummy1(void) __naked
{
__asm
.globl fs_lshift_a
expf_scale_and_add:
push ar0
push ar1
push ar2
push ar3
mov dpl, _n
lcall ___uchar2fs // turn n into float
lcall ___fsmul // n * scale_factor
dec sp
dec sp
dec sp
dec sp
push dpl
push dph
push b
push acc
mov dpl, (_expf_y + 0)
mov dph, (_expf_y + 1)
mov b, (_expf_y + 2)
mov a, (_expf_y + 3)
lcall ___fsadd // x += (n * scale_factor)
dec sp
dec sp
dec sp
dec sp
ret
__endasm;
}
static void dummy(void) __naked
{
__asm
.globl fs_lshift_a
fs_lshift_a:
jz fs_lshift_done
push ar0
mov r0, a
fs_lshift_loop:
clr c
mov a, r1
rlc a
mov r1, a
mov a, r2
rlc a
mov r2, a
mov a, r3
rlc a
mov r3, a
mov a, r4
rlc a
mov r4, a
djnz r0, fs_lshift_loop
pop ar0
fs_lshift_done:
ret
__endasm;
}
#else // not MATH_ASM_MCS51
#define P0 0.2499999995E+0
#define P1 0.4160288626E-2
#define Q0 0.5000000000E+0
#define Q1 0.4998717877E-1
#define P(z) ((P1*z)+P0)
#define Q(z) ((Q1*z)+Q0)
#define C1 0.693359375
#define C2 -2.1219444005469058277e-4
#define BIGX 88.72283911 /* ln(HUGE_VALF) */
#define EXPEPS 1.0E-7 /* exp(1.0E-7)=0.0000001 */
#define K1 1.4426950409 /* 1/ln(2) */
float expf(float x) _FLOAT_FUNC_REENTRANT
{
int n;
float xn, g, r, z, y;
bool sign;
if(x>=0.0)
{ y=x; sign=0; }
else
{ y=-x; sign=1; }
if(y<EXPEPS) return 1.0;
if(y>BIGX)
{
if(sign)
{
errno=ERANGE;
return HUGE_VALF
;
}
else
{
return 0.0;
}
}
z=y*K1;
n=z;
if(n<0) --n;
if(z-n>=0.5) ++n;
xn=n;
g=((y-xn*C1))-xn*C2;
z=g*g;
r=P(z)*g;
r=0.5+(r/(Q(z)-r));
n++;
z=ldexpf(r, n);
if(sign)
return 1.0/z;
else
return z;
}
#endif