From 8f06a164cb7cff853a0263d35bfaddc298e27147 Mon Sep 17 00:00:00 2001 From: danny0926 <42000puzzles@gmail.com> Date: Mon, 7 Mar 2022 15:14:35 +0800 Subject: [PATCH 01/11] Update 5_practice.c --- 5_practice.c | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/5_practice.c b/5_practice.c index 70d6651..1c46a0b 100644 --- a/5_practice.c +++ b/5_practice.c @@ -9,18 +9,19 @@ int main() { time_t t; // Input the number N - + scanf("%d", &N); // Locate the memory for x, yr, yi; + x = (double*)malloc(sizeof(double)); + yr = (double*)malloc(sizeof(double)); + yi = (double*)malloc(sizeof(double)); // Initial setting for x, for example, x[k] = k - - - + for(k = 0; k < N; k++) x[k] = k; t = clock(); // yr[n]+i*yi[n] = sum(exp(-i 2 Pi k n / N)*x[k], k=0..N-1), n=0..N-1 - + @@ -32,7 +33,9 @@ int main() { printf("%d ms for discrete Fourier Transform of %d elements\n", t, N); // free the memory located by x, yr, yi - + free(x); + free(yr); + free(yi); return 100; From 45a1e97e0c7a5245a8a873812343be98a4a22d4c Mon Sep 17 00:00:00 2001 From: danny0926 <42000puzzles@gmail.com> Date: Mon, 7 Mar 2022 15:48:34 +0800 Subject: [PATCH 02/11] Update 5_practice.c --- 5_practice.c | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/5_practice.c b/5_practice.c index 1c46a0b..9596670 100644 --- a/5_practice.c +++ b/5_practice.c @@ -12,21 +12,24 @@ int main() { scanf("%d", &N); // Locate the memory for x, yr, yi; - x = (double*)malloc(sizeof(double)); - yr = (double*)malloc(sizeof(double)); - yi = (double*)malloc(sizeof(double)); + x = (double*)malloc(N*sizeof(double)); + yr = (double*)malloc(N*sizeof(double)); + yi = (double*)malloc(N*sizeof(double)); // Initial setting for x, for example, x[k] = k for(k = 0; k < N; k++) x[k] = k; t = clock(); // yr[n]+i*yi[n] = sum(exp(-i 2 Pi k n / N)*x[k], k=0..N-1), n=0..N-1 - - - - - - + for(n = 0; n < N; n++){ + yr[n] = 0.0; + yi[n] = 0.0; + for(k = 0; k < N; k++){ + double theta = (-2*M_PI*k*n)/N; + yr[n] = yr[n] + cos(theta)*x[k]; + yi[n] = yi[n] + sin(theta)*x[k]; + } + } // output the results t = clock() - t; From d6ccfa1cc156dc4334567f3b1e7e141d9b14a721 Mon Sep 17 00:00:00 2001 From: danny0926 <42000puzzles@gmail.com> Date: Thu, 17 Mar 2022 20:52:20 +0800 Subject: [PATCH 03/11] Update 5_practice.c MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 加速 --- 5_practice.c | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/5_practice.c b/5_practice.c index 9596670..a9d57f5 100644 --- a/5_practice.c +++ b/5_practice.c @@ -5,7 +5,7 @@ int main() { // Declare all the variables int k, n, N; - double *x, *yr, *yi; + double *x, *yr, *yi, a, b, an, bn, temp, P, c, s; time_t t; // Input the number N @@ -21,14 +21,26 @@ int main() { t = clock(); // yr[n]+i*yi[n] = sum(exp(-i 2 Pi k n / N)*x[k], k=0..N-1), n=0..N-1 + a = cos(2*M_PI/N); b = -sin(2*M_PI/N); for(n = 0; n < N; n++){ yr[n] = 0.0; yi[n] = 0.0; + P = 2*M_PI/N; + an = 1; bn = 0; + c = 1; s = 0; for(k = 0; k < N; k++){ - double theta = (-2*M_PI*k*n)/N; - yr[n] = yr[n] + cos(theta)*x[k]; - yi[n] = yi[n] + sin(theta)*x[k]; + yr[n] += c*x[k]; + yi[n] -= s*x[k]; + temp = c * an - s * bn; + s = c * bn + s * an; + c = temp; +// double theta = (-2*M_PI*k*n)/N; +// yr[n] = yr[n] + cos(theta)*x[k]; +// yi[n] = yi[n] + sin(theta)*x[k]; } + temp = an * a - bn * b; + bn = an * b + bn * a; + an = temp; } // output the results From 7bb44fa276d903564280589d1cabfd5c8500609e Mon Sep 17 00:00:00 2001 From: danny0926 <42000puzzles@gmail.com> Date: Thu, 17 Mar 2022 21:42:16 +0800 Subject: [PATCH 04/11] Create 6_sorting.c --- 6_sorting.c | 91 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 91 insertions(+) create mode 100644 6_sorting.c diff --git a/6_sorting.c b/6_sorting.c new file mode 100644 index 0000000..4beda70 --- /dev/null +++ b/6_sorting.c @@ -0,0 +1,91 @@ +#include // for printf function +#include // for memory allocation +#include // for time calculation +#include // for sine and cosine functions +int main() { + // Declare all the variables + int k, m, n, N; + double *x, *y, z, p; + time_t t; + + // Input the number N + printf("Input N: "); + scanf("%d",&N); + + // Locate the memory for x and y; + x = (double *) malloc(N*sizeof(double)); + y = (double *) malloc(N*sizeof(double)); + + // Initial setting for x, for example, x[k] = 1.0*rand()/RAND_MAX + srand( time(NULL) ); + for(k=0;k x[k]) { + z = x[n]; + x[n] = x[k]; + x[k] = z; + } + } + } + t = clock() - t; + + // print y, x, and time + printf("Sorting %d elements: %f s\n", N, 1.0*t/CLOCKS_PER_SEC); + if(N<10) { + printf("y \t\t x\n"); + for(k=0;k %d\n",k,n); + z = x[n]; + x[n] = x[k]; + x[k] = z; + n++; + } + } + z = x[n]; + x[n] = x[N-1]; + x[N-1] = z; + + // x[0] is at the correct location! + if(N<10) { + printf("y \t\t x\n"); + for(k=0;k Date: Sat, 19 Mar 2022 10:06:36 +0800 Subject: [PATCH 05/11] Update 5_practice.c --- 5_practice.c | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/5_practice.c b/5_practice.c index a9d57f5..6151566 100644 --- a/5_practice.c +++ b/5_practice.c @@ -21,26 +21,28 @@ int main() { t = clock(); // yr[n]+i*yi[n] = sum(exp(-i 2 Pi k n / N)*x[k], k=0..N-1), n=0..N-1 - a = cos(2*M_PI/N); b = -sin(2*M_PI/N); + + a = cos(2*M_PI/N); b = -sin(2*M_PI/N);//***** for(n = 0; n < N; n++){ yr[n] = 0.0; yi[n] = 0.0; - P = 2*M_PI/N; - an = 1; bn = 0; - c = 1; s = 0; + P = 2*M_PI/N;//***** + an = 1; bn = 0;//***** + c = 1; s = 0;//***** for(k = 0; k < N; k++){ - yr[n] += c*x[k]; - yi[n] -= s*x[k]; - temp = c * an - s * bn; - s = c * bn + s * an; - c = temp; + yr[n] += c*x[k];//***** + yi[n] -= s*x[k];//***** + temp = c * an - s * bn;//***** + s = c * bn + s * an;//***** + c = temp;//***** + // double theta = (-2*M_PI*k*n)/N; // yr[n] = yr[n] + cos(theta)*x[k]; // yi[n] = yi[n] + sin(theta)*x[k]; } - temp = an * a - bn * b; - bn = an * b + bn * a; - an = temp; + temp = an * a - bn * b;//***** + bn = an * b + bn * a;//***** + an = temp; //***** } // output the results From 880114e7cf1e101fee51ebe52b5e455f0ace70c5 Mon Sep 17 00:00:00 2001 From: danny0926 <42000puzzles@gmail.com> Date: Sat, 19 Mar 2022 12:20:27 +0800 Subject: [PATCH 06/11] Update 6_sorting.c --- 6_sorting.c | 54 ++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/6_sorting.c b/6_sorting.c index 4beda70..8aaa6e9 100644 --- a/6_sorting.c +++ b/6_sorting.c @@ -2,6 +2,12 @@ #include // for memory allocation #include // for time calculation #include // for sine and cosine functions + + + +int quick_sort(double*, int, int); +void SWAP(double, double); + int main() { // Declare all the variables int k, m, n, N; @@ -87,5 +93,51 @@ int main() { } int quick_sort(double *x, int L, int R){ - + double temp; + if(L < R){ + double pk = x[L]; + int i = L, j = R + 1; + while(i < j){ + while(x[i] < pk) i++; + while(x[j] > pk) j--; + if(i < j){ + temp = x[i]; + x[i] = x[j]; + x[j] = temp; + } + } + temp = x[L]; + x[L] = x[j]; + x[j] = temp; +// SWAP(x[L], x[j]); + quick_sort(x, L, j - 1); + quick_sort(x, j + 1, R); + } } + +//void SWAP(double X, double Y){ +// X ^= Y; +// Y ^= X; +// X ^= Y; +//} + + + + + + + + + + + + + + + + + + + + + From 97a290d2c15a35b4ec8f59b73bbacc40a15d294b Mon Sep 17 00:00:00 2001 From: danny0926 <42000puzzles@gmail.com> Date: Sat, 19 Mar 2022 12:20:31 +0800 Subject: [PATCH 07/11] Create test.c --- test.c | 147 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 147 insertions(+) create mode 100644 test.c diff --git a/test.c b/test.c new file mode 100644 index 0000000..aaa85ff --- /dev/null +++ b/test.c @@ -0,0 +1,147 @@ +#include // for printf function +#include // for memory allocation +#include // for time calculation +#include // for sine and cosine functions + + + +void quick_sort(double*, int, int); +int Partition(double*, int, int); +int main() { + // Declare all the variables + int k, m, n, N; + double *x, *y, z, p; + time_t t; + + // Input the number N + printf("Input N: "); + scanf("%d",&N); + + // Locate the memory for x and y; + x = (double *) malloc(N*sizeof(double)); + y = (double *) malloc(N*sizeof(double)); + + // Initial setting for x, for example, x[k] = 1.0*rand()/RAND_MAX + srand( time(NULL) ); + for(k=0;k x[k]) { + z = x[n]; + x[n] = x[k]; + x[k] = z; + } + } + } + t = clock() - t; + + // print y, x, and time + printf("Sorting %d elements: %f s\n", N, 1.0*t/CLOCKS_PER_SEC); + + printf("bubble sort : \n"); + for(k=0;k %d\n",k,n); +// z = x[n]; +// x[n] = x[k]; +// x[k] = z; +// n++; +// } +// } +// z = x[n]; +// x[n] = x[N-1]; +// x[N-1] = z; + + // x[0] is at the correct location! +// if(N<10) { +// printf("y \t\t x\n"); +// for(k=0;k Date: Sat, 19 Mar 2022 16:56:54 +0800 Subject: [PATCH 08/11] Update test.c --- test.c | 105 +++++++++++++++++++++++++++++++++++++++------------------ 1 file changed, 73 insertions(+), 32 deletions(-) diff --git a/test.c b/test.c index aaa85ff..78c0aff 100644 --- a/test.c +++ b/test.c @@ -7,6 +7,9 @@ void quick_sort(double*, int, int); int Partition(double*, int, int); +double Select(double*, int, int, int); +double Find_Med(double*, int, int); +double QD(double*, int, int); int main() { // Declare all the variables int k, m, n, N; @@ -26,7 +29,11 @@ int main() { for(k=0;k %d\n",k,n); -// z = x[n]; -// x[n] = x[k]; -// x[k] = z; -// n++; -// } -// } -// z = x[n]; -// x[n] = x[N-1]; -// x[N-1] = z; - - // x[0] is at the correct location! -// if(N<10) { -// printf("y \t\t x\n"); -// for(k=0;k Date: Sat, 19 Mar 2022 19:34:49 +0800 Subject: [PATCH 09/11] Update test.c --- test.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test.c b/test.c index 78c0aff..3de08cd 100644 --- a/test.c +++ b/test.c @@ -5,11 +5,11 @@ -void quick_sort(double*, int, int); -int Partition(double*, int, int); -double Select(double*, int, int, int); -double Find_Med(double*, int, int); -double QD(double*, int, int); +void quick_sort(double*, int, int); //quick sort +int Partition(double*, int, int); //sub algo of quick sort +double Select(double*, int, int, int); // find the i-th number +double Find_Med(double*, int, int); //find the median +double QD(double*, int, int); // find the Quartile Deviation int main() { // Declare all the variables int k, m, n, N; From 748126d9033e6f702e16855736cd471344c624cf Mon Sep 17 00:00:00 2001 From: danny0926 <42000puzzles@gmail.com> Date: Sat, 19 Mar 2022 19:35:40 +0800 Subject: [PATCH 10/11] Update 6_sorting.c --- 6_sorting.c | 167 +++++++++++++++++++++++++++++++++------------------- 1 file changed, 106 insertions(+), 61 deletions(-) diff --git a/6_sorting.c b/6_sorting.c index 8aaa6e9..3de08cd 100644 --- a/6_sorting.c +++ b/6_sorting.c @@ -5,9 +5,11 @@ -int quick_sort(double*, int, int); -void SWAP(double, double); - +void quick_sort(double*, int, int); //quick sort +int Partition(double*, int, int); //sub algo of quick sort +double Select(double*, int, int, int); // find the i-th number +double Find_Med(double*, int, int); //find the median +double QD(double*, int, int); // find the Quartile Deviation int main() { // Declare all the variables int k, m, n, N; @@ -27,7 +29,11 @@ int main() { for(k=0;k %d\n",k,n); - z = x[n]; - x[n] = x[k]; - x[k] = z; - n++; - } - } - z = x[n]; - x[n] = x[N-1]; - x[N-1] = z; - // x[0] is at the correct location! - if(N<10) { - printf("y \t\t x\n"); - for(k=0;k pk) j--; - if(i < j){ - temp = x[i]; - x[i] = x[j]; - x[j] = temp; - } - } - temp = x[L]; - x[L] = x[j]; - x[j] = temp; -// SWAP(x[L], x[j]); - quick_sort(x, L, j - 1); - quick_sort(x, j + 1, R); + if(P < r){ + int q = Partition(x, P, r); + quick_sort(x, P, q - 1); + quick_sort(x, q + 1, r); } } -//void SWAP(double X, double Y){ -// X ^= Y; -// Y ^= X; -// X ^= Y; -//} - - +int Partition(double *x, int P, int r){ + double a = x[r]; + int i = P - 1, k; + for(k = P; k < r; k++){ + if(x[k] <= a){ + i++; + double temp; + temp = x[i]; + x[i] = x[k]; + x[k] = temp; + } + } + double temp; + temp = x[i + 1]; + x[i + 1] = x[r]; + x[r] = temp; + return (i + 1); +} +double Select(double *x, int P, int r, int i){ + if(P < r){ + int q = Partition(x, P, r); + int k = q - P + 1; + if(i == k) return x[q]; + else if(i < k) return Select(x, P, q - 1, i); + else return Select(x, q + 1, r, i - k); + } + else{ + return x[r]; + } +} +double Find_Med(double *x, int P, int r){ + int number = r - P + 1; + if(number % 2 == 0){ + double right = Select(x, P, r, number/2 + 1); + double left = Select(x, P, r, number/2); + return (right + left)/2; + } + else{ + return Select(x, P, r, number/2 + 1); + } +} +double QD(double *x, int P, int r){ // ¥|¤À¦ì®t Quartile Deviation, QD + int number = r - P + 1; + if(number/2 == 0){ + double Q1 = Find_Med(x, P, number/2 - 1); + printf("Q1 = %f\n", Q1); + double Q3 = Find_Med(x, number/2, r); + printf("Q3 = %f\n", Q3); + return Q3 - Q1; + } + else{ + double Q1 = Find_Med(x, P, number/2 - 1); + printf("Q1 = %f\n", Q1); + double Q3 = Find_Med(x, (number + 1)/2, r); + printf("Q3 = %f\n", Q3); + return Q3 - Q1; + } +} From 9bb446ce4b5f22a759dbbd3c3bdfb53f5de1a4aa Mon Sep 17 00:00:00 2001 From: danny0926 <42000puzzles@gmail.com> Date: Fri, 30 Sep 2022 14:21:34 +0800 Subject: [PATCH 11/11] 9/30 --- 2-D-Multigrid.c | 266 +++++++++ 2DpoissonSolver.c | 425 ++++++++++++++ ...257\346\266\210\345\216\273\346\263\225.c" | 90 +++ FFT.c | 98 ++++ FFT2.c | 517 +++++++++++++++++ MGsolver.c | 278 +++++++++ MGsolver2D.c | 201 +++++++ Plot.dat | 0 base_p.c | 65 +++ inverseFFT.c | 319 +++++++++++ radix2.c | 64 +++ selcet_number.c | 188 +++++++ test.c | 527 +++++++++++++----- test2.c | 201 +++++++ test3.c | 29 + ...253\351\200\237\344\271\230\346\263\225.c" | 291 ++++++++++ 16 files changed, 3414 insertions(+), 145 deletions(-) create mode 100644 2-D-Multigrid.c create mode 100644 2DpoissonSolver.c create mode 100644 "3\346\242\235\347\232\204\351\253\230\346\226\257\346\266\210\345\216\273\346\263\225.c" create mode 100644 FFT.c create mode 100644 FFT2.c create mode 100644 MGsolver.c create mode 100644 MGsolver2D.c create mode 100644 Plot.dat create mode 100644 base_p.c create mode 100644 inverseFFT.c create mode 100644 radix2.c create mode 100644 selcet_number.c create mode 100644 test2.c create mode 100644 test3.c create mode 100644 "\345\277\253\351\200\237\344\271\230\346\263\225.c" diff --git a/2-D-Multigrid.c b/2-D-Multigrid.c new file mode 100644 index 0000000..8bee687 --- /dev/null +++ b/2-D-Multigrid.c @@ -0,0 +1,266 @@ +// Last change: Abdul Hannan 03 Jan 2019 +#include +#include +#include + +//Global Variables +int prolong(), power(int, int);; +void rest(), add(); +double T[300][300], res[300][300], err[300][300], Tcalc, Tn, h, BE, e=1e-10, er, x, y, x2, y2, r, pi=3.14159265359; +int i, j, k, l, w=2, c, m=256, a, init, cyc, iter=1, n=256, f, f2, ln, lmin, level=0; +float lmax; + +int main() +{ +//Declaration (Local Variables) +float le, b, v; +char bc, sch, choice; +FILE *fp, *fp2, *fp3; + +//File opening +fp = fopen("Poisson Log.csv", "a"); //Data logging +if(fp == NULL){ + printf("Couldn't open file\n"); + return 1; +} + +fp2 = fopen("Plot.dat", "w"); //For plotting results in Tecplot +if(fp2 == NULL){ + printf("Couldn't open file\n"); + return 1; +} + +fp3 = fopen("Residual.csv", "w"); //For recording observations +if(fp3 == NULL){ + printf("Couldn't open file\n"); + return 1; +} + +/*Inputs_________________________________________________________________________________________________________________*/ +x = 1.0/m; +x2 = x*x; +y = 1.0/n; +y2 = y*y; +printf("x = %f", x); +printf("\ty = %f", y); + +x: printf("\n\nChoose an AMG Scheme\nEnter a character: [V]\t[W]\n"); +scanf("%c", &sch); +if (sch == 'V') +{ + fprintf(fp, "V,"); +} +if (sch=='W') +{ + fprintf(fp, "W,"); + printf("Enter the grid level at W point:\t"); + scanf("%d", &lmin); + fprintf(fp, "%d,", lmin); +} +printf("\nEnter the no. of levels: "); +scanf("%d", &ln); +fprintf(fp, "%d,", ln); +printf("\nEnter the no. of inner iterations: "); +scanf("%d", &init); +printf("\nEnter the no. of %c cycles: ", sch); +scanf("%d", &cyc); +fprintf(fp, "%d,", cyc); +/*_________________________________________________________________________________________________________________________*/ + +/*Solution Initialization__________________________________________________________________________________________________*/ +for (i=0; i<=n; i++) +{ + for (j=0; j<=m; j++) + { + T[i][j]=0; + } +} +/*_________________________________________________________________________________________________________________________*/ + +//Solution +clock_t begin = clock(); +while (1) +{ + BE=0; +/*Gauss-Seidel____________________________________________________________________________________________________________*/ + for (i= 1; i0) + { + prolong (); + } + add(); + iter++; + } + } + + else if ((sch == 'W') || (sch == 'w')) + { + while(iter<=cyc) + { + k=n/2; + l=m/2; + level++; + rest(); + k=k*2; + l=l*2; + while (level>lmin) + { + prolong (); + } + k=k/2; + l=l/2; + rest(); + k=k*2; + l=l*2; + while (level>0) + { + prolong (); + } + add(); + iter++; + } + } +/*_________________________________________________________________________________________________________________________*/ + + fprintf(fp3, "%g\n", BE); + + if (BEBE) + BE = fabs(res[i][j]); + } + } +} + +void rest() +{ + while (level<=ln) + { + f = power(w ,level); + f2= f*f; + for (a=1; a +#include +#include +#include +#include +#include + + +void Exact_Discretization(double **A, int N){ + int i,j,k; + for(i=0;i0) A[k][k-(N-1)] = 1; + if(i>0) A[k][k-1] = 1; + if(i= 0; --i){ + for(j = i+1; j < N; ++j){ + x[i] -= M[i][j]*x[j]; + } + x[i] /= M[i][i]; + } + + // free memory + free(M[0]); + free(M); +} + +double Residual(double *r, double **A, double *x, double *b, int N){ + //compute r = b-Ax + int i, j; + double v = 0.0; + for(i = 0; i < N; ++i){ + r[i] = b[i]; + + for(j = 0; j < N; ++j){ + r[i] -= A[i][j]*x[j]; + } + + if(fabs(r[i]) > v) v = fabs(r[i]); + } + return v; +} + +double error(double *x, double *u, int N){ + int i; + double e, v = 0.0; + + for(i = 0; i < N; ++i){ + e = fabs(x[i]-u[i]); + if(e > v) v = e; + } + return v; +} + +void FFT_p(complex *y, complex *x, int N){ + int p; + if(N%2 == 0) p = 2; + else if(N%3 == 0) p = 3; + else if(N%5 == 0) p = 5; + else p = N; + + + complex *Nx, *Ny; + int i, j, k, bias, idx; + + bias = N/p; //bias = N/p : °¾²¾¶q + + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(j = 0; j < p; ++j){ + for(k = 0; k < bias; ++k){ + Nx[k + j*bias] = x[p*k + j]; + } + } + + // F[p][p] + +// complex F[p][p]; + + complex **F; + F = (complex**)malloc(p*sizeof(complex*)); + for(i = 0; i < p; ++i){ + F[i] = (complex*)malloc(p*sizeof(complex)); + } + + for(i = 0; i < p; ++i){ + F[0][i] = 1; + F[i][0] = 1; + } + for(i = 1; i < p/2+1; ++i){ + for(j = 1; j < p; ++j){ + F[i][j] = cexp(-2*M_PI*i*j*I/p); + } + } + for(i = p/2+1; i < p; ++i){ //¦³¤Ï¹ïºÙ + for(j = 1; j < p; ++j){ + F[i][j] = F[p-i][p-j]; + } + } + + + //recursive + if(N == p){ + for(i = 0; i < p; ++i){ + y[i] = 0; + for(j = 0; j < p; ++j){ + y[i] += F[i][j]*x[j]; + } + } + } + else{ + + for(j = 0; j < p; ++j){ + FFT_p(Ny + j*bias, Nx + j*bias, bias); + } + + for(i = 0; i < bias; ++i){ + for(k = 0; k < p; ++k){ + Ny[i + k*bias] *= cexp(-2*M_PI*k*i*I/N); + } + + for(j = 0; j < p; ++j){ + idx = i + j*bias; + y[idx] = 0; + for(k = 0; k < p; ++k){ + y[idx] += F[j][k]*Ny[i + k*bias]; + } + + } + } + + } + + free(Nx); free(Ny); +} + +void DST(double *y, double *x, int N){ + //y is the answer + + int i; + N = 2*N + 2; + complex *Nx, *Ny; + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(i = 0; i < N; ++i){ + Nx[i] = 0.0; + } + + for(i = 0; i < N/2-1; ++i){ + Nx[i+1] = x[i]; + Nx[(i+1) + N/2] = -x[N/2-(i+2)]; + } + + + FFT_p(Ny, Nx, N); + + + for(i = 0; i < N/2-1; ++i){ + y[i] = -0.5*cimag(Ny[i+1]); + } + + free(Nx); free(Ny); + +} + +void iDST(double *y, double *x, int N){ + //y is the answer + int i; + N = 2*N + 2; + complex *Nx, *Ny; + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(i = 0; i < N; ++i){ + Nx[i] = 0.0; + } + + for(i = 0; i < N/2-1; ++i){ + Nx[i+1] = x[i]; + Nx[(i+1) + N/2] = -x[N/2-(i+2)]; + } + + + FFT_p(Ny, Nx, N); + + + for(i = 0; i < N/2-1; ++i){ + y[i] = (-2.0/N)*cimag(Ny[i+1]); + } + + free(Nx); free(Ny); +} + +void Transpose(double **A, int N){ + //transpose a NxN matrix + int i , j; + double temp; + + for(i = 0; i < N; ++i){ + for(j = i+1; j < N; ++j){ + temp = A[i][j]; + A[i][j] = A[j][i]; + A[j][i] = temp; + } + } +} + +void DST2D(double **X, int N){ + int i; + + for(i = 0; i < N; ++i){ + DST(X[i], X[i], N); + } + + Transpose(X, N); + + for(i = 0; i < N; ++i){ + DST(X[i], X[i], N); + } + + Transpose(X, N); +} + +void iDST2D(double **X, int N){ + int i; + + for(i = 0; i < N; ++i){ + iDST(X[i], X[i], N); + } + + Transpose(X, N); + + for(i = 0; i < N; ++i){ + iDST(X[i], X[i], N); + } + + Transpose(X, N); +} + +void FastPoissonSolver(double **F, int N){ + int i, j; + double h; + + h = 1.0/(N+1); + DST2D(F, N); + + + + for(i = 0; i < N; ++i){ + for(j = 0; j < N; ++j){ + F[i][j]= F[i][j]/((2*(cos(M_PI*(i+1)*h)-1))+(2*(cos(M_PI*(j+1)*h)-1))); + } + } + + + iDST2D(F, N); + +} + +int main(){ + //solve AU = F, where r = b - Ax is the residue + // U is unknown, F is known + int i, j, k, N, M; + double **A, *x, *u, **U, *b, **F, *r; + time_t t; + + N = 32; + + M = (N-1)*(N-1); + //create memory, ±N2ºûarray¶}¦¨¤@¦C + A = (double**) malloc(M*sizeof(double*)); + A[0] = (double*) malloc(M*M*sizeof(double)); + for(i = 1; i < M; ++i){ + A[i] = A[i-1] + M; + } + x = (double*) malloc(M*sizeof(double)); + r = (double*) malloc(M*sizeof(double)); + + + b = (double*) malloc(M*sizeof(double)); + F = (double**)malloc((N-1)*sizeof(double*)); + F[0] = b; + for(i = 1; i < N-1; ++i){ + F[i] = F[i-1] + N - 1; + } + + + u = (double*) malloc(M*sizeof(double)); + U = (double**) malloc((N-1)*sizeof(double*)); + U[0] = u; + for(i = 1; i < N-1; ++i){ + U[i] = U[i-1] + N - 1; + } + //end of open memory + + + //initialize + for(i = 0; i < M*M; ++i){ + A[0][i] = 0.0; + } + //end of initialize + + //make the matrix A + Exact_Discretization(A, N); + + + //make F + Source(F, N); + + //case1. use Gauss elimination solve + t = clock(); + BandGaussianElimination(A, x, b, M); + + t = clock() - t; + printf("use Gauss elimination: %f sec\n", 1.0*t/CLOCKS_PER_SEC); + printf("the residual : %e\n", Residual(r, A, x, b, M)); + + //calculate the true answer in order to find the error + ExactSolution(U, N); + printf("error is: %e\n", error(x, u, M)); + + + //case 2. use DST fast poisson solver + t = clock(); + FastPoissonSolver(F, N-1); + t = clock() - t; + + printf("use DST fast poisson solver: %f sec\n", 1.0*t/CLOCKS_PER_SEC); + printf("error is: %e\n", error(b, u, M)); + + + //free memory + free(A[0]); + free(A); + free(x); + free(r); + free(b); + free(u); + free(U); + free(F); + + + return 0; +} diff --git "a/3\346\242\235\347\232\204\351\253\230\346\226\257\346\266\210\345\216\273\346\263\225.c" "b/3\346\242\235\347\232\204\351\253\230\346\226\257\346\266\210\345\216\273\346\263\225.c" new file mode 100644 index 0000000..a7b9c24 --- /dev/null +++ "b/3\346\242\235\347\232\204\351\253\230\346\226\257\346\266\210\345\216\273\346\263\225.c" @@ -0,0 +1,90 @@ +#include +#include +//¼e±a¯x°}ªº°ª´µ®ø¥hªk + +float* BandGaussianElimination(float **A, float *B, int M, int N){ + //begin from A[M][M], end to A[N-1][N-1] , doing in this submatrix + //N by N matrix + //solve Ax = B + int i, j, k, n; + float a, b, c; + a = A[0][0]; + b = A[1][0]; + c = A[0][1]; + //a, b, c ¤£¥þ¬Û¦P + //¥Ñ¤W¦Ó¤U + + for(i = 0; i < N-1; ++i){ + a = A[i][i]; + b = A[i+1][i]; + c = A[i][i+1]; + + + A[i+1][i] = 0; + A[i+1][i+1] -= b*c/a; + + B[i+1] -= B[i]*b/a; + + + + } + + + //¥Ñ¤U¦Ó¤W + for(i = N-1; i >= 1; --i){ + a = A[i][i]; + b = A[i-1][i]; + + + A[i-1][i] = 0; + B[i-1] -= B[i]* b/a;; + } + + for(i = 0; i < N; ++i){ + a = A[i][i]; + A[i][i] = 1; + B[i] /= a; + } + + return B; + +} + + + + + + + +int main(){ + int i, j, N; + N = 3; + + float **A = (float**)malloc(N*sizeof(float*)); + for(i = 0; i < N; i++){ + A[i] = (float*)malloc(N*sizeof(float)); + } + + float *b = (float*)malloc(N*sizeof(float)); + + for(i = 0; i < N; i++){ + for(j = 0; j < N; j++){ + scanf("%f", &A[i][j]); + } + } + + for(i = 0; i < N; i++){ + scanf("%f", &b[i]); + } + + float *x = (float*)malloc(N*sizeof(float)); + + x = BandGaussianElimination(A, b, 0, N); + + for(i = 0; i < N; i++){ + printf("%f ", x[i]); + } + + + return 0; +} diff --git a/FFT.c b/FFT.c new file mode 100644 index 0000000..b12e4b3 --- /dev/null +++ b/FFT.c @@ -0,0 +1,98 @@ +#include +#include +#include +#include +int FFT(double *yr, double *yi, double *xr, double *xi, int N){ + + if(N == 2){ + yr[0] = xr[0] + xr[1]; + yi[0] = xi[0] + xi[1]; + yr[1] = xr[0] - xr[1]; + yi[1] = xi[0] - xi[1]; + return 0; + } + else{ + int n, k, c, s; + double *yEr, *yEi, *yOr, *yOi, *xEr, *xEi, *xOr, *xOi; + + xEr = (double*) malloc((N/2)*sizeof(double)); + xEi = (double*) malloc((N/2)*sizeof(double)); + yEr = (double*) malloc((N/2)*sizeof(double)); + yEi = (double*) malloc((N/2)*sizeof(double)); + xOr = (double*) malloc((N/2)*sizeof(double)); + xOi = (double*) malloc((N/2)*sizeof(double)); + yOr = (double*) malloc((N/2)*sizeof(double)); + yOi = (double*) malloc((N/2)*sizeof(double)); + for(n = 0; n < N/2; n++){ + xEr[n] = xr[2*n]; + xEi[n] = xi[2*n]; + xOr[n] = xr[2*n + 1]; + xOi[n] = xi[2*n + 1]; + } + FFT(yEr, yEi, xEr, xEi, N/2); + FFT(yEr, yEi, xEr, xEi, N/2); + + for(k = 0; k < N/2; k++){ + c = cos(2*M_PI*k/N); + s = -sin(2*M_PI*k/N); + + yr[k] = yEr[k] + (yOr[k]*c - yOi[k]*s); + yi[k] = yEr[k] + (yOr[k]*s + yOi[k]*c); + yr[N/2 + k] = yEr[k] - (yOr[k]*c - yOi[k]*s); + yi[N/2 + k] = yEi[k] - (yOr[k]*s + yOi[k]*c); + } + free(xEr); free(xEi); + free(yEr); free(yEi); + free(xOr); free(xOi); + free(yOr); free(yOi); + } +} + + + + + +int main(){ + FILE *fp; + char ch; + int L = 0, N = 1; + int k, n; + double *xr, *xi, *yr, *yi, *zr, *zi, *ur, *ui; + time_t T; + + + fp = fopen("hw6.csv", "r"); + while(!feof(fp)){ + ch = fgetc(fp); + if(ch == '\n') L++; + } + fclose(fp); + + while(N < L) N = N * 2; + N = N/2; + + xr = (double*) malloc(N*sizeof(double)); + xi = (double*) malloc(N*sizeof(double)); + yr = (double*) malloc(N*sizeof(double)); + yi = (double*) malloc(N*sizeof(double)); + zr = (double*) malloc(N*sizeof(double)); + zi = (double*) malloc(N*sizeof(double)); + ur = (double*) malloc(N*sizeof(double)); + ui = (double*) malloc(N*sizeof(double)); + + fp = fopen("hw6.csv", "r"); + for(k = 0; k < N; k++){ + fscanf(fp, "%d %d", &n, &L); + xr[k] = 1.0*L; + xi[k] = 0.0; + } + fclose(fp); + + T = clock(); + FFT(zr, zi, xr, xi, N); + T = clock() - T; + printf("time = %d ms \n", T); + + + return 0; +} diff --git a/FFT2.c b/FFT2.c new file mode 100644 index 0000000..2be14e3 --- /dev/null +++ b/FFT2.c @@ -0,0 +1,517 @@ +#include +#include +#include +#include +#include + + +// FFT_2(complex *y, complex *x, int N)¡G FFT2 by bit reverse +// FFT_3(complex *y, complex *x, int N)¡G FFT3 by bit reverse +// FFT_5(complex *y, complex *x, int N)¡G FFT5 by bit reverse +// FFT2(complex *y, complex *x, int N)¡G FFT with base 2 +// FFT3(complex *y, complex *x, int N)¡G FFT with base 3 +// FFT5(complex *y, complex *x, int N)¡G FFT with base 5 +// FFTp(complex *y, complex *x, int N)¡G FFT with 2^q * 3^q * 5^r + +// FFT_p(complex *y, complex *x, int N) : FFT with 2^q * 3^q * 5^r , and others are FT +// y: put the result after FFT +// x: the input of FFT +// N: the size of input array + +void FFT_2(complex *x, int N){ + int i, j, k, p, M, L, temp; + complex u, w, t; + + M = N/2; + j = 0; + for(i = 1; i < N-1; i++){ + L = M; + while(j >= L){ + j -= L; + L /= 2; + } + j += L; + if(j > i){ +// printf("%d <-> %d\n",i,j); + t = x[i]; + x[i] = x[j]; + x[j] = t; + } + } + j = 1; + while(j < N){ + w = 1.0; + u = cexp(-I*M_PI/j); + for(i = 0; i < j; i++){ + for(k = 0; k < N; k += 2*j){ + x[k+i+j] *= w; + t = x[k + i]; + x[k + i] = x[k + i] + x[k + i + j]; + x[k+i+j] = t - x[k + i + j]; + } + w *= u; + } + j = j << 1; + } + +} + +void FFT_3(complex *x, int N){ + int i, j, k, p, M, L, temp; + complex u, w, t1, t2, t3, w1, w2; + + M = N/3; + j = 0; + for(i = 1; i < N-1; i++){ + L = M; + while(j >= 2*L){ + j -= 2*L; + L /= 3; + } + j += L; + if(j > i){ + printf("%d <-> %d\n",i,j); + t1 = x[i]; + x[i] = x[j]; + x[j] = t1; + } + } + + j = 1; + w1 = cexp(-I*M_PI*2/3); // w4 == w1 + w2 = cexp(-I*M_PI*4/3); + + while(j < N){ + w = 1.0; + u = cexp(-I*2*M_PI/N); + for(i = 0; i < j; i++){ + for(k = 0; k < N; k += 3*j){ + + t1 = x[k+i]; + t2 = x[k+i+j]*w; + t3 = x[k+i+2*j]*cpow(w, 2); + + x[k+i] = t1 + t2 + t3; + x[k+i+j] = t1 + w1*t2 + w2*t3; + x[k+i+2*j] = t1 + w2*t2 + w1*t3; + } + w *= u; + } + j = j*3; + } + +} + +void FFT_5(complex *x, int N){ + int i, j, k, p, M, L, temp; + complex u, w, t1, t2, t3, t4, t5, w1, w2, w3, w4; + + M = N/5; + j = 0; + for(i = 1; i < N-1; i++){ + L = M; + while(j >= 4*L){ + j -= 4*L; + L /= 5; + } + j += L; + if(j > i){ +// printf("%d <-> %d\n",i,j); + t1 = x[i]; + x[i] = x[j]; + x[j] = t1; + } + } + + j = 1; + w1 = cexp(-I*M_PI*2/5); // w16 == w6 == w1 + w2 = cexp(-I*M_PI*4/5); // w12 == w2 + w3 = cexp(-I*M_PI*6/5); // w8 == w3 + w4 = cexp(-I*M_PI*8/5); // w9 == w4 + + while(j < N){ + w = 1.0; + u = cexp(-I*2*M_PI/N); + for(i = 0; i < j; i++){ + for(k = 0; k < N; k += 5*j){ + + t1 = x[k+i]; + t2 = x[k+i+j]*w; + t3 = x[k+i+2*j]*cpow(w, 2); + t4 = x[k+i+3*j]*cpow(w, 3); + t5 = x[k+i+4*j]*cpow(w, 4); + + x[k+i] = t1 + t2 + t3 + t4 + t5; + x[k+i+j] = t1 + w1*t2 + w2*t3 + w3*t4 + w4*t5; + x[k+i+2*j] = t1 + w2*t2 + w4*t3 + w1*t4 + w3*t5; + x[k+i+3*j] = t1 + w3*t2 + w1*t3 + w4*t4 + w2*t5; + x[k+i+4*j] = t1 + w4*t2 + w3*t3 + w2*t4 + w1*t5; + } + w *= u; + } + j = j*5; + } + +} + +void FFT2(complex *y, complex *x, int N){ + complex *Nx, *Ny; + int i; + complex wn, w; + + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(i = 0; i < N/2; i++){ + Nx[i] = x[2*i]; + Nx[i + N/2] = x[2*i + 1]; + } + + + if(N == 2){ + y[0] = x[0] + x[1]; + y[1] = x[0] - x[1]; + } + else{ + FFT2(Ny, Nx, N/2); + FFT2(Ny + N/2, Nx + N/2, N/2); + + wn = 1; w = cexp(-I*2*M_PI/N); + + for(i = 0; i < N/2; i++){ + y[i] = Ny[i] + wn*Ny[i + N/2]; + y[i + N/2] = Ny[i] - wn*Ny[i + N/2]; + wn *= w; + } + } + free(Nx); free(Ny); +} + +void FFT3(complex *y, complex *x, int N){ + complex *Nx, *Ny, w, wn, w1, w2; + complex wn1, wn2; + int i; + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(i = 0; i < N/3; i++){ + Nx[i] = x[3*i]; + Nx[i + N/3] = x[3*i + 1]; + Nx[i + 2*N/3] = x[3*i + 2]; + } + + w1 = cexp(-I*M_PI*2/3); // w4 == w1 + w2 = cexp(-I*M_PI*4/3); + if(N == 3){ + + y[0] = x[0] + x[1] + x[2]; + y[1] = x[0] + w1*x[1] + w2*x[2]; + y[2] = x[0] + w2*x[1] + w1*x[2]; + + } + else{ + FFT3(Ny, Nx, N/3); + FFT3(Ny + N/3, Nx + N/3, N/3); + FFT3(Ny + 2*N/3, Nx + 2*N/3, N/3); + + wn = 1; w = cexp(-I*2*M_PI/N); + for(i = 0; i < N/3; i++){ + + wn1 = wn*Ny[i + N/3]; + wn2 = cpow(wn, 2)*Ny[i + 2*N/3]; + + y[i] = Ny[i] + wn1 + wn2; + y[i + N/3] = Ny[i] + w1*wn1 + w2*wn2; + y[i + 2*N/3] = Ny[i] + w2*wn1 + w1*wn2; + wn *= w; + } + } + free(Nx); free(Ny); + + +} + +void FFT5(complex *y, complex *x, int N){ + complex *Nx, *Ny, w, wn, w1, w2, w3, w4; + complex wn1, wn2, wn3, wn4; + int i; + + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(i = 0; i < N/5; i++){ + Nx[i] = x[5*i]; + Nx[i + N/5] = x[5*i + 1]; + Nx[i + 2*N/5] = x[5*i + 2]; + Nx[i + 3*N/5] = x[5*i + 3]; + Nx[i + 4*N/5] = x[5*i + 4]; + } + + w1 = cexp(-I*M_PI*2/5); // w16 == w6 == w1 + w2 = cexp(-I*M_PI*4/5); // w12 == w2 + w3 = cexp(-I*M_PI*6/5); // w8 == w3 + w4 = cexp(-I*M_PI*8/5); // w9 == w4 + + if(N == 5){ + + y[0] = x[0] + x[1] + x[2] + x[3] + x[4]; + y[1] = x[0] + w1*x[1] + w2*x[2] + w3*x[3] + w4*x[4]; + y[2] = x[0] + w2*x[1] + w4*x[2] + w1*x[3] + w3*x[4]; + y[3] = x[0] + w3*x[1] + w1*x[2] + w4*x[3] + w2*x[4]; + y[4] = x[0] + w4*x[1] + w3*x[2] + w2*x[3] + w1*x[4]; + + } + else{ + FFT5(Ny, Nx, N/5); + FFT5(Ny + N/5, Nx + N/5, N/5); + FFT5(Ny + 2*N/5, Nx + 2*N/5, N/5); + FFT5(Ny + 3*N/5, Nx + 3*N/5, N/5); + FFT5(Ny + 4*N/5, Nx + 4*N/5, N/5); + + wn = 1; w = cexp(-I*2*M_PI/N); + + + for(i = 0; i < N/5; i++){ + wn1 = wn*Ny[i + N/5]; + wn2 = cpow(wn, 2)*Ny[i + 2*N/5]; + wn3 = cpow(wn, 3)*Ny[i + 3*N/5]; + wn4 = cpow(wn, 4)*Ny[i + 4*N/5]; + + y[i] = Ny[i] + wn1 + wn2 + wn3 + wn4; + y[i + N/5] = Ny[i] + w1*wn1 + w2*wn2 + w3*wn3 + w4*wn4; + y[i + 2*N/5] = Ny[i] + w2*wn1 + w4*wn2 + w1*wn3 + w3*wn4; + y[i + 3*N/5] = Ny[i] + w3*wn1 + w1*wn2 + w4*wn3 + w2*wn4; + y[i + 4*N/5] = Ny[i] + w4*wn1 + w3*wn2 + w2*wn3 + w1*wn4; + wn *= w; + } + } + free(Nx); free(Ny); + + +} + + +void FFTp(complex *y, complex *x, int N){ + complex *Nx, *Ny; + int i; + complex w, w1, w2, w3, w4, wn, wn1, wn2, wn3, wn4; + + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + if(N % 2 == 0){ + for(i = 0; i < N/2; i++){ + Nx[i] = x[2*i]; + Nx[i + N/2] = x[2*i + 1]; + } + + + if(N == 2){ + y[0] = x[0] + x[1]; + y[1] = x[0] - x[1]; + } + else{ + FFTp(Ny, Nx, N/2); + FFTp(Ny + N/2, Nx + N/2, N/2); + + wn = 1; w = cexp(-I*2*M_PI/N); + + for(i = 0; i < N/2; i++){ + y[i] = Ny[i] + wn*Ny[i + N/2]; + y[i + N/2] = Ny[i] - wn*Ny[i + N/2]; + wn *= w; + } + } + } + if(N % 3 == 0){ + for(i = 0; i < N/3; i++){ + Nx[i] = x[3*i]; + Nx[i + N/3] = x[3*i + 1]; + Nx[i + 2*N/3] = x[3*i + 2]; + } + + w1 = cexp(-I*M_PI*2/3); // w4 == w1 + w2 = cexp(-I*M_PI*4/3); + if(N == 3){ + + y[0] = x[0] + x[1] + x[2]; + y[1] = x[0] + w1*x[1] + w2*x[2]; + y[2] = x[0] + w2*x[1] + w1*x[2]; + + } + else{ + FFTp(Ny, Nx, N/3); + FFTp(Ny + N/3, Nx + N/3, N/3); + FFTp(Ny + 2*N/3, Nx + 2*N/3, N/3); + + wn = 1; w = cexp(-I*2*M_PI/N); + for(i = 0; i < N/3; i++){ + wn1 = wn*Ny[i + N/3]; + wn2 = cpow(wn, 2)*Ny[i + 2*N/3]; + y[i] = Ny[i] + wn1 + wn2; + y[i + N/3] = Ny[i] + w1*wn1 + w2*wn2; + y[i + 2*N/3] = Ny[i] + w2*wn1 + w1*wn2; + wn *= w; + } + } + } + if(N % 5 == 0){ + for(i = 0; i < N/5; i++){ + Nx[i] = x[5*i]; + Nx[i + N/5] = x[5*i + 1]; + Nx[i + 2*N/5] = x[5*i + 2]; + Nx[i + 3*N/5] = x[5*i + 3]; + Nx[i + 4*N/5] = x[5*i + 4]; + } + + w1 = cexp(-I*M_PI*2/5); // w16 == w6 == w1 + w2 = cexp(-I*M_PI*4/5); // w12 == w2 + w3 = cexp(-I*M_PI*6/5); // w8 == w3 + w4 = cexp(-I*M_PI*8/5); // w9 == w4 + + if(N == 5){ + + y[0] = x[0] + x[1] + x[2] + x[3] + x[4]; + y[1] = x[0] + w1*x[1] + w2*x[2] + w3*x[3] + w4*x[4]; + y[2] = x[0] + w2*x[1] + w4*x[2] + w1*x[3] + w3*x[4]; + y[3] = x[0] + w3*x[1] + w1*x[2] + w4*x[3] + w2*x[4]; + y[4] = x[0] + w4*x[1] + w3*x[2] + w2*x[3] + w1*x[4]; + + } + else{ + FFTp(Ny, Nx, N/5); + FFTp(Ny + N/5, Nx + N/5, N/5); + FFTp(Ny + 2*N/5, Nx + 2*N/5, N/5); + FFTp(Ny + 3*N/5, Nx + 3*N/5, N/5); + FFTp(Ny + 4*N/5, Nx + 4*N/5, N/5); + + wn = 1; w = cexp(-I*2*M_PI/N); + + + for(i = 0; i < N/5; i++){ + wn1 = wn*Ny[i + N/5]; + wn2 = cpow(wn, 2)*Ny[i + 2*N/5]; + wn3 = cpow(wn, 3)*Ny[i + 3*N/5]; + wn4 = cpow(wn, 4)*Ny[i + 4*N/5]; + + y[i] = Ny[i] + wn1 + wn2 + wn3 + wn4; + y[i + N/5] = Ny[i] + w1*wn1 + w2*wn2 + w3*wn3 + w4*wn4; + y[i + 2*N/5] = Ny[i] + w2*wn1 + w4*wn2 + w1*wn3 + w3*wn4; + y[i + 3*N/5] = Ny[i] + w3*wn1 + w1*wn2 + w4*wn3 + w2*wn4; + y[i + 4*N/5] = Ny[i] + w4*wn1 + w3*wn2 + w2*wn3 + w1*wn4; + wn *= w; + } + } + } +} + + +void FFT_p(complex *y, complex *x, int N){ + int p; + if(N%2 == 0) p = 2; + else if(N%3 == 0) p = 3; + else if(N%5 == 0) p = 5; + else p = N; + + + complex *Nx, *Ny; + int i, j, k, bias, idx; + + bias = N/p; //bias = N/p : °¾²¾¶q + + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(j = 0; j < p; j++){ + for(k = 0; k < bias; k++){ + Nx[k + j*bias] = x[p*k + j]; + } + } + + // F[p][p] + + complex F[p][p]; + for(i = 0; i < p; i++){ + F[0][i] = 1; + F[i][0] = 1; + } + for(i = 1; i < p/2+1; i++){ + for(j = 1; j < p; j++){ + F[i][j] = cexp(-2*M_PI*i*j*I/p); + } + } + for(i = p/2+1; i < p; i++){ //¦³¤Ï¹ïºÙ + for(j = 1; j < p; j++){ + F[i][j] = F[p-i][p-j]; + } + } + + + //recursive + if(N == p){ + for(i = 0; i < p; i++){ + y[i] = 0; + for(j = 0; j < p; j++){ + y[i] += F[i][j]*x[j]; + } + } + } + else{ + + for(j = 0; j < p; j++){ + FFT_p(Ny + j*bias, Nx + j*bias, bias); + } + + for(i = 0; i < bias; i++){ + for(k = 0; k < p; k++){ + Ny[i + k*bias] *= cexp(-2*M_PI*k*i*I/N); + } + + for(j = 0; j < p; j++){ + idx = i + j*bias; + y[idx] = 0; + for(k = 0; k < p; k++){ + y[idx] += F[j][k]*Ny[i + k*bias]; + } + + } + } + + } + + free(Nx); free(Ny); +} + + + +int main(){ + complex *x, *y; + int i , N, temp; + time_t T; +// N = 30; + N = pow(2, 26)*pow(3, 0)*pow(5, 0); + x = (complex*) malloc(N*sizeof(complex)); + y = (complex*) malloc(N*sizeof(complex)); + + printf("N = %d\n", N); + for(i = 0; i < N; i++) x[i] = i; + + T = clock(); + FFT2(y, x, N); + T = clock() - T; + printf("FFT2 time cost ¡G%d ms\n", T); + + + for(i = 0; i < N; i++) x[i] = i; + T = clock(); + FFT_2(x, N); + T = clock() - T; + printf("FFT_2 time cost ¡G%d ms\n", T); + +// printf("after FFT\n"); +// for(i = 0; i < N; i++) printf("%f + %f i\n", creal(y[i]), cimag(y[i])); + + + return 0; +} + diff --git a/MGsolver.c b/MGsolver.c new file mode 100644 index 0000000..c4ea714 --- /dev/null +++ b/MGsolver.c @@ -0,0 +1,278 @@ +#include +#include +#include +#include +#include +#include + + +int Initialize(double *u, int N){ + int i; + for(i = 0; i < N; ++i) u[i] = 0.0; + + return 0; +} + +int Evolution(double *u, double *un, int N){ + int i; + for(i = 0; i < N; ++i) u[i] += un[i]; + return 0; +} + +double Residual(double *r, double *u, double *f, int N){ + //compute r = f-Au + int i; + double v = 0.0; + int m = (int)sqrt(N); + + + for(i = 1; i < N; i++){ + if(i < m){ + r[i] = f[i] - (4*u[i] - u[i-1] - u[i+1] - u[i+m]); + } + else if(i >= m*(m-1)){ + r[i] = f[i] - (4*u[i] - u[i+1] - u[i-1] - u[i-m]); + } + else{ + r[i] = f[i] - (4*u[i] - u[i-1] - u[i-m] - u[i+1] - u[i+m]); + } + if(fabs(r[i]) > v) v = fabs(r[i]); + } + + return v; +} + +int Smooth(double *un, double *u, double *f, double omega, int N){ + int i; + double v = 0.0; + int m = (int)sqrt(N); + + + for(i = 1; i < N; i++){ + if(i < m){ + un[i] = u[i] + omega*(0.5)*(4*u[i] - u[i-1] - u[i+1] - u[i+m] - f[i]); + } + else if(i >= m*(m-1)){ + un[i] = u[i] - omega*(0.5)*(4*u[i] - u[i+1] - u[i-1] - u[i-m] - f[i]); + } + else{ + un[i] = u[i] - omega*(0.5)*(4*u[i] - u[i-1] - u[i-m] - u[i+1] - u[i+m] - f[i]); + } + } + + for(i = 1; i < N; i++){ + u[i] = un[i]; + } + + return 0; +} + +int Restriction(double *r_f, double *r_c, int Nc){ + // go to coarse grid + //r_c:coarse + //r_f:fine + + int i, j; + int m = ((int)sqrt(Nc))/2; + + j = 0; + for(i = m+1; i < Nc; ++i){ + //¸õ¹L²Ä¤@¦C & ³Ì«á¤@¦C + if(i%m == m-1){ + //¸õ¹L¬Û¾Fªº¤@¦C + i = i+m+1; + continue; + } + if(i&m == 0){ + i = i+m; + continue; + } + //i ¬°¤¤¤ßªº8­Ó¤è¦V³£­nºâ ,µM«á°£16 + r_c[j] = 0.0625*(4*r_f[i] + 2*(r_f[i-m] + r_f[i-1] + r_f[i+1] + r_f[i+m]) + r_f[i-m-1] + r_f[i-m+1] + r_f[i+m-1] + r_f[i+m+1]); + ++j; + ++i; + } + + return 0; +} + +int Interprolation(double *u_f, double *u_c, int Nf){ + //go to fine grid + //r_c:coarse + //r_f:fine + + + int i, j; + int mf = (int)sqrt(Nf); + int mc = mf/2; + + if(mf%2 == 0){ + //°¸¼Æ + int row = 1; + + for(i = 0; i < mc; ++i){ + for(j = 0; j < mc; ++j){ + int InFineIdx = mf*row + 2*j + 1; + int InCoaresIdx = mc*i+j; + + u_f[InFineIdx] = u_c[InCoaresIdx];//¤¤¤ß + + if(j != mc-1){ + u_f[InFineIdx+1] += u_c[InCoaresIdx]/2;//¥k + u_f[InFineIdx+1-mf] += u_c[InCoaresIdx]/4;//¥k¤W + u_f[InFineIdx+1+mf] += u_c[InCoaresIdx]/4;//¥k¤U + } + + if(i != mc-1){ + u_f[InFineIdx+mf] += u_c[InCoaresIdx]/2;//¤U + u_f[InFineIdx+1+mf] += u_c[InCoaresIdx]/4;//¥ª¤U + u_f[InFineIdx-1+mf] += u_c[InCoaresIdx]/4;//¥k¤U + } + + u_f[InFineIdx-mf] += u_c[InCoaresIdx]/2;//¤W + u_f[InFineIdx-1] += u_c[InCoaresIdx]/2;//¥ª + u_f[InFineIdx-1-mf] += u_c[InCoaresIdx]/4;//¥ª¤W + + if(i == mc-1 && j == mc-1){ + u_f[InFineIdx-1+mf] -= u_c[InCoaresIdx]/4;//¥k¤U³Q¦hºâ¤@¦¸ + } + } + row += 2; + } + + } + else{ + //©_¼Æ + int row = 1; + for(i = 0; i < mc; ++i){ + for(j = 0; j < mc; ++j){ + int InFineIdx = mf*row + 2*j + 1; + int InCoaresIdx = mc*i+j; + + u_f[InFineIdx-mf-1] = u_c[InCoaresIdx]/4;//¥ª¤W + u_f[InFineIdx-mf] = u_c[InCoaresIdx]/2;//¤W + u_f[InFineIdx-mf+1] = u_c[InCoaresIdx]/4;//¥k¤W + u_f[InFineIdx-1] = u_c[InCoaresIdx]/2;//¥ª + u_f[InFineIdx] = u_c[InCoaresIdx];//¤¤¤ß + u_f[InFineIdx+1] = u_c[InCoaresIdx]/2;//¥k + u_f[InFineIdx+mf-1] = u_c[InCoaresIdx]/4;//¥ª¤U + u_f[InFineIdx+mf] = u_c[InCoaresIdx]/2;//¤U + u_f[InFineIdx+mf+1] = u_c[InCoaresIdx]/4;//¥k¤U + } + } + } + + return 0; +} + +void Source(double **F, int N){ + // make F + + int i, j, k; + double x, y, h; + h = 1.0/N; + + for(i = 0; i < N-1; ++i){ + for(j = 0; j < N-1; ++j){ + x = (i+1)*h; + y = (j+1)*h; + F[i][j] = -(1.0+4.0)*h*h*M_PI*M_PI*sin(M_PI*x)*sin(2*M_PI*y); + } + } +} + + + +int main(){ + int i, j, k, s, *Level, N, M; + double **MEM, **f, **u, **un, **r; + double res, res0, w, h; + + // solving Au = f + N = 4; + M = (N-1)*(N-1); + + int levelSize = 0; + int tempN = N; + while(tempN != 0){ + ++levelSize; + tempN /= 2; + } +// printf("size is %d\n", levelSize); + //open memory + Level = (int*)malloc((levelSize+1)*sizeof(int)); + MEM = (double**)malloc(5*levelSize*sizeof(double*)); + u = MEM + levelSize; + un = u + levelSize; + f = un + levelSize; + r = f + levelSize; + + N = 2; + for(j = 1; j < levelSize; ++j){ + printf("%d\n", N); + // open memory once + MEM[j-1] = (double*)malloc(4*(N+1)*sizeof(double)); + + u[j] = MEM[j-1]; + un[j] = u[j] + N + 1; + f[j] = un[j] + N+1; + r[j] = f[j] + N+1; + + h = 1.0/N; + + + } + + + + + + // end of open memory + + + //make F + + + + +// for(i = 0; i < levelSize; i++) printf("%d\n", Level[i]); + + + //begin the multigrid +// res0 = Residual(r, u, f, Level[levelSize]); +// printf("%f\n", res0); + +// for(i = 0; i < M; i++) printf("%f %f %f\n",r[i], u[i], b[i]); + + +// w = 2.0/3; +// +// for(k = 0; k < 10; ++k){ +// for(j = levelSize; j > 0; --j){ +// Smooth(un, u, f, w, Level[j]); +// Residual(r, u, f, Level[j]); +// +// +// } +// } + + + + + + + + + + + + + + + + + + + + return 0; +} diff --git a/MGsolver2D.c b/MGsolver2D.c new file mode 100644 index 0000000..e3e58f2 --- /dev/null +++ b/MGsolver2D.c @@ -0,0 +1,201 @@ +#include +#include +#include +#include + +double Residual(double **r, double **u, double **f, int N); + +int main(){ + int i, j, k, s, L = 5, *M, N, N2; + double ***Mem, ***u, ***un, ***f, ***r, res, res0, h, w; + time_t t1, t2; + + // [ 0, -1, 0] + // solving Au = f, A = [-1, 4, -1] + // [ 0, -1, 0] + + // for memory allocation! + M = (int *) malloc((L+1)*sizeof(int)); + Mem = (double ***) malloc(5*L*sizeof(double*)); + u = Mem + L - 1; // -1 makes u[0~L-1] to u[1~L] + un = u + L; + f = un + L; + r = f + L; + + N = 2; + for(k = 1; k <= L; ++k) { + printf("%d\n", N); // N = 2^(j+1), j=1..L + N2 = (N+1)*(N+1); + Mem[k-1] = (double **) malloc(4*(N+1)*sizeof(double*)); // open memory once + Mem[k-1][0] = (double*) malloc(4*N2*sizeof(double)); + + u[k] = Mem[k-1]; + un[k] = u[k] + (N+1); + f[k] = un[k] + (N+1); + r[k] = f[k] + (N+1); + + u[k][0] = Mem[k-1][0]; + un[k][0] = u[k][0] + N2; + f[k][0] = un[k][0] + N2; + r[k][0] = f[k][0] + N2; + + for(i = 1; i <= N; ++i){ + u[k][i] = u[k][i-1] + N + 1; + un[k][i] = un[k][i-1] + N + 1; + f[k][i] = f[k][i-1] + N + 1; + r[k][i] = r[k][i-1] + N + 1; + } + + h = 1.0/N; + + for(i = 0; i <= N; ++i){ + for(j = 0; j <= N; ++j){ + u[k][i][j] = 0.0; + un[k][i][j] = 0.0; + r[k][i][j] = 0.0; + f[k][i][j] = -1.0*h*h; + } + + } + M[k] = N; + N = N*2; + } + k = L; + res0 = Residual(r[k],u[k],f[k],M[k]); + w = 0.9; + for(k = 0 ; k < 10; ++k){ + for(j = L; j > 1; --j){ + SmootherG(un[j],u[j],f[j],w,M[j]); + Residual(r[j],u[j],f[j],M[j]); + Restriction(r[j], f[j-1], M[j-1]); + Initial(u[j-1],M[j-1]); + } + SmootherJ(un[j],u[j],f[j],1.0,M[j]); + for(j = 2; j <= L; ++j){ + Interpolation(un[j], u[j-1], M[j]); + Evolution(u[j],un[j],M[j]); + SmootherG(un[j],u[j],f[j],w,M[j]); + } + j = L; + res = Residual(r[j],u[j],f[j],M[j]); + printf("%d Res: %e Res_reduce: %.3f\n",k,res,res/res0); + res0 = res; + } + + + for(i=0;imaxres) maxres = fabs(r[i][j]); + } + + + } + return maxres; +} + +int SmootherJ(double **un, double **u, double **f, double omega, int N){ + int i, j; + + for(i = 1; i < N; ++i){ + for(j = 1; j < N; ++j){ + un[i][j] = u[i][j] + omega*(0.25)*( (u[i+1][j] - 4.0*u[i][j] + u[i-1][j] + u[i][j+1] + u[i][j-1]) - f[i][j] ); + + } + } + for(i = 1; i < N; ++i){ + for(j = 1; j < N; ++j){ + u[i][j] = un[i][j]; + } + } + + return 0; +} + +int SmootherG(double **un, double **u, double **f, double omega, int N){ + int i, j; + omega = 1.0; + + for(i = 1; i < N; ++i){ + for(j = 1; j < N; ++j){ + u[i][j] = u[i][j] + omega*(0.25)*( (u[i+1][j] - 4.0*u[i][j] + u[i-1][j] + u[i][j+1] + u[i][j-1]) - f[i][j] ); + + } + } + + return 0; +} + +int Restriction(double **r_f, double **r_c, int Nc) { + // + int i, j; + for(i = 1; i < Nc; ++i){ + for(j = 1; j < Nc; ++j){ + r_c[i][j] = ( r_f[2*i-1][2*j-1] + 2.0*r_f[2*i][2*j-1] + r_f[2*i+1][2*j-1]+ + 2.0*r_f[2*i-1][2*j] + 4.0*r_f[2*i][2*j] + 2.0*r_f[2*i+1][2*j]+ + r_f[2*i-1][2*j+1] + 2.0*r_f[2*i][2*j+1] + r_f[2*i+1][2*j+1] )/4.0; + + } + + } + return 0; +} + +int Interpolation(double **u_f, double **u_c, int Nf) { + // + int i, j; + + for(i = 1; i < Nf; ++i){ + for(j = 1; j < Nf; ++j){ + + if(i%2 == 0 && j%2 == 0){ + u_f[i][j] = u_c[i/2][j/2]; + } + else if(i%2 == 1 && j%2 == 0){ + u_f[i][j] = 0.5*(u_c[i/2][j/2] + u_c[i/2+1][j/2]); + } + else if(i%2 == 0 && j%2 == 1){ + u_f[i][j] = 0.5*(u_c[i/2][j/2] + u_c[i/2][j/2+1]); + } + else{ + u_f[i][j] = 0.25*( u_c[i/2][j/2] + u_c[i/2][j/2+1] + + u_c[i/2+1][j/2] + u_c[i/2+1][j/2+1]); + } + } + } + + return 0; +} + + diff --git a/Plot.dat b/Plot.dat new file mode 100644 index 0000000..e69de29 diff --git a/base_p.c b/base_p.c new file mode 100644 index 0000000..bbb2721 --- /dev/null +++ b/base_p.c @@ -0,0 +1,65 @@ +#include +#include + + +void PrintBaseP(int a, int p){ // input 1­Ó¾ã¼ÆÂରp¶i¦ì¡A¨Ã¦L¥X¦¹p¶i¦ì¤Î¤Ï§Ç + if(a < p) printf("k = %d\nj = %d\n", a, a); + else{ + int i, cnt, temp, *output, size; + temp = a; + cnt = 0; + while(temp > p-1){ + cnt++; + temp /= p; + } + cnt++; + size = cnt; + printf("cnt = %d\n", cnt); + output = (int*)malloc(cnt*sizeof(int));//0 ~ cnt + + + while(cnt > 0){ + printf("put in : %d\n", a%p); + output[--cnt] = a%p; + a /= p; + } + printf("size = %d\n",sizeof(output)); + + printf("k = "); + for(i = 0; i < size; i++){ + printf("%d", output[i]); + } + printf("\n"); + + printf("j = "); + for(i = size-1; i >= 0; i--){ + printf("%d", output[i]); + } + printf("\n"); + + + + } + + + + + + + +} + + + + + + +int main(){ + base_p_1(100, 4); + + + + + + return 0; +} diff --git a/inverseFFT.c b/inverseFFT.c new file mode 100644 index 0000000..f5ad9c8 --- /dev/null +++ b/inverseFFT.c @@ -0,0 +1,319 @@ +#include +#include +#include +#include +#include +#include + + +void FFT2(complex *y, complex *x, int N){ + complex *Nx, *Ny; + int i; + complex wn, w; + + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(i = 0; i < N/2; i++){ + Nx[i] = x[2*i]; + Nx[i + N/2] = x[2*i + 1]; + } + + + if(N == 2){ + y[0] = x[0] + x[1]; + y[1] = x[0] - x[1]; + } + else{ + FFT2(Ny, Nx, N/2); + FFT2(Ny + N/2, Nx + N/2, N/2); + + wn = 1; w = cexp(-I*2*M_PI/N); + + for(i = 0; i < N/2; i++){ + y[i] = Ny[i] + wn*Ny[i + N/2]; + y[i + N/2] = Ny[i] - wn*Ny[i + N/2]; + wn *= w; + } + } + free(Nx); free(Ny); +} + +void iFFT2(complex *y, complex *x, int N){ + complex *Nx, *Ny; + int i; + complex wn, w; + + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(i = 0; i < N/2; i++){ + Nx[i] = x[2*i]; + Nx[i + N/2] = x[2*i + 1]; + } + + + if(N == 2){ + y[0] = x[0] + x[1]; + y[1] = x[0] - x[1]; + } + else{ + iFFT2(Ny, Nx, N/2); + iFFT2(Ny + N/2, Nx + N/2, N/2); + + wn = 1; w = cexp(I*2*M_PI/N); + + for(i = 0; i < N/2; i++){ + y[i] = Ny[i] + wn*Ny[i + N/2]; + y[i + N/2] = Ny[i] - wn*Ny[i + N/2]; + wn *= w; + } + } + + free(Nx); free(Ny); +} + + +void FFT_p(complex *y, complex *x, int N){ + int p; + if(N%2 == 0) p = 2; + else if(N%3 == 0) p = 3; + else if(N%5 == 0) p = 5; + else p = N; + + + complex *Nx, *Ny; + int i, j, k, bias, idx; + + bias = N/p; //bias = N/p : °¾²¾¶q + + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(j = 0; j < p; j++){ + for(k = 0; k < bias; k++){ + Nx[k + j*bias] = x[p*k + j]; + } + } + + // F[p][p] + + complex F[p][p]; + for(i = 0; i < p; i++){ + F[0][i] = 1; + F[i][0] = 1; + } + for(i = 1; i < p/2+1; i++){ + for(j = 1; j < p; j++){ + F[i][j] = cexp(-2*M_PI*i*j*I/p); + } + } + for(i = p/2+1; i < p; i++){ //¦³¤Ï¹ïºÙ + for(j = 1; j < p; j++){ + F[i][j] = F[p-i][p-j]; + } + } + + + //recursive + if(N == p){ + for(i = 0; i < p; i++){ + y[i] = 0; + for(j = 0; j < p; j++){ + y[i] += F[i][j]*x[j]; + } + } + } + else{ + + for(j = 0; j < p; j++){ + FFT_p(Ny + j*bias, Nx + j*bias, bias); + } + + for(i = 0; i < bias; i++){ + for(k = 0; k < p; k++){ + Ny[i + k*bias] *= cexp(-2*M_PI*k*i*I/N); + } + + for(j = 0; j < p; j++){ + idx = i + j*bias; + y[idx] = 0; + for(k = 0; k < p; k++){ + y[idx] += F[j][k]*Ny[i + k*bias]; + } + + } + } + + } + + free(Nx); free(Ny); +} + + + +void iFFT_p(complex *y, complex *x, int N){ + int p; + if(N%2 == 0) p = 2; + else if(N%3 == 0) p = 3; + else if(N%5 == 0) p = 5; + else p = N; + + + complex *Nx, *Ny; + int i, j, k, bias, idx; + + bias = N/p; //bias = N/p : °¾²¾¶q + + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(j = 0; j < p; j++){ + for(k = 0; k < bias; k++){ + Nx[k + j*bias] = x[p*k + j]; + } + } + + // F[p][p] + + complex F[p][p]; + for(i = 0; i < p; i++){ + F[0][i] = 1; + F[i][0] = 1; + } + for(i = 1; i < p/2+1; i++){ + for(j = 1; j < p; j++){ + F[i][j] = cexp(2*M_PI*i*j*I/p); + } + } + for(i = p/2+1; i < p; i++){ //¦³¤Ï¹ïºÙ + for(j = 1; j < p; j++){ + F[i][j] = F[p-i][p-j]; + } + } + + + //recursive + if(N == p){ + for(i = 0; i < p; i++){ + y[i] = 0; + for(j = 0; j < p; j++){ + y[i] += F[i][j]*x[j]; + } + } + } + else{ + + for(j = 0; j < p; j++){ + iFFT_p(Ny + j*bias, Nx + j*bias, bias); + } + + for(i = 0; i < bias; i++){ + for(k = 0; k < p; k++){ + Ny[i + k*bias] *= cexp(2*M_PI*k*i*I/N); + } + + for(j = 0; j < p; j++){ + idx = i + j*bias; + y[idx] = 0; + for(k = 0; k < p; k++){ + y[idx] += F[j][k]*Ny[i + k*bias]; + } + + } + } + + } + + free(Nx); free(Ny); +} + + + + + + +int main(){ + complex *x, *y, *u, *v, *ans_t; + int i , N; + + N = 3; + + + x = (complex*) malloc(2*N*sizeof(complex)); + y = (complex*) malloc(2*N*sizeof(complex)); + u = (complex*) malloc(2*N*sizeof(complex)); + v = (complex*) malloc(2*N*sizeof(complex)); + ans_t = (complex*) malloc(2*N*sizeof(complex)); + + char temp[N], a; + printf("enter the first number x : \n"); + fgets(temp, N*2, stdin); + for(i = 0; i < N; i++){ + a = temp[i]; + x[N-1-i] = atoi(&a); + } + for(i = N; i < 2*N; i++) x[i] = 0; + + printf("\n"); + printf("enter the second number y : \n"); + fgets(temp, N*2, stdin); + for(i = 0; i < N; i++){ + a = temp[i]; + y[N-1-i] = atoi(&a); + } + for(i = N; i < 2*N; i++) y[i] = 0; + +// for(i = 0; i < 2*N; i++){ +// printf("%f + %f i\n", creal(x[i]), cimag(x[i])); +// } +// printf("\n"); +// +// for(i = 0; i < 2*N; i++){ +// printf("%f + %f i\n", creal(y[i]), cimag(y[i])); +// } + + //u: x after fft + //v: y after fft + + FFT_p(u, x, 2*N); + FFT_p(v, y, 2*N); + +// for(i = 0; i < 2*N; i++){ +// printf("%f + %f i\n", creal(u[i]), cimag(u[i])); +// } +// printf("\n"); +// +// for(i = 0; i < 2*N; i++){ +// printf("%f + %f i\n", creal(v[i]), cimag(v[i])); +// } +// printf("\n"); + + for(i = 0; i < 2*N; i++){ + u[i]*=v[i]; + ans_t[i] = 0; + } + +// for(i = 0; i < 2*N; i++){ +// printf("%f + %f i\n", creal(u[i]), cimag(u[i])); +// } +// printf("\n"); + + iFFT_p(ans_t, u, 2*N); + for(i = 0; i < 2*N; i++) ans_t[i] /= (2*N); + + for(i = 0; i < 2*N; i++) printf("%f\n", creal(ans_t[i])); + + float ans; + ans = 0; + for(i = 0; i < 2*N; i++){ + ans += pow(10, i)*creal(ans_t[i]); + } + + + printf("ans is %f\n", ans); + + + return 0; +} + diff --git a/radix2.c b/radix2.c new file mode 100644 index 0000000..c045bab --- /dev/null +++ b/radix2.c @@ -0,0 +1,64 @@ +#include +#include +#include +#include +time_t timer; +int main() +{ + int j, k, p = 20, n, L, M, N = 1<=L) { + j -= L; + L /= 2; + } + j += L; + if(j>k) { + t = x[k]; + x[k] = x[j]; + x[j] = t; + } + } + + // Butterfly Radix + j = 1; + while(j multiply w + x[k+n+j] *= w; + // FFT2 for k+n and k+n+j + t = x[k+n]; + x[k+n ] = x[k+n]+x[k+n+j]; + x[k+n+j] = t-x[k+n+j]; + } + w *= u; + } + j = j << 1; + } + timer = clock() - timer; +// for(k=0;k // for printf function +#include // for memory allocation +#include // for time calculation +#include // for sine and cosine functions + + + +void quick_sort(double*, int, int); //quick sort +int Partition(double*, int, int); //sub algo of quick sort +double Select(double*, int, int, int); // find the i-th number +double Find_Med(double*, int, int); //find the median +double QD(double*, int, int); // find the Quartile Deviation +int main() { + // Declare all the variables + int k, m, n, N; + double *x, *y, z, p; + time_t t; + + // Input the number N + printf("Input N: "); + scanf("%d",&N); + + // Locate the memory for x and y; + x = (double *) malloc(N*sizeof(double)); + y = (double *) malloc(N*sizeof(double)); + + // Initial setting for x, for example, x[k] = 1.0*rand()/RAND_MAX + srand( time(NULL) ); + for(k=0;k x[k]) { + z = x[n]; + x[n] = x[k]; + x[k] = z; + } + } + } + t = clock() - t; + + // print y, x, and time + printf("Bubble sort %d elements: %f s\n", N, 1.0*t/CLOCKS_PER_SEC); + + printf("bubble sort : \n"); + for(k=0;k // for printf function -#include // for memory allocation -#include // for time calculation -#include // for sine and cosine functions - - - -void quick_sort(double*, int, int); //quick sort -int Partition(double*, int, int); //sub algo of quick sort -double Select(double*, int, int, int); // find the i-th number -double Find_Med(double*, int, int); //find the median -double QD(double*, int, int); // find the Quartile Deviation -int main() { - // Declare all the variables - int k, m, n, N; - double *x, *y, z, p; - time_t t; +#include +#include +#include +#include +#include +#include - // Input the number N - printf("Input N: "); - scanf("%d",&N); - - // Locate the memory for x and y; - x = (double *) malloc(N*sizeof(double)); - y = (double *) malloc(N*sizeof(double)); - - // Initial setting for x, for example, x[k] = 1.0*rand()/RAND_MAX - srand( time(NULL) ); - for(k=0;k x[k]) { - z = x[n]; - x[n] = x[k]; - x[k] = z; - } + +void Exact_Discretization(double **A, int N){ + int i,j,k; + for(i=0;i0) A[k][k-(N-1)] = 1; + if(i>0) A[k][k-1] = 1; + if(i= 0; --i){ + for(j = i+1; j < N; ++j){ + x[i] -= M[i][j]*x[j]; + } + x[i] /= M[i][i]; + } + // free memory + free(M[0]); + free(M); +} - - // free the memory located by x, y - free(x); - free(y); - - return 100; +double Residual(double *r, double **A, double *x, double *b, int N){ + //compute r = b-Ax + int i, j; + double v = 0.0; + for(i = 0; i < N; ++i){ + r[i] = b[i]; + + for(j = 0; j < N; ++j){ + r[i] -= A[i][j]*x[j]; + } + + if(fabs(r[i]) > v) v = fabs(r[i]); + } + return v; } -void quick_sort(double *x, int P, int r){ - double temp; - if(P < r){ - int q = Partition(x, P, r); - quick_sort(x, P, q - 1); - quick_sort(x, q + 1, r); +double error(double *x, double *u, int N){ + int i; + double e, v = 0.0; + + for(i = 0; i < N; ++i){ + e = fabs(x[i]-u[i]); + if(e > v) v = e; } + return v; } -int Partition(double *x, int P, int r){ - double a = x[r]; - int i = P - 1, k; - for(k = P; k < r; k++){ - if(x[k] <= a){ - i++; - double temp; - temp = x[i]; - x[i] = x[k]; - x[k] = temp; +void FFT_p(complex *y, complex *x, int N){ + int p; + if(N%2 == 0) p = 2; + else if(N%3 == 0) p = 3; + else if(N%5 == 0) p = 5; + else p = N; + + + complex *Nx, *Ny; + int i, j, k, bias, idx; + + bias = N/p; //bias = N/p : °¾²¾¶q + + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(j = 0; j < p; ++j){ + for(k = 0; k < bias; ++k){ + Nx[k + j*bias] = x[p*k + j]; } } - double temp; - temp = x[i + 1]; - x[i + 1] = x[r]; - x[r] = temp; - return (i + 1); -} + + // F[p][p] + +// complex F[p][p]; -double Select(double *x, int P, int r, int i){ - if(P < r){ - int q = Partition(x, P, r); - int k = q - P + 1; - if(i == k) return x[q]; - else if(i < k) return Select(x, P, q - 1, i); - else return Select(x, q + 1, r, i - k); + complex **F; + F = (complex**)malloc(p*sizeof(complex*)); + for(i = 0; i < p; ++i){ + F[i] = (complex*)malloc(p*sizeof(complex)); + } + + for(i = 0; i < p; ++i){ + F[0][i] = 1; + F[i][0] = 1; + } + for(i = 1; i < p/2+1; ++i){ + for(j = 1; j < p; ++j){ + F[i][j] = cexp(-2*M_PI*i*j*I/p); + } + } + for(i = p/2+1; i < p; ++i){ //¦³¤Ï¹ïºÙ + for(j = 1; j < p; ++j){ + F[i][j] = F[p-i][p-j]; + } + } + + + //recursive + if(N == p){ + for(i = 0; i < p; ++i){ + y[i] = 0; + for(j = 0; j < p; ++j){ + y[i] += F[i][j]*x[j]; + } + } } else{ - return x[r]; + + for(j = 0; j < p; ++j){ + FFT_p(Ny + j*bias, Nx + j*bias, bias); + } + + for(i = 0; i < bias; ++i){ + for(k = 0; k < p; ++k){ + Ny[i + k*bias] *= cexp(-2*M_PI*k*i*I/N); + } + + for(j = 0; j < p; ++j){ + idx = i + j*bias; + y[idx] = 0; + for(k = 0; k < p; ++k){ + y[idx] += F[j][k]*Ny[i + k*bias]; + } + + } + } + } -} + + free(Nx); free(Ny); +} -double Find_Med(double *x, int P, int r){ - int number = r - P + 1; - if(number % 2 == 0){ - double right = Select(x, P, r, number/2 + 1); - double left = Select(x, P, r, number/2); - return (right + left)/2; +void DST(double *y, double *x, int N){ + //y is the answer + + int i; + N = 2*N + 2; + complex *Nx, *Ny; + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(i = 0; i < N; ++i){ + Nx[i] = 0.0; } - else{ - return Select(x, P, r, number/2 + 1); + + for(i = 0; i < N/2-1; ++i){ + Nx[i+1] = x[i]; + Nx[(i+1) + N/2] = -x[N/2-(i+2)]; } + + + FFT_p(Ny, Nx, N); + + + for(i = 0; i < N/2-1; ++i){ + y[i] = -0.5*cimag(Ny[i+1]); + } + + free(Nx); free(Ny); + } -double QD(double *x, int P, int r){ // ¥|¤À¦ì®t Quartile Deviation, QD - int number = r - P + 1; - if(number/2 == 0){ - double Q1 = Find_Med(x, P, number/2 - 1); - printf("Q1 = %f\n", Q1); - double Q3 = Find_Med(x, number/2, r); - printf("Q3 = %f\n", Q3); - return Q3 - Q1; +void iDST(double *y, double *x, int N){ + //y is the answer + int i; + N = 2*N + 2; + complex *Nx, *Ny; + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(i = 0; i < N; ++i){ + Nx[i] = 0.0; } - else{ - double Q1 = Find_Med(x, P, number/2 - 1); - printf("Q1 = %f\n", Q1); - double Q3 = Find_Med(x, (number + 1)/2, r); - printf("Q3 = %f\n", Q3); - return Q3 - Q1; - } + + for(i = 0; i < N/2-1; ++i){ + Nx[i+1] = x[i]; + Nx[(i+1) + N/2] = -x[N/2-(i+2)]; + } + + + FFT_p(Ny, Nx, N); + + + for(i = 0; i < N/2-1; ++i){ + y[i] = (-2.0/N)*cimag(Ny[i+1]); + } + + free(Nx); free(Ny); } +void Transpose(double **A, int N){ + //transpose a NxN matrix + int i , j; + double temp; + + for(i = 0; i < N; ++i){ + for(j = i+1; j < N; ++j){ + temp = A[i][j]; + A[i][j] = A[j][i]; + A[j][i] = temp; + } + } +} +void DST2D(double **X, int N){ + int i; + + for(i = 0; i < N; ++i){ + DST(X[i], X[i], N); + } + + Transpose(X, N); + + for(i = 0; i < N; ++i){ + DST(X[i], X[i], N); + } + + Transpose(X, N); +} +void iDST2D(double **X, int N){ + int i; + + for(i = 0; i < N; ++i){ + iDST(X[i], X[i], N); + } + + Transpose(X, N); + + for(i = 0; i < N; ++i){ + iDST(X[i], X[i], N); + } + + Transpose(X, N); +} +void FastPoissonSolver(double **F, int N){ + int i, j; + double h; + + h = 1.0/(N+1); + DST2D(F, N); + + + + for(i = 0; i < N; ++i){ + for(j = 0; j < N; ++j){ + F[i][j]= F[i][j]/((2*(cos(M_PI*(i+1)*h)-1))+(2*(cos(M_PI*(j+1)*h)-1))); + } + } + + + iDST2D(F, N); + +} +int main(){ + //solve AU = F, where r = b - Ax is the residue + // U is unknown, F is known + int i, j, k, N, M; + double **A, *x, *u, **U, *b, **F, *r; + time_t t; + + N = 16; + + M = (N-1)*(N-1); + //create memory, ±N2ºûarray¶}¦¨¤@¦C + A = (double**) malloc(M*sizeof(double*)); + A[0] = (double*) malloc(M*M*sizeof(double)); + for(i = 1; i < M; ++i){ + A[i] = A[i-1] + M; + } + x = (double*) malloc(M*sizeof(double)); + r = (double*) malloc(M*sizeof(double)); + + + b = (double*) malloc(M*sizeof(double)); + F = (double**)malloc((N-1)*sizeof(double*)); + F[0] = b; + for(i = 1; i < N-1; ++i){ + F[i] = F[i-1] + N - 1; + } + + + u = (double*) malloc(M*sizeof(double)); + U = (double**) malloc((N-1)*sizeof(double*)); + U[0] = u; + for(i = 1; i < N-1; ++i){ + U[i] = U[i-1] + N - 1; + } + //end of open memory + + + //initialize + for(i = 0; i < M*M; ++i){ + A[0][i] = 0.0; + } + //end of initialize + + //make the matrix A + Exact_Discretization(A, N); + + //make F + Source(F, N); + + //case1. use Gauss elimination solve + t = clock(); + BandGaussianElimination(A, x, b, M); - - - - - - - - - + t = clock() - t; + printf("use Gauss elimination: %f sec\n", 1.0*t/CLOCKS_PER_SEC); + printf("the residual : %e\n", Residual(r, A, x, b, M)); + + //calculate the true answer in order to find the error + ExactSolution(U, N); + printf("error is: %e\n", error(x, u, M)); + + + //case 2. use DST fast poisson solver + t = clock(); + FastPoissonSolver(F, N-1); + t = clock() - t; + + printf("use DST fast poisson solver: %f sec\n", 1.0*t/CLOCKS_PER_SEC); + printf("error is: %e\n", error(b, u, M)); + + + //free memory + free(A[0]); + free(A); + free(x); + free(r); + free(b); + free(u); + free(U); + free(F); + + + return 0; +} diff --git a/test2.c b/test2.c new file mode 100644 index 0000000..e3e58f2 --- /dev/null +++ b/test2.c @@ -0,0 +1,201 @@ +#include +#include +#include +#include + +double Residual(double **r, double **u, double **f, int N); + +int main(){ + int i, j, k, s, L = 5, *M, N, N2; + double ***Mem, ***u, ***un, ***f, ***r, res, res0, h, w; + time_t t1, t2; + + // [ 0, -1, 0] + // solving Au = f, A = [-1, 4, -1] + // [ 0, -1, 0] + + // for memory allocation! + M = (int *) malloc((L+1)*sizeof(int)); + Mem = (double ***) malloc(5*L*sizeof(double*)); + u = Mem + L - 1; // -1 makes u[0~L-1] to u[1~L] + un = u + L; + f = un + L; + r = f + L; + + N = 2; + for(k = 1; k <= L; ++k) { + printf("%d\n", N); // N = 2^(j+1), j=1..L + N2 = (N+1)*(N+1); + Mem[k-1] = (double **) malloc(4*(N+1)*sizeof(double*)); // open memory once + Mem[k-1][0] = (double*) malloc(4*N2*sizeof(double)); + + u[k] = Mem[k-1]; + un[k] = u[k] + (N+1); + f[k] = un[k] + (N+1); + r[k] = f[k] + (N+1); + + u[k][0] = Mem[k-1][0]; + un[k][0] = u[k][0] + N2; + f[k][0] = un[k][0] + N2; + r[k][0] = f[k][0] + N2; + + for(i = 1; i <= N; ++i){ + u[k][i] = u[k][i-1] + N + 1; + un[k][i] = un[k][i-1] + N + 1; + f[k][i] = f[k][i-1] + N + 1; + r[k][i] = r[k][i-1] + N + 1; + } + + h = 1.0/N; + + for(i = 0; i <= N; ++i){ + for(j = 0; j <= N; ++j){ + u[k][i][j] = 0.0; + un[k][i][j] = 0.0; + r[k][i][j] = 0.0; + f[k][i][j] = -1.0*h*h; + } + + } + M[k] = N; + N = N*2; + } + k = L; + res0 = Residual(r[k],u[k],f[k],M[k]); + w = 0.9; + for(k = 0 ; k < 10; ++k){ + for(j = L; j > 1; --j){ + SmootherG(un[j],u[j],f[j],w,M[j]); + Residual(r[j],u[j],f[j],M[j]); + Restriction(r[j], f[j-1], M[j-1]); + Initial(u[j-1],M[j-1]); + } + SmootherJ(un[j],u[j],f[j],1.0,M[j]); + for(j = 2; j <= L; ++j){ + Interpolation(un[j], u[j-1], M[j]); + Evolution(u[j],un[j],M[j]); + SmootherG(un[j],u[j],f[j],w,M[j]); + } + j = L; + res = Residual(r[j],u[j],f[j],M[j]); + printf("%d Res: %e Res_reduce: %.3f\n",k,res,res/res0); + res0 = res; + } + + + for(i=0;imaxres) maxres = fabs(r[i][j]); + } + + + } + return maxres; +} + +int SmootherJ(double **un, double **u, double **f, double omega, int N){ + int i, j; + + for(i = 1; i < N; ++i){ + for(j = 1; j < N; ++j){ + un[i][j] = u[i][j] + omega*(0.25)*( (u[i+1][j] - 4.0*u[i][j] + u[i-1][j] + u[i][j+1] + u[i][j-1]) - f[i][j] ); + + } + } + for(i = 1; i < N; ++i){ + for(j = 1; j < N; ++j){ + u[i][j] = un[i][j]; + } + } + + return 0; +} + +int SmootherG(double **un, double **u, double **f, double omega, int N){ + int i, j; + omega = 1.0; + + for(i = 1; i < N; ++i){ + for(j = 1; j < N; ++j){ + u[i][j] = u[i][j] + omega*(0.25)*( (u[i+1][j] - 4.0*u[i][j] + u[i-1][j] + u[i][j+1] + u[i][j-1]) - f[i][j] ); + + } + } + + return 0; +} + +int Restriction(double **r_f, double **r_c, int Nc) { + // + int i, j; + for(i = 1; i < Nc; ++i){ + for(j = 1; j < Nc; ++j){ + r_c[i][j] = ( r_f[2*i-1][2*j-1] + 2.0*r_f[2*i][2*j-1] + r_f[2*i+1][2*j-1]+ + 2.0*r_f[2*i-1][2*j] + 4.0*r_f[2*i][2*j] + 2.0*r_f[2*i+1][2*j]+ + r_f[2*i-1][2*j+1] + 2.0*r_f[2*i][2*j+1] + r_f[2*i+1][2*j+1] )/4.0; + + } + + } + return 0; +} + +int Interpolation(double **u_f, double **u_c, int Nf) { + // + int i, j; + + for(i = 1; i < Nf; ++i){ + for(j = 1; j < Nf; ++j){ + + if(i%2 == 0 && j%2 == 0){ + u_f[i][j] = u_c[i/2][j/2]; + } + else if(i%2 == 1 && j%2 == 0){ + u_f[i][j] = 0.5*(u_c[i/2][j/2] + u_c[i/2+1][j/2]); + } + else if(i%2 == 0 && j%2 == 1){ + u_f[i][j] = 0.5*(u_c[i/2][j/2] + u_c[i/2][j/2+1]); + } + else{ + u_f[i][j] = 0.25*( u_c[i/2][j/2] + u_c[i/2][j/2+1] + + u_c[i/2+1][j/2] + u_c[i/2+1][j/2+1]); + } + } + } + + return 0; +} + + diff --git a/test3.c b/test3.c new file mode 100644 index 0000000..896cccd --- /dev/null +++ b/test3.c @@ -0,0 +1,29 @@ +#include +#include +#include +#include +#include +#include +#include + + + +int main(){ + int i, j, k, N, M; + double **A, *x, *u, **U, *b, **F, *r, *y, *z; + time_t t; + + N = 4; + + b = (double*) malloc(M*sizeof(double)); + F = (double**)malloc((N-1)*sizeof(double*)); + F[0] = b; + for(i = 1; i < N-1; ++i){ + F[i] = F[i-1] + N - 1; + } + + F[1][0] = 26.0; + printf("%f\n", b[N-1]); + + return 0; +} diff --git "a/\345\277\253\351\200\237\344\271\230\346\263\225.c" "b/\345\277\253\351\200\237\344\271\230\346\263\225.c" new file mode 100644 index 0000000..1eb472d --- /dev/null +++ "b/\345\277\253\351\200\237\344\271\230\346\263\225.c" @@ -0,0 +1,291 @@ +#include +#include +#include +#include +#include +#include +int FinSize = 0; + +void FFT_p(complex *y, complex *x, int N){ + int p; + if(N%2 == 0) p = 2; + else if(N%3 == 0) p = 3; + else if(N%5 == 0) p = 5; + else p = N; + + + complex *Nx, *Ny; + int i, j, k, bias, idx; + + bias = N/p; //bias = N/p : °¾²¾¶q + + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(j = 0; j < p; j++){ + for(k = 0; k < bias; k++){ + Nx[k + j*bias] = x[p*k + j]; + } + } + + // F[p][p] + + complex F[p][p]; + for(i = 0; i < p; i++){ + F[0][i] = 1; + F[i][0] = 1; + } + for(i = 1; i < p/2+1; i++){ + for(j = 1; j < p; j++){ + F[i][j] = cexp(-2*M_PI*i*j*I/p); + } + } + for(i = p/2+1; i < p; i++){ //¦³¤Ï¹ïºÙ + for(j = 1; j < p; j++){ + F[i][j] = F[p-i][p-j]; + } + } + + + //recursive + if(N == p){ + for(i = 0; i < p; i++){ + y[i] = 0; + for(j = 0; j < p; j++){ + y[i] += F[i][j]*x[j]; + } + } + } + else{ + + for(j = 0; j < p; j++){ + FFT_p(Ny + j*bias, Nx + j*bias, bias); + } + + for(i = 0; i < bias; i++){ + for(k = 0; k < p; k++){ + Ny[i + k*bias] *= cexp(-2*M_PI*k*i*I/N); + } + + for(j = 0; j < p; j++){ + idx = i + j*bias; + y[idx] = 0; + for(k = 0; k < p; k++){ + y[idx] += F[j][k]*Ny[i + k*bias]; + } + + } + } + + } + + free(Nx); free(Ny); +} + + + +void iFFT_p(complex *y, complex *x, int N){ + int p; + if(N%2 == 0) p = 2; + else if(N%3 == 0) p = 3; + else if(N%5 == 0) p = 5; + else p = N; + + + complex *Nx, *Ny; + int i, j, k, bias, idx; + + bias = N/p; //bias = N/p : °¾²¾¶q + + Nx = (complex*) malloc(N*sizeof(complex)); + Ny = (complex*) malloc(N*sizeof(complex)); + + for(j = 0; j < p; j++){ + for(k = 0; k < bias; k++){ + Nx[k + j*bias] = x[p*k + j]; + } + } + + // F[p][p] + + complex F[p][p]; + for(i = 0; i < p; i++){ + F[0][i] = 1; + F[i][0] = 1; + } + for(i = 1; i < p/2+1; i++){ + for(j = 1; j < p; j++){ + F[i][j] = cexp(2*M_PI*i*j*I/p); + } + } + for(i = p/2+1; i < p; i++){ //¦³¤Ï¹ïºÙ + for(j = 1; j < p; j++){ + F[i][j] = F[p-i][p-j]; + } + } + + + //recursive + if(N == p){ + for(i = 0; i < p; i++){ + y[i] = 0; + for(j = 0; j < p; j++){ + y[i] += F[i][j]*x[j]; + } + } + } + else{ + + for(j = 0; j < p; j++){ + iFFT_p(Ny + j*bias, Nx + j*bias, bias); + } + + for(i = 0; i < bias; i++){ + for(k = 0; k < p; k++){ + Ny[i + k*bias] *= cexp(2*M_PI*k*i*I/N); + } + + for(j = 0; j < p; j++){ + idx = i + j*bias; + y[idx] = 0; + for(k = 0; k < p; k++){ + y[idx] += F[j][k]*Ny[i + k*bias]; + } + + } + } + + } + + free(Nx); free(Ny); +} + + +int* FastMul(char x[], char y[]){ + complex *Cx, *Cy, *u, *v, *ans_t; + int i, j, size_x = 0, size_y = 0; + + // find the max size of these two numbers + + while(x[size_x] != '\0') size_x++; + while(y[size_y] != '\0') size_y++; + + + // add 0s at the first of the small one + if(size_x > size_y){ + FinSize = size_x; + char new_y[FinSize]; + for(i = 0; i < size_x-size_y; i++) new_y[i] = '0'; + for(i = size_x-size_y; i < size_x; i++) new_y[i] = y[i-size_x+size_y]; + y = new_y; + } + else if(size_x < size_y){ + FinSize = size_y; + char new_x[FinSize]; + for(i = 0; i < size_y-size_x; i++) new_x[i] = '0'; + for(i = size_y-size_x; i < size_y; i++) new_x[i] = x[i-size_y+size_x]; + x = new_x; + } + else FinSize = size_x; + + + // use FFT and iFFT to find the answer + Cx = (complex*) malloc(2*FinSize*sizeof(complex)); + Cy = (complex*) malloc(2*FinSize*sizeof(complex)); + u = (complex*) malloc(2*FinSize*sizeof(complex)); + v = (complex*) malloc(2*FinSize*sizeof(complex)); + ans_t = (complex*) malloc(2*FinSize*sizeof(complex)); + + + char a; //¥Î¨Ó¥N´Àchar¡A¥H«KÂন¾ã¼Æ + for(i = 0; i < FinSize; i++){ + a = x[i]; + Cx[FinSize-1-i] = atoi(&a); // change the char to the int + } + for(i = FinSize; i < 2*FinSize; i++) Cx[i] = 0; + + + for(i = 0; i < FinSize; i++){ + a = y[i]; + Cy[FinSize-1-i] = atoi(&a); + } + for(i = FinSize; i < 2*FinSize; i++) Cy[i] = 0; + + //u: x after fft + //v: y after fft + + //after this , only use the in fact 2*N + FinSize = 2*FinSize; + FFT_p(u, Cx, FinSize); + FFT_p(v, Cy, FinSize); + + for(i = 0; i < FinSize; i++){ + u[i]*=v[i]; + ans_t[i] = 0; + } + + + iFFT_p(ans_t, u, FinSize); + for(i = 0; i < FinSize; i++) ans_t[i] /= FinSize; + + + //±Nµ²ªG­Ó§O©ñ¤Jarray¤¤ + int *res = (int*)malloc(FinSize*sizeof(int)); + for(i = 0; i < FinSize; i++) res[i] = 0; + + char ans[3]; + int temp; + for(i = 0; i < FinSize; i++){ + temp = (int)(creal(ans_t[i])+0.5); //+0.5 ¬°¤F¥i¥HµL±ø¥ó±Ë¥h¤p¼Æ³¡¤À //³o³¡¤À¥i¥[³t + itoa(temp, ans, 10); + + int size_ans; + //size of this number + if(ans[1] == '\0' && ans[2] == '\0') size_ans = 1; + else if(ans[2] == '\0') size_ans = 2; + else size_ans = 3; + + for(j = 0; j < size_ans; j++){ + + a = ans[j]; + int idx = FinSize-size_ans-i+j; + res[idx] += atoi(&a); + while(res[idx] >= 10){ + res[idx] -= 10; + res[idx-1] += 1; + } + + } + + } + + int idx = 0; + while(res[idx] == 0) idx++; + // the number is begin from idx + FinSize = FinSize-idx; + int *res2 = (int*)malloc((FinSize)*sizeof(int)); + for(i = 0; i < (FinSize); i++){ + res2[i] = res[i+idx]; + } + + + + return res2; +} + + +int main(){ + int i ; + + char a[6] = {'1', '2', '3', '4', '5', '6'}; + char b[5] = {'5', '6', '7', '8', '9'}; + + int *ans = FastMul(a, b); + + for(i = 0; i < FinSize; i++) printf("%d", ans[i]); + + + + return 0; +} +