added some utilities

libsrc/arith.cc

@@ -258,6 +258,18 @@ vector<long> sqfreedivs(long a, const vector<long>& plist)
  return dlist;
 }
 
+// gcc division truncates towards 0, while we need rounding, with a
+// consistent behaviour for halves (they go up here).
+//
+// For b>0, rounded_division(a,b) = q such that a/b = q + r/b with -1/2 <= r/b < 1/2
+long rounded_division(long a, long b)
+{
+  std::ldiv_t qr = ldiv(a, b);
+  long r = qr.rem, q = qr.quot;
+  long r2 = r<<1;
+  return (r2<-b? q-1: (r2>=b? q+1: q));
+}
+
 long mod(long a, long b)
 {long c;
  if (b<0) b=-b;
@@ -395,11 +407,34 @@ long chi2(long a)
 { static const long table8[8] = {0,1,0,-1,0,-1,0,1};
   return table8[posmod(a,8)];
 }
-
+
+// set root to rounded sqrt(a) if a>=0, return 1 iff exact
+int isqrt(long a, long& root)
+{
+  if (a<0) {return 0;}
+  root = round(sqrt(a));
+  return a==root*root;
+}
+
+// return rounded sqrt(a) (undefined for a<0)
+long isqrt(const long a)
+{
+  long r;
+  isqrt(a,r);
+  return r;
+}
+
+long squarefree_part(long d)
+{
+  if (d==0) return d;
+  long maxd = sqdivs(d).back();
+  return (d/maxd)/maxd;
+}
+
 long chi4(long a)
 { static const long table4[4] = {0,1,0,-1};
   return table4[posmod(a,4)];
-} 
+}
 
 long hilbert2(long a, long b)
 { static long table44[4][4] = {{0,0,0,0},
@@ -516,4 +551,31 @@ int is_valid_conductor(long n)
   return std::all_of(plist.begin(), plist.end(), [m] (const long& p) {return val(p,m)<=2;});
 }
 
+
+// a=b*q+r, return 1 iff r==0
+int divrem(long a, long b, long& q, long& r)
+{
+  std::ldiv_t qr = ldiv(a, b);
+  r = qr.rem;
+  q = qr.quot;
+  return (r==0);
+}
+
+// a=b*q+r, return 1 iff r==0
+int divrem(int a, int b, int& q, int& r)
+{
+  std::div_t qr = div(a, b);
+  r = qr.rem;
+  q = qr.quot;
+  return (r==0);
+}
+
+// return list of integers from first to last inclusive
+vector<long> range(long first, long last)
+{
+  vector<long> ans(last-first+1);
+  std::iota(ans.begin(), ans.end(), first);
+  return ans;
+}
+
 /* END OF FILE */

libsrc/eclib/arith.h

@@ -130,10 +130,19 @@ inline vector<long> sqfreedivs(long n)
   return sqfreedivs(n, pdivs(n));
 }
 
-
 long mod(long a, long modb); /* modulus in range plus or minus half mod */
 long posmod(long a, long modb); /* ordinary modulus, but definitely positive */
 
+// gcc division truncates towards 0, while we need rounding, with a
+// consistent behaviour for halves (they go up here).
+//
+// For b>0, rounded_division(a,b) = q such that a/b = q + r/b with -1/2 <= r/b < 1/2
+
+// The bigint version of a similar function roundover(a,b) is in
+// marith.h implmented as {bigint a0=(a%b); bigint c = (a-a0)/b;
+// if(2*a0>b) c+=1; return c;} which is not quite the same
+long rounded_division(long a, long b);
+
 long gcd(long, long);
 int gcd(int, int);
 long lcm(long, long);
@@ -149,19 +158,26 @@ inline int is_zero(long n) {return n==0;}
 long val(long factor, long number); // order of factor in number
 
 inline int divides(long factor,long number) 
-{ 
+{
   return (factor==0) ? (number==0) : (number%factor==0);
 }
 
-inline int ndivides(long factor,long number) 
+inline int ndivides(long factor,long number)
 {
   return (factor==0) ? (number!=0) : (number%factor!=0);
-  //  return !::divides(factor,number);
 }
 
+// a=b*q+r, return 1 iff r==0
+int divrem(long a, long b, long& q, long& r);
+int divrem(int a, int b, int& q, int& r);
+
 inline long m1pow(long a) {return (a%2 ?  -1 : +1);}
-inline int sign(long a)   {return (a==0? 0: (a>0? 1: -1));} 
-inline int sign(double a) {return (a==0? 0: (a>0? 1: -1));} 
+inline int sign(long a)   {return (a==0? 0: (a>0? 1: -1));}
+
+// set root to rounded sqrt(a) if a>=0, return 1 iff exact
+int isqrt(long a, long& root);
+// return rounded sqrt(a) (undefined for a<0)
+long isqrt(const long a);
 
 long chi2(long a);
 long chi4(long a);
@@ -173,5 +189,9 @@ int intlog2(long& n, long& e, int roundup);
 int is_squarefree(long n);
 int is_valid_conductor(long n);
 
+long squarefree_part(long d);
+
+// return list of integers from first to last inclusive
+vector<long> range(long first, long last);
 
 #endif

libsrc/eclib/mat.h

@@ -158,6 +158,8 @@ mat echelon(const mat& m, vec& pcols, vec& npcols,
 mat addscalar(const mat&, scalar);
 vec apply(const mat&, const vec&);
 
+mat reduce_modp(const mat& m, const scalar& p);
+
 // Construct an NTL mat_lzz_p (matrix mod p) from a mat mod pr
 mat_zz_p mat_zz_p_from_mat(const mat& M, scalar pr);
 

libsrc/eclib/vec.h

@@ -127,6 +127,8 @@ inline vec operator+(const vec& v1, const vec& v2)
 inline vec addmodp(const vec& v1, const vec& v2, scalar pr)
 { vec ans(v1); ans.addmodp(v2,pr); return ans;}
 
+vec reduce_modp(const vec& v, const scalar& p);
+
 inline vec operator-(const vec& v1, const vec& v2)
 { vec ans(v1); ans-=v2; return ans;}
 

libsrc/mat.cc

@@ -1131,6 +1131,16 @@ mat echelonl(const mat& entries, vec& pc, vec& npc,
   return ans;
 }
 
+mat reduce_modp(const mat& m, const scalar& p)
+{
+  if (p==0) return m;
+  long nr=m.nrows(), nc=m.ncols();
+  mat ans(nr,nc);
+  for(long i=1; i<=nr; i++)
+    for(long j=1; j<=nc; j++)
+      ans(i,j) = mod(m(i,j),p);
+  return ans;
+}
 
 void elimp(const mat& m, long r1, long r2, long pos, scalar pr)
 {

libsrc/vec.cc

@@ -114,6 +114,16 @@ void vec::addmodp(const vec& w, scalar pr)
     cerr << "Incompatible vecs in vec::addmodp"<<endl;
 }
 
+vec reduce_modp(const vec& v, const scalar& p)
+{
+  if (p==0) return v;
+  long d=dim(v);
+  vec ans(d);
+  for(long i=1; i<=d; i++)
+    ans[i] = mod(v[i], p);
+  return ans;
+}
+
 vec& vec::operator-=(const vec& q2)
 {
   scalar* vi=entries; scalar* wi=q2.entries; long i=d;