#ifndef VECTOR_ARITHMETICS
#define VECTOR_ARITHMETICS

#include <vector>
#include <cmath>
using namespace std;


#pragma mark -
#pragma mark Vectors: Extending operations on STL vectors
#pragma mark


template<class TYPE>
inline vector<TYPE> operator*(vector<TYPE> x, const vector<TYPE> &y){
	for(unsigned long i=0; i<x.size(); ++i){
		x[i] *= y[i];
	}
	return x;
}


template<class TYPE>
inline vector<TYPE> operator+(vector<TYPE> x, const vector<TYPE> &y){
	for(unsigned long i=0; i<x.size(); ++i){
		x[i] += y[i];
	}
	return x;
}

template<class TYPE>
inline vector<TYPE> operator-(vector<TYPE> x, const vector<TYPE> &y){
	for(unsigned long i=0; i<x.size(); ++i){
		x[i] -= y[i];
	}
	return x;
}

template<class TYPE>
inline vector<TYPE> operator/(vector<TYPE> x, const vector<TYPE> &y){
	for(unsigned long i=0; i<x.size(); ++i){
		x[i] /= y[i];
	}
	return x;
}

template<class TYPE>
inline vector<TYPE> operator*(vector<TYPE> x, double scalar){
	for(unsigned long i=0; i<x.size(); ++i){
		x[i] *= scalar;
	}
	return x;
}

template<class TYPE>
inline vector<TYPE> operator/(vector<TYPE> x, double scalar){
	for(unsigned long i=0; i<x.size(); ++i){
		x[i] /= scalar;
	}
	return x;
}


inline vector<double> sqrt(vector<double> x){
	for(std::vector<double>::iterator i=x.begin(); i!=x.end(); ++i){
		*i = std::sqrt(*i);
	}
	return x;
}


inline vector<double> exp(vector<double> x){
	for(std::vector<double>::iterator i=x.begin(); i!=x.end(); ++i){
		*i = std::exp(*i);
	}
	return x;
}



inline vector<double> log(vector<double> x){
	for(std::vector<double>::iterator i=x.begin(); i!=x.end(); ++i){
		*i = std::log(*i);
	}
	return x;
}


inline vector<double> pow(vector<double> x, double power){
	for(std::vector<double>::iterator i=x.begin(); i!=x.end(); ++i){
		*i = std::pow(*i,power);
	}
	return x;
}


inline vector<double> abs(vector<double> x){
	for(std::vector<double>::iterator i=x.begin(); i!=x.end(); ++i){
		*i = std::abs(*i);
	}
	return x;
}



//pairwise maxima of vector entries (vectors must have same size)
inline vector<double> max(vector<double> x, const vector<double> &y){
	std::vector<double>::iterator i;
	std::vector<double>::const_iterator j;
	for(i=x.begin(), j=y.begin(); i!=x.end(); ++i, ++j){
		*i = max(*i, *j);
	}
	return x;
}


inline vector<double> max(vector<double> x, double b){
	std::vector<double>::iterator i;
	for(i=x.begin(); i!=x.end(); ++i){
		*i = max(*i, b);
	}
	return x;
}



// tensor product between a scalar and a vector of size DIM2
// resulting vector respresents a 1*DIM2 dimensional matrix stored in row-major order
// That is, matrix[r*DIM2+c] = p1[r]*p2[c]
inline vector<double> TensorProduct(double p1, const vector<double> &p2){
	unsigned int c, DIM2=p2.size();
	vector<double> p(1*DIM2);
	for(c=0; c<DIM2; ++c){
		p[c] = p1 * p2[c];
	}
	return p;
}


// tensor product between two vectors of sizes DIM1 and DIM2, respectively
// resulting vector respresents a DIM1*DIM2 dimensional matrix stored in row-major order
// That is, matrix[r*DIM2+c] = p1[r]*p2[c]
inline vector<double> TensorProduct(const vector<double> &p1, const vector<double> &p2){
	unsigned int r, c, DIM1 = p1.size(), DIM2=p2.size();
	vector<double> p(DIM1*DIM2);
	for(r=0; r<DIM1; ++r){
		for(c=0; c<DIM2; ++c){
			p[r*DIM2 + c] = p1[r] * p2[c];
		}
	}
	return p;
}




inline double TensorProduct(double p1, double p2){
	return p1*p2;
}


inline vector<double> TensorProduct(const vector<double> &p1, double p2){
	vector<double> p(p1.size());
	for(unsigned int i=0; i<p1.size(); ++i){
		p[i] = p1[i] * p2;
	}
	return p;
}



#endif