Post by thread:[PIC]: [MATH] x^y routine implementation

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 .....listservKILLspam@spam@mitvma.mit.edu?body=SET%20PICList%20DIGEST -- http://www.piclist.com hint: To leave the PICList piclist-unsubscribe-requestKILLspammitvma.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 ;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-requestKILLspam.....mitvma.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 ;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_OUTTakeThisOuTmitvma.mit.edu

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}> > 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 -- http://www.piclist.com hint: PICList Posts must start with ONE topic: "[PIC]:" PIC only "[EE]:" engineering "[OT]:" off topic "[AD]:" ad's