La librairie que je suis en train de concevoir pour exécuter, entre autres choses, des calculs sur des polynômes, des vecteurs et des matrices.
#include <stdarg.h>
#include <math.h>
using namespace std;
class QiXMath : public std::exception
{
public:
QiXMath(char* msg) { _whatMsg = msg; }
virtual const char* what() const throw() { return _whatMsg; }
protected:
char* _whatMsg;
};
//Pour représenter des fractions, ce qui est très utile pour les polynôme quand on veut des calculs exacts.
//Note : pour simplifier des nombres décimaux qui ont des séquences répétitives voir http://en.wikipedia.org/wiki/QiFraction_(mathematics)
//TOCHECK : l'opérateur /(int) ne doit pas être nécéssaire puisque j'ai un constructeur qui a un int en entrée, pourquoi dois-je en avoir un ???
template <class T>
class QiFraction
{
public:
//Faire un constructeur avec un seul paramètre : le numérateur.
QiFraction(T n, T d) { _n = n; _d = d; }
QiFraction(T n=0) { _n = n; _d = 1; }
QiFraction(float n) { _n = (T) n; _d = 1; }
QiFraction(int n) { _n = (T) n; _d = 1; }
void dump() const
{
if (_n == 0 || _d == 1) cout << _n;
else {
if (_n == _d) cout << "1";
else cout << _n << "/" << _d;
}
}
QiFraction<T>& operator+=(const QiFraction<T>& a)
{
T P = pgcd(_d, a._d);
T NewD = (_d * a._d) / P;
_n = _n * (NewD / _d) + a._n * (NewD / a._d);
_d = NewD;
return *this;
}
QiFraction<T>& operator-=(const QiFraction<T>& a)
{
T P = pgcd(_d, a._d);
T NewD = (_d * a._d) / P;
_n = _n * (NewD / _d) - a._n * (NewD / a._d);
_d = NewD;
return *this;
}
QiFraction<T>& operator*=(const QiFraction<T>& a)
{
_n *= a._n;
_d *= a._d;
//On simplifie la fraction
T Common = pgcd(_n, _d);
_n /= Common;
_d /= Common;
return *this;
}
QiFraction<T>& operator/=(const QiFraction<T>& a)
{
_n *= a._d;
_d *= a._n;
//On simplifie la fraction
T Common = pgcd(_n, _d);
_n /= Common;
_d /= Common;
return *this;
}
QiFraction<T> operator+(const QiFraction<T>& a) const
{
QiFraction Result(_n, _d);
return Result += a;
}
QiFraction<T> operator-(const QiFraction<T>& a) const
{
QiFraction Result(_n, _d);
return Result -= a;
}
QiFraction<T> operator*(const QiFraction<T>& a) const
{
QiFraction Result(_n, _d);
return Result *= a;
}
QiFraction<T> operator/(const QiFraction<T>& a) const
{
QiFraction Result(_n, _d);
return Result /= a;
}
QiFraction<T> operator/(T a) const
{
QiFraction<T> Result(_n, _d);
return Result /= QiFraction<T>(a);
}
QiFraction<T> operator/(int a) const
{
QiFraction<T> Result(_n, _d);
return Result /= QiFraction<T>(a);
}
bool operator>=(float t) const { return (float)_n/(float)_d > t; }
bool operator==(float t) const { return (float)_n/(float)_d == t; }
operator float() const { return (float)_n/(float)_d; }
QiFraction<T>& inverse()
{
T Tmp = _n;
_n = _d;
_d = _n;
return *this;
}
T pgcd(T u, T v/*, int verbose=false*/)
{
T M = 1;
if (u<0) u = (T)0-u;
if (v<0) v = (T)0-v;
while (u && v) {
// if (verbose) printf("pgcd(%d, %d)\n", u, v);
if (!(u&1) && !(v&1)) { M <<= 1; u>>=1, v>>=1; }
else if (!(u&1) && v&1) u>>=1;
else if (u&1 && !(v&1)) v>>=1;
else { //Les deux nombres sont impairs
if (u>v) u=(u-v)>>1;
else v=(v-u)>>1;
}
}
return M * (u ? u : v);
}
//Note : ppcm(a,b) = (a*b)/pgcd(a,b)
T _n, _d; //Numérateur et dénominateur
};
//###TODO (25%) : Exceptions
//###TODO : cacher l'évaluation du degré d'un polynôme
//###TODO : faire la factorisation simple d'un polynôme
//###TODO : faire des opérateurs pour QiMatrix et QiPolynomial
//###TODO : faire des opérateurs pour QiVector
template <class T>
class QiPolynomial
{
public:
QiPolynomial(int rank)
{
_comps = NULL;
_init(rank);
}
~QiPolynomial()
{
DELETEARRAY(_comps);
}
QiPolynomial(const QiPolynomial<T>& original)
{
if (this != &original) {
_init(original.rank());
int i;
for(i=0; i<=rank(); i++) _comps[i] = original.get(i);
}
}
protected:
void _init(int rank)
{
if (rank<0) return;
_rank = rank;
_comps = new T[_rank+1];
}
void _expand(int nbWantedSlots)
{
int NewSize = nbWantedSlots;
_comps = new T[NewSize];
if (!_comps) throw new QiXMath("QiPolynomial::_expand : NotImplemented");
// memcpy();
}
public:
//
void null()
{
int i;
for(i=0; i<=rank(); i++) _comps[i] = 0.0f;
}
void primitive()
{
int Degree = degree();
if (rank() < Degree+1) _expand(Degree+1);
int i, MaxIndex = degree();
for(i=MaxIndex; i>=0; i--) _comps[i+1] = _comps[i]/(i+1);
_comps[0] = 0.0f;
}
void add(const QiPolynomial<T>& op)
{
int i;
int MaxIndex = op.rank();
if (rank() < MaxIndex) MaxIndex = rank();
for(i=0; i<=MaxIndex; i++) set(i, get(i) + op.get(i));
}
QiPolynomial<T>& operator *= (const QiPolynomial<T>& a)
{
multiply(a);
return *this;
}
bool multiply(const QiPolynomial<T>& op)
{
if (degree() + op.degree() > rank()) return false;
QiPolynomial Temp(degree() + op.degree());
int MaxIt = op.degree(), i, j;
int CurrentDegree = degree();
Temp.null();
for(i=0; i<=MaxIt; i++) {
for(j=0; j<=CurrentDegree; j++) {
T CurrentValue = Temp.get(i+j);
T NewInc = op.get(i) * get(j);
Temp.set(i+j, CurrentValue + NewInc);
}
}
set(Temp);
return true;
}
void multiply(T coeff)
{
int i;
int MaxIndex = degree();
for(i=0; i<=MaxIndex; i++) set(i, get(i)*coeff);
}
QiPolynomial<T>& operator /= (const T& a)
{
int i;
int MaxIndex = degree();
for(i=0; i<=MaxIndex; i++) {
set(i, get(i) / a);
}
return *this;
}
void set(const QiPolynomial& src)
{
int i;
null();
int MaxIndex = rank();
if (MaxIndex>src.rank()) MaxIndex = src.rank();
for(i=0; i<=MaxIndex; i++) set(i, src.get(i));
// _rank = src.rank();
}
void setRaw(int maxDegree, ...)
{
if (maxDegree > rank()) return;
va_list args;
va_start(args, maxDegree);
int i, j;
for(i=maxDegree; i>=0; i--) {
float Value = (float) va_arg(args, double);
_comps[i] = Value;
}
va_end(args);
}
//On remplace la variable x du polynôme par un autre polynôme (x par x+1 par exemple).
//On substitue donc Q(x) à x dans P(x), ce qui donne P(Q(x)), un autre polynôme en x.
//La méthode :
void substitute(const QiPolynomial& newVar)
{
QiPolynomial Value(newVar.degree() + degree()), Total(newVar.degree() + degree());
QiPolynomial Param(newVar.degree() + degree());
int i;
Param.null();
Param.set(0, 1.0f);
for(i=0; i<=degree(); i++) {
Value.set(Param);
Value.multiply(get(i));
Total.add(Value);
Param.multiply(newVar);
}
set(Total);
}
T get(int index) const
{
if (index<0 || index>rank()+1) throw new QiXMath("QiPolynomial::get out of bound");
return _comps[index];
}
void set(int index, T value)
{
if (index<0 || index>rank()+1) throw new QiXMath("QiPolynomial::set out of bound");
_comps[index] = value;
}
float evaluate(T x) const
{
int MaxIndex = degree();
float Current = 1.0f, Total = 0.0f;
int i;
for(i=0; i<=MaxIndex; i++) {
Total += Current * _comps[i];
Current *= x;
}
return Total;
}
int degree() const //Renvoie le degré du polynôme courant : l'indice du premier coefficient non-nul en partant des puissances élevées.
{
int i = rank();
while (get(i) == 0.0f && i>0) i--;
return i;
}
void dump() const
{
int i, Degree = degree();
for(i=Degree; i>=0; i--) {
if (i<Degree && get(i)>=0.0f) printf("+");
printf("%.10f", (float) get(i));
if (!i) continue;
if (i==1) printf("X");
if (i>1) printf("X^%d", i);
}
printf("\n");
}
void dumpObjects() const
{
int i, Degree = degree();
for(i=Degree; i>=0; i--) {
T CurValue = get(i);
if (i<Degree && CurValue>=0.0f) printf("+");
if (CurValue != 0.0f) {
get(i).dump();
if (!i) continue;
if (i==1) printf("X");
if (i>1) printf("X^%d", i);
}
}
printf("\n");
}
int rank() const { return _rank; }
protected:
T* _comps;
int _rank;
};
template <class T>
class QiVector
{
public:
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class ExQiVector : public std::exception
{
public:
ExQiVector(char*) { }
};
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
QiVector(int rank) { _comps = new T[rank]; _rank = rank;}
~QiVector() { DELETEARRAY(_comps) }
void null() {
int i;
for(i=0; i<rank(); i++) _comps[i] = 0.0f;
}
//On crée un vecteur de base
void base(int index)
{
if (index<0 || index>=rank()) return;
null();
set(index, 1.0f);
}
//###TODO : faire des exceptions pour les out of bounds.
bool set(int index, const T& value) {
if (index<0 || index >= rank()) return false;
_comps[index] = value;
return true;
}
bool set(QiVector<T>& newValue)
{
if (rank() != newValue.rank()) return false;
int i;
for(i=0; i<rank(); i++) _comps[i] = newValue.get(i);
return true;
}
bool add(int index, const T& value)
{
if (index<0 || index>=rank()) return false;
return set(index, get(index) + value);
}
T get(int index) const
{
return _comps[index];
}
QiVector<T>& operator=(const QiVector<T>& original)
{
//###TODO : si les rangs sont différents on doit changer le rang de l'objet courant
if (this != &original && rank() == original.rank()) {
int i;
for(i=0; i<rank(); i++) _comps[i] = original.get(i);
}
return *this;
}
int rank() const { return _rank; }
void dump() const
{
int i;
printf("(");
for(i=0; i<rank(); i++) {
if(i) printf(", ");
printf("%f", _comps[i]);
}
printf(")\n");
}
void dumpObjects() const
{
int i;
printf("(");
for(i=0; i<rank(); i++) {
if(i) printf(", ");
_comps[i].dump();
}
printf(")\n");
}
protected:
T* _comps;
int _rank;
};
//Matrice carrée constituée de vecteurs QiVector
template <class T>
class QiMatrix
{
public:
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class ExQiMatrix : public std::exception
{
public:
ExQiMatrix(char*) { }
};
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
QiMatrix(int rank) { _init(rank); }
protected:
bool _init(int rank)
{
_comps = new QiVector<T>*[rank];
int i;
for(i=0; i<rank; i++) _comps[i] = new QiVector<T>(rank);
_rank = rank;
return true;
}
bool _release()
{
int i;
for(i=0; i<rank(); i++) {
DELETESINGLE(_comps[i]);
}
DELETEARRAY(_comps);
}
void _count() const
{
_nbUps = 0, _nbLows = 0, _nbDiags = 0;
int i, j;
for(i=0; i<rank(); i++) {
for(j=0; j<rank(); j++) {
float Value = get(i, j);
if (i>j && Value) _nbLows++;
if (i<j && Value) _nbUps++;
if (i==j && Value) _nbDiags++;
}
}
//
_nbItemsInHalf = _rank*(_rank-1)/2 + _rank;
}
public:
~QiMatrix() { _release(); }
void null() { int i; for(i=0; i<rank(); i++) _comps[i]->null(); }
void identity()
{
null();
int i;
for(i=0; i<rank(); i++) _comps[i]->set(i, 1.0f);
}
T get(int i, int j) const
{
if (i<0 || j<0) return 0.0f;
if (i>=rank() || j>=rank()) return 0.0f;
return _comps[i]->get(j);
}
QiVector<T>& get(int index) const
{
if (index<0 || index>=rank()) return *_comps[0];
return *_comps[index];
}
bool set(int i, int j, T value)
{
if (j<0 || j>=rank()) return false;
_comps[i]->set(j, value);
return true;
}
bool setRaw(int setRank, ...) {
if (setRank != rank()) return false;
va_list args;
va_start(args, setRank);
int i, j;
for(i=0; i<rank(); i++) {
for(j=0; j<rank(); j++) {
if (!_comps[i]->set(j, (float) va_arg(args, double))) {
va_end(args);
return false;
}
}
}
va_end(args);
return true;
}
bool setRow(int index, QiVector<T>& row)
{
if (rank() != row.rank()) return false;
_comps[index]->set(row);
return true;
}
//Multiplication matrice-matrice inplace
QiMatrix<T>& operator *= (const QiMatrix<T>& op)
{
multiply(op);
return *this;
}
bool multiply(const QiMatrix<T>& op)
{
if (rank() != op.rank()) return false;
QiVector<T> Temp(rank());
int i, j, k;
for(i=0; i<rank(); i++) {
for(Temp.null(), j=0; j<rank(); j++) {
for(k=0; k<rank(); k++) {
Temp.add(j, get(i, k) * op.get(k, j));
}
}
setRow(i, Temp);
}
return true;
}
//
QiVector<T> operator * (const QiVector<T>& op) const
{
QiVector<T> Result(op.rank());
multiply(op, Result);
return Result;
}
void multiply(const QiVector<T>& op, QiVector<T>& result) const
{
int i, j;
for(i=0; i<rank(); i++) {
T Total = 0.0f;
for(j=0; j<rank(); j++) {
Total += get(i, j) * op.get(j);
}
result.set(i, Total);
}
}
//Copy constructor
QiMatrix(const QiMatrix<T>& original)
{
_init(original.rank());
int i;
for(i=0; i<rank(); i++) *_comps[i] = original.get(i);
}
bool isLower() const
{
_count();
return _nbLows + _nbDiags == _nbItemsInHalf;
}
bool isDiagonal() const
{
_count();
return _nbDiags == rank();
}
//Inverse la matrice courante "inplace"
void inverse()
{
if (isLower()) {
//Cas pour les matrices triangulaires inférieures
QiVector<T> Temp(rank());
int i, j, k;
T Total;
for(i=0; i<rank(); i++) {
Temp.null();
for(j=0; j<=i; j++) {
Total = 0.0f;
for(k=j; k<i; k++) {
T NormalIK = get(i, k);
T InvKJ = get(k, j);
Total += NormalIK * InvKJ;
}
T DeltaIJ = i == j ? 1.0f : 0.0f;
T MainCoeff = get(i, i);
T InverseCoeff = (DeltaIJ - Total) / MainCoeff;
Temp.set(j, InverseCoeff);
}
setRow(i, Temp);
}
} else if(isDiagonal()) {
int i;
for(i=0; i<rank(); i++) set(i, i, 1.0f/get(i, i));
} else throw new QiXMath("QiMatrix::inverse : not yet implemented.");
}
void dump() const
{
int i;
for(i=0; i<rank(); i++) _comps[i]->dump();
}
void dumpObjects() const
{
int i;
for(i=0; i<rank(); i++) _comps[i]->dumpObjects();
}
int rank() const { return _rank; }
//
bool isUpper() const;
bool isInversible() const;
protected:
QiVector<T>** _comps;
int _rank;
//
mutable int _nbUps, _nbLows, _nbDiags, _nbItemsInHalf;
};
Bien se souvenir qu'ici je cherche plus à être exact que rapide, c'est sûr que pour un moteur 3D j'utiliserai les classes QiMatrix et QiVector avec beaucoup de pincettes !!!
