I am making a 3D renderer that displays a cube based on some given points and edges. It uses this formula for projecting the 3D points onto a 2D screen:

x_proj = (x*f)/(z+f)

y_proj = (y*f)/(z+f)

where f is the focal length and x_proj and y_proj are the coordinates for the points to be placed on the screen. My problem is that the cube I’m displaying looks more like a cuboid.
This is what it looks like:
Needs to be a cube, looks like a cuboid
The points are stored in a list and is projected based of a projection function.

points = [[-50, -50, 50],
          [-50, 50, 50],
          [50, 50, 50],
          [50, -50, 50],
          [-50, -50, -50],
          [-50, 50, -50],
          [50, 50, -50],
          [50, -50, -50]]

edges = [[0, 4],
         [0, 3],
         [0, 1],
         [1, 5],
         [1, 2],
         [2, 3],
         [2, 6],
         [4, 5],
         [4, 7],
         [5, 6],
         [6, 7],
         [3, 7]]

# Other functions

def project(pos, lines, fcl):
    pos_prj = []
    for i in pos:
        x = i[0]
        y = i[1]
        z = i[2]
        try:
            x_prj = (x*fcl)/(z+fcl)
            y_prj = (y*fcl)/(z+fcl)
        except:
            x_prj = x
            y_prj = y
        draw_point([x_prj, y_prj])
        pos_prj.append([x_prj, y_prj])
    for i in lines:
        draw_edge(pos_prj[i[0]], pos_prj[i[1]])

# Projection function is called later in the main loop

I tried changing the formula I used a bit (changing the operators, modifying the formula, etc.) but it didn’t work.