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