Re: bug report

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

[Kiyoshi Ohgishi:]
> 
> I found the following.	
[...]
> ? randomprime()=
> {
> p=random;
> while(isprime(p)!=1,p++);
> p
> }
> ? while(ellap(E=ellinit([0,0,0,random,random]),p=randomprime)!=1,)
>   ***   bug in apell1: doubling?, please report

Indeed, there was a missing case. I had mistakenly decided it couldn't occur
but left the assertion behind, just in case. The fix is small, but the patch
is quite big, since I added/modified some comments while I was investigating.

  Karim.

P.S: randomprime() = nextprime(random);
would be more efficient [not that it matters much...].

P.S2: Anybody feels like implementing SEA ?

Index: src/modules/elliptic.c
===================================================================
RCS file: /home/megrez/cvsroot/pari/src/modules/elliptic.c,v
retrieving revision 1.4
diff -c -r1.4 elliptic.c
*** src/modules/elliptic.c	1999/10/01 14:13:46	1.4
--- src/modules/elliptic.c	1999/10/11 14:45:16
***************
*** 1547,1559 ****
      fh = powsell(cp4,f,h,p);
      s = itos(gceil(gsqrt(gdiv(pordmin,bcon),DEFAULTPREC))) >> 1;
      nb = min(128, s >> 1);
      if (bcon == gun)
      { /* first time: initialize */
        tx = newbloc(s+1); *tx = evaltyp(t_VECSMALL) | evallg(s+1);
        ty = newbloc(s+1);
        ti = newbloc(s+1);
      }
!     else f = powsell(cp4,f,bcon,p);
      if (!fh) goto FOUND;
  
      p1 = gcopy(fh);
--- 1547,1560 ----
      fh = powsell(cp4,f,h,p);
      s = itos(gceil(gsqrt(gdiv(pordmin,bcon),DEFAULTPREC))) >> 1;
      nb = min(128, s >> 1);
+     /* look for h s.t f^h = 0 */
      if (bcon == gun)
      { /* first time: initialize */
        tx = newbloc(s+1); *tx = evaltyp(t_VECSMALL) | evallg(s+1);
        ty = newbloc(s+1);
        ti = newbloc(s+1);
      }
!     else f = powsell(cp4,f,bcon,p); /* F */
      if (!fh) goto FOUND;
  
      p1 = gcopy(fh);
***************
*** 1561,1596 ****
      j = lgefint(p);
      for (i=1; i<=nb; i++)
      { /* baby steps */
!       pts[i] = (long)p1;
        _fix(p1+1, j); tx[i] = _low((GEN)p1[1]);
        _fix(p1+2, j); ty[i] = _low((GEN)p1[2]);
!       p1 = addsell(cp4,p1,f,p); /* f^h * F^nb */
        if (!p1) { h = addii(h, mulsi(i,bcon)); goto FOUND; }
      }
      mfh = dummycopy(fh);
      mfh[2] = lnegi((GEN)mfh[2]);
!     fg = addsell(cp4,p1,mfh,p); /* F^nb */
!     if (!fg) { h = addii(h, mulsi(nb,bcon)); goto FOUND; }
      u = cgetg(nb+1, t_VEC);
      av2 = avma; /* more baby steps, nb points at a time */
      while (i <= s)
      {
        long maxj;
!       for (j=1; j<=nb; j++)
        {
!         p1 = (GEN)pts[j];
          u[j] = lsubii((GEN)fg[1], (GEN)p1[1]);
!         if (u[j] == zero)
          {
!           if (!signe(p1[2]) || !egalii((GEN)p1[2],(GEN)fg[2]))
!             { h = addii(h, mulsi(i+j-1,bcon)); goto FOUND; }
!           /* doubling never occurs */
!           err(bugparier, "apell1: doubling?");
          }
        }
        v = multi_invmod(u, p);
        maxj = (i-1 + nb <= s)? nb: s % nb;
!       for (j=1; j<=maxj; j++,i++)
        {
          p1 = (GEN)pts[j];
          addsell_part2(cp4,p1,fg,p, (GEN)v[j]);
--- 1562,1597 ----
      j = lgefint(p);
      for (i=1; i<=nb; i++)
      { /* baby steps */
!       pts[i] = (long)p1; /* h.f + (i-1).F */
        _fix(p1+1, j); tx[i] = _low((GEN)p1[1]);
        _fix(p1+2, j); ty[i] = _low((GEN)p1[2]);
!       p1 = addsell(cp4,p1,f,p); /* h.f + i.F */
        if (!p1) { h = addii(h, mulsi(i,bcon)); goto FOUND; }
      }
      mfh = dummycopy(fh);
      mfh[2] = lnegi((GEN)mfh[2]);
!     fg = addsell(cp4,p1,mfh,p); /* nb.F */
!     if (!fg) { h = mulsi(nb,bcon); goto FOUND; }
      u = cgetg(nb+1, t_VEC);
      av2 = avma; /* more baby steps, nb points at a time */
      while (i <= s)
      {
        long maxj;
!       for (j=1; j<=nb; j++) /* adding nb.F (part 1) */
        {
!         p1 = (GEN)pts[j]; /* h.f + (i-nb-1+j-1).F */
          u[j] = lsubii((GEN)fg[1], (GEN)p1[1]);
!         if (u[j] == zero) /* sum = 0 or doubling */
          {
!           long k = i+j-2;
!           if (egalii((GEN)p1[2],(GEN)fg[2])) k -= 2*nb; /* fg = p1 */
!           h = addii(h, mulsi(k,bcon));
!           goto FOUND;
          }
        }
        v = multi_invmod(u, p);
        maxj = (i-1 + nb <= s)? nb: s % nb;
!       for (j=1; j<=maxj; j++,i++) /* adding nb.F (part 2) */
        {
          p1 = (GEN)pts[j];
          addsell_part2(cp4,p1,fg,p, (GEN)v[j]);
***************
*** 1617,1623 ****
      gaffect(fg, (GEN)pts[1]);
      for (j=2; j<=nb; j++) /* pts = first nb multiples of fg */
        gaffect(addsell(cp4,(GEN)pts[j-1],fg,p), (GEN)pts[j]);
!     /* replace fg by fg^nb since we do nb points at a time */
      avma = av2;
      fg = gcopy((GEN)pts[nb]);
      av2 = avma;
--- 1618,1624 ----
      gaffect(fg, (GEN)pts[1]);
      for (j=2; j<=nb; j++) /* pts = first nb multiples of fg */
        gaffect(addsell(cp4,(GEN)pts[j-1],fg,p), (GEN)pts[j]);
!     /* replace fg by nb.fg since we do nb points at a time */
      avma = av2;
      fg = gcopy((GEN)pts[nb]);
      av2 = avma;
***************
*** 1647,1655 ****
              p1 = addsell(cp4, powsell(cp4,f,stoi(j2),p),fh,p);
              if (egalii((GEN)p1[1], (GEN)ftest[1]))
              {
-               h  = addii(h, mulsi(j2,bcon));
                if (egalii((GEN)p1[2], (GEN)ftest[2])) i = -i;
!               h = addii(h, mulsi(s, mulsi(i, bcon)));
                goto FOUND;
              }
            }
--- 1648,1655 ----
              p1 = addsell(cp4, powsell(cp4,f,stoi(j2),p),fh,p);
              if (egalii((GEN)p1[1], (GEN)ftest[1]))
              {
                if (egalii((GEN)p1[2], (GEN)ftest[2])) i = -i;
!               h = addii(h, mulii(addis(mulss(s,i), j2), bcon));
                goto FOUND;
              }
            }
__
Karim Belabas                    email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud             Tel: (00 33) 1 69 15 57 48
F-91405 Orsay (France)           Fax: (00 33) 1 69 15 60 19
--
PARI/GP Home Page: http://hasse.mathematik.tu-muenchen.de/ntsw/pari/