Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Tree: 4368a1a84d
Fetching contributors…

Cannot retrieve contributors at this time

276 lines (257 sloc) 10.555 kB
/* Copyright 1990-2011, Jsoftware Inc. All rights reserved. */
/* License in license.txt. */
/* */
/* Verbs: a#"r w */
#include "j.h"
#define REPF(f) A f(J jt,A a,A w,I wf,I wcr)
static REPF(jtrepzdx){A p,q,x;P*wp;
RZ(a&&w);
if(SPARSE&AT(w)){wp=PAV(w); x=SPA(wp,e);}
else x=jt->fill&&AN(jt->fill)?jt->fill:filler(w);
RZ(p=repeat(ravel(rect(a)),ravel(stitch(IX(wcr?*(wf+AS(w)):1),num[-1]))));
RZ(q=irs2(w,x,0L,wcr,0L,jtover));
R irs2(p,q,0L,1L,wcr+!wcr,jtfrom);
} /* (dense complex) # (dense or sparse) */
static REPF(jtrepzsx){A q,x,y;I c,d,j,k=-1,m,p=0,*qv,*xv,*yv;P*ap;
RZ(a&&w);
ap=PAV(a); x=SPA(ap,x); m=AN(x);
if(!AN(SPA(ap,a)))R repzdx(ravel(x),w,wf,wcr);
y=SPA(ap,i); yv=AV(y);
RZ(x=cvt(INT,vec(FL,2*m,AV(x)))); xv=AV(x);
if(equ(zero,SPA(ap,e))){
k=c=*(wf+AS(w));
if(!wf&&SPARSE&AT(w)){A a,y;I m,n,q,*v;P*wp;
wp=PAV(w); a=SPA(wp,a);
if(AN(a)&&!*AV(a)){
y=SPA(wp,i); v=AS(y); m=v[0]; n=v[1]; v=AV(y);
k=m?v[(m-1)*n]+1:0; q=0;
DO(m, if(q==*v)++q; else if(q<*v){k=q; break;} v+=n;);
}}
ASSERT(k<=IMAX-1,EVLIMIT);
if(c==k)RZ(w=irs2(sc(1+k),w,0L,0L,wcr,jttake));
DO(2*m, ASSERT(0<=xv[i],EVDOMAIN); p+=xv[i]; ASSERT(0<=p,EVLIMIT););
GA(q,INT,p,1,0); qv=AV(q);
DO(m, c=*xv++; d=*xv++; j=yv[i]; DO(c, *qv++=j;); DO(d, *qv++=k;););
R irs2(q,w,0L,1L,wcr,jtfrom);
}
ASSERT(0,EVNONCE);
} /* (sparse complex) #"r (dense or sparse) */
#define REPB(T) \
{T*u,*v=(T*)zv; \
for(i=0;i<c;++i){ \
u=i*m+(T*)wv; \
for(j=0,iv=(I*)b;j<q;++j,u+=SZI)switch(*iv++){ \
case B0001: *v++=u[3]; break; \
case B0010: *v++=u[2]; break; \
case B0011: *v++=u[2]; *v++=u[3]; break; \
case B0100: *v++=u[1]; break; \
case B0101: *v++=u[1]; *v++=u[3]; break; \
case B0110: *v++=u[1]; *v++=u[2]; break; \
case B0111: *v++=u[1]; *v++=u[2]; *v++=u[3]; break; \
case B1000: *v++=u[0]; break; \
case B1001: *v++=u[0]; *v++=u[3]; break; \
case B1010: *v++=u[0]; *v++=u[2]; break; \
case B1011: *v++=u[0]; *v++=u[2]; *v++=u[3]; break; \
case B1100: *v++=u[0]; *v++=u[1]; break; \
case B1101: *v++=u[0]; *v++=u[1]; *v++=u[3]; break; \
case B1110: *v++=u[0]; *v++=u[1]; *v++=u[2]; break; \
case B1111: *v++=u[0]; *v++=u[1]; *v++=u[2]; *v++=u[3]; \
} \
if(r){B*c=(B*)iv; DO(r, if(c[i])*v++=u[i];);} \
}}
#if !SY_64 && SY_WIN32
static REPF(jtrepbdx){A z;B*b;C*wv,*zv;I c,i,*iv,j,k,m,p,q,r,zn;
RZ(a&&w);
if(SPARSE&AT(w))R irs2(ifb(AN(a),BAV(a)),w,0L,1L,wcr,jtfrom);
m=AN(a); q=m/SZI; r=m%SZI;
ASSERT(m==*(wf+AS(w)),EVLENGTH);
b=BAV(a); p=bsum(m,b); zn=m?p*(AN(w)/m):0;
ASSERT(0<=zn,EVLIMIT);
GA(z,AT(w),zn,AR(w),AS(w)); *(wf+AS(z))=p;
wv=CAV(w); zv=CAV(z);
RE(c=prod(wf,AS(w)));
if(zn)switch(k=AN(w)/(c*m)*bp(AT(w)),FL&AT(w)||k!=sizeof(D)?k:0){
case sizeof(C): REPB(C); break;
case sizeof(S): REPB(S); break;
case sizeof(I): REPB(I); break;
case sizeof(D): REPB(D); break;
default: {C*u;I k1=k,k2=k*2,k3=k*3,k4=k*4,km=k*m;
for(i=0;i<c;++i){
u=i*km+wv;
for(j=0,iv=(I*)b;j<q;++j,u+=k4)switch(*iv++){
case B0001: MC(zv,k3+u,k1); zv+=k1; break;
case B0010: MC(zv,k2+u,k1); zv+=k1; break;
case B0011: MC(zv,k2+u,k2); zv+=k2; break;
case B0100: MC(zv,k1+u,k1); zv+=k1; break;
case B0101: MC(zv,k1+u,k1); zv+=k1; MC(zv,k3+u,k1); zv+=k1; break;
case B0110: MC(zv,k1+u,k2); zv+=k2; break;
case B0111: MC(zv,k1+u,k3); zv+=k3; break;
case B1000: MC(zv, u,k1); zv+=k1; break;
case B1001: MC(zv, u,k1); zv+=k1; MC(zv,k3+u,k1); zv+=k1; break;
case B1010: MC(zv, u,k1); zv+=k1; MC(zv,k2+u,k1); zv+=k1; break;
case B1011: MC(zv, u,k1); zv+=k1; MC(zv,k2+u,k2); zv+=k2; break;
case B1100: MC(zv, u,k2); zv+=k2; break;
case B1101: MC(zv, u,k2); zv+=k2; MC(zv,k3+u,k1); zv+=k1; break;
case B1110: MC(zv, u,k3); zv+=k3; break;
case B1111: MC(zv, u,k4); zv+=k4;
}
if(r){B*c=(B*)iv; DO(r, if(c[i]){MC(zv,u+i*k,k); zv+=k;});}
}}}
R RELOCATE(w,z);
} /* (dense boolean)#"r (dense or sparse) */
#else
static REPF(jtrepbdx){A z;B*b;C*wv,*zv;I c,k,m,p,zn;
RZ(a&&w);
if(SPARSE&AT(w))R irs2(ifb(AN(a),BAV(a)),w,0L,1L,wcr,jtfrom);
m=AN(a);
b=BAV(a); p=bsum(m,b); zn=m?p*(AN(w)/m):0;
ASSERT(0<=zn,EVLIMIT);
GA(z,AT(w),zn,AR(w),AS(w)); *(wf+AS(z))=p;
wv=CAV(w); zv=CAV(z);
if(!zn)R z;
RE(c=prod(wf,AS(w))); k=AN(w)/(c*m)*bp(AT(w));
DO(c, DO(m, if(b[i]){MC(zv,wv,k); zv+=k;} wv+=k;););
R RELOCATE(w,z);
} /* (dense boolean)#"r (dense or sparse) */
#endif
static REPF(jtrepbsx){A ai,c,d,e,g,q,x,wa,wx,wy,y,y1,z,zy;B*b;I*dv,*gv,j,m,n,*u,*v,*v0;P*ap,*wp,*zp;
RZ(a&&w);
ap=PAV(a); e=SPA(ap,e);
y=SPA(ap,i); u=AV(y);
x=SPA(ap,x); n=AN(x); b=BAV(x);
if(!AN(SPA(ap,a)))R irs2(ifb(n,b),w,0L,1L,wcr,jtfrom);
if(!*BAV(e)){
GA(q,INT,n,1,0); v=v0=AV(q);
DO(n, if(*b++)*v++=u[i];);
AN(q)=*AS(q)=v-v0;
R irs2(q,w,0L,1L,wcr,jtfrom);
}
wp=PAV(w);
if(DENSE&AT(w)||all0(eq(sc(wf),SPA(wp,a)))){RZ(q=denseit(a)); R irs2(ifb(AN(q),BAV(q)),w,0L,1L,wcr,jtfrom);}
wa=SPA(wp,a); wy=SPA(wp,i); wx=SPA(wp,x);
RZ(q=aslash(CPLUS,a));
GA(z,AT(w),1,AR(w),AS(w)); *(wf+AS(z))=m=*AV(q);
RZ(c=indexof(wa,sc(wf)));
RZ(y1=fromr(c,wy));
RZ(q=not(eps(y1,ravel(repeat(not(x),y)))));
m=*AS(a)-m;
GA(ai,INT,m,1,0); v=AV(ai); DO(n, if(!*b++)*v++=u[i];);
RZ(g=grade1(over(ai,repeat(q,y1)))); gv=AV(g);
GA(d,INT,AN(y1),1,0); dv=AV(d); j=0; DO(AN(g), if(m>gv[i])++j; else dv[gv[i]-m]=j;);
RZ(zy=repeat(q,wy)); v=AV(zy)+*AV(c); m=*(1+AS(zy)); DO(*AS(zy), *v-=dv[i]; v+=m;);
zp=PAV(z);
SPB(zp,a,ca(wa));
SPB(zp,e,SPA(wp,e));
SPB(zp,i,zy);
SPB(zp,x,repeat(q,wx));
R z;
} /* (sparse boolean) #"r (dense or sparse) */
static REPF(jtrepidx){A y;I j,m,p=0,*v,*x;
RZ(a&&w);
RZ(a=vi(a)); x=AV(a);
m=*AS(a);
DO(m, ASSERT(0<=x[i],EVDOMAIN); p+=x[i]; ASSERT(0<=p,EVLIMIT););
GA(y,INT,p,1,0); v=AV(y);
DO(m, j=i; DO(x[j], *v++=j;););
R irs2(y,w,0L,1L,wcr,jtfrom);
} /* (dense integer) #"r (dense or sparse) */
static REPF(jtrepisx){A e,q,x,y;I c,j,m,p=0,*qv,*xv,*yv;P*ap;
RZ(a&&w);
ap=PAV(a); e=SPA(ap,e);
y=SPA(ap,i); yv=AV(y);
x=SPA(ap,x); if(!(INT&AT(x)))RZ(x=cvt(INT,x)); xv=AV(x);
if(!AN(SPA(ap,a)))R repidx(ravel(x),w,wf,wcr);
if(!*AV(e)){
m=AN(x);
DO(m, ASSERT(0<=xv[i],EVDOMAIN); p+=xv[i]; ASSERT(0<=p,EVLIMIT););
GA(q,INT,p,1,0); qv=AV(q);
DO(m, c=xv[i]; j=yv[i]; DO(c, *qv++=j;););
R irs2(q,w,0L,1L,wcr,jtfrom);
}
ASSERT(0,EVNONCE);
} /* (sparse integer) #"r (dense or sparse) */
static REPF(jtrep1d){A z;C*wv,*zv;I c,k,m,n,p=0,q,t,*ws,zk,zn;
RZ(a&&w);
t=AT(a); m=AN(a); ws=AS(w); n=wcr?ws[wf]:1;
if(t&CMPX){
if(wcr)R repzdx(from(apv(n,0L,0L),a),w, wf,wcr);
else R repzdx(a,irs2(apv(m,0L,0L),w,0L,1L,0L,jtfrom),wf,1L );
}
if(t&B01){B*x=BAV(a); DO(m,p+=x[i];);}
else{I*x;
RZ(a=vi(a)); x=AV(a);
DO(m, ASSERT(0<=x[i],EVDOMAIN); p+=x[i]; ASSERT(0<=p,EVLIMIT););
}
RE(q=mult(p,n));
RE(zn=n?mult(q,AN(w)/n):0);
GA(z,AT(w),zn,AR(w)+!wcr,ws); *(wf+AS(z))=q;
if(!zn)R z;
wv=CAV(w); zv=CAV(z);
RE(c=prod(wf,ws)); k=AN(w)/(c*n)*bp(AT(w)); zk=p*k;
DO(c*n, mvc(zk,zv,k,wv); zv+=zk; wv+=k;);
R RELOCATE(w,z);
} /* scalar #"r dense or dense #"0 dense */
static B jtrep1sa(J jt,A a,I*c,I*d){A x;B b;I*v;
b=1&&AT(a)&CMPX;
if(b)RZ(x=rect(a)) else x=a;
if(AR(a)){ASSERT(equ(one,aslash(CSTARDOT,le(zero,ravel(x)))),EVDOMAIN); RZ(x=aslash(CPLUS,x));}
if(!(INT&AT(x)))RZ(x=cvt(INT,x));
v=AV(x); *c=v[0]; *d=b?v[1]:0;
ASSERT(0<=*c&&0<=*d,EVDOMAIN);
R 1;
} /* process a in a#"0 w */
static REPF(jtrep1s){A ax,e,x,y,z;B*b;I c,d,cd,j,k,m,n,p,q,*u,*v,wr,*ws;P*wp,*zp;
RZ(a&&w);
if(AT(a)&SCMPX)R rep1d(denseit(a),w,wf,wcr);
RE(rep1sa(a,&c,&d)); cd=c+d;
if(DENSE&AT(w))R rep1d(d?jdot2(sc(c),sc(d)):sc(c),w,wf,wcr);
wr=AR(w); ws=AS(w); n=wcr?*(wf+ws):1; RE(m=mult(n,cd));
wp=PAV(w); e=SPA(wp,e); ax=SPA(wp,a); y=SPA(wp,i); x=SPA(wp,x);
GA(z,AT(w),1,wr+!wcr,ws); *(wf+AS(z))=m; zp=PAV(z);
RE(b=bfi(wr,ax,1));
if(wcr&&b[wf]){ /* along sparse axis */
u=AS(y); p=u[0]; q=u[1]; u=AV(y);
RZ(x=repeat(sc(c),x));
RZ(y=repeat(sc(c),y));
if(p&&1<c){
j=0; DO(wf, j+=b[i];); v=j+AV(y);
if(AN(ax)==1+j){u+=j; DO(p, m=cd**u; u+=q; DO(c, *v=m+i; v+=q;););}
else{A xx;I h,i,j1=1+j,*uu;
GA(xx,INT,j1,1,0); uu=AV(xx);
k=0; DO(j1, uu[i]=u[i];);
for(i=0;i<p;++i,u+=q)
if(ICMP(uu,u,j1)||i==p-1){
h=(i==p-1)+i-k; k=i; m=cd*uu[j];
DO(j1, uu[i]=u[i];);
DO(h, DO(c, *v=m+i; v+=q;););
}
RZ(xx=grade1(y));
RZ(x=from(xx,x));
RZ(y=from(xx,y));
}}}else{A xx; /* along dense axis */
j=0; DO(wcr, j+=!b[wf+i];);
RZ(y=ca(y));
if(d){xx=jt->fill; jt->fill=e;}
x=irs2(AR(a)&&CMPX&AT(a)?a:d?jdot2(sc(c),sc(d)):sc(c),x,0L,1L,j,jtrepeat);
if(d)jt->fill=xx;
RZ(x);
}
SPB(zp,e,e);
SPB(zp,a,ax);
SPB(zp,i,y);
SPB(zp,x,x);
R z;
} /* scalar #"r sparse or sparse #"0 (dense or sparse) */
F2(jtrepeat){B ab,wb;I acr,ar,at,m,wcr,wf,wr,wt;
RZ(a&&w);
ar=AR(a); acr=jt->rank?jt->rank[0]:ar;
wr=AR(w); wcr=jt->rank?jt->rank[1]:wr; wf=wr-wcr; jt->rank=0;
at=AT(a); ab=1&&at&DENSE;
wt=AT(w); wb=1&&wt&DENSE;
if(1<acr||acr<ar)R rank2ex(a,w,0L,MIN(1,acr),wcr,jtrepeat);
ASSERT(!acr||!wcr||(m=*AS(a),m==*(wf+AS(w))),EVLENGTH);
if(!acr||!wcr)R ab&&wb?rep1d(a,w,wf,wcr):rep1s(a,w,wf,wcr);
if(at&CMPX+SCMPX)R ab?repzdx(a,w,wf,wcr):repzsx(a,w,wf,wcr);
if(at&B01 +SB01 )R ab?repbdx(a,w,wf,wcr):repbsx(a,w,wf,wcr);
/* integer */ R ab?repidx(a,w,wf,wcr):repisx(a,w,wf,wcr);
} /* a#"r w main control */
Jump to Line
Something went wrong with that request. Please try again.