Searching \ for '[PIC]: [MATH]: Power function, part 2' in subject line. ()
Make payments with PayPal - it's fast, free and secure! Help us get a faster server
FAQ page:
Search entire site for: '[MATH]: Power function, part 2'.

Exact match. Not showing close matches.
PICList Thread
'[PIC]: [MATH]: Power function, part 2'
2000\09\29@175856 by D. Schouten

picon face
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.


-- hint: PICList Posts must start with ONE topic:
"[PIC]:" PIC only "[EE]:" engineering "[OT]:" off topic "[AD]:" ad's

2000\09\30@043514 by Nikolai Golovchenko

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 =

y is in the range from 0.9993 (x = -1/1024) to 5e-20 (x
= -64 + 1/1024). Okay....


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) * 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
        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
        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
       rrf x2lo, f      ;get next multiplier bit
       bnc exp2loopnext ;skip subtraction if carry is clear

       movf templo, w
       subwf temp, f
       movf temphi, w
        incfsz temphi, w
         subwf ylo, f
        decf yhi, f


       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
       rrf x2hi, f

       test x2hi
        return         ;done
       rrf y2hi, f
       rrf y2lo, f
       decfsz x2hi, f
        goto exp2loop2
       return          ;done

Now we need to calculate the table...
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
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:

       movf index, w
       call exp2tableStart
       movwf templo
       incf index, w
       call exp2tableStart
       movwf temphi

       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'

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 :)

---- Original Message ----
From: D. Schouten <>
Sent: Saturday, September 30, 2000 0:56:03
 To: Nikolai Golovchenko
Subj: [PIC]: [MATH]: Power function, part 2

{Quote hidden}

part 2 8615 bytes content-type:application/octet-stream; name="log2exp2.asm" (decode)

part 3 144 bytes
-- hint: The list server can filter out subtopics
(like ads or off topics) for you. See

'[PIC]: [MATH]: Power function, part 2'
2000\10\03@131101 by D. Schouten
picon face
Hi List / Nikolai,

Thanks for the Exp2 routine Nikolai!

I've also implemented this last piece of the power function. It works
perfectly now but it took a little time to debug. I found two little
bugs in the Exp2 routine in the log2exp2.asm attachment:

1.) the code :

       tstf x2hi
        return         ;done

has to be placed UNDER :

       clrc            ;align integer bits
       rrf x2hi, f     ;to lsb
       rrf x2hi, f

instead of above. This is done right in the explanation posted on the
list but wrong in the attached log2exp2.asm.

2.) in the piece of code :

       movlw 17     ;limit shifts number to 16 maximum
       subwf x2hi, w
       movlw 16
        movwf x2hi

the instruction :


has to be changed in :


to limit shifts number to 16 like Nikolai says in the comment.

Again, thank you very much Nikolai. The routine is not only very
accurate, it's also efficient. The whole power function takes less
then 1000 cycles to calculate and uses less then 300 words of
codespace. WOW!

Again thanks!



-- hint: To leave the PICList

More... (looser matching)
- Last day of these posts
- In 2000 , 2001 only
- Today
- New search...