double
SQRT (double a, double approx)
{
doubleinterpolate;

if (approx <= 0.0)
approx = 1.0;
do {
interpolate = ((a / approx) - approx) / 2;
approx += interpolate;
} while (fabs (interpolate) >= 1.0E-6);

return approx;
}