#include <stdio.h>
#include <complex.h>

double _Complex f1(double _Complex x) {
    return (x - 1)*(x + 1)*(x + _Complex_I)*(x - _Complex_I);
}

double _Complex f(double _Complex x) {
    return csin(x) / x;
}

double _Complex fWrap(double _Complex (*f)(double _Complex), double _Complex roots[], int num, double _Complex x) {
    double _Complex res = f(x);
    for (int i = 0; i < num; i++) {
        res /= (x - roots[i]);
    }
    return res;
}

int Muller(double _Complex (*f)(double _Complex), double _Complex x0, double _Complex x1, double _Complex x2, int maxIter, double argPrecision, double funcPrecision, double _Complex *res) {
    double _Complex h1 = x1 - x0;
    double _Complex h2 = x2 - x1;
    double _Complex delta1 = (f(x1) - f(x0))/h1;
    double _Complex delta2 = (f(x2) - f(x1))/h2;
    double _Complex d = (delta2 - delta1)/(h1 + h2);
    double _Complex E, b, D, h, p;
    int i = 2;

    while (i < maxIter) {
        //printf("%f, %f, %f\n", creal(x0), creal(x1), creal(x2));
        //printf("%f, %f, %f\n", cimag(x0), cimag(x1), cimag(x2));
        if (cabs(x2 - x1) < argPrecision || cabs(f(x2)) <= funcPrecision) {
            *res = x2;
            return i;
        }

        b = delta2 + (h2 * d);
        D = csqrt(b * b - 4 * f(x2) * d);

        if (cabs(b - D) < cabs(b + d)) {
            E = b + D;
        } else {
            E = b - D;
        }

        h = -2 * f(x2) / E;
        p = x2 + h;

        x0 = x1;
        x1 = x2;
        x2 = p;

        h1 = x1 - x0;
        h2 = x2 - x1;
        delta1 = (f(x1) - f(x0))/h1;
        delta2 = (f(x2) - f(x1))/h2;
        d = (delta2 - delta1)/(h1 + h2);
        i++;
    }

    return -1;
}

int MullerFunc(double _Complex (*f)(double _Complex), double _Complex x0, double _Complex x1, int rootNum, int maxIter, double argPrecision, double funcPrecision, double _Complex res[]) {
    int founded = 0;
    int lastres = 1;
    while (founded < rootNum && lastres > 0) {
        double _Complex x2 = (x0 + x1) / 2.0;
        double _Complex h1 = x1 - x0;
        double _Complex h2 = x2 - x1;
        double _Complex delta1 = (fWrap(f, res, founded, x1) - fWrap(f, res, founded, x0))/h1;
        double _Complex delta2 = (fWrap(f, res, founded, x2) - fWrap(f, res, founded, x1))/h2;
        double _Complex d = (delta2 - delta1)/(h1 + h2);
        double _Complex E, b, D, h, p;
        int i = 2;

        while (i < maxIter) {
            //printf("%f, %f, %f\n", creal(x0), creal(x1), creal(x2));
            //printf("%f, %f, %f\n", cimag(x0), cimag(x1), cimag(x2));
            if (cabs(x2 - x1) < argPrecision || cabs(fWrap(f, res, founded, x2)) <= funcPrecision) {
                res[founded] = x2;
                founded++;
                lastres = i;
                break;
            }

            b = delta2 + (h2 * d);
            D = csqrt(b * b - 4 * f(x2) * d);

            if (cabs(b - D) < cabs(b + d)) {
                E = b + D;
            } else {
                E = b - D;
            }

            h = -2 * fWrap(f, res, founded, x2) / E;
            p = x2 + h;

            x0 = x1;
            x1 = x2;
            x2 = p;

            h1 = x1 - x0;
            h2 = x2 - x1;
            delta1 = (fWrap(f, res, founded, x1) - fWrap(f, res, founded, x0))/h1;
            delta2 = (fWrap(f, res, founded, x2) - fWrap(f, res, founded, x1))/h2;
            d = (delta2 - delta1)/(h1 + h2);
            i++;
        }

        if (i == maxIter) lastres = -1;
    }

    return founded;
}

int main() {

    double _Complex x0 = -1.0 + 0.0 * _Complex_I;
    double _Complex x1 = 2.0 + 0.0 * _Complex_I;
    double _Complex x2 = 2.0 + 0.0 * _Complex_I;

    double _Complex res[10];

    //printf("%f\n", f(x2));

    //printf("%d\n", Muller(f, x0, x1, x2, 50, 0.005, 0.005, &res));
    int num = MullerFunc(f1, x0, x1, 8, 500, 0.005, 0.005, res);
    for (int i = 0; i < num; i++) {
        printf("real: %f\tim: %f\n", creal(res[i]), cimag(res[i]));
    }
    //printf("%f\n", creal(res));


    return 0;
}