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