sc/source/core/tool/interpr5.cxx

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