Метод Ньютона для комплексных чисел [закрыт]

Здраствуйте, должен решить уравнение z^3 - 1 = 0, z - комплексное число.

1. Код работает, только вывод один корень из трех, но верный. Что нудно сделать, чтобы выводило три корня? Как мне кажется, проблема в том, что нигде не задается (не передается) стартовая точка и по этому один корень, а так бы для разных точек были бы разные корни, но я в этом не уверен.
2. На основе результатов нужно будет сделать график бассейнов решений, т.е от это будет комплексная плоскость, на которой каждый пиксель будет приводить одному из решений. Как это сделать, может кто-то знает? Как мне кажется, это должно быть два цикла с очень маленьким шагом, но что дальше?

Код :

#include <cmath>
#include <iostream>
using namespace std;
class CComplex
{
public :
    double _re;
    double _im;
    CComplex();
    CComplex(double re, double im);
    CComplex(const CComplex &pCopy);
    CComplex conjugate();
    double complexfabs();
    CComplex operator + (CComplex pAdd);
    CComplex operator - (CComplex pAdd);
    CComplex operator * (CComplex pMul);
    CComplex operator / (CComplex pDiv);
    CComplex uintpower(unsigned int n);
    friend ostream& operator << (ostream &out, const CComplex pVal)
    {
        out<<pVal._re<<(pVal._im < 0 ? " - " : " + ")<<fabs(pVal._im)<<"i";
        return out;
    }
};
CComplex::CComplex()
{
    _re = _im = 0;
}
CComplex::CComplex(double re, double im)
{
    _re = re;
    _im = im;
}
CComplex::CComplex(const CComplex &pCopy)
{
    _re = pCopy._re;
    _im = pCopy._im;
}
CComplex
CComplex::conjugate()
{
    CComplex pThis = (*this);
    pThis._im *= -1;
    return pThis;
}
double CComplex::complexfabs()
{
    return _re*_re + _im*_im;
}
CComplex
CComplex::operator +(CComplex pAdd)
{
    CComplex pThis = (*this);
    pThis._re += pAdd._re;
    pThis._im += pAdd._im;
    return pThis;
}
CComplex
CComplex::operator - (CComplex pAdd)
{
    CComplex pThis = (*this);
    pThis._re -= pAdd._re;
    pThis._im -= pAdd._im;
    return pThis;
}
CComplex
CComplex::operator *(CComplex pMul)
{
    CComplex pThis = (*this);
    pThis._re = _re*pMul._re - _im*pMul._im;
    pThis._im = _re*pMul._im + _im*pMul._re;
    return pThis;
}
CComplex
CComplex::operator / (CComplex pDiv)
{
    CComplex pThis = (*this);
    double divider = pDiv.complexfabs();
    pThis = pThis*pDiv.conjugate();
    pThis._re /= divider;
    pThis._im /= divider;
    return pThis;
}
CComplex
CComplex::uintpower(unsigned int n)
{
    CComplex pThis = CComplex(1, 0);
    for(unsigned int i = 0; i < n; i++)
        pThis = pThis*(*this);
    return pThis;
}
CComplex f(CComplex z, unsigned int n)
{
    return z.uintpower(n) - CComplex(1, 0);
}
CComplex df(CComplex z, unsigned int n)
{
    return CComplex(n, 0)*z.uintpower(n - 1);
}
int main()
{
    unsigned int n = 3;
    unsigned int i;
    double err = 1e-12;
    CComplex z(1,1);

    std::cout << "First root" <<std::endl;
    for(i = 0; f(z, n).complexfabs() > err; i++)
    {
        z = z - f(z, n) / df(z, n);
        cout<<"\riteration "<<i + 1<<" root : "<<z<<std::endl;    
    }
   system("pause");
    return 0;
}