by Nikolai Golovchenko debugged by Daniel Schouten [daniels at XS4ALL.NL]
Well, let's calculate y=2^x! x is Q6.10 notation (always negative with implicit sign). Output y has Q0.16 notation (always less then 1). Therefore x should not be zero (2^0 = 1).
y is in the range from 0.9993 (x = -1/1024) to 5e-20 (x = -64 + 1/1024). Okay....
y=2^x
y=2^-(32*x.15+16*x.14+...1*x.10 + 0.5*x.9+0.25*x.8+...)
\========================/ \===================/
x_int x_frac
x_int = 0:1:63
x_frac = 0:1/1024:(1-1/1024)
y=2^-(x_int+x_frac)
y=2^(-x_int) * 2^(-x_frac)
y= (2^(-x_frac)) >> x_int
=========================
So we first find 2^-x_frac using a table and then shift the result right x_int times.
First we should convert the fraction part through a table to a new fraction part. 2^(-0..-1) = 0.5..1
; y = 2^x2
;Input: x2hi:x2lo
;Output: yhi:ylo
;Form table index from 4 higher bits of x2 fraction part
rlf x2lo, f
rlf x2ho, w
movwf index
rlf x2lo, w
rlf index, f
;here middle 6 bits of x2lo contain multiplier for linear
;interpolation, and index has 4 lower bits of table pointer
rlf index, w ;multiply index by 2 for 16 bit
;entries addressing
andlw 0x1E ;and clear lsb and higher bits
movwf index
;Read first point to temphi:templo
call exp2_table
;clear result and subtract the first point from it
clrf yhi
clrf ylo
movf templo, w
subwf ylo, f
movf temphi, w
skpc
incfsz temphi, w
subwf yhi, f
;Read next point to temphi:templo
incf index, f
incf index, f
call exp2_table
;find difference with previous point to temphi:templo
movf ylo, w
addwf templo, f
movf yhi, w
skpnc
incfsz yhi, w
addwf temphi, f
;multiply the difference by the next 6 fraction bits in
;x2lo and divide by 64. Then subtract it from the current
;result in yhi:ylo
movlw 6
movwf index ;use index as a loop counter
rrf x2lo, f ;align the 6 bits to the lsb
clrf temp ;we will use another temp
;for mul purposes
rrf yhi, f ;rotate y right 2 times
rrf ylo, f ;to make easy multiplication.
rrf temp, f ;lower bits will go to temp
rrf yhi, f
rrf ylo, f
rrf temp, f
exp2loop
rrf x2lo, f ;get next multiplier bit
bnc exp2loopnext ;skip subtraction if carry is clear
movf templo, w
subwf temp, f
movf temphi, w
skpc
incfsz temphi, w
subwf ylo, f
skpc
decf yhi, f
exp2loopnext
rrf yhi, f
rrf ylo, f
rrf temp, f
decfsz index, f
goto exp2loop
;shift the result 8 bits left
movf ylo, w
movwf yhi
movf temp, w
movwf ylo
;Now we will shift y right the number of times
;equal to the integer part of x2
clrc ;align integer bits
rrf x2hi, f ;to lsb
clrc
rrf x2hi, f
test x2hi
skpnz
return ;done
exp2loop2
clrc
rrf y2hi, f
rrf y2lo, f
decfsz x2hi, f
goto exp2loop2
return ;done
Now we need to calculate the table...
1-2.^(0:-1/16:-1)
ans =
Columns 1 through 4
0 0.04239671930143 0.08299595679533 0.12187391981335
Columns 5 through 8
0.15910358474629 0.19475483402537 0.22889458729603 0.26158692703025
Columns 9 through 12
0.29289321881345 0.32287222653155 0.35158022267450 0.37907109396326
Columns 13 through 16
0.40539644249864 0.43060568262165 0.45474613366737 0.47786310878629
Column 17
0.50000000000000
Multiplied by 65536 and rounded
round(65536*(1-2.^(0:-1/16:-1)))
ans =
Columns 1 through 6
0 2779 5439 7987 10427 12763
Columns 7 through 12
15001 17143 19195 21160 23041 24843
Columns 13 through 17
26568 28220 29802 31317 32768
»
So, the table would be:
exp2_table
movf index, w
call exp2tableStart
movwf templo
incf index, w
call exp2tableStart
movwf temphi
return
exp2tableStart
addwf PCL, f
DT 0 & 0xFF, 0 >> 8, 2779 & 0xFF, 2779 >> 8
DT 5439 & 0xFF, 5439 >> 8, 7987 & 0xFF, 7987 >> 8
DT 10427 & 0xFF, 10427 >> 8, 12763 & 0xFF, 12763 >> 8
DT 15001 & 0xFF, 15001 >> 8, 17143 & 0xFF, 17143 >> 8
DT 19195 & 0xFF, 19195 >> 8, 21160 & 0xFF, 21160 >> 8
DT 23041 & 0xFF, 23041 >> 8, 24843 & 0xFF, 24843 >> 8
DT 26568 & 0xFF, 26568 >> 8, 28220 & 0xFF, 28220 >> 8
DT 29802 & 0xFF, 29802 >> 8, 31317 & 0xFF, 31317 >> 8
DT 32768 & 0xFF, 32768 >> 8
IF (($ - 1) >> 8) - ((exp2tableStart + 1) >> 8) != 0
error 'exp2 table crossed 8-bit boundary'
ENDIF
And the cleaned up version is below... Yeah, this one was harder and not really the reverse of log!
Don't forget to check x and y for zero in x^y.
#include <p16f84.inc>
cblock 0x20
xhi, xlo
templo, temphi, temp
index
x2hi, x2lo
endc
#define tval 65300
movlw tval & 0xFF
movwf xlo
movlw tval >> 8
movwf xhi
call log2
nop
call exp2
nop
;************************************************
; EXP2 routine for fixed point unsigned
; 16 bit values in 6.10 notation (implicit minus)
;
; Input: x2hi:x2lo unsigned Q6.10 (implicit minus)
; Output: xhi:xlo unsigned Q0.16
; Temporary: index, templo, temphi, temp
;
; Maximum error: 0.025% (for x2 < 1)
;
; Size: 121 words
; Time: min 2+8+18*2+10+16+10*6-1+6+2=139
; max 2+8+18*2+10+16+17*6-1+6+1+14+6*16-1+2=291
;
; Note: Zero input is illegal! Will result in
; zero output. (2^0 = 1, but x is always less
; than 1)
;
; 30 Sep 2000 by Nikolai Golovchenko
; 02 Oct 2000 debugged by Daniel Schouten [daniels at XS4ALL.NL]
;************************************************
exp2
;Form table index from 4 higher bits of x2 fraction part
rlf x2lo, f
rlf x2hi, w
movwf index
rlf x2lo, w
rlf index, f
;here middle 6 bits of x2lo contain multiplier for linear
;interpolation, and index has 4 lower bits of table pointer
rlf index, w ;multiply index by 2 for 16 bit
;entries addressing
andlw 0x1E ;and clear lsb and higher bits
movwf index
;Read first point to temphi:templo
call exp2_table
;clear result and subtract the first point from it
clrf xhi
clrf xlo
movf templo, w
subwf xlo, f
movf temphi, w
skpc
incfsz temphi, w
subwf xhi, f
;Read next point to temphi:templo
incf index, f
incf index, f
call exp2_table
;find difference with previous point to temphi:templo
movf xlo, w
addwf templo, f
movf xhi, w
skpnc
incfsz xhi, w
addwf temphi, f
;multiply the difference by the next 6 fraction bits in
;x2lo and divide by 64. Then subtract it from the current
;result in xhi:xlo
movlw 6
movwf index ;use index as a loop counter
rrf x2lo, f ;align the 6 bits to the lsb
clrf temp ;we will use another temp
;for mul purposes
rrf xhi, f ;rotate x right 2 times
rrf xlo, f ;to make easy multiplication.
rrf temp, f ;lower bits will go to temp
rrf xhi, f
rrf xlo, f
rrf temp, f
exp2loop
rrf x2lo, f ;get next multiplier bit
bnc exp2loopnext ;skip subtraction if carry is clear
movf templo, w
subwf temp, f
movf temphi, w
skpc
incfsz temphi, w
subwf xlo, f
skpc
decf xhi, f
exp2loopnext
rrf xhi, f
rrf xlo, f
rrf temp, f
decfsz index, f
goto exp2loop
;shift the result 8 bits left
movf xlo, w
movwf xhi
movf temp, w
movwf xlo
;Now we will shift x right the number of times
;equal to the integer part of x2
clrc ;align integer bits
rrf x2hi, f ;to lsb
clrc
rrf x2hi, f
tstf x2hi
skpnz
return ;done
movlw 17 ;limit shifts number to 16 maximum
subwf x2hi, w
movlw 16
skpnc
movwf x2hi
;if fractional part is zero than set carry, else
;clear carry
movf xlo, w
iorwf xhi, w
setc
skpz
clrc
;limit number of si
;perform shifts
exp2loop2
rrf xhi, f
rrf xlo, f
clrc
decfsz x2hi, f
goto exp2loop2
return ;done
exp2_table
movf index, w
call exp2tableStart
movwf templo
incf index, w
call exp2tableStart
movwf temphi
return
exp2tableStart
addwf PCL, f
DT 0 & 0xFF, 0 >> 8, 2779 & 0xFF, 2779 >> 8
DT 5439 & 0xFF, 5439 >> 8, 7987 & 0xFF, 7987 >> 8
DT 10427 & 0xFF, 10427 >> 8, 12763 & 0xFF, 12763 >> 8
DT 15001 & 0xFF, 15001 >> 8, 17143 & 0xFF, 17143 >> 8
DT 19195 & 0xFF, 19195 >> 8, 21160 & 0xFF, 21160 >> 8
DT 23041 & 0xFF, 23041 >> 8, 24843 & 0xFF, 24843 >> 8
DT 26568 & 0xFF, 26568 >> 8, 28220 & 0xFF, 28220 >> 8
DT 29802 & 0xFF, 29802 >> 8, 31317 & 0xFF, 31317 >> 8
DT 32768 & 0xFF, 32768 >> 8
IF (($ - 1) >> 8) - ((exp2tableStart + 1) >> 8) != 0
error 'exp2 table crossed 8-bit boundary'
ENDIF
;************************************************