        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
