TDF Git
Code Review
gerrit.libreoffice.org / core / refs/heads/master / . / sc / source / core / tool / interpr5.cxx
blob: 2f4fe0f19fa1ecd0292642b44f5b44b7168d1c51 [file] [log] [blame]
/* -*- 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 = 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& 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 = pMat1;
if (!pMat)
{
fVal = fVal1;
pMat = 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 = pMat1;
if (!pMat)
{
sStr = sStr1;
pMat = 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 = pMat1;
if (!pMat)
{
fVal = fVal1;
pMat = 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 = pMat1;
if (!pMat)
{
fVal = fVal1;
pMat = 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 = pMat1;
if (!pMat)
{
fVal = fVal1;
pMat = 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,
const ScMatrixRef& pMatY, SCSIZE nN)
{
// ScMatrix matrices are zero based, index access (column,row)
double fDenominator = lcl_GetColumnSumProduct(pMatA, nC, pMatA, nC, nC, nN);
double fNumerator = lcl_GetColumnSumProduct(pMatA, nC, pMatY, 0, nC, nN);
double fFactor = 2.0 * (fNumerator/fDenominator);
for (SCSIZE row = nC; row < nN; row++)
pMatY->PutDouble(
pMatY->GetDouble(row) - fFactor * pMatA->GetDouble(nC,row), row);
}
// Same with transposed matrices A and Y.
void lcl_TApplyHouseholderTransformation(const ScMatrixRef& pMatA, SCSIZE nR,
const ScMatrixRef& pMatY, SCSIZE nN)
{
// ScMatrix matrices are zero based, index access (column,row)
double fDenominator = lcl_TGetColumnSumProduct(pMatA, nR, pMatA, nR, nR, nN);
double fNumerator = lcl_TGetColumnSumProduct(pMatA, nR, pMatY, 0, nR, nN);
double fFactor = 2.0 * (fNumerator/fDenominator);
for (SCSIZE col = nR; col < nN; col++)
pMatY->PutDouble(
pMatY->GetDouble(col) - fFactor * pMatA->GetDouble(col,nR), col);
}
/* Solve for X in R*X=S using back substitution. The solution X overwrites S.
* Uses R from the result of the QR decomposition of a NxK matrix A.
* S is a column vector given as matrix, with at least elements on index
* 0 to K-1; elements on index>=K are ignored. Vector R must not have zero
* elements, no check is done.
*/
void lcl_SolveWithUpperRightTriangle(const ScMatrixRef& pMatA,
::std::vector< double>& pVecR, const ScMatrixRef& pMatS,
SCSIZE nK, bool bIsTransposed)
{
// ScMatrix matrices are zero based, index access (column,row)
SCSIZE row;
// SCSIZE is never negative, therefore test with rowp1=row+1
for (SCSIZE rowp1 = nK; rowp1>0; rowp1--)
{
row = rowp1-1;
KahanSum fSum = pMatS->GetDouble(row);
for (SCSIZE col = rowp1; col<nK ; col++)
if (bIsTransposed)
fSum -= pMatA->GetDouble(row,col) * pMatS->GetDouble(col);
else
fSum -= pMatA->GetDouble(col,row) * pMatS->GetDouble(col);
pMatS->PutDouble( fSum.get() / pVecR[row] , row);
}
}
/* Solve for X in R' * X= T using forward substitution. The solution X
* overwrites T. Uses R from the result of the QR decomposition of a NxK
* matrix A. T is a column vectors given as matrix, with at least elements on
* index 0 to K-1; elements on index>=K are ignored. Vector R must not have
* zero elements, no check is done.
*/
void lcl_SolveWithLowerLeftTriangle(const ScMatrixRef& pMatA,
::std::vector< double>& pVecR, const ScMatrixRef& pMatT,
SCSIZE nK, bool bIsTransposed)
{
// ScMatrix matrices are zero based, index access (column,row)
for (SCSIZE row = 0; row < nK; row++)
{
KahanSum fSum = pMatT -> GetDouble(row);
for (SCSIZE col=0; col < row; col++)
{
if (bIsTransposed)
fSum -= pMatA->GetDouble(col,row) * pMatT->GetDouble(col);
else
fSum -= pMatA->GetDouble(row,col) * pMatT->GetDouble(col);
}
pMatT->PutDouble( fSum.get() / pVecR[row] , row);
}
}
/* Calculates Z = R * B
* R is given in matrix A and vector VecR as obtained from the QR
* decomposition in lcl_CalculateQRdecomposition. B and Z are column vectors
* given as matrix with at least index 0 to K-1; elements on index>=K are
* not used.
*/
void lcl_ApplyUpperRightTriangle(const ScMatrixRef& pMatA,
::std::vector< double>& pVecR, const ScMatrixRef& pMatB,
const ScMatrixRef& pMatZ, SCSIZE nK, bool bIsTransposed)
{
// ScMatrix matrices are zero based, index access (column,row)
for (SCSIZE row = 0; row < nK; row++)
{
KahanSum fSum = pVecR[row] * pMatB->GetDouble(row);
for (SCSIZE col = row+1; col < nK; col++)
if (bIsTransposed)
fSum += pMatA->GetDouble(row,col) * pMatB->GetDouble(col);
else
fSum += pMatA->GetDouble(col,row) * pMatB->GetDouble(col);
pMatZ->PutDouble( fSum.get(), row);
}
}
double lcl_GetMeanOverAll(const ScMatrixRef& pMat, SCSIZE nN)
{
KahanSum fSum = 0.0;
for (SCSIZE i=0 ; i<nN; i++)
fSum += pMat->GetDouble(i);
return fSum.get()/static_cast<double>(nN);
}
// Calculates means of the columns of matrix X. X is a RxC matrix;
// ResMat is a 1xC matrix (=row).
void lcl_CalculateColumnMeans(const ScMatrixRef& pX, const ScMatrixRef& pResMat,
SCSIZE nC, SCSIZE nR)
{
for (SCSIZE i=0; i < nC; i++)
{
KahanSum fSum =0.0;
for (SCSIZE k=0; k < nR; k++)
fSum += pX->GetDouble(i,k); // GetDouble(Column,Row)
pResMat ->PutDouble( fSum.get()/static_cast<double>(nR),i);
}
}
// Calculates means of the rows of matrix X. X is a RxC matrix;
// ResMat is a Rx1 matrix (=column).
void lcl_CalculateRowMeans(const ScMatrixRef& pX, const ScMatrixRef& pResMat,
SCSIZE nC, SCSIZE nR)
{
for (SCSIZE k=0; k < nR; k++)
{
KahanSum fSum = 0.0;
for (SCSIZE i=0; i < nC; i++)
fSum += pX->GetDouble(i,k); // GetDouble(Column,Row)
pResMat ->PutDouble( fSum.get()/static_cast<double>(nC),k);
}
}
void lcl_CalculateColumnsDelta(const ScMatrixRef& pMat, const ScMatrixRef& pColumnMeans,
SCSIZE nC, SCSIZE nR)
{
for (SCSIZE i = 0; i < nC; i++)
for (SCSIZE k = 0; k < nR; k++)
pMat->PutDouble( ::rtl::math::approxSub
(pMat->GetDouble(i,k) , pColumnMeans->GetDouble(i) ) , i, k);
}
void lcl_CalculateRowsDelta(const ScMatrixRef& pMat, const ScMatrixRef& pRowMeans,
SCSIZE nC, SCSIZE nR)
{
for (SCSIZE k = 0; k < nR; k++)
for (SCSIZE i = 0; i < nC; i++)
pMat->PutDouble( ::rtl::math::approxSub
( pMat->GetDouble(i,k) , pRowMeans->GetDouble(k) ) , i, k);
}
// Case1 = simple regression
// MatX = X - MeanX, MatY = Y - MeanY, y - haty = (y - MeanY) - (haty - MeanY)
// = (y-MeanY)-((slope*x+a)-(slope*MeanX+a)) = (y-MeanY)-slope*(x-MeanX)
double lcl_GetSSresid(const ScMatrixRef& pMatX, const ScMatrixRef& pMatY, double fSlope,
SCSIZE nN)
{
KahanSum fSum = 0.0;
for (SCSIZE i=0; i<nN; i++)
{
const double fTemp = pMatY->GetDouble(i) - fSlope * pMatX->GetDouble(i);
fSum += fTemp * fTemp;
}
return fSum.get();
}
}
// Fill default values in matrix X, transform Y to log(Y) in case LOGEST|GROWTH,
// determine sizes of matrices X and Y, determine kind of regression, clone
// Y in case LOGEST|GROWTH, if constant.
bool ScInterpreter::CheckMatrix(bool _bLOG, sal_uInt8& nCase, SCSIZE& nCX,
SCSIZE& nCY, SCSIZE& nRX, SCSIZE& nRY, SCSIZE& M,
SCSIZE& N, ScMatrixRef& pMatX, ScMatrixRef& pMatY)
{
nCX = 0;
nCY = 0;
nRX = 0;
nRY = 0;
M = 0;
N = 0;
pMatY->GetDimensions(nCY, nRY);
const SCSIZE nCountY = nCY * nRY;
for ( SCSIZE i = 0; i < nCountY; i++ )
{
if (!pMatY->IsValue(i))
{
PushIllegalArgument();
return false;
}
}
if ( _bLOG )
{
ScMatrixRef pNewY = pMatY->CloneIfConst();
for (SCSIZE nElem = 0; nElem < nCountY; nElem++)
{
const double fVal = pNewY->GetDouble(nElem);
if (fVal <= 0.0)
{
PushIllegalArgument();
return false;
}
else
pNewY->PutDouble(log(fVal), nElem);
}
pMatY = pNewY;
}
if (pMatX)
{
pMatX->GetDimensions(nCX, nRX);
const SCSIZE nCountX = nCX * nRX;
for ( SCSIZE i = 0; i < nCountX; i++ )
if (!pMatX->IsValue(i))
{
PushIllegalArgument();
return false;
}
if (nCX == nCY && nRX == nRY)
{
nCase = 1; // simple regression
M = 1;
N = nCountY;
}
else if (nCY != 1 && nRY != 1)
{
PushIllegalArgument();
return false;
}
else if (nCY == 1)
{
if (nRX != nRY)
{
PushIllegalArgument();
return false;
}
else
{
nCase = 2; // Y is column
N = nRY;
M = nCX;
}
}
else if (nCX != nCY)
{
PushIllegalArgument();
return false;
}
else
{
nCase = 3; // Y is row
N = nCY;
M = nRX;
}
}
else
{
pMatX = GetNewMat(nCY, nRY, /*bEmpty*/true);
nCX = nCY;
nRX = nRY;
if (!pMatX)
{
PushIllegalArgument();
return false;
}
for ( SCSIZE i = 1; i <= nCountY; i++ )
pMatX->PutDouble(static_cast<double>(i), i-1);
nCase = 1;
N = nCountY;
M = 1;
}
return true;
}
// LINEST
void ScInterpreter::ScLinest()
{
CalculateRGPRKP(false);
}
// LOGEST
void ScInterpreter::ScLogest()
{
CalculateRGPRKP(true);
}
void ScInterpreter::CalculateRGPRKP(bool _bRKP)
{
sal_uInt8 nParamCount = GetByte();
if (!MustHaveParamCount( nParamCount, 1, 4 ))
return;
bool bConstant, bStats;
// optional forth parameter
if (nParamCount == 4)
bStats = GetBool();
else
bStats = false;
// The third parameter may not be missing in ODF, if the forth parameter
// is present. But Excel allows it with default true, we too.
if (nParamCount >= 3)
{
if (IsMissing())
{
Pop();
bConstant = true;
// PushIllegalParameter(); if ODF behavior is desired
// return;
}
else
bConstant = GetBool();
}
else
bConstant = true;
ScMatrixRef pMatX;
if (nParamCount >= 2)
{
if (IsMissing())
{ //In ODF1.2 empty second parameter (which is two ;; ) is allowed
Pop();
pMatX = nullptr;
}
else
{
pMatX = GetMatrix();
}
}
else
pMatX = nullptr;
ScMatrixRef pMatY = GetMatrix();
if (!pMatY)
{
PushIllegalParameter();
return;
}
// 1 = simple; 2 = multiple with Y as column; 3 = multiple with Y as row
sal_uInt8 nCase;
SCSIZE nCX, nCY; // number of columns
SCSIZE nRX, nRY; //number of rows
SCSIZE K = 0, N = 0; // K=number of variables X, N=number of data samples
if (!CheckMatrix(_bRKP,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY))
{
PushIllegalParameter();
return;
}
// Enough data samples?
if ((bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1))
{
PushIllegalParameter();
return;
}
ScMatrixRef pResMat;
if (bStats)
pResMat = GetNewMat(K+1,5, /*bEmpty*/true);
else
pResMat = GetNewMat(K+1,1, /*bEmpty*/true);
if (!pResMat)
{
PushError(FormulaError::CodeOverflow);
return;
}
// Fill unused cells in pResMat; order (column,row)
if (bStats)
{
for (SCSIZE i=2; i<K+1; i++)
{
pResMat->PutError( FormulaError::NotAvailable, i, 2);
pResMat->PutError( FormulaError::NotAvailable, i, 3);
pResMat->PutError( FormulaError::NotAvailable, i, 4);
}
}
// Uses sum(x-MeanX)^2 and not [sum x^2]-N * MeanX^2 in case bConstant.
// Clone constant matrices, so that Mat = Mat - Mean is possible.
double fMeanY = 0.0;
if (bConstant)
{
ScMatrixRef pNewX = pMatX->CloneIfConst();
ScMatrixRef pNewY = pMatY->CloneIfConst();
if (!pNewX || !pNewY)
{
PushError(FormulaError::CodeOverflow);
return;
}
pMatX = pNewX;
pMatY = pNewY;
// DeltaY is possible here; DeltaX depends on nCase, so later
fMeanY = lcl_GetMeanOverAll(pMatY, N);
for (SCSIZE i=0; i<N; i++)
{
pMatY->PutDouble( ::rtl::math::approxSub(pMatY->GetDouble(i),fMeanY), i );
}
}
if (nCase==1)
{
// calculate simple regression
double fMeanX = 0.0;
if (bConstant)
{ // Mat = Mat - Mean
fMeanX = lcl_GetMeanOverAll(pMatX, N);
for (SCSIZE i=0; i<N; i++)
{
pMatX->PutDouble( ::rtl::math::approxSub(pMatX->GetDouble(i),fMeanX), i );
}
}
double fSumXY = lcl_GetSumProduct(pMatX,pMatY,N);
double fSumX2 = lcl_GetSumProduct(pMatX,pMatX,N);
if (fSumX2==0.0)
{
PushNoValue(); // all x-values are identical
return;
}
double fSlope = fSumXY / fSumX2;
double fIntercept = 0.0;
if (bConstant)
fIntercept = fMeanY - fSlope * fMeanX;
pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, 1, 0); //order (column,row)
pResMat->PutDouble(_bRKP ? exp(fSlope) : fSlope, 0, 0);
if (bStats)
{
double fSSreg = fSlope * fSlope * fSumX2;
pResMat->PutDouble(fSSreg, 0, 4);
double fDegreesFreedom =static_cast<double>( bConstant ? N-2 : N-1 );
pResMat->PutDouble(fDegreesFreedom, 1, 3);
double fSSresid = lcl_GetSSresid(pMatX,pMatY,fSlope,N);
pResMat->PutDouble(fSSresid, 1, 4);
if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0)
{ // exact fit; test SSreg too, because SSresid might be
// unequal zero due to round of errors
pResMat->PutDouble(0.0, 1, 4); // SSresid
pResMat->PutError( FormulaError::NotAvailable, 0, 3); // F
pResMat->PutDouble(0.0, 1, 2); // RMSE
pResMat->PutDouble(0.0, 0, 1); // SigmaSlope
if (bConstant)
pResMat->PutDouble(0.0, 1, 1); //SigmaIntercept
else
pResMat->PutError( FormulaError::NotAvailable, 1, 1);
pResMat->PutDouble(1.0, 0, 2); // R^2
}
else
{
double fFstatistic = (fSSreg / static_cast<double>(K))
/ (fSSresid / fDegreesFreedom);
pResMat->PutDouble(fFstatistic, 0, 3);
// standard error of estimate
double fRMSE = sqrt(fSSresid / fDegreesFreedom);
pResMat->PutDouble(fRMSE, 1, 2);
double fSigmaSlope = fRMSE / sqrt(fSumX2);
pResMat->PutDouble(fSigmaSlope, 0, 1);
if (bConstant)
{
double fSigmaIntercept = fRMSE
* sqrt(fMeanX*fMeanX/fSumX2 + 1.0/static_cast<double>(N));
pResMat->PutDouble(fSigmaIntercept, 1, 1);
}
else
{
pResMat->PutError( FormulaError::NotAvailable, 1, 1);
}
double fR2 = fSSreg / (fSSreg + fSSresid);
pResMat->PutDouble(fR2, 0, 2);
}
}
PushMatrix(pResMat);
}
else // calculate multiple regression;
{
// Uses a QR decomposition X = QR. The solution B = (X'X)^(-1) * X' * Y
// becomes B = R^(-1) * Q' * Y
if (nCase ==2) // Y is column
{
::std::vector< double> aVecR(N); // for QR decomposition
// Enough memory for needed matrices?
ScMatrixRef pMeans = GetNewMat(K, 1, /*bEmpty*/true); // mean of each column
ScMatrixRef pMatZ; // for Q' * Y , inter alia
if (bStats)
pMatZ = pMatY->Clone(); // Y is used in statistic, keep it
else
pMatZ = pMatY; // Y can be overwritten
ScMatrixRef pSlopes = GetNewMat(1,K, /*bEmpty*/true); // from b1 to bK
if (!pMeans || !pMatZ || !pSlopes)
{
PushError(FormulaError::CodeOverflow);
return;
}
if (bConstant)
{
lcl_CalculateColumnMeans(pMatX, pMeans, K, N);
lcl_CalculateColumnsDelta(pMatX, pMeans, K, N);
}
if (!lcl_CalculateQRdecomposition(pMatX, aVecR, K, N))
{
PushNoValue();
return;
}
// Later on we will divide by elements of aVecR, so make sure
// that they aren't zero.
bool bIsSingular=false;
for (SCSIZE row=0; row < K && !bIsSingular; row++)
bIsSingular = aVecR[row] == 0.0;
if (bIsSingular)
{
PushNoValue();
return;
}
// Z = Q' Y;
for (SCSIZE col = 0; col < K; col++)
{
lcl_ApplyHouseholderTransformation(pMatX, col, pMatZ, N);
}
// B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
// result Z should have zeros for index>=K; if not, ignore values
for (SCSIZE col = 0; col < K ; col++)
{
pSlopes->PutDouble( pMatZ->GetDouble(col), col);
}
lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, false);
double fIntercept = 0.0;
if (bConstant)
fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K);
// Fill first line in result matrix
pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, K, 0 );
for (SCSIZE i = 0; i < K; i++)
pResMat->PutDouble(_bRKP ? exp(pSlopes->GetDouble(i))
: pSlopes->GetDouble(i) , K-1-i, 0);
if (bStats)
{
double fSSreg = 0.0;
double fSSresid = 0.0;
// re-use memory of Z;
pMatZ->FillDouble(0.0, 0, 0, 0, N-1);
// Z = R * Slopes
lcl_ApplyUpperRightTriangle(pMatX, aVecR, pSlopes, pMatZ, K, false);
// Z = Q * Z, that is Q * R * Slopes = X * Slopes
for (SCSIZE colp1 = K; colp1 > 0; colp1--)
{
lcl_ApplyHouseholderTransformation(pMatX, colp1-1, pMatZ,N);
}
fSSreg =lcl_GetSumProduct(pMatZ, pMatZ, N);
// re-use Y for residuals, Y = Y-Z
for (SCSIZE row = 0; row < N; row++)
pMatY->PutDouble(pMatY->GetDouble(row) - pMatZ->GetDouble(row), row);
fSSresid = lcl_GetSumProduct(pMatY, pMatY, N);
pResMat->PutDouble(fSSreg, 0, 4);
pResMat->PutDouble(fSSresid, 1, 4);
double fDegreesFreedom =static_cast<double>( bConstant ? N-K-1 : N-K );
pResMat->PutDouble(fDegreesFreedom, 1, 3);
if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0)
{ // exact fit; incl. observed values Y are identical
pResMat->PutDouble(0.0, 1, 4); // SSresid
// F = (SSreg/K) / (SSresid/df) = #DIV/0!
pResMat->PutError( FormulaError::NotAvailable, 0, 3); // F
// RMSE = sqrt(SSresid / df) = sqrt(0 / df) = 0
pResMat->PutDouble(0.0, 1, 2); // RMSE
// SigmaSlope[i] = RMSE * sqrt(matrix[i,i]) = 0 * sqrt(...) = 0
for (SCSIZE i=0; i<K; i++)
pResMat->PutDouble(0.0, K-1-i, 1);
// SigmaIntercept = RMSE * sqrt(...) = 0
if (bConstant)
pResMat->PutDouble(0.0, K, 1); //SigmaIntercept
else
pResMat->PutError( FormulaError::NotAvailable, K, 1);
// R^2 = SSreg / (SSreg + SSresid) = 1.0
pResMat->PutDouble(1.0, 0, 2); // R^2
}
else
{
double fFstatistic = (fSSreg / static_cast<double>(K))
/ (fSSresid / fDegreesFreedom);
pResMat->PutDouble(fFstatistic, 0, 3);
// standard error of estimate = root mean SSE
double fRMSE = sqrt(fSSresid / fDegreesFreedom);
pResMat->PutDouble(fRMSE, 1, 2);
// standard error of slopes
// = RMSE * sqrt(diagonal element of (R' R)^(-1) )
// standard error of intercept
// = RMSE * sqrt( Xmean * (R' R)^(-1) * Xmean' + 1/N)
// (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as
// a whole matrix, but iterate over unit vectors.
KahanSum aSigmaIntercept = 0.0;
double fPart; // for Xmean * single column of (R' R)^(-1)
for (SCSIZE col = 0; col < K; col++)
{
//re-use memory of MatZ
pMatZ->FillDouble(0.0,0,0,0,K-1); // Z = unit vector e
pMatZ->PutDouble(1.0, col);
//Solve R' * Z = e
lcl_SolveWithLowerLeftTriangle(pMatX, aVecR, pMatZ, K, false);
// Solve R * Znew = Zold
lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pMatZ, K, false);
// now Z is column col in (R' R)^(-1)
double fSigmaSlope = fRMSE * sqrt(pMatZ->GetDouble(col));
pResMat->PutDouble(fSigmaSlope, K-1-col, 1);
// (R' R) ^(-1) is symmetric
if (bConstant)
{
fPart = lcl_GetSumProduct(pMeans, pMatZ, K);
aSigmaIntercept += fPart * pMeans->GetDouble(col);
}
}
if (bConstant)
{
double fSigmaIntercept = fRMSE
* sqrt( (aSigmaIntercept + 1.0 / static_cast<double>(N) ).get() );
pResMat->PutDouble(fSigmaIntercept, K, 1);
}
else
{
pResMat->PutError( FormulaError::NotAvailable, K, 1);
}
double fR2 = fSSreg / (fSSreg + fSSresid);
pResMat->PutDouble(fR2, 0, 2);
}
}
PushMatrix(pResMat);
}
else // nCase == 3, Y is row, all matrices are transposed
{
::std::vector< double> aVecR(N); // for QR decomposition
// Enough memory for needed matrices?
ScMatrixRef pMeans = GetNewMat(1, K, /*bEmpty*/true); // mean of each row
ScMatrixRef pMatZ; // for Q' * Y , inter alia
if (bStats)
pMatZ = pMatY->Clone(); // Y is used in statistic, keep it
else
pMatZ = pMatY; // Y can be overwritten
ScMatrixRef pSlopes = GetNewMat(K,1, /*bEmpty*/true); // from b1 to bK
if (!pMeans || !pMatZ || !pSlopes)
{
PushError(FormulaError::CodeOverflow);
return;
}
if (bConstant)
{
lcl_CalculateRowMeans(pMatX, pMeans, N, K);
lcl_CalculateRowsDelta(pMatX, pMeans, N, K);
}
if (!lcl_TCalculateQRdecomposition(pMatX, aVecR, K, N))
{
PushNoValue();
return;
}
// Later on we will divide by elements of aVecR, so make sure
// that they aren't zero.
bool bIsSingular=false;
for (SCSIZE row=0; row < K && !bIsSingular; row++)
bIsSingular = aVecR[row] == 0.0;
if (bIsSingular)
{
PushNoValue();
return;
}
// Z = Q' Y
for (SCSIZE row = 0; row < K; row++)
{
lcl_TApplyHouseholderTransformation(pMatX, row, pMatZ, N);
}
// B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
// result Z should have zeros for index>=K; if not, ignore values
for (SCSIZE col = 0; col < K ; col++)
{
pSlopes->PutDouble( pMatZ->GetDouble(col), col);
}
lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, true);
double fIntercept = 0.0;
if (bConstant)
fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K);
// Fill first line in result matrix
pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, K, 0 );
for (SCSIZE i = 0; i < K; i++)
pResMat->PutDouble(_bRKP ? exp(pSlopes->GetDouble(i))
: pSlopes->GetDouble(i) , K-1-i, 0);
if (bStats)
{
double fSSreg = 0.0;
double fSSresid = 0.0;
// re-use memory of Z;
pMatZ->FillDouble(0.0, 0, 0, N-1, 0);
// Z = R * Slopes
lcl_ApplyUpperRightTriangle(pMatX, aVecR, pSlopes, pMatZ, K, true);
// Z = Q * Z, that is Q * R * Slopes = X * Slopes
for (SCSIZE rowp1 = K; rowp1 > 0; rowp1--)
{
lcl_TApplyHouseholderTransformation(pMatX, rowp1-1, pMatZ,N);
}
fSSreg =lcl_GetSumProduct(pMatZ, pMatZ, N);
// re-use Y for residuals, Y = Y-Z
for (SCSIZE col = 0; col < N; col++)
pMatY->PutDouble(pMatY->GetDouble(col) - pMatZ->GetDouble(col), col);
fSSresid = lcl_GetSumProduct(pMatY, pMatY, N);
pResMat->PutDouble(fSSreg, 0, 4);
pResMat->PutDouble(fSSresid, 1, 4);
double fDegreesFreedom =static_cast<double>( bConstant ? N-K-1 : N-K );
pResMat->PutDouble(fDegreesFreedom, 1, 3);
if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0)
{ // exact fit; incl. case observed values Y are identical
pResMat->PutDouble(0.0, 1, 4); // SSresid
// F = (SSreg/K) / (SSresid/df) = #DIV/0!
pResMat->PutError( FormulaError::NotAvailable, 0, 3); // F
// RMSE = sqrt(SSresid / df) = sqrt(0 / df) = 0
pResMat->PutDouble(0.0, 1, 2); // RMSE
// SigmaSlope[i] = RMSE * sqrt(matrix[i,i]) = 0 * sqrt(...) = 0
for (SCSIZE i=0; i<K; i++)
pResMat->PutDouble(0.0, K-1-i, 1);
// SigmaIntercept = RMSE * sqrt(...) = 0
if (bConstant)
pResMat->PutDouble(0.0, K, 1); //SigmaIntercept
else
pResMat->PutError( FormulaError::NotAvailable, K, 1);
// R^2 = SSreg / (SSreg + SSresid) = 1.0
pResMat->PutDouble(1.0, 0, 2); // R^2
}
else
{
double fFstatistic = (fSSreg / static_cast<double>(K))
/ (fSSresid / fDegreesFreedom);
pResMat->PutDouble(fFstatistic, 0, 3);
// standard error of estimate = root mean SSE
double fRMSE = sqrt(fSSresid / fDegreesFreedom);
pResMat->PutDouble(fRMSE, 1, 2);
// standard error of slopes
// = RMSE * sqrt(diagonal element of (R' R)^(-1) )
// standard error of intercept
// = RMSE * sqrt( Xmean * (R' R)^(-1) * Xmean' + 1/N)
// (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as
// a whole matrix, but iterate over unit vectors.
// (R' R) ^(-1) is symmetric
KahanSum aSigmaIntercept = 0.0;
double fPart; // for Xmean * single col of (R' R)^(-1)
for (SCSIZE row = 0; row < K; row++)
{
//re-use memory of MatZ
pMatZ->FillDouble(0.0,0,0,K-1,0); // Z = unit vector e
pMatZ->PutDouble(1.0, row);
//Solve R' * Z = e
lcl_SolveWithLowerLeftTriangle(pMatX, aVecR, pMatZ, K, true);
// Solve R * Znew = Zold
lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pMatZ, K, true);
// now Z is column col in (R' R)^(-1)
double fSigmaSlope = fRMSE * sqrt(pMatZ->GetDouble(row));
pResMat->PutDouble(fSigmaSlope, K-1-row, 1);
if (bConstant)
{
fPart = lcl_GetSumProduct(pMeans, pMatZ, K);
aSigmaIntercept += fPart * pMeans->GetDouble(row);
}
}
if (bConstant)
{
double fSigmaIntercept = fRMSE
* sqrt( (aSigmaIntercept + 1.0 / static_cast<double>(N) ).get() );
pResMat->PutDouble(fSigmaIntercept, K, 1);
}
else
{
pResMat->PutError( FormulaError::NotAvailable, K, 1);
}
double fR2 = fSSreg / (fSSreg + fSSresid);
pResMat->PutDouble(fR2, 0, 2);
}
}
PushMatrix(pResMat);
}
}
}
void ScInterpreter::ScTrend()
{
CalculateTrendGrowth(false);
}
void ScInterpreter::ScGrowth()
{
CalculateTrendGrowth(true);
}
void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
{
sal_uInt8 nParamCount = GetByte();
if (!MustHaveParamCount( nParamCount, 1, 4 ))
return;
// optional forth parameter
bool bConstant;
if (nParamCount == 4)
bConstant = GetBool();
else
bConstant = true;
// The third parameter may be missing in ODF, although the forth parameter
// is present. Default values depend on data not yet read.
ScMatrixRef pMatNewX;
if (nParamCount >= 3)
{
if (IsMissing())
{
Pop();
pMatNewX = nullptr;
}
else
pMatNewX = GetMatrix();
}
else
pMatNewX = nullptr;
//In ODF1.2 empty second parameter (which is two ;; ) is allowed
//Defaults will be set in CheckMatrix
ScMatrixRef pMatX;
if (nParamCount >= 2)
{
if (IsMissing())
{
Pop();
pMatX = nullptr;
}
else
{
pMatX = GetMatrix();
}
}
else
pMatX = nullptr;
ScMatrixRef pMatY = GetMatrix();
if (!pMatY)
{
PushIllegalParameter();
return;
}
// 1 = simple; 2 = multiple with Y as column; 3 = multiple with Y as row
sal_uInt8 nCase;
SCSIZE nCX, nCY; // number of columns
SCSIZE nRX, nRY; //number of rows
SCSIZE K = 0, N = 0; // K=number of variables X, N=number of data samples
if (!CheckMatrix(_bGrowth,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY))
{
PushIllegalParameter();
return;
}
// Enough data samples?
if ((bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1))
{
PushIllegalParameter();
return;
}
// Set default pMatNewX if necessary
SCSIZE nCXN, nRXN;
SCSIZE nCountXN;
if (!pMatNewX)
{
nCXN = nCX;
nRXN = nRX;
nCountXN = nCXN * nRXN;
pMatNewX = pMatX->Clone(); // pMatX will be changed to X-meanX
}
else
{
pMatNewX->GetDimensions(nCXN, nRXN);
if ((nCase == 2 && K != nCXN) || (nCase == 3 && K != nRXN))
{
PushIllegalArgument();
return;
}
nCountXN = nCXN * nRXN;
for (SCSIZE i = 0; i < nCountXN; i++)
if (!pMatNewX->IsValue(i))
{
PushIllegalArgument();
return;
}
}
ScMatrixRef pResMat; // size depends on nCase
if (nCase == 1)
pResMat = GetNewMat(nCXN,nRXN, /*bEmpty*/true);
else
{
if (nCase==2)
pResMat = GetNewMat(1,nRXN, /*bEmpty*/true);
else
pResMat = GetNewMat(nCXN,1, /*bEmpty*/true);
}
if (!pResMat)
{
PushError(FormulaError::CodeOverflow);
return;
}
// Uses sum(x-MeanX)^2 and not [sum x^2]-N * MeanX^2 in case bConstant.
// Clone constant matrices, so that Mat = Mat - Mean is possible.
double fMeanY = 0.0;
if (bConstant)
{
ScMatrixRef pCopyX = pMatX->CloneIfConst();
ScMatrixRef pCopyY = pMatY->CloneIfConst();
if (!pCopyX || !pCopyY)
{
PushError(FormulaError::MatrixSize);
return;
}
pMatX = pCopyX;
pMatY = pCopyY;
// DeltaY is possible here; DeltaX depends on nCase, so later
fMeanY = lcl_GetMeanOverAll(pMatY, N);
for (SCSIZE i=0; i<N; i++)
{
pMatY->PutDouble( ::rtl::math::approxSub(pMatY->GetDouble(i),fMeanY), i );
}
}
if (nCase==1)
{
// calculate simple regression
double fMeanX = 0.0;
if (bConstant)
{ // Mat = Mat - Mean
fMeanX = lcl_GetMeanOverAll(pMatX, N);
for (SCSIZE i=0; i<N; i++)
{
pMatX->PutDouble( ::rtl::math::approxSub(pMatX->GetDouble(i),fMeanX), i );
}
}
double fSumXY = lcl_GetSumProduct(pMatX,pMatY,N);
double fSumX2 = lcl_GetSumProduct(pMatX,pMatX,N);
if (fSumX2==0.0)
{
PushNoValue(); // all x-values are identical
return;
}
double fSlope = fSumXY / fSumX2;
double fHelp;
if (bConstant)
{
double fIntercept = fMeanY - fSlope * fMeanX;
for (SCSIZE i = 0; i < nCountXN; i++)
{
fHelp = pMatNewX->GetDouble(i)*fSlope + fIntercept;
pResMat->PutDouble(_bGrowth ? exp(fHelp) : fHelp, i);
}
}
else
{
for (SCSIZE i = 0; i < nCountXN; i++)
{
fHelp = pMatNewX->GetDouble(i)*fSlope;
pResMat->PutDouble(_bGrowth ? exp(fHelp) : fHelp, i);
}
}
}
else // calculate multiple regression;
{
if (nCase ==2) // Y is column
{
::std::vector< double> aVecR(N); // for QR decomposition
// Enough memory for needed matrices?
ScMatrixRef pMeans = GetNewMat(K, 1, /*bEmpty*/true); // mean of each column
ScMatrixRef pSlopes = GetNewMat(1,K, /*bEmpty*/true); // from b1 to bK
if (!pMeans || !pSlopes)
{
PushError(FormulaError::CodeOverflow);
return;
}
if (bConstant)
{
lcl_CalculateColumnMeans(pMatX, pMeans, K, N);
lcl_CalculateColumnsDelta(pMatX, pMeans, K, N);
}
if (!lcl_CalculateQRdecomposition(pMatX, aVecR, K, N))
{
PushNoValue();
return;
}
// Later on we will divide by elements of aVecR, so make sure
// that they aren't zero.
bool bIsSingular=false;
for (SCSIZE row=0; row < K && !bIsSingular; row++)
bIsSingular = aVecR[row] == 0.0;
if (bIsSingular)
{
PushNoValue();
return;
}
// Z := Q' Y; Y is overwritten with result Z
for (SCSIZE col = 0; col < K; col++)
{
lcl_ApplyHouseholderTransformation(pMatX, col, pMatY, N);
}
// B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
// result Z should have zeros for index>=K; if not, ignore values
for (SCSIZE col = 0; col < K ; col++)
{
pSlopes->PutDouble( pMatY->GetDouble(col), col);
}
lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, false);
// Fill result matrix
lcl_MFastMult(pMatNewX,pSlopes,pResMat,nRXN,K,1);
if (bConstant)
{
double fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K);
for (SCSIZE row = 0; row < nRXN; row++)
pResMat->PutDouble(pResMat->GetDouble(row)+fIntercept, row);
}
if (_bGrowth)
{
for (SCSIZE i = 0; i < nRXN; i++)
pResMat->PutDouble(exp(pResMat->GetDouble(i)), i);
}
}
else
{ // nCase == 3, Y is row, all matrices are transposed
::std::vector< double> aVecR(N); // for QR decomposition
// Enough memory for needed matrices?
ScMatrixRef pMeans = GetNewMat(1, K, /*bEmpty*/true); // mean of each row
ScMatrixRef pSlopes = GetNewMat(K,1, /*bEmpty*/true); // row from b1 to bK
if (!pMeans || !pSlopes)
{
PushError(FormulaError::CodeOverflow);
return;
}
if (bConstant)
{
lcl_CalculateRowMeans(pMatX, pMeans, N, K);
lcl_CalculateRowsDelta(pMatX, pMeans, N, K);
}
if (!lcl_TCalculateQRdecomposition(pMatX, aVecR, K, N))
{
PushNoValue();
return;
}
// Later on we will divide by elements of aVecR, so make sure
// that they aren't zero.
bool bIsSingular=false;
for (SCSIZE row=0; row < K && !bIsSingular; row++)
bIsSingular = aVecR[row] == 0.0;
if (bIsSingular)
{
PushNoValue();
return;
}
// Z := Q' Y; Y is overwritten with result Z
for (SCSIZE row = 0; row < K; row++)
{
lcl_TApplyHouseholderTransformation(pMatX, row, pMatY, N);
}
// B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
// result Z should have zeros for index>=K; if not, ignore values
for (SCSIZE col = 0; col < K ; col++)
{
pSlopes->PutDouble( pMatY->GetDouble(col), col);
}
lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, true);
// Fill result matrix
lcl_MFastMult(pSlopes,pMatNewX,pResMat,1,K,nCXN);
if (bConstant)
{
double fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K);
for (SCSIZE col = 0; col < nCXN; col++)
pResMat->PutDouble(pResMat->GetDouble(col)+fIntercept, col);
}
if (_bGrowth)
{
for (SCSIZE i = 0; i < nCXN; i++)
pResMat->PutDouble(exp(pResMat->GetDouble(i)), i);
}
}
}
PushMatrix(pResMat);
}
void ScInterpreter::ScMatRef()
{
// In case it contains relative references resolve them as usual.
Push( *pCur );
ScAddress aAdr;
PopSingleRef( aAdr );
ScRefCellValue aCell(mrDoc, aAdr);
if (aCell.getType() != CELLTYPE_FORMULA)
{
PushError( FormulaError::NoRef );
return;
}
if (aCell.getFormula()->IsRunning())
{
// Twisted odd corner case where an array element's cell tries to
// access the top left matrix while it is still running, see tdf#88737
// This is a hackish workaround, not a general solution, the matrix
// isn't available anyway and FormulaError::CircularReference would be set.
PushError( FormulaError::RetryCircular );
return;
}
const ScMatrix* pMat = aCell.getFormula()->GetMatrix();
if (pMat)
{
SCSIZE nCols, nRows;
pMat->GetDimensions( nCols, nRows );
SCSIZE nC = static_cast<SCSIZE>(aPos.Col() - aAdr.Col());
SCSIZE nR = static_cast<SCSIZE>(aPos.Row() - aAdr.Row());
#if 0
// XXX: this could be an additional change for tdf#145085 to not
// display the URL in a voluntary entered 2-rows array context.
// However, that might as well be used on purpose to implement a check
// on the URL, which existing documents may have done, the more that
// before the accompanying change of
// ScFormulaCell::GetResultDimensions() the cell array was forced to
// two rows. Do not change without compelling reason. Note that this
// repeating top cell is what Excel implements, but it has no
// additional value so probably isn't used there. Exporting to and
// using in Excel though will lose this capability.
if (aCell.mpFormula->GetCode()->IsHyperLink())
{
// Row 2 element is the URL that is not to be displayed, fake a
// 1-row cell-text-only matrix that is repeated.
assert(nRows == 2);
nR = 0;
}
#endif
if ((nCols <= nC && nCols != 1) || (nRows <= nR && nRows != 1))
PushNA();
else
{
const ScMatrixValue nMatVal = pMat->Get( nC, nR);
ScMatValType nMatValType = nMatVal.nType;
if (ScMatrix::IsNonValueType( nMatValType))
{
if (ScMatrix::IsEmptyPathType( nMatValType))
{ // result of empty false jump path
nFuncFmtType = SvNumFormatType::LOGICAL;
PushInt(0);
}
else if (ScMatrix::IsEmptyType( nMatValType))
{
// Not inherited (really?) and display as empty string, not 0.
PushTempToken( new ScEmptyCellToken( false, true));
}
else
PushString( nMatVal.GetString() );
}
else
{
// Determine nFuncFmtType type before PushDouble().
mrDoc.GetNumberFormatInfo(mrContext, nCurFmtType, nCurFmtIndex, aAdr);
nFuncFmtType = nCurFmtType;
nFuncFmtIndex = nCurFmtIndex;
PushDouble(nMatVal.fVal); // handles DoubleError
}
}
}
else
{
// Determine nFuncFmtType type before PushDouble().
mrDoc.GetNumberFormatInfo(mrContext, nCurFmtType, nCurFmtIndex, aAdr);
nFuncFmtType = nCurFmtType;
nFuncFmtIndex = nCurFmtIndex;
// If not a result matrix, obtain the cell value.
FormulaError nErr = aCell.getFormula()->GetErrCode();
if (nErr != FormulaError::NONE)
PushError( nErr );
else if (aCell.getFormula()->IsValue())
PushDouble(aCell.getFormula()->GetValue());
else
{
svl::SharedString aVal = aCell.getFormula()->GetString();
PushString( aVal );
}
}
}
void ScInterpreter::ScInfo()
{
if( !MustHaveParamCount( GetByte(), 1 ) )
return;
OUString aStr = GetString().getString();
ScCellKeywordTranslator::transKeyword(aStr, &ScGlobal::GetLocale(), ocInfo);
if( aStr == "SYSTEM" )
PushString( u"" SC_INFO_OSVERSION ""_ustr );
else if( aStr == "OSVERSION" )
#if (defined LINUX || defined __FreeBSD__)
PushString(Application::GetOSVersion());
#elif defined MACOSX
// TODO tdf#140286 handle MACOSX version to get result compatible to Excel
PushString("Windows (32-bit) NT 5.01");
#else // handle Windows (WNT, WIN_NT, WIN32, _WIN32)
// TODO tdf#140286 handle Windows version to get a result compatible to Excel
PushString( "Windows (32-bit) NT 5.01" );
#endif
else if( aStr == "RELEASE" )
PushString( ::utl::Bootstrap::getBuildIdData( OUString() ) );
else if( aStr == "NUMFILE" )
PushDouble( 1 );
else if( aStr == "RECALC" )
PushString( ScResId( mrDoc.GetAutoCalc() ? STR_RECALC_AUTO : STR_RECALC_MANUAL ) );
else if (aStr == "DIRECTORY" || aStr == "MEMAVAIL" || aStr == "MEMUSED" || aStr == "ORIGIN" || aStr == "TOTMEM")
PushNA();
else
PushIllegalArgument();
}
/* vim:set shiftwidth=4 softtabstop=4 expandtab: */