	title	'Eratosthenes Sieve for 80x86 Real Mode'
	name	sieve
	page	50,132

;
; Eratosthenes Sieve for 80x86 Real Mode
; Implemented by Ray Duncan, April 1987
; After Gilbreath, Byte September 1981, and January 1983
;

niter   equ     100             ; number of iterations
asize   equ     16380           ; size of array "flags"

cr	equ	0dh			; ASCII carriage return
lf	equ	0ah			; ASCII line feed

stdin	equ 	0			; handle for standard input
stdout	equ	1			; handle for standard output


_TEXT	segment para public 'CODE'

	assume cs:_TEXT,ds:_DATA,es:_DATA

sieve	proc    near

	mov	ax,seg _DATA
	mov	ds,ax
	mov	es,ax

	mov	dx,0			; convert number of iterations
	mov	ax,niter
	mov	cx,10
	mov	si,offset msg1a+3
	call	binasc

	mov	dx,offset msg1		; display message
	mov	cx,msg1_len		; "Starting N iterations of Sieve"
	call	pmsg

	call	getmsec			; get current time in msec
	push	dx			; and save it ...
	push	ax

        mov     counter,niter		; initialize iterations counter

sieve1:					; a sieve iteration starts here...
	mov	di,offset flags		; initialize flags array
	mov	cx,asize		; to all bytes = TRUE
	mov	al,1
	cld
	rep stosb

 	xor	si,si			; SI=  index to flags array
	xor	di,di			; DI = primes counter

sieve2:     				; main loop of sieve
        test	byte ptr flags[si],1	; is this a prime?
	jnz	short sieve4		; jump if prime

sieve3:	inc     si			; bump to next slot in "flags"
        cmp     si,asize		; are we done?
        jle	sieve2			; jump to test another

 	dec     word ptr counter	; more iterations?
        jnz     sieve1			; jump, another iteration needed
	jmp	sieve7

sieve4:					; prime found, zap its multiples
	mov	bx,si			; copy i
 	mov     dx,bx 			; DX = prime = i + i + 3
        add	dx,dx
        add     dx,3
	xor	al,al
	jmp	short sieve6

sieve5:	mov	byte ptr flags[bx],al	; zero this multiple

sieve6:	add     bx,dx			; find next multiple of prime
        cmp     bx,asize		; have we exhausted the array?
        jle	sieve5			; not yet, zap it
	inc     di			; count primes and try next
	jmp	sieve3

sieve7:					; done with all iterations...
	call	getmsec			; get current time
	push	dx			; and save it...
	push	ax

	mov	ax,di			; convert number of primes
	mov	dx,0
	mov	cx,10
	mov	si,offset msg2a+4
	call	binasc

	mov	dx,offset msg2		; display "Number of primes: "
	mov	cx,msg2_len
	call	pmsg

	pop	ax			; stop time:  low word
	pop	dx			;             high word

	pop	bx			; start time: low word
	pop	cx			;             high word

	sub	ax,bx			; DX:AX = stop - start 
	sbb	dx,cx

	mov	cx,niter		; divide by number of iterations
	idiv	cx

	mov	dx,0			; convert msec to ASCII
	mov	cx,10
	mov	si,offset msg3a+4
	call	binasc

	mov	dx,offset msg3		; display "Elapsed time:"
	mov	cx,msg3_len
	call	pmsg

 	mov     ax,04C00h		; exit to DOS with
        int     21H			; return code = 0

sieve	endp


getmsec	proc	near			; DX:AX := current time in msec.

	mov	ah,2ch			; read time from MS-DOS	
	int	21h
	push	dx			; save seconds, hundredths
	mov	al,ch			; AX := hours
	cbw
	mov	bx,60			; hours -> minutes
	imul	bx
	xor	ch,ch			; isolate system minutes
	add	ax,cx			; and find total minutes
	mov	bx,60			; minutes -> seconds
	imul	bx
	pop	cx			; recover seconds, hundredths
	mov	bl,ch			; get seconds
	xor	bh,bh
	add	ax,bx			; find total seconds
	adc	dx,0			; carry if necessary
	xor	ch,ch			; save centisec.
	mov	bp,cx
					; total seconds * 100 the hard way
	mov	bx,ax			; double multiply * 10 
	mov	cx,dx
	add	ax,ax			; * 2
	adc	dx,dx
	add	ax,ax			; * 4
	adc	dx,dx
	add	ax,bx			; * 5
	adc	dx,cx
	add	ax,ax			; * 10
	adc	dx,dx

	mov	bx,ax			; double multiply * 10 
	mov	cx,dx
	add	ax,ax			; * 2
	adc	dx,dx
	add	ax,ax			; * 4
	adc	dx,dx
	add	ax,bx			; * 5
	adc	dx,cx
	add	ax,ax			; * 10
	adc	dx,dx

	add	ax,bp			; add in hundredths of seconds
	adc	dx,0
					; now convert total to msec.
	mov	bx,ax			; double multiply * 10 
	mov	cx,dx
	add	ax,ax			; * 2
	adc	dx,dx
	add	ax,ax			; * 4
	adc	dx,dx
	add	ax,bx			; * 5
	adc	dx,cx
	add	ax,ax			; * 10
	adc	dx,dx

	ret				; return DX:AX = msec.

getmsec	endp


;
; BINASC: Convert 32 bit binary value to ASCII string.
;
; Call with  DX:AX = signed 32 bit value
;	     CX    = radix
;            SI    = last byte of area to store resulting string
;	             (make sure enough room is available to store
;		      the string in the radix you have selected.)
;
; Destroys AX, BX, CX, DX, and SI.
;
binasc 	proc	near			; convert DX:AX to ASCII.

	mov	byte ptr [si],'0' 	; force storage of at least 1 digit.
	or	dx,dx			; test sign of 32 bit value,
	pushf				; and save sign on stack.
	jns	bin1			; jump if it was positive.
	not	dx			; negative, take 2's complement
	not	ax			; of the value. 
	add	ax,1
	adc	dx,0
bin1:					; divide 32 bit value by radix 
					; to extract next digit for the
					; forming string.
	mov	bx,ax			; is the value zero yet?
	or	bx,dx
	jz	bin3			; yes, we are done converting.
	call	divide			; no, divide by radix.
	add	bl,'0'			; convert remainder to ASCII digit.
	cmp	bl,'9'			; might be converting to hex ASCII,
	jle	bin2			; jump if in range 0-9,
	add	bl,'A'-'9'-1		; correct it if in range A-F.
bin2:	mov	[si],bl			; store this character into string.
	dec	si			; back up through string,
	jmp	bin1			; and do it again.
bin3:					; restore sign flag,
	popf				; was original value negative?
	jns	bin4			; no, jump
	mov	byte ptr [si],'-'	; yes,store sign into output.
bin4:	ret				; back to caller.

binasc	endp

;
; General purpose 32 bit by 16 bit unsigned divide.
; This must be used instead of the plain machine unsigned divide
; for cases where the quotient may overflow 16 bits (for example,
; dividing 100,000 by 2).  If called with a zero divisor, this
; routine returns the dividend unchanged and gives no warning.
;
; Call with DX:AX = 32 bit dividend
;           CX    = divisor
;
; Returns   DX:AX = quotient
;           BX    = remainder
;	    CX    = divisor (unchanged)
;
divide	proc	near			; Divide DX:AX by CX

	jcxz	div1			; exit if divide by zero
	push	ax			; 0:dividend_upper/divisor
	mov	ax,dx
	xor	dx,dx
	div	cx
	mov	bx,ax			; BX = quotient1
	pop	ax			; remainder1:dividend_lower/divisor
	div	cx
	xchg	bx,dx			; DX:AX = quotient1:quotient2

div1:	ret				; BX = remainder2

divide	endp



pmsg	proc	near			; print a message on std output
					; call with DS:EDX = address 
					;           ECX    = length 
	mov	ah,40h
	mov	bx,stdout
	int	21h
	ret

pmsg	endp


_TEXT	ends


_DATA	segment para public 'DATA'

flags	db      asize+1 dup (?)
counter	dw      ?

msg1	db	cr,lf,'Starting '
msg1a	db	'     iterations of Sieve...',cr,lf
msg1_len equ $-msg1

msg2	db	cr,lf,'Primes found:  '
msg2a	db	'     ',cr,lf
msg2_len equ $-msg2

msg3	db	cr,lf,'Elapsed time:  '
msg3a	db	'       msec. per iteration',cr,lf
msg3_len equ $-msg3

_DATA	ends


_STACK	segment byte stack 'stack'
	db      4096 dup (?)
_STACK	ends

	end	sieve

