expr.c

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the program and library             */
/*         SCIP --- Solving Constraint Integer Programs                      */
/*                                                                           */
/*    Copyright (C) 2002-2019 Konrad-Zuse-Zentrum                            */
/*                            fuer Informationstechnik Berlin                */
/*                                                                           */
/*  SCIP is distributed under the terms of the ZIB Academic License.         */
/*                                                                           */
/*  You should have received a copy of the ZIB Academic License              */
/*  along with SCIP; see the file COPYING. If not visit scip.zib.de.         */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/**@file   nlpi/expr.c
 * @brief  methods for expressions, expression trees, expression graphs, and related
 * @author Stefan Vigerske
 * @author Thorsten Gellermann
 * @author Ingmar Vierhaus (exprparse)
 */

/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/

#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <ctype.h>

#include "nlpi/pub_expr.h"
#include "nlpi/struct_expr.h"
#include "nlpi/exprinterpret.h"

#include "scip/intervalarith.h"
#include "scip/pub_misc.h"
#include "scip/misc.h"
#include "scip/pub_message.h"


#define SCIP_EXPRESSION_MAXCHILDEST 16       /**< estimate on maximal number of children */

/** sign of a value (-1 or +1)
 *
 *  0.0 has sign +1
 */
#define SIGN(x) ((x) >= 0.0 ? 1.0 : -1.0)

/** ensures that a block memory array has at least a given size
 *
 *  if cursize is 0, then *array1 can be NULL
 */
#define ensureBlockMemoryArraySize(blkmem, array1, cursize, minsize)    \
   do {                                                                 \
      int __newsize;                                                    \
      assert((blkmem)  != NULL);                                        \
      if( *(cursize) >= (minsize) )                                     \
         break;                                                         \
      __newsize = calcGrowSize(minsize);                                \
      assert(__newsize >= (minsize));                                   \
      SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, array1, *(cursize), __newsize) ); \
      *(cursize) = __newsize;                                           \
   } while( FALSE )

#ifdef SCIP_DISABLED_CODE /* this macro is currently not used, which offends lint, so disable it */
/** ensures that two block memory arrays have at least a given size
 *
 *  if cursize is 0, then arrays can be NULL
 */
#define ensureBlockMemoryArraySize2(blkmem, array1, array2, cursize, minsize) \
   do {                                                                 \
      int __newsize;                                                    \
      assert((blkmem)  != NULL);                                        \
      if( *(cursize) >= (minsize) )                                     \
         break;                                                         \
      __newsize = calcGrowSize(minsize);                                \
      assert(__newsize >= (minsize));                                   \
      SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, array1, *(cursize), __newsize) ); \
      SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, array2, *(cursize), __newsize) ); \
      *(cursize) = __newsize;                                           \
   } while( FALSE )
#endif

/** ensures that three block memory arrays have at least a given size
 *
 *  if cursize is 0, then arrays can be NULL
 */
#define ensureBlockMemoryArraySize3(blkmem, array1, array2, array3, cursize, minsize) \
   do {                                                                 \
      int __newsize;                                                    \
      assert((blkmem)  != NULL);                                        \
      if( *(cursize) >= (minsize) )                                     \
         break;                                                         \
      __newsize = calcGrowSize(minsize);                                \
      assert(__newsize >= (minsize));                                   \
      SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, array1, *(cursize), __newsize) ); \
      SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, array2, *(cursize), __newsize) ); \
      SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, array3, *(cursize), __newsize) ); \
      *(cursize) = __newsize;                                           \
   } while( FALSE )

/**@name Miscellaneous private methods */
/**@{ */

/** calculate memory size for dynamically allocated arrays (copied from scip/set.c) */
static
int calcGrowSize(
   int                   num                 /**< minimum number of entries to store */
   )
{
   int size;

   /* calculate the size with this loop, such that the resulting numbers are always the same (-> block memory) */
   size = 4;
   while( size < num )
      size = (int)(1.2 * size + 4);

   return size;
}

/** expression graph nodes comparison to use in sorting methods
 *
 * The nodes need to have been added to the expression graph (depth,pos >= 0).
 * The better node is the one with the lower depth and lower position, if depth is equal.
 */
static
SCIP_DECL_SORTPTRCOMP(exprgraphnodecomp)
{
   SCIP_EXPRGRAPHNODE* node1 = (SCIP_EXPRGRAPHNODE*)elem1;
   SCIP_EXPRGRAPHNODE* node2 = (SCIP_EXPRGRAPHNODE*)elem2;

   assert(node1 != NULL);
   assert(node2 != NULL);
   assert(node1->depth >= 0);
   assert(node1->pos >= 0);
   assert(node2->depth >= 0);
   assert(node2->pos >= 0);

   if( node1->depth != node2->depth )
      return node1->depth - node2->depth;

   /* there should be no two nodes on the same position */
   assert((node1->pos != node2->pos) || (node1 == node2));

   return node1->pos - node2->pos;
}

/** 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) */
static
SCIP_Bool isLbBetter(
   SCIP_Real             minstrength,        /**< minimal relative improvement required to be a better bound */
   SCIP_Real             newlb,              /**< new lower bound */
   SCIP_Real             oldlb,              /**< old lower bound */
   SCIP_Real             oldub               /**< old upper bound */
   )
{
   SCIP_Real eps;

   /* nothing can be tighter than an empty interval */
   if( oldlb > oldub )
      return FALSE;

   eps = REALABS(oldlb);
   eps = MIN(oldub - oldlb, eps);
   return EPSGT(newlb, oldlb, minstrength * MAX(eps, 1e-3));
}

/** 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) */
static
SCIP_Bool isUbBetter(
   SCIP_Real             minstrength,        /**< minimal relative improvement required to be a better bound */
   SCIP_Real             newub,              /**< new upper bound */
   SCIP_Real             oldlb,              /**< old lower bound */
   SCIP_Real             oldub               /**< old upper bound */
   )
{
   SCIP_Real eps;

   /* nothing can be tighter than an empty interval */
   if( oldlb > oldub )
      return FALSE;

   eps = REALABS(oldub);
   eps = MIN(oldub - oldlb, eps);
   return EPSLT(newub, oldub, minstrength * MAX(eps, 1e-3));
}

/**@} */

/**@name Expression curvature methods */
/**@{ */

/** curvature names as strings */
static
const char* curvnames[4] =
   {
      "unknown",
      "convex",
      "concave",
      "linear"
   };

#undef SCIPexprcurvAdd

/** gives curvature for a sum of two functions with given curvature */
SCIP_EXPRCURV SCIPexprcurvAdd(
   SCIP_EXPRCURV         curv1,              /**< curvature of first summand */
   SCIP_EXPRCURV         curv2               /**< curvature of second summand */
   )
{
   return (SCIP_EXPRCURV) (curv1 & curv2);
}

/** gives the curvature for the negation of a function with given curvature */
SCIP_EXPRCURV SCIPexprcurvNegate(
   SCIP_EXPRCURV         curvature           /**< curvature of function */
   )
{
   switch( curvature )
   {
   case SCIP_EXPRCURV_CONCAVE:
      return SCIP_EXPRCURV_CONVEX;

   case SCIP_EXPRCURV_CONVEX:
      return SCIP_EXPRCURV_CONCAVE;

   case SCIP_EXPRCURV_LINEAR:
   case SCIP_EXPRCURV_UNKNOWN:
      /* can return curvature, do this below */
      break;

   default:
      SCIPerrorMessage("unknown curvature status.\n");
      SCIPABORT();
   }

   return curvature;
}

/** gives curvature for a functions with given curvature multiplied by a constant factor */
SCIP_EXPRCURV SCIPexprcurvMultiply(
   SCIP_Real             factor,             /**< constant factor */
   SCIP_EXPRCURV         curvature           /**< curvature of other factor */
   )
{
   if( factor == 0.0 )
      return SCIP_EXPRCURV_LINEAR;
   if( factor > 0.0 )
      return curvature;
   return SCIPexprcurvNegate(curvature);
}

/** gives curvature for base^exponent for given bounds and curvature of base-function and constant exponent */
SCIP_EXPRCURV SCIPexprcurvPower(
   SCIP_INTERVAL         basebounds,         /**< bounds on base function */
   SCIP_EXPRCURV         basecurv,           /**< curvature of base function */
   SCIP_Real             exponent            /**< exponent */
   )
{
   SCIP_Bool expisint;

   assert(basebounds.inf <= basebounds.sup);

   if( exponent == 0.0 )
      return SCIP_EXPRCURV_LINEAR;

   if( exponent == 1.0 )
      return basecurv;

   expisint = EPSISINT(exponent, 0.0); /*lint !e835*/

   /* if exponent is fractional, then power is not defined for a negative base
    * thus, consider only positive part of basebounds
    */
   if( !expisint && basebounds.inf < 0.0 )
   {
      basebounds.inf = 0.0;
      if( basebounds.sup < 0.0 )
         return SCIP_EXPRCURV_LINEAR;
   }

   /* if basebounds contains 0.0, consider negative and positive interval separately, if possible */
   if( basebounds.inf < 0.0 && basebounds.sup > 0.0 )
   {
      SCIP_INTERVAL leftbounds;
      SCIP_INTERVAL rightbounds;

      /* 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 */
      if( exponent < 0.0 )
         return SCIP_EXPRCURV_UNKNOWN;

      SCIPintervalSetBounds(&leftbounds,  basebounds.inf, 0.0);
      SCIPintervalSetBounds(&rightbounds, 0.0, basebounds.sup);

      return (SCIP_EXPRCURV) (SCIPexprcurvPower(leftbounds,  basecurv, exponent) & SCIPexprcurvPower(rightbounds, basecurv, exponent));
   }
   assert(basebounds.inf >= 0.0 || basebounds.sup <= 0.0);

   /* (base^exponent)'' = exponent * ( (exponent-1) base^(exponent-2) (base')^2 + base^(exponent-1) base'' )
    *
    * if base'' is positive, i.e., base is convex, then
    * - for base > 0.0 and exponent > 1.0, the second deriv. is positive -> convex
    * - for base < 0.0 and exponent > 1.0, we can't say (first and second summand opposite signs)
    * - for base > 0.0 and 0.0 < exponent < 1.0, we can't say (first sommand negative, second summand positive)
    * - for base > 0.0 and exponent < 0.0, we can't say (first and second summand opposite signs)
    * - for base < 0.0 and exponent < 0.0 and even, the second deriv. is positive -> convex
    * - for base < 0.0 and exponent < 0.0 and odd, the second deriv. is negative -> concave
    *
    * if base'' is negative, i.e., base is concave, then
    * - for base > 0.0 and exponent > 1.0, we can't say (first summand positive, second summand negative)
    * - for base < 0.0 and exponent > 1.0 and even, the second deriv. is positive -> convex
    * - for base < 0.0 and exponent > 1.0 and odd, the second deriv. is negative -> concave
    * - for base > 0.0 and 0.0 < exponent < 1.0, the second deriv. is negative -> concave
    * - for base > 0.0 and exponent < 0.0, the second deriv. is positive -> convex
    * - for base < 0.0 and exponent < 0.0, we can't say (first and second summand opposite signs)
    *
    * if base'' is zero, i.e., base is linear, then
    *   (base^exponent)'' = exponent * (exponent-1) base^(exponent-2) (base')^2
    * - just multiply signs
    */

   if( basecurv == SCIP_EXPRCURV_LINEAR )
   {
      SCIP_Real sign;

      /* base^(exponent-2) is negative, if base < 0.0 and exponent is odd */
      sign = exponent * (exponent - 1.0);
      assert(basebounds.inf >= 0.0 || expisint);
      if( basebounds.inf < 0.0 && ((int)exponent)%2 != 0 )
         sign *= -1.0;
      assert(sign != 0.0);

      return sign > 0.0 ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
   }

   if( basecurv == SCIP_EXPRCURV_CONVEX )
   {
      if( basebounds.sup <= 0.0 && exponent < 0.0 && expisint )
         return ((int)exponent)%2 == 0 ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
      if( basebounds.inf >= 0.0 && exponent > 1.0 )
         return SCIP_EXPRCURV_CONVEX ;
      return SCIP_EXPRCURV_UNKNOWN;
   }

   if( basecurv == SCIP_EXPRCURV_CONCAVE )
   {
      if( basebounds.sup <= 0.0 && exponent > 1.0 && expisint )
         return ((int)exponent)%2 == 0 ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
      if( basebounds.inf >= 0.0 && exponent < 1.0 )
         return exponent < 0.0 ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
      return SCIP_EXPRCURV_UNKNOWN;
   }

   return SCIP_EXPRCURV_UNKNOWN;
}

/** gives curvature for a monomial with given curvatures and bounds for each factor
 *
 *  See Maranas and Floudas, Finding All Solutions of Nonlinearly Constrained Systems of Equations, JOGO 7, 1995
 *  for the categorization in the case that all factors are linear.
 */
SCIP_EXPRCURV SCIPexprcurvMonomial(
   int                   nfactors,           /**< number of factors in monomial */
   SCIP_Real*            exponents,          /**< exponents in monomial, or NULL if all 1.0 */
   int*                  factoridxs,         /**< indices of factors (but not exponents), or NULL if identity mapping */
   SCIP_EXPRCURV*        factorcurv,         /**< curvature of each factor */
   SCIP_INTERVAL*        factorbounds        /**< bounds of each factor */
   )
{
   SCIP_Real mult;
   SCIP_Real e;
   SCIP_EXPRCURV curv;
   SCIP_EXPRCURV fcurv;
   int nnegative;
   int npositive;
   SCIP_Real sum;
   SCIP_Bool expcurvpos;
   SCIP_Bool expcurvneg;
   int j;
   int f;

   assert(nfactors >= 0);
   assert(factorcurv   != NULL || nfactors == 0);
   assert(factorbounds != NULL || nfactors == 0);

   if( nfactors == 0 )
      return SCIP_EXPRCURV_LINEAR;

   if( nfactors == 1 )
   {
      f = factoridxs != NULL ? factoridxs[0] : 0;
      e = exponents != NULL ? exponents[0] : 1.0;
      /* SCIPdebugMessage("monomial [%g,%g]^%g is %s\n",
         factorbounds[f].inf, factorbounds[f].sup, e,
         SCIPexprcurvGetName(SCIPexprcurvPower(factorbounds[f], factorcurv[f], e))); */
      return SCIPexprcurvPower(factorbounds[f], factorcurv[f], e);  /*lint !e613*/
   }

   mult = 1.0;

   nnegative = 0; /* number of negative exponents */
   npositive = 0; /* number of positive exponents */
   sum = 0.0;     /* sum of exponents */
   expcurvpos = TRUE; /* whether exp_j * f_j''(x) >= 0 for all factors (assuming f_j >= 0) */
   expcurvneg = TRUE; /* whether exp_j * f_j''(x) <= 0 for all factors (assuming f_j >= 0) */

   for( j = 0; j < nfactors; ++j )
   {
      f = factoridxs != NULL ? factoridxs[j] : j;
      if( factorcurv[f] == SCIP_EXPRCURV_UNKNOWN ) /*lint !e613*/
         return SCIP_EXPRCURV_UNKNOWN;
      if( factorbounds[f].inf < 0.0 && factorbounds[f].sup > 0.0 )  /*lint !e613*/
         return SCIP_EXPRCURV_UNKNOWN;

      e = exponents != NULL ? exponents[j] : 1.0;
      if( e < 0.0 )
         ++nnegative;
      else
         ++npositive;
      sum += e;

      if( factorbounds[f].inf < 0.0 )  /*lint !e613*/
      {
         /* if argument is negative, then exponent should be integer */
         assert(EPSISINT(e, 0.0));  /*lint !e835*/

         /* flip j'th argument: (f_j)^(exp_j) = (-1)^(exp_j) (-f_j)^(exp_j) */

         /* -f_j has negated curvature of f_j */
         fcurv = SCIPexprcurvNegate(factorcurv[f]);  /*lint !e613*/

         /* negate monomial, if exponent is odd, i.e., (-1)^(exp_j) = -1 */
         if( (int)e % 2 != 0 )
            mult *= -1.0;
      }
      else
      {
         fcurv = factorcurv[f];  /*lint !e613*/
      }

      /* check if exp_j * fcurv is convex (>= 0) and/or concave */
      fcurv = SCIPexprcurvMultiply(e, fcurv);
      if( !(fcurv & SCIP_EXPRCURV_CONVEX) )
         expcurvpos = FALSE;
      if( !(fcurv & SCIP_EXPRCURV_CONCAVE) )
         expcurvneg = FALSE;
   }

   /* if all factors are linear, then a product f_j^exp_j with f_j >= 0 is convex if
    * - all exponents are negative, or
    * - 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
    * further, the product is concave if
    * - all exponents are positive and the sum of exponents is <= 1.0
    *
    * if factors are nonlinear, then we require additionally, that for convexity
    * - each factor is convex if exp_j >= 0, or concave if exp_j <= 0, i.e., exp_j*f_j'' >= 0
    * and for concavity, we require that
    * - all factors are concave, i.e., exp_j*f_j'' <= 0
    */

   if( nnegative == nfactors && expcurvpos )
      curv = SCIP_EXPRCURV_CONVEX;
   else if( nnegative == nfactors-1 && EPSGE(sum, 1.0, 1e-9) && expcurvpos )
      curv = SCIP_EXPRCURV_CONVEX;
   else if( npositive == nfactors && EPSLE(sum, 1.0, 1e-9) && expcurvneg )
      curv = SCIP_EXPRCURV_CONCAVE;
   else
      curv = SCIP_EXPRCURV_UNKNOWN;
   curv = SCIPexprcurvMultiply(mult, curv);

   return curv;
}

/** gives name as string for a curvature */
const char* SCIPexprcurvGetName(
   SCIP_EXPRCURV         curv                /**< curvature */
   )
{
   assert(curv <= SCIP_EXPRCURV_LINEAR);  /*lint !e685*/

   return curvnames[curv];
}

/**@} */

/**@name Quadratic expression data private methods */
/**@{ */

/** creates SCIP_EXPRDATA_QUADRATIC data structure from given quadratic elements */
static
SCIP_RETCODE quadraticdataCreate(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPRDATA_QUADRATIC** quadraticdata,  /**< buffer to store pointer to quadratic data */
   SCIP_Real             constant,           /**< constant */
   int                   nchildren,          /**< number of children */
   SCIP_Real*            lincoefs,           /**< linear coefficients of children, or NULL if all 0.0 */
   int                   nquadelems,         /**< number of quadratic elements */
   SCIP_QUADELEM*        quadelems           /**< quadratic elements */
   )
{
   assert(blkmem != NULL);
   assert(quadraticdata != NULL);
   assert(quadelems != NULL || nquadelems == 0);
   assert(nchildren >= 0);

   SCIP_ALLOC( BMSallocBlockMemory(blkmem, quadraticdata) );

   (*quadraticdata)->constant   = constant;
   (*quadraticdata)->lincoefs   = NULL;
   (*quadraticdata)->nquadelems = nquadelems;
   (*quadraticdata)->quadelems  = NULL;
   (*quadraticdata)->sorted     = (nquadelems <= 1);

   if( lincoefs != NULL )
   {
      SCIP_ALLOC( BMSduplicateBlockMemoryArray(blkmem, &(*quadraticdata)->lincoefs, lincoefs, nchildren) );
   }

   if( nquadelems > 0 )
   {
      SCIP_ALLOC( BMSduplicateBlockMemoryArray(blkmem, &(*quadraticdata)->quadelems, quadelems, nquadelems) );
   }

   return SCIP_OKAY;
}

/** sorts quadratic elements in a SCIP_EXPRDATA_QUADRATIC data structure */
static
void quadraticdataSort(
   SCIP_EXPRDATA_QUADRATIC* quadraticdata    /**< quadratic data */
   )
{
   assert(quadraticdata != NULL);

   if( quadraticdata->sorted )
   {
#ifndef NDEBUG
      int i;
      for( i = 1; i < quadraticdata->nquadelems; ++i )
      {
         assert(quadraticdata->quadelems[i].idx1 <= quadraticdata->quadelems[i].idx2);
         assert(quadraticdata->quadelems[i-1].idx1 <= quadraticdata->quadelems[i].idx1);
         assert(quadraticdata->quadelems[i-1].idx1 < quadraticdata->quadelems[i].idx1 ||
            quadraticdata->quadelems[i-1].idx2 <= quadraticdata->quadelems[i].idx2);
      }
#endif
      return;
   }

   if( quadraticdata->nquadelems > 0 )
      SCIPquadelemSort(quadraticdata->quadelems, quadraticdata->nquadelems);

   quadraticdata->sorted = TRUE;
}

/**@} */

/**@name Polynomial expression data private methods */
/**@{ */

/** compares two monomials
 *
 *  gives 0 if monomials are equal */
static
SCIP_DECL_SORTPTRCOMP(monomialdataCompare)
{
   SCIP_EXPRDATA_MONOMIAL* monomial1;
   SCIP_EXPRDATA_MONOMIAL* monomial2;

   int i;

   assert(elem1 != NULL);
   assert(elem2 != NULL);

   monomial1 = (SCIP_EXPRDATA_MONOMIAL*)elem1;
   monomial2 = (SCIP_EXPRDATA_MONOMIAL*)elem2;

   /* make sure, both monomials are equal */
   SCIPexprSortMonomialFactors(monomial1);
   SCIPexprSortMonomialFactors(monomial2);

   /* for the first factor where both monomials differ,
    * we return either the difference in the child indices, if children are different
    * or the sign of the difference in the exponents
    */
   for( i = 0; i < monomial1->nfactors && i < monomial2->nfactors; ++i )
   {
      if( monomial1->childidxs[i] != monomial2->childidxs[i] )
         return monomial1->childidxs[i] - monomial2->childidxs[i];
      if( monomial1->exponents[i] > monomial2->exponents[i] )
         return 1;
      else if( monomial1->exponents[i] < monomial2->exponents[i] )
         return -1;
   }

   /* if the factors of one monomial are a proper subset of the factors of the other monomial,
    * we return the difference in the number of monomials
    */
   return monomial1->nfactors - monomial2->nfactors;
}

/** ensures that the factors arrays of a monomial have at least a given size */
static
SCIP_RETCODE monomialdataEnsureFactorsSize(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPRDATA_MONOMIAL*  monomialdata,    /**< monomial data */
   int                   minsize             /**< minimal size of factors arrays */
   )
{
   assert(blkmem != NULL);
   assert(monomialdata != NULL);

   if( minsize > monomialdata->factorssize )
   {
      int newsize;

      newsize = calcGrowSize(minsize);
      SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &monomialdata->childidxs, monomialdata->factorssize, newsize) );
      SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &monomialdata->exponents, monomialdata->factorssize, newsize) );
      monomialdata->factorssize = newsize;
   }
   assert(minsize <= monomialdata->factorssize);

   return SCIP_OKAY;
}

/** creates SCIP_EXPRDATA_POLYNOMIAL data structure from given monomials */
static
SCIP_RETCODE polynomialdataCreate(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPRDATA_POLYNOMIAL** polynomialdata,/**< buffer to store pointer to polynomial data */
   int                   nmonomials,         /**< number of monomials */
   SCIP_EXPRDATA_MONOMIAL** monomials,       /**< monomials */
   SCIP_Real             constant,           /**< constant part */
   SCIP_Bool             copymonomials       /**< whether to copy monomials, or copy only given pointers, in which case polynomialdata assumes ownership of monomial structure */
   )
{
   assert(blkmem != NULL);
   assert(polynomialdata != NULL);
   assert(monomials != NULL || nmonomials == 0);

   SCIP_ALLOC( BMSallocBlockMemory(blkmem, polynomialdata) );

   (*polynomialdata)->constant = constant;
   (*polynomialdata)->nmonomials  = nmonomials;
   (*polynomialdata)->monomialssize = nmonomials;
   (*polynomialdata)->monomials   = NULL;
   (*polynomialdata)->sorted   = (nmonomials <= 1);

   if( nmonomials > 0 )
   {
      int i;

      if( copymonomials )
      {
         SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &(*polynomialdata)->monomials, nmonomials) );

         for( i = 0; i < nmonomials; ++i )
         {
            assert(monomials[i] != NULL);  /*lint !e613*/
            SCIP_CALL( SCIPexprCreateMonomial(blkmem, &(*polynomialdata)->monomials[i],
                  monomials[i]->coef, monomials[i]->nfactors, monomials[i]->childidxs, monomials[i]->exponents) );  /*lint !e613*/
         }
      }
      else
      {
         SCIP_ALLOC( BMSduplicateBlockMemoryArray(blkmem, &(*polynomialdata)->monomials, monomials, nmonomials) );
      }
   }

   return SCIP_OKAY;
}

/** creates a copy of a SCIP_EXPRDATA_POLYNOMIAL data structure */
static
SCIP_RETCODE polynomialdataCopy(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPRDATA_POLYNOMIAL** polynomialdata,/**< buffer to store pointer to polynomial data */
   SCIP_EXPRDATA_POLYNOMIAL* sourcepolynomialdata /**< polynomial data to copy */
   )
{
   assert(blkmem != NULL);
   assert(polynomialdata != NULL);
   assert(sourcepolynomialdata != NULL);

   SCIP_ALLOC( BMSduplicateBlockMemory(blkmem, polynomialdata, sourcepolynomialdata) );

   (*polynomialdata)->monomialssize = sourcepolynomialdata->nmonomials;
   if( sourcepolynomialdata->nmonomials > 0 )
   {
      int i;

      SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &(*polynomialdata)->monomials, (*polynomialdata)->monomialssize) );

      for( i = 0; i < sourcepolynomialdata->nmonomials; ++i )
      {
         assert(sourcepolynomialdata->monomials[i] != NULL);  /*lint !e613*/
         SCIP_CALL( SCIPexprCreateMonomial(blkmem, &(*polynomialdata)->monomials[i], sourcepolynomialdata->monomials[i]->coef,
               sourcepolynomialdata->monomials[i]->nfactors, sourcepolynomialdata->monomials[i]->childidxs, sourcepolynomialdata->monomials[i]->exponents) );
         (*polynomialdata)->monomials[i]->sorted = sourcepolynomialdata->monomials[i]->sorted;
      }
   }
   else
   {
      (*polynomialdata)->monomials = NULL;
   }

   return SCIP_OKAY;
}

/** frees a SCIP_EXPRDATA_POLYNOMIAL data structure */
static
void polynomialdataFree(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPRDATA_POLYNOMIAL** polynomialdata /**< pointer to polynomial data to free */
   )
{
   assert(blkmem != NULL);
   assert(polynomialdata != NULL);
   assert(*polynomialdata != NULL);

   if( (*polynomialdata)->monomialssize > 0 )
   {
      int i;

      for( i = 0; i < (*polynomialdata)->nmonomials; ++i )
      {
         assert((*polynomialdata)->monomials[i] != NULL);
         SCIPexprFreeMonomial(blkmem, &(*polynomialdata)->monomials[i]);
         assert((*polynomialdata)->monomials[i] == NULL);
      }

      BMSfreeBlockMemoryArray(blkmem, &(*polynomialdata)->monomials, (*polynomialdata)->monomialssize);
   }
   assert((*polynomialdata)->monomials == NULL);

   BMSfreeBlockMemory(blkmem, polynomialdata);
}

/** ensures that the monomials array of a polynomial has at least a given size */
static
SCIP_RETCODE polynomialdataEnsureMonomialsSize(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
   int                   minsize             /**< minimal size of monomials array */
   )
{
   assert(blkmem != NULL);
   assert(polynomialdata != NULL);

   ensureBlockMemoryArraySize(blkmem, &polynomialdata->monomials, &polynomialdata->monomialssize, minsize);
   assert(minsize <= polynomialdata->monomialssize);

   return SCIP_OKAY;
}

/** adds an array of monomials to a polynomial */
static
SCIP_RETCODE polynomialdataAddMonomials(
   BMS_BLKMEM*           blkmem,             /**< block memory of expression */
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
   int                   nmonomials,         /**< number of monomials to add */
   SCIP_EXPRDATA_MONOMIAL** monomials,       /**< the monomials to add */
   SCIP_Bool             copymonomials       /**< whether to copy monomials or to assume ownership */
   )
{
   int i;

   assert(blkmem != NULL);
   assert(polynomialdata != NULL);
   assert(monomials != NULL || nmonomials == 0);

   if( nmonomials == 0 )
      return SCIP_OKAY;

   SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, polynomialdata->nmonomials + nmonomials) );
   assert(polynomialdata->monomialssize >= polynomialdata->nmonomials + nmonomials);

   if( copymonomials )
   {
      for( i = 0; i < nmonomials; ++i )
      {
         assert(monomials[i] != NULL);  /*lint !e613*/
         SCIP_CALL( SCIPexprCreateMonomial(blkmem, &polynomialdata->monomials[polynomialdata->nmonomials + i],
               monomials[i]->coef, monomials[i]->nfactors, monomials[i]->childidxs, monomials[i]->exponents) );  /*lint !e613*/
      }
   }
   else
   {
      BMScopyMemoryArray(&polynomialdata->monomials[polynomialdata->nmonomials], monomials, nmonomials);  /*lint !e866*/
   }
   polynomialdata->nmonomials += nmonomials;

   polynomialdata->sorted = (polynomialdata->nmonomials <= 1);

   return SCIP_OKAY;
}

/** ensures that monomials of a polynomial are sorted */
static
void polynomialdataSortMonomials(
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata  /**< polynomial expression */
   )
{
   assert(polynomialdata != NULL);

   if( polynomialdata->sorted )
   {
#ifndef NDEBUG
      int i;

      /* a polynom with more than one monoms can only be sorted if its monoms are sorted */
      for( i = 1; i < polynomialdata->nmonomials; ++i )
      {
         assert(polynomialdata->monomials[i-1]->sorted);
         assert(polynomialdata->monomials[i]->sorted);
         assert(monomialdataCompare(polynomialdata->monomials[i-1], polynomialdata->monomials[i]) <= 0);
      }
#endif
      return;
   }

   if( polynomialdata->nmonomials > 0 )
      SCIPsortPtr((void**)polynomialdata->monomials, monomialdataCompare, polynomialdata->nmonomials);

   polynomialdata->sorted = TRUE;
}

/** merges monomials that differ only in coefficient into a single monomial
 *
 *  Eliminates monomials with coefficient between -eps and eps.
 */
static
void polynomialdataMergeMonomials(
   BMS_BLKMEM*           blkmem,             /**< block memory */
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
   SCIP_Real             eps,                /**< threshold under which numbers are treat as zero */
   SCIP_Bool             mergefactors        /**< whether to merge factors in monomials too */
   )
{
   int i;
   int offset;
   int oldnfactors;

   assert(polynomialdata != NULL);
   assert(eps >= 0.0);

   polynomialdataSortMonomials(polynomialdata);

   /* merge monomials by adding their coefficients, eliminate monomials with no factors or zero coefficient*/
   offset = 0;
   i = 0;
   while( i + offset < polynomialdata->nmonomials )
   {
      if( offset > 0 )
      {
         assert(polynomialdata->monomials[i] == NULL);
         assert(polynomialdata->monomials[i+offset] != NULL);
         polynomialdata->monomials[i] = polynomialdata->monomials[i+offset];
#ifndef NDEBUG
         polynomialdata->monomials[i+offset] = NULL;
#endif
      }

      if( mergefactors )
      {
         oldnfactors = polynomialdata->monomials[i]->nfactors;
         SCIPexprMergeMonomialFactors(polynomialdata->monomials[i], eps);

         /* if monomial has changed, then we cannot assume anymore that polynomial is sorted */
         if( oldnfactors != polynomialdata->monomials[i]->nfactors )
            polynomialdata->sorted = FALSE;
      }

      while( i+offset+1 < polynomialdata->nmonomials )
      {
         assert(polynomialdata->monomials[i+offset+1] != NULL);
         if( mergefactors )
         {
            oldnfactors = polynomialdata->monomials[i+offset+1]->nfactors;
            SCIPexprMergeMonomialFactors(polynomialdata->monomials[i+offset+1], eps);

            /* if monomial has changed, then we cannot assume anymore that polynomial is sorted */
            if( oldnfactors != polynomialdata->monomials[i+offset+1]->nfactors )
               polynomialdata->sorted = FALSE;
         }
         if( monomialdataCompare((void*)polynomialdata->monomials[i], (void*)polynomialdata->monomials[i+offset+1]) != 0 )
            break;
         polynomialdata->monomials[i]->coef += polynomialdata->monomials[i+offset+1]->coef;
         SCIPexprFreeMonomial(blkmem, &polynomialdata->monomials[i+offset+1]);
         ++offset;
      }

      if( polynomialdata->monomials[i]->nfactors == 0 )
      {
         /* constant monomial */
         polynomialdata->constant += polynomialdata->monomials[i]->coef;
         SCIPexprFreeMonomial(blkmem, &polynomialdata->monomials[i]);
         ++offset;
         continue;
      }

      if( EPSZ(polynomialdata->monomials[i]->coef, eps) )
      {
         SCIPexprFreeMonomial(blkmem, &polynomialdata->monomials[i]);
         ++offset;
         continue;
      }

      ++i;
   }

#ifndef NDEBUG
   for( ; i < polynomialdata->nmonomials; ++i )
      assert(polynomialdata->monomials[i] == NULL);
#endif

   polynomialdata->nmonomials -= offset;

   if( EPSZ(polynomialdata->constant, eps) )
      polynomialdata->constant = 0.0;
}

/** multiplies each summand of a polynomial by a given constant */
static
void polynomialdataMultiplyByConstant(
   BMS_BLKMEM*           blkmem,             /**< block memory */
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
   SCIP_Real             factor              /**< constant factor */
   )
{
   int i;

   assert(polynomialdata != NULL);

   if( factor == 1.0 )
      return;

   if( factor == 0.0 )
   {
      for( i = 0; i < polynomialdata->nmonomials; ++i )
         SCIPexprFreeMonomial(blkmem, &polynomialdata->monomials[i]);
      polynomialdata->nmonomials = 0;
   }
   else
   {
      for( i = 0; i < polynomialdata->nmonomials; ++i )
         SCIPexprChgMonomialCoef(polynomialdata->monomials[i], polynomialdata->monomials[i]->coef * factor);
   }

   polynomialdata->constant *= factor;
}

/** multiplies each summand of a polynomial by a given monomial */
static
SCIP_RETCODE polynomialdataMultiplyByMonomial(
   BMS_BLKMEM*           blkmem,             /**< block memory */
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
   SCIP_EXPRDATA_MONOMIAL* factor,           /**< monomial factor */
   int*                  childmap            /**< map children in factor to children in expr, or NULL for 1:1 */
   )
{
   int i;

   assert(blkmem != NULL);
   assert(factor != NULL);
   assert(polynomialdata != NULL);

   if( factor->nfactors == 0 )
   {
      polynomialdataMultiplyByConstant(blkmem, polynomialdata, factor->coef);
      return SCIP_OKAY;
   }

   /* multiply each monomial by factor */
   for( i = 0; i < polynomialdata->nmonomials; ++i )
   {
      SCIP_CALL( SCIPexprMultiplyMonomialByMonomial(blkmem, polynomialdata->monomials[i], factor, childmap) );
   }

   /* add new monomial for constant multiplied by factor */
   if( polynomialdata->constant != 0.0 )
   {
      SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, polynomialdata->nmonomials+1) );
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &polynomialdata->monomials[polynomialdata->nmonomials], polynomialdata->constant, 0, NULL, NULL) );
      SCIP_CALL( SCIPexprMultiplyMonomialByMonomial(blkmem, polynomialdata->monomials[polynomialdata->nmonomials], factor, childmap) );
      ++polynomialdata->nmonomials;
      polynomialdata->sorted = FALSE;
      polynomialdata->constant = 0.0;
   }

   return SCIP_OKAY;
}

/** multiplies a polynomial by a polynomial
 *
 *  Factors need to be different.
 */
static
SCIP_RETCODE polynomialdataMultiplyByPolynomial(
   BMS_BLKMEM*           blkmem,             /**< block memory */
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
   SCIP_EXPRDATA_POLYNOMIAL* factordata,     /**< polynomial factor data */
   int*                  childmap            /**< map children in factor to children in polynomialdata, or NULL for 1:1 */
   )
{
   int i1;
   int i2;
   int orignmonomials;

   assert(blkmem != NULL);
   assert(polynomialdata != NULL);
   assert(factordata != NULL);
   assert(polynomialdata != factordata);

   if( factordata->nmonomials == 0 )
   {
      polynomialdataMultiplyByConstant(blkmem, polynomialdata, factordata->constant);
      return SCIP_OKAY;
   }

   if( factordata->nmonomials == 1 && factordata->constant == 0.0 )
   {
      SCIP_CALL( polynomialdataMultiplyByMonomial(blkmem, polynomialdata, factordata->monomials[0], childmap) );
      return SCIP_OKAY;
   }

   /* turn constant into a monomial, so we can assume below that constant is 0.0 */
   if( polynomialdata->constant != 0.0 )
   {
      SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, polynomialdata->nmonomials+1) );
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &polynomialdata->monomials[polynomialdata->nmonomials], polynomialdata->constant, 0, NULL, NULL) );
      ++polynomialdata->nmonomials;
      polynomialdata->sorted = FALSE;
      polynomialdata->constant = 0.0;
   }

   SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, polynomialdata->nmonomials * (factordata->nmonomials + (factordata->constant == 0.0 ? 0 : 1))) );

   /* for each monomial in factordata (except the last, if factordata->constant is 0),
    * duplicate monomials from polynomialdata and multiply them by the monomial for factordata */
   orignmonomials = polynomialdata->nmonomials;
   for( i2 = 0; i2 < factordata->nmonomials; ++i2 )
   {
      /* add a copy of original monomials to end of polynomialdata's monomials array */
      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 */
      SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, orignmonomials, polynomialdata->monomials, TRUE) );
      assert(polynomialdata->nmonomials == (i2+2) * orignmonomials);

      /* multiply each copied monomial by current monomial from factordata */
      for( i1 = (i2+1) * orignmonomials; i1 < (i2+2) * orignmonomials; ++i1 )
      {
         SCIP_CALL( SCIPexprMultiplyMonomialByMonomial(blkmem, polynomialdata->monomials[i1], factordata->monomials[i2], childmap) );
      }

      if( factordata->constant == 0.0 && i2 == factordata->nmonomials - 2 )
      {
         ++i2;
         break;
      }
   }

   if( factordata->constant != 0.0 )
   {
      assert(i2 == factordata->nmonomials);
      /* multiply original monomials in polynomialdata by constant in factordata */
      for( i1 = 0; i1 < orignmonomials; ++i1 )
         SCIPexprChgMonomialCoef(polynomialdata->monomials[i1], polynomialdata->monomials[i1]->coef * factordata->constant);
   }
   else
   {
      assert(i2 == factordata->nmonomials - 1);
      /* multiply original monomials in polynomialdata by last monomial in factordata */
      for( i1 = 0; i1 < orignmonomials; ++i1 )
      {
         SCIP_CALL( SCIPexprMultiplyMonomialByMonomial(blkmem, polynomialdata->monomials[i1], factordata->monomials[i2], childmap) );
      }
   }

   return SCIP_OKAY;
}

/** takes a power of a polynomial
 *
 *  Exponent needs to be an integer,
 *  polynomial needs to be a monomial, if exponent is negative.
 */
static
SCIP_RETCODE polynomialdataPower(
   BMS_BLKMEM*           blkmem,             /**< block memory */
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
   int                   exponent            /**< exponent of power operation */
   )
{
   SCIP_EXPRDATA_POLYNOMIAL* factor;
   int i;

   assert(blkmem != NULL);
   assert(polynomialdata != NULL);

   if( exponent == 0 )
   {
      /* x^0 = 1, except if x = 0 */
      if( polynomialdata->nmonomials == 0 && polynomialdata->constant == 0.0 )
      {
         polynomialdata->constant = 0.0;
      }
      else
      {
         polynomialdata->constant = 1.0;

         for( i = 0; i < polynomialdata->nmonomials; ++i )
            SCIPexprFreeMonomial(blkmem, &polynomialdata->monomials[i]);
         polynomialdata->nmonomials = 0;
      }

      return SCIP_OKAY;
   }

   if( exponent == 1 )
      return SCIP_OKAY;

   if( polynomialdata->nmonomials == 1 && polynomialdata->constant == 0.0 )
   {
      /* polynomial is a single monomial */
      SCIPexprMonomialPower(polynomialdata->monomials[0], exponent);
      return SCIP_OKAY;
   }

   if( polynomialdata->nmonomials == 0 )
   {
      /* polynomial is a constant */
      polynomialdata->constant = pow(polynomialdata->constant, (SCIP_Real)exponent);
      return SCIP_OKAY;
   }

   assert(exponent >= 2); /* negative exponents not allowed if more than one monom */

   /* todo improve, look how SCIPintervalPowerScalar in intervalarith.c does it */

   /* get copy of our polynomial */
   SCIP_CALL( polynomialdataCopy(blkmem, &factor, polynomialdata) );

   /* do repeated multiplication */
   for( i = 2; i <= exponent; ++i )
   {
      SCIP_CALL( polynomialdataMultiplyByPolynomial(blkmem, polynomialdata, factor, NULL) );
      polynomialdataMergeMonomials(blkmem, polynomialdata, 0.0, TRUE);
   }

   /* free copy again */
   polynomialdataFree(blkmem, &factor);

   return SCIP_OKAY;
}

/** applies a mapping of child indices to the indices used in polynomial monomials */
static
void polynomialdataApplyChildmap(
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data */
   int*                  childmap            /**< mapping of child indices */
   )
{
   SCIP_EXPRDATA_MONOMIAL* monomial;
   int i;
   int j;

   assert(polynomialdata != NULL);

   for( i = 0; i < polynomialdata->nmonomials; ++i )
   {
      monomial = polynomialdata->monomials[i];
      assert(monomial != NULL);

      for( j = 0; j < monomial->nfactors; ++j )
      {
         monomial->childidxs[j] = childmap[monomial->childidxs[j]];
         assert(monomial->childidxs[j] >= 0);
      }
      monomial->sorted = FALSE;
   }

   polynomialdata->sorted = FALSE;
}

/** replaces a factor in a monomial by a polynomial and expands the result */
static
SCIP_RETCODE polynomialdataExpandMonomialFactor(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_MESSAGEHDLR*     messagehdlr,        /**< message handler */
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata, /**< polynomial data where to expand a monomial */
   int                   monomialpos,        /**< position of monomial which factor to expand */
   int                   factorpos,          /**< position of factor in monomial to expand */
   SCIP_EXPRDATA_POLYNOMIAL* factorpolynomial,/**< polynomial that should replace factor */
   int*                  childmap,           /**< map of child indices in factorpolynomial to children of polynomial */
   int                   maxexpansionexponent,/**< maximal exponent for which polynomials (with > 1 summands) are expanded */
   SCIP_Bool*            success             /**< buffer to store whether expansion has been done */
   )
{
   SCIP_EXPRDATA_POLYNOMIAL* factorpolynomialcopy;
   SCIP_EXPRDATA_MONOMIAL* monomial;
   int i;

   assert(blkmem != NULL);
   assert(polynomialdata != NULL);
   assert(factorpolynomial != NULL);
   assert(childmap != NULL || factorpolynomial->nmonomials == 0);
   assert(success != NULL);
   assert(monomialpos >= 0);
   assert(monomialpos < polynomialdata->nmonomials);
   assert(factorpos >= 0);

   monomial = polynomialdata->monomials[monomialpos];
   assert(monomial != NULL);
   assert(factorpos < monomial->nfactors);

   *success = TRUE;

   if( factorpolynomial->nmonomials == 0 )
   {
      /* factorpolynomial is a constant */

      if( !EPSISINT(monomial->exponents[factorpos], 0.0) && factorpolynomial->constant < 0.0 )  /*lint !e835*/
      {
         /* if polynomial is a negative constant and our exponent is not integer, then cannot do expansion */
         SCIPmessagePrintWarning(messagehdlr, "got negative constant %g to the power of a noninteger exponent %g\n", factorpolynomial->constant, monomial->exponents[factorpos]);
         *success = FALSE;
         return SCIP_OKAY;
      }
      monomial->coef *= pow(factorpolynomial->constant, monomial->exponents[factorpos]);

      /* move last factor to position factorpos */
      if( factorpos < monomial->nfactors-1 )
      {
         monomial->exponents[factorpos] = monomial->exponents[monomial->nfactors-1];
         monomial->childidxs[factorpos] = monomial->childidxs[monomial->nfactors-1];
      }
      --monomial->nfactors;
      monomial->sorted = FALSE;
      polynomialdata->sorted = FALSE;

      return SCIP_OKAY;
   }

   if( factorpolynomial->constant == 0.0 && factorpolynomial->nmonomials == 1 )
   {
      /* factorpolynomial is a single monomial */
      SCIP_EXPRDATA_MONOMIAL* factormonomial;
      int childidx;
      SCIP_Real exponent;

      factormonomial = factorpolynomial->monomials[0];
      assert(factormonomial != NULL);

      if( !EPSISINT(monomial->exponents[factorpos], 0.0) )  /*lint !e835*/
      {
         if( factormonomial->coef < 0.0 )
         {
            /* if coefficient of monomial is negative and our exponent is not integer, then do not do expansion
             * @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
             */
            *success = FALSE;
            return SCIP_OKAY;
         }
         if( factormonomial->nfactors > 1 )
         {
            /* @todo if there is an even number of factors in factormonomial that are negative, then they always multiply to something positive
             * however, we cannot expand them as below, since we cannot compute the single powers
             * since we do not have the bounds on the factors here, we skip expansion in this case
             * MINLPLib instances tls2,4,6 are examples where we are loosing here (do not recognize convexity)
             */
            *success = FALSE;
            return SCIP_OKAY;
         }
      }

      SCIP_CALL( monomialdataEnsureFactorsSize(blkmem, monomial, monomial->nfactors + factormonomial->nfactors) );

      for( i = 0; i < factormonomial->nfactors; ++i )
      {
         childidx = childmap[factormonomial->childidxs[i]];  /*lint !e613*/
         /* can do this because monomial->exponents[factorpos] is assumed to be integer or factormonomial has positive coefficient and only one factor
          * thus, if factormonomial->exponents[i] is fractional, then we can assume that it's argument is positive
          */
         exponent = factormonomial->exponents[i] * monomial->exponents[factorpos];
         SCIP_CALL( SCIPexprAddMonomialFactors(blkmem, monomial, 1, &childidx, &exponent) );
      }

      monomial->coef *= pow(factormonomial->coef, monomial->exponents[factorpos]);

      /* move last factor to position factorpos */
      if( factorpos < monomial->nfactors-1 )
      {
         monomial->exponents[factorpos] = monomial->exponents[monomial->nfactors-1];
         monomial->childidxs[factorpos] = monomial->childidxs[monomial->nfactors-1];
      }
      --monomial->nfactors;
      monomial->sorted = FALSE;
      polynomialdata->sorted = FALSE;

      return SCIP_OKAY;
   }

   /* if exponent is negative or fractional and the polynomial is not just a monomial, then we cannot do expansion */
   if( !EPSISINT(monomial->exponents[factorpos], 0.0) || monomial->exponents[factorpos] < 0.0 )  /*lint !e835*/
   {
      *success = FALSE;
      return SCIP_OKAY;
   }

   /* if exponent is too large, skip expansion */
   if( monomial->exponents[factorpos] > maxexpansionexponent )
   {
      *success = FALSE;
      return SCIP_OKAY;
   }

   /* check whether maximal degree of expansion would exceed maxexpansionexponent
    * that is, assume monomial is f1^a1 f2^a2 ... and we want to expand f1 = (g11^beta11 g12^beta12... + g21^beta21 g22^beta22 ... + ...)
    * then we do this only if all ai and all beta are > 0.0 and a1 max(beta11+beta12+..., beta21+beta22+..., ...) + a2 + ... < maxexpansionexponent
    * exception (there need to be one) is if monomial is just f1
    */
   if( maxexpansionexponent < INT_MAX && (monomial->nfactors > 1 || monomial->exponents[factorpos] != 1.0) )
   {
      SCIP_Real restdegree;
      SCIP_Real degree;
      int j;

      restdegree = -monomial->exponents[factorpos];
      for( i = 0; i < monomial->nfactors; ++i )
      {
         if( monomial->exponents[i] < 0.0 )
         {
            /* ai < 0.0 */
            SCIPdebugMessage("skip expansion because factor %d in monomial has negative exponent\n", i);
            *success = FALSE;
            return SCIP_OKAY;
         }
         restdegree += monomial->exponents[i];
      }

      for( i = 0; i < factorpolynomial->nmonomials; ++i )
      {
         degree = 0.0;
         for( j = 0; j < factorpolynomial->monomials[i]->nfactors; ++j )
         {
            if( factorpolynomial->monomials[i]->exponents[j] < 0.0 )
            {
               /* beta_ij < 0.0 */
               SCIPdebugMessage("skip expansion because %d'th factor in %d'th monomial of factorpolynomial is negative\n", i, j);
               *success = FALSE;
               return SCIP_OKAY;
            }
            degree += factorpolynomial->monomials[i]->exponents[j];
         }
         if( degree * monomial->exponents[factorpos] + restdegree > maxexpansionexponent )
         {
            /* (beta_i1+beta_i2+...)*monomial->exponents[factorpos] + rest > maxexpansion */
            SCIPdebugMessage("skip expansion because degree of %d'th monomial would yield degree %g > max = %d in expansion\n",
               i, degree * monomial->exponents[factorpos] + restdegree, maxexpansionexponent);
            *success = FALSE;
            return SCIP_OKAY;
         }
      }
   }

   /* create a copy of factor */
   SCIP_CALL( polynomialdataCopy(blkmem, &factorpolynomialcopy, factorpolynomial) );
   /* apply childmap to copy */
   polynomialdataApplyChildmap(factorpolynomialcopy, childmap);
   /* create power of factor */
   SCIP_CALL( polynomialdataPower(blkmem, factorpolynomialcopy, (int)EPSFLOOR(monomial->exponents[factorpos], 0.0)) );  /*lint !e835*/

   /* remove factor from monomial by moving last factor to position factorpos */
   if( factorpos < monomial->nfactors-1 )
   {
      monomial->exponents[factorpos] = monomial->exponents[monomial->nfactors-1];
      monomial->childidxs[factorpos] = monomial->childidxs[monomial->nfactors-1];
   }
   --monomial->nfactors;
   monomial->sorted = FALSE;

   /* multiply factor with this reduced monomial */
   SCIP_CALL( polynomialdataMultiplyByMonomial(blkmem, factorpolynomialcopy, monomial, NULL) );

   /* remove monomial from polynomial and move last monomial to monomialpos */
   SCIPexprFreeMonomial(blkmem, &polynomialdata->monomials[monomialpos]);
   if( monomialpos < polynomialdata->nmonomials-1 )
      polynomialdata->monomials[monomialpos] = polynomialdata->monomials[polynomialdata->nmonomials-1];
   --polynomialdata->nmonomials;
   polynomialdata->sorted = FALSE;

   /* add factorpolynomialcopy to polynomial */
   SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, factorpolynomialcopy->nmonomials, factorpolynomialcopy->monomials, FALSE) );
   polynomialdata->constant += factorpolynomialcopy->constant;

   factorpolynomialcopy->nmonomials = 0;
   polynomialdataFree(blkmem, &factorpolynomialcopy);

   return SCIP_OKAY;
}

/**@} */

/**@name Expression operand private methods */
/**@{ */

/** a default implementation of expression interval evaluation that always gives a correct result */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntDefault )
{   /*lint --e{715}*/
   SCIPintervalSetEntire(infinity, result);

   return SCIP_OKAY;
}

/** a default implementation of expression curvature check that always gives a correct result */
static
SCIP_DECL_EXPRCURV( exprcurvDefault )
{   /*lint --e{715}*/
   *result = SCIP_EXPRCURV_UNKNOWN;

   return SCIP_OKAY;
}

/** point evaluation for EXPR_VAR */
static
SCIP_DECL_EXPREVAL( exprevalVar )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(varvals != NULL);

   *result = varvals[opdata.intval];

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_VAR */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntVar )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(varvals != NULL);

   *result = varvals[opdata.intval];

   return SCIP_OKAY;
}

/** curvature for EXPR_VAR */
static
SCIP_DECL_EXPRCURV( exprcurvVar )
{   /*lint --e{715}*/
   assert(result  != NULL);

   *result = SCIP_EXPRCURV_LINEAR;

   return SCIP_OKAY;
}

/** point evaluation for EXPR_CONST */
static
SCIP_DECL_EXPREVAL( exprevalConst )
{   /*lint --e{715}*/
   assert(result != NULL);

   *result = opdata.dbl;

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_CONST */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntConst )
{   /*lint --e{715}*/
   assert(result != NULL);

   SCIPintervalSet(result, opdata.dbl);

   return SCIP_OKAY;
}

/** curvature for EXPR_CONST */
static
SCIP_DECL_EXPRCURV( exprcurvConst )
{   /*lint --e{715}*/
   assert(result  != NULL);

   *result = SCIP_EXPRCURV_LINEAR;

   return SCIP_OKAY;
}

/** point evaluation for EXPR_PARAM */
static
SCIP_DECL_EXPREVAL( exprevalParam )
{   /*lint --e{715}*/
   assert(result    != NULL);
   assert(paramvals != NULL );

   *result = paramvals[opdata.intval];

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_PARAM */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntParam )
{   /*lint --e{715}*/
   assert(result    != NULL);
   assert(paramvals != NULL );

   SCIPintervalSet(result, paramvals[opdata.intval]);

   return SCIP_OKAY;
}

/** curvature for EXPR_PARAM */
static
SCIP_DECL_EXPRCURV( exprcurvParam )
{   /*lint --e{715}*/
   assert(result  != NULL);

   *result = SCIP_EXPRCURV_LINEAR;

   return SCIP_OKAY;
}

/** point evaluation for EXPR_PLUS */
static
SCIP_DECL_EXPREVAL( exprevalPlus )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = argvals[0] + argvals[1];

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_PLUS */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntPlus )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalAdd(infinity, result, argvals[0], argvals[1]);

   return SCIP_OKAY;
}

/** curvature for EXPR_PLUS */
static
SCIP_DECL_EXPRCURV( exprcurvPlus )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argcurv != NULL);

   *result = SCIPexprcurvAdd(argcurv[0], argcurv[1]);

   return SCIP_OKAY;
}

/** point evaluation for EXPR_MINUS */
static
SCIP_DECL_EXPREVAL( exprevalMinus )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = argvals[0] - argvals[1];

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_MINUS */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntMinus )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalSub(infinity, result, argvals[0], argvals[1]);

   return SCIP_OKAY;
}

/** curvature for EXPR_MINUS */
static
SCIP_DECL_EXPRCURV( exprcurvMinus )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argcurv != NULL);

   *result = SCIPexprcurvAdd(argcurv[0], SCIPexprcurvNegate(argcurv[1]));

   return SCIP_OKAY;
}

/** point evaluation for EXPR_MUL */
static
SCIP_DECL_EXPREVAL( exprevalMult )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = argvals[0] * argvals[1];

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_MUL */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntMult )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalMul(infinity, result, argvals[0], argvals[1]);

   return SCIP_OKAY;
}

/** curvature for EXPR_MUL */
static
SCIP_DECL_EXPRCURV( exprcurvMult )
{   /*lint --e{715}*/
   assert(result    != NULL);
   assert(argcurv   != NULL);
   assert(argbounds != NULL);

   /* if one factor is constant, then product is
    * - linear, if constant is 0.0
    * - same curvature as other factor, if constant is positive
    * - negated curvature of other factor, if constant is negative
    *
    * if both factors are not constant, then product may not be convex nor concave
    */
   if( argbounds[1].inf == argbounds[1].sup )  /*lint !e777*/
      *result = SCIPexprcurvMultiply(argbounds[1].inf, argcurv[0]);
   else if( argbounds[0].inf == argbounds[0].sup )  /*lint !e777*/
      *result = SCIPexprcurvMultiply(argbounds[0].inf, argcurv[1]);
   else
      *result = SCIP_EXPRCURV_UNKNOWN;

   return SCIP_OKAY;
}

/** point evaluation for EXPR_DIV */
static
#if defined(__GNUC__) && __GNUC__ * 100 + __GNUC_MINOR__ * 10 >= 490 && !defined(__INTEL_COMPILER)
__attribute__((no_sanitize_undefined))
#endif
SCIP_DECL_EXPREVAL( exprevalDiv )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = argvals[0] / argvals[1];

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_DIV */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntDiv )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalDiv(infinity, result, argvals[0], argvals[1]);

   return SCIP_OKAY;
}

/** curvature for EXPR_DIV */
static
SCIP_DECL_EXPRCURV( exprcurvDiv )
{   /*lint --e{715}*/
   assert(result    != NULL);
   assert(argcurv   != NULL);
   assert(argbounds != NULL);

   /* if denominator is constant, then quotient has curvature sign(denominator) * curv(nominator)
    *
    * if nominator is a constant, then quotient is
    * - sign(nominator) * convex, if denominator is concave and positive
    * - sign(nominator) * concave, if denominator is convex and negative
    *
    * if denominator is positive but convex, then we don't know, e.g.,
    *   - 1/x^2 is convex for x>=0
    *   - 1/(1+(x-1)^2) is neither convex nor concave for x >= 0
    *
    * if both nominator and denominator are not constant, then quotient may not be convex nor concave
    */
   if( argbounds[1].inf == argbounds[1].sup )  /*lint !e777*/
   {
      /* denominator is constant */
      *result = SCIPexprcurvMultiply(argbounds[1].inf, argcurv[0]);
   }
   else if( argbounds[0].inf == argbounds[0].sup )  /*lint !e777*/
   {
      /* nominator is constant */
      if( argbounds[1].inf >= 0.0 && (argcurv[1] & SCIP_EXPRCURV_CONCAVE) )
         *result = SCIPexprcurvMultiply(argbounds[0].inf, SCIP_EXPRCURV_CONVEX);
      else if( argbounds[1].sup <= 0.0 && (argcurv[1] & SCIP_EXPRCURV_CONVEX) )
         *result = SCIPexprcurvMultiply(argbounds[0].inf, SCIP_EXPRCURV_CONCAVE);
      else
         *result = SCIP_EXPRCURV_UNKNOWN;
   }
   else
   {
      /* denominator and nominator not constant */
      *result = SCIP_EXPRCURV_UNKNOWN;
   }

   return SCIP_OKAY;
}

/** point evaluation for EXPR_SQUARE */
static
SCIP_DECL_EXPREVAL( exprevalSquare )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = argvals[0] * argvals[0];

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_SQUARE */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntSquare )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalSquare(infinity, result, argvals[0]);

   return SCIP_OKAY;
}

/** curvature for EXPR_SQUARE */
static
SCIP_DECL_EXPRCURV( exprcurvSquare )
{   /*lint --e{715}*/
   assert(result    != NULL);
   assert(argcurv   != NULL);
   assert(argbounds != NULL);

   *result = SCIPexprcurvPower(argbounds[0], argcurv[0], 2.0);

   return SCIP_OKAY;
}

/** point evaluation for EXPR_SQRT */
static
SCIP_DECL_EXPREVAL( exprevalSquareRoot )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = sqrt(argvals[0]);

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_SQRT */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntSquareRoot )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalSquareRoot(infinity, result, argvals[0]);

   return SCIP_OKAY;
}

/** curvature for EXPR_SQRT */
static
SCIP_DECL_EXPRCURV( exprcurvSquareRoot )
{   /*lint --e{715}*/
   assert(result    != NULL);
   assert(argcurv   != NULL);

   /* square-root is concave, if child is concave
    * otherwise, we don't know
    */

   if( argcurv[0] & SCIP_EXPRCURV_CONCAVE )
      *result = SCIP_EXPRCURV_CONCAVE;
   else
      *result = SCIP_EXPRCURV_UNKNOWN;

   return SCIP_OKAY;
}

/** point evaluation for EXPR_REALPOWER */
static
SCIP_DECL_EXPREVAL( exprevalRealPower )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = pow(argvals[0], opdata.dbl);

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_REALPOWER */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntRealPower )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalPowerScalar(infinity, result, argvals[0], opdata.dbl);

   return SCIP_OKAY;
}

/** curvature for EXPR_REALPOWER */
static
SCIP_DECL_EXPRCURV( exprcurvRealPower )
{   /*lint --e{715}*/
   assert(result    != NULL);
   assert(argcurv   != NULL);
   assert(argbounds != NULL);

   *result = SCIPexprcurvPower(argbounds[0], argcurv[0], opdata.dbl);

   return SCIP_OKAY;
}

/** point evaluation for EXPR_INTPOWER */
static
#if defined(__GNUC__) && __GNUC__ * 100 + __GNUC_MINOR__ * 10 >= 490 && !defined(__INTEL_COMPILER)
__attribute__((no_sanitize_undefined))
#endif
SCIP_DECL_EXPREVAL( exprevalIntPower )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   switch( opdata.intval )
   {
   case -1:
      *result = 1.0 / argvals[0];
      return SCIP_OKAY;

   case 0:
      *result = 1.0;
      return SCIP_OKAY;

   case 1:
      *result = argvals[0];
      return SCIP_OKAY;

   case 2:
      *result = argvals[0] * argvals[0];
      return SCIP_OKAY;

   default:
      *result = pow(argvals[0], (SCIP_Real)opdata.intval);
   }

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_INTPOWER */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntIntPower )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalPowerScalar(infinity, result, argvals[0], (SCIP_Real)opdata.intval);

   return SCIP_OKAY;
}

/** curvature for EXPR_INTPOWER */
static
SCIP_DECL_EXPRCURV( exprcurvIntPower )
{   /*lint --e{715}*/
   assert(result    != NULL);
   assert(argcurv   != NULL);
   assert(argbounds != NULL);

   *result = SCIPexprcurvPower(argbounds[0], argcurv[0], (SCIP_Real)opdata.intval);

   return SCIP_OKAY;
}

/** point evaluation for EXPR_SIGNPOWER */
static
SCIP_DECL_EXPREVAL( exprevalSignPower )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   if( argvals[0] > 0 )
      *result =  pow( argvals[0], opdata.dbl);
   else
      *result = -pow(-argvals[0], opdata.dbl);

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_SIGNPOWER */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntSignPower )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalSignPowerScalar(infinity, result, argvals[0], opdata.dbl);

   return SCIP_OKAY;
}

/** curvature for EXPR_SIGNPOWER */
static
SCIP_DECL_EXPRCURV( exprcurvSignPower )
{   /*lint --e{715}*/
   SCIP_INTERVAL tmp;
   SCIP_EXPRCURV left;
   SCIP_EXPRCURV right;

   assert(result    != NULL);
   assert(argcurv   != NULL);
   assert(argbounds != NULL);

   /* for x <= 0, signpower(x,c) = -(-x)^c
    * for x >= 0, signpower(x,c) =  ( x)^c
    *
    * thus, get curvatures for both parts and "intersect" them
    */

   if( argbounds[0].inf < 0 )
   {
      SCIPintervalSetBounds(&tmp, 0.0, -argbounds[0].inf);
      left = SCIPexprcurvNegate(SCIPexprcurvPower(tmp, SCIPexprcurvNegate(argcurv[0]), opdata.dbl));
   }
   else
   {
      left = SCIP_EXPRCURV_LINEAR;
   }

   if( argbounds[0].sup > 0 )
   {
      SCIPintervalSetBounds(&tmp, 0.0,  argbounds[0].sup);
      right = SCIPexprcurvPower(tmp, argcurv[0], opdata.dbl);
   }
   else
   {
      right = SCIP_EXPRCURV_LINEAR;
   }

   *result = (SCIP_EXPRCURV) (left & right);

   return SCIP_OKAY;
}

/** point evaluation for EXPR_EXP */
static
SCIP_DECL_EXPREVAL( exprevalExp )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = exp(argvals[0]);

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_EXP */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntExp )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalExp(infinity, result, argvals[0]);

   return SCIP_OKAY;
}

/** curvature for EXPR_EXP */
static
SCIP_DECL_EXPRCURV( exprcurvExp )
{   /*lint --e{715}*/
   assert(result    != NULL);
   assert(argcurv   != NULL);

   /* expression is convex if child is convex
    * otherwise, we don't know
    */
   if( argcurv[0] & SCIP_EXPRCURV_CONVEX )
      *result = SCIP_EXPRCURV_CONVEX;
   else
      *result = SCIP_EXPRCURV_UNKNOWN;

   return SCIP_OKAY;
}

/** point evaluation for EXPR_LOG */
static
SCIP_DECL_EXPREVAL( exprevalLog )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = log(argvals[0]);

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_LOG */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntLog )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalLog(infinity, result, argvals[0]);

   return SCIP_OKAY;
}

/** curvature for EXPR_LOG */
static
SCIP_DECL_EXPRCURV( exprcurvLog )
{   /*lint --e{715}*/
   assert(result    != NULL);
   assert(argcurv   != NULL);

   /* expression is concave if child is concave
    * otherwise, we don't know
    */
   if( argcurv[0] & SCIP_EXPRCURV_CONCAVE )
      *result = SCIP_EXPRCURV_CONCAVE;
   else
      *result = SCIP_EXPRCURV_UNKNOWN;

   return SCIP_OKAY;
}

/** point evaluation for EXPR_SIN */
static
SCIP_DECL_EXPREVAL( exprevalSin )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = sin(argvals[0]);

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_SIN */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntSin )
{   /*lint --e{715}*/
   assert(result != NULL);
   assert(argvals != NULL);
   assert(nargs == 1);

   SCIPintervalSin(infinity, result, *argvals);

   return SCIP_OKAY;
}

/* @todo implement exprcurvSin */
#define exprcurvSin exprcurvDefault

/** point evaluation for EXPR_COS */
static
SCIP_DECL_EXPREVAL( exprevalCos )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = cos(argvals[0]);

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_COS */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntCos )
{   /*lint --e{715}*/
   assert(result != NULL);
   assert(argvals != NULL);
   assert(nargs == 1);

   SCIPintervalCos(infinity, result, *argvals);

   return SCIP_OKAY;
}

/* @todo implement exprcurvCos */
#define exprcurvCos exprcurvDefault

/** point evaluation for EXPR_TAN */
static
SCIP_DECL_EXPREVAL( exprevalTan )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = tan(argvals[0]);

   return SCIP_OKAY;
}

/* @todo implement SCIPintervalTan */
#define exprevalIntTan exprevalIntDefault

/* @todo implement exprcurvTan */
#define exprcurvTan exprcurvDefault

/* erf and erfi do not seem to exist on every system, and we cannot really handle them anyway, so they are currently disabled */
#ifdef SCIP_DISABLED_CODE
static
SCIP_DECL_EXPREVAL( exprevalErf )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = erf(argvals[0]);

   return SCIP_OKAY;
}

/* @todo implement SCIPintervalErf */
#define exprevalIntErf exprevalIntDefault

/* @todo implement SCIPintervalErf */
#define exprcurvErf exprcurvDefault

static
SCIP_DECL_EXPREVAL( exprevalErfi )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   /* @TODO implement erfi evaluation */
   SCIPerrorMessage("erfi not implemented");

   return SCIP_ERROR;
}

/* @todo implement SCIPintervalErfi */
#define exprevalIntErfi NULL

#define exprcurvErfi exprcurvDefault
#endif

/** point evaluation for EXPR_MIN */
static
SCIP_DECL_EXPREVAL( exprevalMin )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = MIN(argvals[0], argvals[1]);

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_MIN */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntMin )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalMin(infinity, result, argvals[0], argvals[1]);

   return SCIP_OKAY;
}

/** curvature for EXPR_MIN */
static
SCIP_DECL_EXPRCURV( exprcurvMin )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argcurv != NULL);

   /* the minimum of two concave functions is concave
    * otherwise, we don't know
    */

   if( (argcurv[0] & SCIP_EXPRCURV_CONCAVE) && (argcurv[1] & SCIP_EXPRCURV_CONCAVE) )
      *result = SCIP_EXPRCURV_CONCAVE;
   else
      *result = SCIP_EXPRCURV_UNKNOWN;

   return SCIP_OKAY;
}

/** point evaluation for EXPR_MAX */
static
SCIP_DECL_EXPREVAL( exprevalMax )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = MAX(argvals[0], argvals[1]);

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_MAX */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntMax )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalMax(infinity, result, argvals[0], argvals[1]);

   return SCIP_OKAY;
}

/** curvature for EXPR_MAX */
static
SCIP_DECL_EXPRCURV( exprcurvMax )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argcurv != NULL);

   /* the maximum of two convex functions is convex
    * otherwise, we don't know
    */
   if( (argcurv[0] & SCIP_EXPRCURV_CONVEX) && (argcurv[1] & SCIP_EXPRCURV_CONVEX) )
      *result = SCIP_EXPRCURV_CONVEX;
   else
      *result = SCIP_EXPRCURV_UNKNOWN;

   return SCIP_OKAY;
}

/** point evaluation for EXPR_ABS */
static
SCIP_DECL_EXPREVAL( exprevalAbs )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = ABS(argvals[0]);

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_ABS */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntAbs )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalAbs(infinity, result, argvals[0]);

   return SCIP_OKAY;
}

/** curvature for EXPR_ABS */
static
SCIP_DECL_EXPRCURV( exprcurvAbs )
{   /*lint --e{715}*/
   assert(result    != NULL);
   assert(argcurv   != NULL);
   assert(argbounds != NULL);

   /* if child is only negative, then abs(child) = -child
    * if child is only positive, then abs(child) = child
    * if child is both positive and negative, but also linear, then abs(child) is convex
    * otherwise, we don't know
    */
   if( argbounds[0].sup <= 0.0 )
      *result = SCIPexprcurvMultiply(-1.0, argcurv[0]);
   else if( argbounds[0].inf >= 0.0 )
      *result = argcurv[0];
   else if( argcurv[0] == SCIP_EXPRCURV_LINEAR )
      *result = SCIP_EXPRCURV_CONVEX;
   else
      *result = SCIP_EXPRCURV_UNKNOWN;

   return SCIP_OKAY;
}

/** point evaluation for EXPR_SIGN */
static
SCIP_DECL_EXPREVAL( exprevalSign )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   *result = SIGN(argvals[0]);

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_SIGN */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntSign )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalSign(infinity, result, argvals[0]);

   return SCIP_OKAY;
}

/** curvature for EXPR_SIGN */
static
SCIP_DECL_EXPRCURV( exprcurvSign )
{   /*lint --e{715}*/
   assert(result    != NULL);
   assert(argbounds != NULL);

   /* if sign of child is clear, then sign is linear otherwise, we don't know */
   if( argbounds[0].sup <= 0.0 || argbounds[0].inf >= 0.0 )
      *result = SCIP_EXPRCURV_LINEAR;
   else
      *result = SCIP_EXPRCURV_UNKNOWN;

   return SCIP_OKAY;
}

/** point evaluation for EXPR_SUM */
static
SCIP_DECL_EXPREVAL( exprevalSum )
{   /*lint --e{715}*/
   int i;

   assert(result  != NULL);
   assert(argvals != NULL);

   *result = 0.0;
   for( i = 0; i < nargs; ++i )
      *result += argvals[i];

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_SUM */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntSum )
{   /*lint --e{715}*/
   int i;

   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalSet(result, 0.0);

   for( i = 0; i < nargs; ++i )
      SCIPintervalAdd(infinity, result, *result, argvals[i]);

   return SCIP_OKAY;
}

/** curvature for EXPR_SUM */
static
SCIP_DECL_EXPRCURV( exprcurvSum )
{   /*lint --e{715}*/
   int i;

   assert(result  != NULL);
   assert(argcurv != NULL);

   *result = SCIP_EXPRCURV_LINEAR;

   for( i = 0; i < nargs; ++i )
      *result = SCIPexprcurvAdd(*result, argcurv[i]);

   return SCIP_OKAY;
}

/** point evaluation for EXPR_PRODUCT */
static
SCIP_DECL_EXPREVAL( exprevalProduct )
{   /*lint --e{715}*/
   int i;

   assert(result  != NULL);
   assert(argvals != NULL);

   *result = 1.0;
   for( i = 0; i < nargs; ++i )
      *result *= argvals[i];

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_PRODUCT */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntProduct )
{   /*lint --e{715}*/
   int i;

   assert(result  != NULL);
   assert(argvals != NULL);

   SCIPintervalSet(result, 1.0);

   for( i = 0; i < nargs; ++i )
      SCIPintervalMul(infinity, result, *result, argvals[i]);

   return SCIP_OKAY;
}

/** curvature for EXPR_PRODUCT */
static
SCIP_DECL_EXPRCURV( exprcurvProduct )
{   /*lint --e{715}*/
   SCIP_Bool hadnonconst;
   SCIP_Real constants;
   int i;

   assert(result    != NULL);
   assert(argcurv   != NULL);
   assert(argbounds != NULL);

   /* if all factors are constant, then product is linear (even constant)
    * if only one factor is not constant, then product is curvature of this factor, multiplied by sign of product of remaining factors
    */
   *result = SCIP_EXPRCURV_LINEAR;
   hadnonconst = FALSE;
   constants = 1.0;

   for( i = 0; i < nargs; ++i )
   {
      if( argbounds[i].inf == argbounds[i].sup )  /*lint !e777*/
      {
         constants *= argbounds[i].inf;
      }
      else if( !hadnonconst )
      {
         /* first non-constant child */
         *result = argcurv[i];
         hadnonconst = TRUE;
      }
      else
      {
         /* more than one non-constant child, thus don't know curvature */
         *result = SCIP_EXPRCURV_UNKNOWN;
         break;
      }
   }

   *result = SCIPexprcurvMultiply(constants, *result);

   return SCIP_OKAY;
}

/** point evaluation for EXPR_LINEAR */
static
SCIP_DECL_EXPREVAL( exprevalLinear )
{   /*lint --e{715}*/
   SCIP_Real* coef;
   int i;

   assert(result  != NULL);
   assert(argvals != NULL || nargs == 0);
   assert(opdata.data != NULL);

   coef = &((SCIP_Real*)opdata.data)[nargs];

   *result = *coef;
   for( i = nargs-1, --coef; i >= 0; --i, --coef )
      *result += *coef * argvals[i];  /*lint !e613*/

   assert(++coef == (SCIP_Real*)opdata.data);

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_LINEAR */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntLinear )
{   /*lint --e{715}*/
   assert(result  != NULL);
   assert(argvals != NULL || nargs == 0);
   assert(opdata.data != NULL);

   SCIPintervalScalprodScalars(infinity, result, nargs, argvals, (SCIP_Real*)opdata.data);
   SCIPintervalAddScalar(infinity, result, *result, ((SCIP_Real*)opdata.data)[nargs]);

   return SCIP_OKAY;
}

/** curvature for EXPR_LINEAR */
static
SCIP_DECL_EXPRCURV( exprcurvLinear )
{   /*lint --e{715}*/
   SCIP_Real* data;
   int i;

   assert(result  != NULL);
   assert(argcurv != NULL);

   data = (SCIP_Real*)opdata.data;
   assert(data != NULL);

   *result = SCIP_EXPRCURV_LINEAR;

   for( i = 0; i < nargs; ++i )
      *result = SCIPexprcurvAdd(*result, SCIPexprcurvMultiply(data[i], argcurv[i]));

   return SCIP_OKAY;
}

/** expression data copy for EXPR_LINEAR */
static
SCIP_DECL_EXPRCOPYDATA( exprCopyDataLinear )
{  /*lint --e{715}*/
   SCIP_Real* targetdata;

   assert(blkmem != NULL);
   assert(nchildren >= 0);
   assert(opdatatarget != NULL);

   /* for a linear expression, we need to copy the array that holds the coefficients and constant term */
   assert(opdatasource.data != NULL);
   SCIP_ALLOC( BMSduplicateBlockMemoryArray(blkmem, &targetdata, (SCIP_Real*)opdatasource.data, nchildren + 1) );  /*lint !e866*/
   opdatatarget->data = targetdata;

   return SCIP_OKAY;
}

/** expression data free for EXPR_LINEAR */
static
SCIP_DECL_EXPRFREEDATA( exprFreeDataLinear )
{  /*lint --e{715}*/
   SCIP_Real* freedata;

   assert(blkmem != NULL);
   assert(nchildren >= 0);

   freedata = (SCIP_Real*)opdata.data;
   assert(freedata != NULL);

   BMSfreeBlockMemoryArray(blkmem, &freedata, nchildren + 1);  /*lint !e866*/
}

/** point evaluation for EXPR_QUADRATIC */
static
SCIP_DECL_EXPREVAL( exprevalQuadratic )
{   /*lint --e{715}*/
   SCIP_EXPRDATA_QUADRATIC* quaddata;
   SCIP_Real* lincoefs;
   SCIP_QUADELEM* quadelems;
   int nquadelems;
   int i;

   assert(result  != NULL);
   assert(argvals != NULL || nargs == 0);

   quaddata = (SCIP_EXPRDATA_QUADRATIC*)opdata.data;
   assert(quaddata != NULL);

   lincoefs   = quaddata->lincoefs;
   nquadelems = quaddata->nquadelems;
   quadelems  = quaddata->quadelems;

   assert(quadelems != NULL || nquadelems == 0);
   assert(argvals != NULL || nquadelems == 0);

   *result = quaddata->constant;

   if( lincoefs != NULL )
   {
      for( i = nargs-1; i >= 0; --i )
         *result += lincoefs[i] * argvals[i];  /*lint !e613*/
   }

   for( i = 0; i < nquadelems; ++i, ++quadelems )  /*lint !e613*/
      *result += quadelems->coef * argvals[quadelems->idx1] * argvals[quadelems->idx2];  /*lint !e613*/

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_QUADRATIC */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntQuadratic )
{   /*lint --e{715}*/
   SCIP_EXPRDATA_QUADRATIC* quaddata;
   SCIP_Real* lincoefs;
   SCIP_QUADELEM* quadelems;
   int nquadelems;
   int i;
   int argidx;
   SCIP_Real sqrcoef;
   SCIP_INTERVAL lincoef;
   SCIP_INTERVAL tmp;

   assert(result  != NULL);
   assert(argvals != NULL || nargs == 0);

   quaddata = (SCIP_EXPRDATA_QUADRATIC*)opdata.data;
   assert(quaddata != NULL);

   lincoefs   = quaddata->lincoefs;
   nquadelems = quaddata->nquadelems;
   quadelems  = quaddata->quadelems;

   assert(quadelems != NULL || nquadelems == 0);
   assert(argvals   != NULL || nargs == 0);

   /* something fast for case of only one child */
   if( nargs == 1 )
   {
      SCIPintervalSet(&lincoef, lincoefs != NULL ? lincoefs[0] : 0.0);

      sqrcoef = 0.0;
      for( i = 0; i < nquadelems; ++i )
      {
         assert(quadelems[i].idx1 == 0);  /*lint !e613*/
         assert(quadelems[i].idx2 == 0);  /*lint !e613*/
         sqrcoef += quadelems[i].coef;    /*lint !e613*/
      }

      SCIPintervalQuad(infinity, result, sqrcoef, lincoef, argvals[0]);  /*lint !e613*/
      SCIPintervalAddScalar(infinity, result, *result, quaddata->constant);

      return SCIP_OKAY;
   }

   if( nargs == 2 && nquadelems > 0 )
   {
      /* if it's a bivariate quadratic expression with bilinear term, do something special */
      SCIP_Real ax;  /* square coefficient of first  child */
      SCIP_Real ay;  /* square coefficient of second child */
      SCIP_Real axy; /* bilinear coefficient */

      ax = 0.0;
      ay = 0.0;
      axy = 0.0;
      for( i = 0; i < nquadelems; ++i )
         if( quadelems[i].idx1 == 0 && quadelems[i].idx2 == 0 )       /*lint !e613*/
            ax += quadelems[i].coef;                                  /*lint !e613*/
         else if( quadelems[i].idx1 == 1 && quadelems[i].idx2 == 1 )  /*lint !e613*/
            ay += quadelems[i].coef;                                  /*lint !e613*/
         else
            axy += quadelems[i].coef;                                 /*lint !e613*/

      SCIPintervalQuadBivar(infinity, result, ax, ay, axy,
         lincoefs != NULL ? lincoefs[0] : 0.0, lincoefs != NULL ? lincoefs[1] : 0.0,
         argvals[0], argvals[1]);                                     /*lint !e613*/
      SCIPdebugMessage("%g x^2 + %g y^2 + %g x y + %g x + %g y = [%g,%g] for x = [%g,%g], y = [%g,%g]\n",
         ax, ay, axy, lincoefs != NULL ? lincoefs[0] : 0.0, lincoefs != NULL ? lincoefs[1] : 0.0,
         result->inf, result->sup, argvals[0].inf, argvals[0].sup, argvals[1].inf, argvals[1].sup);  /*lint !e613*/

      SCIPintervalAddScalar(infinity, result, *result, quaddata->constant);

      return SCIP_OKAY;
   }

   /* make sure coefficients are sorted */
   quadraticdataSort(quaddata);

   SCIPintervalSet(result, quaddata->constant);

   /* for each argument, we collect it's linear index from lincoefs, it's square coefficients and all factors from bilinear terms
    * then we compute the interval sqrcoef*x^2 + lincoef*x and add it to result
    * @todo split quadratic expression into bivariate quadratic terms and apply the above method
    */
   i = 0;
   for( argidx = 0; argidx < nargs; ++argidx )
   {
      if( i == nquadelems || quadelems[i].idx1 > argidx )  /*lint !e613*/
      {
         /* there are no quadratic terms with argidx in its first argument, that should be easy to handle */
         if( lincoefs != NULL )
         {
            SCIPintervalMulScalar(infinity, &tmp, argvals[argidx], lincoefs[argidx]);  /*lint !e613*/
            SCIPintervalAdd(infinity, result, *result, tmp);
         }
         continue;
      }

      sqrcoef = 0.0;
      SCIPintervalSet(&lincoef, lincoefs != NULL ? lincoefs[argidx] : 0.0);

      assert(i < nquadelems && quadelems[i].idx1 == argidx);  /*lint !e613*/
      do
      {
         if( quadelems[i].idx2 == argidx )  /*lint !e613*/
         {
            sqrcoef += quadelems[i].coef;   /*lint !e613*/
         }
         else
         {
            SCIPintervalMulScalar(infinity, &tmp, argvals[quadelems[i].idx2], quadelems[i].coef);  /*lint !e613*/
            SCIPintervalAdd(infinity, &lincoef, lincoef, tmp);
         }
         ++i;
      }
      while( i < nquadelems && quadelems[i].idx1 == argidx );  /*lint !e613*/
      assert(i == nquadelems || quadelems[i].idx1 > argidx);   /*lint !e613*/

      SCIPintervalQuad(infinity, &tmp, sqrcoef, lincoef, argvals[argidx]);  /*lint !e613*/
      SCIPintervalAdd(infinity, result, *result, tmp);
   }
   assert(i == nquadelems);

   return SCIP_OKAY;
}

/** curvature for EXPR_QUADRATIC */
static
SCIP_DECL_EXPRCURV( exprcurvQuadratic )
{   /*lint --e{715}*/
   SCIP_EXPRDATA_QUADRATIC* data;
   SCIP_QUADELEM* quadelems;
   int nquadelems;
   SCIP_Real* lincoefs;
   int i;

   assert(result    != NULL);
   assert(argcurv   != NULL);
   assert(argbounds != NULL);

   data = (SCIP_EXPRDATA_QUADRATIC*)opdata.data;
   assert(data != NULL);

   lincoefs   = data->lincoefs;
   quadelems  = data->quadelems;
   nquadelems = data->nquadelems;

   *result = SCIP_EXPRCURV_LINEAR;

   if( lincoefs != NULL )
      for( i = 0; i < nargs; ++i )
         *result = SCIPexprcurvAdd(*result, SCIPexprcurvMultiply(lincoefs[i], argcurv[i]));

   /* @todo could try cholesky factorization if all children linear...
    * @todo should then cache the result
    */
   for( i = 0; i < nquadelems && *result != SCIP_EXPRCURV_UNKNOWN; ++i )
   {
      if( quadelems[i].coef == 0.0 )
         continue;

      if( argbounds[quadelems[i].idx1].inf == argbounds[quadelems[i].idx1].sup &&  /*lint !e777*/
         +argbounds[quadelems[i].idx2].inf == argbounds[quadelems[i].idx2].sup
         )  /*lint !e777*/
      {
         /* both factors are constants -> curvature does not change */
         continue;
      }

      if( argbounds[quadelems[i].idx1].inf == argbounds[quadelems[i].idx1].sup )  /*lint !e777*/
      {
         /* first factor is constant, second is not -> add curvature of second */
         *result = SCIPexprcurvAdd(*result, SCIPexprcurvMultiply(quadelems[i].coef * argbounds[quadelems[i].idx1].inf, argcurv[quadelems[i].idx2]));
      }
      else if( argbounds[quadelems[i].idx2].inf == argbounds[quadelems[i].idx2].sup )  /*lint !e777*/
      {
         /* first factor is not constant, second is -> add curvature of first */
         *result = SCIPexprcurvAdd(*result, SCIPexprcurvMultiply(quadelems[i].coef * argbounds[quadelems[i].idx2].inf, argcurv[quadelems[i].idx1]));
      }
      else if( quadelems[i].idx1 == quadelems[i].idx2 )
      {
         /* both factors not constant, but the same (square term) */
         *result = SCIPexprcurvAdd(*result, SCIPexprcurvMultiply(quadelems[i].coef, SCIPexprcurvPower(argbounds[quadelems[i].idx1], argcurv[quadelems[i].idx1], 2.0)));
      }
      else
      {
         /* two different non-constant factors -> can't tell about curvature */
         *result = SCIP_EXPRCURV_UNKNOWN;
      }
   }

   return SCIP_OKAY;
}

/** expression data copy for EXPR_QUADRATIC */
static
SCIP_DECL_EXPRCOPYDATA( exprCopyDataQuadratic )
{  /*lint --e{715}*/
   SCIP_EXPRDATA_QUADRATIC* sourcedata;

   assert(blkmem != NULL);
   assert(opdatatarget != NULL);

   sourcedata = (SCIP_EXPRDATA_QUADRATIC*)opdatasource.data;
   assert(sourcedata != NULL);

   SCIP_CALL( quadraticdataCreate(blkmem, (SCIP_EXPRDATA_QUADRATIC**)&opdatatarget->data,
         sourcedata->constant, nchildren, sourcedata->lincoefs, sourcedata->nquadelems, sourcedata->quadelems) );

   return SCIP_OKAY;
}

/** expression data free for EXPR_QUADRATIC */
static
SCIP_DECL_EXPRFREEDATA( exprFreeDataQuadratic )
{  /*lint --e{715}*/
   SCIP_EXPRDATA_QUADRATIC* quadraticdata;

   assert(blkmem != NULL);
   assert(nchildren >= 0);

   quadraticdata = (SCIP_EXPRDATA_QUADRATIC*)opdata.data;
   assert(quadraticdata != NULL);

   if( quadraticdata->lincoefs != NULL )
   {
      BMSfreeBlockMemoryArray(blkmem, &quadraticdata->lincoefs, nchildren);
   }

   if( quadraticdata->nquadelems > 0 )
   {
      assert(quadraticdata->quadelems != NULL);
      BMSfreeBlockMemoryArray(blkmem, &quadraticdata->quadelems, quadraticdata->nquadelems);
   }

   BMSfreeBlockMemory(blkmem, &quadraticdata);
}

/** point evaluation for EXPR_POLYNOMIAL */
static
SCIP_DECL_EXPREVAL( exprevalPolynomial )
{   /*lint --e{715}*/
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
   SCIP_EXPRDATA_MONOMIAL*   monomialdata;
   SCIP_Real childval;
   SCIP_Real exponent;
   SCIP_Real monomialval;
   int i;
   int j;

   assert(result != NULL);
   assert(argvals != NULL || nargs == 0);
   assert(opdata.data != NULL);

   polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)opdata.data;
   assert(polynomialdata != NULL);

   *result = polynomialdata->constant;

   for( i = 0; i < polynomialdata->nmonomials; ++i )
   {
      monomialdata = polynomialdata->monomials[i];
      assert(monomialdata != NULL);

      monomialval = monomialdata->coef;
      for( j = 0; j < monomialdata->nfactors; ++j )
      {
         assert(monomialdata->childidxs[j] >= 0);
         assert(monomialdata->childidxs[j] < nargs);

         childval = argvals[monomialdata->childidxs[j]];  /*lint !e613*/
         if( childval == 1.0 )  /* 1^anything == 1 */
            continue;

         exponent = monomialdata->exponents[j];

         if( childval == 0.0 )
         {
            if( exponent > 0.0 )
            {
               /* 0^positive == 0 */
               monomialval = 0.0;
               break;
            }
            else if( exponent < 0.0 )
            {
               /* 0^negative = nan (or should it be +inf?, doesn't really matter) */
#ifdef NAN
               *result = NAN;
#else
               /* cppcheck-suppress wrongmathcall */
               *result = pow(0.0, -1.0);
#endif
               return SCIP_OKAY;
            }
            /* 0^0 == 1 */
            continue;
         }

         /* cover some special exponents separately to avoid calling expensive pow function */
         if( exponent == 0.0 )
            continue;
         if( exponent == 1.0 )
         {
            monomialval *= childval;
            continue;
         }
         if( exponent == 2.0 )
         {
            monomialval *= childval * childval;
            continue;
         }
         if( exponent == 0.5 )
         {
            monomialval *= sqrt(childval);
            continue;
         }
         if( exponent == -1.0 )
         {
            monomialval /= childval;
            continue;
         }
         if( exponent == -2.0 )
         {
            monomialval /= childval * childval;
            continue;
         }
         monomialval *= pow(childval, exponent);
      }

      *result += monomialval;
   }

   return SCIP_OKAY;
}

/** interval evaluation for EXPR_POLYNOMIAL */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntPolynomial )
{   /*lint --e{715}*/
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
   SCIP_EXPRDATA_MONOMIAL*   monomialdata;
   SCIP_INTERVAL childval;
   SCIP_INTERVAL monomialval;
   SCIP_Real exponent;
   int i;
   int j;

   assert(result != NULL);
   assert(argvals != NULL || nargs == 0);
   assert(opdata.data != NULL);

   polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)opdata.data;
   assert(polynomialdata != NULL);

   SCIPintervalSet(result, polynomialdata->constant);

   for( i = 0; i < polynomialdata->nmonomials; ++i )
   {
      monomialdata = polynomialdata->monomials[i];
      assert(monomialdata != NULL);

      SCIPintervalSet(&monomialval, monomialdata->coef);
      for( j = 0; j < monomialdata->nfactors && !SCIPintervalIsEntire(infinity, monomialval); ++j )
      {
         assert(monomialdata->childidxs[j] >= 0);
         assert(monomialdata->childidxs[j] < nargs);

         childval = argvals[monomialdata->childidxs[j]];  /*lint !e613*/

         exponent = monomialdata->exponents[j];

         /* cover some special exponents separately to avoid calling expensive pow function */
         if( exponent == 0.0 )
            continue;

         if( exponent == 1.0 )
         {
            SCIPintervalMul(infinity, &monomialval, monomialval, childval);
            continue;
         }

         if( exponent == 2.0 )
         {
            SCIPintervalSquare(infinity, &childval, childval);
            SCIPintervalMul(infinity, &monomialval, monomialval, childval);
            continue;
         }

         if( exponent == 0.5 )
         {
            SCIPintervalSquareRoot(infinity, &childval, childval);
            if( SCIPintervalIsEmpty(infinity, childval) )
            {
               SCIPintervalSetEmpty(result);
               break;
            }
            SCIPintervalMul(infinity, &monomialval, monomialval, childval);
            continue;
         }
         else if( exponent == -1.0 )
         {
            SCIPintervalDiv(infinity, &monomialval, monomialval, childval);
         }
         else if( exponent == -2.0 )
         {
            SCIPintervalSquare(infinity, &childval, childval);
            SCIPintervalDiv(infinity, &monomialval, monomialval, childval);
         }
         else
         {
            SCIPintervalPowerScalar(infinity, &childval, childval, exponent);
            if( SCIPintervalIsEmpty(infinity, childval) )
            {
               SCIPintervalSetEmpty(result);
               return SCIP_OKAY;
            }
            SCIPintervalMul(infinity, &monomialval, monomialval, childval);
         }

         /* the cases in which monomialval gets empty should have been catched */
         assert(!SCIPintervalIsEmpty(infinity, monomialval));
      }

      SCIPintervalAdd(infinity, result, *result, monomialval);
   }

   return SCIP_OKAY;
}

/** curvature for EXPR_POLYNOMIAL */
static
SCIP_DECL_EXPRCURV( exprcurvPolynomial )
{   /*lint --e{715}*/
   SCIP_EXPRDATA_POLYNOMIAL* data;
   SCIP_EXPRDATA_MONOMIAL** monomials;
   SCIP_EXPRDATA_MONOMIAL* monomial;
   int nmonomials;
   int i;

   assert(result    != NULL);
   assert(argcurv   != NULL);
   assert(argbounds != NULL);

   data = (SCIP_EXPRDATA_POLYNOMIAL*)opdata.data;
   assert(data != NULL);

   monomials  = data->monomials;
   nmonomials = data->nmonomials;

   *result = SCIP_EXPRCURV_LINEAR;

   for( i = 0; i < nmonomials && *result != SCIP_EXPRCURV_UNKNOWN; ++i )
   {
      /* we assume that some simplifier was running, so that monomials do not have constants in their factors and such that all factors are different
       * (result would still be correct)
       */
      monomial = monomials[i];
      *result = SCIPexprcurvAdd(*result, SCIPexprcurvMultiply(monomial->coef, SCIPexprcurvMonomial(monomial->nfactors, monomial->exponents, monomial->childidxs, argcurv, argbounds)));
   }

   return SCIP_OKAY;
}

/** expression data copy for EXPR_POLYNOMIAL */
static
SCIP_DECL_EXPRCOPYDATA( exprCopyDataPolynomial )
{  /*lint --e{715}*/
   SCIP_EXPRDATA_POLYNOMIAL* sourcepolynomialdata;
   SCIP_EXPRDATA_POLYNOMIAL* targetpolynomialdata;

   assert(blkmem != NULL);
   assert(opdatatarget != NULL);

   sourcepolynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)opdatasource.data;
   assert(sourcepolynomialdata != NULL);

   SCIP_CALL( polynomialdataCopy(blkmem, &targetpolynomialdata, sourcepolynomialdata) );

   opdatatarget->data = (void*)targetpolynomialdata;

   return SCIP_OKAY;
}

/** expression data free for EXPR_POLYNOMIAL */
static
SCIP_DECL_EXPRFREEDATA( exprFreeDataPolynomial )
{  /*lint --e{715}*/
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;

   assert(blkmem != NULL);

   polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)opdata.data;
   assert(polynomialdata != NULL);

   polynomialdataFree(blkmem, &polynomialdata);
}

/** point evaluation for user expression */
static
SCIP_DECL_EXPREVAL( exprevalUser )
{  /*lint --e{715}*/
   SCIP_EXPRDATA_USER* exprdata;

   exprdata = (SCIP_EXPRDATA_USER*) opdata.data;

   SCIP_CALL( exprdata->eval(exprdata->userdata, nargs, argvals, result, NULL, NULL) );

   return SCIP_OKAY;
}

/** interval evaluation for user expression */
static
SCIP_DECL_EXPRINTEVAL( exprevalIntUser )
{  /*lint --e{715}*/
   SCIP_EXPRDATA_USER* exprdata;

   exprdata = (SCIP_EXPRDATA_USER*) opdata.data;

   if( exprdata->inteval != NULL )
   {
      SCIP_CALL( exprdata->inteval(infinity, exprdata->userdata, nargs, argvals, result, NULL, NULL) );
   }
   else
   {
      /* if user does not provide interval evaluation, then return a result that is always correct */
      SCIPintervalSetEntire(infinity, result);
   }

   return SCIP_OKAY;
}

/** curvature check for user expression */
static
SCIP_DECL_EXPRCURV( exprcurvUser )
{
   SCIP_EXPRDATA_USER* exprdata;

   exprdata = (SCIP_EXPRDATA_USER*) opdata.data;

   if( exprdata->curv != NULL )
   {
      SCIP_CALL( exprdata->curv(infinity, exprdata->userdata, nargs, argbounds, argcurv, result) );
   }
   else
   {
      /* if user does not provide curvature check, then return unknown (which is handled like indefinite) */
      *result = SCIP_EXPRCURV_UNKNOWN;
   }

   return SCIP_OKAY;
}

/** data copy for user expression */
static
SCIP_DECL_EXPRCOPYDATA( exprCopyDataUser )
{
   SCIP_EXPRDATA_USER* exprdatasource;
   SCIP_EXPRDATA_USER* exprdatatarget;

   assert(blkmem != NULL);
   assert(opdatatarget != NULL);

   exprdatasource = (SCIP_EXPRDATA_USER*)opdatasource.data;
   assert(exprdatasource != NULL);

   /* duplicate expression data */
   SCIP_ALLOC( BMSduplicateBlockMemory(blkmem, &exprdatatarget, exprdatasource) );

   /* duplicate user expression data, if any */
   if( exprdatasource->copydata != NULL )
   {
      SCIP_CALL( exprdatasource->copydata(blkmem, nchildren, exprdatasource->userdata, &exprdatatarget->userdata) );
   }
   else
   {
      /* if no copy function for data, then there has to be no data */
      assert(exprdatatarget->userdata == NULL);
   }

   opdatatarget->data = (void*)exprdatatarget;

   return SCIP_OKAY;
}

/** data free for user expression */
static
SCIP_DECL_EXPRFREEDATA( exprFreeDataUser )
{
   SCIP_EXPRDATA_USER* exprdata;

   assert(blkmem != NULL);

   exprdata = (SCIP_EXPRDATA_USER*)opdata.data;

   /* free user expression data, if any */
   if( exprdata->freedata != NULL )
   {
      exprdata->freedata(blkmem, nchildren, exprdata->userdata);
   }
   else
   {
      assert(exprdata->userdata == NULL);
   }

   /* free expression data */
   BMSfreeBlockMemory(blkmem, &exprdata);
}

/** element in table of expression operands */
struct exprOpTableElement
{
   const char*           name;               /**< name of operand (used for printing) */
   int                   nargs;              /**< number of arguments (negative if not fixed) */
   SCIP_DECL_EXPREVAL    ((*eval));          /**< evaluation function */
   SCIP_DECL_EXPRINTEVAL ((*inteval));       /**< interval evaluation function */
   SCIP_DECL_EXPRCURV    ((*curv));          /**< curvature check function */
   SCIP_DECL_EXPRCOPYDATA ((*copydata));     /**< expression data copy function, or NULL to only opdata union */
   SCIP_DECL_EXPRFREEDATA ((*freedata));     /**< expression data free function, or NULL if nothing to free */
};

#define EXPROPEMPTY {NULL, -1, NULL, NULL, NULL, NULL, NULL}

/** table containing for each operand the name, the number of children, and some evaluation functions */
static
struct exprOpTableElement exprOpTable[] =
   {
      EXPROPEMPTY,
      { "variable",          0, exprevalVar,        exprevalIntVar,        exprcurvVar,        NULL, NULL  },
      { "constant",          0, exprevalConst,      exprevalIntConst,      exprcurvConst,      NULL, NULL  },
      { "parameter",         0, exprevalParam,      exprevalIntParam,      exprcurvParam,      NULL, NULL  },
      EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY,
      { "plus",              2, exprevalPlus,       exprevalIntPlus,       exprcurvPlus,       NULL, NULL  },
      { "minus",             2, exprevalMinus,      exprevalIntMinus,      exprcurvMinus,      NULL, NULL  },
      { "mul",               2, exprevalMult,       exprevalIntMult,       exprcurvMult,       NULL, NULL  },
      { "div",               2, exprevalDiv,        exprevalIntDiv,        exprcurvDiv,        NULL, NULL  },
      { "sqr",               1, exprevalSquare,     exprevalIntSquare,     exprcurvSquare,     NULL, NULL  },
      { "sqrt",              1, exprevalSquareRoot, exprevalIntSquareRoot, exprcurvSquareRoot, NULL, NULL  },
      { "realpower",         1, exprevalRealPower,  exprevalIntRealPower,  exprcurvRealPower,  NULL, NULL  },
      { "intpower",          1, exprevalIntPower,   exprevalIntIntPower,   exprcurvIntPower,   NULL, NULL  },
      { "signpower",         1, exprevalSignPower,  exprevalIntSignPower,  exprcurvSignPower,  NULL, NULL  },
      { "exp",               1, exprevalExp,        exprevalIntExp,        exprcurvExp,        NULL, NULL  },
      { "log",               1, exprevalLog,        exprevalIntLog,        exprcurvLog,        NULL, NULL  },
      { "sin",               1, exprevalSin,        exprevalIntSin,        exprcurvSin,        NULL, NULL  },
      { "cos",               1, exprevalCos,        exprevalIntCos,        exprcurvCos,        NULL, NULL  },
      { "tan",               1, exprevalTan,        exprevalIntTan,        exprcurvTan,        NULL, NULL  },
      /* { "erf",               1, exprevalErf,        exprevalIntErf,        exprcurvErf,        NULL, NULL  }, */
      /* { "erfi",              1, exprevalErfi,       exprevalIntErfi        exprcurvErfi,       NULL, NULL  }, */
      EXPROPEMPTY, EXPROPEMPTY,
      { "min",               2, exprevalMin,        exprevalIntMin,        exprcurvMin,        NULL, NULL  },
      { "max",               2, exprevalMax,        exprevalIntMax,        exprcurvMax,        NULL, NULL  },
      { "abs",               1, exprevalAbs,        exprevalIntAbs,        exprcurvAbs,        NULL, NULL  },
      { "sign",              1, exprevalSign,       exprevalIntSign,       exprcurvSign,       NULL, NULL  },
      EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY,
      EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY,
      EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY,
      EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY,
      EXPROPEMPTY, EXPROPEMPTY, EXPROPEMPTY,
      { "sum",              -2, exprevalSum,        exprevalIntSum,        exprcurvSum,        NULL, NULL  },
      { "prod",             -2, exprevalProduct,    exprevalIntProduct,    exprcurvProduct,    NULL, NULL  },
      { "linear",           -2, exprevalLinear,     exprevalIntLinear,     exprcurvLinear,     exprCopyDataLinear,     exprFreeDataLinear     },
      { "quadratic",        -2, exprevalQuadratic,  exprevalIntQuadratic,  exprcurvQuadratic,  exprCopyDataQuadratic,  exprFreeDataQuadratic  },
      { "polynomial",       -2, exprevalPolynomial, exprevalIntPolynomial, exprcurvPolynomial, exprCopyDataPolynomial, exprFreeDataPolynomial },
      { "user",             -2, exprevalUser,       exprevalIntUser,       exprcurvUser,       exprCopyDataUser,       exprFreeDataUser       }
   };

/**@} */

/**@name Expression operand methods */
/**@{ */

/** gives the name of an operand as string */
const char* SCIPexpropGetName(
   SCIP_EXPROP           op                  /**< expression operand */
   )
{
   assert(op < SCIP_EXPR_LAST);

   return exprOpTable[op].name;
}

/** gives the number of children of a simple operand */
int SCIPexpropGetNChildren(
   SCIP_EXPROP           op                  /**< expression operand */
   )
{
   assert(op < SCIP_EXPR_LAST);

   return exprOpTable[op].nargs;
}

/**@} */

/**@name Expressions private methods */
/**@{ */

/** creates an expression
 *
 *  Note, that the expression is allocated but for the children only the pointer is copied.
 */
static
SCIP_RETCODE exprCreate(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPR**           expr,               /**< pointer to buffer for expression address */
   SCIP_EXPROP           op,                 /**< operand of expression */
   int                   nchildren,          /**< number of children */
   SCIP_EXPR**           children,           /**< children */
   SCIP_EXPROPDATA       opdata              /**< operand data */
   )
{
   assert(blkmem != NULL);
   assert(expr   != NULL);
   assert(children != NULL || nchildren == 0);
   assert(children == NULL || nchildren >  0);

   SCIP_ALLOC( BMSallocBlockMemory(blkmem, expr) );

   (*expr)->op        = op;
   (*expr)->nchildren = nchildren;
   (*expr)->children  = children;
   (*expr)->data      = opdata;

   return SCIP_OKAY;
}

/** tries to convert a given (operator,operatordata) pair into a polynomial operator with corresponding data
 *
 *  Does not do this for constants.
 *  If conversion is not possible or operator is already polynomial, *op and *data are
 *  left untouched.
 */
static
SCIP_RETCODE exprConvertToPolynomial(
   BMS_BLKMEM*           blkmem,             /**< block memory */
   SCIP_EXPROP*          op,                 /**< pointer to expression operator */
   SCIP_EXPROPDATA*      data,               /**< pointer to expression data */
   int                   nchildren           /**< number of children of operator */
   )
{
   assert(blkmem != NULL);
   assert(op != NULL);
   assert(data != NULL);

   switch( *op )
   {
   case SCIP_EXPR_VARIDX:
   case SCIP_EXPR_PARAM:
   case SCIP_EXPR_CONST:
      break;

   case SCIP_EXPR_PLUS:
   {
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* monomials[2];
      int childidx;
      SCIP_Real exponent;

      assert(nchildren == 2);

      /* create monomial for first child */
      childidx = 0;
      exponent = 1.0;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomials[0], 1.0, 1, &childidx, &exponent) );

      /* create monomial for second child */
      childidx = 1;
      exponent = 1.0;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomials[1], 1.0, 1, &childidx, &exponent) );

      /* create polynomial for sum of children */
      SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 2, monomials, 0.0, FALSE) );

      *op = SCIP_EXPR_POLYNOMIAL;
      data->data = (void*)polynomialdata;

      break;
   }

   case SCIP_EXPR_MINUS:
   {
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* monomials[2];
      int childidx;
      SCIP_Real exponent;

      assert(nchildren == 2);

      /* create monomial for first child */
      childidx = 0;
      exponent = 1.0;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomials[0],  1.0, 1, &childidx, &exponent) );

      /* create monomial for second child */
      childidx = 1;
      exponent = 1.0;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomials[1], -1.0, 1, &childidx, &exponent) );

      /* create polynomial for difference of children */
      SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 2, monomials, 0.0, FALSE) );

      *op = SCIP_EXPR_POLYNOMIAL;
      data->data = (void*)polynomialdata;

      break;
   }

   case SCIP_EXPR_MUL:
   {
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* monomial;
      int childidx[2];
      SCIP_Real exponent[2];

      assert(nchildren == 2);

      /* create monomial for product of children */
      childidx[0] = 0;
      childidx[1] = 1;
      exponent[0] = 1.0;
      exponent[1] = 1.0;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 2, childidx, exponent) );

      /* create polynomial */
      SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );

      *op = SCIP_EXPR_POLYNOMIAL;
      data->data = (void*)polynomialdata;

      break;
   }

   case SCIP_EXPR_DIV:
   {
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* monomial;
      int childidx[2];
      SCIP_Real exponent[2];

      assert(nchildren == 2);

      /* create monomial for division of children */
      childidx[0] = 0;
      childidx[1] = 1;
      exponent[0] =  1.0;
      exponent[1] = -1.0;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 2, childidx, exponent) );

      /* create polynomial */
      SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );

      *op = SCIP_EXPR_POLYNOMIAL;
      data->data = (void*)polynomialdata;

      break;
   }

   case SCIP_EXPR_SQUARE:
   {
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* monomial;
      int childidx;
      SCIP_Real exponent;

      assert(nchildren == 1);

      /* create monomial for square of child */
      childidx = 0;
      exponent = 2.0;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &exponent) );

      /* create polynomial */
      SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );

      *op = SCIP_EXPR_POLYNOMIAL;
      data->data = (void*)polynomialdata;

      break;
   }

   case SCIP_EXPR_SQRT:
   {
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* monomial;
      int childidx;
      SCIP_Real exponent;

      assert(nchildren == 1);

      /* create monomial for square root of child */
      childidx = 0;
      exponent = 0.5;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &exponent) );

      /* create polynomial */
      SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );

      *op = SCIP_EXPR_POLYNOMIAL;
      data->data = (void*)polynomialdata;

      break;
   }

   case SCIP_EXPR_REALPOWER:
   {
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* monomial;
      int childidx;

      assert(nchildren == 1);

      /* convert to child0 to the power of exponent */

      /* create monomial for power of first child */
      childidx = 0;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &data->dbl) );

      /* create polynomial */
      SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );

      *op = SCIP_EXPR_POLYNOMIAL;
      data->data = (void*)polynomialdata;

      break;
   }

   case SCIP_EXPR_SIGNPOWER:
   {
      SCIP_Real exponent;

      assert(nchildren == 1);

      /* check if exponent is an odd integer */
      exponent = data->dbl;
      if( EPSISINT(exponent, 0.0) && (int)exponent % 2 != 0 )  /*lint !e835*/
      {
         /* convert to child0 to the power of exponent, since sign is kept by taking power */
         SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
         SCIP_EXPRDATA_MONOMIAL* monomial;
         int childidx;

         /* create monomial for power of first child */
         childidx = 0;
         SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &exponent) );

         /* create polynomial */
         SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );

         *op = SCIP_EXPR_POLYNOMIAL;
         data->data = (void*)polynomialdata;
      }
      /* if exponent is not an odd integer constant, then keep it as signpower expression */
      break;
   }

   case SCIP_EXPR_INTPOWER:
   {
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* monomial;
      int childidx;
      SCIP_Real exponent;

      assert(nchildren == 1);

      /* create monomial for power of child */
      childidx = 0;
      exponent = data->intval;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &exponent) );

      /* create polynomial */
      SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );

      *op = SCIP_EXPR_POLYNOMIAL;
      data->data = (void*)polynomialdata;

      break;
   }

   case SCIP_EXPR_EXP:
   case SCIP_EXPR_LOG:
   case SCIP_EXPR_SIN:
   case SCIP_EXPR_COS:
   case SCIP_EXPR_TAN:
      /* case SCIP_EXPR_ERF: */
      /* case SCIP_EXPR_ERFI: */
   case SCIP_EXPR_MIN:
   case SCIP_EXPR_MAX:
   case SCIP_EXPR_ABS:
   case SCIP_EXPR_SIGN:
   case SCIP_EXPR_USER:
      break;

   case SCIP_EXPR_SUM:
   {
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* monomial;
      int childidx;
      int i;
      SCIP_Real exponent;

      /* create empty polynomial */
      SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 0, NULL, 0.0, FALSE) );
      SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, nchildren) );
      assert(polynomialdata->monomialssize >= nchildren);

      /* add summands as monomials */
      childidx = 0;
      exponent = 1.0;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &exponent) );
      for( i = 0; i < nchildren; ++i )
      {
         monomial->childidxs[0] = i;
         SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, 1, &monomial, TRUE) );
      }
      SCIPexprFreeMonomial(blkmem, &monomial);

      *op = SCIP_EXPR_POLYNOMIAL;
      data->data = (void*)polynomialdata;

      break;
   }

   case SCIP_EXPR_PRODUCT:
   {
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* monomial;
      int childidx;
      int i;
      SCIP_Real exponent;

      /* create monomial */
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 0, NULL, NULL) );
      SCIP_CALL( monomialdataEnsureFactorsSize(blkmem, monomial, nchildren) );
      exponent = 1.0;
      for( i = 0; i < nchildren; ++i )
      {
         childidx = i;
         SCIP_CALL( SCIPexprAddMonomialFactors(blkmem, monomial, 1, &childidx, &exponent) );
      }

      /* create polynomial */
      SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 1, &monomial, 0.0, FALSE) );

      *op = SCIP_EXPR_POLYNOMIAL;
      data->data = (void*)polynomialdata;

      break;
   }

   case SCIP_EXPR_LINEAR:
   {
      SCIP_Real* lineardata;
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* monomial;
      int childidx;
      int i;
      SCIP_Real exponent;

      /* get coefficients of linear term */
      lineardata = (SCIP_Real*)data->data;
      assert(lineardata != NULL);

      /* create polynomial consisting of constant from linear term */
      SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 0, NULL, lineardata[nchildren], FALSE) );
      /* ensure space for linear coefficients */
      SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, nchildren) );
      assert(polynomialdata->monomialssize >= nchildren);

      /* add summands as monomials */
      childidx = 0;
      exponent = 1.0;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &monomial, 1.0, 1, &childidx, &exponent) );
      for( i = 0; i < nchildren; ++i )
      {
         monomial->coef = lineardata[i];
         monomial->childidxs[0] = i;
         SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, 1, &monomial, TRUE) );
      }
      SCIPexprFreeMonomial(blkmem, &monomial);

      /* free linear expression data */
      exprFreeDataLinear(blkmem, nchildren, *data);

      *op = SCIP_EXPR_POLYNOMIAL;
      data->data = (void*)polynomialdata;

      break;
   }

   case SCIP_EXPR_QUADRATIC:
   {
      SCIP_EXPRDATA_QUADRATIC* quaddata;
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* squaremonomial;
      SCIP_EXPRDATA_MONOMIAL* bilinmonomial;
      SCIP_EXPRDATA_MONOMIAL* linmonomial;
      int childidx[2];
      SCIP_Real exponent[2];
      int i;

      /* get data of quadratic expression */
      quaddata = (SCIP_EXPRDATA_QUADRATIC*)data->data;
      assert(quaddata != NULL);

      /* create empty polynomial */
      SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 0, NULL, quaddata->constant, FALSE) );
      /* ensure space for linear and quadratic terms */
      SCIP_CALL( polynomialdataEnsureMonomialsSize(blkmem, polynomialdata, (quaddata->lincoefs != NULL ? nchildren : 0) + quaddata->nquadelems) );
      assert(polynomialdata->monomialssize >= quaddata->nquadelems);

      childidx[0] = 0;
      childidx[1] = 0;

      /* create monomial templates */
      exponent[0] = 2.0;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &squaremonomial, 1.0, 1, childidx, exponent) );
      exponent[0] = 1.0;
      exponent[1] = 1.0;
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &bilinmonomial,  1.0, 2, childidx, exponent) );
      SCIP_CALL( SCIPexprCreateMonomial(blkmem, &linmonomial,    1.0, 1, childidx, exponent) );

      /* add linear terms as monomials */
      if( quaddata->lincoefs != NULL )
         for( i = 0; i < nchildren; ++i )
            if( quaddata->lincoefs[i] != 0.0 )
            {
               linmonomial->childidxs[0] = i;
               linmonomial->coef         = quaddata->lincoefs[i];
               SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, 1, &linmonomial, TRUE) );
            }

      /* add quadratic terms as monomials */
      for( i = 0; i < quaddata->nquadelems; ++i )
      {
         if( quaddata->quadelems[i].idx1 == quaddata->quadelems[i].idx2 )
         {
            squaremonomial->childidxs[0] = quaddata->quadelems[i].idx1;
            squaremonomial->coef         = quaddata->quadelems[i].coef;
            SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, 1, &squaremonomial, TRUE) );
         }
         else
         {
            bilinmonomial->childidxs[0] = quaddata->quadelems[i].idx1;
            bilinmonomial->childidxs[1] = quaddata->quadelems[i].idx2;
            bilinmonomial->coef         = quaddata->quadelems[i].coef;
            SCIP_CALL( polynomialdataAddMonomials(blkmem, polynomialdata, 1, &bilinmonomial, TRUE) );
         }
      }
      SCIPexprFreeMonomial(blkmem, &squaremonomial);
      SCIPexprFreeMonomial(blkmem, &bilinmonomial);
      SCIPexprFreeMonomial(blkmem, &linmonomial);

      /* free quadratic expression data */
      exprFreeDataQuadratic(blkmem, nchildren, *data);

      *op = SCIP_EXPR_POLYNOMIAL;
      data->data = (void*)polynomialdata;

      break;
   }

   case SCIP_EXPR_POLYNOMIAL:
   case SCIP_EXPR_LAST:
      break;
   }  /*lint !e788*/

   return SCIP_OKAY;
}

/** converts polynomial expression back into simpler expression, if possible */
static
SCIP_RETCODE exprUnconvertPolynomial(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPROP*          op,                 /**< pointer to expression operator */
   SCIP_EXPROPDATA*      data,               /**< pointer to expression data holding polynomial data */
   int                   nchildren,          /**< number of children of operator */
   void**                children            /**< children array */
   )
{
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
   SCIP_EXPRDATA_MONOMIAL* monomial;
   int maxdegree;
   int nlinmonomials;
   int i;
   int j;

   assert(blkmem != NULL);
   assert(op != NULL);
   assert(*op == SCIP_EXPR_POLYNOMIAL);
   assert(data != NULL);
   assert(children != NULL || nchildren == 0);

   polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)data->data;
   assert(polynomialdata != NULL);

   /* make sure monomials are sorted and merged */
   polynomialdataMergeMonomials(blkmem, polynomialdata, 0.0, TRUE);

   /* if no monomials, then leave as it is */
   if( polynomialdata->nmonomials == 0 )
      return SCIP_OKAY;

   /* check maximal degree of polynomial only - not considering children expressions
    * check number of linear monomials */
   maxdegree = 0;
   nlinmonomials = 0;
   for( i = 0; i < polynomialdata->nmonomials; ++i )
   {
      int monomialdegree;

      monomial = polynomialdata->monomials[i];
      assert(monomial != NULL);

      monomialdegree = 0;
      for(j = 0; j < monomial->nfactors; ++j )
      {
         if( !EPSISINT(monomial->exponents[j], 0.0) || monomial->exponents[j] < 0.0 )  /*lint !e835*/
         {
            monomialdegree = SCIP_EXPR_DEGREEINFINITY;
            break;
         }

         monomialdegree += (int)EPSROUND(monomial->exponents[j], 0.0);  /*lint !e835*/
      }

      if( monomialdegree == SCIP_EXPR_DEGREEINFINITY )
      {
         maxdegree = SCIP_EXPR_DEGREEINFINITY;
         break;
      }

      if( monomialdegree == 1 )
         ++nlinmonomials;

      if( monomialdegree > maxdegree )
         maxdegree = monomialdegree;
   }
   assert(maxdegree > 0 );

   if( maxdegree == 1 )
   {
      /* polynomial is a linear expression in children */

      /* polynomial simplification and monomial merging should ensure that monomial i corresponds to child i and that there are not unused children */
      assert(polynomialdata->nmonomials == nchildren);
      assert(polynomialdata->nmonomials == nlinmonomials);

      if( polynomialdata->constant == 0.0 && polynomialdata->nmonomials == 2 && polynomialdata->monomials[0]->coef == 1.0 && polynomialdata->monomials[1]->coef == 1.0 )
      {
         /* polynomial is addition of two expressions, so turn into SCIP_EXPR_PLUS */
         assert(polynomialdata->monomials[0]->nfactors == 1);
         assert(polynomialdata->monomials[0]->exponents[0] == 1.0);
         assert(polynomialdata->monomials[1]->nfactors == 1);
         assert(polynomialdata->monomials[1]->exponents[0] == 1.0);

         polynomialdataFree(blkmem, &polynomialdata);
         data->data = NULL;

         /* change operator type to PLUS */
         *op = SCIP_EXPR_PLUS;

         return SCIP_OKAY;
      }

      if( polynomialdata->constant == 0.0 && polynomialdata->nmonomials == 2 && polynomialdata->monomials[0]->coef == 1.0 && polynomialdata->monomials[1]->coef == -1.0 )
      {
         /* polynomial is substraction of two expressions, so turn into SCIP_EXPR_MINUS */
         assert(polynomialdata->monomials[0]->nfactors == 1);
         assert(polynomialdata->monomials[0]->exponents[0] == 1.0);
         assert(polynomialdata->monomials[1]->nfactors == 1);
         assert(polynomialdata->monomials[1]->exponents[0] == 1.0);

         polynomialdataFree(blkmem, &polynomialdata);
         data->data = NULL;

         /* change operator type to MINUS */
         *op = SCIP_EXPR_MINUS;

         return SCIP_OKAY;
      }

      if( polynomialdata->constant == 0.0 && polynomialdata->nmonomials == 2 && polynomialdata->monomials[0]->coef == -1.0 && polynomialdata->monomials[1]->coef == 1.0 )
      {
         /* polynomial is substraction of two expressions, so turn into SCIP_EXPR_MINUS */
         void* tmp;

         assert(polynomialdata->monomials[0]->nfactors == 1);
         assert(polynomialdata->monomials[0]->exponents[0] == 1.0);
         assert(polynomialdata->monomials[1]->nfactors == 1);
         assert(polynomialdata->monomials[1]->exponents[0] == 1.0);

         polynomialdataFree(blkmem, &polynomialdata);
         data->data = NULL;

         /* swap children */
         tmp = children[1];           /*lint !e613*/
         children[1] = children[0];   /*lint !e613*/
         children[0] = tmp;           /*lint !e613*/

         /* change operator type to MINUS */
         *op = SCIP_EXPR_MINUS;

         return SCIP_OKAY;
      }

      if( polynomialdata->constant == 0.0 )
      {
         /* check if all monomials have coefficient 1.0 */
         for( i = 0; i < polynomialdata->nmonomials; ++i )
            if( polynomialdata->monomials[i]->coef != 1.0 )
               break;

         if( i == polynomialdata->nmonomials )
         {
            /* polynomial is sum of children, so turn into SCIP_EXPR_SUM */

            polynomialdataFree(blkmem, &polynomialdata);
            data->data = NULL;

            /* change operator type to MINUS */
            *op = SCIP_EXPR_SUM;

            return SCIP_OKAY;
         }
      }

      /* turn polynomial into linear expression */
      {
         SCIP_Real* lindata;

         /* monomial merging should ensure that each child appears in at most one monomial,
          * that monomials are ordered according to the child index, and that constant monomials have been removed
          */

         /* setup data of linear expression */
         SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &lindata, polynomialdata->nmonomials + 1) );

         for( i = 0; i < polynomialdata->nmonomials; ++i )
         {
            assert(polynomialdata->monomials[i]->childidxs[0] == i);
            assert(polynomialdata->monomials[i]->exponents[0] == 1.0);
            lindata[i] = polynomialdata->monomials[i]->coef;  /*lint !e644*/
         }
         lindata[i] = polynomialdata->constant;

         polynomialdataFree(blkmem, &polynomialdata);
         *op = SCIP_EXPR_LINEAR;
         data->data = (void*)lindata;

         return SCIP_OKAY;
      }
   }

   if( maxdegree == 2 && (polynomialdata->nmonomials > 1 || polynomialdata->constant != 0.0 || polynomialdata->monomials[0]->coef != 1.0) )
   {
      /* 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 */
      SCIP_EXPRDATA_QUADRATIC* quaddata;
      int quadelemidx;

      SCIP_ALLOC( BMSallocBlockMemory(blkmem, &quaddata) );
      SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &quaddata->quadelems, polynomialdata->nmonomials - nlinmonomials) );
      quaddata->nquadelems = polynomialdata->nmonomials - nlinmonomials;
      quaddata->constant = polynomialdata->constant;
      quaddata->sorted = FALSE; /* quadratic data is sorted different than polynomials */

      if( nlinmonomials > 0 )
      {
         SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &quaddata->lincoefs, nchildren) );
         BMSclearMemoryArray(quaddata->lincoefs, nchildren);
      }
      else
         quaddata->lincoefs = NULL;

      quadelemidx = 0;
      for( i = 0; i < polynomialdata->nmonomials; ++i )
      {
         assert(polynomialdata->monomials[i]->nfactors == 1 || polynomialdata->monomials[i]->nfactors == 2);
         if( polynomialdata->monomials[i]->nfactors == 1 )
         {
            if( polynomialdata->monomials[i]->exponents[0] == 1.0 )
            {
               /* monomial is a linear term */
               assert(quaddata->lincoefs != NULL);
               quaddata->lincoefs[polynomialdata->monomials[i]->childidxs[0]] += polynomialdata->monomials[i]->coef;
            }
            else
            {
               /* monomial should be a square term */
               assert(polynomialdata->monomials[i]->exponents[0] == 2.0);
               assert(quadelemidx < quaddata->nquadelems);
               quaddata->quadelems[quadelemidx].idx1 = polynomialdata->monomials[i]->childidxs[0];
               quaddata->quadelems[quadelemidx].idx2 = polynomialdata->monomials[i]->childidxs[0];
               quaddata->quadelems[quadelemidx].coef = polynomialdata->monomials[i]->coef;
               ++quadelemidx;
            }
         }
         else
         {
            /* monomial should be a bilinear term */
            assert(polynomialdata->monomials[i]->exponents[0] == 1.0);
            assert(polynomialdata->monomials[i]->exponents[1] == 1.0);
            assert(quadelemidx < quaddata->nquadelems);
            quaddata->quadelems[quadelemidx].idx1 = MIN(polynomialdata->monomials[i]->childidxs[0], polynomialdata->monomials[i]->childidxs[1]);
            quaddata->quadelems[quadelemidx].idx2 = MAX(polynomialdata->monomials[i]->childidxs[0], polynomialdata->monomials[i]->childidxs[1]);
            quaddata->quadelems[quadelemidx].coef = polynomialdata->monomials[i]->coef;
            ++quadelemidx;
         }
      }
      assert(quadelemidx == quaddata->nquadelems);

      polynomialdataFree(blkmem, &polynomialdata);

      *op = SCIP_EXPR_QUADRATIC;
      data->data = (void*)quaddata;

      return SCIP_OKAY;
   }

   if( polynomialdata->constant == 0.0 && polynomialdata->nmonomials == 1 && polynomialdata->monomials[0]->coef == 1.0 )
   {
      /* polynomial is product of children */
      monomial = polynomialdata->monomials[0];
      assert(monomial->nfactors == nchildren);

      if( monomial->nfactors == 1 )
      {
         /* polynomial is x^k for some k */
         assert(monomial->exponents[0] != 1.0); /* should have been handled before */
         assert(monomial->childidxs[0] == 0);

         if( monomial->exponents[0] == 2.0 )
         {
            /* polynomial is x^2, so turn into SCIP_EXPR_SQUARE */

            polynomialdataFree(blkmem, &polynomialdata);
            data->data = NULL;

            *op = SCIP_EXPR_SQUARE;

            return SCIP_OKAY;
         }

         if( EPSISINT(monomial->exponents[0], 0.0) )  /*lint !e835*/
         {
            /* k is an integer, so turn into SCIP_EXPR_INTPOWER */
            int exponent;

            exponent = (int)EPSROUND(monomial->exponents[0], 0.0);  /*lint !e835*/

            polynomialdataFree(blkmem, &polynomialdata);

            *op = SCIP_EXPR_INTPOWER;
            data->intval = exponent;

            return SCIP_OKAY;
         }

         if( monomial->exponents[0] == 0.5 )
         {
            /* polynomial is sqrt(x), so turn into SCIP_EXPR_SQRT */

            polynomialdataFree(blkmem, &polynomialdata);
            data->data = NULL;

            *op = SCIP_EXPR_SQRT;

            return SCIP_OKAY;
         }

         {
            /* polynomial is x^a with a some real number, so turn into SCIP_EXPR_REALPOWER */
            SCIP_Real exponent;

            exponent = monomial->exponents[0];

            polynomialdataFree(blkmem, &polynomialdata);

            *op = SCIP_EXPR_REALPOWER;
            data->dbl = exponent;

            return SCIP_OKAY;
         }
      }

      if( maxdegree == 2 && monomial->nfactors == 2 )
      {
         /* polynomial is product of two children, so turn into SCIP_EXPR_MUL */
         assert(monomial->exponents[0] == 1.0);
         assert(monomial->exponents[1] == 1.0);

         polynomialdataFree(blkmem, &polynomialdata);
         data->data = NULL;

         *op = SCIP_EXPR_MUL;

         return SCIP_OKAY;
      }

      if( maxdegree == monomial->nfactors )
      {
         /* polynomial is a product of n children, so turn into SCIP_EXPR_PRODUCT */

         polynomialdataFree(blkmem, &polynomialdata);
         data->data = NULL;

         *op = SCIP_EXPR_PRODUCT;

         return SCIP_OKAY;
      }

      if( monomial->nfactors == 2 && monomial->exponents[0] == 1.0 && monomial->exponents[1] == -1.0 )
      {
         /* polynomial is x/y, so turn into SCIP_EXPR_DIV */
         assert(monomial->childidxs[0] == 0);
         assert(monomial->childidxs[1] == 1);

         polynomialdataFree(blkmem, &polynomialdata);
         data->data = NULL;

         *op = SCIP_EXPR_DIV;

         return SCIP_OKAY;
      }

      if( monomial->nfactors == 2 && monomial->exponents[0] == -1.0 && monomial->exponents[1] == 1.0 )
      {
         /* polynomial is y/x, so turn into SCIP_EXPR_DIV */
         void* tmp;

         assert(monomial->childidxs[0] == 0);
         assert(monomial->childidxs[1] == 1);

         polynomialdataFree(blkmem, &polynomialdata);
         data->data = NULL;

         /* swap children */
         tmp = children[1];           /*lint !e613*/
         children[1] = children[0];   /*lint !e613*/
         children[0] = tmp;           /*lint !e613*/

         *op = SCIP_EXPR_DIV;

         return SCIP_OKAY;
      }
   }

   return SCIP_OKAY;
}

/** adds copies of expressions to the array of children of a sum, product, linear, quadratic, or polynomial expression
 *
 *  For a sum or product expression, this corresponds to add additional summands and factors, resp.
 *  For a linear expression, this corresponds to add each expression with coefficient 1.0.
 *  For a quadratic or polynomial expression, only the children array may be enlarged, the expression itself remains the same.
 */
static
SCIP_RETCODE exprsimplifyAddChildren(
   BMS_BLKMEM*           blkmem,             /**< block memory */
   SCIP_EXPR*            expr,               /**< quadratic or polynomial expression */
   int                   nexprs,             /**< number of expressions to add */
   SCIP_EXPR**           exprs,              /**< expressions to add */
   SCIP_Bool             comparechildren,    /**< whether to compare expressions with already existing children (no effect for sum and product) */
   SCIP_Real             eps,                /**< which epsilon to use when comparing expressions */
   int*                  childmap            /**< array where to store mapping of indices from exprs to children array in expr, or NULL if not of interest */
   )
{
   int i;

   assert(blkmem != NULL);
   assert(expr != NULL);
   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);
   assert(exprs != NULL || nexprs == 0);

   if( nexprs == 0 )
      return SCIP_OKAY;

   switch( expr->op )
   {
   case SCIP_EXPR_SUM:
   case SCIP_EXPR_PRODUCT:
   {
      SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, expr->nchildren, expr->nchildren + nexprs) );
      for( i = 0; i < nexprs; ++i )
      {
         SCIP_CALL( SCIPexprCopyDeep(blkmem, &expr->children[expr->nchildren + i], exprs[i]) );  /*lint !e613*/
         if( childmap != NULL )
            childmap[i] = expr->nchildren + i;
      }
      expr->nchildren += nexprs;

      break;
   }

   case SCIP_EXPR_LINEAR:
   case SCIP_EXPR_QUADRATIC:
   case SCIP_EXPR_POLYNOMIAL:
   {
      int j;
      int orignchildren;
      SCIP_Bool existsalready;

      orignchildren = expr->nchildren;
      SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, expr->nchildren, expr->nchildren + nexprs) );

      for( i = 0; i < nexprs; ++i )
      {
         existsalready = FALSE;
         if( comparechildren )
            for( j = 0; j < orignchildren; ++j )
               /* during simplification of polynomials, their may be NULL's in children array */
               if( expr->children[j] != NULL && SCIPexprAreEqual(expr->children[j], exprs[i], eps) )  /*lint !e613*/
               {
                  existsalready = TRUE;
                  break;
               }

         if( !existsalready )
         {
            /* add copy of exprs[j] to children array */
            SCIP_CALL( SCIPexprCopyDeep(blkmem, &expr->children[expr->nchildren], exprs[i]) );  /*lint !e613*/
            if( childmap != NULL )
               childmap[i] = expr->nchildren;
            ++expr->nchildren;
         }
         else
         {
            if( childmap != NULL )
               childmap[i] = j;  /*lint !e644*/
            if( expr->op == SCIP_EXPR_LINEAR )
            {
               /* if linear expression, increase coefficient by 1.0 */
               ((SCIP_Real*)expr->data.data)[j] += 1.0;
            }
         }
      }

      /* shrink children array to actually used size */
      assert(comparechildren || expr->nchildren == orignchildren + nexprs);
      SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, orignchildren + nexprs, expr->nchildren) );

      if( expr->op == SCIP_EXPR_LINEAR && expr->nchildren > orignchildren )
      {
         /* if linear expression, then add 1.0 coefficients for new expressions */
         SCIP_Real* data;

         data = (SCIP_Real*)expr->data.data;
         SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &data, orignchildren + 1, expr->nchildren + 1) );
         data[expr->nchildren] = data[orignchildren]; /* move constant from old end to new end */
         for( i = orignchildren; i < expr->nchildren; ++i )
            data[i] = 1.0;
         expr->data.data = (void*)data;
      }
      else if( expr->op == SCIP_EXPR_QUADRATIC && expr->nchildren > orignchildren )
      {
         /* if quadratic expression, then add 0.0 linear coefficients for new expressions */
         SCIP_EXPRDATA_QUADRATIC* data;

         data = (SCIP_EXPRDATA_QUADRATIC*)expr->data.data;
         SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &data->lincoefs, orignchildren, expr->nchildren) );
         BMSclearMemoryArray(&data->lincoefs[orignchildren], expr->nchildren - orignchildren);  /*lint !e866*/
      }

      break;
   }

   default:
      SCIPerrorMessage("exprsimplifyAddChildren cannot be called for operand %d\n", expr->op);
      return SCIP_INVALIDDATA;
   }  /*lint !e788*/

   return SCIP_OKAY;
}

/** converts expressions into polynomials, where possible and obvious */
static
SCIP_RETCODE exprsimplifyConvertToPolynomials(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPR*            expr                /**< expression to convert */
   )
{
   int i;

   assert(expr != NULL);

   for( i = 0; i < expr->nchildren; ++i )
   {
      SCIP_CALL( exprsimplifyConvertToPolynomials(blkmem, expr->children[i]) );
   }

   SCIP_CALL( exprConvertToPolynomial(blkmem, &expr->op, &expr->data, expr->nchildren) );

   return SCIP_OKAY;
}

/** removes duplicate children in a polynomial expression
 *
 *  Leaves NULL's in children array.
 */
static
SCIP_RETCODE exprsimplifyRemoveDuplicatePolynomialChildren(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPR*            expr,               /**< expression */
   SCIP_Real             eps                 /**< threshold for zero */
   )
{
   SCIP_Bool foundduplicates;
   int* childmap;
   int i;
   int j;

   assert(blkmem != NULL);
   assert(expr != NULL);
   assert(SCIPexprGetOperator(expr) == SCIP_EXPR_POLYNOMIAL);

   if( expr->nchildren == 0 )
      return SCIP_OKAY;

   SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &childmap, expr->nchildren) );

   foundduplicates = FALSE;
   for( i = 0; i < expr->nchildren; ++i )
   {
      if( expr->children[i] == NULL )
         continue;
      childmap[i] = i;  /*lint !e644*/

      for( j = i+1; j < expr->nchildren; ++j )
      {
         if( expr->children[j] == NULL )
            continue;

         if( SCIPexprAreEqual(expr->children[i], expr->children[j], eps) )
         {
            /* forget about expr j and remember that is to be replaced by i */
            SCIPexprFreeDeep(blkmem, &expr->children[j]);
            childmap[j] = i;
            foundduplicates = TRUE;
         }
      }
   }

   /* apply childmap to monomials */
   if( foundduplicates )
      polynomialdataApplyChildmap((SCIP_EXPRDATA_POLYNOMIAL*)expr->data.data, childmap);

   /* free childmap */
   BMSfreeBlockMemoryArray(blkmem, &childmap, expr->nchildren);

   return SCIP_OKAY;
}

/** eliminates NULL's in children array and shrinks it to actual size */
static
SCIP_RETCODE exprsimplifyRemovePolynomialNullChildren(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPR*            expr                /**< expression */
   )
{
   int* childmap;
   int lastnonnull;
   int i;

   assert(blkmem != NULL);
   assert(expr != NULL);
   assert(SCIPexprGetOperator(expr) == SCIP_EXPR_POLYNOMIAL);

   if( expr->nchildren == 0 )
      return SCIP_OKAY;

   SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &childmap, expr->nchildren) );

   /* close gaps in children array */
   lastnonnull = expr->nchildren-1;
   while( lastnonnull >= 0 && expr->children[lastnonnull] == NULL )
      --lastnonnull;
   for( i = 0; i <= lastnonnull; ++i )
   {
      if( expr->children[i] != NULL )
      {
         childmap[i] = i; /* child at index i is not moved */  /*lint !e644*/
         continue;
      }
      assert(expr->children[lastnonnull] != NULL);

      /* move child at lastnonnull to position i */
      expr->children[i] = expr->children[lastnonnull];
      expr->children[lastnonnull] = NULL;
      childmap[lastnonnull] = i;

      /* update lastnonnull */
      --lastnonnull;
      while( lastnonnull >= 0 && expr->children[lastnonnull] == NULL )
         --lastnonnull;
   }
   assert(i > lastnonnull);

   /* apply childmap to monomials */
   if( lastnonnull < expr->nchildren-1 )
      polynomialdataApplyChildmap((SCIP_EXPRDATA_POLYNOMIAL*)expr->data.data, childmap);

   BMSfreeBlockMemoryArray(blkmem, &childmap, expr->nchildren);

   /* shrink children array */
   if( lastnonnull >= 0 )
   {
      SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &expr->children, expr->nchildren, lastnonnull+1) );
      expr->nchildren = lastnonnull+1;
   }
   else
   {
      BMSfreeBlockMemoryArray(blkmem, &expr->children, expr->nchildren);
      expr->nchildren = 0;
   }

   return SCIP_OKAY;
}

/** checks which children are still in use and frees those which are not */
static
SCIP_RETCODE exprsimplifyRemovePolynomialUnusedChildren(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPR*            expr                /**< polynomial expression */
   )
{
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
   SCIP_EXPRDATA_MONOMIAL* monomial;
   SCIP_Bool* childinuse;
   int i;
   int j;

   assert(blkmem != NULL);
   assert(expr != NULL);

   if( expr->nchildren == 0 )
      return SCIP_OKAY;

   polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)expr->data.data;
   assert(polynomialdata != NULL);

   /* check which children are still in use */
   SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &childinuse, expr->nchildren) );
   BMSclearMemoryArray(childinuse, expr->nchildren);  /*lint !e644*/
   for( i = 0; i < polynomialdata->nmonomials; ++i )
   {
      monomial = polynomialdata->monomials[i];
      assert(monomial != NULL);

      for( j = 0; j < monomial->nfactors; ++j )
      {
         assert(monomial->childidxs[j] >= 0);
         assert(monomial->childidxs[j] < expr->nchildren);
         childinuse[monomial->childidxs[j]] = TRUE;
      }
   }

   /* free children that are not used in any monomial */
   for( i = 0; i < expr->nchildren; ++i )
      if( expr->children[i] != NULL && !childinuse[i] )
         SCIPexprFreeDeep(blkmem, &expr->children[i]);

   BMSfreeBlockMemoryArray(blkmem, &childinuse, expr->nchildren);

   return SCIP_OKAY;
}

/** flattens polynomials in polynomials, check for constants in non-polynomials expressions
 *
 *  exprsimplifyConvertToPolynomials should have been called before to eliminate simple polynomial operands.
 */
static
SCIP_RETCODE exprsimplifyFlattenPolynomials(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_MESSAGEHDLR*     messagehdlr,        /**< message handler */
   SCIP_EXPR*            expr,               /**< expression */
   SCIP_Real             eps,                /**< threshold, under which values are treat as 0 */
   int                   maxexpansionexponent/**< maximal exponent for which we still expand non-monomial polynomials */
   )
{
   int i;

   assert(expr != NULL);

   for( i = 0; i < expr->nchildren; ++i )
   {
      SCIP_CALL( exprsimplifyFlattenPolynomials(blkmem, messagehdlr, expr->children[i], eps, maxexpansionexponent) );
   }

   switch( SCIPexprGetOperator(expr) )
   {
   case SCIP_EXPR_VARIDX:
   case SCIP_EXPR_CONST:
   case SCIP_EXPR_PARAM:
   case SCIP_EXPR_PLUS:
   case SCIP_EXPR_MINUS:
   case SCIP_EXPR_MUL:
   case SCIP_EXPR_DIV:
   case SCIP_EXPR_SQUARE:
   case SCIP_EXPR_SQRT:
   case SCIP_EXPR_INTPOWER:
   case SCIP_EXPR_REALPOWER:
   case SCIP_EXPR_SIGNPOWER:
      break;

   case SCIP_EXPR_EXP:
   case SCIP_EXPR_LOG:
   case SCIP_EXPR_SIN:
   case SCIP_EXPR_COS:
   case SCIP_EXPR_TAN:
      /* case SCIP_EXPR_ERF: */
      /* case SCIP_EXPR_ERFI: */
   case SCIP_EXPR_ABS:
   case SCIP_EXPR_SIGN:
   {
      /* check if argument is a constant */
      if( (expr->children[0]->op == SCIP_EXPR_POLYNOMIAL && SCIPexprGetNChildren(expr->children[0]) == 0) ||
         expr->children[0]->op == SCIP_EXPR_CONST )
      {
         SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
         SCIP_Real exprval;

         /* since child0 has no children and it's polynomial was flattened, it should have no monomials */
         assert(expr->children[0]->op != SCIP_EXPR_POLYNOMIAL || SCIPexprGetNMonomials(expr->children[0]) == 0);

         /* evaluate expression in constant polynomial */
         SCIP_CALL( SCIPexprEval(expr, NULL, NULL, &exprval) );

         /* create polynomial */
         SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 0, NULL, exprval, FALSE) );

         expr->op = SCIP_EXPR_POLYNOMIAL;
         expr->data.data = (void*)polynomialdata;

         /* forget child */
         SCIPexprFreeDeep(blkmem, &expr->children[0]);
         BMSfreeBlockMemoryArray(blkmem, &expr->children, 1);
         expr->nchildren = 0;
      }

      break;
   }

   case SCIP_EXPR_MIN:
   case SCIP_EXPR_MAX:
   {
      /* check if both arguments are constants */
      if( ((expr->children[0]->op == SCIP_EXPR_POLYNOMIAL && SCIPexprGetNChildren(expr->children[0]) == 0) || expr->children[0]->op == SCIP_EXPR_CONST) &&
         ((expr->children[1]->op == SCIP_EXPR_POLYNOMIAL && SCIPexprGetNChildren(expr->children[1]) == 0) || expr->children[1]->op == SCIP_EXPR_CONST) )
      {
         SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
         SCIP_Real exprval;

         /* since children have no children and it's polynomial was flattened, it should have no monomials */
         assert(expr->children[0]->op != SCIP_EXPR_POLYNOMIAL || SCIPexprGetNMonomials(expr->children[0]) == 0);
         assert(expr->children[1]->op != SCIP_EXPR_POLYNOMIAL || SCIPexprGetNMonomials(expr->children[1]) == 0);

         /* evaluate expression in constants */
         SCIP_CALL( SCIPexprEval(expr, NULL, NULL, &exprval) );

         /* create polynomial */
         SCIP_CALL( polynomialdataCreate(blkmem, &polynomialdata, 0, NULL, exprval, FALSE) );

         expr->op = SCIP_EXPR_POLYNOMIAL;
         expr->data.data = (void*)polynomialdata;

         /* forget children */
         SCIPexprFreeDeep(blkmem, &expr->children[0]);
         SCIPexprFreeDeep(blkmem, &expr->children[1]);
         BMSfreeBlockMemoryArray(blkmem, &expr->children, 2);
         expr->nchildren = 0;
      }

      break;
   }

   case SCIP_EXPR_SUM:
   case SCIP_EXPR_PRODUCT:
   case SCIP_EXPR_LINEAR:
   case SCIP_EXPR_QUADRATIC:
   case SCIP_EXPR_USER:
      break;

   case SCIP_EXPR_POLYNOMIAL:
   {
      SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
      SCIP_EXPRDATA_MONOMIAL* monomial;
      SCIP_Bool removechild;
      int* childmap;
      int childmapsize;
      int j;

      /* simplify current polynomial */
      SCIP_CALL( exprsimplifyRemoveDuplicatePolynomialChildren(blkmem, expr, eps) );
      SCIPexprMergeMonomials(blkmem, expr, eps, TRUE);

      polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)expr->data.data;
      assert(polynomialdata != NULL);

      SCIPdebugMessage("expand factors in expression ");
      SCIPdebug( SCIPexprPrint(expr, messagehdlr, NULL, NULL, NULL, NULL) );
      SCIPdebugPrintf("\n");

      childmap = NULL;
      childmapsize = 0;

      /* resolve children that are constants
       * we do this first, because it reduces the degree and number of factors in the monomials,
       *   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
       */
      for( i = 0; i < expr->nchildren; ++i )
      {
         if( expr->children[i] == NULL )
            continue;

         if( SCIPexprGetOperator(expr->children[i]) != SCIP_EXPR_CONST )
            continue;

         removechild = TRUE; /* we intend to delete children[i] */

         if( childmapsize < expr->children[i]->nchildren )
         {
            int newsize;

            newsize = calcGrowSize(expr->children[i]->nchildren);
            SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &childmap, childmapsize, newsize) );
            childmapsize = newsize;
         }

         /* put constant of child i into every monomial where child i is used */
         for( j = 0; j < polynomialdata->nmonomials; ++j )
         {
            int factorpos;
            SCIP_Bool success;

            monomial = polynomialdata->monomials[j];
            /* if monomial is not sorted, then polynomial should not be sorted either, or have only one monomial */
            assert(monomial->sorted || !polynomialdata->sorted || polynomialdata->nmonomials <= 1);

            if( SCIPexprFindMonomialFactor(monomial, i, &factorpos) )
            {
               assert(factorpos >= 0);
               assert(factorpos < monomial->nfactors);
               /* assert that factors have been merged */
               assert(factorpos == 0 || monomial->childidxs[factorpos-1] != i);
               assert(factorpos == monomial->nfactors-1 || monomial->childidxs[factorpos+1] != i);

               /* SCIPdebugMessage("attempt expanding child %d at monomial %d factor %d\n", i, j, factorpos);
                  SCIPdebug( SCIPexprPrint(expr, NULL, NULL, NULL) ); SCIPdebugPrintf("\n");
                  SCIPdebug( SCIPexprPrint(expr->children[i], NULL, NULL, NULL) ); SCIPdebugPrintf("\n"); */

               if( !EPSISINT(monomial->exponents[factorpos], 0.0) && SCIPexprGetOpReal(expr->children[i]) < 0.0 )  /*lint !e835*/
               {
                  /* if constant is negative and our exponent is not integer, then cannot do expansion */
                  SCIPmessagePrintWarning(messagehdlr, "got negative constant %g to the power of a noninteger exponent %g\n",
                     SCIPexprGetOpReal(expr->children[i]), monomial->exponents[factorpos]);
                  success = FALSE;
               }
               else
               {
                  monomial->coef *= pow(SCIPexprGetOpReal(expr->children[i]), monomial->exponents[factorpos]);

                  /* move last factor to position factorpos */
                  if( factorpos < monomial->nfactors-1 )
                  {
                     monomial->exponents[factorpos] = monomial->exponents[monomial->nfactors-1];
                     monomial->childidxs[factorpos] = monomial->childidxs[monomial->nfactors-1];
                  }
                  --monomial->nfactors;
                  monomial->sorted = FALSE;
                  polynomialdata->sorted = FALSE;

                  success = TRUE;
               }

               if( !success )
                  removechild = FALSE;
            }
         }

         /* forget about child i, if it is not used anymore */
         if( removechild )
            SCIPexprFreeDeep(blkmem, &expr->children[i]);

         /* simplify current polynomial again */
         SCIPexprMergeMonomials(blkmem, expr, eps, TRUE);
      }

      /* try to resolve children that are polynomials itself */
      for( i = 0; i < expr->nchildren; ++i )
      {
         if( expr->children[i] == NULL )
            continue;

         if( SCIPexprGetOperator(expr->children[i]) != SCIP_EXPR_POLYNOMIAL )
            continue;

         removechild = TRUE; /* we intend to delete children[i] */

         if( childmapsize < expr->children[i]->nchildren )
         {
            int newsize;

            newsize = calcGrowSize(expr->children[i]->nchildren);
            SCIP_ALLOC( BMSreallocBlockMemoryArray(blkmem, &childmap, childmapsize, newsize) );
            childmapsize = newsize;
         }

         /* add children of child i */
         SCIP_CALL( exprsimplifyAddChildren(blkmem, expr, expr->children[i]->nchildren, expr->children[i]->children, TRUE, eps, childmap) );

         /* put polynomial of child i into every monomial where child i is used */
         j = 0;
         while( j < polynomialdata->nmonomials )
         {
            int factorpos;
            SCIP_Bool success;

            monomial = polynomialdata->monomials[j];
            /* if monomial is not sorted, then polynomial should not be sorted either, or have only one monomial */
            assert(monomial->sorted || !polynomialdata->sorted || polynomialdata->nmonomials <= 1);

            if( SCIPexprFindMonomialFactor(monomial, i, &factorpos) )
            {
               assert(factorpos >= 0);
               assert(factorpos < monomial->nfactors);
               /* assert that factors have been merged */
               assert(factorpos == 0 || monomial->childidxs[factorpos-1] != i);
               assert(factorpos == monomial->nfactors-1 || monomial->childidxs[factorpos+1] != i);

               /* SCIPdebugMessage("attempt expanding child %d at monomial %d factor %d\n", i, j, factorpos);
                  SCIPdebug( SCIPexprPrint(expr, NULL, NULL, NULL) ); SCIPdebugPrintf("\n");
                  SCIPdebug( SCIPexprPrint(expr->children[i], NULL, NULL, NULL) ); SCIPdebugPrintf("\n"); */

               SCIP_CALL( polynomialdataExpandMonomialFactor(blkmem, messagehdlr, polynomialdata, j, factorpos,
                     (SCIP_EXPRDATA_POLYNOMIAL*)expr->children[i]->data.data, childmap, maxexpansionexponent, &success) );

               if( !success )
               {
                  removechild = FALSE;
                  ++j;
               }
            }
            else
               ++j;

            /* expansion may remove monomials[j], move a monomial from the end to position j, or add new monomials to the end of polynomialdata
             * we thus repeat with index j, if a factor was successfully expanded
             */
         }

         /* forget about child i, if it is not used anymore */
         if( removechild )
            SCIPexprFreeDeep(blkmem, &expr->children[i]);

         /* simplify current polynomial again */
         SCIPexprMergeMonomials(blkmem, expr, eps, TRUE);
      }

      BMSfreeBlockMemoryArrayNull(blkmem, &childmap, childmapsize);

      /* free children that are not in use anymore */
      SCIP_CALL( exprsimplifyRemovePolynomialUnusedChildren(blkmem, expr) );

      /* remove NULLs from children array */
      SCIP_CALL( exprsimplifyRemovePolynomialNullChildren(blkmem, expr) );

      /* if no children left, then it's a constant polynomial -> change into EXPR_CONST */
      if( expr->nchildren == 0 )
      {
         SCIP_Real val;

         /* if no children, then it should also have no monomials */
         assert(polynomialdata->nmonomials == 0);

         val = polynomialdata->constant;
         polynomialdataFree(blkmem, &polynomialdata);

         expr->op = SCIP_EXPR_CONST;
         expr->data.dbl = val;
      }

      SCIPdebugMessage("-> ");
      SCIPdebug( SCIPexprPrint(expr, messagehdlr, NULL, NULL, NULL, NULL) );
      SCIPdebugPrintf("\n");

      break;
   }

   case SCIP_EXPR_LAST:
      break;
   }  /*lint !e788*/

   return SCIP_OKAY;
}

/** separates linear monomials from an expression, if it is a polynomial expression
 *
 *  Separates only those linear terms whose variable is not used otherwise in the expression.
 */
static
SCIP_RETCODE exprsimplifySeparateLinearFromPolynomial(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPR*            expr,               /**< expression */
   SCIP_Real             eps,                /**< threshold, under which positive values are treat as 0 */
   int                   nvars,              /**< number of variables in expression */
   int*                  nlinvars,           /**< buffer to store number of linear variables in linear part */
   int*                  linidxs,            /**< array to store indices of variables in expression tree which belong to linear part */
   SCIP_Real*            lincoefs            /**< array to store coefficients of linear part */
   )
{
   SCIP_EXPRDATA_POLYNOMIAL* polynomialdata;
   SCIP_EXPRDATA_MONOMIAL* monomial;
   int* varsusage;
   int* childusage;
   int childidx;
   int i;
   int j;

   assert(blkmem != NULL);
   assert(expr != NULL);
   assert(nlinvars != NULL);
   assert(linidxs != NULL);
   assert(lincoefs != NULL);

   *nlinvars = 0;

   if( SCIPexprGetOperator(expr) != SCIP_EXPR_POLYNOMIAL )
      return SCIP_OKAY;

   if( SCIPexprGetNChildren(expr) == 0 )
      return SCIP_OKAY;

   polynomialdata = (SCIP_EXPRDATA_POLYNOMIAL*)expr->data.data;
   assert(polynomialdata != NULL);

   /* get variable usage */
   SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &varsusage, nvars) );
   BMSclearMemoryArray(varsusage, nvars);  /*lint !e644*/
   SCIPexprGetVarsUsage(expr, varsusage);

   /* get child usage: how often each child is used in the polynomial */
   SCIP_ALLOC( BMSallocBlockMemoryArray(blkmem, &childusage, expr->nchildren) );
   BMSclearMemoryArray(childusage, expr->nchildren);  /*lint !e644*/
   for( i = 0; i < polynomialdata->nmonomials; ++i )
   {
      monomial = polynomialdata->monomials[i];
      assert(monomial != NULL);
      for( j = 0; j < monomial->nfactors; ++j )
      {
         assert(monomial->childidxs[j] >= 0);
         assert(monomial->childidxs[j] < expr->nchildren);
         ++childusage[monomial->childidxs[j]];
      }
   }

   /* move linear monomials out of polynomial */
   for( i = 0; i < polynomialdata->nmonomials; ++i )
   {
      monomial = polynomialdata->monomials[i];
      assert(monomial != NULL);
      if( monomial->nfactors != 1 )
         continue;
      if( monomial->exponents[0] != 1.0 )
         continue;
      childidx = monomial->childidxs[0];
      if( SCIPexprGetOperator(expr->children[childidx]) != SCIP_EXPR_VARIDX )
         continue;

      /* we are at a linear monomial in a variable */
      assert(SCIPexprGetOpIndex(expr->children[childidx]) < nvars);
      if( childusage[childidx] == 1 && varsusage[SCIPexprGetOpIndex(expr->children[childidx])] == 1 )
      {
         /* if the child expression is not used in another monomial (which would due to merging be not linear)
          * and if the variable is not used somewhere else in the tree,
          * then move this monomial into linear part and free child
          */
         linidxs[*nlinvars]  = SCIPexprGetOpIndex(expr->children[childidx]);
         lincoefs[*nlinvars] = monomial->coef;
         ++*nlinvars;

         SCIPexprFreeDeep(blkmem, &expr->children[childidx]);
         monomial->coef = 0.0;
         monomial->nfactors = 0;
      }
   }

   BMSfreeBlockMemoryArray(blkmem, &varsusage, nvars);
   BMSfreeBlockMemoryArray(blkmem, &childusage, expr->nchildren);

   if( *nlinvars > 0 )
   {
      /* if we did something, cleanup polynomial (e.g., remove monomials with coefficient 0.0) */
      polynomialdataMergeMonomials(blkmem, polynomialdata, eps, FALSE);
      SCIP_CALL( exprsimplifyRemovePolynomialNullChildren(blkmem, expr) );
   }

   return SCIP_OKAY;
}

/** converts polynomial expressions back into simpler expressions, where possible */
static
SCIP_RETCODE exprsimplifyUnconvertPolynomials(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_EXPR*            expr                /**< expression to convert back */
   )
{
   int i;

   assert(blkmem != NULL);
   assert(expr != NULL);

   for( i = 0; i < expr->nchildren; ++i )
   {
      SCIP_CALL( exprsimplifyUnconvertPolynomials(blkmem, expr->children[i]) );
   }

   if( expr->op != SCIP_EXPR_POLYNOMIAL )
      return SCIP_OKAY;

   SCIP_CALL( exprUnconvertPolynomial(blkmem, &expr->op, &expr->data, expr->nchildren, (void**)expr->children) );

   return SCIP_OKAY;
}

static
SCIP_DECL_HASHGETKEY( exprparseVarTableGetKey )
{  /*lint --e{715}*/
   return (void*)((char*)elem + sizeof(int));
}

/** parses a variable name from a string and creates corresponding expression
 *
 *  Creates a new variable index if variable not seen before, updates varnames and vartable structures.
 */
static
SCIP_RETCODE exprparseReadVariable(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   const char**          str,                /**< pointer to the string to be parsed */
   SCIP_EXPR**           expr,               /**< buffer to store pointer to created expression */
   int*                  nvars,              /**< running number of encountered variables so far */
   int**                 varnames,           /**< pointer to buffer to store new variable names */
   int*                  varnameslength,     /**< pointer to length of the varnames buffer array */
   SCIP_HASHTABLE*       vartable,           /**< hash table for variable names and corresponding expression index */
   SCIP_Real             coefficient,        /**< coefficient to be used when creating the expression */
   const char*           varnameendptr       /**< if a <varname> should be parsed, set this to NULL. Then, str points to the '<'
                                                  else, str should point to the first letter of the varname, and varnameendptr should
                                                  point one char behind the last char of the variable name */
   )
{
   int namelength;
   int varidx;
   char varname[SCIP_MAXSTRLEN];
   void* element;

   assert(blkmem != NULL);
   assert(str != NULL);
   assert(expr != NULL);
   assert(nvars != NULL);
   assert(varnames != NULL);
   assert(vartable != NULL);

   if( varnameendptr == NULL )
   {
      ++*str;
      varnameendptr = *str;
      while( varnameendptr[0] != '>' )
         ++varnameendptr;
   }

   namelength = varnameendptr - *str; /*lint !e712*/
   if( namelength >= SCIP_MAXSTRLEN )
   {
      SCIPerrorMessage("Variable name %.*s is too long for buffer in exprparseReadVariable.\n", namelength, str);
      return SCIP_READERROR;
   }

   memcpy(varname, *str, namelength * sizeof(char));
   varname[namelength] = '\0';

   element = SCIPhashtableRetrieve(vartable, varname);
   if( element != NULL )
   {
      /* variable is old friend */
      assert(strcmp((char*)element + sizeof(int), varname) == 0);

      varidx = *(int*)element;
   }
   else
   {
      /* variable is new */
      varidx = *nvars;

      (*varnameslength) -= (int)(1 + (strlen(varname) + 1) / sizeof(int) + 1);
      if( *varnameslength < 0 )
      {
         SCIPerrorMessage("Buffer in exprparseReadVariable is too short for varaible name %.*s.\n", namelength, str);
         return SCIP_READERROR;
      }

      /* store index of variable and variable name in varnames buffer */
      **varnames = varidx;
      (void) SCIPstrncpy((char*)(*varnames + 1), varname, (int)strlen(varname)+1);

      /* insert variable into hashtable */
      SCIP_CALL( SCIPhashtableInsert(vartable, (void*)*varnames) );

      ++*nvars;
      *varnames += 1 + (strlen(varname) + 1) / sizeof(int) + 1;
   }

   /* create VARIDX expression, put into LINEAR expression if we have coefficient != 1 */
   SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_VARIDX, varidx) );  /*lint !e613*/
   if( coefficient != 1.0 )
   {
      SCIP_CALL( SCIPexprCreateLinear(blkmem, expr, 1, expr, &coefficient, 0.0) );
   }

   /* Move pointer to char behind end of variable */
   *str = varnameendptr + 1;

   /* consprint sometimes prints a variable type identifier which we don't need */
   if( (*str)[0] == '[' && (*str)[2] == ']' &&
       ((*str)[1] == SCIP_VARTYPE_BINARY_CHAR  ||
        (*str)[1] == SCIP_VARTYPE_INTEGER_CHAR ||
        (*str)[1] == SCIP_VARTYPE_IMPLINT_CHAR ||
        (*str)[1] == SCIP_VARTYPE_CONTINUOUS_CHAR ) )
      *str += 3;

   return SCIP_OKAY;
}

/** if str[0] points to an opening parenthesis, this function sets endptr to point to the matching closing bracket in str
 *
 *  Searches for at most length characters.
 */
static
SCIP_RETCODE exprparseFindClosingParenthesis(
   const char*           str,                /**< pointer to the string to be parsed */
   const char**          endptr,             /**< pointer to point to the closing parenthesis */
   int                   length              /**< length of the string to be parsed */
   )
{
   int nopenbrackets;

   assert(str[0] == '(');

   *endptr = str;

   /* find the end of this expression */
   nopenbrackets = 0;
   while( (*endptr - str ) < length && !(nopenbrackets == 1 && *endptr[0] == ')') )
   {
      if( *endptr[0] == '(')
         ++nopenbrackets;
      if( *endptr[0] == ')')
         --nopenbrackets;
      ++*endptr;
   }

   if( *endptr[0] != ')' )
   {
      SCIPerrorMessage("unable to find closing parenthesis in unbalanced expression %.*s\n", length, str);
      return SCIP_READERROR;
   }

   return SCIP_OKAY;
}

/** this function sets endptr to point to the next separating comma in str
 *
 *  That is, for a given string like "x+f(x,y),z", endptr will point to the comma before "z"
 *
 *  Searches for at most length characters.
 */
static
SCIP_RETCODE exprparseFindSeparatingComma(
   const char*           str,                /**< pointer to the string to be parsed */
   const char**          endptr,             /**< pointer to point to the comma */
   int                   length              /**< length of the string to be parsed */
   )
{
   int nopenbrackets;

   *endptr = str;

   /* find a comma without open brackets */
   nopenbrackets = 0;
   while( (*endptr - str ) < length && !(nopenbrackets == 0 && *endptr[0] == ',') )
   {
      if( *endptr[0] == '(')
         ++nopenbrackets;
      if( *endptr[0] == ')')
         --nopenbrackets;
      ++*endptr;
   }

   if( *endptr[0] != ',' )
   {
      SCIPerrorMessage("unable to find separating comma in unbalanced expression %.*s\n", length, str);
      return SCIP_READERROR;
   }

   return SCIP_OKAY;
}

/** parses an expression from a string */
static
SCIP_RETCODE exprParse(
   BMS_BLKMEM*           blkmem,             /**< block memory data structure */
   SCIP_MESSAGEHDLR*     messagehdlr,        /**< message handler */
   SCIP_EXPR**           expr,               /**< buffer to store pointer to created expression */
   const char*           str,                /**< pointer to the string to be parsed */
   int                   length,             /**< length of the string to be parsed */
   const char*           lastchar,           /**< pointer to the last char of str that should be parsed */
   int*                  nvars,              /**< running number of encountered variables so far */
   int**                 varnames,           /**< pointer to buffer to store new variable names */
   int*                  varnameslength,     /**< pointer to length of the varnames buffer array */
   SCIP_HASHTABLE*       vartable,           /**< hash table for variable names and corresponding expression index */
   int                   recursiondepth      /**< current recursion depth */
   )
{   /*lint --e{712,747}*/
   SCIP_EXPR* arg1;
   SCIP_EXPR* arg2;
   const char* subexpptr;
   const char* subexpendptr;
   const char* strstart;
   const char* endptr;
   char* nonconstendptr;
   SCIP_Real number;
   int subexplength;
   int nopenbrackets;

   assert(blkmem != NULL);
   assert(expr != NULL);
   assert(str != NULL);
   assert(lastchar >= str);
   assert(nvars != NULL);
   assert(varnames != NULL);
   assert(vartable != NULL);

   assert(recursiondepth < 100);

   strstart = str; /* might be needed for error message... */

   SCIPdebugMessage("exprParse (%i): parsing %.*s\n", recursiondepth, (int) (lastchar-str + 1), str);

   /* ignore whitespace */
   while( isspace((unsigned char)*str) )
      ++str;

   /* look for a sum or difference not contained in brackets */
   subexpptr = str;
   nopenbrackets = 0;

   /* find the end of this expression
    * a '+' right at the beginning indicates a coefficient, not treated here, or a summation
    */
   while( subexpptr != lastchar && !(nopenbrackets == 0 && (subexpptr[0] == '+' || subexpptr[0] == '-') && subexpptr != str) )
   {
      if( subexpptr[0] == '(')
         ++nopenbrackets;
      if( subexpptr[0] == ')')
         --nopenbrackets;
      ++subexpptr;
   }

   if( subexpptr != lastchar )
   {
      SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str, (int) ((subexpptr - 1) - str + 1), subexpptr - 1, nvars,
            varnames, varnameslength, vartable, recursiondepth + 1) );

      if( subexpptr[0] == '+' )
         ++subexpptr;
      SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, subexpptr , (int) (lastchar - (subexpptr ) + 1), lastchar, nvars,
            varnames, varnameslength, vartable, recursiondepth + 1) );

      /* make new expression from two arguments
       * we always use add, because we leave the operator between the found expressions in the second argument
       * this way, we do not have to worry about ''minus brackets'' in the case of more then two summands:
       *   a - b - c = a + (-b -c)
       */
      SCIP_CALL( SCIPexprAdd(blkmem, expr, 1.0, arg1, 1.0, arg2, 0.0) );

      SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
      SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
      SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );

      return SCIP_OKAY;
   }

   /* check for a bracketed subexpression */
   if( str[0] == '(' )
   {
      nopenbrackets = 0;

      subexplength = -1;    /* we do not want the closing bracket in the string */
      subexpptr = str + 1;  /* leave out opening bracket */

      /* find the end of this expression */
      while( subexplength < length && !(nopenbrackets == 1 && str[0] == ')') )
      {
         if( str[0] == '(' )
            ++nopenbrackets;
         if( str[0] == ')' )
            --nopenbrackets;
         ++str;
         ++subexplength;
      }
      subexpendptr = str - 1; /* leave out closing bracket */

      SCIP_CALL( exprParse(blkmem, messagehdlr, expr, subexpptr, subexplength, subexpendptr, nvars, varnames,
            varnameslength, vartable, recursiondepth + 1) );
      ++str;
   }
   else if( isdigit((unsigned char)str[0]) || ((str[0] == '-' || str[0] == '+')
         && (isdigit((unsigned char)str[1]) || str[1] == ' ')) )
   {
      /* check if there is a lonely minus coming, indicating a -1.0 */
      if( str[0] == '-'  && str[1] == ' ' )
      {
         number = -1.0;
         nonconstendptr = (char*) str + 1;
      }
      /* check if there is a number coming */
      else if( !SCIPstrToRealValue(str, &number, &nonconstendptr) )
      {
         SCIPerrorMessage("error parsing number from <%s>\n", str);
         return SCIP_READERROR;
      }
      str = nonconstendptr;

      /* ignore whitespace */
      while( isspace((unsigned char)*str) && str != lastchar )
         ++str;

      if( str[0] != '*' && str[0] != '/' && str[0] != '+' && str[0] != '-' && str[0] != '^' )
      {
         if( str < lastchar )
         {
            SCIP_CALL( exprParse(blkmem, messagehdlr, expr, str, (int)(lastchar - str) + 1, lastchar, nvars, varnames,
                  varnameslength, vartable, recursiondepth + 1) );
            SCIP_CALL( SCIPexprMulConstant(blkmem, expr, *expr, number) );
         }
         else
         {
            SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_CONST, number) );
         }
         str = lastchar + 1;
      }
      else
      {
         SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_CONST, number) );
      }
   }
   else if( str[0] == '<' )
   {
      /* check if expressions begins with a variable */
      SCIP_CALL( exprparseReadVariable(blkmem, &str, expr, nvars, varnames, varnameslength, vartable, 1.0, NULL) );
   }
   /* four character operators */
   else if( strncmp(str, "sqrt", 4) == 0 )
   {
      str += 4;
      SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
      SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, endptr - str - 1, endptr -1, nvars, varnames,
            varnameslength, vartable, recursiondepth + 1) );
      str = endptr + 1;

      SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_SQRT, arg1) );
   }
   /* three character operators with 1 argument */
   else if(
      strncmp(str, "abs", 3) == 0 ||
      strncmp(str, "cos", 3) == 0 ||
      strncmp(str, "exp", 3) == 0 ||
      strncmp(str, "log", 3) == 0 ||
      strncmp(str, "sin", 3) == 0 ||
      strncmp(str, "sqr", 3) == 0 ||
      strncmp(str, "tan", 3) == 0 )
   {
      const char* opname = str;

      str += 3;
      SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
      SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, endptr - str - 1, endptr -1, nvars, varnames,
            varnameslength, vartable, recursiondepth + 1) );
      str = endptr + 1;

      if( strncmp(opname, "abs", 3) == 0 )
      {
         SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_ABS, arg1) );
      }
      else if( strncmp(opname, "cos", 3) == 0 )
      {
         SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_COS, arg1) );
      }
      else if( strncmp(opname, "exp", 3) == 0 )
      {
         SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_EXP, arg1) );
      }
      else if( strncmp(opname, "log", 3) == 0 )
      {
         SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_LOG, arg1) );
      }
      else if( strncmp(opname, "sin", 3) == 0 )
      {
         SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_SIN, arg1) );
      }
      else if( strncmp(opname, "sqr", 3) == 0 )
      {
         SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_SQUARE, arg1) );
      }
      else
      {
         assert(strncmp(opname, "tan", 3) == 0);
         SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_TAN, arg1) );
      }
   }
   /* three character operators with 2 arguments */
   else if(
      strncmp(str, "max", 3) == 0 ||
      strncmp(str, "min", 3) == 0 )
   {
      /* we have a string of the form "min(...,...)" or "max(...,...)"
       * first find the closing parenthesis, then the comma
       */
      const char* comma;
      SCIP_EXPROP op;

      op = (str[1] == 'a' ? SCIP_EXPR_MAX : SCIP_EXPR_MIN);

      str += 3;
      SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );

      SCIP_CALL( exprparseFindSeparatingComma(str+1, &comma, endptr - str - 1) );

      /* parse first argument [str+1..comma-1] */
      SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, comma - str - 1, comma - 1, nvars, varnames,
            varnameslength, vartable, recursiondepth + 1) );

      /* parse second argument [comma+1..endptr] */
      ++comma;
      while( comma < endptr && *comma == ' ' )
         ++comma;

      SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, comma, endptr - comma, endptr - 1, nvars, varnames,
            varnameslength, vartable, recursiondepth + 1) );

      SCIP_CALL( SCIPexprCreate(blkmem, expr, op, arg1, arg2) );

      str = endptr + 1;
   }
   else if( strncmp(str, "power", 5) == 0 )
   {
      /* we have a string of the form "power(...,integer)" (thus, intpower)
       * first find the closing parenthesis, then the comma
       */
      const char* comma;
      int exponent;

      str += 5;
      SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );

      SCIP_CALL( exprparseFindSeparatingComma(str+1, &comma, endptr - str - 1) );

      /* parse first argument [str+1..comma-1] */
      SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, comma - str - 1, comma - 1, nvars, varnames,
            varnameslength, vartable, recursiondepth + 1) );

      ++comma;
      /* parse second argument [comma, endptr-1]: it needs to be an integer */
      while( comma < endptr && *comma == ' ' )
         ++comma;
      if( !isdigit((unsigned char)comma[0]) && !((comma[0] == '-' || comma[0] == '+') && isdigit((unsigned char)comma[1])) )
      {
         SCIPerrorMessage("error parsing integer exponent from <%s>\n", comma);
      }
      if( !SCIPstrToIntValue(comma, &exponent, &nonconstendptr) )
      {
         SCIPerrorMessage("error parsing integer from <%s>\n", comma);
         return SCIP_READERROR;
      }

      SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_INTPOWER, arg1, exponent) );

      str = endptr + 1;
   }
   else if( strncmp(str, "realpower", 9) == 0 || strncmp(str, "signpower", 9) == 0 )
   {
      /* we have a string of the form "realpower(...,double)" or "signpower(...,double)"
       * first find the closing parenthesis, then the comma
       */
      const char* opname = str;
      const char* comma;

      str += 9;
      SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );

      SCIP_CALL( exprparseFindSeparatingComma(str+1, &comma, endptr - str - 1) );

      /* parse first argument [str+1..comma-1] */
      SCIP_CALL( exprParse(blkmem, messagehdlr, &arg1, str + 1, comma - str - 1, comma - 1, nvars, varnames,
            varnameslength, vartable, recursiondepth + 1) );

      ++comma;
      /* parse second argument [comma, endptr-1]: it needs to be an number */
      while( comma < endptr && *comma == ' ' )
         ++comma;
      if( !isdigit((unsigned char)comma[0]) && !((comma[0] == '-' || comma[0] == '+') && isdigit((unsigned char)comma[1])) )
      {
         SCIPerrorMessage("error parsing number exponent from <%s>\n", comma);
      }
      if( !SCIPstrToRealValue(comma, &number, &nonconstendptr) )
      {
         SCIPerrorMessage("error parsing number from <%s>\n", comma);
         return SCIP_READERROR;
      }

      if( strncmp(opname, "realpower", 9) == 0 )
      {
         SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_REALPOWER, arg1, number) );
      }
      else
      {
         assert(strncmp(opname, "signpower", 9) == 0);
         SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_SIGNPOWER, arg1, number) );
      }

      str = endptr + 1;
   }
   else if( isalpha(*str) || *str == '_' || *str == '#' )
   {
      /* check for a variable, that was not recognized earlier because somebody omitted the '<' and '>' we need for
       * SCIPparseVarName, making everyones life harder;
       * we allow only variable names starting with a character or underscore here
       */
      const char* varnamestartptr = str;

      /* allow only variable names containing characters, digits, hash marks, and underscores here */
      while( isalnum(str[0]) || str[0] == '_' || str[0] == '#' )
         ++str;

      SCIP_CALL( exprparseReadVariable(blkmem, &varnamestartptr, expr, nvars, varnames, varnameslength,
            vartable, 1.0, str) );
   }
   else
   {
      SCIPerrorMessage("parsing of invalid expression %.*s.\n", (int) (lastchar - str + 1), str);
      return SCIP_READERROR;
   }

   /* if we are one char behind lastchar, we are done */
   if( str == lastchar + 1)
   {
      SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
      SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
      SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );

      return SCIP_OKAY;
   }

   /* check if we are still in bounds */
   if( str > lastchar + 1)
   {
      SCIPerrorMessage("error finding first expression in \"%.*s\" took us outside of given subexpression length\n", length, strstart);
      return SCIP_READERROR;
   }

   /* ignore whitespace */
   while( isspace((unsigned char)*str) && str != lastchar + 1 )
      ++str;

   /* maybe now we're done? */
   if( str >= lastchar + 1)
   {
      SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
      SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
      SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );

      return SCIP_OKAY;
   }

   if( str[0] == '^' )
   {
      /* a '^' behind the found expression indicates a power */
      SCIP_Real constant;

      arg1 = *expr;
      ++str;
      while( isspace((unsigned char)*str) && str != lastchar + 1 )
         ++str;

      if( isdigit((unsigned char)str[0]) || ((str[0] == '-' || str[0] == '+') && isdigit((unsigned char)str[1])) )
      {
         /* there is a number coming */
         if( !SCIPstrToRealValue(str, &number, &nonconstendptr) )
         {
            SCIPerrorMessage("error parsing number from <%s>\n", str);
            return SCIP_READERROR;
         }

         SCIP_CALL( SCIPexprCreate(blkmem, &arg2, SCIP_EXPR_CONST, number) );
         str = nonconstendptr;

         constant = SCIPexprGetOpReal(arg2);
         SCIPexprFreeDeep(blkmem, &arg2);

         /* expr^number is intpower or realpower */
         if( EPSISINT(constant, 0.0) ) /*lint !e835*/
         {
            SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_INTPOWER, arg1, (int)constant) );
         }
         else
         {
            SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_REALPOWER, arg1, constant) );
         }
      }
      else if( str[0] == '(' )
      {
         /* we use exprParse to evaluate the exponent */

         SCIP_CALL( exprparseFindClosingParenthesis(str, &endptr, length) );
         SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, str + 1, endptr - str - 1, endptr -1, nvars, varnames,
               varnameslength, vartable, recursiondepth + 1) );

         if( SCIPexprGetOperator(arg2) != SCIP_EXPR_CONST )
         {
            /* reformulate arg1^arg2 as exp(arg2 * log(arg1)) */
            SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_LOG, arg1) );
            SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_MUL, *expr, arg2) );
            SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_EXP, *expr) );
         }
         else
         {
            /* expr^number is intpower or realpower */
            constant = SCIPexprGetOpReal(arg2);
            SCIPexprFreeDeep(blkmem, &arg2);
            if( EPSISINT(constant, 0.0) ) /*lint !e835*/
            {
               SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_INTPOWER, arg1, (int)constant) );
            }
            else
            {
               SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_REALPOWER, arg1, constant) );
            }
         }
         str = endptr + 1;
      }
      else
      {
         SCIPerrorMessage("unexpected string following ^ in  %.*s\n", length, str);
         return SCIP_READERROR;
      }

      /* ignore whitespace */
      while( isspace((unsigned char)*str) && str != lastchar + 1 )
         ++str;
   }

   /* check for a two argument operator that is not a multiplication */
   if( str <= lastchar && (str[0] == '+' || str[0] == '-' || str[0] == '/') )
   {
      char op;

      op = str[0];
      arg1 = *expr;

      /* step forward over the operator to go to the beginning of the second argument */
      ++str;

      SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, str, (int) (lastchar - str + 1), lastchar, nvars, varnames,
            varnameslength, vartable, recursiondepth + 1) );
      str = lastchar + 1;

      /* make new expression from two arguments */
      if( op == '+')
      {
         SCIP_CALL( SCIPexprAdd(blkmem, expr, 1.0, arg1, 1.0, arg2, 0.0) );
      }
      else if( op == '-')
      {
         SCIP_CALL( SCIPexprAdd(blkmem, expr, 1.0, arg1, -1.0, arg2, 0.0) );
      }
      else if( op == '*' )
      {
         if( SCIPexprGetOperator(arg1) == SCIP_EXPR_CONST )
         {
            SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg2, SCIPexprGetOpReal(arg1)) );
         }
         else if( SCIPexprGetOperator(arg2) == SCIP_EXPR_CONST )
         {
            SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg1, SCIPexprGetOpReal(arg2)) );
         }
         else
         {
            SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_MUL, arg1, arg2) );
         }
      }
      else
      {
         assert(op == '/');

         if( SCIPexprGetOperator(arg2) == SCIP_EXPR_CONST )
         {
            SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg1, 1.0 / SCIPexprGetOpReal(arg2)) );
            SCIPexprFreeShallow(blkmem, &arg2);
         }
         else
         {
            SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_DIV, arg1, arg2) );
         }
      }
   }

   /* ignore whitespace */
   while( isspace((unsigned char)*str) )
      ++str;

   /* we are either done or we have a multiplication? */
   if( str >= lastchar + 1 )
   {
      SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
      SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
      SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );

      return SCIP_OKAY;
   }

   /* if there is a part of the string left to be parsed, we assume that this as a multiplication */
   arg1 = *expr;

   /* stepping over multiplication operator if needed */
   if( str[0] == '*' )
   {
      ++str;
   }
   else if( str[0] != '(' )
   {
      SCIPdebugMessage("No operator found, assuming a multiplication before %.*s\n", (int) (lastchar - str + 1), str);
   }

   SCIP_CALL( exprParse(blkmem, messagehdlr, &arg2, str, (int) (lastchar - str + 1), lastchar, nvars, varnames,
         varnameslength, vartable, recursiondepth + 1) );

   if( SCIPexprGetOperator(arg1) == SCIP_EXPR_CONST )
   {
      SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg2, SCIPexprGetOpReal(arg1)) );
      SCIPexprFreeDeep(blkmem, &arg1);
   }
   else if( SCIPexprGetOperator(arg2) == SCIP_EXPR_CONST )
   {
      SCIP_CALL( SCIPexprMulConstant(blkmem, expr, arg1, SCIPexprGetOpReal(arg2)) );
      SCIPexprFreeDeep(blkmem, &arg2);
   }
   else
   {
      SCIP_CALL( SCIPexprCreate(blkmem, expr, SCIP_EXPR_MUL, arg1, arg2) );
   }

   SCIPdebugMessage("exprParse (%i): returns expression ", recursiondepth);
   SCIPdebug( SCIPexprPrint(*expr, messagehdlr, NULL, NULL, NULL, NULL) );
   SCIPdebug( SCIPmessagePrintInfo(messagehdlr, "\n") );

   return SCIP_OKAY;
}

/**@} */

/**@name Expression methods */
/**@{ */

/* In debug mode, the following methods are implemented as function calls to ensure
 * type validity.
 * In optimized mode, the methods are implemented as defines to improve performance.
 * However, we want to have them in the library anyways, so we have to undef the defines.
 */

#undef SCIPexprGetOperator
#undef SCIPexprGetNChildren
#undef SCIPexprGetChildren
#undef SCIPexprGetOpIndex
#undef SCIPexprGetOpReal
#undef SCIPexprGetOpData
#undef SCIPexprGetRealPowerExponent
#undef SCIPexprGetIntPowerExponent
#undef SCIPexprGetSignPowerExponent
#undef SCIPexprGetLinearCoefs
#undef SCIPexprGetLinearConstant
#undef SCIPexprGetQuadElements
#undef SCIPexprGetQuadConstant
#undef SCIPexprGetQuadLinearCoefs
#undef SCIPexprGetNQuadElements
#undef SCIPexprGetMonomials
#undef SCIPexprGetNMonomials
#undef SCIPexprGetPolynomialConstant
#undef SCIPexprGetMonomialCoef
#undef SCIPexprGetMonomialNFactors
#undef SCIPexprGetMonomialChildIndices
#undef SCIPexprGetMonomialExponents
#undef SCIPexprGetUserData
#undef SCIPexprHasUserEstimator
#undef SCIPexprGetUserEvalCapability

/** gives operator of expression */
SCIP_EXPROP SCIPexprGetOperator(
   SCIP_EXPR*            expr                /**< expression */
   )
{
   assert(expr != NULL);

   return expr->op;
}

/** gives number of children of an expression */
int SCIPexprGetNChildren(
   SCIP_EXPR*            expr                /**< expression */
   )
{
   assert(expr != NULL);

   return expr->nchildren;
}

/** gives pointer to array with children of an expression */
SCIP_EXPR** SCIPexprGetChildren(
   SCIP_EXPR*            expr                /**< expression */
   )
{
   assert(expr != NULL);

   return expr->children;
}

/** gives index belonging to a SCIP_EXPR_VARIDX or SCIP_EXPR_PARAM operand */
int SCIPexprGetOpIndex(
   SCIP_EXPR*            expr                /**< expression */
   )
{
   assert(expr != NULL);
   assert(expr->op == SCIP_EXPR_VARIDX || expr->op == SCIP_EXPR_PARAM);

   return expr->data.intval;
}

/** gives real belonging to a SCIP_EXPR_CONST operand */
SCIP_Real SCIPexprGetOpReal(
   SCIP_EXPR*            expr                /**< expression */
   )
{
   assert(expr != NULL);
   assert(expr->op == SCIP_EXPR_CONST);

   return expr->data.dbl;
}

/** gives void* belonging to a complex operand */
void* SCIPexprGetOpData(
   SCIP_EXPR*            expr                /**< expression */
   )
{
   assert(expr != NULL);
   assert(expr->op >= SCIP_EXPR_SUM); /* only complex operands store their data as void* */

   return expr->data.data;
}

/** gives exponent belonging to a SCIP_EXPR_REALPOWER expression */
SCIP_Real SCIPexprGetRealPowerExponent(
   SCIP_EXPR*            expr                /**< expression */
   )
{
   assert(expr != NULL);
   assert(expr->op == SCIP_EXPR_REALPOWER);

   return expr->data.dbl;
}

/** gives exponent belonging to a SCIP_EXPR_INTPOWER expression */
int SCIPexprGetIntPowerExponent(
   SCIP_EXPR*            expr                /**< expression */
   )
{
   assert(expr != NULL);
   assert(expr->op == SCIP_EXPR_INTPOWER);

   return expr->data.intval;
}

/** gives exponent belonging to a SCIP_EXPR_SIGNPOWER expression */
SCIP_Real SCIPexprGetSignPowerExponent(
   SCIP_EXPR*            expr                /**< expression */
   )
{
   assert(expr != NULL);
   assert(expr->op == SCIP_EXPR_SIGNPOWER);

   return expr->data.dbl;
}

/** gives linear coefficients belonging to a SCIP_EXPR_LINEAR expression */
SCIP_Real* SCIPexprGetLinearCoefs(
   SCIP_EXPR*            expr                /**< expression */
   )
{
   assert(expr != NULL);
   assert(expr->op == SCIP_EXPR_LINEAR);
   assert(expr->data.data != NULL);

   /* the coefficients are stored in the first nchildren elements of the array stored as expression data */
   return (SCIP_Real*)expr->data.data;
}

/** gives constant belonging to a SCIP_EXPR_LINEAR expression */
SCIP_Real SCIPexprGetLinearConstant(
   SCIP_EXPR*            expr                /**< expression */
   )
{
   assert(expr != NULL);
   assert(expr->op == SCIP_EXPR_LINEAR);
   assert(expr->data.data != NULL);

   /* the constant is stored in the nchildren's element of the array stored as expression data */
   return ((SCIP_Real*)expr->data.data)[expr->nchildren];
}

/** gives quadratic elements belonging to a SCIP_EXPR_QUADRATIC expression */
SCIP_QUADELEM* SCIPexprGetQuadElements(
   SCIP_EXPR*            expr                /**< quadratic expression */
   )
{
   assert(expr != NULL);
   assert(expr->op == SCIP_EXPR_QUADRATIC);
   assert(expr->data.data != NULL);

   return ((SCIP_EXPRDATA_QUADRATIC*)expr->data.data)->quadelems;
}

/** gives constant belonging to a SCIP_EXPR_QUADRATIC expression */
SCIP_Real SCIPexprGetQuadConstant(
   SCIP_EXPR*            expr                /**< quadratic expression */
   )
{
   assert(expr != NULL);
   assert(expr->op == SCIP_EXPR_QUADRATIC);
   assert(expr->data.data != NULL);

   return ((SCIP_EXPRDATA_QUADRATIC*)expr->data.data)->constant;
}

/** gives linear coefficients belonging to a SCIP_EXPR_QUADRATIC expression
 * can be NULL if all coefficients are 0.0 */
SCIP_Real* SCIPexprGetQuadLinearCoefs(
   SCIP_EXPR*            expr                /**< quadratic expression */
   )
{
   assert(expr != NULL);
   assert(expr->op == SCIP_EXPR_QUADRATIC);
   assert(expr->data.data != NULL);

   return ((SCIP_EXPRDATA_QUADRATIC*)expr->data.data)->lincoefs;
}

/** gives number of quadratic elements belonging to a SCIP_EXPR_QUADRATIC expression */
int SCIPexprGetNQuadElements(
   SCIP_EXPR*            expr                /**< quadratic expression */
   )
{
   assert(expr != NULL);
   assert(expr->op == SCIP_EXPR_QUADRATIC);
   assert(expr->data.data != NULL);

   return ((SCIP_EXPRDATA_QUADRATIC*)expr->data.data)->nquadelems;
}

/** gives the monomials belonging to a SCIP_EXPR_POLYNOMIAL expression */
SCIP_EXPRDATA_MONOMIAL** SCIPexprGetMonomials(
   SCIP_EXPR*            expr                /**< expression */
   )
{
   assert(expr != NULL);
   assert(expr->op == SCIP_EXPR_POLYNOMIAL);
   assert(expr->data.data != NULL);