C++语言的复数与矩阵运算库

风花雪月

总共包含三个文件:
COMPLEX.H
matrix.h
text.cpp
本程序打包下载见本贴4楼


COMPLEX.H如下：

1. /* -------------------------------------------------------------------- */
   
2. /* Z++ Version 1.10               complex.h       Last revised 10/10/92 */
   
3. /*                                                                      */
   
4. /* Complex number class for Turbo C++/Borland C++.                      */
   
5. /* Copyright 1992 by Carl W. Moreland                                   */
   
6. /* -------------------------------------------------------------------- */
   
7. 
   
8. #ifndef COMPLEXdotH
   
9. #define COMPLEXdotH
   
10. 
    
11. #include <math.h>
    
12. #include <iostream.h>
    
13. 
    
14. const unsigned char Z_RADIANS = 0;
    
15. const unsigned char Z_DEGREES = 1;
    
16. const unsigned char Z_COMMA   = 0;        // (x, y)
    
17. const unsigned char Z_LETTER  = 1;        // x + iy
    
18. 
    
19. class complex
    
20. {
    
21. public:
    
22.   double re, im;
    
23. 
    
24. private:
    
25.   static unsigned char zArgMode;
    
26.   static unsigned char zPrintMode;
    
27.   static unsigned char zLetter;
    
28. 
    
29. public:
    
30.   complex(void): re(0), im(0) {}
    
31.   complex(const double real, const double imag=0): re(real), im(imag) {}
    
32.   complex(const complex& z): re(z.re), im(z.im) {}
    
33. 
    
34.   friend double    re(const complex& z) {        // real part
    
35.     return z.re;
    
36.   }
    
37.   friend double    im(const complex& z) {        // imaginary part
    
38.     return z.im;
    
39.   }
    
40.   friend double  real(const complex& z) {        // real part
    
41.     return z.re;
    
42.   }
    
43.   friend double  imag(const complex& z) {        // imaginary part
    
44.     return z.im;
    
45.   }
    
46.   friend double   mag(const complex& z) {        // magnitude |z|
    
47.     return sqrt(z.re*z.re + z.im*z.im);
    
48.   }
    
49.   friend double   arg(const complex& z);        // argument
    
50.   friend double   ang(const complex& z) {        // angle
    
51.     return arg(z);
    
52.   }
    
53.   friend double    ph(const complex& z) {        // phase
    
54.     return arg(z);
    
55.   }
    
56.   friend complex conj(const complex& z) {        // complex conjugate
    
57.     return complex(z.re, -z.im);
    
58.   }
    
59.   friend double  norm(const complex& z) {        // norm
    
60.     return z.re*z.re + z.im*z.im;
    
61.   }
    
62. 
    
63.   friend complex rtop(double x,   double y=0);
    
64.   friend complex ptor(double mag, double angle=0);
    
65.   complex& topolar(void);
    
66.   complex& torect(void);
    
67. 
    
68.   void operator = (const complex& z) {                // z1 = z2
    
69.     re = z.re;
    
70.     im = z.im;
    
71.   }
    
72.   complex& operator += (const complex& z) {        // z1 += z2
    
73.     re += z.re;
    
74.     im += z.im;
    
75.     return *this;
    
76.   }
    
77.   complex& operator -= (const complex& z) {        // z1 -= z2
    
78.     re -= z.re;
    
79.     im -= z.im;
    
80.     return *this;
    
81.   }
    
82.   complex& operator *= (const complex& z) {        // z1 *= z2
    
83.     *this = *this * z;
    
84.     return *this;
    
85.   }
    
86.   complex& operator /= (const complex& z) {        // z1 /= z2
    
87.     *this = *this / z;
    
88.     return *this;
    
89.   }
    
90.   complex operator + (void) const {                // +z1
    
91.     return *this;
    
92.   }
    
93.   complex operator - (void) const {                // -z1
    
94.     return complex(-re, -im);
    
95.   }
    
96. 
    
97.   friend complex operator + (const complex& z1, const complex& z2) {
    
98.     return complex(z1.re + z2.re, z1.im + z2.im);
    
99.   }
    
100.   friend complex operator + (const complex& z, const double x) {
     
101.     return complex(z.re+x, z.im);
     
102.   }
     
103.   friend complex operator + (const double x, const complex& z) {
     
104.     return complex(z.re+x, z.im);
     
105.   }
     
106.   friend complex operator - (const complex& z1, const complex& z2) {
     
107.     return complex(z1.re - z2.re, z1.im - z2.im);
     
108.   }
     
109.   friend complex operator - (const complex&, const double);
     
110.   friend complex operator - (const double x, const complex& z) {
     
111.     return complex(x-z.re, -z.im);
     
112.   }
     
113.   friend complex operator * (const complex& z1, const complex& z2) {
     
114.     double re = z1.re*z2.re - z1.im*z2.im;
     
115.     double im = z1.re*z2.im + z1.im*z2.re;
     
116.     return complex(re, im);
     
117.   }
     
118.   friend complex operator * (const complex& z, const double x) {
     
119.     return complex(z.re*x, z.im*x);
     
120.   }
     
121.   friend complex operator * (const double x, const complex& z) {
     
122.     return complex(z.re*x, z.im*x);
     
123.   }
     
124.   friend complex operator / (const complex&, const complex&);
     
125.   friend complex operator / (const complex& z, const double x) {
     
126.     return complex(z.re/x, z.im/x);
     
127.   }
     
128.   friend complex operator / (const double, const complex&);
     
129.   friend complex operator ^ (const complex& z1, const complex& z2) {
     
130.     return pow(z1, z2);
     
131.   }
     
132. 
     
133.   friend int operator == (const complex& z1, const complex& z2) {
     
134.     return (z1.re == z2.re) && (z1.im == z2.im);
     
135.   }
     
136.   friend int operator != (const complex& z1, const complex& z2) {
     
137.     return (z1.re != z2.re) || (z1.im != z2.im);
     
138.   }
     
139. 
     
140.   friend double   abs(const complex& z);
     
141.   friend complex sqrt(const complex& z);
     
142.   friend complex pow(const complex& base, const complex& exp);
     
143.   friend complex pow(const complex& base, const double   exp);
     
144.   friend complex pow(const double   base, const complex& exp);
     
145. 
     
146.   friend complex   exp(const complex& z);
     
147.   friend complex   log(const complex& z);
     
148.   friend complex    ln(const complex& z);
     
149.   friend complex log10(const complex& z);
     
150. 
     
151.   friend complex  cos(const complex& z);
     
152.   friend complex  sin(const complex& z);
     
153.   friend complex  tan(const complex& z);
     
154. 
     
155.   friend complex asin(const complex& z);
     
156.   friend complex acos(const complex& z);
     
157.   friend complex atan(const complex& z);
     
158. 
     
159.   friend complex sinh(const complex& z);
     
160.   friend complex cosh(const complex& z);
     
161.   friend complex tanh(const complex& z);
     
162. 
     
163.   void SetArgMode(unsigned char mode) const {
     
164.     if(mode == Z_RADIANS || mode == Z_DEGREES)
     
165.       zArgMode = mode;
     
166.   }
     
167.   void SetPrintMode(unsigned char mode) const {
     
168.     if(mode == Z_COMMA || mode == Z_LETTER)
     
169.       zPrintMode = mode;
     
170.   }
     
171.   void SetLetter(unsigned char letter) const {
     
172.     zLetter = letter;
     
173.   }
     
174. 
     
175.   friend ostream& operator<<(ostream&, const complex&);
     
176.   friend istream& operator>>(istream&, const complex&);
     
177. };
     
178. 
     
179. static const complex Z0(0, 0);                // complex number 0
     
180. static const complex Z1(1, 0);                // complex number 1
     
181. static const complex Zi(0, 1);                // complex number i
     
182. static const complex Zinf(HUGE_VAL, HUGE_VAL); // complex number infinity
     
183. static const complex Complex;
     
184. 
     
185. /* -------------------------------------------------------------------- */
     
186. /* Here is the same class with the name capitalized. If you prefer this */
     
187. /* simply remove the comment delimiters below and comment out the 5     */
     
188. /* static globals above.                                                */
     
189. /* -------------------------------------------------------------------- */
     
190. /*
     
191. class Complex: public complex
     
192. {
     
193. public:
     
194.   Complex(void): re(0), im(0) {}
     
195.   Complex(const double real, const double imag=0): re(real), im(imag) {}
     
196.   Complex(const complex& z): re(z.re), im(z.im) {}
     
197. };
     
198. 
     
199. static const Complex Z0(0, 0);                // complex number 0
     
200. static const Complex Z1(1, 0);                // complex number 1
     
201. static const Complex Zi(0, 1);                // complex number i
     
202. static const Complex Zinf(HUGE_VAL, HUGE_VAL); // complex number infinity
     
203. static const Complex complex;
     
204. */
     
205. 
     
206. #endif

复制代码

[ 本帖最后由 风花雪月 于 2006-10-10 00:44 编辑 ]

风花雪月

matrix.h如下：

1. #include "afxwin.h"
   
2. #include "strstrea.h"
   
3. #include "fstream.h"
   
4. #include "iomanip.h"
   
5. #include "stdlib.h"
   
6. #include "malloc.h"
   
7. #include "math.h"
   
8. 
   
9. #define CDMatrix (CMatrix<double>)
   
10. #define CIMatrix (CMatrix<int>)
    
11. 
    
12. template<class Type>
    
13. class CMatrix :public CObject
    
14. {
    
15.         protected:               
    
16.                         Type    **f;
    
17.                         int     row;
    
18.                         int     col;
    
19.         public:
    
20.                 int  GetRow(){return row;}
    
21.                 int  GetCol(){return col;}
    
22.                 void SetRow(int r){row=r;}
    
23.                 void SetCol(int l){col=l;}
    
24.                 Type** GetF(){return f;}               
    
25.                 CMatrix();
    
26.                 CMatrix(int,int);
    
27.                 CMatrix(CMatrix &copy);
    
28.                 ~CMatrix();
    
29.                 BOOL Ones(int);
    
30.                 BOOL Ones(int,int);
    
31.                 BOOL Zeros(int);
    
32.                 BOOL Zeros(int,int);
    
33.                 CMatrix Turn();
    
34.                 CMatrix Inverse();
    
35.                 CMatrix Construct(CMatrix &mat,CString style);
    
36.                 CMatrix GetChild(int,int,int,int);
    
37.                 BOOL Diag(Type *,int);
    
38.                 BOOL Malloc(int,int);
    
39.                 CMatrix   operator*(CMatrix  &m);
    
40.                 friend CMatrix  operator*(double,CMatrix);
    
41.                 friend CMatrix  operator*(CMatrix,double);
    
42.                 CMatrix   operator/(Type x);
    
43.                 CMatrix   operator+(CMatrix  &m);
    
44.                 CMatrix   operator+(Type x);
    
45.                 CMatrix   operator-(CMatrix  &m);
    
46.                 CMatrix   operator-(Type x);
    
47.                 CMatrix   operator-();
    
48.                 Type* &operator[](int);
    
49.                 friend istream  &operator>>(istream  &,CMatrix  &);
    
50.                 friend ostream  &operator<<(ostream  &,CMatrix  &);
    
51.                 void operator=(const CMatrix  &m);
    
52.                 void operator=(Type);
    
53.                 Type ** MatrixAlloc(int,int);
    
54.                 void    MatrixFree(Type **);
    
55.                 void    Hessenberg();
    
56.                 BOOL    Cholesky();
    
57.                 CMatrix QREigen(double,int);// return eigenvalue in two row matrix of a real matrix                
    
58.                 CMatrix Jacobi(double eps,int jt);
    
59. };
    
60. 
    
61. template<class Type>
    
62. Type* &CMatrix<Type>::operator[](int i)
    
63. {
    
64.         Type *tp;
    
65.         tp=(Type *)f[i];
    
66.         return tp;
    
67. }
    
68. template<class Type>
    
69. CMatrix<Type>::CMatrix()
    
70. {
    
71.         f=NULL;
    
72.         row=0;
    
73.         col=0;
    
74. }
    
75. template<class Type>
    
76. CMatrix<Type>::CMatrix(int m,int n)
    
77. {       
    
78.         f=MatrixAlloc(m,n);
    
79.         row=m;
    
80.         col=n;  
    
81. }
    
82. template<class Type>
    
83. CMatrix<Type>::~CMatrix()
    
84. {
    
85.         if(f!=NULL) MatrixFree(f);
    
86. }
    
87. 
    
88. template<class Type>
    
89. BOOL CMatrix<Type>::Ones(int m)
    
90. {        
    
91.         f=MatrixAlloc(m,m);
    
92.         for(int i=0;i<m;i++)
    
93.         for(int j=0;j<m;j++)
    
94.         {
    
95.                 if(i!=j)f[i][j]=0;
    
96.                 else f[i][j]=1;
    
97.         }
    
98.         row=m;
    
99.         col=m;
    
100.         return TRUE;
     
101. }
     
102. template<class Type>
     
103. BOOL CMatrix<Type>::Ones(int m,int n)
     
104. {        
     
105.         f=MatrixAlloc(m,n);
     
106.         for(int i=0;i<m;i++)
     
107.         for(int j=0;j<n;j++)
     
108.         {
     
109.                 f[i][j]=1;
     
110.         }
     
111.         row=m;
     
112.         col=n;
     
113.         return TRUE;
     
114. }
     
115. template<class Type>
     
116. BOOL CMatrix<Type>::Zeros(int m)
     
117. {        
     
118.         f=MatrixAlloc(m,m);
     
119.         for(int i=0;i<m;i++)
     
120.         for(int j=0;j<m;j++)
     
121.         {
     
122.                 f[i][j]=0;
     
123.         }
     
124.         row=m;
     
125.         col=m;
     
126.         return TRUE;
     
127. }
     
128. template<class Type>
     
129. BOOL CMatrix<Type>::Zeros(int m,int n)
     
130. {        
     
131.         f=MatrixAlloc(m,n);
     
132.         if(f==NULL) return FALSE;
     
133.         for(int i=0;i<m;i++)
     
134.         for(int j=0;j<n;j++)
     
135.         {
     
136.                 f[i][j]=0;
     
137.         }
     
138.         row=m;
     
139.         col=n;
     
140.         return TRUE;
     
141. }
     
142. 
     
143. template<class Type>
     
144. CMatrix<Type> CMatrix<Type>::operator*(CMatrix<Type>  &mat)
     
145. {
     
146.         int l,m,n;
     
147.         l=GetRow();
     
148.         m=GetCol();
     
149.         n=mat.GetCol();
     
150.         CMatrix tmp(l,n);
     
151.         for(int k=0;k<l;k++)
     
152.         for(int j=0;j<n;j++)
     
153.         {
     
154.                 for(int i=0;i<m;i++)
     
155.                 tmp.f[k][j]+=f[k][i]*mat.f[i][j];
     
156.         }       
     
157.         return tmp;
     
158. }
     
159. 
     
160. template<class Type>
     
161. CMatrix<Type> CMatrix<Type>::operator/(Type x)
     
162. {
     
163.         int l,n;
     
164.         l=GetRow();
     
165.         n=GetCol();
     
166.         CMatrix<Type> tmp(l,n);
     
167.         for(int j=0;j<l;j++)
     
168.         for(int i=0;i<n;i++)
     
169.                 tmp.f[j][i]=f[j][i]/x;
     
170.         return tmp;
     
171. }
     
172. 
     
173. template<class Type>
     
174. CMatrix<Type> CMatrix<Type>::operator+( CMatrix<Type>  &mat)
     
175. {
     
176.         int l1,m1;
     
177.         l1=GetRow();
     
178.         m1=GetCol();
     
179.         CMatrix<Type> tmp(l1,m1);
     
180.         for(int j=0;j<l1;j++)
     
181.         for(int i=0;i<m1;i++)
     
182.                 tmp.f[j][i]=f[j][i]+mat.f[j][i];
     
183.         return tmp;
     
184. }
     
185. 
     
186. template<class Type>
     
187. CMatrix<Type> CMatrix<Type>::operator+(Type x)
     
188. {
     
189.         int l,m;
     
190.         l=GetRow();
     
191.         m=GetCol();
     
192.         CMatrix<Type> tmp(l,m);
     
193.         for(int j=0;j<l;j++)
     
194.         for(int i=0;i<m;i++)
     
195.                 tmp.f[j][i]=f[j][i]+x;
     
196.         return tmp;
     
197. }
     
198. 
     
199. template<class Type>
     
200. CMatrix<Type> CMatrix<Type>::operator-(CMatrix<Type>  &mat)
     
201. {
     
202.         int l1,m1;
     
203.         l1=GetRow();
     
204.         m1=GetCol();
     
205.         CMatrix<Type> tmp(l1,m1);
     
206.         for(int j=0;j<l1;j++)
     
207.         for(int i=0;i<m1;i++)
     
208.                 tmp.f[j][i]=f[j][i]-mat.f[j][i];
     
209.         return tmp;
     
210. }
     
211. template<class Type>
     
212. CMatrix<Type> CMatrix<Type>::operator-()
     
213. {
     
214.         for(int i=0;i<row;i++)
     
215.         for(int j=0;j<col;j++)
     
216.                 f[i][j]=-f[i][j];
     
217.         return *this;
     
218. }       
     
219. 
     
220. template<class Type>
     
221. CMatrix<Type> CMatrix<Type>::operator-( Type x)
     
222. {
     
223.         int l,m;
     
224.         l=GetRow();
     
225.         m=GetCol();
     
226.         CMatrix<Type> temp(l,m);
     
227.         for(int j=0;j<l;j++)
     
228.         for(int i=0;i<m;i++)
     
229.                 temp.f[j][i]=f[j][i]-x;
     
230.         return temp;
     
231. }
     
232. 
     
233. template<class Type>
     
234. istream &operator>>(istream &in,CMatrix<Type> &mat)
     
235. {   
     
236.         int i,j;
     
237.         for(i=0;i<mat.row;i++)
     
238.         for(j=0;j<mat.col;j++)
     
239.                  in>>mat.f[i][j];
     
240.         cout<<endl;
     
241.         return in;
     
242. };
     
243. 
     
244. template<class Type>
     
245. ostream  & operator<<(ostream& out,CMatrix<Type> &v1)
     
246. {
     
247.         if(v1.GetF()==NULL)
     
248.         {
     
249.                 out<<"This Matrix cannot be output!";
     
250.                 return out<<endl;
     
251.         }
     
252.         out<<setiosflags(ios::right||ios::fixed);
     
253.         out<<setprecision(5);
     
254.         out<<"["<<endl;
     
255.         int mr=v1.GetRow();
     
256.         int mc=v1.GetCol();
     
257.         for(int j=0;j<mr;j++)
     
258.         {
     
259.                 for(int i=0;i<mc-1;i++)
     
260.                         out<<setw(12)<<v1.f[j][i];
     
261.                 out<<setw(12)<<v1.f[j][mc-1]<<";"<<endl;
     
262.         }
     
263.         return out<<"]"<<endl;
     
264.        
     
265. }
     
266. template<class Type>
     
267. CMatrix<Type>::CMatrix(CMatrix<Type> &copy)
     
268. {
     
269.         f=NULL;
     
270.         *this=copy;
     
271. }
     
272.        
     
273. template<class Type>
     
274. void CMatrix<Type>::operator=(const CMatrix<Type> &copy)
     
275. {
     
276.        
     
277.         if(this==&copy) return;
     
278.         if(f!=NULL) MatrixFree(f);
     
279.         int m,n;
     
280.         m=copy.row;
     
281.         n=copy.col;
     
282.         f=MatrixAlloc(m,n);
     
283.           for(int i=0;i<m;i++)
     
284.         for(int j=0;j<n;j++)
     
285.                 f[i][j]=copy.f[i][j];
     
286.         row=m;
     
287.         col=n;  
     
288. }
     
289. 
     
290. template<class Type>
     
291. void CMatrix<Type>::operator=(Type d)
     
292. {
     
293.        
     
294.         if(f==NULL)
     
295.         {
     
296.                 cout<<"The Matrix is not be allociated"<<endl;
     
297.                 return;
     
298.         }
     
299.         for(int i=0;i<row;i++)
     
300.         for(int j=0;j<col;j++)
     
301.                 f[i][j]=d;
     
302. }
     
303. 
     
304. template<class Type>
     
305. Type ** CMatrix<Type>::MatrixAlloc(int r,int c)
     
306. {
     
307.     Type *x,**y;
     
308.     int   n;
     
309.     x=(Type *)calloc(r*c,sizeof(Type));
     
310.     y=(Type **)calloc(r,sizeof(Type *));
     
311.         for(n=0;n<=r-1;++n)
     
312.         y[n]=&x[c*n];
     
313.         return(y);
     
314. }
     
315. 
     
316. template<class Type>
     
317. void CMatrix<Type>::MatrixFree(Type **x)
     
318. {
     
319.    free(x[0]);
     
320.    free(x);
     
321. }
     
322. template<class Type>
     
323. BOOL CMatrix<Type>::Malloc(int m,int n)
     
324. {
     
325.         f=MatrixAlloc(m,n);
     
326.         if(f==NULL) return FALSE;
     
327.         row=m;
     
328.         col=n;
     
329.         return TRUE;
     
330. }
     
331. 
     
332.        
     
333. template<class Type>
     
334. CMatrix<Type> CMatrix<Type>::Turn()
     
335. {
     
336.         CMatrix<Type> tmp(col,row);
     
337.         for(int i=0;i<row;i++)
     
338.         for(int j=0;j<col;j++)
     
339.                 tmp.f[j][i]=f[i][j];
     
340.         return tmp;
     
341. }
     
342. 
     
343. template<class Type>
     
344. BOOL CMatrix<Type>::Diag(Type *array,int m)
     
345. {
     
346.         f=MatrixAlloc(m,m);
     
347.         for(int i=0;i<m;i++)
     
348.         for(int j=0;j<m;j++)
     
349.                 if(i==j) f[i][j]=array[i];
     
350.         row=m;
     
351.         col=m;
     
352.         return TRUE;
     
353. }
     
354. 
     
355. template<class Type>
     
356. CMatrix<Type> CMatrix<Type>::Construct(CMatrix<Type> &mat,CString style)
     
357. {
     
358.         int i,j;
     
359.         CMatrix<Type> tmp;
     
360.         if(style=="LR"||style=="lr")
     
361.         {
     
362.                 if(row!=mat.row)    return tmp;
     
363.                 if(!tmp.Malloc(row,col+mat.col)) return tmp;
     
364.                 for(i=0;i<tmp.row;i++)
     
365.                 {
     
366.                         for(j=0;j<tmp.col;j++)
     
367.                         {
     
368.                                 if(j<col) tmp.f[i][j]=f[i][j];
     
369.                                 else      tmp.f[i][j]=mat.f[i][j-col];
     
370.                         }
     
371.                 }
     
372.         }
     
373.         return tmp;
     
374. }
     
375. 
     
376. template<class Type>       
     
377. CMatrix<Type> CMatrix<Type>::GetChild(int sr,int sc,int er,int ec)
     
378. {
     
379.         int i,j;
     
380.         CMatrix<Type> tmp(er-sr+1,ec-sc+1);
     
381.         for(i=0;i<tmp.row;i++)
     
382.         for(j=0;j<tmp.col;j++)
     
383.                 tmp[i][j]=f[sr+i][sc+j];
     
384.         return tmp;
     
385. }
     
386. template<class Type>       
     
387. CMatrix<Type> CMatrix<Type>::Inverse()
     
388. {  
     
389.         Type d,p;
     
390.     int *is,*js,i,j,k,v;
     
391.     if(row!=col)  
     
392.         {  
     
393.                 cout<<"\nrow!=column,this matrix can't be inversed";
     
394.                 return *this;   
     
395.         }
     
396.     is=new int[row];
     
397.     js=new int[row];
     
398.     for(k=0;k<=row-1;k++)
     
399.         {
     
400.                 d=0.0;
     
401.             for(i=k;i<=row-1;i++)
     
402.             for(j=k;j<=row-1;j++)
     
403.             {  
     
404.                         p=fabs(f[i][j]);
     
405.                     if(p>d)
     
406.                         {  
     
407.                                 d=p;
     
408.                                 is[k]=i;
     
409.                                 js[k]=j;  
     
410.                         }
     
411.                 }
     
412.             if(d+1.0==1.0)     //singular
     
413.             {  
     
414.                         delete is,js;
     
415.                     cerr<<"\nerror*****not inv,be careful your matrix had been changed\n";
     
416.                     return *this;   
     
417.                 }
     
418.             if(is[k]!=k)
     
419.             for(j=0;j<=row-1;j++)
     
420.                 {   
     
421.                         v=is[k];
     
422.                     p=f[k][j];
     
423.                         f[k][j]=f[v][j];
     
424.                         f[v][j]=p;
     
425.                 }
     
426.             if(js[k]!=k)
     
427.             for(i=0;i<=row-1;i++)
     
428.                 {  
     
429.                         v=js[k];
     
430.                     p=f[i][k];
     
431.                         f[i][k]=f[i][v];
     
432.                         f[i][v]=p;
     
433.                 }
     
434.             f[k][k]=1.0/f[k][k];
     
435.             for(j=0;j<=row-1;j++)
     
436.             if(j!=k)
     
437.                f[k][j]*=f[k][k];
     
438.             for(i=0;i<=row-1;i++)
     
439.                         if(i!=k)
     
440.                 for(j=0;j<=row-1;j++)
     
441.                                 if(j!=k)
     
442.                              f[i][j]-=f[i][k]*f[k][j];
     
443.                 for(i=0;i<=row-1;i++)
     
444.             if(i!=k)
     
445.                          f[i][k]=-f[i][k]*f[k][k];
     
446.         }
     
447.     // change row and column after inverse
     
448.         for(k=row-1;k>=0;k--)
     
449.         {  
     
450.                 if(js[k]!=k)
     
451.                 for(j=0;j<=row-1;j++)
     
452.                 {  
     
453.                         v=js[k];
     
454.                     p=f[k][j];
     
455.                         f[k][j]=f[v][j];
     
456.                         f[v][j]=p;
     
457.                 }
     
458.             if(is[k]!=k)
     
459.                          for(i=0;i<=row-1;i++)
     
460.                          {  
     
461.                                  v=is[k];
     
462.                                  p=f[i][k];  
     
463.                                  f[i][k]=f[i][v];
     
464.                                  f[i][v]=p;
     
465.                      }
     
466.       }
     
467.         delete is,js;
     
468.         return *this;   
     
469. }
     
470. 
     
471. template<class Type>       
     
472. CMatrix<Type> operator * (double num,CMatrix<Type> right)
     
473. {
     
474.     CMatrix<Type> tmp(right.row,right.col);
     
475.         int i,j;
     
476.         for(i=0;i<right.row;i++)
     
477.         for(j=0;j<right.col;j++)
     
478.                  tmp.f[i][j]=num*right.f[i][j];
     
479.         return tmp;
     
480. }
     
481. 
     
482. template<class Type>       
     
483. CMatrix<Type> operator * (CMatrix<Type> left,double num)
     
484. {   
     
485.         CMatrix<Type> tmp(left.row,left.col);
     
486.         int i,j;
     
487.         for(i=0;i<left.row;i++)
     
488.         for(j=0;j<left.col;j++)
     
489.                  tmp.f[i][j]=num*left.f[i][j];
     
490.         return tmp;
     
491. }
     
492. template<class Type>       
     
493. void CMatrix<Type>::Hessenberg()
     
494. {
     
495.         int i,j,k;
     
496.         double d,t;
     
497.         if(row==0)
     
498.         {
     
499.             cerr<<"\na null matrix,can't use hessenberg transformation";
     
500.             exit(1);
     
501.          }
     
502.         if(row!=col)
     
503.         {
     
504.             cerr<<"\nnot a square matrix,can't use hessenberg transformation";
     
505.             exit(1);
     
506.          }
     
507.         for(k=1;k<=row-2;k++)
     
508.         {
     
509.             d=0.0;
     
510.             for(j=k;j<=row-1;j++)
     
511.             {
     
512.                  t=f[j][k-1];
     
513.                  if(fabs(t)>fabs(d))      { d=t;  i=j;  }
     
514.              }
     
515.             if(fabs(d)+1.0!=1.0)
     
516.             {
     
517.                  if(i!=k)
     
518.                  {      for(j=k-1;j<=row-1;j++)
     
519.                         {     t=f[i][j];  f[i][j]=f[k][j];  f[k][j]=t;
     
520.                          }
     
521.                         for(j=0;j<=row-1;j++)
     
522.                         {      t=f[j][i];  f[j][i]=f[j][k];  f[j][k]=t;
     
523.                          }
     
524.                   }
     
525.                  for(i=k+1;i<=row-1;i++)
     
526.                  {       t=f[i][k-1]/d;
     
527.                          f[i][k-1]=0.0;
     
528.                          for(j=k;j<=row-1;j++)
     
529.                          {      f[i][j]-=t*f[k][j];
     
530.                           }
     
531.                          for(j=0;j<=row-1;j++)
     
532.                          {      f[j][k]+=t*f[j][i];
     
533.                           }
     
534.                   }
     
535.              }
     
536.         }
     
537. }
     
538. 
     
539. template<class Type>       
     
540. BOOL CMatrix<Type>::Cholesky()
     
541. {
     
542.         int i,j,k;
     
543.         double d;
     
544.         if(row!=col)
     
545.         {        cerr<<"\nnot a squre matrix, can't use Cholesky disintegration";
     
546.                 return FALSE;
     
547.         }
     
548.         if((f[0][0]+1.0==1.0)||(f[0][0]<0.0))
     
549.         {        cerr<<"\nnot a Zhengdin matrix, can't use Cholesky disintegration";
     
550.                 return FALSE;
     
551.         }
     
552.         f[0][0]=sqrt(f[0][0]);
     
553.         d=f[0][0];
     
554.         for(i=1;i<=row-1;i++)
     
555.         {        f[i][0]/=f[0][0];
     
556.         }
     
557.         for(j=1;j<=row-1;j++)
     
558.         {        for(k=0;k<=j-1;k++)
     
559.                 {   f[j][j]-=f[j][k]*f[j][k];
     
560.                 }
     
561.         if((f[j][j]+1.0==1.0)||(f[j][j]<0.0))
     
562.                 {   cerr<<"\nnot a zhendin matrix, can't use Cholesky disintegration";
     
563.                     return FALSE;
     
564.                 }
     
565.                 f[j][j]=sqrt(f[j][j]);
     
566.                 d*=f[j][j];
     
567.             for(i=j+1;i<=row-1;i++)
     
568.                 {        for(k=0;k<=j-1;k++)
     
569.                 f[i][j]-=f[i][k]*f[j][k];
     
570.                     f[i][j]/=f[j][j];
     
571.                 }
     
572.         }
     
573.         for(i=0;i<=row-2;i++)
     
574.                 for(j=i+1;j<=row-1;j++)
     
575.                         f[i][j]=0.0;
     
576.     return TRUE;
     
577. }
     
578. template<class Type>       
     
579. CMatrix<Type> CMatrix<Type>::Jacobi(double eps,int jt)
     
580. {
     
581.                 CMatrix<Type> v;
     
582.                 v.Zeros(row,col);
     
583.                 if(row!=col) return v;
     
584.                 int i,j,p,q,l,n;
     
585.                 n=row;
     
586.                 for(i=0;i<=n-1;i++)    // check the matrix's symmetrition
     
587.                         for(j=0;j<=n-1;j++)
     
588.                                 if((f[i][j]-f[j][i])>0.001)
     
589.                                         return v;
     
590.                 double fm,cn,sn,omega,x,y,d;
     
591.                 l=1;
     
592.                 for (i=0; i<=n-1; i++)
     
593.                 {
     
594.                         v[i][i]=1.0;
     
595.                         for (j=0; j<=n-1; j++)
     
596.                         if (i!=j) v[i][j]=0.0;
     
597.                 }
     
598.     while (1==1)//????????????????????????????????????
     
599.       { fm=0.0;
     
600.         for (i=0; i<=n-1; i++)
     
601.         for (j=0; j<=n-1; j++)
     
602.           { d=fabs(f[i][j]);
     
603.             if ((i!=j)&&(d>fm))
     
604.               { fm=d; p=i; q=j;}
     
605.           }
     
606.         if (fm<eps)  return v;
     
607.         if (l>jt)  return v;
     
608.         l=l+1;
     
609.         x=-f[p][q]; y=(f[q][q]-f[p][p])/2.0;
     
610.         omega=x/sqrt(x*x+y*y);
     
611.         if (y<0.0) omega=-omega;
     
612.         sn=1.0+sqrt(1.0-omega*omega);
     
613.         sn=omega/sqrt(2.0*sn);
     
614.         cn=sqrt(1.0-sn*sn);
     
615.         fm=f[p][p];
     
616.         f[p][p]=fm*cn*cn+f[q][q]*sn*sn+f[p][q]*omega;
     
617.         f[q][q]=fm*sn*sn+f[q][q]*cn*cn-f[p][q]*omega;
     
618.         f[p][q]=0.0; f[q][p]=0.0;
     
619.         for (j=0; j<=n-1; j++)
     
620.         if ((j!=p)&&(j!=q))
     
621.         {
     
622.             fm=f[p][j];
     
623.             f[p][j]=fm*cn+f[q][j]*sn;
     
624.             f[q][j]=-fm*sn+f[q][j]*cn;
     
625.           }
     
626.         for (i=0; i<=n-1; i++)
     
627.           if ((i!=p)&&(i!=q))
     
628.           {
     
629.               fm=f[i][p];
     
630.               f[i][p]=fm*cn+f[i][q]*sn;
     
631.               f[i][q]=-fm*sn+f[i][q]*cn;
     
632.             }
     
633.         for (i=0; i<=n-1; i++)
     
634.           {
     
635.             fm=v[i][p];
     
636.             v[i][p]=fm*cn+v[i][q]*sn;
     
637.             v[i][q]=-fm*sn+v[i][q]*cn;
     
638.           }
     
639.       }
     
640.     return v;
     
641. }
     
642. 
     
643. 
     
644. 
     
645. 
     
646. template<class Type>
     
647. CMatrix<double> CMatrix<Type>::QREigen(double eps,int jt)
     
648. {              // return eigenvalue in two row matrix of a real matrix
     
649.         int m,it,i,j,k,l,ii,jj,kk,ll;
     
650.     CMatrix<double> uv;
     
651.         uv.Zeros(row,2);
     
652.         if(row!=col) return uv;
     
653.         int n=row;
     
654.     double b,c,w,g,xy,p,q,r,x,s,e,ff,z,y;
     
655.     double *u,*v,*a;
     
656.     u=(double*)calloc(n,sizeof(double));
     
657.     v=(double*)calloc(n,sizeof(double));
     
658.     a=(double*)calloc(n*n,sizeof(double));
     
659.     this->Hessenberg();
     
660.         for(i=0;i<n;i++)
     
661.                 for(j=0;j<n;j++)
     
662.                         a[i*n+j]=f[i][j];
     
663.         it=0; m=n;
     
664.     while (m!=0)
     
665.       { l=m-1;
     
666.         while ((l>0)&&(fabs(a[l*n+l-1])>eps*
     
667.               (fabs(a[(l-1)*n+l-1])+fabs(a[l*n+l])))) l=l-1;
     
668.         ii=(m-1)*n+m-1; jj=(m-1)*n+m-2;
     
669.         kk=(m-2)*n+m-1; ll=(m-2)*n+m-2;
     
670.         if (l==m-1)
     
671.           { u[m-1]=a[(m-1)*n+m-1]; v[m-1]=0.0;
     
672.             m=m-1; it=0;
     
673.           }
     
674.         else if (l==m-2)
     
675.           { b=-(a[ii]+a[ll]);
     
676.             c=a[ii]*a[ll]-a[jj]*a[kk];
     
677.             w=b*b-4.0*c;
     
678.             y=sqrt(fabs(w));
     
679.             if (w>0.0)
     
680.               { xy=1.0;
     
681.                 if (b<0.0) xy=-1.0;
     
682.                 u[m-1]=(-b-xy*y)/2.0;
     
683.                 u[m-2]=c/u[m-1];
     
684.                 v[m-1]=0.0; v[m-2]=0.0;
     
685.               }
     
686.             else
     
687.               { u[m-1]=-b/2.0; u[m-2]=u[m-1];
     
688.                 v[m-1]=y/2.0; v[m-2]=-v[m-1];
     
689.               }
     
690.             m=m-2; it=0;
     
691.           }
     
692.         else
     
693.           { if (it>=jt)
     
694.               { cout<<"fail\n";
     
695.                 return uv;
     
696.               }
     
697.             it=it+1;
     
698.             for (j=l+2; j<=m-1; j++)
     
699.               a[j*n+j-2]=0.0;
     
700.             for (j=l+3; j<=m-1; j++)
     
701.               a[j*n+j-3]=0.0;
     
702.             for (k=l; k<=m-2; k++)
     
703.               { if (k!=l)
     
704.                   { p=a[k*n+k-1]; q=a[(k+1)*n+k-1];
     
705.                     r=0.0;
     
706.                     if (k!=m-2) r=a[(k+2)*n+k-1];
     
707.                   }
     
708.                 else
     
709.                   { x=a[ii]+a[ll];
     
710.                     y=a[ll]*a[ii]-a[kk]*a[jj];
     
711.                     ii=l*n+l; jj=l*n+l+1;
     
712.                     kk=(l+1)*n+l; ll=(l+1)*n+l+1;
     
713.                     p=a[ii]*(a[ii]-x)+a[jj]*a[kk]+y;
     
714.                     q=a[kk]*(a[ii]+a[ll]-x);
     
715.                     r=a[kk]*a[(l+2)*n+l+1];
     
716.                   }
     
717.                 if ((fabs(p)+fabs(q)+fabs(r))!=0.0)
     
718.                   { xy=1.0;
     
719.                     if (p<0.0) xy=-1.0;
     
720.                     s=xy*sqrt(p*p+q*q+r*r);
     
721.                     if (k!=l) a[k*n+k-1]=-s;
     
722.                     e=-q/s; ff=-r/s; x=-p/s;
     
723.                     y=-x-ff*r/(p+s);
     
724.                     g=e*r/(p+s);
     
725.                     z=-x-e*q/(p+s);
     
726.                     for (j=k; j<=m-1; j++)
     
727.                       { ii=k*n+j; jj=(k+1)*n+j;
     
728.                         p=x*a[ii]+e*a[jj];
     
729.                         q=e*a[ii]+y*a[jj];
     
730.                         r=ff*a[ii]+g*a[jj];
     
731.                         if (k!=m-2)
     
732.                           { kk=(k+2)*n+j;
     
733.                             p=p+ff*a[kk];
     
734.                             q=q+g*a[kk];
     
735.                             r=r+z*a[kk]; a[kk]=r;
     
736.                           }
     
737.                         a[jj]=q; a[ii]=p;
     
738.                       }
     
739.                     j=k+3;
     
740.                     if (j>=m-1) j=m-1;
     
741.                     for (i=l; i<=j; i++)
     
742.                       { ii=i*n+k; jj=i*n+k+1;
     
743.                         p=x*a[ii]+e*a[jj];
     
744.                         q=e*a[ii]+y*a[jj];
     
745.                         r=ff*a[ii]+g*a[jj];
     
746.                         if (k!=m-2)
     
747.                           { kk=i*n+k+2;
     
748.                             p=p+ff*a[kk];
     
749.                             q=q+g*a[kk];
     
750.                             r=r+z*a[kk]; a[kk]=r;
     
751.                           }
     
752.                         a[jj]=q; a[ii]=p;
     
753.                       }
     
754.                   }
     
755.               }
     
756.           }
     
757.       }
     
758. 
     
759.     for(i=0;i<n;i++)
     
760.                 uv[i][0]=u[i];
     
761.         for(i=0;i<n;i++)
     
762.                 uv[i][1]=v[i];
     
763.         return uv;
     
764.   }

复制代码

风花雪月

text.cpp如下：

1. #include "matrix.h"
   
2. #include <math.h>
   
3. #include <fstream.h>
   
4. #include <malloc.h>
   
5. 
   
6. const unsigned char Z_RADIANS = 0;
   
7. const unsigned char Z_DEGREES = 1;
   
8. const unsigned char Z_COMMA   = 0;        // (x, y)
   
9. const unsigned char Z_LETTER  = 1;        // x + iy
   
10. const double M_PI  = 1.1415926;       
    
11. class complex
    
12. {
    
13. public:
    
14.   double re, im;
    
15. 
    
16. private:
    
17.   static unsigned char zArgMode;
    
18.   static unsigned char zPrintMode;
    
19.   static unsigned char zLetter;
    
20. 
    
21. public:
    
22.   complex(void): re(0), im(0) {}
    
23.   complex(const double real, const double imag=0): re(real), im(imag) {}
    
24.   complex(const complex& z): re(z.re), im(z.im) {}
    
25. 
    
26.   friend double    re(const complex& z) {        // real part
    
27.     return z.re;
    
28.   }
    
29.   friend double    im(const complex& z) {        // imaginary part
    
30.     return z.im;
    
31.   }
    
32.   friend double  real(const complex& z) {        // real part
    
33.     return z.re;
    
34.   }
    
35.   friend double  imag(const complex& z) {        // imaginary part
    
36.     return z.im;
    
37.   }
    
38.   friend double   mag(const complex& z) {        // magnitude |z|
    
39.     return sqrt(z.re*z.re + z.im*z.im);
    
40.   }
    
41.   friend double   arg(const complex& z);        // argument
    
42.   friend double   ang(const complex& z) {        // angle
    
43.     return arg(z);
    
44.   }
    
45.   friend double    ph(const complex& z) {        // phase
    
46.     return arg(z);
    
47.   }
    
48.   friend complex conj(const complex& z) {        // complex conjugate
    
49.     return complex(z.re, -z.im);
    
50.   }
    
51.   friend double  norm(const complex& z) {        // norm
    
52.     return z.re*z.re + z.im*z.im;
    
53.   }
    
54. 
    
55.   friend complex rtop(double x,   double y=0);
    
56.   friend complex ptor(double mag, double angle=0);
    
57.   complex& topolar(void);
    
58.   complex& torect(void);
    
59. 
    
60.   void operator = (const complex& z) {                // z1 = z2
    
61.     re = z.re;
    
62.     im = z.im;
    
63.   }
    
64.   complex& operator += (const complex& z) {        // z1 += z2
    
65.     re += z.re;
    
66.     im += z.im;
    
67.     return *this;
    
68.   }
    
69.   complex& operator -= (const complex& z) {        // z1 -= z2
    
70.     re -= z.re;
    
71.     im -= z.im;
    
72.     return *this;
    
73.   }
    
74.   complex& operator *= (const complex& z) {        // z1 *= z2
    
75.     *this = *this * z;
    
76.     return *this;
    
77.   }
    
78.   complex& operator /= (const complex& z) {        // z1 /= z2
    
79.     *this = *this / z;
    
80.     return *this;
    
81.   }
    
82.   complex operator + (void) const {                // +z1
    
83.     return *this;
    
84.   }
    
85.   complex operator - (void) const {                // -z1
    
86.     return complex(-re, -im);
    
87.   }
    
88. 
    
89.   friend complex operator + (const complex& z1, const complex& z2) {
    
90.     return complex(z1.re + z2.re, z1.im + z2.im);
    
91.   }
    
92.   friend complex operator + (const complex& z, const double x) {
    
93.     return complex(z.re+x, z.im);
    
94.   }
    
95.   friend complex operator + (const double x, const complex& z) {
    
96.     return complex(z.re+x, z.im);
    
97.   }
    
98.   friend complex operator - (const complex& z1, const complex& z2) {
    
99.     return complex(z1.re - z2.re, z1.im - z2.im);
    
100.   }
     
101.   friend complex operator - (const complex&, const double);
     
102.   friend complex operator - (const double x, const complex& z) {
     
103.     return complex(x-z.re, -z.im);
     
104.   }
     
105.   friend complex operator * (const complex& z1, const complex& z2) {
     
106.     double re = z1.re*z2.re - z1.im*z2.im;
     
107.     double im = z1.re*z2.im + z1.im*z2.re;
     
108.     return complex(re, im);
     
109.   }
     
110.   friend complex operator * (const complex& z, const double x) {
     
111.     return complex(z.re*x, z.im*x);
     
112.   }
     
113.   friend complex operator * (const double x, const complex& z) {
     
114.     return complex(z.re*x, z.im*x);
     
115.   }
     
116.   friend complex operator / (const complex&, const complex&);
     
117.   friend complex operator / (const complex& z, const double x) {
     
118.     return complex(z.re/x, z.im/x);
     
119.   }
     
120.   friend complex operator / (const double, const complex&);
     
121.   friend complex operator ^ (const complex& z1, const complex& z2) {
     
122.     return pow(z1, z2);
     
123.   }
     
124. 
     
125.   friend int operator == (const complex& z1, const complex& z2) {
     
126.     return (z1.re == z2.re) && (z1.im == z2.im);
     
127.   }
     
128.   friend int operator != (const complex& z1, const complex& z2) {
     
129.     return (z1.re != z2.re) || (z1.im != z2.im);
     
130.   }
     
131. 
     
132.   friend double   abs(const complex& z);
     
133.   friend complex sqrt(const complex& z);
     
134.   friend complex pow(const complex& base, const complex& exp);
     
135.   friend complex pow(const complex& base, const double   exp);
     
136.   friend complex pow(const double   base, const complex& exp);
     
137. 
     
138.   friend complex   exp(const complex& z);
     
139.   friend complex   log(const complex& z);
     
140.   friend complex    ln(const complex& z);
     
141.   friend complex log10(const complex& z);
     
142. 
     
143.   friend complex  cos(const complex& z);
     
144.   friend complex  sin(const complex& z);
     
145.   friend complex  tan(const complex& z);
     
146. 
     
147.   friend complex asin(const complex& z);
     
148.   friend complex acos(const complex& z);
     
149.   friend complex atan(const complex& z);
     
150. 
     
151.   friend complex sinh(const complex& z);
     
152.   friend complex cosh(const complex& z);
     
153.   friend complex tanh(const complex& z);
     
154. 
     
155.   void SetArgMode(unsigned char mode) const {
     
156.     if(mode == Z_RADIANS || mode == Z_DEGREES)
     
157.       zArgMode = mode;
     
158.   }
     
159.   void SetPrintMode(unsigned char mode) const {
     
160.     if(mode == Z_COMMA || mode == Z_LETTER)
     
161.       zPrintMode = mode;
     
162.   }
     
163.   void SetLetter(unsigned char letter) const {
     
164.     zLetter = letter;
     
165.   }
     
166. 
     
167.   friend ostream& operator<<(ostream&, const complex&);
     
168.   friend istream& operator>>(istream&, const complex&);
     
169. };
     
170. 
     
171. static const complex Z0(0, 0);                // complex number 0
     
172. static const complex Z1(1, 0);                // complex number 1
     
173. static const complex Zi(0, 1);                // complex number i
     
174. static const complex Zinf(HUGE_VAL, HUGE_VAL); // complex number infinity
     
175. static const complex Complex;
     
176. ///////////////////////////////////////////////
     
177. /////////////////////////////////////////////////
     
178. ///////////////////////////////////////////////
     
179. unsigned char complex::zArgMode   = Z_RADIANS;
     
180. unsigned char complex::zPrintMode = Z_COMMA;
     
181. unsigned char complex::zLetter    = 'i';
     
182. 
     
183. /* -------------------------------------------------------------------- */
     
184. 
     
185. double arg(const complex& z)                // argument (angle)
     
186. {
     
187.   if(z.re == 0 && z.im == 0)                // this is actually a domain error
     
188.     return 0;
     
189.   if(complex::zArgMode == Z_RADIANS)
     
190.     return atan2(z.im, z.re);
     
191.   return atan2(z.im, z.re)/M_PI*180;
     
192. }
     
193. 
     
194. complex ptor(double mag, double angle)  // polar-to-rectangular
     
195. {
     
196.   if(complex::zArgMode == Z_RADIANS)
     
197.     return complex(mag*cos(angle), mag*sin(angle));
     
198.   return complex(mag*cos(angle/180*M_PI), mag*sin(angle/180*M_PI));
     
199. }
     
200. 
     
201. complex rtop(double x, double y)        // rectangular-to-polar
     
202. {
     
203.   if(x == 0 && y == 0)                        // this is actually a domain error
     
204.     return Z0;
     
205.   if(complex::zArgMode == Z_RADIANS)
     
206.     return complex(sqrt(x*x + y*y), atan2(y, x));
     
207.   return complex(sqrt(x*x + y*y), atan2(y, x)*180/M_PI);
     
208. }
     
209. 
     
210. complex& complex::topolar(void)
     
211. {
     
212.   double re_tmp = re;
     
213. 
     
214.   if(re != 0 || im != 0)                // z = (0,0) is a domain error
     
215.   {
     
216.     re = sqrt(re*re + im*im);
     
217.     im = atan2(im, re_tmp);
     
218.   }
     
219. 
     
220.   if(complex::zArgMode == Z_DEGREES)
     
221.     im *= (180/M_PI);
     
222. 
     
223.   return *this;
     
224. }
     
225. 
     
226. complex& complex::torect(void)
     
227. {
     
228.   double re_tmp = re;
     
229. 
     
230.   re = re_tmp*cos(im);
     
231.   im = re_tmp*sin(im);
     
232. 
     
233.   return *this;
     
234. }
     
235. 
     
236. /* ----- Operators ---------------------------------------------------- */
     
237. 
     
238. complex operator/(const complex& z1, const complex& z2)
     
239. {
     
240.   if(z2 == Z0)
     
241.     return complex(Zinf);                // z2 = Z0 is an error!
     
242. 
     
243.   double denom = z2.re*z2.re + z2.im*z2.im;
     
244.   double re = (z1.re*z2.re + z1.im*z2.im)/denom;
     
245.   double im = (z2.re*z1.im - z2.im*z1.re)/denom;
     
246.   return complex(re, im);
     
247. }
     
248. 
     
249. complex operator/(const double x, const complex& z)
     
250. {
     
251.   if(z == Z0)
     
252.     return complex(Zinf);                // z = Z0 is an error!
     
253. 
     
254.   double denom = z.re*z.re + z.im*z.im;
     
255.   return complex(x*z.re/denom, -z.im*x/denom);
     
256. }
     
257. 
     
258. /* ----- Math functions ----------------------------------------------- */
     
259. 
     
260. double abs(const complex& z)
     
261. {
     
262.   return sqrt(z.re*z.re + z.im*z.im);
     
263. }
     
264. 
     
265. complex sqrt(const complex& z)
     
266. {
     
267.   return ptor(sqrt(abs(z)), arg(z)/2);
     
268. }
     
269. 
     
270. complex pow(const complex& base, const double exponent)
     
271. {
     
272.   if(base != Z0 && exponent == 0.0)
     
273.     return complex(1,0);
     
274. 
     
275.   if (base == Z0 && exponent > 0)
     
276.     return Z0;
     
277. 
     
278.   // base == Z0 && exponent == 0 is undefined!
     
279. 
     
280.   return ptor(pow(abs(base), exponent), exponent*arg(base));
     
281. }
     
282. 
     
283. complex pow(const double base, const complex& exponent)
     
284. {
     
285.   if(base != 0.0 && exponent == Z0)
     
286.     return complex(1,0);
     
287. 
     
288.   if (base == 0 && re(exponent) > 0)
     
289.     return complex(0,0);
     
290. 
     
291.   // base == 0 && re(exponent) == 0 is undefined!
     
292. 
     
293.   if(base > 0.0)
     
294.     return exp(exponent * log(fabs(base)));
     
295. 
     
296.   return exp(exponent * complex(log(fabs(base)), M_PI));
     
297. }
     
298. 
     
299. complex pow(const complex& base, const complex& exponent)
     
300. {
     
301.   if(base != Z0 && exponent == Z0)
     
302.     return complex(1,0);
     
303. 
     
304.   if(base == Z0 && re(exponent) > 0)
     
305.     return complex(0,0);
     
306. 
     
307.   // base == Z0 && re(exponent) == 0 is undefined!
     
308. 
     
309.   return exp(exponent * log(base));
     
310. }
     
311. 
     
312. /* ----- Trig functions ----------------------------------------------- */
     
313. 
     
314. complex exp(const complex& z)
     
315. {
     
316.   double mag = exp(z.re);
     
317.   return complex(mag*cos(z.im), mag*sin(z.im));
     
318. }
     
319. 
     
320. complex log(const complex& z)
     
321. {
     
322.   return complex(log(mag(z)), atan2(z.im, z.re));
     
323. }
     
324. /*
     
325. complex log10(const complex& z)
     
326. {
     
327.   return log(z) * M_LOG10E;
     
328. }
     
329. */
     
330. complex sin(const complex& z)
     
331. {
     
332.   return (exp(Zi*z) - exp(-Zi*z))/(2*Zi);
     
333. }
     
334. 
     
335. complex cos(const complex& z)
     
336. {
     
337.   return (exp(Zi*z) + exp(-Zi*z))/2;
     
338. }
     
339. 
     
340. complex tan(const complex& z)
     
341. {
     
342.   return sin(z)/cos(z);
     
343. }
     
344. 
     
345. complex asin(const complex& z)
     
346. {
     
347.   return -Zi*log(Zi*z+sqrt(1-z*z));
     
348. }
     
349. 
     
350. complex acos(const complex& z)
     
351. {
     
352.   return -Zi*log(z+Zi*sqrt(1-z*z));
     
353. }
     
354. 
     
355. complex atan(const complex& z)
     
356. {
     
357.   return -0.5*Zi * log((1+Zi*z)/(1-Zi*z));
     
358. }
     
359. 
     
360. complex sinh(const complex& z)
     
361. {
     
362.   return (exp(z) - exp(-z))/2;
     
363. }
     
364. 
     
365. complex cosh(const complex& z)
     
366. {
     
367.   return (exp(z) + exp(-z))/2;
     
368. }
     
369. 
     
370. complex tanh(const complex& z)
     
371. {
     
372.   return sinh(z)/cosh(z);
     
373. }
     
374. 
     
375. /* ----- Stream I/O --------------------------------------------------- */
     
376. 
     
377. ostream& operator<<(ostream& s, const complex& z)
     
378. {
     
379.   char sign[] = "   ";
     
380. 
     
381.   if(complex::zPrintMode == Z_COMMA)
     
382.     return s << "(" << z.re << ", " << z.im << ")";
     
383. 
     
384.   if(z.im == 0 || z.im/fabs(z.im) == 1)
     
385.     sign[1] = '+';
     
386.   else
     
387.     sign[1] = '-';
     
388.   return s << z.re << sign << complex::zLetter << fabs(z.im);
     
389. 
     
390. }
     
391. 
     
392. istream& operator>>(istream& s, const complex& z)
     
393. {
     
394.   char ch;
     
395. 
     
396.   s >> ch;
     
397.   if(ch == '(')
     
398.   {
     
399.     s >> z.re >> ch;
     
400.     if(ch == ',')
     
401.       s >> z.im >> ch;
     
402.     if(ch != ')')
     
403.       s.clear(ios::badbit | s.rdstate());
     
404.   }
     
405.   else
     
406.   {
     
407.     s.putback(ch);
     
408.     s >> z.re;
     
409.   }
     
410. 
     
411.   return s;
     
412. }
     
413. ///////////////////////////////////////
     
414. ////////////////////fft
     
415. void Fft(complex CA[],int n,int k,complex CB[],int l,int il)
     
416. {
     
417.     int it,m,is,i,j,nv,l0;
     
418.     double p,q,s,vr,vi,poddr,poddi;
     
419.         double *pr,*pi,*fr,*fi;
     
420.         pr=(double*)malloc(n*sizeof(double));
     
421.     pi=(double*)malloc(n*sizeof(double));
     
422.         fr=(double*)malloc(n*sizeof(double));
     
423.         fi=(double*)malloc(n*sizeof(double));
     
424.         for(i=0;i<n;i++)
     
425.         {
     
426.                 pr[i]=CA[i].re;
     
427.                 pi[i]=CA[i].im;
     
428.         }
     
429.     for (it=0; it<=n-1; it++)
     
430.       { m=it; is=0;
     
431.         for (i=0; i<=k-1; i++)
     
432.           { j=m/2; is=2*is+(m-2*j); m=j;}
     
433.         fr[it]=pr[is]; fi[it]=pi[is];
     
434.       }
     
435.     pr[0]=1.0; pi[0]=0.0;
     
436.     p=6.283185306/(1.0*n);
     
437.     pr[1]=cos(p); pi[1]=-sin(p);
     
438.     if (l!=0) pi[1]=-pi[1];
     
439.     for (i=2; i<=n-1; i++)
     
440.       { p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];
     
441.         s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
     
442.         pr[i]=p-q; pi[i]=s-p-q;
     
443.       }
     
444.     for (it=0; it<=n-2; it=it+2)
     
445.       { vr=fr[it]; vi=fi[it];
     
446.         fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1];
     
447.         fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1];
     
448.       }
     
449.     m=n/2; nv=2;
     
450.     for (l0=k-2; l0>=0; l0--)
     
451.       { m=m/2; nv=2*nv;
     
452.         for (it=0; it<=(m-1)*nv; it=it+nv)
     
453.           for (j=0; j<=(nv/2)-1; j++)
     
454.             { p=pr[m*j]*fr[it+j+nv/2];
     
455.               q=pi[m*j]*fi[it+j+nv/2];
     
456.               s=pr[m*j]+pi[m*j];
     
457.               s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
     
458.               poddr=p-q; poddi=s-p-q;
     
459.               fr[it+j+nv/2]=fr[it+j]-poddr;
     
460.               fi[it+j+nv/2]=fi[it+j]-poddi;
     
461.               fr[it+j]=fr[it+j]+poddr;
     
462.               fi[it+j]=fi[it+j]+poddi;
     
463.             }
     
464.       }
     
465.     if (l!=0)
     
466.       for (i=0; i<=n-1; i++)
     
467.         { fr[i]=fr[i]/(1.0*n);
     
468.           fi[i]=fi[i]/(1.0*n);
     
469.         }
     
470.     if (il!=0)
     
471.       for (i=0; i<=n-1; i++)
     
472.         { pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
     
473.           if (fabs(fr[i])<0.000001*fabs(fi[i]))
     
474.             { if ((fi[i]*fr[i])>0) pi[i]=90.0;
     
475.               else pi[i]=-90.0;
     
476.             }
     
477.           else
     
478.             pi[i]=atan(fi[i]/fr[i])*360.0/6.283185306;
     
479.         }
     
480.         for(i=0;i<n;i++)
     
481.         {
     
482.                 CA[i].re=pr[i];
     
483.                 CA[i].im=pi[i];
     
484.         }
     
485.    
     
486.           for(i=0;i<n;i++)
     
487.           {
     
488.                   CB[i].re=fr[i];
     
489.                   CB[i].im=fi[i];
     
490.           }
     
491.     return;
     
492.   }
     
493. /////////////////////////////////////
     
494. void  main()
     
495. {
     
496.     complex a,b;
     
497.         a.re=2;
     
498.         a.im=3;
     
499.         b=a;
     
500.         b=exp(a);
     
501.         cout<<b;
     
502.         int n=64,i,j;
     
503.         double deltaT=0.1,R;
     
504.         ifstream ifile("d:\\work\\tj.fft");
     
505.         ofstream ofile("d:\\work\\tj.out");
     
506.         /*for(i=0;i<n;i++)
     
507.        ofile<<exp(-i*deltaT)<<endl;
     
508.        
     
509.         complex q;
     
510.         double L;
     
511.         int j;
     
512.         for(L=0;L<n;L++)
     
513.         {
     
514.                 ifile.seekg(ios::beg);
     
515.                 q=0;
     
516.                 for(j=0;j<n;j++)
     
517.                 {
     
518.                         ifile>>R;
     
519.                         q+=R*exp(-2*M_PI*L*Zi*L*j/n);
     
520.                 }
     
521.                 q*=deltaT;
     
522.                 ofile<<q<<endl;
     
523.         }
     
524.         CMatrix<int> A;
     
525.         A.Ones(3);
     
526.         cout<<A;
     
527.     CMatrix<complex> B;
     
528.         B.Ones(3);
     
529.         cout<<B;
     
530. */
     
531.         complex CA[64];
     
532.         complex p[64],f[64];
     
533.     for(i=0;i<64;i++)
     
534.         {
     
535.                 ifile>>p[i];
     
536.         }
     
537.         Fft(p,64,6,f,1,1);
     
538.         for(i=0;i<=15;i++)
     
539.         {
     
540.                 for(j=0;j<=3;j++)
     
541.                         ofile<<f[4*i+j]<<"    ";
     
542.                 ofile<<endl;
     
543.         }
     
544.         ofile<<endl<<endl<<endl;
     
545. 
     
546. for(i=0;i<=15;i++)
     
547.         {
     
548.                 for(j=0;j<=3;j++)
     
549.                         ofile<<p[4*i+j]<<"    ";
     
550.                 ofile<<endl;
     
551.         }
     
552. 
     
553. 
     
554. 
     
555.         return ;
     
556. };

复制代码

天堂泪

:lol

assist

有人测试过了吗？

风花雪月

原帖由 assist 于 2007-9-21 09:39 发表

                               
登录/注册后可看大图

有人测试过了吗？

有人测试过，不过会不会因为转贴造成错误，暂时没办法判断