Optimise AdjformEx to use int32_t instead of mpq_class

Author: Dominic Price

core/Adjform.cc

@@ -359,20 +359,22 @@ namespace cadabra {
 		return -(Adjform::value_type)data->begin().number_of_children();
 	}
 
-	AdjformEx::rational_type AdjformEx::zero = 0;
+	AdjformEx::integer_type AdjformEx::zero = 0;
 
 	AdjformEx::AdjformEx()
 	{
 
 	}
 
-	AdjformEx::AdjformEx(const Adjform& adjform, const AdjformEx::rational_type& value, const Ex& prefactor)
+	AdjformEx::AdjformEx(const Adjform& adjform, const AdjformEx::integer_type& value, const Ex& prefactor)
 		: prefactor(prefactor)
 	{
+		sizeof(mpq_class);
+		sizeof(int64_t);
 		set(adjform, value);
 	}
 
-	AdjformEx::AdjformEx(const Adjform& adjform, const rational_type& value, Ex::iterator prefactor)
+	AdjformEx::AdjformEx(const Adjform& adjform, const integer_type& value, Ex::iterator prefactor)
 		: prefactor(prefactor)
 	{
 		set(adjform, value);
@@ -412,7 +414,7 @@ namespace cadabra {
 		one(tensor_begin->multiplier);
 	}
 
-	AdjformEx::rational_type AdjformEx::compare(const AdjformEx& other) const
+	AdjformEx::integer_type AdjformEx::compare(const AdjformEx& other) const
 	{
 		// Early failure checks
 		if (data.empty() || data.size() != other.data.size())
@@ -421,7 +423,7 @@ namespace cadabra {
 		// Find the numeric factor between the first two terms, then loop over all
 		// other terms checking that the factor is the same. If not, return 0
 		auto a_it = data.begin(), b_it = other.data.begin(), a_end = data.end();
-		rational_type factor = a_it->second / b_it->second;
+		integer_type factor = a_it->second / b_it->second;
 		while (a_it != a_end) {
 			if (a_it->first != b_it->first)
 				return 0;
@@ -438,13 +440,13 @@ namespace cadabra {
 			add(kv.first, kv.second);
 	}
 
-	void AdjformEx::combine(const AdjformEx& other, AdjformEx::rational_type factor)
+	void AdjformEx::combine(const AdjformEx& other, AdjformEx::integer_type factor)
 	{
 		for (const auto& kv : other.data)
 			add(kv.first, kv.second * factor);
 	}
 
-	void AdjformEx::multiply(const AdjformEx::rational_type& k)
+	void AdjformEx::multiply(const AdjformEx::integer_type& k)
 	{
 		for (auto& kv : data)
 			kv.second *= k;
@@ -519,33 +521,33 @@ namespace cadabra {
 		return tensor;
 	}
 
-	const AdjformEx::rational_type& AdjformEx::get(const Adjform& adjform) const
+	const AdjformEx::integer_type& AdjformEx::get(const Adjform& adjform) const
 	{
 		auto pos = data.find(adjform);
 		return (pos == data.end()) ? zero : pos->second;
 	}
 
-	void AdjformEx::set(const Adjform& term, const AdjformEx::rational_type& value)
+	void AdjformEx::set(const Adjform& term, const AdjformEx::integer_type& value)
 	{
 		if (!term.empty())
 			set_(term, value);
 	}
 
-	void AdjformEx::set_(const Adjform& term, const AdjformEx::rational_type& value)
+	void AdjformEx::set_(const Adjform& term, const AdjformEx::integer_type& value)
 	{
 		if (value != 0)
 			data[term] = value;
 		else
 			data.erase(term);
 	}
 
-	void AdjformEx::add(const Adjform& term, const AdjformEx::rational_type& value)
+	void AdjformEx::add(const Adjform& term, const AdjformEx::integer_type& value)
 	{
 		if (!term.empty())
 			add_(term, value);
 	}
 
-	void AdjformEx::add_(const Adjform& term, const AdjformEx::rational_type& value)
+	void AdjformEx::add_(const Adjform& term, const AdjformEx::integer_type& value)
 	{
 		auto elem = data.find(term);
 		if (elem == data.end() && value != 0) {

core/Adjform.hh

@@ -84,23 +84,23 @@ namespace cadabra {
 	class AdjformEx
 	{
 	public:
-		using rational_type = mpq_class;
-		using map_t = std::map<Adjform, rational_type>;
+		using integer_type = int32_t;
+		using map_t = std::map<Adjform, integer_type>;
 		using iterator = map_t::iterator;
 		using const_iterator = map_t::const_iterator;
 
 		AdjformEx();
-		AdjformEx(const Adjform& adjform, const rational_type& value = 1, const Ex& prefactor = Ex());
-		AdjformEx(const Adjform& adjform, const rational_type& value, Ex::iterator prefactor);
+		AdjformEx(const Adjform& adjform, const integer_type& value = 1, const Ex& prefactor = Ex());
+		AdjformEx(const Adjform& adjform, const integer_type& value, Ex::iterator prefactor);
 		AdjformEx(Ex& tr, Ex::iterator it, IndexMap& index_map, const Kernel& kernel);
 
 		// Check if 'other' is a linear multiple of 'this' and return
 		// the numeric factor if so, otherwise returns 0
-		rational_type compare(const AdjformEx& other) const;
+		integer_type compare(const AdjformEx& other) const;
 
 		void combine(const AdjformEx& other); // Add all contributions from 'other' into 'this'
-		void combine(const AdjformEx& other, rational_type factor);
-		void multiply(const rational_type& k); // Multiply all terms by a constant factor
+		void combine(const AdjformEx& other, integer_type factor);
+		void multiply(const integer_type& k); // Multiply all terms by a constant factor
 
 		iterator begin();
 		const_iterator begin() const;
@@ -119,11 +119,11 @@ namespace cadabra {
 		Ex& get_tensor_ex();
 
 		// Get the value of the term, or zero if it doesn't exist
-		const rational_type& get(const Adjform& adjform) const;
+		const integer_type& get(const Adjform& adjform) const;
 		// Sets the given term to value, creating/removing the term if required
-		void set(const Adjform& adjform, const rational_type& value = 1);
+		void set(const Adjform& adjform, const integer_type& value = 1);
 		// Adds value to the given term, creating/removing the term if required
-		void add(const Adjform& adjform, const rational_type& value = 1);
+		void add(const Adjform& adjform, const integer_type& value = 1);
 
 		// Symmetrize in the given indices
 		// e.g. if the only term is abcd then
@@ -143,12 +143,12 @@ namespace cadabra {
 
 	private:
 		// Unsafe (but faster) versions of the public functions
-		void set_(const Adjform& adjform, const rational_type& value = 1);
-		void add_(const Adjform& adjform, const rational_type& value = 1);
+		void set_(const Adjform& adjform, const integer_type& value = 1);
+		void add_(const Adjform& adjform, const integer_type& value = 1);
 		map_t data;
 		Ex prefactor;
 		Ex tensor;
-		static rational_type zero;
+		static integer_type zero;
 	};
 }
 

core/algorithms/meld.cc

@@ -342,8 +342,8 @@ bool similar_tensor_form(const Kernel& kernel, Ex& tr, Ex::iterator it1, Ex::ite
 meld::result_t meld::apply_tableaux(iterator it)
 {
 	using namespace boost::numeric::ublas;
-	using matrix_type = matrix<AdjformEx::rational_type>;
-	using vector_type = vector<AdjformEx::rational_type>;
+	using matrix_type = matrix<mpq_class>;
+	using vector_type = vector<mpq_class>;
 
 	result_t res = result_t::l_no_action;
 
@@ -368,7 +368,7 @@ meld::result_t meld::apply_tableaux(iterator it)
 		// ensure that when solving for linear dependence there are as many
 		// unknowns as equations
 		matrix_type coeffs;
-		LinearSolver<AdjformEx::rational_type> solver;
+		LinearSolver<mpq_class> solver;
 
 		// The adjform in position 'i' of 'mapping' represents the term corresponding
 		// to the 'i'th row of 'coeffs'
@@ -425,7 +425,7 @@ meld::result_t meld::apply_tableaux(iterator it)
 						// Ensure that the next term is bigger than any other term to begin with
 						Adjform next_term;
 						next_term.push_back(std::numeric_limits<Adjform::value_type>::max());
-						AdjformEx::rational_type sum = 0;
+						mpq_class sum = 0;
 						for (size_t i = 0; i < adjforms.size(); ++i) {
 							if (lhs_its[i] != adjforms[i].end() && lhs_its[i]->first == cur_term) {
 								sum += x(i) * lhs_its[i]->second;
@@ -437,7 +437,7 @@ meld::result_t meld::apply_tableaux(iterator it)
 								next_term = lhs_its[i]->first;
 						}
 
-						AdjformEx::rational_type rhs_sum;
+						mpq_class rhs_sum;
 						if (rhs_it == cur_adjform.end() || rhs_it->first != cur_term) {
 							rhs_sum = 0;
 						}