Midpoint ellipse algorithm in 8086 Assembly

I'm trying to implement an ellipse drawing algorithm with 8086 Assembly in VGA (320x200) mode, but the shape is not drawing properly for bigger radii.

As you can see the shape is concave or straight not elliptical. I get acceptable result when the axis have small lengths i.e (10, 15) but any more then that and it stops resembling an ellipse. I used an FPU to calculate all values.

    .387

SCREEN_WIDTH = 320
SCREEN_HEIGHT = 200
HALF_SCREEN_WIDTH = SCREEN_WIDTH / 2
HALF_SCREEN_HEIGHT = SCREEN_HEIGHT / 2

stack1 segment stack
            dw 300 dup(?)
stack_top   dw ?    ;stack_top pointer

stack1 ends

code1 segment
start1:
    mov ax,seg stack1 
	mov ss,ax								; initialize stack
	mov sp, offset stack_top

    mov al, 13h     
    mov ah, 0       ; change mode
    int 10h         ; enter graphics mode

    mov byte ptr cs:[color], 10
    
start_draw:
    call clear_screen
    call ellipse

wait_key:
    in  al, 60h ;take keyboard input
    cmp al, 1   ;if ESC end
    je end_program
    jmp wait_key

;-------------------------------------------------------------
;-------------------------------------------------------------
x_radius            dw 40   ;width of the ellipse
y_radius            dw 50   ;height of the ellipse

x_derivative        dw ?    ;derivates calculated from ellipse equation
y_derivative        dw ?

x_squared           dw ?    ;squared values for convenience
y_squared           dw ?

decision_param      dw ?    

four                dw 4    ;int values used in some of the expressions
two                 dw 2    ;storing them so I can load them to the FPU
;----------------------
calc_x_derivative:
;x_derivative = 2 * y_squared * x
;RPN: 2 y_squared x * *
    finit
    fild    word ptr cs:[y_squared]
    fimul   word ptr cs:[two] ;instead of multiplying by 2 
    fimul   word ptr cs:[x]
    fistp   word ptr cs:[x_derivative]
    ret
    
calc_y_derivative:
;y_derivative = 2 * x_squared * y
;RPN: 2 x_squared y * *
    finit
    fild    word ptr cs:[x_squared]
    fimul   word ptr cs:[two] 
    fimul   word ptr cs:[y]
    fistp   word ptr cs:[y_derivative]
    ret

ellipse:
    ;initialize x as 0 and y as y_radius
    mov     word ptr cs:[x], 0
    mov     ax, word ptr cs:[y_radius]
    mov     word ptr cs:[y], ax

x_radius_squared:
    finit 
    fild    word ptr cs:[y_radius]
    fimul   word ptr cs:[y_radius]
    fistp   word ptr cs:[y_squared]

y_radius_squared:
    finit
    fild    word ptr cs:[x_radius]
    fimul   word ptr cs:[x_radius]
    fistp   word ptr cs:[x_squared]

derivatives:
    call calc_x_derivative
    call calc_y_derivative
    
calc_decision_param_1:
;d1 = y_squared - (x_squared * ry) + (1/4 * x_squared)  
;RPN: y_squared x_squared ry * - x_squared 4 / +
    finit
    fild    word ptr cs:[y_squared]
    fild    word ptr cs:[x_squared]
    fild    word ptr cs:[y_radius]
    fmul
    fsub
    fild    word ptr cs:[x_squared]
    fidiv   word ptr cs:[four]
    fadd
    fistp   word ptr cs:[decision_param]

region_1_loop:
    call color_all_symmetries

    inc word ptr cs:[x] 
    call calc_x_derivative

    cmp word ptr cs:[decision_param], 0 ;compare decision_param to 0
    jge region_1_alternative_next       ;if its lower decrement y and calculate alternative decision_param

region_1_next:
;decision_param = decision_param + x_derivative + (y_squared)
;RPN: decision_param x_derivative + y_squared +
    finit
    fild word ptr cs:[decision_param]
    fild word ptr cs:[x_derivative]
    fadd
    fild word ptr cs:[y_squared]
    fadd
    fistp word ptr cs:[decision_param]
    jmp region_1_loop_end

region_1_alternative_next:
;decision_param = decision_param + x_derivative - y_derivative + (y_squared)
;RPN: decision_param x_derivative + y_derivative - y_squared +
    dec word ptr cs:[y]
    call calc_y_derivative
    finit
    fild word ptr cs:[decision_param]
    fild word ptr cs:[x_derivative]
    fadd
    fild word ptr cs:[y_derivative]
    fsub
    fild word ptr cs:[y_squared]
    fadd
    fistp word ptr cs:[decision_param]

region_1_loop_end:
    mov ax, word ptr cs:[x_derivative]
    cmp ax, word ptr cs:[y_derivative]     ;check if x_derivative < y_derivative
    jge calc_decision_param_2              ;if it is we are done with region 1
    jmp region_1_loop

calc_decision_param_2:
;decision_param = (y_squared * (x + 1/2) * (x + 1/2)) + (x_squared * (y - 1) * (y - 1)) - (x_squared * y_squared)
;RPN: y_squared x 1 2 / + * x 1 2 / + * x_squared y 1 - * y 1 - * + x_squared y_squared * -
    fild word ptr cs:[y_squared]
    fild word ptr cs:[x]
    fld1
    fidiv word ptr cs:[two]
    fadd
    fmul
    fild word ptr cs:[x]
    fld1
    fidiv word ptr cs:[two]
    fadd
    fmul
    fild word ptr cs:[x_squared]
    fild word ptr cs:[y]
    fld1
    fsub
    fmul
    fild word ptr cs:[y]
    fld1
    fsub
    fmul
    fadd
    fild word ptr cs:[x_squared]
    fimul word ptr cs:[y_squared]
    fsub
    fistp word ptr cs:[decision_param]

region_2_loop:
    call color_all_symmetries

    dec word ptr cs:[y]
    call calc_y_derivative

    cmp word ptr cs:[decision_param], 0 ;check if decision_param < 0
    jge alternative_region_2_next

region_2_next:
;decision_param = decision_param + x_squared - y_derivative
;RPN: decision_param x_squared + y_derivative -
    finit
    fild word ptr cs:[decision_param]
    fild word ptr cs:[x_squared]
    fadd
    fild word ptr cs:[y_derivative]
    fsub
    fistp word ptr cs:[decision_param]
    jmp region_2_loop_end


alternative_region_2_next:
    inc word ptr cs:[x]
    call calc_x_derivative

region_2_alternative_decision_param:
;decision_param = decision_param + x_derivative - y_derivative + (x_squared)
;RPN: decision_param x_derivative + y_derivative - x_squared +
    finit
    fild word ptr cs:[decision_param]
    fild word ptr cs:[x_derivative]
    fadd
    fild word ptr cs:[y_derivative]
    fsub
    fild word ptr cs:[x_squared]
    fadd
    fistp word ptr cs:[decision_param]

region_2_loop_end:
;y_derivative = y_derivative - (2 * x_squared)
;RPN: y_derivative 2 x_squared * -
    finit
    fild word ptr cs:[y_derivative]
    fild word ptr cs:[two]
    fild word ptr cs:[x_squared]
    fmul
    fsub
    fistp word ptr cs:[y_derivative]

    cmp word ptr cs:[y], 0      ; check if y <= 0
    jg region_2_loop       

end_ellipse:
    ret

;-------------------------------------------------------------
;---------------------  color_pixel  -------------------------

x       dw ?
y       dw ?

temp_y  dw ?
temp_x  dw ?

color   db ?
;---------------------
color_all_symmetries:
    ; Coloring points (xc + x, yc + y), (xc - x, yc + y), (xc + x, yc - y), (xc - x, yc - y)
    push ax

    ;Point (xc + x, yc - y)
    mov ax, HALF_SCREEN_HEIGHT
    sub ax, word ptr cs:[y]
    mov word ptr cs:[temp_y], ax
    mov ax, word ptr cs:[x]
    add ax, HALF_SCREEN_WIDTH
    mov word ptr cs:[temp_x], ax
    call color_pixel

    ;Point (xc - x, yc - y)
    mov ax, HALF_SCREEN_HEIGHT
    sub ax, word ptr cs:[y]
    mov word ptr cs:[temp_y], ax
    mov ax, HALF_SCREEN_WIDTH
    sub ax, word ptr cs:[x]
    mov word ptr cs:[temp_x], ax
    call color_pixel

    ;Point (xc + x, yc + y)
    mov ax, HALF_SCREEN_HEIGHT
    add ax, word ptr cs:[y]
    mov word ptr cs:[temp_y], ax
    mov ax, word ptr cs:[x]
    add ax, HALF_SCREEN_WIDTH
    mov word ptr cs:[temp_x], ax
    call color_pixel

    ;Point (xc - x, yc - y)
    mov ax, HALF_SCREEN_HEIGHT
    add ax, word ptr cs:[y]
    mov word ptr cs:[temp_y], ax
    mov ax, HALF_SCREEN_WIDTH
    sub ax, word ptr cs:[x]
    mov word ptr cs:[temp_x], ax
    call color_pixel

    pop ax
    ret

color_pixel:
    push ax
    push es
    ;position needs to be in a format y*320 + x
    mov ax, 0a000h
    mov es, ax
    mov ax, word ptr cs:[temp_y]    ;selecting a row
    mov bx, SCREEN_WIDTH            
    mul bx                          ;ax multiplying y by 320 to get proper position
    mov bx, word ptr cs:[temp_x]    ;selecting a column
    add bx, ax                      ;adding to the pixel position
    mov al, byte ptr cs:[color]     ;adding a color
    mov byte ptr es:[bx], al

    pop es
    pop ax
    ret
;-------------------------------------------------------------


end_program:
    mov al, 3h
    mov ah, 0
    int 10h

    mov ax, 4c00h
    int 21h

clear_screen:
    mov ax, 0a000h
    mov es, ax
    xor ax, ax
    mov di, ax
    mov cx, SCREEN_WIDTH * SCREEN_HEIGHT
    mov al, 0
    cld
    rep stosb   ;while cx != 0, mov byte ptr es:[di], al; di++; cx--
    ret

code1 ends

end start1

I tried many different version of the same algorithm all to no success. I think it might be an issue with my variables(especially the derivatives and decision parameter) overflowing I tried fixing it by using double words instead of words to declare variables but I had trouble comparing them later in the algorithm, and when I did get the program to compile it still did not work as intended.

> I think it might be an issue with my variables(especially the derivatives and decision parameter) overflowing

Your suspicion about variables overflowing is correct.
As an example, consider the computation for calc_decision_param_1:

> ;d1 = y_squared - (x_squared * ry) + (1/4 * x_squared)

The result for this is (50 * 50) - (40 * 40 * 50) + (40 * 40 / 4) = -77100
This doesn't fit a signed word anymore and therefore the fistp word ptr cs:[decision_param] instruction can't but store the default 'IntegerIndefinite value 8000h' in the word-sized destination.

What you should do is defining some of your variables as dword integers: x_derivative, y_derivative, and decision_param. All other variables will fit their current word size.

> but I had trouble comparing them later in the algorithm

The code that checks whether decision_param < 0 should simply look at the sign bit of the dword variable:

test byte ptr cs:[decision_param+3], 80h
jz   alternative_region_2_next

And the code that checks whether x_derivative < y_derivative could use ficomp to compare these signed dword integers. Based on the (conflicting) comments (;check if x_derivative < y_derivative and ;if it is we are done with region 1), it would look like:

  fninit
  fild    dword ptr cs:[x_derivative]
  ficomp  dword ptr cs:[y_derivative]
  fnstsw  ax
  sahf
  jb      calc_decision_param_2   ; done if x_derivative < y_derivative
  jmp     region_1_loop

calc_decision_param_2:

It would seem that x_radius_squared: and y_radius_squared: do exactly the opposite from what their names tell us!