share/algebra/gcdex.mac

   1 /* Implementation of the Maxima user functions gcdex and igcdex
   2
   3    gcdex(f,g) yields a list [a,b,u]
   4    where u is the GCD of f and g, and u = f*a+g*b.
   5
   6    F and G should be univariate polynomials, or a main variable
   7    should be supplied so that we are working in the polynomial ring over
   8    the rational functions in the other variables, and thus in a
   9    principal ideal domain, and the ideal generated by (f,g) is indeed
  10    generated by the gcd of f and g.
  11
  12    igcdex(n,k) yields a list [a,b,u]
  13    where u is the GCD of the integers n and k, and u = n*a+k*b.
  14
  15    igcdex is called from the function gcdex for integer arguments.
  16
  17    This file is part of Maxima -- GPL CAS based on DOE-MACSYMA
  18 */
  19
  20 gcdex(f, g, [var]):=
  21    block([q0:[], q1:[], ok:true, lis1:[], lis2:[], lis3:[], q:[]],
  22       if integerp(f) and integerp(g) then return(igcdex(f, g)),
  23       if var = [] then var:listofvars([f,g]),
  24       if not (length(var) = 1) then
  25       (
  26          if (f=0 and g=0) then return([0,0,0]),
  27          error("Specify one main variable or give univariate polynomials")
  28       ),
  29       var:var[1],
  30       f:rat(f),
  31       g:rat(g),
  32       /* Do not divide by zero, but return the result */
  33       if g = 0 then return([1, 0, f]),
  34       /* divide(f,g) => [q:quotient(f,g), remainder:f-g*q] */
  35       q0:divide(f, g, var),
  36       /* if f/g is 0 then we reverse them */
  37       if first(q0) = 0 then
  38       (  /* the deg f < deg g */
  39          lis2:gcdex(g, f, var),
  40          return([second(lis2), first(lis2), third(lis2)])
  41       ),
  42       if second(q0) = 0 then
  43          lis2:[0, 1, g]
  44       else
  45       (
  46          q1:divide(g, second(q0), var),
  47          lis1:[1, -first(q0), second(q0)],
  48          if second(q1) = 0 then
  49             lis2:lis1
  50          else
  51          (  /* lisi are always perpendicular to [f,g,-1] */
  52             lis2:[-first(q1), 1+first(q0)*first(q1), second(q1)],
  53             while (ok) do
  54             (
  55                q:divide(third(lis1), third(lis2), var),
  56                lis3:lis1-lis2*first(q),
  57                if third(lis3) = 0 then
  58                   ok:false
  59                else
  60                (
  61                   lis1:lis2,
  62                   lis2:lis3/first(content(third(lis3), var))
  63                )
  64             )
  65          )
  66       ),
  67       if freeof(var, third(lis2)) then
  68          lis2/third(lis2)
  69       else
  70          lis2
  71    );
  72
  73 igcdex(a,b) :=
  74    block([x:0, lx:1, y:1, ly:0, q, r],
  75       if not (integerp(a) and integerp(b)) then
  76          error("igcdex: The arguments must be integers.")
  77       else
  78       (
  79          while b # 0 do
  80          (
  81             [q,r]:divide(a,b),
  82             [a,b]:[b,r],
  83             [x,lx]:[lx-q*x,x],
  84             [y,ly]:[ly-q*y,y]
  85          ),
  86          [lx,ly,a]
  87       )
  88    )$