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