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