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