1b3f79822SAndrew Rist /**************************************************************
2cdf0e10cSrcweir *
3b3f79822SAndrew Rist * Licensed to the Apache Software Foundation (ASF) under one
4b3f79822SAndrew Rist * or more contributor license agreements. See the NOTICE file
5b3f79822SAndrew Rist * distributed with this work for additional information
6b3f79822SAndrew Rist * regarding copyright ownership. The ASF licenses this file
7b3f79822SAndrew Rist * to you under the Apache License, Version 2.0 (the
8b3f79822SAndrew Rist * "License"); you may not use this file except in compliance
9b3f79822SAndrew Rist * with the License. You may obtain a copy of the License at
10b3f79822SAndrew Rist *
11b3f79822SAndrew Rist * http://www.apache.org/licenses/LICENSE-2.0
12b3f79822SAndrew Rist *
13b3f79822SAndrew Rist * Unless required by applicable law or agreed to in writing,
14b3f79822SAndrew Rist * software distributed under the License is distributed on an
15b3f79822SAndrew Rist * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
16b3f79822SAndrew Rist * KIND, either express or implied. See the License for the
17b3f79822SAndrew Rist * specific language governing permissions and limitations
18b3f79822SAndrew Rist * under the License.
19b3f79822SAndrew Rist *
20b3f79822SAndrew Rist *************************************************************/
21b3f79822SAndrew Rist
22b3f79822SAndrew Rist
23cdf0e10cSrcweir
24cdf0e10cSrcweir // MARKER(update_precomp.py): autogen include statement, do not remove
25cdf0e10cSrcweir #include "precompiled_sc.hxx"
26cdf0e10cSrcweir
27cdf0e10cSrcweir // INCLUDE ---------------------------------------------------------------
28cdf0e10cSrcweir
29cdf0e10cSrcweir #ifndef INCLUDED_RTL_MATH_HXX
30cdf0e10cSrcweir #include <rtl/math.hxx>
31cdf0e10cSrcweir #endif
32cdf0e10cSrcweir #include <rtl/logfile.hxx>
33cdf0e10cSrcweir #include <string.h>
34cdf0e10cSrcweir #include <math.h>
35cdf0e10cSrcweir #include <stdio.h>
36cdf0e10cSrcweir
37cdf0e10cSrcweir #if OSL_DEBUG_LEVEL > 1
38cdf0e10cSrcweir #include <stdio.h>
39cdf0e10cSrcweir #endif
40cdf0e10cSrcweir #include <unotools/bootstrap.hxx>
41cdf0e10cSrcweir #include <svl/zforlist.hxx>
42cdf0e10cSrcweir
43cdf0e10cSrcweir #include "interpre.hxx"
44cdf0e10cSrcweir #include "global.hxx"
45cdf0e10cSrcweir #include "compiler.hxx"
46cdf0e10cSrcweir #include "cell.hxx"
47cdf0e10cSrcweir #include "document.hxx"
48cdf0e10cSrcweir #include "dociter.hxx"
49cdf0e10cSrcweir #include "scmatrix.hxx"
50cdf0e10cSrcweir #include "globstr.hrc"
51cdf0e10cSrcweir #include "cellkeytranslator.hxx"
52cdf0e10cSrcweir
53cdf0e10cSrcweir #include <string.h>
54cdf0e10cSrcweir #include <math.h>
55cdf0e10cSrcweir #include <vector>
56cdf0e10cSrcweir
57cdf0e10cSrcweir using ::std::vector;
58cdf0e10cSrcweir using namespace formula;
59cdf0e10cSrcweir
60cdf0e10cSrcweir const double fInvEpsilon = 1.0E-7;
61cdf0e10cSrcweir
62cdf0e10cSrcweir // -----------------------------------------------------------------------
63cdf0e10cSrcweir struct MatrixAdd : public ::std::binary_function<double,double,double>
64cdf0e10cSrcweir {
operator ()MatrixAdd65cdf0e10cSrcweir inline double operator() (const double& lhs, const double& rhs) const
66cdf0e10cSrcweir {
67cdf0e10cSrcweir return ::rtl::math::approxAdd( lhs,rhs);
68cdf0e10cSrcweir }
69cdf0e10cSrcweir };
70cdf0e10cSrcweir struct MatrixSub : public ::std::binary_function<double,double,double>
71cdf0e10cSrcweir {
operator ()MatrixSub72cdf0e10cSrcweir inline double operator() (const double& lhs, const double& rhs) const
73cdf0e10cSrcweir {
74cdf0e10cSrcweir return ::rtl::math::approxSub( lhs,rhs);
75cdf0e10cSrcweir }
76cdf0e10cSrcweir };
77cdf0e10cSrcweir struct MatrixMul : public ::std::binary_function<double,double,double>
78cdf0e10cSrcweir {
operator ()MatrixMul79cdf0e10cSrcweir inline double operator() (const double& lhs, const double& rhs) const
80cdf0e10cSrcweir {
81cdf0e10cSrcweir return lhs * rhs;
82cdf0e10cSrcweir }
83cdf0e10cSrcweir };
84cdf0e10cSrcweir struct MatrixDiv : public ::std::binary_function<double,double,double>
85cdf0e10cSrcweir {
operator ()MatrixDiv86cdf0e10cSrcweir inline double operator() (const double& lhs, const double& rhs) const
87cdf0e10cSrcweir {
88cdf0e10cSrcweir return ScInterpreter::div( lhs,rhs);
89cdf0e10cSrcweir }
90cdf0e10cSrcweir };
91cdf0e10cSrcweir struct MatrixPow : public ::std::binary_function<double,double,double>
92cdf0e10cSrcweir {
operator ()MatrixPow93cdf0e10cSrcweir inline double operator() (const double& lhs, const double& rhs) const
94cdf0e10cSrcweir {
95cdf0e10cSrcweir return ::pow( lhs,rhs);
96cdf0e10cSrcweir }
97cdf0e10cSrcweir };
98cdf0e10cSrcweir
ScGetGCD(double fx,double fy)99cdf0e10cSrcweir double ScInterpreter::ScGetGCD(double fx, double fy)
100cdf0e10cSrcweir {
101cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::div" );
102cdf0e10cSrcweir // By ODFF definition GCD(0,a) => a. This is also vital for the code in
103cdf0e10cSrcweir // ScGCD() to work correctly with a preset fy=0.0
104cdf0e10cSrcweir if (fy == 0.0)
105cdf0e10cSrcweir return fx;
106cdf0e10cSrcweir else if (fx == 0.0)
107cdf0e10cSrcweir return fy;
108cdf0e10cSrcweir else
109cdf0e10cSrcweir {
110cdf0e10cSrcweir double fz = fmod(fx, fy);
111cdf0e10cSrcweir while (fz > 0.0)
112cdf0e10cSrcweir {
113cdf0e10cSrcweir fx = fy;
114cdf0e10cSrcweir fy = fz;
115cdf0e10cSrcweir fz = fmod(fx, fy);
116cdf0e10cSrcweir }
117cdf0e10cSrcweir return fy;
118cdf0e10cSrcweir }
119cdf0e10cSrcweir }
120cdf0e10cSrcweir
ScGCD()121cdf0e10cSrcweir void ScInterpreter::ScGCD()
122cdf0e10cSrcweir {
123cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScGCD" );
124cdf0e10cSrcweir short nParamCount = GetByte();
125cdf0e10cSrcweir if ( MustHaveParamCountMin( nParamCount, 1 ) )
126cdf0e10cSrcweir {
127cdf0e10cSrcweir double fx, fy = 0.0;
128cdf0e10cSrcweir ScRange aRange;
129cdf0e10cSrcweir size_t nRefInList = 0;
130cdf0e10cSrcweir while (!nGlobalError && nParamCount-- > 0)
131cdf0e10cSrcweir {
132cdf0e10cSrcweir switch (GetStackType())
133cdf0e10cSrcweir {
134cdf0e10cSrcweir case svDouble :
135cdf0e10cSrcweir case svString:
136cdf0e10cSrcweir case svSingleRef:
137cdf0e10cSrcweir {
138cdf0e10cSrcweir fx = ::rtl::math::approxFloor( GetDouble());
139cdf0e10cSrcweir if (fx < 0.0)
140cdf0e10cSrcweir {
141cdf0e10cSrcweir PushIllegalArgument();
142cdf0e10cSrcweir return;
143cdf0e10cSrcweir }
144cdf0e10cSrcweir fy = ScGetGCD(fx, fy);
145cdf0e10cSrcweir }
146cdf0e10cSrcweir break;
147cdf0e10cSrcweir case svDoubleRef :
148cdf0e10cSrcweir case svRefList :
149cdf0e10cSrcweir {
150cdf0e10cSrcweir sal_uInt16 nErr = 0;
151cdf0e10cSrcweir PopDoubleRef( aRange, nParamCount, nRefInList);
152cdf0e10cSrcweir double nCellVal;
153cdf0e10cSrcweir ScValueIterator aValIter(pDok, aRange, glSubTotal);
154cdf0e10cSrcweir if (aValIter.GetFirst(nCellVal, nErr))
155cdf0e10cSrcweir {
156cdf0e10cSrcweir do
157cdf0e10cSrcweir {
158cdf0e10cSrcweir fx = ::rtl::math::approxFloor( nCellVal);
159cdf0e10cSrcweir if (fx < 0.0)
160cdf0e10cSrcweir {
161cdf0e10cSrcweir PushIllegalArgument();
162cdf0e10cSrcweir return;
163cdf0e10cSrcweir }
164cdf0e10cSrcweir fy = ScGetGCD(fx, fy);
165cdf0e10cSrcweir } while (nErr == 0 && aValIter.GetNext(nCellVal, nErr));
166cdf0e10cSrcweir }
167cdf0e10cSrcweir SetError(nErr);
168cdf0e10cSrcweir }
169cdf0e10cSrcweir break;
170cdf0e10cSrcweir case svMatrix :
171cdf0e10cSrcweir {
172cdf0e10cSrcweir ScMatrixRef pMat = PopMatrix();
173cdf0e10cSrcweir if (pMat)
174cdf0e10cSrcweir {
175cdf0e10cSrcweir SCSIZE nC, nR;
176cdf0e10cSrcweir pMat->GetDimensions(nC, nR);
177cdf0e10cSrcweir if (nC == 0 || nR == 0)
178cdf0e10cSrcweir SetError(errIllegalArgument);
179cdf0e10cSrcweir else
180cdf0e10cSrcweir {
181cdf0e10cSrcweir SCSIZE nCount = nC * nR;
182cdf0e10cSrcweir for ( SCSIZE j = 0; j < nCount; j++ )
183cdf0e10cSrcweir {
184cdf0e10cSrcweir if (!pMat->IsValue(j))
185cdf0e10cSrcweir {
186cdf0e10cSrcweir PushIllegalArgument();
187cdf0e10cSrcweir return;
188cdf0e10cSrcweir }
189cdf0e10cSrcweir fx = ::rtl::math::approxFloor( pMat->GetDouble(j));
190cdf0e10cSrcweir if (fx < 0.0)
191cdf0e10cSrcweir {
192cdf0e10cSrcweir PushIllegalArgument();
193cdf0e10cSrcweir return;
194cdf0e10cSrcweir }
195cdf0e10cSrcweir fy = ScGetGCD(fx, fy);
196cdf0e10cSrcweir }
197cdf0e10cSrcweir }
198cdf0e10cSrcweir }
199cdf0e10cSrcweir }
200cdf0e10cSrcweir break;
201cdf0e10cSrcweir default : SetError(errIllegalParameter); break;
202cdf0e10cSrcweir }
203cdf0e10cSrcweir }
204cdf0e10cSrcweir PushDouble(fy);
205cdf0e10cSrcweir }
206cdf0e10cSrcweir }
207cdf0e10cSrcweir
ScLCM()208cdf0e10cSrcweir void ScInterpreter:: ScLCM()
209cdf0e10cSrcweir {
210cdf0e10cSrcweir short nParamCount = GetByte();
211cdf0e10cSrcweir if ( MustHaveParamCountMin( nParamCount, 1 ) )
212cdf0e10cSrcweir {
213cdf0e10cSrcweir double fx, fy = 1.0;
214cdf0e10cSrcweir ScRange aRange;
215cdf0e10cSrcweir size_t nRefInList = 0;
216cdf0e10cSrcweir while (!nGlobalError && nParamCount-- > 0)
217cdf0e10cSrcweir {
218cdf0e10cSrcweir switch (GetStackType())
219cdf0e10cSrcweir {
220cdf0e10cSrcweir case svDouble :
221cdf0e10cSrcweir case svString:
222cdf0e10cSrcweir case svSingleRef:
223cdf0e10cSrcweir {
224cdf0e10cSrcweir fx = ::rtl::math::approxFloor( GetDouble());
225cdf0e10cSrcweir if (fx < 0.0)
226cdf0e10cSrcweir {
227cdf0e10cSrcweir PushIllegalArgument();
228cdf0e10cSrcweir return;
229cdf0e10cSrcweir }
230cdf0e10cSrcweir if (fx == 0.0 || fy == 0.0)
231cdf0e10cSrcweir fy = 0.0;
232cdf0e10cSrcweir else
233cdf0e10cSrcweir fy = fx * fy / ScGetGCD(fx, fy);
234cdf0e10cSrcweir }
235cdf0e10cSrcweir break;
236cdf0e10cSrcweir case svDoubleRef :
237cdf0e10cSrcweir case svRefList :
238cdf0e10cSrcweir {
239cdf0e10cSrcweir sal_uInt16 nErr = 0;
240cdf0e10cSrcweir PopDoubleRef( aRange, nParamCount, nRefInList);
241cdf0e10cSrcweir double nCellVal;
242cdf0e10cSrcweir ScValueIterator aValIter(pDok, aRange, glSubTotal);
243cdf0e10cSrcweir if (aValIter.GetFirst(nCellVal, nErr))
244cdf0e10cSrcweir {
245cdf0e10cSrcweir do
246cdf0e10cSrcweir {
247cdf0e10cSrcweir fx = ::rtl::math::approxFloor( nCellVal);
248cdf0e10cSrcweir if (fx < 0.0)
249cdf0e10cSrcweir {
250cdf0e10cSrcweir PushIllegalArgument();
251cdf0e10cSrcweir return;
252cdf0e10cSrcweir }
253cdf0e10cSrcweir if (fx == 0.0 || fy == 0.0)
254cdf0e10cSrcweir fy = 0.0;
255cdf0e10cSrcweir else
256cdf0e10cSrcweir fy = fx * fy / ScGetGCD(fx, fy);
257cdf0e10cSrcweir } while (nErr == 0 && aValIter.GetNext(nCellVal, nErr));
258cdf0e10cSrcweir }
259cdf0e10cSrcweir SetError(nErr);
260cdf0e10cSrcweir }
261cdf0e10cSrcweir break;
262cdf0e10cSrcweir case svMatrix :
263cdf0e10cSrcweir {
264cdf0e10cSrcweir ScMatrixRef pMat = PopMatrix();
265cdf0e10cSrcweir if (pMat)
266cdf0e10cSrcweir {
267cdf0e10cSrcweir SCSIZE nC, nR;
268cdf0e10cSrcweir pMat->GetDimensions(nC, nR);
269cdf0e10cSrcweir if (nC == 0 || nR == 0)
270cdf0e10cSrcweir SetError(errIllegalArgument);
271cdf0e10cSrcweir else
272cdf0e10cSrcweir {
273cdf0e10cSrcweir SCSIZE nCount = nC * nR;
274cdf0e10cSrcweir for ( SCSIZE j = 0; j < nCount; j++ )
275cdf0e10cSrcweir {
276cdf0e10cSrcweir if (!pMat->IsValue(j))
277cdf0e10cSrcweir {
278cdf0e10cSrcweir PushIllegalArgument();
279cdf0e10cSrcweir return;
280cdf0e10cSrcweir }
281cdf0e10cSrcweir fx = ::rtl::math::approxFloor( pMat->GetDouble(j));
282cdf0e10cSrcweir if (fx < 0.0)
283cdf0e10cSrcweir {
284cdf0e10cSrcweir PushIllegalArgument();
285cdf0e10cSrcweir return;
286cdf0e10cSrcweir }
287cdf0e10cSrcweir if (fx == 0.0 || fy == 0.0)
288cdf0e10cSrcweir fy = 0.0;
289cdf0e10cSrcweir else
290cdf0e10cSrcweir fy = fx * fy / ScGetGCD(fx, fy);
291cdf0e10cSrcweir }
292cdf0e10cSrcweir }
293cdf0e10cSrcweir }
294cdf0e10cSrcweir }
295cdf0e10cSrcweir break;
296cdf0e10cSrcweir default : SetError(errIllegalParameter); break;
297cdf0e10cSrcweir }
298cdf0e10cSrcweir }
299cdf0e10cSrcweir PushDouble(fy);
300cdf0e10cSrcweir }
301cdf0e10cSrcweir }
302cdf0e10cSrcweir
GetNewMat(SCSIZE nC,SCSIZE nR)303cdf0e10cSrcweir ScMatrixRef ScInterpreter::GetNewMat(SCSIZE nC, SCSIZE nR)
304cdf0e10cSrcweir {
305cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetNewMat" );
306cdf0e10cSrcweir ScMatrix* pMat = new ScMatrix( nC, nR);
307cdf0e10cSrcweir pMat->SetErrorInterpreter( this);
308cdf0e10cSrcweir // A temporary matrix is mutable and ScMatrix::CloneIfConst() returns the
309cdf0e10cSrcweir // very matrix.
310cdf0e10cSrcweir pMat->SetImmutable( false);
311cdf0e10cSrcweir SCSIZE nCols, nRows;
312cdf0e10cSrcweir pMat->GetDimensions( nCols, nRows);
313cdf0e10cSrcweir if ( nCols != nC || nRows != nR )
314c2202f61Smseidel { // arbitrary limit of elements exceeded
315cdf0e10cSrcweir SetError( errStackOverflow);
316cdf0e10cSrcweir pMat->Delete();
317cdf0e10cSrcweir pMat = NULL;
318cdf0e10cSrcweir }
319cdf0e10cSrcweir return pMat;
320cdf0e10cSrcweir }
321cdf0e10cSrcweir
CreateMatrixFromDoubleRef(const FormulaToken * pToken,SCCOL nCol1,SCROW nRow1,SCTAB nTab1,SCCOL nCol2,SCROW nRow2,SCTAB nTab2)322cdf0e10cSrcweir ScMatrixRef ScInterpreter::CreateMatrixFromDoubleRef( const FormulaToken* pToken,
323cdf0e10cSrcweir SCCOL nCol1, SCROW nRow1, SCTAB nTab1,
324cdf0e10cSrcweir SCCOL nCol2, SCROW nRow2, SCTAB nTab2 )
325cdf0e10cSrcweir {
326cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CreateMatrixFromDoubleRef" );
327cdf0e10cSrcweir ScMatrixRef pMat = NULL;
328cdf0e10cSrcweir if (nTab1 == nTab2 && !nGlobalError)
329cdf0e10cSrcweir {
330cdf0e10cSrcweir ScTokenMatrixMap::const_iterator aIter;
331cdf0e10cSrcweir if ( static_cast<SCSIZE>(nRow2 - nRow1 + 1) *
332cdf0e10cSrcweir static_cast<SCSIZE>(nCol2 - nCol1 + 1) >
333cdf0e10cSrcweir ScMatrix::GetElementsMax() )
334cdf0e10cSrcweir SetError(errStackOverflow);
335cdf0e10cSrcweir else if (pTokenMatrixMap && ((aIter = pTokenMatrixMap->find( pToken))
336cdf0e10cSrcweir != pTokenMatrixMap->end()))
337cdf0e10cSrcweir pMat = static_cast<ScToken*>((*aIter).second.get())->GetMatrix();
338cdf0e10cSrcweir else
339cdf0e10cSrcweir {
340cdf0e10cSrcweir SCSIZE nMatCols = static_cast<SCSIZE>(nCol2 - nCol1 + 1);
341cdf0e10cSrcweir SCSIZE nMatRows = static_cast<SCSIZE>(nRow2 - nRow1 + 1);
342cdf0e10cSrcweir pMat = GetNewMat( nMatCols, nMatRows);
343cdf0e10cSrcweir if (pMat && !nGlobalError)
344cdf0e10cSrcweir {
345cdf0e10cSrcweir // Set position where the next entry is expected.
346cdf0e10cSrcweir SCROW nNextRow = nRow1;
347cdf0e10cSrcweir SCCOL nNextCol = nCol1;
348cdf0e10cSrcweir // Set last position as if there was a previous entry.
349cdf0e10cSrcweir SCROW nThisRow = nRow2;
350cdf0e10cSrcweir SCCOL nThisCol = nCol1 - 1;
351cdf0e10cSrcweir ScCellIterator aCellIter( pDok, nCol1, nRow1, nTab1, nCol2,
352cdf0e10cSrcweir nRow2, nTab2);
353cdf0e10cSrcweir for (ScBaseCell* pCell = aCellIter.GetFirst(); pCell; pCell =
354cdf0e10cSrcweir aCellIter.GetNext())
355cdf0e10cSrcweir {
356cdf0e10cSrcweir nThisCol = aCellIter.GetCol();
357cdf0e10cSrcweir nThisRow = aCellIter.GetRow();
358cdf0e10cSrcweir if (nThisCol != nNextCol || nThisRow != nNextRow)
359cdf0e10cSrcweir {
360cdf0e10cSrcweir // Fill empty between iterator's positions.
361cdf0e10cSrcweir for ( ; nNextCol <= nThisCol; ++nNextCol)
362cdf0e10cSrcweir {
363cdf0e10cSrcweir SCSIZE nC = nNextCol - nCol1;
364cdf0e10cSrcweir SCSIZE nMatStopRow = ((nNextCol < nThisCol) ?
365cdf0e10cSrcweir nMatRows : nThisRow - nRow1);
366cdf0e10cSrcweir for (SCSIZE nR = nNextRow - nRow1; nR <
367cdf0e10cSrcweir nMatStopRow; ++nR)
368cdf0e10cSrcweir {
369cdf0e10cSrcweir pMat->PutEmpty( nC, nR);
370cdf0e10cSrcweir }
371cdf0e10cSrcweir nNextRow = nRow1;
372cdf0e10cSrcweir }
373cdf0e10cSrcweir }
374cdf0e10cSrcweir if (nThisRow == nRow2)
375cdf0e10cSrcweir {
376cdf0e10cSrcweir nNextCol = nThisCol + 1;
377cdf0e10cSrcweir nNextRow = nRow1;
378cdf0e10cSrcweir }
379cdf0e10cSrcweir else
380cdf0e10cSrcweir {
381cdf0e10cSrcweir nNextCol = nThisCol;
382cdf0e10cSrcweir nNextRow = nThisRow + 1;
383cdf0e10cSrcweir }
384cdf0e10cSrcweir if (HasCellEmptyData(pCell))
385cdf0e10cSrcweir {
386cdf0e10cSrcweir pMat->PutEmpty( static_cast<SCSIZE>(nThisCol-nCol1),
387cdf0e10cSrcweir static_cast<SCSIZE>(nThisRow-nRow1));
388cdf0e10cSrcweir }
389cdf0e10cSrcweir else if (HasCellValueData(pCell))
390cdf0e10cSrcweir {
391cdf0e10cSrcweir ScAddress aAdr( nThisCol, nThisRow, nTab1);
392cdf0e10cSrcweir double fVal = GetCellValue( aAdr, pCell);
393cdf0e10cSrcweir if ( nGlobalError )
394cdf0e10cSrcweir {
395cdf0e10cSrcweir fVal = CreateDoubleError( nGlobalError);
396cdf0e10cSrcweir nGlobalError = 0;
397cdf0e10cSrcweir }
398cdf0e10cSrcweir pMat->PutDouble( fVal,
399cdf0e10cSrcweir static_cast<SCSIZE>(nThisCol-nCol1),
400cdf0e10cSrcweir static_cast<SCSIZE>(nThisRow-nRow1));
401cdf0e10cSrcweir }
402cdf0e10cSrcweir else
403cdf0e10cSrcweir {
404cdf0e10cSrcweir String aStr;
405cdf0e10cSrcweir GetCellString( aStr, pCell);
406cdf0e10cSrcweir if ( nGlobalError )
407cdf0e10cSrcweir {
408cdf0e10cSrcweir double fVal = CreateDoubleError( nGlobalError);
409cdf0e10cSrcweir nGlobalError = 0;
410cdf0e10cSrcweir pMat->PutDouble( fVal,
411cdf0e10cSrcweir static_cast<SCSIZE>(nThisCol-nCol1),
412cdf0e10cSrcweir static_cast<SCSIZE>(nThisRow-nRow1));
413cdf0e10cSrcweir }
414cdf0e10cSrcweir else
415cdf0e10cSrcweir pMat->PutString( aStr,
416cdf0e10cSrcweir static_cast<SCSIZE>(nThisCol-nCol1),
417cdf0e10cSrcweir static_cast<SCSIZE>(nThisRow-nRow1));
418cdf0e10cSrcweir }
419cdf0e10cSrcweir }
420cdf0e10cSrcweir // Fill empty if iterator's last position wasn't the end.
421cdf0e10cSrcweir if (nThisCol != nCol2 || nThisRow != nRow2)
422cdf0e10cSrcweir {
423cdf0e10cSrcweir for ( ; nNextCol <= nCol2; ++nNextCol)
424cdf0e10cSrcweir {
425cdf0e10cSrcweir SCSIZE nC = nNextCol - nCol1;
426cdf0e10cSrcweir for (SCSIZE nR = nNextRow - nRow1; nR < nMatRows; ++nR)
427cdf0e10cSrcweir {
428cdf0e10cSrcweir pMat->PutEmpty( nC, nR);
429cdf0e10cSrcweir }
430cdf0e10cSrcweir nNextRow = nRow1;
431cdf0e10cSrcweir }
432cdf0e10cSrcweir }
433cdf0e10cSrcweir if (pTokenMatrixMap)
434cdf0e10cSrcweir pTokenMatrixMap->insert( ScTokenMatrixMap::value_type(
435cdf0e10cSrcweir pToken, new ScMatrixToken( pMat)));
436cdf0e10cSrcweir }
437cdf0e10cSrcweir }
438cdf0e10cSrcweir }
439cdf0e10cSrcweir else // not a 2D matrix
440cdf0e10cSrcweir SetError(errIllegalParameter);
441cdf0e10cSrcweir return pMat;
442cdf0e10cSrcweir }
443cdf0e10cSrcweir
444cdf0e10cSrcweir
GetMatrix()445cdf0e10cSrcweir ScMatrixRef ScInterpreter::GetMatrix()
446cdf0e10cSrcweir {
447cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::GetMatrix" );
448cdf0e10cSrcweir ScMatrixRef pMat = NULL;
449cdf0e10cSrcweir switch (GetRawStackType())
450cdf0e10cSrcweir {
451cdf0e10cSrcweir case svSingleRef :
452cdf0e10cSrcweir {
453cdf0e10cSrcweir ScAddress aAdr;
454cdf0e10cSrcweir PopSingleRef( aAdr );
455cdf0e10cSrcweir pMat = GetNewMat(1, 1);
456cdf0e10cSrcweir if (pMat)
457cdf0e10cSrcweir {
458cdf0e10cSrcweir ScBaseCell* pCell = GetCell( aAdr );
459cdf0e10cSrcweir if (HasCellEmptyData(pCell))
460cdf0e10cSrcweir pMat->PutEmpty( 0 );
461cdf0e10cSrcweir else if (HasCellValueData(pCell))
462cdf0e10cSrcweir pMat->PutDouble(GetCellValue(aAdr, pCell), 0);
463cdf0e10cSrcweir else
464cdf0e10cSrcweir {
465cdf0e10cSrcweir String aStr;
466cdf0e10cSrcweir GetCellString(aStr, pCell);
467cdf0e10cSrcweir pMat->PutString(aStr, 0);
468cdf0e10cSrcweir }
469cdf0e10cSrcweir }
470cdf0e10cSrcweir }
471cdf0e10cSrcweir break;
472cdf0e10cSrcweir case svDoubleRef:
473cdf0e10cSrcweir {
474cdf0e10cSrcweir SCCOL nCol1, nCol2;
475cdf0e10cSrcweir SCROW nRow1, nRow2;
476cdf0e10cSrcweir SCTAB nTab1, nTab2;
477cdf0e10cSrcweir const ScToken* p = sp ? static_cast<const ScToken*>(pStack[sp-1]) : NULL;
478cdf0e10cSrcweir PopDoubleRef(nCol1, nRow1, nTab1, nCol2, nRow2, nTab2);
479cdf0e10cSrcweir pMat = CreateMatrixFromDoubleRef( p, nCol1, nRow1, nTab1,
480cdf0e10cSrcweir nCol2, nRow2, nTab2);
481cdf0e10cSrcweir }
482cdf0e10cSrcweir break;
483cdf0e10cSrcweir case svMatrix:
484cdf0e10cSrcweir pMat = PopMatrix();
485cdf0e10cSrcweir break;
486cdf0e10cSrcweir case svError :
487cdf0e10cSrcweir case svMissing :
488cdf0e10cSrcweir case svDouble :
489cdf0e10cSrcweir {
490cdf0e10cSrcweir double fVal = GetDouble();
491cdf0e10cSrcweir pMat = GetNewMat( 1, 1);
492cdf0e10cSrcweir if ( pMat )
493cdf0e10cSrcweir {
494cdf0e10cSrcweir if ( nGlobalError )
495cdf0e10cSrcweir {
496cdf0e10cSrcweir fVal = CreateDoubleError( nGlobalError);
497cdf0e10cSrcweir nGlobalError = 0;
498cdf0e10cSrcweir }
499cdf0e10cSrcweir pMat->PutDouble( fVal, 0);
500cdf0e10cSrcweir }
501cdf0e10cSrcweir }
502cdf0e10cSrcweir break;
503cdf0e10cSrcweir case svString :
504cdf0e10cSrcweir {
505cdf0e10cSrcweir String aStr = GetString();
506cdf0e10cSrcweir pMat = GetNewMat( 1, 1);
507cdf0e10cSrcweir if ( pMat )
508cdf0e10cSrcweir {
509cdf0e10cSrcweir if ( nGlobalError )
510cdf0e10cSrcweir {
511cdf0e10cSrcweir double fVal = CreateDoubleError( nGlobalError);
512cdf0e10cSrcweir pMat->PutDouble( fVal, 0);
513cdf0e10cSrcweir nGlobalError = 0;
514cdf0e10cSrcweir }
515cdf0e10cSrcweir else
516cdf0e10cSrcweir pMat->PutString( aStr, 0);
517cdf0e10cSrcweir }
518cdf0e10cSrcweir }
519cdf0e10cSrcweir break;
520cdf0e10cSrcweir default:
521cdf0e10cSrcweir PopError();
522cdf0e10cSrcweir SetError( errIllegalArgument);
523cdf0e10cSrcweir break;
524cdf0e10cSrcweir }
525cdf0e10cSrcweir return pMat;
526cdf0e10cSrcweir }
527cdf0e10cSrcweir
ScMatValue()528cdf0e10cSrcweir void ScInterpreter::ScMatValue()
529cdf0e10cSrcweir {
530cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatValue" );
531cdf0e10cSrcweir if ( MustHaveParamCount( GetByte(), 3 ) )
532cdf0e10cSrcweir {
533cdf0e10cSrcweir // 0 to count-1
534cdf0e10cSrcweir SCSIZE nR = static_cast<SCSIZE>(::rtl::math::approxFloor(GetDouble()));
535cdf0e10cSrcweir SCSIZE nC = static_cast<SCSIZE>(::rtl::math::approxFloor(GetDouble()));
536cdf0e10cSrcweir switch (GetStackType())
537cdf0e10cSrcweir {
538cdf0e10cSrcweir case svSingleRef :
539cdf0e10cSrcweir {
540cdf0e10cSrcweir ScAddress aAdr;
541cdf0e10cSrcweir PopSingleRef( aAdr );
542cdf0e10cSrcweir ScBaseCell* pCell = GetCell( aAdr );
543cdf0e10cSrcweir if (pCell && pCell->GetCellType() == CELLTYPE_FORMULA)
544cdf0e10cSrcweir {
545cdf0e10cSrcweir sal_uInt16 nErrCode = ((ScFormulaCell*)pCell)->GetErrCode();
546cdf0e10cSrcweir if (nErrCode != 0)
547cdf0e10cSrcweir PushError( nErrCode);
548cdf0e10cSrcweir else
549cdf0e10cSrcweir {
550cdf0e10cSrcweir const ScMatrix* pMat = ((ScFormulaCell*)pCell)->GetMatrix();
551cdf0e10cSrcweir CalculateMatrixValue(pMat,nC,nR);
552cdf0e10cSrcweir }
553cdf0e10cSrcweir }
554cdf0e10cSrcweir else
555cdf0e10cSrcweir PushIllegalParameter();
556cdf0e10cSrcweir }
557cdf0e10cSrcweir break;
558cdf0e10cSrcweir case svDoubleRef :
559cdf0e10cSrcweir {
560cdf0e10cSrcweir SCCOL nCol1;
561cdf0e10cSrcweir SCROW nRow1;
562cdf0e10cSrcweir SCTAB nTab1;
563cdf0e10cSrcweir SCCOL nCol2;
564cdf0e10cSrcweir SCROW nRow2;
565cdf0e10cSrcweir SCTAB nTab2;
566cdf0e10cSrcweir PopDoubleRef(nCol1, nRow1, nTab1, nCol2, nRow2, nTab2);
567cdf0e10cSrcweir if (nCol2 - nCol1 >= static_cast<SCCOL>(nR) &&
568cdf0e10cSrcweir nRow2 - nRow1 >= static_cast<SCROW>(nC) &&
569cdf0e10cSrcweir nTab1 == nTab2)
570cdf0e10cSrcweir {
571cdf0e10cSrcweir ScAddress aAdr( sal::static_int_cast<SCCOL>( nCol1 + nR ),
572cdf0e10cSrcweir sal::static_int_cast<SCROW>( nRow1 + nC ), nTab1 );
573cdf0e10cSrcweir ScBaseCell* pCell = GetCell( aAdr );
574cdf0e10cSrcweir if (HasCellValueData(pCell))
575cdf0e10cSrcweir PushDouble(GetCellValue( aAdr, pCell ));
576cdf0e10cSrcweir else
577cdf0e10cSrcweir {
578cdf0e10cSrcweir String aStr;
579cdf0e10cSrcweir GetCellString(aStr, pCell);
580cdf0e10cSrcweir PushString(aStr);
581cdf0e10cSrcweir }
582cdf0e10cSrcweir }
583cdf0e10cSrcweir else
584cdf0e10cSrcweir PushNoValue();
585cdf0e10cSrcweir }
586cdf0e10cSrcweir break;
587cdf0e10cSrcweir case svMatrix:
588cdf0e10cSrcweir {
589cdf0e10cSrcweir ScMatrixRef pMat = PopMatrix();
590cdf0e10cSrcweir CalculateMatrixValue(pMat,nC,nR);
591cdf0e10cSrcweir }
592cdf0e10cSrcweir break;
593cdf0e10cSrcweir default:
594cdf0e10cSrcweir PopError();
595cdf0e10cSrcweir PushIllegalParameter();
596cdf0e10cSrcweir break;
597cdf0e10cSrcweir }
598cdf0e10cSrcweir }
599cdf0e10cSrcweir }
CalculateMatrixValue(const ScMatrix * pMat,SCSIZE nC,SCSIZE nR)600cdf0e10cSrcweir void ScInterpreter::CalculateMatrixValue(const ScMatrix* pMat,SCSIZE nC,SCSIZE nR)
601cdf0e10cSrcweir {
602cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CalculateMatrixValue" );
603cdf0e10cSrcweir if (pMat)
604cdf0e10cSrcweir {
605cdf0e10cSrcweir SCSIZE nCl, nRw;
606cdf0e10cSrcweir pMat->GetDimensions(nCl, nRw);
607cdf0e10cSrcweir if (nC < nCl && nR < nRw)
608cdf0e10cSrcweir {
609cdf0e10cSrcweir ScMatValType nMatValType;
610cdf0e10cSrcweir const ScMatrixValue* pMatVal = pMat->Get( nC, nR,nMatValType);
611cdf0e10cSrcweir if (ScMatrix::IsNonValueType( nMatValType))
612cdf0e10cSrcweir PushString( pMatVal->GetString() );
613cdf0e10cSrcweir else
614cdf0e10cSrcweir PushDouble(pMatVal->fVal);
615cdf0e10cSrcweir // also handles DoubleError
616cdf0e10cSrcweir }
617cdf0e10cSrcweir else
618cdf0e10cSrcweir PushNoValue();
619cdf0e10cSrcweir }
620cdf0e10cSrcweir else
621cdf0e10cSrcweir PushNoValue();
622cdf0e10cSrcweir }
623cdf0e10cSrcweir
ScEMat()624cdf0e10cSrcweir void ScInterpreter::ScEMat()
625cdf0e10cSrcweir {
626cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScEMat" );
627cdf0e10cSrcweir if ( MustHaveParamCount( GetByte(), 1 ) )
628cdf0e10cSrcweir {
629cdf0e10cSrcweir SCSIZE nDim = static_cast<SCSIZE>(::rtl::math::approxFloor(GetDouble()));
630cdf0e10cSrcweir if ( nDim * nDim > ScMatrix::GetElementsMax() || nDim == 0)
631cdf0e10cSrcweir PushIllegalArgument();
632cdf0e10cSrcweir else
633cdf0e10cSrcweir {
634cdf0e10cSrcweir ScMatrixRef pRMat = GetNewMat(nDim, nDim);
635cdf0e10cSrcweir if (pRMat)
636cdf0e10cSrcweir {
637cdf0e10cSrcweir MEMat(pRMat, nDim);
638cdf0e10cSrcweir PushMatrix(pRMat);
639cdf0e10cSrcweir }
640cdf0e10cSrcweir else
641cdf0e10cSrcweir PushIllegalArgument();
642cdf0e10cSrcweir }
643cdf0e10cSrcweir }
644cdf0e10cSrcweir }
645cdf0e10cSrcweir
MEMat(ScMatrix * mM,SCSIZE n)646cdf0e10cSrcweir void ScInterpreter::MEMat(ScMatrix* mM, SCSIZE n)
647cdf0e10cSrcweir {
648cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::MEMat" );
649cdf0e10cSrcweir mM->FillDouble(0.0, 0, 0, n-1, n-1);
650cdf0e10cSrcweir for (SCSIZE i = 0; i < n; i++)
651cdf0e10cSrcweir mM->PutDouble(1.0, i, i);
652cdf0e10cSrcweir }
653cdf0e10cSrcweir
654cdf0e10cSrcweir /* Matrix LUP decomposition according to the pseudocode of "Introduction to
655cdf0e10cSrcweir * Algorithms" by Cormen, Leiserson, Rivest, Stein.
656cdf0e10cSrcweir *
657cdf0e10cSrcweir * Added scaling for numeric stability.
658cdf0e10cSrcweir *
659cdf0e10cSrcweir * Given an n x n nonsingular matrix A, find a permutation matrix P, a unit
660cdf0e10cSrcweir * lower-triangular matrix L, and an upper-triangular matrix U such that PA=LU.
661cdf0e10cSrcweir * Compute L and U "in place" in the matrix A, the original content is
662cdf0e10cSrcweir * destroyed. Note that the diagonal elements of the U triangular matrix
663cdf0e10cSrcweir * replace the diagonal elements of the L-unit matrix (that are each ==1). The
664cdf0e10cSrcweir * permutation matrix P is an array, where P[i]=j means that the i-th row of P
665cdf0e10cSrcweir * contains a 1 in column j. Additionally keep track of the number of
666cdf0e10cSrcweir * permutations (row exchanges).
667cdf0e10cSrcweir *
668cdf0e10cSrcweir * Returns 0 if a singular matrix is encountered, else +1 if an even number of
66986e1cf34SPedro Giffuni * permutations occurred, or -1 if odd, which is the sign of the determinant.
670cdf0e10cSrcweir * This may be used to calculate the determinant by multiplying the sign with
671cdf0e10cSrcweir * the product of the diagonal elements of the LU matrix.
672cdf0e10cSrcweir */
lcl_LUP_decompose(ScMatrix * mA,const SCSIZE n,::std::vector<SCSIZE> & P)673cdf0e10cSrcweir static int lcl_LUP_decompose( ScMatrix* mA, const SCSIZE n,
674cdf0e10cSrcweir ::std::vector< SCSIZE> & P )
675cdf0e10cSrcweir {
676cdf0e10cSrcweir int nSign = 1;
677cdf0e10cSrcweir // Find scale of each row.
678cdf0e10cSrcweir ::std::vector< double> aScale(n);
679cdf0e10cSrcweir for (SCSIZE i=0; i < n; ++i)
680cdf0e10cSrcweir {
681cdf0e10cSrcweir double fMax = 0.0;
682cdf0e10cSrcweir for (SCSIZE j=0; j < n; ++j)
683cdf0e10cSrcweir {
684cdf0e10cSrcweir double fTmp = fabs( mA->GetDouble( j, i));
685cdf0e10cSrcweir if (fMax < fTmp)
686cdf0e10cSrcweir fMax = fTmp;
687cdf0e10cSrcweir }
688cdf0e10cSrcweir if (fMax == 0.0)
689cdf0e10cSrcweir return 0; // singular matrix
690cdf0e10cSrcweir aScale[i] = 1.0 / fMax;
691cdf0e10cSrcweir }
692cdf0e10cSrcweir // Represent identity permutation, P[i]=i
693cdf0e10cSrcweir for (SCSIZE i=0; i < n; ++i)
694cdf0e10cSrcweir P[i] = i;
695cdf0e10cSrcweir // "Recursion" on the diagonale.
696cdf0e10cSrcweir SCSIZE l = n - 1;
697cdf0e10cSrcweir for (SCSIZE k=0; k < l; ++k)
698cdf0e10cSrcweir {
699cdf0e10cSrcweir // Implicit pivoting. With the scale found for a row, compare values of
700cdf0e10cSrcweir // a column and pick largest.
701cdf0e10cSrcweir double fMax = 0.0;
702cdf0e10cSrcweir double fScale = aScale[k];
703cdf0e10cSrcweir SCSIZE kp = k;
704cdf0e10cSrcweir for (SCSIZE i = k; i < n; ++i)
705cdf0e10cSrcweir {
706cdf0e10cSrcweir double fTmp = fScale * fabs( mA->GetDouble( k, i));
707cdf0e10cSrcweir if (fMax < fTmp)
708cdf0e10cSrcweir {
709cdf0e10cSrcweir fMax = fTmp;
710cdf0e10cSrcweir kp = i;
711cdf0e10cSrcweir }
712cdf0e10cSrcweir }
713cdf0e10cSrcweir if (fMax == 0.0)
714cdf0e10cSrcweir return 0; // singular matrix
715cdf0e10cSrcweir // Swap rows. The pivot element will be at mA[k,kp] (row,col notation)
716cdf0e10cSrcweir if (k != kp)
717cdf0e10cSrcweir {
718cdf0e10cSrcweir // permutations
719cdf0e10cSrcweir SCSIZE nTmp = P[k];
720cdf0e10cSrcweir P[k] = P[kp];
721cdf0e10cSrcweir P[kp] = nTmp;
722cdf0e10cSrcweir nSign = -nSign;
723cdf0e10cSrcweir // scales
724cdf0e10cSrcweir double fTmp = aScale[k];
725cdf0e10cSrcweir aScale[k] = aScale[kp];
726cdf0e10cSrcweir aScale[kp] = fTmp;
727cdf0e10cSrcweir // elements
728cdf0e10cSrcweir for (SCSIZE i=0; i < n; ++i)
729cdf0e10cSrcweir {
730cdf0e10cSrcweir double fMatTmp = mA->GetDouble( i, k);
731cdf0e10cSrcweir mA->PutDouble( mA->GetDouble( i, kp), i, k);
732cdf0e10cSrcweir mA->PutDouble( fMatTmp, i, kp);
733cdf0e10cSrcweir }
734cdf0e10cSrcweir }
735cdf0e10cSrcweir // Compute Schur complement.
736cdf0e10cSrcweir for (SCSIZE i = k+1; i < n; ++i)
737cdf0e10cSrcweir {
738cdf0e10cSrcweir double fTmp = mA->GetDouble( k, i) / mA->GetDouble( k, k);
739cdf0e10cSrcweir mA->PutDouble( fTmp, k, i);
740cdf0e10cSrcweir for (SCSIZE j = k+1; j < n; ++j)
741cdf0e10cSrcweir mA->PutDouble( mA->GetDouble( j, i) - fTmp * mA->GetDouble( j,
742cdf0e10cSrcweir k), j, i);
743cdf0e10cSrcweir }
744cdf0e10cSrcweir }
745cdf0e10cSrcweir #if OSL_DEBUG_LEVEL > 1
746cdf0e10cSrcweir fprintf( stderr, "\n%s\n", "lcl_LUP_decompose(): LU");
747cdf0e10cSrcweir for (SCSIZE i=0; i < n; ++i)
748cdf0e10cSrcweir {
749cdf0e10cSrcweir for (SCSIZE j=0; j < n; ++j)
750cdf0e10cSrcweir fprintf( stderr, "%8.2g ", mA->GetDouble( j, i));
751cdf0e10cSrcweir fprintf( stderr, "\n%s\n", "");
752cdf0e10cSrcweir }
753cdf0e10cSrcweir fprintf( stderr, "\n%s\n", "lcl_LUP_decompose(): P");
754cdf0e10cSrcweir for (SCSIZE j=0; j < n; ++j)
755cdf0e10cSrcweir fprintf( stderr, "%5u ", (unsigned)P[j]);
756cdf0e10cSrcweir fprintf( stderr, "\n%s\n", "");
757cdf0e10cSrcweir #endif
758cdf0e10cSrcweir
759cdf0e10cSrcweir bool bSingular=false;
760cdf0e10cSrcweir for (SCSIZE i=0; i < n && !bSingular; i++)
761cdf0e10cSrcweir bSingular = (mA->GetDouble(i,i) == 0.0);
762cdf0e10cSrcweir if (bSingular)
763cdf0e10cSrcweir nSign = 0;
764cdf0e10cSrcweir
765cdf0e10cSrcweir return nSign;
766cdf0e10cSrcweir }
767cdf0e10cSrcweir
768cdf0e10cSrcweir
769cdf0e10cSrcweir /* Solve a LUP decomposed equation Ax=b. LU is a combined matrix of L and U
770cdf0e10cSrcweir * triangulars and P the permutation vector as obtained from
771cdf0e10cSrcweir * lcl_LUP_decompose(). B is the right-hand side input vector, X is used to
772cdf0e10cSrcweir * return the solution vector.
773cdf0e10cSrcweir */
lcl_LUP_solve(const ScMatrix * mLU,const SCSIZE n,const::std::vector<SCSIZE> & P,const::std::vector<double> & B,::std::vector<double> & X)774cdf0e10cSrcweir static void lcl_LUP_solve( const ScMatrix* mLU, const SCSIZE n,
775cdf0e10cSrcweir const ::std::vector< SCSIZE> & P, const ::std::vector< double> & B,
776cdf0e10cSrcweir ::std::vector< double> & X )
777cdf0e10cSrcweir {
778cdf0e10cSrcweir SCSIZE nFirst = SCSIZE_MAX;
779cdf0e10cSrcweir // Ax=b => PAx=Pb, with decomposition LUx=Pb.
780cdf0e10cSrcweir // Define y=Ux and solve for y in Ly=Pb using forward substitution.
781cdf0e10cSrcweir for (SCSIZE i=0; i < n; ++i)
782cdf0e10cSrcweir {
783cdf0e10cSrcweir double fSum = B[P[i]];
784cdf0e10cSrcweir // Matrix inversion comes with a lot of zeros in the B vectors, we
785cdf0e10cSrcweir // don't have to do all the computing with results multiplied by zero.
786cdf0e10cSrcweir // Until then, simply lookout for the position of the first nonzero
787cdf0e10cSrcweir // value.
788cdf0e10cSrcweir if (nFirst != SCSIZE_MAX)
789cdf0e10cSrcweir {
790cdf0e10cSrcweir for (SCSIZE j = nFirst; j < i; ++j)
791cdf0e10cSrcweir fSum -= mLU->GetDouble( j, i) * X[j]; // X[j] === y[j]
792cdf0e10cSrcweir }
793cdf0e10cSrcweir else if (fSum)
794cdf0e10cSrcweir nFirst = i;
795cdf0e10cSrcweir X[i] = fSum; // X[i] === y[i]
796cdf0e10cSrcweir }
797cdf0e10cSrcweir // Solve for x in Ux=y using back substitution.
798cdf0e10cSrcweir for (SCSIZE i = n; i--; )
799cdf0e10cSrcweir {
800cdf0e10cSrcweir double fSum = X[i]; // X[i] === y[i]
801cdf0e10cSrcweir for (SCSIZE j = i+1; j < n; ++j)
802cdf0e10cSrcweir fSum -= mLU->GetDouble( j, i) * X[j]; // X[j] === x[j]
803cdf0e10cSrcweir X[i] = fSum / mLU->GetDouble( i, i); // X[i] === x[i]
804cdf0e10cSrcweir }
805cdf0e10cSrcweir #if OSL_DEBUG_LEVEL >1
806cdf0e10cSrcweir fprintf( stderr, "\n%s\n", "lcl_LUP_solve():");
807cdf0e10cSrcweir for (SCSIZE i=0; i < n; ++i)
808cdf0e10cSrcweir fprintf( stderr, "%8.2g ", X[i]);
809cdf0e10cSrcweir fprintf( stderr, "%s\n", "");
810cdf0e10cSrcweir #endif
811cdf0e10cSrcweir }
812cdf0e10cSrcweir
813cdf0e10cSrcweir
ScMatDet()814cdf0e10cSrcweir void ScInterpreter::ScMatDet()
815cdf0e10cSrcweir {
816cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatDet" );
817cdf0e10cSrcweir if ( MustHaveParamCount( GetByte(), 1 ) )
818cdf0e10cSrcweir {
819cdf0e10cSrcweir ScMatrixRef pMat = GetMatrix();
820cdf0e10cSrcweir if (!pMat)
821cdf0e10cSrcweir {
822cdf0e10cSrcweir PushIllegalParameter();
823cdf0e10cSrcweir return;
824cdf0e10cSrcweir }
825cdf0e10cSrcweir if ( !pMat->IsNumeric() )
826cdf0e10cSrcweir {
827cdf0e10cSrcweir PushNoValue();
828cdf0e10cSrcweir return;
829cdf0e10cSrcweir }
830cdf0e10cSrcweir SCSIZE nC, nR;
831cdf0e10cSrcweir pMat->GetDimensions(nC, nR);
832cdf0e10cSrcweir if ( nC != nR || nC == 0 || (sal_uLong) nC * nC > ScMatrix::GetElementsMax() )
833cdf0e10cSrcweir PushIllegalArgument();
834cdf0e10cSrcweir else
835cdf0e10cSrcweir {
836cdf0e10cSrcweir // LUP decomposition is done inplace, use copy.
837cdf0e10cSrcweir ScMatrixRef xLU = pMat->Clone();
838cdf0e10cSrcweir if (!xLU)
839cdf0e10cSrcweir PushError( errCodeOverflow);
840cdf0e10cSrcweir else
841cdf0e10cSrcweir {
842cdf0e10cSrcweir ::std::vector< SCSIZE> P(nR);
843cdf0e10cSrcweir int nDetSign = lcl_LUP_decompose( xLU, nR, P);
844cdf0e10cSrcweir if (!nDetSign)
845cdf0e10cSrcweir PushInt(0); // singular matrix
846cdf0e10cSrcweir else
847cdf0e10cSrcweir {
848cdf0e10cSrcweir // In an LU matrix the determinant is simply the product of
849cdf0e10cSrcweir // all diagonal elements.
850cdf0e10cSrcweir double fDet = nDetSign;
851cdf0e10cSrcweir ScMatrix* pLU = xLU;
852cdf0e10cSrcweir for (SCSIZE i=0; i < nR; ++i)
853cdf0e10cSrcweir fDet *= pLU->GetDouble( i, i);
854cdf0e10cSrcweir PushDouble( fDet);
855cdf0e10cSrcweir }
856cdf0e10cSrcweir }
857cdf0e10cSrcweir }
858cdf0e10cSrcweir }
859cdf0e10cSrcweir }
860cdf0e10cSrcweir
ScMatInv()861cdf0e10cSrcweir void ScInterpreter::ScMatInv()
862cdf0e10cSrcweir {
863cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatInv" );
864cdf0e10cSrcweir if ( MustHaveParamCount( GetByte(), 1 ) )
865cdf0e10cSrcweir {
866cdf0e10cSrcweir ScMatrixRef pMat = GetMatrix();
867cdf0e10cSrcweir if (!pMat)
868cdf0e10cSrcweir {
869cdf0e10cSrcweir PushIllegalParameter();
870cdf0e10cSrcweir return;
871cdf0e10cSrcweir }
872cdf0e10cSrcweir if ( !pMat->IsNumeric() )
873cdf0e10cSrcweir {
874cdf0e10cSrcweir PushNoValue();
875cdf0e10cSrcweir return;
876cdf0e10cSrcweir }
877cdf0e10cSrcweir SCSIZE nC, nR;
878cdf0e10cSrcweir pMat->GetDimensions(nC, nR);
879cdf0e10cSrcweir if ( nC != nR || nC == 0 || (sal_uLong) nC * nC > ScMatrix::GetElementsMax() )
880cdf0e10cSrcweir PushIllegalArgument();
881cdf0e10cSrcweir else
882cdf0e10cSrcweir {
883cdf0e10cSrcweir // LUP decomposition is done inplace, use copy.
884cdf0e10cSrcweir ScMatrixRef xLU = pMat->Clone();
885cdf0e10cSrcweir // The result matrix.
886cdf0e10cSrcweir ScMatrixRef xY = GetNewMat( nR, nR);
887cdf0e10cSrcweir if (!xLU || !xY)
888cdf0e10cSrcweir PushError( errCodeOverflow);
889cdf0e10cSrcweir else
890cdf0e10cSrcweir {
891cdf0e10cSrcweir ::std::vector< SCSIZE> P(nR);
892cdf0e10cSrcweir int nDetSign = lcl_LUP_decompose( xLU, nR, P);
893cdf0e10cSrcweir if (!nDetSign)
894cdf0e10cSrcweir PushIllegalArgument();
895cdf0e10cSrcweir else
896cdf0e10cSrcweir {
897cdf0e10cSrcweir // Solve equation for each column.
898cdf0e10cSrcweir ScMatrix* pY = xY;
899cdf0e10cSrcweir ::std::vector< double> B(nR);
900cdf0e10cSrcweir ::std::vector< double> X(nR);
901cdf0e10cSrcweir for (SCSIZE j=0; j < nR; ++j)
902cdf0e10cSrcweir {
903cdf0e10cSrcweir for (SCSIZE i=0; i < nR; ++i)
904cdf0e10cSrcweir B[i] = 0.0;
905cdf0e10cSrcweir B[j] = 1.0;
906cdf0e10cSrcweir lcl_LUP_solve( xLU, nR, P, B, X);
907cdf0e10cSrcweir for (SCSIZE i=0; i < nR; ++i)
908cdf0e10cSrcweir pY->PutDouble( X[i], j, i);
909cdf0e10cSrcweir }
910cdf0e10cSrcweir #if 0
911cdf0e10cSrcweir /* Possible checks for ill-condition:
912cdf0e10cSrcweir * 1. Scale matrix, invert scaled matrix. If there are
913cdf0e10cSrcweir * elements of the inverted matrix that are several
914cdf0e10cSrcweir * orders of magnitude greater than 1 =>
915cdf0e10cSrcweir * ill-conditioned.
916cdf0e10cSrcweir * Just how much is "several orders"?
917cdf0e10cSrcweir * 2. Invert the inverted matrix and assess whether the
918cdf0e10cSrcweir * result is sufficiently close to the original matrix.
919cdf0e10cSrcweir * If not => ill-conditioned.
920cdf0e10cSrcweir * Just what is sufficient?
921cdf0e10cSrcweir * 3. Multiplying the inverse by the original matrix should
922cdf0e10cSrcweir * produce a result sufficiently close to the identity
923cdf0e10cSrcweir * matrix.
924cdf0e10cSrcweir * Just what is sufficient?
925cdf0e10cSrcweir *
926cdf0e10cSrcweir * The following is #3.
927cdf0e10cSrcweir */
928cdf0e10cSrcweir ScMatrixRef xR = GetNewMat( nR, nR);
929cdf0e10cSrcweir if (xR)
930cdf0e10cSrcweir {
931cdf0e10cSrcweir ScMatrix* pR = xR;
932cdf0e10cSrcweir lcl_MFastMult( pMat, pY, pR, nR, nR, nR);
933cdf0e10cSrcweir #if OSL_DEBUG_LEVEL > 1
934cdf0e10cSrcweir fprintf( stderr, "\n%s\n", "ScMatInv(): mult-identity");
935cdf0e10cSrcweir #endif
936cdf0e10cSrcweir for (SCSIZE i=0; i < nR; ++i)
937cdf0e10cSrcweir {
938cdf0e10cSrcweir for (SCSIZE j=0; j < nR; ++j)
939cdf0e10cSrcweir {
940cdf0e10cSrcweir double fTmp = pR->GetDouble( j, i);
941cdf0e10cSrcweir #if OSL_DEBUG_LEVEL > 1
942cdf0e10cSrcweir fprintf( stderr, "%8.2g ", fTmp);
943cdf0e10cSrcweir #endif
944cdf0e10cSrcweir if (fabs( fTmp - (i == j)) > fInvEpsilon)
945cdf0e10cSrcweir SetError( errIllegalArgument);
946cdf0e10cSrcweir }
947cdf0e10cSrcweir #if OSL_DEBUG_LEVEL > 1
948cdf0e10cSrcweir fprintf( stderr, "\n%s\n", "");
949cdf0e10cSrcweir #endif
950cdf0e10cSrcweir }
951cdf0e10cSrcweir }
952cdf0e10cSrcweir #endif
953cdf0e10cSrcweir if (nGlobalError)
954cdf0e10cSrcweir PushError( nGlobalError);
955cdf0e10cSrcweir else
956cdf0e10cSrcweir PushMatrix( pY);
957cdf0e10cSrcweir }
958cdf0e10cSrcweir }
959cdf0e10cSrcweir }
960cdf0e10cSrcweir }
961cdf0e10cSrcweir }
962cdf0e10cSrcweir
ScMatMult()963cdf0e10cSrcweir void ScInterpreter::ScMatMult()
964cdf0e10cSrcweir {
965cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatMult" );
966cdf0e10cSrcweir if ( MustHaveParamCount( GetByte(), 2 ) )
967cdf0e10cSrcweir {
968cdf0e10cSrcweir ScMatrixRef pMat2 = GetMatrix();
969cdf0e10cSrcweir ScMatrixRef pMat1 = GetMatrix();
970cdf0e10cSrcweir ScMatrixRef pRMat;
971cdf0e10cSrcweir if (pMat1 && pMat2)
972cdf0e10cSrcweir {
973cdf0e10cSrcweir if ( pMat1->IsNumeric() && pMat2->IsNumeric() )
974cdf0e10cSrcweir {
975cdf0e10cSrcweir SCSIZE nC1, nC2;
976cdf0e10cSrcweir SCSIZE nR1, nR2;
977cdf0e10cSrcweir pMat1->GetDimensions(nC1, nR1);
978cdf0e10cSrcweir pMat2->GetDimensions(nC2, nR2);
979cdf0e10cSrcweir if (nC1 != nR2)
980cdf0e10cSrcweir PushIllegalArgument();
981cdf0e10cSrcweir else
982cdf0e10cSrcweir {
983cdf0e10cSrcweir pRMat = GetNewMat(nC2, nR1);
984cdf0e10cSrcweir if (pRMat)
985cdf0e10cSrcweir {
986cdf0e10cSrcweir double sum;
987cdf0e10cSrcweir for (SCSIZE i = 0; i < nR1; i++)
988cdf0e10cSrcweir {
989cdf0e10cSrcweir for (SCSIZE j = 0; j < nC2; j++)
990cdf0e10cSrcweir {
991cdf0e10cSrcweir sum = 0.0;
992cdf0e10cSrcweir for (SCSIZE k = 0; k < nC1; k++)
993cdf0e10cSrcweir {
994cdf0e10cSrcweir sum += pMat1->GetDouble(k,i)*pMat2->GetDouble(j,k);
995cdf0e10cSrcweir }
996cdf0e10cSrcweir pRMat->PutDouble(sum, j, i);
997cdf0e10cSrcweir }
998cdf0e10cSrcweir }
999cdf0e10cSrcweir PushMatrix(pRMat);
1000cdf0e10cSrcweir }
1001cdf0e10cSrcweir else
1002cdf0e10cSrcweir PushIllegalArgument();
1003cdf0e10cSrcweir }
1004cdf0e10cSrcweir }
1005cdf0e10cSrcweir else
1006cdf0e10cSrcweir PushNoValue();
1007cdf0e10cSrcweir }
1008cdf0e10cSrcweir else
1009cdf0e10cSrcweir PushIllegalParameter();
1010cdf0e10cSrcweir }
1011cdf0e10cSrcweir }
1012cdf0e10cSrcweir
ScMatTrans()1013cdf0e10cSrcweir void ScInterpreter::ScMatTrans()
1014cdf0e10cSrcweir {
1015cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatTrans" );
1016cdf0e10cSrcweir if ( MustHaveParamCount( GetByte(), 1 ) )
1017cdf0e10cSrcweir {
1018cdf0e10cSrcweir ScMatrixRef pMat = GetMatrix();
1019cdf0e10cSrcweir ScMatrixRef pRMat;
1020cdf0e10cSrcweir if (pMat)
1021cdf0e10cSrcweir {
1022cdf0e10cSrcweir SCSIZE nC, nR;
1023cdf0e10cSrcweir pMat->GetDimensions(nC, nR);
1024cdf0e10cSrcweir pRMat = GetNewMat(nR, nC);
1025cdf0e10cSrcweir if ( pRMat )
1026cdf0e10cSrcweir {
1027cdf0e10cSrcweir pMat->MatTrans(*pRMat);
1028cdf0e10cSrcweir PushMatrix(pRMat);
1029cdf0e10cSrcweir }
1030cdf0e10cSrcweir else
1031cdf0e10cSrcweir PushIllegalArgument();
1032cdf0e10cSrcweir }
1033cdf0e10cSrcweir else
1034cdf0e10cSrcweir PushIllegalParameter();
1035cdf0e10cSrcweir }
1036cdf0e10cSrcweir }
1037cdf0e10cSrcweir
1038cdf0e10cSrcweir
1039cdf0e10cSrcweir /** Minimum extent of one result matrix dimension.
1040cdf0e10cSrcweir For a row or column vector to be replicated the larger matrix dimension is
1041cdf0e10cSrcweir returned, else the smaller dimension.
1042cdf0e10cSrcweir */
lcl_GetMinExtent(SCSIZE n1,SCSIZE n2)1043cdf0e10cSrcweir inline SCSIZE lcl_GetMinExtent( SCSIZE n1, SCSIZE n2 )
1044cdf0e10cSrcweir {
1045cdf0e10cSrcweir if (n1 == 1)
1046cdf0e10cSrcweir return n2;
1047cdf0e10cSrcweir else if (n2 == 1)
1048cdf0e10cSrcweir return n1;
1049cdf0e10cSrcweir else if (n1 < n2)
1050cdf0e10cSrcweir return n1;
1051cdf0e10cSrcweir else
1052cdf0e10cSrcweir return n2;
1053cdf0e10cSrcweir }
1054cdf0e10cSrcweir
1055cdf0e10cSrcweir template<class _Function>
lcl_MatrixCalculation(const _Function & _pOperation,ScMatrix * pMat1,ScMatrix * pMat2,ScInterpreter * _pIterpreter)1056cdf0e10cSrcweir ScMatrixRef lcl_MatrixCalculation(const _Function& _pOperation,ScMatrix* pMat1, ScMatrix* pMat2,ScInterpreter* _pIterpreter)
1057cdf0e10cSrcweir {
1058cdf0e10cSrcweir SCSIZE nC1, nC2, nMinC;
1059cdf0e10cSrcweir SCSIZE nR1, nR2, nMinR;
1060cdf0e10cSrcweir SCSIZE i, j;
1061cdf0e10cSrcweir pMat1->GetDimensions(nC1, nR1);
1062cdf0e10cSrcweir pMat2->GetDimensions(nC2, nR2);
1063cdf0e10cSrcweir nMinC = lcl_GetMinExtent( nC1, nC2);
1064cdf0e10cSrcweir nMinR = lcl_GetMinExtent( nR1, nR2);
1065cdf0e10cSrcweir ScMatrixRef xResMat = _pIterpreter->GetNewMat(nMinC, nMinR);
1066cdf0e10cSrcweir if (xResMat)
1067cdf0e10cSrcweir {
1068cdf0e10cSrcweir ScMatrix* pResMat = xResMat;
1069cdf0e10cSrcweir for (i = 0; i < nMinC; i++)
1070cdf0e10cSrcweir {
1071cdf0e10cSrcweir for (j = 0; j < nMinR; j++)
1072cdf0e10cSrcweir {
1073cdf0e10cSrcweir if (pMat1->IsValueOrEmpty(i,j) && pMat2->IsValueOrEmpty(i,j))
1074cdf0e10cSrcweir {
1075cdf0e10cSrcweir double d = _pOperation(pMat1->GetDouble(i,j),pMat2->GetDouble(i,j));
1076cdf0e10cSrcweir pResMat->PutDouble( d, i, j);
1077cdf0e10cSrcweir }
1078cdf0e10cSrcweir else
1079cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i, j);
1080cdf0e10cSrcweir }
1081cdf0e10cSrcweir }
1082cdf0e10cSrcweir }
1083cdf0e10cSrcweir return xResMat;
1084cdf0e10cSrcweir }
1085cdf0e10cSrcweir
MatConcat(ScMatrix * pMat1,ScMatrix * pMat2)1086cdf0e10cSrcweir ScMatrixRef ScInterpreter::MatConcat(ScMatrix* pMat1, ScMatrix* pMat2)
1087cdf0e10cSrcweir {
1088cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::MatConcat" );
1089cdf0e10cSrcweir SCSIZE nC1, nC2, nMinC;
1090cdf0e10cSrcweir SCSIZE nR1, nR2, nMinR;
1091cdf0e10cSrcweir SCSIZE i, j;
1092cdf0e10cSrcweir pMat1->GetDimensions(nC1, nR1);
1093cdf0e10cSrcweir pMat2->GetDimensions(nC2, nR2);
1094cdf0e10cSrcweir nMinC = lcl_GetMinExtent( nC1, nC2);
1095cdf0e10cSrcweir nMinR = lcl_GetMinExtent( nR1, nR2);
1096cdf0e10cSrcweir ScMatrixRef xResMat = GetNewMat(nMinC, nMinR);
1097cdf0e10cSrcweir if (xResMat)
1098cdf0e10cSrcweir {
1099cdf0e10cSrcweir ScMatrix* pResMat = xResMat;
1100cdf0e10cSrcweir for (i = 0; i < nMinC; i++)
1101cdf0e10cSrcweir {
1102cdf0e10cSrcweir for (j = 0; j < nMinR; j++)
1103cdf0e10cSrcweir {
1104cdf0e10cSrcweir sal_uInt16 nErr = pMat1->GetErrorIfNotString( i, j);
1105cdf0e10cSrcweir if (!nErr)
1106cdf0e10cSrcweir nErr = pMat2->GetErrorIfNotString( i, j);
1107cdf0e10cSrcweir if (nErr)
1108cdf0e10cSrcweir pResMat->PutError( nErr, i, j);
1109cdf0e10cSrcweir else
1110cdf0e10cSrcweir {
1111cdf0e10cSrcweir String aTmp( pMat1->GetString( *pFormatter, i, j));
1112cdf0e10cSrcweir aTmp += pMat2->GetString( *pFormatter, i, j);
1113cdf0e10cSrcweir pResMat->PutString( aTmp, i, j);
1114cdf0e10cSrcweir }
1115cdf0e10cSrcweir }
1116cdf0e10cSrcweir }
1117cdf0e10cSrcweir }
1118cdf0e10cSrcweir return xResMat;
1119cdf0e10cSrcweir }
1120cdf0e10cSrcweir
1121cdf0e10cSrcweir
1122cdf0e10cSrcweir // fuer DATE, TIME, DATETIME
lcl_GetDiffDateTimeFmtType(short & nFuncFmt,short nFmt1,short nFmt2)1123cdf0e10cSrcweir void lcl_GetDiffDateTimeFmtType( short& nFuncFmt, short nFmt1, short nFmt2 )
1124cdf0e10cSrcweir {
1125cdf0e10cSrcweir if ( nFmt1 != NUMBERFORMAT_UNDEFINED || nFmt2 != NUMBERFORMAT_UNDEFINED )
1126cdf0e10cSrcweir {
1127cdf0e10cSrcweir if ( nFmt1 == nFmt2 )
1128cdf0e10cSrcweir {
1129cdf0e10cSrcweir if ( nFmt1 == NUMBERFORMAT_TIME || nFmt1 == NUMBERFORMAT_DATETIME )
1130cdf0e10cSrcweir nFuncFmt = NUMBERFORMAT_TIME; // Zeiten ergeben Zeit
1131cdf0e10cSrcweir // else: nichts besonderes, Zahl (Datum - Datum := Tage)
1132cdf0e10cSrcweir }
1133cdf0e10cSrcweir else if ( nFmt1 == NUMBERFORMAT_UNDEFINED )
1134cdf0e10cSrcweir nFuncFmt = nFmt2; // z.B. Datum + Tage := Datum
1135cdf0e10cSrcweir else if ( nFmt2 == NUMBERFORMAT_UNDEFINED )
1136cdf0e10cSrcweir nFuncFmt = nFmt1;
1137cdf0e10cSrcweir else
1138cdf0e10cSrcweir {
1139cdf0e10cSrcweir if ( nFmt1 == NUMBERFORMAT_DATE || nFmt2 == NUMBERFORMAT_DATE ||
1140cdf0e10cSrcweir nFmt1 == NUMBERFORMAT_DATETIME || nFmt2 == NUMBERFORMAT_DATETIME )
1141cdf0e10cSrcweir {
1142cdf0e10cSrcweir if ( nFmt1 == NUMBERFORMAT_TIME || nFmt2 == NUMBERFORMAT_TIME )
1143cdf0e10cSrcweir nFuncFmt = NUMBERFORMAT_DATETIME; // Datum + Zeit
1144cdf0e10cSrcweir }
1145cdf0e10cSrcweir }
1146cdf0e10cSrcweir }
1147cdf0e10cSrcweir }
1148cdf0e10cSrcweir
1149cdf0e10cSrcweir
ScAdd()1150cdf0e10cSrcweir void ScInterpreter::ScAdd()
1151cdf0e10cSrcweir {
1152cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScAdd" );
1153cdf0e10cSrcweir CalculateAddSub(sal_False);
1154cdf0e10cSrcweir }
CalculateAddSub(sal_Bool _bSub)1155cdf0e10cSrcweir void ScInterpreter::CalculateAddSub(sal_Bool _bSub)
1156cdf0e10cSrcweir {
1157cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CalculateAddSub" );
1158cdf0e10cSrcweir ScMatrixRef pMat1 = NULL;
1159cdf0e10cSrcweir ScMatrixRef pMat2 = NULL;
1160cdf0e10cSrcweir double fVal1 = 0.0, fVal2 = 0.0;
1161cdf0e10cSrcweir short nFmt1, nFmt2;
1162cdf0e10cSrcweir nFmt1 = nFmt2 = NUMBERFORMAT_UNDEFINED;
1163cdf0e10cSrcweir short nFmtCurrencyType = nCurFmtType;
1164cdf0e10cSrcweir sal_uLong nFmtCurrencyIndex = nCurFmtIndex;
1165cdf0e10cSrcweir short nFmtPercentType = nCurFmtType;
1166cdf0e10cSrcweir if ( GetStackType() == svMatrix )
1167cdf0e10cSrcweir pMat2 = GetMatrix();
1168cdf0e10cSrcweir else
1169cdf0e10cSrcweir {
1170cdf0e10cSrcweir fVal2 = GetDouble();
1171cdf0e10cSrcweir switch ( nCurFmtType )
1172cdf0e10cSrcweir {
1173cdf0e10cSrcweir case NUMBERFORMAT_DATE :
1174cdf0e10cSrcweir case NUMBERFORMAT_TIME :
1175cdf0e10cSrcweir case NUMBERFORMAT_DATETIME :
1176cdf0e10cSrcweir nFmt2 = nCurFmtType;
1177cdf0e10cSrcweir break;
1178cdf0e10cSrcweir case NUMBERFORMAT_CURRENCY :
1179cdf0e10cSrcweir nFmtCurrencyType = nCurFmtType;
1180cdf0e10cSrcweir nFmtCurrencyIndex = nCurFmtIndex;
1181cdf0e10cSrcweir break;
1182cdf0e10cSrcweir case NUMBERFORMAT_PERCENT :
1183cdf0e10cSrcweir nFmtPercentType = NUMBERFORMAT_PERCENT;
1184cdf0e10cSrcweir break;
1185cdf0e10cSrcweir }
1186cdf0e10cSrcweir }
1187cdf0e10cSrcweir if ( GetStackType() == svMatrix )
1188cdf0e10cSrcweir pMat1 = GetMatrix();
1189cdf0e10cSrcweir else
1190cdf0e10cSrcweir {
1191cdf0e10cSrcweir fVal1 = GetDouble();
1192cdf0e10cSrcweir switch ( nCurFmtType )
1193cdf0e10cSrcweir {
1194cdf0e10cSrcweir case NUMBERFORMAT_DATE :
1195cdf0e10cSrcweir case NUMBERFORMAT_TIME :
1196cdf0e10cSrcweir case NUMBERFORMAT_DATETIME :
1197cdf0e10cSrcweir nFmt1 = nCurFmtType;
1198cdf0e10cSrcweir break;
1199cdf0e10cSrcweir case NUMBERFORMAT_CURRENCY :
1200cdf0e10cSrcweir nFmtCurrencyType = nCurFmtType;
1201cdf0e10cSrcweir nFmtCurrencyIndex = nCurFmtIndex;
1202cdf0e10cSrcweir break;
1203cdf0e10cSrcweir case NUMBERFORMAT_PERCENT :
1204cdf0e10cSrcweir nFmtPercentType = NUMBERFORMAT_PERCENT;
1205cdf0e10cSrcweir break;
1206cdf0e10cSrcweir }
1207cdf0e10cSrcweir }
1208cdf0e10cSrcweir if (pMat1 && pMat2)
1209cdf0e10cSrcweir {
1210cdf0e10cSrcweir ScMatrixRef pResMat;
1211cdf0e10cSrcweir if ( _bSub )
1212cdf0e10cSrcweir {
1213cdf0e10cSrcweir MatrixSub aSub;
1214cdf0e10cSrcweir pResMat = lcl_MatrixCalculation(aSub ,pMat1, pMat2,this);
1215cdf0e10cSrcweir }
1216cdf0e10cSrcweir else
1217cdf0e10cSrcweir {
1218cdf0e10cSrcweir MatrixAdd aAdd;
1219cdf0e10cSrcweir pResMat = lcl_MatrixCalculation(aAdd ,pMat1, pMat2,this);
1220cdf0e10cSrcweir }
1221cdf0e10cSrcweir
1222cdf0e10cSrcweir if (!pResMat)
1223cdf0e10cSrcweir PushNoValue();
1224cdf0e10cSrcweir else
1225cdf0e10cSrcweir PushMatrix(pResMat);
1226cdf0e10cSrcweir }
1227cdf0e10cSrcweir else if (pMat1 || pMat2)
1228cdf0e10cSrcweir {
1229cdf0e10cSrcweir double fVal;
1230cdf0e10cSrcweir sal_Bool bFlag;
1231cdf0e10cSrcweir ScMatrixRef pMat = pMat1;
1232cdf0e10cSrcweir if (!pMat)
1233cdf0e10cSrcweir {
1234cdf0e10cSrcweir fVal = fVal1;
1235cdf0e10cSrcweir pMat = pMat2;
1236cdf0e10cSrcweir bFlag = sal_True; // double - Matrix
1237cdf0e10cSrcweir }
1238cdf0e10cSrcweir else
1239cdf0e10cSrcweir {
1240cdf0e10cSrcweir fVal = fVal2;
1241cdf0e10cSrcweir bFlag = sal_False; // Matrix - double
1242cdf0e10cSrcweir }
1243cdf0e10cSrcweir SCSIZE nC, nR;
1244cdf0e10cSrcweir pMat->GetDimensions(nC, nR);
1245cdf0e10cSrcweir ScMatrixRef pResMat = GetNewMat(nC, nR);
1246cdf0e10cSrcweir if (pResMat)
1247cdf0e10cSrcweir {
1248cdf0e10cSrcweir SCSIZE nCount = nC * nR;
1249cdf0e10cSrcweir if (bFlag || !_bSub )
1250cdf0e10cSrcweir {
1251cdf0e10cSrcweir for ( SCSIZE i = 0; i < nCount; i++ )
1252cdf0e10cSrcweir {
1253cdf0e10cSrcweir if (pMat->IsValue(i))
1254cdf0e10cSrcweir pResMat->PutDouble( _bSub ? ::rtl::math::approxSub( fVal, pMat->GetDouble(i)) : ::rtl::math::approxAdd( pMat->GetDouble(i), fVal), i);
1255cdf0e10cSrcweir else
1256cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i);
1257cdf0e10cSrcweir } // for ( SCSIZE i = 0; i < nCount; i++ )
1258cdf0e10cSrcweir } // if (bFlag || !_bSub )
1259cdf0e10cSrcweir else
1260cdf0e10cSrcweir {
1261cdf0e10cSrcweir for ( SCSIZE i = 0; i < nCount; i++ )
1262cdf0e10cSrcweir { if (pMat->IsValue(i))
1263cdf0e10cSrcweir pResMat->PutDouble( ::rtl::math::approxSub( pMat->GetDouble(i), fVal), i);
1264cdf0e10cSrcweir else
1265cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i);
1266cdf0e10cSrcweir } // for ( SCSIZE i = 0; i < nCount; i++ )
1267cdf0e10cSrcweir }
1268cdf0e10cSrcweir PushMatrix(pResMat);
1269cdf0e10cSrcweir }
1270cdf0e10cSrcweir else
1271cdf0e10cSrcweir PushIllegalArgument();
1272cdf0e10cSrcweir }
1273cdf0e10cSrcweir else if ( _bSub )
1274cdf0e10cSrcweir PushDouble( ::rtl::math::approxSub( fVal1, fVal2 ) );
1275cdf0e10cSrcweir else
1276cdf0e10cSrcweir PushDouble( ::rtl::math::approxAdd( fVal1, fVal2 ) );
1277cdf0e10cSrcweir if ( nFmtCurrencyType == NUMBERFORMAT_CURRENCY )
1278cdf0e10cSrcweir {
1279cdf0e10cSrcweir nFuncFmtType = nFmtCurrencyType;
1280cdf0e10cSrcweir nFuncFmtIndex = nFmtCurrencyIndex;
1281cdf0e10cSrcweir }
1282cdf0e10cSrcweir else
1283cdf0e10cSrcweir {
1284cdf0e10cSrcweir lcl_GetDiffDateTimeFmtType( nFuncFmtType, nFmt1, nFmt2 );
1285cdf0e10cSrcweir if ( nFmtPercentType == NUMBERFORMAT_PERCENT && nFuncFmtType == NUMBERFORMAT_NUMBER )
1286cdf0e10cSrcweir nFuncFmtType = NUMBERFORMAT_PERCENT;
1287cdf0e10cSrcweir }
1288cdf0e10cSrcweir }
1289cdf0e10cSrcweir
ScAmpersand()1290cdf0e10cSrcweir void ScInterpreter::ScAmpersand()
1291cdf0e10cSrcweir {
1292cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScAmpersand" );
1293cdf0e10cSrcweir ScMatrixRef pMat1 = NULL;
1294cdf0e10cSrcweir ScMatrixRef pMat2 = NULL;
1295cdf0e10cSrcweir String sStr1, sStr2;
1296cdf0e10cSrcweir if ( GetStackType() == svMatrix )
1297cdf0e10cSrcweir pMat2 = GetMatrix();
1298cdf0e10cSrcweir else
1299cdf0e10cSrcweir sStr2 = GetString();
1300cdf0e10cSrcweir if ( GetStackType() == svMatrix )
1301cdf0e10cSrcweir pMat1 = GetMatrix();
1302cdf0e10cSrcweir else
1303cdf0e10cSrcweir sStr1 = GetString();
1304cdf0e10cSrcweir if (pMat1 && pMat2)
1305cdf0e10cSrcweir {
1306cdf0e10cSrcweir ScMatrixRef pResMat = MatConcat(pMat1, pMat2);
1307cdf0e10cSrcweir if (!pResMat)
1308cdf0e10cSrcweir PushNoValue();
1309cdf0e10cSrcweir else
1310cdf0e10cSrcweir PushMatrix(pResMat);
1311cdf0e10cSrcweir }
1312cdf0e10cSrcweir else if (pMat1 || pMat2)
1313cdf0e10cSrcweir {
1314cdf0e10cSrcweir String sStr;
1315cdf0e10cSrcweir sal_Bool bFlag;
1316cdf0e10cSrcweir ScMatrixRef pMat = pMat1;
1317cdf0e10cSrcweir if (!pMat)
1318cdf0e10cSrcweir {
1319cdf0e10cSrcweir sStr = sStr1;
1320cdf0e10cSrcweir pMat = pMat2;
1321cdf0e10cSrcweir bFlag = sal_True; // double - Matrix
1322cdf0e10cSrcweir }
1323cdf0e10cSrcweir else
1324cdf0e10cSrcweir {
1325cdf0e10cSrcweir sStr = sStr2;
1326cdf0e10cSrcweir bFlag = sal_False; // Matrix - double
1327cdf0e10cSrcweir }
1328cdf0e10cSrcweir SCSIZE nC, nR;
1329cdf0e10cSrcweir pMat->GetDimensions(nC, nR);
1330cdf0e10cSrcweir ScMatrixRef pResMat = GetNewMat(nC, nR);
1331cdf0e10cSrcweir if (pResMat)
1332cdf0e10cSrcweir {
1333cdf0e10cSrcweir SCSIZE nCount = nC * nR;
1334cdf0e10cSrcweir if (nGlobalError)
1335cdf0e10cSrcweir {
1336cdf0e10cSrcweir for ( SCSIZE i = 0; i < nCount; i++ )
1337cdf0e10cSrcweir pResMat->PutError( nGlobalError, i);
1338cdf0e10cSrcweir }
1339cdf0e10cSrcweir else if (bFlag)
1340cdf0e10cSrcweir {
1341cdf0e10cSrcweir for ( SCSIZE i = 0; i < nCount; i++ )
1342cdf0e10cSrcweir {
1343cdf0e10cSrcweir sal_uInt16 nErr = pMat->GetErrorIfNotString( i);
1344cdf0e10cSrcweir if (nErr)
1345cdf0e10cSrcweir pResMat->PutError( nErr, i);
1346cdf0e10cSrcweir else
1347cdf0e10cSrcweir {
1348cdf0e10cSrcweir String aTmp( sStr);
1349cdf0e10cSrcweir aTmp += pMat->GetString( *pFormatter, i);
1350cdf0e10cSrcweir pResMat->PutString( aTmp, i);
1351cdf0e10cSrcweir }
1352cdf0e10cSrcweir }
1353cdf0e10cSrcweir }
1354cdf0e10cSrcweir else
1355cdf0e10cSrcweir {
1356cdf0e10cSrcweir for ( SCSIZE i = 0; i < nCount; i++ )
1357cdf0e10cSrcweir {
1358cdf0e10cSrcweir sal_uInt16 nErr = pMat->GetErrorIfNotString( i);
1359cdf0e10cSrcweir if (nErr)
1360cdf0e10cSrcweir pResMat->PutError( nErr, i);
1361cdf0e10cSrcweir else
1362cdf0e10cSrcweir {
1363cdf0e10cSrcweir String aTmp( pMat->GetString( *pFormatter, i));
1364cdf0e10cSrcweir aTmp += sStr;
1365cdf0e10cSrcweir pResMat->PutString( aTmp, i);
1366cdf0e10cSrcweir }
1367cdf0e10cSrcweir }
1368cdf0e10cSrcweir }
1369cdf0e10cSrcweir PushMatrix(pResMat);
1370cdf0e10cSrcweir }
1371cdf0e10cSrcweir else
1372cdf0e10cSrcweir PushIllegalArgument();
1373cdf0e10cSrcweir }
1374cdf0e10cSrcweir else
1375cdf0e10cSrcweir {
1376cdf0e10cSrcweir if ( CheckStringResultLen( sStr1, sStr2 ) )
1377cdf0e10cSrcweir sStr1 += sStr2;
1378cdf0e10cSrcweir PushString(sStr1);
1379cdf0e10cSrcweir }
1380cdf0e10cSrcweir }
1381cdf0e10cSrcweir
ScSub()1382cdf0e10cSrcweir void ScInterpreter::ScSub()
1383cdf0e10cSrcweir {
1384cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSub" );
1385cdf0e10cSrcweir CalculateAddSub(sal_True);
1386cdf0e10cSrcweir }
1387cdf0e10cSrcweir
ScMul()1388cdf0e10cSrcweir void ScInterpreter::ScMul()
1389cdf0e10cSrcweir {
1390cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMul" );
1391cdf0e10cSrcweir ScMatrixRef pMat1 = NULL;
1392cdf0e10cSrcweir ScMatrixRef pMat2 = NULL;
1393cdf0e10cSrcweir double fVal1 = 0.0, fVal2 = 0.0;
1394cdf0e10cSrcweir short nFmtCurrencyType = nCurFmtType;
1395cdf0e10cSrcweir sal_uLong nFmtCurrencyIndex = nCurFmtIndex;
1396cdf0e10cSrcweir if ( GetStackType() == svMatrix )
1397cdf0e10cSrcweir pMat2 = GetMatrix();
1398cdf0e10cSrcweir else
1399cdf0e10cSrcweir {
1400cdf0e10cSrcweir fVal2 = GetDouble();
1401cdf0e10cSrcweir switch ( nCurFmtType )
1402cdf0e10cSrcweir {
1403cdf0e10cSrcweir case NUMBERFORMAT_CURRENCY :
1404cdf0e10cSrcweir nFmtCurrencyType = nCurFmtType;
1405cdf0e10cSrcweir nFmtCurrencyIndex = nCurFmtIndex;
1406cdf0e10cSrcweir break;
1407cdf0e10cSrcweir }
1408cdf0e10cSrcweir }
1409cdf0e10cSrcweir if ( GetStackType() == svMatrix )
1410cdf0e10cSrcweir pMat1 = GetMatrix();
1411cdf0e10cSrcweir else
1412cdf0e10cSrcweir {
1413cdf0e10cSrcweir fVal1 = GetDouble();
1414cdf0e10cSrcweir switch ( nCurFmtType )
1415cdf0e10cSrcweir {
1416cdf0e10cSrcweir case NUMBERFORMAT_CURRENCY :
1417cdf0e10cSrcweir nFmtCurrencyType = nCurFmtType;
1418cdf0e10cSrcweir nFmtCurrencyIndex = nCurFmtIndex;
1419cdf0e10cSrcweir break;
1420cdf0e10cSrcweir }
1421cdf0e10cSrcweir }
1422cdf0e10cSrcweir if (pMat1 && pMat2)
1423cdf0e10cSrcweir {
1424cdf0e10cSrcweir MatrixMul aMul;
1425cdf0e10cSrcweir ScMatrixRef pResMat = lcl_MatrixCalculation(aMul,pMat1, pMat2,this);
1426cdf0e10cSrcweir if (!pResMat)
1427cdf0e10cSrcweir PushNoValue();
1428cdf0e10cSrcweir else
1429cdf0e10cSrcweir PushMatrix(pResMat);
1430cdf0e10cSrcweir }
1431cdf0e10cSrcweir else if (pMat1 || pMat2)
1432cdf0e10cSrcweir {
1433cdf0e10cSrcweir double fVal;
1434cdf0e10cSrcweir ScMatrixRef pMat = pMat1;
1435cdf0e10cSrcweir if (!pMat)
1436cdf0e10cSrcweir {
1437cdf0e10cSrcweir fVal = fVal1;
1438cdf0e10cSrcweir pMat = pMat2;
1439cdf0e10cSrcweir }
1440cdf0e10cSrcweir else
1441cdf0e10cSrcweir fVal = fVal2;
1442cdf0e10cSrcweir SCSIZE nC, nR;
1443cdf0e10cSrcweir pMat->GetDimensions(nC, nR);
1444cdf0e10cSrcweir ScMatrixRef pResMat = GetNewMat(nC, nR);
1445cdf0e10cSrcweir if (pResMat)
1446cdf0e10cSrcweir {
1447cdf0e10cSrcweir SCSIZE nCount = nC * nR;
1448cdf0e10cSrcweir for ( SCSIZE i = 0; i < nCount; i++ )
1449cdf0e10cSrcweir if (pMat->IsValue(i))
1450cdf0e10cSrcweir pResMat->PutDouble(pMat->GetDouble(i)*fVal, i);
1451cdf0e10cSrcweir else
1452cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i);
1453cdf0e10cSrcweir PushMatrix(pResMat);
1454cdf0e10cSrcweir }
1455cdf0e10cSrcweir else
1456cdf0e10cSrcweir PushIllegalArgument();
1457cdf0e10cSrcweir }
1458cdf0e10cSrcweir else
1459cdf0e10cSrcweir PushDouble(fVal1 * fVal2);
1460cdf0e10cSrcweir if ( nFmtCurrencyType == NUMBERFORMAT_CURRENCY )
1461cdf0e10cSrcweir {
1462cdf0e10cSrcweir nFuncFmtType = nFmtCurrencyType;
1463cdf0e10cSrcweir nFuncFmtIndex = nFmtCurrencyIndex;
1464cdf0e10cSrcweir }
1465cdf0e10cSrcweir }
1466cdf0e10cSrcweir
ScDiv()1467cdf0e10cSrcweir void ScInterpreter::ScDiv()
1468cdf0e10cSrcweir {
1469cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScDiv" );
1470cdf0e10cSrcweir ScMatrixRef pMat1 = NULL;
1471cdf0e10cSrcweir ScMatrixRef pMat2 = NULL;
1472cdf0e10cSrcweir double fVal1 = 0.0, fVal2 = 0.0;
1473cdf0e10cSrcweir short nFmtCurrencyType = nCurFmtType;
1474cdf0e10cSrcweir sal_uLong nFmtCurrencyIndex = nCurFmtIndex;
1475cdf0e10cSrcweir short nFmtCurrencyType2 = NUMBERFORMAT_UNDEFINED;
1476cdf0e10cSrcweir if ( GetStackType() == svMatrix )
1477cdf0e10cSrcweir pMat2 = GetMatrix();
1478cdf0e10cSrcweir else
1479cdf0e10cSrcweir {
1480cdf0e10cSrcweir fVal2 = GetDouble();
1481cdf0e10cSrcweir // hier kein Currency uebernehmen, 123kg/456DM sind nicht DM
1482cdf0e10cSrcweir nFmtCurrencyType2 = nCurFmtType;
1483cdf0e10cSrcweir }
1484cdf0e10cSrcweir if ( GetStackType() == svMatrix )
1485cdf0e10cSrcweir pMat1 = GetMatrix();
1486cdf0e10cSrcweir else
1487cdf0e10cSrcweir {
1488cdf0e10cSrcweir fVal1 = GetDouble();
1489cdf0e10cSrcweir switch ( nCurFmtType )
1490cdf0e10cSrcweir {
1491cdf0e10cSrcweir case NUMBERFORMAT_CURRENCY :
1492cdf0e10cSrcweir nFmtCurrencyType = nCurFmtType;
1493cdf0e10cSrcweir nFmtCurrencyIndex = nCurFmtIndex;
1494cdf0e10cSrcweir break;
1495cdf0e10cSrcweir }
1496cdf0e10cSrcweir }
1497cdf0e10cSrcweir if (pMat1 && pMat2)
1498cdf0e10cSrcweir {
1499cdf0e10cSrcweir MatrixDiv aDiv;
1500cdf0e10cSrcweir ScMatrixRef pResMat = lcl_MatrixCalculation(aDiv,pMat1, pMat2,this);
1501cdf0e10cSrcweir if (!pResMat)
1502cdf0e10cSrcweir PushNoValue();
1503cdf0e10cSrcweir else
1504cdf0e10cSrcweir PushMatrix(pResMat);
1505cdf0e10cSrcweir }
1506cdf0e10cSrcweir else if (pMat1 || pMat2)
1507cdf0e10cSrcweir {
1508cdf0e10cSrcweir double fVal;
1509cdf0e10cSrcweir sal_Bool bFlag;
1510cdf0e10cSrcweir ScMatrixRef pMat = pMat1;
1511cdf0e10cSrcweir if (!pMat)
1512cdf0e10cSrcweir {
1513cdf0e10cSrcweir fVal = fVal1;
1514cdf0e10cSrcweir pMat = pMat2;
1515cdf0e10cSrcweir bFlag = sal_True; // double - Matrix
1516cdf0e10cSrcweir }
1517cdf0e10cSrcweir else
1518cdf0e10cSrcweir {
1519cdf0e10cSrcweir fVal = fVal2;
1520cdf0e10cSrcweir bFlag = sal_False; // Matrix - double
1521cdf0e10cSrcweir }
1522cdf0e10cSrcweir SCSIZE nC, nR;
1523cdf0e10cSrcweir pMat->GetDimensions(nC, nR);
1524cdf0e10cSrcweir ScMatrixRef pResMat = GetNewMat(nC, nR);
1525cdf0e10cSrcweir if (pResMat)
1526cdf0e10cSrcweir {
1527cdf0e10cSrcweir SCSIZE nCount = nC * nR;
1528cdf0e10cSrcweir if (bFlag)
1529cdf0e10cSrcweir { for ( SCSIZE i = 0; i < nCount; i++ )
1530cdf0e10cSrcweir if (pMat->IsValue(i))
1531cdf0e10cSrcweir pResMat->PutDouble( div( fVal, pMat->GetDouble(i)), i);
1532cdf0e10cSrcweir else
1533cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i);
1534cdf0e10cSrcweir }
1535cdf0e10cSrcweir else
1536cdf0e10cSrcweir { for ( SCSIZE i = 0; i < nCount; i++ )
1537cdf0e10cSrcweir if (pMat->IsValue(i))
1538cdf0e10cSrcweir pResMat->PutDouble( div( pMat->GetDouble(i), fVal), i);
1539cdf0e10cSrcweir else
1540cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i);
1541cdf0e10cSrcweir }
1542cdf0e10cSrcweir PushMatrix(pResMat);
1543cdf0e10cSrcweir }
1544cdf0e10cSrcweir else
1545cdf0e10cSrcweir PushIllegalArgument();
1546cdf0e10cSrcweir }
1547cdf0e10cSrcweir else
1548cdf0e10cSrcweir {
1549cdf0e10cSrcweir PushDouble( div( fVal1, fVal2) );
1550cdf0e10cSrcweir }
1551cdf0e10cSrcweir if ( nFmtCurrencyType == NUMBERFORMAT_CURRENCY && nFmtCurrencyType2 != NUMBERFORMAT_CURRENCY )
1552cdf0e10cSrcweir { // auch DM/DM ist nicht DM bzw. DEM/EUR nicht DEM
1553cdf0e10cSrcweir nFuncFmtType = nFmtCurrencyType;
1554cdf0e10cSrcweir nFuncFmtIndex = nFmtCurrencyIndex;
1555cdf0e10cSrcweir }
1556cdf0e10cSrcweir }
1557cdf0e10cSrcweir
ScPower()1558cdf0e10cSrcweir void ScInterpreter::ScPower()
1559cdf0e10cSrcweir {
1560cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScPower" );
1561cdf0e10cSrcweir if ( MustHaveParamCount( GetByte(), 2 ) )
1562cdf0e10cSrcweir ScPow();
1563cdf0e10cSrcweir }
1564cdf0e10cSrcweir
ScPow()1565cdf0e10cSrcweir void ScInterpreter::ScPow()
1566cdf0e10cSrcweir {
1567cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScPow" );
1568cdf0e10cSrcweir ScMatrixRef pMat1 = NULL;
1569cdf0e10cSrcweir ScMatrixRef pMat2 = NULL;
1570cdf0e10cSrcweir double fVal1 = 0.0, fVal2 = 0.0;
1571cdf0e10cSrcweir if ( GetStackType() == svMatrix )
1572cdf0e10cSrcweir pMat2 = GetMatrix();
1573cdf0e10cSrcweir else
1574cdf0e10cSrcweir fVal2 = GetDouble();
1575cdf0e10cSrcweir if ( GetStackType() == svMatrix )
1576cdf0e10cSrcweir pMat1 = GetMatrix();
1577cdf0e10cSrcweir else
1578cdf0e10cSrcweir fVal1 = GetDouble();
1579cdf0e10cSrcweir if (pMat1 && pMat2)
1580cdf0e10cSrcweir {
1581cdf0e10cSrcweir MatrixPow aPow;
1582cdf0e10cSrcweir ScMatrixRef pResMat = lcl_MatrixCalculation(aPow,pMat1, pMat2,this);
1583cdf0e10cSrcweir if (!pResMat)
1584cdf0e10cSrcweir PushNoValue();
1585cdf0e10cSrcweir else
1586cdf0e10cSrcweir PushMatrix(pResMat);
1587cdf0e10cSrcweir }
1588cdf0e10cSrcweir else if (pMat1 || pMat2)
1589cdf0e10cSrcweir {
1590cdf0e10cSrcweir double fVal;
1591cdf0e10cSrcweir sal_Bool bFlag;
1592cdf0e10cSrcweir ScMatrixRef pMat = pMat1;
1593cdf0e10cSrcweir if (!pMat)
1594cdf0e10cSrcweir {
1595cdf0e10cSrcweir fVal = fVal1;
1596cdf0e10cSrcweir pMat = pMat2;
1597cdf0e10cSrcweir bFlag = sal_True; // double - Matrix
1598cdf0e10cSrcweir }
1599cdf0e10cSrcweir else
1600cdf0e10cSrcweir {
1601cdf0e10cSrcweir fVal = fVal2;
1602cdf0e10cSrcweir bFlag = sal_False; // Matrix - double
1603cdf0e10cSrcweir }
1604cdf0e10cSrcweir SCSIZE nC, nR;
1605cdf0e10cSrcweir pMat->GetDimensions(nC, nR);
1606cdf0e10cSrcweir ScMatrixRef pResMat = GetNewMat(nC, nR);
1607cdf0e10cSrcweir if (pResMat)
1608cdf0e10cSrcweir {
1609cdf0e10cSrcweir SCSIZE nCount = nC * nR;
1610cdf0e10cSrcweir if (bFlag)
1611cdf0e10cSrcweir { for ( SCSIZE i = 0; i < nCount; i++ )
1612cdf0e10cSrcweir if (pMat->IsValue(i))
1613a067bd65SRob Weir pResMat->PutDouble(pow(fVal,pMat->GetDouble(i)), i);
1614cdf0e10cSrcweir else
1615cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i);
1616cdf0e10cSrcweir }
1617cdf0e10cSrcweir else
1618cdf0e10cSrcweir { for ( SCSIZE i = 0; i < nCount; i++ )
1619cdf0e10cSrcweir if (pMat->IsValue(i))
1620a067bd65SRob Weir pResMat->PutDouble(pow(pMat->GetDouble(i),fVal), i);
1621cdf0e10cSrcweir else
1622cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), i);
1623cdf0e10cSrcweir }
1624cdf0e10cSrcweir PushMatrix(pResMat);
1625cdf0e10cSrcweir }
1626cdf0e10cSrcweir else
1627cdf0e10cSrcweir PushIllegalArgument();
1628cdf0e10cSrcweir }
1629cdf0e10cSrcweir else
1630a067bd65SRob Weir PushDouble(pow(fVal1,fVal2));
1631cdf0e10cSrcweir }
1632cdf0e10cSrcweir
ScSumProduct()1633cdf0e10cSrcweir void ScInterpreter::ScSumProduct()
1634cdf0e10cSrcweir {
1635cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSumProduct" );
1636cdf0e10cSrcweir sal_uInt8 nParamCount = GetByte();
1637cdf0e10cSrcweir if ( !MustHaveParamCount( nParamCount, 1, 30 ) )
1638cdf0e10cSrcweir return;
1639cdf0e10cSrcweir
1640cdf0e10cSrcweir ScMatrixRef pMat1 = NULL;
1641cdf0e10cSrcweir ScMatrixRef pMat2 = NULL;
1642cdf0e10cSrcweir ScMatrixRef pMat = NULL;
1643cdf0e10cSrcweir pMat2 = GetMatrix();
1644cdf0e10cSrcweir if (!pMat2)
1645cdf0e10cSrcweir {
1646cdf0e10cSrcweir PushIllegalParameter();
1647cdf0e10cSrcweir return;
1648cdf0e10cSrcweir }
1649cdf0e10cSrcweir SCSIZE nC, nC1;
1650cdf0e10cSrcweir SCSIZE nR, nR1;
1651cdf0e10cSrcweir pMat2->GetDimensions(nC, nR);
1652cdf0e10cSrcweir pMat = pMat2;
1653cdf0e10cSrcweir MatrixMul aMul;
1654cdf0e10cSrcweir for (sal_uInt16 i = 1; i < nParamCount; i++)
1655cdf0e10cSrcweir {
1656cdf0e10cSrcweir pMat1 = GetMatrix();
1657cdf0e10cSrcweir if (!pMat1)
1658cdf0e10cSrcweir {
1659cdf0e10cSrcweir PushIllegalParameter();
1660cdf0e10cSrcweir return;
1661cdf0e10cSrcweir }
1662cdf0e10cSrcweir pMat1->GetDimensions(nC1, nR1);
1663cdf0e10cSrcweir if (nC1 != nC || nR1 != nR)
1664cdf0e10cSrcweir {
1665cdf0e10cSrcweir PushNoValue();
1666cdf0e10cSrcweir return;
1667cdf0e10cSrcweir }
1668cdf0e10cSrcweir ScMatrixRef pResMat = lcl_MatrixCalculation(aMul,pMat1, pMat,this);
1669cdf0e10cSrcweir if (!pResMat)
1670cdf0e10cSrcweir {
1671cdf0e10cSrcweir PushNoValue();
1672cdf0e10cSrcweir return;
1673cdf0e10cSrcweir }
1674cdf0e10cSrcweir else
1675cdf0e10cSrcweir pMat = pResMat;
1676cdf0e10cSrcweir }
1677cdf0e10cSrcweir double fSum = 0.0;
1678cdf0e10cSrcweir SCSIZE nCount = pMat->GetElementCount();
1679cdf0e10cSrcweir for (SCSIZE j = 0; j < nCount; j++)
1680cdf0e10cSrcweir {
1681cdf0e10cSrcweir if (!pMat->IsString(j))
1682cdf0e10cSrcweir fSum += pMat->GetDouble(j);
1683cdf0e10cSrcweir }
1684cdf0e10cSrcweir PushDouble(fSum);
1685cdf0e10cSrcweir }
1686cdf0e10cSrcweir
ScSumX2MY2()1687cdf0e10cSrcweir void ScInterpreter::ScSumX2MY2()
1688cdf0e10cSrcweir {
1689cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSumX2MY2" );
1690cdf0e10cSrcweir CalculateSumX2MY2SumX2DY2(sal_False);
1691cdf0e10cSrcweir }
CalculateSumX2MY2SumX2DY2(sal_Bool _bSumX2DY2)1692cdf0e10cSrcweir void ScInterpreter::CalculateSumX2MY2SumX2DY2(sal_Bool _bSumX2DY2)
1693cdf0e10cSrcweir {
1694cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CalculateSumX2MY2SumX2DY2" );
1695cdf0e10cSrcweir if ( !MustHaveParamCount( GetByte(), 2 ) )
1696cdf0e10cSrcweir return;
1697cdf0e10cSrcweir
1698cdf0e10cSrcweir ScMatrixRef pMat1 = NULL;
1699cdf0e10cSrcweir ScMatrixRef pMat2 = NULL;
1700cdf0e10cSrcweir SCSIZE i, j;
1701cdf0e10cSrcweir pMat2 = GetMatrix();
1702cdf0e10cSrcweir pMat1 = GetMatrix();
1703cdf0e10cSrcweir if (!pMat2 || !pMat1)
1704cdf0e10cSrcweir {
1705cdf0e10cSrcweir PushIllegalParameter();
1706cdf0e10cSrcweir return;
1707cdf0e10cSrcweir }
1708cdf0e10cSrcweir SCSIZE nC1, nC2;
1709cdf0e10cSrcweir SCSIZE nR1, nR2;
1710cdf0e10cSrcweir pMat2->GetDimensions(nC2, nR2);
1711cdf0e10cSrcweir pMat1->GetDimensions(nC1, nR1);
1712cdf0e10cSrcweir if (nC1 != nC2 || nR1 != nR2)
1713cdf0e10cSrcweir {
1714cdf0e10cSrcweir PushNoValue();
1715cdf0e10cSrcweir return;
1716cdf0e10cSrcweir }
1717cdf0e10cSrcweir double fVal, fSum = 0.0;
1718cdf0e10cSrcweir for (i = 0; i < nC1; i++)
1719cdf0e10cSrcweir for (j = 0; j < nR1; j++)
1720cdf0e10cSrcweir if (!pMat1->IsString(i,j) && !pMat2->IsString(i,j))
1721cdf0e10cSrcweir {
1722cdf0e10cSrcweir fVal = pMat1->GetDouble(i,j);
1723cdf0e10cSrcweir fSum += fVal * fVal;
1724cdf0e10cSrcweir fVal = pMat2->GetDouble(i,j);
1725cdf0e10cSrcweir if ( _bSumX2DY2 )
1726cdf0e10cSrcweir fSum += fVal * fVal;
1727cdf0e10cSrcweir else
1728cdf0e10cSrcweir fSum -= fVal * fVal;
1729cdf0e10cSrcweir }
1730cdf0e10cSrcweir PushDouble(fSum);
1731cdf0e10cSrcweir }
1732cdf0e10cSrcweir
ScSumX2DY2()1733cdf0e10cSrcweir void ScInterpreter::ScSumX2DY2()
1734cdf0e10cSrcweir {
1735cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSumX2DY2" );
1736cdf0e10cSrcweir CalculateSumX2MY2SumX2DY2(sal_True);
1737cdf0e10cSrcweir }
1738cdf0e10cSrcweir
ScSumXMY2()1739cdf0e10cSrcweir void ScInterpreter::ScSumXMY2()
1740cdf0e10cSrcweir {
1741cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScSumXMY2" );
1742cdf0e10cSrcweir if ( !MustHaveParamCount( GetByte(), 2 ) )
1743cdf0e10cSrcweir return;
1744cdf0e10cSrcweir
1745cdf0e10cSrcweir ScMatrixRef pMat1 = NULL;
1746cdf0e10cSrcweir ScMatrixRef pMat2 = NULL;
1747cdf0e10cSrcweir pMat2 = GetMatrix();
1748cdf0e10cSrcweir pMat1 = GetMatrix();
1749cdf0e10cSrcweir if (!pMat2 || !pMat1)
1750cdf0e10cSrcweir {
1751cdf0e10cSrcweir PushIllegalParameter();
1752cdf0e10cSrcweir return;
1753cdf0e10cSrcweir }
1754cdf0e10cSrcweir SCSIZE nC1, nC2;
1755cdf0e10cSrcweir SCSIZE nR1, nR2;
1756cdf0e10cSrcweir pMat2->GetDimensions(nC2, nR2);
1757cdf0e10cSrcweir pMat1->GetDimensions(nC1, nR1);
1758cdf0e10cSrcweir if (nC1 != nC2 || nR1 != nR2)
1759cdf0e10cSrcweir {
1760cdf0e10cSrcweir PushNoValue();
1761cdf0e10cSrcweir return;
1762cdf0e10cSrcweir } // if (nC1 != nC2 || nR1 != nR2)
1763cdf0e10cSrcweir MatrixSub aSub;
1764cdf0e10cSrcweir ScMatrixRef pResMat = lcl_MatrixCalculation(aSub,pMat1, pMat2,this);
1765cdf0e10cSrcweir if (!pResMat)
1766cdf0e10cSrcweir {
1767cdf0e10cSrcweir PushNoValue();
1768cdf0e10cSrcweir }
1769cdf0e10cSrcweir else
1770cdf0e10cSrcweir {
1771cdf0e10cSrcweir double fVal, fSum = 0.0;
1772cdf0e10cSrcweir SCSIZE nCount = pResMat->GetElementCount();
1773cdf0e10cSrcweir for (SCSIZE i = 0; i < nCount; i++)
1774cdf0e10cSrcweir if (!pResMat->IsString(i))
1775cdf0e10cSrcweir {
1776cdf0e10cSrcweir fVal = pResMat->GetDouble(i);
1777cdf0e10cSrcweir fSum += fVal * fVal;
1778cdf0e10cSrcweir }
1779cdf0e10cSrcweir PushDouble(fSum);
1780cdf0e10cSrcweir }
1781cdf0e10cSrcweir }
1782cdf0e10cSrcweir
ScFrequency()1783cdf0e10cSrcweir void ScInterpreter::ScFrequency()
1784cdf0e10cSrcweir {
1785cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScFrequency" );
1786cdf0e10cSrcweir if ( !MustHaveParamCount( GetByte(), 2 ) )
1787cdf0e10cSrcweir return;
1788cdf0e10cSrcweir
1789cdf0e10cSrcweir vector<double> aBinArray;
1790cdf0e10cSrcweir vector<long> aBinIndexOrder;
1791cdf0e10cSrcweir
1792cdf0e10cSrcweir GetSortArray(1, aBinArray, &aBinIndexOrder);
1793cdf0e10cSrcweir SCSIZE nBinSize = aBinArray.size();
1794cdf0e10cSrcweir if (nGlobalError)
1795cdf0e10cSrcweir {
1796cdf0e10cSrcweir PushNoValue();
1797cdf0e10cSrcweir return;
1798cdf0e10cSrcweir }
1799cdf0e10cSrcweir
1800cdf0e10cSrcweir vector<double> aDataArray;
1801cdf0e10cSrcweir GetSortArray(1, aDataArray);
1802cdf0e10cSrcweir SCSIZE nDataSize = aDataArray.size();
1803cdf0e10cSrcweir
1804cdf0e10cSrcweir if (aDataArray.empty() || nGlobalError)
1805cdf0e10cSrcweir {
1806cdf0e10cSrcweir PushNoValue();
1807cdf0e10cSrcweir return;
1808cdf0e10cSrcweir }
1809cdf0e10cSrcweir ScMatrixRef pResMat = GetNewMat(1, nBinSize+1);
1810cdf0e10cSrcweir if (!pResMat)
1811cdf0e10cSrcweir {
1812cdf0e10cSrcweir PushIllegalArgument();
1813cdf0e10cSrcweir return;
1814cdf0e10cSrcweir }
1815cdf0e10cSrcweir
1816cdf0e10cSrcweir if (nBinSize != aBinIndexOrder.size())
1817cdf0e10cSrcweir {
1818cdf0e10cSrcweir PushIllegalArgument();
1819cdf0e10cSrcweir return;
1820cdf0e10cSrcweir }
1821cdf0e10cSrcweir
1822cdf0e10cSrcweir SCSIZE j;
1823cdf0e10cSrcweir SCSIZE i = 0;
1824cdf0e10cSrcweir for (j = 0; j < nBinSize; ++j)
1825cdf0e10cSrcweir {
1826cdf0e10cSrcweir SCSIZE nCount = 0;
1827cdf0e10cSrcweir while (i < nDataSize && aDataArray[i] <= aBinArray[j])
1828cdf0e10cSrcweir {
1829cdf0e10cSrcweir ++nCount;
1830cdf0e10cSrcweir ++i;
1831cdf0e10cSrcweir }
1832cdf0e10cSrcweir pResMat->PutDouble(static_cast<double>(nCount), aBinIndexOrder[j]);
1833cdf0e10cSrcweir }
1834cdf0e10cSrcweir pResMat->PutDouble(static_cast<double>(nDataSize-i), j);
1835cdf0e10cSrcweir PushMatrix(pResMat);
1836cdf0e10cSrcweir }
1837cdf0e10cSrcweir
1838cdf0e10cSrcweir // -----------------------------------------------------------------------------
1839cdf0e10cSrcweir // Helper methods for LINEST/LOGEST and TREND/GROWTH
1840cdf0e10cSrcweir // All matrices must already exist and have the needed size, no control tests
1841cdf0e10cSrcweir // done. Those methodes, which names start with lcl_T, are adapted to case 3,
1842cdf0e10cSrcweir // where Y (=observed values) is given as row.
1843cdf0e10cSrcweir // Remember, ScMatrix matrices are zero based, index access (column,row).
1844cdf0e10cSrcweir // -----------------------------------------------------------------------------
1845cdf0e10cSrcweir
1846cdf0e10cSrcweir // Multiply n x m Mat A with m x l Mat B to n x l Mat R
lcl_MFastMult(ScMatrixRef pA,ScMatrixRef pB,ScMatrixRef pR,SCSIZE n,SCSIZE m,SCSIZE l)1847cdf0e10cSrcweir void lcl_MFastMult( ScMatrixRef pA, ScMatrixRef pB, ScMatrixRef pR, SCSIZE n, SCSIZE m, SCSIZE l )
1848cdf0e10cSrcweir {
1849cdf0e10cSrcweir double sum;
1850cdf0e10cSrcweir for (SCSIZE row = 0; row < n; row++)
1851cdf0e10cSrcweir {
1852cdf0e10cSrcweir for (SCSIZE col = 0; col < l; col++)
1853cdf0e10cSrcweir { // result element(col, row) =sum[ (row of A) * (column of B)]
1854cdf0e10cSrcweir sum = 0.0;
1855cdf0e10cSrcweir for (SCSIZE k = 0; k < m; k++)
1856cdf0e10cSrcweir sum += pA->GetDouble(k,row) * pB->GetDouble(col,k);
1857cdf0e10cSrcweir pR->PutDouble(sum, col, row);
1858cdf0e10cSrcweir }
1859cdf0e10cSrcweir }
1860cdf0e10cSrcweir }
1861cdf0e10cSrcweir
1862cdf0e10cSrcweir // <A;B> over all elements; uses the matrices as vectors of length M
lcl_GetSumProduct(ScMatrixRef pMatA,ScMatrixRef pMatB,SCSIZE nM)1863cdf0e10cSrcweir double lcl_GetSumProduct( ScMatrixRef pMatA, ScMatrixRef pMatB, SCSIZE nM )
1864cdf0e10cSrcweir {
1865cdf0e10cSrcweir double fSum = 0.0;
1866cdf0e10cSrcweir for (SCSIZE i=0; i<nM; i++)
1867cdf0e10cSrcweir fSum += pMatA->GetDouble(i) * pMatB->GetDouble(i);
1868cdf0e10cSrcweir return fSum;
1869cdf0e10cSrcweir }
1870cdf0e10cSrcweir
1871cdf0e10cSrcweir // Special version for use within QR decomposition.
1872cdf0e10cSrcweir // Euclidean norm of column index C starting in row index R;
1873cdf0e10cSrcweir // matrix A has count N rows.
lcl_GetColumnEuclideanNorm(ScMatrixRef pMatA,SCSIZE nC,SCSIZE nR,SCSIZE nN)1874cdf0e10cSrcweir double lcl_GetColumnEuclideanNorm( ScMatrixRef pMatA, SCSIZE nC, SCSIZE nR, SCSIZE nN )
1875cdf0e10cSrcweir {
1876cdf0e10cSrcweir double fNorm = 0.0;
1877cdf0e10cSrcweir for (SCSIZE row=nR; row<nN; row++)
1878cdf0e10cSrcweir fNorm += (pMatA->GetDouble(nC,row)) * (pMatA->GetDouble(nC,row));
1879cdf0e10cSrcweir return sqrt(fNorm);
1880cdf0e10cSrcweir }
1881cdf0e10cSrcweir
1882cdf0e10cSrcweir // Euclidean norm of row index R starting in column index C;
1883cdf0e10cSrcweir // matrix A has count N columns.
lcl_TGetColumnEuclideanNorm(ScMatrixRef pMatA,SCSIZE nR,SCSIZE nC,SCSIZE nN)1884cdf0e10cSrcweir double lcl_TGetColumnEuclideanNorm( ScMatrixRef pMatA, SCSIZE nR, SCSIZE nC, SCSIZE nN )
1885cdf0e10cSrcweir {
1886cdf0e10cSrcweir double fNorm = 0.0;
1887cdf0e10cSrcweir for (SCSIZE col=nC; col<nN; col++)
1888cdf0e10cSrcweir fNorm += (pMatA->GetDouble(col,nR)) * (pMatA->GetDouble(col,nR));
1889cdf0e10cSrcweir return sqrt(fNorm);
1890cdf0e10cSrcweir }
1891cdf0e10cSrcweir
1892cdf0e10cSrcweir // Special version for use within QR decomposition.
1893cdf0e10cSrcweir // Maximum norm of column index C starting in row index R;
1894cdf0e10cSrcweir // matrix A has count N rows.
lcl_GetColumnMaximumNorm(ScMatrixRef pMatA,SCSIZE nC,SCSIZE nR,SCSIZE nN)1895cdf0e10cSrcweir double lcl_GetColumnMaximumNorm( ScMatrixRef pMatA, SCSIZE nC, SCSIZE nR, SCSIZE nN )
1896cdf0e10cSrcweir {
1897cdf0e10cSrcweir double fNorm = 0.0;
1898cdf0e10cSrcweir for (SCSIZE row=nR; row<nN; row++)
1899cdf0e10cSrcweir if (fNorm < fabs(pMatA->GetDouble(nC,row)))
1900cdf0e10cSrcweir fNorm = fabs(pMatA->GetDouble(nC,row));
1901cdf0e10cSrcweir return fNorm;
1902cdf0e10cSrcweir }
1903cdf0e10cSrcweir
1904cdf0e10cSrcweir // Maximum norm of row index R starting in col index C;
1905cdf0e10cSrcweir // matrix A has count N columns.
lcl_TGetColumnMaximumNorm(ScMatrixRef pMatA,SCSIZE nR,SCSIZE nC,SCSIZE nN)1906cdf0e10cSrcweir double lcl_TGetColumnMaximumNorm( ScMatrixRef pMatA, SCSIZE nR, SCSIZE nC, SCSIZE nN )
1907cdf0e10cSrcweir {
1908cdf0e10cSrcweir double fNorm = 0.0;
1909cdf0e10cSrcweir for (SCSIZE col=nC; col<nN; col++)
1910cdf0e10cSrcweir if (fNorm < fabs(pMatA->GetDouble(col,nR)))
1911cdf0e10cSrcweir fNorm = fabs(pMatA->GetDouble(col,nR));
1912cdf0e10cSrcweir return fNorm;
1913cdf0e10cSrcweir }
1914cdf0e10cSrcweir
1915cdf0e10cSrcweir // Special version for use within QR decomposition.
1916cdf0e10cSrcweir // <A(Ca);B(Cb)> starting in row index R;
1917cdf0e10cSrcweir // Ca and Cb are indices of columns, matrices A and B have count N rows.
lcl_GetColumnSumProduct(ScMatrixRef pMatA,SCSIZE nCa,ScMatrixRef pMatB,SCSIZE nCb,SCSIZE nR,SCSIZE nN)1918cdf0e10cSrcweir double lcl_GetColumnSumProduct( ScMatrixRef pMatA, SCSIZE nCa, ScMatrixRef pMatB,
1919cdf0e10cSrcweir SCSIZE nCb, SCSIZE nR, SCSIZE nN )
1920cdf0e10cSrcweir {
1921cdf0e10cSrcweir double fResult = 0.0;
1922cdf0e10cSrcweir for (SCSIZE row=nR; row<nN; row++)
1923cdf0e10cSrcweir fResult += pMatA->GetDouble(nCa,row) * pMatB->GetDouble(nCb,row);
1924cdf0e10cSrcweir return fResult;
1925cdf0e10cSrcweir }
1926cdf0e10cSrcweir
1927cdf0e10cSrcweir // <A(Ra);B(Rb)> starting in column index C;
1928cdf0e10cSrcweir // Ra and Rb are indices of rows, matrices A and B have count N columns.
lcl_TGetColumnSumProduct(ScMatrixRef pMatA,SCSIZE nRa,ScMatrixRef pMatB,SCSIZE nRb,SCSIZE nC,SCSIZE nN)1929cdf0e10cSrcweir double lcl_TGetColumnSumProduct( ScMatrixRef pMatA, SCSIZE nRa,
1930cdf0e10cSrcweir ScMatrixRef pMatB, SCSIZE nRb, SCSIZE nC, SCSIZE nN )
1931cdf0e10cSrcweir {
1932cdf0e10cSrcweir double fResult = 0.0;
1933cdf0e10cSrcweir for (SCSIZE col=nC; col<nN; col++)
1934cdf0e10cSrcweir fResult += pMatA->GetDouble(col,nRa) * pMatB->GetDouble(col,nRb);
1935cdf0e10cSrcweir return fResult;
1936cdf0e10cSrcweir }
1937cdf0e10cSrcweir
193886e1cf34SPedro Giffuni // #118029# no mathematical signum, but used to switch between adding and subtracting
lcl_GetSign(double fValue)1939cdf0e10cSrcweir double lcl_GetSign(double fValue)
1940cdf0e10cSrcweir {
194186e1cf34SPedro Giffuni return (fValue >= 0.0 ? 1.0 : -1.0 );
1942cdf0e10cSrcweir }
1943cdf0e10cSrcweir
1944cdf0e10cSrcweir /* Calculates a QR decomposition with Householder reflection.
1945cdf0e10cSrcweir * For each NxK matrix A exists a decomposition A=Q*R with an orthogonal
1946cdf0e10cSrcweir * NxN matrix Q and a NxK matrix R.
1947cdf0e10cSrcweir * Q=H1*H2*...*Hk with Householder matrices H. Such a householder matrix can
1948cdf0e10cSrcweir * be build from a vector u by H=I-(2/u'u)*(u u'). This vectors u are returned
1949cdf0e10cSrcweir * in the columns of matrix A, overwriting the old content.
1950cdf0e10cSrcweir * The matrix R has a quadric upper part KxK with values in the upper right
1951cdf0e10cSrcweir * triangle and zeros in all other elements. Here the diagonal elements of R
1952cdf0e10cSrcweir * are stored in the vector R and the other upper right elements in the upper
1953cdf0e10cSrcweir * right of the matrix A.
1954cdf0e10cSrcweir * The function returns false, if calculation breaks. But because of round-off
1955cdf0e10cSrcweir * errors singularity is often not detected.
1956cdf0e10cSrcweir */
lcl_CalculateQRdecomposition(ScMatrixRef pMatA,::std::vector<double> & pVecR,SCSIZE nK,SCSIZE nN)1957cdf0e10cSrcweir bool lcl_CalculateQRdecomposition(ScMatrixRef pMatA,
1958cdf0e10cSrcweir ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN)
1959cdf0e10cSrcweir {
1960cdf0e10cSrcweir double fScale ;
1961cdf0e10cSrcweir double fEuclid ;
1962cdf0e10cSrcweir double fFactor ;
1963cdf0e10cSrcweir double fSignum ;
1964cdf0e10cSrcweir double fSum ;
1965cdf0e10cSrcweir // ScMatrix matrices are zero based, index access (column,row)
1966cdf0e10cSrcweir for (SCSIZE col = 0; col <nK; col++)
1967cdf0e10cSrcweir {
1968cdf0e10cSrcweir // calculate vector u of the householder transformation
1969cdf0e10cSrcweir fScale = lcl_GetColumnMaximumNorm(pMatA, col, col, nN);
1970cdf0e10cSrcweir if (fScale == 0.0)
1971cdf0e10cSrcweir // A is singular
1972cdf0e10cSrcweir return false;
1973cdf0e10cSrcweir
1974cdf0e10cSrcweir for (SCSIZE row = col; row <nN; row++)
1975cdf0e10cSrcweir pMatA->PutDouble( pMatA->GetDouble(col,row)/fScale, col, row);
1976cdf0e10cSrcweir
1977cdf0e10cSrcweir fEuclid = lcl_GetColumnEuclideanNorm(pMatA, col, col, nN);
1978cdf0e10cSrcweir fFactor = 1.0/fEuclid/(fEuclid + fabs(pMatA->GetDouble(col,col)));
1979cdf0e10cSrcweir fSignum = lcl_GetSign(pMatA->GetDouble(col,col));
1980cdf0e10cSrcweir pMatA->PutDouble( pMatA->GetDouble(col,col) + fSignum*fEuclid, col,col);
1981cdf0e10cSrcweir pVecR[col] = -fSignum * fScale * fEuclid;
1982cdf0e10cSrcweir
1983cdf0e10cSrcweir // apply Householder transformation to A
1984cdf0e10cSrcweir for (SCSIZE c=col+1; c<nK; c++)
1985cdf0e10cSrcweir {
1986cdf0e10cSrcweir fSum =lcl_GetColumnSumProduct(pMatA, col, pMatA, c, col, nN);
1987cdf0e10cSrcweir for (SCSIZE row = col; row <nN; row++)
1988cdf0e10cSrcweir pMatA->PutDouble( pMatA->GetDouble(c,row)
1989cdf0e10cSrcweir - fSum * fFactor * pMatA->GetDouble(col,row), c, row);
1990cdf0e10cSrcweir }
1991cdf0e10cSrcweir }
1992cdf0e10cSrcweir return true;
1993cdf0e10cSrcweir }
1994cdf0e10cSrcweir
1995cdf0e10cSrcweir // same with transposed matrix A, N is count of columns, K count of rows
lcl_TCalculateQRdecomposition(ScMatrixRef pMatA,::std::vector<double> & pVecR,SCSIZE nK,SCSIZE nN)1996cdf0e10cSrcweir bool lcl_TCalculateQRdecomposition(ScMatrixRef pMatA,
1997cdf0e10cSrcweir ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN)
1998cdf0e10cSrcweir {
1999cdf0e10cSrcweir double fScale ;
2000cdf0e10cSrcweir double fEuclid ;
2001cdf0e10cSrcweir double fFactor ;
2002cdf0e10cSrcweir double fSignum ;
2003cdf0e10cSrcweir double fSum ;
2004cdf0e10cSrcweir // ScMatrix matrices are zero based, index access (column,row)
2005cdf0e10cSrcweir for (SCSIZE row = 0; row <nK; row++)
2006cdf0e10cSrcweir {
2007cdf0e10cSrcweir // calculate vector u of the householder transformation
2008cdf0e10cSrcweir fScale = lcl_TGetColumnMaximumNorm(pMatA, row, row, nN);
2009cdf0e10cSrcweir if (fScale == 0.0)
2010cdf0e10cSrcweir // A is singular
2011cdf0e10cSrcweir return false;
2012cdf0e10cSrcweir
2013cdf0e10cSrcweir for (SCSIZE col = row; col <nN; col++)
2014cdf0e10cSrcweir pMatA->PutDouble( pMatA->GetDouble(col,row)/fScale, col, row);
2015cdf0e10cSrcweir
2016cdf0e10cSrcweir fEuclid = lcl_TGetColumnEuclideanNorm(pMatA, row, row, nN);
2017cdf0e10cSrcweir fFactor = 1.0/fEuclid/(fEuclid + fabs(pMatA->GetDouble(row,row)));
2018cdf0e10cSrcweir fSignum = lcl_GetSign(pMatA->GetDouble(row,row));
2019cdf0e10cSrcweir pMatA->PutDouble( pMatA->GetDouble(row,row) + fSignum*fEuclid, row,row);
2020cdf0e10cSrcweir pVecR[row] = -fSignum * fScale * fEuclid;
2021cdf0e10cSrcweir
2022cdf0e10cSrcweir // apply Householder transformation to A
2023cdf0e10cSrcweir for (SCSIZE r=row+1; r<nK; r++)
2024cdf0e10cSrcweir {
2025cdf0e10cSrcweir fSum =lcl_TGetColumnSumProduct(pMatA, row, pMatA, r, row, nN);
2026cdf0e10cSrcweir for (SCSIZE col = row; col <nN; col++)
2027cdf0e10cSrcweir pMatA->PutDouble( pMatA->GetDouble(col,r)
2028cdf0e10cSrcweir - fSum * fFactor * pMatA->GetDouble(col,row), col, r);
2029cdf0e10cSrcweir }
2030cdf0e10cSrcweir }
2031cdf0e10cSrcweir return true;
2032cdf0e10cSrcweir }
2033cdf0e10cSrcweir
2034cdf0e10cSrcweir
2035cdf0e10cSrcweir /* Applies a Householder transformation to a column vector Y with is given as
2036cdf0e10cSrcweir * Nx1 Matrix. The Vektor u, from which the Householder transformation is build,
2037cdf0e10cSrcweir * is the column part in matrix A, with column index C, starting with row
2038cdf0e10cSrcweir * index C. A is the result of the QR decomposition as obtained from
2039cdf0e10cSrcweir * lcl_CaluclateQRdecomposition.
2040cdf0e10cSrcweir */
lcl_ApplyHouseholderTransformation(ScMatrixRef pMatA,SCSIZE nC,ScMatrixRef pMatY,SCSIZE nN)2041cdf0e10cSrcweir void lcl_ApplyHouseholderTransformation(ScMatrixRef pMatA, SCSIZE nC,
2042cdf0e10cSrcweir ScMatrixRef pMatY, SCSIZE nN)
2043cdf0e10cSrcweir {
2044cdf0e10cSrcweir // ScMatrix matrices are zero based, index access (column,row)
2045cdf0e10cSrcweir double fDenominator = lcl_GetColumnSumProduct(pMatA, nC, pMatA, nC, nC, nN);
2046cdf0e10cSrcweir double fNumerator = lcl_GetColumnSumProduct(pMatA, nC, pMatY, 0, nC, nN);
2047cdf0e10cSrcweir double fFactor = 2.0 * (fNumerator/fDenominator);
2048cdf0e10cSrcweir for (SCSIZE row = nC; row < nN; row++)
2049cdf0e10cSrcweir pMatY->PutDouble(
2050cdf0e10cSrcweir pMatY->GetDouble(row) - fFactor * pMatA->GetDouble(nC,row), row);
2051cdf0e10cSrcweir }
2052cdf0e10cSrcweir
2053cdf0e10cSrcweir // Same with transposed matrices A and Y.
lcl_TApplyHouseholderTransformation(ScMatrixRef pMatA,SCSIZE nR,ScMatrixRef pMatY,SCSIZE nN)2054cdf0e10cSrcweir void lcl_TApplyHouseholderTransformation(ScMatrixRef pMatA, SCSIZE nR,
2055cdf0e10cSrcweir ScMatrixRef pMatY, SCSIZE nN)
2056cdf0e10cSrcweir {
2057cdf0e10cSrcweir // ScMatrix matrices are zero based, index access (column,row)
2058cdf0e10cSrcweir double fDenominator = lcl_TGetColumnSumProduct(pMatA, nR, pMatA, nR, nR, nN);
2059cdf0e10cSrcweir double fNumerator = lcl_TGetColumnSumProduct(pMatA, nR, pMatY, 0, nR, nN);
2060cdf0e10cSrcweir double fFactor = 2.0 * (fNumerator/fDenominator);
2061cdf0e10cSrcweir for (SCSIZE col = nR; col < nN; col++)
2062cdf0e10cSrcweir pMatY->PutDouble(
2063cdf0e10cSrcweir pMatY->GetDouble(col) - fFactor * pMatA->GetDouble(col,nR), col);
2064cdf0e10cSrcweir }
2065cdf0e10cSrcweir
2066cdf0e10cSrcweir /* Solve for X in R*X=S using back substitution. The solution X overwrites S.
2067cdf0e10cSrcweir * Uses R from the result of the QR decomposition of a NxK matrix A.
2068cdf0e10cSrcweir * S is a column vector given as matrix, with at least elements on index
2069cdf0e10cSrcweir * 0 to K-1; elements on index>=K are ignored. Vector R must not have zero
2070cdf0e10cSrcweir * elements, no check is done.
2071cdf0e10cSrcweir */
lcl_SolveWithUpperRightTriangle(ScMatrixRef pMatA,::std::vector<double> & pVecR,ScMatrixRef pMatS,SCSIZE nK,bool bIsTransposed)2072cdf0e10cSrcweir void lcl_SolveWithUpperRightTriangle(ScMatrixRef pMatA,
2073cdf0e10cSrcweir ::std::vector< double>& pVecR, ScMatrixRef pMatS,
2074cdf0e10cSrcweir SCSIZE nK, bool bIsTransposed)
2075cdf0e10cSrcweir {
2076cdf0e10cSrcweir // ScMatrix matrices are zero based, index access (column,row)
2077cdf0e10cSrcweir double fSum;
2078cdf0e10cSrcweir SCSIZE row;
2079cdf0e10cSrcweir // SCSIZE is never negative, therefore test with rowp1=row+1
2080cdf0e10cSrcweir for (SCSIZE rowp1 = nK; rowp1>0; rowp1--)
2081cdf0e10cSrcweir {
2082cdf0e10cSrcweir row = rowp1-1;
2083cdf0e10cSrcweir fSum = pMatS->GetDouble(row);
2084cdf0e10cSrcweir for (SCSIZE col = rowp1; col<nK ; col++)
2085cdf0e10cSrcweir if (bIsTransposed)
2086cdf0e10cSrcweir fSum -= pMatA->GetDouble(row,col) * pMatS->GetDouble(col);
2087cdf0e10cSrcweir else
2088cdf0e10cSrcweir fSum -= pMatA->GetDouble(col,row) * pMatS->GetDouble(col);
2089cdf0e10cSrcweir pMatS->PutDouble( fSum / pVecR[row] , row);
2090cdf0e10cSrcweir }
2091cdf0e10cSrcweir }
2092cdf0e10cSrcweir
2093cdf0e10cSrcweir /* Solve for X in R' * X= T using forward substitution. The solution X
2094cdf0e10cSrcweir * overwrites T. Uses R from the result of the QR decomposition of a NxK
2095cdf0e10cSrcweir * matrix A. T is a column vectors given as matrix, with at least elements on
2096cdf0e10cSrcweir * index 0 to K-1; elements on index>=K are ignored. Vector R must not have
2097cdf0e10cSrcweir * zero elements, no check is done.
2098cdf0e10cSrcweir */
lcl_SolveWithLowerLeftTriangle(ScMatrixRef pMatA,::std::vector<double> & pVecR,ScMatrixRef pMatT,SCSIZE nK,bool bIsTransposed)2099cdf0e10cSrcweir void lcl_SolveWithLowerLeftTriangle(ScMatrixRef pMatA,
2100cdf0e10cSrcweir ::std::vector< double>& pVecR, ScMatrixRef pMatT,
2101cdf0e10cSrcweir SCSIZE nK, bool bIsTransposed)
2102cdf0e10cSrcweir {
2103cdf0e10cSrcweir // ScMatrix matrices are zero based, index access (column,row)
2104cdf0e10cSrcweir double fSum;
2105cdf0e10cSrcweir for (SCSIZE row = 0; row < nK; row++)
2106cdf0e10cSrcweir {
2107cdf0e10cSrcweir fSum = pMatT -> GetDouble(row);
2108cdf0e10cSrcweir for (SCSIZE col=0; col < row; col++)
2109cdf0e10cSrcweir {
2110cdf0e10cSrcweir if (bIsTransposed)
2111cdf0e10cSrcweir fSum -= pMatA->GetDouble(col,row) * pMatT->GetDouble(col);
2112cdf0e10cSrcweir else
2113cdf0e10cSrcweir fSum -= pMatA->GetDouble(row,col) * pMatT->GetDouble(col);
2114cdf0e10cSrcweir }
2115cdf0e10cSrcweir pMatT->PutDouble( fSum / pVecR[row] , row);
2116cdf0e10cSrcweir }
2117cdf0e10cSrcweir }
2118cdf0e10cSrcweir
2119cdf0e10cSrcweir /* Calculates Z = R * B
2120cdf0e10cSrcweir * R is given in matrix A and vector VecR as obtained from the QR
2121*ebe15e47SJohn Bampton * decomposition in lcl_CalculateQRdecomposition. B and Z are column vectors
2122cdf0e10cSrcweir * given as matrix with at least index 0 to K-1; elements on index>=K are
2123cdf0e10cSrcweir * not used.
2124cdf0e10cSrcweir */
lcl_ApplyUpperRightTriangle(ScMatrixRef pMatA,::std::vector<double> & pVecR,ScMatrixRef pMatB,ScMatrixRef pMatZ,SCSIZE nK,bool bIsTransposed)2125cdf0e10cSrcweir void lcl_ApplyUpperRightTriangle(ScMatrixRef pMatA,
2126cdf0e10cSrcweir ::std::vector< double>& pVecR, ScMatrixRef pMatB,
2127cdf0e10cSrcweir ScMatrixRef pMatZ, SCSIZE nK, bool bIsTransposed)
2128cdf0e10cSrcweir {
2129cdf0e10cSrcweir // ScMatrix matrices are zero based, index access (column,row)
2130cdf0e10cSrcweir double fSum;
2131cdf0e10cSrcweir for (SCSIZE row = 0; row < nK; row++)
2132cdf0e10cSrcweir {
2133cdf0e10cSrcweir fSum = pVecR[row] * pMatB->GetDouble(row);
2134cdf0e10cSrcweir for (SCSIZE col = row+1; col < nK; col++)
2135cdf0e10cSrcweir if (bIsTransposed)
2136cdf0e10cSrcweir fSum += pMatA->GetDouble(row,col) * pMatB->GetDouble(col);
2137cdf0e10cSrcweir else
2138cdf0e10cSrcweir fSum += pMatA->GetDouble(col,row) * pMatB->GetDouble(col);
2139cdf0e10cSrcweir pMatZ->PutDouble( fSum, row);
2140cdf0e10cSrcweir }
2141cdf0e10cSrcweir }
2142cdf0e10cSrcweir
2143cdf0e10cSrcweir
2144cdf0e10cSrcweir
lcl_GetMeanOverAll(ScMatrixRef pMat,SCSIZE nN)2145cdf0e10cSrcweir double lcl_GetMeanOverAll(ScMatrixRef pMat, SCSIZE nN)
2146cdf0e10cSrcweir {
2147cdf0e10cSrcweir double fSum = 0.0;
2148cdf0e10cSrcweir for (SCSIZE i=0 ; i<nN; i++)
2149cdf0e10cSrcweir fSum += pMat->GetDouble(i);
2150cdf0e10cSrcweir return fSum/static_cast<double>(nN);
2151cdf0e10cSrcweir }
2152cdf0e10cSrcweir
2153cdf0e10cSrcweir // Calculates means of the columns of matrix X. X is a RxC matrix;
2154cdf0e10cSrcweir // ResMat is a 1xC matrix (=row).
lcl_CalculateColumnMeans(ScMatrixRef pX,ScMatrixRef pResMat,SCSIZE nC,SCSIZE nR)2155cdf0e10cSrcweir void lcl_CalculateColumnMeans(ScMatrixRef pX, ScMatrixRef pResMat,
2156cdf0e10cSrcweir SCSIZE nC, SCSIZE nR)
2157cdf0e10cSrcweir {
2158cdf0e10cSrcweir double fSum = 0.0;
2159cdf0e10cSrcweir for (SCSIZE i=0; i < nC; i++)
2160cdf0e10cSrcweir {
2161cdf0e10cSrcweir fSum =0.0;
2162cdf0e10cSrcweir for (SCSIZE k=0; k < nR; k++)
2163cdf0e10cSrcweir fSum += pX->GetDouble(i,k); // GetDouble(Column,Row)
2164cdf0e10cSrcweir pResMat ->PutDouble( fSum/static_cast<double>(nR),i);
2165cdf0e10cSrcweir }
2166cdf0e10cSrcweir }
2167cdf0e10cSrcweir
2168cdf0e10cSrcweir // Calculates means of the rows of matrix X. X is a RxC matrix;
2169cdf0e10cSrcweir // ResMat is a Rx1 matrix (=column).
lcl_CalculateRowMeans(ScMatrixRef pX,ScMatrixRef pResMat,SCSIZE nC,SCSIZE nR)2170cdf0e10cSrcweir void lcl_CalculateRowMeans(ScMatrixRef pX, ScMatrixRef pResMat,
2171cdf0e10cSrcweir SCSIZE nC, SCSIZE nR)
2172cdf0e10cSrcweir {
2173cdf0e10cSrcweir double fSum = 0.0;
2174cdf0e10cSrcweir for (SCSIZE k=0; k < nR; k++)
2175cdf0e10cSrcweir {
2176cdf0e10cSrcweir fSum =0.0;
2177cdf0e10cSrcweir for (SCSIZE i=0; i < nC; i++)
2178cdf0e10cSrcweir fSum += pX->GetDouble(i,k); // GetDouble(Column,Row)
2179cdf0e10cSrcweir pResMat ->PutDouble( fSum/static_cast<double>(nC),k);
2180cdf0e10cSrcweir }
2181cdf0e10cSrcweir }
2182cdf0e10cSrcweir
lcl_CalculateColumnsDelta(ScMatrixRef pMat,ScMatrixRef pColumnMeans,SCSIZE nC,SCSIZE nR)2183cdf0e10cSrcweir void lcl_CalculateColumnsDelta(ScMatrixRef pMat, ScMatrixRef pColumnMeans,
2184cdf0e10cSrcweir SCSIZE nC, SCSIZE nR)
2185cdf0e10cSrcweir {
2186cdf0e10cSrcweir for (SCSIZE i = 0; i < nC; i++)
2187cdf0e10cSrcweir for (SCSIZE k = 0; k < nR; k++)
2188cdf0e10cSrcweir pMat->PutDouble( ::rtl::math::approxSub
2189cdf0e10cSrcweir (pMat->GetDouble(i,k) , pColumnMeans->GetDouble(i) ) , i, k);
2190cdf0e10cSrcweir }
2191cdf0e10cSrcweir
lcl_CalculateRowsDelta(ScMatrixRef pMat,ScMatrixRef pRowMeans,SCSIZE nC,SCSIZE nR)2192cdf0e10cSrcweir void lcl_CalculateRowsDelta(ScMatrixRef pMat, ScMatrixRef pRowMeans,
2193cdf0e10cSrcweir SCSIZE nC, SCSIZE nR)
2194cdf0e10cSrcweir {
2195cdf0e10cSrcweir for (SCSIZE k = 0; k < nR; k++)
2196cdf0e10cSrcweir for (SCSIZE i = 0; i < nC; i++)
2197cdf0e10cSrcweir pMat->PutDouble( ::rtl::math::approxSub
2198cdf0e10cSrcweir ( pMat->GetDouble(i,k) , pRowMeans->GetDouble(k) ) , i, k);
2199cdf0e10cSrcweir }
2200cdf0e10cSrcweir
2201cdf0e10cSrcweir // Case1 = simple regression
2202cdf0e10cSrcweir // MatX = X - MeanX, MatY = Y - MeanY, y - haty = (y - MeanY) - (haty - MeanY)
2203cdf0e10cSrcweir // = (y-MeanY)-((slope*x+a)-(slope*MeanX+a)) = (y-MeanY)-slope*(x-MeanX)
lcl_GetSSresid(ScMatrixRef pMatX,ScMatrixRef pMatY,double fSlope,SCSIZE nN)2204cdf0e10cSrcweir double lcl_GetSSresid(ScMatrixRef pMatX, ScMatrixRef pMatY, double fSlope,
2205cdf0e10cSrcweir SCSIZE nN)
2206cdf0e10cSrcweir {
2207cdf0e10cSrcweir double fSum = 0.0;
2208cdf0e10cSrcweir double fTemp = 0.0;
2209cdf0e10cSrcweir for (SCSIZE i=0; i<nN; i++)
2210cdf0e10cSrcweir {
2211cdf0e10cSrcweir fTemp = pMatY->GetDouble(i) - fSlope * pMatX->GetDouble(i);
2212cdf0e10cSrcweir fSum += fTemp * fTemp;
2213cdf0e10cSrcweir }
2214cdf0e10cSrcweir return fSum;
2215cdf0e10cSrcweir }
2216cdf0e10cSrcweir
2217cdf0e10cSrcweir // Fill default values in matrix X, transform Y to log(Y) in case LOGEST|GROWTH,
2218cdf0e10cSrcweir // determine sizes of matrices X and Y, determine kind of regression, clone
2219cdf0e10cSrcweir // Y in case LOGEST|GROWTH, if constant.
CheckMatrix(bool _bLOG,sal_uInt8 & nCase,SCSIZE & nCX,SCSIZE & nCY,SCSIZE & nRX,SCSIZE & nRY,SCSIZE & M,SCSIZE & N,ScMatrixRef & pMatX,ScMatrixRef & pMatY)2220cdf0e10cSrcweir bool ScInterpreter::CheckMatrix(bool _bLOG, sal_uInt8& nCase, SCSIZE& nCX,
2221cdf0e10cSrcweir SCSIZE& nCY, SCSIZE& nRX, SCSIZE& nRY, SCSIZE& M,
2222cdf0e10cSrcweir SCSIZE& N, ScMatrixRef& pMatX, ScMatrixRef& pMatY)
2223cdf0e10cSrcweir {
2224cdf0e10cSrcweir
2225cdf0e10cSrcweir nCX = 0;
2226cdf0e10cSrcweir nCY = 0;
2227cdf0e10cSrcweir nRX = 0;
2228cdf0e10cSrcweir nRY = 0;
2229cdf0e10cSrcweir M = 0;
2230cdf0e10cSrcweir N = 0;
2231cdf0e10cSrcweir pMatY->GetDimensions(nCY, nRY);
2232cdf0e10cSrcweir const SCSIZE nCountY = nCY * nRY;
2233cdf0e10cSrcweir for ( SCSIZE i = 0; i < nCountY; i++ )
2234cdf0e10cSrcweir {
2235cdf0e10cSrcweir if (!pMatY->IsValue(i))
2236cdf0e10cSrcweir {
2237cdf0e10cSrcweir PushIllegalArgument();
2238cdf0e10cSrcweir return false;
2239cdf0e10cSrcweir }
2240cdf0e10cSrcweir }
2241cdf0e10cSrcweir
2242cdf0e10cSrcweir if ( _bLOG )
2243cdf0e10cSrcweir {
2244cdf0e10cSrcweir ScMatrixRef pNewY = pMatY->CloneIfConst();
2245cdf0e10cSrcweir for (SCSIZE nElem = 0; nElem < nCountY; nElem++)
2246cdf0e10cSrcweir {
2247cdf0e10cSrcweir const double fVal = pNewY->GetDouble(nElem);
2248cdf0e10cSrcweir if (fVal <= 0.0)
2249cdf0e10cSrcweir {
2250cdf0e10cSrcweir PushIllegalArgument();
2251cdf0e10cSrcweir return false;
2252cdf0e10cSrcweir }
2253cdf0e10cSrcweir else
2254cdf0e10cSrcweir pNewY->PutDouble(log(fVal), nElem);
2255cdf0e10cSrcweir }
2256cdf0e10cSrcweir pMatY = pNewY;
2257cdf0e10cSrcweir }
2258cdf0e10cSrcweir
2259cdf0e10cSrcweir if (pMatX)
2260cdf0e10cSrcweir {
2261cdf0e10cSrcweir pMatX->GetDimensions(nCX, nRX);
2262cdf0e10cSrcweir const SCSIZE nCountX = nCX * nRX;
2263cdf0e10cSrcweir for ( SCSIZE i = 0; i < nCountX; i++ )
2264cdf0e10cSrcweir if (!pMatX->IsValue(i))
2265cdf0e10cSrcweir {
2266cdf0e10cSrcweir PushIllegalArgument();
2267cdf0e10cSrcweir return false;
2268cdf0e10cSrcweir }
2269cdf0e10cSrcweir if (nCX == nCY && nRX == nRY)
2270cdf0e10cSrcweir {
2271cdf0e10cSrcweir nCase = 1; // simple regression
2272cdf0e10cSrcweir M = 1;
2273cdf0e10cSrcweir N = nCountY;
2274cdf0e10cSrcweir }
2275cdf0e10cSrcweir else if (nCY != 1 && nRY != 1)
2276cdf0e10cSrcweir {
2277cdf0e10cSrcweir PushIllegalArgument();
2278cdf0e10cSrcweir return false;
2279cdf0e10cSrcweir }
2280cdf0e10cSrcweir else if (nCY == 1)
2281cdf0e10cSrcweir {
2282cdf0e10cSrcweir if (nRX != nRY)
2283cdf0e10cSrcweir {
2284cdf0e10cSrcweir PushIllegalArgument();
2285cdf0e10cSrcweir return false;
2286cdf0e10cSrcweir }
2287cdf0e10cSrcweir else
2288cdf0e10cSrcweir {
2289cdf0e10cSrcweir nCase = 2; // Y is column
2290cdf0e10cSrcweir N = nRY;
2291cdf0e10cSrcweir M = nCX;
2292cdf0e10cSrcweir }
2293cdf0e10cSrcweir }
2294cdf0e10cSrcweir else if (nCX != nCY)
2295cdf0e10cSrcweir {
2296cdf0e10cSrcweir PushIllegalArgument();
2297cdf0e10cSrcweir return false;
2298cdf0e10cSrcweir }
2299cdf0e10cSrcweir else
2300cdf0e10cSrcweir {
2301cdf0e10cSrcweir nCase = 3; // Y is row
2302cdf0e10cSrcweir N = nCY;
2303cdf0e10cSrcweir M = nRX;
2304cdf0e10cSrcweir }
2305cdf0e10cSrcweir }
2306cdf0e10cSrcweir else
2307cdf0e10cSrcweir {
2308cdf0e10cSrcweir pMatX = GetNewMat(nCY, nRY);
2309cdf0e10cSrcweir nCX = nCY;
2310cdf0e10cSrcweir nRX = nRY;
2311cdf0e10cSrcweir if (!pMatX)
2312cdf0e10cSrcweir {
2313cdf0e10cSrcweir PushIllegalArgument();
2314cdf0e10cSrcweir return false;
2315cdf0e10cSrcweir }
2316cdf0e10cSrcweir for ( SCSIZE i = 1; i <= nCountY; i++ )
2317cdf0e10cSrcweir pMatX->PutDouble(static_cast<double>(i), i-1);
2318cdf0e10cSrcweir nCase = 1;
2319cdf0e10cSrcweir N = nCountY;
2320cdf0e10cSrcweir M = 1;
2321cdf0e10cSrcweir }
2322cdf0e10cSrcweir return true;
2323cdf0e10cSrcweir }
2324cdf0e10cSrcweir // -----------------------------------------------------------------------------
2325cdf0e10cSrcweir
2326cdf0e10cSrcweir // LINEST
ScRGP()2327cdf0e10cSrcweir void ScInterpreter::ScRGP()
2328cdf0e10cSrcweir {
2329cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRGP" );
2330cdf0e10cSrcweir CalulateRGPRKP(false);
2331cdf0e10cSrcweir }
2332cdf0e10cSrcweir
2333cdf0e10cSrcweir // LOGEST
ScRKP()2334cdf0e10cSrcweir void ScInterpreter::ScRKP()
2335cdf0e10cSrcweir {
2336cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRKP" );
2337cdf0e10cSrcweir CalulateRGPRKP(true);
2338cdf0e10cSrcweir }
2339cdf0e10cSrcweir
CalulateRGPRKP(bool _bRKP)2340cdf0e10cSrcweir void ScInterpreter::CalulateRGPRKP(bool _bRKP)
2341cdf0e10cSrcweir {
2342cdf0e10cSrcweir sal_uInt8 nParamCount = GetByte();
2343cdf0e10cSrcweir if ( !MustHaveParamCount( nParamCount, 1, 4 ) )
2344cdf0e10cSrcweir return;
2345cdf0e10cSrcweir bool bConstant, bStats;
2346cdf0e10cSrcweir
2347cdf0e10cSrcweir // optional forth parameter
2348cdf0e10cSrcweir if (nParamCount == 4)
2349cdf0e10cSrcweir bStats = GetBool();
2350cdf0e10cSrcweir else
2351cdf0e10cSrcweir bStats = false;
2352cdf0e10cSrcweir
2353cdf0e10cSrcweir // The third parameter may not be missing in ODF, if the forth parameter
2354cdf0e10cSrcweir // is present. But Excel allows it with default true, we too.
2355cdf0e10cSrcweir if (nParamCount >= 3)
2356cdf0e10cSrcweir {
2357cdf0e10cSrcweir if (IsMissing())
2358cdf0e10cSrcweir {
2359cdf0e10cSrcweir Pop();
2360cdf0e10cSrcweir bConstant = true;
2361cdf0e10cSrcweir // PushIllegalParameter(); if ODF behavior is desired
2362cdf0e10cSrcweir // return;
2363cdf0e10cSrcweir }
2364cdf0e10cSrcweir else
2365cdf0e10cSrcweir bConstant = GetBool();
2366cdf0e10cSrcweir }
2367cdf0e10cSrcweir else
2368cdf0e10cSrcweir bConstant = true;
2369cdf0e10cSrcweir
2370cdf0e10cSrcweir ScMatrixRef pMatX;
2371cdf0e10cSrcweir if (nParamCount >= 2)
2372cdf0e10cSrcweir {
2373cdf0e10cSrcweir if (IsMissing())
2374cdf0e10cSrcweir {
2375cdf0e10cSrcweir // In ODF1.2 empty second parameter (which is two ;; ) is allowed
2376cdf0e10cSrcweir Pop();
2377cdf0e10cSrcweir pMatX = NULL;
2378cdf0e10cSrcweir }
2379cdf0e10cSrcweir else
2380cdf0e10cSrcweir {
2381cdf0e10cSrcweir pMatX = GetMatrix();
2382cdf0e10cSrcweir }
2383cdf0e10cSrcweir }
2384cdf0e10cSrcweir else
2385cdf0e10cSrcweir pMatX = NULL;
2386cdf0e10cSrcweir
2387cdf0e10cSrcweir ScMatrixRef pMatY;
2388cdf0e10cSrcweir pMatY = GetMatrix();
2389cdf0e10cSrcweir if (!pMatY)
2390cdf0e10cSrcweir {
2391cdf0e10cSrcweir PushIllegalParameter();
2392cdf0e10cSrcweir return;
2393cdf0e10cSrcweir }
2394cdf0e10cSrcweir
2395cdf0e10cSrcweir // 1 = simple; 2 = multiple with Y as column; 3 = multiple with Y as row
2396cdf0e10cSrcweir sal_uInt8 nCase;
2397cdf0e10cSrcweir
2398cdf0e10cSrcweir SCSIZE nCX, nCY; // number of columns
2399cdf0e10cSrcweir SCSIZE nRX, nRY; // number of rows
2400cdf0e10cSrcweir SCSIZE K = 0, N = 0; // K=number of variables X, N=number of data samples
2401cdf0e10cSrcweir if ( !CheckMatrix(_bRKP,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY) )
2402cdf0e10cSrcweir {
2403cdf0e10cSrcweir PushIllegalParameter();
2404cdf0e10cSrcweir return;
2405cdf0e10cSrcweir }
2406cdf0e10cSrcweir
2407cdf0e10cSrcweir // Enough data samples?
2408cdf0e10cSrcweir if ( (bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1) )
2409cdf0e10cSrcweir {
2410cdf0e10cSrcweir PushIllegalParameter();
2411cdf0e10cSrcweir return;
2412cdf0e10cSrcweir }
2413cdf0e10cSrcweir
2414cdf0e10cSrcweir ScMatrixRef pResMat;
2415cdf0e10cSrcweir if (bStats)
2416cdf0e10cSrcweir pResMat = GetNewMat(K+1,5);
2417cdf0e10cSrcweir else
2418cdf0e10cSrcweir pResMat = GetNewMat(K+1,1);
2419cdf0e10cSrcweir if (!pResMat)
2420cdf0e10cSrcweir {
2421cdf0e10cSrcweir PushError(errCodeOverflow);
2422cdf0e10cSrcweir return;
2423cdf0e10cSrcweir }
2424cdf0e10cSrcweir // Fill unused cells in pResMat; order (column,row)
2425cdf0e10cSrcweir if (bStats)
2426cdf0e10cSrcweir {
2427cdf0e10cSrcweir for (SCSIZE i=2; i<K+1; i++)
2428cdf0e10cSrcweir {
2429cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 2 );
2430cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 3 );
2431cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 4 );
2432cdf0e10cSrcweir }
2433cdf0e10cSrcweir }
2434cdf0e10cSrcweir
2435cdf0e10cSrcweir // Uses sum(x-MeanX)^2 and not [sum x^2]-N * MeanX^2 in case bConstant.
2436cdf0e10cSrcweir // Clone constant matrices, so that Mat = Mat - Mean is possible.
2437cdf0e10cSrcweir double fMeanY = 0.0;
2438cdf0e10cSrcweir if (bConstant)
2439cdf0e10cSrcweir {
2440cdf0e10cSrcweir ScMatrixRef pNewX = pMatX->CloneIfConst();
2441cdf0e10cSrcweir ScMatrixRef pNewY = pMatY->CloneIfConst();
2442cdf0e10cSrcweir if (!pNewX || !pNewY)
2443cdf0e10cSrcweir {
2444cdf0e10cSrcweir PushError(errCodeOverflow);
2445cdf0e10cSrcweir return;
2446cdf0e10cSrcweir }
2447cdf0e10cSrcweir pMatX = pNewX;
2448cdf0e10cSrcweir pMatY = pNewY;
2449cdf0e10cSrcweir // DeltaY is possible here; DeltaX depends on nCase, so later
2450cdf0e10cSrcweir fMeanY = lcl_GetMeanOverAll(pMatY, N);
2451cdf0e10cSrcweir for (SCSIZE i=0; i<N; i++)
2452cdf0e10cSrcweir {
2453cdf0e10cSrcweir pMatY->PutDouble( ::rtl::math::approxSub(pMatY->GetDouble(i),fMeanY), i );
2454cdf0e10cSrcweir }
2455cdf0e10cSrcweir }
2456cdf0e10cSrcweir
2457cdf0e10cSrcweir if (nCase==1)
2458cdf0e10cSrcweir {
2459cdf0e10cSrcweir // calculate simple regression
2460cdf0e10cSrcweir double fMeanX = 0.0;
2461cdf0e10cSrcweir if (bConstant)
2462cdf0e10cSrcweir { // Mat = Mat - Mean
2463cdf0e10cSrcweir fMeanX = lcl_GetMeanOverAll(pMatX, N);
2464cdf0e10cSrcweir for (SCSIZE i=0; i<N; i++)
2465cdf0e10cSrcweir {
2466cdf0e10cSrcweir pMatX->PutDouble( ::rtl::math::approxSub(pMatX->GetDouble(i),fMeanX), i );
2467cdf0e10cSrcweir }
2468cdf0e10cSrcweir }
2469cdf0e10cSrcweir double fSumXY = lcl_GetSumProduct(pMatX,pMatY,N);
2470cdf0e10cSrcweir double fSumX2 = lcl_GetSumProduct(pMatX,pMatX,N);
2471cdf0e10cSrcweir if (fSumX2==0.0)
2472cdf0e10cSrcweir {
2473cdf0e10cSrcweir PushNoValue(); // all x-values are identical
2474cdf0e10cSrcweir return;
2475cdf0e10cSrcweir }
2476cdf0e10cSrcweir double fSlope = fSumXY / fSumX2;
2477cdf0e10cSrcweir double fIntercept = 0.0;
2478cdf0e10cSrcweir if (bConstant)
2479cdf0e10cSrcweir fIntercept = fMeanY - fSlope * fMeanX;
2480cdf0e10cSrcweir pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, 1, 0); //order (column,row)
2481cdf0e10cSrcweir pResMat->PutDouble(_bRKP ? exp(fSlope) : fSlope, 0, 0);
2482cdf0e10cSrcweir
2483cdf0e10cSrcweir if (bStats)
2484cdf0e10cSrcweir {
2485cdf0e10cSrcweir double fSSreg = fSlope * fSlope * fSumX2;
2486cdf0e10cSrcweir pResMat->PutDouble(fSSreg, 0, 4);
2487cdf0e10cSrcweir
2488cdf0e10cSrcweir double fDegreesFreedom =static_cast<double>( (bConstant) ? N-2 : N-1 );
2489cdf0e10cSrcweir pResMat->PutDouble(fDegreesFreedom, 1, 3);
2490cdf0e10cSrcweir
2491cdf0e10cSrcweir double fSSresid = lcl_GetSSresid(pMatX,pMatY,fSlope,N);
2492cdf0e10cSrcweir pResMat->PutDouble(fSSresid, 1, 4);
2493cdf0e10cSrcweir
2494cdf0e10cSrcweir if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0)
2495cdf0e10cSrcweir { // exact fit; test SSreg too, because SSresid might be
2496cdf0e10cSrcweir // unequal zero due to round of errors
2497cdf0e10cSrcweir pResMat->PutDouble(0.0, 1, 4); // SSresid
2498cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 0, 3); // F
2499cdf0e10cSrcweir pResMat->PutDouble(0.0, 1, 2); // RMSE
2500cdf0e10cSrcweir pResMat->PutDouble(0.0, 0, 1); // SigmaSlope
2501cdf0e10cSrcweir if (bConstant)
2502cdf0e10cSrcweir pResMat->PutDouble(0.0, 1, 1); //SigmaIntercept
2503cdf0e10cSrcweir else
2504cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 1, 1);
2505cdf0e10cSrcweir pResMat->PutDouble(1.0, 0, 2); // R^2
2506cdf0e10cSrcweir }
2507cdf0e10cSrcweir else
2508cdf0e10cSrcweir {
2509cdf0e10cSrcweir double fFstatistic = (fSSreg / static_cast<double>(K))
2510cdf0e10cSrcweir / (fSSresid / fDegreesFreedom);
2511cdf0e10cSrcweir pResMat->PutDouble(fFstatistic, 0, 3);
2512cdf0e10cSrcweir
2513cdf0e10cSrcweir // standard error of estimate
2514cdf0e10cSrcweir double fRMSE = sqrt(fSSresid / fDegreesFreedom);
2515cdf0e10cSrcweir pResMat->PutDouble(fRMSE, 1, 2);
2516cdf0e10cSrcweir
2517cdf0e10cSrcweir double fSigmaSlope = fRMSE / sqrt(fSumX2);
2518cdf0e10cSrcweir pResMat->PutDouble(fSigmaSlope, 0, 1);
2519cdf0e10cSrcweir
2520cdf0e10cSrcweir if (bConstant)
2521cdf0e10cSrcweir {
2522cdf0e10cSrcweir double fSigmaIntercept = fRMSE
2523cdf0e10cSrcweir * sqrt(fMeanX*fMeanX/fSumX2 + 1.0/static_cast<double>(N));
2524cdf0e10cSrcweir pResMat->PutDouble(fSigmaIntercept, 1, 1);
2525cdf0e10cSrcweir }
2526cdf0e10cSrcweir else
2527cdf0e10cSrcweir {
2528cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 1, 1);
2529cdf0e10cSrcweir }
2530cdf0e10cSrcweir
2531cdf0e10cSrcweir double fR2 = fSSreg / (fSSreg + fSSresid);
2532cdf0e10cSrcweir pResMat->PutDouble(fR2, 0, 2);
2533cdf0e10cSrcweir }
2534cdf0e10cSrcweir }
2535cdf0e10cSrcweir PushMatrix(pResMat);
2536cdf0e10cSrcweir }
2537cdf0e10cSrcweir else // calculate multiple regression;
2538cdf0e10cSrcweir {
2539cdf0e10cSrcweir // Uses a QR decomposition X = QR. The solution B = (X'X)^(-1) * X' * Y
2540cdf0e10cSrcweir // becomes B = R^(-1) * Q' * Y
2541cdf0e10cSrcweir if (nCase ==2) // Y is column
2542cdf0e10cSrcweir {
2543cdf0e10cSrcweir ::std::vector< double> aVecR(N); // for QR decomposition
2544cdf0e10cSrcweir // Enough memory for needed matrices?
2545cdf0e10cSrcweir ScMatrixRef pMeans = GetNewMat(K, 1); // mean of each column
2546cdf0e10cSrcweir ScMatrixRef pMatZ; // for Q' * Y , inter alia
2547cdf0e10cSrcweir if (bStats)
2548cdf0e10cSrcweir pMatZ = pMatY->Clone(); // Y is used in statistic, keep it
2549cdf0e10cSrcweir else
2550cdf0e10cSrcweir pMatZ = pMatY; // Y can be overwritten
2551cdf0e10cSrcweir ScMatrixRef pSlopes = GetNewMat(1,K); // from b1 to bK
2552cdf0e10cSrcweir if (!pMeans || !pMatZ || !pSlopes)
2553cdf0e10cSrcweir {
2554cdf0e10cSrcweir PushError(errCodeOverflow);
2555cdf0e10cSrcweir return;
2556cdf0e10cSrcweir }
2557cdf0e10cSrcweir if (bConstant)
2558cdf0e10cSrcweir {
2559cdf0e10cSrcweir lcl_CalculateColumnMeans(pMatX, pMeans, K, N);
2560cdf0e10cSrcweir lcl_CalculateColumnsDelta(pMatX, pMeans, K, N);
2561cdf0e10cSrcweir }
2562cdf0e10cSrcweir if (!lcl_CalculateQRdecomposition(pMatX, aVecR, K, N))
2563cdf0e10cSrcweir {
2564cdf0e10cSrcweir PushNoValue();
2565cdf0e10cSrcweir return;
2566cdf0e10cSrcweir }
2567cdf0e10cSrcweir // Later on we will divide by elements of aVecR, so make sure
2568cdf0e10cSrcweir // that they aren't zero.
2569cdf0e10cSrcweir bool bIsSingular=false;
2570cdf0e10cSrcweir for (SCSIZE row=0; row < K && !bIsSingular; row++)
2571cdf0e10cSrcweir bIsSingular = bIsSingular || aVecR[row]==0.0;
2572cdf0e10cSrcweir if (bIsSingular)
2573cdf0e10cSrcweir {
2574cdf0e10cSrcweir PushNoValue();
2575cdf0e10cSrcweir return;
2576cdf0e10cSrcweir }
2577cdf0e10cSrcweir // Z = Q' Y;
2578cdf0e10cSrcweir for (SCSIZE col = 0; col < K; col++)
2579cdf0e10cSrcweir {
2580cdf0e10cSrcweir lcl_ApplyHouseholderTransformation(pMatX, col, pMatZ, N);
2581cdf0e10cSrcweir }
2582cdf0e10cSrcweir // B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
2583cdf0e10cSrcweir // result Z should have zeros for index>=K; if not, ignore values
2584cdf0e10cSrcweir for (SCSIZE col = 0; col < K ; col++)
2585cdf0e10cSrcweir {
2586cdf0e10cSrcweir pSlopes->PutDouble( pMatZ->GetDouble(col), col);
2587cdf0e10cSrcweir }
2588cdf0e10cSrcweir lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, false);
2589cdf0e10cSrcweir double fIntercept = 0.0;
2590cdf0e10cSrcweir if (bConstant)
2591cdf0e10cSrcweir fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K);
2592cdf0e10cSrcweir // Fill first line in result matrix
2593cdf0e10cSrcweir pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, K, 0 );
2594cdf0e10cSrcweir for (SCSIZE i = 0; i < K; i++)
2595cdf0e10cSrcweir pResMat->PutDouble(_bRKP ? exp(pSlopes->GetDouble(i))
2596cdf0e10cSrcweir : pSlopes->GetDouble(i) , K-1-i, 0);
2597cdf0e10cSrcweir
2598cdf0e10cSrcweir
2599cdf0e10cSrcweir if (bStats)
2600cdf0e10cSrcweir {
2601cdf0e10cSrcweir double fSSreg = 0.0;
2602cdf0e10cSrcweir double fSSresid = 0.0;
2603cdf0e10cSrcweir // re-use memory of Z;
2604cdf0e10cSrcweir pMatZ->FillDouble(0.0, 0, 0, 0, N-1);
2605cdf0e10cSrcweir // Z = R * Slopes
2606cdf0e10cSrcweir lcl_ApplyUpperRightTriangle(pMatX, aVecR, pSlopes, pMatZ, K, false);
2607cdf0e10cSrcweir // Z = Q * Z, that is Q * R * Slopes = X * Slopes
2608cdf0e10cSrcweir for (SCSIZE colp1 = K; colp1 > 0; colp1--)
2609cdf0e10cSrcweir {
2610cdf0e10cSrcweir lcl_ApplyHouseholderTransformation(pMatX, colp1-1, pMatZ,N);
2611cdf0e10cSrcweir }
2612cdf0e10cSrcweir fSSreg =lcl_GetSumProduct(pMatZ, pMatZ, N);
2613cdf0e10cSrcweir // re-use Y for residuals, Y = Y-Z
2614cdf0e10cSrcweir for (SCSIZE row = 0; row < N; row++)
2615cdf0e10cSrcweir pMatY->PutDouble(pMatY->GetDouble(row) - pMatZ->GetDouble(row), row);
2616cdf0e10cSrcweir fSSresid = lcl_GetSumProduct(pMatY, pMatY, N);
2617cdf0e10cSrcweir pResMat->PutDouble(fSSreg, 0, 4);
2618cdf0e10cSrcweir pResMat->PutDouble(fSSresid, 1, 4);
2619cdf0e10cSrcweir
2620cdf0e10cSrcweir double fDegreesFreedom =static_cast<double>( (bConstant) ? N-K-1 : N-K );
2621cdf0e10cSrcweir pResMat->PutDouble(fDegreesFreedom, 1, 3);
2622cdf0e10cSrcweir
2623cdf0e10cSrcweir if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0)
2624cdf0e10cSrcweir { // exact fit; incl. observed values Y are identical
2625cdf0e10cSrcweir pResMat->PutDouble(0.0, 1, 4); // SSresid
2626cdf0e10cSrcweir // F = (SSreg/K) / (SSresid/df) = #DIV/0!
2627cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 0, 3); // F
2628cdf0e10cSrcweir // RMSE = sqrt(SSresid / df) = sqrt(0 / df) = 0
2629cdf0e10cSrcweir pResMat->PutDouble(0.0, 1, 2); // RMSE
2630cdf0e10cSrcweir // SigmaSlope[i] = RMSE * sqrt(matrix[i,i]) = 0 * sqrt(...) = 0
2631cdf0e10cSrcweir for (SCSIZE i=0; i<K; i++)
2632cdf0e10cSrcweir pResMat->PutDouble(0.0, K-1-i, 1);
2633cdf0e10cSrcweir
2634cdf0e10cSrcweir // SigmaIntercept = RMSE * sqrt(...) = 0
2635cdf0e10cSrcweir if (bConstant)
2636cdf0e10cSrcweir pResMat->PutDouble(0.0, K, 1); //SigmaIntercept
2637cdf0e10cSrcweir else
2638cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
2639cdf0e10cSrcweir
2640cdf0e10cSrcweir // R^2 = SSreg / (SSreg + SSresid) = 1.0
2641cdf0e10cSrcweir pResMat->PutDouble(1.0, 0, 2); // R^2
2642cdf0e10cSrcweir }
2643cdf0e10cSrcweir else
2644cdf0e10cSrcweir {
2645cdf0e10cSrcweir double fFstatistic = (fSSreg / static_cast<double>(K))
2646cdf0e10cSrcweir / (fSSresid / fDegreesFreedom);
2647cdf0e10cSrcweir pResMat->PutDouble(fFstatistic, 0, 3);
2648cdf0e10cSrcweir
2649cdf0e10cSrcweir // standard error of estimate = root mean SSE
2650cdf0e10cSrcweir double fRMSE = sqrt(fSSresid / fDegreesFreedom);
2651cdf0e10cSrcweir pResMat->PutDouble(fRMSE, 1, 2);
2652cdf0e10cSrcweir
2653cdf0e10cSrcweir // standard error of slopes
2654cdf0e10cSrcweir // = RMSE * sqrt(diagonal element of (R' R)^(-1) )
2655cdf0e10cSrcweir // standard error of intercept
2656cdf0e10cSrcweir // = RMSE * sqrt( Xmean * (R' R)^(-1) * Xmean' + 1/N)
2657cdf0e10cSrcweir // (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as
2658cdf0e10cSrcweir // a whole matrix, but iterate over unit vectors.
2659cdf0e10cSrcweir double fSigmaSlope = 0.0;
2660cdf0e10cSrcweir double fSigmaIntercept = 0.0;
2661cdf0e10cSrcweir double fPart; // for Xmean * single column of (R' R)^(-1)
2662cdf0e10cSrcweir for (SCSIZE col = 0; col < K; col++)
2663cdf0e10cSrcweir {
2664cdf0e10cSrcweir //re-use memory of MatZ
2665cdf0e10cSrcweir pMatZ->FillDouble(0.0,0,0,0,K-1); // Z = unit vector e
2666cdf0e10cSrcweir pMatZ->PutDouble(1.0, col);
2667cdf0e10cSrcweir //Solve R' * Z = e
2668cdf0e10cSrcweir lcl_SolveWithLowerLeftTriangle(pMatX, aVecR, pMatZ, K, false);
2669cdf0e10cSrcweir // Solve R * Znew = Zold
2670cdf0e10cSrcweir lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pMatZ, K, false);
2671cdf0e10cSrcweir // now Z is column col in (R' R)^(-1)
2672cdf0e10cSrcweir fSigmaSlope = fRMSE * sqrt(pMatZ->GetDouble(col));
2673cdf0e10cSrcweir pResMat->PutDouble(fSigmaSlope, K-1-col, 1);
2674cdf0e10cSrcweir // (R' R) ^(-1) is symmetric
2675cdf0e10cSrcweir if (bConstant)
2676cdf0e10cSrcweir {
2677cdf0e10cSrcweir fPart = lcl_GetSumProduct(pMeans, pMatZ, K);
2678cdf0e10cSrcweir fSigmaIntercept += fPart * pMeans->GetDouble(col);
2679cdf0e10cSrcweir }
2680cdf0e10cSrcweir }
2681cdf0e10cSrcweir if (bConstant)
2682cdf0e10cSrcweir {
2683cdf0e10cSrcweir fSigmaIntercept = fRMSE
2684cdf0e10cSrcweir * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N));
2685cdf0e10cSrcweir pResMat->PutDouble(fSigmaIntercept, K, 1);
2686cdf0e10cSrcweir }
2687cdf0e10cSrcweir else
2688cdf0e10cSrcweir {
2689cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
2690cdf0e10cSrcweir }
2691cdf0e10cSrcweir
2692cdf0e10cSrcweir double fR2 = fSSreg / (fSSreg + fSSresid);
2693cdf0e10cSrcweir pResMat->PutDouble(fR2, 0, 2);
2694cdf0e10cSrcweir }
2695cdf0e10cSrcweir }
2696cdf0e10cSrcweir PushMatrix(pResMat);
2697cdf0e10cSrcweir }
2698cdf0e10cSrcweir else // nCase == 3, Y is row, all matrices are transposed
2699cdf0e10cSrcweir {
2700cdf0e10cSrcweir ::std::vector< double> aVecR(N); // for QR decomposition
2701cdf0e10cSrcweir // Enough memory for needed matrices?
2702cdf0e10cSrcweir ScMatrixRef pMeans = GetNewMat(1, K); // mean of each row
2703cdf0e10cSrcweir ScMatrixRef pMatZ; // for Q' * Y , inter alia
2704cdf0e10cSrcweir if (bStats)
2705cdf0e10cSrcweir pMatZ = pMatY->Clone(); // Y is used in statistic, keep it
2706cdf0e10cSrcweir else
2707cdf0e10cSrcweir pMatZ = pMatY; // Y can be overwritten
2708cdf0e10cSrcweir ScMatrixRef pSlopes = GetNewMat(K,1); // from b1 to bK
2709cdf0e10cSrcweir if (!pMeans || !pMatZ || !pSlopes)
2710cdf0e10cSrcweir {
2711cdf0e10cSrcweir PushError(errCodeOverflow);
2712cdf0e10cSrcweir return;
2713cdf0e10cSrcweir }
2714cdf0e10cSrcweir if (bConstant)
2715cdf0e10cSrcweir {
2716cdf0e10cSrcweir lcl_CalculateRowMeans(pMatX, pMeans, N, K);
2717cdf0e10cSrcweir lcl_CalculateRowsDelta(pMatX, pMeans, N, K);
2718cdf0e10cSrcweir }
2719cdf0e10cSrcweir
2720cdf0e10cSrcweir if (!lcl_TCalculateQRdecomposition(pMatX, aVecR, K, N))
2721cdf0e10cSrcweir {
2722cdf0e10cSrcweir PushNoValue();
2723cdf0e10cSrcweir return;
2724cdf0e10cSrcweir }
2725cdf0e10cSrcweir
2726cdf0e10cSrcweir // Later on we will divide by elements of aVecR, so make sure
2727cdf0e10cSrcweir // that they aren't zero.
2728cdf0e10cSrcweir bool bIsSingular=false;
2729cdf0e10cSrcweir for (SCSIZE row=0; row < K && !bIsSingular; row++)
2730cdf0e10cSrcweir bIsSingular = bIsSingular || aVecR[row]==0.0;
2731cdf0e10cSrcweir if (bIsSingular)
2732cdf0e10cSrcweir {
2733cdf0e10cSrcweir PushNoValue();
2734cdf0e10cSrcweir return;
2735cdf0e10cSrcweir }
2736cdf0e10cSrcweir // Z = Q' Y
2737cdf0e10cSrcweir for (SCSIZE row = 0; row < K; row++)
2738cdf0e10cSrcweir {
2739cdf0e10cSrcweir lcl_TApplyHouseholderTransformation(pMatX, row, pMatZ, N);
2740cdf0e10cSrcweir }
2741cdf0e10cSrcweir // B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
2742cdf0e10cSrcweir // result Z should have zeros for index>=K; if not, ignore values
2743cdf0e10cSrcweir for (SCSIZE col = 0; col < K ; col++)
2744cdf0e10cSrcweir {
2745cdf0e10cSrcweir pSlopes->PutDouble( pMatZ->GetDouble(col), col);
2746cdf0e10cSrcweir }
2747cdf0e10cSrcweir lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, true);
2748cdf0e10cSrcweir double fIntercept = 0.0;
2749cdf0e10cSrcweir if (bConstant)
2750cdf0e10cSrcweir fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K);
2751cdf0e10cSrcweir // Fill first line in result matrix
2752cdf0e10cSrcweir pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, K, 0 );
2753cdf0e10cSrcweir for (SCSIZE i = 0; i < K; i++)
2754cdf0e10cSrcweir pResMat->PutDouble(_bRKP ? exp(pSlopes->GetDouble(i))
2755cdf0e10cSrcweir : pSlopes->GetDouble(i) , K-1-i, 0);
2756cdf0e10cSrcweir
2757cdf0e10cSrcweir
2758cdf0e10cSrcweir if (bStats)
2759cdf0e10cSrcweir {
2760cdf0e10cSrcweir double fSSreg = 0.0;
2761cdf0e10cSrcweir double fSSresid = 0.0;
2762cdf0e10cSrcweir // re-use memory of Z;
2763cdf0e10cSrcweir pMatZ->FillDouble(0.0, 0, 0, N-1, 0);
2764cdf0e10cSrcweir // Z = R * Slopes
2765cdf0e10cSrcweir lcl_ApplyUpperRightTriangle(pMatX, aVecR, pSlopes, pMatZ, K, true);
2766cdf0e10cSrcweir // Z = Q * Z, that is Q * R * Slopes = X * Slopes
2767cdf0e10cSrcweir for (SCSIZE rowp1 = K; rowp1 > 0; rowp1--)
2768cdf0e10cSrcweir {
2769cdf0e10cSrcweir lcl_TApplyHouseholderTransformation(pMatX, rowp1-1, pMatZ,N);
2770cdf0e10cSrcweir }
2771cdf0e10cSrcweir fSSreg =lcl_GetSumProduct(pMatZ, pMatZ, N);
2772cdf0e10cSrcweir // re-use Y for residuals, Y = Y-Z
2773cdf0e10cSrcweir for (SCSIZE col = 0; col < N; col++)
2774cdf0e10cSrcweir pMatY->PutDouble(pMatY->GetDouble(col) - pMatZ->GetDouble(col), col);
2775cdf0e10cSrcweir fSSresid = lcl_GetSumProduct(pMatY, pMatY, N);
2776cdf0e10cSrcweir pResMat->PutDouble(fSSreg, 0, 4);
2777cdf0e10cSrcweir pResMat->PutDouble(fSSresid, 1, 4);
2778cdf0e10cSrcweir
2779cdf0e10cSrcweir double fDegreesFreedom =static_cast<double>( (bConstant) ? N-K-1 : N-K );
2780cdf0e10cSrcweir pResMat->PutDouble(fDegreesFreedom, 1, 3);
2781cdf0e10cSrcweir
2782cdf0e10cSrcweir if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0)
2783cdf0e10cSrcweir { // exact fit; incl. case observed values Y are identical
2784cdf0e10cSrcweir pResMat->PutDouble(0.0, 1, 4); // SSresid
2785cdf0e10cSrcweir // F = (SSreg/K) / (SSresid/df) = #DIV/0!
2786cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 0, 3); // F
2787cdf0e10cSrcweir // RMSE = sqrt(SSresid / df) = sqrt(0 / df) = 0
2788cdf0e10cSrcweir pResMat->PutDouble(0.0, 1, 2); // RMSE
2789cdf0e10cSrcweir // SigmaSlope[i] = RMSE * sqrt(matrix[i,i]) = 0 * sqrt(...) = 0
2790cdf0e10cSrcweir for (SCSIZE i=0; i<K; i++)
2791cdf0e10cSrcweir pResMat->PutDouble(0.0, K-1-i, 1);
2792cdf0e10cSrcweir
2793cdf0e10cSrcweir // SigmaIntercept = RMSE * sqrt(...) = 0
2794cdf0e10cSrcweir if (bConstant)
2795cdf0e10cSrcweir pResMat->PutDouble(0.0, K, 1); //SigmaIntercept
2796cdf0e10cSrcweir else
2797cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
2798cdf0e10cSrcweir
2799cdf0e10cSrcweir // R^2 = SSreg / (SSreg + SSresid) = 1.0
2800cdf0e10cSrcweir pResMat->PutDouble(1.0, 0, 2); // R^2
2801cdf0e10cSrcweir }
2802cdf0e10cSrcweir else
2803cdf0e10cSrcweir {
2804cdf0e10cSrcweir double fFstatistic = (fSSreg / static_cast<double>(K))
2805cdf0e10cSrcweir / (fSSresid / fDegreesFreedom);
2806cdf0e10cSrcweir pResMat->PutDouble(fFstatistic, 0, 3);
2807cdf0e10cSrcweir
2808cdf0e10cSrcweir // standard error of estimate = root mean SSE
2809cdf0e10cSrcweir double fRMSE = sqrt(fSSresid / fDegreesFreedom);
2810cdf0e10cSrcweir pResMat->PutDouble(fRMSE, 1, 2);
2811cdf0e10cSrcweir
2812cdf0e10cSrcweir // standard error of slopes
2813cdf0e10cSrcweir // = RMSE * sqrt(diagonal element of (R' R)^(-1) )
2814cdf0e10cSrcweir // standard error of intercept
2815cdf0e10cSrcweir // = RMSE * sqrt( Xmean * (R' R)^(-1) * Xmean' + 1/N)
2816cdf0e10cSrcweir // (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as
2817cdf0e10cSrcweir // a whole matrix, but iterate over unit vectors.
2818cdf0e10cSrcweir // (R' R) ^(-1) is symmetric
2819cdf0e10cSrcweir double fSigmaSlope = 0.0;
2820cdf0e10cSrcweir double fSigmaIntercept = 0.0;
2821cdf0e10cSrcweir double fPart; // for Xmean * single col of (R' R)^(-1)
2822cdf0e10cSrcweir for (SCSIZE row = 0; row < K; row++)
2823cdf0e10cSrcweir {
2824cdf0e10cSrcweir //re-use memory of MatZ
2825cdf0e10cSrcweir pMatZ->FillDouble(0.0,0,0,K-1,0); // Z = unit vector e
2826cdf0e10cSrcweir pMatZ->PutDouble(1.0, row);
2827cdf0e10cSrcweir //Solve R' * Z = e
2828cdf0e10cSrcweir lcl_SolveWithLowerLeftTriangle(pMatX, aVecR, pMatZ, K, true);
2829cdf0e10cSrcweir // Solve R * Znew = Zold
2830cdf0e10cSrcweir lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pMatZ, K, true);
2831cdf0e10cSrcweir // now Z is column col in (R' R)^(-1)
2832cdf0e10cSrcweir fSigmaSlope = fRMSE * sqrt(pMatZ->GetDouble(row));
2833cdf0e10cSrcweir pResMat->PutDouble(fSigmaSlope, K-1-row, 1);
2834cdf0e10cSrcweir if (bConstant)
2835cdf0e10cSrcweir {
2836cdf0e10cSrcweir fPart = lcl_GetSumProduct(pMeans, pMatZ, K);
2837cdf0e10cSrcweir fSigmaIntercept += fPart * pMeans->GetDouble(row);
2838cdf0e10cSrcweir }
2839cdf0e10cSrcweir }
2840cdf0e10cSrcweir if (bConstant)
2841cdf0e10cSrcweir {
2842cdf0e10cSrcweir fSigmaIntercept = fRMSE
2843cdf0e10cSrcweir * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N));
2844cdf0e10cSrcweir pResMat->PutDouble(fSigmaIntercept, K, 1);
2845cdf0e10cSrcweir }
2846cdf0e10cSrcweir else
2847cdf0e10cSrcweir {
2848cdf0e10cSrcweir pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
2849cdf0e10cSrcweir }
2850cdf0e10cSrcweir
2851cdf0e10cSrcweir double fR2 = fSSreg / (fSSreg + fSSresid);
2852cdf0e10cSrcweir pResMat->PutDouble(fR2, 0, 2);
2853cdf0e10cSrcweir }
2854cdf0e10cSrcweir }
2855cdf0e10cSrcweir PushMatrix(pResMat);
2856cdf0e10cSrcweir }
2857cdf0e10cSrcweir }
2858cdf0e10cSrcweir return;
2859cdf0e10cSrcweir }
2860cdf0e10cSrcweir
ScTrend()2861cdf0e10cSrcweir void ScInterpreter::ScTrend()
2862cdf0e10cSrcweir {
2863cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScTrend" );
2864cdf0e10cSrcweir CalculateTrendGrowth(false);
2865cdf0e10cSrcweir }
2866cdf0e10cSrcweir
ScGrowth()2867cdf0e10cSrcweir void ScInterpreter::ScGrowth()
2868cdf0e10cSrcweir {
2869cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScGrowth" );
2870cdf0e10cSrcweir CalculateTrendGrowth(true);
2871cdf0e10cSrcweir }
2872cdf0e10cSrcweir
CalculateTrendGrowth(bool _bGrowth)2873cdf0e10cSrcweir void ScInterpreter::CalculateTrendGrowth(bool _bGrowth)
2874cdf0e10cSrcweir {
2875cdf0e10cSrcweir sal_uInt8 nParamCount = GetByte();
2876cdf0e10cSrcweir if ( !MustHaveParamCount( nParamCount, 1, 4 ) )
2877cdf0e10cSrcweir return;
2878cdf0e10cSrcweir
2879cdf0e10cSrcweir // optional fourth parameter
2880cdf0e10cSrcweir bool bConstant;
2881cdf0e10cSrcweir if (nParamCount == 4)
2882cdf0e10cSrcweir bConstant = GetBool();
2883cdf0e10cSrcweir else
2884cdf0e10cSrcweir bConstant = true;
2885cdf0e10cSrcweir
2886cdf0e10cSrcweir // The third parameter may be missing in ODF, although the fourth parameter
2887cdf0e10cSrcweir // is present. Default values depend on data not yet read.
2888cdf0e10cSrcweir ScMatrixRef pMatNewX;
2889cdf0e10cSrcweir if (nParamCount >= 3)
2890cdf0e10cSrcweir {
2891cdf0e10cSrcweir if (IsMissing())
2892cdf0e10cSrcweir {
2893cdf0e10cSrcweir Pop();
2894cdf0e10cSrcweir pMatNewX = NULL;
2895cdf0e10cSrcweir }
2896cdf0e10cSrcweir else
2897cdf0e10cSrcweir pMatNewX = GetMatrix();
2898cdf0e10cSrcweir }
2899cdf0e10cSrcweir else
2900cdf0e10cSrcweir pMatNewX = NULL;
2901cdf0e10cSrcweir
2902cdf0e10cSrcweir // In ODF1.2 empty second parameter (which is two ;; ) is allowed.
2903cdf0e10cSrcweir // Defaults will be set in CheckMatrix.
2904cdf0e10cSrcweir ScMatrixRef pMatX;
2905cdf0e10cSrcweir if (nParamCount >= 2)
2906cdf0e10cSrcweir {
2907cdf0e10cSrcweir if (IsMissing())
2908cdf0e10cSrcweir {
2909cdf0e10cSrcweir Pop();
2910cdf0e10cSrcweir pMatX = NULL;
2911cdf0e10cSrcweir }
2912cdf0e10cSrcweir else
2913cdf0e10cSrcweir {
2914cdf0e10cSrcweir pMatX = GetMatrix();
2915cdf0e10cSrcweir }
2916cdf0e10cSrcweir }
2917cdf0e10cSrcweir else
2918cdf0e10cSrcweir pMatX = NULL;
2919cdf0e10cSrcweir
2920cdf0e10cSrcweir ScMatrixRef pMatY;
2921cdf0e10cSrcweir pMatY = GetMatrix();
2922cdf0e10cSrcweir if (!pMatY)
2923cdf0e10cSrcweir {
2924cdf0e10cSrcweir PushIllegalParameter();
2925cdf0e10cSrcweir return;
2926cdf0e10cSrcweir }
2927cdf0e10cSrcweir
2928cdf0e10cSrcweir // 1 = simple; 2 = multiple with Y as column; 3 = multiple with Y as row
2929cdf0e10cSrcweir sal_uInt8 nCase;
2930cdf0e10cSrcweir
2931cdf0e10cSrcweir SCSIZE nCX, nCY; // number of columns
2932cdf0e10cSrcweir SCSIZE nRX, nRY; // number of rows
2933cdf0e10cSrcweir SCSIZE K = 0, N = 0; // K=number of variables X, N=number of data samples
2934cdf0e10cSrcweir if ( !CheckMatrix(_bGrowth,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY) )
2935cdf0e10cSrcweir {
2936cdf0e10cSrcweir PushIllegalParameter();
2937cdf0e10cSrcweir return;
2938cdf0e10cSrcweir }
2939cdf0e10cSrcweir
2940cdf0e10cSrcweir // Enough data samples?
2941cdf0e10cSrcweir if ( (bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1) )
2942cdf0e10cSrcweir {
2943cdf0e10cSrcweir PushIllegalParameter();
2944cdf0e10cSrcweir return;
2945cdf0e10cSrcweir }
2946cdf0e10cSrcweir
2947cdf0e10cSrcweir // Set default pMatNewX if necessary
2948cdf0e10cSrcweir SCSIZE nCXN, nRXN;
2949cdf0e10cSrcweir SCSIZE nCountXN;
2950cdf0e10cSrcweir if (!pMatNewX)
2951cdf0e10cSrcweir {
2952cdf0e10cSrcweir nCXN = nCX;
2953cdf0e10cSrcweir nRXN = nRX;
2954cdf0e10cSrcweir nCountXN = nCXN * nRXN;
2955cdf0e10cSrcweir pMatNewX = pMatX->Clone(); // pMatX will be changed to X-meanX
2956cdf0e10cSrcweir }
2957cdf0e10cSrcweir else
2958cdf0e10cSrcweir {
2959cdf0e10cSrcweir pMatNewX->GetDimensions(nCXN, nRXN);
2960cdf0e10cSrcweir if ((nCase == 2 && K != nCXN) || (nCase == 3 && K != nRXN))
2961cdf0e10cSrcweir {
2962cdf0e10cSrcweir PushIllegalArgument();
2963cdf0e10cSrcweir return;
2964cdf0e10cSrcweir }
2965cdf0e10cSrcweir nCountXN = nCXN * nRXN;
2966cdf0e10cSrcweir for ( SCSIZE i = 0; i < nCountXN; i++ )
2967cdf0e10cSrcweir if (!pMatNewX->IsValue(i))
2968cdf0e10cSrcweir {
2969cdf0e10cSrcweir PushIllegalArgument();
2970cdf0e10cSrcweir return;
2971cdf0e10cSrcweir }
2972cdf0e10cSrcweir }
2973cdf0e10cSrcweir ScMatrixRef pResMat; // size depends on nCase
2974cdf0e10cSrcweir if (nCase == 1)
2975cdf0e10cSrcweir pResMat = GetNewMat(nCXN,nRXN);
2976cdf0e10cSrcweir else
2977cdf0e10cSrcweir {
2978cdf0e10cSrcweir if (nCase==2)
2979cdf0e10cSrcweir pResMat = GetNewMat(1,nRXN);
2980cdf0e10cSrcweir else
2981cdf0e10cSrcweir pResMat = GetNewMat(nCXN,1);
2982cdf0e10cSrcweir }
2983cdf0e10cSrcweir if (!pResMat)
2984cdf0e10cSrcweir {
2985cdf0e10cSrcweir PushError(errCodeOverflow);
2986cdf0e10cSrcweir return;
2987cdf0e10cSrcweir }
2988cdf0e10cSrcweir // Uses sum(x-MeanX)^2 and not [sum x^2]-N * MeanX^2 in case bConstant.
2989cdf0e10cSrcweir // Clone constant matrices, so that Mat = Mat - Mean is possible.
2990cdf0e10cSrcweir double fMeanY = 0.0;
2991cdf0e10cSrcweir if (bConstant)
2992cdf0e10cSrcweir {
2993cdf0e10cSrcweir ScMatrixRef pCopyX = pMatX->CloneIfConst();
2994cdf0e10cSrcweir ScMatrixRef pCopyY = pMatY->CloneIfConst();
2995cdf0e10cSrcweir if (!pCopyX || !pCopyY)
2996cdf0e10cSrcweir {
2997cdf0e10cSrcweir PushError(errStackOverflow);
2998cdf0e10cSrcweir return;
2999cdf0e10cSrcweir }
3000cdf0e10cSrcweir pMatX = pCopyX;
3001cdf0e10cSrcweir pMatY = pCopyY;
3002cdf0e10cSrcweir // DeltaY is possible here; DeltaX depends on nCase, so later
3003cdf0e10cSrcweir fMeanY = lcl_GetMeanOverAll(pMatY, N);
3004cdf0e10cSrcweir for (SCSIZE i=0; i<N; i++)
3005cdf0e10cSrcweir {
3006cdf0e10cSrcweir pMatY->PutDouble( ::rtl::math::approxSub(pMatY->GetDouble(i),fMeanY), i );
3007cdf0e10cSrcweir }
3008cdf0e10cSrcweir }
3009cdf0e10cSrcweir
3010cdf0e10cSrcweir if (nCase==1)
3011cdf0e10cSrcweir {
3012cdf0e10cSrcweir // calculate simple regression
3013cdf0e10cSrcweir double fMeanX = 0.0;
3014cdf0e10cSrcweir if (bConstant)
3015cdf0e10cSrcweir { // Mat = Mat - Mean
3016cdf0e10cSrcweir fMeanX = lcl_GetMeanOverAll(pMatX, N);
3017cdf0e10cSrcweir for (SCSIZE i=0; i<N; i++)
3018cdf0e10cSrcweir {
3019cdf0e10cSrcweir pMatX->PutDouble( ::rtl::math::approxSub(pMatX->GetDouble(i),fMeanX), i );
3020cdf0e10cSrcweir }
3021cdf0e10cSrcweir }
3022cdf0e10cSrcweir double fSumXY = lcl_GetSumProduct(pMatX,pMatY,N);
3023cdf0e10cSrcweir double fSumX2 = lcl_GetSumProduct(pMatX,pMatX,N);
3024cdf0e10cSrcweir if (fSumX2==0.0)
3025cdf0e10cSrcweir {
3026cdf0e10cSrcweir PushNoValue(); // all x-values are identical
3027cdf0e10cSrcweir return;
3028cdf0e10cSrcweir }
3029cdf0e10cSrcweir double fSlope = fSumXY / fSumX2;
3030cdf0e10cSrcweir double fIntercept = 0.0;
3031cdf0e10cSrcweir double fHelp;
3032cdf0e10cSrcweir if (bConstant)
3033cdf0e10cSrcweir {
3034cdf0e10cSrcweir fIntercept = fMeanY - fSlope * fMeanX;
3035cdf0e10cSrcweir for (SCSIZE i = 0; i < nCountXN; i++)
3036cdf0e10cSrcweir {
3037cdf0e10cSrcweir fHelp = pMatNewX->GetDouble(i)*fSlope + fIntercept;
3038cdf0e10cSrcweir pResMat->PutDouble(_bGrowth ? exp(fHelp) : fHelp, i);
3039cdf0e10cSrcweir }
3040cdf0e10cSrcweir }
3041cdf0e10cSrcweir else
3042cdf0e10cSrcweir {
3043cdf0e10cSrcweir for (SCSIZE i = 0; i < nCountXN; i++)
3044cdf0e10cSrcweir {
3045cdf0e10cSrcweir fHelp = pMatNewX->GetDouble(i)*fSlope;
3046cdf0e10cSrcweir pResMat->PutDouble(_bGrowth ? exp(fHelp) : fHelp, i);
3047cdf0e10cSrcweir }
3048cdf0e10cSrcweir }
3049cdf0e10cSrcweir }
3050cdf0e10cSrcweir else // calculate multiple regression;
3051cdf0e10cSrcweir {
3052cdf0e10cSrcweir if (nCase ==2) // Y is column
3053cdf0e10cSrcweir {
3054cdf0e10cSrcweir ::std::vector< double> aVecR(N); // for QR decomposition
3055cdf0e10cSrcweir // Enough memory for needed matrices?
3056cdf0e10cSrcweir ScMatrixRef pMeans = GetNewMat(K, 1); // mean of each column
3057cdf0e10cSrcweir ScMatrixRef pSlopes = GetNewMat(1,K); // from b1 to bK
3058cdf0e10cSrcweir if (!pMeans || !pSlopes)
3059cdf0e10cSrcweir {
3060cdf0e10cSrcweir PushError(errCodeOverflow);
3061cdf0e10cSrcweir return;
3062cdf0e10cSrcweir }
3063cdf0e10cSrcweir if (bConstant)
3064cdf0e10cSrcweir {
3065cdf0e10cSrcweir lcl_CalculateColumnMeans(pMatX, pMeans, K, N);
3066cdf0e10cSrcweir lcl_CalculateColumnsDelta(pMatX, pMeans, K, N);
3067cdf0e10cSrcweir }
3068cdf0e10cSrcweir if (!lcl_CalculateQRdecomposition(pMatX, aVecR, K, N))
3069cdf0e10cSrcweir {
3070cdf0e10cSrcweir PushNoValue();
3071cdf0e10cSrcweir return;
3072cdf0e10cSrcweir }
3073cdf0e10cSrcweir // Later on we will divide by elements of aVecR, so make sure
3074cdf0e10cSrcweir // that they aren't zero.
3075cdf0e10cSrcweir bool bIsSingular=false;
3076cdf0e10cSrcweir for (SCSIZE row=0; row < K && !bIsSingular; row++)
3077cdf0e10cSrcweir bIsSingular = bIsSingular || aVecR[row]==0.0;
3078cdf0e10cSrcweir if (bIsSingular)
3079cdf0e10cSrcweir {
3080cdf0e10cSrcweir PushNoValue();
3081cdf0e10cSrcweir return;
3082cdf0e10cSrcweir }
3083cdf0e10cSrcweir // Z := Q' Y; Y is overwritten with result Z
3084cdf0e10cSrcweir for (SCSIZE col = 0; col < K; col++)
3085cdf0e10cSrcweir {
3086cdf0e10cSrcweir lcl_ApplyHouseholderTransformation(pMatX, col, pMatY, N);
3087cdf0e10cSrcweir }
3088cdf0e10cSrcweir // B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
3089cdf0e10cSrcweir // result Z should have zeros for index>=K; if not, ignore values
3090cdf0e10cSrcweir for (SCSIZE col = 0; col < K ; col++)
3091cdf0e10cSrcweir {
3092cdf0e10cSrcweir pSlopes->PutDouble( pMatY->GetDouble(col), col);
3093cdf0e10cSrcweir }
3094cdf0e10cSrcweir lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, false);
3095cdf0e10cSrcweir
3096cdf0e10cSrcweir // Fill result matrix
3097cdf0e10cSrcweir lcl_MFastMult(pMatNewX,pSlopes,pResMat,nRXN,K,1);
3098cdf0e10cSrcweir if (bConstant)
3099cdf0e10cSrcweir {
3100cdf0e10cSrcweir double fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K);
3101cdf0e10cSrcweir for (SCSIZE row = 0; row < nRXN; row++)
3102cdf0e10cSrcweir pResMat->PutDouble(pResMat->GetDouble(row)+fIntercept, row);
3103cdf0e10cSrcweir }
3104cdf0e10cSrcweir if (_bGrowth)
3105cdf0e10cSrcweir {
3106cdf0e10cSrcweir for (SCSIZE i = 0; i < nRXN; i++)
3107cdf0e10cSrcweir pResMat->PutDouble(exp(pResMat->GetDouble(i)), i);
3108cdf0e10cSrcweir }
3109cdf0e10cSrcweir }
3110cdf0e10cSrcweir else
3111cdf0e10cSrcweir { // nCase == 3, Y is row, all matrices are transposed
3112cdf0e10cSrcweir
3113cdf0e10cSrcweir ::std::vector< double> aVecR(N); // for QR decomposition
3114cdf0e10cSrcweir // Enough memory for needed matrices?
3115cdf0e10cSrcweir ScMatrixRef pMeans = GetNewMat(1, K); // mean of each row
3116cdf0e10cSrcweir ScMatrixRef pSlopes = GetNewMat(K,1); // row from b1 to bK
3117cdf0e10cSrcweir if (!pMeans || !pSlopes)
3118cdf0e10cSrcweir {
3119cdf0e10cSrcweir PushError(errCodeOverflow);
3120cdf0e10cSrcweir return;
3121cdf0e10cSrcweir }
3122cdf0e10cSrcweir if (bConstant)
3123cdf0e10cSrcweir {
3124cdf0e10cSrcweir lcl_CalculateRowMeans(pMatX, pMeans, N, K);
3125cdf0e10cSrcweir lcl_CalculateRowsDelta(pMatX, pMeans, N, K);
3126cdf0e10cSrcweir }
3127cdf0e10cSrcweir if (!lcl_TCalculateQRdecomposition(pMatX, aVecR, K, N))
3128cdf0e10cSrcweir {
3129cdf0e10cSrcweir PushNoValue();
3130cdf0e10cSrcweir return;
3131cdf0e10cSrcweir }
3132cdf0e10cSrcweir // Later on we will divide by elements of aVecR, so make sure
3133cdf0e10cSrcweir // that they aren't zero.
3134cdf0e10cSrcweir bool bIsSingular=false;
3135cdf0e10cSrcweir for (SCSIZE row=0; row < K && !bIsSingular; row++)
3136cdf0e10cSrcweir bIsSingular = bIsSingular || aVecR[row]==0.0;
3137cdf0e10cSrcweir if (bIsSingular)
3138cdf0e10cSrcweir {
3139cdf0e10cSrcweir PushNoValue();
3140cdf0e10cSrcweir return;
3141cdf0e10cSrcweir }
3142cdf0e10cSrcweir // Z := Q' Y; Y is overwritten with result Z
3143cdf0e10cSrcweir for (SCSIZE row = 0; row < K; row++)
3144cdf0e10cSrcweir {
3145cdf0e10cSrcweir lcl_TApplyHouseholderTransformation(pMatX, row, pMatY, N);
3146cdf0e10cSrcweir }
3147cdf0e10cSrcweir // B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
3148cdf0e10cSrcweir // result Z should have zeros for index>=K; if not, ignore values
3149cdf0e10cSrcweir for (SCSIZE col = 0; col < K ; col++)
3150cdf0e10cSrcweir {
3151cdf0e10cSrcweir pSlopes->PutDouble( pMatY->GetDouble(col), col);
3152cdf0e10cSrcweir }
3153cdf0e10cSrcweir lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, true);
3154cdf0e10cSrcweir
3155cdf0e10cSrcweir // Fill result matrix
3156cdf0e10cSrcweir lcl_MFastMult(pSlopes,pMatNewX,pResMat,1,K,nCXN);
3157cdf0e10cSrcweir if (bConstant)
3158cdf0e10cSrcweir {
3159cdf0e10cSrcweir double fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K);
3160cdf0e10cSrcweir for (SCSIZE col = 0; col < nCXN; col++)
3161cdf0e10cSrcweir pResMat->PutDouble(pResMat->GetDouble(col)+fIntercept, col);
3162cdf0e10cSrcweir }
3163cdf0e10cSrcweir if (_bGrowth)
3164cdf0e10cSrcweir {
3165cdf0e10cSrcweir for (SCSIZE i = 0; i < nCXN; i++)
3166cdf0e10cSrcweir pResMat->PutDouble(exp(pResMat->GetDouble(i)), i);
3167cdf0e10cSrcweir }
3168cdf0e10cSrcweir }
3169cdf0e10cSrcweir }
3170cdf0e10cSrcweir PushMatrix(pResMat);
3171cdf0e10cSrcweir return;
3172cdf0e10cSrcweir }
3173cdf0e10cSrcweir
ScMatRef()3174cdf0e10cSrcweir void ScInterpreter::ScMatRef()
3175cdf0e10cSrcweir {
3176cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatRef" );
3177cdf0e10cSrcweir // Falls Deltarefs drin sind...
3178cdf0e10cSrcweir Push( (FormulaToken&)*pCur );
3179cdf0e10cSrcweir ScAddress aAdr;
3180cdf0e10cSrcweir PopSingleRef( aAdr );
3181cdf0e10cSrcweir ScFormulaCell* pCell = (ScFormulaCell*) GetCell( aAdr );
3182cdf0e10cSrcweir if( pCell && pCell->GetCellType() == CELLTYPE_FORMULA )
3183cdf0e10cSrcweir {
3184cdf0e10cSrcweir const ScMatrix* pMat = pCell->GetMatrix();
3185cdf0e10cSrcweir if( pMat )
3186cdf0e10cSrcweir {
3187cdf0e10cSrcweir SCSIZE nCols, nRows;
3188cdf0e10cSrcweir pMat->GetDimensions( nCols, nRows );
3189cdf0e10cSrcweir SCSIZE nC = static_cast<SCSIZE>(aPos.Col() - aAdr.Col());
3190cdf0e10cSrcweir SCSIZE nR = static_cast<SCSIZE>(aPos.Row() - aAdr.Row());
3191cdf0e10cSrcweir if ((nCols <= nC && nCols != 1) || (nRows <= nR && nRows != 1))
3192cdf0e10cSrcweir PushNA();
3193cdf0e10cSrcweir else
3194cdf0e10cSrcweir {
3195cdf0e10cSrcweir ScMatValType nMatValType;
3196cdf0e10cSrcweir const ScMatrixValue* pMatVal = pMat->Get( nC, nR, nMatValType);
3197cdf0e10cSrcweir if (ScMatrix::IsNonValueType( nMatValType))
3198cdf0e10cSrcweir {
3199cdf0e10cSrcweir if (ScMatrix::IsEmptyPathType( nMatValType))
3200cdf0e10cSrcweir { // result of empty sal_False jump path
3201cdf0e10cSrcweir nFuncFmtType = NUMBERFORMAT_LOGICAL;
3202cdf0e10cSrcweir PushInt(0);
3203cdf0e10cSrcweir }
3204cdf0e10cSrcweir else if (ScMatrix::IsEmptyType( nMatValType))
3205cdf0e10cSrcweir {
3206cdf0e10cSrcweir // Not inherited (really?) and display as empty string, not 0.
3207cdf0e10cSrcweir PushTempToken( new ScEmptyCellToken( false, true));
3208cdf0e10cSrcweir }
3209cdf0e10cSrcweir else
3210cdf0e10cSrcweir PushString( pMatVal->GetString() );
3211cdf0e10cSrcweir }
3212cdf0e10cSrcweir else
3213cdf0e10cSrcweir {
3214cdf0e10cSrcweir PushDouble(pMatVal->fVal); // handles DoubleError
3215cdf0e10cSrcweir pDok->GetNumberFormatInfo( nCurFmtType, nCurFmtIndex, aAdr, pCell );
3216cdf0e10cSrcweir nFuncFmtType = nCurFmtType;
3217cdf0e10cSrcweir nFuncFmtIndex = nCurFmtIndex;
3218cdf0e10cSrcweir }
3219cdf0e10cSrcweir }
3220cdf0e10cSrcweir }
3221cdf0e10cSrcweir else
3222cdf0e10cSrcweir {
3223cdf0e10cSrcweir // If not a result matrix, obtain the cell value.
3224cdf0e10cSrcweir sal_uInt16 nErr = pCell->GetErrCode();
3225cdf0e10cSrcweir if (nErr)
3226cdf0e10cSrcweir PushError( nErr );
3227cdf0e10cSrcweir else if( pCell->IsValue() )
3228cdf0e10cSrcweir PushDouble( pCell->GetValue() );
3229cdf0e10cSrcweir else
3230cdf0e10cSrcweir {
3231cdf0e10cSrcweir String aVal;
3232cdf0e10cSrcweir pCell->GetString( aVal );
3233cdf0e10cSrcweir PushString( aVal );
3234cdf0e10cSrcweir }
3235cdf0e10cSrcweir pDok->GetNumberFormatInfo( nCurFmtType, nCurFmtIndex, aAdr, pCell );
3236cdf0e10cSrcweir nFuncFmtType = nCurFmtType;
3237cdf0e10cSrcweir nFuncFmtIndex = nCurFmtIndex;
3238cdf0e10cSrcweir }
3239cdf0e10cSrcweir }
3240cdf0e10cSrcweir else
3241cdf0e10cSrcweir PushError( errNoRef );
3242cdf0e10cSrcweir }
3243cdf0e10cSrcweir
ScInfo()3244cdf0e10cSrcweir void ScInterpreter::ScInfo()
3245cdf0e10cSrcweir {
3246cdf0e10cSrcweir RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScInfo" );
3247cdf0e10cSrcweir if( MustHaveParamCount( GetByte(), 1 ) )
3248cdf0e10cSrcweir {
3249cdf0e10cSrcweir String aStr = GetString();
3250cdf0e10cSrcweir ScCellKeywordTranslator::transKeyword(aStr, ScGlobal::GetLocale(), ocInfo);
3251cdf0e10cSrcweir if( aStr.EqualsAscii( "SYSTEM" ) )
3252cdf0e10cSrcweir PushString( String( RTL_CONSTASCII_USTRINGPARAM( SC_INFO_OSVERSION ) ) );
3253cdf0e10cSrcweir else if( aStr.EqualsAscii( "OSVERSION" ) )
3254cdf0e10cSrcweir PushString( String( RTL_CONSTASCII_USTRINGPARAM( "Windows (32-bit) NT 5.01" ) ) );
3255cdf0e10cSrcweir else if( aStr.EqualsAscii( "RELEASE" ) )
3256cdf0e10cSrcweir PushString( ::utl::Bootstrap::getBuildIdData( ::rtl::OUString() ) );
3257cdf0e10cSrcweir else if( aStr.EqualsAscii( "NUMFILE" ) )
3258cdf0e10cSrcweir PushDouble( 1 );
3259cdf0e10cSrcweir else if( aStr.EqualsAscii( "RECALC" ) )
3260cdf0e10cSrcweir PushString( ScGlobal::GetRscString( pDok->GetAutoCalc() ? STR_RECALC_AUTO : STR_RECALC_MANUAL ) );
3261cdf0e10cSrcweir else
3262cdf0e10cSrcweir PushIllegalArgument();
3263cdf0e10cSrcweir }
3264cdf0e10cSrcweir }
3265