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