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