# PICMicrocontollerMathMethod

## 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
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
movf x2hi, w
skpnc
incf x2hi, 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
movf x2hi, w
skpnc
incf x2hi, 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
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
;************************************************

```

 file: /Techref/microchip/math/power/16lr-ng.htm, 11KB, , updated: 2004/4/6 19:52, local time: 2018/12/10 12:56, owner: NG--944, TOP NEW HELP FIND:  54.221.147.93: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! /  PIC Microcontoller Math Method

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.

Attn spammers: All posts are reviewed before being made visible to anyone other than the poster.
 Did you find what you needed? "No. I'm looking for: " "No. Take me to the search page." "No. Take me to the top so I can drill down by catagory" "No. I'm willing to pay for help, please refer me to a qualified consultant" "No. But I'm interested. me at when this page is expanded."

 PICList 2018 contributors: o List host: MIT, Site host massmind.org, Top posters @20181210 RussellMc, Van Horn, David, Sean Breheny, David C Brown, Neil, Isaac M. Bavaresco, Bob Blick, Harold Hallikainen, AB Pearce - UKRI STFC, 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.

.