vector.h

矩阵运算的模板类

	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 order
template <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 A
template <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