00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "config.h"
00022
00023
00024 #include <stdlib.h>
00025 #include <stdio.h>
00026 #include <math.h>
00027
00028
00029 #ifdef HAVE_COMPLEX_H
00030 #include <complex.h>
00031 #endif
00032
00033 #include <fftw3.h>
00034
00035
00036 #include "nfft3.h"
00037 #include "nfft3util.h"
00038
00039
00040 #define KPI2 6.2831853071795864769252867665590057683943387987502
00041
00042 int main(void)
00043 {
00044
00045
00046
00047
00048
00049
00050 const int N = 8;
00051
00052
00053
00054
00055
00056 fpt_set set = fpt_init(1,lrint(ceil(log2((double)N))),0U);
00057
00058
00059 double *alpha = malloc((N+2)*sizeof(double)),
00060 *beta = malloc((N+2)*sizeof(double)),
00061 *gamma = malloc((N+2)*sizeof(double));
00062
00063
00064 alpha[0] = beta[0] = 0.0;
00065
00066 gamma[0] = 1.0;
00067
00068
00069 {
00070 int k;
00071 for (k = 0; k <= N; k++)
00072 {
00073 alpha[k+1] = ((double)(2*k+1))/((double)(k+1));
00074 beta[k+1] = 0.0;
00075 gamma[k+1] = -((double)(k))/((double)(k+1));
00076 }
00077 }
00078
00079 printf(
00080 "Computing a fast polynomial transform (FPT) and a fast discrete cosine \n"
00081 "transform (DCT) to evaluate\n\n"
00082 " f_j = a_0 P_0(x_j) + a_1 P_1(x_j) + ... + a_N P_N(x_j), j=0,1,...,N,\n\n"
00083 "with N=%d, x_j = cos(j*pi/N), j=0,1,...N, the Chebyshev nodes, a_k,\n"
00084 "k=0,1,...,N, random Fourier coefficients in [-1,1]x[-1,1]*I, and P_k,\n"
00085 "k=0,1,...,N, the Legendre polynomials.",N
00086 );
00087
00088
00089 nfft_srand48(314);
00090
00091
00092
00093
00094
00095
00096
00097
00098 fpt_precompute(set,0,alpha,beta,gamma,0,1000.0);
00099
00100
00101 {
00102
00103 double _Complex *a = malloc((N+1)*sizeof(double _Complex));
00104 double _Complex *b = malloc((N+1)*sizeof(double _Complex));
00105 double *f = malloc((N+1)*sizeof(double _Complex));
00106
00107
00108 const int NP1 = N + 1;
00109 fftw_r2r_kind kind = FFTW_REDFT00;
00110 fftw_plan p = fftw_plan_many_r2r(1, &NP1, 1, (double*)b, NULL, 2, 1,
00111 (double*)f, NULL, 1, 1, &kind, 0U);
00112
00113
00114 {
00115 int k;
00116 printf("\n2) Random Fourier coefficients a_k, k=0,1,...,N:\n");
00117 for (k = 0; k <= N; k++)
00118 {
00119 a[k] = 2.0*nfft_drand48() - 1.0;
00120 printf(" a_%-2d = %+5.3lE\n",k,creal(a[k]));
00121 }
00122 }
00123
00124
00125 fpt_trafo(set,0,a,b,N,0U);
00126
00127
00128
00129
00130
00131 {
00132 int j;
00133 for (j = 1; j < N; j++)
00134 b[j] *= 0.5;
00135 }
00136
00137
00138 fftw_execute(p);
00139
00140 {
00141 int j;
00142 printf("\n3) Function values f_j, j=1,1,...,M:\n");
00143 for (j = 0; j <= N; j++)
00144 printf(" f_%-2d = %+5.3lE\n",j,f[j]);
00145 }
00146
00147
00148 free(a);
00149 free(b);
00150 free(f);
00151
00152
00153 fftw_destroy_plan(p);
00154 }
00155
00156
00157 fpt_finalize(set);
00158 free(alpha);
00159 free(beta);
00160 free(gamma);
00161
00162 return EXIT_SUCCESS;
00163 }