Searching \ for '[PIC]: [MATH] x^y routine implementation' in subject line. ()
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.
'[PIC]: [MATH] x^y routine implementation'
2000\09\27@130355 by

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 listservmitvma.mit.edu?body=SET%20PICList%20DIGEST

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

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
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
movf x2hi, w
skpnc
incf x2hi, 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 listservmitvma.mit.edu?body=SET%20PICList%20DIGEST

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

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
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
movf x2hi, w
skpnc
incf x2hi, 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
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-requestmitvma.mit.edu

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
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
movf x2hi, w
skpnc
incf x2hi, 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
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-requestmitvma.mit.edu

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!

--

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}
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.

Scott Dattalo wrote:
{Quote hidden}

--

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: rallcfltd.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"
================================================================

--