forqfvec (iterative qfminim)

[Bill Allombert <Bill.Allombert@math.u-bordeaux1.fr>]

Hello PARI developers,

Please find an experimental patch that provide a function forqfvec
that allow to iterate on the vectors that qfminim would returns.
This was suggested by Ariel Pacetti.

This is the synopsis:

forqfvec(v,q,b,expr): q being a square and symmetric matrix representing a
positive definite quadratic form, evaluate expr for all vector v such that
q(v)<b.

Some example of use:
my(s);forqfvecsmall(v,matid(8),30,s+=Vec(v));s
##
my(s);forqfvecsmall(v,matid(8),,s+=Vec(v));s
##
my(s);forqfvecsmall(v,matid(8),30,if(s++<10,print(v)));
##
forqfvecsmall(v,matid(8),,return(Vec(v)))

Comments welcome!

Cheers,
Bill.

commit da0037ad04511d5bed3666842297eb3fd7cc792f
Author: Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>
Date:   Sun Feb 7 19:16:54 2010 +0100

    Add function forqfvec.

diff --git a/src/basemath/bibli1.c b/src/basemath/bibli1.c
index 02e68f0..0a33c92 100644
--- a/src/basemath/bibli1.c
+++ b/src/basemath/bibli1.c
@@ -1647,6 +1647,99 @@ addcolumntomatrix(GEN V, GEN invp, GEN L)
   return 1;
 }
 
+void
+forqfvec(void *E, void (*fun)(void *, GEN), GEN a, GEN BORNE)
+{
+  GEN x,u,r;
+  long n = lg(a), i, j, k;
+  pari_sp av1, av=avma, lim;
+  double p,BOUND,*v,*y,*z,**q;
+  const double eps = 0.0001;
+
+  if (typ(a) != t_MAT || !RgM_is_ZM(a)) pari_err(typeer,"forqfvec");
+  if (!BORNE) BORNE = gen_0;
+  else
+  {
+    BORNE = gfloor(BORNE);
+    if (typ(BORNE) != t_INT) pari_err(typeer, "minim0");
+  }
+
+  minim_alloc(n, &q, &x, &y, &z, &v);
+  av1 = avma;
+
+  u = lllgramint(a);
+  if (lg(u) != n) pari_err(talker,"not a definite form in minim0");
+  a = qf_apply_ZM(a,u);
+
+  n--;
+  a = RgM_gtofp(a, DEFAULTPREC);
+  r = qfgaussred_positive(a);
+  if (!r) pari_err(precer, "minim0");
+  for (j=1; j<=n; j++)
+  {
+    v[j] = rtodbl(gcoeff(r,j,j));
+    for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j));
+  }
+
+  if (gequal0(BORNE))
+  {
+    double c;
+    p = rtodbl(gcoeff(a,1,1));
+    for (i=2; i<=n; i++) { c = rtodbl(gcoeff(a,i,i)); if (c < p) p = c; }
+    BORNE = roundr(dbltor(p));
+  }
+  else
+    p = gtodouble(BORNE);
+  BOUND = p * (1 + eps);
+  if (BOUND == p) pari_err(precer, "minim0");
+
+  av1 = avma; lim = stack_lim(av1,1);
+  k = n; y[n] = z[n] = 0;
+  x[n] = (long)sqrt(BOUND/v[n]);
+  for(;;x[1]--)
+  {
+    do
+    {
+      if (k>1)
+      {
+        long l = k-1;
+        z[l] = 0;
+        for (j=k; j<=n; j++) z[l] += q[l][j]*x[j];
+        p = (double)x[k] + z[k];
+        y[l] = y[k] + p*p*v[k];
+        x[l] = (long)floor(sqrt((BOUND-y[l])/v[l])-z[l]);
+        k = l;
+      }
+      for(;;)
+      {
+        p = (double)x[k] + z[k];
+        if (y[k] + p*p*v[k] <= BOUND) break;
+        k++; x[k]--;
+      }
+    } while (k > 1);
+    if (! x[1] && y[1]<=eps) break;
+
+    p = (double)x[1] + z[1]; p = y[1] + p*p*v[1]; /* norm(x) */
+    fun(E, x); if (loop_break()) break;
+  }
+  avma = av;
+}
+
+static void
+gp_evalvoid(void *E, GEN x)
+{
+  set_lex(-1,x);
+  closure_evalvoid((GEN)E);
+}
+
+void
+forqfvec0(GEN a, GEN BORNE, GEN code)
+{
+  push_lex(gen_0, code);
+  forqfvec((void*) code, &gp_evalvoid, a, BORNE);
+  pop_lex(1);
+}
+
 /* Minimal vectors for the integral definite quadratic form: a.
  * Result u:
  *   u[1]= Number of vectors of square norm <= BORNE
diff --git a/src/functions/linear_algebra/forqfvec b/src/functions/linear_algebra/forqfvec
new file mode 100644
index 0000000..cfb559e
--- /dev/null
+++ b/src/functions/linear_algebra/forqfvec
@@ -0,0 +1,9 @@
+Function: forqfvec
+Section: linear_algebra
+C-Name: forqfvec0
+Prototype: vVGDGI
+Help:forqfvec(v,q,b,expr): q being a square and symmetric matrix
+ representing a positive definite quadratic form, evaluate expr for all 
+ vector v such that q(v)<b.
+Doc: $q$ being a square and symmetric matrix representing a positive definite
+ quadratic form, evaluate \kbd{expr} for all vector $v$ such that $q(v)<b$.
diff --git a/src/headers/paridecl.h b/src/headers/paridecl.h
index f55771b..6d6b1b1 100644
--- a/src/headers/paridecl.h
+++ b/src/headers/paridecl.h
@@ -1075,6 +1075,8 @@ GEN     Q_from_QR(GEN x, long prec);
 GEN     R_from_QR(GEN x, long prec);
 GEN     algdep(GEN x, long n);
 GEN     algdep0(GEN x, long n, long bit);
+void    forqfvec(void *E, void (*fun)(void *, GEN), GEN a, GEN BORNE);
+void    forqfvec0(GEN a, GEN BORNE, GEN code);
 GEN     gram_matrix(GEN M);
 GEN     lindep(GEN x);
 GEN     lindep0(GEN x, long flag);