Make the type of empty arithmetic result the natural type given the arguments (instead of the type for boolean args)

Author: HenryHRich

jsrc/ca.c

@@ -59,7 +59,8 @@ static UI imodpow(UI x,I n,UI m){
  if(unlikely(m==1))R 0; UI z=1; // if n=0 result is 1 unless m=1, then 0
  UI mrecip=((UI)(-IMIN)/m); mrecip=(mrecip<<1)+((((UI)(-IMIN)-mrecip*m)<<1)>=m);  // 2^64%m, possibly low by as much as 2^-64
 #define modm(x) ({UI t; doubletype tt; tt=(doubletype)(x)*(doubletype)mrecip; t=tt>>BW; t=(x)-t*m; if(unlikely(t>=m))t-=m; t;})
- // x%m using mrecip.  mrecip may be low by 2^-64; x*mrecip is truncated.  Resulting modulus may be high by at most m
+ // x%m using mrecip.  x*mrecip is truncated, which gives the remainder.  mrecip may be 2^-64 low, so the remainder may be x*m/2^64 high.  This is never more than m too high, so a single correction suffices
+ // we expect the correction to be rare so we use a branch
 // #define modm(x) (x)%m
  while(n){UI zz=z*x; zz=modm(zz); x=x*x; x=modm(x); z=n&1?zz:z; n>>=1;}  //  repeated square/mod
  R z;

jsrc/va2.c

@@ -476,22 +476,29 @@ static A jtva2(J jt,AD * RESTRICT a,AD * RESTRICT w,AD * RESTRICT self,UI allran
    aadocv=&vainfo->p2[(at>>(INTX-1))+((at+wt)>>INTX)];
   }else{
 
-   // If an operand is empty, turn it to Boolean, and if the OTHER operand is non-numeric, turn that to Boolean too (leaving
-   //  rank and shape untouched).  This change to the other operand is notional only - we won't actually convert
-   // when there is an empty - but it guarantees that execution on an empty never fails.
    // If we switch a sparse nonnumeric matrix to boolean, that may be a space problem; but we don't
    // support nonnumeric sparse now
    // if an operand is sparse, replace its type with the corresponding non-sparse type, for purposes of testing operand precisions
    if(unlikely(ISSPARSE(at|wt))){
     at&=~SPARSE; wt&=~SPARSE;
     jtinplace=0;  // We use jtinplace==0 as a flag meaning 'sparse'
    }
+   // Get the result type and routine
+   adocv=var(self,at,wt);
+   if(unlikely(adocv.f==0)){
+    // There is no routine for these argument types.
+    // If an operand is empty, turn it to Boolean, and if the OTHER operand is non-numeric, turn that to Boolean too (leaving
+    //  rank and shape untouched).  This change to the other operand is notional only - we won't actually convert
+    // when there is an empty - but it guarantees that execution on an empty never fails.
+    at=((-AN(a)&(-AN(w)|-(at&NUMERIC)))>=0)?B01:at;
+    wt=((-AN(w)&(-AN(a)|-(wt&NUMERIC)))>=0)?B01:wt;
+// obsolete     if(AN(a)==0){at=B01;if(!(wt&NUMERIC))wt=B01;}  // switch empty arg to Boolean & ensure compatibility with other arg
+// obsolete     if(AN(w)==0){wt=B01;if(!(at&NUMERIC))at=B01;}
+    adocv=var(self,at,wt);   // rerun the decode with safer types
+   }
+   aadocv=&adocv;  // we save the address of the struct
    // We could allocate the result block here & avoid the test after the allocation later.  But we would have to check for agreement etc
-   if(AN(a)==0){at=B01;if(!(wt&NUMERIC))wt=B01;}  // switch empty arg to Boolean & ensure compatibility with other arg
-   if(AN(w)==0){wt=B01;if(!(at&NUMERIC))at=B01;}
-
-   // Figure out the result type.  Don't signal the error from it yet, because domain has lower priority than agreement
-   adocv=var(self,at,wt); aadocv=&adocv;
+   // Don't signal domain error on the types yet, because domain has lower priority than agreement
   }
  }
 

test/g0a.ijs

@@ -56,7 +56,41 @@ NB. handling empty arrays -----------------------------------------------
 (":0 2 3$s:@<"0 'a') -: ":0 2 3$s:@<"0&.> <'a'
 (":0 2 3$s:@<"0 'a') -: ":0 2 3$<"0@s: <'a'
 
-4!:55 ;:'n'
+NB. type of empty arithmetic arrays
+empties =. (0$0);'';(0$00);(0$0.);(0$0.j0.);(0$a:);(0$1x);(0$1r2)
+nonempties =. 0;' ';00;0.;0.j0.;a:;1x;1r2
+(".;._2 (0 : 0)) -: empties 3!:0@:+&>/ nonempties
+  4 4   4  8 16 4  64 128
+  4 4   4  8 16 4  64 128
+  4 4   4  8 16 4  64 128
+  8 4   8  8 16 4   8   8
+ 16 4  16 16 16 4  16  16
+  4 4   4  8 16 4  64 128
+ 64 4  64  8 16 4  64 128
+128 4 128  8 16 4 128 128
+)
+(".;._2 (0 : 0)) -: empties 3!:0@:*&>/ nonempties
+  1 1   4  8 16 1  64 128
+  1 1   4  8 16 1  64 128
+  4 1   4  8 16 1  64 128
+  8 1   8  8 16 1   8   8
+ 16 1  16 16 16 1  16  16
+  1 1   4  8 16 1  64 128
+ 64 1  64  8 16 1  64 128
+128 1 128  8 16 1 128 128
+)
+(".;._2 (0 : 0)) -: empties 3!:0@:-&>/ nonempties
+  4 4   4  8 16 4  64 128
+  4 4   4  8 16 4  64 128
+  4 4   4  8 16 4  64 128
+  8 4   8  8 16 4   8   8
+ 16 4  16 16 16 4  16  16
+  4 4   4  8 16 4  64 128
+ 64 4  64  8 16 4  64 128
+128 4 128  8 16 4 128 128
+)
+
+4!:55 ;:'n empties nonempties'
 
 
 

test/g520.ijs

@@ -137,6 +137,10 @@ for_t. i. 4 do.
   assert. (_6 ]\ 0 0 1 0 0 0  0 0 2e_20 0 0 0) ('' run128_9) 2;(,."1 (_2) ]\ 00 0);(0$00);(0$0.0);M;0.0;Yt
   assert. (_6 ]\ 1 2 0 0 0 0  0 2e_20 0 0 0 0) ('' run128_9) 6;(,."1 (_2) ]\ 0 2);(00 1);(1. 2.);M;0.0;Yt
   assert. (_6 ]\ 0 0 1 2 0 3  0 0 2e_20 6e_20 0 15e_20) ('' run128_9) 8;(,."1 (_2) ]\ (4$0) , 3 3);((3$0) , 2 3 5);((3$0), 1. 2. 3.);M;0.0;Yt
+  NB. change sign of result when Yt = 8 & fetch is from a Slack
+  assert. (- _6 ]\ 1 0 0 0 0 0  0 0 0 0 0 0) ('' run128_9) 00;(,."1 (_2) ]\ 00 0);(0$00);(0$0.0);M;0.0;(8 (0}) Yt)
+  assert. (- _6 ]\ 0 0 1 0 0 0  0 0 2e_20 0 0 0) ('' run128_9) 2;(,."1 (_2) ]\ 00 0);(0$00);(0$0.0);M;0.0;(8 (2}) Yt)
+  assert. (_6 ]\ 1 2 0 0 0 0  0 2e_20 0 0 0 0) ('' run128_9) 6;(,."1 (_2) ]\ 0 2);(00 1);(1. 2.);M;0.0;(8 (6}) Yt)
   NB. clamping to threshold
   M =. 1 1e_20 (<1;3;2 3)} 2 {. ,: |: _4 ]\ _1. 0 1 2  _1e_5 _2e_9 0 1e_9   0 0 0 0 0 0 0 1e_10   NB. input by columns
   assert. (2 {. ,: _1. 0 1 2) ('' run128_9) 00;(,."1 (_2) ]\ 00 0);(0$00);(0$0.0);M;0.5;Yt