/*************************************************************************
 *
 *  $RCSfile: hommatrixtemplate.hxx,v $
 *
 *  $Revision: 1.1 $
 *
 *  last change: $Author: aw $ $Date: 2003-04-15 16:02:53 $
 *
 *  The Contents of this file are made available subject to the terms of
 *  either of the following licenses
 *
 *         - GNU Lesser General Public License Version 2.1
 *         - Sun Industry Standards Source License Version 1.1
 *
 *  Sun Microsystems Inc., October, 2000
 *
 *  GNU Lesser General Public License Version 2.1
 *  =============================================
 *  Copyright 2000 by Sun Microsystems, Inc.
 *  901 San Antonio Road, Palo Alto, CA 94303, USA
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License version 2.1, as published by the Free Software Foundation.
 *
 *  This library 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
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
 *  MA  02111-1307  USA
 *
 *
 *  Sun Industry Standards Source License Version 1.1
 *  =================================================
 *  The contents of this file are subject to the Sun Industry Standards
 *  Source License Version 1.1 (the "License"); You may not use this file
 *  except in compliance with the License. You may obtain a copy of the
 *  License at http://www.openoffice.org/license.html.
 *
 *  Software provided under this License is provided on an "AS IS" basis,
 *  WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING,
 *  WITHOUT LIMITATION, WARRANTIES THAT THE SOFTWARE IS FREE OF DEFECTS,
 *  MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE, OR NON-INFRINGING.
 *  See the License for the specific provisions governing your rights and
 *  obligations concerning the Software.
 *
 *  The Initial Developer of the Original Code is: Sun Microsystems, Inc.
 *
 *  Copyright: 2000 by Sun Microsystems, Inc.
 *
 *  All Rights Reserved.
 *
 *  Contributor(s): _______________________________________
 *
 *
 ************************************************************************/

#ifndef _HOMMATRIX_TEMPLATE_HXX
#define _HOMMATRIX_TEMPLATE_HXX

#ifndef _SAL_TYPES_H_
#include <sal/types.h>
#endif

#include <math.h>
#include <string.h>

inline double ImplGetDefaultValue(sal_uInt16 nRow, sal_uInt16 nColumn)
{
    if(nRow == nColumn)
        return 1.0;
    return 0.0;
}

template < int _RowSize > class ImplMatLine
{
    enum { RowSize = _RowSize };

    double                                          mfValue[RowSize];

public:
    ImplMatLine(sal_uInt16 nRow = 0L, ImplMatLine< RowSize >* pToBeCopied = 0L)
    {
        if(pToBeCopied)
        {
            memcpy(&mfValue, pToBeCopied, sizeof(double) * RowSize);
        }
        else
        {
            for(sal_uInt16 a(0); a < RowSize; a++)
            {
                mfValue[a] = ImplGetDefaultValue(nRow, a);
            }
        }
    }

    double Get(sal_uInt16 nColumn) const
    {
        return mfValue[nColumn];
    }

    void Set(sal_uInt16 nColumn, const double& rValue)
    {
        mfValue[nColumn] = rValue;
    }
};

template < int _RowSize > class ImplHomMatrixTemplate
{
    enum { RowSize = _RowSize };

    sal_uInt32                                      mnRefCount;

    ImplMatLine< RowSize >                          maLine[RowSize - 1];
    ImplMatLine< RowSize >*                         mpLine;

public:
    // Is last line used?
    sal_Bool IsLastLineDefault() const
    {
        if(!mpLine)
            return sal_True;

        for(sal_uInt16 a(0); a < RowSize; a++)
        {
            if(ImplGetDefaultValue((RowSize - 1), a) != mpLine->Get(a))
            {
                return sal_False;
            }
        }

        // reset last line, it equals default
        delete ((ImplHomMatrixTemplate< RowSize >*)this)->mpLine;
        ((ImplHomMatrixTemplate< RowSize >*)this)->mpLine = 0L;

        return sal_True;
    }

    // This constructor is only used form the static identity matrix, thus
    // the RefCount is set to 1 to never 'delete' this static incarnation.
    ImplHomMatrixTemplate()
    :   mnRefCount(1),
        mpLine(0L)
    {
        // complete initialization with identity matrix, all lines
        // were initialized with a trailing 1 followed by 0's.
        for(sal_uInt16 a(1); a < (RowSize - 1); a++)
        {
            maLine[a].Set(0, 0.0);
            maLine[a].Set(a, 1.0);
        }
    }

    ImplHomMatrixTemplate(const ImplHomMatrixTemplate& rToBeCopied)
    :   mnRefCount(0),
        mpLine(0L)
    {
        // complete initialization using copy
        for(sal_uInt16 a(0); a < (RowSize - 1); a++)
        {
            memcpy(&maLine[a], &rToBeCopied.maLine[a], sizeof(ImplMatLine< RowSize >));
        }

        if(rToBeCopied.mpLine)
        {
            mpLine = new ImplMatLine< RowSize >((RowSize - 1), rToBeCopied.mpLine);
        }
    }

    ~ImplHomMatrixTemplate()
    {
        if(mpLine)
        {
            delete mpLine;
        }
    }

    const sal_uInt32 GetRefCount() const { return mnRefCount; }
    void IncRefCount() { mnRefCount++; }
    void DecRefCount() { mnRefCount--; }

    sal_uInt16 GetEdgeLength() const { return RowSize; }

    double Get(sal_uInt16 nRow, sal_uInt16 nColumn) const
    {
        if(nRow < (RowSize - 1))
            return maLine[nRow].Get(nColumn);
        if(mpLine)
            return mpLine->Get(nColumn);
        return ImplGetDefaultValue((RowSize - 1), nColumn);
    }

    void Set(sal_uInt16 nRow, sal_uInt16 nColumn, const double& rValue)
    {
        if(nRow < (RowSize - 1))
        {
            maLine[nRow].Set(nColumn, rValue);
        }
        else if(mpLine)
        {
            mpLine->Set(nColumn, rValue);
        }
        else
        {
            if(ImplGetDefaultValue((RowSize - 1), nColumn) != rValue)
            {
                mpLine = new ImplMatLine< RowSize >((RowSize - 1), 0L);
                mpLine->Set(nColumn, rValue);
            }
        }
    }

    void TestLastLine()
    {
        if(mpLine)
        {
            sal_Bool bNecessary(sal_False);

            for(sal_uInt16 a(0);!bNecessary && a < RowSize; a++)
            {
                if(ImplGetDefaultValue((RowSize - 1), a) != mpLine->Get(a))
                {
                    bNecessary = sal_True;
                }
            }

            if(!bNecessary)
            {
                delete mpLine;
                mpLine = 0L;
            }
        }
    }

    // Left-upper decompositon
    sal_Bool Ludcmp(sal_uInt16 nIndex[], sal_Int16& nParity)
    {
        double fBig, fSum, fDum;
        double fStorage[RowSize];
        sal_uInt16 a, b, c, nAMax;

        nParity = 1;

        // Calc the max of each line. If a line is empty,
        // stop immediately since matrix is not invertible then.
        for(a = 0; a < RowSize; a++)
        {
            fBig = 0.0;

            for(b = 0; b < RowSize; b++)
            {
                double fTemp(fabs(Get(a, b)));

                if(fTemp > fBig)
                    fBig = fTemp;
            }

            if(0.0 == fBig)
                return sal_False;

            fStorage[a] = 1.0 / fBig;
        }

        // start normalizing
        for(b = 0; b < RowSize; b++)
        {
            for(a = 0; a < b; a++)
            {
                fSum = Get(a, b);

                for(c = 0; c < a; c++)
                {
                    fSum -= Get(a, c) * Get(c, b);
                }

                Set(a, b, fSum);
            }

            fBig = 0.0;

            for(a = b; a < RowSize; a++)
            {
                fSum = Get(a, b);

                for(c = 0; c < b; c++)
                {
                    fSum -= Get(a, c) * Get(c, b);
                }

                Set(a, b, fSum);
                fDum = fStorage[a] * fabs(fSum);

                if(fDum >= fBig)
                {
                    fBig = fDum;
                    nAMax = a;
                }
            }

            if(b != nAMax)
            {
                for(c = 0; c < RowSize; c++)
                {
                    fDum = Get(nAMax, c);
                    Set(nAMax, c, Get(b, c));
                    Set(b, c, fDum);
                }

                nParity = -nParity;
                fStorage[nAMax] = fStorage[b];
            }

            nIndex[b] = nAMax;

            // here the failure of precision occurs
            if(0.0 == fabs(Get(b, b)))
                return sal_False;

            if(b != (RowSize - 1))
            {
                fDum = 1.0 / Get(b, b);

                for(a = b + 1; a < RowSize; a++)
                {
                    Set(a, b, Get(a, b) * fDum);
                }
            }
        }

        return sal_True;
    }

    void Lubksb(const sal_uInt16 nIndex[], double fRow[]) const
    {
        sal_uInt16 b, ip;
        sal_Int16 a, a2 = -1;
        double fSum;

        for(a = 0; a < RowSize; a++)
        {
            ip = nIndex[a];
            fSum = fRow[ip];
            fRow[ip] = fRow[a];

            if(a2 >= 0)
            {
                for(b = a2; b < a; b++)
                {
                    fSum -= Get(a, b) * fRow[b];
                }
            }
            else if(0.0 != fSum)
            {
                a2 = a;
            }

            fRow[a] = fSum;
        }

        for(a = (RowSize - 1); a >= 0; a--)
        {
            fSum = fRow[a];

            for(b = a + 1; b < RowSize; b++)
            {
                fSum -= Get(a, b) * fRow[b];
            }

            if(0.0 != Get(a, a))
            {
                fRow[a] = fSum / Get(a, a);
            }
        }
    }

    sal_Bool IsIdentity() const
    {
        // last line needs no testing if not existing
        const sal_uInt16 nMaxLine = (mpLine) ? RowSize : (RowSize - 1);

        for(sal_uInt16 a(0); a < nMaxLine; a++)
        {
            for(sal_uInt16 b(0); b < RowSize; b++)
            {
                if(Get(a, b) != ImplGetDefaultValue(a, b))
                    return sal_False;
            }
        }

        return sal_True;
    }

    sal_Bool IsInvertible() const
    {
        ImplHomMatrixTemplate aWork(*this);
        sal_uInt16 nIndex[RowSize];
        sal_Int16 nParity;

        return aWork.Ludcmp(nIndex, nParity);
    }

    sal_Bool IsNormalized() const
    {
        if(!mpLine)
            return sal_True;

        const double fHomValue(Get((RowSize - 1), (RowSize - 1)));

        if(1.0 != fHomValue && 0.0 != fHomValue)
            return sal_False;

        return sal_True;
    }

    void DoInvert(const ImplHomMatrixTemplate& rWork, const sal_uInt16 nIndex[])
    {
        double fArray[RowSize];

        for(sal_uInt16 a(0); a < RowSize; a++)
        {
            // prepare line
            for(sal_uInt16 b(0); b < RowSize; b++)
                fArray[b] = ImplGetDefaultValue(a, b);

            // expand line
            rWork.Lubksb(nIndex, fArray);

            // copy line transposed to this matrix
            for(b = 0; b < RowSize; b++)
                Set(b, a, fArray[b]);
        }

        // evtl. get rid of last matrix line
        TestLastLine();
    }

    void DoNormalize()
    {
        if(mpLine)
        {
            const double fHomValue(Get((RowSize - 1), (RowSize - 1)));

            for(sal_uInt16 a(0); a < RowSize; a++)
            {
                for(sal_uInt16 b(0); b < RowSize; b++)
                {
                    Set(a, b, Get(a, b) / fHomValue);
                }
            }

            // evtl. get rid of last matrix line
            TestLastLine();
        }
    }

    double DoDeterminant() const
    {
        ImplHomMatrixTemplate aWork(*this);
        sal_uInt16 nIndex[RowSize];
        sal_Int16 nParity;
        double fRetval(0.0);

        if(aWork.Ludcmp(nIndex, nParity))
        {
            fRetval = (double)nParity;

            // last line needs no multiply if not existing; default value would be 1.
            const sal_uInt16 nMaxLine = (aWork.mpLine) ? RowSize : (RowSize - 1);

            for(sal_uInt16 a(0); a < nMaxLine; a++)
            {
                fRetval *= aWork.Get(a, a);
            }
        }

        return fRetval;
    }

    double DoTrace() const
    {
        double fTrace = (mpLine) ? 0.0 : 1.0;
        const sal_uInt16 nMaxLine = (mpLine) ? RowSize : (RowSize - 1);

        for(sal_uInt16 a(0); a < nMaxLine; a++)
        {
            fTrace += Get(a, a);
        }

        return fTrace;
    }

    void DoTranspose()
    {
        for(sal_uInt16 a(0); a < (RowSize - 1); a++)
        {
            for(sal_uInt16 b(a + 1); b < RowSize; b++)
            {
                const double fTemp(Get(a, b));
                Set(a, b, Get(b, a));
                Set(b, a, fTemp);
            }
        }

        TestLastLine();
    }

    void DoAddMatrix(const ImplHomMatrixTemplate& rMat)
    {
        for(sal_uInt16 a(0); a < RowSize; a++)
        {
            for(sal_uInt16 b(0); b < RowSize; b++)
            {
                Set(a, b, Get(a, b) + rMat.Get(a, b));
            }
        }

        TestLastLine();
    }

    void DoSubMatrix(const ImplHomMatrixTemplate& rMat)
    {
        for(sal_uInt16 a(0); a < RowSize; a++)
        {
            for(sal_uInt16 b(0); b < RowSize; b++)
            {
                Set(a, b, Get(a, b) - rMat.Get(a, b));
            }
        }

        TestLastLine();
    }

    void DoMulMatrix(const double& rfValue)
    {
        for(sal_uInt16 a(0); a < RowSize; a++)
        {
            for(sal_uInt16 b(0); b < RowSize; b++)
            {
                Set(a, b, Get(a, b) * rfValue);
            }
        }

        TestLastLine();
    }

    void DoMulMatrix(const ImplHomMatrixTemplate& rMat)
    {
        // create a copy as source for the original values
        const ImplHomMatrixTemplate aCopy(*this);

        if(rMat.mpLine)
        {
            // complex first lines
            for(sal_uInt16 a(0); a < (RowSize - 1); a++)
            {
                for(sal_uInt16 b(0); b < RowSize; b++)
                {
                    double fValue(0.0);

                    for(sal_uInt16 c(0); c < RowSize; c++)
                    {
                        fValue += aCopy.Get(c, b) * rMat.Get(a, c);
                    }

                    Set(a, b, fValue);
                }
            }

            if(mpLine)
            {
                // complex last line
                for(sal_uInt16 b(0); b < RowSize; b++)
                {
                    double fValue(0.0);

                    for(sal_uInt16 c(0); c < RowSize; c++)
                    {
                        fValue += aCopy.Get(c, b) * rMat.Get((RowSize - 1), c);
                    }

                    Set((RowSize - 1), b, fValue);
                }

                TestLastLine();
            }
            else
            {
                // last line is rMat line, copy
                mpLine = new ImplMatLine< RowSize >((RowSize - 1), rMat.mpLine);
            }
        }
        else
        {
            // simplified first lines
            for(sal_uInt16 a(0); a < (RowSize - 1); a++)
            {
                for(sal_uInt16 b(0); b < RowSize; b++)
                {
                    double fValue(0.0);

                    for(sal_uInt16 c(0); c < (RowSize - 1); c++)
                    {
                        fValue += aCopy.Get(c, b) * rMat.Get(a, c);
                    }

                    if(b == (RowSize -1))
                    {
                        fValue += aCopy.Get(a, (RowSize - 1));
                    }

                    Set(a, b, fValue);
                }
            }

            if(mpLine)
            {
                // simplified last line
                for(sal_uInt16 b(0); b < RowSize; b++)
                {
                    double fValue(0.0);

                    for(sal_uInt16 c(0); c < (RowSize - 1); c++)
                    {
                        fValue += aCopy.Get(c, b) * rMat.Get((RowSize - 1), c);
                    }

                    if(b == (RowSize -1))
                    {
                        fValue += aCopy.Get((RowSize - 1), (RowSize - 1));
                    }

                    Set((RowSize - 1), b, fValue);
                }

                TestLastLine();
            }
            else
            {
                // no last line
            }
        }
    }

    sal_Bool IsEqual(const ImplHomMatrixTemplate& rMat)
    {
        const sal_uInt16 nMaxLine = (mpLine || rMat.mpLine) ? RowSize : (RowSize - 1);

        for(sal_uInt16 a(0); a < nMaxLine; a++)
        {
            for(sal_uInt16 b(0); b < RowSize; b++)
            {
                if(Get(a, b) != rMat.Get(a, b))
                    return sal_False;
            }
        }

        return sal_True;
    }
};

#endif // _HOMMATRIX_TEMPLATE_HXX