Mercurial > hg > jgplsrc

comparison ar.c @ 0:e0bbaa717f41 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
lol J
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Mon, 25 Nov 2013 11:56:30 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:e0bbaa717f41
1 /* Copyright 1990-2011, Jsoftware Inc. All rights reserved. */
2 /* License in license.txt. */
3 /* */
4 /* Adverbs: Reduce (Insert) and Outer Product */
5
6 #include "j.h"
7 #include "vasm.h"
8 #include "ve.h"
9 #include "vcomp.h"
10 #include "ar.h"
11
12 static DF1(jtreduce);
13
14
15 #define PARITY2 u=(UC*)&s; b=0; b^=*u++; b^=*u++;
16 #define PARITY4 u=(UC*)&s; b=0; b^=*u++; b^=*u++; b^=*u++; b^=*u++;
17 #define PARITY8 u=(UC*)&s; b=0; b^=*u++; b^=*u++; b^=*u++; b^=*u++; b^=*u++; b^=*u++; b^=*u++; b^=*u++;
18
19 #if SY_64
20 #define PARITYW PARITY8
21 #else
22 #define PARITYW PARITY4
23 #endif
24
25 #if SY_ALIGN
26 #define VDONE(T,PAR) \
27 {I q=n/sizeof(T);T s,*y=(T*)x; DO(m, s=0; DO(q, s^=*y++;); PAR; *z++=b==pc;);}
28
29 static void vdone(I m,I n,B*x,B*z,B pc){B b,*u;
30 if(1==m){UI s,*xi;
31 s=0; b=0;
32 xi=(I*)x; DO(n/SZI, s^=*xi++;);
33 u=(B*)xi; DO(n%SZI, b^=*u++;);
34 u=(B*)&s; DO(SZI, b^=*u++;);
35 *z=b==pc;
36 }else if(0==n%sizeof(UI ))VDONE(UI, PARITYW)
37 else if(0==n%sizeof(UINT))VDONE(UINT,PARITY4)
38 else if(0==n%sizeof(US ))VDONE(US, PARITY2)
39 else DO(m, b=0; DO(n, b^=*x++;); *z++=b==pc;);
40 }
41 #else
42 static void vdone(I m,I n,B*x,B*z,B pc){B b;I q,r;UC*u;UI s,*y;
43 q=n/SZI; r=n%SZI; y=(UI*)x;
44 switch((r?2:0)+pc){
45 case 0: DO(m, s=0; DO(q, s^=*y++;); PARITYW; *z++=!b;); break;
46 case 1: DO(m, s=0; DO(q, s^=*y++;); PARITYW; *z++= b;); break;
47 case 2: DO(m, s=0; DO(q, s^=*y++;); PARITYW; u=(UC*)y; DO(r, b^=*u++;); *z++=!b; x+=n; y=(UI*)x;); break;
48 case 3: DO(m, s=0; DO(q, s^=*y++;); PARITYW; u=(UC*)y; DO(r, b^=*u++;); *z++= b; x+=n; y=(UI*)x;); break;
49 }}
50 #endif
51
52 #define RBFXLOOP(T,pfx) \
53 {T*xx=(T*)x,*yy,*z0,*zz=(T*)z; \
54 q=d/sizeof(T); \
55 for(j=0;j<m;++j){ \
56 yy=xx; xx-=q; z0=zz; DO(q, --xx; --yy; --zz; *zz=pfx(*xx,*yy);); \
57 DO(n-2, zz=z0; DO(q, --xx; --zz; *zz=pfx(*xx,*zz););); \
58 }} /* non-commutative */
59
60 #define RCFXLOOP(T,pfx) \
61 {T*xx=(T*)x,*yy,*z0,*zz=(T*)z; \
62 q=d/sizeof(T); \
63 for(j=0;j<m;++j){ \
64 yy=xx; xx+=q; z0=zz; DO(q, *zz++=pfx(*yy,*xx); ++xx; ++yy;); \
65 DO(n-2, zz=z0; DO(q, *zz++=pfx(*zz,*xx); ++xx; );); \
66 }} /* commutative */
67
68 #if SY_ALIGN
69 #define RBFXODDSIZE(pfx,bpfx) RBFXLOOP(C,bpfx)
70 #define REDUCECFX REDUCEBFX
71 #else
72 #define RBFXODDSIZE(pfx,bpfx) \
73 {B*zz;I r,t,*xi,*yi,*zi; \
74 q=d/SZI; r=d%SZI; xi=(I*)x; zz=z; \
75 for(j=0;j<m;++j,zz-=d){ \
76 yi=xi; xi=(I*)((B*)xi-d); zi=(I*)zz; \
77 DO(q, --xi; --yi; *--zi=pfx(*xi,*yi);); \
78 xi=(I*)((B*)xi-r); yi=(I*)((B*)yi-r); t=pfx(*xi,*yi); MC((B*)zi-r,&t,r); \
79 DO(n-2, zi=(I*)zz; DO(q, --xi; --zi; *zi=pfx(*xi,*zi);); \
80 xi=(I*)((B*)xi-r); zi=(I*)((B*)zi-r); t=pfx(*xi,*zi); MC( zi, &t,r);); \
81 }} /* non-commutative */
82
83 #define RCFXODDSIZE(pfx,bpfx) \
84 {I r,t,*xi,*yi,*zi; \
85 q=d/SZI; r=d%SZI; \
86 for(j=0;j<m;++j,x+=d,z+=d){ \
87 yi=(I*)x; x+=d; xi=(I*)x; zi=(I*)z; DO(q, *zi++=pfx(*yi,*xi); ++xi; ++yi;); t=pfx(*yi,*xi); MC(zi,&t,r); \
88 DO(n-2, x+=d; xi=(I*)x; zi=(I*)z; DO(q, *zi++=pfx(*zi,*xi); ++xi; ); t=pfx(*zi,*xi); MC(zi,&t,r);); \
89 }} /* commutative */
90
91 #define REDUCECFX(f,pfx,ipfx,spfx,bpfx,vdo) \
92 AHDRP(f,B,B){B*y=0;I d,j,q; \
93 if(c==n){vdo; R;} \
94 d=c/n; \
95 if(1==n)DO(d, *z++=*x++;) \
96 else if(0==d%sizeof(UI ))RCFXLOOP(UI, pfx) \
97 else if(0==d%sizeof(UINT))RCFXLOOP(UINT,ipfx) \
98 else if(0==d%sizeof(US ))RCFXLOOP(US, spfx) \
99 else RCFXODDSIZE(pfx,bpfx) \
100 } /* commutative */
101
102 #endif
103
104 #define REDUCEBFX(f,pfx,ipfx,spfx,bpfx,vdo) \
105 AHDRP(f,B,B){B*y=0;I d,j,q; \
106 if(c==n){vdo; R;} \
107 d=c/n; x+=m*c; z+=m*d; \
108 if(1==n)DO(d, *--z=*--x;) \
109 else if(0==d%sizeof(UI ))RBFXLOOP(UI, pfx) \
110 else if(0==d%sizeof(UINT))RBFXLOOP(UINT,ipfx) \
111 else if(0==d%sizeof(US ))RBFXLOOP(US, spfx) \
112 else RBFXODDSIZE(pfx,bpfx) \
113 } /* non-commutative */
114
115 REDUCECFX( eqinsB, EQ, IEQ, SEQ, BEQ, vdone(m,n,x,z,(B)(n%2)))
116 REDUCECFX( neinsB, NE, INE, SNE, BNE, vdone(m,n,x,z,1 ))
117 REDUCECFX( orinsB, OR, IOR, SOR, BOR, DO(m, *z++=1&&memchr(x,C1,n); x+=c;))
118 REDUCECFX( andinsB, AND, IAND, SAND, BAND, DO(m, *z++=! memchr(x,C0,n); x+=c;))
119 REDUCEBFX( ltinsB, LT, ILT, SLT, BLT, DO(m, *z++= *(x+n-1)&&!memchr(x,C1,n-1)?1:0; x+=c;))
120 REDUCEBFX( leinsB, LE, ILE, SLE, BLE, DO(m, *z++=!*(x+n-1)&&!memchr(x,C0,n-1)?0:1; x+=c;))
121 REDUCEBFX( gtinsB, GT, IGT, SGT, BGT, DO(m, y=memchr(x,C0,n); *z++=1&&(y?1&(y-x):1&n); x+=c;))
122 REDUCEBFX( geinsB, GE, IGE, SGE, BGE, DO(m, y=memchr(x,C1,n); *z++=! (y?1&(y-x):1&n); x+=c;))
123 REDUCEBFX( norinsB, NOR, INOR, SNOR, BNOR, DO(m, y=memchr(x,C1,n); d=y?y-x:n; *z++=(1&d)==d<n-1; x+=c;))
124 REDUCEBFX(nandinsB, NAND,INAND,SNAND,BNAND,DO(m, y=memchr(x,C0,n); d=y?y-x:n; *z++=(1&d)!=d<n-1; x+=c;))
125
126
127 #if SY_ALIGN
128 REDUCEPFX(plusinsB,I,B,PLUS)
129 #else
130 AHDRR(plusinsB,I,B){I d,dw,i,p,q,r,r1,s;UC*tu;UI*v;
131 if(c==n&&n<SZI)DO(m, s=0; DO(n, s+=*x++;); *z++=s;)
132 else if(c==n){UI t;
133 p=n/SZI; q=p/255; r=p%255; r1=n%SZI; tu=(UC*)&t;
134 for(i=0;i<m;++i){
135 s=0; v=(UI*)x;
136 DO(q, t=0; DO(255, t+=*v++;); DO(SZI, s+=tu[i];));
137 t=0; DO(r, t+=*v++;); DO(SZI, s+=tu[i];);
138 x=(B*)v; DO(r1, s+=*x++;);
139 *z++=s;
140 }}else{A t;UI*tv;
141 d=c/n; dw=(d+SZI-1)/SZI; p=dw*SZI; memset(z,C0,m*d*SZI);
142 q=n/255; r=n%255;
143 t=ga(INT,dw,1,0); if(!t)R;
144 tu=(UC*)AV(t); tv=(UI*)tu; v=(UI*)x;
145 for(i=0;i<m;++i,z+=d){
146 DO(q, memset(tv,C0,p); DO(255, DO(dw,tv[i]+=v[i];); x+=d; v=(UI*)x;); DO(d,z[i]+=tu[i];));
147 memset(tv,C0,p); DO(r, DO(dw,tv[i]+=v[i];); x+=d; v=(UI*)x;); DO(d,z[i]+=tu[i];) ;
148 }}} /* +/"r w on boolean w, originally by Roger Moore */
149 #endif
150
151
152 REDUCEOVF( plusinsI, I, I, PLUSR, PLUSVV, PLUSRV)
153 REDUCEOVF(minusinsI, I, I, MINUSR,MINUSVV,MINUSRV)
154 REDUCEOVF(tymesinsI, I, I, TYMESR,TYMESVV,TYMESRV)
155
156 REDUCCPFX( plusinsO, D, I, PLUSO)
157 REDUCCPFX(minusinsO, D, I, MINUSO)
158 REDUCCPFX(tymesinsO, D, I, TYMESO)
159
160 REDUCENAN( plusinsD, D, D, PLUS )
161 REDUCENAN( plusinsZ, Z, Z, zplus )
162 REDUCEPFX( plusinsX, X, X, xplus )
163
164 REDUCEPFX(minusinsB, I, B, MINUS )
165 REDUCENAN(minusinsD, D, D, MINUS )
166 REDUCENAN(minusinsZ, Z, Z, zminus)
167
168 REDUCEPFX(tymesinsD, D, D, TYMES )
169 REDUCEPFX(tymesinsZ, Z, Z, ztymes)
170
171 REDUCENAN( divinsD, D, D, DIV )
172 REDUCENAN( divinsZ, Z, Z, zdiv )
173
174 REDUCEPFX( maxinsI, I, I, MAX )
175 REDUCEPFX( maxinsD, D, D, MAX )
176 REDUCEPFX( maxinsX, X, X, XMAX )
177 REDUCEPFX( maxinsS, SB,SB,SBMAX )
178
179 REDUCEPFX( mininsI, I, I, MIN )
180 REDUCEPFX( mininsD, D, D, MIN )
181 REDUCEPFX( mininsX, X, X, XMIN )
182 REDUCEPFX( mininsS, SB,SB,SBMIN )
183
184
185 static DF1(jtred0){DECLF;A x;I f,r,wr,*s;
186 wr=AR(w); r=jt->rank?jt->rank[1]:wr; f=wr-r; s=AS(w);
187 jt->rank=0;
188 GA(x,AT(w),0L,r,f+s);
189 R reitem(vec(INT,f,s),lamin1(df1(x,iden(fs))));
190 } /* f/"r w identity case */
191
192 static DF1(jtredg){PROLOG;DECLF;A y,z;B p;C*u,*v;I i,k,n,old,r,wr,yn,yr,*ys,yt;
193 RZ(w);
194 ASSERT(DENSE&AT(w),EVNONCE);
195 wr=AR(w); r=jt->rank?jt->rank[1]:wr; jt->rank=0;
196 if(r<wr)R rank1ex(w,self,r,jtredg);
197 n=IC(w); p=ARELATIVE(w);
198 RZ(z=tail(w)); yt=AT(z); yn=AN(z); yr=AR(z); ys=1+AS(w);
199 k=yn*bp(yt); v=CAV(w)+k*(n-1);
200 old=jt->tbase+jt->ttop;
201 for(i=1;i<n;++i){
202 v-=k;
203 GA(y,yt,yn,yr,ys); u=CAV(y);
204 if(p){A1*wv=(A1*)v,*yv=(A1*)u;I d=(I)w-(I)y; AFLAG(y)=AFREL; DO(yn, yv[i]=d+wv[i];);}else MC(u,v,k);
205 RZ(z=CALL2(f2,y,z,fs));
206 gc(z,old);
207 }
208 EPILOG(z);
209 } /* f/"r w for general f and 1<(-r){$w */
210
211
212 static C fca[]={CSTAR, CPLUS, CEQ, CMIN, CMAX, CPLUSDOT, CSTARDOT, CNE, 0}; /* commutative & associative */
213
214 static A jtredsp1a(J jt,C id,A z,A e,I n,I r,I*s){A t;B b,p=0;D d=1;
215 switch(id){
216 default: ASSERT(0,EVNONCE);
217 case CPLUSDOT: R n?gcd(z,e):ca(e);
218 case CSTARDOT: R n?lcm(z,e):ca(e);
219 case CMIN: R n?minimum(z,e):ca(e);
220 case CMAX: R n?maximum(z,e):ca(e);
221 case CPLUS: if(n&&equ(e,zero))R z; DO(r, d*=s[i];); t=tymes(e,d>IMAX?scf(d-n):sc((I)d-n)); R n?plus (z,t):t;
222 case CSTAR: if(n&&equ(e,one ))R z; DO(r, d*=s[i];); t=expn2(e,d>IMAX?scf(d-n):sc((I)d-n)); R n?tymes(z,t):t;
223 case CEQ: p=1; /* fall thru */
224 case CNE:
225 ASSERT(B01&AT(e),EVNONCE);
226 if(!n)*BAV(z)=p;
227 b=1; DO(r, if(!(s[i]%2)){b=0; break;});
228 R !p==*BAV(e)&&b!=n%2?not(z):z;
229 }} /* f/w on sparse vector w, post processing */
230
231 static A jtredsp1(J jt,A w,A self,C id,VF ado,I cv,I f,I r,I zt){A e,x,z;I m,n;P*wp;
232 RZ(w);
233 wp=PAV(w); e=SPA(wp,e); x=SPA(wp,x); n=AN(x); m=*AS(w);
234 GA(z,zt,1,0,0);
235 if(n){ado(jt,1L,n,n,AV(z),AV(x)); RE(0); if(m==n)R z;}
236 R redsp1a(id,z,e,n,AR(w),AS(w));
237 } /* f/"r w for sparse vector w */
238
239 DF1(jtredravel){A f,x,z;C id;I cv,n;P*wp;VF ado;
240 RZ(w);
241 f=VAV(self)->f;
242 if(!(SPARSE&AT(w)))R reduce(AN(w)?gah(1L,w):mtv,f);
243 wp=PAV(w); x=SPA(wp,x); n=AN(x);
244 id=vaid(VAV(f)->f);
245 while(1){
246 vains(id,AT(x),&ado,&cv);
247 ASSERT(ado,EVNONCE);
248 GA(z,rtype(cv),1,0,0);
249 if(n)ado(jt,1L,n,n,AV(z),AV(x));
250 if(jt->jerr!=EWOV)R redsp1a(id,z,SPA(wp,e),n,AR(w),AS(w));;
251 }} /* f/@, w */
252
253 static A jtredspd(J jt,A w,A self,C id,VF ado,I cv,I f,I r,I zt){A a,e,x,z,zx;I c,m,n,*s,t,*v,wr,*ws,xf,xr;P*wp,*zp;
254 RZ(w);
255 ASSERT(strchr(fca,id),EVNONCE);
256 wp=PAV(w); a=SPA(wp,a); e=SPA(wp,e); x=SPA(wp,x); s=AS(x);
257 xr=r; v=AV(a); DO(AN(a), if(f<v[i])--xr;); xf=AR(x)-xr;
258 m=prod(xf,s); c=m?AN(x)/m:0; n=s[xf];
259 GA(zx,zt,AN(x)/n,AR(x)-1,s); ICPY(xf+AS(zx),1+xf+s,xr-1);
260 ado(jt,m,c,n,AV(zx),AV(x)); RE(0);
261 switch(id){
262 case CPLUS: if(!equ(e,zero))RZ(e=tymes(e,sc(n))); break;
263 case CSTAR: if(!equ(e,one )&&!equ(e,zero))RZ(e=expn2(e,sc(n))); break;
264 case CEQ: ASSERT(B01&AT(x),EVNONCE); if(!*BAV(e)&&0==n%2)e=one; break;
265 case CNE: ASSERT(B01&AT(x),EVNONCE); if( *BAV(e)&&1==n%2)e=zero;
266 }
267 if(AT(e)!=AT(zx)){t=maxtype(AT(e),AT(zx)); if(t!=AT(zx))RZ(zx=cvt(t,zx));}
268 wr=AR(w); ws=AS(w);
269 GA(z,STYPE(AT(zx)),1,wr-1,ws); if(1<wr)ICPY(f+AS(z),f+1+ws,wr-1);
270 zp=PAV(z);
271 RZ(a=ca(a)); v=AV(a); DO(AN(a), if(f<v[i])--v[i];);
272 SPB(zp,a,a);
273 SPB(zp,e,cvt(AT(zx),e));
274 SPB(zp,i,SPA(wp,i));
275 SPB(zp,x,zx);
276 R z;
277 } /* f/"r w for sparse w, rank > 1, dense axis */
278
279 static B jtredspsprep(J jt,C id,I f,I zt,A a,A e,A x,A y,I*zm,I**zdv,B**zpv,I**zqv,C**zxxv,A*zsn){
280 A d,p,q,sn=0,xx;B*pv;C*xxv;I*dv,j,k,m,mm,*qv=0,*u,*v,yc,yr,yr1,*yv;
281 v=AS(y); yr=v[0]; yc=v[1]; yr1=yr-1;
282 RZ(d=grade1(eq(a,sc(f)))); dv=AV(d);
283 DO(AN(a), if(i!=dv[i]){RZ(q=grade1p(d,y)); qv=AV(q); break;});
284 GA(p,B01,yr,1,0); pv=BAV(p); memset(pv,C0,yr);
285 u=yv=AV(y); m=mm=0; j=-1; if(qv)v=yv+yc*qv[0];
286 for(k=0;k<yr1;++k){
287 if(qv){u=v; v=yv+yc*qv[1+k];}else v=u+yc;
288 DO(yc-1, if(u[dv[i]]!=v[dv[i]]){++m; pv[k]=1; mm=MAX(mm,k-j); j=k; break;});
289 if(!qv)u+=yc;
290 }
291 if(yr){++m; pv[yr1]=1; mm=MAX(mm,yr1-j);}
292 if(qv){j=mm*aii(x); GA(xx,AT(x),j,1,0); xxv=CAV(xx);}
293 switch(id){
294 case CPLUS: case CPLUSDOT: j=!equ(e,zero); break;
295 case CSTAR: case CSTARDOT: j=!equ(e,one); break;
296 case CMIN: j=!equ(e,zt&B01?one :zt&INT?sc(IMAX):ainf ); break;
297 case CMAX: j=!equ(e,zt&B01?zero:zt&INT?sc(IMIN):scf(infm)); break;
298 case CEQ: j=!*BAV(e); break;
299 case CNE: j= *BAV(e); break;
300 }
301 if(j)GA(sn,INT,m,1,0);
302 *zm=m; *zdv=dv; *zpv=pv; *zqv=qv; *zxxv=xxv; *zsn=sn;
303 R 1;
304 }
305
306 static B jtredspse(J jt,C id,I wm,I xt,A e,A zx,A sn,A*ze,A*zzx){A b;B nz;I t,zt;
307 RZ(b=ne(zero,sn)); nz=!all0(b); zt=AT(zx);
308 switch(id){
309 case CPLUS: if(nz)RZ(zx=plus (zx, tymes(e,sn) )); RZ(e= tymes(e,sc(wm)) ); break;
310 case CSTAR: if(nz)RZ(zx=tymes(zx,bcvt(1,expn2(e,sn)))); RZ(e=bcvt(1,expn2(e,sc(wm)))); break;
311 case CPLUSDOT: if(nz)RZ(zx=gcd(zx,from(b,over(zero,e)))); break;
312 case CSTARDOT: if(nz)RZ(zx=lcm(zx,from(b,over(one ,e)))); break;
313 case CMIN: if(nz)RZ(zx=minimum(zx,from(b,over(zt&B01?one: zt&INT?sc(IMAX):ainf, e)))); break;
314 case CMAX: if(nz)RZ(zx=maximum(zx,from(b,over(zt&B01?zero:zt&INT?sc(IMIN):scf(infm),e)))); break;
315 case CEQ: ASSERT(B01&xt,EVNONCE); if(nz)RZ(zx=eq(zx,eq(zero,residue(num[2],sn)))); if(!(wm%2))e=one; break;
316 case CNE: ASSERT(B01&xt,EVNONCE); if(nz)RZ(zx=ne(zx,eq(one, residue(num[2],sn)))); if(!(wm%2))e=zero; break;
317 }
318 if(AT(e)!=AT(zx)){t=maxtype(AT(e),AT(zx)); if(t!=AT(zx))RZ(zx=cvt(t,zx));}
319 *ze=e; *zzx=zx;
320 R 1;
321 }
322
323 static A jtredsps(J jt,A w,A self,C id,VF ado,I cv,I f,I r,I zt){A a,a1,e,sn,x,x1=0,y,z,zx,zy;B*pv;
324 C*xv,*xxv,*zv;I*dv,i,m,n,*qv,*sv,*v,wr,xk,xt,wm,*ws,xc,yc,yr,*yu,*yv,zk;P*wp,*zp;
325 RZ(w);
326 ASSERT(strchr(fca,id),EVNONCE);
327 wr=AR(w); ws=AS(w); wm=ws[f];
328 wp=PAV(w); a=SPA(wp,a); e=SPA(wp,e);
329 y=SPA(wp,i); v=AS(y); yr=v[0]; yc=v[1]; yv=AV(y);
330 x=SPA(wp,x); xt=AT(x); xc=aii(x);
331 RZ(redspsprep(id,f,zt,a,e,x,y,&m,&dv,&pv,&qv,&xxv,&sn));
332 xv=CAV(x); xk=xc*bp(xt);
333 GA(zx,zt,m*xc,AR(x),AS(x)); *AS(zx)=m; zv=CAV(zx); zk=xc*bp(zt);
334 GA(zy,INT,m*(yc-1),2,0); v=AS(zy); v[0]=m; v[1]=yc-1; yu=AV(zy);
335 v=qv; if(sn)sv=AV(sn);
336 for(i=0;i<m;++i){A y;B*p1;C*u;I*vv;
337 p1=1+(B*)memchr(pv,C1,yr); n=p1-pv; if(sn)sv[i]=wm-n; pv=p1;
338 vv=qv?yv+yc**v:yv; DO(yc-1, *yu++=vv[dv[i]];);
339 if(1<n){if(qv){u=xxv; DO(n, MC(u,xv+xk*v[i],xk); u+=xk;);} ado(jt,1L,n*xc,n,zv,qv?xxv:xv); RE(0);}
340 else if(zk==xk)MC(zv,qv?xv+xk**v:xv,xk);
341 else {if(!x1)GA(x1,xt,xc,1,0); MC(AV(x1),qv?xv+xk**v:xv,xk); RZ(y=cvt(zt,x1)); MC(zv,AV(y),zk);}
342 zv+=zk; if(qv)v+=n; else{xv+=n*xk; yv+=n*yc;}
343 }
344 if(sn)RZ(redspse(id,wm,xt,e,zx,sn,&e,&zx));
345 RZ(a1=ca(a)); v=AV(a1); n=0; DO(AN(a), if(f!=v[i])v[n++]=v[i]-(f<v[i]););
346 GA(z,STYPE(AT(zx)),1,wr-1,ws); if(1<r)ICPY(f+AS(z),f+1+ws,r-1);
347 zp=PAV(z);
348 SPB(zp,a,vec(INT,n,v));
349 SPB(zp,e,cvt(AT(zx),e));
350 SPB(zp,x,zx);
351 SPB(zp,i,zy);
352 R z;
353 } /* f/"r w for sparse w, rank > 1, sparse axis */
354
355 static DF1(jtreducesp){A a,g,x,y,z;B b;C id;I cv,f,n,r,rr[2],*v,wn,wr,*ws,wt,zt;P*wp;VF ado;
356 RZ(w);
357 wr=AR(w); r=jt->rank?jt->rank[1]:wr; f=wr-r;
358 wn=AN(w); ws=AS(w); n=r?ws[f]:1;
359 wt=AT(w); wt=wn?DTYPE(wt):B01;
360 g=VAV(self)->f; id=vaid(g);
361 if(!n)R red0(w,self);
362 vains(id,wt,&ado,&cv);
363 if(2==n&&!(ado&&strchr(fca,id))){
364 rr[0]=0; rr[1]=r;
365 jt->rank=rr; x=from(zero,w);
366 jt->rank=rr; y=from(one, w);
367 R df2(x,y,g);
368 }
369 if(!ado)R redg(w,self);
370 if(1==n)R tail(w);
371 zt=rtype(cv);
372 jt->rank=0;
373 if(1==wr)z=redsp1(w,self,id,ado,cv,f,r,zt);
374 else{
375 wp=PAV(w); a=SPA(wp,a); v=AV(a);
376 b=0; DO(AN(a), if(f==v[i]){b=1; break;});
377 z=b?redsps(w,self,id,ado,cv,f,r,zt):redspd(w,self,id,ado,cv,f,r,zt);
378 }
379 R jt->jerr==EWOV?(rr[1]=r,jt->rank=rr,reducesp(w,self)):z;
380 } /* f/"r for sparse w */
381
382 #define BR2IFX(T,F) {T*u=(T*)wv,*v=u+d,x,y; \
383 GA(z,B01,wn/2,wr-1,ws); zv=BAV(z); \
384 if(1<d)DO(m, DO(d, x=*u++; y=*v++; *zv++=x F y; ); u+=d; v+=d;) \
385 else DO(m, x=*u++; y=*u++; *zv++=x F y; ); \
386 }
387 #define BR2PFX(T,F) {T*u=(T*)wv,*v=u+d,x,y; \
388 GA(z,B01,wn/2,wr-1,ws); zv=BAV(z); \
389 if(1<d)DO(m, DO(d, x=*u++; y=*v++; *zv++=F(x,y);); u+=d; v+=d;) \
390 else DO(m, x=*u++; y=*u++; *zv++=F(x,y); ); \
391 }
392 #define BTABIFX(F) {btab[0 ]=0 F 0; \
393 btab[SYS&SYS_LILENDIAN?256: 1]=0 F 1; \
394 btab[SYS&SYS_LILENDIAN? 1:256]=1 F 0; \
395 btab[257 ]=1 F 1; \
396 }
397 #define BTABPFX(F) {btab[0 ]=F(0,0); \
398 btab[SYS&SYS_LILENDIAN?256: 1]=F(0,1); \
399 btab[SYS&SYS_LILENDIAN? 1:256]=F(1,0); \
400 btab[257 ]=F(1,1); \
401 }
402 #define BR2CASE(t,id) ((id)+256*(t))
403
404 static B jtreduce2(J jt,A w,C id,I f,I r,A*zz){A z=0;B b=0,btab[258],*zv;I c,d,m,wn,wr,*ws,*wv;
405 wn=AN(w); wr=AR(w); ws=AS(w); wv=AV(w);
406 m=prod(f,ws); c=m?wn/m:prod(r,f+ws); d=c/2;
407 switch(BR2CASE(AT(w),id)){
408 case BR2CASE(B01,CEQ ): if(b=1==r)BTABIFX(== ); break;
409 case BR2CASE(B01,CNE ): if(b=1==r)BTABIFX(!= ); break;
410 case BR2CASE(B01,CLT ): if(b=1==r)BTABIFX(< ); break;
411 case BR2CASE(B01,CLE ): if(b=1==r)BTABIFX(<= ); break;
412 case BR2CASE(B01,CGT ): if(b=1==r)BTABIFX(> ); break;
413 case BR2CASE(B01,CGE ): if(b=1==r)BTABIFX(>= ); break;
414 case BR2CASE(B01,CMAX ):
415 case BR2CASE(B01,CPLUSDOT): if(b=1==r)BTABIFX(|| ); break;
416 case BR2CASE(B01,CPLUSCO ): if(b=1==r)BTABPFX(BNOR ); break;
417 case BR2CASE(B01,CMIN ):
418 case BR2CASE(B01,CSTAR ):
419 case BR2CASE(B01,CSTARDOT): if(b=1==r)BTABIFX(&& ); break;
420 case BR2CASE(B01,CSTARCO ): if(b=1==r)BTABPFX(BNAND); break;
421 case BR2CASE(LIT,CEQ ): BR2IFX(C,== ); break;
422 case BR2CASE(LIT,CNE ): BR2IFX(C,!= ); break;
423 case BR2CASE(INT,CEQ ): BR2IFX(I,== ); break;
424 case BR2CASE(INT,CLT ): BR2IFX(I,< ); break;
425 case BR2CASE(INT,CLE ): BR2IFX(I,<= ); break;
426 case BR2CASE(INT,CGT ): BR2IFX(I,> ); break;
427 case BR2CASE(INT,CGE ): BR2IFX(I,>= ); break;
428 case BR2CASE(INT,CNE ): BR2IFX(I,!= ); break;
429 case BR2CASE(FL, CEQ ): BR2PFX(D,TEQ); break;
430 case BR2CASE(FL, CLT ): BR2PFX(D,TLT); break;
431 case BR2CASE(FL, CLE ): BR2PFX(D,TLE); break;
432 case BR2CASE(FL, CGT ): BR2PFX(D,TGT); break;
433 case BR2CASE(FL, CGE ): BR2PFX(D,TGE); break;
434 case BR2CASE(FL, CNE ): BR2PFX(D,TNE); break;
435 }
436 if(b){S*u=(S*)wv; GA(z,B01,wn/2,wr-1,ws); zv=BAV(z); DO(m, *zv++=btab[*u++];);}
437 if(z&&1<r){I*u=f+AS(z),*v=f+1+ws; DO(r-1, *u++=*v++;);}
438 *zz=z;
439 R 1;
440 } /* f/"r for dense w over an axis of length 2 */
441
442 static DF1(jtreduce){A z;C id;I c,cv,f,m,n,r,rr[2],t,wn,wr,*ws,wt,zt;VF ado;
443 RZ(w);
444 wr=AR(w); r=jt->rank?jt->rank[1]:wr; f=wr-r;
445 wn=AN(w); ws=AS(w); n=r?ws[f]:1;
446 if(SPARSE&AT(w))R reducesp(w,self);
447 wt=AT(w); wt=wn?wt:B01;
448 id=vaid(VAV(self)->f);
449 switch(n){
450 case 0: R red0(w,self);
451 case 1: R tail(w);
452 case 2: RZ(reduce2(w,id,f,r,&z)); if(z)R z;
453 }
454 vains(id,wt,&ado,&cv);
455 if(!ado)R redg(w,self);
456 zt=rtype(cv); jt->rank=0;
457 GA(z,zt,wn/n,MAX(0,wr-1),ws); if(1<r)ICPY(f+AS(z),f+1+ws,r-1);
458 if((t=atype(cv))&&t!=wt)RZ(w=cvt(t,w));
459 m=prod(f,ws); c=m?wn/m:prod(r,f+ws);
460 ado(jt,m,c,n,AV(z),AV(w));
461 if(jt->jerr)R jt->jerr==EWOV?(rr[1]=r,jt->rank=rr,reduce(w,self)):0; else R cv&VRI+VRD?cvz(cv,z):z;
462 } /* f/"r w main control */
463
464 static A jtredcatsp(J jt,A w,A z,I r){A a,q,x,y;B*b;I c,d,e,f,j,k,m,n,n1,p,*u,*v,wr,*ws,xr;P*wp,*zp;
465 ws=AS(w); wr=AR(w); f=wr-r; p=ws[1+f];
466 wp=PAV(w); x=SPA(wp,x); y=SPA(wp,i); a=SPA(wp,a); v=AV(a);
467 m=*AS(y); n=AN(a); n1=n-1; xr=AR(x);
468 RZ(b=bfi(wr,a,1));
469 c=b[f]; d=b[1+f]; if(c&&d)b[f]=0; e=f+!c;
470 j=0; DO(n, if(e==v[i]){j=i; break;});
471 k=1; DO(f, if(!b[i])++k;);
472 zp=PAV(z); SPB(zp,e,ca(SPA(wp,e)));
473 GA(q,INT,n-(c&&d),1,0); v=AV(q); DO(wr, if(b[i])*v++=i-(i>f);); SPB(zp,a,q);
474 if(c&&d){ /* sparse sparse */
475 SPB(zp,x,ca(x));
476 SPB(zp,i,q=repeatr(ne(a,sc(f)),y));
477 v=j+AV(q); u=j+AV(y);
478 DO(m, *v=p*u[0]+u[1]; v+=n1; u+=n;);
479 }else if(!c&&!d){ /* dense dense */
480 u=AS(x); GA(q,AT(x),AN(x),xr-1,u); v=k+AS(q); *v=u[k]*u[1+k]; ICPY(1+v,2+k+u,xr-k-2);
481 MC(AV(q),AV(x),AN(x)*bp(AT(x)));
482 SPB(zp,x,q); SPB(zp,i,ca(y));
483 }else{ /* other */
484 GA(q,INT,xr,1,0); v=AV(q);
485 if(1!=k){*v++=0; *v++=k; e=0; DO(xr-1, ++e; if(e!=k)*v++=e;); RZ(x=cant2(q,x));}
486 v=AV(q); u=AS(x); *v=u[0]*u[1]; ICPY(1+v,2+u,xr-1); RZ(x=reshape(vec(INT,xr-1,v),x));
487 e=ws[f+c]; RZ(y=repeat(sc(e),y)); v=j+AV(y);
488 if(c)DO(m, k=p**v; DO(e, *v=k+ i; v+=n;);)
489 else DO(m, k= *v; DO(e, *v=k+p*i; v+=n;););
490 RZ(q=grade1(y)); RZ(y=from(q,y)); RZ(x=from(q,x));
491 SPB(zp,x,x); SPB(zp,i,y);
492 }
493 R z;
494 } /* ,/"r w for sparse w, 2<r */
495
496 static DF1(jtredcat){A z;B b;I f,r,*s,*v,wr;
497 RZ(w);
498 wr=AR(w); r=jt->rank?jt->rank[1]:wr; f=wr-r; s=AS(w); jt->rank=0;
499 b=1==r&&1==s[f];
500 if(2>r&&!b)R ca(w);
501 GA(z,AT(w),AN(w),wr-1,s);
502 if(!b){v=f+AS(z); RE(*v=mult(s[f],s[1+f])); ICPY(1+v,2+f+s,r-2);}
503 if(SPARSE&AT(w))R redcatsp(w,z,r);
504 MC(AV(z),AV(w),AN(w)*bp(AT(w)));
505 if(ARELATIVE(w)){AFLAG(z)=AFREL; z=relocate((I)w-(I)z,z);}
506 R z;
507 } /* ,/"r w */
508
509 static DF1(jtredsemi){I f,n,r,*s,wr;
510 RZ(w);
511 wr=AR(w); r=jt->rank?jt->rank[1]:wr; f=wr-r; s=AS(w); n=r?s[f]:1;
512 if(2>n){ASSERT(n,EVDOMAIN); R tail(w);}
513 if(BOX&AT(w))R irs2(rank1ex(curtail(w),0L,r-1,jtbox),tail(w),0L,r,r-1,jtover);
514 else R irs1(w,0L,r-1,jtbox);
515 } /* ;/"r w */
516
517 static DF1(jtredstitch){A c,y;I f,n,r,*s,*v,wr;
518 RZ(w);
519 wr=AR(w); r=jt->rank?jt->rank[1]:wr; f=wr-r; jt->rank=0;
520 s=AS(w); n=r?s[f]:1;
521 ASSERT(n,EVDOMAIN);
522 if(1==n)R irs1(w,0L,r,jthead);
523 if(1==r)R 2==n?ca(w):irs2(irs2(num[-2],w,0L,0L,1L,jtdrop),irs2(num[-2],w,0L,0L,1L,jttake),0L,1L,0L,jtover);
524 if(2==r)R irs1(w,0L,2L,jtcant1);
525 RZ(c=IX(wr)); v=AV(c); v[f]=f+1; v[f+1]=f; RZ(y=cant2(c,w));
526 if(SPARSE&AT(w)){A x;
527 GA(x,INT,f+r-1,1,0); v=AV(x); ICPY(v,AS(y),f+1);
528 RE(v[f+1]=mult(s[f],s[f+2])); ICPY(v+f+2,s+3+f,r-3);
529 R reshape(x,y);
530 }else{
531 v=AS(y);
532 RE(v[f+1]=mult(s[f],s[f+2])); ICPY(v+f+2,s+3+f,r-3);
533 --AR(y);
534 R y;
535 }} /* ,./"r w */
536
537 static DF1(jtredstiteach){A*wv,y;I n,p,r,t,wd;
538 RZ(w);
539 n=AN(w);
540 if(!(2<n&&1==AR(w)&&BOX&AT(w)))R reduce(w,self);
541 wv=AAV(w); wd=(I)w*ARELATIVE(w); y=WVR(0); p=IC(y); t=AT(y);
542 DO(n, y=WVR(i); r=AR(y); if(!(r&&r<=2&&p==IC(y)&&t==AT(y)))R reduce(w,self););
543 R box(razeh(w));
544 } /* ,.&.>/ w */
545
546 static DF1(jtredcateach){A*u,*v,*wv,x,*xv,z,*zv;I f,m,mn,n,r,wd,wr,*ws,zm,zn;
547 RZ(w);
548 wr=AR(w); ws=AS(w); r=jt->rank?jt->rank[1]:wr; f=wr-r; jt->rank=0;
549 n=r?ws[f]:1;
550 if(!r||1>=n)R reshape(repeat(ne(sc(f),IX(wr)),shape(w)),n?w:ace);
551 if(!(BOX&AT(w)))R df1(cant2(sc(f),w),qq(ds(CBOX),one));
552 zn=AN(w)/n; zm=prod(f,ws); m=zm?AN(w)/(zm*n):prod(r-1,ws+f+1); mn=m*n;
553 GA(z,BOX,zn,wr-1,ws); ICPY(AS(z)+f,ws+f+1,r-1);
554 GA(x,BOX,n,1,0); xv=AAV(x);
555 zv=AAV(z); wv=AAV(w); wd=(I)w*ARELATIVE(w);
556 DO(zm, u=wv; DO(m, v=u++; DO(n, xv[i]=AADR(wd,*v); v+=m;); RZ(*zv++=raze(x));); wv+=mn;);
557 R z;
558 } /* ,&.>/"r w */
559
560 static DF2(jtoprod){R df2(a,w,VAV(self)->h);}
561
562
563 F1(jtslash){A h;AF f1=jtreduce;C c;V*v;
564 RZ(w);
565 if(NOUN&AT(w))R evger(w,sc(GINSERT));
566 v=VAV(w);
567 switch(v->id){
568 case CCOMMA: f1=jtredcat; break;
569 case CCOMDOT: f1=jtredstitch; break;
570 case CSEMICO: f1=jtredsemi; break;
571 case CUNDER: if(COPE==ID(v->g)){c=ID(v->f); if(c==CCOMMA)f1=jtredcateach; else if(c==CCOMDOT)f1=jtredstiteach;}
572 }
573 RZ(h=qq(w,v2(lr(w),RMAX)));
574 R fdef(CSLASH,VERB, f1,jtoprod, w,0L,h, 0L, RMAX,RMAX,RMAX);
575 }
576
577 A jtaslash (J jt,C c, A w){RZ( w); R df1( w, slash(ds(c)) );}
578 A jtaslash1(J jt,C c, A w){RZ( w); R df1( w,qq(slash(ds(c)),one));}
579 A jtatab (J jt,C c,A a,A w){RZ(a&&w); R df2(a,w, slash(ds(c)) );}
580
581
582 static AHDRR(jtmeanD,D,D){I d,i;D*y;D v,*zz;
583 d=c/n;
584 NAN0;
585 if(1==d)DO(m, v= *x++; DO(n-1, v+=*x++;); *z++=v/n;)
586 else for(i=0;i<m;++i){
587 y=x; x+=d; zz=z; DO(d, *z++ =*x+++ *y++;);
588 DO(n-3, z=zz; DO(d, *z+++=*x++; ));
589 z=zz; DO(d, *z =(*z+*x++)/n; ++z;);
590 }
591 NAN1V;
592 } /* based on REDUCEPFX; 2<n */
593
594 static AHDRR(jtmeanI,D,I){I d,i;I*y;D v,*zz;
595 d=c/n;
596 if(1==d)DO(m, v=(D)*x++; DO(n-1, v+=*x++;); *z++=v/n;)
597 else for(i=0;i<m;++i){
598 y=x; x+=d; zz=z; DO(d, *z++ =*x+++(D)*y++;);
599 DO(n-3, z=zz; DO(d, *z+++=*x++; ));
600 z=zz; DO(d, *z =(*z+*x++)/n; ++z;);
601 }} /* based on REDUCEPFX; 2<n */
602
603 DF1(jtmean){A z;I c,f,m,n,r,wn,wr,*ws,wt;
604 RZ(w);
605 wr=AR(w); r=jt->rank?jt->rank[1]:wr; f=wr-r; jt->rank=0;
606 wt=AT(w); wn=AN(w); ws=AS(w); n=r?ws[f]:1;
607 if(!(wn&&2<n&&wt&INT+FL))R divide(df1(w,qq(slash(ds(CPLUS)),sc(r))),sc(n));
608 GA(z,FL,wn/n,MAX(0,wr-1),ws); if(1<r)ICPY(f+AS(z),f+1+ws,r-1);
609 m=prod(f,ws); c=m?wn/m:prod(r,f+ws);
610 if(wt&INT)meanI(m,c,n,DAV(z), AV(w));
611 else meanD(m,c,n,DAV(z),DAV(w));
612 RE(0); R z;
613 } /* (+/%#)"r w */