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