Hello, please include the following snippet. Its application domain depends on that of quo() which is univariate, currently. It is tested but not heavily so. It is not optimized. Thanks, ralf Input: x,y rational polynomials Computes u,v,d such that ux+vy = gcd(x,y). /////////// Extended GCD ex xgcd (const ex& a, const ex& b, ex& u, ex& v, const symbol& s) { const ex& x = a.normal(); const ex& y = b.normal(); if (is_exactly_a<numeric>(x) && is_exactly_a<numeric>(y) && x.info(info_flags::integer) && y.info(info_flags::integer)) { cln::cl_I cu, cv, cr; cr = cln::xgcd (cln::the<cln::cl_I>(ex_to<numeric>(x).to_cl_N()), cln::the<cln::cl_I>(ex_to<numeric>(y).to_cl_N()), &cu, &cv); u = numeric(cu); v = numeric(cv); return numeric(cr); } if (is_exactly_a<numeric>(y) && y.info(info_flags::integer)) { u = 0; v = numeric(1)/y; return 1; } if (is_exactly_a<numeric>(x) && x.info(info_flags::integer)) { v = 0; u = numeric(1)/x; return 1; } if (x.info(info_flags::rational_polynomial) && y.info(info_flags::rational_polynomial)) { ex d, p, q; d = gcd (x, y); int dx = degree(x, s), dy = degree(y, s); if (dx < dy) { p = y; q = x; } else { p = x; q = y; } ex u1,u2,u3,v1,v2,v3,w1,w2,w3; u1 = p; u2 = 1; u3 = 0; v1 = q; v2 = 0; v3 = 1; while (!v1.is_zero()) { ex Q = quo (u1, v1, s); w1 = u1 - Q*v1; w2 = u2 - Q*v2; w3 = u3 - Q*v3; u1 = v1.expand(); u2 = v2.expand(); u3 = v3.expand(); v1 = w1.expand(); v2 = w2.expand(); v3 = w3.expand(); } u = (u2 / ((u1/d).normal())).expand(); v = (u3 / ((u1/d).normal())).expand(); //cerr << "CHK: " << (u*x+v*y).expand() << endl << flush; return d; } throw invalid_argument("xgcd(): not a rational polynomial"); }
Hi! On Sun, 3 Oct 2004, Ralf Stephan wrote:
please include the following snippet. Its application domain depends on that of quo() which is univariate, currently. It is tested but not heavily so. It is not optimized.
Thanks, ralf
Input: x,y rational polynomials Computes u,v,d such that ux+vy = gcd(x,y).
/////////// Extended GCD ex xgcd (const ex& a, const ex& b, ex& u, ex& v, const symbol& s)
What's the point in having this twice? In normal.cpp:1275 we have: /** Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X) * and b(X) in Z[X]. * * @param a first multivariate polynomial * @param b second multivariate polynomial * @param ca pointer to expression that will receive the cofactor of a, or NULL * @param cb pointer to expression that will receive the cofactor of b, or NULL * @param check_args check whether a and b are polynomials with rational * coefficients (defaults to "true") * @return the GCD as a new expression */ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) { // ... Is the root of the problem the missing documentation of cofactors in the tutorial or am I missing something fundamental?! Confused -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Is the root of the problem the missing documentation of cofactors in the tutorial or am I missing something fundamental?!
Actually, yes, both, sorry. This code: ex e1, e2; symbol x("x"); cout << gcd(x-1, x*x*x, &e1, &e2, true) << "," << e1 << "," << e2 << endl; gives 1,-1+x,x^3 where my function and Pari's bezout() gives 1,-1-x-x^2,1 (or a different order but you get the idea) since -(x-1)(1+x+x^2) + x^3 = 1. Best regards, ralf
Hello, In case you did not understand my last message on this, here again: You wrote
On Sun, 3 Oct 2004, Ralf Stephan wrote:
ex xgcd (const ex& a, const ex& b, ex& u, ex& v, const symbol& s)
What's the point in having this twice? In normal.cpp:1275 we have:
/** Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X) * and b(X) in Z[X]. * @param ca pointer to expression that will receive the cofactor of a, or NULL * @param cb pointer to expression that will receive the cofactor of b, or NULL
To declare something as cofactor does not mean it is. What the function returns does NOT satisfy the usual cofactor equation gcd(a,b) = ca*a + cb*b. Vide my earlier example, or a = (x-1)(x+1), b = (x-1)^3 which should give ca = (3-x)/4 and cb = 1/4. But it gives ca = (1-x)^2, cb = 1+x. Of course, if your gcd() WOULD return the cofactors, normal() would be broken, as it depends on the odd behaviour which I have no explanation for. It would be likewise helpful for me if you just tell your intention to look into this or not, so I can either go on or try to fix your code. Certainly, you are more familiar with the code. Best, ralf
Hi! On Mon, Nov 08, 2004 at 12:45:41PM +0100, Ralf Stephan wrote:
To declare something as cofactor does not mean it is. What the function returns does NOT satisfy the usual cofactor equation
gcd(a,b) = ca*a + cb*b.
These are not the cofactors (at least, I've never seen the quantities from this equation referred to as "cofactors"). The cofactors in GiNaC's gcd() are defined by a = ca*gcd(a,b) -> ca = a/gcd(a,b) b = cb*gcd(a,b) -> cb = b/gcd(a,b) Like many other things in GiNaC, this follows the behavior of Maple. What you want (the solutions to s,t in gcd(a,b)=s*a+t*b) is indeed something entirely different, and computed by a different algorithm (gcdex() in Maple).
Vide my earlier example, or a = (x-1)(x+1), b = (x-1)^3 which should give ca = (3-x)/4 and cb = 1/4. But it gives ca = (1-x)^2, cb = 1+x.
(It actually gives ca=1+x, cb=(1-x)^2.) Bye, Christian -- / Physics is an algorithm \/ http://www.uni-mainz.de/~bauec002/
Hello,
gcd(a,b) = ca*a + cb*b.
These are not the cofactors (at least, I've never seen the quantities from this equation referred to as "cofactors").
Heh. See C. Yap, Algorithmic Algebra, Ch. 2, §2, p. 5 which is online: http://www.cs.nyu.edu/yap/book/ He agrees with you in that the algorithm to compute them is handled separately from the treatise of gcd computation. Voilà, extgcd. Now, we are back at square one where I just sent the xgcd() function to the list. But it had a bug. If you want to have the corrected version in GiNaC, I will send it gladly and directly to you. ralf
Hi! On Mon, Nov 08, 2004 at 04:15:54PM +0100, Ralf Stephan wrote:
Now, we are back at square one where I just sent the xgcd() function to the list. But it had a bug. If you want to have the corrected version in GiNaC, I will send it gladly and directly to you.
Yep, keep 'em coming. What I find a bit suspicious about your earlier xgcd() is its generous use of normal(). Maybe it would be better to use multiply_lcm() and divide() instead (cf. sqrfree() and frac_cancel()). @Richy: Are we still using Maple naming conventions? Should this be called gcdex()? Bye, Christian -- / Physics is an algorithm \/ http://www.uni-mainz.de/~bauec002/
On Tue, 16 Nov 2004, Christian Bauer wrote:
@Richy: Are we still using Maple naming conventions? Should this be called gcdex()?
I don't see a reason why Maple naming conventions should be any better than others. If we've been following Maple naming conventions it's been because we didn't know better. Maple's argument ordering doesn't follow any conventions either and fortunately we chose not to adopt it. Well, let's see how others call it: Maple: gcdex gcdex(A,B,x,'s','t'). Note that the signature is different, anyways. CLN: xgcd There is cl_I xgcd(const cl_I& a, const cl_I& b, cl_I* u, cl_I* v). Integers only! Pari: bezout After the B'ezout Identity. The signature is: bezout(x,y): gives a 3-dimensional row vector [u,v,d] such that d=gcd(x,y) and u*x+v*y=d. Mma: ExtendedGCD That long name follows Mathematica's unique verbose naming convention. MuPAD: gcdex They copied all the conventions from Maple when they started with the whole system. Singular: extgcd extgcd(x,y) returns a list of 3 objects of the same type as the type of the arguments. Lidia: xgcd But it's not the same as in CLN: they reverse (a,b) <-> (u,v): bigint xgcd(bigint& u, bigint& v, const bigint& a, const bigint& b) [shudder] Integers only! ... That should be enough to give up looking for a consensus. @Ralf: Can you say anything about the uniqueness of the cofactors computed by your xgcd function? I mean, when you have g==u*a+v*b you also have g==(u+b)*v+(v-a)*b and so on. With integers it is trivial to normalize u and v, e.g. to the smallest possible absolute values. Regards -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Maple: gcdex CLN: xgcd Pari: bezout Mma: ExtendedGCD MuPAD: gcdex Singular: extgcd Lidia: xgcd
@Ralf: Can you say anything about the uniqueness of the cofactors computed by your xgcd function? I mean, when you have g==u*a+v*b you also have g==(u+b)*v+(v-a)*b and so on. With integers it is trivial to normalize u and v, e.g. to the smallest possible absolute values.
You mean g==(u+b)*a+(v-a)*b, no I can't say anything about this, the literature is mostly about integers, and I have not the background to tackle it. BTW the Magma manual says Ch. 44 UNIVARIATE POLYNOMIAL RINGS 187 Xgcd(f, g) XGCD(f, g) The extended greatest common divisor of polynomials f and g in a univariate polynomial ring P : the function returns polynomials c, a and b in P such that c is the GCD f and g (as defined in the function GreatestCommonDivisor), and a * f + b * g = c. The coefficient ring must be a field. Since the GCD c is unique, the multipliers a and b are unique if f and g are both non-zero. For polynomials over the rational field, a modular algorithm due to Allan Steel(unpublished) is used; over other fields the basic Euclidean algorithm is used. So there is another algorithm to ponder... ralf
Hi! On Wed, 17 Nov 2004, Ralf Stephan wrote:
@Ralf: Can you say anything about the uniqueness of the cofactors computed by your xgcd function? I mean, when you have g==u*a+v*b you also have g==(u+b)*v+(v-a)*b and so on. With integers it is trivial to normalize u and v, e.g. to the smallest possible absolute values.
You mean g==(u+b)*a+(v-a)*b, no I can't say anything about this,
If that turns out to be undetermined, it should just be documented as such.
the literature is mostly about integers, and I have not the background to tackle it. BTW the Magma manual says
Ch. 44 UNIVARIATE POLYNOMIAL RINGS 187
Xgcd(f, g) XGCD(f, g)
The extended greatest common divisor of polynomials f and g in a univariate polynomial ring
P : the function returns polynomials c, a and b in P such that c is the GCD f and g (as defined in the function GreatestCommonDivisor), and a * f + b * g = c. The coefficient ring must be a field. Since the GCD c is unique, the multipliers a and b are unique if f and g are both non-zero.
As it stands, that last sentence is not correct. To see that, just drop any two elements of any integral domain in a and b and try shifting them according to above equation. Maybe there is some other way to "normalize" u and v to the most simple form? Or maybe that is a non-issue for some reason we don't see?
For polynomials over the rational field, a modular algorithm due to Allan Steel(unpublished) is used; over other fields the basic Euclidean algorithm is used.
So there is another algorithm to ponder...
...not unless Allan shares some of his insights with us. :-) Regards -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Richard B. Kreckel wrote:
You mean g==(u+b)*a+(v-a)*b, no I can't say anything about this,
If that turns out to be undetermined, it should just be documented as such.
If a and b are coprime, then there is a unique solution if g of degree less than deg(a)+deg(b), the normalization is done by euclidean division of u and v by b and a
For polynomials over the rational field, a modular algorithm due to Allan Steel(unpublished) is used; over other fields the basic Euclidean algorithm is used.
So there is another algorithm to ponder...
...not unless Allan shares some of his insights with us. :-)
It's probably done by solving the same Bezout identity for a few primes and reconstruction by Chinese remaindering, just like for (in my opinion fastest, except for small inputs) modular GCD algorithm. The main difference being that the upper bound for coefficients in u and v is larger (bounded by the resultant of a and v instead of by the Landau-Mignotte bound). If my recollection is correct, von zur Gathen & Gerhard book describes a modular Bezout algorithm. Bernard Parisse
participants (4)
-
bernard.parisse
-
Christian Bauer
-
Ralf Stephan
-
Richard B. Kreckel