00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <math.h>
00022 #include <stdio.h>
00023 #include "infft.h"
00024 #include "nfft3util.h"
00025 #include "legendre.h"
00026
00027
00028 DK(KSQRTPII,0.56418958354775628694807945156077258584405062932900);
00029
00030 static inline R alpha_al(const int k, const int n)
00031 {
00032 if (k > 0)
00033 {
00034 if (k < n)
00035 return IF(k%2,K(1.0),K(-1.0));
00036 else
00037 return SQRT(((R)(2*k+1))/((R)(k-n+1)))*SQRT((((R)(2*k+1))/((R)(k+n+1))));
00038 }
00039 else if (k == 0)
00040 {
00041 if (n == 0)
00042 return K(1.0);
00043 else
00044 return IF(n%2,K(0.0),K(-1.0));
00045 }
00046 return K(0.0);
00047 }
00048
00049 static inline R beta_al(const int k, const int n)
00050 {
00051 if (0 <= k && k < n)
00052 return K(1.0);
00053 else
00054 return K(0.0);
00055 }
00056
00057 static inline R gamma_al(const int k, const int n)
00058 {
00059 if (k == -1)
00060 return SQRT(KSQRTPII*nfft_lambda((R)(n),K(0.5)));
00061 else if (k <= n)
00062 return K(0.0);
00063 else
00064 return -SQRT(((R)(k-n))/((R)(k-n+1))*((R)(k+n))/((R)(k+n+1)));
00065 }
00066
00067 void alpha_al_row(R *alpha, const int N, const int n)
00068 {
00069 int j;
00070 R *p = alpha;
00071 for (j = -1; j <= N; j++)
00072 *p++ = alpha_al(j,n);
00073 }
00074
00075 void beta_al_row(R *beta, const int N, const int n)
00076 {
00077 int j;
00078 R *p = beta;
00079 for (j = -1; j <= N; j++)
00080 *p++ = beta_al(j,n);
00081 }
00082
00083 void gamma_al_row(R *gamma, const int N, const int n)
00084 {
00085 int j;
00086 R *p = gamma;
00087 for (j = -1; j <= N; j++)
00088 *p++ = gamma_al(j,n);
00089 }
00090
00091 inline void alpha_al_all(R *alpha, const int N)
00092 {
00093 int i,j;
00094 R *p = alpha;
00095 for (i = 0; i <= N; i++)
00096 for (j = -1; j <= N; j++)
00097 *p++ = alpha_al(j,i);
00098 }
00099
00100 inline void beta_al_all(R *alpha, const int N)
00101 {
00102 int i,j;
00103 R *p = alpha;
00104 for (i = 0; i <= N; i++)
00105 for (j = -1; j <= N; j++)
00106 *p++ = beta_al(j,i);
00107 }
00108
00109 inline void gamma_al_all(R *alpha, const int N)
00110 {
00111 int i,j;
00112 R *p = alpha;
00113 for (i = 0; i <= N; i++)
00114 for (j = -1; j <= N; j++)
00115 *p++ = gamma_al(j,i);
00116 }
00117
00118 void eval_al(R *x, R *y, const int size, const int k, R *alpha,
00119 R *beta, R *gamma)
00120 {
00121
00122
00123
00124 int i,j;
00125 R a,b,x_val_act,a_old;
00126 R *x_act, *y_act;
00127 R *alpha_act, *beta_act, *gamma_act;
00128
00129
00130 x_act = x;
00131 y_act = y;
00132 for (i = 0; i < size; i++)
00133 {
00134 a = 1.0;
00135 b = 0.0;
00136 x_val_act = *x_act;
00137
00138 if (k == 0)
00139 {
00140 *y_act = 1.0;
00141 }
00142 else
00143 {
00144 alpha_act = &(alpha[k]);
00145 beta_act = &(beta[k]);
00146 gamma_act = &(gamma[k]);
00147 for (j = k; j > 1; j--)
00148 {
00149 a_old = a;
00150 a = b + a_old*((*alpha_act)*x_val_act+(*beta_act));
00151 b = a_old*(*gamma_act);
00152 alpha_act--;
00153 beta_act--;
00154 gamma_act--;
00155 }
00156 *y_act = (a*((*alpha_act)*x_val_act+(*beta_act))+b);
00157 }
00158 x_act++;
00159 y_act++;
00160 }
00161 }
00162
00163 int eval_al_thresh(R *x, R *y, const int size, const int k, R *alpha,
00164 R *beta, R *gamma, R threshold)
00165 {
00166
00167
00168
00169 int i,j;
00170 R a,b,x_val_act,a_old;
00171 R *x_act, *y_act;
00172 R *alpha_act, *beta_act, *gamma_act;
00173
00174
00175 x_act = x;
00176 y_act = y;
00177 for (i = 0; i < size; i++)
00178 {
00179 a = 1.0;
00180 b = 0.0;
00181 x_val_act = *x_act;
00182
00183 if (k == 0)
00184 {
00185 *y_act = 1.0;
00186 }
00187 else
00188 {
00189 alpha_act = &(alpha[k]);
00190 beta_act = &(beta[k]);
00191 gamma_act = &(gamma[k]);
00192 for (j = k; j > 1; j--)
00193 {
00194 a_old = a;
00195 a = b + a_old*((*alpha_act)*x_val_act+(*beta_act));
00196 b = a_old*(*gamma_act);
00197 alpha_act--;
00198 beta_act--;
00199 gamma_act--;
00200 }
00201 *y_act = (a*((*alpha_act)*x_val_act+(*beta_act))+b);
00202 if (fabs(*y_act) > threshold)
00203 {
00204 return 1;
00205 }
00206 }
00207 x_act++;
00208 y_act++;
00209 }
00210 return 0;
00211 }