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