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