Pour ceux qui en veulent encore plus, je vous invite à aller faire un tour sur le forum "Autres langages" qui, grâce à bluestorm (et d'autres) propose depuis aujourd'hui un atelier très intéressant. Il s'agit d'écrire un interpréteur pour un mini langage de programmation.
Il y a différentes difficultés et si vous allez jusqu'au bout, vous aurez déjà un vrai langage bien conséquent et capable de beaucoup plus que ce qu'on pourrait croire au premier abord.
Si vous avez des questions d'ordre C++, n'hésitez pas à les poser sur ce forum. Vous recevrez certainement plus d'aide. Pour le reste, postez là-bas pour garder la cohérence du sujet.
Un des outils souvent utilisés par les administateurs de machines sous Linux est le petit programme appelé "diff". Il permet d'afficher la différence entre deux fichiers en spécifiant quelles lignes sont différentes dans les 2 fichiers passés en argument. Ce qui donne par exemple (sans aucune option) avec les deux fichiers suivants,
#include <iostream>
using namespace std;
int main()
{
cout << "Hello, World !" << endl;
return 0;
}
#include <iostream>
int main()
{
std::cout << "Hello, World !" << std::endl;
return 0;
}
//Fin du programme
le résultat suivant:
diff Fichier1 Fichier2
2c2
< using namespace std;
---
>
6c6
< cout << "Hello, World !" << endl;
---
> std::cout << "Hello, World !" << std::endl;
8a9,10
>
> //Fin du programme
Il faut lire le résultat de la manière suivante:
2c2 veut dire que la ligne 2 de chaque fichier est différente et ce qui suit correspond aux lignes dans chacun des fichiers séparées par "---".
Idem pour 6c6 qui indique que la ligne 6 a été changée. La ligne 8a9,10 indique que les lignes 9 et 10 ont été ajoutées.
Si on appelle le programme avec les fichiers dans l'ordre opposé, on obtient:
diff fichier2 fichier1
2c2
<
---
> using namespace std;
6c6
< std::cout << "Hello, World !" << std::endl;
---
> cout << "Hello, World !" << endl;
9,10d8
<
< //Fin du programme
Le début est le même, mais cette fois on voit que les lignes 9 et 10 ont été détruites dans le fichier de droite.
Vous remarquerez que les lignes identiques ne sont jamais affichées.
L'exercice
Le but de cet exercice est de réaliser un programme de diff. Fonctionnant sur le même principe que celui de Linux.
Votre programme devra recevoir deux fichiers (les chemins) en argument et montrer les différences entre les deux fichiers en indiquant si les lignes ont été ajoutées (a), changées (c) ou détruites (d) avec bien sûr le numéro des lignes différentes et le contenu de la ligne dans chacun des 2 fichiers en utilisant la notation < et > pour indiquer de quel fichier provient la ligne affichée.
Vous avez jusqu'au 31 mai pour soumettre vos réponses à Réponse_Exercices.
Ah, même si j'ai mis du temps à comprendre, c'est plus facile que les polynômes . Par contre, j'ai pas vraiment compris le truc des lignes détruites... c'est quand le fichier entré en premier a moins de lignes que l'autre ? J'ai pas vraiment compris l'affichage... pareil que 8a9,10...
Si les fichiers n'ont aucune différence, il n'y a rien d'affiché, si ?
Pour ce qui est de savoir les lignes qui ont été ajoutées/détruites, je pense qu'il faut compter le nombre de lignes dans chacun des deux fichiers : s'ils n'en ont pas le même nombre, c'est qu'il y a des lignes qui ont été ajoutées ou détruites, sinon pas de souci.
Je n'ai pas reçu beaucoup de réponses (3) pour l'exercice sur les polynômes.
La base de la classe
Comment re présenter un polynôme ? Le plus simple est certainement de stocker les coefficients. Le plus facile pour cela est d'utiliser un vector<double>.
On peut également utiliser une map ou une list si on veut optimiser certains cas particuliers du type <math>\(x^{42} + 2\)</math> où il y a beaucoup de termes intermédiaires nuls.
Nous avons donc la base suivante:
class Polynome{
//...
private:
std::vector<double> m_termes;
};
Constructeurs
Les constructeurs prenant en argument un paramètre de type double et le constructeur par défaut peuvent être créés ensembles.
Polynome(double a = 0.)
:m_termes(1,a) //On remplit le vector avec un terme valant a
{}
Un constructeur recevant un vector peut également s'avérer utile. Comme l'utilisateur peut avoir entré un tableau de type 3,4,0,5,0,0,0 avecdes 0s inutiles, il faut les simplifier. On crée donc une fonction simplifie() qui fait cela.
Le plus simple ici est d'utiliser les fonctions membres du vector. On a ainsi les versions [] et at() qui vérifie l'indice. Comme toujours, il y a une version const et non-const.
Vérifier si deux polynomes sont égaux consiste à vérifier l'égalité des termes. Pour se simplifier la vie, on utilise l'opérateur de comparaison de vector.
Comme toujours, on déclare les opérateurs +=,-= et *= dans la classe puis les versions +,- et * à l'extérieur. Au niveau des calculs, on utilise la définition (par exemple de Wikipédia).
La division euclidienne donne accès à 2 informations, le quotient et le reste. On peut ainsi créer un opérateur % et / en créant une seule fonction.
A nouveau, on utilise l'algorithme usuel. Comme la fonction n'a pas besoin d'un objet pour fonctionner, on peut la mettre en statique.
typedef std::pair<Polynome, Polynome > Paire;
static Paire division(const Polynome& dividende, const Polynome& diviseur)
{
if (diviseur == 0.)
throw std::runtime_error("Polynome: Division par 0");
if (dividende.degre() < diviseur.degre())
return Paire(0, dividende);
if (dividende == diviseur)
return Paire(1,0);
Polynome reste(dividende);
Polynome quotient;
while (reste.degre() > diviseur.degre())
{
std::vector<double> vec(reste.degre() - diviseur.degre() + 1,0);
vec.back() = reste.m_termes.back() / diviseur.m_termes.back();
Polynome p(vec);
reste -= p*diviseur;
reste.simplifie();
quotient += p;
}
return Paire(quotient,reste);
}
Avec cela, on a terminé les niveaux 1 et 2. Voici le code complet avec quelques fioritures en plus.
Pour le niveau 3, il s'agissait d'ajouter des templates. Cela ne fait pas grand chose de nouveau. Voici une version complète avec des commentaires Doxygen (en anglais).
#ifndef __POLYNOMIAL_HPP_INCLUDED__
#define __POLYNOMIAL_HPP_INCLUDED__
/**
* \file Polynom.hpp
* \author Nanoc
* \brief A Polynomial class
* \date 19 April 2009
* \version 1.1
*
* The operators are implemented naively. Could be improved
* \todo Add a solve() function
*/
#include <vector>
#include <cmath>
#include <stdexcept>
#include <ostream>
/**
* \class Polynomial<>
*
* \brief A class to represent polynomials of any degree
*
* All the usual operators are overloaded. One can evaluate the polynomial at any point or compute its derivative
* Access to all terms is also granted
*
*/
template<typename T>
class Polynomial
{
public:
// Constructors -------------------------------------------------------------------------------
/**
* \brief Constructor. Buils a Polynomial from a single value
*/
Polynomial(const T& a = static_cast<T>(0))
:m_terms(1,a)
{}
/**
* \brief Constructor. Builds a Polynomial from a vector of terms
*/
Polynomial(const std::vector<T>& tab)
:m_terms(tab)
{
simplify();
}
/**
* \brief Constructor. Buils a Polynomial from a "C-like" array
*/
template<size_t N>
Polynomial(const T (&tab)[N])
:m_terms(&tab[0], &tab[N])
{
simplify();
}
// Term access ------------------------------------------------------------------------------
/**
* \brief Access to a term (read only)
*/
inline const T& operator[](size_t n) const
{
return m_terms[n];
}
/**
* \brief Access to a term
*/
inline T& operator[](size_t n)
{
return m_terms[n];
}
/**
* \brief Access to a term with rangecheck (read only)
*/
inline const T& at(size_t n) const
{
return m_terms.at(n);
}
/**
* \brief Access to a term with rangecheck
*/
inline T& at(size_t n)
{
return m_terms.at(n);
}
// Degree -----------------------------------------------------------------------------------
/**
* \brief Returns the degree of the polynomial
*/
inline size_t degree() const
{
return m_terms.size() - 1;
}
/**
* \brief Returns the length of the internal array (should be degree()+1)
*/
inline size_t length() const
{
return m_terms.size();
}
// Evaluation -----------------------------------------------------------------------------------
/**
* \brief Evaluates the polynom at a given point x
*/
template <typename S>
S evaluate(const S& x) const
{
S result = static_cast<S>(0);
for (size_t i(0); i<m_terms.size(); ++i)
{
result += m_terms[i] * power(x,i);
}
return result;
}
/**
* \brief Evaluates the polynom at a given point x
*/
template <typename S>
inline S operator()(const S& x) const
{
return evaluate(x);
}
//Operators ----------------------------------------------------------------------------
/**
* \brief Sum term by term of two polynomials
*/
Polynomial<T>& operator+=(const Polynomial<T>& q)
{
size_t i = 0;
for (; i<q.length() && i<m_terms.size(); ++i)
{
m_terms[i] += q[i];
}
for (; i<q.length() ; ++i)
m_terms.push_back(q[i]);
return *this;
}
/**
* \brief Difference term by term of two polynomials
*/
Polynomial<T>& operator-=(const Polynomial<T>& q)
{
size_t i = 0;
for (; i<q.length() && i<m_terms.size(); ++i)
{
m_terms[i] -= q[i];
}
for (; i<q.length() ; ++i)
m_terms.push_back(-q[i]);
return *this;
}
/**
* \brief Product of two Polynomials
*/
Polynomial<T>& operator*=(const Polynomial<T>& q)
{
if (q == Polynomial<T>(0) || *this == Polynomial<T>(0))
{
m_terms.assign(1,static_cast<T>(0));
return *this;
}
const size_t len(q.degree() + degree() + 1);
std::vector<T> vec(len,static_cast<T>(0));
for (size_t k(0); k<len; ++k)
for (size_t i(0); i< this->length();++i)
{
const size_t j = k-i;
if (j<q.length())
vec[k] += m_terms[i]*q[j];
}
*this = Polynomial<T>(vec);
return *this;
}
/**
* \brief Quotient of the Euclidian division
*/
inline Polynomial<T>& operator/=(const Polynomial<T>& q)
{
*this = division(*this,q).first;
return *this;
}
/**
* \brief Remainder of the Euclidian division
*/
inline Polynomial<T>& operator%=(const Polynomial<T>& q)
{
*this = division(*this,q).second;
return *this;
}
//Derivative -------------------------------------------------------------------------
/**
* \brief Returns the derivative of the polynom
*/
Polynomial<T> derivative() const
{
if (degree() == 0)
return 0;
Polynomial<T> ret(*this);
for (size_t i(0); i<ret.degree(); ++i)
ret[i] = ret[i+1] * (i+1);
ret.m_terms.pop_back();
return ret;
}
private:
//Simplification -----------------------------------------------------------------------
/**
* \brief Simplifies the polynom by removing all 0-terms after the highest non-zero term
*/
void simplify()
{
size_t counter = 0;
typedef typename std::vector<T>::reverse_iterator Iterator;
for (Iterator it = m_terms.rbegin(); it != m_terms.rend() && *it == static_cast<T>(0) ; ++it)
{
++counter;
}
m_terms.resize(m_terms.size() - counter);
}
/**
* \brief Notation simplification
*/
typedef std::pair<Polynomial<T>, Polynomial<T> > Pair;
/**
* \brief Returns a Pair containing the quotient and the remainder of a division
*/
static Pair division(const Polynomial& dividend, const Polynomial& divisor)
{
if (divisor == Polynomial<T>(static_cast<T>(0)))
throw std::runtime_error("Polnomial: Division by 0");
if (dividend.degree() < divisor.degree())
return Pair(0, dividend);
if (dividend == divisor)
return Pair(1,0);
Polynomial<T> remainder(dividend);
Polynomial quotient;
while (remainder.degree() > divisor.degree())
{
std::vector<T> vec(remainder.degree() - divisor.degree() + 1,0);
vec.back() = remainder.m_terms.back() / divisor.m_terms.back();
Polynomial<T> p(vec);
remainder -= p*divisor;
remainder.simplify();
quotient += p;
}
return Pair(quotient,remainder);
}
//Attributes-----------------------------------------------------------------------------
std::vector<T> m_terms; /*!< The coefficients */
//Friends--------------------------------------------------------------------------------
friend bool operator==(const Polynomial<T>& p, const Polynomial<T>& q);
};
// Stream -----------------------------------------------------------------------------------
/**
* \brief Writes the polynom in a stream
*/
template <typename T>
std::ostream& operator<<(std::ostream& stream, const Polynomial<T>& poly)
{
if (poly.degree() == 0 && poly[0] == 0)
{
stream << static_cast<T>(0);
return stream;
}
if (poly[poly.length() - 1] < static_cast<T>(0))
stream << "- ";
for (size_t i(1);i<poly.length();++i)
{
if (poly[poly.length() - i] != static_cast<T>(0))
stream << std::abs(poly[poly.length() - i]) << "x^" << poly.length() - i;
if (poly[poly.length() - (i+1)] > static_cast<T>(0))
stream << " + ";
else if (poly[poly.length() - (i+1)] < static_cast<T>(0))
stream << " - ";
}
if (poly[0] != static_cast<T>(0))
stream << std::abs(poly[0]);
return stream;
}
//Comparison operators ------------------------------------------------------------------------
/**
* \brief Compares two polynoms to see if they are equal
*/
template <typename T>
bool operator==(const Polynomial<T>& p, const Polynomial<T>& q)
{
return p.m_terms == q.m_terms;
}
/**
* \brief Compares two polynoms to see if they are different
*/
template <typename T>
inline bool operator!=(const Polynomial<T>& p, const Polynomial<T>& q)
{
return !(p == q);
}
// Binary Operators ---------------------------------------------------------------------------
template<typename T>
inline Polynomial<T> operator+(const Polynomial<T>& p, const Polynomial<T>& q)
{
return Polynomial<T>(p) += q;
}
template<typename T>
inline Polynomial<T> operator-(const Polynomial<T>& p, const Polynomial<T>& q)
{
return Polynomial<T>(p) -= q;
}
template<typename T>
inline Polynomial<T> operator*(const Polynomial<T>& p, const Polynomial<T>& q)
{
return Polynomial<T>(p) *= q;
}
template<typename T>
inline Polynomial<T> operator/(const Polynomial<T>& p, const Polynomial<T>& q)
{
return Polynomial<T>(p) /= q;
}
template<typename T>
inline Polynomial<T> operator%(const Polynomial<T>& p, const Polynomial<T>& q)
{
return Polynomial<T>(p) %= q;
}
//Unary operators -------------------------------------------------------------------------------
template<typename T>
Polynomial<T> operator-(const Polynomial<T>& p)
{
Polynomial<T> result(p);
for (size_t i(0); i<result.length(); ++i)
result[i] = -result[i];
return result;
}
template<typename T>
Polynomial<T> operator+(const Polynomial<T>& p)
{
Polynomial<T> result(p);
for (size_t i(0); i<result.length(); ++i)
result[i] = +result[i];
return result;
}
//Logical operators ----------------------------------------------------------------------------
template<typename T>
inline bool operator!(const Polynomial<T>& p)
{
return p == Polynomial<T>(0);
}
//Derivative -----------------------------------------------------------------------------------
/**
* \brief Returns the derivative of the Polynomial passed
*/
template<typename T>
inline Polynomial<T> derivative(const Polynomial<T>& p)
{
return p.derivative();
}
//Useful functions -----------------------------------------------------------------------------
/**
* \brief Returns the n-th power of a given x
*/
template<typename T>
T power(const T& x, size_t n)
{
T result = static_cast<T>(1);
for (size_t i(0); i<n ; ++i)
result *= x;
return result;
}
//Simplification typedef -----------------------------------------------------------------------
/**
* \typedef Poly
* A simpler name for standard polynomial using double as coefficients
*/
typedef Polynomial<double> Poly;
/**
* \typedef Polyd
* A simpler name for standard polynomial using double as coefficients
*/
typedef Polynomial<double> Polyd;
/**
* \typedef Polyf
* A simpler name for standard polynomial using float as coefficients
*/
typedef Polynomial<float> Polyf;
#endif //__POLYNOM_HPP_INCLUDED__
Je ne t'avais pas envoyé ce que j'avais fait, j'avais beaucoup de choses que je voulais rédiger (expliquer comment on écrit les opérateurs), plus une expérience à mener au sujet des expressions template. Mais j'ai eu d'autres choses à faire en court de route.
Je peux t'envoyer le code, ou le poster ici, comme tu préfères.
PS:
- polynomes et matrices sont deux "exceptions", où il est préférable d'implémenter *= en fonction de *.
- il suffit d'utiliser operator==(vector, vector) pour la comparaison.
Tiens, j'avais pas vu cet exo ! Je suis un peu dans le même cas que lmghs, j'ai réalisé avec un ami un programme de résolution d'(in)équation ( enfin, il manque l'IHM en php ... ) . Si ça intéresse : http://gmath.servhome.org/index.php .
3 choses :
a) Le logiciel en ligne n'est pas activé donc chargement infini. Une autre version est en développement mais aucune date n'est annoncée pour sa mise en utilisation libre.
b) Dans la doc, ce qui est utile ici se réduit à deux classes sauf erreur de ma part ( ça fait un bail que je ne l'ai pas retouché ) : BasePolynome et Polynome.
c) La division polynomiale est prise en charge.
@hiura: intéressant, mais ça ne marche pas chez moi ;-) Ou alors j'ai pas compris comment enrter l'équation. Quel algorithme utilisez-vous ? Newton + Horner ? Bézout ? Tschirnhausen ?
@hiura: intéressant, mais ça ne marche pas chez moi ;-)
Oui c'est normal, on n'a pas les droits pour exécuter l'application actuellement.
Citation : Nanoc
Quel algorithme utilisez-vous ? Newton + Horner ? Bézout ? Tschirnhausen ?
Je ne suis pas responsable de cette partie de l'application mais je dirais Newton. ( En trois ans de développement, et 6 mois d'inactivisme j'ai oublié certains détails :s . ) Je vais demander. Je peux simplement affirmé que pour les degrés 0 ( x € R ), et 1 ( x = -b / a ) et 2 ( Viète ) c'est d'autres méthodes plus simples qui sont employés.
[ Comme ça je dirais Newton car c'est le seul nom que je reconnais avec Bézout, mais comme quand on l'a codé on ne connaissait même pas les dérivées ... ]
PS : la conversion string -> polynôme n'est pas terrible, il aurait fallu qu'on connaisse les techniques de grammaire pour faire mieux.
C'est en trois parties, et je vais juste vous donner la solution template, la non template est sans grand intérêt (et celle reposant sur les expressions template tout juste démarrée).
Pour pouvoir écrire: "const Polynomial A = (X^5) - (X^4) - (X^3) + (3*X^2) - (2*X);", je suis passé par une classe utilitaire pour les monômes. Et je continue avec deux très vilains hacks:
- le premier avec l'opérateur xor que je détourne pour les puissances -- opérateur dont la priorité est inadaptée à ce que l'on veut faire, d'où le parenthésage systématique car nécessaire dans l'expression ci-dessus.
- un second pour utiliser juste X sans avoir de .cpp, et en jouant avec de très vilains #define -- bref, les enfants ne faites pas ça au travail, ni même à la maison.
- l'évaluation repose sur la méthode de Hörner comme il se doit.
Ensuite, on peut remarquer que :
- pour des raisons de DRY et de symétrie, + et - sont définis en termes de += et -=, tandis que *= et /= sont définis en terme de * et /.
- les retours par valeur se font constants (pour refuser (2*X^2)*=42 -- à l'image de ce que fait le compilo pour les types entiers)
- je ne fournis pas de version non constante de l'opérateur[], ainsi nul besoin de sortir des proxys pour différencier l'accès RW (read-write) de l'accès RO (read-only) et s'enquiquiner à faire des trim à la volée.
Dans les dépendances, il faut rajouter array_helpers.hpp:
Bon bah je m'installe , mon code n'est pas aussi pro que celui de lmghs, ne gère pas les templates (bah oui j'en suis pas encore là dans le cours ) mais je le préfère pour les outils mathématiques élémentaires qu'il apporte.
#include "polynome.h"
#include <cmath>
#include <iomanip>
using namespace std;
/********************************************************
Constructeurs
*******************************************************/
// Constructeur de copie
Polynome::Polynome(const Polynome &P){
m_coeffs = P.m_coeffs;
sturm_computed = P.sturm_computed;
Sturm = P.Sturm;
}
// Par défault : Créer le polynome nul
Polynome::Polynome(){
m_coeffs.push_back(0.);
sturm_computed = false;
}
// Via un vecteur : les coefficients du polynôme sont recopié du vecteur
Polynome::Polynome(const std::vector<double> &coeffs){
m_coeffs = coeffs;
sturm_computed = false;
}
// Via une liste variable : Polynome(nb de coeffs, a0, a1, a2, ...)
// ex : Polynome(3, 1., -2., 4.) crée le polynôme 4x^2-2x+1
// La liste variable est constituée de double !!! Polynome(2, 1, 2) ne renverra pas 2x+1 !!!
Polynome::Polynome(int n, const double x, ...){
sturm_computed = false;
if (n <= 0)
return;
va_list varg; // Variable stockant la liste
va_start(varg, x); // On initialise la variable
m_coeffs.push_back(x);
n--;
do {
m_coeffs.push_back(va_arg(varg, double));
}while(--n!=0);
va_end(varg); // On détruit la liste variable
}
// Via un tableau de double, ceci permet de profiter de l'auto-cast int -> double
// et donc de pouvoir écrire Polynome(3, (double[]) {1, 0, 1.}) pour le polynôme x^2+1
// ou bien double coef[] = {1, 0, 1.}; Polynome(3, coef)
Polynome::Polynome(const int n, const double *liste){
sturm_computed = false;
if (n <= 0)
return;
for(int i = 0 ; i < n ; i++)
m_coeffs.push_back(liste[i]);
}
/********************************************************
******************* Destructeur **********************
*******************************************************/
Polynome::~Polynome(){
// Rien de bien palpitant
}
/********************************************************
****************** Opérateurs ************************
*******************************************************/
// Opérateur d'affectation
Polynome Polynome::operator=(const Polynome &P){
m_coeffs = P.m_coeffs;
sturm_computed = P.sturm_computed;
Sturm = P.Sturm;
return *this;
}
Polynome Polynome::operator+=(const Polynome &P)
{
int n = this->degre();
int m = P.degre();
int i;
if (m <= n)
{
// On rajoute un polynome de degré plus petit ou égal
for(i = 0 ; i <= m ; i++)
m_coeffs[i] += P.m_coeffs[i];
}
else
{
// On rajoute un polynôme de degré plus grand
m_coeffs.reserve(m+1);
for(i = 0 ; i <= n ; i++)
m_coeffs[i] += P.m_coeffs[i];
for(i = n + 1 ; i<= m ; i++)
m_coeffs[i] = P.m_coeffs[i];
}
return *this;
}
Polynome Polynome::operator-=(const Polynome &P)
{
return (*this)+=-P;
}
Polynome Polynome::operator*=(const Polynome &P)
{
*this = *this * P;
return *this;
}
Polynome Polynome::operator/=(const Polynome &P) throw(std::domain_error)
{
if(P.degre() == -1)
throw std::domain_error("Division par zéro");
*this = *this / P;
return *this;
}
Polynome Polynome::operator%=(const Polynome &P) throw(std::domain_error)
{
if(P.degre() == -1)
throw std::domain_error("Division par zéro");
*this = *this % P;
return *this;
}
// Ajoute 2 polynômes entre eux
Polynome Polynome::operator+(const Polynome &P) const{
return Polynome(*this)+=P;
}
// Retranche 2 polynômes
Polynome Polynome::operator-(const Polynome &P) const{
return Polynome(*this)+=-P;
}
// Renvoie l'opposée du polynôme
Polynome Polynome::operator-() const{
int n = m_coeffs.size();
Polynome res;
res.m_coeffs.resize(n);
for(int i = 0 ; i < n ; i++)
res.m_coeffs[i] = - m_coeffs[i];
return res;
}
// Multiplie 2 polynômes entre eux : X^i * X^j = X^(i+j)
Polynome Polynome::operator*(const Polynome &P) const{
int n = this->degre();
int m = P.degre();
if (n < 0 || m < 0)
return Polynome();
vector<double> M(n+m+1,0);
for(int i = 0 ; i<= n ; i++)
for(int j = 0 ; j <= m ; j++)
M[i+j] += m_coeffs[i]*P.m_coeffs[j];
return Polynome(M);
}
// Renvoie le quotient de la division euclidienne
Polynome Polynome::operator/(const Polynome &P) const throw(std::domain_error){
if(P.degre() == -1)
throw std::domain_error("Division par zéro");
Polynome Q,R;
div_eucli(*this, P, Q, R);
return Q;
}
// Renvoie le reste de la division euclidienne
Polynome Polynome::operator%(const Polynome &P) const throw(std::domain_error){
if(P.degre() == -1)
throw std::domain_error("Division par zéro");
Polynome Q;
Polynome R;
div_eucli(*this, P, Q, R);
return R;
}
// Test l'égalité de 2 polynômes
bool Polynome::operator==(const Polynome &P) const{
int n = this->degre();
int m = P.degre();
if (n != m)
return false;
for(int i = 0 ; i <=n ; i++)
if(m_coeffs[i] != P.m_coeffs[i])
return false;
return true;
}
// Test la différence de 2 polynômes
bool Polynome::operator!=(const Polynome &P) const{
return !(*this == P);
}
// Confie à Poly.affiche la tache cout << P
ostream &operator<<(ostream &out, const Polynome &Poly){
Poly.affiche(out);
return out;
}
// Evalue un polynôme en 1 point
double Polynome::operator()(const double x) const{
return this->eval(x);
}
/********************************************************
****************** Arithmétique **********************
*******************************************************/
// Effectue la division euclidienne de A par B, le quotient va dans Q et le reste dans R
void Polynome::div_eucli(const Polynome &A, const Polynome &B, Polynome &Q, Polynome &R) throw(std::domain_error){
if(B.degre() == -1)
throw std::domain_error("Division par zéro");
int n = A.degre();
int m = B.degre();
if (m < 0)
return;
if (n-m < 0){
Q = Polynome();
R = A;
return;
}
R.m_coeffs = A.m_coeffs;
Q.m_coeffs.resize(n-m+1);
double pivot = 1/B.m_coeffs[m], r;
for(int i = n-m ; i >= 0 ; i--){
r = R.m_coeffs[i+m] * pivot;
Q.m_coeffs[i] = r;
if(r != 0)
for(int j = i ; j<=i+m ; j++){
R.m_coeffs[j] -= r * B.m_coeffs[j-i];
if ( abs(R.m_coeffs[j]) < POLYNOME_BIG_EPSILON)
R.m_coeffs[j] = 0;
}
R.m_coeffs[i+m] = 0;
}
}
// Renvoie le pgcd de 2 polynômes
Polynome Polynome::pgcd(const Polynome &A, const Polynome &B){
Polynome A1 = A, B1 = B, Q, R;
do{
R = A1 % B1;
A1 = B1;
B1 = R;
} while(R.degre() != -1);
return A1;
}
// Renvoie la dérivée du polynôme
Polynome::Polynome Polynome::derive(unsigned int d) const{
int n = this->degre();
if (n+1-(int)d <= 0)
return Polynome();
vector<double> res(n+1-d);
int p = 1;
for(unsigned int i = 1 ; i <= d ; i++)
p *= i;
for(int i = 0; i <= n-(int)d ; i++){
res[i] = p*m_coeffs[i+d];
p = p * (d+i+1)/(i+1);
}
return Polynome(res);
}
// Renvoie la primitive du polynôme nulle en 0
Polynome::Polynome Polynome::integre(unsigned int d) const{
int n = this->degre();
if (n+1+d <= 0)
return Polynome();
vector<double> res(n+1+d);
int p = 1;
for(unsigned int i = 1 ; i <= d ; i++)
p *= i;
for(unsigned int i = 0 ; i < d ; i++)
res[i] = 0;
for(int i = 0 ; i <= n ; i++){
res[i+d] = m_coeffs[i]/p;
p = p * (d+i+1)/(i+1);
}
return Polynome(res);
}
// Rend le polynôme unitaire (son coefficient de plus haut degré vaudra 1)
void Polynome::unitaire(){
int n = this->degre();
if (n == -1)
return;
double div = 1 / m_coeffs[n];
for(int i = 0 ; i <= n ; i++)
m_coeffs[i] *= div;
return;
}
/********************************************************
********************* Racines *************************
*******************************************************/
// Renvoie le nombre de racines réelles du polynômes (Théorème de Sturm)
int Polynome::nb_roots(){
if( !sturm_computed )
sturm_compute();
return abs(nb_chg_signe(-M) - nb_chg_signe(M));
}
// Renvoie le nombre de racines réelles du polynômes dans l'intervalle [a,b] (Théorème de Sturm)
int Polynome::nb_roots(double a, double b){
if( !sturm_computed )
sturm_compute();
return abs(nb_chg_signe(a) - nb_chg_signe(b));
}
// Calcul les racines réelles d'un polynôme
vector<double> Polynome::roots(){
int nb_racines = nb_roots();
if (nb_racines == 0)
return vector<double>();
// Recherches les intervalles sur lesquels les racines sont isolés
// Cette recherche se fait via une dichotomie utilisée sur le théorème de Sturm
double *bornes = new double[nb_racines + 1]; // Pour n racines, il y a n+1 intervalles
bool *borne_fixe = new bool[nb_racines + 1];
for(int i=0;i<nb_racines+1;i++)
borne_fixe[i] = false;
vector<double> racines;
bornes[0] = -M; bornes[nb_racines] = M; // On fixe les bornes limites [-M,M]
borne_fixe[0] = true; borne_fixe[nb_racines] = true;
for(int i = 1 ; i < nb_racines ; i++){
// On recherche le premier intervalle qui contient plusieurs racines pour le scindé en 2
// Derniere borne fixe
int gauche = 0, droite, rac_gauche, rac_droite, rac_milieu, nb_rac;
double x_g, x_d, milieu;
bool continuer;
while(borne_fixe[gauche])
gauche++;
gauche--;
// Cette erreur ne devrait jamais se produire
if (gauche > nb_racines){
std::cerr << "Erreur lors du scindage de l'intervalle (recherche de la borne gauche\n";
std::cerr << left << setw(20) << "Valeur des bornes" << "Fixé\n";
for(int i = 0 ; i <= nb_racines ; i++)
std::cerr << left << setw(20) << bornes[i] << borne_fixe[i] << "\n";
std::cerr << flush;
return vector<double>(0);
}
// Prochaine borne fixe
droite = gauche + 1;
while(!borne_fixe[droite])
droite++;
// Cette erreur ne devrait jamais se produire
if (droite > nb_racines){
std::cerr << "Erreur lors du scindage de l'intervalle (recherche de la borne droite\n";
std::cerr << left << setw(20) << "Valeur des bornes" << "Fixé\n";
for(int i = 0 ; i <= nb_racines ; i++)
std::cerr << left << setw(20) << bornes[i] << borne_fixe[i] << "\n";
std::cerr << flush;
return vector<double>(0);
}
x_g = bornes[gauche];
x_d = bornes[droite];
rac_gauche = nb_chg_signe(bornes[gauche]);
rac_droite = nb_chg_signe(bornes[droite]);
// On cherche à couper l'intervalle [x_g, x_d] en 2 intervalles [x_g, milieu] et [x_d, milieu]
// où chacun des 2 contient au moins une racine
continuer = true;
while(continuer){
milieu = (x_g + x_d)/2;
rac_milieu = nb_chg_signe(milieu);
if( rac_milieu == rac_gauche) // Il n'y a pas de racines entre x_g et milieu
x_g = milieu; // donc on décale la borne gauche
else if( rac_milieu == rac_droite) // Il n'y a pas de racines entre milieu et x_d
x_d = milieu; // donc on décale la borne droite
else
continuer = false; // On a réussi à couper l'intervalle
}
// On scinde l'intervalle
nb_rac = abs(rac_gauche - rac_milieu);
bornes[gauche + nb_rac] = milieu;
borne_fixe[gauche + nb_rac] = true;
// Réduit la taille de l'intervalle total [-M, M] pour accélerer la convergence de la méthode de newton
if (gauche == 0)
bornes[gauche] = x_g;
if (droite == nb_racines)
bornes[droite] = x_d;
}
Polynome calc = *this / Polynome::pgcd(*this, this->derive()); // Permet de transformer toutes les racines de multiplicité > 1 en racines de multiplicité 1
for(int i = 0 ; i < nb_racines ; i++){ // Ce qui permet une convergence quadratique de la méthode de newton
int iter = 0; // Nombre d'itérations pour le calcul d'une racine
double rac;
double borne_gauche = bornes[i], borne_droite = bornes[i+1]; // Bornes de recherche de la racine
bool calcul = true; // On effectue le calcul
while(calcul){
try{
iter++;
rac = calc.newton((borne_gauche + borne_droite)/2);
if (rac >= bornes[i] && rac <= bornes[i+1]){ // Si la racine se situe dans l'intervalle de recherche
racines.push_back(rac); // On l'empile sur le vecteur de solutions
calcul = false; // Le calcul s'est effectué jusqu'au bout
calc = calc / Polynome(2, -racines[i], 1.); // Divise le polynôme de recherche par (X-rac) pour supprimer la racine (enlève un bassin d'attraction pour la méthode de newton)
}
}
catch(const std::domain_error &e){
// La dérivée s'est annulée, laisse calcul à true
}
catch(const std::runtime_error &e){
// L'algorithme tourne en boucle, laisse calcul à true
}
if (calcul){ // Si l'on doit continuer à chercher la racine
calc.dicho(borne_gauche, borne_droite, POLYNOME_NB_ITERATIONS_DICHO); // On divise l'intervalle de recherche par 2^POLYNOME_NB_ITERATIONS_DICHO
}
if (iter > POLYNOME_NB_ITERATIONS_RECHERCHE) // Si l'on a effectué trop d'itérations, on abandonne la recherche
calcul = false;
}
}
delete[] bornes;
delete[] borne_fixe;
return racines;
}
// Construit la suite de Sturm
void Polynome::sturm_compute(){
Polynome A = *this;
Polynome B = A.derive();
Polynome R;
A.unitaire();
Sturm.clear();
Sturm.push_back(A);
int r = 1;
do{
Sturm.push_back(B);
r++;
R = - (A % B);
A = B;
B = R;
} while(R.degre() != -1);
// Si le polynôme a des racines multiples on réduit leur multiplicité à 1
if (Sturm[r-1].degre() != 0)
for(int i = 0 ; i < r ; i++)
Sturm[i] = Sturm[i] / Sturm[r-1];
// Toutes les racines seront dans [-M, M]
int degre = this->degre();
M = 0;
double coeff = m_coeffs[degre];
for(int i = 0 ; i <= degre ; i++)
M += abs(m_coeffs[i] / coeff);
if (M < 1)
M = 1;
sturm_computed = true;
}
// Indique si la suite de Sturm est construite ou non
bool Polynome::get_sturmcomputed(){
return sturm_computed;
}
// Calcul le nombre de changement de signe utilisé dans le théorème de Sturm
int Polynome::nb_chg_signe(double x){
int nb = 0;
int n = Sturm.size();
// On recherche le premier terme non nul
int i = 0;
while(i < n && Sturm[i].eval(x) == 0)
i++;
if (i == n)
return 0; // Tout les termes sont nuls
double terme_precedent = Sturm[i].eval(x), calc;
for(i = i+1 ; i < n ; i++){
calc = Sturm[i].eval(x);
if (terme_precedent * calc < 0){
nb++;
terme_precedent = calc;
}
}
return nb;
}
// Approxime la racine via l'algorithme de Newton
double Polynome::newton(double x) throw(std::domain_error, std::runtime_error){
Polynome DP = this->derive();
int iter = 0;
double f_x = eval(x), fp_x;
while(abs(f_x) > POLYNOME_EPSILON){
fp_x = DP.eval(x);
// Exceptions
if (abs(fp_x) < POLYNOME_EPSILON) // La dérivée est nulle
throw std::domain_error("Division par zéro");
if (iter > POLYNOME_NB_ITERATIONS_NEWTON){ // On a effectué trop d'itérations
if (abs(f_x) < POLYNOME_BIG_EPSILON) // Mais on est presque sur un zéro
return x;
throw std::runtime_error("Trop d'itérations");
}
x = x - f_x/fp_x;
f_x = eval(x);
iter++;
}
return x;
}
// Effectue un dichotomie sur un nombre d'itérations fixés
// Elle réduit l'intervalle [b_g, b_d] autour de la racine
void Polynome::dicho(double &b_g, double &b_d, int iter){
double val_g = eval(b_g), val_d = eval(b_d);
double milieu, val_m;
if (val_g * val_d > 0) // On ne peut pas lancer la dichotomie
return;
while (iter > 0){
milieu = (b_g + b_d) / 2;
val_m = eval(milieu);
if (val_g * val_m <= 0){ // La racine est située dans [b_g, milieu]
b_d = milieu; // On se place sur le segment [b_g, milieu]
}
else{ // La racine est située dans [milieu, b_d]
b_g = milieu; // On se place sur le segment [milieu, b_d]
val_g = val_m;
}
iter--;
}
}
/********************************************************
********************** Autres *************************
*******************************************************/
// Evalue le polynôme en un point, méthode de Horner
double Polynome::eval(double x) const{
double res = 0;
int n = this->degre();
for(int i = n ; i >= 0 ; i--){
res*= x;
res+= m_coeffs[i];
}
return res;
}
// Renvoie le degré du polynôme. Par convention le degré du polynôme nul sera -1 et non pas - l'infini
int Polynome::degre() const{
int n = m_coeffs.size() - 1;
while(n >= 0 && m_coeffs[n] == 0)
n--;
return n;
}
// Envoie l'expression litérale du polynôme sur out (utiliser par cout << P)
void Polynome::affiche(ostream &out) const{
int i;
double c;
bool poly_nul = true;
for(i=m_coeffs.size()-1 ; i!=-1 ; i--){
c = m_coeffs[i];
if (i == 0){
if(poly_nul)
out << noshowpos << c;
else
if(c != 0)
out << showpos << c;
}
else{
if(poly_nul && c != 0){
if(c == -1)
out << '-';
else if(c != 1)
out << noshowpos << c;
poly_nul = false;
}
else if(c != 0){
if( c == 1)
out << '+';
else if(c == -1)
out << '-';
else
out << showpos << c;
}
}
if(c != 0){
if(i==1)
out << 'X';
else if(i > 1)
out << "X^(" << noshowpos << i <<')';
}
}
out << noshowpos;
}
et un petit main pour tester ça
#include <iostream>
#include "polynome.h"
#include <windows.h>
using namespace std;
// Affiche un vector<double>
std::ostream &operator<<(std::ostream &out, const std::vector<double> &vec){
int n = vec.size();
for(int i = 0 ; i < n ; i++)
out << "[" << i << "] = " << vec[i] << endl;
return out;
}
int main(int argc, char *argv[])
{
// SetConsoleOutputCP( 1252 );
Polynome A(6, (double[]) {1, 0, 0, -1, 4., 1.});
cout << "A(X) = " << A << endl;
Polynome B(3, 1., 0., 1.); // Création via une liste variable de double
cout << "B(X) = " << B << endl;
cout << endl << "Division euclidienne : A = B * Q + R" << endl;
Polynome Q, R;
Polynome::div_eucli(A,B,Q,R);
cout << "Q(X) = " << Q << endl;
cout << "R(X) = " << R << endl;
cout << "Vérification (on teste en même temps l'opérateur d'égalité) : ";
if (A == B * Q + R)
cout << "ok" << endl;
else
cout << "erreur" << endl;
cout << endl << "A(X) / B(X) = " << A / B << endl;
cout << "A(X) % B(X) = " << A % B << endl;
cout << endl << "Calcul de PGCD" << endl;
double coeffs_A[] = {7, -7, 0, 8, -7, 0, 1};
A = Polynome(7, coeffs_A);
B = Polynome(6, (double[]) {-7, 0, 3, -7, 0, 3});
cout << "A(X) = " << A << endl << "B(X) = " << B << endl;
R = Polynome::pgcd(B,A);
cout << "PGCD(A,B) : " << R << endl;
R.unitaire();
cout << "PGCD(A,B) unitarisé : " << R << endl;
A = Polynome(6, -12., 20., -1., -9., 1., 1.);
//A = Polynome(3, -1., 0., 1.);
A = Polynome(11, 30., 30., -11., 35., -13., 36., 36., 20., 42., 15., -34.);
cout << endl << "A(X) = " << A << endl;
cout << "Nombre de racines réelles : " << A.nb_roots() << endl;
cout << "Racines : " << endl << A.roots();
return 0;
}
ce qui donne au final :
A(X) = X^(5)+4X^(4)-X^(3)+1
B(X) = X^(2)+1
Division euclidienne : A = B * Q + R
Q(X) = X^(3)+4X^(2)-2X-4
R(X) = 2X+5
Vérification (on teste en même temps l'opérateur d'égalité) : ok
A(X) / B(X) = X^(3)+4X^(2)-2X-4
A(X) % B(X) = 2X+5
Calcul de PGCD
A(X) = X^(6)-7X^(4)+8X^(3)-7X+7
B(X) = 3X^(5)-7X^(3)+3X^(2)-7
PGCD(A,B) : -0.25X^(3)-0.25
PGCD(A,B) unitarisé : X^(3)+1
A(X) = -34X^(10)+15X^(9)+42X^(8)+20X^(7)+36X^(6)+36X^(5)-13X^(4)+35X^(3)-11X^(2)
+30X+30
Nombre de racines réelles : 2
Racines :
[0] = -0.574217
[1] = 1.71694
Process returned 0 (0x0) execution time : 0.044 s
Press any key to continue.
Rapidement,
a- ton += n'est pas bon, tu as utilisé reserve au lieu de resize (c'est reserve+push_back, ou resize + op[])
b- utilise la symétrie, tu as tout à y gagner
c- pas de iostream dans un .h (iosfwd si tu passes des stream en paramètre, ostream/istream pour implémenter des opérateurs sur les flux, et iostream pour cout/cin/clog/cerr -- en tirant moins de choses, on y gagne sur les temps de compilation)
d- les spécifications d'exception sont des attrape-couillons qu'il veut mieux bannir de nos codes
e- pas besoin de stdarg dans le .h non plus.
f- il n'y a pas besoin de destructeur, ni de constructeur de copie, ni d'opérateur d'affectation.
g- les opérateurs x= devraient renvoyer des références et non des copies de *this.
h- pour recevoir le tableau de doubles, cf mon constructeur ligne 104 qui permet de recevoir directement un tableau _statique_ de doubles.
i- Pour les racines, j'aurais géré ça de manière complètement externe : les racines ne font pas parti du cœur du polynôme -- tu aurais pû profiter de l'info pour la multiplication, voire même la division, mais cela n'a pas été le cas.
h- donne des noms à tes paramètres dans les signatures des fonctions que tu déclares. En particulier pour la division (
Les problèmes de symétrie ont à plusieurs reprises été évoqués ici (ce fil) et ailleurs (le forum, et encore plus ailleurs). Je te laisse fouiller.
Pour les spécifications d'exception, je te laisse lire cette entrée du GOTW: http://www.gotw.ca/gotw/082.htm ...
Hum j'espérai que cela soit plus verbeux que cela.
Rapidement une spec d'exception n'est pas vérifiée à la compilation, mais à l'exécution. L'effet est que si une exception du mauvais type cherche à sortir, alors le programme sera arrêté. Oh on peut intervenir et faire 2-3 trucs au moment de l'arrêt, mais arrêt il y aura.
Et ce n'est pas du tout un comportement acceptable -- et en plus cela ralentit l'exécution vu que des vérifications sont réalisées.
Pas que je sache.
Enfin, il y a toujours ces histoires de conversions implicites qui doivent pouvoir s'appliquer.
Genre s'il existe un type Toto dépourvu du moins unaire, et un constructeur implicite qui prend un Toto. Je m'avance peut-être. C'est à tester.
Pair<A,B> est une structure contenant deux valeurs de types différents A et B. A et B peuvent aussi être le même type. Dans notre cas, c'est long à écrire, je crée donc un typedef pour simplifier la notation. Le typedef sert à créer un alias sur un nom de type.
La ligne qui te pose problème définit donc un alias nommé "Paire" qui remplacera donc les std::pair<Polynome,Polynome>, autrement dit les couples composés de deux polynomes. Ce que l'on veut pour la division puisqu'il y a le quotient et le reste.