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 byteThen 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.00000000000000Multiplied 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: 2020/10/30 05:47,
owner: NG--944,
34.201.9.19:LOG IN |
| ©2020 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? <A HREF="http://www.piclist.com/techref/microchip/math/power/16lr-ng.htm"> PIC Microcontoller Math Method </A> |
| Did you find what you needed? |
|
PICList 2020 contributors:
o List host: MIT, Site host massmind.org, Top posters @20201030 Bob Blick, RussellMc, madscientistatlarge, Jason White, Neil, David C Brown, Sean Breheny, Harold Hallikainen, James Cameron, Allen Mulvey, * 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! |
.