理解FFT之一 復數

最近在看FFT的內容。參考了2個資料，一個是CLRS，一個是一本不知名的算法書。

都是以多項式為切入口接入FFT。

思路大致是，多項式有2中表達方式：系數，以及值對。

值對的表達方式非常有利於多項式之間的相互運算，同時也需要二者之間的相互轉化。

其中一個簡化計算的思路就是巧妙的選取x的值，減少計算量。

這就涉及到復數的使用。遂開始寫一個復數的類。用來實現和驗證fft算法。

算是自己第一個動手造的輪子。

代碼比較丑陋，以后再慢慢改吧。

#ifndef MYCOMPLEX_H
#define MYCOMPLEX_H

#include <math.h>
#include <iostream>
using namespace std;

namespace MyComplex
{

#define PI 3.1415926

    struct mycomplex
    {
        double r, theta;
        mycomplex(const double& rin =0, const double& thetain=0) : r(rin), theta(thetain){};
        double re() const;
        double im() const;
        void fromReIm(double re, double im);
        bool operator==(const mycomplex& c) const;
        const mycomplex& operator=(const mycomplex& c);
        const mycomplex operator+(const mycomplex& c) const;
        const mycomplex operator-(const mycomplex& c) const;
        const mycomplex operator*(const mycomplex& c) const;
        const mycomplex operator/(const mycomplex& c) const;
        const mycomplex& operator+=(const mycomplex& c);
        const mycomplex& operator-=(const mycomplex& c);
        const mycomplex& operator*=(const mycomplex& c);
        const mycomplex& operator/=(const mycomplex& c);
        const mycomplex operator-() const;
        const mycomplex pow(double p) const;
        const mycomplex conj()const;

        ostream& print_p(ostream& os);
    private:
        void validate();
        //ostream& operator<<(ostream &os);
    };

    double mycomplex::re() const
    {
        return r*cos(theta);
    }

    double mycomplex::im() const
    {
        return r*sin(theta);
    }

    void mycomplex::fromReIm(double re, double im)
    {
        r = sqrt(re*re + im*im);
        if (re == 0.0)
        {
            if (im ==0.0)
            {
                theta = 0;
            }
            else
            {
                theta = (im > 0) ? (PI / 2) : (PI / 2 * 3);
            }
        }
        else
        {
            if (re<0)
            {
                theta = PI + atan(im / re);
            }
            if (re>0 && im<=0)
            {
                theta = 2 * PI + atan(im / re);
            }
            if (re>0 && im>=0)
            {
                theta = atan(im / re);
            }
            
        }
    }

    bool mycomplex::operator==(const mycomplex& c) const
    {
        return r == c.r && theta == c.theta;
    }

    const mycomplex mycomplex::operator+(const mycomplex& c) const
    {
        mycomplex ret;
        ret.fromReIm(r*cos(theta) + c.r*cos(c.theta), r*sin(theta) + c.r*sin(c.theta));
        return ret;
    }

    const mycomplex mycomplex::operator*(const mycomplex& c) const
    {
        mycomplex ret;
        ret.r = r*c.r;
        ret.theta = theta + c.theta;
        ret.validate();
        return ret;
    }

    const mycomplex mycomplex::operator/(const mycomplex& c) const
    {
        mycomplex ret;
        ret.r = r / c.r;
        ret.theta = theta - c.theta;
        ret.validate();
        return ret;
    }

    const mycomplex& mycomplex::operator*=(const mycomplex& c)
    {
        *this = *this * c;
        return *this;
    }

    const mycomplex& mycomplex::operator/=(const mycomplex& c)
    {
        *this = *this / c;
        return *this;
    }

    const mycomplex& mycomplex::operator-=(const mycomplex& c)
    {
        *this = *this - c;
        return *this;
    }

    const mycomplex mycomplex::pow(double p) const
    {
        mycomplex ret;
        ret.r = ::pow(r, p);
        ret.theta = theta * p;
        ret.validate();
        return ret;
    }

    const mycomplex mycomplex::conj() const
    {
        mycomplex ret;
        ret.r = r;
        ret.theta = 2 * PI - theta;
        return ret;
    }

    ostream& mycomplex::print_p(ostream& os)
    {
        os << r << "*exp(" << theta << "i" << ")";
        return os;
    }

    void mycomplex::validate()
    {
        if (theta > 0)
        {
            while (theta > 2 *PI)
            {
                theta -= 2 * PI;
            }
        }
        else
        {
            while (theta / PI < 0)
            {
                theta += 2 * PI;
            }
        }
    }

    const mycomplex& mycomplex::operator+=(const mycomplex& c)
    {
        *this = *this + c;
        return *this;
    }

    const mycomplex mycomplex::operator-() const
    {
        mycomplex ret(r, ((theta <= PI) ? (theta + PI) : (theta - PI)));
        return ret;
    }

    const mycomplex mycomplex::operator-(const mycomplex& c) const
    {
        mycomplex ret;
        ret = *this + (-c);
        return ret;
    }

    const mycomplex& mycomplex::operator=(const mycomplex& c)
    {
        r=c.r;
        theta = c.theta;
        return *this;
    }

    ostream& operator<<(ostream &os, const mycomplex& c)
    {
        os << c.re();
        if (c.im()<0)
        {
            os << c.im() << "i";
        }
        if (c.im() > 0)
        {
            os << "+" << c.im() << "i";
        }
        return os;
    }

}    //namespace MyComplex

#endif