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!