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
;************************************************

  .