Originally posted by MrFun
head hurts . . . . . reading this thread . . . . . . .
head hurts . . . . . reading this thread . . . . . . .
/************************************************************************/ /*** In the relation P_{L+r} = A_r^L P_L + B_r^L P_{L-1} this routine generates B_r^L . Initial conditions are B_0^L = 0, B_{-1}^L = 1. m = order of the problem r = how much to shift l = where to start shift (i.e. l = L in above *) n = number of points desired result is double array of size n workspace is double array of size 4 * n ***/ static void Bshifted( int m, int r, int l, int n, double *result, double *workspace ) { double *cospts, *pi, *piminus1, *pcurr; double di, dn, aval, cval; int i, j; /* allocate work space */ cospts = workspace; pi = cospts + n; piminus1 = pi + n; pcurr = piminus1 + n; /* generate sample points */ dn = (double) n; for (i=0; i<n; i++) { di = (double) i; cospts[i] = cos( (((2.0 * di)+1.0) * PI)/(2.0 * dn)); } if (r < 0) { for (i=0; i<n; i++) result[i] = 1.0; } /* piminus1 starts off as all zeros */ for (i=0; i<n; i++) piminus1[i] = 0.0; if (r == 0) { for (i=0; i<n; i++) result[i] = 0.0; } /* return (piminus1); */ /* generate initial value of pi */ for (i=0; i<n; i++) pi[i] = L2_cn(m,m+l); if (r == 1) { for (i=0; i<n; i++) result[i] = pi[i]; } /* return(pi); */ /* main loop */ if (r > 1) { for (i=m+1; i<m+r; i++) { aval = L2_an(m,i+l); cval = L2_cn(m,i+l); for (j=0; j<n; j++) pcurr[j] = aval * cospts[j] * pi[j] + cval * piminus1[j]; /* need to copy values */ memcpy(piminus1, pi, sizeof(double) * n); memcpy(pi, pcurr, sizeof(double) * n); } memcpy(result, pcurr, sizeof(double) * n); } /* return (pcurr); */ } static void A_revshifted( int m, int r, int l, int n, double *result, double *workspace ) { double *cospts, *pi, *piplus1, *pcurr; double di, dn, aval, cval; int i, j; /* allocate work space */ cospts = workspace; pi = cospts + n; piplus1 = pi + n; pcurr = piplus1 + n; /* generate sample points */ dn = (double) n; for (i=0; i<n; i++) { di = (double) i; cospts[i] = cos((((2.0 * di)+1.0) * PI)/(2.0 * dn)); } if (r < 0) { for (i=0; i<n; i++) result[i] = 0.0; } /* return (zeros(n)); */ /* piplus1 starts off as all ones */ for (i=0; i<n; i++) piplus1[i] = 1.0; if (r == 0) { for (i=0; i<n; i++) result[i] = 1.0; } /* return (piplus1); */ /* generate initial value of pi */ /* aval = -L2_an(m,m+l); */ /* cval = 1.0 / L2_cn(m,m+l); */ aval = L2_ancn(m,m+l); for (i=0; i<n; i++) pi[i] = cospts[i] * aval; if (r == 1) { for (i=0; i<n; i++) result[i] = pi[i]; } /* return(pi); */ /* main loop */ if (r > 1) { for(i = 1 ; i < r ; i ++){ /* aval = -L2_an(m, m + l - i); */ aval = L2_ancn(m, m + l - i); /* cval = 1.0/L2_cn(m, m + l - i); */ cval = L2_cn_inv(m,m+l-i); for(j = 0 ; j < n ; j ++) pcurr[j] = aval * cospts[j] * pi[j] + cval * piplus1[j]; /* need to copy values */ memcpy(piplus1, pi, sizeof(double) * n); memcpy(pi, pcurr, sizeof(double) * n); } /*** now copy the final result ***/ memcpy(result, pcurr, sizeof(double) * n); } /* return (pcurr); */ }
Comment