1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the program and library             */
4    	/*         SCIP --- Solving Constraint Integer Programs                      */
5    	/*                                                                           */
6    	/*  Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB)                      */
7    	/*                                                                           */
8    	/*  Licensed under the Apache License, Version 2.0 (the "License");          */
9    	/*  you may not use this file except in compliance with the License.         */
10   	/*  You may obtain a copy of the License at                                  */
11   	/*                                                                           */
12   	/*      http://www.apache.org/licenses/LICENSE-2.0                           */
13   	/*                                                                           */
14   	/*  Unless required by applicable law or agreed to in writing, software      */
15   	/*  distributed under the License is distributed on an "AS IS" BASIS,        */
16   	/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17   	/*  See the License for the specific language governing permissions and      */
18   	/*  limitations under the License.                                           */
19   	/*                                                                           */
20   	/*  You should have received a copy of the Apache-2.0 license                */
21   	/*  along with SCIP; see the file LICENSE. If not visit scipopt.org.         */
22   	/*                                                                           */
23   	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24   	
25   	/**@file   expr_pow.c
26   	 * @ingroup DEFPLUGINS_EXPR
27   	 * @brief  power expression handler
28   	 * @author Benjamin Mueller
29   	 * @author Ksenia Bestuzheva
30   	 *
31   	 * @todo signpower for exponent < 1 ?
32   	 */
33   	
34   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
35   	
36   	/*lint --e{835}*/
37   	/*lint -e777*/
38   	
39   	#include <string.h>
40   	
41   	#include "scip/expr_pow.h"
42   	#include "scip/pub_expr.h"
43   	#include "scip/expr_value.h"
44   	#include "scip/expr_product.h"
45   	#include "scip/expr_sum.h"
46   	#include "scip/expr_exp.h"
47   	#include "scip/expr_abs.h"
48   	#include "symmetry/struct_symmetry.h"
49   	
50   	#define POWEXPRHDLR_NAME         "pow"
51   	#define POWEXPRHDLR_DESC         "power expression"
52   	#define POWEXPRHDLR_PRECEDENCE   55000
53   	#define POWEXPRHDLR_HASHKEY      SCIPcalcFibHash(21163.0)
54   	
55   	#define SIGNPOWEXPRHDLR_NAME       "signpower"
56   	#define SIGNPOWEXPRHDLR_DESC       "signed power expression"
57   	#define SIGNPOWEXPRHDLR_PRECEDENCE 56000
58   	#define SIGNPOWEXPRHDLR_HASHKEY    SCIPcalcFibHash(21163.1)
59   	
60   	#define INITLPMAXPOWVAL          1e+06 /**< maximal allowed absolute value of power expression at bound,
61   	                                        *   used for adjusting bounds in the convex case in initestimates */
62   	
63   	/*
64   	 * Data structures
65   	 */
66   	
67   	/** sign of a value (-1 or +1)
68   	 *
69   	 * 0.0 has sign +1 here (shouldn't matter, though)
70   	 */
71   	#define SIGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
72   	
73   	#define SIGNPOW_ROOTS_KNOWN 10   /**< up to which (integer) exponents precomputed roots have been stored */
74   	
75   	/** The positive root of the polynomial (n-1) y^n + n y^(n-1) - 1 is needed in separation.
76   	 *  Here we store these roots for small integer values of n.
77   	 */
78   	static
79   	SCIP_Real signpow_roots[SIGNPOW_ROOTS_KNOWN+1] = {
80   	   -1.0,                     /* no root for n=0 */
81   	   -1.0,                     /* no root for n=1 */
82   	   0.41421356237309504880,   /* root for n=2 (-1+sqrt(2)) */
83   	   0.5,                      /* root for n=3 */
84   	   0.56042566045031785945,   /* root for n=4 */
85   	   0.60582958618826802099,   /* root for n=5 */
86   	   0.64146546982884663257,   /* root for n=6 */
87   	   0.67033204760309682774,   /* root for n=7 */
88   	   0.69428385661425826738,   /* root for n=8 */
89   	   0.71453772716733489700,   /* root for n=9 */
90   	   0.73192937842370733350    /* root for n=10 */
91   	};
92   	
93   	/** expression handler data */
94   	struct SCIP_ExprhdlrData
95   	{
96   	   SCIP_Real             minzerodistance;    /**< minimal distance from zero to enforce for child in bound tightening */
97   	   int                   expandmaxexponent;  /**< maximal exponent when to expand power of sum in simplify */
98   	   SCIP_Bool             distribfracexponent;/**< whether a fractional exponent is distributed onto factors on power of product */
99   	
100  	   SCIP_Bool             warnedonpole;       /**< whether we warned on enforcing a minimal distance from zero for child */
101  	};
102  	
103  	/** expression data */
104  	struct SCIP_ExprData
105  	{
106  	   SCIP_Real             exponent;           /**< exponent */
107  	   SCIP_Real             root;               /**< positive root of (n-1) y^n + n y^(n-1) - 1, or
108  	                                                  SCIP_INVALID if not computed yet */
109  	};
110  	
111  	/*
112  	 * Local methods
113  	 */
114  	
115  	/** computes positive root of the polynomial (n-1) y^n + n y^(n-1) - 1 for n > 1 */
116  	static
117  	SCIP_RETCODE computeSignpowerRoot(
118  	   SCIP*                 scip,               /**< SCIP data structure */
119  	   SCIP_Real*            root,               /**< buffer where to store computed root */
120  	   SCIP_Real             exponent            /**< exponent n */
121  	   )
122  	{
123  	   SCIP_Real polyval;
124  	   SCIP_Real gradval;
125  	   int iter;
126  	
127  	   assert(scip != NULL);
128  	   assert(exponent > 1.0);
129  	   assert(root != NULL);
130  	
131  	   /* lookup for popular integer exponent */
132  	   if( SCIPisIntegral(scip, exponent) && exponent-0.5 < SIGNPOW_ROOTS_KNOWN )
133  	   {
134  	      *root = signpow_roots[(int)SCIPfloor(scip, exponent+0.5)];
135  	      return SCIP_OKAY;
136  	   }
137  	
138  	   /* lookup for weymouth exponent */
139  	   if( SCIPisEQ(scip, exponent, 1.852) )
140  	   {
141  	      *root = 0.39821689389382575186;
142  	      return SCIP_OKAY;
143  	   }
144  	
145  	   /* search for a positive root of (n-1) y^n + n y^(n-1) - 1
146  	    * use the closest precomputed root as starting value
147  	    */
148  	   if( exponent >= SIGNPOW_ROOTS_KNOWN )
149  	      *root = signpow_roots[SIGNPOW_ROOTS_KNOWN];
150  	   else if( exponent <= 2.0 )
151  	      *root = signpow_roots[2];
152  	   else
153  	      *root = signpow_roots[(int)SCIPfloor(scip, exponent)];
154  	
155  	   for( iter = 0; iter < 1000; ++iter )
156  	   {
157  	      polyval = (exponent - 1.0) * pow(*root, exponent) + exponent * pow(*root, exponent - 1.0) - 1.0;
158  	      if( fabs(polyval) < 1e-12 && SCIPisZero(scip, polyval) )
159  	         break;
160  	
161  	      /* gradient of (n-1) y^n + n y^(n-1) - 1 is n(n-1)y^(n-1) + n(n-1)y^(n-2) */
162  	      gradval = (exponent - 1.0) * exponent * (pow(*root, exponent - 1.0) + pow(*root, exponent - 2.0));
163  	      if( SCIPisZero(scip, gradval) )
164  	         break;
165  	
166  	      /* update root by adding -polyval/gradval (Newton's method) */
167  	      *root -= polyval / gradval;
168  	      if( *root < 0.0 )
169  	         *root = 0.0;
170  	   }
171  	
172  	   if( !SCIPisZero(scip, polyval) )
173  	   {
174  	      SCIPerrorMessage("failed to compute root for exponent %g\n", exponent);
175  	      return SCIP_ERROR;
176  	   }
177  	   SCIPdebugMsg(scip, "root for %g is %.20g, certainty = %g\n", exponent, *root, polyval);
178  	   /* @todo cache root value for other expressions (an exponent seldom comes alone)?? (they are actually really fast to compute...) */
179  	
180  	   return SCIP_OKAY;
181  	}
182  	
183  	/** computes negative root of the polynomial (n-1) y^n - n y^(n-1) + 1 for n < -1 */
184  	static
185  	SCIP_RETCODE computeHyperbolaRoot(
186  	   SCIP*                 scip,               /**< SCIP data structure */
187  	   SCIP_Real*            root,               /**< buffer where to store computed root */
188  	   SCIP_Real             exponent            /**< exponent n */
189  	   )
190  	{
191  	   SCIP_Real polyval;
192  	   SCIP_Real gradval;
193  	   int iter;
194  	
195  	   assert(scip != NULL);
196  	   assert(exponent < -1.0);
197  	   assert(root != NULL);
198  	
199  	   *root = -2.0;  /* that's the solution for n=-2 */
200  	
201  	   for( iter = 0; iter < 1000; ++iter )
202  	   {
203  	      polyval = (exponent - 1.0) * pow(*root, exponent) - exponent * pow(*root, exponent - 1.0) + 1.0;
204  	      if( fabs(polyval) < 1e-12 && SCIPisZero(scip, polyval) )
205  	         break;
206  	
207  	      /* gradient of (n-1) y^n - n y^(n-1) + 1 is n(n-1)y^(n-1) - n(n-1)y^(n-2) */
208  	      gradval = (exponent - 1.0) * exponent * (pow(*root, exponent - 1.0) - pow(*root, exponent - 2.0));
209  	      if( SCIPisZero(scip, gradval) )
210  	         break;
211  	
212  	      /* update root by adding -polyval/gradval (Newton's method) */
213  	      *root -= polyval / gradval;
214  	      if( *root >= 0.0 )
215  	         *root = -1;
216  	   }
217  	
218  	   if( !SCIPisZero(scip, polyval) )
219  	   {
220  	      SCIPerrorMessage("failed to compute root for exponent %g\n", exponent);
221  	      return SCIP_ERROR;
222  	   }
223  	   SCIPdebugMsg(scip, "root for %g is %.20g, certainty = %g\n", exponent, *root, polyval);
224  	   /* @todo cache root value for other expressions (an exponent seldom comes alone)?? (they are actually really fast to compute...) */
225  	
226  	   return SCIP_OKAY;
227  	}
228  	
229  	/** creates expression data */
230  	static
231  	SCIP_RETCODE createData(
232  	   SCIP*                 scip,               /**< SCIP data structure */
233  	   SCIP_EXPRDATA**       exprdata,           /**< pointer where to store expression data */
234  	   SCIP_Real             exponent            /**< exponent of the power expression */
235  	   )
236  	{
237  	   assert(exprdata != NULL);
238  	
239  	   SCIP_CALL( SCIPallocBlockMemory(scip, exprdata) );
240  	
241  	   (*exprdata)->exponent = exponent;
242  	   (*exprdata)->root = SCIP_INVALID;
243  	
244  	   return SCIP_OKAY;
245  	}
246  	
247  	/** computes a tangent at a reference point by linearization
248  	 *
249  	 * for a normal power, linearization in xref is xref^exponent + exponent * xref^(exponent-1) (x - xref)
250  	 * = (1-exponent) * xref^exponent + exponent * xref^(exponent-1) * x
251  	 *
252  	 * for a signpower, linearization is the same if xref is positive
253  	 * for xref negative it is -(-xref)^exponent + exponent * (-xref)^(exponent-1) (x-xref)
254  	 * = (1-exponent) * (-xref)^(exponent-1) * xref + exponent * (-xref)^(exponent-1) * x
255  	 */
256  	static
257  	void computeTangent(
258  	   SCIP*                 scip,               /**< SCIP data structure */
259  	   SCIP_Bool             signpower,          /**< are we signpower or normal power */
260  	   SCIP_Real             exponent,           /**< exponent */
261  	   SCIP_Real             xref,               /**< reference point where to linearize */
262  	   SCIP_Real*            constant,           /**< buffer to store constant term of secant */
263  	   SCIP_Real*            slope,              /**< buffer to store slope of secant */
264  	   SCIP_Bool*            success             /**< buffer to store whether secant could be computed */
265  	   )
266  	{
267  	   SCIP_Real xrefpow;
268  	
269  	   assert(scip != NULL);
270  	   assert(constant != NULL);
271  	   assert(slope != NULL);
272  	   assert(success != NULL);
273  	   assert(xref != 0.0 || exponent > 0.0);
274  	   /* non-integral exponent -> reference point must be >= 0 or we do signpower */
275  	   assert(EPSISINT(exponent, 0.0) || signpower || !SCIPisNegative(scip, xref));
276  	
277  	   /* TODO power is not differentiable at 0.0 for exponent < 0
278  	    * should we forbid here that xref > 0, do something smart here, or just return success=FALSE?
279  	    */
280  	   /* assert(exponent >= 1.0 || xref > 0.0); */
281  	
282  	   if( !EPSISINT(exponent, 0.0) && !signpower && xref < 0.0 )
283  	      xref = 0.0;
284  	
285  	   xrefpow = pow(signpower ? REALABS(xref) : xref, exponent - 1.0);
286  	
287  	   /* if huge xref and/or exponent too large, then pow may overflow */
288  	   if( !SCIPisFinite(xrefpow) )
289  	   {
290  	      *success = FALSE;
291  	      return;
292  	   }
293  	
294  	   *constant = (1.0 - exponent) * xrefpow * xref;
295  	   *slope = exponent * xrefpow;
296  	   *success = TRUE;
297  	}
298  	
299  	/** computes a secant between lower and upper bound
300  	 *
301  	 * secant is xlb^exponent + (xub^exponent - xlb^exponent) / (xub - xlb) * (x - xlb)
302  	 * = xlb^exponent - slope * xlb + slope * x  with slope = (xub^exponent - xlb^exponent) / (xub - xlb)
303  	 * same if signpower
304  	 */
305  	static
306  	void computeSecant(
307  	   SCIP*                 scip,               /**< SCIP data structure */
308  	   SCIP_Bool             signpower,          /**< are we signpower or normal power */
309  	   SCIP_Real             exponent,           /**< exponent */
310  	   SCIP_Real             xlb,                /**< lower bound on x */
311  	   SCIP_Real             xub,                /**< upper bound on x */
312  	   SCIP_Real*            constant,           /**< buffer to store constant term of secant */
313  	   SCIP_Real*            slope,              /**< buffer to store slope of secant */
314  	   SCIP_Bool*            success             /**< buffer to store whether secant could be computed */
315  	   )
316  	{
317  	   assert(scip != NULL);
318  	   assert(constant != NULL);
319  	   assert(slope != NULL);
320  	   assert(success != NULL);
321  	   assert(xlb >= 0.0 || EPSISINT(exponent, 0.0) || signpower);
322  	   assert(xub >= 0.0 || EPSISINT(exponent, 0.0) || signpower);
323  	   assert(exponent != 1.0);
324  	
325  	   *success = FALSE;
326  	
327  	   /* infinite bounds will not work */
328  	   if( SCIPisInfinity(scip, -xlb) || SCIPisInfinity(scip, xub) )
329  	      return;
330  	
331  	   /* first handle some special cases */
332  	   if( xlb == xub )
333  	   {
334  	      /* usually taken care of in separatePointPow already, but we might be called with different bounds here,
335  	       * e.g., when handling odd or signed power
336  	       */
337  	      *slope = 0.0;
338  	      *constant = pow(xlb, exponent);
339  	   }
340  	   else if( EPSISINT(exponent / 2.0, 0.0) && !signpower && xub > 0.1 && SCIPisFeasEQ(scip, xlb, -xub) )
341  	   {
342  	      /* for normal power with even exponents with xlb ~ -xub the slope would be very close to 0
343  	       * since xub^n - xlb^n is prone to cancellation here, we omit computing this secant (it's probably useless)
344  	       * unless the bounds are close to 0 as well (xub <= 0.1 in the "if" above)
345  	       * or we have exactly xlb=-xub, where we can return a clean 0.0 (though it's probably useless)
346  	       */
347  	      if( xlb == -xub )
348  	      {
349  	         *slope = 0.0;
350  	         *constant = pow(xlb, exponent);
351  	      }
352  	      else
353  	      {
354  	         return;
355  	      }
356  	   }
357  	   else if( xlb == 0.0 && exponent > 0.0 )
358  	   {
359  	      assert(xub >= 0.0);
360  	      *slope = pow(xub, exponent-1.0);
361  	      *constant = 0.0;
362  	   }
363  	   else if( xub == 0.0 && exponent > 0.0 )
364  	   {
365  	      /* normal pow: slope = - xlb^exponent / (-xlb) = xlb^(exponent-1)
366  	       * signpower:  slope = (-xlb)^exponent / (-xlb) = (-xlb)^(exponent-1)
367  	       */
368  	      assert(xlb <= 0.0);  /* so signpower or exponent is integral */
369  	      if( signpower )
370  	         *slope = pow(-xlb, exponent-1.0);
371  	      else
372  	         *slope = pow(xlb, exponent-1.0);
373  	      *constant = 0.0;
374  	   }
375  	   else if( SCIPisEQ(scip, xlb, xub) && (!signpower || xlb >= 0.0 || xub <= 0.0) )
376  	   {
377  	      /* Computing the slope as (xub^n - xlb^n)/(xub-xlb) can lead to cancellation.
378  	       * To avoid this, we replace xub^n by a Taylor expansion of pow at xlb:
379  	       *   xub^n = xlb^n + n xlb^(n-1) (xub-xlb) + 0.5 n*(n-1) xlb^(n-2) (xub-xlb)^2 + 1/6 n*(n-1)*(n-2) xi^(n-3) (xub-xlb)^3  for some xlb < xi < xub
380  	       * Dropping the last term, the slope is  (with an error of O((xub-xlb)^2) = 1e-18)
381  	       *   n*xlb^(n-1) + 0.5 n*(n-1) xlb^(n-2)*(xub-xlb)
382  	       * = n*xlb^(n-1) (1 - 0.5*(n-1)) + 0.5 n*(n-1) xlb^(n-2)*xub
383  	       * = 0.5*n*((3-n)*xlb^(n-1) + (n-1) xlb^(n-2)*xub)
384  	       *
385  	       * test n=2: 0.5*2*((3-2)*xlb + (2-1) 1*xub) = xlb + xub    ok
386  	       *      n=3: 0.5*3*((3-3)*xlb + (3-1) xlb*xub) = 3*xlb*xub ~ xlb^2 + xlb*xub + xub^2  ok
387  	       *
388  	       * The constant is
389  	       *   xlb^n - 0.5*n*((3-n) xlb^(n-1) + (n-1) xlb^(n-2)*xub) * xlb
390  	       * = xlb^n - 0.5*n*(3-n) xlb^n - 0.5*n*(n-1) xlb^(n-1)*xub
391  	       * = (1-0.5*n*(3-n)) xlb^n - 0.5 n*(n-1) xlb^(n-1) xub
392  	       *
393  	       * test n=2: (1-0.5*2*(3-2)) xlb^2 - 0.5 2*(2-1) xlb xub = -xlb*xub
394  	       *           old formula: xlb^2 - (xlb+xub) * xlb = -xlb*xub  ok
395  	       *
396  	       * For signpower with xub <= 0, we can negate xlb and xub:
397  	       *   slope: (sign(xub)|xub|^n - sign(xlb)*|xlb|^n) / (xub-xlb) = -((-xub)^n - (-xlb)^n) / (xub - xlb) = ((-xub)^n - (-xlb)^n) / (-xub - (-xlb))
398  	       *   constant: sign(xlb)|xlb|^n + slope * (xub - xlb) = -((-xlb)^n - slope * (xub - xlb)) = -((-xlb)^n + slope * ((-xub) - (-xlb)))
399  	       */
400  	      SCIP_Real xlb_n;   /* xlb^n */
401  	      SCIP_Real xlb_n1;  /* xlb^(n-1) */
402  	      SCIP_Real xlb_n2;  /* xlb^(n-2) */
403  	
404  	      if( signpower && xub <= 0.0 )
405  	      {
406  	         xlb *= -1.0;
407  	         xub *= -1.0;
408  	      }
409  	
410  	      xlb_n = pow(xlb, exponent);
411  	      xlb_n1 = pow(xlb, exponent - 1.0);
412  	      xlb_n2 = pow(xlb, exponent - 2.0);
413  	
414  	      *slope = 0.5*exponent * ((3.0-exponent) * xlb_n1 + (exponent-1.0) * xlb_n2 * xub);
415  	      *constant = (1.0 - 0.5*exponent*(3.0-exponent)) * xlb_n - 0.5*exponent*(exponent-1.0) * xlb_n1 * xub;
416  	
417  	      if( signpower && xub <= 0.0 )
418  	         *constant *= -1.0;
419  	   }
420  	   else
421  	   {
422  	      SCIP_Real lbval;
423  	      SCIP_Real ubval;
424  	
425  	      if( signpower )
426  	         lbval = SIGN(xlb) * pow(REALABS(xlb), exponent);
427  	      else
428  	         lbval = pow(xlb, exponent);
429  	      if( !SCIPisFinite(lbval) )
430  	         return;
431  	
432  	      if( signpower )
433  	         ubval = SIGN(xub) * pow(REALABS(xub), exponent);
434  	      else
435  	         ubval = pow(xub, exponent);
436  	      if( !SCIPisFinite(ubval) )
437  	         return;
438  	
439  	      /* we still can have bad numerics when xlb^exponent and xub^exponent are very close, but xlb and xub are not
440  	       * for now, only check that things did not cancel out completely
441  	       */
442  	      if( lbval == ubval )
443  	         return;
444  	
445  	      *slope = (ubval - lbval) / (xub - xlb);
446  	      *constant = lbval - *slope * xlb;
447  	   }
448  	
449  	   /* check whether we had overflows */
450  	   if( !SCIPisFinite(*slope) || !SCIPisFinite(*constant) )
451  	      return;
452  	
453  	   *success = TRUE;
454  	}
455  	
456  	/** Separation for parabola
457  	 *
458  	 * - even positive powers: x^2, x^4, x^6 with x arbitrary, or
459  	 * - positive powers > 1: x^1.5, x^2.5 with x >= 0
460  	 <pre>
461  	  100 +--------------------------------------------------------------------+
462  	      |*               +                 +                +               *|
463  	   90 |**                                                     x**2 ********|
464  	      |  *                                                              *  |
465  	   80 |-+*                                                              *+-|
466  	      |   **                                                          **   |
467  	   70 |-+   *                                                        *   +-|
468  	      |     **                                                      **     |
469  	   60 |-+     *                                                    *     +-|
470  	      |       **                                                  **       |
471  	   50 |-+       *                                                *       +-|
472  	      |          **                                            **          |
473  	   40 |-+          *                                          *          +-|
474  	      |             **                                      **             |
475  	   30 |-+            **                                    **            +-|
476  	      |                **                                **                |
477  	   20 |-+                **                            **                +-|
478  	      |                   ***                        ***                   |
479  	   10 |-+                   ***                    ***                   +-|
480  	      |                +       *****     +    *****       +                |
481  	    0 +--------------------------------------------------------------------+
482  	     -10              -5                 0                5                10
483  	 </pre>
484  	 */
485  	static
486  	void estimateParabola(
487  	   SCIP*                 scip,               /**< SCIP data structure */
488  	   SCIP_Real             exponent,           /**< exponent */
489  	   SCIP_Bool             overestimate,       /**< should the power be overestimated? */
490  	   SCIP_Real             xlb,                /**< lower bound on x */
491  	   SCIP_Real             xub,                /**< upper bound on x */
492  	   SCIP_Real             xref,               /**< reference point (where to linearize) */
493  	   SCIP_Real*            constant,           /**< buffer to store constant term of estimator */
494  	   SCIP_Real*            slope,              /**< buffer to store slope of estimator */
495  	   SCIP_Bool*            islocal,            /**< buffer to store whether estimator only locally valid, that is,
496  	                                              *   it depends on given bounds */
497  	   SCIP_Bool*            success             /**< buffer to store whether estimator could be computed */
498  	   )
499  	{
500  	   assert(scip != NULL);
501  	   assert(constant != NULL);
502  	   assert(slope != NULL);
503  	   assert(islocal != NULL);
504  	   assert(success != NULL);
505  	   assert((exponent >= 0.0 && EPSISINT(exponent/2.0, 0.0)) || (exponent > 1.0 && xlb >= 0.0));
506  	
507  	   if( !overestimate )
508  	   {
509  	      computeTangent(scip, FALSE, exponent, xref, constant, slope, success);
510  	      *islocal = FALSE;
511  	   }
512  	   else
513  	   {
514  	      /* overestimation -> secant */
515  	      computeSecant(scip, FALSE, exponent, xlb, xub, constant, slope, success);
516  	      *islocal = TRUE;
517  	   }
518  	}
519  	
520  	
521  	/** Separation for signpower
522  	 *
523  	 * - odd positive powers, x^3, x^5, x^7
524  	 * - sign(x)|x|^n for n > 1
525  	 * - lower bound on x is negative (otherwise one should use separation for parabola)
526  	 <pre>
527  	  100 +--------------------------------------------------------------------+
528  	      |                +                 +                +              **|
529  	      |                                                   x*abs(x) ******* |
530  	      |                                                              **    |
531  	      |                                                            **      |
532  	   50 |-+                                                       ***      +-|
533  	      |                                                      ***           |
534  	      |                                                   ***              |
535  	      |                                               *****                |
536  	      |                                          *****                     |
537  	    0 |-+                        ****************                        +-|
538  	      |                     *****                                          |
539  	      |                *****                                               |
540  	      |              ***                                                   |
541  	      |           ***                                                      |
542  	  -50 |-+      ***                                                       +-|
543  	      |      **                                                            |
544  	      |    **                                                              |
545  	      |  **                                                                |
546  	      |**              +                 +                +                |
547  	 -100 +--------------------------------------------------------------------+
548  	     -10              -5                 0                5                10
549  	 </pre>
550  	 */
551  	static
552  	void estimateSignedpower(
553  	   SCIP*                 scip,               /**< SCIP data structure */
554  	   SCIP_Real             exponent,           /**< exponent */
555  	   SCIP_Real             root,               /**< positive root of the polynomial (n-1) y^n + n y^(n-1) - 1,
556  	                                              *   if xubglobal > 0 */
557  	   SCIP_Bool             overestimate,       /**< should the power be overestimated? */
558  	   SCIP_Real             xlb,                /**< lower bound on x, assumed to be non-positive */
559  	   SCIP_Real             xub,                /**< upper bound on x */
560  	   SCIP_Real             xref,               /**< reference point (where to linearize) */
561  	   SCIP_Real             xlbglobal,          /**< global lower bound on x */
562  	   SCIP_Real             xubglobal,          /**< global upper bound on x */
563  	   SCIP_Real*            constant,           /**< buffer to store constant term of estimator */
564  	   SCIP_Real*            slope,              /**< buffer to store slope of estimator */
565  	   SCIP_Bool*            islocal,            /**< buffer to store whether estimator only locally valid, that is,
566  	                                              *   it depends on given bounds */
567  	   SCIP_Bool*            branchcand,         /**< buffer to indicate whether estimator would improve by branching
568  	                                              *   on it */
569  	   SCIP_Bool*            success             /**< buffer to store whether estimator could be computed */
570  	   )
571  	{
572  	   assert(scip != NULL);
573  	   assert(constant != NULL);
574  	   assert(slope != NULL);
575  	   assert(islocal != NULL);
576  	   assert(branchcand != NULL);
577  	   assert(*branchcand == TRUE);  /* the default */
578  	   assert(success != NULL);
579  	   assert(exponent >= 1.0);
580  	   assert(xlb < 0.0); /* otherwise estimateParabola should have been called */
581  	   assert(xubglobal <= 0.0 || (root > 0.0 && root < 1.0));
582  	
583  	   *success = FALSE;
584  	
585  	   if( !SCIPisPositive(scip, xub) )
586  	   {
587  	      /* easy case */
588  	      if( !overestimate )
589  	      {
590  	         /* underestimator is secant */
591  	         computeSecant(scip, TRUE, exponent, xlb, xub, constant, slope, success);
592  	         *islocal = TRUE;
593  	      }
594  	      else
595  	      {
596  	         /* overestimator is tangent */
597  	
598  	         /* we must linearize left of 0 */
599  	         if( xref > 0.0 )
600  	            xref = 0.0;
601  	
602  	         computeTangent(scip, TRUE, exponent, xref, constant, slope, success);
603  	
604  	         /* if global upper bound is > 0, then the tangent is only valid locally if the reference point is right of
605  	          * -root*xubglobal
606  	          */
607  	         *islocal = SCIPisPositive(scip, xubglobal) && xref > -root * xubglobal;
608  	
609  	         /* tangent doesn't move after branching */
610  	         *branchcand = FALSE;
611  	      }
612  	   }
613  	   else
614  	   {
615  	      SCIP_Real c;
616  	
617  	      if( !overestimate )
618  	      {
619  	         /* compute the special point which decides between secant and tangent */
620  	         c = -xlb * root;
621  	
622  	         if( xref < c )
623  	         {
624  	            /* underestimator is secant between xlb and c */
625  	            computeSecant(scip, TRUE, exponent, xlb, c, constant, slope, success);
626  	            *islocal = TRUE;
627  	         }
628  	         else
629  	         {
630  	            /* underestimator is tangent */
631  	            computeTangent(scip, TRUE, exponent, xref, constant, slope, success);
632  	
633  	            /* if reference point is left of -root*xlbglobal (c w.r.t. global bounds),
634  	             * then tangent is not valid w.r.t. global bounds
635  	             */
636  	            *islocal = xref < -root * xlbglobal;
637  	
638  	            /* tangent doesn't move after branching */
639  	            *branchcand = FALSE;
640  	         }
641  	      }
642  	      else
643  	      {
644  	         /* compute the special point which decides between secant and tangent */
645  	         c = -xub * root;
646  	
647  	         if( xref <= c )
648  	         {
649  	            /* overestimator is tangent */
650  	            computeTangent(scip, TRUE, exponent, xref, constant, slope, success);
651  	
652  	            /* if reference point is right of -root*xubglobal (c w.r.t. global bounds),
653  	             * then tangent is not valid w.r.t. global bounds
654  	             */
655  	            *islocal = xref > -root * xubglobal;
656  	
657  	            /* tangent doesn't move after branching */
658  	            *branchcand = FALSE;
659  	         }
660  	         else
661  	         {
662  	            /* overestimator is secant */
663  	            computeSecant(scip, TRUE, exponent, c, xub, constant, slope, success);
664  	            *islocal = TRUE;
665  	         }
666  	      }
667  	   }
668  	}
669  	
670  	/** Separation for positive hyperbola
671  	 *
672  	 * - x^-2, x^-4 with x arbitrary
673  	 * - x^-0.5, x^-1, x^-1.5, x^-3, x^-5 with x >= 0
674  	 <pre>
675  	  5 +----------------------------------------------------------------------+
676  	    |                 +               * +*               +                 |
677  	    |                                 *  *                 x**(-2) ******* |
678  	  4 |-+                               *  *                               +-|
679  	    |                                 *  *                                 |
680  	    |                                 *  *                                 |
681  	    |                                 *  *                                 |
682  	  3 |-+                               *   *                              +-|
683  	    |                                *    *                                |
684  	    |                                *    *                                |
685  	  2 |-+                              *    *                              +-|
686  	    |                                *    *                                |
687  	    |                               *      *                               |
688  	  1 |-+                             *      *                             +-|
689  	    |                               *      *                               |
690  	    |                             **        **                             |
691  	    |                   **********            **********                   |
692  	  0 |*******************                                *******************|
693  	    |                                                                      |
694  	    |                 +                 +                +                 |
695  	 -1 +----------------------------------------------------------------------+
696  	   -10               -5                 0                5                 10
697  	 </pre>
698  	 */
699  	static
700  	void estimateHyperbolaPositive(
701  	   SCIP*                 scip,               /**< SCIP data structure */
702  	   SCIP_Real             exponent,           /**< exponent */
703  	   SCIP_Real             root,               /**< negative root of the polynomial (n-1) y^n - n y^(n-1) + 1,
704  	                                              *   if x has mixed sign (w.r.t. global bounds?) and underestimating */
705  	   SCIP_Bool             overestimate,       /**< should the power be overestimated? */
706  	   SCIP_Real             xlb,                /**< lower bound on x */
707  	   SCIP_Real             xub,                /**< upper bound on x */
708  	   SCIP_Real             xref,               /**< reference point (where to linearize) */
709  	   SCIP_Real             xlbglobal,          /**< global lower bound on x */
710  	   SCIP_Real             xubglobal,          /**< global upper bound on x */
711  	   SCIP_Real*            constant,           /**< buffer to store constant term of estimator */
712  	   SCIP_Real*            slope,              /**< buffer to store slope of estimator */
713  	   SCIP_Bool*            islocal,            /**< buffer to store whether estimator only locally valid, that is,
714  	                                              *   it depends on given bounds */
715  	   SCIP_Bool*            branchcand,         /**< buffer to indicate whether estimator would improve by branching
716  	                                              *   on it */
717  	   SCIP_Bool*            success             /**< buffer to store whether estimator could be computed */
718  	   )
719  	{
720  	   assert(scip != NULL);
721  	   assert(constant != NULL);
722  	   assert(slope != NULL);
723  	   assert(islocal != NULL);
724  	   assert(branchcand != NULL);
725  	   assert(*branchcand == TRUE);  /* the default */
726  	   assert(success != NULL);
727  	   assert(exponent < 0.0);
728  	   assert(EPSISINT(exponent/2.0, 0.0) || xlb >= 0.0);
729  	
730  	   *success = FALSE;
731  	
732  	   if( !overestimate )
733  	   {
734  	      if( xlb >= 0.0 || xub <= 0.0 )
735  	      {
736  	         /* underestimate and fixed sign -> tangent */
737  	
738  	         /* make sure xref has the same sign as xlb,xub */
739  	         if( xref < 0.0 && xlb >= 0.0 )
740  	            xref = xlb;
741  	         else if( xref > 0.0 && xub <= 0.0 )
742  	            xref = xub;
743  	
744  	         if( SCIPisZero(scip, xref) )
745  	         {
746  	            /* estimator would need to have an (essentially) infinite scope
747  	             * first try to make up a better refpoint
748  	             */
749  	            if( xub > 0.0 )
750  	            {
751  	               /* thus xlb >= 0.0; stay close to xlb (probably = 0) */
752  	               if( !SCIPisInfinity(scip, xub) )
753  	                  xref = 0.9 * xlb + 0.1 * xub;
754  	               else
755  	                  xref = 0.1;
756  	            }
757  	            else
758  	            {
759  	               /* xub <= 0.0; stay close to xub (probably = 0) */
760  	               if( !SCIPisInfinity(scip, -xlb) )
761  	                  xref = 0.1 * xlb + 0.9 * xub;
762  	               else
763  	                  xref = 0.1;
764  	            }
765  	
766  	            /* if still close to 0, then also bounds are close to 0, then just give up */
767  	            if( SCIPisZero(scip, xref) )
768  	               return;
769  	         }
770  	
771  	         computeTangent(scip, FALSE, exponent, xref, constant, slope, success);
772  	
773  	         /* tangent will not change if branching on x (even if only locally valid, see checks below) */
774  	         *branchcand = FALSE;
775  	
776  	         if( EPSISINT(exponent/2.0, 0.0) )
777  	         {
778  	            /* for even exponents (as in the picture):
779  	             * if x has fixed sign globally, then our tangent is also globally valid
780  	             * however, if x has mixed sign, then it depends on the constellation between reference point and global
781  	             * bounds, whether the tangent is globally valid (see also the longer discussion for the mixed-sign
782  	             * underestimator below )
783  	             */
784  	            if( xref > 0.0 && xlbglobal < 0.0 )
785  	            {
786  	               assert(xubglobal > 0.0);  /* since xref > 0.0 */
787  	               assert(root < 0.0); /* root needs to be given */
788  	               /* if on right side, then tangent is only locally valid if xref is too much to the left */
789  	               *islocal = xref < xlbglobal * root;
790  	            }
791  	            else if( xref < 0.0 && xubglobal > 0.0 )
792  	            {
793  	               assert(xlbglobal < 0.0);  /* since xref < 0.0 */
794  	               assert(root < 0.0); /* root needs to be given */
795  	               /* if on left side, then tangent is only locally valid if xref is too much to the right */
796  	               *islocal = xref > xubglobal * root;
797  	            }
798  	            else
799  	               *islocal = FALSE;
800  	         }
801  	         else
802  	         {
803  	            /* for odd exponents, the tangent is only locally valid if the sign of x is not fixed globally */
804  	            *islocal = xlbglobal * xubglobal < 0.0;
805  	         }
806  	      }
807  	      else
808  	      {
809  	         /* underestimate but mixed sign */
810  	         if( SCIPisInfinity(scip, -xlb) )
811  	         {
812  	            if( SCIPisInfinity(scip, xub) )
813  	            {
814  	               /* underestimator is constant 0, but that is globally valid */
815  	               *constant = 0.0;
816  	               *slope = 0.0;
817  	               *islocal = FALSE;
818  	               *success = TRUE;
819  	               return;
820  	            }
821  	
822  	            /* switch sign of x (mirror on ordinate) to make left bound finite and use its estimator */
823  	            estimateHyperbolaPositive(scip, exponent, root, overestimate, -xub, -xlb, -xref, -xubglobal, -xlbglobal,
824  	                                      constant, slope, islocal, branchcand, success);
825  	            if( *success )
826  	               *slope = -*slope;
827  	         }
828  	         else
829  	         {
830  	            /* The convex envelope of x^exponent for x in [xlb, infinity] is a line (secant) between xlb and some positive
831  	             * coordinate xhat, and x^exponent for x > xhat.
832  	             * Further, on [xlb,xub] with xub < xhat, the convex envelope is the secant between xlb and xub.
833  	             *
834  	             * To find xhat, consider the affine-linear function  l(x) = xlb^n + c * (x - xlb)   where n = exponent
835  	             * we look for a value of x such that f(x) and l(x) coincide and such that l(x) will be tangent to f(x) on that
836  	             * point, that is
837  	             * xhat > 0 such that f(xhat) = l(xhat) and f'(xhat) = l'(xhat)
838  	             * => xhat^n = xlb^n + c * (xhat - xlb)   and   n * xhat^(n-1) = c
839  	             * => xhat^n = xlb^n + n * xhat^n - n * xhat^(n-1) * xlb
840  	             * => 0 = xlb^n + (n-1) * xhat^n - n * xhat^(n-1) * xlb
841  	             *
842  	             * Divide by xlb^n, one gets a polynomial that looks very much like the one for signpower, but a sign is
843  	             * different (since this is *not signed* power):
844  	             * 0 = 1 + (n-1) * y^n - n * y^(n-1)  where y = xhat/xlb
845  	             *
846  	             * The solution y < 0 (because xlb < 0 and we want xhat > 0) is what we expect to be given as "root".
847  	             */
848  	            assert(root < 0.0); /* root needs to be given */
849  	            if( xref <= xlb * root )
850  	            {
851  	               /* If the reference point is left of xhat (=xlb*root), then we can take the
852  	                * secant between xlb and root*xlb (= tangent at root*xlb).
853  	                * However, if xub < root*xlb, then we can tilt the estimator to be the secant between xlb and xub.
854  	                */
855  	               computeSecant(scip, FALSE, exponent, xlb, MIN(xlb * root, xub), constant, slope, success);
856  	               *islocal = TRUE;
857  	            }
858  	            else
859  	            {
860  	               /* If reference point is right of xhat, then take the tangent at xref.
861  	                * This will still be underestimating for x in [xlb,0], too.
862  	                * The tangent is globally valid, if we had also generated w.r.t. global bounds.
863  	                */
864  	               computeTangent(scip, FALSE, exponent, xref, constant, slope, success);
865  	               *islocal = xref < xlbglobal * root;
866  	               *branchcand = FALSE;
867  	            }
868  	         }
869  	      }
870  	   }
871  	   else
872  	   {
873  	      /* overestimate and mixed sign -> pole is within domain -> cannot overestimate */
874  	      if( xlb < 0.0 && xub > 0.0 )
875  	         return;
876  	
877  	      /* overestimate and fixed sign -> secant */
878  	      computeSecant(scip, FALSE, exponent, xlb, xub, constant, slope, success);
879  	      *islocal = TRUE;
880  	   }
881  	
882  	}
883  	
884  	/** Separation for mixed-sign hyperbola
885  	 *
886  	 * - x^-1, x^-3, x^-5 without x >= 0 (either x arbitrary or x negative)
887  	 <pre>
888  	    +----------------------------------------------------------------------+
889  	    |                 +                 *                +                 |
890  	  4 |-+                                  *                 x**(-1) *******-|
891  	    |                                    *                                 |
892  	    |                                    *                                 |
893  	    |                                    *                                 |
894  	  2 |-+                                  *                               +-|
895  	    |                                     *                                |
896  	    |                                      **                              |
897  	    |                                        *********                     |
898  	  0 |*********************                            *********************|
899  	    |                     *********                                        |
900  	    |                              **                                      |
901  	    |                                *                                     |
902  	 -2 |-+                               *                                  +-|
903  	    |                                 *                                    |
904  	    |                                 *                                    |
905  	    |                                 *                                    |
906  	 -4 |-+                               *                                  +-|
907  	    |                 +                *+                +                 |
908  	    +----------------------------------------------------------------------+
909  	   -10               -5                 0                5                 10
910  	 </pre>
911  	 */
912  	static
913  	void estimateHyperbolaMixed(
914  	   SCIP*                 scip,               /**< SCIP data structure */
915  	   SCIP_Real             exponent,           /**< exponent */
916  	   SCIP_Bool             overestimate,       /**< should the power be overestimated? */
917  	   SCIP_Real             xlb,                /**< lower bound on x */
918  	   SCIP_Real             xub,                /**< upper bound on x */
919  	   SCIP_Real             xref,               /**< reference point (where to linearize) */
920  	   SCIP_Real             xlbglobal,          /**< global lower bound on x */
921  	   SCIP_Real             xubglobal,          /**< global upper bound on x */
922  	   SCIP_Real*            constant,           /**< buffer to store constant term of estimator */
923  	   SCIP_Real*            slope,              /**< buffer to store slope of estimator */
924  	   SCIP_Bool*            islocal,            /**< buffer to store whether estimator only locally valid, that is,
925  	                                                  it depends on given bounds */
926  	   SCIP_Bool*            branchcand,         /**< buffer to indicate whether estimator would improve by branching
927  	                                                  on it */
928  	   SCIP_Bool*            success             /**< buffer to store whether estimator could be computed */
929  	   )
930  	{
931  	   assert(scip != NULL);
932  	   assert(constant != NULL);
933  	   assert(slope != NULL);
934  	   assert(islocal != NULL);
935  	   assert(branchcand != NULL);
936  	   assert(*branchcand == TRUE);  /* the default */
937  	   assert(success != NULL);
938  	   assert(exponent < 0.0);
939  	   assert(EPSISINT((exponent-1.0)/2.0, 0.0));
940  	   assert(xlb < 0.0);
941  	
942  	   *success = FALSE;
943  	
944  	   if( xub <= 0.0 )
945  	   {
946  	      /* x is negative */
947  	      if( !overestimate )
948  	      {
949  	         /* underestimation -> secant */
950  	         computeSecant(scip, FALSE, exponent, xlb, xub, constant, slope, success);
951  	         *islocal = TRUE;
952  	      }
953  	      else if( !SCIPisZero(scip, xlb/10.0) )
954  	      {
955  	         /* overestimation -> tangent */
956  	
957  	         /* need to linearize left of 0 */
958  	         if( xref > 0.0 )
959  	            xref = 0.0;
960  	
961  	         if( SCIPisZero(scip, xref) )
962  	         {
963  	            /* if xref is very close to 0.0, then slope would be infinite
964  	             * try to move closer to lower bound (if xlb < -10*eps)
965  	             */
966  	            if( !SCIPisInfinity(scip, -xlb) )
967  	               xref = 0.1*xlb + 0.9*xub;
968  	            else
969  	               xref = 0.1;
970  	         }
971  	
972  	         computeTangent(scip, FALSE, exponent, xref, constant, slope, success);
973  	
974  	         /* if x does not have a fixed sign globally, then our tangent is not globally valid
975  	          * (power is not convex on global domain)
976  	          */
977  	         *islocal = xlbglobal * xubglobal < 0.0;
978  	
979  	         /* tangent doesn't move by branching */
980  	         *branchcand = FALSE;
981  	      }
982  	      /* else: xlb is very close to zero, xub is <= 0, so slope would be infinite
983  	       * (for any reference point in [xlb, xub]) -> do not estimate
984  	       */
985  	   }
986  	   /* else: x has mixed sign -> pole is within domain -> cannot estimate */
987  	}
988  	
989  	/** builds an estimator for a power function */
990  	static
991  	SCIP_RETCODE buildPowEstimator(
992  	   SCIP*                 scip,               /**< SCIP data structure */
993  	   SCIP_EXPRDATA*        exprdata,           /**< expression data */
994  	   SCIP_Bool             overestimate,       /**< is this an overestimator? */
995  	   SCIP_Real             childlb,            /**< local lower bound on the child */
996  	   SCIP_Real             childub,            /**< local upper bound on the child */
997  	   SCIP_Real             childglb,           /**< global lower bound on the child */
998  	   SCIP_Real             childgub,           /**< global upper bound on the child */
999  	   SCIP_Bool             childintegral,      /**< whether child is integral */
1000 	   SCIP_Real             refpoint,           /**< reference point */
1001 	   SCIP_Real             exponent,           /**< esponent */
1002 	   SCIP_Real*            coef,               /**< pointer to store the coefficient of the estimator */
1003 	   SCIP_Real*            constant,           /**< pointer to store the constant of the estimator */
1004 	   SCIP_Bool*            success,            /**< pointer to store whether the estimator was built successfully */
1005 	   SCIP_Bool*            islocal,            /**< pointer to store whether the estimator is valid w.r.t. local bounds
1006 	                                              *   only */
1007 	   SCIP_Bool*            branchcand          /**< pointer to indicate whether to consider child for branching
1008 	                                              *   (initialized to TRUE) */
1009 	   )
1010 	{
1011 	   SCIP_Bool isinteger;
1012 	   SCIP_Bool iseven;
1013 	
1014 	   assert(scip != NULL);
1015 	   assert(exprdata != NULL);
1016 	   assert(coef != NULL);
1017 	   assert(constant != NULL);
1018 	   assert(success != NULL);
1019 	   assert(islocal != NULL);
1020 	   assert(branchcand != NULL);
1021 	
1022 	   isinteger = EPSISINT(exponent, 0.0);
1023 	   iseven = isinteger && EPSISINT(exponent / 2.0, 0.0);
1024 	
1025 	   if( exponent == 2.0 )
1026 	   {
1027 	      /* important special case: quadratic case */
1028 	      /* initialize, because SCIPaddSquareXyz only adds to existing values */
1029 	      *success = TRUE;
1030 	      *coef = 0.0;
1031 	      *constant = 0.0;
1032 	
1033 	      if( overestimate )
1034 	      {
1035 	         SCIPaddSquareSecant(scip, 1.0, childlb, childub, coef, constant, success);
1036 	         *islocal = TRUE; /* secants are only valid locally */
1037 	      }
1038 	      else
1039 	      {
1040 	         SCIPaddSquareLinearization(scip, 1.0, refpoint, childintegral, coef, constant, success);
1041 	         *islocal = FALSE; /* linearizations are globally valid */
1042 	         *branchcand = FALSE;  /* there is no improvement due to branching */
1043 	      }
1044 	   }
1045 	   else if( exponent > 0.0 && iseven )
1046 	   {
1047 	      estimateParabola(scip, exponent, overestimate, childlb, childub, refpoint, constant, coef, islocal, success);
1048 	      /* if estimate is locally valid, then we computed a secant and so branching can improve it */
1049 	      *branchcand = *islocal;
1050 	   }
1051 	   else if( exponent > 1.0 && childlb >= 0.0 )
1052 	   {
1053 	      /* make sure we linearize in convex region (if we will linearize) */
1054 	      if( refpoint < 0.0 )
1055 	         refpoint = 0.0;
1056 	
1057 	      estimateParabola(scip, exponent, overestimate, childlb, childub, refpoint, constant, coef, islocal, success);
1058 	
1059 	      /* if estimate is locally valid, then we computed a secant and so branching can improve it */
1060 	      *branchcand = *islocal;
1061 	
1062 	      /* if odd power, then check whether tangent on parabola is also globally valid, that is reference point is
1063 	       * right of -root*global-lower-bound
1064 	       */
1065 	      if( !*islocal && !iseven && childglb < 0.0 )
1066 	      {
1067 	         if( SCIPisInfinity(scip, -childglb) )
1068 	            *islocal = TRUE;
1069 	         else
1070 	         {
1071 	            if( exprdata->root == SCIP_INVALID )
1072 	            {
1073 	               SCIP_CALL( computeSignpowerRoot(scip, &exprdata->root, exponent) );
1074 	            }
1075 	            *islocal = refpoint < exprdata->root * (-childglb);
1076 	         }
1077 	      }
1078 	   }
1079 	   else if( exponent > 1.0 )  /* and !iseven && childlb < 0.0 due to previous if */
1080 	   {
1081 	      /* compute root if not known yet; only needed if mixed sign (global child ub > 0) */
1082 	      if( exprdata->root == SCIP_INVALID && childgub > 0.0 )
1083 	      {
1084 	         SCIP_CALL( computeSignpowerRoot(scip, &exprdata->root, exponent) );
1085 	      }
1086 	      estimateSignedpower(scip, exponent, exprdata->root, overestimate, childlb, childub, refpoint,
1087 	                          -childglb, childgub, constant, coef, islocal, branchcand, success);
1088 	   }
1089 	   else if( exponent < 0.0 && (iseven || childlb >= 0.0) )
1090 	   {
1091 	      /* compute root if not known yet; only needed if mixed sign (globally) and iseven */
1092 	      if( exprdata->root == SCIP_INVALID && iseven )
1093 	      {
1094 	         SCIP_CALL( computeHyperbolaRoot(scip, &exprdata->root, exponent) );
1095 	      }
1096 	      estimateHyperbolaPositive(scip, exponent, exprdata->root, overestimate, childlb, childub, refpoint,
1097 	                                childglb, childgub, constant, coef, islocal, branchcand, success);
1098 	   }
1099 	   else if( exponent < 0.0 )
1100 	   {
1101 	      assert(!iseven); /* should hold due to previous if */
1102 	      assert(childlb < 0.0); /* should hold due to previous if */
1103 	      assert(isinteger); /* should hold because childlb < 0.0 (same as assert above) */
1104 	
1105 	      estimateHyperbolaMixed(scip, exponent, overestimate, childlb, childub, refpoint, childglb, childgub,
1106 	                             constant, coef, islocal, branchcand, success);
1107 	   }
1108 	   else
1109 	   {
1110 	      assert(exponent < 1.0); /* the only case that should be left */
1111 	      assert(exponent > 0.0); /* should hold due to previous if */
1112 	
1113 	      SCIPestimateRoot(scip, exponent, overestimate, childlb, childub, refpoint, constant, coef, islocal, success);
1114 	
1115 	      /* if estimate is locally valid, then we computed a secant, and so branching can improve it */
1116 	      *branchcand = *islocal;
1117 	   }
1118 	
1119 	   return SCIP_OKAY;
1120 	}
1121 	
1122 	/** fills an array of reference points for estimating on the convex side */
1123 	static
1124 	void addTangentRefpoints(
1125 	   SCIP*                 scip,               /**< SCIP data structure */
1126 	   SCIP_Real             exponent,           /**< exponent of the power expression */
1127 	   SCIP_Real             lb,                 /**< lower bound on the child variable */
1128 	   SCIP_Real             ub,                 /**< upper bound on the child variable */
1129 	   SCIP_Real*            refpoints           /**< array to store the reference points */
1130 	   )
1131 	{
1132 	   SCIP_Real maxabsbnd;
1133 	
1134 	   assert(refpoints != NULL);
1135 	
1136 	   maxabsbnd = pow(INITLPMAXPOWVAL, 1 / exponent);
1137 	
1138 	   /* make sure the absolute values of bounds are not too large */
1139 	   if( ub > -maxabsbnd )
1140 	      lb = MAX(lb, -maxabsbnd);
1141 	   if( lb < maxabsbnd )
1142 	      ub = MIN(ub,  maxabsbnd);
1143 	
1144 	   /* in the case when ub < -maxabsbnd or lb > maxabsbnd, we still want to at least make bounds finite */
1145 	   if( SCIPisInfinity(scip, -lb) )
1146 	      lb = MIN(-10.0, ub - 0.1*REALABS(ub));  /*lint !e666 */
1147 	   if( SCIPisInfinity(scip,  ub) )
1148 	      ub = MAX( 10.0, lb + 0.1*REALABS(lb));  /*lint !e666 */
1149 	
1150 	   refpoints[0] = (7.0 * lb + ub) / 8.0;
1151 	   refpoints[1] = (lb + ub) / 2.0;
1152 	   refpoints[2] = (lb + 7.0 * ub) / 8.0;
1153 	}
1154 	
1155 	/** fills an array of reference points for sign(x)*abs(x)^n or x^n (n odd), where x has mixed signs
1156 	 *
1157 	 *  The reference points are: the lower and upper bounds (one for secant and one for tangent);
1158 	 *  and for the second tangent, the point on the convex part of the function between the point
1159 	 *  deciding between tangent and secant, and the corresponding bound
1160 	 */
1161 	static
1162 	SCIP_RETCODE addSignpowerRefpoints(
1163 	   SCIP*                 scip,               /**< SCIP data structure */
1164 	   SCIP_EXPRDATA*        exprdata,           /**< expression data */
1165 	   SCIP_Real             lb,                 /**< lower bound on the child variable */
1166 	   SCIP_Real             ub,                 /**< upper bound on the child variable */
1167 	   SCIP_Real             exponent,           /**< exponent */
1168 	   SCIP_Bool             underestimate,      /**< are the refpoints for an underestimator */
1169 	   SCIP_Real*            refpoints           /**< array to store the reference points */
1170 	   )
1171 	{
1172 	   assert(refpoints != NULL);
1173 	
1174 	   if( (underestimate && SCIPisInfinity(scip, -lb)) || (!underestimate && SCIPisInfinity(scip, ub)) )
1175 	      return SCIP_OKAY;
1176 	
1177 	   if( exprdata->root == SCIP_INVALID )
1178 	   {
1179 	      SCIP_CALL( computeSignpowerRoot(scip, &exprdata->root, exponent) );
1180 	   }
1181 	
1182 	   /* make bounds finite (due to a previous if, only one can be infinite here) */
1183 	   if( SCIPisInfinity(scip, -lb) )
1184 	      lb = -ub * exprdata->root - 1.0;
1185 	   else if( SCIPisInfinity(scip,  ub) )
1186 	      ub = -lb * exprdata->root + 1.0;
1187 	
1188 	   if( underestimate )
1189 	   {
1190 	      /* secant point */
1191 	      refpoints[0] = lb;
1192 	
1193 	      /* tangent points, depending on the special point */
1194 	      if( -lb * exprdata->root < ub - 2.0 )
1195 	         refpoints[2] = ub;
1196 	      if( -lb * exprdata->root < ub - 4.0 )
1197 	         refpoints[1] = (-lb * exprdata->root + ub) / 2.0;
1198 	   }
1199 	
1200 	   if( !underestimate )
1201 	   {
1202 	      /* secant point */
1203 	      refpoints[2] = ub;
1204 	
1205 	      /* tangent points, depending on the special point */
1206 	      if( -ub * exprdata->root > lb + 2.0 )
1207 	         refpoints[0] = lb;
1208 	      if( -ub * exprdata->root > lb + 4.0 )
1209 	         refpoints[1] = (lb - ub * exprdata->root) / 2.0;
1210 	   }
1211 	
1212 	   return SCIP_OKAY;
1213 	}
1214 	
1215 	/** choose reference points for adding initestimates cuts for a power expression */
1216 	static
1217 	SCIP_RETCODE chooseRefpointsPow(
1218 	   SCIP*                 scip,               /**< SCIP data structure */
1219 	   SCIP_EXPRDATA*        exprdata,           /**< expression data */
1220 	   SCIP_Real             lb,                 /**< lower bound on the child variable */
1221 	   SCIP_Real             ub,                 /**< upper bound on the child variable */
1222 	   SCIP_Real*            refpointsunder,     /**< array to store reference points for underestimators */
1223 	   SCIP_Real*            refpointsover,      /**< array to store reference points for overestimators */
1224 	   SCIP_Bool             underestimate,      /**< whether refpoints for underestimation are needed */
1225 	   SCIP_Bool             overestimate        /**< whether refpoints for overestimation are needed */
1226 	   )
1227 	{
1228 	   SCIP_Bool convex;
1229 	   SCIP_Bool concave;
1230 	   SCIP_Bool mixedsign;
1231 	   SCIP_Bool even;
1232 	   SCIP_Real exponent;
1233 	
1234 	   assert(scip != NULL);
1235 	   assert(exprdata != NULL);
1236 	   assert(refpointsunder != NULL && refpointsover != NULL);
1237 	
1238 	   exponent = exprdata->exponent;
1239 	   even = EPSISINT(exponent, 0.0) && EPSISINT(exponent / 2.0, 0.0);
1240 	
1241 	   convex = FALSE;
1242 	   concave = FALSE;
1243 	   mixedsign = lb < 0.0 && ub > 0.0;
1244 	
1245 	   /* convex case:
1246 	    * - parabola with an even degree or positive domain
1247 	    * - hyperbola with a positive domain
1248 	    * - even hyperbola with a negative domain
1249 	    */
1250 	   if( (exponent > 1.0 && (lb >= 0 || even)) || (exponent < 0.0 && lb >= 0) || (exponent < 0.0 && even && ub <= 0.0) )
1251 	      convex = TRUE;
1252 	   /* concave case:
1253 	    * - parabola or hyperbola with a negative domain and (due to previous if) an uneven degree
1254 	    * - root
1255 	    */
1256 	   else if( ub <= 0 || (exponent > 0.0 && exponent < 1.0) )
1257 	      concave = TRUE;
1258 	
1259 	   if( underestimate )
1260 	   {
1261 	      if( convex )
1262 	         addTangentRefpoints(scip, exponent, lb, ub, refpointsunder);
1263 	      else if( (concave && !SCIPisInfinity(scip, -lb) && !SCIPisInfinity(scip, ub)) ||
1264 	               (exponent < 0.0 && even && mixedsign) ) /* concave with finite bounds or mixed even hyperbola */
1265 	      {
1266 	         /* for secant, refpoint doesn't matter, but we add it to signal that the corresponding cut should be created */
1267 	         refpointsunder[0] = (lb + ub) / 2.0;
1268 	      }
1269 	      else if( exponent > 1.0 && !even && mixedsign ) /* mixed signpower */
1270 	      {
1271 	         SCIP_CALL( addSignpowerRefpoints(scip, exprdata, lb, ub, exponent, TRUE, refpointsunder) );
1272 	      }
1273 	      else /* mixed odd hyperbola or an infinite bound */
1274 	         assert((exponent < 0.0 && !even && mixedsign) || SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub));
1275 	   }
1276 	
1277 	   if( overestimate )
1278 	   {
1279 	      if( convex && !SCIPisInfinity(scip, -lb) && !SCIPisInfinity(scip, ub) )
1280 	         refpointsover[0] = (lb + ub) / 2.0;
1281 	      else if( concave )
1282 	         addTangentRefpoints(scip, exponent, lb, ub, refpointsover);
1283 	      else if( exponent > 1.0 && !even && mixedsign ) /* mixed signpower */
1284 	      {
1285 	         SCIP_CALL( addSignpowerRefpoints(scip, exprdata, lb, ub, exponent, FALSE, refpointsover) );
1286 	      }
1287 	      else /* mixed hyperbola or an infinite bound */
1288 	         assert((exponent < 0.0 && mixedsign) || SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub));
1289 	   }
1290 	
1291 	   return SCIP_OKAY;
1292 	}
1293 	
1294 	
1295 	/*
1296 	 * Callback methods of expression handler
1297 	 */
1298 	
1299 	/** compares two power expressions
1300 	 *
1301 	 * the order of two power (normal or signed) is base_1^expo_1 < base_2^expo_2 if and only if
1302 	 * base_1 < base2 or, base_1 = base_2 and expo_1 < expo_2
1303 	 */
1304 	static
1305 	SCIP_DECL_EXPRCOMPARE(comparePow)
1306 	{  /*lint --e{715}*/
1307 	   SCIP_Real expo1;
1308 	   SCIP_Real expo2;
1309 	   int compareresult;
1310 	
1311 	   /**! [SnippetExprComparePow] */
1312 	   compareresult = SCIPcompareExpr(scip, SCIPexprGetChildren(expr1)[0], SCIPexprGetChildren(expr2)[0]);
1313 	   if( compareresult != 0 )
1314 	      return compareresult;
1315 	
1316 	   expo1 = SCIPgetExponentExprPow(expr1);
1317 	   expo2 = SCIPgetExponentExprPow(expr2);
1318 	
1319 	   return expo1 == expo2 ? 0 : expo1 < expo2 ? -1 : 1;
1320 	   /**! [SnippetExprComparePow] */
1321 	}
1322 	
1323 	/** simplifies a pow expression
1324 	 *
1325 	 * Evaluates the power function when its child is a value expression
1326 	 */
1327 	static
1328 	SCIP_DECL_EXPRSIMPLIFY(simplifyPow)
1329 	{  /*lint --e{715}*/
1330 	   SCIP_EXPRHDLR* exprhdlr;
1331 	   SCIP_EXPRHDLRDATA* exprhdlrdata;
1332 	   SCIP_EXPR* base;
1333 	   SCIP_Real exponent;
1334 	
1335 	   assert(scip != NULL);
1336 	   assert(expr != NULL);
1337 	   assert(simplifiedexpr != NULL);
1338 	   assert(SCIPexprGetNChildren(expr) == 1);
1339 	
1340 	   exprhdlr = SCIPexprGetHdlr(expr);
1341 	   assert(exprhdlr != NULL);
1342 	
1343 	   exprhdlrdata = SCIPexprhdlrGetData(exprhdlr);
1344 	   assert(exprhdlrdata != NULL);
1345 	
1346 	   base = SCIPexprGetChildren(expr)[0];
1347 	   assert(base != NULL);
1348 	
1349 	   exponent = SCIPgetExponentExprPow(expr);
1350 	
1351 	   SCIPdebugPrintf("[simplifyPow] simplifying power with expo %g\n", exponent);
1352 	
1353 	   /* enforces POW1 */
1354 	   if( exponent == 0.0 )
1355 	   {
1356 	      SCIPdebugPrintf("[simplifyPow] POW1\n");
1357 	      /* TODO: more checks? */
1358 	      assert(!SCIPisExprValue(scip, base) || SCIPgetValueExprValue(base) != 0.0);
1359 	      SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, 1.0, ownercreate, ownercreatedata) );
1360 	      return SCIP_OKAY;
1361 	   }
1362 	
1363 	   /* enforces POW2 */
1364 	   if( exponent == 1.0 )
1365 	   {
1366 	      SCIPdebugPrintf("[simplifyPow] POW2\n");
1367 	      *simplifiedexpr = base;
1368 	      SCIPcaptureExpr(*simplifiedexpr);
1369 	      return SCIP_OKAY;
1370 	   }
1371 	
1372 	   /* enforces POW3 */
1373 	   if( SCIPisExprValue(scip, base) )
1374 	   {
1375 	      SCIP_Real baseval;
1376 	
1377 	      SCIPdebugPrintf("[simplifyPow] POW3\n");
1378 	      baseval = SCIPgetValueExprValue(base);
1379 	
1380 	      /* the assert below was failing on st_e35 for baseval=-1e-15 and fractional exponent
1381 	       * in the subNLP heuristic; I assume that this was because baseval was evaluated after
1382 	       * variable fixings and that there were just floating-point inaccuracies and 0 was meant,
1383 	       * so I treat -1e-15 as 0 here
1384 	       */
1385 	      if( baseval < 0.0 && fmod(exponent, 1.0) != 0.0 && baseval > -SCIPepsilon(scip) )
1386 	         baseval = 0.0;
1387 	
1388 	      /* TODO check if those are all important asserts */
1389 	      assert(baseval >= 0.0 || fmod(exponent, 1.0) == 0.0);
1390 	      assert(baseval != 0.0 || exponent != 0.0);
1391 	
1392 	      if( baseval != 0.0 || exponent > 0.0 )
1393 	      {
1394 	         SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, pow(baseval, exponent), ownercreate, ownercreatedata) );
1395 	         return SCIP_OKAY;
1396 	      }
1397 	   }
1398 	
1399 	   /* enforces POW11 (exp(x)^n = exp(n*x)) */
1400 	   if( SCIPisExprExp(scip, base) )
1401 	   {
1402 	      SCIP_EXPR* child;
1403 	      SCIP_EXPR* prod;
1404 	      SCIP_EXPR* exponential;
1405 	      SCIP_EXPR* simplifiedprod;
1406 	
1407 	      SCIPdebugPrintf("[simplifyPow] POW11\n");
1408 	      child = SCIPexprGetChildren(base)[0];
1409 	
1410 	      /* multiply child of exponential with exponent */
1411 	      SCIP_CALL( SCIPcreateExprProduct(scip, &prod, 1, &child, exponent, ownercreate, ownercreatedata) );
1412 	
1413 	      /* simplify product */
1414 	      SCIP_CALL( SCIPcallExprSimplify(scip, prod, &simplifiedprod, ownercreate, ownercreatedata) );
1415 	      SCIP_CALL( SCIPreleaseExpr(scip, &prod) );
1416 	
1417 	      /* create exponential with new child */
1418 	      SCIP_CALL( SCIPcreateExprExp(scip, &exponential, simplifiedprod, ownercreate, ownercreatedata) );
1419 	      SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedprod) );
1420 	
1421 	      /* the final simplified expression is the simplification of the just created exponential */
1422 	      SCIP_CALL( SCIPcallExprSimplify(scip, exponential, simplifiedexpr, ownercreate, ownercreatedata) );
1423 	      SCIP_CALL( SCIPreleaseExpr(scip, &exponential) );
1424 	
1425 	      return SCIP_OKAY;
1426 	   }
1427 	
1428 	   /* enforces POW10 */
1429 	   if( SCIPisExprVar(scip, base) )
1430 	   {
1431 	      SCIP_VAR* basevar;
1432 	
1433 	      SCIPdebugPrintf("[simplifyPow] POW10\n");
1434 	      basevar = SCIPgetVarExprVar(base);
1435 	
1436 	      assert(basevar != NULL);
1437 	
1438 	      /* TODO: if exponent is negative, we could fix the binary variable to 1. However, this is a bit tricky because
1439 	       * variables can not be tighten in EXITPRE, where the simplify is also called
1440 	       */
1441 	      if( SCIPvarIsBinary(basevar) && exponent > 0.0 )
1442 	      {
1443 	         *simplifiedexpr = base;
1444 	         SCIPcaptureExpr(*simplifiedexpr);
1445 	         return SCIP_OKAY;
1446 	      }
1447 	   }
1448 	
1449 	   if( EPSISINT(exponent, 0.0) )
1450 	   {
1451 	      SCIP_EXPR* aux;
1452 	      SCIP_EXPR* simplifiedaux;
1453 	
1454 	      /* enforces POW12 (abs(x)^n = x^n if n is even) */
1455 	      if( SCIPisExprAbs(scip, base) && (int)exponent % 2 == 0 )
1456 	      {
1457 	         SCIP_EXPR* newpow;
1458 	
1459 	         SCIPdebugPrintf("[simplifyPow] POWXX\n");
1460 	
1461 	         SCIP_CALL( SCIPcreateExprPow(scip, &newpow, SCIPexprGetChildren(base)[0], exponent, ownercreate, ownercreatedata) );
1462 	         SCIP_CALL( simplifyPow(scip, newpow, simplifiedexpr, ownercreate, ownercreatedata) );
1463 	         SCIP_CALL( SCIPreleaseExpr(scip, &newpow) );
1464 	
1465 	         return SCIP_OKAY;
1466 	      }
1467 	
1468 	      /* enforces POW5
1469 	       * given (pow n (prod 1.0 expr_1 ... expr_k)) we distribute the exponent:
1470 	       * -> (prod 1.0 (pow n expr_1) ... (pow n expr_k))
1471 	       * notes: - since base is simplified and its coefficient is 1.0 (SP8)
1472 	       *        - n is an integer (excluding 1 and 0; see POW1-2 above)
1473 	       */
1474 	      if( SCIPisExprProduct(scip, base) )
1475 	      {
1476 	         SCIP_EXPR* auxproduct;
1477 	         int i;
1478 	
1479 	         /* create empty product */
1480 	         SCIP_CALL( SCIPcreateExprProduct(scip, &auxproduct, 0, NULL, 1.0, ownercreate, ownercreatedata) );
1481 	
1482 	         for( i = 0; i < SCIPexprGetNChildren(base); ++i )
1483 	         {
1484 	            /* create (pow n expr_i) and simplify */
1485 	            SCIP_CALL( SCIPcreateExprPow(scip, &aux, SCIPexprGetChildren(base)[i], exponent, ownercreate, ownercreatedata) );
1486 	            SCIP_CALL( simplifyPow(scip, aux, &simplifiedaux, ownercreate, ownercreatedata) );
1487 	            SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1488 	
1489 	            /* append (pow n expr_i) to product */
1490 	            SCIP_CALL( SCIPappendExprChild(scip, auxproduct, simplifiedaux) );
1491 	            SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedaux) );
1492 	         }
1493 	
1494 	         /* simplify (prod 1.0 (pow n expr_1) ... (pow n expr_k))
1495 	          * this calls simplifyProduct directly, since we know its children are simplified */
1496 	         SCIP_CALL( SCIPcallExprSimplify(scip, auxproduct, simplifiedexpr, ownercreate, ownercreatedata) );
1497 	         SCIP_CALL( SCIPreleaseExpr(scip, &auxproduct) );
1498 	         return SCIP_OKAY;
1499 	      }
1500 	
1501 	      /* enforces POW6
1502 	       * given (pow n (sum 0.0 coef expr)) we can move `pow` inside `sum`:
1503 	       * (pow n (sum 0.0 coef expr) ) -> (sum 0.0 coef^n (pow n expr))
1504 	       * notes: - since base is simplified and its constant is 0, then coef != 1.0 (SS7)
1505 	       *        - n is an integer (excluding 1 and 0; see POW1-2 above)
1506 	       */
1507 	      if( SCIPisExprSum(scip, base) && SCIPexprGetNChildren(base) == 1 && SCIPgetConstantExprSum(base) == 0.0 )
1508 	      {
1509 	         SCIP_Real newcoef;
1510 	
1511 	         SCIPdebugPrintf("[simplifyPow] seeing a sum with one term, exponent %g\n", exponent);
1512 	
1513 	         /* assert SS7 holds */
1514 	         assert(SCIPgetCoefsExprSum(base)[0] != 1.0);
1515 	
1516 	         /* create (pow n expr) and simplify it
1517 	          * note: we call simplifyPow directly, since we know that `expr` is simplified */
1518 	         newcoef = pow(SCIPgetCoefsExprSum(base)[0], exponent);
1519 	         SCIP_CALL( SCIPcreateExprPow(scip, &aux, SCIPexprGetChildren(base)[0], exponent, ownercreate, ownercreatedata) );
1520 	         SCIP_CALL( simplifyPow(scip, aux, &simplifiedaux, ownercreate, ownercreatedata) );
1521 	         SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1522 	
1523 	         /* create (sum (pow n expr)) and simplify it
1524 	          * this calls simplifySum directly, since we know its children are simplified */
1525 	         SCIP_CALL( SCIPcreateExprSum(scip, &aux, 1, &simplifiedaux, &newcoef, 0.0, ownercreate, ownercreatedata) );
1526 	         SCIP_CALL( SCIPcallExprSimplify(scip, aux, simplifiedexpr, ownercreate, ownercreatedata) );
1527 	         SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1528 	         SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedaux) );
1529 	         return SCIP_OKAY;
1530 	      }
1531 	
1532 	      /* enforces POW7 for exponent 2
1533 	       * (const + sum alpha_i expr_i)^2 = sum alpha_i^2 expr_i^2
1534 	       * + sum_{j < i} 2*alpha_i alpha_j expr_i expr_j
1535 	       * + sum const alpha_i expr_i
1536 	       * TODO: put some limits on the number of children of the sum being expanded
1537 	       */
1538 	      if( SCIPisExprSum(scip, base) && exponent == 2.0 && exprhdlrdata->expandmaxexponent >= 2 )
1539 	      {
1540 	         int i;
1541 	         int nchildren;
1542 	         int nexpandedchildren;
1543 	         SCIP_EXPR* expansion;
1544 	         SCIP_EXPR** expandedchildren;
1545 	         SCIP_Real* coefs;
1546 	         SCIP_Real constant;
1547 	
1548 	         SCIPdebugPrintf("[simplifyPow] expanding sum^%g\n", exponent);
1549 	
1550 	         nchildren = SCIPexprGetNChildren(base);
1551 	         nexpandedchildren = nchildren * (nchildren + 1) / 2 + nchildren;
1552 	         SCIP_CALL( SCIPallocBufferArray(scip, &coefs, nexpandedchildren) );
1553 	         SCIP_CALL( SCIPallocBufferArray(scip, &expandedchildren, nexpandedchildren) );
1554 	
1555 	         for( i = 0; i < nchildren; ++i )
1556 	         {
1557 	            int j;
1558 	            SCIP_EXPR* expansionchild;
1559 	            SCIP_EXPR* prodchildren[2];
1560 	            prodchildren[0] = SCIPexprGetChildren(base)[i];
1561 	
1562 	            /* create and simplify expr_i * expr_j */
1563 	            for( j = 0; j < i; ++j )
1564 	            {
1565 	               prodchildren[1] = SCIPexprGetChildren(base)[j];
1566 	               coefs[i*(i+1)/2 + j] = 2 * SCIPgetCoefsExprSum(base)[i] * SCIPgetCoefsExprSum(base)[j];
1567 	
1568 	               SCIP_CALL( SCIPcreateExprProduct(scip, &expansionchild, 2, prodchildren, 1.0, ownercreate,
1569 	                          ownercreatedata) );
1570 	               SCIP_CALL( SCIPcallExprSimplify(scip, expansionchild, &expandedchildren[i*(i+1)/2 + j],
1571 	                          ownercreate, ownercreatedata) ); /* this calls simplifyProduct */
1572 	               SCIP_CALL( SCIPreleaseExpr(scip, &expansionchild) );
1573 	            }
1574 	            /* create and simplify expr_i * expr_i */
1575 	            prodchildren[1] = SCIPexprGetChildren(base)[i];
1576 	            coefs[i*(i+1)/2 + i] = SCIPgetCoefsExprSum(base)[i] * SCIPgetCoefsExprSum(base)[i];
1577 	
1578 	            SCIP_CALL( SCIPcreateExprProduct(scip, &expansionchild, 2, prodchildren, 1.0, ownercreate,
1579 	                       ownercreatedata) );
1580 	            SCIP_CALL( SCIPcallExprSimplify(scip, expansionchild, &expandedchildren[i*(i+1)/2 + i], ownercreate,
1581 	                       ownercreatedata) ); /* this calls simplifyProduct */
1582 	            SCIP_CALL( SCIPreleaseExpr(scip, &expansionchild) );
1583 	         }
1584 	         /* create const * alpha_i expr_i */
1585 	         for( i = 0; i < nchildren; ++i )
1586 	         {
1587 	            coefs[i + nexpandedchildren - nchildren] = 2 * SCIPgetConstantExprSum(base) * SCIPgetCoefsExprSum(base)[i];
1588 	            expandedchildren[i + nexpandedchildren - nchildren] = SCIPexprGetChildren(base)[i];
1589 	         }
1590 	
1591 	         constant = SCIPgetConstantExprSum(base);
1592 	         constant *= constant;
1593 	         /* create sum of all the above and simplify it with simplifySum since all of its children are simplified! */
1594 	         SCIP_CALL( SCIPcreateExprSum(scip, &expansion, nexpandedchildren, expandedchildren, coefs, constant,
1595 	                    ownercreate, ownercreatedata) );
1596 	         SCIP_CALL( SCIPcallExprSimplify(scip, expansion, simplifiedexpr, ownercreate,
1597 	                    ownercreatedata) ); /* this calls simplifySum */
1598 	
1599 	         /* release everything */
1600 	         SCIP_CALL( SCIPreleaseExpr(scip, &expansion) );
1601 	         /* release the *created* expanded children */
1602 	         for( i = 0; i < nexpandedchildren - nchildren; ++i )
1603 	         {
1604 	            SCIP_CALL( SCIPreleaseExpr(scip, &expandedchildren[i]) );
1605 	         }
1606 	         SCIPfreeBufferArray(scip, &expandedchildren);
1607 	         SCIPfreeBufferArray(scip, &coefs);
1608 	
1609 	         return SCIP_OKAY;
1610 	      }
1611 	
1612 	      /* enforces POW7 for exponent > 2 */
1613 	      if( SCIPisExprSum(scip, base) && exponent > 2.0 && exponent <= exprhdlrdata->expandmaxexponent )
1614 	      {
1615 	         SCIPdebugPrintf("[simplifyPow] expanding sum^%g\n", exponent);
1616 	
1617 	         SCIP_CALL( SCIPpowerExprSum(scip, simplifiedexpr, base, (int)exponent, TRUE, ownercreate, ownercreatedata) );
1618 	
1619 	         return SCIP_OKAY;
1620 	      }
1621 	
1622 	   }
1623 	   else
1624 	   {
1625 	      /* enforces POW9
1626 	       *
1627 	       * FIXME code of POW6 is very similar
1628 	       */
1629 	      if( SCIPexprGetNChildren(base) == 1
1630 	         && SCIPisExprSum(scip, base)
1631 	         && SCIPgetConstantExprSum(base) == 0.0
1632 	         && SCIPgetCoefsExprSum(base)[0] >= 0.0 )
1633 	      {
1634 	         SCIP_EXPR* simplifiedaux;
1635 	         SCIP_EXPR* aux;
1636 	         SCIP_Real newcoef;
1637 	
1638 	         SCIPdebugPrintf("[simplifyPow] seeing a sum with one term, exponent %g\n", exponent);
1639 	         /* assert SS7 holds */
1640 	         assert(SCIPgetCoefsExprSum(base)[0] != 1.0);
1641 	
1642 	         /* create (pow n expr) and simplify it
1643 	          * note: we call simplifyPow directly, since we know that `expr` is simplified */
1644 	         SCIP_CALL( SCIPcreateExprPow(scip, &aux, SCIPexprGetChildren(base)[0], exponent, ownercreate,
1645 	                    ownercreatedata) );
1646 	         SCIP_CALL( simplifyPow(scip, aux, &simplifiedaux, ownercreate, ownercreatedata) );
1647 	         SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1648 	
1649 	         /* create (sum (pow n expr)) and simplify it
1650 	          * this calls simplifySum directly, since we know its child is simplified! */
1651 	         newcoef = pow(SCIPgetCoefsExprSum(base)[0], exponent);
1652 	         SCIP_CALL( SCIPcreateExprSum(scip, &aux, 1, &simplifiedaux, &newcoef, 0.0, ownercreate,
1653 	                    ownercreatedata) );
1654 	         SCIP_CALL( SCIPcallExprSimplify(scip, aux, simplifiedexpr, ownercreate, ownercreatedata) );
1655 	         SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1656 	         SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedaux) );
1657 	
1658 	         return SCIP_OKAY;
1659 	      }
1660 	
1661 	      /* enforces POW5a
1662 	       * given (pow n (prod 1.0 expr_1 ... expr_k)) we distribute the exponent:
1663 	       * -> (prod 1.0 (pow n expr_1) ... (pow n expr_k))
1664 	       * notes: - since base is simplified and its coefficient is 1.0 (SP8)
1665 	       * TODO we can enable this more often by default when simplify makes use of bounds on factors
1666 	       */
1667 	      if( exprhdlrdata->distribfracexponent && SCIPisExprProduct(scip, base) )
1668 	      {
1669 	         SCIP_EXPR* aux;
1670 	         SCIP_EXPR* simplifiedaux;
1671 	         SCIP_EXPR* auxproduct;
1672 	         int i;
1673 	
1674 	         /* create empty product */
1675 	         SCIP_CALL( SCIPcreateExprProduct(scip, &auxproduct, 0, NULL, 1.0, ownercreate, ownercreatedata) );
1676 	
1677 	         for( i = 0; i < SCIPexprGetNChildren(base); ++i )
1678 	         {
1679 	            /* create (pow n expr_i) and simplify */
1680 	            SCIP_CALL( SCIPcreateExprPow(scip, &aux, SCIPexprGetChildren(base)[i], exponent, ownercreate, ownercreatedata) );
1681 	            SCIP_CALL( simplifyPow(scip, aux, &simplifiedaux, ownercreate, ownercreatedata) );
1682 	            SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1683 	
1684 	            /* append (pow n expr_i) to product */
1685 	            SCIP_CALL( SCIPappendExprChild(scip, auxproduct, simplifiedaux) );
1686 	            SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedaux) );
1687 	         }
1688 	
1689 	         /* simplify (prod 1.0 (pow n expr_1) ... (pow n expr_k))
1690 	          * this calls simplifyProduct directly, since we know its children are simplified */
1691 	         SCIP_CALL( SCIPcallExprSimplify(scip, auxproduct, simplifiedexpr, ownercreate, ownercreatedata) );
1692 	         SCIP_CALL( SCIPreleaseExpr(scip, &auxproduct) );
1693 	         return SCIP_OKAY;
1694 	      }
1695 	   }
1696 	
1697 	   /* enforces POW8
1698 	    * given (pow n (pow expo expr)) we distribute the exponent:
1699 	    * -> (pow n*expo expr)
1700 	    * notes: n is not 1 or 0, see POW1-2 above
1701 	    */
1702 	   if( SCIPisExprPower(scip, base) )
1703 	   {
1704 	      SCIP_Real newexponent;
1705 	      SCIP_Real baseexponent;
1706 	
1707 	      baseexponent = SCIPgetExponentExprPow(base);
1708 	      newexponent = baseexponent * exponent;
1709 	
1710 	      /* some checks (see POW8 definition in scip_expr.h) to make sure we don't loose an
1711 	       * implicit SCIPexprGetChildren(base)[0] >= 0 constraint
1712 	       *
1713 	       * if newexponent is fractional, then we will still need expr >= 0
1714 	       * if both exponents were integer, then we never required and will not require expr >= 0
1715 	       * if base exponent was an even integer, then we did not require expr >= 0
1716 	       * (but may need to use |expr|^newexponent)
1717 	       */
1718 	      if( !EPSISINT(newexponent, 0.0) ||
1719 	          (EPSISINT(baseexponent, 0.0) && EPSISINT(exponent, 0.0)) ||
1720 	          (EPSISINT(baseexponent, 0.0) && ((int)baseexponent) % 2 == 0) )
1721 	      {
1722 	         SCIP_EXPR* aux;
1723 	
1724 	         if( EPSISINT(baseexponent, 0.0) && ((int)baseexponent) % 2 == 0 &&
1725 	             (!EPSISINT(newexponent, 0.0) || ((int)newexponent) % 2 == 1) )
1726 	         {
1727 	            /* If base exponent was even integer and new exponent will be fractional,
1728 	             * then simplify to |expr|^newexponent to allow eval for expr < 0.
1729 	             * If base exponent was even integer and new exponent will be odd integer,
1730 	             * then simplify to |expr|^newexponent to preserve value for expr < 0.
1731 	             */
1732 	            SCIP_EXPR* simplifiedaux;
1733 	
1734 	            SCIP_CALL( SCIPcreateExprAbs(scip, &aux, SCIPexprGetChildren(base)[0], ownercreate, ownercreatedata) );
1735 	            SCIP_CALL( SCIPcallExprSimplify(scip, aux, &simplifiedaux, ownercreate, ownercreatedata) );
1736 	            SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1737 	            SCIP_CALL( SCIPcreateExprPow(scip, &aux, simplifiedaux, newexponent, ownercreate, ownercreatedata) );
1738 	            SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedaux) );
1739 	         }
1740 	         else
1741 	         {
1742 	            SCIP_CALL( SCIPcreateExprPow(scip, &aux, SCIPexprGetChildren(base)[0], newexponent, ownercreate,
1743 	                       ownercreatedata) );
1744 	         }
1745 	
1746 	         SCIP_CALL( simplifyPow(scip, aux, simplifiedexpr, ownercreate, ownercreatedata) );
1747 	         SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1748 	
1749 	         return SCIP_OKAY;
1750 	      }
1751 	   }
1752 	
1753 	   SCIPdebugPrintf("[simplifyPow] power is simplified\n");
1754 	   *simplifiedexpr = expr;
1755 	
1756 	   /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */
1757 	   SCIPcaptureExpr(*simplifiedexpr);
1758 	
1759 	   return SCIP_OKAY;
1760 	}
1761 	
1762 	/** expression callback to get information for symmetry detection */
1763 	static
1764 	SCIP_DECL_EXPRGETSYMDATA(getSymDataPow)
1765 	{  /*lint --e{715}*/
1766 	   SCIP_EXPRDATA* exprdata;
1767 	
1768 	   assert(scip != NULL);
1769 	   assert(expr != NULL);
1770 	
1771 	   exprdata = SCIPexprGetData(expr);
1772 	   assert(exprdata != NULL);
1773 	
1774 	   SCIP_CALL( SCIPallocBlockMemory(scip, symdata) );
1775 	
1776 	   (*symdata)->nconstants = 1;
1777 	   (*symdata)->ncoefficients = 0;
1778 	
1779 	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*symdata)->constants, 1) );
1780 	   (*symdata)->constants[0] = exprdata->exponent;
1781 	
1782 	   return SCIP_OKAY;
1783 	}
1784 	
1785 	/** expression handler copy callback */
1786 	static
1787 	SCIP_DECL_EXPRCOPYHDLR(copyhdlrPow)
1788 	{  /*lint --e{715}*/
1789 	   SCIP_CALL( SCIPincludeExprhdlrPow(scip) );
1790 	
1791 	   return SCIP_OKAY;
1792 	}
1793 	
1794 	/** expression handler free callback */
1795 	static
1796 	SCIP_DECL_EXPRFREEHDLR(freehdlrPow)
1797 	{  /*lint --e{715}*/
1798 	   assert(exprhdlrdata != NULL);
1799 	   assert(*exprhdlrdata != NULL);
1800 	
1801 	   SCIPfreeBlockMemory(scip, exprhdlrdata);
1802 	
1803 	   return SCIP_OKAY;
1804 	}
1805 	
1806 	/** expression data copy callback */
1807 	static
1808 	SCIP_DECL_EXPRCOPYDATA(copydataPow)
1809 	{  /*lint --e{715}*/
1810 	   SCIP_EXPRDATA* sourceexprdata;
1811 	
1812 	   assert(targetexprdata != NULL);
1813 	   assert(sourceexpr != NULL);
1814 	
1815 	   sourceexprdata = SCIPexprGetData(sourceexpr);
1816 	   assert(sourceexprdata != NULL);
1817 	
1818 	   *targetexprdata = NULL;
1819 	
1820 	   SCIP_CALL( createData(targetscip, targetexprdata, sourceexprdata->exponent) );
1821 	
1822 	   return SCIP_OKAY;
1823 	}
1824 	
1825 	/** expression data free callback */
1826 	static
1827 	SCIP_DECL_EXPRFREEDATA(freedataPow)
1828 	{  /*lint --e{715}*/
1829 	   SCIP_EXPRDATA* exprdata;
1830 	
1831 	   assert(expr != NULL);
1832 	
1833 	   exprdata = SCIPexprGetData(expr);
1834 	   assert(exprdata != NULL);
1835 	
1836 	   SCIPfreeBlockMemory(scip, &exprdata);
1837 	   SCIPexprSetData(expr, NULL);
1838 	
1839 	   return SCIP_OKAY;
1840 	}
1841 	
1842 	/** expression print callback */
1843 	/** @todo: use precedence for better printing */
1844 	static
1845 	SCIP_DECL_EXPRPRINT(printPow)
1846 	{  /*lint --e{715}*/
1847 	   assert(expr != NULL);
1848 	
1849 	   /**! [SnippetExprPrintPow] */
1850 	   switch( stage )
1851 	   {
1852 	      case SCIP_EXPRITER_ENTEREXPR :
1853 	      {
1854 	         /* print function with opening parenthesis */
1855 	         SCIPinfoMessage(scip, file, "(");
1856 	         break;
1857 	      }
1858 	
1859 	      case SCIP_EXPRITER_VISITINGCHILD :
1860 	      {
1861 	         assert(currentchild == 0);
1862 	         break;
1863 	      }
1864 	
1865 	      case SCIP_EXPRITER_LEAVEEXPR :
1866 	      {
1867 	         SCIP_Real exponent = SCIPgetExponentExprPow(expr);
1868 	
1869 	         /* print closing parenthesis */
1870 	         if( exponent >= 0.0 )
1871 	            SCIPinfoMessage(scip, file, ")^%g", exponent);
1872 	         else
1873 	            SCIPinfoMessage(scip, file, ")^(%g)", exponent);
1874 	
1875 	         break;
1876 	      }
1877 	
1878 	      case SCIP_EXPRITER_VISITEDCHILD :
1879 	      default:
1880 	         break;
1881 	   }
1882 	   /**! [SnippetExprPrintPow] */
1883 	
1884 	   return SCIP_OKAY;
1885 	}
1886 	
1887 	/** expression point evaluation callback */
1888 	static
1889 	SCIP_DECL_EXPREVAL(evalPow)
1890 	{  /*lint --e{715}*/
1891 	   SCIP_Real exponent;
1892 	   SCIP_Real base;
1893 	
1894 	   assert(expr != NULL);
1895 	   assert(SCIPexprGetNChildren(expr) == 1);
1896 	   assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]) != SCIP_INVALID);
1897 	
1898 	   exponent = SCIPgetExponentExprPow(expr);
1899 	   base = SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]);
1900 	
1901 	   *val = pow(base, exponent);
1902 	
1903 	   /* if there is a domain, pole, or range error, pow() should return some kind of NaN, infinity, or HUGE_VAL
1904 	    * we could also work with floating point exceptions or errno, but I am not sure this would be thread-safe
1905 	    */
1906 	   if( !SCIPisFinite(*val) || *val == HUGE_VAL || *val == -HUGE_VAL )
1907 	      *val = SCIP_INVALID;
1908 	
1909 	   return SCIP_OKAY;
1910 	}
1911 	
1912 	/** derivative evaluation callback
1913 	 *
1914 	 * computes <gradient, children.dot>
1915 	 * if expr is child^p, then computes
1916 	 * p child^(p-1) dot(child)
1917 	 */
1918 	static
1919 	SCIP_DECL_EXPRFWDIFF(fwdiffPow)
1920 	{  /*lint --e{715}*/
1921 	   SCIP_EXPR* child;
1922 	   SCIP_Real exponent;
1923 	
1924 	   assert(expr != NULL);
1925 	   assert(SCIPexprGetData(expr) != NULL);
1926 	   assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID);
1927 	
1928 	   child = SCIPexprGetChildren(expr)[0];
1929 	   assert(child != NULL);
1930 	   assert(!SCIPisExprValue(scip, child));
1931 	   assert(SCIPexprGetDot(child) != SCIP_INVALID);
1932 	
1933 	   exponent = SCIPgetExponentExprPow(expr);
1934 	   assert(exponent != 0.0);
1935 	
1936 	   /* x^exponent is not differentiable for x = 0 and exponent in ]0,1[ */
1937 	   if( exponent > 0.0 && exponent < 1.0 && SCIPexprGetEvalValue(child) == 0.0 )
1938 	      *dot = SCIP_INVALID;
1939 	   else
1940 	      *dot = exponent * pow(SCIPexprGetEvalValue(child), exponent - 1.0) * SCIPexprGetDot(child);
1941 	
1942 	   return SCIP_OKAY;
1943 	}
1944 	
1945 	/** expression backward forward derivative evaluation callback
1946 	 *
1947 	 * computes partial/partial child ( <gradient, children.dot> )
1948 	 * if expr is child^n, then computes
1949 	 * n * (n - 1) child^(n-2) dot(child)
1950 	 */
1951 	static
1952 	SCIP_DECL_EXPRBWFWDIFF(bwfwdiffPow)
1953 	{  /*lint --e{715}*/
1954 	   SCIP_EXPR* child;
1955 	   SCIP_Real exponent;
1956 	
1957 	   assert(expr != NULL);
1958 	   assert(SCIPexprGetData(expr) != NULL);
1959 	   assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID);
1960 	   assert(childidx == 0);
1961 	
1962 	   child = SCIPexprGetChildren(expr)[0];
1963 	   assert(child != NULL);
1964 	   assert(!SCIPisExprValue(scip, child));
1965 	   assert(SCIPexprGetDot(child) != SCIP_INVALID);
1966 	
1967 	   exponent = SCIPgetExponentExprPow(expr);
1968 	   assert(exponent != 0.0);
1969 	
1970 	   /* x^exponent is not twice differentiable for x = 0 and exponent in ]0,1[ u ]1,2[ */
1971 	   if( exponent > 0.0 && exponent < 2.0 && SCIPexprGetEvalValue(child) == 0.0 && exponent != 1.0 )
1972 	      *bardot = SCIP_INVALID;
1973 	   else
1974 	      *bardot = exponent * (exponent - 1.0) * pow(SCIPexprGetEvalValue(child), exponent - 2.0) * SCIPexprGetDot(child);
1975 	
1976 	   return SCIP_OKAY;
1977 	}
1978 	
1979 	/** expression derivative evaluation callback */
1980 	static
1981 	SCIP_DECL_EXPRBWDIFF(bwdiffPow)
1982 	{  /*lint --e{715}*/
1983 	   SCIP_EXPR* child;
1984 	   SCIP_Real childval;
1985 	   SCIP_Real exponent;
1986 	
1987 	   assert(expr != NULL);
1988 	   assert(SCIPexprGetData(expr) != NULL);
1989 	   assert(childidx == 0);
1990 	
1991 	   child = SCIPexprGetChildren(expr)[0];
1992 	   assert(child != NULL);
1993 	   assert(!SCIPisExprValue(scip, child));
1994 	
1995 	   childval = SCIPexprGetEvalValue(child);
1996 	   assert(childval != SCIP_INVALID);
1997 	
1998 	   exponent = SCIPgetExponentExprPow(expr);
1999 	   assert(exponent != 0.0);
2000 	
2001 	   /* x^exponent is not differentiable for x = 0 and exponent in ]0,1[ */
2002 	   if( exponent > 0.0 && exponent < 1.0 && childval == 0.0 )
2003 	      *val = SCIP_INVALID;
2004 	   else
2005 	      *val = exponent * pow(childval, exponent - 1.0);
2006 	
2007 	   return SCIP_OKAY;
2008 	}
2009 	
2010 	/** expression interval evaluation callback */
2011 	static
2012 	SCIP_DECL_EXPRINTEVAL(intevalPow)
2013 	{  /*lint --e{715}*/
2014 	   SCIP_INTERVAL childinterval;
2015 	   SCIP_Real exponent;
2016 	
2017 	   assert(expr != NULL);
2018 	   assert(SCIPexprGetNChildren(expr) == 1);
2019 	
2020 	   childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]);
2021 	
2022 	   exponent = SCIPgetExponentExprPow(expr);
2023 	
2024 	   if( exponent < 0.0 )
2025 	   {
2026 	      SCIP_EXPRHDLRDATA* exprhdlrdata;
2027 	      exprhdlrdata = SCIPexprhdlrGetData(SCIPexprGetHdlr(expr));
2028 	      assert(exprhdlrdata != NULL);
2029 	
2030 	      if( exprhdlrdata->minzerodistance > 0.0 )
2031 	      {
2032 	         /* avoid small interval around 0 if possible, see also reversepropPow */
2033 	         if( childinterval.inf > -exprhdlrdata->minzerodistance && childinterval.inf < exprhdlrdata->minzerodistance )
2034 	         {
2035 	            if( !exprhdlrdata->warnedonpole && SCIPgetVerbLevel(scip) > SCIP_VERBLEVEL_NONE )
2036 	            {
2037 	               SCIPinfoMessage(scip, NULL, "Changing lower bound for child of pow(.,%g) from %g to %g.\n"
2038 	                  "Check your model formulation or use option expr/" POWEXPRHDLR_NAME "/minzerodistance to avoid this warning.\n",
2039 	                  exponent, childinterval.inf, exprhdlrdata->minzerodistance);
2040 	               SCIPinfoMessage(scip, NULL, "Expression: ");
2041 	               SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
2042 	               SCIPinfoMessage(scip, NULL, "\n");
2043 	               exprhdlrdata->warnedonpole = TRUE;
2044 	            }
2045 	            childinterval.inf = exprhdlrdata->minzerodistance;
2046 	         }
2047 	         else if( childinterval.sup < exprhdlrdata->minzerodistance
2048 	                  && childinterval.sup > -exprhdlrdata->minzerodistance )
2049 	         {
2050 	            if( !exprhdlrdata->warnedonpole && SCIPgetVerbLevel(scip) > SCIP_VERBLEVEL_NONE )
2051 	            {
2052 	               SCIPinfoMessage(scip, NULL, "Changing upper bound for child of pow(.,%g) from %g to %g.\n"
2053 	                  "Check your model formulation or use option expr/" POWEXPRHDLR_NAME "/minzerodistance to avoid this warning.\n",
2054 	                  exponent, childinterval.sup, -exprhdlrdata->minzerodistance);
2055 	               SCIPinfoMessage(scip, NULL, "Expression: ");
2056 	               SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
2057 	               SCIPinfoMessage(scip, NULL, "\n");
2058 	               exprhdlrdata->warnedonpole = TRUE;
2059 	            }
2060 	            childinterval.sup = -exprhdlrdata->minzerodistance;
2061 	         }
2062 	      }
2063 	   }
2064 	
2065 	   if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) )
2066 	   {
2067 	      SCIPintervalSetEmpty(interval);
2068 	      return SCIP_OKAY;
2069 	   }
2070 	
2071 	   SCIPintervalPowerScalar(SCIP_INTERVAL_INFINITY, interval, childinterval, exponent);
2072 	
2073 	   /* make sure 0^negative is an empty interval, as some other codes do not handle intervals like [inf,inf] well
2074 	    * TODO maybe change SCIPintervalPowerScalar?
2075 	    */
2076 	   if( exponent < 0.0 && childinterval.inf == 0.0 && childinterval.sup == 0.0 )
2077 	      SCIPintervalSetEmpty(interval);
2078 	
2079 	   return SCIP_OKAY;
2080 	}
2081 	
2082 	/** expression estimator callback */
2083 	static
2084 	SCIP_DECL_EXPRESTIMATE(estimatePow)
2085 	{  /*lint --e{715}*/
2086 	   SCIP_EXPRDATA* exprdata;
2087 	   SCIP_EXPR* child;
2088 	   SCIP_Real childlb;
2089 	   SCIP_Real childub;
2090 	   SCIP_Real exponent;
2091 	   SCIP_Bool isinteger;
2092 	
2093 	   assert(scip != NULL);
2094 	   assert(expr != NULL);
2095 	   assert(SCIPexprGetNChildren(expr) == 1);
2096 	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), POWEXPRHDLR_NAME) == 0);
2097 	   assert(refpoint != NULL);
2098 	   assert(coefs != NULL);
2099 	   assert(constant != NULL);
2100 	   assert(islocal != NULL);
2101 	   assert(branchcand != NULL);
2102 	   assert(*branchcand == TRUE);  /* the default */
2103 	   assert(success != NULL);
2104 	
2105 	   *success = FALSE;
2106 	
2107 	   /* get aux variables: we over- or underestimate childvar^exponent  */
2108 	   child = SCIPexprGetChildren(expr)[0];
2109 	   assert(child != NULL);
2110 	
2111 	   SCIPdebugMsg(scip, "%sestimation of x^%g at x=%.15g\n",
2112 	      overestimate ? "over" : "under", SCIPgetExponentExprPow(expr), *refpoint);
2113 	
2114 	   /* we can not generate a cut at +/- infinity */
2115 	   if( SCIPisInfinity(scip, REALABS(*refpoint)) )
2116 	      return SCIP_OKAY;
2117 	
2118 	   childlb = localbounds[0].inf;
2119 	   childub = localbounds[0].sup;
2120 	
2121 	   exprdata = SCIPexprGetData(expr);
2122 	   exponent = exprdata->exponent;
2123 	   assert(exponent != 1.0 && exponent != 0.0); /* this should have been simplified */
2124 	
2125 	   /* if child is constant, then return a constant estimator
2126 	    * this can help with small infeasibilities if boundtightening is relaxing bounds too much
2127 	    */
2128 	   if( childlb == childub )
2129 	   {
2130 	      *coefs = 0.0;
2131 	      *constant = pow(childlb, exponent);
2132 	      *success = TRUE;
2133 	      *islocal = globalbounds[0].inf != globalbounds[0].sup;
2134 	      *branchcand = FALSE;
2135 	      return SCIP_OKAY;
2136 	   }
2137 	
2138 	   isinteger = EPSISINT(exponent, 0.0);
2139 	
2140 	   /* if exponent is not integral, then child must be non-negative */
2141 	   if( !isinteger && childlb < 0.0 )
2142 	   {
2143 	      /* somewhere we should have tightened the bound on x, but small tightening are not always applied by SCIP
2144 	       * it is ok to do this tightening here, but let's assert that we were close to 0.0 already
2145 	       */
2146 	      assert(SCIPisFeasZero(scip, childlb));
2147 	      childlb = 0.0;
2148 	   }
2149 	   assert(isinteger || childlb >= 0.0);
2150 	
2151 	   SCIP_CALL( buildPowEstimator(scip, exprdata, overestimate, childlb, childub, globalbounds[0].inf,
2152 	         globalbounds[0].sup, SCIPexprIsIntegral(child), MAX(childlb, *refpoint), exponent, coefs,
2153 	         constant, success, islocal, branchcand) );
2154 	
2155 	   return SCIP_OKAY;
2156 	}
2157 	
2158 	/** expression reverse propagaton callback */
2159 	static
2160 	SCIP_DECL_EXPRREVERSEPROP(reversepropPow)
2161 	{  /*lint --e{715}*/
2162 	   SCIP_INTERVAL child;
2163 	   SCIP_INTERVAL interval;
2164 	   SCIP_Real exponent;
2165 	
2166 	   assert(scip != NULL);
2167 	   assert(expr != NULL);
2168 	   assert(SCIPexprGetNChildren(expr) == 1);
2169 	
2170 	   exponent = SCIPgetExponentExprPow(expr);
2171 	   child = childrenbounds[0];
2172 	
2173 	   SCIPdebugMsg(scip, "reverseprop x^%g in [%.15g,%.15g], x = [%.15g,%.15g]", exponent, bounds.inf, bounds.sup,
2174 	         child.inf, child.sup);
2175 	
2176 	   if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, child) )
2177 	   {
2178 	      *infeasible = TRUE;
2179 	      return SCIP_OKAY;
2180 	   }
2181 	
2182 	   if( SCIPintervalIsEntire(SCIP_INTERVAL_INFINITY, bounds) )
2183 	   {
2184 	      /* if exponent is not integral, then make sure that child is non-negative */
2185 	      if( !EPSISINT(exponent, 0.0) && child.inf < 0.0 )
2186 	      {
2187 	         SCIPintervalSetBounds(&interval, 0.0, child.sup);
2188 	      }
2189 	      else
2190 	      {
2191 	         SCIPdebugMsgPrint(scip, "-> no improvement\n");
2192 	         return SCIP_OKAY;
2193 	      }
2194 	   }
2195 	   else
2196 	   {
2197 	      /* f = pow(c0, alpha) -> c0 = pow(f, 1/alpha) */
2198 	      SCIPintervalPowerScalarInverse(SCIP_INTERVAL_INFINITY, &interval, child, exponent, bounds);
2199 	   }
2200 	
2201 	   if( exponent < 0.0 )
2202 	   {
2203 	      SCIP_EXPRHDLRDATA* exprhdlrdata;
2204 	
2205 	      exprhdlrdata = SCIPexprhdlrGetData(SCIPexprGetHdlr(expr));
2206 	      assert(exprhdlrdata != NULL);
2207 	
2208 	      if( exprhdlrdata->minzerodistance > 0.0 )
2209 	      {
2210 	         /* push lower bound from >= -epsilon to >=  epsilon to avoid pole at 0 (domain error)
2211 	          * push upper bound from <=  epsilon to <= -epsilon to avoid pole at 0 (domain error)
2212 	          * this can lead to a cutoff if domain would otherwise be very close around 0
2213 	          */
2214 	         if( interval.inf > -exprhdlrdata->minzerodistance && interval.inf < exprhdlrdata->minzerodistance )
2215 	         {
2216 	            if( !exprhdlrdata->warnedonpole && SCIPgetVerbLevel(scip) > SCIP_VERBLEVEL_NONE )
2217 	            {
2218 	               SCIPinfoMessage(scip, NULL, "Changing lower bound for child of pow(.,%g) from %g to %g.\n"
2219 	                  "Check your model formulation or use option expr/" POWEXPRHDLR_NAME "/minzerodistance to avoid this warning.\n",
2220 	                  exponent, interval.inf, exprhdlrdata->minzerodistance);
2221 	               SCIPinfoMessage(scip, NULL, "Expression: ");
2222 	               SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
2223 	               SCIPinfoMessage(scip, NULL, "\n");
2224 	               exprhdlrdata->warnedonpole = TRUE;
2225 	            }
2226 	            interval.inf = exprhdlrdata->minzerodistance;
2227 	         }
2228 	         else if( interval.sup < exprhdlrdata->minzerodistance && interval.sup > -exprhdlrdata->minzerodistance )
2229 	         {
2230 	            if( !exprhdlrdata->warnedonpole && SCIPgetVerbLevel(scip) > SCIP_VERBLEVEL_NONE )
2231 	            {
2232 	               SCIPinfoMessage(scip, NULL, "Changing lower bound for child of pow(.,%g) from %g to %g.\n"
2233 	                  "Check your model formulation or use option expr/" POWEXPRHDLR_NAME "/minzerodistance to avoid this warning.\n",
2234 	                  exponent, interval.sup, -exprhdlrdata->minzerodistance);
2235 	               SCIPinfoMessage(scip, NULL, "Expression: ");
2236 	               SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
2237 	               SCIPinfoMessage(scip, NULL, "\n");
2238 	               exprhdlrdata->warnedonpole = TRUE;
2239 	            }
2240 	            interval.sup = -exprhdlrdata->minzerodistance;
2241 	         }
2242 	      }
2243 	   }
2244 	
2245 	   SCIPdebugMsgPrint(scip, " -> [%.15g,%.15g]\n", interval.inf, interval.sup);
2246 	
2247 	   childrenbounds[0] = interval;
2248 	
2249 	   return SCIP_OKAY;
2250 	}
2251 	
2252 	/** initial estimates callback for a power expression */
2253 	static
2254 	SCIP_DECL_EXPRINITESTIMATES(initestimatesPow)
2255 	{
2256 	   SCIP_EXPRDATA* exprdata;
2257 	   SCIP_EXPR* child;
2258 	   SCIP_Real childlb;
2259 	   SCIP_Real childub;
2260 	   SCIP_Real exponent;
2261 	   SCIP_Bool isinteger;
2262 	   SCIP_Bool branchcand;
2263 	   SCIP_Bool success;
2264 	   SCIP_Bool islocal;
2265 	   SCIP_Real refpointsunder[3] = {SCIP_INVALID, SCIP_INVALID, SCIP_INVALID};
2266 	   SCIP_Real refpointsover[3] = {SCIP_INVALID, SCIP_INVALID, SCIP_INVALID};
2267 	   SCIP_Bool overest[6] = {FALSE, FALSE, FALSE, TRUE, TRUE, TRUE};
2268 	   int i;
2269 	
2270 	   assert(scip != NULL);
2271 	   assert(expr != NULL);
2272 	
2273 	   child = SCIPexprGetChildren(expr)[0];
2274 	   assert(child != NULL);
2275 	
2276 	   childlb = bounds[0].inf;
2277 	   childub = bounds[0].sup;
2278 	
2279 	   /* if child is essentially constant, then there should be no point in separation */
2280 	   if( SCIPisEQ(scip, childlb, childub) )
2281 	   {
2282 	      SCIPdebugMsg(scip, "skip initestimates as child seems essentially fixed [%.15g,%.15g]\n", childlb, childub);
2283 	      return SCIP_OKAY;
2284 	   }
2285 	
2286 	   exprdata = SCIPexprGetData(expr);
2287 	   exponent = exprdata->exponent;
2288 	   assert(exponent != 1.0 && exponent != 0.0); /* this should have been simplified */
2289 	
2290 	   isinteger = EPSISINT(exponent, 0.0);
2291 	
2292 	   /* if exponent is not integral, then child must be non-negative */
2293 	   if( !isinteger && childlb < 0.0 )
2294 	   {
2295 	      /* somewhere we should have tightened the bound on x, but small tightening are not always applied by SCIP
2296 	       * it is ok to do this tightening here, but let's assert that we were close to 0.0 already
2297 	       */
2298 	      assert(SCIPisFeasZero(scip, childlb));
2299 	      childlb = 0.0;
2300 	   }
2301 	   assert(isinteger || childlb >= 0.0);
2302 	
2303 	   /* TODO simplify to get 3 refpoints for either under- or overestimate */
2304 	   SCIP_CALL( chooseRefpointsPow(scip, exprdata, childlb, childub, refpointsunder, refpointsover, !overestimate,
2305 	         overestimate) );
2306 	
2307 	   for( i = 0; i < 6 && *nreturned < SCIP_EXPR_MAXINITESTIMATES; ++i )
2308 	   {
2309 	      SCIP_Real refpoint;
2310 	
2311 	      if( (overest[i] && !overestimate) || (!overest[i] && overestimate) )
2312 	         continue;
2313 	
2314 	      assert(overest[i] || i < 3); /* make sure that no out-of-bounds array access will be attempted */
2315 	      refpoint = overest[i] ? refpointsover[i % 3] : refpointsunder[i % 3];
2316 	
2317 	      if( refpoint == SCIP_INVALID )
2318 	         continue;
2319 	
2320 	      assert(SCIPisLE(scip, refpoint, childub) && SCIPisGE(scip, refpoint, childlb));
2321 	
2322 	      branchcand = TRUE;
2323 	      SCIP_CALL( buildPowEstimator(scip, exprdata, overest[i], childlb, childub, childlb, childub,
2324 	            SCIPexprIsIntegral(child), refpoint, exponent, coefs[*nreturned], &constant[*nreturned],
2325 	            &success, &islocal, &branchcand) );
2326 	
2327 	      if( success )
2328 	      {
2329 	         SCIPdebugMsg(scip, "initestimate x^%g for base in [%g,%g] at ref=%g, over:%u -> %g*x+%g\n", exponent,
2330 	               childlb, childub, refpoint, overest[i], coefs[*nreturned][0], constant[*nreturned]);
2331 	         ++*nreturned;
2332 	      }
2333 	   }
2334 	
2335 	   return SCIP_OKAY;
2336 	}
2337 	
2338 	/** expression hash callback */
2339 	static
2340 	SCIP_DECL_EXPRHASH(hashPow)
2341 	{  /*lint --e{715}*/
2342 	   assert(scip != NULL);
2343 	   assert(expr != NULL);
2344 	   assert(SCIPexprGetNChildren(expr) == 1);
2345 	   assert(hashkey != NULL);
2346 	   assert(childrenhashes != NULL);
2347 	
2348 	   /* TODO include exponent into hashkey */
2349 	   *hashkey = POWEXPRHDLR_HASHKEY;
2350 	   *hashkey ^= childrenhashes[0];
2351 	
2352 	   return SCIP_OKAY;
2353 	}
2354 	
2355 	/** expression curvature detection callback */
2356 	static
2357 	SCIP_DECL_EXPRCURVATURE(curvaturePow)
2358 	{  /*lint --e{715}*/
2359 	   SCIP_EXPR* child;
2360 	   SCIP_INTERVAL childinterval;
2361 	   SCIP_Real exponent;
2362 	
2363 	   assert(scip != NULL);
2364 	   assert(expr != NULL);
2365 	   assert(exprcurvature != SCIP_EXPRCURV_UNKNOWN);
2366 	   assert(childcurv != NULL);
2367 	   assert(success != NULL);
2368 	   assert(SCIPexprGetNChildren(expr) == 1);
2369 	
2370 	   exponent = SCIPgetExponentExprPow(expr);
2371 	   child = SCIPexprGetChildren(expr)[0];
2372 	   assert(child != NULL);
2373 	
2374 	   SCIP_CALL( SCIPevalExprActivity(scip, child) );
2375 	   childinterval = SCIPexprGetActivity(child);
2376 	
2377 	   *childcurv = SCIPexprcurvPowerInv(childinterval, exponent, exprcurvature);
2378 	   /* SCIPexprcurvPowerInv return unknown actually means that curv cannot be obtained */
2379 	   *success = *childcurv != SCIP_EXPRCURV_UNKNOWN;
2380 	
2381 	   return SCIP_OKAY;
2382 	}
2383 	
2384 	/** expression monotonicity detection callback */
2385 	static
2386 	SCIP_DECL_EXPRMONOTONICITY(monotonicityPow)
2387 	{  /*lint --e{715}*/
2388 	   SCIP_INTERVAL interval;
2389 	   SCIP_Real exponent;
2390 	   SCIP_Real inf;
2391 	   SCIP_Real sup;
2392 	   SCIP_Bool expisint;
2393 	
2394 	   assert(scip != NULL);
2395 	   assert(expr != NULL);
2396 	   assert(result != NULL);
2397 	   assert(SCIPexprGetNChildren(expr) == 1);
2398 	   assert(childidx == 0);
2399 	
2400 	   assert(SCIPexprGetChildren(expr)[0] != NULL);
2401 	   SCIP_CALL( SCIPevalExprActivity(scip, SCIPexprGetChildren(expr)[0]) );
2402 	   interval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]);
2403 	
2404 	   *result = SCIP_MONOTONE_UNKNOWN;
2405 	   inf = SCIPintervalGetInf(interval);
2406 	   sup = SCIPintervalGetSup(interval);
2407 	   exponent = SCIPgetExponentExprPow(expr);
2408 	   expisint = EPSISINT(exponent, 0.0); /*lint !e835*/
2409 	
2410 	   if( expisint )
2411 	   {
2412 	      SCIP_Bool expisodd = ceil(exponent/2) != exponent/2;
2413 	
2414 	      if( expisodd )
2415 	      {
2416 	         /* x^1, x^3, ... */
2417 	         if( exponent >= 0.0 )
2418 	            *result = SCIP_MONOTONE_INC;
2419 	
2420 	         /* ..., x^-3, x^-1 are decreasing if 0 is not in ]inf,sup[ */
2421 	         else if( inf >= 0.0 || sup <= 0.0 )
2422 	            *result = SCIP_MONOTONE_DEC;
2423 	      }
2424 	      /* ..., x^-4, x^-2, x^2, x^4, ... */
2425 	      else
2426 	      {
2427 	         /* function is not monotone if 0 is in ]inf,sup[ */
2428 	         if( inf >= 0.0 )
2429 	            *result = exponent >= 0.0 ? SCIP_MONOTONE_INC : SCIP_MONOTONE_DEC;
2430 	         else if( sup <= 0.0 )
2431 	            *result = exponent >= 0.0 ? SCIP_MONOTONE_DEC : SCIP_MONOTONE_INC;
2432 	      }
2433 	   }
2434 	   else
2435 	   {
2436 	      /* note that the expression is not defined for negative input values
2437 	       *
2438 	       * - increasing iff exponent >= 0
2439 	       * - decreasing iff exponent <= 0
2440 	       */
2441 	      *result = exponent >= 0.0 ? SCIP_MONOTONE_INC : SCIP_MONOTONE_DEC;
2442 	   }
2443 	
2444 	   return SCIP_OKAY;
2445 	}
2446 	
2447 	/** expression integrality detection callback */
2448 	static
2449 	SCIP_DECL_EXPRINTEGRALITY(integralityPow)
2450 	{  /*lint --e{715}*/
2451 	   SCIP_EXPR* child;
2452 	   SCIP_Real exponent;
2453 	   SCIP_Bool expisint;
2454 	
2455 	   assert(scip != NULL);
2456 	   assert(expr != NULL);
2457 	   assert(isintegral != NULL);
2458 	   assert(SCIPexprGetNChildren(expr) == 1);
2459 	
2460 	   *isintegral = FALSE;
2461 	
2462 	   child = SCIPexprGetChildren(expr)[0];
2463 	   assert(child != NULL);
2464 	
2465 	   /* expression can not be integral if child is not */
2466 	   if( !SCIPexprIsIntegral(child) )
2467 	      return SCIP_OKAY;
2468 	
2469 	   exponent = SCIPgetExponentExprPow(expr);
2470 	   assert(exponent != 0.0);
2471 	   expisint = EPSISINT(exponent, 0.0); /*lint !e835*/
2472 	
2473 	   /* expression is integral if and only if exponent non-negative and integral */
2474 	   *isintegral = expisint && exponent >= 0.0;
2475 	
2476 	   return SCIP_OKAY;
2477 	}
2478 	
2479 	/** simplifies a signpower expression
2480 	 */
2481 	static
2482 	SCIP_DECL_EXPRSIMPLIFY(simplifySignpower)
2483 	{  /*lint --e{715}*/
2484 	   SCIP_EXPR* base;
2485 	   SCIP_Real exponent;
2486 	
2487 	   assert(scip != NULL);
2488 	   assert(expr != NULL);
2489 	   assert(simplifiedexpr != NULL);
2490 	   assert(SCIPexprGetNChildren(expr) == 1);
2491 	
2492 	   base = SCIPexprGetChildren(expr)[0];
2493 	   assert(base != NULL);
2494 	
2495 	   exponent = SCIPgetExponentExprPow(expr);
2496 	   SCIPdebugPrintf("[simplifySignpower] simplifying power with expo %g\n", exponent);
2497 	   assert(exponent >= 1.0);
2498 	
2499 	   /* enforces SPOW2 */
2500 	   if( exponent == 1.0 )
2501 	   {
2502 	      SCIPdebugPrintf("[simplifySignpower] POW2\n");
2503 	      *simplifiedexpr = base;
2504 	      SCIPcaptureExpr(*simplifiedexpr);
2505 	      return SCIP_OKAY;
2506 	   }
2507 	
2508 	   /* enforces SPOW3 */
2509 	   if( SCIPisExprValue(scip, base) )
2510 	   {
2511 	      SCIP_Real baseval;
2512 	
2513 	      SCIPdebugPrintf("[simplifySignpower] POW3\n");
2514 	      baseval = SCIPgetValueExprValue(base);
2515 	
2516 	      SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, SIGN(baseval) * pow(REALABS(baseval), exponent),
2517 	                 ownercreate, ownercreatedata) );
2518 	
2519 	      return SCIP_OKAY;
2520 	   }
2521 	
2522 	   /* enforces SPOW11 (exp(x)^n = exp(n*x))
2523 	    * since exp() is always nonnegative, we can treat signpower as normal power here
2524 	    */
2525 	   if( SCIPisExprExp(scip, base) )
2526 	   {
2527 	      SCIP_EXPR* child;
2528 	      SCIP_EXPR* prod;
2529 	      SCIP_EXPR* exponential;
2530 	      SCIP_EXPR* simplifiedprod;
2531 	
2532 	      SCIPdebugPrintf("[simplifySignpower] POW11\n");
2533 	      child = SCIPexprGetChildren(base)[0];
2534 	
2535 	      /* multiply child of exponential with exponent */
2536 	      SCIP_CALL( SCIPcreateExprProduct(scip, &prod, 1, &child, exponent, ownercreate, ownercreatedata) );
2537 	
2538 	      /* simplify product */
2539 	      SCIP_CALL( SCIPcallExprSimplify(scip, prod, &simplifiedprod, ownercreate, ownercreatedata) );
2540 	      SCIP_CALL( SCIPreleaseExpr(scip, &prod) );
2541 	
2542 	      /* create exponential with new child */
2543 	      SCIP_CALL( SCIPcreateExprExp(scip, &exponential, simplifiedprod, ownercreate, ownercreatedata) );
2544 	      SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedprod) );
2545 	
2546 	      /* the final simplified expression is the simplification of the just created exponential */
2547 	      SCIP_CALL( SCIPcallExprSimplify(scip, exponential, simplifiedexpr, ownercreate, ownercreatedata) );
2548 	      SCIP_CALL( SCIPreleaseExpr(scip, &exponential) );
2549 	
2550 	      return SCIP_OKAY;
2551 	   }
2552 	
2553 	   /* enforces SPOW6 */
2554 	   if( EPSISINT(exponent, 0.0) && ((int)exponent) % 2 == 1 )
2555 	   {
2556 	      SCIP_EXPR* aux;
2557 	
2558 	      /* we do not just change the expression data of expression to say it is a normal power, since, at the moment,
2559 	       * simplify identifies that expressions changed by checking that the pointer of the input expression is
2560 	       * different from the returned (simplified) expression
2561 	      */
2562 	      SCIP_CALL( SCIPcreateExprPow(scip, &aux, base, exponent, ownercreate, ownercreatedata) );
2563 	
2564 	      SCIP_CALL( simplifyPow(scip, aux, simplifiedexpr, ownercreate, ownercreatedata) );
2565 	      SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
2566 	
2567 	      return SCIP_OKAY;
2568 	   }
2569 	
2570 	   /* enforces SPOW10 */
2571 	   if( SCIPisExprVar(scip, base) )
2572 	   {
2573 	      SCIP_VAR* basevar;
2574 	
2575 	      SCIPdebugPrintf("[simplifySignpower] POW10\n");
2576 	      basevar = SCIPgetVarExprVar(base);
2577 	
2578 	      assert(basevar != NULL);
2579 	
2580 	      if( SCIPvarIsBinary(basevar) )
2581 	      {
2582 	         *simplifiedexpr = base;
2583 	         SCIPcaptureExpr(*simplifiedexpr);
2584 	         return SCIP_OKAY;
2585 	      }
2586 	   }
2587 	
2588 	   /* TODO if( SCIPisExprSignpower(scip, base) ... */
2589 	
2590 	   /* enforces SPOW8
2591 	    * given (signpow n (pow expo expr)) we distribute the exponent:
2592 	    * -> (signpow n*expo expr) for even n  (i.e., sign(x^n) * |x|^n = 1 * x^n)
2593 	    * notes: n is an even integer (see SPOW6 above)
2594 	    * FIXME: doesn't this extend to any exponent?
2595 	    * If (pow expo expr) can be negative, it should mean that (-1)^expo = -1
2596 	    * then (signpow n (pow expo expr)) = sign(expr^expo) * |expr^expo|^n
2597 	    * then sign(expr^expo) = sign(expr) and |expr^expo| = |expr|^expo and so
2598 	    * (signpow n (pow expo expr)) = sign(expr^expo) * |expr^expo|^n = sign(expr) * |expr|^(expo*n) = signpow n*expo expr
2599 	    */
2600 	   if( EPSISINT(exponent, 0.0) && SCIPisExprPower(scip, base) )
2601 	   {
2602 	      SCIP_EXPR* aux;
2603 	      SCIP_Real newexponent;
2604 	
2605 	      assert(((int)exponent) % 2 == 0 ); /* odd case should have been handled by SPOW6 */
2606 	
2607 	      newexponent = SCIPgetExponentExprPow(base) * exponent;
2608 	      SCIP_CALL( SCIPcreateExprSignpower(scip, &aux, SCIPexprGetChildren(base)[0], newexponent,
2609 	                 ownercreate, ownercreatedata) );
2610 	      SCIP_CALL( simplifySignpower(scip, aux, simplifiedexpr, ownercreate, ownercreatedata) );
2611 	
2612 	      SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
2613 	
2614 	      return SCIP_OKAY;
2615 	   }
2616 	
2617 	   /* enforces SPOW9 */
2618 	   if( SCIPisExprSum(scip, base)
2619 	      && SCIPexprGetNChildren(base) == 1
2620 	      && SCIPgetConstantExprSum(base) == 0.0 )
2621 	   {
2622 	      SCIP_EXPR* simplifiedaux;
2623 	      SCIP_EXPR* aux;
2624 	      SCIP_Real newcoef;
2625 	
2626 	      SCIPdebugPrintf("[simplifySignpower] seeing a sum with one term, exponent %g\n", exponent);
2627 	      /* assert SS7 holds */
2628 	      assert(SCIPgetCoefsExprSum(base)[0] != 1.0);
2629 	
2630 	      /* create (signpow n expr) and simplify it
2631 	       * note: we call simplifySignpower directly, since we know that `expr` is simplified */
2632 	      SCIP_CALL( SCIPcreateExprSignpower(scip, &aux, SCIPexprGetChildren(base)[0], exponent,
2633 	                 ownercreate, ownercreatedata) );
2634 	      newcoef = SIGN(SCIPgetCoefsExprSum(base)[0]) * pow(REALABS(SCIPgetCoefsExprSum(base)[0]), exponent);
2635 	      SCIP_CALL( simplifySignpower(scip, aux, &simplifiedaux, ownercreate, ownercreatedata) );
2636 	      SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
2637 	
2638 	      /* create (sum (signpow n expr)) and simplify it
2639 	       * this calls simplifySum directly, since we know its child is simplified */
2640 	      SCIP_CALL( SCIPcreateExprSum(scip, &aux, 1, &simplifiedaux, &newcoef, 0.0, ownercreate, ownercreatedata) );
2641 	      SCIP_CALL( SCIPcallExprSimplify(scip, aux, simplifiedexpr, ownercreate, ownercreatedata) );
2642 	      SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
2643 	      SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedaux) );
2644 	      return SCIP_OKAY;
2645 	   }
2646 	
2647 	   SCIPdebugPrintf("[simplifySignpower] signpower is simplified\n");
2648 	   *simplifiedexpr = expr;
2649 	
2650 	   /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */
2651 	   SCIPcaptureExpr(*simplifiedexpr);
2652 	
2653 	   return SCIP_OKAY;
2654 	}
2655 	
2656 	/** expression handler copy callback */
2657 	static
2658 	SCIP_DECL_EXPRCOPYHDLR(copyhdlrSignpower)
2659 	{  /*lint --e{715}*/
2660 	   SCIP_CALL( SCIPincludeExprhdlrSignpower(scip) );
2661 	
2662 	   return SCIP_OKAY;
2663 	}
2664 	
2665 	/** expression print callback */
2666 	static
2667 	SCIP_DECL_EXPRPRINT(printSignpower)
2668 	{  /*lint --e{715}*/
2669 	   assert(expr != NULL);
2670 	
2671 	   switch( stage )
2672 	   {
2673 	      case SCIP_EXPRITER_ENTEREXPR :
2674 	      {
2675 	         SCIPinfoMessage(scip, file, "signpower(");
2676 	         break;
2677 	      }
2678 	
2679 	      case SCIP_EXPRITER_VISITINGCHILD :
2680 	      {
2681 	         assert(currentchild == 0);
2682 	         break;
2683 	      }
2684 	
2685 	      case SCIP_EXPRITER_LEAVEEXPR :
2686 	      {
2687 	         SCIPinfoMessage(scip, file, ",%g)", SCIPgetExponentExprPow(expr));
2688 	         break;
2689 	      }
2690 	
2691 	      case SCIP_EXPRITER_VISITEDCHILD :
2692 	      default:
2693 	         break;
2694 	   }
2695 	
2696 	   return SCIP_OKAY;
2697 	}
2698 	
2699 	/** expression parse callback */
2700 	static
2701 	SCIP_DECL_EXPRPARSE(parseSignpower)
2702 	{  /*lint --e{715}*/
2703 	   SCIP_EXPR* childexpr;
2704 	   SCIP_Real exponent;
2705 	
2706 	   assert(expr != NULL);
2707 	
2708 	   /**! [SnippetExprParseSignpower] */
2709 	   /* parse child expression string */
2710 	   SCIP_CALL( SCIPparseExpr(scip, &childexpr, string, endstring, ownercreate, ownercreatedata) );
2711 	   assert(childexpr != NULL);
2712 	
2713 	   string = *endstring;
2714 	   while( *string == ' ' )
2715 	      ++string;
2716 	
2717 	   if( *string != ',' )
2718 	   {
2719 	      SCIPerrorMessage("Expected comma after first argument of signpower().\n");
2720 	      return SCIP_READERROR;
2721 	   }
2722 	   ++string;
2723 	
2724 	   if( !SCIPparseReal(scip, string, &exponent, (char**)endstring) )
2725 	   {
2726 	      SCIPerrorMessage("Expected numeric exponent for second argument of signpower().\n");
2727 	      return SCIP_READERROR;
2728 	   }
2729 	
2730 	   if( exponent <= 1.0 || SCIPisInfinity(scip, exponent) )
2731 	   {
2732 	      SCIPerrorMessage("Expected finite exponent >= 1.0 for signpower().\n");
2733 	      return SCIP_READERROR;
2734 	   }
2735 	
2736 	   /* create signpower expression */
2737 	   SCIP_CALL( SCIPcreateExprSignpower(scip, expr, childexpr, exponent, ownercreate, ownercreatedata) );
2738 	   assert(*expr != NULL);
2739 	
2740 	   /* release child expression since it has been captured by the signpower expression */
2741 	   SCIP_CALL( SCIPreleaseExpr(scip, &childexpr) );
2742 	
2743 	   *success = TRUE;
2744 	   /**! [SnippetExprParseSignpower] */
2745 	
2746 	   return SCIP_OKAY;
2747 	}
2748 	
2749 	/** expression point evaluation callback */
2750 	static
2751 	SCIP_DECL_EXPREVAL(evalSignpower)
2752 	{  /*lint --e{715}*/
2753 	   SCIP_Real exponent;
2754 	   SCIP_Real base;
2755 	
2756 	   assert(expr != NULL);
2757 	   assert(SCIPexprGetNChildren(expr) == 1);
2758 	   assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]) != SCIP_INVALID);
2759 	
2760 	   exponent = SCIPgetExponentExprPow(expr);
2761 	   base = SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]);
2762 	
2763 	   *val = SIGN(base) * pow(REALABS(base), exponent);
2764 	
2765 	   /* if there is a range error, pow() should return some kind of infinity, or HUGE_VAL
2766 	    * we could also work with floating point exceptions or errno, but I am not sure this would be thread-safe
2767 	    */
2768 	   if( !SCIPisFinite(*val) || *val == HUGE_VAL || *val == -HUGE_VAL )
2769 	      *val = SCIP_INVALID;
2770 	
2771 	   return SCIP_OKAY;
2772 	}
2773 	
2774 	/** expression derivative evaluation callback */
2775 	static
2776 	SCIP_DECL_EXPRBWDIFF(bwdiffSignpower)
2777 	{  /*lint --e{715}*/
2778 	   SCIP_EXPR* child;
2779 	   SCIP_Real childval;
2780 	   SCIP_Real exponent;
2781 	
2782 	   assert(expr != NULL);
2783 	   assert(SCIPexprGetData(expr) != NULL);
2784 	   assert(childidx == 0);
2785 	
2786 	   child = SCIPexprGetChildren(expr)[0];
2787 	   assert(child != NULL);
2788 	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(child)), "val") != 0);
2789 	
2790 	   childval = SCIPexprGetEvalValue(child);
2791 	   assert(childval != SCIP_INVALID);
2792 	
2793 	   exponent = SCIPgetExponentExprPow(expr);
2794 	   assert(exponent >= 1.0);
2795 	
2796 	   *val = exponent * pow(REALABS(childval), exponent - 1.0);
2797 	
2798 	   return SCIP_OKAY;
2799 	}
2800 	
2801 	/** expression interval evaluation callback */
2802 	static
2803 	SCIP_DECL_EXPRINTEVAL(intevalSignpower)
2804 	{  /*lint --e{715}*/
2805 	   SCIP_INTERVAL childinterval;
2806 	
2807 	   assert(expr != NULL);
2808 	   assert(SCIPexprGetNChildren(expr) == 1);
2809 	
2810 	   childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]);
2811 	   if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) )
2812 	   {
2813 	      SCIPintervalSetEmpty(interval);
2814 	      return SCIP_OKAY;
2815 	   }
2816 	
2817 	   SCIPintervalSignPowerScalar(SCIP_INTERVAL_INFINITY, interval, childinterval, SCIPgetExponentExprPow(expr));
2818 	
2819 	   return SCIP_OKAY;
2820 	}
2821 	
2822 	/** expression estimator callback */
2823 	static
2824 	SCIP_DECL_EXPRESTIMATE(estimateSignpower)
2825 	{  /*lint --e{715}*/
2826 	   SCIP_EXPRDATA* exprdata;
2827 	   SCIP_Real childlb;
2828 	   SCIP_Real childub;
2829 	   SCIP_Real childglb;
2830 	   SCIP_Real childgub;
2831 	   SCIP_Real exponent;
2832 	
2833 	   assert(scip != NULL);
2834 	   assert(expr != NULL);
2835 	   assert(SCIPexprGetNChildren(expr) == 1);
2836 	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), SIGNPOWEXPRHDLR_NAME) == 0);
2837 	   assert(refpoint != NULL);
2838 	   assert(coefs != NULL);
2839 	   assert(constant != NULL);
2840 	   assert(islocal != NULL);
2841 	   assert(branchcand != NULL);
2842 	   assert(*branchcand == TRUE);  /* the default */
2843 	   assert(success != NULL);
2844 	
2845 	   *success = FALSE;
2846 	
2847 	   SCIPdebugMsg(scip, "%sestimation of signed x^%g at x=%g\n", overestimate ? "over" : "under",
2848 	         SCIPgetExponentExprPow(expr), *refpoint);
2849 	
2850 	   /* we can not generate a cut at +/- infinity */
2851 	   if( SCIPisInfinity(scip, REALABS(*refpoint)) )
2852 	      return SCIP_OKAY;
2853 	
2854 	   childlb = localbounds[0].inf;
2855 	   childub = localbounds[0].sup;
2856 	
2857 	   childglb = globalbounds[0].inf;
2858 	   childgub = globalbounds[0].sup;
2859 	
2860 	   exprdata = SCIPexprGetData(expr);
2861 	   exponent = exprdata->exponent;
2862 	   assert(exponent > 1.0); /* exponent == 1 should have been simplified */
2863 	
2864 	   /* if child is constant, then return a constant estimator
2865 	    * this can help with small infeasibilities if boundtightening is relaxing bounds too much
2866 	    */
2867 	   if( childlb == childub )
2868 	   {
2869 	      *coefs = 0.0;
2870 	      *constant = SIGN(childlb)*pow(REALABS(childlb), exponent);
2871 	      *success = TRUE;
2872 	      *islocal = childglb != childgub;
2873 	      *branchcand = FALSE;
2874 	      return SCIP_OKAY;
2875 	   }
2876 	
2877 	   if( childlb >= 0.0 )
2878 	   {
2879 	      estimateParabola(scip, exponent, overestimate, childlb, childub, MAX(0.0, *refpoint), constant, coefs,
2880 	            islocal, success);
2881 	
2882 	      *branchcand = *islocal;
2883 	
2884 	      /* if odd or signed power, then check whether tangent on parabola is also globally valid, that is
2885 	       * reference point is right of -root*global-lower-bound
2886 	       */
2887 	      if( !*islocal && childglb < 0.0 )
2888 	      {
2889 	         if( SCIPisInfinity(scip, -childglb) )
2890 	            *islocal = TRUE;
2891 	         else
2892 	         {
2893 	            if( exprdata->root == SCIP_INVALID )
2894 	            {
2895 	               SCIP_CALL( computeSignpowerRoot(scip, &exprdata->root, exponent) );
2896 	            }
2897 	            *islocal = *refpoint < exprdata->root * (-childglb);
2898 	         }
2899 	      }
2900 	   }
2901 	   else  /* and childlb < 0.0 due to previous if */
2902 	   {
2903 	      /* compute root if not known yet; only needed if mixed sign (global child ub > 0) */
2904 	      if( exprdata->root == SCIP_INVALID && childgub > 0.0 )
2905 	      {
2906 	         SCIP_CALL( computeSignpowerRoot(scip, &exprdata->root, exponent) );
2907 	      }
2908 	      estimateSignedpower(scip, exponent, exprdata->root, overestimate, childlb, childub, *refpoint,
2909 	         childglb, childgub, constant, coefs, islocal, branchcand, success);
2910 	   }
2911 	
2912 	   return SCIP_OKAY;
2913 	}
2914 	
2915 	/** initial estimates callback for a signpower expression */
2916 	static
2917 	SCIP_DECL_EXPRINITESTIMATES(initestimatesSignpower)
2918 	{
2919 	   SCIP_EXPRDATA* exprdata;
2920 	   SCIP_Real childlb;
2921 	   SCIP_Real childub;
2922 	   SCIP_Real exponent;
2923 	   SCIP_Bool branchcand;
2924 	   SCIP_Bool success;
2925 	   SCIP_Bool islocal;
2926 	   SCIP_Real refpointsunder[3] = {SCIP_INVALID, SCIP_INVALID, SCIP_INVALID};
2927 	   SCIP_Real refpointsover[3] = {SCIP_INVALID, SCIP_INVALID, SCIP_INVALID};
2928 	   SCIP_Bool overest[6] = {FALSE, FALSE, FALSE, TRUE, TRUE, TRUE};
2929 	   SCIP_Real refpoint;
2930 	   int i;
2931 	
2932 	   assert(scip != NULL);
2933 	   assert(expr != NULL);
2934 	   assert(SCIPexprGetNChildren(expr) == 1);
2935 	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), SIGNPOWEXPRHDLR_NAME) == 0);
2936 	
2937 	   childlb = bounds[0].inf;
2938 	   childub = bounds[0].sup;
2939 	
2940 	   /* if child is essentially constant, then there should be no point in separation */
2941 	   if( SCIPisEQ(scip, childlb, childub) )
2942 	   {
2943 	      SCIPdebugMsg(scip, "skip initestimates as child seems essentially fixed [%.15g,%.15g]\n", childlb, childub);
2944 	      return SCIP_OKAY;
2945 	   }
2946 	
2947 	   exprdata = SCIPexprGetData(expr);
2948 	   exponent = exprdata->exponent;
2949 	   assert(exponent > 1.0); /* this should have been simplified */
2950 	
2951 	   if( childlb >= 0.0 )
2952 	   {
2953 	      if( !overestimate )
2954 	         addTangentRefpoints(scip, exponent, childlb, childub, refpointsunder);
2955 	      if( overestimate && !SCIPisInfinity(scip, childub) )
2956 	         refpointsover[0] = (childlb + childub) / 2.0;
2957 	   }
2958 	   else if( childub <= 0.0 )
2959 	   {
2960 	      if( !overestimate && !SCIPisInfinity(scip, -childlb) )
2961 	         refpointsunder[0] = (childlb + childub) / 2.0;
2962 	      if( overestimate )
2963 	         addTangentRefpoints(scip, exponent, childlb, childub, refpointsunder);
2964 	   }
2965 	   else
2966 	   {
2967 	      SCIP_CALL( addSignpowerRefpoints(scip, exprdata, childlb, childub, exponent, !overestimate, refpointsunder) );
2968 	   }
2969 	
2970 	   /* add cuts for all refpoints */
2971 	   for( i = 0; i < 6 && *nreturned < SCIP_EXPR_MAXINITESTIMATES; ++i )
2972 	   {
2973 	      if( (overest[i] && !overestimate) || (!overest[i] && overestimate) )
2974 	         continue;
2975 	
2976 	      assert(overest[i] || i < 3); /* make sure that no out-of-bounds array access will be attempted */
2977 	      refpoint = overest[i] ? refpointsover[i % 3] : refpointsunder[i % 3];
2978 	      if( refpoint == SCIP_INVALID )
2979 	         continue;
2980 	      assert(SCIPisLE(scip, refpoint, childub) && SCIPisGE(scip, refpoint, childlb));
2981 	
2982 	      if( childlb >= 0 )
2983 	      {
2984 	         estimateParabola(scip, exponent, overest[i], childlb, childub, refpoint, &constant[*nreturned], coefs[*nreturned],
2985 	               &islocal, &success);
2986 	      }
2987 	      else
2988 	      {
2989 	         /* compute root if not known yet; only needed if mixed sign (global child ub > 0) */
2990 	         if( exprdata->root == SCIP_INVALID && childub > 0.0 )
2991 	         {
2992 	            SCIP_CALL( computeSignpowerRoot(scip, &exprdata->root, exponent) );
2993 	         }
2994 	         branchcand = TRUE;
2995 	         estimateSignedpower(scip, exponent, exprdata->root, overest[i], childlb, childub, refpoint,
2996 	               childlb, childub, &constant[*nreturned], coefs[*nreturned], &islocal,
2997 	               &branchcand, &success);
2998 	      }
2999 	
3000 	      if( success )
3001 	         ++*nreturned;
3002 	   }
3003 	
3004 	   return SCIP_OKAY;
3005 	}
3006 	
3007 	/** expression reverse propagaton callback */
3008 	static
3009 	SCIP_DECL_EXPRREVERSEPROP(reversepropSignpower)
3010 	{  /*lint --e{715}*/
3011 	   SCIP_INTERVAL interval;
3012 	   SCIP_INTERVAL exprecip;
3013 	   SCIP_Real exponent;
3014 	
3015 	   assert(scip != NULL);
3016 	   assert(expr != NULL);
3017 	   assert(SCIPexprGetNChildren(expr) == 1);
3018 	
3019 	   exponent = SCIPgetExponentExprPow(expr);
3020 	
3021 	   SCIPdebugMsg(scip, "reverseprop signpow(x,%g) in [%.15g,%.15g]", exponent, bounds.inf, bounds.sup);
3022 	
3023 	   if( SCIPintervalIsEntire(SCIP_INTERVAL_INFINITY, bounds) )
3024 	   {
3025 	      SCIPdebugMsgPrint(scip, "-> no improvement\n");
3026 	      return SCIP_OKAY;
3027 	   }
3028 	
3029 	   /* f = pow(c0, alpha) -> c0 = pow(f, 1/alpha) */
3030 	   SCIPintervalSet(&exprecip, exponent);
3031 	   SCIPintervalReciprocal(SCIP_INTERVAL_INFINITY, &exprecip, exprecip);
3032 	   if( exprecip.inf == exprecip.sup )
3033 	   {
3034 	      SCIPintervalSignPowerScalar(SCIP_INTERVAL_INFINITY, &interval, bounds, exprecip.inf);
3035 	   }
3036 	   else
3037 	   {
3038 	      SCIP_INTERVAL interval1, interval2;
3039 	      SCIPintervalSignPowerScalar(SCIP_INTERVAL_INFINITY, &interval1, bounds, exprecip.inf);
3040 	      SCIPintervalSignPowerScalar(SCIP_INTERVAL_INFINITY, &interval2, bounds, exprecip.sup);
3041 	      SCIPintervalUnify(&interval, interval1, interval2);
3042 	   }
3043 	
3044 	   SCIPdebugMsgPrint(scip, " -> [%.15g,%.15g]\n", interval.inf, interval.sup);
3045 	
3046 	   childrenbounds[0] = interval;
3047 	
3048 	   return SCIP_OKAY;
3049 	}
3050 	
3051 	/** expression hash callback */
3052 	static
3053 	SCIP_DECL_EXPRHASH(hashSignpower)
3054 	{  /*lint --e{715}*/
3055 	   assert(scip != NULL);
3056 	   assert(expr != NULL);
3057 	   assert(SCIPexprGetNChildren(expr) == 1);
3058 	   assert(hashkey != NULL);
3059 	   assert(childrenhashes != NULL);
3060 	
3061 	   /* TODO include exponent into hashkey */
3062 	   *hashkey = SIGNPOWEXPRHDLR_HASHKEY;
3063 	   *hashkey ^= childrenhashes[0];
3064 	
3065 	   return SCIP_OKAY;
3066 	}
3067 	
3068 	/** expression curvature detection callback */
3069 	static
3070 	SCIP_DECL_EXPRCURVATURE(curvatureSignpower)
3071 	{  /*lint --e{715}*/
3072 	   SCIP_EXPR* child;
3073 	   SCIP_INTERVAL childinterval;
3074 	
3075 	   assert(scip != NULL);
3076 	   assert(expr != NULL);
3077 	   assert(exprcurvature != SCIP_EXPRCURV_UNKNOWN);
3078 	   assert(childcurv != NULL);
3079 	   assert(success != NULL);
3080 	   assert(SCIPexprGetNChildren(expr) == 1);
3081 	
3082 	   child = SCIPexprGetChildren(expr)[0];
3083 	   assert(child != NULL);
3084 	
3085 	   SCIP_CALL( SCIPevalExprActivity(scip, child) );
3086 	   childinterval = SCIPexprGetActivity(child);
3087 	
3088 	   if( exprcurvature == SCIP_EXPRCURV_CONVEX )
3089 	   {
3090 	      /* signpower is only convex if argument is convex and non-negative */
3091 	      *childcurv = SCIP_EXPRCURV_CONVEX;
3092 	      *success = childinterval.inf >= 0.0;
3093 	   }
3094 	   else if( exprcurvature == SCIP_EXPRCURV_CONCAVE )
3095 	   {
3096 	      /* signpower is only concave if argument is concave and non-positive */
3097 	      *childcurv = SCIP_EXPRCURV_CONCAVE;
3098 	      *success = childinterval.sup <= 0.0;
3099 	   }
3100 	   else
3101 	      *success = FALSE;
3102 	
3103 	   return SCIP_OKAY;
3104 	}
3105 	
3106 	/** expression monotonicity detection callback */
3107 	static
3108 	SCIP_DECL_EXPRMONOTONICITY(monotonicitySignpower)
3109 	{  /*lint --e{715}*/
3110 	   assert(scip != NULL);
3111 	   assert(expr != NULL);
3112 	   assert(result != NULL);
3113 	
3114 	   *result = SCIP_MONOTONE_INC;
3115 	   return SCIP_OKAY;
3116 	}
3117 	
3118 	/** creates the handler for power expression and includes it into SCIP */
3119 	SCIP_RETCODE SCIPincludeExprhdlrPow(
3120 	   SCIP*                 scip                /**< SCIP data structure */
3121 	   )
3122 	{
3123 	   SCIP_EXPRHDLR* exprhdlr;
3124 	   SCIP_EXPRHDLRDATA* exprhdlrdata;
3125 	
3126 	   SCIP_CALL( SCIPallocClearBlockMemory(scip, &exprhdlrdata) );
3127 	
3128 	   SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, POWEXPRHDLR_NAME, POWEXPRHDLR_DESC, POWEXPRHDLR_PRECEDENCE,
3129 	              evalPow, exprhdlrdata) );
3130 	   assert(exprhdlr != NULL);
3131 	
3132 	   SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrPow, freehdlrPow);
3133 	   SCIPexprhdlrSetCopyFreeData(exprhdlr, copydataPow, freedataPow);
3134 	   SCIPexprhdlrSetSimplify(exprhdlr, simplifyPow);
3135 	   SCIPexprhdlrSetPrint(exprhdlr, printPow);
3136 	   SCIPexprhdlrSetIntEval(exprhdlr, intevalPow);
3137 	   SCIPexprhdlrSetEstimate(exprhdlr, initestimatesPow, estimatePow);
3138 	   SCIPexprhdlrSetReverseProp(exprhdlr, reversepropPow);
3139 	   SCIPexprhdlrSetHash(exprhdlr, hashPow);
3140 	   SCIPexprhdlrSetCompare(exprhdlr, comparePow);
3141 	   SCIPexprhdlrSetDiff(exprhdlr, bwdiffPow, fwdiffPow, bwfwdiffPow);
3142 	   SCIPexprhdlrSetCurvature(exprhdlr, curvaturePow);
3143 	   SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicityPow);
3144 	   SCIPexprhdlrSetIntegrality(exprhdlr, integralityPow);
3145 	   SCIPexprhdlrSetGetSymdata(exprhdlr, getSymDataPow);
3146 	
3147 	   SCIP_CALL( SCIPaddRealParam(scip, "expr/" POWEXPRHDLR_NAME "/minzerodistance",
3148 	      "minimal distance from zero to enforce for child in bound tightening",
3149 	      &exprhdlrdata->minzerodistance, FALSE, SCIPepsilon(scip), 0.0, 1.0, NULL, NULL) );
3150 	
3151 	   SCIP_CALL( SCIPaddIntParam(scip, "expr/" POWEXPRHDLR_NAME "/expandmaxexponent",
3152 	      "maximal exponent when to expand power of sum in simplify",
3153 	      &exprhdlrdata->expandmaxexponent, FALSE, 2, 1, INT_MAX, NULL, NULL) );
3154 	
3155 	   SCIP_CALL( SCIPaddBoolParam(scip, "expr/" POWEXPRHDLR_NAME "/distribfracexponent",
3156 	      "whether a fractional exponent is distributed onto factors on power of product",
3157 	      &exprhdlrdata->distribfracexponent, FALSE, FALSE, NULL, NULL) );
3158 	
3159 	   return SCIP_OKAY;
3160 	}
3161 	
3162 	/** creates the handler for signed power expression and includes it into SCIP */
3163 	SCIP_RETCODE SCIPincludeExprhdlrSignpower(
3164 	   SCIP*                 scip                /**< SCIP data structure */
3165 	   )
3166 	{
3167 	   SCIP_EXPRHDLR* exprhdlr;
3168 	
3169 	   SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, SIGNPOWEXPRHDLR_NAME, SIGNPOWEXPRHDLR_DESC,
3170 	              SIGNPOWEXPRHDLR_PRECEDENCE, evalSignpower, NULL) );
3171 	   assert(exprhdlr != NULL);
3172 	
3173 	   SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrSignpower, NULL);
3174 	   SCIPexprhdlrSetCopyFreeData(exprhdlr, copydataPow, freedataPow);
3175 	   SCIPexprhdlrSetSimplify(exprhdlr, simplifySignpower);
3176 	   SCIPexprhdlrSetPrint(exprhdlr, printSignpower);
3177 	   SCIPexprhdlrSetParse(exprhdlr, parseSignpower);
3178 	   SCIPexprhdlrSetIntEval(exprhdlr, intevalSignpower);
3179 	   SCIPexprhdlrSetEstimate(exprhdlr, initestimatesSignpower, estimateSignpower);
3180 	   SCIPexprhdlrSetReverseProp(exprhdlr, reversepropSignpower);
3181 	   SCIPexprhdlrSetHash(exprhdlr, hashSignpower);
3182 	   SCIPexprhdlrSetCompare(exprhdlr, comparePow);
3183 	   SCIPexprhdlrSetDiff(exprhdlr, bwdiffSignpower, NULL, NULL);
3184 	   SCIPexprhdlrSetCurvature(exprhdlr, curvatureSignpower);
3185 	   SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicitySignpower);
3186 	   SCIPexprhdlrSetIntegrality(exprhdlr, integralityPow);
3187 	   SCIPexprhdlrSetGetSymdata(exprhdlr, getSymDataPow);
3188 	
3189 	   return SCIP_OKAY;
3190 	}
3191 	
3192 	/** creates a power expression */
3193 	SCIP_RETCODE SCIPcreateExprPow(
3194 	   SCIP*                 scip,               /**< SCIP data structure */
3195 	   SCIP_EXPR**           expr,               /**< pointer where to store expression */
3196 	   SCIP_EXPR*            child,              /**< single child */
3197 	   SCIP_Real             exponent,           /**< exponent of the power expression */
3198 	   SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
3199 	   void*                 ownercreatedata     /**< data to pass to ownercreate */
3200 	   )
3201 	{
3202 	   SCIP_EXPRDATA* exprdata;
3203 	
3204 	   assert(expr != NULL);
3205 	   assert(child != NULL);
3206 	
3207 	   SCIP_CALL( createData(scip, &exprdata, exponent) );
3208 	   assert(exprdata != NULL);
3209 	
3210 	   SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPgetExprhdlrPower(scip), exprdata, 1, &child, ownercreate,
3211 	              ownercreatedata) );
3212 	
3213 	   return SCIP_OKAY;
3214 	}
3215 	
3216 	/** creates a signpower expression */
3217 	SCIP_RETCODE SCIPcreateExprSignpower(
3218 	   SCIP*                 scip,               /**< SCIP data structure */
3219 	   SCIP_EXPR**           expr,               /**< pointer where to store expression */
3220 	   SCIP_EXPR*            child,              /**< single child */
3221 	   SCIP_Real             exponent,           /**< exponent of the power expression */
3222 	   SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
3223 	   void*                 ownercreatedata     /**< data to pass to ownercreate */
3224 	   )
3225 	{
3226 	   SCIP_EXPRDATA* exprdata;
3227 	
3228 	   assert(expr != NULL);
3229 	   assert(child != NULL);
3230 	   assert(SCIPfindExprhdlr(scip, SIGNPOWEXPRHDLR_NAME) != NULL);
3231 	
3232 	   SCIP_CALL( createData(scip, &exprdata, exponent) );
3233 	   assert(exprdata != NULL);
3234 	
3235 	   SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPfindExprhdlr(scip, SIGNPOWEXPRHDLR_NAME), exprdata, 1, &child,
3236 	              ownercreate, ownercreatedata) );
3237 	
3238 	   return SCIP_OKAY;
3239 	}
3240 	
3241 	/** indicates whether expression is of signpower-type */  /*lint -e{715}*/
3242 	SCIP_Bool SCIPisExprSignpower(
3243 	   SCIP*                 scip,               /**< SCIP data structure */
3244 	   SCIP_EXPR*            expr                /**< expression */
3245 	   )
3246 	{ /*lint --e{715}*/
3247 	   assert(expr != NULL);
3248 	
3249 	   return strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), SIGNPOWEXPRHDLR_NAME) == 0;
3250 	}
3251 	
3252 	/** computes coefficients of linearization of a square term in a reference point */
3253 	void SCIPaddSquareLinearization(
3254 	   SCIP*                 scip,               /**< SCIP data structure */
3255 	   SCIP_Real             sqrcoef,            /**< coefficient of square term */
3256 	   SCIP_Real             refpoint,           /**< point where to linearize */
3257 	   SCIP_Bool             isint,              /**< whether corresponding variable is a discrete variable, and thus linearization could be moved */
3258 	   SCIP_Real*            lincoef,            /**< buffer to add coefficient of linearization */
3259 	   SCIP_Real*            linconstant,        /**< buffer to add constant of linearization */
3260 	   SCIP_Bool*            success             /**< buffer to set to FALSE if linearization has failed due to large numbers */
3261 	   )
3262 	{
3263 	   assert(scip != NULL);
3264 	   assert(lincoef != NULL);
3265 	   assert(linconstant != NULL);
3266 	   assert(success != NULL);
3267 	
3268 	   if( sqrcoef == 0.0 )
3269 	      return;
3270 	
3271 	   if( SCIPisInfinity(scip, REALABS(refpoint)) )
3272 	   {
3273 	      *success = FALSE;
3274 	      return;
3275 	   }
3276 	
3277 	   if( !isint || SCIPisIntegral(scip, refpoint) )
3278 	   {
3279 	      SCIP_Real tmp;
3280 	
3281 	      /* sqrcoef * x^2  ->  tangent in refpoint = sqrcoef * 2 * refpoint * (x - refpoint) */
3282 	
3283 	      tmp = sqrcoef * refpoint;
3284 	
3285 	      if( SCIPisInfinity(scip, 2.0 * REALABS(tmp)) )
3286 	      {
3287 	         *success = FALSE;
3288 	         return;
3289 	      }
3290 	
3291 	      *lincoef += 2.0 * tmp;
3292 	      tmp *= refpoint;
3293 	      *linconstant -= tmp;
3294 	   }
3295 	   else
3296 	   {
3297 	      /* sqrcoef * x^2 -> secant between f=floor(refpoint) and f+1 = sqrcoef * (f^2 + ((f+1)^2 - f^2) * (x-f))
3298 	       * = sqrcoef * (-f*(f+1) + (2*f+1)*x)
3299 	       */
3300 	      SCIP_Real f;
3301 	      SCIP_Real coef;
3302 	      SCIP_Real constant;
3303 	
3304 	      f = SCIPfloor(scip, refpoint);
3305 	
3306 	      coef     =  sqrcoef * (2.0 * f + 1.0);
3307 	      constant = -sqrcoef * f * (f + 1.0);
3308 	
3309 	      if( SCIPisInfinity(scip, REALABS(coef)) || SCIPisInfinity(scip, REALABS(constant)) )
3310 	      {
3311 	         *success = FALSE;
3312 	         return;
3313 	      }
3314 	
3315 	      *lincoef     += coef;
3316 	      *linconstant += constant;
3317 	   }
3318 	}
3319 	
3320 	/** computes coefficients of secant of a square term */
3321 	void SCIPaddSquareSecant(
3322 	   SCIP*                 scip,               /**< SCIP data structure */
3323 	   SCIP_Real             sqrcoef,            /**< coefficient of square term */
3324 	   SCIP_Real             lb,                 /**< lower bound on variable */
3325 	   SCIP_Real             ub,                 /**< upper bound on variable */
3326 	   SCIP_Real*            lincoef,            /**< buffer to add coefficient of secant */
3327 	   SCIP_Real*            linconstant,        /**< buffer to add constant of secant */
3328 	   SCIP_Bool*            success             /**< buffer to set to FALSE if secant has failed due to large numbers or unboundedness */
3329 	   )
3330 	{
3331 	   SCIP_Real coef;
3332 	   SCIP_Real constant;
3333 	
3334 	   assert(scip != NULL);
3335 	   assert(!SCIPisInfinity(scip,  lb));
3336 	   assert(!SCIPisInfinity(scip, -ub));
3337 	   assert(SCIPisLE(scip, lb, ub));
3338 	   assert(lincoef != NULL);
3339 	   assert(linconstant != NULL);
3340 	   assert(success != NULL);
3341 	
3342 	   if( sqrcoef == 0.0 )
3343 	      return;
3344 	
3345 	   if( SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub) )
3346 	   {
3347 	      /* unboundedness */
3348 	      *success = FALSE;
3349 	      return;
3350 	   }
3351 	
3352 	   /* sqrcoef * x^2 -> sqrcoef * (lb * lb + (ub*ub - lb*lb)/(ub-lb) * (x-lb)) = sqrcoef * (lb*lb + (ub+lb)*(x-lb))
3353 	    *  = sqrcoef * ((lb+ub)*x - lb*ub)
3354 	    */
3355 	   coef     =  sqrcoef * (lb + ub);
3356 	   constant = -sqrcoef * lb * ub;
3357 	   if( SCIPisInfinity(scip, REALABS(coef)) || SCIPisInfinity(scip, REALABS(constant)) )
3358 	   {
3359 	      *success = FALSE;
3360 	      return;
3361 	   }
3362 	
3363 	   *lincoef     += coef;
3364 	   *linconstant += constant;
3365 	}
3366 	
3367 	/** Separation for roots with exponent in [0,1]
3368 	 *
3369 	 * - x^0.5 with x >= 0
3370 	 <pre>
3371 	  8 +----------------------------------------------------------------------+
3372 	    |             +             +              +             +             |
3373 	  7 |-+                                                     x**0.5 ********|
3374 	    |                                                             *********|
3375 	    |                                                      ********        |
3376 	  6 |-+                                             ********             +-|
3377 	    |                                         ******                       |
3378 	  5 |-+                                 ******                           +-|
3379 	    |                             ******                                   |
3380 	    |                        *****                                         |
3381 	  4 |-+                  ****                                            +-|
3382 	    |               *****                                                  |
3383 	  3 |-+          ****                                                    +-|
3384 	    |         ***                                                          |
3385 	    |      ***                                                             |
3386 	  2 |-+  **                                                              +-|
3387 	    |  **                                                                  |
3388 	  1 |**                                                                  +-|
3389 	    |*                                                                     |
3390 	    |*            +             +              +             +             |
3391 	  0 +----------------------------------------------------------------------+
3392 	    0             10            20             30            40            50
3393 	 </pre>
3394 	 */
3395 	void SCIPestimateRoot(
3396 	   SCIP*                 scip,               /**< SCIP data structure */
3397 	   SCIP_Real             exponent,           /**< exponent */
3398 	   SCIP_Bool             overestimate,       /**< should the power be overestimated? */
3399 	   SCIP_Real             xlb,                /**< lower bound on x */
3400 	   SCIP_Real             xub,                /**< upper bound on x */
3401 	   SCIP_Real             xref,               /**< reference point (where to linearize) */
3402 	   SCIP_Real*            constant,           /**< buffer to store constant term of estimator */
3403 	   SCIP_Real*            slope,              /**< buffer to store slope of estimator */
3404 	   SCIP_Bool*            islocal,            /**< buffer to store whether estimator only locally valid, that is,
3405 	                                                  it depends on given bounds */
3406 	   SCIP_Bool*            success             /**< buffer to store whether estimator could be computed */
3407 	   )
3408 	{
3409 	   assert(scip != NULL);
3410 	   assert(constant != NULL);
3411 	   assert(slope != NULL);
3412 	   assert(islocal != NULL);
3413 	   assert(success != NULL);
3414 	   assert(exponent > 0.0);
3415 	   assert(exponent < 1.0);
3416 	   assert(xlb >= 0.0);
3417 	
3418 	   if( !overestimate )
3419 	   {
3420 	      /* underestimate -> secant */
3421 	      computeSecant(scip, FALSE, exponent, xlb, xub, constant, slope, success);
3422 	      *islocal = TRUE;
3423 	   }
3424 	   else
3425 	   {
3426 	      /* overestimate -> tangent */
3427 	
3428 	      /* need to linearize right of 0 */
3429 	      if( xref < 0.0 )
3430 	         xref = 0.0;
3431 	
3432 	      if( SCIPisZero(scip, xref) )
3433 	      {
3434 	         if( SCIPisZero(scip, xub) )
3435 	         {
3436 	            *success = FALSE;
3437 	            *islocal = FALSE;
3438 	            return;
3439 	         }
3440 	
3441 	         /* if xref is 0 (then xlb=0 probably), then slope is infinite, then try to move away from 0 */
3442 	         if( xub < 0.2 )
3443 	            xref = 0.5 * xlb + 0.5 * xub;
3444 	         else
3445 	            xref = 0.1;
3446 	      }
3447 	
3448 	      computeTangent(scip, FALSE, exponent, xref, constant, slope, success);
3449 	      *islocal = FALSE;
3450 	   }
3451 	}
3452 	
3453 	/* from pub_expr.h */
3454 	
3455 	/** gets the exponent of a power or signed power expression */  /*lint -e{715}*/
3456 	SCIP_Real SCIPgetExponentExprPow(
3457 	   SCIP_EXPR*            expr                /**< expression */
3458 	   )
3459 	{
3460 	   SCIP_EXPRDATA* exprdata;
3461 	
3462 	   assert(expr != NULL);
3463 	
3464 	   exprdata = SCIPexprGetData(expr);
3465 	   assert(exprdata != NULL);
3466 	
3467 	   return exprdata->exponent;
3468 	}
3469