orthog(M,T,n)/* routine returns in M an orthogonal matrix Oand in T an upper triangular matrix T such that M=O*T *//* if M non-singular, T non-singular */register double *M,*T;register int n;{	register int i,j,k;	register double *s,*c,*sT,*s1T,*sM;	double sqrt(),dot();	for(k=0,s1T=T,c=M;k<n;k++,s1T+=n+1,c++)		{		if(k>0)			{			for(j=0,sT=T+k,sM=M;j<k;j++,sT+=n,sM++)				{				*sT=dot(sM,n,c,n,n);				sp1(*sT,sM,c,n);				}			}		(*s1T)=dot(c,n,c,n,n);		if((*s1T)<=0)			return(1);		*s1T=sqrt(*s1T);		for(i=0,s=c;i<n;i++,s+=n)			*s/= *s1T;		}	return(0);}sp1(p1,p2,p3,n)register double p1,*p2,*p3;register int n;{	register int i;	for(i=0;i<n;i++,p2+=n,p3+=n)		*p3 -= p1 * *p2;	return;}