[GiNaC-devel] Problems with dummy summation

[Vladimir Kisil]

		Dear All,

		I recently noticed that GiNaC incorrectly expands products with
  dummy index summations, which is illustrated by the following program:

#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;
#include <stdexcept>

int main(){
	realsymbol a("a"), b("b"), c("c"), d("d");
	varidx mu(symbol("mu", "\\mu"), 2);
	try {
		ex e;
		e = (indexed(a, mu)*indexed(b, mu.toggle_variance())+1)*indexed(c, mu)*indexed(d, mu.toggle_variance());
		cout << e.expand() << endl;
		// -> d.mu*c~mu*b.mu*a~mu+d.mu*c~mu

		e = pow(indexed(a, mu)*indexed(b, mu.toggle_variance())+1, 2);
		cout << e.expand() << endl;
		// -> 1+2*b.mu*a~mu+(b.mu)^2*(a~mu)^2
	} catch  (exception &p) {
			cerr << "Got problem: " <<  p.what() << endl;
	}
}

  I made a patch to mul.cpp, ncmul.cpp, power.cpp which _partially_ fix
  this problem. The bigger test looks now like that:

#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;
#include <stdexcept>

int main(){
	realsymbol a("a"), b("b"), c("c"), d("d");
	varidx mu(symbol("mu", "\\mu"), 2),  nu(symbol("nu", "\\nu"), 2);
	try {
		ex e;
		e = (indexed(a, mu)*indexed(b, mu.toggle_variance())+1)*indexed(c, mu)*indexed(d, mu.toggle_variance());
		cout << e.expand() << endl << endl;
		// -> d.mu*c~mu*b.symbol6*a~symbol6+d.mu*c~mu

		e = indexed(c, mu)*indexed(d, mu.toggle_variance())*(indexed(a, mu)*indexed(b, mu.toggle_variance())+1);
		cout << e.expand() << endl << endl;
		// -> d.mu*c~mu+d.mu*a~symbol7*b.symbol7*c~mu

		e = (1+indexed(c, mu)*indexed(d, mu.toggle_variance()))*(indexed(a, mu)*indexed(b, mu.toggle_variance())+1);
		cout << e.expand() << endl << endl;
		// -> 1+b.mu*a~mu+b.mu*c~symbol8*a~mu*d.symbol8+d.mu*c~mu

		e = pow(2*indexed(a, mu)*indexed(b, mu.toggle_variance())+1, 2);
		cout << e.expand() << endl << endl;
		// -> 1+4*b.mu*a~mu*b.symbol9*a~symbol9+4*b.mu*a~mu

		e = pow(2*indexed(a, mu)*indexed(b, mu.toggle_variance())+1, 3);
		cout << e.expand() << endl << endl;
		// -> 1+6*b.mu*a~mu+12*b.mu*a~mu*a~symbol10*b.symbol10+8*b.symbol12*a~symbol12*b.mu*a~symbol11*b.symbol11*a~mu

		e = (indexed(a, mu, nu)*dirac_gamma(mu.toggle_variance()) * dirac_gamma(nu.toggle_variance()) + dirac_ONE()) 
			* indexed(b, mu, nu)*dirac_gamma(mu.toggle_variance()) * dirac_gamma(nu.toggle_variance());
		cout << e.expand() << endl << endl;
		// -> (gamma.mu*gamma.nu*gamma.mu*gamma.nu)*a~mu~nu*b~mu~nu+(gamma.mu*gamma.nu)*b~mu~nu

		e = (indexed(a, mu, nu)*dirac_gamma(mu.toggle_variance()) * dirac_gamma(nu.toggle_variance()) + dirac_ONE()) 
			*( indexed(b, mu, nu)*dirac_gamma(mu.toggle_variance()) * dirac_gamma(nu.toggle_variance()) + dirac_ONE());
		cout << e.expand() << endl << endl;
		// -> b~symbol16~symbol17*a~mu~nu*(gamma.mu*gamma.nu*gamma.symbol16*gamma.symbol17)+ONE+(gamma.mu*gamma.nu)*a~mu~nu+(gamma.mu*gamma.nu)*b~mu~nu

		e = (indexed(a, mu, nu)*dirac_gamma(mu.toggle_variance()) * dirac_gamma(nu.toggle_variance()) + dirac_ONE()) 
			* indexed(b, mu.toggle_variance(), nu.toggle_variance())*dirac_gamma(mu) * dirac_gamma(nu);
		cout << e.expand() << endl << endl;
		// -> (gamma.mu*gamma.nu*gamma~mu*gamma~nu)*a~mu~nu*b.mu.nu+(gamma~mu*gamma~nu)*b.mu.nu

		e = pow(indexed(a, mu, nu)*dirac_gamma(mu.toggle_variance()) * dirac_gamma(nu.toggle_variance()), 2);
		cout << e.expand() << endl << endl;
		// -> (gamma.mu*gamma.nu*gamma.symbol18*gamma.symbol19)*a~mu~nu*a~symbol18~symbol19

		e = pow(2*indexed(a, mu, nu)*dirac_gamma(mu.toggle_variance()) * dirac_gamma(nu.toggle_variance())+dirac_ONE(), 2);
		cout << e.expand() << endl << endl;
		// -> 4*(gamma.mu*gamma.nu*gamma.symbol18*gamma.symbol19)*a~mu~nu*a~symbol18~symbol19+ONE+4*(gamma.mu*gamma.nu)*a~mu~nu

		e = pow(2*indexed(a, mu, nu)*dirac_gamma(mu.toggle_variance()) * dirac_gamma(nu.toggle_variance())+dirac_ONE(), 3);
		cout << e.expand() << endl << endl;
		// -> ONE+8*a~symbol26~symbol27*a~mu~nu*a~symbol28~symbol29*(gamma.mu*gamma.nu*gamma.symbol26*gamma.symbol27*gamma.symbol28*gamma.symbol29)+6*(gamma.mu*gamma.nu)*a~mu~nu+4*(gamma.mu*gamma.nu*gamma.symbol24*gamma.symbol25)*a~mu~nu*a~symbol24~symbol25+4*a~symbol22~symbol23*a~mu~nu*(gamma.mu*gamma.nu*gamma.symbol22*gamma.symbol23)+4*a~symbol20~symbol21*a~mu~nu*(gamma.mu*gamma.nu*gamma.symbol20*gamma.symbol21)
	} catch  (exception &p) {
			cerr << "Got problem: " <<  p.what() << endl;
	}
}
  
  So the problem remains in the cases when one dummy index belongs to
  non-commutative part of product and other---to a commutative part. I do
  not know how to fix this situation in a simple way.

  I sent my patch as an attachment to a separate letter to this list
  (which will be likely hold till moderator's  approval). It includes as
  well a slightly redesigned version of my previous patch for indexed.cpp
  and clifford.cpp (from the Thu Apr 28 18:17:55 CEST 2005).

  Best regards,
  Vladimir

-- 
Vladimir V. Kisil     email: kisilv at maths.leeds.ac.uk
--                      www: http://maths.leeds.ac.uk/~kisilv/