(*:Name: Algebra`PolynomialExtendedGCD` *)

(*:Author: Matthew Markert *)

(*:Context: Algebra`PolynomialExtendedGCD` *)

(*:Package Version: 1.1 *)

(*:Copyright: Copyright 1993-1999,  Wolfram Research, Inc.*)

(*:History: Version 1.0 by Matthew Markert, 1993
	    Version 1.1 by ECM, 1998 *)

(*:Keywords: gcd, finite fields, Euclidean *)

(*:Requirements: none. *)

(*:Warnings: none. *)

(*:Source: any graduate-level introductory algebra text. *)

(*:Limitations: This function only works for monovariate polynomials. *)

(*:Summary: 
This package defines the function PolynomialExtendedGCD.  This works for
univariate polynomials with coefficients in the rationals, the Gaussian 
rationals, or in a finite field if the finite fields package is loaded.
*)

(*:Discussion: 
	The GCD is only properly defined up to a factor of any field 
	element, so we make the GCD a monic polynomial (i.e., having a leading
	nonzero coefficient of 1) where possible.
   
	The kernel has GCD and ExtendedGCD for integers, and it has 
	PolynomialGCD for multivariate polynomials.  Besides the options 
	used here, one could also deal with polynomials over the over the 
	reals or complexes.  Perhaps GCD should have a general mechanism 
	for determining the base field. 
	
	This package uses the Euclidean algorithm.
*)

BeginPackage["Algebra`PolynomialExtendedGCD`"]

(* Usage Messages *)

PolynomialExtendedGCD::usage =
  "PolynomialExtendedGCD[poly1, poly2] gives the extended greatest
common divisor of two monovariate polynomials over the rationals.
PolynomialExtendedGCD[poly1, poly2, Modulus->p] gives the extended
GCD over the integers mod prime p.";

polyextGCD::usage = "debug"

Begin["`Private`"]
Unprotect[PolynomialExtendedGCD]

(* PolynomialExtendedGCD for polynomials over various fields *)

(* Note that Options[PolynomialGCD] = 
	{Extension -> None, Modulus -> 0, Trig -> False}.
*)
Options[PolynomialExtendedGCD] = {Modulus -> 0}


(* Messages *)

PolynomialExtendedGCD::onevar =
  "PolynomialExtendedGCD is defined only for polynomials in one variable.";

PolynomialExtendedGCD::pegcdz =
  "PolynomialExtendedGCD[0, 0] has no unique solution"

PolynomialExtendedGCD::pegcdmz =
  "PolynomialExtendedGCD[0, 0, Modulus -> ``] has no unique solution"

PolynomialExtendedGCD::modp =
  "Value of option Modulus -> `` should be a prime number or 0."

(* ***** integers mod a prime ***** *)

PolynomialExtendedGCD[f_, g_, opts___?OptionQ] :=
  Module[ {fmp, gmp, y, temp, result, inputok=True},
   (
        result
   ) /;
    (p = Modulus /. {opts} /. Options[PolynomialExtendedGCD];
     p =!= 0 &&
     If[!PrimeQ[p],
	Message[PolynomialExtendedGCD::modp, p];
	False,
	True]) &&
    FreeQ[fmp=PolynomialMod[f,p], PolynomialMod] &&
    FreeQ[gmp=PolynomialMod[g,p], PolynomialMod] &&
    (
    	 y = Union[Variables[fmp],Variables[gmp]];
         If[Length[y] > 1,
            Message[PolynomialExtendedGCD::onevar];
            inputok = False;
         ];
         If[ Length[y] === 0,

            (* =========== arguments are polynomials of order 0 ============ *)
            If[fmp == gmp == 0,
               Message[PolynomialExtendedGCD::pegcdmz, p];
               inputok = False,
               (* fmp and gmp are not both 0 *)
               temp = extendedGCD[fmp,gmp];
               If[ Head[temp] === List,
                  (* divide list by gcd (mod p), because when working with
		     univariate polynomials we want to return a
		     monic polynomial (i.e., leading nonzero coeff is 1) *)
                  result = Mod[PowerMod[temp[[1]], -1, p] * temp, p],
                  (* unevaluated extendedGCD *)
                  inputok = False
               ]
            ],

            (* else Length[y] === 1 *)
            (* =========== arguments are polynomials of order 1 ============ *)
            y = First[y];
            result = If[ Exponent[fmp,y]<Exponent[gmp,y],
                   	temp = polyextGCD[gmp,fmp,y,p];
                   	{temp[[1]], Reverse[ temp[[2]] ]},
                   	polyextGCD[fmp,gmp,y,p]
                     ]
         ]; (* end Length[y] === 0 *)
         inputok
    )
  ]



(* ***** Rationals and finite fields ***** *)

PolynomialExtendedGCD[f_, g_, opts___?OptionQ] :=
Module[{fmp = f, gmp = g, temp, y, result, inputok=True},
  (
        result
  ) /; 
       (
       p = Modulus /. {opts} /. Options[PolynomialExtendedGCD];
       p === 0
       ) &&
       (
    (*** A finite field element is not a variable! ***)
    y = Select[
      Union[Variables[fmp],Variables[gmp]],
      (Head[Head[#]]=!=Algebra`FiniteFields`GF)&];
    If[ Length[y] > 1,
      Message[PolynomialExtendedGCD::onevar];
      inputok = False;
    ];
    If[ Length[y] === 0,

      (* =========== arguments are polynomials of order 0 =============== *)
      If[fmp == gmp == 0,
            Message[PolynomialExtendedGCD::pegcdz];
            inputok = False,
            (* fmp and gmp are not both 0 *)
            temp = extendedGCD[fmp,gmp];
            If[ Head[temp] === List,
                (* divide list by gcd (mod p), because when working with
		   univariate polynomials we want to return a
		   monic polynomial (i.e., leading nonzero coeff is 1) *)
                result = temp/temp[[1]],
                (* unevaluated extendedGCD *)
                inputok = False
            ]
      ],

      (* =========== arguments are polynomials of order 1 =============== *)
      (* else Length[y] === 1 *)
      y = First[y];
      If[ Exponent[fmp,y] < Exponent[gmp,y],
        temp = polyextGCD[gmp,fmp,y];
        result = {temp[[1]], Reverse[ temp[[2]] ]},
        result = polyextGCD[fmp,gmp,y]
      ]
    ];
    inputok
       )
] /; (Precision[{f,g}] === Infinity)



(* =================================================================== *)
(* 
The function polyextGCD strongly assumes its input has been error-checked
by PolynomialExtendedGCD.  The degree of f must not be greater than the 
degree of g. The degree of g must be at least 1. 
*)

polyextGCD[f_, 0, y_, p_] :=
  Module[
    {monic = PowerMod[Coefficient[f, y, Exponent[f, y]], -1, p]},
    {PolynomialMod[monic f, p], Mod[monic {1,0}, p]}
  ]

polyextGCD[f_, g_, y_, p_] :=
  Module[
    { fm = f, gm = g, fv = {1, 0}, gv = {0, 1},
      q, r = y, leading, monic},
    While[Exponent[r, y] > 0,
      leading = Coefficient[gm, y, Exponent[gm, y]];
      monic = PowerMod[leading, -1, p];
      q = PolynomialQuotient[
        PolynomialMod[monic fm, p],
        PolynomialMod[monic gm, p], y];
      r = leading *
        PolynomialRemainder[
          PolynomialMod[monic fm, p],
          PolynomialMod[monic gm, p], y];
      {fm, gm, fv, gv}  =
        PolynomialMod[#, p] & /@ {gm, r, gv, fv - q gv};
    ];
    If[gm === 0,
      monic = PowerMod[Coefficient[fm, y, Exponent[fm, y]], -1, p];
      {PolynomialMod[monic fm,p], PolynomialMod[monic fv, p]},
      {1, PowerMod[gm, -1, p] gv}
    ]
  ] /; (f =!= 0)

polyextGCD[ f_, 0, y_] :=
  (
  {f,{1,0}}
  )


polyextGCD[f_, g_, y_] :=
Module[
  { fm = f, gm = g, fv = {1, 0}, gv = {0, 1},
    q, r = y, leading},
  While[Exponent[r, y] > 0,
    q = PolynomialQuotient[fm, gm, y];
    r = PolynomialRemainder[fm, gm, y];
    {fm, gm, fv, gv}  = {gm, r, gv, fv - q gv};
  ];
  If[gm === 0,
    leading = Coefficient[fm, y, Exponent[fm, y]];
    {Expand[fm/leading], Expand[fv/leading]},
    {1, Expand[gv/gm]}
  ]
] /; (f =!= 0)

(* =================================================================== *)

(* Assuming that not both fmp and gmp are zero, this gives
	ExtendedGCD for rational arguments.
 *)

extendedGCD[fmp_?IntegerQ, gmp_?IntegerQ] := ExtendedGCD[fmp, gmp] 

extendedGCD[fmp_, gmp_] :=
   Module[{gcd = GCD[fmp, gmp]},
         If[fmp != 0,
               {gcd, {gcd/fmp, 0}},
               {gcd, {0, gcd/gmp}}
         ]
   ]




End[] (* `Private` context *)

Protect[PolynomialExtendedGCD]

EndPackage[]   (* Algebra`PolynomialExtendedGCD` *)
