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