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