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