?login_element?

Subversion Repositories NedoOS

Rev

Blame | Last modification | View Log | Download

  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
  252.