Numerical recipes in C++

				#include "nr.h"
				
				void NR::spread(const DP y, Vec_IO_DP &yy, const DP x, const int m)
				{
					static int nfac[11]={0,1,1,2,6,24,120,720,5040,40320,362880};
					int ihi,ilo,ix,j,nden;
					DP fac;
				
					int n=yy.size();
					if (m > 10) nrerror("factorial table too small in spread");
					ix=int(x);
					if (x == DP(ix)) yy[ix-1] += y;
					else {
						ilo=MIN(MAX(int(x-0.5*m),0),int(n-m));
						ihi=ilo+m;
						nden=nfac[m];
						fac=x-ilo-1;
						for (j=ilo+1;j						yy[ihi-1] += y*fac/(nden*(x-ihi));
						for (j=ihi-1;j>ilo;j--) {
							nden=(nden/(j-ilo))*(j-ihi);
							yy[j-1] += y*fac/(nden*(x-j));
						}
					}
				}