00001
00002
00003 #include <math.h>
00004
00005 typedef struct FCOMPLEX {float r,i;} fcomplex;
00006
00007 fcomplex Cadd(a,b)
00008 fcomplex a,b;
00009 { fcomplex c;
00010 c.r=a.r+b.r;
00011 c.i=a.i+b.i;
00012 return c;
00013 }
00014
00015 fcomplex Csub(a,b)
00016 fcomplex a,b;
00017 { fcomplex c;
00018 c.r=a.r-b.r;
00019 c.i=a.i-b.i;
00020 return c;
00021 }
00022
00023 fcomplex Cmul(a,b)
00024 fcomplex a,b;
00025 { fcomplex c;
00026 c.r=a.r*b.r-a.i*b.i;
00027 c.i=a.i*b.r+a.r*b.i;
00028 return c;
00029 }
00030
00031 fcomplex Complex(re,im)
00032 float re,im;
00033 { fcomplex c;
00034 c.r=re;
00035 c.i=im;
00036 return c;
00037 }
00038
00039 fcomplex Conjg(z)
00040 fcomplex z;
00041 { fcomplex c;
00042 c.r=z.r;
00043 c.i = -z.i;
00044 return c;
00045 }
00046
00047 fcomplex Cdiv(a,b)
00048 fcomplex a,b;
00049 { fcomplex c;
00050 float r,den;
00051 if (fabs(b.r) >= fabs(b.i)) {
00052 r=b.i/b.r;
00053 den=b.r+r*b.i;
00054 c.r=(a.r+r*a.i)/den;
00055 c.i=(a.i-r*a.r)/den;
00056 } else {
00057 r=b.r/b.i;
00058 den=b.i+r*b.r;
00059 c.r=(a.r*r+a.i)/den;
00060 c.i=(a.i*r-a.r)/den;
00061 }
00062 return c;
00063 }
00064
00065 float Cabs(z)
00066 fcomplex z;
00067 { float x,y,ans,temp;
00068 x=fabs(z.r);
00069 y=fabs(z.i);
00070 if (x == 0.0)
00071 ans=y;
00072 else if (y == 0.0)
00073 ans=x;
00074 else if (x > y) {
00075 temp=y/x;
00076 ans=x*sqrt(1.0+temp*temp);
00077 } else {
00078 temp=x/y;
00079 ans=y*sqrt(1.0+temp*temp);
00080 }
00081 return ans;
00082 }
00083
00084 fcomplex Csqrt(z)
00085 fcomplex z;
00086 { fcomplex c;
00087 float x,y,w,r;
00088 if ((z.r == 0.0) && (z.i == 0.0)) {
00089 c.r=0.0;
00090 c.i=0.0;
00091 return c;
00092 } else {
00093 x=fabs(z.r);
00094 y=fabs(z.i);
00095 if (x >= y) {
00096 r=y/x;
00097 w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
00098 } else {
00099 r=x/y;
00100 w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
00101 }
00102 if (z.r >= 0.0) {
00103 c.r=w;
00104 c.i=z.i/(2.0*w);
00105 } else {
00106 c.i=(z.i >= 0) ? w : -w;
00107 c.r=z.i/(2.0*c.i);
00108 }
00109 return c;
00110 }
00111 }
00112
00113 fcomplex RCmul(x,a)
00114 float x;
00115 fcomplex a;
00116 { fcomplex c;
00117 c.r=x*a.r;
00118 c.i=x*a.i;
00119 return c;
00120 }
