scilab/modules/polynomials/sci_gateway/cpp/sci_sfact.cpp

   1 /*
   2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3  * Copyright (C) 2012 - Scilab Enterprises - Cedric DELAMARRE
   4  *
   5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
   6  *
   7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
   8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
   9  * This file was originally licensed under the terms of the CeCILL v2.1,
  10  * and continues to be available under such terms.
  11  * For more information, see the COPYING file which you should have received
  12  * along with this program.
  13  *
  14  */
  15 /*--------------------------------------------------------------------------*/
  16 #include "polynomials_gw.hxx"
  17 #include "function.hxx"
  18 #include "double.hxx"
  19 #include "polynom.hxx"
  20 #include "overload.hxx"
  21
  22 extern "C"
  23 {
  24 #include "Scierror.h"
  25 #include "localization.h"
  26 #include "sciprint.h"
  27 #include "elem_common.h"
  28     extern int C2F(sfact1)(double*, int*, double*, int*, int*);
  29     extern int C2F(sfact2)(double*, int*, int*, double*, int*, int*);
  30 }
  31 /*--------------------------------------------------------------------------*/
  32 types::Function::ReturnValue sci_sfact(types::typed_list &in, int _iRetCount, types::typed_list &out)
  33 {
  34     types::Polynom* pPolyIn  = NULL;
  35     types::Polynom* pPolyOut = NULL;
  36
  37     int iMaxIt = 100;
  38     int iErr   = 0;
  39     int iOne   = 1;
  40
  41     if (in.size() != 1)
  42     {
  43         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "sfact", 1);
  44         return types::Function::Error;
  45     }
  46
  47     if (_iRetCount > 1)
  48     {
  49         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "sfact", 1);
  50         return types::Function::Error;
  51     }
  52
  53     if (in[0]->isPoly() == false)
  54     {
  55         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_sfact";
  56         return Overload::call(wstFuncName, in, _iRetCount, out);
  57     }
  58
  59     pPolyIn = in[0]->getAs<types::Polynom>();
  60
  61     if (pPolyIn->isComplex())
  62     {
  63         Scierror(999, _("%s: Wrong value for input argument #%d: A real matrix expected.\n"), "sfact", 1);
  64         return types::Function::Error;
  65     }
  66
  67     if (pPolyIn->isScalar())
  68     {
  69         double* pdblCoef = pPolyIn->get(0)->get();
  70
  71         // check symmetry
  72         int iDegree = pPolyIn->get(0)->getRank();
  73         int iDegD2  = (int)(iDegree / 2);
  74         int n       = 1 + iDegD2;
  75
  76         if (2 * iDegD2 != iDegree)
  77         {
  78             Scierror(999, _("%s: Wrong value for input argument #%d: A symmetric polynomial expected.\n"), "sfact", 1);
  79             return types::Function::Error;
  80         }
  81
  82         for (int i = 0; i < n; i++)
  83         {
  84             if (pdblCoef[i] != pdblCoef[iDegree - i])
  85             {
  86                 Scierror(999, _("%s: Wrong value for input argument #%d: A symmetric polynomial expected.\n"), "sfact", 1);
  87                 return types::Function::Error;
  88             }
  89         }
  90
  91         // create result
  92         pPolyOut = new types::Polynom(pPolyIn->getVariableName(), 1, 1);
  93         double* pdblCoefOut = NULL;
  94         types::SinglePoly* pSP = new types::SinglePoly(&pdblCoefOut, iDegD2);
  95         C2F(dcopy)(&n, pdblCoef, &iOne, pdblCoefOut, &iOne);
  96
  97         // perform operation
  98         double* pdblWork = new double[7 * n];
  99         C2F(sfact1)(pdblCoefOut, &iDegD2, pdblWork, &iMaxIt, &iErr);
 100         delete[] pdblWork;
 101         if (iErr == 2)
 102         {
 103             delete pSP;
 104             delete pPolyOut;
 105             Scierror(999, _("%s: Wrong value for input argument #%d: Non negative or null value expected at degree %d.\n"), "sfact", 1, n);
 106             return types::Function::Error;
 107         }
 108         else if (iErr == 1)
 109         {
 110             delete pSP;
 111             delete pPolyOut;
 112             Scierror(999, _("%s: Wrong value for input argument #%d: Convergence problem.\n"), "sfact", 1);
 113             return types::Function::Error;
 114         }
 115         else if (iErr < 0)
 116         {
 117             sciprint("%s: warning: Convergence at 10^%d near.\n", "sfact", iErr);
 118         }
 119
 120         // return result
 121         pPolyOut->set(0, pSP);
 122         delete pSP;
 123     }
 124     else
 125     {
 126         if (pPolyIn->getRows() != pPolyIn->getCols())
 127         {
 128             Scierror(999, _("%s: Wrong size for input argument #%d: Square matrix expected.\n"), "sfact", 1);
 129             return types::Function::Error;
 130         }
 131
 132         int iSize   = pPolyIn->getSize();
 133         int iRows   = pPolyIn->getRows();
 134         int iDegree = pPolyIn->getMaxRank();
 135         int iDegD2  = (int)(iDegree / 2);
 136         int n       = 1 + iDegD2;
 137
 138         double* pdblOut  = new double[iSize * n];
 139         double* pdblWork = new double[(n + 1) * iRows * ((n + 1)*iRows) + 1];
 140
 141         memset(pdblOut, 0x00, iSize * n * sizeof(double));
 142
 143         for (int i = 0; i < iSize; i++)
 144         {
 145             double* pdblIn = pPolyIn->get(i)->get();
 146             int iSizeToCpy = 2 + pPolyIn->get(i)->getSize() - 1 - n;
 147             if (iSizeToCpy > 0)
 148             {
 149                 C2F(dcopy)(&iSizeToCpy, pdblIn + n - 1, &iOne, pdblOut + i, &iSize);
 150             }
 151         }
 152
 153         int nm1 = n - 1;
 154         iMaxIt += n;
 155         C2F(sfact2)(pdblOut, &iRows, &nm1, pdblWork, &iMaxIt, &iErr);
 156         delete[] pdblWork;
 157
 158         if (iErr < 0)
 159         {
 160             delete[] pdblOut;
 161             Scierror(999, _("%s: Wrong value for input argument #%d: Convergence problem.\n"), "sfact", 1);
 162             return types::Function::Error;
 163         }
 164         else if (iErr > 0)
 165         {
 166             delete[] pdblOut;
 167             Scierror(999, _("%s: Wrong value for input argument #%d: singular or asymmetric problem.\n"), "sfact", 1);
 168             return types::Function::Error;
 169         }
 170
 171         pPolyOut = new types::Polynom(pPolyIn->getVariableName(), pPolyIn->getDims(), pPolyIn->getDimsArray());
 172         for (int i = 0; i < iSize; i++)
 173         {
 174             double* pdblCoefOut = NULL;
 175             types::SinglePoly* pSP = new types::SinglePoly(&pdblCoefOut, nm1);
 176             C2F(dcopy)(&n, pdblOut + i, &iSize, pdblCoefOut, &iOne);
 177             pPolyOut->set(i, pSP);
 178             delete pSP;
 179         }
 180
 181         delete[] pdblOut;
 182     }
 183
 184     out.push_back(pPolyOut);
 185     return types::Function::OK;
 186 }
 187 /*--------------------------------------------------------------------------*/
 188