C++语言的复数与矩阵运算库
总共包含三个文件:COMPLEX.H
matrix.h
text.cpp
本程序打包下载见本贴4楼
COMPLEX.H如下：
/* -------------------------------------------------------------------- */
/* Z++ Version 1.10 complex.h Last revised 10/10/92 */
/* */
/* Complex number class for Turbo C++/Borland C++. */
/* Copyright 1992 by Carl W. Moreland */
/* -------------------------------------------------------------------- */
#ifndef COMPLEXdotH
#define COMPLEXdotH
#include <math.h>
#include <iostream.h>
const unsigned char Z_RADIANS = 0;
const unsigned char Z_DEGREES = 1;
const unsigned char Z_COMMA = 0; // (x, y)
const unsigned char Z_LETTER= 1; // x + iy
class complex
{
public:
double re, im;
private:
static unsigned char zArgMode;
static unsigned char zPrintMode;
static unsigned char zLetter;
public:
complex(void): re(0), im(0) {}
complex(const double real, const double imag=0): re(real), im(imag) {}
complex(const complex& z): re(z.re), im(z.im) {}
friend double re(const complex& z) { // real part
return z.re;
}
friend double im(const complex& z) { // imaginary part
return z.im;
}
friend doublereal(const complex& z) { // real part
return z.re;
}
friend doubleimag(const complex& z) { // imaginary part
return z.im;
}
friend double mag(const complex& z) { // magnitude |z|
return sqrt(z.re*z.re + z.im*z.im);
}
friend double arg(const complex& z); // argument
friend double ang(const complex& z) { // angle
return arg(z);
}
friend double ph(const complex& z) { // phase
return arg(z);
}
friend complex conj(const complex& z) { // complex conjugate
return complex(z.re, -z.im);
}
friend doublenorm(const complex& z) { // norm
return z.re*z.re + z.im*z.im;
}
friend complex rtop(double x, double y=0);
friend complex ptor(double mag, double angle=0);
complex& topolar(void);
complex& torect(void);
void operator = (const complex& z) { // z1 = z2
re = z.re;
im = z.im;
}
complex& operator += (const complex& z) { // z1 += z2
re += z.re;
im += z.im;
return *this;
}
complex& operator -= (const complex& z) { // z1 -= z2
re -= z.re;
im -= z.im;
return *this;
}
complex& operator *= (const complex& z) { // z1 *= z2
*this = *this * z;
return *this;
}
complex& operator /= (const complex& z) { // z1 /= z2
*this = *this / z;
return *this;
}
complex operator + (void) const { // +z1
return *this;
}
complex operator - (void) const { // -z1
return complex(-re, -im);
}
friend complex operator + (const complex& z1, const complex& z2) {
return complex(z1.re + z2.re, z1.im + z2.im);
}
friend complex operator + (const complex& z, const double x) {
return complex(z.re+x, z.im);
}
friend complex operator + (const double x, const complex& z) {
return complex(z.re+x, z.im);
}
friend complex operator - (const complex& z1, const complex& z2) {
return complex(z1.re - z2.re, z1.im - z2.im);
}
friend complex operator - (const complex&, const double);
friend complex operator - (const double x, const complex& z) {
return complex(x-z.re, -z.im);
}
friend complex operator * (const complex& z1, const complex& z2) {
double re = z1.re*z2.re - z1.im*z2.im;
double im = z1.re*z2.im + z1.im*z2.re;
return complex(re, im);
}
friend complex operator * (const complex& z, const double x) {
return complex(z.re*x, z.im*x);
}
friend complex operator * (const double x, const complex& z) {
return complex(z.re*x, z.im*x);
}
friend complex operator / (const complex&, const complex&);
friend complex operator / (const complex& z, const double x) {
return complex(z.re/x, z.im/x);
}
friend complex operator / (const double, const complex&);
friend complex operator ^ (const complex& z1, const complex& z2) {
return pow(z1, z2);
}
friend int operator == (const complex& z1, const complex& z2) {
return (z1.re == z2.re) && (z1.im == z2.im);
}
friend int operator != (const complex& z1, const complex& z2) {
return (z1.re != z2.re) || (z1.im != z2.im);
}
friend double abs(const complex& z);
friend complex sqrt(const complex& z);
friend complex pow(const complex& base, const complex& exp);
friend complex pow(const complex& base, const double exp);
friend complex pow(const double base, const complex& exp);
friend complex exp(const complex& z);
friend complex log(const complex& z);
friend complex ln(const complex& z);
friend complex log10(const complex& z);
friend complexcos(const complex& z);
friend complexsin(const complex& z);
friend complextan(const complex& z);
friend complex asin(const complex& z);
friend complex acos(const complex& z);
friend complex atan(const complex& z);
friend complex sinh(const complex& z);
friend complex cosh(const complex& z);
friend complex tanh(const complex& z);
void SetArgMode(unsigned char mode) const {
if(mode == Z_RADIANS || mode == Z_DEGREES)
zArgMode = mode;
}
void SetPrintMode(unsigned char mode) const {
if(mode == Z_COMMA || mode == Z_LETTER)
zPrintMode = mode;
}
void SetLetter(unsigned char letter) const {
zLetter = letter;
}
friend ostream& operator<<(ostream&, const complex&);
friend istream& operator>>(istream&, const complex&);
};
static const complex Z0(0, 0); // complex number 0
static const complex Z1(1, 0); // complex number 1
static const complex Zi(0, 1); // complex number i
static const complex Zinf(HUGE_VAL, HUGE_VAL); // complex number infinity
static const complex Complex;
/* -------------------------------------------------------------------- */
/* Here is the same class with the name capitalized. If you prefer this */
/* simply remove the comment delimiters below and comment out the 5 */
/* static globals above. */
/* -------------------------------------------------------------------- */
/*
class Complex: public complex
{
public:
Complex(void): re(0), im(0) {}
Complex(const double real, const double imag=0): re(real), im(imag) {}
Complex(const complex& z): re(z.re), im(z.im) {}
};
static const Complex Z0(0, 0); // complex number 0
static const Complex Z1(1, 0); // complex number 1
static const Complex Zi(0, 1); // complex number i
static const Complex Zinf(HUGE_VAL, HUGE_VAL); // complex number infinity
static const Complex complex;
*/
#endif
[ 本帖最后由 风花雪月 于 2006-10-10 00:44 编辑 ] matrix.h如下：
#include "afxwin.h"
#include "strstrea.h"
#include "fstream.h"
#include "iomanip.h"
#include "stdlib.h"
#include "malloc.h"
#include "math.h"
#define CDMatrix (CMatrix<double>)
#define CIMatrix (CMatrix<int>)
template<class Type>
class CMatrix :public CObject
{
protected:
Type **f;
int row;
int col;
public:
intGetRow(){return row;}
intGetCol(){return col;}
void SetRow(int r){row=r;}
void SetCol(int l){col=l;}
Type** GetF(){return f;}
CMatrix();
CMatrix(int,int);
CMatrix(CMatrix ©);
~CMatrix();
BOOL Ones(int);
BOOL Ones(int,int);
BOOL Zeros(int);
BOOL Zeros(int,int);
CMatrix Turn();
CMatrix Inverse();
CMatrix Construct(CMatrix &mat,CString style);
CMatrix GetChild(int,int,int,int);
BOOL Diag(Type *,int);
BOOL Malloc(int,int);
CMatrix operator*(CMatrix&m);
friend CMatrixoperator*(double,CMatrix);
friend CMatrixoperator*(CMatrix,double);
CMatrix operator/(Type x);
CMatrix operator+(CMatrix&m);
CMatrix operator+(Type x);
CMatrix operator-(CMatrix&m);
CMatrix operator-(Type x);
CMatrix operator-();
Type* &operator[](int);
friend istream&operator>>(istream&,CMatrix&);
friend ostream&operator<<(ostream&,CMatrix&);
void operator=(const CMatrix&m);
void operator=(Type);
Type ** MatrixAlloc(int,int);
void MatrixFree(Type **);
void Hessenberg();
BOOL Cholesky();
CMatrix QREigen(double,int);// return eigenvalue in two row matrix of a real matrix
CMatrix Jacobi(double eps,int jt);
};
template<class Type>
Type* &CMatrix<Type>::operator[](int i)
{
Type *tp;
tp=(Type *)f;
return tp;
}
template<class Type>
CMatrix<Type>::CMatrix()
{
f=NULL;
row=0;
col=0;
}
template<class Type>
CMatrix<Type>::CMatrix(int m,int n)
{
f=MatrixAlloc(m,n);
row=m;
col=n;
}
template<class Type>
CMatrix<Type>::~CMatrix()
{
if(f!=NULL) MatrixFree(f);
}
template<class Type>
BOOL CMatrix<Type>::Ones(int m)
{
f=MatrixAlloc(m,m);
for(int i=0;i<m;i++)
for(int j=0;j<m;j++)
{
if(i!=j)f=0;
else f=1;
}
row=m;
col=m;
return TRUE;
}
template<class Type>
BOOL CMatrix<Type>::Ones(int m,int n)
{
f=MatrixAlloc(m,n);
for(int i=0;i<m;i++)
for(int j=0;j<n;j++)
{
f=1;
}
row=m;
col=n;
return TRUE;
}
template<class Type>
BOOL CMatrix<Type>::Zeros(int m)
{
f=MatrixAlloc(m,m);
for(int i=0;i<m;i++)
for(int j=0;j<m;j++)
{
f=0;
}
row=m;
col=m;
return TRUE;
}
template<class Type>
BOOL CMatrix<Type>::Zeros(int m,int n)
{
f=MatrixAlloc(m,n);
if(f==NULL) return FALSE;
for(int i=0;i<m;i++)
for(int j=0;j<n;j++)
{
f=0;
}
row=m;
col=n;
return TRUE;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::operator*(CMatrix<Type>&mat)
{
int l,m,n;
l=GetRow();
m=GetCol();
n=mat.GetCol();
CMatrix tmp(l,n);
for(int k=0;k<l;k++)
for(int j=0;j<n;j++)
{
for(int i=0;i<m;i++)
tmp.f+=f*mat.f;
}
return tmp;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::operator/(Type x)
{
int l,n;
l=GetRow();
n=GetCol();
CMatrix<Type> tmp(l,n);
for(int j=0;j<l;j++)
for(int i=0;i<n;i++)
tmp.f=f/x;
return tmp;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::operator+( CMatrix<Type>&mat)
{
int l1,m1;
l1=GetRow();
m1=GetCol();
CMatrix<Type> tmp(l1,m1);
for(int j=0;j<l1;j++)
for(int i=0;i<m1;i++)
tmp.f=f+mat.f;
return tmp;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::operator+(Type x)
{
int l,m;
l=GetRow();
m=GetCol();
CMatrix<Type> tmp(l,m);
for(int j=0;j<l;j++)
for(int i=0;i<m;i++)
tmp.f=f+x;
return tmp;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::operator-(CMatrix<Type>&mat)
{
int l1,m1;
l1=GetRow();
m1=GetCol();
CMatrix<Type> tmp(l1,m1);
for(int j=0;j<l1;j++)
for(int i=0;i<m1;i++)
tmp.f=f-mat.f;
return tmp;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::operator-()
{
for(int i=0;i<row;i++)
for(int j=0;j<col;j++)
f=-f;
return *this;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::operator-( Type x)
{
int l,m;
l=GetRow();
m=GetCol();
CMatrix<Type> temp(l,m);
for(int j=0;j<l;j++)
for(int i=0;i<m;i++)
temp.f=f-x;
return temp;
}
template<class Type>
istream &operator>>(istream &in,CMatrix<Type> &mat)
{
int i,j;
for(i=0;i<mat.row;i++)
for(j=0;j<mat.col;j++)
in>>mat.f;
cout<<endl;
return in;
};
template<class Type>
ostream& operator<<(ostream& out,CMatrix<Type> &v1)
{
if(v1.GetF()==NULL)
{
out<<"This Matrix cannot be output!";
return out<<endl;
}
out<<setiosflags(ios::right||ios::fixed);
out<<setprecision(5);
out<<"["<<endl;
int mr=v1.GetRow();
int mc=v1.GetCol();
for(int j=0;j<mr;j++)
{
for(int i=0;i<mc-1;i++)
out<<setw(12)<<v1.f;
out<<setw(12)<<v1.f<<";"<<endl;
}
return out<<"]"<<endl;
}
template<class Type>
CMatrix<Type>::CMatrix(CMatrix<Type> ©)
{
f=NULL;
*this=copy;
}
template<class Type>
void CMatrix<Type>::operator=(const CMatrix<Type> ©)
{
if(this==©) return;
if(f!=NULL) MatrixFree(f);
int m,n;
m=copy.row;
n=copy.col;
f=MatrixAlloc(m,n);
for(int i=0;i<m;i++)
for(int j=0;j<n;j++)
f=copy.f;
row=m;
col=n;
}
template<class Type>
void CMatrix<Type>::operator=(Type d)
{
if(f==NULL)
{
cout<<"The Matrix is not be allociated"<<endl;
return;
}
for(int i=0;i<row;i++)
for(int j=0;j<col;j++)
f=d;
}
template<class Type>
Type ** CMatrix<Type>::MatrixAlloc(int r,int c)
{
Type *x,**y;
int n;
x=(Type *)calloc(r*c,sizeof(Type));
y=(Type **)calloc(r,sizeof(Type *));
for(n=0;n<=r-1;++n)
y=&x;
return(y);
}
template<class Type>
void CMatrix<Type>::MatrixFree(Type **x)
{
free(x);
free(x);
}
template<class Type>
BOOL CMatrix<Type>::Malloc(int m,int n)
{
f=MatrixAlloc(m,n);
if(f==NULL) return FALSE;
row=m;
col=n;
return TRUE;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::Turn()
{
CMatrix<Type> tmp(col,row);
for(int i=0;i<row;i++)
for(int j=0;j<col;j++)
tmp.f=f;
return tmp;
}
template<class Type>
BOOL CMatrix<Type>::Diag(Type *array,int m)
{
f=MatrixAlloc(m,m);
for(int i=0;i<m;i++)
for(int j=0;j<m;j++)
if(i==j) f=array;
row=m;
col=m;
return TRUE;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::Construct(CMatrix<Type> &mat,CString style)
{
int i,j;
CMatrix<Type> tmp;
if(style=="LR"||style=="lr")
{
if(row!=mat.row) return tmp;
if(!tmp.Malloc(row,col+mat.col)) return tmp;
for(i=0;i<tmp.row;i++)
{
for(j=0;j<tmp.col;j++)
{
if(j<col) tmp.f=f;
else tmp.f=mat.f;
}
}
}
return tmp;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::GetChild(int sr,int sc,int er,int ec)
{
int i,j;
CMatrix<Type> tmp(er-sr+1,ec-sc+1);
for(i=0;i<tmp.row;i++)
for(j=0;j<tmp.col;j++)
tmp=f;
return tmp;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::Inverse()
{
Type d,p;
int *is,*js,i,j,k,v;
if(row!=col)
{
cout<<"\nrow!=column,this matrix can't be inversed";
return *this;
}
is=new int;
js=new int;
for(k=0;k<=row-1;k++)
{
d=0.0;
for(i=k;i<=row-1;i++)
for(j=k;j<=row-1;j++)
{
p=fabs(f);
if(p>d)
{
d=p;
is=i;
js=j;
}
}
if(d+1.0==1.0) //singular
{
delete is,js;
cerr<<"\nerror*****not inv,be careful your matrix had been changed\n";
return *this;
}
if(is!=k)
for(j=0;j<=row-1;j++)
{
v=is;
p=f;
f=f;
f=p;
}
if(js!=k)
for(i=0;i<=row-1;i++)
{
v=js;
p=f;
f=f;
f=p;
}
f=1.0/f;
for(j=0;j<=row-1;j++)
if(j!=k)
f*=f;
for(i=0;i<=row-1;i++)
if(i!=k)
for(j=0;j<=row-1;j++)
if(j!=k)
f-=f*f;
for(i=0;i<=row-1;i++)
if(i!=k)
f=-f*f;
}
// change row and column after inverse
for(k=row-1;k>=0;k--)
{
if(js!=k)
for(j=0;j<=row-1;j++)
{
v=js;
p=f;
f=f;
f=p;
}
if(is!=k)
for(i=0;i<=row-1;i++)
{
v=is;
p=f;
f=f;
f=p;
}
}
delete is,js;
return *this;
}
template<class Type>
CMatrix<Type> operator * (double num,CMatrix<Type> right)
{
CMatrix<Type> tmp(right.row,right.col);
int i,j;
for(i=0;i<right.row;i++)
for(j=0;j<right.col;j++)
tmp.f=num*right.f;
return tmp;
}
template<class Type>
CMatrix<Type> operator * (CMatrix<Type> left,double num)
{
CMatrix<Type> tmp(left.row,left.col);
int i,j;
for(i=0;i<left.row;i++)
for(j=0;j<left.col;j++)
tmp.f=num*left.f;
return tmp;
}
template<class Type>
void CMatrix<Type>::Hessenberg()
{
int i,j,k;
double d,t;
if(row==0)
{
cerr<<"\na null matrix,can't use hessenberg transformation";
exit(1);
}
if(row!=col)
{
cerr<<"\nnot a square matrix,can't use hessenberg transformation";
exit(1);
}
for(k=1;k<=row-2;k++)
{
d=0.0;
for(j=k;j<=row-1;j++)
{
t=f;
if(fabs(t)>fabs(d)) { d=t;i=j;}
}
if(fabs(d)+1.0!=1.0)
{
if(i!=k)
{ for(j=k-1;j<=row-1;j++)
{ t=f;f=f;f=t;
}
for(j=0;j<=row-1;j++)
{ t=f;f=f;f=t;
}
}
for(i=k+1;i<=row-1;i++)
{ t=f/d;
f=0.0;
for(j=k;j<=row-1;j++)
{ f-=t*f;
}
for(j=0;j<=row-1;j++)
{ f+=t*f;
}
}
}
}
}
template<class Type>
BOOL CMatrix<Type>::Cholesky()
{
int i,j,k;
double d;
if(row!=col)
{ cerr<<"\nnot a squre matrix, can't use Cholesky disintegration";
return FALSE;
}
if((f+1.0==1.0)||(f<0.0))
{ cerr<<"\nnot a Zhengdin matrix, can't use Cholesky disintegration";
return FALSE;
}
f=sqrt(f);
d=f;
for(i=1;i<=row-1;i++)
{ f/=f;
}
for(j=1;j<=row-1;j++)
{ for(k=0;k<=j-1;k++)
{ f-=f*f;
}
if((f+1.0==1.0)||(f<0.0))
{ cerr<<"\nnot a zhendin matrix, can't use Cholesky disintegration";
return FALSE;
}
f=sqrt(f);
d*=f;
for(i=j+1;i<=row-1;i++)
{ for(k=0;k<=j-1;k++)
f-=f*f;
f/=f;
}
}
for(i=0;i<=row-2;i++)
for(j=i+1;j<=row-1;j++)
f=0.0;
return TRUE;
}
template<class Type>
CMatrix<Type> CMatrix<Type>::Jacobi(double eps,int jt)
{
CMatrix<Type> v;
v.Zeros(row,col);
if(row!=col) return v;
int i,j,p,q,l,n;
n=row;
for(i=0;i<=n-1;i++) // check the matrix's symmetrition
for(j=0;j<=n-1;j++)
if((f-f)>0.001)
return v;
double fm,cn,sn,omega,x,y,d;
l=1;
for (i=0; i<=n-1; i++)
{
v=1.0;
for (j=0; j<=n-1; j++)
if (i!=j) v=0.0;
}
while (1==1)//????????????????????????????????????
{ fm=0.0;
for (i=0; i<=n-1; i++)
for (j=0; j<=n-1; j++)
{ d=fabs(f);
if ((i!=j)&&(d>fm))
{ fm=d; p=i; q=j;}
}
if (fm<eps)return v;
if (l>jt)return v;
l=l+1;
x=-f; y=(f-f)/2.0;
omega=x/sqrt(x*x+y*y);
if (y<0.0) omega=-omega;
sn=1.0+sqrt(1.0-omega*omega);
sn=omega/sqrt(2.0*sn);
cn=sqrt(1.0-sn*sn);
fm=f;
f=fm*cn*cn+f*sn*sn+f*omega;
f=fm*sn*sn+f*cn*cn-f*omega;
f=0.0; f=0.0;
for (j=0; j<=n-1; j++)
if ((j!=p)&&(j!=q))
{
fm=f;
f=fm*cn+f*sn;
f=-fm*sn+f*cn;
}
for (i=0; i<=n-1; i++)
if ((i!=p)&&(i!=q))
{
fm=f;
f=fm*cn+f*sn;
f=-fm*sn+f*cn;
}
for (i=0; i<=n-1; i++)
{
fm=v;
v=fm*cn+v*sn;
v=-fm*sn+v*cn;
}
}
return v;
}
template<class Type>
CMatrix<double> CMatrix<Type>::QREigen(double eps,int jt)
{ // return eigenvalue in two row matrix of a real matrix
int m,it,i,j,k,l,ii,jj,kk,ll;
CMatrix<double> uv;
uv.Zeros(row,2);
if(row!=col) return uv;
int n=row;
double b,c,w,g,xy,p,q,r,x,s,e,ff,z,y;
double *u,*v,*a;
u=(double*)calloc(n,sizeof(double));
v=(double*)calloc(n,sizeof(double));
a=(double*)calloc(n*n,sizeof(double));
this->Hessenberg();
for(i=0;i<n;i++)
for(j=0;j<n;j++)
a=f;
it=0; m=n;
while (m!=0)
{ l=m-1;
while ((l>0)&&(fabs(a)>eps*
(fabs(a[(l-1)*n+l-1])+fabs(a)))) l=l-1;
ii=(m-1)*n+m-1; jj=(m-1)*n+m-2;
kk=(m-2)*n+m-1; ll=(m-2)*n+m-2;
if (l==m-1)
{ u=a[(m-1)*n+m-1]; v=0.0;
m=m-1; it=0;
}
else if (l==m-2)
{ b=-(a+a);
c=a*a-a*a;
w=b*b-4.0*c;
y=sqrt(fabs(w));
if (w>0.0)
{ xy=1.0;
if (b<0.0) xy=-1.0;
u=(-b-xy*y)/2.0;
u=c/u;
v=0.0; v=0.0;
}
else
{ u=-b/2.0; u=u;
v=y/2.0; v=-v;
}
m=m-2; it=0;
}
else
{ if (it>=jt)
{ cout<<"fail\n";
return uv;
}
it=it+1;
for (j=l+2; j<=m-1; j++)
a=0.0;
for (j=l+3; j<=m-1; j++)
a=0.0;
for (k=l; k<=m-2; k++)
{ if (k!=l)
{ p=a; q=a[(k+1)*n+k-1];
r=0.0;
if (k!=m-2) r=a[(k+2)*n+k-1];
}
else
{ x=a+a;
y=a*a-a*a;
ii=l*n+l; jj=l*n+l+1;
kk=(l+1)*n+l; ll=(l+1)*n+l+1;
p=a*(a-x)+a*a+y;
q=a*(a+a-x);
r=a*a[(l+2)*n+l+1];
}
if ((fabs(p)+fabs(q)+fabs(r))!=0.0)
{ xy=1.0;
if (p<0.0) xy=-1.0;
s=xy*sqrt(p*p+q*q+r*r);
if (k!=l) a=-s;
e=-q/s; ff=-r/s; x=-p/s;
y=-x-ff*r/(p+s);
g=e*r/(p+s);
z=-x-e*q/(p+s);
for (j=k; j<=m-1; j++)
{ ii=k*n+j; jj=(k+1)*n+j;
p=x*a+e*a;
q=e*a+y*a;
r=ff*a+g*a;
if (k!=m-2)
{ kk=(k+2)*n+j;
p=p+ff*a;
q=q+g*a;
r=r+z*a; a=r;
}
a=q; a=p;
}
j=k+3;
if (j>=m-1) j=m-1;
for (i=l; i<=j; i++)
{ ii=i*n+k; jj=i*n+k+1;
p=x*a+e*a;
q=e*a+y*a;
r=ff*a+g*a;
if (k!=m-2)
{ kk=i*n+k+2;
p=p+ff*a;
q=q+g*a;
r=r+z*a; a=r;
}
a=q; a=p;
}
}
}
}
}
for(i=0;i<n;i++)
uv=u;
for(i=0;i<n;i++)
uv=v;
return uv;
} text.cpp如下：
#include "matrix.h"
#include <math.h>
#include <fstream.h>
#include <malloc.h>
const unsigned char Z_RADIANS = 0;
const unsigned char Z_DEGREES = 1;
const unsigned char Z_COMMA = 0; // (x, y)
const unsigned char Z_LETTER= 1; // x + iy
const double M_PI= 1.1415926;
class complex
{
public:
double re, im;
private:
static unsigned char zArgMode;
static unsigned char zPrintMode;
static unsigned char zLetter;
public:
complex(void): re(0), im(0) {}
complex(const double real, const double imag=0): re(real), im(imag) {}
complex(const complex& z): re(z.re), im(z.im) {}
friend double re(const complex& z) { // real part
return z.re;
}
friend double im(const complex& z) { // imaginary part
return z.im;
}
friend doublereal(const complex& z) { // real part
return z.re;
}
friend doubleimag(const complex& z) { // imaginary part
return z.im;
}
friend double mag(const complex& z) { // magnitude |z|
return sqrt(z.re*z.re + z.im*z.im);
}
friend double arg(const complex& z); // argument
friend double ang(const complex& z) { // angle
return arg(z);
}
friend double ph(const complex& z) { // phase
return arg(z);
}
friend complex conj(const complex& z) { // complex conjugate
return complex(z.re, -z.im);
}
friend doublenorm(const complex& z) { // norm
return z.re*z.re + z.im*z.im;
}
friend complex rtop(double x, double y=0);
friend complex ptor(double mag, double angle=0);
complex& topolar(void);
complex& torect(void);
void operator = (const complex& z) { // z1 = z2
re = z.re;
im = z.im;
}
complex& operator += (const complex& z) { // z1 += z2
re += z.re;
im += z.im;
return *this;
}
complex& operator -= (const complex& z) { // z1 -= z2
re -= z.re;
im -= z.im;
return *this;
}
complex& operator *= (const complex& z) { // z1 *= z2
*this = *this * z;
return *this;
}
complex& operator /= (const complex& z) { // z1 /= z2
*this = *this / z;
return *this;
}
complex operator + (void) const { // +z1
return *this;
}
complex operator - (void) const { // -z1
return complex(-re, -im);
}
friend complex operator + (const complex& z1, const complex& z2) {
return complex(z1.re + z2.re, z1.im + z2.im);
}
friend complex operator + (const complex& z, const double x) {
return complex(z.re+x, z.im);
}
friend complex operator + (const double x, const complex& z) {
return complex(z.re+x, z.im);
}
friend complex operator - (const complex& z1, const complex& z2) {
return complex(z1.re - z2.re, z1.im - z2.im);
}
friend complex operator - (const complex&, const double);
friend complex operator - (const double x, const complex& z) {
return complex(x-z.re, -z.im);
}
friend complex operator * (const complex& z1, const complex& z2) {
double re = z1.re*z2.re - z1.im*z2.im;
double im = z1.re*z2.im + z1.im*z2.re;
return complex(re, im);
}
friend complex operator * (const complex& z, const double x) {
return complex(z.re*x, z.im*x);
}
friend complex operator * (const double x, const complex& z) {
return complex(z.re*x, z.im*x);
}
friend complex operator / (const complex&, const complex&);
friend complex operator / (const complex& z, const double x) {
return complex(z.re/x, z.im/x);
}
friend complex operator / (const double, const complex&);
friend complex operator ^ (const complex& z1, const complex& z2) {
return pow(z1, z2);
}
friend int operator == (const complex& z1, const complex& z2) {
return (z1.re == z2.re) && (z1.im == z2.im);
}
friend int operator != (const complex& z1, const complex& z2) {
return (z1.re != z2.re) || (z1.im != z2.im);
}
friend double abs(const complex& z);
friend complex sqrt(const complex& z);
friend complex pow(const complex& base, const complex& exp);
friend complex pow(const complex& base, const double exp);
friend complex pow(const double base, const complex& exp);
friend complex exp(const complex& z);
friend complex log(const complex& z);
friend complex ln(const complex& z);
friend complex log10(const complex& z);
friend complexcos(const complex& z);
friend complexsin(const complex& z);
friend complextan(const complex& z);
friend complex asin(const complex& z);
friend complex acos(const complex& z);
friend complex atan(const complex& z);
friend complex sinh(const complex& z);
friend complex cosh(const complex& z);
friend complex tanh(const complex& z);
void SetArgMode(unsigned char mode) const {
if(mode == Z_RADIANS || mode == Z_DEGREES)
zArgMode = mode;
}
void SetPrintMode(unsigned char mode) const {
if(mode == Z_COMMA || mode == Z_LETTER)
zPrintMode = mode;
}
void SetLetter(unsigned char letter) const {
zLetter = letter;
}
friend ostream& operator<<(ostream&, const complex&);
friend istream& operator>>(istream&, const complex&);
};
static const complex Z0(0, 0); // complex number 0
static const complex Z1(1, 0); // complex number 1
static const complex Zi(0, 1); // complex number i
static const complex Zinf(HUGE_VAL, HUGE_VAL); // complex number infinity
static const complex Complex;
///////////////////////////////////////////////
/////////////////////////////////////////////////
///////////////////////////////////////////////
unsigned char complex::zArgMode = Z_RADIANS;
unsigned char complex::zPrintMode = Z_COMMA;
unsigned char complex::zLetter = 'i';
/* -------------------------------------------------------------------- */
double arg(const complex& z) // argument (angle)
{
if(z.re == 0 && z.im == 0) // this is actually a domain error
return 0;
if(complex::zArgMode == Z_RADIANS)
return atan2(z.im, z.re);
return atan2(z.im, z.re)/M_PI*180;
}
complex ptor(double mag, double angle)// polar-to-rectangular
{
if(complex::zArgMode == Z_RADIANS)
return complex(mag*cos(angle), mag*sin(angle));
return complex(mag*cos(angle/180*M_PI), mag*sin(angle/180*M_PI));
}
complex rtop(double x, double y) // rectangular-to-polar
{
if(x == 0 && y == 0) // this is actually a domain error
return Z0;
if(complex::zArgMode == Z_RADIANS)
return complex(sqrt(x*x + y*y), atan2(y, x));
return complex(sqrt(x*x + y*y), atan2(y, x)*180/M_PI);
}
complex& complex::topolar(void)
{
double re_tmp = re;
if(re != 0 || im != 0) // z = (0,0) is a domain error
{
re = sqrt(re*re + im*im);
im = atan2(im, re_tmp);
}
if(complex::zArgMode == Z_DEGREES)
im *= (180/M_PI);
return *this;
}
complex& complex::torect(void)
{
double re_tmp = re;
re = re_tmp*cos(im);
im = re_tmp*sin(im);
return *this;
}
/* ----- Operators ---------------------------------------------------- */
complex operator/(const complex& z1, const complex& z2)
{
if(z2 == Z0)
return complex(Zinf); // z2 = Z0 is an error!
double denom = z2.re*z2.re + z2.im*z2.im;
double re = (z1.re*z2.re + z1.im*z2.im)/denom;
double im = (z2.re*z1.im - z2.im*z1.re)/denom;
return complex(re, im);
}
complex operator/(const double x, const complex& z)
{
if(z == Z0)
return complex(Zinf); // z = Z0 is an error!
double denom = z.re*z.re + z.im*z.im;
return complex(x*z.re/denom, -z.im*x/denom);
}
/* ----- Math functions ----------------------------------------------- */
double abs(const complex& z)
{
return sqrt(z.re*z.re + z.im*z.im);
}
complex sqrt(const complex& z)
{
return ptor(sqrt(abs(z)), arg(z)/2);
}
complex pow(const complex& base, const double exponent)
{
if(base != Z0 && exponent == 0.0)
return complex(1,0);
if (base == Z0 && exponent > 0)
return Z0;
// base == Z0 && exponent == 0 is undefined!
return ptor(pow(abs(base), exponent), exponent*arg(base));
}
complex pow(const double base, const complex& exponent)
{
if(base != 0.0 && exponent == Z0)
return complex(1,0);
if (base == 0 && re(exponent) > 0)
return complex(0,0);
// base == 0 && re(exponent) == 0 is undefined!
if(base > 0.0)
return exp(exponent * log(fabs(base)));
return exp(exponent * complex(log(fabs(base)), M_PI));
}
complex pow(const complex& base, const complex& exponent)
{
if(base != Z0 && exponent == Z0)
return complex(1,0);
if(base == Z0 && re(exponent) > 0)
return complex(0,0);
// base == Z0 && re(exponent) == 0 is undefined!
return exp(exponent * log(base));
}
/* ----- Trig functions ----------------------------------------------- */
complex exp(const complex& z)
{
double mag = exp(z.re);
return complex(mag*cos(z.im), mag*sin(z.im));
}
complex log(const complex& z)
{
return complex(log(mag(z)), atan2(z.im, z.re));
}
/*
complex log10(const complex& z)
{
return log(z) * M_LOG10E;
}
*/
complex sin(const complex& z)
{
return (exp(Zi*z) - exp(-Zi*z))/(2*Zi);
}
complex cos(const complex& z)
{
return (exp(Zi*z) + exp(-Zi*z))/2;
}
complex tan(const complex& z)
{
return sin(z)/cos(z);
}
complex asin(const complex& z)
{
return -Zi*log(Zi*z+sqrt(1-z*z));
}
complex acos(const complex& z)
{
return -Zi*log(z+Zi*sqrt(1-z*z));
}
complex atan(const complex& z)
{
return -0.5*Zi * log((1+Zi*z)/(1-Zi*z));
}
complex sinh(const complex& z)
{
return (exp(z) - exp(-z))/2;
}
complex cosh(const complex& z)
{
return (exp(z) + exp(-z))/2;
}
complex tanh(const complex& z)
{
return sinh(z)/cosh(z);
}
/* ----- Stream I/O --------------------------------------------------- */
ostream& operator<<(ostream& s, const complex& z)
{
char sign[] = " ";
if(complex::zPrintMode == Z_COMMA)
return s << "(" << z.re << ", " << z.im << ")";
if(z.im == 0 || z.im/fabs(z.im) == 1)
sign = '+';
else
sign = '-';
return s << z.re << sign << complex::zLetter << fabs(z.im);
}
istream& operator>>(istream& s, const complex& z)
{
char ch;
s >> ch;
if(ch == '(')
{
s >> z.re >> ch;
if(ch == ',')
s >> z.im >> ch;
if(ch != ')')
s.clear(ios::badbit | s.rdstate());
}
else
{
s.putback(ch);
s >> z.re;
}
return s;
}
///////////////////////////////////////
////////////////////fft
void Fft(complex CA[],int n,int k,complex CB[],int l,int il)
{
int it,m,is,i,j,nv,l0;
double p,q,s,vr,vi,poddr,poddi;
double *pr,*pi,*fr,*fi;
pr=(double*)malloc(n*sizeof(double));
pi=(double*)malloc(n*sizeof(double));
fr=(double*)malloc(n*sizeof(double));
fi=(double*)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{
pr=CA.re;
pi=CA.im;
}
for (it=0; it<=n-1; it++)
{ m=it; is=0;
for (i=0; i<=k-1; i++)
{ j=m/2; is=2*is+(m-2*j); m=j;}
fr=pr; fi=pi;
}
pr=1.0; pi=0.0;
p=6.283185306/(1.0*n);
pr=cos(p); pi=-sin(p);
if (l!=0) pi=-pi;
for (i=2; i<=n-1; i++)
{ p=pr*pr; q=pi*pi;
s=(pr+pi)*(pr+pi);
pr=p-q; pi=s-p-q;
}
for (it=0; it<=n-2; it=it+2)
{ vr=fr; vi=fi;
fr=vr+fr; fi=vi+fi;
fr=vr-fr; fi=vi-fi;
}
m=n/2; nv=2;
for (l0=k-2; l0>=0; l0--)
{ m=m/2; nv=2*nv;
for (it=0; it<=(m-1)*nv; it=it+nv)
for (j=0; j<=(nv/2)-1; j++)
{ p=pr*fr;
q=pi*fi;
s=pr+pi;
s=s*(fr+fi);
poddr=p-q; poddi=s-p-q;
fr=fr-poddr;
fi=fi-poddi;
fr=fr+poddr;
fi=fi+poddi;
}
}
if (l!=0)
for (i=0; i<=n-1; i++)
{ fr=fr/(1.0*n);
fi=fi/(1.0*n);
}
if (il!=0)
for (i=0; i<=n-1; i++)
{ pr=sqrt(fr*fr+fi*fi);
if (fabs(fr)<0.000001*fabs(fi))
{ if ((fi*fr)>0) pi=90.0;
else pi=-90.0;
}
else
pi=atan(fi/fr)*360.0/6.283185306;
}
for(i=0;i<n;i++)
{
CA.re=pr;
CA.im=pi;
}
for(i=0;i<n;i++)
{
CB.re=fr;
CB.im=fi;
}
return;
}
/////////////////////////////////////
voidmain()
{
complex a,b;
a.re=2;
a.im=3;
b=a;
b=exp(a);
cout<<b;
int n=64,i,j;
double deltaT=0.1,R;
ifstream ifile("d:\\work\\tj.fft");
ofstream ofile("d:\\work\\tj.out");
/*for(i=0;i<n;i++)
ofile<<exp(-i*deltaT)<<endl;
complex q;
double L;
int j;
for(L=0;L<n;L++)
{
ifile.seekg(ios::beg);
q=0;
for(j=0;j<n;j++)
{
ifile>>R;
q+=R*exp(-2*M_PI*L*Zi*L*j/n);
}
q*=deltaT;
ofile<<q<<endl;
}
CMatrix<int> A;
A.Ones(3);
cout<<A;
CMatrix<complex> B;
B.Ones(3);
cout<<B;
*/
complex CA;
complex p,f;
for(i=0;i<64;i++)
{
ifile>>p;
}
Fft(p,64,6,f,1,1);
for(i=0;i<=15;i++)
{
for(j=0;j<=3;j++)
ofile<<f<<" ";
ofile<<endl;
}
ofile<<endl<<endl<<endl;
for(i=0;i<=15;i++)
{
for(j=0;j<=3;j++)
ofile<<p<<" ";
ofile<<endl;
}
return ;
}; :lol 有人测试过了吗？ 原帖由 assist 于 2007-9-21 09:39 发表 http://www.chinavib.com/forum/images/common/back.gif
有人测试过了吗？
有人测试过，不过会不会因为转贴造成错误，暂时没办法判断
页:
[1]