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   [65000], #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
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   [65000], 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

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... ;-)
 
Top