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