PIC Microcontoller Math Method for square root

16 bit integer square root

From: Scott Dattalo

I know at least a few people are interested in this. I finally got around to applying the 18cxxx instruction set to the sqrt algorithm (see: http://www.dattalo.com/technical/software/pic/picsqrt.html and http://www.dattalo.com/technical/theory/sqrt.html ).

It takes only 25 instructions and executes between 68 and 108 (excluding the call). I have no idea what the average is. Although, if I made the test code isochronous, it'd be easy to find that out. I still think there's room for optimization too (but after spending 4 hours fixing an elusive bug in gpsim, I'm not up to it [but I see two places that can be tweaked]). BTW, gpsim runs through all 2^16 test cases in just 2 seconds!

Scott

       list    p=18c242,t=ON,c=132,n=80
        radix   dec


include "p18c242.inc"

  cblock  0
        N_hi,N_lo
        mask
        x,y
        s_hi,s_lo
  endc


        ORG     0               ;Reset Vector


Main

        CLRF    status
        CLRF    x
        CLRF    s_hi
        CLRF    s_lo
        INCF    x,f

; Here's some test code that goes through all 2^16 cases
;
; psuedo code:
:
;  int  s=0;
;  int  x=0;
;
;  do {
;
;    if(x*x == s)
;      x++;
;
;    n = s;
;    w = sqrt(n);
;    if(w+1 != x)
;      failed();
;
;    s++;
;
;  } while (x<256);

;  for(i=0; i< 256; i++)
xxx
        MOVF    x,W             ;W = x*x
        MULWF   x

        movf    prodh,w,0       ;if x*x is equal to the current
        cpfseq  s_hi            ;test case, s, then we need to advance
         bra    L1              ;x. Note that x will always be 1 greater
        movf    prodl,w,0       ;than the sqrt(s)
        cpfseq  s_lo
         bra    L1

        infsnz  x,f
here     bra    here

L1:
  ; N = s
        movf    s_hi,w
        movwf   N_hi
        movf    s_lo,w
        movwf   N_lo


  ; find the sqrt of N (N_hi:N_lo)
        RCALL    sqrt

        infsnz  s_lo,f          ;s++ advance for the next test case
         incf   s_hi,f

        addlw   1               ;if w (sqrt(N)) is one less than x
        cpfseq  x               ; then we've passed
         bra    failed


        bra     xxx

failed:
        bra     xxx             ;Set a break point here to catch failures




;----------------------------------------------------------
;sqrt
;
; The purpose of this routine is take the square root of a 16 bit
;unsigned integer.
;Inputs:  N_hi - High byte of the 16 bit number to be square rooted
;         N_lo - Low byte     "                  "             "
;Outputs: W register returned with the 8 bit result
;
;Memory used:
;  mask
;
;Minimum 68 cycles
;Maximum 108 cycles

sqrt
        MOVLW   0xc0            ;Initialize value for mask
        MOVWF   mask
        MOVLW   0x40            ;Initial value for the root

sq1     CPFSLT   N_hi
         BRA     sq6
          ;Subtract the root developed so far

sq3     RLCF    N_lo,F          ;Shift N left one position
        RLCF    N_hi,F
        BC      sq4

        RRCF    mask,F          ;mov the 2-bit mask down 1
        XORWF   mask,W          ;clear the last and set the next bit

        BNC     sq1

   ;We are almost done. In the last iteration, we have 7 bits of the root. When
   ;"01" is appended to it, we will have a 9-bit number that must be subtracted
   ;from N.

        SUBWF   N_hi,F          ;
        SKPC                    ;If the upper 7 bits cause a borrow, then
          RETURN                ;the appended "01" will as well: We're done.

        SKPNZ                   ;If the result of the subtraction is zero
          BTFSC N_lo,7          ;AND the msb of N_lo is set then the LSB of the
                                ;root is zero.
           XORLW  1             ;Otherwise, it is one.

        RETURN                  ;

sq4     BTFSC   mask,0          ;Need to unconditionally set the current bit of
the root.
          RETURN                ;However, if we're through iterating, then
leave.

        RRNCF   mask,F
        XORWF   mask,W          ;Append "01" to the root developed so far.
sq6     SUBWF   N_hi,F
        IORWF   mask,W          ;Set the current bit
        bra     sq3             ;Go unconditionally set the current bit.



        END

(gpsim is a full-featured software simulator for Microchip PIC microcontrollers distributed under the GNU General Public License. http://dattalo.com/gnupic/gpsim.html )

Archive:

• PICList post "square roots for the 18F"

  .