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 extern void C2F(recbez)(double*, int*, double*, int*, double*, int*, double*, double*);
27 }
28 /*--------------------------------------------------------------------------*/
sci_bezout(types::typed_list & in,int _iRetCount,types::typed_list & out)29 types::Function::ReturnValue sci_bezout(types::typed_list &in, int _iRetCount, types::typed_list &out)
30 {
31 bool bComplex = false;
32 double* pdblIn[2] = {NULL, NULL};
33 int piDegree[2] = {0, 0};
34 double dblErr = 0;
35 int iOne = 1;
36
37 std::wstring wstrName = L"";
38
39 if (in.size() != 2)
40 {
41 Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "bezout", 2);
42 return types::Function::Error;
43 }
44
45 if (_iRetCount > 3)
46 {
47 Scierror(78, _("%s: Wrong number of output argument(s): %d to %d expected.\n"), "bezout", 1, 3);
48 return types::Function::Error;
49 // _iRetCount==3 is undocumented. May be it is used in an internal way
50 }
51
52 // get input arguments
53 for (int i = 0; i < in.size(); i++)
54 {
55 if (in[i]->isPoly() == false && in[i]->isDouble() == false)
56 {
57 std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_bezout";
58 return Overload::call(wstFuncName, in, _iRetCount, out);
59 }
60
61 types::GenericType* pGT = in[i]->getAs<types::GenericType>();
62
63 if (pGT->isComplex())
64 {
65 Scierror(999, _("%s: Wrong type for input argument #%d: A real scalar expected.\n"), "bezout", i + 1);
66 return types::Function::Error;
67 }
68
69 if (pGT->isScalar() == false)
70 {
71 Scierror(999, _("%s: Wrong size for input argument #%d: A real scalar expected.\n"), "bezout", i + 1);
72 return types::Function::Error;
73 }
74
75 if (in[i]->isDouble())
76 {
77 pdblIn[i] = in[i]->getAs<types::Double>()->get();
78 }
79 else // polynom
80 {
81 types::Polynom* pPolyIn = in[i]->getAs<types::Polynom>();
82 if (wstrName != L"" && wstrName != pPolyIn->getVariableName())
83 {
84 Scierror(999, _("%s: Wrong value for input argument #%d: A polynomial '%ls' expected.\n"), "bezout", i + 1, wstrName.c_str());
85 return types::Function::Error;
86 }
87
88 wstrName = pPolyIn->getVariableName();
89 pdblIn[i] = pPolyIn->get(0)->get();
90 piDegree[i] = pPolyIn->get(0)->getRank();
91 }
92 }
93
94 // perform operation
95 int iMaxRank = std::max(piDegree[0], piDegree[1]) + 1;
96 int iMinRank = std::min(piDegree[0], piDegree[1]) + 1;
97 double* pdblWork = new double[10 * iMaxRank + 3 * iMaxRank * iMaxRank];
98 double* pdblOut = new double[2 * (piDegree[0] + piDegree[1] + 2) + iMinRank + 3];
99 int ipb[6];
100
101 C2F(recbez)(pdblIn[0], piDegree, pdblIn[1], piDegree + 1, pdblOut, ipb, pdblWork, &dblErr);
102 delete[] pdblWork;
103
104 // create result
105 int np = ipb[1] - ipb[0];
106 double* pdblSP = NULL;
107 types::SinglePoly* pSP = new types::SinglePoly(&pdblSP, np - 1);
108 memcpy(pdblSP, pdblOut + ipb[0] - 1, np * sizeof(double));
109
110 if (wstrName == L"")
111 {
112 wstrName = L"s";
113 }
114
115 types::Polynom* pPolyGCD = new types::Polynom(wstrName, 1, 1);
116 pPolyGCD->set(0, pSP);
117 delete pSP;
118
119 // return result
120 out.push_back(pPolyGCD);
121
122 if (_iRetCount > 1)
123 {
124 types::Polynom* pPolyU = new types::Polynom(wstrName, 2, 2);
125 for (int i = 0; i < 4; i++)
126 {
127 int ii = i + 1;
128 int iRankU = ipb[ii + 1] - ipb[ii];
129 double* pdblU = NULL;
130 types::SinglePoly* pSPU = new types::SinglePoly(&pdblU, iRankU - 1);
131 memcpy(pdblU, pdblOut + ipb[ii] - 1, iRankU * sizeof(double));
132 pPolyU->set(i, pSPU);
133 delete pSPU;
134 }
135 out.push_back(pPolyU);
136 }
137
138 delete[] pdblOut;
139
140
141 if (_iRetCount == 3)
142 {
143 out.push_back(new types::Double(dblErr));
144 }
145
146 return types::Function::OK;
147 }
148 /*--------------------------------------------------------------------------*/
149