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(pow(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(pow(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(pow(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 double lcl_GetSign(double fValue) 1940 { 1941 if (fValue < 0.0) 1942 return -1.0; 1943 else if (fValue > 0.0) 1944 return 1.0; 1945 else 1946 return 0.0; 1947 } 1948 1949 /* Calculates a QR decomposition with Householder reflection. 1950 * For each NxK matrix A exists a decomposition A=Q*R with an orthogonal 1951 * NxN matrix Q and a NxK matrix R. 1952 * Q=H1*H2*...*Hk with Householder matrices H. Such a householder matrix can 1953 * be build from a vector u by H=I-(2/u'u)*(u u'). This vectors u are returned 1954 * in the columns of matrix A, overwriting the old content. 1955 * The matrix R has a quadric upper part KxK with values in the upper right 1956 * triangle and zeros in all other elements. Here the diagonal elements of R 1957 * are stored in the vector R and the other upper right elements in the upper 1958 * right of the matrix A. 1959 * The function returns false, if calculation breaks. But because of round-off 1960 * errors singularity is often not detected. 1961 */ 1962 bool lcl_CalculateQRdecomposition(ScMatrixRef pMatA, 1963 ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN) 1964 { 1965 double fScale ; 1966 double fEuclid ; 1967 double fFactor ; 1968 double fSignum ; 1969 double fSum ; 1970 // ScMatrix matrices are zero based, index access (column,row) 1971 for (SCSIZE col = 0; col <nK; col++) 1972 { 1973 // calculate vector u of the householder transformation 1974 fScale = lcl_GetColumnMaximumNorm(pMatA, col, col, nN); 1975 if (fScale == 0.0) 1976 // A is singular 1977 return false; 1978 1979 for (SCSIZE row = col; row <nN; row++) 1980 pMatA->PutDouble( pMatA->GetDouble(col,row)/fScale, col, row); 1981 1982 fEuclid = lcl_GetColumnEuclideanNorm(pMatA, col, col, nN); 1983 fFactor = 1.0/fEuclid/(fEuclid + fabs(pMatA->GetDouble(col,col))); 1984 fSignum = lcl_GetSign(pMatA->GetDouble(col,col)); 1985 pMatA->PutDouble( pMatA->GetDouble(col,col) + fSignum*fEuclid, col,col); 1986 pVecR[col] = -fSignum * fScale * fEuclid; 1987 1988 // apply Householder transformation to A 1989 for (SCSIZE c=col+1; c<nK; c++) 1990 { 1991 fSum =lcl_GetColumnSumProduct(pMatA, col, pMatA, c, col, nN); 1992 for (SCSIZE row = col; row <nN; row++) 1993 pMatA->PutDouble( pMatA->GetDouble(c,row) 1994 - fSum * fFactor * pMatA->GetDouble(col,row), c, row); 1995 } 1996 } 1997 return true; 1998 } 1999 2000 // same with transposed matrix A, N is count of columns, K count of rows 2001 bool lcl_TCalculateQRdecomposition(ScMatrixRef pMatA, 2002 ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN) 2003 { 2004 double fScale ; 2005 double fEuclid ; 2006 double fFactor ; 2007 double fSignum ; 2008 double fSum ; 2009 // ScMatrix matrices are zero based, index access (column,row) 2010 for (SCSIZE row = 0; row <nK; row++) 2011 { 2012 // calculate vector u of the householder transformation 2013 fScale = lcl_TGetColumnMaximumNorm(pMatA, row, row, nN); 2014 if (fScale == 0.0) 2015 // A is singular 2016 return false; 2017 2018 for (SCSIZE col = row; col <nN; col++) 2019 pMatA->PutDouble( pMatA->GetDouble(col,row)/fScale, col, row); 2020 2021 fEuclid = lcl_TGetColumnEuclideanNorm(pMatA, row, row, nN); 2022 fFactor = 1.0/fEuclid/(fEuclid + fabs(pMatA->GetDouble(row,row))); 2023 fSignum = lcl_GetSign(pMatA->GetDouble(row,row)); 2024 pMatA->PutDouble( pMatA->GetDouble(row,row) + fSignum*fEuclid, row,row); 2025 pVecR[row] = -fSignum * fScale * fEuclid; 2026 2027 // apply Householder transformation to A 2028 for (SCSIZE r=row+1; r<nK; r++) 2029 { 2030 fSum =lcl_TGetColumnSumProduct(pMatA, row, pMatA, r, row, nN); 2031 for (SCSIZE col = row; col <nN; col++) 2032 pMatA->PutDouble( pMatA->GetDouble(col,r) 2033 - fSum * fFactor * pMatA->GetDouble(col,row), col, r); 2034 } 2035 } 2036 return true; 2037 } 2038 2039 2040 /* Applies a Householder transformation to a column vector Y with is given as 2041 * Nx1 Matrix. The Vektor u, from which the Householder transformation is build, 2042 * is the column part in matrix A, with column index C, starting with row 2043 * index C. A is the result of the QR decomposition as obtained from 2044 * lcl_CaluclateQRdecomposition. 2045 */ 2046 void lcl_ApplyHouseholderTransformation(ScMatrixRef pMatA, SCSIZE nC, 2047 ScMatrixRef pMatY, SCSIZE nN) 2048 { 2049 // ScMatrix matrices are zero based, index access (column,row) 2050 double fDenominator = lcl_GetColumnSumProduct(pMatA, nC, pMatA, nC, nC, nN); 2051 double fNumerator = lcl_GetColumnSumProduct(pMatA, nC, pMatY, 0, nC, nN); 2052 double fFactor = 2.0 * (fNumerator/fDenominator); 2053 for (SCSIZE row = nC; row < nN; row++) 2054 pMatY->PutDouble( 2055 pMatY->GetDouble(row) - fFactor * pMatA->GetDouble(nC,row), row); 2056 } 2057 2058 // Same with transposed matrices A and Y. 2059 void lcl_TApplyHouseholderTransformation(ScMatrixRef pMatA, SCSIZE nR, 2060 ScMatrixRef pMatY, SCSIZE nN) 2061 { 2062 // ScMatrix matrices are zero based, index access (column,row) 2063 double fDenominator = lcl_TGetColumnSumProduct(pMatA, nR, pMatA, nR, nR, nN); 2064 double fNumerator = lcl_TGetColumnSumProduct(pMatA, nR, pMatY, 0, nR, nN); 2065 double fFactor = 2.0 * (fNumerator/fDenominator); 2066 for (SCSIZE col = nR; col < nN; col++) 2067 pMatY->PutDouble( 2068 pMatY->GetDouble(col) - fFactor * pMatA->GetDouble(col,nR), col); 2069 } 2070 2071 /* Solve for X in R*X=S using back substitution. The solution X overwrites S. 2072 * Uses R from the result of the QR decomposition of a NxK matrix A. 2073 * S is a column vector given as matrix, with at least elements on index 2074 * 0 to K-1; elements on index>=K are ignored. Vector R must not have zero 2075 * elements, no check is done. 2076 */ 2077 void lcl_SolveWithUpperRightTriangle(ScMatrixRef pMatA, 2078 ::std::vector< double>& pVecR, ScMatrixRef pMatS, 2079 SCSIZE nK, bool bIsTransposed) 2080 { 2081 // ScMatrix matrices are zero based, index access (column,row) 2082 double fSum; 2083 SCSIZE row; 2084 // SCSIZE is never negative, therefore test with rowp1=row+1 2085 for (SCSIZE rowp1 = nK; rowp1>0; rowp1--) 2086 { 2087 row = rowp1-1; 2088 fSum = pMatS->GetDouble(row); 2089 for (SCSIZE col = rowp1; col<nK ; col++) 2090 if (bIsTransposed) 2091 fSum -= pMatA->GetDouble(row,col) * pMatS->GetDouble(col); 2092 else 2093 fSum -= pMatA->GetDouble(col,row) * pMatS->GetDouble(col); 2094 pMatS->PutDouble( fSum / pVecR[row] , row); 2095 } 2096 } 2097 2098 /* Solve for X in R' * X= T using forward substitution. The solution X 2099 * overwrites T. Uses R from the result of the QR decomposition of a NxK 2100 * matrix A. T is a column vectors given as matrix, with at least elements on 2101 * index 0 to K-1; elements on index>=K are ignored. Vector R must not have 2102 * zero elements, no check is done. 2103 */ 2104 void lcl_SolveWithLowerLeftTriangle(ScMatrixRef pMatA, 2105 ::std::vector< double>& pVecR, ScMatrixRef pMatT, 2106 SCSIZE nK, bool bIsTransposed) 2107 { 2108 // ScMatrix matrices are zero based, index access (column,row) 2109 double fSum; 2110 for (SCSIZE row = 0; row < nK; row++) 2111 { 2112 fSum = pMatT -> GetDouble(row); 2113 for (SCSIZE col=0; col < row; col++) 2114 { 2115 if (bIsTransposed) 2116 fSum -= pMatA->GetDouble(col,row) * pMatT->GetDouble(col); 2117 else 2118 fSum -= pMatA->GetDouble(row,col) * pMatT->GetDouble(col); 2119 } 2120 pMatT->PutDouble( fSum / pVecR[row] , row); 2121 } 2122 } 2123 2124 /* Calculates Z = R * B 2125 * R is given in matrix A and vector VecR as obtained from the QR 2126 * decompostion in lcl_CalculateQRdecomposition. B and Z are column vectors 2127 * given as matrix with at least index 0 to K-1; elements on index>=K are 2128 * not used. 2129 */ 2130 void lcl_ApplyUpperRightTriangle(ScMatrixRef pMatA, 2131 ::std::vector< double>& pVecR, ScMatrixRef pMatB, 2132 ScMatrixRef pMatZ, SCSIZE nK, bool bIsTransposed) 2133 { 2134 // ScMatrix matrices are zero based, index access (column,row) 2135 double fSum; 2136 for (SCSIZE row = 0; row < nK; row++) 2137 { 2138 fSum = pVecR[row] * pMatB->GetDouble(row); 2139 for (SCSIZE col = row+1; col < nK; col++) 2140 if (bIsTransposed) 2141 fSum += pMatA->GetDouble(row,col) * pMatB->GetDouble(col); 2142 else 2143 fSum += pMatA->GetDouble(col,row) * pMatB->GetDouble(col); 2144 pMatZ->PutDouble( fSum, row); 2145 } 2146 } 2147 2148 2149 2150 double lcl_GetMeanOverAll(ScMatrixRef pMat, SCSIZE nN) 2151 { 2152 double fSum = 0.0; 2153 for (SCSIZE i=0 ; i<nN; i++) 2154 fSum += pMat->GetDouble(i); 2155 return fSum/static_cast<double>(nN); 2156 } 2157 2158 // Calculates means of the columns of matrix X. X is a RxC matrix; 2159 // ResMat is a 1xC matrix (=row). 2160 void lcl_CalculateColumnMeans(ScMatrixRef pX, ScMatrixRef pResMat, 2161 SCSIZE nC, SCSIZE nR) 2162 { 2163 double fSum = 0.0; 2164 for (SCSIZE i=0; i < nC; i++) 2165 { 2166 fSum =0.0; 2167 for (SCSIZE k=0; k < nR; k++) 2168 fSum += pX->GetDouble(i,k); // GetDouble(Column,Row) 2169 pResMat ->PutDouble( fSum/static_cast<double>(nR),i); 2170 } 2171 } 2172 2173 // Calculates means of the rows of matrix X. X is a RxC matrix; 2174 // ResMat is a Rx1 matrix (=column). 2175 void lcl_CalculateRowMeans(ScMatrixRef pX, ScMatrixRef pResMat, 2176 SCSIZE nC, SCSIZE nR) 2177 { 2178 double fSum = 0.0; 2179 for (SCSIZE k=0; k < nR; k++) 2180 { 2181 fSum =0.0; 2182 for (SCSIZE i=0; i < nC; i++) 2183 fSum += pX->GetDouble(i,k); // GetDouble(Column,Row) 2184 pResMat ->PutDouble( fSum/static_cast<double>(nC),k); 2185 } 2186 } 2187 2188 void lcl_CalculateColumnsDelta(ScMatrixRef pMat, ScMatrixRef pColumnMeans, 2189 SCSIZE nC, SCSIZE nR) 2190 { 2191 for (SCSIZE i = 0; i < nC; i++) 2192 for (SCSIZE k = 0; k < nR; k++) 2193 pMat->PutDouble( ::rtl::math::approxSub 2194 (pMat->GetDouble(i,k) , pColumnMeans->GetDouble(i) ) , i, k); 2195 } 2196 2197 void lcl_CalculateRowsDelta(ScMatrixRef pMat, ScMatrixRef pRowMeans, 2198 SCSIZE nC, SCSIZE nR) 2199 { 2200 for (SCSIZE k = 0; k < nR; k++) 2201 for (SCSIZE i = 0; i < nC; i++) 2202 pMat->PutDouble( ::rtl::math::approxSub 2203 ( pMat->GetDouble(i,k) , pRowMeans->GetDouble(k) ) , i, k); 2204 } 2205 2206 // Case1 = simple regression 2207 // MatX = X - MeanX, MatY = Y - MeanY, y - haty = (y - MeanY) - (haty - MeanY) 2208 // = (y-MeanY)-((slope*x+a)-(slope*MeanX+a)) = (y-MeanY)-slope*(x-MeanX) 2209 double lcl_GetSSresid(ScMatrixRef pMatX, ScMatrixRef pMatY, double fSlope, 2210 SCSIZE nN) 2211 { 2212 double fSum = 0.0; 2213 double fTemp = 0.0; 2214 for (SCSIZE i=0; i<nN; i++) 2215 { 2216 fTemp = pMatY->GetDouble(i) - fSlope * pMatX->GetDouble(i); 2217 fSum += fTemp * fTemp; 2218 } 2219 return fSum; 2220 } 2221 2222 // Fill default values in matrix X, transform Y to log(Y) in case LOGEST|GROWTH, 2223 // determine sizes of matrices X and Y, determine kind of regression, clone 2224 // Y in case LOGEST|GROWTH, if constant. 2225 bool ScInterpreter::CheckMatrix(bool _bLOG, sal_uInt8& nCase, SCSIZE& nCX, 2226 SCSIZE& nCY, SCSIZE& nRX, SCSIZE& nRY, SCSIZE& M, 2227 SCSIZE& N, ScMatrixRef& pMatX, ScMatrixRef& pMatY) 2228 { 2229 2230 nCX = 0; 2231 nCY = 0; 2232 nRX = 0; 2233 nRY = 0; 2234 M = 0; 2235 N = 0; 2236 pMatY->GetDimensions(nCY, nRY); 2237 const SCSIZE nCountY = nCY * nRY; 2238 for ( SCSIZE i = 0; i < nCountY; i++ ) 2239 { 2240 if (!pMatY->IsValue(i)) 2241 { 2242 PushIllegalArgument(); 2243 return false; 2244 } 2245 } 2246 2247 if ( _bLOG ) 2248 { 2249 ScMatrixRef pNewY = pMatY->CloneIfConst(); 2250 for (SCSIZE nElem = 0; nElem < nCountY; nElem++) 2251 { 2252 const double fVal = pNewY->GetDouble(nElem); 2253 if (fVal <= 0.0) 2254 { 2255 PushIllegalArgument(); 2256 return false; 2257 } 2258 else 2259 pNewY->PutDouble(log(fVal), nElem); 2260 } 2261 pMatY = pNewY; 2262 } 2263 2264 if (pMatX) 2265 { 2266 pMatX->GetDimensions(nCX, nRX); 2267 const SCSIZE nCountX = nCX * nRX; 2268 for ( SCSIZE i = 0; i < nCountX; i++ ) 2269 if (!pMatX->IsValue(i)) 2270 { 2271 PushIllegalArgument(); 2272 return false; 2273 } 2274 if (nCX == nCY && nRX == nRY) 2275 { 2276 nCase = 1; // simple regression 2277 M = 1; 2278 N = nCountY; 2279 } 2280 else if (nCY != 1 && nRY != 1) 2281 { 2282 PushIllegalArgument(); 2283 return false; 2284 } 2285 else if (nCY == 1) 2286 { 2287 if (nRX != nRY) 2288 { 2289 PushIllegalArgument(); 2290 return false; 2291 } 2292 else 2293 { 2294 nCase = 2; // Y is column 2295 N = nRY; 2296 M = nCX; 2297 } 2298 } 2299 else if (nCX != nCY) 2300 { 2301 PushIllegalArgument(); 2302 return false; 2303 } 2304 else 2305 { 2306 nCase = 3; // Y is row 2307 N = nCY; 2308 M = nRX; 2309 } 2310 } 2311 else 2312 { 2313 pMatX = GetNewMat(nCY, nRY); 2314 nCX = nCY; 2315 nRX = nRY; 2316 if (!pMatX) 2317 { 2318 PushIllegalArgument(); 2319 return false; 2320 } 2321 for ( SCSIZE i = 1; i <= nCountY; i++ ) 2322 pMatX->PutDouble(static_cast<double>(i), i-1); 2323 nCase = 1; 2324 N = nCountY; 2325 M = 1; 2326 } 2327 return true; 2328 } 2329 // ----------------------------------------------------------------------------- 2330 2331 // LINEST 2332 void ScInterpreter::ScRGP() 2333 { 2334 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRGP" ); 2335 CalulateRGPRKP(false); 2336 } 2337 2338 // LOGEST 2339 void ScInterpreter::ScRKP() 2340 { 2341 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRKP" ); 2342 CalulateRGPRKP(true); 2343 } 2344 2345 void ScInterpreter::CalulateRGPRKP(bool _bRKP) 2346 { 2347 sal_uInt8 nParamCount = GetByte(); 2348 if ( !MustHaveParamCount( nParamCount, 1, 4 ) ) 2349 return; 2350 bool bConstant, bStats; 2351 2352 // optional forth parameter 2353 if (nParamCount == 4) 2354 bStats = GetBool(); 2355 else 2356 bStats = false; 2357 2358 // The third parameter may not be missing in ODF, if the forth parameter 2359 // is present. But Excel allows it with default true, we too. 2360 if (nParamCount >= 3) 2361 { 2362 if (IsMissing()) 2363 { 2364 Pop(); 2365 bConstant = true; 2366 // PushIllegalParameter(); if ODF behavior is desired 2367 // return; 2368 } 2369 else 2370 bConstant = GetBool(); 2371 } 2372 else 2373 bConstant = true; 2374 2375 ScMatrixRef pMatX; 2376 if (nParamCount >= 2) 2377 { 2378 if (IsMissing()) 2379 { 2380 // In ODF1.2 empty second parameter (which is two ;; ) is allowed 2381 Pop(); 2382 pMatX = NULL; 2383 } 2384 else 2385 { 2386 pMatX = GetMatrix(); 2387 } 2388 } 2389 else 2390 pMatX = NULL; 2391 2392 ScMatrixRef pMatY; 2393 pMatY = GetMatrix(); 2394 if (!pMatY) 2395 { 2396 PushIllegalParameter(); 2397 return; 2398 } 2399 2400 // 1 = simple; 2 = multiple with Y as column; 3 = multiple with Y as row 2401 sal_uInt8 nCase; 2402 2403 SCSIZE nCX, nCY; // number of columns 2404 SCSIZE nRX, nRY; // number of rows 2405 SCSIZE K = 0, N = 0; // K=number of variables X, N=number of data samples 2406 if ( !CheckMatrix(_bRKP,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY) ) 2407 { 2408 PushIllegalParameter(); 2409 return; 2410 } 2411 2412 // Enough data samples? 2413 if ( (bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1) ) 2414 { 2415 PushIllegalParameter(); 2416 return; 2417 } 2418 2419 ScMatrixRef pResMat; 2420 if (bStats) 2421 pResMat = GetNewMat(K+1,5); 2422 else 2423 pResMat = GetNewMat(K+1,1); 2424 if (!pResMat) 2425 { 2426 PushError(errCodeOverflow); 2427 return; 2428 } 2429 // Fill unused cells in pResMat; order (column,row) 2430 if (bStats) 2431 { 2432 for (SCSIZE i=2; i<K+1; i++) 2433 { 2434 pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 2 ); 2435 pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 3 ); 2436 pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 4 ); 2437 } 2438 } 2439 2440 // Uses sum(x-MeanX)^2 and not [sum x^2]-N * MeanX^2 in case bConstant. 2441 // Clone constant matrices, so that Mat = Mat - Mean is possible. 2442 double fMeanY = 0.0; 2443 if (bConstant) 2444 { 2445 ScMatrixRef pNewX = pMatX->CloneIfConst(); 2446 ScMatrixRef pNewY = pMatY->CloneIfConst(); 2447 if (!pNewX || !pNewY) 2448 { 2449 PushError(errCodeOverflow); 2450 return; 2451 } 2452 pMatX = pNewX; 2453 pMatY = pNewY; 2454 // DeltaY is possible here; DeltaX depends on nCase, so later 2455 fMeanY = lcl_GetMeanOverAll(pMatY, N); 2456 for (SCSIZE i=0; i<N; i++) 2457 { 2458 pMatY->PutDouble( ::rtl::math::approxSub(pMatY->GetDouble(i),fMeanY), i ); 2459 } 2460 } 2461 2462 if (nCase==1) 2463 { 2464 // calculate simple regression 2465 double fMeanX = 0.0; 2466 if (bConstant) 2467 { // Mat = Mat - Mean 2468 fMeanX = lcl_GetMeanOverAll(pMatX, N); 2469 for (SCSIZE i=0; i<N; i++) 2470 { 2471 pMatX->PutDouble( ::rtl::math::approxSub(pMatX->GetDouble(i),fMeanX), i ); 2472 } 2473 } 2474 double fSumXY = lcl_GetSumProduct(pMatX,pMatY,N); 2475 double fSumX2 = lcl_GetSumProduct(pMatX,pMatX,N); 2476 if (fSumX2==0.0) 2477 { 2478 PushNoValue(); // all x-values are identical 2479 return; 2480 } 2481 double fSlope = fSumXY / fSumX2; 2482 double fIntercept = 0.0; 2483 if (bConstant) 2484 fIntercept = fMeanY - fSlope * fMeanX; 2485 pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, 1, 0); //order (column,row) 2486 pResMat->PutDouble(_bRKP ? exp(fSlope) : fSlope, 0, 0); 2487 2488 if (bStats) 2489 { 2490 double fSSreg = fSlope * fSlope * fSumX2; 2491 pResMat->PutDouble(fSSreg, 0, 4); 2492 2493 double fDegreesFreedom =static_cast<double>( (bConstant) ? N-2 : N-1 ); 2494 pResMat->PutDouble(fDegreesFreedom, 1, 3); 2495 2496 double fSSresid = lcl_GetSSresid(pMatX,pMatY,fSlope,N); 2497 pResMat->PutDouble(fSSresid, 1, 4); 2498 2499 if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0) 2500 { // exact fit; test SSreg too, because SSresid might be 2501 // unequal zero due to round of errors 2502 pResMat->PutDouble(0.0, 1, 4); // SSresid 2503 pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 0, 3); // F 2504 pResMat->PutDouble(0.0, 1, 2); // RMSE 2505 pResMat->PutDouble(0.0, 0, 1); // SigmaSlope 2506 if (bConstant) 2507 pResMat->PutDouble(0.0, 1, 1); //SigmaIntercept 2508 else 2509 pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 1, 1); 2510 pResMat->PutDouble(1.0, 0, 2); // R^2 2511 } 2512 else 2513 { 2514 double fFstatistic = (fSSreg / static_cast<double>(K)) 2515 / (fSSresid / fDegreesFreedom); 2516 pResMat->PutDouble(fFstatistic, 0, 3); 2517 2518 // standard error of estimate 2519 double fRMSE = sqrt(fSSresid / fDegreesFreedom); 2520 pResMat->PutDouble(fRMSE, 1, 2); 2521 2522 double fSigmaSlope = fRMSE / sqrt(fSumX2); 2523 pResMat->PutDouble(fSigmaSlope, 0, 1); 2524 2525 if (bConstant) 2526 { 2527 double fSigmaIntercept = fRMSE 2528 * sqrt(fMeanX*fMeanX/fSumX2 + 1.0/static_cast<double>(N)); 2529 pResMat->PutDouble(fSigmaIntercept, 1, 1); 2530 } 2531 else 2532 { 2533 pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 1, 1); 2534 } 2535 2536 double fR2 = fSSreg / (fSSreg + fSSresid); 2537 pResMat->PutDouble(fR2, 0, 2); 2538 } 2539 } 2540 PushMatrix(pResMat); 2541 } 2542 else // calculate multiple regression; 2543 { 2544 // Uses a QR decomposition X = QR. The solution B = (X'X)^(-1) * X' * Y 2545 // becomes B = R^(-1) * Q' * Y 2546 if (nCase ==2) // Y is column 2547 { 2548 ::std::vector< double> aVecR(N); // for QR decomposition 2549 // Enough memory for needed matrices? 2550 ScMatrixRef pMeans = GetNewMat(K, 1); // mean of each column 2551 ScMatrixRef pMatZ; // for Q' * Y , inter alia 2552 if (bStats) 2553 pMatZ = pMatY->Clone(); // Y is used in statistic, keep it 2554 else 2555 pMatZ = pMatY; // Y can be overwritten 2556 ScMatrixRef pSlopes = GetNewMat(1,K); // from b1 to bK 2557 if (!pMeans || !pMatZ || !pSlopes) 2558 { 2559 PushError(errCodeOverflow); 2560 return; 2561 } 2562 if (bConstant) 2563 { 2564 lcl_CalculateColumnMeans(pMatX, pMeans, K, N); 2565 lcl_CalculateColumnsDelta(pMatX, pMeans, K, N); 2566 } 2567 if (!lcl_CalculateQRdecomposition(pMatX, aVecR, K, N)) 2568 { 2569 PushNoValue(); 2570 return; 2571 } 2572 // Later on we will divide by elements of aVecR, so make sure 2573 // that they aren't zero. 2574 bool bIsSingular=false; 2575 for (SCSIZE row=0; row < K && !bIsSingular; row++) 2576 bIsSingular = bIsSingular || aVecR[row]==0.0; 2577 if (bIsSingular) 2578 { 2579 PushNoValue(); 2580 return; 2581 } 2582 // Z = Q' Y; 2583 for (SCSIZE col = 0; col < K; col++) 2584 { 2585 lcl_ApplyHouseholderTransformation(pMatX, col, pMatZ, N); 2586 } 2587 // B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z 2588 // result Z should have zeros for index>=K; if not, ignore values 2589 for (SCSIZE col = 0; col < K ; col++) 2590 { 2591 pSlopes->PutDouble( pMatZ->GetDouble(col), col); 2592 } 2593 lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, false); 2594 double fIntercept = 0.0; 2595 if (bConstant) 2596 fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K); 2597 // Fill first line in result matrix 2598 pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, K, 0 ); 2599 for (SCSIZE i = 0; i < K; i++) 2600 pResMat->PutDouble(_bRKP ? exp(pSlopes->GetDouble(i)) 2601 : pSlopes->GetDouble(i) , K-1-i, 0); 2602 2603 2604 if (bStats) 2605 { 2606 double fSSreg = 0.0; 2607 double fSSresid = 0.0; 2608 // re-use memory of Z; 2609 pMatZ->FillDouble(0.0, 0, 0, 0, N-1); 2610 // Z = R * Slopes 2611 lcl_ApplyUpperRightTriangle(pMatX, aVecR, pSlopes, pMatZ, K, false); 2612 // Z = Q * Z, that is Q * R * Slopes = X * Slopes 2613 for (SCSIZE colp1 = K; colp1 > 0; colp1--) 2614 { 2615 lcl_ApplyHouseholderTransformation(pMatX, colp1-1, pMatZ,N); 2616 } 2617 fSSreg =lcl_GetSumProduct(pMatZ, pMatZ, N); 2618 // re-use Y for residuals, Y = Y-Z 2619 for (SCSIZE row = 0; row < N; row++) 2620 pMatY->PutDouble(pMatY->GetDouble(row) - pMatZ->GetDouble(row), row); 2621 fSSresid = lcl_GetSumProduct(pMatY, pMatY, N); 2622 pResMat->PutDouble(fSSreg, 0, 4); 2623 pResMat->PutDouble(fSSresid, 1, 4); 2624 2625 double fDegreesFreedom =static_cast<double>( (bConstant) ? N-K-1 : N-K ); 2626 pResMat->PutDouble(fDegreesFreedom, 1, 3); 2627 2628 if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0) 2629 { // exact fit; incl. observed values Y are identical 2630 pResMat->PutDouble(0.0, 1, 4); // SSresid 2631 // F = (SSreg/K) / (SSresid/df) = #DIV/0! 2632 pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 0, 3); // F 2633 // RMSE = sqrt(SSresid / df) = sqrt(0 / df) = 0 2634 pResMat->PutDouble(0.0, 1, 2); // RMSE 2635 // SigmaSlope[i] = RMSE * sqrt(matrix[i,i]) = 0 * sqrt(...) = 0 2636 for (SCSIZE i=0; i<K; i++) 2637 pResMat->PutDouble(0.0, K-1-i, 1); 2638 2639 // SigmaIntercept = RMSE * sqrt(...) = 0 2640 if (bConstant) 2641 pResMat->PutDouble(0.0, K, 1); //SigmaIntercept 2642 else 2643 pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1); 2644 2645 // R^2 = SSreg / (SSreg + SSresid) = 1.0 2646 pResMat->PutDouble(1.0, 0, 2); // R^2 2647 } 2648 else 2649 { 2650 double fFstatistic = (fSSreg / static_cast<double>(K)) 2651 / (fSSresid / fDegreesFreedom); 2652 pResMat->PutDouble(fFstatistic, 0, 3); 2653 2654 // standard error of estimate = root mean SSE 2655 double fRMSE = sqrt(fSSresid / fDegreesFreedom); 2656 pResMat->PutDouble(fRMSE, 1, 2); 2657 2658 // standard error of slopes 2659 // = RMSE * sqrt(diagonal element of (R' R)^(-1) ) 2660 // standard error of intercept 2661 // = RMSE * sqrt( Xmean * (R' R)^(-1) * Xmean' + 1/N) 2662 // (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as 2663 // a whole matrix, but iterate over unit vectors. 2664 double fSigmaSlope = 0.0; 2665 double fSigmaIntercept = 0.0; 2666 double fPart; // for Xmean * single column of (R' R)^(-1) 2667 for (SCSIZE col = 0; col < K; col++) 2668 { 2669 //re-use memory of MatZ 2670 pMatZ->FillDouble(0.0,0,0,0,K-1); // Z = unit vector e 2671 pMatZ->PutDouble(1.0, col); 2672 //Solve R' * Z = e 2673 lcl_SolveWithLowerLeftTriangle(pMatX, aVecR, pMatZ, K, false); 2674 // Solve R * Znew = Zold 2675 lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pMatZ, K, false); 2676 // now Z is column col in (R' R)^(-1) 2677 fSigmaSlope = fRMSE * sqrt(pMatZ->GetDouble(col)); 2678 pResMat->PutDouble(fSigmaSlope, K-1-col, 1); 2679 // (R' R) ^(-1) is symmetric 2680 if (bConstant) 2681 { 2682 fPart = lcl_GetSumProduct(pMeans, pMatZ, K); 2683 fSigmaIntercept += fPart * pMeans->GetDouble(col); 2684 } 2685 } 2686 if (bConstant) 2687 { 2688 fSigmaIntercept = fRMSE 2689 * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N)); 2690 pResMat->PutDouble(fSigmaIntercept, K, 1); 2691 } 2692 else 2693 { 2694 pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1); 2695 } 2696 2697 double fR2 = fSSreg / (fSSreg + fSSresid); 2698 pResMat->PutDouble(fR2, 0, 2); 2699 } 2700 } 2701 PushMatrix(pResMat); 2702 } 2703 else // nCase == 3, Y is row, all matrices are transposed 2704 { 2705 ::std::vector< double> aVecR(N); // for QR decomposition 2706 // Enough memory for needed matrices? 2707 ScMatrixRef pMeans = GetNewMat(1, K); // mean of each row 2708 ScMatrixRef pMatZ; // for Q' * Y , inter alia 2709 if (bStats) 2710 pMatZ = pMatY->Clone(); // Y is used in statistic, keep it 2711 else 2712 pMatZ = pMatY; // Y can be overwritten 2713 ScMatrixRef pSlopes = GetNewMat(K,1); // from b1 to bK 2714 if (!pMeans || !pMatZ || !pSlopes) 2715 { 2716 PushError(errCodeOverflow); 2717 return; 2718 } 2719 if (bConstant) 2720 { 2721 lcl_CalculateRowMeans(pMatX, pMeans, N, K); 2722 lcl_CalculateRowsDelta(pMatX, pMeans, N, K); 2723 } 2724 2725 if (!lcl_TCalculateQRdecomposition(pMatX, aVecR, K, N)) 2726 { 2727 PushNoValue(); 2728 return; 2729 } 2730 2731 // Later on we will divide by elements of aVecR, so make sure 2732 // that they aren't zero. 2733 bool bIsSingular=false; 2734 for (SCSIZE row=0; row < K && !bIsSingular; row++) 2735 bIsSingular = bIsSingular || aVecR[row]==0.0; 2736 if (bIsSingular) 2737 { 2738 PushNoValue(); 2739 return; 2740 } 2741 // Z = Q' Y 2742 for (SCSIZE row = 0; row < K; row++) 2743 { 2744 lcl_TApplyHouseholderTransformation(pMatX, row, pMatZ, N); 2745 } 2746 // B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z 2747 // result Z should have zeros for index>=K; if not, ignore values 2748 for (SCSIZE col = 0; col < K ; col++) 2749 { 2750 pSlopes->PutDouble( pMatZ->GetDouble(col), col); 2751 } 2752 lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, true); 2753 double fIntercept = 0.0; 2754 if (bConstant) 2755 fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K); 2756 // Fill first line in result matrix 2757 pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, K, 0 ); 2758 for (SCSIZE i = 0; i < K; i++) 2759 pResMat->PutDouble(_bRKP ? exp(pSlopes->GetDouble(i)) 2760 : pSlopes->GetDouble(i) , K-1-i, 0); 2761 2762 2763 if (bStats) 2764 { 2765 double fSSreg = 0.0; 2766 double fSSresid = 0.0; 2767 // re-use memory of Z; 2768 pMatZ->FillDouble(0.0, 0, 0, N-1, 0); 2769 // Z = R * Slopes 2770 lcl_ApplyUpperRightTriangle(pMatX, aVecR, pSlopes, pMatZ, K, true); 2771 // Z = Q * Z, that is Q * R * Slopes = X * Slopes 2772 for (SCSIZE rowp1 = K; rowp1 > 0; rowp1--) 2773 { 2774 lcl_TApplyHouseholderTransformation(pMatX, rowp1-1, pMatZ,N); 2775 } 2776 fSSreg =lcl_GetSumProduct(pMatZ, pMatZ, N); 2777 // re-use Y for residuals, Y = Y-Z 2778 for (SCSIZE col = 0; col < N; col++) 2779 pMatY->PutDouble(pMatY->GetDouble(col) - pMatZ->GetDouble(col), col); 2780 fSSresid = lcl_GetSumProduct(pMatY, pMatY, N); 2781 pResMat->PutDouble(fSSreg, 0, 4); 2782 pResMat->PutDouble(fSSresid, 1, 4); 2783 2784 double fDegreesFreedom =static_cast<double>( (bConstant) ? N-K-1 : N-K ); 2785 pResMat->PutDouble(fDegreesFreedom, 1, 3); 2786 2787 if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0) 2788 { // exact fit; incl. case observed values Y are identical 2789 pResMat->PutDouble(0.0, 1, 4); // SSresid 2790 // F = (SSreg/K) / (SSresid/df) = #DIV/0! 2791 pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 0, 3); // F 2792 // RMSE = sqrt(SSresid / df) = sqrt(0 / df) = 0 2793 pResMat->PutDouble(0.0, 1, 2); // RMSE 2794 // SigmaSlope[i] = RMSE * sqrt(matrix[i,i]) = 0 * sqrt(...) = 0 2795 for (SCSIZE i=0; i<K; i++) 2796 pResMat->PutDouble(0.0, K-1-i, 1); 2797 2798 // SigmaIntercept = RMSE * sqrt(...) = 0 2799 if (bConstant) 2800 pResMat->PutDouble(0.0, K, 1); //SigmaIntercept 2801 else 2802 pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1); 2803 2804 // R^2 = SSreg / (SSreg + SSresid) = 1.0 2805 pResMat->PutDouble(1.0, 0, 2); // R^2 2806 } 2807 else 2808 { 2809 double fFstatistic = (fSSreg / static_cast<double>(K)) 2810 / (fSSresid / fDegreesFreedom); 2811 pResMat->PutDouble(fFstatistic, 0, 3); 2812 2813 // standard error of estimate = root mean SSE 2814 double fRMSE = sqrt(fSSresid / fDegreesFreedom); 2815 pResMat->PutDouble(fRMSE, 1, 2); 2816 2817 // standard error of slopes 2818 // = RMSE * sqrt(diagonal element of (R' R)^(-1) ) 2819 // standard error of intercept 2820 // = RMSE * sqrt( Xmean * (R' R)^(-1) * Xmean' + 1/N) 2821 // (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as 2822 // a whole matrix, but iterate over unit vectors. 2823 // (R' R) ^(-1) is symmetric 2824 double fSigmaSlope = 0.0; 2825 double fSigmaIntercept = 0.0; 2826 double fPart; // for Xmean * single col of (R' R)^(-1) 2827 for (SCSIZE row = 0; row < K; row++) 2828 { 2829 //re-use memory of MatZ 2830 pMatZ->FillDouble(0.0,0,0,K-1,0); // Z = unit vector e 2831 pMatZ->PutDouble(1.0, row); 2832 //Solve R' * Z = e 2833 lcl_SolveWithLowerLeftTriangle(pMatX, aVecR, pMatZ, K, true); 2834 // Solve R * Znew = Zold 2835 lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pMatZ, K, true); 2836 // now Z is column col in (R' R)^(-1) 2837 fSigmaSlope = fRMSE * sqrt(pMatZ->GetDouble(row)); 2838 pResMat->PutDouble(fSigmaSlope, K-1-row, 1); 2839 if (bConstant) 2840 { 2841 fPart = lcl_GetSumProduct(pMeans, pMatZ, K); 2842 fSigmaIntercept += fPart * pMeans->GetDouble(row); 2843 } 2844 } 2845 if (bConstant) 2846 { 2847 fSigmaIntercept = fRMSE 2848 * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N)); 2849 pResMat->PutDouble(fSigmaIntercept, K, 1); 2850 } 2851 else 2852 { 2853 pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1); 2854 } 2855 2856 double fR2 = fSSreg / (fSSreg + fSSresid); 2857 pResMat->PutDouble(fR2, 0, 2); 2858 } 2859 } 2860 PushMatrix(pResMat); 2861 } 2862 } 2863 return; 2864 } 2865 2866 void ScInterpreter::ScTrend() 2867 { 2868 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScTrend" ); 2869 CalculateTrendGrowth(false); 2870 } 2871 2872 void ScInterpreter::ScGrowth() 2873 { 2874 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScGrowth" ); 2875 CalculateTrendGrowth(true); 2876 } 2877 2878 void ScInterpreter::CalculateTrendGrowth(bool _bGrowth) 2879 { 2880 sal_uInt8 nParamCount = GetByte(); 2881 if ( !MustHaveParamCount( nParamCount, 1, 4 ) ) 2882 return; 2883 2884 // optional fourth parameter 2885 bool bConstant; 2886 if (nParamCount == 4) 2887 bConstant = GetBool(); 2888 else 2889 bConstant = true; 2890 2891 // The third parameter may be missing in ODF, although the fourth parameter 2892 // is present. Default values depend on data not yet read. 2893 ScMatrixRef pMatNewX; 2894 if (nParamCount >= 3) 2895 { 2896 if (IsMissing()) 2897 { 2898 Pop(); 2899 pMatNewX = NULL; 2900 } 2901 else 2902 pMatNewX = GetMatrix(); 2903 } 2904 else 2905 pMatNewX = NULL; 2906 2907 // In ODF1.2 empty second parameter (which is two ;; ) is allowed. 2908 // Defaults will be set in CheckMatrix. 2909 ScMatrixRef pMatX; 2910 if (nParamCount >= 2) 2911 { 2912 if (IsMissing()) 2913 { 2914 Pop(); 2915 pMatX = NULL; 2916 } 2917 else 2918 { 2919 pMatX = GetMatrix(); 2920 } 2921 } 2922 else 2923 pMatX = NULL; 2924 2925 ScMatrixRef pMatY; 2926 pMatY = GetMatrix(); 2927 if (!pMatY) 2928 { 2929 PushIllegalParameter(); 2930 return; 2931 } 2932 2933 // 1 = simple; 2 = multiple with Y as column; 3 = multiple with Y as row 2934 sal_uInt8 nCase; 2935 2936 SCSIZE nCX, nCY; // number of columns 2937 SCSIZE nRX, nRY; // number of rows 2938 SCSIZE K = 0, N = 0; // K=number of variables X, N=number of data samples 2939 if ( !CheckMatrix(_bGrowth,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY) ) 2940 { 2941 PushIllegalParameter(); 2942 return; 2943 } 2944 2945 // Enough data samples? 2946 if ( (bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1) ) 2947 { 2948 PushIllegalParameter(); 2949 return; 2950 } 2951 2952 // Set default pMatNewX if necessary 2953 SCSIZE nCXN, nRXN; 2954 SCSIZE nCountXN; 2955 if (!pMatNewX) 2956 { 2957 nCXN = nCX; 2958 nRXN = nRX; 2959 nCountXN = nCXN * nRXN; 2960 pMatNewX = pMatX->Clone(); // pMatX will be changed to X-meanX 2961 } 2962 else 2963 { 2964 pMatNewX->GetDimensions(nCXN, nRXN); 2965 if ((nCase == 2 && K != nCXN) || (nCase == 3 && K != nRXN)) 2966 { 2967 PushIllegalArgument(); 2968 return; 2969 } 2970 nCountXN = nCXN * nRXN; 2971 for ( SCSIZE i = 0; i < nCountXN; i++ ) 2972 if (!pMatNewX->IsValue(i)) 2973 { 2974 PushIllegalArgument(); 2975 return; 2976 } 2977 } 2978 ScMatrixRef pResMat; // size depends on nCase 2979 if (nCase == 1) 2980 pResMat = GetNewMat(nCXN,nRXN); 2981 else 2982 { 2983 if (nCase==2) 2984 pResMat = GetNewMat(1,nRXN); 2985 else 2986 pResMat = GetNewMat(nCXN,1); 2987 } 2988 if (!pResMat) 2989 { 2990 PushError(errCodeOverflow); 2991 return; 2992 } 2993 // Uses sum(x-MeanX)^2 and not [sum x^2]-N * MeanX^2 in case bConstant. 2994 // Clone constant matrices, so that Mat = Mat - Mean is possible. 2995 double fMeanY = 0.0; 2996 if (bConstant) 2997 { 2998 ScMatrixRef pCopyX = pMatX->CloneIfConst(); 2999 ScMatrixRef pCopyY = pMatY->CloneIfConst(); 3000 if (!pCopyX || !pCopyY) 3001 { 3002 PushError(errStackOverflow); 3003 return; 3004 } 3005 pMatX = pCopyX; 3006 pMatY = pCopyY; 3007 // DeltaY is possible here; DeltaX depends on nCase, so later 3008 fMeanY = lcl_GetMeanOverAll(pMatY, N); 3009 for (SCSIZE i=0; i<N; i++) 3010 { 3011 pMatY->PutDouble( ::rtl::math::approxSub(pMatY->GetDouble(i),fMeanY), i ); 3012 } 3013 } 3014 3015 if (nCase==1) 3016 { 3017 // calculate simple regression 3018 double fMeanX = 0.0; 3019 if (bConstant) 3020 { // Mat = Mat - Mean 3021 fMeanX = lcl_GetMeanOverAll(pMatX, N); 3022 for (SCSIZE i=0; i<N; i++) 3023 { 3024 pMatX->PutDouble( ::rtl::math::approxSub(pMatX->GetDouble(i),fMeanX), i ); 3025 } 3026 } 3027 double fSumXY = lcl_GetSumProduct(pMatX,pMatY,N); 3028 double fSumX2 = lcl_GetSumProduct(pMatX,pMatX,N); 3029 if (fSumX2==0.0) 3030 { 3031 PushNoValue(); // all x-values are identical 3032 return; 3033 } 3034 double fSlope = fSumXY / fSumX2; 3035 double fIntercept = 0.0; 3036 double fHelp; 3037 if (bConstant) 3038 { 3039 fIntercept = fMeanY - fSlope * fMeanX; 3040 for (SCSIZE i = 0; i < nCountXN; i++) 3041 { 3042 fHelp = pMatNewX->GetDouble(i)*fSlope + fIntercept; 3043 pResMat->PutDouble(_bGrowth ? exp(fHelp) : fHelp, i); 3044 } 3045 } 3046 else 3047 { 3048 for (SCSIZE i = 0; i < nCountXN; i++) 3049 { 3050 fHelp = pMatNewX->GetDouble(i)*fSlope; 3051 pResMat->PutDouble(_bGrowth ? exp(fHelp) : fHelp, i); 3052 } 3053 } 3054 } 3055 else // calculate multiple regression; 3056 { 3057 if (nCase ==2) // Y is column 3058 { 3059 ::std::vector< double> aVecR(N); // for QR decomposition 3060 // Enough memory for needed matrices? 3061 ScMatrixRef pMeans = GetNewMat(K, 1); // mean of each column 3062 ScMatrixRef pSlopes = GetNewMat(1,K); // from b1 to bK 3063 if (!pMeans || !pSlopes) 3064 { 3065 PushError(errCodeOverflow); 3066 return; 3067 } 3068 if (bConstant) 3069 { 3070 lcl_CalculateColumnMeans(pMatX, pMeans, K, N); 3071 lcl_CalculateColumnsDelta(pMatX, pMeans, K, N); 3072 } 3073 if (!lcl_CalculateQRdecomposition(pMatX, aVecR, K, N)) 3074 { 3075 PushNoValue(); 3076 return; 3077 } 3078 // Later on we will divide by elements of aVecR, so make sure 3079 // that they aren't zero. 3080 bool bIsSingular=false; 3081 for (SCSIZE row=0; row < K && !bIsSingular; row++) 3082 bIsSingular = bIsSingular || aVecR[row]==0.0; 3083 if (bIsSingular) 3084 { 3085 PushNoValue(); 3086 return; 3087 } 3088 // Z := Q' Y; Y is overwritten with result Z 3089 for (SCSIZE col = 0; col < K; col++) 3090 { 3091 lcl_ApplyHouseholderTransformation(pMatX, col, pMatY, N); 3092 } 3093 // B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z 3094 // result Z should have zeros for index>=K; if not, ignore values 3095 for (SCSIZE col = 0; col < K ; col++) 3096 { 3097 pSlopes->PutDouble( pMatY->GetDouble(col), col); 3098 } 3099 lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, false); 3100 3101 // Fill result matrix 3102 lcl_MFastMult(pMatNewX,pSlopes,pResMat,nRXN,K,1); 3103 if (bConstant) 3104 { 3105 double fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K); 3106 for (SCSIZE row = 0; row < nRXN; row++) 3107 pResMat->PutDouble(pResMat->GetDouble(row)+fIntercept, row); 3108 } 3109 if (_bGrowth) 3110 { 3111 for (SCSIZE i = 0; i < nRXN; i++) 3112 pResMat->PutDouble(exp(pResMat->GetDouble(i)), i); 3113 } 3114 } 3115 else 3116 { // nCase == 3, Y is row, all matrices are transposed 3117 3118 ::std::vector< double> aVecR(N); // for QR decomposition 3119 // Enough memory for needed matrices? 3120 ScMatrixRef pMeans = GetNewMat(1, K); // mean of each row 3121 ScMatrixRef pSlopes = GetNewMat(K,1); // row from b1 to bK 3122 if (!pMeans || !pSlopes) 3123 { 3124 PushError(errCodeOverflow); 3125 return; 3126 } 3127 if (bConstant) 3128 { 3129 lcl_CalculateRowMeans(pMatX, pMeans, N, K); 3130 lcl_CalculateRowsDelta(pMatX, pMeans, N, K); 3131 } 3132 if (!lcl_TCalculateQRdecomposition(pMatX, aVecR, K, N)) 3133 { 3134 PushNoValue(); 3135 return; 3136 } 3137 // Later on we will divide by elements of aVecR, so make sure 3138 // that they aren't zero. 3139 bool bIsSingular=false; 3140 for (SCSIZE row=0; row < K && !bIsSingular; row++) 3141 bIsSingular = bIsSingular || aVecR[row]==0.0; 3142 if (bIsSingular) 3143 { 3144 PushNoValue(); 3145 return; 3146 } 3147 // Z := Q' Y; Y is overwritten with result Z 3148 for (SCSIZE row = 0; row < K; row++) 3149 { 3150 lcl_TApplyHouseholderTransformation(pMatX, row, pMatY, N); 3151 } 3152 // B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z 3153 // result Z should have zeros for index>=K; if not, ignore values 3154 for (SCSIZE col = 0; col < K ; col++) 3155 { 3156 pSlopes->PutDouble( pMatY->GetDouble(col), col); 3157 } 3158 lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, true); 3159 3160 // Fill result matrix 3161 lcl_MFastMult(pSlopes,pMatNewX,pResMat,1,K,nCXN); 3162 if (bConstant) 3163 { 3164 double fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K); 3165 for (SCSIZE col = 0; col < nCXN; col++) 3166 pResMat->PutDouble(pResMat->GetDouble(col)+fIntercept, col); 3167 } 3168 if (_bGrowth) 3169 { 3170 for (SCSIZE i = 0; i < nCXN; i++) 3171 pResMat->PutDouble(exp(pResMat->GetDouble(i)), i); 3172 } 3173 } 3174 } 3175 PushMatrix(pResMat); 3176 return; 3177 } 3178 3179 void ScInterpreter::ScMatRef() 3180 { 3181 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScMatRef" ); 3182 // Falls Deltarefs drin sind... 3183 Push( (FormulaToken&)*pCur ); 3184 ScAddress aAdr; 3185 PopSingleRef( aAdr ); 3186 ScFormulaCell* pCell = (ScFormulaCell*) GetCell( aAdr ); 3187 if( pCell && pCell->GetCellType() == CELLTYPE_FORMULA ) 3188 { 3189 const ScMatrix* pMat = pCell->GetMatrix(); 3190 if( pMat ) 3191 { 3192 SCSIZE nCols, nRows; 3193 pMat->GetDimensions( nCols, nRows ); 3194 SCSIZE nC = static_cast<SCSIZE>(aPos.Col() - aAdr.Col()); 3195 SCSIZE nR = static_cast<SCSIZE>(aPos.Row() - aAdr.Row()); 3196 if ((nCols <= nC && nCols != 1) || (nRows <= nR && nRows != 1)) 3197 PushNA(); 3198 else 3199 { 3200 ScMatValType nMatValType; 3201 const ScMatrixValue* pMatVal = pMat->Get( nC, nR, nMatValType); 3202 if (ScMatrix::IsNonValueType( nMatValType)) 3203 { 3204 if (ScMatrix::IsEmptyPathType( nMatValType)) 3205 { // result of empty sal_False jump path 3206 nFuncFmtType = NUMBERFORMAT_LOGICAL; 3207 PushInt(0); 3208 } 3209 else if (ScMatrix::IsEmptyType( nMatValType)) 3210 { 3211 // Not inherited (really?) and display as empty string, not 0. 3212 PushTempToken( new ScEmptyCellToken( false, true)); 3213 } 3214 else 3215 PushString( pMatVal->GetString() ); 3216 } 3217 else 3218 { 3219 PushDouble(pMatVal->fVal); // handles DoubleError 3220 pDok->GetNumberFormatInfo( nCurFmtType, nCurFmtIndex, aAdr, pCell ); 3221 nFuncFmtType = nCurFmtType; 3222 nFuncFmtIndex = nCurFmtIndex; 3223 } 3224 } 3225 } 3226 else 3227 { 3228 // If not a result matrix, obtain the cell value. 3229 sal_uInt16 nErr = pCell->GetErrCode(); 3230 if (nErr) 3231 PushError( nErr ); 3232 else if( pCell->IsValue() ) 3233 PushDouble( pCell->GetValue() ); 3234 else 3235 { 3236 String aVal; 3237 pCell->GetString( aVal ); 3238 PushString( aVal ); 3239 } 3240 pDok->GetNumberFormatInfo( nCurFmtType, nCurFmtIndex, aAdr, pCell ); 3241 nFuncFmtType = nCurFmtType; 3242 nFuncFmtIndex = nCurFmtIndex; 3243 } 3244 } 3245 else 3246 PushError( errNoRef ); 3247 } 3248 3249 void ScInterpreter::ScInfo() 3250 { 3251 RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScInfo" ); 3252 if( MustHaveParamCount( GetByte(), 1 ) ) 3253 { 3254 String aStr = GetString(); 3255 ScCellKeywordTranslator::transKeyword(aStr, ScGlobal::GetLocale(), ocInfo); 3256 if( aStr.EqualsAscii( "SYSTEM" ) ) 3257 PushString( String( RTL_CONSTASCII_USTRINGPARAM( SC_INFO_OSVERSION ) ) ); 3258 else if( aStr.EqualsAscii( "OSVERSION" ) ) 3259 PushString( String( RTL_CONSTASCII_USTRINGPARAM( "Windows (32-bit) NT 5.01" ) ) ); 3260 else if( aStr.EqualsAscii( "RELEASE" ) ) 3261 PushString( ::utl::Bootstrap::getBuildIdData( ::rtl::OUString() ) ); 3262 else if( aStr.EqualsAscii( "NUMFILE" ) ) 3263 PushDouble( 1 ); 3264 else if( aStr.EqualsAscii( "RECALC" ) ) 3265 PushString( ScGlobal::GetRscString( pDok->GetAutoCalc() ? STR_RECALC_AUTO : STR_RECALC_MANUAL ) ); 3266 else 3267 PushIllegalArgument(); 3268 } 3269 } 3270