Problem 42580. Conic equation

A conic of revolution (around the z axis) can be defined by the equation

   s^2 – 2*R*z + (k+1)*z^2 = 0

where s^2=x^2+y^2, R is the vertex radius of curvature, and k is the conic constant: k<-1 for a hyperbola, k=-1 for a parabola, -1<k<0 for a tall ellipse, k=0 for a sphere, and k>0 for a short ellipse.

Write a function z=conic(s,R,k) to calculate height z as a function of radius s for given R and k. Choose the branch of the solution that gives z=s^2/(2*R)+... for small values of s. This defines a concave surface for R>0 and a convex surface for R<0.

The trick is to get full machine precision for all values of s and R. The test suite will require a relative error less than 4*eps, where eps is the machine precision.

Hint (added 2015/09/03): the straightforward solution is

   z = (R-sqrt(R^2-(k+1)*s^2))/(k+1), 

but this does not work if k=-1, gives the wrong branch of the solution if R<0, and is subject to severe roundoff error if s^2 is small compared to R^2. It is possible, however, to find a mathematically equivalent form of the solution that solves all three problems at once.