double quwl ( double rad )
{
	if ( rad < 0 ) return - 1 ;
 
	auto double xN = 1 ;
 
	if ( xN  < sqrt ( DBL_MIN )                 ||
		 xN  > sqrt ( DBL_MAX )                 ||
		 rad > DBL_MAX - xN * xN ) return - 1 ;
 
	auto double xNplus1 = 1 / ( 2 * xN ) * ( xN * xN + rad ) ;
 
	while ( fabs ( xNplus1 - xN ) > rad * 1e-14                       &&
						  xNplus1 > sqrt ( DBL_MIN ) - 1              &&
						  xNplus1 < sqrt ( DBL_MAX ) + 1              &&
			    			  rad < DBL_MAX - xNplus1 * xNplus1 + 1 )
	{
		xN = xNplus1 ;
 
		xNplus1 = 1 / ( 2 * xN ) * ( xN * xN + rad ) ;
	}
 
	return xNplus1 ;
}

while ( fabs ( xNplus1 - sqrt ( rad ) ) < fabs ( xN - sqrt ( rad ) )        &&
						    	    xNplus1 >= sqrt ( DBL_MIN )                 &&
						            xNplus1 <= sqrt ( DBL_MAX )                 &&
			    			            rad <= DBL_MAX - xNplus1 * xNplus1 )
	{
		xN = xNplus1 ;
 
		xNplus1 = 1 / ( 2 * xN ) * ( xN * xN + rad ) ;
	}
 
	return xN ;

double quwl ( double rad )
{
	if ( rad < 0 ) return - 1 ;
 
	auto double xN = rad > 2   ? log ( rad ) / log ( 2 ) :
					 rad > 0.5 ? rad / 2                 : 
					 rad > - 1 ? rad * 2                 : - 1 ;
 
	if ( xN  < sqrt ( DBL_MIN )                 ||
		 xN  > sqrt ( DBL_MAX )                 ||
		 rad > DBL_MAX - xN * xN ) return - 1 ;
 
	auto double xNplus1 = 1 / ( 2 * xN ) * ( xN * xN + rad ) ;
 
	while ( xNplus1 >= sqrt ( DBL_MIN )          &&
			xNplus1 <= sqrt ( DBL_MAX )          &&
			rad <= DBL_MAX - xNplus1 * xNplus1 )
	{
		if ( fabs ( xNplus1 - sqrt ( rad ) ) >  fabs ( xN - sqrt ( rad ) ) )
			 return xN ;
 
		if ( fabs ( xNplus1 - sqrt ( rad ) ) == fabs ( xN - sqrt ( rad ) ) )
			 return xNplus1 > xN ? xNplus1 : xN ;
 
		xN = xNplus1 ;
 
		xNplus1 = 1 / ( 2 * xN ) * ( xN * xN + rad ) ;
	}
 
	return xNplus1 ;
}