new file mode 100644
--- /dev/null
+++ b/vgauss.c
@@ -0,0 +1,113 @@
+/* Copyright 1990-2011, Jsoftware Inc.  All rights reserved. */
+/* License in license.txt.                                   */
+/*                                                                         */
+/* Verbs: Gaussian Elimination                                             */
+
+#include "j.h"
+
+                    
+F1(jtgausselm){A t;C*tv;I c,e,i,j,m,old,r,r1,*s;Q p,*u,*v,*x;
+ F1RANK(2,jtgausselm,0);
+ ASSERT(RAT&AT(w),EVNONCE);
+ ASSERT(2==AR(w),EVRANK);
+ s=AS(w); r=s[0]; c=s[1]; r1=MIN(r,c); v=QAV(w); m=c*bp(AT(w));
+ GA(t,LIT,m,1,0); tv=CAV(t);
+ old=jt->tbase+jt->ttop;
+ for(j=0;j<r1;++j){
+  e=-1; u=v+c*j+j; DO(r-j, if(XDIG(u->n)){e=i+j; break;} u+=c;);  /* find pivot row */
+  ASSERT(0<=e,EVDOMAIN);
+  x=v+c*j; 
+  if(j!=e){u=v+c*e; MC(tv,x,m); MC(x,u,m); MC(u,tv,m);} /* interchange rows e and j */
+  p=x[j]; DO(c, x[i]=qdiv(x[i],p););
+  for(i=0;i<r;++i){
+   if(i==j)continue;
+   u=v+c*i; p=u[j];  /* pivot */
+   DO(c, u[i]=qminus(u[i],qtymes(p,x[i])););
+  }
+  gc(w,old);
+ }
+ R w;
+}    /* Gaussian elimination in place */
+
+static F1(jtdetr){A t,z;C*tv;I c,e,g=1,i,j,k,m,old,r,*s;Q d,p,*u,*v,*x;
+ RZ(w);
+ s=AS(w); r=s[0]; c=s[1];
+ v=QAV(w); 
+ m=c*sizeof(Q); GA(t,LIT,m,1,0); tv=CAV(t);
+ old=jt->tbase+jt->ttop;
+ for(j=0;j<r;++j){
+  e=-1; u=v+c*j+j; DO(r-j, if(XDIG(u->n)){e=i+j; break;} u+=c;);  /* find pivot row */
+  if(0>e)R cvt(RAT,zero);
+  x=v+c*j;
+  if(j!=e){u=v+c*e; MC(tv,x,m); MC(x,u,m); MC(u,tv,m); g=-g;}  /* interchange rows e and j */
+  i=XDIG(x[j].n); if(i==XPINF||i==XNINF)R mark;
+  for(i=j+1;i<r;++i){
+   u=v+c*i;
+   if(XDIG(u[j].n)){p=qdiv(u[j],x[j]); for(k=j+1;k<r;++k)u[k]=qminus(u[k],qtymes(p,x[k]));}
+  }
+  gc(w,old);
+ }
+ d=0<g?*v:qminus(zeroQ,*v); u=v+1+c; DO(r-1, d=qtymes(d,*u); u+=1+c;);
+ RE(0);
+ GA(z,RAT,1,0,0); *QAV(z)=d; R z;
+}    /* determinant on rational matrix; works in place */
+
+static F1(jtdetd){D g,h,p,q,*u,*v,*x,*y,z=1.0;I c,d,e,i,j,k,r,*s;
+ RZ(w);
+ s=AS(w); r=s[0]; c=s[1]; v=DAV(w);
+ NAN0;
+ for(j=0;j<r;++j){
+  x=v+c*j; u=x+j; h=0.0; 
+  DO(r-j, k=i; DO(c-j, g=ABS(*u); if(h<g){h=g; d=j+k; e=j+i;} ++u;); u+=j;);  /* find pivot, maximum abs element */
+  if(h==inf)R mark;
+  if(0==h)R scf(0.0);
+  if(j!=d){u=v+c*d+j; y=x+j; DO(c-j, q=*u; *u=*y; *y=q; ++u;  ++y; ); z=-z;}  /* interchange rows j and d */
+  if(j!=e){u=x+e;     y=x+j; DO(r-j, q=*u; *u=*y; *y=q; u+=c; y+=c;); z=-z;}  /* interchange cols j and e */
+  q=x[j]; z*=q; JBREAK0;
+  for(i=j+1;i<r;++i){
+   u=v+c*i;
+   if(u[j]){p=u[j]/q; for(k=j+1;k<r;++k)u[k]-=p*x[k];}
+ }}
+ NAN1;
+ R scf(z);
+}    /* determinant on real     matrix; works in place */
+
+#define ZABT(v)         ((v).re*(v).re+(v).im*(v).im)
+
+static F1(jtdetz){A t;D g,h;I c,d,e,i,j,k,r,*s;Z p,q,*u,*v,*x,*y,z;
+ RZ(w);
+ z.re=1.0; z.im=0.0;
+ s=AS(w); r=s[0]; c=s[1]; v=ZAV(w);
+ NAN0;
+ for(j=0;j<r;++j){
+  x=v+c*j; u=x+j; h=0.0; 
+  DO(r-j, k=i; DO(c-j, g=ZABT(*u); if(h<g){h=g; d=j+k; e=j+i;} ++u;); u+=j;);  /* find pivot, maximum abs element */
+  if(h==inf)R mark;
+  if(0==h)R scf(0.0);
+  if(j!=d){u=v+c*d;              DO(c-j, q=u[j+i]; u[j+i]=x[j+i]; x[j+i]=q;); z=zminus(zeroZ,z);}  /* interchange rows j and d */
+  if(j!=e){u=v+c*j+e; y=v+c*j+j; DO(r-j, q=*u; *u=*y; *y=q; u+=c; y+=c;);     z=zminus(zeroZ,z);}  /* interchange cols j and e */
+  q=x[j]; z=ztymes(z,q);
+  for(i=j+1;i<r;++i){
+   u=v+c*i;
+   if(ZNZ(u[j])){p=zdiv(u[j],q); for(k=j+1;k<r;++k)u[k]=zminus(u[k],ztymes(p,x[k]));}
+ }}
+ NAN1; RE(0);
+ GA(t,CMPX,1,0,0); *ZAV(t)=z; R t;
+}    /* determinant on complex  matrix; works in place */
+
+F1(jtgaussdet){A z;I*s;
+ RZ(w);
+ ASSERT(2==AR(w),EVRANK);
+ s=AS(w);
+ ASSERT(s[0]==s[1],EVLENGTH);
+ switch(AT(w)){
+  default:   ASSERT(0,EVDOMAIN);
+  case B01:
+  case INT:  R detd(cvt(FL,w));
+  case FL:   z=detd(ca(w));      break;
+  case CMPX: z=detz(ca(w));      break;
+  case XNUM: z=detr(cvt(RAT,w)); break;
+  case RAT:  z=detr(ca(w));
+ }
+ R z==mark?detxm(w,eval("-/ .*")):z;
+}    /* determinant on square matrix */