read recden; read PGCD; p := 2: alias(alpha = RootOf(x^2+x+1)): X := [alpha,x,y]: g := expand(randpoly(X,degree=5,dense)) mod p: a := expand(randpoly(X,degree=10,dense)*g) mod p: b := expand(randpoly(X,degree=10,dense)*g) mod p: X := [x,y,w]: A := convert(a,POLYNOMIAL,X,w=alpha,p); B := convert(b,POLYNOMIAL,X,w=alpha,p); G := PGCD(A,B); divrpoly(A,G); divrpoly(B,G); A,B,C,Z := 'A,B,C,Z': z := (A^2*B^2*Z^2+A^2*B^2+B^2*Z^2+Z^2+B^2)/Z: a := (B^2*C^2*A^2+B^2*C^2+C^2*A^2+A^2+C^2)/A: b := (C^2*Z*B^2+C^2*Z+Z*B^2+B^2+Z)/B: c := (Z*A*C^2+Z*A+A*C^2+C^2+A)/C: f := Z^2*A^3*B^2*C^2*(z^2*b^2*c^2+(a+z+1)^2*(a+c^2+1)): p := expand(f) mod 2: q := diff(p,A) mod 2: X := [A,B,C,Z]: P := convert(p,POLYNOMIAL,X,2); Q := convert(q,POLYNOMIAL,X,2); G := PGCD(P,Q);