#5659 The default gcd algorithm can give wrong results or hang, especially with a small modulus

[David Scherfgen]

By Claude:

The default gcd algorithm (gcd: spmod) can return wrong results, signal errors or hang. It works modulo a prime at random evaluation points. With a small modulus the result changes from call to call, and without a modulus it goes wrong for some inputs whose coefficients involve its first prime, 2305843009213693921 in SBCL:

(%i1) display2d: false$
(%i2) modulus: 3$
(%i3) f: expand((x+1)*(y+1))$
(%i4) g: expand((x+1)*(x*y^2+1))$
(%i5) makelist(gcd(f, g), 12);
(%o5) [x+1,x+1,(x+1)*y+x+1,x+1,x+1,(x+1)*y+x+1,x+1,(x+1)*y+x+1,x+1,
       (x+1)*y+x+1,x+1,(x+1)*y+x+1]
(%i6) modulus: false$
(%i7) P: 2305843009213693921$
(%i8) gcd(x^2+1, x^2+1+P);
(%o8) x^2+2305843009213693922
(%i9) gcd(expand((P*x^2+1)*(x+2)), expand((P*x^2+1)*(x+3)));
(%o9) 1

(x+1)*y+x+1 is f itself, which does not divide g. The correct result is x+1, and the last two results should be 1 and 2305843009213693921*x^2+1.

With modulus: 2 the following returns 1 instead of x*y+x+y, and the third identical call never returns:

(%i1) display2d: false$
(%i2) modulus: 2$
(%i3) f: expand((x+y)*(x*y+x+y))$
(%i4) g: expand(x*(x*y+x+y))$
(%i5) gcd(f, g);
(%o5) 1
(%i6) gcd(f, g);
(%o6) 1
(%i7) gcd(f, g);

With modulus: 13, gcd(expand(5*u*(-6*u-6)), expand(5*u*(2*z^2+u*(-6*z-4*y)-6))) sometimes signals PTPTQUOTIENT: Polynomial quotient is not exact instead of returning u.

Setting gcd: subres gives the correct results in all these examples.