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