#3832 Bug in spmod gcd algorithm for polynomials

[Marius Gerbershagen]

The following input

trigsimp(trigexpand(1/(16*(a^8*m^4*cos(theta)^4+2*a^6*m^4*r^2*cos(theta)^2+a^4*m^4*r^4)*sin(theta)^8+8*(a^10*m^2*r^2*cos(theta)^10+2*(2*a^8*m^2*r^4-a^8*m^3*r^2)*cos(theta)^8+6*(a^6*m^2*r^6-a^6*m^3*r^4)*cos(theta)^6+2*(2*a^4*m^2*r^8-3*a^4*m^3*r^6)*cos(theta)^4+(a^2*m^2*r^10-2*a^2*m^3*r^8)*cos(theta)^2)*sin(theta)^4+a^12*r^4*cos(theta)^16+2*(3*a^10*r^6-2*a^10*m*r^4)*cos(theta)^14+(15*a^8*r^8-20*a^8*m*r^6+4*a^8*m^2*r^4)*cos(theta)^12+4*(5*a^6*r^10-10*a^6*m*r^8+4*a^6*m^2*r^6)*cos(theta)^10+(15*a^4*r^12-40*a^4*m*r^10+24*a^4*m^2*r^8)*cos(theta)^8+2*(3*a^2*r^14-10*a^2*m*r^12+8*a^2*m^2*r^10)*cos(theta)^6+(r^16-4*m*r^14+4*m^2*r^12)*cos(theta)^4)+1/(4*((a^6*m^2*r^2-2*a^6*m^3+a^8*m^2)*cos(theta)^4+2*(a^4*m^2*r^4+(a^6*m^2-2*a^4*m^3)*r^2)*cos(theta)^2+a^2*m^2*r^6+(a^4*m^2-2*a^2*m^3)*r^4)*sin(theta)^4+(a^8*r^4+(a^10-2*a^8*m)*r^2)*cos(theta)^10+2*(2*a^6*r^6+(2*a^8-5*a^6*m)*r^4-(a^8*m-2*a^6*m^2)*r^2)*cos(theta)^8+6*(a^4*r^8+(a^6-3*a^4*m)*r^6-(a^6*m-2*a^4*m^2)*r^4)*cos(theta)^6+2*(2*a^2*r^10+(2*a^4-7*a^2*m)*r^8-3*(a^4*m-2*a^2*m^2)*r^6)*cos(theta)^4+(r^12+(a^2-4*m)*r^10-2*(a^2*m-2*m^2)*r^8)*cos(theta)^2)));

gives an error

PQUOTIENT: Quotient by a polynomial of higher degree (case 2a)
 -- an error. To debug this try: debugmode(true);

I have traced this down to a bug in the zgcd function which computes the gcd of

sin(theta)^4*(cos(theta)^4*(r^2*m^2*a^6*4+(m^3*a^6*(-8)+m^2*a^8*4))+cos(theta)^2*(r^4*m^2*a^4*8+r^2*(m^3*a^4*(-16)+m^2*a^6*8))+(r^6*m^2*a^2*4+r^4*(m^3*a^2*(-8)+m^2*a^4*4)))+(cos(theta)^10*(r^4*a^8*1+r^2*(m*a^8*(-2)+a^10*1))+cos(theta)^8*(r^6*a^6*4+r^4*(m*a^6*(-10)+a^8*4)+r^2*(m^2*a^6*4+m*a^8*(-2)))+cos(theta)^6*(r^8*a^4*6+r^6*(m*a^4*(-18)+a^6*6)+r^4*(m^2*a^4*12+m*a^6*(-6)))+cos(theta)^4*(r^10*a^2*4+r^8*(m*a^2*(-14)+a^4*4)+r^6*(m^2*a^2*12+m*a^4*(-6)))+cos(theta)^2*(r^12*1+r^10*(m*(-4)+a^2*1)+r^8*(m^2*4+m*a^2*(-2))))

and

sin(theta)^8*(cos(theta)^4*m^4*a^8*16+cos(theta)^2*r^2*m^4*a^6*32+r^4*m^4*a^4*16)+sin(theta)^4*(cos(theta)^10*r^2*m^2*a^10*8+cos(theta)^8*(r^4*m^2*a^8*32+r^2*m^3*a^8*(-16))+cos(theta)^6*(r^6*m^2*a^6*48+r^4*m^3*a^6*(-48))+cos(theta)^4*(r^8*m^2*a^4*32+r^6*m^3*a^4*(-48))+cos(theta)^2*(r^10*m^2*a^2*8+r^8*m^3*a^2*(-16)))+(cos(theta)^16*r^4*a^12*1+cos(theta)^14*(r^6*a^10*6+r^4*m*a^10*(-4))+cos(theta)^12*(r^8*a^8*15+r^6*m*a^8*(-20)+r^4*m^2*a^8*4)+cos(theta)^10*(r^10*a^6*20+r^8*m*a^6*(-40)+r^6*m^2*a^6*16)+cos(theta)^8*(r^12*a^4*15+r^10*m*a^4*(-40)+r^8*m^2*a^4*24)+cos(theta)^6*(r^14*a^2*6+r^12*m*a^2*(-20)+r^10*m^2*a^2*16)+cos(theta)^4*(r^16*1+r^14*m*(-4)+r^12*m^2*4))

to be

sin(theta)^4*(cos(theta)^4*(r^2*m^2*a^6*4+(m^3*a^6*(-8)+m^2*a^8*4))+cos(theta)^2*(r^4*m^2*a^4*8+r^2*(m^3*a^4*(-16)+m^2*a^6*8))+(r^6*m^2*a^2*4+r^4*(m^3*a^2*(-8)+m^2*a^4*4)))+(cos(theta)^10*(r^4*a^8*1+r^2*(m*a^8*(-2)+a^10*1))+cos(theta)^8*(r^6*a^6*4+r^4*(m*a^6*(-10)+a^8*4)+r^2*(m^2*a^6*4+m*a^8*(-2)))+cos(theta)^6*(r^8*a^4*6+r^6*(m*a^4*(-18)+a^6*6)+r^4*(m^2*a^4*12+m*a^6*(-6)))+cos(theta)^4*(r^10*a^2*4+r^8*(m*a^2*(-14)+a^4*4)+r^6*(m^2*a^2*12+m*a^4*(-6)))+cos(theta)^2*(r^12*1+r^10*(m*(-4)+a^2*1)+r^8*(m^2*4+m*a^2*(-2))))

Other gcd algorithms seem to work fine, for instance the oldgcd function gives

sin(theta)^4*(cos(theta)^4*m^2*a^6*4+cos(theta)^2*r^2*m^2*a^4*8+r^4*m^2*a^2*4)+(cos(theta)^10*r^2*a^8*1+cos(theta)^8*(r^4*a^6*4+r^2*m*a^6*(-2))+cos(theta)^6*(r^6*a^4*6+r^4*m*a^4*(-6))+cos(theta)^4*(r^8*a^2*4+r^6*m*a^2*(-6))+cos(theta)^2*(r^10*1+r^8*m*(-2)))

Unfortunately, I haven't been able to reduce the example to a simpler one.

[Robert Dodier]

This bug appears to be fixed by commit 73105cf (from merge request #25), but before closing this bug report, I'd like to put together some test cases. If anyone has some test cases, they could post them here, otherwise, I'll put something together.