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