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