calculus.hpp
1 /*
2  * This file is part of CasADi.
3  *
4  * CasADi -- A symbolic framework for dynamic optimization.
5  * Copyright (C) 2010-2023 Joel Andersson, Joris Gillis, Moritz Diehl,
6  * KU 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 
39 
40 namespace casadi {
41 #ifndef SWIG
43 #ifdef M_PI
44  const double pi = M_PI;
45 #else
46  const double pi = 3.14159265358979323846;
47 #endif
48 
50  const double inf = std::numeric_limits<double>::infinity();
51 
53  const double nan = std::numeric_limits<double>::quiet_NaN();
54 
56  const double eps = std::numeric_limits<double>::epsilon();
57 #endif // SWIG
58 
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  // Parametric 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  // Sparsity cast
200  OP_SPARSITY_CAST,
201 
202  OP_LOG1P,
203 
204  OP_EXPM1,
205 
206  OP_HYPOT,
207 
208  OP_LOGSUMEXP,
209 
210  OP_REMAINDER,
211 
212  OP_DUMP,
213 
214  OP_KRON,
215 
216  OP_KRON_CONTRACT,
217 
218  };
219  #define NUM_BUILT_IN_OPS (OP_KRON_CONTRACT+1)
220 
221  #define OP_
222 
223 #ifndef SWIG
224 
226 
229  using std::isfinite;
230  using std::sqrt;
231  using std::sin;
232  using std::cos;
233  using std::tan;
234  using std::atan;
235  using std::asin;
236  using std::acos;
237  using std::sinh;
238  using std::cosh;
239  using std::tanh;
240  using std::exp;
241  using std::log;
242  using std::log10;
243  using std::abs;
244  using std::fabs;
245  using std::floor;
246  using std::ceil;
247  using std::pow;
248  using std::fmod;
249  using std::remainder;
250  using std::atan2;
251  using std::erf;
252  using std::fmin;
253  using std::fmax;
254  using std::fabs;
255  using std::atanh;
256  using std::asinh;
257  using std::acosh;
258  using std::isnan;
259  using std::isinf;
260  using std::log1p;
261  using std::expm1;
262  using std::hypot;
263  using std::copysign;
265 
267  // Implement "missing" operations
268 
270  inline double sign(double x) { return x<0 ? -1 : x>0 ? 1 : x;}
272 
274 
276 
277  inline double simplify(double x) { return x;}
278  inline double constpow(double x, double y) { return pow(x, y);}
279  inline double printme(double x, double y) {
280  std::ios::fmtflags f(uout().flags());
281  uout() << "|> " << y << " : ";
282  uout() << std::setprecision(std::numeric_limits<double>::digits10 + 1) << std::scientific;
283  uout() << x << std::endl;
284  uout().flags(f);
285  return x;
286  }
287  inline bool is_equal(double x, double y, casadi_int depth=0) { return x==y;}
288 
289 
290  // Integer maximum and minimum
291  inline casadi_int casadi_max(casadi_int x, casadi_int y) { return std::max(x, y);}
292  inline casadi_int casadi_min(casadi_int x, casadi_int y) { return std::min(x, y);}
293 
295  inline double if_else_zero(double x, double y) { return x==0 ? 0 : y;}
296  inline double if_else(double x, double y, double z) { return x==0 ? z : y;}
297 #ifdef HAS_ERFINV
298  using ::erfinv;
299 #else // HAS ERFINV
300  inline double erfinv(double x) throw() {
301  // Approximation found in Sourceforge and modified: Not very efficient
302  if (x>=1) {
303  return x==1 ? inf : nan;
304  } else if (x<=-1) {
305  return x==-1 ? -inf : nan;
306  } else if (x<-0.7) {
307  double z = sqrt(-log((1.0+x)/2.0));
308  return -(((1.641345311*z+3.429567803)*z-1.624906493)*z-1.970840454)/
309  ((1.637067800*z+3.543889200)*z+1.0);
310  } else {
311  double y;
312  if (x<0.7) {
313  double z = x*x;
314  y = x*(((-0.140543331*z+0.914624893)*z-1.645349621)*z+0.886226899)/
315  ((((-0.329097515*z+0.012229801)*z+1.442710462)*z-2.118377725)*z+1.0);
316  } else {
317  double z = sqrt(-log((1.0-x)/2.0));
318  y = (((1.641345311*z+3.429567803)*z-1.624906493)*z-1.970840454)/
319  ((1.637067800*z+3.543889200)*z+1.0);
320  }
321 
322  //polish x to full accuracy
323  y = y - (erf(y) - x) / (2.0/sqrt(pi) * exp(-y*y));
324  y = y - (erf(y) - x) / (2.0/sqrt(pi) * exp(-y*y));
325  return y;
326  }
327  }
328 #endif // HAS_ERFINV
330 
331  template<typename T>
332  T twice(const T& x) {
333  return x+x;
334  }
335 
336  template<typename T>
337  T sq(const T& x) {
338  return x*x;
339  }
340 
341  template<casadi_int I>
342  struct UnaryOperation {
344  template<typename T> static inline void fcn(const T& x, T& f);
345 
347  template<typename T> static inline void der(const T& x, const T& f, T* d);
348  };
349 
350  template<casadi_int I>
351  struct BinaryOperation {
353  template<typename T> static inline void fcn(const T& x, const T& y, T& f) {
354  UnaryOperation<I>::fcn(x, f);}
355 
357  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
358  UnaryOperation<I>::der(x, f, d); d[1]=0; }
359  };
360 
361  template<casadi_int I>
362  struct BinaryOperationE {
364  template<typename T> static inline T fcn(const T& x, const T& y) {
365  T ret;
366  BinaryOperation<I>::fcn(x, y, ret);
367  return ret;
368  }
369  };
370 
372  template<casadi_int I>
373  struct DerBinaryOperation {
375  template<typename T> static inline void derf(const T& x, const T& y, T& f, T* d) {
376 
380  T tmp;
381 
383  BinaryOperation<I>::fcn(x, y, tmp);
384 
386  BinaryOperation<I>::der(x, y, tmp, d);
387 
389  f = tmp;
390  }
391  };
392 
394  template<casadi_int I>
395  struct BinaryOperationSS {
397  template<typename T> static inline void fcn(const T& x, const T& y, T& f, casadi_int n) {
398  BinaryOperation<I>::fcn(x, y, f);
399  }
400 
402  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d,
403  casadi_int n) {
404  BinaryOperation<I>::der(x, y, f, d);
405  }
406  };
407 
408 
410  template<casadi_int I>
411  struct BinaryOperationVV {
413  template<typename T> static inline void fcn(const T* x, const T* y, T* f, casadi_int n) {
414  for (casadi_int i=0; i<n; ++i) {
415  BinaryOperation<I>::fcn(*x++, *y++, *f++);
416  }
417  }
418 
420  template<typename T> static inline void der(const T* x, const T* y,
421  const T* f, T* d, casadi_int n) {
422  for (casadi_int i=0; i<n; ++i, d+=2) {
423  BinaryOperation<I>::der(*x++, *y++, *f++, d);
424  }
425  }
426  };
427 
429  template<casadi_int I>
430  struct BinaryOperationVS {
432  template<typename T> static inline void fcn(const T* x, const T& y, T* f, casadi_int n) {
433  for (casadi_int i=0; i<n; ++i) {
434  BinaryOperation<I>::fcn(*x++, y, *f++);
435  }
436  }
437 
439  template<typename T> static inline void der(const T* x, const T& y,
440  const T* f, T* d, casadi_int n) {
441  for (casadi_int i=0; i<n; ++i, d+=2) {
442  BinaryOperation<I>::der(*x++, y, *f++, d);
443  }
444  }
445  };
446 
448  template<casadi_int I>
449  struct BinaryOperationSV {
451  template<typename T> static inline void fcn(const T& x, const T* y, T* f, casadi_int n) {
452  for (casadi_int i=0; i<n; ++i) {
453  BinaryOperation<I>::fcn(x, *y++, *f++);
454  }
455  }
456 
458  template<typename T> static inline void der(const T& x, const T* y,
459  const T* f, T* d, casadi_int n) {
460  for (casadi_int i=0; i<n; ++i, d+=2) {
461  BinaryOperation<I>::der(x, *y++, *f++, d);
462  }
463  }
464  };
465 
468  template<casadi_int I> struct SmoothChecker { static const bool check=true;};
469  template<> struct SmoothChecker<OP_LT>{ static const bool check=false;};
470  template<> struct SmoothChecker<OP_LE>{ static const bool check=false;};
471  template<> struct SmoothChecker<OP_FLOOR>{ static const bool check=false;};
472  template<> struct SmoothChecker<OP_CEIL>{ static const bool check=false;};
473  template<> struct SmoothChecker<OP_FMOD>{ static const bool check=false;};
474  template<> struct SmoothChecker<OP_REMAINDER>{ static const bool check=false;};
475  template<> struct SmoothChecker<OP_EQ>{ static const bool check=false;};
476  template<> struct SmoothChecker<OP_NE>{ static const bool check=false;};
477  template<> struct SmoothChecker<OP_SIGN>{ static const bool check=false;};
478  template<> struct SmoothChecker<OP_COPYSIGN>{ static const bool check=false;};
479  template<> struct SmoothChecker<OP_NOT>{ static const bool check=false;};
480  template<> struct SmoothChecker<OP_AND>{ static const bool check=false;};
481  template<> struct SmoothChecker<OP_OR>{ static const bool check=false;};
482  template<> struct SmoothChecker<OP_IF_ELSE_ZERO>{ static const bool check=false;};
484 
487  template<casadi_int I> struct F0XChecker { static const bool check=false;};
488  template<> struct F0XChecker<OP_ASSIGN>{ static const bool check=true;};
489  template<> struct F0XChecker<OP_MUL>{ static const bool check=true;};
490  template<> struct F0XChecker<OP_DIV>{ static const bool check=true;};
491  template<> struct F0XChecker<OP_NEG>{ static const bool check=true;};
492  template<> struct F0XChecker<OP_SQRT>{ static const bool check=true;};
493  template<> struct F0XChecker<OP_SQ>{ static const bool check=true;};
494  template<> struct F0XChecker<OP_TWICE>{ static const bool check=true;};
495  template<> struct F0XChecker<OP_SIN>{ static const bool check=true;};
496  template<> struct F0XChecker<OP_TAN>{ static const bool check=true;};
497  template<> struct F0XChecker<OP_ATAN>{ static const bool check=true;};
498  template<> struct F0XChecker<OP_ASIN>{ static const bool check=true;};
499  template<> struct F0XChecker<OP_FLOOR>{ static const bool check=true;};
500  template<> struct F0XChecker<OP_CEIL>{ static const bool check=true;};
501  template<> struct F0XChecker<OP_FMOD>{ static const bool check=true;};
502  template<> struct F0XChecker<OP_REMAINDER>{ static const bool check=true;};
503  template<> struct F0XChecker<OP_FABS>{ static const bool check=true;};
504  template<> struct F0XChecker<OP_SIGN>{ static const bool check=true;};
505  template<> struct F0XChecker<OP_COPYSIGN>{ static const bool check=true;};
506  template<> struct F0XChecker<OP_ERF>{ static const bool check=true;};
507  template<> struct F0XChecker<OP_SINH>{ static const bool check=true;};
508  template<> struct F0XChecker<OP_TANH>{ static const bool check=true;};
509  template<> struct F0XChecker<OP_ASINH>{ static const bool check=true;};
510  template<> struct F0XChecker<OP_ATANH>{ static const bool check=true;};
511  template<> struct F0XChecker<OP_ERFINV>{ static const bool check=true;};
512  template<> struct F0XChecker<OP_AND>{ static const bool check=true;};
513  template<> struct F0XChecker<OP_IF_ELSE_ZERO>{ static const bool check=true;};
514  template<> struct F0XChecker<OP_LOG1P>{ static const bool check=true;};
515  template<> struct F0XChecker<OP_EXPM1>{ static const bool check=true;};
517 
520  template<casadi_int I> struct FX0Checker { static const bool check=false;};
521  template<> struct FX0Checker<OP_MUL>{ static const bool check=true;};
522  template<> struct FX0Checker<OP_AND>{ static const bool check=true;};
523  template<> struct FX0Checker<OP_IF_ELSE_ZERO>{ static const bool check=true;};
525 
528  template<casadi_int I> struct F00Checker {
529  static const bool check=F0XChecker<I>::check || FX0Checker<I>::check;
530  };
531  template<> struct F00Checker<OP_ADD>{ static const bool check=true;};
532  template<> struct F00Checker<OP_SUB>{ static const bool check=true;};
533  template<> struct F00Checker<OP_FMIN>{ static const bool check=true;};
534  template<> struct F00Checker<OP_FMAX>{ static const bool check=true;};
535  template<> struct F00Checker<OP_AND>{ static const bool check=true;};
536  template<> struct F00Checker<OP_OR>{ static const bool check=true;};
537  template<> struct F00Checker<OP_COPYSIGN>{ static const bool check=true;};
538  template<> struct F00Checker<OP_LT>{ static const bool check=true;};
539  template<> struct F00Checker<OP_HYPOT>{ static const bool check=true;};
541 
544  template<casadi_int I> struct CommChecker { static const bool check=false;};
545  template<> struct CommChecker<OP_ADD>{ static const bool check=true;};
546  template<> struct CommChecker<OP_MUL>{ static const bool check=true;};
547  template<> struct CommChecker<OP_EQ>{ static const bool check=true;};
548  template<> struct CommChecker<OP_NE>{ static const bool check=true;};
549  template<> struct CommChecker<OP_AND>{ static const bool check=true;};
550  template<> struct CommChecker<OP_OR>{ static const bool check=true;};
551  template<> struct CommChecker<OP_HYPOT>{ static const bool check=true;};
552  template<> struct CommChecker<OP_FMIN>{ static const bool check=true;};
553  template<> struct CommChecker<OP_FMAX>{ static const bool check=true;};
555 
558  template<casadi_int I> struct NonnegativeChecker { static const bool check=false;};
559  template<> struct NonnegativeChecker<OP_SQRT>{ static const bool check=true;};
560  template<> struct NonnegativeChecker<OP_SQ>{ static const bool check=true;};
561  template<> struct NonnegativeChecker<OP_EXP>{ static const bool check=true;};
562  template<> struct NonnegativeChecker<OP_LT>{ static const bool check=true;};
563  template<> struct NonnegativeChecker<OP_LE>{ static const bool check=true;};
564  template<> struct NonnegativeChecker<OP_EQ>{ static const bool check=true;};
565  template<> struct NonnegativeChecker<OP_NE>{ static const bool check=true;};
566  template<> struct NonnegativeChecker<OP_NOT>{ static const bool check=true;};
567  template<> struct NonnegativeChecker<OP_AND>{ static const bool check=true;};
568  template<> struct NonnegativeChecker<OP_OR>{ static const bool check=true;};
569  template<> struct NonnegativeChecker<OP_HYPOT>{ static const bool check=true;};
570  template<> struct NonnegativeChecker<OP_FABS>{ static const bool check=true;};
572 
575  template<casadi_int I> struct BooleanChecker { static const bool check=false;};
576  template<> struct BooleanChecker<OP_LT>{ static const bool check=true;};
577  template<> struct BooleanChecker<OP_LE>{ static const bool check=true;};
578  template<> struct BooleanChecker<OP_EQ>{ static const bool check=true;};
579  template<> struct BooleanChecker<OP_NE>{ static const bool check=true;};
580  template<> struct BooleanChecker<OP_NOT>{ static const bool check=true;};
581  template<> struct BooleanChecker<OP_AND>{ static const bool check=true;};
582  template<> struct BooleanChecker<OP_OR>{ static const bool check=true;};
584 
587  template<casadi_int I> struct NargChecker { static const casadi_int check=1;};
588  template<> struct NargChecker<OP_ADD>{ static const casadi_int check=2;};
589  template<> struct NargChecker<OP_SUB>{ static const casadi_int check=2;};
590  template<> struct NargChecker<OP_MUL>{ static const casadi_int check=2;};
591  template<> struct NargChecker<OP_DIV>{ static const casadi_int check=2;};
592  template<> struct NargChecker<OP_POW>{ static const casadi_int check=2;};
593  template<> struct NargChecker<OP_CONSTPOW>{ static const casadi_int check=2;};
594  template<> struct NargChecker<OP_EQ>{ static const casadi_int check=2;};
595  template<> struct NargChecker<OP_NE>{ static const casadi_int check=2;};
596  template<> struct NargChecker<OP_AND>{ static const casadi_int check=2;};
597  template<> struct NargChecker<OP_OR>{ static const casadi_int check=2;};
598  template<> struct NargChecker<OP_FMIN>{ static const casadi_int check=2;};
599  template<> struct NargChecker<OP_FMAX>{ static const casadi_int check=2;};
600  template<> struct NargChecker<OP_PRINTME>{ static const casadi_int check=2;};
601  template<> struct NargChecker<OP_ATAN2>{ static const casadi_int check=2;};
602  template<> struct NargChecker<OP_IF_ELSE_ZERO>{ static const casadi_int check=2;};
603  template<> struct NargChecker<OP_FMOD>{ static const casadi_int check=2;};
604  template<> struct NargChecker<OP_REMAINDER>{ static const casadi_int check=2;};
605  template<> struct NargChecker<OP_COPYSIGN>{ static const casadi_int check=2;};
606  template<> struct NargChecker<OP_CONST>{ static const casadi_int check=0;};
607  template<> struct NargChecker<OP_PARAMETER>{ static const casadi_int check=0;};
608  template<> struct NargChecker<OP_INPUT>{ static const casadi_int check=0;};
609  template<> struct NargChecker<OP_HYPOT>{ static const casadi_int check=2;};
611 
613  template<>
614  struct UnaryOperation<OP_ASSIGN>{
615  public:
616  template<typename T> static inline void fcn(const T& x, T& f) { f = x;}
617  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 1; }
618  };
619 
621  template<>
622  struct BinaryOperation<OP_ADD>{
623  public:
624  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x+y;}
625  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
626  d[0]=d[1]=1;}
627  };
628 
630  template<>
631  struct BinaryOperation<OP_SUB>{
632  public:
633  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x-y;}
634  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
635  d[0]=1; d[1]=-1;}
636  };
637 
639  template<>
640  struct BinaryOperation<OP_MUL>{
641  public:
642  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x*y;}
643  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
644  d[0]=y; d[1]=x;}
645  };
646 
648  template<>
649  struct BinaryOperation<OP_DIV>{
650  public:
651  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x/y;}
652  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
653  d[0]=1/y; d[1]=-f/y;}
654  };
655 
657  template<>
658  struct UnaryOperation<OP_NEG>{
659  public:
660  template<typename T> static inline void fcn(const T& x, T& f) { f = -x;}
661  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=-1;}
662  };
663 
665  template<>
666  struct UnaryOperation<OP_EXP>{
667  public:
668  template<typename T> static inline void fcn(const T& x, T& f) { f = exp(x);}
669  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=f;}
670  };
671 
673  template<>
674  struct UnaryOperation<OP_LOG>{
675  public:
676  template<typename T> static inline void fcn(const T& x, T& f) { f = log(x);}
677  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=1/x;}
678  };
679 
681  template<>
682  struct BinaryOperation<OP_POW>{
683  public:
684  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = pow(x, y);}
685  // See issue #104 why d[0] is no longer y*f/x
686  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
687  d[0]=y*pow(x, y-1); d[1]=log(x)*f;}
688  };
689 
691  template<>
692  struct BinaryOperation<OP_CONSTPOW>{
693  public:
694  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = pow(x, y);}
695  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
696  d[0]=y*pow(x, y-1); d[1]=0;}
697  };
698 
700  template<>
701  struct UnaryOperation<OP_SQRT>{
702  public:
703  template<typename T> static inline void fcn(const T& x, T& f) { f = sqrt(x);}
704  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=1/(twice(f));}
705  };
706 
708  template<>
709  struct UnaryOperation<OP_SQ>{
710  public:
711  template<typename T> static inline void fcn(const T& x, T& f) { f = sq(x);}
712  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=twice(x);}
713  };
714 
716  template<>
717  struct UnaryOperation<OP_TWICE>{
718  template<typename T> static inline void fcn(const T& x, T& f) { f = 2.*x;}
719  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 2; }
720  };
721 
723  template<>
724  struct UnaryOperation<OP_SIN>{
725  public:
726  template<typename T> static inline void fcn(const T& x, T& f) { f = sin(x);}
727  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=cos(x);}
728  };
729 
731  template<>
732  struct UnaryOperation<OP_COS>{
733  public:
734  template<typename T> static inline void fcn(const T& x, T& f) { f = cos(x);}
735  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=-sin(x);}
736  };
737 
739  template<>
740  struct UnaryOperation<OP_TAN>{
741  public:
742  template<typename T> static inline void fcn(const T& x, T& f) { f = tan(x);}
743  template<typename T> static inline void der(const T& x, const T& f, T* d)
744  { d[0] = 1/sq(cos(x));}
745  };
746 
748  template<>
749  struct UnaryOperation<OP_ASIN>{
750  public:
751  template<typename T> static inline void fcn(const T& x, T& f) { f = asin(x);}
752  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=1/sqrt(1-x*x);}
753  };
754 
756  template<>
757  struct UnaryOperation<OP_ACOS>{
758  public:
759  template<typename T> static inline void fcn(const T& x, T& f) { f = acos(x);}
760  template<typename T> static inline void der(const T& x, const T& f, T* d)
761  { d[0]=-1/sqrt(1-x*x);}
762  };
763 
765  template<>
766  struct UnaryOperation<OP_ATAN>{
767  public:
768  template<typename T> static inline void fcn(const T& x, T& f) { f = atan(x);}
769  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 1/(1+x*x);}
770  };
771 
773  template<>
774  struct BinaryOperation<OP_LT>{
775  public:
776  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x < y;}
777  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
778  d[0]=d[1]=0;}
779  };
780 
782  template<>
783  struct BinaryOperation<OP_LE>{
784  public:
785  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x <= y;}
786  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
787  d[0]=d[1]=0;}
788  };
789 
791  template<>
792  struct UnaryOperation<OP_FLOOR>{
793  public:
794  template<typename T> static inline void fcn(const T& x, T& f) { f = floor(x);}
795  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 0;}
796  };
797 
799  template<>
800  struct UnaryOperation<OP_CEIL>{
801  public:
802  template<typename T> static inline void fcn(const T& x, T& f) { f = ceil(x);}
803  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 0;}
804  };
805 
807  template<>
808  struct BinaryOperation<OP_FMOD>{
809  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = fmod(x, y);}
810  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
811  d[0]=1; d[1]=(f-x)/y;}
812  };
813 
815  template<>
816  struct BinaryOperation<OP_REMAINDER>{
817  template<typename T> static inline void fcn(const T& x, const T& y, T& f) {
818  f = remainder(x, y);}
819  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
820  d[0]=1; d[1]=(f-x)/y;}
821  };
822 
824  template<>
825  struct BinaryOperation<OP_EQ>{
826  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x==y;}
827  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
828  d[0]=d[1]=0;}
829  };
830 
832  template<>
833  struct BinaryOperation<OP_NE>{
834  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x!=y;}
835  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
836  d[0]=d[1]=0;}
837  };
838 
840  template<>
841  struct UnaryOperation<OP_NOT>{
842  public:
843  template<typename T> static inline void fcn(const T& x, T& f) { f = !x;}
844  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 0;}
845  };
846 
848  template<>
849  struct BinaryOperation<OP_AND>{
850  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x && y;}
851  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
852  d[0]=d[1]=0;}
853  };
854 
856  template<>
857  struct BinaryOperation<OP_OR>{
858  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x || y;}
859  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
860  d[0]=d[1]=0;}
861  };
862 
864  template<>
865  struct UnaryOperation<OP_ERF>{
866  template<typename T> static inline void fcn(const T& x, T& f) { f = erf(x);}
867  template<typename T> static inline void der(const T& x, const T& f, T* d) {
868  d[0] = (2/sqrt(pi))*exp(-x*x);}
869  };
870 
872  template<>
873  struct UnaryOperation<OP_FABS>{
874  template<typename T> static inline void fcn(const T& x, T& f) { f = fabs(x);}
875  template<typename T> static inline void der(const T& x, const T& f, T* d) {
876  d[0]=sign(x);}
877  };
878 
880  template<>
881  struct UnaryOperation<OP_SIGN>{
882  template<typename T> static inline void fcn(const T& x, T& f) { f = sign(x);}
883  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=0;}
884  };
885 
887  template<>
888  struct BinaryOperation<OP_COPYSIGN>{
889  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = copysign(x, y);}
890  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
891  T e = 1; d[0]=copysign(e, y)*sign(x); d[1]=0;}
892  };
893 
895  template<>
896  struct BinaryOperation<OP_FMIN>{
897  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = fmin(x, y);}
898  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
899  T a = x<=y;
900  T b = y<=x;
901  T c = a+b;
902  d[0]=a/c; d[1]=b/c;}
903  };
904 
906  template<>
907  struct BinaryOperation<OP_FMAX>{
908  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = fmax(x, y);}
909  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
910  T a = y<=x;
911  T b = x<=y;
912  T c = a+b;
913  d[0]=a/c; d[1]=b/c;}
914  };
915 
917  template<>
918  struct UnaryOperation<OP_INV>{
919  template<typename T> static inline void fcn(const T& x, T& f) { f = 1./x;}
920  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = -f*f; }
921  };
922 
924  template<>
925  struct UnaryOperation<OP_SINH>{
926  template<typename T> static inline void fcn(const T& x, T& f) { f = sinh(x);}
927  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = cosh(x); }
928  };
929 
931  template<>
932  struct UnaryOperation<OP_COSH>{
933  template<typename T> static inline void fcn(const T& x, T& f) { f = cosh(x);}
934  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = sinh(x); }
935  };
936 
938  template<>
939  struct UnaryOperation<OP_TANH>{
940  template<typename T> static inline void fcn(const T& x, T& f) { f = tanh(x);}
941  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 1-f*f; }
942  };
943 
945  template<>
946  struct UnaryOperation<OP_ASINH>{
947  template<typename T> static inline void fcn(const T& x, T& f) { f = asinh(x);}
948  template<typename T> static inline void der(const T& x, const T& f, T* d) {
949  d[0] = 1/sqrt(1+x*x); }
950  };
951 
953  template<>
954  struct UnaryOperation<OP_ACOSH>{
955  template<typename T> static inline void fcn(const T& x, T& f) { f = acosh(x);}
956  template<typename T> static inline void der(const T& x, const T& f, T* d) {
957  d[0] = 1/sqrt(x-1)/sqrt(x+1); }
958  };
959 
961  template<>
962  struct UnaryOperation<OP_ATANH>{
963  template<typename T> static inline void fcn(const T& x, T& f) { f = atanh(x);}
964  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 1/(1-x*x); }
965  };
966 
968  template<>
969  struct UnaryOperation<OP_ERFINV>{
970  template<typename T> static inline void fcn(const T& x, T& f) { f = erfinv(x);}
971  template<typename T> static inline void der(const T& x, const T& f, T* d) {
972  d[0] = (sqrt(pi)/2)*exp(f*f); }
973  };
974 
976  template<>
977  struct UnaryOperation<OP_LOG1P>{
978  template<typename T> static inline void fcn(const T& x, T& f) { f = log1p(x);}
979  template<typename T> static inline void der(const T& x, const T& f, T* d) {
980  d[0] = 1/(1+x);}
981  };
982 
984  template<>
985  struct UnaryOperation<OP_EXPM1>{
986  template<typename T> static inline void fcn(const T& x, T& f) { f = expm1(x);}
987  template<typename T> static inline void der(const T& x, const T& f, T* d) {
988  d[0] = exp(x); }
989  };
990 
992  template<>
993  struct BinaryOperation<OP_PRINTME>{
994  template<typename T> static inline void fcn(const T& x, const T& y, T& f) {f = printme(x, y); }
995  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
996  d[0]=1; d[1]=0;}
997  };
998 
1000  template<>
1001  struct BinaryOperation<OP_ATAN2>{
1002  public:
1003  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = atan2(x, y);}
1004  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
1005  T t = x*x+y*y; d[0]=y/t; d[1]=-x/t;}
1006  };
1007 
1009  template<>
1010  struct BinaryOperation<OP_IF_ELSE_ZERO>{
1011  public:
1012  template<typename T> static inline void fcn(const T& x, const T& y, T& f) {
1013  f = if_else_zero(x, y);}
1014  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
1015  d[0]=0; d[1]=x;}
1016  };
1017 
1019  template<>
1020  struct BinaryOperation<OP_LIFT>{
1021  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x;}
1022  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
1023  d[0] = 1; d[1] = 0; }
1024  };
1025 
1027  template<>
1028  struct BinaryOperation<OP_HYPOT>{
1029  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = hypot(x, y);}
1030  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
1031  d[0] = x/f; d[1] = y/f; }
1032  };
1033 
1034  template<template<casadi_int> class F, typename T>
1035  T operation_getter(casadi_int op) {
1036  switch (static_cast<Operation>(op)) {
1037  case OP_ASSIGN: return F<OP_ASSIGN>::check;
1038  case OP_ADD: return F<OP_ADD>::check;
1039  case OP_SUB: return F<OP_SUB>::check;
1040  case OP_MUL: return F<OP_MUL>::check;
1041  case OP_DIV: return F<OP_DIV>::check;
1042  case OP_NEG: return F<OP_NEG>::check;
1043  case OP_EXP: return F<OP_EXP>::check;
1044  case OP_LOG: return F<OP_LOG>::check;
1045  case OP_POW: return F<OP_POW>::check;
1046  case OP_CONSTPOW: return F<OP_CONSTPOW>::check;
1047  case OP_SQRT: return F<OP_SQRT>::check;
1048  case OP_SQ: return F<OP_SQ>::check;
1049  case OP_TWICE: return F<OP_TWICE>::check;
1050  case OP_SIN: return F<OP_SIN>::check;
1051  case OP_COS: return F<OP_COS>::check;
1052  case OP_TAN: return F<OP_TAN>::check;
1053  case OP_ASIN: return F<OP_ASIN>::check;
1054  case OP_ACOS: return F<OP_ACOS>::check;
1055  case OP_ATAN: return F<OP_ATAN>::check;
1056  case OP_LT: return F<OP_LT>::check;
1057  case OP_LE: return F<OP_LE>::check;
1058  case OP_EQ: return F<OP_EQ>::check;
1059  case OP_NE: return F<OP_NE>::check;
1060  case OP_NOT: return F<OP_NOT>::check;
1061  case OP_AND: return F<OP_AND>::check;
1062  case OP_OR: return F<OP_OR>::check;
1063  case OP_FLOOR: return F<OP_FLOOR>::check;
1064  case OP_CEIL: return F<OP_CEIL>::check;
1065  case OP_FMOD: return F<OP_FMOD>::check;
1066  case OP_REMAINDER: return F<OP_REMAINDER>::check;
1067  case OP_FABS: return F<OP_FABS>::check;
1068  case OP_SIGN: return F<OP_SIGN>::check;
1069  case OP_COPYSIGN: return F<OP_COPYSIGN>::check;
1070  case OP_IF_ELSE_ZERO: return F<OP_IF_ELSE_ZERO>::check;
1071  case OP_ERF: return F<OP_ERF>::check;
1072  case OP_FMIN: return F<OP_FMIN>::check;
1073  case OP_FMAX: return F<OP_FMAX>::check;
1074  case OP_INV: return F<OP_INV>::check;
1075  case OP_SINH: return F<OP_SINH>::check;
1076  case OP_COSH: return F<OP_COSH>::check;
1077  case OP_TANH: return F<OP_TANH>::check;
1078  case OP_ASINH: return F<OP_ASINH>::check;
1079  case OP_ACOSH: return F<OP_ACOSH>::check;
1080  case OP_ATANH: return F<OP_ATANH>::check;
1081  case OP_ATAN2: return F<OP_ATAN2>::check;
1082  case OP_CONST: return F<OP_CONST>::check;
1083  case OP_INPUT: return F<OP_INPUT>::check;
1084  case OP_OUTPUT: return F<OP_OUTPUT>::check;
1085  case OP_PARAMETER: return F<OP_PARAMETER>::check;
1086  case OP_CALL: return F<OP_CALL>::check;
1087  case OP_FIND: return F<OP_FIND>::check;
1088  case OP_LOW: return F<OP_LOW>::check;
1089  case OP_MAP: return F<OP_MAP>::check;
1090  case OP_MTIMES: return F<OP_MTIMES>::check;
1091  case OP_SOLVE: return F<OP_SOLVE>::check;
1092  case OP_TRANSPOSE: return F<OP_TRANSPOSE>::check;
1093  case OP_DETERMINANT: return F<OP_DETERMINANT>::check;
1094  case OP_INVERSE: return F<OP_INVERSE>::check;
1095  case OP_DOT: return F<OP_DOT>::check;
1096  case OP_BILIN: return F<OP_BILIN>::check;
1097  case OP_RANK1: return F<OP_RANK1>::check;
1098  case OP_HORZCAT: return F<OP_HORZCAT>::check;
1099  case OP_VERTCAT: return F<OP_VERTCAT>::check;
1100  case OP_DIAGCAT: return F<OP_DIAGCAT>::check;
1101  case OP_HORZSPLIT: return F<OP_HORZSPLIT>::check;
1102  case OP_VERTSPLIT: return F<OP_VERTSPLIT>::check;
1103  case OP_DIAGSPLIT: return F<OP_DIAGSPLIT>::check;
1104  case OP_RESHAPE: return F<OP_RESHAPE>::check;
1105  case OP_SPARSITY_CAST: return F<OP_SPARSITY_CAST>::check;
1106  case OP_SUBREF: return F<OP_SUBREF>::check;
1107  case OP_SUBASSIGN: return F<OP_SUBASSIGN>::check;
1108  case OP_GETNONZEROS: return F<OP_GETNONZEROS>::check;
1109  case OP_GETNONZEROS_PARAM: return F<OP_GETNONZEROS_PARAM>::check;
1110  case OP_ADDNONZEROS: return F<OP_ADDNONZEROS>::check;
1111  case OP_ADDNONZEROS_PARAM: return F<OP_ADDNONZEROS>::check;
1112  case OP_SETNONZEROS: return F<OP_SETNONZEROS>::check;
1113  case OP_SETNONZEROS_PARAM: return F<OP_SETNONZEROS>::check;
1114  case OP_PROJECT: return F<OP_PROJECT>::check;
1115  case OP_ASSERTION: return F<OP_ASSERTION>::check;
1116  case OP_MONITOR: return F<OP_MONITOR>::check;
1117  case OP_DUMP: return F<OP_DUMP>::check;
1118  case OP_NORM2: return F<OP_NORM2>::check;
1119  case OP_NORM1: return F<OP_NORM1>::check;
1120  case OP_NORMINF: return F<OP_NORMINF>::check;
1121  case OP_NORMF: return F<OP_NORMF>::check;
1122  case OP_MMIN: return F<OP_MMIN>::check;
1123  case OP_MMAX: return F<OP_MMAX>::check;
1124  case OP_HORZREPMAT: return F<OP_HORZREPMAT>::check;
1125  case OP_HORZREPSUM: return F<OP_HORZREPSUM>::check;
1126  case OP_ERFINV: return F<OP_ERFINV>::check;
1127  case OP_PRINTME: return F<OP_PRINTME>::check;
1128  case OP_LIFT: return F<OP_LIFT>::check;
1129  case OP_EINSTEIN: return F<OP_EINSTEIN>::check;
1130  case OP_BSPLINE: return F<OP_BSPLINE>::check;
1131  case OP_CONVEXIFY: return F<OP_CONVEXIFY>::check;
1132  case OP_LOG1P: return F<OP_LOG1P>::check;
1133  case OP_EXPM1: return F<OP_EXPM1>::check;
1134  case OP_HYPOT: return F<OP_HYPOT>::check;
1135  case OP_LOGSUMEXP: return F<OP_LOGSUMEXP>::check;
1136  case OP_KRON: return F<OP_KRON>::check;
1137  case OP_KRON_CONTRACT: return F<OP_KRON_CONTRACT>::check;
1138  }
1139  return T();
1140  }
1141 
1142  template<template<casadi_int> class F>
1143  bool operation_checker(casadi_int op) {
1144  return operation_getter<F, bool>(op);
1145  }
1146 
1148  template<typename T>
1149  struct casadi_math {
1150 
1154  static inline void fun(unsigned char op, const T& x, const T& y, T& f);
1155 
1159  static inline void fun(unsigned char op, const T* x, const T* y, T* f, casadi_int n);
1160 
1164  static inline void fun(unsigned char op, const T* x, const T& y, T* f, casadi_int n);
1165 
1169  static inline void fun(unsigned char op, const T& x, const T* y, T* f, casadi_int n);
1170 
1174  static inline void der(unsigned char op, const T& x, const T& y, const T& f, T* d);
1175 
1179  static inline void derF(unsigned char op, const T& x, const T& y, T& f, T* d);
1180 
1184  static inline void fun_linear(unsigned char op, const T*x, const T* y, T* f);
1185 
1189  static inline bool is_binary(unsigned char op);
1190 
1194  static inline bool is_unary(unsigned char op);
1195 
1199  static inline casadi_int ndeps(unsigned char op);
1200 
1204  static inline std::string print(unsigned char op, const std::string& x,
1205  const std::string& y);
1206  static inline std::string print(unsigned char op, const std::string& x);
1207  static inline std::string name(unsigned char op);
1208  static inline std::string pre(unsigned char op);
1209  static inline std::string sep(unsigned char op);
1210  static inline std::string post(unsigned char op);
1211  };
1212 
1214  template<>
1215  struct casadi_math<casadi_int>{
1216 
1220  static inline void fun(unsigned char op, const casadi_int& x,
1221  const casadi_int& y, casadi_int& f) {
1222  double ff(0);
1223  casadi_math<double>::fun(op, static_cast<double>(x), static_cast<double>(y), ff);
1224  f = static_cast<casadi_int>(ff);
1225  }
1226 
1227  static inline void fun(unsigned char op, const casadi_int* x, const casadi_int* y,
1228  casadi_int* f, casadi_int n) {
1229  for (casadi_int i=0; i<n; ++i) {
1230  double ff(0);
1231  casadi_math<double>::fun(op, static_cast<double>(*x++), static_cast<double>(*y++), ff);
1232  *f++ = static_cast<casadi_int>(ff);
1233  }
1234  }
1235 
1236  static inline void fun(unsigned char op, const casadi_int* x, const casadi_int& y,
1237  casadi_int* f, casadi_int n) {
1238  for (casadi_int i=0; i<n; ++i) {
1239  double ff;
1240  casadi_math<double>::fun(op, static_cast<double>(*x++), static_cast<double>(y), ff);
1241  *f++ = static_cast<casadi_int>(ff);
1242  }
1243  }
1244 
1245  static inline void fun(unsigned char op, const casadi_int& x, const casadi_int* y,
1246  casadi_int* f, casadi_int n) {
1247  for (casadi_int i=0; i<n; ++i) {
1248  double ff;
1249  casadi_math<double>::fun(op, static_cast<double>(x), static_cast<double>(*y++), ff);
1250  *f++ = static_cast<casadi_int>(ff);
1251  }
1252  }
1253 
1257  static inline void der(unsigned char op, const casadi_int& x, const casadi_int& y,
1258  const casadi_int& f, casadi_int* d) {
1259  double d_real[2] = {static_cast<double>(d[0]), static_cast<double>(d[1])};
1260  casadi_math<double>::der(op, static_cast<double>(x), static_cast<double>(y),
1261  static_cast<double>(f), d_real);
1262  d[0] = static_cast<casadi_int>(d_real[0]);
1263  d[1] = static_cast<casadi_int>(d_real[1]);
1264  }
1265 
1269  static inline void derF(unsigned char op, const casadi_int& x, const casadi_int& y,
1270  casadi_int& f, casadi_int* d) {
1271  double d_real[2] = {static_cast<double>(d[0]), static_cast<double>(d[1])};
1272  double f_real = static_cast<double>(f);
1273  casadi_math<double>::derF(op, static_cast<double>(x), static_cast<double>(y), f_real, d_real);
1274  f = static_cast<casadi_int>(f_real);
1275  d[0] = static_cast<casadi_int>(d_real[0]);
1276  d[1] = static_cast<casadi_int>(d_real[1]);
1277  }
1278 
1282  static inline casadi_int ndeps(unsigned char op) {
1283  return casadi_math<double>::ndeps(op);
1284  }
1285 
1289  static inline std::string print(unsigned char op, const std::string& x,
1290  const std::string& y) {
1291  return casadi_math<double>::print(op, x, y);
1292  }
1293  static inline std::string print(unsigned char op, const std::string& x) {
1294  return casadi_math<double>::print(op, x);
1295  }
1296  static inline std::string pre(unsigned char op) {
1297  return casadi_math<double>::pre(op);
1298  }
1299  static inline std::string name(unsigned char op) {
1300  return casadi_math<double>::name(op);
1301  }
1302  static inline std::string sep(unsigned char op) {
1303  return casadi_math<double>::sep(op);
1304  }
1305  static inline std::string post(unsigned char op) {
1306  return casadi_math<double>::post(op);
1307  }
1308  };
1309 
1310  // Template implementations
1311 
1312  template<typename T>
1313  inline void casadi_math<T>::fun(unsigned char op, const T& x, const T& y, T& f) {
1314  // NOTE: We define the implementation in a preprocessor macro to be able to force inlining,
1315  // and to allow extensions in the VM
1316 #define CASADI_MATH_FUN_BUILTIN_GEN(CNAME, X, Y, F, N) \
1317  case OP_ASSIGN: CNAME<OP_ASSIGN>::fcn(X, Y, F, N); break; \
1318  case OP_ADD: CNAME<OP_ADD>::fcn(X, Y, F, N); break; \
1319  case OP_SUB: CNAME<OP_SUB>::fcn(X, Y, F, N); break; \
1320  case OP_MUL: CNAME<OP_MUL>::fcn(X, Y, F, N); break; \
1321  case OP_DIV: CNAME<OP_DIV>::fcn(X, Y, F, N); break; \
1322  case OP_NEG: CNAME<OP_NEG>::fcn(X, Y, F, N); break; \
1323  case OP_EXP: CNAME<OP_EXP>::fcn(X, Y, F, N); break; \
1324  case OP_LOG: CNAME<OP_LOG>::fcn(X, Y, F, N); break; \
1325  case OP_POW: CNAME<OP_POW>::fcn(X, Y, F, N); break; \
1326  case OP_CONSTPOW: CNAME<OP_CONSTPOW>::fcn(X, Y, F, N); break; \
1327  case OP_SQRT: CNAME<OP_SQRT>::fcn(X, Y, F, N); break; \
1328  case OP_SQ: CNAME<OP_SQ>::fcn(X, Y, F, N); break; \
1329  case OP_TWICE: CNAME<OP_TWICE>::fcn(X, Y, F, N); break; \
1330  case OP_SIN: CNAME<OP_SIN>::fcn(X, Y, F, N); break; \
1331  case OP_COS: CNAME<OP_COS>::fcn(X, Y, F, N); break; \
1332  case OP_TAN: CNAME<OP_TAN>::fcn(X, Y, F, N); break; \
1333  case OP_ASIN: CNAME<OP_ASIN>::fcn(X, Y, F, N); break; \
1334  case OP_ACOS: CNAME<OP_ACOS>::fcn(X, Y, F, N); break; \
1335  case OP_ATAN: CNAME<OP_ATAN>::fcn(X, Y, F, N); break; \
1336  case OP_LT: CNAME<OP_LT>::fcn(X, Y, F, N); break; \
1337  case OP_LE: CNAME<OP_LE>::fcn(X, Y, F, N); break; \
1338  case OP_EQ: CNAME<OP_EQ>::fcn(X, Y, F, N); break; \
1339  case OP_NE: CNAME<OP_NE>::fcn(X, Y, F, N); break; \
1340  case OP_NOT: CNAME<OP_NOT>::fcn(X, Y, F, N); break; \
1341  case OP_AND: CNAME<OP_AND>::fcn(X, Y, F, N); break; \
1342  case OP_OR: CNAME<OP_OR>::fcn(X, Y, F, N); break; \
1343  case OP_IF_ELSE_ZERO: CNAME<OP_IF_ELSE_ZERO>::fcn(X, Y, F, N); break; \
1344  case OP_FLOOR: CNAME<OP_FLOOR>::fcn(X, Y, F, N); break; \
1345  case OP_CEIL: CNAME<OP_CEIL>::fcn(X, Y, F, N); break; \
1346  case OP_FMOD: CNAME<OP_FMOD>::fcn(X, Y, F, N); break; \
1347  case OP_REMAINDER: CNAME<OP_REMAINDER>::fcn(X, Y, F, N); break; \
1348  case OP_FABS: CNAME<OP_FABS>::fcn(X, Y, F, N); break; \
1349  case OP_SIGN: CNAME<OP_SIGN>::fcn(X, Y, F, N); break; \
1350  case OP_COPYSIGN: CNAME<OP_COPYSIGN>::fcn(X, Y, F, N); break; \
1351  case OP_ERF: CNAME<OP_ERF>::fcn(X, Y, F, N); break; \
1352  case OP_FMIN: CNAME<OP_FMIN>::fcn(X, Y, F, N); break; \
1353  case OP_FMAX: CNAME<OP_FMAX>::fcn(X, Y, F, N); break; \
1354  case OP_INV: CNAME<OP_INV>::fcn(X, Y, F, N); break; \
1355  case OP_SINH: CNAME<OP_SINH>::fcn(X, Y, F, N); break; \
1356  case OP_COSH: CNAME<OP_COSH>::fcn(X, Y, F, N); break; \
1357  case OP_TANH: CNAME<OP_TANH>::fcn(X, Y, F, N); break; \
1358  case OP_ASINH: CNAME<OP_ASINH>::fcn(X, Y, F, N); break; \
1359  case OP_ACOSH: CNAME<OP_ACOSH>::fcn(X, Y, F, N); break; \
1360  case OP_ATANH: CNAME<OP_ATANH>::fcn(X, Y, F, N); break; \
1361  case OP_ATAN2: CNAME<OP_ATAN2>::fcn(X, Y, F, N); break; \
1362  case OP_ERFINV: CNAME<OP_ERFINV>::fcn(X, Y, F, N); break; \
1363  case OP_LIFT: CNAME<OP_LIFT>::fcn(X, Y, F, N); break; \
1364  case OP_PRINTME: CNAME<OP_PRINTME>::fcn(X, Y, F, N); break; \
1365  case OP_LOG1P: CNAME<OP_LOG1P>::fcn(X, Y, F, N); break; \
1366  case OP_EXPM1: CNAME<OP_EXPM1>::fcn(X, Y, F, N); break; \
1367  case OP_HYPOT: CNAME<OP_HYPOT>::fcn(X, Y, F, N); break;
1368 
1369 #define CASADI_MATH_FUN_BUILTIN(X, Y, F) CASADI_MATH_FUN_BUILTIN_GEN(BinaryOperationSS, X, Y, F, 1)
1370 
1371  switch (op) {
1372  CASADI_MATH_FUN_BUILTIN(x, y, f)
1373  }
1374  }
1375 
1376  template<typename T>
1377  inline void casadi_math<T>::fun(unsigned char op, const T* x, const T* y, T* f, casadi_int n) {
1378  switch (op) {
1379  CASADI_MATH_FUN_BUILTIN_GEN(BinaryOperationVV, x, y, f, n)
1380  }
1381  }
1382 
1383  template<typename T>
1384  inline void casadi_math<T>::fun(unsigned char op, const T* x, const T& y, T* f, casadi_int n) {
1385  switch (op) {
1386  CASADI_MATH_FUN_BUILTIN_GEN(BinaryOperationVS, x, y, f, n)
1387  }
1388  }
1389 
1390  template<typename T>
1391  inline void casadi_math<T>::fun(unsigned char op, const T& x, const T* y, T* f, casadi_int n) {
1392  switch (op) {
1393  CASADI_MATH_FUN_BUILTIN_GEN(BinaryOperationSV, x, y, f, n)
1394  }
1395  }
1396 
1397 
1398  template<typename T>
1399  inline void casadi_math<T>::der(unsigned char op, const T& x, const T& y, const T& f, T* d) {
1400  // NOTE: We define the implementation in a preprocessor macro to be able to force inlining,
1401  // and to allow extensions in the VM
1402 #define CASADI_MATH_DER_BUILTIN(X, Y, F, D) \
1403  case OP_ASSIGN: BinaryOperation<OP_ASSIGN>::der(X, Y, F, D); break; \
1404  case OP_ADD: BinaryOperation<OP_ADD>::der(X, Y, F, D); break; \
1405  case OP_SUB: BinaryOperation<OP_SUB>::der(X, Y, F, D); break; \
1406  case OP_MUL: BinaryOperation<OP_MUL>::der(X, Y, F, D); break; \
1407  case OP_DIV: BinaryOperation<OP_DIV>::der(X, Y, F, D); break; \
1408  case OP_NEG: BinaryOperation<OP_NEG>::der(X, Y, F, D); break; \
1409  case OP_EXP: BinaryOperation<OP_EXP>::der(X, Y, F, D); break; \
1410  case OP_LOG: BinaryOperation<OP_LOG>::der(X, Y, F, D); break; \
1411  case OP_POW: BinaryOperation<OP_POW>::der(X, Y, F, D); break; \
1412  case OP_CONSTPOW: BinaryOperation<OP_CONSTPOW>::der(X, Y, F, D); break; \
1413  case OP_SQRT: BinaryOperation<OP_SQRT>::der(X, Y, F, D); break; \
1414  case OP_SQ: BinaryOperation<OP_SQ>::der(X, Y, F, D); break; \
1415  case OP_TWICE: BinaryOperation<OP_TWICE>::der(X, Y, F, D); break; \
1416  case OP_SIN: BinaryOperation<OP_SIN>::der(X, Y, F, D); break; \
1417  case OP_COS: BinaryOperation<OP_COS>::der(X, Y, F, D); break; \
1418  case OP_TAN: BinaryOperation<OP_TAN>::der(X, Y, F, D); break; \
1419  case OP_ASIN: BinaryOperation<OP_ASIN>::der(X, Y, F, D); break; \
1420  case OP_ACOS: BinaryOperation<OP_ACOS>::der(X, Y, F, D); break; \
1421  case OP_ATAN: BinaryOperation<OP_ATAN>::der(X, Y, F, D); break; \
1422  case OP_LT: BinaryOperation<OP_LT>::der(X, Y, F, D); break; \
1423  case OP_LE: BinaryOperation<OP_LE>::der(X, Y, F, D); break; \
1424  case OP_EQ: BinaryOperation<OP_EQ>::der(X, Y, F, D); break; \
1425  case OP_NE: BinaryOperation<OP_NE>::der(X, Y, F, D); break; \
1426  case OP_NOT: BinaryOperation<OP_NOT>::der(X, Y, F, D); break; \
1427  case OP_AND: BinaryOperation<OP_AND>::der(X, Y, F, D); break; \
1428  case OP_OR: BinaryOperation<OP_OR>::der(X, Y, F, D); break; \
1429  case OP_IF_ELSE_ZERO: BinaryOperation<OP_IF_ELSE_ZERO>::der(X, Y, F, D); break; \
1430  case OP_FLOOR: BinaryOperation<OP_FLOOR>::der(X, Y, F, D); break; \
1431  case OP_CEIL: BinaryOperation<OP_CEIL>::der(X, Y, F, D); break; \
1432  case OP_FMOD: BinaryOperation<OP_FMOD>::der(X, Y, F, D); break; \
1433  case OP_REMAINDER: BinaryOperation<OP_REMAINDER>::der(X, Y, F, D); break; \
1434  case OP_FABS: BinaryOperation<OP_FABS>::der(X, Y, F, D); break; \
1435  case OP_SIGN: BinaryOperation<OP_SIGN>::der(X, Y, F, D); break; \
1436  case OP_COPYSIGN: BinaryOperation<OP_COPYSIGN>::der(X, Y, F, D); break; \
1437  case OP_ERF: BinaryOperation<OP_ERF>::der(X, Y, F, D); break; \
1438  case OP_FMIN: BinaryOperation<OP_FMIN>::der(X, Y, F, D); break; \
1439  case OP_FMAX: BinaryOperation<OP_FMAX>::der(X, Y, F, D); break; \
1440  case OP_INV: BinaryOperation<OP_INV>::der(X, Y, F, D); break; \
1441  case OP_SINH: BinaryOperation<OP_SINH>::der(X, Y, F, D); break; \
1442  case OP_COSH: BinaryOperation<OP_COSH>::der(X, Y, F, D); break; \
1443  case OP_TANH: BinaryOperation<OP_TANH>::der(X, Y, F, D); break; \
1444  case OP_ASINH: BinaryOperation<OP_ASINH>::der(X, Y, F, D); break; \
1445  case OP_ACOSH: BinaryOperation<OP_ACOSH>::der(X, Y, F, D); break; \
1446  case OP_ATANH: BinaryOperation<OP_ATANH>::der(X, Y, F, D); break; \
1447  case OP_ATAN2: BinaryOperation<OP_ATAN2>::der(X, Y, F, D); break; \
1448  case OP_ERFINV: BinaryOperation<OP_ERFINV>::der(X, Y, F, D); break; \
1449  case OP_LIFT: BinaryOperation<OP_LIFT>::der(X, Y, F, D); break; \
1450  case OP_PRINTME: BinaryOperation<OP_PRINTME>::der(X, Y, F, D); break; \
1451  case OP_LOG1P: BinaryOperation<OP_LOG1P>::der(X, Y, F, D); break; \
1452  case OP_EXPM1: BinaryOperation<OP_EXPM1>::der(X, Y, F, D); break; \
1453  case OP_HYPOT: BinaryOperation<OP_HYPOT>::der(X, Y, F, D); break;
1454  switch (op) {
1455  CASADI_MATH_DER_BUILTIN(x, y, f, d)
1456  }
1457  }
1458 
1459 
1460  template<typename T>
1461  inline void casadi_math<T>::derF(unsigned char op, const T& x, const T& y, T& f, T* d) {
1462  // NOTE: We define the implementation in a preprocessor macro to be able to force inlining,
1463  // and to allow extensions in the VM
1464 #define CASADI_MATH_DERF_BUILTIN(X, Y, F, D) \
1465 case OP_ASSIGN: DerBinaryOperation<OP_ASSIGN>::derf(X, Y, F, D); break; \
1466 case OP_ADD: DerBinaryOperation<OP_ADD>::derf(X, Y, F, D); break; \
1467 case OP_SUB: DerBinaryOperation<OP_SUB>::derf(X, Y, F, D); break; \
1468 case OP_MUL: DerBinaryOperation<OP_MUL>::derf(X, Y, F, D); break; \
1469 case OP_DIV: DerBinaryOperation<OP_DIV>::derf(X, Y, F, D); break; \
1470 case OP_NEG: DerBinaryOperation<OP_NEG>::derf(X, Y, F, D); break; \
1471 case OP_EXP: DerBinaryOperation<OP_EXP>::derf(X, Y, F, D); break; \
1472 case OP_LOG: DerBinaryOperation<OP_LOG>::derf(X, Y, F, D); break; \
1473 case OP_POW: DerBinaryOperation<OP_POW>::derf(X, Y, F, D); break; \
1474 case OP_CONSTPOW: DerBinaryOperation<OP_CONSTPOW>::derf(X, Y, F, D); break; \
1475 case OP_SQRT: DerBinaryOperation<OP_SQRT>::derf(X, Y, F, D); break; \
1476 case OP_SQ: DerBinaryOperation<OP_SQ>::derf(X, Y, F, D); break; \
1477 case OP_TWICE: DerBinaryOperation<OP_TWICE>::derf(X, Y, F, D); break; \
1478 case OP_SIN: DerBinaryOperation<OP_SIN>::derf(X, Y, F, D); break; \
1479 case OP_COS: DerBinaryOperation<OP_COS>::derf(X, Y, F, D); break; \
1480 case OP_TAN: DerBinaryOperation<OP_TAN>::derf(X, Y, F, D); break; \
1481 case OP_ASIN: DerBinaryOperation<OP_ASIN>::derf(X, Y, F, D); break; \
1482 case OP_ACOS: DerBinaryOperation<OP_ACOS>::derf(X, Y, F, D); break; \
1483 case OP_ATAN: DerBinaryOperation<OP_ATAN>::derf(X, Y, F, D); break; \
1484 case OP_LT: DerBinaryOperation<OP_LT>::derf(X, Y, F, D); break; \
1485 case OP_LE: DerBinaryOperation<OP_LE>::derf(X, Y, F, D); break; \
1486 case OP_EQ: DerBinaryOperation<OP_EQ>::derf(X, Y, F, D); break; \
1487 case OP_NE: DerBinaryOperation<OP_NE>::derf(X, Y, F, D); break; \
1488 case OP_NOT: DerBinaryOperation<OP_NOT>::derf(X, Y, F, D); break; \
1489 case OP_AND: DerBinaryOperation<OP_AND>::derf(X, Y, F, D); break; \
1490 case OP_OR: DerBinaryOperation<OP_OR>::derf(X, Y, F, D); break; \
1491 case OP_IF_ELSE_ZERO: DerBinaryOperation<OP_IF_ELSE_ZERO>::derf(X, Y, F, D); break; \
1492 case OP_FLOOR: DerBinaryOperation<OP_FLOOR>::derf(X, Y, F, D); break; \
1493 case OP_CEIL: DerBinaryOperation<OP_CEIL>::derf(X, Y, F, D); break; \
1494 case OP_FMOD: DerBinaryOperation<OP_FMOD>::derf(X, Y, F, D); break; \
1495 case OP_REMAINDER: DerBinaryOperation<OP_REMAINDER>::derf(X, Y, F, D); break; \
1496 case OP_FABS: DerBinaryOperation<OP_FABS>::derf(X, Y, F, D); break; \
1497 case OP_SIGN: DerBinaryOperation<OP_SIGN>::derf(X, Y, F, D); break; \
1498 case OP_COPYSIGN: DerBinaryOperation<OP_COPYSIGN>::derf(X, Y, F, D); break; \
1499 case OP_ERF: DerBinaryOperation<OP_ERF>::derf(X, Y, F, D); break; \
1500 case OP_FMIN: DerBinaryOperation<OP_FMIN>::derf(X, Y, F, D); break; \
1501 case OP_FMAX: DerBinaryOperation<OP_FMAX>::derf(X, Y, F, D); break; \
1502 case OP_INV: DerBinaryOperation<OP_INV>::derf(X, Y, F, D); break; \
1503 case OP_SINH: DerBinaryOperation<OP_SINH>::derf(X, Y, F, D); break; \
1504 case OP_COSH: DerBinaryOperation<OP_COSH>::derf(X, Y, F, D); break; \
1505 case OP_TANH: DerBinaryOperation<OP_TANH>::derf(X, Y, F, D); break; \
1506 case OP_ASINH: DerBinaryOperation<OP_ASINH>::derf(X, Y, F, D); break; \
1507 case OP_ACOSH: DerBinaryOperation<OP_ACOSH>::derf(X, Y, F, D); break; \
1508 case OP_ATANH: DerBinaryOperation<OP_ATANH>::derf(X, Y, F, D); break; \
1509 case OP_ATAN2: DerBinaryOperation<OP_ATAN2>::derf(X, Y, F, D); break; \
1510 case OP_ERFINV: DerBinaryOperation<OP_ERFINV>::derf(X, Y, F, D); break; \
1511 case OP_LIFT: DerBinaryOperation<OP_LIFT>::derf(X, Y, F, D); break; \
1512 case OP_PRINTME: DerBinaryOperation<OP_PRINTME>::derf(X, Y, F, D); break; \
1513 case OP_LOG1P: DerBinaryOperation<OP_LOG1P>::derf(X, Y, F, D); break; \
1514 case OP_EXPM1: DerBinaryOperation<OP_EXPM1>::derf(X, Y, F, D); break; \
1515 case OP_HYPOT: DerBinaryOperation<OP_HYPOT>::derf(X, Y, F, D); break;
1516  switch (op) {
1517  CASADI_MATH_DERF_BUILTIN(x, y, f, d)
1518  }
1519  }
1520 
1521  #define CASADI_MATH_BINARY_BUILTIN \
1522  case OP_ADD: \
1523  case OP_SUB: \
1524  case OP_MUL: \
1525  case OP_DIV: \
1526  case OP_POW: \
1527  case OP_CONSTPOW: \
1528  case OP_LT: \
1529  case OP_LE: \
1530  case OP_EQ: \
1531  case OP_NE: \
1532  case OP_AND: \
1533  case OP_OR: \
1534  case OP_COPYSIGN: \
1535  case OP_FMOD: \
1536  case OP_REMAINDER: \
1537  case OP_FMIN: \
1538  case OP_FMAX: \
1539  case OP_ATAN2: \
1540  case OP_PRINTME: \
1541  case OP_LIFT: \
1542  case OP_HYPOT:
1543 
1544  #define CASADI_MATH_UNARY_BUILTIN \
1545  case OP_ASSIGN: \
1546  case OP_NEG: \
1547  case OP_EXP: \
1548  case OP_LOG: \
1549  case OP_SQRT: \
1550  case OP_SQ: \
1551  case OP_TWICE: \
1552  case OP_SIN: \
1553  case OP_COS: \
1554  case OP_TAN: \
1555  case OP_ASIN: \
1556  case OP_ACOS: \
1557  case OP_ATAN: \
1558  case OP_FLOOR: \
1559  case OP_CEIL: \
1560  case OP_NOT: \
1561  case OP_ERF: \
1562  case OP_FABS: \
1563  case OP_SIGN: \
1564  case OP_INV: \
1565  case OP_SINH: \
1566  case OP_COSH: \
1567  case OP_TANH: \
1568  case OP_ASINH: \
1569  case OP_ACOSH: \
1570  case OP_ATANH: \
1571  case OP_ERFINV: \
1572  case OP_LOG1P: \
1573  case OP_EXPM1:
1574 
1575  template<typename T>
1576  inline void casadi_math<T>::fun_linear(unsigned char op, const T* x, const T* y, T* f) {
1577  if (op==OP_ADD || op==OP_SUB) {
1578  for (int i=0;i<3;++i) {
1579  f[i] = T::binary(op, x[i], y[i]);
1580  }
1581  } else if (op==OP_TWICE || op==OP_NEG) {
1582  for (int i=0;i<3;++i) {
1583  f[i] = T::unary(op, x[i]);
1584  }
1585  } else if (op==OP_MUL) {
1586  f[0] += x[0]*y[0];
1587  f[1] += x[0]*y[1];
1588  f[2] += x[0]*y[2];
1589  f[1] += x[1]*y[0];
1590  f[2] += x[1]*y[1];
1591  f[2] += x[1]*y[2];
1592  f[2] += x[2]*y[0];
1593  f[2] += x[2]*y[1];
1594  f[2] += x[2]*y[2];
1595  } else if (op==OP_DIV) {
1596  bool const_argy = y[1].is_zero() && y[2].is_zero();
1597  if (const_argy) {
1598  f[0] = x[0]/y[0];
1599  f[1] = x[1]/y[0];
1600  f[2] = x[2]/y[0];
1601  } else {
1602  f[2] = (x[0]+x[1]+x[2])/(y[0]+y[1]+y[2]);
1603  }
1604  } else if (casadi_math<T>::is_unary(op)) {
1605  bool const_arg = x[1].is_zero() && x[2].is_zero();
1606  if (const_arg) {
1607  f[0] = T::unary(op, x[0]);
1608  } else {
1609  f[2] = T::unary(op, x[0]+x[1]+x[2]);
1610  }
1611 
1612  } else if (casadi_math<T>::is_binary(op)) {
1613  bool const_argx = x[1].is_zero() && x[2].is_zero();
1614  bool const_argy = y[1].is_zero() && y[2].is_zero();
1615  if (const_argx && const_argy) {
1616  f[0] = T::binary(op, x[0], y[0]);
1617  } else {
1618  f[2] = T::binary(op, x[0]+x[1]+x[2], y[0]+y[1]+y[2]);
1619  }
1620  } else {
1621  casadi_error("Not implemented");
1622  }
1623  }
1624 
1625  template<typename T>
1626  bool casadi_math<T>::is_binary(unsigned char op) {
1627  switch (op) {
1628  CASADI_MATH_BINARY_BUILTIN
1629  case OP_IF_ELSE_ZERO:
1630  return true;
1631  default:
1632  return false;
1633  }
1634  }
1635 
1636  template<typename T>
1637  bool casadi_math<T>::is_unary(unsigned char op) {
1638  switch (op) {
1639  CASADI_MATH_UNARY_BUILTIN
1640  return true;
1641  default:
1642  return false;
1643  }
1644  }
1645 
1646  template<typename T>
1647  inline casadi_int casadi_math<T>::ndeps(unsigned char op) {
1648  switch (op) {
1649  case OP_CONST:
1650  case OP_PARAMETER:
1651  case OP_INPUT:
1652  return 0;
1653  CASADI_MATH_BINARY_BUILTIN
1654  case OP_IF_ELSE_ZERO:
1655  return 2;
1656  case OP_CALL:
1657  return -1;
1658  default:
1659  return 1;
1660  }
1661  }
1662 
1663  template<typename T>
1664  inline std::string
1665  casadi_math<T>::print(unsigned char op,
1666  const std::string& x, const std::string& y) {
1667  casadi_assert_dev(ndeps(op)==2);
1668  return pre(op) + x + sep(op) + y + post(op);
1669  }
1670 
1671  template<typename T>
1672  inline std::string
1673  casadi_math<T>::print(unsigned char op, const std::string& x) {
1674  casadi_assert_dev(ndeps(op)==1);
1675  return pre(op) + x + post(op);
1676  }
1677 
1678  template<typename T>
1679  inline std::string casadi_math<T>::name(unsigned char op) {
1680  switch (op) {
1681  case OP_ASSIGN: return "assign";
1682  case OP_ADD: return "add";
1683  case OP_SUB: return "sub";
1684  case OP_MUL: return "mul";
1685  case OP_DIV: return "div";
1686  case OP_NEG: return "neg";
1687  case OP_EXP: return "exp";
1688  case OP_LOG: return "log";
1689  case OP_CONSTPOW:
1690  case OP_POW: return "pow";
1691  case OP_SQRT: return "sqrt";
1692  case OP_SQ: return "sq";
1693  case OP_TWICE: return "twice";
1694  case OP_SIN: return "sin";
1695  case OP_COS: return "cos";
1696  case OP_TAN: return "tan";
1697  case OP_ASIN: return "asin";
1698  case OP_ACOS: return "acos";
1699  case OP_ATAN: return "atan";
1700  case OP_LT: return "lt";
1701  case OP_LE: return "le";
1702  case OP_EQ: return "eq";
1703  case OP_NE: return "ne";
1704  case OP_NOT: return "not";
1705  case OP_AND: return "and";
1706  case OP_OR: return "or";
1707  case OP_FLOOR: return "floor";
1708  case OP_CEIL: return "ceil";
1709  case OP_FMOD: return "fmod";
1710  case OP_REMAINDER: return "remainder";
1711  case OP_FABS: return "fabs";
1712  case OP_SIGN: return "sign";
1713  case OP_COPYSIGN: return "copysign";
1714  case OP_IF_ELSE_ZERO: return "if_else_zero";
1715  case OP_ERF: return "erf";
1716  case OP_FMIN: return "fmin";
1717  case OP_FMAX: return "fmax";
1718  case OP_INV: return "inv";
1719  case OP_SINH: return "sinh";
1720  case OP_COSH: return "cosh";
1721  case OP_TANH: return "tanh";
1722  case OP_ASINH: return "asinh";
1723  case OP_ACOSH: return "acosh";
1724  case OP_ATANH: return "atanh";
1725  case OP_ATAN2: return "atan2";
1726  case OP_CONST: return "const";
1727  case OP_INPUT: return "input";
1728  case OP_OUTPUT: return "output";
1729  case OP_PARAMETER: return "parameter";
1730  case OP_CALL: return "call";
1731  case OP_MTIMES: return "mtimes";
1732  case OP_SOLVE: return "solve";
1733  case OP_TRANSPOSE: return "transpose";
1734  case OP_DETERMINANT: return "determinant";
1735  case OP_INVERSE: return "inverse";
1736  case OP_DOT: return "dot";
1737  case OP_HORZCAT: return "horzcat";
1738  case OP_VERTCAT: return "vertcat";
1739  case OP_DIAGCAT: return "diagcat";
1740  case OP_HORZSPLIT: return "horzsplit";
1741  case OP_VERTSPLIT: return "vertsplit";
1742  case OP_DIAGSPLIT: return "diagsplit";
1743  case OP_RESHAPE: return "reshape";
1744  case OP_SPARSITY_CAST: return "sparsity_cast";
1745  case OP_SUBREF: return "subref";
1746  case OP_SUBASSIGN: return "subassign";
1747  case OP_GETNONZEROS: return "getnonzeros";
1748  case OP_GETNONZEROS_PARAM: return "getnonzeros_param";
1749  case OP_ADDNONZEROS: return "addnonzeros";
1750  case OP_ADDNONZEROS_PARAM: return "addnonzeros_param";
1751  case OP_SETNONZEROS: return "setnonzeros";
1752  case OP_SETNONZEROS_PARAM: return "setnonzeros_param";
1753  case OP_PROJECT: return "project";
1754  case OP_ASSERTION: return "assertion";
1755  case OP_NORM2: return "norm2";
1756  case OP_NORM1: return "norm1";
1757  case OP_NORMINF: return "norminf";
1758  case OP_NORMF: return "normf";
1759  case OP_ERFINV: return "erfinv";
1760  case OP_PRINTME: return "printme";
1761  case OP_LIFT: return "lift";
1762  case OP_EINSTEIN: return "einstein";
1763  case OP_BSPLINE: return "bspline";
1764  case OP_CONVEXIFY: return "convexify";
1765  case OP_LOG1P: return "log1p";
1766  case OP_EXPM1: return "expm1";
1767  case OP_HYPOT: return "hypot";
1768  case OP_LOGSUMEXP: return "logsumexp";
1769  case OP_KRON: return "kron";
1770  case OP_KRON_CONTRACT: return "kron_contract";
1771  }
1772  return "<invalid-op>";
1773  }
1774 
1775  template<typename T>
1776  inline std::string casadi_math<T>::pre(unsigned char op) {
1777  switch (op) {
1778  case OP_ASSIGN: return "";
1779  case OP_ADD: return "(";
1780  case OP_SUB: return "(";
1781  case OP_MUL: return "(";
1782  case OP_DIV: return "(";
1783  case OP_NEG: return "(-";
1784  case OP_TWICE: return "(2.*";
1785  case OP_LT: return "(";
1786  case OP_LE: return "(";
1787  case OP_EQ: return "(";
1788  case OP_NE: return "(";
1789  case OP_NOT: return "(!";
1790  case OP_AND: return "(";
1791  case OP_OR: return "(";
1792  case OP_IF_ELSE_ZERO: return "(";
1793  case OP_INV: return "(1./";
1794  default: return name(op) + "(";
1795  }
1796  }
1797 
1798  template<typename T>
1799  inline std::string casadi_math<T>::sep(unsigned char op) {
1800  switch (op) {
1801  case OP_ADD: return "+";
1802  case OP_SUB: return "-";
1803  case OP_MUL: return "*";
1804  case OP_DIV: return "/";
1805  case OP_LT: return "<";
1806  case OP_LE: return "<=";
1807  case OP_EQ: return "==";
1808  case OP_NE: return "!=";
1809  case OP_AND: return "&&";
1810  case OP_OR: return "||";
1811  case OP_IF_ELSE_ZERO: return "?";
1812  default: return ",";
1813  }
1814  }
1815 
1816  template<typename T>
1817  inline std::string casadi_math<T>::post(unsigned char op) {
1818  switch (op) {
1819  case OP_ASSIGN: return "";
1820  case OP_IF_ELSE_ZERO: return ":0)";
1821  default: return ")";
1822  }
1823  }
1824 
1825 
1826  template<class T, typename SU>
1827  T common_simp_unary(casadi_int op, const T& x, casadi_int depth,
1828  SU&& gen_unary,
1829  bool unique,
1830  bool& hit) {
1831  hit = true;
1832  switch (op) {
1833  case OP_TWICE:
1834  if (x.is_op(OP_MUL) && x.dep(0).is_constant() &&
1835  x.dep(0).is_half())
1836  return x.dep(1); // 2*(0.5*x) = x
1837  else if (x.is_op(OP_MUL) && x.dep(1).is_constant() &&
1838  x.dep(1).is_half())
1839  return x.dep(0); // 2*(x*0.5) = x
1840  break;
1841  case OP_SQ:
1842  if (x.is_op(OP_SQRT))
1843  return x.dep(); // sqrt(x)^2 = x
1844  else if (x.is_op(OP_NEG))
1845  return sq(x.dep()); // (-x)^2 = x^2
1846  else if (x.is_op(OP_FABS))
1847  return sq(x.dep()); // |x|^2 = x^2
1848  break;
1849  case OP_FABS:
1850  if (x.is_nonnegative())
1851  return x;
1852  else if (x.is_op(OP_NEG))
1853  return fabs(x.dep()); // fabs(-x) = fabs(x)
1854  break;
1855  case OP_LOG:
1856  if (x.is_op(OP_EXP))
1857  return x.dep(); // log(exp(x)) = x
1858  break;
1859  case OP_INV:
1860  if (x.is_op(OP_INV))
1861  return x.dep(); // 1/(1/x) = x
1862  break;
1863  case OP_SQRT:
1864  if (x.is_op(OP_SQ))
1865  return fabs(x.dep()); // sqrt(x^2) = x
1866  break;
1867  case OP_COS:
1868  if (x.is_op(OP_NEG))
1869  return cos(x.dep()); // cos(-x) = cos(x)
1870  else if (x.is_op(OP_FABS))
1871  return cos(x.dep()); // cos(|x|) = cos(x)
1872  break;
1873  case OP_NEG:
1874  if (x.is_op(OP_NEG))
1875  return x.dep(); // -(-x) = x
1876  else if (unique && x.is_op(OP_SUB))
1877  return x.dep(1) - x.dep(0); // -(x-y) = y-x
1878  }
1879  hit = false;
1880  return 0;
1881  }
1882 
1883 
1884  // Is the expression an inequality
1885  template<class T>
1886  bool is_ineq(const T& x) {
1887  return x.is_op(OP_LE) || x.is_op(OP_LT);
1888  }
1889 
1890  template<class T, typename SU, typename SB>
1891  T common_simp_binary(casadi_int op, const T& x, const T& y, casadi_int depth,
1892  SU&& gen_unary,
1893  SB&& gen_binary,
1894  bool unique_x,
1895  bool unique_y,
1896  bool& hit) {
1897  hit = true;
1898  switch (op) {
1899  case OP_ADD:
1900  if (x.is_zero())
1901  return y;
1902  else if (y->is_zero()) // term2 is zero
1903  return x;
1904  else if (y.is_op(OP_NEG)) // x + (-y) -> x - y
1905  return x - (-y);
1906  else if (x.is_op(OP_NEG)) // (-x) + y -> y - x
1907  return y - x.dep();
1908  else if (is_equal(x, y, depth))
1909  return 2*x; // x+x -> 2*x
1910  else if (x.is_op(OP_MUL) && y.is_op(OP_MUL) &&
1911  x.dep(0).is_constant() && static_cast<double>(x.dep(0))==0.5 &&
1912  y.dep(0).is_constant() && static_cast<double>(y.dep(0))==0.5 &&
1913  is_equal(y.dep(1), x.dep(1), depth)) // 0.5x+0.5x = x
1914  return x.dep(1);
1915  else if (x.is_op(OP_DIV) && y.is_op(OP_DIV) &&
1916  x.dep(1).is_constant() && static_cast<double>(x.dep(1))==2 &&
1917  y.dep(1).is_constant() && static_cast<double>(y.dep(1))==2 &&
1918  is_equal(y.dep(0), x.dep(0), depth)) // x/2+x/2 = x
1919  return x.dep(0);
1920  else if (x.is_op(OP_MUL) && y.is_op(OP_MUL) &&
1921  x.dep(0).is_constant() && y.dep(0).is_constant() &&
1922  static_cast<double>(x.dep(0))+static_cast<double>(y.dep(0))==1 &&
1923  is_equal(y.dep(1), x.dep(1), depth)
1924  )
1925  return x.dep(1);
1926  else if (x.is_op(OP_SUB) && is_equal(x.dep(1), y, depth))
1927  return x.dep(0);
1928  else if (y.is_op(OP_SUB) && is_equal(x, y.dep(1), depth))
1929  return y.dep(0);
1930  else if (x.is_op(OP_SQ) && y.is_op(OP_SQ) &&
1931  ((x.dep().is_op(OP_SIN) && y.dep().is_op(OP_COS))
1932  || (x.dep().is_op(OP_COS) && y.dep().is_op(OP_SIN)))
1933  && is_equal(x.dep().dep(), y.dep().dep(), depth))
1934  return 1; // sin^2 + cos^2 -> 1
1935  else if (x.is_doubled() && is_equal(x.dep(0), y, depth))
1936  return 3*x.dep(0);
1937  else if (y.is_doubled() && is_equal(y.dep(0), x, depth))
1938  return 3*y.dep(0);
1939  else if (x.is_op(OP_MUL) && y.is_op(OP_MUL) &&
1940  is_equal(x.dep(0), y.dep(0), depth) && (x.dep(1)+y.dep(1)).is_one())
1941  return x.dep(0);
1942  break;
1943  case OP_SUB:
1944  if (y.is_zero()) // term2 is zero
1945  return x;
1946  if (x.is_zero()) // term1 is zero
1947  return -y;
1948  if (is_equal(x, y, depth)) // the terms are equal
1949  return 0;
1950  else if (y.is_op(OP_NEG)) // x - (-y) -> x + y
1951  return x + y.dep();
1952  else if (x.is_op(OP_ADD) && is_equal(x.dep(1), y, depth))
1953  return x.dep(0);
1954  else if (x.is_op(OP_ADD) && is_equal(x.dep(0), y, depth))
1955  return x.dep(1);
1956  else if (y.is_op(OP_ADD) && is_equal(x, y.dep(1), depth))
1957  return -y.dep(0);
1958  else if (y.is_op(OP_ADD) && is_equal(x, y.dep(0), depth))
1959  return -y.dep(1);
1960  else if (x.is_op(OP_NEG))
1961  return -(x.dep() + y);
1962  else if (x.is_op(OP_MUL) && y.is_op(OP_MUL) &&
1963  x.dep(0).is_constant() && y.dep(0).is_constant() &&
1964  static_cast<double>(x.dep(0))-static_cast<double>(y.dep(0))==1 &&
1965  is_equal(y.dep(1), x.dep(1), depth)
1966  )
1967  return x.dep(1);
1968  else if (x.is_doubled() && is_equal(x.dep(0), y, depth))
1969  return y;
1970  else if (y.is_doubled() && is_equal(y.dep(0), x, depth))
1971  return -x;
1972  break;
1973  case OP_MUL:
1974  if (is_equal(y, x, depth))
1975  return sq(x);
1976  else if (!x.is_constant() && y.is_constant())
1977  return y * x;
1978  else if (x.is_zero() || y->is_zero()) // one of the terms is zero
1979  return 0;
1980  else if (x.is_one()) // term1 is one
1981  return y;
1982  else if (y.is_one()) // term2 is one
1983  return x;
1984  else if (y.is_minus_one())
1985  return -x;
1986  else if (x.is_minus_one())
1987  return -y;
1988  else if (y.is_op(OP_INV))
1989  return x/y.inv();
1990  else if (x.is_op(OP_INV))
1991  return y / x.inv();
1992  else if (x.is_constant() && y.is_op(OP_MUL) && y.dep(0).is_constant() &&
1993  static_cast<double>(x)*static_cast<double>(y.dep(0))==1) // 5*(0.2*x) = x
1994  return y.dep(1);
1995  else if (x.is_constant() && y.is_doubled() &&
1996  static_cast<double>(x)*2==1) // 0.5*(2*x) = x
1997  return y.dep(0);
1998  else if (x.is_constant() && y.is_op(OP_DIV) && y.dep(1).is_constant() &&
1999  static_cast<double>(x)==static_cast<double>(y.dep(1))) // 5*(x/5) = x
2000  return y.dep(0);
2001  else if (x.is_constant() && static_cast<double>(x)==2)
2002  return gen_unary(OP_TWICE, y);
2003  else if (y.is_constant() && static_cast<double>(y)==2)
2004  return gen_unary(OP_TWICE, x);
2005  else if (x.is_op(OP_DIV) && is_equal(x.dep(1), y, depth)) // ((2/x)*x)
2006  return x.dep(0);
2007  else if (y.is_op(OP_DIV) &&
2008  is_equal(y.dep(1), x, depth)) // ((2/x)*x)
2009  return y.dep(0);
2010  else if (x.is_op(OP_NEG))
2011  return -(x.dep() * y);
2012  else if (y.is_op(OP_NEG))
2013  return -(x * y.dep());
2014  else if (unique_x && x.is_op(OP_ADD) && x.dep(1).is_op(OP_DIV) &&
2015  is_equal(x.dep(1).dep(1), y, depth))
2016  return x.dep(0)*y + x.dep(1).dep(0); // (a + b/c)*c -> a*c + b
2017  else if (unique_x && x.is_op(OP_ADD) && x.dep(0).is_op(OP_DIV) &&
2018  is_equal(x.dep(0).dep(1), y, depth))
2019  return x.dep(0).dep(0)+x.dep(1)*y; // (a/c + b)*c -> a*c + b
2020  else if (unique_y && y.is_op(OP_ADD) && y.dep(1).is_op(OP_DIV) &&
2021  is_equal(y.dep(1).dep(1), x, depth))
2022  return x*y.dep(0) + y.dep(1).dep(0); // c*(a + b/c) -> c*a + b
2023  else if (unique_y && y.is_op(OP_ADD) && y.dep(0).is_op(OP_DIV) &&
2024  is_equal(y.dep(0).dep(1), x, depth))
2025  return y.dep(0).dep(0) + x*y.dep(1); // c*(a/c + b) -> a + c*b
2026  else if (unique_x && x.is_op(OP_SUB) && x.dep(1).is_op(OP_DIV) &&
2027  is_equal(x.dep(1).dep(1), y, depth))
2028  return x.dep(0)*y - x.dep(1).dep(0); // (a - b/c)*c -> a*c - b
2029  else if (unique_x && x.is_op(OP_SUB) && x.dep(0).is_op(OP_DIV) &&
2030  is_equal(x.dep(0).dep(1), y, depth))
2031  return x.dep(0).dep(0)-x.dep(1)*y; // (a/c - b)*c -> a*c - b
2032  else if (unique_y && y.is_op(OP_SUB) && y.dep(1).is_op(OP_DIV) &&
2033  is_equal(y.dep(1).dep(1), x, depth))
2034  return x*y.dep(0) - y.dep(1).dep(0); // c*(a - b/c) -> c*a - b
2035  else if (unique_y && y.is_op(OP_SUB) && y.dep(0).is_op(OP_DIV) &&
2036  is_equal(y.dep(0).dep(1), x, depth))
2037  return y.dep(0).dep(0) - x*y.dep(1); // c*(a/c - b) -> a - c*b
2038  else if (x.is_constant() && y.is_op(OP_MUL) && y.dep(0).is_constant())
2039  return x*y.dep(0)*y.dep(1); // a*(b*x) -> (a*b)*x
2040  else if (y.is_constant() && x.is_op(OP_MUL) && x.dep(0).is_constant())
2041  return y*x.dep(0)*x.dep(1); // (a*x)*b -> (a*b)*x
2042  break;
2043  case OP_DIV:
2044  if (y.is_zero()) // term2 is zero
2045  return std::numeric_limits<double>::infinity();
2046  else if (x.is_zero()) // term1 is zero
2047  return 0;
2048  else if (y.is_one()) // term2 is one
2049  return x;
2050  else if (y.is_minus_one())
2051  return -x;
2052  else if (y.is_half())
2053  return 2*x;
2054  else if (is_equal(x, y, depth)) // terms are equal
2055  return 1;
2056  else if (x.is_doubled() && y.is_constant() && static_cast<double>(y)==2)
2057  return x.dep(0);
2058  else if (x.is_doubled() && is_equal(y, x.dep(0), depth))
2059  return 2;
2060  else if (x.is_op(OP_MUL) && is_equal(y, x.dep(0), depth))
2061  return x.dep(1);
2062  else if (x.is_op(OP_MUL) && is_equal(y, x.dep(1), depth))
2063  return x.dep(0);
2064  else if (x.is_one())
2065  return y.inv();
2066  else if (y.is_op(OP_INV))
2067  return x*y.inv();
2068  else if (x.is_doubled() && y.is_doubled())
2069  return x.dep(0) / y->dep(0);
2070  else if (y.is_constant() && x.is_op(OP_DIV) && x.dep(1).is_constant() &&
2071  static_cast<double>(y)*static_cast<double>(x.dep(1))==1) // (x/5)/0.2
2072  return x.dep(0);
2073  else if (y.is_op(OP_MUL) &&
2074  is_equal(y.dep(1), x, depth)) // x/(2*x) = 1/2
2075  return 1/y.dep(0); //BinarySX::create(OP_DIV, 1, y.dep(0));
2076  else if (y.is_doubled() &&
2077  is_equal(y.dep(0), x, depth)) // x/(2*x) = 1/2
2078  return 0.5;
2079  else if (x.is_op(OP_NEG) &&
2080  is_equal(x.dep(0), y, depth)) // (-x)/x = -1
2081  return -1;
2082  else if (y.is_op(OP_NEG) &&
2083  is_equal(y.dep(0), x, depth)) // x/(-x) = 1
2084  return -1;
2085  else if (y.is_op(OP_NEG) && x.is_op(OP_NEG) &&
2086  is_equal(x.dep(0), y.dep(0), depth)) // (-x)/(-x) = 1
2087  return 1;
2088  else if (x.is_op(OP_DIV) && is_equal(y, x.dep(0), depth))
2089  return x.dep(1).inv();
2090  else if (x.is_op(OP_NEG))
2091  return -(x.dep() / y);
2092  else if (y.is_op(OP_NEG))
2093  return -(x / y.dep());
2094  break;
2095  case OP_POW:
2096  if (y.is_constant()) {
2097  if (y.is_integer()) {
2098  casadi_int nn = y->to_int();
2099  if (nn == 0) {
2100  return 1;
2101  } else if (nn>100 || nn<-100) { // maximum depth
2102  return gen_binary(OP_CONSTPOW, x, nn);
2103  } else if (nn<0) { // negative power
2104  return 1/pow(x, -nn);
2105  } else if (nn%2 == 1) { // odd power
2106  return x*pow(x, nn-1);
2107  } else { // even power
2108  T rt = pow(x, static_cast<casadi_int>(nn/2));
2109  return rt*rt;
2110  }
2111  } else if (y->to_double()==0.5) {
2112  return sqrt(x);
2113  } else {
2114  return gen_binary(OP_CONSTPOW, x, y);
2115  }
2116  }
2117  break;
2118  case OP_LE:
2119  // Keep chained inequalities intact unless operands are unique
2120  if (!unique_x && !unique_y && (is_ineq(x) || is_ineq(y))) break;
2121  if ((y-x).is_nonnegative())
2122  return 1;
2123  break;
2124  case OP_FMIN:
2125  if (x.is_inf()) return y;
2126  if (y.is_inf()) return x;
2127  if (x.is_minus_inf() || y.is_minus_inf()) return -std::numeric_limits<double>::infinity();
2128  if (is_equal(x, y, depth)) return x;
2129  break;
2130  case OP_FMAX:
2131  if (x.is_minus_inf()) return y;
2132  if (y.is_minus_inf()) return x;
2133  if (x.is_inf() || y.is_inf()) return std::numeric_limits<double>::infinity();
2134  if (is_equal(x, y, depth)) return x;
2135  break;
2136  case OP_LT:
2137  if (!unique_x && !unique_y && (is_ineq(x) || is_ineq(y))) break;
2138  if (((x)-y).is_nonnegative())
2139  return 0;
2140  break;
2141  case OP_EQ:
2142  if (is_equal(x, y, depth))
2143  return 1;
2144  break;
2145  case OP_NE:
2146  if (is_equal(x, y, depth))
2147  return 0;
2148  break;
2149  case OP_IF_ELSE_ZERO:
2150  if (y.is_zero()) {
2151  return y;
2152  } else if (x.is_constant()) {
2153  if (static_cast<double>(x)!=0) {
2154  return y;
2155  } else {
2156  return 0;
2157  }
2158  }
2159  }
2160  hit = false;
2161  return 0;
2162  }
2163 
2164 #endif // SWIG
2165 
2166 } // namespace casadi
2167 
2169 
2170 #endif // CASADI_CALCULUS_HPP
casadi::rosenbrock_tableau::b
const long double b[8]
Definition: rosenbrock_tableau.hpp:109
casadi
The casadi namespace.
Definition: archiver.hpp:32
casadi::uout
CASADI_EXPORT std::ostream & uout()