Merge pull request #524 from kvedala/numerical-methods

Numerical methods

DIRECTORY.md

@@ -224,6 +224,11 @@
   * [Tower Of Hanoi](https://github.com/TheAlgorithms/C/blob/master/misc/Tower_Of_Hanoi.c)
   * [Union Find](https://github.com/TheAlgorithms/C/blob/master/misc/union_Find.c)
 
+## Numerical Methods
+  * [Durand Kerner Roots](https://github.com/TheAlgorithms/C/blob/master/numerical_methods/durand_kerner_roots.c)
+  * [Qr Decomposition](https://github.com/TheAlgorithms/C/blob/master/numerical_methods/qr_decomposition.c)
+  * [Qr Eigen Values](https://github.com/TheAlgorithms/C/blob/master/numerical_methods/qr_eigen_values.c)
+
 ## Project Euler
   * Problem 01
     * [Sol1](https://github.com/TheAlgorithms/C/blob/master/project_euler/Problem%2001/sol1.c)
@@ -303,4 +308,5 @@
   * [Selection Sort](https://github.com/TheAlgorithms/C/blob/master/sorting/Selection_Sort.c)
   * [Shaker Sort](https://github.com/TheAlgorithms/C/blob/master/sorting/shaker_sort.c)
   * [Shell Sort](https://github.com/TheAlgorithms/C/blob/master/sorting/shell_Sort.c)
+  * [Shell Sort2](https://github.com/TheAlgorithms/C/blob/master/sorting/shell_sort2.c)
   * [Stooge Sort](https://github.com/TheAlgorithms/C/blob/master/sorting/Stooge_Sort.c)

README.md

@@ -4,6 +4,8 @@ C
 
 For a full list of all algorithms, please see: [DIRECTORY.md](https://github.com/TheAlgorithms/C/blob/master/DIRECTORY.md)
 
+All the code can be executed and tested online: [![using Google Colab Notebook](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/gist/kvedala/27f1b0b6502af935f6917673ec43bcd7/plot-durand_kerner-log.ipynb)
+
 ## LeetCode Algorithm
 
 - [Solution](https://github.com/TheAlgorithms/C/tree/master/leetcode) for [LeetCode](https://leetcode.com/problemset/all/)

numerical_methods/durand_kerner_roots.c

@@ -0,0 +1,204 @@
+#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;
+}

numerical_methods/qr_decomposition.c

@@ -0,0 +1,153 @@
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+#define ROWS 4
+#define COLUMNS 3
+
+double A[ROWS][COLUMNS] = {
+    {1, 2, 3},
+    {3, 6, 5},
+    {5, 2, 8},
+    {8, 9, 3}};
+
+void print_matrix(double A[][COLUMNS], int M, int N)
+{
+    for (int row = 0; row < M; row++)
+    {
+        for (int col = 0; col < N; col++)
+            printf("% 9.3g\t", A[row][col]);
+        putchar('\n');
+    }
+    putchar('\n');
+}
+
+void print_2d(double **A, int M, int N)
+{
+    for (int row = 0; row < M; row++)
+    {
+        for (int col = 0; col < N; col++)
+            printf("% 9.3g\t", A[row][col]);
+        putchar('\n');
+    }
+    putchar('\n');
+}
+
+double vector_dot(double *a, double *b, int L)
+{
+    double mag = 0.f;
+    for (int i = 0; i < L; i++)
+        mag += a[i] * b[i];
+
+    return mag;
+}
+
+double vector_mag(double *vector, int L)
+{
+    double dot = vector_dot(vector, vector, L);
+    return sqrt(dot);
+}
+
+double *vector_proj(double *a, double *b, double *out, int L)
+{
+    double num = vector_dot(a, b, L);
+    double deno = vector_dot(b, b, L);
+    for (int i = 0; i < L; i++)
+        out[i] = num * b[i] / deno;
+
+    return out;
+}
+
+double *vector_sub(double *a, double *b, double *out, int L)
+{
+    for (int i = 0; i < L; i++)
+        out[i] = a[i] - b[i];
+
+    return out;
+}
+
+void qr_decompose(double A[][COLUMNS], double **Q, double **R, int M, int N)
+{
+    double *col_vector = (double *)malloc(M * sizeof(double));
+    double *col_vector2 = (double *)malloc(M * sizeof(double));
+    double *tmp_vector = (double *)malloc(M * sizeof(double));
+    for (int i = 0; i < N; i++) /* for each column => R is a square matrix of NxN */
+    {
+        for (int j = 0; j < i; j++) /* second dimension of column */
+            R[i][j] = 0.;           /* make R upper triangular */
+
+        /* get corresponding Q vector */
+        for (int j = 0; j < M; j++)
+        {
+            tmp_vector[j] = A[j][i]; /* accumulator for uk */
+            col_vector[j] = A[j][i];
+        }
+        for (int j = 0; j < i; j++)
+        {
+            for (int k = 0; k < M; k++)
+                col_vector2[k] = Q[k][j];
+            vector_proj(col_vector, col_vector2, col_vector2, M);
+            vector_sub(tmp_vector, col_vector2, tmp_vector, M);
+        }
+        double mag = vector_mag(tmp_vector, M);
+        for (int j = 0; j < M; j++)
+            Q[j][i] = tmp_vector[j] / mag;
+
+        /* compute upper triangular values of R */
+        for (int kk = 0; kk < M; kk++)
+            col_vector[kk] = Q[kk][i];
+        for (int k = i; k < N; k++)
+        {
+            for (int kk = 0; kk < M; kk++)
+                col_vector2[kk] = A[kk][k];
+            R[i][k] = vector_dot(col_vector, col_vector2, M);
+        }
+    }
+
+    free(col_vector);
+    free(col_vector2);
+    free(tmp_vector);
+}
+
+int main(void)
+{
+    // double A[][COLUMNS] = {
+    //     {1, -1, 4},
+    //     {1, 4, -2},
+    //     {1, 4, 2},
+    //     {1, -1, 0}};
+
+    print_matrix(A, ROWS, COLUMNS);
+
+    double **R = (double **)malloc(sizeof(double) * COLUMNS * COLUMNS);
+    double **Q = (double **)malloc(sizeof(double) * ROWS * COLUMNS);
+    if (!Q || !R)
+    {
+        perror("Unable to allocate memory for Q & R!");
+        return -1;
+    }
+    for (int i = 0; i < ROWS; i++)
+    {
+        R[i] = (double *)malloc(sizeof(double) * COLUMNS);
+        Q[i] = (double *)malloc(sizeof(double) * COLUMNS);
+        if (!Q[i] || !R[i])
+        {
+            perror("Unable to allocate memory for Q & R.");
+            return -1;
+        }
+    }
+
+    qr_decompose(A, Q, R, ROWS, COLUMNS);
+
+    print_2d(R, ROWS, COLUMNS);
+    print_2d(Q, ROWS, COLUMNS);
+
+    for (int i = 0; i < ROWS; i++)
+    {
+        free(R[i]);
+        free(Q[i]);
+    }
+    free(R);
+    free(Q);
+    return 0;
+}

numerical_methods/qr_eigen_values.c

@@ -0,0 +1,222 @@
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+
+#define LIMS 9
+
+void create_matrix(double **A, int N)
+{
+    int i, j, tmp, lim2 = LIMS >> 1;
+    srand(time(NULL));
+
+    for (i = 0; i < N; i++)
+    {
+        A[i][i] = (rand() % LIMS) - lim2;
+        for (j = i + 1; j < N; j++)
+        {
+            tmp = (rand() % LIMS) - lim2;
+            A[i][j] = tmp;
+            A[j][i] = tmp;
+        }
+    }
+}
+
+void print_matrix(double **A, int M, int N)
+{
+    for (int row = 0; row < M; row++)
+    {
+        for (int col = 0; col < N; col++)
+            printf("% 9.3g\t", A[row][col]);
+        putchar('\n');
+    }
+    putchar('\n');
+}
+
+double vector_dot(double *a, double *b, int L)
+{
+    double mag = 0.f;
+    for (int i = 0; i < L; i++)
+        mag += a[i] * b[i];
+
+    return mag;
+}
+
+double vector_mag(double *vector, int L)
+{
+    double dot = vector_dot(vector, vector, L);
+    return sqrt(dot);
+}
+
+double *vector_proj(double *a, double *b, double *out, int L)
+{
+    double num = vector_dot(a, b, L);
+    double deno = vector_dot(b, b, L);
+    for (int i = 0; i < L; i++)
+        out[i] = num * b[i] / deno;
+
+    return out;
+}
+
+double *vector_sub(double *a, double *b, double *out, int L)
+{
+    for (int i = 0; i < L; i++)
+        out[i] = a[i] - b[i];
+
+    return out;
+}
+
+void qr_decompose(double **A, double **Q, double **R, int M, int N)
+{
+    double *col_vector = (double *)malloc(M * sizeof(double));
+    double *col_vector2 = (double *)malloc(M * sizeof(double));
+    double *tmp_vector = (double *)malloc(M * sizeof(double));
+    for (int i = 0; i < N; i++) /* for each column => R is a square matrix of NxN */
+    {
+        for (int j = 0; j < i; j++) /* second dimension of column */
+            R[i][j] = 0.;           /* make R upper triangular */
+
+        /* get corresponding Q vector */
+        for (int j = 0; j < M; j++)
+        {
+            tmp_vector[j] = A[j][i]; /* accumulator for uk */
+            col_vector[j] = A[j][i];
+        }
+        for (int j = 0; j < i; j++)
+        {
+            for (int k = 0; k < M; k++)
+                col_vector2[k] = Q[k][j];
+            vector_proj(col_vector, col_vector2, col_vector2, M);
+            vector_sub(tmp_vector, col_vector2, tmp_vector, M);
+        }
+        double mag = vector_mag(tmp_vector, M);
+        for (int j = 0; j < M; j++)
+            Q[j][i] = tmp_vector[j] / mag;
+
+        /* compute upper triangular values of R */
+        for (int kk = 0; kk < M; kk++)
+            col_vector[kk] = Q[kk][i];
+        for (int k = i; k < N; k++)
+        {
+            for (int kk = 0; kk < M; kk++)
+                col_vector2[kk] = A[kk][k];
+            R[i][k] = vector_dot(col_vector, col_vector2, M);
+        }
+    }
+
+    free(col_vector);
+    free(col_vector2);
+    free(tmp_vector);
+}
+
+double **mat_mul(double **A, double **B, double **OUT, int R1, int C1, int R2, int C2)
+{
+    if (C1 != R2)
+    {
+        perror("Matrix dimensions mismatch!");
+        return OUT;
+    }
+    for (int i = 0; i < R1; i++)
+        for (int j = 0; j < C2; j++)
+        {
+            OUT[i][j] = 0.f;
+            for (int k = 0; k < C1; k++)
+                OUT[i][j] += A[i][k] * B[k][j];
+        }
+    return OUT;
+}
+
+int main(int argc, char **argv)
+{
+    int mat_size = 5;
+    if (argc == 2)
+        mat_size = atoi(argv[1]);
+
+    if (mat_size < 2)
+    {
+        fprintf(stderr, "Matrix size should be > 2\n");
+        return -1;
+    }
+
+    int i, rows = mat_size, columns = mat_size;
+
+    double **A = (double **)malloc(sizeof(double) * mat_size);
+    double **R = (double **)malloc(sizeof(double) * mat_size);
+    double **Q = (double **)malloc(sizeof(double) * mat_size);
+    double *eigen_vals = (double *)malloc(sizeof(double) * mat_size);
+    if (!Q || !R || !eigen_vals)
+    {
+        perror("Unable to allocate memory for Q & R!");
+        return -1;
+    }
+    for (i = 0; i < mat_size; i++)
+    {
+        A[i] = (double *)malloc(sizeof(double) * mat_size);
+        R[i] = (double *)malloc(sizeof(double) * mat_size);
+        Q[i] = (double *)malloc(sizeof(double) * mat_size);
+        if (!Q[i] || !R[i])
+        {
+            perror("Unable to allocate memory for Q & R.");
+            return -1;
+        }
+    }
+
+    create_matrix(A, mat_size);
+
+    print_matrix(A, mat_size, mat_size);
+
+    int counter = 0, num_eigs = rows - 1;
+    double last_eig = 0;
+
+    while (num_eigs > 0)
+    {
+        while (fabs(A[num_eigs][num_eigs - 1]) > 1e-10)
+        {
+            last_eig = A[num_eigs][num_eigs];
+            for (int i = 0; i < rows; i++)
+                A[i][i] -= last_eig; /* A - cI */
+            qr_decompose(A, Q, R, rows, columns);
+
+#if defined(DEBUG) || !defined(NDEBUG)
+            print_matrix(A, rows, columns);
+            print_matrix(Q, rows, columns);
+            print_matrix(R, columns, columns);
+            printf("-------------------- %d ---------------------\n", ++counter);
+#endif
+            mat_mul(R, Q, A, columns, columns, rows, columns);
+            for (int i = 0; i < rows; i++)
+                A[i][i] += last_eig; /* A + cI */
+        }
+
+        eigen_vals[num_eigs] = A[num_eigs][num_eigs];
+#if defined(DEBUG) || !defined(NDEBUG)
+        printf("========================\n");
+        printf("Eigen value: % g,\n", last_eig);
+        printf("========================\n");
+#endif
+        num_eigs--;
+        rows--;
+        columns--;
+    }
+    eigen_vals[0] = A[0][0];
+
+#if defined(DEBUG) || !defined(NDEBUG)
+    print_matrix(R, mat_size, mat_size);
+    print_matrix(Q, mat_size, mat_size);
+#endif
+    printf("Eigen vals: ");
+    for (i = 0; i < mat_size; i++)
+        printf("% 9.4g\t", eigen_vals[i]);
+
+    for (int i = 0; i < mat_size; i++)
+    {
+        free(A[i]);
+        free(R[i]);
+        free(Q[i]);
+    }
+    free(A);
+    free(R);
+    free(Q);
+    free(eigen_vals);
+    return 0;
+}

sorting/shell_sort2.c

@@ -0,0 +1,78 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+#define ELEMENT_NR 20000
+#define ARRAY_LEN(x) (sizeof(x) / sizeof((x)[0]))
+
+void show_data(int arr[], int len)
+{
+    int i;
+
+    for (i = 0; i < len; i++)
+        printf("%3d ", arr[i]);
+    printf("\n");
+}
+
+void swap(int *a, int *b)
+{
+    int tmp;
+
+    tmp = *a;
+    *a = *b;
+    *b = tmp;
+}
+
+/**
+ * Optimized algorithm - takes half the time as other
+ **/
+void shell_sort(int array[], int LEN)
+{
+    const int gaps[] = {701, 301, 132, 57, 23, 10, 4, 1};
+    const int gap_len = 8;
+    int i, j, g;
+
+    for (g = 0; g < gap_len; g++)
+    {
+        int gap = gaps[g];
+        for (i = gap; i < LEN; i++)
+        {
+            int tmp = array[i];
+
+            for (j = i; j >= gap && (array[j - gap] - tmp) > 0; j -= gap)
+                array[j] = array[j - gap];
+            array[j] = tmp;
+        }
+    }
+#ifdef DEBUG
+    for (i = 0; i < LEN; i++)
+        printf("%s\t", data[i]);
+#endif
+}
+
+int main(int argc, char *argv[])
+{
+    int i;
+    int array[ELEMENT_NR];
+    int range = 500;
+    int size;
+    double time_spent;
+
+    srand(time(NULL));
+    for (i = 0; i < ELEMENT_NR; i++)
+        array[i] = rand() % range + 1;
+
+    size = ARRAY_LEN(array);
+
+    show_data(array, size);
+    clock_t t1 = clock();
+    shell_sort(array, size);
+    clock_t t2 = clock();
+
+    printf("Data Sorted\n");
+    show_data(array, size);
+
+    printf("Time spent sorting: %.4g s\n", (t2 - t1) / CLOCKS_PER_SEC);
+
+    return 0;
+}