[GiNaC-list] [PATCH] gcd: avoid expanding large expressions

[Sheplyakov Alexei]

Hello!

gcd(a, b, &ca, &cb) seems to [almost] always return gcd (and 
co-factors) in expanded form. Program below demonstrates this:

/**
 * gcd_demo.cpp Demonstrate bug/feature of GiNaC's gcd
 */
#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

int main(int argc, char** argv)
{
	symbol x("x");
	symbol y("y");
	ex ca, cb;

	ex a = pow(pow(y,2)-pow(x,2), 10);
	ex b = pow(pow(x,4)-pow(y,4), 12);
	ex g = gcd(a, b, &ca, &cb, false);

	cout << "----------------------------------------------------------" << endl;
	cout << "a = " << a << endl;
	cout << "b = " << b << endl;
	cout << "GCD(a, b) = " << g << endl;
	cout << "a co-factor = " << ca << endl;
	cout << "b co-factor = " << cb << endl;
	cout << "---------------------------------------------------------" << endl;

	return 0;
}

The output of this program is

----------------------------------------------------------
a = (y^2-x^2)^10
b = (-y^4+x^4)^12
GCD(a, b) = -120*y^14*x^6+45*y^4*x^16-10*y^2*x^18+45*y^16*x^4-252*y^10*x^10+x^20+y^20+210*y^8*x^12-10*y^18*x^2-120*y^6*x^14+210*y^12*x^8
a co-factor = 1
b co-factor = 22*y^10*x^18+y^28+10*y^26*x^2+x^28+22*y^18*x^10+10*y^2*x^26-264*y^14*x^14+121*y^8*x^20-165*y^16*x^12+43*y^4*x^24+121*y^20*x^8+100*y^6*x^22-165*y^12*x^16+100*y^22*x^6+43*y^24*x^4
---------------------------------------------------------

Sometimes this [mis]feature causes .normal() to do a lot of unnecessary
job (in particular, if the expression consists of terms with denominators
like \prod(m_i^2-m_j^2)^n), so even with small expressions (about 500
terms) .normal() takes *really* long.

Attached patch fixes this behaviour, thus the program above will print

----------------------------------------------------------
a = (y^2-x^2)^10
b = (-y^4+x^4)^12
GCD(a, b) = (-y^2+x^2)^10
a co-factor = 1
b co-factor = (-y^2+x^2)^2*(y^2+x^2)^12
---------------------------------------------------------


-- 
Best regards,
  Alexei

-------------- next part --------------
Index: normal.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/normal.cpp,v
retrieving revision 1.94.2.5
diff -r1.94.2.5 normal.cpp
1346a1347
> 		const ex& exp_a = a.op(1);
1348c1349,1351
< 			if (p.is_equal(b.op(0))) {
---
> 			ex pb = b.op(0);
> 			const ex& exp_b = b.op(1);
> 			if (p.is_equal(pb)) {
1350d1352
< 				ex exp_a = a.op(1), exp_b = b.op(1);
1364c1366,1391
< 			}
---
> 			} else {
> 				ex p_co, pb_co;
> 				ex p_gcd = gcd(p, pb, &p_co, &pb_co, check_args);
> 				if (p_gcd.is_equal(_ex1)) {
> 					// a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==>
> 					// gcd(a,b) = 1
> 					if (ca)
> 						*ca = a;
> 					if (cb)
> 						*cb = b;
> 					return _ex1;
> 					// XXX: do I need to check for p_gcd = -1?
> 				} else {
> 					// there are common factors:
> 					// a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==>
> 					// gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m
> 					if (exp_a < exp_b) {
> 						return power(p_gcd, exp_a)*
> 							gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false);
> 					} else {
> 						return power(p_gcd, exp_b)*
> 							gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false);
> 					}
> 				} // p_gcd.is_equal(_ex1)
> 			} // p.is_equal(pb)
> 
1372a1400,1414
> 			} 
> 
> 			ex p_co, bpart_co;
> 			ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
> 
> 			if (p_gcd.is_equal(_ex1)) {
> 				// a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
> 				if (ca)
> 					*ca = a;
> 				if (cb)
> 					*cb = b;
> 				return _ex1;
> 			} else {
> 				// a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
> 				return p_gcd*gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false);
1374c1416,1417
< 		}
---
> 		} // is_exactly_a<power>(b)
> 
1384a1428,1444
> 
> 		ex p_co, apart_co;
> 		const ex& exp_b(b.op(1));
> 		ex p_gcd = gcd(a, p, &apart_co, &p_co, false);
> 		if (p_gcd.is_equal(_ex1)) {
> 			// b=p(x)^n, gcd(a, p) = 1 ==> gcd(a, b) == 1
> 			if (ca)
> 				*ca = a;
> 			if (cb)
> 				*cb = b;
> 			return _ex1;
> 		} else {
> 			// there are common factors:
> 			// a(x) = g(x) A(x), b(x) = g(x)^n B(x)^n ==> gcd = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
> 
> 			return p_gcd*gcd(apart_co, power(p_gcd, exp_b-1)*power(p_co, exp_b), ca, cb, false);
> 		} // p_gcd.is_equal(_ex1)
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: Digital signature
Url : http://www.cebix.net/pipermail/ginac-list/attachments/20050506/71cdb9de/attachment.pgp