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;
}