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