2.8 Vandermonde Matrices and Toeplitz Matrices 95 #include "nrutil.h" #define FREERETURN {free_vector(h,1,n);free_vector(g,1,n);return;} void toeplz(float r[],float x[],float y],int n) Solves the Toeplitz systemR(N+)=(i=1,....N).The Toeplitz matrix need not be symmetric.y[1..n]and r[1..2*n-1]are input arrays;x[1..n]is the output array. int j,k,m,m1,m2 float pp,pt1,pt2,qq,qt1,qt2,sd,sgd,sgn,shn,sxn; float *g,*h; http://www.nr. readable files if (r[n]=0.0)nrerror("toeplz-1 singular principal minor"); Permission is opyright(C Sample page g=vector(1,n); h-vector(1,n); x[1]=y[1]/r[n]; Initialize for the recursion. 83g if (n=1)FREERETURN granted fori g[1]=r[n-1]/x[n]; h[1]=r[n+1]/x[n]; for (m=1;m<=n;m++){ Main loop over the recursion. 11-800-872 (including this one) m1=m+1; sxn =-y[m1]; Compute numerator and denominator for x, sd =-r[n] 7423 (North America to any server computer,is tusers to make one paper 1988-1992 by Cambridge University Press.Programs from NUMERICAL RECIPES IN C: for(j=1;j<=m;j++){ sxn +r[n+m1-j]*x[j]; sd +r[n+m1-j]*g[m-j+1]: 2 THE if (sd ==0.0)nrerror("toeplz-2 singular principal minor"); x[m1]=sxn/sd; wnence t. for (j=1;j<=m;j++)x[j] -=x[m1]*g[m-j+1] only). if (m1=n)FREERETURN sgn -r[n-m1]; Compute numerator and denominator for G and H, shn =-r[n+m1]; send sgd =-r[n]; for(j=1;j<=m;j++)[ email sgn +r[n+j-m1]*g[j]; shn +r[n+m1-j]*h[j]; sgd +=r[n+j-m1]*h[m-j+1]; if (sgd ==0.0)nrerror("toeplz-3 singular principal minor"); g[m1】=sgm/sgd; whence G and H. ART OF SCIENTIFIC COMPUTING(ISBN 0-521- h[m1]=shn/sd; k=m; m2=(m+1)>1; pp=g[m1]: Ito directcustserv@cambridge.org 88 personal use.Further reproduction,or 1988-1992 by Numerical Recipes (j=1;j<=m2;j+)( -431085 pt1=g[j]; pt2=g [k]: qti=h[j]; qt2=h [k]; Software. g[j]=pt1-pp*qt2; gk幻=pt2-Pp*qt1; h[j]=qt1-qq*pt2; (outside North America) h[k--]=qt2-qq*pt1 Back for another recurrence. nrerror("toeplz should not arrive here!"); If you are in the business of solving very large Toeplitz systems,you should find out about so-called "new,fast"algorithms,which require only on the order of N(log N)operations compared to N-for Levinson's method.These methods are too complicated to include here2.8 Vandermonde Matrices and Toeplitz Matrices 95 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). #include "nrutil.h" #define FREERETURN {free_vector(h,1,n);free_vector(g,1,n);return;} void toeplz(float r[], float x[], float y[], int n) Solves the Toeplitz system N j=1 R(N+i−j)xj = yi (i = 1,...,N). The Toeplitz matrix need not be symmetric. y[1..n] and r[1..2*n-1] are input arrays; x[1..n] is the output array. { int j,k,m,m1,m2; float pp,pt1,pt2,qq,qt1,qt2,sd,sgd,sgn,shn,sxn; float *g,*h; if (r[n] == 0.0) nrerror("toeplz-1 singular principal minor"); g=vector(1,n); h=vector(1,n); x[1]=y[1]/r[n]; Initialize for the recursion. if (n == 1) FREERETURN g[1]=r[n-1]/r[n]; h[1]=r[n+1]/r[n]; for (m=1;m<=n;m++) { Main loop over the recursion. m1=m+1; sxn = -y[m1]; Compute numerator and denominator for x, sd = -r[n]; for (j=1;j<=m;j++) { sxn += r[n+m1-j]*x[j]; sd += r[n+m1-j]*g[m-j+1]; } if (sd == 0.0) nrerror("toeplz-2 singular principal minor"); x[m1]=sxn/sd; whence x. for (j=1;j<=m;j++) x[j] -= x[m1]*g[m-j+1]; if (m1 == n) FREERETURN sgn = -r[n-m1]; Compute numerator and denominator for G and H, shn = -r[n+m1]; sgd = -r[n]; for (j=1;j<=m;j++) { sgn += r[n+j-m1]*g[j]; shn += r[n+m1-j]*h[j]; sgd += r[n+j-m1]*h[m-j+1]; } if (sgd == 0.0) nrerror("toeplz-3 singular principal minor"); g[m1]=sgn/sgd; whence G and H. h[m1]=shn/sd; k=m; m2=(m+1) >> 1; pp=g[m1]; qq=h[m1]; for (j=1;j<=m2;j++) { pt1=g[j]; pt2=g[k]; qt1=h[j]; qt2=h[k]; g[j]=pt1-pp*qt2; g[k]=pt2-pp*qt1; h[j]=qt1-qq*pt2; h[k--]=qt2-qq*pt1; } } Back for another recurrence. nrerror("toeplz - should not arrive here!"); } If you are in the business of solving very large Toeplitz systems, you should find out about so-called “new, fast” algorithms, which require only on the order of N(log N) 2 operations, compared to N2 for Levinson’s method. These methods are too complicated to include here