Quelle  interpr5.cxx   Sprache: C

/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
/*
 * This file is part of the LibreOffice project.
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
 *
 * This file incorporates work covered by the following license notice:
 *
 *   Licensed to the Apache Software Foundation (ASF) under one or more
 *   contributor license agreements. See the NOTICE file distributed
 *   with this work for additional information regarding copyright
 *   ownership. The ASF licenses this file to you under the Apache
 *   License, Version 2.0 (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.apache.org/licenses/LICENSE-2.0 .
 */

#include <rtl/math.hxx>
#include <string.h>
#include <math.h>

#ifdef DEBUG_SC_LUP_DECOMPOSITION
#include <stdio.h>
#endif

#include <unotools/bootstrap.hxx>
#include <svl/zforlist.hxx>
#include <tools/duration.hxx>

#include <interpre.hxx>
#include <global.hxx>
#include <formulacell.hxx>
#include <document.hxx>
#include <dociter.hxx>
#include <scmatrix.hxx>
#include <globstr.hrc>
#include <scresid.hxx>
#include <cellkeytranslator.hxx>
#include <formulagroup.hxx>
#include <vcl/svapp.hxx> //Application::

#include <vector>

using ::std::vector;
using namespace formula;

namespace {

double MatrixAdd(const double& lhs, const double& rhs)
{
    return ::rtl::math::approxAdd( lhs,rhs);
}

double MatrixSub(const double& lhs, const double& rhs)
{
    return ::rtl::math::approxSub( lhs,rhs);
}

double MatrixMul(const double& lhs, const double& rhs)
{
    return lhs * rhs;
}

double MatrixDiv(const double& lhs, const double& rhs)
{
    return ScInterpreter::div( lhs,rhs);
}

double MatrixPow(const double& lhs, const double& rhs)
{
    return ::pow( lhs,rhs);
}

// Multiply n x m Mat A with m x l Mat B to n x l Mat R
void lcl_MFastMult(const ScMatrixRef& pA, const ScMatrixRef& pB, const ScMatrixRef& pR,
                   SCSIZE n, SCSIZE m, SCSIZE l)
{
    for (SCSIZE row = 0; row < n; row++)
    {
        for (SCSIZE col = 0; col < l; col++)
        {   // result element(col, row) =sum[ (row of A) * (column of B)]
            KahanSum fSum = 0.0;
            for (SCSIZE k = 0; k < m; k++)
                fSum += pA->GetDouble(k,row) * pB->GetDouble(col,k);
            pR->PutDouble(fSum.get(), col, row);
        }
    }
}

}

double ScInterpreter::ScGetGCD(double fx, double fy)
{
    // By ODFF definition GCD(0,a) => a. This is also vital for the code in
    // ScGCD() to work correctly with a preset fy=0.0
    if (fy == 0.0)
        return fx;
    else if (fx == 0.0)
        return fy;
    else
    {
        double fz = fmod(fx, fy);
        while (fz > 0.0)
        {
            fx = fy;
            fy = fz;
            fz = fmod(fx, fy);
        }
        return fy;
    }
}

void ScInterpreter::ScGCD()
{
    short nParamCount = GetByte();
    if ( !MustHaveParamCountMin( nParamCount, 1 ) )
        return;

    double fx, fy = 0.0;
    ScRange aRange;
    size_t nRefInList = 0;
    while (nGlobalError == FormulaError::NONE && nParamCount-- > 0)
    {
        switch (GetStackType())
        {
            case svDouble :
            case svString:
            case svSingleRef:
            {
                fx = ::rtl::math::approxFloor( GetDouble());
                if (fx < 0.0)
                {
                    PushIllegalArgument();
                    return;
                }
                fy = ScGetGCD(fx, fy);
            }
            break;
            case svDoubleRef :
            case svRefList :
            {
                FormulaError nErr = FormulaError::NONE;
                PopDoubleRef( aRange, nParamCount, nRefInList);
                double nCellVal;
                ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
                if (aValIter.GetFirst(nCellVal, nErr))
                {
                    do
                    {
                        fx = ::rtl::math::approxFloor( nCellVal);
                        if (fx < 0.0)
                        {
                            PushIllegalArgument();
                            return;
                        }
                        fy = ScGetGCD(fx, fy);
                    } while (nErr == FormulaError::NONE && aValIter.GetNext(nCellVal, nErr));
                }
                SetError(nErr);
            }
            break;
            case svMatrix :
            case svExternalSingleRef:
            case svExternalDoubleRef:
            {
                ScMatrixRef pMat = GetMatrix();
                if (pMat)
                {
                    SCSIZE nC, nR;
                    pMat->GetDimensions(nC, nR);
                    if (nC == 0 || nR == 0)
                        SetError(FormulaError::IllegalArgument);
                    else
                    {
                     double nVal = pMat->GetGcd();
                     fy = ScGetGCD(nVal,fy);
                    }
                }
            }
            break;
            default : SetError(FormulaError::IllegalParameter); break;
        }
    }
    PushDouble(fy);
}

void ScInterpreter:: ScLCM()
{
    short nParamCount = GetByte();
    if ( !MustHaveParamCountMin( nParamCount, 1 ) )
        return;

    double fx, fy = 1.0;
    ScRange aRange;
    size_t nRefInList = 0;
    while (nGlobalError == FormulaError::NONE && nParamCount-- > 0)
    {
        switch (GetStackType())
        {
            case svDouble :
            case svString:
            case svSingleRef:
            {
                fx = ::rtl::math::approxFloor( GetDouble());
                if (fx < 0.0)
                {
                    PushIllegalArgument();
                    return;
                }
                if (fx == 0.0 || fy == 0.0)
                    fy = 0.0;
                else
                    fy = fx * fy / ScGetGCD(fx, fy);
            }
            break;
            case svDoubleRef :
            case svRefList :
            {
                FormulaError nErr = FormulaError::NONE;
                PopDoubleRef( aRange, nParamCount, nRefInList);
                double nCellVal;
                ScValueIterator aValIter( mrContext, aRange, mnSubTotalFlags );
                if (aValIter.GetFirst(nCellVal, nErr))
                {
                    do
                    {
                        fx = ::rtl::math::approxFloor( nCellVal);
                        if (fx < 0.0)
                        {
                            PushIllegalArgument();
                            return;
                        }
                        if (fx == 0.0 || fy == 0.0)
                            fy = 0.0;
                        else
                            fy = fx * fy / ScGetGCD(fx, fy);
                    } while (nErr == FormulaError::NONE && aValIter.GetNext(nCellVal, nErr));
                }
                SetError(nErr);
            }
            break;
            case svMatrix :
            case svExternalSingleRef:
            case svExternalDoubleRef:
            {
                ScMatrixRef pMat = GetMatrix();
                if (pMat)
                {
                    SCSIZE nC, nR;
                    pMat->GetDimensions(nC, nR);
                    if (nC == 0 || nR == 0)
                        SetError(FormulaError::IllegalArgument);
                    else
                    {
                     double nVal = pMat->GetLcm();
                     fy = (nVal * fy ) / ScGetGCD(nVal, fy);
                    }
                }
            }
            break;
            default : SetError(FormulaError::IllegalParameter); break;
        }
    }
    PushDouble(fy);
}

void ScInterpreter::MakeMatNew(ScMatrixRef& rMat, SCSIZE nC, SCSIZE nR)
{
    rMat->SetErrorInterpreter( this);
    // A temporary matrix is mutable and ScMatrix::CloneIfConst() returns the
    // very matrix.
    rMat->SetMutable();
    SCSIZE nCols, nRows;
    rMat->GetDimensions( nCols, nRows);
    if ( nCols != nC || nRows != nR )
    {   // arbitrary limit of elements exceeded
        SetError( FormulaError::MatrixSize);
        rMat.reset();
    }
}

ScMatrixRef ScInterpreter::GetNewMat(SCSIZE nC, SCSIZE nR, bool bEmpty)
{
    ScMatrixRef pMat;
    if (bEmpty)
        pMat = new ScMatrix(nC, nR);
    else
        pMat = new ScMatrix(nC, nR, 0.0);
    MakeMatNew(pMat, nC, nR);
    return pMat;
}

ScMatrixRef ScInterpreter::GetNewMat(SCSIZE nC, SCSIZE nR, const std::vector<double>& ;rValues)
{
    ScMatrixRef pMat(new ScMatrix(nC, nR, rValues));
    MakeMatNew(pMat, nC, nR);
    return pMat;
}

ScMatrixRef ScInterpreter::CreateMatrixFromDoubleRef( const FormulaToken* pToken,
        SCCOL nCol1, SCROW nRow1, SCTAB nTab1,
        SCCOL nCol2, SCROW nRow2, SCTAB nTab2 )
{
    if (nTab1 != nTab2 || nGlobalError != FormulaError::NONE)
    {
        // Not a 2D matrix.
        SetError(FormulaError::IllegalParameter);
        return nullptr;
    }

    if (nTab1 == nTab2 && pToken)
    {
        const ScComplexRefData& rCRef = *pToken->GetDoubleRef();
        if (rCRef.IsTrimToData())
        {
            // Clamp the size of the matrix area to rows which actually contain data.
            // For e.g. SUM(IF over an entire column, this can make a big difference,
            // But let's not trim the empty space from the top or left as this matters
            // at least in matrix formulas involving IF().
            // Refer ScCompiler::AnnotateTrimOnDoubleRefs() where double-refs are
            // flagged for trimming.
            SCCOL nTempCol = nCol1;
            SCROW nTempRow = nRow1;
            mrDoc.ShrinkToDataArea(nTab1, nTempCol, nTempRow, nCol2, nRow2);
        }
    }

    SCSIZE nMatCols = static_cast<SCSIZE>(nCol2 - nCol1 + 1);
    SCSIZE nMatRows = static_cast<SCSIZE>(nRow2 - nRow1 + 1);

    if (!ScMatrix::IsSizeAllocatable( nMatCols, nMatRows))
    {
        SetError(FormulaError::MatrixSize);
        return nullptr;
    }

    ScTokenMatrixMap::const_iterator aIter;
    if (pToken && ((aIter = maTokenMatrixMap.find( pToken)) != maTokenMatrixMap.end()))
    {
        /* XXX casting const away here is ugly; ScMatrixToken (to which the
         * result of this function usually is assigned) should not be forced to
         * carry a ScConstMatrixRef though.
         * TODO: a matrix already stored in pTokenMatrixMap should be
         * read-only and have a copy-on-write mechanism. Previously all tokens
         * were modifiable so we're already better than before ... */
        return const_cast<FormulaToken*>((*aIter).second.get())->GetMatrix();
    }

    ScMatrixRef pMat = GetNewMat( nMatCols, nMatRows, true);
    if (!pMat || nGlobalError != FormulaError::NONE)
        return nullptr;

    if (!bCalcAsShown)
    {
        // Use fast array fill.
        mrDoc.FillMatrix(*pMat, nTab1, nCol1, nRow1, nCol2, nRow2);
    }
    else
    {
        // Use slower ScCellIterator to round values.

        // TODO: this probably could use CellBucket for faster storage, see
        // sc/source/core/data/column2.cxx and FillMatrixHandler, and then be
        // moved to a function on its own, and/or squeeze the rounding into a
        // similar FillMatrixHandler that would need to keep track of the cell
        // position then.

        // Set position where the next entry is expected.
        SCROW nNextRow = nRow1;
        SCCOL nNextCol = nCol1;
        // Set last position as if there was a previous entry.
        SCROW nThisRow = nRow2;
        SCCOL nThisCol = nCol1 - 1;

        ScCellIterator aCellIter( mrDoc, ScRange( nCol1, nRow1, nTab1, nCol2, nRow2, nTab2));
        for (bool bHas = aCellIter.first(); bHas; bHas = aCellIter.next())
        {
            nThisCol = aCellIter.GetPos().Col();
            nThisRow = aCellIter.GetPos().Row();
            if (nThisCol != nNextCol || nThisRow != nNextRow)
            {
                // Fill empty between iterator's positions.
                for ( ; nNextCol <= nThisCol; ++nNextCol)
                {
                    const SCSIZE nC = nNextCol - nCol1;
                    const SCSIZE nMatStopRow = ((nNextCol < nThisCol) ? nMatRows : nThisRow - nRow1);
                    for (SCSIZE nR = nNextRow - nRow1; nR < nMatStopRow; ++nR)
                    {
                        pMat->PutEmpty( nC, nR);
                    }
                    nNextRow = nRow1;
                }
            }
            if (nThisRow == nRow2)
            {
                nNextCol = nThisCol + 1;
                nNextRow = nRow1;
            }
            else
            {
                nNextCol = nThisCol;
                nNextRow = nThisRow + 1;
            }

            const SCSIZE nMatCol = static_cast<SCSIZE>(nThisCol - nCol1);
            const SCSIZE nMatRow = static_cast<SCSIZE>(nThisRow - nRow1);
            ScRefCellValue aCell( aCellIter.getRefCellValue());
            if (aCellIter.isEmpty() || aCell.hasEmptyValue())
            {
                pMat->PutEmpty( nMatCol, nMatRow);
            }
            else if (aCell.hasError())
            {
                pMat->PutError( aCell.getFormula()->GetErrCode(), nMatCol, nMatRow);
            }
            else if (aCell.hasNumeric())
            {
                double fVal = aCell.getValue();
                // CELLTYPE_FORMULA already stores the rounded value.
                if (aCell.getType() == CELLTYPE_VALUE)
                {
                    // TODO: this could be moved to ScCellIterator to take
                    // advantage of the faster ScAttrArray_IterGetNumberFormat.
                    const ScAddress aAdr( nThisCol, nThisRow, nTab1);
                    const sal_uInt32 nNumFormat = mrDoc.GetNumberFormat( mrContext, aAdr);
                    fVal = mrDoc.RoundValueAsShown( fVal, nNumFormat, &mrContext);
                }
                pMat->PutDouble( fVal, nMatCol, nMatRow);
            }
            else if (aCell.hasString())
            {
                pMat->PutString( mrStrPool.intern( aCell.getString(mrDoc)), nMatCol, nMatRow);
            }
            else
            {
                assert(!"aCell.what?");
                pMat->PutEmpty( nMatCol, nMatRow);
            }
        }

        // Fill empty if iterator's last position wasn't the end.
        if (nThisCol != nCol2 || nThisRow != nRow2)
        {
            for ( ; nNextCol <= nCol2; ++nNextCol)
            {
                SCSIZE nC = nNextCol - nCol1;
                for (SCSIZE nR = nNextRow - nRow1; nR < nMatRows; ++nR)
                {
                    pMat->PutEmpty( nC, nR);
                }
                nNextRow = nRow1;
            }
        }
    }

    if (pToken)
        maTokenMatrixMap.emplace(pToken, new ScMatrixToken( pMat));

    return pMat;
}

ScMatrixRef ScInterpreter::GetMatrix()
{
    ScMatrixRef pMat = nullptr;
    switch (GetRawStackType())
    {
        case svSingleRef :
        {
            ScAddress aAdr;
            PopSingleRef( aAdr );
            pMat = GetNewMat(1, 1);
            if (pMat)
            {
                ScRefCellValue aCell(mrDoc, aAdr);
                if (aCell.hasEmptyValue())
                    pMat->PutEmpty(0, 0);
                else if (aCell.hasNumeric())
                    pMat->PutDouble(GetCellValue(aAdr, aCell), 0);
                else
                {
                    svl::SharedString aStr;
                    GetCellString(aStr, aCell);
                    pMat->PutString(aStr, 0);
                }
            }
        }
        break;
        case svDoubleRef:
        {
            SCCOL nCol1, nCol2;
            SCROW nRow1, nRow2;
            SCTAB nTab1, nTab2;
            const formula::FormulaToken* p = sp ? pStack[sp-1] : nullptr;
            PopDoubleRef(nCol1, nRow1, nTab1, nCol2, nRow2, nTab2);
            pMat = CreateMatrixFromDoubleRef( p, nCol1, nRow1, nTab1,
                    nCol2, nRow2, nTab2);
        }
        break;
        case svMatrix:
            pMat = PopMatrix();
        break;
        case svError :
        case svMissing :
        case svDouble :
        {
            double fVal = GetDouble();
            pMat = GetNewMat( 1, 1);
            if ( pMat )
            {
                if ( nGlobalError != FormulaError::NONE )
                {
                    fVal = CreateDoubleError( nGlobalError);
                    nGlobalError = FormulaError::NONE;
                }
                pMat->PutDouble( fVal, 0);
            }
        }
        break;
        case svString :
        {
            svl::SharedString aStr = GetString();
            pMat = GetNewMat( 1, 1);
            if ( pMat )
            {
                if ( nGlobalError != FormulaError::NONE )
                {
                    double fVal = CreateDoubleError( nGlobalError);
                    pMat->PutDouble( fVal, 0);
                    nGlobalError = FormulaError::NONE;
                }
                else
                    pMat->PutString(aStr, 0);
            }
        }
        break;
        case svExternalSingleRef:
        {
            ScExternalRefCache::TokenRef pToken;
            PopExternalSingleRef(pToken);
            pMat = GetNewMat( 1, 1, true);
            if (!pMat)
                break;
            if (nGlobalError != FormulaError::NONE)
            {
                pMat->PutError( nGlobalError, 0, 0);
                nGlobalError = FormulaError::NONE;
                break;
            }
            switch (pToken->GetType())
            {
                case svError:
                    pMat->PutError( pToken->GetError(), 0, 0);
                break;
                case svDouble:
                    pMat->PutDouble( pToken->GetDouble(), 0, 0);
                break;
                case svString:
                    pMat->PutString( pToken->GetString(), 0, 0);
                break;
                default:
                    ;   // nothing, empty element matrix
            }
        }
        break;
        case svExternalDoubleRef:
            PopExternalDoubleRef(pMat);
        break;
        default:
            PopError();
            SetError( FormulaError::IllegalArgument);
        break;
    }
    return pMat;
}

ScMatrixRef ScInterpreter::GetMatrix( short & rParam, size_t & rRefInList )
{
    switch (GetRawStackType())
    {
        case svRefList:
            {
                ScRange aRange( ScAddress::INITIALIZE_INVALID );
                PopDoubleRef( aRange, rParam, rRefInList);
                if (nGlobalError != FormulaError::NONE)
                    return nullptr;

                SCCOL nCol1, nCol2;
                SCROW nRow1, nRow2;
                SCTAB nTab1, nTab2;
                aRange.GetVars( nCol1, nRow1, nTab1, nCol2, nRow2, nTab2);
                return CreateMatrixFromDoubleRef( nullptr, nCol1, nRow1, nTab1, nCol2, nRow2, nTab2);
            }
        break;
        default:
            return GetMatrix();
    }
}

sc::RangeMatrix ScInterpreter::GetRangeMatrix()
{
    sc::RangeMatrix aRet;
    switch (GetRawStackType())
    {
        case svMatrix:
            aRet = PopRangeMatrix();
        break;
        default:
            aRet.mpMat = GetMatrix();
    }
    return aRet;
}

void ScInterpreter::ScMatValue()
{
    if ( !MustHaveParamCount( GetByte(), 3 ) )
        return;

    // 0 to count-1
    // Theoretically we could have GetSize() instead of GetUInt32(), but
    // really, practically ...
    SCSIZE nR = static_cast<SCSIZE>(GetUInt32());
    SCSIZE nC = static_cast<SCSIZE>(GetUInt32());
    if (nGlobalError != FormulaError::NONE)
    {
        PushError( nGlobalError);
        return;
    }
    switch (GetStackType())
    {
        case svSingleRef :
        {
            ScAddress aAdr;
            PopSingleRef( aAdr );
            ScRefCellValue aCell(mrDoc, aAdr);
            if (aCell.getType() == CELLTYPE_FORMULA)
            {
                FormulaError nErrCode = aCell.getFormula()->GetErrCode();
                if (nErrCode != FormulaError::NONE)
                    PushError( nErrCode);
                else
                {
                    const ScMatrix* pMat = aCell.getFormula()->GetMatrix();
                    CalculateMatrixValue(pMat,nC,nR);
                }
            }
            else
                PushIllegalParameter();
        }
        break;
        case svDoubleRef :
        {
            SCCOL nCol1;
            SCROW nRow1;
            SCTAB nTab1;
            SCCOL nCol2;
            SCROW nRow2;
            SCTAB nTab2;
            PopDoubleRef(nCol1, nRow1, nTab1, nCol2, nRow2, nTab2);
            if (nCol2 - nCol1 >= static_cast<SCCOL>(nR) &&
                    nRow2 - nRow1 >= static_cast<SCROW>(nC) &&
                    nTab1 == nTab2)
            {
                ScAddress aAdr( sal::static_int_cast<SCCOL>( nCol1 + nR ),
                                sal::static_int_cast<SCROW>( nRow1 + nC ), nTab1 );
                ScRefCellValue aCell(mrDoc, aAdr);
                if (aCell.hasNumeric())
                    PushDouble(GetCellValue(aAdr, aCell));
                else
                {
                    svl::SharedString aStr;
                    GetCellString(aStr, aCell);
                    PushString(aStr);
                }
            }
            else
                PushNoValue();
        }
        break;
        case svMatrix:
        {
            ScMatrixRef pMat = PopMatrix();
            CalculateMatrixValue(pMat.get(),nC,nR);
        }
        break;
        default:
            PopError();
            PushIllegalParameter();
        break;
    }
}
void ScInterpreter::CalculateMatrixValue(const ScMatrix* pMat,SCSIZE nC,SCSIZE nR)
{
    if (pMat)
    {
        SCSIZE nCl, nRw;
        pMat->GetDimensions(nCl, nRw);
        if (nC < nCl && nR < nRw)
        {
            const ScMatrixValue nMatVal = pMat->Get( nC, nR);
            ScMatValType nMatValType = nMatVal.nType;
            if (ScMatrix::IsNonValueType( nMatValType))
                PushString( nMatVal.GetString() );
            else
                PushDouble(nMatVal.fVal);
                // also handles DoubleError
        }
        else
            PushNoValue();
    }
    else
        PushNoValue();
}

void ScInterpreter::ScEMat()
{
    if ( !MustHaveParamCount( GetByte(), 1 ) )
        return;

    SCSIZE nDim = static_cast<SCSIZE>(GetUInt32());
    if (nGlobalError != FormulaError::NONE || nDim == 0)
        PushIllegalArgument();
    else if (!ScMatrix::IsSizeAllocatable( nDim, nDim))
        PushError( FormulaError::MatrixSize);
    else
    {
        ScMatrixRef pRMat = GetNewMat(nDim, nDim, /*bEmpty*/true);
        if (pRMat)
        {
            MEMat(pRMat, nDim);
            PushMatrix(pRMat);
        }
        else
            PushIllegalArgument();
    }
}

void ScInterpreter::MEMat(const ScMatrixRef& mM, SCSIZE n)
{
    mM->FillDouble(0.0, 0, 0, n-1, n-1);
    for (SCSIZE i = 0; i < n; i++)
        mM->PutDouble(1.0, i, i);
}

/* Matrix LUP decomposition according to the pseudocode of "Introduction to
 * Algorithms" by Cormen, Leiserson, Rivest, Stein.
 *
 * Added scaling for numeric stability.
 *
 * Given an n x n nonsingular matrix A, find a permutation matrix P, a unit
 * lower-triangular matrix L, and an upper-triangular matrix U such that PA=LU.
 * Compute L and U "in place" in the matrix A, the original content is
 * destroyed. Note that the diagonal elements of the U triangular matrix
 * replace the diagonal elements of the L-unit matrix (that are each ==1). The
 * permutation matrix P is an array, where P[i]=j means that the i-th row of P
 * contains a 1 in column j. Additionally keep track of the number of
 * permutations (row exchanges).
 *
 * Returns 0 if a singular matrix is encountered, else +1 if an even number of
 * permutations occurred, or -1 if odd, which is the sign of the determinant.
 * This may be used to calculate the determinant by multiplying the sign with
 * the product of the diagonal elements of the LU matrix.
 */
static int lcl_LUP_decompose( ScMatrix* mA, const SCSIZE n,
        ::std::vector< SCSIZE> & P )
{
    int nSign = 1;
    // Find scale of each row.
    ::std::vector< double> aScale(n);
    for (SCSIZE i=0; i < n; ++i)
    {
        double fMax = 0.0;
        for (SCSIZE j=0; j < n; ++j)
        {
            double fTmp = fabs( mA->GetDouble( j, i));
            if (fMax < fTmp)
                fMax = fTmp;
        }
        if (fMax == 0.0)
            return 0;       // singular matrix
        aScale[i] = 1.0 / fMax;
    }
    // Represent identity permutation, P[i]=i
    for (SCSIZE i=0; i < n; ++i)
        P[i] = i;
    // "Recursion" on the diagonal.
    SCSIZE l = n - 1;
    for (SCSIZE k=0; k < l; ++k)
    {
        // Implicit pivoting. With the scale found for a row, compare values of
        // a column and pick largest.
        double fMax = 0.0;
        double fScale = aScale[k];
        SCSIZE kp = k;
        for (SCSIZE i = k; i < n; ++i)
        {
            double fTmp = fScale * fabs( mA->GetDouble( k, i));
            if (fMax < fTmp)
            {
                fMax = fTmp;
                kp = i;
            }
        }
        if (fMax == 0.0)
            return 0;       // singular matrix
        // Swap rows. The pivot element will be at mA[k,kp] (row,col notation)
        if (k != kp)
        {
            // permutations
            SCSIZE nTmp = P[k];
            P[k]        = P[kp];
            P[kp]       = nTmp;
            nSign       = -nSign;
            // scales
            double fTmp = aScale[k];
            aScale[k]   = aScale[kp];
            aScale[kp]  = fTmp;
            // elements
            for (SCSIZE i=0; i < n; ++i)
            {
                double fMatTmp = mA->GetDouble( i, k);
                mA->PutDouble( mA->GetDouble( i, kp), i, k);
                mA->PutDouble( fMatTmp, i, kp);
            }
        }
        // Compute Schur complement.
        for (SCSIZE i = k+1; i < n; ++i)
        {
            double fNum = mA->GetDouble( k, i);
            double fDen = mA->GetDouble( k, k);
            mA->PutDouble( fNum/fDen, k, i);
            for (SCSIZE j = k+1; j < n; ++j)
                mA->PutDouble( ( mA->GetDouble( j, i) * fDen  -
                            fNum * mA->GetDouble( j, k) ) / fDen, j, i);
        }
    }
#ifdef DEBUG_SC_LUP_DECOMPOSITION
    fprintf( stderr, "\n%s\n", "lcl_LUP_decompose(): LU");
    for (SCSIZE i=0; i < n; ++i)
    {
        for (SCSIZE j=0; j < n; ++j)
            fprintf( stderr, "%8.2g ", mA->GetDouble( j, i));
        fprintf( stderr, "\n%s\n", "");
    }
    fprintf( stderr, "\n%s\n", "lcl_LUP_decompose(): P");
    for (SCSIZE j=0; j < n; ++j)
        fprintf( stderr, "%5u ", (unsigned)P[j]);
    fprintf( stderr, "\n%s\n", "");
#endif

    bool bSingular=false;
    for (SCSIZE i=0; i<n && !bSingular; i++)
        bSingular = (mA->GetDouble(i,i)) == 0.0;
    if (bSingular)
        nSign = 0;

    return nSign;
}

/* Solve a LUP decomposed equation Ax=b. LU is a combined matrix of L and U
 * triangulars and P the permutation vector as obtained from
 * lcl_LUP_decompose(). B is the right-hand side input vector, X is used to
 * return the solution vector.
 */
static void lcl_LUP_solve( const ScMatrix* mLU, const SCSIZE n,
        const ::std::vector< SCSIZE> & P, const ::std::vector< double> & B,
        ::std::vector< double> & X )
{
    SCSIZE nFirst = SCSIZE_MAX;
    // Ax=b => PAx=Pb, with decomposition LUx=Pb.
    // Define y=Ux and solve for y in Ly=Pb using forward substitution.
    for (SCSIZE i=0; i < n; ++i)
    {
        KahanSum fSum = B[P[i]];
        // Matrix inversion comes with a lot of zeros in the B vectors, we
        // don't have to do all the computing with results multiplied by zero.
        // Until then, simply lookout for the position of the first nonzero
        // value.
        if (nFirst != SCSIZE_MAX)
        {
            for (SCSIZE j = nFirst; j < i; ++j)
                fSum -= mLU->GetDouble( j, i) * X[j];         // X[j] === y[j]
        }
        else if (fSum != 0)
            nFirst = i;
        X[i] = fSum.get();                                    // X[i] === y[i]
    }
    // Solve for x in Ux=y using back substitution.
    for (SCSIZE i = n; i--; )
    {
        KahanSum fSum = X[i];                                 // X[i] === y[i]
        for (SCSIZE j = i+1; j < n; ++j)
            fSum -= mLU->GetDouble( j, i) * X[j];             // X[j] === x[j]
        X[i] = fSum.get() / mLU->GetDouble( i, i);            // X[i] === x[i]
    }
#ifdef DEBUG_SC_LUP_DECOMPOSITION
    fprintf( stderr, "\n%s\n", "lcl_LUP_solve():");
    for (SCSIZE i=0; i < n; ++i)
        fprintf( stderr, "%8.2g ", X[i]);
    fprintf( stderr, "%s\n", "");
#endif
}

void ScInterpreter::ScMatDet()
{
    if ( !MustHaveParamCount( GetByte(), 1 ) )
        return;

    ScMatrixRef pMat = GetMatrix();
    if (!pMat)
    {
        PushIllegalParameter();
        return;
    }
    if ( !pMat->IsNumeric() )
    {
        PushNoValue();
        return;
    }
    SCSIZE nC, nR;
    pMat->GetDimensions(nC, nR);
    if ( nC != nR || nC == 0 )
        PushIllegalArgument();
    else if (!ScMatrix::IsSizeAllocatable( nC, nR))
        PushError( FormulaError::MatrixSize);
    else
    {
        // LUP decomposition is done inplace, use copy.
        ScMatrixRef xLU = pMat->Clone();
        if (!xLU)
            PushError( FormulaError::CodeOverflow);
        else
        {
            ::std::vector< SCSIZE> P(nR);
            int nDetSign = lcl_LUP_decompose( xLU.get(), nR, P);
            if (!nDetSign)
                PushInt(0);     // singular matrix
            else
            {
                // In an LU matrix the determinant is simply the product of
                // all diagonal elements.
                double fDet = nDetSign;
                for (SCSIZE i=0; i < nR; ++i)
                    fDet *= xLU->GetDouble( i, i);
                PushDouble( fDet);
            }
        }
    }
}

void ScInterpreter::ScMatInv()
{
    if ( !MustHaveParamCount( GetByte(), 1 ) )
        return;

    ScMatrixRef pMat = GetMatrix();
    if (!pMat)
    {
        PushIllegalParameter();
        return;
    }
    if ( !pMat->IsNumeric() )
    {
        PushNoValue();
        return;
    }
    SCSIZE nC, nR;
    pMat->GetDimensions(nC, nR);

    if (ScCalcConfig::isOpenCLEnabled())
    {
        sc::FormulaGroupInterpreter *pInterpreter = sc::FormulaGroupInterpreter::getStatic();
        if (pInterpreter != nullptr)
        {
            ScMatrixRef xResMat = pInterpreter->inverseMatrix(*pMat);
            if (xResMat)
            {
                PushMatrix(xResMat);
                return;
            }
        }
    }

    if ( nC != nR || nC == 0 )
        PushIllegalArgument();
    else if (!ScMatrix::IsSizeAllocatable( nC, nR))
        PushError( FormulaError::MatrixSize);
    else
    {
        // LUP decomposition is done inplace, use copy.
        ScMatrixRef xLU = pMat->Clone();
        // The result matrix.
        ScMatrixRef xY = GetNewMat( nR, nR, /*bEmpty*/true );
        if (!xLU || !xY)
            PushError( FormulaError::CodeOverflow);
        else
        {
            ::std::vector< SCSIZE> P(nR);
            int nDetSign = lcl_LUP_decompose( xLU.get(), nR, P);
            if (!nDetSign)
                PushIllegalArgument();
            else
            {
                // Solve equation for each column.
                ::std::vector< double> B(nR);
                ::std::vector< double> X(nR);
                for (SCSIZE j=0; j < nR; ++j)
                {
                    for (SCSIZE i=0; i < nR; ++i)
                        B[i] = 0.0;
                    B[j] = 1.0;
                    lcl_LUP_solve( xLU.get(), nR, P, B, X);
                    for (SCSIZE i=0; i < nR; ++i)
                        xY->PutDouble( X[i], j, i);
                }
#ifdef DEBUG_SC_LUP_DECOMPOSITION
                /* Possible checks for ill-condition:
                 * 1. Scale matrix, invert scaled matrix. If there are
                 *    elements of the inverted matrix that are several
                 *    orders of magnitude greater than 1 =>
                 *    ill-conditioned.
                 *    Just how much is "several orders"?
                 * 2. Invert the inverted matrix and assess whether the
                 *    result is sufficiently close to the original matrix.
                 *    If not => ill-conditioned.
                 *    Just what is sufficient?
                 * 3. Multiplying the inverse by the original matrix should
                 *    produce a result sufficiently close to the identity
                 *    matrix.
                 *    Just what is sufficient?
                 *
                 * The following is #3.
                 */
                const double fInvEpsilon = 1.0E-7;
                ScMatrixRef xR = GetNewMat( nR, nR);
                if (xR)
                {
                    ScMatrix* pR = xR.get();
                    lcl_MFastMult( pMat, xY.get(), pR, nR, nR, nR);
                    fprintf( stderr, "\n%s\n", "ScMatInv(): mult-identity");
                    for (SCSIZE i=0; i < nR; ++i)
                    {
                        for (SCSIZE j=0; j < nR; ++j)
                        {
                            double fTmp = pR->GetDouble( j, i);
                            fprintf( stderr, "%8.2g ", fTmp);
                            if (fabs( fTmp - (i == j)) > fInvEpsilon)
                                SetError( FormulaError::IllegalArgument);
                        }
                    fprintf( stderr, "\n%s\n", "");
                    }
                }
#endif
                if (nGlobalError != FormulaError::NONE)
                    PushError( nGlobalError);
                else
                    PushMatrix( xY);
            }
        }
    }
}

void ScInterpreter::ScMatMult()
{
    if ( !MustHaveParamCount( GetByte(), 2 ) )
        return;

    ScMatrixRef pMat2 = GetMatrix();
    ScMatrixRef pMat1 = GetMatrix();
    ScMatrixRef pRMat;
    if (pMat1 && pMat2)
    {
        if ( pMat1->IsNumeric() && pMat2->IsNumeric() )
        {
            SCSIZE nC1, nC2;
            SCSIZE nR1, nR2;
            pMat1->GetDimensions(nC1, nR1);
            pMat2->GetDimensions(nC2, nR2);
            if (nC1 != nR2)
                PushIllegalArgument();
            else
            {
                pRMat = GetNewMat(nC2, nR1, /*bEmpty*/true);
                if (pRMat)
                {
                    KahanSum fSum;
                    for (SCSIZE i = 0; i < nR1; i++)
                    {
                        for (SCSIZE j = 0; j < nC2; j++)
                        {
                            fSum = 0.0;
                            for (SCSIZE k = 0; k < nC1; k++)
                            {
                                fSum += pMat1->GetDouble(k,i)*pMat2->GetDouble(j,k);
                            }
                            pRMat->PutDouble(fSum.get(), j, i);
                        }
                    }
                    PushMatrix(pRMat);
                }
                else
                    PushIllegalArgument();
            }
        }
        else
            PushNoValue();
    }
    else
        PushIllegalParameter();
}

void ScInterpreter::ScMatSequence()
{
    sal_uInt8 nParamCount = GetByte();
    if (!MustHaveParamCount(nParamCount, 1, 4))
        return;

    // 4th argument is the step number (optional)
    double nSteps = 1.0;
    if (nParamCount == 4)
        nSteps = GetDoubleWithDefault(nSteps);

    // 3d argument is the start number (optional)
    double nStart = 1.0;
    if (nParamCount >= 3)
        nStart = GetDoubleWithDefault(nStart);

    // 2nd argument is the col number (optional)
    sal_Int32 nColumns = 1;
    if (nParamCount >= 2)
    {
        nColumns = GetInt32WithDefault(nColumns);
        if (nColumns < 1)
        {
            PushIllegalArgument();
            return;
        }
    }

    // 1st argument is the row number (required)
    sal_Int32 nRows = GetInt32WithDefault(1);
    if (nRows < 1)
    {
        PushIllegalArgument();
        return;
    }

    if (nGlobalError != FormulaError::NONE)
    {
        PushError(nGlobalError);
        return;
    }

    size_t nMatrixSize = static_cast<size_t>(nColumns) * nRows;
    ScMatrixRef pResMat = GetNewMat(nColumns, nRows, /*bEmpty*/true);
    for (size_t iPos = 0; iPos < nMatrixSize; iPos++)
    {
        pResMat->PutDoubleTrans(nStart, iPos);
        nStart = nStart + nSteps;
    }

    if (pResMat)
    {
        PushMatrix(pResMat);
    }
    else
    {
        PushIllegalParameter();
        return;
    }
}

void ScInterpreter::ScMatTrans()
{
    if ( !MustHaveParamCount( GetByte(), 1 ) )
        return;

    ScMatrixRef pMat = GetMatrix();
    ScMatrixRef pRMat;
    if (pMat)
    {
        SCSIZE nC, nR;
        pMat->GetDimensions(nC, nR);
        pRMat = GetNewMat(nR, nC, /*bEmpty*/true);
        if ( pRMat )
        {
            pMat->MatTrans(*pRMat);
            PushMatrix(pRMat);
        }
        else
            PushIllegalArgument();
    }
    else
        PushIllegalParameter();
}

/** Minimum extent of one result matrix dimension.
    For a row or column vector to be replicated the larger matrix dimension is
    returned, else the smaller dimension.
 */
static SCSIZE lcl_GetMinExtent( SCSIZE n1, SCSIZE n2 )
{
    if (n1 == 1)
        return n2;
    else if (n2 == 1)
        return n1;
    else if (n1 < n2)
        return n1;
    else
        return n2;
}

static ScMatrixRef lcl_MatrixCalculation(
    const ScMatrix& rMat1, const ScMatrix& rMat2, ScInterpreter* pInterpreter, const ScMatrix::CalculateOpFunction& Op)
{
    SCSIZE nC1, nC2, nMinC;
    SCSIZE nR1, nR2, nMinR;
    rMat1.GetDimensions(nC1, nR1);
    rMat2.GetDimensions(nC2, nR2);
    nMinC = lcl_GetMinExtent( nC1, nC2);
    nMinR = lcl_GetMinExtent( nR1, nR2);
    ScMatrixRef xResMat = pInterpreter->GetNewMat(nMinC, nMinR, /*bEmpty*/true);
    if (xResMat)
        xResMat->ExecuteBinaryOp(nMinC, nMinR, rMat1, rMat2, pInterpreter, Op);
    return xResMat;
}

ScMatrixRef ScInterpreter::MatConcat(const ScMatrixRef& pMat1, const ScMatrixRef&&nbsp;pMat2)
{
    SCSIZE nC1, nC2, nMinC;
    SCSIZE nR1, nR2, nMinR;
    pMat1->GetDimensions(nC1, nR1);
    pMat2->GetDimensions(nC2, nR2);
    nMinC = lcl_GetMinExtent( nC1, nC2);
    nMinR = lcl_GetMinExtent( nR1, nR2);
    ScMatrixRef xResMat = GetNewMat(nMinC, nMinR, /*bEmpty*/true);
    if (xResMat)
    {
        xResMat->MatConcat(nMinC, nMinR, pMat1, pMat2, mrContext, mrDoc.GetSharedStringPool());
    }
    return xResMat;
}

// for DATE, TIME, DATETIME, DURATION
static void lcl_GetDiffDateTimeFmtType( SvNumFormatType& nFuncFmt, SvNumFormatType nFmt1, SvNumFormatType nFmt2 )
{
    if ( nFmt1 == SvNumFormatType::UNDEFINED && nFmt2 == SvNumFormatType::UNDEFINED )
        return;

    if ( nFmt1 == nFmt2 )
    {
        if ( nFmt1 == SvNumFormatType::TIME || nFmt1 == SvNumFormatType::DATETIME
                || nFmt1 == SvNumFormatType::DURATION )
            nFuncFmt = SvNumFormatType::DURATION;   // times result in time duration
        // else: nothing special, number (date - date := days)
    }
    else if ( nFmt1 == SvNumFormatType::UNDEFINED )
        nFuncFmt = nFmt2;   // e.g. date + days := date
    else if ( nFmt2 == SvNumFormatType::UNDEFINED )
        nFuncFmt = nFmt1;
    else
    {
        if ( nFmt1 == SvNumFormatType::DATE || nFmt2 == SvNumFormatType::DATE ||
            nFmt1 == SvNumFormatType::DATETIME || nFmt2 == SvNumFormatType::DATETIME )
        {
            if ( nFmt1 == SvNumFormatType::TIME || nFmt2 == SvNumFormatType::TIME )
                nFuncFmt = SvNumFormatType::DATETIME;   // date + time
        }
    }
}

void ScInterpreter::ScAdd()
{
    CalculateAddSub(false);
}

void ScInterpreter::CalculateAddSub(bool _bSub)
{
    ScMatrixRef pMat1 = nullptr;
    ScMatrixRef pMat2 = nullptr;
    double fVal1 = 0.0, fVal2 = 0.0;
    SvNumFormatType nFmt1, nFmt2;
    nFmt1 = nFmt2 = SvNumFormatType::UNDEFINED;
    bool bDuration = false;
    SvNumFormatType nFmtCurrencyType = nCurFmtType;
    sal_uLong nFmtCurrencyIndex = nCurFmtIndex;
    SvNumFormatType nFmtPercentType = nCurFmtType;
    if ( GetStackType() == svMatrix )
        pMat2 = GetMatrix();
    else
    {
        fVal2 = GetDouble();
        switch ( nCurFmtType )
        {
            case SvNumFormatType::DATE :
            case SvNumFormatType::TIME :
            case SvNumFormatType::DATETIME :
            case SvNumFormatType::DURATION :
                nFmt2 = nCurFmtType;
                bDuration = true;
            break;
            case SvNumFormatType::CURRENCY :
                nFmtCurrencyType = nCurFmtType;
                nFmtCurrencyIndex = nCurFmtIndex;
            break;
            case SvNumFormatType::PERCENT :
                nFmtPercentType = SvNumFormatType::PERCENT;
            break;
            default: break;
        }
    }
    if ( GetStackType() == svMatrix )
        pMat1 = GetMatrix();
    else
    {
        fVal1 = GetDouble();
        switch ( nCurFmtType )
        {
            case SvNumFormatType::DATE :
            case SvNumFormatType::TIME :
            case SvNumFormatType::DATETIME :
            case SvNumFormatType::DURATION :
                nFmt1 = nCurFmtType;
                bDuration = true;
            break;
            case SvNumFormatType::CURRENCY :
                nFmtCurrencyType = nCurFmtType;
                nFmtCurrencyIndex = nCurFmtIndex;
            break;
            case SvNumFormatType::PERCENT :
                nFmtPercentType = SvNumFormatType::PERCENT;
            break;
            default: break;
        }
    }
    if (pMat1 && pMat2)
    {
        ScMatrixRef pResMat;
        if ( _bSub )
        {
            pResMat = lcl_MatrixCalculation( *pMat1, *pMat2, this, MatrixSub);
        }
        else
        {
            pResMat = lcl_MatrixCalculation( *pMat1, *pMat2, this, MatrixAdd);
        }

        if (!pResMat)
            PushNoValue();
        else
            PushMatrix(pResMat);
    }
    else if (pMat1 || pMat2)
    {
        double fVal;
        bool bFlag;
        ScMatrixRef pMat = std::move(pMat1);
        if (!pMat)
        {
            fVal = fVal1;
            pMat = std::move(pMat2);
            bFlag = true;           // double - Matrix
        }
        else
        {
            fVal = fVal2;
            bFlag = false;          // Matrix - double
        }
        SCSIZE nC, nR;
        pMat->GetDimensions(nC, nR);
        ScMatrixRef pResMat = GetNewMat(nC, nR, true);
        if (pResMat)
        {
            if (_bSub)
            {
                pMat->SubOp( bFlag, fVal, *pResMat);
            }
            else
            {
                pMat->AddOp( fVal, *pResMat);
            }
            PushMatrix(pResMat);
        }
        else
            PushIllegalArgument();
    }
    else
    {
        // Determine nFuncFmtType type before PushDouble().
        if ( nFmtCurrencyType == SvNumFormatType::CURRENCY )
        {
            nFuncFmtType = nFmtCurrencyType;
            nFuncFmtIndex = nFmtCurrencyIndex;
        }
        else
        {
            lcl_GetDiffDateTimeFmtType( nFuncFmtType, nFmt1, nFmt2 );
            if (nFmtPercentType == SvNumFormatType::PERCENT && nFuncFmtType == SvNumFormatType::NUMBER)
                nFuncFmtType = SvNumFormatType::PERCENT;
        }
        if ((nFuncFmtType == SvNumFormatType::DURATION || bDuration)
                && ((_bSub && std::fabs(fVal1 - fVal2) <= SAL_MAX_INT32)
                    || (!_bSub && std::fabs(fVal1 + fVal2) <= SAL_MAX_INT32)))
        {
            // Limit to microseconds resolution on date inflicted or duration
            // values of 24 hours or more.
            const sal_uInt64 nEpsilon = ((std::fabs(fVal1) >= 1.0 || std::fabs(fVal2) >= 1.0) ?
                    ::tools::Duration::kAccuracyEpsilonNanosecondsMicroseconds :
                    ::tools::Duration::kAccuracyEpsilonNanoseconds);
            if (_bSub)
                PushDouble( ::tools::Duration( fVal1 - fVal2, nEpsilon).GetInDays());
            else
                PushDouble( ::tools::Duration( fVal1 + fVal2, nEpsilon).GetInDays());
        }
        else
        {
            if (_bSub)
                PushDouble( ::rtl::math::approxSub( fVal1, fVal2 ) );
            else
                PushDouble( ::rtl::math::approxAdd( fVal1, fVal2 ) );
        }
    }
}

void ScInterpreter::ScAmpersand()
{
    ScMatrixRef pMat1 = nullptr;
    ScMatrixRef pMat2 = nullptr;
    OUString sStr1, sStr2;
    if ( GetStackType() == svMatrix )
        pMat2 = GetMatrix();
    else
        sStr2 = GetString().getString();
    if ( GetStackType() == svMatrix )
        pMat1 = GetMatrix();
    else
        sStr1 = GetString().getString();
    if (pMat1 && pMat2)
    {
        ScMatrixRef pResMat = MatConcat(pMat1, pMat2);
        if (!pResMat)
            PushNoValue();
        else
            PushMatrix(pResMat);
    }
    else if (pMat1 || pMat2)
    {
        OUString sStr;
        bool bFlag;
        ScMatrixRef pMat = std::move(pMat1);
        if (!pMat)
        {
            sStr = sStr1;
            pMat = std::move(pMat2);
            bFlag = true;           // double - Matrix
        }
        else
        {
            sStr = sStr2;
            bFlag = false;          // Matrix - double
        }
        SCSIZE nC, nR;
        pMat->GetDimensions(nC, nR);
        ScMatrixRef pResMat = GetNewMat(nC, nR, /*bEmpty*/true);
        if (pResMat)
        {
            if (nGlobalError != FormulaError::NONE)
            {
                for (SCSIZE i = 0; i < nC; ++i)
                    for (SCSIZE j = 0; j < nR; ++j)
                        pResMat->PutError( nGlobalError, i, j);
            }
            else if (bFlag)
            {
                for (SCSIZE i = 0; i < nC; ++i)
                    for (SCSIZE j = 0; j < nR; ++j)
                    {
                        FormulaError nErr = pMat->GetErrorIfNotString( i, j);
                        if (nErr != FormulaError::NONE)
                            pResMat->PutError( nErr, i, j);
                        else
                        {
                            OUString aTmp = sStr + pMat->GetString(mrContext, i, j).getString();
                            pResMat->PutString(mrStrPool.intern(aTmp), i, j);
                        }
                    }
            }
            else
            {
                for (SCSIZE i = 0; i < nC; ++i)
                    for (SCSIZE j = 0; j < nR; ++j)
                    {
                        FormulaError nErr = pMat->GetErrorIfNotString( i, j);
                        if (nErr != FormulaError::NONE)
                            pResMat->PutError( nErr, i, j);
                        else
                        {
                            OUString aTmp = pMat->GetString(mrContext, i, j).getString() + sStr;
                            pResMat->PutString(mrStrPool.intern(aTmp), i, j);
                        }
                    }
            }
            PushMatrix(pResMat);
        }
        else
            PushIllegalArgument();
    }
    else
    {
        if ( CheckStringResultLen( sStr1, sStr2.getLength() ) )
            sStr1 += sStr2;
        PushString(sStr1);
    }
}

void ScInterpreter::ScSub()
{
    CalculateAddSub(true);
}

void ScInterpreter::ScMul()
{
    ScMatrixRef pMat1 = nullptr;
    ScMatrixRef pMat2 = nullptr;
    double fVal1 = 0.0, fVal2 = 0.0;
    SvNumFormatType nFmtCurrencyType = nCurFmtType;
    sal_uLong nFmtCurrencyIndex = nCurFmtIndex;
    if ( GetStackType() == svMatrix )
        pMat2 = GetMatrix();
    else
    {
        fVal2 = GetDouble();
        switch ( nCurFmtType )
        {
            case SvNumFormatType::CURRENCY :
                nFmtCurrencyType = nCurFmtType;
                nFmtCurrencyIndex = nCurFmtIndex;
            break;
            default: break;
        }
    }
    if ( GetStackType() == svMatrix )
        pMat1 = GetMatrix();
    else
    {
        fVal1 = GetDouble();
        switch ( nCurFmtType )
        {
            case SvNumFormatType::CURRENCY :
                nFmtCurrencyType = nCurFmtType;
                nFmtCurrencyIndex = nCurFmtIndex;
            break;
            default: break;
        }
    }
    if (pMat1 && pMat2)
    {
        ScMatrixRef pResMat = lcl_MatrixCalculation( *pMat1, *pMat2, this, MatrixMul);
        if (!pResMat)
            PushNoValue();
        else
            PushMatrix(pResMat);
    }
    else if (pMat1 || pMat2)
    {
        double fVal;
        ScMatrixRef pMat = std::move(pMat1);
        if (!pMat)
        {
            fVal = fVal1;
            pMat = std::move(pMat2);
        }
        else
            fVal = fVal2;
        SCSIZE nC, nR;
        pMat->GetDimensions(nC, nR);
        ScMatrixRef pResMat = GetNewMat(nC, nR, /*bEmpty*/true);
        if (pResMat)
        {
            pMat->MulOp( fVal, *pResMat);
            PushMatrix(pResMat);
        }
        else
            PushIllegalArgument();
    }
    else
    {
        // Determine nFuncFmtType type before PushDouble().
        if ( nFmtCurrencyType == SvNumFormatType::CURRENCY )
        {
            nFuncFmtType = nFmtCurrencyType;
            nFuncFmtIndex = nFmtCurrencyIndex;
        }
        PushDouble(fVal1 * fVal2);
    }
}

void ScInterpreter::ScDiv()
{
    ScMatrixRef pMat1 = nullptr;
    ScMatrixRef pMat2 = nullptr;
    double fVal1 = 0.0, fVal2 = 0.0;
    SvNumFormatType nFmtCurrencyType = nCurFmtType;
    sal_uLong nFmtCurrencyIndex = nCurFmtIndex;
    SvNumFormatType nFmtCurrencyType2 = SvNumFormatType::UNDEFINED;
    if ( GetStackType() == svMatrix )
        pMat2 = GetMatrix();
    else
    {
        fVal2 = GetDouble();
        // do not take over currency, 123kg/456USD is not USD
        nFmtCurrencyType2 = nCurFmtType;
    }
    if ( GetStackType() == svMatrix )
        pMat1 = GetMatrix();
    else
    {
        fVal1 = GetDouble();
        switch ( nCurFmtType )
        {
            case SvNumFormatType::CURRENCY :
                nFmtCurrencyType = nCurFmtType;
                nFmtCurrencyIndex = nCurFmtIndex;
            break;
            default: break;
        }
    }
    if (pMat1 && pMat2)
    {
        ScMatrixRef pResMat = lcl_MatrixCalculation( *pMat1, *pMat2, this, MatrixDiv);
        if (!pResMat)
            PushNoValue();
        else
            PushMatrix(pResMat);
    }
    else if (pMat1 || pMat2)
    {
        double fVal;
        bool bFlag;
        ScMatrixRef pMat = std::move(pMat1);
        if (!pMat)
        {
            fVal = fVal1;
            pMat = std::move(pMat2);
            bFlag = true;           // double - Matrix
        }
        else
        {
            fVal = fVal2;
            bFlag = false;          // Matrix - double
        }
        SCSIZE nC, nR;
        pMat->GetDimensions(nC, nR);
        ScMatrixRef pResMat = GetNewMat(nC, nR, /*bEmpty*/true);
        if (pResMat)
        {
            pMat->DivOp( bFlag, fVal, *pResMat);
            PushMatrix(pResMat);
        }
        else
            PushIllegalArgument();
    }
    else
    {
        // Determine nFuncFmtType type before PushDouble().
        if (    nFmtCurrencyType  == SvNumFormatType::CURRENCY &&
                nFmtCurrencyType2 != SvNumFormatType::CURRENCY)
        {   // even USD/USD is not USD
            nFuncFmtType = nFmtCurrencyType;
            nFuncFmtIndex = nFmtCurrencyIndex;
        }
        PushDouble( div( fVal1, fVal2) );
    }
}

void ScInterpreter::ScPower()
{
    if ( MustHaveParamCount( GetByte(), 2 ) )
        ScPow();
}

void ScInterpreter::ScPow()
{
    ScMatrixRef pMat1 = nullptr;
    ScMatrixRef pMat2 = nullptr;
    double fVal1 = 0.0, fVal2 = 0.0;
    if ( GetStackType() == svMatrix )
        pMat2 = GetMatrix();
    else
        fVal2 = GetDouble();
    if ( GetStackType() == svMatrix )
        pMat1 = GetMatrix();
    else
        fVal1 = GetDouble();
    if (pMat1 && pMat2)
    {
        ScMatrixRef pResMat = lcl_MatrixCalculation( *pMat1, *pMat2, this, MatrixPow);
        if (!pResMat)
            PushNoValue();
        else
            PushMatrix(pResMat);
    }
    else if (pMat1 || pMat2)
    {
        double fVal;
        bool bFlag;
        ScMatrixRef pMat = std::move(pMat1);
        if (!pMat)
        {
            fVal = fVal1;
            pMat = std::move(pMat2);
            bFlag = true;           // double - Matrix
        }
        else
        {
            fVal = fVal2;
            bFlag = false;          // Matrix - double
        }
        SCSIZE nC, nR;
        pMat->GetDimensions(nC, nR);
        ScMatrixRef pResMat = GetNewMat(nC, nR, /*bEmpty*/true);
        if (pResMat)
        {
            pMat->PowOp( bFlag, fVal, *pResMat);
            PushMatrix(pResMat);
        }
        else
            PushIllegalArgument();
    }
    else
    {
        PushDouble( sc::power( fVal1, fVal2));
    }
}

void ScInterpreter::ScSumProduct()
{
    short nParamCount = GetByte();
    if ( !MustHaveParamCountMin( nParamCount, 1) )
        return;

    // XXX NOTE: Excel returns #VALUE! for reference list and 0 (why?) for
    // array of references. We calculate the proper individual arrays if sizes
    // match.

    size_t nInRefList = 0;
    ScMatrixRef pMatLast;
    ScMatrixRef pMat;

    pMatLast = GetMatrix( --nParamCount, nInRefList);
    if (!pMatLast)
    {
        PushIllegalParameter();
        return;
    }

    SCSIZE nC, nCLast, nR, nRLast;
    pMatLast->GetDimensions(nCLast, nRLast);
    std::vector<double> aResArray;
    pMatLast->GetDoubleArray(aResArray);

    while (nParamCount--)
    {
        pMat = GetMatrix( nParamCount, nInRefList);
        if (!pMat)
        {
            PushIllegalParameter();
            return;
        }
        pMat->GetDimensions(nC, nR);
        if (nC != nCLast || nR != nRLast)
        {
            PushNoValue();
            return;
        }

        pMat->MergeDoubleArrayMultiply(aResArray);
    }

    KahanSum fSum = 0.0;
    for( double fPosArray : aResArray )
    {
        FormulaError nErr = GetDoubleErrorValue(fPosArray);
        if (nErr == FormulaError::NONE)
            fSum += fPosArray;
        else if (nErr != FormulaError::ElementNaN)
        {
            // Propagate the first error encountered, ignore "this is not a number" elements.
            PushError(nErr);
            return;
        }
    }

    PushDouble(fSum.get());
}

void ScInterpreter::ScSumX2MY2()
{
    CalculateSumX2MY2SumX2DY2(false);
}
void ScInterpreter::CalculateSumX2MY2SumX2DY2(bool _bSumX2DY2)
{
    if ( !MustHaveParamCount( GetByte(), 2 ) )
        return;

    ScMatrixRef pMat1 = nullptr;
    ScMatrixRef pMat2 = nullptr;
    SCSIZE i, j;
    pMat2 = GetMatrix();
    pMat1 = GetMatrix();
    if (!pMat2 || !pMat1)
    {
        PushIllegalParameter();
        return;
    }
    SCSIZE nC1, nC2;
    SCSIZE nR1, nR2;
    pMat2->GetDimensions(nC2, nR2);
    pMat1->GetDimensions(nC1, nR1);
    if (nC1 != nC2 || nR1 != nR2)
    {
        PushNoValue();
        return;
    }
    double fVal;
    KahanSum fSum = 0.0;
    for (i = 0; i < nC1; i++)
        for (j = 0; j < nR1; j++)
            if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
            {
                fVal = pMat1->GetDouble(i,j);
                fSum += fVal * fVal;
                fVal = pMat2->GetDouble(i,j);
                if ( _bSumX2DY2 )
                    fSum += fVal * fVal;
                else
                    fSum -= fVal * fVal;
            }
    PushDouble(fSum.get());
}

void ScInterpreter::ScSumX2DY2()
{
    CalculateSumX2MY2SumX2DY2(true);
}

void ScInterpreter::ScSumXMY2()
{
    if ( !MustHaveParamCount( GetByte(), 2 ) )
        return;

    ScMatrixRef pMat2 = GetMatrix();
    ScMatrixRef pMat1 = GetMatrix();
    if (!pMat2 || !pMat1)
    {
        PushIllegalParameter();
        return;
    }
    SCSIZE nC1, nC2;
    SCSIZE nR1, nR2;
    pMat2->GetDimensions(nC2, nR2);
    pMat1->GetDimensions(nC1, nR1);
    if (nC1 != nC2 || nR1 != nR2)
    {
        PushNoValue();
        return;
    } // if (nC1 != nC2 || nR1 != nR2)
    ScMatrixRef pResMat = lcl_MatrixCalculation( *pMat1, *pMat2, this, MatrixSub);
    if (!pResMat)
    {
        PushNoValue();
    }
    else
    {
        PushDouble(pResMat->SumSquare(false).maAccumulator.get());
    }
}

void ScInterpreter::ScFrequency()
{
    if ( !MustHaveParamCount( GetByte(), 2 ) )
        return;

    vector<double>  aBinArray;
    vector<tools::Long>    aBinIndexOrder;

    GetSortArray( 1, aBinArray, &aBinIndexOrder, false, false );
    SCSIZE nBinSize = aBinArray.size();
    if (nGlobalError != FormulaError::NONE)
    {
        PushNoValue();
        return;
    }

    vector<double>  aDataArray;
    GetSortArray( 1, aDataArray, nullptr, false, false );
    SCSIZE nDataSize = aDataArray.size();

    if (aDataArray.empty() || nGlobalError != FormulaError::NONE)
    {
        PushNoValue();
        return;
    }
    ScMatrixRef pResMat = GetNewMat(1, nBinSize+1, /*bEmpty*/true);
    if (!pResMat)
    {
        PushIllegalArgument();
        return;
    }

    if (nBinSize != aBinIndexOrder.size())
    {
        PushIllegalArgument();
        return;
    }

    SCSIZE j;
    SCSIZE i = 0;
    for (j = 0; j < nBinSize; ++j)
    {
        SCSIZE nCount = 0;
        while (i < nDataSize && aDataArray[i] <= aBinArray[j])
        {
            ++nCount;
            ++i;
        }
        pResMat->PutDouble(static_cast<double>(nCount), aBinIndexOrder[j]);
    }
    pResMat->PutDouble(static_cast<double>(nDataSize-i), j);
    PushMatrix(pResMat);
}

namespace {

// Helper methods for LINEST/LOGEST and TREND/GROWTH
// All matrices must already exist and have the needed size, no control tests
// done. Those methods, which names start with lcl_T, are adapted to case 3,
// where Y (=observed values) is given as row.
// Remember, ScMatrix matrices are zero based, index access (column,row).

// <A;B> over all elements; uses the matrices as vectors of length M
double lcl_GetSumProduct(const ScMatrixRef& pMatA, const ScMatrixRef& pMatB, SCSIZE nM)
{
    KahanSum fSum = 0.0;
    for (SCSIZE i=0; i<nM; i++)
        fSum += pMatA->GetDouble(i) * pMatB->GetDouble(i);
    return fSum.get();
}

// Special version for use within QR decomposition.
// Euclidean norm of column index C starting in row index R;
// matrix A has count N rows.
double lcl_GetColumnEuclideanNorm(const ScMatrixRef& pMatA, SCSIZE nC, SCSIZE nR, SCSIZE nN)
{
    KahanSum fNorm = 0.0;
    for (SCSIZE row=nR; row<nN; row++)
        fNorm  += (pMatA->GetDouble(nC,row)) * (pMatA->GetDouble(nC,row));
    return sqrt(fNorm.get());
}

// Euclidean norm of row index R starting in column index C;
// matrix A has count N columns.
double lcl_TGetColumnEuclideanNorm(const ScMatrixRef& pMatA, SCSIZE nR, SCSIZE nC, SCSIZE nN)
{
    KahanSum fNorm = 0.0;
    for (SCSIZE col=nC; col<nN; col++)
        fNorm  += (pMatA->GetDouble(col,nR)) * (pMatA->GetDouble(col,nR));
    return sqrt(fNorm.get());
}

// Special version for use within QR decomposition.
// Maximum norm of column index C starting in row index R;
// matrix A has count N rows.
double lcl_GetColumnMaximumNorm(const ScMatrixRef& pMatA, SCSIZE nC, SCSIZE nR, SCSIZE nN)
{
    double fNorm = 0.0;
    for (SCSIZE row=nR; row<nN; row++)
    {
        double fVal = fabs(pMatA->GetDouble(nC,row));
        if (fNorm < fVal)
            fNorm = fVal;
    }
    return fNorm;
}

// Maximum norm of row index R starting in col index C;
// matrix A has count N columns.
double lcl_TGetColumnMaximumNorm(const ScMatrixRef& pMatA, SCSIZE nR, SCSIZE nC, SCSIZE nN)
{
    double fNorm = 0.0;
    for (SCSIZE col=nC; col<nN; col++)
    {
        double fVal = fabs(pMatA->GetDouble(col,nR));
        if (fNorm < fVal)
            fNorm = fVal;
    }
    return fNorm;
}

// Special version for use within QR decomposition.
// <A(Ca);B(Cb)> starting in row index R;
// Ca and Cb are indices of columns, matrices A and B have count N rows.
double lcl_GetColumnSumProduct(const ScMatrixRef& pMatA, SCSIZE nCa,
                               const ScMatrixRef& pMatB, SCSIZE nCb, SCSIZE nR, SCSIZE nN)
{
    KahanSum fResult = 0.0;
    for (SCSIZE row=nR; row<nN; row++)
        fResult += pMatA->GetDouble(nCa,row) * pMatB->GetDouble(nCb,row);
    return fResult.get();
}

// <A(Ra);B(Rb)> starting in column index C;
// Ra and Rb are indices of rows, matrices A and B have count N columns.
double lcl_TGetColumnSumProduct(const ScMatrixRef& pMatA, SCSIZE nRa,
                                const ScMatrixRef& pMatB, SCSIZE nRb, SCSIZE nC, SCSIZE nN)
{
    KahanSum fResult = 0.0;
    for (SCSIZE col=nC; col<nN; col++)
        fResult += pMatA->GetDouble(col,nRa) * pMatB->GetDouble(col,nRb);
    return fResult.get();
}

// no mathematical signum, but used to switch between adding and subtracting
double lcl_GetSign(double fValue)
{
    return (fValue >= 0.0 ? 1.0 : -1.0 );
}

/* Calculates a QR decomposition with Householder reflection.
 * For each NxK matrix A exists a decomposition A=Q*R with an orthogonal
 * NxN matrix Q and a NxK matrix R.
 * Q=H1*H2*...*Hk with Householder matrices H. Such a householder matrix can
 * be build from a vector u by H=I-(2/u'u)*(u u'). This vectors u are returned
 * in the columns of matrix A, overwriting the old content.
 * The matrix R has a quadric upper part KxK with values in the upper right
 * triangle and zeros in all other elements. Here the diagonal elements of R
 * are stored in the vector R and the other upper right elements in the upper
 * right of the matrix A.
 * The function returns false, if calculation breaks. But because of round-off
 * errors singularity is often not detected.
 */
bool lcl_CalculateQRdecomposition(const ScMatrixRef& pMatA,
                                  ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN)
{
    // ScMatrix matrices are zero based, index access (column,row)
    for (SCSIZE col = 0; col <nK; col++)
    {
        // calculate vector u of the householder transformation
        const double fScale = lcl_GetColumnMaximumNorm(pMatA, col, col, nN);
        if (fScale == 0.0)
        {
            // A is singular
            return false;
        }
        for (SCSIZE row = col; row <nN; row++)
            pMatA->PutDouble( pMatA->GetDouble(col,row)/fScale, col, row);

        const double fEuclid = lcl_GetColumnEuclideanNorm(pMatA, col, col, nN);
        const double fFactor = 1.0/fEuclid/(fEuclid + fabs(pMatA->GetDouble(col,col)));
        const double fSignum = lcl_GetSign(pMatA->GetDouble(col,col));
        pMatA->PutDouble( pMatA->GetDouble(col,col) + fSignum*fEuclid, col,col);
        pVecR[col] = -fSignum * fScale * fEuclid;

        // apply Householder transformation to A
        for (SCSIZE c=col+1; c<nK; c++)
        {
            const double fSum =lcl_GetColumnSumProduct(pMatA, col, pMatA, c, col, nN);
            for (SCSIZE row = col; row <nN; row++)
                pMatA->PutDouble( pMatA->GetDouble(c,row) - fSum * fFactor * pMatA->GetDouble(col,row), c, row);
        }
    }
    return true;
}

// same with transposed matrix A, N is count of columns, K count of rows
bool lcl_TCalculateQRdecomposition(const ScMatrixRef& pMatA,
                                   ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN)
{
    double fSum ;
    // ScMatrix matrices are zero based, index access (column,row)
    for (SCSIZE row = 0; row <nK; row++)
    {
        // calculate vector u of the householder transformation
        const double fScale = lcl_TGetColumnMaximumNorm(pMatA, row, row, nN);
        if (fScale == 0.0)
        {
            // A is singular
            return false;
        }
        for (SCSIZE col = row; col <nN; col++)
            pMatA->PutDouble( pMatA->GetDouble(col,row)/fScale, col, row);

        const double fEuclid = lcl_TGetColumnEuclideanNorm(pMatA, row, row, nN);
        const double fFactor = 1.0/fEuclid/(fEuclid + fabs(pMatA->GetDouble(row,row)));
        const double fSignum = lcl_GetSign(pMatA->GetDouble(row,row));
        pMatA->PutDouble( pMatA->GetDouble(row,row) + fSignum*fEuclid, row,row);
        pVecR[row] = -fSignum * fScale * fEuclid;

        // apply Householder transformation to A
        for (SCSIZE r=row+1; r<nK; r++)
        {
            fSum =lcl_TGetColumnSumProduct(pMatA, row, pMatA, r, row, nN);
            for (SCSIZE col = row; col <nN; col++)
                pMatA->PutDouble(
                    pMatA->GetDouble(col,r) - fSum * fFactor * pMatA->GetDouble(col,row), col, r);
        }
    }
    return true;
}

/* Applies a Householder transformation to a column vector Y with is given as
 * Nx1 Matrix. The vector u, from which the Householder transformation is built,
 * is the column part in matrix A, with column index C, starting with row
 * index C. A is the result of the QR decomposition as obtained from
 * lcl_CalculateQRdecomposition.
 */
void lcl_ApplyHouseholderTransformation(const ScMatrixRef& pMatA, SCSIZE nC,
--> --------------------

--> maximum size reached

--> --------------------