Hi Dan,

Try the logarithm ruler approach.

   z =3D x^y
log2 z =3D log2 x^y
log2 z =3D y * log2 x

   z =3D 2^(y*log2(x))
   =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D

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 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=3DLOG2(X), 0 < X < 1, 0.16 notation
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
Note that X should be greater than zero, because log2(0)=3D-Inf.
Actually zero x is the easiest case, because 0^y=3D0. (y!=3D0)

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.24792751344=
359
  Columns 5 through 8=20
  -0.32192809488736  -0.39231742277876  -0.45943161863730  -0.52356195605=
701
  Columns 9 through 12=20
  -0.58496250072116  -0.64385618977472  -0.70043971814109  -0.75488750216=
347
  Columns 13 through 16=20
  -0.80735492205760  -0.85798099512757  -0.90689059560852  -0.95419631038=
688
  Column 17=20
  -1.00000000000000

Multiplied by -1024 and rounded (10 bits):
=BB round(-1024*(log2(0.5)-log2(0.5:0.5/16:1)))
ans =3D
  Columns 1 through 6=20
           0          90         174         254         330         402
  Columns 7 through 12=20
         470         536         599         659         717         773
  Columns 13 through 17=20
         827         879         929         977        1024
(note, the table should have 16+1 entries!)
        =20
;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=3D256

;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)!
       =20
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=20
By: D. Schouten .=20

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 =3D 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 mailto:listserv@mitvma.mit.edu?body=3DSET%20PICList%20DIGEST

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