durand kerner roots Algorithm

The Durand-Kerner algorithm is an efficient numerical method for finding all the complex roots of a given polynomial simultaneously. Named after its inventors, Robert Durand and Eduard-Wilhelm Kerner, it was introduced in 1960 as a simple yet powerful technique to solve polynomial equations with complex coefficients. The algorithm is based on the principle of simultaneous iterations, in which an initial approximation of the roots is improved iteratively until the desired level of accuracy is achieved. The Durand-Kerner method is particularly advantageous for solving high-degree polynomials, as it converges rapidly and has a low computational complexity compared to other root-finding techniques. The Durand-Kerner algorithm starts by constructing initial approximations of the roots using complex numbers evenly spaced around the unit circle in the complex plane. These approximations are then iteratively refined by dividing the polynomial by the differences between the current approximations and the original root estimates. Convergence is ensured by the fact that the initial approximations lie in the so-called "basins of attraction" of the true roots, meaning that the iterates are drawn closer to the actual roots with every iteration. The algorithm terminates when the absolute differences between successive approximations fall below a predefined threshold, indicating that the estimated roots are accurate enough. The Durand-Kerner method is particularly well-suited for parallel and distributed computing, as the calculations for each root can be performed independently and simultaneously, allowing for substantial speedup in the overall computation time.
#include <math.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <complex.h>

/**
 * Test the algorithm online:
 * https://gist.github.com/kvedala/27f1b0b6502af935f6917673ec43bcd7
 **/

/***
 * Try the highly unstable Wilkinson's polynomial:
 * ./numerical_methods/durand_kerner_roots.c 1 -210 20615 -1256850 53327946 -1672280820 40171771630 -756111184500 11310276995381 -135585182899530 1307535010540395 -10142299865511450 63030812099294896 -311333643161390640 1206647803780373360 -3599979517947607200 8037811822645051776 -12870931245150988800 13803759753640704000 -8752948036761600000 2432902008176640000
 * */

#define ACCURACY 1e-10

/**
 * define polynomial function
 **/
long double complex function(double *coeffs, unsigned int degree, long double complex x)
{
    long double complex out = 0.;
    unsigned int n;

    for (n = 0; n < degree; n++)
        out += coeffs[n] * cpow(x, degree - n - 1);

    return out;
}

static inline char *complex_str(long double complex x)
{
    static char msg[50];
    double r = creal(x);
    double c = cimag(x);

    sprintf(msg, "% 7.04g%+7.04gj", r, c);

    return msg;
}

char check_termination(long double delta)
{
    static long double past_delta = INFINITY;
    if (fabsl(past_delta - delta) <= ACCURACY || delta < ACCURACY)
        return 1;
    past_delta = delta;
    return 0;
}

/***
 * the comandline inputs are taken as coeffiecients of a polynomial
 **/
int main(int argc, char **argv)
{
    double *coeffs = NULL;
    long double complex *s0 = NULL;
    unsigned int degree = 0;
    unsigned int n, i;

    if (argc < 2)
    {
        printf("Please pass the coefficients of the polynomial as commandline arguments.\n");
        return 0;
    }

    degree = argc - 1;                                                              /*< detected polynomial degree */
    coeffs = (double *)malloc(degree * sizeof(double));                             /**< store all input coefficients */
    s0 = (long double complex *)malloc((degree - 1) * sizeof(long double complex)); /**< number of roots = degree-1 */

    /* initialize random seed: */
    srand(time(NULL));

    if (!coeffs || !s0)
    {
        perror("Unable to allocate memory!");
        if (coeffs)
            free(coeffs);
        if (s0)
            free(s0);
        return EXIT_FAILURE;
    }

#if defined(DEBUG) || !defined(NDEBUG)
    /**
     * store intermediate values to a CSV file
     **/
    FILE *log_file = fopen("durand_kerner.log.csv", "wt");
    if (!log_file)
    {
        perror("Unable to create a storage log file!");
        free(coeffs);
        free(s0);
        return EXIT_FAILURE;
    }
    fprintf(log_file, "iter#,");
#endif

    printf("Computing the roots for:\n\t");
    for (n = 0; n < degree; n++)
    {
        coeffs[n] = strtod(argv[n + 1], NULL);
        if (n < degree - 1 && coeffs[n] != 0)
            printf("(%g) x^%d + ", coeffs[n], degree - n - 1);
        else if (coeffs[n] != 0)
            printf("(%g) x^%d = 0\n", coeffs[n], degree - n - 1);

        double tmp;
        if (n > 0)
            coeffs[n] /= tmp; /* numerical errors less when the first coefficient is "1" */
        else
        {
            tmp = coeffs[0];
            coeffs[0] = 1;
        }

        /* initialize root approximations with random values */
        if (n < degree - 1)
        {
            s0[n] = (long double)rand() + (long double)rand() * I;
#if defined(DEBUG) || !defined(NDEBUG)
            fprintf(log_file, "root_%d,", n);
#endif
        }
    }

#if defined(DEBUG) || !defined(NDEBUG)
    fprintf(log_file, "avg. correction");
    fprintf(log_file, "\n0,");
    for (n = 0; n < degree - 1; n++)
        fprintf(log_file, "%s,", complex_str(s0[n]));
#endif

    double tol_condition = 1;
    unsigned long iter = 0;

    while (!check_termination(tol_condition) && iter < INT_MAX)
    {
        long double complex delta = 0;
        tol_condition = 0;
        iter++;

#if defined(DEBUG) || !defined(NDEBUG)
        fprintf(log_file, "\n%ld,", iter);
#endif

        for (n = 0; n < degree - 1; n++)
        {
            long double complex numerator = function(coeffs, degree, s0[n]);
            long double complex denominator = 1.0;
            for (i = 0; i < degree - 1; i++)
                if (i != n)
                    denominator *= s0[n] - s0[i];

            delta = numerator / denominator;

            if (isnan(cabsl(delta)) || isinf(cabsl(delta)))
            {
                printf("\n\nOverflow/underrun error - got value = %Lg", cabsl(delta));
                goto end;
            }

            s0[n] -= delta;

            tol_condition = fmaxl(tol_condition, fabsl(cabsl(delta)));

#if defined(DEBUG) || !defined(NDEBUG)
            fprintf(log_file, "%s,", complex_str(s0[n]));
#endif
        }
        // tol_condition /= (degree - 1);

#if defined(DEBUG) || !defined(NDEBUG)
        if (iter % 500 == 0)
        {
            printf("Iter: %lu\t", iter);
            for (n = 0; n < degree - 1; n++)
                printf("\t%s", complex_str(s0[n]));
            printf("\t\tabsolute average change: %.4g\n", tol_condition);
        }

        fprintf(log_file, "%.4g", tol_condition);
#endif
    }
end:

#if defined(DEBUG) || !defined(NDEBUG)
    fclose(log_file);
#endif

    printf("\nIterations: %lu\n", iter);
    for (n = 0; n < degree - 1; n++)
        printf("\t%s\n", complex_str(s0[n]));
    printf("absolute average change: %.4g\n", tol_condition);

    free(coeffs);
    free(s0);

    return 0;
}

LANGUAGE:

DARK MODE: