diff -rupN code/cholesky.h code302/cholesky.h
--- code/cholesky.h	Thu May 10 14:24:12 2007
+++ code302/cholesky.h	Tue Sep 16 17:25:12 2008
@@ -3,7 +3,6 @@ struct Cholesky{
 	MatDoub el;
 	Cholesky(MatDoub_I &a) : n(a.nrows()), el(a) {
 		Int i,j,k;
-		VecDoub tmp;
 		Doub sum;
 		if (el.ncols() != n) throw("need square matrix");
 		for (i=0;i<n;i++) {
diff -rupN code/deviates.h code302/deviates.h
--- code/deviates.h	Thu May 10 14:24:12 2007
+++ code302/deviates.h	Tue Sep 16 17:25:12 2008
@@ -99,7 +99,6 @@ struct Gammadev : Normaldev {
 struct Poissondev : Ran {
 	Doub lambda, sqlam, loglam, lamexp, lambold;
 	VecDoub logfact;
-	Int swch;
 	Poissondev(Doub llambda, Ullong i) : Ran(i), lambda(llambda),
 		logfact(1024,-1.), lambold(-1.) {}
 	Int dev() {
diff -rupN code/dlinmin.h code302/dlinmin.h
--- code/dlinmin.h	Mon Apr 16 17:39:41 2007
+++ code302/dlinmin.h	Wed Dec 31 17:00:00 1969
@@ -1,25 +0,0 @@
-template <class T>
-struct Dlinmin {
-	VecDoub &p;
-	VecDoub &xi;
-	Doub &fret;
-	Df1dim<T> df1dim;
-	Mnbrak mnbrak;
-	Dbrent dbrent;
-	Dlinmin(VecDoub_IO &pp, VecDoub_IO &xii, Doub &frett, T &funcd) :
-		p(pp), xi(xii), fret(frett), df1dim(p,xi,funcd) {}
-	void min()
-	{
-		const Doub TOL=1.0e-8;
-		Doub xx,xmin,fx,fb,fa,bx,ax;
-		Int n=p.size();
-		ax=0.0;
-		xx=1.0;
-		mnbrak.bracket(ax,xx,bx,fa,fx,fb,df1dim);
-		fret=dbrent.min(ax,xx,bx,df1dim,TOL,xmin);
-		for (Int j=0;j<n;j++) {
-			xi[j] *= xmin;
-			p[j] += xi[j];
-		}
-	}
-};
diff -rupN code/factorial.h code302/factorial.h
--- code/factorial.h	Mon Apr 16 17:39:41 2007
+++ code302/factorial.h	Wed Dec 31 17:00:00 1969
@@ -1,34 +0,0 @@
-struct Factorial {
-	static const Int NLTAB = 10000;
-	VecDoub tab;
-	VecDoub ltab;
-
-	Factorial() : tab(171), ltab(NLTAB) {
-		Int i;
-		tab[0] = 1.;
-		ltab[0] = 0.;
-		for (i=1;i<171;i++) {
-			tab[i] = i*tab[i-1];
-			ltab[i] = log(tab[i]);
-		}
-		for (i=171;i<NLTAB;i++) ltab[i] = gammln(i+1.);
-	}
-
-	inline Doub fac(const Int i) {
-		if (i<0 || i>170) throw("factorial overflows");
-		return tab[i];
-	}
-
-	Doub facln(const Int i) {
-		if (i<0) throw("negative facln arg");
-		if (i<NLTAB) return ltab[i];
-		return gammln(i+1.);
-	}
-
-	Doub bico(const Int n, const Int k) {
-		if (n<0 || k<0 || k>n) throw("bico bad args");
-		if (n<171) return floor(0.5+tab[n]/(tab[k]*tab[n-k]));
-		if (n<NLTAB) return floor(0.5+exp(ltab[n]-ltab[k]-ltab[n-k]));
-		return floor(0.5+exp(gammln(n+1.)-gammln(k+1.)-gammln(n-k+1.)));
-	}
-};
diff -rupN code/hash.h code302/hash.h
--- code/hash.h	Thu May 10 14:24:13 2007
+++ code302/hash.h	Tue Sep 16 17:25:12 2008
@@ -10,7 +10,7 @@ struct Hashfn1 {
 				k = (Uint *)key;
 				return hasher.int64(*k);
 			case 8:
-				kk = (Ullong *)key;\
+				kk = (Ullong *)key;
 				return hasher.int64(*kk);
 			default:
 				throw("Hashfn1 is for 4 or 8 byte keys only.");
@@ -51,7 +51,7 @@ template<class keyT, class hfnT> struct 
 	VecInt htable, next, garbg;
 	VecUllong thehash;
 	hfnT hash;
-	Hashtable(Int nh, Int nv);
+	Hashtable(Int nh, Int nm);
 	Int iget(const keyT &key);
 	Int iset(const keyT &key);
 	Int ierase(const keyT &key);
@@ -60,9 +60,9 @@ template<class keyT, class hfnT> struct 
 };
 
 template<class keyT, class hfnT>
-Hashtable<keyT,hfnT>::Hashtable(Int nh, Int nv):
-	hash(sizeof(keyT)), nhash(nh), nmax(nv), nn(0), ng(0),
-	htable(nh), next(nv), garbg(nv), thehash(nv) {
+Hashtable<keyT,hfnT>::Hashtable(Int nh, Int nm):
+	hash(sizeof(keyT)), nhash(nh), nmax(nm), nn(0), ng(0),
+	htable(nh), next(nm), garbg(nm), thehash(nm) {
 	for (Int j=0; j<nh; j++) { htable[j] = -1; }
 }
 template<class keyT, class hfnT>
diff -rupN code/kdtree.h code302/kdtree.h
--- code/kdtree.h	Thu May 10 14:24:13 2007
+++ code302/kdtree.h	Tue Sep 16 17:25:12 2008
@@ -53,7 +53,6 @@ template<Int DIM> struct KDtree {
 	Doub disti(Int jpt, Int kpt);
 	Int locate(Point<DIM> pt);
 	Int locate(Int jpt);
-	Int nearest(Int jpt);
 	Int nearest(Point<DIM> pt);
 	void nnearest(Int jpt, Int *nn, Doub *dn, Int n);
 	static void sift_down(Doub *heap, Int *ndx, Int nn);
diff -rupN code/linmin.h code302/linmin.h
--- code/linmin.h	Mon Apr 16 17:39:41 2007
+++ code302/linmin.h	Wed Dec 31 17:00:00 1969
@@ -1,25 +0,0 @@
-template <class T>
-struct Linmin {
-	VecDoub &p;
-	VecDoub &xi;
-	Doub &fret;
-	F1dim<T> f1dim;
-	Mnbrak mnbrak;
-	Brent brent;
-	Linmin(VecDoub_IO &pp, VecDoub_IO &xii, Doub &frett, T &func) :
-		p(pp), xi(xii), fret(frett), f1dim(p,xi,func) {}
-	void min()
-	{
-		const Doub TOL=1.0e-8;
-		Doub xx,xmin,fx,fb,fa,bx,ax;
-		Int n=p.size();
-		ax=0.0;
-		xx=1.0;
-		mnbrak.bracket(ax,xx,bx,fa,fx,fb,f1dim);
-		fret=brent.min(ax,xx,bx,f1dim,TOL,xmin);
-		for (Int j=0;j<n;j++) {
-			xi[j] *= xmin;
-			p[j] += xi[j];
-		}
-	}
-};
diff -rupN code/ludcmp.h code302/ludcmp.h
--- code/ludcmp.h	Thu May 10 14:24:13 2007
+++ code302/ludcmp.h	Tue Sep 16 17:25:12 2008
@@ -27,6 +27,7 @@ LUdcmp::LUdcmp(MatDoub_I &a) : n(a.nrows
 	}
 	for (k=0;k<n;k++) {
 		big=0.0;
+		imax=k;
 		for (i=k;i<n;i++) {
 			temp=vv[i]*abs(lu[i][k]);
 			if (temp > big) {
diff -rupN code/mins_ndim.h code302/mins_ndim.h
--- code/mins_ndim.h	Thu May 10 14:24:13 2007
+++ code302/mins_ndim.h	Tue Sep 16 17:25:12 2008
@@ -187,7 +187,7 @@ struct Frprmn : Linemethod<T> {
 			fp=fret;
 			func.df(p,xi);
 			Doub test=0.0;
-			Doub den=MAX(fp,1.0);
+			Doub den=MAX(abs(fp),1.0);
 			for (Int j=0;j<n;j++) {
 				Doub temp=abs(xi[j])*MAX(abs(p[j]),1.0)/den;
 				if (temp > test) test=temp;
diff -rupN code/mnbrak.h code302/mnbrak.h
--- code/mnbrak.h	Mon Apr 16 17:39:41 2007
+++ code302/mnbrak.h	Wed Dec 31 17:00:00 1969
@@ -1,60 +0,0 @@
-struct Mnbrak {
-	inline void shft3(Doub &a, Doub &b, Doub &c, const Doub d)
-	{
-		a=b;
-		b=c;
-		c=d;
-	}
-	template <class T>
-	void bracket(Doub &ax, Doub &bx, Doub &cx, Doub &fa, Doub &fb, Doub &fc,
-		T &func)
-	{
-		const Doub GOLD=1.618034,GLIMIT=100.0,TINY=1.0e-20;
-		Doub fu;
-		fa=func(ax);
-		fb=func(bx);
-		if (fb > fa) {
-			SWAP(ax,bx);
-			SWAP(fb,fa);
-		}
-		cx=bx+GOLD*(bx-ax);
-		fc=func(cx);
-		while (fb > fc) {
-			Doub r=(bx-ax)*(fb-fc);
-			Doub q=(bx-cx)*(fb-fa);
-			Doub u=bx-((bx-cx)*q-(bx-ax)*r)/
-				(2.0*SIGN(MAX(abs(q-r),TINY),q-r));
-			Doub ulim=bx+GLIMIT*(cx-bx);
-			if ((bx-u)*(u-cx) > 0.0) {
-				fu=func(u);
-				if (fu < fc) {
-					ax=bx;
-					bx=u;
-					fa=fb;
-					fb=fu;
-					return;
-				} else if (fu > fb) {
-					cx=u;
-					fc=fu;
-					return;
-				}
-				u=cx+GOLD*(cx-bx);
-				fu=func(u);
-			} else if ((cx-u)*(u-ulim) > 0.0) {
-				fu=func(u);
-				if (fu < fc) {
-					shft3(bx,cx,u,u+GOLD*(u-cx));
-					shft3(fb,fc,fu,func(u));
-				}
-			} else if ((u-ulim)*(ulim-cx) >= 0.0) {
-				u=ulim;
-				fu=func(u);
-			} else {
-				u=cx+GOLD*(cx-bx);
-				fu=func(u);
-			}
-			shft3(ax,bx,cx,u);
-			shft3(fa,fb,fc,fu);
-		}
-	}
-};
diff -rupN code/mparith.h code302/mparith.h
--- code/mparith.h	Thu May 10 14:24:13 2007
+++ code302/mparith.h	Tue Sep 16 17:25:12 2008
@@ -228,7 +228,7 @@ void MParith::mp2dfr(VecUchar_IO &a, str
 	}
 }
 string MParith::mppi(const Int np) {
-	const Uint IAOFF=48,MACC=2;
+	const Uint MACC=2;
 	Int ir,j,n=np+MACC;
 	Uchar mm;
 	string s;
diff -rupN code/pade.h code302/pade.h
--- code/pade.h	Thu May 10 14:24:13 2007
+++ code302/pade.h	Tue Sep 16 17:25:12 2008
@@ -1,10 +1,8 @@
 Ratfn pade(VecDoub_I &cof)
 {
-	const Doub BIG=1.0e99;
 	Int j,k,n=(cof.size()-1)/2;
 	Doub sum;
 	MatDoub q(n,n),qlu(n,n);
-	VecInt indx(n);
 	VecDoub x(n),y(n),num(n+1),denom(n+1);
 	for (j=0;j<n;j++) {
 		y[j]=cof[n+j+1];
diff -rupN code/psplot.h code302/psplot.h
--- code/psplot.h	Thu May 10 14:24:13 2007
+++ code302/psplot.h	Tue Sep 16 17:25:12 2008
@@ -1,7 +1,7 @@
 struct PSpage {
 
-	FILE *PLT;
-	char *file;
+	static FILE *PLT;
+	static char *file;
 	char fontname[128];
 	Doub fontsize;
 
@@ -17,7 +17,7 @@ struct PSpage {
 		setfont("Times-Roman",12.);
 		setlinewidth(0.5);
 	}
-	PSpage(FILE *PPLT, char *ffile) : PLT(PPLT), file(ffile){}
+	PSpage() {}
 	~PSpage() {if (PLT) close();}
 
 	void setfont(char *fontnam, Doub size) {
@@ -97,8 +97,7 @@ struct PSplot : PSpage {
 	Doub majticsz,minticsz;
 
 	PSplot(PSpage &page, Doub ppll, Doub ppur, Doub qqll, Doub qqur)
-		: PSpage(page.PLT,page.file),
-		pll(ppll), qll(qqll), pur(ppur), qur(qqur),
+		: pll(ppll), qll(qqll), pur(ppur), qur(qqur),
 		xll(ppll), yll(qqll), xur(ppur), yur(qqur), xbox(4), ybox(4),
 		majticsz(8.), minticsz(4.) {
 			strcpy(fontname,page.fontname);
@@ -200,3 +199,5 @@ struct PSplot : PSpage {
 		scales(xmajd,xmind,ymajd,ymind);
 	}
 };
+FILE *PSpage::PLT;
+char *PSpage::file;
diff -rupN code/quasinewton.h code302/quasinewton.h
--- code/quasinewton.h	Thu May 10 14:24:13 2007
+++ code302/quasinewton.h	Tue Sep 16 17:25:12 2008
@@ -36,7 +36,7 @@ void dfpmin(VecDoub_IO &p, const Doub gt
 		for (Int i=0;i<n;i++) dg[i]=g[i];
 		funcd.df(p,g);
 		test=0.0;
-		den=MAX(fret,1.0);
+		den=MAX(abs(fret),1.0);
 		for (Int i=0;i<n;i++) {
 			temp=abs(g[i])*MAX(abs(p[i]),1.0)/den;
 			if (temp > test) test=temp;
diff -rupN code/ran.h code302/ran.h
--- code/ran.h	Thu May 10 14:24:13 2007
+++ code302/ran.h	Tue Sep 16 17:25:12 2008
@@ -132,4 +132,3 @@ struct Ranlim32 {
 		2.32830643653869629E-10 * int32() );
 	}
 };
-
diff -rupN code/rlftfrag.h code302/rlftfrag.h
--- code/rlftfrag.h	Mon Apr 16 17:39:41 2007
+++ code302/rlftfrag.h	Wed Dec 31 17:00:00 1969
@@ -1,60 +0,0 @@
-Int main(void)		// example1
-{
-	const Int N2=256,N3=256;
-	NRMat3d<Doub> data(1,N2,N3);
-	MatDoub speq(1,2*N2);
-
-//	...
-	rlft3(data,speq,1);
-//	...
-	rlft3(data,speq,-1);
-//	...
-	return 0;
-}
-
-Int main(void)		// example2
-{
-	const Int N1=32,N2=64,N3=16;
-	NRMat3d<Doub> data(N1,N2,N3);
-	MatDoub speq(N1,2*N2);
-//	...
-	rlft3(data,speq,1);
-//	...
-	return 0;
-}
-
-Int main(void)		// example3
-{
-	Int j;
-	Doub fac,r,i,*sp1,*sp2;
-	const Int N=32;
-	NRMat3d<Doub> data1(N,N,N),data2(N,N,N);
-	MatDoub speq1(N,2*N),speq2(N,2*N);
-//	...
-	rlft3(data1,speq1,1);
-	rlft3(data2,speq2,1);
-	fac=2.0/(N*N*N);
-	sp1 = &data1[0][0][0];
-	sp2 = &data1[0][0][0];
-	for (j=0;j<N*N*N/2;j++) {
-		r = sp1[0]*sp2[0] - sp1[1]*sp2[1];
-		i = sp1[0]*sp2[1] + sp1[1]*sp2[0];
-		sp1[0] = fac*r;
-		sp1[1] = fac*i;
-		sp1 += 2;
-		sp2 += 2;
-	}
-	sp1 = &speq1[0][0];
-	sp2 = &speq2[0][0];
-	for (j=0;j<N*N;j++) {
-		r = sp1[0]*sp2[0] - sp1[1]*sp2[1];
-		i = sp1[0]*sp2[1] + sp1[1]*sp2[0];
-		sp1[0] = fac*r;
-		sp1[1] = fac*i;
-		sp1 += 2;
-		sp2 += 2;
-	}
-	rlft3(data1,speq1,-1);
-//	...
-	return 0;
-}
diff -rupN code/roots.h code302/roots.h
--- code/roots.h	Thu May 10 14:24:13 2007
+++ code302/roots.h	Tue Sep 16 17:25:12 2008
@@ -272,8 +272,8 @@ Doub rtsafe(T &funcd, const Doub x1, con
 			if (temp == rts) return rts;
 		}
 		if (abs(dx) < xacc) return rts;
-		Doub f=funcd(rts);
-		Doub df=funcd.df(rts);
+		f=funcd(rts);
+		df=funcd.df(rts);
 		if (f < 0.0)
 			xl=rts;
 		else
diff -rupN code/roots_multidim.h code302/roots_multidim.h
--- code/roots_multidim.h	Thu May 10 14:24:13 2007
+++ code302/roots_multidim.h	Tue Sep 16 17:25:12 2008
@@ -180,7 +180,12 @@ void broydn(VecDoub_IO &x, Bool &check, 
 	for (its=1;its<=MAXITS;its++) {
 		if (restrt) {
 			qr=new QRdcmp(fdjac(x,fvec));
-			if (qr->sing) throw("singular Jacobian in broydn");
+			if (qr->sing) {
+				MatDoub one(n,n,0.0);
+				for (i=0;i<n;i++) one[i][i]=1.0;
+				delete qr;
+				qr=new QRdcmp(one);
+			}
 		} else {
 			for (i=0;i<n;i++) s[i]=x[i]-xold[i];
 			for (i=0;i<n;i++) {
@@ -215,6 +220,12 @@ void broydn(VecDoub_IO &x, Bool &check, 
 		}
 		fold=f;
 		qr->rsolve(p,p);
+		Doub slope=0.0;
+		for (i=0;i<n;i++) slope += g[i]*p[i];
+		if (slope >= 0.0) {
+			restrt=true;
+			continue;
+		}
 		lnsrch(xold,fold,g,p,x,f,stpmax,check,fmin);
 		test=0.0;
 		for (i=0;i<n;i++)
diff -rupN code/savgol.h code302/savgol.h
--- code/savgol.h	Thu May 10 14:24:13 2007
+++ code302/savgol.h	Tue Sep 16 17:25:12 2008
@@ -5,7 +5,6 @@ void savgol(VecDoub_O &c, const Int np, 
 	Doub fac,sum;
 	if (np < nl+nr+1 || nl < 0 || nr < 0 || ld > m || nl+nr < m)
 		throw("bad args in savgol");
-	VecInt indx(m+1);
 	MatDoub a(m+1,m+1);
 	VecDoub b(m+1);
 	for (ipj=0;ipj<=(m << 1);ipj++) {
diff -rupN code/series.h code302/series.h
--- code/series.h	Thu May 10 14:24:13 2007
+++ code302/series.h	Tue Sep 16 17:25:12 2008
@@ -106,4 +106,3 @@ struct Levin {
 		return (lastval = val);
 	}
 };
-
diff -rupN code/simplex.h code302/simplex.h
--- code/simplex.h	Thu May 10 14:24:13 2007
+++ code302/simplex.h	Tue Sep 16 17:25:12 2008
@@ -161,16 +161,25 @@ void Simplex::solve()
 	if (verbose)
 		cout << "    at end of phase0,iter= " << nsteps << endl;
 	if (ierr != 0) {
+		delete lu;
+		for (Int i=0;i<n;i++)
+			delete a[i+1];
 		return;
 	}
 	phase1();
 	if (verbose)
 		cout << "    at end of phase1,iter= " << nsteps << endl;
 	if (ierr != 0) {
+		delete lu;
+		for (Int i=0;i<n;i++)
+			delete a[i+1];
 		return;
 	}
 	phase2();
 	prepare_output();
+	delete lu;
+	for (Int i=0;i<n;i++)
+		delete a[i+1];
 }
 
 void Simplex::initialize() {
@@ -290,7 +299,7 @@ void Simplex::phase1()
 				v[i]=1.0/scale[ord[i]];
 				done=false;
 			}
-  			else
+			else
 				v[i]=0.0;
 		}
 		if (done)
@@ -331,7 +340,7 @@ void Simplex::phase2()
 		for (Int i=1;i<=m;i++) {
 			if (ord[i] > 0)
 				v[i]=-obj[ord[i]];
-  			else
+			else
 				v[i]=-obj[ord[i]+nm1];
 		}
 		kp=colpiv(v,2,piv);
diff -rupN code/stiel.h code302/stiel.h
--- code/stiel.h	Thu May 10 14:24:13 2007
+++ code302/stiel.h	Tue Sep 16 17:25:12 2008
@@ -108,4 +108,3 @@ void Stiel::get_weights(VecDoub_O &x, Ve
 	amu0=b[0];
 	gaucof(a,b,amu0,x,w);
 }
-
diff -rupN code/svd.h code302/svd.h
--- code/svd.h	Thu May 10 14:24:13 2007
+++ code302/svd.h	Tue Sep 16 17:25:12 2008
@@ -49,13 +49,13 @@ void SVD::solve(VecDoub_I &b, VecDoub_O 
 
 void SVD::solve(MatDoub_I &b, MatDoub_O &x, Doub thresh = -1.)
 {
-	int i,j,m=b.ncols();
-	if (b.nrows() != n || x.nrows() != n || b.ncols() != x.ncols())
+	int i,j,p=b.ncols();
+	if (b.nrows() != m || x.nrows() != n || x.ncols() != p)
 		throw("SVD::solve bad sizes");
-	VecDoub xx(n);
-	for (j=0;j<m;j++) {
-		for (i=0;i<n;i++) xx[i] = b[i][j];
-		solve(xx,xx,thresh);
+	VecDoub xx(n),bcol(m);
+	for (j=0;j<p;j++) {
+		for (i=0;i<m;i++) bcol[i] = b[i][j];
+		solve(bcol,xx,thresh);
 		for (i=0;i<n;i++) x[i][j] = xx[i];
 	}
 }