# Calculate Square Root to two decimal places

#### westaust55

##### Moderator
Having noted comments/requests for an "improved Square-root function by forum members, and not having spied a routine posted to date,
below is a basic routine for calculating a square root of an integer number from 1 to 65535 to 2 decimal places.

I stopped at two decimal places for 2 reasons:
(a) the integer part can be up to 3 digits (255), which in a word variable leaves 2 for the mantissa/decimal part
(b) the divisors become larger than 65535 after 2 decimal places which would further complicate the maths

Nevertheless the 2 decimal place format and multiplied by 100 is in keeping with the PICAXE inbuilt Sin and COS functions.

Even to 2 decimal places, some internal math can exceed 65535 (numbers like ~500,000 occur) so overflow checking is required.
The code can no doubt be compressed but left as posted so that others may be able to understand how it functions.

The code did become shorter after I realised that the remainder determined from
Remaind = x - ( sqr(x) * sqr(x) )​
was always the same value as the remainder calculated by the math procedure and thus could be calculated in a couple of lines rather than many steps.

Code:
``````; Calculate square root of an integer number from 1 to 65535
;
; The result ia a value to 2 decimal places multiplied by 100
;
; Thus:
; sqr (2)     = 2.14   and is given as   214
; sqr (12358  = 111.16 and is given as 11116
; sqr (65535) = 255.99 and is given as 25599

; The routine is based on the methods described on various
; internet sources such as:
; www.solving-math-problems.com/calculate square root
; www.homeschoolmath.net/teaching/square-root-algorithm
;
;
; Picaxe program routine created by
; Westaust55 (Perth, Western Australia)
;
;
; Date: 30 Oct 2011
; Rev: 0
;
#PICAXE 40X1
#NO_DATA
#NO_TABLE

SYMBOL       x = w4
SYMBOL Remaind = w5
SYMBOL Intern1 = w6
SYMBOL OvrFlow = w7
SYMBOl Factor  = w8
SYMBOL Result  = w9

SYMBOL Counter = b24
SYMBOL Digit   = b25
SYMBOL Hunths  = b26
SYMBOL Charact = b27

; start of the program example
Main:
; Entry of a value
DataEntry:
SERTXD  ( "Input an value:")
SERRXD   , #x		; early 40X1 parts Must have the timeout - PE tells you so!
SERTXD (cr,lf,"Entered Value: ", #x, cr,lf)

Charact = SQR x					; Find the Characterist/Integer part using inbuilt function

Remaind = Charact * Charact			; determine the Remainder
Remaind = x - Remaind

Result = 0						; Clear the mantissa/demical part in case multiple calls

If Remaind = 0 Then Finalise			; If Remainder is 0 then nothing more to do

; Calculate the mantissa/decimal part to 2 decimal places
FOR Counter = 1 TO 2 				; one pass for each decimal place required - STOP AT 2

Factor = Charact 				; here we start to determine a new divisor
FOR Digit = 1 to Counter		; calculate 20 * Sqrt including current decimal places
Factor = Factor * 10
NEXT Digit
Factor = Result * 10 + Factor * 2 	; close, but we have still to add the units part

Intern1 = Remaind * 100			; we temporarily use var Intern1 in case we need to adjust

OvrFlow = Remaind ** 100		; check if we overflow the 65535 word variable limit
IF OvrFlow > 0 THEN
Factor = Factor / 10
Intern1 = Remaind + 5 /10 * 100 ; 5 is to round-up
ENDIF
Remaind = Intern1

Digit = Remaind / Factor		; a first pass but need to check may be too high

DO
Intern1 = Factor + Digit * Digit
IF Intern1 > Remaind THEN
DEC Digit			; too high so decrement the next digit of the square root
ENDIF
LOOP UNTIL Remaind >= Intern1
Remaind = Remaind - Intern1		; calculate the new remainder
Result = Result * 10 + Digit		; merge the next decimal digit into the result

NEXT Counter

Finalise:
Result = Charact * 100 + Result	; add in the Characteristic (ie Integer) part
SERTXD ("Sqrt = ", #Result, cr, lf)
GOTO Main``````

#### westaust55

##### Moderator
Having generated a decimal math based version felt there had to be a more compact method.
An example is given in Wikipedia: http://en.wikipedia.org/wiki/Methods_of_computing_square_roots
Will look soon to generating another version which could be more compact and faster . . .

As the old adage says: "more than one way to skin a cat"

#### westaust55

##### Moderator
After a few night working on a Square Root function using in effect binary math as opposed to decimal math, the results were not as good as hoped.

The core methodology is the same for either base 10 or base 2 math even if it does look a little different.
The next divisor part is computed from:
present square root value * 2 * base
hence the multiply by 20 in the decimal version and multiply by 4 in the binary version.

The code size is roughly the same for the core routine by both methods (182 bytes for decimal method versus 162 bytes for binary method).
Although in effect using 9 bits for the fractional (mantissa) part, the result is accurate to decimal places but a little variable thereafter.
Accuracy by observation compared with MS Excel indicates that results are typically within 0.0015 but occasionally may be out by ~0.005
so not really 3 decimal points accuracy as was hoped.

Again as per the decimal based version posted above, there was a need to monitor math overflow when multiplying large values.

Thanks also to hippy for spotting a bug when I could not see the forest for the trees late in the day with an overflow programming "hiccup".

Code:
``````; Calculate the Square Root for an integer value from 1 to 65535
;
;
; Picaxe routine created by Westaust55 (Perth, Western Australia)
;
; Accuracy by observation compared with MS Excel is that results
; are typically within 0.0015 but occasionally may be out by ~0.005
;
; Date: 4 Nov 2011
; Rev: 0
;
;
#PICAXE 40X1
#NO_DATA
#NO_TABLE

; define variable alias names
SYMBOL       x = w4  ; w4 =  b9:b8
SYMBOL Remaind = w5  ; w5 = b11:b10
SYMBOL SqrMant = w6  ; w6 = b13:b12
SYMBOL SqrRoot = w7  ; w7 = b15:b14
SYMBOL RemOFlw = b16
SYMBOL DivOFlw = b17
SYMBOL Charact = b18
SYMBOL Counter = b19

; start of the program example
Main:
; Entry of a value
SERTXD  ( "Input an value:")
SERRXD   , x		; early 40X1 parts Must have the timeout - PE tells you so!
SERTXD (cr,lf,"Entered Value: ", #x, cr,lf)
;
; Calculate the Integer Part (aka Characteristic) of the Square Root and the "Remainder"
Charact = SQR x				; Use the PICAXE routine to deterine the Characteristic/Integer part
SqrRoot = Charact				; we still need this value for the calculations
Remaind = Charact * Charact
Remaind = x - Remaind			; calc the remainder from which the Mantissa/decimal part is determined

RemOFlw = 0					; clear in case of repeat calls to this routine
DivOFlw = 0
SqrRoot = Charact * 2

For Counter = 0 to 8
RemOFlw = Remaind ** 4 + RemOFlw
Remaind = Remaind * 4

DivOFlw = SqrRoot ** 2
SqrRoot = SqrRoot AND \$FFFE * 2 OR 1

IF SqrRoot > Remaind AND RemOFlw <= DivOFlw THEN SkipOver
Remaind = Remaind - SqrRoot
SqrRoot = SqrRoot OR 2
IF RemOFlw = 0 THEN SkipOver
Remaind = Remaind + 65535
DEC RemOFlw
RemOFlw = REMOFlw * 4

SkipOver:
NEXT Counter
;
; For the mantissa we now have a value where bit15 = 0.50, bit14 = 0.250, bit13 = 0.125, etc
; Move the value 4 decimal places to the left to create an integer value that we can understand
; so for example, 0.3010 appears as %0101100101110111 (=dec 119) but after "shift" will become 3010
SqrMant = SqrRoot >> 1 << 7 ** 10000 				;

; now print out the result for Log10(x) for the entered value for x
SERTXD ("Sqrt(",#x,") = ",#Charact,".")
IF SqrMant <10 THEN
SERTXD ("000")
ELSEIF SqrMant <100 THEN
SERTXD ("00")
ELSEIF SqrMant <1000 THEN
SERTXD ("0")
ENDIF
SERTXD (#SqrMant, cr,lf,cr,lf)

GOTO Main``````
EDIT:
Part of the reason for the level of error is due to only working to 9 bits.
As can be seen from the attached table, when using only 9 bits, the 9th bit has a weighting of 0.00195 so typically the value could be almost 0.002 low since further bits are not included.
Adding a 10th bit would improve the resolution by adding 0.00098 and it would need 11 bits resolution (0.00049) to even try claim a "real" 3 decimal place accuracy.
Will be away again next week but may have another look later.
All that said, 2 decimal places is consistent with the accuracy available for the inbuilt trig funcitons.

#### Attachments

• 27.4 KB Views: 37
Last edited:

#### westaust55

##### Moderator
The code provided in posts 1 and 3 are more for the X1 and X2 PICAXE parts as they use the inbult SQR function to quickly determine the whole/integer part of the result.

For those using a M or M2 part, below is some code that will replicate the SQR function from the X1 and X2 parts providing the integer/whole number part. This could be extended using the same procedure for a decimal place if desired.

Code:
``````; Calculate square root of an integer number from 1 to 65535 giving an integer answer
;
;
; Picaxe program routine created by
; Westaust55 (Perth, Western Australia)
;
;
; Date: 31 Aug 2013
; Rev: 0
SYMBOL Entry	= w1
SYMBOL Intermed	= w2
SYMBOL Tempry2	= w3
SYMBOL Result 	= b8
SYMBOL Tempry1	= b9

Main:
; below are a few test cases - just remove the remark/comment mark (;) from one line and run in the PE simulator to watch the action
; Entry = 64	; sqrt = 8
; Entry = 625	; sqrt = 25
; Entry = 400 	; sqrt = 20
; Entry = 1000 	; sqrt=  31
; Entry = 2000	; sqrt = 44
; Entry = 23456	; sqrt = 153
; Entry = 42356	; sqrt = 205

Tempry1 = Entry / 10000
LOOKUP Tempry1, (0,1,1,1,2,2,2), Result
Intermed = Result * Result
Intermed = Tempry1 - Intermed * 100
Tempry1 = Tempry1 * 100
Tempry1 = Entry / 100 - Tempry1
Intermed = Intermed + Tempry1

Tempry1 = 255 ; = -1 so first increment = 0
DO
INC Tempry1
Tempry2 = Result * 20 + Tempry1 * Tempry1
LOOP UNTIL Tempry2 >= Intermed

IF Tempry2 > Intermed THEN : DEC Tempry1 : ENDIF
Tempry2 = Result * 20 + Tempry1 * Tempry1
Result = Result * 10 + Tempry1

Intermed = Intermed - Tempry2 * 100
Tempry2 = Entry // 100
Intermed = Intermed + Tempry2

Tempry1 = 255 ; = -1 so first increment = 0
DO
INC Tempry1
IF result <10 THEN
Tempry2 = Result * 20 + Tempry1 * Tempry1
ELSE
Tempry2 = Result/10*10
Tempry2 = Result + Tempry2 * 10 + Tempry1 * Tempry1
ENDIF
LOOP UNTIL Tempry2 >= Intermed

IF Tempry2 > Intermed THEN : DEC Tempry1 : ENDIF
; next line only needed if you want to go further to decimal places
; Tempry2 = Result * 20 + Tempry1 * Tempry1
Result = Result * 10 + Tempry1

SERTXD (#RESULT,cr,lf)
PAUSE 1000``````
The below version has an outer looping structure which saves around 20 bytes but uses an extra byte variable.
Code:
``````; Calculate square root of an integer number from 1 to 65535 giving an integer answer
;
;
; Picaxe program routine created by
; Westaust55 (Perth, Western Australia)
;
;
; Date: 1 Sept 2013
; Rev: 0

SYMBOL Entry	= w1
SYMBOL Intermed	= w2
SYMBOL Tempry2	= w3
SYMBOL Result 	= b8
SYMBOL Tempry1	= b9
SYMBOL Counter	= b10

Main:

; Entry = 16	; sqrt = 4
; Entry = 20	; sqrt = 4
; Entry = 25	; sqrt = 5
; Entry = 64	; sqrt = 8
Entry = 625	; sqrt = 25
; Entry = 400 	; sqrt = 20
; Entry = 1000 	; sqrt=  31
; Entry = 2000	; sqrt = 44
; Entry = 23456	; sqrt = 153
; Entry = 42356	; sqrt = 205

Tempry1 = Entry / 10000
LOOKUP Tempry1, (0,1,1,1,2,2,2), Result ; find hunderds digit
Intermed = Result * Result
Intermed = Tempry1 - Intermed * 100
Tempry1 = Tempry1 * 100
Tempry1 = Entry / 100 - Tempry1
Intermed = Intermed  + Tempry1

For Counter = 0 to 1 ; for tens and untis digits
IF counter = 1 THEN
Intermed = Intermed - Tempry2 * 100
Intermed = Entry // 100 + Intermed
ENDIF
Tempry1 = 255 ; = -1 so first increment = 0
DO
INC Tempry1
IF result <10 THEN
Tempry2 = Result * 20 + Tempry1 * Tempry1
ELSE
Tempry2 = Result/10*10
Tempry2 = Result + Tempry2 * 10 + Tempry1 * Tempry1
ENDIF
LOOP UNTIL Tempry2 >= Intermed
IF Tempry2 > Intermed THEN : DEC Tempry1 : ENDIF
Tempry2 = Result * 20 + Tempry1 * Tempry1
Result = Result * 10 + Tempry1
Next Counter
SERTXD (#RESULT,cr,lf)
PAUSE 1000``````

Last edited:

#### womai

##### Senior Member
For integer numbers, the iterative Babylonian method just uses a few lines of code. Note: Due to integer arithmetic it works only for numbers up to 32767, not 65535.

Code:
``````w0 = 32767 ' input
w1 = 1 ' will contain result (square root of w0)

for b4 = 1 to 12
w1 = w0 / w1 + w1 / 2
next b4``````
The choice of 12 iterations guarantees convergence to the closest integer for all values up to 32767. One could limit the number of iterations by checking whether the result has changed compared to the previous iteration, but that requires storing the previous result - needs one more word variable.

#### AllyCat

##### Senior Member
Hi,

A forum search for a "high resolution" Square Root algorithm for a current project found this thread, which seems worth resurrecting. As usual, Westy's code seems to run impeccably, although I haven't examined how it all works. But it now appears that for the present requirements Bill Wann may have devised an elegantly simple solution in post #8 here.

For integer numbers, the iterative Babylonian method just uses a few lines of code..........The choice of 12 iterations guarantees convergence to the closest integer for all values up to 32767. One could limit the number of iterations by checking whether the result has changed compared to the previous iteration, but that requires storing the previous result - needs one more word variable.

The Babelonian method does look appealingly simple, but PICaxe's integer maths produces a few "issues" which I thought are worth reporting here.

Firstly, it doesn't generally produce the "closest" integer to the square root, but the "truncated" (lower) value (as does the PICaxe (X2 family only) SQRT function and the "Digit by Digit" method). It does appear to work for numbers above 32767, but the choice of the "initial estimate" value may be important. For a 16-bits calculation, a starting value of about 50 seems optimum (around 6 iterations for the smallest and largest numbers) but 255 may be needed to avoid overflow errors when the number approaches 65535).

However, an issue that I haven't seen reported before is that for input values of the form N * N - 1 (e.g. 8 and 15) the integer solution oscillates between N and (N - 1), so that if the algorithm checks (only) for equality of the "previous" and "present" roots, then it never finishes and the program hangs. Code:
``````#picaxe 08m2                  ; And most others
for w3 = 0 to 65535         ; Test range
w0 = w3                    ; Test value - Edit as required
;    w0 = w3 * w3 - 1         ; "Problem" values
w1 = 50                 ; Will contain result (50 gives max iterations = 6)
for b8 = 1 to 12
w2 = w1           ; "Previous" value
w1 = w0 / w1 + w1 / 2
if w2 = w1 then exit
next b8
sertxd("Root(",#w0,")= ",#w1," iterations= ",#b8,cr)
next``````
One solution is to exit when the Previous and Present values differ by one unit, which may also speed up the algorithmslightly (by avoiding an "unnecessary" iteration that just leads to the same value again). However, in this case (or whenever an "arbitrary" number of iterations are used), the program does not know whether the lower (N-1) or higher (N) solution has been found. This probably doesn't matter if the result is not to be rounded to the nearest value, but it is a little disconcerting that the result can depend on arbitrary initial conditions.

A rule to determine when the result should be "rounded up" is if the "remainder" (i.e. number - root2) exceeds the found Root. Therefore, I have refined the algorithm to not only terminate with the "minimum" number of iterations, but optionally to either round up to the nearest root, or at least to determine the roots consistently.

Code:
``````#picaxe 08m2            ; And most others
#no_data
#define roundup        ; Root = nearest square root to number
symbol test   = w0
symbol number = w1
symbol root   = w2
symbol last   = w3
for test = 0 to 65535
number = test                       ; All numbers
;    number = test * test - 1             ; "Problem" values
Gosub FindSquareRoot
SerTxd(cr,lf, "SQRT(", #number, ") = ",#root)
next
end

FindSquareRoot:             ; ~28 bytes without rounding (46 with rounding)
root = 50                       ; First estimate (50 gives max ~6 iterations)
do
last = root
root = number / root + root / 2        ; New estimate
last = last - root                       ; Compare previous and present root values
loop until last < 2                      ; Loop until equal or 1 digit higher
#ifdef roundup
last = root + 1 * root - number       ; Root - Remainder (or +1 if already rounded up)
if last > 32767 then                   ; Remainder > Root (=negative value) so round up
inc root
#else
last = root * root
if last > number then                 ; Already rounded up
dec root
#endif
endif
Return``````
Cheers, Alan.

#### premelec

##### Senior Member
Not sure if it's of use - I used to use a routine on a Friden mechanical calculator [60 years ago] which involved subtracting 1-3-5-7-9 counting subtractions and shifting... Also taught my kids about squares with Y by Y marble arrays having them see how many marbles needed to be added to get next square of marbles... 1-3-5-7 etc +2 each time more than last time... My uncle - a retired engineer who died at age 104 - said he used to put himself to sleep doing square roots in his head - I don't know which process he used... ;-)