Searching \ for '[PIC]: [MATH] x^y routine implementation' in subject line. ()
Make payments with PayPal - it's fast, free and secure! Help us get a faster server
FAQ page: www.piclist.com/techref/microchip/math/index.htm?key=math
Search entire site for: '[MATH] x^y routine implementation'.

Exact match. Not showing close matches.
PICList Thread
'[PIC]: [MATH] x^y routine implementation'
2000\09\27@130355 by D. Schouten

picon face
Hi All,


I'm in desparate need of a routine that calculates :

X raised to the power of Y (like the pow(x,y) function in C++)

where 0 < X < 1  (in a 0.16 fixed point notation)

and Y = 1.00 to 1.50 (in a 1.7 fixed point notation)

for a PIC16C73B.

I have about 2K of code-space and about 5000 clock-cycles available. I
already have implemented 8*8, 16*16, 24*16, 16/16, 32/16, 32-bit shift
left and 32-bit shift right routines (all unsigned) routines that can
be used if necessary.

I can't find such a routine on the web and the microchip application
note doesn't bring me any further since they don't have the
implementation for the 16Cxx series and besides that, they use 32-bit
precision wich is too high for my application.

I'm looking forward to any reply.

Thanks!

Daniel...

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
use spam_OUTlistservTakeThisOuTspammitvma.mit.edu?body=SET%20PICList%20DIGEST


2000\09\28@041709 by Nikolai Golovchenko

flavicon
face
Hi Dan,

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 http://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 byte

Then 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.00000000000000

Multiplied 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




Post piclist\2000\09\27\130355a
[PIC]: [MATH] x^y routine implementation By: D. Schouten .
Hi All,


I'm in desparate need of a routine that calculates :

X raised to the power of Y (like the pow(x,y) function in C++)

where 0 < X < 1  (in a 0.16 fixed point notation)

and Y = 1.00 to 1.50 (in a 1.7 fixed point notation)

for a PIC16C73B.

I have about 2K of code-space and about 5000 clock-cycles available. I
already have implemented 8*8, 16*16, 24*16, 16/16, 32/16, 32-bit shift
left and 32-bit shift right routines (all unsigned) routines that can
be used if necessary.

I can't find such a routine on the web and the microchip application
note doesn't bring me any further since they don't have the
implementation for the 16Cxx series and besides that, they use 32-bit
precision wich is too high for my application.

I'm looking forward to any reply.

Thanks!

Daniel...

--
http://www.piclist.com#nomail Going offline? Don't AutoReply us!
use .....listservKILLspamspam@spam@mitvma.mit.edu?body=SET%20PICList%20DIGEST

--
http://www.piclist.com hint: To leave the PICList
piclist-unsubscribe-requestspamKILLspammitvma.mit.edu


2000\09\28@075821 by Nikolai Golovchenko

flavicon
face
I found some time to debug the log2 routine. It actually
works! Error is very small, less than 0.05%.

The whole x^y thing should take about 319*2+160 ~= 800
cycles and 105*2+20=220 words. Not so bad :)

 Nikolai

;************************************************
; LOG2 routine for fixed point unsigned
; 16 bit values in 0.16 notation
;
; Input: xhi:xlo unsigned Q0.16 (modified)
; Output: x2hi:x2lo unsigned Q6.10 (implicit minus)
; Temporary: index, templo, temphi
;
; Maximum error: 0.05%
;
; Size: 105 words
; Time: min 2+1+7*1 -1+10+18*2+8+17+10*8-1+4+2=165
;       max 2+1+7*15-1+10+18*2+8+17+17*8-1+4+2=319
;
; Note: Zero input is illegal! Will result in
;       infinite loop.
;
; 28 Sep 2000 by Nikolai Golovchenko
;************************************************
log2
       clrf x2hi
log2loop
       clrc
       rlf xlo, f
       rlf xhi, 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 byte

;prepare the index, i.e. shift xh right 3 times (4 bit index
;shifted left once)
       rrf xhi, 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. (xlo 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
       rrf 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)!
       return

log2_table
       movf index, w
       call log2tableStart
       movwf templo
       incf index, w
       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
IF (($ - 1) >> 8) - ((log2tableStart + 1) >> 8) != 0
       error 'log2 table crossed 8-bit boundary'
ENDIF
;************************************************

--
http://www.piclist.com hint: To leave the PICList
.....piclist-unsubscribe-requestKILLspamspam.....mitvma.mit.edu


2000\09\28@080926 by Nikolai Golovchenko

flavicon
face
Oh s**t, the table should have 17 entries! Or a roll-over...
Corrected routine below.



----------------------------------------------------

I found some time to debug the log2 routine. It actually
works! Error is very small, less than 0.05%.

The whole x^y thing should take about 319*2+160 ~= 800
cycles and 105*2+20=220 words. Not so bad :)

 Nikolai

;************************************************
; LOG2 routine for fixed point unsigned
; 16 bit values in 0.16 notation
;
; Input: xhi:xlo unsigned Q0.16 (modified)
; Output: x2hi:x2lo unsigned Q6.10 (implicit minus)
; Temporary: index, templo, temphi
;
; Maximum error: 0.05%
;
; Size: 105 words
; Time: min 2+1+7*1 -1+10+20*2+8+17+10*8-1+4+2=169
;       max 2+1+7*15-1+10+20*2+8+17+17*8-1+4+2=323
;
; Note: Zero input is illegal! Will result in
;       infinite loop.
;
; 28 Sep 2000 by Nikolai Golovchenko
;************************************************
log2
       clrf x2hi
log2loop
       clrc
       rlf xlo, f
       rlf xhi, 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 byte

;prepare the index, i.e. shift xh right 3 times (4 bit index
;shifted left once)
       rrf xhi, 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. (xlo 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
       rrf 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)!
       return

log2_table
       movf index, w
       andlw 0x1F
       call log2tableStart
       movwf templo
       incf index, w
       andlw 0x1F
       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
IF (($ - 1) >> 8) - ((log2tableStart + 1) >> 8) != 0
       error 'log2 table crossed 8-bit boundary'
ENDIF
;************************************************

--
http://www.piclist.com hint: To leave the PICList
EraseMEpiclist-unsubscribe-requestspam_OUTspamTakeThisOuTmitvma.mit.edu


2000\09\29@003352 by Scott Dattalo
face
flavicon
face
On Wed, 27 Sep 2000, D. Schouten wrote:

> Hi All,
>
>
> I'm in desparate need of a routine that calculates :
>
>  X raised to the power of Y (like the pow(x,y) function in C++)
>
>  where 0 < X < 1  (in a 0.16 fixed point notation)
>
>  and Y = 1.00 to 1.50 (in a 1.7 fixed point notation)
>
> for a PIC16C73B.


In addition to the methods mentioned by Nikolai, there is also Feynman's Power
Algorithm. I haven't implemented this with pic code, but did test it out in
basic (a long time ago). The algorithm is described in a homework problem in
chapter 1 of Knuth's "The Art of Computer Programming: Fundamental Algorithms".
BTW, if you don't have the three volumes of Knuth's books, then get them!

--
http://www.piclist.com hint: PICList Posts must start with ONE topic:
"[PIC]:" PIC only "[EE]:" engineering "[OT]:" off topic "[AD]:" ad's


2000\09\29@085829 by daan xs4all

picon face
Scott,

thanks for the Knuth's books tip! I will check these out.

I'm currently working on the algorithm. I will keep you
up to date when I've made any progress.

Daniel...


{Original Message removed}

2000\09\29@100606 by M. Adam Davis

flavicon
face
I wish I could find a cheap source for them!  (or the original TEX files...  I
think he should release all of them to the public domain, though I suspect some
publisher still holds the rights...)

He is still, apparently, working on the fourth volume though, so when(if) it
comes out I imagine there will be a reprinting of all of them.

My local school library has the first one, but not the other two  <sigh>, but
occasionally you'll see them on ebay and at other places.

-Adam

Scott Dattalo wrote:
{Quote hidden}

--
http://www.piclist.com hint: PICList Posts must start with ONE topic:
"[PIC]:" PIC only "[EE]:" engineering "[OT]:" off topic "[AD]:" ad's


2000\09\29@100809 by Robert A. LaBudde

flavicon
face
At 02:55 PM 9/29/00 +0100, Daniel wrote:
>Scott,
>
>thanks for the Knuth's books tip! I will check these out.
>
>I'm currently working on the algorithm. I will keep you
>up to date when I've made any progress.

The standard method for x^y (x, y real) in high-level languages is via
scaling to a number between 0 and 1, then using a rational Chebyshev
expansions to calculate exp(a) and log(b). (x^y = e^(y ln x)

If x and y are integer, or y is integer, much simpler methods exist.

================================================================
Robert A. LaBudde, PhD, PAS, Dpl. ACAFS  e-mail: ralspamspam_OUTlcfltd.com
Least Cost Formulations, Ltd.            URL: http://lcfltd.com/
824 Timberlake Drive                     Tel: 757-467-0954
Virginia Beach, VA 23464-3239            Fax: 757-467-2947

"Vere scire est per causas scire"
================================================================

--
http://www.piclist.com hint: PICList Posts must start with ONE topic:
"[PIC]:" PIC only "[EE]:" engineering "[OT]:" off topic "[AD]:" ad's


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