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