part 1 6460 bytes content-type:text/plain; charset=ISO-8859-1 (decoded quoted-printable)
Hi Daniel,
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 attached... 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.
Hope it works for you :)
Nikolai
---- Original Message ----
From: D. Schouten <spam_OUTdanielsTakeThisOuT
xs4all.nl>
Sent: Saturday, September 30, 2000 0:56:03
To: Nikolai Golovchenko
Subj: [PIC]: [MATH]: Power function, part 2
{Quote hidden}> Hi List / Nikolai,
> Thanks again for the log2() function Nikolai.
> I've implemented it in my program and it works perfectly now (took a
> little time to get it work because of a really stupid mistake of me).
> I have also implemented the multiplication by y (compared to the log2
> function, this is a piece of cake of course), so the the
> power-function is now finished till :
> y * log2(x)
> where the result is in an unsigned Q6.10 (implicit minus) notation.
> Now I only have to take care of the exp() function. But the problem is
> that I'm not a math-wizard like Nikolai. So I would really appreciate
> it if you could help me out with that routine also. The output of the
> exp() function has to be in the unsigned Q0.16 format again.
> Again thanks a lot.
> Daniel...
part 2 8615 bytes content-type:application/octet-stream; name="log2exp2.asm" (decode)
part 3 144 bytes
--
http://www.piclist.com hint: The list server can filter out subtopics
(like ads or off topics) for you. See http://www.piclist.com/#topics
