Re: qfbprimeform(-3,1)--> stack corruption

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

On Thu, 2 Oct 2003, Bill Allombert wrote:
> ? qfbprimeform(-3,1)
> %1 = Qfb(1, 1, 33554435)
>
> This looks like a stack corruption do to an abuse of setsigne(y[3],1).

Actually data corruption. gzero's sign was set to 1 [ and gzero is not part
of the stack ]

> Unfortunately, there is no warranty that a PARI function return a value
> on the stack rather than gzero or gun. If we can restrain ourself to use
> such hack, we should write a macro that do it correctly.

I agree with the general comment. Although I don't see how to write a general
macro which would not compute a possibly unnecessary signe(), as in your patch
in fact.

Here is a(n ever so slightly) more efficient patch, which saves a few bitmask
operations [ committed to CVS ].

This original formula is

  av = avma; y[3] = lpileupto(av, shifti(subsi(1, D), -2));

which is wasteful since computing 1 - D, and copying the result (gerepileupto)
is unnecessary when D is even, making the routine more than twice slower than
it should be in half the cases. I made a mistake when "optimizing" this.

     Karim.

Index: src/basemath/arith2.c
===================================================================
RCS file: /home/megrez/cvsroot/pari/src/basemath/arith2.c,v
retrieving revision 1.73
diff -u -w -r1.73 arith2.c
--- src/basemath/arith2.c       2003/09/08 08:20:58     1.73
+++ src/basemath/arith2.c       2003/10/03 09:44:55
@@ -1669,22 +1669,25 @@
   GEN y = cgetg(4,t_QFI);
   long isodd;

-  if (typ(D) != t_INT || signe(D) >= 0) err(typeer,"real_unit_form_by_disc");
-  switch(4 - mod4(D))
+  if (typ(D) != t_INT || signe(D) >= 0) err(typeer,"imag_unit_form_by_disc");
+  switch(mod4(D))
   {
-    case 2:
-    case 3: err(funder2,"imag_unit_form_by_disc");
+    case 0: isodd = 0; break;
+    case 3: isodd = 1; break;
+    default: err(funder2,"imag_unit_form_by_disc");
+             return NULL; /* not reached */
   }
-  y[1] = un; isodd = mpodd(D);
+  y[1] = un;
   y[2] = isodd? un: zero;
-  /* y[3] = (1-D) / 4 or -D / 4, whichever is an integer */
-  y[3] = lshifti(D,-2); setsigne(y[3],1);
+  /* upon return, y[3] = (1-D) / 4 or -D / 4, whichever is an integer */
+  y[3] = lshifti(D,-2);
   if (isodd)
   {
     pari_sp av = avma;
-    y[3] = lpileuptoint(av, addis((GEN)y[3],1));
+    y[3] = lpileuptoint(av, addis((GEN)y[3],-1));
   }
-  return y;
+  /* at this point y[3] < 0 */
+  setsigne(y[3], 1); return y;
 }

 GEN

-- 
Karim Belabas                     Tel: (+33) (0)1 69 15 57 48
Dép. de Mathématiques, Bât. 425   Fax: (+33) (0)1 69 15 60 19
Université Paris-Sud              http://www.math.u-psud.fr/~belabas/
F-91405 Orsay (France)            http://www.parigp-home.de/  [PARI/GP]