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