git.gag.com Git - debian/cc1111/blob - device/lib/logf.c

   1 /*  logf.c: Computes the natural log of a 32 bit float as outlined in [1].
   2
   3     Copyright (C) 2001, 2002  Jesus Calvino-Fraga, jesusc@ieee.org
   4
   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.
   9
  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.
  14
  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 */
  18
  19 /* [1] William James Cody and W.  M.  Waite.  _Software manual for the
  20    elementary functions_, Englewood Cliffs, N.J.:Prentice-Hall, 1980. */
  21
  22 /* Version 1.0 - Initial release */
  23
  24 #define SDCC_MATH_LIB
  25 #include <math.h>
  26 #include <errno.h>
  27
  28
  29 #ifdef MATH_ASM_MCS51
  30
  31 #define SDCC_FLOAT_LIB
  32 #include <float.h>
  33
  34 // TODO: share with other temps
  35 static __data unsigned char logf_tmp[4];
  36
  37 float logf(float x)
  38 {
  39         x;
  40         __asm
  41
  42         // extract the two input, placing it into:
  43         //      sign     exponent   mantiassa
  44         //      ----     --------   ---------
  45         //  x:  sign_a   exp_a     r4/r3/r2
  46
  47         lcall   fsgetarg
  48 logf_neg_check:
  49         jnb     sign_a, logf_zero_check
  50         // TODO: set errno to EDOM (negative numbers not allowed)
  51         ljmp    fs_return_nan
  52
  53 logf_zero_check:
  54         cjne    r4, #0, logf_ok
  55         // TODO: set errno to ERANGE (zero not allowed)
  56         setb    sign_a
  57         ljmp    fs_return_inf
  58
  59 logf_ok:
  60         push    exp_a
  61         mov     a, #3
  62         mov     r1, #0
  63         lcall   fs_rshift_a
  64
  65         clr     a
  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
  71         mov     r0, a
  72 logf_cordic_loop:
  73         mov     ar7, r4         // r7/r6/r5 = x >> i
  74         mov     ar6, r3
  75         mov     ar5, r2
  76         mov     b, r1
  77         mov     a, r0
  78         lcall   __fs_cordic_rshift_r765_unsigned
  79         mov     a, r1           // check if x + (x >> i) > 1.0
  80         add     a, b
  81         mov     a, r2
  82         addc    a, r5
  83         mov     a, r3
  84         addc    a, r6
  85         mov     a, r4
  86         addc    a, r7
  87         anl     a, #0xE0
  88         jnz     logf_cordic_skip
  89         mov     a, r1           // x = x + (x >> i)
  90         add     a, b
  91         mov     r1, a
  92         mov     a, r2
  93         addc    a, r5
  94         mov     r2, a
  95         mov     a, r3
  96         addc    a, r6
  97         mov     r3, a
  98         mov     a, r4
  99         addc    a, r7
 100         mov     r4, a
 101         clr     a               // y = y + log_table[i]
 102         movc    a, @a+dptr
 103         add     a, (_logf_tmp + 0)
 104         mov     (_logf_tmp + 0), a
 105         mov     a, #1
 106         movc    a, @a+dptr
 107         addc    a, (_logf_tmp + 1)
 108         mov     (_logf_tmp + 1), a
 109         mov     a, #2
 110         movc    a, @a+dptr
 111         addc    a, (_logf_tmp + 2)
 112         mov     (_logf_tmp + 2), a
 113         mov     a, #3
 114         movc    a, @a+dptr
 115         addc    a, (_logf_tmp + 3)
 116         mov     (_logf_tmp + 3), a
 117 logf_cordic_skip:
 118         inc     dptr
 119         inc     dptr
 120         inc     dptr
 121         inc     dptr
 122         inc     r0
 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)
 130         mov     exp_a, #129
 131         setb    sign_a
 132         lcall   fs_normalize_a
 133         pop     acc
 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
 137
 138         // TODO: which of these gives best accuracy???
 139         ljmp    fs_zerocheck_return
 140         //ljmp  fs_round_and_return
 141 logf_exponent:
 142         jc      logf_exp_neg
 143         // the input exponent was greater than 126
 144 logf_exp_pos:
 145         add     a, #130
 146         clr     sign_b
 147         sjmp    logf_exp_scale
 148 logf_exp_neg:
 149         // the input exponent was less than 126
 150         cpl     a
 151         add     a, #127
 152         setb    sign_b
 153 logf_exp_scale:
 154         // r0 has abs(exp - 126)
 155         mov     r0, a
 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
 158         lcall   fs_swap_a_b
 159         // multiply r0 by log(2), or r0 * 0xB17218
 160         mov     a, #0x18
 161         mov     b, r0
 162         mul     ab
 163         mov     r1, a
 164         mov     r2, b
 165         mov     a, #0xB1
 166         mov     b, r0
 167         mul     ab
 168         mov     r3, a
 169         mov     r4, b
 170         mov     a, #0x72
 171         mov     b, r0
 172         mul     ab
 173         add     a, r2
 174         mov     r2, a
 175         mov     a, b
 176         addc    a, r3
 177         mov     r3, a
 178         clr     a
 179         addc    a, r4
 180         mov     r4, a
 181         // turn r0 * log(2) into a proper float
 182         mov     exp_a, #134
 183         lcall   fs_normalize_a
 184         // now just add log(fractional) +/- log(2) * abs(exp - 126)
 185         ljmp    fsadd_direct_entry
 186         __endasm;
 187 #pragma less_pedantic
 188 }
 189
 190 #else // not MATH_ASM_MCS51
 191
 192 /*Constants for 24 bits or less (8 decimal digits)*/
 193 #define A0 -0.5527074855E+0
 194 #define B0 -0.6632718214E+1
 195 #define A(w) (A0)
 196 #define B(w) (w+B0)
 197
 198 #define C0  0.70710678118654752440
 199 #define C1  0.693359375 /*355.0/512.0*/
 200 #define C2 -2.121944400546905827679E-4
 201
 202 float logf(const float x) _FLOAT_FUNC_REENTRANT
 203 {
 204 #if     defined(SDCC_mcs51) && defined(SDCC_MODEL_SMALL) \
 205     && !defined(SDCC_NOOVERLAY)
 206     volatile
 207 #endif
 208     float Rz;
 209     float f, z, w, znum, zden, xn;
 210     int n;
 211
 212     if (x<=0.0)
 213     {
 214         errno=EDOM;
 215         return 0.0;
 216     }
 217     f=frexpf(x, &n);
 218     znum=f-0.5;
 219     if (f>C0)
 220     {
 221         znum-=0.5;
 222         zden=(f*0.5)+0.5;
 223     }
 224     else
 225     {
 226         n--;
 227         zden=znum*0.5+0.5;
 228     }
 229     z=znum/zden;
 230     w=z*z;
 231
 232     Rz=z+z*(w*A(w)/B(w));
 233     xn=n;
 234     return ((xn*C2+Rz)+xn*C1);
 235 }
 236
 237 #endif