please dont rip this site

PIC Microcontoller Math Method

X ^ Y Logarithm ruler approach

Nikolai Golovchenko says:

[to implement x^y try the] logarithm ruler approach.
   z = x^y
log2 z = log2 x^y
log2 z = y * log2 x

   z = 2^(y*log2(x))
   ================

You have a couple of options to calculate that expression:

1) approximate calculation of log2 and 2^x with a small table and linear approximation (probably just 16-32 entry table would give a less than 0.5% error).

That would take under 1000 cycles and about 500 words.

See www.dattalo.com for an 8 bit log routine. Exp routine is just the log routine reversed.

2) Using CORDIC method for log2 and 2^x. This is probably a more precise method, but will take ~5000 cycles and as much words.

So, I guess, the way to go is the option 1.

Let's see...

x2=LOG2(X), 0 < X < 1, 0.16 notation
=====================================
Note that X should be greater than zero, because log2(0)=-Inf. Actually zero x is the easiest case, because 0^y=0. (y!=0)

In our case, log2(x) is in the range from -16 to -2.2e-5. So we can use 6.10 notation (to reserve one bit for multiplication by y) with implicit sign (it's always negative) for log result.

To find the integer part of log, we rotate x left until carry is set. The number of rotations is the approximated log result (will give us an interpolation point).

log2
        clrf x2hi
log2loop
        clrc
        rlf xl, f
        rlf xh, f
        incf x2hi, f
        skpc
         goto log2loop
;x2hi contains approximated log (integer part + error)
        clrc            ;normalize x2
        rlf x2hi, f     ;for 6.10 notation
        rlf x2hi, f     ;
        clrf x2lo       ;clear lower byte

Then we use 4 higher bits of the rotated x and get a fraction number from a table of log2(0.5)-log2(0.5:0.5/16:1), that is a table of values that we add to the already found integer part (note they are negative, so we can use subtraction instead).

 Columns 1 through 4
                  0  -0.08746284125034  -0.16992500144231  -0.24792751344359
  Columns 5 through 8 
  -0.32192809488736  -0.39231742277876  -0.45943161863730  -0.52356195605701
  Columns 9 through 12 
  -0.58496250072116  -0.64385618977472  -0.70043971814109  -0.75488750216347
  Columns 13 through 16 
  -0.80735492205760  -0.85798099512757  -0.90689059560852  -0.95419631038688
  Column 17 
  -1.00000000000000

Multiplied by -1024 and rounded (10 bits):
» round(-1024*(log2(0.5)-log2(0.5:0.5/16:1)))
ans =

  Columns 1 through 6 
           0          90         174         254         330         402
  Columns 7 through 12 
         470         536         599         659         717         773
  Columns 13 through 17 
         827         879         929         977        1024
(note, the table should have 16+1 entries!)
         
;prepare the index, i.e. shift xh right 3 times (4 bit index
;shifted left once)
        rrf xh, w
        movwf index
        rrf index, f
        rrf index, w
        andlw 0x1E
        movwf index
        call log2_table ;read the table entry to temphi:templo
;subtract it from current x2
        movf templo, w
        subwf x2lo, f
        movf temphi, w
        skpc
         incfsz temphi, w
          subwf x2hi, f
;read next point
        incf index, f
        incf index, f
        call log2_table ;read the table entry to temphi:templo
;find the difference with the previous point
        movf x2lo, w
        addwf templo, f
        movf x2hi, w
        skpnc
         incf x2hi, w
        addwf temphi, w
;leave only two lsb's of temphi (difference is 10 bit long)
        andlw 0x03
        movwf temphi
;now the difference is in temp
;multiply it by next 8 bits of the rotated x and divide by
;2^8=256

;shift xhi:xlo left by 4 bits to get 8 bits of the multiplier
;in xh. (xl is garbage)
        swapf   xhi, w
        andlw   0xF0
        movwf   xhi
        swapf   xlo, w
        andlw   0x0F
        iorwf   xhi, f
;we will simultaneously multiply and subtract from x2
        clrf xlo        ;use xlo as a temp
        movlw 8
        movwf index     ;use index as a counter
log2loop_mul
        rlf xhi, f      ;shift next multiplier bit to carry
        bnc log2loop_mul2 ;skip if no carry

        movf templo, w  ;subtract temp from x2hi:x2lo:xlo
        subwf xlo, f
        movf temphi, w
        skpc
         incfsz temphi, w
          subwf x2lo, f
        skpc
         decf x2hi, f

log2loop_mul2
        rrf x2hi, f     ;x2hi becomes garbage, but we'll
                        ;overwrite it
        rrf x2lo, f
        rrf xlo, f
        decfsz index, f
         goto log2loop_mul
;multiply x2 by 256 and overwrite x2hi to compensate 8 right
;rotations of x2 during the multiplication above
        movf x2lo, w
        movwf x2hi
        movf xlo, w
        movwf x2lo
;now x2 contains log2(x)!

I didn't debug that, but I hope it helps. Please consider publishing your results for us. :) Nikolai

Here is the finished, debugged, full routine:

;************************************************
;
; LOG2 routine (approximate, using a lookup table)
;
; Input: xhi:xlo - unsigned fractional number in the 0.16
;        format - 0 integer bits and 16 integer bits.
;        Note, the input value is destroyed.
;
; Output: x2hi:x2lo = -log2(xhi:xlo); the result is in
;        the 6.10 format (6 integer and 10 fractional bits).
;        Note that the logarithm of a fractional number
;        is negative, however the returned result is
;        positive. The negative sign is implicit.
;
; Temps:
;   temphi, templo, index
;
; 28 Sep 2000 Nikolai Golovchenko - initial version
; 17 Nov 2003 Nikolai Golovchenko - fixed a bug in the table and
;                                   optimized the routine.
;************************************************
log2

; First, find the number of cleared most significant
; bits in xhi:xlo. This will be the integer portion
; of the result. For example:
; log2([0.5 0.25 0.125 .. 1/65536]) = [-1 -2 -3 .. -16]
;
; If all bits are zero, return 16 (equivalent to
; replacing a zero input value with 0x0001).

        clrf x2hi
log2loop
        clrc
        rlf xlo, f
        rlf xhi, f
        incf x2hi, f
        btfsc x2hi, 4   ; limit the number of iterations to 16
         goto log2IntDone
        skpc
         goto log2loop
log2IntDone

; x2hi contains the integer portion of the result.
;
; xhi:xlo is normalized so that it's a fractional number
; between 0 and 1-1/32768 ~= 0.99997
;
; It's almost equivalent to the fractional portion of 
; the result, but we can minimize the error by using a small,
; 16-element look-up table with linear interpolation.
;

; Now normalize x2 so it matches the 6.10 notation of the
; result, and clear the lower 10 bits.
        clrc
        rlf x2hi, f
        rlf x2hi, f
        clrf x2lo

; Since the look-up table has 16 elements, use the higher
; 4 bits of xhi:xlo as an index to the first point. The
; the next element in the table will be used as the second point.
; The other 12 fractional bits in xhi:xlo can be used for linear 
; interpolation between the two points. Note that we
; actually need one extra element in the table when 
; xhi:xlo is has all ones in the higher 4 bits.
;

; load index with the xhi[7:4] to look up the first point.
        swapf xhi, w
        andlw 0x0F
        movwf index
        call log2_table 

; save the first point. To save memory, we can do
; it in place of x2hi:x2lo since we know the lower 10 
; bits are zero. 
        movf templo, w
        subwf x2lo, f
        movf temphi, w
        skpc
         incfsz temphi, w
          subwf x2hi, f

; look up the second point
        incf index, f
        call log2_table 

; find the difference between the second and the first
; points. temp = temp + x2.
; The difference will be in the lower 10 bits of temp.
        movf x2lo, w
        addwf templo, f
        movf x2hi, w
        skpnc
         incf x2hi, w
        addwf temphi, w
        andlw 0x03              ; ignore all integer bits
        movwf temphi

; Perform the linear interpolation. Basically,
; multiply the difference in temp by the lower
; 12 bits in xhi:xlow. We'll actually use only 8 bits.
;

; shift the 8-bit multiplier into xhi and clear xlo
        swapf   xhi, f
        swapf   xlo, w
        xorwf   xhi, w
        andlw   0x0F
        xorwf   xhi, f
        clrf xlo

; init the loop counter
        movlw 8
        movwf index     
log2loop_mul
        rrf xhi, f      
        bnc log2loop_mul2
        movf templo, w  
        subwf xlo, f
        movf temphi, w
        skpc
         incfsz temphi, w
          subwf x2lo, f
        skpc
         decf x2hi, f
log2loop_mul2
        rrf x2hi, f     
        rrf x2lo, f
        rrf xlo, f
        decfsz index, f
         goto log2loop_mul

; we are done with the linear interpolation:
;
; x2hi:x2lo:xlo = x2hi:x2lo * 256 - xhi * temphi:templo =
;               = 256 * (x2hi:x2lo - xhi / 256 * temphi:templo)

; now shift the lower two bytes up by 8 bits to 
; get the final result:
;
; x2hi:x2lo = x2hi:x2lo - xhi / 256 * temphi:templo
;
        movf x2lo, w
        movwf x2hi
        movf xlo, w
        movwf x2lo
        return

;
; log2_table
;
; Input: index - the number of the element to look up,
;               between 0 and 16. No range checking 
;               is performed.
; Output: temphi:templo = table[index]
;
log2_table
        ; look up the lower byte first
        clrc
        rlf index, w            ; w = index * 2
        call log2tableStart
        movwf templo

        ; look up the higher byte
        setc
        rlf index, w            ; w = index * 2 + 1
        call log2tableStart
        movwf temphi
        return

log2tableStart
        addwf PCL, f
        DT   0 & 0xFF,   0 >> 8,  90 & 0xFF,  90 >> 8
        DT 174 & 0xFF, 174 >> 8, 254 & 0xFF, 254 >> 8
        DT 330 & 0xFF, 330 >> 8, 402 & 0xFF, 402 >> 8
        DT 470 & 0xFF, 470 >> 8, 536 & 0xFF, 536 >> 8
        DT 599 & 0xFF, 599 >> 8, 659 & 0xFF, 659 >> 8
        DT 717 & 0xFF, 717 >> 8, 773 & 0xFF, 773 >> 8
        DT 827 & 0xFF, 827 >> 8, 879 & 0xFF, 879 >> 8
        DT 929 & 0xFF, 929 >> 8, 977 & 0xFF, 977 >> 8
        DT 1024 & 0xFF, 1024 >> 8
 IF (($ - 1) >> 8) - ((log2tableStart + 1) >> 8) != 0
        error 'log2 table crossed 8-bit boundary'
 ENDIF
;************************************************

See also:


file: /Techref/microchip/math/power/16lr-ng.htm, 11KB, , updated: 2004/4/6 20:52, local time: 2018/9/23 02:19, owner: NG--944,
TOP NEW HELP FIND: 
54.158.52.166:LOG IN

 ©2018 These pages are served without commercial sponsorship. (No popup ads, etc...).Bandwidth abuse increases hosting cost forcing sponsorship or shutdown. This server aggressively defends against automated copying for any reason including offline viewing, duplication, etc... Please respect this requirement and DO NOT RIP THIS SITE. Questions?
Please DO link to this page! Digg it! / MAKE! / 

<A HREF="http://www.piclist.com/techref/microchip/math/power/16lr-ng.htm"> PIC Microcontoller Math Method </A>

After you find an appropriate page, you are invited to your to this massmind site! (posts will be visible only to you before review) Just type in the box and press the Post button. (HTML welcomed, but not the <A tag: Instead, use the link box to link to another page. A tutorial is available Members can login to post directly, become page editors, and be credited for their posts.


Link? Put it here: 
if you want a response, please enter your email address: 
Attn spammers: All posts are reviewed before being made visible to anyone other than the poster.
Did you find what you needed?

  PICList 2018 contributors:
o List host: MIT, Site host massmind.org, Top posters @20180923 RussellMc, Van Horn, David, Sean Breheny, David C Brown, Bob Blick, Isaac M. Bavaresco, Neil, Denny Esterline, Harold Hallikainen, John Gardner,
* Page Editors: James Newton, David Cary, and YOU!
* Roman Black of Black Robotics donates from sales of Linistep stepper controller kits.
* Ashley Roll of Digital Nemesis donates from sales of RCL-1 RS232 to TTL converters.
* Monthly Subscribers: Gregg Rew. on-going support is MOST appreciated!
* Contributors: Richard Seriani, Sr.
 

Welcome to www.piclist.com!

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

  .