Re: Bug in arithmetic?

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

[Francois Morain:]
> o.osf1-alpha/gp-sta 
>                    GP/PARI CALCULATOR Version 2.0.13 (alpha)
>                               Alpha 64-bit version
>                 (readline disabled, extended help not available)
> 
>                            Copyright (C) 1989-1998 by
>           C. Batut, K. Belabas, D. Bernardi, H. Cohen and M. Olivier.
> 
> Type ? for help, \q to quit.
> Type ?12 for how to get moral (and possibly technical) support.
> 
>    realprecision = 38 significant digits
>    seriesprecision = 16 significant terms
>    format = g0.38
> 
> parisize = 8000000, primelimit = 500000
> ? n=10^2000+4561;
> 
> For n, we should have 2^(n-1) = 1 mod n, but the answer is not:

This was corrected in patch 1 (due to a mistake in a stupid script, I only
numbered the first 2 patches in the 2.0.13 series...), msg 351 on this list.
I'm reposting it:

[me:]
> This patch corrects a nasty hidden bug in Karatsuba integer multiplication
> (could occur when both operands have more than 155 digits and one of them
> contains many consecutive 0 words...)

*** src/kernel/none/mp.c.orig   Tue Dec 15 16:30:00 1998
--- src/kernel/none/mp.c        Thu Dec 17 17:12:20 1998
***************
*** 1833,1853 ****
    GEN z,z0,y0,yd, zd = (GEN)avma;
    long a,lz,ly = lgefint(y);
  
!   (void)new_chunk(d);
!   if (ly!=2) /* y != 0 */
    {
!     yd = y+ly; y0 = yd-d;
!     while (yd > y0) *--zd = *--yd; /* copy last d words of y */
    }
!   a = ly-2 - d;
!   if (a <= 0)
    {
!     if (a) { z0 = (GEN)avma; while (zd >= z0) *--zd = 0; }
      z = icopy(x);
    }
!   else
!     z = addiispec(x+2, y+2, lgefint(x)-2, a);
!   lz=lgefint(z)+d;
    z[1] = evalsigne(1) | evallgefint(lz);
    z[0] = evaltyp(t_INT) | evallg(lz); return z;
  }
--- 1833,1856 ----
    GEN z,z0,y0,yd, zd = (GEN)avma;
    long a,lz,ly = lgefint(y);
  
!   z0 = new_chunk(d);
!   a = ly-2; yd = y+ly; 
!   if (a >= d)
    {
!     y0 = yd-d; while (yd > y0) *--zd = *--yd; /* copy last d words of y */
!     a -= d;
!     if (a)
!       z = addiispec(x+2, y+2, lgefint(x)-2, a);
!     else
!       z = icopy(x);
    }
!   else
    {
!     y0 = yd-a; while (yd > y0) *--zd = *--yd; /* copy last a words of y */
!     while (zd >= z0) *--zd = 0;    /* complete with 0s */
      z = icopy(x);
    }
!   lz = lgefint(z)+d;
    z[1] = evalsigne(1) | evallgefint(lz);
    z[0] = evaltyp(t_INT) | evallg(lz); return z;
  }
--
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://pari.home.ml.org