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