Valgrind errors with std::vector

Greetings.

I've been working on a polynomial calculator. I was very happy to have finally finished my implementation of Newton's Tangent Method as well as Ruffini's Rule. Together, I intended to use them to find all the zeros of a function.

However, as I started combining the two, I realized that removing a cout statement or adding it in a function would cause the said function to produce different outputs. This made me think that I was seeing either stack or heap corruption.

So I opened up valgrind and I was astonished with the amount of errors I had.

The first error seems the most important one, and the one that leads to several others. But to show it to you, I have to place part of the code here.

I first have a class named Term:

Code:

---

class Term {
public:
 // Methods...
  unsigned int exponent;
  double quocient;
protected:
private:


};

---

The Term class doesn't seem to have any problems at all and those are its only members: a uint32 named exponent and a double named quocient.

Next is the Polynomial class:

Code:

---

class Polynomial {
public:
  Polynomial();
  bool addTerm(Term term);
  bool add0Term(unsigned int exponent);
  //Some other methods
  std::vector<Term> terms;
protected:
private:

};

---

This is where everything gets ruined. I have the member terms which is std::vector<Term>.

During the implementation of addTerm and add0Term (both with the same code, the difference lies just in adding a zero-quocient Term), I check for an empty list as follows:

Code:

---

bool Polynomial::addTerm(Term term) {
  bool anybigger = false;

  ////
  // If this is the first, just push it.
  //
  if(terms.size()==0) {
      terms.push_back(term);
      return true;
  }

//More code...

---

Now, the first valgrind error happens at exactly terms.push_back(term), here's the error:

Quote:

---

==30949== Invalid read of size 4
==30949== at 0x804990F: Polynomial::addTerm(Term) (main.cxx:165)
==30949== by 0x804AD11: main (main.cxx:517)
==30949== Address 0x4037474 is 0 bytes after a block of size 12 alloc'd
==30949== at 0x4008390: operator new(unsigned int) (vg_replace_malloc.c:214)
==30949== by 0x804C4A1: __gnu_cxx::new_allocator<Term>::allocate(unsigned int, void const*) (new_allocator.h:89)
==30949== by 0x804C211: std::_Vector_base<Term, std::allocator<Term> >::_M_allocate(unsigned int) (stl_vector.h:140)
==30949== by 0x804BBC6: std::vector<Term, std::allocator<Term> >::_M_insert_aux(__gnu_cxx::__normal_iterator<Term *, std::vector<Term, std::allocator<Term> > >, Term const&) (vector.tcc:322)
==30949== by 0x804B494: std::vector<Term, std::allocator<Term> >::push_back(Term const&) (stl_vector.h:741)
==30949== by 0x80498A2: Polynomial::addTerm(Term) (main.cxx:153)
==30949== by 0x804ACE7: main (main.cxx:516)

---

And the sequence of calls is presented here:

Code:

---

//main.cxx:516
Term term1(4,-5);
//... more code
Polynomial poly;
//... more code
poly.addTerm(term1);

---

Code:

---

//main.cxx:153
bool Polynomial::addTerm(Term term) {
  //More code...

  ////
  // If this is the first, just push it.
  //
  if(terms.size()==0) {
      terms.push_back(term);
      return true;
  }

//More code...

---

I even ported my application to use the simple vector class I made some time ago and which didn't implement iterators or standard allocators. Even with it, I reached this error, so I think it is safe to admit that *I* am doing something wrong...but what is it? What's wrong with pushing back a term? I should note that even adding a terms.reserve(1); before would make valgrind complain at that same line, so the problem isn't with me adding the term, but it seems somehow related to how terms was created/constructed?

Can somebody help me understand this, please? Otherwise I'll just dump vectors and good with good old arrays, memcpy(), memmove(), new and delete.

Maybe I forgot to implement some kind of important constructor and it defaulted to a buggy one? I doubt that, though...

Thanks,

Jorl17

There is nothing amiss in the code you've posted.
However I certainly believe that there is a problems, but it is just caused by other parts of the code.
The problem is not with vectors themselves.
More code?

Well, I'll certainly post the full code here. Mind any bugs or mis-assumptions, just check to see if I don't do anything nuts. Also, the whole code is a bit messy, since I intend to clean it up and split it across classes (as well as make members private, use strings instead of qstring, etc).

Here it is, and thank you for your help:

Code:

---

/***
 *
 *  This file is part of PolyCalc
 *  Copyright (C) 2007-2010 João Ricardo Lourenço ([email protected])
 *  ====================================================================
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ***/
#include <QString>
#include <iostream>
#include <vector>
#include <math.h>
class Term {
public:
  Term(void);
  Term(unsigned int exponent, double quocient);
  Term(const Term& o);
  void set(unsigned int exponent, double quocient);
  QString getAsString(void);
  QString getInfoAsString(void);
  QString getQuocientAsString(void);
  Term getDerivative(void);
  Term getIntegral(void);
  double getImage(double object);
  unsigned int exponent;
  double quocient;
protected:
private:


};

void Term::set(unsigned int exponent, double quocient) {
  this->exponent=exponent; this->quocient=quocient;
}
Term::Term(unsigned int exponent, double quocient) {
  this->exponent=exponent; this->quocient=quocient;
}
Term::Term(void) {
  this->exponent=0; this->quocient=0;
}
Term::Term(const Term& o) {
  this->exponent=o.exponent; this->quocient=o.quocient;
}

QString Term::getAsString(void) {
  if(quocient==0) {
      return "NaN";
  }
  QString tmp;
  if (exponent==1) {
      if(quocient==1.0f)
        tmp = "x";
      else
        tmp = QString::number(quocient, 'g', 15) + "x";
  } else if (exponent==0) {
      tmp = QString::number(quocient, 'g', 15);
  } else {
      if(quocient==1.0f)
        tmp = "x^" + QString::number(exponent);
      else
        tmp = QString::number(quocient, 'g', 15) + "x^" + QString::number(exponent);
  }
  return tmp;
}

QString Term::getQuocientAsString(void) {
  if(quocient==0) {
      return "NaN";
  }
  return QString::number(quocient, 'g', 15);

}

QString Term::getInfoAsString(void) {
  if(quocient==0) {
      return "NaN";
  }

  QString tmp = "This is a degree "+QString::number(exponent)+" term of quocient "+QString::number(quocient, 'g', 15)+".";

  return tmp;
}

Term Term::getDerivative(void) {
  Term tmp = (*this);

  if(tmp.exponent==0)
      return Term(0,0);
  tmp.quocient*=((double)(tmp.exponent--));
  return tmp;
}
Term Term::getIntegral(void) {
  Term tmp = (*this);

  if(tmp.exponent==0)
      return Term(1,quocient);

  tmp.quocient/=((double)(++tmp.exponent));
  return tmp;
}

double Term::getImage(double object) {

  if(exponent==0) {
      return quocient;
  } else if(exponent==1) {
      return quocient*object;
  } else {
      return quocient*pow(object, exponent);
  }
}

class Polynomial {
public:
  Polynomial();
  bool addTerm(Term term);
  bool add0Term(unsigned int exponent);
  QString getAsString(void);
  Polynomial getDerivative(void);
  Polynomial getIntegral(void);
  double getImage(double object);
  double findZerobyNewtons(double startX, bool& success);
  void printZeros(void);
  std::vector<Term> terms;
protected:
private:

};

Polynomial::Polynomial() {

}

bool Polynomial::addTerm(Term term) {
  bool anybigger = false;

  ////
  // If this is the first, just push it.
  //
  if(terms.size()==0) {
      terms.push_back(term);
      return true;
  }

  if(term.quocient==0)
      return false;

  ////
  // Check to see if this term has a higher degree than all others.
  // If so, then just put it in the beginning.
  //
  unsigned int first_exponent = terms.begin()->exponent;
  unsigned int last_exponent = terms.end()->exponent;
  if(term.exponent>first_exponent) {
      terms.insert(terms.begin(), term);
      return true;
  } else if(term.exponent<last_exponent) {
      terms.push_back(term);
      return true;
  }



  ////
  // See where it fits.
  //
  for (std::vector<Term>::iterator it = terms.begin(); it != terms.end(); it++) {

      ////
      // If there is already a term of such degree, bail out.
      //
      if(it->exponent==term.exponent) {
        std::cout<<"Found a match for "<<term.exponent<<std::endl;
        return false;
      }
      ////
      // If we found a term exactly before (1 degree above) this one,
      // then insert the term after it.
      //
      else if(it->exponent==(term.exponent+1)) {
        if ((it+1)->exponent==term.exponent)
            return false;

        terms.insert(it+1, term);
        return true;
      }

      ////
      // Let's toggle the "awareness" of the iteration. Now we know
      // that we have a term which is bigger in terms of degree.
      // Let's continue iterating to see if we find one that is of a mintor
      // degree.
      //
      else if (it->exponent>term.exponent)
        anybigger=true;

      ////
      // If we found a term of a lower degree and we know that there are terms
      // of bigger degree, then put it right before this one.
      //
      else if (it->exponent<term.exponent && anybigger==true) {
        it--;
        terms.insert(it, term);
        return true;
      }
  }

  ////
  // We are minor than all others, so we can insert ourselves at the end.
  // FIXME: This could be done using terms.last()->exponent>term.exponent
  //
  terms.push_back(term);
  return true;
}

bool Polynomial::add0Term(unsigned int exponent) {
  bool anybigger = false;

  Term term (exponent, 0);
  ////
  // If this is the first, just push it.
  //
  if(terms.size()==0) {
      terms.push_back(term);
      return true;
  }

  ////
  // Check to see if this term has a higher degree than all others.
  // If so, then just put it in the beginning.
  //
  unsigned int first_exponent = terms.begin()->exponent;
  unsigned int last_exponent = terms.end()->exponent;
  if(term.exponent>first_exponent) {
      terms.insert(terms.begin(), term);
      return true;
  } else if(term.exponent<last_exponent) {
      terms.push_back(term);
      return true;
  }



  ////
  // See where it fits.
  //
  for (std::vector<Term>::iterator it = terms.begin(); it != terms.end(); it++) {

      ////
      // If there is already a term of such degree, bail out.
      //
      if(it->exponent==term.exponent) {
        std::cout<<"Found a match for "<<term.exponent<<std::endl;
        return false;
      }
      ////
      // If we found a term exactly before (1 degree above) this one,
      // then insert the term after it.
      //
      else if(it->exponent==(term.exponent+1)) {
        if ((it+1)->exponent==term.exponent)
            return false;

        terms.insert(it+1, term);
        return true;
      }

      ////
      // Let's toggle the "awareness" of the iteration. Now we know
      // that we have a term which is bigger in terms of degree.
      // Let's continue iterating to see if we find one that is of a mintor
      // degree.
      //
      else if (it->exponent>term.exponent)
        anybigger=true;

      ////
      // If we found a term of a lower degree and we know that there are terms
      // of bigger degree, then put it right before this one.
      //
      else if (it->exponent<term.exponent && anybigger==true) {
        it--;
        terms.insert(it, term);
        return true;
      }
  }

  ////
  // We are minor than all others, so we can insert ourselves at the end.
  // FIXME: This could be done using terms.last()->exponent>term.exponent
  //
  terms.push_back(term);
  return true;
}

QString Polynomial::getAsString(void) {
  if(terms.size()==0) {
    return "NaN";
  }
  QString tmp=(terms.begin())->getAsString();
  for (std::vector<Term>::iterator it = terms.begin()+1; it != terms.end(); it++) {
      tmp+=" + "+it->getAsString();
  }
  return tmp;
}

Polynomial Polynomial::getDerivative(void) {
  if(terms.size()==0) {
    return (*this);
  }

  Polynomial tmp;
  Term tmp_term;

    for (std::vector<Term>::iterator it = terms.begin(); it != terms.end(); it++) {
      tmp_term = it->getDerivative();
      if(tmp_term.quocient!=0)
          tmp.terms.push_back(tmp_term);
    }

    return tmp;
}
Polynomial Polynomial::getIntegral(void) {
  if(terms.size()==0) {
    return (*this);
  }

  Polynomial tmp;
  Term tmp_term;

    for (std::vector<Term>::iterator it = terms.begin(); it != terms.end(); it++) {
      tmp_term = it->getIntegral();
      if(tmp_term.quocient!=0)
          tmp.terms.push_back(tmp_term);
    }

    return tmp;
}

double Polynomial::getImage(double object) {
  if(terms.size()==0) {
    return 0;
  }
  double sum=0;
  for (std::vector<Term>::iterator it = terms.begin(); it != terms.end(); it++) {
      sum+=it->getImage(object);
  }

  return sum;
}

double Polynomial::findZerobyNewtons(double startX, bool& success) {

  double currX=startX, oldX=0;
  Polynomial derivative = getDerivative();
  if(derivative.getImage(startX)==0) {
      success=false;
      return 0;
  }
  int n=0;
  long double e=5E-11;

  while (fabs(oldX-currX)>=e) {
      if(n>10000000) {
        success=false;
        return 0;
      }
      oldX=currX;
      //std::cout<<"Iterating with x_"<<n<<"="<<QString::number(currX, 'g', 15).toStdString()<<std::endl;
      currX=currX-getImage(currX)/derivative.getImage(currX);
      n++;
  }
  success=true;

  if(fabs(0-currX)<5E-11) {
      return 0; /** That's too close to zero not to be zero. FIXME: Should we do this? **/
  }
  return currX;
}

Polynomial doRuffiniDivision(long double a, Polynomial p) {
  /** int num=p.terms.size();
  Polynomial tmp;
  Term cur_term;
  cur_term=p.terms.at(0);
  tmp.terms.push_back(Term(cur_term.exponent-1, cur_term.quocient));

  unsigned int last_exponent=cur_term.exponent;
  long double carry=cur_term.quocient;



  for(int i=1; i<num; i++) {
      std::cout<<"We carry: "<<QString::number(carry, 'g', 15).toStdString()<<"\n";
      cur_term=p.terms.at(i);
      if(cur_term.exponent!=(last_exponent-1)) {
        last_exponent--;
        while(last_exponent>cur_term.exponent) {
            std::cout<<"iterating with last_exponent="<<last_exponent<<" and texp="<<cur_term.exponent<<"\n";
            if(last_exponent!=0) {
              std::cout<<"Adding0: "<<Term(last_exponent-1, (carry*a+0)).getAsString().toStdString()<<std::endl;
              tmp.addTerm(Term(last_exponent-1, carry=(carry*a+0)));
            }
            else {
              break; //FIXME: Is this ok? It tries to prevent pushing remainders. Seems unnecessary
              std::cout<<"Breaking shouldn't happen\n";
            }
            --last_exponent;
        }

      }

      std::cout<<"Adding: "<<Term((last_exponent=cur_term.exponent)-1, (carry*a+cur_term.quocient)).getAsString().toStdString()<<std::endl;
      tmp.addTerm(Term((last_exponent=cur_term.exponent)-1, fabs(0-(carry=(carry*a+cur_term.quocient)))<5E-11 ? 0 : carry));




  }

  ////
  // If the difference is below 5E-11, then we can admit that's 0
  // Also, it is a zero for sure if there's not term not multiplying by x and
  // a==0. That is, any polynomial whose terms are of degree >=1 will have
  // 0 as a root.
  //
  if(carry==0 || fabs(0-carry)<5E-11 || (a==0 && cur_term.exponent>=1) )
      std::cout<<"We have a zero!\n";
  else {
      std::cout<<"We have a remainder of: "<<QString::number(carry, 'g', 15).toStdString()<<"\n";
  }

  **/

  Polynomial tmp;

  long double carry=0;
  unsigned int last_exponent=p.terms.begin()->exponent;

  ////
  // Make this a full polynomial
  //
  for (std::vector<Term>::iterator it = p.terms.begin()+1; it != p.terms.end(); it++) {
      //std::cout<<"Iterating with it->exponent="<<it->exponent <<" and last_exponent="<<last_exponent<<std::endl;
      if(it->exponent!=last_exponent-1) {
        p.add0Term(last_exponent-1);
        last_exponent=(p.terms.begin()+1)->exponent;
        it=p.terms.begin()+1;
      }
      else
        last_exponent=it->exponent;

  }

  //std::cout<<p.getAsString().toStdString()<<std::endl;


  last_exponent=p.terms.begin()->exponent;
  tmp.addTerm(Term(last_exponent-1, carry=p.terms.begin()->quocient));
  for (std::vector<Term>::iterator it = p.terms.begin()+1; it != p.terms.end(); it++) {
      //std::cout<<"Adding: "<<Term((last_exponent=it->exponent)-1, (carry*a+it->quocient)).getAsString().toStdString()<<std::endl;
      tmp.addTerm(Term((last_exponent=it->exponent)-1, fabs(0-(carry=(carry*a+it->quocient)))<5E-11 ? 0 : carry));
  }

  ////
  // If the difference is below 5E-11, then we can admit that's 0
  // Also, it is a zero for sure if there's not term not multiplying by x and
  // a==0. That is, any polynomial whose terms are of degree >=1 will have
  // 0 as a root.
  //
  if(carry==0 || fabs(0-carry)<5E-11 || (a==0 && p.terms.end()->exponent>=1) ) {
      //std::cout<<"We have a zero!\n";
  } else {
      //std::cout<<"We have a remainder of: "<<QString::number(carry, 'g', 15).toStdString()<<"\n";
  }
  return tmp;
}

void Polynomial::printZeros(void) {
  Polynomial tmp=(*this);
  int n=0;
  for(;;) {
      bool success=false;
      long double zero=tmp.findZerobyNewtons(-100, success);
      if(!success)
        break;

      //else
      n++;
      std::cout<<"Zero "<<n<<": "<<QString::number(zero, 'g', 15).toStdString()<<std::endl;
      tmp=doRuffiniDivision(zero, tmp);
      std::cout<<"Now working with "<<tmp.getAsString().toStdString()<<std::endl;
  }
}

int main(void) {
  Term term1(4,-5);
  Term term2(2,-3);
  Term term3(5,3);
  Term term4(0,5); //Check case of Term term4(0,-2);  at end (and of adding term4 0,5 in the beggining and getting nans)
  Term term5(6,3);

  Polynomial poly, poly2;
  poly.addTerm(term1);
  poly.addTerm(term2);
  poly.addTerm(term3);
  poly.addTerm(term4);
  std::cout<<"The image of -1.01775 in "<<poly.getAsString().toStdString()<<" is "<<QString::number(poly.getImage(-1.01775105728273), 'g', 15).toStdString()<<std::endl;
  //bool success=false;
  /**long double zero= poly.findZerobyNewtons(-0.0003, success);
  std::cout<<"Success="<<success<<std::endl;
  std::cout<<"One of the zeroes is:"<<QString::number(zero, 'g', 15).toStdString()<<std::endl;
  poly2=doRuffiniDivision(zero, poly);
  std::cout<<"The division outputs: "<<poly2.getAsString().toStdString()<<std::endl;
  **/
  poly.printZeros();
  //std::cout<<term1.getAsString().toStdString()<<std::endl;
  //std::cout<<term1.getInfoAsString().toStdString()<<std::endl;
  //std::cout<<"The derivative of "<<poly.getAsString().toStdString()<<" is: "<<(poly2=poly.getDerivative()).getAsString().toStdString()<<std::endl;
  //std::cout<<"For that to be true, we must make sure that the integral of the derivative is the same as the original. See for yourself:\n";
  //std::cout<<"The integral of "<<poly2.getAsString().toStdString()<< " is: "<<poly2.getIntegral().getAsString().toStdString()<<std::endl;
  //std::cin.get();
  return 0;
}

---

Only glanced it, don't have much time, but this is wrong:

Code:

---

unsigned int last_exponent = terms.end()->exponent;

---

terms.end() returns the iterator PASSED the last element, not to the last element. So terms.end() doesn't point to a valid element....

Thank you very much, that was probably it. My knowledge of iterators is still very limited, as I've never really had a need to use it in my small projects. It was quite stupid of me not to properly read the documentation, though.

I've just recompiled it all and debugged it. It works like a charm. Now to optimize everything :)

Thank you very, very much,

Jorl17