gmp-snap-20040203/mul_fft.c

[Kevin Ryde]

Jason <jasonmoxham at btclick.com> writes:
>
> on athlon linux 2.4.24 gcc-3.2.3
> ./configure --enable-assert --disable-shared
> make
> gcc -Wall -W jay.c .libs/libgmp.a
> ./a.out
>
> we get assertion failure
> mul_fft.c:996: GNU MP assertion failed: (pad_op +pl -pl3)[__i] == 0
>
>
> #include <stdio.h>
> #include <stdlib.h>
> #define WANT_ASSERT 1
> #include "gmp.h"
>
> int	main(void)
> {mpz_t y,z;
>
> mpz_init(y);
> mpz_init(z);
> mpz_set_ui(y,120);
> mpz_mul_2exp(y,y,100000);
> mpz_add_ui(y,y,1);
> mpz_pow_ui(z,y,12);
> return 0;}
>
> can you confirm ?

I think this originates in mpn_fft_add_modF.  If it only happens in
certain configurations then I suspect that's only different mul
thresholds.  It shows up in my --host=none debug build of the
development code (with gcc 3.3).

If I understand all that fft modulo stuff, values are supposed to be
maintained with a[n]<=1 and the lower limbs arbitrary.  The add_modF
code tests for its result r[n] > 1 and reduces, but I think it needs
to do a second reduction sometimes.

Eg. if a and b are both 1,FF,...,FF, then the sum can have r[n]==3,
and a single subtract of 2^n+1 doesn't suffice to get it to r[n]<=1.

The program you posted ends up with such values,

	1 FF FF ... FF ...
	1 00 00 ... A8 ...

which is enough to give r[n]==3.

This looks the same in gmp 4.1.2 and may be a bug there too, though I
haven't tried to make a failing case for there.

I'm thinking to change the "if" to a "while", something like

	r[n] = a[n] + b[n] + mpn_add_n (r, a, b, n);
	while (r[n] > 1)
	  {
	    MPN_DECR_U (r, n+1, CNST_LIMB(1));
	    r[n]--;
	  }

Paul wrote all this code, hopefully he can confirm or refute.