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-2020 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  quaddata->lincoefs[polynomialdata->monomials[i]->childidxs[0]] += polynomialdata->monomials[i]->coef;
3976  }
3977  else
3978  {
3979  /* monomial should be a square term */
3980  assert(polynomialdata->monomials[i]->exponents[0] == 2.0);
3981  assert(quadelemidx < quaddata->nquadelems);
3982  quaddata->quadelems[quadelemidx].idx1 = polynomialdata->monomials[i]->childidxs[0];
3983  quaddata->quadelems[quadelemidx].idx2 = polynomialdata->monomials[i]->childidxs[0];
3984  quaddata->quadelems[quadelemidx].coef = polynomialdata->monomials[i]->coef;
3985  ++quadelemidx;
3986  }
3987  }
3988  else
3989  {
3990  /* monomial should be a bilinear term */
3991  assert(polynomialdata->monomials[i]->exponents[0] == 1.0);
3992  assert(polynomialdata->monomials[i]->exponents[1] == 1.0);
3993  assert(quadelemidx < quaddata->nquadelems);
3994  quaddata->quadelems[quadelemidx].idx1 = MIN(polynomialdata->monomials[i]->childidxs[0], polynomialdata->monomials[i]->childidxs[1]);
3995  quaddata->quadelems[quadelemidx].idx2 = MAX(polynomialdata->monomials[i]->childidxs[0], polynomialdata->monomials[i]->childidxs[1]);
3996  quaddata->quadelems[quadelemidx].coef = polynomialdata->monomials[i]->coef;
3997  ++quadelemidx;
3998  }
3999  }
4000  assert(quadelemidx == quaddata->nquadelems);
4001 
4002  polynomialdataFree(blkmem, &polynomialdata);
4003 
4004  *op = SCIP_EXPR_QUADRATIC;
4005  data->data = (void*)quaddata;
4006 
4007  return SCIP_OKAY;
4008  }
4009 
4010  if( polynomialdata->constant == 0.0 && polynomialdata->nmonomials == 1 && polynomialdata->monomials[0]->coef == 1.0 )
4011  {
4012  /* polynomial is product of children */
4013  monomial = polynomialdata->monomials[0];
4014  assert(monomial->nfactors == nchildren);
4015 
4016  if( monomial->nfactors == 1 )
4017  {
4018  /* polynomial is x^k for some k */
4019  assert(monomial->exponents[0] != 1.0); /* should have been handled before */
4020  assert(monomial->childidxs[0] == 0);
4021 
4022  if( monomial->exponents[0] == 2.0 )
4023  {
4024  /* polynomial is x^2, so turn into SCIP_EXPR_SQUARE */
4025 
4026  polynomialdataFree(blkmem, &polynomialdata);
4027  data->data = NULL;
4028 
4029  *op = SCIP_EXPR_SQUARE;
4030 
4031  return SCIP_OKAY;
4032  }
4033 
4034  if( EPSISINT(monomial->exponents[0], 0.0) ) /*lint !e835*/
4035  {
4036  /* k is an integer, so turn into SCIP_EXPR_INTPOWER */
4037  int exponent;
4038 
4039  exponent = (int)EPSROUND(monomial->exponents[0], 0.0); /*lint !e835*/
4040 
4041  polynomialdataFree(blkmem, &polynomialdata);
4042 
4043  *op = SCIP_EXPR_INTPOWER;
4044  data->intval = exponent;
4045 
4046  return SCIP_OKAY;
4047  }
4048 
4049  if( monomial->exponents[0] == 0.5 )
4050  {
4051  /* polynomial is sqrt(x), so turn into SCIP_EXPR_SQRT */
4052 
4053  polynomialdataFree(blkmem, &polynomialdata);
4054  data->data = NULL;
4055 
4056  *op = SCIP_EXPR_SQRT;
4057 
4058  return SCIP_OKAY;
4059  }
4060 
4061  {
4062  /* polynomial is x^a with a some real number, so turn into SCIP_EXPR_REALPOWER */
4063  SCIP_Real exponent;
4064 
4065  exponent = monomial->exponents[0];
4066 
4067  polynomialdataFree(blkmem, &polynomialdata);
4068 
4069  *op = SCIP_EXPR_REALPOWER;
4070  data->dbl = exponent;
4071 
4072  return SCIP_OKAY;
4073  }
4074  }
4075 
4076  if( maxdegree == 2 && monomial->nfactors == 2 )
4077  {
4078  /* polynomial is product of two children, so turn into SCIP_EXPR_MUL */
4079  assert(monomial->exponents[0] == 1.0);
4080  assert(monomial->exponents[1] == 1.0);
4081 
4082  polynomialdataFree(blkmem, &polynomialdata);
4083  data->data = NULL;
4084 
4085  *op = SCIP_EXPR_MUL;
4086 
4087  return SCIP_OKAY;
4088  }
4089 
4090  if( maxdegree == monomial->nfactors )
4091  {
4092  /* polynomial is a product of n children, so turn into SCIP_EXPR_PRODUCT */
4093 
4094  polynomialdataFree(blkmem, &polynomialdata);
4095  data->data = NULL;
4096 
4097  *op = SCIP_EXPR_PRODUCT;
4098 
4099  return SCIP_OKAY;
4100  }
4101 
4102  if( monomial->nfactors == 2 && monomial->exponents[0] == 1.0 && monomial->exponents[1] == -1.0 )
4103  {
4104  /* polynomial is x/y, so turn into SCIP_EXPR_DIV */
4105  assert(monomial->childidxs[0] == 0);
4106  assert(monomial->childidxs[1] == 1);
4107 
4108  polynomialdataFree(blkmem, &polynomialdata);
4109  data->data = NULL;
4110 
4111  *op = SCIP_EXPR_DIV;
4112 
4113  return SCIP_OKAY;
4114  }
4115 
4116  if( monomial->nfactors == 2 && monomial->exponents[0] == -1.0 && monomial->exponents[1] == 1.0 )
4117  {
4118  /* polynomial is y/x, so turn into SCIP_EXPR_DIV */
4119  void* tmp;
4120 
4121  assert(monomial->childidxs[0] == 0);
4122  assert(monomial->childidxs[1] == 1);
4123 
4124  polynomialdataFree(blkmem, &polynomialdata);
4125  data->data = NULL;
4126 
4127  /* swap children */
4128  tmp = children[1]; /*lint !e613*/
4129  children[1] = children[0]; /*lint !e613*/
4130  children[0] = tmp; /*lint !e613*/
4131 
4132  *op = SCIP_EXPR_DIV;
4133 
4134  return SCIP_OKAY;
4135  }
4136  }
4137 
4138  return SCIP_OKAY;
4139 }
4140 
4141 /** adds copies of expressions to the array of children of a sum, product, linear, quadratic, or polynomial expression
4142  *
4143  * For a sum or product expression, this corresponds to add additional summands and factors, resp.
4144  * For a linear expression, this corresponds to add each expression with coefficient 1.0.
4145  * For a quadratic or polynomial expression, only the children array may be enlarged, the expression itself remains the same.
4146  */
4147 static
4149  BMS_BLKMEM* blkmem, /**< block memory */
4150  SCIP_EXPR* expr, /**< quadratic or polynomial expression */
4151  int nexprs, /**< number of expressions to add */
4152  SCIP_EXPR** exprs, /**< expressions to add */
4153  SCIP_Bool comparechildren, /**< whether to compare expressions with already existing children (no effect for sum and product) */
4154  SCIP_Real eps, /**< which epsilon to use when comparing expressions */
4155  int* childmap /**< array where to store mapping of indices from exprs to children array in expr, or NULL if not of interest */
4156  )
4157 {
4158  int i;
4159 
4160  assert(blkmem != NULL);
4161  assert(expr != NULL);
4162  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);
4163  assert(exprs != NULL || nexprs == 0);
4164 
4165  if( nexprs == 0 )
4166  return SCIP_OKAY;
4167 
4168  switch( expr->op )
4169  {
4170  case SCIP_EXPR_SUM:
4171  case SCIP_EXPR_PRODUCT:
4172  {
4173  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, expr->nchildren, expr->nchildren + nexprs) );
4174  for( i = 0; i < nexprs; ++i )
4175  {
4176  SCIP_CALL( SCIPexprCopyDeep(blkmem, &expr->children[expr->nchildren + i], exprs[i]) ); /*lint !e613*/
4177  if( childmap != NULL )
4178  childmap[i] = expr->nchildren + i;
4179  }
4180  expr->nchildren += nexprs;
4181 
4182  break;
4183  }
4184 
4185  case SCIP_EXPR_LINEAR:
4186  case SCIP_EXPR_QUADRATIC:
4187  case SCIP_EXPR_POLYNOMIAL:
4188  {
4189  int j;
4190  int orignchildren;
4191  SCIP_Bool existsalready;
4192 
4193  orignchildren = expr->nchildren;
4194  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, expr->nchildren, expr->nchildren + nexprs) );
4195 
4196  for( i = 0; i < nexprs; ++i )
4197  {
4198  existsalready = FALSE;
4199  if( comparechildren )
4200  for( j = 0; j < orignchildren; ++j )
4201  /* during simplification of polynomials, their may be NULL's in children array */
4202  if( expr->children[j] != NULL && SCIPexprAreEqual(expr->children[j], exprs[i], eps) ) /*lint !e613*/
4203  {
4204  existsalready = TRUE;
4205  break;
4206  }
4207 
4208  if( !existsalready )
4209  {
4210  /* add copy of exprs[j] to children array */
4211  SCIP_CALL( SCIPexprCopyDeep(blkmem, &expr->children[expr->nchildren], exprs[i]) ); /*lint !e613*/
4212  if( childmap != NULL )
4213  childmap[i] = expr->nchildren;
4214  ++expr->nchildren;
4215  }
4216  else
4217  {
4218  if( childmap != NULL )
4219  childmap[i] = j; /*lint !e644*/
4220  if( expr->op == SCIP_EXPR_LINEAR )
4221  {
4222  /* if linear expression, increase coefficient by 1.0 */
4223  ((SCIP_Real*)expr->data.data)[j] += 1.0;
4224  }
4225  }
4226  }
4227 
4228  /* shrink children array to actually used size */
4229  assert(comparechildren || expr->nchildren == orignchildren + nexprs);
4230  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, orignchildren + nexprs, expr->nchildren) );
4231 
4232  if( expr->op == SCIP_EXPR_LINEAR && expr->nchildren > orignchildren )
4233  {
4234  /* if linear expression, then add 1.0 coefficients for new expressions */
4235  SCIP_Real* data;
4236 
4237  data = (SCIP_Real*)expr->data.data;
4238  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &data, orignchildren + 1, expr->nchildren + 1) );
4239  data[expr->nchildren] = data[orignchildren]; /* move constant from old end to new end */
4240  for( i = orignchildren; i < expr->nchildren; ++i )
4241  data[i] = 1.0;
4242  expr->data.data = (void*)data;
4243  }
4244  else if( expr->op == SCIP_EXPR_QUADRATIC && expr->nchildren > orignchildren )
4245  {
4246  /* if quadratic expression, then add 0.0 linear coefficients for new expressions */
4248 
4249  data = (SCIP_EXPRDATA_QUADRATIC*)expr->data.data;
4250  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &data->lincoefs, orignchildren, expr->nchildren) );
4251  BMSclearMemoryArray(&data->lincoefs[orignchildren], expr->nchildren - orignchildren); /*lint !e866*/
4252  }
4253 
4254  break;
4255  }
4256 
4257  default:
4258  SCIPerrorMessage("exprsimplifyAddChildren cannot be called for operand %d\n", expr->op);
4259  return SCIP_INVALIDDATA;
4260  } /*lint !e788*/
4261 
4262  return SCIP_OKAY;
4263 }
4264 
4265 /** converts expressions into polynomials, where possible and obvious */
4266 static
4268  BMS_BLKMEM* blkmem, /**< block memory data structure */
4269  SCIP_EXPR* expr /**< expression to convert */
4270  )
4271 {
4272  int i;
4273 
4274  assert(expr != NULL);
4275 
4276  for( i = 0; i < expr->nchildren; ++i )
4277  {
4279  }
4280 
4281  SCIP_CALL( exprConvertToPolynomial(blkmem, &expr->op, &expr->data, expr->nchildren) );
4282 
4283  return SCIP_OKAY;
4284 }
4285 
4286 /** removes duplicate children in a polynomial expression
4287  *
4288  * Leaves NULL's in children array.
4289  */
4290 static
4292  BMS_BLKMEM* blkmem, /**< block memory data structure */
4293  SCIP_EXPR* expr, /**< expression */
4294  SCIP_Real eps /**< threshold for zero */
4295  )
4296 {
4297  SCIP_Bool foundduplicates;
4298  int* childmap;
4299  int i;
4300  int j;
4301 
4302  assert(blkmem != NULL);
4303  assert(expr != NULL);
4304  assert(SCIPexprGetOperator(expr) == SCIP_EXPR_POLYNOMIAL);
4305 
4306  if( expr->nchildren == 0 )
4307  return SCIP_OKAY;
4308 
4309  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &childmap, expr->nchildren) );
4310 
4311  foundduplicates = FALSE;
4312  for( i = 0; i < expr->nchildren; ++i )
4313  {
4314  if( expr->children[i] == NULL )
4315  continue;
4316  childmap[i] = i; /*lint !e644*/
4317 
4318  for( j = i+1; j < expr->nchildren; ++j )
4319  {
4320  if( expr->children[j] == NULL )
4321  continue;
4322 
4323  if( SCIPexprAreEqual(expr->children[i], expr->children[j], eps) )
4324  {
4325  /* forget about expr j and remember that is to be replaced by i */
4326  SCIPexprFreeDeep(blkmem, &expr->children[j]);
4327  childmap[j] = i;
4328  foundduplicates = TRUE;
4329  }
4330  }
4331  }
4332 
4333  /* apply childmap to monomials */
4334  if( foundduplicates )
4336 
4337  /* free childmap */
4338  BMSfreeBlockMemoryArray(blkmem, &childmap, expr->nchildren);
4339 
4340  return SCIP_OKAY;
4341 }
4342 
4343 /** eliminates NULL's in children array and shrinks it to actual size */
4344 static
4346  BMS_BLKMEM* blkmem, /**< block memory data structure */
4347  SCIP_EXPR* expr /**< expression */
4348  )
4349 {
4350  int* childmap;
4351  int lastnonnull;
4352  int i;
4353 
4354  assert(blkmem != NULL);
4355  assert(expr != NULL);
4356  assert(SCIPexprGetOperator(expr) == SCIP_EXPR_POLYNOMIAL);
4357 
4358  if( expr->nchildren == 0 )
4359  return SCIP_OKAY;
4360 
4361  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &childmap, expr->nchildren) );
4362 
4363  /* close gaps in children array */
4364  lastnonnull = expr->nchildren-1;
4365  while( lastnonnull >= 0 && expr->children[lastnonnull] == NULL )
4366  --lastnonnull;
4367  for( i = 0; i <= lastnonnull; ++i )
4368  {
4369  if( expr->children[i] != NULL )
4370  {
4371  childmap[i] = i; /* child at index i is not moved */ /*lint !e644*/
4372  continue;
4373  }
4374  assert(expr->children[lastnonnull] != NULL);
4375 
4376  /* move child at lastnonnull to position i */
4377  expr->children[i] = expr->children[lastnonnull];
4378  expr->children[lastnonnull] = NULL;
4379  childmap[lastnonnull] = i;
4380 
4381  /* update lastnonnull */
4382  --lastnonnull;
4383  while( lastnonnull >= 0 && expr->children[lastnonnull] == NULL )
4384  --lastnonnull;
4385  }
4386  assert(i > lastnonnull);
4387 
4388  /* apply childmap to monomials */
4389  if( lastnonnull < expr->nchildren-1 )
4391 
4392  BMSfreeBlockMemoryArray(blkmem, &childmap, expr->nchildren);
4393 
4394  /* shrink children array */
4395  if( lastnonnull >= 0 )
4396  {
4397  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, expr->nchildren, lastnonnull+1) );
4398  expr->nchildren = lastnonnull+1;
4399  }
4400  else
4401  {
4402  BMSfreeBlockMemoryArray(blkmem, &expr->children, expr->nchildren);
4403  expr->nchildren = 0;
4404  }
4405 
4406  return SCIP_OKAY;
4407 }
4408 
4409 /** checks which children are still in use and frees those which are not */
4410 static
4412  BMS_BLKMEM* blkmem, /**< block memory data structure */
4413  SCIP_EXPR* expr /**< polynomial expression */
4414  )
4415 {
4416  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
4417  SCIP_EXPRDATA_MONOMIAL* monomial;
4418  SCIP_Bool* childinuse;
4419  int i;
4420  int j;
4421 
4422  assert(blkmem != NULL);
4423  assert(expr != NULL);
4424 
4425  if( expr->nchildren == 0 )
4426  return SCIP_OKAY;
4427 
4428  polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)expr->data.data;
4429  assert(polynomialdata != NULL);
4430 
4431  /* check which children are still in use */
4432  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &childinuse, expr->nchildren) );
4433  BMSclearMemoryArray(childinuse, expr->nchildren); /*lint !e644*/
4434  for( i = 0; i < polynomialdata->nmonomials; ++i )
4435  {
4436  monomial = polynomialdata->monomials[i];
4437  assert(monomial != NULL);
4438 
4439  for( j = 0; j < monomial->nfactors; ++j )
4440  {
4441  assert(monomial->childidxs[j] >= 0);
4442  assert(monomial->childidxs[j] < expr->nchildren);
4443  childinuse[monomial->childidxs[j]] = TRUE;
4444  }
4445  }
4446 
4447  /* free children that are not used in any monomial */
4448  for( i = 0; i < expr->nchildren; ++i )
4449  if( expr->children[i] != NULL && !childinuse[i] )
4450  SCIPexprFreeDeep(blkmem, &expr->children[i]);
4451 
4452  BMSfreeBlockMemoryArray(blkmem, &childinuse, expr->nchildren);
4453 
4454  return SCIP_OKAY;
4455 }
4456 
4457 /** flattens polynomials in polynomials, check for constants in non-polynomials expressions
4458  *
4459  * exprsimplifyConvertToPolynomials should have been called before to eliminate simple polynomial operands.
4460  */
4461 static
4463  BMS_BLKMEM* blkmem, /**< block memory data structure */
4464  SCIP_MESSAGEHDLR* messagehdlr, /**< message handler */
4465  SCIP_EXPR* expr, /**< expression */
4466  SCIP_Real eps, /**< threshold, under which values are treat as 0 */
4467  int maxexpansionexponent/**< maximal exponent for which we still expand non-monomial polynomials */
4468  )
4469 {
4470  int i;
4471 
4472  assert(expr != NULL);
4473 
4474  for( i = 0; i < expr->nchildren; ++i )
4475  {
4476  SCIP_CALL( exprsimplifyFlattenPolynomials(blkmem, messagehdlr, expr->children[i], eps, maxexpansionexponent) );
4477  }
4478 
4479  switch( SCIPexprGetOperator(expr) )
4480  {
4481  case SCIP_EXPR_VARIDX:
4482  case SCIP_EXPR_CONST:
4483  case SCIP_EXPR_PARAM:
4484  case SCIP_EXPR_PLUS:
4485  case SCIP_EXPR_MINUS:
4486  case SCIP_EXPR_MUL:
4487  case SCIP_EXPR_DIV:
4488  case SCIP_EXPR_SQUARE:
4489  case SCIP_EXPR_SQRT:
4490  case SCIP_EXPR_INTPOWER:
4491  case SCIP_EXPR_REALPOWER:
4492  case SCIP_EXPR_SIGNPOWER:
4493  break;
4494 
4495  case SCIP_EXPR_EXP:
4496  case SCIP_EXPR_LOG:
4497  case SCIP_EXPR_SIN:
4498  case SCIP_EXPR_COS:
4499  case SCIP_EXPR_TAN:
4500  /* case SCIP_EXPR_ERF: */
4501  /* case SCIP_EXPR_ERFI: */
4502  case SCIP_EXPR_ABS:
4503  case SCIP_EXPR_SIGN:
4504  {
4505  /* check if argument is a constant */
4506  if( (expr->children[0]->op == SCIP_EXPR_POLYNOMIAL && SCIPexprGetNChildren(expr->children[0]) == 0) ||
4507  expr->children[0]->op == SCIP_EXPR_CONST )
4508  {
4509  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
4510  SCIP_Real exprval;
4511 
4512  /* since child0 has no children and it's polynomial was flattened, it should have no monomials */
4513  assert(expr->children[0]->op != SCIP_EXPR_POLYNOMIAL || SCIPexprGetNMonomials(expr->children[0]) == 0);
4514 
4515  /* evaluate expression in constant polynomial */
4516  SCIP_CALL( SCIPexprEval(expr, NULL, NULL, &exprval) );
4517 
4518  /* create polynomial */
4519  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 0, NULL, exprval, FALSE) );
4520 
4521  expr->op = SCIP_EXPR_POLYNOMIAL;
4522  expr->data.data = (void*)polynomialdata;
4523 
4524  /* forget child */
4525  SCIPexprFreeDeep(blkmem, &expr->children[0]);
4526  BMSfreeBlockMemoryArray(blkmem, &expr->children, 1);
4527  expr->nchildren = 0;
4528  }
4529 
4530  break;
4531  }
4532 
4533  case SCIP_EXPR_MIN:
4534  case SCIP_EXPR_MAX:
4535  {
4536  /* check if both arguments are constants */
4537  if( ((expr->children[0]->op == SCIP_EXPR_POLYNOMIAL && SCIPexprGetNChildren(expr->children[0]) == 0) || expr->children[0]->op == SCIP_EXPR_CONST) &&
4538  ((expr->children[1]->op == SCIP_EXPR_POLYNOMIAL && SCIPexprGetNChildren(expr->children[1]) == 0) || expr->children[1]->op == SCIP_EXPR_CONST) )
4539  {
4540  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
4541  SCIP_Real exprval;
4542 
4543  /* since children have no children and it's polynomial was flattened, it should have no monomials */
4544  assert(expr->children[0]->op != SCIP_EXPR_POLYNOMIAL || SCIPexprGetNMonomials(expr->children[0]) == 0);
4545  assert(expr->children[1]->op != SCIP_EXPR_POLYNOMIAL || SCIPexprGetNMonomials(expr->children[1]) == 0);
4546 
4547  /* evaluate expression in constants */
4548  SCIP_CALL( SCIPexprEval(expr, NULL, NULL, &exprval) );
4549 
4550  /* create polynomial */
4551  SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 0, NULL, exprval, FALSE) );
4552 
4553  expr->op = SCIP_EXPR_POLYNOMIAL;
4554  expr->data.data = (void*)polynomialdata;
4555 
4556  /* forget children */
4557  SCIPexprFreeDeep(blkmem, &expr->children[0]);
4558  SCIPexprFreeDeep(blkmem, &expr->children[1]);
4559  BMSfreeBlockMemoryArray(blkmem, &expr->children, 2);
4560  expr->nchildren = 0;
4561  }
4562 
4563  break;
4564  }
4565 
4566  case SCIP_EXPR_SUM:
4567  case SCIP_EXPR_PRODUCT:
4568  case SCIP_EXPR_LINEAR:
4569  case SCIP_EXPR_QUADRATIC:
4570  case SCIP_EXPR_USER:
4571  break;
4572 
4573  case SCIP_EXPR_POLYNOMIAL:
4574  {
4575  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
4576  SCIP_EXPRDATA_MONOMIAL* monomial;
4577  SCIP_Bool removechild;
4578  int* childmap;
4579  int childmapsize;
4580  int j;
4581 
4582  /* simplify current polynomial */
4584  SCIPexprMergeMonomials(blkmem, expr, eps, TRUE);
4585 
4586  polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)expr->data.data;
4587  assert(polynomialdata != NULL);
4588 
4589  SCIPdebugMessage("expand factors in expression ");
4590  SCIPdebug( SCIPexprPrint(expr, messagehdlr, NULL, NULL, NULL, NULL) );
4591  SCIPdebugPrintf("\n");
4592 
4593  childmap = NULL;
4594  childmapsize = 0;
4595 
4596  /* resolve children that are constants
4597  * we do this first, because it reduces the degree and number of factors in the monomials,
4598  * 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
4599  */
4600  for( i = 0; i < expr->nchildren; ++i )
4601  {
4602  if( expr->children[i] == NULL )
4603  continue;
4604 
4605  if( SCIPexprGetOperator(expr->children[i]) != SCIP_EXPR_CONST )
4606  continue;
4607 
4608  removechild = TRUE; /* we intend to delete children[i] */
4609 
4610  if( childmapsize < expr->children[i]->nchildren )
4611  {
4612  int newsize;
4613 
4614  newsize = calcGrowSize(expr->children[i]->nchildren);
4615  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &childmap, childmapsize, newsize) );
4616  childmapsize = newsize;
4617  }
4618 
4619  /* put constant of child i into every monomial where child i is used */
4620  for( j = 0; j < polynomialdata->nmonomials; ++j )
4621  {
4622  int factorpos;
4623  SCIP_Bool success;
4624 
4625  monomial = polynomialdata->monomials[j];
4626  /* if monomial is not sorted, then polynomial should not be sorted either, or have only one monomial */
4627  assert(monomial->sorted || !polynomialdata->sorted || polynomialdata->nmonomials <= 1);
4628 
4629  if( SCIPexprFindMonomialFactor(monomial, i, &factorpos) )
4630  {
4631  assert(factorpos >= 0);
4632  assert(factorpos < monomial->nfactors);
4633  /* assert that factors have been merged */
4634  assert(factorpos == 0 || monomial->childidxs[factorpos-1] != i);
4635  assert(factorpos == monomial->nfactors-1 || monomial->childidxs[factorpos+1] != i);
4636 
4637  /* SCIPdebugMessage("attempt expanding child %d at monomial %d factor %d\n", i, j, factorpos);
4638  SCIPdebug( SCIPexprPrint(expr, NULL, NULL, NULL) ); SCIPdebugPrintf("\n");
4639  SCIPdebug( SCIPexprPrint(expr->children[i], NULL, NULL, NULL) ); SCIPdebugPrintf("\n"); */
4640 
4641  if( !EPSISINT(monomial->exponents[factorpos], 0.0) && SCIPexprGetOpReal(expr->children[i]) < 0.0 ) /*lint !e835*/
4642  {
4643  /* if constant is negative and our exponent is not integer, then cannot do expansion */
4644  SCIPmessagePrintWarning(messagehdlr, "got negative constant %g to the power of a noninteger exponent %g\n",
4645  SCIPexprGetOpReal(expr->children[i]), monomial->exponents[factorpos]);
4646  success = FALSE;
4647  }
4648  else
4649  {
4650  monomial->coef *= pow(SCIPexprGetOpReal(expr->children[i]), monomial->exponents[factorpos]);
4651 
4652  /* move last factor to position factorpos */
4653  if( factorpos < monomial->nfactors-1 )
4654  {
4655  monomial->exponents[factorpos] = monomial->exponents[monomial->nfactors-1];
4656  monomial->childidxs[factorpos] = monomial->childidxs[monomial->nfactors-1];
4657  }
4658  --monomial->nfactors;
4659  monomial->sorted = FALSE;
4660  polynomialdata->sorted = FALSE;
4661 
4662  success = TRUE;
4663  }
4664 
4665  if( !success )
4666  removechild = FALSE;
4667  }
4668  }
4669 
4670  /* forget about child i, if it is not used anymore */
4671  if( removechild )
4672  SCIPexprFreeDeep(blkmem, &expr->children[i]);
4673 
4674  /* simplify current polynomial again */
4675  SCIPexprMergeMonomials(blkmem, expr, eps, TRUE);
4676  }
4677 
4678  /* try to resolve children that are polynomials itself */
4679  for( i = 0; i < expr->nchildren; ++i )
4680  {
4681  if( expr->children[i] == NULL )
4682  continue;
4683 
4685  continue;
4686 
4687  removechild = TRUE; /* we intend to delete children[i] */
4688 
4689  if( childmapsize < expr->children[i]->nchildren )
4690  {
4691  int newsize;
4692 
4693  newsize = calcGrowSize(expr->children[i]->nchildren);
4694  SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &childmap, childmapsize, newsize) );
4695  childmapsize = newsize;
4696  }
4697 
4698  /* add children of child i */
4699  SCIP_CALL( exprsimplifyAddChildren(blkmem, expr, expr->children[i]->nchildren, expr->children[i]->children, TRUE, eps, childmap) );
4700 
4701  /* put polynomial of child i into every monomial where child i is used */
4702  j = 0;
4703  while( j < polynomialdata->nmonomials )
4704  {
4705  int factorpos;
4706  SCIP_Bool success;
4707 
4708  monomial = polynomialdata->monomials[j];
4709  /* if monomial is not sorted, then polynomial should not be sorted either, or have only one monomial */
4710  assert(monomial->sorted || !polynomialdata->sorted || polynomialdata->nmonomials <= 1);
4711 
4712  if( SCIPexprFindMonomialFactor(monomial, i, &factorpos) )
4713  {
4714  assert(factorpos >= 0);
4715  assert(factorpos < monomial->nfactors);
4716  /* assert that factors have been merged */
4717  assert(factorpos == 0 || monomial->childidxs[factorpos-1] != i);
4718  assert(factorpos == monomial->nfactors-1 || monomial->childidxs[factorpos+1] != i);
4719 
4720  /* SCIPdebugMessage("attempt expanding child %d at monomial %d factor %d\n", i, j, factorpos);
4721  SCIPdebug( SCIPexprPrint(expr, NULL, NULL, NULL) ); SCIPdebugPrintf("\n");
4722  SCIPdebug( SCIPexprPrint(expr->children[i], NULL, NULL, NULL) ); SCIPdebugPrintf("\n"); */
4723 
4724  SCIP_CALL( polynomialdataExpandMonomialFactor(blkmem, messagehdlr, polynomialdata, j, factorpos,
4725  (SCIP_EXPRDATA_POLYNOMIAL*)expr->children[i]->data.data, childmap, maxexpansionexponent, &success) );
4726 
4727  if( !success )
4728  {
4729  removechild = FALSE;
4730  ++j;
4731  }
4732  }
4733  else
4734  ++j;
4735 
4736  /* expansion may remove monomials[j], move a monomial from the end to position j, or add new monomials to the end of polynomialdata
4737  * we thus repeat with index j, if a factor was successfully expanded
4738  */
4739  }
4740 
4741  /* forget about child i, if it is not used anymore */
4742  if( removechild )
4743  SCIPexprFreeDeep(blkmem, &expr->children[i]);
4744 
4745  /* simplify current polynomial again */
4746  SCIPexprMergeMonomials(blkmem, expr, eps, TRUE);
4747  }
4748 
4749  BMSfreeBlockMemoryArrayNull(blkmem, &childmap, childmapsize);
4750 
4751  /* free children that are not in use anymore */
4753 
4754  /* remove NULLs from children array */
4756 
4757  /* if no children left, then it's a constant polynomial -> change into EXPR_CONST */
4758  if( expr->nchildren == 0 )
4759  {
4760  SCIP_Real val;
4761 
4762  /* if no children, then it should also have no monomials */
4763  assert(polynomialdata->nmonomials == 0);
4764 
4765  val = polynomialdata->constant;
4766  polynomialdataFree(blkmem, &polynomialdata);
4767 
4768  expr->op = SCIP_EXPR_CONST;
4769  expr->data.dbl = val;
4770  }
4771 
4772  SCIPdebugMessage("-> ");
4773  SCIPdebug( SCIPexprPrint(expr, messagehdlr, NULL, NULL, NULL, NULL) );
4774  SCIPdebugPrintf("\n");
4775 
4776  break;
4777  }
4778 
4779  case SCIP_EXPR_LAST:
4780  break;
4781  } /*lint !e788*/
4782 
4783  return SCIP_OKAY;
4784 }
4785 
4786 /** separates linear monomials from an expression, if it is a polynomial expression
4787  *
4788  * Separates only those linear terms whose variable is not used otherwise in the expression.
4789  */
4790 static
4792  BMS_BLKMEM* blkmem, /**< block memory data structure */
4793  SCIP_EXPR* expr, /**< expression */
4794  SCIP_Real eps, /**< threshold, under which positive values are treat as 0 */
4795  int nvars, /**< number of variables in expression */
4796  int* nlinvars, /**< buffer to store number of linear variables in linear part */
4797  int* linidxs, /**< array to store indices of variables in expression tree which belong to linear part */
4798  SCIP_Real* lincoefs /**< array to store coefficients of linear part */
4799  )
4800 {
4801  SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
4802  SCIP_EXPRDATA_MONOMIAL* monomial;
4803  int* varsusage;
4804  int* childusage;
4805  int childidx;
4806  int i;
4807  int j;
4808 
4809  assert(blkmem != NULL);
4810  assert(expr != NULL);
4811  assert(nlinvars != NULL);
4812  assert(linidxs != NULL);
4813  assert(lincoefs != NULL);
4814 
4815  *nlinvars = 0;
4816 
4818  return SCIP_OKAY;
4819 
4820  if( SCIPexprGetNChildren(expr) == 0 )
4821  return SCIP_OKAY;
4822 
4823  polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)expr->data.data;
4824  assert(polynomialdata != NULL);
4825 
4826  /* get variable usage */
4827  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &varsusage, nvars) );
4828  BMSclearMemoryArray(varsusage, nvars); /*lint !e644*/
4829  SCIPexprGetVarsUsage(expr, varsusage);
4830 
4831  /* get child usage: how often each child is used in the polynomial */
4832  SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &childusage, expr->nchildren) );
4833  BMSclearMemoryArray(childusage, expr->nchildren); /*lint !e644*/
4834  for( i = 0; i < polynomialdata->nmonomials; ++i )
4835  {
4836  monomial = polynomialdata->monomials[i];
4837  assert(monomial != NULL);
4838  for( j = 0; j < monomial->nfactors; ++j )
4839  {
4840  assert(monomial->childidxs[j] >= 0);
4841  assert(monomial->childidxs[j] < expr->nchildren);
4842  ++childusage[monomial->childidxs[j]];
4843  }
4844  }
4845 
4846  /* move linear monomials out of polynomial */
4847  for( i = 0; i < polynomialdata->nmonomials; ++i )
4848  {
4849  monomial = polynomialdata->monomials[i];
4850  assert(monomial != NULL);
4851  if( monomial->nfactors != 1 )
4852  continue;
4853  if( monomial->exponents[0] != 1.0 )
4854  continue;
4855  childidx = monomial->childidxs[0];
4856  if( SCIPexprGetOperator(expr->children[childidx]) != SCIP_EXPR_VARIDX )
4857  continue;
4858 
4859  /* we are at a linear monomial in a variable */
4860  assert(SCIPexprGetOpIndex(expr->children[childidx]) < nvars);
4861  if( childusage[childidx] == 1 && varsusage[SCIPexprGetOpIndex(expr->children[childidx])] == 1 )
4862  {
4863  /* if the child expression is not used in another monomial (which would due to merging be not linear)
4864  * and if the variable is not used somewhere else in the tree,
4865  * then move this monomial into linear part and free child
4866  */
4867  linidxs[*nlinvars] = SCIPexprGetOpIndex(expr->children[childidx]);
4868  lincoefs[*nlinvars] = monomial->coef;
4869  ++*nlinvars;
4870 
4871  SCIPexprFreeDeep(blkmem, &expr->children[childidx]);
4872  monomial->coef = 0.0;
4873  monomial->nfactors = 0;
4874  }
4875  }
4876 
4877  BMSfreeBlockMemoryArray(blkmem, &varsusage, nvars);
4878  BMSfreeBlockMemoryArray(blkmem, &childusage, expr->nchildren);
4879 
4880  if( *nlinvars > 0 )
4881  {
4882  /* if we did something, cleanup polynomial (e.g., remove monomials with coefficient 0.0) */
4883  polynomialdataMergeMonomials(blkmem, polynomialdata, eps, FALSE);
4885  }
4886 
4887  return SCIP_OKAY;
4888 }
4889 
4890 /** converts polynomial expressions back into simpler expressions, where possible */
4891 static
4893  BMS_BLKMEM* blkmem, /**< block memory data structure */
4894  SCIP_EXPR* expr /**< expression to convert back */
4895  )
4896 {
4897  int i;
4898 
4899  assert(blkmem != NULL);
4900  assert(expr != NULL);
4901 
4902  for( i = 0; i < expr->nchildren; ++i )
4903  {
4905  }
4906 
4907  if( expr->op != SCIP_EXPR_POLYNOMIAL )
4908  return SCIP_OKAY;
4909 
4910  SCIP_CALL( exprUnconvertPolynomial(blkmem, &expr->op, &expr->data, expr->nchildren, (void**)expr->children) );
4911 
4912  return SCIP_OKAY;
4913 }
4914 
4915 static
4916 SCIP_DECL_HASHGETKEY( exprparseVarTableGetKey )
4917 { /*lint --e{715}*/
4918  return (void*)((char*)elem + sizeof(int));
4919 }
4920 
4921 /** parses a variable name from a string and creates corresponding expression
4922  *
4923  * Creates a new variable index if variable not seen before, updates varnames and vartable structures.
4924  */
4925 static
4927  BMS_BLKMEM* blkmem, /**< block memory data structure */
4928  const char** str, /**< pointer to the string to be parsed */
4929  SCIP_EXPR** expr, /**< buffer to store pointer to created expression */
4930  int* nvars, /**< running number of encountered variables so far */
4931  int** varnames, /**< pointer to buffer to store new variable names */
4932  int* varnameslength, /**< pointer to length of the varnames buffer array */
4933  SCIP_HASHTABLE* vartable, /**< hash table for variable names and corresponding expression index */
4934  SCIP_Real coefficient, /**< coefficient to be used when creating the expression */
4935  const char* varnameendptr /**< if a <varname> should be parsed, set this to NULL. Then, str points to the '<'
4936  else, str should point to the first letter of the varname, and varnameendptr should
4937  point one char behind the last char of the variable name */
4938  )
4939 {
4940  int namelength;
4941  int varidx;
4942  char varname[SCIP_MAXSTRLEN];
4943  void* element;
4944 
4945  assert(blkmem != NULL);
4946  assert(str != NULL);
4947  assert(expr != NULL);
4948  assert(nvars != NULL);
4949  assert(varnames != NULL);
4950  assert(vartable != NULL);
4951 
4952  if( varnameendptr == NULL )
4953  {
4954  ++*str;
4955  varnameendptr = *str;
4956  while( varnameendptr[0] != '>' )
4957  ++varnameendptr;
4958  }
4959 
4960  namelength = varnameendptr - *str; /*lint !e712*/
4961  if( namelength >= SCIP_MAXSTRLEN )
4962  {
4963  SCIPerrorMessage("Variable name %.*s is too long for buffer in exprparseReadVariable.\n", namelength, str);
4964  return SCIP_READERROR;
4965  }
4966 
4967  memcpy(varname, *str, namelength * sizeof(char));
4968  varname[namelength] = '\0';
4969 
4970  element = SCIPhashtableRetrieve(vartable, varname);
4971  if( element != NULL )
4972  {
4973  /* variable is old friend */
4974  assert(strcmp((char*)element + sizeof(int), varname) == 0);
4975 
4976  varidx = *(int*)element;
4977  }
4978  else
4979  {
4980  /* variable is new */
4981  varidx = *nvars;
4982 
4983  (*varnameslength) -= (int)(1 + (strlen(varname) + 1) / sizeof(int) + 1);
4984  if( *varnameslength < 0 )
4985  {
4986  SCIPerrorMessage("Buffer in exprparseReadVariable is too short for varaible name %.*s.\n", namelength, str);
4987  return SCIP_READERROR;
4988  }
4989 
4990  /* store index of variable and variable name in varnames buffer */
4991  **varnames = varidx;
4992  (void) SCIPstrncpy((char*)(*varnames + 1), varname, (int)strlen(varname)+1);
4993 
4994  /* insert variable into hashtable */
4995  SCIP_CALL( SCIPhashtableInsert(vartable, (void*)*varnames) );
4996 
4997  ++*nvars;
4998  *varnames += 1 + (strlen(varname) + 1) / sizeof(int) + 1;
4999  }
5000 
5001  /* create VARIDX expression, put into LINEAR expression if we have coefficient != 1 */
5002  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_VARIDX, varidx) ); /*lint !e613*/
5003  if( coefficient != 1.0 )
5004  {
5005  SCIP_CALL( SCIPexprCreateLinear(blkmem, expr, 1, expr, &coefficient, 0.0) );
5006  }
5007 
5008  /* Move pointer to char behind end of variable */
5009  *str = varnameendptr + 1;
5010 
5011  /* consprint sometimes prints a variable type identifier which we don't need */
5012  if( (*str)[0] == '[' && (*str)[2] == ']' &&
5013  ((*str)[1] == SCIP_VARTYPE_BINARY_CHAR ||
5014  (*str)[1] == SCIP_VARTYPE_INTEGER_CHAR ||
5015  (*str)[1] == SCIP_VARTYPE_IMPLINT_CHAR ||
5016  (*str)[1] == SCIP_VARTYPE_CONTINUOUS_CHAR ) )
5017  *str += 3;
5018 
5019  return SCIP_OKAY;
5020 }
5021 
5022 /** if str[0] points to an opening parenthesis, this function sets endptr to point to the matching closing bracket in str
5023  *
5024  * Searches for at most length characters.
5025  */
5026 static
5028  const char* str, /**< pointer to the string to be parsed */
5029  const char** endptr, /**< pointer to point to the closing parenthesis */
5030  int length /**< length of the string to be parsed */
5031  )
5032 {
5033  int nopenbrackets;
5034 
5035  assert(str[0] == '(');
5036 
5037  *endptr = str;
5038 
5039  /* find the end of this expression */
5040  nopenbrackets = 0;
5041  while( (*endptr - str ) < length && !(nopenbrackets == 1 && *endptr[0] == ')') )
5042  {
5043  if( *endptr[0] == '(')
5044  ++nopenbrackets;
5045  if( *endptr[0] == ')')
5046  --nopenbrackets;
5047  ++*endptr;
5048  }
5049 
5050  if( *endptr[0] != ')' )
5051  {
5052  SCIPerrorMessage("unable to find closing parenthesis in unbalanced expression %.*s\n", length, str);
5053  return SCIP_READERROR;
5054  }
5055 
5056  return SCIP_OKAY;
5057 }
5058 
5059 /** this function sets endptr to point to the next separating comma in str
5060  *
5061  * That is, for a given string like "x+f(x,y),z", endptr will point to the comma before "z"
5062  *
5063  * Searches for at most length characters.
5064  */
5065 static
5067  const char* str, /**< pointer to the string to be parsed */
5068  const char** endptr, /**< pointer to point to the comma */
5069  int length /**< length of the string to be parsed */
5070  )
5071 {
5072  int nopenbrackets;
5073 
5074  *endptr = str;
5075 
5076  /* find a comma without open brackets */
5077  nopenbrackets = 0;
5078  while( (*endptr - str ) < length && !(nopenbrackets == 0 && *endptr[0] == ',') )
5079  {
5080  if( *endptr[0] == '(')
5081  ++nopenbrackets;
5082  if( *endptr[0] == ')')
5083  --nopenbrackets;
5084  ++*endptr;
5085  }
5086 
5087  if( *endptr[0] != ',' )
5088  {
5089  SCIPerrorMessage("unable to find separating comma in unbalanced expression %.*s\n", length, str);
5090  return SCIP_READERROR;
5091  }
5092 
5093  return SCIP_OKAY;
5094 }
5095 
5096 /** parses an expression from a string */
5097 static
5099  BMS_BLKMEM* blkmem, /**< block memory data structure */
5100  SCIP_MESSAGEHDLR* messagehdlr, /**< message handler */
5101  SCIP_EXPR** expr, /**< buffer to store pointer to created expression */
5102  const char* str, /**< pointer to the string to be parsed */
5103  int length, /**< length of the string to be parsed */
5104  const char* lastchar, /**< pointer to the last char of str that should be parsed */
5105  int* nvars, /**< running number of encountered variables so far */
5106  int** varnames, /**< pointer to buffer to store new variable names */
5107  int* varnameslength, /**< pointer to length of the varnames buffer array */
5108  SCIP_HASHTABLE* vartable, /**< hash table for variable names and corresponding expression index */
5109  int recursiondepth /**< current recursion depth */
5110  )
5111 { /*lint --e{712,747}*/
5112  SCIP_EXPR* arg1;
5113  SCIP_EXPR* arg2;
5114  const char* subexpptr;
5115  const char* subexpendptr;
5116  const char* strstart;
5117  const char* endptr;
5118  char* nonconstendptr;
5119  SCIP_Real number;
5120  int subexplength;
5121  int nopenbrackets;
5122 
5123  assert(blkmem != NULL);
5124  assert(expr != NULL);
5125  assert(str != NULL);
5126  assert(lastchar >= str);
5127  assert(nvars != NULL);
5128  assert(varnames != NULL);
5129  assert(vartable != NULL);
5130 
5131  assert(recursiondepth < 100);
5132 
5133  strstart = str; /* might be needed for error message... */
5134 
5135  SCIPdebugMessage("exprParse (%i): parsing %.*s\n", recursiondepth, (int) (lastchar-str + 1), str);
5136 
5137  /* ignore whitespace */
5138  while( isspace((unsigned char)*str) )
5139  ++str;
5140 
5141  /* look for a sum or difference not contained in brackets */
5142  subexpptr = str;
5143  nopenbrackets = 0;
5144 
5145  /* find the end of this expression
5146  * a '+' right at the beginning indicates a coefficient, not treated here, or a summation
5147  */
5148  while( subexpptr != lastchar && !(nopenbrackets == 0 && (subexpptr[0] == '+' || subexpptr[0] == '-') && subexpptr != str) )
5149  {
5150  if( subexpptr[0] == '(')
5151  ++nopenbrackets;
5152  if( subexpptr[0] == ')')
5153  --nopenbrackets;
5154  ++subexpptr;
5155  }
5156 
5157  if( subexpptr != lastchar )
5158  {
5159  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str, (int) ((subexpptr - 1) - str + 1), subexpptr - 1, nvars,
5160  varnames, varnameslength, vartable, recursiondepth + 1) );
5161 
5162  if( subexpptr[0] == '+' )
5163  ++subexpptr;
5164  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, subexpptr , (int) (lastchar - (subexpptr ) + 1), lastchar, nvars,
5165  varnames, varnameslength, vartable, recursiondepth + 1) );
5166 
5167  /* make new expression from two arguments
5168  * we always use add, because we leave the operator between the found expressions in the second argument
5169  * this way, we do not have to worry about ''minus brackets'' in the case of more then two summands:
5170  * a - b - c = a + (-b -c)
5171  */
5172  SCIP_CALL( SCIPexprAdd(blkmem, expr, 1.0, arg1, 1.0, arg2, 0.0) );
5173 
5174  SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
5175  SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
5176  SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );
5177 
5178  return SCIP_OKAY;
5179  }
5180 
5181  /* check for a bracketed subexpression */
5182  if( str[0] == '(' )
5183  {
5184  nopenbrackets = 0;
5185 
5186  subexplength = -1; /* we do not want the closing bracket in the string */
5187  subexpptr = str + 1; /* leave out opening bracket */
5188 
5189  /* find the end of this expression */
5190  while( subexplength < length && !(nopenbrackets == 1 && str[0] == ')') )
5191  {
5192  if( str[0] == '(' )
5193  ++nopenbrackets;
5194  if( str[0] == ')' )
5195  --nopenbrackets;
5196  ++str;
5197  ++subexplength;
5198  }
5199  subexpendptr = str - 1; /* leave out closing bracket */
5200 
5201  SCIP_CALL( exprParse(blkmem, messagehdlr, expr, subexpptr, subexplength, subexpendptr, nvars, varnames,
5202  varnameslength, vartable, recursiondepth + 1) );
5203  ++str;
5204  }
5205  else if( isdigit((unsigned char)str[0]) || ((str[0] == '-' || str[0] == '+')
5206  && (isdigit((unsigned char)str[1]) || str[1] == ' ')) )
5207  {
5208  /* check if there is a lonely minus coming, indicating a -1.0 */
5209  if( str[0] == '-' && str[1] == ' ' )
5210  {
5211  number = -1.0;
5212  nonconstendptr = (char*) str + 1;
5213  }
5214  /* check if there is a number coming */
5215  else if( !SCIPstrToRealValue(str, &number, &nonconstendptr) )
5216  {
5217  SCIPerrorMessage("error parsing number from <%s>\n", str);
5218  return SCIP_READERROR;
5219  }
5220  str = nonconstendptr;
5221 
5222  /* ignore whitespace */
5223  while( isspace((unsigned char)*str) && str != lastchar )
5224  ++str;
5225 
5226  if( str[0] != '*' && str[0] != '/' && str[0] != '+' && str[0] != '-' && str[0] != '^' )
5227  {
5228  if( str < lastchar )
5229  {
5230  SCIP_CALL( exprParse(blkmem, messagehdlr, expr, str, (int)(lastchar - str) + 1, lastchar, nvars, varnames,
5231  varnameslength, vartable, recursiondepth + 1) );
5232  SCIP_CALL( SCIPexprMulConstant(blkmem, expr, *expr, number) );
5233  }
5234  else
5235  {
5236  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_CONST, number) );
5237  }
5238  str = lastchar + 1;
5239  }
5240  else
5241  {
5242  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_CONST, number) );
5243  }
5244  }
5245  else if( str[0] == '<' )
5246  {
5247  /* check if expressions begins with a variable */
5248  SCIP_CALL( exprparseReadVariable(blkmem, &str, expr, nvars, varnames, varnameslength, vartable, 1.0, NULL) );
5249  }
5250  /* four character operators */
5251  else if( strncmp(str, "sqrt", 4) == 0 )
5252  {
5253  str += 4;
5254  SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
5255  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, endptr - str - 1, endptr -1, nvars, varnames,
5256  varnameslength, vartable, recursiondepth + 1) );
5257  str = endptr + 1;
5258 
5259  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_SQRT, arg1) );
5260  }
5261  /* three character operators with 1 argument */
5262  else if(
5263  strncmp(str, "abs", 3) == 0 ||
5264  strncmp(str, "cos", 3) == 0 ||
5265  strncmp(str, "exp", 3) == 0 ||
5266  strncmp(str, "log", 3) == 0 ||
5267  strncmp(str, "sin", 3) == 0 ||
5268  strncmp(str, "sqr", 3) == 0 ||
5269  strncmp(str, "tan", 3) == 0 )
5270  {
5271  const char* opname = str;
5272 
5273  str += 3;
5274  SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
5275  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, endptr - str - 1, endptr -1, nvars, varnames,
5276  varnameslength, vartable, recursiondepth + 1) );
5277  str = endptr + 1;
5278 
5279  if( strncmp(opname, "abs", 3) == 0 )
5280  {
5281  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_ABS, arg1) );
5282  }
5283  else if( strncmp(opname, "cos", 3) == 0 )
5284  {
5285  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_COS, arg1) );
5286  }
5287  else if( strncmp(opname, "exp", 3) == 0 )
5288  {
5289  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_EXP, arg1) );
5290  }
5291  else if( strncmp(opname, "log", 3) == 0 )
5292  {
5293  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_LOG, arg1) );
5294  }
5295  else if( strncmp(opname, "sin", 3) == 0 )
5296  {
5297  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_SIN, arg1) );
5298  }
5299  else if( strncmp(opname, "sqr", 3) == 0 )
5300  {
5301  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_SQUARE, arg1) );
5302  }
5303  else
5304  {
5305  assert(strncmp(opname, "tan", 3) == 0);
5306  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_TAN, arg1) );
5307  }
5308  }
5309  /* three character operators with 2 arguments */
5310  else if(
5311  strncmp(str, "max", 3) == 0 ||
5312  strncmp(str, "min", 3) == 0 )
5313  {
5314  /* we have a string of the form "min(...,...)" or "max(...,...)"
5315  * first find the closing parenthesis, then the comma
5316  */
5317  const char* comma;
5318  SCIP_EXPROP op;
5319 
5320  op = (str[1] == 'a' ? SCIP_EXPR_MAX : SCIP_EXPR_MIN);
5321 
5322  str += 3;
5323  SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
5324 
5325  SCIP_CALL( exprparseFindSeparatingComma(str+1, &comma, endptr - str - 1) );
5326 
5327  /* parse first argument [str+1..comma-1] */
5328  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, comma - str - 1, comma - 1, nvars, varnames,
5329  varnameslength, vartable, recursiondepth + 1) );
5330 
5331  /* parse second argument [comma+1..endptr] */
5332  ++comma;
5333  while( comma < endptr && *comma == ' ' )
5334  ++comma;
5335 
5336  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, comma, endptr - comma, endptr - 1, nvars, varnames,
5337  varnameslength, vartable, recursiondepth + 1) );
5338 
5339  SCIP_CALL( SCIPexprCreate(blkmem, expr, op, arg1, arg2) );
5340 
5341  str = endptr + 1;
5342  }
5343  else if( strncmp(str, "power", 5) == 0 )
5344  {
5345  /* we have a string of the form "power(...,integer)" (thus, intpower)
5346  * first find the closing parenthesis, then the comma
5347  */
5348  const char* comma;
5349  int exponent;
5350 
5351  str += 5;
5352  SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
5353 
5354  SCIP_CALL( exprparseFindSeparatingComma(str+1, &comma, endptr - str - 1) );
5355 
5356  /* parse first argument [str+1..comma-1] */
5357  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, comma - str - 1, comma - 1, nvars, varnames,
5358  varnameslength, vartable, recursiondepth + 1) );
5359 
5360  ++comma;
5361  /* parse second argument [comma, endptr-1]: it needs to be an integer */
5362  while( comma < endptr && *comma == ' ' )
5363  ++comma;
5364  if( !isdigit((unsigned char)comma[0]) && !((comma[0] == '-' || comma[0] == '+') && isdigit((unsigned char)comma[1])) )
5365  {
5366  SCIPerrorMessage("error parsing integer exponent from <%s>\n", comma);
5367  }
5368  if( !SCIPstrToIntValue(comma, &exponent, &nonconstendptr) )
5369  {
5370  SCIPerrorMessage("error parsing integer from <%s>\n", comma);
5371  return SCIP_READERROR;
5372  }
5373 
5374  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_INTPOWER, arg1, exponent) );
5375 
5376  str = endptr + 1;
5377  }
5378  else if( strncmp(str, "realpower", 9) == 0 || strncmp(str, "signpower", 9) == 0 )
5379  {
5380  /* we have a string of the form "realpower(...,double)" or "signpower(...,double)"
5381  * first find the closing parenthesis, then the comma
5382  */
5383  const char* opname = str;
5384  const char* comma;
5385 
5386  str += 9;
5387  SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
5388 
5389  SCIP_CALL( exprparseFindSeparatingComma(str+1, &comma, endptr - str - 1) );
5390 
5391  /* parse first argument [str+1..comma-1] */
5392  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, comma - str - 1, comma - 1, nvars, varnames,
5393  varnameslength, vartable, recursiondepth + 1) );
5394 
5395  ++comma;
5396  /* parse second argument [comma, endptr-1]: it needs to be an number */
5397  while( comma < endptr && *comma == ' ' )
5398  ++comma;
5399  if( !isdigit((unsigned char)comma[0]) && !((comma[0] == '-' || comma[0] == '+') && isdigit((unsigned char)comma[1])) )
5400  {
5401  SCIPerrorMessage("error parsing number exponent from <%s>\n", comma);
5402  }
5403  if( !SCIPstrToRealValue(comma, &number, &nonconstendptr) )
5404  {
5405  SCIPerrorMessage("error parsing number from <%s>\n", comma);
5406  return SCIP_READERROR;
5407  }
5408 
5409  if( strncmp(opname, "realpower", 9) == 0 )
5410  {
5411  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_REALPOWER, arg1, number) );
5412  }
5413  else
5414  {
5415  assert(strncmp(opname, "signpower", 9) == 0);
5416  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_SIGNPOWER, arg1, number) );
5417  }
5418 
5419  str = endptr + 1;
5420  }
5421  else if( isalpha(*str) || *str == '_' || *str == '#' )
5422  {
5423  /* check for a variable, that was not recognized earlier because somebody omitted the '<' and '>' we need for
5424  * SCIPparseVarName, making everyones life harder;
5425  * we allow only variable names starting with a character or underscore here
5426  */
5427  const char* varnamestartptr = str;
5428 
5429  /* allow only variable names containing characters, digits, hash marks, and underscores here */
5430  while( isalnum(str[0]) || str[0] == '_' || str[0] == '#' )
5431  ++str;
5432 
5433  SCIP_CALL( exprparseReadVariable(blkmem, &varnamestartptr, expr, nvars, varnames, varnameslength,
5434  vartable, 1.0, str) );
5435  }
5436  else
5437  {
5438  SCIPerrorMessage("parsing of invalid expression %.*s.\n", (int) (lastchar - str + 1), str);
5439  return SCIP_READERROR;
5440  }
5441 
5442  /* if we are one char behind lastchar, we are done */
5443  if( str == lastchar + 1)
5444  {
5445  SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
5446  SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
5447  SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );
5448 
5449  return SCIP_OKAY;
5450  }
5451 
5452  /* check if we are still in bounds */
5453  if( str > lastchar + 1)
5454  {
5455  SCIPerrorMessage("error finding first expression in \"%.*s\" took us outside of given subexpression length\n", length, strstart);
5456  return SCIP_READERROR;
5457  }
5458 
5459  /* ignore whitespace */
5460  while( isspace((unsigned char)*str) && str != lastchar + 1 )
5461  ++str;
5462 
5463  /* maybe now we're done? */
5464  if( str >= lastchar + 1)
5465  {
5466  SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
5467  SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
5468  SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );
5469 
5470  return SCIP_OKAY;
5471  }
5472 
5473  if( str[0] == '^' )
5474  {
5475  /* a '^' behind the found expression indicates a power */
5476  SCIP_Real constant;
5477 
5478  arg1 = *expr;
5479  ++str;
5480  while( isspace((unsigned char)*str) && str != lastchar + 1 )
5481  ++str;
5482 
5483  if( isdigit((unsigned char)str[0]) || ((str[0] == '-' || str[0] == '+') && isdigit((unsigned char)str[1])) )
5484  {
5485  /* there is a number coming */
5486  if( !SCIPstrToRealValue(str, &number, &nonconstendptr) )
5487  {
5488  SCIPerrorMessage("error parsing number from <%s>\n", str);
5489  return SCIP_READERROR;
5490  }
5491 
5492  SCIP_CALL( SCIPexprCreate(blkmem, &arg2, SCIP_EXPR_CONST, number) );
5493  str = nonconstendptr;
5494 
5495  constant = SCIPexprGetOpReal(arg2);
5496  SCIPexprFreeDeep(blkmem, &arg2);
5497 
5498  /* expr^number is intpower or realpower */
5499  if( EPSISINT(constant, 0.0) ) /*lint !e835*/
5500  {
5501  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_INTPOWER, arg1, (int)constant) );
5502  }
5503  else
5504  {
5505  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_REALPOWER, arg1, constant) );
5506  }
5507  }
5508  else if( str[0] == '(' )
5509  {
5510  /* we use exprParse to evaluate the exponent */
5511 
5512  SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
5513  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, str + 1, endptr - str - 1, endptr -1, nvars, varnames,
5514  varnameslength, vartable, recursiondepth + 1) );
5515 
5516  if( SCIPexprGetOperator(arg2) != SCIP_EXPR_CONST )
5517  {
5518  /* reformulate arg1^arg2 as exp(arg2 * log(arg1)) */
5519  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_LOG, arg1) );
5520  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_MUL, *expr, arg2) );
5521  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_EXP, *expr) );
5522  }
5523  else
5524  {
5525  /* expr^number is intpower or realpower */
5526  constant = SCIPexprGetOpReal(arg2);
5527  SCIPexprFreeDeep(blkmem, &arg2);
5528  if( EPSISINT(constant, 0.0) ) /*lint !e835*/
5529  {
5530  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_INTPOWER, arg1, (int)constant) );
5531  }
5532  else
5533  {
5534  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_REALPOWER, arg1, constant) );
5535  }
5536  }
5537  str = endptr + 1;
5538  }
5539  else
5540  {
5541  SCIPerrorMessage("unexpected string following ^ in %.*s\n", length, str);
5542  return SCIP_READERROR;
5543  }
5544 
5545  /* ignore whitespace */
5546  while( isspace((unsigned char)*str) && str != lastchar + 1 )
5547  ++str;
5548  }
5549 
5550  /* check for a two argument operator that is not a multiplication */
5551  if( str <= lastchar && (str[0] == '+' || str[0] == '-' || str[0] == '/') )
5552  {
5553  char op;
5554 
5555  op = str[0];
5556  arg1 = *expr;
5557 
5558  /* step forward over the operator to go to the beginning of the second argument */
5559  ++str;
5560 
5561  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, str, (int) (lastchar - str + 1), lastchar, nvars, varnames,
5562  varnameslength, vartable, recursiondepth + 1) );
5563  str = lastchar + 1;
5564 
5565  /* make new expression from two arguments */
5566  if( op == '+')
5567  {
5568  SCIP_CALL( SCIPexprAdd(blkmem, expr, 1.0, arg1, 1.0, arg2, 0.0) );
5569  }
5570  else if( op == '-')
5571  {
5572  SCIP_CALL( SCIPexprAdd(blkmem, expr, 1.0, arg1, -1.0, arg2, 0.0) );
5573  }
5574  else if( op == '*' )
5575  {
5576  if( SCIPexprGetOperator(arg1) == SCIP_EXPR_CONST )
5577  {
5578  SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg2, SCIPexprGetOpReal(arg1)) );
5579  }
5580  else if( SCIPexprGetOperator(arg2) == SCIP_EXPR_CONST )
5581  {
5582  SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg1, SCIPexprGetOpReal(arg2)) );
5583  }
5584  else
5585  {
5586  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_MUL, arg1, arg2) );
5587  }
5588  }
5589  else
5590  {
5591  assert(op == '/');
5592 
5593  if( SCIPexprGetOperator(arg2) == SCIP_EXPR_CONST )
5594  {
5595  SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg1, 1.0 / SCIPexprGetOpReal(arg2)) );
5596  SCIPexprFreeShallow(blkmem, &arg2);
5597  }
5598  else
5599  {
5600  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_DIV, arg1, arg2) );
5601  }
5602  }
5603  }
5604 
5605  /* ignore whitespace */
5606  while( isspace((unsigned char)*str) )
5607  ++str;
5608 
5609  /* we are either done or we have a multiplication? */
5610  if( str >= lastchar + 1 )
5611  {
5612  SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
5613  SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
5614  SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );
5615 
5616  return SCIP_OKAY;
5617  }
5618 
5619  /* if there is a part of the string left to be parsed, we assume that this as a multiplication */
5620  arg1 = *expr;
5621 
5622  /* stepping over multiplication operator if needed */
5623  if( str[0] == '*' )
5624  {
5625  ++str;
5626  }
5627  else if( str[0] != '(' )
5628  {
5629  SCIPdebugMessage("No operator found, assuming a multiplication before %.*s\n", (int) (lastchar - str + 1), str);
5630  }
5631 
5632  SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, str, (int) (lastchar - str + 1), lastchar, nvars, varnames,
5633  varnameslength, vartable, recursiondepth + 1) );
5634 
5635  if( SCIPexprGetOperator(arg1) == SCIP_EXPR_CONST )
5636  {
5637  SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg2, SCIPexprGetOpReal(arg1)) );
5638  SCIPexprFreeDeep(blkmem, &arg1);
5639  }
5640  else if( SCIPexprGetOperator(arg2) == SCIP_EXPR_CONST )
5641  {
5642  SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg1, SCIPexprGetOpReal(arg2)) );
5643  SCIPexprFreeDeep(blkmem, &arg2);
5644  }
5645  else
5646  {
5647  SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_MUL, arg1, arg2) );
5648  }
5649 
5650  SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
5651  SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
5652  SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );
5653 
5654  return SCIP_OKAY;
5655 }
5656 
5657 /**@} */
5658 
5659 /**@name Expression methods */
5660 /**@{ */
5661 
5662 /* In debug mode, the following methods are implemented as function calls to ensure
5663  * type validity.
5664  * In optimized mode, the methods are implemented as defines to improve performance.
5665  * However, we want to have them in the library anyways, so we have to undef the defines.
5666  */
5667 
5668 #undef SCIPexprGetOperator
5669 #undef SCIPexprGetNChildren
5670 #undef SCIPexprGetChildren
5671 #undef SCIPexprGetOpIndex
5672 #undef SCIPexprGetOpReal
5673 #undef SCIPexprGetOpData
5674 #undef SCIPexprGetRealPowerExponent
5675 #undef SCIPexprGetIntPowerExponent
5676 #undef SCIPexprGetSignPowerExponent
5677 #undef SCIPexprGetLinearCoefs
5678 #undef SCIPexprGetLinearConstant
5679 #undef SCIPexprGetQuadElements
5680 #undef SCIPexprGetQuadConstant
5681 #undef SCIPexprGetQuadLinearCoefs
5682 #undef SCIPexprGetNQuadElements
5683 #undef SCIPexprGetMonomials
5684 #undef SCIPexprGetNMonomials
5685 #undef SCIPexprGetPolynomialConstant
5686 #undef SCIPexprGetMonomialCoef
5687 #undef SCIPexprGetMonomialNFactors
5688 #undef SCIPexprGetMonomialChildIndices
5689 #undef SCIPexprGetMonomialExponents
5690 #undef SCIPexprGetUserData
5691 #undef SCIPexprHasUserEstimator
5692 #undef SCIPexprGetUserEvalCapability
5693 
5694 /** gives operator of expression */
5696  SCIP_EXPR* expr /**< expression */
5697  )
5698 {
5699  assert(expr != NULL);
5700 
5701  return expr->op;
5702 }
5703 
5704 /** gives number of children of an expression */
5706  SCIP_EXPR* expr /**< expression */
5707  )
5708 {
5709  assert(expr != NULL);
5710 
5711  return expr->nchildren;
5712 }
5713 
5714 /** gives pointer to array with children of an expression */
5716  SCIP_EXPR* expr /**< expression */
5717  )
5718 {
5719  assert(expr != NULL);
5720 
5721  return expr->children;
5722 }
5723 
5724 /** gives index belonging to a SCIP_EXPR_VARIDX or SCIP_EXPR_PARAM operand */
5726  SCIP_EXPR* expr /**< expression */
5727  )
5728 {
5729  assert(expr != NULL);
5730  assert(expr->op == SCIP_EXPR_VARIDX || expr->op == SCIP_EXPR_PARAM);
5731 
5732  return expr->data.intval;
5733 }
5734 
5735 /** gives real belonging to a SCIP_EXPR_CONST operand */
5737  SCIP_EXPR* expr /**< expression */
5738  )
5739 {
5740  assert(expr != NULL);
5741  assert(expr->op == SCIP_EXPR_CONST);
5742 
5743  return expr->data.dbl;
5744 }
5745 
5746 /** gives void* belonging to a complex operand */
5748  SCIP_EXPR* expr /**< expression */
5749  )
5750 {
5751  assert(expr != NULL);
5752  assert(expr->op >= SCIP_EXPR_SUM); /* only complex operands store their data as void* */
5753 
5754  return expr->data.data;
5755 }
5756 
5757 /** gives exponent belonging to a SCIP_EXPR_REALPOWER expression */
5759  SCIP_EXPR* expr /**< expression */
5760  )
5761 {
5762  assert(expr != NULL);
5763  assert(expr->op == SCIP_EXPR_REALPOWER);
5764 
5765  return expr->data.dbl;
5766 }
5767 
5768 /** gives exponent belonging to a SCIP_EXPR_INTPOWER expression */
5770  SCIP_EXPR* expr /**< expression */
5771  )
5772 {
5773  assert(expr != NULL);
5774  assert(expr->op == SCIP_EXPR_INTPOWER);
5775 
5776  return expr->data.intval;
5777 }
5778 
5779 /** gives exponent belonging to a SCIP_EXPR_SIGNPOWER expression */
5781  SCIP_EXPR* expr /**< expression */
5782  )
5783 {
5784  assert(expr != NULL);
5785  assert(expr->op == SCIP_EXPR_SIGNPOWER);
5786 
5787  return expr->data.dbl;
5788 }
5789 
5790 /** gives linear coefficients belonging to a SCIP_EXPR_LINEAR expression */
5792  SCIP_EXPR* expr /**< expression */
5793  )
5794 {
5795  assert(expr != NULL);
5796  assert(expr->op == SCIP_EXPR_LINEAR);
5797  assert(expr->data.data != NULL);
5798 
5799  /* the coefficients are stored in the first nchildren elements of the array stored as expression data */
5800  return (SCIP_Real*)expr->data.data;
5801 }
5802 
5803 /** gives constant belonging to a SCIP_EXPR_LINEAR expression */
5805  SCIP_EXPR* expr /**< expression */
5806  )
5807 {
5808  assert(expr != NULL);
5809  assert(expr->op == SCIP_EXPR_LINEAR);
5810  assert(expr->data.data != NULL);
5811 
5812  /* the constant is stored in the nchildren's element of the array stored as expression data */
5813  return ((SCIP_Real*)expr->data.data)[expr->nchildren];
5814 }
5815 
5816 /** gives quadratic elements belonging to a SCIP_EXPR_QUADRATIC expression */
5818  SCIP_EXPR* expr /**< quadratic expression */
5819  )
5820 {
5821  assert(expr != NULL);
5822  assert(expr->op == SCIP_EXPR_QUADRATIC);
5823  assert(expr->data.data != NULL);
5824 
5825  return ((SCIP_EXPRDATA_QUADRATIC*)expr->data.data)->quadelems;
5826 }
5827 
5828 /** gives constant belonging to a SCIP_EXPR_QUADRATIC expression */
5830  SCIP_EXPR* expr /**< quadratic expression */
5831  )
5832 {
5833  assert(expr != NULL);
5834  assert(expr->op == SCIP_EXPR_QUADRATIC);
5835  assert(expr->data.data != NULL);
5836 
5837  return ((SCIP_EXPRDATA_QUADRATIC*)expr->data.data)->constant;
5838 }
5839 
5840 /** gives linear coefficients belonging to a SCIP_EXPR_QUADRATIC expression
5841  * can be NULL if all coefficients are 0.0 */
5843  SCIP_EXPR* expr /**< quadratic expression */
5844  )
5845 {
5846  assert(expr != NULL);
5847  assert(expr->op == SCIP_EXPR_QUADRATIC);
5848  assert(expr->data.data != NULL);
5849 
5850  return ((SCIP_EXPRDATA_QUADRATIC*)expr->data.data)->lincoefs;
5851 }
5852 
5853 /** gives number of quadratic elements belonging to a SCIP_EXPR_QUADRATIC expression */
5855