bin/47301: Miscompilation by bundled GCC

[tg%gmplib.org@localhost]

>Number:         47301
>Category:       bin
>Synopsis:       Miscompilation by bundled GCC
>Confidential:   no
>Severity:       serious
>Priority:       low
>Responsible:    bin-bug-people
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Sun Dec 09 16:00:01 +0000 2012
>Originator:     Torbjorn Granlund
>Release:        NetBSD 6.0
>Organization:
KTH, Sweden
>Environment:
NetBSD slug.gmplib.org 6.0 NetBSD 6.0 (GENERIC) vax
>Description:
The current GMP source code is miscompiled by /usr/bin/gcc.
The result is crashes in n-factorial computations.
>How-To-Repeat:
Either get ftp://gmplib.org/pub/snapshot/gmp-5.1.0-RC1.tar.xz
and isolate the problem of oddfac_1.c in mpz_oddfac_1.  Expect a crash
for n <= ODD_FACTORIAL_TABLE_LIMIT.

I have created a smaller test case which I have not executed, but
where it is very obvious that there is a miscompile.

Source:

typedef unsigned long int mp_limb_t;
typedef struct
{
int _mp_alloc;
int _mp_size;
mp_limb_t *_mp_d;
} __mpz_struct;
typedef __mpz_struct mpz_t[1];
typedef mp_limb_t * mp_ptr;
typedef long int mp_size_t;
typedef __mpz_struct *mpz_ptr;
extern const mp_limb_t __gmp_oddfac_table[];
extern const mp_limb_t __gmp_odd2fac_table[];
static mp_limb_t id_to_n (mp_limb_t id) { return id*3+1+(id&1); }
static mp_limb_t n_to_bit (mp_limb_t n) { return ((n-5)|1)/3U; }
static mp_limb_t limb_apprsqrt (mp_limb_t x) { return x; }
static void
mpz_2multiswing_1 (mpz_ptr x, mp_limb_t n, mp_ptr sieve, mp_ptr factors)
{
  mp_limb_t prod, max_prod;
  mp_size_t j;
  j = 0;
  prod = -(n & 1);
  n &= ~1L;
  prod = (prod & n) + 1;
  max_prod = ((~ ((mp_limb_t) (0))) >> 0) / (n-1);
  mp_limb_t __q, __prime;
  __prime = 3;
  if (prod > max_prod)
    {
      factors[j++] = prod;
      prod = 1;
    }
  __q = n;
  do {
    __q /= __prime;
    if ((__q & 1) != 0)
      prod *= __prime;
  } while (__q >= __prime);

  mp_limb_t s;
  mp_limb_t prime;
  s = limb_apprsqrt(n);
  s = n_to_bit (s);
  mp_limb_t __mask, __index, __max_i, __i;
  __i = n_to_bit (5);
  __index = __i / 32;
  __mask = ((mp_limb_t) 1L) << (__i % 32);
  __i += 0;
  __max_i = s;
  do {
    ++__i;
    if ((sieve[__index] & __mask) == 0)
      {
        prime = id_to_n(__i);
        mp_limb_t __q, __prime;
        __prime = prime;
        if (prod > max_prod)
          {
            factors[j++] = prod;
            prod = 1;
          }
        __q = n;
        do {
            __q /= __prime;
            if ((__q & 1) != 0)
              prod *= __prime;
        } while (__q >= __prime);
      }
    __mask = __mask << 1 | __mask >> 31;
    __index += __mask & 1;
  } while (__i <= __max_i);
  s++;

  if (j != 0)
    {
      factors[j++] = prod;
      __gmpz_prodlimbs (x, factors, j);
    }
}
void
miscompiled_function (mpz_ptr x, mp_limb_t n, unsigned flag)
{
  if (n <= 16)
    {
      /* This is miscompiled.  The store address is completely bogus.  */
      mp_limb_t xxx = __gmp_oddfac_table[n];
      (x->_mp_d)[0] = xxx;
      (x->_mp_size) = 1;
    }
  else
    {
      unsigned s;
      mp_ptr factors;
      mp_limb_t tn;
      mp_limb_t prod, max_prod, i;
      mp_size_t j;
      s = 0;
      for (tn = n; tn >= 400; s++)
        tn >>= 1;
      j = 0;
      factors = (mp_limb_t *) __builtin_alloca(199 * sizeof (mp_limb_t));
      prod = 1;
      max_prod = 98723491;
      do {
        i = 21;
        factors[j++] = 0x27065f73L;
        do {
          if (prod > max_prod)
            {
              factors[j++] = prod;
              prod = i;
            }
          else
            prod *= i;
          i += 2;
        } while (i <= tn);
        max_prod <<= 1;
        tn >>= 1;
      } while (tn > 20);
      factors[j++] = prod;
      factors[j++] = __gmp_odd2fac_table[(tn - 1) >> 1];
      factors[j++] = __gmp_oddfac_table[tn >> 1];
      __gmpz_prodlimbs (x, factors, j);

      if (s != 0)
        {
          mpz_t mswing;
          mp_ptr sieve;
          mp_size_t size;
          size = n / (32 - 0) + 4;
          mpz_ptr __x = mswing;
          __x->_mp_alloc = size;
          __x->_mp_d = (mp_limb_t *) __builtin_alloca(size * sizeof 
(mp_limb_t));
          sieve = (mswing->_mp_d) + size / 2 + 1;
          size = (__gmp_primesieve (sieve, n - 1) + 1) / n + 1;
          factors = (mp_limb_t *) __builtin_alloca(size * sizeof (mp_limb_t));
          do
            {
              mp_ptr square, px;
              mp_size_t nx, ns;
              mp_limb_t cy;
              s--;
              mpz_2multiswing_1 (mswing, n >> s, sieve, factors);
              nx = (x->_mp_size);
              ns = (mswing->_mp_size);
              nx = size + ns;
              px = (nx) > (x->_mp_alloc) != 0 ? (mp_ptr) __gmpz_realloc(x,nx) : 
x->_mp_d;
              cy = __gmpn_mul (px, square, size, mswing->_mp_d, ns);
              (x->_mp_size) = nx - (cy == 0);
            }
          while (s != 0);
        }
    }
}

Compile with the default compiler (claimed as 4.1.3) using -O
and no other flags.  This code is generated for the first part
of miscompiled_function:

miscompiled_function:
        .word 0xfc0
        subl2 $60,%sp
        cmpl 8(%ap),$16
        jlequ .L4
        cmpl 8(%ap),$399
        jgtru .L6
        movl 8(%ap),%r4
        clrl -60(%fp)
        jbr .L8
.L4:
        movl 8(%ap),%r0
        movl __gmp_oddfac_table[%r0],*8(%r0)               BAD
        movl 4(%ap),%r1
        movl $1,4(%r1)
        ret

The line annotated BAD stores in a completely invalid place,
related to the (non-pointer!) argument n.  Note that the source
operand is valud, using n properly.

It is fairly easy to make the compiler generate valid code by
making innocent changes to this file.

>Fix:
I have not made any attempt at fixing this compiler bug, since
I feel gcc 4.1 is obsolete.

If I were you, I'd upgrade to a current GCC version.  I notice
that your fear of GPL 3 seems limited given that the binutils
of this install is quite new, so GPL 3 must not be a valid
reason for staying obsolete with GCC.