scilab/modules/polynomials/sci_gateway/cpp/sci_roots.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
  17 #include <cmath>
  18
  19 #include "polynomials_gw.hxx"
  20 #include "function.hxx"
  21 #include "double.hxx"
  22 #include "string.hxx"
  23 #include "polynom.hxx"
  24 #include "overload.hxx"
  25 #include "context.hxx"
  26 #include "execvisitor.hxx"
  27
  28 extern "C"
  29 {
  30 #include "Scierror.h"
  31 #include "localization.h"
  32 #include "elem_common.h"
  33     extern int C2F(rpoly)(double*, int*, double*, double*, int*);
  34     extern int C2F(wscal)(int*, double*, double*, double*, double*, int*);
  35 }
  36 /*--------------------------------------------------------------------------*/
  37 types::Function::ReturnValue sci_roots(types::typed_list &in, int _iRetCount, types::typed_list &out)
  38 {
  39     std::wstring wstrAlgo   = L"e"; // e = eigen (default), f = fast
  40     types::Double* pDblIn   = NULL;
  41     types::Double* pDblOut  = NULL;
  42     types::Polynom* pPolyIn = NULL;
  43
  44     double* pdblInReal   = NULL;
  45     double* pdblInImg    = NULL;
  46     double* pdblTempReal = NULL;
  47     double* pdblTempImg  = NULL;
  48
  49     int iOne = 1;
  50     int imOne = -1;
  51     int iSize = 0;
  52     bool bComplex = false;
  53     types::Function::ReturnValue ret = types::Function::OK;
  54
  55     if (in.size() < 1 || in.size() > 2)
  56     {
  57         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "roots", 1, 2);
  58         return types::Function::Error;
  59     }
  60
  61     if (_iRetCount > 1)
  62     {
  63         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "roots", 1);
  64         return types::Function::Error;
  65     }
  66
  67     // get algo type
  68     if (in.size() == 2)
  69     {
  70         if (in[1]->isString() == false)
  71         {
  72             Scierror(999, _("%s: Wrong type for input argument #%d : string expected.\n"), "roots", 2);
  73             return types::Function::Error;
  74         }
  75
  76         types::String* pStrAlgo = in[1]->getAs<types::String>();
  77         if (pStrAlgo->isScalar() == false)
  78         {
  79             Scierror(999, _("%s: Wrong size for input argument #%d : A scalar expected.\n"), "roots", 2);
  80             return types::Function::Error;
  81         }
  82
  83         wstrAlgo = pStrAlgo->get(0);
  84         if (wstrAlgo != L"e" && wstrAlgo != L"f")
  85         {
  86             Scierror(999, _("%s: Wrong value for input argument #%d : ""%s"" or ""%s"" expected.\n"), "roots", 2, "e", "f");
  87             return types::Function::Error;
  88         }
  89     }
  90
  91     if (in[0]->isDouble())
  92     {
  93         // for Matlab compatibility root of the vector of coefficients
  94         pDblIn = in[0]->getAs<types::Double>();
  95         if (pDblIn->isEmpty())
  96         {
  97             out.push_back(types::Double::Empty());
  98             return types::Function::OK;
  99         }
 100
 101         iSize = pDblIn->getSize();
 102
 103         // old fortran function dtild
 104         // switch elements of a vector. [1 2 3] => [3 2 1]
 105         pdblInReal = new double[iSize];
 106         C2F(dcopy)(&iSize, pDblIn->get(), &iOne, pdblInReal, &imOne);
 107         if (pDblIn->isComplex())
 108         {
 109             bComplex = true;
 110             pdblInImg = new double[iSize];
 111             C2F(dcopy)(&iSize, pDblIn->getImg(), &iOne, pdblInImg, &imOne);
 112         }
 113     }
 114     else if (in[0]->isPoly())
 115     {
 116         pPolyIn = in[0]->getAs<types::Polynom>();
 117
 118         if (pPolyIn->isScalar() == false)
 119         {
 120             Scierror(999, _("%s: Wrong size for input argument #%d : A scalar expected.\n"), "roots", 1);
 121             return types::Function::Error;
 122         }
 123
 124         iSize      = pPolyIn->getMaxRank() + 1;
 125         pdblInReal = pPolyIn->get(0)->get();
 126
 127         if (pPolyIn->isComplex())
 128         {
 129             bComplex = true;
 130             pdblInImg = pPolyIn->get(0)->getImg();
 131         }
 132     }
 133     else
 134     {
 135         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_roots";
 136         return Overload::call(wstFuncName, in, _iRetCount, out);
 137     }
 138
 139     // If "fast" algo was chosen and polynomial is complex,
 140     // then produce an error.
 141     if (wstrAlgo == L"f" && bComplex)
 142     {
 143         Scierror(999, _("%s: Wrong value for input argument #%d : If algo is ""%s"", a real is expected.\n"), "roots", 2, "f");
 144         return types::Function::Error;
 145     }
 146
 147     double t = 0;
 148     while (t == 0)
 149     {
 150         iSize--;
 151         if (iSize < 0)
 152         {
 153             out.push_back(types::Double::Empty());
 154             return types::Function::OK;
 155         }
 156
 157         t = std::fabs(pdblInReal[iSize]);
 158         if (bComplex)
 159         {
 160             t += std::fabs(pdblInImg[iSize]);
 161         }
 162     }
 163
 164     // If "fast" algo was chosen and polynomial has degree greater than 100,
 165     // then produce an error.
 166     if (wstrAlgo == L"f" && iSize > 100)
 167     {
 168         Scierror(999, _("%s: Wrong value for input argument #%d : If algo is ""%s"", a degree less than %d expected.\n"), "roots", 2, "f", 100);
 169         return types::Function::Error;
 170     }
 171
 172     if (wstrAlgo == L"f")
 173     {
 174         // real polynomial: rpoly algorithm
 175         // this alg is much more speedy, but it may happens that it gives
 176         // erroneous results without messages : example
 177         // roots(%s^31-8*%s^30+9*%s^29+0.995) should have two real roots near
 178         // 1.355 and 6.65 and the other ones inside a circle centered in 0
 179         // with radius 1
 180
 181         pDblOut = new types::Double(iSize, 1, true);
 182         int iFail = 0;
 183         int iSizeP1 = iSize + 1;
 184
 185         double* pdblTempReal = new double[iSize + 1];
 186         C2F(dcopy)(&iSizeP1, pdblInReal, &iOne, pdblTempReal, &imOne);
 187         C2F(rpoly)(pdblTempReal, &iSize, pDblOut->get(), pDblOut->getImg(), &iFail);
 188         delete[] pdblTempReal;
 189
 190         if (iFail)
 191         {
 192             if (iFail == 1)
 193             {
 194                 Scierror(999, _("%s: Convergence problem...\n"), "roots");
 195             }
 196             else if (iFail == 2)
 197             {
 198                 Scierror(999, _("%s: Leading coefficient is zero.\n"), "roots");
 199             }
 200             else if (iFail == 3)
 201             {
 202                 Scierror(999, _("%s: Too high degree (max 100).\n"), "roots");
 203             }
 204
 205             return types::Function::Error;
 206         }
 207
 208         out.push_back(pDblOut);
 209     }
 210     else // wstrAlgo == L"e"
 211     {
 212         // Companion matrix method
 213         int iSizeM1 = iSize - 1;
 214         int iSizeP1 = iSize + 1;
 215         double dblOne = 1;
 216         double* pdblTempReal = new double[iSize];
 217         double* pdblTempImg = NULL;
 218         double sr = pdblInReal[iSize];
 219
 220         C2F(dcopy)(&iSize, pdblInReal, &iOne, pdblTempReal, &imOne);
 221
 222         if (bComplex)
 223         {
 224             pdblTempImg = new double[iSize];
 225             C2F(dcopy)(&iSize, pdblInImg, &iOne, pdblTempImg, &imOne);
 226
 227             double si = pdblInImg[iSize];
 228             double t = sr * sr + si * si;
 229             sr = -sr / t;
 230             si = si / t;
 231
 232             C2F(wscal)(&iSize, &sr, &si, pdblTempReal, pdblTempImg, &iOne);
 233         }
 234         else
 235         {
 236             double dbl = -1 / sr;
 237             C2F(dscal)(&iSize, &dbl, pdblTempReal, &iOne);
 238         }
 239
 240         pDblOut = new types::Double(iSize, iSize, bComplex);
 241         double* pdblOutReal = pDblOut->get();
 242         double* pdblOutImg = NULL;
 243
 244         memset(pdblOutReal, 0x00, pDblOut->getSize() * sizeof(double));
 245         C2F(dset)(&iSizeM1, &dblOne, &pdblOutReal[iSize], &iSizeP1);
 246         C2F(dcopy)(&iSize, pdblTempReal, &iOne, pdblOutReal, &iOne);
 247         delete[] pdblTempReal;
 248
 249         if (bComplex)
 250         {
 251             pdblOutImg = pDblOut->getImg();
 252             memset(pdblOutImg, 0x00, pDblOut->getSize() * sizeof(double));
 253             C2F(dcopy)(&iSize, pdblTempImg, &iOne, pdblOutImg, &iOne);
 254             delete[] pdblTempImg;
 255         }
 256
 257         //call spec
 258         types::InternalType* pSpec = symbol::Context::getInstance()->get(symbol::Symbol(L"spec"));
 259         if (pSpec && pSpec->isFunction())
 260         {
 261             types::Function *funcSpec = pSpec->getAs<types::Function>();
 262             types::typed_list tlInput;
 263             types::optional_list tlOpt;
 264             tlInput.push_back(pDblOut);
 265
 266             ret = funcSpec->call(tlInput, tlOpt, 1, out);
 267             pDblOut->killMe();
 268         }
 269         else
 270         {
 271             Scierror(999, _("%s: unable to find spec function\n"), "roots");
 272             return types::Function::Error;
 273         }
 274     }
 275
 276     if (pDblIn)
 277     {
 278         delete[] pdblInReal;
 279         if (bComplex)
 280         {
 281             delete[] pdblInImg;
 282         }
 283     }
 284
 285     return ret;
 286 }
 287 /*--------------------------------------------------------------------------*/
 288