[elementary_functions] linspace() c++ gateway gives tremendous speedup
[scilab.git] / scilab / modules / elementary_functions / sci_gateway / cpp / sci_linspace.cpp
1 /*
2  *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  *  Copyright (C) 2018- St├ęphane MOTTELET
4  *
5  * This file is hereby licensed under the terms of the GNU GPL v2.0,
6  * For more information, see the COPYING file which you should have received
7  * along with this program.
8  *
9  */
10
11 #include "double.hxx"
12 #include "function.hxx"
13 #include "int.hxx"
14 #include "overload.hxx"
15
16 extern "C"
17 {
18 #include "Scierror.h"
19 #include "localization.h"
20 }
21
22 bool fillRange(double* pdblOut, double* pdblMin, double* pdblMax, int iRows, int iCols);
23
24 /* ==================================================================== */
25 types::Function::ReturnValue sci_linspace(types::typed_list &in, int _iRetCount, types::typed_list &out)
26 {
27     int iCols = 100;
28     types::Double* pDblOut;
29
30     if (in.size() != 2 & in.size() != 3)
31     {
32         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "linspace", 2,3);
33         return types::Function::Error;
34     }
35
36     if (_iRetCount > 1)
37     {
38         Scierror(78, _("%s: Wrong number of output argument(s): %d expected."), "linspace", 1);
39         return types::Function::Error;
40     }
41
42     types::Double* pDbl[2];
43     for (int i=0; i<2; i++)
44     {
45         if (in[i]->isDouble())
46         {
47             pDbl[i] = in[i]->getAs<types::Double>();
48         }
49         else
50         {
51             // other types -> overload
52             std::wstring wstFuncName = L"%" + in[i]->getShortTypeStr() + L"_linspace";
53             return Overload::call(wstFuncName, in, _iRetCount, out);
54         }
55     }
56
57     // Check dimensions are the same
58     int  iDims0  = pDbl[0]->getDims();
59     int* piDims0 = pDbl[0]->getDimsArray();
60     int  iDims1  = pDbl[1]->getDims();
61     int* piDims1 = pDbl[1]->getDimsArray();
62     if (iDims0 != iDims1)
63     {
64         Scierror(999, _("%s: Arguments %d and %d must have same dimensions.\n"), "linspace",1,2);
65         return types::Function::Error;
66     }
67     for (int i = 0; i < iDims0; i++)
68     {
69         if (piDims0[i] != piDims1[i])
70         {
71             Scierror(999, _("%s: Arguments %d and %d must have same dimensions.\n"), "linspace",1,2);
72             return types::Function::Error;
73         }
74     }
75
76     if (in.size() == 3)
77     {
78         if (in[2]->isDouble() && in[2]->getAs<types::Double>()->isComplex() == false
79             && in[2]->getAs<types::Double>()->getSize() == 1)
80         {
81             double dblCols = in[2]->getAs<types::Double>()->get(0);
82             if (std::floor(dblCols) != dblCols)
83             {
84                 Scierror(999, _("%s: Argument #%d: An integer value expected.\n"), "linspace",3);
85                 return types::Function::Error;
86             }
87             if (dblCols <= 0)
88             {
89                 // empty matrix case
90                 out.push_back(types::Double::Empty());
91                 return types::Function::OK;
92             }
93             iCols = (int) dblCols;
94         }
95         else
96         {
97             Scierror(999, _("%s: Argument #%d: An integer value expected.\n"), "linspace",3);
98             return types::Function::Error;
99         }
100     }
101
102     // yet another empty matrix case
103     if (pDbl[0]->getSize() == 0)
104     {
105         out.push_back(types::Double::Empty());
106         return types::Function::OK;
107     }
108
109     // generation is done by considering array as a column vector
110     int iRows =  pDbl[0]->getSize();
111     // pDblOut is resized later
112     pDblOut = new types::Double(iRows,iCols);
113
114     if (!fillRange(pDblOut->get(), pDbl[0]->get(), pDbl[1]->get(), iRows, iCols))
115     {
116         // if Infs or NaNs
117         pDblOut->killMe();
118         return types::Function::Error;
119     }
120
121     if (pDbl[0]->isComplex() || pDbl[1]->isComplex())
122     {
123         int iReal;
124         for (iReal = 0; iReal < 2; iReal++) {
125             if (!pDbl[iReal]->isComplex())
126             {
127                 pDbl[iReal] = pDbl[iReal]->clone();
128                 pDbl[iReal]->setComplex(true);
129                 break;
130             }
131         }
132         // Complexify pDblOut
133         pDblOut->setComplex(true);
134         bool status = fillRange(pDblOut->getImg(), pDbl[0]->getImg(), pDbl[1]->getImg(), iRows, iCols);
135         if (iReal < 2)
136         {
137               pDbl[iReal]->killMe();
138         }
139         if (status != true) // if Infs or NaNs
140         {
141             pDblOut->killMe();
142             return types::Function::Error;
143         }
144     }
145
146     int *piNewDims = new int[iDims0+1];
147     // keep the first dimension unchanged
148     piNewDims[0] = piDims0[0];
149     int iDim = 1;
150     for (int i=1; i<iDims0; i++)
151     {
152         // squeeze subsequent single dimensions
153         if (piDims0[i]>1)
154         {
155             piNewDims[iDim++] = piDims0[i];
156         }
157     }
158     // add the suplementary dimension
159     piNewDims[iDim++] = iCols;
160     // reshape the matrix:
161     pDblOut->reshape(piNewDims, iDim);
162     out.push_back(pDblOut);
163
164     return types::Function::OK;
165 }
166
167 bool fillRange(double* pdblOut, double* pdblMin, double* pdblMax, int iRows, int iCols)
168 {
169     double* step = new double[iRows];
170     for (int j = 0, k = (iCols-1)*iRows; j < iRows; j++)
171     {
172         step[j] = (pdblMax[j]-pdblMin[j])/(iCols-1);
173         // checking Infs and NaNs
174         int indInfOrNan = std::isinf(pdblMin[j]) || std::isnan(pdblMin[j]) ? 1 : 0;
175         indInfOrNan = indInfOrNan == 0 ? (std::isinf(pdblMax[j]) || std::isnan(pdblMax[j]) ? 2 : 0) : indInfOrNan;
176         if (indInfOrNan > 0)
177         {
178             delete[] step;
179             Scierror(999, _("%s: Argument #%d: %%nan and %%inf values are forbidden.\n"), "linspace",indInfOrNan);
180             return false;
181         }
182         // last column is enforced (http://bugzilla.scilab.org/10966)
183         pdblOut[k++] = pdblMax[j];
184     }
185     // doing the linear range generation
186     for (int i = 0; i<iCols-1; i++)
187     {
188         for (int j = 0; j < iRows; j++)
189         {
190               *(pdblOut++) = pdblMin[j]+i*step[j];
191         }
192     }
193     delete[] step;
194
195     return true;
196 }
197
198