Merge branch 'master' of jsoftware.com:jsource

# Conflicts:
#	jsrc/vfrom.c

jlibrary/addons/dev/modular/manifest.ijs

@@ -6,7 +6,7 @@ DESCRIPTION=: 0 : 0
 modular arithmetic
 )
 
-VERSION=: '1.0.3'
+VERSION=: '1.0.4'
 
 FILES=: 0 : 0
 modular.ijs

jlibrary/addons/dev/modular/modular.ijs

@@ -74,3 +74,72 @@ NB. errors:
 %. m. M i. 3 3
 
 )
+
+NB. modular determinant
+NB. modification of gaussian elimination from
+NB. math/misc/linear.ijs
+
+NB. failsafe case; never seen in practice
+mdetx=: 4 : 'x|-/ . * x: y'
+
+NB. main adverb
+modular_det=: 1 : 0
+mult=.* m. m
+diff=. - m. m
+div=.% m. m
+A=. y
+r =. #A
+i=. 0
+p=. 1
+while. i<r-1 do.
+  k=. 1 i.~ 1=m +. tc=.i}. i{"1 A   NB. need an invertible pivot
+  if. r=i+k do.  NB. if so, no invertible pivot in lower column
+     if. 0=#tc-.0 do. 0 return. NB. level 1: lower column all 0, det is 0
+     end.
+     p=. p mult mult/(<0 1)|:(i,i){.A    NB. reduce to the unprocessd block
+     B=.(i,i)}.A
+     if. 2=#B do.               NB. level 2: do 2 by 2 det
+       'a b c d'=.,B
+       p mult (a mult d)diff b mult c return.
+     end.
+     if. 1 e. m +. {.B do.      NB. level 3: look for pivot in  row
+       p mult m modular_det_j_ |:B return.
+     end.
+     if. 1 e. t=.,m +. B do.   NB. level 4: look for a pivot anywhere (won't be in 0th row or col by above levels)
+       'I J'=.($B)#: t i. 1
+       B=.(<0,I)C. B
+       B=.(<0,J) C."1 B
+       p mult m modular_det_j_ B return.
+     end.
+     if. 1 < >./g=.+./B do.     NB. level 5: try to factor gcds out of columns
+      p mult (mult/g) mult m modular_det_j_ B %"1 g return.
+     end.
+     if. 1 < >./g=.+./"1 B do.  NB. level 6: try to factor gcds out of rows
+      p mult (mult/g) mult m modular_det_j_ B %"1 0 g return.
+     end.
+     if. 1 e. ,m +. S=.({. , }. diff"1 {.)B do.  NB. level 7: try row shear
+     p mult m modular_det_j_ S  return.
+     end.
+     if. 1 e. ,m +. S=.({."1 ,. }."1 diff"1 0 {."1)B do.  NB. level 8: try column shear
+     p mult m modular_det_j_ S  return.
+     end.
+     p mult m mdetx_j_ B     return.       NB. level 9: use extended
+  end.
+NB. end of dealing with no easy pivot
+NB. this is where the main work is done and most of the above is usually skipped:
+    if. k do.     NB. move a pivot from lower in column
+      A=. (<i,i+k) C. A
+      p=. p mult _1
+    end.
+NB. Time to partial pivot
+    r1=.>:i
+    col=.i{"1 A
+    A=.(r1{.A),(r1}.A) diff (r1}.col) mult/ (i{A) div i{col
+    i=. >:i
+end.
+p mult mult/(<0 1)|:A
+)
+
+Mdet=:1 : 0
+(m modular_det_j_)"2 y
+)

jlibrary/system/main/stdlib.ijs

@@ -1,7 +1,7 @@
 18!:4 <'z'
 3 : 0 ''
 
-JLIB=: '9.5.1'
+JLIB=: '9.5.2'
 
 notdef=. 0: ~: 4!:0 @ <
 hostpathsep=: ('/\'{~6=9!:12'')&(I. @ (e.&'/\')@] })

jsrc/jgmp.h

@@ -296,6 +296,7 @@ EXTERN Q Q__;            // x: __ NB. _1r0 internal form
 #define jmpz_bin_ui __gmpz_bin_ui            // https://gmplib.org/manual/Number-Theoretic-Functions
 #define jmpz_cdiv_q __gmpz_cdiv_q            // https://gmplib.org/manual/Integer-Division
 #define jmpz_fdiv_q __gmpz_fdiv_q            // https://gmplib.org/manual/Integer-Division
+#define jmpz_fits_slong_p __gmpz_fits_slong_p // https://gmplib.org/manual/Miscellaneous-Integer-Functions
 #define jmpz_clear __gmpz_clear              // https://gmplib.org/manual/Initializing-Integers
 #define jmpz_cmp __gmpz_cmp                  // https://gmplib.org/manual/Integer-Comparisons
 #define jmpz_cmpabs_ui __gmpz_cmpabs_ui      // https://gmplib.org/manual/Integer-Comparisons
@@ -305,6 +306,7 @@ EXTERN Q Q__;            // x: __ NB. _1r0 internal form
 #define jmpz_fdiv_qr __gmpz_fdiv_qr          // https://gmplib.org/manual/Integer-Division
 #define jmpz_fdiv_qr_ui __gmpz_fdiv_qr_ui    // https://gmplib.org/manual/Integer-Division
 #define jmpz_fdiv_r __gmpz_fdiv_r            // https://gmplib.org/manual/Integer-Division
+#define jmpz_fits_slong_p __gmpz_fits_slong_p // https://gmplib.org/manual/Miscellaneous-Integer-Functions
 #define jmpz_gcd __gmpz_gcd                  // https://gmplib.org/manual/Number-Theoretic-Functions
 #define jmpz_get_d __gmpz_get_d              // https://gmplib.org/manual/Converting-Integers
 #define jmpz_get_d_2exp __gmpz_get_d_2exp    // https://gmplib.org/manual/Converting-Integers
@@ -356,6 +358,7 @@ EXTERN void (*jmpz_fdiv_q)(mpz_t, const mpz_t, const mpz_t);
 EXTERN void (*jmpz_fdiv_qr)(mpz_t, mpz_t, const mpz_t, const mpz_t);
 EXTERN void (*jmpz_fdiv_qr_ui)(mpz_t, mpz_t, const mpz_t, mpir_ui);
 EXTERN void (*jmpz_fdiv_r)(mpz_t, const mpz_t, const mpz_t);
+EXTERN int  (*jmpz_fits_slong_p)(const mpz_t);
 EXTERN void (*jmpz_gcd)(mpz_t, const mpz_t, const mpz_t);
 EXTERN D    (*jmpz_get_d)(const mpz_t);
 EXTERN D    (*jmpz_get_d_2exp)(long int*, const mpz_t);
@@ -501,6 +504,7 @@ extern void jfree4gmp(void*, size_t);
 #define icmpXX(x,y) shimXX(jmpz_cmp, x,y)           // *x-y
 #define icmpXI(x,y) shimXI(jmpz_cmp_si, x,y)
 #define Xfdiv_qXX(x, y) XshimXX(jmpz_fdiv_q, x, y)      // <.x%y
+#define ifits_slong_pX(x) shimX(jmpz_fits_slong_p, x)   // x = IMIN<.IMAX>.x (nonzero result if true)
 #define XgcdXX(x, y) XshimXX(jmpz_gcd, x, y)            // x+.y
 #define DgetX(y) shimX(jmpz_get_d, y)    // y+0.0
 #define IgetX(y) shimX(jmpz_get_si, y)   // (I)y  NB. UINT_MAX&|&.(-&INT_MIN) y

jsrc/jgmpinit.c

@@ -419,6 +419,7 @@ void jgmpinit(C*libpath) {
  jgmpfn(mpz_bin_ui);       // https://gmplib.org/manual/Number-Theoretic-Functions
  jgmpfn(mpz_cdiv_q);       // https://gmplib.org/manual/Integer-Division
  jgmpfn(mpz_fdiv_q);       // https://gmplib.org/manual/Integer-Division
+ jgmpfn(mpz_fits_slong_p); // https://gmplib.org/manual/Miscellaneous-Integer-Functions
  jgmpfn(mpz_clear);        // https://gmplib.org/manual/Initializing-Integers
  jgmpfn(mpz_cmp);          // https://gmplib.org/manual/Integer-Comparisons
  jgmpfn(mpz_cmpabs_ui);    // https://gmplib.org/manual/Integer-Comparisons

jsrc/u.c

@@ -235,8 +235,17 @@ I jti0(J jt,A w){ARGCHK1(w);
   }
   R cval;  // too-large values don't convert, handle separately
  }
+ if (AT(w)&XNUM) {
+  X x= XAV(w)[0];
+  if (ifits_slong_pX(x)) {
+   I r= IgetX(x);
+   if (unlikely(IMIN==r)) R -IMAX;
+   R r;
+  }
+  if (XSGN(x)>0) R IMAX; else R -IMAX;
+ }
  if(!(w=vi(w)))R 0; ASSERT(!AR(w),EVRANK);
- if (AT(w)&(RAT+XNUM) || ISGMP(w)) SEGFAULT;
+ if (ISGMP(w)||AT(w)&RAT) SEGFAULT; // this should never happen
  R IAV(w)[0];
 }  // can't move the ASSERT earlier without breaking a lot of tests
 

jsrc/vfrom.c

@@ -869,6 +869,8 @@ A jtcompidx(J jt,I axislen,A ind){
  // to hold a bitmask of the values that have not been crossed off
  RZ(ind=likely(ISDENSETYPE(AT(ind),INT))?ind:cvt(INT,ind));  // ind is now an INT vector, possibly the input argument
  I allolen=MAX(axislen,1);  // number of words needed.  There must be at least one value crossed off, but we always need at least 1 word for bitmask
+                            // We also reserve one word of buffer between the indices and the bitmask since, in the case of (for example) (<<<_1) { i.65, we could end up overwriting the last word of the bitmask before reading it
+                            // An alternative would be to pipeline the loop, loading the mask for the next iteration before completing the previous iteration, but this would be annoying and bloat
  A z; GATV0(z,INT,allolen,1) I *zv0=IAV1(z), *zv=zv0;   // allocate the result/temp block.  
  I bwds=(axislen+(BW-1))>>LGBW;  // number of words needed: one bit for each valid index vallue
  I *bv=zv+allolen-bwds; mvc(bwds*SZI,bv,SY_64?4*SZI:2*SZI,validitymask); bv[bwds-1]=~((~1ll)<<((axislen-1)&(BW-1)));  // fill the block with 1s to indicate we need to write; clear ending 0s

jsrc/vrand.c

@@ -620,6 +620,9 @@ static F2(jtrollksub){A z;I an,*av,k,m1,n,p,q,r,sh;UI m,mk,s,t,*u,x=jt->rngdata-
   // Output the rest, one bit at a time
   t=NEXT;  // Get random # for bits
   B*c=(B*)u; DQ(r&(SZI-1), *c++=1&t; t>>=1;);
+ }else if(unlikely(IMAX==m&&XNUM&AT(w))){
+  // extended output
+  X*xv=XAV(z); RZ(xv); AT(z)= XNUM; mvc(n*SZI, xv, SZI, MEMSET00); DO(n,(*xv++)=XAV(roll(w))[0];)
  }else{
   // integer output
   r=n; s=GMOF(m,x); if(s==x)s=0;
@@ -799,6 +802,7 @@ DF2(jtdeal){A z;I at,j,k,m,n,wt,*zv;UI c,s,t,x=jt->rngdata->rngparms[jt->rngdata
   // calculate the number of values to deal: m, plus a factor times the expected number of collisions, plus 2 for good measure.  Will never exceed n.  Repeats a little less than 1% of the time for n between 30 and 300
   A h=sc(m+4+(I)((n<1000?2.4:2.2)*((D)m+(D)n*(pow((((D)(n-1))/(D)n),(D)m)-1)))); NOUNROLL do{RZ(z=nub(rollksub(h,w)));}while(AN(z)<m); RZ(z=jttake(JTIPW,a,z));
 #else
+  if(unlikely(n==IMAX&&XNUM&AT(w))){A h= plus(num(2),a); NOUNROLL do{RZ(z=nub(rollksub(h,w)));}while(AN(z)<m);R take(a,z);}
   A h,y; I d,*hv,i,i1,p,q,*v,*yv;
   FULLHASHSIZE(2*m,INTSIZE,1,0,p);
   GATV0(h,INT,p,1); hv=AV(h); DO(p, hv[i]=0;);

test/g520.ijs

@@ -1510,6 +1510,11 @@ NB. for '' { i. 0 3 bug fix (3 2 0) -: $ (2 0$4) {"2 i. 3 0 3
 
 NB. x{y complementary indexing ------------------------------------------
 
+(i.64) -: (<<<_1) { i.65
+(i.64) -: (<<<64) { i.65
+(i.128) -: (<<<_1) { i.129
+(i.128) -: (<<<128) { i.129
+
 jot=:<$0
 
 NB. Boolean

test/g600.ijs

@@ -1015,6 +1015,9 @@ f =: ^
 0 (ger"0 -: appcyc"0) i. 0 [ cycinit ger=:$`*`f
 0 (ger"0 -: appcyc"0) i. 0 [ cycinit ger=:(5 5 5"_)`*`f
 
+5 = # 5 ? 2980293480239480239480239480239480239482039x
+63 < >./ 2 ^. 5 ? 2980293480239480239480239480239480239482039x
+
 a=:_1 4 6 8 3 5 8 _1 7 4
 'value error' -: undefname`0:`*: "0 etx a