Scippy

SCIP

Solving Constraint Integer Programs

expr.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2019 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file nlpi/expr.c
17  * @brief methods for expressions, expression trees, expression graphs, and related
18  * @author Stefan Vigerske
19  * @author Thorsten Gellermann
20  * @author Ingmar Vierhaus (exprparse)
21  */
22 
23 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
24 
25 #include <stdarg.h>
26 #include <string.h>
27 #include <math.h>
28 #include <ctype.h>
29 
30 #include "nlpi/pub_expr.h"
31 #include "nlpi/struct_expr.h"
32 #include "nlpi/exprinterpret.h"
33 
34 #include "scip/intervalarith.h"
35 #include "scip/pub_misc.h"
36 #include "scip/misc.h"
37 #include "scip/pub_message.h"
38 
39 
40 #define SCIP_EXPRESSION_MAXCHILDEST 16 /**< estimate on maximal number of children */
41 
42 /** sign of a value (-1 or +1)
43  *
44  * 0.0 has sign +1
45  */
46 #define SIGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
47 
48 /** ensures that a block memory array has at least a given size
49  *
50  * if cursize is 0, then *array1 can be NULL
51  */
52 #define ensureBlockMemoryArraySize(blkmem, array1, cursize, minsize) \
53  do { \
54  int __newsize; \
55  assert((blkmem) != NULL); \
56  if( *(cursize) >= (minsize) ) \
57  break; \
58  __newsize = calcGrowSize(minsize); \
59  assert(__newsize >= (minsize)); \
60  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, array1, *(cursize), __newsize) ); \
61  *(cursize) = __newsize; \
62  } while( FALSE )
63 
64 #ifdef SCIP_DISABLED_CODE /* this macro is currently not used, which offends lint, so disable it */
65 /** ensures that two block memory arrays have at least a given size
66  *
67  * if cursize is 0, then arrays can be NULL
68  */
69 #define ensureBlockMemoryArraySize2(blkmem, array1, array2, cursize, minsize) \
70  do { \
71  int __newsize; \
72  assert((blkmem) != NULL); \
73  if( *(cursize) >= (minsize) ) \
74  break; \
75  __newsize = calcGrowSize(minsize); \
76  assert(__newsize >= (minsize)); \
77  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, array1, *(cursize), __newsize) ); \
78  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, array2, *(cursize), __newsize) ); \
79  *(cursize) = __newsize; \
80  } while( FALSE )
81 #endif
82 
83 /** ensures that three block memory arrays have at least a given size
84  *
85  * if cursize is 0, then arrays can be NULL
86  */
87 #define ensureBlockMemoryArraySize3(blkmem, array1, array2, array3, cursize, minsize) \
88  do { \
89  int __newsize; \
90  assert((blkmem) != NULL); \
91  if( *(cursize) >= (minsize) ) \
92  break; \
93  __newsize = calcGrowSize(minsize); \
94  assert(__newsize >= (minsize)); \
95  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, array1, *(cursize), __newsize) ); \
96  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, array2, *(cursize), __newsize) ); \
97  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, array3, *(cursize), __newsize) ); \
98  *(cursize) = __newsize; \
99  } while( FALSE )
100 
101 /**@name Miscellaneous private methods */
102 /**@{ */
103 
104 /** calculate memory size for dynamically allocated arrays (copied from scip/set.c) */
105 static
107  int num /**< minimum number of entries to store */
108  )
109 {
110  int size;
111 
112  /* calculate the size with this loop, such that the resulting numbers are always the same (-> block memory) */
113  size = 4;
114  while( size < num )
115  size = (int)(1.2 * size + 4);
116 
117  return size;
118 }
119 
120 /** expression graph nodes comparison to use in sorting methods
121  *
122  * The nodes need to have been added to the expression graph (depth,pos >= 0).
123  * The better node is the one with the lower depth and lower position, if depth is equal.
124  */
125 static
126 SCIP_DECL_SORTPTRCOMP(exprgraphnodecomp)
127 {
128  SCIP_EXPRGRAPHNODE* node1 = (SCIP_EXPRGRAPHNODE*)elem1;
129  SCIP_EXPRGRAPHNODE* node2 = (SCIP_EXPRGRAPHNODE*)elem2;
130 
131  assert(node1 != NULL);
132  assert(node2 != NULL);
133  assert(node1->depth >= 0);
134  assert(node1->pos >= 0);
135  assert(node2->depth >= 0);
136  assert(node2->pos >= 0);
137 
138  if( node1->depth != node2->depth )
139  return node1->depth - node2->depth;
140 
141  /* there should be no two nodes on the same position */
142  assert((node1->pos != node2->pos) || (node1 == node2));
143 
144  return node1->pos - node2->pos;
145 }
146 
147 /** checks if a given new lower bound is tighter (w.r.t. given bound strengthening epsilon) than the old one (copied from scip/set.c) */
148 static
150  SCIP_Real minstrength, /**< minimal relative improvement required to be a better bound */
151  SCIP_Real newlb, /**< new lower bound */
152  SCIP_Real oldlb, /**< old lower bound */
153  SCIP_Real oldub /**< old upper bound */
154  )
155 {
156  SCIP_Real eps;
157 
158  /* nothing can be tighter than an empty interval */
159  if( oldlb > oldub )
160  return FALSE;
161 
162  eps = REALABS(oldlb);
163  eps = MIN(oldub - oldlb, eps);
164  return EPSGT(newlb, oldlb, minstrength * MAX(eps, 1e-3));
165 }
166 
167 /** checks if a given new upper bound is tighter (w.r.t. given bound strengthening epsilon) than the old one (copied from scip/set.c) */
168 static
170  SCIP_Real minstrength, /**< minimal relative improvement required to be a better bound */
171  SCIP_Real newub, /**< new upper bound */
172  SCIP_Real oldlb, /**< old lower bound */
173  SCIP_Real oldub /**< old upper bound */
174  )
175 {
176  SCIP_Real eps;
177 
178  /* nothing can be tighter than an empty interval */
179  if( oldlb > oldub )
180  return FALSE;
181 
182  eps = REALABS(oldub);
183  eps = MIN(oldub - oldlb, eps);
184  return EPSLT(newub, oldub, minstrength * MAX(eps, 1e-3));
185 }
186 
187 /**@} */
188 
189 /**@name Expression curvature methods */
190 /**@{ */
191 
192 /** curvature names as strings */
193 static
194 const char* curvnames[4] =
195  {
196  "unknown",
197  "convex",
198  "concave",
199  "linear"
200  };
201 
202 #undef SCIPexprcurvAdd
203 
204 /** gives curvature for a sum of two functions with given curvature */
206  SCIP_EXPRCURV curv1, /**< curvature of first summand */
207  SCIP_EXPRCURV curv2 /**< curvature of second summand */
208  )
209 {
210  return (SCIP_EXPRCURV) (curv1 & curv2);
211 }
212 
213 /** gives the curvature for the negation of a function with given curvature */
215  SCIP_EXPRCURV curvature /**< curvature of function */
216  )
217 {
218  switch( curvature )
219  {
221  return SCIP_EXPRCURV_CONVEX;
222 
224  return SCIP_EXPRCURV_CONCAVE;
225 
228  /* can return curvature, do this below */
229  break;
230 
231  default:
232  SCIPerrorMessage("unknown curvature status.\n");
233  SCIPABORT();
234  }
235 
236  return curvature;
237 }
238 
239 /** gives curvature for a functions with given curvature multiplied by a constant factor */
241  SCIP_Real factor, /**< constant factor */
242  SCIP_EXPRCURV curvature /**< curvature of other factor */
243  )
244 {
245  if( factor == 0.0 )
246  return SCIP_EXPRCURV_LINEAR;
247  if( factor > 0.0 )
248  return curvature;
249  return SCIPexprcurvNegate(curvature);
250 }
251 
252 /** gives curvature for base^exponent for given bounds and curvature of base-function and constant exponent */
254  SCIP_INTERVAL basebounds, /**< bounds on base function */
255  SCIP_EXPRCURV basecurv, /**< curvature of base function */
256  SCIP_Real exponent /**< exponent */
257  )
258 {
259  SCIP_Bool expisint;
260 
261  assert(basebounds.inf <= basebounds.sup);
262 
263  if( exponent == 0.0 )
264  return SCIP_EXPRCURV_LINEAR;
265 
266  if( exponent == 1.0 )
267  return basecurv;
268 
269  expisint = EPSISINT(exponent, 0.0); /*lint !e835*/
270 
271  /* if exponent is fractional, then power is not defined for a negative base
272  * thus, consider only positive part of basebounds
273  */
274  if( !expisint && basebounds.inf < 0.0 )
275  {
276  basebounds.inf = 0.0;
277  if( basebounds.sup < 0.0 )
278  return SCIP_EXPRCURV_LINEAR;
279  }
280 
281  /* if basebounds contains 0.0, consider negative and positive interval separately, if possible */
282  if( basebounds.inf < 0.0 && basebounds.sup > 0.0 )
283  {
284  SCIP_INTERVAL leftbounds;
285  SCIP_INTERVAL rightbounds;
286 
287  /* something like x^(-2) may look convex on each side of zero, but is not convex on the whole interval due to the singularity at 0.0 */
288  if( exponent < 0.0 )
289  return SCIP_EXPRCURV_UNKNOWN;
290 
291  SCIPintervalSetBounds(&leftbounds, basebounds.inf, 0.0);
292  SCIPintervalSetBounds(&rightbounds, 0.0, basebounds.sup);
293 
294  return (SCIP_EXPRCURV) (SCIPexprcurvPower(leftbounds, basecurv, exponent) & SCIPexprcurvPower(rightbounds, basecurv, exponent));
295  }
296  assert(basebounds.inf >= 0.0 || basebounds.sup <= 0.0);
297 
298  /* (base^exponent)'' = exponent * ( (exponent-1) base^(exponent-2) (base')^2 + base^(exponent-1) base'' )
299  *
300  * if base'' is positive, i.e., base is convex, then
301  * - for base > 0.0 and exponent > 1.0, the second deriv. is positive -> convex
302  * - for base < 0.0 and exponent > 1.0, we can't say (first and second summand opposite signs)
303  * - for base > 0.0 and 0.0 < exponent < 1.0, we can't say (first sommand negative, second summand positive)
304  * - for base > 0.0 and exponent < 0.0, we can't say (first and second summand opposite signs)
305  * - for base < 0.0 and exponent < 0.0 and even, the second deriv. is positive -> convex
306  * - for base < 0.0 and exponent < 0.0 and odd, the second deriv. is negative -> concave
307  *
308  * if base'' is negative, i.e., base is concave, then
309  * - for base > 0.0 and exponent > 1.0, we can't say (first summand positive, second summand negative)
310  * - for base < 0.0 and exponent > 1.0 and even, the second deriv. is positive -> convex
311  * - for base < 0.0 and exponent > 1.0 and odd, the second deriv. is negative -> concave
312  * - for base > 0.0 and 0.0 < exponent < 1.0, the second deriv. is negative -> concave
313  * - for base > 0.0 and exponent < 0.0, the second deriv. is positive -> convex
314  * - for base < 0.0 and exponent < 0.0, we can't say (first and second summand opposite signs)
315  *
316  * if base'' is zero, i.e., base is linear, then
317  * (base^exponent)'' = exponent * (exponent-1) base^(exponent-2) (base')^2
318  * - just multiply signs
319  */
320 
321  if( basecurv == SCIP_EXPRCURV_LINEAR )
322  {
323  SCIP_Real sign;
324 
325  /* base^(exponent-2) is negative, if base < 0.0 and exponent is odd */
326  sign = exponent * (exponent - 1.0);
327  assert(basebounds.inf >= 0.0 || expisint);
328  if( basebounds.inf < 0.0 && ((int)exponent)%2 != 0 )
329  sign *= -1.0;
330  assert(sign != 0.0);
331 
332  return sign > 0.0 ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
333  }
334 
335  if( basecurv == SCIP_EXPRCURV_CONVEX )
336  {
337  if( basebounds.sup <= 0.0 && exponent < 0.0 && expisint )
338  return ((int)exponent)%2 == 0 ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
339  if( basebounds.inf >= 0.0 && exponent > 1.0 )
340  return SCIP_EXPRCURV_CONVEX ;
341  return SCIP_EXPRCURV_UNKNOWN;
342  }
343 
344  if( basecurv == SCIP_EXPRCURV_CONCAVE )
345  {
346  if( basebounds.sup <= 0.0 && exponent > 1.0 && expisint )
347  return ((int)exponent)%2 == 0 ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
348  if( basebounds.inf >= 0.0 && exponent < 1.0 )
349  return exponent < 0.0 ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
350  return SCIP_EXPRCURV_UNKNOWN;
351  }
352 
353  return SCIP_EXPRCURV_UNKNOWN;
354 }
355 
356 /** gives curvature for a monomial with given curvatures and bounds for each factor
357  *
358  * See Maranas and Floudas, Finding All Solutions of Nonlinearly Constrained Systems of Equations, JOGO 7, 1995
359  * for the categorization in the case that all factors are linear.
360  */
362  int nfactors, /**< number of factors in monomial */
363  SCIP_Real* exponents, /**< exponents in monomial, or NULL if all 1.0 */
364  int* factoridxs, /**< indices of factors (but not exponents), or NULL if identity mapping */
365  SCIP_EXPRCURV* factorcurv, /**< curvature of each factor */
366  SCIP_INTERVAL* factorbounds /**< bounds of each factor */
367  )
368 {
369  SCIP_Real mult;
370  SCIP_Real e;
371  SCIP_EXPRCURV curv;
372  SCIP_EXPRCURV fcurv;
373  int nnegative;
374  int npositive;
375  SCIP_Real sum;
376  SCIP_Bool expcurvpos;
377  SCIP_Bool expcurvneg;
378  int j;
379  int f;
380 
381  assert(nfactors >= 0);
382  assert(factorcurv != NULL || nfactors == 0);
383  assert(factorbounds != NULL || nfactors == 0);
384 
385  if( nfactors == 0 )
386  return SCIP_EXPRCURV_LINEAR;
387 
388  if( nfactors == 1 )
389  {
390  f = factoridxs != NULL ? factoridxs[0] : 0;
391  e = exponents != NULL ? exponents[0] : 1.0;
392  /* SCIPdebugMessage("monomial [%g,%g]^%g is %s\n",
393  factorbounds[f].inf, factorbounds[f].sup, e,
394  SCIPexprcurvGetName(SCIPexprcurvPower(factorbounds[f], factorcurv[f], e))); */
395  return SCIPexprcurvPower(factorbounds[f], factorcurv[f], e); /*lint !e613*/
396  }
397 
398  mult = 1.0;
399 
400  nnegative = 0; /* number of negative exponents */
401  npositive = 0; /* number of positive exponents */
402  sum = 0.0; /* sum of exponents */
403  expcurvpos = TRUE; /* whether exp_j * f_j''(x) >= 0 for all factors (assuming f_j >= 0) */
404  expcurvneg = TRUE; /* whether exp_j * f_j''(x) <= 0 for all factors (assuming f_j >= 0) */
405 
406  for( j = 0; j < nfactors; ++j )
407  {
408  f = factoridxs != NULL ? factoridxs[j] : j;
409  if( factorcurv[f] == SCIP_EXPRCURV_UNKNOWN ) /*lint !e613*/
410  return SCIP_EXPRCURV_UNKNOWN;
411  if( factorbounds[f].inf < 0.0 && factorbounds[f].sup > 0.0 ) /*lint !e613*/
412  return SCIP_EXPRCURV_UNKNOWN;
413 
414  e = exponents != NULL ? exponents[j] : 1.0;
415  if( e < 0.0 )
416  ++nnegative;
417  else
418  ++npositive;
419  sum += e;
420 
421  if( factorbounds[f].inf < 0.0 ) /*lint !e613*/
422  {
423  /* if argument is negative, then exponent should be integer */
424  assert(EPSISINT(e, 0.0)); /*lint !e835*/
425 
426  /* flip j'th argument: (f_j)^(exp_j) = (-1)^(exp_j) (-f_j)^(exp_j) */
427 
428  /* -f_j has negated curvature of f_j */
429  fcurv = SCIPexprcurvNegate(factorcurv[f]); /*lint !e613*/
430 
431  /* negate monomial, if exponent is odd, i.e., (-1)^(exp_j) = -1 */
432  if( (int)e % 2 != 0 )
433  mult *= -1.0;
434  }
435  else
436  {
437  fcurv = factorcurv[f]; /*lint !e613*/
438  }
439 
440  /* check if exp_j * fcurv is convex (>= 0) and/or concave */
441  fcurv = SCIPexprcurvMultiply(e, fcurv);
442  if( !(fcurv & SCIP_EXPRCURV_CONVEX) )
443  expcurvpos = FALSE;
444  if( !(fcurv & SCIP_EXPRCURV_CONCAVE) )
445  expcurvneg = FALSE;
446  }
447 
448  /* if all factors are linear, then a product f_j^exp_j with f_j >= 0 is convex if
449  * - all exponents are negative, or
450  * - all except one exponent j* are negative and exp_j* >= 1 - sum_{j!=j*}exp_j, but the latter is equivalent to sum_j exp_j >= 1
451  * further, the product is concave if
452  * - all exponents are positive and the sum of exponents is <= 1.0
453  *
454  * if factors are nonlinear, then we require additionally, that for convexity
455  * - each factor is convex if exp_j >= 0, or concave if exp_j <= 0, i.e., exp_j*f_j'' >= 0
456  * and for concavity, we require that
457  * - all factors are concave, i.e., exp_j*f_j'' <= 0
458  */
459 
460  if( nnegative == nfactors && expcurvpos )
461  curv = SCIP_EXPRCURV_CONVEX;
462  else if( nnegative == nfactors-1 && EPSGE(sum, 1.0, 1e-9) && expcurvpos )
463  curv = SCIP_EXPRCURV_CONVEX;
464  else if( npositive == nfactors && EPSLE(sum, 1.0, 1e-9) && expcurvneg )
465  curv = SCIP_EXPRCURV_CONCAVE;
466  else
467  curv = SCIP_EXPRCURV_UNKNOWN;
468  curv = SCIPexprcurvMultiply(mult, curv);
469 
470  return curv;
471 }
472 
473 /** gives name as string for a curvature */
475  SCIP_EXPRCURV curv /**< curvature */
476  )
477 {
478  assert(curv <= SCIP_EXPRCURV_LINEAR); /*lint !e685*/
479 
480  return curvnames[curv];
481 }
482 
483 /**@} */
484 
485 /**@name Quadratic expression data private methods */
486 /**@{ */
487 
488 /** creates SCIP_EXPRDATA_QUADRATIC data structure from given quadratic elements */
489 static
491  BMS_BLKMEM* blkmem, /**< block memory data structure */
492  SCIP_EXPRDATA_QUADRATIC** quadraticdata, /**< buffer to store pointer to quadratic data */
493  SCIP_Real constant, /**< constant */
494  int nchildren, /**< number of children */
495  SCIP_Real* lincoefs, /**< linear coefficients of children, or NULL if all 0.0 */
496  int nquadelems, /**< number of quadratic elements */
497  SCIP_QUADELEM* quadelems /**< quadratic elements */
498  )
499 {
500  assert(blkmem != NULL);
501  assert(quadraticdata != NULL);
502  assert(quadelems != NULL || nquadelems == 0);
503  assert(nchildren >= 0);
504 
505  SCIP_ALLOC( BMSallocBlockMemory(blkmem, quadraticdata) );
506 
507  (*quadraticdata)->constant = constant;
508  (*quadraticdata)->lincoefs = NULL;
509  (*quadraticdata)->nquadelems = nquadelems;
510  (*quadraticdata)->quadelems = NULL;
511  (*quadraticdata)->sorted = (nquadelems <= 1);
512 
513  if( lincoefs != NULL )
514  {
515  SCIP_ALLOC( BMSduplicateBlockMemoryArray(blkmem, &(*quadraticdata)->lincoefs, lincoefs, nchildren) );
516  }
517 
518  if( nquadelems > 0 )
519  {
520  SCIP_ALLOC( BMSduplicateBlockMemoryArray(blkmem, &(*quadraticdata)->quadelems, quadelems, nquadelems) );
521  }
522 
523  return SCIP_OKAY;
524 }
525 
526 /** sorts quadratic elements in a SCIP_EXPRDATA_QUADRATIC data structure */
527 static
529  SCIP_EXPRDATA_QUADRATIC* quadraticdata /**< quadratic data */
530  )
531 {
532  assert(quadraticdata != NULL);
533 
534  if( quadraticdata->sorted )
535  {
536 #ifndef NDEBUG
537  int i;
538  for( i = 1; i < quadraticdata->nquadelems; ++i )
539  {
540  assert(quadraticdata->quadelems[i].idx1 <= quadraticdata->quadelems[i].idx2);
541  assert(quadraticdata->quadelems[i-1].idx1 <= quadraticdata->quadelems[i].idx1);
542  assert(quadraticdata->quadelems[i-1].idx1 < quadraticdata->quadelems[i].idx1 ||
543  quadraticdata->quadelems[i-1].idx2 <= quadraticdata->quadelems[i].idx2);
544  }
545 #endif
546  return;
547  }
548 
549  if( quadraticdata->nquadelems > 0 )
550  SCIPquadelemSort(quadraticdata->quadelems, quadraticdata->nquadelems);
551 
552  quadraticdata->sorted = TRUE;
553 }
554 
555 /**@} */
556 
557 /**@name Polynomial expression data private methods */
558 /**@{ */
559 
560 /** compares two monomials
561  *
562  * gives 0 if monomials are equal */
563 static
564 SCIP_DECL_SORTPTRCOMP(monomialdataCompare)
565 {
566  SCIP_EXPRDATA_MONOMIAL* monomial1;
567  SCIP_EXPRDATA_MONOMIAL* monomial2;
568 
569  int i;
570 
571  assert(elem1 != NULL);
572  assert(elem2 != NULL);
573 
574  monomial1 = (SCIP_EXPRDATA_MONOMIAL*)elem1;
575  monomial2 = (SCIP_EXPRDATA_MONOMIAL*)elem2;
576 
577  /* make sure, both monomials are equal */
578  SCIPexprSortMonomialFactors(monomial1);
579  SCIPexprSortMonomialFactors(monomial2);
580 
581  /* for the first factor where both monomials differ,
582  * we return either the difference in the child indices, if children are different
583  * or the sign of the difference in the exponents
584  */
585  for( i = 0; i < monomial1->nfactors && i < monomial2->nfactors; ++i )
586  {
587  if( monomial1->childidxs[i] != monomial2->childidxs[i] )
588  return monomial1->childidxs[i] - monomial2->childidxs[i];
589  if( monomial1->exponents[i] > monomial2->exponents[i] )
590  return 1;
591  else if( monomial1->exponents[i] < monomial2->exponents[i] )
592  return -1;
593  }
594 
595  /* if the factors of one monomial are a proper subset of the factors of the other monomial,
596  * we return the difference in the number of monomials
597  */
598  return monomial1->nfactors - monomial2->nfactors;
599 }
600 
601 /** ensures that the factors arrays of a monomial have at least a given size */
602 static
604  BMS_BLKMEM* blkmem, /**< block memory data structure */
605  SCIP_EXPRDATA_MONOMIAL* monomialdata, /**< monomial data */
606  int minsize /**< minimal size of factors arrays */
607  )
608 {
609  assert(blkmem != NULL);
610  assert(monomialdata != NULL);
611 
612  if( minsize > monomialdata->factorssize )
613  {
614  int newsize;
615 
616  newsize = calcGrowSize(minsize);
617  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &monomialdata->childidxs, monomialdata->factorssize, newsize) );
618  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &monomialdata->exponents, monomialdata->factorssize, newsize) );
619  monomialdata->factorssize = newsize;
620  }
621  assert(minsize <= monomialdata->factorssize);
622 
623  return SCIP_OKAY;
624 }
625 
626 /** creates SCIP_EXPRDATA_POLYNOMIAL data structure from given monomials */
627 static
629  BMS_BLKMEM* blkmem, /**< block memory data structure */
630  SCIP_EXPRDATA_POLYNOMIAL** polynomialdata,/**< buffer to store pointer to polynomial data */
631  int nmonomials, /**< number of monomials */
632  SCIP_EXPRDATA_MONOMIAL** monomials, /**< monomials */
633  SCIP_Real constant, /**< constant part */
634  SCIP_Bool copymonomials /**< whether to copy monomials, or copy only given pointers, in which case polynomialdata assumes ownership of monomial structure */
635  )
636 {
637  assert(blkmem != NULL);
638  assert(polynomialdata != NULL);
639  assert(monomials != NULL || nmonomials == 0);
640 
641  SCIP_ALLOC( BMSallocBlockMemory(blkmem, polynomialdata) );
642 
643  (*polynomialdata)->constant = constant;
644  (*polynomialdata)->nmonomials = nmonomials;
645  (*polynomialdata)->monomialssize = nmonomials;
646  (*polynomialdata)->monomials = NULL;
647  (*polynomialdata)->sorted = (nmonomials <= 1);
648 
649  if( nmonomials > 0 )
650  {
651  int i;
652 
653  if( copymonomials )
654  {
655  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &(*polynomialdata)->monomials, nmonomials) );
656 
657  for( i = 0; i < nmonomials; ++i )
658  {
659  assert(monomials[i] != NULL); /*lint !e613*/
660  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &(*polynomialdata)->monomials[i],
661  monomials[i]->coef, monomials[i]->nfactors, monomials[i]->childidxs, monomials[i]->exponents) ); /*lint !e613*/
662  }
663  }
664  else
665  {
666  SCIP_ALLOC( BMSduplicateBlockMemoryArray(blkmem, &(*polynomialdata)->monomials, monomials, nmonomials) );
667  }
668  }
669 
670  return SCIP_OKAY;
671 }
672 
673 /** creates a copy of a SCIP_EXPRDATA_POLYNOMIAL data structure */
674 static
676  BMS_BLKMEM* blkmem, /**< block memory data structure */
677  SCIP_EXPRDATA_POLYNOMIAL** polynomialdata,/**< buffer to store pointer to polynomial data */
678  SCIP_EXPRDATA_POLYNOMIAL* sourcepolynomialdata /**< polynomial data to copy */
679  )
680 {
681  assert(blkmem != NULL);
682  assert(polynomialdata != NULL);
683  assert(sourcepolynomialdata != NULL);
684 
685  SCIP_ALLOC( BMSduplicateBlockMemory(blkmem, polynomialdata, sourcepolynomialdata) );
686 
687  (*polynomialdata)->monomialssize = sourcepolynomialdata->nmonomials;
688  if( sourcepolynomialdata->nmonomials > 0 )
689  {
690  int i;
691 
692  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &(*polynomialdata)->monomials, (*polynomialdata)->monomialssize) );
693 
694  for( i = 0; i < sourcepolynomialdata->nmonomials; ++i )
695  {
696  assert(sourcepolynomialdata->monomials[i] != NULL); /*lint !e613*/
697  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &(*polynomialdata)->monomials[i], sourcepolynomialdata->monomials[i]->coef,
698  sourcepolynomialdata->monomials[i]->nfactors, sourcepolynomialdata->monomials[i]->childidxs, sourcepolynomialdata->monomials[i]->exponents) );
699  (*polynomialdata)->monomials[i]->sorted = sourcepolynomialdata->monomials[i]->sorted;
700  }
701  }
702  else
703  {
704  (*polynomialdata)->monomials = NULL;
705  }
706 
707  return SCIP_OKAY;
708 }
709 
710 /** frees a SCIP_EXPRDATA_POLYNOMIAL data structure */
711 static
713  BMS_BLKMEM* blkmem, /**< block memory data structure */
714  SCIP_EXPRDATA_POLYNOMIAL** polynomialdata /**< pointer to polynomial data to free */
715  )
716 {
717  assert(blkmem != NULL);
718  assert(polynomialdata != NULL);
719  assert(*polynomialdata != NULL);
720 
721  if( (*polynomialdata)->monomialssize > 0 )
722  {
723  int i;
724 
725  for( i = 0; i < (*polynomialdata)->nmonomials; ++i )
726  {
727  assert((*polynomialdata)->monomials[i] != NULL);
728  SCIPexprFreeMonomial(blkmem, &(*polynomialdata)->monomials[i]);
729  assert((*polynomialdata)->monomials[i] == NULL);
730  }
731 
732  BMSfreeBlockMemoryArray(blkmem, &(*polynomialdata)->monomials, (*polynomialdata)->monomialssize);
733  }
734  assert((*polynomialdata)->monomials == NULL);
735 
736  BMSfreeBlockMemory(blkmem, polynomialdata);
737 }
738 
739 /** ensures that the monomials array of a polynomial has at least a given size */
740 static
742  BMS_BLKMEM* blkmem, /**< block memory data structure */
743  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
744  int minsize /**< minimal size of monomials array */
745  )
746 {
747  assert(blkmem != NULL);
748  assert(polynomialdata != NULL);
749 
750  ensureBlockMemoryArraySize(blkmem, &polynomialdata->monomials, &polynomialdata->monomialssize, minsize);
751  assert(minsize <= polynomialdata->monomialssize);
752 
753  return SCIP_OKAY;
754 }
755 
756 /** adds an array of monomials to a polynomial */
757 static
759  BMS_BLKMEM* blkmem, /**< block memory of expression */
760  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
761  int nmonomials, /**< number of monomials to add */
762  SCIP_EXPRDATA_MONOMIAL** monomials, /**< the monomials to add */
763  SCIP_Bool copymonomials /**< whether to copy monomials or to assume ownership */
764  )
765 {
766  int i;
767 
768  assert(blkmem != NULL);
769  assert(polynomialdata != NULL);
770  assert(monomials != NULL || nmonomials == 0);
771 
772  if( nmonomials == 0 )
773  return SCIP_OKAY;
774 
775  SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, polynomialdata->nmonomials + nmonomials) );
776  assert(polynomialdata->monomialssize >= polynomialdata->nmonomials + nmonomials);
777 
778  if( copymonomials )
779  {
780  for( i = 0; i < nmonomials; ++i )
781  {
782  assert(monomials[i] != NULL); /*lint !e613*/
783  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &polynomialdata->monomials[polynomialdata->nmonomials + i],
784  monomials[i]->coef, monomials[i]->nfactors, monomials[i]->childidxs, monomials[i]->exponents) ); /*lint !e613*/
785  }
786  }
787  else
788  {
789  BMScopyMemoryArray(&polynomialdata->monomials[polynomialdata->nmonomials], monomials, nmonomials); /*lint !e866*/
790  }
791  polynomialdata->nmonomials += nmonomials;
792 
793  polynomialdata->sorted = (polynomialdata->nmonomials <= 1);
794 
795  return SCIP_OKAY;
796 }
797 
798 /** ensures that monomials of a polynomial are sorted */
799 static
801  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata /**< polynomial expression */
802  )
803 {
804  assert(polynomialdata != NULL);
805 
806  if( polynomialdata->sorted )
807  {
808 #ifndef NDEBUG
809  int i;
810 
811  /* a polynom with more than one monoms can only be sorted if its monoms are sorted */
812  for( i = 1; i < polynomialdata->nmonomials; ++i )
813  {
814  assert(polynomialdata->monomials[i-1]->sorted);
815  assert(polynomialdata->monomials[i]->sorted);
816  assert(monomialdataCompare(polynomialdata->monomials[i-1], polynomialdata->monomials[i]) <= 0);
817  }
818 #endif
819  return;
820  }
821 
822  if( polynomialdata->nmonomials > 0 )
823  SCIPsortPtr((void**)polynomialdata->monomials, monomialdataCompare, polynomialdata->nmonomials);
824 
825  polynomialdata->sorted = TRUE;
826 }
827 
828 /** merges monomials that differ only in coefficient into a single monomial
829  *
830  * Eliminates monomials with coefficient between -eps and eps.
831  */
832 static
834  BMS_BLKMEM* blkmem, /**< block memory */
835  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
836  SCIP_Real eps, /**< threshold under which numbers are treat as zero */
837  SCIP_Bool mergefactors /**< whether to merge factors in monomials too */
838  )
839 {
840  int i;
841  int offset;
842  int oldnfactors;
843 
844  assert(polynomialdata != NULL);
845  assert(eps >= 0.0);
846 
847  polynomialdataSortMonomials(polynomialdata);
848 
849  /* merge monomials by adding their coefficients, eliminate monomials with no factors or zero coefficient*/
850  offset = 0;
851  i = 0;
852  while( i + offset < polynomialdata->nmonomials )
853  {
854  if( offset > 0 )
855  {
856  assert(polynomialdata->monomials[i] == NULL);
857  assert(polynomialdata->monomials[i+offset] != NULL);
858  polynomialdata->monomials[i] = polynomialdata->monomials[i+offset];
859 #ifndef NDEBUG
860  polynomialdata->monomials[i+offset] = NULL;
861 #endif
862  }
863 
864  if( mergefactors )
865  {
866  oldnfactors = polynomialdata->monomials[i]->nfactors;
867  SCIPexprMergeMonomialFactors(polynomialdata->monomials[i], eps);
868 
869  /* if monomial has changed, then we cannot assume anymore that polynomial is sorted */
870  if( oldnfactors != polynomialdata->monomials[i]->nfactors )
871  polynomialdata->sorted = FALSE;
872  }
873 
874  while( i+offset+1 < polynomialdata->nmonomials )
875  {
876  assert(polynomialdata->monomials[i+offset+1] != NULL);
877  if( mergefactors )
878  {
879  oldnfactors = polynomialdata->monomials[i+offset+1]->nfactors;
880  SCIPexprMergeMonomialFactors(polynomialdata->monomials[i+offset+1], eps);
881 
882  /* if monomial has changed, then we cannot assume anymore that polynomial is sorted */
883  if( oldnfactors != polynomialdata->monomials[i+offset+1]->nfactors )
884  polynomialdata->sorted = FALSE;
885  }
886  if( monomialdataCompare((void*)polynomialdata->monomials[i], (void*)polynomialdata->monomials[i+offset+1]) != 0 )
887  break;
888  polynomialdata->monomials[i]->coef += polynomialdata->monomials[i+offset+1]->coef;
889  SCIPexprFreeMonomial(blkmem, &polynomialdata->monomials[i+offset+1]);
890  ++offset;
891  }
892 
893  if( polynomialdata->monomials[i]->nfactors == 0 )
894  {
895  /* constant monomial */
896  polynomialdata->constant += polynomialdata->monomials[i]->coef;
897  SCIPexprFreeMonomial(blkmem, &polynomialdata->monomials[i]);
898  ++offset;
899  continue;
900  }
901 
902  if( EPSZ(polynomialdata->monomials[i]->coef, eps) )
903  {
904  SCIPexprFreeMonomial(blkmem, &polynomialdata->monomials[i]);
905  ++offset;
906  continue;
907  }
908 
909  ++i;
910  }
911 
912 #ifndef NDEBUG
913  for( ; i < polynomialdata->nmonomials; ++i )
914  assert(polynomialdata->monomials[i] == NULL);
915 #endif
916 
917  polynomialdata->nmonomials -= offset;
918 
919  if( EPSZ(polynomialdata->constant, eps) )
920  polynomialdata->constant = 0.0;
921 }
922 
923 /** multiplies each summand of a polynomial by a given constant */
924 static
926  BMS_BLKMEM* blkmem, /**< block memory */
927  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
928  SCIP_Real factor /**< constant factor */
929  )
930 {
931  int i;
932 
933  assert(polynomialdata != NULL);
934 
935  if( factor == 1.0 )
936  return;
937 
938  if( factor == 0.0 )
939  {
940  for( i = 0; i < polynomialdata->nmonomials; ++i )
941  SCIPexprFreeMonomial(blkmem, &polynomialdata->monomials[i]);
942  polynomialdata->nmonomials = 0;
943  }
944  else
945  {
946  for( i = 0; i < polynomialdata->nmonomials; ++i )
947  SCIPexprChgMonomialCoef(polynomialdata->monomials[i], polynomialdata->monomials[i]->coef * factor);
948  }
949 
950  polynomialdata->constant *= factor;
951 }
952 
953 /** multiplies each summand of a polynomial by a given monomial */
954 static
956  BMS_BLKMEM* blkmem, /**< block memory */
957  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
958  SCIP_EXPRDATA_MONOMIAL* factor, /**< monomial factor */
959  int* childmap /**< map children in factor to children in expr, or NULL for 1:1 */
960  )
961 {
962  int i;
963 
964  assert(blkmem != NULL);
965  assert(factor != NULL);
966  assert(polynomialdata != NULL);
967 
968  if( factor->nfactors == 0 )
969  {
970  polynomialdataMultiplyByConstant(blkmem, polynomialdata, factor->coef);
971  return SCIP_OKAY;
972  }
973 
974  /* multiply each monomial by factor */
975  for( i = 0; i < polynomialdata->nmonomials; ++i )
976  {
977  SCIP_CALL( SCIPexprMultiplyMonomialByMonomial(blkmem, polynomialdata->monomials[i], factor, childmap) );
978  }
979 
980  /* add new monomial for constant multiplied by factor */
981  if( polynomialdata->constant != 0.0 )
982  {
983  SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, polynomialdata->nmonomials+1) );
984  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &polynomialdata->monomials[polynomialdata->nmonomials], polynomialdata->constant, 0, NULL, NULL) );
985  SCIP_CALL( SCIPexprMultiplyMonomialByMonomial(blkmem, polynomialdata->monomials[polynomialdata->nmonomials], factor, childmap) );
986  ++polynomialdata->nmonomials;
987  polynomialdata->sorted = FALSE;
988  polynomialdata->constant = 0.0;
989  }
990 
991  return SCIP_OKAY;
992 }
993 
994 /** multiplies a polynomial by a polynomial
995  *
996  * Factors need to be different.
997  */
998 static
1000  BMS_BLKMEM* blkmem, /**< block memory */
1001  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
1002  SCIP_EXPRDATA_POLYNOMIAL* factordata, /**< polynomial factor data */
1003  int* childmap /**< map children in factor to children in polynomialdata, or NULL for 1:1 */
1004  )
1005 {
1006  int i1;
1007  int i2;
1008  int orignmonomials;
1009 
1010  assert(blkmem != NULL);
1011  assert(polynomialdata != NULL);
1012  assert(factordata != NULL);
1013  assert(polynomialdata != factordata);
1014 
1015  if( factordata->nmonomials == 0 )
1016  {
1017  polynomialdataMultiplyByConstant(blkmem, polynomialdata, factordata->constant);
1018  return SCIP_OKAY;
1019  }
1020 
1021  if( factordata->nmonomials == 1 && factordata->constant == 0.0 )
1022  {
1023  SCIP_CALL( polynomialdataMultiplyByMonomial(blkmem, polynomialdata, factordata->monomials[0], childmap) );
1024  return SCIP_OKAY;
1025  }
1026 
1027  /* turn constant into a monomial, so we can assume below that constant is 0.0 */
1028  if( polynomialdata->constant != 0.0 )
1029  {
1030  SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, polynomialdata->nmonomials+1) );
1031  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &polynomialdata->monomials[polynomialdata->nmonomials], polynomialdata->constant, 0, NULL, NULL) );
1032  ++polynomialdata->nmonomials;
1033  polynomialdata->sorted = FALSE;
1034  polynomialdata->constant = 0.0;
1035  }
1036 
1037  SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, polynomialdata->nmonomials * (factordata->nmonomials + (factordata->constant == 0.0 ? 0 : 1))) );
1038 
1039  /* for each monomial in factordata (except the last, if factordata->constant is 0),
1040  * duplicate monomials from polynomialdata and multiply them by the monomial for factordata */
1041  orignmonomials = polynomialdata->nmonomials;
1042  for( i2 = 0; i2 < factordata->nmonomials; ++i2 )
1043  {
1044  /* add a copy of original monomials to end of polynomialdata's monomials array */
1045  assert(polynomialdata->nmonomials + orignmonomials <= polynomialdata->monomialssize); /* reallocating in polynomialdataAddMonomials would make the polynomialdata->monomials invalid, so assert that above the monomials array was made large enough */
1046  SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, orignmonomials, polynomialdata->monomials, TRUE) );
1047  assert(polynomialdata->nmonomials == (i2+2) * orignmonomials);
1048 
1049  /* multiply each copied monomial by current monomial from factordata */
1050  for( i1 = (i2+1) * orignmonomials; i1 < (i2+2) * orignmonomials; ++i1 )
1051  {
1052  SCIP_CALL( SCIPexprMultiplyMonomialByMonomial(blkmem, polynomialdata->monomials[i1], factordata->monomials[i2], childmap) );
1053  }
1054 
1055  if( factordata->constant == 0.0 && i2 == factordata->nmonomials - 2 )
1056  {
1057  ++i2;
1058  break;
1059  }
1060  }
1061 
1062  if( factordata->constant != 0.0 )
1063  {
1064  assert(i2 == factordata->nmonomials);
1065  /* multiply original monomials in polynomialdata by constant in factordata */
1066  for( i1 = 0; i1 < orignmonomials; ++i1 )
1067  SCIPexprChgMonomialCoef(polynomialdata->monomials[i1], polynomialdata->monomials[i1]->coef * factordata->constant);
1068  }
1069  else
1070  {
1071  assert(i2 == factordata->nmonomials - 1);
1072  /* multiply original monomials in polynomialdata by last monomial in factordata */
1073  for( i1 = 0; i1 < orignmonomials; ++i1 )
1074  {
1075  SCIP_CALL( SCIPexprMultiplyMonomialByMonomial(blkmem, polynomialdata->monomials[i1], factordata->monomials[i2], childmap) );
1076  }
1077  }
1078 
1079  return SCIP_OKAY;
1080 }
1081 
1082 /** takes a power of a polynomial
1083  *
1084  * Exponent needs to be an integer,
1085  * polynomial needs to be a monomial, if exponent is negative.
1086  */
1087 static
1089  BMS_BLKMEM* blkmem, /**< block memory */
1090  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
1091  int exponent /**< exponent of power operation */
1092  )
1093 {
1094  SCIP_EXPRDATA_POLYNOMIAL* factor;
1095  int i;
1096 
1097  assert(blkmem != NULL);
1098  assert(polynomialdata != NULL);
1099 
1100  if( exponent == 0 )
1101  {
1102  /* x^0 = 1, except if x = 0 */
1103  if( polynomialdata->nmonomials == 0 && polynomialdata->constant == 0.0 )
1104  {
1105  polynomialdata->constant = 0.0;
1106  }
1107  else
1108  {
1109  polynomialdata->constant = 1.0;
1110 
1111  for( i = 0; i < polynomialdata->nmonomials; ++i )
1112  SCIPexprFreeMonomial(blkmem, &polynomialdata->monomials[i]);
1113  polynomialdata->nmonomials = 0;
1114  }
1115 
1116  return SCIP_OKAY;
1117  }
1118 
1119  if( exponent == 1 )
1120  return SCIP_OKAY;
1121 
1122  if( polynomialdata->nmonomials == 1 && polynomialdata->constant == 0.0 )
1123  {
1124  /* polynomial is a single monomial */
1125  SCIPexprMonomialPower(polynomialdata->monomials[0], exponent);
1126  return SCIP_OKAY;
1127  }
1128 
1129  if( polynomialdata->nmonomials == 0 )
1130  {
1131  /* polynomial is a constant */
1132  polynomialdata->constant = pow(polynomialdata->constant, (SCIP_Real)exponent);
1133  return SCIP_OKAY;
1134  }
1135 
1136  assert(exponent >= 2); /* negative exponents not allowed if more than one monom */
1137 
1138  /* todo improve, look how SCIPintervalPowerScalar in intervalarith.c does it */
1139 
1140  /* get copy of our polynomial */
1141  SCIP_CALL( polynomialdataCopy(blkmem, &factor, polynomialdata) );
1142 
1143  /* do repeated multiplication */
1144  for( i = 2; i <= exponent; ++i )
1145  {
1146  SCIP_CALL( polynomialdataMultiplyByPolynomial(blkmem, polynomialdata, factor, NULL) );
1147  polynomialdataMergeMonomials(blkmem, polynomialdata, 0.0, TRUE);
1148  }
1149 
1150  /* free copy again */
1151  polynomialdataFree(blkmem, &factor);
1152 
1153  return SCIP_OKAY;
1154 }
1155 
1156 /** applies a mapping of child indices to the indices used in polynomial monomials */
1157 static
1159  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
1160  int* childmap /**< mapping of child indices */
1161  )
1162 {
1163  SCIP_EXPRDATA_MONOMIAL* monomial;
1164  int i;
1165  int j;
1166 
1167  assert(polynomialdata != NULL);
1168 
1169  for( i = 0; i < polynomialdata->nmonomials; ++i )
1170  {
1171  monomial = polynomialdata->monomials[i];
1172  assert(monomial != NULL);
1173 
1174  for( j = 0; j < monomial->nfactors; ++j )
1175  {
1176  monomial->childidxs[j] = childmap[monomial->childidxs[j]];
1177  assert(monomial->childidxs[j] >= 0);
1178  }
1179  monomial->sorted = FALSE;
1180  }
1181 
1182  polynomialdata->sorted = FALSE;
1183 }
1184 
1185 /** replaces a factor in a monomial by a polynomial and expands the result */
1186 static
1188  BMS_BLKMEM* blkmem, /**< block memory data structure */
1189  SCIP_MESSAGEHDLR* messagehdlr, /**< message handler */
1190  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data where to expand a monomial */
1191  int monomialpos, /**< position of monomial which factor to expand */
1192  int factorpos, /**< position of factor in monomial to expand */
1193  SCIP_EXPRDATA_POLYNOMIAL* factorpolynomial,/**< polynomial that should replace factor */
1194  int* childmap, /**< map of child indices in factorpolynomial to children of polynomial */
1195  int maxexpansionexponent,/**< maximal exponent for which polynomials (with > 1 summands) are expanded */
1196  SCIP_Bool* success /**< buffer to store whether expansion has been done */
1197  )
1198 {
1199  SCIP_EXPRDATA_POLYNOMIAL* factorpolynomialcopy;
1200  SCIP_EXPRDATA_MONOMIAL* monomial;
1201  int i;
1202 
1203  assert(blkmem != NULL);
1204  assert(polynomialdata != NULL);
1205  assert(factorpolynomial != NULL);
1206  assert(childmap != NULL || factorpolynomial->nmonomials == 0);
1207  assert(success != NULL);
1208  assert(monomialpos >= 0);
1209  assert(monomialpos < polynomialdata->nmonomials);
1210  assert(factorpos >= 0);
1211 
1212  monomial = polynomialdata->monomials[monomialpos];
1213  assert(monomial != NULL);
1214  assert(factorpos < monomial->nfactors);
1215 
1216  *success = TRUE;
1217 
1218  if( factorpolynomial->nmonomials == 0 )
1219  {
1220  /* factorpolynomial is a constant */
1221 
1222  if( !EPSISINT(monomial->exponents[factorpos], 0.0) && factorpolynomial->constant < 0.0 ) /*lint !e835*/
1223  {
1224  /* if polynomial is a negative constant and our exponent is not integer, then cannot do expansion */
1225  SCIPmessagePrintWarning(messagehdlr, "got negative constant %g to the power of a noninteger exponent %g\n", factorpolynomial->constant, monomial->exponents[factorpos]);
1226  *success = FALSE;
1227  return SCIP_OKAY;
1228  }
1229  monomial->coef *= pow(factorpolynomial->constant, monomial->exponents[factorpos]);
1230 
1231  /* move last factor to position factorpos */
1232  if( factorpos < monomial->nfactors-1 )
1233  {
1234  monomial->exponents[factorpos] = monomial->exponents[monomial->nfactors-1];
1235  monomial->childidxs[factorpos] = monomial->childidxs[monomial->nfactors-1];
1236  }
1237  --monomial->nfactors;
1238  monomial->sorted = FALSE;
1239  polynomialdata->sorted = FALSE;
1240 
1241  return SCIP_OKAY;
1242  }
1243 
1244  if( factorpolynomial->constant == 0.0 && factorpolynomial->nmonomials == 1 )
1245  {
1246  /* factorpolynomial is a single monomial */
1247  SCIP_EXPRDATA_MONOMIAL* factormonomial;
1248  int childidx;
1249  SCIP_Real exponent;
1250 
1251  factormonomial = factorpolynomial->monomials[0];
1252  assert(factormonomial != NULL);
1253 
1254  if( !EPSISINT(monomial->exponents[factorpos], 0.0) ) /*lint !e835*/
1255  {
1256  if( factormonomial->coef < 0.0 )
1257  {
1258  /* if coefficient of monomial is negative and our exponent is not integer, then do not do expansion
1259  * @todo the only case where this could make sense is if the factors can be negative, i.e., when we have negative arguments with an odd exponent: (-x^a)^b = (-x)^(ab) for a odd
1260  */
1261  *success = FALSE;
1262  return SCIP_OKAY;
1263  }
1264  if( factormonomial->nfactors > 1 )
1265  {
1266  /* @todo if there is an even number of factors in factormonomial that are negative, then they always multiply to something positive
1267  * however, we cannot expand them as below, since we cannot compute the single powers
1268  * since we do not have the bounds on the factors here, we skip expansion in this case
1269  * MINLPLib instances tls2,4,6 are examples where we are loosing here (do not recognize convexity)
1270  */
1271  *success = FALSE;
1272  return SCIP_OKAY;
1273  }
1274  }
1275 
1276  SCIP_CALL( monomialdataEnsureFactorsSize(blkmem, monomial, monomial->nfactors + factormonomial->nfactors) );
1277 
1278  for( i = 0; i < factormonomial->nfactors; ++i )
1279  {
1280  childidx = childmap[factormonomial->childidxs[i]]; /*lint !e613*/
1281  /* can do this because monomial->exponents[factorpos] is assumed to be integer or factormonomial has positive coefficient and only one factor
1282  * thus, if factormonomial->exponents[i] is fractional, then we can assume that it's argument is positive
1283  */
1284  exponent = factormonomial->exponents[i] * monomial->exponents[factorpos];
1285  SCIP_CALL( SCIPexprAddMonomialFactors(blkmem, monomial, 1, &childidx, &exponent) );
1286  }
1287 
1288  monomial->coef *= pow(factormonomial->coef, monomial->exponents[factorpos]);
1289 
1290  /* move last factor to position factorpos */
1291  if( factorpos < monomial->nfactors-1 )
1292  {
1293  monomial->exponents[factorpos] = monomial->exponents[monomial->nfactors-1];
1294  monomial->childidxs[factorpos] = monomial->childidxs[monomial->nfactors-1];
1295  }
1296  --monomial->nfactors;
1297  monomial->sorted = FALSE;
1298  polynomialdata->sorted = FALSE;
1299 
1300  return SCIP_OKAY;
1301  }
1302 
1303  /* if exponent is negative or fractional and the polynomial is not just a monomial, then we cannot do expansion */
1304  if( !EPSISINT(monomial->exponents[factorpos], 0.0) || monomial->exponents[factorpos] < 0.0 ) /*lint !e835*/
1305  {
1306  *success = FALSE;
1307  return SCIP_OKAY;
1308  }
1309 
1310  /* if exponent is too large, skip expansion */
1311  if( monomial->exponents[factorpos] > maxexpansionexponent )
1312  {
1313  *success = FALSE;
1314  return SCIP_OKAY;
1315  }
1316 
1317  /* check whether maximal degree of expansion would exceed maxexpansionexponent
1318  * that is, assume monomial is f1^a1 f2^a2 ... and we want to expand f1 = (g11^beta11 g12^beta12... + g21^beta21 g22^beta22 ... + ...)
1319  * then we do this only if all ai and all beta are > 0.0 and a1 max(beta11+beta12+..., beta21+beta22+..., ...) + a2 + ... < maxexpansionexponent
1320  * exception (there need to be one) is if monomial is just f1
1321  */
1322  if( maxexpansionexponent < INT_MAX && (monomial->nfactors > 1 || monomial->exponents[factorpos] != 1.0) )
1323  {
1324  SCIP_Real restdegree;
1325  SCIP_Real degree;
1326  int j;
1327 
1328  restdegree = -monomial->exponents[factorpos];
1329  for( i = 0; i < monomial->nfactors; ++i )
1330  {
1331  if( monomial->exponents[i] < 0.0 )
1332  {
1333  /* ai < 0.0 */
1334  SCIPdebugMessage("skip expansion because factor %d in monomial has negative exponent\n", i);
1335  *success = FALSE;
1336  return SCIP_OKAY;
1337  }
1338  restdegree += monomial->exponents[i];
1339  }
1340 
1341  for( i = 0; i < factorpolynomial->nmonomials; ++i )
1342  {
1343  degree = 0.0;
1344  for( j = 0; j < factorpolynomial->monomials[i]->nfactors; ++j )
1345  {
1346  if( factorpolynomial->monomials[i]->exponents[j] < 0.0 )
1347  {
1348  /* beta_ij < 0.0 */
1349  SCIPdebugMessage("skip expansion because %d'th factor in %d'th monomial of factorpolynomial is negative\n", i, j);
1350  *success = FALSE;
1351  return SCIP_OKAY;
1352  }
1353  degree += factorpolynomial->monomials[i]->exponents[j];
1354  }
1355  if( degree * monomial->exponents[factorpos] + restdegree > maxexpansionexponent )
1356  {
1357  /* (beta_i1+beta_i2+...)*monomial->exponents[factorpos] + rest > maxexpansion */
1358  SCIPdebugMessage("skip expansion because degree of %d'th monomial would yield degree %g > max = %d in expansion\n",
1359  i, degree * monomial->exponents[factorpos] + restdegree, maxexpansionexponent);
1360  *success = FALSE;
1361  return SCIP_OKAY;
1362  }
1363  }
1364  }
1365 
1366  /* create a copy of factor */
1367  SCIP_CALL( polynomialdataCopy(blkmem, &factorpolynomialcopy, factorpolynomial) );
1368  /* apply childmap to copy */
1369  polynomialdataApplyChildmap(factorpolynomialcopy, childmap);
1370  /* create power of factor */
1371  SCIP_CALL( polynomialdataPower(blkmem, factorpolynomialcopy, (int)EPSFLOOR(monomial->exponents[factorpos], 0.0)) ); /*lint !e835*/
1372 
1373  /* remove factor from monomial by moving last factor to position factorpos */
1374  if( factorpos < monomial->nfactors-1 )
1375  {
1376  monomial->exponents[factorpos] = monomial->exponents[monomial->nfactors-1];
1377  monomial->childidxs[factorpos] = monomial->childidxs[monomial->nfactors-1];
1378  }
1379  --monomial->nfactors;
1380  monomial->sorted = FALSE;
1381 
1382  /* multiply factor with this reduced monomial */
1383  SCIP_CALL( polynomialdataMultiplyByMonomial(blkmem, factorpolynomialcopy, monomial, NULL) );
1384 
1385  /* remove monomial from polynomial and move last monomial to monomialpos */
1386  SCIPexprFreeMonomial(blkmem, &polynomialdata->monomials[monomialpos]);
1387  if( monomialpos < polynomialdata->nmonomials-1 )
1388  polynomialdata->monomials[monomialpos] = polynomialdata->monomials[polynomialdata->nmonomials-1];
1389  --polynomialdata->nmonomials;
1390  polynomialdata->sorted = FALSE;
1391 
1392  /* add factorpolynomialcopy to polynomial */
1393  SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, factorpolynomialcopy->nmonomials, factorpolynomialcopy->monomials, FALSE) );
1394  polynomialdata->constant += factorpolynomialcopy->constant;
1395 
1396  factorpolynomialcopy->nmonomials = 0;
1397  polynomialdataFree(blkmem, &factorpolynomialcopy);
1398 
1399  return SCIP_OKAY;
1400 }
1401 
1402 /**@} */
1403 
1404 /**@name Expression operand private methods */
1405 /**@{ */
1406 
1407 /** a default implementation of expression interval evaluation that always gives a correct result */
1408 static
1409 SCIP_DECL_EXPRINTEVAL( exprevalIntDefault )
1410 { /*lint --e{715}*/
1412 
1413  return SCIP_OKAY;
1414 }
1415 
1416 /** a default implementation of expression curvature check that always gives a correct result */
1417 static
1418 SCIP_DECL_EXPRCURV( exprcurvDefault )
1419 { /*lint --e{715}*/
1420  *result = SCIP_EXPRCURV_UNKNOWN;
1421 
1422  return SCIP_OKAY;
1423 }
1424 
1425 /** point evaluation for EXPR_VAR */
1426 static
1427 SCIP_DECL_EXPREVAL( exprevalVar )
1428 { /*lint --e{715}*/
1429  assert(result != NULL);
1430  assert(varvals != NULL);
1431 
1432  *result = varvals[opdata.intval];
1433 
1434  return SCIP_OKAY;
1435 }
1436 
1437 /** interval evaluation for EXPR_VAR */
1438 static
1439 SCIP_DECL_EXPRINTEVAL( exprevalIntVar )
1440 { /*lint --e{715}*/
1441  assert(result != NULL);
1442  assert(varvals != NULL);
1443 
1444  *result = varvals[opdata.intval];
1445 
1446  return SCIP_OKAY;
1447 }
1448 
1449 /** curvature for EXPR_VAR */
1450 static
1451 SCIP_DECL_EXPRCURV( exprcurvVar )
1452 { /*lint --e{715}*/
1453  assert(result != NULL);
1454 
1455  *result = SCIP_EXPRCURV_LINEAR;
1456 
1457  return SCIP_OKAY;
1458 }
1459 
1460 /** point evaluation for EXPR_CONST */
1461 static
1462 SCIP_DECL_EXPREVAL( exprevalConst )
1463 { /*lint --e{715}*/
1464  assert(result != NULL);
1465 
1466  *result = opdata.dbl;
1467 
1468  return SCIP_OKAY;
1469 }
1470 
1471 /** interval evaluation for EXPR_CONST */
1472 static
1473 SCIP_DECL_EXPRINTEVAL( exprevalIntConst )
1474 { /*lint --e{715}*/
1475  assert(result != NULL);
1476 
1477  SCIPintervalSet(result, opdata.dbl);
1478 
1479  return SCIP_OKAY;
1480 }
1481 
1482 /** curvature for EXPR_CONST */
1483 static
1484 SCIP_DECL_EXPRCURV( exprcurvConst )
1485 { /*lint --e{715}*/
1486  assert(result != NULL);
1487 
1488  *result = SCIP_EXPRCURV_LINEAR;
1489 
1490  return SCIP_OKAY;
1491 }
1492 
1493 /** point evaluation for EXPR_PARAM */
1494 static
1495 SCIP_DECL_EXPREVAL( exprevalParam )
1496 { /*lint --e{715}*/
1497  assert(result != NULL);
1498  assert(paramvals != NULL );
1499 
1500  *result = paramvals[opdata.intval];
1501 
1502  return SCIP_OKAY;
1503 }
1504 
1505 /** interval evaluation for EXPR_PARAM */
1506 static
1507 SCIP_DECL_EXPRINTEVAL( exprevalIntParam )
1508 { /*lint --e{715}*/
1509  assert(result != NULL);
1510  assert(paramvals != NULL );
1511 
1512  SCIPintervalSet(result, paramvals[opdata.intval]);
1513 
1514  return SCIP_OKAY;
1515 }
1516 
1517 /** curvature for EXPR_PARAM */
1518 static
1519 SCIP_DECL_EXPRCURV( exprcurvParam )
1520 { /*lint --e{715}*/
1521  assert(result != NULL);
1522 
1523  *result = SCIP_EXPRCURV_LINEAR;
1524 
1525  return SCIP_OKAY;
1526 }
1527 
1528 /** point evaluation for EXPR_PLUS */
1529 static
1530 SCIP_DECL_EXPREVAL( exprevalPlus )
1531 { /*lint --e{715}*/
1532  assert(result != NULL);
1533  assert(argvals != NULL);
1534 
1535  *result = argvals[0] + argvals[1];
1536 
1537  return SCIP_OKAY;
1538 }
1539 
1540 /** interval evaluation for EXPR_PLUS */
1541 static
1542 SCIP_DECL_EXPRINTEVAL( exprevalIntPlus )
1543 { /*lint --e{715}*/
1544  assert(result != NULL);
1545  assert(argvals != NULL);
1546 
1547  SCIPintervalAdd(infinity, result, argvals[0], argvals[1]);
1548 
1549  return SCIP_OKAY;
1550 }
1551 
1552 /** curvature for EXPR_PLUS */
1553 static
1554 SCIP_DECL_EXPRCURV( exprcurvPlus )
1555 { /*lint --e{715}*/
1556  assert(result != NULL);
1557  assert(argcurv != NULL);
1558 
1559  *result = SCIPexprcurvAdd(argcurv[0], argcurv[1]);
1560 
1561  return SCIP_OKAY;
1562 }
1563 
1564 /** point evaluation for EXPR_MINUS */
1565 static
1566 SCIP_DECL_EXPREVAL( exprevalMinus )
1567 { /*lint --e{715}*/
1568  assert(result != NULL);
1569  assert(argvals != NULL);
1570 
1571  *result = argvals[0] - argvals[1];
1572 
1573  return SCIP_OKAY;
1574 }
1575 
1576 /** interval evaluation for EXPR_MINUS */
1577 static
1578 SCIP_DECL_EXPRINTEVAL( exprevalIntMinus )
1579 { /*lint --e{715}*/
1580  assert(result != NULL);
1581  assert(argvals != NULL);
1582 
1583  SCIPintervalSub(infinity, result, argvals[0], argvals[1]);
1584 
1585  return SCIP_OKAY;
1586 }
1587 
1588 /** curvature for EXPR_MINUS */
1589 static
1590 SCIP_DECL_EXPRCURV( exprcurvMinus )
1591 { /*lint --e{715}*/
1592  assert(result != NULL);
1593  assert(argcurv != NULL);
1594 
1595  *result = SCIPexprcurvAdd(argcurv[0], SCIPexprcurvNegate(argcurv[1]));
1596 
1597  return SCIP_OKAY;
1598 }
1599 
1600 /** point evaluation for EXPR_MUL */
1601 static
1602 SCIP_DECL_EXPREVAL( exprevalMult )
1603 { /*lint --e{715}*/
1604  assert(result != NULL);
1605  assert(argvals != NULL);
1606 
1607  *result = argvals[0] * argvals[1];
1608 
1609  return SCIP_OKAY;
1610 }
1611 
1612 /** interval evaluation for EXPR_MUL */
1613 static
1614 SCIP_DECL_EXPRINTEVAL( exprevalIntMult )
1615 { /*lint --e{715}*/
1616  assert(result != NULL);
1617  assert(argvals != NULL);
1618 
1619  SCIPintervalMul(infinity, result, argvals[0], argvals[1]);
1620 
1621  return SCIP_OKAY;
1622 }
1623 
1624 /** curvature for EXPR_MUL */
1625 static
1626 SCIP_DECL_EXPRCURV( exprcurvMult )
1627 { /*lint --e{715}*/
1628  assert(result != NULL);
1629  assert(argcurv != NULL);
1630  assert(argbounds != NULL);
1631 
1632  /* if one factor is constant, then product is
1633  * - linear, if constant is 0.0
1634  * - same curvature as other factor, if constant is positive
1635  * - negated curvature of other factor, if constant is negative
1636  *
1637  * if both factors are not constant, then product may not be convex nor concave
1638  */
1639  if( argbounds[1].inf == argbounds[1].sup ) /*lint !e777*/
1640  *result = SCIPexprcurvMultiply(argbounds[1].inf, argcurv[0]);
1641  else if( argbounds[0].inf == argbounds[0].sup ) /*lint !e777*/
1642  *result = SCIPexprcurvMultiply(argbounds[0].inf, argcurv[1]);
1643  else
1644  *result = SCIP_EXPRCURV_UNKNOWN;
1645 
1646  return SCIP_OKAY;
1647 }
1648 
1649 /** point evaluation for EXPR_DIV */
1650 static
1651 #if defined(__GNUC__) && __GNUC__ * 100 + __GNUC_MINOR__ * 10 >= 490 && !defined(__INTEL_COMPILER)
1652 __attribute__((no_sanitize_undefined))
1653 #endif
1654 SCIP_DECL_EXPREVAL( exprevalDiv )
1655 { /*lint --e{715}*/
1656  assert(result != NULL);
1657  assert(argvals != NULL);
1658 
1659  *result = argvals[0] / argvals[1];
1660 
1661  return SCIP_OKAY;
1662 }
1663 
1664 /** interval evaluation for EXPR_DIV */
1665 static
1666 SCIP_DECL_EXPRINTEVAL( exprevalIntDiv )
1667 { /*lint --e{715}*/
1668  assert(result != NULL);
1669  assert(argvals != NULL);
1670 
1671  SCIPintervalDiv(infinity, result, argvals[0], argvals[1]);
1672 
1673  return SCIP_OKAY;
1674 }
1675 
1676 /** curvature for EXPR_DIV */
1677 static
1678 SCIP_DECL_EXPRCURV( exprcurvDiv )
1679 { /*lint --e{715}*/
1680  assert(result != NULL);
1681  assert(argcurv != NULL);
1682  assert(argbounds != NULL);
1683 
1684  /* if denominator is constant, then quotient has curvature sign(denominator) * curv(nominator)
1685  *
1686  * if nominator is a constant, then quotient is
1687  * - sign(nominator) * convex, if denominator is concave and positive
1688  * - sign(nominator) * concave, if denominator is convex and negative
1689  *
1690  * if denominator is positive but convex, then we don't know, e.g.,
1691  * - 1/x^2 is convex for x>=0
1692  * - 1/(1+(x-1)^2) is neither convex nor concave for x >= 0
1693  *
1694  * if both nominator and denominator are not constant, then quotient may not be convex nor concave
1695  */
1696  if( argbounds[1].inf == argbounds[1].sup ) /*lint !e777*/
1697  {
1698  /* denominator is constant */
1699  *result = SCIPexprcurvMultiply(argbounds[1].inf, argcurv[0]);
1700  }
1701  else if( argbounds[0].inf == argbounds[0].sup ) /*lint !e777*/
1702  {
1703  /* nominator is constant */
1704  if( argbounds[1].inf >= 0.0 && (argcurv[1] & SCIP_EXPRCURV_CONCAVE) )
1705  *result = SCIPexprcurvMultiply(argbounds[0].inf, SCIP_EXPRCURV_CONVEX);
1706  else if( argbounds[1].sup <= 0.0 && (argcurv[1] & SCIP_EXPRCURV_CONVEX) )
1707  *result = SCIPexprcurvMultiply(argbounds[0].inf, SCIP_EXPRCURV_CONCAVE);
1708  else
1709  *result = SCIP_EXPRCURV_UNKNOWN;
1710  }
1711  else
1712  {
1713  /* denominator and nominator not constant */
1714  *result = SCIP_EXPRCURV_UNKNOWN;
1715  }
1716 
1717  return SCIP_OKAY;
1718 }
1719 
1720 /** point evaluation for EXPR_SQUARE */
1721 static
1722 SCIP_DECL_EXPREVAL( exprevalSquare )
1723 { /*lint --e{715}*/
1724  assert(result != NULL);
1725  assert(argvals != NULL);
1726 
1727  *result = argvals[0] * argvals[0];
1728 
1729  return SCIP_OKAY;
1730 }
1731 
1732 /** interval evaluation for EXPR_SQUARE */
1733 static
1734 SCIP_DECL_EXPRINTEVAL( exprevalIntSquare )
1735 { /*lint --e{715}*/
1736  assert(result != NULL);
1737  assert(argvals != NULL);
1738 
1739  SCIPintervalSquare(infinity, result, argvals[0]);
1740 
1741  return SCIP_OKAY;
1742 }
1743 
1744 /** curvature for EXPR_SQUARE */
1745 static
1746 SCIP_DECL_EXPRCURV( exprcurvSquare )
1747 { /*lint --e{715}*/
1748  assert(result != NULL);
1749  assert(argcurv != NULL);
1750  assert(argbounds != NULL);
1751 
1752  *result = SCIPexprcurvPower(argbounds[0], argcurv[0], 2.0);
1753 
1754  return SCIP_OKAY;
1755 }
1756 
1757 /** point evaluation for EXPR_SQRT */
1758 static
1759 SCIP_DECL_EXPREVAL( exprevalSquareRoot )
1760 { /*lint --e{715}*/
1761  assert(result != NULL);
1762  assert(argvals != NULL);
1763 
1764  *result = sqrt(argvals[0]);
1765 
1766  return SCIP_OKAY;
1767 }
1768 
1769 /** interval evaluation for EXPR_SQRT */
1770 static
1771 SCIP_DECL_EXPRINTEVAL( exprevalIntSquareRoot )
1772 { /*lint --e{715}*/
1773  assert(result != NULL);
1774  assert(argvals != NULL);
1775 
1776  SCIPintervalSquareRoot(infinity, result, argvals[0]);
1777 
1778  return SCIP_OKAY;
1779 }
1780 
1781 /** curvature for EXPR_SQRT */
1782 static
1783 SCIP_DECL_EXPRCURV( exprcurvSquareRoot )
1784 { /*lint --e{715}*/
1785  assert(result != NULL);
1786  assert(argcurv != NULL);
1787 
1788  /* square-root is concave, if child is concave
1789  * otherwise, we don't know
1790  */
1791 
1792  if( argcurv[0] & SCIP_EXPRCURV_CONCAVE )
1793  *result = SCIP_EXPRCURV_CONCAVE;
1794  else
1795  *result = SCIP_EXPRCURV_UNKNOWN;
1796 
1797  return SCIP_OKAY;
1798 }
1799 
1800 /** point evaluation for EXPR_REALPOWER */
1801 static
1802 SCIP_DECL_EXPREVAL( exprevalRealPower )
1803 { /*lint --e{715}*/
1804  assert(result != NULL);
1805  assert(argvals != NULL);
1806 
1807  *result = pow(argvals[0], opdata.dbl);
1808 
1809  return SCIP_OKAY;
1810 }
1811 
1812 /** interval evaluation for EXPR_REALPOWER */
1813 static
1814 SCIP_DECL_EXPRINTEVAL( exprevalIntRealPower )
1815 { /*lint --e{715}*/
1816  assert(result != NULL);
1817  assert(argvals != NULL);
1818 
1819  SCIPintervalPowerScalar(infinity, result, argvals[0], opdata.dbl);
1820 
1821  return SCIP_OKAY;
1822 }
1823 
1824 /** curvature for EXPR_REALPOWER */
1825 static
1826 SCIP_DECL_EXPRCURV( exprcurvRealPower )
1827 { /*lint --e{715}*/
1828  assert(result != NULL);
1829  assert(argcurv != NULL);
1830  assert(argbounds != NULL);
1831 
1832  *result = SCIPexprcurvPower(argbounds[0], argcurv[0], opdata.dbl);
1833 
1834  return SCIP_OKAY;
1835 }
1836 
1837 /** point evaluation for EXPR_INTPOWER */
1838 static
1839 #if defined(__GNUC__) && __GNUC__ * 100 + __GNUC_MINOR__ * 10 >= 490 && !defined(__INTEL_COMPILER)
1840 __attribute__((no_sanitize_undefined))
1841 #endif
1842 SCIP_DECL_EXPREVAL( exprevalIntPower )
1843 { /*lint --e{715}*/
1844  assert(result != NULL);
1845  assert(argvals != NULL);
1846 
1847  switch( opdata.intval )
1848  {
1849  case -1:
1850  *result = 1.0 / argvals[0];
1851  return SCIP_OKAY;
1852 
1853  case 0:
1854  *result = 1.0;
1855  return SCIP_OKAY;
1856 
1857  case 1:
1858  *result = argvals[0];
1859  return SCIP_OKAY;
1860 
1861  case 2:
1862  *result = argvals[0] * argvals[0];
1863  return SCIP_OKAY;
1864 
1865  default:
1866  *result = pow(argvals[0], (SCIP_Real)opdata.intval);
1867  }
1868 
1869  return SCIP_OKAY;
1870 }
1871 
1872 /** interval evaluation for EXPR_INTPOWER */
1873 static
1874 SCIP_DECL_EXPRINTEVAL( exprevalIntIntPower )
1875 { /*lint --e{715}*/
1876  assert(result != NULL);
1877  assert(argvals != NULL);
1878 
1879  SCIPintervalPowerScalar(infinity, result, argvals[0], (SCIP_Real)opdata.intval);
1880 
1881  return SCIP_OKAY;
1882 }
1883 
1884 /** curvature for EXPR_INTPOWER */
1885 static
1886 SCIP_DECL_EXPRCURV( exprcurvIntPower )
1887 { /*lint --e{715}*/
1888  assert(result != NULL);
1889  assert(argcurv != NULL);
1890  assert(argbounds != NULL);
1891 
1892  *result = SCIPexprcurvPower(argbounds[0], argcurv[0], (SCIP_Real)opdata.intval);
1893 
1894  return SCIP_OKAY;
1895 }
1896 
1897 /** point evaluation for EXPR_SIGNPOWER */
1898 static
1899 SCIP_DECL_EXPREVAL( exprevalSignPower )
1900 { /*lint --e{715}*/
1901  assert(result != NULL);
1902  assert(argvals != NULL);
1903 
1904  if( argvals[0] > 0 )
1905  *result = pow( argvals[0], opdata.dbl);
1906  else
1907  *result = -pow(-argvals[0], opdata.dbl);
1908 
1909  return SCIP_OKAY;
1910 }
1911 
1912 /** interval evaluation for EXPR_SIGNPOWER */
1913 static
1914 SCIP_DECL_EXPRINTEVAL( exprevalIntSignPower )
1915 { /*lint --e{715}*/
1916  assert(result != NULL);
1917  assert(argvals != NULL);
1918 
1919  SCIPintervalSignPowerScalar(infinity, result, argvals[0], opdata.dbl);
1920 
1921  return SCIP_OKAY;
1922 }
1923 
1924 /** curvature for EXPR_SIGNPOWER */
1925 static
1926 SCIP_DECL_EXPRCURV( exprcurvSignPower )
1927 { /*lint --e{715}*/
1928  SCIP_INTERVAL tmp;
1929  SCIP_EXPRCURV left;
1930  SCIP_EXPRCURV right;
1931 
1932  assert(result != NULL);
1933  assert(argcurv != NULL);
1934  assert(argbounds != NULL);
1935 
1936  /* for x <= 0, signpower(x,c) = -(-x)^c
1937  * for x >= 0, signpower(x,c) = ( x)^c
1938  *
1939  * thus, get curvatures for both parts and "intersect" them
1940  */
1941 
1942  if( argbounds[0].inf < 0 )
1943  {
1944  SCIPintervalSetBounds(&tmp, 0.0, -argbounds[0].inf);
1945  left = SCIPexprcurvNegate(SCIPexprcurvPower(tmp, SCIPexprcurvNegate(argcurv[0]), opdata.dbl));
1946  }
1947  else
1948  {
1949  left = SCIP_EXPRCURV_LINEAR;
1950  }
1951 
1952  if( argbounds[0].sup > 0 )
1953  {
1954  SCIPintervalSetBounds(&tmp, 0.0, argbounds[0].sup);
1955  right = SCIPexprcurvPower(tmp, argcurv[0], opdata.dbl);
1956  }
1957  else
1958  {
1959  right = SCIP_EXPRCURV_LINEAR;
1960  }
1961 
1962  *result = (SCIP_EXPRCURV) (left & right);
1963 
1964  return SCIP_OKAY;
1965 }
1966 
1967 /** point evaluation for EXPR_EXP */
1968 static
1969 SCIP_DECL_EXPREVAL( exprevalExp )
1970 { /*lint --e{715}*/
1971  assert(result != NULL);
1972  assert(argvals != NULL);
1973 
1974  *result = exp(argvals[0]);
1975 
1976  return SCIP_OKAY;
1977 }
1978 
1979 /** interval evaluation for EXPR_EXP */
1980 static
1981 SCIP_DECL_EXPRINTEVAL( exprevalIntExp )
1982 { /*lint --e{715}*/
1983  assert(result != NULL);
1984  assert(argvals != NULL);
1985 
1986  SCIPintervalExp(infinity, result, argvals[0]);
1987 
1988  return SCIP_OKAY;
1989 }
1990 
1991 /** curvature for EXPR_EXP */
1992 static
1993 SCIP_DECL_EXPRCURV( exprcurvExp )
1994 { /*lint --e{715}*/
1995  assert(result != NULL);
1996  assert(argcurv != NULL);
1997 
1998  /* expression is convex if child is convex
1999  * otherwise, we don't know
2000  */
2001  if( argcurv[0] & SCIP_EXPRCURV_CONVEX )
2002  *result = SCIP_EXPRCURV_CONVEX;
2003  else
2004  *result = SCIP_EXPRCURV_UNKNOWN;
2005 
2006  return SCIP_OKAY;
2007 }
2008 
2009 /** point evaluation for EXPR_LOG */
2010 static
2011 SCIP_DECL_EXPREVAL( exprevalLog )
2012 { /*lint --e{715}*/
2013  assert(result != NULL);
2014  assert(argvals != NULL);
2015 
2016  *result = log(argvals[0]);
2017 
2018  return SCIP_OKAY;
2019 }
2020 
2021 /** interval evaluation for EXPR_LOG */
2022 static
2023 SCIP_DECL_EXPRINTEVAL( exprevalIntLog )
2024 { /*lint --e{715}*/
2025  assert(result != NULL);
2026  assert(argvals != NULL);
2027 
2028  SCIPintervalLog(infinity, result, argvals[0]);
2029 
2030  return SCIP_OKAY;
2031 }
2032 
2033 /** curvature for EXPR_LOG */
2034 static
2035 SCIP_DECL_EXPRCURV( exprcurvLog )
2036 { /*lint --e{715}*/
2037  assert(result != NULL);
2038  assert(argcurv != NULL);
2039 
2040  /* expression is concave if child is concave
2041  * otherwise, we don't know
2042  */
2043  if( argcurv[0] & SCIP_EXPRCURV_CONCAVE )
2044  *result = SCIP_EXPRCURV_CONCAVE;
2045  else
2046  *result = SCIP_EXPRCURV_UNKNOWN;
2047 
2048  return SCIP_OKAY;
2049 }
2050 
2051 /** point evaluation for EXPR_SIN */
2052 static
2053 SCIP_DECL_EXPREVAL( exprevalSin )
2054 { /*lint --e{715}*/
2055  assert(result != NULL);
2056  assert(argvals != NULL);
2057 
2058  *result = sin(argvals[0]);
2059 
2060  return SCIP_OKAY;
2061 }
2062 
2063 /** interval evaluation for EXPR_SIN */
2064 static
2065 SCIP_DECL_EXPRINTEVAL( exprevalIntSin )
2066 { /*lint --e{715}*/
2067  assert(result != NULL);
2068  assert(argvals != NULL);
2069  assert(nargs == 1);
2070 
2071  SCIPintervalSin(infinity, result, *argvals);
2072 
2073  return SCIP_OKAY;
2074 }
2075 
2076 /* @todo implement exprcurvSin */
2077 #define exprcurvSin exprcurvDefault
2078 
2079 /** point evaluation for EXPR_COS */
2080 static
2081 SCIP_DECL_EXPREVAL( exprevalCos )
2082 { /*lint --e{715}*/
2083  assert(result != NULL);
2084  assert(argvals != NULL);
2085 
2086  *result = cos(argvals[0]);
2087 
2088  return SCIP_OKAY;
2089 }
2090 
2091 /** interval evaluation for EXPR_COS */
2092 static
2093 SCIP_DECL_EXPRINTEVAL( exprevalIntCos )
2094 { /*lint --e{715}*/
2095  assert(result != NULL);
2096  assert(argvals != NULL);
2097  assert(nargs == 1);
2098 
2099  SCIPintervalCos(infinity, result, *argvals);
2100 
2101  return SCIP_OKAY;
2102 }
2103 
2104 /* @todo implement exprcurvCos */
2105 #define exprcurvCos exprcurvDefault
2106 
2107 /** point evaluation for EXPR_TAN */
2108 static
2109 SCIP_DECL_EXPREVAL( exprevalTan )
2110 { /*lint --e{715}*/
2111  assert(result != NULL);
2112  assert(argvals != NULL);
2113 
2114  *result = tan(argvals[0]);
2115 
2116  return SCIP_OKAY;
2117 }
2118 
2119 /* @todo implement SCIPintervalTan */
2120 #define exprevalIntTan exprevalIntDefault
2121 
2122 /* @todo implement exprcurvTan */
2123 #define exprcurvTan exprcurvDefault
2124 
2125 /* erf and erfi do not seem to exist on every system, and we cannot really handle them anyway, so they are currently disabled */
2126 #ifdef SCIP_DISABLED_CODE
2127 static
2128 SCIP_DECL_EXPREVAL( exprevalErf )
2129 { /*lint --e{715}*/
2130  assert(result != NULL);
2131  assert(argvals != NULL);
2132 
2133  *result = erf(argvals[0]);
2134 
2135  return SCIP_OKAY;
2136 }
2137 
2138 /* @todo implement SCIPintervalErf */
2139 #define exprevalIntErf exprevalIntDefault
2140 
2141 /* @todo implement SCIPintervalErf */
2142 #define exprcurvErf exprcurvDefault
2143 
2144 static
2145 SCIP_DECL_EXPREVAL( exprevalErfi )
2146 { /*lint --e{715}*/
2147  assert(result != NULL);
2148  assert(argvals != NULL);
2149 
2150  /* @TODO implement erfi evaluation */
2151  SCIPerrorMessage("erfi not implemented");
2152 
2153  return SCIP_ERROR;
2154 }
2155 
2156 /* @todo implement SCIPintervalErfi */
2157 #define exprevalIntErfi NULL
2158 
2159 #define exprcurvErfi exprcurvDefault
2160 #endif
2161 
2162 /** point evaluation for EXPR_MIN */
2163 static
2164 SCIP_DECL_EXPREVAL( exprevalMin )
2165 { /*lint --e{715}*/
2166  assert(result != NULL);
2167  assert(argvals != NULL);
2168 
2169  *result = MIN(argvals[0], argvals[1]);
2170 
2171  return SCIP_OKAY;
2172 }
2173 
2174 /** interval evaluation for EXPR_MIN */
2175 static
2176 SCIP_DECL_EXPRINTEVAL( exprevalIntMin )
2177 { /*lint --e{715}*/
2178  assert(result != NULL);
2179  assert(argvals != NULL);
2180 
2181  SCIPintervalMin(infinity, result, argvals[0], argvals[1]);
2182 
2183  return SCIP_OKAY;
2184 }
2185 
2186 /** curvature for EXPR_MIN */
2187 static
2188 SCIP_DECL_EXPRCURV( exprcurvMin )
2189 { /*lint --e{715}*/
2190  assert(result != NULL);
2191  assert(argcurv != NULL);
2192 
2193  /* the minimum of two concave functions is concave
2194  * otherwise, we don't know
2195  */
2196 
2197  if( (argcurv[0] & SCIP_EXPRCURV_CONCAVE) && (argcurv[1] & SCIP_EXPRCURV_CONCAVE) )
2198  *result = SCIP_EXPRCURV_CONCAVE;
2199  else
2200  *result = SCIP_EXPRCURV_UNKNOWN;
2201 
2202  return SCIP_OKAY;
2203 }
2204 
2205 /** point evaluation for EXPR_MAX */
2206 static
2207 SCIP_DECL_EXPREVAL( exprevalMax )
2208 { /*lint --e{715}*/
2209  assert(result != NULL);
2210  assert(argvals != NULL);
2211 
2212  *result = MAX(argvals[0], argvals[1]);
2213 
2214  return SCIP_OKAY;
2215 }
2216 
2217 /** interval evaluation for EXPR_MAX */
2218 static
2219 SCIP_DECL_EXPRINTEVAL( exprevalIntMax )
2220 { /*lint --e{715}*/
2221  assert(result != NULL);
2222  assert(argvals != NULL);
2223 
2224  SCIPintervalMax(infinity, result, argvals[0], argvals[1]);
2225 
2226  return SCIP_OKAY;
2227 }
2228 
2229 /** curvature for EXPR_MAX */
2230 static
2231 SCIP_DECL_EXPRCURV( exprcurvMax )
2232 { /*lint --e{715}*/
2233  assert(result != NULL);
2234  assert(argcurv != NULL);
2235 
2236  /* the maximum of two convex functions is convex
2237  * otherwise, we don't know
2238  */
2239  if( (argcurv[0] & SCIP_EXPRCURV_CONVEX) && (argcurv[1] & SCIP_EXPRCURV_CONVEX) )
2240  *result = SCIP_EXPRCURV_CONVEX;
2241  else
2242  *result = SCIP_EXPRCURV_UNKNOWN;
2243 
2244  return SCIP_OKAY;
2245 }
2246 
2247 /** point evaluation for EXPR_ABS */
2248 static
2249 SCIP_DECL_EXPREVAL( exprevalAbs )
2250 { /*lint --e{715}*/
2251  assert(result != NULL);
2252  assert(argvals != NULL);
2253 
2254  *result = ABS(argvals[0]);
2255 
2256  return SCIP_OKAY;
2257 }
2258 
2259 /** interval evaluation for EXPR_ABS */
2260 static
2261 SCIP_DECL_EXPRINTEVAL( exprevalIntAbs )
2262 { /*lint --e{715}*/
2263  assert(result != NULL);
2264  assert(argvals != NULL);
2265 
2266  SCIPintervalAbs(infinity, result, argvals[0]);
2267 
2268  return SCIP_OKAY;
2269 }
2270 
2271 /** curvature for EXPR_ABS */
2272 static
2273 SCIP_DECL_EXPRCURV( exprcurvAbs )
2274 { /*lint --e{715}*/
2275  assert(result != NULL);
2276  assert(argcurv != NULL);
2277  assert(argbounds != NULL);
2278 
2279  /* if child is only negative, then abs(child) = -child
2280  * if child is only positive, then abs(child) = child
2281  * if child is both positive and negative, but also linear, then abs(child) is convex
2282  * otherwise, we don't know
2283  */
2284  if( argbounds[0].sup <= 0.0 )
2285  *result = SCIPexprcurvMultiply(-1.0, argcurv[0]);
2286  else if( argbounds[0].inf >= 0.0 )
2287  *result = argcurv[0];
2288  else if( argcurv[0] == SCIP_EXPRCURV_LINEAR )
2289  *result = SCIP_EXPRCURV_CONVEX;
2290  else
2291  *result = SCIP_EXPRCURV_UNKNOWN;
2292 
2293  return SCIP_OKAY;
2294 }
2295 
2296 /** point evaluation for EXPR_SIGN */
2297 static
2298 SCIP_DECL_EXPREVAL( exprevalSign )
2299 { /*lint --e{715}*/
2300  assert(result != NULL);
2301  assert(argvals != NULL);
2302 
2303  *result = SIGN(argvals[0]);
2304 
2305  return SCIP_OKAY;
2306 }
2307 
2308 /** interval evaluation for EXPR_SIGN */
2309 static
2310 SCIP_DECL_EXPRINTEVAL( exprevalIntSign )
2311 { /*lint --e{715}*/
2312  assert(result != NULL);
2313  assert(argvals != NULL);
2314 
2315  SCIPintervalSign(infinity, result, argvals[0]);
2316 
2317  return SCIP_OKAY;
2318 }
2319 
2320 /** curvature for EXPR_SIGN */
2321 static
2322 SCIP_DECL_EXPRCURV( exprcurvSign )
2323 { /*lint --e{715}*/
2324  assert(result != NULL);
2325  assert(argbounds != NULL);
2326 
2327  /* if sign of child is clear, then sign is linear otherwise, we don't know */
2328  if( argbounds[0].sup <= 0.0 || argbounds[0].inf >= 0.0 )
2329  *result = SCIP_EXPRCURV_LINEAR;
2330  else
2331  *result = SCIP_EXPRCURV_UNKNOWN;
2332 
2333  return SCIP_OKAY;
2334 }
2335 
2336 /** point evaluation for EXPR_SUM */
2337 static
2338 SCIP_DECL_EXPREVAL( exprevalSum )
2339 { /*lint --e{715}*/
2340  int i;
2341 
2342  assert(result != NULL);
2343  assert(argvals != NULL);
2344 
2345  *result = 0.0;
2346  for( i = 0; i < nargs; ++i )
2347  *result += argvals[i];
2348 
2349  return SCIP_OKAY;
2350 }
2351 
2352 /** interval evaluation for EXPR_SUM */
2353 static
2354 SCIP_DECL_EXPRINTEVAL( exprevalIntSum )
2355 { /*lint --e{715}*/
2356  int i;
2357 
2358  assert(result != NULL);
2359  assert(argvals != NULL);
2360 
2361  SCIPintervalSet(result, 0.0);
2362 
2363  for( i = 0; i < nargs; ++i )
2364  SCIPintervalAdd(infinity, result, *result, argvals[i]);
2365 
2366  return SCIP_OKAY;
2367 }
2368 
2369 /** curvature for EXPR_SUM */
2370 static
2371 SCIP_DECL_EXPRCURV( exprcurvSum )
2372 { /*lint --e{715}*/
2373  int i;
2374 
2375  assert(result != NULL);
2376  assert(argcurv != NULL);
2377 
2378  *result = SCIP_EXPRCURV_LINEAR;
2379 
2380  for( i = 0; i < nargs; ++i )
2381  *result = SCIPexprcurvAdd(*result, argcurv[i]);
2382 
2383  return SCIP_OKAY;
2384 }
2385 
2386 /** point evaluation for EXPR_PRODUCT */
2387 static
2388 SCIP_DECL_EXPREVAL( exprevalProduct )
2389 { /*lint --e{715}*/
2390  int i;
2391 
2392  assert(result != NULL);
2393  assert(argvals != NULL);
2394 
2395  *result = 1.0;
2396  for( i = 0; i < nargs; ++i )
2397  *result *= argvals[i];
2398 
2399  return SCIP_OKAY;
2400 }
2401 
2402 /** interval evaluation for EXPR_PRODUCT */
2403 static
2404 SCIP_DECL_EXPRINTEVAL( exprevalIntProduct )
2405 { /*lint --e{715}*/
2406  int i;
2407 
2408  assert(result != NULL);
2409  assert(argvals != NULL);
2410 
2411  SCIPintervalSet(result, 1.0);
2412 
2413  for( i = 0; i < nargs; ++i )
2414  SCIPintervalMul(infinity, result, *result, argvals[i]);
2415 
2416  return SCIP_OKAY;
2417 }
2418 
2419 /** curvature for EXPR_PRODUCT */
2420 static
2421 SCIP_DECL_EXPRCURV( exprcurvProduct )
2422 { /*lint --e{715}*/
2423  SCIP_Bool hadnonconst;
2424  SCIP_Real constants;
2425  int i;
2426 
2427  assert(result != NULL);
2428  assert(argcurv != NULL);
2429  assert(argbounds != NULL);
2430 
2431  /* if all factors are constant, then product is linear (even constant)
2432  * if only one factor is not constant, then product is curvature of this factor, multiplied by sign of product of remaining factors
2433  */
2434  *result = SCIP_EXPRCURV_LINEAR;
2435  hadnonconst = FALSE;
2436  constants = 1.0;
2437 
2438  for( i = 0; i < nargs; ++i )
2439  {
2440  if( argbounds[i].inf == argbounds[i].sup ) /*lint !e777*/
2441  {
2442  constants *= argbounds[i].inf;
2443  }
2444  else if( !hadnonconst )
2445  {
2446  /* first non-constant child */
2447  *result = argcurv[i];
2448  hadnonconst = TRUE;
2449  }
2450  else
2451  {
2452  /* more than one non-constant child, thus don't know curvature */
2453  *result = SCIP_EXPRCURV_UNKNOWN;
2454  break;
2455  }
2456  }
2457 
2458  *result = SCIPexprcurvMultiply(constants, *result);
2459 
2460  return SCIP_OKAY;
2461 }
2462 
2463 /** point evaluation for EXPR_LINEAR */
2464 static
2465 SCIP_DECL_EXPREVAL( exprevalLinear )
2466 { /*lint --e{715}*/
2467  SCIP_Real* coef;
2468  int i;
2469 
2470  assert(result != NULL);
2471  assert(argvals != NULL || nargs == 0);
2472  assert(opdata.data != NULL);
2473 
2474  coef = &((SCIP_Real*)opdata.data)[nargs];
2475 
2476  *result = *coef;
2477  for( i = nargs-1, --coef; i >= 0; --i, --coef )
2478  *result += *coef * argvals[i]; /*lint !e613*/
2479 
2480  assert(++coef == (SCIP_Real*)opdata.data);
2481 
2482  return SCIP_OKAY;
2483 }
2484 
2485 /** interval evaluation for EXPR_LINEAR */
2486 static
2487 SCIP_DECL_EXPRINTEVAL( exprevalIntLinear )
2488 { /*lint --e{715}*/
2489  assert(result != NULL);
2490  assert(argvals != NULL || nargs == 0);
2491  assert(opdata.data != NULL);
2492 
2493  SCIPintervalScalprodScalars(infinity, result, nargs, argvals, (SCIP_Real*)opdata.data);
2494  SCIPintervalAddScalar(infinity, result, *result, ((SCIP_Real*)opdata.data)[nargs]);
2495 
2496  return SCIP_OKAY;
2497 }
2498 
2499 /** curvature for EXPR_LINEAR */
2500 static
2501 SCIP_DECL_EXPRCURV( exprcurvLinear )
2502 { /*lint --e{715}*/
2503  SCIP_Real* data;
2504  int i;
2505 
2506  assert(result != NULL);
2507  assert(argcurv != NULL);
2508 
2509  data = (SCIP_Real*)opdata.data;
2510  assert(data != NULL);
2511 
2512  *result = SCIP_EXPRCURV_LINEAR;
2513 
2514  for( i = 0; i < nargs; ++i )
2515  *result = SCIPexprcurvAdd(*result, SCIPexprcurvMultiply(data[i], argcurv[i]));
2516 
2517  return SCIP_OKAY;
2518 }
2519 
2520 /** expression data copy for EXPR_LINEAR */
2521 static
2522 SCIP_DECL_EXPRCOPYDATA( exprCopyDataLinear )
2523 { /*lint --e{715}*/
2524  SCIP_Real* targetdata;
2525 
2526  assert(blkmem != NULL);
2527  assert(nchildren >= 0);
2528  assert(opdatatarget != NULL);
2529 
2530  /* for a linear expression, we need to copy the array that holds the coefficients and constant term */
2531  assert(opdatasource.data != NULL);
2532  SCIP_ALLOC( BMSduplicateBlockMemoryArray(blkmem, &targetdata, (SCIP_Real*)opdatasource.data, nchildren + 1) ); /*lint !e866*/
2533  opdatatarget->data = targetdata;
2534 
2535  return SCIP_OKAY;
2536 }
2537 
2538 /** expression data free for EXPR_LINEAR */
2539 static
2540 SCIP_DECL_EXPRFREEDATA( exprFreeDataLinear )
2541 { /*lint --e{715}*/
2542  SCIP_Real* freedata;
2543 
2544  assert(blkmem != NULL);
2545  assert(nchildren >= 0);
2546 
2547  freedata = (SCIP_Real*)opdata.data;
2548  assert(freedata != NULL);
2549 
2550  BMSfreeBlockMemoryArray(blkmem, &freedata, nchildren + 1); /*lint !e866*/
2551 }
2552 
2553 /** point evaluation for EXPR_QUADRATIC */
2554 static
2555 SCIP_DECL_EXPREVAL( exprevalQuadratic )
2556 { /*lint --e{715}*/
2557  SCIP_EXPRDATA_QUADRATIC* quaddata;
2558  SCIP_Real* lincoefs;
2559  SCIP_QUADELEM* quadelems;
2560  int nquadelems;
2561  int i;
2562 
2563  assert(result != NULL);
2564  assert(argvals != NULL || nargs == 0);
2565 
2566  quaddata = (SCIP_EXPRDATA_QUADRATIC*)opdata.data;
2567  assert(quaddata != NULL);
2568 
2569  lincoefs = quaddata->lincoefs;
2570  nquadelems = quaddata->nquadelems;
2571  quadelems = quaddata->quadelems;
2572 
2573  assert(quadelems != NULL || nquadelems == 0);
2574  assert(argvals != NULL || nquadelems == 0);
2575 
2576  *result = quaddata->constant;
2577 
2578  if( lincoefs != NULL )
2579  {
2580  for( i = nargs-1; i >= 0; --i )
2581  *result += lincoefs[i] * argvals[i]; /*lint !e613*/
2582  }
2583 
2584  for( i = 0; i < nquadelems; ++i, ++quadelems ) /*lint !e613*/
2585  *result += quadelems->coef * argvals[quadelems->idx1] * argvals[quadelems->idx2]; /*lint !e613*/
2586 
2587  return SCIP_OKAY;
2588 }
2589 
2590 /** interval evaluation for EXPR_QUADRATIC */
2591 static
2592 SCIP_DECL_EXPRINTEVAL( exprevalIntQuadratic )
2593 { /*lint --e{715}*/
2594  SCIP_EXPRDATA_QUADRATIC* quaddata;
2595  SCIP_Real* lincoefs;
2596  SCIP_QUADELEM* quadelems;
2597  int nquadelems;
2598  int i;
2599  int argidx;
2600  SCIP_Real sqrcoef;
2601  SCIP_INTERVAL lincoef;
2602  SCIP_INTERVAL tmp;
2603 
2604  assert(result != NULL);
2605  assert(argvals != NULL || nargs == 0);
2606 
2607  quaddata = (SCIP_EXPRDATA_QUADRATIC*)opdata.data;
2608  assert(quaddata != NULL);
2609 
2610  lincoefs = quaddata->lincoefs;
2611  nquadelems = quaddata->nquadelems;
2612  quadelems = quaddata->quadelems;
2613 
2614  assert(quadelems != NULL || nquadelems == 0);
2615  assert(argvals != NULL || nargs == 0);
2616 
2617  /* something fast for case of only one child */
2618  if( nargs == 1 )
2619  {
2620  SCIPintervalSet(&lincoef, lincoefs != NULL ? lincoefs[0] : 0.0);
2621 
2622  sqrcoef = 0.0;
2623  for( i = 0; i < nquadelems; ++i )
2624  {
2625  assert(quadelems[i].idx1 == 0); /*lint !e613*/
2626  assert(quadelems[i].idx2 == 0); /*lint !e613*/
2627  sqrcoef += quadelems[i].coef; /*lint !e613*/
2628  }
2629 
2630  SCIPintervalQuad(infinity, result, sqrcoef, lincoef, argvals[0]); /*lint !e613*/
2631  SCIPintervalAddScalar(infinity, result, *result, quaddata->constant);
2632 
2633  return SCIP_OKAY;
2634  }
2635 
2636  if( nargs == 2 && nquadelems > 0 )
2637  {
2638  /* if it's a bivariate quadratic expression with bilinear term, do something special */
2639  SCIP_Real ax; /* square coefficient of first child */
2640  SCIP_Real ay; /* square coefficient of second child */
2641  SCIP_Real axy; /* bilinear coefficient */
2642 
2643  ax = 0.0;
2644  ay = 0.0;
2645  axy = 0.0;
2646  for( i = 0; i < nquadelems; ++i )
2647  if( quadelems[i].idx1 == 0 && quadelems[i].idx2 == 0 ) /*lint !e613*/
2648  ax += quadelems[i].coef; /*lint !e613*/
2649  else if( quadelems[i].idx1 == 1 && quadelems[i].idx2 == 1 ) /*lint !e613*/
2650  ay += quadelems[i].coef; /*lint !e613*/
2651  else
2652  axy += quadelems[i].coef; /*lint !e613*/
2653 
2654  SCIPintervalQuadBivar(infinity, result, ax, ay, axy,
2655  lincoefs != NULL ? lincoefs[0] : 0.0, lincoefs != NULL ? lincoefs[1] : 0.0,
2656  argvals[0], argvals[1]); /*lint !e613*/
2657  SCIPdebugMessage("%g x^2 + %g y^2 + %g x y + %g x + %g y = [%g,%g] for x = [%g,%g], y = [%g,%g]\n",
2658  ax, ay, axy, lincoefs != NULL ? lincoefs[0] : 0.0, lincoefs != NULL ? lincoefs[1] : 0.0,
2659  result->inf, result->sup, argvals[0].inf, argvals[0].sup, argvals[1].inf, argvals[1].sup); /*lint !e613*/
2660 
2661  SCIPintervalAddScalar(infinity, result, *result, quaddata->constant);
2662 
2663  return SCIP_OKAY;
2664  }
2665 
2666  /* make sure coefficients are sorted */
2667  quadraticdataSort(quaddata);
2668 
2669  SCIPintervalSet(result, quaddata->constant);
2670 
2671  /* for each argument, we collect it's linear index from lincoefs, it's square coefficients and all factors from bilinear terms
2672  * then we compute the interval sqrcoef*x^2 + lincoef*x and add it to result
2673  * @todo split quadratic expression into bivariate quadratic terms and apply the above method
2674  */
2675  i = 0;
2676  for( argidx = 0; argidx < nargs; ++argidx )
2677  {
2678  if( i == nquadelems || quadelems[i].idx1 > argidx ) /*lint !e613*/
2679  {
2680  /* there are no quadratic terms with argidx in its first argument, that should be easy to handle */
2681  if( lincoefs != NULL )
2682  {
2683  SCIPintervalMulScalar(infinity, &tmp, argvals[argidx], lincoefs[argidx]); /*lint !e613*/
2684  SCIPintervalAdd(infinity, result, *result, tmp);
2685  }
2686  continue;
2687  }
2688 
2689  sqrcoef = 0.0;
2690  SCIPintervalSet(&lincoef, lincoefs != NULL ? lincoefs[argidx] : 0.0);
2691 
2692  assert(i < nquadelems && quadelems[i].idx1 == argidx); /*lint !e613*/
2693  do
2694  {
2695  if( quadelems[i].idx2 == argidx ) /*lint !e613*/
2696  {
2697  sqrcoef += quadelems[i].coef; /*lint !e613*/
2698  }
2699  else
2700  {
2701  SCIPintervalMulScalar(infinity, &tmp, argvals[quadelems[i].idx2], quadelems[i].coef); /*lint !e613*/
2702  SCIPintervalAdd(infinity, &lincoef, lincoef, tmp);
2703  }
2704  ++i;
2705  }
2706  while( i < nquadelems && quadelems[i].idx1 == argidx ); /*lint !e613*/
2707  assert(i == nquadelems || quadelems[i].idx1 > argidx); /*lint !e613*/
2708 
2709  SCIPintervalQuad(infinity, &tmp, sqrcoef, lincoef, argvals[argidx]); /*lint !e613*/
2710  SCIPintervalAdd(infinity, result, *result, tmp);
2711  }
2712  assert(i == nquadelems);
2713 
2714  return SCIP_OKAY;
2715 }
2716 
2717 /** curvature for EXPR_QUADRATIC */
2718 static
2719 SCIP_DECL_EXPRCURV( exprcurvQuadratic )
2720 { /*lint --e{715}*/
2722  SCIP_QUADELEM* quadelems;
2723  int nquadelems;
2724  SCIP_Real* lincoefs;
2725  int i;
2726 
2727  assert(result != NULL);
2728  assert(argcurv != NULL);
2729  assert(argbounds != NULL);
2730 
2731  data = (SCIP_EXPRDATA_QUADRATIC*)opdata.data;
2732  assert(data != NULL);
2733 
2734  lincoefs = data->lincoefs;
2735  quadelems = data->quadelems;
2736  nquadelems = data->nquadelems;
2737 
2738  *result = SCIP_EXPRCURV_LINEAR;
2739 
2740  if( lincoefs != NULL )
2741  for( i = 0; i < nargs; ++i )
2742  *result = SCIPexprcurvAdd(*result, SCIPexprcurvMultiply(lincoefs[i], argcurv[i]));
2743 
2744  /* @todo could try cholesky factorization if all children linear...
2745  * @todo should then cache the result
2746  */
2747  for( i = 0; i < nquadelems && *result != SCIP_EXPRCURV_UNKNOWN; ++i )
2748  {
2749  if( quadelems[i].coef == 0.0 )
2750  continue;
2751 
2752  if( argbounds[quadelems[i].idx1].inf == argbounds[quadelems[i].idx1].sup && /*lint !e777*/
2753  +argbounds[quadelems[i].idx2].inf == argbounds[quadelems[i].idx2].sup
2754  ) /*lint !e777*/
2755  {
2756  /* both factors are constants -> curvature does not change */
2757  continue;
2758  }
2759 
2760  if( argbounds[quadelems[i].idx1].inf == argbounds[quadelems[i].idx1].sup ) /*lint !e777*/
2761  {
2762  /* first factor is constant, second is not -> add curvature of second */
2763  *result = SCIPexprcurvAdd(*result, SCIPexprcurvMultiply(quadelems[i].coef * argbounds[quadelems[i].idx1].inf, argcurv[quadelems[i].idx2]));
2764  }
2765  else if( argbounds[quadelems[i].idx2].inf == argbounds[quadelems[i].idx2].sup ) /*lint !e777*/
2766  {
2767  /* first factor is not constant, second is -> add curvature of first */
2768  *result = SCIPexprcurvAdd(*result, SCIPexprcurvMultiply(quadelems[i].coef * argbounds[quadelems[i].idx2].inf, argcurv[quadelems[i].idx1]));
2769  }
2770  else if( quadelems[i].idx1 == quadelems[i].idx2 )
2771  {
2772  /* both factors not constant, but the same (square term) */
2773  *result = SCIPexprcurvAdd(*result, SCIPexprcurvMultiply(quadelems[i].coef, SCIPexprcurvPower(argbounds[quadelems[i].idx1], argcurv[quadelems[i].idx1], 2.0)));
2774  }
2775  else
2776  {
2777  /* two different non-constant factors -> can't tell about curvature */
2778  *result = SCIP_EXPRCURV_UNKNOWN;
2779  }
2780  }
2781 
2782  return SCIP_OKAY;
2783 }
2784 
2785 /** expression data copy for EXPR_QUADRATIC */
2786 static
2787 SCIP_DECL_EXPRCOPYDATA( exprCopyDataQuadratic )
2788 { /*lint --e{715}*/
2789  SCIP_EXPRDATA_QUADRATIC* sourcedata;
2790 
2791  assert(blkmem != NULL);
2792  assert(opdatatarget != NULL);
2793 
2794  sourcedata = (SCIP_EXPRDATA_QUADRATIC*)opdatasource.data;
2795  assert(sourcedata != NULL);
2796 
2797  SCIP_CALL( quadraticdataCreate(blkmem, (SCIP_EXPRDATA_QUADRATIC**)&opdatatarget->data,
2798  sourcedata->constant, nchildren, sourcedata->lincoefs, sourcedata->nquadelems, sourcedata->quadelems) );
2799 
2800  return SCIP_OKAY;
2801 }
2802 
2803 /** expression data free for EXPR_QUADRATIC */
2804 static
2805 SCIP_DECL_EXPRFREEDATA( exprFreeDataQuadratic )
2806 { /*lint --e{715}*/
2807  SCIP_EXPRDATA_QUADRATIC* quadraticdata;
2808 
2809  assert(blkmem != NULL);
2810  assert(nchildren >= 0);
2811 
2812  quadraticdata = (SCIP_EXPRDATA_QUADRATIC*)opdata.data;
2813  assert(quadraticdata != NULL);
2814 
2815  if( quadraticdata->lincoefs != NULL )
2816  {
2817  BMSfreeBlockMemoryArray(blkmem, &quadraticdata->lincoefs, nchildren);
2818  }
2819 
2820  if( quadraticdata->nquadelems > 0 )
2821  {
2822  assert(quadraticdata->quadelems != NULL);
2823  BMSfreeBlockMemoryArray(blkmem, &quadraticdata->quadelems, quadraticdata->nquadelems);
2824  }
2825 
2826  BMSfreeBlockMemory(blkmem, &quadraticdata);
2827 }
2828 
2829 /** point evaluation for EXPR_POLYNOMIAL */
2830 static
2831 SCIP_DECL_EXPREVAL( exprevalPolynomial )
2832 { /*lint --e{715}*/
2833  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
2834  SCIP_EXPRDATA_MONOMIAL* monomialdata;
2835  SCIP_Real childval;
2836  SCIP_Real exponent;
2837  SCIP_Real monomialval;
2838  int i;
2839  int j;
2840 
2841  assert(result != NULL);
2842  assert(argvals != NULL || nargs == 0);
2843  assert(opdata.data != NULL);
2844 
2845  polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)opdata.data;
2846  assert(polynomialdata != NULL);
2847 
2848  *result = polynomialdata->constant;
2849 
2850  for( i = 0; i < polynomialdata->nmonomials; ++i )
2851  {
2852  monomialdata = polynomialdata->monomials[i];
2853  assert(monomialdata != NULL);
2854 
2855  monomialval = monomialdata->coef;
2856  for( j = 0; j < monomialdata->nfactors; ++j )
2857  {
2858  assert(monomialdata->childidxs[j] >= 0);
2859  assert(monomialdata->childidxs[j] < nargs);
2860 
2861  childval = argvals[monomialdata->childidxs[j]]; /*lint !e613*/
2862  if( childval == 1.0 ) /* 1^anything == 1 */
2863  continue;
2864 
2865  exponent = monomialdata->exponents[j];
2866 
2867  if( childval == 0.0 )
2868  {
2869  if( exponent > 0.0 )
2870  {
2871  /* 0^positive == 0 */
2872  monomialval = 0.0;
2873  break;
2874  }
2875  else if( exponent < 0.0 )
2876  {
2877  /* 0^negative = nan (or should it be +inf?, doesn't really matter) */
2878 #ifdef NAN
2879  *result = NAN;
2880 #else
2881  /* cppcheck-suppress wrongmathcall */
2882  *result = pow(0.0, -1.0);
2883 #endif
2884  return SCIP_OKAY;
2885  }
2886  /* 0^0 == 1 */
2887  continue;
2888  }
2889 
2890  /* cover some special exponents separately to avoid calling expensive pow function */
2891  if( exponent == 0.0 )
2892  continue;
2893  if( exponent == 1.0 )
2894  {
2895  monomialval *= childval;
2896  continue;
2897  }
2898  if( exponent == 2.0 )
2899  {
2900  monomialval *= childval * childval;
2901  continue;
2902  }
2903  if( exponent == 0.5 )
2904  {
2905  monomialval *= sqrt(childval);
2906  continue;
2907  }
2908  if( exponent == -1.0 )
2909  {
2910  monomialval /= childval;
2911  continue;
2912  }
2913  if( exponent == -2.0 )
2914  {
2915  monomialval /= childval * childval;
2916  continue;
2917  }
2918  monomialval *= pow(childval, exponent);
2919  }
2920 
2921  *result += monomialval;
2922  }
2923 
2924  return SCIP_OKAY;
2925 }
2926 
2927 /** interval evaluation for EXPR_POLYNOMIAL */
2928 static
2929 SCIP_DECL_EXPRINTEVAL( exprevalIntPolynomial )
2930 { /*lint --e{715}*/
2931  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
2932  SCIP_EXPRDATA_MONOMIAL* monomialdata;
2933  SCIP_INTERVAL childval;
2934  SCIP_INTERVAL monomialval;
2935  SCIP_Real exponent;
2936  int i;
2937  int j;
2938 
2939  assert(result != NULL);
2940  assert(argvals != NULL || nargs == 0);
2941  assert(opdata.data != NULL);
2942 
2943  polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)opdata.data;
2944  assert(polynomialdata != NULL);
2945 
2946  SCIPintervalSet(result, polynomialdata->constant);
2947 
2948  for( i = 0; i < polynomialdata->nmonomials; ++i )
2949  {
2950  monomialdata = polynomialdata->monomials[i];
2951  assert(monomialdata != NULL);
2952 
2953  SCIPintervalSet(&monomialval, monomialdata->coef);
2954  for( j = 0; j < monomialdata->nfactors && !SCIPintervalIsEntire(infinity, monomialval); ++j )
2955  {
2956  assert(monomialdata->childidxs[j] >= 0);
2957  assert(monomialdata->childidxs[j] < nargs);
2958 
2959  childval = argvals[monomialdata->childidxs[j]]; /*lint !e613*/
2960 
2961  exponent = monomialdata->exponents[j];
2962 
2963  /* cover some special exponents separately to avoid calling expensive pow function */
2964  if( exponent == 0.0 )
2965  continue;
2966 
2967  if( exponent == 1.0 )
2968  {
2969  SCIPintervalMul(infinity, &monomialval, monomialval, childval);
2970  continue;
2971  }
2972 
2973  if( exponent == 2.0 )
2974  {
2975  SCIPintervalSquare(infinity, &childval, childval);
2976  SCIPintervalMul(infinity, &monomialval, monomialval, childval);
2977  continue;
2978  }
2979 
2980  if( exponent == 0.5 )
2981  {
2982  SCIPintervalSquareRoot(infinity, &childval, childval);
2983  if( SCIPintervalIsEmpty(infinity, childval) )
2984  {
2985  SCIPintervalSetEmpty(result);
2986  break;
2987  }
2988  SCIPintervalMul(infinity, &monomialval, monomialval, childval);
2989  continue;
2990  }
2991  else if( exponent == -1.0 )
2992  {
2993  SCIPintervalDiv(infinity, &monomialval, monomialval, childval);
2994  }
2995  else if( exponent == -2.0 )
2996  {
2997  SCIPintervalSquare(infinity, &childval, childval);
2998  SCIPintervalDiv(infinity, &monomialval, monomialval, childval);
2999  }
3000  else
3001  {
3002  SCIPintervalPowerScalar(infinity, &childval, childval, exponent);
3003  if( SCIPintervalIsEmpty(infinity, childval) )
3004  {
3005  SCIPintervalSetEmpty(result);
3006  return SCIP_OKAY;
3007  }
3008  SCIPintervalMul(infinity, &monomialval, monomialval, childval);
3009  }
3010 
3011  /* the cases in which monomialval gets empty should have been catched */
3012  assert(!SCIPintervalIsEmpty(infinity, monomialval));
3013  }
3014 
3015  SCIPintervalAdd(infinity, result, *result, monomialval);
3016  }
3017 
3018  return SCIP_OKAY;
3019 }
3020 
3021 /** curvature for EXPR_POLYNOMIAL */
3022 static
3023 SCIP_DECL_EXPRCURV( exprcurvPolynomial )
3024 { /*lint --e{715}*/
3026  SCIP_EXPRDATA_MONOMIAL** monomials;
3027  SCIP_EXPRDATA_MONOMIAL* monomial;
3028  int nmonomials;
3029  int i;
3030 
3031  assert(result != NULL);
3032  assert(argcurv != NULL);
3033  assert(argbounds != NULL);
3034 
3035  data = (SCIP_EXPRDATA_POLYNOMIAL*)opdata.data;
3036  assert(data != NULL);
3037 
3038  monomials = data->monomials;
3039  nmonomials = data->nmonomials;
3040 
3041  *result = SCIP_EXPRCURV_LINEAR;
3042 
3043  for( i = 0; i < nmonomials && *result != SCIP_EXPRCURV_UNKNOWN; ++i )
3044  {
3045  /* we assume that some simplifier was running, so that monomials do not have constants in their factors and such that all factors are different
3046  * (result would still be correct)
3047  */
3048  monomial = monomials[i];
3049  *result = SCIPexprcurvAdd(*result, SCIPexprcurvMultiply(monomial->coef, SCIPexprcurvMonomial(monomial->nfactors, monomial->exponents, monomial->childidxs, argcurv, argbounds)));
3050  }
3051 
3052  return SCIP_OKAY;
3053 }
3054 
3055 /** expression data copy for EXPR_POLYNOMIAL */
3056 static
3057 SCIP_DECL_EXPRCOPYDATA( exprCopyDataPolynomial )
3058 { /*lint --e{715}*/
3059  SCIP_EXPRDATA_POLYNOMIAL* sourcepolynomialdata;
3060  SCIP_EXPRDATA_POLYNOMIAL* targetpolynomialdata;
3061 
3062  assert(blkmem != NULL);
3063  assert(opdatatarget != NULL);
3064 
3065  sourcepolynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)opdatasource.data;
3066  assert(sourcepolynomialdata != NULL);
3067 
3068  SCIP_CALL( polynomialdataCopy(blkmem, &targetpolynomialdata, sourcepolynomialdata) );
3069 
3070  opdatatarget->data = (void*)targetpolynomialdata;
3071 
3072  return SCIP_OKAY;
3073 }
3074 
3075 /** expression data free for EXPR_POLYNOMIAL */
3076 static
3077 SCIP_DECL_EXPRFREEDATA( exprFreeDataPolynomial )
3078 { /*lint --e{715}*/
3079  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3080 
3081  assert(blkmem != NULL);
3082 
3083  polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)opdata.data;
3084  assert(polynomialdata != NULL);
3085 
3086  polynomialdataFree(blkmem, &polynomialdata);
3087 }
3088 
3089 /** point evaluation for user expression */
3090 static
3091 SCIP_DECL_EXPREVAL( exprevalUser )
3092 { /*lint --e{715}*/
3093  SCIP_EXPRDATA_USER* exprdata;
3094 
3095  exprdata = (SCIP_EXPRDATA_USER*) opdata.data;
3096 
3097  SCIP_CALL( exprdata->eval(exprdata->userdata, nargs, argvals, result, NULL, NULL) );
3098 
3099  return SCIP_OKAY;
3100 }
3101 
3102 /** interval evaluation for user expression */
3103 static
3104 SCIP_DECL_EXPRINTEVAL( exprevalIntUser )
3105 { /*lint --e{715}*/
3106  SCIP_EXPRDATA_USER* exprdata;
3107 
3108  exprdata = (SCIP_EXPRDATA_USER*) opdata.data;
3109 
3110  if( exprdata->inteval != NULL )
3111  {
3112  SCIP_CALL( exprdata->inteval(infinity, exprdata->userdata, nargs, argvals, result, NULL, NULL) );
3113  }
3114  else
3115  {
3116  /* if user does not provide interval evaluation, then return a result that is always correct */
3118  }
3119 
3120  return SCIP_OKAY;
3121 }
3122 
3123 /** curvature check for user expression */
3124 static
3125 SCIP_DECL_EXPRCURV( exprcurvUser )
3126 {
3127  SCIP_EXPRDATA_USER* exprdata;
3128 
3129  exprdata = (SCIP_EXPRDATA_USER*) opdata.data;
3130 
3131  if( exprdata->curv != NULL )
3132  {
3133  SCIP_CALL( exprdata->curv(infinity, exprdata->userdata, nargs, argbounds, argcurv, result) );
3134  }
3135  else
3136  {
3137  /* if user does not provide curvature check, then return unknown (which is handled like indefinite) */
3138  *result = SCIP_EXPRCURV_UNKNOWN;
3139  }
3140 
3141  return SCIP_OKAY;
3142 }
3143 
3144 /** data copy for user expression */
3145 static
3146 SCIP_DECL_EXPRCOPYDATA( exprCopyDataUser )
3147 {
3148  SCIP_EXPRDATA_USER* exprdatasource;
3149  SCIP_EXPRDATA_USER* exprdatatarget;
3150 
3151  assert(blkmem != NULL);
3152  assert(opdatatarget != NULL);
3153 
3154  exprdatasource = (SCIP_EXPRDATA_USER*)opdatasource.data;
3155  assert(exprdatasource != NULL);
3156 
3157  /* duplicate expression data */
3158  SCIP_ALLOC( BMSduplicateBlockMemory(blkmem, &exprdatatarget, exprdatasource) );
3159 
3160  /* duplicate user expression data, if any */
3161  if( exprdatasource->copydata != NULL )
3162  {
3163  SCIP_CALL( exprdatasource->copydata(blkmem, nchildren, exprdatasource->userdata, &exprdatatarget->userdata) );
3164  }
3165  else
3166  {
3167  /* if no copy function for data, then there has to be no data */
3168  assert(exprdatatarget->userdata == NULL);
3169  }
3170 
3171  opdatatarget->data = (void*)exprdatatarget;
3172 
3173  return SCIP_OKAY;
3174 }
3175 
3176 /** data free for user expression */
3177 static
3178 SCIP_DECL_EXPRFREEDATA( exprFreeDataUser )
3179 {
3180  SCIP_EXPRDATA_USER* exprdata;
3181 
3182  assert(blkmem != NULL);
3183 
3184  exprdata = (SCIP_EXPRDATA_USER*)opdata.data;
3185 
3186  /* free user expression data, if any */
3187  if( exprdata->freedata != NULL )
3188  {
3189  exprdata->freedata(blkmem, nchildren, exprdata->userdata);
3190  }
3191  else
3192  {
3193  assert(exprdata->userdata == NULL);
3194  }
3195 
3196  /* free expression data */
3197  BMSfreeBlockMemory(blkmem, &exprdata);
3198 }
3199 
3200 /** element in table of expression operands */
3201 struct exprOpTableElement
3202 {
3203  const char* name; /**< name of operand (used for printing) */
3204  int nargs; /**< number of arguments (negative if not fixed) */
3205  SCIP_DECL_EXPREVAL ((*eval)); /**< evaluation function */
3206  SCIP_DECL_EXPRINTEVAL ((*inteval)); /**< interval evaluation function */
3207  SCIP_DECL_EXPRCURV ((*curv)); /**< curvature check function */
3208  SCIP_DECL_EXPRCOPYDATA ((*copydata)); /**< expression data copy function, or NULL to only opdata union */
3209  SCIP_DECL_EXPRFREEDATA ((*freedata)); /**< expression data free function, or NULL if nothing to free */
3210 };
3211 
3212 #define EXPROPEMPTY {NULL, -1, NULL, NULL, NULL, NULL, NULL}
3213 
3214 /** table containing for each operand the name, the number of children, and some evaluation functions */
3215 static
3216 struct exprOpTableElement exprOpTable[] =
3217  {
3218  EXPROPEMPTY,
3219  { "variable", 0, exprevalVar, exprevalIntVar, exprcurvVar, NULL, NULL },
3220  { "constant", 0, exprevalConst, exprevalIntConst, exprcurvConst, NULL, NULL },
3221  { "parameter", 0, exprevalParam, exprevalIntParam, exprcurvParam, NULL, NULL },
3223  { "plus", 2, exprevalPlus, exprevalIntPlus, exprcurvPlus, NULL, NULL },
3224  { "minus", 2, exprevalMinus, exprevalIntMinus, exprcurvMinus, NULL, NULL },
3225  { "mul", 2, exprevalMult, exprevalIntMult, exprcurvMult, NULL, NULL },
3226  { "div", 2, exprevalDiv, exprevalIntDiv, exprcurvDiv, NULL, NULL },
3227  { "sqr", 1, exprevalSquare, exprevalIntSquare, exprcurvSquare, NULL, NULL },
3228  { "sqrt", 1, exprevalSquareRoot, exprevalIntSquareRoot, exprcurvSquareRoot, NULL, NULL },
3229  { "realpower", 1, exprevalRealPower, exprevalIntRealPower, exprcurvRealPower, NULL, NULL },
3230  { "intpower", 1, exprevalIntPower, exprevalIntIntPower, exprcurvIntPower, NULL, NULL },
3231  { "signpower", 1, exprevalSignPower, exprevalIntSignPower, exprcurvSignPower, NULL, NULL },
3232  { "exp", 1, exprevalExp, exprevalIntExp, exprcurvExp, NULL, NULL },
3233  { "log", 1, exprevalLog, exprevalIntLog, exprcurvLog, NULL, NULL },
3234  { "sin", 1, exprevalSin, exprevalIntSin, exprcurvSin, NULL, NULL },
3235  { "cos", 1, exprevalCos, exprevalIntCos, exprcurvCos, NULL, NULL },
3236  { "tan", 1, exprevalTan, exprevalIntTan, exprcurvTan, NULL, NULL },
3237  /* { "erf", 1, exprevalErf, exprevalIntErf, exprcurvErf, NULL, NULL }, */
3238  /* { "erfi", 1, exprevalErfi, exprevalIntErfi exprcurvErfi, NULL, NULL }, */
3240  { "min", 2, exprevalMin, exprevalIntMin, exprcurvMin, NULL, NULL },
3241  { "max", 2, exprevalMax, exprevalIntMax, exprcurvMax, NULL, NULL },
3242  { "abs", 1, exprevalAbs, exprevalIntAbs, exprcurvAbs, NULL, NULL },
3243  { "sign", 1, exprevalSign, exprevalIntSign, exprcurvSign, NULL, NULL },
3249  { "sum", -2, exprevalSum, exprevalIntSum, exprcurvSum, NULL, NULL },
3250  { "prod", -2, exprevalProduct, exprevalIntProduct, exprcurvProduct, NULL, NULL },
3251  { "linear", -2, exprevalLinear, exprevalIntLinear, exprcurvLinear, exprCopyDataLinear, exprFreeDataLinear },
3252  { "quadratic", -2, exprevalQuadratic, exprevalIntQuadratic, exprcurvQuadratic, exprCopyDataQuadratic, exprFreeDataQuadratic },
3253  { "polynomial", -2, exprevalPolynomial, exprevalIntPolynomial, exprcurvPolynomial, exprCopyDataPolynomial, exprFreeDataPolynomial },
3254  { "user", -2, exprevalUser, exprevalIntUser, exprcurvUser, exprCopyDataUser, exprFreeDataUser }
3255  };
3256 
3257 /**@} */
3258 
3259 /**@name Expression operand methods */
3260 /**@{ */
3261 
3262 /** gives the name of an operand as string */
3263 const char* SCIPexpropGetName(
3264  SCIP_EXPROP op /**< expression operand */
3265  )
3266 {
3267  assert(op < SCIP_EXPR_LAST);
3268 
3269  return exprOpTable[op].name;
3270 }
3271 
3272 /** gives the number of children of a simple operand */
3274  SCIP_EXPROP op /**< expression operand */
3275  )
3276 {
3277  assert(op < SCIP_EXPR_LAST);
3278 
3279  return exprOpTable[op].nargs;
3280 }
3281 
3282 /**@} */
3283 
3284 /**@name Expressions private methods */
3285 /**@{ */
3286 
3287 /** creates an expression
3288  *
3289  * Note, that the expression is allocated but for the children only the pointer is copied.
3290  */
3291 static
3293  BMS_BLKMEM* blkmem, /**< block memory data structure */
3294  SCIP_EXPR** expr, /**< pointer to buffer for expression address */
3295  SCIP_EXPROP op, /**< operand of expression */
3296  int nchildren, /**< number of children */
3297  SCIP_EXPR** children, /**< children */
3298  SCIP_EXPROPDATA opdata /**< operand data */
3299  )
3300 {
3301  assert(blkmem != NULL);
3302  assert(expr != NULL);
3303  assert(children != NULL || nchildren == 0);
3304  assert(children == NULL || nchildren > 0);
3305 
3306  SCIP_ALLOC( BMSallocBlockMemory(blkmem, expr) );
3307 
3308  (*expr)->op = op;
3309  (*expr)->nchildren = nchildren;
3310  (*expr)->children = children;
3311  (*expr)->data = opdata;
3312 
3313  return SCIP_OKAY;
3314 }
3315 
3316 /** tries to convert a given (operator,operatordata) pair into a polynomial operator with corresponding data
3317  *
3318  * Does not do this for constants.
3319  * If conversion is not possible or operator is already polynomial, *op and *data are
3320  * left untouched.
3321  */
3322 static
3324  BMS_BLKMEM* blkmem, /**< block memory */
3325  SCIP_EXPROP* op, /**< pointer to expression operator */
3326  SCIP_EXPROPDATA* data, /**< pointer to expression data */
3327  int nchildren /**< number of children of operator */
3328  )
3329 {
3330  assert(blkmem != NULL);
3331  assert(op != NULL);
3332  assert(data != NULL);
3333 
3334  switch( *op )
3335  {
3336  case SCIP_EXPR_VARIDX:
3337  case SCIP_EXPR_PARAM:
3338  case SCIP_EXPR_CONST:
3339  break;
3340 
3341  case SCIP_EXPR_PLUS:
3342  {
3343  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3344  SCIP_EXPRDATA_MONOMIAL* monomials[2];
3345  int childidx;
3346  SCIP_Real exponent;
3347 
3348  assert(nchildren == 2);
3349 
3350  /* create monomial for first child */
3351  childidx = 0;
3352  exponent = 1.0;
3353  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomials[0], 1.0, 1, &childidx, &exponent) );
3354 
3355  /* create monomial for second child */
3356  childidx = 1;
3357  exponent = 1.0;
3358  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomials[1], 1.0, 1, &childidx, &exponent) );
3359 
3360  /* create polynomial for sum of children */
3361  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 2, monomials, 0.0, FALSE) );
3362 
3363  *op = SCIP_EXPR_POLYNOMIAL;
3364  data->data = (void*)polynomialdata;
3365 
3366  break;
3367  }
3368 
3369  case SCIP_EXPR_MINUS:
3370  {
3371  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3372  SCIP_EXPRDATA_MONOMIAL* monomials[2];
3373  int childidx;
3374  SCIP_Real exponent;
3375 
3376  assert(nchildren == 2);
3377 
3378  /* create monomial for first child */
3379  childidx = 0;
3380  exponent = 1.0;
3381  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomials[0], 1.0, 1, &childidx, &exponent) );
3382 
3383  /* create monomial for second child */
3384  childidx = 1;
3385  exponent = 1.0;
3386  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomials[1], -1.0, 1, &childidx, &exponent) );
3387 
3388  /* create polynomial for difference of children */
3389  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 2, monomials, 0.0, FALSE) );
3390 
3391  *op = SCIP_EXPR_POLYNOMIAL;
3392  data->data = (void*)polynomialdata;
3393 
3394  break;
3395  }
3396 
3397  case SCIP_EXPR_MUL:
3398  {
3399  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3400  SCIP_EXPRDATA_MONOMIAL* monomial;
3401  int childidx[2];
3402  SCIP_Real exponent[2];
3403 
3404  assert(nchildren == 2);
3405 
3406  /* create monomial for product of children */
3407  childidx[0] = 0;
3408  childidx[1] = 1;
3409  exponent[0] = 1.0;
3410  exponent[1] = 1.0;
3411  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 2, childidx, exponent) );
3412 
3413  /* create polynomial */
3414  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );
3415 
3416  *op = SCIP_EXPR_POLYNOMIAL;
3417  data->data = (void*)polynomialdata;
3418 
3419  break;
3420  }
3421 
3422  case SCIP_EXPR_DIV:
3423  {
3424  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3425  SCIP_EXPRDATA_MONOMIAL* monomial;
3426  int childidx[2];
3427  SCIP_Real exponent[2];
3428 
3429  assert(nchildren == 2);
3430 
3431  /* create monomial for division of children */
3432  childidx[0] = 0;
3433  childidx[1] = 1;
3434  exponent[0] = 1.0;
3435  exponent[1] = -1.0;
3436  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 2, childidx, exponent) );
3437 
3438  /* create polynomial */
3439  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );
3440 
3441  *op = SCIP_EXPR_POLYNOMIAL;
3442  data->data = (void*)polynomialdata;
3443 
3444  break;
3445  }
3446 
3447  case SCIP_EXPR_SQUARE:
3448  {
3449  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3450  SCIP_EXPRDATA_MONOMIAL* monomial;
3451  int childidx;
3452  SCIP_Real exponent;
3453 
3454  assert(nchildren == 1);
3455 
3456  /* create monomial for square of child */
3457  childidx = 0;
3458  exponent = 2.0;
3459  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &exponent) );
3460 
3461  /* create polynomial */
3462  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );
3463 
3464  *op = SCIP_EXPR_POLYNOMIAL;
3465  data->data = (void*)polynomialdata;
3466 
3467  break;
3468  }
3469 
3470  case SCIP_EXPR_SQRT:
3471  {
3472  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3473  SCIP_EXPRDATA_MONOMIAL* monomial;
3474  int childidx;
3475  SCIP_Real exponent;
3476 
3477  assert(nchildren == 1);
3478 
3479  /* create monomial for square root of child */
3480  childidx = 0;
3481  exponent = 0.5;
3482  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &exponent) );
3483 
3484  /* create polynomial */
3485  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );
3486 
3487  *op = SCIP_EXPR_POLYNOMIAL;
3488  data->data = (void*)polynomialdata;
3489 
3490  break;
3491  }
3492 
3493  case SCIP_EXPR_REALPOWER:
3494  {
3495  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3496  SCIP_EXPRDATA_MONOMIAL* monomial;
3497  int childidx;
3498 
3499  assert(nchildren == 1);
3500 
3501  /* convert to child0 to the power of exponent */
3502 
3503  /* create monomial for power of first child */
3504  childidx = 0;
3505  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &data->dbl) );
3506 
3507  /* create polynomial */
3508  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );
3509 
3510  *op = SCIP_EXPR_POLYNOMIAL;
3511  data->data = (void*)polynomialdata;
3512 
3513  break;
3514  }
3515 
3516  case SCIP_EXPR_SIGNPOWER:
3517  {
3518  SCIP_Real exponent;
3519 
3520  assert(nchildren == 1);
3521 
3522  /* check if exponent is an odd integer */
3523  exponent = data->dbl;
3524  if( EPSISINT(exponent, 0.0) && (int)exponent % 2 != 0 ) /*lint !e835*/
3525  {
3526  /* convert to child0 to the power of exponent, since sign is kept by taking power */
3527  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3528  SCIP_EXPRDATA_MONOMIAL* monomial;
3529  int childidx;
3530 
3531  /* create monomial for power of first child */
3532  childidx = 0;
3533  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &exponent) );
3534 
3535  /* create polynomial */
3536  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );
3537 
3538  *op = SCIP_EXPR_POLYNOMIAL;
3539  data->data = (void*)polynomialdata;
3540  }
3541  /* if exponent is not an odd integer constant, then keep it as signpower expression */
3542  break;
3543  }
3544 
3545  case SCIP_EXPR_INTPOWER:
3546  {
3547  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3548  SCIP_EXPRDATA_MONOMIAL* monomial;
3549  int childidx;
3550  SCIP_Real exponent;
3551 
3552  assert(nchildren == 1);
3553 
3554  /* create monomial for power of child */
3555  childidx = 0;
3556  exponent = data->intval;
3557  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &exponent) );
3558 
3559  /* create polynomial */
3560  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );
3561 
3562  *op = SCIP_EXPR_POLYNOMIAL;
3563  data->data = (void*)polynomialdata;
3564 
3565  break;
3566  }
3567 
3568  case SCIP_EXPR_EXP:
3569  case SCIP_EXPR_LOG:
3570  case SCIP_EXPR_SIN:
3571  case SCIP_EXPR_COS:
3572  case SCIP_EXPR_TAN:
3573  /* case SCIP_EXPR_ERF: */
3574  /* case SCIP_EXPR_ERFI: */
3575  case SCIP_EXPR_MIN:
3576  case SCIP_EXPR_MAX:
3577  case SCIP_EXPR_ABS:
3578  case SCIP_EXPR_SIGN:
3579  case SCIP_EXPR_USER:
3580  break;
3581 
3582  case SCIP_EXPR_SUM:
3583  {
3584  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3585  SCIP_EXPRDATA_MONOMIAL* monomial;
3586  int childidx;
3587  int i;
3588  SCIP_Real exponent;
3589 
3590  /* create empty polynomial */
3591  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 0, NULL, 0.0, FALSE) );
3592  SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, nchildren) );
3593  assert(polynomialdata->monomialssize >= nchildren);
3594 
3595  /* add summands as monomials */
3596  childidx = 0;
3597  exponent = 1.0;
3598  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &exponent) );
3599  for( i = 0; i < nchildren; ++i )
3600  {
3601  monomial->childidxs[0] = i;
3602  SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, 1, &monomial, TRUE) );
3603  }
3604  SCIPexprFreeMonomial(blkmem, &monomial);
3605 
3606  *op = SCIP_EXPR_POLYNOMIAL;
3607  data->data = (void*)polynomialdata;
3608 
3609  break;
3610  }
3611 
3612  case SCIP_EXPR_PRODUCT:
3613  {
3614  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3615  SCIP_EXPRDATA_MONOMIAL* monomial;
3616  int childidx;
3617  int i;
3618  SCIP_Real exponent;
3619 
3620  /* create monomial */
3621  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 0, NULL, NULL) );
3622  SCIP_CALL( monomialdataEnsureFactorsSize(blkmem, monomial, nchildren) );
3623  exponent = 1.0;
3624  for( i = 0; i < nchildren; ++i )
3625  {
3626  childidx = i;
3627  SCIP_CALL( SCIPexprAddMonomialFactors(blkmem, monomial, 1, &childidx, &exponent) );
3628  }
3629 
3630  /* create polynomial */
3631  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );
3632 
3633  *op = SCIP_EXPR_POLYNOMIAL;
3634  data->data = (void*)polynomialdata;
3635 
3636  break;
3637  }
3638 
3639  case SCIP_EXPR_LINEAR:
3640  {
3641  SCIP_Real* lineardata;
3642  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3643  SCIP_EXPRDATA_MONOMIAL* monomial;
3644  int childidx;
3645  int i;
3646  SCIP_Real exponent;
3647 
3648  /* get coefficients of linear term */
3649  lineardata = (SCIP_Real*)data->data;
3650  assert(lineardata != NULL);
3651 
3652  /* create polynomial consisting of constant from linear term */
3653  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 0, NULL, lineardata[nchildren], FALSE) );
3654  /* ensure space for linear coefficients */
3655  SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, nchildren) );
3656  assert(polynomialdata->monomialssize >= nchildren);
3657 
3658  /* add summands as monomials */
3659  childidx = 0;
3660  exponent = 1.0;
3661  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &exponent) );
3662  for( i = 0; i < nchildren; ++i )
3663  {
3664  monomial->coef = lineardata[i];
3665  monomial->childidxs[0] = i;
3666  SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, 1, &monomial, TRUE) );
3667  }
3668  SCIPexprFreeMonomial(blkmem, &monomial);
3669 
3670  /* free linear expression data */
3671  exprFreeDataLinear(blkmem, nchildren, *data);
3672 
3673  *op = SCIP_EXPR_POLYNOMIAL;
3674  data->data = (void*)polynomialdata;
3675 
3676  break;
3677  }
3678 
3679  case SCIP_EXPR_QUADRATIC:
3680  {
3681  SCIP_EXPRDATA_QUADRATIC* quaddata;
3682  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3683  SCIP_EXPRDATA_MONOMIAL* squaremonomial;
3684  SCIP_EXPRDATA_MONOMIAL* bilinmonomial;
3685  SCIP_EXPRDATA_MONOMIAL* linmonomial;
3686  int childidx[2];
3687  SCIP_Real exponent[2];
3688  int i;
3689 
3690  /* get data of quadratic expression */
3691  quaddata = (SCIP_EXPRDATA_QUADRATIC*)data->data;
3692  assert(quaddata != NULL);
3693 
3694  /* create empty polynomial */
3695  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 0, NULL, quaddata->constant, FALSE) );
3696  /* ensure space for linear and quadratic terms */
3697  SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, (quaddata->lincoefs != NULL ? nchildren : 0) + quaddata->nquadelems) );
3698  assert(polynomialdata->monomialssize >= quaddata->nquadelems);
3699 
3700  childidx[0] = 0;
3701  childidx[1] = 0;
3702 
3703  /* create monomial templates */
3704  exponent[0] = 2.0;
3705  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &squaremonomial, 1.0, 1, childidx, exponent) );
3706  exponent[0] = 1.0;
3707  exponent[1] = 1.0;
3708  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &bilinmonomial, 1.0, 2, childidx, exponent) );
3709  SCIP_CALL( SCIPexprCreateMonomial(blkmem, &linmonomial, 1.0, 1, childidx, exponent) );
3710 
3711  /* add linear terms as monomials */
3712  if( quaddata->lincoefs != NULL )
3713  for( i = 0; i < nchildren; ++i )
3714  if( quaddata->lincoefs[i] != 0.0 )
3715  {
3716  linmonomial->childidxs[0] = i;
3717  linmonomial->coef = quaddata->lincoefs[i];
3718  SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, 1, &linmonomial, TRUE) );
3719  }
3720 
3721  /* add quadratic terms as monomials */
3722  for( i = 0; i < quaddata->nquadelems; ++i )
3723  {
3724  if( quaddata->quadelems[i].idx1 == quaddata->quadelems[i].idx2 )
3725  {
3726  squaremonomial->childidxs[0] = quaddata->quadelems[i].idx1;
3727  squaremonomial->coef = quaddata->quadelems[i].coef;
3728  SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, 1, &squaremonomial, TRUE) );
3729  }
3730  else
3731  {
3732  bilinmonomial->childidxs[0] = quaddata->quadelems[i].idx1;
3733  bilinmonomial->childidxs[1] = quaddata->quadelems[i].idx2;
3734  bilinmonomial->coef = quaddata->quadelems[i].coef;
3735  SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, 1, &bilinmonomial, TRUE) );
3736  }
3737  }
3738  SCIPexprFreeMonomial(blkmem, &squaremonomial);
3739  SCIPexprFreeMonomial(blkmem, &bilinmonomial);
3740  SCIPexprFreeMonomial(blkmem, &linmonomial);
3741 
3742  /* free quadratic expression data */
3743  exprFreeDataQuadratic(blkmem, nchildren, *data);
3744 
3745  *op = SCIP_EXPR_POLYNOMIAL;
3746  data->data = (void*)polynomialdata;
3747 
3748  break;
3749  }
3750 
3751  case SCIP_EXPR_POLYNOMIAL:
3752  case SCIP_EXPR_LAST:
3753  break;
3754  } /*lint !e788*/
3755 
3756  return SCIP_OKAY;
3757 }
3758 
3759 /** converts polynomial expression back into simpler expression, if possible */
3760 static
3762  BMS_BLKMEM* blkmem, /**< block memory data structure */
3763  SCIP_EXPROP* op, /**< pointer to expression operator */
3764  SCIP_EXPROPDATA* data, /**< pointer to expression data holding polynomial data */
3765  int nchildren, /**< number of children of operator */
3766  void** children /**< children array */
3767  )
3768 {
3769  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
3770  SCIP_EXPRDATA_MONOMIAL* monomial;
3771  int maxdegree;
3772  int nlinmonomials;
3773  int i;
3774  int j;
3775 
3776  assert(blkmem != NULL);
3777  assert(op != NULL);
3778  assert(*op == SCIP_EXPR_POLYNOMIAL);
3779  assert(data != NULL);
3780  assert(children != NULL || nchildren == 0);
3781 
3782  polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)data->data;
3783  assert(polynomialdata != NULL);
3784 
3785  /* make sure monomials are sorted and merged */
3786  polynomialdataMergeMonomials(blkmem, polynomialdata, 0.0, TRUE);
3787 
3788  /* if no monomials, then leave as it is */
3789  if( polynomialdata->nmonomials == 0 )
3790  return SCIP_OKAY;
3791 
3792  /* check maximal degree of polynomial only - not considering children expressions
3793  * check number of linear monomials */
3794  maxdegree = 0;
3795  nlinmonomials = 0;
3796  for( i = 0; i < polynomialdata->nmonomials; ++i )
3797  {
3798  int monomialdegree;
3799 
3800  monomial = polynomialdata->monomials[i];
3801  assert(monomial != NULL);
3802 
3803  monomialdegree = 0;
3804  for(j = 0; j < monomial->nfactors; ++j )
3805  {
3806  if( !EPSISINT(monomial->exponents[j], 0.0) || monomial->exponents[j] < 0.0 ) /*lint !e835*/
3807  {
3808  monomialdegree = SCIP_EXPR_DEGREEINFINITY;
3809  break;
3810  }
3811 
3812  monomialdegree += (int)EPSROUND(monomial->exponents[j], 0.0); /*lint !e835*/
3813  }
3814 
3815  if( monomialdegree == SCIP_EXPR_DEGREEINFINITY )
3816  {
3817  maxdegree = SCIP_EXPR_DEGREEINFINITY;
3818  break;
3819  }
3820 
3821  if( monomialdegree == 1 )
3822  ++nlinmonomials;
3823 
3824  if( monomialdegree > maxdegree )
3825  maxdegree = monomialdegree;
3826  }
3827  assert(maxdegree > 0 );
3828 
3829  if( maxdegree == 1 )
3830  {
3831  /* polynomial is a linear expression in children */
3832 
3833  /* polynomial simplification and monomial merging should ensure that monomial i corresponds to child i and that there are not unused children */
3834  assert(polynomialdata->nmonomials == nchildren);
3835  assert(polynomialdata->nmonomials == nlinmonomials);
3836 
3837  if( polynomialdata->constant == 0.0 && polynomialdata->nmonomials == 2 && polynomialdata->monomials[0]->coef == 1.0 && polynomialdata->monomials[1]->coef == 1.0 )
3838  {
3839  /* polynomial is addition of two expressions, so turn into SCIP_EXPR_PLUS */
3840  assert(polynomialdata->monomials[0]->nfactors == 1);
3841  assert(polynomialdata->monomials[0]->exponents[0] == 1.0);
3842  assert(polynomialdata->monomials[1]->nfactors == 1);
3843  assert(polynomialdata->monomials[1]->exponents[0] == 1.0);
3844 
3845  polynomialdataFree(blkmem, &polynomialdata);
3846  data->data = NULL;
3847 
3848  /* change operator type to PLUS */
3849  *op = SCIP_EXPR_PLUS;
3850 
3851  return SCIP_OKAY;
3852  }
3853 
3854  if( polynomialdata->constant == 0.0 && polynomialdata->nmonomials == 2 && polynomialdata->monomials[0]->coef == 1.0 && polynomialdata->monomials[1]->coef == -1.0 )
3855  {
3856  /* polynomial is substraction of two expressions, so turn into SCIP_EXPR_MINUS */
3857  assert(polynomialdata->monomials[0]->nfactors == 1);
3858  assert(polynomialdata->monomials[0]->exponents[0] == 1.0);
3859  assert(polynomialdata->monomials[1]->nfactors == 1);
3860  assert(polynomialdata->monomials[1]->exponents[0] == 1.0);
3861 
3862  polynomialdataFree(blkmem, &polynomialdata);
3863  data->data = NULL;
3864 
3865  /* change operator type to MINUS */
3866  *op = SCIP_EXPR_MINUS;
3867 
3868  return SCIP_OKAY;
3869  }
3870 
3871  if( polynomialdata->constant == 0.0 && polynomialdata->nmonomials == 2 && polynomialdata->monomials[0]->coef == -1.0 && polynomialdata->monomials[1]->coef == 1.0 )
3872  {
3873  /* polynomial is substraction of two expressions, so turn into SCIP_EXPR_MINUS */
3874  void* tmp;
3875 
3876  assert(polynomialdata->monomials[0]->nfactors == 1);
3877  assert(polynomialdata->monomials[0]->exponents[0] == 1.0);
3878  assert(polynomialdata->monomials[1]->nfactors == 1);
3879  assert(polynomialdata->monomials[1]->exponents[0] == 1.0);
3880 
3881  polynomialdataFree(blkmem, &polynomialdata);
3882  data->data = NULL;
3883 
3884  /* swap children */
3885  tmp = children[1]; /*lint !e613*/
3886  children[1] = children[0]; /*lint !e613*/
3887  children[0] = tmp; /*lint !e613*/
3888 
3889  /* change operator type to MINUS */
3890  *op = SCIP_EXPR_MINUS;
3891 
3892  return SCIP_OKAY;
3893  }
3894 
3895  if( polynomialdata->constant == 0.0 )
3896  {
3897  /* check if all monomials have coefficient 1.0 */
3898  for( i = 0; i < polynomialdata->nmonomials; ++i )
3899  if( polynomialdata->monomials[i]->coef != 1.0 )
3900  break;
3901 
3902  if( i == polynomialdata->nmonomials )
3903  {
3904  /* polynomial is sum of children, so turn into SCIP_EXPR_SUM */
3905 
3906  polynomialdataFree(blkmem, &polynomialdata);
3907  data->data = NULL;
3908 
3909  /* change operator type to MINUS */
3910  *op = SCIP_EXPR_SUM;
3911 
3912  return SCIP_OKAY;
3913  }
3914  }
3915 
3916  /* turn polynomial into linear expression */
3917  {
3918  SCIP_Real* lindata;
3919 
3920  /* monomial merging should ensure that each child appears in at most one monomial,
3921  * that monomials are ordered according to the child index, and that constant monomials have been removed
3922  */
3923 
3924  /* setup data of linear expression */
3925  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &lindata, polynomialdata->nmonomials + 1) );
3926 
3927  for( i = 0; i < polynomialdata->nmonomials; ++i )
3928  {
3929  assert(polynomialdata->monomials[i]->childidxs[0] == i);
3930  assert(polynomialdata->monomials[i]->exponents[0] == 1.0);
3931  lindata[i] = polynomialdata->monomials[i]->coef; /*lint !e644*/
3932  }
3933  lindata[i] = polynomialdata->constant;
3934 
3935  polynomialdataFree(blkmem, &polynomialdata);
3936  *op = SCIP_EXPR_LINEAR;
3937  data->data = (void*)lindata;
3938 
3939  return SCIP_OKAY;
3940  }
3941  }
3942 
3943  if( maxdegree == 2 && (polynomialdata->nmonomials > 1 || polynomialdata->constant != 0.0 || polynomialdata->monomials[0]->coef != 1.0) )
3944  {
3945  /* polynomial is quadratic expression with more than one summand or with a constant or a square or bilinear term with coefficient != 1.0, so turn into SCIP_EXPR_QUADRATIC */
3946  SCIP_EXPRDATA_QUADRATIC* quaddata;
3947  int quadelemidx;
3948 
3949  SCIP_ALLOC( BMSallocBlockMemory(blkmem, &quaddata) );
3950  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &quaddata->quadelems, polynomialdata->nmonomials - nlinmonomials) );
3951  quaddata->nquadelems = polynomialdata->nmonomials - nlinmonomials;
3952  quaddata->constant = polynomialdata->constant;
3953  quaddata->sorted = FALSE; /* quadratic data is sorted different than polynomials */
3954 
3955  if( nlinmonomials > 0 )
3956  {
3957  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &quaddata->lincoefs, nchildren) );
3958  BMSclearMemoryArray(quaddata->lincoefs, nchildren);
3959  }
3960  else
3961  quaddata->lincoefs = NULL;
3962 
3963  quadelemidx = 0;
3964  for( i = 0; i < polynomialdata->nmonomials; ++i )
3965  {
3966  assert(polynomialdata->monomials[i]->nfactors == 1 || polynomialdata->monomials[i]->nfactors == 2);
3967  if( polynomialdata->monomials[i]->nfactors == 1 )
3968  {
3969  if( polynomialdata->monomials[i]->exponents[0] == 1.0 )
3970  {
3971  /* monomial is a linear term */
3972  assert(quaddata->lincoefs != NULL);
3973  quaddata->lincoefs[polynomialdata->monomials[i]->childidxs[0]] += polynomialdata->monomials[i]->coef;
3974  }
3975  else
3976  {
3977  /* monomial should be a square term */
3978  assert(polynomialdata->monomials[i]->exponents[0] == 2.0);
3979  assert(quadelemidx < quaddata->nquadelems);
3980  quaddata->quadelems[quadelemidx].idx1 = polynomialdata->monomials[i]->childidxs[0];
3981  quaddata->quadelems[quadelemidx].idx2 = polynomialdata->monomials[i]->childidxs[0];
3982  quaddata->quadelems[quadelemidx].coef = polynomialdata->monomials[i]->coef;
3983  ++quadelemidx;
3984  }
3985  }
3986  else
3987  {
3988  /* monomial should be a bilinear term */
3989  assert(polynomialdata->monomials[i]->exponents[0] == 1.0);
3990  assert(polynomialdata->monomials[i]->exponents[1] == 1.0);
3991  assert(quadelemidx < quaddata->nquadelems);
3992  quaddata->quadelems[quadelemidx].idx1 = MIN(polynomialdata->monomials[i]->childidxs[0], polynomialdata->monomials[i]->childidxs[1]);
3993  quaddata->quadelems[quadelemidx].idx2 = MAX(polynomialdata->monomials[i]->childidxs[0], polynomialdata->monomials[i]->childidxs[1]);
3994  quaddata->quadelems[quadelemidx].coef = polynomialdata->monomials[i]->coef;
3995  ++quadelemidx;
3996  }
3997  }
3998  assert(quadelemidx == quaddata->nquadelems);
3999 
4000  polynomialdataFree(blkmem, &polynomialdata);
4001 
4002  *op = SCIP_EXPR_QUADRATIC;
4003  data->data = (void*)quaddata;
4004 
4005  return SCIP_OKAY;
4006  }
4007 
4008  if( polynomialdata->constant == 0.0 && polynomialdata->nmonomials == 1 && polynomialdata->monomials[0]->coef == 1.0 )
4009  {
4010  /* polynomial is product of children */
4011  monomial = polynomialdata->monomials[0];
4012  assert(monomial->nfactors == nchildren);
4013 
4014  if( monomial->nfactors == 1 )
4015  {
4016  /* polynomial is x^k for some k */
4017  assert(monomial->exponents[0] != 1.0); /* should have been handled before */
4018  assert(monomial->childidxs[0] == 0);
4019 
4020  if( monomial->exponents[0] == 2.0 )
4021  {
4022  /* polynomial is x^2, so turn into SCIP_EXPR_SQUARE */
4023 
4024  polynomialdataFree(blkmem, &polynomialdata);
4025  data->data = NULL;
4026 
4027  *op = SCIP_EXPR_SQUARE;
4028 
4029  return SCIP_OKAY;
4030  }
4031 
4032  if( EPSISINT(monomial->exponents[0], 0.0) ) /*lint !e835*/
4033  {
4034  /* k is an integer, so turn into SCIP_EXPR_INTPOWER */
4035  int exponent;
4036 
4037  exponent = (int)EPSROUND(monomial->exponents[0], 0.0); /*lint !e835*/
4038 
4039  polynomialdataFree(blkmem, &polynomialdata);
4040 
4041  *op = SCIP_EXPR_INTPOWER;
4042  data->intval = exponent;
4043 
4044  return SCIP_OKAY;
4045  }
4046 
4047  if( monomial->exponents[0] == 0.5 )
4048  {
4049  /* polynomial is sqrt(x), so turn into SCIP_EXPR_SQRT */
4050 
4051  polynomialdataFree(blkmem, &polynomialdata);
4052  data->data = NULL;
4053 
4054  *op = SCIP_EXPR_SQRT;
4055 
4056  return SCIP_OKAY;
4057  }
4058 
4059  {
4060  /* polynomial is x^a with a some real number, so turn into SCIP_EXPR_REALPOWER */
4061  SCIP_Real exponent;
4062 
4063  exponent = monomial->exponents[0];
4064 
4065  polynomialdataFree(blkmem, &polynomialdata);
4066 
4067  *op = SCIP_EXPR_REALPOWER;
4068  data->dbl = exponent;
4069 
4070  return SCIP_OKAY;
4071  }
4072  }
4073 
4074  if( maxdegree == 2 && monomial->nfactors == 2 )
4075  {
4076  /* polynomial is product of two children, so turn into SCIP_EXPR_MUL */
4077  assert(monomial->exponents[0] == 1.0);
4078  assert(monomial->exponents[1] == 1.0);
4079 
4080  polynomialdataFree(blkmem, &polynomialdata);
4081  data->data = NULL;
4082 
4083  *op = SCIP_EXPR_MUL;
4084 
4085  return SCIP_OKAY;
4086  }
4087 
4088  if( maxdegree == monomial->nfactors )
4089  {
4090  /* polynomial is a product of n children, so turn into SCIP_EXPR_PRODUCT */
4091 
4092  polynomialdataFree(blkmem, &polynomialdata);
4093  data->data = NULL;
4094 
4095  *op = SCIP_EXPR_PRODUCT;
4096 
4097  return SCIP_OKAY;
4098  }
4099 
4100  if( monomial->nfactors == 2 && monomial->exponents[0] == 1.0 && monomial->exponents[1] == -1.0 )
4101  {
4102  /* polynomial is x/y, so turn into SCIP_EXPR_DIV */
4103  assert(monomial->childidxs[0] == 0);
4104  assert(monomial->childidxs[1] == 1);
4105 
4106  polynomialdataFree(blkmem, &polynomialdata);
4107  data->data = NULL;
4108 
4109  *op = SCIP_EXPR_DIV;
4110 
4111  return SCIP_OKAY;
4112  }
4113 
4114  if( monomial->nfactors == 2 && monomial->exponents[0] == -1.0 && monomial->exponents[1] == 1.0 )
4115  {
4116  /* polynomial is y/x, so turn into SCIP_EXPR_DIV */
4117  void* tmp;
4118 
4119  assert(monomial->childidxs[0] == 0);
4120  assert(monomial->childidxs[1] == 1);
4121 
4122  polynomialdataFree(blkmem, &polynomialdata);
4123  data->data = NULL;
4124 
4125  /* swap children */
4126  tmp = children[1]; /*lint !e613*/
4127  children[1] = children[0]; /*lint !e613*/
4128  children[0] = tmp; /*lint !e613*/
4129 
4130  *op = SCIP_EXPR_DIV;
4131 
4132  return SCIP_OKAY;
4133  }
4134  }
4135 
4136  return SCIP_OKAY;
4137 }
4138 
4139 /** adds copies of expressions to the array of children of a sum, product, linear, quadratic, or polynomial expression
4140  *
4141  * For a sum or product expression, this corresponds to add additional summands and factors, resp.
4142  * For a linear expression, this corresponds to add each expression with coefficient 1.0.
4143  * For a quadratic or polynomial expression, only the children array may be enlarged, the expression itself remains the same.
4144  */
4145 static
4147  BMS_BLKMEM* blkmem, /**< block memory */
4148  SCIP_EXPR* expr, /**< quadratic or polynomial expression */
4149  int nexprs, /**< number of expressions to add */
4150  SCIP_EXPR** exprs, /**< expressions to add */
4151  SCIP_Bool comparechildren, /**< whether to compare expressions with already existing children (no effect for sum and product) */
4152  SCIP_Real eps, /**< which epsilon to use when comparing expressions */
4153  int* childmap /**< array where to store mapping of indices from exprs to children array in expr, or NULL if not of interest */
4154  )
4155 {
4156  int i;
4157 
4158  assert(blkmem != NULL);
4159  assert(expr != NULL);
4160  assert(expr->op == SCIP_EXPR_SUM || expr->op == SCIP_EXPR_PRODUCT || expr->op == SCIP_EXPR_LINEAR || expr->op == SCIP_EXPR_QUADRATIC || expr->op == SCIP_EXPR_POLYNOMIAL);
4161  assert(exprs != NULL || nexprs == 0);
4162 
4163  if( nexprs == 0 )
4164  return SCIP_OKAY;
4165 
4166  switch( expr->op )
4167  {
4168  case SCIP_EXPR_SUM:
4169  case SCIP_EXPR_PRODUCT:
4170  {
4171  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, expr->nchildren, expr->nchildren + nexprs) );
4172  for( i = 0; i < nexprs; ++i )
4173  {
4174  SCIP_CALL( SCIPexprCopyDeep(blkmem, &expr->children[expr->nchildren + i], exprs[i]) ); /*lint !e613*/
4175  if( childmap != NULL )
4176  childmap[i] = expr->nchildren + i;
4177  }
4178  expr->nchildren += nexprs;
4179 
4180  break;
4181  }
4182 
4183  case SCIP_EXPR_LINEAR:
4184  case SCIP_EXPR_QUADRATIC:
4185  case SCIP_EXPR_POLYNOMIAL:
4186  {
4187  int j;
4188  int orignchildren;
4189  SCIP_Bool existsalready;
4190 
4191  orignchildren = expr->nchildren;
4192  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, expr->nchildren, expr->nchildren + nexprs) );
4193 
4194  for( i = 0; i < nexprs; ++i )
4195  {
4196  existsalready = FALSE;
4197  if( comparechildren )
4198  for( j = 0; j < orignchildren; ++j )
4199  /* during simplification of polynomials, their may be NULL's in children array */
4200  if( expr->children[j] != NULL && SCIPexprAreEqual(expr->children[j], exprs[i], eps) ) /*lint !e613*/
4201  {
4202  existsalready = TRUE;
4203  break;
4204  }
4205 
4206  if( !existsalready )
4207  {
4208  /* add copy of exprs[j] to children array */
4209  SCIP_CALL( SCIPexprCopyDeep(blkmem, &expr->children[expr->nchildren], exprs[i]) ); /*lint !e613*/
4210  if( childmap != NULL )
4211  childmap[i] = expr->nchildren;
4212  ++expr->nchildren;
4213  }
4214  else
4215  {
4216  if( childmap != NULL )
4217  childmap[i] = j; /*lint !e644*/
4218  if( expr->op == SCIP_EXPR_LINEAR )
4219  {
4220  /* if linear expression, increase coefficient by 1.0 */
4221  ((SCIP_Real*)expr->data.data)[j] += 1.0;
4222  }
4223  }
4224  }
4225 
4226  /* shrink children array to actually used size */
4227  assert(comparechildren || expr->nchildren == orignchildren + nexprs);
4228  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, orignchildren + nexprs, expr->nchildren) );
4229 
4230  if( expr->op == SCIP_EXPR_LINEAR && expr->nchildren > orignchildren )
4231  {
4232  /* if linear expression, then add 1.0 coefficients for new expressions */
4233  SCIP_Real* data;
4234 
4235  data = (SCIP_Real*)expr->data.data;
4236  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &data, orignchildren + 1, expr->nchildren + 1) );
4237  data[expr->nchildren] = data[orignchildren]; /* move constant from old end to new end */
4238  for( i = orignchildren; i < expr->nchildren; ++i )
4239  data[i] = 1.0;
4240  expr->data.data = (void*)data;
4241  }
4242  else if( expr->op == SCIP_EXPR_QUADRATIC && expr->nchildren > orignchildren )
4243  {
4244  /* if quadratic expression, then add 0.0 linear coefficients for new expressions */
4246 
4247  data = (SCIP_EXPRDATA_QUADRATIC*)expr->data.data;
4248  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &data->lincoefs, orignchildren, expr->nchildren) );
4249  BMSclearMemoryArray(&data->lincoefs[orignchildren], expr->nchildren - orignchildren); /*lint !e866*/
4250  }
4251 
4252  break;
4253  }
4254 
4255  default:
4256  SCIPerrorMessage("exprsimplifyAddChildren cannot be called for operand %d\n", expr->op);
4257  return SCIP_INVALIDDATA;
4258  } /*lint !e788*/
4259 
4260  return SCIP_OKAY;
4261 }
4262 
4263 /** converts expressions into polynomials, where possible and obvious */
4264 static
4266  BMS_BLKMEM* blkmem, /**< block memory data structure */
4267  SCIP_EXPR* expr /**< expression to convert */
4268  )
4269 {
4270  int i;
4271 
4272  assert(expr != NULL);
4273 
4274  for( i = 0; i < expr->nchildren; ++i )
4275  {
4277  }
4278 
4279  SCIP_CALL( exprConvertToPolynomial(blkmem, &expr->op, &expr->data, expr->nchildren) );
4280 
4281  return SCIP_OKAY;
4282 }
4283 
4284 /** removes duplicate children in a polynomial expression
4285  *
4286  * Leaves NULL's in children array.
4287  */
4288 static
4290  BMS_BLKMEM* blkmem, /**< block memory data structure */
4291  SCIP_EXPR* expr, /**< expression */
4292  SCIP_Real eps /**< threshold for zero */
4293  )
4294 {
4295  SCIP_Bool foundduplicates;
4296  int* childmap;
4297  int i;
4298  int j;
4299 
4300  assert(blkmem != NULL);
4301  assert(expr != NULL);
4302  assert(SCIPexprGetOperator(expr) == SCIP_EXPR_POLYNOMIAL);
4303 
4304  if( expr->nchildren == 0 )
4305  return SCIP_OKAY;
4306 
4307  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &childmap, expr->nchildren) );
4308 
4309  foundduplicates = FALSE;
4310  for( i = 0; i < expr->nchildren; ++i )
4311  {
4312  if( expr->children[i] == NULL )
4313  continue;
4314  childmap[i] = i; /*lint !e644*/
4315 
4316  for( j = i+1; j < expr->nchildren; ++j )
4317  {
4318  if( expr->children[j] == NULL )
4319  continue;
4320 
4321  if( SCIPexprAreEqual(expr->children[i], expr->children[j], eps) )
4322  {
4323  /* forget about expr j and remember that is to be replaced by i */
4324  SCIPexprFreeDeep(blkmem, &expr->children[j]);
4325  childmap[j] = i;
4326  foundduplicates = TRUE;
4327  }
4328  }
4329  }
4330 
4331  /* apply childmap to monomials */
4332  if( foundduplicates )
4334 
4335  /* free childmap */
4336  BMSfreeBlockMemoryArray(blkmem, &childmap, expr->nchildren);
4337 
4338  return SCIP_OKAY;
4339 }
4340 
4341 /** eliminates NULL's in children array and shrinks it to actual size */
4342 static
4344  BMS_BLKMEM* blkmem, /**< block memory data structure */
4345  SCIP_EXPR* expr /**< expression */
4346  )
4347 {
4348  int* childmap;
4349  int lastnonnull;
4350  int i;
4351 
4352  assert(blkmem != NULL);
4353  assert(expr != NULL);
4354  assert(SCIPexprGetOperator(expr) == SCIP_EXPR_POLYNOMIAL);
4355 
4356  if( expr->nchildren == 0 )
4357  return SCIP_OKAY;
4358 
4359  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &childmap, expr->nchildren) );
4360 
4361  /* close gaps in children array */
4362  lastnonnull = expr->nchildren-1;
4363  while( lastnonnull >= 0 && expr->children[lastnonnull] == NULL )
4364  --lastnonnull;
4365  for( i = 0; i <= lastnonnull; ++i )
4366  {
4367  if( expr->children[i] != NULL )
4368  {
4369  childmap[i] = i; /* child at index i is not moved */ /*lint !e644*/
4370  continue;
4371  }
4372  assert(expr->children[lastnonnull] != NULL);
4373 
4374  /* move child at lastnonnull to position i */
4375  expr->children[i] = expr->children[lastnonnull];
4376  expr->children[lastnonnull] = NULL;
4377  childmap[lastnonnull] = i;
4378 
4379  /* update lastnonnull */
4380  --lastnonnull;
4381  while( lastnonnull >= 0 && expr->children[lastnonnull] == NULL )
4382  --lastnonnull;
4383  }
4384  assert(i > lastnonnull);
4385 
4386  /* apply childmap to monomials */
4387  if( lastnonnull < expr->nchildren-1 )
4389 
4390  BMSfreeBlockMemoryArray(blkmem, &childmap, expr->nchildren);
4391 
4392  /* shrink children array */
4393  if( lastnonnull >= 0 )
4394  {
4395  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, expr->nchildren, lastnonnull+1) );
4396  expr->nchildren = lastnonnull+1;
4397  }
4398  else
4399  {
4400  BMSfreeBlockMemoryArray(blkmem, &expr->children, expr->nchildren);
4401  expr->nchildren = 0;
4402  }
4403 
4404  return SCIP_OKAY;
4405 }
4406 
4407 /** checks which children are still in use and frees those which are not */
4408 static
4410  BMS_BLKMEM* blkmem, /**< block memory data structure */
4411  SCIP_EXPR* expr /**< polynomial expression */
4412  )
4413 {
4414  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
4415  SCIP_EXPRDATA_MONOMIAL* monomial;
4416  SCIP_Bool* childinuse;
4417  int i;
4418  int j;
4419 
4420  assert(blkmem != NULL);
4421  assert(expr != NULL);
4422 
4423  if( expr->nchildren == 0 )
4424  return SCIP_OKAY;
4425 
4426  polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)expr->data.data;
4427  assert(polynomialdata != NULL);
4428 
4429  /* check which children are still in use */
4430  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &childinuse, expr->nchildren) );
4431  BMSclearMemoryArray(childinuse, expr->nchildren); /*lint !e644*/
4432  for( i = 0; i < polynomialdata->nmonomials; ++i )
4433  {
4434  monomial = polynomialdata->monomials[i];
4435  assert(monomial != NULL);
4436 
4437  for( j = 0; j < monomial->nfactors; ++j )
4438  {
4439  assert(monomial->childidxs[j] >= 0);
4440  assert(monomial->childidxs[j] < expr->nchildren);
4441  childinuse[monomial->childidxs[j]] = TRUE;
4442  }
4443  }
4444 
4445  /* free children that are not used in any monomial */
4446  for( i = 0; i < expr->nchildren; ++i )
4447  if( expr->children[i] != NULL && !childinuse[i] )
4448  SCIPexprFreeDeep(blkmem, &expr->children[i]);
4449 
4450  BMSfreeBlockMemoryArray(blkmem, &childinuse, expr->nchildren);
4451 
4452  return SCIP_OKAY;
4453 }
4454 
4455 /** flattens polynomials in polynomials, check for constants in non-polynomials expressions
4456  *
4457  * exprsimplifyConvertToPolynomials should have been called before to eliminate simple polynomial operands.
4458  */
4459 static
4461  BMS_BLKMEM* blkmem, /**< block memory data structure */
4462  SCIP_MESSAGEHDLR* messagehdlr, /**< message handler */
4463  SCIP_EXPR* expr, /**< expression */
4464  SCIP_Real eps, /**< threshold, under which values are treat as 0 */
4465  int maxexpansionexponent/**< maximal exponent for which we still expand non-monomial polynomials */
4466  )
4467 {
4468  int i;
4469 
4470  assert(expr != NULL);
4471 
4472  for( i = 0; i < expr->nchildren; ++i )
4473  {
4474  SCIP_CALL( exprsimplifyFlattenPolynomials(blkmem, messagehdlr, expr->children[i], eps, maxexpansionexponent) );
4475  }
4476 
4477  switch( SCIPexprGetOperator(expr) )
4478  {
4479  case SCIP_EXPR_VARIDX:
4480  case SCIP_EXPR_CONST:
4481  case SCIP_EXPR_PARAM:
4482  case SCIP_EXPR_PLUS:
4483  case SCIP_EXPR_MINUS:
4484  case SCIP_EXPR_MUL:
4485  case SCIP_EXPR_DIV:
4486  case SCIP_EXPR_SQUARE:
4487  case SCIP_EXPR_SQRT:
4488  case SCIP_EXPR_INTPOWER:
4489  case SCIP_EXPR_REALPOWER:
4490  case SCIP_EXPR_SIGNPOWER:
4491  break;
4492 
4493  case SCIP_EXPR_EXP:
4494  case SCIP_EXPR_LOG:
4495  case SCIP_EXPR_SIN:
4496  case SCIP_EXPR_COS:
4497  case SCIP_EXPR_TAN:
4498  /* case SCIP_EXPR_ERF: */
4499  /* case SCIP_EXPR_ERFI: */
4500  case SCIP_EXPR_ABS:
4501  case SCIP_EXPR_SIGN:
4502  {
4503  /* check if argument is a constant */
4504  if( (expr->children[0]->op == SCIP_EXPR_POLYNOMIAL && SCIPexprGetNChildren(expr->children[0]) == 0) ||
4505  expr->children[0]->op == SCIP_EXPR_CONST )
4506  {
4507  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
4508  SCIP_Real exprval;
4509 
4510  /* since child0 has no children and it's polynomial was flattened, it should have no monomials */
4511  assert(expr->children[0]->op != SCIP_EXPR_POLYNOMIAL || SCIPexprGetNMonomials(expr->children[0]) == 0);
4512 
4513  /* evaluate expression in constant polynomial */
4514  SCIP_CALL( SCIPexprEval(expr, NULL, NULL, &exprval) );
4515 
4516  /* create polynomial */
4517  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 0, NULL, exprval, FALSE) );
4518 
4519  expr->op = SCIP_EXPR_POLYNOMIAL;
4520  expr->data.data = (void*)polynomialdata;
4521 
4522  /* forget child */
4523  SCIPexprFreeDeep(blkmem, &expr->children[0]);
4524  BMSfreeBlockMemoryArray(blkmem, &expr->children, 1);
4525  expr->nchildren = 0;
4526  }
4527 
4528  break;
4529  }
4530 
4531  case SCIP_EXPR_MIN:
4532  case SCIP_EXPR_MAX:
4533  {
4534  /* check if both arguments are constants */
4535  if( ((expr->children[0]->op == SCIP_EXPR_POLYNOMIAL && SCIPexprGetNChildren(expr->children[0]) == 0) || expr->children[0]->op == SCIP_EXPR_CONST) &&
4536  ((expr->children[1]->op == SCIP_EXPR_POLYNOMIAL && SCIPexprGetNChildren(expr->children[1]) == 0) || expr->children[1]->op == SCIP_EXPR_CONST) )
4537  {
4538  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
4539  SCIP_Real exprval;
4540 
4541  /* since children have no children and it's polynomial was flattened, it should have no monomials */
4542  assert(expr->children[0]->op != SCIP_EXPR_POLYNOMIAL || SCIPexprGetNMonomials(expr->children[0]) == 0);
4543  assert(expr->children[1]->op != SCIP_EXPR_POLYNOMIAL || SCIPexprGetNMonomials(expr->children[1]) == 0);
4544 
4545  /* evaluate expression in constants */
4546  SCIP_CALL( SCIPexprEval(expr, NULL, NULL, &exprval) );
4547 
4548  /* create polynomial */
4549  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 0, NULL, exprval, FALSE) );
4550 
4551  expr->op = SCIP_EXPR_POLYNOMIAL;
4552  expr->data.data = (void*)polynomialdata;
4553 
4554  /* forget children */
4555  SCIPexprFreeDeep(blkmem, &expr->children[0]);
4556  SCIPexprFreeDeep(blkmem, &expr->children[1]);
4557  BMSfreeBlockMemoryArray(blkmem, &expr->children, 2);
4558  expr->nchildren = 0;
4559  }
4560 
4561  break;
4562  }
4563 
4564  case SCIP_EXPR_SUM:
4565  case SCIP_EXPR_PRODUCT:
4566  case SCIP_EXPR_LINEAR:
4567  case SCIP_EXPR_QUADRATIC:
4568  case SCIP_EXPR_USER:
4569  break;
4570 
4571  case SCIP_EXPR_POLYNOMIAL:
4572  {
4573  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
4574  SCIP_EXPRDATA_MONOMIAL* monomial;
4575  SCIP_Bool removechild;
4576  int* childmap;
4577  int childmapsize;
4578  int j;
4579 
4580  /* simplify current polynomial */
4582  SCIPexprMergeMonomials(blkmem, expr, eps, TRUE);
4583 
4584  polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)expr->data.data;
4585  assert(polynomialdata != NULL);
4586 
4587  SCIPdebugMessage("expand factors in expression ");
4588  SCIPdebug( SCIPexprPrint(expr, messagehdlr, NULL, NULL, NULL, NULL) );
4589  SCIPdebugPrintf("\n");
4590 
4591  childmap = NULL;
4592  childmapsize = 0;
4593 
4594  /* resolve children that are constants
4595  * we do this first, because it reduces the degree and number of factors in the monomials,
4596  * thereby allowing some expansions of polynomials that may not be possible otherwise, e.g., turning c0*c1 with c0=quadratic and c1=constant into a single monomial
4597  */
4598  for( i = 0; i < expr->nchildren; ++i )
4599  {
4600  if( expr->children[i] == NULL )
4601  continue;
4602 
4603  if( SCIPexprGetOperator(expr->children[i]) != SCIP_EXPR_CONST )
4604  continue;
4605 
4606  removechild = TRUE; /* we intend to delete children[i] */
4607 
4608  if( childmapsize < expr->children[i]->nchildren )
4609  {
4610  int newsize;
4611 
4612  newsize = calcGrowSize(expr->children[i]->nchildren);
4613  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &childmap, childmapsize, newsize) );
4614  childmapsize = newsize;
4615  }
4616 
4617  /* put constant of child i into every monomial where child i is used */
4618  for( j = 0; j < polynomialdata->nmonomials; ++j )
4619  {
4620  int factorpos;
4621  SCIP_Bool success;
4622 
4623  monomial = polynomialdata->monomials[j];
4624  /* if monomial is not sorted, then polynomial should not be sorted either, or have only one monomial */
4625  assert(monomial->sorted || !polynomialdata->sorted || polynomialdata->nmonomials <= 1);
4626 
4627  if( SCIPexprFindMonomialFactor(monomial, i, &factorpos) )
4628  {
4629  assert(factorpos >= 0);
4630  assert(factorpos < monomial->nfactors);
4631  /* assert that factors have been merged */
4632  assert(factorpos == 0 || monomial->childidxs[factorpos-1] != i);
4633  assert(factorpos == monomial->nfactors-1 || monomial->childidxs[factorpos+1] != i);
4634 
4635  /* SCIPdebugMessage("attempt expanding child %d at monomial %d factor %d\n", i, j, factorpos);
4636  SCIPdebug( SCIPexprPrint(expr, NULL, NULL, NULL) ); SCIPdebugPrintf("\n");
4637  SCIPdebug( SCIPexprPrint(expr->children[i], NULL, NULL, NULL) ); SCIPdebugPrintf("\n"); */
4638 
4639  if( !EPSISINT(monomial->exponents[factorpos], 0.0) && SCIPexprGetOpReal(expr->children[i]) < 0.0 ) /*lint !e835*/
4640  {
4641  /* if constant is negative and our exponent is not integer, then cannot do expansion */
4642  SCIPmessagePrintWarning(messagehdlr, "got negative constant %g to the power of a noninteger exponent %g\n",
4643  SCIPexprGetOpReal(expr->children[i]), monomial->exponents[factorpos]);
4644  success = FALSE;
4645  }
4646  else
4647  {
4648  monomial->coef *= pow(SCIPexprGetOpReal(expr->children[i]), monomial->exponents[factorpos]);
4649 
4650  /* move last factor to position factorpos */
4651  if( factorpos < monomial->nfactors-1 )
4652  {
4653  monomial->exponents[factorpos] = monomial->exponents[monomial->nfactors-1];
4654  monomial->childidxs[factorpos] = monomial->childidxs[monomial->nfactors-1];
4655  }
4656  --monomial->nfactors;
4657  monomial->sorted = FALSE;
4658  polynomialdata->sorted = FALSE;
4659 
4660  success = TRUE;
4661  }
4662 
4663  if( !success )
4664  removechild = FALSE;
4665  }
4666  }
4667 
4668  /* forget about child i, if it is not used anymore */
4669  if( removechild )
4670  SCIPexprFreeDeep(blkmem, &expr->children[i]);
4671 
4672  /* simplify current polynomial again */
4673  SCIPexprMergeMonomials(blkmem, expr, eps, TRUE);
4674  }
4675 
4676  /* try to resolve children that are polynomials itself */
4677  for( i = 0; i < expr->nchildren; ++i )
4678  {
4679  if( expr->children[i] == NULL )
4680  continue;
4681 
4683  continue;
4684 
4685  removechild = TRUE; /* we intend to delete children[i] */
4686 
4687  if( childmapsize < expr->children[i]->nchildren )
4688  {
4689  int newsize;
4690 
4691  newsize = calcGrowSize(expr->children[i]->nchildren);
4692  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &childmap, childmapsize, newsize) );
4693  childmapsize = newsize;
4694  }
4695 
4696  /* add children of child i */
4697  SCIP_CALL( exprsimplifyAddChildren(blkmem, expr, expr->children[i]->nchildren, expr->children[i]->children, TRUE, eps, childmap) );
4698 
4699  /* put polynomial of child i into every monomial where child i is used */
4700  j = 0;
4701  while( j < polynomialdata->nmonomials )
4702  {
4703  int factorpos;
4704  SCIP_Bool success;
4705 
4706  monomial = polynomialdata->monomials[j];
4707  /* if monomial is not sorted, then polynomial should not be sorted either, or have only one monomial */
4708  assert(monomial->sorted || !polynomialdata->sorted || polynomialdata->nmonomials <= 1);
4709 
4710  if( SCIPexprFindMonomialFactor(monomial, i, &factorpos) )
4711  {
4712  assert(factorpos >= 0);
4713  assert(factorpos < monomial->nfactors);
4714  /* assert that factors have been merged */
4715  assert(factorpos == 0 || monomial->childidxs[factorpos-1] != i);
4716  assert(factorpos == monomial->nfactors-1 || monomial->childidxs[factorpos+1] != i);
4717 
4718  /* SCIPdebugMessage("attempt expanding child %d at monomial %d factor %d\n", i, j, factorpos);
4719  SCIPdebug( SCIPexprPrint(expr, NULL, NULL, NULL) ); SCIPdebugPrintf("\n");
4720  SCIPdebug( SCIPexprPrint(expr->children[i], NULL, NULL, NULL) ); SCIPdebugPrintf("\n"); */
4721 
4722  SCIP_CALL( polynomialdataExpandMonomialFactor(blkmem, messagehdlr, polynomialdata, j, factorpos,
4723  (SCIP_EXPRDATA_POLYNOMIAL*)expr->children[i]->data.data, childmap, maxexpansionexponent, &success) );
4724 
4725  if( !success )
4726  {
4727  removechild = FALSE;
4728  ++j;
4729  }
4730  }
4731  else
4732  ++j;
4733 
4734  /* expansion may remove monomials[j], move a monomial from the end to position j, or add new monomials to the end of polynomialdata
4735  * we thus repeat with index j, if a factor was successfully expanded
4736  */
4737  }
4738 
4739  /* forget about child i, if it is not used anymore */
4740  if( removechild )
4741  SCIPexprFreeDeep(blkmem, &expr->children[i]);
4742 
4743  /* simplify current polynomial again */
4744  SCIPexprMergeMonomials(blkmem, expr, eps, TRUE);
4745  }
4746 
4747  BMSfreeBlockMemoryArrayNull(blkmem, &childmap, childmapsize);
4748 
4749  /* free children that are not in use anymore */
4751 
4752  /* remove NULLs from children array */
4754 
4755  /* if no children left, then it's a constant polynomial -> change into EXPR_CONST */
4756  if( expr->nchildren == 0 )
4757  {
4758  SCIP_Real val;
4759 
4760  /* if no children, then it should also have no monomials */
4761  assert(polynomialdata->nmonomials == 0);
4762 
4763  val = polynomialdata->constant;
4764  polynomialdataFree(blkmem, &polynomialdata);
4765 
4766  expr->op = SCIP_EXPR_CONST;
4767  expr->data.dbl = val;
4768  }
4769 
4770  SCIPdebugMessage("-> ");
4771  SCIPdebug( SCIPexprPrint(expr, messagehdlr, NULL, NULL, NULL, NULL) );
4772  SCIPdebugPrintf("\n");
4773 
4774  break;
4775  }
4776 
4777  case SCIP_EXPR_LAST:
4778  break;
4779  } /*lint !e788*/
4780 
4781  return SCIP_OKAY;
4782 }
4783 
4784 /** separates linear monomials from an expression, if it is a polynomial expression
4785  *
4786  * Separates only those linear terms whose variable is not used otherwise in the expression.
4787  */
4788 static
4790  BMS_BLKMEM* blkmem, /**< block memory data structure */
4791  SCIP_EXPR* expr, /**< expression */
4792  SCIP_Real eps, /**< threshold, under which positive values are treat as 0 */
4793  int nvars, /**< number of variables in expression */
4794  int* nlinvars, /**< buffer to store number of linear variables in linear part */
4795  int* linidxs, /**< array to store indices of variables in expression tree which belong to linear part */
4796  SCIP_Real* lincoefs /**< array to store coefficients of linear part */
4797  )
4798 {
4799  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
4800  SCIP_EXPRDATA_MONOMIAL* monomial;
4801  int* varsusage;
4802  int* childusage;
4803  int childidx;
4804  int i;
4805  int j;
4806 
4807  assert(blkmem != NULL);
4808  assert(expr != NULL);
4809  assert(nlinvars != NULL);
4810  assert(linidxs != NULL);
4811  assert(lincoefs != NULL);
4812 
4813  *nlinvars = 0;
4814 
4816  return SCIP_OKAY;
4817 
4818  if( SCIPexprGetNChildren(expr) == 0 )
4819  return SCIP_OKAY;
4820 
4821  polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)expr->data.data;
4822  assert(polynomialdata != NULL);
4823 
4824  /* get variable usage */
4825  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &varsusage, nvars) );
4826  BMSclearMemoryArray(varsusage, nvars); /*lint !e644*/
4827  SCIPexprGetVarsUsage(expr, varsusage);
4828 
4829  /* get child usage: how often each child is used in the polynomial */
4830  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &childusage, expr->nchildren) );
4831  BMSclearMemoryArray(childusage, expr->nchildren); /*lint !e644*/
4832  for( i = 0; i < polynomialdata->nmonomials; ++i )
4833  {
4834  monomial = polynomialdata->monomials[i];
4835  assert(monomial != NULL);
4836  for( j = 0; j < monomial->nfactors; ++j )
4837  {
4838  assert(monomial->childidxs[j] >= 0);
4839  assert(monomial->childidxs[j] < expr->nchildren);
4840  ++childusage[monomial->childidxs[j]];
4841  }
4842  }
4843 
4844  /* move linear monomials out of polynomial */
4845  for( i = 0; i < polynomialdata->nmonomials; ++i )
4846  {
4847  monomial = polynomialdata->monomials[i];
4848  assert(monomial != NULL);
4849  if( monomial->nfactors != 1 )
4850  continue;
4851  if( monomial->exponents[0] != 1.0 )
4852  continue;
4853  childidx = monomial->childidxs[0];
4854  if( SCIPexprGetOperator(expr->children[childidx]) != SCIP_EXPR_VARIDX )
4855  continue;
4856 
4857  /* we are at a linear monomial in a variable */
4858  assert(SCIPexprGetOpIndex(expr->children[childidx]) < nvars);
4859  if( childusage[childidx] == 1 && varsusage[SCIPexprGetOpIndex(expr->children[childidx])] == 1 )
4860  {
4861  /* if the child expression is not used in another monomial (which would due to merging be not linear)
4862  * and if the variable is not used somewhere else in the tree,
4863  * then move this monomial into linear part and free child
4864  */
4865  linidxs[*nlinvars] = SCIPexprGetOpIndex(expr->children[childidx]);
4866  lincoefs[*nlinvars] = monomial->coef;
4867  ++*nlinvars;
4868 
4869  SCIPexprFreeDeep(blkmem, &expr->children[childidx]);
4870  monomial->coef = 0.0;
4871  monomial->nfactors = 0;
4872  }
4873  }
4874 
4875  BMSfreeBlockMemoryArray(blkmem, &varsusage, nvars);
4876  BMSfreeBlockMemoryArray(blkmem, &childusage, expr->nchildren);
4877 
4878  if( *nlinvars > 0 )
4879  {
4880  /* if we did something, cleanup polynomial (e.g., remove monomials with coefficient 0.0) */
4881  polynomialdataMergeMonomials(blkmem, polynomialdata, eps, FALSE);
4883  }
4884 
4885  return SCIP_OKAY;
4886 }
4887 
4888 /** converts polynomial expressions back into simpler expressions, where possible */
4889 static
4891  BMS_BLKMEM* blkmem, /**< block memory data structure */
4892  SCIP_EXPR* expr /**< expression to convert back */
4893  )
4894 {
4895  int i;
4896 
4897  assert(blkmem != NULL);
4898  assert(expr != NULL);
4899 
4900  for( i = 0; i < expr->nchildren; ++i )
4901  {
4903  }
4904 
4905  if( expr->op != SCIP_EXPR_POLYNOMIAL )
4906  return SCIP_OKAY;
4907 
4908  SCIP_CALL( exprUnconvertPolynomial(blkmem, &expr->op, &expr->data, expr->nchildren, (void**)expr->children) );
4909 
4910  return SCIP_OKAY;
4911 }
4912 
4913 static
4914 SCIP_DECL_HASHGETKEY( exprparseVarTableGetKey )
4915 { /*lint --e{715}*/
4916  return (void*)((char*)elem + sizeof(int));
4917 }
4918 
4919 /** parses a variable name from a string and creates corresponding expression
4920  *
4921  * Creates a new variable index if variable not seen before, updates varnames and vartable structures.
4922  */
4923 static
4925  BMS_BLKMEM* blkmem, /**< block memory data structure */
4926  const char** str, /**< pointer to the string to be parsed */
4927  SCIP_EXPR** expr, /**< buffer to store pointer to created expression */
4928  int* nvars, /**< running number of encountered variables so far */
4929  int** varnames, /**< pointer to buffer to store new variable names */
4930  int* varnameslength, /**< pointer to length of the varnames buffer array */
4931  SCIP_HASHTABLE* vartable, /**< hash table for variable names and corresponding expression index */
4932  SCIP_Real coefficient, /**< coefficient to be used when creating the expression */
4933  const char* varnameendptr /**< if a <varname> should be parsed, set this to NULL. Then, str points to the '<'
4934  else, str should point to the first letter of the varname, and varnameendptr should
4935  point one char behind the last char of the variable name */
4936  )
4937 {
4938  int namelength;
4939  int varidx;
4940  char varname[SCIP_MAXSTRLEN];
4941  void* element;
4942 
4943  assert(blkmem != NULL);
4944  assert(str != NULL);
4945  assert(expr != NULL);
4946  assert(nvars != NULL);
4947  assert(varnames != NULL);
4948  assert(vartable != NULL);
4949 
4950  if( varnameendptr == NULL )
4951  {
4952  ++*str;
4953  varnameendptr = *str;
4954  while( varnameendptr[0] != '>' )
4955  ++varnameendptr;
4956  }
4957 
4958  namelength = varnameendptr - *str; /*lint !e712*/
4959  if( namelength >= SCIP_MAXSTRLEN )
4960  {
4961  SCIPerrorMessage("Variable name %.*s is too long for buffer in exprparseReadVariable.\n", namelength, str);
4962  return SCIP_READERROR;
4963  }
4964 
4965  memcpy(varname, *str, namelength * sizeof(char));
4966  varname[namelength] = '\0';
4967 
4968  element = SCIPhashtableRetrieve(vartable, varname);
4969  if( element != NULL )
4970  {
4971  /* variable is old friend */
4972  assert(strcmp((char*)element + sizeof(int), varname) == 0);
4973 
4974  varidx = *(int*)element;
4975  }
4976  else
4977  {
4978  /* variable is new */
4979  varidx = *nvars;
4980 
4981  (*varnameslength) -= (int)(1 + (strlen(varname) + 1) / sizeof(int) + 1);
4982  if( *varnameslength < 0 )
4983  {
4984  SCIPerrorMessage("Buffer in exprparseReadVariable is too short for varaible name %.*s.\n", namelength, str);
4985  return SCIP_READERROR;
4986  }
4987 
4988  /* store index of variable and variable name in varnames buffer */
4989  **varnames = varidx;
4990  (void) SCIPstrncpy((char*)(*varnames + 1), varname, (int)strlen(varname)+1);
4991 
4992  /* insert variable into hashtable */
4993  SCIP_CALL( SCIPhashtableInsert(vartable, (void*)*varnames) );
4994 
4995  ++*nvars;
4996  *varnames += 1 + (strlen(varname) + 1) / sizeof(int) + 1;
4997  }
4998 
4999  /* create VARIDX expression, put into LINEAR expression if we have coefficient != 1 */
5000  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_VARIDX, varidx) ); /*lint !e613*/
5001  if( coefficient != 1.0 )
5002  {
5003  SCIP_CALL( SCIPexprCreateLinear(blkmem, expr, 1, expr, &coefficient, 0.0) );
5004  }
5005 
5006  /* Move pointer to char behind end of variable */
5007  *str = varnameendptr + 1;
5008 
5009  /* consprint sometimes prints a variable type identifier which we don't need */
5010  if( (*str)[0] == '[' && (*str)[2] == ']' &&
5011  ((*str)[1] == SCIP_VARTYPE_BINARY_CHAR ||
5012  (*str)[1] == SCIP_VARTYPE_INTEGER_CHAR ||
5013  (*str)[1] == SCIP_VARTYPE_IMPLINT_CHAR ||
5014  (*str)[1] == SCIP_VARTYPE_CONTINUOUS_CHAR ) )
5015  *str += 3;
5016 
5017  return SCIP_OKAY;
5018 }
5019 
5020 /** if str[0] points to an opening parenthesis, this function sets endptr to point to the matching closing bracket in str
5021  *
5022  * Searches for at most length characters.
5023  */
5024 static
5026  const char* str, /**< pointer to the string to be parsed */
5027  const char** endptr, /**< pointer to point to the closing parenthesis */
5028  int length /**< length of the string to be parsed */
5029  )
5030 {
5031  int nopenbrackets;
5032 
5033  assert(str[0] == '(');
5034 
5035  *endptr = str;
5036 
5037  /* find the end of this expression */
5038  nopenbrackets = 0;
5039  while( (*endptr - str ) < length && !(nopenbrackets == 1 && *endptr[0] == ')') )
5040  {
5041  if( *endptr[0] == '(')
5042  ++nopenbrackets;
5043  if( *endptr[0] == ')')
5044  --nopenbrackets;
5045  ++*endptr;
5046  }
5047 
5048  if( *endptr[0] != ')' )
5049  {
5050  SCIPerrorMessage("unable to find closing parenthesis in unbalanced expression %.*s\n", length, str);
5051  return SCIP_READERROR;
5052  }
5053 
5054  return SCIP_OKAY;
5055 }
5056 
5057 /** this function sets endptr to point to the next separating comma in str
5058  *
5059  * That is, for a given string like "x+f(x,y),z", endptr will point to the comma before "z"
5060  *
5061  * Searches for at most length characters.
5062  */
5063 static
5065  const char* str, /**< pointer to the string to be parsed */
5066  const char** endptr, /**< pointer to point to the comma */
5067  int length /**< length of the string to be parsed */
5068  )
5069 {
5070  int nopenbrackets;
5071 
5072  *endptr = str;
5073 
5074  /* find a comma without open brackets */
5075  nopenbrackets = 0;
5076  while( (*endptr - str ) < length && !(nopenbrackets == 0 && *endptr[0] == ',') )
5077  {
5078  if( *endptr[0] == '(')
5079  ++nopenbrackets;
5080  if( *endptr[0] == ')')
5081  --nopenbrackets;
5082  ++*endptr;
5083  }
5084 
5085  if( *endptr[0] != ',' )
5086  {
5087  SCIPerrorMessage("unable to find separating comma in unbalanced expression %.*s\n", length, str);
5088  return SCIP_READERROR;
5089  }
5090 
5091  return SCIP_OKAY;
5092 }
5093 
5094 /** parses an expression from a string */
5095 static
5097  BMS_BLKMEM* blkmem, /**< block memory data structure */
5098  SCIP_MESSAGEHDLR* messagehdlr, /**< message handler */
5099  SCIP_EXPR** expr, /**< buffer to store pointer to created expression */
5100  const char* str, /**< pointer to the string to be parsed */
5101  int length, /**< length of the string to be parsed */
5102  const char* lastchar, /**< pointer to the last char of str that should be parsed */
5103  int* nvars, /**< running number of encountered variables so far */
5104  int** varnames, /**< pointer to buffer to store new variable names */
5105  int* varnameslength, /**< pointer to length of the varnames buffer array */
5106  SCIP_HASHTABLE* vartable, /**< hash table for variable names and corresponding expression index */
5107  int recursiondepth /**< current recursion depth */
5108  )
5109 { /*lint --e{712,747}*/
5110  SCIP_EXPR* arg1;
5111  SCIP_EXPR* arg2;
5112  const char* subexpptr;
5113  const char* subexpendptr;
5114  const char* strstart;
5115  const char* endptr;
5116  char* nonconstendptr;
5117  SCIP_Real number;
5118  int subexplength;
5119  int nopenbrackets;
5120 
5121  assert(blkmem != NULL);
5122  assert(expr != NULL);
5123  assert(str != NULL);
5124  assert(lastchar >= str);
5125  assert(nvars != NULL);
5126  assert(varnames != NULL);
5127  assert(vartable != NULL);
5128 
5129  assert(recursiondepth < 100);
5130 
5131  strstart = str; /* might be needed for error message... */
5132 
5133  SCIPdebugMessage("exprParse (%i): parsing %.*s\n", recursiondepth, (int) (lastchar-str + 1), str);
5134 
5135  /* ignore whitespace */
5136  while( isspace((unsigned char)*str) )
5137  ++str;
5138 
5139  /* look for a sum or difference not contained in brackets */
5140  subexpptr = str;
5141  nopenbrackets = 0;
5142 
5143  /* find the end of this expression
5144  * a '+' right at the beginning indicates a coefficient, not treated here, or a summation
5145  */
5146  while( subexpptr != lastchar && !(nopenbrackets == 0 && (subexpptr[0] == '+' || subexpptr[0] == '-') && subexpptr != str) )
5147  {
5148  if( subexpptr[0] == '(')
5149  ++nopenbrackets;
5150  if( subexpptr[0] == ')')
5151  --nopenbrackets;
5152  ++subexpptr;
5153  }
5154 
5155  if( subexpptr != lastchar )
5156  {
5157  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str, (int) ((subexpptr - 1) - str + 1), subexpptr - 1, nvars,
5158  varnames, varnameslength, vartable, recursiondepth + 1) );
5159 
5160  if( subexpptr[0] == '+' )
5161  ++subexpptr;
5162  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, subexpptr , (int) (lastchar - (subexpptr ) + 1), lastchar, nvars,
5163  varnames, varnameslength, vartable, recursiondepth + 1) );
5164 
5165  /* make new expression from two arguments
5166  * we always use add, because we leave the operator between the found expressions in the second argument
5167  * this way, we do not have to worry about ''minus brackets'' in the case of more then two summands:
5168  * a - b - c = a + (-b -c)
5169  */
5170  SCIP_CALL( SCIPexprAdd(blkmem, expr, 1.0, arg1, 1.0, arg2, 0.0) );
5171 
5172  SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
5173  SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
5174  SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );
5175 
5176  return SCIP_OKAY;
5177  }
5178 
5179  /* check for a bracketed subexpression */
5180  if( str[0] == '(' )
5181  {
5182  nopenbrackets = 0;
5183 
5184  subexplength = -1; /* we do not want the closing bracket in the string */
5185  subexpptr = str + 1; /* leave out opening bracket */
5186 
5187  /* find the end of this expression */
5188  while( subexplength < length && !(nopenbrackets == 1 && str[0] == ')') )
5189  {
5190  if( str[0] == '(' )
5191  ++nopenbrackets;
5192  if( str[0] == ')' )
5193  --nopenbrackets;
5194  ++str;
5195  ++subexplength;
5196  }
5197  subexpendptr = str - 1; /* leave out closing bracket */
5198 
5199  SCIP_CALL( exprParse(blkmem, messagehdlr, expr, subexpptr, subexplength, subexpendptr, nvars, varnames,
5200  varnameslength, vartable, recursiondepth + 1) );
5201  ++str;
5202  }
5203  else if( isdigit((unsigned char)str[0]) || ((str[0] == '-' || str[0] == '+')
5204  && (isdigit((unsigned char)str[1]) || str[1] == ' ')) )
5205  {
5206  /* check if there is a lonely minus coming, indicating a -1.0 */
5207  if( str[0] == '-' && str[1] == ' ' )
5208  {
5209  number = -1.0;
5210  nonconstendptr = (char*) str + 1;
5211  }
5212  /* check if there is a number coming */
5213  else if( !SCIPstrToRealValue(str, &number, &nonconstendptr) )
5214  {
5215  SCIPerrorMessage("error parsing number from <%s>\n", str);
5216  return SCIP_READERROR;
5217  }
5218  str = nonconstendptr;
5219 
5220  /* ignore whitespace */
5221  while( isspace((unsigned char)*str) && str != lastchar )
5222  ++str;
5223 
5224  if( str[0] != '*' && str[0] != '/' && str[0] != '+' && str[0] != '-' && str[0] != '^' )
5225  {
5226  if( str < lastchar )
5227  {
5228  SCIP_CALL( exprParse(blkmem, messagehdlr, expr, str, (int)(lastchar - str) + 1, lastchar, nvars, varnames,
5229  varnameslength, vartable, recursiondepth + 1) );
5230  SCIP_CALL( SCIPexprMulConstant(blkmem, expr, *expr, number) );
5231  }
5232  else
5233  {
5234  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_CONST, number) );
5235  }
5236  str = lastchar + 1;
5237  }
5238  else
5239  {
5240  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_CONST, number) );
5241  }
5242  }
5243  else if( str[0] == '<' )
5244  {
5245  /* check if expressions begins with a variable */
5246  SCIP_CALL( exprparseReadVariable(blkmem, &str, expr, nvars, varnames, varnameslength, vartable, 1.0, NULL) );
5247  }
5248  /* four character operators */
5249  else if( strncmp(str, "sqrt", 4) == 0 )
5250  {
5251  str += 4;
5252  SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
5253  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, endptr - str - 1, endptr -1, nvars, varnames,
5254  varnameslength, vartable, recursiondepth + 1) );
5255  str = endptr + 1;
5256 
5257  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_SQRT, arg1) );
5258  }
5259  /* three character operators with 1 argument */
5260  else if(
5261  strncmp(str, "abs", 3) == 0 ||
5262  strncmp(str, "cos", 3) == 0 ||
5263  strncmp(str, "exp", 3) == 0 ||
5264  strncmp(str, "log", 3) == 0 ||
5265  strncmp(str, "sin", 3) == 0 ||
5266  strncmp(str, "sqr", 3) == 0 ||
5267  strncmp(str, "tan", 3) == 0 )
5268  {
5269  const char* opname = str;
5270 
5271  str += 3;
5272  SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
5273  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, endptr - str - 1, endptr -1, nvars, varnames,
5274  varnameslength, vartable, recursiondepth + 1) );
5275  str = endptr + 1;
5276 
5277  if( strncmp(opname, "abs", 3) == 0 )
5278  {
5279  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_ABS, arg1) );
5280  }
5281  else if( strncmp(opname, "cos", 3) == 0 )
5282  {
5283  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_COS, arg1) );
5284  }
5285  else if( strncmp(opname, "exp", 3) == 0 )
5286  {
5287  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_EXP, arg1) );
5288  }
5289  else if( strncmp(opname, "log", 3) == 0 )
5290  {
5291  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_LOG, arg1) );
5292  }
5293  else if( strncmp(opname, "sin", 3) == 0 )
5294  {
5295  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_SIN, arg1) );
5296  }
5297  else if( strncmp(opname, "sqr", 3) == 0 )
5298  {
5299  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_SQUARE, arg1) );
5300  }
5301  else
5302  {
5303  assert(strncmp(opname, "tan", 3) == 0);
5304  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_TAN, arg1) );
5305  }
5306  }
5307  /* three character operators with 2 arguments */
5308  else if(
5309  strncmp(str, "max", 3) == 0 ||
5310  strncmp(str, "min", 3) == 0 )
5311  {
5312  /* we have a string of the form "min(...,...)" or "max(...,...)"
5313  * first find the closing parenthesis, then the comma
5314  */
5315  const char* comma;
5316  SCIP_EXPROP op;
5317 
5318  op = (str[1] == 'a' ? SCIP_EXPR_MAX : SCIP_EXPR_MIN);
5319 
5320  str += 3;
5321  SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
5322 
5323  SCIP_CALL( exprparseFindSeparatingComma(str+1, &comma, endptr - str - 1) );
5324 
5325  /* parse first argument [str+1..comma-1] */
5326  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, comma - str - 1, comma - 1, nvars, varnames,
5327  varnameslength, vartable, recursiondepth + 1) );
5328 
5329  /* parse second argument [comma+1..endptr] */
5330  ++comma;
5331  while( comma < endptr && *comma == ' ' )
5332  ++comma;
5333 
5334  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, comma, endptr - comma, endptr - 1, nvars, varnames,
5335  varnameslength, vartable, recursiondepth + 1) );
5336 
5337  SCIP_CALL( SCIPexprCreate(blkmem, expr, op, arg1, arg2) );
5338 
5339  str = endptr + 1;
5340  }
5341  else if( strncmp(str, "power", 5) == 0 )
5342  {
5343  /* we have a string of the form "power(...,integer)" (thus, intpower)
5344  * first find the closing parenthesis, then the comma
5345  */
5346  const char* comma;
5347  int exponent;
5348 
5349  str += 5;
5350  SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
5351 
5352  SCIP_CALL( exprparseFindSeparatingComma(str+1, &comma, endptr - str - 1) );
5353 
5354  /* parse first argument [str+1..comma-1] */
5355  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, comma - str - 1, comma - 1, nvars, varnames,
5356  varnameslength, vartable, recursiondepth + 1) );
5357 
5358  ++comma;
5359  /* parse second argument [comma, endptr-1]: it needs to be an integer */
5360  while( comma < endptr && *comma == ' ' )
5361  ++comma;
5362  if( !isdigit((unsigned char)comma[0]) && !((comma[0] == '-' || comma[0] == '+') && isdigit((unsigned char)comma[1])) )
5363  {
5364  SCIPerrorMessage("error parsing integer exponent from <%s>\n", comma);
5365  }
5366  if( !SCIPstrToIntValue(comma, &exponent, &nonconstendptr) )
5367  {
5368  SCIPerrorMessage("error parsing integer from <%s>\n", comma);
5369  return SCIP_READERROR;
5370  }
5371 
5372  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_INTPOWER, arg1, exponent) );
5373 
5374  str = endptr + 1;
5375  }
5376  else if( strncmp(str, "realpower", 9) == 0 || strncmp(str, "signpower", 9) == 0 )
5377  {
5378  /* we have a string of the form "realpower(...,double)" or "signpower(...,double)"
5379  * first find the closing parenthesis, then the comma
5380  */
5381  const char* opname = str;
5382  const char* comma;
5383 
5384  str += 9;
5385  SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
5386 
5387  SCIP_CALL( exprparseFindSeparatingComma(str+1, &comma, endptr - str - 1) );
5388 
5389  /* parse first argument [str+1..comma-1] */
5390  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, comma - str - 1, comma - 1, nvars, varnames,
5391  varnameslength, vartable, recursiondepth + 1) );
5392 
5393  ++comma;
5394  /* parse second argument [comma, endptr-1]: it needs to be an number */
5395  while( comma < endptr && *comma == ' ' )
5396  ++comma;
5397  if( !isdigit((unsigned char)comma[0]) && !((comma[0] == '-' || comma[0] == '+') && isdigit((unsigned char)comma[1])) )
5398  {
5399  SCIPerrorMessage("error parsing number exponent from <%s>\n", comma);
5400  }
5401  if( !SCIPstrToRealValue(comma, &number, &nonconstendptr) )
5402  {
5403  SCIPerrorMessage("error parsing number from <%s>\n", comma);
5404  return SCIP_READERROR;
5405  }
5406 
5407  if( strncmp(opname, "realpower", 9) == 0 )
5408  {
5409  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_REALPOWER, arg1, number) );
5410  }
5411  else
5412  {
5413  assert(strncmp(opname, "signpower", 9) == 0);
5414  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_SIGNPOWER, arg1, number) );
5415  }
5416 
5417  str = endptr + 1;
5418  }
5419  else if( isalpha(*str) || *str == '_' || *str == '#' )
5420  {
5421  /* check for a variable, that was not recognized earlier because somebody omitted the '<' and '>' we need for
5422  * SCIPparseVarName, making everyones life harder;
5423  * we allow only variable names starting with a character or underscore here
5424  */
5425  const char* varnamestartptr = str;
5426 
5427  /* allow only variable names containing characters, digits, hash marks, and underscores here */
5428  while( isalnum(str[0]) || str[0] == '_' || str[0] == '#' )
5429  ++str;
5430 
5431  SCIP_CALL( exprparseReadVariable(blkmem, &varnamestartptr, expr, nvars, varnames, varnameslength,
5432  vartable, 1.0, str) );
5433  }
5434  else
5435  {
5436  SCIPerrorMessage("parsing of invalid expression %.*s.\n", (int) (lastchar - str + 1), str);
5437  return SCIP_READERROR;
5438  }
5439 
5440  /* if we are one char behind lastchar, we are done */
5441  if( str == lastchar + 1)
5442  {
5443  SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
5444  SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
5445  SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );
5446 
5447  return SCIP_OKAY;
5448  }
5449 
5450  /* check if we are still in bounds */
5451  if( str > lastchar + 1)
5452  {
5453  SCIPerrorMessage("error finding first expression in \"%.*s\" took us outside of given subexpression length\n", length, strstart);
5454  return SCIP_READERROR;
5455  }
5456 
5457  /* ignore whitespace */
5458  while( isspace((unsigned char)*str) && str != lastchar + 1 )
5459  ++str;
5460 
5461  /* maybe now we're done? */
5462  if( str >= lastchar + 1)
5463  {
5464  SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
5465  SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
5466  SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );
5467 
5468  return SCIP_OKAY;
5469  }
5470 
5471  if( str[0] == '^' )
5472  {
5473  /* a '^' behind the found expression indicates a power */
5474  SCIP_Real constant;
5475 
5476  arg1 = *expr;
5477  ++str;
5478  while( isspace((unsigned char)*str) && str != lastchar + 1 )
5479  ++str;
5480 
5481  if( isdigit((unsigned char)str[0]) || ((str[0] == '-' || str[0] == '+') && isdigit((unsigned char)str[1])) )
5482  {
5483  /* there is a number coming */
5484  if( !SCIPstrToRealValue(str, &number, &nonconstendptr) )
5485  {
5486  SCIPerrorMessage("error parsing number from <%s>\n", str);
5487  return SCIP_READERROR;
5488  }
5489 
5490  SCIP_CALL( SCIPexprCreate(blkmem, &arg2, SCIP_EXPR_CONST, number) );
5491  str = nonconstendptr;
5492 
5493  constant = SCIPexprGetOpReal(arg2);
5494  SCIPexprFreeDeep(blkmem, &arg2);
5495 
5496  /* expr^number is intpower or realpower */
5497  if( EPSISINT(constant, 0.0) ) /*lint !e835*/
5498  {
5499  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_INTPOWER, arg1, (int)constant) );
5500  }
5501  else
5502  {
5503  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_REALPOWER, arg1, constant) );
5504  }
5505  }
5506  else if( str[0] == '(' )
5507  {
5508  /* we use exprParse to evaluate the exponent */
5509 
5510  SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
5511  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, str + 1, endptr - str - 1, endptr -1, nvars, varnames,
5512  varnameslength, vartable, recursiondepth + 1) );
5513 
5514  if( SCIPexprGetOperator(arg2) != SCIP_EXPR_CONST )
5515  {
5516  /* reformulate arg1^arg2 as exp(arg2 * log(arg1)) */
5517  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_LOG, arg1) );
5518  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_MUL, *expr, arg2) );
5519  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_EXP, *expr) );
5520  }
5521  else
5522  {
5523  /* expr^number is intpower or realpower */
5524  constant = SCIPexprGetOpReal(arg2);
5525  SCIPexprFreeDeep(blkmem, &arg2);
5526  if( EPSISINT(constant, 0.0) ) /*lint !e835*/
5527  {
5528  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_INTPOWER, arg1, (int)constant) );
5529  }
5530  else
5531  {
5532  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_REALPOWER, arg1, constant) );
5533  }
5534  }
5535  str = endptr + 1;
5536  }
5537  else
5538  {
5539  SCIPerrorMessage("unexpected string following ^ in %.*s\n", length, str);
5540  return SCIP_READERROR;
5541  }
5542 
5543  /* ignore whitespace */
5544  while( isspace((unsigned char)*str) && str != lastchar + 1 )
5545  ++str;
5546  }
5547 
5548  /* check for a two argument operator that is not a multiplication */
5549  if( str <= lastchar && (str[0] == '+' || str[0] == '-' || str[0] == '/') )
5550  {
5551  char op;
5552 
5553  op = str[0];
5554  arg1 = *expr;
5555 
5556  /* step forward over the operator to go to the beginning of the second argument */
5557  ++str;
5558 
5559  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, str, (int) (lastchar - str + 1), lastchar, nvars, varnames,
5560  varnameslength, vartable, recursiondepth + 1) );
5561  str = lastchar + 1;
5562 
5563  /* make new expression from two arguments */
5564  if( op == '+')
5565  {
5566  SCIP_CALL( SCIPexprAdd(blkmem, expr, 1.0, arg1, 1.0, arg2, 0.0) );
5567  }
5568  else if( op == '-')
5569  {
5570  SCIP_CALL( SCIPexprAdd(blkmem, expr, 1.0, arg1, -1.0, arg2, 0.0) );
5571  }
5572  else if( op == '*' )
5573  {
5574  if( SCIPexprGetOperator(arg1) == SCIP_EXPR_CONST )
5575  {
5576  SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg2, SCIPexprGetOpReal(arg1)) );
5577  }
5578  else if( SCIPexprGetOperator(arg2) == SCIP_EXPR_CONST )
5579  {
5580  SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg1, SCIPexprGetOpReal(arg2)) );
5581  }
5582  else
5583  {
5584  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_MUL, arg1, arg2) );
5585  }
5586  }
5587  else
5588  {
5589  assert(op == '/');
5590 
5591  if( SCIPexprGetOperator(arg2) == SCIP_EXPR_CONST )
5592  {
5593  SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg1, 1.0 / SCIPexprGetOpReal(arg2)) );
5594  SCIPexprFreeShallow(blkmem, &arg2);
5595  }
5596  else
5597  {
5598  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_DIV, arg1, arg2) );
5599  }
5600  }
5601  }
5602 
5603  /* ignore whitespace */
5604  while( isspace((unsigned char)*str) )
5605  ++str;
5606 
5607  /* we are either done or we have a multiplication? */
5608  if( str >= lastchar + 1 )
5609  {
5610  SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
5611  SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
5612  SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );
5613 
5614  return SCIP_OKAY;
5615  }
5616 
5617  /* if there is a part of the string left to be parsed, we assume that this as a multiplication */
5618  arg1 = *expr;
5619 
5620  /* stepping over multiplication operator if needed */
5621  if( str[0] == '*' )
5622  {
5623  ++str;
5624  }
5625  else if( str[0] != '(' )
5626  {
5627  SCIPdebugMessage("No operator found, assuming a multiplication before %.*s\n", (int) (lastchar - str + 1), str);
5628  }
5629 
5630  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, str, (int) (lastchar - str + 1), lastchar, nvars, varnames,
5631  varnameslength, vartable, recursiondepth + 1) );
5632 
5633  if( SCIPexprGetOperator(arg1) == SCIP_EXPR_CONST )
5634  {
5635  SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg2, SCIPexprGetOpReal(arg1)) );
5636  SCIPexprFreeDeep(blkmem, &arg1);
5637  }
5638  else if( SCIPexprGetOperator(arg2) == SCIP_EXPR_CONST )
5639  {
5640  SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg1, SCIPexprGetOpReal(arg2)) );
5641  SCIPexprFreeDeep(blkmem, &arg2);
5642  }
5643  else
5644  {
5645  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_MUL, arg1, arg2) );
5646  }
5647 
5648  SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
5649  SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
5650  SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );
5651 
5652  return SCIP_OKAY;
5653 }
5654 
5655 /**@} */
5656 
5657 /**@name Expression methods */
5658 /**@{ */
5659 
5660 /* In debug mode, the following methods are implemented as function calls to ensure
5661  * type validity.
5662  * In optimized mode, the methods are implemented as defines to improve performance.
5663  * However, we want to have them in the library anyways, so we have to undef the defines.
5664  */
5665 
5666 #undef SCIPexprGetOperator
5667 #undef SCIPexprGetNChildren
5668 #undef SCIPexprGetChildren
5669 #undef SCIPexprGetOpIndex
5670 #undef SCIPexprGetOpReal
5671 #undef SCIPexprGetOpData
5672 #undef SCIPexprGetRealPowerExponent
5673 #undef SCIPexprGetIntPowerExponent
5674 #undef SCIPexprGetSignPowerExponent
5675 #undef SCIPexprGetLinearCoefs
5676 #undef SCIPexprGetLinearConstant
5677 #undef SCIPexprGetQuadElements
5678 #undef SCIPexprGetQuadConstant
5679 #undef SCIPexprGetQuadLinearCoefs
5680 #undef SCIPexprGetNQuadElements
5681 #undef SCIPexprGetMonomials
5682 #undef SCIPexprGetNMonomials
5683 #undef SCIPexprGetPolynomialConstant
5684 #undef SCIPexprGetMonomialCoef
5685 #undef SCIPexprGetMonomialNFactors
5686 #undef SCIPexprGetMonomialChildIndices
5687 #undef SCIPexprGetMonomialExponents
5688 #undef SCIPexprGetUserData
5689 #undef SCIPexprHasUserEstimator
5690 #undef SCIPexprGetUserEvalCapability
5691 
5692 /** gives operator of expression */
5694  SCIP_EXPR* expr /**< expression */
5695  )
5696 {
5697  assert(expr != NULL);
5698 
5699  return expr->op;
5700 }
5701 
5702 /** gives number of children of an expression */
5704  SCIP_EXPR* expr /**< expression */
5705  )
5706 {
5707  assert(expr != NULL);
5708 
5709  return expr->nchildren;
5710 }
5711 
5712 /** gives pointer to array with children of an expression */
5714  SCIP_EXPR* expr /**< expression */
5715  )
5716 {
5717  assert(expr != NULL);
5718 
5719  return expr->children;
5720 }
5721 
5722 /** gives index belonging to a SCIP_EXPR_VARIDX or SCIP_EXPR_PARAM operand */
5724  SCIP_EXPR* expr /**< expression */
5725  )
5726 {
5727  assert(expr != NULL);
5728  assert(expr->op == SCIP_EXPR_VARIDX || expr->op == SCIP_EXPR_PARAM);
5729 
5730  return expr->data.intval;
5731 }
5732 
5733 /** gives real belonging to a SCIP_EXPR_CONST operand */
5735  SCIP_EXPR* expr /**< expression */
5736  )
5737 {
5738  assert(expr != NULL);
5739  assert(expr->op == SCIP_EXPR_CONST);
5740 
5741  return expr->data.dbl;
5742 }
5743 
5744 /** gives void* belonging to a complex operand */
5746  SCIP_EXPR* expr /**< expression */
5747  )
5748 {
5749  assert(expr != NULL);
5750  assert(expr->op >= SCIP_EXPR_SUM); /* only complex operands store their data as void* */
5751 
5752  return expr->data.data;
5753 }
5754 
5755 /** gives exponent belonging to a SCIP_EXPR_REALPOWER expression */
5757  SCIP_EXPR* expr /**< expression */
5758  )
5759 {
5760  assert(expr != NULL);
5761  assert(expr->op == SCIP_EXPR_REALPOWER);
5762 
5763  return expr->data.dbl;
5764 }
5765 
5766 /** gives exponent belonging to a SCIP_EXPR_INTPOWER expression */
5768  SCIP_EXPR* expr /**< expression */
5769  )
5770 {
5771  assert(expr != NULL);
5772  assert(expr->op == SCIP_EXPR_INTPOWER);
5773 
5774  return expr->data.intval;
5775 }
5776 
5777 /** gives exponent belonging to a SCIP_EXPR_SIGNPOWER expression */
5779  SCIP_EXPR* expr /**< expression */
5780  )
5781 {
5782  assert(expr != NULL);
5783  assert(expr->op == SCIP_EXPR_SIGNPOWER);
5784 
5785  return expr->data.dbl;
5786 }
5787 
5788 /** gives linear coefficients belonging to a SCIP_EXPR_LINEAR expression */
5790  SCIP_EXPR* expr /**< expression */
5791  )
5792 {
5793  assert(expr != NULL);
5794  assert(expr->op == SCIP_EXPR_LINEAR);
5795  assert(expr->data.data != NULL);
5796 
5797  /* the coefficients are stored in the first nchildren elements of the array stored as expression data */
5798  return (SCIP_Real*)expr->data.data;
5799 }
5800 
5801 /** gives constant belonging to a SCIP_EXPR_LINEAR expression */
5803  SCIP_EXPR* expr /**< expression */
5804  )
5805 {
5806  assert(expr != NULL);
5807  assert(expr->op == SCIP_EXPR_LINEAR);
5808  assert(expr->data.data != NULL);
5809 
5810  /* the constant is stored in the nchildren's element of the array stored as expression data */
5811  return ((SCIP_Real*)expr->data.data)[expr->nchildren];
5812 }
5813 
5814 /** gives quadratic elements belonging to a SCIP_EXPR_QUADRATIC expression */
5816  SCIP_EXPR* expr /**< quadratic expression */
5817  )
5818 {
5819  assert(expr != NULL);
5820  assert(expr->op == SCIP_EXPR_QUADRATIC);
5821  assert(expr->data.data != NULL);
5822 
5823  return ((SCIP_EXPRDATA_QUADRATIC*)expr->data.data)->quadelems;
5824 }
5825 
5826 /** gives constant belonging to a SCIP_EXPR_QUADRATIC expression */
5828  SCIP_EXPR* expr /**< quadratic expression */
5829  )
5830 {
5831  assert(expr != NULL);
5832  assert(expr->op == SCIP_EXPR_QUADRATIC);
5833  assert(expr->data.data != NULL);
5834 
5835  return ((SCIP_EXPRDATA_QUADRATIC*)expr->data.data)->constant;
5836 }
5837 
5838 /** gives linear coefficients belonging to a SCIP_EXPR_QUADRATIC expression
5839  * can be NULL if all coefficients are 0.0 */
5841  SCIP_EXPR* expr /**< quadratic expression */
5842  )
5843 {
5844  assert(expr != NULL);
5845  assert(expr->op == SCIP_EXPR_QUADRATIC);
5846  assert(expr->data.data != NULL);
5847 
5848  return ((SCIP_EXPRDATA_QUADRATIC*)expr->data.data)->lincoefs;
5849 }
5850 
5851 /** gives number of quadratic elements belonging to a SCIP_EXPR_QUADRATIC expression */
5853  SCIP_EXPR* expr /**< quadratic expression */
5854  )
5855 {
5856  assert(expr != NULL);
5857  assert(expr->op == SCIP_EXPR_QUADRATIC);
5858  assert(expr->data.data != NULL);
5859 
5860  return ((SCIP_EXPRDATA_QUADRATIC*)expr->data.data)->nquadelems;
5861 }
5862 
5863 /** gives the monomials belonging to a SCIP_EXPR_POLYNOMIAL expression */
5865  SCIP_EXPR* expr /**< expression */
5866  )
5867 {
5868  assert(expr != NULL);
5869  assert(expr->op == SCIP_EXPR_POLYNOMIAL);
5870  assert(expr->data.data != NULL);
5871 
5872  return ((SC