Code

int EllipsoidSolver(Vec3 E, Vec3 y,
                    float& error, float& resid, Vec3& x, Vec3& n,
                    float tol=1e-4, float rtol=1e-4,
                    int maxiter=10, bool useX=false)
{
    const float Rmin = 1e-4;
    const Vec3 One = Vec3(1.0f);

    const Vec3 y2 = y*y;
    const Vec3 E2 = E*E;
    const float yn = norm(y);
    const float R = min(E);

    float t;
    int iter = 0;

    // Check if y is too close to center of the ellipsoid.
    if (yn < Rmin)
        return -maxiter;

    // Start guess.
    if (useX)
        t = norm(y-x)/norm(x/E2);
    else
        t = -R + norm(y);

    while (iter < maxiter)
    {
        x = (E2/(Vec3(t)+E2))*y;
        Vec3 xp = -x/(Vec3(t)+E2);
        float F = dot(One, x*x/E2) - 1.0f;
        float Fp = dot(One, 2.0f*x*xp/E2);
        float dt = F/Fp;
        t = t - dt;
        iter = iter + 1;
        error = std::fabs(dt);
        resid = std::fabs(F);
        if (error < tol && resid < rtol)
        {
            n = x/E2;
            float rn = 1.0f/norm(n);
            n = rn*n;
            return iter;
        }
    }
    
    return -maxiter;
}