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;};
553 
556  template<casadi_int I> struct NonnegativeChecker { static const bool check=false;};
557  template<> struct NonnegativeChecker<OP_SQRT>{ static const bool check=true;};
558  template<> struct NonnegativeChecker<OP_SQ>{ static const bool check=true;};
559  template<> struct NonnegativeChecker<OP_EXP>{ static const bool check=true;};
560  template<> struct NonnegativeChecker<OP_LT>{ static const bool check=true;};
561  template<> struct NonnegativeChecker<OP_LE>{ static const bool check=true;};
562  template<> struct NonnegativeChecker<OP_EQ>{ static const bool check=true;};
563  template<> struct NonnegativeChecker<OP_NE>{ static const bool check=true;};
564  template<> struct NonnegativeChecker<OP_NOT>{ static const bool check=true;};
565  template<> struct NonnegativeChecker<OP_AND>{ static const bool check=true;};
566  template<> struct NonnegativeChecker<OP_OR>{ static const bool check=true;};
567  template<> struct NonnegativeChecker<OP_HYPOT>{ static const bool check=true;};
568  template<> struct NonnegativeChecker<OP_FABS>{ static const bool check=true;};
570 
573  template<casadi_int I> struct NargChecker { static const casadi_int check=1;};
574  template<> struct NargChecker<OP_ADD>{ static const casadi_int check=2;};
575  template<> struct NargChecker<OP_SUB>{ static const casadi_int check=2;};
576  template<> struct NargChecker<OP_MUL>{ static const casadi_int check=2;};
577  template<> struct NargChecker<OP_DIV>{ static const casadi_int check=2;};
578  template<> struct NargChecker<OP_POW>{ static const casadi_int check=2;};
579  template<> struct NargChecker<OP_CONSTPOW>{ static const casadi_int check=2;};
580  template<> struct NargChecker<OP_EQ>{ static const casadi_int check=2;};
581  template<> struct NargChecker<OP_NE>{ static const casadi_int check=2;};
582  template<> struct NargChecker<OP_AND>{ static const casadi_int check=2;};
583  template<> struct NargChecker<OP_OR>{ static const casadi_int check=2;};
584  template<> struct NargChecker<OP_FMIN>{ static const casadi_int check=2;};
585  template<> struct NargChecker<OP_FMAX>{ static const casadi_int check=2;};
586  template<> struct NargChecker<OP_PRINTME>{ static const casadi_int check=2;};
587  template<> struct NargChecker<OP_ATAN2>{ static const casadi_int check=2;};
588  template<> struct NargChecker<OP_IF_ELSE_ZERO>{ static const casadi_int check=2;};
589  template<> struct NargChecker<OP_FMOD>{ static const casadi_int check=2;};
590  template<> struct NargChecker<OP_REMAINDER>{ static const casadi_int check=2;};
591  template<> struct NargChecker<OP_COPYSIGN>{ static const casadi_int check=2;};
592  template<> struct NargChecker<OP_CONST>{ static const casadi_int check=0;};
593  template<> struct NargChecker<OP_PARAMETER>{ static const casadi_int check=0;};
594  template<> struct NargChecker<OP_INPUT>{ static const casadi_int check=0;};
595  template<> struct NargChecker<OP_HYPOT>{ static const casadi_int check=2;};
597 
599  template<>
600  struct UnaryOperation<OP_ASSIGN>{
601  public:
602  template<typename T> static inline void fcn(const T& x, T& f) { f = x;}
603  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 1; }
604  };
605 
607  template<>
608  struct BinaryOperation<OP_ADD>{
609  public:
610  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x+y;}
611  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
612  d[0]=d[1]=1;}
613  };
614 
616  template<>
617  struct BinaryOperation<OP_SUB>{
618  public:
619  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x-y;}
620  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
621  d[0]=1; d[1]=-1;}
622  };
623 
625  template<>
626  struct BinaryOperation<OP_MUL>{
627  public:
628  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x*y;}
629  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
630  d[0]=y; d[1]=x;}
631  };
632 
634  template<>
635  struct BinaryOperation<OP_DIV>{
636  public:
637  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x/y;}
638  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
639  d[0]=1/y; d[1]=-f/y;}
640  };
641 
643  template<>
644  struct UnaryOperation<OP_NEG>{
645  public:
646  template<typename T> static inline void fcn(const T& x, T& f) { f = -x;}
647  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=-1;}
648  };
649 
651  template<>
652  struct UnaryOperation<OP_EXP>{
653  public:
654  template<typename T> static inline void fcn(const T& x, T& f) { f = exp(x);}
655  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=f;}
656  };
657 
659  template<>
660  struct UnaryOperation<OP_LOG>{
661  public:
662  template<typename T> static inline void fcn(const T& x, T& f) { f = log(x);}
663  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=1/x;}
664  };
665 
667  template<>
668  struct BinaryOperation<OP_POW>{
669  public:
670  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = pow(x, y);}
671  // See issue #104 why d[0] is no longer y*f/x
672  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
673  d[0]=y*pow(x, y-1); d[1]=log(x)*f;}
674  };
675 
677  template<>
678  struct BinaryOperation<OP_CONSTPOW>{
679  public:
680  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = pow(x, y);}
681  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
682  d[0]=y*pow(x, y-1); d[1]=0;}
683  };
684 
686  template<>
687  struct UnaryOperation<OP_SQRT>{
688  public:
689  template<typename T> static inline void fcn(const T& x, T& f) { f = sqrt(x);}
690  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=1/(twice(f));}
691  };
692 
694  template<>
695  struct UnaryOperation<OP_SQ>{
696  public:
697  template<typename T> static inline void fcn(const T& x, T& f) { f = sq(x);}
698  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=twice(x);}
699  };
700 
702  template<>
703  struct UnaryOperation<OP_TWICE>{
704  template<typename T> static inline void fcn(const T& x, T& f) { f = 2.*x;}
705  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 2; }
706  };
707 
709  template<>
710  struct UnaryOperation<OP_SIN>{
711  public:
712  template<typename T> static inline void fcn(const T& x, T& f) { f = sin(x);}
713  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=cos(x);}
714  };
715 
717  template<>
718  struct UnaryOperation<OP_COS>{
719  public:
720  template<typename T> static inline void fcn(const T& x, T& f) { f = cos(x);}
721  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=-sin(x);}
722  };
723 
725  template<>
726  struct UnaryOperation<OP_TAN>{
727  public:
728  template<typename T> static inline void fcn(const T& x, T& f) { f = tan(x);}
729  template<typename T> static inline void der(const T& x, const T& f, T* d)
730  { d[0] = 1/sq(cos(x));}
731  };
732 
734  template<>
735  struct UnaryOperation<OP_ASIN>{
736  public:
737  template<typename T> static inline void fcn(const T& x, T& f) { f = asin(x);}
738  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=1/sqrt(1-x*x);}
739  };
740 
742  template<>
743  struct UnaryOperation<OP_ACOS>{
744  public:
745  template<typename T> static inline void fcn(const T& x, T& f) { f = acos(x);}
746  template<typename T> static inline void der(const T& x, const T& f, T* d)
747  { d[0]=-1/sqrt(1-x*x);}
748  };
749 
751  template<>
752  struct UnaryOperation<OP_ATAN>{
753  public:
754  template<typename T> static inline void fcn(const T& x, T& f) { f = atan(x);}
755  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 1/(1+x*x);}
756  };
757 
759  template<>
760  struct BinaryOperation<OP_LT>{
761  public:
762  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x < y;}
763  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
764  d[0]=d[1]=0;}
765  };
766 
768  template<>
769  struct BinaryOperation<OP_LE>{
770  public:
771  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x <= y;}
772  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
773  d[0]=d[1]=0;}
774  };
775 
777  template<>
778  struct UnaryOperation<OP_FLOOR>{
779  public:
780  template<typename T> static inline void fcn(const T& x, T& f) { f = floor(x);}
781  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 0;}
782  };
783 
785  template<>
786  struct UnaryOperation<OP_CEIL>{
787  public:
788  template<typename T> static inline void fcn(const T& x, T& f) { f = ceil(x);}
789  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 0;}
790  };
791 
793  template<>
794  struct BinaryOperation<OP_FMOD>{
795  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = fmod(x, y);}
796  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
797  d[0]=1; d[1]=(f-x)/y;}
798  };
799 
801  template<>
802  struct BinaryOperation<OP_REMAINDER>{
803  template<typename T> static inline void fcn(const T& x, const T& y, T& f) {
804  f = remainder(x, y);}
805  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
806  d[0]=1; d[1]=(f-x)/y;}
807  };
808 
810  template<>
811  struct BinaryOperation<OP_EQ>{
812  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x==y;}
813  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
814  d[0]=d[1]=0;}
815  };
816 
818  template<>
819  struct BinaryOperation<OP_NE>{
820  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x!=y;}
821  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
822  d[0]=d[1]=0;}
823  };
824 
826  template<>
827  struct UnaryOperation<OP_NOT>{
828  public:
829  template<typename T> static inline void fcn(const T& x, T& f) { f = !x;}
830  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 0;}
831  };
832 
834  template<>
835  struct BinaryOperation<OP_AND>{
836  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x && y;}
837  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
838  d[0]=d[1]=0;}
839  };
840 
842  template<>
843  struct BinaryOperation<OP_OR>{
844  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x || y;}
845  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
846  d[0]=d[1]=0;}
847  };
848 
850  template<>
851  struct UnaryOperation<OP_ERF>{
852  template<typename T> static inline void fcn(const T& x, T& f) { f = erf(x);}
853  template<typename T> static inline void der(const T& x, const T& f, T* d) {
854  d[0] = (2/sqrt(pi))*exp(-x*x);}
855  };
856 
858  template<>
859  struct UnaryOperation<OP_FABS>{
860  template<typename T> static inline void fcn(const T& x, T& f) { f = fabs(x);}
861  template<typename T> static inline void der(const T& x, const T& f, T* d) {
862  d[0]=sign(x);}
863  };
864 
866  template<>
867  struct UnaryOperation<OP_SIGN>{
868  template<typename T> static inline void fcn(const T& x, T& f) { f = sign(x);}
869  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0]=0;}
870  };
871 
873  template<>
874  struct BinaryOperation<OP_COPYSIGN>{
875  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = copysign(x, y);}
876  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
877  T e = 1; d[0]=copysign(e, y); d[1]=0;}
878  };
879 
881  template<>
882  struct BinaryOperation<OP_FMIN>{
883  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = fmin(x, y);}
884  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
885  T a = x<=y;
886  T b = y<=x;
887  T c = a+b;
888  d[0]=a/c; d[1]=b/c;}
889  };
890 
892  template<>
893  struct BinaryOperation<OP_FMAX>{
894  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = fmax(x, y);}
895  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
896  T a = y<=x;
897  T b = x<=y;
898  T c = a+b;
899  d[0]=a/c; d[1]=b/c;}
900  };
901 
903  template<>
904  struct UnaryOperation<OP_INV>{
905  template<typename T> static inline void fcn(const T& x, T& f) { f = 1./x;}
906  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = -f*f; }
907  };
908 
910  template<>
911  struct UnaryOperation<OP_SINH>{
912  template<typename T> static inline void fcn(const T& x, T& f) { f = sinh(x);}
913  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = cosh(x); }
914  };
915 
917  template<>
918  struct UnaryOperation<OP_COSH>{
919  template<typename T> static inline void fcn(const T& x, T& f) { f = cosh(x);}
920  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = sinh(x); }
921  };
922 
924  template<>
925  struct UnaryOperation<OP_TANH>{
926  template<typename T> static inline void fcn(const T& x, T& f) { f = tanh(x);}
927  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 1-f*f; }
928  };
929 
931  template<>
932  struct UnaryOperation<OP_ASINH>{
933  template<typename T> static inline void fcn(const T& x, T& f) { f = asinh(x);}
934  template<typename T> static inline void der(const T& x, const T& f, T* d) {
935  d[0] = 1/sqrt(1+x*x); }
936  };
937 
939  template<>
940  struct UnaryOperation<OP_ACOSH>{
941  template<typename T> static inline void fcn(const T& x, T& f) { f = acosh(x);}
942  template<typename T> static inline void der(const T& x, const T& f, T* d) {
943  d[0] = 1/sqrt(x-1)/sqrt(x+1); }
944  };
945 
947  template<>
948  struct UnaryOperation<OP_ATANH>{
949  template<typename T> static inline void fcn(const T& x, T& f) { f = atanh(x);}
950  template<typename T> static inline void der(const T& x, const T& f, T* d) { d[0] = 1/(1-x*x); }
951  };
952 
954  template<>
955  struct UnaryOperation<OP_ERFINV>{
956  template<typename T> static inline void fcn(const T& x, T& f) { f = erfinv(x);}
957  template<typename T> static inline void der(const T& x, const T& f, T* d) {
958  d[0] = (sqrt(pi)/2)*exp(f*f); }
959  };
960 
962  template<>
963  struct UnaryOperation<OP_LOG1P>{
964  template<typename T> static inline void fcn(const T& x, T& f) { f = log1p(x);}
965  template<typename T> static inline void der(const T& x, const T& f, T* d) {
966  d[0] = 1/(1+x);}
967  };
968 
970  template<>
971  struct UnaryOperation<OP_EXPM1>{
972  template<typename T> static inline void fcn(const T& x, T& f) { f = expm1(x);}
973  template<typename T> static inline void der(const T& x, const T& f, T* d) {
974  d[0] = exp(x); }
975  };
976 
978  template<>
979  struct BinaryOperation<OP_PRINTME>{
980  template<typename T> static inline void fcn(const T& x, const T& y, T& f) {f = printme(x, y); }
981  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
982  d[0]=1; d[1]=0;}
983  };
984 
986  template<>
987  struct BinaryOperation<OP_ATAN2>{
988  public:
989  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = atan2(x, y);}
990  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
991  T t = x*x+y*y; d[0]=y/t; d[1]=-x/t;}
992  };
993 
995  template<>
996  struct BinaryOperation<OP_IF_ELSE_ZERO>{
997  public:
998  template<typename T> static inline void fcn(const T& x, const T& y, T& f) {
999  f = if_else_zero(x, y);}
1000  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
1001  d[0]=0; d[1]=x;}
1002  };
1003 
1005  template<>
1006  struct BinaryOperation<OP_LIFT>{
1007  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = x;}
1008  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
1009  d[0] = 1; d[1] = 0; }
1010  };
1011 
1013  template<>
1014  struct BinaryOperation<OP_HYPOT>{
1015  template<typename T> static inline void fcn(const T& x, const T& y, T& f) { f = hypot(x, y);}
1016  template<typename T> static inline void der(const T& x, const T& y, const T& f, T* d) {
1017  d[0] = x/f; d[1] = y/f; }
1018  };
1019 
1020  template<template<casadi_int> class F, typename T>
1021  T operation_getter(casadi_int op) {
1022  switch (static_cast<Operation>(op)) {
1023  case OP_ASSIGN: return F<OP_ASSIGN>::check;
1024  case OP_ADD: return F<OP_ADD>::check;
1025  case OP_SUB: return F<OP_SUB>::check;
1026  case OP_MUL: return F<OP_MUL>::check;
1027  case OP_DIV: return F<OP_DIV>::check;
1028  case OP_NEG: return F<OP_NEG>::check;
1029  case OP_EXP: return F<OP_EXP>::check;
1030  case OP_LOG: return F<OP_LOG>::check;
1031  case OP_POW: return F<OP_POW>::check;
1032  case OP_CONSTPOW: return F<OP_CONSTPOW>::check;
1033  case OP_SQRT: return F<OP_SQRT>::check;
1034  case OP_SQ: return F<OP_SQ>::check;
1035  case OP_TWICE: return F<OP_TWICE>::check;
1036  case OP_SIN: return F<OP_SIN>::check;
1037  case OP_COS: return F<OP_COS>::check;
1038  case OP_TAN: return F<OP_TAN>::check;
1039  case OP_ASIN: return F<OP_ASIN>::check;
1040  case OP_ACOS: return F<OP_ACOS>::check;
1041  case OP_ATAN: return F<OP_ATAN>::check;
1042  case OP_LT: return F<OP_LT>::check;
1043  case OP_LE: return F<OP_LE>::check;
1044  case OP_EQ: return F<OP_EQ>::check;
1045  case OP_NE: return F<OP_NE>::check;
1046  case OP_NOT: return F<OP_NOT>::check;
1047  case OP_AND: return F<OP_AND>::check;
1048  case OP_OR: return F<OP_OR>::check;
1049  case OP_FLOOR: return F<OP_FLOOR>::check;
1050  case OP_CEIL: return F<OP_CEIL>::check;
1051  case OP_FMOD: return F<OP_FMOD>::check;
1052  case OP_REMAINDER: return F<OP_REMAINDER>::check;
1053  case OP_FABS: return F<OP_FABS>::check;
1054  case OP_SIGN: return F<OP_SIGN>::check;
1055  case OP_COPYSIGN: return F<OP_COPYSIGN>::check;
1056  case OP_IF_ELSE_ZERO: return F<OP_IF_ELSE_ZERO>::check;
1057  case OP_ERF: return F<OP_ERF>::check;
1058  case OP_FMIN: return F<OP_FMIN>::check;
1059  case OP_FMAX: return F<OP_FMAX>::check;
1060  case OP_INV: return F<OP_INV>::check;
1061  case OP_SINH: return F<OP_SINH>::check;
1062  case OP_COSH: return F<OP_COSH>::check;
1063  case OP_TANH: return F<OP_TANH>::check;
1064  case OP_ASINH: return F<OP_ASINH>::check;
1065  case OP_ACOSH: return F<OP_ACOSH>::check;
1066  case OP_ATANH: return F<OP_ATANH>::check;
1067  case OP_ATAN2: return F<OP_ATAN2>::check;
1068  case OP_CONST: return F<OP_CONST>::check;
1069  case OP_INPUT: return F<OP_INPUT>::check;
1070  case OP_OUTPUT: return F<OP_OUTPUT>::check;
1071  case OP_PARAMETER: return F<OP_PARAMETER>::check;
1072  case OP_CALL: return F<OP_CALL>::check;
1073  case OP_FIND: return F<OP_FIND>::check;
1074  case OP_LOW: return F<OP_LOW>::check;
1075  case OP_MAP: return F<OP_MAP>::check;
1076  case OP_MTIMES: return F<OP_MTIMES>::check;
1077  case OP_SOLVE: return F<OP_SOLVE>::check;
1078  case OP_TRANSPOSE: return F<OP_TRANSPOSE>::check;
1079  case OP_DETERMINANT: return F<OP_DETERMINANT>::check;
1080  case OP_INVERSE: return F<OP_INVERSE>::check;
1081  case OP_DOT: return F<OP_DOT>::check;
1082  case OP_BILIN: return F<OP_BILIN>::check;
1083  case OP_RANK1: return F<OP_RANK1>::check;
1084  case OP_HORZCAT: return F<OP_HORZCAT>::check;
1085  case OP_VERTCAT: return F<OP_VERTCAT>::check;
1086  case OP_DIAGCAT: return F<OP_DIAGCAT>::check;
1087  case OP_HORZSPLIT: return F<OP_HORZSPLIT>::check;
1088  case OP_VERTSPLIT: return F<OP_VERTSPLIT>::check;
1089  case OP_DIAGSPLIT: return F<OP_DIAGSPLIT>::check;
1090  case OP_RESHAPE: return F<OP_RESHAPE>::check;
1091  case OP_SPARSITY_CAST: return F<OP_SPARSITY_CAST>::check;
1092  case OP_SUBREF: return F<OP_SUBREF>::check;
1093  case OP_SUBASSIGN: return F<OP_SUBASSIGN>::check;
1094  case OP_GETNONZEROS: return F<OP_GETNONZEROS>::check;
1095  case OP_GETNONZEROS_PARAM: return F<OP_GETNONZEROS_PARAM>::check;
1096  case OP_ADDNONZEROS: return F<OP_ADDNONZEROS>::check;
1097  case OP_ADDNONZEROS_PARAM: return F<OP_ADDNONZEROS>::check;
1098  case OP_SETNONZEROS: return F<OP_SETNONZEROS>::check;
1099  case OP_SETNONZEROS_PARAM: return F<OP_SETNONZEROS>::check;
1100  case OP_PROJECT: return F<OP_PROJECT>::check;
1101  case OP_ASSERTION: return F<OP_ASSERTION>::check;
1102  case OP_MONITOR: return F<OP_MONITOR>::check;
1103  case OP_DUMP: return F<OP_DUMP>::check;
1104  case OP_NORM2: return F<OP_NORM2>::check;
1105  case OP_NORM1: return F<OP_NORM1>::check;
1106  case OP_NORMINF: return F<OP_NORMINF>::check;
1107  case OP_NORMF: return F<OP_NORMF>::check;
1108  case OP_MMIN: return F<OP_MMIN>::check;
1109  case OP_MMAX: return F<OP_MMAX>::check;
1110  case OP_HORZREPMAT: return F<OP_HORZREPMAT>::check;
1111  case OP_HORZREPSUM: return F<OP_HORZREPSUM>::check;
1112  case OP_ERFINV: return F<OP_ERFINV>::check;
1113  case OP_PRINTME: return F<OP_PRINTME>::check;
1114  case OP_LIFT: return F<OP_LIFT>::check;
1115  case OP_EINSTEIN: return F<OP_EINSTEIN>::check;
1116  case OP_BSPLINE: return F<OP_BSPLINE>::check;
1117  case OP_CONVEXIFY: return F<OP_CONVEXIFY>::check;
1118  case OP_LOG1P: return F<OP_LOG1P>::check;
1119  case OP_EXPM1: return F<OP_EXPM1>::check;
1120  case OP_HYPOT: return F<OP_HYPOT>::check;
1121  case OP_LOGSUMEXP: return F<OP_LOGSUMEXP>::check;
1122  case OP_KRON: return F<OP_KRON>::check;
1123  case OP_KRON_CONTRACT: return F<OP_KRON_CONTRACT>::check;
1124  }
1125  return T();
1126  }
1127 
1128  template<template<casadi_int> class F>
1129  bool operation_checker(casadi_int op) {
1130  return operation_getter<F, bool>(op);
1131  }
1132 
1134  template<typename T>
1135  struct casadi_math {
1136 
1140  static inline void fun(unsigned char op, const T& x, const T& y, T& f);
1141 
1145  static inline void fun(unsigned char op, const T* x, const T* y, T* f, casadi_int n);
1146 
1150  static inline void fun(unsigned char op, const T* x, const T& y, T* f, casadi_int n);
1151 
1155  static inline void fun(unsigned char op, const T& x, const T* y, T* f, casadi_int n);
1156 
1160  static inline void der(unsigned char op, const T& x, const T& y, const T& f, T* d);
1161 
1165  static inline void derF(unsigned char op, const T& x, const T& y, T& f, T* d);
1166 
1170  static inline void fun_linear(unsigned char op, const T*x, const T* y, T* f);
1171 
1175  static inline bool is_binary(unsigned char op);
1176 
1180  static inline bool is_unary(unsigned char op);
1181 
1185  static inline casadi_int ndeps(unsigned char op);
1186 
1190  static inline std::string print(unsigned char op, const std::string& x,
1191  const std::string& y);
1192  static inline std::string print(unsigned char op, const std::string& x);
1193  static inline std::string name(unsigned char op);
1194  static inline std::string pre(unsigned char op);
1195  static inline std::string sep(unsigned char op);
1196  static inline std::string post(unsigned char op);
1197  };
1198 
1200  template<>
1201  struct casadi_math<casadi_int>{
1202 
1206  static inline void fun(unsigned char op, const casadi_int& x,
1207  const casadi_int& y, casadi_int& f) {
1208  double ff(0);
1209  casadi_math<double>::fun(op, static_cast<double>(x), static_cast<double>(y), ff);
1210  f = static_cast<casadi_int>(ff);
1211  }
1212 
1213  static inline void fun(unsigned char op, const casadi_int* x, const casadi_int* y,
1214  casadi_int* f, casadi_int n) {
1215  for (casadi_int i=0; i<n; ++i) {
1216  double ff(0);
1217  casadi_math<double>::fun(op, static_cast<double>(*x++), static_cast<double>(*y++), ff);
1218  *f++ = static_cast<casadi_int>(ff);
1219  }
1220  }
1221 
1222  static inline void fun(unsigned char op, const casadi_int* x, const casadi_int& y,
1223  casadi_int* f, casadi_int n) {
1224  for (casadi_int i=0; i<n; ++i) {
1225  double ff;
1226  casadi_math<double>::fun(op, static_cast<double>(*x++), static_cast<double>(y), ff);
1227  *f++ = static_cast<casadi_int>(ff);
1228  }
1229  }
1230 
1231  static inline void fun(unsigned char op, const casadi_int& x, const casadi_int* y,
1232  casadi_int* f, casadi_int n) {
1233  for (casadi_int i=0; i<n; ++i) {
1234  double ff;
1235  casadi_math<double>::fun(op, static_cast<double>(x), static_cast<double>(*y++), ff);
1236  *f++ = static_cast<casadi_int>(ff);
1237  }
1238  }
1239 
1243  static inline void der(unsigned char op, const casadi_int& x, const casadi_int& y,
1244  const casadi_int& f, casadi_int* d) {
1245  double d_real[2] = {static_cast<double>(d[0]), static_cast<double>(d[1])};
1246  casadi_math<double>::der(op, static_cast<double>(x), static_cast<double>(y),
1247  static_cast<double>(f), d_real);
1248  d[0] = static_cast<casadi_int>(d_real[0]);
1249  d[1] = static_cast<casadi_int>(d_real[1]);
1250  }
1251 
1255  static inline void derF(unsigned char op, const casadi_int& x, const casadi_int& y,
1256  casadi_int& f, casadi_int* d) {
1257  double d_real[2] = {static_cast<double>(d[0]), static_cast<double>(d[1])};
1258  double f_real = static_cast<double>(f);
1259  casadi_math<double>::derF(op, static_cast<double>(x), static_cast<double>(y), f_real, d_real);
1260  f = static_cast<casadi_int>(f_real);
1261  d[0] = static_cast<casadi_int>(d_real[0]);
1262  d[1] = static_cast<casadi_int>(d_real[1]);
1263  }
1264 
1268  static inline casadi_int ndeps(unsigned char op) {
1269  return casadi_math<double>::ndeps(op);
1270  }
1271 
1275  static inline std::string print(unsigned char op, const std::string& x,
1276  const std::string& y) {
1277  return casadi_math<double>::print(op, x, y);
1278  }
1279  static inline std::string print(unsigned char op, const std::string& x) {
1280  return casadi_math<double>::print(op, x);
1281  }
1282  static inline std::string pre(unsigned char op) {
1283  return casadi_math<double>::pre(op);
1284  }
1285  static inline std::string name(unsigned char op) {
1286  return casadi_math<double>::name(op);
1287  }
1288  static inline std::string sep(unsigned char op) {
1289  return casadi_math<double>::sep(op);
1290  }
1291  static inline std::string post(unsigned char op) {
1292  return casadi_math<double>::post(op);
1293  }
1294  };
1295 
1296  // Template implementations
1297 
1298  template<typename T>
1299  inline void casadi_math<T>::fun(unsigned char op, const T& x, const T& y, T& f) {
1300  // NOTE: We define the implementation in a preprocessor macro to be able to force inlining,
1301  // and to allow extensions in the VM
1302 #define CASADI_MATH_FUN_BUILTIN_GEN(CNAME, X, Y, F, N) \
1303  case OP_ASSIGN: CNAME<OP_ASSIGN>::fcn(X, Y, F, N); break; \
1304  case OP_ADD: CNAME<OP_ADD>::fcn(X, Y, F, N); break; \
1305  case OP_SUB: CNAME<OP_SUB>::fcn(X, Y, F, N); break; \
1306  case OP_MUL: CNAME<OP_MUL>::fcn(X, Y, F, N); break; \
1307  case OP_DIV: CNAME<OP_DIV>::fcn(X, Y, F, N); break; \
1308  case OP_NEG: CNAME<OP_NEG>::fcn(X, Y, F, N); break; \
1309  case OP_EXP: CNAME<OP_EXP>::fcn(X, Y, F, N); break; \
1310  case OP_LOG: CNAME<OP_LOG>::fcn(X, Y, F, N); break; \
1311  case OP_POW: CNAME<OP_POW>::fcn(X, Y, F, N); break; \
1312  case OP_CONSTPOW: CNAME<OP_CONSTPOW>::fcn(X, Y, F, N); break; \
1313  case OP_SQRT: CNAME<OP_SQRT>::fcn(X, Y, F, N); break; \
1314  case OP_SQ: CNAME<OP_SQ>::fcn(X, Y, F, N); break; \
1315  case OP_TWICE: CNAME<OP_TWICE>::fcn(X, Y, F, N); break; \
1316  case OP_SIN: CNAME<OP_SIN>::fcn(X, Y, F, N); break; \
1317  case OP_COS: CNAME<OP_COS>::fcn(X, Y, F, N); break; \
1318  case OP_TAN: CNAME<OP_TAN>::fcn(X, Y, F, N); break; \
1319  case OP_ASIN: CNAME<OP_ASIN>::fcn(X, Y, F, N); break; \
1320  case OP_ACOS: CNAME<OP_ACOS>::fcn(X, Y, F, N); break; \
1321  case OP_ATAN: CNAME<OP_ATAN>::fcn(X, Y, F, N); break; \
1322  case OP_LT: CNAME<OP_LT>::fcn(X, Y, F, N); break; \
1323  case OP_LE: CNAME<OP_LE>::fcn(X, Y, F, N); break; \
1324  case OP_EQ: CNAME<OP_EQ>::fcn(X, Y, F, N); break; \
1325  case OP_NE: CNAME<OP_NE>::fcn(X, Y, F, N); break; \
1326  case OP_NOT: CNAME<OP_NOT>::fcn(X, Y, F, N); break; \
1327  case OP_AND: CNAME<OP_AND>::fcn(X, Y, F, N); break; \
1328  case OP_OR: CNAME<OP_OR>::fcn(X, Y, F, N); break; \
1329  case OP_IF_ELSE_ZERO: CNAME<OP_IF_ELSE_ZERO>::fcn(X, Y, F, N); break; \
1330  case OP_FLOOR: CNAME<OP_FLOOR>::fcn(X, Y, F, N); break; \
1331  case OP_CEIL: CNAME<OP_CEIL>::fcn(X, Y, F, N); break; \
1332  case OP_FMOD: CNAME<OP_FMOD>::fcn(X, Y, F, N); break; \
1333  case OP_REMAINDER: CNAME<OP_REMAINDER>::fcn(X, Y, F, N); break; \
1334  case OP_FABS: CNAME<OP_FABS>::fcn(X, Y, F, N); break; \
1335  case OP_SIGN: CNAME<OP_SIGN>::fcn(X, Y, F, N); break; \
1336  case OP_COPYSIGN: CNAME<OP_COPYSIGN>::fcn(X, Y, F, N); break; \
1337  case OP_ERF: CNAME<OP_ERF>::fcn(X, Y, F, N); break; \
1338  case OP_FMIN: CNAME<OP_FMIN>::fcn(X, Y, F, N); break; \
1339  case OP_FMAX: CNAME<OP_FMAX>::fcn(X, Y, F, N); break; \
1340  case OP_INV: CNAME<OP_INV>::fcn(X, Y, F, N); break; \
1341  case OP_SINH: CNAME<OP_SINH>::fcn(X, Y, F, N); break; \
1342  case OP_COSH: CNAME<OP_COSH>::fcn(X, Y, F, N); break; \
1343  case OP_TANH: CNAME<OP_TANH>::fcn(X, Y, F, N); break; \
1344  case OP_ASINH: CNAME<OP_ASINH>::fcn(X, Y, F, N); break; \
1345  case OP_ACOSH: CNAME<OP_ACOSH>::fcn(X, Y, F, N); break; \
1346  case OP_ATANH: CNAME<OP_ATANH>::fcn(X, Y, F, N); break; \
1347  case OP_ATAN2: CNAME<OP_ATAN2>::fcn(X, Y, F, N); break; \
1348  case OP_ERFINV: CNAME<OP_ERFINV>::fcn(X, Y, F, N); break; \
1349  case OP_LIFT: CNAME<OP_LIFT>::fcn(X, Y, F, N); break; \
1350  case OP_PRINTME: CNAME<OP_PRINTME>::fcn(X, Y, F, N); break; \
1351  case OP_LOG1P: CNAME<OP_LOG1P>::fcn(X, Y, F, N); break; \
1352  case OP_EXPM1: CNAME<OP_EXPM1>::fcn(X, Y, F, N); break; \
1353  case OP_HYPOT: CNAME<OP_HYPOT>::fcn(X, Y, F, N); break;
1354 
1355 #define CASADI_MATH_FUN_BUILTIN(X, Y, F) CASADI_MATH_FUN_BUILTIN_GEN(BinaryOperationSS, X, Y, F, 1)
1356 
1357  switch (op) {
1358  CASADI_MATH_FUN_BUILTIN(x, y, f)
1359  }
1360  }
1361 
1362  template<typename T>
1363  inline void casadi_math<T>::fun(unsigned char op, const T* x, const T* y, T* f, casadi_int n) {
1364  switch (op) {
1365  CASADI_MATH_FUN_BUILTIN_GEN(BinaryOperationVV, x, y, f, n)
1366  }
1367  }
1368 
1369  template<typename T>
1370  inline void casadi_math<T>::fun(unsigned char op, const T* x, const T& y, T* f, casadi_int n) {
1371  switch (op) {
1372  CASADI_MATH_FUN_BUILTIN_GEN(BinaryOperationVS, x, y, f, n)
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(BinaryOperationSV, x, y, f, n)
1380  }
1381  }
1382 
1383 
1384  template<typename T>
1385  inline void casadi_math<T>::der(unsigned char op, const T& x, const T& y, const T& f, T* d) {
1386  // NOTE: We define the implementation in a preprocessor macro to be able to force inlining,
1387  // and to allow extensions in the VM
1388 #define CASADI_MATH_DER_BUILTIN(X, Y, F, D) \
1389  case OP_ASSIGN: BinaryOperation<OP_ASSIGN>::der(X, Y, F, D); break; \
1390  case OP_ADD: BinaryOperation<OP_ADD>::der(X, Y, F, D); break; \
1391  case OP_SUB: BinaryOperation<OP_SUB>::der(X, Y, F, D); break; \
1392  case OP_MUL: BinaryOperation<OP_MUL>::der(X, Y, F, D); break; \
1393  case OP_DIV: BinaryOperation<OP_DIV>::der(X, Y, F, D); break; \
1394  case OP_NEG: BinaryOperation<OP_NEG>::der(X, Y, F, D); break; \
1395  case OP_EXP: BinaryOperation<OP_EXP>::der(X, Y, F, D); break; \
1396  case OP_LOG: BinaryOperation<OP_LOG>::der(X, Y, F, D); break; \
1397  case OP_POW: BinaryOperation<OP_POW>::der(X, Y, F, D); break; \
1398  case OP_CONSTPOW: BinaryOperation<OP_CONSTPOW>::der(X, Y, F, D); break; \
1399  case OP_SQRT: BinaryOperation<OP_SQRT>::der(X, Y, F, D); break; \
1400  case OP_SQ: BinaryOperation<OP_SQ>::der(X, Y, F, D); break; \
1401  case OP_TWICE: BinaryOperation<OP_TWICE>::der(X, Y, F, D); break; \
1402  case OP_SIN: BinaryOperation<OP_SIN>::der(X, Y, F, D); break; \
1403  case OP_COS: BinaryOperation<OP_COS>::der(X, Y, F, D); break; \
1404  case OP_TAN: BinaryOperation<OP_TAN>::der(X, Y, F, D); break; \
1405  case OP_ASIN: BinaryOperation<OP_ASIN>::der(X, Y, F, D); break; \
1406  case OP_ACOS: BinaryOperation<OP_ACOS>::der(X, Y, F, D); break; \
1407  case OP_ATAN: BinaryOperation<OP_ATAN>::der(X, Y, F, D); break; \
1408  case OP_LT: BinaryOperation<OP_LT>::der(X, Y, F, D); break; \
1409  case OP_LE: BinaryOperation<OP_LE>::der(X, Y, F, D); break; \
1410  case OP_EQ: BinaryOperation<OP_EQ>::der(X, Y, F, D); break; \
1411  case OP_NE: BinaryOperation<OP_NE>::der(X, Y, F, D); break; \
1412  case OP_NOT: BinaryOperation<OP_NOT>::der(X, Y, F, D); break; \
1413  case OP_AND: BinaryOperation<OP_AND>::der(X, Y, F, D); break; \
1414  case OP_OR: BinaryOperation<OP_OR>::der(X, Y, F, D); break; \
1415  case OP_IF_ELSE_ZERO: BinaryOperation<OP_IF_ELSE_ZERO>::der(X, Y, F, D); break; \
1416  case OP_FLOOR: BinaryOperation<OP_FLOOR>::der(X, Y, F, D); break; \
1417  case OP_CEIL: BinaryOperation<OP_CEIL>::der(X, Y, F, D); break; \
1418  case OP_FMOD: BinaryOperation<OP_FMOD>::der(X, Y, F, D); break; \
1419  case OP_REMAINDER: BinaryOperation<OP_REMAINDER>::der(X, Y, F, D); break; \
1420  case OP_FABS: BinaryOperation<OP_FABS>::der(X, Y, F, D); break; \
1421  case OP_SIGN: BinaryOperation<OP_SIGN>::der(X, Y, F, D); break; \
1422  case OP_COPYSIGN: BinaryOperation<OP_COPYSIGN>::der(X, Y, F, D); break; \
1423  case OP_ERF: BinaryOperation<OP_ERF>::der(X, Y, F, D); break; \
1424  case OP_FMIN: BinaryOperation<OP_FMIN>::der(X, Y, F, D); break; \
1425  case OP_FMAX: BinaryOperation<OP_FMAX>::der(X, Y, F, D); break; \
1426  case OP_INV: BinaryOperation<OP_INV>::der(X, Y, F, D); break; \
1427  case OP_SINH: BinaryOperation<OP_SINH>::der(X, Y, F, D); break; \
1428  case OP_COSH: BinaryOperation<OP_COSH>::der(X, Y, F, D); break; \
1429  case OP_TANH: BinaryOperation<OP_TANH>::der(X, Y, F, D); break; \
1430  case OP_ASINH: BinaryOperation<OP_ASINH>::der(X, Y, F, D); break; \
1431  case OP_ACOSH: BinaryOperation<OP_ACOSH>::der(X, Y, F, D); break; \
1432  case OP_ATANH: BinaryOperation<OP_ATANH>::der(X, Y, F, D); break; \
1433  case OP_ATAN2: BinaryOperation<OP_ATAN2>::der(X, Y, F, D); break; \
1434  case OP_ERFINV: BinaryOperation<OP_ERFINV>::der(X, Y, F, D); break; \
1435  case OP_LIFT: BinaryOperation<OP_LIFT>::der(X, Y, F, D); break; \
1436  case OP_PRINTME: BinaryOperation<OP_PRINTME>::der(X, Y, F, D); break; \
1437  case OP_LOG1P: BinaryOperation<OP_LOG1P>::der(X, Y, F, D); break; \
1438  case OP_EXPM1: BinaryOperation<OP_EXPM1>::der(X, Y, F, D); break; \
1439  case OP_HYPOT: BinaryOperation<OP_HYPOT>::der(X, Y, F, D); break;
1440  switch (op) {
1441  CASADI_MATH_DER_BUILTIN(x, y, f, d)
1442  }
1443  }
1444 
1445 
1446  template<typename T>
1447  inline void casadi_math<T>::derF(unsigned char op, const T& x, const T& y, T& f, T* d) {
1448  // NOTE: We define the implementation in a preprocessor macro to be able to force inlining,
1449  // and to allow extensions in the VM
1450 #define CASADI_MATH_DERF_BUILTIN(X, Y, F, D) \
1451 case OP_ASSIGN: DerBinaryOperation<OP_ASSIGN>::derf(X, Y, F, D); break; \
1452 case OP_ADD: DerBinaryOperation<OP_ADD>::derf(X, Y, F, D); break; \
1453 case OP_SUB: DerBinaryOperation<OP_SUB>::derf(X, Y, F, D); break; \
1454 case OP_MUL: DerBinaryOperation<OP_MUL>::derf(X, Y, F, D); break; \
1455 case OP_DIV: DerBinaryOperation<OP_DIV>::derf(X, Y, F, D); break; \
1456 case OP_NEG: DerBinaryOperation<OP_NEG>::derf(X, Y, F, D); break; \
1457 case OP_EXP: DerBinaryOperation<OP_EXP>::derf(X, Y, F, D); break; \
1458 case OP_LOG: DerBinaryOperation<OP_LOG>::derf(X, Y, F, D); break; \
1459 case OP_POW: DerBinaryOperation<OP_POW>::derf(X, Y, F, D); break; \
1460 case OP_CONSTPOW: DerBinaryOperation<OP_CONSTPOW>::derf(X, Y, F, D); break; \
1461 case OP_SQRT: DerBinaryOperation<OP_SQRT>::derf(X, Y, F, D); break; \
1462 case OP_SQ: DerBinaryOperation<OP_SQ>::derf(X, Y, F, D); break; \
1463 case OP_TWICE: DerBinaryOperation<OP_TWICE>::derf(X, Y, F, D); break; \
1464 case OP_SIN: DerBinaryOperation<OP_SIN>::derf(X, Y, F, D); break; \
1465 case OP_COS: DerBinaryOperation<OP_COS>::derf(X, Y, F, D); break; \
1466 case OP_TAN: DerBinaryOperation<OP_TAN>::derf(X, Y, F, D); break; \
1467 case OP_ASIN: DerBinaryOperation<OP_ASIN>::derf(X, Y, F, D); break; \
1468 case OP_ACOS: DerBinaryOperation<OP_ACOS>::derf(X, Y, F, D); break; \
1469 case OP_ATAN: DerBinaryOperation<OP_ATAN>::derf(X, Y, F, D); break; \
1470 case OP_LT: DerBinaryOperation<OP_LT>::derf(X, Y, F, D); break; \
1471 case OP_LE: DerBinaryOperation<OP_LE>::derf(X, Y, F, D); break; \
1472 case OP_EQ: DerBinaryOperation<OP_EQ>::derf(X, Y, F, D); break; \
1473 case OP_NE: DerBinaryOperation<OP_NE>::derf(X, Y, F, D); break; \
1474 case OP_NOT: DerBinaryOperation<OP_NOT>::derf(X, Y, F, D); break; \
1475 case OP_AND: DerBinaryOperation<OP_AND>::derf(X, Y, F, D); break; \
1476 case OP_OR: DerBinaryOperation<OP_OR>::derf(X, Y, F, D); break; \
1477 case OP_IF_ELSE_ZERO: DerBinaryOperation<OP_IF_ELSE_ZERO>::derf(X, Y, F, D); break; \
1478 case OP_FLOOR: DerBinaryOperation<OP_FLOOR>::derf(X, Y, F, D); break; \
1479 case OP_CEIL: DerBinaryOperation<OP_CEIL>::derf(X, Y, F, D); break; \
1480 case OP_FMOD: DerBinaryOperation<OP_FMOD>::derf(X, Y, F, D); break; \
1481 case OP_REMAINDER: DerBinaryOperation<OP_REMAINDER>::derf(X, Y, F, D); break; \
1482 case OP_FABS: DerBinaryOperation<OP_FABS>::derf(X, Y, F, D); break; \
1483 case OP_SIGN: DerBinaryOperation<OP_SIGN>::derf(X, Y, F, D); break; \
1484 case OP_COPYSIGN: DerBinaryOperation<OP_COPYSIGN>::derf(X, Y, F, D); break; \
1485 case OP_ERF: DerBinaryOperation<OP_ERF>::derf(X, Y, F, D); break; \
1486 case OP_FMIN: DerBinaryOperation<OP_FMIN>::derf(X, Y, F, D); break; \
1487 case OP_FMAX: DerBinaryOperation<OP_FMAX>::derf(X, Y, F, D); break; \
1488 case OP_INV: DerBinaryOperation<OP_INV>::derf(X, Y, F, D); break; \
1489 case OP_SINH: DerBinaryOperation<OP_SINH>::derf(X, Y, F, D); break; \
1490 case OP_COSH: DerBinaryOperation<OP_COSH>::derf(X, Y, F, D); break; \
1491 case OP_TANH: DerBinaryOperation<OP_TANH>::derf(X, Y, F, D); break; \
1492 case OP_ASINH: DerBinaryOperation<OP_ASINH>::derf(X, Y, F, D); break; \
1493 case OP_ACOSH: DerBinaryOperation<OP_ACOSH>::derf(X, Y, F, D); break; \
1494 case OP_ATANH: DerBinaryOperation<OP_ATANH>::derf(X, Y, F, D); break; \
1495 case OP_ATAN2: DerBinaryOperation<OP_ATAN2>::derf(X, Y, F, D); break; \
1496 case OP_ERFINV: DerBinaryOperation<OP_ERFINV>::derf(X, Y, F, D); break; \
1497 case OP_LIFT: DerBinaryOperation<OP_LIFT>::derf(X, Y, F, D); break; \
1498 case OP_PRINTME: DerBinaryOperation<OP_PRINTME>::derf(X, Y, F, D); break; \
1499 case OP_LOG1P: DerBinaryOperation<OP_LOG1P>::derf(X, Y, F, D); break; \
1500 case OP_EXPM1: DerBinaryOperation<OP_EXPM1>::derf(X, Y, F, D); break; \
1501 case OP_HYPOT: DerBinaryOperation<OP_HYPOT>::derf(X, Y, F, D); break;
1502  switch (op) {
1503  CASADI_MATH_DERF_BUILTIN(x, y, f, d)
1504  }
1505  }
1506 
1507  #define CASADI_MATH_BINARY_BUILTIN \
1508  case OP_ADD: \
1509  case OP_SUB: \
1510  case OP_MUL: \
1511  case OP_DIV: \
1512  case OP_POW: \
1513  case OP_CONSTPOW: \
1514  case OP_LT: \
1515  case OP_LE: \
1516  case OP_EQ: \
1517  case OP_NE: \
1518  case OP_AND: \
1519  case OP_OR: \
1520  case OP_COPYSIGN: \
1521  case OP_FMOD: \
1522  case OP_REMAINDER: \
1523  case OP_FMIN: \
1524  case OP_FMAX: \
1525  case OP_ATAN2: \
1526  case OP_PRINTME: \
1527  case OP_LIFT: \
1528  case OP_HYPOT:
1529 
1530  #define CASADI_MATH_UNARY_BUILTIN \
1531  case OP_ASSIGN: \
1532  case OP_NEG: \
1533  case OP_EXP: \
1534  case OP_LOG: \
1535  case OP_SQRT: \
1536  case OP_SQ: \
1537  case OP_TWICE: \
1538  case OP_SIN: \
1539  case OP_COS: \
1540  case OP_TAN: \
1541  case OP_ASIN: \
1542  case OP_ACOS: \
1543  case OP_ATAN: \
1544  case OP_FLOOR: \
1545  case OP_CEIL: \
1546  case OP_NOT: \
1547  case OP_ERF: \
1548  case OP_FABS: \
1549  case OP_SIGN: \
1550  case OP_INV: \
1551  case OP_SINH: \
1552  case OP_COSH: \
1553  case OP_TANH: \
1554  case OP_ASINH: \
1555  case OP_ACOSH: \
1556  case OP_ATANH: \
1557  case OP_ERFINV: \
1558  case OP_LOG1P: \
1559  case OP_EXPM1:
1560 
1561  template<typename T>
1562  inline void casadi_math<T>::fun_linear(unsigned char op, const T* x, const T* y, T* f) {
1563  if (op==OP_ADD || op==OP_SUB) {
1564  for (int i=0;i<3;++i) {
1565  f[i] = T::binary(op, x[i], y[i]);
1566  }
1567  } else if (op==OP_TWICE || op==OP_NEG) {
1568  for (int i=0;i<3;++i) {
1569  f[i] = T::unary(op, x[i]);
1570  }
1571  } else if (op==OP_MUL) {
1572  f[0] += x[0]*y[0];
1573  f[1] += x[0]*y[1];
1574  f[2] += x[0]*y[2];
1575  f[1] += x[1]*y[0];
1576  f[2] += x[1]*y[1];
1577  f[2] += x[1]*y[2];
1578  f[2] += x[2]*y[0];
1579  f[2] += x[2]*y[1];
1580  f[2] += x[2]*y[2];
1581  } else if (op==OP_DIV) {
1582  bool const_argy = y[1].is_zero() && y[2].is_zero();
1583  if (const_argy) {
1584  f[0] = x[0]/y[0];
1585  f[1] = x[1]/y[0];
1586  f[2] = x[2]/y[0];
1587  } else {
1588  f[2] = (x[0]+x[1]+x[2])/(y[0]+y[1]+y[2]);
1589  }
1590  } else if (casadi_math<T>::is_unary(op)) {
1591  bool const_arg = x[1].is_zero() && x[2].is_zero();
1592  if (const_arg) {
1593  f[0] = T::unary(op, x[0]);
1594  } else {
1595  f[2] = T::unary(op, x[0]+x[1]+x[2]);
1596  }
1597 
1598  } else if (casadi_math<T>::is_binary(op)) {
1599  bool const_argx = x[1].is_zero() && x[2].is_zero();
1600  bool const_argy = y[1].is_zero() && y[2].is_zero();
1601  if (const_argx && const_argy) {
1602  f[0] = T::binary(op, x[0], y[0]);
1603  } else {
1604  f[2] = T::binary(op, x[0]+x[1]+x[2], y[0]+y[1]+y[2]);
1605  }
1606  } else {
1607  casadi_error("Not implemented");
1608  }
1609  }
1610 
1611  template<typename T>
1612  bool casadi_math<T>::is_binary(unsigned char op) {
1613  switch (op) {
1614  CASADI_MATH_BINARY_BUILTIN
1615  case OP_IF_ELSE_ZERO:
1616  return true;
1617  default:
1618  return false;
1619  }
1620  }
1621 
1622  template<typename T>
1623  bool casadi_math<T>::is_unary(unsigned char op) {
1624  switch (op) {
1625  CASADI_MATH_UNARY_BUILTIN
1626  return true;
1627  default:
1628  return false;
1629  }
1630  }
1631 
1632  template<typename T>
1633  inline casadi_int casadi_math<T>::ndeps(unsigned char op) {
1634  switch (op) {
1635  case OP_CONST:
1636  case OP_PARAMETER:
1637  case OP_INPUT:
1638  return 0;
1639  CASADI_MATH_BINARY_BUILTIN
1640  case OP_IF_ELSE_ZERO:
1641  return 2;
1642  case OP_CALL:
1643  return -1;
1644  default:
1645  return 1;
1646  }
1647  }
1648 
1649  template<typename T>
1650  inline std::string
1651  casadi_math<T>::print(unsigned char op,
1652  const std::string& x, const std::string& y) {
1653  casadi_assert_dev(ndeps(op)==2);
1654  return pre(op) + x + sep(op) + y + post(op);
1655  }
1656 
1657  template<typename T>
1658  inline std::string
1659  casadi_math<T>::print(unsigned char op, const std::string& x) {
1660  casadi_assert_dev(ndeps(op)==1);
1661  return pre(op) + x + post(op);
1662  }
1663 
1664  template<typename T>
1665  inline std::string casadi_math<T>::name(unsigned char op) {
1666  switch (op) {
1667  case OP_ASSIGN: return "assign";
1668  case OP_ADD: return "add";
1669  case OP_SUB: return "sub";
1670  case OP_MUL: return "mul";
1671  case OP_DIV: return "div";
1672  case OP_NEG: return "neg";
1673  case OP_EXP: return "exp";
1674  case OP_LOG: return "log";
1675  case OP_CONSTPOW:
1676  case OP_POW: return "pow";
1677  case OP_SQRT: return "sqrt";
1678  case OP_SQ: return "sq";
1679  case OP_TWICE: return "twice";
1680  case OP_SIN: return "sin";
1681  case OP_COS: return "cos";
1682  case OP_TAN: return "tan";
1683  case OP_ASIN: return "asin";
1684  case OP_ACOS: return "acos";
1685  case OP_ATAN: return "atan";
1686  case OP_LT: return "lt";
1687  case OP_LE: return "le";
1688  case OP_EQ: return "eq";
1689  case OP_NE: return "ne";
1690  case OP_NOT: return "not";
1691  case OP_AND: return "and";
1692  case OP_OR: return "or";
1693  case OP_FLOOR: return "floor";
1694  case OP_CEIL: return "ceil";
1695  case OP_FMOD: return "fmod";
1696  case OP_REMAINDER: return "remainder";
1697  case OP_FABS: return "fabs";
1698  case OP_SIGN: return "sign";
1699  case OP_COPYSIGN: return "copysign";
1700  case OP_IF_ELSE_ZERO: return "if_else_zero";
1701  case OP_ERF: return "erf";
1702  case OP_FMIN: return "fmin";
1703  case OP_FMAX: return "fmax";
1704  case OP_INV: return "inv";
1705  case OP_SINH: return "sinh";
1706  case OP_COSH: return "cosh";
1707  case OP_TANH: return "tanh";
1708  case OP_ASINH: return "asinh";
1709  case OP_ACOSH: return "acosh";
1710  case OP_ATANH: return "atanh";
1711  case OP_ATAN2: return "atan2";
1712  case OP_CONST: return "const";
1713  case OP_INPUT: return "input";
1714  case OP_OUTPUT: return "output";
1715  case OP_PARAMETER: return "parameter";
1716  case OP_CALL: return "call";
1717  case OP_MTIMES: return "mtimes";
1718  case OP_SOLVE: return "solve";
1719  case OP_TRANSPOSE: return "transpose";
1720  case OP_DETERMINANT: return "determinant";
1721  case OP_INVERSE: return "inverse";
1722  case OP_DOT: return "dot";
1723  case OP_HORZCAT: return "horzcat";
1724  case OP_VERTCAT: return "vertcat";
1725  case OP_DIAGCAT: return "diagcat";
1726  case OP_HORZSPLIT: return "horzsplit";
1727  case OP_VERTSPLIT: return "vertsplit";
1728  case OP_DIAGSPLIT: return "diagsplit";
1729  case OP_RESHAPE: return "reshape";
1730  case OP_SPARSITY_CAST: return "sparsity_cast";
1731  case OP_SUBREF: return "subref";
1732  case OP_SUBASSIGN: return "subassign";
1733  case OP_GETNONZEROS: return "getnonzeros";
1734  case OP_GETNONZEROS_PARAM: return "getnonzeros_param";
1735  case OP_ADDNONZEROS: return "addnonzeros";
1736  case OP_ADDNONZEROS_PARAM: return "addnonzeros_param";
1737  case OP_SETNONZEROS: return "setnonzeros";
1738  case OP_SETNONZEROS_PARAM: return "setnonzeros_param";
1739  case OP_PROJECT: return "project";
1740  case OP_ASSERTION: return "assertion";
1741  case OP_NORM2: return "norm2";
1742  case OP_NORM1: return "norm1";
1743  case OP_NORMINF: return "norminf";
1744  case OP_NORMF: return "normf";
1745  case OP_ERFINV: return "erfinv";
1746  case OP_PRINTME: return "printme";
1747  case OP_LIFT: return "lift";
1748  case OP_EINSTEIN: return "einstein";
1749  case OP_BSPLINE: return "bspline";
1750  case OP_CONVEXIFY: return "convexify";
1751  case OP_LOG1P: return "log1p";
1752  case OP_EXPM1: return "expm1";
1753  case OP_HYPOT: return "hypot";
1754  case OP_LOGSUMEXP: return "logsumexp";
1755  case OP_KRON: return "kron";
1756  case OP_KRON_CONTRACT: return "kron_contract";
1757  }
1758  return "<invalid-op>";
1759  }
1760 
1761  template<typename T>
1762  inline std::string casadi_math<T>::pre(unsigned char op) {
1763  switch (op) {
1764  case OP_ASSIGN: return "";
1765  case OP_ADD: return "(";
1766  case OP_SUB: return "(";
1767  case OP_MUL: return "(";
1768  case OP_DIV: return "(";
1769  case OP_NEG: return "(-";
1770  case OP_TWICE: return "(2.*";
1771  case OP_LT: return "(";
1772  case OP_LE: return "(";
1773  case OP_EQ: return "(";
1774  case OP_NE: return "(";
1775  case OP_NOT: return "(!";
1776  case OP_AND: return "(";
1777  case OP_OR: return "(";
1778  case OP_IF_ELSE_ZERO: return "(";
1779  case OP_INV: return "(1./";
1780  default: return name(op) + "(";
1781  }
1782  }
1783 
1784  template<typename T>
1785  inline std::string casadi_math<T>::sep(unsigned char op) {
1786  switch (op) {
1787  case OP_ADD: return "+";
1788  case OP_SUB: return "-";
1789  case OP_MUL: return "*";
1790  case OP_DIV: return "/";
1791  case OP_LT: return "<";
1792  case OP_LE: return "<=";
1793  case OP_EQ: return "==";
1794  case OP_NE: return "!=";
1795  case OP_AND: return "&&";
1796  case OP_OR: return "||";
1797  case OP_IF_ELSE_ZERO: return "?";
1798  default: return ",";
1799  }
1800  }
1801 
1802  template<typename T>
1803  inline std::string casadi_math<T>::post(unsigned char op) {
1804  switch (op) {
1805  case OP_ASSIGN: return "";
1806  case OP_IF_ELSE_ZERO: return ":0)";
1807  default: return ")";
1808  }
1809  }
1810 
1811 
1812  template<class T, typename SU>
1813  T common_simp_unary(casadi_int op, const T& x, casadi_int depth,
1814  SU&& gen_unary,
1815  bool unique,
1816  bool& hit) {
1817  hit = true;
1818  switch (op) {
1819  case OP_TWICE:
1820  if (x.is_op(OP_MUL) && x.dep(0).is_constant() &&
1821  x.dep(0).is_half())
1822  return x.dep(1); // 2*(0.5*x) = x
1823  else if (x.is_op(OP_MUL) && x.dep(1).is_constant() &&
1824  x.dep(1).is_half())
1825  return x.dep(0); // 2*(x*0.5) = x
1826  break;
1827  case OP_SQ:
1828  if (x.is_op(OP_SQRT))
1829  return x.dep(); // sqrt(x)^2 = x
1830  else if (x.is_op(OP_NEG))
1831  return sq(x.dep()); // (-x)^2 = x^2
1832  else if (x.is_op(OP_FABS))
1833  return sq(x.dep()); // |x|^2 = x^2
1834  break;
1835  case OP_FABS:
1836  if (x.is_nonnegative())
1837  return x;
1838  else if (x.is_op(OP_NEG))
1839  return fabs(x.dep()); // fabs(-x) = fabs(x)
1840  break;
1841  case OP_LOG:
1842  if (x.is_op(OP_EXP))
1843  return x.dep(); // log(exp(x)) = x
1844  break;
1845  case OP_INV:
1846  if (x.is_op(OP_INV))
1847  return x.dep(); // 1/(1/x) = x
1848  break;
1849  case OP_SQRT:
1850  if (x.is_op(OP_SQ))
1851  return fabs(x.dep()); // sqrt(x^2) = x
1852  break;
1853  case OP_COS:
1854  if (x.is_op(OP_NEG))
1855  return cos(x.dep()); // cos(-x) = cos(x)
1856  else if (x.is_op(OP_FABS))
1857  return cos(x.dep()); // cos(|x|) = cos(x)
1858  break;
1859  case OP_NEG:
1860  if (x.is_op(OP_NEG))
1861  return x.dep(); // -(-x) = x
1862  else if (unique && x.is_op(OP_SUB))
1863  return x.dep(1) - x.dep(0); // -(x-y) = y-x
1864  }
1865  hit = false;
1866  return 0;
1867  }
1868 
1869 
1870  template<class T, typename SU, typename SB>
1871  T common_simp_binary(casadi_int op, const T& x, const T& y, casadi_int depth,
1872  SU&& gen_unary,
1873  SB&& gen_binary,
1874  bool unique_x,
1875  bool unique_y,
1876  bool& hit) {
1877  hit = true;
1878  switch (op) {
1879  case OP_ADD:
1880  if (x.is_zero())
1881  return y;
1882  else if (y->is_zero()) // term2 is zero
1883  return x;
1884  else if (y.is_op(OP_NEG)) // x + (-y) -> x - y
1885  return x - (-y);
1886  else if (x.is_op(OP_NEG)) // (-x) + y -> y - x
1887  return y - x.dep();
1888  else if (is_equal(x, y, depth))
1889  return 2*x; // x+x -> 2*x
1890  else if (x.is_op(OP_MUL) && y.is_op(OP_MUL) &&
1891  x.dep(0).is_constant() && static_cast<double>(x.dep(0))==0.5 &&
1892  y.dep(0).is_constant() && static_cast<double>(y.dep(0))==0.5 &&
1893  is_equal(y.dep(1), x.dep(1), depth)) // 0.5x+0.5x = x
1894  return x.dep(1);
1895  else if (x.is_op(OP_DIV) && y.is_op(OP_DIV) &&
1896  x.dep(1).is_constant() && static_cast<double>(x.dep(1))==2 &&
1897  y.dep(1).is_constant() && static_cast<double>(y.dep(1))==2 &&
1898  is_equal(y.dep(0), x.dep(0), depth)) // x/2+x/2 = x
1899  return x.dep(0);
1900  else if (x.is_op(OP_MUL) && y.is_op(OP_MUL) &&
1901  x.dep(0).is_constant() && y.dep(0).is_constant() &&
1902  static_cast<double>(x.dep(0))+static_cast<double>(y.dep(0))==1 &&
1903  is_equal(y.dep(1), x.dep(1), depth)
1904  )
1905  return x.dep(1);
1906  else if (x.is_op(OP_SUB) && is_equal(x.dep(1), y, depth))
1907  return x.dep(0);
1908  else if (y.is_op(OP_SUB) && is_equal(x, y.dep(1), depth))
1909  return y.dep(0);
1910  else if (x.is_op(OP_SQ) && y.is_op(OP_SQ) &&
1911  ((x.dep().is_op(OP_SIN) && y.dep().is_op(OP_COS))
1912  || (x.dep().is_op(OP_COS) && y.dep().is_op(OP_SIN)))
1913  && is_equal(x.dep().dep(), y.dep().dep(), depth))
1914  return 1; // sin^2 + cos^2 -> 1
1915  else if (x.is_doubled() && is_equal(x.dep(0), y, depth))
1916  return 3*x.dep(0);
1917  else if (y.is_doubled() && is_equal(y.dep(0), x, depth))
1918  return 3*y.dep(0);
1919  else if (x.is_op(OP_MUL) && y.is_op(OP_MUL) &&
1920  is_equal(x.dep(0), y.dep(0), depth) && (x.dep(1)+y.dep(1)).is_one())
1921  return x.dep(0);
1922  break;
1923  case OP_SUB:
1924  if (y.is_zero()) // term2 is zero
1925  return x;
1926  if (x.is_zero()) // term1 is zero
1927  return -y;
1928  if (is_equal(x, y, depth)) // the terms are equal
1929  return 0;
1930  else if (y.is_op(OP_NEG)) // x - (-y) -> x + y
1931  return x + y.dep();
1932  else if (x.is_op(OP_ADD) && is_equal(x.dep(1), y, depth))
1933  return x.dep(0);
1934  else if (x.is_op(OP_ADD) && is_equal(x.dep(0), y, depth))
1935  return x.dep(1);
1936  else if (y.is_op(OP_ADD) && is_equal(x, y.dep(1), depth))
1937  return -y.dep(0);
1938  else if (y.is_op(OP_ADD) && is_equal(x, y.dep(0), depth))
1939  return -y.dep(1);
1940  else if (x.is_op(OP_NEG))
1941  return -(x.dep() + y);
1942  else if (x.is_op(OP_MUL) && y.is_op(OP_MUL) &&
1943  x.dep(0).is_constant() && y.dep(0).is_constant() &&
1944  static_cast<double>(x.dep(0))-static_cast<double>(y.dep(0))==1 &&
1945  is_equal(y.dep(1), x.dep(1), depth)
1946  )
1947  return x.dep(1);
1948  else if (x.is_doubled() && is_equal(x.dep(0), y, depth))
1949  return y;
1950  else if (y.is_doubled() && is_equal(y.dep(0), x, depth))
1951  return -x;
1952  break;
1953  case OP_MUL:
1954  if (is_equal(y, x, depth))
1955  return sq(x);
1956  else if (!x.is_constant() && y.is_constant())
1957  return y * x;
1958  else if (x.is_zero() || y->is_zero()) // one of the terms is zero
1959  return 0;
1960  else if (x.is_one()) // term1 is one
1961  return y;
1962  else if (y.is_one()) // term2 is one
1963  return x;
1964  else if (y.is_minus_one())
1965  return -x;
1966  else if (x.is_minus_one())
1967  return -y;
1968  else if (y.is_op(OP_INV))
1969  return x/y.inv();
1970  else if (x.is_op(OP_INV))
1971  return y / x.inv();
1972  else if (x.is_constant() && y.is_op(OP_MUL) && y.dep(0).is_constant() &&
1973  static_cast<double>(x)*static_cast<double>(y.dep(0))==1) // 5*(0.2*x) = x
1974  return y.dep(1);
1975  else if (x.is_constant() && y.is_doubled() &&
1976  static_cast<double>(x)*2==1) // 0.5*(2*x) = x
1977  return y.dep(0);
1978  else if (x.is_constant() && y.is_op(OP_DIV) && y.dep(1).is_constant() &&
1979  static_cast<double>(x)==static_cast<double>(y.dep(1))) // 5*(x/5) = x
1980  return y.dep(0);
1981  else if (x.is_constant() && static_cast<double>(x)==2)
1982  return gen_unary(OP_TWICE, y);
1983  else if (y.is_constant() && static_cast<double>(y)==2)
1984  return gen_unary(OP_TWICE, x);
1985  else if (x.is_op(OP_DIV) && is_equal(x.dep(1), y, depth)) // ((2/x)*x)
1986  return x.dep(0);
1987  else if (y.is_op(OP_DIV) &&
1988  is_equal(y.dep(1), x, depth)) // ((2/x)*x)
1989  return y.dep(0);
1990  else if (x.is_op(OP_NEG))
1991  return -(x.dep() * y);
1992  else if (y.is_op(OP_NEG))
1993  return -(x * y.dep());
1994  else if (unique_x && x.is_op(OP_ADD) && x.dep(1).is_op(OP_DIV) &&
1995  is_equal(x.dep(1).dep(1), y, depth))
1996  return x.dep(0)*y + x.dep(1).dep(0); // (a + b/c)*c -> a*c + b
1997  else if (unique_x && x.is_op(OP_ADD) && x.dep(0).is_op(OP_DIV) &&
1998  is_equal(x.dep(0).dep(1), y, depth))
1999  return x.dep(0).dep(0)+x.dep(1)*y; // (a/c + b)*c -> a*c + b
2000  else if (unique_y && y.is_op(OP_ADD) && y.dep(1).is_op(OP_DIV) &&
2001  is_equal(y.dep(1).dep(1), x, depth))
2002  return x*y.dep(0) + y.dep(1).dep(0); // c*(a + b/c) -> c*a + b
2003  else if (unique_y && y.is_op(OP_ADD) && y.dep(0).is_op(OP_DIV) &&
2004  is_equal(y.dep(0).dep(1), x, depth))
2005  return y.dep(0).dep(0) + x*y.dep(1); // c*(a/c + b) -> a + c*b
2006  else if (unique_x && x.is_op(OP_SUB) && x.dep(1).is_op(OP_DIV) &&
2007  is_equal(x.dep(1).dep(1), y, depth))
2008  return x.dep(0)*y - x.dep(1).dep(0); // (a - b/c)*c -> a*c - b
2009  else if (unique_x && x.is_op(OP_SUB) && x.dep(0).is_op(OP_DIV) &&
2010  is_equal(x.dep(0).dep(1), y, depth))
2011  return x.dep(0).dep(0)-x.dep(1)*y; // (a/c - b)*c -> a*c - b
2012  else if (unique_y && y.is_op(OP_SUB) && y.dep(1).is_op(OP_DIV) &&
2013  is_equal(y.dep(1).dep(1), x, depth))
2014  return x*y.dep(0) - y.dep(1).dep(0); // c*(a - b/c) -> c*a - b
2015  else if (unique_y && y.is_op(OP_SUB) && y.dep(0).is_op(OP_DIV) &&
2016  is_equal(y.dep(0).dep(1), x, depth))
2017  return y.dep(0).dep(0) - x*y.dep(1); // c*(a/c - b) -> a - c*b
2018  else if (x.is_constant() && y.is_op(OP_MUL) && y.dep(0).is_constant())
2019  return x*y.dep(0)*y.dep(1); // a*(b*x) -> (a*b)*x
2020  else if (y.is_constant() && x.is_op(OP_MUL) && x.dep(0).is_constant())
2021  return y*x.dep(0)*x.dep(1); // (a*x)*b -> (a*b)*x
2022  break;
2023  case OP_DIV:
2024  if (y.is_zero()) // term2 is zero
2025  return std::numeric_limits<double>::infinity();
2026  else if (x.is_zero()) // term1 is zero
2027  return 0;
2028  else if (y.is_one()) // term2 is one
2029  return x;
2030  else if (y.is_minus_one())
2031  return -x;
2032  else if (y.is_half())
2033  return 2*x;
2034  else if (is_equal(x, y, depth)) // terms are equal
2035  return 1;
2036  else if (x.is_doubled() && y.is_constant() && static_cast<double>(y)==2)
2037  return x.dep(0);
2038  else if (x.is_doubled() && is_equal(y, x.dep(0), depth))
2039  return 2;
2040  else if (x.is_op(OP_MUL) && is_equal(y, x.dep(0), depth))
2041  return x.dep(1);
2042  else if (x.is_op(OP_MUL) && is_equal(y, x.dep(1), depth))
2043  return x.dep(0);
2044  else if (x.is_one())
2045  return y.inv();
2046  else if (y.is_op(OP_INV))
2047  return x*y.inv();
2048  else if (x.is_doubled() && y.is_doubled())
2049  return x.dep(0) / y->dep(0);
2050  else if (y.is_constant() && x.is_op(OP_DIV) && x.dep(1).is_constant() &&
2051  static_cast<double>(y)*static_cast<double>(x.dep(1))==1) // (x/5)/0.2
2052  return x.dep(0);
2053  else if (y.is_op(OP_MUL) &&
2054  is_equal(y.dep(1), x, depth)) // x/(2*x) = 1/2
2055  return 1/y.dep(0); //BinarySX::create(OP_DIV, 1, y.dep(0));
2056  else if (y.is_doubled() &&
2057  is_equal(y.dep(0), x, depth)) // x/(2*x) = 1/2
2058  return 0.5;
2059  else if (x.is_op(OP_NEG) &&
2060  is_equal(x.dep(0), y, depth)) // (-x)/x = -1
2061  return -1;
2062  else if (y.is_op(OP_NEG) &&
2063  is_equal(y.dep(0), x, depth)) // x/(-x) = 1
2064  return -1;
2065  else if (y.is_op(OP_NEG) && x.is_op(OP_NEG) &&
2066  is_equal(x.dep(0), y.dep(0), depth)) // (-x)/(-x) = 1
2067  return 1;
2068  else if (x.is_op(OP_DIV) && is_equal(y, x.dep(0), depth))
2069  return x.dep(1).inv();
2070  else if (x.is_op(OP_NEG))
2071  return -(x.dep() / y);
2072  else if (y.is_op(OP_NEG))
2073  return -(x / y.dep());
2074  break;
2075  case OP_POW:
2076  if (y.is_constant()) {
2077  if (y.is_integer()) {
2078  casadi_int nn = y->to_int();
2079  if (nn == 0) {
2080  return 1;
2081  } else if (nn>100 || nn<-100) { // maximum depth
2082  return gen_binary(OP_CONSTPOW, x, nn);
2083  } else if (nn<0) { // negative power
2084  return 1/pow(x, -nn);
2085  } else if (nn%2 == 1) { // odd power
2086  return x*pow(x, nn-1);
2087  } else { // even power
2088  T rt = pow(x, static_cast<casadi_int>(nn/2));
2089  return rt*rt;
2090  }
2091  } else if (y->to_double()==0.5) {
2092  return sqrt(x);
2093  } else {
2094  return gen_binary(OP_CONSTPOW, x, y);
2095  }
2096  }
2097  break;
2098  case OP_LE:
2099  if ((y-x).is_nonnegative())
2100  return 1;
2101  break;
2102  case OP_FMIN:
2103  if (x.is_inf()) return y;
2104  if (y.is_inf()) return x;
2105  if (x.is_minus_inf() || y.is_minus_inf()) return -std::numeric_limits<double>::infinity();
2106  if (is_equal(x, y, depth)) return x;
2107  break;
2108  case OP_FMAX:
2109  if (x.is_minus_inf()) return y;
2110  if (y.is_minus_inf()) return x;
2111  if (x.is_inf() || y.is_inf()) return std::numeric_limits<double>::infinity();
2112  if (is_equal(x, y, depth)) return x;
2113  break;
2114  case OP_LT:
2115  if (((x)-y).is_nonnegative())
2116  return 0;
2117  break;
2118  case OP_EQ:
2119  if (is_equal(x, y, depth))
2120  return 1;
2121  break;
2122  case OP_NE:
2123  if (is_equal(x, y, depth))
2124  return 0;
2125  break;
2126  case OP_IF_ELSE_ZERO:
2127  if (y.is_zero()) {
2128  return y;
2129  } else if (x.is_constant()) {
2130  if (static_cast<double>(x)!=0) {
2131  return y;
2132  } else {
2133  return 0;
2134  }
2135  }
2136  }
2137  hit = false;
2138  return 0;
2139  }
2140 
2141 #endif // SWIG
2142 
2143 } // namespace casadi
2144 
2146 
2147 #endif // CASADI_CALCULUS_HPP
casadi
The casadi namespace.
Definition: archiver.hpp:32
casadi::uout
CASADI_EXPORT std::ostream & uout()