1 /*
2  *    This file is part of CasADi.
3  *
4  *    CasADi -- A symbolic framework for dynamic optimization.
5  *    Copyright (C) 2010-2014 Joel Andersson, Joris Gillis, Moritz Diehl,
6  *                            K.U. Leuven. All rights reserved.
7  *    Copyright (C) 2011-2014 Greg Horn
8  *
9  *    CasADi is free software; you can redistribute it and/or
10  *    modify it under the terms of the GNU Lesser General Public
11  *    License as published by the Free Software Foundation; either
12  *    version 3 of the License, or (at your option) any later version.
13  *
14  *    CasADi is distributed in the hope that it will be useful,
15  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
16  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  *    Lesser General Public License for more details.
18  *
19  *    You should have received a copy of the GNU Lesser General Public
20  *    License along with CasADi; if not, write to the Free Software
21  *    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22  *
23  */
24 
25 
26 #ifndef CASADI_CALCULUS_HPP
27 #define CASADI_CALCULUS_HPP
28 
29 #include <iostream>
30 #include <string>
31 #include <cmath>
32 #include <limits>
33 #include <algorithm>
34 #include "casadi_common.hpp"
35 
36 // Define pi if the compiler fails to do so
37 
38 /// \cond INTERNAL
39 
40 namespace casadi {
41 #ifndef SWIG
42   /// Define pi
43 #ifdef M_PI
44   const double pi = M_PI;
45 #else
46   const double pi = 3.14159265358979323846;
47 #endif
48 
49   /// infinity
50   const double inf = std::numeric_limits<double>::infinity();
51 
52   /// Not a number
53   const double nan = std::numeric_limits<double>::quiet_NaN();
54 
55   /// Machine epsilon
56   const double eps = std::numeric_limits<double>::epsilon();
57 #endif // SWIG
58 
59   /// Enum for quick access to any node
60   enum Operation {
61     // Simple assignment
62     OP_ASSIGN,
63 
64     // Standard unary and binary functions
65     OP_ADD,  OP_SUB,  OP_MUL,  OP_DIV,
66     OP_NEG,  OP_EXP,  OP_LOG,  OP_POW, OP_CONSTPOW,
67     OP_SQRT,  OP_SQ,  OP_TWICE,
68     OP_SIN,  OP_COS,  OP_TAN,
69     OP_ASIN,  OP_ACOS,  OP_ATAN,
70     OP_LT, OP_LE, OP_EQ, OP_NE, OP_NOT, OP_AND, OP_OR,
71     OP_FLOOR,  OP_CEIL,  OP_FMOD, OP_FABS, OP_SIGN, OP_COPYSIGN, OP_IF_ELSE_ZERO,
72     OP_ERF,  OP_FMIN,  OP_FMAX,
73     OP_INV,
74     OP_SINH,  OP_COSH,  OP_TANH,
75     OP_ASINH, OP_ACOSH, OP_ATANH,
76     OP_ATAN2,
77 
78     // Double constant
79     OP_CONST,
80 
81     // Function input and output
82     OP_INPUT, OP_OUTPUT,
83 
84     // Free parameter
85     OP_PARAMETER,
86 
87     // Embedded function call
88     OP_CALL,
89 
90     // Find first nonzero in a vector
91     OP_FIND,
92 
93     // Find first nonzero in a vector
94     OP_LOW,
95 
96     // Embedded function call in parallel
97     OP_MAP,
98 
99     // Matrix multiplication
100     OP_MTIMES,
101 
102     // Solve linear system of equations
103     OP_SOLVE,
104 
105     // Matrix transpose
106     OP_TRANSPOSE,
107 
108     // Matrix determinant
109     OP_DETERMINANT,
110 
111     // Matrix inverse
112     OP_INVERSE,
113 
114     // Inner product
115     OP_DOT,
116 
117     // Bilinear form
118     OP_BILIN,
119 
120     // Rank-1 update
121     OP_RANK1,
122 
123     // Horizontal concatenation
124     OP_HORZCAT,
125 
126     // Vertical concatenation of vectors
127     OP_VERTCAT,
128 
129     // Diagonal concatenation
130     OP_DIAGCAT,
131 
132     // Horizontal split
133     OP_HORZSPLIT,
134 
135     // Vertical split of vectors
136     OP_VERTSPLIT,
137 
138     // Diagonal split
139     OP_DIAGSPLIT,
140 
141     // Reshape an expression
142     OP_RESHAPE,
143 
144     // Submatrix reference
145     OP_SUBREF,
146 
147     // Submatrix assignment
148     OP_SUBASSIGN,
149 
150     // Nonzero reference
151     OP_GETNONZEROS,
152 
153     // Parameteric nonzero reference
154     OP_GETNONZEROS_PARAM,
155 
156     // Nonzero addition
157     OP_ADDNONZEROS,
158 
159     // parametric nonzero addition
160     OP_ADDNONZEROS_PARAM,
161 
162     // Nonzero assignment
163     OP_SETNONZEROS,
164 
165     // Parametric nonzero assignment
166     OP_SETNONZEROS_PARAM,
167 
168     // Set sparse
169     OP_PROJECT,
170 
171     // Assertion
172     OP_ASSERTION,
173 
174     // Monitor
175     OP_MONITOR,
176 
177     // Norms
178     OP_NORM2, OP_NORM1, OP_NORMINF, OP_NORMF,
179 
180     // min/max
181     OP_MMIN, OP_MMAX,
182 
183     // Horizontal repeat
184     OP_HORZREPMAT,
185 
186     // Horizontal repeat sum
187     OP_HORZREPSUM,
188 
189     OP_ERFINV,
190     OP_PRINTME,
191     OP_LIFT,
192 
193     OP_EINSTEIN,
194 
195     OP_BSPLINE,
196 
197     OP_CONVEXIFY,
198   };
199   #define NUM_BUILT_IN_OPS (OP_CONVEXIFY+1)
200 
201   #define OP_
202 
203 #ifndef SWIG
204 
205   ///@{
206   /** \brief Enable using elementary numerical operations without std:: prefix */
207   using std::isfinite;
208   using std::sqrt;
209   using std::sin;
210   using std::cos;
211   using std::tan;
212   using std::atan;
213   using std::asin;
214   using std::acos;
215   using std::sinh;
216   using std::cosh;
217   using std::tanh;
218   using std::exp;
219   using std::log;
220   using std::log10;
221   using std::abs;
222   using std::fabs;
223   using std::floor;
224   using std::ceil;
225   using std::pow;
226   using std::fmod;
227   using std::atan2;
228   using std::erf;
229   using std::fmin;
230   using std::fmax;
231   using std::fabs;
232   using std::atanh;
233   using std::asinh;
234   using std::acosh;
235   using std::isnan;
236   using std::isinf;
237   ///@}
238 
239   ///@{
240   // Implement "missing" operations
241 
242   /// Sign function, note that sign(nan) == nan
sign(double x)243   inline double sign(double x) { return x<0 ? -1 : x>0 ? 1 : x;}
244   ///@}
245 
246   ///@}
247 
248   ///@{
249   /** \brief  CasADi additions */
simplify(double x)250   inline double simplify(double x) { return x;}
constpow(double x,double y)251   inline double constpow(double x, double y) { return pow(x, y);}
printme(double x,double y)252   inline double printme(double x, double y) {
253     std::ios::fmtflags f(uout().flags());
254     uout() << "|> " << y << " : ";
255     uout() << std::setprecision(std::numeric_limits<double>::digits10 + 1) << std::scientific;
256     uout() << x << std::endl;
257     uout().flags(f);
258     return x;
259   }
is_equal(double x,double y,casadi_int depth=0)260   inline bool is_equal(double x, double y, casadi_int depth=0) { return x==y;}
261 
262   #ifdef HAS_COPYSIGN
263   using std::copysign;
264   #else
265   /// copysign function
copysign(double x,double y)266   inline double copysign(double x, double y) { return y>=0 ? fabs(x) : -fabs(x);}
267   #endif //HAS_COPYSIGN
268 
269   // Integer maximum and minimum
casadi_max(casadi_int x,casadi_int y)270   inline casadi_int casadi_max(casadi_int x, casadi_int y) { return std::max(x, y);}
casadi_min(casadi_int x,casadi_int y)271   inline casadi_int casadi_min(casadi_int x, casadi_int y) { return std::min(x, y);}
272 
273   /// Conditional assignment
if_else_zero(double x,double y)274   inline double if_else_zero(double x, double y) { return x==0 ? 0 : y;}
if_else(double x,double y,double z)275   inline double if_else(double x, double y, double z) { return x==0 ? z : y;}
276 #ifdef HAS_ERFINV
277   using ::erfinv;
278 #else // HAS ERFINV
erfinv(double x)279   inline double erfinv(double x) throw() {
280     // Approximation found in Sourceforge and modified: Not very efficient
281     if (x>=1) {
282       return x==1 ? inf : nan;
283     } else if (x<=-1) {
284       return x==-1 ? -inf : nan;
285     } else if (x<-0.7) {
286       double z = sqrt(-log((1.0+x)/2.0));
287       return -(((1.641345311*z+3.429567803)*z-1.624906493)*z-1.970840454)/
288           ((1.637067800*z+3.543889200)*z+1.0);
289     } else {
290       double y;
291       if (x<0.7) {
292         double z = x*x;
293         y = x*(((-0.140543331*z+0.914624893)*z-1.645349621)*z+0.886226899)/
294             ((((-0.329097515*z+0.012229801)*z+1.442710462)*z-2.118377725)*z+1.0);
295       } else {
296         double z = sqrt(-log((1.0-x)/2.0));
297         y = (((1.641345311*z+3.429567803)*z-1.624906493)*z-1.970840454)/
298             ((1.637067800*z+3.543889200)*z+1.0);
299       }
300 
301       //polish x to full accuracy
302       y = y - (erf(y) - x) / (2.0/sqrt(pi) * exp(-y*y));
303       y = y - (erf(y) - x) / (2.0/sqrt(pi) * exp(-y*y));
304       return y;
305     }
306   }
307 #endif // HAS_ERFINV
308   ///@}
309 
310   template<typename T>
twice(const T & x)311   T twice(const T& x) {
312     return x+x;
313   }
314 
315   template<typename T>
sq(const T & x)316   T sq(const T& x) {
317     return x*x;
318   }
319 
320   template<casadi_int I>
321   struct UnaryOperation {
322     /// Function evaluation
323     template<typename T> static inline void fcn(const T& x, T& f);
324 
325     /// Partial derivatives
326     template<typename T> static inline void der(const T& x, const T& f, T* d);
327   };
328 
329   template<casadi_int I>
330   struct BinaryOperation {
331     /// Function evaluation
fcncasadi::BinaryOperation332     template<typename T> static inline void fcn(const T& x, const T& y, T& f) {
333         UnaryOperation<I>::fcn(x, f);}
334 
335     /// Partial derivatives - binary function
dercasadi::BinaryOperation336     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
337         UnaryOperation<I>::der(x, f, d); d[1]=0; }
338   };
339 
340   template<casadi_int I>
341   struct BinaryOperationE {
342     /// Function evaluation
fcncasadi::BinaryOperationE343     template<typename T> static inline T fcn(const T& x, const T& y) {
344       T ret;
345       BinaryOperation<I>::fcn(x, y, ret);
346       return ret;
347     }
348   };
349 
350   /// Calculate function and derivative
351   template<casadi_int I>
352   struct DerBinaryOpertion {
353     /// Perform the operation
derfcasadi::DerBinaryOpertion354     template<typename T> static inline void derf(const T& x, const T& y, T& f, T* d) {
355 
356       /** First save to temp since f might have the same address as x or y,
357       * in which case it will be incorrect in the second call
358       */
359       T tmp;
360 
361       /// Evaluate the function
362       BinaryOperation<I>::fcn(x, y, tmp);
363 
364       /// Evaluate the partial derivatives
365       BinaryOperation<I>::der(x, y, tmp, d);
366 
367       /// Now save f
368       f = tmp;
369     }
370   };
371 
372   /// Perform a binary operation on two scalars
373   template<casadi_int I>
374   struct BinaryOperationSS {
375     /// Function evaluation
fcncasadi::BinaryOperationSS376     template<typename T> static inline void fcn(const T& x, const T& y, T& f, casadi_int n) {
377       BinaryOperation<I>::fcn(x, y, f);
378     }
379 
380     /// Partial derivatives - binary function
dercasadi::BinaryOperationSS381     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d,
382         casadi_int n) {
383       BinaryOperation<I>::der(x, y, f, d);
384     }
385   };
386 
387 
388   /// Perform a binary operation on two vectors
389   template<casadi_int I>
390   struct BinaryOperationVV {
391     /// Function evaluation
fcncasadi::BinaryOperationVV392     template<typename T> static inline void fcn(const T* x, const T* y, T* f, casadi_int n) {
393       for (casadi_int i=0; i<n; ++i) {
394         BinaryOperation<I>::fcn(*x++, *y++, *f++);
395       }
396     }
397 
398     /// Partial derivatives - binary function
dercasadi::BinaryOperationVV399     template<typename T> static inline void der(const T* x, const T* y,
400         const T* f, T* d, casadi_int n) {
401       for (casadi_int i=0; i<n; ++i, d+=2) {
402         BinaryOperation<I>::der(*x++, *y++, *f++, d);
403       }
404     }
405   };
406 
407   /// Perform a binary operation on a vector and a scalar
408   template<casadi_int I>
409   struct BinaryOperationVS {
410     /// Function evaluation
fcncasadi::BinaryOperationVS411     template<typename T> static inline void fcn(const T* x, const T& y, T* f, casadi_int n) {
412       for (casadi_int i=0; i<n; ++i) {
413         BinaryOperation<I>::fcn(*x++, y, *f++);
414       }
415     }
416 
417     /// Partial derivatives - binary function
dercasadi::BinaryOperationVS418     template<typename T> static inline void der(const T* x, const T& y,
419         const T* f, T* d, casadi_int n) {
420       for (casadi_int i=0; i<n; ++i, d+=2) {
421         BinaryOperation<I>::der(*x++, y, *f++, d);
422       }
423     }
424   };
425 
426   /// Perform a binary operation on a scalar and a vector
427   template<casadi_int I>
428   struct BinaryOperationSV {
429     /// Function evaluation
fcncasadi::BinaryOperationSV430     template<typename T> static inline void fcn(const T& x, const T* y, T* f, casadi_int n) {
431       for (casadi_int i=0; i<n; ++i) {
432         BinaryOperation<I>::fcn(x, *y++, *f++);
433       }
434     }
435 
436     /// Partial derivatives - binary function
dercasadi::BinaryOperationSV437     template<typename T> static inline void der(const T& x, const T* y,
438         const T* f, T* d, casadi_int n) {
439       for (casadi_int i=0; i<n; ++i, d+=2) {
440         BinaryOperation<I>::der(x, *y++, *f++, d);
441       }
442     }
443   };
444 
445   ///@{
446   /// Smoothness (by default true)
447   template<casadi_int I> struct SmoothChecker { static const bool check=true;};
448   template<>      struct SmoothChecker<OP_LT>{ static const bool check=false;};
449   template<>      struct SmoothChecker<OP_LE>{ static const bool check=false;};
450   template<>      struct SmoothChecker<OP_FLOOR>{ static const bool check=false;};
451   template<>      struct SmoothChecker<OP_CEIL>{ static const bool check=false;};
452   template<>      struct SmoothChecker<OP_FMOD>{ static const bool check=false;};
453   template<>      struct SmoothChecker<OP_EQ>{ static const bool check=false;};
454   template<>      struct SmoothChecker<OP_NE>{ static const bool check=false;};
455   template<>      struct SmoothChecker<OP_SIGN>{ static const bool check=false;};
456   template<>      struct SmoothChecker<OP_COPYSIGN>{ static const bool check=false;};
457   template<>      struct SmoothChecker<OP_NOT>{ static const bool check=false;};
458   template<>      struct SmoothChecker<OP_AND>{ static const bool check=false;};
459   template<>      struct SmoothChecker<OP_OR>{ static const bool check=false;};
460   template<>      struct SmoothChecker<OP_IF_ELSE_ZERO>{ static const bool check=false;};
461   ///@}
462 
463   ///@{
464   /// If evaluated with the first argument zero, is the result zero?
465   template<casadi_int I> struct F0XChecker { static const bool check=false;};
466   template<>      struct F0XChecker<OP_ASSIGN>{ static const bool check=true;};
467   template<>      struct F0XChecker<OP_MUL>{ static const bool check=true;};
468   template<>      struct F0XChecker<OP_DIV>{ static const bool check=true;};
469   template<>      struct F0XChecker<OP_NEG>{ static const bool check=true;};
470   template<>      struct F0XChecker<OP_POW>{ static const bool check=true;};
471   template<>      struct F0XChecker<OP_CONSTPOW>{ static const bool check=true;};
472   template<>      struct F0XChecker<OP_SQRT>{ static const bool check=true;};
473   template<>      struct F0XChecker<OP_SQ>{ static const bool check=true;};
474   template<>      struct F0XChecker<OP_TWICE>{ static const bool check=true;};
475   template<>      struct F0XChecker<OP_SIN>{ static const bool check=true;};
476   template<>      struct F0XChecker<OP_TAN>{ static const bool check=true;};
477   template<>      struct F0XChecker<OP_ATAN>{ static const bool check=true;};
478   template<>      struct F0XChecker<OP_ASIN>{ static const bool check=true;};
479   template<>      struct F0XChecker<OP_FLOOR>{ static const bool check=true;};
480   template<>      struct F0XChecker<OP_CEIL>{ static const bool check=true;};
481   template<>      struct F0XChecker<OP_FMOD>{ static const bool check=true;};
482   template<>      struct F0XChecker<OP_FABS>{ static const bool check=true;};
483   template<>      struct F0XChecker<OP_SIGN>{ static const bool check=true;};
484   template<>      struct F0XChecker<OP_COPYSIGN>{ static const bool check=true;};
485   template<>      struct F0XChecker<OP_ERF>{ static const bool check=true;};
486   template<>      struct F0XChecker<OP_SINH>{ static const bool check=true;};
487   template<>      struct F0XChecker<OP_TANH>{ static const bool check=true;};
488   template<>      struct F0XChecker<OP_ASINH>{ static const bool check=true;};
489   template<>      struct F0XChecker<OP_ATANH>{ static const bool check=true;};
490   template<>      struct F0XChecker<OP_ERFINV>{ static const bool check=true;};
491   template<>      struct F0XChecker<OP_AND>{ static const bool check=true;};
492   template<>      struct F0XChecker<OP_IF_ELSE_ZERO>{ static const bool check=true;};
493   ///@}
494 
495   ///@{
496   /// If evaluated with the second argument zero, is the result zero?
497   template<casadi_int I> struct FX0Checker { static const bool check=false;};
498   template<>      struct FX0Checker<OP_MUL>{ static const bool check=true;};
499   template<>      struct FX0Checker<OP_AND>{ static const bool check=true;};
500   template<>      struct FX0Checker<OP_IF_ELSE_ZERO>{ static const bool check=true;};
501   ///@}
502 
503   ///@{
504   /// If evaluated with both arguments zero, is the result zero?
505   template<casadi_int I> struct F00Checker {
506     static const bool check=F0XChecker<I>::check || FX0Checker<I>::check;
507   };
508   template<>      struct F00Checker<OP_ADD>{ static const bool check=true;};
509   template<>      struct F00Checker<OP_SUB>{ static const bool check=true;};
510   template<>      struct F00Checker<OP_FMIN>{ static const bool check=true;};
511   template<>      struct F00Checker<OP_FMAX>{ static const bool check=true;};
512   template<>      struct F00Checker<OP_AND>{ static const bool check=true;};
513   template<>      struct F00Checker<OP_OR>{ static const bool check=true;};
514   template<>      struct F00Checker<OP_COPYSIGN>{ static const bool check=true;};
515   template<>      struct F00Checker<OP_LT>{ static const bool check=true;};
516   ///@}
517 
518   ///@{
519   /// Is commutative
520   template<casadi_int I> struct CommChecker { static const bool check=false;};
521   template<>      struct CommChecker<OP_ADD>{ static const bool check=true;};
522   template<>      struct CommChecker<OP_MUL>{ static const bool check=true;};
523   template<>      struct CommChecker<OP_EQ>{ static const bool check=true;};
524   template<>      struct CommChecker<OP_NE>{ static const bool check=true;};
525   template<>      struct CommChecker<OP_AND>{ static const bool check=true;};
526   template<>      struct CommChecker<OP_OR>{ static const bool check=true;};
527   ///@}
528 
529   ///@{
530   /// Always non-negative (false by default)
531   template<casadi_int I> struct NonnegativeChecker { static const bool check=false;};
532   template<>      struct NonnegativeChecker<OP_SQRT>{ static const bool check=true;};
533   template<>      struct NonnegativeChecker<OP_SQ>{ static const bool check=true;};
534   template<>      struct NonnegativeChecker<OP_EXP>{ static const bool check=true;};
535   template<>      struct NonnegativeChecker<OP_LT>{ static const bool check=true;};
536   template<>      struct NonnegativeChecker<OP_LE>{ static const bool check=true;};
537   template<>      struct NonnegativeChecker<OP_EQ>{ static const bool check=true;};
538   template<>      struct NonnegativeChecker<OP_NE>{ static const bool check=true;};
539   template<>      struct NonnegativeChecker<OP_NOT>{ static const bool check=true;};
540   template<>      struct NonnegativeChecker<OP_AND>{ static const bool check=true;};
541   template<>      struct NonnegativeChecker<OP_OR>{ static const bool check=true;};
542   ///@}
543 
544   ///@{
545   /// Is the operation binary as opposed to unary
546   template<casadi_int I> struct NargChecker { static const casadi_int check=1;};
547   template<>      struct NargChecker<OP_ADD>{ static const casadi_int check=2;};
548   template<>      struct NargChecker<OP_SUB>{ static const casadi_int check=2;};
549   template<>      struct NargChecker<OP_MUL>{ static const casadi_int check=2;};
550   template<>      struct NargChecker<OP_DIV>{ static const casadi_int check=2;};
551   template<>      struct NargChecker<OP_POW>{ static const casadi_int check=2;};
552   template<>      struct NargChecker<OP_CONSTPOW>{ static const casadi_int check=2;};
553   template<>      struct NargChecker<OP_EQ>{ static const casadi_int check=2;};
554   template<>      struct NargChecker<OP_NE>{ static const casadi_int check=2;};
555   template<>      struct NargChecker<OP_AND>{ static const casadi_int check=2;};
556   template<>      struct NargChecker<OP_OR>{ static const casadi_int check=2;};
557   template<>      struct NargChecker<OP_FMIN>{ static const casadi_int check=2;};
558   template<>      struct NargChecker<OP_FMAX>{ static const casadi_int check=2;};
559   template<>      struct NargChecker<OP_PRINTME>{ static const casadi_int check=2;};
560   template<>      struct NargChecker<OP_ATAN2>{ static const casadi_int check=2;};
561   template<>      struct NargChecker<OP_IF_ELSE_ZERO>{ static const casadi_int check=2;};
562   template<>      struct NargChecker<OP_FMOD>{ static const casadi_int check=2;};
563   template<>      struct NargChecker<OP_COPYSIGN>{ static const casadi_int check=2;};
564   template<>      struct NargChecker<OP_CONST>{ static const casadi_int check=0;};
565   template<>      struct NargChecker<OP_PARAMETER>{ static const casadi_int check=0;};
566   template<>      struct NargChecker<OP_INPUT>{ static const casadi_int check=0;};
567   ///@}
568 
569   /// Simple assignment
570   template<>
571   struct UnaryOperation<OP_ASSIGN>{
572   public:
fcncasadi::UnaryOperation573     template<typename T> static inline void fcn(const T& x, T& f) { f = x;}
dercasadi::UnaryOperation574     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 1; }
575   };
576 
577   /// Addition
578   template<>
579   struct BinaryOperation<OP_ADD>{
580   public:
fcncasadi::BinaryOperation581     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x+y;}
dercasadi::BinaryOperation582     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
583         d[0]=d[1]=1;}
584   };
585 
586   /// Subtraction
587   template<>
588   struct BinaryOperation<OP_SUB>{
589   public:
fcncasadi::BinaryOperation590     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x-y;}
dercasadi::BinaryOperation591     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
592         d[0]=1; d[1]=-1;}
593   };
594 
595   /// Multiplication
596   template<>
597   struct BinaryOperation<OP_MUL>{
598   public:
fcncasadi::BinaryOperation599     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x*y;}
dercasadi::BinaryOperation600     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
601         d[0]=y; d[1]=x;}
602   };
603 
604   /// Division
605   template<>
606   struct BinaryOperation<OP_DIV>{
607   public:
fcncasadi::BinaryOperation608     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x/y;}
dercasadi::BinaryOperation609     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
610         d[0]=1/y; d[1]=-f/y;}
611   };
612 
613   /// Negation
614   template<>
615   struct UnaryOperation<OP_NEG>{
616   public:
fcncasadi::UnaryOperation617     template<typename T> static inline void fcn(const T& x, T& f) { f = -x;}
dercasadi::UnaryOperation618     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=-1;}
619   };
620 
621   /// Natural exponent
622   template<>
623   struct UnaryOperation<OP_EXP>{
624   public:
fcncasadi::UnaryOperation625     template<typename T> static inline void fcn(const T& x, T& f) { f = exp(x);}
dercasadi::UnaryOperation626     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=f;}
627   };
628 
629   /// Natural logarithm
630   template<>
631   struct UnaryOperation<OP_LOG>{
632   public:
fcncasadi::UnaryOperation633     template<typename T> static inline void fcn(const T& x, T& f) { f = log(x);}
dercasadi::UnaryOperation634     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=1/x;}
635   };
636 
637   /// Power, defined only for x>=0
638   template<>
639   struct BinaryOperation<OP_POW>{
640   public:
fcncasadi::BinaryOperation641     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = pow(x, y);}
642     // See issue #104 why d[0] is no longer y*f/x
dercasadi::BinaryOperation643     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
644         d[0]=y*pow(x, y-1); d[1]=log(x)*f;}
645   };
646 
647   /// Power, defined only for y constant
648   template<>
649   struct BinaryOperation<OP_CONSTPOW>{
650   public:
fcncasadi::BinaryOperation651     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = pow(x, y);}
dercasadi::BinaryOperation652     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
653         d[0]=y*pow(x, y-1); d[1]=0;}
654   };
655 
656   /// Square root
657   template<>
658   struct UnaryOperation<OP_SQRT>{
659   public:
fcncasadi::UnaryOperation660     template<typename T> static inline void fcn(const T& x, T& f) { f = sqrt(x);}
dercasadi::UnaryOperation661     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=1/(twice(f));}
662   };
663 
664   /// Square
665   template<>
666   struct UnaryOperation<OP_SQ>{
667   public:
fcncasadi::UnaryOperation668     template<typename T> static inline void fcn(const T& x, T& f) { f = sq(x);}
dercasadi::UnaryOperation669     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=twice(x);}
670   };
671 
672   /// Times two
673   template<>
674   struct UnaryOperation<OP_TWICE>{
fcncasadi::UnaryOperation675     template<typename T> static inline void fcn(const T& x, T& f) { f = 2.*x;}
dercasadi::UnaryOperation676     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 2; }
677   };
678 
679   /// Sine
680   template<>
681   struct UnaryOperation<OP_SIN>{
682   public:
fcncasadi::UnaryOperation683     template<typename T> static inline void fcn(const T& x, T& f) { f = sin(x);}
dercasadi::UnaryOperation684     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=cos(x);}
685   };
686 
687   /// Cosine
688   template<>
689   struct UnaryOperation<OP_COS>{
690   public:
fcncasadi::UnaryOperation691     template<typename T> static inline void fcn(const T& x, T& f) { f = cos(x);}
dercasadi::UnaryOperation692     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=-sin(x);}
693   };
694 
695   /// Tangent
696   template<>
697   struct UnaryOperation<OP_TAN>{
698   public:
fcncasadi::UnaryOperation699     template<typename T> static inline void fcn(const T& x, T& f) { f = tan(x);}
dercasadi::UnaryOperation700     template<typename T> static inline void der(const T& x, const T& f, T* d)
701     { d[0] = 1/sq(cos(x));}
702   };
703 
704   /// Arcus sine
705   template<>
706   struct UnaryOperation<OP_ASIN>{
707   public:
fcncasadi::UnaryOperation708     template<typename T> static inline void fcn(const T& x, T& f) { f = asin(x);}
dercasadi::UnaryOperation709     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=1/sqrt(1-x*x);}
710   };
711 
712   /// Arcus cosine
713   template<>
714   struct UnaryOperation<OP_ACOS>{
715   public:
fcncasadi::UnaryOperation716     template<typename T> static inline void fcn(const T& x, T& f) { f = acos(x);}
dercasadi::UnaryOperation717     template<typename T> static inline void der(const T& x, const T& f, T* d)
718     { d[0]=-1/sqrt(1-x*x);}
719   };
720 
721   /// Arcus tangent
722   template<>
723   struct UnaryOperation<OP_ATAN>{
724   public:
fcncasadi::UnaryOperation725     template<typename T> static inline void fcn(const T& x, T& f) { f = atan(x);}
dercasadi::UnaryOperation726     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 1/(1+x*x);}
727   };
728 
729   /// Less than
730   template<>
731   struct BinaryOperation<OP_LT>{
732   public:
fcncasadi::BinaryOperation733     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x < y;}
dercasadi::BinaryOperation734     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
735         d[0]=d[1]=0;}
736   };
737 
738   /// Less or equal to
739   template<>
740   struct BinaryOperation<OP_LE>{
741   public:
fcncasadi::BinaryOperation742     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x <= y;}
dercasadi::BinaryOperation743     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
744         d[0]=d[1]=0;}
745   };
746 
747   /// Floor function
748   template<>
749   struct UnaryOperation<OP_FLOOR>{
750   public:
fcncasadi::UnaryOperation751     template<typename T> static inline void fcn(const T& x, T& f) { f = floor(x);}
dercasadi::UnaryOperation752     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 0;}
753   };
754 
755   /// Ceil function
756   template<>
757   struct UnaryOperation<OP_CEIL>{
758   public:
fcncasadi::UnaryOperation759     template<typename T> static inline void fcn(const T& x, T& f) { f = ceil(x);}
dercasadi::UnaryOperation760     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 0;}
761   };
762 
763   /// Remainder of division
764   template<>
765   struct BinaryOperation<OP_FMOD>{
fcncasadi::BinaryOperation766     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = fmod(x, y);}
dercasadi::BinaryOperation767     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
768       d[0]=1; d[1]=(f-x)/y;}
769   };
770 
771   /// Equal to
772   template<>
773   struct BinaryOperation<OP_EQ>{
fcncasadi::BinaryOperation774     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x==y;}
dercasadi::BinaryOperation775     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
776         d[0]=d[1]=0;}
777   };
778 
779   /// Not equal to
780   template<>
781   struct BinaryOperation<OP_NE>{
fcncasadi::BinaryOperation782     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x!=y;}
dercasadi::BinaryOperation783     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
784         d[0]=d[1]=0;}
785   };
786 
787   /// Logical not
788   template<>
789   struct UnaryOperation<OP_NOT>{
790   public:
fcncasadi::UnaryOperation791     template<typename T> static inline void fcn(const T& x, T& f) { f = !x;}
dercasadi::UnaryOperation792     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 0;}
793   };
794 
795   /// Logical and
796   template<>
797   struct BinaryOperation<OP_AND>{
fcncasadi::BinaryOperation798     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x && y;}
dercasadi::BinaryOperation799     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
800         d[0]=d[1]=0;}
801   };
802 
803   /// Logical or
804   template<>
805   struct BinaryOperation<OP_OR>{
fcncasadi::BinaryOperation806     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x || y;}
dercasadi::BinaryOperation807     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
808         d[0]=d[1]=0;}
809   };
810 
811   /// Error function
812   template<>
813   struct UnaryOperation<OP_ERF>{
fcncasadi::UnaryOperation814     template<typename T> static inline void fcn(const T& x, T& f) { f = erf(x);}
dercasadi::UnaryOperation815     template<typename T> static inline void der(const T& x, const T& f, T* d) {
816         d[0] = (2/sqrt(pi))*exp(-x*x);}
817   };
818 
819   /// Absolute value
820   template<>
821   struct UnaryOperation<OP_FABS>{
fcncasadi::UnaryOperation822     template<typename T> static inline void fcn(const T& x, T& f) { f = fabs(x);}
dercasadi::UnaryOperation823     template<typename T> static inline void der(const T& x, const T& f, T* d) {
824         d[0]=sign(x);}
825   };
826 
827   /// Sign
828   template<>
829   struct UnaryOperation<OP_SIGN>{
fcncasadi::UnaryOperation830     template<typename T> static inline void fcn(const T& x, T& f) { f = sign(x);}
dercasadi::UnaryOperation831     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=0;}
832   };
833 
834   /// Copysign
835   template<>
836   struct BinaryOperation<OP_COPYSIGN>{
fcncasadi::BinaryOperation837     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = copysign(x, y);}
dercasadi::BinaryOperation838     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
839         T e = 1; d[0]=copysign(e, y); d[1]=0;}
840   };
841 
842   /// Minimum
843   template<>
844   struct BinaryOperation<OP_FMIN>{
fcncasadi::BinaryOperation845     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = fmin(x, y);}
dercasadi::BinaryOperation846     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
847         d[0]=x<=y; d[1]=!d[0];}
848   };
849 
850   /// Maximum
851   template<>
852   struct BinaryOperation<OP_FMAX>{
fcncasadi::BinaryOperation853     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = fmax(x, y);}
dercasadi::BinaryOperation854     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
855         d[0]=x>=y; d[1]=!d[0];}
856   };
857 
858   /// Elementwise inverse
859   template<>
860   struct UnaryOperation<OP_INV>{
fcncasadi::UnaryOperation861     template<typename T> static inline void fcn(const T& x, T& f) { f = 1./x;}
dercasadi::UnaryOperation862     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = -f*f; }
863   };
864 
865   /// Hyperbolic sine
866   template<>
867   struct UnaryOperation<OP_SINH>{
fcncasadi::UnaryOperation868     template<typename T> static inline void fcn(const T& x, T& f) { f = sinh(x);}
dercasadi::UnaryOperation869     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = cosh(x); }
870   };
871 
872   /// Hyperbolic cosine
873   template<>
874   struct UnaryOperation<OP_COSH>{
fcncasadi::UnaryOperation875     template<typename T> static inline void fcn(const T& x, T& f) { f = cosh(x);}
dercasadi::UnaryOperation876     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = sinh(x); }
877   };
878 
879   /// Hyperbolic tangent
880   template<>
881   struct UnaryOperation<OP_TANH>{
fcncasadi::UnaryOperation882     template<typename T> static inline void fcn(const T& x, T& f) { f = tanh(x);}
dercasadi::UnaryOperation883     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 1-f*f; }
884   };
885 
886   /// Inverse hyperbolic sine
887   template<>
888   struct UnaryOperation<OP_ASINH>{
fcncasadi::UnaryOperation889     template<typename T> static inline void fcn(const T& x, T& f) { f = asinh(x);}
dercasadi::UnaryOperation890     template<typename T> static inline void der(const T& x, const T& f, T* d) {
891         d[0] = 1/sqrt(1+x*x); }
892   };
893 
894   /// Inverse hyperbolic cosine
895   template<>
896   struct UnaryOperation<OP_ACOSH>{
fcncasadi::UnaryOperation897     template<typename T> static inline void fcn(const T& x, T& f) { f = acosh(x);}
dercasadi::UnaryOperation898     template<typename T> static inline void der(const T& x, const T& f, T* d) {
899         d[0] = 1/sqrt(x-1)/sqrt(x+1); }
900   };
901 
902   /// Inverse hyperbolic tangent
903   template<>
904   struct UnaryOperation<OP_ATANH>{
fcncasadi::UnaryOperation905     template<typename T> static inline void fcn(const T& x, T& f) { f = atanh(x);}
dercasadi::UnaryOperation906     template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 1/(1-x*x); }
907   };
908 
909   /// Inverse of error function
910   template<>
911   struct UnaryOperation<OP_ERFINV>{
fcncasadi::UnaryOperation912     template<typename T> static inline void fcn(const T& x, T& f) { f = erfinv(x);}
dercasadi::UnaryOperation913     template<typename T> static inline void der(const T& x, const T& f, T* d) {
914         d[0] = (sqrt(pi)/2)*exp(f*f); }
915   };
916 
917   /// Identity operator with the side effect of printing
918   template<>
919   struct BinaryOperation<OP_PRINTME>{
fcncasadi::BinaryOperation920     template<typename T> static inline void fcn(const T& x, const T& y, T& f) {f = printme(x, y); }
dercasadi::BinaryOperation921     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
922         d[0]=1; d[1]=0;}
923   };
924 
925   /// Arctan2
926   template<>
927   struct BinaryOperation<OP_ATAN2>{
928   public:
fcncasadi::BinaryOperation929     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = atan2(x, y);}
dercasadi::BinaryOperation930     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
931         T t = x*x+y*y; d[0]=y/t; d[1]=-x/t;}
932   };
933 
934   /// Conditional assignment
935   template<>
936   struct BinaryOperation<OP_IF_ELSE_ZERO>{
937   public:
fcncasadi::BinaryOperation938     template<typename T> static inline void fcn(const T& x, const T& y, T& f) {
939         f = if_else_zero(x, y);}
dercasadi::BinaryOperation940     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
941         d[0]=0; d[1]=if_else_zero(x, T(1));}
942   };
943 
944   /// Inverse of error function
945   template<>
946   struct BinaryOperation<OP_LIFT>{
fcncasadi::BinaryOperation947     template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x;}
dercasadi::BinaryOperation948     template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
949         d[0] = 1; d[1] = 0; }
950   };
951 
952   template<template<casadi_int> class F, typename T>
operation_getter(casadi_int op)953   T operation_getter(casadi_int op) {
954     switch (static_cast<Operation>(op)) {
955     case OP_ASSIGN:        return F<OP_ASSIGN>::check;
956     case OP_ADD:           return F<OP_ADD>::check;
957     case OP_SUB:           return F<OP_SUB>::check;
958     case OP_MUL:           return F<OP_MUL>::check;
959     case OP_DIV:           return F<OP_DIV>::check;
960     case OP_NEG:           return F<OP_NEG>::check;
961     case OP_EXP:           return F<OP_EXP>::check;
962     case OP_LOG:           return F<OP_LOG>::check;
963     case OP_POW:           return F<OP_POW>::check;
964     case OP_CONSTPOW:      return F<OP_CONSTPOW>::check;
965     case OP_SQRT:          return F<OP_SQRT>::check;
966     case OP_SQ:            return F<OP_SQ>::check;
967     case OP_TWICE:         return F<OP_TWICE>::check;
968     case OP_SIN:           return F<OP_SIN>::check;
969     case OP_COS:           return F<OP_COS>::check;
970     case OP_TAN:           return F<OP_TAN>::check;
971     case OP_ASIN:          return F<OP_ASIN>::check;
972     case OP_ACOS:          return F<OP_ACOS>::check;
973     case OP_ATAN:          return F<OP_ATAN>::check;
974     case OP_LT:            return F<OP_LT>::check;
975     case OP_LE:            return F<OP_LE>::check;
976     case OP_EQ:            return F<OP_EQ>::check;
977     case OP_NE:            return F<OP_NE>::check;
978     case OP_NOT:           return F<OP_NOT>::check;
979     case OP_AND:           return F<OP_AND>::check;
980     case OP_OR:            return F<OP_OR>::check;
981     case OP_FLOOR:         return F<OP_FLOOR>::check;
982     case OP_CEIL:          return F<OP_CEIL>::check;
983     case OP_FMOD:          return F<OP_FMOD>::check;
984     case OP_FABS:          return F<OP_FABS>::check;
985     case OP_SIGN:          return F<OP_SIGN>::check;
986     case OP_COPYSIGN:      return F<OP_COPYSIGN>::check;
987     case OP_IF_ELSE_ZERO:  return F<OP_IF_ELSE_ZERO>::check;
988     case OP_ERF:           return F<OP_ERF>::check;
989     case OP_FMIN:          return F<OP_FMIN>::check;
990     case OP_FMAX:          return F<OP_FMAX>::check;
991     case OP_INV:           return F<OP_INV>::check;
992     case OP_SINH:          return F<OP_SINH>::check;
993     case OP_COSH:          return F<OP_COSH>::check;
994     case OP_TANH:          return F<OP_TANH>::check;
995     case OP_ASINH:         return F<OP_ASINH>::check;
996     case OP_ACOSH:         return F<OP_ACOSH>::check;
997     case OP_ATANH:         return F<OP_ATANH>::check;
998     case OP_ATAN2:         return F<OP_ATAN2>::check;
999     case OP_CONST:         return F<OP_CONST>::check;
1000     case OP_INPUT:         return F<OP_INPUT>::check;
1001     case OP_OUTPUT:        return F<OP_OUTPUT>::check;
1002     case OP_PARAMETER:     return F<OP_PARAMETER>::check;
1003     case OP_CALL:          return F<OP_CALL>::check;
1004     case OP_FIND:          return F<OP_FIND>::check;
1005     case OP_LOW:           return F<OP_LOW>::check;
1006     case OP_MAP:           return F<OP_MAP>::check;
1007     case OP_MTIMES:        return F<OP_MTIMES>::check;
1008     case OP_SOLVE:         return F<OP_SOLVE>::check;
1009     case OP_TRANSPOSE:     return F<OP_TRANSPOSE>::check;
1010     case OP_DETERMINANT:   return F<OP_DETERMINANT>::check;
1011     case OP_INVERSE:       return F<OP_INVERSE>::check;
1012     case OP_DOT:           return F<OP_DOT>::check;
1013     case OP_BILIN:         return F<OP_BILIN>::check;
1014     case OP_RANK1:         return F<OP_RANK1>::check;
1015     case OP_HORZCAT:       return F<OP_HORZCAT>::check;
1016     case OP_VERTCAT:       return F<OP_VERTCAT>::check;
1017     case OP_DIAGCAT:       return F<OP_DIAGCAT>::check;
1018     case OP_HORZSPLIT:     return F<OP_HORZSPLIT>::check;
1019     case OP_VERTSPLIT:     return F<OP_VERTSPLIT>::check;
1020     case OP_DIAGSPLIT:     return F<OP_DIAGSPLIT>::check;
1021     case OP_RESHAPE:       return F<OP_RESHAPE>::check;
1022     case OP_SUBREF:        return F<OP_SUBREF>::check;
1023     case OP_SUBASSIGN:     return F<OP_SUBASSIGN>::check;
1024     case OP_GETNONZEROS:   return F<OP_GETNONZEROS>::check;
1025     case OP_GETNONZEROS_PARAM:   return F<OP_GETNONZEROS_PARAM>::check;
1026     case OP_ADDNONZEROS:   return F<OP_ADDNONZEROS>::check;
1027     case OP_ADDNONZEROS_PARAM:   return F<OP_ADDNONZEROS>::check;
1028     case OP_SETNONZEROS:   return F<OP_SETNONZEROS>::check;
1029     case OP_SETNONZEROS_PARAM:   return F<OP_SETNONZEROS>::check;
1030     case OP_PROJECT:       return F<OP_PROJECT>::check;
1031     case OP_ASSERTION:     return F<OP_ASSERTION>::check;
1032     case OP_MONITOR:       return F<OP_MONITOR>::check;
1033     case OP_NORM2:         return F<OP_NORM2>::check;
1034     case OP_NORM1:         return F<OP_NORM1>::check;
1035     case OP_NORMINF:       return F<OP_NORMINF>::check;
1036     case OP_NORMF:         return F<OP_NORMF>::check;
1037     case OP_MMIN:          return F<OP_MMIN>::check;
1038     case OP_MMAX:          return F<OP_MMAX>::check;
1039     case OP_HORZREPMAT:    return F<OP_HORZREPMAT>::check;
1040     case OP_HORZREPSUM:    return F<OP_HORZREPSUM>::check;
1041     case OP_ERFINV:        return F<OP_ERFINV>::check;
1042     case OP_PRINTME:       return F<OP_PRINTME>::check;
1043     case OP_LIFT:          return F<OP_LIFT>::check;
1044     case OP_EINSTEIN:      return F<OP_EINSTEIN>::check;
1045     case OP_BSPLINE:       return F<OP_BSPLINE>::check;
1046     case OP_CONVEXIFY:     return F<OP_CONVEXIFY>::check;
1047     }
1048     return T();
1049   }
1050 
1051   template<template<casadi_int> class F>
operation_checker(casadi_int op)1052   bool operation_checker(casadi_int op) {
1053     return operation_getter<F, bool>(op);
1054   }
1055 
1056   /// Easy access to all the functions for a particular type
1057   template<typename T>
1058   struct casadi_math {
1059 
1060     /** \brief Evaluate a built in function (scalar-scalar) */
1061     static inline void fun(unsigned char op, const T& x, const T& y, T& f);
1062 
1063     /** \brief Evaluate a built in function (vector-vector) */
1064     static inline void fun(unsigned char op, const T* x, const T* y, T* f, casadi_int n);
1065 
1066     /** \brief Evaluate a built in function (vector-scalar) */
1067     static inline void fun(unsigned char op, const T* x, const T& y, T* f, casadi_int n);
1068 
1069     /** \brief Evaluate a built in function (scalar-vector) */
1070     static inline void fun(unsigned char op, const T& x, const T* y, T* f, casadi_int n);
1071 
1072     /** \brief Evaluate a built in derivative function */
1073     static inline void der(unsigned char op, const T& x, const T& y, const T& f, T* d);
1074 
1075     /** \brief Evaluate the function and the derivative function */
1076     static inline void derF(unsigned char op, const T& x, const T& y, T& f, T* d);
1077 
1078     /** \brief Is binary operation? */
1079     static inline bool is_binary(unsigned char op);
1080 
1081     /** \brief Is unary operation? */
1082     static inline bool is_unary(unsigned char op);
1083 
1084     /** \brief Number of dependencies */
1085     static inline casadi_int ndeps(unsigned char op);
1086 
1087     /** \brief Print */
1088     static inline std::string print(unsigned char op, const std::string& x,
1089                              const std::string& y);
1090     static inline std::string print(unsigned char op, const std::string& x);
1091     static inline std::string name(unsigned char op);
1092     static inline std::string pre(unsigned char op);
1093     static inline std::string sep(unsigned char op);
1094     static inline std::string post(unsigned char op);
1095   };
1096 
1097   /// Specialize the class so that it can be used with integer type
1098   template<>
1099   struct casadi_math<casadi_int>{
1100 
1101     /** \brief Evaluate a built in function */
funcasadi::casadi_math1102     static inline void fun(unsigned char op, const casadi_int& x,
1103         const casadi_int& y, casadi_int& f) {
1104       double ff(0);
1105       casadi_math<double>::fun(op, static_cast<double>(x), static_cast<double>(y), ff);
1106       f = static_cast<casadi_int>(ff);
1107     }
1108 
funcasadi::casadi_math1109     static inline void fun(unsigned char op, const casadi_int* x, const casadi_int* y,
1110         casadi_int* f, casadi_int n) {
1111       for (casadi_int i=0; i<n; ++i) {
1112         double ff(0);
1113         casadi_math<double>::fun(op, static_cast<double>(*x++), static_cast<double>(*y++), ff);
1114         *f++ = static_cast<casadi_int>(ff);
1115       }
1116     }
1117 
funcasadi::casadi_math1118     static inline void fun(unsigned char op, const casadi_int* x, const casadi_int& y,
1119         casadi_int* f, casadi_int n) {
1120       for (casadi_int i=0; i<n; ++i) {
1121         double ff;
1122         casadi_math<double>::fun(op, static_cast<double>(*x++), static_cast<double>(y), ff);
1123         *f++ = static_cast<casadi_int>(ff);
1124       }
1125     }
1126 
funcasadi::casadi_math1127     static inline void fun(unsigned char op, const casadi_int& x, const casadi_int* y,
1128         casadi_int* f, casadi_int n) {
1129       for (casadi_int i=0; i<n; ++i) {
1130         double ff;
1131         casadi_math<double>::fun(op, static_cast<double>(x), static_cast<double>(*y++), ff);
1132         *f++ = static_cast<casadi_int>(ff);
1133       }
1134     }
1135 
1136     /** \brief Evaluate a built in derivative function */
dercasadi::casadi_math1137     static inline void der(unsigned char op, const casadi_int& x, const casadi_int& y,
1138         const casadi_int& f, casadi_int* d) {
1139       double d_real[2] = {static_cast<double>(d[0]), static_cast<double>(d[1])};
1140       casadi_math<double>::der(op, static_cast<double>(x), static_cast<double>(y),
1141                                static_cast<double>(f), d_real);
1142       d[0] = static_cast<casadi_int>(d_real[0]);
1143       d[1] = static_cast<casadi_int>(d_real[1]);
1144     }
1145 
1146     /** \brief Evaluate the function and the derivative function */
derFcasadi::casadi_math1147     static inline void derF(unsigned char op, const casadi_int& x, const casadi_int& y,
1148         casadi_int& f, casadi_int* d) {
1149       double d_real[2] = {static_cast<double>(d[0]), static_cast<double>(d[1])};
1150       double f_real = static_cast<double>(f);
1151       casadi_math<double>::derF(op, static_cast<double>(x), static_cast<double>(y), f_real, d_real);
1152       f = static_cast<casadi_int>(f_real);
1153       d[0] = static_cast<casadi_int>(d_real[0]);
1154       d[1] = static_cast<casadi_int>(d_real[1]);
1155     }
1156 
1157     /** \brief Number of dependencies */
ndepscasadi::casadi_math1158     static inline casadi_int ndeps(unsigned char op) {
1159       return casadi_math<double>::ndeps(op);
1160     }
1161 
1162     /** \brief Print */
printcasadi::casadi_math1163     static inline std::string print(unsigned char op, const std::string& x,
1164                                     const std::string& y) {
1165       return casadi_math<double>::print(op, x, y);
1166     }
printcasadi::casadi_math1167     static inline std::string print(unsigned char op, const std::string& x) {
1168       return casadi_math<double>::print(op, x);
1169     }
precasadi::casadi_math1170     static inline std::string pre(unsigned char op) {
1171       return casadi_math<double>::pre(op);
1172     }
namecasadi::casadi_math1173     static inline std::string name(unsigned char op) {
1174       return casadi_math<double>::name(op);
1175     }
sepcasadi::casadi_math1176     static inline std::string sep(unsigned char op) {
1177       return casadi_math<double>::sep(op);
1178     }
postcasadi::casadi_math1179     static inline std::string post(unsigned char op) {
1180       return casadi_math<double>::post(op);
1181     }
1182   };
1183 
1184   // Template implementations
1185 
1186   template<typename T>
fun(unsigned char op,const T & x,const T & y,T & f)1187   inline void casadi_math<T>::fun(unsigned char op, const T& x, const T& y, T& f) {
1188     // NOTE: We define the implementation in a preprocessor macro to be able to force inlining,
1189     //  and to allow extensions in the VM
1190 #define CASADI_MATH_FUN_BUILTIN_GEN(CNAME, X, Y, F, N)                  \
1191     case OP_ASSIGN:    CNAME<OP_ASSIGN>::fcn(X, Y, F, N);        break; \
1192   case OP_ADD:       CNAME<OP_ADD>::fcn(X, Y, F, N);           break;   \
1193   case OP_SUB:       CNAME<OP_SUB>::fcn(X, Y, F, N);           break;   \
1194   case OP_MUL:       CNAME<OP_MUL>::fcn(X, Y, F, N);           break;   \
1195   case OP_DIV:       CNAME<OP_DIV>::fcn(X, Y, F, N);           break;   \
1196   case OP_NEG:       CNAME<OP_NEG>::fcn(X, Y, F, N);           break;   \
1197   case OP_EXP:       CNAME<OP_EXP>::fcn(X, Y, F, N);           break;   \
1198   case OP_LOG:       CNAME<OP_LOG>::fcn(X, Y, F, N);           break;   \
1199   case OP_POW:       CNAME<OP_POW>::fcn(X, Y, F, N);           break;   \
1200   case OP_CONSTPOW:  CNAME<OP_CONSTPOW>::fcn(X, Y, F, N);      break;   \
1201   case OP_SQRT:      CNAME<OP_SQRT>::fcn(X, Y, F, N);          break;   \
1202   case OP_SQ:        CNAME<OP_SQ>::fcn(X, Y, F, N);            break;   \
1203   case OP_TWICE:     CNAME<OP_TWICE>::fcn(X, Y, F, N);         break;   \
1204   case OP_SIN:       CNAME<OP_SIN>::fcn(X, Y, F, N);           break;   \
1205   case OP_COS:       CNAME<OP_COS>::fcn(X, Y, F, N);           break;   \
1206   case OP_TAN:       CNAME<OP_TAN>::fcn(X, Y, F, N);           break;   \
1207   case OP_ASIN:      CNAME<OP_ASIN>::fcn(X, Y, F, N);          break;   \
1208   case OP_ACOS:      CNAME<OP_ACOS>::fcn(X, Y, F, N);          break;   \
1209   case OP_ATAN:      CNAME<OP_ATAN>::fcn(X, Y, F, N);          break;   \
1210   case OP_LT:        CNAME<OP_LT>::fcn(X, Y, F, N);            break;   \
1211   case OP_LE:        CNAME<OP_LE>::fcn(X, Y, F, N);            break;   \
1212   case OP_EQ:        CNAME<OP_EQ>::fcn(X, Y, F, N);            break;   \
1213   case OP_NE:        CNAME<OP_NE>::fcn(X, Y, F, N);            break;   \
1214   case OP_NOT:       CNAME<OP_NOT>::fcn(X, Y, F, N);           break;   \
1215   case OP_AND:       CNAME<OP_AND>::fcn(X, Y, F, N);           break;   \
1216   case OP_OR:        CNAME<OP_OR>::fcn(X, Y, F, N);            break;   \
1217   case OP_IF_ELSE_ZERO: CNAME<OP_IF_ELSE_ZERO>::fcn(X, Y, F, N); break; \
1218   case OP_FLOOR:     CNAME<OP_FLOOR>::fcn(X, Y, F, N);         break;   \
1219   case OP_CEIL:      CNAME<OP_CEIL>::fcn(X, Y, F, N);          break;   \
1220   case OP_FMOD:      CNAME<OP_FMOD>::fcn(X, Y, F, N);          break;   \
1221   case OP_FABS:      CNAME<OP_FABS>::fcn(X, Y, F, N);          break;   \
1222   case OP_SIGN:      CNAME<OP_SIGN>::fcn(X, Y, F, N);          break;   \
1223   case OP_COPYSIGN:  CNAME<OP_COPYSIGN>::fcn(X, Y, F, N);      break;   \
1224   case OP_ERF:       CNAME<OP_ERF>::fcn(X, Y, F, N);           break;   \
1225   case OP_FMIN:      CNAME<OP_FMIN>::fcn(X, Y, F, N);          break;   \
1226   case OP_FMAX:      CNAME<OP_FMAX>::fcn(X, Y, F, N);          break;   \
1227   case OP_INV:       CNAME<OP_INV>::fcn(X, Y, F, N);           break;   \
1228   case OP_SINH:      CNAME<OP_SINH>::fcn(X, Y, F, N);          break;   \
1229   case OP_COSH:      CNAME<OP_COSH>::fcn(X, Y, F, N);          break;   \
1230   case OP_TANH:      CNAME<OP_TANH>::fcn(X, Y, F, N);          break;   \
1231   case OP_ASINH:     CNAME<OP_ASINH>::fcn(X, Y, F, N);         break;   \
1232   case OP_ACOSH:     CNAME<OP_ACOSH>::fcn(X, Y, F, N);         break;   \
1233   case OP_ATANH:     CNAME<OP_ATANH>::fcn(X, Y, F, N);         break;   \
1234   case OP_ATAN2:     CNAME<OP_ATAN2>::fcn(X, Y, F, N);         break;   \
1235   case OP_ERFINV:    CNAME<OP_ERFINV>::fcn(X, Y, F, N);        break;   \
1236   case OP_LIFT:      CNAME<OP_LIFT>::fcn(X, Y, F, N);          break;   \
1237   case OP_PRINTME:   CNAME<OP_PRINTME>::fcn(X, Y, F, N);       break;
1238 
1239 #define CASADI_MATH_FUN_BUILTIN(X, Y, F) CASADI_MATH_FUN_BUILTIN_GEN(BinaryOperationSS, X, Y, F, 1)
1240 
1241     switch (op) {
1242     CASADI_MATH_FUN_BUILTIN(x, y, f)
1243       }
1244   }
1245 
1246   template<typename T>
fun(unsigned char op,const T * x,const T * y,T * f,casadi_int n)1247   inline void casadi_math<T>::fun(unsigned char op, const T* x, const T* y, T* f, casadi_int n) {
1248     switch (op) {
1249       CASADI_MATH_FUN_BUILTIN_GEN(BinaryOperationVV, x, y, f, n)
1250         }
1251   }
1252 
1253   template<typename T>
fun(unsigned char op,const T * x,const T & y,T * f,casadi_int n)1254   inline void casadi_math<T>::fun(unsigned char op, const T* x, const T& y, T* f, casadi_int n) {
1255     switch (op) {
1256       CASADI_MATH_FUN_BUILTIN_GEN(BinaryOperationVS, x, y, f, n)
1257         }
1258   }
1259 
1260   template<typename T>
fun(unsigned char op,const T & x,const T * y,T * f,casadi_int n)1261   inline void casadi_math<T>::fun(unsigned char op, const T& x, const T* y, T* f, casadi_int n) {
1262     switch (op) {
1263       CASADI_MATH_FUN_BUILTIN_GEN(BinaryOperationSV, x, y, f, n)
1264         }
1265   }
1266 
1267 
1268   template<typename T>
der(unsigned char op,const T & x,const T & y,const T & f,T * d)1269   inline void casadi_math<T>::der(unsigned char op, const T& x, const T& y, const T& f, T* d) {
1270     // NOTE: We define the implementation in a preprocessor macro to be able to force inlining,
1271     // and to allow extensions in the VM
1272 #define CASADI_MATH_DER_BUILTIN(X, Y, F, D)                             \
1273     case OP_ASSIGN:    BinaryOperation<OP_ASSIGN>::der(X, Y, F, D);     break; \
1274   case OP_ADD:       BinaryOperation<OP_ADD>::der(X, Y, F, D);        break; \
1275   case OP_SUB:       BinaryOperation<OP_SUB>::der(X, Y, F, D);        break; \
1276   case OP_MUL:       BinaryOperation<OP_MUL>::der(X, Y, F, D);        break; \
1277   case OP_DIV:       BinaryOperation<OP_DIV>::der(X, Y, F, D);        break; \
1278   case OP_NEG:       BinaryOperation<OP_NEG>::der(X, Y, F, D);        break; \
1279   case OP_EXP:       BinaryOperation<OP_EXP>::der(X, Y, F, D);        break; \
1280   case OP_LOG:       BinaryOperation<OP_LOG>::der(X, Y, F, D);        break; \
1281   case OP_POW:       BinaryOperation<OP_POW>::der(X, Y, F, D);        break; \
1282   case OP_CONSTPOW:  BinaryOperation<OP_CONSTPOW>::der(X, Y, F, D);   break; \
1283   case OP_SQRT:      BinaryOperation<OP_SQRT>::der(X, Y, F, D);       break; \
1284   case OP_SQ:        BinaryOperation<OP_SQ>::der(X, Y, F, D);         break; \
1285   case OP_TWICE:     BinaryOperation<OP_TWICE>::der(X, Y, F, D);      break; \
1286   case OP_SIN:       BinaryOperation<OP_SIN>::der(X, Y, F, D);        break; \
1287   case OP_COS:       BinaryOperation<OP_COS>::der(X, Y, F, D);        break; \
1288   case OP_TAN:       BinaryOperation<OP_TAN>::der(X, Y, F, D);        break; \
1289   case OP_ASIN:      BinaryOperation<OP_ASIN>::der(X, Y, F, D);       break; \
1290   case OP_ACOS:      BinaryOperation<OP_ACOS>::der(X, Y, F, D);       break; \
1291   case OP_ATAN:      BinaryOperation<OP_ATAN>::der(X, Y, F, D);       break; \
1292   case OP_LT:        BinaryOperation<OP_LT>::der(X, Y, F, D);         break; \
1293   case OP_LE:        BinaryOperation<OP_LE>::der(X, Y, F, D);         break; \
1294   case OP_EQ:        BinaryOperation<OP_EQ>::der(X, Y, F, D);         break; \
1295   case OP_NE:        BinaryOperation<OP_NE>::der(X, Y, F, D);         break; \
1296   case OP_NOT:       BinaryOperation<OP_NOT>::der(X, Y, F, D);        break; \
1297   case OP_AND:       BinaryOperation<OP_AND>::der(X, Y, F, D);        break; \
1298   case OP_OR:        BinaryOperation<OP_OR>::der(X, Y, F, D);         break; \
1299   case OP_IF_ELSE_ZERO: BinaryOperation<OP_IF_ELSE_ZERO>::der(X, Y, F, D);         break; \
1300   case OP_FLOOR:     BinaryOperation<OP_FLOOR>::der(X, Y, F, D);      break; \
1301   case OP_CEIL:      BinaryOperation<OP_CEIL>::der(X, Y, F, D);       break; \
1302   case OP_FMOD:      BinaryOperation<OP_FMOD>::der(X, Y, F, D);       break; \
1303   case OP_FABS:      BinaryOperation<OP_FABS>::der(X, Y, F, D);       break; \
1304   case OP_SIGN:      BinaryOperation<OP_SIGN>::der(X, Y, F, D);       break; \
1305   case OP_COPYSIGN:  BinaryOperation<OP_COPYSIGN>::der(X, Y, F, D);   break; \
1306   case OP_ERF:       BinaryOperation<OP_ERF>::der(X, Y, F, D);        break; \
1307   case OP_FMIN:      BinaryOperation<OP_FMIN>::der(X, Y, F, D);       break; \
1308   case OP_FMAX:      BinaryOperation<OP_FMAX>::der(X, Y, F, D);       break; \
1309   case OP_INV:       BinaryOperation<OP_INV>::der(X, Y, F, D);        break; \
1310   case OP_SINH:      BinaryOperation<OP_SINH>::der(X, Y, F, D);       break; \
1311   case OP_COSH:      BinaryOperation<OP_COSH>::der(X, Y, F, D);       break; \
1312   case OP_TANH:      BinaryOperation<OP_TANH>::der(X, Y, F, D);       break; \
1313   case OP_ASINH:     BinaryOperation<OP_ASINH>::der(X, Y, F, D);      break; \
1314   case OP_ACOSH:     BinaryOperation<OP_ACOSH>::der(X, Y, F, D);      break; \
1315   case OP_ATANH:     BinaryOperation<OP_ATANH>::der(X, Y, F, D);      break; \
1316   case OP_ATAN2:     BinaryOperation<OP_ATAN2>::der(X, Y, F, D);      break; \
1317   case OP_ERFINV:    BinaryOperation<OP_ERFINV>::der(X, Y, F, D);     break; \
1318   case OP_LIFT:      BinaryOperation<OP_LIFT>::der(X, Y, F, D);       break; \
1319   case OP_PRINTME:   BinaryOperation<OP_PRINTME>::der(X, Y, F, D);    break;
1320 
1321     switch (op) {
1322     CASADI_MATH_DER_BUILTIN(x, y, f, d)
1323       }
1324   }
1325 
1326 
1327     template<typename T>
derF(unsigned char op,const T & x,const T & y,T & f,T * d)1328       inline void casadi_math<T>::derF(unsigned char op, const T& x, const T& y, T& f, T* d) {
1329     // NOTE: We define the implementation in a preprocessor macro to be able to force inlining,
1330     // and to allow extensions in the VM
1331 #define CASADI_MATH_DERF_BUILTIN(X, Y, F, D)                            \
1332     case OP_ASSIGN:    DerBinaryOpertion<OP_ASSIGN>::derf(X, Y, F, D);        break; \
1333   case OP_ADD:       DerBinaryOpertion<OP_ADD>::derf(X, Y, F, D);        break; \
1334   case OP_SUB:       DerBinaryOpertion<OP_SUB>::derf(X, Y, F, D);        break; \
1335   case OP_MUL:       DerBinaryOpertion<OP_MUL>::derf(X, Y, F, D);        break; \
1336   case OP_DIV:       DerBinaryOpertion<OP_DIV>::derf(X, Y, F, D);        break; \
1337   case OP_NEG:       DerBinaryOpertion<OP_NEG>::derf(X, Y, F, D);        break; \
1338   case OP_EXP:       DerBinaryOpertion<OP_EXP>::derf(X, Y, F, D);        break; \
1339   case OP_LOG:       DerBinaryOpertion<OP_LOG>::derf(X, Y, F, D);        break; \
1340   case OP_POW:       DerBinaryOpertion<OP_POW>::derf(X, Y, F, D);        break; \
1341   case OP_CONSTPOW:  DerBinaryOpertion<OP_CONSTPOW>::derf(X, Y, F, D);   break; \
1342   case OP_SQRT:      DerBinaryOpertion<OP_SQRT>::derf(X, Y, F, D);       break; \
1343   case OP_SQ:        DerBinaryOpertion<OP_SQ>::derf(X, Y, F, D);         break; \
1344   case OP_TWICE:     DerBinaryOpertion<OP_TWICE>::derf(X, Y, F, D);      break; \
1345   case OP_SIN:       DerBinaryOpertion<OP_SIN>::derf(X, Y, F, D);        break; \
1346   case OP_COS:       DerBinaryOpertion<OP_COS>::derf(X, Y, F, D);        break; \
1347   case OP_TAN:       DerBinaryOpertion<OP_TAN>::derf(X, Y, F, D);        break; \
1348   case OP_ASIN:      DerBinaryOpertion<OP_ASIN>::derf(X, Y, F, D);       break; \
1349   case OP_ACOS:      DerBinaryOpertion<OP_ACOS>::derf(X, Y, F, D);       break; \
1350   case OP_ATAN:      DerBinaryOpertion<OP_ATAN>::derf(X, Y, F, D);       break; \
1351   case OP_LT:        DerBinaryOpertion<OP_LT>::derf(X, Y, F, D);         break; \
1352   case OP_LE:        DerBinaryOpertion<OP_LE>::derf(X, Y, F, D);         break; \
1353   case OP_EQ:        DerBinaryOpertion<OP_EQ>::derf(X, Y, F, D);         break; \
1354   case OP_NE:        DerBinaryOpertion<OP_NE>::derf(X, Y, F, D);         break; \
1355   case OP_NOT:       DerBinaryOpertion<OP_NOT>::derf(X, Y, F, D);        break; \
1356   case OP_AND:       DerBinaryOpertion<OP_AND>::derf(X, Y, F, D);        break; \
1357   case OP_OR:        DerBinaryOpertion<OP_OR>::derf(X, Y, F, D);         break; \
1358   case OP_IF_ELSE_ZERO: DerBinaryOpertion<OP_IF_ELSE_ZERO>::derf(X, Y, F, D);         break; \
1359   case OP_FLOOR:     DerBinaryOpertion<OP_FLOOR>::derf(X, Y, F, D);      break; \
1360   case OP_CEIL:      DerBinaryOpertion<OP_CEIL>::derf(X, Y, F, D);       break; \
1361   case OP_FMOD:      DerBinaryOpertion<OP_FMOD>::derf(X, Y, F, D);       break; \
1362   case OP_FABS:      DerBinaryOpertion<OP_FABS>::derf(X, Y, F, D);       break; \
1363   case OP_SIGN:      DerBinaryOpertion<OP_SIGN>::derf(X, Y, F, D);       break; \
1364   case OP_COPYSIGN:  DerBinaryOpertion<OP_COPYSIGN>::derf(X, Y, F, D);   break; \
1365   case OP_ERF:       DerBinaryOpertion<OP_ERF>::derf(X, Y, F, D);        break; \
1366   case OP_FMIN:      DerBinaryOpertion<OP_FMIN>::derf(X, Y, F, D);       break; \
1367   case OP_FMAX:      DerBinaryOpertion<OP_FMAX>::derf(X, Y, F, D);       break; \
1368   case OP_INV:       DerBinaryOpertion<OP_INV>::derf(X, Y, F, D);        break; \
1369   case OP_SINH:      DerBinaryOpertion<OP_SINH>::derf(X, Y, F, D);       break; \
1370   case OP_COSH:      DerBinaryOpertion<OP_COSH>::derf(X, Y, F, D);       break; \
1371   case OP_TANH:      DerBinaryOpertion<OP_TANH>::derf(X, Y, F, D);       break; \
1372   case OP_ASINH:     DerBinaryOpertion<OP_ASINH>::derf(X, Y, F, D);      break; \
1373   case OP_ACOSH:     DerBinaryOpertion<OP_ACOSH>::derf(X, Y, F, D);      break; \
1374   case OP_ATANH:     DerBinaryOpertion<OP_ATANH>::derf(X, Y, F, D);      break; \
1375   case OP_ATAN2:     DerBinaryOpertion<OP_ATAN2>::derf(X, Y, F, D);      break; \
1376   case OP_ERFINV:    DerBinaryOpertion<OP_ERFINV>::derf(X, Y, F, D);     break; \
1377   case OP_LIFT:      DerBinaryOpertion<OP_LIFT>::derf(X, Y, F, D);       break; \
1378   case OP_PRINTME:   DerBinaryOpertion<OP_PRINTME>::derf(X, Y, F, D);    break;
1379 
1380     switch (op) {
1381       CASADI_MATH_DERF_BUILTIN(x, y, f, d)
1382         }
1383   }
1384 
1385   #define CASADI_MATH_BINARY_BUILTIN              \
1386     case OP_ADD:                                  \
1387     case OP_SUB:                                  \
1388     case OP_MUL:                                  \
1389     case OP_DIV:                                  \
1390     case OP_POW:                                  \
1391     case OP_CONSTPOW:                             \
1392     case OP_LT:                                   \
1393     case OP_LE:                                   \
1394     case OP_EQ:                                   \
1395     case OP_NE:                                   \
1396     case OP_AND:                                  \
1397     case OP_OR:                                   \
1398     case OP_IF_ELSE_ZERO:                         \
1399     case OP_COPYSIGN:                             \
1400     case OP_FMOD:                                 \
1401     case OP_FMIN:                                 \
1402     case OP_FMAX:                                 \
1403     case OP_ATAN2:                                \
1404     case OP_PRINTME:                              \
1405     case OP_LIFT:
1406 
1407   #define CASADI_MATH_UNARY_BUILTIN              \
1408     case OP_ASSIGN:                              \
1409     case OP_NEG:                                 \
1410     case OP_EXP:                                 \
1411     case OP_LOG:                                 \
1412     case OP_SQRT:                                \
1413     case OP_SQ:                                  \
1414     case OP_TWICE:                               \
1415     case OP_SIN:                                 \
1416     case OP_COS:                                 \
1417     case OP_TAN:                                 \
1418     case OP_ASIN:                                \
1419     case OP_ACOS:                                \
1420     case OP_ATAN:                                \
1421     case OP_FLOOR:                               \
1422     case OP_CEIL:                                \
1423     case OP_NOT:                                 \
1424     case OP_ERF:                                 \
1425     case OP_FABS:                                \
1426     case OP_SIGN:                                \
1427     case OP_INV:                                 \
1428     case OP_SINH:                                \
1429     case OP_COSH:                                \
1430     case OP_TANH:                                \
1431     case OP_ASINH:                               \
1432     case OP_ACOSH:                               \
1433     case OP_ATANH:                               \
1434     case OP_ERFINV:
1435 
1436   template<typename T>
is_binary(unsigned char op)1437   bool casadi_math<T>::is_binary(unsigned char op) {
1438     switch (op) {
1439       CASADI_MATH_BINARY_BUILTIN
1440       return true;
1441     default:
1442       return false;
1443     }
1444   }
1445 
1446   template<typename T>
is_unary(unsigned char op)1447   bool casadi_math<T>::is_unary(unsigned char op) {
1448     switch (op) {
1449       CASADI_MATH_UNARY_BUILTIN
1450       return true;
1451     default:
1452       return false;
1453     }
1454   }
1455 
1456   template<typename T>
ndeps(unsigned char op)1457   inline casadi_int casadi_math<T>::ndeps(unsigned char op) {
1458     switch (op) {
1459       case OP_CONST:
1460       case OP_PARAMETER:
1461       case OP_INPUT:
1462         return 0;
1463         CASADI_MATH_BINARY_BUILTIN
1464           return 2;
1465       default:
1466         return 1;
1467     }
1468   }
1469 
1470   template<typename T>
1471   inline std::string
print(unsigned char op,const std::string & x,const std::string & y)1472   casadi_math<T>::print(unsigned char op,
1473                         const std::string& x, const std::string& y) {
1474     casadi_assert_dev(ndeps(op)==2);
1475     return pre(op) + x + sep(op) + y + post(op);
1476   }
1477 
1478   template<typename T>
1479   inline std::string
print(unsigned char op,const std::string & x)1480   casadi_math<T>::print(unsigned char op, const std::string& x) {
1481     casadi_assert_dev(ndeps(op)==1);
1482     return pre(op) + x + post(op);
1483   }
1484 
1485   template<typename T>
name(unsigned char op)1486   inline std::string casadi_math<T>::name(unsigned char op) {
1487     switch (op) {
1488     case OP_ASSIGN:         return "assign";
1489     case OP_ADD:            return "add";
1490     case OP_SUB:            return "sub";
1491     case OP_MUL:            return "mul";
1492     case OP_DIV:            return "div";
1493     case OP_NEG:            return "neg";
1494     case OP_EXP:            return "exp";
1495     case OP_LOG:            return "log";
1496     case OP_CONSTPOW:
1497     case OP_POW:            return "pow";
1498     case OP_SQRT:           return "sqrt";
1499     case OP_SQ:             return "sq";
1500     case OP_TWICE:          return "twice";
1501     case OP_SIN:            return "sin";
1502     case OP_COS:            return "cos";
1503     case OP_TAN:            return "tan";
1504     case OP_ASIN:           return "asin";
1505     case OP_ACOS:           return "acos";
1506     case OP_ATAN:           return "atan";
1507     case OP_LT:             return "lt";
1508     case OP_LE:             return "le";
1509     case OP_EQ:             return "eq";
1510     case OP_NE:             return "ne";
1511     case OP_NOT:            return "not";
1512     case OP_AND:            return "and";
1513     case OP_OR:             return "or";
1514     case OP_FLOOR:          return "floor";
1515     case OP_CEIL:           return "ceil";
1516     case OP_FMOD:           return "fmod";
1517     case OP_FABS:           return "fabs";
1518     case OP_SIGN:           return "sign";
1519     case OP_COPYSIGN:       return "copysign";
1520     case OP_IF_ELSE_ZERO:   return "if_else_zero";
1521     case OP_ERF:            return "erf";
1522     case OP_FMIN:           return "fmin";
1523     case OP_FMAX:           return "fmax";
1524     case OP_INV:            return "inv";
1525     case OP_SINH:           return "sinh";
1526     case OP_COSH:           return "cosh";
1527     case OP_TANH:           return "tanh";
1528     case OP_ASINH:          return "asinh";
1529     case OP_ACOSH:          return "acosh";
1530     case OP_ATANH:          return "atanh";
1531     case OP_ATAN2:          return "atan2";
1532     case OP_CONST:          return "const";
1533     case OP_INPUT:          return "input";
1534     case OP_OUTPUT:         return "output";
1535     case OP_PARAMETER:      return "parameter";
1536     case OP_CALL:           return "call";
1537     case OP_MTIMES:         return "mtimes";
1538     case OP_SOLVE:          return "solve";
1539     case OP_TRANSPOSE:      return "transpose";
1540     case OP_DETERMINANT:    return "determinant";
1541     case OP_INVERSE:        return "inverse";
1542     case OP_DOT:            return "dot";
1543     case OP_HORZCAT:        return "horzcat";
1544     case OP_VERTCAT:        return "vertcat";
1545     case OP_DIAGCAT:        return "diagcat";
1546     case OP_HORZSPLIT:      return "horzsplit";
1547     case OP_VERTSPLIT:      return "vertsplit";
1548     case OP_DIAGSPLIT:      return "diagsplit";
1549     case OP_RESHAPE:        return "reshape";
1550     case OP_SUBREF:         return "subref";
1551     case OP_SUBASSIGN:      return "subassign";
1552     case OP_GETNONZEROS:    return "getnonzeros";
1553     case OP_GETNONZEROS_PARAM:    return "getnonzeros_param";
1554     case OP_ADDNONZEROS:    return "addnonzeros";
1555     case OP_ADDNONZEROS_PARAM:    return "addnonzeros_param";
1556     case OP_SETNONZEROS:    return "setnonzeros";
1557     case OP_SETNONZEROS_PARAM:    return "setnonzeros_param";
1558     case OP_PROJECT:        return "project";
1559     case OP_ASSERTION:      return "assertion";
1560     case OP_NORM2:          return "norm2";
1561     case OP_NORM1:          return "norm1";
1562     case OP_NORMINF:        return "norminf";
1563     case OP_NORMF:          return "normf";
1564     case OP_ERFINV:         return "erfinv";
1565     case OP_PRINTME:        return "printme";
1566     case OP_LIFT:           return "lift";
1567     case OP_EINSTEIN:       return "einstein";
1568     case OP_BSPLINE:        return "bspline";
1569     case OP_CONVEXIFY:      return "convexify";
1570     }
1571     return nullptr;
1572   }
1573 
1574   template<typename T>
pre(unsigned char op)1575   inline std::string casadi_math<T>::pre(unsigned char op) {
1576     switch (op) {
1577     case OP_ASSIGN:    return "";
1578     case OP_ADD:       return "(";
1579     case OP_SUB:       return "(";
1580     case OP_MUL:       return "(";
1581     case OP_DIV:       return "(";
1582     case OP_NEG:       return "(-";
1583     case OP_TWICE:     return "(2.*";
1584     case OP_LT:        return "(";
1585     case OP_LE:        return "(";
1586     case OP_EQ:        return "(";
1587     case OP_NE:        return "(";
1588     case OP_NOT:       return "(!";
1589     case OP_AND:       return "(";
1590     case OP_OR:        return "(";
1591     case OP_IF_ELSE_ZERO: return "(";
1592     case OP_INV:       return "(1./";
1593     default: return name(op) + "(";
1594     }
1595   }
1596 
1597   template<typename T>
sep(unsigned char op)1598   inline std::string casadi_math<T>::sep(unsigned char op) {
1599     switch (op) {
1600     case OP_ADD:       return "+";
1601     case OP_SUB:       return "-";
1602     case OP_MUL:       return "*";
1603     case OP_DIV:       return "/";
1604     case OP_LT:        return "<";
1605     case OP_LE:        return "<=";
1606     case OP_EQ:        return "==";
1607     case OP_NE:        return "!=";
1608     case OP_AND:       return "&&";
1609     case OP_OR:        return "||";
1610     case OP_IF_ELSE_ZERO: return "?";
1611     default:           return ",";
1612     }
1613   }
1614 
1615   template<typename T>
post(unsigned char op)1616   inline std::string casadi_math<T>::post(unsigned char op) {
1617     switch (op) {
1618     case OP_ASSIGN:       return "";
1619     case OP_IF_ELSE_ZERO: return ":0)";
1620     default:              return ")";
1621     }
1622   }
1623 
1624 #endif // SWIG
1625 
1626 } // namespace casadi
1627 
1628 /// \endcond
1629 
1630 #endif // CASADI_CALCULUS_HPP
1631