diff --git a/DIRECTORY.md b/DIRECTORY.md index 5066685c..d2d91054 100644 --- a/DIRECTORY.md +++ b/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) diff --git a/README.md b/README.md index 21e09760..c0a4398b 100644 --- a/README.md +++ b/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/) diff --git a/numerical_methods/durand_kerner_roots.c b/numerical_methods/durand_kerner_roots.c new file mode 100644 index 00000000..c88e7c8c --- /dev/null +++ b/numerical_methods/durand_kerner_roots.c @@ -0,0 +1,204 @@ +#include +#include +#include +#include +#include +#include +#include + +/** + * 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; +} \ No newline at end of file diff --git a/numerical_methods/qr_decomposition.c b/numerical_methods/qr_decomposition.c new file mode 100644 index 00000000..a681bbf3 --- /dev/null +++ b/numerical_methods/qr_decomposition.c @@ -0,0 +1,153 @@ +#include +#include +#include + +#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; +} \ No newline at end of file diff --git a/numerical_methods/qr_eigen_values.c b/numerical_methods/qr_eigen_values.c new file mode 100644 index 00000000..42183751 --- /dev/null +++ b/numerical_methods/qr_eigen_values.c @@ -0,0 +1,222 @@ +#include +#include +#include +#include + +#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; +} \ No newline at end of file diff --git a/sorting/shell_sort2.c b/sorting/shell_sort2.c new file mode 100644 index 00000000..7802edf1 --- /dev/null +++ b/sorting/shell_sort2.c @@ -0,0 +1,78 @@ +#include +#include +#include + +#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; +}