Page MenuHomeHEPForge
Files F11951693

Polynomial.cc
No OneTemporary
Actions

Polynomial.cc
View Options

// -*- C++ -*-
/* Polynomial.cc
* Contains definition of the class Polynomial and associated types and operators
* Created on: Jul 7, 2010
* Author: malin
*/
#include "Polynomial.h"
#include <limits>
using namespace std;
namespace ColorFull {
// Returning Monomial at place i.
const Monomial& Polynomial::at(int i) const {
return poly.at(i);
}
// Returning Monomial at place i.
Monomial& Polynomial::at(int i) {
return poly.at(i);
}
// Returning number of terms in Polynomial
int Polynomial::size() const {
return poly.size();
}
// Adding a Monomial
void Polynomial::push_back(Monomial Mon) {
poly.push_back(Mon);
}
// Erase the Monomial at place i
void Polynomial::erase(int i) {
poly.erase(poly.begin() + i);
}
// Is the polynomial empty?
bool Polynomial::empty() const {
if (poly.empty())
return true;
else
return false;
}
// Erase info in polynomial
void Polynomial::clear() {
poly.clear();
}
// Take complex conjugate of the polynomial
void Polynomial::conjugate(){
// Loop over Monomials and complex conjugate them
for(uint i=0; i< poly.size(); i++) poly.at(i).conjugate();
}
// Collects terms with same power of 2 and N and Cf
void Polynomial::simplify() {
// If empty Polynomial or only one Monomial, do nothing
if (poly.size() <= 1)
return;
// Collect different Monomials, keep 0th Monomial
polynomial terms;
terms.push_back(poly.at(0));
poly.erase(poly.begin());
// Move terms from Polynomial to unique set of terms
// as long as terms remain in Polynomial
while (poly.size() > 0) {
bool was_found = false;
// Check if term already there
for (uint i = 0; (i < terms.size() && !was_found); i++) {
// If same powers of 2, N and cf
if (poly.at(0).pow_2 == terms.at(i).pow_2 && poly.at(0).pow_N
== terms.at(i).pow_N && poly.at(0).pow_cf
== terms.at(i).pow_cf) {
was_found = true;
// Change the n_times (if no numerical factor)
// or the num factor of term to be added
cnum for_testing = 1.0;
// If already stored term, or part to be added, has a numeric part,
// store info in numeric part
if ( poly.at(0).num != for_testing or terms.at(i).num != for_testing) {
// Float to add to num, put numerical factor, and n_times,
// and sign here
// cnum c_to_add=to_add*poly.at(0).num;
cnum c_to_add;
c_to_add = poly.at(0).num;
c_to_add *= poly.at(0).n_times;
// Move int-part of existing Col_str to float part
terms.at(i).num *= terms.at(i).n_times;
terms.at(i).n_times = 1;
// Add cnum from the term
terms.at(i).num += c_to_add;
} else { //if no numeric part change the int-part instead
terms.at(i).n_times += poly.at(0).n_times;
}
}
}
// If term not already there, add term
if (!was_found)
terms.push_back(poly.at(0));
// Erase the term in poly
poly.erase(poly.begin());
}
poly = terms;
}
// Define the operator << for Polynomial
std::ostream& operator<<(std::ostream& out, const Polynomial & Poly) {
if (Poly.size() == 0)
out << "1";
else if (Poly.size() == 1) out << Poly.at(0);//changed 10 06 11
// out << "size 1";
else {
if (Poly.size() > 1)
out << "(";
// Write out 0th element only if !=0, or only element
if ( Poly.size() == 1 )
out << Poly.at(0);
else if(!Poly.at(0).n_times == 0 ){
out << Poly.at(0);
out << "+";
}
for (int i = 1; i < Poly.size(); i++) {
out << Poly.at(i);
if (i != Poly.size() - 1)
out << "+";
}
if (Poly.size() > 1)
out << ")";
}
return out;
}
// Define the operator == for Polynomial
// By def, each TERM has to be identical, this is against intuition
// as order matters 1+2 is not 2+1 (TODO, change this (define normal order?) )
bool operator==( const Polynomial & Poly1, const Polynomial & Poly2){
if(Poly1.size() != Poly2.size() ) return false;
// All terms should be the same in order for Polynomials to be the same
for (int i=0; i < Poly1.size(); i++ ){
if( Poly1.at(i) != Poly2.at(i ) ) return false;
}
return true;
}
// Define the operator != for Polynomial
bool operator!=( const Polynomial & Poly1, const Polynomial & Poly2){
if( Poly1==Poly2) return false;
else return true;
}
// Define the operator + for Polynomial and a single Monomial
Polynomial operator+(const Polynomial & Poly, const Monomial & Mon){
//cout << "operator+(Poly, Mon): Incoming Polynomial " << Poly << endl;
//cout << "operator+(Poly, Mon): Incoming Monomial " << Mon << endl;
Polynomial out_Poly(Poly);
// If original Polynomial was empty=1, push_back dummy Monomial=1
// (such that after the addition the Polynomial has indeed two terms)
if(out_Poly.empty()){
//cout << "operator+(Poly, Mon): Empty Polynomial " << Poly << endl;
Monomial dummy_Mon;
out_Poly.push_back(dummy_Mon);
//cout << "operator+(Poly, Mon): Appended dummy Mon " << Poly << endl;
}
// Add term unless it's 0
if(Mon.n_times != 0) out_Poly.push_back( Mon );
//cout << "operator+(Poly, Mon): After addition " << Poly << endl;
return out_Poly;
}
// Define the operator - for Polynomial and a single Monomial
Polynomial operator-( const Polynomial & Poly, const Monomial & Mon){
Monomial out_Mon(Mon);
// Change sing of monomial
out_Mon.n_times*=(-1); // Corrected sign in unused operator 12 09 18, was (1)
// Add term unless it's 0
return Poly+out_Mon;
}
// Define the operator - for Polynomial and a single Monomial
Polynomial operator-(const Monomial & Mon, const Polynomial & Poly){
Polynomial out_Poly(Poly);
// If Poly is empty=1, add a default Monomial=1, which can change sign
if (out_Poly.empty()) {
Monomial dummy_Mon;
out_Poly.push_back(dummy_Mon);
}
// Change sing of Polynomial, by looping over terms
for (uint m=0; m < out_Poly.poly.size(); m++)
out_Poly.poly.at(m).n_times*=(-1);
// Add Polynomial and Monomial
return out_Poly+Mon;
}
// Define the operator + Polynomial and single Monomial
Polynomial operator+(const Monomial & Mon, const Polynomial & Poly ){
return Poly + Mon;
}
// Define the operator + for Polynomials
// Recently changed 101118
Polynomial operator+( const Polynomial & Poly1, const Polynomial & Poly2) {
Polynomial out_Poly1(Poly1);
Polynomial out_Poly2(Poly2);
// If original Poly2 was empty=1, push_back dummy Monomial=1
// after adding empty Monomial the Polynomial is still 1
if (out_Poly2.empty()) {
Monomial dummy_Mon;
out_Poly2.push_back(dummy_Mon);
}
// Add terms from Poly2 to Poly1
for (int i = 0; i < out_Poly2.size(); i++) {
out_Poly1 = out_Poly1 + out_Poly2.at(i);
}
return out_Poly1;
}
// Define the operator - for Polynomials
Polynomial operator-(const Polynomial & Poly1, const Polynomial & Poly2) {
Polynomial out_Poly2(Poly2);
// If Poly2 is empty=1 make it contain a Monomial,
// as a Monomial by construction is 1
if(out_Poly2.empty()) {
Monomial dummy_Mon;
out_Poly2.push_back(dummy_Mon);
}
// Change sing of Poly2, by looping over terms
for (uint m=0; m < out_Poly2.poly.size(); m++)
out_Poly2.poly.at(m).n_times*=(-1);
// Now the subtraction is equal to an addition
return Poly1+out_Poly2;
}
// Define the operator * for Polynomial and int
Polynomial operator*(const Polynomial & Poly, const int in){
Polynomial out_Poly(Poly);
// If initially empty Polynomial=1, append empty Monomial=1, to have something to multiply with
Monomial dummy_Mon;
if (out_Poly.empty())
out_Poly.push_back(dummy_Mon);
// Loop over terms in Polynomial
for (int i=0; i< out_Poly.size(); i++){
out_Poly.at(i).n_times=out_Poly.at(i).n_times*in;
}
return out_Poly;
}
// Define the operator * for int and Polynomial
Polynomial operator*(const int in, const Polynomial & Poly) {
return Poly * in;
}
// Define the operator * for Polynomial and cnum
Polynomial operator*(const Polynomial & Poly, const cnum c) {
Polynomial out_Poly(Poly);
// If initially empty Polynomial=1, append empty Monomial=1, to have something to multiply with
Monomial dummy_Mon;
if (out_Poly.empty())
out_Poly.push_back(dummy_Mon);
// Loop over terms in Polynomial
for (int i = 0; i < out_Poly.size(); i++) {
//cout << "in loop with i=" << i;
out_Poly.at(i) = out_Poly.at(i) * c;
}
return out_Poly;
}
Polynomial operator*(const cnum c, const Polynomial & Poly) {
return Poly * c;
}
// Define the operator * for Polynomial and double
Polynomial operator*(const Polynomial & Poly, const double d) {
Polynomial out_Poly(Poly);
// If initially empty Polynomial=1, append empty Monomial=1, to have something to multiply with
Monomial dummy_Mon;
if (out_Poly.empty())
out_Poly.push_back(dummy_Mon);
// Loop over terms in Polynomial
for (int i = 0; i < out_Poly.size(); i++) {
out_Poly.poly.at(i).num = out_Poly.at(i).num * d;
}
return out_Poly;
}
// Define the operator * for Polynomial and double
Polynomial operator*(const double d, const Polynomial & Poly) {
return Poly * d;
}
// Define the operator * for Polynomial and Monomial
Polynomial operator*(const Polynomial & Poly, const Monomial & Mon){
// To contain the result
Polynomial Poly_res;
// Special case that the Polynomial is empty=1
if( Poly.empty() ) {
Poly_res.push_back(Mon);
return Poly_res;
}
else{
for (int i1=0; i1< Poly.size(); i1++){
Poly_res.push_back( Poly.at(i1)*Mon);
}
}
return Poly_res;
}
// Define the operator * for Monomial and Polynomial
Polynomial operator*(const Monomial & Mon, const Polynomial & Poly){
return Poly*Mon;
}
// Define the operator * for Polynomials
Polynomial operator*( const Polynomial & Poly1, const Polynomial & Poly2){
// To contain the result
Polynomial Poly_res;
// Special case that the vectors are empty=1, needed for special treatment of sum of 0-monomials
// Poly_res is also empty, and thus =1, 1*1=1
if(Poly1.empty() && Poly2.empty()){
//Monomial Mon_tmp;
//Poly_res.push_back( Mon_tmp);
return Poly_res;
}
// If one Poly is empty, =1, return other
if(!Poly1.empty() && Poly2.empty() ) return Poly1;
else if( Poly1.empty() && !Poly2.empty()) return Poly2;
// If both vectors contain info
else {
for (int i1=0; i1< Poly1.size(); i1++){
for (int i2=0; i2< Poly2.size(); i2++){
// Add terms if non of them is 0, else just don't add
if( Poly1.at(i1).n_times!=0 && Poly2.at(i2).n_times!=0 )
Poly_res.push_back(Poly1.at(i1)*Poly2.at(i2));
//Poly_res= Poly_res+ ( Poly1.at(i1)*Poly2.at(i2) ); // if used initialize poly to 0 first
//cout << "operator*(Poly1, Poly2): for " << i1 <<" " << i2 << " " << Poly_res <<endl;
}
//cout << "operator*(Poly1, Poly2): Looped over " << i1 <<endl;
}
// If, by now, the Poly_res is empty that's because all terms where 0
if( Poly_res.empty() ){
Monomial Mon_tmp;
Mon_tmp.n_times=0;
Poly_res.push_back( Mon_tmp );
}
}
return Poly_res;
}
// Define the operator << for vector of Polynomial
std::ostream& operator<<(std::ostream& out, const Poly_vec & Poly_v){
out <<"{";
// Loop over entries
for( uint i=0; i< Poly_v.size(); i++ ){
out << Poly_v.at(i);
// If not last element print ","
if (i<Poly_v.size()-1 ) out << ", ";
}
out <<"}";
return out;
}
// Define the operator << for Poly_matr
std::ostream& operator<<(std::ostream& out, const Poly_matr & Poly_m){
out <<"{" <<endl;
// Loop over rows
for(uint i=0; i< Poly_m.size(); i++ ){
out <<"{";
// Loop over columns
for(uint j=0; j< Poly_m.at(i).size(); j++ ){
// Print element
cout.width( 60 );
ostringstream outstr;
outstr << Poly_m.at(i).at(j);
// If not last element print ","
if (j<Poly_m.at(i).size()-1 ) outstr << ",";
out << outstr.str();
//out << Poly_m.at(i).at(j) << "";
}
out <<"}";
// If not last row, print ","
if (i<Poly_m.at(i).size()-1 ) out << ",";
out << endl;
}
out <<"}" <<endl;
return out;
}
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Jun 3, 4:34 PM (9 h, 19 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5323897
Default Alt Text
Polynomial.cc (11 KB)

Event Timeline