#include "abstractinterp.h"
#include <QString>

#include <math.h>

AbstractInterp::AbstractInterp( const QVector<double>& x,
                                const QVector<double>& y,
                                int m )
    : m_n( x.size() ),
      m_mm( m ),
      m_jsav( 0 ),
      m_correlated( false ),
      m_x( x ),
      m_y( y )
{
    m_dj = qMax( 1, int( pow( double( m_n ), 0.25 ) ) );
}

AbstractInterp::~AbstractInterp()
{
}

int AbstractInterp::locate( const double x )
{
    int ju, jm, jl;
    if ( m_n < 2 || m_mm < 2 || m_mm > m_n )
        throw QString( "locate size error" );

    bool ascnd = ( m_x.at( m_n - 1 ) >= m_x.at( 0 ) );
    jl = 0;
    ju = m_n - 1;
    while ( ju - jl > 1 )
    {
        jm = ( ju + jl ) >> 1;
        if ( ( x >= m_x.at( jm ) ) == ascnd )
            jl = jm;
        else
            ju = jm;
    }
    m_correlated = abs( jl - m_jsav ) > m_dj ? false : true;
    m_jsav = jl;
    return qMax( 0, qMin( m_n - m_mm, jl - ( ( m_mm - 2 ) >> 1 ) ) );
}

int AbstractInterp::hunt( const double x )
{
    int jl = m_jsav, jm, ju, inc = 1;
    if ( m_n < 2 || m_mm < 2 || m_mm > m_n )
        throw QString( "hunt size error" );

    bool ascnd = ( m_x.at( m_n - 1 ) >= m_x.at( 0 ) );
    if ( jl < 0 || jl > ( m_n - 1 ) )
    {
        jl = 0;
        ju = m_n - 1;
    }
    else
    {
        if ( ( x >= m_x.at( jl ) ) == ascnd )
        {
            for (;;)
            {
                ju = jl + inc;
                if ( ju >= m_n - 1 )
                {
                    ju = m_n - 1;
                    break;
                }
                else if ( ( x < m_x.at( ju ) ) == ascnd )
                    break;
                else
                {
                    jl = ju;
                    inc += inc;
                }
            }
        }
        else
        {
            ju = jl;
            for (;;)
            {
                jl = jl - inc;
                if ( jl <= 0 )
                {
                    jl = 0;
                    break;
                }
                else if ( ( x >= m_x.at( jl ) ) == ascnd )
                    break;
                else
                {
                    ju = jl;
                    inc += inc;
                }
            }
        }
    }

    while ( ju - jl > 1 )
    {
        jm = ( ju + jl ) >> 1;
        if ( ( x >= m_x.at( jm ) ) == ascnd )
            jl = jm;
        else
            ju = jm;
    }
    m_correlated = abs( jl - m_jsav ) > m_dj ? false : true;
    m_jsav = jl;
    return qMax( 0, qMin( m_n - m_mm, jl - ( ( m_mm - 2 ) >> 1 ) ) );
}

double AbstractInterp::interp( double x )
{
    int jlo = m_correlated ? hunt( x ) : locate( x );
    return rawInterp( jlo, x );
}
