;Variables:
Square:4	;32 bits source (=0 on return)
Root:2		;16 bits result	(may get updated after rounding)
Remain:2	;16 bits remainder (becomes meaningless after rounding)
Counter		;loop counter (=0 on return)

;Routine to square 32 bit integer with.optional rounding to nearest root (by Antonie Daas)
; Input: Square:4,  Output: Root:2, Rmndr:2,  Local: Counter
; 9 data registers including input and output,
; 45 instructions (+21 for rounding), 631 cycles max. incl. return (+ 19 for rounding)
Sqroot	clrf	Root		;clear initial root
	clrf	Root+1
	clrf	Remain		;clear initial remainder
	clrf	Remain+1
	movlw	16		;set number of bits to process
	movwf	Counter		;set bit counter
;Test for next bit of Root
sqrt_lp	movlw	0x40		;load test mask
	subwf	Square+3,W	;compare mask with msb of Square
	movf	Root,W		;get low byte of Root so far
	skpc			;if borrow, then:
	 incfsz	Root,W		;add borrow bit
	 subwf	Remain,W	;compare Root low so far with remainder low so far
	movf	Root+1,W	;get high byte of Root so far
	skpc			;if borrow, then:
	 incfsz	Root+1,W	;add borrow bit
	 subwf	Remain+1,W	;compare Root high so far with remainder high so far
	skpc			;if borrow (C=0), then:
	 goto	sqrt_rl		;skip subtraction and append root with 0
;Subtract root so far from Sumsqr and Rmdr
	movlw	0x40		;load test mask again
	subwf	Square+3,F	;subtract mask from msb Square
	movf	Root,W		;get low byte of root so far
	skpc			;if borrow, then:
	 incfsz	Root,W		;add borrow bit
	 subwf	Remain,F	;subtract root low so far from remainder low so far
	movf	Root+1,W	;get high byte of root so far
	skpc			;if borrow, then:
	 incfsz	Root+1,W	;add borrow bit
	 subwf	Remain+1,F	;subtract root high so far from remainder high so far (C=1)
sqrt_rl	rlf	Root,F		;append new root bit
	rlf	Root+1,F	;and propagate previous bits, C=0
	rlf	Square,F	;shift Square up 1 bit
	rlf	Square+1,F	;and second byte,
	rlf	Square+2,F	;and third byte
	rlf	Square+3,F	;and fourth byte, shift msb into Remainder
	rlf	Remain,F	;append remainder
	rlf	Remain+1,F	;shift up other bits, C=0
	rlf	Square,F	;shift Square up 2nd bit
	rlf	Square+1,F	;and second byte
	rlf	Square+2,F	;and third byte
	rlf	Square+3,F	;and fourth byte, shift into Remainder
	rlf	Remain,F	;append remainder second bit
	rlf	Remain+1,F	;shift high byte up, C=0 (except for last itteration)
	decfsz	Counter,F	;if not all 16 bits done, then:
	 goto	sqrt_lp		;repeat procedure
;Round result to nearest integer
	skpnc			;if Carry shifted out of Remain, then:
	 goto	sqrt_ru		;skip test as remainder > root
	movlw	0x40		;load test mask
	subwf	Square+3,W	;compare mask with msb of Square fraction
	movf	Root,W		;get low byte of Root so far
	skpc			;if borrow, then:
	 incfsz	Root,W		;add borrow bit
	 subwf	Remain,W	;compare Root low so far with remainder low so far
	movf	Root+1,W	;get high byte of Root so far
	skpc			;if borrow, then:
	 incfsz	Root+1,W	;add borrow bit
	 subwf	Remain+1,W	;compare Root high so far with remainder high so far
	skpc			;if borrow, then:
	 return			;done, keep rounded down
sqrt_ru	incf	Root,F		;round up
	skpnz			;if overflow, then:
	 incfsz	Root+1,F	;propagate carry, if result is not 0, then:
	 return			;done, if result 0, then:
	decf	Root,F		;undo rounding up
	decf	Root+1,F	;and its high byte
	return