#include <iostream>
#include <cmath>
#include <cfloat>
using namespace std;

int const N = 1597463175;         // 0.00175128778162252   HARDER
// int const N = 1597463174;      // 0.001751301557861324  LOMONT

float inv_sqrt_newton( float x ) {
	float xhalf = 0.5f*x;
	int xi = *reinterpret_cast<int *>( &x );
	xi = N - (xi >> 1);
	x = *reinterpret_cast<float *>( &xi );
	return x*(1.5f - xhalf*x*x); 
}

int main() {
	cout.precision( 16 );

	double sumerror, sumsqrerror, maxerror;

	for ( float x = 1.0; x < 2.0; x += FLT_EPSILON ) {
		double y = 1.0/std::sqrt( static_cast<double>( x ) );
		double xx = static_cast<double>( inv_sqrt_newton( x ) );
		double rel = ( y - xx );
		sumerror += std::fabs( rel );
		sumsqrerror += ( rel )*( rel );
		maxerror = std::max( maxerror, std::fabs( rel ) );
	}

	for ( float x = 2.0; x < 4.0; x += 2.0*FLT_EPSILON ) {
		double y = 1.0/std::sqrt( static_cast<double>( x ) );
		double xx = static_cast<double>( inv_sqrt_newton( x ) );
		double rel = ( y - xx )/y;
		sumerror += std::fabs( rel );
		sumsqrerror += ( rel )*( rel );
		maxerror = std::max( maxerror, std::fabs( rel ) );
	}

	cout << "Sum of errors  " << "\t" << sumerror    << "\t" << sumerror << endl;
	cout << "Sum of squares " << "\t" << sumsqrerror << "\t" << sumsqrerror << endl;
	cout << "Maximum error  " << "\t" << maxerror    << "\t" << maxerror << endl;

	return 0;
}