Fix +/\ +/\. (+/%#)/\ with NaNs

Author: HenryHRich

jsrc/ap.c

@@ -734,17 +734,22 @@ static DF2(jtinfixd){F12IP;A z;C*x,*y;I c=0,d,k,m,n,p,q,r,*s,wr,*ws,wt,zc;
 #define SETZ    {s=yv; DQ(c, *zv++=*s++;  );}
 #define SETZD   {s=yv; DQ(c, *zv++=*s++/d;);}
 
+// check for infinity/NaN, EVNAN if so
+#define CKINFB(addr)
+#define CKINFI(addr)
+#define CKINFD(addr) ASSERT(((US*)(addr))[3]&0x7ffc!=0x7ffc,EVNAN)
+
 #define MOVSUMAVG(Tw,Ty,ty,Tz,tz,xd,SET)  \
  {Tw*u,*v;Ty*s,x=0,*yv;Tz*zv;                                  \
   GATVS(z,tz,c*(1+p),AR(w),AS(w),tz##SIZE,GACOPYSHAPE,R 0); AS(z)[0]=1+p;                    \
   zv=(Tz*)AV(z); u=v=(Tw*)AV(w);                               \
   if(1==c){                                                    \
-   DQ(m, x+=*v++;); *zv++=xd;                                  \
-   DQ(p, x+=(Ty)*v++-(Ty)*u++; *zv++=xd;);                     \
+   DQ(m, CKINF##Tw(v) x+=*v++;); *zv++=xd;                                  \
+   DQ(p, CKINF##Tw(v) x+=(Ty)*v++-(Ty)*u++; *zv++=xd;);                     \
   }else{                                                       \
    GATVS(y,ty,c,1,0,ty##SIZE,GACOPYSHAPE0,R 0); s=yv=(Ty*)AV1(y); DQ(c, *s++=0;);            \
-   DQ(m, s=yv; DQ(c, *s+++=*v++;);); SET;                      \
-   DQ(p, s=yv; DQ(c, x=*s+++=(Ty)*v++-(Ty)*u++; *zv++=xd;););  \
+   DQ(m, s=yv; DQ(c, CKINF##Tw(v) *s+++=*v++;);); SET;                      \
+   DQ(p, s=yv; DQ(c, CKINF##Tw(v) x=*s+++=(Ty)*v++-(Ty)*u++; *zv++=xd;););  \
  }}
 
 static A jtmovsumavg1(J jt,I m,A w,A fs,B avg){A y,z;D d=(D)m;I c,p,wt;

jsrc/as.c

@@ -354,8 +354,9 @@ static DF2(jtofxinv){F12IP;A f,fs,z;C c;I t;V*v;
  F2RANKW(0,RMAX,jtofxinv,self);
  fs=FAV(self)->fgh[0]; f=FAV(fs)->fgh[0]; v=FAV(f); c=v->id; t=AT(w);  // self = f/\. fs = f/  f = f  v = verb info for f
  if(!(c==CPLUS||c==CBDOT&&t&INT||((c&-2)==CEQ)&&t&B01))R outfix(a,w,self);  // if not +/\. or m b./\. or =/\. or ~:/\.
- A z0,z1; z=irs2(dfv1(z0,w,fs),dfv2(z1,a,w,bslash(fs)),c==CPLUS?ds(CMINUS):f, RMAX,-1L,jtatomic2);
- if(jt->jerr==EVNAN){RESETERR; R outfix(a,w,self);}else R z;
+ A z0; RZ(dfv1(z0,w,fs)) if(t&FL&&unlikely(!all0(jtisnan(jt,z0)))){RESETERR; R outfix(a,w,self);}  // take u/ y; if FL, see if result contains NaN, failover if so to slow version
+ A z1; z=irs2(z0,dfv2(z1,a,w,bslash(fs)),c==CPLUS?ds(CMINUS):f, RMAX,-1L,jtatomic2);  // x u/\. y = (u/ y) uinv x u/\ y
+ if(jt->jerr==EVNAN){RESETERR; R outfix(a,w,self);}else R z;  // failover if NaN
 }    /* a f/\. w where f has an "undo" */
 
 static DF2(jtofxassoc){F12IP;A f,i,j,p,s,x,z;C id,*zv;I c,d,k,kc,m,r,t;V*v;VA2 adocv;

jsrc/j.h

@@ -2084,6 +2084,7 @@ if(likely(type _i<3)){z=(type _i<1)?1:(type _i==1)?_zzt[0]:_zzt[0]*_zzt[1];}else
 // PRODX takes the product of init and v[0..n-1], generating error if overflow, but waiting till the end so no error if there is a 0 in the product
 // overflow sets z to the error value of 0; if we see a multiplicand of 0 we stop right away so we can skip the error
 // This is written to be branchless for rank < 3
+// n may be negative
 #if SY_64
 // I have been unable to make clang produce a simple loop that doesn't end with a backward branch.  So I am going to handle ranks 0-2 here and call a subroutine for the rest
 #define PRODX(z,n,v,init) \

jsrc/ja.h

@@ -920,6 +920,7 @@ extern void jfree4gmp(void*,size_t);
 #define raincr(x) __atomic_fetch_add(&AC(x),1,__ATOMIC_ACQ_REL)
 #define rasv(x)   {I c=AC(x); if(likely(!ACISPERM(c))){if(c<0)AC(x)=(I)((UI)c+(ACINPLACE+ACUC1));else raincr(x); radescend(x,sv)}}  // better a misbranch than an atomic instruction if c<0.  Could avoid recur check if AC>1
 #define ra(x)   {I c=AC(x); if(likely(!ACISPERM(c))){if(c<0)AC(x)=(I)((UI)c+(ACINPLACE+ACUC1));else raincr(x); radescend(x)}}  // better a misbranch than an atomic instruction if c<0.  Could avoid recur check if AC>1
+#define raN(x,n)   {I c=AC(x); if(likely(!ACISPERM(c))){if(c<0)AC(x)=(I)((UI)c+(ACINPLACE+(n)));else __atomic_fetch_add(&AC(x),(n),__ATOMIC_ACQ_REL); radescend(x)}}  // add N to usecount all at once
 #define racontents(x)   {I c=AC(x); if(MEMAUDIT!=0&&c<0)SEGFAULT; if(likely(!ACISPERM(c))){raincr(x); radescend(x)}}  // Used on contents of box, which cannot have AC<0
 #define rareccontents(x)   {I c=AC(x); if(MEMAUDIT!=0&&c<0)SEGFAULT; if(likely(!ACISPERM(c))){raincr(x);}}  // Used on contents of recursive box, which cannot have AC<0 and does not need recursion
 #define raname(x) {if(likely(!ACISPERM(AC(x))))raincr(x);}  // NAME is not inplaceable, seldom local; just add 1.  No traverse needed on ra

jsrc/vcat.c

@@ -22,7 +22,7 @@ F2(jtsetfv){F12IP;A q=jt->fill;
   jt->fillv=CAV(q); jt->fillvlen=bpnoun(t);  // jt->fillv points to the fill atom
  }else{fillv0(t);}    // empty or no fill.  create std fill in fillv0 and point jt->fillv at it
  w=TYPESEQ(t,AT(w))?w:cvt(t,w);  // convert w to dest type.  note if w is boxed and nonempty this won't change it
- if(w==0&&jt->jerr==EVDOMAIN)jt->jerr=EVINHOMO; // if we got an error here (always called DOMAIN), show it as EVHOMO when we eformat
+ if(unlikely(w==0)&&jt->jerr==EVDOMAIN)jt->jerr=EVINHOMO; // if we got an error here (always called DOMAIN), show it as EVHOMO when we eformat
  R w;
 }
 // simpler version when there is only w.  We make AT(w) an argument to save a fetch

jsrc/vf.c

@@ -10,16 +10,20 @@ F1(jtcurtail){F12IP; R jtdrop(jtfg,num(-1),w);}
 
 F1(jtshift1){F12IP;R drop(num(-1),over(num(1),w));}
 
-static A jttk0(J jt,B b,A a,A w){A z;I k,m=0,n,p,r,*s,*u;
+// take empty or from empty, a=take shape, w=takefrom. b=1 if a DOES NOT contain a 0
+// scaf* pass n and u as args
+static INLINE A jttk0(J jt,B b,A a,A w){A z;I k,m=0,n,p,r,*s,*u;
  r=AR(w); n=AN(a); u=AV(a); 
  if(!b){PRODX(m,n,u,1) ASSERT(m>IMIN,EVLIMIT); PRODX(m,r-n,n+AS(w),ABS(m))}
  GA(z,AT(w),m,r,AS(w)); 
+// scaf* if take not empty, set fill; if take empty, return
  s=AS(z); DO(n, p=u[i]; ASSERT(p>IMIN,EVLIMIT); *s++=ABS(p););
  if(m){k=bpnoun(AT(w)); mvc(k*m,AVn(r,z),k,jt->fillv);}
  R z;
 }
 
-static F2(jttks){F12IP;PROLOG(0092);A a1,q,x,y,z;B b,c;I an,m,r,*s,*u,*v;P*wp,*zp;
+// scaf* pass an and u as args
+static INLINE F2(jttks){F12IP;PROLOG(0092);A a1,q,x,y,z;B b,c;I an,m,r,*s,*u,*v;P*wp,*zp;
  an=AN(a); u=AV(a); r=AR(w); s=AS(w); 
  GASPARSE(z,AT(w),1,r,s); v=AS(z); DO(an, v[i]=ABS(u[i]););
  zp=PAV(z); wp=PAV(w);
@@ -48,12 +52,16 @@ static F2(jttks){F12IP;PROLOG(0092);A a1,q,x,y,z;B b,c;I an,m,r,*s,*u,*v;P*wp,*z
  EPILOG(z);
 }    /* take on sparse array w */
 
-// general take routine.  a is result frame followed by signed take values i. e. shape of result, w is array
-static F2(jttk){F12IP;PROLOG(0093);A y,z;B b=0;C*yv,*zv;I c,d,dy,dz,e,i,k,m,n,p,q,r,*s,t,*u;
+// general take/drop routine.  a is result frame followed by signed take values i. e. shape of result, w is array
+// scaf* pass n and u as args
+static F2(jttk){F12IP;PROLOG(0093);A y,z;C*yv,*zv;I c,d,dy,dz,e,i,k,m,n,p,q,r,*s,t,*u;
  n=AN(a); u=AV(a); r=AR(w); s=AS(w); t=AT(w);
  if(unlikely(ISSPARSE(t)))R tks(a,w);
- DO(n, if(!u[i]){b=1; break;}); if(!b)DO(r-n, if(!s[n+i]){b=1; break;});  // if empty take, or take from empty cell, set b
- if(((b-1)&AN(w))==0)R tk0(b,a,w);   // this handles empty w, so PROD OK below   b||!AN(w)
+// obsolete  DO(n, if(!u[i]){b=1; break;}); if(!b)DO(r-n, if(!s[n+i]){b=1; break;});  // if empty take, or take from empty cell, set b
+ DO(n, if(!u[i])goto emptytake;); DO(r-n, if(!s[n+i])goto emptytake;);  // if empty take, or take from empty cell, set b
+ if(unlikely(AN(w)==0)){B b=0; if(0){emptytake: b=1;} R tk0(b,a,w);}
+// obsolete  if(((b-1)&AN(w))==0)R tk0(b,a,w);   // this handles empty w, so PROD OK below   b||!AN(w)
+#if 1   // obsolete 
  k=bpnoun(t); z=w; c=q=1;  // c will be #cells for this axis
  // process take one axis at a time
  for(i=0;i<n;++i){I itemsize;
@@ -62,23 +70,25 @@ static F2(jttk){F12IP;PROLOG(0093);A y,z;B b=0;C*yv,*zv;I c,d,dy,dz,e,i,k,m,n,p,
    PROD(itemsize,r-i-1,s+i+1);  // size of item of cell
    DPMULDE(c*itemsize,q,d); GA(y,t,d,r,AS(z)); AS(y)[i]=q;  // this catches q=IMIN: mult error or GA error   d=#cells*itemsize*#taken items
    if(q>m)mvc(k*AN(y),CAVn(r,y),k,jt->fillv);   // overtake - fill the whole area
-   itemsize *= k; e=itemsize*MIN(m,q);  //  itemsize=in bytes; e=total bytes moved per item
+   itemsize*=k; e=itemsize*MIN(m,q);  //  itemsize=in bytes; e=total bytes moved per item
    dy=itemsize*q; yv=CAV(y);
    dz=itemsize*m; zv=CAV(z);
    m-=q; I yzdiff=dy-dz; yv+=REPSGN(p&m)&yzdiff; zv-=REPSGN(p&-m)&yzdiff;
    DQ(c, MC(yv,zv,e); yv+=dy; zv+=dz;);
    z=y;
   }
  }
+#else
+// scaf* set fill if needed
+#endif
  EPILOG(z);
 }
 
 F2(jttake){F12IP;A s;I acr,af,ar,n,*v,wcr,wf,wr;
- 
  ARGCHK2(a,w); I wt=AT(w);  // wt=type of w
  acr=jt->ranks>>RANKTX; wcr=(RANKT)jt->ranks; RESETRANK;  // save ranks before they are destroyed 
  if(unlikely(ISSPARSE(AT(a))))RZ(a=denseit(a));  //if a is empty this destroys jt->ranks
- if(likely(!ISSPARSE(wt)))RZ(w=jtsetfv1(jt,w,AT(w)));   // pity to do this before we know we need fill
+ if(likely(!ISSPARSE(wt)))RZ(w=jtsetfv1(jt,w,AT(w)));   // scaf* wrong to do this before we know we need fill
  ar=AR(a); acr=ar<acr?ar:acr; af=ar-acr;  // ?r=rank, ?cr=cell rank, ?f=length of frame
  wr=AR(w); wcr=wr<wcr?wr:wcr; wf=wr-wcr;
  if(((af-1)&(acr-2))>=0){
@@ -99,6 +109,7 @@ F2(jttake){F12IP;A s;I acr,af,ar,n,*v,wcr,wf,wr;
   if(i<AN(s)){
    s=ca(s); if(!ISDENSETYPE(AT(a),FL))RZ(a=ccvt(FL,a,0));  // copy area we are going to change; put a in a form where we can recognize infinity
    for(;i<AN(s);++i){if(DAV(a)[i]==IMIN)IAV(s)[i]=IMIN;else if(INF(DAV(a)[i]))IAV(s)[i]=wcr?ws[wf+i]:1;}  // kludge.  The problem is which huge values to consider infinite.  This is how it was done
+// scaf what's this IMIN business?  32-bit? avoid all the conversions
   }
  }
  a=s;
@@ -107,7 +118,7 @@ F2(jttake){F12IP;A s;I acr,af,ar,n,*v,wcr,wf,wr;
   I tklen = IAV(a)[0];  // get the one number in a, the take amount
   I tkasign = REPSGN(tklen);  // 0 if tklen nonneg, ~0 if neg
   I nitems = ws[0];  // number of items of w
-  I tkabs = (tklen^tkasign)-tkasign;  // ABS(tklen)
+  I tkabs = (tklen^tkasign)-tkasign;  // (UI)ABS(tklen)
   if((UI)tkabs<=(UI)nitems) {  // if this is not an overtake...  (unsigned to handle overflow, if tklen=IMIN).
    // calculate offset
    I woffset = tkasign&(tklen + nitems);   // x+#y if x neg, 0 if x pos
@@ -191,7 +202,7 @@ static F1(jtrsh0){F12IP;A x,y;I wcr,wf,wr,*ws;
 }
 
 F1(jthead){F12IP;I wcr,wf,wr;
- 
+
  ARGCHK1(w);
  wr=AR(w); wcr=(RANKT)jt->ranks; wcr=wr<wcr?wr:wcr; wf=wr-wcr;  // no RESETRANK so that we can pass rank into other code
  if(unlikely(!wcr)){RETF(RETARG(w))  // {."0, a NOP

jsrc/vt.c

@@ -664,6 +664,34 @@ xx =: i. 1e6  1 1
 NB. fill not used if not needed
 1: (i. 2 3 2 5 6) ,!.'a'"3 i. 2 5 6
 
+NB. x fill t v y, checking ,!.fill"v
+t =: {{ xx =: x [ yy =: y [ mm =: m [ nn =: n
+        if. 0=#m do. assert. (x {{ x , y }}"n y) -: x ,"n y else. assert. (x m {{ x ,!.m y }}"n y) -: x ,!.m"n y end. }}
+
+NB. try all ranks of x/y, and with a 0 at each position
+tt =: {{
+xxx =: x [ yyy =: y
+cr =: x ,&#&$ y  NB. x and y are always a cell of result
+x (m t _) y  NB. check the single cell
+for_r. ,/ ,"0/~ 0 1 2 do.
+  'prefx prefy' =: pref =: {.&(?5 5)&.> r
+  xs =. (prefx , $x) ($!.'',) y [ ys =. (prefy , $y) ($!.'',) x 
+  xs (m t cr) ys  NB. check rank without 0 insertion
+  for_s.  ,/ x ,"0/&i.&#&$ y do.  NB. insert 0 in each position
+    'sx sy' =: s
+    ((($!.'',)~  prefx , 0 sx} $) x) (m t cr) ((($!.'',)~   prefy , 0 sy} $) y)
+  end.
+end.
+1
+}}
+
+(>: ((>: 0.2 0.5 0.8 I. ? 0) ?@$ 8) ?@$ 100) ('' tt) (>: ((>: 0.2 0.5 0.8 I. ? 0) ?@$ 8) ?@$ 100) 
+(>: ((>: 0.2 0.5 0.8 I. ? 0) ?@$ 8) ?@$ 100) (100 tt) (>: ((>: 0.2 0.5 0.8 I. ? 0) ?@$ 8) ?@$ 100) 
+(>: ((>: 0.2 0.5 0.8 I. ? 0) ?@$ 8) ?@$ 100x) ('' tt) (>: ((>: 0.2 0.5 0.8 I. ? 0) ?@$ 8) ?@$ 100) 
+(>: ((>: 0.2 0.5 0.8 I. ? 0) ?@$ 8) ?@$ 100x) (100x tt) (>: ((>: 0.2 0.5 0.8 I. ? 0) ?@$ 8) ?@$ 100) 
+(<"0 >: ((>: 0.2 0.5 0.8 I. ? 0) ?@$ 8) ?@$ 100) ('' tt) (<"0 >: ((>: 0.2 0.5 0.8 I. ? 0) ?@$ 8) ?@$ 100) 
+(<"0 >: ((>: 0.2 0.5 0.8 I. ? 0) ?@$ 8) ?@$ 100) ((<5) tt) (<"0 >: ((>: 0.2 0.5 0.8 I. ? 0) ?@$ 8) ?@$ 100)
+
 
 randfini''
 

test/g320.ijs

@@ -581,6 +581,9 @@ a =: i. 1000 1000
 (7!:2 '; 2 (<@(]@]/)\) a') > 1.5 * 7!:2 '2 ;@:(<@(]@]/)\) a'
 (; 2 (<@(]@]/)\) a) -: 2 ;@:(<@(]@]/)\) a
 
+6 _. _. _. 15 (= +. *.&(128!:5))"0 (3) +/\ 1 2 3 _. 4 5 6
+2 _. _. _. 5 (= +. *.&(128!:5))"0 (3) (+/%#)\ 1 2 3 _. 4 5 6
+
 
 randfini''