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