sdcc/device/lib/logf.c

   1 /*-------------------------------------------------------------------------
   2    logf.c - Computes the natural log of a 32 bit float as outlined in [1].
   3
   4    Copyright (C) 2001, 2002, Jesus Calvino-Fraga, jesusc@ieee.org
   5
   6    This library is free software; you can redistribute it and/or modify it
   7    under the terms of the GNU General Public License as published by the
   8    Free Software Foundation; either version 2, or (at your option) any
   9    later version.
  10
  11    This library is distributed in the hope that it will be useful,
  12    but WITHOUT ANY WARRANTY; without even the implied warranty of
  13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14    GNU General Public License for more details.
  15
  16    You should have received a copy of the GNU General Public License
  17    along with this library; see the file COPYING. If not, write to the
  18    Free Software Foundation, 51 Franklin Street, Fifth Floor, Boston,
  19    MA 02110-1301, USA.
  20
  21    As a special exception, if you link this library with other files,
  22    some of which are compiled with SDCC, to produce an executable,
  23    this library does not by itself cause the resulting executable to
  24    be covered by the GNU General Public License. This exception does
  25    not however invalidate any other reasons why the executable file
  26    might be covered by the GNU General Public License.
  27 -------------------------------------------------------------------------*/
  28
  29 /* [1] William James Cody and W.  M.  Waite.  _Software manual for the
  30    elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */
  31
  32 /* Version 1.0 - Initial release */
  33
  34 #define __SDCC_MATH_LIB
  35 #include <math.h>
  36 #include <errno.h>
  37
  38
  39 #ifdef MATH_ASM_MCS51
  40
  41 #define __SDCC_FLOAT_LIB
  42 #include <float.h>
  43
  44 // TODO: share with other temps
  45 static __data unsigned char logf_tmp[4];
  46
  47 float logf(float x)
  48 {
  49         x;
  50         __asm
  51
  52         // extract the two input, placing it into:
  53         //      sign     exponent   mantiassa
  54         //      ----     --------   ---------
  55         //  x:  sign_a   exp_a     r4/r3/r2
  56
  57         lcall   fsgetarg
  58 logf_neg_check:
  59         jnb     sign_a, logf_zero_check
  60         // TODO: set errno to EDOM (negative numbers not allowed)
  61         lcall   fs_return_nan
  62         ljmp    logf_exit
  63
  64 logf_zero_check:
  65         cjne    r4, #0, logf_ok
  66         // TODO: set errno to ERANGE (zero not allowed)
  67         setb    sign_a
  68         lcall   fs_return_inf
  69         ljmp    logf_exit
  70
  71 logf_ok:
  72         push    exp_a
  73         mov     a, #3
  74         mov     r1, #0
  75         lcall   fs_rshift_a
  76
  77         clr     a
  78         mov     (_logf_tmp + 0), a      // y = 0
  79         mov     (_logf_tmp + 1), a
  80         mov     (_logf_tmp + 2), a
  81         mov     (_logf_tmp + 3), a
  82         mov     dptr, #__fs_natural_log_table
  83         mov     r0, a
  84 logf_cordic_loop:
  85         mov     ar7, r4         // r7/r6/r5 = x >> i
  86         mov     ar6, r3
  87         mov     ar5, r2
  88         mov     b, r1
  89         mov     a, r0
  90         lcall   __fs_cordic_rshift_r765_unsigned
  91         mov     a, r1           // check if x + (x >> i) > 1.0
  92         add     a, b
  93         mov     a, r2
  94         addc    a, r5
  95         mov     a, r3
  96         addc    a, r6
  97         mov     a, r4
  98         addc    a, r7
  99         anl     a, #0xE0
 100         jnz     logf_cordic_skip
 101         mov     a, r1           // x = x + (x >> i)
 102         add     a, b
 103         mov     r1, a
 104         mov     a, r2
 105         addc    a, r5
 106         mov     r2, a
 107         mov     a, r3
 108         addc    a, r6
 109         mov     r3, a
 110         mov     a, r4
 111         addc    a, r7
 112         mov     r4, a
 113         clr     a               // y = y + log_table[i]
 114         movc    a, @a+dptr
 115         add     a, (_logf_tmp + 0)
 116         mov     (_logf_tmp + 0), a
 117         mov     a, #1
 118         movc    a, @a+dptr
 119         addc    a, (_logf_tmp + 1)
 120         mov     (_logf_tmp + 1), a
 121         mov     a, #2
 122         movc    a, @a+dptr
 123         addc    a, (_logf_tmp + 2)
 124         mov     (_logf_tmp + 2), a
 125         mov     a, #3
 126         movc    a, @a+dptr
 127         addc    a, (_logf_tmp + 3)
 128         mov     (_logf_tmp + 3), a
 129 logf_cordic_skip:
 130         inc     dptr
 131         inc     dptr
 132         inc     dptr
 133         inc     dptr
 134         inc     r0
 135         cjne    r0, #30, logf_cordic_loop
 136         // at this point, _logf_tmp has the natural log of the positive
 137         // normalized fractional part (0.5 to 1.0 -> 0.693 to 0.0)
 138         mov     r4, (_logf_tmp + 3)
 139         mov     r3, (_logf_tmp + 2)
 140         mov     r2, (_logf_tmp + 1)
 141         mov     r1, (_logf_tmp + 0)
 142         mov     exp_a, #129
 143         setb    sign_a
 144         lcall   fs_normalize_a
 145         pop     acc
 146         cjne    a, #126, logf_exponent
 147         // if the input exponent was 126, then we don't need to add
 148         // anything and we can just return the with we have already
 149
 150         // TODO: which of these gives best accuracy???
 151         lcall   fs_zerocheck_return
 152         //lcall fs_round_and_return
 153         sjmp    logf_exit
 154 logf_exponent:
 155         jc      logf_exp_neg
 156         // the input exponent was greater than 126
 157 logf_exp_pos:
 158         add     a, #130
 159         clr     sign_b
 160         sjmp    logf_exp_scale
 161 logf_exp_neg:
 162         // the input exponent was less than 126
 163         cpl     a
 164         add     a, #127
 165         setb    sign_b
 166 logf_exp_scale:
 167         // r0 has abs(exp - 126)
 168         mov     r0, a
 169         // put the log of faction into b, so we can use a to compute
 170         // the offset to be added to it or subtracted from it
 171         lcall   fs_swap_a_b
 172         // multiply r0 by log(2), or r0 * 0xB17218
 173         mov     a, #0x18
 174         mov     b, r0
 175         mul     ab
 176         mov     r1, a
 177         mov     r2, b
 178         mov     a, #0xB1
 179         mov     b, r0
 180         mul     ab
 181         mov     r3, a
 182         mov     r4, b
 183         mov     a, #0x72
 184         mov     b, r0
 185         mul     ab
 186         add     a, r2
 187         mov     r2, a
 188         mov     a, b
 189         addc    a, r3
 190         mov     r3, a
 191         clr     a
 192         addc    a, r4
 193         mov     r4, a
 194         // turn r0 * log(2) into a proper float
 195         mov     exp_a, #134
 196         lcall   fs_normalize_a
 197         // now just add log(fractional) +/- log(2) * abs(exp - 126)
 198         lcall   fsadd_direct_entry
 199 logf_exit:
 200         __endasm;
 201 #pragma less_pedantic
 202 }
 203
 204 #else // not MATH_ASM_MCS51
 205
 206 /*Constants for 24 bits or less (8 decimal digits)*/
 207 #define A0 -0.5527074855E+0
 208 #define B0 -0.6632718214E+1
 209 #define A(w) (A0)
 210 #define B(w) (w+B0)
 211
 212 #define C0  0.70710678118654752440
 213 #define C1  0.693359375 /*355.0/512.0*/
 214 #define C2 -2.121944400546905827679E-4
 215
 216 float logf(float x) _FLOAT_FUNC_REENTRANT
 217 {
 218 #if     defined(__SDCC_mcs51) && defined(__SDCC_MODEL_SMALL) \
 219     && !defined(__SDCC_NOOVERLAY)
 220     volatile
 221 #endif
 222     float Rz;
 223     float f, z, w, znum, zden, xn;
 224     int n;
 225
 226     if (x<=0.0)
 227     {
 228         errno=EDOM;
 229         return 0.0;
 230     }
 231     f=frexpf(x, &n);
 232     znum=f-0.5;
 233     if (f>C0)
 234     {
 235         znum-=0.5;
 236         zden=(f*0.5)+0.5;
 237     }
 238     else
 239     {
 240         n--;
 241         zden=znum*0.5+0.5;
 242     }
 243     z=znum/zden;
 244     w=z*z;
 245
 246     Rz=z+z*(w*A(w)/B(w));
 247     xn=n;
 248     return ((xn*C2+Rz)+xn*C1);
 249 }
 250
 251 #endif