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