251 lines
3.2 KiB
C
251 lines
3.2 KiB
C
#include <math.h>
|
|
|
|
typedef struct FCOMPLEX {float r,i;} fcomplex;
|
|
|
|
#if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */
|
|
|
|
fcomplex Cadd(fcomplex a, fcomplex b)
|
|
{
|
|
fcomplex c;
|
|
c.r=a.r+b.r;
|
|
c.i=a.i+b.i;
|
|
return c;
|
|
}
|
|
|
|
fcomplex Csub(fcomplex a, fcomplex b)
|
|
{
|
|
fcomplex c;
|
|
c.r=a.r-b.r;
|
|
c.i=a.i-b.i;
|
|
return c;
|
|
}
|
|
|
|
|
|
fcomplex Cmul(fcomplex a, fcomplex b)
|
|
{
|
|
fcomplex c;
|
|
c.r=a.r*b.r-a.i*b.i;
|
|
c.i=a.i*b.r+a.r*b.i;
|
|
return c;
|
|
}
|
|
|
|
fcomplex Complex(float re, float im)
|
|
{
|
|
fcomplex c;
|
|
c.r=re;
|
|
c.i=im;
|
|
return c;
|
|
}
|
|
|
|
fcomplex Conjg(fcomplex z)
|
|
{
|
|
fcomplex c;
|
|
c.r=z.r;
|
|
c.i = -z.i;
|
|
return c;
|
|
}
|
|
|
|
fcomplex Cdiv(fcomplex a, fcomplex b)
|
|
{
|
|
fcomplex c;
|
|
float r,den;
|
|
if (fabs(b.r) >= fabs(b.i)) {
|
|
r=b.i/b.r;
|
|
den=b.r+r*b.i;
|
|
c.r=(a.r+r*a.i)/den;
|
|
c.i=(a.i-r*a.r)/den;
|
|
} else {
|
|
r=b.r/b.i;
|
|
den=b.i+r*b.r;
|
|
c.r=(a.r*r+a.i)/den;
|
|
c.i=(a.i*r-a.r)/den;
|
|
}
|
|
return c;
|
|
}
|
|
|
|
float Cabs(fcomplex z)
|
|
{
|
|
float x,y,ans,temp;
|
|
x=fabs(z.r);
|
|
y=fabs(z.i);
|
|
if (x == 0.0)
|
|
ans=y;
|
|
else if (y == 0.0)
|
|
ans=x;
|
|
else if (x > y) {
|
|
temp=y/x;
|
|
ans=x*sqrt(1.0+temp*temp);
|
|
} else {
|
|
temp=x/y;
|
|
ans=y*sqrt(1.0+temp*temp);
|
|
}
|
|
return ans;
|
|
}
|
|
|
|
fcomplex Csqrt(fcomplex z)
|
|
{
|
|
fcomplex c;
|
|
float x,y,w,r;
|
|
if ((z.r == 0.0) && (z.i == 0.0)) {
|
|
c.r=0.0;
|
|
c.i=0.0;
|
|
return c;
|
|
} else {
|
|
x=fabs(z.r);
|
|
y=fabs(z.i);
|
|
if (x >= y) {
|
|
r=y/x;
|
|
w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
|
|
} else {
|
|
r=x/y;
|
|
w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
|
|
}
|
|
if (z.r >= 0.0) {
|
|
c.r=w;
|
|
c.i=z.i/(2.0*w);
|
|
} else {
|
|
c.i=(z.i >= 0) ? w : -w;
|
|
c.r=z.i/(2.0*c.i);
|
|
}
|
|
return c;
|
|
}
|
|
}
|
|
|
|
fcomplex RCmul(float x, fcomplex a)
|
|
{
|
|
fcomplex c;
|
|
c.r=x*a.r;
|
|
c.i=x*a.i;
|
|
return c;
|
|
}
|
|
|
|
#else /* ANSI */
|
|
/* traditional - K&R */
|
|
|
|
fcomplex Cadd(a,b)
|
|
fcomplex a,b;
|
|
{
|
|
fcomplex c;
|
|
c.r=a.r+b.r;
|
|
c.i=a.i+b.i;
|
|
return c;
|
|
}
|
|
|
|
fcomplex Csub(a,b)
|
|
fcomplex a,b;
|
|
{
|
|
fcomplex c;
|
|
c.r=a.r-b.r;
|
|
c.i=a.i-b.i;
|
|
return c;
|
|
}
|
|
|
|
|
|
fcomplex Cmul(a,b)
|
|
fcomplex a,b;
|
|
{
|
|
fcomplex c;
|
|
c.r=a.r*b.r-a.i*b.i;
|
|
c.i=a.i*b.r+a.r*b.i;
|
|
return c;
|
|
}
|
|
|
|
fcomplex Complex(re,im)
|
|
float im,re;
|
|
{
|
|
fcomplex c;
|
|
c.r=re;
|
|
c.i=im;
|
|
return c;
|
|
}
|
|
|
|
fcomplex Conjg(z)
|
|
fcomplex z;
|
|
{
|
|
fcomplex c;
|
|
c.r=z.r;
|
|
c.i = -z.i;
|
|
return c;
|
|
}
|
|
|
|
fcomplex Cdiv(a,b)
|
|
fcomplex a,b;
|
|
{
|
|
fcomplex c;
|
|
float r,den;
|
|
if (fabs(b.r) >= fabs(b.i)) {
|
|
r=b.i/b.r;
|
|
den=b.r+r*b.i;
|
|
c.r=(a.r+r*a.i)/den;
|
|
c.i=(a.i-r*a.r)/den;
|
|
} else {
|
|
r=b.r/b.i;
|
|
den=b.i+r*b.r;
|
|
c.r=(a.r*r+a.i)/den;
|
|
c.i=(a.i*r-a.r)/den;
|
|
}
|
|
return c;
|
|
}
|
|
|
|
float Cabs(z)
|
|
fcomplex z;
|
|
{
|
|
float x,y,ans,temp;
|
|
x=fabs(z.r);
|
|
y=fabs(z.i);
|
|
if (x == 0.0)
|
|
ans=y;
|
|
else if (y == 0.0)
|
|
ans=x;
|
|
else if (x > y) {
|
|
temp=y/x;
|
|
ans=x*sqrt(1.0+temp*temp);
|
|
} else {
|
|
temp=x/y;
|
|
ans=y*sqrt(1.0+temp*temp);
|
|
}
|
|
return ans;
|
|
}
|
|
|
|
fcomplex Csqrt(z)
|
|
fcomplex z;
|
|
{
|
|
fcomplex c;
|
|
float x,y,w,r;
|
|
if ((z.r == 0.0) && (z.i == 0.0)) {
|
|
c.r=0.0;
|
|
c.i=0.0;
|
|
return c;
|
|
} else {
|
|
x=fabs(z.r);
|
|
y=fabs(z.i);
|
|
if (x >= y) {
|
|
r=y/x;
|
|
w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
|
|
} else {
|
|
r=x/y;
|
|
w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
|
|
}
|
|
if (z.r >= 0.0) {
|
|
c.r=w;
|
|
c.i=z.i/(2.0*w);
|
|
} else {
|
|
c.i=(z.i >= 0) ? w : -w;
|
|
c.r=z.i/(2.0*c.i);
|
|
}
|
|
return c;
|
|
}
|
|
}
|
|
|
|
fcomplex RCmul(x,a)
|
|
fcomplex a;
|
|
float x;
|
|
{
|
|
fcomplex c;
|
|
c.r=x*a.r;
|
|
c.i=x*a.i;
|
|
return c;
|
|
}
|
|
|
|
#endif /* ANSI */
|