#include <vector>
#include <arithmetic_iterators/arithmetic_iterator.h>
#include <arithmetic_iterators/ai_square.h>
#include <arithmetic_iterators/linalg.h>



namespace traditional {
    template <class P1, class P2, class O1>
    inline void subtract(const P1 &x, const P2 &y, O1 &diff) {
	int l = x.size();
	for(int i=0; i<l; ++i)
	    diff[i] = x[i]-y[i];
    }
    template <class P1, class P2>
    inline vector<typename P1::value_type> subtract(const P1 &x, const P2 &y) {
	vector<typename P1::value_type> diff(x.size());
	subtract(x,y,diff);
	return diff;
    }

    template <class P1, class O1>
    inline void square(const P1 &x, O1 &sq) {
	int l = x.size();
	for(int i=0; i<l; ++i)
	    sq[i] = x[i]*x[i];
    }
    template <class P1>
    inline vector<typename P1::value_type> square(const P1 &x) {
	vector<typename P1::value_type> sq(x.size());
	square(x, sq);
	return sq;
    }
}




template <class P1, class P2>
typename P1::value_type test_lazy(const P1 &x, const P2 &y)
{
    return sum(square(x - y));
}


template <class P1, class P2>
typename P1::value_type test_hand(const P1 &x, const P2 &y)
{
    typename P1::const_iterator xit = x.begin();
    typename P2::const_iterator yit = y.begin();
    typename P1::value_type acc = 0;

    for(; xit != x.end(); ++xit, ++yit)
	acc += (*xit-*yit)*(*xit-*yit);
    return acc;
}


template <class P1, class P2, class VEC1, class VEC2>
typename P1::value_type test_hand2(const P1 &x, const P2 &y, VEC1 &diff, VEC2 &sq)
{
    int l = x.size();
    for(int i=0; i<l; ++i)
	diff[i] = x[i] - y[i];
    for(int i=0; i<l; ++i)
	sq[i] = diff[i]*diff[i];
    typename P1::value_type acc = 0;
    for(int i=0; i<l; ++i)
	acc += sq[i];
    return acc;
}



template <class P1, class P2>
typename P1::value_type test_traditional(const P1 &x, const P2 &y)
{
    return sum(traditional::square(traditional::subtract(x,y)));
}


template <class P1, class P2, class VEC1, class VEC2>
typename P1::value_type test_traditional2(const P1 &x, const P2 &y, VEC1 &diff, VEC2 &sq)
{
    traditional::subtract(x,y, diff);
    traditional::square(diff, sq);
    return sum(sq);
}


using namespace NTPM;


int
main()
{
    const int L = 10000;
    const int ITERS = 1000;
    float x[L];
    float y[L];


    for(int i=0; i<L; ++i)
    {
	x[i] = (float)i*0.3;
	y[i] = (float)i*0.2;
    }

    vector<float> x_v(x,x+L);
    vector<float> y_v(y,y+L);
    vector<float> diff(L);
    vector<float> sq(L);

    /* Build objects that promise the data. */
    NTPMai::Promise<float*> x_l = make_promise(x,x+L);
    NTPMai::Promise<float*> y_l = make_promise(y,y+L);

    cout << "lazy = " << test_lazy(x_l, y_l) << endl;
    cout << "hand = " << test_hand(x_v, y_v) << endl;
    cout << "hand = " << test_hand2(x_v, y_v, diff, sq) << endl;
    cout << "traditional = " << test_traditional(x_v, y_v) << endl;
    cout << "traditional = " << test_traditional2(x_v, y_v, diff, sq) << endl;


    /* Evaluate sum of squared difference of x and y. */

    for(int i=0; i<ITERS; ++i)
	test_hand2(x_v, y_v, diff, sq);

    for(int i=0; i<ITERS; ++i)
	test_lazy(x_l, y_l);

    for(int i=0; i<ITERS; ++i)
	test_hand(x_v, y_v);

    for(int i=0; i<ITERS; ++i)
	test_traditional(x_v, y_v);

    for(int i=0; i<ITERS; ++i)
	test_traditional2(x_v, y_v, diff, sq);
}