Euler 10

;----------------------------------------------------------
; Project Euler - Problem 10
;
; Find the sum of all the primes below two million.
;
; Requirements:
;
;   * at least a 386 with math coprocessor
;   * about 192 KiB free RAM
;
; Usage: euler10.com
;
; Prints the sum to standard output.
;
; Return codes (%errorlevel%):
;
;   * 0 Okay.
;   * 1 Not enough memory for the sieve data left.
;
    cpu 386
    org 100h
    [map all]

; TODO Move the stack down to save some space.
PROGRAM_SEGMENT_SIZE   equ 1000h   ; in paragraphs.

;----------------------------------------------------------
segment .data
; The limit up to which all prime numbers will be added.
; It is a variable instead of a constant which means the
; program may be extended to let the user enter this value
; at runtime.
limit:  dd 2000000

segment .bss
; Segment of the sieve data.  The data may span more than one
; 64 KiB segment.  In the case of the 2 million limit given in
; the Project Euler problem #10 it will take up about 128 KiB.
sieve_segment: resw 1

;----------------------------------------------------------
start:

segment .data
.banner_txt:
    db "Project Euler #10",13,10,13,10
    db "Sum all prime numbers below $"
.banner_txt2:
    db '.',13,10,13,10,'$'

segment .text
    mov ah,9
    mov dx,.banner_txt
    int 21h
    xor edx,edx
    mov eax,[limit]
    call print_u64
    mov ah,9
    mov dx,.banner_txt2
    int 21h

    call memory_init
    call sieve
    call add_primes
    call print_u64

    mov ax,4c00h
    int 21h

;----------------------------------------------------------
; Check if enough memory is available and calculate the
; base segment for the sieve data behind the program memory.
;
; In: limit The number up to which the prime numbers should
;           be added up.
; Out: sieve_segment The start segment of the sieve data.
memory_init:

segment .data
.mem_txt:
    db "Not enough memory.",13,10,'$'

segment .text
    mov eax,[limit]
    ; limit / 128 (/2 = only odd numbers, /8 = 8 bits per byte, /16 = paragraphs)
    shr eax,8
    inc ax      ; ("round" up.)
    jz .not_enough_memory_error
    add ax,PROGRAM_SEGMENT_SIZE
    jc .not_enough_memory_error
    cmp [2],ax
    jb .not_enough_memory_error

    ; Calculate and store the first segment address of the sieve data.
    mov ax,cs
    add ax,PROGRAM_SEGMENT_SIZE
    mov [sieve_segment],ax

    ret

.not_enough_memory_error:
    ; TODO Print amount of needed vs. available memory.
    mov ah,9
    mov dx,.mem_txt
    int 21h
    mov ax,4c01h
    int 21h

;----------------------------------------------------------
; Sieve of Eratosthenes implementation.
;
; In:   sieve_segment Start segment of the sieve data.
;       limit Up to this number the sieve determines prime numbers.
; Out:  [sieve_segment] Lookup bitarray with set bit for every odd prime number.
sieve:

segment .bss
.limit_sqrt:    resd 1

segment .text
    call fill_sieve

    ; Zeroeth bit of sieve data = 0
    mov bx,[sieve_segment]
    mov es,bx
    mov byte [es:0],11111110b

    ; Bounds for outer and inner loop.
    fild dword [limit]
    fsqrt
    fistp dword [.limit_sqrt]
    mov esi,[limit]
    shr esi,1

    mov edi,3           ; i (EDI)
.L1:
    mov eax,edi
    shr eax,1
    call get_bit
    jnc .not_prime

    mov edx,edi         ; j (EDX)
    shr edx,1
    add edx,edi
.L2:
    mov eax,edx
    call get_bit        ; Called to initialise AX, CX, and ES:BX.
    btr ax,cx
    mov [es:bx],al

    add edx,edi
    cmp edx,esi
    jb .L2

.not_prime:
    add edi,2
    cmp edi,[.limit_sqrt]
    jbe .L1

    ret

;----------------------------------------------------------
; Fills the sieve data memory with set bits.
;
; As it fills with 32 bit values it may fill up to three
; excess bytes with set bits.
;
; In:   sieve_segment
;       limit
; Out:  [sieve_segmet] All set bits.
fill_sieve:
    mov bx,[sieve_segment]
    mov es,bx

; Fill whole 64 KiB segments.
    mov eax,[limit]
    shr eax,20
    mov dx,ax
    mov eax,0ffffffffh
.L1:
    or dx,dx
    jz .whole_segments_finished
    xor di,di
    mov cx,4000h
    rep stosd

    mov bx,es
    add bx,1000h
    mov es,bx

    dec dx
    jmp .L1

.whole_segments_finished:
; Fill part of last 64 KiB segment.
    mov ecx,[limit]
    shr ecx,4
    inc ecx
    and ecx,0ffffh
    xor di,di
    rep stosd

    ret

;----------------------------------------------------------
; Adds all the prime numbers below [limit].
;
; In:   [limit]
; Out:  EDX:EAX     64 bit integer with the sum.
add_primes:

    xor edi,edi ; sum_hi (EDI)
    mov esi,2   ; sum_lo (ESI)
    mov edx,3   ; i (EDX)
.L1:
    mov eax,edx
    shr eax,1
    call get_bit
    jnc .not_prime

    add esi,edx
    jnc .skip
    inc edi
.skip:

.not_prime:
    add edx,2
    cmp edx,[limit]
    jb .L1

    mov eax,esi
    mov edx,edi
    ret

;----------------------------------------------------------
; Gets a bit from the bitarray (and a far pointer, and the byte value, and
; the bit offset within the byte value).
;
; in:   EAX bit offset.
;       sieve_segment
; out:  CF = bit value.
;       AL byte value.
;       ES:BX address of byte value.
;       CX bit offset within byte.
get_bit:
    mov cx,ax
    and cx,00000111b
    shr eax,3
    mov ebx,eax
    xor bx,bx
    shr ebx,4
    add bx,[sieve_segment]
    mov es,bx
    mov bx,ax
    mov al,[es:bx]
    bt ax,cx
    ret

;----------------------------------------------------------
; Prints an unsigned 64 bit integer value.
;
; In:   EDX:EAX Unsigned 64 bit integer value.
print_u64:

    mov edi,eax
    mov esi,edx
    mov ebx,10
    xor cx,cx           ; CX is digit counter.
.L1:
    mov eax,esi
    xor edx,edx
    div ebx
    mov esi,eax

    mov eax,edi
    div ebx
    mov edi,eax
    push dx
    inc cl

    mov eax,edi
    or eax,esi
    jnz .L1

.print:
    mov ah,2
    pop dx
    or dl,'0'
    int 21h
    loop .print

    ret

;----------------------------------------------------------