Faster Square Roots - Printable Version

Code:

FUNCTION FASTCALL SQRT (radicand as FLOAT) as FLOAT
ASM
; FLOATS arrive in A ED CB
;A is the exponent.
          
          AND   A               ; Test for zero argument 
          RET   Z               ; Return with zero.
          
          ;Strictly we should test the number for being negative and quit if it is.
          ;But let's assume we like imaginary numbers, hmm?
          ; If you'd rather break it change to a jump to an error below.
          ;BIT   7,(HL)          ; Test the bit.
          ;JR    NZ,REPORT       ; back to REPORT_A 
                                ; 'Invalid argument'
          RES 7,E               ; Now it's a positive number, no matter what.
          
          call __FPSTACK_PUSH   ; Okay, We put it on the calc stack. Stack contains ABS(x)
                    
          ;   Halve the exponent to achieve a good guess.(accurate with .25 16 64 etc.)

                                ; Remember, A is the exponent.
          XOR   $80             ; toggle sign of exponent
          SRA   A               ; shift right, bit 7 unchanged.
          INC   A               ;
          JR    Z,ASIS          ; forward with say .25 -> .5
          JP    P,ASIS          ; leave increment if value > .5
          DEC   A               ; restore to shift only.
ASIS:     XOR   $80             ; restore sign.
          
          call __FPSTACK_PUSH   ; Okay, NOW we put the guess on the stack
          rst  28h	; ROM CALC    ;;guess,x
          DEFB $C3              ;;st-mem-3              x,guess
          DEFB $02              ;;delete                x 

SLOOP:    DEFB  $31             ;;duplicate             x,x.
          DEFB  $E3             ;;get-mem-3             x,x,guess
          DEFB  $C4             ;;st-mem-4              x,x,guess
          DEFB  $05             ;;div                   x,x/guess.
          DEFB  $E3             ;;get-mem-3             x,x/guess,guess          
          DEFB  $0F             ;;addition              x,x/guess+guess          
          DEFB  $A2             ;;stk-half              x,x/guess+guess,.5
          DEFB  $04             ;;multiply              x,(x/guess+guess)*.5
          DEFB  $C3             ;;st-mem-3              x,newguess
          DEFB  $E4             ;;get-mem-4             x,newguess,oldguess
          DEFB  $03             ;;subtract              x,newguess-oldguess
          DEFB  $2A             ;;abs                   x,difference.
          DEFB  $37             ;;greater-0             x,(0/1).
          DEFB  $00             ;;jump-true             x.

          DEFB  SLOOP - $       ;;to sloop              x.

          DEFB  $02             ;;delete                .
          DEFB  $E3             ;;get-mem-3             retrieve final guess.
          DEFB  $38             ;;end-calc              sqr x.

          jp __FPSTACK_POP
          
END ASM
END FUNCTION


FUNCTION FASTCALL SQRT16(radicand as uInteger) as uByte
asm
    XOR A
    AND A
    
    ld  a,l
    ld  l,h
    ld	de,0040h	; 40h appends "01" to D
	ld	h,d

    ld b,7
sqrt16loop:
	sbc	hl,de		; IF speed is critical, and you don't mind spending the extra bytes, you could unroll this loop 7 times instead of DJNZ.
	jr	nc,$+3	
	add	hl,de		
	ccf			
	rl	d		
	rla			
	adc	hl,hl		
	rla			
	adc	hl,hl		
    DJNZ sqrt16loop
	
    sbc	hl,de		; optimised last iteration
	ccf
	rl	d

    ld a,d
end asm
END FUNCTION


FUNCTION t AS ULONG
   RETURN INT((65536 * PEEK (23674) + 256 * PEEK(23673) + PEEK (23672)))
END FUNCTION

CLS

DIM a,b as float
DIM i as uInteger
DIM time as long

PRINT "ROM","FAST"

REM show it's as accurate
for i=1 to 15
    LET a=rnd * 32768
    PRINT SQR(a),SQRT(a)
next i

PRINT
PRINT "Over 500 Cycles:"
PRINT

REM ROM version 500 times.
LET time=t()
for i=1 to 500
    b=SQR(a)
next i
PRINT "Rom routine: ";t()-time;" Frames."

REM MY version 500 times.
LET time=t()
for i=1 to 500
    b=SQRT(a)
next i
PRINT "Fast routine: ";t()-time;" Frames."

PRINT
PRINT "PRESS A KEY"

PAUSE 1: PAUSE 0

CLS

PRINT "NUM   FAST      INTEGER"

REM show it's as accurate
for i=1 to 15
    LET a=INT(rnd * 32768)
    PRINT a;TAB 6;SQRT(a);TAB 16;SQRT16(a)
next i

PRINT
PRINT "Over 500 Cycles:"
PRINT

REM MY version 500 times.
LET time=t()
for i=1 to 500
    b=SQRT(a)
next i
PRINT "Fast routine: ";t()-time;" Frames."

REM MY Integer version 500 times.
LET time=t()
for i=1 to 500
    b=SQRT16(a)
next i
PRINT "Integer routine: ";t()-time;" Frames."

Code:

FUNCTION FASTCALL sqrtLF (num as uLong) as uInteger
REM incoming is DEHL
REM outbound is HL

asm
LD A,D
OR E
JP Z, sqrtLF16bit ; we're inside a 16 bit number. We can use the faster version.

LD b,16 ; b times round
EXX ; Out to root and rem - we're doing most of this in alternate registers.
 LD DE,0
 LD HL,0 ; DEHL = remainder
 LD BC,0 ; BC = root
EXX   ;back to num and loop
sqrtLFasmloop:
EXX  ; out to root and rem

SLA  C ; root <<= 1
RL  B   ;

SLA L ; rem=rem<<1
RL  H  ;
RL  E    ;
RL  D     ;
NOP
  
SLA L ; rem=rem<<1
RL  H  ;
RL  E    ;
RL  D     ;
EXX  ; back to Num and loop 

LD a,d    ; A = inputnum>>30
AND 192
RLCA
RLCA

SLA  L ; num <<= 1
RL  H
RL  E
RL  D

SLA  L ; num <<= 1
RL  H
RL  E
RL  D

EXX  ; out to root and rem

  ADD A,L     ; a=a+L              ; REM=REM+num>>30
  LD L,A      ; a-> L               ;
  JR NC, sqrtLFasmloophop1           ;
  INC H
  JR NC, sqrtLFasmloophop1
  INC DE           ;                               ;

sqrtLFasmloophop1:
  INC BC                       ; root=root+1
sqrtLFasmloophop2:

  ; DEHL = Remainder
  ; BC = root
  
  ; if rem >= root then
  LD A,D
  OR E
  JR NZ, sqrtLFasmthen ; if rem > 65535 then rem is definitely > root and we go to true
  
  LD A, H
  CP B
  JR C, sqrtLFasmelse ; H<B - that is rem<root so rem>=root is false and we go to else 
  JR NZ, sqrtLFasmthen ; H isn't zero though, so we could do a carry from it, so we're good to say HL is larger.
  
  ; if h is out, then it's down to L and C 
  LD A,L
  CP C
  JR C, sqrtLFasmelse ; L<C - that is rem<root so rem>=root is false and we go to else
  ; must be true - go to true.
      
  sqrtLFasmthen:
  ;remainder=remainder-root
  AND A ; clear carry flag
  SBC HL,BC ; take root away from the lower half of rem.
  JP NC, sqrtLFasmhop3 ; we didn't take away too much, so we're okay to loop round.

  ; if we're here, we did take away too much. We need to borrow from DE
  DEC DE ; borrow off DE
  
  sqrtLFasmhop3:
  INC BC ;root=root+1
  JP sqrtLFasmloopend
                  
  ;else 
  sqrtLFasmelse:
  DEC BC ;root=root-1
  ;end if
  
sqrtLFasmloopend:
  
EXX  ; back to num
  DJNZ sqrtLFasmloop
  
EXX ; out to root and rem
  PUSH BC
EXX ; back to normal
POP HL
   
    
SRA  H
RES 7,H
RR  L       ; Hl=HL/2 - root/2 is the answer.
jr sqrtLFexitFunction

sqrtLF16bit:

ld  a,l
ld  l,h
ld   de,0040h   ; 40h appends "01" to D
ld   h,d
ld b,7
sqrtLFsqrt16loop:
  sbc   hl,de      ; IF speed is critical, and you don't mind spending the extra bytes, you could unroll this loop 7 times instead of DJNZ.
  jr    nc,$+3   
  add   hl,de      
  ccf         
  rl    d      
  rla         
  adc   hl,hl      
  rla         
  adc   hl,hl      
DJNZ sqrtLFsqrt16loop
      
sbc   hl,de      ; optimised last iteration
ccf
rl   d
ld h,0
ld l,d 
ld de,0
sqrtLFexitFunction:   
end asm

END FUNCTION