GT Chen's Blog

這是Soccer所寫的函式

注意：
  1:結構定義與自定義(typedef)的變數名字不可一樣
        否則會出現編譯錯誤
  2:只供C使用   C++本身數學函式有內建複數
----------以下開始------------
//--------------the operations of complex------------------------//
struct complex{
                 double Re;
                 double Im;
               };

typedef struct complex _complex;


_complex Cadd(_complex a,_complex b)
{
    _complex c;

    c.Re=a.Re+b.Re;
    c.Im=a.Im+b.Im;

    return c;
}

_complex Csub(_complex a,_complex b)
{
    _complex c;

    c.Re=a.Re-b.Re;
    c.Im=a.Im-b.Im;

    return c;
}


_complex Cmulti(_complex a,_complex b)
{
    _complex c;

    c.Re=a.Re*b.Re-a.Im*b.Im;
    c.Im=a.Im*b.Re+a.Re*b.Im;

    return c;
}


_complex Cdiv(_complex a,_complex b)
{
    _complex c;
    double r,den;

    if(fabs(b.Re)>=fabs(b.Im))
    {
       r=b.Im/b.Re;
       den=b.Re+r*b.Im;
       c.Re=(a.Re+r*a.Im)/den;
       c.Im=(a.Im-r*a.Re)/den;
     }
    else
    {
        r=b.Re/b.Im;
        den=b.Im+r*b.Re;
        c.Re=(a.Re*r+a.Im)/den;
        c.Im=(a.Im*r-a.Re)/den;
     }

     return c;
}


double Ccabs(_complex z)      //the norm of complex
{
     double x,y,ans,temp;
     x=fabs(z.Re);
     y=fabs(z.Im);

     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;
}


_complex Cinner(_complex a,_complex b)   //the inner product of complex
{
   _complex c;

   c.Re=a.Re*b.Re+a.Im*b.Im;
   c.Im=a.Re*b.Im-a.Im*b.Re;

   return c;
}


_complex Cconst(double a)  //transform constant to complex  real part
{
   _complex c;

   c.Re=a;
   c.Im=0.0;

   return c;
}


_complex Cimg(double a)  //transform constant to complex    img part
{
   _complex c;

    c.Re=0.0;
    c.Im=a;

    return c;
}


_complex Csqrt(_complex z)  //the squre of complex
{
    _complex c;
    double x,y,w,r;

    if((z.Re==0.0)&&(z.Im==0.0))
    {
        c.Re=0.0;
        c.Im=0.0;
        return c;
    }
    else
    {
        x=fabs(z.Re);
        y=fabs(z.Im);

        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.Re>=0.0)
        {
           c.Re=w;
           c.Im=z.Im/(2.0*w);
        }
        else
        {
           c.Im=(z.Im>=0.0) ? w: -w;
           c.Re=z.Im/(2.0*c.Im);
        }
    return c;
    }
}


_complex Cconju(_complex z)   //the conjugate of complex
{

    _complex c;

    c.Re=z.Re;
    c.Im=-z.Im;

    return c;
}
-----------以上為止-------------

存成complex.h
只要在程式開頭引入即可使用

EX:

  _complex G1,G2,G3;    //令變數G1,G2,G3為complex variable

  G1.Re=1.0;         //real part of G1=1.0
  G1.Im=4.0;         //imaginary part of G1=4.0

  G2=Cconst(3.0);    //令G2=3+0i
  G3=Cimg(2.0);      //令G3=0+2i

  Cadd(G1,G2);       //G1+G2
  Csub(G1,G2);       //G1-G2
  Cmulti(G1,G2);     //G1*G2
  Cdiv(G1,G2);       //G1/G2
  Ccabs(G1);         //G1的絕對值 (不再是複數)
  Cinner(G1,G2);     //G1,G2的內積
  Csqrt(G1);         //G1的平方
  Cconju(G1);        //G1的共軛數

  printf("G1=%.16e %.16e\n",G1.Re,G1.Im);   //印出G1