/************************************************************** * * 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 * * Unless required by applicable law or agreed to in writing, * software distributed under the License is distributed on an * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY * KIND, either express or implied. See the License for the * specific language governing permissions and limitations * under the License. * *************************************************************/ // MARKER(update_precomp.py): autogen include statement, do not remove #include "precompiled_sc.hxx" // INCLUDE --------------------------------------------------------------- #ifndef INCLUDED_RTL_MATH_HXX #include #endif #include #include #include #include #if OSL_DEBUG_LEVEL > 1 #include #endif #include #include #include "interpre.hxx" #include "global.hxx" #include "compiler.hxx" #include "cell.hxx" #include "document.hxx" #include "dociter.hxx" #include "scmatrix.hxx" #include "globstr.hrc" #include "cellkeytranslator.hxx" #include "osversiondef.hxx" #include #include #include using ::std::vector; using namespace formula; const double fInvEpsilon = 1.0E-7; // ----------------------------------------------------------------------- struct MatrixAdd : public ::std::binary_function { inline double operator() (const double& lhs, const double& rhs) const { return ::rtl::math::approxAdd( lhs,rhs); } }; struct MatrixSub : public ::std::binary_function { inline double operator() (const double& lhs, const double& rhs) const { return ::rtl::math::approxSub( lhs,rhs); } }; struct MatrixMul : public ::std::binary_function { inline double operator() (const double& lhs, const double& rhs) const { return lhs * rhs; } }; struct MatrixDiv : public ::std::binary_function { inline double operator() (const double& lhs, const double& rhs) const { return ScInterpreter::div( lhs,rhs); } }; struct MatrixPow : public ::std::binary_function { inline double operator() (const double& lhs, const double& rhs) const { return ::pow( lhs,rhs); } }; double ScInterpreter::ScGetGCD(double fx, double fy) { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::div" ); // 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() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScGCD" ); short nParamCount = GetByte(); if ( MustHaveParamCountMin( nParamCount, 1 ) ) { double fx, fy = 0.0; ScRange aRange; size_t nRefInList = 0; while (!nGlobalError && 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 : { sal_uInt16 nErr = 0; PopDoubleRef( aRange, nParamCount, nRefInList); double nCellVal; ScValueIterator aValIter(pDok, aRange, glSubTotal); if (aValIter.GetFirst(nCellVal, nErr)) { do { fx = ::rtl::math::approxFloor( nCellVal); if (fx < 0.0) { PushIllegalArgument(); return; } fy = ScGetGCD(fx, fy); } while (nErr == 0 && aValIter.GetNext(nCellVal, nErr)); } SetError(nErr); } break; case svMatrix : { ScMatrixRef pMat = PopMatrix(); if (pMat) { SCSIZE nC, nR; pMat->GetDimensions(nC, nR); if (nC == 0 || nR == 0) SetError(errIllegalArgument); else { SCSIZE nCount = nC * nR; for ( SCSIZE j = 0; j < nCount; j++ ) { if (!pMat->IsValue(j)) { PushIllegalArgument(); return; } fx = ::rtl::math::approxFloor( pMat->GetDouble(j)); if (fx < 0.0) { PushIllegalArgument(); return; } fy = ScGetGCD(fx, fy); } } } } break; default : SetError(errIllegalParameter); break; } } PushDouble(fy); } } void ScInterpreter:: ScLCM() { short nParamCount = GetByte(); if ( MustHaveParamCountMin( nParamCount, 1 ) ) { double fx, fy = 1.0; ScRange aRange; size_t nRefInList = 0; while (!nGlobalError && 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 : { sal_uInt16 nErr = 0; PopDoubleRef( aRange, nParamCount, nRefInList); double nCellVal; ScValueIterator aValIter(pDok, aRange, glSubTotal); 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 == 0 && aValIter.GetNext(nCellVal, nErr)); } SetError(nErr); } break; case svMatrix : { ScMatrixRef pMat = PopMatrix(); if (pMat) { SCSIZE nC, nR; pMat->GetDimensions(nC, nR); if (nC == 0 || nR == 0) SetError(errIllegalArgument); else { SCSIZE nCount = nC * nR; for ( SCSIZE j = 0; j < nCount; j++ ) { if (!pMat->IsValue(j)) { PushIllegalArgument(); return; } fx = ::rtl::math::approxFloor( pMat->GetDouble(j)); if (fx < 0.0) { PushIllegalArgument(); return; } if (fx == 0.0 || fy == 0.0) fy = 0.0; else fy = fx * fy / ScGetGCD(fx, fy); } } } } break; default : SetError(errIllegalParameter); break; } } PushDouble(fy); } } ScMatrixRef ScInterpreter::GetNewMat(SCSIZE nC, SCSIZE nR) { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetNewMat" ); ScMatrix* pMat = new ScMatrix( nC, nR); pMat->SetErrorInterpreter( this); // A temporary matrix is mutable and ScMatrix::CloneIfConst() returns the // very matrix. pMat->SetImmutable( false); SCSIZE nCols, nRows; pMat->GetDimensions( nCols, nRows); if ( nCols != nC || nRows != nR ) { // arbitray limit of elements exceeded SetError( errStackOverflow); pMat->Delete(); pMat = NULL; } return pMat; } ScMatrixRef ScInterpreter::CreateMatrixFromDoubleRef( const FormulaToken* pToken, SCCOL nCol1, SCROW nRow1, SCTAB nTab1, SCCOL nCol2, SCROW nRow2, SCTAB nTab2 ) { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CreateMatrixFromDoubleRef" ); ScMatrixRef pMat = NULL; if (nTab1 == nTab2 && !nGlobalError) { ScTokenMatrixMap::const_iterator aIter; if ( static_cast(nRow2 - nRow1 + 1) * static_cast(nCol2 - nCol1 + 1) > ScMatrix::GetElementsMax() ) SetError(errStackOverflow); else if (pTokenMatrixMap && ((aIter = pTokenMatrixMap->find( pToken)) != pTokenMatrixMap->end())) pMat = static_cast((*aIter).second.get())->GetMatrix(); else { SCSIZE nMatCols = static_cast(nCol2 - nCol1 + 1); SCSIZE nMatRows = static_cast(nRow2 - nRow1 + 1); pMat = GetNewMat( nMatCols, nMatRows); if (pMat && !nGlobalError) { // 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( pDok, nCol1, nRow1, nTab1, nCol2, nRow2, nTab2); for (ScBaseCell* pCell = aCellIter.GetFirst(); pCell; pCell = aCellIter.GetNext()) { nThisCol = aCellIter.GetCol(); nThisRow = aCellIter.GetRow(); if (nThisCol != nNextCol || nThisRow != nNextRow) { // Fill empty between iterator's positions. for ( ; nNextCol <= nThisCol; ++nNextCol) { SCSIZE nC = nNextCol - nCol1; 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; } if (HasCellEmptyData(pCell)) { pMat->PutEmpty( static_cast(nThisCol-nCol1), static_cast(nThisRow-nRow1)); } else if (HasCellValueData(pCell)) { ScAddress aAdr( nThisCol, nThisRow, nTab1); double fVal = GetCellValue( aAdr, pCell); if ( nGlobalError ) { fVal = CreateDoubleError( nGlobalError); nGlobalError = 0; } pMat->PutDouble( fVal, static_cast(nThisCol-nCol1), static_cast(nThisRow-nRow1)); } else { String aStr; GetCellString( aStr, pCell); if ( nGlobalError ) { double fVal = CreateDoubleError( nGlobalError); nGlobalError = 0; pMat->PutDouble( fVal, static_cast(nThisCol-nCol1), static_cast(nThisRow-nRow1)); } else pMat->PutString( aStr, static_cast(nThisCol-nCol1), static_cast(nThisRow-nRow1)); } } // 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 (pTokenMatrixMap) pTokenMatrixMap->insert( ScTokenMatrixMap::value_type( pToken, new ScMatrixToken( pMat))); } } } else // not a 2D matrix SetError(errIllegalParameter); return pMat; } ScMatrixRef ScInterpreter::GetMatrix() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetMatrix" ); ScMatrixRef pMat = NULL; switch (GetRawStackType()) { case svSingleRef : { ScAddress aAdr; PopSingleRef( aAdr ); pMat = GetNewMat(1, 1); if (pMat) { ScBaseCell* pCell = GetCell( aAdr ); if (HasCellEmptyData(pCell)) pMat->PutEmpty( 0 ); else if (HasCellValueData(pCell)) pMat->PutDouble(GetCellValue(aAdr, pCell), 0); else { String aStr; GetCellString(aStr, pCell); pMat->PutString(aStr, 0); } } } break; case svDoubleRef: { SCCOL nCol1, nCol2; SCROW nRow1, nRow2; SCTAB nTab1, nTab2; const ScToken* p = sp ? static_cast(pStack[sp-1]) : NULL; 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 ) { fVal = CreateDoubleError( nGlobalError); nGlobalError = 0; } pMat->PutDouble( fVal, 0); } } break; case svString : { String aStr = GetString(); pMat = GetNewMat( 1, 1); if ( pMat ) { if ( nGlobalError ) { double fVal = CreateDoubleError( nGlobalError); pMat->PutDouble( fVal, 0); nGlobalError = 0; } else pMat->PutString( aStr, 0); } } break; default: PopError(); SetError( errIllegalArgument); break; } return pMat; } void ScInterpreter::ScMatValue() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatValue" ); if ( MustHaveParamCount( GetByte(), 3 ) ) { // 0 to count-1 SCSIZE nR = static_cast(::rtl::math::approxFloor(GetDouble())); SCSIZE nC = static_cast(::rtl::math::approxFloor(GetDouble())); switch (GetStackType()) { case svSingleRef : { ScAddress aAdr; PopSingleRef( aAdr ); ScBaseCell* pCell = GetCell( aAdr ); if (pCell && pCell->GetCellType() == CELLTYPE_FORMULA) { sal_uInt16 nErrCode = ((ScFormulaCell*)pCell)->GetErrCode(); if (nErrCode != 0) PushError( nErrCode); else { const ScMatrix* pMat = ((ScFormulaCell*)pCell)->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(nR) && nRow2 - nRow1 >= static_cast(nC) && nTab1 == nTab2) { ScAddress aAdr( sal::static_int_cast( nCol1 + nR ), sal::static_int_cast( nRow1 + nC ), nTab1 ); ScBaseCell* pCell = GetCell( aAdr ); if (HasCellValueData(pCell)) PushDouble(GetCellValue( aAdr, pCell )); else { String aStr; GetCellString(aStr, pCell); PushString(aStr); } } else PushNoValue(); } break; case svMatrix: { ScMatrixRef pMat = PopMatrix(); CalculateMatrixValue(pMat,nC,nR); } break; default: PopError(); PushIllegalParameter(); break; } } } void ScInterpreter::CalculateMatrixValue(const ScMatrix* pMat,SCSIZE nC,SCSIZE nR) { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CalculateMatrixValue" ); if (pMat) { SCSIZE nCl, nRw; pMat->GetDimensions(nCl, nRw); if (nC < nCl && nR < nRw) { ScMatValType nMatValType; const ScMatrixValue* pMatVal = pMat->Get( nC, nR,nMatValType); if (ScMatrix::IsNonValueType( nMatValType)) PushString( pMatVal->GetString() ); else PushDouble(pMatVal->fVal); // also handles DoubleError } else PushNoValue(); } else PushNoValue(); } void ScInterpreter::ScEMat() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScEMat" ); if ( MustHaveParamCount( GetByte(), 1 ) ) { SCSIZE nDim = static_cast(::rtl::math::approxFloor(GetDouble())); if ( nDim * nDim > ScMatrix::GetElementsMax() || nDim == 0) PushIllegalArgument(); else { ScMatrixRef pRMat = GetNewMat(nDim, nDim); if (pRMat) { MEMat(pRMat, nDim); PushMatrix(pRMat); } else PushIllegalArgument(); } } } void ScInterpreter::MEMat(ScMatrix* mM, SCSIZE n) { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::MEMat" ); 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 occured, 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 diagonale. 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 fTmp = mA->GetDouble( k, i) / mA->GetDouble( k, k); mA->PutDouble( fTmp, k, i); for (SCSIZE j = k+1; j < n; ++j) mA->PutDouble( mA->GetDouble( j, i) - fTmp * mA->GetDouble( j, k), j, i); } } #if OSL_DEBUG_LEVEL > 1 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) { double 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) nFirst = i; X[i] = fSum; // X[i] === y[i] } // Solve for x in Ux=y using back substitution. for (SCSIZE i = n; i--; ) { double 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 / mLU->GetDouble( i, i); // X[i] === x[i] } #if OSL_DEBUG_LEVEL >1 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() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatDet" ); if ( MustHaveParamCount( GetByte(), 1 ) ) { ScMatrixRef pMat = GetMatrix(); if (!pMat) { PushIllegalParameter(); return; } if ( !pMat->IsNumeric() ) { PushNoValue(); return; } SCSIZE nC, nR; pMat->GetDimensions(nC, nR); if ( nC != nR || nC == 0 || (sal_uLong) nC * nC > ScMatrix::GetElementsMax() ) PushIllegalArgument(); else { // LUP decomposition is done inplace, use copy. ScMatrixRef xLU = pMat->Clone(); if (!xLU) PushError( errCodeOverflow); else { ::std::vector< SCSIZE> P(nR); int nDetSign = lcl_LUP_decompose( xLU, 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; ScMatrix* pLU = xLU; for (SCSIZE i=0; i < nR; ++i) fDet *= pLU->GetDouble( i, i); PushDouble( fDet); } } } } } void ScInterpreter::ScMatInv() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatInv" ); if ( MustHaveParamCount( GetByte(), 1 ) ) { ScMatrixRef pMat = GetMatrix(); if (!pMat) { PushIllegalParameter(); return; } if ( !pMat->IsNumeric() ) { PushNoValue(); return; } SCSIZE nC, nR; pMat->GetDimensions(nC, nR); if ( nC != nR || nC == 0 || (sal_uLong) nC * nC > ScMatrix::GetElementsMax() ) PushIllegalArgument(); else { // LUP decomposition is done inplace, use copy. ScMatrixRef xLU = pMat->Clone(); // The result matrix. ScMatrixRef xY = GetNewMat( nR, nR); if (!xLU || !xY) PushError( errCodeOverflow); else { ::std::vector< SCSIZE> P(nR); int nDetSign = lcl_LUP_decompose( xLU, nR, P); if (!nDetSign) PushIllegalArgument(); else { // Solve equation for each column. ScMatrix* pY = xY; ::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, nR, P, B, X); for (SCSIZE i=0; i < nR; ++i) pY->PutDouble( X[i], j, i); } #if 0 /* 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. */ ScMatrixRef xR = GetNewMat( nR, nR); if (xR) { ScMatrix* pR = xR; lcl_MFastMult( pMat, pY, pR, nR, nR, nR); #if OSL_DEBUG_LEVEL > 1 fprintf( stderr, "\n%s\n", "ScMatInv(): mult-identity"); #endif for (SCSIZE i=0; i < nR; ++i) { for (SCSIZE j=0; j < nR; ++j) { double fTmp = pR->GetDouble( j, i); #if OSL_DEBUG_LEVEL > 1 fprintf( stderr, "%8.2g ", fTmp); #endif if (fabs( fTmp - (i == j)) > fInvEpsilon) SetError( errIllegalArgument); } #if OSL_DEBUG_LEVEL > 1 fprintf( stderr, "\n%s\n", ""); #endif } } #endif if (nGlobalError) PushError( nGlobalError); else PushMatrix( pY); } } } } } void ScInterpreter::ScMatMult() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatMult" ); if ( MustHaveParamCount( GetByte(), 2 ) ) { 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); if (pRMat) { double sum; for (SCSIZE i = 0; i < nR1; i++) { for (SCSIZE j = 0; j < nC2; j++) { sum = 0.0; for (SCSIZE k = 0; k < nC1; k++) { sum += pMat1->GetDouble(k,i)*pMat2->GetDouble(j,k); } pRMat->PutDouble(sum, j, i); } } PushMatrix(pRMat); } else PushIllegalArgument(); } } else PushNoValue(); } else PushIllegalParameter(); } } void ScInterpreter::ScMatTrans() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatTrans" ); if ( MustHaveParamCount( GetByte(), 1 ) ) { ScMatrixRef pMat = GetMatrix(); ScMatrixRef pRMat; if (pMat) { SCSIZE nC, nR; pMat->GetDimensions(nC, nR); pRMat = GetNewMat(nR, nC); 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. */ inline 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; } template ScMatrixRef lcl_MatrixCalculation(const _Function& _pOperation,ScMatrix* pMat1, ScMatrix* pMat2,ScInterpreter* _pIterpreter) { SCSIZE nC1, nC2, nMinC; SCSIZE nR1, nR2, nMinR; SCSIZE i, j; pMat1->GetDimensions(nC1, nR1); pMat2->GetDimensions(nC2, nR2); nMinC = lcl_GetMinExtent( nC1, nC2); nMinR = lcl_GetMinExtent( nR1, nR2); ScMatrixRef xResMat = _pIterpreter->GetNewMat(nMinC, nMinR); if (xResMat) { ScMatrix* pResMat = xResMat; for (i = 0; i < nMinC; i++) { for (j = 0; j < nMinR; j++) { if (pMat1->IsValueOrEmpty(i,j) && pMat2->IsValueOrEmpty(i,j)) { double d = _pOperation(pMat1->GetDouble(i,j),pMat2->GetDouble(i,j)); pResMat->PutDouble( d, i, j); } else pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i, j); } } } return xResMat; } ScMatrixRef ScInterpreter::MatConcat(ScMatrix* pMat1, ScMatrix* pMat2) { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::MatConcat" ); SCSIZE nC1, nC2, nMinC; SCSIZE nR1, nR2, nMinR; SCSIZE i, j; pMat1->GetDimensions(nC1, nR1); pMat2->GetDimensions(nC2, nR2); nMinC = lcl_GetMinExtent( nC1, nC2); nMinR = lcl_GetMinExtent( nR1, nR2); ScMatrixRef xResMat = GetNewMat(nMinC, nMinR); if (xResMat) { ScMatrix* pResMat = xResMat; for (i = 0; i < nMinC; i++) { for (j = 0; j < nMinR; j++) { sal_uInt16 nErr = pMat1->GetErrorIfNotString( i, j); if (!nErr) nErr = pMat2->GetErrorIfNotString( i, j); if (nErr) pResMat->PutError( nErr, i, j); else { String aTmp( pMat1->GetString( *pFormatter, i, j)); aTmp += pMat2->GetString( *pFormatter, i, j); pResMat->PutString( aTmp, i, j); } } } } return xResMat; } // fuer DATE, TIME, DATETIME void lcl_GetDiffDateTimeFmtType( short& nFuncFmt, short nFmt1, short nFmt2 ) { if ( nFmt1 != NUMBERFORMAT_UNDEFINED || nFmt2 != NUMBERFORMAT_UNDEFINED ) { if ( nFmt1 == nFmt2 ) { if ( nFmt1 == NUMBERFORMAT_TIME || nFmt1 == NUMBERFORMAT_DATETIME ) nFuncFmt = NUMBERFORMAT_TIME; // Zeiten ergeben Zeit // else: nichts besonderes, Zahl (Datum - Datum := Tage) } else if ( nFmt1 == NUMBERFORMAT_UNDEFINED ) nFuncFmt = nFmt2; // z.B. Datum + Tage := Datum else if ( nFmt2 == NUMBERFORMAT_UNDEFINED ) nFuncFmt = nFmt1; else { if ( nFmt1 == NUMBERFORMAT_DATE || nFmt2 == NUMBERFORMAT_DATE || nFmt1 == NUMBERFORMAT_DATETIME || nFmt2 == NUMBERFORMAT_DATETIME ) { if ( nFmt1 == NUMBERFORMAT_TIME || nFmt2 == NUMBERFORMAT_TIME ) nFuncFmt = NUMBERFORMAT_DATETIME; // Datum + Zeit } } } } void ScInterpreter::ScAdd() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScAdd" ); CalculateAddSub(sal_False); } void ScInterpreter::CalculateAddSub(sal_Bool _bSub) { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CalculateAddSub" ); ScMatrixRef pMat1 = NULL; ScMatrixRef pMat2 = NULL; double fVal1 = 0.0, fVal2 = 0.0; short nFmt1, nFmt2; nFmt1 = nFmt2 = NUMBERFORMAT_UNDEFINED; short nFmtCurrencyType = nCurFmtType; sal_uLong nFmtCurrencyIndex = nCurFmtIndex; short nFmtPercentType = nCurFmtType; if ( GetStackType() == svMatrix ) pMat2 = GetMatrix(); else { fVal2 = GetDouble(); switch ( nCurFmtType ) { case NUMBERFORMAT_DATE : case NUMBERFORMAT_TIME : case NUMBERFORMAT_DATETIME : nFmt2 = nCurFmtType; break; case NUMBERFORMAT_CURRENCY : nFmtCurrencyType = nCurFmtType; nFmtCurrencyIndex = nCurFmtIndex; break; case NUMBERFORMAT_PERCENT : nFmtPercentType = NUMBERFORMAT_PERCENT; break; } } if ( GetStackType() == svMatrix ) pMat1 = GetMatrix(); else { fVal1 = GetDouble(); switch ( nCurFmtType ) { case NUMBERFORMAT_DATE : case NUMBERFORMAT_TIME : case NUMBERFORMAT_DATETIME : nFmt1 = nCurFmtType; break; case NUMBERFORMAT_CURRENCY : nFmtCurrencyType = nCurFmtType; nFmtCurrencyIndex = nCurFmtIndex; break; case NUMBERFORMAT_PERCENT : nFmtPercentType = NUMBERFORMAT_PERCENT; break; } } if (pMat1 && pMat2) { ScMatrixRef pResMat; if ( _bSub ) { MatrixSub aSub; pResMat = lcl_MatrixCalculation(aSub ,pMat1, pMat2,this); } else { MatrixAdd aAdd; pResMat = lcl_MatrixCalculation(aAdd ,pMat1, pMat2,this); } if (!pResMat) PushNoValue(); else PushMatrix(pResMat); } else if (pMat1 || pMat2) { double fVal; sal_Bool bFlag; ScMatrixRef pMat = pMat1; if (!pMat) { fVal = fVal1; pMat = pMat2; bFlag = sal_True; // double - Matrix } else { fVal = fVal2; bFlag = sal_False; // Matrix - double } SCSIZE nC, nR; pMat->GetDimensions(nC, nR); ScMatrixRef pResMat = GetNewMat(nC, nR); if (pResMat) { SCSIZE nCount = nC * nR; if (bFlag || !_bSub ) { for ( SCSIZE i = 0; i < nCount; i++ ) { if (pMat->IsValue(i)) pResMat->PutDouble( _bSub ? ::rtl::math::approxSub( fVal, pMat->GetDouble(i)) : ::rtl::math::approxAdd( pMat->GetDouble(i), fVal), i); else pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i); } // for ( SCSIZE i = 0; i < nCount; i++ ) } // if (bFlag || !_bSub ) else { for ( SCSIZE i = 0; i < nCount; i++ ) { if (pMat->IsValue(i)) pResMat->PutDouble( ::rtl::math::approxSub( pMat->GetDouble(i), fVal), i); else pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i); } // for ( SCSIZE i = 0; i < nCount; i++ ) } PushMatrix(pResMat); } else PushIllegalArgument(); } else if ( _bSub ) PushDouble( ::rtl::math::approxSub( fVal1, fVal2 ) ); else PushDouble( ::rtl::math::approxAdd( fVal1, fVal2 ) ); if ( nFmtCurrencyType == NUMBERFORMAT_CURRENCY ) { nFuncFmtType = nFmtCurrencyType; nFuncFmtIndex = nFmtCurrencyIndex; } else { lcl_GetDiffDateTimeFmtType( nFuncFmtType, nFmt1, nFmt2 ); if ( nFmtPercentType == NUMBERFORMAT_PERCENT && nFuncFmtType == NUMBERFORMAT_NUMBER ) nFuncFmtType = NUMBERFORMAT_PERCENT; } } void ScInterpreter::ScAmpersand() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScAmpersand" ); ScMatrixRef pMat1 = NULL; ScMatrixRef pMat2 = NULL; String sStr1, sStr2; if ( GetStackType() == svMatrix ) pMat2 = GetMatrix(); else sStr2 = GetString(); if ( GetStackType() == svMatrix ) pMat1 = GetMatrix(); else sStr1 = GetString(); if (pMat1 && pMat2) { ScMatrixRef pResMat = MatConcat(pMat1, pMat2); if (!pResMat) PushNoValue(); else PushMatrix(pResMat); } else if (pMat1 || pMat2) { String sStr; sal_Bool bFlag; ScMatrixRef pMat = pMat1; if (!pMat) { sStr = sStr1; pMat = pMat2; bFlag = sal_True; // double - Matrix } else { sStr = sStr2; bFlag = sal_False; // Matrix - double } SCSIZE nC, nR; pMat->GetDimensions(nC, nR); ScMatrixRef pResMat = GetNewMat(nC, nR); if (pResMat) { SCSIZE nCount = nC * nR; if (nGlobalError) { for ( SCSIZE i = 0; i < nCount; i++ ) pResMat->PutError( nGlobalError, i); } else if (bFlag) { for ( SCSIZE i = 0; i < nCount; i++ ) { sal_uInt16 nErr = pMat->GetErrorIfNotString( i); if (nErr) pResMat->PutError( nErr, i); else { String aTmp( sStr); aTmp += pMat->GetString( *pFormatter, i); pResMat->PutString( aTmp, i); } } } else { for ( SCSIZE i = 0; i < nCount; i++ ) { sal_uInt16 nErr = pMat->GetErrorIfNotString( i); if (nErr) pResMat->PutError( nErr, i); else { String aTmp( pMat->GetString( *pFormatter, i)); aTmp += sStr; pResMat->PutString( aTmp, i); } } } PushMatrix(pResMat); } else PushIllegalArgument(); } else { if ( CheckStringResultLen( sStr1, sStr2 ) ) sStr1 += sStr2; PushString(sStr1); } } void ScInterpreter::ScSub() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSub" ); CalculateAddSub(sal_True); } void ScInterpreter::ScMul() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMul" ); ScMatrixRef pMat1 = NULL; ScMatrixRef pMat2 = NULL; double fVal1 = 0.0, fVal2 = 0.0; short nFmtCurrencyType = nCurFmtType; sal_uLong nFmtCurrencyIndex = nCurFmtIndex; if ( GetStackType() == svMatrix ) pMat2 = GetMatrix(); else { fVal2 = GetDouble(); switch ( nCurFmtType ) { case NUMBERFORMAT_CURRENCY : nFmtCurrencyType = nCurFmtType; nFmtCurrencyIndex = nCurFmtIndex; break; } } if ( GetStackType() == svMatrix ) pMat1 = GetMatrix(); else { fVal1 = GetDouble(); switch ( nCurFmtType ) { case NUMBERFORMAT_CURRENCY : nFmtCurrencyType = nCurFmtType; nFmtCurrencyIndex = nCurFmtIndex; break; } } if (pMat1 && pMat2) { MatrixMul aMul; ScMatrixRef pResMat = lcl_MatrixCalculation(aMul,pMat1, pMat2,this); 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); if (pResMat) { SCSIZE nCount = nC * nR; for ( SCSIZE i = 0; i < nCount; i++ ) if (pMat->IsValue(i)) pResMat->PutDouble(pMat->GetDouble(i)*fVal, i); else pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i); PushMatrix(pResMat); } else PushIllegalArgument(); } else PushDouble(fVal1 * fVal2); if ( nFmtCurrencyType == NUMBERFORMAT_CURRENCY ) { nFuncFmtType = nFmtCurrencyType; nFuncFmtIndex = nFmtCurrencyIndex; } } void ScInterpreter::ScDiv() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScDiv" ); ScMatrixRef pMat1 = NULL; ScMatrixRef pMat2 = NULL; double fVal1 = 0.0, fVal2 = 0.0; short nFmtCurrencyType = nCurFmtType; sal_uLong nFmtCurrencyIndex = nCurFmtIndex; short nFmtCurrencyType2 = NUMBERFORMAT_UNDEFINED; if ( GetStackType() == svMatrix ) pMat2 = GetMatrix(); else { fVal2 = GetDouble(); // hier kein Currency uebernehmen, 123kg/456DM sind nicht DM nFmtCurrencyType2 = nCurFmtType; } if ( GetStackType() == svMatrix ) pMat1 = GetMatrix(); else { fVal1 = GetDouble(); switch ( nCurFmtType ) { case NUMBERFORMAT_CURRENCY : nFmtCurrencyType = nCurFmtType; nFmtCurrencyIndex = nCurFmtIndex; break; } } if (pMat1 && pMat2) { MatrixDiv aDiv; ScMatrixRef pResMat = lcl_MatrixCalculation(aDiv,pMat1, pMat2,this); if (!pResMat) PushNoValue(); else PushMatrix(pResMat); } else if (pMat1 || pMat2) { double fVal; sal_Bool bFlag; ScMatrixRef pMat = pMat1; if (!pMat) { fVal = fVal1; pMat = pMat2; bFlag = sal_True; // double - Matrix } else { fVal = fVal2; bFlag = sal_False; // Matrix - double } SCSIZE nC, nR; pMat->GetDimensions(nC, nR); ScMatrixRef pResMat = GetNewMat(nC, nR); if (pResMat) { SCSIZE nCount = nC * nR; if (bFlag) { for ( SCSIZE i = 0; i < nCount; i++ ) if (pMat->IsValue(i)) pResMat->PutDouble( div( fVal, pMat->GetDouble(i)), i); else pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i); } else { for ( SCSIZE i = 0; i < nCount; i++ ) if (pMat->IsValue(i)) pResMat->PutDouble( div( pMat->GetDouble(i), fVal), i); else pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i); } PushMatrix(pResMat); } else PushIllegalArgument(); } else { PushDouble( div( fVal1, fVal2) ); } if ( nFmtCurrencyType == NUMBERFORMAT_CURRENCY && nFmtCurrencyType2 != NUMBERFORMAT_CURRENCY ) { // auch DM/DM ist nicht DM bzw. DEM/EUR nicht DEM nFuncFmtType = nFmtCurrencyType; nFuncFmtIndex = nFmtCurrencyIndex; } } void ScInterpreter::ScPower() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScPower" ); if ( MustHaveParamCount( GetByte(), 2 ) ) ScPow(); } void ScInterpreter::ScPow() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScPow" ); ScMatrixRef pMat1 = NULL; ScMatrixRef pMat2 = NULL; 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) { MatrixPow aPow; ScMatrixRef pResMat = lcl_MatrixCalculation(aPow,pMat1, pMat2,this); if (!pResMat) PushNoValue(); else PushMatrix(pResMat); } else if (pMat1 || pMat2) { double fVal; sal_Bool bFlag; ScMatrixRef pMat = pMat1; if (!pMat) { fVal = fVal1; pMat = pMat2; bFlag = sal_True; // double - Matrix } else { fVal = fVal2; bFlag = sal_False; // Matrix - double } SCSIZE nC, nR; pMat->GetDimensions(nC, nR); ScMatrixRef pResMat = GetNewMat(nC, nR); if (pResMat) { SCSIZE nCount = nC * nR; if (bFlag) { for ( SCSIZE i = 0; i < nCount; i++ ) if (pMat->IsValue(i)) pResMat->PutDouble(pow(fVal,pMat->GetDouble(i)), i); else pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i); } else { for ( SCSIZE i = 0; i < nCount; i++ ) if (pMat->IsValue(i)) pResMat->PutDouble(pow(pMat->GetDouble(i),fVal), i); else pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i); } PushMatrix(pResMat); } else PushIllegalArgument(); } else PushDouble(pow(fVal1,fVal2)); } void ScInterpreter::ScSumProduct() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSumProduct" ); sal_uInt8 nParamCount = GetByte(); if ( !MustHaveParamCount( nParamCount, 1, 30 ) ) return; ScMatrixRef pMat1 = NULL; ScMatrixRef pMat2 = NULL; ScMatrixRef pMat = NULL; pMat2 = GetMatrix(); if (!pMat2) { PushIllegalParameter(); return; } SCSIZE nC, nC1; SCSIZE nR, nR1; pMat2->GetDimensions(nC, nR); pMat = pMat2; MatrixMul aMul; for (sal_uInt16 i = 1; i < nParamCount; i++) { pMat1 = GetMatrix(); if (!pMat1) { PushIllegalParameter(); return; } pMat1->GetDimensions(nC1, nR1); if (nC1 != nC || nR1 != nR) { PushNoValue(); return; } ScMatrixRef pResMat = lcl_MatrixCalculation(aMul,pMat1, pMat,this); if (!pResMat) { PushNoValue(); return; } else pMat = pResMat; } double fSum = 0.0; SCSIZE nCount = pMat->GetElementCount(); for (SCSIZE j = 0; j < nCount; j++) { if (!pMat->IsString(j)) fSum += pMat->GetDouble(j); } PushDouble(fSum); } void ScInterpreter::ScSumX2MY2() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSumX2MY2" ); CalculateSumX2MY2SumX2DY2(sal_False); } void ScInterpreter::CalculateSumX2MY2SumX2DY2(sal_Bool _bSumX2DY2) { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CalculateSumX2MY2SumX2DY2" ); if ( !MustHaveParamCount( GetByte(), 2 ) ) return; ScMatrixRef pMat1 = NULL; ScMatrixRef pMat2 = NULL; 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, fSum = 0.0; for (i = 0; i < nC1; i++) for (j = 0; j < nR1; j++) if (!pMat1->IsString(i,j) && !pMat2->IsString(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); } void ScInterpreter::ScSumX2DY2() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSumX2DY2" ); CalculateSumX2MY2SumX2DY2(sal_True); } void ScInterpreter::ScSumXMY2() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSumXMY2" ); if ( !MustHaveParamCount( GetByte(), 2 ) ) return; ScMatrixRef pMat1 = NULL; ScMatrixRef pMat2 = NULL; 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; } // if (nC1 != nC2 || nR1 != nR2) MatrixSub aSub; ScMatrixRef pResMat = lcl_MatrixCalculation(aSub,pMat1, pMat2,this); if (!pResMat) { PushNoValue(); } else { double fVal, fSum = 0.0; SCSIZE nCount = pResMat->GetElementCount(); for (SCSIZE i = 0; i < nCount; i++) if (!pResMat->IsString(i)) { fVal = pResMat->GetDouble(i); fSum += fVal * fVal; } PushDouble(fSum); } } void ScInterpreter::ScFrequency() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScFrequency" ); if ( !MustHaveParamCount( GetByte(), 2 ) ) return; vector aBinArray; vector aBinIndexOrder; GetSortArray(1, aBinArray, &aBinIndexOrder); SCSIZE nBinSize = aBinArray.size(); if (nGlobalError) { PushNoValue(); return; } vector aDataArray; GetSortArray(1, aDataArray); SCSIZE nDataSize = aDataArray.size(); if (aDataArray.empty() || nGlobalError) { PushNoValue(); return; } ScMatrixRef pResMat = GetNewMat(1, nBinSize+1); 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(nCount), aBinIndexOrder[j]); } pResMat->PutDouble(static_cast(nDataSize-i), j); PushMatrix(pResMat); } // ----------------------------------------------------------------------------- // Helper methods for LINEST/LOGEST and TREND/GROWTH // All matrices must already exist and have the needed size, no control tests // done. Those methodes, 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). // ----------------------------------------------------------------------------- // Multiply n x m Mat A with m x l Mat B to n x l Mat R void lcl_MFastMult( ScMatrixRef pA, ScMatrixRef pB, ScMatrixRef pR, SCSIZE n, SCSIZE m, SCSIZE l ) { double sum; 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)] sum = 0.0; for (SCSIZE k = 0; k < m; k++) sum += pA->GetDouble(k,row) * pB->GetDouble(col,k); pR->PutDouble(sum, col, row); } } } // over all elements; uses the matrices as vectors of length M double lcl_GetSumProduct( ScMatrixRef pMatA, ScMatrixRef pMatB, SCSIZE nM ) { double fSum = 0.0; for (SCSIZE i=0; iGetDouble(i) * pMatB->GetDouble(i); return fSum; } // 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( ScMatrixRef pMatA, SCSIZE nC, SCSIZE nR, SCSIZE nN ) { double fNorm = 0.0; for (SCSIZE row=nR; rowGetDouble(nC,row)) * (pMatA->GetDouble(nC,row)); return sqrt(fNorm); } // Euclidean norm of row index R starting in column index C; // matrix A has count N columns. double lcl_TGetColumnEuclideanNorm( ScMatrixRef pMatA, SCSIZE nR, SCSIZE nC, SCSIZE nN ) { double fNorm = 0.0; for (SCSIZE col=nC; colGetDouble(col,nR)) * (pMatA->GetDouble(col,nR)); return sqrt(fNorm); } // 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( ScMatrixRef pMatA, SCSIZE nC, SCSIZE nR, SCSIZE nN ) { double fNorm = 0.0; for (SCSIZE row=nR; rowGetDouble(nC,row))) fNorm = fabs(pMatA->GetDouble(nC,row)); return fNorm; } // Maximum norm of row index R starting in col index C; // matrix A has count N columns. double lcl_TGetColumnMaximumNorm( ScMatrixRef pMatA, SCSIZE nR, SCSIZE nC, SCSIZE nN ) { double fNorm = 0.0; for (SCSIZE col=nC; colGetDouble(col,nR))) fNorm = fabs(pMatA->GetDouble(col,nR)); return fNorm; } // Special version for use within QR decomposition. // starting in row index R; // Ca and Cb are indices of columns, matrices A and B have count N rows. double lcl_GetColumnSumProduct( ScMatrixRef pMatA, SCSIZE nCa, ScMatrixRef pMatB, SCSIZE nCb, SCSIZE nR, SCSIZE nN ) { double fResult = 0.0; for (SCSIZE row=nR; rowGetDouble(nCa,row) * pMatB->GetDouble(nCb,row); return fResult; } // starting in column index C; // Ra and Rb are indices of rows, matrices A and B have count N columns. double lcl_TGetColumnSumProduct( ScMatrixRef pMatA, SCSIZE nRa, ScMatrixRef pMatB, SCSIZE nRb, SCSIZE nC, SCSIZE nN ) { double fResult = 0.0; for (SCSIZE col=nC; colGetDouble(col,nRa) * pMatB->GetDouble(col,nRb); return fResult; } // #118029# 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(ScMatrixRef pMatA, ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN) { double fScale ; double fEuclid ; double fFactor ; double fSignum ; double fSum ; // ScMatrix matrices are zero based, index access (column,row) for (SCSIZE col = 0; col PutDouble( pMatA->GetDouble(col,row)/fScale, col, row); fEuclid = lcl_GetColumnEuclideanNorm(pMatA, col, col, nN); fFactor = 1.0/fEuclid/(fEuclid + fabs(pMatA->GetDouble(col,col))); 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; cPutDouble( 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(ScMatrixRef pMatA, ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN) { double fScale ; double fEuclid ; double fFactor ; double fSignum ; double fSum ; // ScMatrix matrices are zero based, index access (column,row) for (SCSIZE row = 0; row PutDouble( pMatA->GetDouble(col,row)/fScale, col, row); fEuclid = lcl_TGetColumnEuclideanNorm(pMatA, row, row, nN); fFactor = 1.0/fEuclid/(fEuclid + fabs(pMatA->GetDouble(row,row))); 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; rPutDouble( 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 Vektor u, from which the Householder transformation is build, * 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_CaluclateQRdecomposition. */ void lcl_ApplyHouseholderTransformation(ScMatrixRef pMatA, SCSIZE nC, 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(ScMatrixRef pMatA, SCSIZE nR, 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(ScMatrixRef pMatA, ::std::vector< double>& pVecR, ScMatrixRef pMatS, SCSIZE nK, bool bIsTransposed) { // ScMatrix matrices are zero based, index access (column,row) double fSum; SCSIZE row; // SCSIZE is never negative, therefore test with rowp1=row+1 for (SCSIZE rowp1 = nK; rowp1>0; rowp1--) { row = rowp1-1; fSum = pMatS->GetDouble(row); for (SCSIZE col = rowp1; colGetDouble(row,col) * pMatS->GetDouble(col); else fSum -= pMatA->GetDouble(col,row) * pMatS->GetDouble(col); pMatS->PutDouble( fSum / 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(ScMatrixRef pMatA, ::std::vector< double>& pVecR, ScMatrixRef pMatT, SCSIZE nK, bool bIsTransposed) { // ScMatrix matrices are zero based, index access (column,row) double fSum; for (SCSIZE row = 0; row < nK; row++) { 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 / pVecR[row] , row); } } /* Calculates Z = R * B * R is given in matrix A and vector VecR as obtained from the QR * decompostion 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(ScMatrixRef pMatA, ::std::vector< double>& pVecR, ScMatrixRef pMatB, ScMatrixRef pMatZ, SCSIZE nK, bool bIsTransposed) { // ScMatrix matrices are zero based, index access (column,row) double fSum; for (SCSIZE row = 0; row < nK; row++) { 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, row); } } double lcl_GetMeanOverAll(ScMatrixRef pMat, SCSIZE nN) { double fSum = 0.0; for (SCSIZE i=0 ; iGetDouble(i); return fSum/static_cast(nN); } // Calculates means of the columns of matrix X. X is a RxC matrix; // ResMat is a 1xC matrix (=row). void lcl_CalculateColumnMeans(ScMatrixRef pX, ScMatrixRef pResMat, SCSIZE nC, SCSIZE nR) { double fSum = 0.0; for (SCSIZE i=0; i < nC; i++) { fSum =0.0; for (SCSIZE k=0; k < nR; k++) fSum += pX->GetDouble(i,k); // GetDouble(Column,Row) pResMat ->PutDouble( fSum/static_cast(nR),i); } } // Calculates means of the rows of matrix X. X is a RxC matrix; // ResMat is a Rx1 matrix (=column). void lcl_CalculateRowMeans(ScMatrixRef pX, ScMatrixRef pResMat, SCSIZE nC, SCSIZE nR) { double fSum = 0.0; for (SCSIZE k=0; k < nR; k++) { fSum =0.0; for (SCSIZE i=0; i < nC; i++) fSum += pX->GetDouble(i,k); // GetDouble(Column,Row) pResMat ->PutDouble( fSum/static_cast(nC),k); } } void lcl_CalculateColumnsDelta(ScMatrixRef pMat, 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(ScMatrixRef pMat, 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(ScMatrixRef pMatX, ScMatrixRef pMatY, double fSlope, SCSIZE nN) { double fSum = 0.0; double fTemp = 0.0; for (SCSIZE i=0; iGetDouble(i) - fSlope * pMatX->GetDouble(i); fSum += fTemp * fTemp; } return fSum; } // 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); nCX = nCY; nRX = nRY; if (!pMatX) { PushIllegalArgument(); return false; } for ( SCSIZE i = 1; i <= nCountY; i++ ) pMatX->PutDouble(static_cast(i), i-1); nCase = 1; N = nCountY; M = 1; } return true; } // ----------------------------------------------------------------------------- // LINEST void ScInterpreter::ScRGP() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRGP" ); CalulateRGPRKP(false); } // LOGEST void ScInterpreter::ScRKP() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRKP" ); CalulateRGPRKP(true); } void ScInterpreter::CalulateRGPRKP(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 = NULL; } else { pMatX = GetMatrix(); } } else pMatX = NULL; ScMatrixRef pMatY; 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 && (NPutString(ScGlobal::GetRscString(STR_NV_STR), i, 2 ); pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 3 ); pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 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(errCodeOverflow); return; } pMatX = pNewX; pMatY = pNewY; // DeltaY is possible here; DeltaX depends on nCase, so later fMeanY = lcl_GetMeanOverAll(pMatY, N); for (SCSIZE i=0; iPutDouble( ::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; iPutDouble( ::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( (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->PutString(ScGlobal::GetRscString(STR_NV_STR), 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->PutString(ScGlobal::GetRscString(STR_NV_STR), 1, 1); pResMat->PutDouble(1.0, 0, 2); // R^2 } else { double fFstatistic = (fSSreg / static_cast(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(N)); pResMat->PutDouble(fSigmaIntercept, 1, 1); } else { pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 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); // 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); // from b1 to bK if (!pMeans || !pMatZ || !pSlopes) { PushError(errCodeOverflow); 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 = 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( (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->PutString(ScGlobal::GetRscString(STR_NV_STR), 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; iPutDouble(0.0, K-1-i, 1); // SigmaIntercept = RMSE * sqrt(...) = 0 if (bConstant) pResMat->PutDouble(0.0, K, 1); //SigmaIntercept else pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1); // R^2 = SSreg / (SSreg + SSresid) = 1.0 pResMat->PutDouble(1.0, 0, 2); // R^2 } else { double fFstatistic = (fSSreg / static_cast(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. double fSigmaSlope = 0.0; double fSigmaIntercept = 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) 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); fSigmaIntercept += fPart * pMeans->GetDouble(col); } } if (bConstant) { fSigmaIntercept = fRMSE * sqrt(fSigmaIntercept + 1.0 / static_cast(N)); pResMat->PutDouble(fSigmaIntercept, K, 1); } else { pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 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); // 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); // from b1 to bK if (!pMeans || !pMatZ || !pSlopes) { PushError(errCodeOverflow); 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 = 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( (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->PutString(ScGlobal::GetRscString(STR_NV_STR), 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; iPutDouble(0.0, K-1-i, 1); // SigmaIntercept = RMSE * sqrt(...) = 0 if (bConstant) pResMat->PutDouble(0.0, K, 1); //SigmaIntercept else pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1); // R^2 = SSreg / (SSreg + SSresid) = 1.0 pResMat->PutDouble(1.0, 0, 2); // R^2 } else { double fFstatistic = (fSSreg / static_cast(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 double fSigmaSlope = 0.0; double fSigmaIntercept = 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) fSigmaSlope = fRMSE * sqrt(pMatZ->GetDouble(row)); pResMat->PutDouble(fSigmaSlope, K-1-row, 1); if (bConstant) { fPart = lcl_GetSumProduct(pMeans, pMatZ, K); fSigmaIntercept += fPart * pMeans->GetDouble(row); } } if (bConstant) { fSigmaIntercept = fRMSE * sqrt(fSigmaIntercept + 1.0 / static_cast(N)); pResMat->PutDouble(fSigmaIntercept, K, 1); } else { pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1); } double fR2 = fSSreg / (fSSreg + fSSresid); pResMat->PutDouble(fR2, 0, 2); } } PushMatrix(pResMat); } } return; } void ScInterpreter::ScTrend() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScTrend" ); CalculateTrendGrowth(false); } void ScInterpreter::ScGrowth() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScGrowth" ); CalculateTrendGrowth(true); } void ScInterpreter::CalculateTrendGrowth(bool _bGrowth) { sal_uInt8 nParamCount = GetByte(); if ( !MustHaveParamCount( nParamCount, 1, 4 ) ) return; // optional fourth parameter bool bConstant; if (nParamCount == 4) bConstant = GetBool(); else bConstant = true; // The third parameter may be missing in ODF, although the fourth parameter // is present. Default values depend on data not yet read. ScMatrixRef pMatNewX; if (nParamCount >= 3) { if (IsMissing()) { Pop(); pMatNewX = NULL; } else pMatNewX = GetMatrix(); } else pMatNewX = NULL; // 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 = NULL; } else { pMatX = GetMatrix(); } } else pMatX = NULL; ScMatrixRef pMatY; 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 && (NClone(); // 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); else { if (nCase==2) pResMat = GetNewMat(1,nRXN); else pResMat = GetNewMat(nCXN,1); } if (!pResMat) { PushError(errCodeOverflow); 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(errStackOverflow); return; } pMatX = pCopyX; pMatY = pCopyY; // DeltaY is possible here; DeltaX depends on nCase, so later fMeanY = lcl_GetMeanOverAll(pMatY, N); for (SCSIZE i=0; iPutDouble( ::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; iPutDouble( ::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; double fHelp; if (bConstant) { 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); // mean of each column ScMatrixRef pSlopes = GetNewMat(1,K); // from b1 to bK if (!pMeans || !pSlopes) { PushError(errCodeOverflow); 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 = 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); // mean of each row ScMatrixRef pSlopes = GetNewMat(K,1); // row from b1 to bK if (!pMeans || !pSlopes) { PushError(errCodeOverflow); 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 = 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); return; } void ScInterpreter::ScMatRef() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatRef" ); // Falls Deltarefs drin sind... Push( (FormulaToken&)*pCur ); ScAddress aAdr; PopSingleRef( aAdr ); ScFormulaCell* pCell = (ScFormulaCell*) GetCell( aAdr ); if( pCell && pCell->GetCellType() == CELLTYPE_FORMULA ) { const ScMatrix* pMat = pCell->GetMatrix(); if( pMat ) { SCSIZE nCols, nRows; pMat->GetDimensions( nCols, nRows ); SCSIZE nC = static_cast(aPos.Col() - aAdr.Col()); SCSIZE nR = static_cast(aPos.Row() - aAdr.Row()); if ((nCols <= nC && nCols != 1) || (nRows <= nR && nRows != 1)) PushNA(); else { ScMatValType nMatValType; const ScMatrixValue* pMatVal = pMat->Get( nC, nR, nMatValType); if (ScMatrix::IsNonValueType( nMatValType)) { if (ScMatrix::IsEmptyPathType( nMatValType)) { // result of empty sal_False jump path nFuncFmtType = NUMBERFORMAT_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( pMatVal->GetString() ); } else { PushDouble(pMatVal->fVal); // handles DoubleError pDok->GetNumberFormatInfo( nCurFmtType, nCurFmtIndex, aAdr, pCell ); nFuncFmtType = nCurFmtType; nFuncFmtIndex = nCurFmtIndex; } } } else { // If not a result matrix, obtain the cell value. sal_uInt16 nErr = pCell->GetErrCode(); if (nErr) PushError( nErr ); else if( pCell->IsValue() ) PushDouble( pCell->GetValue() ); else { String aVal; pCell->GetString( aVal ); PushString( aVal ); } pDok->GetNumberFormatInfo( nCurFmtType, nCurFmtIndex, aAdr, pCell ); nFuncFmtType = nCurFmtType; nFuncFmtIndex = nCurFmtIndex; } } else PushError( errNoRef ); } void ScInterpreter::ScInfo() { RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScInfo" ); if( MustHaveParamCount( GetByte(), 1 ) ) { String aStr = GetString(); ScCellKeywordTranslator::transKeyword(aStr, ScGlobal::GetLocale(), ocInfo); if( aStr.EqualsAscii( "SYSTEM" ) ) PushString( String( RTL_CONSTASCII_USTRINGPARAM( SC_INFO_OSVERSION ) ) ); else if( aStr.EqualsAscii( "OSVERSION" ) ) PushString( String( RTL_CONSTASCII_USTRINGPARAM( "Windows (32-bit) NT 5.01" ) ) ); else if( aStr.EqualsAscii( "RELEASE" ) ) PushString( ::utl::Bootstrap::getBuildIdData( ::rtl::OUString() ) ); else if( aStr.EqualsAscii( "NUMFILE" ) ) PushDouble( 1 ); else if( aStr.EqualsAscii( "RECALC" ) ) PushString( ScGlobal::GetRscString( pDok->GetAutoCalc() ? STR_RECALC_AUTO : STR_RECALC_MANUAL ) ); else PushIllegalArgument(); } }