[GiNaC-list] code: extended gcd

[Ralf Stephan]

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");
}