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