Difference between revisions of "Square roots in 6502 machine code"

From BeebWiki
Jump to: navigation, search
m (1 revision)
m (1 revision)
 
(No difference)

Latest revision as of 19:13, 8 March 2015

A simple way of calculating integer square roots is to count how many times increasing odd numbers can be subtracted from the starting value. For example:

32-1=31 -> 31-3=28 -> 28-5=23 -> 23-7=16 -> 16-9=7 -> 7-11<0

 1          2          3          4          5

-> SQR(30)=5, remainder 7

This can be done in 6502 machine code as follows:

   \ Calculate 16-bit square root
   \ ----------------------------
   \ On entry, (in,in+1)=input value
   \ On exit,  Y=(out)=root, X=(out+1)=remainder
   \ If out=in+2 this lets you use !in=value
   
   .sqr                            :\ On entry, !in=input value
   LDY #1:STY out+0:DEY:STY out+1  :\ Initialise out to first subtrand
   .sqr_loop                       :\ Repeatedly subtract increasing
   SEC                             :\   odd numbers until in<0
   LDA in+0:TAX:SBC out+0:STA in+0 :\ in=in-subtrand, remainder in X
   LDA in+1:SBC out+1:STA in+1
   BCC sqr_done                    :\ in<0, all done
   INY                             :\
   LDA out+0:ADC #1:STA out+0      :\ step +2 to next odd number
   BCC sqr_loop                    :\ no overflow, subtract again
   INC out+1:BNE sqr_loop          :\ INC high byte and subtract again
   .sqr_done                       :\ Y=root, X=remainder
   STY out+0:STX out+1             :\ out?0=root, out?1=remainder
   RTS

You can test this with:

   FOR A%=1 to 65535:!in=A%:CALL sqr:P.A%,out?0,out?1:NEXT

or with:

   FOR A%=1 to 65535:!in=A%:B%=USR sqr:P.A%,(B%AND&FF00)DIV256,(B%AND&FF0000)DIV65536:NEXT

It is simple to reduce the code to calculate roots of 8-bit numbers by removing the code to subtract in+1. The code uses an 8-bit counter for the root, so to calculate roots of numbers large than 16 bits different code is needed, as the square root of &10000 (a 17-bit number) is &100 (a 9-bit number).

Jgharston 21:37, 30 August 2007 (BST)