Superquadric surfaces play an important role in many different research fields due to their ability to model a variety of shapes with a small number of parameters. One of their core applications is the estimation of object shape characteristics from a set of discrete samples from the surface of an object. However, the corresponding optimization problem is prone to numerical instabilities in some regions of the parameter space. To mitigate this problem, lower bound constraints to the shape parameters are applied during the optimization thus limiting the range of shapes which can be accurately represented by superquadrics. Therefore, the exact modeling of very common shapes like cuboids and cylinders is error-prone in practice. In this article we investigate this problem and provide a numerically stable formulation for the evaluation of the superquadric surface function and for its gradient. This new formulation enables numerically stable fitting of superquadrics in the previously constrained region, i.e., in its full range including cuboids and cylinders. In addition, the new formulation also leads to faster convergence speed. The theoretical contributions are substantiated by experiments on synthetic as well as real data.