Vector<T> C=A;
	for(int i=1;i<=C.Size();i++) C.vec[i]*=B.vec[i];
	return C;
}


template <typename T> Vector<T> mtl_divide(const Vector<T>& A, const Vector<T>& B){
	assert(A.Size()==B.Size());
	assert(A.type==B.type);
	Vector<T> C=A;
	for(int i=1;i<=C.Size();i++) C.vec[i]/=B.vec[i];
	return C;
}

template <typename T> ostream& mtl_ostream(ostream& s, const Vector<T>& A) {	s << "  ";	for(int i=1; i<=A.n; i++) {
		s << Num2str(A[i]) << " ";
	}	s << endl;	return s;}

template <typename T> istream& mtl_istream(istream& s, Vector<T>& A) {	for (int i=1; i<=A.n; i++)		s >> A[i];	return s;}

template <typename T> inline int Length(const Vector<T>& A) {
	return A.Size(); 
}


template <typename T> inline T Sum(const Vector<T>& A) {
	return A.Sum_(); 
}


template <typename T> inline T Max(const Vector<T>& A,int* idx=0){
	return A.Max_(idx); 
}


template <typename T> Matrix<T> Diag(const Vector<T> A){	Matrix<T> C(A.Size(),A.Size());	C=0;	for(int i=1;i<=NumRows(C);i++) C[i][i]=A[i];	return C;}

template <typename T> double Norm(const Vector<T>& A, int p=2){
	if(p==mtl_numeric_limits<int>::max()) {
		double maxA=0;
		for(int i=1;i<=A.Size();i++) if(maxA<Abs(A[i])) maxA=Abs(A[i]);
		return maxA;		
	}
	else if(p==mtl_numeric_limits<int>::min()) {
		double minA=mtl_numeric_limits<double>::max();
		for(int i=1;i<=A.Size();i++) if(minA>Abs(A[i])) minA=Abs(A[i]);
		return minA;
	}
	else return Vector<T>::normp(A,p);}


template <typename T> double Vector<T>::normp(const Vector<T>& A, int p) {
	double z=0, sum=0;
	int i;
	for(i=1;i<=A.Size();i++) if(z<Abs(A.vec[i])) z=Abs(A.vec[i]);
	if(z==0) return 0;
	for(i=1;i<=A.Size();i++) sum+=Ipow(Abs(A.vec[i])/z,p);
	return (z*pow(sum,1.0/p));
}

template <typename T> inline T Mean(const Vector<T>& A){	return A.Sum_()/A.n;}


template <typename T> Vector<T> Cumsum(const Vector<T>& A){
	Vector<T> C(A.n);
	T s=0;
	for(int i=1;i<=A.n;i++) { C(i)=s+A(i); s+=A(i); }
	return C;
}


template <typename T> Vector<T> Cumprod(const Vector<T>& A){
	Vector<T> C(A.n);
	T p=1;
	for(int i=1;i<=A.n;i++) { C(i)=p*A(i); p*=A(i); }
	return C;
}

template <typename T> string Num2str(const Vector<T>& A){	string s;	for(int i=1;i<=A.Size();i++){		if(i>1) s += " ";		s += Num2str(A[i]);	}	return s;}

// order=1: ascending order// order=-1: descending ordertemplate <typename T> int VecSort(const Vector<T>& A, Vector<T>& v, Vector<int>& orderA, int order){	int i;	int n_=A.Size();	v=A;	orderA.Resize(n_);	for(i=1;i<=n_;i++) orderA[i]=i;	for(i=1;i<=n_;i++){		for(int j=i+1;j<=n_;j++){			if(order*v[i] > order*v[j]){				T a=v[i];				v[i]=v[j];				v[j]=a;				int b=orderA[i];				orderA[i]=orderA[j];				orderA[j]=b;			}		}	}	return 1;}

template <typename T> inline Vector<T> Transpose(const Vector<T>& A){
	return A.ct();}

template <typename T> Matrix<T> OuterProd(const Vector<T>& A, const Vector<T>& B){	Matrix<T> C(A.Size(),B.Size());	for(int i=1;i<=A.Size();i++) for(int j=1;j<=B.Size();j++) C[i][j]=A.vec[i]*B.vec[j];	return C;}

template <typename T> T InnerProd(const Vector<T>& A, const Vector<T>& B){	assert(A.Size()==B.Size());	T sum=T();	for(int i=1;i<=A.Size();i++) sum+=A.vec[i]*B.vec[i];	return sum;}


template <typename T> Vector<T> Cross(const Vector<T>& A, const Vector<T>& B){
	assert(A.Size()==3 && B.Size()==3);
	Vector<T> C(3,A.Type());
	C[1]=A[2]*B[3]-A[3]*B[2];
	C[2]=A[3]*B[1]-A[1]*B[3];
	C[3]=A[1]*B[2]-A[2]*B[1];
	return C;
}


template <typename T> inline T Dot(const Vector<T>& A, const Vector<T>& B){
	return InnerProd(A,B);
}

template <typename T> inline Vector<T> ArrayMultiply(const Vector<T>& A, const Vector<T>& B){	return mtl_times(A,B);}

template <typename T> inline Vector<T> ArrayDivide(const Vector<T>& A, const Vector<T>& B){	return mtl_divide(A,B);}

// ToMatrix returns mxn matrix whose rows are repetitions of vector Atemplate <typename T> Matrix<T> ToMatrix(const Vector<T>& A, int m=1){
	assert(A.Type()==ROW_VECTOR);	Matrix<T> C(m,A.Size());	for(int i=1;i<=m;i++) C[i]=A;	return C;}


template <typename T> void Vector<T>::FromInt(const Vector<int>& A){
	Resize(A.Size());
	for(int i=1;i<=A.Size();i++) (*this)[i]=T(A[i]);
	return;
}


template <typename T> void Vector<T>::FromFloat(const Vector<float>& A){
	Resize(A.Size());
	for(int i=1;i<=A.Size();i++) (*this)[i]=T(A[i]);
	return;
}


template <typename T> void Vector<T>::FromDouble(const Vector<double>& A){
	Resize(A.Size());
	for(int i=1;i<=A.Size();i++) (*this)[i]=T(A[i]);
	return;
}


template <typename T> Vector<int> Vector<T>::ToInt(){
	Vector<int> C(Size());
	for(int i=1;i<=A.Size();i++) C[i]=(double)(*this)[i];
	return C;
}


template <typename T> Vector<float> Vector<T>::ToFloat(){
	Vector<float> C(Size());
	for(int i=1;i<=A.Size();i++) C[i]=(float)(*this)[i];
	return C;
}


template <typename T> Vector<double> Vector<T>::ToDouble(){
	Vector<double> C(Size());
	for(int i=1;i<=A.Size();i++) C[i]=(double)(*this)[i];
	return C;
}


template <typename T> Vector<int> Int(const Vector<T>& A){
	Vector<int> C(A.Size());
	for(int i=1;i<=A.Size();i++) C[i]=(int)A[i];
	return C;
}


template <typename T> Vector<float> Float(const Vector<T>& A){
	Vector<float> C(A.Size());
	for(int i=1;i<=A.Size();i++) C[i]=(float)A[i];
	return C;
}

template <typename T> Vector<double> Double(const Vector<T>& A){	Vector<double> C(A.Size());	for(int i=1;i<=A.Size();i++) C[i]=(double)A[i];	return C;}


template <typename T> Vector<T> Abs(const Vector<T>& A){
	Vector<T> C(A.Size());
	for(int i=1;i<=A.Size();i++) C.vec[i]=T(Abs(A.vec[i]));
	return C;
}


template <typename T> Vector<T> Sign(const Vector<T>& A){
	Vector<T> C(A.Size());
	for(int i=1;i<=A.Size();i++) C.vec[i]=T( ( (A.vec[i]>0) ? 1 : ((A.vec[i]<0) ? (-1) : 0) ) );
	return C;
}

template <typename T> Vector<T> Pow(const Vector<T>& A, double b){
	Vector<T> C(A.Size());
	for(int i=1;i<=A.Size();i++) C.vec[i]=pow((double)A.vec[i],b);
	return C;
}


template <typename T> Vector<T> Pow(const Vector<T>& A, int b){
	Vector<T> C(A.Size());
	for(int i=1;i<=A.Size();i++) C.vec[i]=Ipow(A.vec[i],b);
	return C;
}


template <typename T> Vector<T> Exp(const Vector<T>& A){	Vector<T> C(A.Size());	for(int i=1;i<=A.Size();i++) C.vec[i]=exp(A.vec[i]);	return C;}

template <typename T> Vector<T> Log(const Vector<T>& A){	Vector<T> C(A.Size());	for(int i=1;i<=A.Size();i++) C.vec[i]=Log(A.vec[i]);	return C;}

template <typename T> Vector<T> Digamma(const Vector<T>& A){	Vector<T> C(A.Size());	for(int i=1;i<=A.Size();i++) C.vec[i]=Digamma(A.vec[i]);	return C;}

template <typename T> Vector<T> Gammaln(const Vector<T>& A){	Vector<T> C(A.Size());	for(int i=1;i<=A.Size();i++) C.vec[i]=Gammaln(A.vec[i]);	return C;}


template <typename T> Vector<T> FFT(const Vector<T>& A,int N){
	Vector<T> C;
	Vector<T>::FFTBase(A,N,1,C);
	return C;
}


template <typename T> Vector<T> IFFT(const Vector<T>& A,int N){
	Vector<T> C;
	Vector<T>::FFTBase(A,N,-1,C);
	return C;
}


template <typename T> int Vector<T>::FFTBase(const Vector<T>& A,int N,int isign,Vector<T>& C){
	int log2n=0; int pow2n=1;
	for(log2n=0;pow2n<N;log2n++) pow2n*=2;
	C.Resize(2*pow2n);
	for(int i=1;i<=A.Size();i++) C[i]=A[i];
	Vector<T>::FFTCore(C,pow2n,isign);
	if(isign<0) C /= T(pow2n);
	return pow2n;
}


template <typename T> inline void FFT_SWAP(T& a,T& b){
	T temp=a;
	a=b;
	b=temp;
}


template <typename T> int Vector<T>::FFTCore(Vector<T>& data,int nn,int isign)
{
	int n,mmax,m,j,istep,i;
	double wtemp,wr,wpr,wpi,wi,theta;
	double tempr,tempi;

	n=nn << 1;
	j=1;
	for (i=1;i<n;i+=2) {
		if (j > i) {
			FFT_SWAP(data[j],data[i]);
			FFT_SWAP(data[j+1],data[i+1]);
		}
		m=n >> 1;
		while (m >= 2 && j > m) {
			j -= m;
			m >>= 1;
		}
		j += m;
	}
	mmax=2;
	while (n > mmax) {
		istep=2*mmax;
		theta=6.28318530717959/(isign*mmax);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		for (m=1;m<mmax;m+=2) {
			for (i=m;i<=n;i+=istep) {
				j=i+mmax;
				tempr=wr*data[j]-wi*data[j+1];
				tempi=wr*data[j+1]+wi*data[j];
				data[j]=data[i]-tempr;
				data[j+1]=data[i+1]-tempi;
				data[i] += tempr;
				data[i+1] += tempi;
			}
			wr=(wtemp=wr)*wpr-wi*wpi+wr;
			wi=wi*wpr+wtemp*wpi+wi;
		}
		mmax=istep;
	}
	return 1;
}


template <typename T> Vector<T> Merge(const Vector<T>& A, const Vector<T>& B)
{
	Vector<T> C(A.Size()+B.Size());
	int i=1,n=1;
	for(i=1;i<=A.Size();i++,n++) C[n]=A[i];
	for(i=1;i<=B.Size();i++,n++) C[n]=B[i];
	return C;
}


// ------------------------------------------------------------------------
// Vector utilities
// ------------------------------------------------------------------------
template <typename T> Vector<T> Vector<T>::Rand(int m)
{
	Vector<T> C(m);
	for(int i=1;i<=m;i++) C[i]=(float)RANDF();
	return C;
}


template <typename T> Vector<T> Vector<T>::Randn(int m)
{
	Vector<T> C(m);
	for(int i=1;i<=m;i++) C[i]=T(GaussianRandom());
	return C;
}


template <typename T> Vector<T> Vector<T>::CRandn(int m)
{
	Vector<T> C(m);
	for(int i=1;i<=m;i++) {
		C[i]=T(GaussianRandom(),GaussianRandom());
	}
	return C;
}


template <typename T> void Vector<T>::SetPermIdentity()
{
	for(int i=1;i<=n;i++) (*this)[i]=T(i);
}


#ifndef DISABLE_COMPLEX
#ifdef __GNUC__template <typename T> ostream& operator<<(ostream& s, const complex<T>& c) {   s << "(" << c.real() << "," << c.imag() << ")";   return s;}#endif
#endif#undef local_max#undef local_min
#endif