This code is written in HD64180 Assembler and runs under CP/M. I have a MicroMint SB180 single board computer which I run the Z-System on. Great stuff.


; MODULE NAME: SQRT.ASM
; FUNCTION: This program calculates the square root of 13 out to
;	 	some ungodly number of places.
;
; The method used here is:
;
;	given x and x^2 evaluate (x+1)(x+1) = x^2 + 2x + 1
;
; The x in this case is the part of the sqrt already known, the
; x^2 is known also.  The values x=36055 and x^2=1299963025 are
; compiled into the program as starting values.  The next guess
; at the value of the sqrt will be 36056 (ie x+1 = 36056).  Thus,
; to find the square of 36056, we calculate 2x (=36055+36055), add
; this to x^2 (=1299963025), and then add one.  If the '12' in x^2
; changes to a '13' we then know that 36056 is too big, and we should
; back up and try 360550.  If it does not change to '13' we know that
; this (x+1) is not large enough, so we try 36057.  In this case
; we can use the value of (x^2 + 2x + 1) as the new value of x^2.
; If the '12'->'13' we append a 0 to the end of x and append a 00
; to the end of x^2 and start the process again.
;
; Note that checking for a '12' or a '13' is sufficient to evaluate
; (x^2 + 2x +1).  Since the sqrt of 13 is endless the value of 13^2
; will be either less than 13 (12.999999...) or greater than 13
; (13.000...0000674829287...), but never _exactly_ 13.
;
; I use BCD here to pack two digits into a byte.  This allows the
; calculation of lots (relatively :)) of digits but is probably slower
; than say storing one digit per byte and being satisfied with fewer
; digits in the answer.
;
;===============================================================
;	.Z80
;
digits	EQU	13013		; number of digits wanted
DPL	EQU	65		; # of digits per line on console
BDOS	EQU	00005H
;
BOTTOM	EQU	00300H		; absolute bottom of available memory
TOP	EQU	0D800H		; absolute top of available memory
SIZE	EQU	TOP-BOTTOM	; max # of bytes of available memory
;
; Memory layout:
;
; Memory is allocated in quantities of 'allocation units'.  This
; is equal to one-fourth of the available memory.
;
aunit	equ	(size / 4)
;
; 2x is first in memory.  It is here because it can be written
; over by the 'savesqrt' program without hurting x or x^2.
; The least significant two digits (LS2D) are stored in the first byte.
; It grows up into memory.
;
ptr3	equ	bottom		; storage for 2x
;
; There are two copies of x^2.  One for the current value of x^2,
; and one for the 'being calculated' x^2.  Note that x^2 occupies
; twice as much space as x.  x^2 is stored most significant two digits
; (MS2D) in the lowest bytes of its allocated space and grows up into
; memory.  The MS2D never move, a new LSD is simply appended at the
; next available byte higher in memory.
;
ptr2a	equ	ptr3 + aunit	; 1st copy of x^2
aend	equ	ptr2a + 4	; LS2D of x^2 after init
ptr2b	equ	ptr2a + aunit	; 2nd copy of x^2
bend	equ	ptr2b + 4	; LS2D of x^2 after init
;
; x is stored LS2D in the top of its allocated space and grows down
; into memory.  The digits are moved all moved down and a new LSD
; is appended.
;
ptr1	equ	ptr2b + (2 * aunit) - 1
;
; storage for pointers to start and end of x and x^2, for use by the
; program 'savesqrt'.
;
spot1	equ	ptr2b + aunit	; spot for ptr to start of x
spot2	equ	spot1 + 2	; spot for ptr to end of x
spot3	equ	spot2 + 2	; spot for ptr to start of x^2
spot4	equ	spot3 + 2	; spot for ptr to end of x^2
;
	cseg
;
; Start by clearing memory and initializing x and x^2.
;
SQRT0:	LD	HL,BOTTOM	;Start of data memory
	LD	DE,BOTTOM+1	;Next byte
	LD	BC,SIZE-01H	;Amount available
	XOR	A		;Clear a
	LD	(HL),A		;Clear first
	LDIR			;Clear the rest
;
	LD	HL,DATA		;Move x into position
	LD	DE,PTR1-2	;(to here)
	LD	BC,3		;# of 'sqrt' data bytes
	LDIR
	LD	HL,DATA1	;Move x^2 into position
	LD	DE,PTR2A	;This is copy #1 of x^2
	LD	BC,5		;# of 'sqr' data bytes
	LDIR
	LD	HL,DATA1	;Move x^2 into position
	LD	DE,PTR2B	;This is copy #2 of x^2
	LD	BC,5		;# of 'sqr' data bytes
	LDIR
;
; Print what we know so far onto the console...
;
	CALL	LNNUM		;Display line number
	LD	A,'3'		;Display 3
	CALL	OUTT
	LD	A,'6'		;Display 6
	CALL	OUTT
	LD	A,'0'		;Display 0
	CALL	OUTT
	LD	A,'5'		;Display 5
	CALL	OUTT
	JP	SQRT3		;No need to check 1st time
;
; Check to see if (x^2 + 2x + 1) is > 13.
;
SEARCH:	LD	HL,(SRCHP1)	;Storage for current
	LD	A,(HL)		; 2nd byte of sqr
	CP	013H		;Was sqr > 13....?
	JP	Z,OVER		;Yes, sqrt too big, back up...
;
; It is < 13, so use the new value of (x^2 + 2x + 1) as the next
; x^2.  Do this by swapping some pointers...
;
	LD	DE,(XSQR1)	;Get 1st pointer to end of sqr
	LD	HL,(XSQR2)	;Get second
	LD	(XSQR2),DE	;Switch places
	LD	(XSQR1),HL
	LD	DE,(SRCHP1)	;Pointer to search byte 1
	LD	HL,(SRCHP2)	;Pointer to search byte 2
	LD	(SRCHP2),DE	;Switch places
	LD	(SRCHP1),HL
;
; Here is where 2x is calculated.
;
SQRT3:	LD	DE,PTR1		;tail-end of x
	LD	HL,PTR3		;ptr to storage for 2x
	LD	BC,(COUNT)	;compute this many places
	sra	b		; divide bc by 2 (BCD!)
	rr	c
	LD	A,B		;swap B and C
	LD	B,C		; B is inner loop counter for djnz
	ld	c,a		; C is outer loop counter now
	INC	B		; inc to negate 1st dec in loop
	INC	C		; inc to negate 1st dec in loop
	XOR	A		;Clear carry flag
;
SQRT1:	LD	A,(DE)		;Get 2 digits of x
	ADC	A,A		;Add so can use daa
	DAA			;Adjust the two digits
	LD	(HL),A		;Save 2x result
	INC	HL		;Point to next 2x spot
	DEC	DE		;Point to next two x digits
	DJNZ	SQRT1		; (decrements reg b)
	DEC	C
	JR	NZ,SQRT1
;
; This is where x^2, 2x, and 1 are added together.
;
	LD	DE,(XSQR1)	; ptr to LS2D of x^2
	LD	IY,(XSQR2)	;Put result in other x^2 string
	LD	HL,PTR3		; ptr to LS2D of 2x
	LD	BC,(COUNT)	;Compute this many places
	LD	A,B		; swap B and C
	LD	B,C		; B is inner loop counter for djnz
	ld	c,a		; C is outer loop counter now
	INC	B		; inc to negate 1st dec in loop
	INC	C		; inc to negate 1st dec in loop
	SCF			; the 1 in x^2 + 2x + 1
;
SQRT2:	LD	A,(DE)		;Get 2 digits of x^2
	ADC	A,(HL)		;Add to 2x
	DAA			;Adjust
	LD	(IY+0),A	;Save as new x^2
	INC	HL		;Next 2x pair
	DEC	DE		;Next x^2 pair
	DEC	IY		;Next store spot
	DJNZ	SQRT2		; (decrements B)
	DEC	C
	JR	NZ,SQRT2
;
; Here is where x is made into x+1.  If the LSD of x is a 9
; it will be incremented to hex digit 'a', but since we are
; actually calculating the value of x with a LSD of '10' this
; will be corrected by a decrement back to '9' in 'over:'.
;
	LD	HL,PTR1		;Point to 1st pair of x
	INC	(HL)		;Next estimate of sqrt
	JP	SEARCH		;How was last estimate?
;
; (x^2 + 2x + 1) is > 13, so we need to back up...
;
; First we check to see if the required number of digits of x
; have been calculated.  Checking here ensures that x^2 will
; be 13.(lots of zeros), rather than 12.(lots of nines).
; Personal preference...
;
; Since the new x^2 is > 13 the new value if x is too big, so
; the previous last digit of x is the one we want, so the COUNT
; value is incremented and we roll back to the previous x.
;
OVER:	XOR	A		;Clear carry
	LD	DE,(COUNT)	;Get number of places found
	LD	HL,digits	;How many wanted
	SBC	HL,DE		;This many?
	JP	C,SQEND		;Yes, done
;
	LD	HL,(COUNT)	;Point to digit count
	INC	HL		;Have found one place of x
	LD	(COUNT),HL	;Save
;
	LD	HL,(XSQR1)	;Pointer to LS2D of x^2 #1
	INC	HL		; append new LS2D
	LD	(XSQR1),HL	;Save
	LD	HL,(XSQR2)	;Pointer to LS2D of x^2 #2
	INC	HL		; append new LS2D
	LD	(XSQR2),HL	;Save
;
	LD	HL,PTR1		;Point to 1st pair of x
	DEC	(HL)		;Last est too large
;
; Uncommnet this line to skip the console printouts...
;
;	JR	OVERA
;
; Display the digit just found.  Display 'DPL' digits to a line,
; preceeded by the line number (in hex).
;
	LD	A,(HL)		; byte with the latest digit
	AND	0FH		; mask for it...
	OR	030H		; change to ASCII character
	CALL	OUTT		; display digit found
	LD	HL,VDCNT1	; # of digits on this line so far
	DEC	(HL)
	JR	NZ,OVERA	; jp if still more to go...
	LD	A,DPL		; this many digits per line
	LD	(HL),A
	LD	A,0DH		;Cr
	CALL	OUTT
	LD	A,0AH		;Lf
	CALL	OUTT
	CALL	LNNUM		;Display line number
;
; Here we shift the value of x so a '0' nibble can be inserted
; at the end.
;
OVERA:	LD	HL,PTR1		;Point to 1st pair of x
	LD	BC,(COUNT)	;Shift all bytes
	srl	b		; divide bc by 2 (BCD!)
	rr	c
	LD	A,B		; swap B and C
	LD	B,C		; B is inner loop counter for djnz
	ld	c,a		; C is outer loop counter now
	INC	B		; inc to negate 1st dec in loop
	INC	C		; inc to negate 1st dec in loop
	XOR	A		;Clear reg a
;
ROTATE:	RLD			;Move sqrt digits up one nibble
	DEC	HL		;Point ot next pair to shift
	DJNZ	ROTATE		; (decrements reg b)
	DEC	C
	JR	NZ,ROTATE
;
	JP	SQRT3		; ok, go do next digit...
;
; Increment and print the current line number on the console.
;
LNNUM:	LD	HL,VDCNT2	; storage spot for line #
	INC	(HL)		; inc it
	LD	C,(HL)
	CALL	OUTHX		; print it
	LD	A,':'
	CALL	OUTT
	LD	A,' '
	CALL	OUTT
	RET
;
; Print a BCD byte in reg C as 2 digits.
;
OUTHX:	LD	A,C
	RRA			;Get upper nibble
	RRA
	RRA
	RRA
	CALL	HEX1		;Convert upper nibble
	LD	A,C		;For lower nib
HEX1:	AND	0FH		;Mask out upper nib
	ADD	A,090H		;Convert nib to ascii
	DAA
	ADC	A,040H
	DAA
	CALL	OUTT		;Send out
	RET
;
; Print character in reg A on console.
;
OUTT:	PUSH	HL		;Save pointer
	PUSH	BC		;Save
	LD	E,A		;Byte to print
	LD	C,06H		;Output
	CALL	BDOS		;Send to CON:
	POP	BC		;Restore
	POP	HL		;Restore pointer
	RET	
;
; All done, print the last digit, terminate the line on the
; console, and exit.
;
SQEND:	LD	HL,PTR1		;Point to 1st pair of x
	LD	A,(HL)
	AND	0FH
	OR	030H		; it is a decimal digit...
	CALL	OUTT		;Display last digit
	LD	A,0AH
	CALL	OUTT
	LD	A,0DH
	CALL	OUTT
	LD	C,00H
;
	ld	hl,ptr1		; address of end of x
	inc	hl		; adjustment needed here...
	ld	(spot2),hl	; save for program 'savesqrt'
	ld	de,(count)	; # of digits
	srl	d		; divide de by two (BCD!)
	rr	e		; with carry from d
	xor	a		; clear carry flag
	sbc	hl,de		; calc start of x
	ld	(spot1),hl	; save for 'savesqrt'
;
	ld	hl,(srchp1)	; get ptr to start of x^2
	ld	(spot3),hl	; save for 'savesqrt'
	ld	hl,(xsqr2)	; get ptr to end of x^2
	inc	hl		; adjustment needed here...
	ld	(spot4),hl	; save for 'savesqrt'
;
	JP	BDOS		; exit
;
;Storage and pointers
;
COUNT:	DEFW	0006H		;Spot for # of places found so far
XSQR1:	DEFW	AEND		; spot for ptr to LS2D of x^2
XSQR2:	DEFW	BEND		; spot for ptr to LS2D of x^2
SRCHP1:	DEFW	PTR2B		; spot for ptr to MS2D of x^2
SRCHP2:	DEFW	PTR2A		; spot for ptr to MS2D of x^2
VDCNT1:	DEFB	9		;Char count ('01: 36055' is 9 char)
VDCNT2:	DEFB	0FFH		;Line count
;
DATA:	DEFB	003H		;At start x = 36055
	DEFB	060H
	DEFB	055H
;
DATA1:	DEFB	012H		;1st digit of x^2
	DEFB	099H		;36055 squared = 1299963025
	DEFB	096H
	DEFB	030H
	DEFB	025H
;
	END