Re: forqfvec (iterative qfminim)

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

Hello PARI-dev,

Ariel told me forqfvec did not work with non LLL reduced matrices,
so here a new patch that fix that, and a new example:

? my(z=vector(30));forqfvec(v,matid(8),30,z[v~*v]++);z
%1 = [8, 56, 224, 568, 1008, 1568, 2752, 4664, 6056, 7056, 10656,
15904, 17584, 19264, 28224, 37432, 39312, 42392, 54880, 71568, 77056,
74592, 97344, 130592, 126008, 123088, 163520, 195392, 195120, 197568]

Cheers,
Bill.

diff --git a/src/basemath/bibli1.c b/src/basemath/bibli1.c
index 69882e0..3b48975 100644
--- a/src/basemath/bibli1.c
+++ b/src/basemath/bibli1.c
@@ -1648,6 +1648,107 @@ addcolumntomatrix(GEN V, GEN invp, GEN L)
   return 1;
 }
 
+struct qfvec
+{
+  GEN a, r, u;
+};
+
+static void
+forqfvec_init(struct qfvec *qv, GEN a)
+{
+  if (typ(a) != t_MAT || !RgM_is_ZM(a)) pari_err(typeer, "forqfvec");
+  qv->u = lllgramint(a);
+  if (lg(qv->u) != lg(a)) pari_err(talker, "not a definite form in minim0");
+  qv->a = RgM_gtofp(qf_apply_ZM(a, qv->u), DEFAULTPREC);
+  qv->r = qfgaussred_positive(qv->a);
+  if (!qv->r) pari_err(precer, "minim0");
+}
+ 
+static void
+forqfvec(void *E, long (*fun)(void *, GEN, GEN, double), struct qfvec *qv, GEN BORNE)
+{
+  GEN x, a = qv->a, r = qv->r, u = qv->u;
+  long n = lg(a), i, j, k;
+  double p,BOUND,*v,*y,*z,**q;
+  const double eps = 0.0001;
+  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);
+  n--;
+  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");
+
+  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) */
+    if (fun(E, u, x, p)) break;
+  }
+}
+
+static long
+_gp_forqf(void *E, GEN u, GEN x, double p)
+{
+  pari_sp av = avma;
+  set_lex(-1, ZM_zc_mul(u, x));
+  closure_evalvoid((GEN)E);
+  avma = av;
+  return loop_break();
+}
+
+void
+forqfvec0(GEN a, GEN BORNE, GEN code)
+{
+  pari_sp av = avma;
+  struct qfvec qv;
+  forqfvec_init(&qv, a);
+  push_lex(gen_0, code);
+  forqfvec((void*) code, &_gp_forqf, &qv, BORNE);
+  pop_lex(1);
+  avma = av;
+}
+
 /* 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..c3ee950
--- /dev/null
+++ b/src/functions/linear_algebra/forqfvec
@@ -0,0 +1,15 @@
+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$.
+Variant: The following function is also available:
+ \fun{void}{forqfvec}{void *E, long (*fun)(void *, GEN, double), GEN q, GEN b}:
+ Evaluate \kbd{fun(E,v,m)} on all $v$ such that $q(v)<b$, where $v$ is a
+ \typ{VECSMALL} and $m=q(v)$ is a C double. The function \kbd{fun} must
+ return $0$, unless \kbd{forqfvec} should stop, in which case, it should
+ return $1$.
diff --git a/src/headers/paridecl.h b/src/headers/paridecl.h
index 2708168..9f83d4a 100644
--- a/src/headers/paridecl.h
+++ b/src/headers/paridecl.h
@@ -1084,6 +1084,7 @@ 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    forqfvec0(GEN a, GEN BORNE, GEN code);
 GEN     gram_matrix(GEN M);
 GEN     lindep(GEN x);
 GEN     lindep0(GEN x, long flag);