Librairie de Math

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