//
// === FAST INVERSE SQUARE ROOT  ===
// From "Graphics Gems V", Alan Paeth (Editor)
// ISBN 0125434553/9649-1547332-306386
// Published by Ap Profession, 1995
//

#include <iostream>
#include <iomanip>
#include <cmath>

#include <sys/time.h>
#include <time.h>

#include "frsqrt.hh"

#define ITER 300

static double
second ()
{
    struct timeval  tvd;

    gettimeofday( & tvd, NULL );

    return tvd.tv_sec + (tvd.tv_usec * 1e-6);
}

int
main ( int argc, char ** argv )
{
    init_rsqrt();

    const double  res3 = 1.0 / std::sqrt(3.0);
    const double  res5 = 1.0 / std::sqrt(5.0);
    
    std::cout << std::setprecision( 16 ) << std::scientific
              << "first example: 1 / sqrt(3)" << std::endl
              << "  exact  = "  << res3 << std::endl
              << "  float  = " << rsqrtf( 3.0f )
              << std::setprecision( 4 ) << std::scientific
              << ", error = " << std::abs( res3 - rsqrtf( 3.0f ) ) / res3 << std::endl
              << std::setprecision( 16 ) << std::scientific
              << "  double = " << rsqrt( 3.0 )
              << std::setprecision( 4 ) << std::scientific
              << ", error = " << std::abs( res3 - rsqrt( 3.0 ) ) / res3 << std::endl;
    
    std::cout << std::setprecision( 16 ) << std::scientific
              << "second example: 1 / sqrt(5)" << std::endl
              << "  exact  = " << res5 << std::endl
              << "  float  = " << rsqrtf( 5.0f )
              << std::setprecision( 4 ) << std::scientific
              << ", error = " << std::abs( res5 - rsqrtf( 5.0f ) ) / res5 << std::endl
              << std::setprecision( 16 ) << std::scientific
              << "  double = " << rsqrt( 5.0 )
              << std::setprecision( 4 ) << std::scientific
              << ", error = " << std::abs( res5 - rsqrt( 5.0 ) ) / res5  << std::endl;

    std::cout << std::endl << "Benchmark" << std::endl << std::endl;
    
    {
        float    d1, d2, d;
        double   start, end;
        
        start = second();
        
        d1 = 0.0;
        for ( unsigned i = 0; i < ITER; i++ )
            for ( d = 2.0; d < 1000000.0; d += 4.0 )
            {
                d1 += 1.0f / std::sqrt( d   );
                d1 += 1.0f / std::sqrt( d+1.0f );
                d1 += 1.0f / std::sqrt( d+2.0f );
                d1 += 1.0f / std::sqrt( d+3.0f );
            }// for
        
        end = second();
        std::cout << "(float)  time for 1.0 / sqrt = "
                  << std::setprecision( 2 ) << std::fixed << end - start << " sec"
                  << std::setprecision( 16 ) << std::scientific << " (res = " << d1 << ")" << std::endl;
        
        start = second();
        
        d2 = 0.0;
        for ( unsigned i = 0; i < ITER; i++ )
            for ( d = 2.0; d < 1000000.0; d += 4.0 )
            {
                d2 += rsqrtf( d );
                d2 += rsqrtf( d+1.0f );
                d2 += rsqrtf( d+2.0f );
                d2 += rsqrtf( d+3.0f );
            }// for
        
        end = second();
        std::cout << "(float)  time for      rsqrt = "
                  << std::setprecision( 2 ) << std::fixed << end - start << " sec"
                  << std::setprecision( 16 ) << std::scientific << " (res = " << d2 << ")" << std::endl;
    }

    {
        double   d1, d2, d;
        double   start, end;
        
        start = second();
        
        d1 = 0.0;
        for ( unsigned i = 0; i < ITER; i++ )
            for ( d = 2.0; d < 1000000.0; d += 4.0 )
            {
                d1 += 1.0 / std::sqrt( d   );
                d1 += 1.0 / std::sqrt( d+1.0 );
                d1 += 1.0 / std::sqrt( d+2.0 );
                d1 += 1.0 / std::sqrt( d+3.0 );
            }// for
        
        end = second();
        std::cout << "(double)  time for 1.0 / sqrt = "
                  << std::setprecision( 2 ) << std::fixed << end - start << " sec"
                  << std::setprecision( 16 ) << std::scientific << " (res = " << d1 << ")" << std::endl;
        
        start = second();
        
        d2 = 0.0;
        for ( unsigned i = 0; i < ITER; i++ )
            for ( d = 2.0; d < 1000000.0; d += 4.0 )
            {
                d2 += rsqrt( d );
                d2 += rsqrt( d+1.0 );
                d2 += rsqrt( d+2.0 );
                d2 += rsqrt( d+3.0 );
            }// for
        
        end = second();
        std::cout << "(double)  time for      rsqrt = "
                  << std::setprecision( 2 ) << std::fixed << end - start << " sec"
                  << std::setprecision( 16 ) << std::scientific << " (res = " << d2 << ")" << std::endl;
    }
    
    return 0;
}