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