Skip to content

Commit 0538404

Browse files
StepfenShawnStepfenShawn
authored
Merge pull request TheAlgorithms#524 from kvedala/numerical-methods
Numerical methods
2 parents d1f6643 + dbcd87a commit 0538404

File tree

6 files changed

+665
-0
lines changed

6 files changed

+665
-0
lines changed

‎DIRECTORY.md

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -224,6 +224,11 @@
224224
* [Tower Of Hanoi](https://github.com/TheAlgorithms/C/blob/master/misc/Tower_Of_Hanoi.c)
225225
* [Union Find](https://github.com/TheAlgorithms/C/blob/master/misc/union_Find.c)
226226

227+
## Numerical Methods
228+
* [Durand Kerner Roots](https://github.com/TheAlgorithms/C/blob/master/numerical_methods/durand_kerner_roots.c)
229+
* [Qr Decomposition](https://github.com/TheAlgorithms/C/blob/master/numerical_methods/qr_decomposition.c)
230+
* [Qr Eigen Values](https://github.com/TheAlgorithms/C/blob/master/numerical_methods/qr_eigen_values.c)
231+
227232
## Project Euler
228233
* Problem 01
229234
* [Sol1](https://github.com/TheAlgorithms/C/blob/master/project_euler/Problem%2001/sol1.c)
@@ -303,4 +308,5 @@
303308
* [Selection Sort](https://github.com/TheAlgorithms/C/blob/master/sorting/Selection_Sort.c)
304309
* [Shaker Sort](https://github.com/TheAlgorithms/C/blob/master/sorting/shaker_sort.c)
305310
* [Shell Sort](https://github.com/TheAlgorithms/C/blob/master/sorting/shell_Sort.c)
311+
* [Shell Sort2](https://github.com/TheAlgorithms/C/blob/master/sorting/shell_sort2.c)
306312
* [Stooge Sort](https://github.com/TheAlgorithms/C/blob/master/sorting/Stooge_Sort.c)

‎README.md

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,8 @@ C
44

55
For a full list of all algorithms, please see: [DIRECTORY.md](https://github.com/TheAlgorithms/C/blob/master/DIRECTORY.md)
66

7+
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)
8+
79
## LeetCode Algorithm
810

911
- [Solution](https://github.com/TheAlgorithms/C/tree/master/leetcode) for [LeetCode](https://leetcode.com/problemset/all/)

‎numerical_methods/durand_kerner_roots.c

Lines changed: 204 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,204 @@
1+
#include <math.h>
2+
#include <time.h>
3+
#include <stdio.h>
4+
#include <stdlib.h>
5+
#include <string.h>
6+
#include <limits.h>
7+
#include <complex.h>
8+
9+
/**
10+
* Test the algorithm online:
11+
* https://gist.github.com/kvedala/27f1b0b6502af935f6917673ec43bcd7
12+
**/
13+
14+
/***
15+
* Try the highly unstable Wilkinson's polynomial:
16+
* ./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
17+
* */
18+
19+
#define ACCURACY 1e-10
20+
21+
/**
22+
* define polynomial function
23+
**/
24+
long double complex function(double *coeffs, unsigned int degree, long double complex x)
25+
{
26+
long double complex out = 0.;
27+
unsigned int n;
28+
29+
for (n = 0; n < degree; n++)
30+
out += coeffs[n] * cpow(x, degree - n - 1);
31+
32+
return out;
33+
}
34+
35+
static inline char *complex_str(long double complex x)
36+
{
37+
static char msg[50];
38+
double r = creal(x);
39+
double c = cimag(x);
40+
41+
sprintf(msg, "% 7.04g%+7.04gj", r, c);
42+
43+
return msg;
44+
}
45+
46+
char check_termination(long double delta)
47+
{
48+
static long double past_delta = INFINITY;
49+
if (fabsl(past_delta - delta) <= ACCURACY || delta < ACCURACY)
50+
return 1;
51+
past_delta = delta;
52+
return 0;
53+
}
54+
55+
/***
56+
* the comandline inputs are taken as coeffiecients of a polynomial
57+
**/
58+
int main(int argc, char **argv)
59+
{
60+
double *coeffs = NULL;
61+
long double complex *s0 = NULL;
62+
unsigned int degree = 0;
63+
unsigned int n, i;
64+
65+
if (argc < 2)
66+
{
67+
printf("Please pass the coefficients of the polynomial as commandline arguments.\n");
68+
return 0;
69+
}
70+
71+
degree = argc - 1; /*< detected polynomial degree */
72+
coeffs = (double *)malloc(degree * sizeof(double)); /**< store all input coefficients */
73+
s0 = (long double complex *)malloc((degree - 1) * sizeof(long double complex)); /**< number of roots = degree-1 */
74+
75+
/* initialize random seed: */
76+
srand(time(NULL));
77+
78+
if (!coeffs || !s0)
79+
{
80+
perror("Unable to allocate memory!");
81+
if (coeffs)
82+
free(coeffs);
83+
if (s0)
84+
free(s0);
85+
return EXIT_FAILURE;
86+
}
87+
88+
#if defined(DEBUG) || !defined(NDEBUG)
89+
/**
90+
* store intermediate values to a CSV file
91+
**/
92+
FILE *log_file = fopen("durand_kerner.log.csv", "wt");
93+
if (!log_file)
94+
{
95+
perror("Unable to create a storage log file!");
96+
free(coeffs);
97+
free(s0);
98+
return EXIT_FAILURE;
99+
}
100+
fprintf(log_file, "iter#,");
101+
#endif
102+
103+
printf("Computing the roots for:\n\t");
104+
for (n = 0; n < degree; n++)
105+
{
106+
coeffs[n] = strtod(argv[n + 1], NULL);
107+
if (n < degree - 1 && coeffs[n] != 0)
108+
printf("(%g) x^%d + ", coeffs[n], degree - n - 1);
109+
else if (coeffs[n] != 0)
110+
printf("(%g) x^%d = 0\n", coeffs[n], degree - n - 1);
111+
112+
double tmp;
113+
if (n > 0)
114+
coeffs[n] /= tmp; /* numerical errors less when the first coefficient is "1" */
115+
else
116+
{
117+
tmp = coeffs[0];
118+
coeffs[0] = 1;
119+
}
120+
121+
/* initialize root approximations with random values */
122+
if (n < degree - 1)
123+
{
124+
s0[n] = (long double)rand() + (long double)rand() * I;
125+
#if defined(DEBUG) || !defined(NDEBUG)
126+
fprintf(log_file, "root_%d,", n);
127+
#endif
128+
}
129+
}
130+
131+
#if defined(DEBUG) || !defined(NDEBUG)
132+
fprintf(log_file, "avg. correction");
133+
fprintf(log_file, "\n0,");
134+
for (n = 0; n < degree - 1; n++)
135+
fprintf(log_file, "%s,", complex_str(s0[n]));
136+
#endif
137+
138+
double tol_condition = 1;
139+
unsigned long iter = 0;
140+
141+
while (!check_termination(tol_condition) && iter < INT_MAX)
142+
{
143+
long double complex delta = 0;
144+
tol_condition = 0;
145+
iter++;
146+
147+
#if defined(DEBUG) || !defined(NDEBUG)
148+
fprintf(log_file, "\n%ld,", iter);
149+
#endif
150+
151+
for (n = 0; n < degree - 1; n++)
152+
{
153+
long double complex numerator = function(coeffs, degree, s0[n]);
154+
long double complex denominator = 1.0;
155+
for (i = 0; i < degree - 1; i++)
156+
if (i != n)
157+
denominator *= s0[n] - s0[i];
158+
159+
delta = numerator / denominator;
160+
161+
if (isnan(cabsl(delta)) || isinf(cabsl(delta)))
162+
{
163+
printf("\n\nOverflow/underrun error - got value = %Lg", cabsl(delta));
164+
goto end;
165+
}
166+
167+
s0[n] -= delta;
168+
169+
tol_condition = fmaxl(tol_condition, fabsl(cabsl(delta)));
170+
171+
#if defined(DEBUG) || !defined(NDEBUG)
172+
fprintf(log_file, "%s,", complex_str(s0[n]));
173+
#endif
174+
}
175+
// tol_condition /= (degree - 1);
176+
177+
#if defined(DEBUG) || !defined(NDEBUG)
178+
if (iter % 500 == 0)
179+
{
180+
printf("Iter: %lu\t", iter);
181+
for (n = 0; n < degree - 1; n++)
182+
printf("\t%s", complex_str(s0[n]));
183+
printf("\t\tabsolute average change: %.4g\n", tol_condition);
184+
}
185+
186+
fprintf(log_file, "%.4g", tol_condition);
187+
#endif
188+
}
189+
end:
190+
191+
#if defined(DEBUG) || !defined(NDEBUG)
192+
fclose(log_file);
193+
#endif
194+
195+
printf("\nIterations: %lu\n", iter);
196+
for (n = 0; n < degree - 1; n++)
197+
printf("\t%s\n", complex_str(s0[n]));
198+
printf("absolute average change: %.4g\n", tol_condition);
199+
200+
free(coeffs);
201+
free(s0);
202+
203+
return 0;
204+
}

‎numerical_methods/qr_decomposition.c

Lines changed: 153 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,153 @@
1+
#include <stdio.h>
2+
#include <math.h>
3+
#include <stdlib.h>
4+
5+
#define ROWS 4
6+
#define COLUMNS 3
7+
8+
double A[ROWS][COLUMNS] = {
9+
{1, 2, 3},
10+
{3, 6, 5},
11+
{5, 2, 8},
12+
{8, 9, 3}};
13+
14+
void print_matrix(double A[][COLUMNS], int M, int N)
15+
{
16+
for (int row = 0; row < M; row++)
17+
{
18+
for (int col = 0; col < N; col++)
19+
printf("% 9.3g\t", A[row][col]);
20+
putchar('\n');
21+
}
22+
putchar('\n');
23+
}
24+
25+
void print_2d(double **A, int M, int N)
26+
{
27+
for (int row = 0; row < M; row++)
28+
{
29+
for (int col = 0; col < N; col++)
30+
printf("% 9.3g\t", A[row][col]);
31+
putchar('\n');
32+
}
33+
putchar('\n');
34+
}
35+
36+
double vector_dot(double *a, double *b, int L)
37+
{
38+
double mag = 0.f;
39+
for (int i = 0; i < L; i++)
40+
mag += a[i] * b[i];
41+
42+
return mag;
43+
}
44+
45+
double vector_mag(double *vector, int L)
46+
{
47+
double dot = vector_dot(vector, vector, L);
48+
return sqrt(dot);
49+
}
50+
51+
double *vector_proj(double *a, double *b, double *out, int L)
52+
{
53+
double num = vector_dot(a, b, L);
54+
double deno = vector_dot(b, b, L);
55+
for (int i = 0; i < L; i++)
56+
out[i] = num * b[i] / deno;
57+
58+
return out;
59+
}
60+
61+
double *vector_sub(double *a, double *b, double *out, int L)
62+
{
63+
for (int i = 0; i < L; i++)
64+
out[i] = a[i] - b[i];
65+
66+
return out;
67+
}
68+
69+
void qr_decompose(double A[][COLUMNS], double **Q, double **R, int M, int N)
70+
{
71+
double *col_vector = (double *)malloc(M * sizeof(double));
72+
double *col_vector2 = (double *)malloc(M * sizeof(double));
73+
double *tmp_vector = (double *)malloc(M * sizeof(double));
74+
for (int i = 0; i < N; i++) /* for each column => R is a square matrix of NxN */
75+
{
76+
for (int j = 0; j < i; j++) /* second dimension of column */
77+
R[i][j] = 0.; /* make R upper triangular */
78+
79+
/* get corresponding Q vector */
80+
for (int j = 0; j < M; j++)
81+
{
82+
tmp_vector[j] = A[j][i]; /* accumulator for uk */
83+
col_vector[j] = A[j][i];
84+
}
85+
for (int j = 0; j < i; j++)
86+
{
87+
for (int k = 0; k < M; k++)
88+
col_vector2[k] = Q[k][j];
89+
vector_proj(col_vector, col_vector2, col_vector2, M);
90+
vector_sub(tmp_vector, col_vector2, tmp_vector, M);
91+
}
92+
double mag = vector_mag(tmp_vector, M);
93+
for (int j = 0; j < M; j++)
94+
Q[j][i] = tmp_vector[j] / mag;
95+
96+
/* compute upper triangular values of R */
97+
for (int kk = 0; kk < M; kk++)
98+
col_vector[kk] = Q[kk][i];
99+
for (int k = i; k < N; k++)
100+
{
101+
for (int kk = 0; kk < M; kk++)
102+
col_vector2[kk] = A[kk][k];
103+
R[i][k] = vector_dot(col_vector, col_vector2, M);
104+
}
105+
}
106+
107+
free(col_vector);
108+
free(col_vector2);
109+
free(tmp_vector);
110+
}
111+
112+
int main(void)
113+
{
114+
// double A[][COLUMNS] = {
115+
// {1, -1, 4},
116+
// {1, 4, -2},
117+
// {1, 4, 2},
118+
// {1, -1, 0}};
119+
120+
print_matrix(A, ROWS, COLUMNS);
121+
122+
double **R = (double **)malloc(sizeof(double) * COLUMNS * COLUMNS);
123+
double **Q = (double **)malloc(sizeof(double) * ROWS * COLUMNS);
124+
if (!Q || !R)
125+
{
126+
perror("Unable to allocate memory for Q & R!");
127+
return -1;
128+
}
129+
for (int i = 0; i < ROWS; i++)
130+
{
131+
R[i] = (double *)malloc(sizeof(double) * COLUMNS);
132+
Q[i] = (double *)malloc(sizeof(double) * COLUMNS);
133+
if (!Q[i] || !R[i])
134+
{
135+
perror("Unable to allocate memory for Q & R.");
136+
return -1;
137+
}
138+
}
139+
140+
qr_decompose(A, Q, R, ROWS, COLUMNS);
141+
142+
print_2d(R, ROWS, COLUMNS);
143+
print_2d(Q, ROWS, COLUMNS);
144+
145+
for (int i = 0; i < ROWS; i++)
146+
{
147+
free(R[i]);
148+
free(Q[i]);
149+
}
150+
free(R);
151+
free(Q);
152+
return 0;
153+
}

0 commit comments

Comments
 (0)