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
20gcdex(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
73igcdex(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   )$
89